1 //////////////////////////////////////////////////////
3 // PID performance checker of the TRD
5 // Performs checks of ESD information for TRD-PID and recalculate PID based on OCDB information
6 // Also provides performance plots on detector based on the PID information - for the moment only
7 // MC source is used but V0 is also possible. Here is a list of detector performance checks
8 // - Integrated dE/dx per chamber
9 // - <PH> as function of time bin and local radial position
10 // - number of clusters/tracklet
11 // - number of tracklets/track
13 // Author : Alex Wilk <wilka@uni-muenster.de>
14 // Alex Bercuci <A.Bercuci@gsi.de>
15 // Markus Fasel <M.Fasel@gsi.de>
17 ///////////////////////////////////////////////////////
28 #include "TProfile2D.h"
30 #include "TGraphErrors.h"
33 #include <TClonesArray.h>
34 #include <TObjArray.h>
37 #include "AliESDEvent.h"
38 #include "AliESDInputHandler.h"
39 #include "AliTrackReference.h"
41 #include "AliAnalysisTask.h"
43 #include "AliTRDtrackerV1.h"
44 #include "AliTRDtrackV1.h"
45 #include "AliTRDcluster.h"
46 #include "AliTRDReconstructor.h"
47 #include "AliCDBManager.h"
48 #include "AliTRDpidUtil.h"
50 #include "Cal/AliTRDCalPID.h"
51 #include "Cal/AliTRDCalPIDNN.h"
52 #include "AliTRDcheckPID.h"
53 #include "info/AliTRDtrackInfo.h"
54 #include "info/AliTRDpidInfo.h"
56 Char_t const * AliTRDcheckPID::fgMethod[3] = {"LQ", "NN", "ESD"};
58 ClassImp(AliTRDcheckPID)
60 //________________________________________________________________________
61 AliTRDcheckPID::AliTRDcheckPID()
62 :AliTRDrecoTask("checkPID", "TRD PID checker")
68 ,fMinNTracklets(AliTRDgeometry::kNlayer)
69 ,fMaxNTracklets(AliTRDgeometry::kNlayer)
72 // Default constructor
75 fReconstructor = new AliTRDReconstructor();
76 fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
78 // Initialize momentum axis with default values
79 Double_t defaultMomenta[AliTRDCalPID::kNMom+1];
80 for(Int_t imom = 0; imom < AliTRDCalPID::kNMom+1; imom++)
81 defaultMomenta[imom] = AliTRDCalPID::GetMomentumBinning(imom);
82 SetMomentumBinning(AliTRDCalPID::kNMom, defaultMomenta);
84 memset(fEfficiency, 0, AliPID::kSPECIES*sizeof(TObjArray*));
86 fUtil = new AliTRDpidUtil();
88 DefineOutput(1, TObjArray::Class()); // PID for reference maker
92 //________________________________________________________________________
93 AliTRDcheckPID::~AliTRDcheckPID()
95 if(fPID){fPID->Delete(); delete fPID;}
96 if(fGraph){fGraph->Delete(); delete fGraph;}
97 if(fReconstructor) delete fReconstructor;
98 if(fUtil) delete fUtil;
102 //________________________________________________________________________
103 void AliTRDcheckPID::CreateOutputObjects()
108 OpenFile(0, "RECREATE");
109 fContainer = Histos();
111 fPID = new TObjArray();
112 fPID->SetOwner(kTRUE);
115 //________________________________________________________
116 void AliTRDcheckPID::Exec(Option_t *opt)
124 AliTRDrecoTask::Exec(opt);
130 //_______________________________________________________
131 TObjArray * AliTRDcheckPID::Histos(){
134 // Create QA histograms
136 if(fContainer) return fContainer;
138 Int_t xBins = AliPID::kSPECIES*fMomentumAxis->GetNbins();
139 fContainer = new TObjArray(); fContainer->Expand(kNPlots);
141 const Float_t epsilon = 1./(2*(AliTRDpidUtil::kBins-1)); // get nice histos with bin center at 0 and 1
144 const Int_t kNmethodsPID=Int_t(sizeof(fgMethod)/sizeof(Char_t*));
145 // histos with posterior probabilities for all particle species
146 for(Int_t is=AliPID::kSPECIES; is--;){
147 fEfficiency[is] = new TObjArray(kNmethodsPID);
148 fEfficiency[is]->SetOwner();
149 fEfficiency[is]->SetName(Form("Eff_%s", AliPID::ParticleShortName(is)));
150 fContainer->AddAt(fEfficiency[is], is?kEfficiencyMu+is-1:kEfficiency);
151 for(Int_t im=kNmethodsPID; im--;){
152 if(!(h = (TH2F*)gROOT->FindObject(Form("PID_%s_%s", fgMethod[im], AliPID::ParticleShortName(is))))) {
153 h = new TH2F(Form("PID_%s_%s", fgMethod[im], AliPID::ParticleShortName(is)), "", xBins, -0.5, xBins - 0.5,
154 AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon);
156 fEfficiency[is]->AddAt(h, im);
160 // histos of the dE/dx distribution for all 5 particle species and 11 momenta
161 if(!(h = (TH2F*)gROOT->FindObject("dEdx"))){
162 h = new TH2F("dEdx", "",
163 xBins, -0.5, xBins - 0.5,
167 fContainer->AddAt(h, kdEdx);
169 // histos of the dE/dx slices for all 5 particle species and 11 momenta
170 if(!(h = (TH2F*)gROOT->FindObject("dEdxSlice"))){
171 h = new TH2F("dEdxSlice", "",
172 xBins*AliTRDpidUtil::kLQslices, -0.5, xBins*AliTRDpidUtil::kLQslices - 0.5,
175 fContainer->AddAt(h, kdEdxSlice);
177 // histos of the pulse height distribution for all 5 particle species and 11 momenta
178 TObjArray *fPH = new TObjArray(2);
179 fPH->SetOwner(); fPH->SetName("PH");
180 fContainer->AddAt(fPH, kPH);
181 if(!(h = (TProfile2D*)gROOT->FindObject("PHT"))){
182 h = new TProfile2D("PHT", "",
183 xBins, -0.5, xBins - 0.5,
184 AliTRDtrackerV1::GetNTimeBins(), -0.5, AliTRDtrackerV1::GetNTimeBins() - 0.5);
187 if(!(h = (TProfile2D*)gROOT->FindObject("PHX"))){
188 h = new TProfile2D("PHX", "",
189 xBins, -0.5, xBins - 0.5,
190 AliTRDtrackerV1::GetNTimeBins(), 0., .5*AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght());
194 // histos of the number of clusters distribution for all 5 particle species and 11 momenta
195 if(!(h = (TH2F*)gROOT->FindObject("NClus"))){
196 h = new TH2F("NClus", "",
197 xBins, -0.5, xBins - 0.5,
200 fContainer->AddAt(h, kNClus);
203 // momentum distributions - absolute and in momentum bins
204 if(!(h = (TH1F*)gROOT->FindObject("hMom"))){
205 h = new TH1F("hMom", "momentum distribution", fMomentumAxis->GetNbins(), fMomentumAxis->GetXmin(), fMomentumAxis->GetXmax());
207 fContainer->AddAt(h, kMomentum);
209 if(!(h = (TH1F*)gROOT->FindObject("hMomBin"))){
210 h = new TH1F("hMomBin", "momentum distribution in momentum bins", fMomentumAxis->GetNbins(), fMomentumAxis->GetXmin(), fMomentumAxis->GetXmax());
212 fContainer->AddAt(h, kMomentumBin);
214 // Number of tracklets per track
215 if(!(h = (TH2F*)gROOT->FindObject("nTracklets"))){
216 h = new TH2F("nTracklets", "",
217 xBins, -0.5, xBins - 0.5,
218 AliTRDgeometry::kNlayer, 0.5, AliTRDgeometry::kNlayer+.5);
220 fContainer->AddAt(h, kNTracklets);
226 //________________________________________________________________________
227 Bool_t AliTRDcheckPID::CheckTrackQuality(const AliTRDtrackV1* track) const
230 // Check if the track is ok for PID
233 Int_t ntracklets = track->GetNumberOfTracklets();
234 if(ntracklets >= fMinNTracklets && ntracklets <= fMaxNTracklets) return 1;
241 //________________________________________________________________________
242 Int_t AliTRDcheckPID::CalcPDG(AliTRDtrackV1* track)
244 // Documentation to come
246 track -> SetReconstructor(fReconstructor);
247 (const_cast<AliTRDrecoParam*>(fReconstructor->GetRecoParam()))->SetPIDNeuralNetwork();
250 if(track -> GetPID(AliPID::kElectron) > track -> GetPID(AliPID::kMuon) + track -> GetPID(AliPID::kPion) + track -> GetPID(AliPID::kKaon) + track -> GetPID(AliPID::kProton)){
253 else if(track -> GetPID(kProton) > track -> GetPID(AliPID::kPion) && track -> GetPID(AliPID::kProton) > track -> GetPID(AliPID::kKaon) && track -> GetPID(AliPID::kProton) > track -> GetPID(AliPID::kMuon)){
256 else if(track -> GetPID(AliPID::kKaon) > track -> GetPID(AliPID::kMuon) && track -> GetPID(AliPID::kKaon) > track -> GetPID(AliPID::kPion)){
259 else if(track -> GetPID(AliPID::kMuon) > track -> GetPID(AliPID::kPion)){
268 //_______________________________________________________
269 TH1 *AliTRDcheckPID::PlotLQ(const AliTRDtrackV1 *track)
272 // Plot the probabilities for electrons using 2-dim LQ
276 AliDebug(2, "No ESD info available.");
280 if(track) fkTrack = track;
282 AliDebug(2, "No Track defined.");
287 status = fkESD -> GetStatus();
288 if(!(status&AliESDtrack::kTRDin)) return NULL;
290 if(!CheckTrackQuality(fkTrack)) return NULL;
292 AliTRDtrackV1 cTrack(*fkTrack);
293 cTrack.SetReconstructor(fReconstructor);
296 Float_t momentum = 0.;
299 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
300 pdg = fkMC->GetPDG();
302 //AliWarning("No MC info available!");
303 momentum = cTrack.GetMomentum(0);
304 pdg = CalcPDG(&cTrack);
306 if(!IsInRange(momentum)) return NULL;
308 (const_cast<AliTRDrecoParam*>(fReconstructor->GetRecoParam()))->SetPIDNeuralNetwork(kFALSE);
310 if(cTrack.GetNumberOfTrackletsPID() < fMinNTracklets) return NULL;
311 Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
314 TObjArray *eff(NULL);
315 for(Int_t is=AliPID::kSPECIES; is--;){
316 if(!(eff = dynamic_cast<TObjArray *>(fContainer->At(is?kEfficiencyMu+is-1:kEfficiency)))){
317 AliWarning("No Histogram List defined.");
320 if(!(hPIDLQ = dynamic_cast<TH2F *>(eff->At(AliTRDpidUtil::kLQ)))){
321 AliWarning("No Histogram defined.");
325 hPIDLQ -> Fill(FindBin(species, momentum), cTrack.GetPID(is));
333 //_______________________________________________________
334 TH1 *AliTRDcheckPID::PlotNN(const AliTRDtrackV1 *track)
337 // Plot the probabilities for electrons using 2-dim LQ
341 AliDebug(2, "No ESD info available.");
345 if(track) fkTrack = track;
347 AliDebug(2, "No Track defined.");
352 status = fkESD -> GetStatus();
353 if(!(status&AliESDtrack::kTRDin)) return NULL;
355 if(!CheckTrackQuality(fkTrack)) return NULL;
357 AliTRDtrackV1 cTrack(*fkTrack);
358 cTrack.SetReconstructor(fReconstructor);
361 Float_t momentum = 0.;
363 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
364 pdg = fkMC->GetPDG();
366 //AliWarning("No MC info available!");
367 momentum = cTrack.GetMomentum(0);
368 pdg = CalcPDG(&cTrack);
370 if(!IsInRange(momentum)) return NULL;
372 (const_cast<AliTRDrecoParam*>(fReconstructor->GetRecoParam()))->SetPIDNeuralNetwork();
374 if(cTrack.GetNumberOfTrackletsPID() < fMinNTracklets) return NULL;
376 Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
380 TObjArray *eff(NULL);
381 for(Int_t is=AliPID::kSPECIES; is--;){
382 if(!(eff = dynamic_cast<TObjArray *>(fContainer->At(is?kEfficiencyMu+is-1:kEfficiency)))){
383 AliWarning("No Histogram List defined.");
386 if(!(hPIDNN = dynamic_cast<TH2F *>(eff->At(AliTRDpidUtil::kNN)))){
387 AliWarning("No Histogram defined.");
391 hPIDNN->Fill(FindBin(species, momentum), cTrack.GetPID(is));
398 //_______________________________________________________
399 TH1 *AliTRDcheckPID::PlotESD(const AliTRDtrackV1 *track)
402 // Plot the probabilities for electrons using 2-dim LQ
406 AliDebug(2, "No ESD info available.");
410 if(track) fkTrack = track;
412 AliDebug(2, "No Track defined.");
417 status = fkESD -> GetStatus();
418 if(!(status&AliESDtrack::kTRDin)) return NULL;
420 if(!CheckTrackQuality(fkTrack)) return NULL;
421 if(fkTrack->GetNumberOfTrackletsPID() < fMinNTracklets) return NULL;
424 Float_t momentum = 0.;
426 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
427 pdg = fkMC->GetPDG();
429 //AliWarning("No MC info available!");
430 AliTRDtrackV1 cTrack(*fkTrack);
431 cTrack.SetReconstructor(fReconstructor);
432 momentum = cTrack.GetMomentum(0);
433 pdg = CalcPDG(&cTrack);
435 if(!IsInRange(momentum)) return NULL;
437 const Double32_t *pidESD = fkESD->GetResponseIter();
438 Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
441 TObjArray *eff(NULL);
442 for(Int_t is=AliPID::kSPECIES; is--;){
443 if(!(eff = dynamic_cast<TObjArray *>(fContainer->At(is?kEfficiencyMu+is-1:kEfficiency)))){
444 AliWarning("No Histogram List defined.");
447 if(!(hPID = dynamic_cast<TH2F *>(eff->At(AliTRDpidUtil::kESD)))){
448 AliWarning("No Histogram defined.");
452 AliTRDpidInfo *pid = new AliTRDpidInfo();
454 hPID->Fill(FindBin(species, momentum), pidESD[is]);
461 //_______________________________________________________
462 TH1 *AliTRDcheckPID::PlotdEdx(const AliTRDtrackV1 *track)
465 // Plot the probabilities for electrons using 2-dim LQ
468 if(track) fkTrack = track;
470 AliDebug(2, "No Track defined.");
474 if(!CheckTrackQuality(fkTrack)) return NULL;
477 if(!(hdEdx = dynamic_cast<TH2F *>(fContainer->At(kdEdx)))){
478 AliWarning("No Histogram defined.");
482 AliTRDtrackV1 cTrack(*fkTrack);
483 cTrack.SetReconstructor(fReconstructor);
485 Float_t momentum = 0.;
487 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
488 pdg = fkMC->GetPDG();
490 //AliWarning("No MC info available!");
491 momentum = cTrack.GetMomentum(0);
492 pdg = CalcPDG(&cTrack);
494 Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
495 if(!IsInRange(momentum)) return NULL;
497 (const_cast<AliTRDrecoParam*>(fReconstructor->GetRecoParam()))->SetPIDNeuralNetwork(kTRUE);
498 // (const_cast<AliTRDrecoParam*>(fReconstructor->GetRecoParam()))->SetPIDNeuralNetwork(kFALSE);
500 Int_t iBin = FindBin(species, momentum);
501 AliTRDseedV1 *tracklet = NULL;
502 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
504 tracklet = cTrack.GetTracklet(iChamb);
505 if(!tracklet) continue;
506 tracklet -> CookdEdx(AliTRDpidUtil::kNNslices);
507 for(Int_t i = AliTRDpidUtil::kNNslices; i--;) sumdEdx += tracklet->GetdEdx()[i];
508 // tracklet -> CookdEdx(AliTRDpidUtil::kLQslices);
509 // for(Int_t i = 3; i--;) sumdEdx += tracklet->GetdEdx()[i];
510 sumdEdx /= AliTRDCalPIDNN::kMLPscale;
511 hdEdx -> Fill(iBin, sumdEdx);
517 //_______________________________________________________
518 TH1 *AliTRDcheckPID::PlotdEdxSlice(const AliTRDtrackV1 *track)
521 // Plot the probabilities for electrons using 2-dim LQ
524 if(track) fkTrack = track;
526 AliDebug(2, "No Track defined.");
530 if(!CheckTrackQuality(fkTrack)) return NULL;
533 if(!(hdEdxSlice = dynamic_cast<TH2F *>(fContainer->At(kdEdxSlice)))){
534 AliWarning("No Histogram defined.");
538 AliTRDtrackV1 cTrack(*fkTrack);
539 cTrack.SetReconstructor(fReconstructor);
541 Float_t momentum = 0.;
543 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
544 pdg = fkMC->GetPDG();
546 //AliWarning("No MC info available!");
547 momentum = cTrack.GetMomentum(0);
548 pdg = CalcPDG(&cTrack);
550 if(!IsInRange(momentum)) return NULL;
552 (const_cast<AliTRDrecoParam*>(fReconstructor->GetRecoParam()))->SetPIDNeuralNetwork(kFALSE);
553 Int_t iMomBin = fMomentumAxis->FindBin(momentum);
554 Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
556 AliTRDseedV1 *tracklet = NULL;
557 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
558 tracklet = cTrack.GetTracklet(iChamb);
559 if(!tracklet) continue;
560 tracklet -> CookdEdx(AliTRDpidUtil::kLQslices);
561 fdEdx = const_cast<Float_t *>(tracklet->GetdEdx());
562 for(Int_t iSlice = 0; iSlice < AliTRDpidUtil::kLQslices; iSlice++){
563 hdEdxSlice -> Fill(species * fMomentumAxis->GetNbins() * AliTRDpidUtil::kLQslices + (iMomBin-1) * AliTRDpidUtil::kLQslices + iSlice, fdEdx[iSlice]);
571 //_______________________________________________________
572 TH1 *AliTRDcheckPID::PlotPH(const AliTRDtrackV1 *track)
575 // Plot the probabilities for electrons using 2-dim LQ
578 if(track) fkTrack = track;
580 AliDebug(2, "No Track defined.");
584 if(!CheckTrackQuality(fkTrack)) return NULL;
586 TObjArray *arr = NULL;
587 TProfile2D *hPHX, *hPHT;
588 if(!(arr = dynamic_cast<TObjArray *>(fContainer->At(kPH)))){
589 AliWarning("No Histogram defined.");
592 hPHT = (TProfile2D*)arr->At(0);
593 hPHX = (TProfile2D*)arr->At(1);
596 Float_t momentum = 0.;
598 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
599 pdg = fkMC->GetPDG();
601 //AliWarning("No MC info available!");
602 AliTRDtrackV1 cTrack(*fkTrack);
603 cTrack.SetReconstructor(fReconstructor);
604 momentum = cTrack.GetMomentum(0);
605 pdg = CalcPDG(&cTrack);
607 if(!IsInRange(momentum)) return NULL;;
609 AliTRDseedV1 *tracklet = NULL;
610 AliTRDcluster *cluster = NULL;
611 Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
612 Int_t iBin = FindBin(species, momentum);
613 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
614 tracklet = fkTrack->GetTracklet(iChamb);
615 if(!tracklet) continue;
616 Float_t x0 = tracklet->GetX0();
617 for(Int_t ic = 0; ic < AliTRDtrackerV1::GetNTimeBins(); ic++){
618 if(!(cluster = tracklet->GetClusters(ic))) continue;
619 hPHT -> Fill(iBin, cluster->GetLocalTimeBin(), TMath::Abs(cluster->GetQ()));
620 hPHX -> Fill(iBin, x0 - cluster->GetX(), tracklet->GetdQdl(ic));
627 //_______________________________________________________
628 TH1 *AliTRDcheckPID::PlotNClus(const AliTRDtrackV1 *track)
631 // Plot the probabilities for electrons using 2-dim LQ
634 if(track) fkTrack = track;
636 AliDebug(2, "No Track defined.");
640 if(!CheckTrackQuality(fkTrack)) return NULL;
643 if(!(hNClus = dynamic_cast<TH2F *>(fContainer->At(kNClus)))){
644 AliWarning("No Histogram defined.");
650 Float_t momentum = 0.;
652 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
653 pdg = fkMC->GetPDG();
655 //AliWarning("No MC info available!");
656 AliTRDtrackV1 cTrack(*fkTrack);
657 cTrack.SetReconstructor(fReconstructor);
658 momentum = cTrack.GetMomentum(0);
659 pdg = CalcPDG(&cTrack);
661 if(!IsInRange(momentum)) return NULL;
663 Int_t species = AliTRDpidUtil::AliTRDpidUtil::Pdg2Pid(pdg);
664 Int_t iBin = FindBin(species, momentum);
665 AliTRDseedV1 *tracklet = NULL;
666 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
667 tracklet = fkTrack->GetTracklet(iChamb);
668 if(!tracklet) continue;
669 hNClus -> Fill(iBin, tracklet->GetN());
675 //_______________________________________________________
676 TH1 *AliTRDcheckPID::PlotNTracklets(const AliTRDtrackV1 *track)
679 // Plot the probabilities for electrons using 2-dim LQ
682 if(track) fkTrack = track;
684 AliDebug(2, "No Track defined.");
689 if(!(hTracklets = dynamic_cast<TH2F *>(fContainer->At(kNTracklets)))){
690 AliWarning("No Histogram defined.");
694 AliTRDtrackV1 cTrack(*fkTrack);
695 cTrack.SetReconstructor(fReconstructor);
697 Float_t momentum = 0.;
699 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
700 pdg = fkMC->GetPDG();
702 //AliWarning("No MC info available!");
703 momentum = cTrack.GetMomentum(0);
704 pdg = CalcPDG(&cTrack);
706 Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
707 if(!IsInRange(momentum)) return NULL;
709 Int_t iBin = FindBin(species, momentum);
710 hTracklets -> Fill(iBin, cTrack.GetNumberOfTracklets());
714 //_______________________________________________________
715 TH1 *AliTRDcheckPID::PlotMom(const AliTRDtrackV1 *track)
718 // Plot the probabilities for electrons using 2-dim LQ
721 if(track) fkTrack = track;
723 AliDebug(2, "No Track defined.");
727 if(!CheckTrackQuality(fkTrack)) return NULL;
730 if(!(hMom = dynamic_cast<TH1F *>(fContainer->At(kMomentum)))){
731 AliWarning("No Histogram defined.");
737 Float_t momentum = 0.;
739 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
740 pdg = fkMC->GetPDG();
742 //AliWarning("No MC info available!");
743 AliTRDtrackV1 cTrack(*fkTrack);
744 cTrack.SetReconstructor(fReconstructor);
745 momentum = cTrack.GetMomentum(0);
746 pdg = CalcPDG(&cTrack);
748 if(IsInRange(momentum)) hMom -> Fill(momentum);
753 //_______________________________________________________
754 TH1 *AliTRDcheckPID::PlotMomBin(const AliTRDtrackV1 *track)
757 // Plot the probabilities for electrons using 2-dim LQ
760 if(track) fkTrack = track;
762 AliDebug(2, "No Track defined.");
766 if(!CheckTrackQuality(fkTrack)) return NULL;
769 if(!(hMomBin = dynamic_cast<TH1F *>(fContainer->At(kMomentumBin)))){
770 AliWarning("No Histogram defined.");
776 Float_t momentum = 0.;
779 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
780 pdg = fkMC->GetPDG();
782 //AliWarning("No MC info available!");
783 AliTRDtrackV1 cTrack(*fkTrack);
784 cTrack.SetReconstructor(fReconstructor);
785 momentum = cTrack.GetMomentum(0);
787 if(IsInRange(momentum)) hMomBin -> Fill(fMomentumAxis->FindBin(momentum));
792 //________________________________________________________
793 Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig)
795 // Steering function to retrieve performance plots
797 Bool_t kFIRST = kTRUE;
798 TGraphErrors *g = NULL;
800 TObjArray *arr = NULL;
801 TH1 *h1 = NULL, *h=NULL;
803 TList *content = NULL;
806 gPad->Divide(2, 1, 1.e-5, 1.e-5);
807 TList *l=gPad->GetListOfPrimitives();
808 TVirtualPad *pad = ((TVirtualPad*)l->At(0));pad->cd();
809 pad->SetMargin(0.1, 0.01, 0.1, 0.01);
811 TLegend *legEff = new TLegend(.64, .84, .98, .98);
812 legEff->SetBorderSize(1);legEff->SetTextSize(0.03255879);
813 legEff->SetFillColor(0);
814 h=new TH1S("hEff", "", 1, .5, 11.);
815 h->SetLineColor(1);h->SetLineWidth(1);
817 ax->SetTitle("p [GeV/c]");
818 ax->SetRangeUser(.5, 11.);
819 ax->SetMoreLogLabels();
821 ax->SetTitle("#epsilon_{#pi} [%]");
823 ax->SetRangeUser(1.e-2, 10.);
825 content = (TList *)fGraph->FindObject(Form("Eff_%s", AliTRDCalPID::GetPartName(AliPID::kPion)));
826 if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kLQ))) break;
827 if(!g->GetN()) break;
828 legEff->SetHeader("PID Method [PION]");
829 g->Draw("pc"); legEff->AddEntry(g, "LQ 2D", "pl");
830 if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kNN))) break;
831 g->Draw("pc"); legEff->AddEntry(g, "NN", "pl");
832 if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kESD))) break;
833 g->Draw("p"); legEff->AddEntry(g, "ESD", "pl");
841 pad = ((TVirtualPad*)l->At(1));pad->cd();
842 pad->SetMargin(0.1, 0.01, 0.1, 0.01);
843 h=new TH1S("hThr", "", 1, .5, 11.);
844 h->SetLineColor(1);h->SetLineWidth(1);
846 ax->SetTitle("p [GeV/c]");
847 ax->SetMoreLogLabels();
849 ax->SetTitle("Threshold [%]");
850 ax->SetRangeUser(5.e-2, 1.);
852 content = (TList *)fGraph->FindObject("Thres");
853 if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kLQ))) break;
854 if(!g->GetN()) break;
856 if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kNN))) break;
858 if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kESD))) break;
866 gPad->Divide(2, 1, 1.e-5, 1.e-5);
867 TList *l=gPad->GetListOfPrimitives();
868 TVirtualPad *pad = ((TVirtualPad*)l->At(0));pad->cd();
869 pad->SetMargin(0.1, 0.01, 0.1, 0.01);
871 TLegend *legEff = new TLegend(.64, .84, .98, .98);
872 legEff->SetBorderSize(1);legEff->SetTextSize(0.03255879);
873 legEff->SetFillColor(0);
874 h = (TH1S*)gROOT->FindObject("hEff");
875 h=(TH1S*)h->Clone("hEff_K");
876 h->SetYTitle("#epsilon_{K} [%]");
877 h->GetYaxis()->SetRangeUser(1.e-2, 1.e2);
879 content = (TList *)fGraph->FindObject(Form("Eff_%s", AliTRDCalPID::GetPartName(AliPID::kKaon)));
880 if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kLQ))) break;
881 if(!g->GetN()) break;
882 legEff->SetHeader("PID Method [KAON]");
883 g->Draw("pc"); legEff->AddEntry(g, "LQ 2D", "pl");
884 if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kNN))) break;
885 g->Draw("pc"); legEff->AddEntry(g, "NN", "pl");
886 if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kESD))) break;
887 g->Draw("p"); legEff->AddEntry(g, "ESD", "pl");
894 TLegend *legEff2 = new TLegend(.64, .84, .98, .98);
895 legEff2->SetBorderSize(1);legEff2->SetTextSize(0.03255879);
896 legEff2->SetFillColor(0);
897 pad = ((TVirtualPad*)l->At(1));pad->cd();
898 pad->SetMargin(0.1, 0.01, 0.1, 0.01);
899 h=(TH1S*)h->Clone("hEff_p");
900 h->SetYTitle("#epsilon_{p} [%]");
901 h->GetYaxis()->SetRangeUser(1.e-2, 1.e2);
903 content = (TList *)fGraph->FindObject(Form("Eff_%s", AliTRDCalPID::GetPartName(AliPID::kProton)));
904 if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kLQ))) break;
905 if(!g->GetN()) break;
906 legEff2->SetHeader("PID Method [PROTON]");
907 g->Draw("pc"); legEff2->AddEntry(g, "LQ 2D", "pl");
908 if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kNN))) break;
909 g->Draw("pc"); legEff2->AddEntry(g, "NN", "pl");
910 if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kESD))) break;
911 g->Draw("p"); legEff2->AddEntry(g, "ESD", "pl");
920 // save 2.0 GeV projection as reference
921 TLegend *legdEdx = new TLegend(.7, .7, .98, .98);
922 legdEdx->SetBorderSize(1);
924 if(!(h2 = (TH2F*)(fContainer->At(kdEdx)))) break;
925 legdEdx->SetHeader("Particle Species");
926 gPad->SetMargin(0.1, 0.01, 0.1, 0.01);
927 for(Int_t is = AliPID::kSPECIES-1; is>=0; is--){
928 Int_t bin = FindBin(is, 2.);
929 h1 = h2->ProjectionY(Form("px%d", is), bin, bin);
930 if(!h1->GetEntries()) continue;
931 h1->Scale(1./h1->Integral());
932 h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
934 h1->GetXaxis()->SetTitle("dE/dx (a.u.)");
935 h1->GetYaxis()->SetTitle("<Entries>");
937 h = (TH1F*)h1->DrawClone(kFIRST ? "c" : "samec");
938 legdEdx->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "l");
952 gPad->Divide(2, 1, 1.e-5, 1.e-5);
953 TList *l=gPad->GetListOfPrimitives();
955 // save 2.0 GeV projection as reference
956 TLegend *legPH = new TLegend(.4, .7, .68, .98);
957 legPH->SetBorderSize(1);legPH->SetFillColor(0);
958 legPH->SetHeader("Particle Species");
959 if(!(arr = (TObjArray*)(fContainer->At(kPH)))) break;
960 if(!(h2 = (TProfile2D*)(arr->At(0)))) break;
962 TVirtualPad *pad = ((TVirtualPad*)l->At(0));pad->cd();
963 pad->SetMargin(0.1, 0.01, 0.1, 0.01);
965 for(Int_t is=0; is<AliPID::kSPECIES; is++){
966 Int_t bin = FindBin(is, 2.);
967 h1 = h2->ProjectionY(Form("pyt%d", is), bin, bin);
968 if(!h1->GetEntries()) continue;
969 h1->SetMarkerStyle(24);
970 h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
971 h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
973 h1->GetXaxis()->SetTitle("t_{drift} [100*ns]");
974 h1->GetYaxis()->SetTitle("<dQ/dt> [a.u.]");
976 h = (TH1F*)h1->DrawClone(kFIRST ? "c" : "samec");
977 legPH->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "pl");
981 pad = ((TVirtualPad*)l->At(1));pad->cd();
982 pad->SetMargin(0.1, 0.01, 0.1, 0.01);
983 if(!(h2 = (TProfile2D*)(arr->At(1)))) break;
985 for(Int_t is=0; is<AliPID::kSPECIES; is++){
986 Int_t bin = FindBin(is, 2.);
987 h1 = h2->ProjectionY(Form("pyx%d", is), bin, bin);
988 if(!h1->GetEntries()) continue;
989 h1->SetMarkerStyle(24);
990 h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
991 h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
993 h1->GetXaxis()->SetTitle("x_{drift} [cm]");
994 h1->GetYaxis()->SetTitle("<dQ/dl> [a.u./cm]");
996 h = (TH1F*)h1->DrawClone(kFIRST ? "c" : "samec");
1009 // save 2.0 GeV projection as reference
1010 TLegend *legNClus = new TLegend(.13, .7, .4, .98);
1011 legNClus->SetBorderSize(1);
1012 legNClus->SetFillColor(0);
1015 if(!(h2 = (TH2F*)(fContainer->At(kNClus)))) break;
1016 legNClus->SetHeader("Particle Species");
1017 for(Int_t is=0; is<AliPID::kSPECIES; is++){
1018 Int_t bin = FindBin(is, 2.);
1019 h1 = h2->ProjectionY(Form("pyNClus%d", is), bin, bin);
1020 if(!h1->GetEntries()) continue;
1021 h1->Scale(100./h1->Integral());
1022 //h1->SetMarkerStyle(24);
1023 //h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
1024 h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
1026 h1->GetXaxis()->SetTitle("N^{cl}/tracklet");
1027 h1->GetYaxis()->SetTitle("Prob. [%]");
1028 h = (TH1F*)h1->DrawClone("c");
1030 h->GetXaxis()->SetRangeUser(0., 35.);
1032 } else h = (TH1F*)h1->DrawClone("samec");
1034 legNClus->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "l");
1048 TLegend *legNClus = new TLegend(.4, .7, .68, .98);
1049 legNClus->SetBorderSize(1);
1051 if(!(h2 = (TH2F*)(fContainer->At(kNTracklets)))) break;
1052 legNClus->SetHeader("Particle Species");
1053 for(Int_t is=0; is<AliPID::kSPECIES; is++){
1054 Int_t bin = FindBin(is, 2.);
1055 h1 = h2->ProjectionY(Form("pyNTracklets%d", is), bin, bin);
1056 if(!h1->GetEntries()) continue;
1057 h1->Scale(100./h1->Integral());
1058 //h1->SetMarkerStyle(24);
1059 //h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
1060 h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
1062 h1->GetXaxis()->SetTitle("N^{trklt}/track");
1063 h1->GetXaxis()->SetRangeUser(1.,6.);
1064 h1->GetYaxis()->SetTitle("Prob. [%]");
1066 h = (TH1F*)h1->DrawClone(kFIRST ? "c" : "samec");
1067 legNClus->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "l");
1079 AliInfo(Form("Reference plot [%d] missing result", ifig));
1083 //________________________________________________________________________
1084 Bool_t AliTRDcheckPID::PostProcess()
1086 // Draw result to the screen
1087 // Called once at the end of the query
1090 Printf("ERROR: list not available");
1094 TObjArray *eff(NULL);
1096 fGraph = new TObjArray(2*AliPID::kSPECIES);
1099 if(!(eff = dynamic_cast<TObjArray *>(fContainer->At(kEfficiency)))){
1100 AliError("Efficiency container for Electrons missing.");
1103 EvaluateEfficiency(eff, fGraph, AliPID::kPion, 0.9);
1104 EvaluateEfficiency(eff, fGraph, AliPID::kKaon, 0.9);
1105 EvaluateEfficiency(eff, fGraph, AliPID::kProton, 0.9);
1111 //________________________________________________________________________
1112 void AliTRDcheckPID::EvaluateEfficiency(TObjArray * const histoContainer, TObjArray *results, Int_t species, Float_t electronEfficiency){
1113 // Process PID information for pion efficiency
1115 fUtil->SetElectronEfficiency(electronEfficiency);
1118 const Int_t kNmethodsPID=Int_t(sizeof(fgMethod)/sizeof(Char_t*));
1119 Color_t colors[kNmethodsPID] = {kBlue, kGreen+2, kRed};
1120 Int_t markerStyle[kNmethodsPID] = {7, 7, 24};
1121 // efficiency graphs
1122 TGraphErrors *g(NULL);
1123 TObjArray *eff = new TObjArray(kNmethodsPID); eff->SetOwner(); eff->SetName(Form("Eff_%s", AliTRDCalPID::GetPartName(species)));
1124 results->AddAt(eff, species);
1125 for(Int_t iMethod = 0; iMethod < kNmethodsPID; iMethod++){
1126 eff->AddAt(g = new TGraphErrors(), iMethod);
1127 g->SetName(Form("%s", fgMethod[iMethod]));
1128 g->SetLineColor(colors[iMethod]);
1129 g->SetMarkerColor(colors[iMethod]);
1130 g->SetMarkerStyle(markerStyle[iMethod]);
1133 // Threshold graphs if not already
1134 TObjArray *thres(NULL);
1135 if(!(results->At(AliPID::kSPECIES))){
1136 thres = new TObjArray(kNmethodsPID); thres->SetOwner();
1137 thres->SetName("Thres");
1138 results->AddAt(thres, AliPID::kSPECIES);
1139 for(Int_t iMethod = 0; iMethod < kNmethodsPID; iMethod++){
1140 thres->AddAt(g = new TGraphErrors(), iMethod);
1141 g->SetName(Form("%s", fgMethod[iMethod]));
1142 g->SetLineColor(colors[iMethod]);
1143 g->SetMarkerColor(colors[iMethod]);
1144 g->SetMarkerStyle(markerStyle[iMethod]);
1148 TH2F *hPtr[kNmethodsPID]={
1149 (TH2F*)histoContainer->At(AliTRDpidUtil::kLQ),
1150 (TH2F*)histoContainer->At(AliTRDpidUtil::kNN),
1151 (TH2F*)histoContainer->At(AliTRDpidUtil::kESD)
1153 for(Int_t iMom = 0; iMom < fMomentumAxis->GetNbins(); iMom++){
1154 Float_t mom(fMomentumAxis->GetBinCenter(iMom+1));
1156 Int_t binEl(fMomentumAxis->GetNbins() * AliPID::kElectron + iMom + 1),
1157 binXX(fMomentumAxis->GetNbins() * species + iMom + 1);
1159 for(Int_t iMethod = 0; iMethod < kNmethodsPID; iMethod++){
1160 // Calculate the Species Efficiency at electronEfficiency% electron efficiency for each Method
1162 TH1D *histo1 = hPtr[iMethod] -> ProjectionY(Form("%s_el", fgMethod[iMethod]), binEl, binEl);
1163 TH1D *histo2 = hPtr[iMethod] -> ProjectionY(Form("%s_%s", fgMethod[iMethod], AliTRDCalPID::GetPartName(species)), binXX, binXX);
1165 if(!fUtil->CalculatePionEffi(histo1, histo2)) continue;
1167 g=(TGraphErrors*)eff->At(iMethod);
1168 g->SetPoint(iMom, mom, 1.e2*fUtil->GetPionEfficiency());
1169 g->SetPointError(iMom, 0., 1.e2*fUtil->GetError());
1170 AliDebug(2, Form("%s Efficiency for %s is : %f +/- %f", AliTRDCalPID::GetPartName(species), fgMethod[iMethod], fUtil->GetPionEfficiency(), fUtil->GetError()));
1172 if(!thres) continue;
1173 g=(TGraphErrors*)thres->At(iMethod);
1174 g->SetPoint(iMom, mom, fUtil->GetThreshold());
1175 g->SetPointError(iMom, 0., 0.);