PWGHF/hfe/AliHFEtrdPIDqaV1.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 // Class AliHFEtrdPIDqaV1
  17 // Monitoring TRD PID in the HFE PID montioring framework. The following
  18 // quantities are monitored:
  19 //   TRD electron likelihood
  20 //   TRD dE/dx (Absolute values)
  21 //   TPC dE/dx (Number of sigmas, control histogram)
  22 // (Always as function of momentum, particle species and centrality
  23 // before and after cut)
  24 // More information about the PID monitoring framework can be found in
  25 // AliHFEpidQAmanager.cxx and AliHFEdetPIDqa.cxx
  26 //
  27 // Author:
  28 //    Markus Fasel <M.Fasel@gsi.de>
  29 //
  30
  31 #include <TAxis.h>
  32 #include <TBrowser.h>
  33 #include <TH2.h>
  34 #include <THnSparse.h>
  35 #include <TString.h>
  36
  37 #include "AliESDtrack.h"
  38 #include "AliLog.h"
  39 #include "AliPID.h"
  40 #include "AliPIDResponse.h"
  41
  42 #include "AliHFEcollection.h"
  43 #include "AliHFEpidBase.h"
  44 #include "AliHFEpidQAmanager.h"
  45 #include "AliHFEpidTPC.h"
  46 #include "AliHFEpidTRD.h"
  47 #include "AliHFEtrdPIDqaV1.h"
  48
  49 ClassImp(AliHFEtrdPIDqaV1)
  50
  51 //____________________________________________________________
  52 AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1():
  53     AliHFEdetPIDqa(),
  54     fHistos(NULL)
  55 {
  56   //
  57   // Dummy constructor
  58   //
  59 }
  60
  61 //____________________________________________________________
  62 AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1(const Char_t *name):
  63     AliHFEdetPIDqa(name, "QA for TRD"),
  64     fHistos(NULL)
  65 {
  66   //
  67   // Default constructor
  68   //
  69 }
  70
  71 //____________________________________________________________
  72 AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1(const AliHFEtrdPIDqaV1 &o):
  73     AliHFEdetPIDqa(o),
  74     fHistos(NULL)
  75 {
  76   //
  77   // Copy constructor
  78   //
  79 }
  80
  81 //____________________________________________________________
  82 AliHFEtrdPIDqaV1 &AliHFEtrdPIDqaV1::operator=(const AliHFEtrdPIDqaV1 &o){
  83   //
  84   // Make assignment
  85   //
  86   if(this == &o) return *this;
  87   AliHFEdetPIDqa::operator=(o);
  88   fHistos = o.fHistos;
  89
  90   return *this;
  91 }
  92 //____________________________________________________________
  93 AliHFEtrdPIDqaV1::~AliHFEtrdPIDqaV1(){
  94   //
  95   // Deconstructor
  96   //
  97   if(fHistos) delete fHistos;
  98
  99 }
 100 //_________________________________________________________
 101 Long64_t AliHFEtrdPIDqaV1::Merge(TCollection *coll){
 102   //
 103   // Merge with other objects
 104   //
 105   if(!coll) return 0;
 106   if(coll->IsEmpty()) return 1;
 107
 108   TIter it(coll);
 109   AliHFEtrdPIDqaV1 *refQA = NULL;
 110   TObject *o = NULL;
 111   Long64_t count = 0;
 112   TList listHistos;
 113   while((o = it())){
 114     refQA = dynamic_cast<AliHFEtrdPIDqaV1 *>(o);
 115     if(!refQA) continue;
 116
 117     listHistos.Add(refQA->fHistos);
 118     count++;
 119   }
 120   fHistos->Merge(&listHistos);
 121   return count + 1;
 122 }
 123
 124 //_________________________________________________________
 125 void AliHFEtrdPIDqaV1::Browse(TBrowser *b){
 126   //
 127   // Browse the PID QA
 128   //
 129   if(b){
 130     if(fHistos){
 131       b->Add(fHistos, fHistos->GetName());
 132
 133       // Make Projections of the dE/dx Spectra and add them to a new Folder
 134       TString specnames[4] = {"All", "Electrons", "Pions", "Protons"};
 135       Int_t specind[4] = {-1, AliPID::kElectron, AliPID::kPion, AliPID::kProton};
 136       TList *listTM = new TList;
 137       listTM->SetOwner();
 138       TList *listLike = new TList;
 139       listLike->SetOwner();
 140       TList *listCharge = new TList;
 141       listCharge->SetOwner();
 142       TList *listTPCnsigma = new TList;
 143       listTPCnsigma->SetOwner();
 144
 145       TH2 *hptr = NULL;
 146       TList *lctm, *lclike, *lccharge, *lcTPCsig;
 147       for(Int_t icent = -1; icent < 11; icent++){
 148         lctm = new TList;
 149         lctm->SetName(icent < 0 ? "MinBias" : Form("Centrality%d", icent));
 150         lctm->SetOwner();
 151         listTM->Add(lctm);
 152         lclike = new TList;
 153         lclike->SetName(icent < 0 ? "MinBias" : Form("Centrality%d", icent));
 154         lclike->SetOwner();
 155         listLike->Add(lclike);
 156         lccharge = new TList;
 157         lccharge->SetName(icent < 0 ? "MinBias" : Form("Centrality%d", icent));
 158         lccharge->SetOwner();
 159         listCharge->Add(lccharge);
 160         lcTPCsig = new TList;
 161         lcTPCsig->SetName(icent < 0 ? "MinBias" : Form("Centrality%d", icent));
 162         lcTPCsig->SetOwner();
 163         listTPCnsigma->Add(lcTPCsig);
 164         for(Int_t ispec = 0; ispec < 4; ispec++){
 165           for(Int_t istep = 0; istep < 2; istep++){
 166             hptr = MakeTRDspectrumTM(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], icent);
 167             if(hptr){
 168               hptr->SetName(Form("hTRDtm%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
 169               lctm->Add(hptr);
 170             }
 171             hptr = MakeTRDlikelihoodDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], icent);
 172             hptr->SetName(Form("hTRDlike%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
 173             lclike->Add(hptr);
 174             hptr = MakeTRDchargeDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], icent);
 175             hptr->SetName(Form("hTRDcharge%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
 176             lccharge->Add(hptr);
 177             hptr = MakeTPCspectrumNsigma(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], icent);
 178             hptr->SetName(Form("hTPCspectrum%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
 179             lcTPCsig->Add(hptr);
 180           }
 181         }
 182       }
 183
 184       b->Add(listTM, "Projections Truncated Mean");
 185       b->Add(listLike, "Projections Likelihood distribution");
 186       b->Add(listCharge, "Projections Tracklet Charge");
 187       b->Add(listTPCnsigma, "Projections TPC spectra");
 188     }
 189   }
 190 }
 191
 192 //____________________________________________________________
 193 void AliHFEtrdPIDqaV1::Initialize(){
 194   //
 195   // Initialize QA histos for TRD PID
 196   //
 197
 198   AliDebug(1, "Initializing PID QA for TRD");
 199   // Make common binning
 200   const Int_t kPIDbins = AliPID::kSPECIES + 1;
 201   const Int_t kSteps = 2;
 202   const Double_t kMinPID = -1;
 203   const Double_t kMinP = 0.;
 204   const Double_t kMaxPID = (Double_t)AliPID::kSPECIES;
 205   const Double_t kMaxP = 20.;
 206   const Int_t kCentralityBins = 11;
 207
 208   // Define number of bins
 209   Int_t kPbins = fQAmanager->HasHighResolutionHistos() ? 1000 : 100;
 210   Int_t tpcSigmaBins = fQAmanager->HasHighResolutionHistos() ? 1400 : 140;
 211   Int_t trdLikelihoodBins = fQAmanager->HasHighResolutionHistos() ? 200 : 100;
 212
 213   fHistos = new AliHFEcollection("trdqahistos", "Collection of TRD QA histograms");
 214
 215   // Create Control Histogram monitoring the TPC sigma between the selection steps
 216   Int_t nBinsTPCSigma[5] = {kPIDbins, kPbins, tpcSigmaBins, kSteps, kCentralityBins};
 217   Double_t minTPCSigma[5] = {kMinPID, kMinP, -12., 0., 0.};
 218   Double_t maxTPCSigma[5] = {kMaxPID, kMaxP, 12., 2., 11.};
 219   fHistos->CreateTHnSparse("hTPCsigma", "TPC sigma; species p [GeV/c]; TPC dEdx - <dE/dx>|_{el} [#sigma]; selection step", 5, nBinsTPCSigma, minTPCSigma, maxTPCSigma);
 220   fHistos->Sumw2("hTPCsigma");
 221   // Create Monitoring histogram for the Likelihood distribution
 222   Int_t nBinsTRDlike[5] = {kPIDbins, kPbins, trdLikelihoodBins, kSteps, kCentralityBins};
 223   Double_t minTRDlike[5] = {kMinPID, kMinP, 0., 0., 0.};
 224   Double_t maxTRDlike[5] = {kMaxPID, kMaxP, 1., 2., 11.};
 225   fHistos->CreateTHnSparse("hTRDlikelihood", "TRD Likelihood Distribution; species; p [GeV/c]; TRD electron Likelihood; selection step", 5, nBinsTRDlike, minTRDlike, maxTRDlike);
 226   fHistos->Sumw2("hTRDlikelihood");
 227   // Create Monitoring histogram for the TRD total charge
 228   const Int_t kTRDchargeBins = 1000;
 229   Int_t nBinsTRDcharge[5] = {kPIDbins, kPbins, kTRDchargeBins, kSteps, kCentralityBins};
 230   Double_t minTRDcharge[5] = {kMinPID, kMinP, 0., 0., 0.};
 231   Double_t maxTRDcharge[5] = {kMaxPID, kMaxP, 100000., 2., 11.};
 232   fHistos->CreateTHnSparse("hTRDcharge", "Total TRD charge; species; p [GeV/c]; TRD charge [a.u.]; selection step", 5, nBinsTRDcharge, minTRDcharge, maxTRDcharge);
 233   fHistos->Sumw2("hTRDcharge");
 234   // Monitoring of the TRD truncated mean according to version 1
 235   const Int_t kTRDtmBins = 1000;
 236   Int_t nBinsTRDtm[5] = {kPIDbins, kPbins, kTRDtmBins, kSteps, kCentralityBins};
 237   Double_t minTRDtm[5] = {kMinPID, kMinP, 0., 0., 0.};
 238   Double_t maxTRDtm[5] = {kMaxPID, kMaxP, 20000., 2., 11.};
 239   fHistos->CreateTHnSparse("hTRDtruncatedMean", "TRD truncated Mean; species; p [GeV/c]; TRD signal [a.u.]; selection step", 5, nBinsTRDtm, minTRDtm, maxTRDtm);
 240   fHistos->Sumw2("hTRDtruncatedMean");
 241
 242   // Monitoring of the number of tracklets
 243   fHistos->CreateTH2F("hNtrackletsBefore", "Number of tracklets before PID; species; p (GeV/c), Number of Tracklets", kPbins, kMinP, kMaxP, 7, 0., 7.);
 244   fHistos->Sumw2("hNtrackletsBefore");
 245   fHistos->CreateTH2F("hNtrackletsAfter", "Number of tracklets after PID; species; p (GeV/c), Number of Tracklets", kPbins, kMinP, kMaxP, 7, 0., 7.);
 246   fHistos->Sumw2("hNtrackletsAfter");
 247 }
 248
 249 //____________________________________________________________
 250 void AliHFEtrdPIDqaV1::ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step){
 251   //
 252   // Process the track, fill the containers
 253   //
 254   AliDebug(1, Form("QA started for TRD PID for step %d", (Int_t)step));
 255   Int_t species = track->GetAbInitioPID();
 256   if(species >= AliPID::kSPECIES) species = -1;
 257   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
 258
 259   AliHFEpidTRD *trdpid = dynamic_cast<AliHFEpidTRD *>(fQAmanager->GetDetectorPID(AliHFEpid::kTRDpid));
 260   const AliPIDResponse *pidResponse = trdpid ? trdpid->GetPIDResponse() : NULL;
 261
 262   Double_t container[5];
 263   container[0] = species;
 264   container[1] = trdpid ? trdpid->GetP(track->GetRecTrack(), anatype) : 0.;
 265   container[2] = pidResponse ? pidResponse->NumberOfSigmasTPC(track->GetRecTrack(), AliPID::kElectron) : 0.;
 266   container[3] = step;
 267   container[4] = track->GetCentrality();
 268   AliDebug(1, Form("Species %d, p %f, number of sigma TPC %f, step %f, centrality %f", (Int_t)species,(Float_t) container[1],(Float_t) container[2],(Float_t) container[3],(Float_t) container[4]));
 269   fHistos->Fill("hTPCsigma", container);
 270   AliDebug(1, "Filled  TPC sigma\n");
 271
 272
 273
 274   container[2] = trdpid ? trdpid->GetElectronLikelihood(static_cast<const AliVTrack*>(track->GetRecTrack()), anatype) : 0;
 275   fHistos->Fill("hTRDlikelihood", container);
 276   AliDebug(1, "Filled likelihood\n");
 277
 278   if(track->IsESDanalysis()){
 279     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track->GetRecTrack());
 280     if(esdtrack){
 281       container[2] = trdpid ? trdpid->GetTRDSignalV1(esdtrack) : 0;
 282       fHistos->Fill("hTRDtruncatedMean", container);
 283     }
 284   }
 285   for(UInt_t ily = 0; ily < 6; ily++){
 286     container[2] = trdpid ? trdpid->GetChargeLayer(track->GetRecTrack(), ily, anatype) : 0;
 287     if(container[2] < 1e-3) continue; // Filter out 0 entries
 288     fHistos->Fill("hTRDcharge", container);
 289     AliDebug(1, "Filled TRD charge\n");
 290   }
 291
 292   Int_t ntracklets = track->GetRecTrack()->GetTRDntrackletsPID();
 293   if(step == AliHFEdetPIDqa::kBeforePID)
 294     fHistos->Fill("hNtrackletsBefore", trdpid ? trdpid->GetP(track->GetRecTrack(), anatype) : 0., ntracklets);
 295   else
 296     fHistos->Fill("hNtrackletsAfter", trdpid ? trdpid->GetP(track->GetRecTrack(), anatype) : 0., ntracklets);
 297   AliDebug(1, "Filled tracklet\n");
 298 }
 299
 300 //_________________________________________________________
 301 TH2 *AliHFEtrdPIDqaV1::MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
 302   //
 303   // Get the TPC control histogram for the TRD selection step (either before or after PID)
 304   //
 305   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTPCsigma"));
 306   if(!histo){
 307     AliError("QA histogram monitoring TPC nSigma not available");
 308     return NULL;
 309   }
 310   if(species > -1 && species < AliPID::kSPECIES){
 311     // cut on species (if available)
 312     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
 313   }
 314   TString centname, centtitle;
 315   Bool_t hasCentralityInfo = kTRUE;
 316   if(centralityClass > -1){
 317     if(histo->GetNdimensions() < 5){
 318       AliError("Centrality Information not available");
 319       centname = centtitle = "MinBias";
 320       hasCentralityInfo = kFALSE;
 321     } else {
 322       // Project centrality classes
 323       // -1 is Min. Bias
 324       histo->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
 325       centname = Form("Cent%d", centralityClass);
 326       centtitle = Form("Centrality %d", centralityClass);
 327     }
 328   } else {
 329     histo->GetAxis(4)->SetRange(1, histo->GetAxis(4)->GetNbins()-1);
 330     centname = centtitle = "MinBias";
 331     hasCentralityInfo = kTRUE;
 332   }
 333   histo->GetAxis(3)->SetRange(step + 1, step + 1);
 334
 335   TH2 *hSpec = histo->Projection(2, 1);
 336   // construct title and name
 337   TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
 338   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
 339   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
 340   TString histname = Form("hSigmaTPC%s%s%s", specID.Data(), stepname.Data(), centname.Data());
 341   TString histtitle = Form("TPC Sigma for %s %s PID %s", speciesname.Data(), stepname.Data(), centtitle.Data());
 342   hSpec->SetName(histname.Data());
 343   hSpec->SetTitle(histtitle.Data());
 344
 345   // Unset range on the original histogram
 346   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
 347   histo->GetAxis(3)->SetRange(0, histo->GetAxis(3)->GetNbins());
 348   if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());
 349   return hSpec;
 350 }
 351
 352 //_________________________________________________________
 353 TH2 *AliHFEtrdPIDqaV1::MakeTRDspectrumTM(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
 354   //
 355   // Get the TPC control histogram for the TRD selection step (either before or after PID)
 356   //
 357   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDtruncatedMean"));
 358   if(!histo){
 359     AliError("QA histogram monitoring TPC nSigma not available");
 360     return NULL;
 361   }
 362   if(species > -1 && species < AliPID::kSPECIES){
 363     // cut on species (if available)
 364     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
 365   }
 366   TString centname, centtitle;
 367   Bool_t hasCentralityInfo = kTRUE;
 368   if(centralityClass > -1){
 369      if(histo->GetNdimensions() < 5){
 370       AliError("Centrality Information not available");
 371       centname = centtitle = "MinBias";
 372       hasCentralityInfo= kFALSE;
 373     } else {
 374       // Project centrality classes
 375       // -1 is Min. Bias
 376       histo->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
 377       centname = Form("Cent%d", centralityClass);
 378       centtitle = Form("Centrality %d", centralityClass);
 379     }
 380   } else {
 381     histo->GetAxis(4)->SetRange(1, histo->GetAxis(4)->GetNbins()-1);
 382     centname = centtitle = "MinBias";
 383     hasCentralityInfo= kFALSE;
 384   }
 385   histo->GetAxis(3)->SetRange(step + 1, step + 1);
 386
 387   TH2 *hSpec = histo->Projection(2, 1);
 388   // construct title and name
 389   TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
 390   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
 391   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
 392   TString histname = Form("hTMTRD%s%s%s", specID.Data(), stepname.Data(), centname.Data());
 393   TString histtitle = Form("TRD Truncated Mean for %s %s PID %s", speciesname.Data(), stepname.Data(), centtitle.Data());
 394   hSpec->SetName(histname.Data());
 395   hSpec->SetTitle(histtitle.Data());
 396
 397   // Unset range on the original histogram
 398   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
 399   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
 400   if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());
 401   return hSpec;
 402 }
 403
 404 //_________________________________________________________
 405 TH2 *AliHFEtrdPIDqaV1::MakeTRDlikelihoodDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
 406   //
 407   // Make Histogram for TRD Likelihood distribution
 408   //
 409   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDlikelihood"));
 410   if(!histo){
 411     AliError("QA histogram monitoring TRD Electron Likelihood not available");
 412     return NULL;
 413   }
 414   if(species > -1 && species < AliPID::kSPECIES){
 415     // cut on species (if available)
 416     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
 417   }
 418   TString centname, centtitle;
 419   Bool_t hasCentralityInfo = kTRUE;
 420   if(centralityClass > -1){
 421     if(histo->GetNdimensions() < 5){
 422       AliError("Centrality Information not available");
 423       centname = centtitle = "MinBias";
 424       hasCentralityInfo= kFALSE;
 425     } else {
 426       // Project centrality classes
 427       // -1 is Min. Bias
 428       histo->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
 429       centname = Form("Cent%d", centralityClass);
 430       centtitle = Form("Centrality %d", centralityClass);
 431     }
 432   } else {
 433     histo->GetAxis(4)->SetRange(1, histo->GetAxis(4)->GetNbins()-1);
 434     centname = centtitle = "MinBias";
 435     hasCentralityInfo= kTRUE;
 436   }
 437   histo->GetAxis(3)->SetRange(step + 1, step + 1);
 438
 439   TH2 *hSpec = histo->Projection(2, 1);
 440   // construct title and name
 441   TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
 442   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
 443   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
 444   TString histname = Form("hLikeElTRD%s%s%s", specID.Data(), stepname.Data(),centname.Data());
 445   TString histtitle = Form("TRD electron Likelihood for %s %s PID %s", speciesname.Data(), stepname.Data(),centtitle.Data());
 446   hSpec->SetName(histname.Data());
 447   hSpec->SetTitle(histtitle.Data());
 448
 449   // Unset range on the original histogram
 450   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
 451   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
 452   if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());
 453   return hSpec;
 454 }
 455
 456 //_________________________________________________________
 457 TH2 *AliHFEtrdPIDqaV1::MakeTRDchargeDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
 458   //
 459   // Make Histogram for TRD Likelihood distribution
 460   //
 461   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDcharge"));
 462   if(!histo){
 463     AliError("QA histogram monitoring TRD total charge not available");
 464     return NULL;
 465   }
 466   if(species > -1 && species < AliPID::kSPECIES){
 467     // cut on species (if available)
 468     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
 469   }
 470   TString centname, centtitle;
 471   Bool_t hasCentralityInfo = kTRUE;
 472   if(centralityClass > -1){
 473     if(histo->GetNdimensions() < 5){
 474       AliError("Centrality Information not available");
 475       centname = centtitle = "MinBias";
 476       hasCentralityInfo= kFALSE;
 477     } else {
 478       // Project centrality classes
 479       // -1 is Min. Bias
 480       histo->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
 481       centname = Form("Cent%d", centralityClass);
 482       centtitle = Form("Centrality %d", centralityClass);
 483     }
 484   } else {
 485     histo->GetAxis(4)->SetRange(1, histo->GetAxis(4)->GetNbins()-1);
 486     centname = centtitle = "MinBias";
 487     hasCentralityInfo= kTRUE;
 488   }
 489   histo->GetAxis(3)->SetRange(step + 1, step + 1);
 490
 491   TH2 *hSpec = histo->Projection(2, 1);
 492   // construct title and name
 493   TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
 494   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
 495   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
 496   TString histname = Form("hChargeTRD%s%s%s", specID.Data(), stepname.Data(), centname.Data());
 497   TString histtitle = Form("TRD total charge for %s %s PID %s", speciesname.Data(), stepname.Data(), centtitle.Data());
 498   hSpec->SetName(histname.Data());
 499   hSpec->SetTitle(histtitle.Data());
 500
 501   // Unset range on the original histogram
 502   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
 503   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
 504   if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());
 505   return hSpec;
 506 }
 507