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