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