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       if (fhNClustersForITSPID[i])
 613         fhNClustersForITSPID[i]->Add(entry->fhNClustersForITSPID[i]);
 614       if (fhNMissingITSPoints[i])
 615         fhNMissingITSPoints[i] ->Add(entry->fhNMissingITSPoints[i]);
 616
 617       fhC11[i]               ->Add(entry->fhC11[i]              );
 618       fhC22[i]               ->Add(entry->fhC22[i]              );
 619       fhC33[i]               ->Add(entry->fhC33[i]              );
 620       fhC44[i]               ->Add(entry->fhC44[i]              );
 621       fhC55[i]               ->Add(entry->fhC55[i]              );
 622
 623       fhRel1PtUncertainty[i] ->Add(entry->fhRel1PtUncertainty[i]);
 624
 625       fhDXY[i]               ->Add(entry->fhDXY[i]              );
 626       fhDZ[i]                ->Add(entry->fhDZ[i]               );
 627       fhDXYDZ[i]             ->Add(entry->fhDXYDZ[i]          );
 628       fhDXYvsDZ[i]           ->Add(entry->fhDXYvsDZ[i]          );
 629
 630       fhDXYNormalized[i]     ->Add(entry->fhDXYNormalized[i]    );
 631       fhDZNormalized[i]      ->Add(entry->fhDZNormalized[i]     );
 632       fhDXYvsDZNormalized[i] ->Add(entry->fhDXYvsDZNormalized[i]);
 633       fhNSigmaToVertex[i]    ->Add(entry->fhNSigmaToVertex[i]);
 634
 635       fhPt[i]                ->Add(entry->fhPt[i]);
 636       fhEta[i]               ->Add(entry->fhEta[i]);
 637     }
 638
 639     fhCutStatistics  ->Add(entry->fhCutStatistics);
 640     fhCutCorrelation ->Add(entry->fhCutCorrelation);
 641
 642     count++;
 643   }
 644   return count+1;
 645 }
 646
 647 //____________________________________________________________________
 648 AliESDtrackCuts* AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
 649 {
 650   // creates an AliESDtrackCuts object and fills it with standard (pre data-taking) values for TPC-only cuts
 651
 652   Printf("AliESDtrackCuts::GetStandardTPCOnlyTrackCuts: Creating track cuts for TPC-only.");
 653
 654   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 655
 656   esdTrackCuts->SetMinNClustersTPC(50);
 657   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
 658   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 659
 660   esdTrackCuts->SetMaxDCAToVertexZ(3.2);
 661   esdTrackCuts->SetMaxDCAToVertexXY(2.4);
 662   esdTrackCuts->SetDCAToVertex2D(kTRUE);
 663
 664   return esdTrackCuts;
 665 }
 666
 667 //____________________________________________________________________
 668 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(Bool_t selPrimaries)
 669 {
 670   // creates an AliESDtrackCuts object and fills it with standard values for ITS-TPC cuts for pp 2009 data
 671
 672   Printf("AliESDtrackCuts::GetStandardITSTPCTrackCuts: Creating track cuts for ITS+TPC.");
 673
 674   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 675
 676   // TPC
 677   esdTrackCuts->SetRequireTPCStandAlone(kTRUE); // to get chi2 and ncls of kTPCin
 678   esdTrackCuts->SetMinNClustersTPC(70);
 679   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
 680   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 681   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 682   // ITS
 683   esdTrackCuts->SetRequireITSRefit(kTRUE);
 684   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 685                                          AliESDtrackCuts::kAny);
 686   if(selPrimaries) {
 687     // 7*(0.0050+0.0060/pt^0.9)
 688     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0350+0.0420/pt^0.9");
 689   }
 690   esdTrackCuts->SetMaxDCAToVertexZ(1.e6);
 691   esdTrackCuts->SetDCAToVertex2D(kFALSE);
 692   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 693   //esdTrackCuts->SetEtaRange(-0.8,+0.8);
 694
 695   return esdTrackCuts;
 696 }
 697
 698 //____________________________________________________________________
 699 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(Bool_t selPrimaries,Int_t clusterCut)
 700 {
 701   // creates an AliESDtrackCuts object and fills it with standard values for ITS-TPC cuts for pp 2010 data
 702   // if clusterCut = 1, the cut on the number of clusters is replaced by
 703   // a cut on the number of crossed rows and on the ration crossed
 704   // rows/findable clusters
 705
 706   Printf("AliESDtrackCuts::GetStandardITSTPCTrackCuts: Creating track cuts for ITS+TPC.");
 707
 708   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 709
 710   // TPC
 711   if(clusterCut == 0)  esdTrackCuts->SetMinNClustersTPC(70);
 712   else if (clusterCut == 1) {
 713     esdTrackCuts->SetMinNCrossedRowsTPC(70);
 714     esdTrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
 715   }
 716   else {
 717     Printf("Wrong value of the clusterCut parameter (%d), using cut on Nclusters",clusterCut);
 718     esdTrackCuts->SetMinNClustersTPC(70);
 719   }
 720   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
 721   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 722   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 723   // ITS
 724   esdTrackCuts->SetRequireITSRefit(kTRUE);
 725   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 726                                          AliESDtrackCuts::kAny);
 727   if(selPrimaries) {
 728     // 7*(0.0026+0.0050/pt^1.01)
 729     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
 730   }
 731   esdTrackCuts->SetMaxDCAToVertexZ(2);
 732   esdTrackCuts->SetDCAToVertex2D(kFALSE);
 733   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 734
 735   return esdTrackCuts;
 736 }
 737
 738 //____________________________________________________________________
 739
 740
 741
 742 //____________________________________________________________________
 743 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSPureSATrackCuts2009(Bool_t selPrimaries, Bool_t useForPid)
 744 {
 745   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks
 746
 747   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 748   esdTrackCuts->SetRequireITSPureStandAlone(kTRUE);
 749   esdTrackCuts->SetRequireITSRefit(kTRUE);
 750   esdTrackCuts->SetMinNClustersITS(4);
 751   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 752                                          AliESDtrackCuts::kAny);
 753   esdTrackCuts->SetMaxChi2PerClusterITS(1.);
 754
 755   if(selPrimaries) {
 756     // 7*(0.0085+0.0026/pt^1.55)
 757     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0595+0.0182/pt^1.55");
 758   }
 759   if(useForPid){
 760     esdTrackCuts->SetRequireITSPid(kTRUE);
 761   }
 762   return esdTrackCuts;
 763 }
 764
 765 //____________________________________________________________________
 766 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSPureSATrackCuts2010(Bool_t selPrimaries, Bool_t useForPid)
 767 {
 768   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks - pp 2010
 769
 770   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 771   esdTrackCuts->SetRequireITSPureStandAlone(kTRUE);
 772   esdTrackCuts->SetRequireITSRefit(kTRUE);
 773   esdTrackCuts->SetMinNClustersITS(4);
 774   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 775                                          AliESDtrackCuts::kAny);
 776   esdTrackCuts->SetMaxChi2PerClusterITS(2.5);
 777
 778   if(selPrimaries) {
 779     // 7*(0.0033+0.0045/pt^1.3)
 780     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0231+0.0315/pt^1.3");
 781   }
 782   if(useForPid){
 783     esdTrackCuts->SetRequireITSPid(kTRUE);
 784   }
 785   return esdTrackCuts;
 786 }
 787
 788 //____________________________________________________________________
 789 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSSATrackCuts2009(Bool_t selPrimaries, Bool_t useForPid)
 790 {
 791   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks
 792
 793   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 794   esdTrackCuts->SetRequireITSStandAlone(kTRUE);
 795   esdTrackCuts->SetRequireITSPureStandAlone(kFALSE);
 796   esdTrackCuts->SetRequireITSRefit(kTRUE);
 797   esdTrackCuts->SetMinNClustersITS(4);
 798   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 799                                          AliESDtrackCuts::kAny);
 800   esdTrackCuts->SetMaxChi2PerClusterITS(1.);
 801
 802   if(selPrimaries) {
 803     // 7*(0.0085+0.0026/pt^1.55)
 804     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0595+0.0182/pt^1.55");
 805   }
 806   if(useForPid){
 807     esdTrackCuts->SetRequireITSPid(kTRUE);
 808   }
 809   return esdTrackCuts;
 810 }
 811
 812 //____________________________________________________________________
 813 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSSATrackCuts2010(Bool_t selPrimaries, Bool_t useForPid)
 814 {
 815   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks --pp 2010
 816
 817   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts;
 818   esdTrackCuts->SetRequireITSStandAlone(kTRUE);
 819   esdTrackCuts->SetRequireITSPureStandAlone(kFALSE);
 820   esdTrackCuts->SetRequireITSRefit(kTRUE);
 821   esdTrackCuts->SetMinNClustersITS(4);
 822   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 823                                          AliESDtrackCuts::kAny);
 824   esdTrackCuts->SetMaxChi2PerClusterITS(2.5);
 825
 826   if(selPrimaries) {
 827     // 7*(0.0033+0.0045/pt^1.3)
 828     esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0231+0.0315/pt^1.3");
 829   }
 830   if(useForPid){
 831     esdTrackCuts->SetRequireITSPid(kTRUE);
 832   }
 833   return esdTrackCuts;
 834 }
 835
 836 //____________________________________________________________________
 837 AliESDtrackCuts* AliESDtrackCuts::GetStandardITSSATrackCutsPbPb2010(Bool_t selPrimaries, Bool_t useForPid)
 838 {
 839   // creates an AliESDtrackCuts object and fills it with standard values for ITS pure SA tracks -- PbPb 2010
 840
 841   AliESDtrackCuts* esdTrackCuts = GetStandardITSSATrackCuts2010(selPrimaries, useForPid);
 842   esdTrackCuts->SetMaxNOfMissingITSPoints(1);
 843
 844   return esdTrackCuts;
 845 }
 846
 847 //____________________________________________________________________
 848 Int_t AliESDtrackCuts::GetReferenceMultiplicity(const AliESDEvent* esd, Bool_t tpcOnly)
 849 {
 850   // Gets reference multiplicity following the standard cuts and a defined fiducial volume
 851   // tpcOnly = kTRUE -> consider TPC-only tracks
 852   //         = kFALSE -> consider global tracks
 853
 854   if (!tpcOnly)
 855   {
 856     Printf("AliESDtrackCuts::GetReferenceMultiplicity: Not implemented for global tracks!");
 857     return -1;
 858   }
 859
 860   static AliESDtrackCuts* esdTrackCuts = 0;
 861   if (!esdTrackCuts)
 862   {
 863     esdTrackCuts = GetStandardTPCOnlyTrackCuts();
 864     esdTrackCuts->SetEtaRange(-0.8, 0.8);
 865     esdTrackCuts->SetPtRange(0.15);
 866   }
 867
 868   Int_t nTracks = esdTrackCuts->CountAcceptedTracks(esd);
 869
 870   return nTracks;
 871 }
 872
 873 //____________________________________________________________________
 874 Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* const esdTrack)
 875 {
 876   // Calculates the number of sigma to the vertex.
 877
 878   Float_t b[2];
 879   Float_t bRes[2];
 880   Float_t bCov[3];
 881   esdTrack->GetImpactParameters(b,bCov);
 882
 883   if (bCov[0]<=0 || bCov[2]<=0) {
 884     AliDebugClass(1, "Estimated b resolution lower or equal zero!");
 885     bCov[0]=0; bCov[2]=0;
 886   }
 887   bRes[0] = TMath::Sqrt(bCov[0]);
 888   bRes[1] = TMath::Sqrt(bCov[2]);
 889
 890   // -----------------------------------
 891   // How to get to a n-sigma cut?
 892   //
 893   // The accumulated statistics from 0 to d is
 894   //
 895   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 896   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 897   //
 898   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
 899   // Can this be expressed in a different way?
 900
 901   if (bRes[0] == 0 || bRes[1] ==0)
 902     return -1;
 903
 904   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 905
 906   // work around precision problem
 907   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 908   // 1e-15 corresponds to nsigma ~ 7.7
 909   if (TMath::Exp(-d * d / 2) < 1e-15)
 910     return 1000;
 911
 912   Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 913   return nSigma;
 914 }
 915
 916 void AliESDtrackCuts::EnableNeededBranches(TTree* tree)
 917 {
 918   // enables the branches needed by AcceptTrack, for a list see comment of AcceptTrack
 919
 920   tree->SetBranchStatus("fTracks.fFlags", 1);
 921   tree->SetBranchStatus("fTracks.fITSncls", 1);
 922   tree->SetBranchStatus("fTracks.fTPCncls", 1);
 923   tree->SetBranchStatus("fTracks.fITSchi2", 1);
 924   tree->SetBranchStatus("fTracks.fTPCchi2", 1);
 925   tree->SetBranchStatus("fTracks.fC*", 1);
 926   tree->SetBranchStatus("fTracks.fD", 1);
 927   tree->SetBranchStatus("fTracks.fZ", 1);
 928   tree->SetBranchStatus("fTracks.fCdd", 1);
 929   tree->SetBranchStatus("fTracks.fCdz", 1);
 930   tree->SetBranchStatus("fTracks.fCzz", 1);
 931   tree->SetBranchStatus("fTracks.fP*", 1);
 932   tree->SetBranchStatus("fTracks.fR*", 1);
 933   tree->SetBranchStatus("fTracks.fKinkIndexes*", 1);
 934 }
 935
 936 //____________________________________________________________________
 937 Bool_t AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack)
 938 {
 939   //
 940   // figure out if the tracks survives all the track cuts defined
 941   //
 942   // the different quality parameter and kinematic values are first
 943   // retrieved from the track. then it is found out what cuts the
 944   // track did not survive and finally the cuts are imposed.
 945
 946   // this function needs the following branches:
 947   // fTracks.fFlags
 948   // fTracks.fITSncls
 949   // fTracks.fTPCncls
 950   // fTracks.fITSchi2
 951   // fTracks.fTPCchi2
 952   // fTracks.fC   //GetExternalCovariance
 953   // fTracks.fD   //GetImpactParameters
 954   // fTracks.fZ   //GetImpactParameters
 955   // fTracks.fCdd //GetImpactParameters
 956   // fTracks.fCdz //GetImpactParameters
 957   // fTracks.fCzz //GetImpactParameters
 958   // fTracks.fP   //GetPxPyPz
 959   // fTracks.fR   //GetMass
 960   // fTracks.fP   //GetMass
 961   // fTracks.fKinkIndexes
 962
 963   UInt_t status = esdTrack->GetStatus();
 964
 965   // getting quality parameters from the ESD track
 966   Int_t nClustersITS = esdTrack->GetITSclusters(0);
 967   Int_t nClustersTPC = -1;
 968   if(fCutRequireTPCStandAlone) {
 969     nClustersTPC = esdTrack->GetTPCNclsIter1();
 970   }
 971   else {
 972     nClustersTPC = esdTrack->GetTPCclusters(0);
 973   }
 974   Float_t nCrossedRowsTPC = esdTrack->GetTPCClusterInfo(2,1);
 975   Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 976   if (esdTrack->GetTPCNclsF()>0) {
 977     ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC / esdTrack->GetTPCNclsF();
 978   }
 979
 980   Int_t nClustersTPCShared = esdTrack->GetTPCnclsS();
 981   Float_t fracClustersTPCShared = -1.;
 982
 983   Float_t chi2PerClusterITS = -1;
 984   Float_t chi2PerClusterTPC = -1;
 985   if (nClustersITS!=0)
 986     chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
 987   if (nClustersTPC!=0) {
 988     if(fCutRequireTPCStandAlone) {
 989       chi2PerClusterTPC = esdTrack->GetTPCchi2Iter1()/Float_t(nClustersTPC);
 990     } else {
 991       chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
 992     }
 993     fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClustersTPC);
 994   }
 995
 996   Double_t extCov[15];
 997   esdTrack->GetExternalCovariance(extCov);
 998
 999   // getting the track to vertex parameters
