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();
75 //printf("AliTRDdEdxBaseUtils::BinLogX warning xmin < epsilon! nothing done, axis not set. %e\n", from); exit(1);
79 Double_t *new_bins = new Double_t[bins + 1];
82 const Double_t factor = pow(to/from, 1./bins);
84 for (int i = 1; i <= bins; i++) {
85 new_bins[i] = factor * new_bins[i-1];
87 axis->Set(bins, new_bins);
91 void AliTRDdEdxBaseUtils::GetCDFCuts(const TH1D *hh, const Int_t ncut, Double_t cuts[], const Double_t cdfs[], const Double_t thres)
94 //counts of hh is sorted
97 for(Int_t ii=0; ii<ncut; ii++){
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);
117 TMath::Sort(nsel, datas, id, kFALSE);
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);
124 cuts[ii] = datas[id[Int_t(icdf*nsel)]];
128 Double_t AliTRDdEdxBaseUtils::GetMeanRMS(const Double_t nn, const Double_t sum, const Double_t w2s, Double_t * grms, Double_t * gerr)
131 //calculate mean (with error) and rms from sum, w2s, nn
132 //if nn=0, mean, error, and rms all = 0
135 Double_t tmean = 0, trms = 0, terr = 0;
140 const Double_t arg = w2s/nn-tmean*tmean;
141 if(TMath::Abs(arg)<EPSILON){
146 printf("AliTRDdEdxBaseUtils::GetMeanRMS error negative sqrt argument!! %e -- %e %e %f\n", arg, w2s, sum, nn); exit(1);
149 trms = TMath::Sqrt(arg);
152 terr = trms/TMath::Sqrt(nn);
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)
169 //calculate truncated mean
170 //return <x*w>_{low-high according to x}
175 for(Int_t ii=0; ii<nx; ii++){
176 printf("test %d/%d %f\n", ii, nx, xdata[ii]);
182 TMath::Sort(nx, xdata, index, kFALSE);
187 const Int_t istart = Int_t (nx*lowfrac);
188 const Int_t istop = Int_t (nx*highfrac);
190 //=,< correct, because when low=0, high=1 it is correct
191 for(Int_t ii=istart; ii<istop; ii++){
194 weight = wws[index[ii]];
196 const Double_t sx = xdata[index[ii]]*weight;
202 //printf("test in loop %d/%d %f %f %f\n", ii, nused, sx, sum, w2s);
206 return GetMeanRMS(nused, sum, w2s, grms, gerr);
209 Double_t AliTRDdEdxBaseUtils::TruncatedMean(const TH1 *hh, const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr)
212 //do truncation on histogram
214 //if hh is scaled, be sure Sumw2 is called before scaling!! then mean, rms and err will all be correct
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);
223 printf("AliTRDdEdxBaseUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
227 npreTrunc += bc*bc/be/be;
230 const Double_t nstart = npreTrunc*lowfrac;
231 const Double_t nstop = npreTrunc*highfrac;
233 //with Double_t this should also handle normalized hist
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);
243 printf("AliTRDdEdxBaseUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
247 const Double_t weight = bc*bc/be/be;
249 //<= correct, because when high=1, nstop has to be included
250 if(ntot>nstart && ntot<=nstop){
252 const Double_t val = hh->GetBinCenter(itmp);
254 w2s += weight*val*val;
258 //printf("test %d %f %f --- %f %f -- %f %f\n", itmp, weight, val, sum, w2s, nused, nsample);
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);
270 return GetMeanRMS(nused, sum, w2s, grms, gerr);
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)
276 //fit slices of hh using truncation
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();
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();
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);
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);
302 for(Int_t ix=x0; ix<=x1; ix++){
304 const Double_t rawcount = hh->Integral(ix,ix,0, ny+1);
305 if(rawcount<EPSILON){
309 TH1D *htmp = new TH1D(Form("FitSlicesY_%s_%d", hh->GetName(), ix),"",ny, ymin, ymax);
310 htmp->GetXaxis()->Set(ny, newbins);
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);
318 printf("AliTRDdEdxBaseUtils::FitSlicesY error %d %d %e %e\n", ix, iy, be, bc); exit(1);
323 htmp->SetBinContent(iy, bc);
324 htmp->SetBinError(iy, be);
326 ntot += (bc/be)*(bc/be);
328 //if(be) printf("test %d %d : %f %f %f\n", ix, iy, bc, be, pow(bc/be,2));
331 hnor->SetBinContent(ix, ntot);
332 hnor->SetBinError( ix, 0);
334 if(ntot<thres || htmp->GetRMS()<EPSILON){
340 Double_t trms = -999, terr = -999;
341 const Double_t tmean = TruncatedMean(htmp, lowfrac, highfrac, &trms, &terr);
343 hmpv->SetBinContent(ix, tmean);
344 hmpv->SetBinError( ix, terr);
346 hwid->SetBinContent(ix, trms);
347 hwid->SetBinError( ix, 0);
349 hres->SetBinContent(ix, tmean>EPSILON ? trms/tmean:0);
350 hres->SetBinError( ix, 0);
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());
366 //===================================================================================
367 // TRD Analysis Fast Tool
368 //===================================================================================
370 Int_t AliTRDdEdxBaseUtils::GetNtracklet(const AliESDEvent *esd)
373 //number of trd tracklet in one esd event
375 const Int_t ntrk0 = esd->GetNumberOfTracks();
377 for(Int_t ii=0; ii<ntrk0; ii++){
378 ntrdv1 += esd->GetTrack(ii)->GetTRDntracklets();
383 AliTRDtrackV1 * AliTRDdEdxBaseUtils::GetTRDtrackV1(const AliESDtrack * esdtrack)
386 //Get TRD friend track
389 AliESDfriendTrack * friendtrk = (AliESDfriendTrack *)esdtrack->GetFriendTrack();
391 //printf("xlulog AliAnalysisTaskCosmicTRD::GetTRDtrack no friend!!\n"); exit(1);
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)) )
405 Bool_t AliTRDdEdxBaseUtils::IsInSameStack(const AliTRDtrackV1 *trdtrack)
408 // to check if all tracklets are in the same stack, useful in cosmic
411 TVectorD secs(18), stks(5);
413 for(Int_t ilayer = 0; ilayer < 6; ilayer++){
414 AliTRDseedV1 *tracklet=trdtrack->GetTracklet(ilayer);
418 const Int_t det = tracklet->GetDetector();
419 const Int_t isector = AliTRDgeometry::GetSector(det);
420 const Int_t istack = AliTRDgeometry::GetStack(det);
426 if(secs.Sum()!=1 || stks.Sum()!=1){
433 Double_t AliTRDdEdxBaseUtils::GetRedefinedPhi(Double_t phi)
436 //redefine phi in pi/2 - pi/2 + 2pi
439 if(phi>0 && phi < TMath::Pi()/2){
440 phi += 2*TMath::Pi();
446 Double_t AliTRDdEdxBaseUtils::Getdydx(const AliTRDseedV1 *tracklet)
452 return tracklet->GetYref(1);
457 Double_t AliTRDdEdxBaseUtils::Getdzdx(const AliTRDseedV1 *tracklet)
463 return tracklet->GetZref(1);
468 Double_t AliTRDdEdxBaseUtils::Getdldx(const AliTRDseedV1 *tracklet)
471 //return angular normalization factor = dldx
474 return TMath::Sqrt(1+tracklet->GetYref(1)*tracklet->GetYref(1)+tracklet->GetZref(1)*tracklet->GetZref(1));
479 AliTRDseedV1 * AliTRDdEdxBaseUtils::GetFirstTracklet(const AliTRDtrackV1 *trdtrack)
485 AliTRDseedV1 *tracklet = 0x0;
486 for(Int_t ilayer = 0; ilayer < 6; ilayer++){
487 tracklet=trdtrack->GetTracklet(ilayer);
497 AliTRDseedV1 * AliTRDdEdxBaseUtils::GetLastTracklet(const AliTRDtrackV1 *trdtrack)
502 AliTRDseedV1 *tracklet = 0x0;
504 for(Int_t ilayer = 5; ilayer >= 0; ilayer--){
505 tracklet=trdtrack->GetTracklet(ilayer);
515 void AliTRDdEdxBaseUtils::GetFirstSectorStackMomentum(const AliTRDtrackV1 *trdtrack, Int_t & isec, Int_t & istk, Double_t & mom)
523 AliTRDseedV1 *tracklet = GetFirstTracklet(trdtrack);
525 const Int_t det = tracklet->GetDetector();
526 isec = AliTRDgeometry::GetSector(det);
527 istk = AliTRDgeometry::GetStack(det);
529 mom = tracklet->GetMomentum();
535 //===================================================================================
536 // Detector and Data Constant
537 //===================================================================================
538 Int_t AliTRDdEdxBaseUtils::ToDetector(const Int_t gtb)
543 return gtb/AliTRDseedV1::kNtb;
546 Int_t AliTRDdEdxBaseUtils::ToTimeBin(const Int_t gtb)
551 return gtb%AliTRDseedV1::kNtb;
554 Int_t AliTRDdEdxBaseUtils::ToSector(const Int_t gtb)
559 return AliTRDgeometry::GetSector(ToDetector(gtb));
562 Int_t AliTRDdEdxBaseUtils::ToStack(const Int_t gtb)
567 return AliTRDgeometry::GetStack(ToDetector(gtb));
570 Int_t AliTRDdEdxBaseUtils::ToLayer(const Int_t gtb)
575 return AliTRDgeometry::GetLayer(ToDetector(gtb));
578 void AliTRDdEdxBaseUtils::CheckRunB(const TString listrun1kg, const Int_t run, TString & type)
583 if(listrun1kg.Contains(Form("%d",run))){
591 TString AliTRDdEdxBaseUtils::GetRunType(const Int_t run)
598 if(run>=121527 && run<= 126460)//LHC10d
600 else if(run>=126461 && run<= 130930)//LHC10e
602 else if(run>=136782 && run <= 139846)//LHC10h
603 type="2010PbPbLHC10h";
604 else if(run>= 143224 && run<= 143237)//2011Feb
605 type="2011cosmicFeb";
606 else if(run>= 150587 && run<= 154930){
607 type="2011cosmicMayJun";
609 CheckRunB("154601 154602 154629 154634 154636 154639 154643", run, type);
611 else if(run>=173047 && run<=180164){
614 CheckRunB("173047 173049 173050 173065 173070 173092 173095 173113 173131 173160 174154 174156 174192 174194", run, type);
624 void AliTRDdEdxBaseUtils::PrintControl()
627 //print out control variable
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());
632 //===================================================================================
633 // dEdx Parameterization
634 //===================================================================================
635 void AliTRDdEdxBaseUtils::FastFitdEdxTR(TH1 * hh)
638 //fit dedx tr from mean
641 TF1 *ff=new TF1("ff", MeandEdxTRLogx, -0.5, 4.5, 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;
652 ff->SetParameters(par);
655 for(int ii=0; ii<8; ii++){
656 printf("par[%d]=%e;\n", ii, ff->GetParameter(ii));
659 TFile *fout=new TFile("fastfit.root","recreate");
669 Double_t AliTRDdEdxBaseUtils::Q0MeanTRDpp(const Double_t bg)
672 //truncated Mean Q_{xx} in TRD
678 //2012-0605 /u/xlu/.task/CommondEdx/myAnaData/newTune_r56595/inuse/plot/fastFit
685 par[5]=-1.495764e-06;
686 par[6]=-2.536119e+00;
689 return scale*MeandEdxTR(&bg, par);
692 Double_t AliTRDdEdxBaseUtils::Q0MeanTRDPbPb(const Double_t bg)
695 //truncated Mean Q_{xx} in TRD
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;
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);
715 Double_t AliTRDdEdxBaseUtils::Q1MeanTRDpp(const Double_t bg)
718 //truncated Mean 1/(1/Q)_{xx} in TRD
723 //So 4. Mär 13:30:51 CET 2012
725 par[0]= 2.434646e-01;
726 par[1]= 1.400211e+00;
727 par[2]= 6.937471e+00;
728 par[3]= 7.758118e-02;
729 par[4]= 1.097372e+01;
730 par[5]= 4.297518e-04;
731 par[6]= -1.806266e+00;
732 par[7]= 1.543811e+00;
734 //hhtype2Q1b2c2 scale 0.418629 at ltbg 0.650000
736 return 0.418629*MeandEdxTR(&bg, par);
739 Double_t AliTRDdEdxBaseUtils::Q1MeanTRDPbPb(const Double_t bg)
742 //truncated Mean 1/(1/Q)_{xx} in TRD
747 //So 4. Mär 13:30:52 CET 2012
749 par[0]= 2.193660e-01;
750 par[1]= 2.051864e+00;
751 par[2]= 6.825112e+00;
752 par[3]= 6.151693e-02;
753 par[4]= 1.390343e+01;
754 par[5]= 6.010032e-05;
755 par[6]= -1.676324e+00;
756 par[7]= 1.838873e+00;
758 //hhtype4Q1b2c2 scale 0.457988 at ltbg 0.650000
760 return 0.457988*MeandEdxTR(&bg, par);
763 Double_t AliTRDdEdxBaseUtils::QMeanTPC(const Double_t bg)
770 //Mi 15. Feb 14:48:05 CET 2012
771 //train_2012-02-13_1214.12001, tpcsig
778 return MeandEdx(&bg, par);
781 Double_t AliTRDdEdxBaseUtils::MeandEdxTR(const Double_t * xx, const Double_t * pin)
784 //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
788 for(int ii=0; ii<3; ii++){
791 return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
794 Double_t AliTRDdEdxBaseUtils::MeanTR(const Double_t * xx, const Double_t * par)
797 //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
801 const Double_t bg = xx[0];
802 const Double_t gamma = sqrt(1+bg*bg);
804 const Double_t p0 = TMath::Abs(par[1]);
805 const Double_t p1 = TMath::Abs(par[2]);
806 const Double_t p2 = TMath::Abs(par[3]);
808 const Double_t zz = TMath::Log(gamma);
809 const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
811 return par[0]+tryield;
814 Double_t AliTRDdEdxBaseUtils::MeandEdx(const Double_t * xx, const Double_t * par)
816 return ALEPH(xx, par);
819 Double_t AliTRDdEdxBaseUtils::ALEPH(const Double_t * xx, const Double_t * par)
822 //ALEPH parametrization for dEdx
826 const Double_t bg = xx[0];
827 const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
829 const Double_t p0 = TMath::Abs(par[0]);
830 const Double_t p1 = TMath::Abs(par[1]);
832 //after redefining p2 (<0) -> ln P2
833 //const Double_t p2 = par[2];
836 const Double_t p2 = TMath::Abs(par[2]);
838 const Double_t p3 = TMath::Abs(par[3]);
839 const Double_t p4 = TMath::Abs(par[4]);
841 const Double_t aa = TMath::Power(beta, p3);
843 //numerically very bad when p2~1e-15, bg~1e3, p4~5.
844 const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
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!!
852 //----> fail for 10h, not used, restore back!
853 //const Double_t lbg = TMath::Log(bg);
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) );
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);
861 return (p1-aa-bb)*p0/aa;
864 Double_t AliTRDdEdxBaseUtils::ToLogx(FFunc func, const Double_t * xx, const Double_t * par)
867 //f(x)-> f(y) with y=log10(x)
869 const Double_t x2[]={TMath::Power(10, xx[0])};
870 return func(x2, par);