]>
Commit | Line | Data |
---|---|---|
77ea1c6f | 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 | |
40 | ||
41 | // --- ROOT system --- | |
11f2ec15 | 42 | #include "TList.h" |
77ea1c6f | 43 | #include "TMath.h" |
44 | #include "TMinuit.h" | |
45 | ||
11f2ec15 | 46 | #include "TCanvas.h" |
47 | #include "TH1.h" | |
48 | #include "TH2.h" | |
49 | #include "TF1.h" | |
50 | #include "TROOT.h" | |
51 | ||
77ea1c6f | 52 | // --- AliRoot header files --- |
11f2ec15 | 53 | //#include "AliLog.h" |
77ea1c6f | 54 | #include "AliPHOSRawDecoderv1.h" |
55 | #include "AliPHOSPulseGenerator.h" | |
56 | ||
57 | ||
58 | ClassImp(AliPHOSRawDecoderv1) | |
59 | ||
60 | //----------------------------------------------------------------------------- | |
11f2ec15 | 61 | AliPHOSRawDecoderv1::AliPHOSRawDecoderv1():AliPHOSRawDecoder(), |
62 | fSampleParamsLow(0x0),fSampleParamsHigh(0x0),fToFit(0x0) | |
77ea1c6f | 63 | { |
64 | //Default constructor. | |
65 | } | |
66 | ||
67 | //----------------------------------------------------------------------------- | |
68 | AliPHOSRawDecoderv1::AliPHOSRawDecoderv1(AliRawReader* rawReader, AliAltroMapping **mapping): | |
11f2ec15 | 69 | AliPHOSRawDecoder(rawReader,mapping), |
70 | fSampleParamsLow(0x0),fSampleParamsHigh(0x0),fToFit(0x0) | |
77ea1c6f | 71 | { |
72 | //Construct a decoder object. | |
73 | //Is is user responsibility to provide next raw event | |
74 | //using AliRawReader::NextEvent(). | |
75 | ||
76 | if(!gMinuit) | |
77 | gMinuit = new TMinuit(100); | |
11f2ec15 | 78 | fSampleParamsHigh =new TArrayD(5) ; |
79 | fSampleParamsHigh->AddAt(4.25,0) ; | |
80 | fSampleParamsHigh->AddAt(0.094,1) ; | |
81 | fSampleParamsHigh->AddAt(0.0151,2) ; | |
82 | fSampleParamsHigh->AddAt(0.0384,3) ; | |
83 | fSampleParamsLow=new TArrayD(5) ; | |
84 | fSampleParamsLow->AddAt(5.14,0) ; | |
85 | fSampleParamsLow->AddAt(0.0970,1) ; | |
86 | fSampleParamsLow->AddAt(0.0088,2) ; | |
87 | fSampleParamsLow->AddAt(0.0346,3) ; | |
88 | fToFit = new TList() ; | |
77ea1c6f | 89 | } |
90 | ||
91 | //----------------------------------------------------------------------------- | |
92 | AliPHOSRawDecoderv1::~AliPHOSRawDecoderv1() | |
93 | { | |
94 | //Destructor. | |
11f2ec15 | 95 | if(fSampleParamsLow){ |
96 | delete fSampleParamsLow ; | |
97 | fSampleParamsLow=0 ; | |
98 | } | |
99 | if(fSampleParamsHigh){ | |
100 | delete fSampleParamsHigh ; | |
101 | fSampleParamsHigh=0; | |
102 | } | |
103 | if(fToFit){ | |
104 | delete fToFit ; | |
105 | fToFit=0 ; | |
106 | } | |
77ea1c6f | 107 | } |
108 | ||
109 | //----------------------------------------------------------------------------- | |
110 | AliPHOSRawDecoderv1::AliPHOSRawDecoderv1(const AliPHOSRawDecoderv1 &phosDecoder ): | |
11f2ec15 | 111 | AliPHOSRawDecoder(phosDecoder), |
112 | fSampleParamsLow(0x0),fSampleParamsHigh(0x0),fToFit(0x0) | |
77ea1c6f | 113 | { |
114 | //Copy constructor. | |
11f2ec15 | 115 | fToFit = new TList() ; |
116 | fSampleParamsLow =new TArrayD(*(phosDecoder.fSampleParamsLow)) ; | |
117 | fSampleParamsHigh=new TArrayD(*(phosDecoder.fSampleParamsHigh)) ; | |
77ea1c6f | 118 | } |
119 | ||
120 | //----------------------------------------------------------------------------- | |
11f2ec15 | 121 | AliPHOSRawDecoderv1& AliPHOSRawDecoderv1::operator = (const AliPHOSRawDecoderv1 &phosDecoder) |
77ea1c6f | 122 | { |
123 | //Assignment operator. | |
124 | ||
11f2ec15 | 125 | // if(this != &phosDecoder) { |
77ea1c6f | 126 | // } |
11f2ec15 | 127 | fToFit = new TList() ; |
128 | if(fSampleParamsLow){ | |
129 | fSampleParamsLow = phosDecoder.fSampleParamsLow ; | |
130 | fSampleParamsHigh= phosDecoder.fSampleParamsHigh ; | |
131 | } | |
132 | else{ | |
133 | fSampleParamsLow =new TArrayD(*(phosDecoder.fSampleParamsLow)) ; | |
134 | fSampleParamsHigh=new TArrayD(*(phosDecoder.fSampleParamsHigh)) ; | |
135 | } | |
77ea1c6f | 136 | return *this; |
137 | } | |
138 | ||
139 | //----------------------------------------------------------------------------- | |
140 | Bool_t AliPHOSRawDecoderv1::NextDigit() | |
141 | { | |
142 | //Extract an energy deposited in the crystal, | |
143 | //crystal' position (module,column,row), | |
144 | //time and gain (high or low). | |
145 | //First collects sample, then evaluates it and if it has | |
146 | //reasonable shape, fits it with Gamma2 function and extracts | |
147 | //energy and time. | |
11f2ec15 | 148 | |
149 | // TCanvas * c = (TCanvas *)gROOT->FindObjectAny("canvMy") ; | |
150 | // TH1S * h = new TH1S("s","",200,0.5,200.5) ; | |
151 | // TF1 * fff = new TF1("fff","[0]+[1]*((x-[2])+[3]*(x-[2])*(x-[2]))*(exp(-(x-[2])*[4])+[5]*exp(-(x-[2])*[6]))",0.,1000.) ; | |
77ea1c6f | 152 | |
153 | AliCaloRawStream* in = fCaloStream; | |
154 | ||
155 | Int_t iBin = 0; | |
156 | Int_t tLength = 0; | |
157 | fEnergy = -111; | |
158 | Float_t pedMean = 0; | |
159 | Int_t nPed = 0; | |
160 | Float_t baseLine = 1.0; | |
161 | const Float_t nPreSamples = 10; | |
162 | ||
163 | while ( in->Next() ) { | |
164 | ||
165 | if(!tLength) { | |
166 | tLength = in->GetTimeLength(); | |
167 | if(tLength!=fSamples->GetSize()) { | |
168 | delete fSamples ; | |
169 | fSamples = new TArrayI(tLength); | |
170 | } | |
171 | else{ | |
172 | for(Int_t i=0; i<fSamples->GetSize(); i++){ | |
173 | fSamples->AddAt(0,i) ; | |
174 | } | |
175 | } | |
176 | } | |
177 | ||
178 | // Fit the full sample | |
179 | if(in->IsNewHWAddress() && iBin>0) { | |
180 | ||
181 | Double_t pedestal =0. ; | |
182 | if(fPedSubtract){ | |
183 | if (nPed > 0) | |
184 | pedestal = (Double_t)(pedMean/nPed); | |
185 | else | |
186 | return kFALSE; | |
187 | } | |
188 | ||
189 | //calculate energy | |
190 | //first estimate if this sample looks like gamma2 function | |
191 | Double_t aMean=0. ; | |
192 | Double_t aRMS=0. ; | |
193 | Int_t tStart = 0 ; | |
194 | Int_t cnts=0 ; | |
195 | for(Int_t i=0; i<fSamples->GetSize(); i++){ | |
196 | if(fSamples->At(i)>0){ | |
197 | Double_t de=fSamples->At(i)-pedestal ; | |
198 | aMean+=de ; | |
199 | aRMS+=de*de ; | |
200 | cnts++; | |
201 | if(de>2 && tStart==0) | |
202 | tStart=i ; | |
203 | if(fSamples->At(i)>fEnergy) | |
204 | fEnergy=fSamples->At(i) ; | |
205 | } | |
206 | } | |
207 | if(cnts>0){ | |
208 | aMean/=cnts; | |
209 | aRMS=aRMS/cnts-aMean*aMean; | |
210 | } | |
211 | ||
77ea1c6f | 212 | //IF sample has reasonable mean and RMS, try to fit it with gamma2 |
213 | if(fEnergy>2.&& cnts >20 && aMean>0. && aRMS>2.){ //more or less reasonable sample | |
214 | ||
215 | gMinuit->mncler(); // Reset Minuit's list of paramters | |
216 | gMinuit->SetPrintLevel(-1) ; // No Printout | |
217 | gMinuit->SetFCN(AliPHOSRawDecoderv1::UnfoldingChiSquare) ; | |
218 | // To set the address of the minimization function | |
11f2ec15 | 219 | |
220 | fToFit->Clear() ; | |
221 | if(fLowGainFlag){ | |
222 | fSampleParamsLow->AddAt(pedestal,4) ; | |
223 | fToFit->AddFirst((TObject*)fSampleParamsLow) ; | |
224 | } | |
225 | else{ | |
226 | fSampleParamsHigh->AddAt(pedestal,4) ; | |
227 | fToFit->AddFirst((TObject*)fSampleParamsHigh) ; | |
228 | } | |
229 | fToFit->AddLast((TObject*)fSamples) ; | |
230 | ||
231 | gMinuit->SetObjectFit((TObject*)fToFit) ; // To tranfer pointer to UnfoldingChiSquare | |
77ea1c6f | 232 | Int_t ierflg ; |
11f2ec15 | 233 | gMinuit->mnparm(0, "t0", 1.*tStart, 0.1, 0, 0, ierflg) ; |
77ea1c6f | 234 | if(ierflg != 0){ |
11f2ec15 | 235 | // AliWarning(Form("Unable to set initial value for fit procedure : t0=%e\n",1.*tStart) ) ; |
77ea1c6f | 236 | fEnergy=0. ; |
237 | fTime=-999. ; | |
238 | return kTRUE ; //will scan further | |
239 | } | |
11f2ec15 | 240 | Double_t amp0=(fEnergy-pedestal)*0.0032; |
241 | ||
242 | gMinuit->mnparm(1, "Energy", amp0 , 0.001*amp0, 0, 0, ierflg) ; | |
77ea1c6f | 243 | if(ierflg != 0){ |
11f2ec15 | 244 | // AliWarning(Form("Unable to set initial value for fit procedure : E=%e\n", amp0)) ; |
77ea1c6f | 245 | fEnergy=0. ; |
246 | fTime=-999. ; | |
247 | return kTRUE ; //will scan further | |
248 | } | |
249 | ||
11f2ec15 | 250 | Double_t p0 = 0.0001 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly |
77ea1c6f | 251 | // depends on it. |
252 | Double_t p1 = 1.0 ; | |
253 | Double_t p2 = 0.0 ; | |
254 | gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls | |
255 | gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient | |
256 | // gMinuit->SetMaxIterations(100); | |
257 | gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings | |
258 | ||
259 | gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize | |
260 | ||
11f2ec15 | 261 | Double_t err,t0err ; |
262 | Double_t t0,efit ; | |
263 | gMinuit->GetParameter(0,t0, t0err) ; | |
264 | gMinuit->GetParameter(1,efit, err) ; | |
77ea1c6f | 265 | |
11f2ec15 | 266 | Double_t a,alpha ; |
267 | if(fLowGainFlag){ | |
268 | a=fSampleParamsLow->At(0) ; | |
269 | alpha=fSampleParamsLow->At(1) ; | |
270 | } | |
271 | else{ | |
272 | a=fSampleParamsHigh->At(0) ; | |
273 | alpha=fSampleParamsHigh->At(1) ; | |
274 | } | |
275 | ||
276 | // c->cd() ; | |
277 | // h->Draw() ; | |
278 | // if(fLowGainFlag){ | |
279 | // fff->SetParameters(pedestal,efit,t0,a,alpha,fSampleParamsLow->At(2),fSampleParamsLow->At(3)) ; | |
280 | // } | |
281 | // else{ | |
282 | // fff->SetParameters(pedestal,efit,t0,a,alpha,fSampleParamsHigh->At(2),fSampleParamsHigh->At(3)) ; | |
283 | // } | |
284 | // fff->Draw("same") ; | |
285 | // c->Update(); | |
286 | ||
287 | efit*=(2.*a+TMath::Sqrt(4.*a*a+alpha*alpha))/alpha/alpha*TMath::Exp(-1.+(alpha-TMath::Sqrt(4.*a*a+alpha*alpha))/2./a) ; | |
288 | //printf("efit=%f, t0=%e +- %e, ped=%f \n",efit,t0,t0err,pedestal) ; | |
77ea1c6f | 289 | Double_t fmin,fedm,errdef ; |
290 | Int_t npari,nparx,istat; | |
291 | gMinuit->mnstat(fmin,fedm,errdef,npari,nparx,istat) ; | |
11f2ec15 | 292 | |
293 | //if(fLowGainFlag) | |
294 | // printf("LowGain \n") ; | |
295 | //else | |
296 | // printf("highGain \n") ; | |
297 | ||
298 | //printf("fmin=%e \n",fmin) ; | |
299 | //getchar() ; | |
300 | ||
301 | if(1){ //fmin < 3.+0.3*efit ){ //Chi^2 of a good sample | |
302 | if(efit>0.){ | |
77ea1c6f | 303 | fEnergy=efit ; |
304 | fTime=t0 ; | |
305 | } | |
306 | } | |
307 | else{ | |
308 | fEnergy=0 ; //bad sample | |
309 | fTime=-999.; | |
310 | } | |
77ea1c6f | 311 | |
11f2ec15 | 312 | fTime*=fPulseGenerator->GetRawFormatTimeTrigger() ; |
313 | ||
77ea1c6f | 314 | if (fEnergy < baseLine) fEnergy = 0; |
315 | } | |
316 | else{ //bad sample | |
317 | fEnergy=0. ; | |
318 | fTime=-999. ; | |
319 | } | |
320 | ||
321 | return kTRUE; | |
322 | } | |
323 | ||
324 | fLowGainFlag = in->IsLowGain(); | |
325 | fTime = fPulseGenerator->GetRawFormatTimeTrigger() * in->GetTime(); | |
326 | fModule = in->GetModule()+1; | |
327 | fRow = in->GetRow() +1; | |
328 | fColumn = in->GetColumn()+1; | |
329 | ||
330 | ||
331 | // Fill array with samples | |
332 | iBin++; | |
333 | if(tLength-iBin < nPreSamples) { | |
334 | pedMean += in->GetSignal(); | |
335 | nPed++; | |
336 | } | |
337 | fSamples->AddAt(in->GetSignal(),tLength-iBin); | |
11f2ec15 | 338 | // h->SetBinContent(tLength-iBin+1,in->GetSignal()) ; |
77ea1c6f | 339 | |
340 | } // in.Next() | |
341 | ||
342 | return kFALSE; | |
343 | } | |
344 | //_____________________________________________________________________________ | |
345 | void AliPHOSRawDecoderv1::UnfoldingChiSquare(Int_t & /*nPar*/, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag) | |
346 | { | |
347 | // Calculates the Chi square for the samples minimization | |
348 | // Number of parameters, Gradient, Chi squared, parameters, what to do | |
349 | ||
11f2ec15 | 350 | TList * toFit= (TList*)gMinuit->GetObjectFit() ; |
351 | TArrayD * params=(TArrayD*)toFit->At(0) ; | |
352 | TArrayI * samples = (TArrayI*)toFit->At(1) ; | |
77ea1c6f | 353 | |
354 | fret = 0. ; | |
355 | if(iflag == 2) | |
11f2ec15 | 356 | for(Int_t iparam = 0 ; iparam < 2 ; iparam++) |
77ea1c6f | 357 | Grad[iparam] = 0 ; // Will evaluate gradient |
358 | ||
11f2ec15 | 359 | Int_t nSamples=samples->GetSize() ; //Math::Min(70,samples->GetSize()) ; |
360 | Double_t t0=x[0] ; | |
361 | Double_t en=x[1] ; | |
362 | Double_t a=params->At(0) ; | |
363 | Double_t alpha=params->At(1) ; | |
364 | Double_t b=params->At(2) ; | |
365 | Double_t beta=params->At(3) ; | |
77ea1c6f | 366 | |
367 | for(Int_t i = 0 ; i < nSamples ; i++) { | |
11f2ec15 | 368 | if(samples->At(i)==0 || samples->At(i)==1023) //zero or overflow |
77ea1c6f | 369 | continue ; |
370 | Double_t dt=i*1.-t0 ; | |
11f2ec15 | 371 | Double_t diff=float(samples->At(i))-Gamma2(dt,en,params) ; |
372 | Double_t w=0.1+0.005*i ; //Mean Pedestal RMS + rising modulation | |
77ea1c6f | 373 | // if(w==0)w=1. ; |
374 | diff/=w ; | |
375 | if(iflag == 2){ // calculate gradient | |
77ea1c6f | 376 | if(dt>=0.){ |
11f2ec15 | 377 | Grad[0] += -2.*en*diff*((alpha*dt*(1.+a*dt)-1.-2.*a*dt)*TMath::Exp(-alpha*dt)+ |
378 | b*(beta*dt*(1.+a*dt)-1.-2.*a*dt)*TMath::Exp(-beta*dt)) /w ; //derivative over t0 | |
379 | Grad[1] += -2.*diff*(dt+a*dt*dt)*(TMath::Exp(-alpha*dt)+ | |
380 | b*TMath::Exp(-dt*beta))/w ; | |
77ea1c6f | 381 | } |
382 | } | |
383 | fret += diff*diff ; | |
384 | } | |
77ea1c6f | 385 | if(nSamples){ |
386 | fret/=nSamples ; | |
387 | if(iflag == 2){ | |
11f2ec15 | 388 | for(Int_t iparam = 0 ; iparam < 2 ; iparam++) |
77ea1c6f | 389 | Grad[iparam] /= nSamples ; |
390 | } | |
391 | } | |
77ea1c6f | 392 | |
393 | } | |
394 | //----------------------------------------------------------------------------- | |
11f2ec15 | 395 | Double_t AliPHOSRawDecoderv1::Gamma2(Double_t dt,Double_t en,TArrayD * params){ //Function for fitting samples |
77ea1c6f | 396 | //parameters: |
11f2ec15 | 397 | //dt-time after start |
77ea1c6f | 398 | //en-amplutude |
11f2ec15 | 399 | //function parameters |
400 | ||
401 | Double_t ped=params->At(4) ; | |
77ea1c6f | 402 | if(dt<0.) |
11f2ec15 | 403 | return ped ; //pedestal |
77ea1c6f | 404 | else |
11f2ec15 | 405 | return ped+en*(dt+params->At(0)*dt*dt)*(TMath::Exp(-dt*params->At(1))+params->At(2)*TMath::Exp(-dt*params->At(3))) ; |
77ea1c6f | 406 | } |
407 |