source: trunk/MagicSoft/Simulation/Detector/include-MTrigger/MTrigger.cxx@ 393

/////////////////////////////////////////////////////////////////
//
//  MTrigger
//
//  
#include "MTrigger.hxx"

#include "TROOT.h"
#include "TFile.h"
#include "TH1.h"
#include "TObjArray.h"
#include "MGTriggerSignal.hxx"


MTrigger::MTrigger() { 
  // ============================================================
  //
  //  default constructor 
  //  
  //  The procedure is the following: 
  //
  //  1. Allocation of some memory needed
  //  2. some parameters of the trigger are set to default.   
  //  3. if a File MTrigger.card exists in the current directory, 
  //     this parameters of the trigger may be changed
  //  4. Then the all signals are set to zero

  FILE *unit_mtrig ; 
  Int_t endflag = 1  ; 
  Int_t bthresholdpixel = FALSE;
  char   datac[256] ; 
  char   dummy[50] ; 
  char   input_thres[50];
  Int_t  i, ii ; 

  Float_t threshold ; 
 
  // 
  //   allocate the memory for the 2dim arrays (a_sig, d_sig ) 
  //

  for( Int_t j=0; j<TRIGGER_PIXELS; j++ ) { 

    a_sig[j] = new Float_t[TRIGGER_TIME_SLICES] ; 

    d_sig[j] = new Float_t[TRIGGER_TIME_SLICES] ; 
  } 

  //
  //   set the values for the standard response pulse
  //

  fwhm_resp = RESPONSE_FWHM       ; 
  ampl_resp = RESPONSE_AMPLITUDE  ; 
    
  threshold = CHANNEL_THRESHOLD  ; 
  

  gate_leng        = TRIGGER_GATE       ; 
  trigger_multi    = TRIGGER_MULTI      ; 
  trigger_geometry = TRIGGER_GEOM ; 

  //
  //  check if the file MTrigger.card exists
  //

  if ( (unit_mtrig = fopen ("MTrigger.card", "r")) != 0 ) {
    cout << "[MTrigger]  use the values from MTrigger.card "<< endl ; 
    
    while ( endflag == 1  ) {
      //
      //
      fgets (datac, 255, unit_mtrig) ;
      //      printf ("--> %s <--", datac ) ;
      
      //
      //   now compare the line with controlcard words
      //
      
      if (      strncmp (datac, "channel_threshold", 17 ) == 0 ) {
        sscanf (datac, "%s %f", dummy, &threshold ) ;
      }
      else if ( strncmp (datac, "gate_length", 11 ) == 0 ) {
        sscanf (datac, "%s %f", dummy, &gate_leng ) ;
      }
      else if ( strncmp (datac, "response_fwhm", 13 ) == 0 ) {
        sscanf (datac, "%s %f", dummy, &fwhm_resp ) ;
      }
      else if ( strncmp (datac, "response_ampl", 13 ) == 0 ) {
        sscanf (datac, "%s %f", dummy, &ampl_resp ) ;
      }
      else if ( strncmp (datac, "multiplicity", 12 ) == 0 ) {
        sscanf (datac, "%s %f", dummy, &trigger_multi ) ;
      }
      else if ( strncmp (datac, "topology", 8 ) == 0 ) {
        sscanf (datac, "%s %i", dummy, &trigger_geometry ) ;
      }
      else if ( strncmp (datac, "threshold_file", 14 ) == 0 ) {
        sscanf (datac, "%s %s", dummy, input_thres ) ;
        bthresholdpixel=TRUE;
      }

      if ( feof(unit_mtrig) != 0 ) {
        endflag = 0 ;
      }
      
    }

    fclose ( unit_mtrig ) ; 
  }
  else {
    cout << "[MTrigger]  use the standard values for MTrigger "<< endl ; 
  }

  cout << endl
       << "[MTrigger]  Setting up the MTrigger with this values "<< endl ; 
  if(bthresholdpixel){
    cout<<endl
        << "[MTrigger]    ChannelThreshold from file:   "<<input_thres
        <<endl;
  }
  else{
    cout << endl 
         << "[MTrigger]    ChannelThreshold:   " << threshold << " mV" 
         << endl ; 
  }
  cout << "[MTrigger]    Gate Length:        " << gate_leng  << " ns"
       << endl ; 
  cout << "[MTrigger]    Response FWHM:      " << fwhm_resp  << " ns"
       << endl ; 
  cout << "[MTrigger]    Response Amplitude: " << ampl_resp  << " mV"
       << endl ; 
  cout << "[MTrigger]    Trigger Multiplicity:      " << trigger_multi  << " pixels"
       << endl ; 
  cout << "[MTrigger]    Trigger Topology: " << trigger_geometry 
       << endl ; 
  
  cout << endl ; 


  //
  //   we have introduced individual thresholds for all pixels
  //
  FILE *unit_thres;

  if (bthresholdpixel == TRUE) {
    if ((unit_thres=fopen(input_thres, "r"))==0){
      cout<<"WARNING: not able to read ..."<<input_thres<<endl;
      cout<<"Threshold will be set to "<<threshold<<" for all pixels"<<endl;
      for (Int_t k=0; k<TRIGGER_PIXELS; k++ ) {
        chan_thres[k] = threshold ; 
      }
    }
    else {
      for (i=0;i<TRIGGER_PIXELS;i++){
        fscanf(unit_thres, "%f",&chan_thres[i]);
      }
      fclose (unit_thres);
    }
  }
  else {
    for (Int_t k=0; k<TRIGGER_PIXELS; k++ ) {
      chan_thres[k] = threshold ; 
    }
  }


  //
  //    set up the response shape
  // 
     
  Float_t   sigma ; 
  Float_t   x, x0 ; 

  sigma = fwhm_resp / 2.35 ; 
  x0 = 3*sigma ; 
  
  for (i=0; i< RESPONSE_SLICES ; i++ ) {  

    x = i * (1./((Float_t)SLICES_PER_NSEC)) 
      + (1./( 2 * (Float_t)SLICES_PER_NSEC ))  ; 
    
    sing_resp[i] = 
      ampl_resp * expf(-0.5 * (x-x0)*(x-x0) / (sigma*sigma) ) ; 
  
  }

  //
  //    look for the time between start of response function and the 
  //    maximum value of the response function. This is needed by the
  //    member functions FillNSB() and FillStar() 
  //  

  Int_t imax  = 0  ;
  Float_t max = 0. ; 
  for (i=0; i< RESPONSE_SLICES ; i++ ) {  
    if ( sing_resp[i] > max ) {
      imax = i ; 
      max  = sing_resp[i] ; 
    }
  }
 
  peak_time = ( (Float_t) imax )  / ( (Float_t) SLICES_PER_NSEC ) ;
 

  //
  //   the amplitude of one single photo electron is not a constant. 
  //   There exists a measured distribution from Razmik. This distribution
  //   is used to simulate the noise of the amplitude. 
  //   For this a histogramm (histPmt) is created and filled with the
  //   values. 
  // 

  histPmt = new TH1F ("histPmt","Noise of PMT", 40, 0., 40.) ;
  
  Stat_t ValRazmik[41] = { 0., 2.14, 2.06, 2.05, 2.05, 2.06, 2.07, 2.08,  2.15,
                           2.27, 2.40, 2.48, 2.55, 2.50, 2.35, 2.20,  2.10,
                           1.90, 1.65, 1.40, 1.25, 1.00, 0.80, 0.65,  0.50,
                           0.35, 0.27, 0.20, 0.18, 0.16, 0.14, 0.12,  0.10, 
                           0.08, 0.06, 0.04, 0.02, 0.01, 0.005,0.003, 0.001} ; 

  histMean =  histPmt->GetMean() ;   
  
  histPmt->SetContent( ValRazmik) ; 

  histMean =  histPmt->GetMean() ; 

  //
  //   create the random generator for the Electronic Noise
  // 

  GenElec = new TRandom() ; 

  //
  //  Read in the lookup table for NN trigger
  //
  
  FILE *unit ;
  int id ;

  i = 0 ; 

  if ( (unit = fopen("../include-MTrigger/TABLE_NEXT_NEIGHBOUR", "r" )) == 0 ) { 
    cout << "ERROR: not able to read ../include-MTrigger/TABLE_NEXT_NEIGHBOUR"
         << endl ; 
    exit(123) ; 
  }  
  else { 
    while ( i < TRIGGER_PIXELS )
      {
        fscanf ( unit, " %d", &id ) ;
        
        for ( Int_t k=0; k<6; k++ ) { 
          fscanf ( unit, "%d ", &NN[i][k]  ) ; 
        }
        i++ ;
      }
    
    fclose (unit) ;
  }


  //
  //
  //  set all the booleans used to FALSE, indicating that the pixel is not 
  //  used in this event. 
  //
  
  for ( i =0 ; i <TRIGGER_PIXELS ; i++ ) { 
    used [i] = FALSE ;
    dknt [i] = FALSE ;

    nphotshow[i] = 0 ;
    nphotnsb [i] = 0 ;
    nphotstar[i] = 0 ;

    baseline[i] = 0 ; 
  }

  for ( ii=0 ; ii<TRIGGER_TIME_SLICES; ii++ ) { 
    sum_d_sig[ii] = 0. ;  
  }
  
  //
  //   set the information about the Different Level Triggers to zero
  //

  nZero = nFirst = nSecond = 0 ; 

  for (ii=0 ; ii<TRIGGER_TIME_SLICES; ii++ ) { 
    SlicesZero[ii] = FALSE;
  }
  
  for ( i = 0 ; i < 5 ; i++) {
    SlicesFirst[i]  = -50.0 ; 
    SlicesSecond[i] = -50.0 ;
    PixelsFirst[i]  = -1;
    PixelsSecond[i] = -1;
  }
  cout << " end of MTrigger::MTrigger()" << endl ; 
} 

