correct previous non intended commit, but leaving some new things
[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
d468be51 43
e37e3c84 44ClassImp( AliCaloRawAnalyzerPeakFinder )
45
bab48e74 46
47AliCaloRawAnalyzerPeakFinder::AliCaloRawAnalyzerPeakFinder() :AliCaloRawAnalyzer("Peak-Finder", "PF"),
fc7cd737 48 fPeakFinderVectors(0),
168c7b3c 49 fRunOnAlien(false),
50 fIsInitialized(false)
d655d7dd 51{
fc7cd737 52 //Comment
168c7b3c 53 fAlgo= Algo::kPeakFinder;
bab48e74 54 fPeakFinderVectors = new AliCaloPeakFinderVectors() ;
55 ResetVectors();
56 LoadVectorsOCDB();
bab48e74 57}
58
bab48e74 59void
60AliCaloRawAnalyzerPeakFinder::ResetVectors()
61{
fc7cd737 62 //As name implies
168c7b3c 63 for(int i=0; i < PF::MAXSTART; i++)
d655d7dd 64 {
168c7b3c 65 for(int j=0; j < PF::SAMPLERANGE; j++ )
d655d7dd 66 {
d655d7dd 67 for(int k=0; k < 100; k++ )
68 {
69 fPFAmpVectors[i][j][k] = 0;
70 fPFTofVectors[i][j][k] = 0;
48a2e3eb 71 fPFAmpVectorsCoarse[i][j][k] = 0;
72 fPFTofVectorsCoarse[i][j][k] = 0;
d655d7dd 73 }
74 }
75 }
d655d7dd 76}
77
78
57839add 79AliCaloRawAnalyzerPeakFinder::~AliCaloRawAnalyzerPeakFinder()
d655d7dd 80{
81 //comment
d655d7dd 82}
83
84
48a2e3eb 85Double_t
86AliCaloRawAnalyzerPeakFinder::ScanCoarse(const Double_t *const array, const int length ) const
87{
bab48e74 88 // Fisrt (coarce) estimate of Amplitude using the Peak-Finder.
89 // The output of the first iteration is sued to select vectors
90 // for the second iteration.
91
48a2e3eb 92 Double_t tmpTof = 0;
93 Double_t tmpAmp= 0;
94
95 for(int i=0; i < length; i++)
96 {
97 tmpTof += fPFTofVectorsCoarse[0][length][i]*array[i];
98 tmpAmp += fPFAmpVectorsCoarse[0][length][i]*array[i];
99 }
100
101 tmpTof = tmpTof / tmpAmp ;
102 return tmpTof;
103}
104
105
57839add 106AliCaloFitResults
107AliCaloRawAnalyzerPeakFinder::Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1, const UInt_t altrocfg2 )
d655d7dd 108{
1f847d9f 109 // Evaluation of amplitude and TOF
168c7b3c 110 if( fIsInitialized == false )
111 {
112 cout << __FILE__ << ":" << __LINE__ << "ERROR, peakfinder vectors not loaded" << endl;
113 return AliCaloFitResults(kInvalid, kInvalid);
114 }
115
d655d7dd 116 // Extracting the amplitude using the Peak-Finder algorithm
117 // The amplitude is a weighted sum of the samples using
118 // optimum weights.
119
120 short maxampindex; //index of maximum amplitude
121 short maxamp; //Maximum amplitude
bab48e74 122 fAmp = 0;
d655d7dd 123 int index = SelectBunch( bunchvector, &maxampindex, &maxamp );
168c7b3c 124
d655d7dd 125 if( index >= 0)
126 {
127 Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index)) , altrocfg1, altrocfg2, fReversed );
128 Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed );
f57baa2d 129 short timebinOffset = maxampindex - (bunchvector.at( index ).GetLength()-1);
168c7b3c 130
2cd0ffda 131 if( maxf < fAmpCut || ( maxamp - ped) > fOverflowCut ) // (maxamp - ped) > fOverflowCut = Close to saturation (use low gain then)
d655d7dd 132 {
168c7b3c 133 return AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset);
2cd0ffda 134 }
135 else if ( maxf >= fAmpCut )
136 {
137 int first = 0;
138 int last = 0;
139 short maxrev = maxampindex - bunchvector.at(index).GetStartBin();
92d9f317 140 SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxrev, &first, &last, fFitArrayCut);
d655d7dd 141 int nsamples = last - first;
2cd0ffda 142
143 if( ( nsamples ) >= fNsampleCut ) // no if statement needed really; keep for readability
d655d7dd 144 {
145 int startbin = bunchvector.at(index).GetStartBin();
48a2e3eb 146 int n = last - first;
d655d7dd 147 int pfindex = n - fNsampleCut;
168c7b3c 148 pfindex = pfindex > PF::SAMPLERANGE ? PF::SAMPLERANGE : pfindex;
8ffb0474 149 int dt = maxampindex - startbin -2;
8ffb0474 150 int tmpindex = 0;
48a2e3eb 151 Float_t tmptof = ScanCoarse( &fReversed[dt] , n );
152
153 if( tmptof < -1 )
154 {
155 tmpindex = 0;
156 }
157 else
158 if( tmptof > -1 && tmptof < 100 )
159 {
bab48e74 160 tmpindex = 1;
48a2e3eb 161 }
162 else
163 {
164 tmpindex = 2;
165 }
bab48e74 166
8ffb0474 167 double tof = 0;
168c7b3c 168 for(int k=0; k < PF::SAMPLERANGE; k++ )
8ffb0474 169 {
170 tof += fPFTofVectors[0][pfindex][k]*fReversed[ dt +k + tmpindex -1 ];
171 }
168c7b3c 172 for( int i=0; i < PF::SAMPLERANGE; i++ )
bab48e74 173 {
174 {
f67ce0df 175
bab48e74 176 fAmp += fPFAmpVectors[0][pfindex][i]*fReversed[ dt +i +tmpindex -1 ];
177 }
178 }
f67ce0df 179
bab48e74 180 if( TMath::Abs( (maxf - fAmp )/maxf ) > 0.1 )
48a2e3eb 181 {
bab48e74 182 fAmp = maxf;
48a2e3eb 183 }
48a2e3eb 184
92d9f317 185 tof = timebinOffset - 0.01*tof/fAmp - fL1Phase/TIMEBINWITH; // clock
2cd0ffda 186 Float_t chi2 = CalculateChi2(fAmp, tof-timebinOffset+maxrev, first, last);
187 Int_t ndf = last - first - 1; // nsamples - 2
168c7b3c 188 return AliCaloFitResults( maxamp, ped , Ret::kFitPar, fAmp, tof,
2cd0ffda 189 timebinOffset, chi2, ndf,
168c7b3c 190 Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
d655d7dd 191 }
192 else
193 {
2cd0ffda 194 Float_t chi2 = CalculateChi2(maxf, maxrev, first, last);
195 Int_t ndf = last - first - 1; // nsamples - 2
168c7b3c 196 return AliCaloFitResults( maxamp, ped , Ret::kCrude, maxf, timebinOffset,
197 timebinOffset, chi2, ndf, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
d655d7dd 198 }
2cd0ffda 199 } // ampcut
d655d7dd 200 }
168c7b3c 201 return AliCaloFitResults(kInvalid, kInvalid);
d655d7dd 202}
203
204
bab48e74 205void
168c7b3c 206AliCaloRawAnalyzerPeakFinder::CopyVectors( const AliCaloPeakFinderVectors *const pfv )
bab48e74 207{
208 // As name implies
bab48e74 209 if ( pfv != 0)
210 {
168c7b3c 211 for(int i = 0; i < PF::MAXSTART ; i++)
bab48e74 212 {
168c7b3c 213 for( int j=0; j < PF::SAMPLERANGE; j++)
bab48e74 214 {
0cb95545 215 pfv->GetVector( i, j, fPFAmpVectors[i][j] , fPFTofVectors[i][j],
216 fPFAmpVectorsCoarse[i][j] , fPFTofVectorsCoarse[i][j] );
168c7b3c 217
bab48e74 218 fPeakFinderVectors->SetVector( i, j, fPFAmpVectors[i][j], fPFTofVectors[i][j],
219 fPFAmpVectorsCoarse[i][j], fPFTofVectorsCoarse[i][j] );
220 }
221 }
222 }
223 else
224 {
225 AliFatal( "pfv = ZERO !!!!!!!");
bab48e74 226 }
227}
228
229
bab48e74 230void
231AliCaloRawAnalyzerPeakFinder::LoadVectorsOCDB()
232{
233 //Loading of Peak-Finder vectors from the
234 //Offline Condition Database (OCDB)
bab48e74 235 AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/PeakFinder/");
11bcc677 236
bab48e74 237 if( entry != 0 )
11bcc677 238 {
f73babef 239 //cout << __FILE__ << ":" << __LINE__ << ": Printing metadata !! " << endl;
240 //entry->PrintMetaData();
11bcc677 241 AliCaloPeakFinderVectors *pfv = (AliCaloPeakFinderVectors *)entry->GetObject();
242 if( pfv == 0 )
bab48e74 243 {
11bcc677 244 cout << __FILE__ << ":" << __LINE__ << "_ ERRROR " << endl;
245 }
11bcc677 246 CopyVectors( pfv );
1f847d9f 247
11bcc677 248 if( pfv != 0 )
249 {
250 fIsInitialized = true;
168c7b3c 251 }
11bcc677 252 }
bab48e74 253}
254
255
bab48e74 256void
257AliCaloRawAnalyzerPeakFinder::WriteRootFile() const
afae9650 258{ // Utility function to write Peak-Finder vectors to an root file
bab48e74 259 // The output is used to create an OCDB entry.
bab48e74 260 fPeakFinderVectors->PrintVectors();
bab48e74 261 TFile *f = new TFile("peakfindervectors2.root", "recreate" );
262 fPeakFinderVectors->Write();
263 f->Close();
264 delete f;
48a2e3eb 265}
168c7b3c 266
267
268void
269AliCaloRawAnalyzerPeakFinder::PrintVectors()
afae9650 270{ // Utility function to write Peak-Finder vectors
168c7b3c 271 for(int i=0; i < 20; i++)
272 {
273 for( int j = 0; j < PF::MAXSTART; j ++ )
274 {
275 for( int k=0; k < PF::SAMPLERANGE; k++ )
276 {
277 cout << fPFAmpVectors[j][k][i] << "\t" ;
278 }
279 }
280 cout << endl;
281 }
282 cout << __FILE__ << ":" << __LINE__ << ":.... DONE !!" << endl;
283}