PWG2/AliESDtrackCuts.cxx

   1 #include "AliESDtrackCuts.h"
   2
   3
   4 #include <AliESDtrack.h>
   5 #include <AliESD.h>
   6 #include <AliLog.h>
   7 #include <TTree.h>
   8 #include <TFile.h>
   9
  10 //____________________________________________________________________
  11 ClassImp(AliESDtrackCuts)
  12
  13 // Cut names
  14 const Char_t* AliESDtrackCuts::fgkCutNames[kNCuts] = {
  15  "require TPC refit",
  16  "require ITS refit",
  17  "n clusters TPC",
  18  "n clusters ITS",
  19  "#Chi^{2}/clusters TPC",
  20  "#Chi^{2}/clusters ITS",
  21  "cov 11",
  22  "cov 22",
  23  "cov 33",
  24  "cov 44",
  25  "cov 55",
  26  "trk-to-vtx",
  27  "trk-to-vtx failed",
  28  "kink daughters",
  29  "p",
  30  "p_{T}",
  31  "p_{x}",
  32  "p_{y}",
  33  "p_{z}",
  34  "y",
  35  "eta"
  36 };
  37
  38 //____________________________________________________________________
  39 AliESDtrackCuts::AliESDtrackCuts() : TNamed(),
  40   fCutMinNClusterTPC(0),
  41   fCutMinNClusterITS(0),
  42   fCutMaxChi2PerClusterTPC(0),
  43   fCutMaxChi2PerClusterITS(0),
  44   fCutMaxC11(0),
  45   fCutMaxC22(0),
  46   fCutMaxC33(0),
  47   fCutMaxC44(0),
  48   fCutMaxC55(0),
  49   fCutAcceptKinkDaughters(0),
  50   fCutRequireTPCRefit(0),
  51   fCutRequireITSRefit(0),
  52   fCutNsigmaToVertex(0),
  53   fCutSigmaToVertexRequired(0),
  54   fPMin(0),
  55   fPMax(0),
  56   fPtMin(0),
  57   fPtMax(0),
  58   fPxMin(0),
  59   fPxMax(0),
  60   fPyMin(0),
  61   fPyMax(0),
  62   fPzMin(0),
  63   fPzMax(0),
  64   fEtaMin(0),
  65   fEtaMax(0),
  66   fRapMin(0),
  67   fRapMax(0),
  68   fHistogramsOn(0),
  69   ffDTheoretical(0),
  70   fhCutStatistics(0),
  71   fhCutCorrelation(0)
  72 {
  73   //
  74   // default constructor
  75   //
  76
  77   Init();
  78 }
  79
  80 //____________________________________________________________________
  81 AliESDtrackCuts::AliESDtrackCuts(Char_t* name, Char_t* title) : TNamed(name,title),
  82   fCutMinNClusterTPC(0),
  83   fCutMinNClusterITS(0),
  84   fCutMaxChi2PerClusterTPC(0),
  85   fCutMaxChi2PerClusterITS(0),
  86   fCutMaxC11(0),
  87   fCutMaxC22(0),
  88   fCutMaxC33(0),
  89   fCutMaxC44(0),
  90   fCutMaxC55(0),
  91   fCutAcceptKinkDaughters(0),
  92   fCutRequireTPCRefit(0),
  93   fCutRequireITSRefit(0),
  94   fCutNsigmaToVertex(0),
  95   fCutSigmaToVertexRequired(0),
  96   fPMin(0),
  97   fPMax(0),
  98   fPtMin(0),
  99   fPtMax(0),
 100   fPxMin(0),
 101   fPxMax(0),
 102   fPyMin(0),
 103   fPyMax(0),
 104   fPzMin(0),
 105   fPzMax(0),
 106   fEtaMin(0),
 107   fEtaMax(0),
 108   fRapMin(0),
 109   fRapMax(0),
 110   fHistogramsOn(0),
 111   fhCutStatistics(0),
 112   fhCutCorrelation(0)
 113 {
 114   //
 115   // constructor
 116   //
 117
 118   Init();
 119
 120   //##############################################################################
 121   // setting default cuts
 122   SetMinNClustersTPC();
 123   SetMinNClustersITS();
 124   SetMaxChi2PerClusterTPC();
 125   SetMaxChi2PerClusterITS();
 126   SetMaxCovDiagonalElements();
 127   SetRequireTPCRefit();
 128   SetRequireITSRefit();
 129   SetAcceptKingDaughters();
 130   SetMinNsigmaToVertex();
 131   SetRequireSigmaToVertex();
 132   SetPRange();
 133   SetPtRange();
 134   SetPxRange();
 135   SetPyRange();
 136   SetPzRange();
 137   SetEtaRange();
 138   SetRapRange();
 139
 140   SetHistogramsOn();
 141 }
 142
 143 //_____________________________________________________________________________
 144 AliESDtrackCuts::AliESDtrackCuts(const AliESDtrackCuts &c) : TNamed(c),
 145   fCutMinNClusterTPC(0),
 146   fCutMinNClusterITS(0),
 147   fCutMaxChi2PerClusterTPC(0),
 148   fCutMaxChi2PerClusterITS(0),
 149   fCutMaxC11(0),
 150   fCutMaxC22(0),
 151   fCutMaxC33(0),
 152   fCutMaxC44(0),
 153   fCutMaxC55(0),
 154   fCutAcceptKinkDaughters(0),
 155   fCutRequireTPCRefit(0),
 156   fCutRequireITSRefit(0),
 157   fCutNsigmaToVertex(0),
 158   fCutSigmaToVertexRequired(0),
 159   fPMin(0),
 160   fPMax(0),
 161   fPtMin(0),
 162   fPtMax(0),
 163   fPxMin(0),
 164   fPxMax(0),
 165   fPyMin(0),
 166   fPyMax(0),
 167   fPzMin(0),
 168   fPzMax(0),
 169   fEtaMin(0),
 170   fEtaMax(0),
 171   fRapMin(0),
 172   fRapMax(0),
 173   fHistogramsOn(0),
 174   ffDTheoretical(0),
 175   fhCutStatistics(0),
 176   fhCutCorrelation(0)
 177 {
 178   //
 179   // copy constructor
 180   //
 181
 182   ((AliESDtrackCuts &) c).Copy(*this);
 183 }
 184
 185 AliESDtrackCuts::~AliESDtrackCuts()
 186 {
 187   //
 188   // destructor
 189   //
 190
 191   for (Int_t i=0; i<2; i++) {
 192
 193     if (fhNClustersITS[i])
 194       delete fhNClustersITS[i];
 195     if (fhNClustersTPC[i])
 196       delete fhNClustersTPC[i];
 197     if (fhChi2PerClusterITS[i])
 198       delete fhChi2PerClusterITS[i];
 199     if (fhChi2PerClusterTPC[i])
 200       delete fhChi2PerClusterTPC[i];
 201     if (fhC11[i])
 202       delete fhC11[i];
 203     if (fhC22[i])
 204       delete fhC22[i];
 205     if (fhC33[i])
 206       delete fhC33[i];
 207     if (fhC44[i])
 208       delete fhC44[i];
 209     if (fhC55[i])
 210     delete fhC55[i];
 211
 212     if (fhDXY[i])
 213       delete fhDXY[i];
 214     if (fhDZ[i])
 215       delete fhDZ[i];
 216     if (fhDXYvsDZ[i])
 217       delete fhDXYvsDZ[i];
 218
 219     if (fhDXYNormalized[i])
 220       delete fhDXYNormalized[i];
 221     if (fhDZNormalized[i])
 222       delete fhDZNormalized[i];
 223     if (fhDXYvsDZNormalized[i])
 224       delete fhDXYvsDZNormalized[i];
 225     if (fhNSigmaToVertex[i])
 226       delete fhNSigmaToVertex[i];
 227   }
 228
 229   if (ffDTheoretical)
 230     delete ffDTheoretical;
 231
 232   if (fhCutStatistics)
 233     delete fhCutStatistics;
 234   if (fhCutCorrelation)
 235     delete fhCutCorrelation;
 236 }
 237
 238 void AliESDtrackCuts::Init()
 239 {
 240   //
 241   // sets everything to zero
 242   //
 243
 244   fCutMinNClusterTPC = 0;
 245   fCutMinNClusterITS = 0;
 246
 247   fCutMaxChi2PerClusterTPC = 0;
 248   fCutMaxChi2PerClusterITS = 0;
 249
 250   fCutMaxC11 = 0;
 251   fCutMaxC22 = 0;
 252   fCutMaxC33 = 0;
 253   fCutMaxC44 = 0;
 254   fCutMaxC55 = 0;
 255
 256   fCutAcceptKinkDaughters = 0;
 257   fCutRequireTPCRefit = 0;
 258   fCutRequireITSRefit = 0;
 259
 260   fCutNsigmaToVertex = 0;
 261   fCutSigmaToVertexRequired = 0;
 262
 263   fPMin = 0;
 264   fPMax = 0;
 265   fPtMin = 0;
 266   fPtMax = 0;
 267   fPxMin = 0;
 268   fPxMax = 0;
 269   fPyMin = 0;
 270   fPyMax = 0;
 271   fPzMin = 0;
 272   fPzMax = 0;
 273   fEtaMin = 0;
 274   fEtaMax = 0;
 275   fRapMin = 0;
 276   fRapMax = 0;
 277
 278   fHistogramsOn = kFALSE;
 279
 280   for (Int_t i=0; i<2; ++i)
 281   {
 282     fhNClustersITS[i] = 0;
 283     fhNClustersTPC[i] = 0;
 284
 285     fhChi2PerClusterITS[i] = 0;
 286     fhChi2PerClusterTPC[i] = 0;
 287
 288     fhC11[i] = 0;
 289     fhC22[i] = 0;
 290     fhC33[i] = 0;
 291     fhC44[i] = 0;
 292     fhC55[i] = 0;
 293
 294     fhDXY[i] = 0;
 295     fhDZ[i] = 0;
 296     fhDXYvsDZ[i] = 0;
 297
 298     fhDXYNormalized[i] = 0;
 299     fhDZNormalized[i] = 0;
 300     fhDXYvsDZNormalized[i] = 0;
 301     fhNSigmaToVertex[i] = 0;
 302   }
 303   ffDTheoretical = 0;
 304
 305   fhCutStatistics = 0;
 306   fhCutCorrelation = 0;
 307 }
 308
 309 //_____________________________________________________________________________
 310 AliESDtrackCuts &AliESDtrackCuts::operator=(const AliESDtrackCuts &c)
 311 {
 312   //
 313   // Assignment operator
 314   //
 315
 316   if (this != &c) ((AliESDtrackCuts &) c).Copy(*this);
 317   return *this;
 318 }
 319
 320 //_____________________________________________________________________________
 321 void AliESDtrackCuts::Copy(TObject &c) const
 322 {
 323   //
 324   // Copy function
 325   //
 326
 327   AliESDtrackCuts& target = (AliESDtrackCuts &) c;
 328
 329   target.Init();
 330
 331   target.fCutMinNClusterTPC = fCutMinNClusterTPC;
 332   target.fCutMinNClusterITS = fCutMinNClusterITS;
 333
 334   target.fCutMaxChi2PerClusterTPC = fCutMaxChi2PerClusterTPC;
 335   target.fCutMaxChi2PerClusterITS = fCutMaxChi2PerClusterITS;
 336
 337   target.fCutMaxC11 = fCutMaxC11;
 338   target.fCutMaxC22 = fCutMaxC22;
 339   target.fCutMaxC33 = fCutMaxC33;
 340   target.fCutMaxC44 = fCutMaxC44;
 341   target.fCutMaxC55 = fCutMaxC55;
 342
 343   target.fCutAcceptKinkDaughters = fCutAcceptKinkDaughters;
 344   target.fCutRequireTPCRefit = fCutRequireTPCRefit;
 345   target.fCutRequireITSRefit = fCutRequireITSRefit;
 346
 347   target.fCutNsigmaToVertex = fCutNsigmaToVertex;
 348   target.fCutSigmaToVertexRequired = fCutSigmaToVertexRequired;
 349
 350   target.fPMin = fPMin;
 351   target.fPMax = fPMax;
 352   target.fPtMin = fPtMin;
 353   target.fPtMax = fPtMax;
 354   target.fPxMin = fPxMin;
 355   target.fPxMax = fPxMax;
 356   target.fPyMin = fPyMin;
 357   target.fPyMax = fPyMax;
 358   target.fPzMin = fPzMin;
 359   target.fPzMax = fPzMax;
 360   target.fEtaMin = fEtaMin;
 361   target.fEtaMax = fEtaMax;
 362   target.fRapMin = fRapMin;
 363   target.fRapMax = fRapMax;
 364
 365   target.fHistogramsOn = fHistogramsOn;
 366
 367   for (Int_t i=0; i<2; ++i)
 368   {
 369     if (fhNClustersITS[i]) target.fhNClustersITS[i] = (TH1F*) fhNClustersITS[i]->Clone();
 370     if (fhNClustersTPC[i]) target.fhNClustersTPC[i] = (TH1F*) fhNClustersTPC[i]->Clone();
 371
 372     if (fhChi2PerClusterITS[i]) target.fhChi2PerClusterITS[i] = (TH1F*) fhChi2PerClusterITS[i]->Clone();
 373     if (fhChi2PerClusterTPC[i]) target.fhChi2PerClusterTPC[i] = (TH1F*) fhChi2PerClusterTPC[i]->Clone();
 374
 375     if (fhC11[i]) target.fhC11[i] = (TH1F*) fhC11[i]->Clone();
 376     if (fhC22[i]) target.fhC22[i] = (TH1F*) fhC22[i]->Clone();
 377     if (fhC33[i]) target.fhC33[i] = (TH1F*) fhC33[i]->Clone();
 378     if (fhC44[i]) target.fhC44[i] = (TH1F*) fhC44[i]->Clone();
 379     if (fhC55[i]) target.fhC55[i] = (TH1F*) fhC55[i]->Clone();
 380
 381     if (fhDXY[i]) target.fhDXY[i] = (TH1F*) fhDXY[i]->Clone();
 382     if (fhDZ[i]) target.fhDZ[i] = (TH1F*) fhDZ[i]->Clone();
 383     if (fhDXYvsDZ[i]) target.fhDXYvsDZ[i] = (TH2F*) fhDXYvsDZ[i]->Clone();
 384
 385     if (fhDXYNormalized[i]) target.fhDXYNormalized[i] = (TH1F*) fhDXYNormalized[i]->Clone();
 386     if (fhDZNormalized[i]) target.fhDZNormalized[i] = (TH1F*) fhDZNormalized[i]->Clone();
 387     if (fhDXYvsDZNormalized[i]) target.fhDXYvsDZNormalized[i] = (TH2F*) fhDXYvsDZNormalized[i]->Clone();
 388     if (fhNSigmaToVertex[i]) target.fhNSigmaToVertex[i] = (TH1F*) fhNSigmaToVertex[i]->Clone();
 389   }
 390   if (ffDTheoretical) target.ffDTheoretical = (TF1*) ffDTheoretical->Clone();
 391
 392   if (fhCutStatistics) target.fhCutStatistics = (TH1F*) fhCutStatistics->Clone();
 393   if (fhCutCorrelation) target.fhCutCorrelation = (TH2F*) fhCutCorrelation->Clone();
 394
 395   TNamed::Copy(c);
 396 }
 397
 398 //_____________________________________________________________________________
 399 Long64_t AliESDtrackCuts::Merge(TCollection* list) {
 400   // Merge a list of AliESDtrackCuts objects with this (needed for PROOF)
 401   // Returns the number of merged objects (including this)
 402
 403   if (!list)
 404     return 0;
 405
 406   if (list->IsEmpty())
 407     return 1;
 408
 409   if (!fHistogramsOn)
 410     return 0;
 411
 412   TIterator* iter = list->MakeIterator();
 413   TObject* obj;
 414
 415
 416   // collection of measured and generated histograms
 417   Int_t count = 0;
 418   while ((obj = iter->Next())) {
 419
 420     AliESDtrackCuts* entry = dynamic_cast<AliESDtrackCuts*>(obj);
 421     if (entry == 0)
 422       continue;
 423
 424     if (!entry->fHistogramsOn)
 425       continue;
 426
 427     for (Int_t i=0; i<2; i++) {
 428
 429       fhNClustersITS[i]      ->Add(entry->fhNClustersITS[i]     );
 430       fhNClustersTPC[i]      ->Add(entry->fhNClustersTPC[i]     );
 431
 432       fhChi2PerClusterITS[i] ->Add(entry->fhChi2PerClusterITS[i]);
 433       fhChi2PerClusterTPC[i] ->Add(entry->fhChi2PerClusterTPC[i]);
 434
 435       fhC11[i]               ->Add(entry->fhC11[i]              );
 436       fhC22[i]               ->Add(entry->fhC22[i]              );
 437       fhC33[i]               ->Add(entry->fhC33[i]              );
 438       fhC44[i]               ->Add(entry->fhC44[i]              );
 439       fhC55[i]               ->Add(entry->fhC55[i]              );
 440
 441       fhDXY[i]               ->Add(entry->fhDXY[i]              );
 442       fhDZ[i]                ->Add(entry->fhDZ[i]               );
 443       fhDXYvsDZ[i]           ->Add(entry->fhDXYvsDZ[i]          );
 444
 445       fhDXYNormalized[i]     ->Add(entry->fhDXYNormalized[i]    );
 446       fhDZNormalized[i]      ->Add(entry->fhDZNormalized[i]     );
 447       fhDXYvsDZNormalized[i] ->Add(entry->fhDXYvsDZNormalized[i]);
 448       fhNSigmaToVertex[i]    ->Add(entry->fhNSigmaToVertex[i]);
 449
 450     }
 451
 452     fhCutStatistics  ->Add(entry->fhCutStatistics);
 453     fhCutCorrelation ->Add(entry->fhCutCorrelation);
 454
 455     count++;
 456   }
 457
 458   return count+1;
 459 }
 460
 461
 462 //____________________________________________________________________
 463 Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)
 464 {
 465   // Calculates the number of sigma to the vertex.
 466
 467   Float_t b[2];
 468   Float_t bRes[2];
 469   Float_t bCov[3];
 470   esdTrack->GetImpactParameters(b,bCov);
 471   if (bCov[0]<=0 || bCov[2]<=0) {
 472     AliDebug(1, "Estimated b resolution lower or equal zero!");
 473     bCov[0]=0; bCov[2]=0;
 474   }
 475   bRes[0] = TMath::Sqrt(bCov[0]);
 476   bRes[1] = TMath::Sqrt(bCov[2]);
 477
 478   // -----------------------------------
 479   // How to get to a n-sigma cut?
 480   //
 481   // The accumulated statistics from 0 to d is
 482   //
 483   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 484   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 485   //
 486   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
 487   // Can this be expressed in a different way?
 488
 489   if (bRes[0] == 0 || bRes[1] ==0)
 490     return -1;
 491
 492   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 493
 494   // stupid rounding problem screws up everything:
 495   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 496   if (TMath::Exp(-d * d / 2) < 1e-10)
 497     return 1000;
 498
 499   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 500   return d;
 501 }
 502
 503 void AliESDtrackCuts::EnableNeededBranches(TTree* tree)
 504 {
 505   // enables the branches needed by AcceptTrack, for a list see comment of AcceptTrack
 506
 507   tree->SetBranchStatus("fTracks.fFlags", 1);
 508   tree->SetBranchStatus("fTracks.fITSncls", 1);
 509   tree->SetBranchStatus("fTracks.fTPCncls", 1);
 510   tree->SetBranchStatus("fTracks.fITSchi2", 1);
 511   tree->SetBranchStatus("fTracks.fTPCchi2", 1);
 512   tree->SetBranchStatus("fTracks.fC*", 1);
 513   tree->SetBranchStatus("fTracks.fD", 1);
 514   tree->SetBranchStatus("fTracks.fZ", 1);
 515   tree->SetBranchStatus("fTracks.fCdd", 1);
 516   tree->SetBranchStatus("fTracks.fCdz", 1);
 517   tree->SetBranchStatus("fTracks.fCzz", 1);
 518   tree->SetBranchStatus("fTracks.fP*", 1);
 519   tree->SetBranchStatus("fTracks.fR*", 1);
 520   tree->SetBranchStatus("fTracks.fKinkIndexes*", 1);
 521 }
 522
 523 //____________________________________________________________________
 524 Bool_t
 525 AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
 526   //
 527   // figure out if the tracks survives all the track cuts defined
 528   //
 529   // the different quality parameter and kinematic values are first
 530   // retrieved from the track. then it is found out what cuts the
 531   // track did not survive and finally the cuts are imposed.
 532
 533   // this function needs the following branches:
 534   // fTracks.fFlags
 535   // fTracks.fITSncls
 536   // fTracks.fTPCncls
 537   // fTracks.fITSchi2
 538   // fTracks.fTPCchi2
 539   // fTracks.fC   //GetExternalCovariance
 540   // fTracks.fD   //GetImpactParameters
 541   // fTracks.fZ   //GetImpactParameters
 542   // fTracks.fCdd //GetImpactParameters
 543   // fTracks.fCdz //GetImpactParameters
 544   // fTracks.fCzz //GetImpactParameters
 545   // fTracks.fP   //GetPxPyPz
 546   // fTracks.fR   //GetMass
 547   // fTracks.fP   //GetMass
 548   // fTracks.fKinkIndexes
 549
 550   UInt_t status = esdTrack->GetStatus();
 551
 552   // dummy array
 553   Int_t  fIdxInt[200];
 554
 555   // getting quality parameters from the ESD track
 556   Int_t nClustersITS = esdTrack->GetITSclusters(fIdxInt);
 557   Int_t nClustersTPC = esdTrack->GetTPCclusters(fIdxInt);
 558
 559
 560
 561   Float_t chi2PerClusterITS = -1;
 562   Float_t chi2PerClusterTPC = -1;
 563   if (nClustersITS!=0)
 564     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
 565   if (nClustersTPC!=0)
 566     chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
 567
 568   Double_t extCov[15];
 569   esdTrack->GetExternalCovariance(extCov);
 570
 571   // getting the track to vertex parameters
 572   Float_t nSigmaToVertex = GetSigmaToVertex(esdTrack);
 573
 574   // getting the kinematic variables of the track
 575   // (assuming the mass is known)
 576   Double_t p[3];
 577   esdTrack->GetPxPyPz(p);
 578   Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
 579   Float_t pt       = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2));
 580   Float_t energy   = TMath::Sqrt(TMath::Power(esdTrack->GetMass(),2) + TMath::Power(momentum,2));
 581
 582
 583   //y-eta related calculations
 584   Float_t eta = -100.;
 585   Float_t y   = -100.;
 586   if((momentum != TMath::Abs(p[2]))&&(momentum != 0))
 587     eta = 0.5*TMath::Log((momentum + p[2])/(momentum - p[2]));
 588   if((energy != TMath::Abs(p[2]))&&(momentum != 0))
 589     y = 0.5*TMath::Log((energy + p[2])/(energy - p[2]));
 590
 591
 592   //########################################################################
 593   // cut the track?
 594
 595   Bool_t cuts[kNCuts];
 596   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
 597
 598   // track quality cuts
 599   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
 600     cuts[0]=kTRUE;
 601   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
 602     cuts[1]=kTRUE;
 603   if (nClustersTPC<fCutMinNClusterTPC)
 604     cuts[2]=kTRUE;
 605   if (nClustersITS<fCutMinNClusterITS)
 606     cuts[3]=kTRUE;
 607   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
 608     cuts[4]=kTRUE;
 609   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
 610     cuts[5]=kTRUE;
 611   if (extCov[0]  > fCutMaxC11)
 612     cuts[6]=kTRUE;
 613   if (extCov[2]  > fCutMaxC22)
 614     cuts[7]=kTRUE;
 615   if (extCov[5]  > fCutMaxC33)
 616     cuts[8]=kTRUE;
 617   if (extCov[9]  > fCutMaxC44)
 618     cuts[9]=kTRUE;
 619   if (extCov[14]  > fCutMaxC55)
 620     cuts[10]=kTRUE;
 621   if (nSigmaToVertex > fCutNsigmaToVertex)
 622     cuts[11] = kTRUE;
 623   // if n sigma could not be calculated
 624   if (nSigmaToVertex<0 && fCutSigmaToVertexRequired)
 625     cuts[12]=kTRUE;
 626   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
 627     cuts[13]=kTRUE;
 628   // track kinematics cut
 629   if((momentum < fPMin) || (momentum > fPMax))
 630     cuts[14]=kTRUE;
 631   if((pt < fPtMin) || (pt > fPtMax))
 632     cuts[15] = kTRUE;
 633   if((p[0] < fPxMin) || (p[0] > fPxMax))
 634     cuts[16] = kTRUE;
 635   if((p[1] < fPyMin) || (p[1] > fPyMax))
 636     cuts[17] = kTRUE;
 637   if((p[2] < fPzMin) || (p[2] > fPzMax))
 638     cuts[18] = kTRUE;
 639   if((eta < fEtaMin) || (eta > fEtaMax))
 640     cuts[19] = kTRUE;
 641   if((y < fRapMin) || (y > fRapMax))
 642     cuts[20] = kTRUE;
 643
 644   Bool_t cut=kFALSE;
 645   for (Int_t i=0; i<kNCuts; i++)
 646     if (cuts[i]) cut = kTRUE;
 647
 648   //########################################################################
 649   // filling histograms
 650   if (fHistogramsOn) {
 651     fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n tracks")));
 652
 653     if (cut)
 654       fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut tracks")));
 655
 656     for (Int_t i=0; i<kNCuts; i++) {
 657       if (cuts[i])
 658         fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
 659
 660       for (Int_t j=i; j<kNCuts; j++) {
 661         if (cuts[i] && cuts[j]) {
 662           Float_t x = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
 663           Float_t y = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
 664           fhCutCorrelation->Fill(x,y);
 665         }
 666       }
 667     }
 668
 669
 670     fhNClustersITS[0]->Fill(nClustersITS);
 671     fhNClustersTPC[0]->Fill(nClustersTPC);
 672     fhChi2PerClusterITS[0]->Fill(chi2PerClusterITS);
 673     fhChi2PerClusterTPC[0]->Fill(chi2PerClusterTPC);
 674
 675     fhC11[0]->Fill(extCov[0]);
 676     fhC22[0]->Fill(extCov[2]);
 677     fhC33[0]->Fill(extCov[5]);
 678     fhC44[0]->Fill(extCov[9]);
 679     fhC55[0]->Fill(extCov[14]);
 680
 681     Float_t b[2];
 682     Float_t bRes[2];
 683     Float_t bCov[3];
 684     esdTrack->GetImpactParameters(b,bCov);
 685     if (bCov[0]<=0 || bCov[2]<=0) {
 686       AliDebug(1, "Estimated b resolution lower or equal zero!");
 687       bCov[0]=0; bCov[2]=0;
 688     }
 689     bRes[0] = TMath::Sqrt(bCov[0]);
 690     bRes[1] = TMath::Sqrt(bCov[2]);
 691
 692     fhDZ[0]->Fill(b[1]);
 693     fhDXY[0]->Fill(b[0]);
 694     fhDXYvsDZ[0]->Fill(b[1],b[0]);
 695
 696     if (bRes[0]!=0 && bRes[1]!=0) {
 697       fhDZNormalized[0]->Fill(b[1]/bRes[1]);
 698       fhDXYNormalized[0]->Fill(b[0]/bRes[0]);
 699       fhDXYvsDZNormalized[0]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
 700       fhNSigmaToVertex[0]->Fill(nSigmaToVertex);
 701     }
 702   }
 703
 704   //########################################################################
 705   // cut the track!
 706   if (cut) return kFALSE;
 707
 708   //########################################################################
 709   // filling histograms after cut
 710   if (fHistogramsOn) {
 711     fhNClustersITS[1]->Fill(nClustersITS);
 712     fhNClustersTPC[1]->Fill(nClustersTPC);
 713     fhChi2PerClusterITS[1]->Fill(chi2PerClusterITS);
 714     fhChi2PerClusterTPC[1]->Fill(chi2PerClusterTPC);
 715
 716     fhC11[1]->Fill(extCov[0]);
 717     fhC22[1]->Fill(extCov[2]);
 718     fhC33[1]->Fill(extCov[5]);
 719     fhC44[1]->Fill(extCov[9]);
 720     fhC55[1]->Fill(extCov[14]);
 721
 722     Float_t b[2];
 723     Float_t bRes[2];
 724     Float_t bCov[3];
 725     esdTrack->GetImpactParameters(b,bCov);
 726     if (bCov[0]<=0 || bCov[2]<=0) {
 727       AliDebug(1, "Estimated b resolution lower or equal zero!");
 728       bCov[0]=0; bCov[2]=0;
 729     }
 730     bRes[0] = TMath::Sqrt(bCov[0]);
 731     bRes[1] = TMath::Sqrt(bCov[2]);
 732
 733     fhDZ[1]->Fill(b[1]);
 734     fhDXY[1]->Fill(b[0]);
 735     fhDXYvsDZ[1]->Fill(b[1],b[0]);
 736
 737     if (bRes[0]!=0 && bRes[1]!=0)
 738     {
 739       fhDZNormalized[1]->Fill(b[1]/bRes[1]);
 740       fhDXYNormalized[1]->Fill(b[0]/bRes[0]);
 741       fhDXYvsDZNormalized[1]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
 742       fhNSigmaToVertex[1]->Fill(nSigmaToVertex);
 743     }
 744   }
 745
 746   return kTRUE;
 747 }
 748
 749 //____________________________________________________________________
 750 TObjArray*
 751 AliESDtrackCuts::GetAcceptedTracks(AliESD* esd)
 752 {
 753   //
 754   // returns an array of all tracks that pass the cuts
 755   //
 756
 757   TObjArray* acceptedTracks = new TObjArray();
 758
 759   // loop over esd tracks
 760   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
 761     AliESDtrack* track = esd->GetTrack(iTrack);
 762
 763     if (AcceptTrack(track))
 764       acceptedTracks->Add(track);
 765   }
 766
 767   return acceptedTracks;
 768 }
 769
 770 //____________________________________________________________________
 771 Int_t
 772 AliESDtrackCuts::CountAcceptedTracks(AliESD* esd)
 773 {
 774   //
 775   // returns an the number of tracks that pass the cuts
 776   //
 777
 778   Int_t count = 0;
 779
 780   // loop over esd tracks
 781   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
 782     AliESDtrack* track = esd->GetTrack(iTrack);
 783
 784     if (AcceptTrack(track))
 785       count++;
 786   }
 787
 788   return count;
 789 }
 790
 791 //____________________________________________________________________
 792  void AliESDtrackCuts::DefineHistograms(Int_t color) {
 793    //
 794    // diagnostics histograms are defined
 795    //
 796
 797    fHistogramsOn=kTRUE;
 798
 799    //###################################################################################
 800    // defining histograms
 801
 802    fhCutStatistics = new TH1F("cut_statistics","cut statistics",kNCuts+4,-0.5,kNCuts+3.5);
 803
 804    fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
 805    fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
 806
 807    fhCutCorrelation = new TH2F("cut_correlation","cut correlation",kNCuts,-0.5,kNCuts-0.5,kNCuts,-0.5,kNCuts-0.5);;
 808
 809    for (Int_t i=0; i<kNCuts; i++) {
 810      fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fgkCutNames[i]);
 811      fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fgkCutNames[i]);
 812      fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fgkCutNames[i]);
 813    }
 814
 815   fhCutStatistics  ->SetLineColor(color);
 816   fhCutCorrelation ->SetLineColor(color);
 817   fhCutStatistics  ->SetLineWidth(2);
 818   fhCutCorrelation ->SetLineWidth(2);
 819
 820   Char_t str[256];
 821   for (Int_t i=0; i<2; i++) {
 822     if (i==0) sprintf(str," ");
 823     else sprintf(str,"_cut");
 824
 825     fhNClustersITS[i]        = new TH1F(Form("nClustersITS%s",str),"",8,-0.5,7.5);
 826     fhNClustersTPC[i]        = new TH1F(Form("nClustersTPC%s",str),"",165,-0.5,164.5);
 827     fhChi2PerClusterITS[i]   = new TH1F(Form("chi2PerClusterITS%s",str),"",500,0,10);
 828     fhChi2PerClusterTPC[i]   = new TH1F(Form("chi2PerClusterTPC%s",str),"",500,0,10);
 829
 830     fhC11[i]                 = new TH1F(Form("covMatrixDiagonal11%s",str),"",2000,0,20);
 831     fhC22[i]                 = new  TH1F(Form("covMatrixDiagonal22%s",str),"",2000,0,20);
 832     fhC33[i]                 = new  TH1F(Form("covMatrixDiagonal33%s",str),"",1000,0,1);
 833     fhC44[i]                 = new  TH1F(Form("covMatrixDiagonal44%s",str),"",1000,0,5);
 834     fhC55[i]                 = new  TH1F(Form("covMatrixDiagonal55%s",str),"",1000,0,5);
 835
 836     fhDXY[i]                 = new  TH1F(Form("dXY%s",str),"",500,-10,10);
 837     fhDZ[i]                  = new  TH1F(Form("dZ%s",str),"",500,-10,10);
 838     fhDXYvsDZ[i]             = new  TH2F(Form("dXYvsDZ%s",str),"",200,-10,10,200,-10,10);
 839
 840     fhDXYNormalized[i]       = new  TH1F(Form("dXYNormalized%s",str),"",500,-10,10);
 841     fhDZNormalized[i]        = new  TH1F(Form("dZNormalized%s",str),"",500,-10,10);
 842     fhDXYvsDZNormalized[i]   = new  TH2F(Form("dXYvsDZNormalized%s",str),"",200,-10,10,200,-10,10);
 843
 844     fhNSigmaToVertex[i]         = new  TH1F(Form("nSigmaToVertex%s",str),"",500,0,50);
 845
 846     fhNClustersITS[i]->SetXTitle("n ITS clusters");
 847     fhNClustersTPC[i]->SetXTitle("n TPC clusters");
 848     fhChi2PerClusterITS[i]->SetXTitle("#Chi^{2} per ITS cluster");
 849     fhChi2PerClusterTPC[i]->SetXTitle("#Chi^{2} per TPC cluster");
 850
 851     fhC11[i]->SetXTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
 852     fhC22[i]->SetXTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
 853     fhC33[i]->SetXTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
 854     fhC44[i]->SetXTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
 855     fhC55[i]->SetXTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
 856
 857     fhDXY[i]->SetXTitle("transverse impact parameter");
 858     fhDZ[i]->SetXTitle("longitudinal impact parameter");
 859     fhDXYvsDZ[i]->SetXTitle("longitudinal impact parameter");
 860     fhDXYvsDZ[i]->SetYTitle("transverse impact parameter");
 861
 862     fhDXYNormalized[i]->SetXTitle("normalized trans impact par");
 863     fhDZNormalized[i]->SetXTitle("normalized long impact par");
 864     fhDXYvsDZNormalized[i]->SetXTitle("normalized long impact par");
 865     fhDXYvsDZNormalized[i]->SetYTitle("normalized trans impact par");
 866     fhNSigmaToVertex[i]->SetXTitle("n #sigma to vertex");
 867
 868     fhNClustersITS[i]->SetLineColor(color);   fhNClustersITS[i]->SetLineWidth(2);
 869     fhNClustersTPC[i]->SetLineColor(color);   fhNClustersTPC[i]->SetLineWidth(2);
 870     fhChi2PerClusterITS[i]->SetLineColor(color);   fhChi2PerClusterITS[i]->SetLineWidth(2);
 871     fhChi2PerClusterTPC[i]->SetLineColor(color);   fhChi2PerClusterTPC[i]->SetLineWidth(2);
 872
 873     fhC11[i]->SetLineColor(color);   fhC11[i]->SetLineWidth(2);
 874     fhC22[i]->SetLineColor(color);   fhC22[i]->SetLineWidth(2);
 875     fhC33[i]->SetLineColor(color);   fhC33[i]->SetLineWidth(2);
 876     fhC44[i]->SetLineColor(color);   fhC44[i]->SetLineWidth(2);
 877     fhC55[i]->SetLineColor(color);   fhC55[i]->SetLineWidth(2);
 878
 879     fhDXY[i]->SetLineColor(color);   fhDXY[i]->SetLineWidth(2);
 880     fhDZ[i]->SetLineColor(color);   fhDZ[i]->SetLineWidth(2);
 881
 882     fhDXYNormalized[i]->SetLineColor(color);   fhDXYNormalized[i]->SetLineWidth(2);
 883     fhDZNormalized[i]->SetLineColor(color);    fhDZNormalized[i]->SetLineWidth(2);
 884     fhNSigmaToVertex[i]->SetLineColor(color);  fhNSigmaToVertex[i]->SetLineWidth(2);
 885   }
 886
 887   // The number of sigmas to the vertex is per definition gaussian
 888   ffDTheoretical = new TF1("nSigmaToVertexTheoretical","([0]/2.506628274)*exp(-(x**2)/2)",0,50);
 889   ffDTheoretical->SetParameter(0,1);
 890 }
 891
 892
 893
 894 //____________________________________________________________________
 895 void
 896 AliESDtrackCuts::Print(const Option_t*) const {
 897   //
 898   // print method - still to be implemented
 899   //
 900
 901   AliInfo("AliESDtrackCuts...");
 902 }
 903
 904
 905 //____________________________________________________________________
 906 void AliESDtrackCuts::SaveHistograms(Char_t* dir) {
 907   //
 908   // saves the histograms in a directory (dir)
 909   //
 910
 911
 912   if (!fHistogramsOn) {
 913     AliDebug(0, "Histograms not on - cannot save histograms!!!");
 914     return;
 915   }
 916
 917   gDirectory->mkdir(dir);
 918   gDirectory->cd(dir);
 919   TString outfile = 0;
 920   outfile = dir; outfile += ".root";
 921   TFile *f = TFile::Open(outfile.Data(),"recreate");
 922
 923   gDirectory->mkdir("before_cuts");
 924   gDirectory->mkdir("after_cuts");
 925
 926   // a factor of 2 is needed since n sigma is positive
 927   ffDTheoretical->SetParameter(0,2*fhNSigmaToVertex[0]->Integral("width"));
 928   ffDTheoretical->Write("nSigmaToVertexTheory");
 929
 930   fhCutStatistics->Write();
 931   fhCutCorrelation->Write();
 932
 933   for (Int_t i=0; i<2; i++) {
 934     if (i==0)
 935       gDirectory->cd("before_cuts");
 936     else
 937       gDirectory->cd("after_cuts");
 938     fhNClustersITS[i]        ->Write();
 939     fhNClustersTPC[i]        ->Write();
 940     fhChi2PerClusterITS[i]   ->Write();
 941     fhChi2PerClusterTPC[i]   ->Write();
 942
 943     fhC11[i]                 ->Write();
 944     fhC22[i]                 ->Write();
 945     fhC33[i]                 ->Write();
 946     fhC44[i]                 ->Write();
 947     fhC55[i]                 ->Write();
 948
 949     fhDXY[i]                 ->Write();
 950     fhDZ[i]                  ->Write();
 951     fhDXYvsDZ[i]             ->Write();
 952
 953     fhDXYNormalized[i]       ->Write();
 954     fhDZNormalized[i]        ->Write();
 955     fhDXYvsDZNormalized[i]   ->Write();
 956     fhNSigmaToVertex[i]      ->Write();
 957
 958     gDirectory->cd("../");
 959   }
 960
 961   f->Close();
 962
 963   gDirectory->cd("../");
 964 }
 965
 966
 967