MTrigger::~MTrigger() {
  // ============================================================//
  //   destructor
  //
  int i;
  // delete histPmt ; 
  for(i=0;i<TRIGGER_PIXELS;i++){
    delete [] a_sig[i];
    delete [] d_sig[i];
  }
  delete GenElec;
}

  
void MTrigger::Reset() { 
  // ============================================================
  //
  //  reset all values of the signals to zero
  //
  Int_t  i, ii ; 
  
  for ( i =0 ; i <TRIGGER_PIXELS ; i++ ) {
    used [i] = FALSE ; 
    dknt [i] = FALSE ; 
    
    nphotshow[i] = 0 ; 
    nphotnsb [i] = 0 ; 
    nphotstar[i] = 0 ; 
    } 
  
  for ( ii=0 ; ii<TRIGGER_TIME_SLICES; ii++ ) { 
    sum_d_sig[ii] = 0. ;  
  }    
}

void MTrigger::ClearZero() {
  //
  //   set the information about the Zero Level Trigger to zero
  //

  Int_t i;

  nZero = 0 ; 

  for (i=0 ; i<TRIGGER_TIME_SLICES; i++ ) { 
    SlicesZero[i] = FALSE;
  }
  
}

void MTrigger::ClearFirst() {
  //
  //   set the information about the First Level Trigger to zero
  //

  Int_t i;

  nFirst =  0 ; 
  
  for ( i = 0 ; i < 5 ; i++) {
    SlicesFirst[i]  = -50.0 ; 
    PixelsFirst[i]  = -1;
  }
}

Float_t  MTrigger::FillShow(Int_t iPix, Float_t time) { 
  // ============================================================
  //
  //     Fills the information of one single Phe electron that 
  //     comes from the shower 
  //
  
  //
  //  First check the time 
  //
  
  if ( time < 0. || time > TOTAL_TRIGGER_TIME ) {
    cout << "     WARNING: time of phe out of time range: " << time << endl; 
    return 0. ;      
  }
  else {
    return ( Fill( iPix, time, CASE_SHOW ) )  ; 
  }
}

Float_t  MTrigger::FillNSB(Int_t iPix, Float_t time) { 
  // ============================================================
  //
  //     Fills the information of one single Phe electron that 
  //     comes from the shower 
  //

  //
  //  First check the time 
  //

  if ( time < 0. || time > TOTAL_TRIGGER_TIME ) {
    cout << "     WARNING: time of phe out of time range: " << time << endl; 
    return 0. ;      
  }
  else {
    return ( Fill( iPix, time - peak_time, CASE_NSB ) )  ; 
  }
}

Float_t  MTrigger::FillStar(Int_t iPix, Float_t time) { 
  // ============================================================
  //
  //     Fills the information of one single Phe electron that 
  //     comes from the shower 
  //

  //
  //  First check the time 
  //

  if ( time < 0. || time > TOTAL_TRIGGER_TIME ) {
    cout << "     WARNING: time of phe out of time range: " << time << endl; 
    return 0. ;      
  }
  else {
    return ( Fill( iPix, time - peak_time, CASE_STAR ) )  ; 
  }
}

