]>
Commit | Line | Data |
---|---|---|
6951a056 | 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 | // TRD dEdx base utils | |
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 | #include "TF1.h" | |
30 | #include "TFile.h" | |
31 | #include "TH1D.h" | |
32 | #include "TH2D.h" | |
33 | #include "THnSparse.h" | |
34 | #include "TMath.h" | |
35 | #include "TMatrixD.h" | |
36 | #include "TMinuit.h" | |
37 | #include "TObjArray.h" | |
38 | #include "TRandom3.h" | |
39 | #include "TStopwatch.h" | |
40 | #include "TVectorD.h" | |
41 | ||
42 | #include "TTreeStream.h" | |
43 | ||
6951a056 | 44 | #include "AliESDEvent.h" |
45 | #include "AliESDfriendTrack.h" | |
46 | #include "AliESDtrack.h" | |
6951a056 | 47 | #include "AliTRDtrackV1.h" |
48 | ||
49 | #include "AliTRDdEdxBaseUtils.h" | |
50 | ||
6fa94550 | 51 | #define EPSILON 1e-8 |
6951a056 | 52 | |
687aa844 | 53 | Double_t AliTRDdEdxBaseUtils::fgQ0Frac = 0.55; |
6951a056 | 54 | Double_t AliTRDdEdxBaseUtils::fgQ1Frac = 0.5; |
55 | Double_t AliTRDdEdxBaseUtils::fgTimeBinCountCut = 0.0; | |
56 | Int_t AliTRDdEdxBaseUtils::fgCalibTPCnclsCut = 70; | |
57 | Bool_t AliTRDdEdxBaseUtils::fgExBOn = kTRUE; | |
58 | Bool_t AliTRDdEdxBaseUtils::fgPadGainOn = kTRUE; | |
687aa844 | 59 | Double_t AliTRDdEdxBaseUtils::fgQScale = 50; |
6951a056 | 60 | |
61 | //=================================================================================== | |
62 | // Math and Histogram | |
63 | //=================================================================================== | |
6856b6a8 | 64 | void AliTRDdEdxBaseUtils::BinLogX(TAxis *axis) |
65 | { | |
66 | // | |
67 | // Method for the correct logarithmic binning of histograms | |
68 | // copied and modified from AliTPCcalibBase | |
69 | ||
70 | const Int_t bins = axis->GetNbins(); | |
71 | ||
72 | const Double_t from = axis->GetXmin(); | |
73 | const Double_t to = axis->GetXmax(); | |
ca66da8f | 74 | if (from<EPSILON){ |
396eba11 | 75 | //printf("AliTRDdEdxBaseUtils::BinLogX warning xmin < epsilon! nothing done, axis not set. %e\n", from); exit(1); |
76 | return; | |
ca66da8f | 77 | } |
78 | ||
6856b6a8 | 79 | Double_t *new_bins = new Double_t[bins + 1]; |
80 | ||
81 | new_bins[0] = from; | |
82 | const Double_t factor = pow(to/from, 1./bins); | |
83 | ||
84 | for (int i = 1; i <= bins; i++) { | |
85 | new_bins[i] = factor * new_bins[i-1]; | |
86 | } | |
87 | axis->Set(bins, new_bins); | |
88 | delete [] new_bins; | |
89 | } | |
90 | ||
6951a056 | 91 | void AliTRDdEdxBaseUtils::GetCDFCuts(const TH1D *hh, const Int_t ncut, Double_t cuts[], const Double_t cdfs[], const Double_t thres) |
92 | { | |
93 | // | |
94 | //counts of hh is sorted | |
95 | // | |
96 | ||
97 | for(Int_t ii=0; ii<ncut; ii++){ | |
98 | cuts[ii] = -999; | |
99 | } | |
100 | ||
101 | Int_t nsel = 0; | |
102 | const Int_t nbin = hh->GetNbinsX(); | |
103 | Double_t datas[nbin]; | |
104 | for(Int_t ii=1; ii<=nbin; ii++){ | |
105 | const Double_t res = hh->GetBinContent(ii); | |
106 | if(res<thres){ | |
107 | continue; | |
108 | } | |
109 | ||
110 | datas[nsel] = res; | |
111 | nsel++; | |
112 | } | |
113 | if(!nsel) | |
114 | return; | |
115 | ||
116 | Int_t id[nsel]; | |
117 | TMath::Sort(nsel, datas, id, kFALSE); | |
118 | ||
119 | for(Int_t ii=0; ii<ncut; ii++){ | |
120 | const Double_t icdf = cdfs[ii]; | |
121 | if(icdf<0 || icdf>1){ | |
122 | printf("AliTRDdEdxBaseUtils::GetCDFCuts error cdfs[%d] %15f out of range!\n", ii, icdf); exit(1); | |
123 | } | |
124 | cuts[ii] = datas[id[Int_t(icdf*nsel)]]; | |
125 | } | |
126 | } | |
127 | ||
128 | Double_t AliTRDdEdxBaseUtils::GetMeanRMS(const Double_t nn, const Double_t sum, const Double_t w2s, Double_t * grms, Double_t * gerr) | |
129 | { | |
130 | // | |
131 | //calculate mean (with error) and rms from sum, w2s, nn | |
132 | //if nn=0, mean, error, and rms all = 0 | |
133 | // | |
134 | ||
135 | Double_t tmean = 0, trms = 0, terr = 0; | |
136 | ||
137 | if(nn>EPSILON){ | |
138 | tmean = sum/nn; | |
139 | ||
140 | const Double_t arg = w2s/nn-tmean*tmean; | |
141 | if(TMath::Abs(arg)<EPSILON){ | |
142 | trms = 0; | |
143 | } | |
144 | else{ | |
145 | if( arg <0 ){ | |
146 | printf("AliTRDdEdxBaseUtils::GetMeanRMS error negative sqrt argument!! %e -- %e %e %f\n", arg, w2s, sum, nn); exit(1); | |
147 | } | |
148 | ||
149 | trms = TMath::Sqrt(arg); | |
150 | } | |
151 | ||
152 | terr = trms/TMath::Sqrt(nn); | |
153 | } | |
154 | ||
155 | if(grms){ | |
156 | (*grms) = trms; | |
157 | } | |
158 | ||
159 | if(gerr){ | |
160 | (*gerr) = terr; | |
161 | } | |
162 | ||
163 | return tmean; | |
164 | } | |
165 | ||
166 | Double_t AliTRDdEdxBaseUtils::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) | |
167 | { | |
168 | // | |
169 | //calculate truncated mean | |
170 | //return <x*w>_{low-high according to x} | |
171 | // | |
172 | ||
173 | /* | |
174 | //test-> | |
175 | for(Int_t ii=0; ii<nx; ii++){ | |
176 | printf("test %d/%d %f\n", ii, nx, xdata[ii]); | |
177 | } | |
178 | //<--test | |
179 | */ | |
180 | ||
181 | Int_t index[nx]; | |
182 | TMath::Sort(nx, xdata, index, kFALSE); | |
183 | ||
184 | Int_t nused = 0; | |
185 | Double_t sum = 0; | |
186 | Double_t w2s = 0; | |
187 | const Int_t istart = Int_t (nx*lowfrac); | |
188 | const Int_t istop = Int_t (nx*highfrac); | |
189 | ||
190 | //=,< correct, because when low=0, high=1 it is correct | |
191 | for(Int_t ii=istart; ii<istop; ii++){ | |
192 | Double_t weight = 1; | |
193 | if(wws){ | |
194 | weight = wws[index[ii]]; | |
195 | } | |
196 | const Double_t sx = xdata[index[ii]]*weight; | |
197 | ||
198 | sum += sx; | |
199 | w2s += sx*sx; | |
200 | ||
201 | nused++; | |
202 | //printf("test in loop %d/%d %f %f %f\n", ii, nused, sx, sum, w2s); | |
203 | ||
204 | } | |
205 | ||
206 | return GetMeanRMS(nused, sum, w2s, grms, gerr); | |
207 | } | |
208 | ||
209 | Double_t AliTRDdEdxBaseUtils::TruncatedMean(const TH1 *hh, const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr) | |
210 | { | |
211 | // | |
212 | //do truncation on histogram | |
213 | // | |
214 | //if hh is scaled, be sure Sumw2 is called before scaling!! then mean, rms and err will all be correct | |
215 | ||
216 | //with under-/over-flow | |
217 | Double_t npreTrunc = 0; | |
218 | for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){ | |
219 | const Double_t be = hh->GetBinError(itmp); | |
220 | const Double_t bc = hh->GetBinContent(itmp); | |
221 | if(be<EPSILON){ | |
222 | if(bc>EPSILON){ | |
223 | printf("AliTRDdEdxBaseUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1); | |
224 | } | |
225 | continue; | |
226 | } | |
227 | npreTrunc += bc*bc/be/be; | |
228 | } | |
229 | ||
230 | const Double_t nstart = npreTrunc*lowfrac; | |
231 | const Double_t nstop = npreTrunc*highfrac; | |
232 | ||
233 | //with Double_t this should also handle normalized hist | |
234 | Double_t ntot = 0; | |
235 | Double_t nused = 0; | |
236 | Double_t sum = 0; | |
237 | Double_t w2s = 0; | |
238 | for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){ | |
239 | const Double_t be = hh->GetBinError(itmp); | |
240 | const Double_t bc = hh->GetBinContent(itmp); | |
241 | if(be<EPSILON){ | |
242 | if(bc>EPSILON){ | |
243 | printf("AliTRDdEdxBaseUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1); | |
244 | } | |
245 | continue; | |
246 | } | |
247 | const Double_t weight = bc*bc/be/be; | |
248 | ntot+=weight; | |
249 | //<= correct, because when high=1, nstop has to be included | |
250 | if(ntot>nstart && ntot<=nstop){ | |
251 | ||
252 | const Double_t val = hh->GetBinCenter(itmp); | |
253 | sum += weight*val; | |
254 | w2s += weight*val*val; | |
255 | ||
256 | nused += weight; | |
257 | ||
258 | //printf("test %d %f %f --- %f %f -- %f %f\n", itmp, weight, val, sum, w2s, nused, nsample); | |
259 | } | |
260 | else if(ntot>nstop){ | |
261 | if(itmp>=hh->GetNbinsX()){ | |
262 | printf("AliTRDdEdxBaseUtils::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); | |
263 | nused = 0; | |
264 | w2s = sum = 0; | |
265 | } | |
266 | break; | |
267 | } | |
268 | } | |
269 | ||
270 | return GetMeanRMS(nused, sum, w2s, grms, gerr); | |
271 | } | |
272 | ||
273 | void AliTRDdEdxBaseUtils::FitSlicesY(const TH2D *hh, TH1D *&hnor, TH1D *&hmpv, TH1D *&hwid, TH1D *&hres, const Double_t thres, const Double_t lowfrac, const Double_t highfrac) | |
274 | { | |
275 | // | |
276 | //fit slices of hh using truncation | |
277 | // | |
278 | ||
279 | const Int_t x0 = hh->GetXaxis()->GetFirst(); | |
280 | const Int_t x1 = hh->GetXaxis()->GetLast(); | |
281 | const Int_t y0 = hh->GetYaxis()->GetFirst(); | |
282 | const Int_t y1 = hh->GetYaxis()->GetLast(); | |
283 | ||
284 | const Int_t nx = hh->GetNbinsX(); | |
285 | const Int_t ny = hh->GetNbinsY(); | |
286 | const Double_t xmin = hh->GetXaxis()->GetXmin(); | |
287 | const Double_t xmax = hh->GetXaxis()->GetXmax(); | |
288 | const Double_t ymin = hh->GetYaxis()->GetXmin(); | |
289 | const Double_t ymax = hh->GetYaxis()->GetXmax(); | |
290 | ||
291 | hnor = new TH1D(Form("%s_amp",hh->GetName()), "", nx, xmin, xmax); | |
292 | hmpv = new TH1D(Form("%s_mpv",hh->GetName()), "", nx, xmin, xmax); | |
293 | hwid = new TH1D(Form("%s_wid",hh->GetName()), "", nx, xmin, xmax); | |
294 | hres = new TH1D(Form("%s_res",hh->GetName()), "", nx, xmin, xmax); | |
295 | ||
687aa844 | 296 | Double_t newbins[ny+1]; |
297 | for(Int_t ii=1; ii<=ny+1; ii++){ | |
298 | const Double_t lowx= hh->GetYaxis()->GetBinLowEdge(ii); | |
299 | newbins[ii-1]=lowx; | |
300 | } | |
301 | ||
6951a056 | 302 | for(Int_t ix=x0; ix<=x1; ix++){ |
303 | //to speed up | |
304 | const Double_t rawcount = hh->Integral(ix,ix,0, ny+1); | |
305 | if(rawcount<EPSILON){ | |
306 | continue; | |
307 | } | |
308 | ||
309 | TH1D *htmp = new TH1D(Form("FitSlicesY_%s_%d", hh->GetName(), ix),"",ny, ymin, ymax); | |
687aa844 | 310 | htmp->GetXaxis()->Set(ny, newbins); |
6951a056 | 311 | Double_t ntot = 0; |
312 | for(Int_t iy=y0; iy<=y1; iy++){ | |
313 | const Double_t be = hh->GetBinError(ix,iy); | |
314 | const Double_t bc = hh->GetBinContent(ix, iy); | |
315 | ||
316 | if(be<EPSILON){ | |
317 | if(bc>EPSILON){ | |
318 | printf("AliTRDdEdxBaseUtils::FitSlicesY error %d %d %e %e\n", ix, iy, be, bc); exit(1); | |
319 | } | |
320 | continue; | |
321 | } | |
322 | ||
323 | htmp->SetBinContent(iy, bc); | |
324 | htmp->SetBinError(iy, be); | |
325 | ||
326 | ntot += (bc/be)*(bc/be); | |
327 | ||
328 | //if(be) printf("test %d %d : %f %f %f\n", ix, iy, bc, be, pow(bc/be,2)); | |
329 | } | |
330 | ||
331 | hnor->SetBinContent(ix, ntot); | |
332 | hnor->SetBinError( ix, 0); | |
333 | ||
334 | if(ntot<thres || htmp->GetRMS()<EPSILON){ | |
335 | delete htmp; | |
336 | continue; | |
337 | } | |
338 | ||
339 | //test htmp->Draw(); | |
340 | Double_t trms = -999, terr = -999; | |
341 | const Double_t tmean = TruncatedMean(htmp, lowfrac, highfrac, &trms, &terr); | |
342 | ||
343 | hmpv->SetBinContent(ix, tmean); | |
344 | hmpv->SetBinError( ix, terr); | |
345 | ||
346 | hwid->SetBinContent(ix, trms); | |
347 | hwid->SetBinError( ix, 0); | |
348 | ||
349 | hres->SetBinContent(ix, tmean>EPSILON ? trms/tmean:0); | |
350 | hres->SetBinError( ix, 0); | |
351 | ||
352 | delete htmp; | |
353 | } | |
354 | ||
355 | TH1 *hhs[]={hnor, hmpv, hwid, hres}; | |
356 | const TString yt[]={"N", "MPV", "#sigma", "#sigma/MPV"}; | |
357 | const Int_t nh = sizeof(hhs)/sizeof(TH1*); | |
358 | for(Int_t ii=0; ii<nh; ii++){ | |
359 | hhs[ii]->SetYTitle(Form("%s of %s", yt[ii].Data(), hh->GetYaxis()->GetTitle())); | |
360 | hhs[ii]->SetXTitle(hh->GetXaxis()->GetTitle()); | |
361 | hhs[ii]->GetYaxis()->SetTitleOffset(hh->GetYaxis()->GetTitleOffset()); | |
362 | hhs[ii]->SetTitle(hh->GetTitle()); | |
363 | } | |
364 | } | |
365 | ||
366 | //=================================================================================== | |
367 | // TRD Analysis Fast Tool | |
368 | //=================================================================================== | |
369 | ||
370 | Int_t AliTRDdEdxBaseUtils::GetNtracklet(const AliESDEvent *esd) | |
371 | { | |
372 | // | |
373 | //number of trd tracklet in one esd event | |
374 | // | |
375 | const Int_t ntrk0 = esd->GetNumberOfTracks(); | |
376 | Int_t ntrdv1=0; | |
377 | for(Int_t ii=0; ii<ntrk0; ii++){ | |
378 | ntrdv1 += esd->GetTrack(ii)->GetTRDntracklets(); | |
379 | } | |
380 | return ntrdv1; | |
381 | } | |
382 | ||
383 | AliTRDtrackV1 * AliTRDdEdxBaseUtils::GetTRDtrackV1(const AliESDtrack * esdtrack) | |
384 | { | |
385 | // | |
386 | //Get TRD friend track | |
387 | // | |
388 | ||
389 | AliESDfriendTrack * friendtrk = (AliESDfriendTrack *)esdtrack->GetFriendTrack(); | |
390 | if(!friendtrk){ | |
391 | //printf("xlulog AliAnalysisTaskCosmicTRD::GetTRDtrack no friend!!\n"); exit(1); | |
392 | return 0x0; | |
393 | } | |
394 | ||
395 | TObject *calibObject=0x0; | |
396 | AliTRDtrackV1 * trdtrack=0x0; | |
397 | for(Int_t l=0; (calibObject=friendtrk->GetCalibObject(l)); l++) { | |
398 | if( (trdtrack=dynamic_cast<AliTRDtrackV1*>(calibObject)) ) | |
399 | break; | |
400 | } | |
401 | ||
402 | return trdtrack; | |
403 | } | |
404 | ||
405 | Bool_t AliTRDdEdxBaseUtils::IsInSameStack(const AliTRDtrackV1 *trdtrack) | |
406 | { | |
407 | // | |
408 | // to check if all tracklets are in the same stack, useful in cosmic | |
409 | // | |
410 | ||
411 | TVectorD secs(18), stks(5); | |
412 | ||
413 | for(Int_t ilayer = 0; ilayer < 6; ilayer++){ | |
414 | AliTRDseedV1 *tracklet=trdtrack->GetTracklet(ilayer); | |
415 | if(!tracklet) | |
416 | continue; | |
417 | ||
418 | const Int_t det = tracklet->GetDetector(); | |
419 | const Int_t isector = AliTRDgeometry::GetSector(det); | |
420 | const Int_t istack = AliTRDgeometry::GetStack(det); | |
421 | ||
422 | secs[isector] = 1; | |
423 | stks[istack] = 1; | |
424 | } | |
425 | ||
426 | if(secs.Sum()!=1 || stks.Sum()!=1){ | |
427 | return kFALSE; | |
428 | } | |
429 | else | |
430 | return kTRUE; | |
431 | } | |
432 | ||
ca66da8f | 433 | Double_t AliTRDdEdxBaseUtils::GetRedefinedPhi(Double_t phi) |
434 | { | |
435 | // | |
436 | //redefine phi in pi/2 - pi/2 + 2pi | |
437 | // | |
438 | ||
439 | if(phi>0 && phi < TMath::Pi()/2){ | |
440 | phi += 2*TMath::Pi(); | |
441 | } | |
442 | ||
443 | return phi; | |
444 | } | |
445 | ||
446 | Double_t AliTRDdEdxBaseUtils::Getdydx(const AliTRDseedV1 *tracklet) | |
fe829aa0 | 447 | { |
448 | // | |
ca66da8f | 449 | //get dydx |
fe829aa0 | 450 | // |
451 | if(tracklet) | |
ca66da8f | 452 | return tracklet->GetYref(1); |
fe829aa0 | 453 | else |
454 | return -999; | |
455 | } | |
456 | ||
ca66da8f | 457 | Double_t AliTRDdEdxBaseUtils::Getdzdx(const AliTRDseedV1 *tracklet) |
458 | { | |
459 | // | |
460 | //get dzdx | |
461 | // | |
462 | if(tracklet) | |
463 | return tracklet->GetZref(1); | |
464 | else | |
465 | return -999; | |
466 | } | |
467 | ||
468 | Double_t AliTRDdEdxBaseUtils::Getdldx(const AliTRDseedV1 *tracklet) | |
469 | { | |
470 | // | |
471 | //return angular normalization factor = dldx | |
472 | // | |
473 | if(tracklet) | |
474 | return TMath::Sqrt(1+tracklet->GetYref(1)*tracklet->GetYref(1)+tracklet->GetZref(1)*tracklet->GetZref(1)); | |
475 | else | |
476 | return -999; | |
477 | } | |
478 | ||
479 | AliTRDseedV1 * AliTRDdEdxBaseUtils::GetFirstTracklet(const AliTRDtrackV1 *trdtrack) | |
480 | { | |
481 | // | |
482 | //as function name | |
483 | // | |
484 | ||
485 | AliTRDseedV1 *tracklet = 0x0; | |
486 | for(Int_t ilayer = 0; ilayer < 6; ilayer++){ | |
487 | tracklet=trdtrack->GetTracklet(ilayer); | |
488 | if(!tracklet) | |
489 | continue; | |
490 | ||
491 | break; | |
492 | } | |
493 | ||
494 | return tracklet; | |
495 | } | |
496 | ||
fe829aa0 | 497 | AliTRDseedV1 * AliTRDdEdxBaseUtils::GetLastTracklet(const AliTRDtrackV1 *trdtrack) |
498 | { | |
499 | // | |
500 | //get last tracklet | |
501 | // | |
502 | AliTRDseedV1 *tracklet = 0x0; | |
503 | ||
504 | for(Int_t ilayer = 5; ilayer >= 0; ilayer--){ | |
505 | tracklet=trdtrack->GetTracklet(ilayer); | |
506 | if(!tracklet) | |
507 | continue; | |
508 | ||
509 | break; | |
510 | } | |
511 | ||
512 | return tracklet; | |
513 | } | |
514 | ||
ca66da8f | 515 | void AliTRDdEdxBaseUtils::GetFirstSectorStackMomentum(const AliTRDtrackV1 *trdtrack, Int_t & isec, Int_t & istk, Double_t & mom) |
6951a056 | 516 | { |
517 | // | |
518 | //as function name | |
519 | // | |
520 | isec = istk = -999; | |
521 | mom = -999; | |
522 | ||
ca66da8f | 523 | AliTRDseedV1 *tracklet = GetFirstTracklet(trdtrack); |
524 | if(tracklet){ | |
6951a056 | 525 | const Int_t det = tracklet->GetDetector(); |
526 | isec = AliTRDgeometry::GetSector(det); | |
527 | istk = AliTRDgeometry::GetStack(det); | |
ca66da8f | 528 | |
6951a056 | 529 | mom = tracklet->GetMomentum(); |
6951a056 | 530 | } |
531 | ||
ca66da8f | 532 | return; |
6951a056 | 533 | } |
534 | ||
535 | //=================================================================================== | |
536 | // Detector and Data Constant | |
537 | //=================================================================================== | |
538 | Int_t AliTRDdEdxBaseUtils::ToDetector(const Int_t gtb) | |
539 | { | |
540 | // | |
541 | //gtb = det*Ntb+itb | |
542 | // | |
543 | return gtb/AliTRDseedV1::kNtb; | |
544 | } | |
545 | ||
546 | Int_t AliTRDdEdxBaseUtils::ToTimeBin(const Int_t gtb) | |
547 | { | |
548 | // | |
549 | //gtb = det*Ntb+itb | |
550 | // | |
551 | return gtb%AliTRDseedV1::kNtb; | |
552 | } | |
553 | ||
554 | Int_t AliTRDdEdxBaseUtils::ToSector(const Int_t gtb) | |
555 | { | |
556 | // | |
557 | //return sector | |
558 | // | |
559 | return AliTRDgeometry::GetSector(ToDetector(gtb)); | |
560 | } | |
561 | ||
562 | Int_t AliTRDdEdxBaseUtils::ToStack(const Int_t gtb) | |
563 | { | |
564 | // | |
565 | //return stack | |
566 | // | |
567 | return AliTRDgeometry::GetStack(ToDetector(gtb)); | |
568 | } | |
569 | ||
570 | Int_t AliTRDdEdxBaseUtils::ToLayer(const Int_t gtb) | |
571 | { | |
572 | // | |
573 | //return layer | |
574 | // | |
575 | return AliTRDgeometry::GetLayer(ToDetector(gtb)); | |
576 | } | |
577 | ||
ca66da8f | 578 | void AliTRDdEdxBaseUtils::CheckRunB(const TString listrun1kg, const Int_t run, TString & type) |
579 | { | |
580 | // | |
581 | //check run b field | |
582 | // | |
583 | if(listrun1kg.Contains(Form("%d",run))){ | |
584 | type+="1kG"; | |
585 | } | |
586 | else { | |
587 | type+="5kG"; | |
588 | } | |
589 | } | |
590 | ||
6951a056 | 591 | TString AliTRDdEdxBaseUtils::GetRunType(const Int_t run) |
592 | { | |
593 | // | |
594 | //return run type | |
595 | // | |
596 | ||
597 | TString type; | |
598 | if(run>=121527 && run<= 126460)//LHC10d | |
ca66da8f | 599 | type="2010ppLHC10d"; |
6951a056 | 600 | else if(run>=126461 && run<= 130930)//LHC10e |
ca66da8f | 601 | type="2010ppLHC10e"; |
6951a056 | 602 | else if(run>=136782 && run <= 139846)//LHC10h |
ca66da8f | 603 | type="2010PbPbLHC10h"; |
6951a056 | 604 | else if(run>= 143224 && run<= 143237)//2011Feb |
ca66da8f | 605 | type="2011cosmicFeb"; |
6951a056 | 606 | else if(run>= 150587 && run<= 154930){ |
ca66da8f | 607 | type="2011cosmicMayJun"; |
6951a056 | 608 | |
ca66da8f | 609 | CheckRunB("154601 154602 154629 154634 154636 154639 154643", run, type); |
610 | } | |
611 | else if(run>=173047 && run<=180164){ | |
612 | type="2012cosmic"; | |
613 | ||
614 | CheckRunB("173047 173049 173050 173065 173070 173092 173095 173113 173131 173160 174154 174156 174192 174194", run, type); | |
6951a056 | 615 | } |
616 | else{ | |
617 | type="unknown"; | |
618 | } | |
619 | ||
620 | type.ToUpper(); | |
621 | return type; | |
622 | } | |
623 | ||
624 | void AliTRDdEdxBaseUtils::PrintControl() | |
625 | { | |
626 | // | |
627 | //print out control variable | |
628 | // | |
687aa844 | 629 | printf("\nAliTRDdEdxBaseUtils::PrintControl Q0Frac %.2f, Q1Frac %.2f, TimeBinCountCut %.2f, CalibTPCnclsCut %d, IsExBOn %d, IsPadGainOn %d, QScale %.2f\n\n", Q0Frac(), Q1Frac(), TimeBinCountCut(), CalibTPCnclsCut(), IsExBOn(), IsPadGainOn(), QScale()); |
6951a056 | 630 | } |
631 | ||
632 | //=================================================================================== | |
633 | // dEdx Parameterization | |
634 | //=================================================================================== | |
fe829aa0 | 635 | void AliTRDdEdxBaseUtils::FastFitdEdxTR(TH1 * hh) |
636 | { | |
637 | // | |
638 | //fit dedx tr from mean | |
639 | // | |
640 | ||
641 | TF1 *ff=new TF1("ff", MeandEdxTRLogx, -0.5, 4.5, 8); | |
642 | Double_t par[8]; | |
643 | par[0]= 2.397001e-01; | |
644 | par[1]= 1.334697e+00; | |
645 | par[2]= 6.967470e+00; | |
646 | par[3]= 9.055289e-02; | |
647 | par[4]= 9.388760e+00; | |
648 | par[5]= 9.452742e-04; | |
649 | par[6]= -1.866091e+00; | |
650 | par[7]= 1.403545e+00; | |
651 | ||
652 | ff->SetParameters(par); | |
653 | hh->Fit(ff,"N"); | |
654 | ||
655 | for(int ii=0; ii<8; ii++){ | |
656 | printf("par[%d]=%e;\n", ii, ff->GetParameter(ii)); | |
657 | } | |
658 | ||
659 | TFile *fout=new TFile("fastfit.root","recreate"); | |
660 | hh->Write(); | |
661 | ff->Write(); | |
662 | fout->Save(); | |
663 | fout->Close(); | |
664 | ||
665 | delete fout; | |
666 | delete ff; | |
667 | } | |
6951a056 | 668 | |
669 | Double_t AliTRDdEdxBaseUtils::Q0MeanTRDpp(const Double_t bg) | |
670 | { | |
671 | // | |
672 | //truncated Mean Q_{xx} in TRD | |
673 | // | |
674 | ||
675 | Double_t par[8]; | |
fe829aa0 | 676 | Double_t scale = 1; |
677 | ||
678 | //2012-0605 /u/xlu/.task/CommondEdx/myAnaData/newTune_r56595/inuse/plot/fastFit | |
679 | scale = 0.9257; | |
680 | par[0]=8.316476e-01; | |
681 | par[1]=1.600907e+00; | |
682 | par[2]=7.728447e+00; | |
683 | par[3]=6.235622e-02; | |
684 | par[4]=1.136209e+01; | |
685 | par[5]=-1.495764e-06; | |
686 | par[6]=-2.536119e+00; | |
687 | par[7]=3.416031e+00; | |
688 | ||
689 | return scale*MeandEdxTR(&bg, par); | |
6951a056 | 690 | } |
691 | ||
692 | Double_t AliTRDdEdxBaseUtils::Q0MeanTRDPbPb(const Double_t bg) | |
693 | { | |
694 | // | |
695 | //truncated Mean Q_{xx} in TRD | |
696 | // | |
697 | ||
698 | Double_t par[8]; | |
699 | ||
700 | //03132012161259 | |
701 | //opt: PbPbQ0 | |
702 | par[0]= 1.844912e-01; | |
703 | par[1]= 2.509702e+00; | |
704 | par[2]= 6.744031e+00; | |
705 | par[3]= 7.355123e-02; | |
706 | par[4]= 1.166023e+01; | |
707 | par[5]= 1.736186e-04; | |
708 | par[6]= -1.716063e+00; | |
709 | par[7]= 1.611366e+00; | |
710 | ||
711 | ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see4.log:hhtype4Q0b2c2 scale 0.460994 at ltbg 0.650000 | |
712 | return 0.460994*MeandEdxTR(&bg, par); | |
713 | } | |
714 | ||
715 | Double_t AliTRDdEdxBaseUtils::Q1MeanTRDpp(const Double_t bg) | |
716 | { | |
717 | // | |
718 | //truncated Mean 1/(1/Q)_{xx} in TRD | |
719 | // | |
720 | ||
721 | Double_t par[8]; | |
722 | ||
723 |