Factorization of the different raw fitting algorithms in EMCAL (Per Thomas)
[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 // fFp(0), 
50   fL1Phase(0),
51   fAmp(0),
52   fTof(0)
53 {
54   //Comment 
55   snprintf(fName, 256,"%s", name);
56   snprintf(fNameShort,256, "%s", nameshort);
57     
58   for(int i=0; i < ALTROMAXSAMPLES; i++ )
59     {
60       fReversed[i] = 0;
61     }
62
63   // fFp = fopen("amp2.txt", "w");
64
65 }
66
67 AliCaloRawAnalyzer::~AliCaloRawAnalyzer()
68 {
69
70 }
71
72
73 void 
74 AliCaloRawAnalyzer::SetTimeConstraint(const int min, const int max ) 
75 {
76   //Require that the bin if the maximum ADC value is between min and max (timebin)
77   if(  ( min > max ) || min > ALTROMAXSAMPLES  || max > ALTROMAXSAMPLES  )
78     {
79       AliWarning( Form( "Attempt to set Invalid time bin range (Min , Max) = (%d, %d), Ingored",  min, max ) ); 
80     }
81   else
82     {
83       fMinTimeIndex  = min;
84       fMaxTimeIndex  = max;
85     }
86 }
87
88
89 UShort_t 
90 AliCaloRawAnalyzer::Max(const UShort_t *data, const int length )
91 {
92   //------------
93   UShort_t tmpmax  = data[0];
94
95   for(int i=0; i < length; i++)
96     {
97       if( tmpmax  <  data[i] )
98         {
99           tmpmax = data[i];
100         }
101     }
102   return tmpmax;
103 }
104
105
106 void 
107 AliCaloRawAnalyzer::SelectSubarray( const Double_t *data, const int length, const short maxindex, int *const first,  int *const last, const int cut)
108 {
109   //Selection of subset of data from one bunch that will be used for fitting or
110   //Peak finding.  Go to the left and right of index of the maximum time bin
111   //Until the ADC value is less that fFitArrayCut, or derivative changes sign (data jump)
112   int tmpfirst =  maxindex;
113   int tmplast  =  maxindex;
114   Double_t prevFirst =  data[maxindex];
115   Double_t prevLast  =  data[maxindex];  
116   bool firstJump = false;
117   bool lastJump = false;
118
119   while( (tmpfirst >= 0) && (data[tmpfirst] >= cut ) && (!firstJump) ) 
120     {
121       // jump check:
122       if (tmpfirst != maxindex) { // neighbor to maxindex can share peak with maxindex
123         if ( data[tmpfirst] >= prevFirst) {
124           firstJump = true;
125         }
126       }
127       prevFirst = data[tmpfirst];
128       tmpfirst -- ;
129     }
130   
131   while( (tmplast < length) && (data[tmplast] >=  cut ) && (!lastJump) ) 
132     {
133       // jump check:
134       if (tmplast != maxindex) { // neighbor to maxindex can share peak with maxindex
135         if ( data[tmplast] >= prevLast) {
136           lastJump = true;
137         }
138       }
139       prevLast = data[tmplast];
140       tmplast ++;
141     }
142
143   // we keep one pre- or post- sample if we can (as in online)
144   // check though if we ended up on a 'jump', or out of bounds: if so, back up
145   if (firstJump || tmpfirst<0) tmpfirst ++;
146   if (lastJump || tmplast>=length) tmplast --;
147
148   *first = tmpfirst;
149   *last =  tmplast;
150   return;
151 }
152
153
154
155 Float_t 
156 AliCaloRawAnalyzer::ReverseAndSubtractPed( const AliCaloBunchInfo *bunch, const UInt_t /*altrocfg1*/,  const UInt_t /*altrocfg2*/, double *outarray ) const
157 {
158   //Time sample comes in reversed order, revers them back
159   //Subtract the baseline based on content of altrocfg1 and altrocfg2.
160   Int_t length =  bunch->GetLength(); 
161   const UShort_t *sig = bunch->GetData();
162
163   double ped = EvaluatePedestal( sig , length);
164
165   for( int i=0; i < length; i++ )
166     {
167       outarray[i] = sig[length -i -1] - ped;
168     }
169
170   return ped;
171 }
172
173
174
175 Float_t 
176 AliCaloRawAnalyzer::EvaluatePedestal(const UShort_t * const data, const int /*length*/ ) const
177 {
178   //  double ped = 0;
179   double tmp = 0;
180
181   if( fIsZerosupressed == false ) 
182     {
183       for(int i=0; i < fNsamplePed; i++ )
184         {
185           tmp += data[i];
186         }
187    }
188
189   return  tmp/fNsamplePed;
190  }
191
192
193 short  
194 AliCaloRawAnalyzer::Max( const AliCaloBunchInfo *const bunch , int *const maxindex )
195 {
196   //comment
197   short tmpmax = -1;
198   int tmpindex = -1;
199   const UShort_t *sig = bunch->GetData();
200
201   for(int i=0; i < bunch->GetLength(); i++ )
202     {
203       if( sig[i] > tmpmax )
204         {
205           tmpmax   =  sig[i];
206           tmpindex =  i; 
207         }
208     }
209   
210   if(maxindex != 0 )
211     {
212       //   *maxindex =  bunch->GetLength() -1 - tmpindex + bunch->GetStartBin(); 
213        *maxindex =  bunch->GetLength() -1 - tmpindex + bunch->GetStartBin(); 
214     }
215   
216   return  tmpmax;
217 }
218
219
220 bool  
221 AliCaloRawAnalyzer::CheckBunchEdgesForMax( const AliCaloBunchInfo *const bunch ) const
222 {
223   // a bunch is considered invalid if the maximum is in the first or last time-bin
224   short tmpmax = -1;
225   int tmpindex = -1;
226   const UShort_t *sig = bunch->GetData();
227
228   for(int i=0; i < bunch->GetLength(); i++ )
229     {
230       if( sig[i] > tmpmax )
231         {
232           tmpmax   =  sig[i];
233           tmpindex =  i; 
234         }
235     }
236   
237   bool bunchOK = true;
238   if (tmpindex == 0 || tmpindex == (bunch->GetLength()-1) )
239     {
240       bunchOK = false;
241     }
242   
243   return  bunchOK;
244 }
245
246
247 int 
248 AliCaloRawAnalyzer::SelectBunch( const vector<AliCaloBunchInfo> &bunchvector,short *const maxampbin, short *const maxamplitude )
249 {
250   //We select the bunch with the highest amplitude unless any time constraints is set
251   short max = -1;
252   short bunchindex = -1;
253   short maxall = -1;
254   int indx = -1;
255
256   for(unsigned int i=0; i < bunchvector.size(); i++ )
257     {
258       max = Max(  &bunchvector.at(i), &indx ); // CRAP PTH, bug fix, trouble if more than one bunches  
259       if( IsInTimeRange( indx, fMaxTimeIndex, fMinTimeIndex) )
260         {
261           if( max > maxall )
262             {
263               maxall = max;
264               bunchindex = i;
265               *maxampbin     = indx;
266               *maxamplitude  = max;
267             }
268         }
269     }
270  
271   if (bunchindex >= 0) {
272     bool bunchOK = CheckBunchEdgesForMax( &bunchvector.at(bunchindex) );
273     if (! bunchOK) { 
274       bunchindex = -1; 
275     }
276   }
277
278   return  bunchindex;
279 }
280
281
282 bool 
283 AliCaloRawAnalyzer::IsInTimeRange( const int maxindex, const int maxtindx, const int mintindx )
284 {
285   // Ckeck if the index of the max ADC vaue is consistent with trigger.
286   if( ( mintindx  < 0 && maxtindx   < 0) ||maxtindx  < 0 )
287     {
288       return true; 
289     }
290   return ( maxindex < maxtindx ) && ( maxindex > mintindx  ) ? true : false;
291 }
292
293
294 void 
295 AliCaloRawAnalyzer::PrintBunches( const vector<AliCaloBunchInfo> &bvctr ) const
296 {
297   //comment
298   cout << __FILE__ << __LINE__<< "*************** Printing Bunches *******************" << endl;
299   cout << __FILE__ << __LINE__<< "*** There are " << bvctr.size() << ", bunches" << endl;
300
301   for(unsigned int i=0; i < bvctr.size() ; i++ )
302     {
303       PrintBunch( bvctr.at(i) );
304       cout << " bunch = "  <<  i  << endl;
305     }
306   cout << __FILE__ << __LINE__<< "*************** Done ... *******************" << endl;
307 }
308
309
310 void 
311 AliCaloRawAnalyzer::PrintBunch( const AliCaloBunchInfo &bunch ) const
312 {
313   //comment
314   cout << __FILE__ << ":" << __LINE__ << endl;
315   cout << __FILE__ << __LINE__   << ", startimebin = " << bunch.GetStartBin() << ", length =" <<  bunch.GetLength() << endl;
316   const UShort_t *sig =  bunch.GetData();  
317   
318   for ( Int_t j = 0; j <  bunch.GetLength();  j++) 
319     {
320       printf("%d\t", sig[j] );
321     }
322   cout << endl; 
323 }
324
325
326 Double_t
327 AliCaloRawAnalyzer::CalculateChi2(const Double_t amp, const Double_t time,
328                                   const Int_t first, const Int_t last,
329                                   const Double_t adcErr, 
330                                   const Double_t tau)
331 {
332   //   Input:
333   //   amp   - max amplitude;
334   //   time    - time of max amplitude; 
335   //   first, last - sample array indices to be used
336   //   adcErr   - nominal error of amplitude measurement (one value for all channels)
337   //           if adcErr<0 that mean adcErr=1.
338   //   tau   - filter time response (in timebin units)
339   // Output:
340   //   chi2 - chi2
341
342   if (first == last || first<0 ) { // signal consists of single sample, chi2 estimate (0) not too well defined.. 
343     // or, first is negative, the indices are not valid
344     return Ret::kDummy;
345   }
346
347   int nsamples =  last - first + 1;
348   // printf(" AliCaloRawAnalyzer::CalculateChi2 : first %i last %i : nsamples %i : amp %3.2f time %3.2f \n", first, last, nsamples, amp, time); 
349
350   Int_t x = 0;
351   Double_t chi2 = 0;
352   Double_t dy = 0.0, xx = 0.0, f=0.0;
353
354   for (Int_t i=0; i<nsamples; i++) {
355     x     = first + i; // timebin
356     xx    = (x - time + tau) / tau; // help variable
357     f     = 0;
358     if (xx > 0) {
359       f = amp * xx*xx * TMath::Exp(2 * (1 - xx )) ;
360     }
361     dy    = fReversed[x] - f; 
362     chi2 += dy*dy;
363     // printf(" AliCaloRawAnalyzer::CalculateChi2 : %i : y %f -> f %f : dy %f \n", i, fReversed[first+i], f, dy); 
364   }
365
366   if (adcErr>0.0) { // weight chi2
367     chi2 /= (adcErr*adcErr);
368   }
369   return chi2;
370 }
371
372
373 void
374 AliCaloRawAnalyzer::CalculateMeanAndRMS(const Int_t first, const Int_t last,
375                                         Double_t & mean, Double_t & rms)
376 {
377   //   Input:
378   //   first, last - sample array indices to be used
379   // Output:
380   //   mean and RMS of samples 
381   //
382   // To possibly be used to differentiate good signals from bad before fitting
383   // 
384   mean = Ret::kDummy;
385   rms =  Ret::kDummy;
386
387   if (first == last || first<0 ) { // signal consists of single sample, chi2 estimate (0) not too well defined.. 
388     // or, first is negative, the indices are not valid
389     return;
390   }
391
392   int nsamples =  last - first + 1;
393   //  printf(" AliCaloRawAnalyzer::CalculateMeanAndRMS : first %i last %i : nsamples %i \n", first, last, nsamples); 
394
395   int x = 0;
396   Double_t sampleSum = 0; // sum of samples
397   Double_t squaredSampleSum = 0; // sum of samples squared
398
399   for (Int_t i=0; i<nsamples; i++) {
400     x = first + i;
401     sampleSum += fReversed[x];
402     squaredSampleSum += (fReversed[x] * fReversed[x]);
403   }
404
405   mean = sampleSum / nsamples;   
406   Double_t squaredMean = squaredSampleSum / nsamples;    
407   // The variance (rms squared) is equal to the mean of the squares minus the square of the mean..       
408   rms = sqrt(squaredMean - mean*mean);
409
410   return;
411 }
412
413
414
415 // AliCaloFitResults
416 // AliCaloRawAnalyzer::Evaluate( const vector<AliCaloBunchInfo>  &/*bunchvector*/, const UInt_t /*altrocfg1*/,  const UInt_t /*altrocfg2*/)
417 // { // method to do the selection of what should possibly be fitted
418 //   // not implemented for base class
419 //   cout << __FILE__ << ":" << __LINE__ << " " << endl;
420   
421 //   return AliCaloFitResults( 0, 0 );
422 // }
423
424
425
426 int
427 AliCaloRawAnalyzer::PreFitEvaluateSamples( const vector<AliCaloBunchInfo>  &bunchvector, const UInt_t altrocfg1,  
428                                            const UInt_t altrocfg2, Int_t & index, Float_t & maxf, short & maxamp, 
429                                            short & maxrev, Float_t & ped, int & first, int & last,const int acut )
430 { // method to do the selection of what should possibly be fitted
431   int nsamples  = 0;
432   short maxampindex = 0;
433   index = SelectBunch( bunchvector,  &maxampindex,  &maxamp ); // select the bunch with the highest amplitude unless any time constraints is set
434
435   
436   if( index >= 0 && maxamp >= acut ) // something valid was found, and non-zero amplitude
437     {
438       // use more convenient numbering and possibly subtract pedestal
439       ped  = ReverseAndSubtractPed( &(bunchvector.at(index)),  altrocfg1, altrocfg2, fReversed  );
440       maxf = TMath::MaxElement( bunchvector.at(index).GetLength(),  fReversed );
441       
442       if ( maxf >= acut  ) // possibly significant signal
443         {
444           // select array around max to possibly be used in fit
445           maxrev = maxampindex - bunchvector.at(index).GetStartBin(); 
446           SelectSubarray( fReversed,  bunchvector.at(index).GetLength(),  maxrev, &first, &last, acut );
447
448           // sanity check: maximum should not be in first or last bin
449           // if we should do a fit
450           if (first!=maxrev && last!=maxrev) 
451             {
452               // calculate how many samples we have 
453               nsamples =  last - first + 1;       
454             }
455         }
456     }
457
458   return nsamples;
459 }
460