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 <TH2.h>
  33 #include <THnSparse.h>
  34 #include <TString.h>
  35
  36 #include "AliLog.h"
  37 #include "AliPID.h"
  38
  39 #include "AliHFEcollection.h"
  40 #include "AliHFEpidBase.h"
  41 #include "AliHFEpidQAmanager.h"
  42 #include "AliHFEpidTPC.h"
  43 #include "AliHFEpidTRD.h"
  44 #include "AliHFEtrdPIDqaV1.h"
  45
  46 ClassImp(AliHFEtrdPIDqaV1)
  47
  48 //____________________________________________________________
  49 AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1():
  50     AliHFEdetPIDqa(),
  51     fHistos(NULL)
  52 {
  53   //
  54   // Dummy constructor
  55   //
  56 }
  57
  58 //____________________________________________________________
  59 AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1(const Char_t *name):
  60     AliHFEdetPIDqa(name, "QA for TRD"),
  61     fHistos(NULL)
  62 {
  63   //
  64   // Default constructor
  65   //
  66 }
  67
  68 //____________________________________________________________
  69 AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1(const AliHFEtrdPIDqaV1 &o):
  70     AliHFEdetPIDqa(o),
  71     fHistos(NULL)
  72 {
  73   //
  74   // Copy constructor
  75   //
  76 }
  77
  78 //____________________________________________________________
  79 AliHFEtrdPIDqaV1 &AliHFEtrdPIDqaV1::operator=(const AliHFEtrdPIDqaV1 &o){
  80   //
  81   // Make assignment
  82   //
  83   AliHFEdetPIDqa::operator=(o);
  84   fHistos = o.fHistos;
  85
  86   return *this;
  87 }
  88
  89 //_________________________________________________________
  90 Long64_t AliHFEtrdPIDqaV1::Merge(TCollection *coll){
  91   //
  92   // Merge with other objects
  93   //
  94   if(!coll) return 0;
  95   if(coll->IsEmpty()) return 1;
  96
  97   TIter it(coll);
  98   AliHFEtrdPIDqaV1 *refQA = NULL;
  99   TObject *o = NULL;
 100   Long64_t count = 0;
 101   TList listHistos;
 102   while((o = it())){
 103     refQA = dynamic_cast<AliHFEtrdPIDqaV1 *>(o);
 104     if(!refQA) continue;
 105
 106     listHistos.Add(refQA->fHistos);
 107     count++;
 108   }
 109   fHistos->Merge(&listHistos);
 110   return count + 1;
 111 }
 112
 113 //____________________________________________________________
 114 void AliHFEtrdPIDqaV1::Initialize(){
 115   //
 116   // Initialize QA histos for TRD PID
 117   //
 118
 119   AliDebug(1, "Initializing PID QA for TRD");
 120   // Make common binning
 121   const Int_t kPIDbins = AliPID::kSPECIES + 1;
 122   const Int_t kPbins = 100;
 123   const Int_t kSteps = 2;
 124   const Double_t kMinPID = -1;
 125   const Double_t kMinP = 0.;
 126   const Double_t kMaxPID = (Double_t)AliPID::kSPECIES;
 127   const Double_t kMaxP = 20.;
 128
 129   fHistos = new AliHFEcollection("trdqahistos", "Collection of TRD QA histograms");
 130
 131   // Create Control Histogram monitoring the TPC sigma between the selection steps
 132   const Int_t kTPCSigmaBins = 140;
 133   Int_t nBinsTPCSigma[4] = {kPIDbins, kPbins, kTPCSigmaBins, kSteps};
 134   Double_t minTPCSigma[4] = {kMinPID, kMinP, -12., 0};
 135   Double_t maxTPCSigma[4] = {kMaxPID, kMaxP, 12., 2.};
 136   fHistos->CreateTHnSparse("hTPCsigma", "TPC sigma; species p [GeV/c]; TPC dEdx - <dE/dx>|_{el} [#sigma]; selection step", 4, nBinsTPCSigma, minTPCSigma, maxTPCSigma);
 137   // Create Monitoring histogram for the Likelihood distribution
 138   const Int_t kTRDLikelihoodBins = 100;
 139   Int_t nBinsTRDlike[4] = {kPIDbins, kPbins, kTRDLikelihoodBins, kSteps};
 140   Double_t minTRDlike[4] = {kMinPID, kMinP, 0., 0.};
 141   Double_t maxTRDlike[4] = {kMaxPID, kMaxP, 1., 2.};
 142   fHistos->CreateTHnSparse("hTRDlikelihood", "TRD Likelihood Distribution; species; p [GeV/c]; TRD electron Likelihood; selection step", 4, nBinsTRDlike, minTRDlike, maxTRDlike);
 143   // Create Monitoring histogram for the TRD total charge
 144   const Int_t kTRDchargeBins = 1000;
 145   Int_t nBinsTRDcharge[4] = {kPIDbins, kPbins, kTRDchargeBins, kSteps};
 146   Double_t minTRDcharge[4] = {kMinPID, kMinP, 0., 0.};
 147   Double_t maxTRDcharge[4] = {kMaxPID, kMaxP, 100000., 2.};
 148   fHistos->CreateTHnSparse("hTRDcharge", "Total TRD charge; species; p [GeV/c]; TRD charge [a.u.]; selection step", 4, nBinsTRDcharge, minTRDcharge, maxTRDcharge);
 149 }
 150
 151 //____________________________________________________________
 152 void AliHFEtrdPIDqaV1::ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step){
 153   //
 154   // Process the track, fill the containers
 155   //
 156   AliDebug(1, Form("QA started for TRD PID for step %d", (Int_t)step));
 157   Int_t species = track->GetAbInitioPID();
 158   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
 159
 160   AliHFEpidTRD *trdpid = dynamic_cast<AliHFEpidTRD *>(fQAmanager->GetDetectorPID(AliHFEpid::kTRDpid));
 161   AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fQAmanager->GetDetectorPID(AliHFEpid::kTPCpid));
 162
 163   Double_t container[4];
 164   container[0] = species;
 165   container[1] = trdpid->GetP(track->GetRecTrack(), anatype);
 166   container[2] = tpcpid->NumberOfSigmas(track->GetRecTrack(), AliPID::kElectron, anatype);
 167   container[3] = step;
 168   fHistos->Fill("hTPCsigma", container);
 169
 170   container[2] = trdpid->GetElectronLikelihood(track->GetRecTrack(), anatype);
 171   fHistos->Fill("hTRDlikelihood", container);
 172
 173   for(UInt_t ily = 0; ily < 6; ily++){
 174     container[2] = trdpid->GetChargeLayer(track->GetRecTrack(), ily, anatype);
 175     fHistos->Fill("hTRDcharge", container);
 176   }
 177  if(species >= AliPID::kSPECIES) species = -1;
 178 }
 179
 180 //_________________________________________________________
 181 TH2 *AliHFEtrdPIDqaV1::MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t species){
 182   //
 183   // Get the TPC control histogram for the TRD selection step (either before or after PID)
 184   //
 185   printf("histos :%p\n", fHistos);
 186   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTPCsigma"));
 187   if(!histo){
 188     AliError("QA histogram monitoring TPC nSigma not available");
 189     return NULL;
 190   }
 191   if(species > -1 && species < AliPID::kSPECIES){
 192     // cut on species (if available)
 193     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
 194   }
 195   histo->GetAxis(3)->SetRange(step + 1, step + 1);
 196
 197   TH2 *hSpec = histo->Projection(2, 1);
 198   // construct title and name
 199   TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
 200   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
 201   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
 202   TString histname = Form("hSigmaTPC%s%s", specID.Data(), stepname.Data());
 203   TString histtitle = Form("TPC Sigma for %s %s PID", speciesname.Data(), stepname.Data());
 204   hSpec->SetName(histname.Data());
 205   hSpec->SetTitle(histtitle.Data());
 206
 207   // Unset range on the original histogram
 208   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
 209   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
 210   return hSpec;
 211 }
 212
 213 //_________________________________________________________
 214 TH2 *AliHFEtrdPIDqaV1::MakeTRDlikelihoodDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species){
 215   //
 216   // Make Histogram for TRD Likelihood distribution
 217   //
 218   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDlikelihood"));
 219   if(!histo){
 220     AliError("QA histogram monitoring TRD Electron Likelihood not available");
 221     return NULL;
 222   }
 223   if(species > -1 && species < AliPID::kSPECIES){
 224     // cut on species (if available)
 225     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
 226   }
 227   histo->GetAxis(3)->SetRangeUser(step + 1, step + 1);
 228
 229   TH2 *hSpec = histo->Projection(2, 1);
 230   // construct title and name
 231   TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
 232   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
 233   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
 234   TString histname = Form("hLikeElTRD%s%s", specID.Data(), stepname.Data());
 235   TString histtitle = Form("TRD electron Likelihood for %s %s PID", speciesname.Data(), stepname.Data());
 236   hSpec->SetName(histname.Data());
 237   hSpec->SetTitle(histtitle.Data());
 238
 239   // Unset range on the original histogram
 240   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
 241   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
 242   return hSpec;
 243 }
 244
 245 //_________________________________________________________
 246 TH2 *AliHFEtrdPIDqaV1::MakeTRDchargeDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species){
 247   //
 248   // Make Histogram for TRD Likelihood distribution
 249   //
 250   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDcharge"));
 251   if(!histo){
 252     AliError("QA histogram monitoring TRD total charge not available");
 253     return NULL;
 254   }
 255   if(species > -1 && species < AliPID::kSPECIES){
 256     // cut on species (if available)
 257     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
 258   }
 259   histo->GetAxis(3)->SetRange(step + 1, step + 1);
 260
 261   TH2 *hSpec = histo->Projection(2, 1);
 262   // construct title and name
 263   TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
 264   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
 265   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
 266   TString histname = Form("hChargeTRD%s%s", specID.Data(), stepname.Data());
 267   TString histtitle = Form("TRD total charge for %s %s PID", speciesname.Data(), stepname.Data());
 268   hSpec->SetName(histname.Data());
 269   hSpec->SetTitle(histtitle.Data());
 270
 271   // Unset range on the original histogram
 272   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
 273   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
 274   return hSpec;
 275 }
 276