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