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