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