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