Line data Source code
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 :
44 : #include "AliESDEvent.h"
45 : #include "AliESDfriendTrack.h"
46 : #include "AliESDtrack.h"
47 : #include "AliTRDtrackV1.h"
48 :
49 : #include "AliTRDdEdxBaseUtils.h"
50 :
51 : #define EPSILON 1e-8
52 :
53 : Double_t AliTRDdEdxBaseUtils::fgQ0Frac = 0.55;
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;
59 : Double_t AliTRDdEdxBaseUtils::fgQScale = 50;
60 :
61 : //===================================================================================
62 : // Math and Histogram
63 : //===================================================================================
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 0 : const Int_t bins = axis->GetNbins();
71 :
72 0 : const Double_t from = axis->GetXmin();
73 0 : const Double_t to = axis->GetXmax();
74 0 : if (from<EPSILON){
75 : //printf("AliTRDdEdxBaseUtils::BinLogX warning xmin < epsilon! nothing done, axis not set. %e\n", from); exit(1);
76 0 : return;
77 : }
78 :
79 0 : Double_t *new_bins = new Double_t[bins + 1];
80 :
81 0 : new_bins[0] = from;
82 0 : const Double_t factor = pow(to/from, 1./bins);
83 :
84 0 : for (int i = 1; i <= bins; i++) {
85 0 : new_bins[i] = factor * new_bins[i-1];
86 : }
87 0 : axis->Set(bins, new_bins);
88 0 : delete [] new_bins;
89 0 : }
90 :
91 : void AliTRDdEdxBaseUtils::GetCDFCuts(const TH1D *hh, Int_t ncut, Double_t cuts[], const Double_t cdfs[], Double_t thres)
92 : {
93 : //
94 : //counts of hh is sorted
95 : //
96 :
97 0 : for(Int_t ii=0; ii<ncut; ii++){
98 0 : cuts[ii] = -999;
99 : }
100 :
101 : Int_t nsel = 0;
102 0 : const Int_t nbin = hh->GetNbinsX();
103 0 : Double_t datas[nbin];
104 0 : for(Int_t ii=1; ii<=nbin; ii++){
105 0 : const Double_t res = hh->GetBinContent(ii);
106 0 : if(res<thres){
107 0 : continue;
108 : }
109 :
110 0 : datas[nsel] = res;
111 0 : nsel++;
112 0 : }
113 0 : if(!nsel)
114 0 : return;
115 :
116 0 : Int_t id[nsel];
117 0 : TMath::Sort(nsel, datas, id, kFALSE);
118 :
119 0 : for(Int_t ii=0; ii<ncut; ii++){
120 0 : const Double_t icdf = cdfs[ii];
121 0 : if(icdf<0 || icdf>1){
122 0 : printf("AliTRDdEdxBaseUtils::GetCDFCuts error cdfs[%d] %15f out of range!\n", ii, icdf); exit(1);
123 : }
124 0 : cuts[ii] = datas[id[Int_t(icdf*nsel)]];
125 : }
126 0 : }
127 :
128 : Double_t AliTRDdEdxBaseUtils::GetMeanRMS(Double_t nn, Double_t sum, 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 208 : if(nn>EPSILON){
138 42 : tmean = sum/nn;
139 :
140 42 : const Double_t arg = w2s/nn-tmean*tmean;
141 42 : if(TMath::Abs(arg)<EPSILON){
142 : trms = 0;
143 0 : }
144 : else{
145 42 : if( arg <0 ){
146 0 : printf("AliTRDdEdxBaseUtils::GetMeanRMS error negative sqrt argument!! %e -- %e %e %f\n", arg, w2s, sum, nn); exit(1);
147 : }
148 :
149 42 : trms = TMath::Sqrt(arg);
150 : }
151 :
152 42 : terr = trms/TMath::Sqrt(nn);
153 42 : }
154 :
155 104 : if(grms){
156 0 : (*grms) = trms;
157 0 : }
158 :
159 104 : if(gerr){
160 0 : (*gerr) = terr;
161 0 : }
162 :
163 104 : return tmean;
164 : }
165 :
166 : Double_t AliTRDdEdxBaseUtils::TruncatedMean(Int_t nx, const Double_t xdata[], Double_t lowfrac, 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 208 : Int_t index[nx];
182 104 : TMath::Sort(nx, xdata, index, kFALSE);
183 :
184 : Int_t nused = 0;
185 : Double_t sum = 0;
186 : Double_t w2s = 0;
187 104 : const Int_t istart = Int_t (nx*lowfrac);
188 104 : const Int_t istop = Int_t (nx*highfrac);
189 :
190 : //=,< correct, because when low=0, high=1 it is correct
191 4994 : for(Int_t ii=istart; ii<istop; ii++){
192 : Double_t weight = 1;
193 2393 : if(wws){
194 0 : weight = wws[index[ii]];
195 0 : }
196 2393 : const Double_t sx = xdata[index[ii]]*weight;
197 :
198 2393 : sum += sx;
199 2393 : w2s += sx*sx;
200 :
201 2393 : nused++;
202 : //printf("test in loop %d/%d %f %f %f\n", ii, nused, sx, sum, w2s);
203 :
204 : }
205 :
206 104 : return GetMeanRMS(nused, sum, w2s, grms, gerr);
207 104 : }
208 :
209 : Double_t AliTRDdEdxBaseUtils::TruncatedMean(const TH1 *hh, Double_t lowfrac, 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 0 : for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
219 0 : const Double_t be = hh->GetBinError(itmp);
220 0 : const Double_t bc = hh->GetBinContent(itmp);
221 0 : if(be<EPSILON){
222 0 : if(bc>EPSILON){
223 0 : printf("AliTRDdEdxBaseUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
224 : }
225 0 : continue;
226 : }
227 0 : npreTrunc += bc*bc/be/be;
228 0 : }
229 :
230 0 : const Double_t nstart = npreTrunc*lowfrac;
231 0 : 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 0 : for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
239 0 : const Double_t be = hh->GetBinError(itmp);
240 0 : const Double_t bc = hh->GetBinContent(itmp);
241 0 : if(be<EPSILON){
242 0 : if(bc>EPSILON){
243 0 : printf("AliTRDdEdxBaseUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
244 : }
245 0 : continue;
246 : }
247 0 : const Double_t weight = bc*bc/be/be;
248 0 : ntot+=weight;
249 : //<= correct, because when high=1, nstop has to be included
250 0 : if(ntot>nstart && ntot<=nstop){
251 :
252 0 : const Double_t val = hh->GetBinCenter(itmp);
253 0 : sum += weight*val;
254 0 : w2s += weight*val*val;
255 :
256 0 : nused += weight;
257 :
258 : //printf("test %d %f %f --- %f %f -- %f %f\n", itmp, weight, val, sum, w2s, nused, nsample);
259 0 : }
260 0 : else if(ntot>nstop){
261 0 : if(itmp>=hh->GetNbinsX()){
262 0 : 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 0 : }
266 0 : break;
267 : }
268 0 : }
269 :
270 0 : return GetMeanRMS(nused, sum, w2s, grms, gerr);
271 : }
272 :
273 : void AliTRDdEdxBaseUtils::FitSlicesY(const TH2D *hh, TH1D *&hnor, TH1D *&hmpv, TH1D *&hwid, TH1D *&hres, Double_t thres, Double_t lowfrac, Double_t highfrac)
274 : {
275 : //
276 : //fit slices of hh using truncation
277 : //
278 :
279 0 : const Int_t x0 = hh->GetXaxis()->GetFirst();
280 0 : const Int_t x1 = hh->GetXaxis()->GetLast();
281 0 : const Int_t y0 = hh->GetYaxis()->GetFirst();
282 0 : const Int_t y1 = hh->GetYaxis()->GetLast();
283 :
284 0 : const Int_t nx = hh->GetNbinsX();
285 0 : const Int_t ny = hh->GetNbinsY();
286 0 : const Double_t xmin = hh->GetXaxis()->GetXmin();
287 0 : const Double_t xmax = hh->GetXaxis()->GetXmax();
288 0 : const Double_t ymin = hh->GetYaxis()->GetXmin();
289 0 : const Double_t ymax = hh->GetYaxis()->GetXmax();
290 :
291 0 : hnor = new TH1D(Form("%s_amp",hh->GetName()), "", nx, xmin, xmax);
292 0 : hmpv = new TH1D(Form("%s_mpv",hh->GetName()), "", nx, xmin, xmax);
293 0 : hwid = new TH1D(Form("%s_wid",hh->GetName()), "", nx, xmin, xmax);
294 0 : hres = new TH1D(Form("%s_res",hh->GetName()), "", nx, xmin, xmax);
295 :
296 0 : Double_t newbins[ny+1];
297 0 : for(Int_t ii=1; ii<=ny+1; ii++){
298 0 : const Double_t lowx= hh->GetYaxis()->GetBinLowEdge(ii);
299 0 : newbins[ii-1]=lowx;
300 : }
301 :
302 0 : for(Int_t ix=x0; ix<=x1; ix++){
303 : //to speed up
304 0 : const Double_t rawcount = hh->Integral(ix,ix,0, ny+1);
305 0 : if(rawcount<EPSILON){
306 0 : continue;
307 : }
308 :
309 0 : TH1D *htmp = new TH1D(Form("FitSlicesY_%s_%d", hh->GetName(), ix),"",ny, ymin, ymax);
310 0 : htmp->GetXaxis()->Set(ny, newbins);
311 : Double_t ntot = 0;
312 0 : for(Int_t iy=y0; iy<=y1; iy++){
313 0 : const Double_t be = hh->GetBinError(ix,iy);
314 0 : const Double_t bc = hh->GetBinContent(ix, iy);
315 :
316 0 : if(be<EPSILON){
317 0 : if(bc>EPSILON){
318 0 : printf("AliTRDdEdxBaseUtils::FitSlicesY error %d %d %e %e\n", ix, iy, be, bc); exit(1);
319 : }
320 0 : continue;
321 : }
322 :
323 0 : htmp->SetBinContent(iy, bc);
324 0 : htmp->SetBinError(iy, be);
325 :
326 0 : 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 0 : }
330 :
331 0 : hnor->SetBinContent(ix, ntot);
332 0 : hnor->SetBinError( ix, 0);
333 :
334 0 : if(ntot<thres || htmp->GetRMS()<EPSILON){
335 0 : delete htmp;
336 0 : continue;
337 : }
338 :
339 : //test htmp->Draw();
340 0 : Double_t trms = -999, terr = -999;
341 0 : const Double_t tmean = TruncatedMean(htmp, lowfrac, highfrac, &trms, &terr);
342 :
343 0 : hmpv->SetBinContent(ix, tmean);
344 0 : hmpv->SetBinError( ix, terr);
345 :
346 0 : hwid->SetBinContent(ix, trms);
347 0 : hwid->SetBinError( ix, 0);
348 :
349 0 : hres->SetBinContent(ix, tmean>EPSILON ? trms/tmean:0);
350 0 : hres->SetBinError( ix, 0);
351 :
352 0 : delete htmp;
353 0 : }
354 :
355 0 : TH1 *hhs[]={hnor, hmpv, hwid, hres};
356 0 : const TString yt[]={"N", "MPV", "#sigma", "#sigma/MPV"};
357 : const Int_t nh = sizeof(hhs)/sizeof(TH1*);
358 0 : for(Int_t ii=0; ii<nh; ii++){
359 0 : hhs[ii]->SetYTitle(Form("%s of %s", yt[ii].Data(), hh->GetYaxis()->GetTitle()));
360 0 : hhs[ii]->SetXTitle(hh->GetXaxis()->GetTitle());
361 0 : hhs[ii]->GetYaxis()->SetTitleOffset(hh->GetYaxis()->GetTitleOffset());
362 0 : hhs[ii]->SetTitle(hh->GetTitle());
363 : }
364 0 : }
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 0 : const Int_t ntrk0 = esd->GetNumberOfTracks();
376 : Int_t ntrdv1=0;
377 0 : for(Int_t ii=0; ii<ntrk0; ii++){
378 0 : ntrdv1 += esd->GetTrack(ii)->GetTRDntracklets();
379 : }
380 0 : return ntrdv1;
381 : }
382 :
383 : AliTRDtrackV1 * AliTRDdEdxBaseUtils::GetTRDtrackV1(const AliESDtrack * esdtrack)
384 : {
385 : //
386 : //Get TRD friend track
387 : //
388 :
389 0 : AliESDfriendTrack * friendtrk = (AliESDfriendTrack *)esdtrack->GetFriendTrack();
390 0 : if(!friendtrk){
391 : //printf("xlulog AliAnalysisTaskCosmicTRD::GetTRDtrack no friend!!\n"); exit(1);
392 0 : return 0x0;
393 : }
394 :
395 : TObject *calibObject=0x0;
396 : AliTRDtrackV1 * trdtrack=0x0;
397 0 : for(Int_t l=0; (calibObject=friendtrk->GetCalibObject(l)); l++) {
398 0 : if( (trdtrack=dynamic_cast<AliTRDtrackV1*>(calibObject)) )
399 0 : break;
400 : }
401 :
402 : return trdtrack;
403 0 : }
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 0 : TVectorD secs(18), stks(5);
412 :
413 0 : for(Int_t ilayer = 0; ilayer < 6; ilayer++){
414 0 : AliTRDseedV1 *tracklet=trdtrack->GetTracklet(ilayer);
415 0 : if(!tracklet)
416 0 : continue;
417 :
418 0 : const Int_t det = tracklet->GetDetector();
419 0 : const Int_t isector = AliTRDgeometry::GetSector(det);
420 0 : const Int_t istack = AliTRDgeometry::GetStack(det);
421 :
422 0 : secs[isector] = 1;
423 0 : stks[istack] = 1;
424 0 : }
425 :
426 0 : if(secs.Sum()!=1 || stks.Sum()!=1){
427 0 : return kFALSE;
428 : }
429 : else
430 0 : return kTRUE;
431 0 : }
432 :
433 : Double_t AliTRDdEdxBaseUtils::GetRedefinedPhi(Double_t phi)
434 : {
435 : //
436 : //redefine phi in pi/2 - pi/2 + 2pi
437 : //
438 :
439 0 : if(phi>0 && phi < TMath::Pi()/2){
440 0 : phi += 2*TMath::Pi();
441 0 : }
442 :
443 0 : return phi;
444 : }
445 :
446 : Double_t AliTRDdEdxBaseUtils::Getdydx(const AliTRDseedV1 *tracklet)
447 : {
448 : //
449 : //get dydx
450 : //
451 0 : if(tracklet)
452 0 : return tracklet->GetYref(1);
453 : else
454 0 : return -999;
455 0 : }
456 :
457 : Double_t AliTRDdEdxBaseUtils::Getdzdx(const AliTRDseedV1 *tracklet)
458 : {
459 : //
460 : //get dzdx
461 : //
462 0 : if(tracklet)
463 0 : return tracklet->GetZref(1);
464 : else
465 0 : return -999;
466 0 : }
467 :
468 : Double_t AliTRDdEdxBaseUtils::Getdldx(const AliTRDseedV1 *tracklet)
469 : {
470 : //
471 : //return angular normalization factor = dldx
472 : //
473 12772 : if(tracklet)
474 6386 : return TMath::Sqrt(1+tracklet->GetYref(1)*tracklet->GetYref(1)+tracklet->GetZref(1)*tracklet->GetZref(1));
475 : else
476 0 : return -999;
477 6386 : }
478 :
479 : AliTRDseedV1 * AliTRDdEdxBaseUtils::GetFirstTracklet(const AliTRDtrackV1 *trdtrack)
480 : {
481 : //
482 : //as function name
483 : //
484 :
485 : AliTRDseedV1 *tracklet = 0x0;
486 0 : for(Int_t ilayer = 0; ilayer < 6; ilayer++){
487 0 : tracklet=trdtrack->GetTracklet(ilayer);
488 0 : if(!tracklet)
489 : continue;
490 :
491 0 : break;
492 : }
493 :
494 0 : return tracklet;
495 : }
496 :
497 : AliTRDseedV1 * AliTRDdEdxBaseUtils::GetLastTracklet(const AliTRDtrackV1 *trdtrack)
498 : {
499 : //
500 : //get last tracklet
501 : //
502 : AliTRDseedV1 *tracklet = 0x0;
503 :
504 0 : for(Int_t ilayer = 5; ilayer >= 0; ilayer--){
505 0 : tracklet=trdtrack->GetTracklet(ilayer);
506 0 : if(!tracklet)
507 : continue;
508 :
509 0 : break;
510 : }
511 :
512 0 : return tracklet;
513 : }
514 :
515 : void AliTRDdEdxBaseUtils::GetFirstSectorStackMomentum(const AliTRDtrackV1 *trdtrack, Int_t & isec, Int_t & istk, Double_t & mom)
516 : {
517 : //
518 : //as function name
519 : //
520 0 : isec = istk = -999;
521 0 : mom = -999;
522 :
523 0 : AliTRDseedV1 *tracklet = GetFirstTracklet(trdtrack);
524 0 : if(tracklet){
525 0 : const Int_t det = tracklet->GetDetector();
526 0 : isec = AliTRDgeometry::GetSector(det);
527 0 : istk = AliTRDgeometry::GetStack(det);
528 :
529 0 : mom = tracklet->GetMomentum();
530 0 : }
531 :
532 : return;
533 0 : }
534 :
535 : //===================================================================================
536 : // Detector and Data Constant
537 : //===================================================================================
538 : Int_t AliTRDdEdxBaseUtils::ToDetector(Int_t gtb)
539 : {
540 : //
541 : //gtb = det*Ntb+itb
542 : //
543 8888 : return gtb/AliTRDseedV1::kNtb;
544 : }
545 :
546 : Int_t AliTRDdEdxBaseUtils::ToTimeBin(Int_t gtb)
547 : {
548 : //
549 : //gtb = det*Ntb+itb
550 : //
551 8888 : return gtb%AliTRDseedV1::kNtb;
552 : }
553 :
554 : Int_t AliTRDdEdxBaseUtils::ToSector(Int_t gtb)
555 : {
556 : //
557 : //return sector
558 : //
559 0 : return AliTRDgeometry::GetSector(ToDetector(gtb));
560 : }
561 :
562 : Int_t AliTRDdEdxBaseUtils::ToStack(Int_t gtb)
563 : {
564 : //
565 : //return stack
566 : //
567 0 : return AliTRDgeometry::GetStack(ToDetector(gtb));
568 : }
569 :
570 : Int_t AliTRDdEdxBaseUtils::ToLayer(Int_t gtb)
571 : {
572 : //
573 : //return layer
574 : //
575 0 : return AliTRDgeometry::GetLayer(ToDetector(gtb));
576 : }
577 :
578 : void AliTRDdEdxBaseUtils::CheckRunB(TString listrun1kg, Int_t run, TString & type)
579 : {
580 : //
581 : //check run b field
582 : //
583 0 : if(listrun1kg.Contains(Form("%d",run))){
584 0 : type+="1kG";
585 0 : }
586 : else {
587 0 : type+="5kG";
588 : }
589 0 : }
590 :
591 : TString AliTRDdEdxBaseUtils::GetRunType(Int_t run)
592 : {
593 : //
594 : //return run type
595 : //
596 :
597 0 : TString type;
598 0 : if(run>=121527 && run<= 126460)//LHC10d
599 0 : type="2010ppLHC10d";
600 0 : else if(run>=126461 && run<= 130930)//LHC10e
601 0 : type="2010ppLHC10e";
602 0 : else if(run>=136782 && run <= 139846)//LHC10h
603 0 : type="2010PbPbLHC10h";
604 0 : else if(run>= 143224 && run<= 143237)//2011Feb
605 0 : type="2011cosmicFeb";
606 0 : else if(run>= 150587 && run<= 154930){
607 0 : type="2011cosmicMayJun";
608 :
609 0 : CheckRunB("154601 154602 154629 154634 154636 154639 154643", run, type);
610 0 : }
611 0 : else if(run>=173047 && run<=180164){
612 0 : type="2012cosmic";
613 :
614 0 : CheckRunB("173047 173049 173050 173065 173070 173092 173095 173113 173131 173160 174154 174156 174192 174194", run, type);
615 0 : }
616 : else{
617 0 : type="unknown";
618 : }
619 :
620 0 : type.ToUpper();
621 : return type;
622 0 : }
623 :
624 : void AliTRDdEdxBaseUtils::PrintControl()
625 : {
626 : //
627 : //print out control variable
628 : //
629 4 : 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());
630 2 : }
631 :
632 : //===================================================================================
633 : // dEdx Parameterization
634 : //===================================================================================
635 : void AliTRDdEdxBaseUtils::FastFitdEdxTR(TH1 * hh)
636 : {
637 : //
638 : //fit dedx tr from mean
639 : //
640 :
641 0 : TF1 *ff=new TF1("ff", MeandEdxTRLogx, -0.5, 4.5, 8);
642 0 : Double_t par[8];
643 0 : par[0]= 2.397001e-01;
644 0 : par[1]= 1.334697e+00;
645 0 : par[2]= 6.967470e+00;
646 0 : par[3]= 9.055289e-02;
647 0 : par[4]= 9.388760e+00;
648 0 : par[5]= 9.452742e-04;
649 0 : par[6]= -1.866091e+00;
650 0 : par[7]= 1.403545e+00;
651 :
652 0 : ff->SetParameters(par);
653 0 : hh->Fit(ff,"N");
654 :
655 0 : for(int ii=0; ii<8; ii++){
656 0 : printf("par[%d]=%e;\n", ii, ff->GetParameter(ii));
657 : }
658 :
659 0 : TFile *fout=new TFile("fastfit.root","recreate");
660 0 : hh->Write();
661 0 : ff->Write();
662 0 : fout->Save();
663 0 : fout->Close();
664 :
665 0 : delete fout;
666 0 : delete ff;
667 0 : }
668 :
669 : Double_t AliTRDdEdxBaseUtils::Q0MeanTRDpp(Double_t bg)
670 : {
671 : //
672 : //truncated Mean Q_{xx} in TRD
673 : //
674 :
675 0 : Double_t par[8];
676 : Double_t scale = 1;
677 :
678 : //2012-0605 /u/xlu/.task/CommondEdx/myAnaData/newTune_r56595/inuse/plot/fastFit
679 : scale = 0.9257;
680 0 : par[0]=8.316476e-01;
681 0 : par[1]=1.600907e+00;
682 0 : par[2]=7.728447e+00;
683 0 : par[3]=6.235622e-02;
684 0 : par[4]=1.136209e+01;
685 0 : par[5]=-1.495764e-06;
686 0 : par[6]=-2.536119e+00;
687 0 : par[7]=3.416031e+00;
688 :
689 0 : return scale*MeandEdxTR(&bg, par);
690 0 : }
691 :
692 : Double_t AliTRDdEdxBaseUtils::Q0MeanTRDPbPb(Double_t bg)
693 : {
694 : //
695 : //truncated Mean Q_{xx} in TRD
696 : //
697 :
698 0 : Double_t par[8];
699 :
700 : //03132012161259
701 : //opt: PbPbQ0
702 0 : par[0]= 1.844912e-01;
703 0 : par[1]= 2.509702e+00;
704 0 : par[2]= 6.744031e+00;
705 0 : par[3]= 7.355123e-02;
706 0 : par[4]= 1.166023e+01;
707 0 : par[5]= 1.736186e-04;
708 0 : par[6]= -1.716063e+00;
709 0 : 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 0 : return 0.460994*MeandEdxTR(&bg, par);
713 0 : }
714 :
715 : Double_t AliTRDdEdxBaseUtils::Q1MeanTRDpp(Double_t bg)
716 : {
717 : //
718 : //truncated Mean 1/(1/Q)_{xx} in TRD
719 : //
720 :
721 0 : Double_t par[8];
722 :
723 : //So 4. Mär 13:30:51 CET 2012
724 : //opt: trdppQ1
725 0 : par[0]= 2.434646e-01;
726 0 : par[1]= 1.400211e+00;
727 0 : par[2]= 6.937471e+00;
728 0 : par[3]= 7.758118e-02;
729 0 : par[4]= 1.097372e+01;
730 0 : par[5]= 4.297518e-04;
731 0 : par[6]= -1.806266e+00;
732 0 : par[7]= 1.543811e+00;
733 :
734 : //hhtype2Q1b2c2 scale 0.418629 at ltbg 0.650000
735 :
736 0 : return 0.418629*MeandEdxTR(&bg, par);
737 0 : }
738 :
739 : Double_t AliTRDdEdxBaseUtils::Q1MeanTRDPbPb(Double_t bg)
740 : {
741 : //
742 : //truncated Mean 1/(1/Q)_{xx} in TRD
743 : //
744 :
745 0 : Double_t par[8];
746 :
747 : //So 4. Mär 13:30:52 CET 2012
748 : //opt: trdPbPbQ1
749 0 : par[0]= 2.193660e-01;
750 0 : par[1]= 2.051864e+00;
751 0 : par[2]= 6.825112e+00;
752 0 : par[3]= 6.151693e-02;
753 0 : par[4]= 1.390343e+01;
754 0 : par[5]= 6.010032e-05;
755 0 : par[6]= -1.676324e+00;
756 0 : par[7]= 1.838873e+00;
757 :
758 : //hhtype4Q1b2c2 scale 0.457988 at ltbg 0.650000
759 :
760 0 : return 0.457988*MeandEdxTR(&bg, par);
761 0 : }
762 :
763 : Double_t AliTRDdEdxBaseUtils::QMeanTPC(Double_t bg)
764 : {
765 : //
766 : //bethe bloch in TPC
767 : //
768 :
769 0 : Double_t par[5];
770 : //Mi 15. Feb 14:48:05 CET 2012
771 : //train_2012-02-13_1214.12001, tpcsig
772 0 : par[0]= 4.401269;
773 0 : par[1]= 9.725370;
774 0 : par[2]= 0.000178;
775 0 : par[3]= 1.904962;
776 0 : par[4]= 1.426576;
777 :
778 0 : return MeandEdx(&bg, par);
779 0 : }
780 :
781 : Double_t AliTRDdEdxBaseUtils::MeandEdxTR(const Double_t * xx, const Double_t * pin)
782 : {
783 : //
784 : //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
785 : //npar = 8 = 3+5
786 : //
787 0 : Double_t ptr[4]={0};
788 0 : for(int ii=0; ii<3; ii++){
789 0 : ptr[ii+1]=pin[ii];
790 : }
791 0 : return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
792 0 : }
793 :
794 : Double_t AliTRDdEdxBaseUtils::MeanTR(const Double_t * xx, const Double_t * par)
795 : {
796 : //
797 : //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
798 : //npar = 4
799 : //
800 :
801 0 : const Double_t bg = xx[0];
802 0 : const Double_t gamma = sqrt(1+bg*bg);
803 :
804 0 : const Double_t p0 = TMath::Abs(par[1]);
805 0 : const Double_t p1 = TMath::Abs(par[2]);
806 0 : const Double_t p2 = TMath::Abs(par[3]);
807 :
808 0 : const Double_t zz = TMath::Log(gamma);
809 0 : const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
810 :
811 0 : return par[0]+tryield;
812 : }
813 :
814 : Double_t AliTRDdEdxBaseUtils::MeandEdx(const Double_t * xx, const Double_t * par)
815 : {
816 0 : return ALEPH(xx, par);
817 : }
818 :
819 : Double_t AliTRDdEdxBaseUtils::ALEPH(const Double_t * xx, const Double_t * par)
820 : {
821 : //
822 : //ALEPH parametrization for dEdx
823 : //npar = 5
824 : //
825 :
826 0 : const Double_t bg = xx[0];
827 0 : const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
828 :
829 0 : const Double_t p0 = TMath::Abs(par[0]);
830 0 : const Double_t p1 = TMath::Abs(par[1]);
831 :
832 : //after redefining p2 (<0) -> ln P2
833 : //const Double_t p2 = par[2];
834 :
835 : //restore back
836 0 : const Double_t p2 = TMath::Abs(par[2]);
837 :
838 0 : const Double_t p3 = TMath::Abs(par[3]);
839 0 : const Double_t p4 = TMath::Abs(par[4]);
840 :
841 0 : const Double_t aa = TMath::Power(beta, p3);
842 :
843 : //numerically very bad when p2~1e-15, bg~1e3, p4~5.
844 0 : const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
845 :
846 : //+1 to avoid the singularity when bg<1 and p4>>1
847 : //very^inf important! with this hesse NOTPOS->OK and no nan or inf in migrad
848 : //i.e. numerically very stable
849 : //the reason is when bg<1 ln blows up very fast with p4>>1 and nan/inf appears in migrad search
850 : //the fitting will adjust the parameters as if bg is not shifted, the fitted curve looks fine!!
851 : //-- 2012 Aug 8
852 : //----> fail for 10h, not used, restore back!
853 : //const Double_t lbg = TMath::Log(bg);
854 :
855 : //redefine p2(<0) -> ln P2
856 : //corresponds to alephFit.C fAlephOPt = 11
857 : //const Double_t bb = p2 + TMath::Log( 1 + TMath::Exp(-p2-p4*lbg) );
858 :
859 : //printf("test----- %f %f -- %f %f %f %f %f --- %f %f %f\n", bg, beta, p0, p1, p2, p3, p4, p0/aa, aa, bb);
860 :
861 0 : return (p1-aa-bb)*p0/aa;
862 : }
863 :
864 : Double_t AliTRDdEdxBaseUtils::ToLogx(FFunc func, const Double_t * xx, const Double_t * par)
865 : {
866 : //
867 : //f(x)-> f(y) with y=log10(x)
868 : //
869 0 : const Double_t x2[]={TMath::Power(10, xx[0])};
870 0 : return func(x2, par);
871 0 : }
872 :
|
