Adding include path to allow compilation of CleanGeom task
[u/mrichter/AliRoot.git] / EMCAL / AliCaloRawAnalyzer.cxx
1 // -*- mode: c++ -*-
2 /**************************************************************************
3  * This file is property of and copyright by                              *
4  * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009     *
5  *                                                                        *
6  * Primary Author: Per Thomas Hille <perthomas.hille@yale.edu>            *
7  *                                                                        *
8  * Contributors are mentioned in the code where appropriate.              *
9  * Please report bugs to p.t.hille@fys.uio.no                             *
10  *                                                                        *
11  * Permission to use, copy, modify and distribute this software and its   *
12  * documentation strictly for non-commercial purposes is hereby granted   *
13  * without fee, provided that the above copyright notice appears in all   *
14  * copies and that both the copyright notice and this permission notice   *
15  * appear in the supporting documentation. The authors make no claims     *
16  * about the suitability of this software for any purpose. It is          *
17  * provided "as is" without express or implied warranty.                  *
18  **************************************************************************/
19
20
21 // Base class for extraction 
22 // of signal amplitude and peak position
23 // From CALO Calorimeter RAW data (from the RCU)
24 // Contains some utilities for preparing / selecting
25 // Signals suitable for signal extraction
26 // By derived classes
27
28 #include "AliLog.h"
29 #include "AliCaloRawAnalyzer.h"
30 #include "AliCaloBunchInfo.h"
31 #include "AliCaloFitResults.h"
32 #include "TMath.h"
33 #include <iostream>
34 using namespace std;
35
36 ClassImp(AliCaloRawAnalyzer)  
37
38 AliCaloRawAnalyzer::AliCaloRawAnalyzer(const char *name, const char *nameshort) :  TObject(),
39   fMinTimeIndex(-1),
40   fMaxTimeIndex(-1),
41   fFitArrayCut(5),
42   fAmpCut(4),
43   fNsampleCut(5),
44   fOverflowCut(950),
45   fNsamplePed(3),
46   fIsZerosupressed( false ),
47   fVerbose( false ),
48   fAlgo(Algo::kNONE), 
49   fL1Phase(0),
50   fAmp(0),
51   fTof(0),
52   fTau( EMCAL::TAU )
53 {
54   //Comment 
55   snprintf(fName, 256,"%s", name);
56   snprintf(fNameShort,256, "%s", nameshort);
57   //  sprintf(fName ,"%s", name);
58   // sprintf(fNameShort, "%s", nameshort);
59   
60   for(int i=0; i < ALTROMAXSAMPLES; i++ )
61     {
62       fReversed[i] = 0;
63     }
64 }
65
66 AliCaloRawAnalyzer::~AliCaloRawAnalyzer()
67 {
68
69 }
70
71
72 void 
73 AliCaloRawAnalyzer::SetTimeConstraint(const int min, const int max ) 
74 {
75   //Require that the bin if the maximum ADC value is between min and max (timebin)
76   if(  ( min > max ) || min > ALTROMAXSAMPLES  || max > ALTROMAXSAMPLES  )
77     {
78       AliWarning( Form( "Attempt to set Invalid time bin range (Min , Max) = (%d, %d), Ingored",  min, max ) ); 
79     }
80   else
81     {
82       fMinTimeIndex  = min;
83       fMaxTimeIndex  = max;
84     }
85 }
86
87
88 UShort_t 
89 AliCaloRawAnalyzer::Max(const UShort_t *data, const int length ) const
90 {
91   //------------
92   UShort_t tmpmax  = data[0];
93
94   for(int i=0; i < length; i++)
95     {
96       if( tmpmax  <  data[i] )
97         {
98           tmpmax = data[i];
99         }
100     }
101   return tmpmax;
102 }
103
104
105 void 
106 AliCaloRawAnalyzer::SelectSubarray( const Double_t *data, const int length, const short maxindex,int *const first,  int *const last, const int cut) const
107 {
108   //Selection of subset of data from one bunch that will be used for fitting or
109   //Peak finding.  Go to the left and right of index of the maximum time bin
110   //Until the ADC value is less that fFitArrayCut, or derivative changes sign (data jump)
111   int tmpfirst =  maxindex;
112   int tmplast  =  maxindex;
113   Double_t prevFirst =  data[maxindex];
114   Double_t prevLast  =  data[maxindex];  
115   bool firstJump = false;
116   bool lastJump = false;
117
118   while( (tmpfirst >= 0) && (data[tmpfirst] >= cut ) && (!firstJump) ) 
119     {
120       // jump check:
121       if (tmpfirst != maxindex) { // neighbor to maxindex can share peak with maxindex
122         if ( data[tmpfirst] >= prevFirst) {
123           firstJump = true;
124         }
125       }
126       prevFirst = data[tmpfirst];
127       tmpfirst -- ;
128     }
129   
130   while( (tmplast < length) && (data[tmplast] >=  cut ) && (!lastJump) ) 
131     {
132       // jump check:
133       if (tmplast != maxindex) { // neighbor to maxindex can share peak with maxindex
134         if ( data[tmplast] >= prevLast) {
135           lastJump = true;
136         }
137       }
138       prevLast = data[tmplast];
139       tmplast ++;
140     }
141
142   // we keep one pre- or post- sample if we can (as in online)
143   // check though if we ended up on a 'jump', or out of bounds: if so, back up
144   if (firstJump || tmpfirst<0) tmpfirst ++;
145   if (lastJump || tmplast>=length) tmplast --;
146
147   *first = tmpfirst;
148   *last =  tmplast;
149   return;
150 }
151
152
153
154 Float_t 
155 AliCaloRawAnalyzer::ReverseAndSubtractPed( const AliCaloBunchInfo *bunch, const UInt_t /*altrocfg1*/,  const UInt_t /*altrocfg2*/, double *outarray ) const
156 {
157   //Time sample comes in reversed order, revers them back
158   //Subtract the baseline based on content of altrocfg1 and altrocfg2.
159   Int_t length =  bunch->GetLength(); 
160   const UShort_t *sig = bunch->GetData();
161
162   double ped = EvaluatePedestal( sig , length);
163
164   for( int i=0; i < length; i++ )
165     {
166       outarray[i] = sig[length -i -1] - ped;
167     }
168
169   return ped;
170 }
171
172
173
174 Float_t 
175 AliCaloRawAnalyzer::EvaluatePedestal(const UShort_t * const data, const int /*length*/ ) const
176 {
177   //  double ped = 0;
178   double tmp = 0;
179
180   if( fIsZerosupressed == false ) 
181     {
182       for(int i=0; i < fNsamplePed; i++ )
183         {
184           tmp += data[i];
185         }
186    }
187
188   return  tmp/fNsamplePed;
189  }
190
191
192 short  
193 AliCaloRawAnalyzer::Max( const AliCaloBunchInfo *const bunch , int *const maxindex ) const
194 {
195   //comment
196   short tmpmax = -1;
197   int tmpindex = -1;
198   const UShort_t *sig = bunch->GetData();
199
200   for(int i=0; i < bunch->GetLength(); i++ )
201     {
202       if( sig[i] > tmpmax )
203         {
204           tmpmax   =  sig[i];
205           tmpindex =  i; 
206         }
207     }
208   
209   if(maxindex != 0 )
210     {
211       //   *maxindex =  bunch->GetLength() -1 - tmpindex + bunch->GetStartBin(); 
212        *maxindex =  bunch->GetLength() -1 - tmpindex + bunch->GetStartBin(); 
213     }
214   
215   return  tmpmax;
216 }
217
218
219 bool  
220 AliCaloRawAnalyzer::CheckBunchEdgesForMax( const AliCaloBunchInfo *const bunch ) const
221 {
222   // a bunch is considered invalid if the maximum is in the first or last time-bin
223   short tmpmax = -1;
224   int tmpindex = -1;
225   const UShort_t *sig = bunch->GetData();
226
227   for(int i=0; i < bunch->GetLength(); i++ )
228     {
229       if( sig[i] > tmpmax )
230         {
231           tmpmax   =  sig[i];
232           tmpindex =  i; 
233         }
234     }
235   
236   bool bunchOK = true;
237   if (tmpindex == 0 || tmpindex == (bunch->GetLength()-1) )
238     {
239       bunchOK = false;
240     }
241   
242   return  bunchOK;
243 }
244
245
246 int 
247 AliCaloRawAnalyzer::SelectBunch( const vector<AliCaloBunchInfo> &bunchvector, short *const maxampbin, short *const maxamplitude )
248 {
249   //We select the bunch with the highest amplitude unless any time constraints is set
250   short max = -1;
251   short bunchindex = -1;
252   short maxall = -1;
253   int indx = -1;
254
255   for(unsigned int i=0; i < bunchvector.size(); i++ )
256     {
257       max = Max(  &bunchvector.at(i), &indx ); // CRAP PTH, bug fix, trouble if more than one bunches  
258       if( IsInTimeRange( indx, fMaxTimeIndex, fMinTimeIndex) )
259         {
260           if( max > maxall )
261             {
262               maxall = max;
263               bunchindex = i;
264               *maxampbin     = indx;
265               *maxamplitude  = max;
266             }
267         }
268     }
269  
270   if (bunchindex >= 0) {
271     bool bunchOK = CheckBunchEdgesForMax( &bunchvector.at(bunchindex) );
272     if (! bunchOK) { 
273       bunchindex = -1; 
274     }
275   }
276
277   return  bunchindex;
278 }
279
280
281 bool 
282 AliCaloRawAnalyzer::IsInTimeRange( const int maxindex, const int maxtindx, const int mintindx ) const
283 {
284   // Ckeck if the index of the max ADC vaue is consistent with trigger.
285   if( ( mintindx  < 0 && maxtindx   < 0) ||maxtindx  < 0 )
286     {
287       return true; 
288     }
289   return ( maxindex < maxtindx ) && ( maxindex > mintindx  ) ? true : false;
290 }
291
292
293 void 
294 AliCaloRawAnalyzer::PrintBunches( const vector<AliCaloBunchInfo> &bvctr ) 
295 {
296   //comment
297   cout << __FILE__ << __LINE__<< "*************** Printing Bunches *******************" << endl;
298   cout << __FILE__ << __LINE__<< "*** There are " << bvctr.size() << ", bunches" << endl;
299
300   for(unsigned int i=0; i < bvctr.size() ; i++ )
301     {
302       PrintBunch( bvctr.at(i) );
303       cout << " bunch = "  <<  i  << endl;
304     }
305   cout << __FILE__ << __LINE__<< "*************** Done ... *******************" << endl;
306 }
307
308
309 void 
310 AliCaloRawAnalyzer::PrintBunch( const AliCaloBunchInfo &bunch )
311 {
312   //comment
313   cout << __FILE__ << ":" << __LINE__ << endl;
314   cout << __FILE__ << __LINE__   << ", startimebin = " << bunch.GetStartBin() << ", length =" <<  bunch.GetLength() << endl;
315   const UShort_t *sig =  bunch.GetData();  
316   
317   for ( Int_t j = 0; j <  bunch.GetLength();  j++) 
318     {
319       printf("%d\t", sig[j] );
320     }
321   cout << endl; 
322 }
323
324
325 Double_t
326 AliCaloRawAnalyzer::CalculateChi2(const Double_t amp, const Double_t time,
327                                   const Int_t first, const Int_t last,
328                                   const Double_t adcErr, 
329                                   const Double_t tau) const
330 {
331   //   Input:
332   //   amp   - max amplitude;
333   //   time    - time of max amplitude; 
334   //   first, last - sample array indices to be used
335   //   adcErr   - nominal error of amplitude measurement (one value for all channels)
336   //           if adcErr<0 that mean adcErr=1.
337   //   tau   - filter time response (in timebin units)
338   // Output:
339   //   chi2 - chi2
340
341   if (first == last || first<0 ) { // signal consists of single sample, chi2 estimate (0) not too well defined.. 
342     // or, first is negative, the indices are not valid
343     return Ret::kDummy;
344   }
345
346   int nsamples =  last - first + 1;
347   // printf(" AliCaloRawAnalyzer::CalculateChi2 : first %i last %i : nsamples %i : amp %3.2f time %3.2f \n", first, last, nsamples, amp, time); 
348
349   Int_t x = 0;
350   Double_t chi2 = 0;
351   Double_t dy = 0.0, xx = 0.0, f=0.0;
352
353   for (Int_t i=0; i<nsamples; i++) {
354     x     = first + i; // timebin
355     xx    = (x - time + tau) / tau; // help variable
356     f     = 0;
357     if (xx > 0) {
358       f = amp * xx*xx * TMath::Exp(2 * (1 - xx )) ;
359     }
360     dy    = fReversed[x] - f; 
361     chi2 += dy*dy;
362     // printf(" AliCaloRawAnalyzer::CalculateChi2 : %i : y %f -> f %f : dy %f \n", i, fReversed[first+i], f, dy); 
363   }
364
365   if (adcErr>0.0) { // weight chi2
366     chi2 /= (adcErr*adcErr);
367   }
368   return chi2;
369 }
370
371
372 void
373 AliCaloRawAnalyzer::CalculateMeanAndRMS(const Int_t first, const Int_t last,
374                                         Double_t & mean, Double_t & rms)
375 {
376   //   Input:
377   //   first, last - sample array indices to be used
378   // Output:
379   //   mean and RMS of samples 
380   //
381   // To possibly be used to differentiate good signals from bad before fitting
382   // 
383   mean = Ret::kDummy;
384   rms =  Ret::kDummy;
385
386   if (first == last || first<0 ) { // signal consists of single sample, chi2 estimate (0) not too well defined.. 
387     // or, first is negative, the indices are not valid
388     return;
389   }
390
391   int nsamples =  last - first + 1;
392   //  printf(" AliCaloRawAnalyzer::CalculateMeanAndRMS : first %i last %i : nsamples %i \n", first, last, nsamples); 
393
394   int x = 0;
395   Double_t sampleSum = 0; // sum of samples
396   Double_t squaredSampleSum = 0; // sum of samples squared
397
398   for (Int_t i=0; i<nsamples; i++) {
399     x = first + i;
400     sampleSum += fReversed[x];
401     squaredSampleSum += (fReversed[x] * fReversed[x]);
402   }
403
404   mean = sampleSum / nsamples;   
405   Double_t squaredMean = squaredSampleSum / nsamples;    
406   // The variance (rms squared) is equal to the mean of the squares minus the square of the mean..       
407   rms = sqrt(squaredMean - mean*mean);
408
409   return;
410 }
411
412
413
414 int
415 AliCaloRawAnalyzer::PreFitEvaluateSamples( const vector<AliCaloBunchInfo>  &bunchvector, const UInt_t altrocfg1,  
416                                            const UInt_t altrocfg2, Int_t & index, Float_t & maxf, short & maxamp, 
417                                            short & maxrev, Float_t & ped, int & first, int & last,const int acut )
418 { // method to do the selection of what should possibly be fitted
419   int nsamples  = 0;
420   short maxampindex = 0;
421   index = SelectBunch( bunchvector,  &maxampindex,  &maxamp ); // select the bunch with the highest amplitude unless any time constraints is set
422
423   
424   if( index >= 0 && maxamp >= acut ) // something valid was found, and non-zero amplitude
425     {
426       // use more convenient numbering and possibly subtract pedestal
427       ped  = ReverseAndSubtractPed( &(bunchvector.at(index)),  altrocfg1, altrocfg2, fReversed  );
428       maxf = TMath::MaxElement( bunchvector.at(index).GetLength(),  fReversed );
429       
430       if ( maxf >= acut  ) // possibly significant signal
431         {
432           // select array around max to possibly be used in fit
433           maxrev = maxampindex - bunchvector.at(index).GetStartBin(); 
434           SelectSubarray( fReversed,  bunchvector.at(index).GetLength(),  maxrev, &first, &last, acut );
435
436           // sanity check: maximum should not be in first or last bin
437           // if we should do a fit
438           if (first!=maxrev && last!=maxrev) 
439             {
440               // calculate how many samples we have 
441               nsamples =  last - first + 1;       
442             }
443         }
444     }
445
446   return nsamples;
447 }
448