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 "AliAODTrack.h"
  37 #include "AliESDtrack.h"
  38 #include "AliESDpid.h"
  39 #include "AliExternalTrackParam.h"
  40 #include "AliPID.h"
  41
  42 #include "AliHFEcollection.h"
  43 #include "AliHFEpidBase.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(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   const AliVParticle *rectrack = track->GetRecTrack();
 159   if(species >= AliPID::kSPECIES) species = -1;
 160   if(!TString(rectrack->IsA()->GetName()).CompareTo("AliESDtrack")) ProcessESDtrack(dynamic_cast<const AliESDtrack *>(rectrack), step, species);
 161   else if(!TString(rectrack->IsA()->GetName()).CompareTo("AliAODTrack")) ProcessAODtrack(dynamic_cast<const AliAODTrack *>(rectrack), step, species);
 162   else  AliWarning(Form("Object type %s not supported\n", rectrack->IsA()->GetName()));
 163 }
 164
 165 //____________________________________________________________
 166 void AliHFEtrdPIDqaV1::ProcessESDtrack(const AliESDtrack *track, AliHFEdetPIDqa::EStep_t step, Int_t species){
 167   //
 168   // Fill QA histograms
 169   //
 170   Double_t container[4];
 171   container[0] = species;
 172   container[1] = track->GetOuterParam() ? track->GetOuterParam()->P() : track->P();
 173   container[2] = fESDpid->NumberOfSigmasTPC(track, AliPID::kElectron);
 174   container[3] = step;
 175   fHistos->Fill("hTPCsigma", container);
 176
 177   Double_t pid[5];
 178   track->GetTRDpid(pid);
 179   container[2] = pid[AliPID::kElectron];
 180   fHistos->Fill("hTRDlikelihood", container);
 181
 182   Double_t charge = 0.;
 183   for(Int_t ily = 0; ily < 6; ily++){
 184     charge = 0.;
 185     for(Int_t isl = 0; isl < track->GetNumberOfTRDslices(); isl++){
 186       charge += track->GetTRDslice(ily, isl);
 187     }
 188     fHistos->Fill("hTRDcharge", container);
 189   }
 190 }
 191
 192 //_________________________________________________________
 193 void AliHFEtrdPIDqaV1::ProcessAODtrack(const AliAODTrack * /*track*/, AliHFEdetPIDqa::EStep_t /*step*/, Int_t /*species*/){
 194   //
 195   // To be implemented
 196   //
 197 }
 198
 199 //_________________________________________________________
 200 TH2 *AliHFEtrdPIDqaV1::MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t species){
 201   //
 202   // Get the TPC control histogram for the TRD selection step (either before or after PID)
 203   //
 204   printf("histos :%p\n", fHistos);
 205   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTPCsigma"));
 206   if(!histo){
 207     AliError("QA histogram monitoring TPC nSigma not available");
 208     return NULL;
 209   }
 210   if(species > -1 && species < AliPID::kSPECIES){
 211     // cut on species (if available)
 212     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
 213   }
 214   histo->GetAxis(3)->SetRangeUser(species + 1, species + 1);
 215
 216   TH2 *hSpec = histo->Projection(2, 1);
 217   // construct title and name
 218   TString stepname = step ? "before" : "after";
 219   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
 220   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
 221   TString histname = Form("hSigmaTPC%s%s", specID.Data(), stepname.Data());
 222   TString histtitle = Form("TPC Sigma for %s %s PID", speciesname.Data(), stepname.Data());
 223   hSpec->SetName(histname.Data());
 224   hSpec->SetTitle(histtitle.Data());
 225
 226   // Unset range on the original histogram
 227   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
 228   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
 229   return hSpec;
 230 }
 231
 232 //_________________________________________________________
 233 TH2 *AliHFEtrdPIDqaV1::MakeTRDlikelihoodDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species){
 234   //
 235   // Make Histogram for TRD Likelihood distribution
 236   //
 237   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDlikelihood"));
 238   if(!histo){
 239     AliError("QA histogram monitoring TRD Electron Likelihood not available");
 240     return NULL;
 241   }
 242   if(species > -1 && species < AliPID::kSPECIES){
 243     // cut on species (if available)
 244     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
 245   }
 246   histo->GetAxis(3)->SetRangeUser(species + 1, species + 1);
 247
 248   TH2 *hSpec = histo->Projection(2, 1);
 249   // construct title and name
 250   TString stepname = step ? "before" : "after";
 251   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
 252   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
 253   TString histname = Form("hLikeElTRD%s%s", specID.Data(), stepname.Data());
 254   TString histtitle = Form("TRD electron Likelihood for %s %s PID", speciesname.Data(), stepname.Data());
 255   hSpec->SetName(histname.Data());
 256   hSpec->SetTitle(histtitle.Data());
 257
 258   // Unset range on the original histogram
 259   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
 260   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
 261   return hSpec;
 262 }
 263
 264 //_________________________________________________________
 265 TH2 *AliHFEtrdPIDqaV1::MakeTRDchargeDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species){
 266   //
 267   // Make Histogram for TRD Likelihood distribution
 268   //
 269   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDcharge"));
 270   if(!histo){
 271     AliError("QA histogram monitoring TRD total charge not available");
 272     return NULL;
 273   }
 274   if(species > -1 && species < AliPID::kSPECIES){
 275     // cut on species (if available)
 276     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
 277   }
 278   histo->GetAxis(3)->SetRangeUser(species + 1, species + 1);
 279
 280   TH2 *hSpec = histo->Projection(2, 1);
 281   // construct title and name
 282   TString stepname = step ? "before" : "after";
 283   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
 284   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
 285   TString histname = Form("hChargeTRD%s%s", specID.Data(), stepname.Data());
 286   TString histtitle = Form("TRD total charge for %s %s PID", speciesname.Data(), stepname.Data());
 287   hSpec->SetName(histname.Data());
 288   hSpec->SetTitle(histtitle.Data());
 289
 290   // Unset range on the original histogram
 291   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
 292   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
 293   return hSpec;
 294 }
 295