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