1 /**************************************************************************
2 * This file is property of and copyright by *
3 * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009 *
5 * Primary Author: Per Thomas Hille <perthomas.hille@yale.edu> *
7 * Contributors are mentioned in the code where appropriate. *
8 * Please report bugs to p.t.hille@fys.uio.no *
10 * Permission to use, copy, modify and distribute this software and its *
11 * documentation strictly for non-commercial purposes is hereby granted *
12 * without fee, provided that the above copyright notice appears in all *
13 * copies and that both the copyright notice and this permission notice *
14 * appear in the supporting documentation. The authors make no claims *
15 * about the suitability of this software for any purpose. It is *
16 * provided "as is" without express or implied warranty. *
17 **************************************************************************/
19 // The Peak-Finder algorithm
20 // The amplitude is extracted as a
21 // weighted sum of the samples using the
22 // best possible weights.
23 // The wights is calculated only once and the
24 // Actual extraction of amplitude and peak position
25 // Is done with a simple vector multiplication, allowing for
26 // Extreemely fast computations.
28 #include "AliCaloRawAnalyzerPeakFinder.h"
29 #include "AliCaloBunchInfo.h"
30 #include "AliCaloFitResults.h"
38 ClassImp( AliCaloRawAnalyzerPeakFinder )
40 AliCaloRawAnalyzerPeakFinder::AliCaloRawAnalyzerPeakFinder() :AliCaloRawAnalyzer("Peak-Finder", "PF")
48 for(int i=0; i < MAXSTART; i++)
50 for(int j=0; j < SAMPLERANGE; j++ )
52 fPFAmpVectors[i][j] = new double[100];
53 fPFTofVectors[i][j] = new double[100];
54 fPFAmpVectorsCoarse[i][j] = new double[100];
55 fPFTofVectorsCoarse[i][j] = new double[100];
57 for(int k=0; k < 100; k++ )
59 fPFAmpVectors[i][j][k] = 0;
60 fPFTofVectors[i][j][k] = 0;
61 fPFAmpVectorsCoarse[i][j][k] = 0;
62 fPFTofVectorsCoarse[i][j][k] = 0;
72 AliCaloRawAnalyzerPeakFinder::~AliCaloRawAnalyzerPeakFinder()
75 for(int i=0; i < MAXSTART; i++)
77 for(int j=0; j < SAMPLERANGE; j++ )
79 delete[] fPFAmpVectors[i][j];
80 delete[] fPFTofVectors[i][j];
81 delete[] fPFAmpVectorsCoarse[i][j];
82 delete[] fPFTofVectorsCoarse[i][j];
89 AliCaloRawAnalyzerPeakFinder::ScanCoarse(const Double_t *const array, const int length ) const
94 for(int i=0; i < length; i++)
96 tmpTof += fPFTofVectorsCoarse[0][length][i]*array[i];
97 tmpAmp += fPFAmpVectorsCoarse[0][length][i]*array[i];
100 tmpTof = tmpTof / tmpAmp ;
106 AliCaloRawAnalyzerPeakFinder::Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1, const UInt_t altrocfg2 )
108 // Extracting the amplitude using the Peak-Finder algorithm
109 // The amplitude is a weighted sum of the samples using
112 short maxampindex; //index of maximum amplitude
113 short maxamp; //Maximum amplitude
122 // cout << __FILE__ << __LINE__ << "\tendbin = " << bunchvector.at(index).GetEndBin() << "\tstartbin = " << bunchvector.at(index).GetStartBin() << endl;
124 int index = SelectBunch( bunchvector, &maxampindex, &maxamp );
128 Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index)) , altrocfg1, altrocfg2, fReversed );
129 Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed );
131 if( maxf < fAmpCut || ( maxamp - ped) > 900 ) // (maxamp - ped) > 900 = Close to saturation (use low gain then)
133 // cout << __FILE__ << __LINE__ <<":, maxamp = " << maxamp << ", ped = "<< ped << ",. maxf = "<< maxf << ", maxampindex = "<< maxampindex << endl;
134 return AliCaloFitResults( maxamp, ped, -1, maxf, maxampindex, -1, -1 );
140 if ( maxf > fAmpCut )
144 SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxampindex - bunchvector.at(index).GetStartBin(), &first, &last);
145 int nsamples = last - first;
146 if( ( nsamples ) >= fNsampleCut )
148 int startbin = bunchvector.at(index).GetStartBin();
149 int n = last - first;
150 int pfindex = n - fNsampleCut;
151 pfindex = pfindex > SAMPLERANGE ? SAMPLERANGE : pfindex;
153 short timebinOffset = maxampindex - (bunchvector.at( index ).GetLength()-1);
155 int dt = maxampindex - startbin -2;
157 // cout << __FILE__ << __LINE__ <<"\t The coarse estimated t0 is " << ScanCoarse( &fReversed[dt] , n ) << endl;
160 // Float_t tmptof = ScanCoarse( &fReversed[dt] , n );
162 // cout << __FILE__ << __LINE__ << ", dt = " << dt << ",\tmaxamindex = " << maxampindex << "\tstartbin = "<< startbin << endl;
164 for( int i=0; i < SAMPLERANGE; i++ )
166 for( int j = 0; j < 3; j++ )
168 // fAmpA[j] += fPFAmpVectors[0][pfindex][i]*tmp[j];
169 fAmpA[j] += fPFAmpVectors[0][pfindex][i]*fReversed[ dt +i +j -1 ];
176 for(int k=0; k < 3; k ++)
178 // cout << __FILE__ << __LINE__ << "amp[="<< k <<"] = " << fAmpA[k] << endl;
179 if( TMath::Abs(fAmpA[k] - ( maxamp - ped) ) < diff)
181 diff = TMath::Abs(fAmpA[k] - ( maxamp - ped));
186 Float_t tmptof = ScanCoarse( &fReversed[dt] , n );
193 if( tmptof > -1 && tmptof < 100 )
203 for(int k=0; k < SAMPLERANGE; k++ )
205 tof += fPFTofVectors[0][pfindex][k]*fReversed[ dt +k + tmpindex -1 ];
208 // cout << __FILE__ << __LINE__ << "tofRaw = "<< tof / fAmpA[tmpindex] << endl;
210 // tof = tof / fAmpA[tmpindex] + (dt + startbin)*100;
212 if( TMath::Abs( (maxf - fAmpA[tmpindex])/maxf ) > 0.1 )
214 fAmpA[tmpindex] = maxf;
219 // tof = (dt + startbin + tmpindex )*100 - tof/fAmpA[tmpindex];
220 // tof = ( timebinOffset )*100 - tof/fAmpA[tmpindex]; // ns
221 tof = timebinOffset - 0.01*tof/fAmpA[tmpindex]; // clock ticks
223 // tof = tof/fAmpA[tmpindex];
226 return AliCaloFitResults( maxamp, ped , -1, fAmpA[tmpindex], tof, -2, -3 );
230 return AliCaloFitResults( maxamp, ped , -5, maxf, -6, -7, -8 );
234 // cout << __FILE__ << __LINE__ << "WARNING, returning amp = -1 " << endl;
235 return AliCaloFitResults(-1, -1);
240 AliCaloRawAnalyzerPeakFinder::LoadVectors()
242 //Read in the Peak finder vecors from file
243 for(int i = 0; i < MAXSTART ; i++)
245 for( int j=0; j < SAMPLERANGE; j++)
247 char filenameCoarse[256];
250 int n = j+fNsampleCut;
252 // double start = (double)i+0.5;
253 double start = (double)i+0;
255 sprintf(filename, "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt1.0.txt", getenv("ALICE_ROOT"), start, n);
256 sprintf(filenameCoarse, "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt3.0.txt", getenv("ALICE_ROOT"), start, n);
258 FILE *fp = fopen(filename, "r");
259 FILE *fpc = fopen(filenameCoarse, "r");
263 AliFatal( Form( "could not open file: %s", filename ) );
268 AliFatal( Form( "could not open file: %s", filenameCoarse ) );
272 for(int m = 0; m < n ; m++ )
274 cout << __FILE__ << __LINE__ << "i="<<i <<"\tj=" <<j << "\tm=" << m << endl;
276 fscanf(fp, "%lf\t", &fPFAmpVectors[i][j][m] );
277 // fPFAmpVectorsCoarse[i][j][m] = 1;
278 fscanf(fpc, "%lf\t", &fPFAmpVectorsCoarse[i][j][m] );
284 for(int m = 0; m < n ; m++ )
286 // fPFTofVectors[i][j][m] = 1;
288 fscanf(fp, "%lf\t", &fPFTofVectors[i][j][m] );
289 fscanf(fpc, "%lf\t", &fPFTofVectorsCoarse[i][j][m] );
290 // fPFTofVectorsCoarse[i][j][m] = 1;
305 AliCaloRawAnalyzerPeakFinder::PolTof( const double rectof ) const
308 static Double_t p0 = -55.69;
309 static Double_t p1 = 3.178;
310 static Double_t p2 = -0.05587;
311 static Double_t p3 = 0.0003185;
312 static Double_t p4 = -7.91E-7;
313 static Double_t p5 = 7.576E-10;