]>
Commit | Line | Data |
---|---|---|
baff65a9 | 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 decodes the stream of ALTRO samples to extract | |
19 | // the PHOS "digits" of current event. Uses fitting procedure | |
20 | // to separate reasonable samples | |
21 | // | |
22 | // Typical use case: | |
23 | // AliRawReader* rf = new AliRawReaderDate("2006run2211.raw"); | |
24 | // AliPHOSRawDecoder dc(rf); | |
25 | // while (rf->NextEvent()) { | |
baff65a9 | 26 | // dc.SubtractPedestals(kTRUE); |
27 | // while ( dc.NextDigit() ) { | |
28 | // Int_t module = dc.GetModule(); | |
29 | // Int_t column = dc.GetColumn(); | |
30 | // Int_t row = dc.GetRow(); | |
31 | // Double_t amplitude = dc.GetEnergy(); | |
32 | // Double_t time = dc.GetTime(); | |
33 | // Bool_t IsLowGain = dc.IsLowGain(); | |
34 | // .......... | |
35 | // } | |
36 | // } | |
37 | ||
38 | // Author: Dmitri Peressounko using idea of Y.Kucheryaev | |
39 | ||
40 | // --- ROOT system --- | |
41 | #include "TList.h" | |
42 | #include "TMath.h" | |
43 | #include "TMinuit.h" | |
44 | ||
45 | #include "TCanvas.h" | |
46 | #include "TH1.h" | |
47 | #include "TH2.h" | |
48 | #include "TF1.h" | |
49 | #include "TROOT.h" | |
50 | ||
51 | // --- AliRoot header files --- | |
52 | //#include "AliLog.h" | |
53 | #include "AliPHOSRawDecoderv2.h" | |
54 | #include "AliPHOSPulseGenerator.h" | |
55 | ||
56 | ||
57 | ClassImp(AliPHOSRawDecoderv2) | |
58 | ||
59 | //----------------------------------------------------------------------------- | |
60 | AliPHOSRawDecoderv2::AliPHOSRawDecoderv2():AliPHOSRawDecoder(), | |
61 | fNtimeSamples(25),fRMScut(11.) | |
62 | { | |
63 | //Default constructor. | |
64 | fLGpar[0]=0.971 ; | |
65 | fLGpar[1]=0.0465; | |
66 | fLGpar[2]=1.56 ; | |
67 | fHGpar[0]=0.941 ; | |
68 | fHGpar[1]=0.0436; | |
69 | fHGpar[2]=1.65 ; | |
70 | } | |
71 | ||
72 | //----------------------------------------------------------------------------- | |
73 | AliPHOSRawDecoderv2::AliPHOSRawDecoderv2(AliRawReader* rawReader, AliAltroMapping **mapping): | |
74 | AliPHOSRawDecoder(rawReader,mapping), | |
75 | fNtimeSamples(25),fRMScut(11.) | |
76 | { | |
77 | //Construct a decoder object. | |
78 | //Is is user responsibility to provide next raw event | |
79 | //using AliRawReader::NextEvent(). | |
80 | fLGpar[0]=0.971 ; | |
81 | fLGpar[1]=0.0465; | |
82 | fLGpar[2]=1.56 ; | |
83 | fHGpar[0]=0.941 ; | |
84 | fHGpar[1]=0.0436; | |
85 | fHGpar[2]=1.65 ; | |
86 | } | |
87 | ||
88 | //----------------------------------------------------------------------------- | |
89 | AliPHOSRawDecoderv2::~AliPHOSRawDecoderv2() | |
90 | { | |
91 | //Destructor. | |
92 | //Nothing to delete | |
93 | } | |
94 | ||
95 | //----------------------------------------------------------------------------- | |
96 | AliPHOSRawDecoderv2::AliPHOSRawDecoderv2(const AliPHOSRawDecoderv2 &phosDecoder ): | |
97 | AliPHOSRawDecoder(phosDecoder), | |
98 | fNtimeSamples(25),fRMScut(11.) | |
99 | { | |
100 | //Copy constructor. | |
101 | fNtimeSamples=phosDecoder.fNtimeSamples ; | |
102 | for(Int_t i=0; i<3;i++){ | |
103 | fLGpar[i]=phosDecoder.fLGpar[i] ; | |
104 | fHGpar[i]=phosDecoder.fHGpar[i] ; | |
105 | } | |
106 | fRMScut=phosDecoder.fRMScut ; | |
107 | } | |
108 | ||
109 | //----------------------------------------------------------------------------- | |
110 | AliPHOSRawDecoderv2& AliPHOSRawDecoderv2::operator = (const AliPHOSRawDecoderv2 &phosDecoder) | |
111 | { | |
112 | //Assignment operator. | |
113 | ||
114 | fNtimeSamples=phosDecoder.fNtimeSamples ; | |
115 | for(Int_t i=0; i<3;i++){ | |
116 | fLGpar[i]=phosDecoder.fLGpar[i] ; | |
117 | fHGpar[i]=phosDecoder.fHGpar[i] ; | |
118 | } | |
119 | fRMScut=phosDecoder.fRMScut ; | |
120 | return *this; | |
121 | } | |
122 | ||
123 | //----------------------------------------------------------------------------- | |
124 | Bool_t AliPHOSRawDecoderv2::NextDigit() | |
125 | { | |
126 | //Extract an energy deposited in the crystal, | |
127 | //crystal' position (module,column,row), | |
128 | //time and gain (high or low). | |
129 | //First collects sample, then evaluates it and if it has | |
130 | //reasonable shape, fits it with Gamma2 function and extracts | |
131 | //energy and time. | |
132 | ||
39569d36 | 133 | TCanvas * cs = (TCanvas*)gROOT->FindObjectAny("CSample") ; |
134 | if(!cs) | |
135 | cs = new TCanvas("CSample","CSample") ; | |
136 | ||
137 | TH1D * h = (TH1D*)gROOT->FindObjectAny("hSample") ; | |
138 | if(!h) | |
139 | h=new TH1D("hSample","",200,0.,200.) ; | |
baff65a9 | 140 | |
141 | ||
142 | AliCaloRawStream* in = fCaloStream; | |
143 | ||
144 | Int_t iBin = fSamples->GetSize(); | |
145 | fEnergy = 0; | |
146 | Double_t pedMean = 0; | |
39569d36 | 147 | Double_t pedRMS = 0; |
baff65a9 | 148 | Int_t nPed = 0; |
149 | Double_t baseLine = 1.0; | |
150 | const Int_t nPreSamples = 10; | |
39569d36 | 151 | fQuality = 0. ; |
baff65a9 | 152 | |
153 | while ( in->Next() ) { | |
154 | ||
155 | // Fit the full sample | |
156 | if(in->IsNewHWAddress() && iBin!=fSamples->GetSize()) { | |
157 | ||
158 | Double_t pedestal =0. ; | |
159 | if(fPedSubtract){ | |
160 | if (nPed > 0) | |
161 | pedestal = (Double_t)(pedMean/nPed); | |
162 | else | |
163 | return kFALSE; | |
164 | } | |
39569d36 | 165 | for(Int_t i=0; i<fSamples->GetSize(); i++){ |
166 | h->SetBinContent(i+1,fSamples->At(i)) ; | |
167 | } | |
baff65a9 | 168 | |
169 | //calculate time and energy | |
170 | Int_t maxBin=0 ; | |
171 | Int_t maxAmp=0 ; | |
172 | for(Int_t i=iBin; i<fSamples->GetSize(); i++){ | |
39569d36 | 173 | if(maxAmp<fSamples->At(i)){ |
baff65a9 | 174 | maxBin=i ; |
175 | maxAmp=fSamples->At(i) ; | |
176 | } | |
177 | } | |
39569d36 | 178 | if(maxBin==fSamples->GetSize()-1){//bad sample |
179 | fEnergy=0. ; | |
180 | fTime=-999.; | |
181 | return kTRUE ; | |
182 | } | |
183 | fEnergy=Double_t(maxAmp)-pedestal ; | |
baff65a9 | 184 | fOverflow =0 ; //look for plato on the top of sample |
39569d36 | 185 | if(fEnergy>500 && //this is not fluctuation of soft sample |
baff65a9 | 186 | maxBin<fSamples->GetSize()-1 && fSamples->At(maxBin+1)==maxAmp){ //and there is a plato |
187 | fOverflow = kTRUE ; | |
188 | } | |
189 | ||
39569d36 | 190 | // if(fEnergy>500.){ |
191 | // if(fRow==54 && fColumn==24){ | |
192 | // printf("fE=%f, ped=%f, row=%d, col=%d \n",fEnergy,pedestal,fRow,fColumn) ; | |
193 | // if(fOverflow)printf(" Overflow \n") ; | |
194 | // else printf("iBin=%d, maxBin=%d, maxAmp=%d,Amp(+1)=%d,Amp(-1)=%d \n",iBin,maxBin,maxAmp,fSamples->At(maxBin+1),fSamples->At(maxBin-1)) ; | |
baff65a9 | 195 | // cs->cd() ; |
196 | // h->Draw() ; | |
197 | // cs->Update() ; | |
198 | // getchar() ; | |
39569d36 | 199 | // } |
baff65a9 | 200 | |
a28333c5 | 201 | if(fOverflow) |
202 | return kTRUE ; //do not calculate energy and time for overflowed channels | |
203 | ||
baff65a9 | 204 | if(fEnergy<baseLine){ //do not evaluate time, drop this sample |
205 | fEnergy=0. ; | |
206 | fTime=-999.; | |
207 | return kTRUE ; | |
208 | } | |
209 | ||
210 | //else calculate time | |
211 | fTime=0. ; | |
212 | Double_t tRMS = 0. ; | |
213 | Double_t tW = 0. ; | |
214 | Int_t cnts=0 ; | |
215 | Double_t a,b,c ; | |
216 | if(fLowGainFlag){ | |
217 | a=fLGpar[0] ; | |
218 | b=fLGpar[1] ; | |
219 | c=fLGpar[2] ; | |
220 | } | |
221 | else{ | |
222 | a=fHGpar[0] ; | |
223 | b=fHGpar[1] ; | |
224 | c=fHGpar[2] ; | |
225 | } | |
226 | ||
227 | for(Int_t i=iBin+1; i<fSamples->GetSize()&& cnts<fNtimeSamples; i++){ | |
228 | if(fSamples->At(i)<pedestal) | |
229 | continue ; | |
230 | //Presently we do not use channels with overflow | |
231 | // if(fOverflow && (fSamples->At(i)==maxAmp ||fSamples->At(i-1)==maxAmp)) //can not calculate time in overflow bin | |
232 | // continue ; | |
233 | if(fTimes->At(i)-fTimes->At(i-1)!=1) //do not use samples with non-consequtive points | |
234 | continue ; | |
235 | Double_t de=fSamples->At(i)-pedestal ; | |
236 | Double_t av = de+fSamples->At(i-1)-pedestal ; | |
237 | if(av<=0.) //this is fluctuation around pedestal, scip | |
238 | continue ; | |
239 | Double_t ds = fSamples->At(i)-fSamples->At(i-1) ; | |
240 | Double_t ti = ds/av ; // calculate log. derivative | |
241 | ti=a/(ti+b)-c*ti ; // and compare with parameterization | |
242 | ti=Double_t(fTimes->At(i))-ti ; | |
243 | Double_t wi = TMath::Abs(ds) ; | |
244 | fTime+=ti*wi ; | |
245 | tW+=wi; | |
246 | tRMS+=ti*ti*wi ; | |
247 | cnts++ ; | |
248 | } | |
249 | ||
250 | if(tW>0.){ | |
251 | fTime/=tW ; | |
39569d36 | 252 | fQuality = tRMS/tW-fTime*fTime ; |
253 | //Normalize quality | |
baff65a9 | 254 | //printf("t0=%f, RMS=%f, cut=%f \n",fTime,tRMS,fRMScut) ; |
39569d36 | 255 | // if(tRMS>=fRMScut){ //bad sample |
256 | // fTime=-999. ; | |
257 | // fEnergy=0. ; | |
258 | // } | |
baff65a9 | 259 | } |
260 | else{ | |
261 | fTime=-999. ; | |
39569d36 | 262 | fQuality=999. ; |
263 | } | |
264 | ||
265 | Bool_t isBad = 0 ; | |
266 | for(Int_t i=iBin+1; i<fSamples->GetSize()-1&&!isBad; i++){ | |
267 | if(fSamples->At(i)>fSamples->At(i-1)+5 && fSamples->At(i)>fSamples->At(i+1)+5) { //single jump | |
268 | isBad=1 ; | |
269 | } | |
baff65a9 | 270 | } |
39569d36 | 271 | if(pedestal<10.) |
272 | isBad=1 ; | |
273 | ||
274 | pedRMS=pedRMS/nPed-pedestal*pedestal ; | |
275 | if(pedRMS>0.1) | |
276 | isBad=1 ; | |
277 | ||
278 | for(Int_t i=iBin+1; i<fSamples->GetSize()-1&&!isBad; i++){ | |
279 | if(fSamples->At(i)<pedestal-1) | |
280 | isBad=1 ; | |
281 | } | |
282 | ||
283 | //two maxima | |
284 | ||
285 | ||
286 | if(fEnergy>10. && !isBad ){ | |
287 | printf("fE=%f, ped=%f, fQuality=%f, pedRMS=%f \n",fEnergy,pedestal,fQuality,pedRMS) ; | |
288 | if(fOverflow)printf(" Overflow \n") ; | |
289 | if(isBad)printf("bad") ; | |
290 | // else printf("iBin=%d, maxBin=%d, maxAmp=%d,Amp(+1)=%d,Amp(-1)=%d \n",iBin,maxBin,maxAmp,fSamples->At(maxBin+1),fSamples->At(maxBin-1)) ; | |
291 | cs->cd() ; | |
292 | h->Draw() ; | |
293 | cs->Update() ; | |
294 | getchar() ; | |
295 | } | |
baff65a9 | 296 | |
baff65a9 | 297 | |
298 | return kTRUE ; //will scan further | |
299 | } | |
300 | ||
301 | fLowGainFlag = in->IsLowGain(); | |
302 | fModule = in->GetModule()+1; | |
303 | fRow = in->GetRow() +1; | |
304 | fColumn = in->GetColumn()+1; | |
305 | ||
306 | ||
307 | // Fill array with samples | |
308 | iBin--; | |
309 | if(fPedSubtract && (in->GetTime() < nPreSamples)) { | |
310 | pedMean += in->GetSignal(); | |
39569d36 | 311 | pedRMS += in->GetSignal()*in->GetSignal(); |
baff65a9 | 312 | nPed++; |
313 | } | |
314 | fSamples->AddAt(in->GetSignal(),iBin); | |
315 | fTimes->AddAt(in->GetTime(),iBin); | |
316 | } // in.Next() | |
317 | ||
318 | return kFALSE; | |
319 | } |