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
 702   //y-eta related calculations
 703   Float_t eta = -100.;
 704   Float_t y   = -100.;
 705   if((momentum != TMath::Abs(p[2]))&&(momentum != 0))
 706     eta = 0.5*TMath::Log((momentum + p[2])/(momentum - p[2]));
 707   if((energy != TMath::Abs(p[2]))&&(momentum != 0))
 708     y = 0.5*TMath::Log((energy + p[2])/(energy - p[2]));
 709
 710   Float_t relUncertainty1Pt = TMath::Sqrt(extCov[14])*pt;
 711
 712   //########################################################################
 713   // cut the track?
 714
 715   Bool_t cuts[kNCuts];
 716   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
 717
 718   // track quality cuts
 719   if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
 720     cuts[0]=kTRUE;
 721   if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
 722     cuts[1]=kTRUE;
 723   if (nClustersTPC<fCutMinNClusterTPC)
 724     cuts[2]=kTRUE;
 725   if (nClustersITS<fCutMinNClusterITS)
 726     cuts[3]=kTRUE;
 727   if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
 728     cuts[4]=kTRUE;
 729   if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
 730     cuts[5]=kTRUE;
 731   if (extCov[0]  > fCutMaxC11)
 732     cuts[6]=kTRUE;
 733   if (extCov[2]  > fCutMaxC22)
 734     cuts[7]=kTRUE;
 735   if (extCov[5]  > fCutMaxC33)
 736     cuts[8]=kTRUE;
 737   if (extCov[9]  > fCutMaxC44)
 738     cuts[9]=kTRUE;
 739   if (extCov[14]  > fCutMaxC55)
 740     cuts[10]=kTRUE;
 741   if (nSigmaToVertex > fCutNsigmaToVertex && fCutSigmaToVertexRequired)
 742     cuts[11] = kTRUE;
 743   // if n sigma could not be calculated
 744   if (nSigmaToVertex<0 && fCutSigmaToVertexRequired)
 745     cuts[12]=kTRUE;
 746   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
 747     cuts[13]=kTRUE;
 748   // track kinematics cut
 749   if((momentum < fPMin) || (momentum > fPMax))
 750     cuts[14]=kTRUE;
 751   if((pt < fPtMin) || (pt > fPtMax))
 752     cuts[15] = kTRUE;
 753   if((p[0] < fPxMin) || (p[0] > fPxMax))
 754     cuts[16] = kTRUE;
 755   if((p[1] < fPyMin) || (p[1] > fPyMax))
 756     cuts[17] = kTRUE;
 757   if((p[2] < fPzMin) || (p[2] > fPzMax))
 758     cuts[18] = kTRUE;
 759   if((eta < fEtaMin) || (eta > fEtaMax))
 760     cuts[19] = kTRUE;
 761   if((y < fRapMin) || (y > fRapMax))
 762     cuts[20] = kTRUE;
 763   if (fCutDCAToVertex2D && dcaToVertex > 1)
 764     cuts[21] = kTRUE;
 765   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
 766     cuts[22] = kTRUE;
 767   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
 768     cuts[23] = kTRUE;
 769   if (fCutDCAToVertex2D && fCutMinDCAToVertexXY > 0 && fCutMinDCAToVertexZ > 0 && dcaToVertexXY*dcaToVertexXY/fCutMinDCAToVertexXY/fCutMinDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMinDCAToVertexZ/fCutMinDCAToVertexZ < 1)
 770     cuts[24] = kTRUE;
 771   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) < fCutMinDCAToVertexXY)
 772     cuts[25] = kTRUE;
 773   if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) < fCutMinDCAToVertexZ)
 774     cuts[26] = kTRUE;
 775
 776   for (Int_t i = 0; i < 3; i++)
 777     cuts[27+i] = !CheckITSClusterRequirement(fCutClusterRequirementITS[i], esdTrack->HasPointOnITSLayer(i*2), esdTrack->HasPointOnITSLayer(i*2+1));
 778
 779   if (fCutRequireITSStandAlone && ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin)))
 780     cuts[30]=kTRUE;
 781
 782   if (relUncertainty1Pt > fCutMaxRel1PtUncertainty)
 783      cuts[31]=kTRUE;
 784
 785   Bool_t cut=kFALSE;
 786   for (Int_t i=0; i<kNCuts; i++)
 787     if (cuts[i]) {cut = kTRUE;}
 788
 789
 790
 791   //########################################################################
 792   // filling histograms
 793   if (fHistogramsOn) {
 794     fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n tracks")));
 795     if (cut)
 796       fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut tracks")));
 797
 798     for (Int_t i=0; i<kNCuts; i++) {
 799       if (cuts[i])
 800         fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
 801
 802       for (Int_t j=i; j<kNCuts; j++) {
 803         if (cuts[i] && cuts[j]) {
 804           Float_t xC = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
 805           Float_t yC = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
 806           fhCutCorrelation->Fill(xC, yC);
 807         }
 808       }
 809     }
 810   }
 811
 812   // now we loop over the filling of the histograms twice: once "before" the cut, once "after"
 813   // the code is not in a function due to too many local variables that would need to be passed
 814
 815   for (Int_t id = 0; id < 2; id++)
 816   {
 817     // id = 0 --> before cut
 818     // id = 1 --> after cut
 819
 820     if (fHistogramsOn)
 821     {
 822       fhNClustersITS[id]->Fill(nClustersITS);
 823       fhNClustersTPC[id]->Fill(nClustersTPC);
 824       fhChi2PerClusterITS[id]->Fill(chi2PerClusterITS);
 825       fhChi2PerClusterTPC[id]->Fill(chi2PerClusterTPC);
 826
 827       fhC11[id]->Fill(extCov[0]);
 828       fhC22[id]->Fill(extCov[2]);
 829       fhC33[id]->Fill(extCov[5]);
 830       fhC44[id]->Fill(extCov[9]);
 831       fhC55[id]->Fill(extCov[14]);
 832
 833       fhRel1PtUncertainty[id]->Fill(relUncertainty1Pt);
 834
 835       fhPt[id]->Fill(pt);
 836       fhEta[id]->Fill(eta);
 837
 838       Float_t bRes[2];
 839       bRes[0] = TMath::Sqrt(bCov[0]);
 840       bRes[1] = TMath::Sqrt(bCov[2]);
 841
 842       fhDZ[id]->Fill(b[1]);
 843       fhDXY[id]->Fill(b[0]);
 844       fhDXYDZ[id]->Fill(dcaToVertex);
 845       fhDXYvsDZ[id]->Fill(b[1],b[0]);
 846
 847       if (bRes[0]!=0 && bRes[1]!=0) {
 848         fhDZNormalized[id]->Fill(b[1]/bRes[1]);
 849         fhDXYNormalized[id]->Fill(b[0]/bRes[0]);
 850         fhDXYvsDZNormalized[id]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
 851         fhNSigmaToVertex[id]->Fill(nSigmaToVertex);
 852       }
 853     }
 854
 855     // cut the track
 856     if (cut)
 857       return kFALSE;
 858   }
 859
 860   return kTRUE;
 861 }
 862
 863 //____________________________________________________________________
 864 Bool_t AliESDtrackCuts::CheckITSClusterRequirement(ITSClusterRequirement req, Bool_t clusterL1, Bool_t clusterL2)
 865 {
 866   // checks if the cluster requirement is fullfilled (in this case: return kTRUE)
 867
 868   switch (req)
 869   {
 870         case kOff:        return kTRUE;
 871         case kNone:       return !clusterL1 && !clusterL2;
 872         case kAny:        return clusterL1 || clusterL2;
 873         case kFirst:      return clusterL1;
 874         case kOnlyFirst:  return clusterL1 && !clusterL2;
 875         case kSecond:     return clusterL2;
 876         case kOnlySecond: return clusterL2 && !clusterL1;
 877         case kBoth:       return clusterL1 && clusterL2;
 878   }
 879
 880   return kFALSE;
 881 }
 882
 883 //____________________________________________________________________
 884 AliESDtrack* AliESDtrackCuts::GetTPCOnlyTrack(AliESDEvent* esd, Int_t iTrack)
 885 {
 886   // creates a TPC only track from the given esd track
 887   // the track has to be deleted by the user
 888   //
 889   // NB. most of the functionality to get a TPC only track from an ESD track is in AliESDtrack, where it should be
 890   // there are only missing propagations here that are needed for old data
 891   // this function will therefore become obsolete
 892   //
 893   // adapted from code provided by CKB
 894
 895   if (!esd->GetPrimaryVertexTPC())
 896     return 0; // No TPC vertex no TPC tracks
 897
 898   if(!esd->GetPrimaryVertexTPC()->GetStatus())
 899     return 0; // TPC Vertex is created by default in AliESDEvent, do not use in this case
 900
 901   AliESDtrack* track = esd->GetTrack(iTrack);
 902   if (!track)
 903     return 0;
 904
 905   AliESDtrack *tpcTrack = new AliESDtrack();
 906
 907   // This should have been done during the reconstruction
 908   // fixed by Juri in r26675
 909   // but recalculate for older data CKB
 910   Float_t p[2],cov[3];
 911   track->GetImpactParametersTPC(p,cov);
 912   if(p[0]==0&&p[1]==0)
 913     track->RelateToVertexTPC(esd->GetPrimaryVertexTPC(),esd->GetMagneticField(),kVeryBig);
 914   // BKC
 915
 916   // only true if we have a tpc track
 917   if (!track->FillTPCOnlyTrack(*tpcTrack))
 918   {
 919     delete tpcTrack;
 920     return 0;
 921   }
 922
 923   // propagate to Vertex
 924   // not needed for normal reconstructed ESDs...
 925   // Double_t pTPC[2],covTPC[3];
 926   // tpcTrack->PropagateToDCA(esd->GetPrimaryVertexTPC(), esd->GetMagneticField(), 10000,  pTPC, covTPC);
 927
 928   return tpcTrack;
 929 }
 930
 931 //____________________________________________________________________
 932 TObjArray* AliESDtrackCuts::GetAcceptedTracks(AliESDEvent* esd,Bool_t bTPC)
 933 {
 934   //
 935   // returns an array of all tracks that pass the cuts
 936   // or an array of TPC only tracks (propagated to the TPC vertex during reco)
 937   // tracks that pass the cut
 938
 939   TObjArray* acceptedTracks = new TObjArray();
 940
 941   // loop over esd tracks
 942   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
 943     if(bTPC){
 944       if(!esd->GetPrimaryVertexTPC())return acceptedTracks; // No TPC vertex no TPC tracks
 945       if(!esd->GetPrimaryVertexTPC()->GetStatus())return acceptedTracks; // No proper TPC vertex, only the default
 946
 947       AliESDtrack *tpcTrack = GetTPCOnlyTrack(esd, iTrack);
 948       if (!tpcTrack)
 949         continue;
 950
 951       if (AcceptTrack(tpcTrack)) {
 952         acceptedTracks->Add(tpcTrack);
 953       }
 954       else
 955         delete tpcTrack;
 956     }
 957     else
 958     {
 959       AliESDtrack* track = esd->GetTrack(iTrack);
 960       if(AcceptTrack(track))
 961         acceptedTracks->Add(track);
 962     }
 963   }
 964   if(bTPC)acceptedTracks->SetOwner(kTRUE);
 965   return acceptedTracks;
 966 }
 967
 968 //____________________________________________________________________
 969 Int_t AliESDtrackCuts::CountAcceptedTracks(AliESDEvent* esd)
 970 {
 971   //
 972   // returns an the number of tracks that pass the cuts
 973   //
 974
 975   Int_t count = 0;
 976
 977   // loop over esd tracks
 978   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
 979     AliESDtrack* track = esd->GetTrack(iTrack);
 980     if (AcceptTrack(track))
 981       count++;
 982   }
 983
 984   return count;
 985 }
 986
 987 //____________________________________________________________________
 988  void AliESDtrackCuts::DefineHistograms(Int_t color) {
 989    //
 990    // diagnostics histograms are defined
 991    //
 992
 993    fHistogramsOn=kTRUE;
 994
 995    Bool_t oldStatus = TH1::AddDirectoryStatus();
 996    TH1::AddDirectory(kFALSE);
 997
 998    //###################################################################################
 999    // defining histograms
1000
1001    fhCutStatistics = new TH1F("cut_statistics","cut statistics",kNCuts+4,-0.5,kNCuts+3.5);
1002
1003    fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
1004    fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
1005
1006    fhCutCorrelation = new TH2F("cut_correlation","cut correlation",kNCuts,-0.5,kNCuts-0.5,kNCuts,-0.5,kNCuts-0.5);;
1007
1008    for (Int_t i=0; i<kNCuts; i++) {
1009      fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fgkCutNames[i]);
1010      fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fgkCutNames[i]);
1011      fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fgkCutNames[i]);
1012    }
1013
1014   fhCutStatistics  ->SetLineColor(color);
1015   fhCutCorrelation ->SetLineColor(color);
1016   fhCutStatistics  ->SetLineWidth(2);
1017   fhCutCorrelation ->SetLineWidth(2);
1018
1019   for (Int_t i=0; i<2; i++) {
1020     fhNClustersITS[i]        = new TH1F("nClustersITS"    ,"",8,-0.5,7.5);
1021     fhNClustersTPC[i]        = new TH1F("nClustersTPC"     ,"",165,-0.5,164.5);
1022     fhChi2PerClusterITS[i]   = new TH1F("chi2PerClusterITS","",500,0,10);
1023     fhChi2PerClusterTPC[i]   = new TH1F("chi2PerClusterTPC","",500,0,10);
1024
1025     fhC11[i]                 = new TH1F("covMatrixDiagonal11","",2000,0,20);
1026     fhC22[i]                 = new TH1F("covMatrixDiagonal22","",2000,0,20);
1027     fhC33[i]                 = new TH1F("covMatrixDiagonal33","",1000,0,0.1);
1028     fhC44[i]                 = new TH1F("covMatrixDiagonal44","",1000,0,0.1);
1029     fhC55[i]                 = new TH1F("covMatrixDiagonal55","",1000,0,5);
1030
1031     fhRel1PtUncertainty[i]   = new TH1F("rel1PtUncertainty","",1000,0,5);
1032
1033     fhDXY[i]                 = new TH1F("dXY"    ,"",500,-10,10);
1034     fhDZ[i]                  = new TH1F("dZ"     ,"",500,-10,10);
1035     fhDXYDZ[i]               = new TH1F("dXYDZ"  ,"",500,0,10);
1036     fhDXYvsDZ[i]             = new TH2F("dXYvsDZ","",200,-10,10,200,-10,10);
1037
1038     fhDXYNormalized[i]       = new TH1F("dXYNormalized"    ,"",500,-10,10);
1039     fhDZNormalized[i]        = new TH1F("dZNormalized"     ,"",500,-10,10);
1040     fhDXYvsDZNormalized[i]   = new TH2F("dXYvsDZNormalized","",200,-10,10,200,-10,10);
1041
1042     fhNSigmaToVertex[i]      = new TH1F("nSigmaToVertex","",500,0,10);
1043
1044     fhPt[i]                  = new TH1F("pt"     ,"p_{T} distribution;p_{T} (GeV/c)", 800, 0.0, 10.0);
1045     fhEta[i]                 = new TH1F("eta"     ,"#eta distribution;#eta",40,-2.0,2.0);
1046
1047     fhNClustersITS[i]->SetTitle("n ITS clusters");
1048     fhNClustersTPC[i]->SetTitle("n TPC clusters");
1049     fhChi2PerClusterITS[i]->SetTitle("#Chi^{2} per ITS cluster");
1050     fhChi2PerClusterTPC[i]->SetTitle("#Chi^{2} per TPC cluster");
1051
1052     fhC11[i]->SetTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
1053     fhC22[i]->SetTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
1054     fhC33[i]->SetTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
1055     fhC44[i]->SetTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
1056     fhC55[i]->SetTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
1057
1058     fhRel1PtUncertainty[i]->SetTitle("rel. uncertainty of 1/p_{T}");
1059
1060     fhDXY[i]->SetXTitle("transverse impact parameter (cm)");
1061     fhDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
1062     fhDXYDZ[i]->SetTitle("absolute impact parameter;sqrt(dXY**2 + dZ**2)  (cm)");
1063     fhDXYvsDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
1064     fhDXYvsDZ[i]->SetYTitle("transverse impact parameter (cm)");
1065
1066     fhDXYNormalized[i]->SetTitle("normalized trans impact par (n#sigma)");
1067     fhDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
1068     fhDXYvsDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
1069     fhDXYvsDZNormalized[i]->SetYTitle("normalized trans impact par (n#sigma)");
1070     fhNSigmaToVertex[i]->SetTitle("n #sigma to vertex");
1071
1072     fhNClustersITS[i]->SetLineColor(color);   fhNClustersITS[i]->SetLineWidth(2);
1073     fhNClustersTPC[i]->SetLineColor(color);   fhNClustersTPC[i]->SetLineWidth(2);
1074     fhChi2PerClusterITS[i]->SetLineColor(color);   fhChi2PerClusterITS[i]->SetLineWidth(2);
1075     fhChi2PerClusterTPC[i]->SetLineColor(color);   fhChi2PerClusterTPC[i]->SetLineWidth(2);
1076
1077     fhC11[i]->SetLineColor(color);   fhC11[i]->SetLineWidth(2);
1078     fhC22[i]->SetLineColor(color);   fhC22[i]->SetLineWidth(2);
1079     fhC33[i]->SetLineColor(color);   fhC33[i]->SetLineWidth(2);
1080     fhC44[i]->SetLineColor(color);   fhC44[i]->SetLineWidth(2);
1081     fhC55[i]->SetLineColor(color);   fhC55[i]->SetLineWidth(2);
1082
1083     fhRel1PtUncertainty[i]->SetLineColor(color); fhRel1PtUncertainty[i]->SetLineWidth(2);
1084
1085     fhDXY[i]->SetLineColor(color);   fhDXY[i]->SetLineWidth(2);
1086     fhDZ[i]->SetLineColor(color);    fhDZ[i]->SetLineWidth(2);
1087     fhDXYDZ[i]->SetLineColor(color); fhDXYDZ[i]->SetLineWidth(2);
1088
1089     fhDXYNormalized[i]->SetLineColor(color);   fhDXYNormalized[i]->SetLineWidth(2);
1090     fhDZNormalized[i]->SetLineColor(color);    fhDZNormalized[i]->SetLineWidth(2);
1091     fhNSigmaToVertex[i]->SetLineColor(color);  fhNSigmaToVertex[i]->SetLineWidth(2);
1092   }
1093
1094   // The number of sigmas to the vertex is per definition gaussian
1095   ffDTheoretical = new TF1("nSigmaToVertexTheoretical","([0]/2.506628274)*exp(-(x**2)/2)",0,50);
1096   ffDTheoretical->SetParameter(0,1);
1097
1098   TH1::AddDirectory(oldStatus);
1099 }
1100
1101 //____________________________________________________________________
1102 Bool_t AliESDtrackCuts::LoadHistograms(const Char_t* dir)
1103 {
1104   //
1105   // loads the histograms from a file
1106   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
1107   //
1108
1109   if (!dir)
1110     dir = GetName();
1111
1112   if (!gDirectory->cd(dir))
1113     return kFALSE;
1114
1115   ffDTheoretical = dynamic_cast<TF1*> (gDirectory->Get("nSigmaToVertexTheory"));
1116
1117   fhCutStatistics = dynamic_cast<TH1F*> (gDirectory->Get("cut_statistics"));
1118   fhCutCorrelation = dynamic_cast<TH2F*> (gDirectory->Get("cut_correlation"));
1119
1120   for (Int_t i=0; i<2; i++) {
1121     if (i==0)
1122     {
1123       gDirectory->cd("before_cuts");
1124     }
1125     else
1126       gDirectory->cd("after_cuts");
1127
1128     fhNClustersITS[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersITS"     ));
1129     fhNClustersTPC[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersTPC"     ));
1130     fhChi2PerClusterITS[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterITS"));
1131     fhChi2PerClusterTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterTPC"));
1132
1133     fhC11[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal11"));
1134     fhC22[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal22"));
1135     fhC33[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal33"));
1136     fhC44[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal44"));
1137     fhC55[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal55"));
1138
1139     fhRel1PtUncertainty[i] = dynamic_cast<TH1F*> (gDirectory->Get("rel1PtUncertainty"));
1140
1141     fhDXY[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXY"    ));
1142     fhDZ[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZ"     ));
1143     fhDXYDZ[i] =   dynamic_cast<TH1F*> (gDirectory->Get("dXYDZ"));
1144     fhDXYvsDZ[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZ"));
1145
1146     fhDXYNormalized[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXYNormalized"    ));
1147     fhDZNormalized[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZNormalized"     ));
1148     fhDXYvsDZNormalized[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZNormalized"));
1149     fhNSigmaToVertex[i] = dynamic_cast<TH1F*> (gDirectory->Get("nSigmaToVertex"));
1150
1151     fhPt[i] = dynamic_cast<TH1F*> (gDirectory->Get("pt"));
1152     fhEta[i] = dynamic_cast<TH1F*> (gDirectory->Get("eta"));
1153
1154     gDirectory->cd("../");
1155   }
1156
1157   gDirectory->cd("..");
1158
1159   return kTRUE;
1160 }
1161
1162 //____________________________________________________________________
1163 void AliESDtrackCuts::SaveHistograms(const Char_t* dir) {
1164   //
1165   // saves the histograms in a directory (dir)
1166   //
1167
1168   if (!fHistogramsOn) {
1169     AliDebug(0, "Histograms not on - cannot save histograms!!!");
1170     return;
1171   }
1172
1173   if (!dir)
1174     dir = GetName();
1175
1176   gDirectory->mkdir(dir);
1177   gDirectory->cd(dir);
1178
1179   gDirectory->mkdir("before_cuts");
1180   gDirectory->mkdir("after_cuts");
1181
1182   // a factor of 2 is needed since n sigma is positive
1183   ffDTheoretical->SetParameter(0,2*fhNSigmaToVertex[0]->Integral("width"));
1184   ffDTheoretical->Write("nSigmaToVertexTheory");
1185
1186   fhCutStatistics->Write();
1187   fhCutCorrelation->Write();
1188
1189   for (Int_t i=0; i<2; i++) {
1190     if (i==0)
1191       gDirectory->cd("before_cuts");
1192     else
1193       gDirectory->cd("after_cuts");
1194
1195     fhNClustersITS[i]        ->Write();
1196     fhNClustersTPC[i]        ->Write();
1197     fhChi2PerClusterITS[i]   ->Write();
1198     fhChi2PerClusterTPC[i]   ->Write();
1199
1200     fhC11[i]                 ->Write();
1201     fhC22[i]                 ->Write();
1202     fhC33[i]                 ->Write();
1203     fhC44[i]                 ->Write();
1204     fhC55[i]                 ->Write();
1205
1206     fhRel1PtUncertainty[i]   ->Write();
1207
1208     fhDXY[i]                 ->Write();
1209     fhDZ[i]                  ->Write();
1210     fhDXYDZ[i]               ->Write();
1211     fhDXYvsDZ[i]             ->Write();
1212
1213     fhDXYNormalized[i]       ->Write();
1214     fhDZNormalized[i]        ->Write();
1215     fhDXYvsDZNormalized[i]   ->Write();
1216     fhNSigmaToVertex[i]      ->Write();
1217
1218     fhPt[i]                  ->Write();
1219     fhEta[i]                 ->Write();
1220
1221     gDirectory->cd("../");
1222   }
1223
1224   gDirectory->cd("../");
1225 }
1226
1227 //____________________________________________________________________
1228 void AliESDtrackCuts::DrawHistograms()
1229 {
1230   // draws some histograms
1231
1232   TCanvas* canvas1 = new TCanvas(Form("%s_1", GetName()), "Track Quality Results1", 800, 800);
1233   canvas1->Divide(2, 2);
1234
1235   canvas1->cd(1);
1236   fhNClustersTPC[0]->SetStats(kFALSE);
1237   fhNClustersTPC[0]->Draw();
1238
1239   canvas1->cd(2);
1240   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1241   fhChi2PerClusterTPC[0]->Draw();
1242
1243   canvas1->cd(3);
1244   fhNSigmaToVertex[0]->SetStats(kFALSE);
1245   fhNSigmaToVertex[0]->GetXaxis()->SetRangeUser(0, 10);
1246   fhNSigmaToVertex[0]->Draw();
1247
1248   canvas1->SaveAs(Form("%s_%s.gif", GetName(), canvas1->GetName()));
1249
1250   TCanvas* canvas2 = new TCanvas(Form("%s_2", GetName()), "Track Quality Results2", 1200, 800);
1251   canvas2->Divide(3, 2);
1252
1253   canvas2->cd(1);
1254   fhC11[0]->SetStats(kFALSE);
1255   gPad->SetLogy();
1256   fhC11[0]->Draw();
1257
1258   canvas2->cd(2);
1259   fhC22[0]->SetStats(kFALSE);
1260   gPad->SetLogy();
1261   fhC22[0]->Draw();
1262
1263   canvas2->cd(3);
1264   fhC33[0]->SetStats(kFALSE);
1265   gPad->SetLogy();
1266   fhC33[0]->Draw();
1267
1268   canvas2->cd(4);
1269   fhC44[0]->SetStats(kFALSE);
1270   gPad->SetLogy();
1271   fhC44[0]->Draw();
1272
1273   canvas2->cd(5);
1274   fhC55[0]->SetStats(kFALSE);
1275   gPad->SetLogy();
1276   fhC55[0]->Draw();
1277
1278   canvas2->cd(6);
1279   fhRel1PtUncertainty[0]->SetStats(kFALSE);
1280   gPad->SetLogy();
1281   fhRel1PtUncertainty[0]->Draw();
1282
1283   canvas2->SaveAs(Form("%s_%s.gif", GetName(), canvas2->GetName()));
1284
1285   TCanvas* canvas3 = new TCanvas(Form("%s_3", GetName()), "Track Quality Results3", 1200, 800);
1286   canvas3->Divide(3, 2);
1287
1288   canvas3->cd(1);
1289   fhDXY[0]->SetStats(kFALSE);
1290   gPad->SetLogy();
1291   fhDXY[0]->Draw();
1292
1293   canvas3->cd(2);
1294   fhDZ[0]->SetStats(kFALSE);
1295   gPad->SetLogy();
1296   fhDZ[0]->Draw();
1297
1298   canvas3->cd(3);
1299   fhDXYvsDZ[0]->SetStats(kFALSE);
1300   gPad->SetLogz();
1301   gPad->SetRightMargin(0.15);
1302   fhDXYvsDZ[0]->Draw("COLZ");
1303
1304   canvas3->cd(4);
1305   fhDXYNormalized[0]->SetStats(kFALSE);
1306   gPad->SetLogy();
1307   fhDXYNormalized[0]->Draw();
1308
1309   canvas3->cd(5);
1310   fhDZNormalized[0]->SetStats(kFALSE);
1311   gPad->SetLogy();
1312   fhDZNormalized[0]->Draw();
1313
1314   canvas3->cd(6);
1315   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1316   gPad->SetLogz();
1317   gPad->SetRightMargin(0.15);
1318   fhDXYvsDZNormalized[0]->Draw("COLZ");
1319
1320   canvas3->SaveAs(Form("%s_%s.gif", GetName(), canvas3->GetName()));
1321
1322   TCanvas* canvas4 = new TCanvas(Form("%s_4", GetName()), "Track Quality Results4", 800, 500);
1323   canvas4->Divide(2, 1);
1324
1325   canvas4->cd(1);
1326   fhCutStatistics->SetStats(kFALSE);
1327   fhCutStatistics->LabelsOption("v");
1328   gPad->SetBottomMargin(0.3);
1329   fhCutStatistics->Draw();
1330
1331   canvas4->cd(2);
1332   fhCutCorrelation->SetStats(kFALSE);
1333   fhCutCorrelation->LabelsOption("v");
1334   gPad->SetBottomMargin(0.3);
1335   gPad->SetLeftMargin(0.3);
1336   fhCutCorrelation->Draw("COLZ");
1337
1338   canvas4->SaveAs(Form("%s_%s.gif", GetName(), canvas4->GetName()));
1339
1340   /*canvas->cd(1);
1341   fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1342   fhDXYvsDZNormalized[0]->DrawCopy("COLZ");
1343
1344   canvas->cd(2);
1345   fhNClustersTPC[0]->SetStats(kFALSE);
1346   fhNClustersTPC[0]->DrawCopy();
1347
1348   canvas->cd(3);
1349   fhChi2PerClusterITS[0]->SetStats(kFALSE);
1350   fhChi2PerClusterITS[0]->DrawCopy();
1351   fhChi2PerClusterITS[1]->SetLineColor(2);
1352   fhChi2PerClusterITS[1]->DrawCopy("SAME");
1353
1354   canvas->cd(4);
1355   fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1356   fhChi2PerClusterTPC[0]->DrawCopy();
1357   fhChi2PerClusterTPC[1]->SetLineColor(2);
1358   fhChi2PerClusterTPC[1]->DrawCopy("SAME");*/
1359 }
1360