Float_t MTrigger::Fill( Int_t iPix, Float_t time, Int_t fall ) {
  // ============================================================
  //
  //     Fills the information in the array for the analog signal
  //

  Float_t PmtAmp = 0 ;   //  Amplitude of the PMT signal (results from noise)

  if ( iPix < 0 ) {
    cout << "     ERROR: in MTrigger::Fill()    " << endl ; 
    cout << "     ERROR: Pixel Id < 0 ---> Exit " << endl ; 
    exit (1) ; 
  }
  else if ( iPix >= CAMERA_PIXELS ) {
    cout << "     ERROR: in MTrigger::Fill()    " << endl ; 
    cout << "     ERROR: Pixel Id > CAMERA_PIXELS ---> Exit " << endl ; 
    exit (1) ; 
  }
  else if ( iPix >= TRIGGER_PIXELS ) {
    //
    //  We have not to fill information in the trigger part, 
    //  but we must create the height of the puls going into 
    //  the FADC simulation
    //
    PmtAmp =  (histPmt->GetRandom()/histMean) ; 
    
    //
    //  But we fill the information in the counters of phe's
    //
    
    if ( fall == CASE_SHOW )
      nphotshow[iPix]++ ; 
    else if ( fall == CASE_NSB ) 
      nphotshow[iPix]++ ; 
    else if ( fall == CASE_STAR ) 
      nphotstar[iPix]++ ; 
    
      
  }
  else {  
    //
    // we have a trigger pixel and we fill it 
    //
    Int_t i ; 

    //
    //  but at the beginning we must check if this pixel is
    //  hitted the first time
    //
    
    if ( used[iPix] == FALSE ) {
      used [iPix] = TRUE ; 
      //      baseline[iPix] = 0. ; 
      
      for (i=0; i < TRIGGER_TIME_SLICES; i++ ) {
        a_sig[iPix][i] = 0. ; 
        d_sig[iPix][i] = 0. ; 
      }
    }
    
    //
    // get the randomized amplitude
    //
    PmtAmp =  (histPmt->GetRandom()/histMean) ; 

    //
    // select the first slice to fill
    //

    Int_t ichan = (Int_t) ( time * ((Float_t) SLICES_PER_NSEC) ) ; 

    //
    //  look over the response signal and put it in the signal line
    //

    for ( i = 0 ; i<RESPONSE_SLICES; i++ ) {
      
      if ( (ichan+i) >= 0  &&  
           (ichan+i)  < TRIGGER_TIME_SLICES ) {  
        a_sig[iPix][ichan+i] += PmtAmp * sing_resp[i] ;
      } 
    }

    //
    //  we fill the information in the counters of phe's
    //
    
    if ( fall == CASE_SHOW )
      nphotshow[iPix]++ ; 
    else if ( fall == CASE_NSB ) 
      nphotshow[iPix]++ ; 
    else if ( fall == CASE_STAR ) 
      nphotstar[iPix]++ ; 

    //
    //
    return PmtAmp ;
  }
  return PmtAmp ; 
}



void MTrigger::ElecNoise() {
  // ============================================================
  //
  //    adds the noise due to optronic and electronic 
  //    to the signal
  //
  Float_t rausch ; 

  rausch = RESPONSE_AMPLITUDE * 0.3 ; 
 
  for ( Int_t i=0 ; i < TRIGGER_PIXELS; i++  ) {
    if ( used [i] == TRUE ) {
       
      for ( Int_t ii=1 ; ii<TRIGGER_TIME_SLICES; ii++ ) { 

        a_sig [i][ii] +=  GenElec->Gaus(0., rausch  ) ; 

      }
    }
  }
}

void MTrigger::SetMultiplicity(Int_t multi){
  //=============================================================
  //
  //    It sets the private member trigger_multi

  trigger_multi=multi;
}

void MTrigger::SetTopology(Int_t topo){
  //=============================================================
  //
  //    It sets the private member trigger_geometry

  trigger_geometry=topo;
}

void MTrigger::SetThreshold(Float_t thres[]){
  //=============================================================
  //
  //    It sets the private member chan_thres[TRIGGER_PIXELS]

  Int_t i;

  for(i=0;i<TRIGGER_PIXELS;i++){
    chan_thres[i]=thres[i];
  }
}

