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