1 /**************************************************************************
2 * This file is property of and copyright by *
3 * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009 *
5 * Primary Author: Per Thomas Hille <p.t.hille@fys.uio.no> *
7 * Contributors are mentioned in the code where appropriate. *
8 * Please report bugs to p.t.hille@fys.uio.no *
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 **************************************************************************/
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)
27 #include "AliCaloRawAnalyzerKStandard.h"
28 #include "AliCaloBunchInfo.h"
29 #include "AliCaloFitResults.h"
40 ClassImp( AliCaloRawAnalyzerKStandard )
43 AliCaloRawAnalyzerKStandard::AliCaloRawAnalyzerKStandard() : AliCaloRawAnalyzerFitter("Chi Square ( kStandard )", "KStandard")
46 fAlgo = Algo::kStandard;
50 AliCaloRawAnalyzerKStandard::~AliCaloRawAnalyzerKStandard()
57 AliCaloRawAnalyzerKStandard::Evaluate( const vector<AliCaloBunchInfo> &bunchlist, const UInt_t altrocfg1, const UInt_t altrocfg2 )
59 Float_t pedEstimate = 0;
64 Float_t ampEstimate = 0;
65 short timeEstimate = 0;
70 Bool_t fitDone = kFALSE;
71 int nsamples = PreFitEvaluateSamples( bunchlist, altrocfg1, altrocfg2, bunchIndex, ampEstimate,
72 maxADC, timeEstimate, pedEstimate, first, last, (int)fAmpCut );
75 if (ampEstimate >= fAmpCut )
78 Int_t timebinOffset = bunchlist.at(bunchIndex).GetStartBin() - (bunchlist.at(bunchIndex).GetLength()-1);
81 if ( nsamples > 1 && maxADC< OVERFLOWCUT )
83 FitRaw(first, last, amp, time, chi2, fitDone);
84 time += timebinOffset;
85 timeEstimate += timebinOffset;
91 Float_t ampAsymm = (amp - ampEstimate)/(amp + ampEstimate);
92 Float_t timeDiff = time - timeEstimate;
94 if ( (TMath::Abs(ampAsymm) > 0.1) || (TMath::Abs(timeDiff) > 2) )
105 amp += (0.5 - gRandom->Rndm());
107 time = time * TIMEBINWITH;
110 return AliCaloFitResults( -99, -99, fAlgo , amp, time,
111 (int)time, chi2, ndf, Ret::kDummy );
113 return AliCaloFitResults( Ret::kInvalid, Ret::kInvalid );
117 //____________________________________________________________________________
119 AliCaloRawAnalyzerKStandard::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Float_t & chi2, Bool_t & fitDone) const
121 // Fits the raw signal time distribution
122 int nsamples = lastTimeBin - firstTimeBin + 1;
124 if (nsamples < 3) { return; }
126 TGraph *gSig = new TGraph( nsamples);
128 for (int i=0; i<nsamples; i++)
130 Int_t timebin = firstTimeBin + i;
131 gSig->SetPoint(i, timebin, GetReversed(timebin));
134 TF1 * signalF = new TF1("signal", RawResponseFunction, 0, TIMEBINS , 5);
135 signalF->SetParameters(10.,5., TAU ,ORDER,0.); //set all defaults once, just to be safe
136 signalF->SetParNames("amp","t0","tau","N","ped");
137 signalF->FixParameter(2,TAU);
138 signalF->FixParameter(3,ORDER);
139 signalF->FixParameter(4, 0);
140 signalF->SetParameter(1, time);
141 signalF->SetParameter(0, amp);
142 signalF->SetParLimits(0, 0.5*amp, 2*amp );
143 signalF->SetParLimits(1, time - 4, time + 4);
146 gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points
147 amp = signalF->GetParameter(0);
148 time = signalF->GetParameter(1);
149 chi2 = signalF->GetChisquare();
152 catch (const std::exception & e) {
153 AliError( Form("TGraph Fit exception %s", e.what()) );
154 // stay with default amp and time in case of exception, i.e. no special action required
159 delete gSig; // delete TGraph
164 //__________________________________________________________________
166 AliCaloRawAnalyzerKStandard::FitParabola(const TGraph *gSig, Float_t & amp) const
168 //BEG YS alternative methods to calculate the amplitude
169 Double_t * ymx = gSig->GetX() ;
170 Double_t * ymy = gSig->GetY() ;
172 Double_t ymMaxX[kN] = {0., 0., 0.} ;
173 Double_t ymMaxY[kN] = {0., 0., 0.} ;
175 // find the maximum amplitude
177 for (Int_t ymi = 0; ymi < gSig->GetN(); ymi++) {
178 if (ymy[ymi] > ymMaxY[0] ) {
179 ymMaxY[0] = ymy[ymi] ; //<========== This is the maximum amplitude
180 ymMaxX[0] = ymx[ymi] ;
184 // find the maximum by fitting a parabola through the max and the two adjacent samples
185 if ( ymiMax < gSig->GetN()-1 && ymiMax > 0) {
186 ymMaxY[1] = ymy[ymiMax+1] ;
187 ymMaxY[2] = ymy[ymiMax-1] ;
188 ymMaxX[1] = ymx[ymiMax+1] ;
189 ymMaxX[2] = ymx[ymiMax-1] ;
190 if (ymMaxY[0]*ymMaxY[1]*ymMaxY[2] > 0) {
191 //fit a parabola through the 3 points y= a+bx+x*x*x
199 for (Int_t i = 0; i < kN ; i++) {
202 sx2 += ymMaxX[i]*ymMaxX[i] ;
203 sx3 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ;
204 sx4 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ;
205 sxy += ymMaxX[i]*ymMaxY[i] ;
206 sx2y += ymMaxX[i]*ymMaxX[i]*ymMaxY[i] ;
208 Double_t cN = (sx2y*kN-sy*sx2)*(sx3*sx-sx2*sx2)-(sx2y*sx-sxy*sx2)*(sx3*kN-sx*sx2);
209 Double_t cD = (sx4*kN-sx2*sx2)*(sx3*sx-sx2*sx2)-(sx4*sx-sx3*sx2)*(sx3*kN-sx*sx2) ;
210 Double_t c = cN / cD ;
211 Double_t b = ((sx2y*kN-sy*sx2)-c*(sx4*kN-sx2*sx2))/(sx3*kN-sx*sx2) ;
212 Double_t a = (sy-b*sx-c*sx2)/kN ;
213 Double_t xmax = -b/(2*c) ;
214 ymax = a + b*xmax + c*xmax*xmax ;//<========== This is the maximum amplitude
219 Double_t diff = TMath::Abs(1-ymMaxY[0]/amp) ;
229 //__________________________________________________________________
231 AliCaloRawAnalyzerKStandard::RawResponseFunction(Double_t *x, Double_t *par)
233 Double_t signal = 0.;
234 Double_t tau = par[2];
236 Double_t ped = par[4];
237 Double_t xx = ( x[0] - par[1] + tau ) / tau ;
242 signal = ped + par[0] * TMath::Power(xx , n) * TMath::Exp(n * (1 - xx )) ;