1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 // TRD dEdx base utils
24 // lu@physi.uni-heidelberg.de
33 #include "THnSparse.h"
37 #include "TObjArray.h"
39 #include "TStopwatch.h"
42 #include "TTreeStream.h"
44 #include "AliESDEvent.h"
45 #include "AliESDfriendTrack.h"
46 #include "AliESDtrack.h"
47 #include "AliTRDtrackV1.h"
49 #include "AliTRDdEdxBaseUtils.h"
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;
61 //===================================================================================
63 //===================================================================================
64 void AliTRDdEdxBaseUtils::BinLogX(TAxis *axis)
67 // Method for the correct logarithmic binning of histograms
68 // copied and modified from AliTPCcalibBase
70 const Int_t bins = axis->GetNbins();
72 const Double_t from = axis->GetXmin();
73 const Double_t to = axis->GetXmax();
74 if (from<EPSILON) return;
75 Double_t *new_bins = new Double_t[bins + 1];
78 const Double_t factor = pow(to/from, 1./bins);
80 for (int i = 1; i <= bins; i++) {
81 new_bins[i] = factor * new_bins[i-1];
83 axis->Set(bins, new_bins);
87 void AliTRDdEdxBaseUtils::GetCDFCuts(const TH1D *hh, const Int_t ncut, Double_t cuts[], const Double_t cdfs[], const Double_t thres)
90 //counts of hh is sorted
93 for(Int_t ii=0; ii<ncut; ii++){
98 const Int_t nbin = hh->GetNbinsX();
100 for(Int_t ii=1; ii<=nbin; ii++){
101 const Double_t res = hh->GetBinContent(ii);
113 TMath::Sort(nsel, datas, id, kFALSE);
115 for(Int_t ii=0; ii<ncut; ii++){
116 const Double_t icdf = cdfs[ii];
117 if(icdf<0 || icdf>1){
118 printf("AliTRDdEdxBaseUtils::GetCDFCuts error cdfs[%d] %15f out of range!\n", ii, icdf); exit(1);
120 cuts[ii] = datas[id[Int_t(icdf*nsel)]];
124 Double_t AliTRDdEdxBaseUtils::GetMeanRMS(const Double_t nn, const Double_t sum, const Double_t w2s, Double_t * grms, Double_t * gerr)
127 //calculate mean (with error) and rms from sum, w2s, nn
128 //if nn=0, mean, error, and rms all = 0
131 Double_t tmean = 0, trms = 0, terr = 0;
136 const Double_t arg = w2s/nn-tmean*tmean;
137 if(TMath::Abs(arg)<EPSILON){
142 printf("AliTRDdEdxBaseUtils::GetMeanRMS error negative sqrt argument!! %e -- %e %e %f\n", arg, w2s, sum, nn); exit(1);
145 trms = TMath::Sqrt(arg);
148 terr = trms/TMath::Sqrt(nn);
162 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)
165 //calculate truncated mean
166 //return <x*w>_{low-high according to x}
171 for(Int_t ii=0; ii<nx; ii++){
172 printf("test %d/%d %f\n", ii, nx, xdata[ii]);
178 TMath::Sort(nx, xdata, index, kFALSE);
183 const Int_t istart = Int_t (nx*lowfrac);
184 const Int_t istop = Int_t (nx*highfrac);
186 //=,< correct, because when low=0, high=1 it is correct
187 for(Int_t ii=istart; ii<istop; ii++){
190 weight = wws[index[ii]];
192 const Double_t sx = xdata[index[ii]]*weight;
198 //printf("test in loop %d/%d %f %f %f\n", ii, nused, sx, sum, w2s);
202 return GetMeanRMS(nused, sum, w2s, grms, gerr);
205 Double_t AliTRDdEdxBaseUtils::TruncatedMean(const TH1 *hh, const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr)
208 //do truncation on histogram
210 //if hh is scaled, be sure Sumw2 is called before scaling!! then mean, rms and err will all be correct
212 //with under-/over-flow
213 Double_t npreTrunc = 0;
214 for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
215 const Double_t be = hh->GetBinError(itmp);
216 const Double_t bc = hh->GetBinContent(itmp);
219 printf("AliTRDdEdxBaseUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
223 npreTrunc += bc*bc/be/be;
226 const Double_t nstart = npreTrunc*lowfrac;
227 const Double_t nstop = npreTrunc*highfrac;
229 //with Double_t this should also handle normalized hist
234 for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
235 const Double_t be = hh->GetBinError(itmp);
236 const Double_t bc = hh->GetBinContent(itmp);
239 printf("AliTRDdEdxBaseUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
243 const Double_t weight = bc*bc/be/be;
245 //<= correct, because when high=1, nstop has to be included
246 if(ntot>nstart && ntot<=nstop){
248 const Double_t val = hh->GetBinCenter(itmp);
250 w2s += weight*val*val;
254 //printf("test %d %f %f --- %f %f -- %f %f\n", itmp, weight, val, sum, w2s, nused, nsample);
257 if(itmp>=hh->GetNbinsX()){
258 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);
266 return GetMeanRMS(nused, sum, w2s, grms, gerr);
269 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)
272 //fit slices of hh using truncation
275 const Int_t x0 = hh->GetXaxis()->GetFirst();
276 const Int_t x1 = hh->GetXaxis()->GetLast();
277 const Int_t y0 = hh->GetYaxis()->GetFirst();
278 const Int_t y1 = hh->GetYaxis()->GetLast();
280 const Int_t nx = hh->GetNbinsX();
281 const Int_t ny = hh->GetNbinsY();
282 const Double_t xmin = hh->GetXaxis()->GetXmin();
283 const Double_t xmax = hh->GetXaxis()->GetXmax();
284 const Double_t ymin = hh->GetYaxis()->GetXmin();
285 const Double_t ymax = hh->GetYaxis()->GetXmax();
287 hnor = new TH1D(Form("%s_amp",hh->GetName()), "", nx, xmin, xmax);
288 hmpv = new TH1D(Form("%s_mpv",hh->GetName()), "", nx, xmin, xmax);
289 hwid = new TH1D(Form("%s_wid",hh->GetName()), "", nx, xmin, xmax);
290 hres = new TH1D(Form("%s_res",hh->GetName()), "", nx, xmin, xmax);
292 Double_t newbins[ny+1];
293 for(Int_t ii=1; ii<=ny+1; ii++){
294 const Double_t lowx= hh->GetYaxis()->GetBinLowEdge(ii);
298 for(Int_t ix=x0; ix<=x1; ix++){
300 const Double_t rawcount = hh->Integral(ix,ix,0, ny+1);
301 if(rawcount<EPSILON){
305 TH1D *htmp = new TH1D(Form("FitSlicesY_%s_%d", hh->GetName(), ix),"",ny, ymin, ymax);
306 htmp->GetXaxis()->Set(ny, newbins);
308 for(Int_t iy=y0; iy<=y1; iy++){
309 const Double_t be = hh->GetBinError(ix,iy);
310 const Double_t bc = hh->GetBinContent(ix, iy);
314 printf("AliTRDdEdxBaseUtils::FitSlicesY error %d %d %e %e\n", ix, iy, be, bc); exit(1);
319 htmp->SetBinContent(iy, bc);
320 htmp->SetBinError(iy, be);
322 ntot += (bc/be)*(bc/be);
324 //if(be) printf("test %d %d : %f %f %f\n", ix, iy, bc, be, pow(bc/be,2));
327 hnor->SetBinContent(ix, ntot);
328 hnor->SetBinError( ix, 0);
330 if(ntot<thres || htmp->GetRMS()<EPSILON){
336 Double_t trms = -999, terr = -999;
337 const Double_t tmean = TruncatedMean(htmp, lowfrac, highfrac, &trms, &terr);
339 hmpv->SetBinContent(ix, tmean);
340 hmpv->SetBinError( ix, terr);
342 hwid->SetBinContent(ix, trms);
343 hwid->SetBinError( ix, 0);
345 hres->SetBinContent(ix, tmean>EPSILON ? trms/tmean:0);
346 hres->SetBinError( ix, 0);
351 TH1 *hhs[]={hnor, hmpv, hwid, hres};
352 const TString yt[]={"N", "MPV", "#sigma", "#sigma/MPV"};
353 const Int_t nh = sizeof(hhs)/sizeof(TH1*);
354 for(Int_t ii=0; ii<nh; ii++){
355 hhs[ii]->SetYTitle(Form("%s of %s", yt[ii].Data(), hh->GetYaxis()->GetTitle()));
356 hhs[ii]->SetXTitle(hh->GetXaxis()->GetTitle());
357 hhs[ii]->GetYaxis()->SetTitleOffset(hh->GetYaxis()->GetTitleOffset());
358 hhs[ii]->SetTitle(hh->GetTitle());
362 //===================================================================================
363 // TRD Analysis Fast Tool
364 //===================================================================================
366 Int_t AliTRDdEdxBaseUtils::GetNtracklet(const AliESDEvent *esd)
369 //number of trd tracklet in one esd event
371 const Int_t ntrk0 = esd->GetNumberOfTracks();
373 for(Int_t ii=0; ii<ntrk0; ii++){
374 ntrdv1 += esd->GetTrack(ii)->GetTRDntracklets();
379 AliTRDtrackV1 * AliTRDdEdxBaseUtils::GetTRDtrackV1(const AliESDtrack * esdtrack)
382 //Get TRD friend track
385 AliESDfriendTrack * friendtrk = (AliESDfriendTrack *)esdtrack->GetFriendTrack();
387 //printf("xlulog AliAnalysisTaskCosmicTRD::GetTRDtrack no friend!!\n"); exit(1);
391 TObject *calibObject=0x0;
392 AliTRDtrackV1 * trdtrack=0x0;
393 for(Int_t l=0; (calibObject=friendtrk->GetCalibObject(l)); l++) {
394 if( (trdtrack=dynamic_cast<AliTRDtrackV1*>(calibObject)) )
401 Bool_t AliTRDdEdxBaseUtils::IsInSameStack(const AliTRDtrackV1 *trdtrack)
404 // to check if all tracklets are in the same stack, useful in cosmic
407 TVectorD secs(18), stks(5);
409 for(Int_t ilayer = 0; ilayer < 6; ilayer++){
410 AliTRDseedV1 *tracklet=trdtrack->GetTracklet(ilayer);
414 const Int_t det = tracklet->GetDetector();
415 const Int_t isector = AliTRDgeometry::GetSector(det);
416 const Int_t istack = AliTRDgeometry::GetStack(det);
422 if(secs.Sum()!=1 || stks.Sum()!=1){
429 Double_t AliTRDdEdxBaseUtils::GetDeltaPhi(const AliTRDseedV1 *tracklet)
432 //get phi difference to normal incidence
435 return TMath::ATan(tracklet->GetYref(1));
440 AliTRDseedV1 * AliTRDdEdxBaseUtils::GetLastTracklet(const AliTRDtrackV1 *trdtrack)
445 AliTRDseedV1 *tracklet = 0x0;
447 for(Int_t ilayer = 5; ilayer >= 0; ilayer--){
448 tracklet=trdtrack->GetTracklet(ilayer);
458 AliTRDseedV1 * AliTRDdEdxBaseUtils::GetFirstSectorStackMomentum(const AliTRDtrackV1 *trdtrack, Int_t & isec, Int_t & istk, Double_t & mom)
466 AliTRDseedV1 *tracklet = 0x0;
467 for(Int_t ilayer = 0; ilayer < 6; ilayer++){
468 tracklet=trdtrack->GetTracklet(ilayer);
472 const Int_t det = tracklet->GetDetector();
473 isec = AliTRDgeometry::GetSector(det);
474 istk = AliTRDgeometry::GetStack(det);
476 mom = tracklet->GetMomentum();
484 //===================================================================================
485 // Detector and Data Constant
486 //===================================================================================
487 Int_t AliTRDdEdxBaseUtils::ToDetector(const Int_t gtb)
492 return gtb/AliTRDseedV1::kNtb;
495 Int_t AliTRDdEdxBaseUtils::ToTimeBin(const Int_t gtb)
500 return gtb%AliTRDseedV1::kNtb;
503 Int_t AliTRDdEdxBaseUtils::ToSector(const Int_t gtb)
508 return AliTRDgeometry::GetSector(ToDetector(gtb));
511 Int_t AliTRDdEdxBaseUtils::ToStack(const Int_t gtb)
516 return AliTRDgeometry::GetStack(ToDetector(gtb));
519 Int_t AliTRDdEdxBaseUtils::ToLayer(const Int_t gtb)
524 return AliTRDgeometry::GetLayer(ToDetector(gtb));
527 TString AliTRDdEdxBaseUtils::GetRunType(const Int_t run)
534 if(run>=121527 && run<= 126460)//LHC10d
536 else if(run>=126461 && run<= 130930)//LHC10e
538 else if(run>=136782 && run <= 139846)//LHC10h
539 type="PbPb2010LHC10h";
540 else if(run>= 143224 && run<= 143237)//2011Feb
541 type="cosmic2011Feb";
542 else if(run>= 150587 && run<= 154930){
543 type="cosmic2011MayJun";
545 TString runstr(Form("%d",run));
546 const TString listrun1kg("154601 154602 154629 154634 154636 154639 154643");
547 if(listrun1kg.Contains(runstr)){
562 void AliTRDdEdxBaseUtils::PrintControl()
565 //print out control variable
567 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());
570 //===================================================================================
571 // dEdx Parameterization
572 //===================================================================================
573 void AliTRDdEdxBaseUtils::FastFitdEdxTR(TH1 * hh)
576 //fit dedx tr from mean
579 TF1 *ff=new TF1("ff", MeandEdxTRLogx, -0.5, 4.5, 8);
581 par[0]= 2.397001e-01;
582 par[1]= 1.334697e+00;
583 par[2]= 6.967470e+00;
584 par[3]= 9.055289e-02;
585 par[4]= 9.388760e+00;
586 par[5]= 9.452742e-04;
587 par[6]= -1.866091e+00;
588 par[7]= 1.403545e+00;
590 ff->SetParameters(par);
593 for(int ii=0; ii<8; ii++){
594 printf("par[%d]=%e;\n", ii, ff->GetParameter(ii));
597 TFile *fout=new TFile("fastfit.root","recreate");
607 Double_t AliTRDdEdxBaseUtils::Q0MeanTRDpp(const Double_t bg)
610 //truncated Mean Q_{xx} in TRD
616 //2012-0605 /u/xlu/.task/CommondEdx/myAnaData/newTune_r56595/inuse/plot/fastFit
623 par[5]=-1.495764e-06;
624 par[6]=-2.536119e+00;
627 return scale*MeandEdxTR(&bg, par);
630 Double_t AliTRDdEdxBaseUtils::Q0MeanTRDPbPb(const Double_t bg)
633 //truncated Mean Q_{xx} in TRD
640 par[0]= 1.844912e-01;
641 par[1]= 2.509702e+00;
642 par[2]= 6.744031e+00;
643 par[3]= 7.355123e-02;
644 par[4]= 1.166023e+01;
645 par[5]= 1.736186e-04;
646 par[6]= -1.716063e+00;
647 par[7]= 1.611366e+00;
649 ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see4.log:hhtype4Q0b2c2 scale 0.460994 at ltbg 0.650000
650 return 0.460994*MeandEdxTR(&bg, par);
653 Double_t AliTRDdEdxBaseUtils::Q1MeanTRDpp(const Double_t bg)
656 //truncated Mean 1/(1/Q)_{xx} in TRD
661 //So 4. Mär 13:30:51 CET 2012
663 par[0]= 2.434646e-01;
664 par[1]= 1.400211e+00;
665 par[2]= 6.937471e+00;
666 par[3]= 7.758118e-02;
667 par[4]= 1.097372e+01;
668 par[5]= 4.297518e-04;
669 par[6]= -1.806266e+00;
670 par[7]= 1.543811e+00;
672 //hhtype2Q1b2c2 scale 0.418629 at ltbg 0.650000
674 return 0.418629*MeandEdxTR(&bg, par);
677 Double_t AliTRDdEdxBaseUtils::Q1MeanTRDPbPb(const Double_t bg)
680 //truncated Mean 1/(1/Q)_{xx} in TRD
685 //So 4. Mär 13:30:52 CET 2012
687 par[0]= 2.193660e-01;
688 par[1]= 2.051864e+00;
689 par[2]= 6.825112e+00;
690 par[3]= 6.151693e-02;
691 par[4]= 1.390343e+01;
692 par[5]= 6.010032e-05;
693 par[6]= -1.676324e+00;
694 par[7]= 1.838873e+00;
696 //hhtype4Q1b2c2 scale 0.457988 at ltbg 0.650000
698 return 0.457988*MeandEdxTR(&bg, par);
701 Double_t AliTRDdEdxBaseUtils::QMeanTPC(const Double_t bg)
708 //Mi 15. Feb 14:48:05 CET 2012
709 //train_2012-02-13_1214.12001, tpcsig
716 return MeandEdx(&bg, par);
719 Double_t AliTRDdEdxBaseUtils::MeandEdxTR(const Double_t * xx, const Double_t * pin)
722 //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
726 for(int ii=0; ii<3; ii++){
729 return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
732 Double_t AliTRDdEdxBaseUtils::MeanTR(const Double_t * xx, const Double_t * par)
735 //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
739 const Double_t bg = xx[0];
740 const Double_t gamma = sqrt(1+bg*bg);
742 const Double_t p0 = TMath::Abs(par[1]);
743 const Double_t p1 = TMath::Abs(par[2]);
744 const Double_t p2 = TMath::Abs(par[3]);
746 const Double_t zz = TMath::Log(gamma);
747 const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
749 return par[0]+tryield;
752 Double_t AliTRDdEdxBaseUtils::MeandEdx(const Double_t * xx, const Double_t * par)
755 //ALEPH parametrization for dEdx
759 const Double_t bg = xx[0];
760 const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
762 const Double_t p0 = TMath::Abs(par[0]);
763 const Double_t p1 = TMath::Abs(par[1]);
764 const Double_t p2 = TMath::Abs(par[2]);
765 const Double_t p3 = TMath::Abs(par[3]);
766 const Double_t p4 = TMath::Abs(par[4]);
768 const Double_t aa = TMath::Power(beta, p3);
769 const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
771 //printf("test----- %f %f -- %f %f %f %f %f --- %f %f %f\n", bg, beta, p0, p1, p2, p3, p4, p0/aa, aa, bb);
773 return (p1-aa-bb)*p0/aa;
776 Double_t AliTRDdEdxBaseUtils::ToLogx(FFunc func, const Double_t * xx, const Double_t * par)
779 //f(x)-> f(y) with y=log10(x)
781 const Double_t x2[]={TMath::Power(10, xx[0])};
782 return func(x2, par);