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

/////////////////////////////////////////////////////////////////
//
//  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  ; 

  overlaping_time = TRIGGER_OVERLAPING;

  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, "overlaping", 10 ) == 0 ) {
        sscanf (datac, "%s %f", dummy, &overlaping_time ) ;
      }
      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]    Overlaping time:    " << overlaping_time << " 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) ;
  }


  //
  //  Read in the lookup table for trigger cells
  //
  
  i = 0 ; 

  if ( (unit = fopen("../include-MTrigger/TABLE_PIXELS_IN_CELLS", "r" )) == 0 ) { 
    cout << "ERROR: not able to read ../include-MTrigger/TABLE_PIXELS_IN_CELLS"
         << endl ; 
    exit(123) ; 
  }  
  else { 
    while ( i < TRIGGER_PIXELS )
      { 
        for ( Int_t k=0; k<TRIGGER_CELLS; k++ ) { 
          TC[k][i]=FALSE; 
        }
        i++ ;
      }
    while ( feof(unit) == 0 ) {
        for ( Int_t k=0; k<TRIGGER_CELLS; k++ ) { 
          fscanf ( unit, "%d ", &i ) ;
          if ((i-1)<TRIGGER_PIXELS)
            TC[k][i-1]=TRUE;
        }
    }
    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 ; 
    SlicesSecond[i] = -50 ;
    PixelsFirst[i]  = -1;
    PixelsSecond[i] = -1;
  }
  cout << " end of MTrigger::MTrigger()" << endl ; 
} 

