]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/EvTrkSelection/AliAnalysisTrackingUncertainties.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
enabling cross checks of em et
[u/mrichter/AliRoot.git] / PWGPP / EvTrkSelection / AliAnalysisTrackingUncertainties.cxx
CommitLineData
b938b7fc 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
16//////////////////////////////////////////////////////////////////////////////
17// //
18// Analysis task for the systematic study of the uncertainties related to //
19// the tracking and ITS-TPC matching efficiency for different particle //
20// species. //
21// //
22//////////////////////////////////////////////////////////////////////////////
23
24
25#include "Riostream.h"
26#include "TH1F.h"
27#include "TH2D.h"
28#include "TH3F.h"
29#include "THn.h"
30#include "TList.h"
31#include "TMath.h"
32#include "TParticlePDG.h"
33//
34#include "AliAnalysisTaskSE.h"
35#include "AliAnalysisManager.h"
36#include "AliPID.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"
46#include "AliStack.h"
47#include "AliLog.h"
48//
49#include "AliAnalysisTrackingUncertainties.h"
50
51
52ClassImp(AliAnalysisTrackingUncertainties)
53
54// histogram constants
55const Int_t kNumberOfAxes = 5;
56
57//________________________________________________________________________
58AliAnalysisTrackingUncertainties::AliAnalysisTrackingUncertainties()
59 : AliAnalysisTaskSE("TaskTestPA"),
60 fESD(0),
61 fESDpid(0),
62 fUtils(0),
63 fMCtrue(0),
64 fListHist(0),
65 fESDtrackCuts(0),
66 fMatchTr(),
67 fMatchChi()
68{
69 // default Constructor
70 /* fast compilation test
71
72 gSystem->Load("libANALYSIS");
73 gSystem->Load("libANALYSISalice");
74 .L AliAnalysisTrackingUncertainties.cxx++
75 */
76}
77
78
79//________________________________________________________________________
80AliAnalysisTrackingUncertainties::AliAnalysisTrackingUncertainties(const char *name)
81 : AliAnalysisTaskSE(name),
82 fESD(0),
83 fESDpid(0),
84 fUtils(0),
85 fMCtrue(0),
86 fListHist(0),
87 fESDtrackCuts(0),
88 fMatchTr(),
89 fMatchChi()
90{
91 //
92 // standard constructur which should be used
93 //
94 fMCtrue = kTRUE;
95 //
96 // create track cuts
97 //
98 fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
99 fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
100 fESDtrackCuts->SetEtaRange(-1., 1.);
101 //
102 // analysis utils if needed
103 //
104 fUtils = new AliAnalysisUtils();
105 //
106
107 // Output slot #0 writes into a TList container
108 DefineOutput(1, TList::Class());
109
110}
111
112
113//________________________________________________________________________
114void AliAnalysisTrackingUncertainties::UserCreateOutputObjects()
115{
116 //
117 // Create histograms
118 // Called once
119 //
120 fListHist = new TList();
121 fListHist->SetOwner(kTRUE);
122 //
a95c2a62 123 // (1.) basic QA and statistics histograms
b938b7fc 124 //
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);
127 //
a95c2a62 128 // (2.) track cut variation histograms
129 //
b938b7fc 130 InitializeTrackCutHistograms();
131 //
a95c2a62 132 // (3.) ITS -> TPC matching histograms
133 //
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};
137 //
138 TString axisNameMatch[kNumberOfAxes] = {"matchChi2","pT","eta","phi","pid"};
139 TString axisTitleMatch[kNumberOfAxes] = {"matchChi2","pT","eta","phi","pid"};
140 //
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]);
145 }
146 BinLogAxis(hBestMatch, 1);
147 fListHist->Add(hBestMatch);
148 //
b938b7fc 149 //
b938b7fc 150 //
151 const int nbPt=40;
152 const double ptMax=5;
153 //
154 TH2F * hNMatch = new TH2F("hNMatch","N Matches",nbPt,0,ptMax,kMaxMatch+1,-0.5,kMaxMatch+0.5);
b938b7fc 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);
b938b7fc 158 fListHist->Add(hAllMatch);
159 fListHist->Add(hAllMatchGlo);
160 //
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);
8ea0d018 169 //add default track cuts in the output list
170 fListHist->Add(fESDtrackCuts);
b938b7fc 171 //
172 // post data
173 //
174 PostData(1, fListHist);
175
176
177}
178
179
180
181//________________________________________________________________________
182void AliAnalysisTrackingUncertainties::UserExec(Option_t *)
183{
184 //
185 // main event loop
186 //
187 fESD = dynamic_cast<AliESDEvent*>( InputEvent() );
188 if (!fESD) {
189 PostData(1, fListHist);
190 return;
191 }
192 //
193 if (!fESDpid) fESDpid = ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid();
194 //
195 // Check Monte Carlo information and other access first:
196 //
197 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
198 if (!eventHandler) {
199 fMCtrue = kFALSE;
200 }
201 //
202 // extract generated particles information
203 //
204 AliMCEvent* mcEvent = 0x0;
205 AliStack* stack = 0x0;
206 if (eventHandler) mcEvent = eventHandler->MCEvent();
207 if (!mcEvent) {
208 if (fMCtrue) return;
209 }
210 if (fMCtrue) {
211 //
212 stack = mcEvent->Stack();
213 if (!stack) return;
214 //
215 for(Int_t i = 0; i < stack->GetNtrack(); i++) {
216 /* at the moment nothing is needed here
217 TParticle * trackMC = stack->Particle(i);
218 Double_t rap = trackMC->Eta();
219 Double_t y = trackMC->Y();
220 Double_t pT = trackMC->Pt();
221 */
222 //
223 }
224
225 }
226 //
227 if (!fESDtrackCuts) {
228 PostData(1, fListHist);
229 return;
230 }
231 //
232 // monitor vertex position and selection
233 //
234 TH2F * histVertexSelection = (TH2F *) fListHist->FindObject("histVertexSelection");
235 //
236 Float_t vertexZ = 0.;
237 if (IsVertexAccepted(fESD, vertexZ)) {
238 histVertexSelection->Fill(vertexZ, 0);
239 } else {
240 histVertexSelection->Fill(vertexZ, 1);
241 return;
242 }
243 //
244 // fill track cut variation histograms
245 //
246 ProcessTrackCutVariation();
247 //
248 // fill ITS->TPC matching histograms
249 //
250 ProcessItsTpcMatching();
251
252 // Post output data
253 PostData(1, fListHist);
254
255}
256
257
258//________________________________________________________________________
259void AliAnalysisTrackingUncertainties::ProcessItsTpcMatching(){
260 //
261 // check how many its-sa tracks get matched to TPC
262 //
263 int ntr = fESD->GetNumberOfTracks();
264 //
265 // initialize histograms
266 //
267 TH2F * hNMatch = (TH2F*) fListHist->FindObject("hNMatch");
a95c2a62 268 THnF * hBestMatch = (THnF*) fListHist->FindObject("hBestMatch");
b938b7fc 269 TH2F * hAllMatch = (TH2F*) fListHist->FindObject("hAllMatch");
270 TH2F * hAllMatchGlo = (TH2F*) fListHist->FindObject("hAllMatchGlo");
271 //
272 TH2F * hNMatchBg = (TH2F*) fListHist->FindObject("hNMatchBg");
273 TH2F * hBestMatchBg = (TH2F*) fListHist->FindObject("hBestMatchBg");
274 TH2F * hAllMatchBg = (TH2F*) fListHist->FindObject("hAllMatchBg");
275 TH2F * hAllMatchGloBg = (TH2F*) fListHist->FindObject("hAllMatchGloBg");
276 //
277 for (int it=0;it<ntr;it++) {
278 AliESDtrack* trSA = fESD->GetTrack(it);
279 if (!trSA->IsOn(AliESDtrack::kITSpureSA) || !trSA->IsOn(AliESDtrack::kITSrefit)) continue;
280 double pt = trSA->Pt();
281 //
282 Int_t nmatch = 0;
48b2ef2f 283 for (int i=kMaxMatch;i--;) {fMatchChi[i]=0; fMatchTr[i]=0;} // reset array
284 for (int it1=0;it1<ntr;it1++) {
b938b7fc 285 if (it1==it) continue;
286 AliESDtrack* trESD = fESD->GetTrack(it1);
287 if (!trESD->IsOn(AliESDtrack::kTPCrefit)) continue;
288 Match(trSA,trESD, nmatch);
289 }
290 //
291 hNMatch->Fill(pt,nmatch);
a95c2a62 292 if (nmatch>0) { // matched tracks
293 Double_t vecHistMatch[kNumberOfAxes] = {fMatchChi[0], pt, trSA->Eta(), trSA->Phi(), 0};
294 hBestMatch->Fill(vecHistMatch);
295 } else { // un-matched tracks -> should be in overflow bin
296 Double_t vecHistMatch[kNumberOfAxes] = {9999999., pt, trSA->Eta(), trSA->Phi(), 0};
297 hBestMatch->Fill(vecHistMatch);
298 }
299 //
b938b7fc 300 for (int imt=nmatch;imt--;) {
301 hAllMatch->Fill(pt,fMatchChi[imt]);
302 if (fMatchTr[imt]->IsOn(AliESDtrack::kITSrefit)) hAllMatchGlo->Fill(pt,fMatchChi[imt]);
303 }
304 //
305 nmatch = 0;
306 for (int i=kMaxMatch;i--;) {fMatchChi[i]=0; fMatchTr[i]=0;}
307 for (int it1=0;it1<ntr;it1++) {
308 if (it1==it) continue;
309 AliESDtrack* trESD = fESD->GetTrack(it1);
310 if (!trESD->IsOn(AliESDtrack::kTPCrefit)) continue;
48b2ef2f 311 Match(trSA,trESD, nmatch, TMath::Pi());
b938b7fc 312 }
313 //
314 hNMatchBg->Fill(pt,nmatch);
315 if (nmatch>0) hBestMatchBg->Fill(pt,fMatchChi[0]);
316 for (int imt=nmatch;imt--;) {
317 hAllMatchBg->Fill(pt,fMatchChi[imt]);
318 if (fMatchTr[imt]->IsOn(AliESDtrack::kITSrefit)) hAllMatchGloBg->Fill(pt,fMatchChi[imt]);
319 }
320 //
321 }
322
323
324}
325
326
48b2ef2f 327void AliAnalysisTrackingUncertainties::Match(const AliESDtrack* tr0, const AliESDtrack* tr1, Int_t &nmatch, Double_t rotate) {
b938b7fc 328 //
329 // check if two tracks are matching, possible rotation for combinatoric backgr.
330 //
331 Float_t bField = fESD->GetMagneticField();
332 //
333 const AliExternalTrackParam* trtpc0 = tr1->GetInnerParam();
334 if (!trtpc0) return;
335 AliExternalTrackParam trtpc(*trtpc0);
336 //
337 if (TMath::Abs(rotate)>1e-5) {
338 const double *par = trtpc.GetParameter();
339 const double *cov = trtpc.GetCovariance();
340 double alp = trtpc.GetAlpha() + rotate;
341 trtpc.Set(trtpc.GetX(),alp,par,cov);
342 }
343 //
344 if (!trtpc.Rotate(tr0->GetAlpha())) return;
345 if (!trtpc.PropagateTo(tr0->GetX(),bField)) return;
346 double chi2 = tr0->GetPredictedChi2(&trtpc);
347 if (chi2>kMaxChi2) return;
348 //
349 int ins;
350 for (ins=0;ins<nmatch;ins++) if (chi2<fMatchChi[ins]) break;
351 if (ins>=kMaxMatch) return;
352
353 for (int imv=nmatch;imv>ins;imv--) {
354 if (imv>=kMaxMatch) continue;
355 fMatchTr[imv] = fMatchTr[imv-1];
356 fMatchChi[imv] = fMatchChi[imv-1];
357 }
358 fMatchTr[ins] = tr1;
359 fMatchChi[ins] = chi2;
360 nmatch++;
361 if (nmatch>=kMaxMatch) nmatch = kMaxMatch;
362 //
363}
364
365
366//________________________________________________________________________
367void AliAnalysisTrackingUncertainties::ProcessTrackCutVariation() {
368 //
369 // fill track cut variation histograms - undo cuts step-by-step and fill histograms
370 //
371 //
372 // initialize histograms
373 //
a95c2a62 374 THnF * histNcl = (THnF *) fListHist->FindObject("histNcl");
375 THnF * histChi2Tpc = (THnF *) fListHist->FindObject("histChi2Tpc");
376 THnF * histDcaZ = (THnF *) fListHist->FindObject("histDcaZ");
377 THnF * histSpd = (THnF *) fListHist->FindObject("histSpd");
378 THnF * histNcr = (THnF *) fListHist->FindObject("histNcr");
379 THnF * histCRoverFC = (THnF *) fListHist->FindObject("histCRoverFC");
380 THnF * histChi2Its = (THnF *) fListHist->FindObject("histChi2Its");
381 THnF * histTpcLength = (THnF *) fListHist->FindObject("histTpcLength");
382 THnF * histTpcItsMatch = (THnF *) fListHist->FindObject("histTpcItsMatch");
b938b7fc 383 //
384 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
385 //
386 for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
387 //
388 AliESDtrack *track =fESD->GetTrack(i);
389 //
390 // relevant variables
391 //
8ea0d018 392 //Double_t pid = Double_t(GetPid(track));
a95c2a62 393 //
394 Int_t nclsTPC = track->GetTPCncls();
395 Float_t pT = track->Pt();
396 Float_t eta = track->Eta();
397 Float_t phi = track->Phi();
398 Float_t chi2TPC = track->GetTPCchi2();
399 Float_t ncrTPC = track->GetTPCCrossedRows();
400 Int_t nclsTPCF = track->GetTPCNclsF();
401 Float_t nCRoverFC = track->GetTPCCrossedRows();
402 Double_t chi2ITS = track->GetITSchi2();
403 Int_t nclsITS = track->GetITSclusters(0);
404 Float_t tpcLength = 0.;
405
406 if (track->GetInnerParam() && track->GetESDEvent()) {
407 tpcLength = track->GetLengthInActiveZone(1, 1.8, 220, track->GetESDEvent()->GetMagneticField());
408 }
409
b938b7fc 410 if (nclsTPC != 0) {
411 chi2TPC /= nclsTPC;
412 } else {
413 chi2TPC = 999.;
414 }
a95c2a62 415
416 if (nclsTPCF !=0) {
417 nCRoverFC /= nclsTPCF;
418 } else {
419 nCRoverFC = 999.;
420 }
421
422 if (nclsITS != 0){
423 chi2ITS /= nclsITS;
424 }else {
425 chi2ITS = 999.;
426 }
b938b7fc 427 //
428 track->GetImpactParameters(dca, cov);
429 //
430 // (1.) fill number of clusters histogram
431 //
432 Int_t minNclsTPC = fESDtrackCuts->GetMinNClusterTPC();
433 fESDtrackCuts->SetMinNClustersTPC(0);
434 if (fESDtrackCuts->AcceptTrack(track)) {
a95c2a62 435 for(Int_t iPid = 0; iPid < 6; iPid++) {
436 Double_t vecHistNcl[kNumberOfAxes] = {nclsTPC, pT, eta, phi, iPid};
437 if (IsConsistentWithPid(iPid, track)) histNcl->Fill(vecHistNcl);
438 }
b938b7fc 439 }
440 fESDtrackCuts->SetMinNClustersTPC(minNclsTPC);
441 //
442 // (2.) fill chi2 TPC histogram
443 //
444 Float_t maxChi2 = fESDtrackCuts->GetMaxChi2PerClusterTPC();
445 fESDtrackCuts->SetMaxChi2PerClusterTPC(999.);
446 if (fESDtrackCuts->AcceptTrack(track)) {
a95c2a62 447 for(Int_t iPid = 0; iPid < 6; iPid++) {
448 Double_t vecHistChi2Tpc[kNumberOfAxes] = {chi2TPC, pT, eta, phi, iPid};
449 if (IsConsistentWithPid(iPid, track)) histChi2Tpc->Fill(vecHistChi2Tpc);
450 }
b938b7fc 451 }
452 fESDtrackCuts->SetMaxChi2PerClusterTPC(maxChi2);
453 //
454 // (3.) fill dca_z histogram
455 //
456 Float_t maxDcaZ = fESDtrackCuts->GetMaxDCAToVertexZ();
457 fESDtrackCuts->SetMaxDCAToVertexZ(999.);
458 if (fESDtrackCuts->AcceptTrack(track)) {
a95c2a62 459 for(Int_t iPid = 0; iPid < 6; iPid++) {
460 Double_t vecHistDcaZ[kNumberOfAxes] = {TMath::Abs(dca[1]), pT, eta, phi, iPid};
461 if (IsConsistentWithPid(iPid, track)) histDcaZ->Fill(vecHistDcaZ);
462 }
b938b7fc 463 }
464 fESDtrackCuts->SetMaxDCAToVertexZ(maxDcaZ);
465 //
466 // (4.) fill hit in SPD histogram
467 //
48b2ef2f 468 fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
b938b7fc 469 if (fESDtrackCuts->AcceptTrack(track)) {
470 Int_t hasPoint = 0;
471 if (track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)) hasPoint = 1;
a95c2a62 472 for(Int_t iPid = 0; iPid < 6; iPid++) {
473 Double_t vecHistSpd[kNumberOfAxes] = {hasPoint, pT, eta, phi, iPid};
474 if (IsConsistentWithPid(iPid, track)) histSpd->Fill(vecHistSpd);
475 }
b938b7fc 476 }
477 fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
a95c2a62 478 //
479 // (5.) fill number of crossed rows histogram
480 //
481 //Int_t minNcrTPC = fESDtrackCuts->GetMinNCrossedRowsTPC(); //wait for getter in ESDtrackCuts
482 Int_t minNcrTPC = 0; //for now use standard value from 2010 !!
483 fESDtrackCuts->SetMinNCrossedRowsTPC(0);
484 if (fESDtrackCuts->AcceptTrack(track)) {
485 for(Int_t iPid = 0; iPid < 6; iPid++) {
486 Double_t vecHistNcr[kNumberOfAxes] = {ncrTPC, pT, eta, phi, iPid};
487 if (IsConsistentWithPid(iPid, track)) histNcr->Fill(vecHistNcr);
488 }
489 }
490 fESDtrackCuts->SetMinNCrossedRowsTPC(minNcrTPC);
491 //
492 // (6.) fill crossed rows over findable clusters histogram
493 //
494 //Int_t minCRoverFC = fESDtrackCuts->GetMinRatioCrossedRowsOverFindableClustersTPC(); //wait for getter in ESDtrackCuts
495 Int_t minCRoverFC = 0.; //for now use standard value from 2010 !!
496 fESDtrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.);
497 if (fESDtrackCuts->AcceptTrack(track)) {
498 for(Int_t iPid = 0; iPid < 6; iPid++) {
499 Double_t vecHistCRoverFC[kNumberOfAxes] = {nCRoverFC, pT, eta, phi, iPid};
500 if (IsConsistentWithPid(iPid, track)) histCRoverFC->Fill(vecHistCRoverFC);
501 }
502 }
503 fESDtrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(minCRoverFC);
504 //
505 // (7.) fill chi2 ITS histogram
506 //
507 Float_t maxChi2ITS = fESDtrackCuts->GetMaxChi2PerClusterITS();
508 fESDtrackCuts->SetMaxChi2PerClusterITS(999.);
509 if (fESDtrackCuts->AcceptTrack(track)) {
510 for(Int_t iPid = 0; iPid < 6; iPid++) {
511 Double_t vecHistChi2ITS[kNumberOfAxes] = {chi2ITS, pT, eta, phi, iPid};
512 if (IsConsistentWithPid(iPid, track)) histChi2Its->Fill(vecHistChi2ITS);
513 }
514 }
515 fESDtrackCuts->SetMaxChi2PerClusterITS(maxChi2ITS);
516 //
517 // (8.) fill active length in TPC histogram
518 //
519 Int_t minTpcLength = fESDtrackCuts->GetMinLengthActiveVolumeTPC();
520 fESDtrackCuts->SetMinLengthActiveVolumeTPC(0);
521 if (fESDtrackCuts->AcceptTrack(track)) {
522 for(Int_t iPid = 0; iPid < 6; iPid++) {
523 Double_t vecHistTpcLength[kNumberOfAxes] = {tpcLength, pT, eta, phi, iPid};
524 if (IsConsistentWithPid(iPid, track)) histTpcLength->Fill(vecHistTpcLength);
525 }
526 }
527 fESDtrackCuts->SetMinLengthActiveVolumeTPC(minTpcLength);
528 //
529 // (9.) fill TPC->ITS matching efficiency histogram
530 //
531 Bool_t isMatched = kFALSE;
532 // remove all ITS requirements
533 //
534 // Leonardo and Emilia:
8ea0d018 535 // -> if MC is available: fill it only for true primaries,
536 // --to be done for every cut?
a95c2a62 537 // -> Postprocessing: plot histogram with 1 divided by histogram with 0 as a function of pT/eta/phi
8ea0d018 538 // -> Do we want to remove the DCA cut?
539 Bool_t refit=fESDtrackCuts->GetRequireITSRefit();
540 Float_t chi2tpc= fESDtrackCuts->GetMaxChi2TPCConstrainedGlobal();
541 Float_t chi2its= fESDtrackCuts->GetMaxChi2PerClusterITS();
542 //TString str = fESDtrackCuts->GetMaxDCAToVertexXYPtDep();
543
a95c2a62 544 fESDtrackCuts->SetRequireITSRefit(kFALSE);
a95c2a62 545 fESDtrackCuts->SetMaxChi2TPCConstrainedGlobal(99999.);
546 fESDtrackCuts->SetMaxChi2PerClusterITS(999999.);
8ea0d018 547 //TString str = fESDtrackCuts->GetMaxDCAToVertexXYPtDep();
548 //fESDtrackCuts->SetMaxDCAToVertexXYPtDep();
549 fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
550
551 if (fESDtrackCuts->AcceptTrack(track)) {
552 for(Int_t iPid = 0; iPid < 6; iPid++) {
553 Double_t vecHistTpcItsMatch[kNumberOfAxes] = {isMatched, pT, eta, phi, iPid};
554 if (IsConsistentWithPid(iPid, track)) histTpcItsMatch->Fill(vecHistTpcItsMatch); // fill with 1 here
555 }
556 }
557 //apply back the cuts
558 fESDtrackCuts->SetRequireITSRefit(refit);
559 fESDtrackCuts->SetMaxChi2TPCConstrainedGlobal(chi2tpc);
560 fESDtrackCuts->SetMaxChi2PerClusterITS(chi2its);
561 //fESDtrackCuts->SetMaxDCAToVertexXYPtDep(str.Data());
562 fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
563 //set is matched
564 isMatched=kTRUE;
a95c2a62 565 if (fESDtrackCuts->AcceptTrack(track)) {
a95c2a62 566 for(Int_t iPid = 0; iPid < 6; iPid++) {
567 Double_t vecHistTpcItsMatch[kNumberOfAxes] = {isMatched, pT, eta, phi, iPid};
568 if (IsConsistentWithPid(iPid, track)) histTpcItsMatch->Fill(vecHistTpcItsMatch); // fill with 0 here
569 }
570 }
8ea0d018 571
b938b7fc 572 } // end of track loop
573
574
575}
576
577
578
579//________________________________________________________________________
580void AliAnalysisTrackingUncertainties::Terminate(Option_t *)
581{
582 // Draw result to the screen
583 // Called once at the end of the query
584
585
586}
587
588
589//________________________________________________________________________
590void AliAnalysisTrackingUncertainties::BinLogAxis(const THn *h, Int_t axisNumber) {
591 //
592 // Method for the correct logarithmic binning of histograms
593 //
594 TAxis *axis = h->GetAxis(axisNumber);
595 int bins = axis->GetNbins();
596
597 Double_t from = axis->GetXmin();
598 Double_t to = axis->GetXmax();
599 Double_t *newBins = new Double_t[bins + 1];
600
601 newBins[0] = from;
602 Double_t factor = pow(to/from, 1./bins);
603
604 for (int i = 1; i <= bins; i++) {
605 newBins[i] = factor * newBins[i-1];
606 }
607 axis->Set(bins, newBins);
608 delete [] newBins;
609
610}
611
612
613//________________________________________________________________________
614Bool_t AliAnalysisTrackingUncertainties::IsVertexAccepted(AliESDEvent * esd, Float_t &vertexZ) {
615 //
616 // function to check if a proper vertex is reconstructed and write z-position in vertexZ
617 //
618 vertexZ = -999.;
619 Bool_t vertexOkay = kFALSE;
620 const AliESDVertex *vertex = esd->GetPrimaryVertexTracks();
621 if (vertex->GetNContributors() < 1) {
622 //
623 vertex = esd->GetPrimaryVertexSPD();
624 if (vertex->GetNContributors() < 1) {
625 vertexOkay = kFALSE; }
626 else {
627 vertexOkay = kTRUE;
628 }
629 //
630 TString vtxTyp = vertex->GetTitle();
631 Double_t cov[6]={0};
632 vertex->GetCovarianceMatrix(cov);
633 Double_t zRes = TMath::Sqrt(cov[5]);
634 if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) vertexOkay = kFALSE;
635 }
636 else {
637 vertexOkay = kTRUE;
638 }
639
640 vertexZ = vertex->GetZ();
641 return vertexOkay;
642
643}
644
4d737435 645//________________________________________________________________________
646AliAnalysisTrackingUncertainties::ESpecies_t AliAnalysisTrackingUncertainties::GetPid(const AliESDtrack * const tr, Bool_t useTPCTOF) const {
647 //
648 // Determine particle species for a given track
649 // Two approaches can be used: As default the selection is done using TPC-only, in addition
650 // the TOF usage is optional. In case of TPC-TOF, a valid TOF signal has to be provided for
651 // the given track. The identification is delegated to helper function for each species.
652 // Tracks which are selected as more than one species (ambiguous decision) are rejected.
653 //
654 // @Return: Particles species (kUndef in case no identification is possible)
655 //
656 if(!fESDpid) return kUndef;
657 if(useTPCTOF && !(tr->GetStatus() & AliVTrack::kTOFpid)) return kUndef;
658
659 Bool_t isElectron(kFALSE), isPion(kFALSE), isKaon(kFALSE), isProton(kFALSE);
660 Int_t nspec(0);
661 if((isElectron = IsElectron(tr, useTPCTOF))) nspec++;
662 if((isPion = IsPion(tr, useTPCTOF))) nspec++;
663 if((isKaon = IsKaon(tr, useTPCTOF))) nspec++;
664 if((isProton = IsProton(tr,useTPCTOF))) nspec++;
665 if(nspec != 1) return kUndef; // No decision or ambiguous decision;
666 if(isElectron) return kSpecElectron;
a95c2a62 667 if(isPion) return kSpecPion;
668 if(isProton) return kSpecProton;
669 if(isKaon) return kSpecKaon;
4d737435 670 return kUndef;
671}
672
673//________________________________________________________________________
674Bool_t AliAnalysisTrackingUncertainties::IsElectron(const AliESDtrack * const tr, Bool_t useTPCTOF) const {
675 //
676 // Selection of electron candidates using the upper half of the TPC sigma band, starting at
677 // the mean ignoring its shift, and going up to 3 sigma above the mean. In case TOF information
678 // is available, tracks which are incompatible with electrons within 3 sigma are rejected. If
679 // no TOF information is used, the momentum regions where the kaon and the proton line cross
680 // the electron line are cut out using a 3 sigma cut around the kaon or proton line.
681 //
682
683 Float_t nsigmaElectronTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kElectron);
684 if(nsigmaElectronTPC < 0 || nsigmaElectronTPC > 3) return kFALSE;
685
686 if(useTPCTOF){
687 Float_t nsigmaElectronTOF = fESDpid->NumberOfSigmasTOF(tr, AliPID::kElectron);
688 if(TMath::Abs(nsigmaElectronTOF) > 3) return kFALSE;
689 else return kTRUE;
690 } else {
691 Float_t nsigmaKaonTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kKaon),
692 nsigmaProtonTPC =fESDpid->NumberOfSigmasTPC(tr, AliPID::kProton);
693 if(TMath::Abs(nsigmaKaonTPC < 3) || TMath::Abs(nsigmaProtonTPC < 3)) return kFALSE;
694 else return kTRUE;
695 }
696}
697
698//________________________________________________________________________
a95c2a62 699Bool_t AliAnalysisTrackingUncertainties::IsConsistentWithPid(Int_t type, const AliESDtrack * const tr) {
700 //
701 // just check if the PID is consistent with a given hypothesis in order to
702 // investigate effects which are only dependent on the energy loss.
703 //
704 if (type == kSpecPion) return IsPion(tr);
705 if (type == kSpecKaon) return IsKaon(tr);
706 if (type == kSpecProton) return IsProton(tr);
707 if (type == kSpecElectron) return IsElectron(tr);
708 if (type == kAll) return kTRUE;
709 return kFALSE;
710
4d737435 711}
712
713//________________________________________________________________________
a95c2a62 714Bool_t AliAnalysisTrackingUncertainties::IsPion(const AliESDtrack * const tr, Bool_t /*useTPCPTOF*/) const{
715 //
716 // Selectron of pion candidates
717 // @TODO: To be implemented
718 //
719 Float_t nsigmaPionTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kPion);
720 if (TMath::Abs(nsigmaPionTPC) < 3) return kTRUE;
721 return kFALSE;
722
4d737435 723}
724
725//________________________________________________________________________
a95c2a62 726Bool_t AliAnalysisTrackingUncertainties::IsKaon(const AliESDtrack * const tr, Bool_t /*useTPCPTOF*/) const {
727 //
728 // Selection of kaon candidates
729 // @TODO: To be implemented
730 //
731 Float_t nsigmaKaonTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kKaon);
732 if (TMath::Abs(nsigmaKaonTPC) < 3) return kTRUE;
733 return kFALSE;
734
735}
736
737//________________________________________________________________________
738Bool_t AliAnalysisTrackingUncertainties::IsProton(const AliESDtrack * const tr, Bool_t /*useTPCPTOF*/) const{
739 //
740 // Selection of proton candidates
741 // @TODO: To be implemented
742 //
743 Float_t nsigmaProtonTPC = fESDpid->NumberOfSigmasTPC(tr, AliPID::kProton);
744 if (TMath::Abs(nsigmaProtonTPC) < 3) return kTRUE;
745 return kFALSE;
746
4d737435 747}
b938b7fc 748
749//________________________________________________________________________
750void AliAnalysisTrackingUncertainties::InitializeTrackCutHistograms() {
751 //
752 // create histograms for the track cut studies
753 //
754 //
755 // (1.) number of clusters
a95c2a62 756 // 0-ncl, 1-pt, 2-eta, 3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
b938b7fc 757 Int_t binsNcl[kNumberOfAxes] = { 40, 50, 20, 18, 6};
758 Double_t minNcl[kNumberOfAxes] = { 20, 0.1, -1, 0, -0.5};
759 Double_t maxNcl[kNumberOfAxes] = {160, 20, +1, 2*TMath::Pi(), 5.5};
760 //
761 TString axisNameNcl[kNumberOfAxes] = {"ncl","pT","eta","phi","pid"};
762 TString axisTitleNcl[kNumberOfAxes] = {"ncl","pT","eta","phi","pid"};
763 //
764 THnF * histNcl = new THnF("histNcl","number of clusters histogram",kNumberOfAxes, binsNcl, minNcl, maxNcl);
765 BinLogAxis(histNcl, 1);
766 fListHist->Add(histNcl);
767 //
768 for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
769 histNcl->GetAxis(iaxis)->SetName(axisNameNcl[iaxis]);
770 histNcl->GetAxis(iaxis)->SetTitle(axisTitleNcl[iaxis]);
771 }
772 //
773 // (2.) chi2/cls-TPC
a95c2a62 774 // 0-chi2, 1-pt, 2-eta, 3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
b938b7fc 775 Int_t binsChi2Tpc[kNumberOfAxes] = { 40, 50, 20, 18, 6};
776 Double_t minChi2Tpc[kNumberOfAxes] = { 0, 0.1, -1, 0, -0.5};
777 Double_t maxChi2Tpc[kNumberOfAxes] = { 8, 20, +1, 2*TMath::Pi(), 5.5};
778 //
779 TString axisNameChi2Tpc[kNumberOfAxes] = {"chi2tpc","pT","eta","phi","pid"};
780 TString axisTitleChi2Tpc[kNumberOfAxes] = {"chi2tpc","pT","eta","phi","pid"};
781 //
a95c2a62 782 THnF * histChi2Tpc = new THnF("histChi2Tpc","chi2 per cls. in TPC",kNumberOfAxes, binsChi2Tpc, minChi2Tpc, maxChi2Tpc);
b938b7fc 783 BinLogAxis(histChi2Tpc, 1);
784 fListHist->Add(histChi2Tpc);
785 //
786 for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
787 histChi2Tpc->GetAxis(iaxis)->SetName(axisNameChi2Tpc[iaxis]);
788 histChi2Tpc->GetAxis(iaxis)->SetTitle(axisTitleChi2Tpc[iaxis]);
789 }
790 //
791 // (3.) dca_z
a95c2a62 792 // 0-dcaZ, 1-pt, 2-eta, 3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
48b2ef2f 793 Int_t binsDcaZ[kNumberOfAxes] = { 20, 50, 20, 18, 6};
b938b7fc 794 Double_t minDcaZ[kNumberOfAxes] = { 0, 0.1, -1, 0, -0.5};
48b2ef2f 795 Double_t maxDcaZ[kNumberOfAxes] = { 4, 20, +1, 2*TMath::Pi(), 5.5};
b938b7fc 796 //
797 TString axisNameDcaZ[kNumberOfAxes] = {"dcaZ","pT","eta","phi","pid"};
798 TString axisTitleDcaZ[kNumberOfAxes] = {"dcaZ","pT","eta","phi","pid"};
799 //
a95c2a62 800 THnF * histDcaZ = new THnF("histDcaZ","dca_z to prim. vtx.",kNumberOfAxes, binsDcaZ, minDcaZ, maxDcaZ);
b938b7fc 801 BinLogAxis(histDcaZ, 1);
802 fListHist->Add(histDcaZ);
803 //
804 for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
805 histDcaZ->GetAxis(iaxis)->SetName(axisNameDcaZ[iaxis]);
806 histDcaZ->GetAxis(iaxis)->SetTitle(axisTitleDcaZ[iaxis]);
807 }
808 //
809 // (4.) hit in SPD layer
a95c2a62 810 // 0-spdHit, 1-pt, 2-eta, 3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
b938b7fc 811 Int_t binsSpd[kNumberOfAxes] = { 2, 50, 20, 18, 6};
812 Double_t minSpd[kNumberOfAxes] = { -0.5, 0.1, -1, 0, -0.5};
813 Double_t maxSpd[kNumberOfAxes] = { 1.5, 20, +1, 2*TMath::Pi(), 5.5};
814 //
815 TString axisNameSpd[kNumberOfAxes] = {"spdHit","pT","eta","phi","pid"};
816 TString axisTitleSpd[kNumberOfAxes] = {"spdHit","pT","eta","phi","pid"};
817 //
a95c2a62 818 THnF * histSpd = new THnF("histSpd","hit in SPD layer or not",kNumberOfAxes, binsSpd, minSpd, maxSpd);
b938b7fc 819 BinLogAxis(histSpd, 1);
820 fListHist->Add(histSpd);
821 //
822 for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
823 histSpd->GetAxis(iaxis)->SetName(axisNameSpd[iaxis]);
824 histSpd->GetAxis(iaxis)->SetTitle(axisTitleSpd[iaxis]);
825 }
a95c2a62 826 //
827 // (5.) number of crossed rows
828 // 0-ncr, 1-pt, 2-eta, 3-phi, 4-pid(0,unid,etc.)
829 Int_t binsNcr[kNumberOfAxes] = { 40, 50, 20, 18, 6};
830 Double_t minNcr[kNumberOfAxes] = { 20, 0.1, -1, 0, -0.5};
831 Double_t maxNcr[kNumberOfAxes] = {160, 20, +1, 2*TMath::Pi(), 5.5};
832 //
833 TString axisNameNcr[kNumberOfAxes] = {"Ncr","pT","eta","phi","pid"};
834 TString axisTitleNcr[kNumberOfAxes] = {"Ncr","pT","eta","phi","pid"};
835 //
836 THnF * histNcr = new THnF("histNcr","number of crossed rows TPC histogram",kNumberOfAxes, binsNcr, minNcr, maxNcr);
837 BinLogAxis(histNcr, 1);
838 fListHist->Add(histNcr);
839 //
840 for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
841 histNcr->GetAxis(iaxis)->SetName(axisNameNcr[iaxis]);
842 histNcr->GetAxis(iaxis)->SetTitle(axisTitleNcr[iaxis]);
843 }
844 //
845 // (6.) ratio crossed rows over findable clusters
846 // 0-CRoverFC, 1-pt, 2-eta, 3-phi, 4-pid(0,unid,etc.)
847 Int_t binsCRoverFC[kNumberOfAxes] = { 26, 50, 20, 18, 6};
848 Double_t minCRoverFC[kNumberOfAxes] = { 0.4, 0.1, -1, 0, -0.5};
849 Double_t maxCRoverFC[kNumberOfAxes] = { 1.8, 20, +1, 2*TMath::Pi(), 5.5};
850 //
851 TString axisNameCRoverFC[kNumberOfAxes] = {"CRoverFC","pT","eta","phi","pid"};
852 TString axisTitleCRoverFC[kNumberOfAxes] = {"CRoverFC","pT","eta","phi","pid"};
853 //
854 THnF * histCRoverFC = new THnF("histCRoverFC","number of crossed rows over findable clusters histogram",kNumberOfAxes, binsCRoverFC, minCRoverFC, maxCRoverFC);
855 BinLogAxis(histCRoverFC, 1);
856 fListHist->Add(histCRoverFC);
857 //
858 for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
859 histCRoverFC->GetAxis(iaxis)->SetName(axisNameCRoverFC[iaxis]);
860 histCRoverFC->GetAxis(iaxis)->SetTitle(axisTitleCRoverFC[iaxis]);
861 }
862 //
863 // (7.) max chi2 / ITS cluster
864 // 0-Chi2Its, 1-pt, 2-eta, 3-phi, 4-pid(0,unid,etc.)
865 Int_t binsChi2Its[kNumberOfAxes] = { 25, 50, 20, 18, 6};
866 Double_t minChi2Its[kNumberOfAxes] = { 0, 0.1, -1, 0, -0.5};
867 Double_t maxChi2Its[kNumberOfAxes] = { 50, 20, +1, 2*TMath::Pi(), 5.5};
868 //
869 TString axisNameChi2Its[kNumberOfAxes] = {"Chi2Its","pT","eta","phi","pid"};
870 TString axisTitleChi2Its[kNumberOfAxes] = {"Chi2Its","pT","eta","phi","pid"};
871 //
872 THnF * histChi2Its = new THnF("histChi2Its","number of crossed rows TPC histogram",kNumberOfAxes, binsChi2Its, minChi2Its, maxChi2Its);
873 BinLogAxis(histChi2Its, 1);
874 fListHist->Add(histChi2Its);
875 //
876 for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
877 histChi2Its->GetAxis(iaxis)->SetName(axisNameChi2Its[iaxis]);
878 histChi2Its->GetAxis(iaxis)->SetTitle(axisTitleChi2Its[iaxis]);
879 }
880 //
881 // (8.) tpc active volume length
882 // 0-TpcLength, 1-pt, 2-eta, 3-phi, 4-pid(0,unid,etc.)
883 Int_t binsTpcLength[kNumberOfAxes] = { 40, 50, 20, 18, 6};
884 Double_t minTpcLength[kNumberOfAxes] = { 0, 0.1, -1, 0, -0.5};
885 Double_t maxTpcLength[kNumberOfAxes] = { 170, 20, +1, 2*TMath::Pi(), 5.5};
886 //
887 TString axisNameTpcLength[kNumberOfAxes] = {"TpcLength","pT","eta","phi","pid"};
888 TString axisTitleTpcLength[kNumberOfAxes] = {"TpcLength","pT","eta","phi","pid"};
889 //
890 THnF * histTpcLength = new THnF("histTpcLength","number of crossed rows TPC histogram",kNumberOfAxes, binsTpcLength, minTpcLength, maxTpcLength);
891 BinLogAxis(histTpcLength, 1);
892 fListHist->Add(histTpcLength);
893 //
894 for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
895 histTpcLength->GetAxis(iaxis)->SetName(axisNameTpcLength[iaxis]);
896 histTpcLength->GetAxis(iaxis)->SetTitle(axisTitleTpcLength[iaxis]);
897 }
898 //
899 // (9.) match TPC->ITS
900 // 0-is matched, 1-pt, 2-eta, 3-phi, 4-pid(0-3 -> electron-proton, 4 -> undef, 5 -> all)
901 Int_t binsTpcItsMatch[kNumberOfAxes] = { 2, 50, 20, 18, 6};
902 Double_t minTpcItsMatch[kNumberOfAxes] = { -0.5, 0.1, -1, 0, -0.5};
903 Double_t maxTpcItsMatch[kNumberOfAxes] = { 1.5, 20, +1, 2*TMath::Pi(), 5.5};
904 //
905 TString axisNameTpcItsMatch[kNumberOfAxes] = {"isMatched","pT","eta","phi","pid"};
906 TString axisTitleTpcItsMatch[kNumberOfAxes] = {"isMatched","pT","eta","phi","pid"};
907 //
908 THnF * histTpcItsMatch = new THnF("histTpcItsMatch","TPC -> ITS matching",kNumberOfAxes, binsTpcItsMatch, minTpcItsMatch, maxTpcItsMatch);
909 BinLogAxis(histTpcItsMatch, 1);
910 fListHist->Add(histTpcItsMatch);
911 //
912 for (Int_t iaxis=0; iaxis<kNumberOfAxes;iaxis++){
913 histTpcItsMatch->GetAxis(iaxis)->SetName(axisNameTpcItsMatch[iaxis]);
914 histTpcItsMatch->GetAxis(iaxis)->SetTitle(axisTitleTpcItsMatch[iaxis]);
915 }
b938b7fc 916
917
918
919
920}
921