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