LoadMyLibs();
TFile f("Output.root");
- //AliPerformanceRes * compObj = (AliPerformanceRes*)f.Get("AliPerformanceRes");
- AliPerformanceRes * compObj = (AliPerformanceRes*)cOutput->FindObject("AliPerformanceRes");
+ AliPerformanceRes * compObj = (AliPerformanceRes*)coutput->FindObject("AliPerformanceRes");
// analyse comparison data
compObj->Analyse();
#include "AliESDEvent.h"
#include "AliESDVertex.h"
#include "AliESDtrack.h"
+#include "AliESDfriendTrack.h"
+#include "AliESDfriend.h"
#include "AliLog.h"
#include "AliMCEvent.h"
#include "AliMCParticle.h"
binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
}
- /*
- Int_t nPtBins = 31;
- Double_t binsPt[32] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.25,2.5,2.75,3.,3.5,4.,5.,6.,8.,10.};
- Double_t ptMin = 0., ptMax = 10.;
-
- if(IsHptGenerator() == kTRUE) {
- nPtBins = 100;
- ptMin = 0.; ptMax = 100.;
- }
- */
-
Double_t yMin = -0.02, yMax = 0.02;
Double_t zMin = -12.0, zMax = 12.0;
if(GetAnalysisMode() == 3) { // TrackRef coordinate system
zMin = -100.; zMax = 100.;
}
- // res_y:res_z:res_phi,res_lambda:res_1pt:y:z:eta:phi:pt
- Int_t binsResolHisto[10]={100,100,100,100,100,25,50,30,144,nPtBins};
- Double_t minResolHisto[10]={-1.,-1.,-0.03,-0.03,-0.2, yMin, zMin,-1.5,0.,ptMin};
- Double_t maxResolHisto[10]={ 1., 1., 0.03, 0.03, 0.2, yMax, zMax, 1.5,2.*TMath::Pi(), ptMax};
+ // res_y:res_z:res_phi,res_lambda:res_pt:y:z:eta:phi:pt
+ Int_t binsResolHisto[10]={100,100,100,100,100,25,50,144,30,nPtBins};
+ Double_t minResolHisto[10]={-1.,-1.,-0.03,-0.03,-0.2, yMin, zMin, 0., -1.5, ptMin};
+ Double_t maxResolHisto[10]={ 1., 1., 0.03, 0.03, 0.2, yMax, zMax, 2.*TMath::Pi(), 1.5, ptMax};
- fResolHisto = new THnSparseF("fResolHisto","res_y:res_z:res_phi:res_lambda:res_1pt:y:z:eta:phi:pt",10,binsResolHisto,minResolHisto,maxResolHisto);
+ fResolHisto = new THnSparseF("fResolHisto","res_y:res_z:res_phi:res_lambda:res_pt:y:z:phi:eta:pt",10,binsResolHisto,minResolHisto,maxResolHisto);
fResolHisto->SetBinEdges(9,binsPt);
fResolHisto->GetAxis(0)->SetTitle("y-y_{mc} (cm)");
fResolHisto->GetAxis(1)->SetTitle("z-z_{mc} (cm)");
fResolHisto->GetAxis(2)->SetTitle("#phi-#phi_{mc} (rad)");
fResolHisto->GetAxis(3)->SetTitle("#lambda-#lambda_{mc} (rad)");
- fResolHisto->GetAxis(4)->SetTitle("(p_{Tmc}/p_{T}-1)");
+ fResolHisto->GetAxis(4)->SetTitle("(p_{T}/p_{Tmc}-1)");
fResolHisto->GetAxis(5)->SetTitle("y_{mc} (cm)");
fResolHisto->GetAxis(6)->SetTitle("z_{mc} (cm)");
- fResolHisto->GetAxis(7)->SetTitle("#eta_{mc}");
- fResolHisto->GetAxis(8)->SetTitle("#phi_{mc} (rad)");
+ fResolHisto->GetAxis(7)->SetTitle("#phi_{mc} (rad)");
+ fResolHisto->GetAxis(8)->SetTitle("#eta_{mc}");
fResolHisto->GetAxis(9)->SetTitle("p_{Tmc} (GeV/c)");
fResolHisto->Sumw2();
- //pull_y:pull_z:pull_lambda:pull_phi:pull_1pt:y:z:eta:phi:pt
- Int_t binsPullHisto[10]={100,100,100,100,100,50,50,30,144,nPtBins};
- Double_t minPullHisto[10]={-5.,-5.,-5.,-5.,-5.,yMin, zMin,-1.5, 0., ptMin};
- Double_t maxPullHisto[10]={ 5., 5., 5., 5., 5., yMax, zMax, 1.5, 2.*TMath::Pi(),ptMax};
+ ////pull_y:pull_z:pull_phi:pull_lambda:pull_1pt:y:z:eta:phi:pt
+ //Int_t binsPullHisto[10]={100,100,100,100,100,50,50,30,144,nPtBins};
+ //Double_t minPullHisto[10]={-5.,-5.,-5.,-5.,-5.,yMin, zMin,-1.5, 0., ptMin};
+ //Double_t maxPullHisto[10]={ 5., 5., 5., 5., 5., yMax, zMax, 1.5, 2.*TMath::Pi(),ptMax};
+ //fPullHisto = new THnSparseF("fPullHisto","pull_y:pull_z:pull_phi:pull_lambda:pull_1pt:y:z:eta:phi:pt",10,binsPullHisto,minPullHisto,maxPullHisto);
- fPullHisto = new THnSparseF("fPullHisto","pull_y:pull_z:pull_phi:pull_lambda:pull_1pt:y:z:eta:phi:pt",10,binsPullHisto,minPullHisto,maxPullHisto);
- if(!IsHptGenerator()) fPullHisto->SetBinEdges(9,binsPt);
+ //pull_y:pull_z:pull_snp:pull_tgl:pull_1pt:y:z:snp:tgl:1pt
+ Int_t binsPullHisto[10]={100,100,100,100,100,50,50,50,50,nPtBins};
+ Double_t minPullHisto[10]={-5.,-5.,-5.,-5.,-5.,yMin, zMin,-1., -2.0, ptMin};
+ Double_t maxPullHisto[10]={ 5., 5., 5., 5., 5., yMax, zMax, 1., 2.0, ptMax};
+ fPullHisto = new THnSparseF("fPullHisto","pull_y:pull_z:pull_y:pull_z:pull_snp:pull_tgl:pull_1pt:y:z:snp:tgl:1pt",10,binsPullHisto,minPullHisto,maxPullHisto);
+ /*
+ if(!IsHptGenerator()) fPullHisto->SetBinEdges(9,bins1Pt);
fPullHisto->GetAxis(0)->SetTitle("(y-y_{mc})/#sigma");
fPullHisto->GetAxis(1)->SetTitle("(z-z_{mc})/#sigma");
fPullHisto->GetAxis(2)->SetTitle("(#phi-#phi_{mc})/#sigma");
fPullHisto->GetAxis(8)->SetTitle("#phi_{mc} (rad)");
fPullHisto->GetAxis(9)->SetTitle("p_{Tmc} (GeV/c)");
fPullHisto->Sumw2();
+ */
+
+ fPullHisto->GetAxis(0)->SetTitle("(y-y_{mc})/#sigma");
+ fPullHisto->GetAxis(1)->SetTitle("(z-z_{mc})/#sigma");
+ fPullHisto->GetAxis(2)->SetTitle("(sin#phi-sin#phi_{mc})/#sigma");
+ fPullHisto->GetAxis(3)->SetTitle("(tan#lambda-tan#lambda_{mc})/#sigma");
+ fPullHisto->GetAxis(4)->SetTitle("(p_{Tmc}/p_{T}-1)/#sigma");
+ fPullHisto->GetAxis(5)->SetTitle("y_{mc} (cm)");
+ fPullHisto->GetAxis(6)->SetTitle("z_{mc} (cm)");
+ fPullHisto->GetAxis(7)->SetTitle("sin#phi_{mc}");
+ fPullHisto->GetAxis(8)->SetTitle("tan#lambda_{mc}");
+ fPullHisto->GetAxis(9)->SetTitle("1/p_{Tmc} (GeV/c)^{-1}");
+ fPullHisto->Sumw2();
// Init cuts
if(!fCutsMC)
//
if(!stack) return;
Int_t label = TMath::Abs(esdTrack->GetLabel());
- if(label == 0) return;
-
TParticle* particle = stack->Particle(label);
if(!particle) return;
if(!particle->GetPDG()) return;
if(particle->GetPDG()->Charge()==0) return;
//printf("charge %d \n",particle->GetPDG()->Charge());
+
+ // Only 5 charged particle species (e,mu,pi,K,p)
+ if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) return;
+
+ // exclude electrons
+ if (fCutsMC->GetEP()==TMath::Abs(particle->GetPdgCode())) return;
- Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
+ Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
Float_t mceta = particle->Eta();
Float_t mcphi = particle->Phi();
if(mcphi<0) mcphi += 2.*TMath::Pi();
Float_t mcpt = particle->Pt();
+ Float_t mcsnp = TMath::Sin(TMath::ATan2(particle->Py(),particle->Px()));
+ Float_t mctgl = TMath::Tan(TMath::ATan2(particle->Pz(),particle->Pt()));
// nb. TPC clusters cut
if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
// select primaries
if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
{
+ if(mcpt == 0) return;
+
deltaYTPC= track->GetY()-particle->Vy();
deltaZTPC = track->GetZ()-particle->Vz();
deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(particle->Py(),particle->Px());
- delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
+ //delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
+ deltaPtTPC = (track->Pt()-mcpt) / mcpt;
pullYTPC= (track->GetY()-particle->Vy()) / TMath::Sqrt(track->GetSigmaY2());
pullZTPC = (track->GetZ()-particle->Vz()) / TMath::Sqrt(track->GetSigmaZ2());
- Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
- Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
- pullLambdaTPC = deltaLambdaTPC / sigma_lambda;
- pullPhiTPC = deltaPhiTPC / sigma_phi;
+ //Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
+ //Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
+ pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2());
+ pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
//pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
//pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
else pull1PtTPC = 0.;
- Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
+ Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,particle->Vy(),particle->Vz(),mcphi,mceta,mcpt};
fResolHisto->Fill(vResolHisto);
- Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
+ Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mcsnp,mctgl,1./mcpt};
fPullHisto->Fill(vPullHisto);
}
}
if(!stack) return;
Int_t label = TMath::Abs(esdTrack->GetLabel());
- if(label == 0) return;
-
TParticle* particle = stack->Particle(label);
if(!particle) return;
if(particle->GetPDG()->Charge()==0) return;
//printf("charge %d \n",particle->GetPDG()->Charge());
+
+ // Only 5 charged particle species (e,mu,pi,K,p)
+ if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) return;
+
+ // exclude electrons
+ if (fCutsMC->GetEP()==TMath::Abs(particle->GetPdgCode())) return;
+
Float_t mceta = particle->Eta();
Float_t mcphi = particle->Phi();
if(mcphi<0) mcphi += 2.*TMath::Pi();
Float_t mcpt = particle->Pt();
+ Float_t mcsnp = TMath::Sin(TMath::ATan2(particle->Py(),particle->Px()));
+ Float_t mctgl = TMath::Tan(TMath::ATan2(particle->Pz(),particle->Pt()));
if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
Int_t clusterITS[200];
if(esdTrack->GetITSclusters(clusterITS)<fCutsRC->GetMinNClustersITS()) return; // min. nb. ITS clusters
- Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
+ Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
// select primaries
if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
{
+ if(mcpt == 0) return;
+
+ deltaYTPC= esdTrack->GetY()-particle->Vy();
+ deltaZTPC = esdTrack->GetZ()-particle->Vz();
+ deltaLambdaTPC = TMath::ATan2(esdTrack->Pz(),esdTrack->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
+ deltaPhiTPC = TMath::ATan2(esdTrack->Py(),esdTrack->Px())-TMath::ATan2(particle->Py(),particle->Px());
+ //delta1PtTPC = (esdTrack->OneOverPt()-1./mcpt)*mcpt;
+ deltaPtTPC = (esdTrack->Pt()-mcpt) / mcpt;
+
+ pullYTPC= (esdTrack->GetY()-particle->Vy()) / TMath::Sqrt(esdTrack->GetSigmaY2());
+ pullZTPC = (esdTrack->GetZ()-particle->Vz()) / TMath::Sqrt(esdTrack->GetSigmaZ2());
+
+ //Double_t sigma_lambda = 1./(1.+esdTrack->GetTgl()*esdTrack->GetTgl()) * TMath::Sqrt(esdTrack->GetSigmaTgl2());
+ //Double_t sigma_phi = 1./TMath::Sqrt(1-esdTrack->GetSnp()*esdTrack->GetSnp()) * TMath::Sqrt(esdTrack->GetSigmaSnp2());
+ pullPhiTPC = (esdTrack->GetSnp() - mcsnp) / TMath::Sqrt(esdTrack->GetSigmaSnp2());
+ pullLambdaTPC = (esdTrack->GetTgl() - mctgl) / TMath::Sqrt(esdTrack->GetSigmaTgl2());
+
+ //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(esdTrack->GetSigmaTgl2());
+ //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(esdTrack->GetSigmaSnp2());
+ if (mcpt) pull1PtTPC = (esdTrack->OneOverPt()-1./mcpt) / TMath::Sqrt(esdTrack->GetSigma1Pt2());
+ else pull1PtTPC = 0.;
+
+ Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,particle->Vy(),particle->Vz(),mcphi,mceta,mcpt};
+ fResolHisto->Fill(vResolHisto);
+
+ Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mcsnp,mctgl,1./mcpt};
+ fPullHisto->Fill(vPullHisto);
+
+
+ /*
deltaYTPC= esdTrack->GetY()-particle->Vy();
deltaZTPC = esdTrack->GetZ()-particle->Vz();
deltaLambdaTPC = TMath::ATan2(esdTrack->Pz(),esdTrack->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
fPullHisto->Fill(vPullHisto);
+ */
}
}
if(!stack) return;
Int_t label = TMath::Abs(esdTrack->GetLabel());
- if(label == 0) return;
-
TParticle* particle = stack->Particle(label);
if(!particle) return;
if(particle->GetPDG()->Charge()==0) return;
//printf("charge %d \n",particle->GetPDG()->Charge());
+ // Only 5 charged particle species (e,mu,pi,K,p)
+ if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) return;
+
+ // exclude electrons
+ if (fCutsMC->GetEP()==TMath::Abs(particle->GetPdgCode())) return;
+
Float_t mceta = particle->Eta();
Float_t mcphi = particle->Phi();
if(mcphi<0) mcphi += 2.*TMath::Pi();
Float_t mcpt = particle->Pt();
+ Float_t mcsnp = TMath::Sin(TMath::ATan2(particle->Py(),particle->Px()));
+ Float_t mctgl = TMath::Tan(TMath::ATan2(particle->Pz(),particle->Pt()));
if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
- Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
+ Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
// select primaries
if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
{
+
+ if(mcpt == 0) return;
+
+ deltaYTPC= track->GetY()-particle->Vy();
+ deltaZTPC = track->GetZ()-particle->Vz();
+ deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
+ deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(particle->Py(),particle->Px());
+ //delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
+ deltaPtTPC = (track->Pt()-mcpt) / mcpt;
+
+ pullYTPC= (track->GetY()-particle->Vy()) / TMath::Sqrt(track->GetSigmaY2());
+ pullZTPC = (track->GetZ()-particle->Vz()) / TMath::Sqrt(track->GetSigmaZ2());
+
+ //Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
+ //Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
+ pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2());
+ pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
+
+ //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
+ //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
+ if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
+ else pull1PtTPC = 0.;
+
+ Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,particle->Vy(),particle->Vz(),mcphi,mceta,mcpt};
+ fResolHisto->Fill(vResolHisto);
+
+ Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mcsnp,mctgl,1./mcpt};
+ fPullHisto->Fill(vPullHisto);
+
+ /*
+
deltaYTPC= track->GetY()-particle->Vy();
deltaZTPC = track->GetZ()-particle->Vz();
deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
fPullHisto->Fill(vPullHisto);
+
+ */
}
}
if(!track) return;
Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
- esdTrack->GetImpactParameters(dca,cov);
+ esdTrack->GetImpactParametersTPC(dca,cov);
//
// Fill rec vs MC information
if(!mcEvent) return;
Int_t label = TMath::Abs(esdTrack->GetLabel());
- if(label == 0) return;
AliMCParticle *mcParticle = mcEvent->GetTrack(label);
if(!mcParticle) return;
AliTrackReference *ref0 = GetFirstTPCTrackRef(mcParticle);
if(!ref0) return;
+ // Only 5 charged particle species (e,mu,pi,K,p)
+ TParticle *particle = mcParticle->Particle();
+ if(!particle) return;
+
+ if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) return;
+
+ // exclude electrons
+ if (fCutsMC->GetEP()==TMath::Abs(particle->GetPdgCode())) return;
+
// calculate alpha angle
Double_t xyz[3] = {ref0->X(),ref0->Y(),ref0->Z()};
Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
Float_t mcphi = ref0->Phi();
if(mcphi<0) mcphi += 2.*TMath::Pi();
Float_t mcpt = ref0->Pt();
+ Float_t mcsnp = TMath::Sin(TMath::ATan2(ref0->Py(),ref0->Px()));
+ Float_t mctgl = TMath::Tan(TMath::ATan2(ref0->Pz(),ref0->Pt()));
if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
- Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
+ Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
// select primaries
if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
{
- //deltaYTPC= track->GetY()-ref0->Y();
- deltaYTPC= track->GetY(); // it corresponds to deltaY after alpha rotation
+ if(mcpt == 0) return;
+
+ deltaYTPC= track->GetY(); // already rotated
deltaZTPC = track->GetZ()-ref0->Z();
deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(ref0->Pz(),ref0->Pt());
deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(ref0->Py(),ref0->Px());
- if(mcpt) delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
- else delta1PtTPC = 0.;
+ //delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
+ deltaPtTPC = (track->Pt()-mcpt) / mcpt;
- //pullYTPC= (track->GetY()-ref0->Y()) / TMath::Sqrt(track->GetSigmaY2());
- pullYTPC= track->GetY() / TMath::Sqrt(track->GetSigmaY2()); //it corresponds to deltaY after alpha rotation
+ pullYTPC= track->GetY() / TMath::Sqrt(track->GetSigmaY2());
pullZTPC = (track->GetZ()-ref0->Z()) / TMath::Sqrt(track->GetSigmaZ2());
+
+ //Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
+ //Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
+ pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2());
+ pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
//pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
//pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
+ if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
+ else pull1PtTPC = 0.;
- Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
- Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
- pullLambdaTPC = deltaLambdaTPC / sigma_lambda;
- pullPhiTPC = deltaPhiTPC / sigma_phi;
+ Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,ref0->Y(),ref0->Z(),mcphi,mceta,mcpt};
+ fResolHisto->Fill(vResolHisto);
+
+ Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,ref0->Y(),ref0->Z(),mcsnp,mctgl,1./mcpt};
+ fPullHisto->Fill(vPullHisto);
+ }
+
+ if(track) delete track;
+}
+
+//_____________________________________________________________________________
+void AliPerformanceRes::ProcessOuterTPC(AliMCEvent *const mcEvent, AliESDtrack *const esdTrack, AliESDfriendTrack *const friendTrack)
+{
+ //
+ // Fill resolution comparison information (outer params at TPC reference point)
+ //
+ if(!friendTrack) return;
+
+ const AliExternalTrackParam * outerParam = friendTrack->GetTPCOut();
+ if(!outerParam) return;
+
+ // create new AliExternalTrackParam
+ AliExternalTrackParam *track = new AliExternalTrackParam(*outerParam);
+ if(!track) return;
+
+ Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
+ esdTrack->GetImpactParametersTPC(dca,cov);
+
+ //
+ // Fill rec vs MC information
+ //
+ if(!mcEvent) return;
+
+ Int_t label = TMath::Abs(esdTrack->GetLabel());
+ AliMCParticle *mcParticle = mcEvent->GetTrack(label);
+ if(!mcParticle) return;
+
+ // get the last TPC track reference
+ AliTrackReference *ref0 = GetLastTPCTrackRef(mcParticle);
+ if(!ref0) return;
+
+ // Only 5 charged particle species (e,mu,pi,K,p)
+ TParticle *particle = mcParticle->Particle();
+ if(!particle) return;
+
+ if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) return;
+
+ // exclude electrons
+ if (fCutsMC->GetEP()==TMath::Abs(particle->GetPdgCode())) return;
+
+ // calculate alpha angle
+ Double_t xyz[3] = {ref0->X(),ref0->Y(),ref0->Z()};
+ Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
+ //if (alpha<0) alpha += TMath::TwoPi();
- if(mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
+ // rotate outer track to local coordinate system
+ // and propagate to the radius of the last track reference of TPC
+ Double_t trRadius = TMath::Sqrt(xyz[1] * xyz[1] + xyz[0] * xyz[0]);
+ Bool_t isOK = track->Propagate(alpha,trRadius,AliTracker::GetBz());
+ if(!isOK) return;
+
+ Float_t mceta = -TMath::Log(TMath::Tan(0.5 * ref0->Theta()));
+ Float_t mcphi = ref0->Phi();
+ if(mcphi<0) mcphi += 2.*TMath::Pi();
+ Float_t mcpt = ref0->Pt();
+ Float_t mcsnp = TMath::Sin(TMath::ATan2(ref0->Py(),ref0->Px()));
+ Float_t mctgl = TMath::Tan(TMath::ATan2(ref0->Pz(),ref0->Pt()));
+
+ if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
+ if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
+
+ Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
+ Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
+
+ // select primaries
+ if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
+ {
+ if(mcpt == 0) return;
+
+ deltaYTPC= track->GetY(); // already rotated
+ deltaZTPC = track->GetZ()-ref0->Z();
+ deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(ref0->Pz(),ref0->Pt());
+ deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(ref0->Py(),ref0->Px());
+ //delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
+ deltaPtTPC = (track->Pt()-mcpt) / mcpt;
+
+ pullYTPC= track->GetY() / TMath::Sqrt(track->GetSigmaY2());
+ pullZTPC = (track->GetZ()-ref0->Z()) / TMath::Sqrt(track->GetSigmaZ2());
+
+ //Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
+ //Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
+ pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2());
+ pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
+
+ //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
+ //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
+ if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
else pull1PtTPC = 0.;
- Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,ref0->Y(),ref0->Z(),mceta,mcphi,mcpt};
+ Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,ref0->Y(),ref0->Z(),mcphi,mceta,mcpt};
fResolHisto->Fill(vResolHisto);
- Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,ref0->Y(),ref0->Z(),mceta,mcphi,mcpt};
+ Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,ref0->Y(),ref0->Z(),mcsnp,mctgl,1./mcpt};
fPullHisto->Fill(vPullHisto);
}
// return first TPC track reference
if(!mcParticle) return 0;
+ // find first track reference
+ // check direction to select proper reference point for looping tracks
Int_t nTrackRef = mcParticle->GetNumberOfTrackReferences();
- AliTrackReference *ref0 = 0;
+ AliTrackReference *ref = 0;
+ AliTrackReference *refIn = 0;
for (Int_t iref = 0; iref < nTrackRef; iref++) {
- ref0 = mcParticle->GetTrackReference(iref);
- if(ref0 && (ref0->DetectorId()==AliTrackReference::kTPC))
+ ref = mcParticle->GetTrackReference(iref);
+ if(ref && (ref->DetectorId()==AliTrackReference::kTPC))
+ {
+ Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
+ if(dir < 0.) break;
+
+ refIn = ref;
break;
+ }
}
-return ref0;
+return refIn;
}
//_____________________________________________________________________________
-void AliPerformanceRes::Exec(AliMCEvent* const mcEvent, AliESDEvent* const esdEvent, const Bool_t bUseMC)
+AliTrackReference * AliPerformanceRes::GetLastTPCTrackRef(AliMCParticle *mcParticle)
+{
+ // return last TPC track reference
+ if(!mcParticle) return 0;
+
+ // find last track reference
+ // check direction to select proper reference point for looping tracks
+ Int_t nTrackRef = mcParticle->GetNumberOfTrackReferences();
+ AliTrackReference *ref = 0;
+ AliTrackReference *refOut = 0;
+ for (Int_t iref = 0; iref < nTrackRef; iref++) {
+ ref = mcParticle->GetTrackReference(iref);
+ if(ref && (ref->DetectorId()==AliTrackReference::kTPC)) {
+ Float_t dir=ref->X()*ref->Px()+ref->Y()*ref->Py();
+ if(dir< 0.0) break;
+ refOut = ref;
+ }
+ }
+
+return refOut;
+}
+
+//_____________________________________________________________________________
+void AliPerformanceRes::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
{
// Process comparison information
//
genHeader->PrimaryVertex(vtxMC);
} // end bUseMC
+
+ // use ESD friends
+ if(bUseESDfriend) {
+ if(!esdFriend) {
+ AliDebug(AliLog::kError, "esdFriend not available");
+ return;
+ }
+ }
// Process events
for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
{
- AliESDtrack*track = esdEvent->GetTrack(iTrack);
+ AliESDtrack *track = esdEvent->GetTrack(iTrack);
if(!track) continue;
+ AliESDfriendTrack *friendTrack=0;
+ if(bUseESDfriend) {
+ friendTrack=esdFriend->GetTrack(iTrack);
+ if(!friendTrack) continue;
+ }
+
if(GetAnalysisMode() == 0) ProcessTPC(stack,track);
else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track);
else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
else if(GetAnalysisMode() == 3) ProcessInnerTPC(mcEvent,track);
+ else if(GetAnalysisMode() == 4) ProcessOuterTPC(mcEvent,track,friendTrack);
else {
printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
return;
}
//_____________________________________________________________________________
-void AliPerformanceRes::Analyse(){
+void AliPerformanceRes::Analyse() {
// Analyse comparison information and store output histograms
// in the folder "folderRes"
//
{
for(Int_t j=5; j<10; j++)
{
- if(j!=7) fResolHisto->GetAxis(7)->SetRangeUser(-0.9,0.9); // eta window
+ if(j!=8) fResolHisto->GetAxis(8)->SetRangeUser(-0.9,0.89); // eta window
fResolHisto->GetAxis(9)->SetRangeUser(0.16,10.); // pt threshold
+ //if(j!=8) fResolHisto->GetAxis(8)->SetRangeUser(0.,0.9); // eta window
+ //fResolHisto->GetAxis(9)->SetRangeUser(0.16,3.); // pt window
if(GetAnalysisMode() == 3) fResolHisto->GetAxis(5)->SetRangeUser(-80.,80.); // y range
h2D = (TH2F*)fResolHisto->Projection(i,j);
h = AliPerformanceRes::MakeResol(h2D,1,0,100);
sprintf(name,"h_res_%d_vs_%d",i,j);
-
- /*
- if(i<2) {
- h->SetMinimum(0.);
- h->SetMaximum(1.);
- } else {
- h->SetMinimum(0.);
- h->SetMaximum(0.2);
- }
- */
h->SetName(name);
h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
h = AliPerformanceRes::MakeResol(h2D,1,1,100);
//h = (TH1F*)arr->At(1);
sprintf(name,"h_mean_res_%d_vs_%d",i,j);
- /*
- h->SetMinimum(-0.05);
- h->SetMaximum(0.05);
- if(i<2) {
- h->SetMinimum(-0.05);
- h->SetMaximum(0.05);
- } else {
- h->SetMinimum(-0.02);
- h->SetMaximum(0.02);
- }
- */
h->SetName(name);
h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
if(j==9) h->SetBit(TH1::kLogX);
aFolderObj->Add(h);
- fResolHisto->GetAxis(7)->SetRangeUser(-1.5,1.5);
+ fResolHisto->GetAxis(8)->SetRangeUser(-1.5,1.5);
fResolHisto->GetAxis(9)->SetRangeUser(0.0,10.);
//
- if(j!=7) fPullHisto->GetAxis(7)->SetRangeUser(-0.9,0.9); // eta window
- fPullHisto->GetAxis(9)->SetRangeUser(0.16,10.); // pt threshold
+ if(j!=8) fPullHisto->GetAxis(8)->SetRangeUser(-1.0,0.99); // theta window
+ else fPullHisto->GetAxis(8)->SetRangeUser(-1.5,1.5); // theta window
+ fPullHisto->GetAxis(9)->SetRangeUser(0.0,10.); // pt threshold
if(GetAnalysisMode() == 3) fPullHisto->GetAxis(5)->SetRangeUser(-80.,80.); // y range
h2D = (TH2F*)fPullHisto->Projection(i,j);
h = AliPerformanceRes::MakeResol(h2D,1,0,100);
sprintf(name,"h_pull_%d_vs_%d",i,j);
- //h->SetMinimum(0.);
- //h->SetMaximum(2.);
h->SetName(name);
h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
sprintf(title,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
h->SetTitle(title);
- if(j==9) h->SetBit(TH1::kLogX);
+ //if(j==9) h->SetBit(TH1::kLogX);
aFolderObj->Add(h);
h = AliPerformanceRes::MakeResol(h2D,1,1,100);
sprintf(name,"h_mean_pull_%d_vs_%d",i,j);
- //h->SetMinimum(-0.2);
- //h->SetMaximum(0.2);
h->SetName(name);
h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
sprintf(title,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
h->SetTitle(title);
- if(j==9) h->SetBit(TH1::kLogX);
+ //if(j==9) h->SetBit(TH1::kLogX);
aFolderObj->Add(h);
-
- fPullHisto->GetAxis(7)->SetRangeUser(-1.5,1.5);
- fPullHisto->GetAxis(9)->SetRangeUser(0.0,10.);
}
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff