Adding Sqitch - Use - Qtot/Qmax
[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>
d655d7dd 39using namespace std;
40
bab48e74 41
e37e3c84 42ClassImp( AliCaloRawAnalyzerPeakFinder )
43
bab48e74 44
45AliCaloRawAnalyzerPeakFinder::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
60void
61AliCaloRawAnalyzerPeakFinder::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 69void
70AliCaloRawAnalyzerPeakFinder::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 89AliCaloRawAnalyzerPeakFinder::~AliCaloRawAnalyzerPeakFinder()
d655d7dd 90{
91 //comment
d655d7dd 92}
93
94
48a2e3eb 95Double_t
96AliCaloRawAnalyzerPeakFinder::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 116AliCaloFitResults
117AliCaloRawAnalyzerPeakFinder::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 219void
168c7b3c 220AliCaloRawAnalyzerPeakFinder::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 244void
245AliCaloRawAnalyzerPeakFinder::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 271void
bab48e74 272AliCaloRawAnalyzerPeakFinder::LoadVectorsASCII()
d655d7dd 273{
bab48e74 274 //Read in the Peak finder vecors from ASCI files
168c7b3c 275 fIsInitialized= true;
276
277 for(int i = 0; i < PF::MAXSTART ; i++)
d655d7dd 278 {
168c7b3c 279 for( int j=0; j < PF::SAMPLERANGE; j++)
d655d7dd 280 {
48a2e3eb 281 char filenameCoarse[256];
d655d7dd 282 char filename[256];
48a2e3eb 283 int n = j+fNsampleCut;
48a2e3eb 284 double start = (double)i+0;
bab48e74 285
48a2e3eb 286 sprintf(filename, "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt1.0.txt", getenv("ALICE_ROOT"), start, n);
287 sprintf(filenameCoarse, "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt3.0.txt", getenv("ALICE_ROOT"), start, n);
bab48e74 288
48a2e3eb 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 }
48a2e3eb 296 if(fpc == 0)
297 {
298 AliFatal( Form( "could not open file: %s", filenameCoarse ) );
299 }
d655d7dd 300 else
301 {
302 for(int m = 0; m < n ; m++ )
303 {
bab48e74 304 fscanf(fp, "%lf\t", &fPFAmpVectors[i][j][m] );
48a2e3eb 305 fscanf(fpc, "%lf\t", &fPFAmpVectorsCoarse[i][j][m] );
d655d7dd 306 }
48a2e3eb 307 fscanf(fp, "\n" );
308 fscanf(fpc, "\n" );
d655d7dd 309 for(int m = 0; m < n ; m++ )
310 {
48a2e3eb 311 fscanf(fp, "%lf\t", &fPFTofVectors[i][j][m] );
312 fscanf(fpc, "%lf\t", &fPFTofVectorsCoarse[i][j][m] );
d655d7dd 313 }
314
bab48e74 315 fPeakFinderVectors->SetVector( i, j, fPFAmpVectors[i][j], fPFTofVectors[i][j],
316 fPFAmpVectorsCoarse[i][j], fPFTofVectorsCoarse[i][j] );
317
d655d7dd 318 fclose (fp);
48a2e3eb 319 fclose (fpc);
d655d7dd 320 }
321 }
322 }
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}