remove exec
[u/mrichter/AliRoot.git] / EMCAL / AliCaloRawAnalyzer.cxx
index 38b317b..bdec18b 100644 (file)
@@ -1,3 +1,4 @@
+// -*- mode: c++ -*-
 /**************************************************************************
  * This file is property of and copyright by                              *
  * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009     *
 #include <iostream>
 using namespace std;
 
-//ClassImp(AliCaloRawAnalyzer)  
+ClassImp(AliCaloRawAnalyzer)  
 
 AliCaloRawAnalyzer::AliCaloRawAnalyzer(const char *name, const char *nameshort) :  TObject(),
-                                                                                  fMinTimeIndex(-1),
-                                                                                  fMaxTimeIndex(-1),
-                                                                                  fFitArrayCut(5),
-                                                                                  fAmpCut(4),
-                                                                                  fNsampleCut(5),
-                                                                                  fNsamplePed(3),
-                                                                                  fIsZerosupressed( false ),
-                                                                                  fVerbose( false )
+  fMinTimeIndex(-1),
+  fMaxTimeIndex(-1),
+  fFitArrayCut(5),
+  fAmpCut(4),
+  fNsampleCut(5),
+  fOverflowCut(950),
+  fNsamplePed(3),
+  fIsZerosupressed( false ),
+  fVerbose( false ),
+  fAlgo(Algo::kNONE), 
+  fL1Phase(0),
+  fAmp(0),
+  fTof(0),
+  fTau( EMCAL::TAU )
 {
   //Comment 
-  sprintf(fName, "%s", name);
-  sprintf(fNameShort, "%s", nameshort);
-    
-  for(int i=0; i < MAXSAMPLES; i++ )
+  snprintf(fName, 256,"%s", name);
+  snprintf(fNameShort,256, "%s", nameshort);
+  //  sprintf(fName ,"%s", name);
+  // sprintf(fNameShort, "%s", nameshort);
+  
+  for(int i=0; i < ALTROMAXSAMPLES; i++ )
     {
       fReversed[i] = 0;
     }
@@ -64,7 +73,7 @@ void
 AliCaloRawAnalyzer::SetTimeConstraint(const int min, const int max ) 
 {
   //Require that the bin if the maximum ADC value is between min and max (timebin)
-  if(  ( min > max ) || min > MAXSAMPLES  || max > MAXSAMPLES  )
+  if(  ( min > max ) || min > ALTROMAXSAMPLES  || max > ALTROMAXSAMPLES  )
     {
       AliWarning( Form( "Attempt to set Invalid time bin range (Min , Max) = (%d, %d), Ingored",  min, max ) ); 
     }
@@ -94,35 +103,50 @@ AliCaloRawAnalyzer::Max(const UShort_t *data, const int length ) const
 
 
 void 
-AliCaloRawAnalyzer::SelectSubarray( const Double_t *fData, const int length, const short maxindex, int *const first,  int *const last ) const
+AliCaloRawAnalyzer::SelectSubarray( const Double_t *data, const int length, const short maxindex,int *const first,  int *const last, const int cut) const
 {
   //Selection of subset of data from one bunch that will be used for fitting or
   //Peak finding.  Go to the left and right of index of the maximum time bin
-  //Untile the ADC value is less that fFitArrayCut
+  //Until the ADC value is less that fFitArrayCut, or derivative changes sign (data jump)
   int tmpfirst =  maxindex;
   int tmplast  =  maxindex;
-  
-  while((  tmpfirst  ) > 0  &&  ( fData[tmpfirst] >  fFitArrayCut   ))  
+  Double_t prevFirst =  data[maxindex];
+  Double_t prevLast  =  data[maxindex];  
+  bool firstJump = false;
+  bool lastJump = false;
+
+  while( (tmpfirst >= 0) && (data[tmpfirst] >= cut ) && (!firstJump) ) 
     {
+      // jump check:
+      if (tmpfirst != maxindex) { // neighbor to maxindex can share peak with maxindex
+       if ( data[tmpfirst] >= prevFirst) {
+         firstJump = true;
+       }
+      }
+      prevFirst = data[tmpfirst];
       tmpfirst -- ;
     }
   
-  while(( tmplast ) <  (length-1)   && ( fData [tmplast] >  fFitArrayCut ))
+  while( (tmplast < length) && (data[tmplast] >=  cut ) && (!lastJump) ) 
     {
+      // jump check:
+      if (tmplast != maxindex) { // neighbor to maxindex can share peak with maxindex
+       if ( data[tmplast] >= prevLast) {
+         lastJump = true;
+       }
+      }
+      prevLast = data[tmplast];
       tmplast ++;
     }
 
+  // we keep one pre- or post- sample if we can (as in online)
+  // check though if we ended up on a 'jump', or out of bounds: if so, back up
+  if (firstJump || tmpfirst<0) tmpfirst ++;
+  if (lastJump || tmplast>=length) tmplast --;
 
-  *first = tmpfirst +1;
-  *last =  tmplast -1;
-
-  int nsamples = *last - *first + 1;
-  if (nsamples < fNsampleCut) {
-    // keep edge bins (below threshold) also, for low # of samples case    
-    *first = tmpfirst;
-    *last =  tmplast; 
-  }
-
+  *first = tmpfirst;
+  *last =  tmplast;
+  return;
 }
 
 
@@ -220,7 +244,7 @@ AliCaloRawAnalyzer::CheckBunchEdgesForMax( const AliCaloBunchInfo *const bunch )
 
 
 int 
-AliCaloRawAnalyzer::SelectBunch( const vector<AliCaloBunchInfo> &bunchvector,short *const maxampbin, short *const maxamplitude ) const
+AliCaloRawAnalyzer::SelectBunch( const vector<AliCaloBunchInfo> &bunchvector, short *const maxampbin, short *const maxamplitude )
 {
   //We select the bunch with the highest amplitude unless any time constraints is set
   short max = -1;
@@ -231,7 +255,7 @@ AliCaloRawAnalyzer::SelectBunch( const vector<AliCaloBunchInfo> &bunchvector,sho
   for(unsigned int i=0; i < bunchvector.size(); i++ )
     {
       max = Max(  &bunchvector.at(i), &indx ); // CRAP PTH, bug fix, trouble if more than one bunches  
-      if( IsInTimeRange( indx) )
+      if( IsInTimeRange( indx, fMaxTimeIndex, fMinTimeIndex) )
        {
          if( max > maxall )
            {
@@ -255,19 +279,19 @@ AliCaloRawAnalyzer::SelectBunch( const vector<AliCaloBunchInfo> &bunchvector,sho
 
 
 bool 
-AliCaloRawAnalyzer::IsInTimeRange( const int maxindex  ) const
+AliCaloRawAnalyzer::IsInTimeRange( const int maxindex, const int maxtindx, const int mintindx ) const
 {
   // Ckeck if the index of the max ADC vaue is consistent with trigger.
-  if( ( fMinTimeIndex  < 0 && fMaxTimeIndex  < 0) ||fMaxTimeIndex  < 0 )
+  if( ( mintindx  < 0 && maxtindx   < 0) ||maxtindx  < 0 )
     {
       return true; 
     }
-  return ( maxindex < fMaxTimeIndex ) && ( maxindex > fMinTimeIndex  ) ? true : false;
+  return ( maxindex < maxtindx ) && ( maxindex > mintindx  ) ? true : false;
 }
 
 
 void 
-AliCaloRawAnalyzer::PrintBunches( const vector<AliCaloBunchInfo> &bvctr ) const
+AliCaloRawAnalyzer::PrintBunches( const vector<AliCaloBunchInfo> &bvctr ) 
 {
   //comment
   cout << __FILE__ << __LINE__<< "*************** Printing Bunches *******************" << endl;
@@ -283,7 +307,7 @@ AliCaloRawAnalyzer::PrintBunches( const vector<AliCaloBunchInfo> &bvctr ) const
 
 
 void 
-AliCaloRawAnalyzer::PrintBunch( const AliCaloBunchInfo &bunch ) const
+AliCaloRawAnalyzer::PrintBunch( const AliCaloBunchInfo &bunch )
 {
   //comment
   cout << __FILE__ << ":" << __LINE__ << endl;
@@ -298,40 +322,124 @@ AliCaloRawAnalyzer::PrintBunch( const AliCaloBunchInfo &bunch ) const
 }
 
 
-AliCaloFitResults
-AliCaloRawAnalyzer::Evaluate( const vector<AliCaloBunchInfo>  &/*bunchvector*/, const UInt_t /*altrocfg1*/,  const UInt_t /*altrocfg2*/)
-{ // method to do the selection of what should possibly be fitted
-  // not implemented for base class
-  return AliCaloFitResults( 0, 0 );
+Double_t
+AliCaloRawAnalyzer::CalculateChi2(const Double_t amp, const Double_t time,
+                                 const Int_t first, const Int_t last,
+                                 const Double_t adcErr, 
+                                 const Double_t tau) const
+{
+  //   Input:
+  //   amp   - max amplitude;
+  //   time    - time of max amplitude; 
+  //   first, last - sample array indices to be used
+  //   adcErr   - nominal error of amplitude measurement (one value for all channels)
+  //           if adcErr<0 that mean adcErr=1.
+  //   tau   - filter time response (in timebin units)
+  // Output:
+  //   chi2 - chi2
+
+  if (first == last || first<0 ) { // signal consists of single sample, chi2 estimate (0) not too well defined.. 
+    // or, first is negative, the indices are not valid
+    return Ret::kDummy;
+  }
+
+  int nsamples =  last - first + 1;
+  // printf(" AliCaloRawAnalyzer::CalculateChi2 : first %i last %i : nsamples %i : amp %3.2f time %3.2f \n", first, last, nsamples, amp, time); 
+
+  Int_t x = 0;
+  Double_t chi2 = 0;
+  Double_t dy = 0.0, xx = 0.0, f=0.0;
+
+  for (Int_t i=0; i<nsamples; i++) {
+    x     = first + i; // timebin
+    xx    = (x - time + tau) / tau; // help variable
+    f     = 0;
+    if (xx > 0) {
+      f = amp * xx*xx * TMath::Exp(2 * (1 - xx )) ;
+    }
+    dy    = fReversed[x] - f; 
+    chi2 += dy*dy;
+    // printf(" AliCaloRawAnalyzer::CalculateChi2 : %i : y %f -> f %f : dy %f \n", i, fReversed[first+i], f, dy); 
+  }
+
+  if (adcErr>0.0) { // weight chi2
+    chi2 /= (adcErr*adcErr);
+  }
+  return chi2;
+}
+
+
+void
+AliCaloRawAnalyzer::CalculateMeanAndRMS(const Int_t first, const Int_t last,
+                                       Double_t & mean, Double_t & rms)
+{
+  //   Input:
+  //   first, last - sample array indices to be used
+  // Output:
+  //   mean and RMS of samples 
+  //
+  // To possibly be used to differentiate good signals from bad before fitting
+  // 
+  mean = Ret::kDummy;
+  rms =  Ret::kDummy;
+
+  if (first == last || first<0 ) { // signal consists of single sample, chi2 estimate (0) not too well defined.. 
+    // or, first is negative, the indices are not valid
+    return;
+  }
+
+  int nsamples =  last - first + 1;
+  //  printf(" AliCaloRawAnalyzer::CalculateMeanAndRMS : first %i last %i : nsamples %i \n", first, last, nsamples); 
+
+  int x = 0;
+  Double_t sampleSum = 0; // sum of samples
+  Double_t squaredSampleSum = 0; // sum of samples squared
+
+  for (Int_t i=0; i<nsamples; i++) {
+    x = first + i;
+    sampleSum += fReversed[x];
+    squaredSampleSum += (fReversed[x] * fReversed[x]);
+  }
+
+  mean = sampleSum / nsamples;          
+  Double_t squaredMean = squaredSampleSum / nsamples;   
+  // The variance (rms squared) is equal to the mean of the squares minus the square of the mean..      
+  rms = sqrt(squaredMean - mean*mean);
+
+  return;
 }
 
 
+
 int
-AliCaloRawAnalyzer::PreFitEvaluateSamples( const vector<AliCaloBunchInfo>  &bunchvector, const UInt_t altrocfg1,  const UInt_t altrocfg2, Int_t & index, Float_t & maxf, short & maxamp, short & maxrev, Float_t & ped, int & first, int & last)
+AliCaloRawAnalyzer::PreFitEvaluateSamples( const vector<AliCaloBunchInfo>  &bunchvector, const UInt_t altrocfg1,  
+                                          const UInt_t altrocfg2, Int_t & index, Float_t & maxf, short & maxamp, 
+                                          short & maxrev, Float_t & ped, int & first, int & last,const int acut )
 { // method to do the selection of what should possibly be fitted
   int nsamples  = 0;
   short maxampindex = 0;
   index = SelectBunch( bunchvector,  &maxampindex,  &maxamp ); // select the bunch with the highest amplitude unless any time constraints is set
 
   
-  if( index >= 0 && maxamp > fAmpCut) // something valid was found, and non-zero amplitude
+  if( index >= 0 && maxamp >= acut ) // something valid was found, and non-zero amplitude
     {
       // use more convenient numbering and possibly subtract pedestal
       ped  = ReverseAndSubtractPed( &(bunchvector.at(index)),  altrocfg1, altrocfg2, fReversed  );
       maxf = TMath::MaxElement( bunchvector.at(index).GetLength(),  fReversed );
       
-      if ( maxf > fAmpCut ) // possibly significant signal
+      if ( maxf >= acut  ) // possibly significant signal
        {
          // select array around max to possibly be used in fit
          maxrev = maxampindex - bunchvector.at(index).GetStartBin(); 
-         SelectSubarray( fReversed,  bunchvector.at(index).GetLength(),  maxrev, &first, &last);
+         SelectSubarray( fReversed,  bunchvector.at(index).GetLength(),  maxrev, &first, &last, acut );
 
          // sanity check: maximum should not be in first or last bin
          // if we should do a fit
-         if (first!=maxrev && last!=maxrev) {
-           // calculate how many samples we have 
-           nsamples =  last - first + 1;         
-         }
+         if (first!=maxrev && last!=maxrev) 
+           {
+             // calculate how many samples we have 
+             nsamples =  last - first + 1;       
+           }
        }
     }