1 /**************************************************************************
2 * Author: Panos Christakoglou. *
3 * Contributors are mentioned in the code where appropriate. *
5 * Permission to use, copy, modify and distribute this software and its *
6 * documentation strictly for non-commercial purposes is hereby granted *
7 * without fee, provided that the above copyright notice appears in all *
8 * copies and that both the copyright notice and this permission notice *
9 * appear in the supporting documentation. The authors make no claims *
10 * about the suitability of this software for any purpose. It is *
11 * provided "as is" without express or implied warranty. *
12 **************************************************************************/
14 /* $Id: AliProtonAnalysisBase.cxx 31056 2009-02-16 14:31:41Z pchrist $ */
16 //-----------------------------------------------------------------
17 // AliProtonAnalysisBase class
18 // This is the class to deal with the proton analysis
19 // Origin: Panos Christakoglou | Panos.Christakoglou@cern.ch
20 //-----------------------------------------------------------------
21 #include <Riostream.h>
28 #include <AliExternalTrackParam.h>
29 #include <AliESDEvent.h>
31 #include <AliVertexerTracks.h>
32 #include <AliTPCpidESD.h>
36 #include "AliProtonAnalysisBase.h"
38 ClassImp(AliProtonAnalysisBase)
40 //____________________________________________________________________//
41 AliProtonAnalysisBase::AliProtonAnalysisBase() :
42 TObject(), fProtonAnalysisLevel("ESD"),
43 fTriggerMode(kMB2), fProtonAnalysisMode(kTPC), fProtonPIDMode(kBayesian),
44 fAnalysisEtaMode(kFALSE),
45 fVxMax(100.), fVyMax(100.), fVzMax(100.),
46 fNBinsX(0), fMinX(0), fMaxX(0),
47 fNBinsY(0), fMinY(0), fMaxY(0),
48 fMinTPCClusters(0), fMinITSClusters(0),
49 fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
50 fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
51 fMaxSigmaToVertex(0), fMaxSigmaToVertexTPC(0),
52 fMaxDCAXY(0), fMaxDCAXYTPC(0),
53 fMaxDCAZ(0), fMaxDCAZTPC(0),
54 fMaxDCA3D(0), fMaxDCA3DTPC(0),
55 fMaxConstrainChi2(0), fMinTPCdEdxPoints(0),
56 fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
57 fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
58 fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE),
59 fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
60 fMaxSigmaToVertexFlag(kFALSE), fMaxSigmaToVertexTPCFlag(kFALSE),
61 fMaxDCAXYFlag(kFALSE), fMaxDCAXYTPCFlag(kFALSE),
62 fMaxDCAZFlag(kFALSE), fMaxDCAZTPCFlag(kFALSE),
63 fMaxDCA3DFlag(kFALSE), fMaxDCA3DTPCFlag(kFALSE),
64 fMaxConstrainChi2Flag(kFALSE),
65 fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
66 fESDpidFlag(kFALSE), fTPCpidFlag(kFALSE), fTOFpidFlag(kFALSE),
67 fPointOnITSLayer1Flag(0), fPointOnITSLayer2Flag(0),
68 fPointOnITSLayer3Flag(0), fPointOnITSLayer4Flag(0),
69 fPointOnITSLayer5Flag(0), fPointOnITSLayer6Flag(0),
70 fMinTPCdEdxPointsFlag(kFALSE),
71 fFunctionProbabilityFlag(kFALSE),
73 fElectronFunction(0), fMuonFunction(0),
74 fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
77 for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
78 for(Int_t i = 0; i < 24; i++) {
82 fListVertexQA = new TList();
83 fListVertexQA->SetName("fListVertexQA");
84 TH1F *gHistVx = new TH1F("gHistVx",
85 "Vx distribution;V_{x} [cm];Entries",
87 gHistVx->SetFillColor(kRed-2);
88 fListVertexQA->Add(gHistVx);
89 TH1F *gHistVxAccepted = new TH1F("gHistVxaccepted",
90 "Vx distribution;V_{x} [cm];Entries",
92 fListVertexQA->Add(gHistVxAccepted);
93 TH1F *gHistVy = new TH1F("gHistVy",
94 "Vy distribution;V_{y} [cm];Entries",
96 gHistVy->SetFillColor(kRed-2);
97 fListVertexQA->Add(gHistVy);
98 TH1F *gHistVyAccepted = new TH1F("gHistVyaccepted",
99 "Vy distribution;V_{y} [cm];Entries",
101 fListVertexQA->Add(gHistVyAccepted);
102 TH1F *gHistVz = new TH1F("gHistVz",
103 "Vz distribution;V_{z} [cm];Entries",
105 gHistVz->SetFillColor(kRed-2);
106 fListVertexQA->Add(gHistVz);
107 TH1F *gHistVzAccepted = new TH1F("gHistVzaccepted",
108 "Vz distribution;V_{z} [cm];Entries",
110 fListVertexQA->Add(gHistVzAccepted);
115 //____________________________________________________________________//
116 AliProtonAnalysisBase::~AliProtonAnalysisBase() {
118 if(fElectronFunction) delete fElectronFunction;
119 if(fMuonFunction) delete fMuonFunction;
120 if(fPionFunction) delete fPionFunction;
121 if(fKaonFunction) delete fKaonFunction;
122 if(fProtonFunction) delete fProtonFunction;
123 if(fListVertexQA) delete fListVertexQA;
126 //____________________________________________________________________//
127 Double_t AliProtonAnalysisBase::GetParticleFraction(Int_t i, Double_t p) {
128 //Return the a priori probs
130 if(fFunctionProbabilityFlag) {
131 if(i == 0) partFrac = fElectronFunction->Eval(p);
132 if(i == 1) partFrac = fMuonFunction->Eval(p);
133 if(i == 2) partFrac = fPionFunction->Eval(p);
134 if(i == 3) partFrac = fKaonFunction->Eval(p);
135 if(i == 4) partFrac = fProtonFunction->Eval(p);
137 else partFrac = fPartFrac[i];
142 //____________________________________________________________________//
143 Bool_t AliProtonAnalysisBase::IsInPhaseSpace(AliESDtrack* const track) {
144 // Checks if the track is outside the analyzed y-Pt phase space
145 Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
148 if((fProtonAnalysisMode == kTPC) || (fProtonAnalysisMode == kHybrid)) {
149 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
151 gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0; eta = -10.0;
154 gPt = tpcTrack->Pt();
155 gPx = tpcTrack->Px();
156 gPy = tpcTrack->Py();
157 gPz = tpcTrack->Pz();
158 eta = tpcTrack->Eta();
160 }//standalone TPC or Hybrid TPC approaches
169 if((gPt < fMinY) || (gPt > fMaxY)) {
171 Printf("IsInPhaseSpace: Track rejected because it has a Pt value of %lf (accepted interval: %lf - %lf)",gPt,fMinY,fMaxY);
174 if(fAnalysisEtaMode) {
175 if((eta < fMinX) || (eta > fMaxX)) {
177 Printf("IsInPhaseSpace: Track rejected because it has an eta value of %lf (accepted interval: %lf - %lf)",eta,fMinX,fMaxX);
182 if((Rapidity(gPx,gPy,gPz) < fMinX) || (Rapidity(gPx,gPy,gPz) > fMaxX)) {
184 Printf("IsInPhaseSpace: Track rejected because it has a y value of %lf (accepted interval: %lf - %lf)",Rapidity(gPx,gPy,gPz),fMinX,fMaxX);
192 //____________________________________________________________________//
193 Bool_t AliProtonAnalysisBase::IsAccepted(AliESDEvent *esd,
194 const AliESDVertex *vertex,
195 AliESDtrack* track) {
196 // Checks if the track is excluded from the cuts
197 Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
198 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0}; //The impact parameters and their covariance.
199 Double_t dca3D = 0.0;
201 if((fProtonAnalysisMode == kTPC)||(fProtonAnalysisMode == kHybrid)) {
202 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
204 gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
205 dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
206 cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
209 gPt = tpcTrack->Pt();
210 gPx = tpcTrack->Px();
211 gPy = tpcTrack->Py();
212 gPz = tpcTrack->Pz();
213 tpcTrack->PropagateToDCA(vertex,
214 esd->GetMagneticField(),
217 }//standalone TPC or hybrid TPC approaches
223 track->PropagateToDCA(vertex,
224 esd->GetMagneticField(),
227 dca3D = TMath::Sqrt(TMath::Power(dca[0],2) +
228 TMath::Power(dca[1],2));
231 Int_t nClustersITS = track->GetITSclusters(fIdxInt);
232 Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
234 Float_t chi2PerClusterITS = -1;
236 chi2PerClusterITS = track->GetITSchi2()/Float_t(nClustersITS);
237 Float_t chi2PerClusterTPC = -1;
239 chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
242 track->GetExternalCovariance(extCov);
244 if(fPointOnITSLayer1Flag) {
245 if(!track->HasPointOnITSLayer(0)) {
247 Printf("IsAccepted: Track rejected because it doesn't have a point on the 1st ITS layer");
251 if(fPointOnITSLayer2Flag) {
252 if(!track->HasPointOnITSLayer(1)) {
254 Printf("IsAccepted: Track rejected because it doesn't have a point on the 2nd ITS layer");
258 if(fPointOnITSLayer3Flag) {
259 if(!track->HasPointOnITSLayer(2)) {
261 Printf("IsAccepted: Track rejected because it doesn't have a point on the 3rd ITS layer");
265 if(fPointOnITSLayer4Flag) {
266 if(!track->HasPointOnITSLayer(3)) {
268 Printf("IsAccepted: Track rejected because it doesn't have a point on the 4th ITS layer");
272 if(fPointOnITSLayer5Flag) {
273 if(!track->HasPointOnITSLayer(4)) {
275 Printf("IsAccepted: Track rejected because it doesn't have a point on the 5th ITS layer");
279 if(fPointOnITSLayer6Flag) {
280 if(!track->HasPointOnITSLayer(5)) {
282 Printf("IsAccepted: Track rejected because it doesn't have a point on the 6th ITS layer");
286 if(fMinITSClustersFlag) {
287 if(nClustersITS < fMinITSClusters) {
289 Printf("IsAccepted: Track rejected because it has %d ITS points (min. requested: %d)",nClustersITS,fMinITSClusters);
293 if(fMaxChi2PerITSClusterFlag) {
294 if(chi2PerClusterITS > fMaxChi2PerITSCluster) {
296 Printf("IsAccepted: Track rejected because it has a chi2 per ITS cluster %lf (max. requested: %lf)",chi2PerClusterITS,fMaxChi2PerITSCluster);
300 if(fMinTPCClustersFlag) {
301 if(nClustersTPC < fMinTPCClusters) {
303 Printf("IsAccepted: Track rejected because it has %d TPC clusters (min. requested: %d)",nClustersTPC,fMinTPCClusters);
307 if(fMaxChi2PerTPCClusterFlag) {
308 if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) {
310 Printf("IsAccepted: Track rejected because it has a chi2 per TPC cluster %lf (max. requested: %lf)",chi2PerClusterTPC,fMaxChi2PerTPCCluster);
315 if(extCov[0] > fMaxCov11) {
317 Printf("IsAccepted: Track rejected because it has a cov11 value of %lf (max. requested: %lf)",extCov[0],fMaxCov11);
322 if(extCov[2] > fMaxCov22) {
324 Printf("IsAccepted: Track rejected because it has a cov22 value of %lf (max. requested: %lf)",extCov[2],fMaxCov22);
329 if(extCov[5] > fMaxCov33) {
331 Printf("IsAccepted: Track rejected because it has a cov33 value of %lf (max. requested: %lf)",extCov[5],fMaxCov33);
336 if(extCov[9] > fMaxCov44) {
338 Printf("IsAccepted: Track rejected because it has a cov44 value of %lf (max. requested: %lf)",extCov[9],fMaxCov44);
343 if(extCov[14] > fMaxCov55) {
345 Printf("IsAccepted: Track rejected because it has a cov55 value of %lf (max. requested: %lf)",extCov[14],fMaxCov55);
349 if(fMaxSigmaToVertexFlag) {
350 if(GetSigmaToVertex(track) > fMaxSigmaToVertex) {
352 Printf("IsAccepted: Track rejected because it has a %lf sigmas to vertex (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertex);
356 if(fMaxSigmaToVertexTPCFlag) {
357 if(GetSigmaToVertex(track) > fMaxSigmaToVertexTPC) {
359 Printf("IsAccepted: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertexTPC);
364 if(TMath::Abs(dca[0]) > fMaxDCAXY) {
366 Printf("IsAccepted: Track rejected because it has a value of dca(xy) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXY);
370 if(fMaxDCAXYTPCFlag) {
371 if(TMath::Abs(dca[0]) > fMaxDCAXYTPC) {
373 Printf("IsAccepted: Track rejected because it has a value of dca(xy) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXYTPC);
378 if(TMath::Abs(dca[1]) > fMaxDCAZ) {
380 Printf("IsAccepted: Track rejected because it has a value of dca(z) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZ);
384 if(fMaxDCAZTPCFlag) {
385 if(TMath::Abs(dca[1]) > fMaxDCAZTPC) {
387 Printf("IsAccepted: Track rejected because it has a value of dca(z) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZTPC);
392 if(TMath::Abs(dca3D) > fMaxDCA3D) {
394 Printf("IsAccepted: Track rejected because it has a value of dca(3D) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3D);
398 if(fMaxDCA3DTPCFlag) {
399 if(TMath::Abs(dca3D) > fMaxDCA3DTPC) {
401 Printf("IsAccepted: Track rejected because it has a value of dca(3D) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3DTPC);
405 if(fMaxConstrainChi2Flag) {
406 if(track->GetConstrainedChi2() > 0)
407 if(TMath::Log(track->GetConstrainedChi2()) > fMaxConstrainChi2) {
409 Printf("IsAccepted: Track rejected because it has a value of the constrained chi2 to the vertex of %lf (max. requested: %lf)",TMath::Log(track->GetConstrainedChi2()),fMaxConstrainChi2);
413 if(fMinTPCdEdxPointsFlag) {
414 if(track->GetTPCsignalN() < fMinTPCdEdxPoints) {
416 Printf("IsAccepted: Track rejected because it has %d TPC points for the calculation of the energy loss (min. requested: %d)",track->GetTPCsignalN(),fMinTPCdEdxPoints);
421 if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) {
423 Printf("IsAccepted: Track rejected because it has no ITS refit flag");
428 if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) {
430 Printf("IsAccepted: Track rejected because it has no TPC refit flag");
435 if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) {
437 Printf("IsAccepted: Track rejected because it has no ESD pid flag");
442 if ((track->GetStatus() & AliESDtrack::kTPCpid) == 0) {
444 Printf("IsAccepted: Track rejected because it has no TPC pid flag");
449 if ((track->GetStatus() & AliESDtrack::kTOFpid) == 0) {
451 Printf("IsAccepted: Track rejected because it has no TOF pid flag");
459 //____________________________________________________________________//
460 Float_t AliProtonAnalysisBase::GetSigmaToVertex(AliESDtrack* esdTrack) const {
461 // Calculates the number of sigma to the vertex.
466 if((fProtonAnalysisMode == kTPC)&&(fProtonAnalysisMode != kHybrid))
467 esdTrack->GetImpactParametersTPC(b,bCov);
469 esdTrack->GetImpactParameters(b,bCov);
471 if (bCov[0]<=0 || bCov[2]<=0) {
472 //AliDebug(1, "Estimated b resolution lower or equal zero!");
473 bCov[0]=0; bCov[2]=0;
475 bRes[0] = TMath::Sqrt(bCov[0]);
476 bRes[1] = TMath::Sqrt(bCov[2]);
478 if (bRes[0] == 0 || bRes[1] ==0) return -1;
480 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
482 if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
484 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
489 //____________________________________________________________________//
490 Double_t AliProtonAnalysisBase::Rapidity(Double_t gPx,
492 Double_t gPz) const {
493 //returns the rapidity of the proton - to be removed
494 Double_t fMass = 9.38270000000000048e-01;
496 Double_t gP = TMath::Sqrt(TMath::Power(gPx,2) +
497 TMath::Power(gPy,2) +
498 TMath::Power(gPz,2));
499 Double_t energy = TMath::Sqrt(gP*gP + fMass*fMass);
502 y = 0.5*TMath::Log((energy + gPz)/(energy - gPz));
507 //________________________________________________________________________
508 const AliESDVertex* AliProtonAnalysisBase::GetVertex(AliESDEvent* esd,
513 // Get the vertex from the ESD and returns it if the vertex is valid
514 // Second argument decides which vertex is used (this selects
515 // also the quality criteria that are applied)
516 const AliESDVertex* vertex = 0;
518 vertex = esd->GetPrimaryVertexSPD();
519 else if(mode == kTPC){
520 Double_t kBz = esd->GetMagneticField();
521 AliVertexerTracks vertexer(kBz);
522 vertexer.SetTPCMode();
523 AliESDVertex *vTPC = vertexer.FindPrimaryVertex(esd);
524 esd->SetPrimaryVertexTPC(vTPC);
525 for (Int_t i=0; i<esd->GetNumberOfTracks(); i++) {
526 AliESDtrack *t = esd->GetTrack(i);
527 t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
530 vertex = esd->GetPrimaryVertexTPC();
532 else if(mode == kGlobal)
533 vertex = esd->GetPrimaryVertex();
535 Printf("GetVertex: ERROR: Invalid second argument %d", mode);
539 Printf("GetVertex: Event rejected because there is no valid vertex object");
543 // check Ncontributors
544 if(vertex->GetNContributors() <= 0) {
546 Printf("GetVertex: Event rejected because the number of contributors for the vertex determination is <= 0");
551 Double_t zRes = vertex->GetZRes();
554 Printf("GetVertex: Event rejected because the value of the vertex resolution in z is 0");
557 ((TH1F *)(fListVertexQA->At(0)))->Fill(vertex->GetXv());
558 ((TH1F *)(fListVertexQA->At(2)))->Fill(vertex->GetYv());
559 ((TH1F *)(fListVertexQA->At(4)))->Fill(vertex->GetZv());
561 if(TMath::Abs(vertex->GetXv()) > gVxMax) {
563 Printf("GetVertex: Event rejected because it has a Vx value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetXv()),gVxMax,gVxMax);
566 if(TMath::Abs(vertex->GetYv()) > gVyMax) {
568 Printf("GetVertex: Event rejected because it has a Vy value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetYv()),gVyMax,gVyMax);
571 if(TMath::Abs(vertex->GetZv()) > gVzMax) {
573 Printf("GetVertex: Event rejected because it has a Vz value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetZv()),gVzMax,gVzMax);
576 ((TH1F *)(fListVertexQA->At(1)))->Fill(vertex->GetXv());
577 ((TH1F *)(fListVertexQA->At(3)))->Fill(vertex->GetYv());
578 ((TH1F *)(fListVertexQA->At(5)))->Fill(vertex->GetZv());
583 //________________________________________________________________________
584 Bool_t AliProtonAnalysisBase::IsEventTriggered(const AliESDEvent *esd,
585 TriggerMode trigger) {
586 // check if the event was triggered
587 ULong64_t triggerMask = esd->GetTriggerMask();
589 // definitions from p-p.cfg
590 ULong64_t spdFO = (1 << 14);
591 ULong64_t v0left = (1 << 11);
592 ULong64_t v0right = (1 << 12);
596 if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
601 if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
606 if (triggerMask & spdFO)
615 //________________________________________________________________________
616 TCanvas *AliProtonAnalysisBase::GetListOfCuts() {
617 // return the list of cuts and their cut values for reference
624 TCanvas *c = new TCanvas("cListOfCuts","List of cuts",0,0,900,600);
625 c->SetFillColor(10); c->SetHighLightColor(41);
628 c->cd(1)->SetFillColor(10);
629 l.DrawLatex(0.3,0.9,"Analysis details: List of cuts\n\n");
631 listOfCuts = "Analysis level: "; listOfCuts += fProtonAnalysisLevel;
632 l.DrawLatex(0.1,0.82,listOfCuts.Data());
633 listOfCuts = "Analysis mode: ";
634 if(fProtonAnalysisMode == kTPC) listOfCuts += "TPC standalone";
635 if(fProtonAnalysisMode == kHybrid) listOfCuts += "Hybrid TPC";
636 if(fProtonAnalysisMode == kGlobal) listOfCuts += "Global tracking";
637 l.DrawLatex(0.1,0.74,listOfCuts.Data());
638 listOfCuts = "Trigger mode: ";
639 if(fTriggerMode == kMB1) listOfCuts += "Minimum bias 1";
640 if(fTriggerMode == kMB2) listOfCuts += "Minimum bias 2";
641 if(fTriggerMode == kSPDFASTOR) listOfCuts += "FastOR (SPD)";
642 l.DrawLatex(0.1,0.66,listOfCuts.Data());
643 listOfCuts = "PID mode: ";
644 if(fProtonPIDMode == kBayesian) listOfCuts += "Bayesian PID";
645 if(fProtonPIDMode == kRatio) listOfCuts += "Z = ln((dE/dx)_{exp.}/(dE/dx)_{theor.})";
646 if(fProtonPIDMode == kSigma1) {
647 listOfCuts += "N_{#sigma}(1) area: "; listOfCuts += fNSigma;
648 listOfCuts += " #sigma";
650 if(fProtonPIDMode == kSigma2) {
651 listOfCuts += "N_{#sigma}(2) area: "; listOfCuts += fNSigma;
652 listOfCuts += " #sigma";
654 l.DrawLatex(0.1,0.58,listOfCuts.Data());
655 listOfCuts = "Accepted vertex diamond: ";
656 l.DrawLatex(0.1,0.5,listOfCuts.Data());
657 listOfCuts = "|V_{x}| < "; listOfCuts += fVxMax; listOfCuts += " [cm]";
658 l.DrawLatex(0.6,0.5,listOfCuts.Data());
659 listOfCuts = "|V_{y}| < "; listOfCuts += fVyMax; listOfCuts += " [cm]";
660 l.DrawLatex(0.6,0.4,listOfCuts.Data());
661 listOfCuts = "|V_{z}| < "; listOfCuts += fVzMax; listOfCuts += " [cm]";
662 l.DrawLatex(0.6,0.3,listOfCuts.Data());
663 listOfCuts = "Phase space: ";
664 l.DrawLatex(0.1,0.2,listOfCuts.Data());
665 if(fAnalysisEtaMode) listOfCuts = "|#eta| < ";
666 else listOfCuts = "|y| < ";
667 listOfCuts += TMath::Abs(fMinX);
668 l.DrawLatex(0.35,0.2,listOfCuts.Data());
669 listOfCuts = "N_{bins} = ";
670 listOfCuts += fNBinsX; listOfCuts += " (binning: ";
671 listOfCuts += (fMaxX-fMinX)/fNBinsX; listOfCuts += ")";
672 l.DrawLatex(0.35,0.15,listOfCuts.Data());
673 listOfCuts = ""; listOfCuts += fMinY; listOfCuts += " < P_{T} < ";
674 listOfCuts += fMaxY; listOfCuts += "GeV/c";
675 l.DrawLatex(0.35,0.1,listOfCuts.Data());
676 listOfCuts = "N_{bins} = ";
677 listOfCuts += fNBinsY; listOfCuts += " (binning: ";
678 listOfCuts += (fMaxY-fMinY)/fNBinsY; listOfCuts += ")";
679 l.DrawLatex(0.35,0.05,listOfCuts.Data());
681 c->cd(2)->SetFillColor(2);
682 l.DrawLatex(0.3,0.9,"ITS related cuts");
683 listOfCuts = "Request a cluster on SPD1: ";
684 listOfCuts += fPointOnITSLayer1Flag;
685 l.DrawLatex(0.1,0.8,listOfCuts.Data());
686 listOfCuts = "Request a cluster on SPD2: ";
687 listOfCuts += fPointOnITSLayer2Flag;
688 l.DrawLatex(0.1,0.7,listOfCuts.Data());
689 listOfCuts = "Request a cluster on SDD1: ";
690 listOfCuts += fPointOnITSLayer3Flag;
691 l.DrawLatex(0.1,0.6,listOfCuts.Data());
692 listOfCuts = "Request a cluster on SDD2: ";
693 listOfCuts += fPointOnITSLayer4Flag;
694 l.DrawLatex(0.1,0.5,listOfCuts.Data());
695 listOfCuts = "Request a cluster on SSD1: ";
696 listOfCuts += fPointOnITSLayer5Flag;
697 l.DrawLatex(0.1,0.4,listOfCuts.Data());
698 listOfCuts = "Request a cluster on SSD2: ";
699 listOfCuts += fPointOnITSLayer6Flag;
700 l.DrawLatex(0.1,0.3,listOfCuts.Data());
701 listOfCuts = "Minimum number of ITS clusters: ";
702 if(fMinITSClustersFlag) listOfCuts += fMinITSClusters;
703 else listOfCuts += "Not used";
704 l.DrawLatex(0.1,0.2,listOfCuts.Data());
705 listOfCuts = "Maximum #chi^{2} per ITS cluster: ";
706 if(fMaxChi2PerITSClusterFlag) listOfCuts += fMaxChi2PerITSCluster;
707 else listOfCuts += "Not used";
708 l.DrawLatex(0.1,0.1,listOfCuts.Data());
710 c->cd(3)->SetFillColor(3);
711 l.DrawLatex(0.3,0.9,"TPC related cuts");
712 listOfCuts = "Minimum number of TPC clusters: ";
713 if(fMinTPCClustersFlag) listOfCuts += fMinTPCClusters;
714 else listOfCuts += "Not used";
715 l.DrawLatex(0.1,0.7,listOfCuts.Data());
716 listOfCuts = "Maximum #chi^{2} per TPC cluster: ";
717 if(fMaxChi2PerTPCClusterFlag) listOfCuts += fMaxChi2PerTPCCluster;
718 else listOfCuts += "Not used";
719 l.DrawLatex(0.1,0.5,listOfCuts.Data());
720 listOfCuts = "Minimum number of TPC points for the dE/dx: ";
721 if(fMinTPCdEdxPointsFlag) listOfCuts += fMinTPCdEdxPoints;
722 else listOfCuts += "Not used";
723 l.DrawLatex(0.1,0.3,listOfCuts.Data());
725 c->cd(4)->SetFillColor(4);
726 l.DrawLatex(0.3,0.9,"Tracking related cuts");
727 listOfCuts = "Maximum cov11: ";
728 if(fMaxCov11Flag) listOfCuts += fMaxCov11;
729 else listOfCuts += "Not used";
730 l.DrawLatex(0.1,0.75,listOfCuts.Data());
731 listOfCuts = "Maximum cov22: ";
732 if(fMaxCov22Flag) listOfCuts += fMaxCov22;
733 else listOfCuts += "Not used";
734 l.DrawLatex(0.1,0.6,listOfCuts.Data());
735 listOfCuts = "Maximum cov33: ";
736 if(fMaxCov33Flag) listOfCuts += fMaxCov33;
737 else listOfCuts += "Not used";
738 l.DrawLatex(0.1,0.45,listOfCuts.Data());
739 listOfCuts = "Maximum cov44: ";
740 if(fMaxCov44Flag) listOfCuts += fMaxCov44;
741 else listOfCuts += "Not used";
742 l.DrawLatex(0.1,0.3,listOfCuts.Data());
743 listOfCuts = "Maximum cov55: ";
744 if(fMaxCov55Flag) listOfCuts += fMaxCov55;
745 else listOfCuts += "Not used";
746 l.DrawLatex(0.1,0.15,listOfCuts.Data());
748 c->cd(5)->SetFillColor(5);
749 l.DrawLatex(0.3,0.9,"DCA related cuts");
750 listOfCuts = "Maximum sigma to vertex: ";
751 if(fMaxSigmaToVertexFlag) listOfCuts += fMaxSigmaToVertex;
752 else listOfCuts += "Not used";
753 l.DrawLatex(0.1,0.8,listOfCuts.Data());
754 listOfCuts = "Maximum sigma to vertex (TPC): ";
755 if(fMaxSigmaToVertexTPCFlag) listOfCuts += fMaxSigmaToVertexTPC;
756 else listOfCuts += "Not used";
757 l.DrawLatex(0.1,0.72,listOfCuts.Data());
758 listOfCuts = "Maximum DCA in xy: ";
759 if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXY;
760 else listOfCuts += "Not used";
761 l.DrawLatex(0.1,0.64,listOfCuts.Data());
762 listOfCuts = "Maximum DCA in z: ";
763 if(fMaxDCAZFlag) listOfCuts += fMaxDCAZ;
764 else listOfCuts += "Not used";
765 l.DrawLatex(0.1,0.56,listOfCuts.Data());
766 listOfCuts = "Maximum DCA in 3D: ";
767 if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3D;
768 else listOfCuts += "Not used";
769 l.DrawLatex(0.1,0.48,listOfCuts.Data());
770 listOfCuts = "Maximum DCA in xy (TPC): ";
771 if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXYTPC;
772 else listOfCuts += "Not used";
773 l.DrawLatex(0.1,0.4,listOfCuts.Data());
774 listOfCuts = "Maximum DCA in z (TPC): ";
775 if(fMaxDCAZFlag) listOfCuts += fMaxDCAZTPC;
776 else listOfCuts += "Not used";
777 l.DrawLatex(0.1,0.32,listOfCuts.Data());
778 listOfCuts = "Maximum DCA in 3D (TPC): ";
779 if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3DTPC;
780 else listOfCuts += "Not used";
781 l.DrawLatex(0.1,0.24,listOfCuts.Data());
782 listOfCuts = "Maximum constrained #chi^{2} (vertex): ";
783 if(fMaxConstrainChi2Flag) listOfCuts += fMaxConstrainChi2;
784 else listOfCuts += "Not used";
785 l.DrawLatex(0.1,0.16,listOfCuts.Data());
787 c->cd(6)->SetFillColor(6);
788 l.DrawLatex(0.3,0.9,"Tracking bits related cuts");
789 listOfCuts = "Request the ITS refit bit: "; listOfCuts += fITSRefitFlag;
790 l.DrawLatex(0.1,0.7,listOfCuts.Data());
791 listOfCuts = "Request the TPC refit bit: "; listOfCuts += fTPCRefitFlag;
792 l.DrawLatex(0.1,0.5,listOfCuts.Data());
793 listOfCuts = "Request the TPC pid bit: "; listOfCuts += fTPCpidFlag;
794 l.DrawLatex(0.1,0.3,listOfCuts.Data());
795 listOfCuts = "Request the ESD pid bit: "; listOfCuts += fESDpidFlag;
796 l.DrawLatex(0.1,0.1,listOfCuts.Data());
801 //________________________________________________________________________
802 Bool_t AliProtonAnalysisBase::IsProton(AliESDtrack *track) {
803 //Function that checks if a track is a proton
804 Double_t probability[5];
805 Double_t gPt = 0.0, gP = 0.0, gEta = 0.0;
806 Long64_t fParticleType = 0;
808 //Bayesian approach for the PID
809 if(fProtonPIDMode == kBayesian) {
810 if((fProtonAnalysisMode == kTPC)||(fProtonAnalysisMode == kHybrid)) {
811 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
813 gPt = tpcTrack->Pt();
815 track->GetTPCpid(probability);
817 }//TPC standalone or Hybrid TPC
818 else if(fProtonAnalysisMode == kGlobal) {
821 track->GetESDpid(probability);
825 for(Int_t i = 0; i < AliPID::kSPECIES; i++)
826 rcc += probability[i]*GetParticleFraction(i,gP);
829 for(Int_t i = 0; i < AliPID::kSPECIES; i++)
830 w[i] = probability[i]*GetParticleFraction(i,gP)/rcc;
831 fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
833 if(fParticleType == 4)
836 //Ratio of the measured over the theoretical dE/dx a la STAR
837 else if(fProtonPIDMode == kRatio) {
838 Printf("The kRatio mode is not implemented yet!!!");
841 //Definition of an N-sigma area around the dE/dx vs P band
842 else if(fProtonPIDMode == kSigma1) {
843 Double_t fAlephParameters[5];
844 fAlephParameters[0] = 4.23232575531564326e+00; //50*0.76176e-1;
845 fAlephParameters[1] = 8.68482806165147636e+00; //10.632;
846 fAlephParameters[2] = 1.34000000000000005e-05; //0.13279e-4;
847 fAlephParameters[3] = 2.30445734159456084e+00; //1.8631;
848 fAlephParameters[4] = 2.25624744086878559e+00; //1.9479;
850 AliTPCpidESD *fPidObject = new AliTPCpidESD();
851 fPidObject->SetBetheBlochParameters(fAlephParameters[0]/50.,fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
853 Double_t nsigma = TMath::Abs(fPidObject->GetNumberOfSigmas(track,AliPID::kProton));
854 if(nsigma <= fNSigma)
856 }//kSigma1 PID method
857 //Another definition of an N-sigma area around the dE/dx vs P band
858 else if(fProtonPIDMode == kSigma2) {
859 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
861 gPt = tpcTrack->Pt();
863 gEta = tpcTrack->Eta();
865 //We start the P slices at >0.3GeV/c with a bining of 50MeV/c ==> Int_t(0.3001/0.05) = 6
866 Int_t nbinP = Int_t(gP/0.05) - 6;
867 Double_t tpcSignal = track->GetTPCsignal();
868 Double_t dEdxTheory = Bethe(gP/9.38270000000000048e-01);
869 Double_t dEdxSigma = fdEdxSigma[nbinP];
870 Double_t nsigma = TMath::Abs(tpcSignal - dEdxTheory)/(tpcSignal*(dEdxSigma/TMath::Sqrt(track->GetTPCsignalN())));
871 if(nsigma <= fNSigma)
878 //________________________________________________________________________
879 void AliProtonAnalysisBase::SetdEdxBandInfo(const char *filename) {
880 // This function is used in case the kSigma1 or kSigma2 PID mode is selected
881 // It takes as an argument the name of the ascii file (for the time being)
882 // that is generated as a prior process.
883 // This ascii file has three columns: The min. P value (bins of 50MeV/c)
884 // the mean and the sigma of the dE/dx distributions for protons coming
885 // from a gaussian fit.
889 Double_t gPtMin = 0.0;
892 in >> gPtMin >> fdEdxMean[iCounter] >> fdEdxSigma[iCounter];
894 Printf("Momentum bin: %d - Min momentum: %lf - mean(dE/dx): %lf - sigma(dE/dx): %lf",iCounter+1,gPtMin,fdEdxMean[iCounter],fdEdxSigma[iCounter]);
899 //________________________________________________________________________
900 Double_t AliProtonAnalysisBase::Bethe(Double_t bg) const {
901 // This is the Bethe-Bloch function normalised to 1 at the minimum
902 // We renormalize it based on the MC information
903 // WARNING: To be revised soon!!!
904 // This is just a temporary fix
905 Double_t normalization = 49.2;
907 Double_t beta2 = bg2/(1.+ bg2);
908 Double_t bb = 8.62702e-2*(9.14550 - beta2 - TMath::Log(3.51000e-5 + 1./bg2))/beta2;
910 const Float_t kmeanCorrection =0.1;
911 Double_t meanCorrection =(1+(bb-1)*kmeanCorrection);
912 bb *= meanCorrection;
914 return normalization*bb;