]>
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 | ||
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.) ; | |
141 | ||
142 | ||
143 | AliCaloRawStream* in = fCaloStream; | |
144 | ||
145 | Int_t iBin = fSamples->GetSize(); | |
146 | fEnergy = 0; | |
147 | Double_t pedMean = 0; | |
148 | Int_t nPed = 0; | |
149 | Double_t baseLine = 1.0; | |
150 | const Int_t nPreSamples = 10; | |
151 | ||
152 | while ( in->Next() ) { | |
153 | ||
154 | // Fit the full sample | |
155 | if(in->IsNewHWAddress() && iBin!=fSamples->GetSize()) { | |
156 | ||
157 | Double_t pedestal =0. ; | |
158 | if(fPedSubtract){ | |
159 | if (nPed > 0) | |
160 | pedestal = (Double_t)(pedMean/nPed); | |
161 | else | |
162 | return kFALSE; | |
163 | } | |
164 | // for(Int_t i=0; i<fSamples->GetSize(); i++){ | |
165 | // h->SetBinContent(i+1,fSamples->At(i)) ; | |
166 | // } | |
167 | ||
168 | //calculate time and energy | |
169 | Int_t maxBin=0 ; | |
170 | Int_t maxAmp=0 ; | |
171 | for(Int_t i=iBin; i<fSamples->GetSize(); i++){ | |
172 | Double_t de=fSamples->At(i)-pedestal ; | |
173 | if(fEnergy<de){ | |
174 | fEnergy=de ; | |
175 | maxBin=i ; | |
176 | maxAmp=fSamples->At(i) ; | |
177 | } | |
178 | } | |
179 | fOverflow =0 ; //look for plato on the top of sample | |
180 | if(fEnergy>700 && //this is not fluctuation of soft sample | |
181 | maxBin<fSamples->GetSize()-1 && fSamples->At(maxBin+1)==maxAmp){ //and there is a plato | |
182 | fOverflow = kTRUE ; | |
183 | } | |
184 | ||
a28333c5 | 185 | // if(fEnergy>1020.){ |
baff65a9 | 186 | // printf("fE=%f \n",fEnergy) ; |
baff65a9 | 187 | // cs->cd() ; |
188 | // h->Draw() ; | |
189 | // cs->Update() ; | |
190 | // getchar() ; | |
191 | // } | |
192 | ||
a28333c5 | 193 | if(fOverflow) |
194 | return kTRUE ; //do not calculate energy and time for overflowed channels | |
195 | ||
baff65a9 | 196 | if(fEnergy<baseLine){ //do not evaluate time, drop this sample |
197 | fEnergy=0. ; | |
198 | fTime=-999.; | |
199 | return kTRUE ; | |
200 | } | |
201 | ||
202 | //else calculate time | |
203 | fTime=0. ; | |
204 | Double_t tRMS = 0. ; | |
205 | Double_t tW = 0. ; | |
206 | Int_t cnts=0 ; | |
207 | Double_t a,b,c ; | |
208 | if(fLowGainFlag){ | |
209 | a=fLGpar[0] ; | |
210 | b=fLGpar[1] ; | |
211 | c=fLGpar[2] ; | |
212 | } | |
213 | else{ | |
214 | a=fHGpar[0] ; | |
215 | b=fHGpar[1] ; | |
216 | c=fHGpar[2] ; | |
217 | } | |
218 | ||
219 | for(Int_t i=iBin+1; i<fSamples->GetSize()&& cnts<fNtimeSamples; i++){ | |
220 | if(fSamples->At(i)<pedestal) | |
221 | continue ; | |
222 | //Presently we do not use channels with overflow | |
223 | // if(fOverflow && (fSamples->At(i)==maxAmp ||fSamples->At(i-1)==maxAmp)) //can not calculate time in overflow bin | |
224 | // continue ; | |
225 | if(fTimes->At(i)-fTimes->At(i-1)!=1) //do not use samples with non-consequtive points | |
226 | continue ; | |
227 | Double_t de=fSamples->At(i)-pedestal ; | |
228 | Double_t av = de+fSamples->At(i-1)-pedestal ; | |
229 | if(av<=0.) //this is fluctuation around pedestal, scip | |
230 | continue ; | |
231 | Double_t ds = fSamples->At(i)-fSamples->At(i-1) ; | |
232 | Double_t ti = ds/av ; // calculate log. derivative | |
233 | ti=a/(ti+b)-c*ti ; // and compare with parameterization | |
234 | ti=Double_t(fTimes->At(i))-ti ; | |
235 | Double_t wi = TMath::Abs(ds) ; | |
236 | fTime+=ti*wi ; | |
237 | tW+=wi; | |
238 | tRMS+=ti*ti*wi ; | |
239 | cnts++ ; | |
240 | } | |
241 | ||
242 | if(tW>0.){ | |
243 | fTime/=tW ; | |
244 | tRMS/=tW ; | |
245 | tRMS-=fTime*fTime ; | |
246 | //printf("t0=%f, RMS=%f, cut=%f \n",fTime,tRMS,fRMScut) ; | |
247 | if(tRMS>=fRMScut){ //bad sample | |
248 | fTime=-999. ; | |
249 | fEnergy=0. ; | |
250 | } | |
251 | } | |
252 | else{ | |
253 | fTime=-999. ; | |
254 | } | |
baff65a9 | 255 | |
256 | fTime*=fPulseGenerator->GetRawFormatTimeTrigger() ; | |
257 | ||
258 | return kTRUE ; //will scan further | |
259 | } | |
260 | ||
261 | fLowGainFlag = in->IsLowGain(); | |
262 | fModule = in->GetModule()+1; | |
263 | fRow = in->GetRow() +1; | |
264 | fColumn = in->GetColumn()+1; | |
265 | ||
266 | ||
267 | // Fill array with samples | |
268 | iBin--; | |
269 | if(fPedSubtract && (in->GetTime() < nPreSamples)) { | |
270 | pedMean += in->GetSignal(); | |
271 | nPed++; | |
272 | } | |
273 | fSamples->AddAt(in->GetSignal(),iBin); | |
274 | fTimes->AddAt(in->GetTime(),iBin); | |
275 | } // in.Next() | |
276 | ||
277 | return kFALSE; | |
278 | } |