add protection against truncated events + coverity - Rachid
[u/mrichter/AliRoot.git] / EMCAL / AliCaloRawAnalyzerPeakFinder.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   perthomas.hille@yale.edu                       *
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 // The Peak-Finder algorithm
21 // The amplitude is extracted  as a
22 // weighted sum of the samples using the 
23 // best possible weights.
24 // The wights is calculated only once and the
25 // Actual extraction of amplitude and peak position
26 // Is done with a simple vector multiplication, allowing for
27 // Extreemely fast computations. 
28
29 #include "AliCaloRawAnalyzerPeakFinder.h"
30 #include "AliCaloBunchInfo.h"
31 #include "AliCaloFitResults.h"
32 #include "TMath.h"
33 #include "AliLog.h"
34 #include "AliCDBEntry.h"
35 #include "AliCDBManager.h"
36 #include "TFile.h"
37 #include "AliCaloPeakFinderVectors.h"
38 #include <iostream>
39
40 using namespace std;
41
42
43 ClassImp( AliCaloRawAnalyzerPeakFinder )
44
45
46 AliCaloRawAnalyzerPeakFinder::AliCaloRawAnalyzerPeakFinder() :AliCaloRawAnalyzer("Peak-Finder", "PF"),  
47                                                               fPeakFinderVectors(0),
48                                                               fRunOnAlien(false),
49                                                               fIsInitialized(false)
50 {
51   //Comment
52   fAlgo= Algo::kPeakFinder;
53   InitOCDB(fRunOnAlien);
54   fPeakFinderVectors = new AliCaloPeakFinderVectors() ;
55   ResetVectors();
56   LoadVectorsOCDB();
57 }
58
59
60 void 
61 AliCaloRawAnalyzerPeakFinder::InitOCDB(bool alien) const
62 {
63   // Setting the default OCDB pathe depending on wether we work locally or on the GRID.
64   if( !AliCDBManager::Instance()->IsDefaultStorageSet ())
65     {
66       AliCDBManager::Instance()->SetDefaultStorage(  alien == true ? "alien://$ALICE_ROOT/OCDB" : "local://$ALICE_ROOT/OCDB" );
67       AliCDBManager::Instance()->SetRun(100);
68     }
69 }
70
71
72 void  
73 AliCaloRawAnalyzerPeakFinder::ResetVectors()
74 {
75   //As name implies
76   for(int i=0; i < PF::MAXSTART; i++)
77     {
78       for(int j=0; j < PF::SAMPLERANGE; j++ )
79         {
80           for(int k=0; k < 100; k++ )
81             {
82               fPFAmpVectors[i][j][k] = 0; 
83               fPFTofVectors[i][j][k] = 0;
84               fPFAmpVectorsCoarse[i][j][k] = 0;
85               fPFTofVectorsCoarse[i][j][k] = 0; 
86             }
87         }
88     }
89 }
90
91
92 AliCaloRawAnalyzerPeakFinder::~AliCaloRawAnalyzerPeakFinder()
93 {
94   //comment
95 }
96
97
98 Double_t  
99 AliCaloRawAnalyzerPeakFinder::ScanCoarse(const Double_t *const array, const int length ) const
100 {
101   // Fisrt (coarce) estimate of Amplitude using the Peak-Finder.
102   // The output of the first iteration is sued to select vectors 
103   // for the second iteration.
104
105   Double_t tmpTof = 0;
106   Double_t tmpAmp= 0;
107
108   for(int i=0; i < length; i++)
109     {
110       tmpTof += fPFTofVectorsCoarse[0][length][i]*array[i]; 
111       tmpAmp += fPFAmpVectorsCoarse[0][length][i]*array[i]; 
112     }
113   
114   tmpTof = tmpTof / tmpAmp ;
115   return tmpTof;
116 }
117
118
119 AliCaloFitResults 
120 AliCaloRawAnalyzerPeakFinder::Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1,  const UInt_t altrocfg2 )
121 {
122   // Evaluation of amplitude and TOF
123   if( fIsInitialized == false )
124     {
125       cout << __FILE__ << ":" << __LINE__ << "ERROR, peakfinder vectors not loaded" << endl;
126       return  AliCaloFitResults(kInvalid, kInvalid);
127     }
128
129   // Extracting the amplitude using the Peak-Finder algorithm
130   // The amplitude is a weighted sum of the samples using 
131   // optimum weights.
132
133   short maxampindex; //index of maximum amplitude
134   short maxamp; //Maximum amplitude
135   fAmp = 0;
136   int index = SelectBunch( bunchvector,  &maxampindex,  &maxamp );
137   
138   if( index >= 0)
139     {
140       Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index))  ,  altrocfg1, altrocfg2, fReversed  );
141       Float_t maxf = TMath::MaxElement(   bunchvector.at(index).GetLength(),  fReversed );
142       short timebinOffset = maxampindex - (bunchvector.at( index ).GetLength()-1); 
143  
144       if(  maxf < fAmpCut  ||  ( maxamp - ped) > fOverflowCut  ) // (maxamp - ped) > fOverflowCut = Close to saturation (use low gain then)
145         {
146           return  AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset);
147         }            
148       else if ( maxf >= fAmpCut )
149         {
150           int first = 0;
151           int last = 0;
152           short maxrev = maxampindex  -  bunchvector.at(index).GetStartBin();     
153           SelectSubarray( fReversed,  bunchvector.at(index).GetLength(), maxrev, &first, &last, fFitArrayCut);
154           int nsamples =  last - first;
155
156           if( ( nsamples  )  >= fNsampleCut ) // no if statement needed really; keep for readability
157             {
158               int startbin = bunchvector.at(index).GetStartBin();  
159               int n = last - first;  
160               int pfindex = n - fNsampleCut; 
161               pfindex = pfindex > PF::SAMPLERANGE ? PF::SAMPLERANGE : pfindex;
162
163               int dt =  maxampindex - startbin -2; 
164               int tmpindex = 0;
165
166
167               Float_t tmptof = ScanCoarse( &fReversed[dt] , n );
168               
169               if( tmptof < -1 )
170                 {
171                   tmpindex = 0;
172                 }
173               else
174                 if( tmptof  > -1 && tmptof < 100 )
175                   {
176                     tmpindex = 1;
177                   }
178                 else
179                   {
180                     tmpindex = 2;
181                   }
182
183               double tof = 0;
184             
185               for(int k=0; k < PF::SAMPLERANGE; k++   )
186                 {
187                   tof +=  fPFTofVectors[0][pfindex][k]*fReversed[ dt  +k + tmpindex -1 ];   
188                 }
189             
190               for( int i=0; i < PF::SAMPLERANGE; i++ )
191                 {
192                   {
193                     fAmp += fPFAmpVectors[0][pfindex][i]*fReversed[ dt  +i  +tmpindex -1 ];
194                   }
195                 }
196               if( TMath::Abs(  (maxf - fAmp  )/maxf )  >   0.1 )
197                 {
198                   fAmp = maxf;
199                 }
200               
201               tof = timebinOffset - 0.01*tof/fAmp - fL1Phase/TIMEBINWITH; // clock
202               
203               // use local-array time for chi2 estimate
204               Float_t chi2 = CalculateChi2(fAmp, tof-timebinOffset+maxrev, first, last);
205               Int_t ndf = last - first - 1; // nsamples - 2
206               return AliCaloFitResults( maxamp, ped , Ret::kFitPar, fAmp, tof, 
207                                         timebinOffset, chi2, ndf,
208                                         Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );  
209             }
210           else
211             {
212               Float_t chi2 = CalculateChi2(maxf, maxrev, first, last);
213               Int_t ndf = last - first - 1; // nsamples - 2
214               return AliCaloFitResults( maxamp, ped , Ret::kCrude, maxf, timebinOffset,
215                                         timebinOffset, chi2, ndf, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) ); 
216             }
217         } // ampcut
218     }
219   return  AliCaloFitResults(kInvalid, kInvalid);
220 }
221
222
223 void   
224 AliCaloRawAnalyzerPeakFinder::CopyVectors( const AliCaloPeakFinderVectors *const pfv )
225 {
226   // As name implies
227   if ( pfv != 0)
228     {
229       for(int i = 0;  i < PF::MAXSTART ; i++)
230         {
231           for( int j=0; j < PF::SAMPLERANGE; j++)  
232             {
233               pfv->GetVector( i, j, fPFAmpVectors[i][j] ,  fPFTofVectors[i][j],    
234                               fPFAmpVectorsCoarse[i][j] , fPFTofVectorsCoarse[i][j]  ); 
235
236               fPeakFinderVectors->SetVector( i, j, fPFAmpVectors[i][j], fPFTofVectors[i][j],    
237                                              fPFAmpVectorsCoarse[i][j], fPFTofVectorsCoarse[i][j] );   
238             }
239         }
240     }
241   else
242     {
243       AliFatal( "pfv = ZERO !!!!!!!");
244     } 
245 }
246
247
248 void   
249 AliCaloRawAnalyzerPeakFinder::LoadVectorsOCDB()
250 {
251   //Loading of Peak-Finder  vectors from the 
252   //Offline Condition Database  (OCDB)
253   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/PeakFinder/");
254   
255   if( entry != 0 )
256   {
257     cout << __FILE__ << ":" << __LINE__ << ": Printing metadata !! " << endl;
258     entry->PrintMetaData();
259     AliCaloPeakFinderVectors  *pfv = (AliCaloPeakFinderVectors *)entry->GetObject(); 
260     if( pfv == 0 )
261     {
262       cout << __FILE__ << ":" << __LINE__ << "_ ERRROR " << endl;
263     }
264     CopyVectors( pfv );
265     
266     if( pfv != 0 )
267     {
268       fIsInitialized = true;
269     }
270   }
271 }
272
273
274 void 
275 AliCaloRawAnalyzerPeakFinder::LoadVectorsASCII()
276 {
277   //Read in the Peak finder vecors from ASCI files
278   fIsInitialized= true;  
279   const Int_t buffersize = 256;
280   for(int i = 0;  i < PF::MAXSTART ; i++)
281   {
282     for( int j=0; j < PF::SAMPLERANGE; j++)
283     {
284       char filenameCoarse[buffersize];
285       char filename[buffersize];
286       int n = j+fNsampleCut;
287       double start = (double)i+0;
288       
289       snprintf(filename, buffersize,       "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt1.0.txt", getenv("ALICE_ROOT"), start, n);
290       snprintf(filenameCoarse, buffersize, "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt3.0.txt", getenv("ALICE_ROOT"), start, n);
291       
292       FILE *fp  =  fopen(filename, "r");
293       FILE *fpc =  fopen(filenameCoarse, "r");
294       
295       if( fp == 0 )
296             {
297               AliFatal( Form( "could not open file: %s", filename ) );
298             }
299       else if(fpc == 0)
300             {
301               AliFatal( Form( "could not open file: %s", filenameCoarse ) );
302             }
303       else
304             {
305               for(int m = 0; m < n ; m++ )
306         {
307           fscanf(fp,  "%lf\t", &fPFAmpVectors[i][j][m] );
308           fscanf(fpc, "%lf\t", &fPFAmpVectorsCoarse[i][j][m] );
309         }
310               fscanf(fp,   "\n" );
311               fscanf(fpc,  "\n" );
312               for(int m = 0; m < n ; m++ )
313         {
314           fscanf(fp, "%lf\t",   &fPFTofVectors[i][j][m]  );
315           fscanf(fpc, "%lf\t",  &fPFTofVectorsCoarse[i][j][m]  );  
316         }
317               
318               fPeakFinderVectors->SetVector( i, j, fPFAmpVectors[i][j], fPFTofVectors[i][j],    
319                                       fPFAmpVectorsCoarse[i][j], fPFTofVectorsCoarse[i][j] );   
320         
321             }
322       
323       if(fp) fclose (fp );
324       if(fpc)fclose (fpc);
325       
326     }
327   }
328 }
329
330
331 void   
332 AliCaloRawAnalyzerPeakFinder::WriteRootFile() const
333 {
334   // Utility function to write Peak-Finder vectors to an root file
335   // The output is used to create an OCDB entry.
336   fPeakFinderVectors->PrintVectors();
337   TFile *f = new TFile("peakfindervectors2.root",  "recreate" );
338   fPeakFinderVectors->Write();
339   f->Close();
340   delete f;
341 }
342
343
344 void 
345 AliCaloRawAnalyzerPeakFinder::PrintVectors()
346 {
347   for(int i=0; i < 20; i++)
348     {
349       for( int j = 0; j < PF::MAXSTART; j ++ )
350         {
351           for( int k=0; k < PF::SAMPLERANGE; k++ )
352             {
353               cout << fPFAmpVectors[j][k][i] << "\t" ;
354             }
355         }
356       cout << endl;
357     }
358   cout << __FILE__ << ":" << __LINE__ << ":.... DONE !!" << endl;
359 }