1000   Float_t nSigmaToVertex = GetSigmaToVertex(esdTrack);
1001
1002   Float_t b[2];
1003   Float_t bCov[3];
1004   esdTrack->GetImpactParameters(b,bCov);
1005   if (bCov[0]<=0 || bCov[2]<=0) {
1006     AliDebug(1, "Estimated b resolution lower or equal zero!");
1007     bCov[0]=0; bCov[2]=0;
1008   }
1009
1010
1011   // set pt-dependent DCA cuts, if requested
1012   SetPtDepDCACuts(esdTrack->Pt());
1013
1014
1015   Float_t dcaToVertexXY = b[0];
1016   Float_t dcaToVertexZ = b[1];
1017
1018   Float_t dcaToVertex = -1;
1019
1020   if (fCutDCAToVertex2D)
1021   {
1022     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1023   }
1024   else
1025     dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1026
1027   // getting the kinematic variables of the track
1028   // (assuming the mass is known)
1029   Double_t p[3];
1030   esdTrack->GetPxPyPz(p);
1031
1032   Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
1033   Float_t pt       = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2));
1034   Float_t energy   = TMath::Sqrt(TMath::Power(esdTrack->GetMass(),2) + TMath::Power(momentum,2));
1035
1036   //y-eta related calculations
1037   Float_t eta = -100.;
1038   Float_t y   = -100.;
1039   if((momentum != TMath::Abs(p[2]))&&(momentum != 0))
1040     eta = 0.5*TMath::Log((momentum + p[2])/(momentum - p[2]));
1041   if((energy != TMath::Abs(p[2]))&&(momentum != 0))
1042     y = 0.5*TMath::Log((energy + p[2])/(energy - p[2]));
1043
1044   if (extCov[14] < 0)
1045   {
1046     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]);
1047     return kFALSE;
1048   }
1049   Float_t relUncertainty1Pt = TMath::Sqrt(extCov[14])*pt;
1050
1051   //########################################################################
1052   // cut the track?
1053
1054   Bool_t cuts[kNCuts];
1055   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1056
1057   // track quality cuts
1058   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1059     cuts[0]=kTRUE;
1060   if (fCutRequireTPCStandAlone && (status&AliESDtrack::kTPCin)==0)
1061     cuts[1]=kTRUE;
1062   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1063     cuts[2]=kTRUE;
1064   if (nClustersTPC<fCutMinNClusterTPC)
1065     cuts[3]=kTRUE;
1066   if (nClustersITS<fCutMinNClusterITS)
1067     cuts[4]=kTRUE;
1068   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1069     cuts[5]=kTRUE;
1070   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1071     cuts[6]=kTRUE;
1072   if (extCov[0]  > fCutMaxC11)
1073     cuts[7]=kTRUE;
1074   if (extCov[2]  > fCutMaxC22)
1075     cuts[8]=kTRUE;
1076   if (extCov[5]  > fCutMaxC33)
1077     cuts[9]=kTRUE;
1078   if (extCov[9]  > fCutMaxC44)
1079     cuts[10]=kTRUE;
1080   if (extCov[14]  > fCutMaxC55)
1081     cuts[11]=kTRUE;
1082   if (nSigmaToVertex > fCutNsigmaToVertex && fCutSigmaToVertexRequired)
1083     cuts[12] = kTRUE;
1084   // if n sigma could not be calculated
1085   if (nSigmaToVertex<0 && fCutSigmaToVertexRequired)
1086     cuts[13]=kTRUE;
1087   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1088     cuts[14]=kTRUE;
1089   // track kinematics cut
1090   if((momentum < fPMin) || (momentum > fPMax))
1091     cuts[15]=kTRUE;
1092   if((pt < fPtMin) || (pt > fPtMax))
1093     cuts[16] = kTRUE;
1094   if((p[0] < fPxMin) || (p[0] > fPxMax))
1095     cuts[17] = kTRUE;
1096   if((p[1] < fPyMin) || (p[1] > fPyMax))
1097     cuts[18] = kTRUE;
1098   if((p[2] < fPzMin) || (p[2] > fPzMax))
1099     cuts[19] = kTRUE;
1100   if((eta < fEtaMin) || (eta > fEtaMax))
1101     cuts[20] = kTRUE;
1102   if((y < fRapMin) || (y > fRapMax))
1103     cuts[21] = kTRUE;
1104   if (fCutDCAToVertex2D && dcaToVertex > 1)
1105     cuts[22] = kTRUE;
1106   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1107     cuts[23] = kTRUE;
1108   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1109     cuts[24] = kTRUE;
1110   if (fCutDCAToVertex2D && fCutMinDCAToVertexXY > 0 && fCutMinDCAToVertexZ > 0 && dcaToVertexXY*dcaToVertexXY/fCutMinDCAToVertexXY/fCutMinDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMinDCAToVertexZ/fCutMinDCAToVertexZ < 1)
1111     cuts[25] = kTRUE;
1112   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) < fCutMinDCAToVertexXY)
1113     cuts[26] = kTRUE;
1114   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) < fCutMinDCAToVertexZ)
1115     cuts[27] = kTRUE;
1116
1117   for (Int_t i = 0; i < 3; i++)
1118     cuts[28+i] = !CheckITSClusterRequirement(fCutClusterRequirementITS[i], esdTrack->HasPointOnITSLayer(i*2), esdTrack->HasPointOnITSLayer(i*2+1));
1119
1120   if(fCutRequireITSStandAlone || fCutRequireITSpureSA){
1121     if ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin)){
1122       // TPC tracks
1123       cuts[31] = kTRUE;
1124     }else{
1125       // ITS standalone tracks
1126       if(fCutRequireITSStandAlone && !fCutRequireITSpureSA){
1127         if(status & AliESDtrack::kITSpureSA) cuts[31] = kTRUE;
1128       }else if(fCutRequireITSpureSA){
1129         if(!(status & AliESDtrack::kITSpureSA)) cuts[31] = kTRUE;
1130       }
1131     }
1132   }
1133
1134   if (relUncertainty1Pt > fCutMaxRel1PtUncertainty)
1135      cuts[32] = kTRUE;
1136
1137   if (!fCutAcceptSharedTPCClusters && nClustersTPCShared!=0)
1138     cuts[33] = kTRUE;
1139
1140   if (fracClustersTPCShared > fCutMaxFractionSharedTPCClusters)
1141     cuts[34] = kTRUE;
1142
1143   Int_t nITSPointsForPid=0;
1144   UChar_t clumap=esdTrack->GetITSClusterMap();
1145   for(Int_t i=2; i<6; i++){
1146     if(clumap&(1<<i)) ++nITSPointsForPid;
1147   }
1148   if(fCutRequireITSPid && nITSPointsForPid<3) cuts[35] = kTRUE;
1149
1150
1151   if (nCrossedRowsTPC<fCutMinNCrossedRowsTPC)
1152     cuts[36]=kTRUE;
1153   if (ratioCrossedRowsOverFindableClustersTPC<fCutMinRatioCrossedRowsOverFindableClustersTPC)
1154     cuts[37]=kTRUE;
1155
1156   Int_t nMissITSpts=0;
1157   Int_t idet,statusLay;
1158   Float_t xloc,zloc;
1159   for(Int_t iLay=0; iLay<6; iLay++){
1160     Bool_t retc=esdTrack->GetITSModuleIndexInfo(iLay,idet,statusLay,xloc,zloc);
1161     if(retc && statusLay==5) ++nMissITSpts;
1162   }
1163   if(nMissITSpts>fCutMaxMissingITSPoints) cuts[38] = kTRUE;
1164
1165   Bool_t cut=kFALSE;
1166   for (Int_t i=0; i<kNCuts; i++)
1167     if (cuts[i]) {cut = kTRUE;}
1168
1169
1170   //########################################################################
1171   // filling histograms
1172   if (fHistogramsOn) {
1173     fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n tracks")));
1174     if (cut)
1175       fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut tracks")));
1176
1177     for (Int_t i=0; i<kNCuts; i++) {
1178       if (fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i]) < 1)
1179         AliFatal(Form("Inconsistency! Cut %d with name %s not found", i, fgkCutNames[i]));
1180
1181       if (cuts[i])
1182         fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
1183
1184       for (Int_t j=i; j<kNCuts; j++) {
1185         if (cuts[i] && cuts[j]) {
1186           Float_t xC = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
1187           Float_t yC = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
1188           fhCutCorrelation->Fill(xC, yC);
1189         }
1190       }
1191     }
1192   }
1193
1194   // now we loop over the filling of the histograms twice: once "before" the cut, once "after"
1195   // the code is not in a function due to too many local variables that would need to be passed
1196
1197   for (Int_t id = 0; id < 2; id++)
1198   {
1199     // id = 0 --> before cut
1200     // id = 1 --> after cut
1201
1202     if (fHistogramsOn)
1203     {
1204       fhNClustersITS[id]->Fill(nClustersITS);
1205       fhNClustersTPC[id]->Fill(nClustersTPC);
1206       fhNSharedClustersTPC[id]->Fill(nClustersTPCShared);
1207       fhNCrossedRowsTPC[id]->Fill(nCrossedRowsTPC);
1208       fhRatioCrossedRowsOverFindableClustersTPC[id]->Fill(ratioCrossedRowsOverFindableClustersTPC);
1209       fhChi2PerClusterITS[id]->Fill(chi2PerClusterITS);
1210       fhChi2PerClusterTPC[id]->Fill(chi2PerClusterTPC);
1211       fhNClustersForITSPID[id]->Fill(nITSPointsForPid);
1212       fhNMissingITSPoints[id]->Fill(nMissITSpts);
1213
1214       fhC11[id]->Fill(extCov[0]);
1215       fhC22[id]->Fill(extCov[2]);
1216       fhC33[id]->Fill(extCov[5]);
1217       fhC44[id]->Fill(extCov[9]);
1218       fhC55[id]->Fill(extCov[14]);
1219
1220       fhRel1PtUncertainty[id]->Fill(relUncertainty1Pt);
1221
1222       fhPt[id]->Fill(pt);
1223       fhEta[id]->Fill(eta);
1224
1225       Float_t bRes[2];
1226       bRes[0] = TMath::Sqrt(bCov[0]);
1227       bRes[1] = TMath::Sqrt(bCov[2]);
1228
1229       fhDZ[id]->Fill(b[1]);
1230       fhDXY[id]->Fill(b[0]);
1231       fhDXYDZ[id]->Fill(dcaToVertex);
1232       fhDXYvsDZ[id]->Fill(b[1],b[0]);
1233
1234       if (bRes[0]!=0 && bRes[1]!=0) {
1235         fhDZNormalized[id]->Fill(b[1]/bRes[1]);
1236         fhDXYNormalized[id]->Fill(b[0]/bRes[0]);
1237         fhDXYvsDZNormalized[id]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
1238         fhNSigmaToVertex[id]->Fill(nSigmaToVertex);
1239       }
1240     }
1241
1242     // cut the track
1243     if (cut)
1244       return kFALSE;
1245   }
1246
1247   return kTRUE;
1248 }
1249
1250 //____________________________________________________________________
1251 Bool_t AliESDtrackCuts::CheckITSClusterRequirement(ITSClusterRequirement req, Bool_t clusterL1, Bool_t clusterL2)
1252 {
1253   // checks if the cluster requirement is fullfilled (in this case: return kTRUE)
1254
1255   switch (req)
1256   {
1257         case kOff:        return kTRUE;
1258         case kNone:       return !clusterL1 && !clusterL2;
1259         case kAny:        return clusterL1 || clusterL2;
1260         case kFirst:      return clusterL1;
1261         case kOnlyFirst:  return clusterL1 && !clusterL2;
1262         case kSecond:     return clusterL2;
1263         case kOnlySecond: return clusterL2 && !clusterL1;
1264         case kBoth:       return clusterL1 && clusterL2;
1265   }
1266
1267   return kFALSE;
1268 }
1269
1270 //____________________________________________________________________
1271 AliESDtrack* AliESDtrackCuts::GetTPCOnlyTrack(AliESDEvent* esd, Int_t iTrack)
1272 {
1273
1274   // Utility function to
1275   // create a TPC only track from the given esd track
1276   //
1277   // IMPORTANT: The track has to be deleted by the user
1278   //
1279   // NB. most of the functionality to get a TPC only track from an ESD track is in AliESDtrack, where it should be
1280   // there are only missing propagations here that are needed for old data
1281   // this function will therefore become obsolete
1282   //
1283   // adapted from code provided by CKB
1284
1285   if (!esd->GetPrimaryVertexTPC())
1286     return 0; // No TPC vertex no TPC tracks
1287
1288   if(!esd->GetPrimaryVertexTPC()->GetStatus())
1289     return 0; // TPC Vertex is created by default in AliESDEvent, do not use in this case
1290
1291   AliESDtrack* track = esd->GetTrack(iTrack);
1292   if (!track)
1293     return 0;
1294
1295   AliESDtrack *tpcTrack = new AliESDtrack();
1296
1297   // only true if we have a tpc track
1298   if (!track->FillTPCOnlyTrack(*tpcTrack))
1299   {
1300     delete tpcTrack;
1301     return 0;
1302   }
1303
1304   // propagate to Vertex
1305   // not needed for normal reconstructed ESDs...
1306   // Double_t pTPC[2],covTPC[3];
1307   // tpcTrack->PropagateToDCA(esd->GetPrimaryVertexTPC(), esd->GetMagneticField(), 10000,  pTPC, covTPC);
1308
1309   return tpcTrack;
1310 }
1311
1312 //____________________________________________________________________
1313 TObjArray* AliESDtrackCuts::GetAcceptedTracks(AliESDEvent* esd,Bool_t bTPC)
1314 {
1315   //
1316   // returns an array of all tracks that pass the cuts
1317   // or an array of TPC only tracks (propagated to the TPC vertex during reco)
1318   // tracks that pass the cut
1319   //
1320   // NOTE: List has to be deleted by the user
1321
1322   TObjArray* acceptedTracks = new TObjArray();
1323
1324   // loop over esd tracks
1325   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
1326     if(bTPC){
1327       if(!esd->GetPrimaryVertexTPC())return acceptedTracks; // No TPC vertex no TPC tracks
1328       if(!esd->GetPrimaryVertexTPC()->GetStatus())return acceptedTracks; // No proper TPC vertex, only the default
1329
1330       AliESDtrack *tpcTrack = GetTPCOnlyTrack(esd, iTrack);
1331       if (!tpcTrack)
1332         continue;
1333
1334       if (AcceptTrack(tpcTrack)) {
1335         acceptedTracks->Add(tpcTrack);
1336       }
1337       else
1338         delete tpcTrack;
1339     }
1340     else
1341     {
1342       AliESDtrack* track = esd->GetTrack(iTrack);
1343       if(AcceptTrack(track))
1344         acceptedTracks->Add(track);
1345     }
1346   }
1347   if(bTPC)acceptedTracks->SetOwner(kTRUE);
1348   return acceptedTracks;
1349 }
1350
1351 //____________________________________________________________________
1352 Int_t AliESDtrackCuts::CountAcceptedTracks(const AliESDEvent* const esd)
1353 {
1354   //
1355   // returns an the number of tracks that pass the cuts
1356   //
1357
1358   Int_t count = 0;
1359
1360   // loop over esd tracks
1361   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
1362     AliESDtrack* track = esd->GetTrack(iTrack);
1363     if (AcceptTrack(track))
1364       count++;
1365   }
1366
1367   return count;
1368 }
1369
1370 //____________________________________________________________________
1371  void AliESDtrackCuts::DefineHistograms(Int_t color) {
1372    //
1373    // diagnostics histograms are defined
1374    //
1375
1376    fHistogramsOn=kTRUE;
1377
1378    Bool_t oldStatus = TH1::AddDirectoryStatus();
1379    TH1::AddDirectory(kFALSE);
1380
1381    //###################################################################################
1382    // defining histograms
1383
1384    fhCutStatistics = new TH1F("cut_statistics","cut statistics",kNCuts+4,-0.5,kNCuts+3.5);
1385
1386    fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
1387    fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
1388
1389    fhCutCorrelation = new TH2F("cut_correlation","cut correlation",kNCuts,-0.5,kNCuts-0.5,kNCuts,-0.5,kNCuts-0.5);;
1390
1391    for (Int_t i=0; i<kNCuts; i++) {
1392      fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fgkCutNames[i]);
1393      fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fgkCutNames[i]);
1394      fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fgkCutNames[i]);
1395    }
1396
1397   fhCutStatistics  ->SetLineColor(color);
1398   fhCutCorrelation ->SetLineColor(color);
1399   fhCutStatistics  ->SetLineWidth(2);
1400   fhCutCorrelation ->SetLineWidth(2);
1401
1402   for (Int_t i=0; i<2; i++) {
1403     fhNClustersITS[i]        = new TH1F("nClustersITS"    ,"",8,-0.5,7.5);
1404     fhNClustersTPC[i]        = new TH1F("nClustersTPC"     ,"",165,-0.5,164.5);
1405     fhNSharedClustersTPC[i]  = new TH1F("nSharedClustersTPC"     ,"",165,-0.5,164.5);
1406     fhNCrossedRowsTPC[i]     = new TH1F("nCrossedRowsTPC"  ,"",165,-0.5,164.5);
1407     fhRatioCrossedRowsOverFindableClustersTPC[i]     = new TH1F("ratioCrossedRowsOverFindableClustersTPC"  ,"",60,0,1.5);
1408     fhChi2PerClusterITS[i]   = new TH1F("chi2PerClusterITS","",500,0,10);
1409     fhChi2PerClusterTPC[i]   = new TH1F("chi2PerClusterTPC","",500,0,10);
1410     fhNClustersForITSPID[i]  = new TH1F("nPointsForITSpid","",5,-0.5,4.5);
1411     fhNMissingITSPoints[i]   = new TH1F("nMissingITSClusters","",7,-0.5,6.5);
1412
1413     fhC11[i]                 = new TH1F("covMatrixDiagonal11","",2000,0,20);
1414     fhC22[i]                 = new TH1F("covMatrixDiagonal22","",2000,0,20);
1415     fhC33[i]                 = new TH1F("covMatrixDiagonal33","",1000,0,0.1);
1416     fhC44[i]                 = new TH1F("covMatrixDiagonal44","",1000,0,0.1);
1417     fhC55[i]                 = new TH1F("covMatrixDiagonal55","",1000,0,5);
1418
1419     fhRel1PtUncertainty[i]   = new TH1F("rel1PtUncertainty","",1000,0,5);
1420
1421     fhDXY[i]                 = new TH1F("dXY"    ,"",500,-10,10);
1422     fhDZ[i]                  = new TH1F("dZ"     ,"",500,-10,10);
1423     fhDXYDZ[i]               = new TH1F("dXYDZ"  ,"",500,0,10);
1424     fhDXYvsDZ[i]             = new TH2F("dXYvsDZ","",200,-10,10,200,-10,10);
1425
1426     fhDXYNormalized[i]       = new TH1F("dXYNormalized"    ,"",500,-10,10);
1427     fhDZNormalized[i]        = new TH1F("dZNormalized"     ,"",500,-10,10);
1428     fhDXYvsDZNormalized[i]   = new TH2F("dXYvsDZNormalized","",200,-10,10,200,-10,10);
1429
1430     fhNSigmaToVertex[i]      = new TH1F("nSigmaToVertex","",500,0,10);
1431
1432     fhPt[i]                  = new TH1F("pt"     ,"p_{T} distribution;p_{T} (GeV/c)", 800, 0.0, 10.0);
1433     fhEta[i]                 = new TH1F("eta"     ,"#eta distribution;#eta",40,-2.0,2.0);
1434
1435     fhNClustersITS[i]->SetTitle("n ITS clusters");
1436     fhNClustersTPC[i]->SetTitle("n TPC clusters");
1437     fhNSharedClustersTPC[i]->SetTitle("n TPC shared clusters");
1438     fhChi2PerClusterITS[i]->SetTitle("#Chi^{2} per ITS cluster");
1439     fhChi2PerClusterTPC[i]->SetTitle("#Chi^{2} per TPC cluster");
1440     fhNClustersForITSPID[i]->SetTitle("n ITS points for PID");
1441     fhNMissingITSPoints[i]->SetTitle("n ITS layers with missing cluster");
1442
1443     fhC11[i]->SetTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
1444     fhC22[i]->SetTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
1445     fhC33[i]->SetTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
1446     fhC44[i]->SetTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
1447     fhC55[i]->SetTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
1448
1449     fhRel1PtUncertainty[i]->SetTitle("rel. uncertainty of 1/p_{T}");
1450
1451     fhDXY[i]->SetXTitle("transverse impact parameter (cm)");
1452     fhDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
1453     fhDXYDZ[i]->SetTitle("absolute impact parameter;sqrt(dXY**2 + dZ**2)  (cm)");
1454     fhDXYvsDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
1455     fhDXYvsDZ[i]->SetYTitle("transverse impact parameter (cm)");
1456
1457     fhDXYNormalized[i]->SetTitle("normalized trans impact par (n#sigma)");
1458     fhDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
1459     fhDXYvsDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
1460     fhDXYvsDZNormalized[i]->SetYTitle("normalized trans impact par (n#sigma)");
1461     fhNSigmaToVertex[i]->SetTitle("n #sigma to vertex");
1462
1463     fhNClustersITS[i]->SetLineColor(color);   fhNClustersITS[i]->SetLineWidth(2);
1464     fhNClustersTPC[i]->SetLineColor(color);   fhNClustersTPC[i]->SetLineWidth(2);
1465     fhNSharedClustersTPC[i]->SetLineColor(color);  fhNSharedClustersTPC[i]->SetLineWidth(2);
1466     fhChi2PerClusterITS[i]->SetLineColor(color);   fhChi2PerClusterITS[i]->SetLineWidth(2);
1467     fhChi2PerClusterTPC[i]->SetLineColor(color);   fhChi2PerClusterTPC[i]->SetLineWidth(2);
1468     fhNClustersForITSPID[i]->SetLineColor(color);  fhNClustersForITSPID[i]->SetLineWidth(2);
1469     fhNMissingITSPoints[i]->SetLineColor(color);   fhNMissingITSPoints[i]->SetLineWidth(2);
1470
1471     fhC11[i]->SetLineColor(color);   fhC11[i]->SetLineWidth(2);
1472     fhC22[i]->SetLineColor(color);   fhC22[i]->SetLineWidth(2);
1473     fhC33[i]->SetLineColor(color);   fhC33[i]->SetLineWidth(2);
1474     fhC44[i]->SetLineColor(color);   fhC44[i]->SetLineWidth(2);
1475     fhC55[i]->SetLineColor(color);   fhC55[i]->SetLineWidth(2);
1476
1477     fhRel1PtUncertainty[i]->SetLineColor(color); fhRel1PtUncertainty[i]->SetLineWidth(2);
1478
1479     fhDXY[i]->SetLineColor(color);   fhDXY[i]->SetLineWidth(2);
1480     fhDZ[i]->SetLineColor(color);    fhDZ[i]->SetLineWidth(2);
1481     fhDXYDZ[i]->SetLineColor(color); fhDXYDZ[i]->SetLineWidth(2);
1482
1483     fhDXYNormalized[i]->SetLineColor(color);   fhDXYNormalized[i]->SetLineWidth(2);
1484     fhDZNormalized[i]->SetLineColor(color);    fhDZNormalized[i]->SetLineWidth(2);
1485     fhNSigmaToVertex[i]->SetLineColor(color);  fhNSigmaToVertex[i]->SetLineWidth(2);
1486   }
1487
1488   // The number of sigmas to the vertex is per definition gaussian
1489   ffDTheoretical = new TF1("nSigmaToVertexTheoretical","([0]/2.506628274)*exp(-(x**2)/2)",0,50);
1490   ffDTheoretical->SetParameter(0,1);
1491
1492   TH1::AddDirectory(oldStatus);
1493 }
1494
1495 //____________________________________________________________________
1496 Bool_t AliESDtrackCuts::LoadHistograms(const Char_t* dir)
1497 {
1498   //
1499   // loads the histograms from a file
1500   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
1501   //
1502
1503   if (!dir)
1504     dir = GetName();
1505
1506   if (!gDirectory->cd(dir))
1507     return kFALSE;
1508
1509   ffDTheoretical = dynamic_cast<TF1*> (gDirectory->Get("nSigmaToVertexTheory"));
1510
1511   fhCutStatistics = dynamic_cast<TH1F*> (gDirectory->Get("cut_statistics"));
1512   fhCutCorrelation = dynamic_cast<TH2F*> (gDirectory->Get("cut_correlation"));
1513
1514   for (Int_t i=0; i<2; i++) {
1515     if (i==0)
1516     {
1517       gDirectory->cd("before_cuts");
1518     }
1519     else
1520       gDirectory->cd("after_cuts");
1521
1522     fhNClustersITS[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersITS"     ));
1523     fhNClustersTPC[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersTPC"     ));
1524     fhNSharedClustersTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get("nSharedClustersTPC"     ));
1525     fhNCrossedRowsTPC[i]   = dynamic_cast<TH1F*> (gDirectory->Get("nCrossedRowsTPC"  ));
1526     fhRatioCrossedRowsOverFindableClustersTPC[i]   = dynamic_cast<TH1F*> (gDirectory->Get("ratioCrossedRowsOverFindableClustersTPC"  ));
1527     fhChi2PerClusterITS[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterITS"));
1528     fhChi2PerClusterTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterTPC"));
1529     fhNClustersForITSPID[i] = dynamic_cast<TH1F*> (gDirectory->Get("nPointsForITSpid"));
1530     fhNMissingITSPoints[i] = dynamic_cast<TH1F*> (gDirectory->Get("nMissingITSClusters"));
1531
1532     fhC11[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal11"));
1533     fhC22[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal22"));
1534     fhC33[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal33"));
1535     fhC44[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal44"));
1536     fhC55[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal55"));
1537
1538     fhRel1PtUncertainty[i] = dynamic_cast<TH1F*> (gDirectory->Get("rel1PtUncertainty"));
1539
1540     fhDXY[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXY"    ));
1541     fhDZ[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZ"     ));
1542     fhDXYDZ[i] =   dynamic_cast<TH1F*> (gDirectory->Get("dXYDZ"));
1543     fhDXYvsDZ[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZ"));
1544
1545     fhDXYNormalized[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXYNormalized"    ));
1546     fhDZNormalized[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZNormalized"     ));
1547     fhDXYvsDZNormalized[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZNormalized"));
1548     fhNSigmaToVertex[i] = dynamic_cast<TH1F*> (gDirectory->Get("nSigmaToVertex"));
1549
1550     fhPt[i] = dynamic_cast<TH1F*> (gDirectory->Get("pt"));
1551     fhEta[i] = dynamic_cast<TH1F*> (gDirectory->Get("eta"));
1552
1553     gDirectory->cd("../");
1554   }
1555
1556   gDirectory->cd("..");
1557
1558   return kTRUE;
1559 }
1560
1561 //____________________________________________________________________
1562 void AliESDtrackCuts::SaveHistograms(const Char_t* dir) {
1563   //
1564   // saves the histograms in a directory (dir)
1565   //
1566
1567   if (!fHistogramsOn) {
1568     AliDebug(0, "Histograms not on - cannot save histograms!!!");
1569     return;
1570   }
1571
1572   if (!dir)
1573     dir = GetName();
1574
1575   gDirectory->mkdir(dir);
1576   gDirectory->cd(dir);
1577
1578   gDirectory->mkdir("before_cuts");
1579   gDirectory->mkdir("after_cuts");
1580
1581   // a factor of 2 is needed since n sigma is positive
1582   ffDTheoretical->SetParameter(0,2*fhNSigmaToVertex[0]->Integral("width"));
1583   ffDTheoretical->Write("nSigmaToVertexTheory");
1584
1585   fhCutStatistics->Write();
1586   fhCutCorrelation->Write();
1587
1588   for (Int_t i=0; i<2; i++) {
1589     if (i==0)
1590       gDirectory->cd("before_cuts");
1591     else
1592       gDirectory->cd("after_cuts");
1593
1594     fhNClustersITS[i]        ->Write();
1595     fhNClustersTPC[i]        ->Write();
1596     fhNSharedClustersTPC[i]  ->Write();
1597     fhNCrossedRowsTPC[i]     ->Write();
1598     fhRatioCrossedRowsOverFindableClustersTPC[i]     ->Write();
1599     fhChi2PerClusterITS[i]   ->Write();
1600     fhChi2PerClusterTPC[i]   ->Write();
1601     fhNClustersForITSPID[i]  ->Write();
1602     fhNMissingITSPoints[i]   ->Write();
1603
1604     fhC11[i]                 ->Write();
1605     fhC22[i]                 ->Write();
1606     fhC33[i]                 ->Write();
1607     fhC44[i]                 ->Write();
1608     fhC55[i]                 ->Write();
1609
1610     fhRel1PtUncertainty[i]   ->Write();
1611
1612     fhDXY[i]                 ->Write();
1613     fhDZ[i]                  ->Write();
1614     fhDXYDZ[i]               ->Write();
1615     fhDXYvsDZ[i]             ->Write();
1616
1617     fhDXYNormalized[i]       ->Write();
1618     fhDZNormalized[i]        ->Write();
1619     fhDXYvsDZNormalized[i]   ->Write();
1620     fhNSigmaToVertex[i]      ->Write();
1621
1622     fhPt[i]                  ->Write();
1623     fhEta[i]                 ->Write();
1624
1625     gDirectory->cd("../");
1626   }
1627
1628   gDirectory->cd("../");
1629 }
1630
1631 //____________________________________________________________________
1632 void AliESDtrackCuts::DrawHistograms()
1633 {
1634   // draws some histograms
1635
1636   TCanvas* canvas1 = new TCanvas(Form("%s_1", GetName()), "Track Quality Results1", 800, 800);
1637   canvas1->Divide(2, 2);
1638
1639   canvas1->cd(1);
1640   fhNClustersTPC[0]->SetStats(kFALSE);
1641   fhNClustersTPC[0]->Draw();
1642
1643   canvas1->cd(2);
1644   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1645   fhChi2PerClusterTPC[0]->Draw();
1646
1647   canvas1->cd(3);
1648   fhNSigmaToVertex[0]->SetStats(kFALSE);
1649   fhNSigmaToVertex[0]->GetXaxis()->SetRangeUser(0, 10);
1650   fhNSigmaToVertex[0]->Draw();
1651
1652   canvas1->SaveAs(Form("%s_%s.gif", GetName(), canvas1->GetName()));
1653
1654   TCanvas* canvas2 = new TCanvas(Form("%s_2", GetName()), "Track Quality Results2", 1200, 800);
1655   canvas2->Divide(3, 2);
1656
1657   canvas2->cd(1);
1658   fhC11[0]->SetStats(kFALSE);
1659   gPad->SetLogy();
1660   fhC11[0]->Draw();
1661
1662   canvas2->cd(2);
1663   fhC22[0]->SetStats(kFALSE);
1664   gPad->SetLogy();
1665   fhC22[0]->Draw();
1666
1667   canvas2->cd(3);
1668   fhC33[0]->SetStats(kFALSE);
1669   gPad->SetLogy();
1670   fhC33[0]->Draw();
1671
1672   canvas2->cd(4);
1673   fhC44[0]->SetStats(kFALSE);
1674   gPad->SetLogy();
1675   fhC44[0]->Draw();
1676
1677   canvas2->cd(5);
1678   fhC55[0]->SetStats(kFALSE);
1679   gPad->SetLogy();
1680   fhC55[0]->Draw();
1681
1682   canvas2->cd(6);
1683   fhRel1PtUncertainty[0]->SetStats(kFALSE);
1684   gPad->SetLogy();
1685   fhRel1PtUncertainty[0]->Draw();
1686
1687   canvas2->SaveAs(Form("%s_%s.gif", GetName(), canvas2->GetName()));
1688
1689   TCanvas* canvas3 = new TCanvas(Form("%s_3", GetName()), "Track Quality Results3", 1200, 800);
1690   canvas3->Divide(3, 2);
1691
1692   canvas3->cd(1);
1693   fhDXY[0]->SetStats(kFALSE);
1694   gPad->SetLogy();
1695   fhDXY[0]->Draw();
1696
1697   canvas3->cd(2);
1698   fhDZ[0]->SetStats(kFALSE);
1699   gPad->SetLogy();
1700   fhDZ[0]->Draw();
1701
1702   canvas3->cd(3);
1703   fhDXYvsDZ[0]->SetStats(kFALSE);
1704   gPad->SetLogz();
1705   gPad->SetRightMargin(0.15);
1706   fhDXYvsDZ[0]->Draw("COLZ");
1707
1708   canvas3->cd(4);
1709   fhDXYNormalized[0]->SetStats(kFALSE);
1710   gPad->SetLogy();
1711   fhDXYNormalized[0]->Draw();
1712
1713   canvas3->cd(5);
1714   fhDZNormalized[0]->SetStats(kFALSE);
1715   gPad->SetLogy();
1716   fhDZNormalized[0]->Draw();
1717
1718   canvas3->cd(6);
1719   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1720   gPad->SetLogz();
1721   gPad->SetRightMargin(0.15);
1722   fhDXYvsDZNormalized[0]->Draw("COLZ");
1723
1724   canvas3->SaveAs(Form("%s_%s.gif", GetName(), canvas3->GetName()));
1725
1726   TCanvas* canvas4 = new TCanvas(Form("%s_4", GetName()), "Track Quality Results4", 800, 500);
1727   canvas4->Divide(2, 1);
1728
1729   canvas4->cd(1);
1730   fhCutStatistics->SetStats(kFALSE);
1731   fhCutStatistics->LabelsOption("v");
1732   gPad->SetBottomMargin(0.3);
1733   fhCutStatistics->Draw();
1734
1735   canvas4->cd(2);
1736   fhCutCorrelation->SetStats(kFALSE);
1737   fhCutCorrelation->LabelsOption("v");
1738   gPad->SetBottomMargin(0.3);
1739   gPad->SetLeftMargin(0.3);
1740   fhCutCorrelation->Draw("COLZ");
1741
1742   canvas4->SaveAs(Form("%s_%s.gif", GetName(), canvas4->GetName()));
1743
1744   /*canvas->cd(1);
1745   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1746   fhDXYvsDZNormalized[0]->DrawCopy("COLZ");
1747
1748   canvas->cd(2);
1749   fhNClustersTPC[0]->SetStats(kFALSE);
1750   fhNClustersTPC[0]->DrawCopy();
1751
1752   canvas->cd(3);
1753   fhChi2PerClusterITS[0]->SetStats(kFALSE);
1754   fhChi2PerClusterITS[0]->DrawCopy();
1755   fhChi2PerClusterITS[1]->SetLineColor(2);
1756   fhChi2PerClusterITS[1]->DrawCopy("SAME");
1757
1758   canvas->cd(4);
1759   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1760   fhChi2PerClusterTPC[0]->DrawCopy();
1761   fhChi2PerClusterTPC[1]->SetLineColor(2);
1762   fhChi2PerClusterTPC[1]->DrawCopy("SAME");*/
1763 }
1764 //--------------------------------------------------------------------------
1765 void AliESDtrackCuts::SetPtDepDCACuts(Double_t pt) {
1766   //
1767   // set the pt-dependent DCA cuts
1768   //
1769
1770   if(f1CutMaxDCAToVertexXYPtDep) {
1771      fCutMaxDCAToVertexXY=f1CutMaxDCAToVertexXYPtDep->Eval(pt);
1772   }
1773
1774   if(f1CutMaxDCAToVertexZPtDep) {
1775     fCutMaxDCAToVertexZ=f1CutMaxDCAToVertexZPtDep->Eval(pt);
1776   }
1777
1778   if(f1CutMinDCAToVertexXYPtDep) {
1779     fCutMinDCAToVertexXY=f1CutMinDCAToVertexXYPtDep->Eval(pt);
1780   }
1781
1782   if(f1CutMinDCAToVertexZPtDep) {
1783     fCutMinDCAToVertexZ=f1CutMinDCAToVertexZPtDep->Eval(pt);
1784   }
1785
1786
1787   return;
1788 }
1789
1790
1791
1792 //--------------------------------------------------------------------------
1793 Bool_t AliESDtrackCuts::CheckPtDepDCA(TString dist,Bool_t print) const {
1794   //
1795   // Check the correctness of the string syntax
1796   //
1797   Bool_t retval=kTRUE;
1798
1799   if(!dist.Contains("pt")) {
1800     if(print) printf("AliESDtrackCuts::CheckPtDepDCA(): string must contain \"pt\"\n");
1801     retval= kFALSE;
1802   }
1803   return retval;
1804 }
1805
1806  void AliESDtrackCuts::SetMaxDCAToVertexXYPtDep(const char *dist){
1807
1808    if(f1CutMaxDCAToVertexXYPtDep){
1809      delete f1CutMaxDCAToVertexXYPtDep;
1810      // resetiing both
1811      f1CutMaxDCAToVertexXYPtDep = 0;
1812      fCutMaxDCAToVertexXYPtDep = "";
1813    }
1814    if(!CheckPtDepDCA(dist,kTRUE)){
1815      return;
1816    }
1817    fCutMaxDCAToVertexXYPtDep = dist;
1818    TString tmp(dist);
1819    tmp.ReplaceAll("pt","x");
1820    f1CutMaxDCAToVertexXYPtDep = new TFormula("f1CutMaxDCAToVertexXYPtDep",tmp.Data());
1821
1822 }
1823
1824  void AliESDtrackCuts::SetMaxDCAToVertexZPtDep(const char *dist){
1825
1826
1827    if(f1CutMaxDCAToVertexZPtDep){
1828      delete f1CutMaxDCAToVertexZPtDep;
1829      // resetiing both
1830      f1CutMaxDCAToVertexZPtDep = 0;
1831      fCutMaxDCAToVertexZPtDep = "";
1832    }
1833    if(!CheckPtDepDCA(dist,kTRUE))return;
1834
1835    fCutMaxDCAToVertexZPtDep = dist;
1836    TString tmp(dist);
1837    tmp.ReplaceAll("pt","x");
1838    f1CutMaxDCAToVertexZPtDep = new TFormula("f1CutMaxDCAToVertexZPtDep",tmp.Data());
1839
1840
1841 }
1842
1843
1844  void AliESDtrackCuts::SetMinDCAToVertexXYPtDep(const char *dist){
1845
1846
1847    if(f1CutMinDCAToVertexXYPtDep){
1848      delete f1CutMinDCAToVertexXYPtDep;
1849      // resetiing both
1850      f1CutMinDCAToVertexXYPtDep = 0;
1851      fCutMinDCAToVertexXYPtDep = "";
1852    }
1853    if(!CheckPtDepDCA(dist,kTRUE))return;
1854
1855    fCutMinDCAToVertexXYPtDep = dist;
1856    TString tmp(dist);
1857    tmp.ReplaceAll("pt","x");
1858    f1CutMinDCAToVertexXYPtDep = new TFormula("f1CutMinDCAToVertexXYPtDep",tmp.Data());
1859
1860 }
1861
1862
1863  void AliESDtrackCuts::SetMinDCAToVertexZPtDep(const char *dist){
1864
1865
1866
1867    if(f1CutMinDCAToVertexZPtDep){
1868      delete f1CutMinDCAToVertexZPtDep;
1869      // resetiing both
1870      f1CutMinDCAToVertexZPtDep = 0;
1871      fCutMinDCAToVertexZPtDep = "";
1872    }
1873    if(!CheckPtDepDCA(dist,kTRUE))return;
1874    fCutMinDCAToVertexZPtDep = dist;
1875    TString tmp(dist);
1876    tmp.ReplaceAll("pt","x");
1877    f1CutMinDCAToVertexZPtDep = new TFormula("f1CutMinDCAToVertexZPtDep",tmp.Data());
1878 }