Removing obsolete histograms. Fix Coverity reports
[u/mrichter/AliRoot.git] / EMCAL / AliCaloRawAnalyzerPeakFinder.cxx
CommitLineData
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>
92d9f317 39
d655d7dd 40using namespace std;
41
bab48e74 42
e37e3c84 43ClassImp( AliCaloRawAnalyzerPeakFinder )
44
bab48e74 45
46AliCaloRawAnalyzerPeakFinder::AliCaloRawAnalyzerPeakFinder() :AliCaloRawAnalyzer("Peak-Finder", "PF"),
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
60void
61AliCaloRawAnalyzerPeakFinder::InitOCDB(bool alien) const
62{
63 // Setting the default OCDB pathe depending on wether we work locally or on the GRID.
b2eed176 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 }
bab48e74 69}
d655d7dd 70
d655d7dd 71
bab48e74 72void
73AliCaloRawAnalyzerPeakFinder::ResetVectors()
74{
fc7cd737 75 //As name implies
168c7b3c 76 for(int i=0; i < PF::MAXSTART; i++)
d655d7dd 77 {
168c7b3c 78 for(int j=0; j < PF::SAMPLERANGE; j++ )
d655d7dd 79 {
d655d7dd 80 for(int k=0; k < 100; k++ )
81 {
82 fPFAmpVectors[i][j][k] = 0;
83 fPFTofVectors[i][j][k] = 0;
48a2e3eb 84 fPFAmpVectorsCoarse[i][j][k] = 0;
85 fPFTofVectorsCoarse[i][j][k] = 0;
d655d7dd 86 }
87 }
88 }
d655d7dd 89}
90
91
57839add 92AliCaloRawAnalyzerPeakFinder::~AliCaloRawAnalyzerPeakFinder()
d655d7dd 93{
94 //comment
d655d7dd 95}
96
97
48a2e3eb 98Double_t
99AliCaloRawAnalyzerPeakFinder::ScanCoarse(const Double_t *const array, const int length ) const
100{
bab48e74 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
48a2e3eb 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
57839add 119AliCaloFitResults
120AliCaloRawAnalyzerPeakFinder::Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1, const UInt_t altrocfg2 )
d655d7dd 121{
1f847d9f 122 // Evaluation of amplitude and TOF
168c7b3c 123 if( fIsInitialized == false )
124 {
125 cout << __FILE__ << ":" << __LINE__ << "ERROR, peakfinder vectors not loaded" << endl;
126 return AliCaloFitResults(kInvalid, kInvalid);
127 }
128
d655d7dd 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
bab48e74 135 fAmp = 0;
d655d7dd 136 int index = SelectBunch( bunchvector, &maxampindex, &maxamp );
168c7b3c 137
d655d7dd 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 );
f57baa2d 142 short timebinOffset = maxampindex - (bunchvector.at( index ).GetLength()-1);
168c7b3c 143
2cd0ffda 144 if( maxf < fAmpCut || ( maxamp - ped) > fOverflowCut ) // (maxamp - ped) > fOverflowCut = Close to saturation (use low gain then)
d655d7dd 145 {
168c7b3c 146 return AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset);
2cd0ffda 147 }
148 else if ( maxf >= fAmpCut )
149 {
150 int first = 0;
151 int last = 0;
152 short maxrev = maxampindex - bunchvector.at(index).GetStartBin();
92d9f317 153 SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxrev, &first, &last, fFitArrayCut);
d655d7dd 154 int nsamples = last - first;
2cd0ffda 155
156 if( ( nsamples ) >= fNsampleCut ) // no if statement needed really; keep for readability
d655d7dd 157 {
158 int startbin = bunchvector.at(index).GetStartBin();
48a2e3eb 159 int n = last - first;
d655d7dd 160 int pfindex = n - fNsampleCut;
168c7b3c 161 pfindex = pfindex > PF::SAMPLERANGE ? PF::SAMPLERANGE : pfindex;
8ffb0474 162 int dt = maxampindex - startbin -2;
8ffb0474 163 int tmpindex = 0;
48a2e3eb 164 Float_t tmptof = ScanCoarse( &fReversed[dt] , n );
165
166 if( tmptof < -1 )
167 {
168 tmpindex = 0;
169 }
170 else
171 if( tmptof > -1 && tmptof < 100 )
172 {
bab48e74 173 tmpindex = 1;
48a2e3eb 174 }
175 else
176 {
177 tmpindex = 2;
178 }
bab48e74 179
8ffb0474 180 double tof = 0;
168c7b3c 181 for(int k=0; k < PF::SAMPLERANGE; k++ )
8ffb0474 182 {
183 tof += fPFTofVectors[0][pfindex][k]*fReversed[ dt +k + tmpindex -1 ];
184 }
168c7b3c 185 for( int i=0; i < PF::SAMPLERANGE; i++ )
bab48e74 186 {
187 {
f67ce0df 188
bab48e74 189 fAmp += fPFAmpVectors[0][pfindex][i]*fReversed[ dt +i +tmpindex -1 ];
190 }
191 }
f67ce0df 192
bab48e74 193 if( TMath::Abs( (maxf - fAmp )/maxf ) > 0.1 )
48a2e3eb 194 {
bab48e74 195 fAmp = maxf;
48a2e3eb 196 }
48a2e3eb 197
92d9f317 198 tof = timebinOffset - 0.01*tof/fAmp - fL1Phase/TIMEBINWITH; // clock
2cd0ffda 199 Float_t chi2 = CalculateChi2(fAmp, tof-timebinOffset+maxrev, first, last);
200 Int_t ndf = last - first - 1; // nsamples - 2
168c7b3c 201 return AliCaloFitResults( maxamp, ped , Ret::kFitPar, fAmp, tof,
2cd0ffda 202 timebinOffset, chi2, ndf,
168c7b3c 203 Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
d655d7dd 204 }
205 else
206 {
2cd0ffda 207 Float_t chi2 = CalculateChi2(maxf, maxrev, first, last);
208 Int_t ndf = last - first - 1; // nsamples - 2
168c7b3c 209 return AliCaloFitResults( maxamp, ped , Ret::kCrude, maxf, timebinOffset,
210 timebinOffset, chi2, ndf, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
d655d7dd 211 }
2cd0ffda 212 } // ampcut
d655d7dd 213 }
168c7b3c 214 return AliCaloFitResults(kInvalid, kInvalid);
d655d7dd 215}
216
217
bab48e74 218void
168c7b3c 219AliCaloRawAnalyzerPeakFinder::CopyVectors( const AliCaloPeakFinderVectors *const pfv )
bab48e74 220{
221 // As name implies
bab48e74 222 if ( pfv != 0)
223 {
168c7b3c 224 for(int i = 0; i < PF::MAXSTART ; i++)
bab48e74 225 {
168c7b3c 226 for( int j=0; j < PF::SAMPLERANGE; j++)
bab48e74 227 {
0cb95545 228 pfv->GetVector( i, j, fPFAmpVectors[i][j] , fPFTofVectors[i][j],
229 fPFAmpVectorsCoarse[i][j] , fPFTofVectorsCoarse[i][j] );
168c7b3c 230
bab48e74 231 fPeakFinderVectors->SetVector( i, j, fPFAmpVectors[i][j], fPFTofVectors[i][j],
232 fPFAmpVectorsCoarse[i][j], fPFTofVectorsCoarse[i][j] );
233 }
234 }
235 }
236 else
237 {
238 AliFatal( "pfv = ZERO !!!!!!!");
bab48e74 239 }
240}
241
242
bab48e74 243void
244AliCaloRawAnalyzerPeakFinder::LoadVectorsOCDB()
245{
246 //Loading of Peak-Finder vectors from the
247 //Offline Condition Database (OCDB)
bab48e74 248 AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/PeakFinder/");
11bcc677 249
bab48e74 250 if( entry != 0 )
11bcc677 251 {
f73babef 252 //cout << __FILE__ << ":" << __LINE__ << ": Printing metadata !! " << endl;
253 //entry->PrintMetaData();
11bcc677 254 AliCaloPeakFinderVectors *pfv = (AliCaloPeakFinderVectors *)entry->GetObject();
255 if( pfv == 0 )
bab48e74 256 {
11bcc677 257 cout << __FILE__ << ":" << __LINE__ << "_ ERRROR " << endl;
258 }
11bcc677 259 CopyVectors( pfv );
1f847d9f 260
11bcc677 261 if( pfv != 0 )
262 {
263 fIsInitialized = true;
168c7b3c 264 }
11bcc677 265 }
bab48e74 266}
267
268
d655d7dd 269void
bab48e74 270AliCaloRawAnalyzerPeakFinder::LoadVectorsASCII()
d655d7dd 271{
bab48e74 272 //Read in the Peak finder vecors from ASCI files
168c7b3c 273 fIsInitialized= true;
a51e676d 274 const Int_t buffersize = 256;
168c7b3c 275 for(int i = 0; i < PF::MAXSTART ; i++)
a51e676d 276 {
277 for( int j=0; j < PF::SAMPLERANGE; j++)
d655d7dd 278 {
a51e676d 279 char filenameCoarse[buffersize];
280 char filename[buffersize];
281 int n = j+fNsampleCut;
282 double start = (double)i+0;
283
284 snprintf(filename, buffersize, "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt1.0.txt", getenv("ALICE_ROOT"), start, n);
285 snprintf(filenameCoarse, buffersize, "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt3.0.txt", getenv("ALICE_ROOT"), start, n);
286
287 FILE *fp = fopen(filename, "r");
288 FILE *fpc = fopen(filenameCoarse, "r");
289
290 if( fp == 0 )
d655d7dd 291 {
292 AliFatal( Form( "could not open file: %s", filename ) );
293 }
a51e676d 294 else if(fpc == 0)
48a2e3eb 295 {
296 AliFatal( Form( "could not open file: %s", filenameCoarse ) );
297 }
a51e676d 298 else
d655d7dd 299 {
300 for(int m = 0; m < n ; m++ )
a51e676d 301 {
2d35e3b8 302 fscanf(fp, "%lf\t", &fPFAmpVectors[i][j][m] );
a51e676d 303 fscanf(fpc, "%lf\t", &fPFAmpVectorsCoarse[i][j][m] );
304 }
48a2e3eb 305 fscanf(fp, "\n" );
306 fscanf(fpc, "\n" );
d655d7dd 307 for(int m = 0; m < n ; m++ )
a51e676d 308 {
309 fscanf(fp, "%lf\t", &fPFTofVectors[i][j][m] );
310 fscanf(fpc, "%lf\t", &fPFTofVectorsCoarse[i][j][m] );
311 }
d655d7dd 312
bab48e74 313 fPeakFinderVectors->SetVector( i, j, fPFAmpVectors[i][j], fPFTofVectors[i][j],
a51e676d 314 fPFAmpVectorsCoarse[i][j], fPFTofVectorsCoarse[i][j] );
315
d655d7dd 316 }
11bcc677 317
2d35e3b8 318 if(fp) fclose (fp );
319 if(fpc)fclose (fpc);
11bcc677 320
d655d7dd 321 }
a51e676d 322 }
d655d7dd 323}
48a2e3eb 324
325
bab48e74 326void
327AliCaloRawAnalyzerPeakFinder::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
339void
340AliCaloRawAnalyzerPeakFinder::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}