PWG0/esdTrackCuts/AliESDtrackCuts.cxx

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