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