void MTrigger::ReadThreshold(char name[]){
  //=============================================================
  //
  //   It reads values for threshold of each pixel from file name

  FILE *unit;
  Int_t i=0;

  if ((unit=fopen(name, "r"))==0){
    cout<<"WARNING: not able to read ..."<<name<<endl;
  }
  else {
    while (i<TRIGGER_PIXELS){
      fscanf(unit, "%f",&chan_thres[i++]);
    }
    fclose (unit);
  }
        
}

void MTrigger::Diskriminate() {
  // ============================================================
  //
  //    Diskriminates the analog signal
  //
  //    one very important part is the calucaltion of the baseline
  //    shift. Because of the AC coupling of the PMT, only the 
  //    fluctuations are interesting. If there are a lot of phe, 
  //    a so-called shift of the baseline occurs.
  //
 
  Int_t  iM = 0 ; 
  Int_t  i, ii  ; 


  Int_t  jmax = (Int_t) (gate_leng * SLICES_PER_NSEC )  ; 

  //
  //    first of all determine the integral of all signals to get
  //    the baseline shift. 
  // 

  for ( i=0 ; i < TRIGGER_PIXELS ; i++ ) {
    if ( used[i] == TRUE ) {
      baseline[i] = 0. ; 

      for ( ii = 0 ;  ii < TRIGGER_TIME_SLICES ; ii++ ) {
        baseline[i] += a_sig[i][ii] ; 
      } 

      baseline[i] = baseline[i] / ( (Float_t ) TRIGGER_TIME_SLICES)  ;

      //
      //  now correct the baseline shift in the analog signal!!
      //
      for ( ii = 0 ;  ii < TRIGGER_TIME_SLICES ; ii++ ) {
        a_sig[i][ii] = a_sig[i][ii] - baseline[i] ; 
      } 
    }
  }
  
  //  
  //  now the diskrimination is coming
  //
  //  take only that pixel which are used
  //

  for ( i=0 ; i < TRIGGER_PIXELS; i++  ) {
    if ( used [i] == TRUE ) {
           
      for ( ii=1 ; ii<TRIGGER_TIME_SLICES; ii++ ) { 
        //
        // first check if the signal is crossing the CHANNEL_THRESHOLD 
        // form low to big signals
        //
        
        if ( a_sig[i][ii-1] <   chan_thres[i]  && 
             a_sig[i][ii]   >=  chan_thres[i]  ) { 
          { 
            if ( dknt[i] == FALSE ) {
              dknt [i] = TRUE ; 
              iM++ ; 
            }
            //    cout << " disk " << ii  ; 
            // 
            //   put the standard diskriminator signal in 
            //   the diskriminated signal 
            //
            for ( Int_t j=0 ; j < jmax ; j++ ) { 
              
              if ( ii+j  < TRIGGER_TIME_SLICES ) { 
                d_sig  [i][ii+j] = 1. ;  
              } 
            } 
            ii = ii + jmax ; 
          }
        }
        else d_sig[i][ii]=0.;
      }
    }
  }  
}


