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, proviyaded 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 purapose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 //////////////////////////////////////////////////////////////////////////////
18 // Analysis task for the systematic study of the uncertainties related to //
19 // the tracking and ITS-TPC matching efficiency for different particle //
22 //////////////////////////////////////////////////////////////////////////////
25 #include "Riostream.h"
32 #include "TParticlePDG.h"
34 #include "AliAnalysisTaskSE.h"
35 #include "AliAnalysisManager.h"
37 #include "AliESDtrackCuts.h"
38 #include "AliESDVertex.h"
39 #include "AliESDEvent.h"
40 #include "AliESDInputHandler.h"
41 #include "AliESDtrack.h"
42 #include "AliESDpid.h"
43 #include "AliESDUtils.h"
44 #include "AliMCEventHandler.h"
45 #include "AliMCEvent.h"
49 #include "AliAnalysisTrackingUncertainties.h"
52 ClassImp(AliAnalysisTrackingUncertainties)
54 // histogram constants
55 const Int_t kNumberOfAxes = 5;
57 //________________________________________________________________________
58 AliAnalysisTrackingUncertainties::AliAnalysisTrackingUncertainties()
59 : AliAnalysisTaskSE("TaskTestPA"),
69 // default Constructor
70 /* fast compilation test
72 gSystem->Load("libANALYSIS");
73 gSystem->Load("libANALYSISalice");
74 .L AliAnalysisTrackingUncertainties.cxx++
79 //________________________________________________________________________
80 AliAnalysisTrackingUncertainties::AliAnalysisTrackingUncertainties(const char *name)
81 : AliAnalysisTaskSE(name),
92 // standard constructur which should be used
98 fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
99 fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
100 fESDtrackCuts->SetEtaRange(-1., 1.);
102 // analysis utils if needed
104 fUtils = new AliAnalysisUtils();
107 // Output slot #0 writes into a TList container
108 DefineOutput(1, TList::Class());
113 //________________________________________________________________________
114 void AliAnalysisTrackingUncertainties::UserCreateOutputObjects()
120 fListHist = new TList();
121 fListHist->SetOwner(kTRUE);
123 // (1.) basic QA and statistics histograms
125 TH2F * histVertexSelection = new TH2F("histVertexSelection", "vertex selection; vertex z (cm); accepted/rejected", 100, -50., 50., 2, -0.5, 1.5);
126 fListHist->Add(histVertexSelection);
128 // (2.) track cut variation histograms
130 InitializeTrackCutHistograms();
132 // (3.) ITS -> TPC matching histograms
134 Int_t binsMatch[kNumberOfAxes] = { 10, 50, 20, 18, 6};
135 Double_t minMatch[kNumberOfAxes] = { 0, 0.1, -1, 0, -0.5};
136 Double_t maxMatch[kNumberOfAxes] = {200, 20, +1, 2*TMath::Pi(), 5.5};
138 TString axisNameMatch[kNumberOfAxes] = {"matchChi2","pT","eta","phi","pid"};
139 TString axisTitleMatch[kNumberOfAxes] = {"matchChi2","pT","eta","phi","pid"};
141 THnF * hBestMatch = new THnF("hBestMatch","ITS -> TPC matching ",kNumberOfAxes, binsMatch, minMatch, maxMatch);
142 for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
143 hBestMatch->GetAxis(iaxis)->SetName(axisNameMatch[iaxis]);
144 hBestMatch->GetAxis(iaxis)->SetTitle(axisTitleMatch[iaxis]);
146 BinLogAxis(hBestMatch, 1);
147 fListHist->Add(hBestMatch);
152 const double ptMax=5;
154 TH2F * hNMatch = new TH2F("hNMatch","N Matches",nbPt,0,ptMax,kMaxMatch+1,-0.5,kMaxMatch+0.5);
155 TH2F * hAllMatch = new TH2F("hAllMatch","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
156 TH2F * hAllMatchGlo = new TH2F("hAllMatchGlo","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
157 fListHist->Add(hNMatch);
158 fListHist->Add(hAllMatch);
159 fListHist->Add(hAllMatchGlo);
161 TH2F * hNMatchBg = new TH2F("hNMatchBg","N Matches",nbPt,0,ptMax,kMaxMatch+1,-0.5,kMaxMatch+0.5);
162 TH2F * hBestMatchBg = new TH2F("hBestMatchBg","Best Match Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
163 TH2F * hAllMatchBg = new TH2F("hAllMatchBg","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
164 TH2F * hAllMatchGloBg = new TH2F("hAllMatchGloBg","All Matches Chi2",nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2);
165 fListHist->Add(hNMatchBg);
166 fListHist->Add(hBestMatchBg);
167 fListHist->Add(hAllMatchBg);
168 fListHist->Add(hAllMatchGloBg);
172 PostData(1, fListHist);
179 //________________________________________________________________________
180 void AliAnalysisTrackingUncertainties::UserExec(Option_t *)
185 fESD = dynamic_cast<AliESDEvent*>( InputEvent() );
187 PostData(1, fListHist);
191 if (!fESDpid) fESDpid = ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid();
193 // Check Monte Carlo information and other access first:
195 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
200 // extract generated particles information
202 AliMCEvent* mcEvent = 0x0;
203 AliStack* stack = 0x0;
204 if (eventHandler) mcEvent = eventHandler->MCEvent();
210 stack = mcEvent->Stack();
213 for(Int_t i = 0; i < stack->GetNtrack(); i++) {
214 /* at the moment nothing is needed here
215 TParticle * trackMC = stack->Particle(i);
216 Double_t rap = trackMC->Eta();
217 Double_t y = trackMC->Y();
218 Double_t pT = trackMC->Pt();
225 if (!fESDtrackCuts) {
226 PostData(1, fListHist);
230 // monitor vertex position and selection
232 TH2F * histVertexSelection = (TH2F *) fListHist->FindObject("histVertexSelection");
234 Float_t vertexZ = 0.;
235 if (IsVertexAccepted(fESD, vertexZ)) {
236 histVertexSelection->Fill(vertexZ, 0);
238 histVertexSelection->Fill(vertexZ, 1);
242 // fill track cut variation histograms
244 ProcessTrackCutVariation();
246 // fill ITS->TPC matching histograms
248 ProcessItsTpcMatching();
251 PostData(1, fListHist);
256 //________________________________________________________________________
257 void AliAnalysisTrackingUncertainties::ProcessItsTpcMatching(){
259 // check how many its-sa tracks get matched to TPC
261 int ntr = fESD->GetNumberOfTracks();
263 // initialize histograms
265 TH2F * hNMatch = (TH2F*) fListHist->FindObject("hNMatch");
266 THnF * hBestMatch = (THnF*) fListHist->FindObject("hBestMatch");
267 TH2F * hAllMatch = (TH2F*) fListHist->FindObject("hAllMatch");
268 TH2F * hAllMatchGlo = (TH2F*) fListHist->FindObject("hAllMatchGlo");
270 TH2F * hNMatchBg = (TH2F*) fListHist->FindObject("hNMatchBg");
271 TH2F * hBestMatchBg = (TH2F*) fListHist->FindObject("hBestMatchBg");
272 TH2F * hAllMatchBg = (TH2F*) fListHist->FindObject("hAllMatchBg");
273 TH2F * hAllMatchGloBg = (TH2F*) fListHist->FindObject("hAllMatchGloBg");
275 for (int it=0;it<ntr;it++) {
276 AliESDtrack* trSA = fESD->GetTrack(it);
277 if (!trSA->IsOn(AliESDtrack::kITSpureSA) || !trSA->IsOn(AliESDtrack::kITSrefit)) continue;
278 double pt = trSA->Pt();
281 for (int i=kMaxMatch;i--;) {fMatchChi[i]=0; fMatchTr[i]=0;} // reset array
282 for (int it1=0;it1<ntr;it1++) {
283 if (it1==it) continue;
284 AliESDtrack* trESD = fESD->GetTrack(it1);
285 if (!trESD->IsOn(AliESDtrack::kTPCrefit)) continue;
286 Match(trSA,trESD, nmatch);
289 hNMatch->Fill(pt,nmatch);
290 if (nmatch>0) { // matched tracks
291 Double_t vecHistMatch[kNumberOfAxes] = {fMatchChi[0], pt, trSA->Eta(), trSA->Phi(), 0};
292 hBestMatch->Fill(vecHistMatch);
293 } else { // un-matched tracks -> should be in overflow bin
294 Double_t vecHistMatch[kNumberOfAxes] = {9999999., pt, trSA->Eta(), trSA->Phi(), 0};
295 hBestMatch->Fill(vecHistMatch);
298 for (int imt=nmatch;imt--;) {
299 hAllMatch->Fill(pt,fMatchChi[imt]);
300 if (fMatchTr[imt]->IsOn(AliESDtrack::kITSrefit)) hAllMatchGlo->Fill(pt,fMatchChi[imt]);
304 for (int i=kMaxMatch;i--;) {fMatchChi[i]=0; fMatchTr[i]=0;}
305 for (int it1=0;it1<ntr;it1++) {
306 if (it1==it) continue;
307 AliESDtrack* trESD = fESD->GetTrack(it1);
308 if (!trESD->IsOn(AliESDtrack::kTPCrefit)) continue;
309 Match(trSA,trESD, nmatch, TMath::Pi());
312 hNMatchBg->Fill(pt,nmatch);
313 if (nmatch>0) hBestMatchBg->Fill(pt,fMatchChi[0]);
314 for (int imt=nmatch;imt--;) {
315 hAllMatchBg->Fill(pt,fMatchChi[imt]);
316 if (fMatchTr[imt]->IsOn(AliESDtrack::kITSrefit)) hAllMatchGloBg->Fill(pt,fMatchChi[imt]);
325 void AliAnalysisTrackingUncertainties::Match(const AliESDtrack* tr0, const AliESDtrack* tr1, Int_t &nmatch, Double_t rotate) {
327 // check if two tracks are matching, possible rotation for combinatoric backgr.
329 Float_t bField = fESD->GetMagneticField();
331 const AliExternalTrackParam* trtpc0 = tr1->GetInnerParam();
333 AliExternalTrackParam trtpc(*trtpc0);
335 if (TMath::Abs(rotate)>1e-5) {
336 const double *par = trtpc.GetParameter();
337 const double *cov = trtpc.GetCovariance();
338 double alp = trtpc.GetAlpha() + rotate;
339 trtpc.Set(trtpc.GetX(),alp,par,cov);
342 if (!trtpc.Rotate(tr0->GetAlpha())) return;
343 if (!trtpc.PropagateTo(tr0->GetX(),bField)) return;
344 double chi2 = tr0->GetPredictedChi2(&trtpc);
345 if (chi2>kMaxChi2) return;
348 for (ins=0;ins<nmatch;ins++) if (chi2<fMatchChi[ins]) break;
349 if (ins>=kMaxMatch) return;
351 for (int imv=nmatch;imv>ins;imv--) {
352 if (imv>=kMaxMatch) continue;
353 fMatchTr[imv] = fMatchTr[imv-1];
354 fMatchChi[imv] = fMatchChi[imv-1];
357 fMatchChi[ins] = chi2;
359 if (nmatch>=kMaxMatch) nmatch = kMaxMatch;
364 //________________________________________________________________________
365 void AliAnalysisTrackingUncertainties::ProcessTrackCutVariation() {
367 // fill track cut variation histograms - undo cuts step-by-step and fill histograms
370 // initialize histograms
372 THnF * histNcl = (THnF *) fListHist->FindObject("histNcl");
373 THnF * histChi2Tpc = (THnF *) fListHist->FindObject("histChi2Tpc");
374 THnF * histDcaZ = (THnF *) fListHist->FindObject("histDcaZ");
375 THnF * histSpd = (THnF *) fListHist->FindObject("histSpd");
376 THnF * histNcr = (THnF *) fListHist->FindObject("histNcr");
377 THnF * histCRoverFC = (THnF *) fListHist->FindObject("histCRoverFC");
378 THnF * histChi2Its = (THnF *) fListHist->FindObject("histChi2Its");
379 THnF * histTpcLength = (THnF *) fListHist->FindObject("histTpcLength");
380 THnF * histTpcItsMatch = (THnF *) fListHist->FindObject("histTpcItsMatch");
382 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
384 for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
386 AliESDtrack *track =fESD->GetTrack(i);
388 // relevant variables
390 Double_t pid = Double_t(GetPid(track));
392 Int_t nclsTPC = track->GetTPCncls();
393 Float_t pT = track->Pt();
394 Float_t eta = track->Eta();
395 Float_t phi = track->Phi();
396 Float_t chi2TPC = track->GetTPCchi2();
397 Float_t ncrTPC = track->GetTPCCrossedRows();
398 Int_t nclsTPCF = track->GetTPCNclsF();
399 Float_t nCRoverFC = track->GetTPCCrossedRows();
400 Double_t chi2ITS = track->GetITSchi2();
401 Int_t nclsITS = track->GetITSclusters(0);
402 Float_t tpcLength = 0.;
404 if (track->GetInnerParam() && track->GetESDEvent()) {
405 tpcLength = track->GetLengthInActiveZone(1, 1.8, 220, track->GetESDEvent()->GetMagneticField());
415 nCRoverFC /= nclsTPCF;
427 track->GetImpactParameters(dca, cov);
429 // (1.) fill number of clusters histogram
431 Int_t minNclsTPC = fESDtrackCuts->GetMinNClusterTPC();
432 fESDtrackCuts->SetMinNClustersTPC(0);
433 if (fESDtrackCuts->AcceptTrack(track)) {
434 for(Int_t iPid = 0; iPid < 6; iPid++) {
435 Double_t vecHistNcl[kNumberOfAxes] = {nclsTPC, pT, eta, phi, iPid};
436 if (IsConsistentWithPid(iPid, track)) histNcl->Fill(vecHistNcl);
439 fESDtrackCuts->SetMinNClustersTPC(minNclsTPC);
441 // (2.) fill chi2 TPC histogram
443 Float_t maxChi2 = fESDtrackCuts->GetMaxChi2PerClusterTPC();
444 fESDtrackCuts->SetMaxChi2PerClusterTPC(999.);
445 if (fESDtrackCuts->AcceptTrack(track)) {
446 for(Int_t iPid = 0; iPid < 6; iPid++) {
447 Double_t vecHistChi2Tpc[kNumberOfAxes] = {chi2TPC, pT, eta, phi, iPid};
448 if (IsConsistentWithPid(iPid, track)) histChi2Tpc->Fill(vecHistChi2Tpc);
451 fESDtrackCuts->SetMaxChi2PerClusterTPC(maxChi2);
453 // (3.) fill dca_z histogram
455 Float_t maxDcaZ = fESDtrackCuts->GetMaxDCAToVertexZ();
456 fESDtrackCuts->SetMaxDCAToVertexZ(999.);
457 if (fESDtrackCuts->AcceptTrack(track)) {
458 for(Int_t iPid = 0; iPid < 6; iPid++) {
459 Double_t vecHistDcaZ[kNumberOfAxes] = {TMath::Abs(dca[1]), pT, eta, phi, iPid};
460 if (IsConsistentWithPid(iPid, track)) histDcaZ->Fill(vecHistDcaZ);
463 fESDtrackCuts->SetMaxDCAToVertexZ(maxDcaZ);
465 // (4.) fill hit in SPD histogram
467 fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
468 if (fESDtrackCuts->AcceptTrack(track)) {
470 if (track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)) hasPoint = 1;
471 for(Int_t iPid = 0; iPid < 6; iPid++) {
472 Double_t vecHistSpd[kNumberOfAxes] = {hasPoint, pT, eta, phi, iPid};
473 if (IsConsistentWithPid(iPid, track)) histSpd->Fill(vecHistSpd);
476 fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
478 // (5.) fill number of crossed rows histogram
480 //Int_t minNcrTPC = fESDtrackCuts->GetMinNCrossedRowsTPC(); //wait for getter in ESDtrackCuts
481 Int_t minNcrTPC = 0; //for now use standard value from 2010 !!
482 fESDtrackCuts->SetMinNCrossedRowsTPC(0);
483 if (fESDtrackCuts->AcceptTrack(track)) {
484 for(Int_t iPid = 0; iPid < 6; iPid++) {
485 Double_t vecHistNcr[kNumberOfAxes] = {ncrTPC, pT, eta, phi, iPid};
486 if (IsConsistentWithPid(iPid, track)) histNcr->Fill(vecHistNcr);
489 fESDtrackCuts->SetMinNCrossedRowsTPC(minNcrTPC);
491 // (6.) fill crossed rows over findable clusters histogram
493 //Int_t minCRoverFC = fESDtrackCuts->GetMinRatioCrossedRowsOverFindableClustersTPC(); //wait for getter in ESDtrackCuts
494 Int_t minCRoverFC = 0.; //for now use standard value from 2010 !!
495 fESDtrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.);
496 if (fESDtrackCuts->AcceptTrack(track)) {
497 for(Int_t iPid = 0; iPid < 6; iPid++) {
498 Double_t vecHistCRoverFC[kNumberOfAxes] = {nCRoverFC, pT, eta, phi, iPid};
499 if (IsConsistentWithPid(iPid, track)) histCRoverFC->Fill(vecHistCRoverFC);
502 fESDtrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(minCRoverFC);
504 // (7.) fill chi2 ITS histogram
506 Float_t maxChi2ITS = fESDtrackCuts->GetMaxChi2PerClusterITS();
507 fESDtrackCuts->SetMaxChi2PerClusterITS(999.);
508 if (fESDtrackCuts->AcceptTrack(track)) {
509 for(Int_t iPid = 0; iPid < 6; iPid++) {
510 Double_t vecHistChi2ITS[kNumberOfAxes] = {chi2ITS, pT, eta, phi, iPid};
511 if (IsConsistentWithPid(iPid, track)) histChi2Its->Fill(vecHistChi2ITS);
514 fESDtrackCuts->SetMaxChi2PerClusterITS(maxChi2ITS);
516 // (8.) fill active length in TPC histogram
518 Int_t minTpcLength = fESDtrackCuts->GetMinLengthActiveVolumeTPC();
519 fESDtrackCuts->SetMinLengthActiveVolumeTPC(0);
520 if (fESDtrackCuts->AcceptTrack(track)) {
521 for(Int_t iPid = 0; iPid < 6; iPid++) {
522 Double_t vecHistTpcLength[kNumberOfAxes] = {tpcLength, pT, eta, phi, iPid};
523 if (IsConsistentWithPid(iPid, track)) histTpcLength->Fill(vecHistTpcLength);
526 fESDtrackCuts->SetMinLengthActiveVolumeTPC(minTpcLength);
528 // (9.) fill TPC->ITS matching efficiency histogram
530 Bool_t isMatched = kFALSE;
531 // remove all ITS requirements
533 // Leonardo and Emilia:
534 // -> beautify this (save the previous cut)
535 // -> check if really all cuts related to the matching are removed
536 // -> histogram with default cuts
537 // -> if MC is available: fill it only for true primaries
538 // -> Postprocessing: plot histogram with 1 divided by histogram with 0 as a function of pT/eta/phi
541 fESDtrackCuts->SetRequireITSRefit(kFALSE);
542 fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
543 fESDtrackCuts->SetMaxChi2TPCConstrainedGlobal(99999.);
544 fESDtrackCuts->SetMaxChi2PerClusterITS(999999.);
546 TString str = fESDtrackCuts->GetMaxDCAToVertexXYPtDep();
547 fESDtrackCuts->SetMaxDCAToVertexXYPtDep();
548 if (fESDtrackCuts->AcceptTrack(track)) {
549 UInt_t status = track->GetStatus();
550 isMatched = status&AliESDtrack::kITSrefit; // this will be obsolete
551 for(Int_t iPid = 0; iPid < 6; iPid++) {
552 Double_t vecHistTpcItsMatch[kNumberOfAxes] = {isMatched, pT, eta, phi, iPid};
553 if (IsConsistentWithPid(iPid, track)) histTpcItsMatch->Fill(vecHistTpcItsMatch); // fill with 0 here
556 fESDtrackCuts->SetRequireITSRefit(kTRUE);
557 fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
558 fESDtrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
559 fESDtrackCuts->SetMaxDCAToVertexXYPtDep(str.Data());
560 fESDtrackCuts->SetMaxChi2PerClusterITS(36);
563 } // end of track loop
570 //________________________________________________________________________
571 void AliAnalysisTrackingUncertainties::Terminate(Option_t *)
573 // Draw result to the screen
574 // Called once at the end of the query
580 //________________________________________________________________________
581 void AliAnalysisTrackingUncertainties::BinLogAxis(const THn *h, Int_t axisNumber) {
583 // Method for the correct logarithmic binning of histograms
585 TAxis *axis = h->GetAxis(axisNumber);
586 int bins = axis->GetNbins();
588 Double_t from = axis->GetXmin();
589 Double_t to = axis->GetXmax();
590 Double_t *newBins = new Double_t[bins + 1];
593 Double_t factor = pow(to/from, 1./bins);
595 for (int i = 1; i <= bins; i++) {
596 newBins[i] = factor * newBins[i-1];
598 axis->Set(bins, newBins);
604 //________________________________________________________________________
605 Bool_t AliAnalysisTrackingUncertainties::IsVertexAccepted(AliESDEvent * esd, Float_t &vertexZ) {
607 // function to check if a proper vertex is reconstructed and write z-position in vertexZ
610 Bool_t vertexOkay = kFALSE;
611 const AliESDVertex *vertex = esd->GetPrimaryVertexTracks();
612 if (vertex->GetNContributors() < 1) {
614 vertex = esd->GetPrimaryVertexSPD();
615 if (vertex->GetNContributors() < 1) {
616 vertexOkay = kFALSE; }
621 TString vtxTyp = vertex->GetTitle();
623 vertex->GetCovarianceMatrix(cov);
624 Double_t zRes = TMath::Sqrt(cov[5]);
625 if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) vertexOkay = kFALSE;
631 vertexZ = vertex->GetZ();
636 //________________________________________________________________________
637 AliAnalysisTrackingUncertainties::ESpecies_t AliAnalysisTrackingUncertainties::GetPid(const AliESDtrack * const tr, Bool_t useTPCTOF) const {
639 // Determine particle species for a given track
640 // Two approaches can be used: As default the selection is done using TPC-only, in addition
641 // the TOF usage is optional. In case of TPC-TOF, a valid TOF signal has to be provided for
642 // the given track. The identification is delegated to helper function for each species.
643 // Tracks which are selected as more than one species (ambiguous decision) are rejected.
645 // @Return: Particles species (kUndef in case no identification is possible)
647 if(!fESDpid) return kUndef;
648 if(useTPCTOF && !(tr->GetStatus() & AliVTrack::kTOFpid)) return kUndef;
650 Bool_t isElectron(kFALSE), isPion(kFALSE), isKaon(kFALSE), isProton(kFALSE);
652 if((isElectron = IsElectron(tr, useTPCTOF))) nspec++;
653 if((isPion = IsPion(tr, useTPCTOF))) nspec++;
654 if((isKaon = IsKaon(tr, useTPCTOF))) nspec++;
655 if((isProton = IsProton(tr,useTPCTOF))) nspec++;
656 if(nspec != 1) return kUndef; // No decision or ambiguous decision;
657 if(isElectron) return kSpecElectron;
658 if(isPion) return kSpecPion;
659 if(isProton) return kSpecProton;
660 if(isKaon) return kSpecKaon;
664 //________________________________________________________________________
665 Bool_t AliAnalysisTrackingUncertainties::IsElectron(const AliESDtrack * const tr, Bool_t useTPCTOF) const {
667 // Selection of electron candidates using the upper half of the TPC sigma band, starting at
668 // the mean ignoring its shift, and going up to 3 sigma above the mean. In case TOF information
669 // is available, tracks which are incompatible with electrons within 3 sigma are rejected. If
670 // no TOF information is used, the momentum regions where the kaon and the proton line cross
671 // the electron line are cut out using a 3 sigma cut around the kaon or proton line.
674 Float_t nsigmaElectronTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kElectron);
675 if(nsigmaElectronTPC < 0 || nsigmaElectronTPC > 3) return kFALSE;
678 Float_t nsigmaElectronTOF = fESDpid->NumberOfSigmasTOF(tr, AliPID::kElectron);
679 if(TMath::Abs(nsigmaElectronTOF) > 3) return kFALSE;
682 Float_t nsigmaKaonTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kKaon),
683 nsigmaProtonTPC =fESDpid->NumberOfSigmasTPC(tr, AliPID::kProton);
684 if(TMath::Abs(nsigmaKaonTPC < 3) || TMath::Abs(nsigmaProtonTPC < 3)) return kFALSE;
689 //________________________________________________________________________
690 Bool_t AliAnalysisTrackingUncertainties::IsConsistentWithPid(Int_t type, const AliESDtrack * const tr) {
692 // just check if the PID is consistent with a given hypothesis in order to
693 // investigate effects which are only dependent on the energy loss.
695 if (type == kSpecPion) return IsPion(tr);
696 if (type == kSpecKaon) return IsKaon(tr);
697 if (type == kSpecProton) return IsProton(tr);
698 if (type == kSpecElectron) return IsElectron(tr);
699 if (type == kAll) return kTRUE;
704 //________________________________________________________________________
705 Bool_t AliAnalysisTrackingUncertainties::IsPion(const AliESDtrack * const tr, Bool_t /*useTPCPTOF*/) const{
707 // Selectron of pion candidates
708 // @TODO: To be implemented
710 Float_t nsigmaPionTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kPion);
711 if (TMath::Abs(nsigmaPionTPC) < 3) return kTRUE;
716 //________________________________________________________________________
717 Bool_t AliAnalysisTrackingUncertainties::IsKaon(const AliESDtrack * const tr, Bool_t /*useTPCPTOF*/) const {
719 // Selection of kaon candidates
720 // @TODO: To be implemented
722 Float_t nsigmaKaonTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kKaon);
723 if (TMath::Abs(nsigmaKaonTPC) < 3) return kTRUE;
728 //________________________________________________________________________
729 Bool_t AliAnalysisTrackingUncertainties::IsProton(const AliESDtrack * const tr, Bool_t /*useTPCPTOF*/) const{
731 // Selection of proton candidates
732 // @TODO: To be implemented
734 Float_t nsigmaProtonTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kProton);
735 if (TMath::Abs(nsigmaProtonTPC) < 3) return kTRUE;
740 //________________________________________________________________________
741 void AliAnalysisTrackingUncertainties::InitializeTrackCutHistograms() {
743 // create histograms for the track cut studies
746 // (1.) number of clusters
747 // 0-ncl, 1-pt, 2-eta, 3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
748 Int_t binsNcl[kNumberOfAxes] = { 40, 50, 20, 18, 6};
749 Double_t minNcl[kNumberOfAxes] = { 20, 0.1, -1, 0, -0.5};
750 Double_t maxNcl[kNumberOfAxes] = {160, 20, +1, 2*TMath::Pi(), 5.5};
752 TString axisNameNcl[kNumberOfAxes] = {"ncl","pT","eta","phi","pid"};
753 TString axisTitleNcl[kNumberOfAxes] = {"ncl","pT","eta","phi","pid"};
755 THnF * histNcl = new THnF("histNcl","number of clusters histogram",kNumberOfAxes, binsNcl, minNcl, maxNcl);
756 BinLogAxis(histNcl, 1);
757 fListHist->Add(histNcl);
759 for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
760 histNcl->GetAxis(iaxis)->SetName(axisNameNcl[iaxis]);
761 histNcl->GetAxis(iaxis)->SetTitle(axisTitleNcl[iaxis]);
765 // 0-chi2, 1-pt, 2-eta, 3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
766 Int_t binsChi2Tpc[kNumberOfAxes] = { 40, 50, 20, 18, 6};
767 Double_t minChi2Tpc[kNumberOfAxes] = { 0, 0.1, -1, 0, -0.5};
768 Double_t maxChi2Tpc[kNumberOfAxes] = { 8, 20, +1, 2*TMath::Pi(), 5.5};
770 TString axisNameChi2Tpc[kNumberOfAxes] = {"chi2tpc","pT","eta","phi","pid"};
771 TString axisTitleChi2Tpc[kNumberOfAxes] = {"chi2tpc","pT","eta","phi","pid"};
773 THnF * histChi2Tpc = new THnF("histChi2Tpc","chi2 per cls. in TPC",kNumberOfAxes, binsChi2Tpc, minChi2Tpc, maxChi2Tpc);
774 BinLogAxis(histChi2Tpc, 1);
775 fListHist->Add(histChi2Tpc);
777 for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
778 histChi2Tpc->GetAxis(iaxis)->SetName(axisNameChi2Tpc[iaxis]);
779 histChi2Tpc->GetAxis(iaxis)->SetTitle(axisTitleChi2Tpc[iaxis]);
783 // 0-dcaZ, 1-pt, 2-eta, 3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
784 Int_t binsDcaZ[kNumberOfAxes] = { 20, 50, 20, 18, 6};
785 Double_t minDcaZ[kNumberOfAxes] = { 0, 0.1, -1, 0, -0.5};
786 Double_t maxDcaZ[kNumberOfAxes] = { 4, 20, +1, 2*TMath::Pi(), 5.5};
788 TString axisNameDcaZ[kNumberOfAxes] = {"dcaZ","pT","eta","phi","pid"};
789 TString axisTitleDcaZ[kNumberOfAxes] = {"dcaZ","pT","eta","phi","pid"};
791 THnF * histDcaZ = new THnF("histDcaZ","dca_z to prim. vtx.",kNumberOfAxes, binsDcaZ, minDcaZ, maxDcaZ);
792 BinLogAxis(histDcaZ, 1);
793 fListHist->Add(histDcaZ);
795 for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
796 histDcaZ->GetAxis(iaxis)->SetName(axisNameDcaZ[iaxis]);
797 histDcaZ->GetAxis(iaxis)->SetTitle(axisTitleDcaZ[iaxis]);
800 // (4.) hit in SPD layer
801 // 0-spdHit, 1-pt, 2-eta, 3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
802 Int_t binsSpd[kNumberOfAxes] = { 2, 50, 20, 18, 6};
803 Double_t minSpd[kNumberOfAxes] = { -0.5, 0.1, -1, 0, -0.5};
804 Double_t maxSpd[kNumberOfAxes] = { 1.5, 20, +1, 2*TMath::Pi(), 5.5};
806 TString axisNameSpd[kNumberOfAxes] = {"spdHit","pT","eta","phi","pid"};
807 TString axisTitleSpd[kNumberOfAxes] = {"spdHit","pT","eta","phi","pid"};
809 THnF * histSpd = new THnF("histSpd","hit in SPD layer or not",kNumberOfAxes, binsSpd, minSpd, maxSpd);
810 BinLogAxis(histSpd, 1);
811 fListHist->Add(histSpd);
813 for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
814 histSpd->GetAxis(iaxis)->SetName(axisNameSpd[iaxis]);
815 histSpd->GetAxis(iaxis)->SetTitle(axisTitleSpd[iaxis]);
818 // (5.) number of crossed rows
819 // 0-ncr, 1-pt, 2-eta, 3-phi, 4-pid(0,unid,etc.)
820 Int_t binsNcr[kNumberOfAxes] = { 40, 50, 20, 18, 6};
821 Double_t minNcr[kNumberOfAxes] = { 20, 0.1, -1, 0, -0.5};
822 Double_t maxNcr[kNumberOfAxes] = {160, 20, +1, 2*TMath::Pi(), 5.5};
824 TString axisNameNcr[kNumberOfAxes] = {"Ncr","pT","eta","phi","pid"};
825 TString axisTitleNcr[kNumberOfAxes] = {"Ncr","pT","eta","phi","pid"};
827 THnF * histNcr = new THnF("histNcr","number of crossed rows TPC histogram",kNumberOfAxes, binsNcr, minNcr, maxNcr);
828 BinLogAxis(histNcr, 1);
829 fListHist->Add(histNcr);
831 for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
832 histNcr->GetAxis(iaxis)->SetName(axisNameNcr[iaxis]);
833 histNcr->GetAxis(iaxis)->SetTitle(axisTitleNcr[iaxis]);
836 // (6.) ratio crossed rows over findable clusters
837 // 0-CRoverFC, 1-pt, 2-eta, 3-phi, 4-pid(0,unid,etc.)
838 Int_t binsCRoverFC[kNumberOfAxes] = { 26, 50, 20, 18, 6};
839 Double_t minCRoverFC[kNumberOfAxes] = { 0.4, 0.1, -1, 0, -0.5};
840 Double_t maxCRoverFC[kNumberOfAxes] = { 1.8, 20, +1, 2*TMath::Pi(), 5.5};
842 TString axisNameCRoverFC[kNumberOfAxes] = {"CRoverFC","pT","eta","phi","pid"};
843 TString axisTitleCRoverFC[kNumberOfAxes] = {"CRoverFC","pT","eta","phi","pid"};
845 THnF * histCRoverFC = new THnF("histCRoverFC","number of crossed rows over findable clusters histogram",kNumberOfAxes, binsCRoverFC, minCRoverFC, maxCRoverFC);
846 BinLogAxis(histCRoverFC, 1);
847 fListHist->Add(histCRoverFC);
849 for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
850 histCRoverFC->GetAxis(iaxis)->SetName(axisNameCRoverFC[iaxis]);
851 histCRoverFC->GetAxis(iaxis)->SetTitle(axisTitleCRoverFC[iaxis]);
854 // (7.) max chi2 / ITS cluster
855 // 0-Chi2Its, 1-pt, 2-eta, 3-phi, 4-pid(0,unid,etc.)
856 Int_t binsChi2Its[kNumberOfAxes] = { 25, 50, 20, 18, 6};
857 Double_t minChi2Its[kNumberOfAxes] = { 0, 0.1, -1, 0, -0.5};
858 Double_t maxChi2Its[kNumberOfAxes] = { 50, 20, +1, 2*TMath::Pi(), 5.5};
860 TString axisNameChi2Its[kNumberOfAxes] = {"Chi2Its","pT","eta","phi","pid"};
861 TString axisTitleChi2Its[kNumberOfAxes] = {"Chi2Its","pT","eta","phi","pid"};
863 THnF * histChi2Its = new THnF("histChi2Its","number of crossed rows TPC histogram",kNumberOfAxes, binsChi2Its, minChi2Its, maxChi2Its);
864 BinLogAxis(histChi2Its, 1);
865 fListHist->Add(histChi2Its);
867 for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
868 histChi2Its->GetAxis(iaxis)->SetName(axisNameChi2Its[iaxis]);
869 histChi2Its->GetAxis(iaxis)->SetTitle(axisTitleChi2Its[iaxis]);
872 // (8.) tpc active volume length
873 // 0-TpcLength, 1-pt, 2-eta, 3-phi, 4-pid(0,unid,etc.)
874 Int_t binsTpcLength[kNumberOfAxes] = { 40, 50, 20, 18, 6};
875 Double_t minTpcLength[kNumberOfAxes] = { 0, 0.1, -1, 0, -0.5};
876 Double_t maxTpcLength[kNumberOfAxes] = { 170, 20, +1, 2*TMath::Pi(), 5.5};
878 TString axisNameTpcLength[kNumberOfAxes] = {"TpcLength","pT","eta","phi","pid"};
879 TString axisTitleTpcLength[kNumberOfAxes] = {"TpcLength","pT","eta","phi","pid"};
881 THnF * histTpcLength = new THnF("histTpcLength","number of crossed rows TPC histogram",kNumberOfAxes, binsTpcLength, minTpcLength, maxTpcLength);
882 BinLogAxis(histTpcLength, 1);
883 fListHist->Add(histTpcLength);
885 for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
886 histTpcLength->GetAxis(iaxis)->SetName(axisNameTpcLength[iaxis]);
887 histTpcLength->GetAxis(iaxis)->SetTitle(axisTitleTpcLength[iaxis]);
890 // (9.) match TPC->ITS
891 // 0-is matched, 1-pt, 2-eta, 3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
892 Int_t binsTpcItsMatch[kNumberOfAxes] = { 2, 50, 20, 18, 6};
893 Double_t minTpcItsMatch[kNumberOfAxes] = { -0.5, 0.1, -1, 0, -0.5};
894 Double_t maxTpcItsMatch[kNumberOfAxes] = { 1.5, 20, +1, 2*TMath::Pi(), 5.5};
896 TString axisNameTpcItsMatch[kNumberOfAxes] = {"isMatched","pT","eta","phi","pid"};
897 TString axisTitleTpcItsMatch[kNumberOfAxes] = {"isMatched","pT","eta","phi","pid"};
899 THnF * histTpcItsMatch = new THnF("histTpcItsMatch","TPC -> ITS matching",kNumberOfAxes, binsTpcItsMatch, minTpcItsMatch, maxTpcItsMatch);
900 BinLogAxis(histTpcItsMatch, 1);
901 fListHist->Add(histTpcItsMatch);
903 for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
904 histTpcItsMatch->GetAxis(iaxis)->SetName(axisNameTpcItsMatch[iaxis]);
905 histTpcItsMatch->GetAxis(iaxis)->SetTitle(axisTitleTpcItsMatch[iaxis]);