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   }
 734   esdTrackCuts->SetMaxDCAToVertexZ(1.e6);
 735   esdTrackCuts->SetDCAToVertex2D(kFALSE);
 736   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 737   //esdTrackCuts->SetEtaRange(-0.8,+0.8);
 738
 739   esdTrackCuts->SetMaxChi2PerClusterITS(36);
 740   esdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
 741
 742   return esdTrackCuts;
 743 }
 744
 745 //____________________________________________________________________
 746 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(Bool_t selPrimaries,Int_t clusterCut)
 747 {
 748   // creates an AliESDtrackCuts object and fills it with standard values for ITS-TPC cuts for pp 2010 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 (2010 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.0026+0.0050/pt^1.01)
 776     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
 777   }
 778   esdTrackCuts->SetMaxDCAToVertexZ(2);
 779   esdTrackCuts->SetDCAToVertex2D(kFALSE);
 780   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 781
 782   esdTrackCuts->SetMaxChi2PerClusterITS(36);
 783   esdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
 784
 785   return esdTrackCuts;
 786 }
 787
 788 //____________________________________________________________________
 789 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSPureSATrackCuts2009(Bool_t selPrimaries, Bool_t useForPid)
 790 {
 791   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks
 792
 793   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 794   esdTrackCuts->SetRequireITSPureStandAlone(kTRUE);
 795   esdTrackCuts->SetRequireITSRefit(kTRUE);
 796   esdTrackCuts->SetMinNClustersITS(4);
 797   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 798                                          AliESDtrackCuts::kAny);
 799   esdTrackCuts->SetMaxChi2PerClusterITS(1.);
 800
 801   if(selPrimaries) {
 802     // 7*(0.0085+0.0026/pt^1.55)
 803     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0595+0.0182/pt^1.55");
 804   }
 805   if(useForPid){
 806     esdTrackCuts->SetRequireITSPid(kTRUE);
 807   }
 808   return esdTrackCuts;
 809 }
 810
 811 //____________________________________________________________________
 812 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSPureSATrackCuts2010(Bool_t selPrimaries, Bool_t useForPid)
 813 {
 814   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks - pp 2010
 815
 816   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 817   esdTrackCuts->SetRequireITSPureStandAlone(kTRUE);
 818   esdTrackCuts->SetRequireITSRefit(kTRUE);
 819   esdTrackCuts->SetMinNClustersITS(4);
 820   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 821                                          AliESDtrackCuts::kAny);
 822   esdTrackCuts->SetMaxChi2PerClusterITS(2.5);
 823
 824   if(selPrimaries) {
 825     // 7*(0.0033+0.0045/pt^1.3)
 826     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0231+0.0315/pt^1.3");
 827   }
 828   if(useForPid){
 829     esdTrackCuts->SetRequireITSPid(kTRUE);
 830   }
 831   return esdTrackCuts;
 832 }
 833
 834 //____________________________________________________________________
 835 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSSATrackCuts2009(Bool_t selPrimaries, Bool_t useForPid)
 836 {
 837   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks
 838
 839   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 840   esdTrackCuts->SetRequireITSStandAlone(kTRUE);
 841   esdTrackCuts->SetRequireITSPureStandAlone(kFALSE);
 842   esdTrackCuts->SetRequireITSRefit(kTRUE);
 843   esdTrackCuts->SetMinNClustersITS(4);
 844   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 845                                          AliESDtrackCuts::kAny);
 846   esdTrackCuts->SetMaxChi2PerClusterITS(1.);
 847
 848   if(selPrimaries) {
 849     // 7*(0.0085+0.0026/pt^1.55)
 850     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0595+0.0182/pt^1.55");
 851   }
 852   if(useForPid){
 853     esdTrackCuts->SetRequireITSPid(kTRUE);
 854   }
 855   return esdTrackCuts;
 856 }
 857
 858 //____________________________________________________________________
 859 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSSATrackCuts2010(Bool_t selPrimaries, Bool_t useForPid)
 860 {
 861   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks --pp 2010
 862
 863   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 864   esdTrackCuts->SetRequireITSStandAlone(kTRUE);
 865   esdTrackCuts->SetRequireITSPureStandAlone(kFALSE);
 866   esdTrackCuts->SetRequireITSRefit(kTRUE);
 867   esdTrackCuts->SetMinNClustersITS(4);
 868   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 869                                          AliESDtrackCuts::kAny);
 870   esdTrackCuts->SetMaxChi2PerClusterITS(2.5);
 871
 872   if(selPrimaries) {
 873     // 7*(0.0033+0.0045/pt^1.3)
 874     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0231+0.0315/pt^1.3");
 875   }
 876   if(useForPid){
 877     esdTrackCuts->SetRequireITSPid(kTRUE);
 878   }
 879   return esdTrackCuts;
 880 }
 881
 882 //____________________________________________________________________
 883 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSSATrackCutsPbPb2010(Bool_t selPrimaries, Bool_t useForPid)
 884 {
 885   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks -- PbPb 2010
 886
 887   AliESDtrackCuts* esdTrackCuts = GetStandardITSSATrackCuts2010(selPrimaries, useForPid);
 888   esdTrackCuts->SetMaxNOfMissingITSPoints(1);
 889
 890   return esdTrackCuts;
 891 }
 892 //____________________________________________________________________
 893
 894 AliESDtrackCuts* AliESDtrackCuts::GetStandardV0DaughterCuts()
 895 {
 896   // creates a AliESDtrackCuts object and fills it with standard cuts for V0 daughters
 897   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 898   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 899   esdTrackCuts->SetMinNClustersTPC(70);
 900   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 901   return esdTrackCuts;
 902 }
 903
 904 //____________________________________________________________________
 905 Int_t AliESDtrackCuts::GetReferenceMultiplicity(const AliESDEvent* esd, Bool_t tpcOnly)
 906 {
 907   // Gets reference multiplicity following the standard cuts and a defined fiducial volume
 908   // tpcOnly = kTRUE -> consider TPC-only tracks
 909   //         = kFALSE -> consider global tracks
 910   //
 911   // DEPRECATED Use GetReferenceMultiplicity with the enum as second argument instead
 912
 913   if (!tpcOnly)
 914   {
 915     AliErrorClass("Not implemented for global tracks!");
 916     return -1;
 917   }
 918
 919   static AliESDtrackCuts* esdTrackCuts = 0;
 920   if (!esdTrackCuts)
 921   {
 922     esdTrackCuts = GetStandardTPCOnlyTrackCuts();
 923     esdTrackCuts->SetEtaRange(-0.8, 0.8);
 924     esdTrackCuts->SetPtRange(0.15);
 925   }
 926
 927   Int_t nTracks = esdTrackCuts->CountAcceptedTracks(esd);
 928
 929   return nTracks;
 930 }
 931
 932 //____________________________________________________________________
 933 Float_t AliESDtrackCuts::GetSigmaToVertex(const AliESDtrack* const esdTrack)
 934 {
 935   // Calculates the number of sigma to the vertex.
 936
 937   Float_t b[2];
 938   Float_t bRes[2];
 939   Float_t bCov[3];
 940   esdTrack->GetImpactParameters(b,bCov);
 941
 942   if (bCov[0]<=0 || bCov[2]<=0) {
 943     AliDebugClass(1, "Estimated b resolution lower or equal zero!");
 944     bCov[0]=0; bCov[2]=0;
 945   }
 946   bRes[0] = TMath::Sqrt(bCov[0]);
 947   bRes[1] = TMath::Sqrt(bCov[2]);
 948
 949   // -----------------------------------
 950   // How to get to a n-sigma cut?
 951   //
 952   // The accumulated statistics from 0 to d is
 953   //
 954   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 955   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 956   //
 957   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
 958   // Can this be expressed in a different way?
 959
 960   if (bRes[0] == 0 || bRes[1] ==0)
 961     return -1;
 962
 963   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 964
 965   // work around precision problem
 966   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 967   // 1e-15 corresponds to nsigma ~ 7.7
 968   if (TMath::Exp(-d * d / 2) < 1e-15)
 969     return 1000;
 970
 971   Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 972   return nSigma;
 973 }
 974
 975 void AliESDtrackCuts::EnableNeededBranches(TTree* tree)
 976 {
 977   // enables the branches needed by AcceptTrack, for a list see comment of AcceptTrack
 978
 979   tree->SetBranchStatus("fTracks.fFlags", 1);
 980   tree->SetBranchStatus("fTracks.fITSncls", 1);
 981   tree->SetBranchStatus("fTracks.fTPCncls", 1);
 982   tree->SetBranchStatus("fTracks.fITSchi2", 1);
 983   tree->SetBranchStatus("fTracks.fTPCchi2", 1);
 984   tree->SetBranchStatus("fTracks.fC*", 1);
 985   tree->SetBranchStatus("fTracks.fD", 1);
 986   tree->SetBranchStatus("fTracks.fZ", 1);
 987   tree->SetBranchStatus("fTracks.fCdd", 1);
 988   tree->SetBranchStatus("fTracks.fCdz", 1);
 989   tree->SetBranchStatus("fTracks.fCzz", 1);
 990   tree->SetBranchStatus("fTracks.fP*", 1);
 991   tree->SetBranchStatus("fTracks.fR*", 1);
 992   tree->SetBranchStatus("fTracks.fKinkIndexes*", 1);
 993 }
 994
 995 //____________________________________________________________________
 996 Bool_t AliESDtrackCuts::AcceptTrack(const AliESDtrack* esdTrack, const AliESDEvent* esdEvent)
 997 {
 998   //
 999   // figure out if the tracks survives all the track cuts defined
1000   //
1001   // the different quality parameter and kinematic values are first
1002   // retrieved from the track. then it is found out what cuts the
1003   // track did not survive and finally the cuts are imposed.
1004
1005   // this function needs the following branches:
1006   // fTracks.fFlags
1007   // fTracks.fITSncls
1008   // fTracks.fTPCncls
1009   // fTracks.fITSchi2
1010   // fTracks.fTPCchi2
1011   // fTracks.fC   //GetExternalCovariance
1012   // fTracks.fD   //GetImpactParameters
1013   // fTracks.fZ   //GetImpactParameters
1014   // fTracks.fCdd //GetImpactParameters
1015   // fTracks.fCdz //GetImpactParameters
1016   // fTracks.fCzz //GetImpactParameters
1017   // fTracks.fP   //GetPxPyPz
1018   // fTracks.fR   //GetMass
1019   // fTracks.fP   //GetMass
1020   // fTracks.fKinkIndexes
1021   //
1022   // esdEvent is only required for the MaxChi2TPCConstrainedVsGlobal
1023
1024   UInt_t status = esdTrack->GetStatus();
1025
1026   // getting quality parameters from the ESD track
1027   Int_t nClustersITS = esdTrack->GetITSclusters(0);
1028   Int_t nClustersTPC = -1;
1029   if(fCutRequireTPCStandAlone) {
1030     nClustersTPC = esdTrack->GetTPCNclsIter1();
1031   }
1032   else {
1033     nClustersTPC = esdTrack->GetTPCclusters(0);
1034   }
1035
1036   //Pt dependent NClusters Cut
1037   if(f1CutMinNClustersTPCPtDep) {
1038     if(esdTrack->Pt()<fCutMaxPtDepNClustersTPC)
1039       fCutMinNClusterTPC = f1CutMinNClustersTPCPtDep->Eval(esdTrack->Pt());
1040     else
1041       fCutMinNClusterTPC = f1CutMinNClustersTPCPtDep->Eval(fCutMaxPtDepNClustersTPC);
1042   }
1043
1044   Float_t nCrossedRowsTPC = esdTrack->GetTPCCrossedRows();
1045   Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1046   if (esdTrack->GetTPCNclsF()>0) {
1047     ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC / esdTrack->GetTPCNclsF();
1048   }
1049
1050   Int_t nClustersTPCShared = esdTrack->GetTPCnclsS();
1051   Float_t fracClustersTPCShared = -1.;
1052
1053   Float_t chi2PerClusterITS = -1;
1054   Float_t chi2PerClusterTPC = -1;
1055   if (nClustersITS!=0)
1056     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1057   if (nClustersTPC!=0) {
1058     if(fCutRequireTPCStandAlone) {
1059       chi2PerClusterTPC = esdTrack->GetTPCchi2Iter1()/Float_t(nClustersTPC);
1060     } else {
1061       chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
1062     }
1063     fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClustersTPC);
1064   }
1065
1066   Double_t extCov[15];
1067   esdTrack->GetExternalCovariance(extCov);
1068
1069   Float_t b[2];
1070   Float_t bCov[3];
1071   esdTrack->GetImpactParameters(b,bCov);
1072   if (bCov[0]<=0 || bCov[2]<=0) {
1073     AliDebug(1, "Estimated b resolution lower or equal zero!");
1074     bCov[0]=0; bCov[2]=0;
1075   }
1076
1077
1078   // set pt-dependent DCA cuts, if requested
1079   SetPtDepDCACuts(esdTrack->Pt());
1080
1081
1082   Float_t dcaToVertexXY = b[0];
1083   Float_t dcaToVertexZ = b[1];
1084
1085   Float_t dcaToVertex = -1;
1086
1087   if (fCutDCAToVertex2D)
1088   {
1089     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1090   }
1091   else
1092     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1093
1094   // getting the kinematic variables of the track
1095   // (assuming the mass is known)
1096   Double_t p[3];
1097   esdTrack->GetPxPyPz(p);
1098
1099   Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
1100   Float_t pt       = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2));
1101   Float_t energy   = TMath::Sqrt(TMath::Power(esdTrack->GetMass(),2) + TMath::Power(momentum,2));
1102
1103   //y-eta related calculations
1104   Float_t eta = -100.;
1105   Float_t y   = -100.;
1106   if((momentum != TMath::Abs(p[2]))&&(momentum != 0))
1107     eta = 0.5*TMath::Log((momentum + p[2])/(momentum - p[2]));
1108   if((energy != TMath::Abs(p[2]))&&(momentum != 0))
1109     y = 0.5*TMath::Log((energy + p[2])/(energy - p[2]));
1110
1111   if (extCov[14] < 0)
1112   {
1113     AliWarning(Form("GetSigma1Pt2() returns negative value for external covariance matrix element fC[14]: %f. Corrupted track information, track will not be accepted!", extCov[14]));
1114     return kFALSE;
1115   }
1116   Float_t relUncertainty1Pt = TMath::Sqrt(extCov[14])*pt;
1117
1118   //########################################################################
1119   // cut the track?
1120
1121   Bool_t cuts[kNCuts];
1122   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1123
1124   // track quality cuts
1125   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1126     cuts[0]=kTRUE;
1127   if (fCutRequireTPCStandAlone && (status&AliESDtrack::kTPCin)==0)
1128     cuts[1]=kTRUE;
1129   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1130     cuts[2]=kTRUE;
1131   if (nClustersTPC<fCutMinNClusterTPC)
1132     cuts[3]=kTRUE;
1133   if (nClustersITS<fCutMinNClusterITS)
1134     cuts[4]=kTRUE;
1135   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1136     cuts[5]=kTRUE;
1137   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1138     cuts[6]=kTRUE;
1139   if (extCov[0]  > fCutMaxC11)
1140     cuts[7]=kTRUE;
1141   if (extCov[2]  > fCutMaxC22)
1142     cuts[8]=kTRUE;
1143   if (extCov[5]  > fCutMaxC33)
1144     cuts[9]=kTRUE;
1145   if (extCov[9]  > fCutMaxC44)
1146     cuts[10]=kTRUE;
1147   if (extCov[14]  > fCutMaxC55)
1148     cuts[11]=kTRUE;
1149
1150   // cut 12 and 13 see below
1151
1152   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1153     cuts[14]=kTRUE;
1154   // track kinematics cut
1155   if((momentum < fPMin) || (momentum > fPMax))
1156     cuts[15]=kTRUE;
1157   if((pt < fPtMin) || (pt > fPtMax))
1158     cuts[16] = kTRUE;
1159   if((p[0] < fPxMin) || (p[0] > fPxMax))
1160     cuts[17] = kTRUE;
1161   if((p[1] < fPyMin) || (p[1] > fPyMax))
1162     cuts[18] = kTRUE;
1163   if((p[2] < fPzMin) || (p[2] > fPzMax))
1164     cuts[19] = kTRUE;
1165   if((eta < fEtaMin) || (eta > fEtaMax))
1166     cuts[20] = kTRUE;
1167   if((y < fRapMin) || (y > fRapMax))
1168     cuts[21] = kTRUE;
1169   if (fCutDCAToVertex2D && dcaToVertex > 1)
1170     cuts[22] = kTRUE;
1171   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1172     cuts[23] = kTRUE;
1173   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1174     cuts[24] = kTRUE;
1175   if (fCutDCAToVertex2D && fCutMinDCAToVertexXY > 0 && fCutMinDCAToVertexZ > 0 && dcaToVertexXY*dcaToVertexXY/fCutMinDCAToVertexXY/fCutMinDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMinDCAToVertexZ/fCutMinDCAToVertexZ < 1)
1176     cuts[25] = kTRUE;
1177   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) < fCutMinDCAToVertexXY)
1178     cuts[26] = kTRUE;
1179   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) < fCutMinDCAToVertexZ)
1180     cuts[27] = kTRUE;
1181
1182   for (Int_t i = 0; i < 3; i++)
1183     cuts[28+i] = !CheckITSClusterRequirement(fCutClusterRequirementITS[i], esdTrack->HasPointOnITSLayer(i*2), esdTrack->HasPointOnITSLayer(i*2+1));
1184
1185   if(fCutRequireITSStandAlone || fCutRequireITSpureSA){
1186     if ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin)){
1187       // TPC tracks
1188       cuts[31] = kTRUE;
1189     }else{
1190       // ITS standalone tracks
1191       if(fCutRequireITSStandAlone && !fCutRequireITSpureSA){
1192         if(status & AliESDtrack::kITSpureSA) cuts[31] = kTRUE;
1193       }else if(fCutRequireITSpureSA){
1194         if(!(status & AliESDtrack::kITSpureSA)) cuts[31] = kTRUE;
1195       }
1196     }
1197   }
1198
1199   if (relUncertainty1Pt > fCutMaxRel1PtUncertainty)
1200      cuts[32] = kTRUE;
1201
1202   if (!fCutAcceptSharedTPCClusters && nClustersTPCShared!=0)
1203     cuts[33] = kTRUE;
1204
1205   if (fracClustersTPCShared > fCutMaxFractionSharedTPCClusters)
1206     cuts[34] = kTRUE;
1207
1208   Int_t nITSPointsForPid=0;
1209   UChar_t clumap=esdTrack->GetITSClusterMap();
1210   for(Int_t i=2; i<6; i++){
1211     if(clumap&(1<<i)) ++nITSPointsForPid;
1212   }
1213   if(fCutRequireITSPid && nITSPointsForPid<3) cuts[35] = kTRUE;
1214
1215
1216   if (nCrossedRowsTPC<fCutMinNCrossedRowsTPC)
1217     cuts[36]=kTRUE;
1218   if (ratioCrossedRowsOverFindableClustersTPC<fCutMinRatioCrossedRowsOverFindableClustersTPC)
1219     cuts[37]=kTRUE;
1220
1221   Int_t nMissITSpts=0;
1222   Int_t idet,statusLay;
1223   Float_t xloc,zloc;
1224   for(Int_t iLay=0; iLay<6; iLay++){
1225     Bool_t retc=esdTrack->GetITSModuleIndexInfo(iLay,idet,statusLay,xloc,zloc);
1226     if(retc && statusLay==5) ++nMissITSpts;
1227   }
1228   if(nMissITSpts>fCutMaxMissingITSPoints) cuts[38] = kTRUE;
1229
1230   Bool_t cut=kFALSE;
1231   for (Int_t i=0; i<kNCuts; i++)
1232     if (cuts[i]) {cut = kTRUE;}
1233
1234   // for performance evaluate the CPU intensive cuts only when the others have passed, and when they are requested
1235   Double_t chi2TPCConstrainedVsGlobal = -2;
1236   Float_t nSigmaToVertex = -2;
1237   if (!cut)
1238   {
1239     // getting the track to vertex parameters
1240     if (fCutSigmaToVertexRequired)
1241     {
1242       nSigmaToVertex = GetSigmaToVertex(esdTrack);
1243       if (nSigmaToVertex > fCutNsigmaToVertex && fCutSigmaToVertexRequired)
1244       {
1245         cuts[12] = kTRUE;
1246         cut = kTRUE;
1247       }
1248       // if n sigma could not be calculated
1249       if (nSigmaToVertex<0 && fCutSigmaToVertexRequired)
1250       {
1251         cuts[13] = kTRUE;
1252         cut = kTRUE;
1253       }
1254     }
1255
1256     // max chi2 TPC constrained vs global track only if track passed the other cut
1257     if (fCutMaxChi2TPCConstrainedVsGlobal < 1e9)
1258     {
1259       if (!esdEvent)
1260         AliFatal("fCutMaxChi2TPCConstrainedVsGlobal set but ESD event not provided.");
1261
1262       // get vertex
1263       const AliESDVertex* vertex = 0;
1264       if (fCutMaxChi2TPCConstrainedVsGlobalVertexType & kVertexTracks)
1265         vertex = esdEvent->GetPrimaryVertexTracks();
1266
1267       if ((!vertex || !vertex->GetStatus()) && fCutMaxChi2TPCConstrainedVsGlobalVertexType & kVertexSPD)
1268         vertex = esdEvent->GetPrimaryVertexSPD();
1269
1270       if ((!vertex || !vertex->GetStatus()) && fCutMaxChi2TPCConstrainedVsGlobalVertexType & kVertexTPC)
1271         vertex = esdEvent->GetPrimaryVertexTPC();
1272
1273       if (vertex->GetStatus())
1274         chi2TPCConstrainedVsGlobal = esdTrack->GetChi2TPCConstrainedVsGlobal(vertex);
1275
1276       if (chi2TPCConstrainedVsGlobal < 0 || chi2TPCConstrainedVsGlobal > fCutMaxChi2TPCConstrainedVsGlobal)
1277       {
1278         cuts[39] = kTRUE;
1279         cut = kTRUE;
1280       }
1281     }
1282   }
1283
1284   //########################################################################
1285   // filling histograms
1286   if (fHistogramsOn) {
1287     fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n tracks")));
1288     if (cut)
1289       fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut tracks")));
1290
1291     for (Int_t i=0; i<kNCuts; i++) {
1292       if (fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i]) < 1)
1293         AliFatal(Form("Inconsistency! Cut %d with name %s not found", i, fgkCutNames[i]));
1294
1295       if (cuts[i])
1296         fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
1297
1298       for (Int_t j=i; j<kNCuts; j++) {
1299         if (cuts[i] && cuts[j]) {
1300           Float_t xC = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
1301           Float_t yC = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
1302           fhCutCorrelation->Fill(xC, yC);
1303         }
1304       }
1305     }
1306   }
1307
1308   // now we loop over the filling of the histograms twice: once "before" the cut, once "after"
1309   // the code is not in a function due to too many local variables that would need to be passed
1310
1311   for (Int_t id = 0; id < 2; id++)
1312   {
1313     // id = 0 --> before cut
1314     // id = 1 --> after cut
1315
1316     if (fHistogramsOn)
1317     {
1318       fhNClustersITS[id]->Fill(nClustersITS);
1319       fhNClustersTPC[id]->Fill(nClustersTPC);
1320       fhNSharedClustersTPC[id]->Fill(nClustersTPCShared);
1321       fhNCrossedRowsTPC[id]->Fill(nCrossedRowsTPC);
1322       fhRatioCrossedRowsOverFindableClustersTPC[id]->Fill(ratioCrossedRowsOverFindableClustersTPC);
1323       fhChi2PerClusterITS[id]->Fill(chi2PerClusterITS);
1324       fhChi2PerClusterTPC[id]->Fill(chi2PerClusterTPC);
1325       fhChi2TPCConstrainedVsGlobal[id]->Fill(chi2TPCConstrainedVsGlobal);
1326       fhNClustersForITSPID[id]->Fill(nITSPointsForPid);
1327       fhNMissingITSPoints[id]->Fill(nMissITSpts);
1328
1329       fhC11[id]->Fill(extCov[0]);
1330       fhC22[id]->Fill(extCov[2]);
1331       fhC33[id]->Fill(extCov[5]);
1332       fhC44[id]->Fill(extCov[9]);
1333       fhC55[id]->Fill(extCov[14]);
1334
1335       fhRel1PtUncertainty[id]->Fill(relUncertainty1Pt);
1336
1337       fhPt[id]->Fill(pt);
1338       fhEta[id]->Fill(eta);
1339
1340       Float_t bRes[2];
1341       bRes[0] = TMath::Sqrt(bCov[0]);
1342       bRes[1] = TMath::Sqrt(bCov[2]);
1343
1344       fhDZ[id]->Fill(b[1]);
1345       fhDXY[id]->Fill(b[0]);
1346       fhDXYDZ[id]->Fill(dcaToVertex);
1347       fhDXYvsDZ[id]->Fill(b[1],b[0]);
1348
1349       if (bRes[0]!=0 && bRes[1]!=0) {
1350         fhDZNormalized[id]->Fill(b[1]/bRes[1]);
1351         fhDXYNormalized[id]->Fill(b[0]/bRes[0]);
1352         fhDXYvsDZNormalized[id]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
1353         fhNSigmaToVertex[id]->Fill(nSigmaToVertex);
1354       }
1355     }
1356
1357     // cut the track
1358     if (cut)
1359       return kFALSE;
1360   }
1361
1362   return kTRUE;
1363 }
1364
1365 //____________________________________________________________________
1366 Bool_t AliESDtrackCuts::CheckITSClusterRequirement(ITSClusterRequirement req, Bool_t clusterL1, Bool_t clusterL2)
1367 {
1368   // checks if the cluster requirement is fullfilled (in this case: return kTRUE)
1369
1370   switch (req)
1371   {
1372         case kOff:        return kTRUE;
1373         case kNone:       return !clusterL1 && !clusterL2;
1374         case kAny:        return clusterL1 || clusterL2;
1375         case kFirst:      return clusterL1;
1376         case kOnlyFirst:  return clusterL1 && !clusterL2;
1377         case kSecond:     return clusterL2;
1378         case kOnlySecond: return clusterL2 && !clusterL1;
1379         case kBoth:       return clusterL1 && clusterL2;
1380   }
1381
1382   return kFALSE;
1383 }
1384
1385 //____________________________________________________________________
1386 AliESDtrack* AliESDtrackCuts::GetTPCOnlyTrack(const AliESDEvent* esd, Int_t iTrack)
1387 {
1388   // Utility function to create a TPC only track from the given esd track
1389   //
1390   // IMPORTANT: The track has to be deleted by the user
1391   //
1392   // NB. most of the functionality to get a TPC only track from an ESD track is in AliESDtrack, where it should be
1393   // there are only missing propagations here that are needed for old data
1394   // this function will therefore become obsolete
1395   //
1396   // adapted from code provided by CKB
1397
1398   if (!esd->GetPrimaryVertexTPC())
1399     return 0; // No TPC vertex no TPC tracks
1400
1401   if(!esd->GetPrimaryVertexTPC()->GetStatus())
1402     return 0; // TPC Vertex is created by default in AliESDEvent, do not use in this case
1403
1404   AliESDtrack* track = esd->GetTrack(iTrack);
1405   if (!track)
1406     return 0;
1407
1408   AliESDtrack *tpcTrack = new AliESDtrack();
1409
1410   // only true if we have a tpc track
1411   if (!track->FillTPCOnlyTrack(*tpcTrack))
1412   {
1413     delete tpcTrack;
1414     return 0;
1415   }
1416
1417   return tpcTrack;
1418 }
1419
1420 //____________________________________________________________________
1421 TObjArray* AliESDtrackCuts::GetAcceptedTracks(const AliESDEvent* esd, Bool_t bTPC)
1422 {
1423   //
1424   // returns an array of all tracks that pass the cuts
1425   // or an array of TPC only tracks (propagated to the TPC vertex during reco)
1426   // tracks that pass the cut
1427   //
1428   // NOTE: List has to be deleted by the user
1429
1430   TObjArray* acceptedTracks = new TObjArray();
1431
1432   // loop over esd tracks
1433   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
1434     if(bTPC){
1435       if(!esd->GetPrimaryVertexTPC())return acceptedTracks; // No TPC vertex no TPC tracks
1436       if(!esd->GetPrimaryVertexTPC()->GetStatus())return acceptedTracks; // No proper TPC vertex, only the default
1437
1438       AliESDtrack *tpcTrack = GetTPCOnlyTrack(esd, iTrack);
1439       if (!tpcTrack)
1440         continue;
1441
1442       if (AcceptTrack(tpcTrack, esd)) {
1443         acceptedTracks->Add(tpcTrack);
1444       }
1445       else
1446         delete tpcTrack;
1447     }
1448     else
1449     {
1450       AliESDtrack* track = esd->GetTrack(iTrack);
1451       if(AcceptTrack(track, esd))
1452         acceptedTracks->Add(track);
1453     }
1454   }
1455   if(bTPC)acceptedTracks->SetOwner(kTRUE);
1456   return acceptedTracks;
1457 }
1458
1459 //____________________________________________________________________
1460 Int_t AliESDtrackCuts::CountAcceptedTracks(const AliESDEvent* const esd)
1461 {
1462   //
1463   // returns an the number of tracks that pass the cuts
1464   //
1465
1466   Int_t count = 0;
1467
1468   // loop over esd tracks
1469   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
1470     AliESDtrack* track = esd->GetTrack(iTrack);
1471     if (AcceptTrack(track, esd))
1472       count++;
1473   }
1474
1475   return count;
1476 }
1477
1478 //____________________________________________________________________
1479  void AliESDtrackCuts::DefineHistograms(Int_t color) {
1480    //
1481    // diagnostics histograms are defined
1482    //
1483
1484    fHistogramsOn=kTRUE;
1485
1486    Bool_t oldStatus = TH1::AddDirectoryStatus();
1487    TH1::AddDirectory(kFALSE);
1488
1489    //###################################################################################
1490    // defining histograms
1491
1492    fhCutStatistics = new TH1F("cut_statistics","cut statistics",kNCuts+4,-0.5,kNCuts+3.5);
1493
1494    fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
1495    fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
1496
1497    fhCutCorrelation = new TH2F("cut_correlation","cut correlation",kNCuts,-0.5,kNCuts-0.5,kNCuts,-0.5,kNCuts-0.5);;
1498
1499    for (Int_t i=0; i<kNCuts; i++) {
1500      fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fgkCutNames[i]);
1501      fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fgkCutNames[i]);
1502      fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fgkCutNames[i]);
1503    }
1504
1505   fhCutStatistics  ->SetLineColor(color);
1506   fhCutCorrelation ->SetLineColor(color);
1507   fhCutStatistics  ->SetLineWidth(2);
1508   fhCutCorrelation ->SetLineWidth(2);
1509
1510   for (Int_t i=0; i<2; i++) {
1511     fhNClustersITS[i]        = new TH1F("nClustersITS"    ,"",8,-0.5,7.5);
1512     fhNClustersTPC[i]        = new TH1F("nClustersTPC"     ,"",165,-0.5,164.5);
1513     fhNSharedClustersTPC[i]  = new TH1F("nSharedClustersTPC"     ,"",165,-0.5,164.5);
1514     fhNCrossedRowsTPC[i]     = new TH1F("nCrossedRowsTPC"  ,"",165,-0.5,164.5);
1515     fhRatioCrossedRowsOverFindableClustersTPC[i]     = new TH1F("ratioCrossedRowsOverFindableClustersTPC"  ,"",60,0,1.5);
1516     fhChi2PerClusterITS[i]   = new TH1F("chi2PerClusterITS","",500,0,10);
1517     fhChi2PerClusterTPC[i]   = new TH1F("chi2PerClusterTPC","",500,0,10);
1518     fhChi2TPCConstrainedVsGlobal[i]   = new TH1F("chi2TPCConstrainedVsGlobal","",600,-2,50);
1519     fhNClustersForITSPID[i]  = new TH1F("nPointsForITSpid","",5,-0.5,4.5);
1520     fhNMissingITSPoints[i]   = new TH1F("nMissingITSClusters","",7,-0.5,6.5);
1521
1522     fhC11[i]                 = new TH1F("covMatrixDiagonal11","",2000,0,20);
1523     fhC22[i]                 = new TH1F("covMatrixDiagonal22","",2000,0,20);
1524     fhC33[i]                 = new TH1F("covMatrixDiagonal33","",1000,0,0.1);
1525     fhC44[i]                 = new TH1F("covMatrixDiagonal44","",1000,0,0.1);
1526     fhC55[i]                 = new TH1F("covMatrixDiagonal55","",1000,0,5);
1527
1528     fhRel1PtUncertainty[i]   = new TH1F("rel1PtUncertainty","",1000,0,5);
1529
1530     fhDXY[i]                 = new TH1F("dXY"    ,"",500,-10,10);
1531     fhDZ[i]                  = new TH1F("dZ"     ,"",500,-10,10);
1532     fhDXYDZ[i]               = new TH1F("dXYDZ"  ,"",500,0,10);
1533     fhDXYvsDZ[i]             = new TH2F("dXYvsDZ","",200,-10,10,200,-10,10);
1534
1535     fhDXYNormalized[i]       = new TH1F("dXYNormalized"    ,"",500,-10,10);
1536     fhDZNormalized[i]        = new TH1F("dZNormalized"     ,"",500,-10,10);
1537     fhDXYvsDZNormalized[i]   = new TH2F("dXYvsDZNormalized","",200,-10,10,200,-10,10);
1538
1539     fhNSigmaToVertex[i]      = new TH1F("nSigmaToVertex","",500,0,10);
1540
1541     fhPt[i]                  = new TH1F("pt"     ,"p_{T} distribution;p_{T} (GeV/c)", 800, 0.0, 10.0);
1542     fhEta[i]                 = new TH1F("eta"     ,"#eta distribution;#eta",40,-2.0,2.0);
1543
1544     fhNClustersITS[i]->SetTitle("n ITS clusters");
1545     fhNClustersTPC[i]->SetTitle("n TPC clusters");
1546     fhNSharedClustersTPC[i]->SetTitle("n TPC shared clusters");
1547     fhChi2PerClusterITS[i]->SetTitle("#Chi^{2} per ITS cluster");
1548     fhChi2PerClusterTPC[i]->SetTitle("#Chi^{2} per TPC cluster");
1549     fhChi2TPCConstrainedVsGlobal[i]->SetTitle("#Chi^{2} TPC constrained track vs global track");
1550     fhNClustersForITSPID[i]->SetTitle("n ITS points for PID");
1551     fhNMissingITSPoints[i]->SetTitle("n ITS layers with missing cluster");
1552
1553     fhC11[i]->SetTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
1554     fhC22[i]->SetTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
1555     fhC33[i]->SetTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
1556     fhC44[i]->SetTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
1557     fhC55[i]->SetTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
1558
1559     fhRel1PtUncertainty[i]->SetTitle("rel. uncertainty of 1/p_{T}");
1560
1561     fhDXY[i]->SetXTitle("transverse impact parameter (cm)");
1562     fhDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
1563     fhDXYDZ[i]->SetTitle("absolute impact parameter;sqrt(dXY**2 + dZ**2)  (cm)");
1564     fhDXYvsDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
1565     fhDXYvsDZ[i]->SetYTitle("transverse impact parameter (cm)");
1566
1567     fhDXYNormalized[i]->SetTitle("normalized trans impact par (n#sigma)");
1568     fhDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
1569     fhDXYvsDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
1570     fhDXYvsDZNormalized[i]->SetYTitle("normalized trans impact par (n#sigma)");
1571     fhNSigmaToVertex[i]->SetTitle("n #sigma to vertex");
1572
1573     fhNClustersITS[i]->SetLineColor(color);   fhNClustersITS[i]->SetLineWidth(2);
1574     fhNClustersTPC[i]->SetLineColor(color);   fhNClustersTPC[i]->SetLineWidth(2);
1575     fhNSharedClustersTPC[i]->SetLineColor(color);  fhNSharedClustersTPC[i]->SetLineWidth(2);
1576     fhChi2PerClusterITS[i]->SetLineColor(color);   fhChi2PerClusterITS[i]->SetLineWidth(2);
1577     fhChi2PerClusterTPC[i]->SetLineColor(color);   fhChi2PerClusterTPC[i]->SetLineWidth(2);
1578     fhChi2TPCConstrainedVsGlobal[i]->SetLineColor(color);   fhChi2TPCConstrainedVsGlobal[i]->SetLineWidth(2);
1579     fhNClustersForITSPID[i]->SetLineColor(color);  fhNClustersForITSPID[i]->SetLineWidth(2);
1580     fhNMissingITSPoints[i]->SetLineColor(color);   fhNMissingITSPoints[i]->SetLineWidth(2);
1581
1582     fhC11[i]->SetLineColor(color);   fhC11[i]->SetLineWidth(2);
1583     fhC22[i]->SetLineColor(color);   fhC22[i]->SetLineWidth(2);
1584     fhC33[i]->SetLineColor(color);   fhC33[i]->SetLineWidth(2);
1585     fhC44[i]->SetLineColor(color);   fhC44[i]->SetLineWidth(2);
1586     fhC55[i]->SetLineColor(color);   fhC55[i]->SetLineWidth(2);
1587
1588     fhRel1PtUncertainty[i]->SetLineColor(color); fhRel1PtUncertainty[i]->SetLineWidth(2);
1589
1590     fhDXY[i]->SetLineColor(color);   fhDXY[i]->SetLineWidth(2);
1591     fhDZ[i]->SetLineColor(color);    fhDZ[i]->SetLineWidth(2);
1592     fhDXYDZ[i]->SetLineColor(color); fhDXYDZ[i]->SetLineWidth(2);
1593
1594     fhDXYNormalized[i]->SetLineColor(color);   fhDXYNormalized[i]->SetLineWidth(2);
1595     fhDZNormalized[i]->SetLineColor(color);    fhDZNormalized[i]->SetLineWidth(2);
1596     fhNSigmaToVertex[i]->SetLineColor(color);  fhNSigmaToVertex[i]->SetLineWidth(2);
1597   }
1598
1599   // The number of sigmas to the vertex is per definition gaussian
1600   ffDTheoretical = new TF1("nSigmaToVertexTheoretical","([0]/2.506628274)*exp(-(x**2)/2)",0,50);
1601   ffDTheoretical->SetParameter(0,1);
1602
1603   TH1::AddDirectory(oldStatus);
1604 }
1605
1606 //____________________________________________________________________
1607 Bool_t AliESDtrackCuts::LoadHistograms(const Char_t* dir)
1608 {
1609   //
1610   // loads the histograms from a file
1611   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
1612   //
1613
1614   if (!dir)
1615     dir = GetName();
1616
1617   if (!gDirectory->cd(dir))
1618     return kFALSE;
1619
1620   ffDTheoretical = dynamic_cast<TF1*> (gDirectory->Get("nSigmaToVertexTheory"));
1621
1622   fhCutStatistics = dynamic_cast<TH1F*> (gDirectory->Get("cut_statistics"));
1623   fhCutCorrelation = dynamic_cast<TH2F*> (gDirectory->Get("cut_correlation"));
1624
1625   for (Int_t i=0; i<2; i++) {
1626     if (i==0)
1627     {
1628       gDirectory->cd("before_cuts");
1629     }
1630     else
1631       gDirectory->cd("after_cuts");
1632
1633     fhNClustersITS[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersITS"     ));
1634     fhNClustersTPC[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersTPC"     ));
1635     fhNSharedClustersTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get("nSharedClustersTPC"     ));
1636     fhNCrossedRowsTPC[i]   = dynamic_cast<TH1F*> (gDirectory->Get("nCrossedRowsTPC"  ));
1637     fhRatioCrossedRowsOverFindableClustersTPC[i]   = dynamic_cast<TH1F*> (gDirectory->Get("ratioCrossedRowsOverFindableClustersTPC"  ));
1638     fhChi2PerClusterITS[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterITS"));
1639     fhChi2PerClusterTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterTPC"));
1640     fhChi2TPCConstrainedVsGlobal[i] = dynamic_cast<TH1F*> (gDirectory->Get("fhChi2TPCConstrainedVsGlobal"));
1641     fhNClustersForITSPID[i] = dynamic_cast<TH1F*> (gDirectory->Get("nPointsForITSpid"));
1642     fhNMissingITSPoints[i] = dynamic_cast<TH1F*> (gDirectory->Get("nMissingITSClusters"));
1643
1644     fhC11[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal11"));
1645     fhC22[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal22"));
1646     fhC33[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal33"));
1647     fhC44[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal44"));
1648     fhC55[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal55"));
1649
1650     fhRel1PtUncertainty[i] = dynamic_cast<TH1F*> (gDirectory->Get("rel1PtUncertainty"));
1651
1652     fhDXY[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXY"    ));
1653     fhDZ[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZ"     ));
1654     fhDXYDZ[i] =   dynamic_cast<TH1F*> (gDirectory->Get("dXYDZ"));
1655     fhDXYvsDZ[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZ"));
1656
1657     fhDXYNormalized[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXYNormalized"    ));
1658     fhDZNormalized[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZNormalized"     ));
1659     fhDXYvsDZNormalized[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZNormalized"));
1660     fhNSigmaToVertex[i] = dynamic_cast<TH1F*> (gDirectory->Get("nSigmaToVertex"));
1661
1662     fhPt[i] = dynamic_cast<TH1F*> (gDirectory->Get("pt"));
1663     fhEta[i] = dynamic_cast<TH1F*> (gDirectory->Get("eta"));
1664
1665     gDirectory->cd("../");
1666   }
1667
1668   gDirectory->cd("..");
1669
1670   return kTRUE;
1671 }
1672
1673 //____________________________________________________________________
1674 void AliESDtrackCuts::SaveHistograms(const Char_t* dir) {
1675   //
1676   // saves the histograms in a directory (dir)
1677   //
1678
1679   if (!fHistogramsOn) {
1680     AliDebug(0, "Histograms not on - cannot save histograms!!!");
1681     return;
1682   }
1683
1684   if (!dir)
1685     dir = GetName();
1686
1687   gDirectory->mkdir(dir);
1688   gDirectory->cd(dir);
1689
1690   gDirectory->mkdir("before_cuts");
1691   gDirectory->mkdir("after_cuts");
1692
1693   // a factor of 2 is needed since n sigma is positive
1694   ffDTheoretical->SetParameter(0,2*fhNSigmaToVertex[0]->Integral("width"));
1695   ffDTheoretical->Write("nSigmaToVertexTheory");
1696
1697   fhCutStatistics->Write();
1698   fhCutCorrelation->Write();
1699
1700   for (Int_t i=0; i<2; i++) {
1701     if (i==0)
1702       gDirectory->cd("before_cuts");
1703     else
1704       gDirectory->cd("after_cuts");
1705
1706     fhNClustersITS[i]        ->Write();
1707     fhNClustersTPC[i]        ->Write();
1708     fhNSharedClustersTPC[i]  ->Write();
1709     fhNCrossedRowsTPC[i]     ->Write();
1710     fhRatioCrossedRowsOverFindableClustersTPC[i]     ->Write();
1711     fhChi2PerClusterITS[i]   ->Write();
1712     fhChi2PerClusterTPC[i]   ->Write();
1713     fhChi2TPCConstrainedVsGlobal[i]   ->Write();
1714     fhNClustersForITSPID[i]  ->Write();
1715     fhNMissingITSPoints[i]   ->Write();
1716
1717     fhC11[i]                 ->Write();
1718     fhC22[i]                 ->Write();
1719     fhC33[i]                 ->Write();
1720     fhC44[i]                 ->Write();
1721     fhC55[i]                 ->Write();
1722
1723     fhRel1PtUncertainty[i]   ->Write();
1724
1725     fhDXY[i]                 ->Write();
1726     fhDZ[i]                  ->Write();
1727     fhDXYDZ[i]               ->Write();
1728     fhDXYvsDZ[i]             ->Write();
1729
1730     fhDXYNormalized[i]       ->Write();
1731     fhDZNormalized[i]        ->Write();
1732     fhDXYvsDZNormalized[i]   ->Write();
1733     fhNSigmaToVertex[i]      ->Write();
1734
1735     fhPt[i]                  ->Write();
1736     fhEta[i]                 ->Write();
1737
1738     gDirectory->cd("../");
1739   }
1740
1741   gDirectory->cd("../");
1742 }
1743
1744 //____________________________________________________________________
1745 void AliESDtrackCuts::DrawHistograms()
1746 {
1747   // draws some histograms
1748
1749   TCanvas* canvas1 = new TCanvas(Form("%s_1", GetName()), "Track Quality Results1", 800, 800);
1750   canvas1->Divide(2, 2);
1751
1752   canvas1->cd(1);
1753   fhNClustersTPC[0]->SetStats(kFALSE);
1754   fhNClustersTPC[0]->Draw();
1755
1756   canvas1->cd(2);
1757   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1758   fhChi2PerClusterTPC[0]->Draw();
1759
1760   canvas1->cd(3);
1761   fhNSigmaToVertex[0]->SetStats(kFALSE);
1762   fhNSigmaToVertex[0]->GetXaxis()->SetRangeUser(0, 10);
1763   fhNSigmaToVertex[0]->Draw();
1764
1765   canvas1->SaveAs(Form("%s_%s.gif", GetName(), canvas1->GetName()));
1766
1767   TCanvas* canvas2 = new TCanvas(Form("%s_2", GetName()), "Track Quality Results2", 1200, 800);
1768   canvas2->Divide(3, 2);
1769
1770   canvas2->cd(1);
1771   fhC11[0]->SetStats(kFALSE);
1772   gPad->SetLogy();
1773   fhC11[0]->Draw();
1774
1775   canvas2->cd(2);
1776   fhC22[0]->SetStats(kFALSE);
1777   gPad->SetLogy();
1778   fhC22[0]->Draw();
1779
1780   canvas2->cd(3);
1781   fhC33[0]->SetStats(kFALSE);
1782   gPad->SetLogy();
1783   fhC33[0]->Draw();
1784
1785   canvas2->cd(4);
1786   fhC44[0]->SetStats(kFALSE);
1787   gPad->SetLogy();
1788   fhC44[0]->Draw();
1789
1790   canvas2->cd(5);
1791   fhC55[0]->SetStats(kFALSE);
1792   gPad->SetLogy();
1793   fhC55[0]->Draw();
1794
1795   canvas2->cd(6);
1796   fhRel1PtUncertainty[0]->SetStats(kFALSE);
1797   gPad->SetLogy();
1798   fhRel1PtUncertainty[0]->Draw();
1799
1800   canvas2->SaveAs(Form("%s_%s.gif", GetName(), canvas2->GetName()));
1801
1802   TCanvas* canvas3 = new TCanvas(Form("%s_3", GetName()), "Track Quality Results3", 1200, 800);
1803   canvas3->Divide(3, 2);
1804
1805   canvas3->cd(1);
1806   fhDXY[0]->SetStats(kFALSE);
1807   gPad->SetLogy();
1808   fhDXY[0]->Draw();
1809
1810   canvas3->cd(2);
1811   fhDZ[0]->SetStats(kFALSE);
1812   gPad->SetLogy();
1813   fhDZ[0]->Draw();
1814
1815   canvas3->cd(3);
1816   fhDXYvsDZ[0]->SetStats(kFALSE);
1817   gPad->SetLogz();
1818   gPad->SetRightMargin(0.15);
1819   fhDXYvsDZ[0]->Draw("COLZ");
1820
1821   canvas3->cd(4);
1822   fhDXYNormalized[0]->SetStats(kFALSE);
1823   gPad->SetLogy();
1824   fhDXYNormalized[0]->Draw();
1825
1826   canvas3->cd(5);
1827   fhDZNormalized[0]->SetStats(kFALSE);
1828   gPad->SetLogy();
1829   fhDZNormalized[0]->Draw();
1830
1831   canvas3->cd(6);
1832   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1833   gPad->SetLogz();
1834   gPad->SetRightMargin(0.15);
1835   fhDXYvsDZNormalized[0]->Draw("COLZ");
1836
1837   canvas3->SaveAs(Form("%s_%s.gif", GetName(), canvas3->GetName()));
1838
1839   TCanvas* canvas4 = new TCanvas(Form("%s_4", GetName()), "Track Quality Results4", 800, 500);
1840   canvas4->Divide(2, 1);
1841
1842   canvas4->cd(1);
1843   fhCutStatistics->SetStats(kFALSE);
1844   fhCutStatistics->LabelsOption("v");
1845   gPad->SetBottomMargin(0.3);
1846   fhCutStatistics->Draw();
1847
1848   canvas4->cd(2);
1849   fhCutCorrelation->SetStats(kFALSE);
1850   fhCutCorrelation->LabelsOption("v");
1851   gPad->SetBottomMargin(0.3);
1852   gPad->SetLeftMargin(0.3);
1853   fhCutCorrelation->Draw("COLZ");
1854
1855   canvas4->SaveAs(Form("%s_%s.gif", GetName(), canvas4->GetName()));
1856
1857   /*canvas->cd(1);
1858   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1859   fhDXYvsDZNormalized[0]->DrawCopy("COLZ");
1860
1861   canvas->cd(2);
1862   fhNClustersTPC[0]->SetStats(kFALSE);
1863   fhNClustersTPC[0]->DrawCopy();
1864
1865   canvas->cd(3);
1866   fhChi2PerClusterITS[0]->SetStats(kFALSE);
1867   fhChi2PerClusterITS[0]->DrawCopy();
1868   fhChi2PerClusterITS[1]->SetLineColor(2);
1869   fhChi2PerClusterITS[1]->DrawCopy("SAME");
1870
1871   canvas->cd(4);
1872   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1873   fhChi2PerClusterTPC[0]->DrawCopy();
1874   fhChi2PerClusterTPC[1]->SetLineColor(2);
1875   fhChi2PerClusterTPC[1]->DrawCopy("SAME");*/
1876 }
1877 //--------------------------------------------------------------------------
1878 void AliESDtrackCuts::SetPtDepDCACuts(Double_t pt) {
1879   //
1880   // set the pt-dependent DCA cuts
1881   //
1882
1883   if(f1CutMaxDCAToVertexXYPtDep) {
1884      fCutMaxDCAToVertexXY=f1CutMaxDCAToVertexXYPtDep->Eval(pt);
1885   }
1886
1887   if(f1CutMaxDCAToVertexZPtDep) {
1888     fCutMaxDCAToVertexZ=f1CutMaxDCAToVertexZPtDep->Eval(pt);
1889   }
1890
1891   if(f1CutMinDCAToVertexXYPtDep) {
1892     fCutMinDCAToVertexXY=f1CutMinDCAToVertexXYPtDep->Eval(pt);
1893   }
1894
1895   if(f1CutMinDCAToVertexZPtDep) {
1896     fCutMinDCAToVertexZ=f1CutMinDCAToVertexZPtDep->Eval(pt);
1897   }
1898
1899
1900   return;
1901 }
1902
1903
1904
1905 //--------------------------------------------------------------------------
1906 Bool_t AliESDtrackCuts::CheckPtDepDCA(TString dist,Bool_t print) const {
1907   //
1908   // Check the correctness of the string syntax
1909   //
1910   Bool_t retval=kTRUE;
1911
1912   if(!dist.Contains("pt")) {
1913     if(print) AliError("string must contain \"pt\"");
1914     retval= kFALSE;
1915   }
1916   return retval;
1917 }
1918
1919  void AliESDtrackCuts::SetMaxDCAToVertexXYPtDep(const char *dist){
1920
1921    if(f1CutMaxDCAToVertexXYPtDep){
1922      delete f1CutMaxDCAToVertexXYPtDep;
1923      // resetiing both
1924      f1CutMaxDCAToVertexXYPtDep = 0;
1925      fCutMaxDCAToVertexXYPtDep = "";
1926    }
1927    if(!CheckPtDepDCA(dist,kTRUE)){
1928      return;
1929    }
1930    fCutMaxDCAToVertexXYPtDep = dist;
1931    TString tmp(dist);
1932    tmp.ReplaceAll("pt","x");
1933    f1CutMaxDCAToVertexXYPtDep = new TFormula("f1CutMaxDCAToVertexXYPtDep",tmp.Data());
1934
1935 }
1936
1937  void AliESDtrackCuts::SetMaxDCAToVertexZPtDep(const char *dist){
1938
1939
1940    if(f1CutMaxDCAToVertexZPtDep){
1941      delete f1CutMaxDCAToVertexZPtDep;
1942      // resetiing both
1943      f1CutMaxDCAToVertexZPtDep = 0;
1944      fCutMaxDCAToVertexZPtDep = "";
1945    }
1946    if(!CheckPtDepDCA(dist,kTRUE))return;
1947
1948    fCutMaxDCAToVertexZPtDep = dist;
1949    TString tmp(dist);
1950    tmp.ReplaceAll("pt","x");
1951    f1CutMaxDCAToVertexZPtDep = new TFormula("f1CutMaxDCAToVertexZPtDep",tmp.Data());
1952
1953
1954 }
1955
1956
1957  void AliESDtrackCuts::SetMinDCAToVertexXYPtDep(const char *dist){
1958
1959
1960    if(f1CutMinDCAToVertexXYPtDep){
1961      delete f1CutMinDCAToVertexXYPtDep;
1962      // resetiing both
1963      f1CutMinDCAToVertexXYPtDep = 0;
1964      fCutMinDCAToVertexXYPtDep = "";
1965    }
1966    if(!CheckPtDepDCA(dist,kTRUE))return;
1967
1968    fCutMinDCAToVertexXYPtDep = dist;
1969    TString tmp(dist);
1970    tmp.ReplaceAll("pt","x");
1971    f1CutMinDCAToVertexXYPtDep = new TFormula("f1CutMinDCAToVertexXYPtDep",tmp.Data());
1972
1973 }
1974
1975
1976  void AliESDtrackCuts::SetMinDCAToVertexZPtDep(const char *dist){
1977
1978
1979
1980    if(f1CutMinDCAToVertexZPtDep){
1981      delete f1CutMinDCAToVertexZPtDep;
1982      // resetiing both
1983      f1CutMinDCAToVertexZPtDep = 0;
1984      fCutMinDCAToVertexZPtDep = "";
1985    }
1986    if(!CheckPtDepDCA(dist,kTRUE))return;
1987    fCutMinDCAToVertexZPtDep = dist;
1988    TString tmp(dist);
1989    tmp.ReplaceAll("pt","x");
1990    f1CutMinDCAToVertexZPtDep = new TFormula("f1CutMinDCAToVertexZPtDep",tmp.Data());
1991 }
1992
1993 AliESDtrackCuts* AliESDtrackCuts::GetMultEstTrackCuts(MultEstTrackCuts cut)
1994 {
1995   // returns the multiplicity estimator track cuts objects to allow for user configuration
1996   // upon first call the objects are created
1997   //
1998   // the cut defined here correspond to GetStandardITSTPCTrackCuts2010 (apart from the one for without SPD)
1999
2000   if (!fgMultEstTrackCuts[kMultEstTrackCutGlobal])
2001   {
2002     // quality cut on ITS+TPC tracks
2003     fgMultEstTrackCuts[kMultEstTrackCutGlobal] = new AliESDtrackCuts();
2004     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetMinNClustersTPC(70);
2005     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetMaxChi2PerClusterTPC(4);
2006     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetAcceptKinkDaughters(kFALSE);
2007     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetRequireTPCRefit(kTRUE);
2008     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetRequireITSRefit(kTRUE);
2009     //multiplicity underestimate if we use global tracks with |eta| > 0.9
2010     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetEtaRange(-0.9, 0.9);
2011
2012     // quality cut on ITS_SA tracks (complementary to ITS+TPC)
2013     fgMultEstTrackCuts[kMultEstTrackCutITSSA] = new AliESDtrackCuts();
2014     fgMultEstTrackCuts[kMultEstTrackCutITSSA]->SetRequireITSRefit(kTRUE);
2015
2016     // primary selection for tracks with SPD hits
2017     fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD] = new AliESDtrackCuts();
2018     fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
2019     fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
2020     fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->SetMaxDCAToVertexZ(2);
2021
2022     // primary selection for tracks w/o SPD hits
2023     fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD] = new AliESDtrackCuts();
2024     fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kNone);
2025     fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->SetMaxDCAToVertexXYPtDep("1.5*(0.0182+0.0350/pt^1.01)");
2026     fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->SetMaxDCAToVertexZ(2);
2027   }
2028
2029   return fgMultEstTrackCuts[cut];
2030 }
2031
2032 Int_t AliESDtrackCuts::GetReferenceMultiplicity(const AliESDEvent* esd, MultEstTrackType trackType, Float_t etaRange)
2033 {
2034   // Get multiplicity estimate based on TPC/ITS tracks and tracklets
2035   // Adapted for AliESDtrackCuts from a version developed by: Ruben Shahoyan, Anton Alkin, Arvinder Palaha
2036   //
2037   // Returns a negative value if no reliable estimate can be provided:
2038   //   -1 SPD vertex missing
2039   //   -2 SPD VertexerZ dispersion too large
2040   //   -3 Track vertex missing (not checked for kTracklets)
2041   //   -4 Distance between SPD and track vertices too large (not checked for kTracklets)
2042   //
2043   // WARNING This functions does not cut on the z vtx. Depending on the eta range requested, you need to restrict your z vertex range!
2044   //
2045   // Strategy for combined estimators
2046   //   1. Count global tracks and flag them
2047   //   2. Count ITSSA as complementaries for ITSTPC+ or as main tracks
2048   //   3. Count the complementary SPD tracklets
2049
2050   const AliESDVertex* vertices[2];
2051   vertices[0] = esd->GetPrimaryVertexSPD();
2052   vertices[1] = esd->GetPrimaryVertexTracks();
2053
2054   if (!vertices[0]->GetStatus())
2055   {
2056     AliDebugClass(AliLog::kDebug, "No SPD vertex. Not able to make a reliable multiplicity estimate.");
2057     return -1;
2058   }
2059
2060   if (vertices[0]->IsFromVertexerZ() && vertices[0]->GetDispersion() > 0.02)
2061   {
2062     AliDebugClass(AliLog::kDebug, "Vertexer z dispersion > 0.02. Not able to make a reliable multiplicity estimate.");
2063     return -2;
2064   }
2065
2066   Int_t multiplicityEstimate = 0;
2067
2068   // SPD tracklet-only estimate
2069   if (trackType == kTracklets)
2070   {
2071     const AliMultiplicity* spdmult = esd->GetMultiplicity();    // spd multiplicity object
2072     for (Int_t i=0; i<spdmult->GetNumberOfTracklets(); ++i)
2073     {
2074       if (TMath::Abs(spdmult->GetEta(i)) > etaRange)
2075         continue; // eta selection for tracklets
2076       multiplicityEstimate++;
2077     }
2078     return multiplicityEstimate;
2079   }
2080
2081   if (!vertices[1]->GetStatus())
2082   {
2083     AliDebugClass(AliLog::kDebug, "No track vertex. Not able to make a reliable multiplicity estimate.");
2084     return -3;
2085   }
2086
2087   // TODO value of displacement to be studied
2088   const Float_t maxDisplacement = 0.5;
2089   //check for displaced vertices
2090   Double_t displacement = TMath::Abs(vertices[0]->GetZ() - vertices[1]->GetZ());
2091   if (displacement > maxDisplacement)
2092   {
2093     AliDebugClass(AliLog::kDebug, Form("Displaced vertices %f > %f",displacement,maxDisplacement));
2094     return -4;
2095   }
2096
2097   // update eta range in track cuts
2098   GetMultEstTrackCuts(kMultEstTrackCutITSSA)->SetEtaRange(-etaRange, etaRange);
2099   GetMultEstTrackCuts(kMultEstTrackCutDCAwSPD)->SetEtaRange(-etaRange, etaRange);
2100   GetMultEstTrackCuts(kMultEstTrackCutDCAwoSPD)->SetEtaRange(-etaRange, etaRange);
2101
2102   //*******************************************************************************************************
2103   //set counters to initial zeros
2104   Int_t tracksITSTPC = 0;     //number of global tracks for a given event
2105   Int_t tracksITSSA = 0;      //number of ITS standalone tracks for a given event
2106   Int_t tracksITSTPCSA_complementary = 0; //number of ITS standalone tracks complementary to TPC for a given event
2107   Int_t trackletsITSTPC_complementary = 0;//number of SPD tracklets complementary to global/ITSSA tracks for a given events
2108   Int_t trackletsITSSA_complementary = 0; //number of SPD tracklets complementary to ITSSA tracks for a given event
2109
2110   const Int_t nESDTracks = esd->GetNumberOfTracks();
2111   Int_t highestID = 0;
2112
2113   // flags for secondary and rejected tracks
2114   const Int_t kRejBit = BIT(15); // set this bit in global tracks if it is rejected by a cut
2115   const Int_t kSecBit = BIT(16); // set this bit in global tracks if it is secondary according to a cut
2116
2117   for(Int_t itracks=0; itracks < nESDTracks; itracks++) {
2118     if(esd->GetTrack(itracks)->GetLabel() > highestID) highestID = esd->GetTrack(itracks)->GetLabel();
2119     esd->GetTrack(itracks)->ResetBit(kSecBit|kRejBit); //reset bits used for flagging secondaries and rejected tracks in case they were changed before this analysis
2120   }
2121   const Int_t maxid = highestID+1; // used to define bool array for check multiple associations of tracklets to one track. array starts at 0.
2122
2123   // bit mask for esd tracks, to check if multiple tracklets are associated to it
2124   Bool_t globalBits[maxid], pureITSBits[maxid];
2125   for(Int_t i=0; i<maxid; i++){ // set all bools to false
2126       globalBits[i]=kFALSE;
2127       pureITSBits[i]=kFALSE;
2128   }
2129
2130   //*******************************************************************************************************
2131   // get multiplicity from global tracks
2132   for(Int_t itracks = 0; itracks < nESDTracks; itracks++) { // flag the tracks
2133     AliESDtrack* track = esd->GetTrack(itracks);
2134
2135     // if track is a secondary from a V0, flag as a secondary
2136     if (track->IsOn(AliESDtrack::kMultInV0)) {
2137       track->SetBit(kSecBit);
2138       continue;
2139     }
2140
2141     //secondary?
2142     if (track->IsOn(AliESDtrack::kMultSec)) {
2143       track->SetBit(kSecBit);
2144       continue;
2145     }
2146
2147     // check tracks with ITS part
2148     //*******************************************************************************************************
2149     if (track->IsOn(AliESDtrack::kITSin) && !track->IsOn(AliESDtrack::kITSpureSA) && trackType == kTrackletsITSTPC) { // track has ITS part but is not an ITS_SA
2150       //*******************************************************************************************************
2151       // TPC+ITS
2152       if (track->IsOn(AliESDtrack::kTPCin)) {  // Global track, has ITS and TPC contributions
2153           if (fgMultEstTrackCuts[kMultEstTrackCutGlobal]->AcceptTrack(track)) { // good ITSTPC track
2154             if (fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->AcceptTrack(track) || fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->AcceptTrack(track)) {
2155               tracksITSTPC++; //global track counted
2156               globalBits[itracks] = kTRUE;
2157             }
2158             else track->SetBit(kSecBit); // large DCA -> secondary, don't count either track not associated tracklet
2159           }
2160           else track->SetBit(kRejBit); // bad quality, don't count the track, but may count tracklet if associated
2161       }
2162       //*******************************************************************************************************
2163       // ITS complementary
2164       else if (fgMultEstTrackCuts[kMultEstTrackCutITSSA]->AcceptTrack(track)) { // good ITS complementary track
2165         if (fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->AcceptTrack(track) || fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->AcceptTrack(track)) {
2166           tracksITSTPCSA_complementary++;
2167           globalBits[itracks] = kTRUE;
2168         }
2169         else track->SetBit(kSecBit); // large DCA -> secondary, don't count either track not associated tracklet
2170       }
2171       else track->SetBit(kRejBit); // bad quality, don't count the track, but may count tracklet if associated
2172     }
2173     //*******************************************************************************************************
2174     // check tracks from ITS_SA_PURE
2175     if (track->IsOn(AliESDtrack::kITSin) && track->IsOn(AliESDtrack::kITSpureSA) && trackType == kTrackletsITSSA){
2176       if (fgMultEstTrackCuts[kMultEstTrackCutITSSA]->AcceptTrack(track)) { // good ITSSA track
2177         if (fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->AcceptTrack(track) || fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->AcceptTrack(track)) {
2178           tracksITSSA++;
2179           pureITSBits[itracks] = kTRUE;
2180         }
2181         else track->SetBit(kRejBit);
2182       }
2183       else track->SetBit(kRejBit);
2184     }
2185   }//ESD tracks counted
2186
2187   //*******************************************************************************************************
2188   // get multiplicity from ITS tracklets to complement TPC+ITS, and ITSpureSA
2189   const AliMultiplicity* spdmult = esd->GetMultiplicity();    // spd multiplicity object
2190   for (Int_t i=0; i<spdmult->GetNumberOfTracklets(); ++i) {
2191     if (TMath::Abs(spdmult->GetEta(i)) > etaRange) continue; // eta selection for tracklets
2192
2193     // if counting tracks+tracklets, check if clusters were already used in tracks
2194     Int_t id1,id2,id3,id4;
2195     spdmult->GetTrackletTrackIDs(i,0,id1,id2); // references for eventual Global/ITS_SA tracks
2196     AliESDtrack* tr1 = id1>=0 ? esd->GetTrack(id1) : 0;
2197     spdmult->GetTrackletTrackIDs(i,1,id3,id4); // references for eventual ITS_SA_pure tracks
2198     AliESDtrack* tr3 = id3>=0 ? esd->GetTrack(id3) : 0;
2199
2200     // are both clusters from the same tracks? If not, skip the tracklet (shouldn't change things much)
2201     if ((id1!=id2 && id1>=0 && id2>=0) || (id3!=id4 && id3>=0 && id4>=0)) continue;
2202
2203     Bool_t bUsedInGlobal = (id1 != -1) ? globalBits[id1] : 0;// has associated global track been associated to a previous tracklet?
2204     Bool_t bUsedInPureITS = (id3 != -1) ? pureITSBits[id3] : 0;// has associated pure ITS track been associated to a previous tracklet?
2205     //*******************************************************************************************************
2206     if (trackType == kTrackletsITSTPC) {
2207       // count tracklets towards global+complementary tracks
2208       if (  (tr1 && !tr1->TestBit(kSecBit)) &&   // reject as secondary
2209         (tr1 &&  tr1->TestBit(kRejBit)) ) {  // count tracklet as bad quality track
2210           if(!bUsedInGlobal){
2211             ++trackletsITSTPC_complementary;
2212             if(id1>0) globalBits[id1] = kTRUE; // mark global track linked to this tracklet as "associated"
2213           }
2214       }
2215       else if(id1<0) {
2216         ++trackletsITSTPC_complementary; // if no associated track, count the tracklet
2217       }
2218     } else {
2219       // count tracklets towards ITS_SA_pure tracks
2220       if (  (tr3 && !tr3->TestBit(kSecBit)) &&   // reject as secondary
2221         (tr3 &&  tr3->TestBit(kRejBit)) ) { // count tracklet as bad quality track
2222         if(!bUsedInPureITS) {
2223           ++trackletsITSSA_complementary;
2224           if(id3>0) pureITSBits[id3] = kTRUE; // mark global track linked to this tracklet as "associated"
2225         }
2226       }
2227       else if(id3<0) {
2228         ++trackletsITSSA_complementary; // if no associated track, count the tracklet
2229       }
2230     }
2231   }
2232
2233   //*******************************************************************************************************
2234   if (trackType == kTrackletsITSTPC)
2235     multiplicityEstimate = tracksITSTPC + tracksITSTPCSA_complementary + trackletsITSTPC_complementary;
2236   else
2237     multiplicityEstimate = tracksITSSA + trackletsITSSA_complementary;
2238
2239   return multiplicityEstimate;
2240 }