void MTrigger::ShowSignal (MMcEvt *McEvt) { 
  // ============================================================
  //
  //  This method is used to book the histogramm to show the signal in 
  //  a special gui frame (class MGTriggerSignal). After the look onto the
  //  signals for a better understanding of the things we will expect 
  //  the gui frame and all histogramms will be destroyed.   
  //
  
  //
  //  first of all create a list of the histograms to show
  //
  //  take only that one with a entry
  
  TH1F *hist ; 
  TH1F *dhist ; 
  Char_t dumm[10]; 
  Char_t name[256]; 
  
  TObjArray  *AList ;
  AList = new TObjArray(10) ; 

  TObjArray  *DList ;
  DList = new TObjArray(10) ; 

  // the list of analog signal histograms
  // at the beginning we initalise 10 elements
  // but this array expand automaticly if neccessay
  
  Int_t ic = 0 ; 
  for ( Int_t i=0 ; i < TRIGGER_PIXELS; i++  ) {
    if ( used [i] == TRUE ) {

      sprintf (dumm, "A_%d", i ) ; 
      sprintf (name, "analog %d", i ) ; 
      
      hist = new TH1F(dumm, name, TRIGGER_TIME_SLICES, 0., TOTAL_TRIGGER_TIME); 
      //
      //  fill the histogram
      //

      for (Int_t ibin=1; ibin <=TRIGGER_TIME_SLICES; ibin++) {
        hist->SetBinContent (ibin, a_sig[i][ibin-1]) ;
      }
      hist->SetMaximum(8.); 
      hist->SetMinimum(-8.); 
      hist->SetStats(kFALSE); 
  
      AList->Add(hist) ; 
      
      sprintf (dumm, "D_%d", i ) ; 
      sprintf (name, "digital %d", i ) ; 
      
      dhist = new TH1F(dumm, name, TRIGGER_TIME_SLICES, 0., TOTAL_TRIGGER_TIME); 
      if ( dknt[i] == TRUE ) {
        //
        //   fill the histogram of digital signal
        //
        for (Int_t ibin=1; ibin <=TRIGGER_TIME_SLICES; ibin++) {
          dhist->SetBinContent (ibin, d_sig[i][ibin-1]) ;
          dhist->SetStats(kFALSE); 
        } 
      }
      dhist->SetMaximum(1.5);
        
      DList->Add(dhist); 
      
      ic++ ; 

    }
  }  

  //
  //   create the Gui Tool
  //
  //

  new MGTriggerSignal(McEvt,
                      AList,
                      DList,
                      gClient->GetRoot(), 
                      gClient->GetRoot(), 
                      400, 400 ) ; 
  
  //
  //   delete the List of histogramms
  //

  AList->Delete() ; 
  DList->Delete() ; 
  
  delete AList ; 
  delete DList ; 
}
  

Int_t MTrigger::ZeroLevel() {
  // ============================================================
  //
  //  This is a level introduced just to speed up the program. 
  //  It makes sense to look for next neighbours only if there
  //  are at least trigger_multi  pixels with a diskriminator 
  //  signal. 
  //
  
  //
  //  first count the pixels with a diskriminator signal
  //
  Int_t iMul = 0 ; 
  for ( Int_t iP =0 ; iP < TRIGGER_PIXELS; iP++ ) {
    //
    //
    if ( dknt[iP] == TRUE ) {
      iMul++ ; 
    }
  }

  //
  //   only if there are at least more pixels than requested 
  //   it make sense to look into details
  if ( iMul >= trigger_multi ) {
    //
    //  fill the sum signal of all diskriminator signals
    // 
    for ( Int_t iP =0 ; iP < TRIGGER_PIXELS; iP++ ) {
      //
      //
      if ( dknt[iP] == TRUE ) {
        //
        // sum it up
        //
        for (Int_t iS=0; iS< TRIGGER_TIME_SLICES; iS++ ) {
          //
          //    
          sum_d_sig [iS] += d_sig[iP][iS] ; 
        } 
      }
    }
    //
    //   run over the sum_d_sig and check each time slice
    //
    Int_t iReturn = 0 ; 
  
    for (Int_t iS=0; iS< TRIGGER_TIME_SLICES; iS++ ) {
      
      if ( sum_d_sig[iS] >= trigger_multi ) {
        iReturn++ ; 
        nZero++;
        SlicesZero[iS] = TRUE ; 
        
      }
      else SlicesZero[iS] = FALSE;
    }

    return ( iReturn ) ; 
  }
  else {
    return 0 ; 
  }
}

