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