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