AliPhysics  master (68a99cd)
AliConversionTrackCuts.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * This file is property of and copyright by the ALICE HLT Project *
3  * ALICE Experiment at CERN, All rights reserved. *
4  * *
5  * Primary Author: Svein Lindal <slindal@fys.uio.no> *
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, provided 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 purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
19 
20 
21 #include "AliConversionTrackCuts.h"
22 //#include "AliAODTrack.h"
23 #include "AliAODEvent.h"
24 #include <TFormula.h>
25 #include <iostream>
26 #include "TH2F.h"
27 #include "AliESDtrackCuts.h"
28 #include "THn.h"
29 
30 using namespace std;
31 ClassImp(AliConversionTrackCuts)
32 
33 
34 const char* AliConversionTrackCuts::fgkCutNames[AliConversionTrackCuts::kNCuts] = {
35  "nClusTPC",
36  "FoundFindable",
37  "Chi2PerNDF",
38  "Kink",
39  "DCA_Z",
40  "DCA_XY",
41  "TPCRefit"
42  "kAccTracks"
43 };
44 
45 
46 
47 //________________________________________________________________________
49  AliAnalysisCuts(),
50  fEsdTrackCuts(NULL),
51  fEsdTrackCutsExtra1(NULL),
52  fEsdTrackCutsExtra2(NULL),
53  fEvent(NULL),
54  fFilterBit(2048),
55  fDCAZmax(3.2*3.2),
56  fDCAXYmax(2.4*2.4),
57  fInitialized(kFALSE),
58  fhPhi(NULL),
59  // fhPt(NULL),
60  //fhPhiPt(NULL),
61  fhdcaxyPt(NULL),
62  fhdcazPt(NULL),
63  fhdca(NULL),
64  fhnclpt(NULL),
65  fhnclsfpt(NULL),
66  fhEtaPhi(NULL),
67  fhTrackEff(NULL),
68  fkCreateTrackEff(kFALSE),
69  fHistograms(NULL)
70 {
71  //Constructor
72 }
73 //________________________________________________________________________
75  AliAnalysisCuts(name, title),
76  fEsdTrackCuts(NULL),
77  fEsdTrackCutsExtra1(NULL),
78  fEsdTrackCutsExtra2(NULL),
79  fEvent(NULL),
80  fFilterBit(2048),
81  fDCAZmax(-1),
82  fDCAXYmax(-1),
83  fInitialized(kFALSE),
84  fhPhi(NULL),
85  //fhPt(NULL),
86  //fhPhiPt(NULL),
87  fhdcaxyPt(NULL),
88  fhdcazPt(NULL),
89  fhdca(NULL),
90  fhnclpt(NULL),
91  fhnclsfpt(NULL),
92  fhEtaPhi(NULL),
93  fhTrackEff(NULL),
94  fkCreateTrackEff(kFALSE),
95  fHistograms(NULL)
96 {
97  //Constructor
98 }
99 
100 
101 //________________________________________________________________________________
104  // if(fHistograms)
105  // delete fHistograms;
106  // fHistograms = NULL;
107 
108  if(fEsdTrackCuts)
109  delete fEsdTrackCuts;
110  fEsdTrackCuts = NULL;
111 
113  delete fEsdTrackCutsExtra1;
114  fEsdTrackCutsExtra1 = NULL;
115 
117  delete fEsdTrackCutsExtra2;
118  fEsdTrackCutsExtra2 = NULL;
119 
120 }
121 
122 //______________________________________________________________________________
124  // Reproduces the cuts of the corresponding bit in the ESD->AOD filtering
125  // (see $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C)
127  const Int_t filterbit = fFilterBit;
128 
129  if (filterbit == 128) {
130  if(!fEsdTrackCuts) {
131  fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
132  fEsdTrackCuts->SetMinNClustersTPC(70);
133  }
134  } else if (filterbit == 256) {
135  if(!fEsdTrackCuts) {
136  // syst study
137  fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
138  fEsdTrackCuts->SetMinNClustersTPC(80);
139  fEsdTrackCuts->SetMaxChi2PerClusterTPC(3);
140  fEsdTrackCuts->SetMaxDCAToVertexZ(2.7);
141  fEsdTrackCuts->SetMaxDCAToVertexXY(1.9);
142  }
143  } else if (filterbit == 512) {
144  if(!fEsdTrackCuts) {
145  // syst study
146  fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
147  fEsdTrackCuts->SetMinNClustersTPC(60);
148  fEsdTrackCuts->SetMaxChi2PerClusterTPC(5);
149  fEsdTrackCuts->SetMaxDCAToVertexZ(3.7);
150  fEsdTrackCuts->SetMaxDCAToVertexXY(2.9);
151  }
152  } else if (filterbit == 1024) {
153  if(!fEsdTrackCuts) {
154  fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
155  fEsdTrackCuts->SetMinNClustersTPC(-1);
156  fEsdTrackCuts->SetMinNCrossedRowsTPC(70);
157  fEsdTrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
158  }
159  } else if (filterbit == 2048) {
160  // mimic hybrid tracks
161  // correspond to esdTrackCutsHTG, but WITHOUT spd constraint. this is checked with the next object
162  if(!fEsdTrackCuts) {
163  fEsdTrackCuts = new AliESDtrackCuts();
164  TFormula *f1NClustersTPCLinearPtDep = new TFormula("f1NClustersTPCLinearPtDep","70.+30./20.*x");
165  fEsdTrackCuts->SetMinNClustersTPCPtDep(f1NClustersTPCLinearPtDep, 100);
166  fEsdTrackCuts->SetMaxChi2PerClusterTPC(4);
167  fEsdTrackCuts->SetRequireTPCStandAlone(kTRUE);
168  fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
169  fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
170  fEsdTrackCuts->SetMaxFractionSharedTPCClusters(0.4);
171 
172  fEsdTrackCuts->SetMaxDCAToVertexXY(2.4);
173  fEsdTrackCuts->SetMaxDCAToVertexZ(3.2);
174  fEsdTrackCuts->SetDCAToVertex2D(kTRUE);
175 
176  fEsdTrackCuts->SetMaxChi2PerClusterITS(36);
177  fEsdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
178 
179  fEsdTrackCuts->SetRequireSigmaToVertex(kFALSE);
180 
181  fEsdTrackCuts->SetEtaRange(-0.9, 0.9);
182  fEsdTrackCuts->SetPtRange(0.1, 1000000.0);
183 
184  fEsdTrackCuts->SetRequireITSRefit(kFALSE); //not here, n
185  }
186  // Add SPD requirement
187  fEsdTrackCutsExtra1 = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");
188  fEsdTrackCutsExtra1->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
189  fEsdTrackCutsExtra1->SetRequireITSRefit(kTRUE);
190  // A track passing fEsdTrackCuts and fEsdTrackCutsExtra1 corresponds to esdTrackCutsHTG
191 
192  fEsdTrackCutsExtra2 = new AliESDtrackCuts("No_SPD", "Reject tracks with cluster in SPD");
193  fEsdTrackCutsExtra2->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kNone);
194  // A track passing fEsdTrackCuts and fEsdTrackCutsExtra2 corresponds to esdTrackCutsHTGC and needs to be constrained
195 
196 
197  }
198 }
199 
200 
201 //______________________________________________________________________________
202 Bool_t AliConversionTrackCuts::AcceptTrack(const AliESDtrack * track) {
203  //Check esd track
204  FillHistograms(kPreCut, track);
205 
206  if( fFilterBit == 256) {
207 
209  const AliExternalTrackParam * param = track->GetConstrainedParam();
210  if(param) {
211  AliESDtrack esdTrack(track);
212  esdTrack.CopyFromVTrack(param);
213 
214  if( !fEsdTrackCuts->IsSelected(&esdTrack)) return kFALSE;
215 
216  FillHistograms(1, track);
217 
218  Double_t dca[2];
219  GetDCA(&esdTrack, dca);
220 
221  FillDCAHist(dca[1], dca[0], &esdTrack);
222  if(fhEtaPhi) fhEtaPhi->Fill(esdTrack.Eta(), esdTrack.Phi());
223  return kTRUE;
224  } else {
225  return kFALSE;
226  }
227 
228  return kFALSE;
229  }
230 
231 
232  if(!fInitialized) {
233  DefineESDCuts();
234  // if(fDCAXYmax > 0) {
235  // if(fEsdTrackCuts) fEsdTrackCuts->SetMaxDCAToVertexXY(fDCAXYmax);
236  // }
237  // if(fDCAZmax > 0) {
238  // if(fEsdTrackCuts) fEsdTrackCuts->SetMaxDCAToVertexZ(fDCAZmax);
239  // }
240 
241  fInitialized = kTRUE;
242  }
243 
244 
245  Double_t dca[2];
246  GetDCA(track, dca);
247 
248 
251  FillHistograms(1, track);
252  FillDCAHist(dca[1], dca[0], track);
253  if(fhEtaPhi) fhEtaPhi->Fill(track->Eta(), track->Phi());
254  return kTRUE;
255  }
256 
258  AliESDtrack internaltrack(track); // local copy used in order to keep constantnes of the input object
259  if (fEsdTrackCutsExtra1 && fEsdTrackCutsExtra1->IsSelected(&internaltrack)) {
260  FillHistograms(1, &internaltrack);
261  FillHistograms(2, &internaltrack);
262 
263  FillDCAHist(dca[1], dca[0], &internaltrack);
264  if(fhEtaPhi) fhEtaPhi->Fill(internaltrack.Eta(), internaltrack.Phi());
265 
266  return kTRUE;
267  }
268 
270  if (fEsdTrackCutsExtra2 && fEsdTrackCutsExtra2->IsSelected(&internaltrack)) {
271  const AliExternalTrackParam * param = internaltrack.GetConstrainedParam();
272  if(param) {
273  AliESDtrack esdTrack(internaltrack);
274  esdTrack.CopyFromVTrack(param);
275 
276  FillHistograms(3, &esdTrack);
277  FillHistograms(1, &esdTrack);
278 
279  FillDCAHist(dca[1], dca[0], &esdTrack);
280  if(fhEtaPhi) fhEtaPhi->Fill(esdTrack.Eta(), esdTrack.Phi());
281 
282  return kTRUE;
283  } else {
284  return kFALSE;
285  }
286  } else {
287  return kFALSE;
288  }
289 
290  cout << "error error, should not be herer!"<<endl;
291  return kFALSE;
292 
293  // FillHistograms(kPreCut + 1, track);
294  // return kTRUE;
295 
296  // fhnclpt->Fill(track->Pt(), track->GetTPCNcls());
297  // if(track->GetTPCNclsF() > 0) fhnclsfpt->Fill(track->Pt(), ((Double_t) track->GetTPCNcls())/track->GetTPCNclsF());
298  // FillHistograms(kPreCut + 1, track);
299 
300  // ///Get impact parameters
301  // Double_t extCov[15];
302  // track->GetExternalCovariance(extCov);
303  // return kTRUE;
304 }
305 
306 Bool_t AliConversionTrackCuts::AcceptTrack(const AliAODTrack * track) {
307  //Check aod track
308 
309  FillHistograms(kPreCut, track);
310 
311  if (fFilterBit == 768) {
312  if(!track->IsHybridGlobalConstrainedGlobal()) return kFALSE;
313 
314  if (!(track->GetStatus() & AliVTrack::kITSrefit)) {
315  return kFALSE;
316  }
317 
318  //The cluster sharing cut can be done with:
319  Double_t frac = Double_t(track->GetTPCnclsS()) / Double_t(track->GetTPCncls());
320  if (frac > 0.4) return kFALSE;
321 
323  FillHistograms(1, track);
324 
326  Double_t dca[2] = { -999, -999};
327  //Bool_t dcaok =
328  GetDCA(track, dca);
329  FillDCAHist(dca[1], dca[0], track);
330 
331 
332  if(track->IsGlobalConstrained()) {
333  FillHistograms(3, track);
334  } else {
335  FillHistograms(2, track);
336  }
337 
338  if(fhEtaPhi) fhEtaPhi->Fill(track->Eta(), track->Phi());
339 
340  return kTRUE;
341 
345  } else if(fFilterBit == 256) {
346  if(!track->IsTPCConstrained()) return kFALSE;
347 
348 
349 
351  Double_t dca[2] = { -999, -999};
352  GetDCA(track, dca);
353 
354  if( (dca[0]*dca[0]/fDCAXYmax + dca[1]*dca[1]/fDCAZmax) > 1 ) {
355  FillHistograms(3, track);
356  return kFALSE;
357  }
358 
359  if(track->GetTPCncls() < 70) {
360  FillHistograms(4, track);
361  return kFALSE;
362  }
363 
364  AliAODVertex * vtx = track->GetProdVertex();
365  if (vtx->GetType() == AliAODVertex::kKink ) {
366  FillHistograms(5, track);
367  return kFALSE;
368  }
369 
370  if(track->Chi2perNDF() > 36) {
371  FillHistograms(6, track);
372  return kFALSE;
373  }
374  if(track->Chi2perNDF() > 26) {
375  FillHistograms(7, track);
376  return kFALSE;
377  }
378  if(track->Chi2perNDF() > 16) {
379  FillHistograms(8, track);
380  return kFALSE;
381  }
382  if(track->Chi2perNDF() > 4) {
383  FillHistograms(9, track);
384  return kFALSE;
385  }
386 
387 
388 
389  FillDCAHist(dca[1], dca[0], track);
390 
391  FillHistograms(2, track);
392  if(fhEtaPhi) fhEtaPhi->Fill(track->Eta(), track->Phi());
393  return kTRUE;
394 
395  }
396  return kFALSE;
397 }
398 
399 
401 Bool_t AliConversionTrackCuts::GetDCA(const AliESDtrack *track, Double_t dcaxyz[2]) {
403  Float_t dca[2];
404  Float_t bCov[3];
405  track->GetImpactParameters(dca,bCov);
406  if (bCov[0]<=0 || bCov[2]<=0) {
407  AliDebug(1, "Estimated b resolution lower or equal zero!");
408  bCov[0]=0; bCov[2]=0;
409  return kFALSE;
410  }
411 
412  dcaxyz[0] = dca[0];
413  dcaxyz[1] = dca[1];
414 
415  return kTRUE;
416 }
417 
419 Bool_t AliConversionTrackCuts::GetDCA(const AliAODTrack *track, Double_t dca[2]) {
421  if(track->TestBit(AliAODTrack::kIsDCA)){
422  dca[0]=track->DCA();
423  dca[1]=track->ZAtDCA();
424  return kTRUE;
425  }
426 
427  Bool_t ok=kFALSE;
428  if(fEvent) {
429  Double_t covdca[3];
430  //AliAODTrack copy(*track);
431  AliExternalTrackParam etp; etp.CopyFromVTrack(track);
432 
433  Float_t xstart = etp.GetX();
434  if(xstart>3.) {
435  dca[0]=-999.;
436  dca[1]=-999.;
437  //printf("This method can be used only for propagation inside the beam pipe \n");
438  return kFALSE;
439  }
440 
441 
442  AliAODVertex *vtx =(AliAODVertex*)(fEvent->GetPrimaryVertex());
443  Double_t fBzkG = fEvent->GetMagneticField(); // z componenent of field in kG
444  ok = etp.PropagateToDCA(vtx,fBzkG,kVeryBig,dca,covdca);
445  //ok = copy.PropagateToDCA(vtx,fBzkG,kVeryBig,dca,covdca);
446  }
447  if(!ok){
448  dca[0]=-999.;
449  dca[1]=-999.;
450  }
451  return ok;
452 }
453 
454 
455 
457  //Create the histograms
458 
459  if(!fHistograms) fHistograms = new TList();
460 
461  fHistograms->SetOwner(kTRUE);
462  fHistograms->SetName("trackCuts");
463 
464  fhPhi = new TH2F(Form("phi_%s", GetName()), Form("phi_%s", GetTitle()), 5, -0.5, 4.5, 32, 0, TMath::TwoPi());
465  // TAxis * xax = fhPhi->GetXaxis();
466  // for(Int_t i = 0; i < kNCuts; i++){
467  // xax->SetBinLabel(xax->FindFixBin(i), fgkCutNames[i]);
468  // }
469  fHistograms->Add(fhPhi);
470 
471 
472  fhEtaPhi = new TH2F(Form("etaphi_%s",GetName()), Form("etaphi_%s", GetTitle()), 36, -0.9, 0.9, 32, 0, TMath::TwoPi());
473  fHistograms->Add(fhEtaPhi);
474 
475  // fhPt = new TH2F("pt", "pt", kNCuts+2, kPreCut -0.5, kNCuts + 0.5,
476  // 20, 0., 20.);
477  // xax = fhPt->GetXaxis();
478  // for(Int_t i = 0; i < kNCuts; i++){
479  // xax->SetBinLabel(xax->FindFixBin(i), fgkCutNames[i]);
480  // }
481  // fHistograms->Add(fhPt);
482 
483  // fhPhiPt = new TH2F("phipt", "phipt", 100, 0, 100, 64, 0, TMath::TwoPi());
484  //fHistograms->Add(fhPhiPt);
485 
486  fhdcaxyPt = new TH2F(Form("dcaxypt_%s", GetName()), Form("dcaxypt_%s", GetTitle()), 20, 0, 20, 50, -2.5, 2.5);
487  fHistograms->Add(fhdcaxyPt);
488 
489  fhdcazPt = new TH2F(Form("dcazpt_%s", GetName()), Form("dcazpt_%s", GetTitle()), 20, 0, 20, 70, -3.5, 3.5);
490  fHistograms->Add(fhdcazPt);
491 
492  fhdca = new TH2F(Form("dca_%s", GetName()), Form("dca_%s", GetTitle()), 70, -3.5, 3.5, 50, -2.5, 2.5);
493  fhdca->SetXTitle("dca z");
494  fhdca->SetYTitle("dca xy");
495 
496 
497  fHistograms->Add(fhdca);
498 
499  // fhnclpt = new TH2F("nclstpcvspt", "nclstpcvspt", 20, 0, 20, 50, 0, 100);
500  // fHistograms->Add(fhnclpt);
501 
502  // fhnclsfpt = new TH2F("nclsfpt", "nclsfpt", 20, 0, 20, 60, 0, 1.2);
503  // fHistograms->Add(fhnclsfpt);
504 
505 
506 
507  if (fkCreateTrackEff) {
508  const Double_t ptbins[23] = {0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4, 1.6, 1.8,
509  2.0, 2.25, 2.5, 2.75, 3.0, 3.25, 3.5, 3.75, 4.0, 4.25,
510  4.5, 4.75, 5.0};
511 
512  const Int_t bins[4] = { 12, 22, 36, 32};
513  const Double_t min[4] = { -12, 0.5, -.9, 0};
514  const Double_t max[4] = { 12, 5, .9, 2*TMath::Pi() };
515 
516  fhTrackEff = new THnF(Form("hTrackEff_%s", GetName()), "hTrackEff", 4, bins, min, max);
517  fhTrackEff->SetBinEdges(1, ptbins);
518  fHistograms->Add(fhTrackEff);
519  }
520 
521 
522  return fHistograms;
523 }
524 
525 void AliConversionTrackCuts::FillHistograms(Int_t cutIndex, const AliVTrack * track) {
526 
527  //Fill histograms
528  if(fhPhi) fhPhi->Fill(cutIndex, track->Phi());
529  // if(fhPt) fhPt->Fill(cutIndex, track->Pt());
530  //if(passed) fhPhiPt->Fill(track->Pt(), track->Phi());
531 
532 }
533 
534 void AliConversionTrackCuts::FillDCAHist(Float_t dcaz, Float_t dcaxy, const AliVTrack * track) {
535  if(fhdcaxyPt) fhdcaxyPt->Fill(track->Pt(), dcaxy);
536  if(fhdcazPt) fhdcazPt->Fill(track->Pt(), dcaz);
537  if(fhdca) fhdca->Fill(dcaz, dcaxy);
538 
539  if(fhTrackEff) {
540  Double_t val[4];
541  val[0] = fEvent->GetPrimaryVertex()->GetZ();
542  val[1] = track->Pt();
543  val[2] = track->Eta();
544  val[3] = track->Phi();
545 
546  fhTrackEff->Fill(val);
547  }
548 }
549 
550 
551 
552 
553 //_________________________________________________________________________________________________
555 {
556 //
557 // Print information on this cut
558 //
559 
560  // printf("Cut name : %s \n", GetName());
561  // printf("Kink daughters are : %s \n", (fRejectKinkDaughters ? "rejected" : "accepted"));
562  // printf("TPC requirements : clusters/findable %f, min. cluster = %d, max chi2 = %f, %s require refit\n", fTPCClusOverFindable, fTPCminNClusters, fTPCmaxChi2, (fRequireTPCRefit) ? "" : "Don't");
563  // printf("ITS requirements : min. cluster = %d (all), %d (SPD), max chi2 = %f \n", fITSminNClusters, fSPDminNClusters, fITSmaxChi2);
564  // printf("DCA z cut : fixed to %f cm \n", fDCAZmax);
565  // printf("DCA xy cut : fixed to %f cm \n", fDCAXYmax)
566  ;
567 }
568 
569 //_________________________________________________________________________________________________
570 
572  AliAODTrack * aodtrack = dynamic_cast<AliAODTrack*>(object);
573  if (aodtrack) {
574  return AcceptTrack(aodtrack);
575  } else {
576  AliESDtrack * track = dynamic_cast<AliESDtrack*>(object);
577  if (track)
578  return AcceptTrack(track);
579  }
580 
581 return kFALSE;
582 }
583 
584 
585 
586 
Bool_t IsSelected(TObject *object)
double Double_t
Definition: External.C:58
Definition: External.C:236
const char * title
Definition: MakeQAPdf.C:27
AliESDtrackCuts * fEsdTrackCuts
int Int_t
Definition: External.C:63
Bool_t AcceptTrack(const AliAODTrack *track)
float Float_t
Definition: External.C:68
void FillDCAHist(Float_t dcaz, Float_t dcaxy, const AliVTrack *track)
void FillHistograms(Int_t cutIndex, const AliVTrack *track)
AliESDtrackCuts * fEsdTrackCutsExtra2
virtual void Print(const Option_t *option="") const
Bool_t GetDCA(const AliAODTrack *track, Double_t dca[2])
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
AliESDtrackCuts * fEsdTrackCutsExtra1