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

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

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


MTrigger::MTrigger() { 
  
  FILE *unit_mtrig ; 
  Int_t endflag = 1  ; 
  char   datac[256] ; 
  char   dummy[50] ; 
  //
  //  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
  
  // 
  //   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  ; 
  
  chan_thres    = CHANNEL_THRESHOLD  ; 
  gate_leng     = TRIGGER_GATE       ; 
  trigger_multi = TRIGGER_MULTI      ; 
 
  //
  //  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, &chan_thres ) ; 
      }
      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 ) ;
      }

      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 ; 
  cout << endl 
       << "[MTrigger]    ChannelThreshold:   " << chan_thres << " 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 << endl ; 
  //
  //    set up the response shape
  // 
  Int_t  i, ii ; 
     
  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) ) ; 

    //      cout << i << "  "
    //           << x << "  "
    //           << sing_resp[i] 
    //           << endl ; 

  } 

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

  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 ;

    nphot[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 ( i = 0 ; i < 5 ; i++) {
    SlicesZero[i]   = 0 ; 
    SlicesFirst[i]  = 0 ; 
    SlicesSecond[i] = 0 ; 
  }
  
  cout << " end of MTrigger::MTrigger()" << endl ; 
} 

MTrigger::~MTrigger() {

  delete histPmt ; 
}
  
void MTrigger::Reset() { 
  //
  //  set all values of the signals to zero
  //
  Int_t  i, ii ; 
  
  for ( i =0 ; i <TRIGGER_PIXELS ; i++ ) {
    used [i] = FALSE ; 
    dknt [i] = FALSE ; 
    
    nphot[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 ( i = 0 ; i < 5 ; i++) {
    SlicesZero[i]   = 0 ; 
    SlicesFirst[i]  = 0 ; 
    SlicesSecond[i] = 0 ; 
  }
    
}


Float_t  MTrigger::Fill( Int_t iPix, Float_t time ) { 
  //
  // fills the information about one single Phe in the Trigger class
  //
  // parameter is the number of the pixel and the time-difference to the
  // first particle
  //
  //

  Int_t i, ichan ; 

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

  //
  //   first we have to check if the pixel iPix is used or not until now
  //   if this is the first use, reset all signal for that pixels
  // 
  if ( iPix >= CAMERA_PIXELS ) {
    // 
    //   the PixelID is greater than the CAMERA  
    //
    cout << " WARNING:  MTrigger::Fill() :  iPix greater than CAMERA_PIXELS"
         << endl ; 
    NoiseAmp = 0. ; 
  }
  else if ( iPix >=TRIGGER_PIXELS ) {
    //
    //   the pixel is inside the Camera, but outside of the TRIGGER-FIELD
    //  
    //   we just scramble the amplitude of the PMT-signal for the FADC
    //
    //   scramble the Amplitude of this single photo electron signal
    //
    NoiseAmp = (histPmt->GetRandom()/histMean) ; 
  } 
  else { 
    //
    //   the photoelectron is contributing to the trigger
    //
    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. ; 
      }
    }
    
    //
    //   then select the time slice to use (ican) 
    // 

    
    if ( time < 0. ) {
      cout << "  WARNING!!   " << time << "  below ZERO!! Very strange!!" 
           << endl ; 
    }
    else if ( time < TOTAL_TRIGGER_TIME ) { 
      nphot[iPix]++ ; 
      //
      ichan = (Int_t) ( time * ((Float_t) SLICES_PER_NSEC) ) ; 
      
      //
      //   scramble the Amplitude of this single photo electron signal
      //
      NoiseAmp = (histPmt->GetRandom()/histMean) ; 
      
      for ( i = 0 ; i<RESPONSE_SLICES; i++ ) {
        
        if ( (ichan+i) < TRIGGER_TIME_SLICES ) {  
          a_sig[iPix][ichan+i] += NoiseAmp * sing_resp[i] ;
 
        } 
      }
    }
    else {
      cout << "  WARNING!!   " << time << "  out of TriggerTimeRange " 
           << TOTAL_TRIGGER_TIME << endl ; 
    }
  }
  
  return NoiseAmp ; 
  
}


Float_t  MTrigger::FillNSB( Int_t iPix, Float_t time ) { 
  //
  // fills the information about one single Phe in the Trigger class
  //
  // parameter is the number of the pixel and the time-difference to the
  // first particle
  //
  //

  Int_t i, ichan ; 

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

  //
  //   first we have to check if the pixel iPix is used or not until now
  //   if this is the first use, reset all signal for that pixels
  // 
  if ( iPix >= CAMERA_PIXELS ) {
    cout << " WARNING:  MTrigger::Fill() :  iPix greater than CAMERA_PIXELS"
         << endl ; 
  }
  else if ( iPix >= TRIGGER_PIXELS ) {
    //
    //   scramble the Amplitude of this single photo electron signal
    //
    NoiseAmp = (histPmt->GetRandom()/histMean) ; 
  } 
  
  else {
    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. ; 
      }
    }

    //
    //   then select the time slice to use (ican) 
    // 
    
    if ( time < 0. ) {
      cout << "  WARNING!!   " << time << "  below ZERO!! Very strange!!" 
           << endl ; 
    }
    else if ( time < TOTAL_TRIGGER_TIME ) { 
      //
      //  FillNSB doesn't add a photon to nphot[iPix] as the method Fill do!!
      //
      
      ichan = (Int_t) ( time * ((Float_t) SLICES_PER_NSEC) ) ; 
      
      //
      //   scramble the Amplitude of this single photo electron signal
      //
      NoiseAmp = (histPmt->GetRandom()/histMean) ; 
      
      for ( i = 0 ; i<RESPONSE_SLICES; i++ ) {
        
        if ( (ichan+i) < TRIGGER_TIME_SLICES ) {  
          a_sig[iPix][ichan+i] += NoiseAmp * sing_resp[i] ; 
        } 
      }
    }
    else {
      cout << "  WARNING!!   " << time << "  out of TriggerTimeRange " 
           << TOTAL_TRIGGER_TIME << endl ; 
    }
  }
  
  return NoiseAmp ; 
  
}

