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::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     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::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)
 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       const AliESDEvent* esdEvent = esdTrack->GetESDEvent();
1260
1261       if (!esdEvent)
1262         AliFatal("fCutMaxChi2TPCConstrainedVsGlobal set but ESD event not set in AliESDTrack. Use AliESDTrack::SetESDEvent before calling AliESDtrackCuts.");
1263
1264       // get vertex
1265       const AliESDVertex* vertex = 0;
1266       if (fCutMaxChi2TPCConstrainedVsGlobalVertexType & kVertexTracks)
1267         vertex = esdEvent->GetPrimaryVertexTracks();
1268
1269       if ((!vertex || !vertex->GetStatus()) && fCutMaxChi2TPCConstrainedVsGlobalVertexType & kVertexSPD)
1270         vertex = esdEvent->GetPrimaryVertexSPD();
1271
1272       if ((!vertex || !vertex->GetStatus()) && fCutMaxChi2TPCConstrainedVsGlobalVertexType & kVertexTPC)
1273         vertex = esdEvent->GetPrimaryVertexTPC();
1274
1275       if (vertex->GetStatus())
1276         chi2TPCConstrainedVsGlobal = esdTrack->GetChi2TPCConstrainedVsGlobal(vertex);
1277
1278       if (chi2TPCConstrainedVsGlobal < 0 || chi2TPCConstrainedVsGlobal > fCutMaxChi2TPCConstrainedVsGlobal)
1279       {
1280         cuts[39] = kTRUE;
1281         cut = kTRUE;
1282       }
1283     }
1284   }
1285
1286   //########################################################################
1287   // filling histograms
1288   if (fHistogramsOn) {
1289     fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n tracks")));
1290     if (cut)
1291       fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut tracks")));
1292
1293     for (Int_t i=0; i<kNCuts; i++) {
1294       if (fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i]) < 1)
1295         AliFatal(Form("Inconsistency! Cut %d with name %s not found", i, fgkCutNames[i]));
1296
1297       if (cuts[i])
1298         fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
1299
1300       for (Int_t j=i; j<kNCuts; j++) {
1301         if (cuts[i] && cuts[j]) {
1302           Float_t xC = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
1303           Float_t yC = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
1304           fhCutCorrelation->Fill(xC, yC);
1305         }
1306       }
1307     }
1308   }
1309
1310   // now we loop over the filling of the histograms twice: once "before" the cut, once "after"
1311   // the code is not in a function due to too many local variables that would need to be passed
1312
1313   for (Int_t id = 0; id < 2; id++)
1314   {
1315     // id = 0 --> before cut
1316     // id = 1 --> after cut
1317
1318     if (fHistogramsOn)
1319     {
1320       fhNClustersITS[id]->Fill(nClustersITS);
1321       fhNClustersTPC[id]->Fill(nClustersTPC);
1322       fhNSharedClustersTPC[id]->Fill(nClustersTPCShared);
1323       fhNCrossedRowsTPC[id]->Fill(nCrossedRowsTPC);
1324       fhRatioCrossedRowsOverFindableClustersTPC[id]->Fill(ratioCrossedRowsOverFindableClustersTPC);
1325       fhChi2PerClusterITS[id]->Fill(chi2PerClusterITS);
1326       fhChi2PerClusterTPC[id]->Fill(chi2PerClusterTPC);
1327       fhChi2TPCConstrainedVsGlobal[id]->Fill(chi2TPCConstrainedVsGlobal);
1328       fhNClustersForITSPID[id]->Fill(nITSPointsForPid);
1329       fhNMissingITSPoints[id]->Fill(nMissITSpts);
1330
1331       fhC11[id]->Fill(extCov[0]);
1332       fhC22[id]->Fill(extCov[2]);
1333       fhC33[id]->Fill(extCov[5]);
1334       fhC44[id]->Fill(extCov[9]);
1335       fhC55[id]->Fill(extCov[14]);
1336
1337       fhRel1PtUncertainty[id]->Fill(relUncertainty1Pt);
1338
1339       fhPt[id]->Fill(pt);
1340       fhEta[id]->Fill(eta);
1341
1342       Float_t bRes[2];
1343       bRes[0] = TMath::Sqrt(bCov[0]);
1344       bRes[1] = TMath::Sqrt(bCov[2]);
1345
1346       fhDZ[id]->Fill(b[1]);
1347       fhDXY[id]->Fill(b[0]);
1348       fhDXYDZ[id]->Fill(dcaToVertex);
1349       fhDXYvsDZ[id]->Fill(b[1],b[0]);
1350
1351       if (bRes[0]!=0 && bRes[1]!=0) {
1352         fhDZNormalized[id]->Fill(b[1]/bRes[1]);
1353         fhDXYNormalized[id]->Fill(b[0]/bRes[0]);
1354         fhDXYvsDZNormalized[id]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
1355         fhNSigmaToVertex[id]->Fill(nSigmaToVertex);
1356       }
1357     }
1358
1359     // cut the track
1360     if (cut)
1361       return kFALSE;
1362   }
1363
1364   return kTRUE;
1365 }
1366
1367 //____________________________________________________________________
1368 Bool_t AliESDtrackCuts::CheckITSClusterRequirement(ITSClusterRequirement req, Bool_t clusterL1, Bool_t clusterL2)
1369 {
1370   // checks if the cluster requirement is fullfilled (in this case: return kTRUE)
1371
1372   switch (req)
1373   {
1374         case kOff:        return kTRUE;
1375         case kNone:       return !clusterL1 && !clusterL2;
1376         case kAny:        return clusterL1 || clusterL2;
1377         case kFirst:      return clusterL1;
1378         case kOnlyFirst:  return clusterL1 && !clusterL2;
1379         case kSecond:     return clusterL2;
1380         case kOnlySecond: return clusterL2 && !clusterL1;
1381         case kBoth:       return clusterL1 && clusterL2;
1382   }
1383
1384   return kFALSE;
1385 }
1386
1387 //____________________________________________________________________
1388 AliESDtrack* AliESDtrackCuts::GetTPCOnlyTrack(const AliESDEvent* esd, Int_t iTrack)
1389 {
1390   // Utility function to create a TPC only track from the given esd track
1391   //
1392   // IMPORTANT: The track has to be deleted by the user
1393   //
1394   // NB. most of the functionality to get a TPC only track from an ESD track is in AliESDtrack, where it should be
1395   // there are only missing propagations here that are needed for old data
1396   // this function will therefore become obsolete
1397   //
1398   // adapted from code provided by CKB
1399
1400   if (!esd->GetPrimaryVertexTPC())
1401     return 0; // No TPC vertex no TPC tracks
1402
1403   if(!esd->GetPrimaryVertexTPC()->GetStatus())
1404     return 0; // TPC Vertex is created by default in AliESDEvent, do not use in this case
1405
1406   AliESDtrack* track = esd->GetTrack(iTrack);
1407   if (!track)
1408     return 0;
1409
1410   AliESDtrack *tpcTrack = new AliESDtrack();
1411
1412   // only true if we have a tpc track
1413   if (!track->FillTPCOnlyTrack(*tpcTrack))
1414   {
1415     delete tpcTrack;
1416     return 0;
1417   }
1418
1419   return tpcTrack;
1420 }
1421
1422 //____________________________________________________________________
1423 TObjArray* AliESDtrackCuts::GetAcceptedTracks(const AliESDEvent* esd, Bool_t bTPC)
1424 {
1425   //
1426   // returns an array of all tracks that pass the cuts
1427   // or an array of TPC only tracks (propagated to the TPC vertex during reco)
1428   // tracks that pass the cut
1429   //
1430   // NOTE: List has to be deleted by the user
1431
1432   TObjArray* acceptedTracks = new TObjArray();
1433
1434   // loop over esd tracks
1435   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
1436     if(bTPC){
1437       if(!esd->GetPrimaryVertexTPC())return acceptedTracks; // No TPC vertex no TPC tracks
1438       if(!esd->GetPrimaryVertexTPC()->GetStatus())return acceptedTracks; // No proper TPC vertex, only the default
1439
1440       AliESDtrack *tpcTrack = GetTPCOnlyTrack(esd, iTrack);
1441       if (!tpcTrack)
1442         continue;
1443
1444       if (AcceptTrack(tpcTrack)) {
1445         acceptedTracks->Add(tpcTrack);
1446       }
1447       else
1448         delete tpcTrack;
1449     }
1450     else
1451     {
1452       AliESDtrack* track = esd->GetTrack(iTrack);
1453       if(AcceptTrack(track))
1454         acceptedTracks->Add(track);
1455     }
1456   }
1457   if(bTPC)acceptedTracks->SetOwner(kTRUE);
1458   return acceptedTracks;
1459 }
1460
1461 //____________________________________________________________________
1462 Int_t AliESDtrackCuts::CountAcceptedTracks(const AliESDEvent* const esd)
1463 {
1464   //
1465   // returns an the number of tracks that pass the cuts
1466   //
1467
1468   Int_t count = 0;
1469
1470   // loop over esd tracks
1471   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
1472     AliESDtrack* track = esd->GetTrack(iTrack);
1473     if (AcceptTrack(track))
1474       count++;
1475   }
1476
1477   return count;
1478 }
1479
1480 //____________________________________________________________________
1481  void AliESDtrackCuts::DefineHistograms(Int_t color) {
1482    //
1483    // diagnostics histograms are defined
1484    //
1485
1486    fHistogramsOn=kTRUE;
1487
1488    Bool_t oldStatus = TH1::AddDirectoryStatus();
1489    TH1::AddDirectory(kFALSE);
1490
1491    //###################################################################################
1492    // defining histograms
1493
1494    fhCutStatistics = new TH1F("cut_statistics","cut statistics",kNCuts+4,-0.5,kNCuts+3.5);
1495
1496    fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
1497    fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
1498
1499    fhCutCorrelation = new TH2F("cut_correlation","cut correlation",kNCuts,-0.5,kNCuts-0.5,kNCuts,-0.5,kNCuts-0.5);;
1500
1501    for (Int_t i=0; i<kNCuts; i++) {
1502      fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fgkCutNames[i]);
1503      fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fgkCutNames[i]);
1504      fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fgkCutNames[i]);
1505    }
1506
1507   fhCutStatistics  ->SetLineColor(color);
1508   fhCutCorrelation ->SetLineColor(color);
1509   fhCutStatistics  ->SetLineWidth(2);
1510   fhCutCorrelation ->SetLineWidth(2);
1511
1512   for (Int_t i=0; i<2; i++) {
1513     fhNClustersITS[i]        = new TH1F("nClustersITS"    ,"",8,-0.5,7.5);
1514     fhNClustersTPC[i]        = new TH1F("nClustersTPC"     ,"",165,-0.5,164.5);
1515     fhNSharedClustersTPC[i]  = new TH1F("nSharedClustersTPC"     ,"",165,-0.5,164.5);
1516     fhNCrossedRowsTPC[i]     = new TH1F("nCrossedRowsTPC"  ,"",165,-0.5,164.5);
1517     fhRatioCrossedRowsOverFindableClustersTPC[i]     = new TH1F("ratioCrossedRowsOverFindableClustersTPC"  ,"",60,0,1.5);
1518     fhChi2PerClusterITS[i]   = new TH1F("chi2PerClusterITS","",500,0,10);
1519     fhChi2PerClusterTPC[i]   = new TH1F("chi2PerClusterTPC","",500,0,10);
1520     fhChi2TPCConstrainedVsGlobal[i]   = new TH1F("chi2TPCConstrainedVsGlobal","",600,-2,50);
1521     fhNClustersForITSPID[i]  = new TH1F("nPointsForITSpid","",5,-0.5,4.5);
1522     fhNMissingITSPoints[i]   = new TH1F("nMissingITSClusters","",7,-0.5,6.5);
1523
1524     fhC11[i]                 = new TH1F("covMatrixDiagonal11","",2000,0,20);
1525     fhC22[i]                 = new TH1F("covMatrixDiagonal22","",2000,0,20);
1526     fhC33[i]                 = new TH1F("covMatrixDiagonal33","",1000,0,0.1);
1527     fhC44[i]                 = new TH1F("covMatrixDiagonal44","",1000,0,0.1);
1528     fhC55[i]                 = new TH1F("covMatrixDiagonal55","",1000,0,5);
1529
1530     fhRel1PtUncertainty[i]   = new TH1F("rel1PtUncertainty","",1000,0,5);
1531
1532     fhDXY[i]                 = new TH1F("dXY"    ,"",500,-10,10);
1533     fhDZ[i]                  = new TH1F("dZ"     ,"",500,-10,10);
1534     fhDXYDZ[i]               = new TH1F("dXYDZ"  ,"",500,0,10);
1535     fhDXYvsDZ[i]             = new TH2F("dXYvsDZ","",200,-10,10,200,-10,10);
1536
1537     fhDXYNormalized[i]       = new TH1F("dXYNormalized"    ,"",500,-10,10);
1538     fhDZNormalized[i]        = new TH1F("dZNormalized"     ,"",500,-10,10);
1539     fhDXYvsDZNormalized[i]   = new TH2F("dXYvsDZNormalized","",200,-10,10,200,-10,10);
1540
1541     fhNSigmaToVertex[i]      = new TH1F("nSigmaToVertex","",500,0,10);
1542
1543     fhPt[i]                  = new TH1F("pt"     ,"p_{T} distribution;p_{T} (GeV/c)", 800, 0.0, 10.0);
1544     fhEta[i]                 = new TH1F("eta"     ,"#eta distribution;#eta",40,-2.0,2.0);
1545
1546     fhNClustersITS[i]->SetTitle("n ITS clusters");
1547     fhNClustersTPC[i]->SetTitle("n TPC clusters");
1548     fhNSharedClustersTPC[i]->SetTitle("n TPC shared clusters");
1549     fhChi2PerClusterITS[i]->SetTitle("#Chi^{2} per ITS cluster");
1550     fhChi2PerClusterTPC[i]->SetTitle("#Chi^{2} per TPC cluster");
1551     fhChi2TPCConstrainedVsGlobal[i]->SetTitle("#Chi^{2} TPC constrained track vs global track");
1552     fhNClustersForITSPID[i]->SetTitle("n ITS points for PID");
1553     fhNMissingITSPoints[i]->SetTitle("n ITS layers with missing cluster");
1554
1555     fhC11[i]->SetTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
1556     fhC22[i]->SetTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
1557     fhC33[i]->SetTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
1558     fhC44[i]->SetTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
1559     fhC55[i]->SetTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
1560
1561     fhRel1PtUncertainty[i]->SetTitle("rel. uncertainty of 1/p_{T}");
1562
1563     fhDXY[i]->SetXTitle("transverse impact parameter (cm)");
1564     fhDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
1565     fhDXYDZ[i]->SetTitle("absolute impact parameter;sqrt(dXY**2 + dZ**2)  (cm)");
1566     fhDXYvsDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
1567     fhDXYvsDZ[i]->SetYTitle("transverse impact parameter (cm)");
1568
1569     fhDXYNormalized[i]->SetTitle("normalized trans impact par (n#sigma)");
1570     fhDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
1571     fhDXYvsDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
1572     fhDXYvsDZNormalized[i]->SetYTitle("normalized trans impact par (n#sigma)");
1573     fhNSigmaToVertex[i]->SetTitle("n #sigma to vertex");
1574
1575     fhNClustersITS[i]->SetLineColor(color);   fhNClustersITS[i]->SetLineWidth(2);
1576     fhNClustersTPC[i]->SetLineColor(color);   fhNClustersTPC[i]->SetLineWidth(2);
1577     fhNSharedClustersTPC[i]->SetLineColor(color);  fhNSharedClustersTPC[i]->SetLineWidth(2);
1578     fhChi2PerClusterITS[i]->SetLineColor(color);   fhChi2PerClusterITS[i]->SetLineWidth(2);
1579     fhChi2PerClusterTPC[i]->SetLineColor(color);   fhChi2PerClusterTPC[i]->SetLineWidth(2);
1580     fhChi2TPCConstrainedVsGlobal[i]->SetLineColor(color);   fhChi2TPCConstrainedVsGlobal[i]->SetLineWidth(2);
1581     fhNClustersForITSPID[i]->SetLineColor(color);  fhNClustersForITSPID[i]->SetLineWidth(2);
1582     fhNMissingITSPoints[i]->SetLineColor(color);   fhNMissingITSPoints[i]->SetLineWidth(2);
1583
1584     fhC11[i]->SetLineColor(color);   fhC11[i]->SetLineWidth(2);
1585     fhC22[i]->SetLineColor(color);   fhC22[i]->SetLineWidth(2);
1586     fhC33[i]->SetLineColor(color);   fhC33[i]->SetLineWidth(2);
1587     fhC44[i]->SetLineColor(color);   fhC44[i]->SetLineWidth(2);
1588     fhC55[i]->SetLineColor(color);   fhC55[i]->SetLineWidth(2);
1589
1590     fhRel1PtUncertainty[i]->SetLineColor(color); fhRel1PtUncertainty[i]->SetLineWidth(2);
1591
1592     fhDXY[i]->SetLineColor(color);   fhDXY[i]->SetLineWidth(2);
1593     fhDZ[i]->SetLineColor(color);    fhDZ[i]->SetLineWidth(2);
1594     fhDXYDZ[i]->SetLineColor(color); fhDXYDZ[i]->SetLineWidth(2);
1595
1596     fhDXYNormalized[i]->SetLineColor(color);   fhDXYNormalized[i]->SetLineWidth(2);
1597     fhDZNormalized[i]->SetLineColor(color);    fhDZNormalized[i]->SetLineWidth(2);
1598     fhNSigmaToVertex[i]->SetLineColor(color);  fhNSigmaToVertex[i]->SetLineWidth(2);
1599   }
1600
1601   // The number of sigmas to the vertex is per definition gaussian
1602   ffDTheoretical = new TF1("nSigmaToVertexTheoretical","([0]/2.506628274)*exp(-(x**2)/2)",0,50);
1603   ffDTheoretical->SetParameter(0,1);
1604
1605   TH1::AddDirectory(oldStatus);
1606 }
1607
1608 //____________________________________________________________________
1609 Bool_t AliESDtrackCuts::LoadHistograms(const Char_t* dir)
1610 {
1611   //
1612   // loads the histograms from a file
1613   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
1614   //
1615
1616   if (!dir)
1617     dir = GetName();
1618
1619   if (!gDirectory->cd(dir))
1620     return kFALSE;
1621
1622   ffDTheoretical = dynamic_cast<TF1*> (gDirectory->Get("nSigmaToVertexTheory"));
1623
1624   fhCutStatistics = dynamic_cast<TH1F*> (gDirectory->Get("cut_statistics"));
1625   fhCutCorrelation = dynamic_cast<TH2F*> (gDirectory->Get("cut_correlation"));
1626
1627   for (Int_t i=0; i<2; i++) {
1628     if (i==0)
1629     {
1630       gDirectory->cd("before_cuts");
1631     }
1632     else
1633       gDirectory->cd("after_cuts");
1634
1635     fhNClustersITS[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersITS"     ));
1636     fhNClustersTPC[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersTPC"     ));
1637     fhNSharedClustersTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get("nSharedClustersTPC"     ));
1638     fhNCrossedRowsTPC[i]   = dynamic_cast<TH1F*> (gDirectory->Get("nCrossedRowsTPC"  ));
1639     fhRatioCrossedRowsOverFindableClustersTPC[i]   = dynamic_cast<TH1F*> (gDirectory->Get("ratioCrossedRowsOverFindableClustersTPC"  ));
1640     fhChi2PerClusterITS[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterITS"));
1641     fhChi2PerClusterTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterTPC"));
1642     fhChi2TPCConstrainedVsGlobal[i] = dynamic_cast<TH1F*> (gDirectory->Get("fhChi2TPCConstrainedVsGlobal"));
1643     fhNClustersForITSPID[i] = dynamic_cast<TH1F*> (gDirectory->Get("nPointsForITSpid"));
1644     fhNMissingITSPoints[i] = dynamic_cast<TH1F*> (gDirectory->Get("nMissingITSClusters"));
1645
1646     fhC11[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal11"));
1647     fhC22[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal22"));
1648     fhC33[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal33"));
1649     fhC44[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal44"));
1650     fhC55[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal55"));
1651
1652     fhRel1PtUncertainty[i] = dynamic_cast<TH1F*> (gDirectory->Get("rel1PtUncertainty"));
1653
1654     fhDXY[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXY"    ));
1655     fhDZ[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZ"     ));
1656     fhDXYDZ[i] =   dynamic_cast<TH1F*> (gDirectory->Get("dXYDZ"));
1657     fhDXYvsDZ[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZ"));
1658
1659     fhDXYNormalized[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXYNormalized"    ));
1660     fhDZNormalized[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZNormalized"     ));
1661     fhDXYvsDZNormalized[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZNormalized"));
1662     fhNSigmaToVertex[i] = dynamic_cast<TH1F*> (gDirectory->Get("nSigmaToVertex"));
1663
1664     fhPt[i] = dynamic_cast<TH1F*> (gDirectory->Get("pt"));
1665     fhEta[i] = dynamic_cast<TH1F*> (gDirectory->Get("eta"));
1666
1667     gDirectory->cd("../");
1668   }
1669
1670   gDirectory->cd("..");
1671
1672   return kTRUE;
1673 }
1674
1675 //____________________________________________________________________
1676 void AliESDtrackCuts::SaveHistograms(const Char_t* dir) {
1677   //
1678   // saves the histograms in a directory (dir)
1679   //
1680
1681   if (!fHistogramsOn) {
1682     AliDebug(0, "Histograms not on - cannot save histograms!!!");
1683     return;
1684   }
1685
1686   if (!dir)
1687     dir = GetName();
1688
1689   gDirectory->mkdir(dir);
1690   gDirectory->cd(dir);
1691
1692   gDirectory->mkdir("before_cuts");
1693   gDirectory->mkdir("after_cuts");
1694
1695   // a factor of 2 is needed since n sigma is positive
1696   ffDTheoretical->SetParameter(0,2*fhNSigmaToVertex[0]->Integral("width"));
1697   ffDTheoretical->Write("nSigmaToVertexTheory");
1698
1699   fhCutStatistics->Write();
1700   fhCutCorrelation->Write();
1701
1702   for (Int_t i=0; i<2; i++) {
1703     if (i==0)
1704       gDirectory->cd("before_cuts");
1705     else
1706       gDirectory->cd("after_cuts");
1707
1708     fhNClustersITS[i]        ->Write();
1709     fhNClustersTPC[i]        ->Write();
1710     fhNSharedClustersTPC[i]  ->Write();
1711     fhNCrossedRowsTPC[i]     ->Write();
1712     fhRatioCrossedRowsOverFindableClustersTPC[i]     ->Write();
1713     fhChi2PerClusterITS[i]   ->Write();
1714     fhChi2PerClusterTPC[i]   ->Write();
1715     fhChi2TPCConstrainedVsGlobal[i]   ->Write();
1716     fhNClustersForITSPID[i]  ->Write();
1717     fhNMissingITSPoints[i]   ->Write();
1718
1719     fhC11[i]                 ->Write();
1720     fhC22[i]                 ->Write();
1721     fhC33[i]                 ->Write();
1722     fhC44[i]                 ->Write();
1723     fhC55[i]                 ->Write();
1724
1725     fhRel1PtUncertainty[i]   ->Write();
1726
1727     fhDXY[i]                 ->Write();
1728     fhDZ[i]                  ->Write();
1729     fhDXYDZ[i]               ->Write();
1730     fhDXYvsDZ[i]             ->Write();
1731
1732     fhDXYNormalized[i]       ->Write();
1733     fhDZNormalized[i]        ->Write();
1734     fhDXYvsDZNormalized[i]   ->Write();
1735     fhNSigmaToVertex[i]      ->Write();
1736
1737     fhPt[i]                  ->Write();
1738     fhEta[i]                 ->Write();
1739
1740     gDirectory->cd("../");
1741   }
1742
1743   gDirectory->cd("../");
1744 }
1745
1746 //____________________________________________________________________
1747 void AliESDtrackCuts::DrawHistograms()
1748 {
1749   // draws some histograms
1750
1751   TCanvas* canvas1 = new TCanvas(Form("%s_1", GetName()), "Track Quality Results1", 800, 800);
1752   canvas1->Divide(2, 2);
1753
1754   canvas1->cd(1);
1755   fhNClustersTPC[0]->SetStats(kFALSE);
1756   fhNClustersTPC[0]->Draw();
1757
1758   canvas1->cd(2);
1759   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1760   fhChi2PerClusterTPC[0]->Draw();
1761
1762   canvas1->cd(3);
1763   fhNSigmaToVertex[0]->SetStats(kFALSE);
1764   fhNSigmaToVertex[0]->GetXaxis()->SetRangeUser(0, 10);
1765   fhNSigmaToVertex[0]->Draw();
1766
1767   canvas1->SaveAs(Form("%s_%s.gif", GetName(), canvas1->GetName()));
1768
1769   TCanvas* canvas2 = new TCanvas(Form("%s_2", GetName()), "Track Quality Results2", 1200, 800);
1770   canvas2->Divide(3, 2);
1771
1772   canvas2->cd(1);
1773   fhC11[0]->SetStats(kFALSE);
1774   gPad->SetLogy();
1775   fhC11[0]->Draw();
1776
1777   canvas2->cd(2);
1778   fhC22[0]->SetStats(kFALSE);
1779   gPad->SetLogy();
1780   fhC22[0]->Draw();
1781
1782   canvas2->cd(3);
1783   fhC33[0]->SetStats(kFALSE);
1784   gPad->SetLogy();
1785   fhC33[0]->Draw();
1786
1787   canvas2->cd(4);
1788   fhC44[0]->SetStats(kFALSE);
1789   gPad->SetLogy();
1790   fhC44[0]->Draw();
1791
1792   canvas2->cd(5);
1793   fhC55[0]->SetStats(kFALSE);
1794   gPad->SetLogy();
1795   fhC55[0]->Draw();
1796
1797   canvas2->cd(6);
1798   fhRel1PtUncertainty[0]->SetStats(kFALSE);
1799   gPad->SetLogy();
1800   fhRel1PtUncertainty[0]->Draw();
1801
1802   canvas2->SaveAs(Form("%s_%s.gif", GetName(), canvas2->GetName()));
1803
1804   TCanvas* canvas3 = new TCanvas(Form("%s_3", GetName()), "Track Quality Results3", 1200, 800);
1805   canvas3->Divide(3, 2);
1806
1807   canvas3->cd(1);
1808   fhDXY[0]->SetStats(kFALSE);
1809   gPad->SetLogy();
1810   fhDXY[0]->Draw();
1811
1812   canvas3->cd(2);
1813   fhDZ[0]->SetStats(kFALSE);
1814   gPad->SetLogy();
1815   fhDZ[0]->Draw();
1816
1817   canvas3->cd(3);
1818   fhDXYvsDZ[0]->SetStats(kFALSE);
1819   gPad->SetLogz();
1820   gPad->SetRightMargin(0.15);
1821   fhDXYvsDZ[0]->Draw("COLZ");
1822
1823   canvas3->cd(4);
1824   fhDXYNormalized[0]->SetStats(kFALSE);
1825   gPad->SetLogy();
1826   fhDXYNormalized[0]->Draw();
1827
1828   canvas3->cd(5);
1829   fhDZNormalized[0]->SetStats(kFALSE);
1830   gPad->SetLogy();
1831   fhDZNormalized[0]->Draw();
1832
1833   canvas3->cd(6);
1834   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1835   gPad->SetLogz();
1836   gPad->SetRightMargin(0.15);
1837   fhDXYvsDZNormalized[0]->Draw("COLZ");
1838
1839   canvas3->SaveAs(Form("%s_%s.gif", GetName(), canvas3->GetName()));
1840
1841   TCanvas* canvas4 = new TCanvas(Form("%s_4", GetName()), "Track Quality Results4", 800, 500);
1842   canvas4->Divide(2, 1);
1843
1844   canvas4->cd(1);
1845   fhCutStatistics->SetStats(kFALSE);
1846   fhCutStatistics->LabelsOption("v");
1847   gPad->SetBottomMargin(0.3);
1848   fhCutStatistics->Draw();
1849
1850   canvas4->cd(2);
1851   fhCutCorrelation->SetStats(kFALSE);
1852   fhCutCorrelation->LabelsOption("v");
1853   gPad->SetBottomMargin(0.3);
1854   gPad->SetLeftMargin(0.3);
1855   fhCutCorrelation->Draw("COLZ");
1856
1857   canvas4->SaveAs(Form("%s_%s.gif", GetName(), canvas4->GetName()));
1858
1859   /*canvas->cd(1);
1860   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1861   fhDXYvsDZNormalized[0]->DrawCopy("COLZ");
1862
1863   canvas->cd(2);
1864   fhNClustersTPC[0]->SetStats(kFALSE);
1865   fhNClustersTPC[0]->DrawCopy();
1866
1867   canvas->cd(3);
1868   fhChi2PerClusterITS[0]->SetStats(kFALSE);
1869   fhChi2PerClusterITS[0]->DrawCopy();
1870   fhChi2PerClusterITS[1]->SetLineColor(2);
1871   fhChi2PerClusterITS[1]->DrawCopy("SAME");
1872
1873   canvas->cd(4);
1874   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1875   fhChi2PerClusterTPC[0]->DrawCopy();
1876   fhChi2PerClusterTPC[1]->SetLineColor(2);
1877   fhChi2PerClusterTPC[1]->DrawCopy("SAME");*/
1878 }
1879 //--------------------------------------------------------------------------
1880 void AliESDtrackCuts::SetPtDepDCACuts(Double_t pt) {
1881   //
1882   // set the pt-dependent DCA cuts
1883   //
1884
1885   if(f1CutMaxDCAToVertexXYPtDep) {
1886      fCutMaxDCAToVertexXY=f1CutMaxDCAToVertexXYPtDep->Eval(pt);
1887   }
1888
1889   if(f1CutMaxDCAToVertexZPtDep) {
1890     fCutMaxDCAToVertexZ=f1CutMaxDCAToVertexZPtDep->Eval(pt);
1891   }
1892
1893   if(f1CutMinDCAToVertexXYPtDep) {
1894     fCutMinDCAToVertexXY=f1CutMinDCAToVertexXYPtDep->Eval(pt);
1895   }
1896
1897   if(f1CutMinDCAToVertexZPtDep) {
1898     fCutMinDCAToVertexZ=f1CutMinDCAToVertexZPtDep->Eval(pt);
1899   }
1900
1901
1902   return;
1903 }
1904
1905
1906
1907 //--------------------------------------------------------------------------
1908 Bool_t AliESDtrackCuts::CheckPtDepDCA(TString dist,Bool_t print) const {
1909   //
1910   // Check the correctness of the string syntax
1911   //
1912   Bool_t retval=kTRUE;
1913
1914   if(!dist.Contains("pt")) {
1915     if(print) AliError("string must contain \"pt\"");
1916     retval= kFALSE;
1917   }
1918   return retval;
1919 }
1920
1921  void AliESDtrackCuts::SetMaxDCAToVertexXYPtDep(const char *dist){
1922
1923    if(f1CutMaxDCAToVertexXYPtDep){
1924      delete f1CutMaxDCAToVertexXYPtDep;
1925      // resetiing both
1926      f1CutMaxDCAToVertexXYPtDep = 0;
1927      fCutMaxDCAToVertexXYPtDep = "";
1928    }
1929    if(!CheckPtDepDCA(dist,kTRUE)){
1930      return;
1931    }
1932    fCutMaxDCAToVertexXYPtDep = dist;
1933    TString tmp(dist);
1934    tmp.ReplaceAll("pt","x");
1935    f1CutMaxDCAToVertexXYPtDep = new TFormula("f1CutMaxDCAToVertexXYPtDep",tmp.Data());
1936
1937 }
1938
1939  void AliESDtrackCuts::SetMaxDCAToVertexZPtDep(const char *dist){
1940
1941
1942    if(f1CutMaxDCAToVertexZPtDep){
1943      delete f1CutMaxDCAToVertexZPtDep;
1944      // resetiing both
1945      f1CutMaxDCAToVertexZPtDep = 0;
1946      fCutMaxDCAToVertexZPtDep = "";
1947    }
1948    if(!CheckPtDepDCA(dist,kTRUE))return;
1949
1950    fCutMaxDCAToVertexZPtDep = dist;
1951    TString tmp(dist);
1952    tmp.ReplaceAll("pt","x");
1953    f1CutMaxDCAToVertexZPtDep = new TFormula("f1CutMaxDCAToVertexZPtDep",tmp.Data());
1954
1955
1956 }
1957
1958
1959  void AliESDtrackCuts::SetMinDCAToVertexXYPtDep(const char *dist){
1960
1961
1962    if(f1CutMinDCAToVertexXYPtDep){
1963      delete f1CutMinDCAToVertexXYPtDep;
1964      // resetiing both
1965      f1CutMinDCAToVertexXYPtDep = 0;
1966      fCutMinDCAToVertexXYPtDep = "";
1967    }
1968    if(!CheckPtDepDCA(dist,kTRUE))return;
1969
1970    fCutMinDCAToVertexXYPtDep = dist;
1971    TString tmp(dist);
1972    tmp.ReplaceAll("pt","x");
1973    f1CutMinDCAToVertexXYPtDep = new TFormula("f1CutMinDCAToVertexXYPtDep",tmp.Data());
1974
1975 }
1976
1977
1978  void AliESDtrackCuts::SetMinDCAToVertexZPtDep(const char *dist){
1979
1980
1981
1982    if(f1CutMinDCAToVertexZPtDep){
1983      delete f1CutMinDCAToVertexZPtDep;
1984      // resetiing both
1985      f1CutMinDCAToVertexZPtDep = 0;
1986      fCutMinDCAToVertexZPtDep = "";
1987    }
1988    if(!CheckPtDepDCA(dist,kTRUE))return;
1989    fCutMinDCAToVertexZPtDep = dist;
1990    TString tmp(dist);
1991    tmp.ReplaceAll("pt","x");
1992    f1CutMinDCAToVertexZPtDep = new TFormula("f1CutMinDCAToVertexZPtDep",tmp.Data());
1993 }
1994
1995 AliESDtrackCuts* AliESDtrackCuts::GetMultEstTrackCuts(MultEstTrackCuts cut)
1996 {
1997   // returns the multiplicity estimator track cuts objects to allow for user configuration
1998   // upon first call the objects are created
1999   //
2000   // the cut defined here correspond to GetStandardITSTPCTrackCuts2010 (apart from the one for without SPD)
2001
2002   if (!fgMultEstTrackCuts[kMultEstTrackCutGlobal])
2003   {
2004     // quality cut on ITS+TPC tracks
2005     fgMultEstTrackCuts[kMultEstTrackCutGlobal] = new AliESDtrackCuts();
2006     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetMinNClustersTPC(70);
2007     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetMaxChi2PerClusterTPC(4);
2008     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetAcceptKinkDaughters(kFALSE);
2009     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetRequireTPCRefit(kTRUE);
2010     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetRequireITSRefit(kTRUE);
2011     //multiplicity underestimate if we use global tracks with |eta| > 0.9
2012     fgMultEstTrackCuts[kMultEstTrackCutGlobal]->SetEtaRange(-0.9, 0.9);
2013
2014     // quality cut on ITS_SA tracks (complementary to ITS+TPC)
2015     fgMultEstTrackCuts[kMultEstTrackCutITSSA] = new AliESDtrackCuts();
2016     fgMultEstTrackCuts[kMultEstTrackCutITSSA]->SetRequireITSRefit(kTRUE);
2017
2018     // primary selection for tracks with SPD hits
2019     fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD] = new AliESDtrackCuts();
2020     fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
2021     fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
2022     fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->SetMaxDCAToVertexZ(2);
2023
2024     // primary selection for tracks w/o SPD hits
2025     fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD] = new AliESDtrackCuts();
2026     fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kNone);
2027     fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->SetMaxDCAToVertexXYPtDep("1.5*(0.0182+0.0350/pt^1.01)");
2028     fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->SetMaxDCAToVertexZ(2);
2029   }
2030
2031   return fgMultEstTrackCuts[cut];
2032 }
2033
2034 Int_t AliESDtrackCuts::GetReferenceMultiplicity(const AliESDEvent* esd, MultEstTrackType trackType, Float_t etaRange)
2035 {
2036   // Get multiplicity estimate based on TPC/ITS tracks and tracklets
2037   // Adapted for AliESDtrackCuts from a version developed by: Ruben Shahoyan, Anton Alkin, Arvinder Palaha
2038   //
2039   // Returns a negative value if no reliable estimate can be provided:
2040   //   -1 SPD vertex missing
2041   //   -2 SPD VertexerZ dispersion too large
2042   //   -3 Track vertex missing (not checked for kTracklets)
2043   //   -4 Distance between SPD and track vertices too large (not checked for kTracklets)
2044   //
2045   // WARNING This functions does not cut on the z vtx. Depending on the eta range requested, you need to restrict your z vertex range!
2046   //
2047   // Strategy for combined estimators
2048   //   1. Count global tracks and flag them
2049   //   2. Count ITSSA as complementaries for ITSTPC+ or as main tracks
2050   //   3. Count the complementary SPD tracklets
2051
2052   const AliESDVertex* vertices[2];
2053   vertices[0] = esd->GetPrimaryVertexSPD();
2054   vertices[1] = esd->GetPrimaryVertexTracks();
2055
2056   if (!vertices[0]->GetStatus())
2057   {
2058     AliDebugClass(AliLog::kDebug, "No SPD vertex. Not able to make a reliable multiplicity estimate.");
2059     return -1;
2060   }
2061
2062   if (vertices[0]->IsFromVertexerZ() && vertices[0]->GetDispersion() > 0.02)
2063   {
2064     AliDebugClass(AliLog::kDebug, "Vertexer z dispersion > 0.02. Not able to make a reliable multiplicity estimate.");
2065     return -2;
2066   }
2067
2068   Int_t multiplicityEstimate = 0;
2069
2070   // SPD tracklet-only estimate
2071   if (trackType == kTracklets)
2072   {
2073     const AliMultiplicity* spdmult = esd->GetMultiplicity();    // spd multiplicity object
2074     for (Int_t i=0; i<spdmult->GetNumberOfTracklets(); ++i)
2075     {
2076       if (TMath::Abs(spdmult->GetEta(i)) > etaRange)
2077         continue; // eta selection for tracklets
2078       multiplicityEstimate++;
2079     }
2080     return multiplicityEstimate;
2081   }
2082
2083   if (!vertices[1]->GetStatus())
2084   {
2085     AliDebugClass(AliLog::kDebug, "No track vertex. Not able to make a reliable multiplicity estimate.");
2086     return -3;
2087   }
2088
2089   // TODO value of displacement to be studied
2090   const Float_t maxDisplacement = 0.5;
2091   //check for displaced vertices
2092   Double_t displacement = TMath::Abs(vertices[0]->GetZ() - vertices[1]->GetZ());
2093   if (displacement > maxDisplacement)
2094   {
2095     AliDebugClass(AliLog::kDebug, Form("Displaced vertices %f > %f",displacement,maxDisplacement));
2096     return -4;
2097   }
2098
2099   // update eta range in track cuts
2100   GetMultEstTrackCuts(kMultEstTrackCutITSSA)->SetEtaRange(-etaRange, etaRange);
2101   GetMultEstTrackCuts(kMultEstTrackCutDCAwSPD)->SetEtaRange(-etaRange, etaRange);
2102   GetMultEstTrackCuts(kMultEstTrackCutDCAwoSPD)->SetEtaRange(-etaRange, etaRange);
2103
2104   //*******************************************************************************************************
2105   //set counters to initial zeros
2106   Int_t tracksITSTPC = 0;     //number of global tracks for a given event
2107   Int_t tracksITSSA = 0;      //number of ITS standalone tracks for a given event
2108   Int_t tracksITSTPCSA_complementary = 0; //number of ITS standalone tracks complementary to TPC for a given event
2109   Int_t trackletsITSTPC_complementary = 0;//number of SPD tracklets complementary to global/ITSSA tracks for a given events
2110   Int_t trackletsITSSA_complementary = 0; //number of SPD tracklets complementary to ITSSA tracks for a given event
2111
2112   const Int_t nESDTracks = esd->GetNumberOfTracks();
2113   Int_t highestID = 0;
2114
2115   // flags for secondary and rejected tracks
2116   const Int_t kRejBit = BIT(15); // set this bit in global tracks if it is rejected by a cut
2117   const Int_t kSecBit = BIT(16); // set this bit in global tracks if it is secondary according to a cut
2118
2119   for(Int_t itracks=0; itracks < nESDTracks; itracks++) {
2120     if(esd->GetTrack(itracks)->GetLabel() > highestID) highestID = esd->GetTrack(itracks)->GetLabel();
2121     esd->GetTrack(itracks)->ResetBit(kSecBit|kRejBit); //reset bits used for flagging secondaries and rejected tracks in case they were changed before this analysis
2122   }
2123   const Int_t maxid = highestID+1; // used to define bool array for check multiple associations of tracklets to one track. array starts at 0.
2124
2125   // bit mask for esd tracks, to check if multiple tracklets are associated to it
2126   Bool_t globalBits[maxid], pureITSBits[maxid];
2127   for(Int_t i=0; i<maxid; i++){ // set all bools to false
2128       globalBits[i]=kFALSE;
2129       pureITSBits[i]=kFALSE;
2130   }
2131
2132   //*******************************************************************************************************
2133   // get multiplicity from global tracks
2134   for(Int_t itracks = 0; itracks < nESDTracks; itracks++) { // flag the tracks
2135     AliESDtrack* track = esd->GetTrack(itracks);
2136
2137     // if track is a secondary from a V0, flag as a secondary
2138     if (track->IsOn(AliESDtrack::kMultInV0)) {
2139       track->SetBit(kSecBit);
2140       continue;
2141     }
2142
2143     //secondary?
2144     if (track->IsOn(AliESDtrack::kMultSec)) {
2145       track->SetBit(kSecBit);
2146       continue;
2147     }
2148
2149     // check tracks with ITS part
2150     //*******************************************************************************************************
2151     if (track->IsOn(AliESDtrack::kITSin) && !track->IsOn(AliESDtrack::kITSpureSA) && trackType == kTrackletsITSTPC) { // track has ITS part but is not an ITS_SA
2152       //*******************************************************************************************************
2153       // TPC+ITS
2154       if (track->IsOn(AliESDtrack::kTPCin)) {  // Global track, has ITS and TPC contributions
2155           if (fgMultEstTrackCuts[kMultEstTrackCutGlobal]->AcceptTrack(track)) { // good ITSTPC track
2156             if (fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->AcceptTrack(track) || fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->AcceptTrack(track)) {
2157               tracksITSTPC++; //global track counted
2158               globalBits[itracks] = kTRUE;
2159             }
2160             else track->SetBit(kSecBit); // large DCA -> secondary, don't count either track not associated tracklet
2161           }
2162           else track->SetBit(kRejBit); // bad quality, don't count the track, but may count tracklet if associated
2163       }
2164       //*******************************************************************************************************
2165       // ITS complementary
2166       else if (fgMultEstTrackCuts[kMultEstTrackCutITSSA]->AcceptTrack(track)) { // good ITS complementary track
2167         if (fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->AcceptTrack(track) || fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->AcceptTrack(track)) {
2168           tracksITSTPCSA_complementary++;
2169           globalBits[itracks] = kTRUE;
2170         }
2171         else track->SetBit(kSecBit); // large DCA -> secondary, don't count either track not associated tracklet
2172       }
2173       else track->SetBit(kRejBit); // bad quality, don't count the track, but may count tracklet if associated
2174     }
2175     //*******************************************************************************************************
2176     // check tracks from ITS_SA_PURE
2177     if (track->IsOn(AliESDtrack::kITSin) && track->IsOn(AliESDtrack::kITSpureSA) && trackType == kTrackletsITSSA){
2178       if (fgMultEstTrackCuts[kMultEstTrackCutITSSA]->AcceptTrack(track)) { // good ITSSA track
2179         if (fgMultEstTrackCuts[kMultEstTrackCutDCAwSPD]->AcceptTrack(track) || fgMultEstTrackCuts[kMultEstTrackCutDCAwoSPD]->AcceptTrack(track)) {
2180           tracksITSSA++;
2181           pureITSBits[itracks] = kTRUE;
2182         }
2183         else track->SetBit(kRejBit);
2184       }
2185       else track->SetBit(kRejBit);
2186     }
2187   }//ESD tracks counted
2188
2189   //*******************************************************************************************************
2190   // get multiplicity from ITS tracklets to complement TPC+ITS, and ITSpureSA
2191   const AliMultiplicity* spdmult = esd->GetMultiplicity();    // spd multiplicity object
2192   for (Int_t i=0; i<spdmult->GetNumberOfTracklets(); ++i) {
2193     if (TMath::Abs(spdmult->GetEta(i)) > etaRange) continue; // eta selection for tracklets
2194
2195     // if counting tracks+tracklets, check if clusters were already used in tracks
2196     Int_t id1,id2,id3,id4;
2197     spdmult->GetTrackletTrackIDs(i,0,id1,id2); // references for eventual Global/ITS_SA tracks
2198     AliESDtrack* tr1 = id1>=0 ? esd->GetTrack(id1) : 0;
2199     spdmult->GetTrackletTrackIDs(i,1,id3,id4); // references for eventual ITS_SA_pure tracks
2200     AliESDtrack* tr3 = id3>=0 ? esd->GetTrack(id3) : 0;
2201
2202     // are both clusters from the same tracks? If not, skip the tracklet (shouldn't change things much)
2203     if ((id1!=id2 && id1>=0 && id2>=0) || (id3!=id4 && id3>=0 && id4>=0)) continue;
2204
2205     Bool_t bUsedInGlobal = (id1 != -1) ? globalBits[id1] : 0;// has associated global track been associated to a previous tracklet?
2206     Bool_t bUsedInPureITS = (id3 != -1) ? pureITSBits[id3] : 0;// has associated pure ITS track been associated to a previous tracklet?
2207     //*******************************************************************************************************
2208     if (trackType == kTrackletsITSTPC) {
2209       // count tracklets towards global+complementary tracks
2210       if (  (tr1 && !tr1->TestBit(kSecBit)) &&   // reject as secondary
2211         (tr1 &&  tr1->TestBit(kRejBit)) ) {  // count tracklet as bad quality track
2212           if(!bUsedInGlobal){
2213             ++trackletsITSTPC_complementary;
2214             if(id1>0) globalBits[id1] = kTRUE; // mark global track linked to this tracklet as "associated"
2215           }
2216       }
2217       else if(id1<0) {
2218         ++trackletsITSTPC_complementary; // if no associated track, count the tracklet
2219       }
2220     } else {
2221       // count tracklets towards ITS_SA_pure tracks
2222       if (  (tr3 && !tr3->TestBit(kSecBit)) &&   // reject as secondary
2223         (tr3 &&  tr3->TestBit(kRejBit)) ) { // count tracklet as bad quality track
2224         if(!bUsedInPureITS) {
2225           ++trackletsITSSA_complementary;
2226           if(id3>0) pureITSBits[id3] = kTRUE; // mark global track linked to this tracklet as "associated"
2227         }
2228       }
2229       else if(id3<0) {
2230         ++trackletsITSSA_complementary; // if no associated track, count the tracklet
2231       }
2232     }
2233   }
2234
2235   //*******************************************************************************************************
2236   if (trackType == kTrackletsITSTPC)
2237     multiplicityEstimate = tracksITSTPC + tracksITSTPCSA_complementary + trackletsITSTPC_complementary;
2238   else
2239     multiplicityEstimate = tracksITSSA + trackletsITSSA_complementary;
2240
2241   return multiplicityEstimate;
2242 }