1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
28 // Markus Fasel <M.Fasel@gsi.de>
34 #include <THnSparse.h>
37 #include "AliESDtrack.h"
40 #include "AliPIDResponse.h"
42 #include "AliHFEcollection.h"
43 #include "AliHFEpidBase.h"
44 #include "AliHFEpidQAmanager.h"
45 #include "AliHFEpidTPC.h"
46 #include "AliHFEpidTRD.h"
47 #include "AliHFEtrdPIDqaV1.h"
49 ClassImp(AliHFEtrdPIDqaV1)
51 //____________________________________________________________
52 AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1():
61 //____________________________________________________________
62 AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1(const Char_t *name):
63 AliHFEdetPIDqa(name, "QA for TRD"),
67 // Default constructor
71 //____________________________________________________________
72 AliHFEtrdPIDqaV1::AliHFEtrdPIDqaV1(const AliHFEtrdPIDqaV1 &o):
81 //____________________________________________________________
82 AliHFEtrdPIDqaV1 &AliHFEtrdPIDqaV1::operator=(const AliHFEtrdPIDqaV1 &o){
86 if(this == &o) return *this;
87 AliHFEdetPIDqa::operator=(o);
92 //____________________________________________________________
93 AliHFEtrdPIDqaV1::~AliHFEtrdPIDqaV1(){
97 if(fHistos) delete fHistos;
100 //_________________________________________________________
101 Long64_t AliHFEtrdPIDqaV1::Merge(TCollection *coll){
103 // Merge with other objects
106 if(coll->IsEmpty()) return 1;
109 AliHFEtrdPIDqaV1 *refQA = NULL;
114 refQA = dynamic_cast<AliHFEtrdPIDqaV1 *>(o);
117 listHistos.Add(refQA->fHistos);
120 fHistos->Merge(&listHistos);
124 //_________________________________________________________
125 void AliHFEtrdPIDqaV1::Browse(TBrowser *b){
131 b->Add(fHistos, fHistos->GetName());
133 // Make Projections of the dE/dx Spectra and add them to a new Folder
134 TString specnames[4] = {"All", "Electrons", "Pions", "Protons"};
135 Int_t specind[4] = {-1, AliPID::kElectron, AliPID::kPion, AliPID::kProton};
136 TList *listTM = new TList;
138 TList *listLike = new TList;
139 listLike->SetOwner();
140 TList *listCharge = new TList;
141 listCharge->SetOwner();
142 TList *listTPCnsigma = new TList;
143 listTPCnsigma->SetOwner();
146 TList *lctm, *lclike, *lccharge, *lcTPCsig;
147 for(Int_t icent = -1; icent < 11; icent++){
149 lctm->SetName(icent < 0 ? "MinBias" : Form("Centrality%d", icent));
153 lclike->SetName(icent < 0 ? "MinBias" : Form("Centrality%d", icent));
155 listLike->Add(lclike);
156 lccharge = new TList;
157 lccharge->SetName(icent < 0 ? "MinBias" : Form("Centrality%d", icent));
158 lccharge->SetOwner();
159 listCharge->Add(lccharge);
160 lcTPCsig = new TList;
161 lcTPCsig->SetName(icent < 0 ? "MinBias" : Form("Centrality%d", icent));
162 lcTPCsig->SetOwner();
163 listTPCnsigma->Add(lcTPCsig);
164 for(Int_t ispec = 0; ispec < 4; ispec++){
165 for(Int_t istep = 0; istep < 2; istep++){
166 hptr = MakeTRDspectrumTM(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], icent);
168 hptr->SetName(Form("hTRDtm%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
171 hptr = MakeTRDlikelihoodDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], icent);
172 hptr->SetName(Form("hTRDlike%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
174 hptr = MakeTRDchargeDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], icent);
175 hptr->SetName(Form("hTRDcharge%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
177 hptr = MakeTPCspectrumNsigma(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], icent);
178 hptr->SetName(Form("hTPCspectrum%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
184 b->Add(listTM, "Projections Truncated Mean");
185 b->Add(listLike, "Projections Likelihood distribution");
186 b->Add(listCharge, "Projections Tracklet Charge");
187 b->Add(listTPCnsigma, "Projections TPC spectra");
192 //____________________________________________________________
193 void AliHFEtrdPIDqaV1::Initialize(){
195 // Initialize QA histos for TRD PID
198 AliDebug(1, "Initializing PID QA for TRD");
199 // Make common binning
200 const Int_t kPIDbins = AliPID::kSPECIES + 1;
201 const Int_t kSteps = 2;
202 const Double_t kMinPID = -1;
203 const Double_t kMinP = 0.;
204 const Double_t kMaxPID = (Double_t)AliPID::kSPECIES;
205 const Double_t kMaxP = 20.;
206 const Int_t kCentralityBins = 11;
208 // Define number of bins
209 Int_t kPbins = fQAmanager->HasHighResolutionHistos() ? 1000 : 100;
210 Int_t tpcSigmaBins = fQAmanager->HasHighResolutionHistos() ? 1400 : 140;
211 Int_t trdLikelihoodBins = fQAmanager->HasHighResolutionHistos() ? 200 : 100;
213 fHistos = new AliHFEcollection("trdqahistos", "Collection of TRD QA histograms");
215 // Create Control Histogram monitoring the TPC sigma between the selection steps
216 Int_t nBinsTPCSigma[5] = {kPIDbins, kPbins, tpcSigmaBins, kSteps, kCentralityBins};
217 Double_t minTPCSigma[5] = {kMinPID, kMinP, -12., 0., 0.};
218 Double_t maxTPCSigma[5] = {kMaxPID, kMaxP, 12., 2., 11.};
219 fHistos->CreateTHnSparse("hTPCsigma", "TPC sigma; species p [GeV/c]; TPC dEdx - <dE/dx>|_{el} [#sigma]; selection step", 5, nBinsTPCSigma, minTPCSigma, maxTPCSigma);
220 fHistos->Sumw2("hTPCsigma");
221 // Create Monitoring histogram for the Likelihood distribution
222 Int_t nBinsTRDlike[5] = {kPIDbins, kPbins, trdLikelihoodBins, kSteps, kCentralityBins};
223 Double_t minTRDlike[5] = {kMinPID, kMinP, 0., 0., 0.};
224 Double_t maxTRDlike[5] = {kMaxPID, kMaxP, 1., 2., 11.};
225 fHistos->CreateTHnSparse("hTRDlikelihood", "TRD Likelihood Distribution; species; p [GeV/c]; TRD electron Likelihood; selection step", 5, nBinsTRDlike, minTRDlike, maxTRDlike);
226 fHistos->Sumw2("hTRDlikelihood");
227 // Create Monitoring histogram for the TRD total charge
228 const Int_t kTRDchargeBins = 1000;
229 Int_t nBinsTRDcharge[5] = {kPIDbins, kPbins, kTRDchargeBins, kSteps, kCentralityBins};
230 Double_t minTRDcharge[5] = {kMinPID, kMinP, 0., 0., 0.};
231 Double_t maxTRDcharge[5] = {kMaxPID, kMaxP, 100000., 2., 11.};
232 fHistos->CreateTHnSparse("hTRDcharge", "Total TRD charge; species; p [GeV/c]; TRD charge [a.u.]; selection step", 5, nBinsTRDcharge, minTRDcharge, maxTRDcharge);
233 fHistos->Sumw2("hTRDcharge");
234 // Monitoring of the TRD truncated mean according to version 1
235 const Int_t kTRDtmBins = 1000;
236 Int_t nBinsTRDtm[5] = {kPIDbins, kPbins, kTRDtmBins, kSteps, kCentralityBins};
237 Double_t minTRDtm[5] = {kMinPID, kMinP, 0., 0., 0.};
238 Double_t maxTRDtm[5] = {kMaxPID, kMaxP, 20000., 2., 11.};
239 fHistos->CreateTHnSparse("hTRDtruncatedMean", "TRD truncated Mean; species; p [GeV/c]; TRD signal [a.u.]; selection step", 5, nBinsTRDtm, minTRDtm, maxTRDtm);
240 fHistos->Sumw2("hTRDtruncatedMean");
242 // Monitoring of the number of tracklets
243 fHistos->CreateTH2F("hNtrackletsBefore", "Number of tracklets before PID; species; p (GeV/c), Number of Tracklets", kPbins, kMinP, kMaxP, 7, 0., 7.);
244 fHistos->Sumw2("hNtrackletsBefore");
245 fHistos->CreateTH2F("hNtrackletsAfter", "Number of tracklets after PID; species; p (GeV/c), Number of Tracklets", kPbins, kMinP, kMaxP, 7, 0., 7.);
246 fHistos->Sumw2("hNtrackletsAfter");
249 //____________________________________________________________
250 void AliHFEtrdPIDqaV1::ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step){
252 // Process the track, fill the containers
254 AliDebug(1, Form("QA started for TRD PID for step %d", (Int_t)step));
255 Int_t species = track->GetAbInitioPID();
256 if(species >= AliPID::kSPECIES) species = -1;
257 AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
259 AliHFEpidTRD *trdpid = dynamic_cast<AliHFEpidTRD *>(fQAmanager->GetDetectorPID(AliHFEpid::kTRDpid));
260 const AliPIDResponse *pidResponse = trdpid ? trdpid->GetPIDResponse() : NULL;
262 Double_t container[5];
263 container[0] = species;
264 container[1] = trdpid ? trdpid->GetP(track->GetRecTrack(), anatype) : 0.;
265 container[2] = pidResponse ? pidResponse->NumberOfSigmasTPC(track->GetRecTrack(), AliPID::kElectron) : 0.;
267 container[4] = track->GetCentrality();
268 AliDebug(1, Form("Species %d, p %f, number of sigma TPC %f, step %f, centrality %f", (Int_t)species,(Float_t) container[1],(Float_t) container[2],(Float_t) container[3],(Float_t) container[4]));
269 fHistos->Fill("hTPCsigma", container);
270 AliDebug(1, "Filled TPC sigma\n");
274 container[2] = trdpid ? trdpid->GetElectronLikelihood(static_cast<const AliVTrack*>(track->GetRecTrack()), anatype) : 0;
275 fHistos->Fill("hTRDlikelihood", container);
276 AliDebug(1, "Filled likelihood\n");
278 if(track->IsESDanalysis()){
279 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track->GetRecTrack());
281 container[2] = trdpid ? trdpid->GetTRDSignalV1(esdtrack) : 0;
282 fHistos->Fill("hTRDtruncatedMean", container);
285 for(UInt_t ily = 0; ily < 6; ily++){
286 container[2] = trdpid ? trdpid->GetChargeLayer(track->GetRecTrack(), ily, anatype) : 0;
287 if(container[2] < 1e-3) continue; // Filter out 0 entries
288 fHistos->Fill("hTRDcharge", container);
289 AliDebug(1, "Filled TRD charge\n");
292 Int_t ntracklets = track->GetRecTrack()->GetTRDntrackletsPID();
293 if(step == AliHFEdetPIDqa::kBeforePID)
294 fHistos->Fill("hNtrackletsBefore", trdpid ? trdpid->GetP(track->GetRecTrack(), anatype) : 0., ntracklets);
296 fHistos->Fill("hNtrackletsAfter", trdpid ? trdpid->GetP(track->GetRecTrack(), anatype) : 0., ntracklets);
297 AliDebug(1, "Filled tracklet\n");
300 //_________________________________________________________
301 TH2 *AliHFEtrdPIDqaV1::MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
303 // Get the TPC control histogram for the TRD selection step (either before or after PID)
305 THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTPCsigma"));
307 AliError("QA histogram monitoring TPC nSigma not available");
310 if(species > -1 && species < AliPID::kSPECIES){
311 // cut on species (if available)
312 histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
314 TString centname, centtitle;
315 Bool_t hasCentralityInfo = kTRUE;
316 if(centralityClass > -1){
317 if(histo->GetNdimensions() < 5){
318 AliError("Centrality Information not available");
319 centname = centtitle = "MinBias";
320 hasCentralityInfo = kFALSE;
322 // Project centrality classes
324 histo->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
325 centname = Form("Cent%d", centralityClass);
326 centtitle = Form("Centrality %d", centralityClass);
329 histo->GetAxis(4)->SetRange(1, histo->GetAxis(4)->GetNbins()-1);
330 centname = centtitle = "MinBias";
331 hasCentralityInfo = kTRUE;
333 histo->GetAxis(3)->SetRange(step + 1, step + 1);
335 TH2 *hSpec = histo->Projection(2, 1);
336 // construct title and name
337 TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
338 TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
339 TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
340 TString histname = Form("hSigmaTPC%s%s%s", specID.Data(), stepname.Data(), centname.Data());
341 TString histtitle = Form("TPC Sigma for %s %s PID %s", speciesname.Data(), stepname.Data(), centtitle.Data());
342 hSpec->SetName(histname.Data());
343 hSpec->SetTitle(histtitle.Data());
345 // Unset range on the original histogram
346 histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
347 histo->GetAxis(3)->SetRange(0, histo->GetAxis(3)->GetNbins());
348 if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());
352 //_________________________________________________________
353 TH2 *AliHFEtrdPIDqaV1::MakeTRDspectrumTM(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
355 // Get the TPC control histogram for the TRD selection step (either before or after PID)
357 THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDtruncatedMean"));
359 AliError("QA histogram monitoring TPC nSigma not available");
362 if(species > -1 && species < AliPID::kSPECIES){
363 // cut on species (if available)
364 histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
366 TString centname, centtitle;
367 Bool_t hasCentralityInfo = kTRUE;
368 if(centralityClass > -1){
369 if(histo->GetNdimensions() < 5){
370 AliError("Centrality Information not available");
371 centname = centtitle = "MinBias";
372 hasCentralityInfo= kFALSE;
374 // Project centrality classes
376 histo->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
377 centname = Form("Cent%d", centralityClass);
378 centtitle = Form("Centrality %d", centralityClass);
381 histo->GetAxis(4)->SetRange(1, histo->GetAxis(4)->GetNbins()-1);
382 centname = centtitle = "MinBias";
383 hasCentralityInfo= kFALSE;
385 histo->GetAxis(3)->SetRange(step + 1, step + 1);
387 TH2 *hSpec = histo->Projection(2, 1);
388 // construct title and name
389 TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
390 TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
391 TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
392 TString histname = Form("hTMTRD%s%s%s", specID.Data(), stepname.Data(), centname.Data());
393 TString histtitle = Form("TRD Truncated Mean for %s %s PID %s", speciesname.Data(), stepname.Data(), centtitle.Data());
394 hSpec->SetName(histname.Data());
395 hSpec->SetTitle(histtitle.Data());
397 // Unset range on the original histogram
398 histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
399 histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
400 if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());
404 //_________________________________________________________
405 TH2 *AliHFEtrdPIDqaV1::MakeTRDlikelihoodDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
407 // Make Histogram for TRD Likelihood distribution
409 THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDlikelihood"));
411 AliError("QA histogram monitoring TRD Electron Likelihood not available");
414 if(species > -1 && species < AliPID::kSPECIES){
415 // cut on species (if available)
416 histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
418 TString centname, centtitle;
419 Bool_t hasCentralityInfo = kTRUE;
420 if(centralityClass > -1){
421 if(histo->GetNdimensions() < 5){
422 AliError("Centrality Information not available");
423 centname = centtitle = "MinBias";
424 hasCentralityInfo= kFALSE;
426 // Project centrality classes
428 histo->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
429 centname = Form("Cent%d", centralityClass);
430 centtitle = Form("Centrality %d", centralityClass);
433 histo->GetAxis(4)->SetRange(1, histo->GetAxis(4)->GetNbins()-1);
434 centname = centtitle = "MinBias";
435 hasCentralityInfo= kTRUE;
437 histo->GetAxis(3)->SetRange(step + 1, step + 1);
439 TH2 *hSpec = histo->Projection(2, 1);
440 // construct title and name
441 TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
442 TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
443 TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
444 TString histname = Form("hLikeElTRD%s%s%s", specID.Data(), stepname.Data(),centname.Data());
445 TString histtitle = Form("TRD electron Likelihood for %s %s PID %s", speciesname.Data(), stepname.Data(),centtitle.Data());
446 hSpec->SetName(histname.Data());
447 hSpec->SetTitle(histtitle.Data());
449 // Unset range on the original histogram
450 histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
451 histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
452 if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());
456 //_________________________________________________________
457 TH2 *AliHFEtrdPIDqaV1::MakeTRDchargeDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
459 // Make Histogram for TRD Likelihood distribution
461 THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDcharge"));
463 AliError("QA histogram monitoring TRD total charge not available");
466 if(species > -1 && species < AliPID::kSPECIES){
467 // cut on species (if available)
468 histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
470 TString centname, centtitle;
471 Bool_t hasCentralityInfo = kTRUE;
472 if(centralityClass > -1){
473 if(histo->GetNdimensions() < 5){
474 AliError("Centrality Information not available");
475 centname = centtitle = "MinBias";
476 hasCentralityInfo= kFALSE;
478 // Project centrality classes
480 histo->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
481 centname = Form("Cent%d", centralityClass);
482 centtitle = Form("Centrality %d", centralityClass);
485 histo->GetAxis(4)->SetRange(1, histo->GetAxis(4)->GetNbins()-1);
486 centname = centtitle = "MinBias";
487 hasCentralityInfo= kTRUE;
489 histo->GetAxis(3)->SetRange(step + 1, step + 1);
491 TH2 *hSpec = histo->Projection(2, 1);
492 // construct title and name
493 TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
494 TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
495 TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
496 TString histname = Form("hChargeTRD%s%s%s", specID.Data(), stepname.Data(), centname.Data());
497 TString histtitle = Form("TRD total charge for %s %s PID %s", speciesname.Data(), stepname.Data(), centtitle.Data());
498 hSpec->SetName(histname.Data());
499 hSpec->SetTitle(histtitle.Data());
501 // Unset range on the original histogram
502 histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
503 histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
504 if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());