]>
Commit | Line | Data |
---|---|---|
379c5c09 | 1 | /************************************************************************** |
2 | * Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /* $Id: $ */ | |
17 | ||
18 | // This class plots samples and qualify their quality according to their shape | |
19 | // | |
20 | // Typical use case: | |
21 | // AliPHOSRawFitter *fitter=new AliPHOSRawFitter(); | |
379c5c09 | 22 | // fitter->SetChannelGeo(module,cellX,cellZ,caloFlag); |
23 | // fitter->SetCalibData(fgCalibData) ; | |
1dfadc32 | 24 | // fitter->Eval(sig,sigStart,sigLength); |
379c5c09 | 25 | // Double_t amplitude = fitter.GetEnergy(); |
26 | // Double_t time = fitter.GetTime(); | |
27 | // Bool_t isLowGain = fitter.GetCaloFlag()==0; | |
28 | ||
29 | // Author: Dmitri Peressounko (Oct.2008) | |
30 | // Modified: Yuri Kharlov (Jul.2009) | |
31 | ||
32 | // --- ROOT system --- | |
33 | #include "TList.h" | |
34 | #include "TMath.h" | |
35 | #include "TMinuit.h" | |
36 | #include "TCanvas.h" | |
37 | #include "TH1.h" | |
38 | #include "TH2.h" | |
39 | #include "TF1.h" | |
40 | #include "TROOT.h" | |
41 | ||
42 | // --- AliRoot header files --- | |
43 | #include "AliLog.h" | |
44 | #include "AliPHOSCalibData.h" | |
45 | #include "AliPHOSRawFitterv2.h" | |
46 | #include "AliPHOSPulseGenerator.h" | |
47 | ||
48 | ClassImp(AliPHOSRawFitterv2) | |
49 | ||
50 | //----------------------------------------------------------------------------- | |
51 | AliPHOSRawFitterv2::AliPHOSRawFitterv2(): | |
52 | AliPHOSRawFitterv0(), | |
53 | fNtimeSamples(25), | |
54 | fRMScut(11.) | |
55 | { | |
56 | //Default constructor. | |
57 | fLGpar[0]=0.971 ; | |
58 | fLGpar[1]=0.0465; | |
59 | fLGpar[2]=1.56 ; | |
60 | fHGpar[0]=0.941 ; | |
61 | fHGpar[1]=0.0436; | |
62 | fHGpar[2]=1.65 ; | |
63 | } | |
64 | ||
65 | //----------------------------------------------------------------------------- | |
66 | AliPHOSRawFitterv2::~AliPHOSRawFitterv2() | |
67 | { | |
68 | //Destructor. | |
69 | } | |
70 | ||
71 | //----------------------------------------------------------------------------- | |
72 | AliPHOSRawFitterv2::AliPHOSRawFitterv2(const AliPHOSRawFitterv2 &phosFitter ): | |
73 | AliPHOSRawFitterv0(phosFitter), | |
74 | fNtimeSamples(25), | |
75 | fRMScut(11.) | |
76 | { | |
77 | //Copy constructor. | |
78 | fNtimeSamples=phosFitter.fNtimeSamples ; | |
79 | for(Int_t i=0; i<3;i++){ | |
80 | fLGpar[i]=phosFitter.fLGpar[i] ; | |
81 | fHGpar[i]=phosFitter.fHGpar[i] ; | |
82 | } | |
83 | fRMScut=phosFitter.fRMScut ; | |
84 | } | |
85 | ||
86 | //----------------------------------------------------------------------------- | |
87 | AliPHOSRawFitterv2& AliPHOSRawFitterv2::operator = (const AliPHOSRawFitterv2 &phosFitter) | |
88 | { | |
89 | //Assignment operator. | |
90 | ||
91 | fNtimeSamples=phosFitter.fNtimeSamples ; | |
92 | for(Int_t i=0; i<3;i++){ | |
93 | fLGpar[i]=phosFitter.fLGpar[i] ; | |
94 | fHGpar[i]=phosFitter.fHGpar[i] ; | |
95 | } | |
96 | fRMScut=phosFitter.fRMScut ; | |
97 | return *this; | |
98 | } | |
99 | ||
100 | //----------------------------------------------------------------------------- | |
92236b27 | 101 | Bool_t AliPHOSRawFitterv2::Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength) |
379c5c09 | 102 | { |
103 | //Extract an energy deposited in the crystal, | |
104 | //crystal' position (module,column,row), | |
105 | //time and gain (high or low). | |
106 | //First collects sample, then evaluates it and if it has | |
107 | //reasonable shape, fits it with Gamma2 function and extracts | |
108 | //energy and time. | |
109 | ||
110 | if (fNBunches > 1) { | |
111 | fQuality = 1000; | |
112 | return kTRUE; | |
113 | } | |
114 | ||
115 | const Float_t kBaseLine = 1.0; | |
116 | const Int_t kPreSamples = 10; | |
117 | ||
118 | Float_t pedMean = 0; | |
119 | Float_t pedRMS = 0; | |
120 | Int_t nPed = 0; | |
121 | UShort_t maxSample = 0; | |
122 | Int_t nMax = 0; | |
123 | ||
124 | TCanvas * cs = (TCanvas*)gROOT->FindObjectAny("CSample") ; | |
125 | if(!cs) | |
126 | cs = new TCanvas("CSample","CSample") ; | |
127 | ||
128 | TH1D * h = (TH1D*)gROOT->FindObjectAny("hSample") ; | |
129 | if(!h) h = new TH1D("hSample","",200,0.,200.) ; | |
130 | ||
131 | Double_t pedestal = 0; | |
92236b27 | 132 | for (Int_t i=0; i<sigLength; i++) { |
379c5c09 | 133 | if (i<kPreSamples) { |
134 | nPed++; | |
92236b27 | 135 | pedMean += signal[i]; |
136 | pedRMS += signal[i]*signal[i] ; | |
379c5c09 | 137 | } |
92236b27 | 138 | if(signal[i] > maxSample) maxSample = signal[i]; |
139 | if(signal[i] == maxSample) nMax++; | |
140 | h->SetBinContent(i+1,signal[i]) ; | |
379c5c09 | 141 | } |
142 | fEnergy = (Double_t)maxSample; | |
143 | if (maxSample > 900 && nMax > 2) fOverflow = kTRUE; | |
144 | ||
145 | if (fPedSubtract) { | |
146 | if (nPed > 0) { | |
147 | fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ; | |
148 | if(fPedestalRMS > 0.) | |
149 | fPedestalRMS = TMath::Sqrt(fPedestalRMS) ; | |
150 | fEnergy -= (Double_t)(pedMean/nPed); // pedestal subtraction | |
151 | } | |
152 | else | |
153 | return kFALSE; | |
154 | } | |
155 | else { | |
156 | //take pedestals from DB | |
157 | pedestal = (Double_t) fAmpOffset ; | |
158 | if (fCalibData) { | |
159 | Float_t truePed = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ; | |
160 | Int_t altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ; | |
161 | pedestal += truePed - altroSettings ; | |
162 | } | |
163 | else{ | |
164 | AliWarning(Form("Can not read data from OCDB")) ; | |
165 | } | |
166 | fEnergy-=pedestal ; | |
167 | } | |
168 | ||
169 | if (fEnergy < kBaseLine) { | |
170 | fEnergy = 0; | |
171 | return kTRUE; | |
172 | } | |
173 | ||
174 | // calculate time | |
175 | fTime=0. ; | |
176 | Double_t tRMS = 0. ; | |
177 | Double_t tW = 0. ; | |
178 | Int_t cnts=0 ; | |
179 | Double_t a=0,b=0,c=0 ; | |
180 | if(fCaloFlag == 0){ // Low gain | |
181 | a=fLGpar[0] ; | |
182 | b=fLGpar[1] ; | |
183 | c=fLGpar[2] ; | |
184 | } | |
185 | else if(fCaloFlag == 1){ // High gain | |
186 | a=fHGpar[0] ; | |
187 | b=fHGpar[1] ; | |
188 | c=fHGpar[2] ; | |
189 | } | |
190 | ||
191 | ||
192 | fQuality = 0. ; | |
193 | ||
92236b27 | 194 | for(Int_t i=1; i<sigLength && cnts<fNtimeSamples; i++){ |
195 | if(signal[i] < pedestal) | |
379c5c09 | 196 | continue ; |
92236b27 | 197 | Double_t de = signal[i] - pedestal ; |
198 | Double_t av = signal[i-1] - pedestal + de; | |
379c5c09 | 199 | if(av<=0.) //this is fluctuation around pedestal, skip it |
200 | continue ; | |
92236b27 | 201 | Double_t ds = signal[i] - signal[i-1] ; |
379c5c09 | 202 | Double_t ti = ds/av ; // calculate log. derivative |
203 | ti = a/(ti+b)-c*ti ; // and compare with parameterization | |
204 | ti = i - ti ; | |
205 | Double_t wi = TMath::Abs(ds) ; | |
206 | fTime += ti*wi ; | |
207 | tW += wi; | |
208 | tRMS += ti*ti*wi ; | |
209 | cnts++ ; | |
210 | } | |
211 | ||
212 | if(tW>0.){ | |
213 | fTime/=tW ; | |
214 | fQuality = tRMS/tW-fTime*fTime ; | |
92236b27 | 215 | fTime+=sigStart; |
379c5c09 | 216 | } |
217 | else{ | |
218 | fTime=-999. ; | |
219 | fQuality=999. ; | |
220 | } | |
221 | ||
222 | Bool_t isBad = 0 ; | |
92236b27 | 223 | for(Int_t i=1; i<sigLength-1&&!isBad; i++){ |
224 | if(signal[i] > signal[i-1]+5 && signal[i] > signal[i+1]+5) { //single jump | |
379c5c09 | 225 | isBad=1 ; |
226 | } | |
227 | } | |
228 | if(pedestal < 10.) | |
229 | isBad=1 ; | |
230 | ||
231 | if(fPedestalRMS > 0.1) | |
232 | isBad=1 ; | |
233 | ||
92236b27 | 234 | for(Int_t i=1; i<sigLength-1&&!isBad; i++){ |
235 | if(signal[i] < pedestal-1) | |
379c5c09 | 236 | isBad=1 ; |
237 | } | |
238 | ||
239 | if(fEnergy>10. && !isBad ){ | |
240 | printf("fE=%f, ped=%f, fQuality=%f, pedRMS=%f \n",fEnergy,pedestal,fQuality,pedRMS) ; | |
241 | if(fOverflow)printf(" Overflow \n") ; | |
242 | if(isBad)printf("bad") ; | |
243 | cs->cd() ; | |
244 | h->Draw() ; | |
245 | cs->Update() ; | |
246 | getchar() ; | |
247 | } | |
248 | ||
249 | return kTRUE; | |
250 | } |