]>
Commit | Line | Data |
---|---|---|
fb3d369e | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, 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 | // | |
17 | // class to calculate TRD dEdx | |
18 | // xx | |
19 | // xx | |
20 | // xx | |
21 | // xx | |
22 | // | |
23 | // Xianguo Lu | |
24 | // lu@physi.uni-heidelberg.de | |
25 | // Xianguo.Lu@cern.ch | |
26 | // | |
27 | // | |
28 | ||
29 | ||
30 | #include "TF1.h" | |
31 | #include "TFile.h" | |
32 | #include "TH1D.h" | |
33 | #include "TH2D.h" | |
34 | #include "THnSparse.h" | |
35 | #include "TMath.h" | |
36 | #include "TMatrixD.h" | |
37 | #include "TMinuit.h" | |
38 | #include "TObjArray.h" | |
39 | #include "TRandom3.h" | |
40 | #include "TStopwatch.h" | |
41 | #include "TVectorD.h" | |
42 | ||
43 | #include "TTreeStream.h" | |
44 | ||
45 | #include "AliCDBId.h" | |
46 | #include "AliCDBMetaData.h" | |
47 | #include "AliCDBStorage.h" | |
48 | #include "AliESDEvent.h" | |
49 | #include "AliESDfriendTrack.h" | |
50 | #include "AliESDtrack.h" | |
d0037ea4 | 51 | #include "AliTRDcalibDB.h" |
52 | #include "AliTRDCalROC.h" | |
fb3d369e | 53 | #include "AliTRDtrackV1.h" |
54 | ||
d0037ea4 | 55 | |
fb3d369e | 56 | #include "AliTRDdEdxUtils.h" |
57 | ||
58 | #define EPSILON 1e-12 | |
59 | ||
60 | THnSparseD * AliTRDdEdxUtils::fgHistGain=0x0; | |
61 | THnSparseD * AliTRDdEdxUtils::fgHistT0=0x0; | |
62 | THnSparseD * AliTRDdEdxUtils::fgHistVd=0x0; | |
63 | TObjArray * AliTRDdEdxUtils::fgHistPHQ=new TObjArray(8); | |
64 | ||
65 | TString AliTRDdEdxUtils::fgCalibFile; | |
66 | TObjArray * AliTRDdEdxUtils::fgObjGain = 0x0; | |
67 | TObjArray * AliTRDdEdxUtils::fgObjT0 = 0x0; | |
68 | TObjArray * AliTRDdEdxUtils::fgObjVd = 0x0; | |
69 | TObjArray * AliTRDdEdxUtils::fgObjPHQ = 0x0; | |
70 | ||
71 | Int_t AliTRDdEdxUtils::fgNchamber = -999; | |
72 | Double_t AliTRDdEdxUtils::fgChamberQ[6]; | |
73 | Double_t AliTRDdEdxUtils::fgChamberTmean[6]; | |
74 | ||
75 | Double_t AliTRDdEdxUtils::fgTrackTmean = -999; | |
76 | ||
d0037ea4 | 77 | Bool_t AliTRDdEdxUtils::fgPadGainOn = kTRUE; |
9a8b0e85 | 78 | Bool_t AliTRDdEdxUtils::fgExBOn = kTRUE; |
d0037ea4 | 79 | Double_t AliTRDdEdxUtils::fgQScale = 45; |
9a8b0e85 | 80 | Double_t AliTRDdEdxUtils::fgQ0Frac = 0.3; |
81 | Double_t AliTRDdEdxUtils::fgQ1Frac = 0.5; | |
82 | Double_t AliTRDdEdxUtils::fgTimeBinCountCut = 0.0; | |
83 | Int_t AliTRDdEdxUtils::fgCalibTPCnclsCut = 70; | |
84 | ||
fb3d369e | 85 | //=================================================================================== |
86 | // Math and Histogram | |
87 | //=================================================================================== | |
88 | void AliTRDdEdxUtils::GetCDFCuts(const TH1D *hh, const Int_t ncut, Double_t cuts[], const Double_t cdfs[], const Double_t thres) | |
89 | { | |
90 | // | |
91 | //counts of hh is sorted | |
92 | // | |
93 | ||
94 | for(Int_t ii=0; ii<ncut; ii++){ | |
95 | cuts[ii] = -999; | |
96 | } | |
97 | ||
98 | Int_t nsel = 0; | |
99 | const Int_t nbin = hh->GetNbinsX(); | |
100 | Double_t datas[nbin]; | |
101 | for(Int_t ii=1; ii<=nbin; ii++){ | |
102 | const Double_t res = hh->GetBinContent(ii); | |
103 | if(res<thres){ | |
104 | continue; | |
105 | } | |
106 | ||
107 | datas[nsel] = res; | |
108 | nsel++; | |
109 | } | |
110 | if(!nsel) | |
111 | return; | |
112 | ||
113 | Int_t id[nsel]; | |
114 | TMath::Sort(nsel, datas, id, kFALSE); | |
115 | ||
116 | for(Int_t ii=0; ii<ncut; ii++){ | |
117 | const Double_t icdf = cdfs[ii]; | |
118 | if(icdf<0 || icdf>1){ | |
119 | printf("AliTRDdEdxUtils::GetCDFCuts error cdfs[%d] %15f out of range!\n", ii, icdf); exit(1); | |
120 | } | |
121 | cuts[ii] = datas[id[Int_t(icdf*nsel)]]; | |
122 | } | |
123 | } | |
124 | ||
125 | Double_t AliTRDdEdxUtils::GetMeanRMS(const Double_t nn, const Double_t sum, const Double_t w2s, Double_t * grms, Double_t * gerr) | |
126 | { | |
127 | // | |
128 | //calculate mean (with error) and rms from sum, w2s, nn | |
129 | //if nn=0, mean, error, and rms all = 0 | |
130 | // | |
131 | ||
132 | Double_t tmean = 0, trms = 0, terr = 0; | |
133 | ||
134 | if(nn>EPSILON){ | |
135 | tmean = sum/nn; | |
136 | ||
137 | const Double_t arg = w2s/nn-tmean*tmean; | |
138 | if(TMath::Abs(arg)<EPSILON){ | |
139 | trms = 0; | |
140 | } | |
141 | else{ | |
142 | if( arg <0 ){ | |
143 | printf("AliTRDdEdxUtils::GetMeanRMS error negative sqrt argument!! %e -- %e %e %f\n", arg, w2s, sum, nn); exit(1); | |
144 | } | |
145 | ||
146 | trms = TMath::Sqrt(arg); | |
147 | } | |
148 | ||
149 | terr = trms/TMath::Sqrt(nn); | |
150 | } | |
151 | ||
152 | if(grms){ | |
153 | (*grms) = trms; | |
154 | } | |
155 | ||
156 | if(gerr){ | |
157 | (*gerr) = terr; | |
158 | } | |
159 | ||
160 | return tmean; | |
161 | } | |
162 | ||
163 | Double_t AliTRDdEdxUtils::TruncatedMean(const Int_t nx, const Double_t xdata[], const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr, Double_t *wws) | |
164 | { | |
165 | // | |
166 | //calculate truncated mean | |
167 | //return <x*w>_{low-high according to x} | |
168 | // | |
169 | ||
170 | /* | |
171 | //test-> | |
172 | for(Int_t ii=0; ii<nx; ii++){ | |
173 | printf("test %d/%d %f\n", ii, nx, xdata[ii]); | |
174 | } | |
175 | //<--test | |
176 | */ | |
177 | ||
178 | Int_t index[nx]; | |
179 | TMath::Sort(nx, xdata, index, kFALSE); | |
180 | ||
181 | Int_t nused = 0; | |
182 | Double_t sum = 0; | |
183 | Double_t w2s = 0; | |
184 | const Int_t istart = Int_t (nx*lowfrac); | |
185 | const Int_t istop = Int_t (nx*highfrac); | |
186 | ||
187 | //=,< correct, because when low=0, high=1 it is correct | |
188 | for(Int_t ii=istart; ii<istop; ii++){ | |
189 | Double_t weight = 1; | |
190 | if(wws){ | |
191 | weight = wws[index[ii]]; | |
192 | } | |
193 | const Double_t sx = xdata[index[ii]]*weight; | |
194 | ||
195 | sum += sx; | |
196 | w2s += sx*sx; | |
197 | ||
198 | nused++; | |
199 | //printf("test in loop %d/%d %f %f %f\n", ii, nused, sx, sum, w2s); | |
200 | ||
201 | } | |
202 | ||
203 | return GetMeanRMS(nused, sum, w2s, grms, gerr); | |
204 | } | |
205 | ||
206 | Double_t AliTRDdEdxUtils::TruncatedMean(const TH1 *hh, const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr) | |
207 | { | |
208 | // | |
209 | //do truncation on histogram | |
210 | // | |
211 | //if hh is scaled, be sure Sumw2 is called before scaling!! then mean, rms and err will all be correct | |
212 | ||
213 | //with under-/over-flow | |
214 | Double_t npreTrunc = 0; | |
215 | for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){ | |
216 | const Double_t be = hh->GetBinError(itmp); | |
217 | const Double_t bc = hh->GetBinContent(itmp); | |
218 | if(be<EPSILON){ | |
219 | if(bc>EPSILON){ | |
220 | printf("AliTRDdEdxUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1); | |
221 | } | |
222 | continue; | |
223 | } | |
224 | npreTrunc += bc*bc/be/be; | |
225 | } | |
226 | ||
227 | const Double_t nstart = npreTrunc*lowfrac; | |
228 | const Double_t nstop = npreTrunc*highfrac; | |
229 | ||
230 | //with Double_t this should also handle normalized hist | |
231 | Double_t ntot = 0; | |
232 | Double_t nused = 0; | |
233 | Double_t sum = 0; | |
234 | Double_t w2s = 0; | |
235 | for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){ | |
236 | const Double_t be = hh->GetBinError(itmp); | |
237 | const Double_t bc = hh->GetBinContent(itmp); | |
238 | if(be<EPSILON){ | |
239 | if(bc>EPSILON){ | |
240 | printf("AliTRDdEdxUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1); | |
241 | } | |
242 | continue; | |
243 | } | |
244 | const Double_t weight = bc*bc/be/be; | |
245 | ntot+=weight; | |
246 | //<= correct, because when high=1, nstop has to be included | |
247 | if(ntot>nstart && ntot<=nstop){ | |
248 | ||
249 | const Double_t val = hh->GetBinCenter(itmp); | |
250 | sum += weight*val; | |
251 | w2s += weight*val*val; | |
252 | ||
253 | nused += weight; | |
254 | ||
255 | //printf("test %d %f %f --- %f %f -- %f %f\n", itmp, weight, val, sum, w2s, nused, nsample); | |
256 | } | |
257 | else if(ntot>nstop){ | |
258 | if(itmp>=hh->GetNbinsX()){ | |
259 | printf("AliTRDdEdxUtils::TruncatedMean warning hist range too small %s %f %f %d %d, %15f %15f %15f; nused w2s sum set to 0\n", hh->GetName(), hh->GetBinLowEdge(1), hh->GetBinLowEdge(itmp), itmp, hh->GetNbinsX(), hh->GetBinContent(hh->GetNbinsX())/hh->Integral(0,hh->GetNbinsX()+1), hh->GetBinContent(hh->GetNbinsX()), hh->Integral(0,hh->GetNbinsX()+1)); //exit(1); | |
260 | nused = 0; | |
261 | w2s = sum = 0; | |
262 | } | |
263 | break; | |
264 | } | |
265 | } | |
266 | ||
267 | return GetMeanRMS(nused, sum, w2s, grms, gerr); | |
268 | } | |
269 | ||
270 | void AliTRDdEdxUtils::FitSlicesY(const TH2D *hh, TH1D *&hnor, TH1D *&hmpv, TH1D *&hwid, TH1D *&hres, const Double_t thres, const Double_t lowfrac, const Double_t highfrac) | |
271 | { | |
272 | // | |
273 | //fit slices of hh using truncation | |
274 | // | |
275 | ||
276 | const Int_t x0 = hh->GetXaxis()->GetFirst(); | |
277 | const Int_t x1 = hh->GetXaxis()->GetLast(); | |
278 | const Int_t y0 = hh->GetYaxis()->GetFirst(); | |
279 | const Int_t y1 = hh->GetYaxis()->GetLast(); | |
280 | ||
281 | const Int_t nx = hh->GetNbinsX(); | |
282 | const Int_t ny = hh->GetNbinsY(); | |
283 | const Double_t xmin = hh->GetXaxis()->GetXmin(); | |
284 | const Double_t xmax = hh->GetXaxis()->GetXmax(); | |
285 | const Double_t ymin = hh->GetYaxis()->GetXmin(); | |
286 | const Double_t ymax = hh->GetYaxis()->GetXmax(); | |
287 | ||
288 | hnor = new TH1D(Form("%s_amp",hh->GetName()), "", nx, xmin, xmax); | |
289 | hmpv = new TH1D(Form("%s_mpv",hh->GetName()), "", nx, xmin, xmax); | |
290 | hwid = new TH1D(Form("%s_wid",hh->GetName()), "", nx, xmin, xmax); | |
291 | hres = new TH1D(Form("%s_res",hh->GetName()), "", nx, xmin, xmax); | |
292 | ||
293 | for(Int_t ix=x0; ix<=x1; ix++){ | |
294 | //to speed up | |
295 | const Double_t rawcount = hh->Integral(ix,ix,0, ny+1); | |
296 | if(rawcount<EPSILON){ | |
297 | continue; | |
298 | } | |
299 | ||
300 | TH1D *htmp = new TH1D(Form("FitSlicesY_%s_%d", hh->GetName(), ix),"",ny, ymin, ymax); | |
301 | Double_t ntot = 0; | |
302 | for(Int_t iy=y0; iy<=y1; iy++){ | |
303 | const Double_t be = hh->GetBinError(ix,iy); | |
304 | const Double_t bc = hh->GetBinContent(ix, iy); | |
305 | ||
306 | if(be<EPSILON){ | |
307 | if(bc>EPSILON){ | |
308 | printf("AliTRDdEdxUtils::FitSlicesY error %d %d %e %e\n", ix, iy, be, bc); exit(1); | |
309 | } | |
310 | continue; | |
311 | } | |
312 | ||
313 | htmp->SetBinContent(iy, bc); | |
314 | htmp->SetBinError(iy, be); | |
315 | ||
316 | ntot += (bc/be)*(bc/be); | |
317 | ||
318 | //if(be) printf("test %d %d : %f %f %f\n", ix, iy, bc, be, pow(bc/be,2)); | |
319 | } | |
320 | ||
321 | hnor->SetBinContent(ix, ntot); | |
322 | hnor->SetBinError( ix, 0); | |
323 | ||
324 | if(ntot<thres || htmp->GetRMS()<EPSILON){ | |
325 | delete htmp; | |
326 | continue; | |
327 | } | |
328 | ||
329 | //test htmp->Draw(); | |
330 | Double_t trms = -999, terr = -999; | |
331 | const Double_t tmean = TruncatedMean(htmp, lowfrac, highfrac, &trms, &terr); | |
332 | ||
333 | hmpv->SetBinContent(ix, tmean); | |
334 | hmpv->SetBinError( ix, terr); | |
335 | ||
336 | hwid->SetBinContent(ix, trms); | |
337 | hwid->SetBinError( ix, 0); | |
338 | ||
339 | hres->SetBinContent(ix, tmean>EPSILON ? trms/tmean:0); | |
340 | hres->SetBinError( ix, 0); | |
341 | ||
342 | delete htmp; | |
343 | } | |
344 | ||
345 | TH1 *hhs[]={hnor, hmpv, hwid, hres}; | |
346 | const TString yt[]={"N", "MPV", "#sigma", "#sigma/MPV"}; | |
347 | const Int_t nh = sizeof(hhs)/sizeof(TH1*); | |
348 | for(Int_t ii=0; ii<nh; ii++){ | |
349 | hhs[ii]->SetYTitle(Form("%s of %s", yt[ii].Data(), hh->GetYaxis()->GetTitle())); | |
350 | hhs[ii]->SetXTitle(hh->GetXaxis()->GetTitle()); | |
351 | hhs[ii]->GetYaxis()->SetTitleOffset(hh->GetYaxis()->GetTitleOffset()); | |
352 | hhs[ii]->SetTitle(hh->GetTitle()); | |
353 | } | |
354 | } | |
355 | ||
356 | //=================================================================================== | |
357 | // TRD Analysis Fast Tool | |
358 | //=================================================================================== | |
359 | ||
360 | Int_t AliTRDdEdxUtils::GetNtracklet(const AliESDEvent *esd) | |
361 | { | |
362 | // | |
363 | //number of trd tracklet in one esd event | |
364 | // | |
365 | const Int_t ntrk0 = esd->GetNumberOfTracks(); | |
366 | Int_t ntrdv1=0; | |
367 | for(Int_t ii=0; ii<ntrk0; ii++){ | |
368 | ntrdv1 += esd->GetTrack(ii)->GetTRDntracklets(); | |
369 | } | |
370 | return ntrdv1; | |
371 | } | |
372 | ||
373 | AliTRDtrackV1 * AliTRDdEdxUtils::GetTRDtrackV1(const AliESDtrack * esdtrack) | |
374 | { | |
375 | // | |
376 | //Get TRD friend track | |
377 | // | |
378 | ||
379 | AliESDfriendTrack * friendtrk = (AliESDfriendTrack *)esdtrack->GetFriendTrack(); | |
380 | if(!friendtrk){ | |
381 | //printf("xlulog AliAnalysisTaskCosmicTRD::GetTRDtrack no friend!!\n"); exit(1); | |
382 | return 0x0; | |
383 | } | |
384 | ||
385 | TObject *calibObject=0x0; | |
386 | AliTRDtrackV1 * trdtrack=0x0; | |
387 | for(Int_t l=0; (calibObject=friendtrk->GetCalibObject(l)); l++) { | |
388 | if( (trdtrack=dynamic_cast<AliTRDtrackV1*>(calibObject)) ) | |
389 | break; | |
390 | } | |
391 | ||
392 | return trdtrack; | |
393 | } | |
394 | ||
395 | Bool_t AliTRDdEdxUtils::IsInSameStack(const AliTRDtrackV1 *trdtrack) | |
396 | { | |
397 | // | |
398 | // to check if all tracklets are in the same stack, useful in cosmic | |
399 | // | |
400 | ||
401 | TVectorD secs(18), stks(5); | |
402 | ||
403 | for(Int_t ilayer = 0; ilayer < 6; ilayer++){ | |
404 | AliTRDseedV1 *tracklet=trdtrack->GetTracklet(ilayer); | |
405 | if(!tracklet) | |
406 | continue; | |
407 | ||
408 | const Int_t det = tracklet->GetDetector(); | |
409 | const Int_t isector = AliTRDgeometry::GetSector(det); | |
410 | const Int_t istack = AliTRDgeometry::GetStack(det); | |
411 | ||
412 | secs[isector] = 1; | |
413 | stks[istack] = 1; | |
414 | } | |
415 | ||
416 | if(secs.Sum()!=1 || stks.Sum()!=1){ | |
417 | return kFALSE; | |
418 | } | |
419 | else | |
420 | return kTRUE; | |
421 | } | |
422 | ||
423 | Bool_t AliTRDdEdxUtils::GetFirstSectorStackMomentum(const AliTRDtrackV1 *trdtrack, Int_t & isec, Int_t & istk, Double_t & mom) | |
424 | { | |
425 | // | |
426 | //as function name | |
427 | // | |
428 | isec = istk = -999; | |
429 | mom = -999; | |
430 | ||
431 | for(Int_t ilayer = 0; ilayer < 6; ilayer++){ | |
432 | AliTRDseedV1 *tracklet=trdtrack->GetTracklet(ilayer); | |
433 | if(!tracklet) | |
434 | continue; | |
435 | ||
436 | const Int_t det = tracklet->GetDetector(); | |
437 | isec = AliTRDgeometry::GetSector(det); | |
438 | istk = AliTRDgeometry::GetStack(det); | |
439 | ||
440 | mom = tracklet->GetMomentum(); | |
441 | ||
442 | break; | |
443 | } | |
444 | ||
445 | if(isec<0) | |
446 | return kFALSE; | |
447 | else | |
448 | return kTRUE; | |
449 | } | |
450 | ||
451 | //=================================================================================== | |
452 | // Calibration | |
453 | //=================================================================================== | |
454 | Double_t AliTRDdEdxUtils::GetCalibTPCscale(const Int_t tpcncls, const Double_t tpcsig) | |
455 | { | |
456 | // | |
457 | //the scale used in calibration | |
458 | // | |
459 | ||
460 | if(tpcncls < CalibTPCnclsCut()) | |
461 | return -999; | |
462 | ||
463 | if(tpcsig<EPSILON) | |
464 | return -999; | |
465 | ||
466 | return tpcsig/120; | |
467 | ||
468 | } | |
469 | ||
470 | Int_t AliTRDdEdxUtils::GetPHQIterator(const Bool_t kinvq, const Double_t mag, const Int_t charge) | |
471 | { | |
472 | // | |
473 | //iterator for calib obj and hist | |
474 | // | |
475 | return kinvq*4 + (mag>0)*2 + (charge>0); | |
476 | } | |
477 | ||
26733ed4 | 478 | TObjArray * AliTRDdEdxUtils::GetObjPHQ() |
479 | { | |
480 | // | |
481 | //return fgObjPHQ, initialized if null | |
482 | // | |
483 | ||
484 | if(!fgObjPHQ){ | |
485 | fgObjPHQ = new TObjArray(8); | |
486 | } | |
487 | ||
488 | return fgObjPHQ; | |
489 | } | |
490 | ||
fb3d369e | 491 | TObjArray * AliTRDdEdxUtils::GetObjPHQ(const Bool_t kinvq, const Double_t mag, const Int_t charge) |
492 | { | |
493 | // | |
494 | //return calib obj | |
495 | // | |
496 | if(!fgObjPHQ){ | |
26733ed4 | 497 | printf("AliTRDdEdxUtils::GetObjPHQ(kinvq, mag, charge) error fgObjPHQ null!!\n"); exit(1); |
fb3d369e | 498 | } |
499 | ||
500 | return (TObjArray*) fgObjPHQ->At(GetPHQIterator(kinvq, mag, charge)); | |
501 | } | |
502 | ||
503 | THnSparseD * AliTRDdEdxUtils::GetHistPHQ(const Bool_t kinvq, const Double_t mag, const Int_t charge) | |
504 | { | |
505 | // | |
506 | //return calib hist | |
507 | // | |
508 | return (THnSparseD*) fgHistPHQ->At(GetPHQIterator(kinvq, mag, charge)); | |
509 | } | |
510 | ||
511 | TString AliTRDdEdxUtils::GetPHQName(const Bool_t kobj, const Int_t iter) | |
512 | { | |
513 | // | |
514 | //get name of calib obj/hist of PHQ | |
515 | // | |
516 | return Form("TRDCalib%sPHQ%d", kobj?"Obj":"Hist", iter); | |
517 | } | |
518 | ||
519 | void AliTRDdEdxUtils::DeleteCalibObj() | |
520 | { | |
521 | // | |
522 | //delete calib obj | |
523 | // | |
524 | delete fgObjGain; | |
525 | delete fgObjT0; | |
526 | delete fgObjVd; | |
527 | ||
528 | fgObjGain = 0x0; | |
529 | fgObjT0 = 0x0; | |
530 | fgObjVd = 0x0; | |
531 | ||
532 | if(fgObjPHQ){ | |
533 | fgObjPHQ->SetOwner(); | |
534 | delete fgObjPHQ; | |
535 | fgObjPHQ = 0x0; | |
536 | } | |
537 | } | |
538 | ||
539 | Bool_t AliTRDdEdxUtils::GenerateDefaultPHQOCDB(const TString path) | |
540 | { | |
541 | // | |
542 | //generate default OCDB object PHQ, do like | |
543 | //AliTRDdEdxUtils::GenerateDefaultPHQOCDB("local://./") | |
544 | // | |
545 | ||
546 | TObjArray * arr8 = new TObjArray(8); | |
547 | arr8->SetOwner(); | |
548 | ||
549 | for(Int_t ii=0; ii<8; ii++){ | |
550 | TObjArray * arr1 = new TObjArray(1); | |
551 | arr1->SetOwner(); | |
552 | arr1->SetName(GetPHQName(1, ii)); | |
553 | ||
554 | const Int_t nbins = NTRDtimebin(); | |
555 | TVectorD * vec = new TVectorD(nbins); | |
556 | for(Int_t jj=0; jj<nbins; jj++){ | |
557 | (*vec)[jj] = 1; | |
558 | } | |
559 | arr1->Add(vec); | |
560 | arr8->Add(arr1); | |
561 | } | |
562 | ||
563 | AliCDBMetaData *metaData= new AliCDBMetaData(); | |
564 | metaData->SetObjectClassName("TObjArray"); | |
565 | metaData->SetResponsible("Raphaelle Bailhache and Xianguo Lu"); | |
566 | ||
567 | AliCDBId id1("TRD/Calib/PHQ", 0, 999999999); | |
568 | AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage(path); | |
569 | gStorage->Put(arr8, id1, metaData); | |
570 | ||
571 | delete metaData; | |
572 | delete arr8; | |
573 | ||
574 | return kTRUE; | |
575 | } | |
576 | ||
577 | void AliTRDdEdxUtils::IniCalibObj() | |
578 | { | |
579 | // | |
580 | //set CalibObj from file, clone to static calib obj | |
581 | // | |
582 | ||
583 | DeleteCalibObj(); | |
584 | ||
585 | TFile *cfile=new TFile(fgCalibFile); | |
586 | if(!cfile->IsOpen()){ | |
587 | printf("AliTRDdEdxUtils::IniCalibObj error fgCalibFile not open! %s\n", fgCalibFile.Data());exit(1); | |
588 | } | |
589 | ||
590 | printf("\nAliTRDdEdxUtils::IniCalibObj file: %s\n", fgCalibFile.Data()); | |
591 | ||
592 | //--- | |
593 | const TString objnames[] ={"TRDCalibObjGain", "TRDCalibObjT0", "TRDCalibObjVd"}; | |
594 | TObjArray ** gobjs[]={ &fgObjGain, &fgObjT0, &fgObjVd}; | |
595 | ||
596 | const Int_t nobj = sizeof(objnames)/sizeof(TString); | |
597 | for(Int_t iobj=0; iobj<nobj; iobj++){ | |
598 | TObjArray *tmpo=0x0; | |
599 | cfile->GetObject(objnames[iobj], tmpo); | |
600 | if(!tmpo){ | |
601 | printf("AliTRDdEdxUtils::IniCalibObj error obj %s not found!\n", objnames[iobj].Data()); exit(1); | |
602 | } | |
603 | ||
604 | (*gobjs[iobj])=(TObjArray*)tmpo->Clone(); | |
605 | (*gobjs[iobj])->SetOwner(); | |
606 | } | |
607 | ||
608 | fgObjPHQ = new TObjArray(8); | |
609 | for(Int_t iter=0; iter<8; iter++){ | |
610 | const TString objn = GetPHQName(1, iter); | |
611 | TObjArray *tmpo=0x0; | |
612 | cfile->GetObject(objn, tmpo); | |
613 | if(!tmpo){ | |
614 | printf("AliTRDdEdxUtils::IniCalibObj error obj %s not found!\n", objn.Data()); exit(1); | |
615 | } | |
616 | ||
617 | TObjArray *obji=(TObjArray*) tmpo->Clone(); | |
618 | obji->SetOwner(); | |
619 | fgObjPHQ->AddAt(obji, iter); | |
620 | } | |
621 | ||
622 | //--- | |
623 | ||
624 | cfile->Close(); | |
625 | delete cfile; | |
626 | } | |
627 | ||
628 | void AliTRDdEdxUtils::DeleteCalibHist() | |
629 | { | |
630 | // | |
631 | //delete calib hist | |
632 | // | |
633 | delete fgHistGain; | |
634 | delete fgHistT0; | |
635 | delete fgHistVd; | |
636 | ||
637 | fgHistGain = 0x0; | |
638 | fgHistT0 = 0x0; | |
639 | fgHistVd = 0x0; | |
640 | ||
641 | //fgHistPHQ owns the hists | |
642 | fgHistPHQ->SetOwner(); | |
643 | fgHistPHQ->Clear(); | |
644 | } | |
645 | ||
646 | void AliTRDdEdxUtils::IniCalibHist(TList *list, const Bool_t kPHQonly) | |
647 | { | |
648 | // | |
649 | //initialize calib hist, list should not own the hist, or list->Clear/delete hist should not be called | |
650 | // | |
651 | ||
652 | DeleteCalibHist(); | |
653 | ||
654 | Int_t nbin[2]; | |
655 | const Double_t xmin[2]={0, 0}; | |
656 | Double_t xmax[2]; | |
657 | ||
658 | nbin[0]=NTRDtimebin(); nbin[1]= 11250; xmax[0]=nbin[0]; xmax[1]=20; | |
659 | for(Int_t iter=0; iter<8; iter++){ | |
660 | const TString hn = GetPHQName(0, iter); | |
661 | THnSparseD *hi = new THnSparseD(hn, "", 2, nbin, xmin, xmax); | |
662 | //fgHistPHQ owns the hists | |
663 | fgHistPHQ->AddAt(hi, iter); | |
664 | list->Add(hi); | |
665 | } | |
666 | ||
667 | if(kPHQonly) | |
668 | return; | |
669 | ||
670 | nbin[0]=NTRDchamber(); nbin[1]= 11250; xmax[0]=nbin[0]; xmax[1]=20; fgHistGain = new THnSparseD("TRDCalibHistGain", "", 2, nbin, xmin, xmax); | |
671 | nbin[0]=NTRDchamber(); nbin[1]= 11250; xmax[0]=nbin[0]; xmax[1]=AliTRDseedV1::kNtb; fgHistT0 = new THnSparseD("TRDCalibHistT0", "", 2, nbin, xmin, xmax); | |
672 | nbin[0]=NTRDchamber(); nbin[1]= 11250; xmax[0]=nbin[0]; xmax[1]=AliTRDseedV1::kNtb; fgHistVd = new THnSparseD("TRDCalibHistVd", "", 2, nbin, xmin, xmax); | |
673 | ||
674 | list->Add(fgHistGain); | |
675 | list->Add(fgHistT0); | |
676 | list->Add(fgHistVd); | |
677 | } | |
678 | ||
679 | Bool_t AliTRDdEdxUtils::ReadCalibHist(const TString filename, const TString listname) | |
680 | { | |
681 | // | |
682 | //used in AliTRDPreprocessorOffline | |
683 | //read in calib hist from file, only for PHQ | |
684 | // | |
685 | DeleteCalibHist(); | |
686 | ||
687 | //maybe already open by others... don't close | |
688 | TFile fcalib(filename); | |
689 | ||
690 | TObjArray * array = (TObjArray*)fcalib.Get(listname); | |
691 | ||
692 | for(Int_t iter=0; iter<8; iter++){ | |
693 | const TString hn = GetPHQName(0, iter); | |
694 | THnSparseD * tmph=0x0; | |
695 | if(array){ | |
696 | tmph = (THnSparseD *) array->FindObject(hn); | |
697 | } | |
698 | else{ | |
699 | tmph = (THnSparseD *) fcalib.Get(hn); | |
700 | } | |
701 | if(!tmph){ | |
702 | printf("AliTRDdEdxUtils::ReadCalibHist warning calib hist not found! %s %s\n", filename.Data(), listname.Data()); | |
703 | fcalib.ls(); | |
879b7730 | 704 | if(array){ |
705 | array->ls(); | |
706 | } | |
fb3d369e | 707 | return kFALSE; |
708 | } | |
709 | THnSparseD *hi = (THnSparseD*)tmph->Clone(); | |
710 | fgHistPHQ->AddAt(hi, iter); | |
711 | } | |
712 | ||
713 | return kTRUE; | |
714 | } | |
715 | ||
716 | void AliTRDdEdxUtils::FillCalibHist(const Int_t ncls, const TVectorD *arrayQ, const TVectorD *arrayX, THnSparseD * hcalib, const Double_t scale) | |
717 | { | |
718 | // | |
719 | //fill calibration hist | |
720 | // | |
721 | if(!hcalib){printf("AliTRDdEdxUtils::FillCalibHist errro hcalib null!!\n"); exit(1);} | |
722 | ||
723 | for(Int_t ii=0; ii<ncls; ii++){ | |
724 | const Double_t dq = (*arrayQ)[ii]; | |
725 | const Double_t xx = (*arrayX)[ii]; | |
726 | ||
727 | const Double_t qmax = hcalib->GetAxis(1)->GetXmax() -0.5 * hcalib->GetAxis(1)->GetBinWidth(1); | |
728 | const Double_t xmin = hcalib->GetAxis(0)->GetXmin(); | |
729 | const Double_t xmax = hcalib->GetAxis(0)->GetXmax(); | |
730 | ||
731 | if(xx>=xmax || xx<xmin){ | |
732 | printf("AliTRDdEdxUtils::FillCalibHist error x overflow or underflow! %s %15f %15f %15f\n", hcalib->GetName(), xx, xmin, xmax); exit(1); | |
733 | } | |
734 | ||
735 | const Double_t var[]={xx, TMath::Min(dq, qmax)/scale}; | |
736 | hcalib->Fill(var); | |
737 | } | |
738 | } | |
739 | ||
740 | void AliTRDdEdxUtils::FillCalibHist(const AliTRDtrackV1 *trdv1, const Bool_t kinvq, const Double_t mag, const Int_t charge, const Double_t scale) | |
741 | { | |
742 | // | |
743 | //get cluster Q and fill calib hist, if kinvq = kTRUE, 1/Q is filled | |
744 | // | |
745 | ||
746 | THnSparseD * hcalib = AliTRDdEdxUtils::GetHistPHQ(kinvq, mag, charge); | |
747 | ||
748 | TVectorD arrayQ(200), arrayX(200); | |
749 | const Int_t ncls = AliTRDdEdxUtils::GetArrayClusterQ(kinvq, &arrayQ, &arrayX, trdv1); | |
750 | FillCalibHist(ncls, &arrayQ, &arrayX, hcalib, kinvq ? 1/scale : scale); | |
751 | ||
752 | static Int_t kprint = 100; | |
753 | if(kprint<0){ | |
754 | printf("\nAliTRDdEdxUtils::FillCalibHist summary: \n"); | |
755 | printf("\nkinvq= %d;\n", kinvq); | |
756 | for(Int_t iq=0; iq<ncls; iq++){ | |
757 | printf("arrayX[%3d] = %15.0f; arrayQ[%3d] = %15f;\n", iq, arrayX[iq], iq, arrayQ[iq]); | |
758 | } | |
759 | printf("\n"); | |
760 | } | |
761 | kprint++; | |
762 | } | |
763 | ||
764 | Int_t AliTRDdEdxUtils::ApplyCalib(const Int_t nc0, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj) | |
765 | { | |
766 | // | |
767 | //apply calibration on arrayQ | |
768 | // | |
769 | if(!cobj){ printf("AliTRDdEdxUtils::ApplyCalib error gain array null!!\n"); exit(1);} | |
770 | ||
771 | TVectorD tmpq(arrayQ->GetNrows()); | |
772 | TVectorD tmpx(arrayX->GetNrows()); | |
773 | Int_t ncls = 0; | |
774 | ||
775 | const TVectorD * gain = (TVectorD*) cobj->At(0); | |
776 | for(Int_t ii=0; ii<nc0; ii++){ | |
777 | const Double_t dq = (*arrayQ)[ii]; | |
778 | const Int_t xx = (Int_t)(*arrayX)[ii]; | |
779 | const Double_t gg = (*gain)[xx]; | |
780 | ||
781 | if(gg<EPSILON){ | |
782 | continue; | |
783 | } | |
784 | ||
785 | tmpq[ncls] = dq*gg; | |
786 | tmpx[ncls] = xx; | |
787 | ncls++; | |
788 | } | |
789 | ||
790 | (*arrayQ)=tmpq; | |
791 | (*arrayX)=tmpx; | |
792 | ||
793 | return ncls; | |
794 | } | |
795 | ||
796 | void AliTRDdEdxUtils::GetPHCountMeanRMS(const TH1D *hnor, TH1D *&hmean) | |
797 | { | |
798 | // | |
799 | //calculate from the ph calib hist the (mean-3sigma) ph-count in the chamber, save in the TH1D output | |
800 | // | |
801 | const Int_t ndet = 540; | |
802 | TObjArray *obj=new TObjArray(ndet); | |
803 | obj->SetOwner(); | |
804 | for(Int_t ii=0; ii<ndet; ii++){ | |
805 | obj->Add(new TVectorD(AliTRDseedV1::kNtb)); | |
806 | } | |
807 | ||
808 | //ibin = binlowedge of bin(ibin+1) = the number fills this bin | |
809 | for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){ | |
810 | const Double_t stat = hnor->GetBinContent(ibin+1); | |
811 | if(stat<EPSILON){ | |
812 | continue; | |
813 | } | |
814 | ||
815 | const Int_t idet = ToDetector(ibin); | |
816 | const Int_t itb = ToTimeBin(ibin); | |
817 | TVectorD *vec=(TVectorD *)obj->At(idet); | |
818 | (*vec)[itb] = stat; | |
819 | } | |
820 | ||
821 | hmean = new TH1D(Form("%sdetmean", hnor->GetName()), "", hnor->GetNbinsX(), hnor->GetXaxis()->GetXmin(), hnor->GetXaxis()->GetXmax()); | |
822 | for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){ | |
823 | const Int_t idet = ToDetector(ibin); | |
824 | const TVectorD *vec=(TVectorD *)obj->At(idet); | |
825 | ||
826 | Int_t nonzero = 0; | |
827 | for(Int_t ii=0; ii<vec->GetNrows(); ii++){ | |
828 | if((*vec)[ii]>EPSILON){ | |
829 | nonzero++; | |
830 | } | |
831 | } | |
832 | ||
833 | Double_t mean = 0; | |
834 | const Double_t lowfrac = 0.6; | |
835 | //if there are too many 0's, reject this chamber by settig mean=rms=0 | |
836 | if(nonzero> (AliTRDseedV1::kNtb*(1-lowfrac)) ){ | |
837 | //only highest (1-lowfrac)*31 timebins are used to estimate the mean and rms! important! otherwise the 0' will make rms very large! | |
838 | mean = TruncatedMean(AliTRDseedV1::kNtb, vec->GetMatrixArray(), lowfrac, 1); | |
839 | } | |
840 | ||
841 | hmean->SetBinContent(ibin+1, mean); | |
842 | } | |
843 | ||
844 | delete obj; | |
845 | } | |
846 | ||
847 | void AliTRDdEdxUtils::CalibOutput(const TList *lin, Int_t run) | |
848 | { | |
849 | // | |
850 | //produce calibration objects | |
851 | // | |
852 | ||
853 | TString objnames("TRDCalibHistGain TRDCalibHistT0 TRDCalibHistVd "); | |
854 | for(Int_t iter=0; iter<8; iter++){ | |
855 | objnames+= GetPHQName(0, iter)+" "; | |
856 | } | |
857 | ||
858 | TList * lout = new TList; | |
859 | lout->SetOwner(); | |
860 | ||
861 | TTreeSRedirector *calibStream = new TTreeSRedirector(Form("TRDCalibStream_%010d.root", run)); | |
862 | ||
863 | const Int_t nh=lin->GetEntries(); | |
864 | for(Int_t ii=0; ii<nh; ii++){ | |
865 | const THnSparseD *hh=(THnSparseD*)lin->At(ii); | |
866 | const TString hname = hh->GetName(); | |
867 | if(!objnames.Contains(hname)) | |
868 | continue; | |
869 | ||
870 | TObjArray * cobj0 = GetCalibObj(hh, run, lout, calibStream); | |
871 | lout->Add(cobj0); | |
872 | } | |
873 | ||
874 | //lout->ls(); | |
875 | ||
876 | //============================================================= | |
877 | //============================================================= | |
878 | ||
879 | TFile *fout=new TFile(Form("TRDCalibObj_%010d.root", run),"recreate"); | |
880 | fout->cd(); | |
881 | const Int_t nout=lout->GetEntries(); | |
882 | for(Int_t ii=0; ii<nout; ii++){ | |
883 | const TString oname = lout->At(ii)->GetName(); | |
884 | if(oname.Contains("Obj")){ | |
885 | TObjArray * cobj = (TObjArray*) lout->At(ii); | |
886 | cobj->Write(oname, TObjArray::kSingleKey); | |
887 | } | |
888 | } | |
889 | fout->Save(); | |
890 | fout->Close(); | |
891 | delete fout; | |
892 | ||
893 | fout=new TFile(Form("TRDCalibList_%010d.root", run),"recreate"); | |
894 | fout->cd(); | |
895 | lin->Write(); | |
896 | lout->Write(); | |
897 | fout->Save(); | |
898 | fout->Close(); | |
899 | delete fout; | |
900 | ||
901 | delete calibStream; | |
902 | ||
903 | /* | |
904 | http://root.cern.ch/root/html/TH1.html | |
905 | When an histogram is created, a reference to it is automatically added to the list of in-memory objects for the current file or directory. This default behaviour can be changed by: | |
906 | ||
907 | h->SetDirectory(0); for the current histogram h | |
908 | TH1::AddDirectory(kFALSE); sets a global switch disabling the reference | |
909 | ||
910 | When the histogram is deleted, the reference to it is removed from the list of objects in memory. When a file is closed, all histograms in memory associated with this file are automatically deleted. | |
911 | */ | |
912 | delete lout; | |
913 | } | |
914 | ||
915 | TObjArray* AliTRDdEdxUtils::GetCalibObj(const THnSparseD *hh, Int_t run, TList *lout, TTreeSRedirector *calibStream) | |
916 | { | |
917 | // | |
918 | //produce calibration objects | |
919 | // | |
920 | ||
921 | const TString hname = hh->GetName(); | |
922 | const Bool_t kinvq = TString(hname(hname.First('Q')+1,1)).Atoi()&4; | |
923 | ||
924 | //---------------------------------------- | |
925 | // Define nbin, tag, cobj0 | |
926 | //---------------------------------------- | |
927 | Int_t nbin =-999; | |
928 | if(hname.Contains("Gain") || hname.Contains("T0") || hname.Contains("Vd")){ | |
929 | nbin = NTRDchamber(); | |
930 | } | |
931 | else if(hname.Contains("PHQ")){ | |
932 | nbin = NTRDtimebin(); | |
933 | } | |
934 | else{ | |
935 | printf("AliTRDdEdxUtils::GetCalibObj error wrong hname!! %s\n", hname.Data()); exit(1); | |
936 | } | |
937 | ||
938 | TString tag(hname); | |
939 | tag.ReplaceAll("Hist","Obj"); | |
940 | ||
941 | TObjArray * cobj0 = new TObjArray(1); | |
942 | cobj0->SetOwner(); | |
943 | cobj0->SetName(tag); | |
944 | cobj0->Add(new TVectorD(nbin)); | |
945 | ||
946 | //---------------------------------------- | |
947 | // Define lowFrac, highFrac | |
948 | //---------------------------------------- | |
949 | Double_t lowFrac = -999, highFrac = -999; | |
950 | if(hname.Contains("Gain") || (hname.Contains("PHQ") && !kinvq) ){ | |
951 | lowFrac = 0.01; highFrac = Q0Frac(); | |
952 | } | |
953 | else if(hname.Contains("PHQ") && kinvq){ | |
954 | lowFrac = Q1Frac(); highFrac = 0.99; | |
955 | } | |
956 | else{ | |
957 | lowFrac = 0.01; | |
958 | highFrac = 0.99; | |
959 | } | |
960 | ||
961 | //---------------------------------------- | |
962 | // Get analysis result | |
963 | //---------------------------------------- | |
964 | TH1::AddDirectory(kFALSE);//important! | |
965 | TH1D *hnor=0x0, *hmpv=0x0, *hres=0x0, *hwid=0x0, *htrdphmean = 0x0;//if(!lout), these have to be deleted | |
966 | TH2D *hpj = hh->Projection(1,0); | |
967 | FitSlicesY(hpj, hnor, hmpv, hwid, hres, 0, lowFrac, highFrac); | |
968 | if(hname.Contains("PHQ")){ | |
969 | GetPHCountMeanRMS(hnor, htrdphmean); | |
970 | if(lout){ | |
971 | lout->Add(htrdphmean); | |
972 | } | |
973 | } | |
974 | delete hpj; | |
975 | ||
976 | if(lout){ | |
977 | lout->Add(hnor); | |
978 | lout->Add(hmpv); | |
979 | lout->Add(hwid); | |
980 | lout->Add(hres); | |
981 | } | |
982 | ||
983 | //---------------------------------------- | |
984 | // Define Counter | |
985 | //---------------------------------------- | |
986 | TVectorD *countDet=0x0; | |
987 | TObjArray *countSSL=0x0; | |
988 | ||
989 | if(hname.Contains("PHQ") && !kinvq){ | |
990 | countDet=new TVectorD(540); | |
991 | countSSL=new TObjArray(90);//SectorStackLayer | |
992 | countSSL->SetOwner(); | |
993 | for(Int_t ii=0; ii<90; ii++){ | |
994 | countSSL->Add(new TVectorD(6)); | |
995 | } | |
996 | } | |
997 | ||
998 | //---------------------------------------- | |
999 | // Fill cobj0 | |
1000 | //---------------------------------------- | |
1001 | ||
1002 | //ibin = binlowedge of bin(ibin+1) = the number fills this bin | |
1003 | for(Int_t ibin=0; ibin<nbin; ibin++){ | |
1004 | Double_t gnor = hnor->GetBinContent(ibin+1); | |
1005 | Double_t gmpv = hmpv->GetBinContent(ibin+1); | |
1006 | Double_t gwid = hwid->GetBinContent(ibin+1); | |
1007 | Double_t gres = hres->GetBinContent(ibin+1); | |
1008 | ||
1009 | //--- set additional cut by kpass | |
1010 | Bool_t kpass = kTRUE; | |
1011 | Double_t gtrdphmean = -999; | |
1012 | if(htrdphmean){ | |
1013 | gtrdphmean = htrdphmean->GetBinContent(ibin+1); | |
1014 | //chamber no statistics (e.g. too many 0's), not usual, not seen in run 143237 | |
1015 | if(gtrdphmean<EPSILON){ | |
1016 | kpass = kFALSE; | |
1017 | } | |
1018 | if(gnor<TimeBinCountCut()*gtrdphmean){ | |
1019 | kpass = kFALSE; | |
1020 | } | |
1021 | } | |
1022 | ||
1023 | //--- set calibration constant p0 | |
1024 | Double_t p0= 0; | |
1025 | ||
1026 | //reason for gmpv=0: | |
1027 | //1)gnor<=3; truncation in hist: (0, 0.6*ntot=1.8 with ntot=3]={1}, in hist entries can pile up so that ntot=2, or 3, and (ntot>nstart && ntot<=nstop) is skipped; | |
1028 | //2)TruncatedMean(hist) out of range (only for Q0, not Q1). | |
1029 | ||
1030 | if(gmpv>EPSILON && kpass){ | |
1031 | if(tag.Contains("T0")){ | |
1032 | p0 = gmpv; | |
1033 | } | |
1034 | else{ | |
1035 | p0 = 1/gmpv; | |
1036 | } | |
1037 | //printf("outcalibobj%s %d %15.6e\n", tag.Data(), ibin, p0); | |
1038 | } | |
1039 | ||
1040 | (*( (TVectorD*)cobj0->At(0) ))[ibin] = p0; | |
1041 | ||
1042 | //--- save optional record | |
1043 | if(p0>EPSILON && countDet && countSSL){ | |
1044 | const Int_t idet = ToDetector(ibin); | |
1045 | (*countDet)[idet]=1; | |
1046 | ||
1047 | const Int_t isector = ToSector(ibin); | |
1048 | const Int_t istack = ToStack(ibin); | |
1049 | const Int_t ilayer = ToLayer(ibin); | |
1050 | TVectorD * vecsectorstack = (TVectorD*)countSSL->At(istack*18+isector); | |
1051 | (*vecsectorstack)[ilayer]=1; | |
1052 | } | |
1053 | ||
1054 | if(calibStream){ | |
1055 | (*calibStream)<<tag<< | |
1056 | "run="<<run<< | |
1057 | "p0="<<p0<< | |
1058 | ||
1059 | "nor="<<gnor<< | |
1060 | "mpv="<<gmpv<< | |
1061 | "wid="<<gwid<< | |
1062 | "res="<<gres<< | |
1063 | "gtrdphmean="<<gtrdphmean<< | |
1064 | ||
1065 | "ibin="<<ibin<< | |
1066 | "\n"; | |
1067 | } | |
1068 | } | |
1069 | ||
1070 | //---------------------------------------- | |
1071 | // Status Report | |
1072 | //---------------------------------------- | |
1073 | if(countDet && countSSL){ | |
1074 | TVectorD count2Dstack(90); | |
1075 | for(Int_t ii=0; ii<90; ii++){ | |
1076 | TVectorD * vecsectorstack = (TVectorD*)countSSL->At(ii); | |
1077 | const Int_t nlayer = (Int_t)vecsectorstack->Sum(); | |
1078 | if(nlayer==6){ | |
1079 | count2Dstack[ii]=1; | |
1080 | } | |
1081 | } | |
1082 | ||
1083 | printf("\nAliTRDdEdxUtils::GetCalibObj Summary run: %d name: %s entries: %.0f ndetector: %03.0f n2dstack %02.0f\n\n", run, hname.Data(), hh->GetEntries(), countDet->Sum(), count2Dstack.Sum()); | |
1084 | } | |
1085 | ||
1086 | //---------------------------------------- | |
1087 | // Clean Up | |
1088 | //---------------------------------------- | |
1089 | ||
1090 | TH1D **hhs[]={&hnor, &hmpv, &hwid, &hres, &htrdphmean}; | |
1091 | const Int_t nhh=sizeof(hhs)/sizeof(TH1D**); | |
1092 | for(Int_t ihh=0; ihh<nhh; ihh++){ | |
1093 | if(!lout){ | |
1094 | delete (*hhs[ihh]); | |
1095 | } | |
1096 | } | |
1097 | ||
1098 | delete countDet; | |
1099 | delete countSSL; | |
1100 | ||
1101 | //---------------------------------------- | |
1102 | ||
1103 | return cobj0; | |
1104 | } | |
1105 | ||
1106 | //=================================================================================== | |
1107 | // dEdx calculation | |
1108 | //=================================================================================== | |
fb3d369e | 1109 | Double_t AliTRDdEdxUtils::ToyCook(const Bool_t kinvq, Int_t &ncluster, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj) |
1110 | { | |
1111 | // | |
1112 | //template for cookdedx | |
1113 | // | |
1114 | if(cobj){ | |
1115 | if(arrayQ && arrayX){ | |
1116 | ncluster = ApplyCalib(ncluster, arrayQ, arrayX, cobj); | |
1117 | } | |
1118 | else{ | |
1119 | printf("AliTRDdEdxUtils::ToyCook arrayQ arrayX null, applycalib can not be applied!\n"); exit(1); | |
1120 | } | |
1121 | } | |
1122 | ||
1123 | Double_t lowFrac =-999, highFrac = -999; | |
1124 | if(kinvq){ | |
1125 | lowFrac = Q1Frac(); highFrac = 0.99; | |
1126 | } | |
1127 | else{ | |
1128 | lowFrac = 0.01; highFrac = Q0Frac(); | |
1129 | } | |
1130 | ||
1131 | Double_t meanQ = TruncatedMean(ncluster, arrayQ->GetMatrixArray(), lowFrac, highFrac); | |
1132 | if(kinvq){ | |
1133 | if(meanQ>EPSILON){ | |
1134 | meanQ = 1/meanQ; | |
1135 | } | |
1136 | } | |
1137 | ||
1138 | return meanQ; | |
1139 | } | |
1140 | ||
1141 | Double_t AliTRDdEdxUtils::CombineddEdx(const Bool_t kinvq, Int_t &concls, TVectorD *coarrayQ, TVectorD *coarrayX, const Int_t tpcncls, const TVectorD *tpcarrayQ, const TVectorD *tpcarrayX, const Int_t trdncls, const TVectorD *trdarrayQ, const TVectorD *trdarrayX) | |
1142 | { | |
1143 | // | |
1144 | //combine tpc and trd dedx | |
1145 | // | |
1146 | ||
1147 | for(Int_t iq=0; iq<tpcncls; iq++){ | |
1148 | (*coarrayQ)[iq]=(*tpcarrayQ)[iq]; | |
1149 | if(tpcarrayX && trdarrayX && coarrayX){ | |
1150 | (*coarrayX)[iq]=(*tpcarrayX)[iq]; | |
1151 | } | |
1152 | } | |
1153 | for(Int_t iq=0; iq<trdncls; iq++){ | |
1154 | (*coarrayQ)[tpcncls+iq]=(*trdarrayQ)[iq]; | |
1155 | if(tpcarrayX && trdarrayX && coarrayX){ | |
1156 | (*coarrayX)[tpcncls+iq]=159+(*trdarrayX)[iq]; | |
1157 | } | |
1158 | } | |
1159 | ||
1160 | concls=trdncls+tpcncls; | |
1161 | ||
1162 | const Double_t coQ = ToyCook(kinvq, concls, coarrayQ, coarrayX); | |
1163 | ||
1164 | return coQ; | |
1165 | } | |
1166 | ||
1167 | ||
1168 | //=================================================================================== | |
1169 | // dEdx Getter and Setter | |
1170 | //=================================================================================== | |
1171 | Double_t AliTRDdEdxUtils::GetAngularCorrection(const AliTRDseedV1 *seed) | |
1172 | { | |
1173 | // | |
1174 | //return angular normalization factor | |
1175 | // | |
1176 | ||
1177 | return TMath::Sqrt(1+seed->GetYref(1)*seed->GetYref(1)+seed->GetZref(1)*seed->GetZref(1)); | |
1178 | } | |
1179 | ||
d0037ea4 | 1180 | Double_t AliTRDdEdxUtils::GetPadGain(const Int_t det, const Int_t icol, const Int_t irow) |
1181 | { | |
1182 | // | |
1183 | //get pad calibration | |
1184 | // | |
1185 | AliTRDcalibDB* calibration = AliTRDcalibDB::Instance(); | |
1186 | if (!calibration) { | |
1187 | printf("AliTRDdEdxUtils::GetPadCalib No AliTRDcalibDB instance available\n"); exit(1); | |
1188 | } | |
1189 | AliTRDCalROC * calGainFactorROC = calibration->GetGainFactorROC(det); | |
1190 | if(!calGainFactorROC){ | |
1191 | printf("AliTRDdEdxUtils::GetPadCalib no calGainFactorROC!\n"); exit(1); | |
1192 | } | |
1193 | ||
1194 | Double_t padgain = -999; | |
1195 | if( icol >= 0 && | |
1196 | icol < calGainFactorROC->GetNcols() && | |
1197 | irow >=0 && | |
1198 | irow < calGainFactorROC->GetNrows()){ | |
1199 | padgain = calGainFactorROC->GetValue(icol, irow); | |
1200 | if(padgain<EPSILON){ | |
1201 | printf("AliTRDdEdxUtils::GetPadGain padgain 0! %f %f -- %d %d %d -- %d %d\n", padgain, EPSILON, det, icol, irow, calGainFactorROC->GetNcols(), calGainFactorROC->GetNrows()); exit(1); | |
1202 | } | |
1203 | } | |
1204 | else{ | |
1205 | //printf("\nAliTRDdEdxUtils::GetPadGain warning!! indices out of range %d %d %d -- %d %d\n\n", det, icol, irow, calGainFactorROC->GetNcols(), calGainFactorROC->GetNrows() ); | |
1206 | } | |
1207 | ||
1208 | return padgain; | |
1209 | } | |
1210 | ||
1211 | Double_t AliTRDdEdxUtils::GetRNDClusterQ(AliTRDcluster *cl) | |
1212 | { | |
1213 | // | |
1214 | //get cluter q from GetRawQ, apply baseline and Kr pad-calibration | |
1215 | // | |
1216 | ||
1217 | const Int_t det = cl->GetDetector(); | |
1218 | const Int_t pad3col = cl->GetPadCol(); | |
1219 | const Int_t padrow = cl->GetPadRow(); | |
1220 | ||
1221 | const Double_t baseline = 10; | |
1222 | ||
1223 | Double_t rndqsum = 0; | |
1224 | for(Int_t ii=0; ii<7; ii++){ | |
25d9c307 | 1225 | if(cl->GetSignals()[ii] < EPSILON){//bad pad marked by electronics |
1226 | continue; | |
1227 | } | |
d0037ea4 | 1228 | |
25d9c307 | 1229 | const Int_t icol = pad3col+(ii-3); |
1230 | const Double_t padgain = GetPadGain(det, icol, padrow); | |
1231 | if(padgain<0){//indices out of range, pad3col near boundary case | |
d0037ea4 | 1232 | continue; |
1233 | } | |
1234 | ||
25d9c307 | 1235 | const Double_t rndsignal = (cl->GetSignals()[ii] - baseline)/(IsPadGainOn()? padgain : 1); |
d0037ea4 | 1236 | |
1237 | //sum it anyway even if signal below baseline, as long as the total is positive | |
1238 | rndqsum += rndsignal; | |
1239 | } | |
1240 | ||
1241 | return rndqsum; | |
1242 | } | |
1243 | ||
fb3d369e | 1244 | Double_t AliTRDdEdxUtils::GetClusterQ(const Bool_t kinvq, const AliTRDseedV1 * seed, const Int_t itb) |
1245 | { | |
1246 | // | |
1247 | //get cluster charge | |
1248 | // | |
1249 | Double_t dq = 0; | |
d0037ea4 | 1250 | AliTRDcluster *cl = 0x0; |
fb3d369e | 1251 | |
25d9c307 | 1252 | //GetRNDClusterQ(cl)>0 ensures that the total sum of q is above baseline*NsignalPhysical. |
d0037ea4 | 1253 | cl = seed->GetClusters(itb); if(cl && GetRNDClusterQ(cl)>0 ) dq+= GetRNDClusterQ(cl);//cl->GetRawQ(); |
1254 | cl = seed->GetClusters(itb+AliTRDseedV1::kNtb); if(cl && GetRNDClusterQ(cl)>0 ) dq+= GetRNDClusterQ(cl);//cl->GetRawQ(); | |
fb3d369e | 1255 | |
1256 | dq /= GetAngularCorrection(seed); | |
1257 | ||
d0037ea4 | 1258 | dq /= QScale(); |
fb3d369e | 1259 | |
1260 | if(kinvq){ | |
1261 | if(dq>EPSILON){ | |
1262 | dq = 1/dq; | |
1263 | } | |
1264 | } | |
1265 | ||
1266 | return dq; | |
1267 | } | |
1268 | ||
1269 | Int_t AliTRDdEdxUtils::GetArrayClusterQ(const Bool_t kinvq, TVectorD *arrayQ, TVectorD *arrayX, const AliTRDtrackV1 *trdtrack, Int_t timeBin0, Int_t timeBin1, Int_t tstep) | |
1270 | { | |
1271 | // | |
1272 | //return nclustter | |
1273 | //(if kinvq, return 1/q array), size of array must be larger than 31*6 | |
1274 | // | |
1275 | if(!arrayQ || arrayQ->GetNrows()< (AliTRDseedV1::kNtb*AliTRDtrackV1::kNplane)){ | |
1276 | printf("AliTRDdEdxUtils::GetArrayClusterQ error arrayQ null or size too small! %d\n", arrayQ? arrayQ->GetNrows() : -999); exit(1); | |
1277 | } | |
879b7730 | 1278 | if(!arrayX || arrayX->GetNrows()< (AliTRDseedV1::kNtb*AliTRDtrackV1::kNplane)){ |
fb3d369e | 1279 | printf("AliTRDdEdxUtils::GetArrayClusterQ error arrayX null or size too small! %d\n", arrayX? arrayX->GetNrows() : -999); exit(1); |
1280 | } | |
1281 | ||
1282 | const Int_t mintb = 0; | |
1283 | const Int_t maxtb = AliTRDseedV1::kNtb-1; | |
1284 | if(timeBin0<mintb) timeBin0=mintb; | |
1285 | if(timeBin1>maxtb) timeBin1=maxtb; | |
1286 | if(tstep<=0) tstep=1; | |
1287 | ||
1288 | //============ | |
1289 | Int_t tbN=0; | |
1290 | Double_t tbQ[200]; | |
1291 | Int_t tbBin[200]; | |
1292 | ||
1293 | for(Int_t ichamber=0; ichamber < AliTRDtrackV1::kNplane; ichamber++){ | |
1294 | const AliTRDseedV1 * seed = trdtrack->GetTracklet(ichamber); | |
1295 | if(!seed) | |
1296 | continue; | |
1297 | ||
1298 | const Int_t det = seed->GetDetector(); | |
1299 | ||
1300 | for(Int_t itb=timeBin0; itb<=timeBin1; itb+=tstep){ | |
1301 | const Double_t dq = GetClusterQ(kinvq, seed, itb); | |
1302 | if(dq<EPSILON) | |
1303 | continue; | |
1304 | ||
1305 | const Int_t gtb = det * AliTRDseedV1::kNtb + itb; | |
1306 | ||
1307 | tbQ[tbN]=dq; | |
1308 | tbBin[tbN]=gtb; | |
1309 | tbN++; | |
1310 | } | |
1311 | } | |
1312 | ||
1313 | Int_t ncls = 0; | |
1314 | for(Int_t iq=0; iq<tbN; iq++){ | |
1315 | if(tbQ[iq]<EPSILON) | |
1316 | continue; | |
1317 | ||
1318 | (*arrayQ)[ncls] = tbQ[iq]; | |
1319 | (*arrayX)[ncls] = tbBin[iq]; | |
1320 | ||
1321 | ncls++; | |
1322 | } | |
1323 | ||
879b7730 | 1324 | static Int_t kprint = 100; |
fb3d369e | 1325 | if(kprint<0){ |
1326 | printf("\nAliTRDdEdxUtils::GetArrayClusterQ raw cluster-Q\n"); | |
1327 | for(Int_t iq=0; iq<ncls; iq++){ | |
1328 | const Int_t ichamber = ToLayer((*arrayX)[iq]); | |
1329 | const AliTRDseedV1 * seed = trdtrack->GetTracklet(ichamber); | |
1330 | if(!seed){ | |
1331 | printf("error seed null!!\n"); exit(1); | |
1332 | } | |
1333 | const Double_t rawq = (*arrayQ)[iq] * 45. * GetAngularCorrection(seed); | |
1334 | printf("esdid=%d; chamber=%d; timebin=%d; rawq= %.3f; myq[%d]= %e;\n", trdtrack->GetESDid(), ichamber, ToTimeBin((*arrayX)[iq]), rawq, iq, (*arrayQ)[iq]); | |
1335 | } | |
1336 | printf("\n"); | |
1337 | } | |
879b7730 | 1338 | kprint++; |
fb3d369e | 1339 | |
1340 | return ncls; | |
1341 | } | |
1342 | ||
1343 | Int_t AliTRDdEdxUtils::UpdateArrayX(const Int_t ncls, TVectorD* arrayX) | |
1344 | { | |
1345 | // | |
1346 | //arrayX det*Ntb+itb -> itb | |
1347 | // | |
1348 | ||
1349 | TVectorD countChamber(6); | |
1350 | for(Int_t ii = 0; ii<ncls; ii++){ | |
1351 | const Int_t xx = (Int_t)(*arrayX)[ii]; | |
1352 | const Int_t idet = ToDetector(xx); | |
1353 | ||
1354 | const Double_t ich = AliTRDgeometry::GetLayer(idet); | |
1355 | const Double_t itb = ToTimeBin(xx); | |
1356 | (*arrayX)[ii] = ich*AliTRDseedV1::kNtb+itb; | |
1357 | ||
1358 | countChamber[ich] = 1; | |
1359 | } | |
1360 | ||
1361 | const Double_t nch = countChamber.Sum(); | |
1362 | return (Int_t) nch; | |
1363 | } | |
1364 | ||
1365 | void AliTRDdEdxUtils::SetChamberQT(const AliTRDtrackV1 *trdtrack, const Int_t kcalib, THnSparseD * hgain, THnSparseD * ht0, THnSparseD * hvd) | |
1366 | { | |
1367 | // | |
1368 | //CookdEdx at TRD track level, use chamber info, related calibrations: chamber-gain; T0, Vd based on raw PH distribution | |
1369 | // | |
1370 | ||
1371 | static Int_t kprint = 100; | |
1372 | ||
1373 | fgNchamber = 0; | |
1374 | for(Int_t ichamber=0; ichamber < AliTRDtrackV1::kNplane; ichamber++){ | |
1375 | //initialize output, default values: 0, so that summation and weighting will automatically discard default quantities | |
1376 | fgChamberQ[ichamber] = fgChamberTmean[ichamber] = 0; | |
1377 | ||
1378 | const AliTRDseedV1 *seed = trdtrack->GetTracklet(ichamber); | |
1379 | if (!seed) | |
1380 | continue; | |
1381 | ||
1382 | const Int_t idet = seed->GetDetector(); | |
1383 | ||
1384 | //------------------------------------------------------------------------- | |
1385 | ||
1386 | Double_t qsum = 0, qtsum = 0, w2sum = 0; | |
1387 | for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++){ | |
1388 | const Double_t dq = GetClusterQ(0, seed, itb); | |
1389 | if(dq<EPSILON) | |
1390 | continue; | |
1391 | ||
1392 | qsum += dq; | |
1393 | qtsum += dq*itb; | |
1394 | w2sum += dq*itb*itb; | |
1395 | } | |
1396 | if(qsum<EPSILON) | |
1397 | continue; | |
1398 | ||
1399 | //------------------------------------------------------------------------- | |
1400 | ||
1401 | Double_t tbm, tbr = 0; | |
1402 | tbm = GetMeanRMS(qsum, qtsum, w2sum, &tbr); | |
1403 | ||
1404 | qsum /= 1.25e3/45.; | |
1405 | ||
1406 | if(hgain){ | |
1407 | const Double_t var[]={idet, qsum}; | |
1408 | hgain->Fill(var); | |
1409 | } | |
1410 | if(ht0){ | |
1411 | const Double_t var[]={idet, tbm}; | |
1412 | ht0->Fill(var); | |
1413 | } | |
1414 | if(hvd){ | |
1415 | const Double_t var[]={idet, tbr}; | |
1416 | hvd->Fill(var); | |
1417 | } | |
1418 | ||
1419 | Double_t gain = 1, t0 = 0, vd = 1; | |
1420 | if(kcalib){ | |
1421 | if(!fgObjGain) {printf("AliTRDdEdxUtils::SetChamberQT error Gain array null!!\n"); exit(1);} | |
1422 | if(! fgObjT0) {printf("AliTRDdEdxUtils::SetChamberQT error T0 array null!!\n"); exit(1);} | |
1423 | if(! fgObjVd) {printf("AliTRDdEdxUtils::SetChamberQT error Vd array null!!\n"); exit(1);} | |
1424 | ||
1425 | const TVectorD * gainvec = (TVectorD*) fgObjGain->At(0); gain = (*gainvec)[idet]; | |
1426 | const TVectorD * t0vec = (TVectorD*) fgObjT0->At(0); t0 = (* t0vec)[idet]; | |
1427 | const TVectorD * vdvec = (TVectorD*) fgObjVd->At(0); vd = (* vdvec)[idet]; | |
1428 | } | |
1429 | if(kprint<0){ | |
1430 | printf("\nAliTRDdEdxUtils::CookdEdxV2\n"); | |
1431 | printf("idet = %d;\n", idet); | |
1432 | printf("gain = %15f; t0 = %15f; vd = %15f;\n", gain, t0, vd); | |
1433 | printf("\n"); | |
1434 | } | |
1435 | ||
1436 | qsum *= gain; | |
1437 | tbm = (tbm-t0)*vd; | |
1438 | ||
1439 | if(qsum<EPSILON) | |
1440 | continue; | |
1441 | ||
1442 | //------------------------------------------------------------------------- | |
1443 | ||
1444 | //should have non-zero value, initialized with default 0 (except for calibrated tbm, may be very close to 0) | |
1445 | fgChamberQ[ichamber] = qsum; | |
1446 | fgChamberTmean[ichamber] = tbm; | |
1447 | fgNchamber++; | |
1448 | } | |
1449 | ||
1450 | if(kprint<0){ | |
1451 | printf("\nAliTRDdEdxUtils::CookdEdxV2 summary:\n"); | |
1452 | ||
1453 | printf("\nfgNchamber = %d\n", fgNchamber); | |
1454 | for(Int_t ich=0; ich<AliTRDtrackV1::kNplane; ich++){ | |
1455 | printf("fgChamberTmean[%d] = %15f; fgChamberQ[%d] = %15f;\n", ich, fgChamberTmean[ich], ich, fgChamberQ[ich]); | |
1456 | } | |
1457 | } | |
1458 | ||
1459 | fgTrackTmean = -999; | |
1460 | if(fgNchamber){ | |
1461 | fgTrackTmean = 0; | |
1462 | for(Int_t ich=0; ich<AliTRDtrackV1::kNplane; ich++){ | |
1463 | fgTrackTmean += fgChamberTmean[ich]; | |
1464 | } | |
1465 | fgTrackTmean /= fgNchamber; | |
1466 | } | |
1467 | ||
1468 | if(kprint<0){ | |
1469 | printf("\nAliTRDdEdxUtils::CookdEdxV2\n"); | |
1470 | printf("GetTrackTmean() %15f\n", GetTrackTmean()); | |
1471 | printf("\n"); | |
1472 | } | |
1473 | kprint++; | |
1474 | ||
1475 | return; | |
1476 | } | |
1477 | ||
1478 | ||
1479 | //=================================================================================== | |
1480 | // dEdx Parameterization | |
1481 | //=================================================================================== | |
1482 | ||
1483 | Double_t AliTRDdEdxUtils::Q0MeanTRDpp(const Double_t bg) | |
1484 | { | |
1485 | // | |
1486 | //truncated Mean Q_{xx} in TRD | |
1487 | // | |
1488 | ||
1489 | Double_t par[8]; | |
1490 | //03132012161150 | |
1491 | //opt: ppQ0 | |
1492 | par[0]= 2.397001e-01; | |
1493 | par[1]= 1.334697e+00; | |
1494 | par[2]= 6.967470e+00; | |
1495 | par[3]= 9.055289e-02; | |
1496 | par[4]= 9.388760e+00; | |
1497 | par[5]= 9.452742e-04; | |
1498 | par[6]= -1.866091e+00; | |
1499 | par[7]= 1.403545e+00; | |
1500 | ||
1501 | ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see2.log:hhtype2Q0b2c2 scale 0.428543 at ltbg 0.650000 | |
1502 | return 0.428543*MeandEdxTR(&bg, par); | |
1503 | } | |
1504 | ||
1505 | Double_t AliTRDdEdxUtils::Q0MeanTRDPbPb(const Double_t bg) | |
1506 | { | |
1507 | // | |
1508 | //truncated Mean Q_{xx} in TRD | |
1509 | // | |
1510 | ||
1511 | Double_t par[8]; | |
1512 | ||
1513 | //03132012161259 | |
1514 | //opt: PbPbQ0 | |
1515 | par[0]= 1.844912e-01; | |
1516 | par[1]= 2.509702e+00; | |
1517 | par[2]= 6.744031e+00; | |
1518 | par[3]= 7.355123e-02; | |
1519 | par[4]= 1.166023e+01; | |
1520 | par[5]= 1.736186e-04; | |
1521 | par[6]= -1.716063e+00; | |
1522 | par[7]= 1.611366e+00; | |
1523 | ||
1524 | ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see4.log:hhtype4Q0b2c2 scale 0.460994 at ltbg 0.650000 | |
1525 | return 0.460994*MeandEdxTR(&bg, par); | |
1526 | } | |
1527 | ||
1528 | Double_t AliTRDdEdxUtils::Q1MeanTRDpp(const Double_t bg) | |
1529 | { | |
1530 | // | |
1531 | //truncated Mean 1/(1/Q)_{xx} in TRD | |
1532 | // | |
1533 | ||
1534 | Double_t par[8]; | |
1535 | ||
1536 |