Int_t MTrigger::FirstLevel() {
  //=================================================
  //
  // This is a level trigger which can look for several 
  // multiplicities (trigger_multi) 
  // and topologies (trigger_geometry)
  //

  Int_t iReturn = 0 ;  // Return value for this function 
    
  //   Definition of needed variables  
  Bool_t Muster[TRIGGER_PIXELS] ; 
  Bool_t Neighb[TRIGGER_PIXELS] ; 
  Int_t  iMulti = 0 ; 

  // We put several wrong topologies which we already know that they
  // are not possible. It can save time.

  if (trigger_geometry==0 &&  trigger_multi>7) {
    cout <<"You are looking for a topology that needs more than six neighbours of the same pixel"<<endl;
    cout <<" Topology   "<<trigger_geometry<<"   Multiplicity   "<<trigger_multi<<endl;;
    return (kFALSE);
  }

  if (trigger_geometry==2 && trigger_multi<3) {
    cout<<"Closed pack geometry with multiplicity "<<trigger_multi<<" does not make sense"<<endl;
    return (kFALSE);
  }
  if (trigger_geometry>2) {
    cout << "This trigger topology is not implemented"<<endl;
    return (kFALSE);
  }
  
  //
  //   loop over all ZeroLevel Trigger
  //
  //   it is only neccessary to look after a ZeroLevel Trigger for
  //   a FirstLevel (NextNeighbour) trigger. 
  //  

  if (nZero) {
    
    //
    //   Then run over all slices 
    //
    
    for ( Int_t iSli = 0;
          iSli < TRIGGER_TIME_SLICES; iSli++ ) {

      //  Check if this time slice has more fired pixels than trigger_multi

      if (SlicesZero[iSli]){
        //
        //  then look in all pixel if the diskriminated signal is 1
        //
        
        for ( Int_t iPix = 0 ; iPix < TRIGGER_PIXELS; iPix++ ) {
          Muster[iPix] = kFALSE ; 
          Neighb[iPix] = kFALSE ;
          if ( used [iPix] == TRUE ) {
            //
            //  now check the diskriminated signal
            //
            if ( d_sig [iPix][iSli] > 0. ) {
              Muster[iPix] = kTRUE ; 
            } 
          }
        } // end of loop over the pixels
        
        //
        //   here we have to look for the topologies
        //
        
        switch(trigger_geometry){
        case 0:{
          
          // It looks for a pixel above threshold which has 
          // trigger_multi-1 neighbour pixels above threshold 

          Bool_t Dummy[TRIGGER_PIXELS] ;

          //  Loop over all pixels
          for (int j=0;j<TRIGGER_PIXELS;j++){

            //
            //  I commented out this line here
            //      Dummy=Muster;

            for (int k=0; k<TRIGGER_PIXELS; k++){
              Neighb[k]=kFALSE;

              Dummy[k]  = Muster[k] ; 
            }
            if(Muster[j]){
              //  If pixel is fired, it checks how many fired neighbours it has
              for (iMulti=1;iMulti<trigger_multi; iMulti++) {
                Neighb[j] = kTRUE ;
                Dummy[j] = kTRUE ;
                if (!PassNextNeighbour(Dummy, &Neighb[0])){
                  break;
                } 
                for (int k=0; k<TRIGGER_PIXELS; k++){
                  if (Neighb[k]){
                    Dummy[k]=kFALSE;
                    Neighb[k]=kFALSE;
                  }
                }
              }
              if (iMulti==trigger_multi ) { 
                //
                //   A NN-Trigger is detected at time Slice 
                //
                PixelsFirst[nFirst] = j;  //  We save pixel that triggers
                SlicesFirst[nFirst++] = iSli ; // We save time when it triggers
                iReturn++ ;
                iSli+=(50*SLICES_PER_NSEC);  // We skip the following 50 ns (dead time)
                break ; 
              }   
            }
          }
          break;
        };
        
        case 1:{

          //   It looks for trigger_multi neighbour pixels above the 
          //   threshold.
          
          for (int j=0;j<TRIGGER_PIXELS;j++){
            if(Muster[j]){
              //  It checks if you can find 
              //  trigger_multi fired neighbour pixels
              Neighb[j] = kTRUE ;
              for (iMulti=1;iMulti<trigger_multi; iMulti++) {
                if (!PassNextNeighbour(Muster, &Neighb[0]))
                  break; 
              }
              if (iMulti==trigger_multi ) { 
                //
                //   A NN-Trigger is detected at time Slice 
                //
                PixelsFirst[nFirst] = j;  //  We save pixel that triggers
                SlicesFirst[nFirst++] = iSli ; //  We save when it triggers
                iReturn++ ;
                iSli+=(50*SLICES_PER_NSEC);  // We skip the following 50 ns (dead time)
                break ; 
              }   
              else {
                //  We put Neighb to kFALSE to check an other pixel
                for (int k=0; k<TRIGGER_PIXELS; k++){
                  if (Neighb[k]){
                    Neighb[k]=kFALSE;
                  }
                }
              }
            }
          }
          break;
        };
        case 2:{
          
          // It looks for trigger_multi closed pack neighbours
          // above threshold
          // Closed pack means that you can take out any pixel
          // and you will still get a trigger for trigger_multi -1

          Int_t closed_pack = 1;

          for (int j=0;j<TRIGGER_PIXELS;j++){
            if(Muster[j]){
              //  It checks if there are trigger_multi
              //  neighbours above threshold
              Neighb[j] = kTRUE ;
              for (iMulti=1;iMulti<trigger_multi; iMulti++){
                if (!PassNextNeighbour(Muster, &Neighb[0]))
                  break; 
              }
              if (iMulti==trigger_multi ) { 
                //
                //   A NN-Trigger is detected at time Slice 
                //

                // Check if there is closed pack topology
                Bool_t Aux1[TRIGGER_PIXELS];
                Bool_t Aux2[TRIGGER_PIXELS];
                for (int jj=0;jj<TRIGGER_PIXELS;jj++)
                  Aux2[jj]=kFALSE;
                
                for (int i=0;i<TRIGGER_PIXELS;i++){
                  if (Neighb[i]) {
                    //  Loop over pixels that achive neighbouring condition
                    
                    // huschel
                    // here also some corrections were neccessary
                    //
                    //Aux1=Neighb;

                    for (int jj=0;jj<TRIGGER_PIXELS;jj++) {
                      Aux1[jj] = Neighb[jj] ;   // huschel

                      Aux2[jj]=kFALSE;
                    }

                    Aux1[i]=kFALSE;
                    Aux2[j]=kTRUE;
                    //  It checks if taking any of the pixels we lose
                    //  neighbouring condition for trigger -1
                    for (iMulti=1;iMulti<(trigger_multi-1);iMulti++){
                      if (!PassNextNeighbour(Aux1, &Aux2[0]))
                        break; 
                    }
                    if (iMulti<(trigger_multi-1)){
                      closed_pack=0;
                      break;
                    }

                  }
                }
                if (closed_pack){
                  PixelsFirst[nFirst] = j;  //  We save pixel that triggers
                  SlicesFirst[nFirst++] = iSli ; //  We save time when it triggers
                  iReturn++ ;
                  iSli+=(50*SLICES_PER_NSEC);  // We skip the following 50 ns (dead time)
                  break ; 
                }
                else {
                  for (int k=0; k<TRIGGER_PIXELS; k++){
                    if (Neighb[k]){
                      Neighb[k]=kFALSE;
                   }
                  }
                }
              }
              else
                for (int k=0; k<TRIGGER_PIXELS; k++)
                  Neighb[k]=kFALSE;               
            }
          }
          break;
        };
        default:{
          cout << "This topology is not implemented yet"<<endl;
          break;
        }
        }
      }
    } // end of loop over the slices  
  }   // end of conditional for a trigger Zero
  
  //
  //   return the Number of FirstLevel Triggers
  //
  return iReturn ; 
}


