PWG3/hfe/AliHFEtpcPIDqa.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 AliHFEtpcPIDqa
  17 // Monitoring TPC PID in the HFE PID montioring framework. The following
  18 // quantities are monitored:
  19 //   TPC dE/dx (Number of sigmas)
  20 //   TPC dE/dx (Absolute values)
  21 // (Always as function of momentum, particle species and centrality
  22 // before and after cut)
  23 // More information about the PID monitoring framework can be found in
  24 // AliHFEpidQAmanager.cxx and AliHFEdetPIDqa.cxx
  25 //
  26 // Author:
  27 //    Markus Fasel <M.Fasel@gsi.de>
  28 //
  29 #include <TClass.h>
  30 #include <TH2.h>
  31 #include <THnSparse.h>
  32 #include <TString.h>
  33
  34 #include "AliAODTrack.h"
  35 #include "AliESDtrack.h"
  36 #include "AliESDpid.h"
  37 #include "AliLog.h"
  38 #include "AliPID.h"
  39
  40 #include "AliHFEcollection.h"
  41 #include "AliHFEpidBase.h"
  42 #include "AliHFEpidQAmanager.h"
  43 #include "AliHFEpidTPC.h"
  44 #include "AliHFEtools.h"
  45 #include "AliHFEtpcPIDqa.h"
  46
  47 ClassImp(AliHFEtpcPIDqa)
  48
  49 //_________________________________________________________
  50 AliHFEtpcPIDqa::AliHFEtpcPIDqa():
  51     AliHFEdetPIDqa()
  52   , fHistos(NULL)
  53 {
  54   //
  55   // Dummy constructor
  56   //
  57 }
  58
  59 //_________________________________________________________
  60 AliHFEtpcPIDqa::AliHFEtpcPIDqa(const char* name):
  61     AliHFEdetPIDqa(name, "QA for TPC")
  62   , fHistos(NULL)
  63 {
  64   //
  65   // Default constructor
  66   //
  67 }
  68
  69 //_________________________________________________________
  70 AliHFEtpcPIDqa::AliHFEtpcPIDqa(const AliHFEtpcPIDqa &o):
  71     AliHFEdetPIDqa(o)
  72 {
  73   //
  74   // Copy constructor
  75   //
  76   o.Copy(*this);
  77 }
  78
  79 //_________________________________________________________
  80 AliHFEtpcPIDqa &AliHFEtpcPIDqa::operator=(const AliHFEtpcPIDqa &o){
  81   //
  82   // Do assignment
  83   //
  84   AliHFEdetPIDqa::operator=(o);
  85   if(&o != this) o.Copy(*this);
  86   return *this;
  87 }
  88
  89 //_________________________________________________________
  90 AliHFEtpcPIDqa::~AliHFEtpcPIDqa(){
  91   //
  92   // Destructor
  93   //
  94   if(fHistos) delete fHistos;
  95 }
  96
  97 //_________________________________________________________
  98 void AliHFEtpcPIDqa::Copy(TObject &o) const {
  99   //
 100   // Make copy
 101   //
 102   AliHFEtpcPIDqa &target = dynamic_cast<AliHFEtpcPIDqa &>(o);
 103   if(target.fHistos){
 104     delete target.fHistos;
 105     target.fHistos = NULL;
 106   }
 107   if(fHistos) target.fHistos = new AliHFEcollection(*fHistos);
 108 }
 109
 110 //_________________________________________________________
 111 Long64_t AliHFEtpcPIDqa::Merge(TCollection *coll){
 112   //
 113   // Merge with other objects
 114   //
 115   if(!coll) return 0;
 116   if(coll->IsEmpty()) return 1;
 117
 118   TIter it(coll);
 119   AliHFEtpcPIDqa *refQA = NULL;
 120   TObject *o = NULL;
 121   Long64_t count = 0;
 122   TList listHistos;
 123   while((o = it())){
 124     refQA = dynamic_cast<AliHFEtpcPIDqa *>(o);
 125     if(!refQA) continue;
 126
 127     listHistos.Add(refQA->fHistos);
 128     count++;
 129   }
 130   fHistos->Merge(&listHistos);
 131   return count + 1;
 132 }
 133
 134 //_________________________________________________________
 135 void AliHFEtpcPIDqa::Initialize(){
 136   //
 137   // Define Histograms
 138   //
 139
 140   fHistos = new AliHFEcollection("tpcqahistos", "Collection of TPC QA histograms");
 141
 142   // Make common binning
 143   const Int_t kNdim = 5;
 144   const Int_t kPIDbins = AliPID::kSPECIES + 1;
 145   const Int_t kPbins = 100;
 146   const Int_t kSteps = 2;
 147   const Int_t kCentralityBins = 11;
 148   const Double_t kMinPID = -1;
 149   const Double_t kMinP = 0.;
 150   const Double_t kMaxPID = (Double_t)AliPID::kSPECIES;
 151   const Double_t kMaxP = 20.;
 152
 153   // 1st histogram: TPC dEdx: (species, p, dEdx, step)
 154   const Int_t kDedxbins = 200;
 155   Int_t nBinsdEdx[kNdim] = {kPIDbins, kPbins, kDedxbins, kSteps, kCentralityBins};
 156   Double_t mindEdx[kNdim] =  {kMinPID, kMinP, 0., 0., 0.};
 157   Double_t maxdEdx[kNdim] =  {kMaxPID, kMaxP, 200, 2., 11.};
 158   fHistos->CreateTHnSparse("tpcDedx", "TPC signal; species; p [GeV/c]; TPC signal [a.u.]; Centrality; Selection Step", kNdim, nBinsdEdx, mindEdx, maxdEdx);
 159   // 2nd histogram: TPC sigmas: (species, p nsigma, step)
 160   const Int_t kSigmaBins = 240;
 161   Int_t nBinsSigma[kNdim] = {kPIDbins, kPbins, kSigmaBins, kSteps, kCentralityBins};
 162   Double_t minSigma[kNdim] = {kMinPID, kMinP, -12., 0., 0.};
 163   Double_t maxSigma[kNdim] = {kMaxPID, kMaxP, 12., 2., 100.};
 164   fHistos->CreateTHnSparse("tpcnSigma", "TPC signal; species; p [GeV/c]; TPC signal [a.u.]; Centrality; Selection Step", kNdim, nBinsSigma, minSigma, maxSigma);
 165
 166   // General TPC QA
 167 }
 168
 169 //_________________________________________________________
 170 void AliHFEtpcPIDqa::ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step){
 171   //
 172   // Fill TPC histograms
 173   //
 174   AliDebug(1, Form("QA started for TPC PID for step %d", (Int_t)step));
 175   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
 176   Float_t centrality = track->GetCentrality();
 177   Int_t species = track->GetAbInitioPID();
 178   if(species >= AliPID::kSPECIES) species = -1;
 179
 180   AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fQAmanager->GetDetectorPID(AliHFEpid::kTPCpid));
 181   Double_t contentSignal[5];
 182   contentSignal[0] = species;
 183   contentSignal[1] = tpcpid->GetP(track->GetRecTrack(), anatype);
 184   contentSignal[2] = GetTPCsignal(track->GetRecTrack(), anatype);
 185   contentSignal[3] = step;
 186   contentSignal[4] = centrality;
 187   (dynamic_cast<THnSparseF *>(fHistos->Get("tpcDedx")))->Fill(contentSignal);
 188
 189   contentSignal[2] = tpcpid->NumberOfSigmas(track->GetRecTrack(), AliPID::kElectron, anatype);
 190   (dynamic_cast<THnSparseF *>(fHistos->Get("tpcnSigma")))->Fill(contentSignal);
 191 }
 192
 193 //_________________________________________________________
 194 Double_t AliHFEtpcPIDqa::GetTPCsignal(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anatype){
 195   //
 196   // Get TPC signal for ESD and AOD track
 197   //
 198   Double_t tpcSignal = 0.;
 199   if(anatype == AliHFEpidObject::kESDanalysis){
 200     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 201     tpcSignal = esdtrack->GetTPCsignal();
 202   } else {
 203     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
 204     tpcSignal = aodtrack->GetDetPid() ? aodtrack->GetDetPid()->GetTPCsignal() : 0.;
 205   }
 206   return tpcSignal;
 207 }
 208
 209 //_________________________________________________________
 210 TH2 *AliHFEtpcPIDqa::MakeSpectrumdEdx(AliHFEdetPIDqa::EStep_t istep, Int_t species){
 211   //
 212   // Plot the Spectrum
 213   //
 214   THnSparseF *hSignal = dynamic_cast<THnSparseF *>(fHistos->Get("tpcDedx"));
 215   hSignal->GetAxis(3)->SetRange(istep + 1, istep + 1);
 216   if(species > 0 && species < AliPID::kSPECIES)
 217     hSignal->GetAxis(0)->SetRange(2 + species, 2 + species);
 218   TH2 *hTmp = hSignal->Projection(2,1);
 219   Char_t hname[256], htitle[256];
 220   sprintf(hname, "hTPCsignal%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
 221   sprintf(htitle, "TPC dE/dx Spectrum %s selection", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
 222   if(species > -1){
 223     strncat(hname, AliPID::ParticleName(species), strlen(AliPID::ParticleName(species)));
 224      Char_t speciesname[256];
 225      sprintf(speciesname, " for %ss", AliPID::ParticleName(species));
 226      strncat(htitle, speciesname, strlen(speciesname));
 227   }
 228   hTmp->SetName(hname);
 229   hTmp->SetTitle(htitle);
 230   hTmp->SetStats(kFALSE);
 231   hTmp->GetXaxis()->SetTitle("p [GeV/c]");
 232   hTmp->GetYaxis()->SetTitle("TPC signal [a.u.]");
 233   hSignal->GetAxis(3)->SetRange(0, hSignal->GetAxis(3)->GetNbins());
 234   hSignal->GetAxis(0)->SetRange(0, hSignal->GetAxis(0)->GetNbins());
 235   return hTmp;
 236 }
 237
 238 //_________________________________________________________
 239 TH2 *AliHFEtpcPIDqa::MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t istep, Int_t species){
 240   //
 241   // Plot the Spectrum
 242   //
 243   THnSparseF *hSignal = dynamic_cast<THnSparseF *>(fHistos->Get("tpcnSigma"));
 244   hSignal->GetAxis(3)->SetRange(istep + 1, istep + 1);
 245   if(species >= 0 && species < AliPID::kSPECIES)
 246     hSignal->GetAxis(0)->SetRange(2 + species, 2 + species);
 247   TH2 *hTmp = hSignal->Projection(2,1);
 248   Char_t hname[256], htitle[256];
 249   sprintf(hname, "hTPCsigma%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
 250   sprintf(htitle, "TPC dE/dx Spectrum[#sigma] %s selection", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
 251   if(species > -1){
 252     strncat(hname, AliPID::ParticleName(species), strlen(AliPID::ParticleName(species)));
 253      Char_t speciesname[256];
 254      sprintf(speciesname, " for %ss", AliPID::ParticleName(species));
 255      strncat(htitle, speciesname, strlen(speciesname));
 256   }
 257   hTmp->SetName(hname);
 258   hTmp->SetTitle(htitle);
 259   hTmp->SetStats(kFALSE);
 260   hTmp->GetXaxis()->SetTitle("p [GeV/c]");
 261   hTmp->GetYaxis()->SetTitle("TPC dE/dx - <dE/dx>|_{el} [#sigma]");
 262   hSignal->GetAxis(3)->SetRange(0, hSignal->GetAxis(3)->GetNbins());
 263   hSignal->GetAxis(0)->SetRange(0, hSignal->GetAxis(0)->GetNbins());
 264   return hTmp;
 265 }
 266