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 "AliTRDinfoGen.h"
54 #include "info/AliTRDtrackInfo.h"
55 #include "info/AliTRDpidInfo.h"
56 #include "info/AliTRDv0Info.h"
58 Char_t const * AliTRDcheckPID::fgMethod[3] = {"LQ", "NN", "ESD"};
60 ClassImp(AliTRDcheckPID)
62 //________________________________________________________________________
63 AliTRDcheckPID::AliTRDcheckPID()
70 ,fMinNTracklets(AliTRDgeometry::kNlayer)
71 ,fMaxNTracklets(AliTRDgeometry::kNlayer)
74 // Default constructor
76 SetNameTitle("TRDcheckPID", "Check TRD PID");
80 //________________________________________________________________________
81 AliTRDcheckPID::AliTRDcheckPID(char* name )
82 :AliTRDrecoTask(name, "Check TRD PID")
88 ,fMinNTracklets(AliTRDgeometry::kNlayer)
89 ,fMaxNTracklets(AliTRDgeometry::kNlayer)
92 // Default constructor
98 DefineInput(3, TObjArray::Class()); // v0 list
99 DefineOutput(2, TObjArray::Class()); // pid info list
103 //________________________________________________________________________
104 void AliTRDcheckPID::LocalInit()
106 // Initialize working data
108 // Initialize momentum axis with default values
109 Double_t defaultMomenta[AliTRDCalPID::kNMom+1];
110 for(Int_t imom = 0; imom < AliTRDCalPID::kNMom+1; imom++)
111 defaultMomenta[imom] = AliTRDCalPID::GetMomentumBinning(imom);
112 SetMomentumBinning(AliTRDCalPID::kNMom, defaultMomenta);
114 memset(fEfficiency, 0, AliPID::kSPECIES*sizeof(TObjArray*));
116 fUtil = new AliTRDpidUtil();
119 //________________________________________________________________________
120 AliTRDcheckPID::~AliTRDcheckPID()
122 if(fPID){fPID->Delete(); delete fPID;}
123 if(fGraph){fGraph->Delete(); delete fGraph;}
124 if(fUtil) delete fUtil;
128 //________________________________________________________________________
129 void AliTRDcheckPID::UserCreateOutputObjects()
134 AliTRDrecoTask::UserCreateOutputObjects();
135 fPID = new TObjArray();
136 fPID->SetOwner(kTRUE);
140 //________________________________________________________
141 void AliTRDcheckPID::UserExec(Option_t *opt)
147 fV0s = dynamic_cast<TObjArray *>(GetInputData(3));
150 AliTRDrecoTask::UserExec(opt);
154 //_______________________________________________________
155 TObjArray * AliTRDcheckPID::Histos(){
158 // Create QA histograms
160 if(fContainer) return fContainer;
162 Int_t xBins = AliPID::kSPECIES*fMomentumAxis->GetNbins();
163 fContainer = new TObjArray(); fContainer->Expand(kNPlots);
165 const Float_t epsilon = 1./(2*(AliTRDpidUtil::kBins-1)); // get nice histos with bin center at 0 and 1
168 const Int_t kNmethodsPID=Int_t(sizeof(fgMethod)/sizeof(Char_t*));
169 // histos with posterior probabilities for all particle species
170 for(Int_t is=AliPID::kSPECIES; is--;){
171 fEfficiency[is] = new TObjArray(kNmethodsPID);
172 fEfficiency[is]->SetOwner();
173 fEfficiency[is]->SetName(Form("Eff_%s", AliPID::ParticleShortName(is)));
174 fContainer->AddAt(fEfficiency[is], is?kEfficiencyMu+is-1:kEfficiency);
175 for(Int_t im=kNmethodsPID; im--;){
176 if(!(h = (TH2F*)gROOT->FindObject(Form("PID_%s_%s", fgMethod[im], AliPID::ParticleShortName(is))))) {
177 h = new TH2F(Form("PID_%s_%s", fgMethod[im], AliPID::ParticleShortName(is)), "", xBins, -0.5, xBins - 0.5,
178 AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon);
180 fEfficiency[is]->AddAt(h, im);
184 // histos of the dE/dx distribution for all 5 particle species and 11 momenta
185 if(!(h = (TH2F*)gROOT->FindObject("dEdx"))){
186 h = new TH2F("dEdx", "",
187 xBins, -0.5, xBins - 0.5,
191 fContainer->AddAt(h, kdEdx);
193 // histos of the dE/dx slices for all 5 particle species and 11 momenta
194 if(!(h = (TH2F*)gROOT->FindObject("dEdxSlice"))){
195 h = new TH2F("dEdxSlice", "",
196 xBins*AliTRDpidUtil::kLQslices, -0.5, xBins*AliTRDpidUtil::kLQslices - 0.5,
199 fContainer->AddAt(h, kdEdxSlice);
201 // histos of the pulse height distribution for all 5 particle species and 11 momenta
202 TObjArray *fPH = new TObjArray(2);
203 fPH->SetOwner(); fPH->SetName("PH");
204 fContainer->AddAt(fPH, kPH);
205 if(!(h = (TProfile2D*)gROOT->FindObject("PHT"))){
206 h = new TProfile2D("PHT", "",
207 xBins, -0.5, xBins - 0.5,
208 AliTRDtrackerV1::GetNTimeBins(), -0.5, AliTRDtrackerV1::GetNTimeBins() - 0.5);
211 if(!(h = (TProfile2D*)gROOT->FindObject("PHX"))){
212 h = new TProfile2D("PHX", "",
213 xBins, -0.5, xBins - 0.5,
214 AliTRDtrackerV1::GetNTimeBins(), 0., .5*AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght());
218 // histos of the number of clusters distribution for all 5 particle species and 11 momenta
219 if(!(h = (TH2F*)gROOT->FindObject("NClus"))){
220 h = new TH2F("NClus", "",
221 xBins, -0.5, xBins - 0.5,
224 fContainer->AddAt(h, kNClus);
227 // momentum distributions - absolute and in momentum bins
228 if(!(h = (TH1F*)gROOT->FindObject("hMom"))){
229 h = new TH1F("hMom", "momentum distribution", fMomentumAxis->GetNbins(), fMomentumAxis->GetXmin(), fMomentumAxis->GetXmax());
231 fContainer->AddAt(h, kMomentum);
233 if(!(h = (TH1F*)gROOT->FindObject("hMomBin"))){
234 h = new TH1F("hMomBin", "momentum distribution in momentum bins", fMomentumAxis->GetNbins(), fMomentumAxis->GetXmin(), fMomentumAxis->GetXmax());
236 fContainer->AddAt(h, kMomentumBin);
238 // Number of tracklets per track
239 if(!(h = (TH2F*)gROOT->FindObject("nTracklets"))){
240 h = new TH2F("nTracklets", "",
241 xBins, -0.5, xBins - 0.5,
242 AliTRDgeometry::kNlayer, 0.5, AliTRDgeometry::kNlayer+.5);
244 fContainer->AddAt(h, kNTracklets);
247 if(!(h = (TH1F*)gROOT->FindObject("nV0"))){
248 h = new TH1F("nV0", "V0s/track;v0/track;entries",
251 fContainer->AddAt(h, kV0);
257 //________________________________________________________________________
258 Bool_t AliTRDcheckPID::CheckTrackQuality(const AliTRDtrackV1* track) const
261 // Check if the track is ok for PID
264 Int_t ntracklets = track->GetNumberOfTracklets();
265 if(ntracklets >= fMinNTracklets && ntracklets <= fMaxNTracklets) return 1;
272 //________________________________________________________________________
273 Int_t AliTRDcheckPID::CalcPDG(AliTRDtrackV1* track)
275 // Documentation to come
277 track -> SetReconstructor(AliTRDinfoGen::Reconstructor());
278 (const_cast<AliTRDrecoParam*>(AliTRDinfoGen::Reconstructor()->GetRecoParam()))->SetPIDNeuralNetwork();
281 if(track -> GetPID(AliPID::kElectron) > track -> GetPID(AliPID::kMuon) + track -> GetPID(AliPID::kPion) + track -> GetPID(AliPID::kKaon) + track -> GetPID(AliPID::kProton)){
284 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)){
287 else if(track -> GetPID(AliPID::kKaon) > track -> GetPID(AliPID::kMuon) && track -> GetPID(AliPID::kKaon) > track -> GetPID(AliPID::kPion)){
290 else if(track -> GetPID(AliPID::kMuon) > track -> GetPID(AliPID::kPion)){
299 //_______________________________________________________
300 TH1 *AliTRDcheckPID::PlotLQ(const AliTRDtrackV1 *track)
303 // Plot the probabilities for electrons using 2-dim LQ
307 AliDebug(2, "No ESD info available.");
311 if(track) fkTrack = track;
313 AliDebug(2, "No Track defined.");
318 status = fkESD -> GetStatus();
319 if(!(status&AliESDtrack::kTRDin)) return NULL;
321 if(!CheckTrackQuality(fkTrack)) return NULL;
323 AliTRDtrackV1 cTrack(*fkTrack);
324 cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
327 Float_t momentum = 0.;
330 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
331 pdg = fkMC->GetPDG();
333 //AliWarning("No MC info available!");
334 momentum = cTrack.GetMomentum(0);
335 pdg = CalcPDG(&cTrack);
337 if(!IsInRange(momentum)) return NULL;
339 (const_cast<AliTRDrecoParam*>(AliTRDinfoGen::Reconstructor()->GetRecoParam()))->SetPIDNeuralNetwork(kFALSE);
341 if(cTrack.GetNumberOfTrackletsPID() < fMinNTracklets) return NULL;
342 Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
345 TObjArray *eff(NULL);
346 for(Int_t is=AliPID::kSPECIES; is--;){
347 if(!(eff = dynamic_cast<TObjArray *>(fContainer->At(is?kEfficiencyMu+is-1:kEfficiency)))){
348 AliWarning("No Histogram List defined.");
351 if(!(hPIDLQ = dynamic_cast<TH2F *>(eff->At(AliTRDpidUtil::kLQ)))){
352 AliWarning("No Histogram defined.");
356 hPIDLQ -> Fill(FindBin(species, momentum), cTrack.GetPID(is));
364 //_______________________________________________________
365 TH1 *AliTRDcheckPID::PlotNN(const AliTRDtrackV1 *track)
368 // Plot the probabilities for electrons using 2-dim LQ
372 AliDebug(2, "No ESD info available.");
376 if(track) fkTrack = track;
378 AliDebug(2, "No Track defined.");
383 status = fkESD -> GetStatus();
384 if(!(status&AliESDtrack::kTRDin)) return NULL;
386 if(!CheckTrackQuality(fkTrack)) return NULL;
388 AliTRDtrackV1 cTrack(*fkTrack);
389 cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
392 Float_t momentum = 0.;
394 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
395 pdg = fkMC->GetPDG();
397 //AliWarning("No MC info available!");
398 momentum = cTrack.GetMomentum(0);
399 pdg = CalcPDG(&cTrack);
401 if(!IsInRange(momentum)) return NULL;
403 (const_cast<AliTRDrecoParam*>(AliTRDinfoGen::Reconstructor()->GetRecoParam()))->SetPIDNeuralNetwork();
405 if(cTrack.GetNumberOfTrackletsPID() < fMinNTracklets) return NULL;
407 Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
411 TObjArray *eff(NULL);
412 for(Int_t is=AliPID::kSPECIES; is--;){
413 if(!(eff = dynamic_cast<TObjArray *>(fContainer->At(is?kEfficiencyMu+is-1:kEfficiency)))){
414 AliWarning("No Histogram List defined.");
417 if(!(hPIDNN = dynamic_cast<TH2F *>(eff->At(AliTRDpidUtil::kNN)))){
418 AliWarning("No Histogram defined.");
422 hPIDNN->Fill(FindBin(species, momentum), cTrack.GetPID(is));
429 //_______________________________________________________
430 TH1 *AliTRDcheckPID::PlotESD(const AliTRDtrackV1 *track)
433 // Plot the probabilities for electrons using 2-dim LQ
437 AliDebug(2, "No ESD info available.");
441 if(track) fkTrack = track;
443 AliDebug(2, "No Track defined.");
448 status = fkESD -> GetStatus();
449 if(!(status&AliESDtrack::kTRDin)) return NULL;
451 if(!CheckTrackQuality(fkTrack)) return NULL;
452 if(fkTrack->GetNumberOfTrackletsPID() < fMinNTracklets) return NULL;
455 Float_t momentum = 0.;
457 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
458 pdg = fkMC->GetPDG();
460 //AliWarning("No MC info available!");
461 AliTRDtrackV1 cTrack(*fkTrack);
462 cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
463 momentum = cTrack.GetMomentum(0);
464 pdg = CalcPDG(&cTrack);
466 if(!IsInRange(momentum)) return NULL;
468 const Double32_t *pidESD = fkESD->GetResponseIter();
469 Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
472 TObjArray *eff(NULL);
473 for(Int_t is=AliPID::kSPECIES; is--;){
474 if(!(eff = dynamic_cast<TObjArray *>(fContainer->At(is?kEfficiencyMu+is-1:kEfficiency)))){
475 AliWarning("No Histogram List defined.");
478 if(!(hPID = dynamic_cast<TH2F *>(eff->At(AliTRDpidUtil::kESD)))){
479 AliWarning("No Histogram defined.");
483 hPID->Fill(FindBin(species, momentum), pidESD[is]);
490 //_______________________________________________________
491 TH1 *AliTRDcheckPID::PlotdEdx(const AliTRDtrackV1 *track)
494 // Plot the probabilities for electrons using 2-dim LQ
497 if(track) fkTrack = track;
499 AliDebug(2, "No Track defined.");
503 if(!CheckTrackQuality(fkTrack)) return NULL;
506 if(!(hdEdx = dynamic_cast<TH2F *>(fContainer->At(kdEdx)))){
507 AliWarning("No Histogram defined.");
511 AliTRDtrackV1 cTrack(*fkTrack);
512 cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
514 Float_t momentum = 0.;
516 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
517 pdg = fkMC->GetPDG();
519 //AliWarning("No MC info available!");
520 momentum = cTrack.GetMomentum(0);
521 pdg = CalcPDG(&cTrack);
523 if(!IsInRange(momentum)) return NULL;
525 // Init exchange container
526 Int_t s(AliTRDpidUtil::Pdg2Pid(pdg));
527 AliTRDpidInfo *pid = new AliTRDpidInfo(s);
529 (const_cast<AliTRDrecoParam*>(AliTRDinfoGen::Reconstructor()->GetRecoParam()))->SetPIDNeuralNetwork(kTRUE);
532 Int_t iBin = FindBin(s, momentum);
533 AliTRDseedV1 *tracklet = NULL;
534 for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
535 tracklet = cTrack.GetTracklet(ily);
536 if(!tracklet) continue;
537 tracklet -> CookdEdx(AliTRDpidUtil::kNNslices);
539 // fill exchange container
540 pid->PushBack(tracklet->GetPlane(),
541 AliTRDpidUtil::GetMomentumBin(tracklet->GetMomentum()), tracklet->GetdEdx());
544 for(Int_t i = AliTRDpidUtil::kNNslices; i--;) sumdEdx += tracklet->GetdEdx()[i];
545 sumdEdx /= AliTRDCalPIDNN::kMLPscale;
546 hdEdx -> Fill(iBin, sumdEdx);
554 //_______________________________________________________
555 TH1 *AliTRDcheckPID::PlotdEdxSlice(const AliTRDtrackV1 *track)
558 // Plot the probabilities for electrons using 2-dim LQ
561 if(track) fkTrack = track;
563 AliDebug(2, "No Track defined.");
567 if(!CheckTrackQuality(fkTrack)) return NULL;
570 if(!(hdEdxSlice = dynamic_cast<TH2F *>(fContainer->At(kdEdxSlice)))){
571 AliWarning("No Histogram defined.");
575 AliTRDtrackV1 cTrack(*fkTrack);
576 cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
578 Float_t momentum = 0.;
580 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
581 pdg = fkMC->GetPDG();
583 //AliWarning("No MC info available!");
584 momentum = cTrack.GetMomentum(0);
585 pdg = CalcPDG(&cTrack);
587 if(!IsInRange(momentum)) return NULL;
589 (const_cast<AliTRDrecoParam*>(AliTRDinfoGen::Reconstructor()->GetRecoParam()))->SetPIDNeuralNetwork(kFALSE);
590 Int_t iMomBin = fMomentumAxis->FindBin(momentum);
591 Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
593 AliTRDseedV1 *tracklet = NULL;
594 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
595 tracklet = cTrack.GetTracklet(iChamb);
596 if(!tracklet) continue;
597 tracklet -> CookdEdx(AliTRDpidUtil::kLQslices);
598 fdEdx = const_cast<Float_t *>(tracklet->GetdEdx());
599 for(Int_t iSlice = 0; iSlice < AliTRDpidUtil::kLQslices; iSlice++){
600 hdEdxSlice -> Fill(species * fMomentumAxis->GetNbins() * AliTRDpidUtil::kLQslices + (iMomBin-1) * AliTRDpidUtil::kLQslices + iSlice, fdEdx[iSlice]);
608 //_______________________________________________________
609 TH1 *AliTRDcheckPID::PlotPH(const AliTRDtrackV1 *track)
612 // Plot the probabilities for electrons using 2-dim LQ
615 if(track) fkTrack = track;
617 AliDebug(2, "No Track defined.");
621 if(!CheckTrackQuality(fkTrack)) return NULL;
623 TObjArray *arr = NULL;
624 TProfile2D *hPHX, *hPHT;
625 if(!(arr = dynamic_cast<TObjArray *>(fContainer->At(kPH)))){
626 AliWarning("No Histogram defined.");
629 hPHT = (TProfile2D*)arr->At(0);
630 hPHX = (TProfile2D*)arr->At(1);
633 Float_t momentum = 0.;
635 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
636 pdg = fkMC->GetPDG();
638 //AliWarning("No MC info available!");
639 AliTRDtrackV1 cTrack(*fkTrack);
640 cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
641 momentum = cTrack.GetMomentum(0);
642 pdg = CalcPDG(&cTrack);
644 if(!IsInRange(momentum)) return NULL;;
646 AliTRDseedV1 *tracklet = NULL;
647 AliTRDcluster *cluster = NULL;
648 Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
649 Int_t iBin = FindBin(species, momentum);
650 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
651 tracklet = fkTrack->GetTracklet(iChamb);
652 if(!tracklet) continue;
653 Float_t x0 = tracklet->GetX0();
654 for(Int_t ic = 0; ic < AliTRDtrackerV1::GetNTimeBins(); ic++){
655 if(!(cluster = tracklet->GetClusters(ic))) continue;
656 hPHT -> Fill(iBin, cluster->GetLocalTimeBin(), TMath::Abs(cluster->GetQ()));
657 hPHX -> Fill(iBin, x0 - cluster->GetX(), tracklet->GetdQdl(ic));
664 //_______________________________________________________
665 TH1 *AliTRDcheckPID::PlotNClus(const AliTRDtrackV1 *track)
668 // Plot the probabilities for electrons using 2-dim LQ
671 if(track) fkTrack = track;
673 AliDebug(2, "No Track defined.");
677 if(!CheckTrackQuality(fkTrack)) return NULL;
680 if(!(hNClus = dynamic_cast<TH2F *>(fContainer->At(kNClus)))){
681 AliWarning("No Histogram defined.");
687 Float_t momentum = 0.;
689 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
690 pdg = fkMC->GetPDG();
692 //AliWarning("No MC info available!");
693 AliTRDtrackV1 cTrack(*fkTrack);
694 cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
695 momentum = cTrack.GetMomentum(0);
696 pdg = CalcPDG(&cTrack);
698 if(!IsInRange(momentum)) return NULL;
700 Int_t species = AliTRDpidUtil::AliTRDpidUtil::Pdg2Pid(pdg);
701 Int_t iBin = FindBin(species, momentum);
702 AliTRDseedV1 *tracklet = NULL;
703 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
704 tracklet = fkTrack->GetTracklet(iChamb);
705 if(!tracklet) continue;
706 hNClus -> Fill(iBin, tracklet->GetN());
712 //_______________________________________________________
713 TH1 *AliTRDcheckPID::PlotNTracklets(const AliTRDtrackV1 *track)
716 // Plot the probabilities for electrons using 2-dim LQ
719 if(track) fkTrack = track;
721 AliDebug(2, "No Track defined.");
726 if(!(hTracklets = dynamic_cast<TH2F *>(fContainer->At(kNTracklets)))){
727 AliWarning("No Histogram defined.");
731 AliTRDtrackV1 cTrack(*fkTrack);
732 cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
734 Float_t momentum = 0.;
736 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
737 pdg = fkMC->GetPDG();
739 //AliWarning("No MC info available!");
740 momentum = cTrack.GetMomentum();
741 pdg = CalcPDG(&cTrack);
743 Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
744 if(!IsInRange(momentum)) return NULL;
746 Int_t iBin = FindBin(species, momentum);
747 hTracklets -> Fill(iBin, cTrack.GetNumberOfTracklets());
751 //_______________________________________________________
752 TH1 *AliTRDcheckPID::PlotMom(const AliTRDtrackV1 *track)
755 // Plot the probabilities for electrons using 2-dim LQ
758 if(track) fkTrack = track;
760 AliDebug(2, "No Track defined.");
764 if(!CheckTrackQuality(fkTrack)) return NULL;
767 if(!(hMom = dynamic_cast<TH1F *>(fContainer->At(kMomentum)))){
768 AliWarning("No Histogram defined.");
774 Float_t momentum = 0.;
776 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
777 pdg = fkMC->GetPDG();
779 //AliWarning("No MC info available!");
780 AliTRDtrackV1 cTrack(*fkTrack);
781 cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
782 momentum = cTrack.GetMomentum(0);
783 pdg = CalcPDG(&cTrack);
785 if(IsInRange(momentum)) hMom -> Fill(momentum);
790 //_______________________________________________________
791 TH1 *AliTRDcheckPID::PlotMomBin(const AliTRDtrackV1 *track)
794 // Plot the probabilities for electrons using 2-dim LQ
797 if(track) fkTrack = track;
799 AliDebug(2, "No Track defined.");
803 if(!CheckTrackQuality(fkTrack)) return NULL;
806 if(!(hMomBin = dynamic_cast<TH1F *>(fContainer->At(kMomentumBin)))){
807 AliWarning("No Histogram defined.");
813 Float_t momentum = 0.;
816 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
817 pdg = fkMC->GetPDG();
819 //AliWarning("No MC info available!");
820 AliTRDtrackV1 cTrack(*fkTrack);
821 cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
822 momentum = cTrack.GetMomentum(0);
824 if(IsInRange(momentum)) hMomBin -> Fill(fMomentumAxis->FindBin(momentum));
828 //_______________________________________________________
829 TH1 *AliTRDcheckPID::PlotV0(const AliTRDtrackV1 *track)
832 // Plot the V0 performance against MC
835 if(track) fkTrack = track;
837 AliDebug(2, "No Track defined.");
840 if(!fkESD->HasV0()) return NULL;
842 AliDebug(1, "No MC defined.");
846 AliWarning("No output container defined.");
849 AliDebug(2, Form("TRACK[%d] species[%s][%d]\n", fkESD->GetId(), fkMC->GetPID()>=0?AliPID::ParticleShortName(fkMC->GetPID()):"none", fkMC->GetPDG()));
852 if(!(h = dynamic_cast<TH1F*>(fContainer->At(kV0)))) return NULL;
853 Int_t sgn(0), n(0); AliTRDv0Info *v0(NULL);
854 for(Int_t iv0(fV0s->GetEntriesFast()); iv0--;){
855 if(!(v0=(AliTRDv0Info*)fV0s->At(iv0))) continue;
856 if(!(sgn = v0->HasTrack(fkESD->GetId()))) continue;
857 //for(Int_t is=AliPID::kSPECIES; is--;) v0->GetPID(is, track);
866 //________________________________________________________
867 Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig)
869 // Steering function to retrieve performance plots
871 Bool_t kFIRST = kTRUE;
872 TGraphErrors *g = NULL;
874 TObjArray *arr = NULL;
875 TH1 *h1 = NULL, *h=NULL;
877 TList *content = NULL;
880 gPad->Divide(2, 1, 1.e-5, 1.e-5);
881 TList *l=gPad->GetListOfPrimitives();
882 TVirtualPad *pad = ((TVirtualPad*)l->At(0));pad->cd();
883 pad->SetMargin(0.1, 0.01, 0.1, 0.01);
885 TLegend *legEff = new TLegend(.64, .84, .98, .98);
886 legEff->SetBorderSize(1);legEff->SetTextSize(0.03255879);
887 legEff->SetFillColor(0);
888 h=new TH1S("hEff", "", 1, .5, 11.);
889 h->SetLineColor(1);h->SetLineWidth(1);
891 ax->SetTitle("p [GeV/c]");
892 ax->SetRangeUser(.5, 11.);
893 ax->SetMoreLogLabels();
895 ax->SetTitle("#epsilon_{#pi} [%]");
897 ax->SetRangeUser(1.e-2, 10.);
899 content = (TList *)fGraph->FindObject(Form("Eff_%s", AliTRDCalPID::GetPartName(AliPID::kPion)));
900 if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kLQ))) break;
901 if(!g->GetN()) break;
902 legEff->SetHeader("PID Method [PION]");
903 g->Draw("pc"); legEff->AddEntry(g, "LQ 2D", "pl");
904 if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kNN))) break;
905 g->Draw("pc"); legEff->AddEntry(g, "NN", "pl");
906 if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kESD))) break;
907 g->Draw("p"); legEff->AddEntry(g, "ESD", "pl");
915 pad = ((TVirtualPad*)l->At(1));pad->cd();
916 pad->SetMargin(0.1, 0.01, 0.1, 0.01);
917 h=new TH1S("hThr", "", 1, .5, 11.);
918 h->SetLineColor(1);h->SetLineWidth(1);
920 ax->SetTitle("p [GeV/c]");
921 ax->SetMoreLogLabels();
923 ax->SetTitle("Threshold [%]");
924 ax->SetRangeUser(5.e-2, 1.);
926 content = (TList *)fGraph->FindObject("Thres");
927 if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kLQ))) break;
928 if(!g->GetN()) break;
930 if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kNN))) break;
932 if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kESD))) break;
940 gPad->Divide(2, 1, 1.e-5, 1.e-5);
941 TList *l=gPad->GetListOfPrimitives();
942 TVirtualPad *pad = ((TVirtualPad*)l->At(0));pad->cd();
943 pad->SetMargin(0.1, 0.01, 0.1, 0.01);
945 TLegend *legEff = new TLegend(.64, .84, .98, .98);
946 legEff->SetBorderSize(1);legEff->SetTextSize(0.03255879);
947 legEff->SetFillColor(0);
948 h = (TH1S*)gROOT->FindObject("hEff");
949 h=(TH1S*)h->Clone("hEff_K");
950 h->SetYTitle("#epsilon_{K} [%]");
951 h->GetYaxis()->SetRangeUser(1.e-2, 1.e2);
953 content = (TList *)fGraph->FindObject(Form("Eff_%s", AliTRDCalPID::GetPartName(AliPID::kKaon)));
954 if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kLQ))) break;
955 if(!g->GetN()) break;
956 legEff->SetHeader("PID Method [KAON]");
957 g->Draw("pc"); legEff->AddEntry(g, "LQ 2D", "pl");
958 if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kNN))) break;
959 g->Draw("pc"); legEff->AddEntry(g, "NN", "pl");
960 if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kESD))) break;
961 g->Draw("p"); legEff->AddEntry(g, "ESD", "pl");
968 TLegend *legEff2 = new TLegend(.64, .84, .98, .98);
969 legEff2->SetBorderSize(1);legEff2->SetTextSize(0.03255879);
970 legEff2->SetFillColor(0);
971 pad = ((TVirtualPad*)l->At(1));pad->cd();
972 pad->SetMargin(0.1, 0.01, 0.1, 0.01);
973 h=(TH1S*)h->Clone("hEff_p");
974 h->SetYTitle("#epsilon_{p} [%]");
975 h->GetYaxis()->SetRangeUser(1.e-2, 1.e2);
977 content = (TList *)fGraph->FindObject(Form("Eff_%s", AliTRDCalPID::GetPartName(AliPID::kProton)));
978 if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kLQ))) break;
979 if(!g->GetN()) break;
980 legEff2->SetHeader("PID Method [PROTON]");
981 g->Draw("pc"); legEff2->AddEntry(g, "LQ 2D", "pl");
982 if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kNN))) break;
983 g->Draw("pc"); legEff2->AddEntry(g, "NN", "pl");
984 if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kESD))) break;
985 g->Draw("p"); legEff2->AddEntry(g, "ESD", "pl");
994 // save 2.0 GeV projection as reference
995 TLegend *legdEdx = new TLegend(.7, .7, .98, .98);
996 legdEdx->SetBorderSize(1);
998 if(!(h2 = (TH2F*)(fContainer->At(kdEdx)))) break;
999 legdEdx->SetHeader("Particle Species");
1000 gPad->SetMargin(0.1, 0.01, 0.1, 0.01);
1001 for(Int_t is = AliPID::kSPECIES-1; is>=0; is--){
1002 Int_t bin = FindBin(is, 2.);
1003 h1 = h2->ProjectionY(Form("px%d", is), bin, bin);
1004 if(!h1->GetEntries()) continue;
1005 h1->Scale(1./h1->Integral());
1006 h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
1008 h1->GetXaxis()->SetTitle("dE/dx (a.u.)");
1009 h1->GetYaxis()->SetTitle("<Entries>");
1011 h = (TH1F*)h1->DrawClone(kFIRST ? "c" : "samec");
1012 legdEdx->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "l");
1026 gPad->Divide(2, 1, 1.e-5, 1.e-5);
1027 TList *l=gPad->GetListOfPrimitives();
1029 // save 2.0 GeV projection as reference
1030 TLegend *legPH = new TLegend(.4, .7, .68, .98);
1031 legPH->SetBorderSize(1);legPH->SetFillColor(0);
1032 legPH->SetHeader("Particle Species");
1033 if(!(arr = (TObjArray*)(fContainer->At(kPH)))) break;
1034 if(!(h2 = (TProfile2D*)(arr->At(0)))) break;
1036 TVirtualPad *pad = ((TVirtualPad*)l->At(0));pad->cd();
1037 pad->SetMargin(0.1, 0.01, 0.1, 0.01);
1039 for(Int_t is=0; is<AliPID::kSPECIES; is++){
1040 Int_t bin = FindBin(is, 2.);
1041 h1 = h2->ProjectionY(Form("pyt%d", is), bin, bin);
1042 if(!h1->GetEntries()) continue;
1043 h1->SetMarkerStyle(24);
1044 h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
1045 h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
1047 h1->GetXaxis()->SetTitle("t_{drift} [100*ns]");
1048 h1->GetYaxis()->SetTitle("<dQ/dt> [a.u.]");
1050 h = (TH1F*)h1->DrawClone(kFIRST ? "c" : "samec");
1051 legPH->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "pl");
1055 pad = ((TVirtualPad*)l->At(1));pad->cd();
1056 pad->SetMargin(0.1, 0.01, 0.1, 0.01);
1057 if(!(h2 = (TProfile2D*)(arr->At(1)))) break;
1059 for(Int_t is=0; is<AliPID::kSPECIES; is++){
1060 Int_t bin = FindBin(is, 2.);
1061 h1 = h2->ProjectionY(Form("pyx%d", is), bin, bin);
1062 if(!h1->GetEntries()) continue;
1063 h1->SetMarkerStyle(24);
1064 h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
1065 h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
1067 h1->GetXaxis()->SetTitle("x_{drift} [cm]");
1068 h1->GetYaxis()->SetTitle("<dQ/dl> [a.u./cm]");
1070 h1->DrawClone(kFIRST ? "c" : "samec");
1083 // save 2.0 GeV projection as reference
1084 TLegend *legNClus = new TLegend(.13, .7, .4, .98);
1085 legNClus->SetBorderSize(1);
1086 legNClus->SetFillColor(0);
1089 if(!(h2 = (TH2F*)(fContainer->At(kNClus)))) break;
1090 legNClus->SetHeader("Particle Species");
1091 for(Int_t is=0; is<AliPID::kSPECIES; is++){
1092 Int_t bin = FindBin(is, 2.);
1093 h1 = h2->ProjectionY(Form("pyNClus%d", is), bin, bin);
1094 if(!h1->GetEntries()) continue;
1095 h1->Scale(100./h1->Integral());
1096 //h1->SetMarkerStyle(24);
1097 //h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
1098 h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
1100 h1->GetXaxis()->SetTitle("N^{cl}/tracklet");
1101 h1->GetYaxis()->SetTitle("Prob. [%]");
1102 h = (TH1F*)h1->DrawClone("c");
1104 h->GetXaxis()->SetRangeUser(0., 35.);
1106 } else h = (TH1F*)h1->DrawClone("samec");
1108 legNClus->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "l");
1122 TLegend *legNClus = new TLegend(.4, .7, .68, .98);
1123 legNClus->SetBorderSize(1);
1125 if(!(h2 = (TH2F*)(fContainer->At(kNTracklets)))) break;
1126 legNClus->SetHeader("Particle Species");
1127 for(Int_t is=0; is<AliPID::kSPECIES; is++){
1128 Int_t bin = FindBin(is, 2.);
1129 h1 = h2->ProjectionY(Form("pyNTracklets%d", is), bin, bin);
1130 if(!h1->GetEntries()) continue;
1131 h1->Scale(100./h1->Integral());
1132 //h1->SetMarkerStyle(24);
1133 //h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
1134 h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
1136 h1->GetXaxis()->SetTitle("N^{trklt}/track");
1137 h1->GetXaxis()->SetRangeUser(1.,6.);
1138 h1->GetYaxis()->SetTitle("Prob. [%]");
1140 h = (TH1F*)h1->DrawClone(kFIRST ? "c" : "samec");
1141 legNClus->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "l");
1153 AliInfo(Form("Reference plot [%d] missing result", ifig));
1157 //________________________________________________________________________
1158 Bool_t AliTRDcheckPID::PostProcess()
1160 // Draw result to the screen
1161 // Called once at the end of the query
1164 Printf("ERROR: list not available");
1168 TObjArray *eff(NULL);
1170 fGraph = new TObjArray(2*AliPID::kSPECIES);
1173 if(!(eff = dynamic_cast<TObjArray *>(fContainer->At(kEfficiency)))){
1174 AliError("Efficiency container for Electrons missing.");
1177 EvaluateEfficiency(eff, fGraph, AliPID::kPion, 0.9);
1178 EvaluateEfficiency(eff, fGraph, AliPID::kKaon, 0.9);
1179 EvaluateEfficiency(eff, fGraph, AliPID::kProton, 0.9);
1185 //________________________________________________________________________
1186 void AliTRDcheckPID::EvaluateEfficiency(const TObjArray * const histoContainer, TObjArray *results, Int_t species, Float_t electronEfficiency){
1187 // Process PID information for pion efficiency
1189 fUtil->SetElectronEfficiency(electronEfficiency);
1192 const Int_t kNmethodsPID=Int_t(sizeof(fgMethod)/sizeof(Char_t*));
1193 Color_t colors[kNmethodsPID] = {kBlue, kGreen+2, kRed};
1194 Int_t markerStyle[kNmethodsPID] = {7, 7, 24};
1195 // efficiency graphs
1196 TGraphErrors *g(NULL);
1197 TObjArray *eff = new TObjArray(kNmethodsPID); eff->SetOwner(); eff->SetName(Form("Eff_%s", AliTRDCalPID::GetPartName(species)));
1198 results->AddAt(eff, species);
1199 for(Int_t iMethod = 0; iMethod < kNmethodsPID; iMethod++){
1200 eff->AddAt(g = new TGraphErrors(), iMethod);
1201 g->SetName(Form("%s", fgMethod[iMethod]));
1202 g->SetLineColor(colors[iMethod]);
1203 g->SetMarkerColor(colors[iMethod]);
1204 g->SetMarkerStyle(markerStyle[iMethod]);
1207 // Threshold graphs if not already
1208 TObjArray *thres(NULL);
1209 if(!(results->At(AliPID::kSPECIES))){
1210 thres = new TObjArray(kNmethodsPID); thres->SetOwner();
1211 thres->SetName("Thres");
1212 results->AddAt(thres, AliPID::kSPECIES);
1213 for(Int_t iMethod = 0; iMethod < kNmethodsPID; iMethod++){
1214 thres->AddAt(g = new TGraphErrors(), iMethod);
1215 g->SetName(Form("%s", fgMethod[iMethod]));
1216 g->SetLineColor(colors[iMethod]);
1217 g->SetMarkerColor(colors[iMethod]);
1218 g->SetMarkerStyle(markerStyle[iMethod]);
1222 TH2F *hPtr[kNmethodsPID]={
1223 (TH2F*)histoContainer->At(AliTRDpidUtil::kLQ),
1224 (TH2F*)histoContainer->At(AliTRDpidUtil::kNN),
1225 (TH2F*)histoContainer->At(AliTRDpidUtil::kESD)
1227 for(Int_t iMom = 0; iMom < fMomentumAxis->GetNbins(); iMom++){
1228 Float_t mom(fMomentumAxis->GetBinCenter(iMom+1));
1230 Int_t binEl(fMomentumAxis->GetNbins() * AliPID::kElectron + iMom + 1),
1231 binXX(fMomentumAxis->GetNbins() * species + iMom + 1);
1233 for(Int_t iMethod = 0; iMethod < kNmethodsPID; iMethod++){
1234 // Calculate the Species Efficiency at electronEfficiency% electron efficiency for each Method
1236 TH1D *histo1 = hPtr[iMethod] -> ProjectionY(Form("%s_el", fgMethod[iMethod]), binEl, binEl);
1237 TH1D *histo2 = hPtr[iMethod] -> ProjectionY(Form("%s_%s", fgMethod[iMethod], AliTRDCalPID::GetPartName(species)), binXX, binXX);
1239 if(!fUtil->CalculatePionEffi(histo1, histo2)) continue;
1241 g=(TGraphErrors*)eff->At(iMethod);
1242 g->SetPoint(iMom, mom, 1.e2*fUtil->GetPionEfficiency());
1243 g->SetPointError(iMom, 0., 1.e2*fUtil->GetError());
1244 AliDebug(2, Form("%s Efficiency for %s is : %f +/- %f", AliTRDCalPID::GetPartName(species), fgMethod[iMethod], fUtil->GetPionEfficiency(), fUtil->GetError()));
1246 if(!thres) continue;
1247 g=(TGraphErrors*)thres->At(iMethod);
1248 g->SetPoint(iMom, mom, fUtil->GetThreshold());
1249 g->SetPointError(iMom, 0., 0.);