#include "AliTRDdEdxBaseUtils.h"
-#define EPSILON 1e-12
+#define EPSILON 1e-8
Double_t AliTRDdEdxBaseUtils::fgQ0Frac = 0.55;
Double_t AliTRDdEdxBaseUtils::fgQ1Frac = 0.5;
const Double_t from = axis->GetXmin();
const Double_t to = axis->GetXmax();
- if (from<EPSILON) return;
+ if (from<EPSILON){
+ //printf("AliTRDdEdxBaseUtils::BinLogX warning xmin < epsilon! nothing done, axis not set. %e\n", from); exit(1);
+ return;
+ }
+
Double_t *new_bins = new Double_t[bins + 1];
new_bins[0] = from;
delete [] new_bins;
}
-void AliTRDdEdxBaseUtils::GetCDFCuts(const TH1D *hh, const Int_t ncut, Double_t cuts[], const Double_t cdfs[], const Double_t thres)
+void AliTRDdEdxBaseUtils::GetCDFCuts(const TH1D *hh, Int_t ncut, Double_t cuts[], const Double_t cdfs[], Double_t thres)
{
//
//counts of hh is sorted
}
}
-Double_t AliTRDdEdxBaseUtils::GetMeanRMS(const Double_t nn, const Double_t sum, const Double_t w2s, Double_t * grms, Double_t * gerr)
+Double_t AliTRDdEdxBaseUtils::GetMeanRMS(Double_t nn, Double_t sum, Double_t w2s, Double_t * grms, Double_t * gerr)
{
//
//calculate mean (with error) and rms from sum, w2s, nn
return tmean;
}
-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)
+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)
{
//
//calculate truncated mean
return GetMeanRMS(nused, sum, w2s, grms, gerr);
}
-Double_t AliTRDdEdxBaseUtils::TruncatedMean(const TH1 *hh, const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr)
+Double_t AliTRDdEdxBaseUtils::TruncatedMean(const TH1 *hh, Double_t lowfrac, Double_t highfrac, Double_t * grms, Double_t * gerr)
{
//
//do truncation on histogram
return GetMeanRMS(nused, sum, w2s, grms, gerr);
}
-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)
+void AliTRDdEdxBaseUtils::FitSlicesY(const TH2D *hh, TH1D *&hnor, TH1D *&hmpv, TH1D *&hwid, TH1D *&hres, Double_t thres, Double_t lowfrac, Double_t highfrac)
{
//
//fit slices of hh using truncation
return kTRUE;
}
-Bool_t AliTRDdEdxBaseUtils::GetFirstSectorStackMomentum(const AliTRDtrackV1 *trdtrack, Int_t & isec, Int_t & istk, Double_t & mom)
+Double_t AliTRDdEdxBaseUtils::GetRedefinedPhi(Double_t phi)
+{
+ //
+ //redefine phi in pi/2 - pi/2 + 2pi
+ //
+
+ if(phi>0 && phi < TMath::Pi()/2){
+ phi += 2*TMath::Pi();
+ }
+
+ return phi;
+}
+
+Double_t AliTRDdEdxBaseUtils::Getdydx(const AliTRDseedV1 *tracklet)
+{
+ //
+ //get dydx
+ //
+ if(tracklet)
+ return tracklet->GetYref(1);
+ else
+ return -999;
+}
+
+Double_t AliTRDdEdxBaseUtils::Getdzdx(const AliTRDseedV1 *tracklet)
+{
+ //
+ //get dzdx
+ //
+ if(tracklet)
+ return tracklet->GetZref(1);
+ else
+ return -999;
+}
+
+Double_t AliTRDdEdxBaseUtils::Getdldx(const AliTRDseedV1 *tracklet)
+{
+ //
+ //return angular normalization factor = dldx
+ //
+ if(tracklet)
+ return TMath::Sqrt(1+tracklet->GetYref(1)*tracklet->GetYref(1)+tracklet->GetZref(1)*tracklet->GetZref(1));
+ else
+ return -999;
+}
+
+AliTRDseedV1 * AliTRDdEdxBaseUtils::GetFirstTracklet(const AliTRDtrackV1 *trdtrack)
{
//
//as function name
//
- isec = istk = -999;
- mom = -999;
+ AliTRDseedV1 *tracklet = 0x0;
for(Int_t ilayer = 0; ilayer < 6; ilayer++){
- AliTRDseedV1 *tracklet=trdtrack->GetTracklet(ilayer);
+ tracklet=trdtrack->GetTracklet(ilayer);
if(!tracklet)
continue;
+ break;
+ }
+
+ return tracklet;
+}
+
+AliTRDseedV1 * AliTRDdEdxBaseUtils::GetLastTracklet(const AliTRDtrackV1 *trdtrack)
+{
+ //
+ //get last tracklet
+ //
+ AliTRDseedV1 *tracklet = 0x0;
+
+ for(Int_t ilayer = 5; ilayer >= 0; ilayer--){
+ tracklet=trdtrack->GetTracklet(ilayer);
+ if(!tracklet)
+ continue;
+
+ break;
+ }
+
+ return tracklet;
+}
+
+void AliTRDdEdxBaseUtils::GetFirstSectorStackMomentum(const AliTRDtrackV1 *trdtrack, Int_t & isec, Int_t & istk, Double_t & mom)
+{
+ //
+ //as function name
+ //
+ isec = istk = -999;
+ mom = -999;
+
+ AliTRDseedV1 *tracklet = GetFirstTracklet(trdtrack);
+ if(tracklet){
const Int_t det = tracklet->GetDetector();
isec = AliTRDgeometry::GetSector(det);
istk = AliTRDgeometry::GetStack(det);
-
+
mom = tracklet->GetMomentum();
-
- break;
}
- if(isec<0)
- return kFALSE;
- else
- return kTRUE;
+ return;
}
//===================================================================================
// Detector and Data Constant
//===================================================================================
-Int_t AliTRDdEdxBaseUtils::ToDetector(const Int_t gtb)
+Int_t AliTRDdEdxBaseUtils::ToDetector(Int_t gtb)
{
//
//gtb = det*Ntb+itb
return gtb/AliTRDseedV1::kNtb;
}
-Int_t AliTRDdEdxBaseUtils::ToTimeBin(const Int_t gtb)
+Int_t AliTRDdEdxBaseUtils::ToTimeBin(Int_t gtb)
{
//
//gtb = det*Ntb+itb
return gtb%AliTRDseedV1::kNtb;
}
-Int_t AliTRDdEdxBaseUtils::ToSector(const Int_t gtb)
+Int_t AliTRDdEdxBaseUtils::ToSector(Int_t gtb)
{
//
//return sector
return AliTRDgeometry::GetSector(ToDetector(gtb));
}
-Int_t AliTRDdEdxBaseUtils::ToStack(const Int_t gtb)
+Int_t AliTRDdEdxBaseUtils::ToStack(Int_t gtb)
{
//
//return stack
return AliTRDgeometry::GetStack(ToDetector(gtb));
}
-Int_t AliTRDdEdxBaseUtils::ToLayer(const Int_t gtb)
+Int_t AliTRDdEdxBaseUtils::ToLayer(Int_t gtb)
{
//
//return layer
return AliTRDgeometry::GetLayer(ToDetector(gtb));
}
-TString AliTRDdEdxBaseUtils::GetRunType(const Int_t run)
+void AliTRDdEdxBaseUtils::CheckRunB(TString listrun1kg, Int_t run, TString & type)
+{
+ //
+ //check run b field
+ //
+ if(listrun1kg.Contains(Form("%d",run))){
+ type+="1kG";
+ }
+ else {
+ type+="5kG";
+ }
+}
+
+TString AliTRDdEdxBaseUtils::GetRunType(Int_t run)
{
//
//return run type
TString type;
if(run>=121527 && run<= 126460)//LHC10d
- type="pp2010LHC10d";
+ type="2010ppLHC10d";
else if(run>=126461 && run<= 130930)//LHC10e
- type="pp2010LHC10e";
+ type="2010ppLHC10e";
else if(run>=136782 && run <= 139846)//LHC10h
- type="PbPb2010LHC10h";
+ type="2010PbPbLHC10h";
else if(run>= 143224 && run<= 143237)//2011Feb
- type="cosmic2011Feb";
+ type="2011cosmicFeb";
else if(run>= 150587 && run<= 154930){
- type="cosmic2011MayJun";
+ type="2011cosmicMayJun";
- TString runstr(Form("%d",run));
- const TString listrun1kg("154601 154602 154629 154634 154636 154639 154643");
- if(listrun1kg.Contains(runstr)){
- type+="1kG";
- }
- else{
- type+="5kG";
- }
+ CheckRunB("154601 154602 154629 154634 154636 154639 154643", run, type);
+ }
+ else if(run>=173047 && run<=180164){
+ type="2012cosmic";
+
+ CheckRunB("173047 173049 173050 173065 173070 173092 173095 173113 173131 173160 174154 174156 174192 174194", run, type);
}
else{
type="unknown";
//===================================================================================
// dEdx Parameterization
//===================================================================================
+void AliTRDdEdxBaseUtils::FastFitdEdxTR(TH1 * hh)
+{
+ //
+ //fit dedx tr from mean
+ //
+
+ TF1 *ff=new TF1("ff", MeandEdxTRLogx, -0.5, 4.5, 8);
+ Double_t par[8];
+ par[0]= 2.397001e-01;
+ par[1]= 1.334697e+00;
+ par[2]= 6.967470e+00;
+ par[3]= 9.055289e-02;
+ par[4]= 9.388760e+00;
+ par[5]= 9.452742e-04;
+ par[6]= -1.866091e+00;
+ par[7]= 1.403545e+00;
+
+ ff->SetParameters(par);
+ hh->Fit(ff,"N");
+
+ for(int ii=0; ii<8; ii++){
+ printf("par[%d]=%e;\n", ii, ff->GetParameter(ii));
+ }
+
+ TFile *fout=new TFile("fastfit.root","recreate");
+ hh->Write();
+ ff->Write();
+ fout->Save();
+ fout->Close();
+
+ delete fout;
+ delete ff;
+}
-Double_t AliTRDdEdxBaseUtils::Q0MeanTRDpp(const Double_t bg)
+Double_t AliTRDdEdxBaseUtils::Q0MeanTRDpp(Double_t bg)
{
//
//truncated Mean Q_{xx} in TRD
//
Double_t par[8];
- //03132012161150
- //opt: ppQ0
-par[0]= 2.397001e-01;
-par[1]= 1.334697e+00;
-par[2]= 6.967470e+00;
-par[3]= 9.055289e-02;
-par[4]= 9.388760e+00;
-par[5]= 9.452742e-04;
-par[6]= -1.866091e+00;
-par[7]= 1.403545e+00;
-
- ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see2.log:hhtype2Q0b2c2 scale 0.428543 at ltbg 0.650000
- return 0.428543*MeandEdxTR(&bg, par);
+ Double_t scale = 1;
+
+ //2012-0605 /u/xlu/.task/CommondEdx/myAnaData/newTune_r56595/inuse/plot/fastFit
+ scale = 0.9257;
+par[0]=8.316476e-01;
+par[1]=1.600907e+00;
+par[2]=7.728447e+00;
+par[3]=6.235622e-02;
+par[4]=1.136209e+01;
+par[5]=-1.495764e-06;
+par[6]=-2.536119e+00;
+par[7]=3.416031e+00;
+
+ return scale*MeandEdxTR(&bg, par);
}
-Double_t AliTRDdEdxBaseUtils::Q0MeanTRDPbPb(const Double_t bg)
+Double_t AliTRDdEdxBaseUtils::Q0MeanTRDPbPb(Double_t bg)
{
//
//truncated Mean Q_{xx} in TRD
return 0.460994*MeandEdxTR(&bg, par);
}
-Double_t AliTRDdEdxBaseUtils::Q1MeanTRDpp(const Double_t bg)
+Double_t AliTRDdEdxBaseUtils::Q1MeanTRDpp(Double_t bg)
{
//
//truncated Mean 1/(1/Q)_{xx} in TRD
return 0.418629*MeandEdxTR(&bg, par);
}
-Double_t AliTRDdEdxBaseUtils::Q1MeanTRDPbPb(const Double_t bg)
+Double_t AliTRDdEdxBaseUtils::Q1MeanTRDPbPb(Double_t bg)
{
//
//truncated Mean 1/(1/Q)_{xx} in TRD
return 0.457988*MeandEdxTR(&bg, par);
}
-Double_t AliTRDdEdxBaseUtils::QMeanTPC(const Double_t bg)
+Double_t AliTRDdEdxBaseUtils::QMeanTPC(Double_t bg)
{
//
//bethe bloch in TPC
return MeandEdx(&bg, par);
}
-Double_t AliTRDdEdxBaseUtils::MeandEdxTR(const Double_t * xx, const Double_t * pin)
+Double_t AliTRDdEdxBaseUtils::MeandEdxTR(const Double_t * xx, const Double_t * pin)
{
//
//ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
}
-Double_t AliTRDdEdxBaseUtils::MeanTR(const Double_t * xx, const Double_t * par)
+Double_t AliTRDdEdxBaseUtils::MeanTR(const Double_t * xx, const Double_t * par)
{
//
//ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
return par[0]+tryield;
}
-Double_t AliTRDdEdxBaseUtils::MeandEdx(const Double_t * xx, const Double_t * par)
+Double_t AliTRDdEdxBaseUtils::MeandEdx(const Double_t * xx, const Double_t * par)
+{
+ return ALEPH(xx, par);
+}
+
+Double_t AliTRDdEdxBaseUtils::ALEPH(const Double_t * xx, const Double_t * par)
{
//
//ALEPH parametrization for dEdx
const Double_t p0 = TMath::Abs(par[0]);
const Double_t p1 = TMath::Abs(par[1]);
+
+ //after redefining p2 (<0) -> ln P2
+ //const Double_t p2 = par[2];
+
+ //restore back
const Double_t p2 = TMath::Abs(par[2]);
+
const Double_t p3 = TMath::Abs(par[3]);
const Double_t p4 = TMath::Abs(par[4]);
const Double_t aa = TMath::Power(beta, p3);
+
+ //numerically very bad when p2~1e-15, bg~1e3, p4~5.
const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
+ //+1 to avoid the singularity when bg<1 and p4>>1
+ //very^inf important! with this hesse NOTPOS->OK and no nan or inf in migrad
+ //i.e. numerically very stable
+ //the reason is when bg<1 ln blows up very fast with p4>>1 and nan/inf appears in migrad search
+ //the fitting will adjust the parameters as if bg is not shifted, the fitted curve looks fine!!
+ //-- 2012 Aug 8
+ //----> fail for 10h, not used, restore back!
+ //const Double_t lbg = TMath::Log(bg);
+
+ //redefine p2(<0) -> ln P2
+ //corresponds to alephFit.C fAlephOPt = 11
+ //const Double_t bb = p2 + TMath::Log( 1 + TMath::Exp(-p2-p4*lbg) );
+
//printf("test----- %f %f -- %f %f %f %f %f --- %f %f %f\n", bg, beta, p0, p1, p2, p3, p4, p0/aa, aa, bb);
return (p1-aa-bb)*p0/aa;
}
-Double_t AliTRDdEdxBaseUtils::ToLogx(FFunc func, const Double_t * xx, const Double_t * par)
+Double_t AliTRDdEdxBaseUtils::ToLogx(FFunc func, const Double_t * xx, const Double_t * par)
{
//
//f(x)-> f(y) with y=log10(x)
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff