ANALYSIS/AliESDtrackCuts.cxx

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