MTrigger::MTrigger(float gate, float overt, float ampl, float fwhm) { 
  // ============================================================
  //
  //  constructor 
  //  
  //  The procedure is the following: 
  //
  //  1. Allocation of some memory needed
  //  2. some parameters of the trigger are set.   
  //  3. Then the all signals are set to zero

  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 = fwhm       ; 
  ampl_resp = ampl  ; 

  overlaping_time = overt;

    
  threshold = CHANNEL_THRESHOLD  ; 
  

  gate_leng        = gate       ; 
  trigger_multi    = TRIGGER_MULTI      ; 
  trigger_geometry = TRIGGER_GEOM ; 

  cout << endl
       << "[MTrigger]  Setting up the MTrigger with this values "<< endl ; 
  cout << "[MTrigger]    Gate Length:        " << gate_leng  << " ns"
       << endl ; 
  cout << "[MTrigger]    Overlaping time:    " << overlaping_time << " ns"
       << endl ; 
  cout << "[MTrigger]    Response FWHM:      " << fwhm_resp  << " ns"
       << endl ; 
  cout << "[MTrigger]    Response Amplitude: " << ampl_resp  << " mV"
       << endl ; 
  cout << endl ; 


  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) ;
  }


  //
  //  Read in the lookup table for trigger cells
  //
  
  i = 0 ; 

  if ( (unit = fopen("../include-MTrigger/TABLE_PIXELS_IN_CELLS", "r" )) == 0 ) { 
    cout << "ERROR: not able to read ../include-MTrigger/TABLE_PIXELS_IN_CELLS"
         << endl ; 
    exit(123) ; 
  }  
  else { 
    while ( i < TRIGGER_PIXELS )
      { 
        for ( Int_t k=0; k<TRIGGER_CELLS; k++ ) { 
          TC[k][i]=FALSE; 
        }
        i++ ;
      }
    while ( feof(unit) == 0 ) {
        for ( Int_t k=0; k<TRIGGER_CELLS; k++ ) { 
          fscanf ( unit, "%d ", &i ) ; 
          if((i-1)<TRIGGER_PIXELS)
            TC[k][i-1]=TRUE;
        }
    }
    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 ; 
    SlicesSecond[i] = -50 ;
    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 ; 
    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::GetResponse(Float_t *resp) {
  // ============================================================
  //
  //    puts the standard response function into the array resp

  for ( Int_t i=0; i< RESPONSE_SLICES; i++ ) {

    resp[i] = sing_resp[i] ; 
  }
  
}

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]){
        //
        //  Loop over trigger cells. It is topology analisy,
        //  therefore it is keep here after multiplicity and 
        //  threshold checks.
        //

        for(Int_t iCell=0; iCell<TRIGGER_CELLS; iCell++){
          //
          //  then look in all pixel of that cell  if the 
          //  diskriminated signal is 1
          //
          for ( Int_t iPix = 0 ; iPix < TRIGGER_PIXELS; iPix++ ) {
            Muster[iPix] = kFALSE ; 
            Neighb[iPix] = kFALSE ;
            // Select pixels which are used and it the current cell
            if ( used [iPix] == TRUE && TC[iCell][iPix]==TRUE) {
              //
              //  now check the diskriminated signal
              //
              if ( d_sig [iPix][iSli] > 0. ) {
                Muster[iPix] = kTRUE ;
              } 
            }
          } // end of loop over the pixels

          //
          //  Here we check which of the "muster" pixels will be fired for
          //  the minimum required overlaping time
          //

          OverlapingTime(Muster, &Muster[0],iSli);

          //
          //   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++){
              
              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)
                  iCell=TRIGGER_CELLS;         // We skip the remaining trigger cells 
                  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)
                  iCell=TRIGGER_CELLS;         // We skip the remaining trigger cells 
                  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
            // The algorithm is not perfect, there still somes cases
            // that are not really well treated

            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 ;
                iMulti=1;

                //while(PassNextNeighbour(Muster, &Neighb[0])) iMulti++;
                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
                                      
                      for (int jj=0;jj<TRIGGER_PIXELS;jj++) {

                        Aux1[jj] = Neighb[jj] ;   // huschel
                        Aux2[jj]=kFALSE;
                      }

                      //  It checks if taking out any of the pixels we lose
                      //  neighbouring condition for trigger_multi -1
                      
                      Aux1[i]=kFALSE;
                      closed_pack=0;
                      for (int jj=0;jj<TRIGGER_PIXELS;jj++) {
                        if (Aux1[jj]==kTRUE){
                          Aux2[jj]=kTRUE;
                          for (iMulti=1;iMulti<(trigger_multi-1);iMulti++){
                            if (!PassNextNeighbour(Aux1, &Aux2[0]))
                              break; 
                          }
                          if (iMulti==(trigger_multi-1)){
                            //  We found a NN trigger for trigger_multi -1
                            //  taking out pixel jj
                            closed_pack=1;
                            break;
                          }
                          Aux2[jj]=kFALSE;
                        }
                      }
                      if (!closed_pack) break;
                      // For some pixell we did not found NN condition
                      // for trigger_multi -1
                    }
                  }
                  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)
                    iCell=TRIGGER_CELLS;         // We skip the remaining trigger cells 
                    break ; 
                  }
                  else {
                    for (int k=0; k<TRIGGER_PIXELS; k++){
                      if (Neighb[k]){
                        Neighb[k]=kFALSE;
                      }
                    }
                  }
                }  //  end if trigger multiplicity achived
                else{
                  for (int k=0; k<TRIGGER_PIXELS; k++)
                    Neighb[k]=kFALSE;
                }                 
              } // end if pixel fired
            } // end loop trigger pixels
            break;
          };  // end case 2:
          default:{
            cout << "This topology is not implemented yet"<<endl;
            break;
          }
          }
        } //end loop over trigger cells.
      }
    } // 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]);
}

void MTrigger::OverlapingTime ( Bool_t m[], Bool_t *n, Int_t ifSli){

  //============================================================
  //
  //  It returns in n the pixels of m that are fired during the
  //  required overlaping time for trigger after ifSli

  int i,j;
  int iNumSli;

  // Translation from ns to slices
  iNumSli=(int) (overlaping_time*SLICES_PER_NSEC);
  if (iNumSli<1) iNumSli=1;

  //  Put pixels that fulfill the requirement in n
  for (i=0;i<TRIGGER_PIXELS;i++){
    if (m[i]==kTRUE){
      for(j=ifSli;j<ifSli+iNumSli;j++){
        if(!d_sig[i][j]){
          n[i]=kFALSE;
          break;
        }
      }
    }
  }

}