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