PWGHF/hfe/AliHFEtrdPIDqa.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 // QA class for TRD PID
  17 // Plot Pion Efficiency at x electron efficiency
  18 // Calculate the threshold parametrisation and save
  19 // them in a root file
  20 //
  21 // Author:
  22 //   Markus Fasel <M.Fasel@gsi.de>
  23 //
  24 #include <TAxis.h>
  25 #include <TBrowser.h>
  26 #include <TClass.h>
  27 #include <TCanvas.h>
  28 #include <TF1.h>
  29 #include <TFile.h>
  30 #include <TGraph.h>
  31 #include <TGraphErrors.h>
  32 #include <THnSparse.h>
  33 #include <TH1.h>
  34 #include <TH2.h>
  35 #include <TIterator.h>
  36 #include <TLegend.h>
  37 #include <TList.h>
  38 #include <TMath.h>
  39 #include <TObjArray.h>
  40 #include <TString.h>
  41
  42 #include "AliAODTrack.h"
  43 #include "AliAODPid.h"
  44 #include "AliESDtrack.h"
  45 #include "AliHFEtrdPIDqa.h"
  46 #include "AliHFEtools.h"
  47 #include "AliHFEpidTRD.h"
  48 #include "AliLog.h"
  49
  50 ClassImp(AliHFEtrdPIDqa)
  51
  52 const Double_t AliHFEtrdPIDqa::fgkElectronEff[kNElectronEffs] = {0.7,0.75, 0.8, 0.85, 0.9, 0.95};
  53 //_______________________________________________________________
  54 // Definition of the common binning
  55 const Int_t AliHFEtrdPIDqa::fgkNBinsCommon[kQuantitiesCommon] = {
  56   AliPID::kSPECIES + 1,         // species
  57   40,                           // p-bins
  58   AliESDtrack::kTRDnPlanes + 1 // tracklets including 0
  59 };
  60 const Double_t AliHFEtrdPIDqa::fgkMinBinCommon[kQuantitiesCommon] = {
  61   -1,      // species
  62   0.1,     // p-bins:
  63   0        // tracklets including 0
  64 };
  65
  66 const Double_t AliHFEtrdPIDqa::fgkMaxBinCommon[kQuantitiesCommon] = {
  67   AliPID::kSPECIES,               // species
  68   10.,                            // p-bins
  69   AliESDtrack::kTRDnPlanes + 1    // tracklets including 0
  70 };
  71 //_______________________________________________________________
  72
  73 //__________________________________________________________________
  74 AliHFEtrdPIDqa::AliHFEtrdPIDqa():
  75   TNamed("trdPIDqa", ""),
  76   fTRDpid(NULL),
  77   fHistos(NULL),
  78   fPionEfficiencies(NULL),
  79   fProtonEfficiencies(NULL),
  80   fKaonEfficiencies(NULL),
  81   fThresholds(NULL),
  82   fShowMessage(kFALSE),
  83   fTotalChargeInSlice0(kFALSE),
  84   fCentralityBin(-1)
  85 {
  86   //
  87   // Default Constructor
  88   //
  89 }
  90
  91 //__________________________________________________________________
  92 AliHFEtrdPIDqa::AliHFEtrdPIDqa(const Char_t *name):
  93   TNamed(name, ""),
  94   fTRDpid(NULL),
  95   fHistos(NULL),
  96   fPionEfficiencies(NULL),
  97   fProtonEfficiencies(NULL),
  98   fKaonEfficiencies(NULL),
  99   fThresholds(NULL),
 100   fShowMessage(kFALSE),
 101   fTotalChargeInSlice0(kFALSE),
 102   fCentralityBin(-1)
 103 {
 104   //
 105   // Main Constructor
 106   //
 107 }
 108
 109 //__________________________________________________________________
 110 AliHFEtrdPIDqa::AliHFEtrdPIDqa(const AliHFEtrdPIDqa &ref):
 111   TNamed(ref),
 112   fTRDpid(NULL),
 113   fHistos(NULL),
 114   fPionEfficiencies(NULL),
 115   fProtonEfficiencies(NULL),
 116   fKaonEfficiencies(NULL),
 117   fThresholds(NULL),
 118   fShowMessage(kFALSE),
 119   fTotalChargeInSlice0(ref.fTotalChargeInSlice0),
 120   fCentralityBin(ref.fCentralityBin)
 121 {
 122   //
 123   // Copy constructor
 124   //
 125   ref.Copy(*this);
 126 }
 127
 128 //__________________________________________________________________
 129 AliHFEtrdPIDqa &AliHFEtrdPIDqa::operator=(const AliHFEtrdPIDqa &ref){
 130   //
 131   // Assignment operator
 132   //
 133   if(this != &ref)
 134     ref.Copy(*this);
 135   return *this;
 136 }
 137
 138 //__________________________________________________________________
 139 AliHFEtrdPIDqa::~AliHFEtrdPIDqa(){
 140   //
 141   // Destructor
 142   //
 143   if(fTRDpid) delete fTRDpid;
 144   if(fHistos) delete fHistos;
 145   if(fPionEfficiencies) delete fPionEfficiencies;
 146   if(fProtonEfficiencies) delete fProtonEfficiencies;
 147   if(fKaonEfficiencies) delete fKaonEfficiencies;
 148 }
 149
 150 //__________________________________________________________________
 151 void AliHFEtrdPIDqa::Copy(TObject &ref) const{
 152   //
 153   // Copies content of this object into object ref
 154   //
 155   TNamed::Copy(ref);
 156
 157   AliHFEtrdPIDqa &target = dynamic_cast<AliHFEtrdPIDqa &>(ref);
 158   target.fTRDpid = fTRDpid;
 159   target.fHistos = dynamic_cast<AliHFEcollection *>(fHistos->Clone());
 160   target.fTotalChargeInSlice0 = fTotalChargeInSlice0;
 161   target.fCentralityBin = fCentralityBin;
 162 }
 163
 164 //__________________________________________________________________
 165 Long64_t AliHFEtrdPIDqa::Merge(TCollection *coll){
 166   //
 167   // Merge objects
 168   //
 169   if(!coll) return 0;
 170   if(coll->IsEmpty()) return 1;
 171
 172   AliHFEtrdPIDqa *refQA = NULL;
 173   TIter it(coll);
 174   TObject *o = NULL;
 175   Long64_t count = 0;
 176   TList listHistos;
 177   while((o = it())){
 178     refQA = dynamic_cast<AliHFEtrdPIDqa *>(o);
 179     if(!refQA) continue;
 180
 181     listHistos.Add(refQA->fHistos);
 182     count++;
 183   }
 184   fHistos->Merge(&listHistos);
 185   return count+1;
 186 }
 187
 188 //__________________________________________________________________
 189 void AliHFEtrdPIDqa::Browse(TBrowser *b){
 190   //
 191   // Enable Browser functionality
 192   //
 193   if(b){
 194     // Add objects to the browser
 195     if(fHistos) b->Add(fHistos, fHistos->GetName());
 196     if(fPionEfficiencies) b->Add(fPionEfficiencies, "Pion Efficiencies");
 197     if(fProtonEfficiencies) b->Add(fProtonEfficiencies, "Proton Efficiencies");
 198     if(fKaonEfficiencies) b->Add(fKaonEfficiencies, "Kaon Efficiencies");
 199     if(fThresholds) b->Add(fThresholds, "Thresholds");
 200   }
 201 }
 202
 203 //__________________________________________________________________
 204 void AliHFEtrdPIDqa::Init(){
 205   //
 206   // Initialize Object
 207   //
 208
 209   fHistos = new AliHFEcollection("TRDqa", "Histos for TRD QA");
 210
 211   CreateLikelihoodHistogram();
 212   CreateQAHistogram();
 213   CreatedEdxHistogram();
 214   CreateHistoTruncatedMean();
 215
 216   fTRDpid = new AliHFEpidTRD("QAtrdPID");
 217   if(fTotalChargeInSlice0) fTRDpid->SetTotalChargeInSlice0();
 218 }
 219
 220 //__________________________________________________________________
 221 void AliHFEtrdPIDqa::CreateLikelihoodHistogram(){
 222   //
 223   // Create Histogram for TRD Likelihood Studies
 224   //
 225   Int_t nbins[kQuantitiesLike]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
 226   nbins[kElectronLike] = 100;
 227   nbins[kNClustersLike] = 200;
 228   nbins[kCentralityBinLike] = 12;
 229   Double_t binMin[kQuantitiesLike]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
 230   Double_t binMax[kQuantitiesLike]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
 231   binMin[kElectronLike] = 0.;
 232   binMin[kNClustersLike] = 0.;
 233   binMin[kCentralityBinLike] = -1.;
 234   binMax[kElectronLike] = 1.;
 235   binMax[kNClustersLike] = 200.;
 236   binMax[kCentralityBinLike] = 11.;
 237
 238   fHistos->CreateTHnSparse("fLikeTRD","TRD Likelihood Studies", kQuantitiesLike, nbins, binMin, binMax);
 239   fHistos->BinLogAxis("fLikeTRD", kP);
 240 }
 241
 242 //__________________________________________________________________
 243 void AliHFEtrdPIDqa::CreateQAHistogram(){
 244   //
 245   // Create Histogram for Basic TRD PID QA
 246   //
 247   AliDebug(1, "Called");
 248   Int_t nbins[kQuantitiesQA]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
 249   nbins[kNonZeroTrackletCharge] = AliESDtrack::kTRDnPlanes + 1;
 250   nbins[kNClusters] = 200;
 251   nbins[kCentralityBinQA] = 12;
 252   Double_t binMin[kQuantitiesQA]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
 253   binMin[kNonZeroTrackletCharge] = 0.;
 254   binMin[kNClusters] = 0.;
 255   binMin[kCentralityBinQA] = -1.;
 256   Double_t binMax[kQuantitiesQA]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
 257   binMax[kNonZeroTrackletCharge] = AliESDtrack::kTRDnPlanes + 1.;
 258   binMax[kNClusters] = 200.;
 259   binMax[kCentralityBinQA] = 11.;
 260   fHistos->CreateTHnSparse("fQAtrack","TRD QA Histogram", kQuantitiesQA, nbins, binMin, binMax);
 261   fHistos->BinLogAxis("fQAtrack", kP);
 262 }
 263
 264 //__________________________________________________________________
 265 void AliHFEtrdPIDqa::CreatedEdxHistogram(){
 266   //
 267   // Create QA histogram for dEdx investigations
 268   //
 269   AliDebug(1, "Called");
 270   Int_t nbins[kQuantitiesdEdx]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
 271   nbins[kdEdx] = 100;
 272   nbins[kNclusters] = 261;
 273   nbins[kNonZeroSlices] = 9;
 274   nbins[kCentralityBindEdx] = 12;
 275   Double_t binMin[kQuantitiesdEdx]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
 276   binMin[kdEdx] = 0.;
 277   binMin[kNclusters] = 0;
 278   binMin[kNonZeroSlices] = 0.;
 279   binMin[kCentralityBindEdx] = -1.;
 280   Double_t binMax[kQuantitiesdEdx]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
 281   binMax[kdEdx] = 10000.;
 282   binMax[kNclusters] = 260.;
 283   binMax[kNonZeroSlices] = 8.;
 284   binMax[kCentralityBindEdx] = 11.;
 285
 286   fHistos->CreateTHnSparse("fQAdEdx","TRD summed dEdx", kQuantitiesdEdx, nbins, binMin, binMax);
 287   fHistos->BinLogAxis("fQAdEdx", kP);
 288   fHistos->Sumw2("fQAdEdx");
 289 }
 290
 291 //__________________________________________________________________
 292 void AliHFEtrdPIDqa::CreateHistoTruncatedMean(){
 293   //
 294   // Create Histogram for Basic TRD PID QA:
 295   //
 296   AliDebug(1, "Called");
 297   Int_t nbins[kQuantitiesTruncMean]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
 298   nbins[kTPCdEdx] = 600;
 299   nbins[kTRDdEdxMethod1] = 1000;
 300   nbins[kTRDdEdxMethod2] = 1000;
 301   nbins[kCentralityBinTruncMean] = 12;
 302   Double_t binMin[kQuantitiesTruncMean]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
 303   binMin[kTPCdEdx] = 0.;
 304   binMin[kTRDdEdxMethod1] = 0.;
 305   binMin[kTRDdEdxMethod2] = 0.;
 306   binMin[kCentralityBinTruncMean] = -1.;
 307   Double_t binMax[kQuantitiesTruncMean]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
 308   binMax[kTPCdEdx] = 600;
 309   binMax[kTRDdEdxMethod1] = 20000.;
 310   binMax[kTRDdEdxMethod2] = 20000.;
 311   binMax[kCentralityBinTruncMean] = 11.;
 312
 313   fHistos->CreateTHnSparse("fTRDtruncMean","TRD TruncatedMean studies", kQuantitiesTruncMean, nbins, binMin, binMax);
 314   fHistos->BinLogAxis("fTRDtruncMean", kP);
 315   fHistos->CreateTH2F("fTRDslicesPions","TRD dEdx per slice for Pions", 8, 0, 8, 2000, 0, 8000);
 316   fHistos->CreateTH2F("fTRDslicesElectrons","TRD dEdx per slice for Electrons", 8, 0, 8, 2000, 0, 8000);
 317 }
 318
 319
 320 //__________________________________________________________________
 321 void AliHFEtrdPIDqa::ProcessTracks(const TObjArray * const tracks, Int_t species){
 322   //
 323   // Process Electron Tracks
 324   //
 325   if(species < -1 || species >= AliPID::kSPECIES) return;
 326   TIter it(tracks);
 327   const AliVTrack *track = NULL;
 328   while((track = dynamic_cast<const AliVTrack *>(it()))){
 329     if(track) ProcessTrack(track, species);
 330   }
 331 }
 332
 333 //__________________________________________________________________
 334 void AliHFEtrdPIDqa::ProcessTrack(const AliVTrack * const track, Int_t species){
 335   //
 336   // Process Single Track
 337   //
 338   if(TString(track->IsA()->GetName()).CompareTo("AliESDtrack") == 0)
 339     ProcessTrackESD(dynamic_cast<const AliESDtrack *>(track), species);
 340   else if(TString(track->IsA()->GetName()).CompareTo("AliAODTrack") == 0)
 341     ProcessTrackAOD(dynamic_cast<const AliAODTrack *>(track), species);
 342 }
 343
 344
 345 //__________________________________________________________________
 346 void AliHFEtrdPIDqa::ProcessTrackESD(const AliESDtrack *track, Int_t species){
 347   //
 348   // Process single ESD track
 349   //
 350   if(!track) return;
 351   if((track->GetStatus() & AliESDtrack::kTRDout) == 0) return;  // require TRD track
 352   FillTRDLikelihoods(track, species);
 353   FillTRDQAplots(track, species);
 354 }
 355
 356 //__________________________________________________________________
 357 void AliHFEtrdPIDqa::ProcessTrackAOD(const AliAODTrack * const track, Int_t /*species*/){
 358   //
 359   // Process single AOD track
 360   // AOD PID object required
 361   //
 362   if(!track) return;
 363   AliAODPid *trackPID = track->GetDetPid();
 364   if(!trackPID) return;
 365
 366 }
 367
 368 //__________________________________________________________________
 369 void AliHFEtrdPIDqa::FillTRDLikelihoods(const AliESDtrack * const track, Int_t species){
 370   //
 371   // Fill TRD Likelihood Histogram
 372   //
 373   Double_t trdLike[AliPID::kSPECIES];
 374   track->GetTRDpid(trdLike);
 375   // Renormalize
 376   Double_t norm =trdLike[AliPID::kElectron]+trdLike[AliPID::kPion];
 377   Double_t likeEle = norm == 0. ? 0. : trdLike[AliPID::kElectron]/norm;
 378   const AliExternalTrackParam *outerPars = track->GetOuterParam();
 379
 380   //Int_t kQuantitiesLike;
 381   Double_t quantities[kQuantitiesLike]; memset(quantities, 0, sizeof(Double_t) * kQuantitiesLike);
 382   // we store:
 383   // species
 384   // p
 385   // ntracklets
 386   // Electron Likelihood
 387   // Centrality
 388   quantities[kSpecies] = species;
 389   quantities[kP] = outerPars ? outerPars->P() : track->P();
 390   quantities[kNTracklets] = track->GetTRDntrackletsPID();
 391   quantities[kElectronLike] = likeEle;
 392   quantities[kNClustersLike] =  track->GetTRDncls();
 393   quantities[kCentralityBinLike] = fCentralityBin;
 394   fHistos->Fill("fLikeTRD", quantities);
 395
 396 }
 397
 398 //__________________________________________________________________
 399 void AliHFEtrdPIDqa::FillTRDQAplots(const AliESDtrack * const track, Int_t species){
 400   //
 401   // Fill QA Plots containing further information
 402   //
 403   const AliExternalTrackParam *outerPars = track->GetOuterParam();
 404
 405   Double_t quantitiesQA[kQuantitiesQA], quantitiesdEdx[kQuantitiesdEdx], quantitiesTruncMean[kQuantitiesTruncMean];
 406   // we store:
 407   // 1. QA
 408   // species
 409   // p
 410   // ntracklets
 411   // Non-zero tracklet charges
 412   // Number of clusters / full track
 413   // 2. dEdx
 414   // species
 415   // p
 416   // ntracklets
 417   // dEdx
 418   // 3. Truncated Mean
 419   // ...
 420   // TPC dEdx
 421   // TRD dEdx Method 1
 422   // TRD dEdx Method 2
 423   quantitiesQA[kSpecies]  = quantitiesdEdx[kSpecies]
 424                           = quantitiesTruncMean[kSpecies]
 425                           = species;
 426   quantitiesQA[kP]  = quantitiesTruncMean[kP]
 427                     = outerPars ? outerPars->P() : track->P();
 428   quantitiesQA[kNTracklets] = quantitiesdEdx[kNTracklets]
 429                             = quantitiesTruncMean[kNTracklets]
 430                             = track->GetTRDntrackletsPID();
 431   quantitiesQA[kNClusters] = quantitiesdEdx[kNclusters] = track->GetTRDncls();
 432
 433
 434   Double_t dEdxSum = 0., qSlice = 0.;
 435   // remove the last slice from the histogram
 436   Int_t ntrackletsNonZero = 0, nSlices = track->GetNumberOfTRDslices(), nSlicesNonZero = 0;
 437   TString speciesname = "pions";
 438   Bool_t selectedForSlicemon = kFALSE;
 439
 440   switch(species){
 441     case AliPID::kElectron: speciesname = "Electrons"; selectedForSlicemon = kTRUE; break;
 442     case AliPID::kPion: speciesname = "Pions"; selectedForSlicemon = kTRUE; break;
 443     default: speciesname = "undefined"; selectedForSlicemon = kFALSE; break;
 444   };
 445   AliDebug(1, Form("species %d, speciesname %s, momentum %f, selected %s", species, speciesname.Data(), track->P(), selectedForSlicemon ? "yes" : "no"));
 446   for(Int_t iplane = 0; iplane < AliESDtrack::kTRDnPlanes; iplane++){
 447     quantitiesdEdx[kP] = track->GetTRDmomentum(iplane);
 448     dEdxSum = 0.;
 449     for(Int_t islice = 0; islice < nSlices; islice++){
 450       if(fTotalChargeInSlice0 && islice >= 7) break;
 451       qSlice = track->GetTRDslice(iplane, fTotalChargeInSlice0 ? islice + 1 : islice);  // hack by mfasel: For data with the new reconstruction, slice 0 is used to store the total charge, the total number of slices is 7 instead of 8
 452       if(qSlice > 1e-1){
 453         // cut out 0 slices
 454         nSlicesNonZero++;
 455         dEdxSum += qSlice;
 456         // Reweighting of the slices for the truncated mean: select all pion tracks above
 457         // 1.5 GeV and monitor the dEdx as function of slice
 458         if(selectedForSlicemon && track->P() > 1.5){
 459           AliDebug(2, Form("Filling Histogram fTRDslices%s", speciesname.Data()));
 460           fHistos->Fill(Form("fTRDslices%s", speciesname.Data()), static_cast<Double_t>(islice), qSlice);
 461         }
 462       }
 463     }
 464     quantitiesdEdx[kNonZeroSlices] = nSlicesNonZero;
 465     quantitiesdEdx[kdEdx] = fTotalChargeInSlice0 ? track->GetTRDslice(iplane, 0) : dEdxSum; // hack by mfasel: In the new reconstruction, the total charge is stored in the first slice, in the old reconstruction it has to be calculated from the slice charges.
 466     if(dEdxSum) ntrackletsNonZero++;
 467     // Fill dEdx histogram
 468     quantitiesdEdx[kCentralityBindEdx] = fCentralityBin;
 469     if(dEdxSum > 1e-1) fHistos->Fill("fQAdEdx", quantitiesdEdx); // Cut out 0 entries
 470   }
 471   quantitiesQA[kNonZeroTrackletCharge] = ntrackletsNonZero;
 472   quantitiesQA[kCentralityBinQA] = fCentralityBin;
 473   fHistos->Fill("fQAtrack", quantitiesQA);
 474
 475   quantitiesTruncMean[kTPCdEdx] = track->GetTPCsignal();
 476   quantitiesTruncMean[kTRDdEdxMethod1] = fTRDpid->GetTRDSignalV1(track, 0.6);
 477   quantitiesTruncMean[kTRDdEdxMethod2] = fTRDpid->GetTRDSignalV2(track, 0.6);
 478   quantitiesTruncMean[kCentralityBinTruncMean] = fCentralityBin;
 479   fHistos->Fill("fTRDtruncMean", quantitiesTruncMean);
 480 }
 481
 482 /////////////////////////////////////////////////////////
 483 //
 484 // Code for Post Processing
 485 //
 486 // //////////////////////////////////////////////////////
 487
 488 //__________________________________________________________________
 489 void AliHFEtrdPIDqa::FinishAnalysis(Int_t fCentrality, Bool_t isGreaterEqual){
 490   //
 491   // Finish Analysis:
 492   // Calculate Electron Efficiency for ntracklets = 4...6
 493   // Calculate thresholds for ntracklets = 4...6
 494   //
 495
 496   if(!fPionEfficiencies){
 497     fPionEfficiencies = new TList;
 498     fPionEfficiencies->SetName("pionEfficiencies");
 499     fPionEfficiencies->SetOwner();
 500   }
 501   if(!fProtonEfficiencies){
 502     fProtonEfficiencies = new TList;
 503     fProtonEfficiencies->SetName("protonEfficiencies");
 504     fProtonEfficiencies->SetOwner();
 505   }
 506   if(!fThresholds){
 507     fThresholds = new TList;
 508     fThresholds->SetName("thresholds");
 509     fThresholds->SetOwner();
 510   }
 511
 512   for(Int_t itr = 4; itr <= 6; itr++){
 513     if(fShowMessage){
 514       printf("========================================\n");
 515       printf("Analysing %d trackltes centrality %i \n", itr, fCentrality);
 516       printf("========================================\n");
 517     }
 518     AnalyseNTracklets(itr, fCentrality, isGreaterEqual);
 519   }
 520 }
 521
 522 //__________________________________________________________________
 523 void AliHFEtrdPIDqa::StoreResults(const Char_t *filename){
 524   //
 525   // Store histos into a root file
 526   //
 527   TFile *outfile = new TFile(filename, "RECREATE");
 528   outfile->cd();
 529   fPionEfficiencies->Clone()->Write(fPionEfficiencies->GetName(), kSingleKey);
 530   fProtonEfficiencies->Clone()->Write(fProtonEfficiencies->GetName(), kSingleKey);
 531   fThresholds->Clone()->Write(fThresholds->GetName(), kSingleKey);
 532   outfile->Close();
 533   delete outfile;
 534 }
 535
 536 //__________________________________________________________________
 537 void AliHFEtrdPIDqa::SaveThresholdParameters(const Char_t *name, Double_t lowerLimit, Double_t upperLimit, Int_t icentrality){
 538   //
 539   // Fit the threshold histograms with the given parametrisation
 540   // and store the TF1 in the file
 541   //
 542
 543   if(!fThresholds){
 544     AliError("Threshold histograms have to be created first");
 545     return;
 546   }
 547
 548     if(fShowMessage){
 549     printf("========================================\n");
 550     printf("Calculating threshold parameters\n");
 551     printf("========================================\n");
 552   }
 553
 554   TList *outlist = new TList;
 555   outlist->SetName("thresholdTRD");
 556
 557   TGraph *threshhist = NULL;
 558   TF1 *threshparam = NULL;
 559   TList *lHistos = NULL, *lFormulas = NULL;
 560   for(Int_t itracklet = 4; itracklet <= 6; itracklet++){
 561
 562     if(fShowMessage){
 563       printf("-------------------------------\n");
 564       printf("Processing %d tracklets\n", itracklet);
 565       printf("-------------------------------\n");
 566     }
 567
 568     Char_t *listname=Form("%dTracklets", itracklet);
 569     if(icentrality!=-1) listname=Form("%dTracklets%dCentrality", itracklet, icentrality);
 570     lHistos = dynamic_cast<TList *>(fThresholds->FindObject(listname));
 571     if(!lHistos){
 572       AliError(Form("Threshold histograms for the case %s not found", listname));
 573       continue;
 574     }
 575     lFormulas = new TList;
 576     lFormulas->SetName(listname);
 577     outlist->Add(lFormulas);
 578
 579     for(Int_t ieff = 0; ieff <  kNElectronEffs; ieff++){
 580
 581       if(fShowMessage){
 582         printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n");
 583         printf("Processing Electron Efficiency %f\n", fgkElectronEff[ieff]);
 584         printf(">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>\n");
 585       }
 586
 587       threshhist = dynamic_cast<TGraph *>(lHistos->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
 588       if(!threshhist) continue;
 589       threshparam = MakeThresholds(threshhist, lowerLimit, upperLimit);
 590       threshparam->SetName(Form("thresh_%d_%d", itracklet, static_cast<Int_t>(fgkElectronEff[ieff] * 100)));
 591       if(icentrality!=-1) threshparam->SetName(Form("thresh_%d_%d_%d", itracklet, icentrality, static_cast<Int_t>(fgkElectronEff[ieff] * 100)));
 592       lFormulas->Add(threshparam);
 593     }
 594   }
 595
 596   // store the output
 597   TFile *outfile = new TFile(name, "RECREATE");
 598   outfile->cd();
 599   outlist->Write(outlist->GetName(), kSingleKey);
 600   outfile->Close();
 601   delete outfile;
 602 }
 603
 604 //__________________________________________________________________
 605 void AliHFEtrdPIDqa::AnalyseNTracklets(Int_t nTracklets, Int_t nCentrality, Bool_t isGreaterEqual){
 606   //
 607   // Calculate Pion Efficiency, Proton Efficiency and Kaon Efficiency at discrete
 608   // elPion Efficiency, Proton Efficiency and Kaon Efficiency at discrete
 609   // electron efficiencies
 610   //
 611   THnSparse *hLikeTRD = dynamic_cast<THnSparseF *>(fHistos->Get("fLikeTRD"));
 612   if(!hLikeTRD){
 613     AliError("Likelihood Histogram not available");
 614     return;
 615   }
 616
 617   Bool_t isPbPb = kFALSE;
 618   if(nCentrality==-1)  isPbPb = kFALSE;
 619   if(nCentrality!=-1)  isPbPb = kTRUE;
 620
 621   Int_t binTracklets = hLikeTRD->GetAxis(kNTracklets)->FindBin(nTracklets);
 622   hLikeTRD->GetAxis(kNTracklets)->SetRange(binTracklets, isGreaterEqual ? 7 : binTracklets);
 623
 624   if(isPbPb){
 625       Int_t binCentrality = hLikeTRD->GetAxis(kCentralityBinLike)->FindBin(nCentrality);
 626       hLikeTRD->GetAxis(kCentralityBinLike)->SetRange(binCentrality, isGreaterEqual ? 11 : binCentrality);
 627       /*
 628        new TCanvas;
 629        TH2 *test = hLikeTRD->Projection(kCentralityBin, kP);
 630        test->Draw("colz");
 631        */
 632   }
 633
 634   Int_t binElectrons = hLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kElectron);
 635   AliDebug(1, Form("BinElectrons %d", binElectrons));
 636   Int_t binPions = hLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kPion);
 637   AliDebug(1, Form("BinPions %d", binPions));
 638   Int_t binProtons =  hLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kProton);
 639   AliDebug(1, Form("BinProtons %d", binProtons));
 640   hLikeTRD->GetAxis(kSpecies)->SetRange(binElectrons, binElectrons);
 641   TH2 *likeElectron = hLikeTRD->Projection(kElectronLike, kP);
 642   likeElectron->SetName("likeElectron");
 643   hLikeTRD->GetAxis(kSpecies)->SetRange(binPions, binPions);
 644   TH2 *likePion = hLikeTRD->Projection(kElectronLike, kP);
 645   likePion->SetName("likePion");
 646   hLikeTRD->GetAxis(kSpecies)->SetRange(binProtons, binProtons);
 647   TH2 *likeProton = hLikeTRD->Projection(kElectronLike, kP);
 648   likeProton->SetName("likeProton");
 649
 650   // Undo ranges
 651   hLikeTRD->GetAxis(kSpecies)->SetRange(0, hLikeTRD->GetAxis(kSpecies)->GetNbins());
 652   hLikeTRD->GetAxis(kNTracklets)->SetRange(0, hLikeTRD->GetAxis(kNTracklets)->GetNbins());
 653   hLikeTRD->GetAxis(kNTracklets)->SetRange(0, hLikeTRD->GetAxis(kCentralityBinLike)->GetNbins());
 654
 655   // Prepare List for output
 656   Char_t *listname=Form("%dTracklets", nTracklets);
 657   if(isPbPb) listname=Form("%dTracklets%dCentrality", nTracklets, nCentrality);
 658
 659
 660   TList *listPions = new TList; listPions->SetName(listname); listPions->SetOwner();
 661   TList *listProtons = new TList; listProtons->SetName(listname); listProtons->SetOwner();
 662   TList *listThresholds = new TList; listThresholds->SetName(listname); listThresholds->SetOwner();
 663   fPionEfficiencies->Add(listPions);
 664   fProtonEfficiencies->Add(listProtons);
 665   fThresholds->Add(listThresholds);
 666
 667   TH1 *probsEl = NULL, *probsPi = NULL, *probsPr = NULL;
 668   TGraphErrors *effPi = NULL, *effPr = NULL, *thresholds = NULL;
 669   Double_t p = 0, dp = 0;
 670   Int_t threshbin = 0;
 671   Double_t eff, error; // value and error
 672   for(Int_t ieff = 0; ieff < kNElectronEffs; ieff++){
 673
 674     if(fShowMessage){
 675       printf("-----------------------------------------\n");
 676       printf("Doing Electron Efficiency %f\n", fgkElectronEff[ieff]);
 677       printf("-----------------------------------------\n");
 678     }
 679     effPi = new TGraphErrors(likeElectron->GetXaxis()->GetNbins());
 680     effPi->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
 681     effPr = new TGraphErrors(likeElectron->GetXaxis()->GetNbins());
 682     effPr->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
 683     thresholds = new TGraphErrors(likeElectron->GetXaxis()->GetNbins());
 684     thresholds->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
 685
 686     // Add to lists
 687     listPions->Add(effPi);
 688     listProtons->Add(effPr);
 689     listThresholds->Add(thresholds);
 690
 691     for(Int_t imom = 1; imom <= likeElectron->GetXaxis()->GetLast(); imom++){
 692       p = likeElectron->GetXaxis()->GetBinCenter(imom);
 693       dp = likeElectron->GetXaxis()->GetBinWidth(imom)/2;
 694
 695       probsEl = likeElectron->ProjectionY("el", imom, imom);
 696       if(!probsEl->GetEntries()) continue;
 697       probsPi = likePion->ProjectionY("pi", imom, imom);
 698       if(!probsPi->GetEntries()) continue;
 699       probsPr = likeProton->ProjectionY("pr", imom, imom);
 700       if(!probsPr->GetEntries()) continue;
 701       AliDebug(1, Form("Calculating Values for p = %f", p));
 702
 703       // Calculate non-electronEfficiency and error
 704       eff = CalculateHadronEfficiency(probsPi, probsEl, fgkElectronEff[ieff], threshbin, error);
 705       thresholds->SetPoint(imom - 1, p, probsEl->GetXaxis()->GetBinCenter(threshbin));
 706       thresholds->SetPointError(imom - 1, dp, EstimateThresholdError(probsEl, threshbin));
 707       AliDebug(1, Form("threshold %d|%f", threshbin, probsEl->GetXaxis()->GetBinCenter(threshbin)));
 708       AliDebug(1, Form("Pion Efficiency %f +- %f", eff, error));
 709       effPi->SetPoint(imom - 1, p, eff);
 710       effPi->SetPointError(imom - 1, dp, error);
 711       eff = CalculateHadronEfficiency(probsPr, probsEl, fgkElectronEff[ieff] , threshbin, error);
 712       AliDebug(1, Form("Proton Efficiency %f", eff));
 713       effPr->SetPoint(imom - 1, p, eff);
 714       effPr->SetPointError(imom - 1, dp, error);
 715
 716       // cleanup
 717       delete probsEl;
 718       delete probsPi;
 719       delete probsPr;
 720     }
 721   }
 722
 723   // remove temporary histograms
 724   delete likeElectron;
 725   delete likePion;
 726   delete likeProton;
 727 }
 728
 729 //__________________________________________________________________
 730 Double_t AliHFEtrdPIDqa::CalculateHadronEfficiency(const TH1 * const hadron, const TH1 *const electron, Double_t eff, Int_t &threshbin, Double_t &error){
 731   //
 732   // Calculate non-electron efficiency
 733   // optionally returns sums as second parameter
 734   //
 735
 736   TH1D eleWorking(*((const TH1D *)electron)), hadronWorking(*((const TH1D *)hadron)); // Leave the original histograms untouched and do calculation including scale on copy
 737   eleWorking.Scale(1./eleWorking.Integral());
 738   hadronWorking.Scale(1./hadronWorking.Integral());
 739
 740   TArrayD sumsEl(eleWorking.GetNbinsX()), sumsHd(eleWorking.GetNbinsX());
 741
 742   // calculate threshold and estimated electron efficiency the threshold was taken
 743   Double_t elEff = 0.;  // estimated electron efficiency at the end
 744   Int_t currentBin = 0, nbins = 0;
 745   for(Int_t ibin = eleWorking.GetXaxis()->GetLast(); ibin >= eleWorking.GetXaxis()->GetFirst(); ibin--){
 746     currentBin = ibin;
 747     nbins++;
 748     elEff += eleWorking.GetBinContent(ibin);
 749     sumsEl[eleWorking.GetXaxis()->GetLast() - ibin] = elEff;
 750     if(elEff >= eff){
 751       // we found the matching bin, break the loop
 752       break;
 753     }
 754   }
 755   threshbin = currentBin;
 756
 757   Double_t hdEff = 0;
 758   for(Int_t ibin = hadronWorking.GetXaxis()->GetLast(); ibin >= threshbin; ibin--) {
 759     hdEff += hadronWorking.GetBinContent(ibin);
 760     sumsHd[hadronWorking.GetXaxis()->GetLast() - ibin] = hdEff;
 761   }
 762
 763   // search sums of electron efficiency for double counts, eliminate in electron and hadron array
 764   TArrayD newsumsEl(100), newsumsHd(100);
 765   Int_t nusable = 0;
 766   for(Int_t ien = 0; ien < nbins; ien++){
 767     if(ien==0){
 768       newsumsEl[0] = sumsEl[0];
 769       nusable++;
 770       continue;
 771     }
 772     Int_t index = TMath::BinarySearch(nusable, newsumsEl.GetArray(), sumsEl[ien]);
 773     if(TMath::Abs(sumsEl[ien] - newsumsEl[index]) < 1e-13){
 774       // element already counted, don't add to the new arrays
 775       continue;
 776     }
 777     newsumsEl[nusable] = sumsEl[ien];
 778     newsumsHd[nusable] = sumsHd[ien];
 779     nusable++;
 780   }
 781
 782   //printf("New array\n");
 783   //for(Int_t ib = 0; ib < nusable; ib++){
 784   //  printf("Electron Efficiency %f, Pion Efficiency %f\n", newsumsEl[ib], newsumsHd[ib]);
 785   //}
 786   //printf("Do Fit\n");
 787
 788   // Calculate error
 789   error = 0;
 790   if(hadronWorking.GetEntries() > 0 && eleWorking.GetEntries() > 0 && nusable > 2){
 791     // Do error calculation in case the bins have enough statistics
 792     TGraph gevh(nusable, newsumsEl.GetArray(), newsumsHd.GetArray());
 793     TF1 evh("evh","pol2", eff-.05, eff+.05);
 794     gevh.Fit(&evh, "Q", "", eff-.05, eff+.05);
 795
 796     // return the error of the pion efficiency
 797     if(((1.-hdEff) < 0) || ((1.- elEff) < 0)){
 798       AliError(" ElEffi or HdEffi > 1. Error can not be calculated. Please increase statistics!");
 799     }   else {
 800       error = TMath::Sqrt(hdEff*(1-hdEff)/hadronWorking.GetEntries()+TMath::Power(evh.Derivative(eff), 2)*elEff*(1-elEff)/eleWorking.GetEntries());
 801     }
 802     AliDebug(2, Form("Pion Effi at [%f] : [%f +/- %f], Threshold[%f]", elEff, hdEff, error, eleWorking.GetBinCenter(threshbin)));
 803     AliDebug(2, Form("Derivative at %4.2f : %f\n", eff, evh.Derivative(eff)));
 804   }
 805
 806   return hdEff;
 807 }
 808
 809 //__________________________________________________________________
 810 Double_t AliHFEtrdPIDqa::EstimateThresholdError(const TH1 * const electron, Int_t threshbin){
 811   //
 812   // Estimate threshold error as sqrt(N_int)/N_ent
 813   // where N_int are the counts integrated under the electron selection band
 814   // and N_ent are all entries in the histogram
 815   //
 816   Int_t nInt = 0;
 817   for(Int_t ibin = threshbin; ibin <= electron->GetXaxis()->GetNbins(); ibin++){
 818     nInt += (Int_t) electron->GetBinContent(ibin);
 819   }
 820   Int_t nEnt = (Int_t) electron->Integral();
 821   if(nEnt == 0) return 1.;
 822   return TMath::Sqrt(static_cast<Double_t>(nInt))/static_cast<Double_t>(nEnt);
 823 }
 824
 825 //__________________________________________________________________
 826 Double_t AliHFEtrdPIDqa::CalculateIntegratedPionEfficiency(UInt_t nTracklets, Double_t electronEff, Double_t pmin, Double_t pmax, Int_t icentrality, Double_t *error){
 827   //
 828   // Calculate Pion Efficiency for a given electron efficiency in the specified momentum range
 829   //
 830   if(nTracklets < 4 || nTracklets > 6){
 831     AliError("Pion Efficiency calculation only available for 4, 5, and 6 tracklets");
 832     return 0.;
 833   }
 834   if(electronEff < 0.6 || electronEff > 1.){
 835     AliError("Pion Efficiency calculation only available in the electron efficiency range 0.6 to 1");
 836     return 0.;
 837   }
 838   if(pmin < 0.1 || pmin > 10 || pmax < 0.1 || pmax > 10.){
 839     AliError("Pion Efficiency calculation only available in the momentum range 0.1 to 10 GeV/c");
 840     return 0.;
 841   }
 842   if(pmax < pmin){
 843     AliError("pmin is expected to be >= pmax");
 844     return 0.;
 845   }
 846
 847   // prerequierements fullfiled
 848   // prepare histos
 849   THnSparse *hLikeTRD = dynamic_cast<THnSparseF *>(fHistos->Get("fLikeTRD"));
 850   if(!hLikeTRD){
 851     AliError("Likelihood Histogram not available");
 852     return 0;
 853   }
 854   Int_t binTracklets = hLikeTRD->GetAxis(kNTracklets)->FindBin(nTracklets);
 855   hLikeTRD->GetAxis(kNTracklets)->SetRange(binTracklets, binTracklets);
 856
 857   if(icentrality!=1){
 858       Int_t binCentrality = hLikeTRD->GetAxis(kCentralityBinLike)->FindBin(icentrality);
 859       hLikeTRD->GetAxis(kCentralityBinLike)->SetRange(binCentrality, binCentrality);
 860   }
 861
 862   Int_t pbinMin = hLikeTRD->GetAxis(kP)->FindBin(pmax),
 863         pbinMax = hLikeTRD->GetAxis(kP)->FindBin(pmax);
 864   hLikeTRD->GetAxis(kP)->SetRange(pbinMin, pbinMax);
 865   Int_t binElectrons = hLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kElectron);
 866   Int_t binPions = hLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kPion);
 867   hLikeTRD->GetAxis(kSpecies)->SetRange(binElectrons, binElectrons);
 868   TH1 *likeElectron = hLikeTRD->Projection(kElectronLike);
 869   likeElectron->Scale(1./likeElectron->Integral());
 870   likeElectron->SetName("likeElectron");
 871   hLikeTRD->GetAxis(kSpecies)->SetRange(binPions, binPions);
 872   TH1 *likePion = hLikeTRD->Projection(kElectronLike);
 873   likePion->Scale(1./likePion->Integral());
 874   likePion->SetName("likePion");
 875
 876   // Undo ranges
 877   hLikeTRD->GetAxis(kSpecies)->SetRange(0, hLikeTRD->GetAxis(kSpecies)->GetNbins());
 878   hLikeTRD->GetAxis(kNTracklets)->SetRange(0, hLikeTRD->GetAxis(kNTracklets)->GetNbins());
 879   hLikeTRD->GetAxis(kP)->SetRange(0, hLikeTRD->GetAxis(kP)->GetNbins());
 880
 881   // Do Calculation
 882   Int_t thresh; Double_t err;
 883   Double_t effpi = CalculateHadronEfficiency(likePion, likeElectron, electronEff, thresh, err);
 884   delete likePion; delete likeElectron;
 885   if(error) *error = err;
 886   return effpi;
 887 }
 888
 889 //__________________________________________________________________
 890 void AliHFEtrdPIDqa::DrawTracklet(Int_t itracklet, Double_t pmin, Double_t pmax, Int_t icentrality, Bool_t doFit){
 891   //
 892   // Draw efficiencies and threshold as function of p
 893   //
 894   if(!(fPionEfficiencies && fProtonEfficiencies && fThresholds)){
 895     AliError("No graphs to draw available");
 896     return;
 897   }
 898
 899   Char_t *listname=Form("%dTracklets", itracklet);
 900   if(icentrality!=-1) listname=Form("%dTracklets%dCentrality", itracklet, icentrality);
 901
 902   TList *lpions = dynamic_cast<TList *>(fPionEfficiencies->FindObject(listname));
 903   TList *lprotons = dynamic_cast<TList *>(fProtonEfficiencies->FindObject(listname));
 904
 905   TList *lthresholds = dynamic_cast<TList *>(fThresholds->FindObject(listname));
 906   if(!(lpions && lprotons && lthresholds)){
 907     AliDebug(1, "Relevant lists not found. Did you forget to run FinishAnalysis()?");
 908     return;
 909   }
 910
 911   TGraphErrors *pi, *pr;
 912   TGraph *tr;
 913   TLegend *leg;
 914   Char_t *canvasname=Form("tracklet%d", itracklet);
 915   if(icentrality!=-1) canvasname=Form("tracklet%dcentrality%d", itracklet, icentrality);
 916   TCanvas *c1 = new TCanvas(canvasname, canvasname, 1024, 768);
 917   c1->Divide(3,2);
 918   TF1 *threshfit = NULL;
 919   for(Int_t ieff = 0; ieff < kNElectronEffs; ieff++){
 920     c1->cd(ieff + 1);
 921     gPad->SetGrid(0,0);
 922     gPad->SetLeftMargin(0.12);
 923     gPad->SetRightMargin(0.08);
 924     pi = dynamic_cast<TGraphErrors *>(lpions->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
 925     pr = dynamic_cast<TGraphErrors *>(lprotons->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
 926     tr = dynamic_cast<TGraph *>(lthresholds->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
 927     if(!(pi && pr && tr)) continue;
 928
 929     // Define nice plot, draw
 930     // Axis Title
 931     pi->GetXaxis()->SetTitle("p / GeV/c");
 932     pi->GetYaxis()->SetTitle("Efficiency");
 933     pr->GetXaxis()->SetTitle("p / GeV/c");
 934     pr->GetYaxis()->SetTitle("Efficiency");
 935     tr->GetXaxis()->SetTitle("p / GeV/c");
 936     tr->GetYaxis()->SetTitle("Efficiency");
 937     pi->GetYaxis()->SetTitleOffset(1.2);
 938     pr->GetYaxis()->SetTitleOffset(1.2);
 939     tr->GetYaxis()->SetTitleOffset(1.2);
 940     pi->GetXaxis()->SetTitleSize(0.045);
 941     pi->GetYaxis()->SetTitleSize(0.045);
 942     pr->GetXaxis()->SetTitleSize(0.045);
 943     pr->GetYaxis()->SetTitleSize(0.045);
 944     tr->GetXaxis()->SetTitleSize(0.045);
 945     tr->GetYaxis()->SetTitleSize(0.045);
 946     // Axis Range
 947     pi->GetYaxis()->SetRangeUser(1e-3, 1.);
 948     pr->GetYaxis()->SetRangeUser(1e-3, 1.);
 949     tr->GetYaxis()->SetRangeUser(1e-3, 1.);
 950     if(pmin >= 0 && pmax >= 0.){
 951       pi->GetXaxis()->SetRangeUser(pmin, pmax);
 952       pr->GetXaxis()->SetRangeUser(pmin, pmax);
 953       tr->GetXaxis()->SetRangeUser(pmin, pmax);
 954     }
 955     // Marker
 956     pi->SetMarkerColor(kRed);
 957     pi->SetMarkerStyle(20);
 958     pr->SetMarkerColor(kBlue);
 959     pr->SetMarkerStyle(21);
 960     tr->SetMarkerColor(kBlack);
 961     tr->SetMarkerStyle(22);
 962     // Title
 963     pi->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
 964     pr->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
 965     tr->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
 966     // Draw
 967     pi->Draw("ape"); pr->Draw("pesame"); tr->Draw("pesame");
 968
 969     // Optionally do Fit
 970     if(doFit){
 971       threshfit = MakeThresholds(tr, pmin, pmax);
 972       threshfit->SetLineColor(kBlack);
 973       threshfit->Draw("same");
 974     }
 975
 976     // Add entries to legend
 977     if(ieff==0){
 978       leg = new TLegend(0.5, 0.65, 0.89, 0.85);
 979       leg->SetBorderSize(0);
 980       leg->SetFillStyle(0);
 981       leg->AddEntry(pi, "Pion Efficiency", "lp");
 982       leg->AddEntry(pr, "Proton Efficiency", "lp");
 983       leg->AddEntry(tr, "Thresholds", "lp");
 984       leg->Draw();
 985       gPad->Update();
 986     }
 987   }
 988   c1->cd();
 989 }
 990
 991 //__________________________________________________________________
 992 TF1 *AliHFEtrdPIDqa::MakeThresholds(TGraph *threshist, Double_t lowerLimit, Double_t upperLimit){
 993   //
 994   // Create TF1 containing the threshold parametrisation
 995   //
 996
 997   TF1 *threshparam = new TF1("thresh", "1-[0]-[1]*x-[2]*TMath::Exp(-[3]*x)", 0.1, 10);
 998   threshist->Fit(threshparam, "NE", "", lowerLimit, upperLimit);
 999   return threshparam;
1000 }
1001
1002 //__________________________________________________________________
1003 void AliHFEtrdPIDqa::ClearLists(){
1004   //
1005   // Clear lists for particle efficiencies and thresholds
1006   //
1007   if(fPionEfficiencies){
1008     fPionEfficiencies->Delete();
1009     delete fPionEfficiencies;
1010     fPionEfficiencies = NULL;
1011   }
1012   if(fProtonEfficiencies){
1013     fProtonEfficiencies->Delete();
1014     delete fProtonEfficiencies;
1015     fProtonEfficiencies = NULL;
1016   }
1017   if(fThresholds){
1018     fThresholds->Delete();
1019     delete fThresholds;
1020     fThresholds = NULL;
1021   }
1022 }
1023
1024 //__________________________________________________________________
1025 Double_t AliHFEtrdPIDqa::EvalPionEfficiency(Int_t ntls, Int_t eEff, Double_t p, Int_t icentrality){
1026   //
1027   // calculate pion efficiency
1028   // Arguments:
1029   //   Number of tracklets
1030   //   Electron Efficiency
1031   //   Momentum
1032   //
1033   TGraphErrors *measurement = GetPionEfficiency(ntls, eEff, icentrality);
1034   if(!measurement) return -1.;
1035   return measurement->Eval(p);
1036 }
1037
1038 //__________________________________________________________________
1039 Double_t AliHFEtrdPIDqa::EvalProtonEfficiency(Int_t ntls, Int_t eEff, Double_t p, Int_t icentrality){
1040   //
1041   // calculate proton efficiency
1042   // Arguments:
1043   //   Number of tracklets
1044   //   Electron Efficiency
1045   //   Momentum
1046   //
1047   TGraphErrors *measurement = GetProtonEfficiency(ntls, eEff, icentrality);
1048   if(!measurement) return -1.;
1049   return measurement->Eval(p);
1050 }
1051
1052 //__________________________________________________________________
1053 Double_t AliHFEtrdPIDqa::EvalThreshold(Int_t ntls, Int_t eEff, Double_t p, Int_t icentrality){
1054   //
1055   // Get the threshold to reach a certain electron efficiency
1056   // Arguments:
1057   //   Number of tracklets
1058   //   Electron Efficiency
1059   //   Momentum
1060   //
1061   TGraph *measurement = GetThreshold(ntls, eEff, icentrality);
1062   if(!measurement) return -1.;
1063   return measurement->Eval(p);
1064 }
1065
1066 //__________________________________________________________________
1067 TGraphErrors *AliHFEtrdPIDqa::GetPionEfficiency(Int_t ntracklets, Int_t eleffpercent, Int_t icentrality){
1068   //
1069   // Get Graph with pion efficiencies
1070   //
1071   Char_t *listname=Form("%dTracklets", ntracklets);
1072   if(icentrality!=-1) listname=Form("%dTracklets%dCentrality", ntracklets, icentrality);
1073   TList *graphs = dynamic_cast<TList *>(fPionEfficiencies->FindObject(listname));
1074   if(!graphs) return NULL;
1075   return dynamic_cast<TGraphErrors *>(graphs->FindObject(Form("eff%d", eleffpercent)));
1076 }
1077
1078 //__________________________________________________________________
1079 TGraphErrors *AliHFEtrdPIDqa::GetProtonEfficiency(Int_t ntracklets, Int_t eleffpercent, Int_t icentrality){
1080   //
1081   // Get Graph with proton efficiencies
1082   //
1083   Char_t *listname=Form("%dTracklets", ntracklets);
1084   if(icentrality!=-1) listname=Form("%dTracklets%dCentrality", ntracklets, icentrality);
1085   TList *graphs = dynamic_cast<TList *>(fProtonEfficiencies->FindObject(listname));
1086   if(!graphs) return NULL;
1087   return dynamic_cast<TGraphErrors *>(graphs->FindObject(Form("eff%d", eleffpercent)));
1088 }
1089
1090 //__________________________________________________________________
1091 TGraph *AliHFEtrdPIDqa::GetThreshold(Int_t ntracklets, Int_t eleffpercent, Int_t icentrality){
1092   //
1093   // Get Graph with threshols
1094   //
1095   Char_t *listname=Form("%dTracklets", ntracklets);
1096   if(icentrality!=-1) listname=Form("%dTracklets%dCentrality", ntracklets, icentrality);
1097   TList *graphs = dynamic_cast<TList *>(fThresholds->FindObject(listname));
1098   if(!graphs) return NULL;
1099   return dynamic_cast<TGraph *>(graphs->FindObject(Form("eff%d", eleffpercent)));
1100 }
1101