void MTrigger::ElecNoise() {

  Float_t rausch ; 

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

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

      }
    }
  }

}


Int_t MTrigger::Diskriminate() {

  cout << " MTrigger::Diskriminate()" << flush ; 

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

      //cout << "Pixel " << i 
      //   << " baseline " << baseline[i] 
      //   <<endl ; 

    }
  }
  
  //
  //  take only that pixel which are used
  //

  for ( i=0 ; i < TRIGGER_PIXELS; i++  ) {
    if ( used [i] == TRUE ) {
      //cout << "Pixel " << i << " used"  ; 
       
      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  && 
             a_sig[i][ii]   >=  chan_thres  ) { 
          { 
            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. ;  
                sum_d_sig [ii+j] += 1. ;  
              } 
            } 
            ii = ii + jmax ; 
          }
        }
      }
      //      cout << endl ; 
    }
  }

  //cout << "**  MTrigger::Diskriminate()  " << iM << endl ; 


  //
  //   determine the number of zero level triggers
  //
  //   zero level trigger = the sum of all diskriminated signals
  //   is above the TRIGGER_MULTI value. 
  //   only for this events it is neccessay to look for next neighbours!!!
  // 
  
  if ( iM > TRIGGER_MULTI ) {
    Int_t iReturn = 0 ; 

    for ( ii=1 ; ii<TRIGGER_TIME_SLICES; ii++ ) { 
      if ( sum_d_sig[ii] > TRIGGER_MULTI ) { 
        iReturn++ ; 

        SlicesZero[nZero++] = ii ; 
        
        //
        //    if a trigger occurs we read out the next 50 nsec
        //
        //    -> don't study the next 50/0.25 = 200 slices
        //
        ii = ii + 200 ; 
      }      
    } 
    
    return ( iReturn ) ; 
  }
  else {
    return ( 0 ) ; 
  }

  return ( 0 ) ; 
  
}  

Int_t MTrigger::FirstLevel() {

  Int_t iReturn = 0 ; 

  Bool_t Muster[TRIGGER_PIXELS] ; 
  Int_t  iMulti = 0 ; 

  // cout << "#### MTrigger::FirstLevel()" << endl ; 
  // cout << nZero << "  " << SlicesZero[0] <<  endl ; 
  
  if ( nZero > 1 ) { 
    cout << " INFORMATION:  more than one Zero Level TRIGGER " << endl ; 
  }

  //
  //   loop over all ZeroLevel Trigger
  //
  //   it is only neccessary to look after a ZeroLevel Trigger for
  //   a FirstLevel (NextNeighbour) trigger. 
  //
  
  for (Int_t iloop = 0; iloop < nZero ; iloop++ ) {
    
    //
    //   Then run over all slices 
    //   start at the ZeroLevelTrigger slices
    //
    
    for ( Int_t iSli = SlicesZero[iloop];
          iSli < SlicesZero[iloop]+ 200; iSli++ ) {

      //
      //  then look in all pixel if the diskriminated signal is 1
      //
      iMulti = 0 ; 

      for ( Int_t iPix = 0 ; iPix < TRIGGER_PIXELS; iPix++ ) {
        Muster[iPix] = kFALSE ; 

        if ( used [iPix] == TRUE ) {
          //
          //  now check the diskriminated signal
          //
          if ( d_sig [iPix][iSli] > 0. ) {

            iMulti++ ; 
            Muster[iPix] = kTRUE ; 
          } 
        }
      } // end of loop over the pixels
      
      //
      //   here we have to look for next neighbours
      //
      
      if ( PassNextNeighbour ( Muster ) ) { 
        //
        //   A NN-Trigger is detected at time Slice 
        //
        SlicesFirst[nFirst++] = iSli ; 
        iReturn++ ;
        break ; 
      }   
    } // end of loop over the slices  
    
  } // end of loop over zerolevelTriggers

  //
  //   return the Number of FirstLevel Triggers
  //
  return iReturn ; 
}


Bool_t MTrigger::PassNextNeighbour ( Bool_t m[] ) {
  //
  //  This method is looking for next neighbour triggers using a 
  //  NNlookup table. This table is builded by the default constructor
  //
  
  Int_t iNN ; 

  //
  //   loop over all trigger pixels
  //
  for ( Int_t i=0; i<TRIGGER_PIXELS; i++) { 
    //
    //  check if this pixel has a diskrminator signal 
    //  (this is inside m[] ) 
    //

    if ( m[i] ) { 
      iNN = 1 ; 
      //      cout << "/ " << i  ; 
      
      //
      //  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 (NN[i][kk] >= TRIGGER_PIXELS ) {

        } 
        // the nextneighbout is inside the TRIGGER_PIXELS
        else { 
          //
          //  look if the boolean of nn pixels is true
          //
          
          if ( m[ NN[i][kk] ] ) { 
            iNN++ ; 
          }  
        }
      } 
      
      //   cout << "   NN  " << iNN ; 
      
      if ( iNN >=4 ) { 
        return ( kTRUE ) ; 
      }  
    } 
  } 
  return ( kFALSE ) ; 
}

Float_t MTrigger::GetFirstLevelTime(Int_t il ) {
  return ( (Float_t)SlicesFirst[il]/ SLICES_PER_NSEC   ) ; 
}



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