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