PWG0/esdTrackCuts/AliESDtrackCuts.cxx

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