PWG0/esdTrackCuts/AliESDtrackCuts.cxx

   1 #include "AliESDtrackCuts.h"
   2
   3 #include <Riostream.h>
   4
   5 //____________________________________________________________________
   6 ClassImp(AliESDtrackCuts);
   7
   8 //____________________________________________________________________
   9 AliESDtrackCuts::AliESDtrackCuts() {
  10
  11   //##############################################################################
  12   // setting default cuts
  13
  14   SetMinNClustersTPC();
  15   SetMinNClustersITS();
  16   SetMaxChi2PerClusterTPC();
  17   SetMaxChi2PerClusterITS();
  18   SetMaxCovDiagonalElements();
  19   SetRequireTPCRefit();
  20   SetRequireITSRefit();
  21   SetAcceptKingDaughters();
  22   SetMinNsigmaToVertex();
  23   SetRequireSigmaToVertex();
  24   SetPRange();
  25   SetPtRange();
  26   SetPxRange();
  27   SetPyRange();
  28   SetPzRange();
  29   SetEtaRange();
  30   SetRapRange();
  31
  32   SetHistogramsOn();
  33
  34   // set the cut names
  35   fCutNames[0]  = "require TPC refit";
  36   fCutNames[1]  = "require ITS refit";
  37   fCutNames[2]  = "n clusters TPC";
  38   fCutNames[3]  = "n clusters ITS";
  39   fCutNames[4]  = "#Chi^{2}/clusters TPC";
  40   fCutNames[5]  = "#Chi^{2}/clusters ITS";
  41   fCutNames[6]  = "cov 11";
  42   fCutNames[7]  = "cov 22";
  43   fCutNames[8]  = "cov 33";
  44   fCutNames[9]  = "cov 44";
  45   fCutNames[10] = "cov 55";
  46   fCutNames[11] = "trk-to-vtx";
  47   fCutNames[12] = "trk-to-vtx failed";
  48   fCutNames[13] = "kink daughters";
  49
  50   fCutNames[14] = "p";
  51   fCutNames[15] = "p_{T}";
  52   fCutNames[16] = "p_{x}";
  53   fCutNames[17] = "p_{y}";
  54   fCutNames[18] = "p_{z}";
  55   fCutNames[19] = "y";
  56   fCutNames[20] = "eta";
  57
  58 }
  59
  60 //____________________________________________________________________
  61 Bool_t
  62 AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack, AliESDVertex* esdVtx, Double_t field) {
  63   //
  64   // re-calculate the track-to-vertex
  65   esdTrack->RelateToVertex(esdVtx, field, 1e99);
  66
  67   return AcceptTrack(esdTrack);
  68 }
  69
  70 //____________________________________________________________________
  71 Bool_t
  72 AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack, Double_t* vtx, Double_t* vtx_res, Double_t field) {
  73
  74   AliESDVertex* esdVtx = new AliESDVertex(vtx, vtx_res,"new vertex");
  75   esdTrack->RelateToVertex(esdVtx, field, 1e99);
  76   return AcceptTrack(esdTrack);
  77 }
  78
  79 //____________________________________________________________________
  80 Bool_t
  81 AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
  82   //
  83   // figure out if the tracks survives all the track cuts defined
  84   //
  85
  86   UInt_t status = esdTrack->GetStatus();
  87
  88   // getting quality parameters from the ESD track
  89   Int_t nClustersITS = esdTrack->GetITSclusters(fIdxInt);
  90   Int_t nClustersTPC = esdTrack->GetTPCclusters(fIdxInt);
  91
  92   Float_t chi2PerClusterITS = -1;
  93   Float_t chi2PerClusterTPC = -1;
  94   if (nClustersITS!=0)
  95     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
  96   if (nClustersTPC!=0)
  97     chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
  98
  99   Double_t extCov[15];
 100   esdTrack->GetExternalCovariance(extCov);
 101
 102   // getting the track to vertex parameters
 103   Float_t b[2];
 104   Float_t bRes[2];
 105   Float_t bCov[3];
 106   esdTrack->GetImpactParameters(b,bCov);
 107   if (bCov[0]<=0 || bCov[2]<=0) {
 108     AliDebug(1, "Estimated b resolution zero!");
 109     bCov[0]=0; bCov[1]=0;
 110   }
 111   bRes[0] = TMath::Sqrt(bCov[0]);
 112   bRes[1] = TMath::Sqrt(bCov[2]);
 113
 114   // FIX !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 115   //
 116   // this is not correct - it will not give n sigma!!!
 117   //
 118   Float_t nSigmaToVertex = -1;
 119   if (bRes[0]!=0 && bRes[1]!=0)
 120     nSigmaToVertex = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 121
 122   // getting the kinematic variables of the track
 123   // (assuming the mass is known)
 124   Double_t p[3];
 125   esdTrack->GetPxPyPz(p);
 126   Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
 127   Float_t pt       = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2));
 128   Float_t energy   = TMath::Sqrt(TMath::Power(esdTrack->GetMass(),2) + TMath::Power(momentum,2));
 129
 130   //y-eta related calculations
 131   Float_t eta = -100.;
 132   Float_t y   = -100.;
 133   if((momentum != TMath::Abs(p[2]))&&(momentum != 0))
 134     eta = 0.5*TMath::Log((momentum + p[2])/(momentum - p[2]));
 135   if((energy != TMath::Abs(p[2]))&&(momentum != 0))
 136     y = 0.5*TMath::Log((energy + p[2])/(energy - p[2]));
 137
 138
 139   //########################################################################
 140   // cut the track?
 141
 142   Bool_t cuts[fNCuts];
 143   for (Int_t i=0; i<fNCuts; i++) cuts[i]=kFALSE;
 144
 145   // track quality cuts
 146   if (fCut_RequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
 147     cuts[0]=kTRUE;
 148   if (fCut_RequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
 149     cuts[1]=kTRUE;
 150   if (nClustersTPC<fCut_MinNClusterTPC)
 151     cuts[2]=kTRUE;
 152   if (nClustersITS<fCut_MinNClusterITS)
 153     cuts[3]=kTRUE;
 154   if (chi2PerClusterTPC>fCut_MaxChi2PerClusterTPC)
 155     cuts[4]=kTRUE;
 156   if (chi2PerClusterITS>fCut_MaxChi2PerClusterITS)
 157     cuts[5]=kTRUE;
 158   if (extCov[0]  > fCut_MaxC11)
 159     cuts[6]=kTRUE;
 160   if (extCov[2]  > fCut_MaxC22)
 161     cuts[7]=kTRUE;
 162   if (extCov[5]  > fCut_MaxC33)
 163     cuts[8]=kTRUE;
 164   if (extCov[9]  > fCut_MaxC44)
 165     cuts[9]=kTRUE;
 166   if (extCov[14]  > fCut_MaxC55)
 167     cuts[10]=kTRUE;
 168   if (nSigmaToVertex > fCut_NsigmaToVertex)
 169     cuts[11] = kTRUE;
 170   // if n sigma could not be calculated
 171   if (nSigmaToVertex<0 && fCut_SigmaToVertexRequired)
 172     cuts[12]=kTRUE;
 173   if (!fCut_AcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
 174     cuts[13]=kTRUE;
 175   // track kinematics cut
 176   if((momentum < fPMin) || (momentum > fPMax))
 177     cuts[14]=kTRUE;
 178   if((pt < fPtMin) || (pt > fPtMax))
 179     cuts[15] = kTRUE;
 180   if((p[0] < fPxMin) || (p[0] > fPxMax))
 181     cuts[16] = kTRUE;
 182   if((p[1] < fPyMin) || (p[1] > fPyMax))
 183     cuts[17] = kTRUE;
 184   if((p[2] < fPzMin) || (p[2] > fPzMax))
 185     cuts[18] = kTRUE;
 186   if((eta < fEtaMin) || (eta > fEtaMax))
 187     cuts[19] = kTRUE;
 188   if((y < fRapMin) || (y > fRapMax))
 189     cuts[20] = kTRUE;
 190
 191   Bool_t cut=kFALSE;
 192   for (Int_t i=0; i<fNCuts; i++)
 193     if (cuts[i]) cut = kTRUE;
 194
 195   //########################################################################
 196   // filling histograms
 197   if (fHistogramsOn) {
 198     hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin("n tracks")));
 199
 200     if (cut)
 201       hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin("n cut tracks")));
 202
 203     for (Int_t i=0; i<fNCuts; i++) {
 204       if (cuts[i])
 205         hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin(fCutNames[i])));
 206
 207       for (Int_t j=i; j<fNCuts; j++) {
 208         if (cuts[i] && cuts[j]) {
 209           Float_t x = hCutCorrelation->GetXaxis()->GetBinCenter(hCutCorrelation->GetXaxis()->FindBin(fCutNames[i]));
 210           Float_t y = hCutCorrelation->GetYaxis()->GetBinCenter(hCutCorrelation->GetYaxis()->FindBin(fCutNames[j]));
 211           hCutCorrelation->Fill(x,y);
 212         }
 213       }
 214     }
 215
 216
 217     hNClustersITS[0]->Fill(nClustersITS);
 218     hNClustersTPC[0]->Fill(nClustersTPC);
 219     hChi2PerClusterITS[0]->Fill(chi2PerClusterITS);
 220     hChi2PerClusterTPC[0]->Fill(chi2PerClusterTPC);
 221
 222     hC11[0]->Fill(extCov[0]);
 223     hC22[0]->Fill(extCov[2]);
 224     hC33[0]->Fill(extCov[5]);
 225     hC44[0]->Fill(extCov[9]);
 226     hC55[0]->Fill(extCov[14]);
 227
 228     hDZ[0]->Fill(b[1]);
 229     hDXY[0]->Fill(b[0]);
 230     hDXYvsDZ[0]->Fill(b[1],b[0]);
 231
 232     if (bRes[0]!=0 && bRes[1]!=0) {
 233       hDZNormalized[0]->Fill(b[1]/bRes[1]);
 234       hDXYNormalized[0]->Fill(b[0]/bRes[0]);
 235       hDXYvsDZNormalized[0]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
 236     }
 237   }
 238
 239   //########################################################################
 240   // cut the track!
 241   if (cut) return kFALSE;
 242
 243   //########################################################################
 244   // filling histograms after cut
 245   if (fHistogramsOn) {
 246     hNClustersITS[1]->Fill(nClustersITS);
 247     hNClustersTPC[1]->Fill(nClustersTPC);
 248     hChi2PerClusterITS[1]->Fill(chi2PerClusterITS);
 249     hChi2PerClusterTPC[1]->Fill(chi2PerClusterTPC);
 250
 251     hC11[1]->Fill(extCov[0]);
 252     hC22[1]->Fill(extCov[2]);
 253     hC33[1]->Fill(extCov[5]);
 254     hC44[1]->Fill(extCov[9]);
 255     hC55[1]->Fill(extCov[14]);
 256
 257     hDZ[1]->Fill(b[1]);
 258     hDXY[1]->Fill(b[0]);
 259     hDXYvsDZ[1]->Fill(b[1],b[0]);
 260
 261     hDZNormalized[1]->Fill(b[1]/bRes[1]);
 262     hDXYNormalized[1]->Fill(b[0]/bRes[0]);
 263     hDXYvsDZNormalized[1]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
 264   }
 265
 266   return kTRUE;
 267 }
 268
 269 //____________________________________________________________________
 270 TObjArray*
 271 AliESDtrackCuts::GetAcceptedTracks(AliESD* esd) {
 272
 273   // returns an array of all tracks that pass the cuts
 274   fAcceptedTracks->Clear();
 275
 276   // loop over esd tracks
 277   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
 278     AliESDtrack* track = esd->GetTrack(iTrack);
 279
 280     if(AcceptTrack(track)) fAcceptedTracks->Add(track);
 281   }
 282
 283   return fAcceptedTracks;
 284 }
 285
 286 //____________________________________________________________________
 287 void
 288 AliESDtrackCuts::DefineHistograms(Int_t color) {
 289
 290   fHistogramsOn=kTRUE;
 291
 292   //###################################################################################
 293   // defining histograms
 294
 295   hCutStatistics = new TH1F("cut_statistics","cut statistics",fNCuts+4,-0.5,fNCuts+3.5);
 296
 297   hCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
 298   hCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
 299
 300   hCutCorrelation = new TH2F("cut_correlation","cut correlation",fNCuts,-0.5,fNCuts-0.5,fNCuts,-0.5,fNCuts-0.5);;
 301
 302   for (Int_t i=0; i<fNCuts; i++) {
 303     hCutStatistics->GetXaxis()->SetBinLabel(i+4,fCutNames[i]);
 304     hCutCorrelation->GetXaxis()->SetBinLabel(i+1,fCutNames[i]);
 305     hCutCorrelation->GetYaxis()->SetBinLabel(i+1,fCutNames[i]);
 306   }
 307
 308   hCutStatistics  ->SetLineColor(color);
 309   hCutCorrelation ->SetLineColor(color);
 310   hCutStatistics  ->SetLineWidth(2);
 311   hCutCorrelation ->SetLineWidth(2);
 312
 313
 314   hNClustersITS        = new TH1F*[2];
 315   hNClustersTPC        = new TH1F*[2];
 316   hChi2PerClusterITS   = new TH1F*[2];
 317   hChi2PerClusterTPC   = new TH1F*[2];
 318
 319   hC11                 = new TH1F*[2];
 320   hC22                 = new TH1F*[2];
 321   hC33                 = new TH1F*[2];
 322   hC44                 = new TH1F*[2];
 323   hC55                 = new TH1F*[2];
 324
 325   hDXY                 = new TH1F*[2];
 326   hDZ                  = new TH1F*[2];
 327   hDXYvsDZ             = new TH2F*[2];
 328
 329   hDXYNormalized       = new TH1F*[2];
 330   hDZNormalized        = new TH1F*[2];
 331   hDXYvsDZNormalized   = new TH2F*[2];
 332
 333
 334   Char_t str[256];
 335   for (Int_t i=0; i<2; i++) {
 336     if (i==0) sprintf(str," ");
 337     else sprintf(str,"_cut");
 338
 339     hNClustersITS[i]        = new TH1F(Form("nClustersITS%s",str),"",8,-0.5,7.5);
 340     hNClustersTPC[i]        = new TH1F(Form("nClustersTPC%s",str),"",165,-0.5,164.5);
 341     hChi2PerClusterITS[i]   = new TH1F(Form("chi2PerClusterITS%s",str),"",500,0,10);
 342     hChi2PerClusterTPC[i]   = new TH1F(Form("chi2PerClusterTPC%s",str),"",500,0,10);
 343
 344     hC11[i]                 = new TH1F(Form("covMatrixDiagonal11%s",str),"",1000,0,5);
 345     hC22[i]                 = new TH1F(Form("covMatrixDiagonal22%s",str),"",1000,0,5);
 346     hC33[i]                 = new TH1F(Form("covMatrixDiagonal33%s",str),"",1000,0,0.5);
 347     hC44[i]                 = new TH1F(Form("covMatrixDiagonal44%s",str),"",1000,0,5);
 348     hC55[i]                 = new TH1F(Form("covMatrixDiagonal55%s",str),"",1000,0,5);
 349
 350     hDXY[i]                 = new TH1F(Form("dXY%s",str),"",500,-10,10);
 351     hDZ[i]                  = new TH1F(Form("dZ%s",str),"",500,-10,10);
 352     hDXYvsDZ[i]             = new TH2F(Form("dXYvsDZ%s",str),"",200,-10,10,200,-10,10);
 353
 354     hDXYNormalized[i]       = new TH1F(Form("dXYNormalized%s",str),"",500,-10,10);
 355     hDZNormalized[i]        = new TH1F(Form("dZNormalized%s",str),"",500,-10,10);
 356     hDXYvsDZNormalized[i]   = new TH2F(Form("dXYvsDZNormalized%s",str),"",200,-10,10,200,-10,10);
 357
 358
 359     hNClustersITS[i]        ->SetXTitle("n ITS clusters");
 360     hNClustersTPC[i]        ->SetXTitle("n TPC clusters");
 361     hChi2PerClusterITS[i]   ->SetXTitle("#Chi^{2} per ITS cluster");
 362     hChi2PerClusterTPC[i]   ->SetXTitle("#Chi^{2} per TPC cluster");
 363
 364     hC11[i]                 ->SetXTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
 365     hC22[i]                 ->SetXTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
 366     hC33[i]                 ->SetXTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
 367     hC44[i]                 ->SetXTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
 368     hC55[i]                 ->SetXTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
 369
 370     hDXY[i]                 ->SetXTitle("transverse impact parameter");
 371     hDZ[i]                  ->SetXTitle("longitudinal impact parameter");
 372     hDXYvsDZ[i]             ->SetXTitle("longitudinal impact parameter");
 373     hDXYvsDZ[i]             ->SetYTitle("transverse impact parameter");
 374
 375     hDXYNormalized[i]       ->SetXTitle("normalized trans impact par");
 376     hDZNormalized[i]        ->SetXTitle("normalized long impact par");
 377     hDXYvsDZNormalized[i]   ->SetXTitle("normalized long impact par");
 378     hDXYvsDZNormalized[i]   ->SetYTitle("normalized trans impact par");
 379
 380     hNClustersITS[i]        ->SetLineColor(color);   hNClustersITS[i]        ->SetLineWidth(2);
 381     hNClustersTPC[i]        ->SetLineColor(color);   hNClustersTPC[i]        ->SetLineWidth(2);
 382     hChi2PerClusterITS[i]   ->SetLineColor(color);   hChi2PerClusterITS[i]   ->SetLineWidth(2);
 383     hChi2PerClusterTPC[i]   ->SetLineColor(color);   hChi2PerClusterTPC[i]   ->SetLineWidth(2);
 384
 385     hC11[i]                 ->SetLineColor(color);   hC11[i]                 ->SetLineWidth(2);
 386     hC22[i]                 ->SetLineColor(color);   hC22[i]                 ->SetLineWidth(2);
 387     hC33[i]                 ->SetLineColor(color);   hC33[i]                 ->SetLineWidth(2);
 388     hC44[i]                 ->SetLineColor(color);   hC44[i]                 ->SetLineWidth(2);
 389     hC55[i]                 ->SetLineColor(color);   hC55[i]                 ->SetLineWidth(2);
 390
 391     hDXY[i]                 ->SetLineColor(color);   hDXY[i]                 ->SetLineWidth(2);
 392     hDZ[i]                  ->SetLineColor(color);   hDZ[i]                  ->SetLineWidth(2);
 393
 394     hDXYNormalized[i]       ->SetLineColor(color);   hDXYNormalized[i]       ->SetLineWidth(2);
 395     hDZNormalized[i]        ->SetLineColor(color);   hDZNormalized[i]        ->SetLineWidth(2);
 396
 397   }
 398 }
 399
 400 //____________________________________________________________________
 401 void
 402 AliESDtrackCuts::Print() {
 403
 404   AliInfo("AliESDtrackCuts...");
 405 }
 406
 407
 408 //____________________________________________________________________
 409 void
 410 AliESDtrackCuts::SaveHistograms(Char_t* dir) {
 411
 412   if (!fHistogramsOn) {
 413     AliDebug(0, "Histograms not on - cannot save histograms!!!");
 414     return;
 415   }
 416
 417   gDirectory->mkdir(dir);
 418   gDirectory->cd(dir);
 419
 420   gDirectory->mkdir("before_cuts");
 421   gDirectory->mkdir("after_cuts");
 422
 423   hCutStatistics->Write();
 424   hCutCorrelation->Write();
 425
 426   for (Int_t i=0; i<2; i++) {
 427     if (i==0)
 428       gDirectory->cd("before_cuts");
 429     else
 430       gDirectory->cd("after_cuts");
 431
 432     hNClustersITS[i]        ->Write();
 433     hNClustersTPC[i]        ->Write();
 434     hChi2PerClusterITS[i]   ->Write();
 435     hChi2PerClusterTPC[i]   ->Write();
 436
 437     hC11[i]                 ->Write();
 438     hC22[i]                 ->Write();
 439     hC33[i]                 ->Write();
 440     hC44[i]                 ->Write();
 441     hC55[i]                 ->Write();
 442
 443     hDXY[i]                 ->Write();
 444     hDZ[i]                  ->Write();
 445     hDXYvsDZ[i]             ->Write();
 446
 447     hDXYNormalized[i]       ->Write();
 448     hDZNormalized[i]        ->Write();
 449     hDXYvsDZNormalized[i]   ->Write();
 450
 451     gDirectory->cd("../");
 452   }
 453
 454   gDirectory->cd("../");
 455 }
 456
 457
 458