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