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