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   // creates a TPC only track from the given esd track
 928   // the track has to be deleted by the user
 929   //
 930   // NB. most of the functionality to get a TPC only track from an ESD track is in AliESDtrack, where it should be
 931   // there are only missing propagations here that are needed for old data
 932   // this function will therefore become obsolete
 933   //
 934   // adapted from code provided by CKB
 935
 936   if (!esd->GetPrimaryVertexTPC())
 937     return 0; // No TPC vertex no TPC tracks
 938
 939   if(!esd->GetPrimaryVertexTPC()->GetStatus())
 940     return 0; // TPC Vertex is created by default in AliESDEvent, do not use in this case
 941
 942   AliESDtrack* track = esd->GetTrack(iTrack);
 943   if (!track)
 944     return 0;
 945
 946   AliESDtrack *tpcTrack = new AliESDtrack();
 947
 948   // This should have been done during the reconstruction
 949   // fixed by Juri in r26675
 950   // but recalculate for older data CKB
 951   Float_t p[2],cov[3];
 952   track->GetImpactParametersTPC(p,cov);
 953   if(p[0]==0&&p[1]==0)
 954     track->RelateToVertexTPC(esd->GetPrimaryVertexTPC(),esd->GetMagneticField(),kVeryBig);
 955   // BKC
 956
 957   // only true if we have a tpc track
 958   if (!track->FillTPCOnlyTrack(*tpcTrack))
 959   {
 960     delete tpcTrack;
 961     return 0;
 962   }
 963
 964   // propagate to Vertex
 965   // not needed for normal reconstructed ESDs...
 966   // Double_t pTPC[2],covTPC[3];
 967   // tpcTrack->PropagateToDCA(esd->GetPrimaryVertexTPC(), esd->GetMagneticField(), 10000,  pTPC, covTPC);
 968
 969   return tpcTrack;
 970 }
 971
 972 //____________________________________________________________________
 973 TObjArray* AliESDtrackCuts::GetAcceptedTracks(AliESDEvent* esd,Bool_t bTPC)
 974 {
 975   //
 976   // returns an array of all tracks that pass the cuts
 977   // or an array of TPC only tracks (propagated to the TPC vertex during reco)
 978   // tracks that pass the cut
 979
 980   TObjArray* acceptedTracks = new TObjArray();
 981
 982   // loop over esd tracks
 983   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
 984     if(bTPC){
 985       if(!esd->GetPrimaryVertexTPC())return acceptedTracks; // No TPC vertex no TPC tracks
 986       if(!esd->GetPrimaryVertexTPC()->GetStatus())return acceptedTracks; // No proper TPC vertex, only the default
 987
 988       AliESDtrack *tpcTrack = GetTPCOnlyTrack(esd, iTrack);
 989       if (!tpcTrack)
 990         continue;
 991
 992       if (AcceptTrack(tpcTrack)) {
 993         acceptedTracks->Add(tpcTrack);
 994       }
 995       else
 996         delete tpcTrack;
 997     }
 998     else
 999     {
1000       AliESDtrack* track = esd->GetTrack(iTrack);
1001       if(AcceptTrack(track))
1002         acceptedTracks->Add(track);
1003     }
1004   }
1005   if(bTPC)acceptedTracks->SetOwner(kTRUE);
1006   return acceptedTracks;
1007 }
1008
1009 //____________________________________________________________________
1010 Int_t AliESDtrackCuts::CountAcceptedTracks(AliESDEvent* esd)
1011 {
1012   //
1013   // returns an the number of tracks that pass the cuts
1014   //
1015
1016   Int_t count = 0;
1017
1018   // loop over esd tracks
1019   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
1020     AliESDtrack* track = esd->GetTrack(iTrack);
1021     if (AcceptTrack(track))
1022       count++;
1023   }
1024
1025   return count;
1026 }
1027
1028 //____________________________________________________________________
1029  void AliESDtrackCuts::DefineHistograms(Int_t color) {
1030    //
1031    // diagnostics histograms are defined
1032    //
1033
1034    fHistogramsOn=kTRUE;
1035
1036    Bool_t oldStatus = TH1::AddDirectoryStatus();
1037    TH1::AddDirectory(kFALSE);
1038
1039    //###################################################################################
1040    // defining histograms
1041
1042    fhCutStatistics = new TH1F("cut_statistics","cut statistics",kNCuts+4,-0.5,kNCuts+3.5);
1043
1044    fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
1045    fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
1046
1047    fhCutCorrelation = new TH2F("cut_correlation","cut correlation",kNCuts,-0.5,kNCuts-0.5,kNCuts,-0.5,kNCuts-0.5);;
1048
1049    for (Int_t i=0; i<kNCuts; i++) {
1050      fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fgkCutNames[i]);
1051      fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fgkCutNames[i]);
1052      fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fgkCutNames[i]);
1053    }
1054
1055   fhCutStatistics  ->SetLineColor(color);
1056   fhCutCorrelation ->SetLineColor(color);
1057   fhCutStatistics  ->SetLineWidth(2);
1058   fhCutCorrelation ->SetLineWidth(2);
1059
1060   for (Int_t i=0; i<2; i++) {
1061     fhNClustersITS[i]        = new TH1F("nClustersITS"    ,"",8,-0.5,7.5);
1062     fhNClustersTPC[i]        = new TH1F("nClustersTPC"     ,"",165,-0.5,164.5);
1063     fhChi2PerClusterITS[i]   = new TH1F("chi2PerClusterITS","",500,0,10);
1064     fhChi2PerClusterTPC[i]   = new TH1F("chi2PerClusterTPC","",500,0,10);
1065
1066     fhC11[i]                 = new TH1F("covMatrixDiagonal11","",2000,0,20);
1067     fhC22[i]                 = new TH1F("covMatrixDiagonal22","",2000,0,20);
1068     fhC33[i]                 = new TH1F("covMatrixDiagonal33","",1000,0,0.1);
1069     fhC44[i]                 = new TH1F("covMatrixDiagonal44","",1000,0,0.1);
1070     fhC55[i]                 = new TH1F("covMatrixDiagonal55","",1000,0,5);
1071
1072     fhRel1PtUncertainty[i]   = new TH1F("rel1PtUncertainty","",1000,0,5);
1073
1074     fhDXY[i]                 = new TH1F("dXY"    ,"",500,-10,10);
1075     fhDZ[i]                  = new TH1F("dZ"     ,"",500,-10,10);
1076     fhDXYDZ[i]               = new TH1F("dXYDZ"  ,"",500,0,10);
1077     fhDXYvsDZ[i]             = new TH2F("dXYvsDZ","",200,-10,10,200,-10,10);
1078
1079     fhDXYNormalized[i]       = new TH1F("dXYNormalized"    ,"",500,-10,10);
1080     fhDZNormalized[i]        = new TH1F("dZNormalized"     ,"",500,-10,10);
1081     fhDXYvsDZNormalized[i]   = new TH2F("dXYvsDZNormalized","",200,-10,10,200,-10,10);
1082
1083     fhNSigmaToVertex[i]      = new TH1F("nSigmaToVertex","",500,0,10);
1084
1085     fhPt[i]                  = new TH1F("pt"     ,"p_{T} distribution;p_{T} (GeV/c)", 800, 0.0, 10.0);
1086     fhEta[i]                 = new TH1F("eta"     ,"#eta distribution;#eta",40,-2.0,2.0);
1087
1088     fhNClustersITS[i]->SetTitle("n ITS clusters");
1089     fhNClustersTPC[i]->SetTitle("n TPC clusters");
1090     fhChi2PerClusterITS[i]->SetTitle("#Chi^{2} per ITS cluster");
1091     fhChi2PerClusterTPC[i]->SetTitle("#Chi^{2} per TPC cluster");
1092
1093     fhC11[i]->SetTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
1094     fhC22[i]->SetTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
1095     fhC33[i]->SetTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
1096     fhC44[i]->SetTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
1097     fhC55[i]->SetTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
1098
1099     fhRel1PtUncertainty[i]->SetTitle("rel. uncertainty of 1/p_{T}");
1100
1101     fhDXY[i]->SetXTitle("transverse impact parameter (cm)");
1102     fhDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
1103     fhDXYDZ[i]->SetTitle("absolute impact parameter;sqrt(dXY**2 + dZ**2)  (cm)");
1104     fhDXYvsDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
1105     fhDXYvsDZ[i]->SetYTitle("transverse impact parameter (cm)");
1106
1107     fhDXYNormalized[i]->SetTitle("normalized trans impact par (n#sigma)");
1108     fhDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
1109     fhDXYvsDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
1110     fhDXYvsDZNormalized[i]->SetYTitle("normalized trans impact par (n#sigma)");
1111     fhNSigmaToVertex[i]->SetTitle("n #sigma to vertex");
1112
1113     fhNClustersITS[i]->SetLineColor(color);   fhNClustersITS[i]->SetLineWidth(2);
1114     fhNClustersTPC[i]->SetLineColor(color);   fhNClustersTPC[i]->SetLineWidth(2);
1115     fhChi2PerClusterITS[i]->SetLineColor(color);   fhChi2PerClusterITS[i]->SetLineWidth(2);
1116     fhChi2PerClusterTPC[i]->SetLineColor(color);   fhChi2PerClusterTPC[i]->SetLineWidth(2);
1117
1118     fhC11[i]->SetLineColor(color);   fhC11[i]->SetLineWidth(2);
1119     fhC22[i]->SetLineColor(color);   fhC22[i]->SetLineWidth(2);
1120     fhC33[i]->SetLineColor(color);   fhC33[i]->SetLineWidth(2);
1121     fhC44[i]->SetLineColor(color);   fhC44[i]->SetLineWidth(2);
1122     fhC55[i]->SetLineColor(color);   fhC55[i]->SetLineWidth(2);
1123
1124     fhRel1PtUncertainty[i]->SetLineColor(color); fhRel1PtUncertainty[i]->SetLineWidth(2);
1125
1126     fhDXY[i]->SetLineColor(color);   fhDXY[i]->SetLineWidth(2);
1127     fhDZ[i]->SetLineColor(color);    fhDZ[i]->SetLineWidth(2);
1128     fhDXYDZ[i]->SetLineColor(color); fhDXYDZ[i]->SetLineWidth(2);
1129
1130     fhDXYNormalized[i]->SetLineColor(color);   fhDXYNormalized[i]->SetLineWidth(2);
1131     fhDZNormalized[i]->SetLineColor(color);    fhDZNormalized[i]->SetLineWidth(2);
1132     fhNSigmaToVertex[i]->SetLineColor(color);  fhNSigmaToVertex[i]->SetLineWidth(2);
1133   }
1134
1135   // The number of sigmas to the vertex is per definition gaussian
1136   ffDTheoretical = new TF1("nSigmaToVertexTheoretical","([0]/2.506628274)*exp(-(x**2)/2)",0,50);
1137   ffDTheoretical->SetParameter(0,1);
1138
1139   TH1::AddDirectory(oldStatus);
1140 }
1141
1142 //____________________________________________________________________
1143 Bool_t AliESDtrackCuts::LoadHistograms(const Char_t* dir)
1144 {
1145   //
1146   // loads the histograms from a file
1147   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
1148   //
1149
1150   if (!dir)
1151     dir = GetName();
1152
1153   if (!gDirectory->cd(dir))
1154     return kFALSE;
1155
1156   ffDTheoretical = dynamic_cast<TF1*> (gDirectory->Get("nSigmaToVertexTheory"));
1157
1158   fhCutStatistics = dynamic_cast<TH1F*> (gDirectory->Get("cut_statistics"));
1159   fhCutCorrelation = dynamic_cast<TH2F*> (gDirectory->Get("cut_correlation"));
1160
1161   for (Int_t i=0; i<2; i++) {
1162     if (i==0)
1163     {
1164       gDirectory->cd("before_cuts");
1165     }
1166     else
1167       gDirectory->cd("after_cuts");
1168
1169     fhNClustersITS[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersITS"     ));
1170     fhNClustersTPC[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersTPC"     ));
1171     fhChi2PerClusterITS[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterITS"));
1172     fhChi2PerClusterTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterTPC"));
1173
1174     fhC11[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal11"));
1175     fhC22[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal22"));
1176     fhC33[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal33"));
1177     fhC44[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal44"));
1178     fhC55[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal55"));
1179
1180     fhRel1PtUncertainty[i] = dynamic_cast<TH1F*> (gDirectory->Get("rel1PtUncertainty"));
1181
1182     fhDXY[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXY"    ));
1183     fhDZ[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZ"     ));
1184     fhDXYDZ[i] =   dynamic_cast<TH1F*> (gDirectory->Get("dXYDZ"));
1185     fhDXYvsDZ[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZ"));
1186
1187     fhDXYNormalized[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXYNormalized"    ));
1188     fhDZNormalized[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZNormalized"     ));
1189     fhDXYvsDZNormalized[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZNormalized"));
1190     fhNSigmaToVertex[i] = dynamic_cast<TH1F*> (gDirectory->Get("nSigmaToVertex"));
1191
1192     fhPt[i] = dynamic_cast<TH1F*> (gDirectory->Get("pt"));
1193     fhEta[i] = dynamic_cast<TH1F*> (gDirectory->Get("eta"));
1194
1195     gDirectory->cd("../");
1196   }
1197
1198   gDirectory->cd("..");
1199
1200   return kTRUE;
1201 }
1202
1203 //____________________________________________________________________
1204 void AliESDtrackCuts::SaveHistograms(const Char_t* dir) {
1205   //
1206   // saves the histograms in a directory (dir)
1207   //
1208
1209   if (!fHistogramsOn) {
1210     AliDebug(0, "Histograms not on - cannot save histograms!!!");
1211     return;
1212   }
1213
1214   if (!dir)
1215     dir = GetName();
1216
1217   gDirectory->mkdir(dir);
1218   gDirectory->cd(dir);
1219
1220   gDirectory->mkdir("before_cuts");
1221   gDirectory->mkdir("after_cuts");
1222
1223   // a factor of 2 is needed since n sigma is positive
1224   ffDTheoretical->SetParameter(0,2*fhNSigmaToVertex[0]->Integral("width"));
1225   ffDTheoretical->Write("nSigmaToVertexTheory");
1226
1227   fhCutStatistics->Write();
1228   fhCutCorrelation->Write();
1229
1230   for (Int_t i=0; i<2; i++) {
1231     if (i==0)
1232       gDirectory->cd("before_cuts");
1233     else
1234       gDirectory->cd("after_cuts");
1235
1236     fhNClustersITS[i]        ->Write();
1237     fhNClustersTPC[i]        ->Write();
1238     fhChi2PerClusterITS[i]   ->Write();
1239     fhChi2PerClusterTPC[i]   ->Write();
1240
1241     fhC11[i]                 ->Write();
1242     fhC22[i]                 ->Write();
1243     fhC33[i]                 ->Write();
1244     fhC44[i]                 ->Write();
1245     fhC55[i]                 ->Write();
1246
1247     fhRel1PtUncertainty[i]   ->Write();
1248
1249     fhDXY[i]                 ->Write();
1250     fhDZ[i]                  ->Write();
1251     fhDXYDZ[i]               ->Write();
1252     fhDXYvsDZ[i]             ->Write();
1253
1254     fhDXYNormalized[i]       ->Write();
1255     fhDZNormalized[i]        ->Write();
1256     fhDXYvsDZNormalized[i]   ->Write();
1257     fhNSigmaToVertex[i]      ->Write();
1258
1259     fhPt[i]                  ->Write();
1260     fhEta[i]                 ->Write();
1261
1262     gDirectory->cd("../");
1263   }
1264
1265   gDirectory->cd("../");
1266 }
1267
1268 //____________________________________________________________________
1269 void AliESDtrackCuts::DrawHistograms()
1270 {
1271   // draws some histograms
1272
1273   TCanvas* canvas1 = new TCanvas(Form("%s_1", GetName()), "Track Quality Results1", 800, 800);
1274   canvas1->Divide(2, 2);
1275
1276   canvas1->cd(1);
1277   fhNClustersTPC[0]->SetStats(kFALSE);
1278   fhNClustersTPC[0]->Draw();
1279
1280   canvas1->cd(2);
1281   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1282   fhChi2PerClusterTPC[0]->Draw();
1283
1284   canvas1->cd(3);
1285   fhNSigmaToVertex[0]->SetStats(kFALSE);
1286   fhNSigmaToVertex[0]->GetXaxis()->SetRangeUser(0, 10);
1287   fhNSigmaToVertex[0]->Draw();
1288
1289   canvas1->SaveAs(Form("%s_%s.gif", GetName(), canvas1->GetName()));
1290
1291   TCanvas* canvas2 = new TCanvas(Form("%s_2", GetName()), "Track Quality Results2", 1200, 800);
1292   canvas2->Divide(3, 2);
1293
1294   canvas2->cd(1);
1295   fhC11[0]->SetStats(kFALSE);
1296   gPad->SetLogy();
1297   fhC11[0]->Draw();
1298
1299   canvas2->cd(2);
1300   fhC22[0]->SetStats(kFALSE);
1301   gPad->SetLogy();
1302   fhC22[0]->Draw();
1303
1304   canvas2->cd(3);
1305   fhC33[0]->SetStats(kFALSE);
1306   gPad->SetLogy();
1307   fhC33[0]->Draw();
1308
1309   canvas2->cd(4);
1310   fhC44[0]->SetStats(kFALSE);
1311   gPad->SetLogy();
1312   fhC44[0]->Draw();
1313
1314   canvas2->cd(5);
1315   fhC55[0]->SetStats(kFALSE);
1316   gPad->SetLogy();
1317   fhC55[0]->Draw();
1318
1319   canvas2->cd(6);
1320   fhRel1PtUncertainty[0]->SetStats(kFALSE);
1321   gPad->SetLogy();
1322   fhRel1PtUncertainty[0]->Draw();
1323
1324   canvas2->SaveAs(Form("%s_%s.gif", GetName(), canvas2->GetName()));
1325
1326   TCanvas* canvas3 = new TCanvas(Form("%s_3", GetName()), "Track Quality Results3", 1200, 800);
1327   canvas3->Divide(3, 2);
1328
1329   canvas3->cd(1);
1330   fhDXY[0]->SetStats(kFALSE);
1331   gPad->SetLogy();
1332   fhDXY[0]->Draw();
1333
1334   canvas3->cd(2);
1335   fhDZ[0]->SetStats(kFALSE);
1336   gPad->SetLogy();
1337   fhDZ[0]->Draw();
1338
1339   canvas3->cd(3);
1340   fhDXYvsDZ[0]->SetStats(kFALSE);
1341   gPad->SetLogz();
1342   gPad->SetRightMargin(0.15);
1343   fhDXYvsDZ[0]->Draw("COLZ");
1344
1345   canvas3->cd(4);
1346   fhDXYNormalized[0]->SetStats(kFALSE);
1347   gPad->SetLogy();
1348   fhDXYNormalized[0]->Draw();
1349
1350   canvas3->cd(5);
1351   fhDZNormalized[0]->SetStats(kFALSE);
1352   gPad->SetLogy();
1353   fhDZNormalized[0]->Draw();
1354
1355   canvas3->cd(6);
1356   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1357   gPad->SetLogz();
1358   gPad->SetRightMargin(0.15);
1359   fhDXYvsDZNormalized[0]->Draw("COLZ");
1360
1361   canvas3->SaveAs(Form("%s_%s.gif", GetName(), canvas3->GetName()));
1362
1363   TCanvas* canvas4 = new TCanvas(Form("%s_4", GetName()), "Track Quality Results4", 800, 500);
1364   canvas4->Divide(2, 1);
1365
1366   canvas4->cd(1);
1367   fhCutStatistics->SetStats(kFALSE);
1368   fhCutStatistics->LabelsOption("v");
1369   gPad->SetBottomMargin(0.3);
1370   fhCutStatistics->Draw();
1371
1372   canvas4->cd(2);
1373   fhCutCorrelation->SetStats(kFALSE);
1374   fhCutCorrelation->LabelsOption("v");
1375   gPad->SetBottomMargin(0.3);
1376   gPad->SetLeftMargin(0.3);
1377   fhCutCorrelation->Draw("COLZ");
1378
1379   canvas4->SaveAs(Form("%s_%s.gif", GetName(), canvas4->GetName()));
1380
1381   /*canvas->cd(1);
1382   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1383   fhDXYvsDZNormalized[0]->DrawCopy("COLZ");
1384
1385   canvas->cd(2);
1386   fhNClustersTPC[0]->SetStats(kFALSE);
1387   fhNClustersTPC[0]->DrawCopy();
1388
1389   canvas->cd(3);
1390   fhChi2PerClusterITS[0]->SetStats(kFALSE);
1391   fhChi2PerClusterITS[0]->DrawCopy();
1392   fhChi2PerClusterITS[1]->SetLineColor(2);
1393   fhChi2PerClusterITS[1]->DrawCopy("SAME");
1394
1395   canvas->cd(4);
1396   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1397   fhChi2PerClusterTPC[0]->DrawCopy();
1398   fhChi2PerClusterTPC[1]->SetLineColor(2);
1399   fhChi2PerClusterTPC[1]->DrawCopy("SAME");*/
1400 }
1401