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