PWG3/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
  41 #include "AliHFEcollection.h"
  42 #include "AliHFEpidBase.h"
  43 #include "AliHFEpidQAmanager.h"
  44 #include "AliHFEpidTPC.h"
  45 #include "AliHFEpidTRD.h"
  46 #include "AliHFEtrdPIDqaV1.h"
  47
  48 ClassImp(AliHFEtrdPIDqaV1)
  49
  50 //____________________________________________________________
  51 AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1():
  52     AliHFEdetPIDqa(),
  53     fHistos(NULL)
  54 {
  55   //
  56   // Dummy constructor
  57   //
  58 }
  59
  60 //____________________________________________________________
  61 AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1(const Char_t *name):
  62     AliHFEdetPIDqa(name, "QA for TRD"),
  63     fHistos(NULL)
  64 {
  65   //
  66   // Default constructor
  67   //
  68 }
  69
  70 //____________________________________________________________
  71 AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1(const AliHFEtrdPIDqaV1 &o):
  72     AliHFEdetPIDqa(o),
  73     fHistos(NULL)
  74 {
  75   //
  76   // Copy constructor
  77   //
  78 }
  79
  80 //____________________________________________________________
  81 AliHFEtrdPIDqaV1 &AliHFEtrdPIDqaV1::operator=(const AliHFEtrdPIDqaV1 &o){
  82   //
  83   // Make assignment
  84   //
  85   AliHFEdetPIDqa::operator=(o);
  86   fHistos = o.fHistos;
  87
  88   return *this;
  89 }
  90
  91 //_________________________________________________________
  92 Long64_t AliHFEtrdPIDqaV1::Merge(TCollection *coll){
  93   //
  94   // Merge with other objects
  95   //
  96   if(!coll) return 0;
  97   if(coll->IsEmpty()) return 1;
  98
  99   TIter it(coll);
 100   AliHFEtrdPIDqaV1 *refQA = NULL;
 101   TObject *o = NULL;
 102   Long64_t count = 0;
 103   TList listHistos;
 104   while((o = it())){
 105     refQA = dynamic_cast<AliHFEtrdPIDqaV1 *>(o);
 106     if(!refQA) continue;
 107
 108     listHistos.Add(refQA->fHistos);
 109     count++;
 110   }
 111   fHistos->Merge(&listHistos);
 112   return count + 1;
 113 }
 114
 115 //_________________________________________________________
 116 void AliHFEtrdPIDqaV1::Browse(TBrowser *b){
 117   //
 118   // Browse the PID QA
 119   //
 120   if(b){
 121     if(fHistos){
 122       b->Add(fHistos, fHistos->GetName());
 123
 124       // Make Projections of the dE/dx Spectra and add them to a new Folder
 125       TString specnames[4] = {"All", "Electrons", "Pions", "Protons"};
 126       Int_t specind[4] = {-1, AliPID::kElectron, AliPID::kPion, AliPID::kProton};
 127       TList *listTM = new TList;
 128       listTM->SetOwner();
 129       TList *listLike = new TList;
 130       listLike->SetOwner();
 131       TList *listCharge = new TList;
 132       listCharge->SetOwner();
 133       TList *listTPCnsigma = new TList;
 134       listTPCnsigma->SetOwner();
 135
 136       TH2 *hptr = NULL;
 137       for(Int_t ispec = 0; ispec < 4; ispec++){
 138         for(Int_t istep = 0; istep < 2; istep++){
 139           hptr = MakeTRDspectrumTM(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
 140           if(hptr){
 141             hptr->SetName(Form("hTRDtm%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
 142             listTM->Add(hptr);
 143           }
 144           hptr = MakeTRDlikelihoodDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
 145           hptr->SetName(Form("hTRDlike%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
 146           listLike->Add(hptr);
 147           hptr = MakeTRDchargeDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
 148           hptr->SetName(Form("hTRDcharge%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
 149           listCharge->Add(hptr);
 150           hptr = MakeTPCspectrumNsigma(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
 151           hptr->SetName(Form("hTPCspectrum%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
 152           listTPCnsigma->Add(hptr);
 153         }
 154       }
 155
 156       b->Add(listTM, "Projections Truncated Mean");
 157       b->Add(listLike, "Projections Likelihood distribution");
 158       b->Add(listCharge, "Projections Tracklet Charge");
 159       b->Add(listTPCnsigma, "Projections TPC spectra");
 160     }
 161   }
 162 }
 163
 164 //____________________________________________________________
 165 void AliHFEtrdPIDqaV1::Initialize(){
 166   //
 167   // Initialize QA histos for TRD PID
 168   //
 169
 170   AliDebug(1, "Initializing PID QA for TRD");
 171   // Make common binning
 172   const Int_t kPIDbins = AliPID::kSPECIES + 1;
 173   const Int_t kPbins = 100;
 174   const Int_t kSteps = 2;
 175   const Double_t kMinPID = -1;
 176   const Double_t kMinP = 0.;
 177   const Double_t kMaxPID = (Double_t)AliPID::kSPECIES;
 178   const Double_t kMaxP = 20.;
 179   const Int_t kCentralityBins = 11;
 180
 181   // Define number of bins
 182   Int_t tpcSigmaBins = fQAmanager->HasHighResolutionHistos() ? 1400 : 140;
 183   Int_t trdLikelihoodBins = fQAmanager->HasHighResolutionHistos() ? 200 : 100;
 184
 185   fHistos = new AliHFEcollection("trdqahistos", "Collection of TRD QA histograms");
 186
 187   // Create Control Histogram monitoring the TPC sigma between the selection steps
 188   Int_t nBinsTPCSigma[5] = {kPIDbins, kPbins, tpcSigmaBins, kSteps, kCentralityBins};
 189   Double_t minTPCSigma[5] = {kMinPID, kMinP, -12., 0., 0.};
 190   Double_t maxTPCSigma[5] = {kMaxPID, kMaxP, 12., 2., 11.};
 191   fHistos->CreateTHnSparse("hTPCsigma", "TPC sigma; species p [GeV/c]; TPC dEdx - <dE/dx>|_{el} [#sigma]; selection step", 4, nBinsTPCSigma, minTPCSigma, maxTPCSigma);
 192   fHistos->Sumw2("hTPCsigma");
 193   // Create Monitoring histogram for the Likelihood distribution
 194   Int_t nBinsTRDlike[5] = {kPIDbins, kPbins, trdLikelihoodBins, kSteps, kCentralityBins};
 195   Double_t minTRDlike[5] = {kMinPID, kMinP, 0., 0., 0.};
 196   Double_t maxTRDlike[5] = {kMaxPID, kMaxP, 1., 2., 11.};
 197   fHistos->CreateTHnSparse("hTRDlikelihood", "TRD Likelihood Distribution; species; p [GeV/c]; TRD electron Likelihood; selection step", 4, nBinsTRDlike, minTRDlike, maxTRDlike);
 198   fHistos->Sumw2("hTRDlikelihood");
 199   // Create Monitoring histogram for the TRD total charge
 200   const Int_t kTRDchargeBins = 1000;
 201   Int_t nBinsTRDcharge[5] = {kPIDbins, kPbins, kTRDchargeBins, kSteps, kCentralityBins};
 202   Double_t minTRDcharge[5] = {kMinPID, kMinP, 0., 0., 0.};
 203   Double_t maxTRDcharge[5] = {kMaxPID, kMaxP, 100000., 2., 11.};
 204   fHistos->CreateTHnSparse("hTRDcharge", "Total TRD charge; species; p [GeV/c]; TRD charge [a.u.]; selection step", 4, nBinsTRDcharge, minTRDcharge, maxTRDcharge);
 205   fHistos->Sumw2("hTRDcharge");
 206   // Monitoring of the TRD truncated mean according to version 1
 207   const Int_t kTRDtmBins = 1000;
 208   Int_t nBinsTRDtm[5] = {kPIDbins, kPbins, kTRDtmBins, kSteps, kCentralityBins};
 209   Double_t minTRDtm[5] = {kMinPID, kMinP, 0., 0., 0.};
 210   Double_t maxTRDtm[5] = {kMaxPID, kMaxP, 20000., 2., 11.};
 211   fHistos->CreateTHnSparse("hTRDtruncatedMean", "TRD truncated Mean; species; p [GeV/c]; TRD signal [a.u.]; selection step", 4, nBinsTRDtm, minTRDtm, maxTRDtm);
 212   fHistos->Sumw2("hTRDtruncatedMean");
 213 }
 214
 215 //____________________________________________________________
 216 void AliHFEtrdPIDqaV1::ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step){
 217   //
 218   // Process the track, fill the containers
 219   //
 220   AliDebug(1, Form("QA started for TRD PID for step %d", (Int_t)step));
 221   Int_t species = track->GetAbInitioPID();
 222   if(species >= AliPID::kSPECIES) species = -1;
 223   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
 224
 225   AliHFEpidTRD *trdpid = dynamic_cast<AliHFEpidTRD *>(fQAmanager->GetDetectorPID(AliHFEpid::kTRDpid));
 226   AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fQAmanager->GetDetectorPID(AliHFEpid::kTPCpid));
 227
 228   Double_t container[5];
 229   container[0] = species;
 230   container[1] = trdpid ? trdpid->GetP(track->GetRecTrack(), anatype) : 0.;
 231   container[2] = tpcpid ? tpcpid->NumberOfSigmas(track->GetRecTrack(), AliPID::kElectron, anatype) : 0.;
 232   container[3] = step;
 233   container[4] = track->GetCentrality();
 234   fHistos->Fill("hTPCsigma", container);
 235
 236   container[2] = trdpid ? trdpid->GetElectronLikelihood(track->GetRecTrack(), anatype) : 0;
 237   fHistos->Fill("hTRDlikelihood", container);
 238
 239   if(track->IsESDanalysis()){
 240     container[2] = trdpid ? trdpid->GetTRDSignalV1(dynamic_cast<const AliESDtrack *>(track->GetRecTrack())) : 0;
 241     fHistos->Fill("hTRDtruncatedMean", container);
 242   }
 243   for(UInt_t ily = 0; ily < 6; ily++){
 244     container[2] = trdpid ? trdpid->GetChargeLayer(track->GetRecTrack(), ily, anatype) : 0;
 245     fHistos->Fill("hTRDcharge", container);
 246   }
 247 }
 248
 249 //_________________________________________________________
 250 TH2 *AliHFEtrdPIDqaV1::MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t species){
 251   //
 252   // Get the TPC control histogram for the TRD selection step (either before or after PID)
 253   //
 254   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTPCsigma"));
 255   if(!histo){
 256     AliError("QA histogram monitoring TPC nSigma not available");
 257     return NULL;
 258   }
 259   if(species > -1 && species < AliPID::kSPECIES){
 260     // cut on species (if available)
 261     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
 262   }
 263   histo->GetAxis(3)->SetRange(step + 1, step + 1);
 264
 265   TH2 *hSpec = histo->Projection(2, 1);
 266   // construct title and name
 267   TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
 268   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
 269   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
 270   TString histname = Form("hSigmaTPC%s%s", specID.Data(), stepname.Data());
 271   TString histtitle = Form("TPC Sigma for %s %s PID", speciesname.Data(), stepname.Data());
 272   hSpec->SetName(histname.Data());
 273   hSpec->SetTitle(histtitle.Data());
 274
 275   // Unset range on the original histogram
 276   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
 277   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
 278   return hSpec;
 279 }
 280
 281 //_________________________________________________________
 282 TH2 *AliHFEtrdPIDqaV1::MakeTRDspectrumTM(AliHFEdetPIDqa::EStep_t step, Int_t species){
 283   //
 284   // Get the TPC control histogram for the TRD selection step (either before or after PID)
 285   //
 286   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDtruncatedMean"));
 287   if(!histo){
 288     AliError("QA histogram monitoring TPC nSigma not available");
 289     return NULL;
 290   }
 291   if(species > -1 && species < AliPID::kSPECIES){
 292     // cut on species (if available)
 293     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
 294   }
 295   histo->GetAxis(3)->SetRange(step + 1, step + 1);
 296
 297   TH2 *hSpec = histo->Projection(2, 1);
 298   // construct title and name
 299   TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
 300   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
 301   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
 302   TString histname = Form("hTMTRD%s%s", specID.Data(), stepname.Data());
 303   TString histtitle = Form("TRD Truncated Mean for %s %s PID", speciesname.Data(), stepname.Data());
 304   hSpec->SetName(histname.Data());
 305   hSpec->SetTitle(histtitle.Data());
 306
 307   // Unset range on the original histogram
 308   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
 309   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
 310   return hSpec;
 311 }
 312
 313 //_________________________________________________________
 314 TH2 *AliHFEtrdPIDqaV1::MakeTRDlikelihoodDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species){
 315   //
 316   // Make Histogram for TRD Likelihood distribution
 317   //
 318   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDlikelihood"));
 319   if(!histo){
 320     AliError("QA histogram monitoring TRD Electron Likelihood not available");
 321     return NULL;
 322   }
 323   if(species > -1 && species < AliPID::kSPECIES){
 324     // cut on species (if available)
 325     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
 326   }
 327   histo->GetAxis(3)->SetRange(step + 1, step + 1);
 328
 329   TH2 *hSpec = histo->Projection(2, 1);
 330   // construct title and name
 331   TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
 332   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
 333   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
 334   TString histname = Form("hLikeElTRD%s%s", specID.Data(), stepname.Data());
 335   TString histtitle = Form("TRD electron Likelihood for %s %s PID", speciesname.Data(), stepname.Data());
 336   hSpec->SetName(histname.Data());
 337   hSpec->SetTitle(histtitle.Data());
 338
 339   // Unset range on the original histogram
 340   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
 341   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
 342   return hSpec;
 343 }
 344
 345 //_________________________________________________________
 346 TH2 *AliHFEtrdPIDqaV1::MakeTRDchargeDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species){
 347   //
 348   // Make Histogram for TRD Likelihood distribution
 349   //
 350   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDcharge"));
 351   if(!histo){
 352     AliError("QA histogram monitoring TRD total charge not available");
 353     return NULL;
 354   }
 355   if(species > -1 && species < AliPID::kSPECIES){
 356     // cut on species (if available)
 357     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
 358   }
 359   histo->GetAxis(3)->SetRange(step + 1, step + 1);
 360
 361   TH2 *hSpec = histo->Projection(2, 1);
 362   // construct title and name
 363   TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
 364   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
 365   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
 366   TString histname = Form("hChargeTRD%s%s", specID.Data(), stepname.Data());
 367   TString histtitle = Form("TRD total charge for %s %s PID", speciesname.Data(), stepname.Data());
 368   hSpec->SetName(histname.Data());
 369   hSpec->SetTitle(histtitle.Data());
 370
 371   // Unset range on the original histogram
 372   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
 373   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
 374   return hSpec;
 375 }
 376