SensorThickness was defined twice. Set inner chip thickness to 250mum to bypass bug...
[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   fPeakFinderVectors = new AliCaloPeakFinderVectors() ;
55   ResetVectors();
56   LoadVectorsOCDB();
57 }
58
59 void  
60 AliCaloRawAnalyzerPeakFinder::ResetVectors()
61 {
62   //As name implies
63   for(int i=0; i < PF::MAXSTART; i++)
64     {
65       for(int j=0; j < PF::SAMPLERANGE; j++ )
66         {
67           for(int k=0; k < 100; k++ )
68             {
69               fPFAmpVectors[i][j][k] = 0; 
70               fPFTofVectors[i][j][k] = 0;
71               fPFAmpVectorsCoarse[i][j][k] = 0;
72               fPFTofVectorsCoarse[i][j][k] = 0; 
73             }
74         }
75     }
76 }
77
78
79 AliCaloRawAnalyzerPeakFinder::~AliCaloRawAnalyzerPeakFinder()
80 {
81   //comment
82 }
83
84
85 Double_t  
86 AliCaloRawAnalyzerPeakFinder::ScanCoarse(const Double_t *const array, const int length ) const
87 {
88   // Fisrt (coarce) estimate of Amplitude using the Peak-Finder.
89   // The output of the first iteration is sued to select vectors 
90   // for the second iteration.
91
92   Double_t tmpTof = 0;
93   Double_t tmpAmp= 0;
94
95   for(int i=0; i < length; i++)
96     {
97       tmpTof += fPFTofVectorsCoarse[0][length][i]*array[i]; 
98       tmpAmp += fPFAmpVectorsCoarse[0][length][i]*array[i]; 
99     }
100   
101   tmpTof = tmpTof / tmpAmp ;
102   return tmpTof;
103 }
104
105
106 AliCaloFitResults 
107 AliCaloRawAnalyzerPeakFinder::Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1,  const UInt_t altrocfg2 )
108 {
109   // Evaluation of amplitude and TOF
110   if( fIsInitialized == false )
111     {
112       cout << __FILE__ << ":" << __LINE__ << "ERROR, peakfinder vectors not loaded" << endl;
113       return  AliCaloFitResults(kInvalid, kInvalid);
114     }
115
116   // Extracting the amplitude using the Peak-Finder algorithm
117   // The amplitude is a weighted sum of the samples using 
118   // optimum weights.
119
120   short maxampindex; //index of maximum amplitude
121   short maxamp; //Maximum amplitude
122   fAmp = 0;
123   int index = SelectBunch( bunchvector,  &maxampindex,  &maxamp );
124   
125   if( index >= 0)
126     {
127       Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index))  ,  altrocfg1, altrocfg2, fReversed  );
128       Float_t maxf = TMath::MaxElement(   bunchvector.at(index).GetLength(),  fReversed );
129       short timebinOffset = maxampindex - (bunchvector.at( index ).GetLength()-1); 
130  
131       if(  maxf < fAmpCut  ||  ( maxamp - ped) > fOverflowCut  ) // (maxamp - ped) > fOverflowCut = Close to saturation (use low gain then)
132         {
133           return  AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset);
134         }            
135       else if ( maxf >= fAmpCut )
136         {
137           int first = 0;
138           int last = 0;
139           short maxrev = maxampindex  -  bunchvector.at(index).GetStartBin();     
140           SelectSubarray( fReversed,  bunchvector.at(index).GetLength(), maxrev, &first, &last, fFitArrayCut);
141           int nsamples =  last - first;
142
143           if( ( nsamples  )  >= fNsampleCut ) // no if statement needed really; keep for readability
144             {
145               int startbin = bunchvector.at(index).GetStartBin();  
146               int n = last - first;  
147               int pfindex = n - fNsampleCut; 
148               pfindex = pfindex > PF::SAMPLERANGE ? PF::SAMPLERANGE : pfindex;
149               int dt =  maxampindex - startbin -2; 
150               int tmpindex = 0;
151               Float_t tmptof = ScanCoarse( &fReversed[dt] , n );
152               
153               if( tmptof < -1 )
154                 {
155                   tmpindex = 0;
156                 }
157               else
158                 if( tmptof  > -1 && tmptof < 100 )
159                   {
160                     tmpindex = 1;
161                   }
162                 else
163                   {
164                     tmpindex = 2;
165                   }
166
167               double tof = 0;
168               for(int k=0; k < PF::SAMPLERANGE; k++   )
169                 {
170                   tof +=  fPFTofVectors[0][pfindex][k]*fReversed[ dt  +k + tmpindex -1 ];   
171                 }
172               for( int i=0; i < PF::SAMPLERANGE; i++ )
173                 {
174                   {
175                     
176                     fAmp += fPFAmpVectors[0][pfindex][i]*fReversed[ dt  +i  +tmpindex -1 ];
177                   }
178                 }
179
180               if( TMath::Abs(  (maxf - fAmp  )/maxf )  >   0.1 )
181                 {
182                   fAmp = maxf;
183                 }
184               
185               tof = timebinOffset - 0.01*tof/fAmp - fL1Phase/TIMEBINWITH; // clock
186               Float_t chi2 = CalculateChi2(fAmp, tof-timebinOffset+maxrev, first, last);
187               Int_t ndf = last - first - 1; // nsamples - 2
188               return AliCaloFitResults( maxamp, ped , Ret::kFitPar, fAmp, tof, 
189                                         timebinOffset, chi2, ndf,
190                                         Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );  
191             }
192           else
193             {
194               Float_t chi2 = CalculateChi2(maxf, maxrev, first, last);
195               Int_t ndf = last - first - 1; // nsamples - 2
196               return AliCaloFitResults( maxamp, ped , Ret::kCrude, maxf, timebinOffset,
197                                         timebinOffset, chi2, ndf, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) ); 
198             }
199         } // ampcut
200     }
201   return  AliCaloFitResults(kInvalid, kInvalid);
202 }
203
204
205 void   
206 AliCaloRawAnalyzerPeakFinder::CopyVectors( const AliCaloPeakFinderVectors *const pfv )
207 {
208   // As name implies
209   if ( pfv != 0)
210     {
211       for(int i = 0;  i < PF::MAXSTART ; i++)
212         {
213           for( int j=0; j < PF::SAMPLERANGE; j++)  
214             {
215               pfv->GetVector( i, j, fPFAmpVectors[i][j] ,  fPFTofVectors[i][j],    
216                               fPFAmpVectorsCoarse[i][j] , fPFTofVectorsCoarse[i][j]  ); 
217
218               fPeakFinderVectors->SetVector( i, j, fPFAmpVectors[i][j], fPFTofVectors[i][j],    
219                                              fPFAmpVectorsCoarse[i][j], fPFTofVectorsCoarse[i][j] );   
220             }
221         }
222     }
223   else
224     {
225       AliFatal( "pfv = ZERO !!!!!!!");
226     } 
227 }
228
229
230 void   
231 AliCaloRawAnalyzerPeakFinder::LoadVectorsOCDB()
232 {
233   //Loading of Peak-Finder  vectors from the 
234   //Offline Condition Database  (OCDB)
235   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/PeakFinder/");
236   
237   if( entry != 0 )
238   {
239     //cout << __FILE__ << ":" << __LINE__ << ": Printing metadata !! " << endl;
240     //entry->PrintMetaData();
241     AliCaloPeakFinderVectors  *pfv = (AliCaloPeakFinderVectors *)entry->GetObject(); 
242     if( pfv == 0 )
243     {
244       cout << __FILE__ << ":" << __LINE__ << "_ ERRROR " << endl;
245     }
246     CopyVectors( pfv );
247     
248     if( pfv != 0 )
249     {
250       fIsInitialized = true;
251     }
252   }
253 }
254
255
256 void   
257 AliCaloRawAnalyzerPeakFinder::WriteRootFile() const
258 { // Utility function to write Peak-Finder vectors to an root file
259   // The output is used to create an OCDB entry.
260   fPeakFinderVectors->PrintVectors();
261   TFile *f = new TFile("peakfindervectors2.root",  "recreate" );
262   fPeakFinderVectors->Write();
263   f->Close();
264   delete f;
265 }
266
267
268 void 
269 AliCaloRawAnalyzerPeakFinder::PrintVectors()
270 { // Utility function to write Peak-Finder vectors 
271   for(int i=0; i < 20; i++)
272     {
273       for( int j = 0; j < PF::MAXSTART; j ++ )
274         {
275           for( int k=0; k < PF::SAMPLERANGE; k++ )
276             {
277               cout << fPFAmpVectors[j][k][i] << "\t" ;
278             }
279         }
280       cout << endl;
281     }
282   cout << __FILE__ << ":" << __LINE__ << ":.... DONE !!" << endl;
283 }