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