PWG3/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 {
  83   //
  84   // Default Constructor
  85   //
  86 }
  87
  88 //__________________________________________________________________
  89 AliHFEtrdPIDqa::AliHFEtrdPIDqa(const Char_t *name):
  90   TNamed(name, ""),
  91   fTRDpid(NULL),
  92   fHistos(NULL),
  93   fPionEfficiencies(NULL),
  94   fProtonEfficiencies(NULL),
  95   fKaonEfficiencies(NULL),
  96   fThresholds(NULL)
  97 {
  98   //
  99   // Main Constructor
 100   //
 101 }
 102
 103 //__________________________________________________________________
 104 AliHFEtrdPIDqa::AliHFEtrdPIDqa(const AliHFEtrdPIDqa &ref):
 105   TNamed(ref),
 106   fTRDpid(NULL),
 107   fHistos(NULL),
 108   fPionEfficiencies(NULL),
 109   fProtonEfficiencies(NULL),
 110   fKaonEfficiencies(NULL),
 111   fThresholds(NULL)
 112 {
 113   //
 114   // Copy constructor
 115   //
 116   ref.Copy(*this);
 117 }
 118
 119 //__________________________________________________________________
 120 AliHFEtrdPIDqa &AliHFEtrdPIDqa::operator=(const AliHFEtrdPIDqa &ref){
 121   //
 122   // Assignment operator
 123   //
 124   if(this != &ref)
 125     ref.Copy(*this);
 126   return *this;
 127 }
 128
 129 //__________________________________________________________________
 130 AliHFEtrdPIDqa::~AliHFEtrdPIDqa(){
 131   //
 132   // Destructor
 133   //
 134   if(fTRDpid) delete fTRDpid;
 135   if(fHistos) delete fHistos;
 136   if(fPionEfficiencies) delete fPionEfficiencies;
 137   if(fProtonEfficiencies) delete fProtonEfficiencies;
 138   if(fKaonEfficiencies) delete fKaonEfficiencies;
 139 }
 140
 141 //__________________________________________________________________
 142 void AliHFEtrdPIDqa::Copy(TObject &ref) const{
 143   //
 144   // Copies content of this object into object ref
 145   //
 146   TNamed::Copy(ref);
 147
 148   AliHFEtrdPIDqa &target = dynamic_cast<AliHFEtrdPIDqa &>(ref);
 149   target.fTRDpid = fTRDpid;
 150   target.fHistos = dynamic_cast<AliHFEcollection *>(fHistos->Clone());
 151 }
 152
 153 //__________________________________________________________________
 154 Long64_t AliHFEtrdPIDqa::Merge(TCollection *coll){
 155   //
 156   // Merge objects
 157   //
 158   if(!coll) return 0;
 159   if(coll->IsEmpty()) return 1;
 160
 161   AliHFEtrdPIDqa *refQA = NULL;
 162   TIter it(coll);
 163   TObject *o = NULL;
 164   Long64_t count = 0;
 165   TList listHistos;
 166   while((o = it())){
 167     refQA = dynamic_cast<AliHFEtrdPIDqa *>(o);
 168     if(!refQA) continue;
 169
 170     listHistos.Add(refQA->fHistos);
 171     count++;
 172   }
 173   fHistos->Merge(&listHistos);
 174   return count+1;
 175 }
 176
 177 //__________________________________________________________________
 178 void AliHFEtrdPIDqa::Browse(TBrowser *b){
 179   //
 180   // Enable Browser functionality
 181   //
 182   if(b){
 183     // Add objects to the browser
 184     if(fHistos) b->Add(fHistos, fHistos->GetName());
 185     if(fPionEfficiencies) b->Add(fPionEfficiencies, "Pion Efficiencies");
 186     if(fProtonEfficiencies) b->Add(fProtonEfficiencies, "Proton Efficiencies");
 187     if(fKaonEfficiencies) b->Add(fKaonEfficiencies, "Kaon Efficiencies");
 188     if(fThresholds) b->Add(fThresholds, "Thresholds");
 189   }
 190 }
 191
 192 //__________________________________________________________________
 193 void AliHFEtrdPIDqa::Init(){
 194   //
 195   // Initialize Object
 196   //
 197
 198   fHistos = new AliHFEcollection("TRDqa", "Histos for TRD QA");
 199
 200   CreateLikelihoodHistogram();
 201   CreateQAHistogram();
 202   CreatedEdxHistogram();
 203   CreateHistoTruncatedMean();
 204
 205   fTRDpid = new AliHFEpidTRD("QAtrdPID");
 206 }
 207
 208 //__________________________________________________________________
 209 void AliHFEtrdPIDqa::CreateLikelihoodHistogram(){
 210   //
 211   // Create Histogram for TRD Likelihood Studies
 212   //
 213   Int_t nbins[kQuantitiesLike]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
 214   nbins[kElectronLike] = 100;
 215   Double_t binMin[kQuantitiesLike]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
 216   Double_t binMax[kQuantitiesLike]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
 217   binMin[kElectronLike] = 0.;
 218   binMax[kElectronLike] = 1.;
 219
 220   fHistos->CreateTHnSparse("fLikeTRD","TRD Likelihood Studies", kQuantitiesLike, nbins, binMin, binMax);
 221   fHistos->BinLogAxis("fLikeTRD", kP);
 222 }
 223
 224 //__________________________________________________________________
 225 void AliHFEtrdPIDqa::CreateQAHistogram(){
 226   //
 227   // Create Histogram for Basic TRD PID QA
 228   //
 229   AliDebug(1, "Called");
 230   Int_t nbins[kQuantitiesQA]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
 231   nbins[kNonZeroTrackletCharge] = AliESDtrack::kTRDnPlanes + 1;
 232   nbins[kNClusters] = 200;
 233   Double_t binMin[kQuantitiesQA]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
 234   binMin[kNonZeroTrackletCharge] = 0.;
 235   binMin[kNClusters] = 0.;
 236   Double_t binMax[kQuantitiesQA]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
 237   binMax[kNonZeroTrackletCharge] = AliESDtrack::kTRDnPlanes + 1.;
 238   binMax[kNClusters] = 200.;
 239
 240   fHistos->CreateTHnSparse("fQAtrack","TRD QA Histogram", kQuantitiesQA, nbins, binMin, binMax);
 241   fHistos->BinLogAxis("fQAtrack", kP);
 242 }
 243
 244 //__________________________________________________________________
 245 void AliHFEtrdPIDqa::CreatedEdxHistogram(){
 246   //
 247   // Create QA histogram for dEdx investigations
 248   //
 249   AliDebug(1, "Called");
 250   Int_t nbins[kQuantitiesdEdx]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
 251   nbins[kdEdx] = 100;
 252   nbins[kNclusters] = 261;
 253   nbins[kNonZeroSlices] = 9;
 254   Double_t binMin[kQuantitiesdEdx]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
 255   binMin[kdEdx] = 0.;
 256   binMin[kNclusters] = 0;
 257   binMin[kNonZeroSlices] = 0.;
 258   Double_t binMax[kQuantitiesdEdx]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
 259   binMax[kdEdx] = 100000.;
 260   binMax[kNclusters] = 260.;
 261   binMax[kNonZeroSlices] = 8.;
 262
 263   fHistos->CreateTHnSparse("fQAdEdx","TRD summed dEdx", kQuantitiesdEdx, nbins, binMin, binMax);
 264   fHistos->BinLogAxis("fQAdEdx", kP);
 265   fHistos->Sumw2("fQAdEdx");
 266 }
 267
 268 //__________________________________________________________________
 269 void AliHFEtrdPIDqa::CreateHistoTruncatedMean(){
 270   //
 271   // Create Histogram for Basic TRD PID QA
 272   //
 273   AliDebug(1, "Called");
 274   Int_t nbins[kQuantitiesTruncMean]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
 275   nbins[kTPCdEdx] = 600;
 276   nbins[kTRDdEdxMethod1] = 1000;
 277   nbins[kTRDdEdxMethod2] = 1000;
 278   Double_t binMin[kQuantitiesTruncMean]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
 279   binMin[kTPCdEdx] = 0.;
 280   binMin[kTRDdEdxMethod1] = 0.;
 281   binMin[kTRDdEdxMethod2] = 0.;
 282   Double_t binMax[kQuantitiesTruncMean]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
 283   binMax[kTPCdEdx] = 600;
 284   binMax[kTRDdEdxMethod1] = 20000.;
 285   binMax[kTRDdEdxMethod2] = 20000.;
 286
 287   fHistos->CreateTHnSparse("fTRDtruncMean","TRD TruncatedMean studies", kQuantitiesTruncMean, nbins, binMin, binMax);
 288   fHistos->BinLogAxis("fTRDtruncMean", kP);
 289   fHistos->CreateTH2F("fTRDslicesPions","TRD dEdx per slice for Pions", 8, 0, 8, 500, 0, 2000);
 290   fHistos->CreateTH2F("fTRDslicesElectrons","TRD dEdx per slice for Electrons", 8, 0, 8, 500, 0, 2000);
 291 }
 292
 293
 294 //__________________________________________________________________
 295 void AliHFEtrdPIDqa::ProcessTracks(TObjArray * const tracks, Int_t species){
 296   //
 297   // Process Electron Tracks
 298   //
 299   if(species < -1 || species >= AliPID::kSPECIES) return;
 300   TIterator *it = tracks->MakeIterator();
 301   AliVTrack *track = NULL;
 302   while((track = dynamic_cast<AliVTrack *>(it->Next()))){
 303     if(track) ProcessTrack(track, species);
 304   }
 305   delete it;
 306 }
 307
 308 //__________________________________________________________________
 309 void AliHFEtrdPIDqa::ProcessTrack(AliVTrack *track, Int_t species){
 310   //
 311   // Process Single Track
 312   //
 313   if(TString(track->IsA()->GetName()).CompareTo("AliESDtrack") == 0)
 314     ProcessTrackESD(dynamic_cast<AliESDtrack *>(track), species);
 315   else if(TString(track->IsA()->GetName()).CompareTo("AliAODTrack") == 0)
 316     ProcessTrackAOD(dynamic_cast<AliAODTrack *>(track), species);
 317 }
 318
 319
 320 //__________________________________________________________________
 321 void AliHFEtrdPIDqa::ProcessTrackESD(AliESDtrack *track, Int_t species){
 322   //
 323   // Process single ESD track
 324   //
 325   if(!track) return;
 326   if((track->GetStatus() & AliESDtrack::kTRDout) == 0) return;  // require TRD track
 327   FillTRDLikelihoods(track, species);
 328   FillTRDQAplots(track, species);
 329 }
 330
 331 //__________________________________________________________________
 332 void AliHFEtrdPIDqa::ProcessTrackAOD(AliAODTrack * const track, Int_t /*species*/){
 333   //
 334   // Process single AOD track
 335   // AOD PID object required
 336   //
 337   if(!track) return;
 338   AliAODPid *trackPID = track->GetDetPid();
 339   if(!trackPID) return;
 340
 341 }
 342
 343 //__________________________________________________________________
 344 void AliHFEtrdPIDqa::FillTRDLikelihoods(AliESDtrack *track, Int_t species){
 345   //
 346   // Fill TRD Likelihood Histogram
 347   //
 348   Double_t trdLike[AliPID::kSPECIES];
 349   track->GetTRDpid(trdLike);
 350   const AliExternalTrackParam *outerPars = track->GetOuterParam();
 351
 352   Double_t quantities[kQuantitiesLike]; memset(quantities, 0, sizeof(Double_t) * kQuantitiesLike);
 353   // we store:
 354   // species
 355   // p
 356   // ntracklets
 357   // Electron Likelihood
 358   quantities[kSpecies] = species;
 359   quantities[kP] = outerPars ? outerPars->P() : track->P();
 360   quantities[kNTracklets] = track->GetTRDntrackletsPID();
 361   quantities[kElectronLike] = trdLike[AliPID::kElectron];
 362
 363   fHistos->Fill("fLikeTRD", quantities);
 364 }
 365
 366 //__________________________________________________________________
 367 void AliHFEtrdPIDqa::FillTRDQAplots(AliESDtrack *track, Int_t species){
 368   //
 369   // Fill QA Plots containing further information
 370   //
 371   const AliExternalTrackParam *outerPars = track->GetOuterParam();
 372
 373   Double_t quantitiesQA[kQuantitiesQA], quantitiesdEdx[kQuantitiesdEdx], quantitiesTruncMean[kQuantitiesTruncMean];
 374   // we store:
 375   // 1. QA
 376   // species
 377   // p
 378   // ntracklets
 379   // Non-zero tracklet charges
 380   // Number of clusters / full track
 381   // 2. dEdx
 382   // species
 383   // p
 384   // ntracklets
 385   // dEdx
 386   // 3. Truncated Mean
 387   // ...
 388   // TPC dEdx
 389   // TRD dEdx Method 1
 390   // TRD dEdx Method 2
 391   quantitiesQA[kSpecies]  = quantitiesdEdx[kSpecies]
 392                           = quantitiesTruncMean[kSpecies]
 393                           = species;
 394   quantitiesQA[kP]  = quantitiesTruncMean[kP]
 395                     = outerPars ? outerPars->P() : track->P();
 396   quantitiesQA[kNTracklets] = quantitiesdEdx[kNTracklets]
 397                             = quantitiesTruncMean[kNTracklets]
 398                             = track->GetTRDntrackletsPID();
 399   quantitiesQA[kNClusters] = quantitiesdEdx[kNclusters] = track->GetTRDncls();
 400
 401
 402   Double_t dEdxSum = 0., qSlice = 0.;
 403   // remove the last slice from the histogram
 404   Int_t ntrackletsNonZero = 0, nSlices = track->GetNumberOfTRDslices(), nSlicesNonZero = 0;
 405   TString speciesname = "pions";
 406   Bool_t selectedForSlicemon = kFALSE;
 407
 408   switch(species){
 409     case AliPID::kElectron: speciesname = "Electrons"; selectedForSlicemon = kTRUE; break;
 410     case AliPID::kPion: speciesname = "Pions"; selectedForSlicemon = kTRUE; break;
 411     default: speciesname = "undefined"; selectedForSlicemon = kFALSE; break;
 412   };
 413   AliDebug(1, Form("species %d, speciesname %s, momentum %f, selected %s", species, speciesname.Data(), track->P(), selectedForSlicemon ? "yes" : "no"));
 414   for(Int_t iplane = 0; iplane < AliESDtrack::kTRDnPlanes; iplane++){
 415     quantitiesdEdx[kP] = track->GetTRDmomentum(iplane);
 416     dEdxSum = 0.;
 417     for(Int_t islice = 0; islice < nSlices; islice++){
 418       qSlice = track->GetTRDslice(iplane, islice);
 419       if(qSlice > 1e-1){
 420         // cut out 0 slices
 421         nSlicesNonZero++;
 422         dEdxSum += qSlice;
 423         // Reweighting of the slices for the truncated mean: select all pion tracks above
 424         // 1.5 GeV and monitor the dEdx as function of slice
 425         if(selectedForSlicemon && track->P() > 1.5){
 426           AliDebug(2, Form("Filling Histogram fTRDslices%s", speciesname.Data()));
 427           fHistos->Fill(Form("fTRDslices%s", speciesname.Data()), static_cast<Double_t>(islice), qSlice);
 428         }
 429       }
 430     }
 431     quantitiesdEdx[kNonZeroSlices] = nSlicesNonZero;
 432     quantitiesdEdx[kdEdx] = dEdxSum;
 433     if(dEdxSum) ntrackletsNonZero++;
 434     // Fill dEdx histogram
 435     if(dEdxSum > 1e-1) fHistos->Fill("fQAdEdx", quantitiesdEdx); // Cut out 0 entries
 436   }
 437   quantitiesQA[kNonZeroTrackletCharge] = ntrackletsNonZero;
 438   fHistos->Fill("fQAtrack", quantitiesQA);
 439
 440   quantitiesTruncMean[kTPCdEdx] = track->GetTPCsignal();
 441   quantitiesTruncMean[kTRDdEdxMethod1] = fTRDpid->GetTRDSignalV1(track, 0.6);
 442   quantitiesTruncMean[kTRDdEdxMethod2] = fTRDpid->GetTRDSignalV2(track, 0.6);
 443   fHistos->Fill("fTRDtruncMean", quantitiesTruncMean);
 444 }
 445
 446 /////////////////////////////////////////////////////////
 447 //
 448 // Code for Post Processing
 449 //
 450 // //////////////////////////////////////////////////////
 451
 452 //__________________________________________________________________
 453 void AliHFEtrdPIDqa::FinishAnalysis(){
 454   //
 455   // Finish Analysis:
 456   // Calculate Electron Efficiency for ntracklets = 4...6
 457   // Calculate thresholds for ntracklets = 4...6
 458   //
 459
 460   if(!fPionEfficiencies){
 461     fPionEfficiencies = new TList;
 462     fPionEfficiencies->SetName("pionEfficiencies");
 463   }
 464   if(!fProtonEfficiencies){
 465     fProtonEfficiencies = new TList;
 466     fProtonEfficiencies->SetName("protonEfficiencies");
 467   }
 468   if(!fThresholds){
 469     fThresholds = new TList;
 470     fThresholds->SetName("thresholds");
 471   }
 472
 473   for(Int_t itr = 4; itr <= 6; itr++){
 474     printf("========================================\n");
 475     printf("Analysing %d trackltes\n", itr);
 476     printf("========================================\n");
 477     AnalyseNTracklets(itr);
 478   }
 479 }
 480
 481 //__________________________________________________________________
 482 void AliHFEtrdPIDqa::StoreResults(const Char_t *filename){
 483   //
 484   // Store histos into a root file
 485   //
 486   TFile *outfile = new TFile(filename, "RECREATE");
 487   outfile->cd();
 488   fPionEfficiencies->Clone()->Write(fPionEfficiencies->GetName(), kSingleKey);
 489   fProtonEfficiencies->Clone()->Write(fProtonEfficiencies->GetName(), kSingleKey);
 490   fThresholds->Clone()->Write(fThresholds->GetName(), kSingleKey);
 491   outfile->Close();
 492   delete outfile;
 493 }
 494
 495 //__________________________________________________________________
 496 void AliHFEtrdPIDqa::SaveThresholdParameters(const Char_t *name){
 497   //
 498   // Fit the threshold histograms with the given parametrisation
 499   // and store the TF1 in the file
 500   //
 501
 502   if(!fThresholds){
 503     AliError("Threshold histograms have to be created first");
 504     return;
 505   }
 506
 507   printf("========================================\n");
 508   printf("Calculating threshold parameters\n");
 509   printf("========================================\n");
 510
 511   TList *outlist = new TList;
 512   outlist->SetName("thresholdTRD");
 513
 514   TGraph *threshhist = NULL;
 515   TF1 *threshparam = NULL;
 516   TList *lHistos = NULL, *lFormulas = NULL;
 517   for(Int_t itracklet = 4; itracklet <= 6; itracklet++){
 518
 519     printf("-------------------------------\n");
 520     printf("Processing %d tracklets\n", itracklet);
 521     printf("-------------------------------\n");
 522
 523     lHistos = dynamic_cast<TList *>(fThresholds->FindObject(Form("%dTracklets", itracklet)));
 524     if(!lHistos){
 525       AliError(Form("Threshold histograms for the case %d tracklets not found", itracklet));
 526       continue;
 527     }
 528     lFormulas = new TList;
 529     lFormulas->SetName(Form("%dTracklets", itracklet));
 530     outlist->Add(lFormulas);
 531
 532     for(Int_t ieff = 0; ieff <  kNElectronEffs; ieff++){
 533
 534       printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n");
 535       printf("Processing Electron Efficiency %f\n", fgkElectronEff[ieff]);
 536       printf(">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>\n");
 537
 538       threshhist = dynamic_cast<TGraph *>(lHistos->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
 539       if(!threshhist) continue;
 540       threshparam = MakeThresholds(threshhist);
 541       threshparam->SetName(Form("thresh_%d_%d", itracklet, static_cast<Int_t>(fgkElectronEff[ieff] * 100)));
 542       lFormulas->Add(threshparam);
 543     }
 544   }
 545
 546   // store the output
 547   TFile *outfile = new TFile(name, "RECREATE");
 548   outfile->cd();
 549   outlist->Write(outlist->GetName(), kSingleKey);
 550   outfile->Close();
 551   delete outfile;
 552 }
 553
 554 //__________________________________________________________________
 555 void AliHFEtrdPIDqa::AnalyseNTracklets(Int_t nTracklets){
 556   //
 557   // Calculate Pion Efficiency, Proton Efficiency and Kaon Efficiency at discrete
 558   // elPion Efficiency, Proton Efficiency and Kaon Efficiency at discrete
 559   // electron efficiencies
 560   //
 561   THnSparse *hLikeTRD = dynamic_cast<THnSparseF *>(fHistos->Get("fLikeTRD"));
 562   if(!hLikeTRD){
 563     AliError("Likelihood Histogram not available");
 564     return;
 565   }
 566   Int_t binTracklets = hLikeTRD->GetAxis(kNTracklets)->FindBin(nTracklets);
 567   hLikeTRD->GetAxis(kNTracklets)->SetRange(binTracklets, binTracklets);
 568
 569   Int_t binElectrons = hLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kElectron);
 570   AliDebug(1, Form("BinElectrons %d", binElectrons));
 571   Int_t binPions = hLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kPion);
 572   AliDebug(1, Form("BinPions %d", binPions));
 573   Int_t binProtons =  hLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kProton);
 574   AliDebug(1, Form("BinProtons %d", binProtons));
 575   hLikeTRD->GetAxis(kSpecies)->SetRange(binElectrons, binElectrons);
 576   TH2 *likeElectron = hLikeTRD->Projection(kElectronLike, kP);
 577   likeElectron->SetName("likeElectron");
 578   hLikeTRD->GetAxis(kSpecies)->SetRange(binPions, binPions);
 579   TH2 *likePion = hLikeTRD->Projection(kElectronLike, kP);
 580   likePion->SetName("likePion");
 581   hLikeTRD->GetAxis(kSpecies)->SetRange(binProtons, binProtons);
 582   TH2 *likeProton = hLikeTRD->Projection(kElectronLike, kP);
 583   likeProton->SetName("likeProton");
 584
 585   // Undo ranges
 586   hLikeTRD->GetAxis(kSpecies)->SetRange(0, hLikeTRD->GetAxis(kSpecies)->GetNbins());
 587   hLikeTRD->GetAxis(kNTracklets)->SetRange(0, hLikeTRD->GetAxis(kNTracklets)->GetNbins());
 588
 589   // Prepare List for output
 590   TList *listPions = new TList; listPions->SetName(Form("%dTracklets", nTracklets));
 591   TList *listProtons = new TList; listProtons->SetName(Form("%dTracklets", nTracklets));
 592   TList *listThresholds = new TList; listThresholds->SetName(Form("%dTracklets", nTracklets));
 593   fPionEfficiencies->Add(listPions);
 594   fProtonEfficiencies->Add(listProtons);
 595   fThresholds->Add(listThresholds);
 596
 597   TH1 *probsEl = NULL, *probsPi = NULL, *probsPr = NULL;
 598   TGraphErrors *effPi = NULL, *effPr = NULL; TGraph *thresholds = NULL;
 599   Double_t p = 0, dp = 0;
 600   Int_t threshbin = 0;
 601   Double_t noElEff[2]; // value and error
 602   for(Int_t ieff = 0; ieff < kNElectronEffs; ieff++){
 603     printf("-----------------------------------------\n");
 604     printf("Doing Electron Efficiency %f\n", fgkElectronEff[ieff]);
 605     printf("-----------------------------------------\n");
 606     effPi = new TGraphErrors(likeElectron->GetXaxis()->GetNbins());
 607     effPi->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
 608     effPr = new TGraphErrors(likeElectron->GetXaxis()->GetNbins());
 609     effPr->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
 610     thresholds = new TGraph(likeElectron->GetXaxis()->GetNbins());
 611     thresholds->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
 612
 613     // Add to lists
 614     listPions->Add(effPi);
 615     listProtons->Add(effPr);
 616     listThresholds->Add(thresholds);
 617
 618     for(Int_t imom = 1; imom <= likeElectron->GetXaxis()->GetLast(); imom++){
 619       p = likeElectron->GetXaxis()->GetBinCenter(imom);
 620       dp = likeElectron->GetXaxis()->GetBinWidth(imom)/2;
 621
 622       probsEl = likeElectron->ProjectionY("el", imom);
 623       if(!probsEl->GetEntries()) continue;
 624       probsEl->Scale(1./probsEl->Integral());
 625       probsPi = likePion->ProjectionY("pi", imom);
 626       if(!probsPi->GetEntries()) continue;
 627       probsPi->Scale(1./probsPi->Integral());
 628       probsPr = likeProton->ProjectionY("pr", imom);
 629       if(!probsPr->GetEntries()) continue;
 630       probsPr->Scale(1./probsPr->Integral());
 631       AliDebug(1, Form("Calculating Values for p = %f", p));
 632
 633       // Calculare threshold we need to achive the x% electron Efficiency
 634       threshbin = GetThresholdBin(probsEl, fgkElectronEff[ieff]);
 635       thresholds->SetPoint(imom - 1, p, probsEl->GetXaxis()->GetBinCenter(threshbin));
 636       AliDebug(1, Form("threshold %d|%f", threshbin, probsEl->GetXaxis()->GetBinCenter(threshbin)));
 637
 638       // Calculate non-electronEfficiency and error
 639       CalculateEfficiency(probsPi, threshbin, noElEff);
 640       AliDebug(1, Form("Pion Efficiency %f", noElEff[0]));
 641       effPi->SetPoint(imom - 1, p, noElEff[0]);
 642       effPi->SetPointError(imom - 1, dp, noElEff[1]);
 643       CalculateEfficiency(probsPr, threshbin, noElEff);
 644       effPr->SetPoint(imom - 1, p, noElEff[0]);
 645       effPr->SetPointError(imom - 1, dp, noElEff[1]);
 646       AliDebug(1, Form("Proton Efficiency %f", noElEff[0]));
 647
 648       // cleanup
 649       delete probsEl;
 650       delete probsPi;
 651       delete probsPr;
 652     }
 653   }
 654
 655   // remove temporary histograms
 656   delete likeElectron;
 657   delete likePion;
 658   delete likeProton;
 659 }
 660
 661 //__________________________________________________________________
 662 Int_t AliHFEtrdPIDqa::GetThresholdBin(TH1 * const input, Double_t eff){
 663   //
 664   // Calculate the threshold bin
 665   //
 666   Double_t integralEff = 0.;
 667   Int_t currentBin = 0;
 668   for(Int_t ibin = input->GetXaxis()->GetLast(); ibin >= input->GetXaxis()->GetFirst(); ibin--){
 669     currentBin = ibin;
 670     integralEff += input->GetBinContent(ibin);
 671     if(integralEff >= eff){
 672       // we found the matching bin, break the loop
 673       break;
 674     }
 675   }
 676   return currentBin;
 677 }
 678
 679 //__________________________________________________________________
 680 Bool_t AliHFEtrdPIDqa::CalculateEfficiency(TH1 * const input, Int_t threshbin, Double_t *par){
 681   //
 682   // Calculate non-electron efficiency
 683   //
 684   Double_t integralEff = 0;
 685   for(Int_t ibin = threshbin; ibin <= input->GetXaxis()->GetLast(); ibin++)
 686     integralEff += input->GetBinContent(ibin);
 687   par[0] = integralEff;
 688
 689   // @TODO: Error calculation
 690   par[1] = 0;
 691
 692   return kTRUE;
 693 }
 694
 695 //__________________________________________________________________
 696 void AliHFEtrdPIDqa::DrawTracklet(Int_t itracklet, Bool_t doFit){
 697   //
 698   // Draw efficiencies and threshold as function of p
 699   //
 700   if(!(fPionEfficiencies && fProtonEfficiencies && fThresholds)){
 701     AliError("No graphs to draw available");
 702     return;
 703   }
 704
 705   TList *lpions = dynamic_cast<TList *>(fPionEfficiencies->FindObject(Form("%dTracklets", itracklet)));
 706   TList *lprotons = dynamic_cast<TList *>(fProtonEfficiencies->FindObject(Form("%dTracklets", itracklet)));
 707
 708   TList *lthresholds = dynamic_cast<TList *>(fThresholds->FindObject(Form("%dTracklets", itracklet)));
 709   if(!(lpions && lprotons && lthresholds)){
 710     AliDebug(1, "Relevant lists not found. Did you forget to run FinishAnalysis()?");
 711     return;
 712   }
 713
 714   TGraphErrors *pi, *pr;
 715   TGraph *tr;
 716   TLegend *leg;
 717   TCanvas *c1 = new TCanvas(Form("tracklet%d", itracklet), Form("Tracklet %d", itracklet), 1024, 768);
 718   c1->Divide(3,2);
 719   TF1 *threshfit = NULL;
 720   for(Int_t ieff = 0; ieff < kNElectronEffs; ieff++){
 721     c1->cd(ieff + 1);
 722     leg = new TLegend(0.6, 0.7, 0.89, 0.89);
 723     leg->SetBorderSize(0);
 724     leg->SetFillStyle(0);
 725     pi = dynamic_cast<TGraphErrors *>(lpions->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
 726     pr = dynamic_cast<TGraphErrors *>(lprotons->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
 727     tr = dynamic_cast<TGraph *>(lthresholds->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
 728     if(!(pi && pr && tr)) continue;
 729
 730     // Define nice plot, draw
 731     // Axis Title
 732     pi->GetXaxis()->SetTitle("p / GeV/c");
 733     pi->GetYaxis()->SetTitle("Efficiency");
 734     pr->GetXaxis()->SetTitle("p / GeV/c");
 735     pr->GetYaxis()->SetTitle("Efficiency");
 736     tr->GetXaxis()->SetTitle("p / GeV/c");
 737     tr->GetYaxis()->SetTitle("Efficiency");
 738     // Axis Range
 739     pi->GetYaxis()->SetRangeUser(0., 1.);
 740     pr->GetYaxis()->SetRangeUser(0., 1.);
 741     tr->GetYaxis()->SetRangeUser(0., 1.);
 742     // Marker
 743     pi->SetMarkerColor(kRed);
 744     pi->SetMarkerStyle(20);
 745     pr->SetMarkerColor(kBlue);
 746     pr->SetMarkerStyle(21);
 747     tr->SetMarkerColor(kBlack);
 748     tr->SetMarkerStyle(22);
 749     // Title
 750     pi->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
 751     pr->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
 752     tr->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
 753     // Draw
 754     pi->Draw("ape"); pr->Draw("pesame"); tr->Draw("psame");
 755
 756     // Optionally do Fit
 757     if(doFit){
 758       threshfit = MakeThresholds(tr);
 759       threshfit->SetLineColor(kBlack);
 760       threshfit->Draw("same");
 761     }
 762
 763     // Add entries to legend
 764     leg->AddEntry(pi, "Pion Efficiency", "lp");
 765     leg->AddEntry(pr, "Proton Efficiency", "lp");
 766     leg->AddEntry(tr, "Thresholds", "lp");
 767     leg->Draw();
 768     c1->Update();
 769   }
 770 }
 771
 772 //__________________________________________________________________
 773 TF1 *AliHFEtrdPIDqa::MakeThresholds(TGraph *threshist){
 774   //
 775   // Create TF1 containing the threshold parametrisation
 776   //
 777
 778   TF1 *threshparam = new TF1("thresh", "1-[0]-[1]*x-[2]*TMath::Exp(-[3]*x)", 0.1, 10);
 779   threshist->Fit(threshparam, "NE", "", 0.1, 3.5);
 780   return threshparam;
 781 }
 782
 783 //__________________________________________________________________
 784 void AliHFEtrdPIDqa::ClearLists(){
 785   //
 786   // Clear lists for particle efficiencies and thresholds
 787   //
 788   if(fPionEfficiencies){
 789     fPionEfficiencies->Delete();
 790     delete fPionEfficiencies;
 791     fPionEfficiencies = NULL;
 792   }
 793   if(fProtonEfficiencies){
 794     fProtonEfficiencies->Delete();
 795     delete fProtonEfficiencies;
 796     fProtonEfficiencies = NULL;
 797   }
 798   if(fThresholds){
 799     fThresholds->Delete();
 800     delete fThresholds;
 801     fThresholds = NULL;
 802   }
 803 }