Bool_t MTrigger::PassNextNeighbour ( Bool_t m[], Bool_t *n) {
  //
  //  This function is looking for a next neighbour of pixels in n[]
  //  above triggers using a NNlookup table. 
  //  This table is builded by the default constructor
  //

  //
  //   loop over all trigger pixels
  //
  
  Bool_t return_val = kFALSE;
  
  for ( Int_t i=0; i<TRIGGER_PIXELS; i++) { 
    //
    //  check if this pixel has a diskrminator signal 
    //  (this is inside n[] ) 
    //
    
    if ( n[i] && !return_val) { 
      
      //
      //  look in the next neighbours from the lookuptable
      //
      
      for ( Int_t kk=0; kk<6; kk++ ) { 
        //
        //  if the nextneighbour is outside the triggerarea do nothing
        //
        if (!return_val){
          if (NN[i][kk] >= TRIGGER_PIXELS ) {
            
          } 
          // the nextneighbour is not inside the TRIGGER_PIXELS
          else { 
            //
            //  look if the boolean of nn pixels is true
            //
            
            if ( m[ NN[i][kk] ] && !n[NN[i][kk]] ) { 
              n[NN[i][kk]]=kTRUE ;
              return_val =kTRUE;
            }  
          }
        }
        else break;
      }
    } 
  }    
  return(return_val);
}

Float_t MTrigger::GetFirstLevelTime( Int_t il ){

  //=============================================================
  //
  //  It gives the time for the il trigger at first level

  return((Float_t) ((Float_t) SlicesFirst[il]/((Float_t) SLICES_PER_NSEC)));
}

Int_t MTrigger::GetFirstLevelPixel( Int_t il ){

  //=============================================================
  //
  //  It gives the pixel that triggers for the il trigger at first level
  return(PixelsFirst[il]);
}