#include "TH2D.h"
#include "TSpline.h"
#include "TF1.h"
-#include "AliTOFGeometry.h"
-#include "AliTOFT0maker.h"
#include "AliCentrality.h"
#include "AliAODPid.h"
#include "AliMCEvent.h"
#include "AliAnalysisManager.h"
#include "AliAODMCHeader.h"
#include "AliAODMCParticle.h"
+#include "TH1D.h"
+#include "TFile.h"
+
ClassImp(AliFlowBayesianPID)
-TH2D* AliFlowBayesianPID::fghPriors[fgkNspecies] = {NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL}; // histo with priors (hardcoded)
+TH2D* AliFlowBayesianPID::fghPriors[fgkNspecies] = {NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL}; // histo with priors (hardcoded)
TSpline3* AliFlowBayesianPID::fgMism = NULL; // function for mismatch
-AliTOFGeometry* AliFlowBayesianPID::fgTofGeo = NULL; // TOF geometry needed to reproduce mismatch shape
+TH1D* AliFlowBayesianPID::fgHtofChannelDist=NULL;
//________________________________________________________________________
AliFlowBayesianPID::AliFlowBayesianPID(AliESDpid *esdpid)
- : AliPIDResponse(), fPIDesd(NULL), fDB(TDatabasePDG::Instance()), fNewTrackParam(0), fTOFresolution(84.0), fTOFResponseF(NULL), fTPCResponseF(NULL), fTOFmaker(NULL),fWTofMism(0.0), fProbTofMism(0.0), fZ(0) ,fMassTOF(0), fBBdata(NULL),fCurrCentrality(100),fPsi(999),fPsiRes(999),fIsMC(kFALSE),fDedx(0.0)
+ : AliPIDResponse(), fPIDesd(NULL), fDB(TDatabasePDG::Instance()), fNewTrackParam(0), fTOFresolution(84.0), fTOFResponseF(NULL), fTPCResponseF(NULL),fWTofMism(0.0), fProbTofMism(0.0), fZ(0) ,fMassTOF(0), fBBdata(NULL),fCurrCentrality(100),fPsi(999),fPsiRes(999),fIsMC(kFALSE),fDedx(0.0)
{
// Constructor
Bool_t redopriors = kFALSE;
fghPriors[7] = new TH2D("fghPriorsHe","Priors as a function of Centrality [p];centrality;p_{t} (GeV/c)",12,-10,110,80,0,20);
redopriors = kTRUE;
}
+ if(! fghPriors[8]){
+ fghPriors[8] = new TH2D("fghPriorsAlfa","Priors as a function of Centrality [p];centrality;p_{t} (GeV/c)",12,-10,110,80,0,20);
+ redopriors = kTRUE;
+ }
if(redopriors) SetPriors();
if(!fgMism) fgMism = GetMismatch();
- if(! fgTofGeo) fgTofGeo = new AliTOFGeometry();
if(! esdpid)
fPIDesd = new AliESDpid();
else
fPIDesd = esdpid;
-
- fTOFmaker = new AliTOFT0maker(fPIDesd);
fProb[0]=0.0;
fProb[1]=0.0;
fProb[5]=0.0;
fProb[6]=0.0;
fProb[7]=0.0;
+ fProb[8]=0.0;
fMass[0] = fDB->GetParticle(11)->Mass(); // e mass
fMass[1] = fDB->GetParticle(13)->Mass(); // mu mass
fMass[5] = fDB->GetParticle(2212)->Mass()+fDB->GetParticle(2112)->Mass(); // p mass
fMass[6] = fDB->GetParticle(2212)->Mass()+fDB->GetParticle(2112)->Mass()*2; // p mass
fMass[7] = (fDB->GetParticle(2212)->Mass()+fDB->GetParticle(2112)->Mass()*2)*0.5; // p mass
+ fMass[8] = (fDB->GetParticle(2212)->Mass()*2+fDB->GetParticle(2112)->Mass()*2)*0.5; // p mass
// TOF response
fTOFResponseF = new TF1("fTOFprob","[0]*TMath::Exp(-(x-[1])*(x-[1])/2/[2]/[2])* (x < [1]+[3]*[2]) + (x > [1]+[3]*[2])*[0]*TMath::Exp(-(x-[1]-[3]*[2]*0.5)*[3]/[2])",-7,7);
fMaskOR[i] = 1; // all dets if available
fMaskCurrent[i] = 0; // current mask
}
+
+ if(!fgHtofChannelDist){
+ TFile *ftofchannel = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
+ fgHtofChannelDist = (TH1D *) ftofchannel->Get("hTOFchanDist");
+ }
+
}
//________________________________________________________________________
AliFlowBayesianPID::~AliFlowBayesianPID(){
// destrucctor that removes all the PID external (non static) objects
if(fTOFResponseF) delete fTOFResponseF;
if(fTPCResponseF) delete fTPCResponseF;
- if(fTOFmaker) delete fTOFmaker;
if(fBBdata) delete fBBdata;
}
//________________________________________________________________________
-void AliFlowBayesianPID::SetDetResponse(AliESDEvent *esd,Float_t centrality,EStartTimeType_t flagStart,Bool_t recomputeT0TOF){
+void AliFlowBayesianPID::SetDetResponse(AliESDEvent *esd,Float_t centrality,EStartTimeType_t flagStart,Bool_t){
// Set the detector responses (including also TPC dE/dx paramterization vs. centrality)
if(!esd){
printf("AliFlowBayesianPID::SetDetResponse -> Error -> No valid esd event");
fBBdata->SetParameter(4, alephParameters[3]);
fBBdata->SetParameter(5, alephParameters[4]);
- if(recomputeT0TOF){
- fTOFmaker->SetTimeResolution(fTOFresolution);
- fTOFmaker->ComputeT0TOF(esd);
- fTOFmaker->WriteInESD(esd);
- }
-
fPIDesd->SetTOFResponse(esd,flagStart);
if(fNewTrackParam){
fPsiRes=999;
}
//________________________________________________________________________
-Float_t AliFlowBayesianPID::GetExpDeDx(const AliESDtrack *t,Int_t iS) const{
+//________________________________________________________________________
+Float_t AliFlowBayesianPID::GetExpDeDx(const AliVTrack *t,Int_t iS) const{
// tuned dE/dx (vs. eta and centrality)
- Double_t ptpc[3];
- t->GetInnerPxPyPz(ptpc);
- Float_t momtpc=TMath::Sqrt(ptpc[0]*ptpc[0] + ptpc[1]*ptpc[1] + ptpc[2]*ptpc[2]);
+ Float_t momtpc=t->GetTPCmomentum();
Float_t dedxExp=0;
if(iS==0) dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kElectron);
else if(iS==5) dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kDeuteron);
else if(iS==6) dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kTriton);
else if(iS==7) dedxExp = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[7])*5;
-
+ else if(iS==8) dedxExp = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[8])*5;
+
Float_t eta = t->Eta();
Float_t etaCorr = 7.98368e-03 - 1.67208e-02 - 1.89776e-01*eta*eta -2.90836e-02*eta*eta + 5.96093e-01*eta*eta*eta*eta + 6.06450e-02*eta*eta*eta*eta - 3.55884e-01*eta*eta*eta*eta*eta*eta;
if(fCurrCentrality < 0){
return dedxExp;
}
//________________________________________________________________________
-Float_t AliFlowBayesianPID::GetExpDeDx(const AliAODTrack *t,Int_t iS) const{
+Float_t AliFlowBayesianPID::GetExpDeDx(const AliVTrack *t,Float_t mass) const{
// tuned dE/dx (vs. eta and centrality)
Float_t momtpc=t->GetTPCmomentum();
- Float_t dedxExp=0;
- if(iS==0) dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kElectron);
- else if(iS==1) dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kMuon);
- else if(iS==2) dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kPion);
- else if(iS==3) dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kKaon);
- else if(iS==4) dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kProton);
- else if(iS==5) dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kDeuteron);
- else if(iS==6) dedxExp = fPIDesd->GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kTriton);
- else if(iS==7) dedxExp = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[7])*5;
+ if(mass < 0.0001) mass = 0.0001;
+
+ Float_t dedxExp = fPIDesd->GetTPCResponse().Bethe(momtpc/mass);
Float_t eta = t->Eta();
Float_t etaCorr = 7.98368e-03 - 1.67208e-02 - 1.89776e-01*eta*eta -2.90836e-02*eta*eta + 5.96093e-01*eta*eta*eta*eta + 6.06450e-02*eta*eta*eta*eta - 3.55884e-01*eta*eta*eta*eta*eta*eta;
+
if(fCurrCentrality < 0){
}
else if(fCurrCentrality < 5) etaCorr += 17E-3;
resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kProton);
dedxExp = GetExpDeDx(t,4);
}
+ else if(iS-1e+9 == 10020){ // d
+ resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kDeuteron);
+ dedxExp = GetExpDeDx(t,5);
+ }
+ else if(iS-1e+9 == 10030){ // t
+ resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kTriton);
+ dedxExp = GetExpDeDx(t,6);
+ }
+ else if(iS-1e+9 == 20030){ // 3He
+ resolutionTPC = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[7])*5*0.07;
+ dedxExp = GetExpDeDx(t,7);
+ }
+ else if(iS-1e+9 == 20040){ // 4He
+ resolutionTPC = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[5])*5*0.07;
+ dedxExp = GetExpDeDx(t,8);
+ }
if(resolutionTPC > -1) dedx = fTPCResponseF->GetRandom()*resolutionTPC + dedxExp;
else dedx = 0;
}
else if(iS==5) resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kDeuteron);
else if(iS==6) resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kTriton);
else if(iS==7) resolutionTPC = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[7])*5*0.07;
+ else if(iS==8) resolutionTPC = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[8])*5*0.07;
if(centr < 0) resolutionTPC *= 0.78;
if(centr < 10) resolutionTPC *= 1.0;
Float_t timeTOF = t->GetTOFsignal() - fPIDesd->GetTOFResponse().GetStartTime(p);
// TOF mismatch weight
- Int_t det[5];
- Float_t length, timeextra, pos[3];
+ Float_t length,timeextra;
/* compute length and expected time */
- fgTofGeo->GetVolumeIndices(t->GetTOFCalChannel(), det);
- fgTofGeo->GetPosPar(det, pos);
- length = 0.;
- for (Int_t i = 0; i < 3; i++) length += pos[i] * pos[i];
- length = TMath::Sqrt(length);
+ length = fgHtofChannelDist->GetBinContent(t->GetTOFCalChannel()%8736);
timeextra = length * 33.3564095198152043;
Float_t dz =t->GetTOFsignalDz();
Float_t dx =t->GetTOFsignalDx();
Float_t mismweight = TMath::Max(fgMism->Eval(timeTOF - timeextra),0.0000001) * ((0.5 + 0.05/pt/pt/pt)*(0.75 + 0.23 * (1.3*dx*dx + 0.7*dz*dz))); // mismatch probabilities
fWTofMism = mismfrac*mismweight;
- Double_t inttimes[8];
+ Double_t inttimes[9];
t->GetIntegratedTimes(inttimes);
inttimes[5] = inttimes[0] / p * fMass[5] * TMath::Sqrt(1+p*p/fMass[5]/fMass[5]);
inttimes[6] = inttimes[0] / p * fMass[6] * TMath::Sqrt(1+p*p/fMass[6]/fMass[6]);
inttimes[7] = inttimes[0] / p * fMass[7] * TMath::Sqrt(1+p*p/fMass[7]/fMass[7]);
+ inttimes[8] = inttimes[0] / p * fMass[8] * TMath::Sqrt(1+p*p/fMass[8]/fMass[8]);
for(Int_t iS=0;iS<fgkNspecies;iS++){
Float_t expsigma = fPIDesd->GetTOFResponse().GetExpectedSigma(p, inttimes[iS], fMass[iS]);
resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kPion);
dedxExp = GetExpDeDx(t,2);
}
- else if(iS==321){
- resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kKaon);
- dedxExp = GetExpDeDx(t,3);
- }
- else if(iS==2212){
- resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kProton);
- dedxExp = GetExpDeDx(t,4);
- }
+ else if(iS==321){
+ resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kKaon);
+ dedxExp = GetExpDeDx(t,3);
+ }
+ else if(iS==2212){
+ resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kProton);
+ dedxExp = GetExpDeDx(t,4);
+ }
+ else if(iS-1e+9 == 10020){ // d
+ resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kDeuteron);
+ dedxExp = GetExpDeDx(t,5);
+ }
+ else if(iS-1e+9 == 10030){ // t
+ resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kTriton);
+ dedxExp = GetExpDeDx(t,6);
+ }
+ else if(iS-1e+9 == 20030){ // 3He
+ resolutionTPC = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[7])*5*0.07;
+ dedxExp = GetExpDeDx(t,7);
+ }
+ else if(iS-1e+9 == 20040){ // 4He
+ resolutionTPC = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[5])*5*0.07;
+ dedxExp = GetExpDeDx(t,8);
+ }
if(resolutionTPC > -1)
dedx = fTPCResponseF->GetRandom()*resolutionTPC + dedxExp;
else dedx = 0;
else if(iS==5) resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kDeuteron);
else if(iS==6) resolutionTPC = fPIDesd->GetTPCResponse().GetExpectedSigma(momtpc,t->GetTPCsignalN(),AliPID::kTriton);
else if(iS==7) resolutionTPC = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[7])*5*0.07;
+ else if(iS==8) resolutionTPC = fPIDesd->GetTPCResponse().Bethe(momtpc/fMass[8])*5*0.07;
if(centr < 0) resolutionTPC *= 0.78;
if(centr < 10) resolutionTPC *= 1.0;
Float_t timeTOF = t->GetTOFsignal() - fPIDesd->GetTOFResponse().GetStartTime(p);
// TOF mismatch weight
- Int_t det[5];
- Float_t length, timeextra, pos[3];
+ Float_t length, timeextra;
/* compute length and expected time */
Float_t etaAbs = TMath::Abs(t->Eta());
Int_t extrapolatedTOFchannel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
if(t->Eta() < 0) extrapolatedTOFchannel = 8736-1- extrapolatedTOFchannel;
extrapolatedTOFchannel = (Int_t(extrapolatedTOFchannel) % 8736);
-
- fgTofGeo->GetVolumeIndices(extrapolatedTOFchannel, det);
- fgTofGeo->GetPosPar(det, pos);
- length = 0.;
- for (Int_t i = 0; i < 3; i++) length += pos[i] * pos[i];
- length = TMath::Sqrt(length);
+ length = fgHtofChannelDist->GetBinContent(extrapolatedTOFchannel);
timeextra = length * 33.3564095198152043;
Float_t mismweight = TMath::Max(fgMism->Eval(timeTOF - timeextra),0.0000001) * ((0.5 + 0.05/pt/pt/pt)); // mismatch probabilities
fWTofMism = mismfrac*mismweight;
- Double_t inttimes[8];
+ Double_t inttimes[9];
t->GetIntegratedTimes(inttimes);
inttimes[5] = inttimes[0] / p * fMass[5] * TMath::Sqrt(1+p*p/fMass[5]/fMass[5]);
inttimes[6] = inttimes[0] / p * fMass[6] * TMath::Sqrt(1+p*p/fMass[6]/fMass[6]);
inttimes[7] = inttimes[0] / p * fMass[7] * TMath::Sqrt(1+p*p/fMass[7]/fMass[7]);
+ inttimes[8] = inttimes[0] / p * fMass[8] * TMath::Sqrt(1+p*p/fMass[8]/fMass[8]);
for(Int_t iS=0;iS<fgkNspecies;iS++){
if(iS==0){
fghPriors[5]->SetBinContent(k1,k2,(fC[icentr][TMath::Max(ipt-5,0)][4]*(1-weight) + 0.2*weight)/deutFact);
fghPriors[6]->SetBinContent(k1,k2,(fC[icentr][TMath::Max(ipt-5,0)][4]*(1-weight) + 0.2*weight)/deutFact/deutFact);
fghPriors[7]->SetBinContent(k1,k2,(fC[icentr][TMath::Max(ipt-5,0)][4]*(1-weight) + 0.2*weight)/deutFact/deutFact/deutFact);
- } // end loop on pt
+ fghPriors[8]->SetBinContent(k1,k2,(fC[icentr][TMath::Max(ipt-5,0)][4]*(1-weight) + 0.2*weight)/deutFact/deutFact/deutFact/deutFact);
+ } // end loop on pt
} // end loop on centrality bins
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff