Line data Source code
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 :
19 : //
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 :
38 : #include "AliEMCALRawResponse.h"
39 :
40 :
41 : using namespace std;
42 :
43 42 : ClassImp( AliCaloRawAnalyzerKStandard )
44 :
45 :
46 12 : AliCaloRawAnalyzerKStandard::AliCaloRawAnalyzerKStandard() : AliCaloRawAnalyzerFitter("Chi Square ( kStandard )", "KStandard")
47 60 : {
48 :
49 12 : fAlgo = Algo::kStandard;
50 24 : }
51 :
52 :
53 : AliCaloRawAnalyzerKStandard::~AliCaloRawAnalyzerKStandard()
54 40 : {
55 : // delete fTf1;
56 40 : }
57 :
58 :
59 : AliCaloFitResults
60 : AliCaloRawAnalyzerKStandard::Evaluate( const vector<AliCaloBunchInfo> &bunchlist, UInt_t altrocfg1, UInt_t altrocfg2 )
61 : {
62 : //Evaluation Amplitude and TOF
63 116 : Float_t pedEstimate = 0;
64 58 : short maxADC = 0;
65 58 : Int_t first = 0;
66 58 : Int_t last = 0;
67 58 : Int_t bunchIndex = 0;
68 58 : Float_t ampEstimate = 0;
69 58 : short timeEstimate = 0;
70 58 : Float_t time = 0;
71 58 : Float_t amp=0;
72 58 : Float_t chi2 = 0;
73 : Int_t ndf = 0;
74 58 : Bool_t fitDone = kFALSE;
75 116 : int nsamples = PreFitEvaluateSamples( bunchlist, altrocfg1, altrocfg2, bunchIndex, ampEstimate,
76 58 : maxADC, timeEstimate, pedEstimate, first, last, (int)fAmpCut );
77 :
78 :
79 58 : if (ampEstimate >= fAmpCut )
80 : {
81 58 : time = timeEstimate;
82 58 : Int_t timebinOffset = bunchlist.at(bunchIndex).GetStartBin() - (bunchlist.at(bunchIndex).GetLength()-1);
83 58 : amp = ampEstimate;
84 :
85 116 : if ( nsamples > 1 && maxADC< OVERFLOWCUT )
86 : {
87 58 : FitRaw(first, last, amp, time, chi2, fitDone);
88 58 : time += timebinOffset;
89 58 : timeEstimate += timebinOffset;
90 58 : ndf = nsamples - 2;
91 58 : }
92 58 : }
93 58 : if ( fitDone )
94 : {
95 58 : Float_t ampAsymm = (amp - ampEstimate)/(amp + ampEstimate);
96 58 : Float_t timeDiff = time - timeEstimate;
97 :
98 116 : if ( (TMath::Abs(ampAsymm) > 0.1) || (TMath::Abs(timeDiff) > 2) )
99 : {
100 0 : amp = ampEstimate;
101 0 : time = timeEstimate;
102 0 : fitDone = kFALSE;
103 0 : }
104 58 : }
105 58 : if (amp >= fAmpCut )
106 : {
107 58 : if ( ! fitDone)
108 : {
109 0 : amp += (0.5 - gRandom->Rndm());
110 0 : }
111 58 : time = time * TIMEBINWITH;
112 58 : time -= fL1Phase;
113 :
114 116 : return AliCaloFitResults( -99, -99, fAlgo , amp, time,
115 58 : (int)time, chi2, ndf, Ret::kDummy );
116 : }
117 0 : return AliCaloFitResults( Ret::kInvalid, Ret::kInvalid );
118 58 : }
119 :
120 :
121 : void
122 : AliCaloRawAnalyzerKStandard::FitRaw( Int_t firstTimeBin, Int_t lastTimeBin,
123 : Float_t & amp, Float_t & time, Float_t & chi2, Bool_t & fitDone) const
124 : {
125 : // Fits the raw signal time distribution
126 116 : int nsamples = lastTimeBin - firstTimeBin + 1;
127 58 : fitDone = kFALSE;
128 58 : if (nsamples < 3) { return; }
129 :
130 58 : TGraph *gSig = new TGraph( nsamples);
131 :
132 954 : for (int i=0; i<nsamples; i++)
133 : {
134 419 : Int_t timebin = firstTimeBin + i;
135 419 : gSig->SetPoint(i, timebin, GetReversed(timebin));
136 : }
137 :
138 58 : TF1 * signalF = new TF1("signal", AliEMCALRawResponse::RawResponseFunction, 0, TIMEBINS , 5);
139 :
140 58 : signalF->SetParameters(10.,5., TAU ,ORDER,0.); //set all defaults once, just to be safe
141 58 : signalF->SetParNames("amp","t0","tau","N","ped");
142 58 : signalF->FixParameter(2,TAU);
143 58 : signalF->FixParameter(3,ORDER);
144 58 : signalF->FixParameter(4, 0);
145 58 : signalF->SetParameter(1, time);
146 58 : signalF->SetParameter(0, amp);
147 58 : signalF->SetParLimits(0, 0.5*amp, 2*amp );
148 58 : signalF->SetParLimits(1, time - 4, time + 4);
149 :
150 : try {
151 116 : gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points
152 116 : amp = signalF->GetParameter(0);
153 116 : time = signalF->GetParameter(1);
154 58 : chi2 = signalF->GetChisquare();
155 58 : fitDone = kTRUE;
156 58 : }
157 : catch (const std::exception & e)
158 : {
159 0 : AliError( Form("TGraph Fit exception %s", e.what()) );
160 : // stay with default amp and time in case of exception, i.e. no special action required
161 0 : fitDone = kFALSE;
162 0 : }
163 :
164 116 : delete signalF;
165 116 : delete gSig; // delete TGraph
166 : return;
167 58 : }
168 :
169 :
|
