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