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