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 // QA class for TRD PID
17 // Plot Pion Efficiency at x electron efficiency
18 // Calculate the threshold parametrisation and save
19 // them in a root file
22 // Markus Fasel <M.Fasel@gsi.de>
31 #include <TGraphErrors.h>
32 #include <THnSparse.h>
35 #include <TIterator.h>
39 #include <TObjArray.h>
42 #include "AliAODTrack.h"
43 #include "AliAODPid.h"
44 #include "AliESDtrack.h"
45 #include "AliHFEtrdPIDqa.h"
46 #include "AliHFEtools.h"
47 #include "AliHFEpidTRD.h"
50 ClassImp(AliHFEtrdPIDqa)
52 const Double_t AliHFEtrdPIDqa::fgkElectronEff[kNElectronEffs] = {0.7,0.75, 0.8, 0.85, 0.9, 0.95};
53 //_______________________________________________________________
54 // Definition of the common binning
55 const Int_t AliHFEtrdPIDqa::fgkNBinsCommon[kQuantitiesCommon] = {
56 AliPID::kSPECIES + 1, // species
58 AliESDtrack::kTRDnPlanes + 1 // tracklets including 0
60 const Double_t AliHFEtrdPIDqa::fgkMinBinCommon[kQuantitiesCommon] = {
63 0 // tracklets including 0
66 const Double_t AliHFEtrdPIDqa::fgkMaxBinCommon[kQuantitiesCommon] = {
67 AliPID::kSPECIES, // species
69 AliESDtrack::kTRDnPlanes + 1 // tracklets including 0
71 //_______________________________________________________________
73 //__________________________________________________________________
74 AliHFEtrdPIDqa::AliHFEtrdPIDqa():
75 TNamed("trdPIDqa", ""),
78 fPionEfficiencies(NULL),
79 fProtonEfficiencies(NULL),
80 fKaonEfficiencies(NULL),
83 fTotalChargeInSlice0(kFALSE),
87 // Default Constructor
91 //__________________________________________________________________
92 AliHFEtrdPIDqa::AliHFEtrdPIDqa(const Char_t *name):
96 fPionEfficiencies(NULL),
97 fProtonEfficiencies(NULL),
98 fKaonEfficiencies(NULL),
100 fShowMessage(kFALSE),
101 fTotalChargeInSlice0(kFALSE),
109 //__________________________________________________________________
110 AliHFEtrdPIDqa::AliHFEtrdPIDqa(const AliHFEtrdPIDqa &ref):
114 fPionEfficiencies(NULL),
115 fProtonEfficiencies(NULL),
116 fKaonEfficiencies(NULL),
118 fShowMessage(kFALSE),
119 fTotalChargeInSlice0(ref.fTotalChargeInSlice0),
120 fCentralityBin(ref.fCentralityBin)
128 //__________________________________________________________________
129 AliHFEtrdPIDqa &AliHFEtrdPIDqa::operator=(const AliHFEtrdPIDqa &ref){
131 // Assignment operator
138 //__________________________________________________________________
139 AliHFEtrdPIDqa::~AliHFEtrdPIDqa(){
143 if(fTRDpid) delete fTRDpid;
144 if(fHistos) delete fHistos;
145 if(fPionEfficiencies) delete fPionEfficiencies;
146 if(fProtonEfficiencies) delete fProtonEfficiencies;
147 if(fKaonEfficiencies) delete fKaonEfficiencies;
150 //__________________________________________________________________
151 void AliHFEtrdPIDqa::Copy(TObject &ref) const{
153 // Copies content of this object into object ref
157 AliHFEtrdPIDqa &target = dynamic_cast<AliHFEtrdPIDqa &>(ref);
158 target.fTRDpid = fTRDpid;
159 target.fHistos = dynamic_cast<AliHFEcollection *>(fHistos->Clone());
160 target.fTotalChargeInSlice0 = fTotalChargeInSlice0;
161 target.fCentralityBin = fCentralityBin;
164 //__________________________________________________________________
165 Long64_t AliHFEtrdPIDqa::Merge(TCollection *coll){
170 if(coll->IsEmpty()) return 1;
172 AliHFEtrdPIDqa *refQA = NULL;
178 refQA = dynamic_cast<AliHFEtrdPIDqa *>(o);
181 listHistos.Add(refQA->fHistos);
184 fHistos->Merge(&listHistos);
188 //__________________________________________________________________
189 void AliHFEtrdPIDqa::Browse(TBrowser *b){
191 // Enable Browser functionality
194 // Add objects to the browser
195 if(fHistos) b->Add(fHistos, fHistos->GetName());
196 if(fPionEfficiencies) b->Add(fPionEfficiencies, "Pion Efficiencies");
197 if(fProtonEfficiencies) b->Add(fProtonEfficiencies, "Proton Efficiencies");
198 if(fKaonEfficiencies) b->Add(fKaonEfficiencies, "Kaon Efficiencies");
199 if(fThresholds) b->Add(fThresholds, "Thresholds");
203 //__________________________________________________________________
204 void AliHFEtrdPIDqa::Init(){
209 fHistos = new AliHFEcollection("TRDqa", "Histos for TRD QA");
211 CreateLikelihoodHistogram();
213 CreatedEdxHistogram();
214 CreateHistoTruncatedMean();
216 fTRDpid = new AliHFEpidTRD("QAtrdPID");
217 if(fTotalChargeInSlice0) fTRDpid->SetTotalChargeInSlice0();
220 //__________________________________________________________________
221 void AliHFEtrdPIDqa::CreateLikelihoodHistogram(){
223 // Create Histogram for TRD Likelihood Studies
225 Int_t nbins[kQuantitiesLike]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
226 nbins[kElectronLike] = 100;
227 nbins[kNClustersLike] = 200;
228 nbins[kCentralityBinLike] = 12;
229 Double_t binMin[kQuantitiesLike]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
230 Double_t binMax[kQuantitiesLike]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
231 binMin[kElectronLike] = 0.;
232 binMin[kNClustersLike] = 0.;
233 binMin[kCentralityBinLike] = -1.;
234 binMax[kElectronLike] = 1.;
235 binMax[kNClustersLike] = 200.;
236 binMax[kCentralityBinLike] = 11.;
238 fHistos->CreateTHnSparse("fLikeTRD","TRD Likelihood Studies", kQuantitiesLike, nbins, binMin, binMax);
239 fHistos->BinLogAxis("fLikeTRD", kP);
242 //__________________________________________________________________
243 void AliHFEtrdPIDqa::CreateQAHistogram(){
245 // Create Histogram for Basic TRD PID QA
247 AliDebug(1, "Called");
248 Int_t nbins[kQuantitiesQA]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
249 nbins[kNonZeroTrackletCharge] = AliESDtrack::kTRDnPlanes + 1;
250 nbins[kNClusters] = 200;
251 nbins[kCentralityBinQA] = 12;
252 Double_t binMin[kQuantitiesQA]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
253 binMin[kNonZeroTrackletCharge] = 0.;
254 binMin[kNClusters] = 0.;
255 binMin[kCentralityBinQA] = -1.;
256 Double_t binMax[kQuantitiesQA]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
257 binMax[kNonZeroTrackletCharge] = AliESDtrack::kTRDnPlanes + 1.;
258 binMax[kNClusters] = 200.;
259 binMax[kCentralityBinQA] = 11.;
260 fHistos->CreateTHnSparse("fQAtrack","TRD QA Histogram", kQuantitiesQA, nbins, binMin, binMax);
261 fHistos->BinLogAxis("fQAtrack", kP);
264 //__________________________________________________________________
265 void AliHFEtrdPIDqa::CreatedEdxHistogram(){
267 // Create QA histogram for dEdx investigations
269 AliDebug(1, "Called");
270 Int_t nbins[kQuantitiesdEdx]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
272 nbins[kNclusters] = 261;
273 nbins[kNonZeroSlices] = 9;
274 nbins[kCentralityBindEdx] = 12;
275 Double_t binMin[kQuantitiesdEdx]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
277 binMin[kNclusters] = 0;
278 binMin[kNonZeroSlices] = 0.;
279 binMin[kCentralityBindEdx] = -1.;
280 Double_t binMax[kQuantitiesdEdx]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
281 binMax[kdEdx] = 10000.;
282 binMax[kNclusters] = 260.;
283 binMax[kNonZeroSlices] = 8.;
284 binMax[kCentralityBindEdx] = 11.;
286 fHistos->CreateTHnSparse("fQAdEdx","TRD summed dEdx", kQuantitiesdEdx, nbins, binMin, binMax);
287 fHistos->BinLogAxis("fQAdEdx", kP);
288 fHistos->Sumw2("fQAdEdx");
291 //__________________________________________________________________
292 void AliHFEtrdPIDqa::CreateHistoTruncatedMean(){
294 // Create Histogram for Basic TRD PID QA:
296 AliDebug(1, "Called");
297 Int_t nbins[kQuantitiesTruncMean]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
298 nbins[kTPCdEdx] = 600;
299 nbins[kTRDdEdxMethod1] = 1000;
300 nbins[kTRDdEdxMethod2] = 1000;
301 nbins[kCentralityBinTruncMean] = 12;
302 Double_t binMin[kQuantitiesTruncMean]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
303 binMin[kTPCdEdx] = 0.;
304 binMin[kTRDdEdxMethod1] = 0.;
305 binMin[kTRDdEdxMethod2] = 0.;
306 binMin[kCentralityBinTruncMean] = -1.;
307 Double_t binMax[kQuantitiesTruncMean]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
308 binMax[kTPCdEdx] = 600;
309 binMax[kTRDdEdxMethod1] = 20000.;
310 binMax[kTRDdEdxMethod2] = 20000.;
311 binMax[kCentralityBinTruncMean] = 11.;
313 fHistos->CreateTHnSparse("fTRDtruncMean","TRD TruncatedMean studies", kQuantitiesTruncMean, nbins, binMin, binMax);
314 fHistos->BinLogAxis("fTRDtruncMean", kP);
315 fHistos->CreateTH2F("fTRDslicesPions","TRD dEdx per slice for Pions", 8, 0, 8, 2000, 0, 8000);
316 fHistos->CreateTH2F("fTRDslicesElectrons","TRD dEdx per slice for Electrons", 8, 0, 8, 2000, 0, 8000);
320 //__________________________________________________________________
321 void AliHFEtrdPIDqa::ProcessTracks(const TObjArray * const tracks, Int_t species){
323 // Process Electron Tracks
325 if(species < -1 || species >= AliPID::kSPECIES) return;
327 const AliVTrack *track = NULL;
328 while((track = dynamic_cast<const AliVTrack *>(it()))){
329 if(track) ProcessTrack(track, species);
333 //__________________________________________________________________
334 void AliHFEtrdPIDqa::ProcessTrack(const AliVTrack * const track, Int_t species){
336 // Process Single Track
338 if(TString(track->IsA()->GetName()).CompareTo("AliESDtrack") == 0)
339 ProcessTrackESD(dynamic_cast<const AliESDtrack *>(track), species);
340 else if(TString(track->IsA()->GetName()).CompareTo("AliAODTrack") == 0)
341 ProcessTrackAOD(dynamic_cast<const AliAODTrack *>(track), species);
345 //__________________________________________________________________
346 void AliHFEtrdPIDqa::ProcessTrackESD(const AliESDtrack *track, Int_t species){
348 // Process single ESD track
351 if((track->GetStatus() & AliESDtrack::kTRDout) == 0) return; // require TRD track
352 FillTRDLikelihoods(track, species);
353 FillTRDQAplots(track, species);
356 //__________________________________________________________________
357 void AliHFEtrdPIDqa::ProcessTrackAOD(const AliAODTrack * const track, Int_t /*species*/){
359 // Process single AOD track
360 // AOD PID object required
363 AliAODPid *trackPID = track->GetDetPid();
364 if(!trackPID) return;
368 //__________________________________________________________________
369 void AliHFEtrdPIDqa::FillTRDLikelihoods(const AliESDtrack * const track, Int_t species){
371 // Fill TRD Likelihood Histogram
373 Double_t trdLike[AliPID::kSPECIES];
374 track->GetTRDpid(trdLike);
376 Double_t norm =trdLike[AliPID::kElectron]+trdLike[AliPID::kPion];
377 Double_t likeEle = norm == 0. ? 0. : trdLike[AliPID::kElectron]/norm;
378 const AliExternalTrackParam *outerPars = track->GetOuterParam();
380 //Int_t kQuantitiesLike;
381 Double_t quantities[kQuantitiesLike]; memset(quantities, 0, sizeof(Double_t) * kQuantitiesLike);
386 // Electron Likelihood
388 quantities[kSpecies] = species;
389 quantities[kP] = outerPars ? outerPars->P() : track->P();
390 quantities[kNTracklets] = track->GetTRDntrackletsPID();
391 quantities[kElectronLike] = likeEle;
392 quantities[kNClustersLike] = track->GetTRDncls();
393 quantities[kCentralityBinLike] = fCentralityBin;
394 fHistos->Fill("fLikeTRD", quantities);
398 //__________________________________________________________________
399 void AliHFEtrdPIDqa::FillTRDQAplots(const AliESDtrack * const track, Int_t species){
401 // Fill QA Plots containing further information
403 const AliExternalTrackParam *outerPars = track->GetOuterParam();
405 Double_t quantitiesQA[kQuantitiesQA], quantitiesdEdx[kQuantitiesdEdx], quantitiesTruncMean[kQuantitiesTruncMean];
411 // Non-zero tracklet charges
412 // Number of clusters / full track
423 quantitiesQA[kSpecies] = quantitiesdEdx[kSpecies]
424 = quantitiesTruncMean[kSpecies]
426 quantitiesQA[kP] = quantitiesTruncMean[kP]
427 = outerPars ? outerPars->P() : track->P();
428 quantitiesQA[kNTracklets] = quantitiesdEdx[kNTracklets]
429 = quantitiesTruncMean[kNTracklets]
430 = track->GetTRDntrackletsPID();
431 quantitiesQA[kNClusters] = quantitiesdEdx[kNclusters] = track->GetTRDncls();
434 Double_t dEdxSum = 0., qSlice = 0.;
435 // remove the last slice from the histogram
436 Int_t ntrackletsNonZero = 0, nSlices = track->GetNumberOfTRDslices(), nSlicesNonZero = 0;
437 TString speciesname = "pions";
438 Bool_t selectedForSlicemon = kFALSE;
441 case AliPID::kElectron: speciesname = "Electrons"; selectedForSlicemon = kTRUE; break;
442 case AliPID::kPion: speciesname = "Pions"; selectedForSlicemon = kTRUE; break;
443 default: speciesname = "undefined"; selectedForSlicemon = kFALSE; break;
445 AliDebug(1, Form("species %d, speciesname %s, momentum %f, selected %s", species, speciesname.Data(), track->P(), selectedForSlicemon ? "yes" : "no"));
446 for(Int_t iplane = 0; iplane < AliESDtrack::kTRDnPlanes; iplane++){
447 quantitiesdEdx[kP] = track->GetTRDmomentum(iplane);
449 for(Int_t islice = 0; islice < nSlices; islice++){
450 if(fTotalChargeInSlice0 && islice >= 7) break;
451 qSlice = track->GetTRDslice(iplane, fTotalChargeInSlice0 ? islice + 1 : islice); // hack by mfasel: For data with the new reconstruction, slice 0 is used to store the total charge, the total number of slices is 7 instead of 8
456 // Reweighting of the slices for the truncated mean: select all pion tracks above
457 // 1.5 GeV and monitor the dEdx as function of slice
458 if(selectedForSlicemon && track->P() > 1.5){
459 AliDebug(2, Form("Filling Histogram fTRDslices%s", speciesname.Data()));
460 fHistos->Fill(Form("fTRDslices%s", speciesname.Data()), static_cast<Double_t>(islice), qSlice);
464 quantitiesdEdx[kNonZeroSlices] = nSlicesNonZero;
465 quantitiesdEdx[kdEdx] = fTotalChargeInSlice0 ? track->GetTRDslice(iplane, 0) : dEdxSum; // hack by mfasel: In the new reconstruction, the total charge is stored in the first slice, in the old reconstruction it has to be calculated from the slice charges.
466 if(dEdxSum) ntrackletsNonZero++;
467 // Fill dEdx histogram
468 quantitiesdEdx[kCentralityBindEdx] = fCentralityBin;
469 if(dEdxSum > 1e-1) fHistos->Fill("fQAdEdx", quantitiesdEdx); // Cut out 0 entries
471 quantitiesQA[kNonZeroTrackletCharge] = ntrackletsNonZero;
472 quantitiesQA[kCentralityBinQA] = fCentralityBin;
473 fHistos->Fill("fQAtrack", quantitiesQA);
475 quantitiesTruncMean[kTPCdEdx] = track->GetTPCsignal();
476 quantitiesTruncMean[kTRDdEdxMethod1] = fTRDpid->GetTRDSignalV1(track, 0.6);
477 quantitiesTruncMean[kTRDdEdxMethod2] = fTRDpid->GetTRDSignalV2(track, 0.6);
478 quantitiesTruncMean[kCentralityBinTruncMean] = fCentralityBin;
479 fHistos->Fill("fTRDtruncMean", quantitiesTruncMean);
482 /////////////////////////////////////////////////////////
484 // Code for Post Processing
486 // //////////////////////////////////////////////////////
488 //__________________________________________________________________
489 void AliHFEtrdPIDqa::FinishAnalysis(Int_t fCentrality, Bool_t isGreaterEqual){
492 // Calculate Electron Efficiency for ntracklets = 4...6
493 // Calculate thresholds for ntracklets = 4...6
496 if(!fPionEfficiencies){
497 fPionEfficiencies = new TList;
498 fPionEfficiencies->SetName("pionEfficiencies");
499 fPionEfficiencies->SetOwner();
501 if(!fProtonEfficiencies){
502 fProtonEfficiencies = new TList;
503 fProtonEfficiencies->SetName("protonEfficiencies");
504 fProtonEfficiencies->SetOwner();
507 fThresholds = new TList;
508 fThresholds->SetName("thresholds");
509 fThresholds->SetOwner();
512 for(Int_t itr = 4; itr <= 6; itr++){
514 printf("========================================\n");
515 printf("Analysing %d trackltes centrality %i \n", itr, fCentrality);
516 printf("========================================\n");
518 AnalyseNTracklets(itr, fCentrality, isGreaterEqual);
522 //__________________________________________________________________
523 void AliHFEtrdPIDqa::StoreResults(const Char_t *filename){
525 // Store histos into a root file
527 TFile *outfile = new TFile(filename, "RECREATE");
529 fPionEfficiencies->Clone()->Write(fPionEfficiencies->GetName(), kSingleKey);
530 fProtonEfficiencies->Clone()->Write(fProtonEfficiencies->GetName(), kSingleKey);
531 fThresholds->Clone()->Write(fThresholds->GetName(), kSingleKey);
536 //__________________________________________________________________
537 void AliHFEtrdPIDqa::SaveThresholdParameters(const Char_t *name, Double_t lowerLimit, Double_t upperLimit, Int_t icentrality){
539 // Fit the threshold histograms with the given parametrisation
540 // and store the TF1 in the file
544 AliError("Threshold histograms have to be created first");
549 printf("========================================\n");
550 printf("Calculating threshold parameters\n");
551 printf("========================================\n");
554 TList *outlist = new TList;
555 outlist->SetName("thresholdTRD");
557 TGraph *threshhist = NULL;
558 TF1 *threshparam = NULL;
559 TList *lHistos = NULL, *lFormulas = NULL;
560 for(Int_t itracklet = 4; itracklet <= 6; itracklet++){
563 printf("-------------------------------\n");
564 printf("Processing %d tracklets\n", itracklet);
565 printf("-------------------------------\n");
568 Char_t *listname=Form("%dTracklets", itracklet);
569 if(icentrality!=-1) listname=Form("%dTracklets%dCentrality", itracklet, icentrality);
570 lHistos = dynamic_cast<TList *>(fThresholds->FindObject(listname));
572 AliError(Form("Threshold histograms for the case %s not found", listname));
575 lFormulas = new TList;
576 lFormulas->SetName(listname);
577 outlist->Add(lFormulas);
579 for(Int_t ieff = 0; ieff < kNElectronEffs; ieff++){
582 printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n");
583 printf("Processing Electron Efficiency %f\n", fgkElectronEff[ieff]);
584 printf(">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>\n");
587 threshhist = dynamic_cast<TGraph *>(lHistos->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
588 if(!threshhist) continue;
589 threshparam = MakeThresholds(threshhist, lowerLimit, upperLimit);
590 threshparam->SetName(Form("thresh_%d_%d", itracklet, static_cast<Int_t>(fgkElectronEff[ieff] * 100)));
591 if(icentrality!=-1) threshparam->SetName(Form("thresh_%d_%d_%d", itracklet, icentrality, static_cast<Int_t>(fgkElectronEff[ieff] * 100)));
592 lFormulas->Add(threshparam);
597 TFile *outfile = new TFile(name, "RECREATE");
599 outlist->Write(outlist->GetName(), kSingleKey);
604 //__________________________________________________________________
605 void AliHFEtrdPIDqa::AnalyseNTracklets(Int_t nTracklets, Int_t nCentrality, Bool_t isGreaterEqual){
607 // Calculate Pion Efficiency, Proton Efficiency and Kaon Efficiency at discrete
608 // elPion Efficiency, Proton Efficiency and Kaon Efficiency at discrete
609 // electron efficiencies
611 THnSparse *hLikeTRD = dynamic_cast<THnSparseF *>(fHistos->Get("fLikeTRD"));
613 AliError("Likelihood Histogram not available");
617 Bool_t isPbPb = kFALSE;
618 if(nCentrality==-1) isPbPb = kFALSE;
619 if(nCentrality!=-1) isPbPb = kTRUE;
621 Int_t binTracklets = hLikeTRD->GetAxis(kNTracklets)->FindBin(nTracklets);
622 hLikeTRD->GetAxis(kNTracklets)->SetRange(binTracklets, isGreaterEqual ? 7 : binTracklets);
625 Int_t binCentrality = hLikeTRD->GetAxis(kCentralityBinLike)->FindBin(nCentrality);
626 hLikeTRD->GetAxis(kCentralityBinLike)->SetRange(binCentrality, isGreaterEqual ? 11 : binCentrality);
629 TH2 *test = hLikeTRD->Projection(kCentralityBin, kP);
634 Int_t binElectrons = hLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kElectron);
635 AliDebug(1, Form("BinElectrons %d", binElectrons));
636 Int_t binPions = hLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kPion);
637 AliDebug(1, Form("BinPions %d", binPions));
638 Int_t binProtons = hLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kProton);
639 AliDebug(1, Form("BinProtons %d", binProtons));
640 hLikeTRD->GetAxis(kSpecies)->SetRange(binElectrons, binElectrons);
641 TH2 *likeElectron = hLikeTRD->Projection(kElectronLike, kP);
642 likeElectron->SetName("likeElectron");
643 hLikeTRD->GetAxis(kSpecies)->SetRange(binPions, binPions);
644 TH2 *likePion = hLikeTRD->Projection(kElectronLike, kP);
645 likePion->SetName("likePion");
646 hLikeTRD->GetAxis(kSpecies)->SetRange(binProtons, binProtons);
647 TH2 *likeProton = hLikeTRD->Projection(kElectronLike, kP);
648 likeProton->SetName("likeProton");
651 hLikeTRD->GetAxis(kSpecies)->SetRange(0, hLikeTRD->GetAxis(kSpecies)->GetNbins());
652 hLikeTRD->GetAxis(kNTracklets)->SetRange(0, hLikeTRD->GetAxis(kNTracklets)->GetNbins());
653 hLikeTRD->GetAxis(kNTracklets)->SetRange(0, hLikeTRD->GetAxis(kCentralityBinLike)->GetNbins());
655 // Prepare List for output
656 Char_t *listname=Form("%dTracklets", nTracklets);
657 if(isPbPb) listname=Form("%dTracklets%dCentrality", nTracklets, nCentrality);
660 TList *listPions = new TList; listPions->SetName(listname); listPions->SetOwner();
661 TList *listProtons = new TList; listProtons->SetName(listname); listProtons->SetOwner();
662 TList *listThresholds = new TList; listThresholds->SetName(listname); listThresholds->SetOwner();
663 fPionEfficiencies->Add(listPions);
664 fProtonEfficiencies->Add(listProtons);
665 fThresholds->Add(listThresholds);
667 TH1 *probsEl = NULL, *probsPi = NULL, *probsPr = NULL;
668 TGraphErrors *effPi = NULL, *effPr = NULL, *thresholds = NULL;
669 Double_t p = 0, dp = 0;
671 Double_t eff, error; // value and error
672 for(Int_t ieff = 0; ieff < kNElectronEffs; ieff++){
675 printf("-----------------------------------------\n");
676 printf("Doing Electron Efficiency %f\n", fgkElectronEff[ieff]);
677 printf("-----------------------------------------\n");
679 effPi = new TGraphErrors(likeElectron->GetXaxis()->GetNbins());
680 effPi->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
681 effPr = new TGraphErrors(likeElectron->GetXaxis()->GetNbins());
682 effPr->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
683 thresholds = new TGraphErrors(likeElectron->GetXaxis()->GetNbins());
684 thresholds->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
687 listPions->Add(effPi);
688 listProtons->Add(effPr);
689 listThresholds->Add(thresholds);
691 for(Int_t imom = 1; imom <= likeElectron->GetXaxis()->GetLast(); imom++){
692 p = likeElectron->GetXaxis()->GetBinCenter(imom);
693 dp = likeElectron->GetXaxis()->GetBinWidth(imom)/2;
695 probsEl = likeElectron->ProjectionY("el", imom, imom);
696 if(!probsEl->GetEntries()) continue;
697 probsPi = likePion->ProjectionY("pi", imom, imom);
698 if(!probsPi->GetEntries()) continue;
699 probsPr = likeProton->ProjectionY("pr", imom, imom);
700 if(!probsPr->GetEntries()) continue;
701 AliDebug(1, Form("Calculating Values for p = %f", p));
703 // Calculate non-electronEfficiency and error
704 eff = CalculateHadronEfficiency(probsPi, probsEl, fgkElectronEff[ieff], threshbin, error);
705 thresholds->SetPoint(imom - 1, p, probsEl->GetXaxis()->GetBinCenter(threshbin));
706 thresholds->SetPointError(imom - 1, dp, EstimateThresholdError(probsEl, threshbin));
707 AliDebug(1, Form("threshold %d|%f", threshbin, probsEl->GetXaxis()->GetBinCenter(threshbin)));
708 AliDebug(1, Form("Pion Efficiency %f +- %f", eff, error));
709 effPi->SetPoint(imom - 1, p, eff);
710 effPi->SetPointError(imom - 1, dp, error);
711 eff = CalculateHadronEfficiency(probsPr, probsEl, fgkElectronEff[ieff] , threshbin, error);
712 AliDebug(1, Form("Proton Efficiency %f", eff));
713 effPr->SetPoint(imom - 1, p, eff);
714 effPr->SetPointError(imom - 1, dp, error);
723 // remove temporary histograms
729 //__________________________________________________________________
730 Double_t AliHFEtrdPIDqa::CalculateHadronEfficiency(const TH1 * const hadron, const TH1 *const electron, Double_t eff, Int_t &threshbin, Double_t &error){
732 // Calculate non-electron efficiency
733 // optionally returns sums as second parameter
736 TH1D eleWorking(*((const TH1D *)electron)), hadronWorking(*((const TH1D *)hadron)); // Leave the original histograms untouched and do calculation including scale on copy
737 eleWorking.Scale(1./eleWorking.Integral());
738 hadronWorking.Scale(1./hadronWorking.Integral());
740 TArrayD sumsEl(eleWorking.GetNbinsX()), sumsHd(eleWorking.GetNbinsX());
742 // calculate threshold and estimated electron efficiency the threshold was taken
743 Double_t elEff = 0.; // estimated electron efficiency at the end
744 Int_t currentBin = 0, nbins = 0;
745 for(Int_t ibin = eleWorking.GetXaxis()->GetLast(); ibin >= eleWorking.GetXaxis()->GetFirst(); ibin--){
748 elEff += eleWorking.GetBinContent(ibin);
749 sumsEl[eleWorking.GetXaxis()->GetLast() - ibin] = elEff;
751 // we found the matching bin, break the loop
755 threshbin = currentBin;
758 for(Int_t ibin = hadronWorking.GetXaxis()->GetLast(); ibin >= threshbin; ibin--) {
759 hdEff += hadronWorking.GetBinContent(ibin);
760 sumsHd[hadronWorking.GetXaxis()->GetLast() - ibin] = hdEff;
763 // search sums of electron efficiency for double counts, eliminate in electron and hadron array
764 TArrayD newsumsEl(100), newsumsHd(100);
766 for(Int_t ien = 0; ien < nbins; ien++){
768 newsumsEl[0] = sumsEl[0];
772 Int_t index = TMath::BinarySearch(nusable, newsumsEl.GetArray(), sumsEl[ien]);
773 if(TMath::Abs(sumsEl[ien] - newsumsEl[index]) < 1e-13){
774 // element already counted, don't add to the new arrays
777 newsumsEl[nusable] = sumsEl[ien];
778 newsumsHd[nusable] = sumsHd[ien];
782 //printf("New array\n");
783 //for(Int_t ib = 0; ib < nusable; ib++){
784 // printf("Electron Efficiency %f, Pion Efficiency %f\n", newsumsEl[ib], newsumsHd[ib]);
786 //printf("Do Fit\n");
790 if(hadronWorking.GetEntries() > 0 && eleWorking.GetEntries() > 0 && nusable > 2){
791 // Do error calculation in case the bins have enough statistics
792 TGraph gevh(nusable, newsumsEl.GetArray(), newsumsHd.GetArray());
793 TF1 evh("evh","pol2", eff-.05, eff+.05);
794 gevh.Fit(&evh, "Q", "", eff-.05, eff+.05);
796 // return the error of the pion efficiency
797 if(((1.-hdEff) < 0) || ((1.- elEff) < 0)){
798 AliError(" ElEffi or HdEffi > 1. Error can not be calculated. Please increase statistics!");
800 error = TMath::Sqrt(hdEff*(1-hdEff)/hadronWorking.GetEntries()+TMath::Power(evh.Derivative(eff), 2)*elEff*(1-elEff)/eleWorking.GetEntries());
802 AliDebug(2, Form("Pion Effi at [%f] : [%f +/- %f], Threshold[%f]", elEff, hdEff, error, eleWorking.GetBinCenter(threshbin)));
803 AliDebug(2, Form("Derivative at %4.2f : %f\n", eff, evh.Derivative(eff)));
809 //__________________________________________________________________
810 Double_t AliHFEtrdPIDqa::EstimateThresholdError(const TH1 * const electron, Int_t threshbin){
812 // Estimate threshold error as sqrt(N_int)/N_ent
813 // where N_int are the counts integrated under the electron selection band
814 // and N_ent are all entries in the histogram
817 for(Int_t ibin = threshbin; ibin <= electron->GetXaxis()->GetNbins(); ibin++){
818 nInt += (Int_t) electron->GetBinContent(ibin);
820 Int_t nEnt = (Int_t) electron->Integral();
821 if(nEnt == 0) return 1.;
822 return TMath::Sqrt(static_cast<Double_t>(nInt))/static_cast<Double_t>(nEnt);
825 //__________________________________________________________________
826 Double_t AliHFEtrdPIDqa::CalculateIntegratedPionEfficiency(UInt_t nTracklets, Double_t electronEff, Double_t pmin, Double_t pmax, Int_t icentrality, Double_t *error){
828 // Calculate Pion Efficiency for a given electron efficiency in the specified momentum range
830 if(nTracklets < 4 || nTracklets > 6){
831 AliError("Pion Efficiency calculation only available for 4, 5, and 6 tracklets");
834 if(electronEff < 0.6 || electronEff > 1.){
835 AliError("Pion Efficiency calculation only available in the electron efficiency range 0.6 to 1");
838 if(pmin < 0.1 || pmin > 10 || pmax < 0.1 || pmax > 10.){
839 AliError("Pion Efficiency calculation only available in the momentum range 0.1 to 10 GeV/c");
843 AliError("pmin is expected to be >= pmax");
847 // prerequierements fullfiled
849 THnSparse *hLikeTRD = dynamic_cast<THnSparseF *>(fHistos->Get("fLikeTRD"));
851 AliError("Likelihood Histogram not available");
854 Int_t binTracklets = hLikeTRD->GetAxis(kNTracklets)->FindBin(nTracklets);
855 hLikeTRD->GetAxis(kNTracklets)->SetRange(binTracklets, binTracklets);
858 Int_t binCentrality = hLikeTRD->GetAxis(kCentralityBinLike)->FindBin(icentrality);
859 hLikeTRD->GetAxis(kCentralityBinLike)->SetRange(binCentrality, binCentrality);
862 Int_t pbinMin = hLikeTRD->GetAxis(kP)->FindBin(pmax),
863 pbinMax = hLikeTRD->GetAxis(kP)->FindBin(pmax);
864 hLikeTRD->GetAxis(kP)->SetRange(pbinMin, pbinMax);
865 Int_t binElectrons = hLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kElectron);
866 Int_t binPions = hLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kPion);
867 hLikeTRD->GetAxis(kSpecies)->SetRange(binElectrons, binElectrons);
868 TH1 *likeElectron = hLikeTRD->Projection(kElectronLike);
869 likeElectron->Scale(1./likeElectron->Integral());
870 likeElectron->SetName("likeElectron");
871 hLikeTRD->GetAxis(kSpecies)->SetRange(binPions, binPions);
872 TH1 *likePion = hLikeTRD->Projection(kElectronLike);
873 likePion->Scale(1./likePion->Integral());
874 likePion->SetName("likePion");
877 hLikeTRD->GetAxis(kSpecies)->SetRange(0, hLikeTRD->GetAxis(kSpecies)->GetNbins());
878 hLikeTRD->GetAxis(kNTracklets)->SetRange(0, hLikeTRD->GetAxis(kNTracklets)->GetNbins());
879 hLikeTRD->GetAxis(kP)->SetRange(0, hLikeTRD->GetAxis(kP)->GetNbins());
882 Int_t thresh; Double_t err;
883 Double_t effpi = CalculateHadronEfficiency(likePion, likeElectron, electronEff, thresh, err);
884 delete likePion; delete likeElectron;
885 if(error) *error = err;
889 //__________________________________________________________________
890 void AliHFEtrdPIDqa::DrawTracklet(Int_t itracklet, Double_t pmin, Double_t pmax, Int_t icentrality, Bool_t doFit){
892 // Draw efficiencies and threshold as function of p
894 if(!(fPionEfficiencies && fProtonEfficiencies && fThresholds)){
895 AliError("No graphs to draw available");
899 Char_t *listname=Form("%dTracklets", itracklet);
900 if(icentrality!=-1) listname=Form("%dTracklets%dCentrality", itracklet, icentrality);
902 TList *lpions = dynamic_cast<TList *>(fPionEfficiencies->FindObject(listname));
903 TList *lprotons = dynamic_cast<TList *>(fProtonEfficiencies->FindObject(listname));
905 TList *lthresholds = dynamic_cast<TList *>(fThresholds->FindObject(listname));
906 if(!(lpions && lprotons && lthresholds)){
907 AliDebug(1, "Relevant lists not found. Did you forget to run FinishAnalysis()?");
911 TGraphErrors *pi, *pr;
914 Char_t *canvasname=Form("tracklet%d", itracklet);
915 if(icentrality!=-1) canvasname=Form("tracklet%dcentrality%d", itracklet, icentrality);
916 TCanvas *c1 = new TCanvas(canvasname, canvasname, 1024, 768);
918 TF1 *threshfit = NULL;
919 for(Int_t ieff = 0; ieff < kNElectronEffs; ieff++){
922 gPad->SetLeftMargin(0.12);
923 gPad->SetRightMargin(0.08);
924 pi = dynamic_cast<TGraphErrors *>(lpions->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
925 pr = dynamic_cast<TGraphErrors *>(lprotons->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
926 tr = dynamic_cast<TGraph *>(lthresholds->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
927 if(!(pi && pr && tr)) continue;
929 // Define nice plot, draw
931 pi->GetXaxis()->SetTitle("p / GeV/c");
932 pi->GetYaxis()->SetTitle("Efficiency");
933 pr->GetXaxis()->SetTitle("p / GeV/c");
934 pr->GetYaxis()->SetTitle("Efficiency");
935 tr->GetXaxis()->SetTitle("p / GeV/c");
936 tr->GetYaxis()->SetTitle("Efficiency");
937 pi->GetYaxis()->SetTitleOffset(1.2);
938 pr->GetYaxis()->SetTitleOffset(1.2);
939 tr->GetYaxis()->SetTitleOffset(1.2);
940 pi->GetXaxis()->SetTitleSize(0.045);
941 pi->GetYaxis()->SetTitleSize(0.045);
942 pr->GetXaxis()->SetTitleSize(0.045);
943 pr->GetYaxis()->SetTitleSize(0.045);
944 tr->GetXaxis()->SetTitleSize(0.045);
945 tr->GetYaxis()->SetTitleSize(0.045);
947 pi->GetYaxis()->SetRangeUser(1e-3, 1.);
948 pr->GetYaxis()->SetRangeUser(1e-3, 1.);
949 tr->GetYaxis()->SetRangeUser(1e-3, 1.);
950 if(pmin >= 0 && pmax >= 0.){
951 pi->GetXaxis()->SetRangeUser(pmin, pmax);
952 pr->GetXaxis()->SetRangeUser(pmin, pmax);
953 tr->GetXaxis()->SetRangeUser(pmin, pmax);
956 pi->SetMarkerColor(kRed);
957 pi->SetMarkerStyle(20);
958 pr->SetMarkerColor(kBlue);
959 pr->SetMarkerStyle(21);
960 tr->SetMarkerColor(kBlack);
961 tr->SetMarkerStyle(22);
963 pi->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
964 pr->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
965 tr->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
967 pi->Draw("ape"); pr->Draw("pesame"); tr->Draw("pesame");
971 threshfit = MakeThresholds(tr, pmin, pmax);
972 threshfit->SetLineColor(kBlack);
973 threshfit->Draw("same");
976 // Add entries to legend
978 leg = new TLegend(0.5, 0.65, 0.89, 0.85);
979 leg->SetBorderSize(0);
980 leg->SetFillStyle(0);
981 leg->AddEntry(pi, "Pion Efficiency", "lp");
982 leg->AddEntry(pr, "Proton Efficiency", "lp");
983 leg->AddEntry(tr, "Thresholds", "lp");
991 //__________________________________________________________________
992 TF1 *AliHFEtrdPIDqa::MakeThresholds(TGraph *threshist, Double_t lowerLimit, Double_t upperLimit){
994 // Create TF1 containing the threshold parametrisation
997 TF1 *threshparam = new TF1("thresh", "1-[0]-[1]*x-[2]*TMath::Exp(-[3]*x)", 0.1, 10);
998 threshist->Fit(threshparam, "NE", "", lowerLimit, upperLimit);
1002 //__________________________________________________________________
1003 void AliHFEtrdPIDqa::ClearLists(){
1005 // Clear lists for particle efficiencies and thresholds
1007 if(fPionEfficiencies){
1008 fPionEfficiencies->Delete();
1009 delete fPionEfficiencies;
1010 fPionEfficiencies = NULL;
1012 if(fProtonEfficiencies){
1013 fProtonEfficiencies->Delete();
1014 delete fProtonEfficiencies;
1015 fProtonEfficiencies = NULL;
1018 fThresholds->Delete();
1024 //__________________________________________________________________
1025 Double_t AliHFEtrdPIDqa::EvalPionEfficiency(Int_t ntls, Int_t eEff, Double_t p, Int_t icentrality){
1027 // calculate pion efficiency
1029 // Number of tracklets
1030 // Electron Efficiency
1033 TGraphErrors *measurement = GetPionEfficiency(ntls, eEff, icentrality);
1034 if(!measurement) return -1.;
1035 return measurement->Eval(p);
1038 //__________________________________________________________________
1039 Double_t AliHFEtrdPIDqa::EvalProtonEfficiency(Int_t ntls, Int_t eEff, Double_t p, Int_t icentrality){
1041 // calculate proton efficiency
1043 // Number of tracklets
1044 // Electron Efficiency
1047 TGraphErrors *measurement = GetProtonEfficiency(ntls, eEff, icentrality);
1048 if(!measurement) return -1.;
1049 return measurement->Eval(p);
1052 //__________________________________________________________________
1053 Double_t AliHFEtrdPIDqa::EvalThreshold(Int_t ntls, Int_t eEff, Double_t p, Int_t icentrality){
1055 // Get the threshold to reach a certain electron efficiency
1057 // Number of tracklets
1058 // Electron Efficiency
1061 TGraph *measurement = GetThreshold(ntls, eEff, icentrality);
1062 if(!measurement) return -1.;
1063 return measurement->Eval(p);
1066 //__________________________________________________________________
1067 TGraphErrors *AliHFEtrdPIDqa::GetPionEfficiency(Int_t ntracklets, Int_t eleffpercent, Int_t icentrality){
1069 // Get Graph with pion efficiencies
1071 Char_t *listname=Form("%dTracklets", ntracklets);
1072 if(icentrality!=-1) listname=Form("%dTracklets%dCentrality", ntracklets, icentrality);
1073 TList *graphs = dynamic_cast<TList *>(fPionEfficiencies->FindObject(listname));
1074 if(!graphs) return NULL;
1075 return dynamic_cast<TGraphErrors *>(graphs->FindObject(Form("eff%d", eleffpercent)));
1078 //__________________________________________________________________
1079 TGraphErrors *AliHFEtrdPIDqa::GetProtonEfficiency(Int_t ntracklets, Int_t eleffpercent, Int_t icentrality){
1081 // Get Graph with proton efficiencies
1083 Char_t *listname=Form("%dTracklets", ntracklets);
1084 if(icentrality!=-1) listname=Form("%dTracklets%dCentrality", ntracklets, icentrality);
1085 TList *graphs = dynamic_cast<TList *>(fProtonEfficiencies->FindObject(listname));
1086 if(!graphs) return NULL;
1087 return dynamic_cast<TGraphErrors *>(graphs->FindObject(Form("eff%d", eleffpercent)));
1090 //__________________________________________________________________
1091 TGraph *AliHFEtrdPIDqa::GetThreshold(Int_t ntracklets, Int_t eleffpercent, Int_t icentrality){
1093 // Get Graph with threshols
1095 Char_t *listname=Form("%dTracklets", ntracklets);
1096 if(icentrality!=-1) listname=Form("%dTracklets%dCentrality", ntracklets, icentrality);
1097 TList *graphs = dynamic_cast<TList *>(fThresholds->FindObject(listname));
1098 if(!graphs) return NULL;
1099 return dynamic_cast<TGraph *>(graphs->FindObject(Form("eff%d", eleffpercent)));