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