PWGHF/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 <TBrowser.h>
  30 #include <TClass.h>
  31 #include <TH2.h>
  32 #include <THnSparse.h>
  33 #include <TString.h>
  34
  35 #include "AliAODTrack.h"
  36 #include "AliESDtrack.h"
  37 #include "AliLog.h"
  38 #include "AliPID.h"
  39 #include "AliPIDResponse.h"
  40
  41 #include "AliHFEcollection.h"
  42 #include "AliHFEpidBase.h"
  43 #include "AliHFEpidQAmanager.h"
  44 #include "AliHFEpidTPC.h"
  45 #include "AliHFEtools.h"
  46 #include "AliHFEtpcPIDqa.h"
  47
  48 ClassImp(AliHFEtpcPIDqa)
  49
  50 //_________________________________________________________
  51 AliHFEtpcPIDqa::AliHFEtpcPIDqa():
  52     AliHFEdetPIDqa()
  53   , fHistos(NULL)
  54   , fBrowseCentrality(-1)
  55 {
  56   //
  57   // Dummy constructor
  58   //
  59 }
  60
  61 //_________________________________________________________
  62 AliHFEtpcPIDqa::AliHFEtpcPIDqa(const char* name):
  63     AliHFEdetPIDqa(name, "QA for TPC")
  64   , fHistos(NULL)
  65   , fBrowseCentrality(-1)
  66 {
  67   //
  68   // Default constructor
  69   //
  70 }
  71
  72 //_________________________________________________________
  73 AliHFEtpcPIDqa::AliHFEtpcPIDqa(const AliHFEtpcPIDqa &o):
  74     AliHFEdetPIDqa(o)
  75   , fHistos(NULL)
  76   , fBrowseCentrality(o.fBrowseCentrality)
  77 {
  78   //
  79   // Copy constructor
  80   //
  81   o.Copy(*this);
  82 }
  83
  84 //_________________________________________________________
  85 AliHFEtpcPIDqa &AliHFEtpcPIDqa::operator=(const AliHFEtpcPIDqa &o){
  86   //
  87   // Do assignment
  88   //
  89   AliHFEdetPIDqa::operator=(o);
  90   if(&o != this) o.Copy(*this);
  91   return *this;
  92 }
  93
  94 //_________________________________________________________
  95 AliHFEtpcPIDqa::~AliHFEtpcPIDqa(){
  96   //
  97   // Destructor
  98   //
  99   if(fHistos) delete fHistos;
 100 }
 101
 102 //_________________________________________________________
 103 void AliHFEtpcPIDqa::Copy(TObject &o) const {
 104   //
 105   // Make copy
 106   //
 107   AliHFEtpcPIDqa &target = dynamic_cast<AliHFEtpcPIDqa &>(o);
 108   if(target.fHistos){
 109     delete target.fHistos;
 110     target.fHistos = NULL;
 111   }
 112   if(fHistos) target.fHistos = new AliHFEcollection(*fHistos);
 113   target.fBrowseCentrality = fBrowseCentrality;
 114 }
 115
 116 //_________________________________________________________
 117 Long64_t AliHFEtpcPIDqa::Merge(TCollection *coll){
 118   //
 119   // Merge with other objects
 120   //
 121   if(!coll) return 0;
 122   if(coll->IsEmpty()) return 1;
 123
 124   TIter it(coll);
 125   AliHFEtpcPIDqa *refQA = NULL;
 126   TObject *o = NULL;
 127   Long64_t count = 0;
 128   TList listHistos;
 129   while((o = it())){
 130     refQA = dynamic_cast<AliHFEtpcPIDqa *>(o);
 131     if(!refQA) continue;
 132
 133     listHistos.Add(refQA->fHistos);
 134     count++;
 135   }
 136   fHistos->Merge(&listHistos);
 137   return count + 1;
 138 }
 139
 140 //_________________________________________________________
 141 void AliHFEtpcPIDqa::Browse(TBrowser *b){
 142   //
 143   // Browse the PID QA
 144   //
 145   if(b){
 146     if(fHistos){
 147       b->Add(fHistos, fHistos->GetName());
 148
 149       // Make Projections of the dE/dx Spectra and add them to a new Folder
 150       TString specnames[AliPID::kSPECIES+1] = {"All", "Electrons", "Muon", "Pions", "Kaon", "Protons"};
 151       Int_t specind[AliPID::kSPECIES+1] = {-1, AliPID::kElectron, AliPID::kMuon, AliPID::kPion, AliPID::kKaon, AliPID::kProton};
 152       TList *listdEdx = new TList;
 153       listdEdx->SetOwner();
 154       TList *listNsigma = new TList;
 155       listNsigma->SetOwner();
 156
 157       TH2 *hptr = NULL;
 158       for(Int_t ispec = 0; ispec < AliPID::kSPECIES+1; ispec++){
 159         for(Int_t istep = 0; istep < 2; istep++){
 160           hptr = MakeSpectrumdEdx(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], fBrowseCentrality);
 161           hptr->SetName(Form("hTPCdEdx%s%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After" , fBrowseCentrality == -1 ? "MinBias" : Form("Cent%d", fBrowseCentrality)));
 162           listdEdx->Add(hptr);
 163           hptr = MakeSpectrumNSigma(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], fBrowseCentrality);
 164           hptr->SetName(Form("hTPCnsigma%s%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After", fBrowseCentrality == -1 ? "MinBias" : Form("Cent%d", fBrowseCentrality)));
 165           listNsigma->Add(hptr);
 166         }
 167       }
 168
 169       b->Add(listdEdx, "Projections dE/dx");
 170       b->Add(listNsigma, "Projections NSigma");
 171     }
 172   }
 173 }
 174
 175 //_________________________________________________________
 176 void AliHFEtpcPIDqa::Initialize(){
 177   //
 178   // Define Histograms
 179   //
 180
 181   fHistos = new AliHFEcollection("tpcqahistos", "Collection of TPC QA histograms");
 182
 183   // Make common binning
 184   const Int_t kNdim = 5;
 185   const Int_t kPIDbins = AliPID::kSPECIES + 1;
 186   const Int_t kSteps = 2;
 187   const Int_t kCentralityBins = 11;
 188   const Double_t kMinPID = -1;
 189   const Double_t kMinP = 0.;
 190   const Double_t kMaxPID = (Double_t)AliPID::kSPECIES;
 191   const Double_t kMaxP = 20.;
 192   // Quantities where one can switch between low and high resolution
 193   Int_t kPbins = fQAmanager->HasHighResolutionHistos() ? 1000 : 100;
 194   Int_t kDedxbins = fQAmanager->HasHighResolutionHistos() ? 400 : 200;
 195   Int_t kSigmaBins = fQAmanager->HasHighResolutionHistos() ? 1400 : 240;
 196
 197   // 1st histogram: TPC dEdx: (species, p, dEdx, step)
 198   Int_t nBinsdEdx[kNdim] = {kPIDbins, kPbins, kDedxbins, kSteps, kCentralityBins};
 199   Double_t mindEdx[kNdim] =  {kMinPID, kMinP, 0., 0., 0.};
 200   Double_t maxdEdx[kNdim] =  {kMaxPID, kMaxP, 200, 2., 11.};
 201   fHistos->CreateTHnSparse("tpcDedx", "TPC signal; species; p [GeV/c]; TPC signal [a.u.]; Selection Step; Centrality", kNdim, nBinsdEdx, mindEdx, maxdEdx);
 202   // 2nd histogram: TPC sigmas: (species, p nsigma, step)
 203   Int_t nBinsSigma[kNdim] = {kPIDbins, kPbins, kSigmaBins, kSteps, kCentralityBins};
 204   Double_t minSigma[kNdim] = {kMinPID, kMinP, -12., 0., 0.};
 205   Double_t maxSigma[kNdim] = {kMaxPID, kMaxP, 12., 2., 11.};
 206   fHistos->CreateTHnSparse("tpcnSigma", "TPC signal; species; p [GeV/c]; TPC signal [a.u.]; Selection Step; Centrality", kNdim, nBinsSigma, minSigma, maxSigma);
 207
 208   // General TPC QA
 209 }
 210
 211 //_________________________________________________________
 212 void AliHFEtpcPIDqa::ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step){
 213   //
 214   // Fill TPC histograms
 215   //
 216   AliDebug(1, Form("QA started for TPC PID for step %d", (Int_t)step));
 217   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
 218   Float_t centrality = track->GetCentrality();
 219   Int_t species = track->GetAbInitioPID();
 220   if(species >= AliPID::kSPECIES) species = -1;
 221
 222   AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fQAmanager->GetDetectorPID(AliHFEpid::kTPCpid));
 223   const AliPIDResponse *pidResponse = tpcpid ? tpcpid->GetPIDResponse() : NULL;
 224   Double_t contentSignal[5];
 225   contentSignal[0] = species;
 226   contentSignal[1] = tpcpid ? tpcpid->GetP(track->GetRecTrack(), anatype) : 0.;
 227   contentSignal[2] = GetTPCsignal(track->GetRecTrack(), anatype);
 228   contentSignal[3] = step;
 229   contentSignal[4] = centrality;
 230   fHistos->Fill("tpcDedx", contentSignal);
 231
 232   contentSignal[2] = pidResponse ? pidResponse->NumberOfSigmasTPC(track->GetRecTrack(), AliPID::kElectron) : 0.;
 233   fHistos->Fill("tpcnSigma", contentSignal);
 234 }
 235
 236 //_________________________________________________________
 237 Double_t AliHFEtpcPIDqa::GetTPCsignal(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anatype){
 238   //
 239   // Get TPC signal for ESD and AOD track
 240   //
 241   Double_t tpcSignal = 0.;
 242   if(anatype == AliHFEpidObject::kESDanalysis){
 243     const AliESDtrack *esdtrack = static_cast<const AliESDtrack *>(track);
 244     tpcSignal = esdtrack->GetTPCsignal();
 245   } else {
 246     const AliAODTrack *aodtrack = static_cast<const AliAODTrack *>(track);
 247     tpcSignal = aodtrack->GetDetPid() ? aodtrack->GetDetPid()->GetTPCsignal() : 0.;
 248   }
 249   return tpcSignal;
 250 }
 251
 252 //_________________________________________________________
 253 TH2 *AliHFEtpcPIDqa::MakeSpectrumdEdx(AliHFEdetPIDqa::EStep_t istep, Int_t species, Int_t centralityClass){
 254   //
 255   // Plot the Spectrum
 256   //
 257   THnSparseF *hSignal = dynamic_cast<THnSparseF *>(fHistos->Get("tpcDedx"));
 258   if(!hSignal) return NULL;
 259   hSignal->GetAxis(3)->SetRange(istep + 1, istep + 1);
 260   if(species >= 0 && species < AliPID::kSPECIES)
 261     hSignal->GetAxis(0)->SetRange(2 + species, 2 + species);
 262   TString hname = Form("hTPCsignal%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after"),
 263           htitle = Form("TPC dE/dx Spectrum %s selection", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
 264   if(species > -1){
 265     hname += AliPID::ParticleName(species);
 266     htitle += Form(" for %ss", AliPID::ParticleName(species));
 267   }
 268   TString centname, centtitle;
 269   Bool_t hasCentralityInfo = kTRUE;
 270   if(centralityClass > -1){
 271     if(hSignal->GetNdimensions() < 5){
 272       AliError("Centrality Information not available");
 273       centname = centtitle = "MinBias";
 274       hasCentralityInfo= kFALSE;
 275     } else {
 276       // Project centrality classes
 277       // -1 is Min. Bias
 278       hSignal->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
 279       centname = Form("Cent%d", centralityClass);
 280       centtitle = Form(" Centrality %d", centralityClass);
 281     }
 282   } else {
 283     centname = centtitle = "MinBias";
 284     hasCentralityInfo= kFALSE;
 285   }
 286   hname += centtitle;
 287   htitle += centtitle;
 288
 289   TH2 *hTmp = hSignal->Projection(2,1);
 290   hTmp->SetName(hname.Data());
 291   hTmp->SetTitle(htitle.Data());
 292   hTmp->SetStats(kFALSE);
 293   hTmp->GetXaxis()->SetTitle("p [GeV/c]");
 294   hTmp->GetYaxis()->SetTitle("TPC signal [a.u.]");
 295   hSignal->GetAxis(3)->SetRange(0, hSignal->GetAxis(3)->GetNbins());
 296   hSignal->GetAxis(0)->SetRange(0, hSignal->GetAxis(0)->GetNbins());
 297   if(hasCentralityInfo) hSignal->GetAxis(4)->SetRange(0, hSignal->GetAxis(4)->GetNbins());
 298   return hTmp;
 299 }
 300
 301 //_________________________________________________________
 302 TH2 *AliHFEtpcPIDqa::MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t istep, Int_t species, Int_t centralityClass){
 303   //
 304   // Plot the Spectrum
 305   //
 306   THnSparseF *hSignal = dynamic_cast<THnSparseF *>(fHistos->Get("tpcnSigma"));
 307   if(!hSignal) return NULL;
 308   hSignal->GetAxis(3)->SetRange(istep + 1, istep + 1);
 309   if(species >= 0 && species < AliPID::kSPECIES)
 310     hSignal->GetAxis(0)->SetRange(2 + species, 2 + species);
 311   TString hname = Form("hTPCsigma%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after"),
 312           htitle = Form("TPC dE/dx Spectrum[#sigma] %s selection", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
 313   if(species > -1){
 314     hname += AliPID::ParticleName(species);
 315     htitle += Form(" for %ss", AliPID::ParticleName(species));
 316   }
 317   TString centname, centtitle;
 318   Bool_t hasCentralityInfo = kTRUE;
 319   if(centralityClass > -1){
 320     if(hSignal->GetNdimensions() < 5){
 321       AliError("Centrality Information not available");
 322       centname = centtitle = "MinBias";
 323       hasCentralityInfo= kFALSE;
 324     } else {
 325       // Project centrality classes
 326       // -1 is Min. Bias
 327       hSignal->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
 328       centname = Form("Cent%d", centralityClass);
 329       centtitle = Form(" Centrality %d", centralityClass);
 330     }
 331   } else {
 332     centname = centtitle = "MinBias";
 333     hasCentralityInfo= kFALSE;
 334   }
 335   hname += centtitle;
 336   htitle += centtitle;
 337
 338   TH2 *hTmp = hSignal->Projection(2,1);
 339   hTmp->SetName(hname.Data());
 340   hTmp->SetTitle(htitle.Data());
 341   hTmp->SetStats(kFALSE);
 342   hTmp->GetXaxis()->SetTitle("p [GeV/c]");
 343   hTmp->GetYaxis()->SetTitle("TPC dE/dx - <dE/dx>|_{el} [#sigma]");
 344   hSignal->GetAxis(3)->SetRange(0, hSignal->GetAxis(3)->GetNbins());
 345   hSignal->GetAxis(0)->SetRange(0, hSignal->GetAxis(0)->GetNbins());
 346   if(hasCentralityInfo) hSignal->GetAxis(4)->SetRange(0, hSignal->GetAxis(4)->GetNbins());
 347   return hTmp;
 348 }
 349