PWGHF/hfe/AliHFEtofPIDqa.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 AliHFEtofPIDqa
  17 //
  18 // Monitoring TOF PID in the HFE PID montioring framework. The following
  19 // quantities are monitored:
  20 //   TOF time distance to electron hypothesis (Number of sigmas)
  21 //   TOF time distance to the pion hypothesis (Absolute value)
  22 //   TPC dE/dx (Number of sigmas, as control histogram)
  23 // (Always as function of momentum, particle species and centrality
  24 // before and after cut)
  25 // More information about the PID monitoring framework can be found in
  26 // AliHFEpidQAmanager.cxx and AliHFEdetPIDqa.cxx
  27 //
  28 // Author:
  29 //
  30 //   Markus Fasel <M.Fasel@gsi.de>
  31 #include <TClass.h>
  32 #include <TArrayD.h>
  33 #include <TH2.h>
  34 #include <THnSparse.h>
  35 #include <TString.h>
  36
  37 #include "AliLog.h"
  38 #include "AliPID.h"
  39 #include "AliPIDResponse.h"
  40 #include "AliVParticle.h"
  41
  42 #include "AliHFEcollection.h"
  43 #include "AliHFEpid.h"
  44 #include "AliHFEpidBase.h"
  45 #include "AliHFEpidQAmanager.h"
  46 #include "AliHFEpidTPC.h"
  47 #include "AliHFEpidTRD.h"
  48 #include "AliHFEpidTOF.h"
  49 #include "AliHFEtools.h"
  50 #include "AliHFEtofPIDqa.h"
  51
  52 ClassImp(AliHFEpidTOF)
  53
  54 //_________________________________________________________
  55 AliHFEtofPIDqa::AliHFEtofPIDqa():
  56     AliHFEdetPIDqa()
  57   , fHistos(NULL)
  58 {
  59   //
  60   // Dummy constructor
  61   //
  62 }
  63
  64 //_________________________________________________________
  65 AliHFEtofPIDqa::AliHFEtofPIDqa(const char* name):
  66     AliHFEdetPIDqa(name, "QA for TOF")
  67   , fHistos(NULL)
  68 {
  69   //
  70   // Default constructor
  71   //
  72 }
  73
  74 //_________________________________________________________
  75 AliHFEtofPIDqa::AliHFEtofPIDqa(const AliHFEtofPIDqa &o):
  76     AliHFEdetPIDqa(o)
  77   , fHistos()
  78 {
  79   //
  80   // Copy constructor
  81   //
  82   o.Copy(*this);
  83 }
  84
  85 //_________________________________________________________
  86 AliHFEtofPIDqa &AliHFEtofPIDqa::operator=(const AliHFEtofPIDqa &o){
  87   //
  88   // Do assignment
  89   //
  90   AliHFEdetPIDqa::operator=(o);
  91   if(&o != this) o.Copy(*this);
  92   return *this;
  93 }
  94
  95 //_________________________________________________________
  96 AliHFEtofPIDqa::~AliHFEtofPIDqa(){
  97   //
  98   // Destructor
  99   //
 100   if(fHistos) delete fHistos;
 101 }
 102
 103 //_________________________________________________________
 104 void AliHFEtofPIDqa::Copy(TObject &o) const {
 105   //
 106   // Make copy
 107   //
 108   AliHFEtofPIDqa &target = dynamic_cast<AliHFEtofPIDqa &>(o);
 109   if(target.fHistos){
 110     delete target.fHistos;
 111     target.fHistos = NULL;
 112   }
 113   if(fHistos) target.fHistos = new AliHFEcollection(*fHistos);
 114 }
 115
 116 //_________________________________________________________
 117 Long64_t AliHFEtofPIDqa::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   AliHFEtofPIDqa *refQA = NULL;
 126   TObject *o = NULL;
 127   Long64_t count = 0;
 128   TList listHistos;
 129   while((o = it())){
 130     refQA = dynamic_cast<AliHFEtofPIDqa *>(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 AliHFEtofPIDqa::Initialize(){
 142   //
 143   // Define Histograms
 144   //
 145
 146   fHistos = new AliHFEcollection("tofqahistos", "Collection of TOF QA histograms");
 147
 148   // Make common binning
 149   const Int_t kPIDbins = AliPID::kSPECIES + 1;
 150   const Int_t kSteps = 2;
 151   const Int_t kCentralityBins = 11;
 152   const Double_t kMinPID = -1;
 153   const Double_t kMinP = 0.;
 154   const Double_t kMaxPID = (Double_t)AliPID::kSPECIES;
 155   const Double_t kMaxP = 20.;
 156   // Quantities where one can switch between low and high resolution
 157   Int_t kPbins = fQAmanager->HasHighResolutionHistos() ?  1000 : 100;
 158   Int_t kSigmaBins = fQAmanager->HasHighResolutionHistos() ? 1400 : 140;
 159
 160   // 1st histogram: TOF sigmas: (species, p nsigma, step)
 161   Int_t nBinsSigma[5] = {kPIDbins, kPbins, kSigmaBins, kSteps, kCentralityBins};
 162   Double_t minSigma[5] = {kMinPID, kMinP, -12., 0, 0};
 163   Double_t maxSigma[5] = {kMaxPID, kMaxP, 12., 2., 11.};
 164   fHistos->CreateTHnSparse("tofnSigma", "TOF signal; species; p [GeV/c]; TOF signal [a.u.]; Selection Step; Centrality", 5, nBinsSigma, minSigma, maxSigma);
 165
 166   // 3rd histogram: TPC sigmas to the electron line: (species, p nsigma, step - only filled if apriori PID information available)
 167   fHistos->CreateTHnSparse("tofMonitorTPC", "TPC signal; species; p [GeV/c]; TPC signal [a.u.]; Selection Step; Centrality", 5, nBinsSigma, minSigma, maxSigma);
 168
 169   // 4th histogram: TOF mismatch template
 170   AliHFEpidTOF *tofpid = dynamic_cast<AliHFEpidTOF *>(fQAmanager->GetDetectorPID(AliHFEpid::kTOFpid));
 171   if(!tofpid) AliDebug(1, "TOF PID not available");
 172   if(tofpid && tofpid->IsGenerateTOFmismatch()){
 173     AliDebug(1, "Prepare histogram for TOF mismatch tracks");
 174     fHistos->CreateTHnSparse("tofMismatch", "TOF signal; species; p [GeV/c]; TOF signal [a.u.]; Selection Step; Centrality", 5, nBinsSigma, minSigma, maxSigma);
 175   }
 176 }
 177
 178 //_________________________________________________________
 179 void AliHFEtofPIDqa::ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step){
 180   //
 181   // Fill TPC histograms
 182   //
 183   //AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
 184   //AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
 185
 186   Float_t centrality = track->GetCentrality();
 187   Int_t species = track->GetAbInitioPID();
 188   if(species >= AliPID::kSPECIES) species = -1;
 189
 190   AliDebug(1, Form("Monitoring particle of type %d for step %d", species, step));
 191   AliHFEpidTOF *tofpid= dynamic_cast<AliHFEpidTOF *>(fQAmanager->GetDetectorPID(AliHFEpid::kTOFpid));
 192
 193   const AliPIDResponse *pidResponse = tofpid ? tofpid->GetPIDResponse() : NULL;
 194   Double_t contentSignal[5];
 195   contentSignal[0] = species;
 196   contentSignal[1] = track->GetRecTrack()->P();
 197   contentSignal[2] = pidResponse ? pidResponse->NumberOfSigmasTOF(track->GetRecTrack(), AliPID::kElectron): 0.;
 198   contentSignal[3] = step;
 199   contentSignal[4] = centrality;
 200   fHistos->Fill("tofnSigma", contentSignal);
 201   contentSignal[2] = pidResponse ? pidResponse->NumberOfSigmasTPC(track->GetRecTrack(), AliPID::kElectron) : 0.;
 202   fHistos->Fill("tofMonitorTPC", contentSignal);
 203
 204   if(tofpid && tofpid->IsGenerateTOFmismatch()){
 205     AliDebug(1,"Filling TOF mismatch histos");
 206     TArrayD nsigmismatch(tofpid->GetNmismatchTracks());
 207     tofpid->GenerateTOFmismatch(track->GetRecTrack(),tofpid->GetNmismatchTracks(), nsigmismatch);
 208     for(int itrk = 0; itrk < tofpid->GetNmismatchTracks(); itrk++){
 209       contentSignal[2] = nsigmismatch[itrk];
 210       fHistos->Fill("tofMismatch",contentSignal);
 211     }
 212   }
 213 }
 214
 215 //_________________________________________________________
 216 TH2 *AliHFEtofPIDqa::MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
 217   //
 218   // Get the TPC control histogram for the TRD selection step (either before or after PID)
 219   //
 220   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("tofMonitorTPC"));
 221   if(!histo){
 222     AliError("QA histogram monitoring TPC nSigma not available");
 223     return NULL;
 224   }
 225   if(species > -1 && species < AliPID::kSPECIES){
 226     // cut on species (if available)
 227     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
 228   }
 229   TString centname, centtitle;
 230   Bool_t hasCentralityInfo = kTRUE;
 231   if(centralityClass > -1){
 232     if(histo->GetNdimensions() < 5){
 233       AliError("Centrality Information not available");
 234       centname = centtitle = "MinBias";
 235       hasCentralityInfo = kFALSE;
 236     } else {
 237       // Project centrality classes
 238       // -1 is Min. Bias
 239       histo->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
 240       centname = Form("Cent%d", centralityClass);
 241       centtitle = Form("Centrality %d", centralityClass);
 242     }
 243   } else {
 244     histo->GetAxis(4)->SetRange(1, histo->GetAxis(4)->GetNbins()-1);
 245     centname = centtitle = "MinBias";
 246     hasCentralityInfo = kTRUE;
 247   }
 248   histo->GetAxis(3)->SetRange(step + 1, step + 1);
 249
 250   TH2 *hSpec = histo->Projection(2, 1);
 251   // construct title and name
 252   TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
 253   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
 254   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
 255   TString histname = Form("hSigmaTPC%s%s%s", specID.Data(), stepname.Data(), centname.Data());
 256   TString histtitle = Form("TPC Sigma for %s %s PID %s", speciesname.Data(), stepname.Data(), centtitle.Data());
 257   hSpec->SetName(histname.Data());
 258   hSpec->SetTitle(histtitle.Data());
 259
 260   // Unset range on the original histogram
 261   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
 262   histo->GetAxis(3)->SetRange(0, histo->GetAxis(3)->GetNbins());
 263   if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());
 264   return hSpec;
 265 }
 266
 267 //_________________________________________________________
 268 TH2 *AliHFEtofPIDqa::MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t istep, Int_t species, Int_t centralityClass){
 269   //
 270   // Plot the Spectrum
 271   //
 272   THnSparseF *hSignal = dynamic_cast<THnSparseF *>(fHistos->Get("tofnSigma"));
 273   if(!hSignal) return NULL;
 274   hSignal->GetAxis(3)->SetRange(istep + 1, istep + 1);
 275   if(species >= 0 && species < AliPID::kSPECIES)
 276     hSignal->GetAxis(0)->SetRange(2 + species, 2 + species);
 277   TString centname, centtitle;
 278   Bool_t hasCentralityInfo = kTRUE;
 279   if(centralityClass > -1){
 280     if(hSignal->GetNdimensions() < 5){
 281       AliError("Centrality Information not available");
 282       centname = centtitle = "MinBias";
 283       hasCentralityInfo = kFALSE;
 284     } else {
 285       // Project centrality classes
 286       // -1 is Min. Bias
 287       AliInfo(Form("Centrality Class: %d", centralityClass));
 288       hSignal->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
 289       centname = Form("Cent%d", centralityClass);
 290       centtitle = Form("Centrality %d", centralityClass);
 291     }
 292   } else {
 293     hSignal->GetAxis(4)->SetRange(1, hSignal->GetAxis(4)->GetNbins()-1);
 294     centname = centtitle = "MinBias";
 295     hasCentralityInfo = kTRUE;
 296   }
 297   TH2 *hTmp = hSignal->Projection(2,1);
 298   TString hname = Form("hTOFsigma%s%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after", centname.Data()),
 299           htitle = Form("Time-of-flight spectrum[#sigma] %s selection %s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after", centtitle.Data());
 300   if(species > -1){
 301     hname += AliPID::ParticleName(species);
 302     htitle += Form(" for %ss", AliPID::ParticleName(species));
 303   }
 304   hTmp->SetName(hname.Data());
 305   hTmp->SetTitle(htitle.Data());
 306   hTmp->SetStats(kFALSE);
 307   hTmp->GetXaxis()->SetTitle("p [GeV/c]");
 308   hTmp->GetYaxis()->SetTitle("TOF time|_{el} - expected time|_{el} [#sigma]");
 309   hSignal->GetAxis(3)->SetRange(0, hSignal->GetAxis(3)->GetNbins());
 310   hSignal->GetAxis(0)->SetRange(0, hSignal->GetAxis(0)->GetNbins());
 311   if(hasCentralityInfo)hSignal->GetAxis(4)->SetRange(0, hSignal->GetAxis(4)->GetNbins());
 312   return hTmp;
 313 }
 314
 315 //_________________________________________________________
 316 TH1 *AliHFEtofPIDqa::GetHistogram(const char *name){
 317   //
 318   // Get the histogram
 319   //
 320   if(!fHistos) return NULL;
 321   TObject *histo = fHistos->Get(name);
 322   if(!histo->InheritsFrom("TH1")) return NULL;
 323   return dynamic_cast<TH1 *>(histo);
 324 }
 325