ANALYSIS/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: AliESDtrackCuts.cxx 24534 2008-03-16 22:22:11Z fca $ */
  17
  18 #include "AliESDtrackCuts.h"
  19
  20 #include <AliESDtrack.h>
  21 #include <AliESDVertex.h>
  22 #include <AliESDEvent.h>
  23 #include <AliMultiplicity.h>
  24 #include <AliLog.h>
  25
  26 #include <TTree.h>
  27 #include <TCanvas.h>
  28 #include <TDirectory.h>
  29 #include <TH2F.h>
  30 #include <TF1.h>
  31 #include <TBits.h>
  32
  33 //____________________________________________________________________
  34 ClassImp(AliESDtrackCuts)
  35
  36 // Cut names
  37 const Char_t* AliESDtrackCuts::fgkCutNames[kNCuts] = {
  38  "require TPC refit",
  39  "require TPC standalone",
  40  "require ITS refit",
  41  "n clusters TPC",
  42  "n clusters ITS",
  43  "#Chi^{2}/cluster TPC",
  44  "#Chi^{2}/cluster ITS",
  45  "cov 11",
  46  "cov 22",
  47  "cov 33",
  48  "cov 44",
  49  "cov 55",
  50  "trk-to-vtx",
  51  "trk-to-vtx failed",
  52  "kink daughters",
  53  "p",
  54  "p_{T}",
  55  "p_{x}",
  56  "p_{y}",
  57  "p_{z}",
  58  "eta",
  59  "y",
  60  "trk-to-vtx max dca 2D absolute",
  61  "trk-to-vtx max dca xy absolute",
  62  "trk-to-vtx max dca z absolute",
  63  "trk-to-vtx min dca 2D absolute",
  64  "trk-to-vtx min dca xy absolute",
  65  "trk-to-vtx min dca z absolute",
  66  "SPD cluster requirement",
  67  "SDD cluster requirement",
  68  "SSD cluster requirement",
  69  "require ITS stand-alone",
  70  "rel 1/pt uncertainty",
  71  "TPC n shared clusters",
  72  "TPC rel shared clusters",
  73  "require ITS Pid",
  74  "n crossed rows TPC",
  75  "n crossed rows / n findable clusters",
  76  "missing ITS points",
  77  "#Chi^{2} TPC constrained vs. global",
  78  "require TOF out",
  79  "TOF Distance cut"
  80 };
  81
  82 AliESDtrackCuts* AliESDtrackCuts::fgMultEstTrackCuts[AliESDtrackCuts::kNMultEstTrackCuts] = { 0, 0, 0, 0 };
  83 Char_t AliESDtrackCuts::fgBeamTypeFlag = -1;
  84
  85 //____________________________________________________________________
  86 AliESDtrackCuts::AliESDtrackCuts(const Char_t* name, const Char_t* title) : AliAnalysisCuts(name,title),
  87   fCutMinNClusterTPC(0),
  88   fCutMinNClusterITS(0),
  89   fCutMinNCrossedRowsTPC(0),
  90   fCutMinRatioCrossedRowsOverFindableClustersTPC(0),
  91   f1CutMinNClustersTPCPtDep(0x0),
  92   fCutMaxPtDepNClustersTPC(0),
  93   fCutMaxChi2PerClusterTPC(0),
  94   fCutMaxChi2PerClusterITS(0),
  95   fCutMaxChi2TPCConstrainedVsGlobal(0),
  96   fCutMaxChi2TPCConstrainedVsGlobalVertexType(kVertexTracks | kVertexSPD),
  97   fCutMaxMissingITSPoints(0),
  98   fCutMaxC11(0),
  99   fCutMaxC22(0),
 100   fCutMaxC33(0),
 101   fCutMaxC44(0),
 102   fCutMaxC55(0),
 103   fCutMaxRel1PtUncertainty(0),
 104   fCutAcceptKinkDaughters(0),
 105   fCutAcceptSharedTPCClusters(0),
 106   fCutMaxFractionSharedTPCClusters(0),
 107   fCutRequireTPCRefit(0),
 108   fCutRequireTPCStandAlone(0),
 109   fCutRequireITSRefit(0),
 110   fCutRequireITSPid(0),
 111   fCutRequireITSStandAlone(0),
 112   fCutRequireITSpureSA(0),
 113   fCutNsigmaToVertex(0),
 114   fCutSigmaToVertexRequired(0),
 115   fCutMaxDCAToVertexXY(0),
 116   fCutMaxDCAToVertexZ(0),
 117   fCutMinDCAToVertexXY(0),
 118   fCutMinDCAToVertexZ(0),
 119   fCutMaxDCAToVertexXYPtDep(""),
 120   fCutMaxDCAToVertexZPtDep(""),
 121   fCutMinDCAToVertexXYPtDep(""),
 122   fCutMinDCAToVertexZPtDep(""),
 123   f1CutMaxDCAToVertexXYPtDep(0x0),
 124   f1CutMaxDCAToVertexZPtDep(0x0),
 125   f1CutMinDCAToVertexXYPtDep(0x0),
 126   f1CutMinDCAToVertexZPtDep(0x0),
 127   fCutDCAToVertex2D(0),
 128   fPMin(0),
 129   fPMax(0),
 130   fPtMin(0),
 131   fPtMax(0),
 132   fPxMin(0),
 133   fPxMax(0),
 134   fPyMin(0),
 135   fPyMax(0),
 136   fPzMin(0),
 137   fPzMax(0),
 138   fEtaMin(0),
 139   fEtaMax(0),
 140   fRapMin(0),
 141   fRapMax(0),
 142   fCutRequireTOFout(kFALSE),
 143   fFlagCutTOFdistance(kFALSE),
 144   fCutTOFdistance(3.),
 145   fHistogramsOn(0),
 146   ffDTheoretical(0),
 147   fhCutStatistics(0),
 148   fhCutCorrelation(0)
 149 {
 150   //
 151   // constructor
 152   //
 153
 154   Init();
 155
 156   //##############################################################################
 157   // setting default cuts
 158   SetMinNClustersTPC();
 159   SetMinNClustersITS();
 160   SetMinNCrossedRowsTPC();
 161   SetMinRatioCrossedRowsOverFindableClustersTPC();
 162   SetMaxChi2PerClusterTPC();
 163   SetMaxChi2PerClusterITS();
 164   SetMaxChi2TPCConstrainedGlobal();
 165   SetMaxChi2TPCConstrainedGlobalVertexType();
 166   SetMaxNOfMissingITSPoints();
 167   SetMaxCovDiagonalElements();
 168   SetMaxRel1PtUncertainty();
 169   SetRequireTPCRefit();
 170   SetRequireTPCStandAlone();
 171   SetRequireITSRefit();
 172   SetRequireITSPid(kFALSE);
 173   SetRequireITSStandAlone(kFALSE);
 174   SetRequireITSPureStandAlone(kFALSE);
 175   SetAcceptKinkDaughters();
 176   SetAcceptSharedTPCClusters();
 177   SetMaxFractionSharedTPCClusters();
 178   SetMaxNsigmaToVertex();
 179   SetMaxDCAToVertexXY();
 180   SetMaxDCAToVertexZ();
 181   SetDCAToVertex2D();
 182   SetMinDCAToVertexXY();
 183   SetMinDCAToVertexZ();
 184   SetPRange();
 185   SetPtRange();
 186   SetPxRange();
 187   SetPyRange();
 188   SetPzRange();
 189   SetEtaRange();
 190   SetRapRange();
 191   SetClusterRequirementITS(kSPD);
 192   SetClusterRequirementITS(kSDD);
 193   SetClusterRequirementITS(kSSD);
 194
 195   SetHistogramsOn();
 196 }
 197
 198 //_____________________________________________________________________________
 199 AliESDtrackCuts::AliESDtrackCuts(const AliESDtrackCuts &c) : AliAnalysisCuts(c),
 200   fCutMinNClusterTPC(0),
 201   fCutMinNClusterITS(0),
 202   fCutMinNCrossedRowsTPC(0),
 203   fCutMinRatioCrossedRowsOverFindableClustersTPC(0),
 204   f1CutMinNClustersTPCPtDep(0x0),
 205   fCutMaxPtDepNClustersTPC(0),
 206   fCutMaxChi2PerClusterTPC(0),
 207   fCutMaxChi2PerClusterITS(0),
 208   fCutMaxChi2TPCConstrainedVsGlobal(0),
 209   fCutMaxChi2TPCConstrainedVsGlobalVertexType(kVertexTracks | kVertexSPD),
 210   fCutMaxMissingITSPoints(0),
 211   fCutMaxC11(0),
 212   fCutMaxC22(0),
 213   fCutMaxC33(0),
 214   fCutMaxC44(0),
 215   fCutMaxC55(0),
 216   fCutMaxRel1PtUncertainty(0),
 217   fCutAcceptKinkDaughters(0),
 218   fCutAcceptSharedTPCClusters(0),
 219   fCutMaxFractionSharedTPCClusters(0),
 220   fCutRequireTPCRefit(0),
 221   fCutRequireTPCStandAlone(0),
 222   fCutRequireITSRefit(0),
 223   fCutRequireITSPid(0),
 224   fCutRequireITSStandAlone(0),
 225   fCutRequireITSpureSA(0),
 226   fCutNsigmaToVertex(0),
 227   fCutSigmaToVertexRequired(0),
 228   fCutMaxDCAToVertexXY(0),
 229   fCutMaxDCAToVertexZ(0),
 230   fCutMinDCAToVertexXY(0),
 231   fCutMinDCAToVertexZ(0),
 232   fCutMaxDCAToVertexXYPtDep(""),
 233   fCutMaxDCAToVertexZPtDep(""),
 234   fCutMinDCAToVertexXYPtDep(""),
 235   fCutMinDCAToVertexZPtDep(""),
 236   f1CutMaxDCAToVertexXYPtDep(0x0),
 237   f1CutMaxDCAToVertexZPtDep(0x0),
 238   f1CutMinDCAToVertexXYPtDep(0x0),
 239   f1CutMinDCAToVertexZPtDep(0x0),
 240   fCutDCAToVertex2D(0),
 241   fPMin(0),
 242   fPMax(0),
 243   fPtMin(0),
 244   fPtMax(0),
 245   fPxMin(0),
 246   fPxMax(0),
 247   fPyMin(0),
 248   fPyMax(0),
 249   fPzMin(0),
 250   fPzMax(0),
 251   fEtaMin(0),
 252   fEtaMax(0),
 253   fRapMin(0),
 254   fRapMax(0),
 255   fCutRequireTOFout(kFALSE),
 256   fFlagCutTOFdistance(kFALSE),
 257   fCutTOFdistance(3.),
 258   fHistogramsOn(0),
 259   ffDTheoretical(0),
 260   fhCutStatistics(0),
 261   fhCutCorrelation(0)
 262 {
 263   //
 264   // copy constructor
 265   //
 266
 267   ((AliESDtrackCuts &) c).Copy(*this);
 268 }
 269
 270 AliESDtrackCuts::~AliESDtrackCuts()
 271 {
 272   //
 273   // destructor
 274   //
 275
 276   for (Int_t i=0; i<2; i++) {
 277
 278     if (fhNClustersITS[i])
 279       delete fhNClustersITS[i];
 280     if (fhNClustersTPC[i])
 281       delete fhNClustersTPC[i];
 282     if (fhNSharedClustersTPC[i])
 283       delete fhNSharedClustersTPC[i];
 284     if (fhNCrossedRowsTPC[i])
 285       delete fhNCrossedRowsTPC[i];
 286     if (fhRatioCrossedRowsOverFindableClustersTPC[i])
 287       delete fhRatioCrossedRowsOverFindableClustersTPC[i];
 288     if (fhChi2PerClusterITS[i])
 289       delete fhChi2PerClusterITS[i];
 290     if (fhChi2PerClusterTPC[i])
 291       delete fhChi2PerClusterTPC[i];
 292     if (fhChi2TPCConstrainedVsGlobal[i])
 293       delete fhChi2TPCConstrainedVsGlobal[i];
 294     if(fhNClustersForITSPID[i])
 295       delete fhNClustersForITSPID[i];
 296     if(fhNMissingITSPoints[i])
 297       delete fhNMissingITSPoints[i];
 298     if (fhC11[i])
 299       delete fhC11[i];
 300     if (fhC22[i])
 301       delete fhC22[i];
 302     if (fhC33[i])
 303       delete fhC33[i];
 304     if (fhC44[i])
 305       delete fhC44[i];
 306     if (fhC55[i])
 307       delete fhC55[i];
 308
 309     if (fhRel1PtUncertainty[i])
 310       delete fhRel1PtUncertainty[i];
 311
 312     if (fhDXY[i])
 313       delete fhDXY[i];
 314     if (fhDZ[i])
 315       delete fhDZ[i];
 316     if (fhDXYDZ[i])
 317       delete fhDXYDZ[i];
 318     if (fhDXYvsDZ[i])
 319       delete fhDXYvsDZ[i];
 320
 321     if (fhDXYNormalized[i])
 322       delete fhDXYNormalized[i];
 323     if (fhDZNormalized[i])
 324       delete fhDZNormalized[i];
 325     if (fhDXYvsDZNormalized[i])
 326       delete fhDXYvsDZNormalized[i];
 327     if (fhNSigmaToVertex[i])
 328       delete fhNSigmaToVertex[i];
 329     if (fhPt[i])
 330       delete fhPt[i];
 331     if (fhEta[i])
 332       delete fhEta[i];
 333     if (fhTOFdistance[i])
 334       delete fhTOFdistance[i];
 335   }
 336
 337   if(f1CutMaxDCAToVertexXYPtDep)delete f1CutMaxDCAToVertexXYPtDep;
 338   f1CutMaxDCAToVertexXYPtDep = 0;
 339   if( f1CutMaxDCAToVertexZPtDep) delete  f1CutMaxDCAToVertexZPtDep;
 340   f1CutMaxDCAToVertexZPtDep = 0;
 341   if( f1CutMinDCAToVertexXYPtDep)delete  f1CutMinDCAToVertexXYPtDep;
 342   f1CutMinDCAToVertexXYPtDep = 0;
 343   if(f1CutMinDCAToVertexZPtDep)delete  f1CutMinDCAToVertexZPtDep;
 344   f1CutMinDCAToVertexZPtDep = 0;
 345
 346
 347   if (ffDTheoretical)
 348     delete ffDTheoretical;
 349
 350   if (fhCutStatistics)
 351     delete fhCutStatistics;
 352   if (fhCutCorrelation)
 353     delete fhCutCorrelation;
 354
 355   if(f1CutMinNClustersTPCPtDep)
 356     delete f1CutMinNClustersTPCPtDep;
 357
 358 }
 359
 360 void AliESDtrackCuts::Init()
 361 {
 362   //
 363   // sets everything to zero
 364   //
 365
 366   fCutMinNClusterTPC = 0;
 367   fCutMinNClusterITS = 0;
 368
 369   fCutMaxChi2PerClusterTPC = 0;
 370   fCutMaxChi2PerClusterITS = 0;
 371   fCutMaxChi2TPCConstrainedVsGlobal = 0;
 372   fCutMaxChi2TPCConstrainedVsGlobalVertexType = kVertexTracks | kVertexSPD;
 373   fCutMaxMissingITSPoints  = 0;
 374
 375   for (Int_t i = 0; i < 3; i++)
 376         fCutClusterRequirementITS[i] = kOff;
 377
 378   fCutMaxC11 = 0;
 379   fCutMaxC22 = 0;
 380   fCutMaxC33 = 0;
 381   fCutMaxC44 = 0;
 382   fCutMaxC55 = 0;
 383
 384   fCutMaxRel1PtUncertainty = 0;
 385
 386   fCutAcceptKinkDaughters = 0;
 387   fCutAcceptSharedTPCClusters = 0;
 388   fCutMaxFractionSharedTPCClusters = 0;
 389   fCutRequireTPCRefit = 0;
 390   fCutRequireTPCStandAlone = 0;
 391   fCutRequireITSRefit = 0;
 392   fCutRequireITSPid = 0;
 393   fCutRequireITSStandAlone = 0;
 394   fCutRequireITSpureSA = 0;
 395
 396   fCutNsigmaToVertex = 0;
 397   fCutSigmaToVertexRequired = 0;
 398   fCutMaxDCAToVertexXY = 0;
 399   fCutMaxDCAToVertexZ = 0;
 400   fCutDCAToVertex2D = 0;
 401   fCutMinDCAToVertexXY = 0;
 402   fCutMinDCAToVertexZ = 0;
 403   fCutMaxDCAToVertexXYPtDep = "";
 404   fCutMaxDCAToVertexZPtDep = "";
 405   fCutMinDCAToVertexXYPtDep = "";
 406   fCutMinDCAToVertexZPtDep = "";
 407
 408   if(f1CutMaxDCAToVertexXYPtDep)delete f1CutMaxDCAToVertexXYPtDep;
 409   f1CutMaxDCAToVertexXYPtDep = 0;
 410   if( f1CutMaxDCAToVertexXYPtDep) delete  f1CutMaxDCAToVertexXYPtDep;
 411   f1CutMaxDCAToVertexXYPtDep = 0;
 412   if( f1CutMaxDCAToVertexZPtDep) delete  f1CutMaxDCAToVertexZPtDep;
 413   f1CutMaxDCAToVertexZPtDep = 0;
 414   if( f1CutMinDCAToVertexXYPtDep)delete  f1CutMinDCAToVertexXYPtDep;
 415   f1CutMinDCAToVertexXYPtDep = 0;
 416   if(f1CutMinDCAToVertexZPtDep)delete f1CutMinDCAToVertexZPtDep;
 417   f1CutMinDCAToVertexZPtDep = 0;
 418
 419
 420   fPMin = 0;
 421   fPMax = 0;
 422   fPtMin = 0;
 423   fPtMax = 0;
 424   fPxMin = 0;
 425   fPxMax = 0;
 426   fPyMin = 0;
 427   fPyMax = 0;
 428   fPzMin = 0;
 429   fPzMax = 0;
 430   fEtaMin = 0;
 431   fEtaMax = 0;
 432   fRapMin = 0;
 433   fRapMax = 0;
 434
 435   fHistogramsOn = kFALSE;
 436
 437   for (Int_t i=0; i<2; ++i)
 438   {
 439     fhNClustersITS[i] = 0;
 440     fhNClustersTPC[i] = 0;
 441     fhNSharedClustersTPC[i] = 0;
 442     fhNCrossedRowsTPC[i] = 0;
 443     fhRatioCrossedRowsOverFindableClustersTPC[i] = 0;
 444
 445     fhChi2PerClusterITS[i] = 0;
 446     fhChi2PerClusterTPC[i] = 0;
 447     fhChi2TPCConstrainedVsGlobal[i] = 0;
 448     fhNClustersForITSPID[i] = 0;
 449     fhNMissingITSPoints[i] = 0;
 450
 451     fhC11[i] = 0;
 452     fhC22[i] = 0;
 453     fhC33[i] = 0;
 454     fhC44[i] = 0;
 455     fhC55[i] = 0;
 456
 457     fhRel1PtUncertainty[i] = 0;
 458
 459     fhDXY[i] = 0;
 460     fhDZ[i] = 0;
 461     fhDXYDZ[i] = 0;
 462     fhDXYvsDZ[i] = 0;
 463
 464     fhDXYNormalized[i] = 0;
 465     fhDZNormalized[i] = 0;
 466     fhDXYvsDZNormalized[i] = 0;
 467     fhNSigmaToVertex[i] = 0;
 468
 469     fhPt[i] = 0;
 470     fhEta[i] = 0;
 471     fhTOFdistance[i] = 0;
 472   }
 473   ffDTheoretical = 0;
 474
 475   fhCutStatistics = 0;
 476   fhCutCorrelation = 0;
 477 }
 478
 479 //_____________________________________________________________________________
 480 AliESDtrackCuts &AliESDtrackCuts::operator=(const AliESDtrackCuts &c)
 481 {
 482   //
 483   // Assignment operator
 484   //
 485
 486   if (this != &c) ((AliESDtrackCuts &) c).Copy(*this);
 487   return *this;
 488 }
 489
 490 //_____________________________________________________________________________
 491 void AliESDtrackCuts::Copy(TObject &c) const
 492 {
 493   //
 494   // Copy function
 495   //
 496
 497   AliESDtrackCuts& target = (AliESDtrackCuts &) c;
 498
 499   target.Init();
 500
 501   target.fCutMinNClusterTPC = fCutMinNClusterTPC;
 502   target.fCutMinNClusterITS = fCutMinNClusterITS;
 503   target.fCutMinNCrossedRowsTPC = fCutMinNCrossedRowsTPC;
 504   target.fCutMinRatioCrossedRowsOverFindableClustersTPC = fCutMinRatioCrossedRowsOverFindableClustersTPC;
 505   if(f1CutMinNClustersTPCPtDep){
 506     target.f1CutMinNClustersTPCPtDep = (TFormula*) f1CutMinNClustersTPCPtDep->Clone("f1CutMinNClustersTPCPtDep");
 507   }
 508   target.fCutMaxPtDepNClustersTPC =   fCutMaxPtDepNClustersTPC;
 509
 510   target.fCutMaxChi2PerClusterTPC = fCutMaxChi2PerClusterTPC;
 511   target.fCutMaxChi2PerClusterITS = fCutMaxChi2PerClusterITS;
 512   target.fCutMaxChi2TPCConstrainedVsGlobal = fCutMaxChi2TPCConstrainedVsGlobal;
 513   target.fCutMaxChi2TPCConstrainedVsGlobalVertexType = fCutMaxChi2TPCConstrainedVsGlobalVertexType;
 514   target.fCutMaxMissingITSPoints = fCutMaxMissingITSPoints;
 515
 516   for (Int_t i = 0; i < 3; i++)
 517     target.fCutClusterRequirementITS[i] = fCutClusterRequirementITS[i];
 518
 519   target.fCutMaxC11 = fCutMaxC11;
 520   target.fCutMaxC22 = fCutMaxC22;
 521   target.fCutMaxC33 = fCutMaxC33;
 522   target.fCutMaxC44 = fCutMaxC44;
 523   target.fCutMaxC55 = fCutMaxC55;
 524
 525   target.fCutMaxRel1PtUncertainty = fCutMaxRel1PtUncertainty;
 526
 527   target.fCutAcceptKinkDaughters = fCutAcceptKinkDaughters;
 528   target.fCutAcceptSharedTPCClusters = fCutAcceptSharedTPCClusters;
 529   target.fCutMaxFractionSharedTPCClusters = fCutMaxFractionSharedTPCClusters;
 530   target.fCutRequireTPCRefit = fCutRequireTPCRefit;
 531   target.fCutRequireTPCStandAlone = fCutRequireTPCStandAlone;
 532   target.fCutRequireITSRefit = fCutRequireITSRefit;
 533   target.fCutRequireITSPid = fCutRequireITSPid;
 534   target.fCutRequireITSStandAlone = fCutRequireITSStandAlone;
 535   target.fCutRequireITSpureSA = fCutRequireITSpureSA;
 536
 537   target.fCutNsigmaToVertex = fCutNsigmaToVertex;
 538   target.fCutSigmaToVertexRequired = fCutSigmaToVertexRequired;
 539   target.fCutMaxDCAToVertexXY = fCutMaxDCAToVertexXY;
 540   target.fCutMaxDCAToVertexZ = fCutMaxDCAToVertexZ;
 541   target.fCutDCAToVertex2D = fCutDCAToVertex2D;
 542   target.fCutMinDCAToVertexXY = fCutMinDCAToVertexXY;
 543   target.fCutMinDCAToVertexZ = fCutMinDCAToVertexZ;
 544
 545   target.fCutMaxDCAToVertexXYPtDep = fCutMaxDCAToVertexXYPtDep;
 546   if(fCutMaxDCAToVertexXYPtDep.Length()>0)target.SetMaxDCAToVertexXYPtDep(fCutMaxDCAToVertexXYPtDep.Data());
 547
 548   target.fCutMaxDCAToVertexZPtDep = fCutMaxDCAToVertexZPtDep;
 549   if(fCutMaxDCAToVertexZPtDep.Length()>0)target.SetMaxDCAToVertexZPtDep(fCutMaxDCAToVertexZPtDep.Data());
 550
 551   target.fCutMinDCAToVertexXYPtDep = fCutMinDCAToVertexXYPtDep;
 552   if(fCutMinDCAToVertexXYPtDep.Length()>0)target.SetMinDCAToVertexXYPtDep(fCutMinDCAToVertexXYPtDep.Data());
 553
 554   target.fCutMinDCAToVertexZPtDep = fCutMinDCAToVertexZPtDep;
 555   if(fCutMinDCAToVertexZPtDep.Length()>0)target.SetMinDCAToVertexZPtDep(fCutMinDCAToVertexZPtDep.Data());
 556
 557   target.fPMin = fPMin;
 558   target.fPMax = fPMax;
 559   target.fPtMin = fPtMin;
 560   target.fPtMax = fPtMax;
 561   target.fPxMin = fPxMin;
 562   target.fPxMax = fPxMax;
 563   target.fPyMin = fPyMin;
 564   target.fPyMax = fPyMax;
 565   target.fPzMin = fPzMin;
 566   target.fPzMax = fPzMax;
 567   target.fEtaMin = fEtaMin;
 568   target.fEtaMax = fEtaMax;
 569   target.fRapMin = fRapMin;
 570   target.fRapMax = fRapMax;
 571
 572   target.fFlagCutTOFdistance = fFlagCutTOFdistance;
 573   target.fCutTOFdistance = fCutTOFdistance;
 574   target.fCutRequireTOFout = fCutRequireTOFout;
 575
 576   target.fHistogramsOn = fHistogramsOn;
 577
 578   for (Int_t i=0; i<2; ++i)
 579   {
 580     if (fhNClustersITS[i]) target.fhNClustersITS[i] = (TH1F*) fhNClustersITS[i]->Clone();
 581     if (fhNClustersTPC[i]) target.fhNClustersTPC[i] = (TH1F*) fhNClustersTPC[i]->Clone();
 582     if (fhNSharedClustersTPC[i]) target.fhNSharedClustersTPC[i] = (TH1F*) fhNSharedClustersTPC[i]->Clone();
 583     if (fhNCrossedRowsTPC[i]) target.fhNCrossedRowsTPC[i] = (TH1F*) fhNCrossedRowsTPC[i]->Clone();
 584     if (fhRatioCrossedRowsOverFindableClustersTPC[i]) target.fhRatioCrossedRowsOverFindableClustersTPC[i] = (TH1F*) fhRatioCrossedRowsOverFindableClustersTPC[i]->Clone();
 585
 586     if (fhChi2PerClusterITS[i]) target.fhChi2PerClusterITS[i] = (TH1F*) fhChi2PerClusterITS[i]->Clone();
 587     if (fhChi2PerClusterTPC[i]) target.fhChi2PerClusterTPC[i] = (TH1F*) fhChi2PerClusterTPC[i]->Clone();
 588     if (fhChi2TPCConstrainedVsGlobal[i]) target.fhChi2TPCConstrainedVsGlobal[i] = (TH1F*) fhChi2TPCConstrainedVsGlobal[i]->Clone();
 589     if (fhNClustersForITSPID[i]) target.fhNClustersForITSPID[i] = (TH1F*) fhNClustersForITSPID[i]->Clone();
 590     if (fhNMissingITSPoints[i]) target.fhNMissingITSPoints[i] = (TH1F*) fhNMissingITSPoints[i]->Clone();
 591
 592     if (fhC11[i]) target.fhC11[i] = (TH1F*) fhC11[i]->Clone();
 593     if (fhC22[i]) target.fhC22[i] = (TH1F*) fhC22[i]->Clone();
 594     if (fhC33[i]) target.fhC33[i] = (TH1F*) fhC33[i]->Clone();
 595     if (fhC44[i]) target.fhC44[i] = (TH1F*) fhC44[i]->Clone();
 596     if (fhC55[i]) target.fhC55[i] = (TH1F*) fhC55[i]->Clone();
 597
 598     if (fhRel1PtUncertainty[i]) target.fhRel1PtUncertainty[i] = (TH1F*) fhRel1PtUncertainty[i]->Clone();
 599
 600     if (fhDXY[i]) target.fhDXY[i] = (TH1F*) fhDXY[i]->Clone();
 601     if (fhDZ[i]) target.fhDZ[i] = (TH1F*) fhDZ[i]->Clone();
 602     if (fhDXYDZ[i]) target.fhDXYDZ[i] = (TH1F*) fhDXYDZ[i]->Clone();
 603     if (fhDXYvsDZ[i]) target.fhDXYvsDZ[i] = (TH2F*) fhDXYvsDZ[i]->Clone();
 604
 605     if (fhDXYNormalized[i]) target.fhDXYNormalized[i] = (TH1F*) fhDXYNormalized[i]->Clone();
 606     if (fhDZNormalized[i]) target.fhDZNormalized[i] = (TH1F*) fhDZNormalized[i]->Clone();
 607     if (fhDXYvsDZNormalized[i]) target.fhDXYvsDZNormalized[i] = (TH2F*) fhDXYvsDZNormalized[i]->Clone();
 608     if (fhNSigmaToVertex[i]) target.fhNSigmaToVertex[i] = (TH1F*) fhNSigmaToVertex[i]->Clone();
 609
 610     if (fhPt[i]) target.fhPt[i] = (TH1F*) fhPt[i]->Clone();
 611     if (fhEta[i]) target.fhEta[i] = (TH1F*) fhEta[i]->Clone();
 612     if (fhTOFdistance[i]) target.fhTOFdistance[i] = (TH2F*) fhTOFdistance[i]->Clone();
 613   }
 614   if (ffDTheoretical) target.ffDTheoretical = (TF1*) ffDTheoretical->Clone();
 615
 616   if (fhCutStatistics) target.fhCutStatistics = (TH1F*) fhCutStatistics->Clone();
 617   if (fhCutCorrelation) target.fhCutCorrelation = (TH2F*) fhCutCorrelation->Clone();
 618
 619   TNamed::Copy(c);
 620 }
 621
 622 //_____________________________________________________________________________
 623 Long64_t AliESDtrackCuts::Merge(TCollection* list) {
 624   // Merge a list of AliESDtrackCuts objects with this (needed for PROOF)
 625   // Returns the number of merged objects (including this)
 626   if (!list)
 627     return 0;
 628   if (list->IsEmpty())
 629     return 1;
 630   if (!fHistogramsOn)
 631     return 0;
 632   TIterator* iter = list->MakeIterator();
 633   TObject* obj;
 634
 635   // collection of measured and generated histograms
 636   Int_t count = 0;
 637   while ((obj = iter->Next())) {
 638
 639     AliESDtrackCuts* entry = dynamic_cast<AliESDtrackCuts*>(obj);
 640     if (entry == 0)
 641       continue;
 642
 643     if (!entry->fHistogramsOn)
 644       continue;
 645
 646     for (Int_t i=0; i<2; i++) {
 647
 648       fhNClustersITS[i]      ->Add(entry->fhNClustersITS[i]     );
 649       fhNClustersTPC[i]      ->Add(entry->fhNClustersTPC[i]     );
 650       if (fhNSharedClustersTPC[i])
 651         fhNSharedClustersTPC[i]      ->Add(entry->fhNSharedClustersTPC[i]     );
 652       if (fhNCrossedRowsTPC[i])
 653         fhNCrossedRowsTPC[i]   ->Add(entry->fhNCrossedRowsTPC[i]     );
 654       if (fhRatioCrossedRowsOverFindableClustersTPC[i])
 655         fhRatioCrossedRowsOverFindableClustersTPC[i]      ->Add(entry->fhRatioCrossedRowsOverFindableClustersTPC[i]     );
 656
 657       fhChi2PerClusterITS[i] ->Add(entry->fhChi2PerClusterITS[i]);
 658       fhChi2PerClusterTPC[i] ->Add(entry->fhChi2PerClusterTPC[i]);
 659       if (fhChi2TPCConstrainedVsGlobal[i])
 660         fhChi2TPCConstrainedVsGlobal[i]->Add(entry->fhChi2TPCConstrainedVsGlobal[i]);
 661       if (fhNClustersForITSPID[i])
 662         fhNClustersForITSPID[i]->Add(entry->fhNClustersForITSPID[i]);
 663       if (fhNMissingITSPoints[i])
 664         fhNMissingITSPoints[i] ->Add(entry->fhNMissingITSPoints[i]);
 665
 666       fhC11[i]               ->Add(entry->fhC11[i]              );
 667       fhC22[i]               ->Add(entry->fhC22[i]              );
 668       fhC33[i]               ->Add(entry->fhC33[i]              );
 669       fhC44[i]               ->Add(entry->fhC44[i]              );
 670       fhC55[i]               ->Add(entry->fhC55[i]              );
 671
 672       fhRel1PtUncertainty[i] ->Add(entry->fhRel1PtUncertainty[i]);
 673
 674       fhDXY[i]               ->Add(entry->fhDXY[i]              );
 675       fhDZ[i]                ->Add(entry->fhDZ[i]               );
 676       fhDXYDZ[i]             ->Add(entry->fhDXYDZ[i]          );
 677       fhDXYvsDZ[i]           ->Add(entry->fhDXYvsDZ[i]          );
 678
 679       fhDXYNormalized[i]     ->Add(entry->fhDXYNormalized[i]    );
 680       fhDZNormalized[i]      ->Add(entry->fhDZNormalized[i]     );
 681       fhDXYvsDZNormalized[i] ->Add(entry->fhDXYvsDZNormalized[i]);
 682       fhNSigmaToVertex[i]    ->Add(entry->fhNSigmaToVertex[i]);
 683
 684       fhPt[i]                ->Add(entry->fhPt[i]);
 685       fhEta[i]               ->Add(entry->fhEta[i]);
 686       fhTOFdistance[i]       ->Add(entry->fhTOFdistance[i]);
 687     }
 688
 689     fhCutStatistics  ->Add(entry->fhCutStatistics);
 690     fhCutCorrelation ->Add(entry->fhCutCorrelation);
 691
 692     count++;
 693   }
 694   return count+1;
 695 }
 696
 697 void AliESDtrackCuts::SetMinNClustersTPCPtDep(TFormula *f1, Float_t ptmax)
 698 {
 699   //
 700   // Sets the pT dependent NClustersTPC cut
 701   //
 702
 703   if(f1){
 704     delete f1CutMinNClustersTPCPtDep;
 705     f1CutMinNClustersTPCPtDep = (TFormula*)f1->Clone("f1CutMinNClustersTPCPtDep");
 706   }
 707   fCutMaxPtDepNClustersTPC=ptmax;
 708 }
 709
 710 //____________________________________________________________________
 711 AliESDtrackCuts* AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
 712 {
 713   // creates an AliESDtrackCuts object and fills it with standard (pre data-taking) values for TPC-only cuts
 714
 715   AliInfoClass("Creating track cuts for TPC-only.");
 716
 717   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 718
 719   esdTrackCuts->SetMinNClustersTPC(50);
 720   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
 721   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 722
 723   esdTrackCuts->SetMaxDCAToVertexZ(3.2);
 724   esdTrackCuts->SetMaxDCAToVertexXY(2.4);
 725   esdTrackCuts->SetDCAToVertex2D(kTRUE);
 726
 727   return esdTrackCuts;
 728 }
 729
 730 //____________________________________________________________________
 731 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(Bool_t selPrimaries)
 732 {
 733   // creates an AliESDtrackCuts object and fills it with standard values for ITS-TPC cuts for pp 2009 data
 734
 735   AliInfoClass("Creating track cuts for ITS+TPC (2009 definition).");
 736
 737   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 738
 739   // TPC
 740   esdTrackCuts->SetRequireTPCStandAlone(kTRUE); // to get chi2 and ncls of kTPCin
 741   esdTrackCuts->SetMinNClustersTPC(70);
 742   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
 743   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 744   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 745   // ITS
 746   esdTrackCuts->SetRequireITSRefit(kTRUE);
 747   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 748                                          AliESDtrackCuts::kAny);
 749   if(selPrimaries) {
 750     // 7*(0.0050+0.0060/pt^0.9)
 751     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0350+0.0420/pt^0.9");
 752     esdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
 753   }
 754   esdTrackCuts->SetMaxDCAToVertexZ(1.e6);
 755   esdTrackCuts->SetDCAToVertex2D(kFALSE);
 756   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 757   //esdTrackCuts->SetEtaRange(-0.8,+0.8);
 758
 759   esdTrackCuts->SetMaxChi2PerClusterITS(36);
 760
 761   return esdTrackCuts;
 762 }
 763
 764 //____________________________________________________________________
 765 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(Bool_t selPrimaries, Int_t clusterCut)
 766 {
 767   // creates an AliESDtrackCuts object and fills it with standard values for ITS-TPC cuts for pp 2011 data
 768   // if clusterCut = 1, the cut on the number of clusters is replaced by
 769   // a cut on the number of crossed rows and on the ration crossed
 770   // rows/findable clusters
 771
 772   AliInfoClass("Creating track cuts for ITS+TPC (2011 definition).");
 773
 774   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 775
 776   // TPC
 777   if(clusterCut == 0)  esdTrackCuts->SetMinNClustersTPC(50);
 778   else if (clusterCut == 1) {
 779     esdTrackCuts->SetMinNCrossedRowsTPC(70);
 780     esdTrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
 781   }
 782   else {
 783     AliWarningClass(Form("Wrong value of the clusterCut parameter (%d), using cut on Nclusters",clusterCut));
 784     esdTrackCuts->SetMinNClustersTPC(50);
 785   }
 786   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
 787   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 788   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 789   // ITS
 790   esdTrackCuts->SetRequireITSRefit(kTRUE);
 791   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 792                                          AliESDtrackCuts::kAny);
 793   if(selPrimaries) {
 794     // 7*(0.0015+0.0050/pt^1.1)
 795     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0105+0.0350/pt^1.1");
 796     esdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
 797   }
 798   esdTrackCuts->SetMaxDCAToVertexZ(2);
 799   esdTrackCuts->SetDCAToVertex2D(kFALSE);
 800   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 801
 802   esdTrackCuts->SetMaxChi2PerClusterITS(36);
 803
 804   return esdTrackCuts;
 805 }
 806
 807 //____________________________________________________________________
 808 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(Bool_t selPrimaries,Int_t clusterCut)
 809 {
 810   // creates an AliESDtrackCuts object and fills it with standard values for ITS-TPC cuts for pp 2010 data
 811   // if clusterCut = 1, the cut on the number of clusters is replaced by
 812   // a cut on the number of crossed rows and on the ration crossed
 813   // rows/findable clusters
 814
 815   AliInfoClass("Creating track cuts for ITS+TPC (2010 definition).");
 816
 817   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 818
 819   // TPC
 820   if(clusterCut == 0)  esdTrackCuts->SetMinNClustersTPC(70);
 821   else if (clusterCut == 1) {
 822     esdTrackCuts->SetMinNCrossedRowsTPC(70);
 823     esdTrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
 824   }
 825   else {
 826     AliWarningClass(Form("Wrong value of the clusterCut parameter (%d), using cut on Nclusters",clusterCut));
 827     esdTrackCuts->SetMinNClustersTPC(70);
 828   }
 829   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
 830   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 831   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 832   // ITS
 833   esdTrackCuts->SetRequireITSRefit(kTRUE);
 834   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 835                                          AliESDtrackCuts::kAny);
 836   if(selPrimaries) {
 837     // 7*(0.0026+0.0050/pt^1.01)
 838     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
 839     esdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
 840   }
 841   esdTrackCuts->SetMaxDCAToVertexZ(2);
 842   esdTrackCuts->SetDCAToVertex2D(kFALSE);
 843   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 844
 845   esdTrackCuts->SetMaxChi2PerClusterITS(36);
 846
 847   return esdTrackCuts;
 848 }
 849
 850 //____________________________________________________________________
 851 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSPureSATrackCuts2009(Bool_t selPrimaries, Bool_t useForPid)
 852 {
 853   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks
 854
 855   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 856   esdTrackCuts->SetRequireITSPureStandAlone(kTRUE);
 857   esdTrackCuts->SetRequireITSRefit(kTRUE);
 858   esdTrackCuts->SetMinNClustersITS(4);
 859   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 860                                          AliESDtrackCuts::kAny);
 861   esdTrackCuts->SetMaxChi2PerClusterITS(1.);
 862
 863   if(selPrimaries) {
 864     // 7*(0.0085+0.0026/pt^1.55)
 865     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0595+0.0182/pt^1.55");
 866   }
 867   if(useForPid){
 868     esdTrackCuts->SetRequireITSPid(kTRUE);
 869   }
 870   return esdTrackCuts;
 871 }
 872
 873 //____________________________________________________________________
 874 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSPureSATrackCuts2010(Bool_t selPrimaries, Bool_t useForPid)
 875 {
 876   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks - pp 2010
 877
 878   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 879   esdTrackCuts->SetRequireITSPureStandAlone(kTRUE);
 880   esdTrackCuts->SetRequireITSRefit(kTRUE);
 881   esdTrackCuts->SetMinNClustersITS(4);
 882   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 883                                          AliESDtrackCuts::kAny);
 884   esdTrackCuts->SetMaxChi2PerClusterITS(2.5);
 885
 886   if(selPrimaries) {
 887     // 7*(0.0033+0.0045/pt^1.3)
 888     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0231+0.0315/pt^1.3");
 889   }
 890   if(useForPid){
 891     esdTrackCuts->SetRequireITSPid(kTRUE);
 892   }
 893   return esdTrackCuts;
 894 }
 895
 896 //____________________________________________________________________
 897 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSSATrackCuts2009(Bool_t selPrimaries, Bool_t useForPid)
 898 {
 899   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks
 900
 901   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 902   esdTrackCuts->SetRequireITSStandAlone(kTRUE);
 903   esdTrackCuts->SetRequireITSPureStandAlone(kFALSE);
 904   esdTrackCuts->SetRequireITSRefit(kTRUE);
 905   esdTrackCuts->SetMinNClustersITS(4);
 906   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 907                                          AliESDtrackCuts::kAny);
 908   esdTrackCuts->SetMaxChi2PerClusterITS(1.);
 909
 910   if(selPrimaries) {
 911     // 7*(0.0085+0.0026/pt^1.55)
 912     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0595+0.0182/pt^1.55");
 913   }
 914   if(useForPid){
 915     esdTrackCuts->SetRequireITSPid(kTRUE);
 916   }
 917   return esdTrackCuts;
 918 }
 919
 920 //____________________________________________________________________
 921 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSSATrackCuts2010(Bool_t selPrimaries, Bool_t useForPid)
 922 {
 923   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks --pp 2010
 924
 925   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 926   esdTrackCuts->SetRequireITSStandAlone(kTRUE);
 927   esdTrackCuts->SetRequireITSPureStandAlone(kFALSE);
 928   esdTrackCuts->SetRequireITSRefit(kTRUE);
 929   esdTrackCuts->SetMinNClustersITS(4);
 930   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 931                                          AliESDtrackCuts::kAny);
 932   esdTrackCuts->SetMaxChi2PerClusterITS(2.5);
 933
 934   if(selPrimaries) {
 935     // 7*(0.0033+0.0045/pt^1.3)
 936     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0231+0.0315/pt^1.3");
 937   }
 938   if(useForPid){
 939     esdTrackCuts->SetRequireITSPid(kTRUE);
 940   }
 941   return esdTrackCuts;
 942 }
 943
 944 //____________________________________________________________________
 945 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSSATrackCutsPbPb2010(Bool_t selPrimaries, Bool_t useForPid)
 946 {
 947   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks -- PbPb 2010
 948
 949   AliESDtrackCuts* esdTrackCuts = GetStandardITSSATrackCuts2010(selPrimaries, useForPid);
 950   esdTrackCuts->SetMaxNOfMissingITSPoints(1);
 951
 952   return esdTrackCuts;
 953 }
 954 //____________________________________________________________________
 955
 956 AliESDtrackCuts* AliESDtrackCuts::GetStandardV0DaughterCuts()
 957 {
 958   // creates a AliESDtrackCuts object and fills it with standard cuts for V0 daughters
 959   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 960   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 961   esdTrackCuts->SetMinNClustersTPC(70);
 962   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 963   return esdTrackCuts;
 964 }
 965
 966 //____________________________________________________________________
 967 Int_t AliESDtrackCuts::GetReferenceMultiplicity(const AliESDEvent* esd, Bool_t tpcOnly)
 968 {
 969   // Gets reference multiplicity following the standard cuts and a defined fiducial volume
 970   // tpcOnly = kTRUE -> consider TPC-only tracks
 971   //         = kFALSE -> consider global tracks
 972   //
 973   // DEPRECATED Use GetReferenceMultiplicity with the enum as second argument instead
 974
 975   if (!tpcOnly)
 976   {
 977     AliErrorClass("Not implemented for global tracks!");
 978     return -1;
 979   }
 980
 981   static AliESDtrackCuts* esdTrackCuts = 0;
 982   if (!esdTrackCuts)
 983   {
 984     esdTrackCuts = GetStandardTPCOnlyTrackCuts();
 985     esdTrackCuts->SetEtaRange(-0.8, 0.8);
 986     esdTrackCuts->SetPtRange(0.15);
 987   }
 988
 989   Int_t nTracks = esdTrackCuts->CountAcceptedTracks(esd);
 990
 991   return nTracks;
 992 }
 993
 994 //____________________________________________________________________
 995 Float_t AliESDtrackCuts::GetSigmaToVertex(const AliESDtrack* const esdTrack)
 996 {
 997   // Calculates the number of sigma to the vertex.
 998
 999   Float_t b[2];
1000   Float_t bRes[2];
1001   Float_t bCov[3];
1002   esdTrack->GetImpactParameters(b,bCov);
1003
1004   if (bCov[0]<=0 || bCov[2]<=0) {
1005     AliDebugClass(1, "Estimated b resolution lower or equal zero!");
1006     bCov[0]=0; bCov[2]=0;
1007   }
1008   bRes[0] = TMath::Sqrt(bCov[0]);
1009   bRes[1] = TMath::Sqrt(bCov[2]);
1010
1011   // -----------------------------------
1012   // How to get to a n-sigma cut?
1013   //
1014   // The accumulated statistics from 0 to d is
1015   //
1016   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
1017   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
1018   //
1019   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
1020   // Can this be expressed in a different way?
1021
1022   if (bRes[0] == 0 || bRes[1] ==0)
1023     return -1;
1024
1025   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
1026
1027   // work around precision problem
1028   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
1029   // 1e-15 corresponds to nsigma ~ 7.7
1030   if (TMath::Exp(-d * d / 2) < 1e-15)
1031     return 1000;
1032
1033   Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
1034   return nSigma;
1035 }
1036
1037 void AliESDtrackCuts::EnableNeededBranches(TTree* tree)
1038 {
1039   // enables the branches needed by AcceptTrack, for a list see comment of AcceptTrack
1040
1041   tree->SetBranchStatus("fTracks.fFlags", 1);
1042   tree->SetBranchStatus("fTracks.fITSncls", 1);
1043   tree->SetBranchStatus("fTracks.fTPCncls", 1);
1044   tree->SetBranchStatus("fTracks.fITSchi2", 1);
1045   tree->SetBranchStatus("fTracks.fTPCchi2", 1);
1046   tree->SetBranchStatus("fTracks.fC*", 1);
1047   tree->SetBranchStatus("fTracks.fD", 1);
1048   tree->SetBranchStatus("fTracks.fZ", 1);
1049   tree->SetBranchStatus("fTracks.fCdd", 1);
1050   tree->SetBranchStatus("fTracks.fCdz", 1);
1051   tree->SetBranchStatus("fTracks.fCzz", 1);
1052   tree->SetBranchStatus("fTracks.fP*", 1);
1053   tree->SetBranchStatus("fTracks.fR*", 1);
1054   tree->SetBranchStatus("fTracks.fKinkIndexes*", 1);
1055 }
1056
1057 //____________________________________________________________________
1058 Bool_t AliESDtrackCuts::AcceptTrack(const AliESDtrack* esdTrack)
1059 {
1060   //
1061   // figure out if the tracks survives all the track cuts defined
1062   //
1063   // the different quality parameter and kinematic values are first
1064   // retrieved from the track. then it is found out what cuts the
1065   // track did not survive and finally the cuts are imposed.
1066
1067   // this function needs the following branches:
1068   // fTracks.fFlags
1069   // fTracks.fITSncls
1070   // fTracks.fTPCncls
1071   // fTracks.fITSchi2
1072   // fTracks.fTPCchi2
1073   // fTracks.fC   //GetExternalCovariance
1074   // fTracks.fD   //GetImpactParameters
1075   // fTracks.fZ   //GetImpactParameters
1076   // fTracks.fCdd //GetImpactParameters
1077   // fTracks.fCdz //GetImpactParameters
1078   // fTracks.fCzz //GetImpactParameters
1079   // fTracks.fP   //GetPxPyPz
1080   // fTracks.fR   //GetMass
1081   // fTracks.fP   //GetMass
1082   // fTracks.fKinkIndexes
1083   //
1084   // esdEvent is only required for the MaxChi2TPCConstrainedVsGlobal
1085
1086   UInt_t status = esdTrack->GetStatus();
1087
1088   // getting quality parameters from the ESD track
1089   Int_t nClustersITS = esdTrack->GetITSclusters(0);
1090   Int_t nClustersTPC = -1;
1091   if(fCutRequireTPCStandAlone) {
1092     nClustersTPC = esdTrack->GetTPCNclsIter1();
1093   }
1094   else {
1095     nClustersTPC = esdTrack->GetTPCclusters(0);
1096   }
1097
1098   //Pt dependent NClusters Cut
1099   if(f1CutMinNClustersTPCPtDep) {
1100     if(esdTrack->Pt()<fCutMaxPtDepNClustersTPC)
1101       fCutMinNClusterTPC = f1CutMinNClustersTPCPtDep->Eval(esdTrack->Pt());
1102     else
1103       fCutMinNClusterTPC = f1CutMinNClustersTPCPtDep->Eval(fCutMaxPtDepNClustersTPC);
1104   }
1105
1106   Float_t nCrossedRowsTPC = esdTrack->GetTPCCrossedRows();
1107   Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1108   if (esdTrack->GetTPCNclsF()>0) {
1109     ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC / esdTrack->GetTPCNclsF();
1110   }
1111
1112   Int_t nClustersTPCShared = esdTrack->GetTPCnclsS();
1113   Float_t fracClustersTPCShared = -1.;
1114
1115   Float_t chi2PerClusterITS = -1;
1116   Float_t chi2PerClusterTPC = -1;
1117   if (nClustersITS!=0)
1118     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1119   if (nClustersTPC!=0) {
1120     if(fCutRequireTPCStandAlone) {
1121       chi2PerClusterTPC = esdTrack->GetTPCchi2Iter1()/Float_t(nClustersTPC);
1122     } else {
1123       chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
1124     }
1125     fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClustersTPC);
1126   }
1127
1128   Double_t extCov[15];
1129   esdTrack->GetExternalCovariance(extCov);
1130
1131   Float_t b[2];
1132   Float_t bCov[3];
1133   esdTrack->GetImpactParameters(b,bCov);
1134   if (bCov[0]<=0 || bCov[2]<=0) {
1135     AliDebug(1, "Estimated b resolution lower or equal zero!");
1136     bCov[0]=0; bCov[2]=0;
1137   }
1138
1139
1140   // set pt-dependent DCA cuts, if requested
1141   SetPtDepDCACuts(esdTrack->Pt());
1142
1143
1144   Float_t dcaToVertexXY = b[0];
1145   Float_t dcaToVertexZ = b[1];
1146
1147   Float_t dcaToVertex = -1;
1148
1149   if (fCutDCAToVertex2D)
1150   {
1151     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1152   }
1153   else
1154     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1155
1156   // getting the kinematic variables of the track
1157   // (assuming the mass is known)
1158   Double_t p[3];
1159   esdTrack->GetPxPyPz(p);
1160
1161   Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
1162   Float_t pt       = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2));
1163   Float_t energy   = TMath::Sqrt(TMath::Power(esdTrack->GetMass(),2) + TMath::Power(momentum,2));
1164
1165   //y-eta related calculations
1166   Float_t eta = -100.;
1167   Float_t y   = -100.;
1168   if((momentum != TMath::Abs(p[2]))&&(momentum != 0))
1169     eta = 0.5*TMath::Log((momentum + p[2])/(momentum - p[2]));
1170   if((energy != TMath::Abs(p[2]))&&(momentum != 0))
1171     y = 0.5*TMath::Log((energy + p[2])/(energy - p[2]));
1172
1173   if (extCov[14] < 0)
1174   {
1175     AliWarning(Form("GetSigma1Pt2() returns negative value for external covariance matrix element fC[14]: %f. Corrupted track information, track will not be accepted!", extCov[14]));
1176     return kFALSE;
1177   }
1178   Float_t relUncertainty1Pt = TMath::Sqrt(extCov[14])*pt;
1179
1180   //########################################################################
1181   // cut the track?
1182
1183   Bool_t cuts[kNCuts];
1184   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1185
1186   // track quality cuts
1187   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1188     cuts[0]=kTRUE;
1189   if (fCutRequireTPCStandAlone && (status&AliESDtrack::kTPCin)==0)
1190     cuts[1]=kTRUE;
1191   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1192     cuts[2]=kTRUE;
1193   if (nClustersTPC<fCutMinNClusterTPC)
1194     cuts[3]=kTRUE;
1195   if (nClustersITS<fCutMinNClusterITS)
1196     cuts[4]=kTRUE;
1197   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1198     cuts[5]=kTRUE;
1199   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1200     cuts[6]=kTRUE;
1201   if (extCov[0]  > fCutMaxC11)
1202     cuts[7]=kTRUE;
1203   if (extCov[2]  > fCutMaxC22)
1204     cuts[8]=kTRUE;
1205   if (extCov[5]  > fCutMaxC33)
1206     cuts[9]=kTRUE;
1207   if (extCov[9]  > fCutMaxC44)
1208     cuts[10]=kTRUE;
1209   if (extCov[14]  > fCutMaxC55)
1210     cuts[11]=kTRUE;
1211
1212   // cut 12 and 13 see below
1213
1214   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1215     cuts[14]=kTRUE;
1216   // track kinematics cut
1217   if((momentum < fPMin) || (momentum > fPMax))
1218     cuts[15]=kTRUE;
1219   if((pt < fPtMin) || (pt > fPtMax))
1220     cuts[16] = kTRUE;
1221   if((p[0] < fPxMin) || (p[0] > fPxMax))
1222     cuts[17] = kTRUE;
1223   if((p[1] < fPyMin) || (p[1] > fPyMax))
1224     cuts[18] = kTRUE;
1225   if((p[2] < fPzMin) || (p[2] > fPzMax))
1226     cuts[19] = kTRUE;
1227   if((eta < fEtaMin) || (eta > fEtaMax))
1228     cuts[20] = kTRUE;
1229   if((y < fRapMin) || (y > fRapMax))
1230     cuts[21] = kTRUE;
1231   if (fCutDCAToVertex2D && dcaToVertex > 1)
1232     cuts[22] = kTRUE;
1233   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1234     cuts[23] = kTRUE;
1235   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1236     cuts[24] = kTRUE;
1237   if (fCutDCAToVertex2D && fCutMinDCAToVertexXY > 0 && fCutMinDCAToVertexZ > 0 && dcaToVertexXY*dcaToVertexXY/fCutMinDCAToVertexXY/fCutMinDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMinDCAToVertexZ/fCutMinDCAToVertexZ < 1)
1238     cuts[25] = kTRUE;
1239   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) < fCutMinDCAToVertexXY)
1240     cuts[26] = kTRUE;
1241   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) < fCutMinDCAToVertexZ)
1242     cuts[27] = kTRUE;
1243
1244   for (Int_t i = 0; i < 3; i++)
1245     cuts[28+i] = !CheckITSClusterRequirement(fCutClusterRequirementITS[i], esdTrack->HasPointOnITSLayer(i*2), esdTrack->HasPointOnITSLayer(i*2+1));
1246
1247   if(fCutRequireITSStandAlone || fCutRequireITSpureSA){
1248     if ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin)){
1249       // TPC tracks
1250       cuts[31] = kTRUE;
1251     }else{
1252       // ITS standalone tracks
1253       if(fCutRequireITSStandAlone && !fCutRequireITSpureSA){
1254         if(status & AliESDtrack::kITSpureSA) cuts[31] = kTRUE;
1255       }else if(fCutRequireITSpureSA){
1256         if(!(status & AliESDtrack::kITSpureSA)) cuts[31] = kTRUE;
1257       }
1258     }
1259   }
1260
1261   if (relUncertainty1Pt > fCutMaxRel1PtUncertainty)
1262      cuts[32] = kTRUE;
1263
1264   if (!fCutAcceptSharedTPCClusters && nClustersTPCShared!=0)
1265     cuts[33] = kTRUE;
1266
1267   if (fracClustersTPCShared > fCutMaxFractionSharedTPCClusters)
1268     cuts[34] = kTRUE;
1269
1270   Int_t nITSPointsForPid=0;
1271   UChar_t clumap=esdTrack->GetITSClusterMap();
1272   for(Int_t i=2; i<6; i++){
1273     if(clumap&(1<<i)) ++nITSPointsForPid;
1274   }
1275   if(fCutRequireITSPid && nITSPointsForPid<3) cuts[35] = kTRUE;
1276
1277
1278   if (nCrossedRowsTPC<fCutMinNCrossedRowsTPC)
1279     cuts[36]=kTRUE;
1280   if (ratioCrossedRowsOverFindableClustersTPC<fCutMinRatioCrossedRowsOverFindableClustersTPC)
1281     cuts[37]=kTRUE;
1282
1283   Int_t nMissITSpts=0;
1284   Int_t idet,statusLay;
1285   Float_t xloc,zloc;
1286   for(Int_t iLay=0; iLay<6; iLay++){
1287     Bool_t retc=esdTrack->GetITSModuleIndexInfo(iLay,idet,statusLay,xloc,zloc);
1288     if(retc && statusLay==5) ++nMissITSpts;
1289   }
1290   if(nMissITSpts>fCutMaxMissingITSPoints) cuts[38] = kTRUE;
1291
1292   //kTOFout
1293   if (fCutRequireTOFout && (status&AliESDtrack::kTOFout)==0)
1294     cuts[40]=kTRUE;
1295
1296   // TOF signal Dz cut
1297   Float_t dxTOF = esdTrack->GetTOFsignalDx();
1298   Float_t dzTOF = esdTrack->GetTOFsignalDz();
1299   if (fFlagCutTOFdistance && (esdTrack->GetStatus() & AliESDtrack::kTOFout) == AliESDtrack::kTOFout){ // applying the TOF distance cut only if requested, and only on tracks that reached the TOF and where associated with a TOF hit
1300           if (fgBeamTypeFlag < 0) {  // the check on the beam type was not done yet
1301                   const AliESDEvent* event = esdTrack->GetESDEvent();
1302                   if (event){
1303                           TString beamTypeESD = event->GetBeamType();
1304                           AliDebug(2,Form("Beam type from ESD event = %s",beamTypeESD.Data()));
1305                           if (beamTypeESD.CompareTo("A-A",TString::kIgnoreCase) == 0){ // we are in PbPb collisions --> fgBeamTypeFlag will be set to 1, to apply the cut on TOF signal Dz
1306                                   fgBeamTypeFlag = 1;
1307                           }
1308                           else { // we are NOT in PbPb collisions --> fgBeamTypeFlag will be set to 0, to NOT apply the cu6 on TOF signal Dz
1309                                   fgBeamTypeFlag = 0;
1310                           }
1311                   }
1312                   else{
1313                           AliFatal("Beam type not available, but it is needed to apply the TOF cut!");
1314                   }
1315           }
1316
1317           if (fgBeamTypeFlag == 1){ // we are in PbPb collisions --> apply the cut on TOF signal Dz
1318                   Float_t radiusTOF = TMath::Sqrt(dxTOF*dxTOF + dzTOF*dzTOF);
1319                   AliDebug(3,Form("TOF check (with fCutTOFdistance = %f) --> dx = %f, dz = %f, radius = %f", fCutTOFdistance, dxTOF, dzTOF, radiusTOF));
1320                   if (radiusTOF > fCutTOFdistance){
1321                           AliDebug(2, Form("************* the radius is outside the range! %f > %f, the track will be skipped", radiusTOF, fCutTOFdistance));
1322                           cuts[41] = kTRUE;
1323                   }
1324           }
1325   }
1326
1327   Bool_t cut=kFALSE;
1328   for (Int_t i=0; i<kNCuts; i++)
1329     if (cuts[i]) {cut = kTRUE;}
1330
1331   // for performance evaluate the CPU intensive cuts only when the others have passed, and when they are requested
1332   Double_t chi2TPCConstrainedVsGlobal = -2;
1333   Float_t nSigmaToVertex = -2;
1334   if (!cut)
1335   {
1336     // getting the track to vertex parameters
1337     if (fCutSigmaToVertexRequired)
1338     {
1339       nSigmaToVertex = GetSigmaToVertex(esdTrack);
1340       if (nSigmaToVertex > fCutNsigmaToVertex && fCutSigmaToVertexRequired)
1341       {
1342         cuts[12] = kTRUE;
1343         cut = kTRUE;
1344       }
1345       // if n sigma could not be calculated
1346       if (nSigmaToVertex<0 && fCutSigmaToVertexRequired)
1347       {
1348         cuts[13] = kTRUE;
1349         cut = kTRUE;
1350       }
1351     }
1352
1353     // max chi2 TPC constrained vs global track only if track passed the other cut
1354     if (fCutMaxChi2TPCConstrainedVsGlobal < 1e9)
1355     {
1356       const AliESDEvent* esdEvent = esdTrack->GetESDEvent();
1357
1358       if (!esdEvent)
1359         AliFatal("fCutMaxChi2TPCConstrainedVsGlobal set but ESD event not set in AliESDTrack. Use AliESDTrack::SetESDEvent before calling AliESDtrackCuts.");
1360
1361       // get vertex
1362       const AliESDVertex* vertex = 0;
1363       if (fCutMaxChi2TPCConstrainedVsGlobalVertexType & kVertexTracks)
1364         vertex = esdEvent->GetPrimaryVertexTracks();
1365
1366       if ((!vertex || !vertex->GetStatus()) && fCutMaxChi2TPCConstrainedVsGlobalVertexType & kVertexSPD)
1367         vertex = esdEvent->GetPrimaryVertexSPD();
1368
1369       if ((!vertex || !vertex->GetStatus()) && fCutMaxChi2TPCConstrainedVsGlobalVertexType & kVertexTPC)
1370         vertex = esdEvent->GetPrimaryVertexTPC();
1371
1372       if (vertex->GetStatus())
1373         chi2TPCConstrainedVsGlobal = esdTrack->GetChi2TPCConstrainedVsGlobal(vertex);
1374
1375       if (chi2TPCConstrainedVsGlobal < 0 || chi2TPCConstrainedVsGlobal > fCutMaxChi2TPCConstrainedVsGlobal)
1376       {
1377         cuts[39] = kTRUE;
1378         cut = kTRUE;
1379       }
1380     }
1381   }
1382
1383   //########################################################################
1384   // filling histograms
1385   if (fHistogramsOn) {
1386     fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n tracks")));
1387     if (cut)
1388       fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut tracks")));
1389
1390     for (Int_t i=0; i<kNCuts; i++) {
1391       if (fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i]) < 1)
1392         AliFatal(Form("Inconsistency! Cut %d with name %s not found", i, fgkCutNames[i]));
1393
1394       if (cuts[i])
1395         fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
1396
1397       for (Int_t j=i; j<kNCuts; j++) {
1398         if (cuts[i] && cuts[j]) {
1399           Float_t xC = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
1400           Float_t yC = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
1401           fhCutCorrelation->Fill(xC, yC);
1402         }
1403       }
1404     }
1405   }
1406
1407   // now we loop over the filling of the histograms twice: once "before" the cut, once "after"
1408   // the code is not in a function due to too many local variables that would need to be passed
1409
1410   for (Int_t id = 0; id < 2; id++)
1411   {
1412     // id = 0 --> before cut
1413     // id = 1 --> after cut
1414
1415     if (fHistogramsOn)
1416     {
1417       fhNClustersITS[id]->Fill(nClustersITS);
1418       fhNClustersTPC[id]->Fill(nClustersTPC);
1419       fhNSharedClustersTPC[id]->Fill(nClustersTPCShared);
1420       fhNCrossedRowsTPC[id]->Fill(nCrossedRowsTPC);
1421       fhRatioCrossedRowsOverFindableClustersTPC[id]->Fill(ratioCrossedRowsOverFindableClustersTPC);
1422       fhChi2PerClusterITS[id]->Fill(chi2PerClusterITS);
1423       fhChi2PerClusterTPC[id]->Fill(chi2PerClusterTPC);
1424       fhChi2TPCConstrainedVsGlobal[id]->Fill(chi2TPCConstrainedVsGlobal);
1425       fhNClustersForITSPID[id]->Fill(nITSPointsForPid);
1426       fhNMissingITSPoints[id]->Fill(nMissITSpts);
1427
1428       fhC11[id]->Fill(extCov[0]);
1429       fhC22[id]->Fill(extCov[2]);
1430       fhC33[id]->Fill(extCov[5]);
1431       fhC44[id]->Fill(extCov[9]);
1432       fhC55[id]->Fill(extCov[14]);
1433
1434       fhRel1PtUncertainty[id]->Fill(relUncertainty1Pt);
1435
1436       fhPt[id]->Fill(pt);
1437       fhEta[id]->Fill(eta);
1438       fhTOFdistance[id]->Fill(dxTOF, dzTOF);
1439
1440       Float_t bRes[2];
1441       bRes[0] = TMath::Sqrt(bCov[0]);
1442       bRes[1] = TMath::Sqrt(bCov[2]);
1443
1444       fhDZ[id]->Fill(b[1]);
1445       fhDXY[id]->Fill(b[0]);
1446       fhDXYDZ[id]->Fill(dcaToVertex);
1447       fhDXYvsDZ[id]->Fill(b[1],b[0]);
1448
1449       if (bRes[0]!=0 && bRes[1]!=0) {
1450         fhDZNormalized[id]->Fill(b[1]/bRes[1]);
1451         fhDXYNormalized[id]->Fill(b[0]/bRes[0]);
1452         fhDXYvsDZNormalized[id]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
1453         fhNSigmaToVertex[id]->Fill(nSigmaToVertex);
1454       }
1455     }
1456
1457     // cut the track
1458     if (cut)
1459       return kFALSE;
1460   }
1461
1462   return kTRUE;
1463 }
1464
1465 //____________________________________________________________________
1466 Bool_t AliESDtrackCuts::CheckITSClusterRequirement(ITSClusterRequirement req, Bool_t clusterL1, Bool_t clusterL2)
1467 {
1468   // checks if the cluster requirement is fullfilled (in this case: return kTRUE)
1469
1470   switch (req)
1471   {
1472         case kOff:        return kTRUE;
1473         case kNone:       return !clusterL1 && !clusterL2;
1474         case kAny:        return clusterL1 || clusterL2;
1475         case kFirst:      return clusterL1;
1476         case kOnlyFirst:  return clusterL1 && !clusterL2;
1477         case kSecond:     return clusterL2;
1478         case kOnlySecond: return clusterL2 && !clusterL1;
1479         case kBoth:       return clusterL1 && clusterL2;
1480   }
1481
1482   return kFALSE;
1483 }
1484
1485 //____________________________________________________________________
1486 AliESDtrack* AliESDtrackCuts::GetTPCOnlyTrack(const AliESDEvent* esd, Int_t iTrack)
1487 {
1488   // Utility function to create a TPC only track from the given esd track
1489   //
1490   // IMPORTANT: The track has to be deleted by the user
1491   //
1492   // NB. most of the functionality to get a TPC only track from an ESD track is in AliESDtrack, where it should be
1493   // there are only missing propagations here that are needed for old data
1494   // this function will therefore become obsolete
1495   //
1496   // adapted from code provided by CKB
1497
1498   if (!esd->GetPrimaryVertexTPC())
1499     return 0; // No TPC vertex no TPC tracks
1500
1501   if(!esd->GetPrimaryVertexTPC()->GetStatus())
1502     return 0; // TPC Vertex is created by default in AliESDEvent, do not use in this case
1503
1504   AliESDtrack* track = esd->GetTrack(iTrack);
1505   if (!track)
1506     return 0;
1507
1508   AliESDtrack *tpcTrack = new AliESDtrack();
1509
1510   // only true if we have a tpc track
1511   if (!track->FillTPCOnlyTrack(*tpcTrack))
1512   {
1513     delete tpcTrack;
1514     return 0;
1515   }
1516
1517   return tpcTrack;
1518 }
1519
1520 //____________________________________________________________________
1521 TObjArray* AliESDtrackCuts::GetAcceptedTracks(const AliESDEvent* esd, Bool_t bTPC)
1522 {
1523   //
1524   // returns an array of all tracks that pass the cuts
1525   // or an array of TPC only tracks (propagated to the TPC vertex during reco)
1526   // tracks that pass the cut
1527   //
1528   // NOTE: List has to be deleted by the user
1529
1530   TObjArray* acceptedTracks = new TObjArray();
1531
1532   // loop over esd tracks
1533   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
1534     if(bTPC){
1535       if(!esd->GetPrimaryVertexTPC())return acceptedTracks; // No TPC vertex no TPC tracks
1536       if(!esd->GetPrimaryVertexTPC()->GetStatus())return acceptedTracks; // No proper TPC vertex, only the default
1537
1538       AliESDtrack *tpcTrack = GetTPCOnlyTrack(esd, iTrack);
1539       if (!tpcTrack)
1540         continue;
1541
1542       if (AcceptTrack(tpcTrack)) {
1543         acceptedTracks->Add(tpcTrack);
1544       }
1545       else
1546         delete tpcTrack;
1547     }
1548     else
1549     {
1550       AliESDtrack* track = esd->GetTrack(iTrack);
1551       if(AcceptTrack(track))
1552         acceptedTracks->Add(track);
1553     }
1554   }
1555   if(bTPC)acceptedTracks->SetOwner(kTRUE);
1556   return acceptedTracks;
1557 }
1558
1559 //____________________________________________________________________
1560 Int_t AliESDtrackCuts::CountAcceptedTracks(const AliESDEvent* const esd)
1561 {
1562   //
1563   // returns an the number of tracks that pass the cuts
1564   //
1565
1566   Int_t count = 0;
1567
1568   // loop over esd tracks
1569   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
1570     AliESDtrack* track = esd->GetTrack(iTrack);
1571     if (AcceptTrack(track))
1572       count++;
1573   }
1574
1575   return count;
1576 }
1577
1578 //____________________________________________________________________
1579  void AliESDtrackCuts::DefineHistograms(Int_t color) {
1580    //
1581    // diagnostics histograms are defined
1582    //
1583
1584    fHistogramsOn=kTRUE;
1585
1586    Bool_t oldStatus = TH1::AddDirectoryStatus();
1587    TH1::AddDirectory(kFALSE);
1588
1589    //###################################################################################
1590    // defining histograms
1591
1592    fhCutStatistics = new TH1F("cut_statistics","cut statistics",kNCuts+4,-0.5,kNCuts+3.5);
1593
1594    fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
1595    fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
1596
1597    fhCutCorrelation = new TH2F("cut_correlation","cut correlation",kNCuts,-0.5,kNCuts-0.5,kNCuts,-0.5,kNCuts-0.5);;
1598
1599    for (Int_t i=0; i<kNCuts; i++) {
1600      fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fgkCutNames[i]);
1601      fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fgkCutNames[i]);
1602      fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fgkCutNames[i]);
1603    }
1604
1605   fhCutStatistics  ->SetLineColor(color);
1606   fhCutCorrelation ->SetLineColor(color);
1607   fhCutStatistics  ->SetLineWidth(2);
1608   fhCutCorrelation ->SetLineWidth(2);
1609
1610   for (Int_t i=0; i<2; i++) {
1611     fhNClustersITS[i]        = new TH1F("nClustersITS"    ,"",8,-0.5,7.5);
1612     fhNClustersTPC[i]        = new TH1F("nClustersTPC"     ,"",165,-0.5,164.5);
1613     fhNSharedClustersTPC[i]  = new TH1F("nSharedClustersTPC"     ,"",165,-0.5,164.5);
1614     fhNCrossedRowsTPC[i]     = new TH1F("nCrossedRowsTPC"  ,"",165,-0.5,164.5);
1615     fhRatioCrossedRowsOverFindableClustersTPC[i]     = new TH1F("ratioCrossedRowsOverFindableClustersTPC"  ,"",60,0,1.5);
1616     fhChi2PerClusterITS[i]   = new TH1F("chi2PerClusterITS","",500,0,10);
1617     fhChi2PerClusterTPC[i]   = new TH1F("chi2PerClusterTPC","",500,0,10);
1618     fhChi2TPCConstrainedVsGlobal[i]   = new TH1F("chi2TPCConstrainedVsGlobal","",600,-2,50);
1619     fhNClustersForITSPID[i]  = new TH1F("nPointsForITSpid","",5,-0.5,4.5);
1620     fhNMissingITSPoints[i]   = new TH1F("nMissingITSClusters","",7,-0.5,6.5);
1621
1622     fhC11[i]                 = new TH1F("covMatrixDiagonal11","",2000,0,20);
1623     fhC22[i]                 = new TH1F("covMatrixDiagonal22","",2000,0,20);
1624     fhC33[i]                 = new TH1F("covMatrixDiagonal33","",1000,0,0.1);
1625     fhC44[i]                 = new TH1F("covMatrixDiagonal44","",1000,0,0.1);
1626     fhC55[i]                 = new TH1F("covMatrixDiagonal55","",1000,0,5);
1627
1628     fhRel1PtUncertainty[i]   = new TH1F("rel1PtUncertainty","",1000,0,5);
1629
1630     fhDXY[i]                 = new TH1F("dXY"    ,"",500,-10,10);
1631     fhDZ[i]                  = new TH1F("dZ"     ,"",500,-10,10);
1632     fhDXYDZ[i]               = new TH1F("dXYDZ"  ,"",500,0,10);
1633     fhDXYvsDZ[i]             = new TH2F("dXYvsDZ","",200,-10,10,200,-10,10);
1634
1635     fhDXYNormalized[i]       = new TH1F("dXYNormalized"    ,"",500,-10,10);
1636     fhDZNormalized[i]        = new TH1F("dZNormalized"     ,"",500,-10,10);
1637     fhDXYvsDZNormalized[i]   = new TH2F("dXYvsDZNormalized","",200,-10,10,200,-10,10);
1638
1639     fhNSigmaToVertex[i]      = new TH1F("nSigmaToVertex","",500,0,10);
1640
1641     fhPt[i]                  = new TH1F("pt"     ,"p_{T} distribution;p_{T} (GeV/c)", 800, 0.0, 10.0);
1642     fhEta[i]                 = new TH1F("eta"     ,"#eta distribution;#eta",40,-2.0,2.0);
1643     fhTOFdistance[i]         = new TH2F("TOFdistance"     ,"TOF distance;dx (cm};dz (cm)", 150, -15, 15, 150, -15, 15);
1644
1645     fhNClustersITS[i]->SetTitle("n ITS clusters");
1646     fhNClustersTPC[i]->SetTitle("n TPC clusters");
1647     fhNSharedClustersTPC[i]->SetTitle("n TPC shared clusters");
1648     fhChi2PerClusterITS[i]->SetTitle("#Chi^{2} per ITS cluster");
1649     fhChi2PerClusterTPC[i]->SetTitle("#Chi^{2} per TPC cluster");
1650     fhChi2TPCConstrainedVsGlobal[i]->SetTitle("#Chi^{2} TPC constrained track vs global track");
1651     fhNClustersForITSPID[i]->SetTitle("n ITS points for PID");
1652     fhNMissingITSPoints[i]->SetTitle("n ITS layers with missing cluster");
1653
1654     fhC11[i]->SetTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
1655     fhC22[i]->SetTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
1656     fhC33[i]->SetTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
1657     fhC44[i]->SetTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
1658     fhC55[i]->SetTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
1659
1660     fhRel1PtUncertainty[i]->SetTitle("rel. uncertainty of 1/p_{T}");
1661
1662     fhDXY[i]->SetXTitle("transverse impact parameter (cm)");
1663     fhDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
1664     fhDXYDZ[i]->SetTitle("absolute impact parameter;sqrt(dXY**2 + dZ**2)  (cm)");
1665     fhDXYvsDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
1666     fhDXYvsDZ[i]->SetYTitle("transverse impact parameter (cm)");
1667
1668     fhDXYNormalized[i]->SetTitle("normalized trans impact par (n#sigma)");
1669     fhDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
1670     fhDXYvsDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
1671     fhDXYvsDZNormalized[i]->SetYTitle("normalized trans impact par (n#sigma)");
1672     fhNSigmaToVertex[i]->SetTitle("n #sigma to vertex");
1673
1674     fhNClustersITS[i]->SetLineColor(color);   fhNClustersITS[i]->SetLineWidth(2);
1675     fhNClustersTPC[i]->SetLineColor(color);   fhNClustersTPC[i]->SetLineWidth(2);
1676     fhNSharedClustersTPC[i]->SetLineColor(color);  fhNSharedClustersTPC[i]->SetLineWidth(2);
1677     fhChi2PerClusterITS[i]->SetLineColor(color);   fhChi2PerClusterITS[i]->SetLineWidth(2);
1678     fhChi2PerClusterTPC[i]->SetLineColor(color);   fhChi2PerClusterTPC[i]->SetLineWidth(2);
1679     fhChi2TPCConstrainedVsGlobal[i]->SetLineColor(color);   fhChi2TPCConstrainedVsGlobal[i]->SetLineWidth(2);
1680     fhNClustersForITSPID[i]->SetLineColor(color);  fhNClustersForITSPID[i]->SetLineWidth(2);
1681     fhNMissingITSPoints[i]->SetLineColor(color);   fhNMissingITSPoints[i]->SetLineWidth(2);
1682
1683     fhC11[i]->SetLineColor(color);   fhC11[i]->SetLineWidth(2);
1684     fhC22[i]->SetLineColor(color);   fhC22[i]->SetLineWidth(2);
1685     fhC33[i]->SetLineColor(color);   fhC33[i]->SetLineWidth(2);
1686     fhC44[i]->SetLineColor(color);   fhC44[i]->SetLineWidth(2);
1687     fhC55[i]->SetLineColor(color);   fhC55[i]->SetLineWidth(2);
1688
1689     fhRel1PtUncertainty[i]->SetLineColor(color); fhRel1PtUncertainty[i]->SetLineWidth(2);
1690
1691     fhDXY[i]->SetLineColor(color);   fhDXY[i]->SetLineWidth(2);
1692     fhDZ[i]->SetLineColor(color);    fhDZ[i]->SetLineWidth(2);
1693     fhDXYDZ[i]->SetLineColor(color); fhDXYDZ[i]->SetLineWidth(2);
1694
1695     fhDXYNormalized[i]->SetLineColor(color);   fhDXYNormalized[i]->SetLineWidth(2);
1696     fhDZNormalized[i]->SetLineColor(color);    fhDZNormalized[i]->SetLineWidth(2);
1697     fhNSigmaToVertex[i]->SetLineColor(color);  fhNSigmaToVertex[i]->SetLineWidth(2);
1698   }
1699
1700   // The number of sigmas to the vertex is per definition gaussian
1701   ffDTheoretical = new TF1("nSigmaToVertexTheoretical","([0]/2.506628274)*exp(-(x**2)/2)",0,50);
1702   ffDTheoretical->SetParameter(0,1);
1703
1704   TH1::AddDirectory(oldStatus);
1705 }
1706
1707 //____________________________________________________________________
1708 Bool_t AliESDtrackCuts::LoadHistograms(const Char_t* dir)
1709 {
1710   //
1711   // loads the histograms from a file
1712   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
1713   //
1714
1715   if (!dir)
1716     dir = GetName();
1717
1718   if (!gDirectory->cd(dir))
1719     return kFALSE;
1720
1721   ffDTheoretical = dynamic_cast<TF1*> (gDirectory->Get("nSigmaToVertexTheory"));
1722
1723   fhCutStatistics = dynamic_cast<TH1F*> (gDirectory->Get("cut_statistics"));
1724   fhCutCorrelation = dynamic_cast<TH2F*> (gDirectory->Get("cut_correlation"));
1725
1726   for (Int_t i=0; i<2; i++) {
1727     if (i==0)
1728     {
1729       gDirectory->cd("before_cuts");
1730     }
1731     else
1732       gDirectory->cd("after_cuts");
1733
1734     fhNClustersITS[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersITS"     ));
1735     fhNClustersTPC[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersTPC"     ));
1736     fhNSharedClustersTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get("nSharedClustersTPC"     ));
1737     fhNCrossedRowsTPC[i]   = dynamic_cast<TH1F*> (gDirectory->Get("nCrossedRowsTPC"  ));
1738     fhRatioCrossedRowsOverFindableClustersTPC[i]   = dynamic_cast<TH1F*> (gDirectory->Get("ratioCrossedRowsOverFindableClustersTPC"  ));
1739     fhChi2PerClusterITS[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterITS"));
1740     fhChi2PerClusterTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterTPC"));
1741     fhChi2TPCConstrainedVsGlobal[i] = dynamic_cast<TH1F*> (gDirectory->Get("fhChi2TPCConstrainedVsGlobal"));
1742     fhNClustersForITSPID[i] = dynamic_cast<TH1F*> (gDirectory->Get("nPointsForITSpid"));
1743     fhNMissingITSPoints[i] = dynamic_cast<TH1F*> (gDirectory->Get("nMissingITSClusters"));
1744
1745     fhC11[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal11"));
1746     fhC22[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal22"));
1747     fhC33[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal33"));
1748     fhC44[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal44"));
1749     fhC55[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal55"));
1750
1751     fhRel1PtUncertainty[i] = dynamic_cast<TH1F*> (gDirectory->Get("rel1PtUncertainty"));
1752
1753     fhDXY[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXY"    ));
1754     fhDZ[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZ"     ));
1755     fhDXYDZ[i] =   dynamic_cast<TH1F*> (gDirectory->Get("dXYDZ"));
1756     fhDXYvsDZ[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZ"));
1757
1758     fhDXYNormalized[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXYNormalized"    ));
1759     fhDZNormalized[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZNormalized"     ));
1760     fhDXYvsDZNormalized[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZNormalized"));
1761     fhNSigmaToVertex[i] = dynamic_cast<TH1F*> (gDirectory->Get("nSigmaToVertex"));
1762
1763     fhPt[i] = dynamic_cast<TH1F*> (gDirectory->Get("pt"));
1764     fhEta[i] = dynamic_cast<TH1F*> (gDirectory->Get("eta"));
1765     fhTOFdistance[i] = dynamic_cast<TH2F*> (gDirectory->Get("TOFdistance"));
1766
1767     gDirectory->cd("../");
1768   }
1769
1770   gDirectory->cd("..");
1771
1772   return kTRUE;
1773 }
1774
1775 //____________________________________________________________________
1776 void AliESDtrackCuts::SaveHistograms(const Char_t* dir) {
1777   //
1778   // saves the histograms in a directory (dir)
1779   //
1780
1781   if (!fHistogramsOn) {
1782     AliDebug(0, "Histograms not on - cannot save histograms!!!");
1783     return;
1784   }
1785
1786   if (!dir)
1787     dir = GetName();
1788
1789   gDirectory->mkdir(dir);
1790   gDirectory->cd(dir);
1791
1792   gDirectory->mkdir("before_cuts");
1793   gDirectory->mkdir("after_cuts");
1794
1795   // a factor of 2 is needed since n sigma is positive
1796   ffDTheoretical->SetParameter(0,2*fhNSigmaToVertex[0]->Integral("width"));
1797   ffDTheoretical->Write("nSigmaToVertexTheory");
1798
1799   fhCutStatistics->Write();
1800   fhCutCorrelation->Write();
1801
1802   for (Int_t i=0; i<2; i++) {
1803     if (i==0)
1804       gDirectory->cd("before_cuts");
1805     else
1806       gDirectory->cd("after_cuts");
1807
1808     fhNClustersITS[i]        ->Write();
1809     fhNClustersTPC[i]        ->Write();
1810     fhNSharedClustersTPC[i]  ->Write();
1811     fhNCrossedRowsTPC[i]     ->Write();
1812     fhRatioCrossedRowsOverFindableClustersTPC[i]     ->Write();
1813     fhChi2PerClusterITS[i]   ->Write();
1814     fhChi2PerClusterTPC[i]   ->Write();
1815     fhChi2TPCConstrainedVsGlobal[i]   ->Write();
1816     fhNClustersForITSPID[i]  ->Write();
1817     fhNMissingITSPoints[i]   ->Write();
1818
1819     fhC11[i]                 ->Write();
1820     fhC22[i]                 ->Write();
1821     fhC33[i]                 ->Write();
1822     fhC44[i]                 ->Write();
1823     fhC55[i]                 ->Write();
1824
1825     fhRel1PtUncertainty[i]   ->Write();
1826
1827     fhDXY[i]                 ->Write();
1828     fhDZ[i]                  ->Write();
1829     fhDXYDZ[i]               ->Write();
1830     fhDXYvsDZ[i]             ->Write();
1831
1832     fhDXYNormalized[i]       ->Write();
1833     fhDZNormalized[i]        ->Write();
1834     fhDXYvsDZNormalized[i]   ->Write();
1835     fhNSigmaToVertex[i]      ->Write();
1836
1837     fhPt[i]                  ->Write();
1838     fhEta[i]                 ->Write();
1839     fhTOFdistance[i]         ->Write();
1840
1841     gDirectory->cd("../");
1842   }
1843
1844   gDirectory->cd("../");
1845 }
1846
1847 //____________________________________________________________________
1848 void AliESDtrackCuts::DrawHistograms()
1849 {
1850   // draws some histograms
1851
1852   TCanvas* canvas1 = new TCanvas(Form("%s_1", GetName()), "Track Quality Results1", 800, 800);
1853   canvas1->Divide(2, 2);
1854
1855   canvas1->cd(1);
1856   fhNClustersTPC[0]->SetStats(kFALSE);
1857   fhNClustersTPC[0]->Draw();
1858
1859   canvas1->cd(2);
1860   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1861   fhChi2PerClusterTPC[0]->Draw();
1862
1863   canvas1->cd(3);
1864   fhNSigmaToVertex[0]->SetStats(kFALSE);
1865   fhNSigmaToVertex[0]->GetXaxis()->SetRangeUser(0, 10);
1866   fhNSigmaToVertex[0]->Draw();
1867
1868   canvas1->SaveAs(Form("%s_%s.gif", GetName(), canvas1->GetName()));
1869
1870   TCanvas* canvas2 = new TCanvas(Form("%s_2", GetName()), "Track Quality Results2", 1200, 800);
1871   canvas2->Divide(3, 2);
1872
1873   canvas2->cd(1);
1874   fhC11[0]->SetStats(kFALSE);
1875   gPad->SetLogy();
1876   fhC11[0]->Draw();
1877
1878   canvas2->cd(2);
1879   fhC22[0]->SetStats(kFALSE);
1880   gPad->SetLogy();
1881   fhC22[0]->Draw();
1882
1883   canvas2->cd(3);
1884   fhC33[0]->SetStats(kFALSE);
1885   gPad->SetLogy();
1886   fhC33[0]->Draw();
1887
1888   canvas2->cd(4);
1889   fhC44[0]->SetStats(kFALSE);
1890   gPad->SetLogy();
1891   fhC44[0]->Draw();
1892
1893   canvas2->cd(5);
1894   fhC55[0]->SetStats(kFALSE);
1895   gPad->SetLogy();
1896   fhC55[0]->Draw();
1897
1898   canvas2->cd(6);
1899   fhRel1PtUncertainty[0]->SetStats(kFALSE);
1900   gPad->SetLogy();
1901   fhRel1PtUncertainty[0]->Draw();
1902
1903   canvas2->SaveAs(Form("%s_%s.gif", GetName(), canvas2->GetName()));
1904
1905   TCanvas* canvas3 = new TCanvas(Form("%s_3", GetName()), "Track Quality Results3", 1200, 800);
1906   canvas3->Divide(3, 2);
1907
1908   canvas3->cd(1);
1909   fhDXY[0]->SetStats(kFALSE);
1910   gPad->SetLogy();
1911   fhDXY[0]->Draw();
1912
1913   canvas3->cd(2);
1914   fhDZ[0]->SetStats(kFALSE);
1915   gPad->SetLogy();
1916   fhDZ[0]->Draw();
1917
1918   canvas3->cd(3);
1919   fhDXYvsDZ[0]->SetStats(kFALSE);
1920   gPad->SetLogz();
1921   gPad->SetRightMargin(0.15);
1922   fhDXYvsDZ[0]->Draw("COLZ");
1923
1924   canvas3->cd(4);
1925   fhDXYNormalized[0]->SetStats(kFALSE);
1926   gPad->SetLogy();
1927   fhDXYNormalized[0]->Draw();
1928
1929   canvas3->cd(5);
1930   fhDZNormalized[0]->SetStats(kFALSE);
1931   gPad->SetLogy();
1932   fhDZNormalized[0]->Draw();
1933
1934   canvas3->cd(6);
1935   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1936   gPad->SetLogz();
1937   gPad->SetRightMargin(0.15);
1938   fhDXYvsDZNormalized[0]->Draw("COLZ");
1939
1940   canvas3->SaveAs(Form("%s_%s.gif", GetName(), canvas3->GetName()));
1941
1942   TCanvas* canvas4 = new TCanvas(Form("%s_4", GetName()), "Track Quality Results4", 800, 500);
1943   canvas4->Divide(2, 1);
1944
1945   canvas4->cd(1);
1946   fhCutStatistics->SetStats(kFALSE);
1947   fhCutStatistics->LabelsOption("v");
1948   gPad->SetBottomMargin(0.3);
1949   fhCutStatistics->Draw();
1950
1951   canvas4->cd(2);
1952   fhCutCorrelation->SetStats(kFALSE);
1953   fhCutCorrelation->LabelsOption("v");
1954   gPad->SetBottomMargin(0.3);
1955   gPad->SetLeftMargin(0.3);
1956   fhCutCorrelation->Draw("COLZ");
1957
1958   canvas4->SaveAs(Form("%s_%s.gif", GetName(), canvas4->GetName()));
1959
1960   /*canvas->cd(1);
1961   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1962   fhDXYvsDZNormalized[0]->DrawCopy("COLZ");
1963
1964   canvas->cd(2);
1965   fhNClustersTPC[0]->SetStats(kFALSE);
1966   fhNClustersTPC[0]->DrawCopy();
1967
1968   canvas->cd(3);
1969   fhChi2PerClusterITS[0]->SetStats(kFALSE);
1970   fhChi2PerClusterITS[0]->DrawCopy();
1971   fhChi2PerClusterITS[1]->SetLineColor(2);
1972   fhChi2PerClusterITS[1]->DrawCopy("SAME");
1973
1974   canvas->cd(4);
1975   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1976   fhChi2PerClusterTPC[0]->DrawCopy();
1977   fhChi2PerClusterTPC[1]->SetLineColor(2);
1978   fhChi2PerClusterTPC[1]->DrawCopy("SAME");*/
1979 }
1980 //--------------------------------------------------------------------------
1981 void AliESDtrackCuts::SetPtDepDCACuts(Double_t pt) {
1982   //
1983   // set the pt-dependent DCA cuts
1984   //
1985
1986   if(f1CutMaxDCAToVertexXYPtDep) {
1987      fCutMaxDCAToVertexXY=f1CutMaxDCAToVertexXYPtDep->Eval(pt);
1988   }
1989
1990   if(f1CutMaxDCAToVertexZPtDep) {
1991     fCutMaxDCAToVertexZ=f1CutMaxDCAToVertexZPtDep->Eval(pt);
1992   }
1993
1994   if(f1CutMinDCAToVertexXYPtDep) {
1995     fCutMinDCAToVertexXY=f1CutMinDCAToVertexXYPtDep->Eval(pt);
1996   }
1997
1998   if(f1CutMinDCAToVertexZPtDep) {
1999     fCutMinDCAToVertexZ=f1CutMinDCAToVertexZPtDep->Eval(pt);
2000   }
2001
2002
2003   return;
2004 }
2005
2006
2007
2008 //--------------------------------------------------------------------------
2009 Bool_t AliESDtrackCuts::CheckPtDepDCA(TString dist,Bool_t print) const {
2010   //
2011   // Check the correctness of the string syntax
2012   //
2013   Bool_t retval=kTRUE;
2014
2015   if(!dist.Contains("pt")) {
2016     if(print) AliError("string must contain \"pt\"");
2017     retval= kFALSE;
2018   }
2019   return retval;
2020 }
2021
2022  void AliESDtrackCuts::SetMaxDCAToVertexXYPtDep(const char *dist){
2023
2024    if(f1CutMaxDCAToVertexXYPtDep){
2025      delete f1CutMaxDCAToVertexXYPtDep;
2026      // resetiing both
2027      f1CutMaxDCAToVertexXYPtDep = 0;
2028      fCutMaxDCAToVertexXYPtDep = "";
2029    }
2030    if(!CheckPtDepDCA(dist,kTRUE)){
2031      return;
2032    }
2033    fCutMaxDCAToVertexXYPtDep = dist;
2034    TString tmp(dist);
2035    tmp.ReplaceAll("pt","x");
2036    f1CutMaxDCAToVertexXYPtDep = new TFormula("f1CutMaxDCAToVertexXYPtDep",tmp.Data());
2037
2038 }
2039
2040  void AliESDtrackCuts::SetMaxDCAToVertexZPtDep(const char *dist){
2041
2042
2043    if(f1CutMaxDCAToVertexZPtDep){
2044      delete f1CutMaxDCAToVertexZPtDep;
2045      // resetiing both
2046      f1CutMaxDCAToVertexZPtDep = 0;
2047      fCutMaxDCAToVertexZPtDep = "";
2048    }
2049    if(!CheckPtDepDCA(dist,kTRUE))return;
2050
2051    fCutMaxDCAToVertexZPtDep = dist;
2052    TString tmp(dist);
2053    tmp.ReplaceAll("pt","x");
2054    f1CutMaxDCAToVertexZPtDep = new TFormula("f1CutMaxDCAToVertexZPtDep",tmp.Data());
2055
2056
2057 }
2058
2059
2060  void AliESDtrackCuts::SetMinDCAToVertexXYPtDep(const char *dist){
2061
2062
2063    if(f1CutMinDCAToVertexXYPtDep){
2064      delete f1CutMinDCAToVertexXYPtDep;
2065      // resetiing both
2066      f1CutMinDCAToVertexXYPtDep = 0;
2067      fCutMinDCAToVertexXYPtDep = "";
2068    }
2069    if(!CheckPtDepDCA(dist,kTRUE))return;
2070
2071    fCutMinDCAToVertexXYPtDep = dist;
2072    TString tmp(dist);
2073    tmp.ReplaceAll("pt","x");
2074    f1CutMinDCAToVertexXYPtDep = new TFormula("f1CutMinDCAToVertexXYPtDep",tmp.Data());
2075
2076 }
2077
2078
2079  void AliESDtrackCuts::SetMinDCAToVertexZPtDep(const char *dist){
2080
2081
2082
2083    if(f1CutMinDCAToVertexZPtDep){
2084      delete f1CutMinDCAToVertexZPtDep;
2085      // resetiing both
2086      f1CutMinDCAToVertexZPtDep = 0;
2087      fCutMinDCAToVertexZPtDep = "";
2088    }
2089    if(!CheckPtDepDCA(dist,kTRUE))return;
2090    fCutMinDCAToVertexZPtDep = dist;
2091    TString tmp(dist);
2092    tmp.ReplaceAll("pt","x");
2093    f1CutMinDCAToVertexZPtDep = new TFormula("f1CutMinDCAToVertexZPtDep",tmp.Data());
2094 }
2095
2096 AliESDtrackCuts* AliESDtrackCuts::GetMultEstTrackCuts(MultEstTrackCuts cut)
2097 {
2098   // returns the multiplicity estimator track cuts objects to allow for user configuration
2099   // upon first call the objects are created
2100   //
2101   // the cut defined here correspond to GetStandardITSTPCTrackCuts2010 (apart from the one for without SPD)
2102
2103   if (!fgMultEstTrackCuts[kMultEstTrackCutGlobal])
2104   {
2105     // quality cut on ITS+TPC tracks
2106     fgMultEstTrackCuts[kMultEstTrackCutGlobal] = new AliESDtrackCuts();
2107     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetMinNClustersTPC(70);
2108     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetMaxChi2PerClusterTPC(4);
2109     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetAcceptKinkDaughters(kFALSE);
2110     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetRequireTPCRefit(kTRUE);
2111     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetRequireITSRefit(kTRUE);
2112     //multiplicity underestimate if we use global tracks with |eta| > 0.9
2113     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetEtaRange(-0.9, 0.9);
2114
2115     // quality cut on ITS_SA tracks (complementary to ITS+TPC)
2116     fgMultEstTrackCuts[kMultEstTrackCutITSSA] = new AliESDtrackCuts();
2117     fgMultEstTrackCuts[kMultEstTrackCutITSSA]->SetRequireITSRefit(kTRUE);
2118
2119     // primary selection for tracks with SPD hits
2120     fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD] = new AliESDtrackCuts();
2121     fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
2122     fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
2123     fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->SetMaxDCAToVertexZ(2);
2124
2125     // primary selection for tracks w/o SPD hits
2126     fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD] = new AliESDtrackCuts();
2127     fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kNone);
2128     fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->SetMaxDCAToVertexXYPtDep("1.5*(0.0182+0.0350/pt^1.01)");
2129     fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->SetMaxDCAToVertexZ(2);
2130   }
2131
2132   return fgMultEstTrackCuts[cut];
2133 }
2134
2135 Int_t AliESDtrackCuts::GetReferenceMultiplicity(const AliESDEvent* esd, MultEstTrackType trackType, Float_t etaRange)
2136 {
2137   // Get multiplicity estimate based on TPC/ITS tracks and tracklets
2138   // Adapted for AliESDtrackCuts from a version developed by: Ruben Shahoyan, Anton Alkin, Arvinder Palaha
2139   //
2140   // Returns a negative value if no reliable estimate can be provided:
2141   //   -1 SPD vertex missing
2142   //   -2 SPD VertexerZ dispersion too large
2143   //   -3 Track vertex missing (not checked for kTracklets)
2144   //   -4 Distance between SPD and track vertices too large (not checked for kTracklets)
2145   //
2146   // WARNING This functions does not cut on the z vtx. Depending on the eta range requested, you need to restrict your z vertex range!
2147   //
2148   // Strategy for combined estimators
2149   //   1. Count global tracks and flag them
2150   //   2. Count ITSSA as complementaries for ITSTPC+ or as main tracks
2151   //   3. Count the complementary SPD tracklets
2152
2153   const AliESDVertex* vertices[2];
2154   vertices[0] = esd->GetPrimaryVertexSPD();
2155   vertices[1] = esd->GetPrimaryVertexTracks();
2156
2157   if (!vertices[0]->GetStatus())
2158   {
2159     AliDebugClass(AliLog::kDebug, "No SPD vertex. Not able to make a reliable multiplicity estimate.");
2160     return -1;
2161   }
2162
2163   if (vertices[0]->IsFromVertexerZ() && vertices[0]->GetDispersion() > 0.02)
2164   {
2165     AliDebugClass(AliLog::kDebug, "Vertexer z dispersion > 0.02. Not able to make a reliable multiplicity estimate.");
2166     return -2;
2167   }
2168
2169   Int_t multiplicityEstimate = 0;
2170
2171   // SPD tracklet-only estimate
2172   if (trackType == kTracklets)
2173   {
2174     const AliMultiplicity* spdmult = esd->GetMultiplicity();    // spd multiplicity object
2175     for (Int_t i=0; i<spdmult->GetNumberOfTracklets(); ++i)
2176     {
2177       if (TMath::Abs(spdmult->GetEta(i)) > etaRange)
2178         continue; // eta selection for tracklets
2179       multiplicityEstimate++;
2180     }
2181     return multiplicityEstimate;
2182   }
2183
2184   if (!vertices[1]->GetStatus())
2185   {
2186     AliDebugClass(AliLog::kDebug, "No track vertex. Not able to make a reliable multiplicity estimate.");
2187     return -3;
2188   }
2189
2190   // TODO value of displacement to be studied
2191   const Float_t maxDisplacement = 0.5;
2192   //check for displaced vertices
2193   Double_t displacement = TMath::Abs(vertices[0]->GetZ() - vertices[1]->GetZ());
2194   if (displacement > maxDisplacement)
2195   {
2196     AliDebugClass(AliLog::kDebug, Form("Displaced vertices %f > %f",displacement,maxDisplacement));
2197     return -4;
2198   }
2199
2200   // update eta range in track cuts
2201   GetMultEstTrackCuts(kMultEstTrackCutITSSA)->SetEtaRange(-etaRange, etaRange);
2202   GetMultEstTrackCuts(kMultEstTrackCutDCAwSPD)->SetEtaRange(-etaRange, etaRange);
2203   GetMultEstTrackCuts(kMultEstTrackCutDCAwoSPD)->SetEtaRange(-etaRange, etaRange);
2204
2205   //*******************************************************************************************************
2206   //set counters to initial zeros
2207   Int_t tracksITSTPC = 0;     //number of global tracks for a given event
2208   Int_t tracksITSSA = 0;      //number of ITS standalone tracks for a given event
2209   Int_t tracksITSTPCSA_complementary = 0; //number of ITS standalone tracks complementary to TPC for a given event
2210   Int_t trackletsITSTPC_complementary = 0;//number of SPD tracklets complementary to global/ITSSA tracks for a given events
2211   Int_t trackletsITSSA_complementary = 0; //number of SPD tracklets complementary to ITSSA tracks for a given event
2212
2213   const Int_t nESDTracks = esd->GetNumberOfTracks();
2214   Int_t highestID = 0;
2215
2216   // flags for secondary and rejected tracks
2217   const Int_t kRejBit = BIT(15); // set this bit in global tracks if it is rejected by a cut
2218   const Int_t kSecBit = BIT(16); // set this bit in global tracks if it is secondary according to a cut
2219
2220   for(Int_t itracks=0; itracks < nESDTracks; itracks++) {
2221     esd->GetTrack(itracks)->ResetBit(kSecBit|kRejBit); //reset bits used for flagging secondaries and rejected tracks in case they were changed before this analysis
2222   }
2223   const Int_t maxid = nESDTracks+1; // used to define bool array for check multiple associations of tracklets to one track. array starts at 0.
2224
2225   // bit mask for esd tracks, to check if multiple tracklets are associated to it
2226   TBits globalBits(maxid), pureITSBits(maxid);
2227   // why labels are used with the data? RS
2228   //  Bool_t globalBits[maxid], pureITSBits[maxid];
2229   //  for(Int_t i=0; i<maxid; i++){ // set all bools to false
2230   //      globalBits[i]=kFALSE;
2231   //      pureITSBits[i]=kFALSE;
2232   //  }
2233
2234   //*******************************************************************************************************
2235   // get multiplicity from global tracks
2236   for(Int_t itracks = 0; itracks < nESDTracks; itracks++) { // flag the tracks
2237     AliESDtrack* track = esd->GetTrack(itracks);
2238
2239     // if track is a secondary from a V0, flag as a secondary
2240     if (track->IsOn(AliESDtrack::kMultInV0)) {
2241       track->SetBit(kSecBit);
2242       continue;
2243     }
2244     /* done via proper DCA cut
2245     //secondary?
2246     if (track->IsOn(AliESDtrack::kMultSec)) {
2247       track->SetBit(kSecBit);
2248       continue;
2249     }
2250     */
2251     // check tracks with ITS part
2252     //*******************************************************************************************************
2253     if (track->IsOn(AliESDtrack::kITSin) && !track->IsOn(AliESDtrack::kITSpureSA) && trackType == kTrackletsITSTPC) { // track has ITS part but is not an ITS_SA
2254       //*******************************************************************************************************
2255       // TPC+ITS
2256       if (track->IsOn(AliESDtrack::kTPCin)) {  // Global track, has ITS and TPC contributions
2257           if (fgMultEstTrackCuts[kMultEstTrackCutGlobal]->AcceptTrack(track)) { // good ITSTPC track
2258             if (fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->AcceptTrack(track) || fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->AcceptTrack(track)) {
2259               tracksITSTPC++; //global track counted
2260               globalBits.SetBitNumber(itracks);
2261             }
2262             else track->SetBit(kSecBit); // large DCA -> secondary, don't count either track not associated tracklet
2263           }
2264           else track->SetBit(kRejBit); // bad quality, don't count the track, but may count tracklet if associated
2265       }
2266       //*******************************************************************************************************
2267       // ITS complementary
2268       else if (fgMultEstTrackCuts[kMultEstTrackCutITSSA]->AcceptTrack(track)) { // good ITS complementary track
2269         if (fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->AcceptTrack(track) || fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->AcceptTrack(track)) {
2270           tracksITSTPCSA_complementary++;
2271           globalBits.SetBitNumber(itracks);
2272         }
2273         else track->SetBit(kSecBit); // large DCA -> secondary, don't count either track not associated tracklet
2274       }
2275       else track->SetBit(kRejBit); // bad quality, don't count the track, but may count tracklet if associated
2276     }
2277     //*******************************************************************************************************
2278     // check tracks from ITS_SA_PURE
2279     if (track->IsOn(AliESDtrack::kITSin) && track->IsOn(AliESDtrack::kITSpureSA) && trackType == kTrackletsITSSA){
2280       if (fgMultEstTrackCuts[kMultEstTrackCutITSSA]->AcceptTrack(track)) { // good ITSSA track
2281         if (fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->AcceptTrack(track) || fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->AcceptTrack(track)) {
2282           tracksITSSA++;
2283           pureITSBits.SetBitNumber(itracks);
2284         }
2285         else track->SetBit(kRejBit);
2286       }
2287       else track->SetBit(kRejBit);
2288     }
2289   }//ESD tracks counted
2290
2291   //*******************************************************************************************************
2292   // get multiplicity from ITS tracklets to complement TPC+ITS, and ITSpureSA
2293   const AliMultiplicity* spdmult = esd->GetMultiplicity();    // spd multiplicity object
2294   for (Int_t i=0; i<spdmult->GetNumberOfTracklets(); ++i) {
2295     if (TMath::Abs(spdmult->GetEta(i)) > etaRange) continue; // eta selection for tracklets
2296
2297     // if counting tracks+tracklets, check if clusters were already used in tracks
2298     Int_t id1,id2,id3,id4;
2299     spdmult->GetTrackletTrackIDs(i,0,id1,id2); // references for eventual Global/ITS_SA tracks
2300     AliESDtrack* tr1 = id1>=0 ? esd->GetTrack(id1) : 0;
2301     spdmult->GetTrackletTrackIDs(i,1,id3,id4); // references for eventual ITS_SA_pure tracks
2302     AliESDtrack* tr3 = id3>=0 ? esd->GetTrack(id3) : 0;
2303
2304     // are both clusters from the same tracks? If not, skip the tracklet (shouldn't change things much)
2305     if ((id1!=id2 && id1>=0 && id2>=0) || (id3!=id4 && id3>=0 && id4>=0)) continue;
2306
2307     Bool_t bUsedInGlobal = (id1 != -1) ? globalBits.TestBitNumber(id1) : 0;// has associated global track been associated to a previous tracklet?
2308     Bool_t bUsedInPureITS = (id3 != -1) ? pureITSBits.TestBitNumber(id3) : 0;// has associated pure ITS track been associated to a previous tracklet?
2309     //*******************************************************************************************************
2310     if (trackType == kTrackletsITSTPC) {
2311       // count tracklets towards global+complementary tracks
2312       if (  (tr1 && !tr1->TestBit(kSecBit)) &&   // reject as secondary
2313         (tr1 &&  tr1->TestBit(kRejBit)) ) {  // count tracklet as bad quality track
2314           if(!bUsedInGlobal){
2315             ++trackletsITSTPC_complementary;
2316             if(id1>0) globalBits.SetBitNumber(id1); // mark global track linked to this tracklet as "associated"
2317           }
2318       }
2319       else if(id1<0) {
2320         ++trackletsITSTPC_complementary; // if no associated track, count the tracklet
2321       }
2322     } else {
2323       // count tracklets towards ITS_SA_pure tracks
2324       if (  (tr3 && !tr3->TestBit(kSecBit)) &&   // reject as secondary
2325         (tr3 &&  tr3->TestBit(kRejBit)) ) { // count tracklet as bad quality track
2326         if(!bUsedInPureITS) {
2327           ++trackletsITSSA_complementary;
2328           if(id3>0) pureITSBits.SetBitNumber(id3); // mark global track linked to this tracklet as "associated"
2329         }
2330       }
2331       else if(id3<0) {
2332         ++trackletsITSSA_complementary; // if no associated track, count the tracklet
2333       }
2334     }
2335   }
2336
2337   //*******************************************************************************************************
2338   if (trackType == kTrackletsITSTPC)
2339     multiplicityEstimate = tracksITSTPC + tracksITSTPCSA_complementary + trackletsITSTPC_complementary;
2340   else
2341     multiplicityEstimate = tracksITSSA + trackletsITSSA_complementary;
2342
2343   return multiplicityEstimate;
2344 }
2345
2346 //____________________________________________________________________
2347 void AliESDtrackCuts::SetRequireStandardTOFmatchCuts(){
2348
2349         // setting the TOF cuts flags (kTOFout = TOF matching distance) to true, to include the selection on the standard TOF matching
2350
2351         SetRequireTOFout(kTRUE);
2352         SetFlagCutTOFdistance(kTRUE);
2353         SetCutTOFdistance(3.);
2354
2355 }
2356