add protection against truncated events + coverity - Rachid
[u/mrichter/AliRoot.git] / EMCAL / AliCaloRawAnalyzerKStandard.cxx
CommitLineData
92d9f317 1/**************************************************************************
2 * This file is property of and copyright by *
3 * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009 *
4 * *
5 * Primary Author: Per Thomas Hille <p.t.hille@fys.uio.no> *
6 * *
7 * Contributors are mentioned in the code where appropriate. *
8 * Please report bugs to p.t.hille@fys.uio.no *
9 * *
10 * Permission to use, copy, modify and distribute this software and its *
11 * documentation strictly for non-commercial purposes is hereby granted *
12 * without fee, provided that the above copyright notice appears in all *
13 * copies and that both the copyright notice and this permission notice *
14 * appear in the supporting documentation. The authors make no claims *
15 * about the suitability of this software for any purpose. It is *
16 * provided "as is" without express or implied warranty. *
17 **************************************************************************/
18
12543482 19//
92d9f317 20// Extraction of amplitude and peak position
21// FRom CALO raw data using
22// least square fit for the
23// Moment assuming identical and
24// independent errors (equivalent with chi square)
25//
26
27#include "AliCaloRawAnalyzerKStandard.h"
28#include "AliCaloBunchInfo.h"
29#include "AliCaloFitResults.h"
30#include "AliLog.h"
31#include "TMath.h"
32#include <stdexcept>
33#include <iostream>
34#include "TF1.h"
35#include "TGraph.h"
36#include "TRandom.h"
37
1f847d9f 38#include "AliEMCALRawResponse.h"
39
40
92d9f317 41using namespace std;
42
92d9f317 43ClassImp( AliCaloRawAnalyzerKStandard )
44
45
396baaf6 46AliCaloRawAnalyzerKStandard::AliCaloRawAnalyzerKStandard() : AliCaloRawAnalyzerFitter("Chi Square ( kStandard )", "KStandard")
92d9f317 47{
48
49 fAlgo = Algo::kStandard;
92d9f317 50}
51
52
53AliCaloRawAnalyzerKStandard::~AliCaloRawAnalyzerKStandard()
54{
d44019b7 55 // delete fTf1;
92d9f317 56}
57
58
92d9f317 59AliCaloFitResults
60AliCaloRawAnalyzerKStandard::Evaluate( const vector<AliCaloBunchInfo> &bunchlist, const UInt_t altrocfg1, const UInt_t altrocfg2 )
61{
4074cc41 62 //Evaluation Amplitude and TOF
92d9f317 63 Float_t pedEstimate = 0;
64 short maxADC = 0;
65 Int_t first = 0;
66 Int_t last = 0;
67 Int_t bunchIndex = 0;
68 Float_t ampEstimate = 0;
69 short timeEstimate = 0;
70 Float_t time = 0;
71 Float_t amp=0;
72 Float_t chi2 = 0;
73 Int_t ndf = 0;
74 Bool_t fitDone = kFALSE;
92d9f317 75 int nsamples = PreFitEvaluateSamples( bunchlist, altrocfg1, altrocfg2, bunchIndex, ampEstimate,
12543482 76 maxADC, timeEstimate, pedEstimate, first, last, (int)fAmpCut );
92d9f317 77
78
79 if (ampEstimate >= fAmpCut )
80 {
81 time = timeEstimate;
82 Int_t timebinOffset = bunchlist.at(bunchIndex).GetStartBin() - (bunchlist.at(bunchIndex).GetLength()-1);
83 amp = ampEstimate;
84
85 if ( nsamples > 1 && maxADC< OVERFLOWCUT )
86 {
87 FitRaw(first, last, amp, time, chi2, fitDone);
88 time += timebinOffset;
89 timeEstimate += timebinOffset;
90 ndf = nsamples - 2;
91 }
92 }
93 if ( fitDone )
94 {
95 Float_t ampAsymm = (amp - ampEstimate)/(amp + ampEstimate);
96 Float_t timeDiff = time - timeEstimate;
97
98 if ( (TMath::Abs(ampAsymm) > 0.1) || (TMath::Abs(timeDiff) > 2) )
99 {
100 amp = ampEstimate;
101 time = timeEstimate;
102 fitDone = kFALSE;
103 }
104 }
105 if (amp >= fAmpCut )
106 {
107 if ( ! fitDone)
108 {
109 amp += (0.5 - gRandom->Rndm());
110 }
92d9f317 111 time = time * TIMEBINWITH;
92d9f317 112 time -= fL1Phase;
113
92d9f317 114 return AliCaloFitResults( -99, -99, fAlgo , amp, time,
12543482 115 (int)time, chi2, ndf, Ret::kDummy );
d44019b7 116 }
92d9f317 117 return AliCaloFitResults( Ret::kInvalid, Ret::kInvalid );
118}
119
92d9f317 120
92d9f317 121void
122 AliCaloRawAnalyzerKStandard::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Float_t & chi2, Bool_t & fitDone) const
d44019b7 123{
124 // Fits the raw signal time distribution
92d9f317 125 int nsamples = lastTimeBin - firstTimeBin + 1;
126 fitDone = kFALSE;
d44019b7 127 if (nsamples < 3) { return; }
92d9f317 128
129 TGraph *gSig = new TGraph( nsamples);
130
131 for (int i=0; i<nsamples; i++)
132 {
133 Int_t timebin = firstTimeBin + i;
134 gSig->SetPoint(i, timebin, GetReversed(timebin));
135 }
1f847d9f 136
1f847d9f 137 TF1 * signalF = new TF1("signal", AliEMCALRawResponse::RawResponseFunction, 0, TIMEBINS , 5);
138
92d9f317 139 signalF->SetParameters(10.,5., TAU ,ORDER,0.); //set all defaults once, just to be safe
140 signalF->SetParNames("amp","t0","tau","N","ped");
d44019b7 141 signalF->FixParameter(2,TAU);
142 signalF->FixParameter(3,ORDER);
143 signalF->FixParameter(4, 0);
92d9f317 144 signalF->SetParameter(1, time);
145 signalF->SetParameter(0, amp);
92d9f317 146 signalF->SetParLimits(0, 0.5*amp, 2*amp );
147 signalF->SetParLimits(1, time - 4, time + 4);
148
149 try {
150 gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points
92d9f317 151 amp = signalF->GetParameter(0);
152 time = signalF->GetParameter(1);
153 chi2 = signalF->GetChisquare();
154 fitDone = kTRUE;
155 }
640f438f 156 catch (const std::exception & e)
157 {
158 AliError( Form("TGraph Fit exception %s", e.what()) );
159 // stay with default amp and time in case of exception, i.e. no special action required
160 fitDone = kFALSE;
92d9f317 161 }
d44019b7 162
92d9f317 163 delete signalF;
92d9f317 164 delete gSig; // delete TGraph
92d9f317 165 return;
166}
167
168