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