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 "AliAnalysisManager.h"
55 #include "info/AliTRDtrackInfo.h"
56 #include "info/AliTRDpidInfo.h"
57 #include "info/AliTRDv0Info.h"
59 Char_t const * AliTRDcheckPID::fgMethod[3] = {"LQ", "NN", "ESD"};
61 ClassImp(AliTRDcheckPID)
63 //________________________________________________________________________
64 AliTRDcheckPID::AliTRDcheckPID()
71 ,fMinNTracklets(AliTRDgeometry::kNlayer)
72 ,fMaxNTracklets(AliTRDgeometry::kNlayer)
75 // Default constructor
77 SetNameTitle("TRDcheckPID", "Check TRD PID");
81 //________________________________________________________________________
82 AliTRDcheckPID::AliTRDcheckPID(char* name )
83 :AliTRDrecoTask(name, "Check TRD PID")
89 ,fMinNTracklets(AliTRDgeometry::kNlayer)
90 ,fMaxNTracklets(AliTRDgeometry::kNlayer)
93 // Default constructor
99 DefineInput(3, TObjArray::Class()); // v0 list
100 DefineOutput(2, TObjArray::Class()); // pid info list
104 //________________________________________________________________________
105 void AliTRDcheckPID::LocalInit()
107 // Initialize working data
109 // Initialize momentum axis with default values
110 Double_t defaultMomenta[AliTRDCalPID::kNMom+1];
111 for(Int_t imom = 0; imom < AliTRDCalPID::kNMom+1; imom++)
112 defaultMomenta[imom] = AliTRDCalPID::GetMomentumBinning(imom);
113 SetMomentumBinning(AliTRDCalPID::kNMom, defaultMomenta);
115 memset(fEfficiency, 0, AliPID::kSPECIES*sizeof(TObjArray*));
117 fUtil = new AliTRDpidUtil();
120 //________________________________________________________________________
121 AliTRDcheckPID::~AliTRDcheckPID()
123 AliAnalysisManager* amg = AliAnalysisManager::GetAnalysisManager();
124 if (amg && amg->IsProofMode()) return;
125 if(fPID){fPID->Delete(); delete fPID;}
126 if(fGraph){fGraph->Delete(); delete fGraph;}
127 if(fUtil) delete fUtil;
131 //________________________________________________________________________
132 void AliTRDcheckPID::UserCreateOutputObjects()
137 AliTRDrecoTask::UserCreateOutputObjects();
138 fPID = new TObjArray();
139 fPID->SetOwner(kTRUE);
143 //________________________________________________________
144 void AliTRDcheckPID::UserExec(Option_t *opt)
150 fV0s = dynamic_cast<TObjArray *>(GetInputData(3));
153 AliTRDrecoTask::UserExec(opt);
157 //_______________________________________________________
158 TObjArray * AliTRDcheckPID::Histos(){
161 // Create QA histograms
163 if(fContainer) return fContainer;
165 Int_t xBins = AliPID::kSPECIES*fMomentumAxis->GetNbins();
166 fContainer = new TObjArray(); fContainer->Expand(kNPlots); fContainer->SetOwner(kTRUE);
168 const Float_t epsilon = 1./(2*(AliTRDpidUtil::kBins-1)); // get nice histos with bin center at 0 and 1
171 const Int_t kNmethodsPID=Int_t(sizeof(fgMethod)/sizeof(Char_t*));
172 // histos with posterior probabilities for all particle species
173 for(Int_t is=AliPID::kSPECIES; is--;){
174 fEfficiency[is] = new TObjArray(kNmethodsPID);
175 fEfficiency[is]->SetOwner();
176 fEfficiency[is]->SetName(Form("Eff_%s", AliPID::ParticleShortName(is)));
177 fContainer->AddAt(fEfficiency[is], is?kEfficiencyMu+is-1:kEfficiency);
178 for(Int_t im=kNmethodsPID; im--;){
179 if(!(h = (TH2F*)gROOT->FindObject(Form("PID_%s_%s", fgMethod[im], AliPID::ParticleShortName(is))))) {
180 h = new TH2F(Form("PID_%s_%s", fgMethod[im], AliPID::ParticleShortName(is)), "", xBins, -0.5, xBins - 0.5,
181 AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon);
183 fEfficiency[is]->AddAt(h, im);
187 // histos of the dE/dx distribution for all 5 particle species and 11 momenta
188 if(!(h = (TH2F*)gROOT->FindObject("dEdx"))){
189 h = new TH2F("dEdx", "",
190 xBins, -0.5, xBins - 0.5,
194 fContainer->AddAt(h, kdEdx);
196 // histos of the dE/dx slices for all 5 particle species and 11 momenta
197 if(!(h = (TH2F*)gROOT->FindObject("dEdxSlice"))){
198 h = new TH2F("dEdxSlice", "",
199 xBins*AliTRDpidUtil::kNNslices, -0.5, xBins*AliTRDpidUtil::kNNslices - 0.5,
202 fContainer->AddAt(h, kdEdxSlice);
204 // histos of the pulse height distribution for all 5 particle species and 11 momenta
205 TObjArray *fPH = new TObjArray(2);
206 fPH->SetOwner(); fPH->SetName("PH");
207 fContainer->AddAt(fPH, kPH);
208 if(!(h = (TProfile2D*)gROOT->FindObject("PHT"))){
209 h = new TProfile2D("PHT", "<PH>(tb);p*species;tb [100*ns];entries",
210 xBins, -0.5, xBins - 0.5,
211 AliTRDseedV1::kNtb, -0.5, AliTRDseedV1::kNtb - 0.5);
214 if(!(h = (TProfile2D*)gROOT->FindObject("PHX"))){
215 h = new TProfile2D("PHX", "<PH>(x);p*species;x_{drift} [cm];entries",
216 xBins, -0.5, xBins - 0.5,
221 // histos of the number of clusters distribution for all 5 particle species and 11 momenta
222 if(!(h = (TH2F*)gROOT->FindObject("NClus"))){
223 h = new TH2F("NClus", "",
224 xBins, -0.5, xBins - 0.5,
227 fContainer->AddAt(h, kNClus);
230 // momentum distributions - absolute and in momentum bins
231 if(!(h = (TH1F*)gROOT->FindObject("hMom"))){
232 h = new TH1F("hMom", "momentum distribution", fMomentumAxis->GetNbins(), fMomentumAxis->GetXmin(), fMomentumAxis->GetXmax());
234 fContainer->AddAt(h, kMomentum);
236 if(!(h = (TH1F*)gROOT->FindObject("hMomBin"))){
237 h = new TH1F("hMomBin", "momentum distribution in momentum bins", fMomentumAxis->GetNbins(), fMomentumAxis->GetXmin(), fMomentumAxis->GetXmax());
239 fContainer->AddAt(h, kMomentumBin);
241 // Number of tracklets per track
242 if(!(h = (TH2F*)gROOT->FindObject("nTracklets"))){
243 h = new TH2F("nTracklets", "",
244 xBins, -0.5, xBins - 0.5,
245 AliTRDgeometry::kNlayer, 0.5, AliTRDgeometry::kNlayer+.5);
247 fContainer->AddAt(h, kNTracklets);
250 if(!(h = (TH1F*)gROOT->FindObject("nV0"))){
251 h = new TH1F("nV0", "V0s/track;v0/track;entries",
254 fContainer->AddAt(h, kV0);
256 // dQ/dl for 1D-Likelihood
257 if(!(h = (TH1F *)gROOT->FindObject("dQdl"))){
258 h = new TH2F("dQdl", "dQ/dl per layer;p*species;dQ/dl [a.u.]", xBins, -0.5, xBins - 0.5, 800, 0., 40000.);
260 fContainer->AddAt(h, kdQdl);
266 //________________________________________________________________________
267 Bool_t AliTRDcheckPID::CheckTrackQuality(const AliTRDtrackV1* track) const
270 // Check if the track is ok for PID
273 Int_t ntracklets = track->GetNumberOfTracklets();
274 if(ntracklets >= fMinNTracklets && ntracklets <= fMaxNTracklets) return 1;
281 //________________________________________________________________________
282 Int_t AliTRDcheckPID::CalcPDG(AliTRDtrackV1* track)
284 // Documentation to come
286 /* track -> SetReconstructor(AliTRDinfoGen::Reconstructor());
287 (const_cast<AliTRDrecoParam*>(AliTRDinfoGen::Reconstructor()->GetRecoParam()))->SetPIDNeuralNetwork();
288 track -> CookPID();*/
290 if(track -> GetPID(AliPID::kElectron) > track -> GetPID(AliPID::kMuon) + track -> GetPID(AliPID::kPion) + track -> GetPID(AliPID::kKaon) + track -> GetPID(AliPID::kProton)){
293 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)){
296 else if(track -> GetPID(AliPID::kKaon) > track -> GetPID(AliPID::kMuon) && track -> GetPID(AliPID::kKaon) > track -> GetPID(AliPID::kPion)){
299 else if(track -> GetPID(AliPID::kMuon) > track -> GetPID(AliPID::kPion)){
308 //_______________________________________________________
309 TH1 *AliTRDcheckPID::PlotLQ(const AliTRDtrackV1 *track)
312 // Plot the probabilities for electrons using 2-dim LQ
316 AliDebug(2, "No ESD info available.");
320 if(track) fkTrack = track;
322 AliDebug(2, "No Track defined.");
327 status = fkESD -> GetStatus();
328 if(!(status&AliESDtrack::kTRDin)) return NULL;
330 if(!CheckTrackQuality(fkTrack)) return NULL;
332 AliTRDtrackV1 cTrack(*fkTrack);
333 cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
336 Float_t momentum = 0.;
339 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
340 pdg = fkMC->GetPDG();
342 //AliWarning("No MC info available!");
343 momentum = cTrack.GetMomentum(0);
344 pdg = CalcPDG(&cTrack);
346 if(!IsInRange(momentum)) return NULL;
348 (const_cast<AliTRDrecoParam*>(AliTRDinfoGen::Reconstructor()->GetRecoParam()))->SetPIDNeuralNetwork(kFALSE);
350 if(cTrack.GetNumberOfTrackletsPID() < fMinNTracklets) return NULL;
351 Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
354 TObjArray *eff(NULL);
355 for(Int_t is=AliPID::kSPECIES; is--;){
356 if(!(eff = dynamic_cast<TObjArray *>(fContainer->At(is?kEfficiencyMu+is-1:kEfficiency)))){
357 AliWarning("No Histogram List defined.");
360 if(!(hPIDLQ = dynamic_cast<TH2F *>(eff->At(AliTRDpidUtil::kLQ)))){
361 AliWarning("No Histogram defined.");
365 hPIDLQ -> Fill(FindBin(species, momentum), cTrack.GetPID(is));
373 //_______________________________________________________
374 TH1 *AliTRDcheckPID::PlotNN(const AliTRDtrackV1 *track)
377 // Plot the probabilities for electrons using 2-dim LQ
381 AliDebug(2, "No ESD info available.");
385 if(track) fkTrack = track;
387 AliDebug(2, "No Track defined.");
392 status = fkESD -> GetStatus();
393 if(!(status&AliESDtrack::kTRDin)) return NULL;
395 if(!CheckTrackQuality(fkTrack)) return NULL;
397 AliTRDtrackV1 cTrack(*fkTrack);
398 cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
401 Float_t momentum = 0.;
403 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
404 pdg = fkMC->GetPDG();
406 //AliWarning("No MC info available!");
407 momentum = cTrack.GetMomentum(0);
408 pdg = CalcPDG(&cTrack);
410 if(!IsInRange(momentum)) return NULL;
412 (const_cast<AliTRDrecoParam*>(AliTRDinfoGen::Reconstructor()->GetRecoParam()))->SetPIDNeuralNetwork();
414 if(cTrack.GetNumberOfTrackletsPID() < fMinNTracklets) return NULL;
416 Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
420 TObjArray *eff(NULL);
421 for(Int_t is=AliPID::kSPECIES; is--;){
422 if(!(eff = dynamic_cast<TObjArray *>(fContainer->At(is?kEfficiencyMu+is-1:kEfficiency)))){
423 AliWarning("No Histogram List defined.");
426 if(!(hPIDNN = dynamic_cast<TH2F *>(eff->At(AliTRDpidUtil::kNN)))){
427 AliWarning("No Histogram defined.");
431 hPIDNN->Fill(FindBin(species, momentum), cTrack.GetPID(is));
438 //_______________________________________________________
439 TH1 *AliTRDcheckPID::PlotESD(const AliTRDtrackV1 *track)
442 // Plot the probabilities for electrons using 2-dim LQ
446 AliDebug(2, "No ESD info available.");
450 if(track) fkTrack = track;
452 AliDebug(2, "No Track defined.");
457 status = fkESD -> GetStatus();
458 if(!(status&AliESDtrack::kTRDin)) return NULL;
460 if(!CheckTrackQuality(fkTrack)) return NULL;
461 if(fkTrack->GetNumberOfTrackletsPID() < fMinNTracklets) return NULL;
464 Float_t momentum = 0.;
466 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
467 pdg = fkMC->GetPDG();
469 //AliWarning("No MC info available!");
470 AliTRDtrackV1 cTrack(*fkTrack);
471 cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
472 momentum = cTrack.GetMomentum(0);
473 pdg = CalcPDG(&cTrack);
475 if(!IsInRange(momentum)) return NULL;
477 const Double32_t *pidESD = fkESD->GetResponseIter();
478 Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
481 TObjArray *eff(NULL);
482 for(Int_t is=AliPID::kSPECIES; is--;){
483 if(!(eff = dynamic_cast<TObjArray *>(fContainer->At(is?kEfficiencyMu+is-1:kEfficiency)))){
484 AliWarning("No Histogram List defined.");
487 if(!(hPID = dynamic_cast<TH2F *>(eff->At(AliTRDpidUtil::kESD)))){
488 AliWarning("No Histogram defined.");
492 hPID->Fill(FindBin(species, momentum), pidESD[is]);
499 //_______________________________________________________
500 TH1 *AliTRDcheckPID::PlotdQdl(const AliTRDtrackV1 *track){
502 // Plot the total charge for the 1D Likelihood method
504 if(track) fkTrack = track;
506 AliDebug(2, "No Track defined");
509 TH2 *hdQdl = dynamic_cast<TH2F *>(fContainer->At(kdQdl));
511 AliWarning("No Histogram defined");
515 if(!CheckTrackQuality(fkTrack)) return NULL;
518 Float_t momentum = 0.;
519 AliTRDtrackV1 cTrack(*fkTrack);
521 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
522 pdg = fkMC->GetPDG();
524 //AliWarning("No MC info available!");
525 momentum = cTrack.GetMomentum(0);
526 pdg = CalcPDG(&cTrack);
528 if(!IsInRange(momentum)) return NULL;
530 // Init exchange container
531 Int_t s(AliTRDpidUtil::Pdg2Pid(pdg));
532 Int_t ibin = FindBin(s, momentum);
534 AliTRDseedV1 *tracklet = NULL;
535 for(Int_t iseed = 0; iseed < 6; iseed++){
536 if(!((tracklet = fkTrack->GetTracklet(iseed)) && tracklet->IsOK())) continue;
537 hdQdl->Fill(ibin, tracklet->GetdQdl());
542 //_______________________________________________________
543 TH1 *AliTRDcheckPID::PlotdEdx(const AliTRDtrackV1 *track)
546 // Plot the probabilities for electrons using 2-dim LQ
549 if(track) fkTrack = track;
551 AliDebug(2, "No Track defined.");
555 if(!CheckTrackQuality(fkTrack)) return NULL;
558 if(!(hdEdx = dynamic_cast<TH2F *>(fContainer->At(kdEdx)))){
559 AliWarning("No Histogram defined.");
563 AliTRDtrackV1 cTrack(*fkTrack);
564 cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
566 Float_t momentum = 0.;
568 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
569 pdg = fkMC->GetPDG();
571 //AliWarning("No MC info available!");
572 momentum = cTrack.GetMomentum(0);
573 pdg = CalcPDG(&cTrack);
575 if(!IsInRange(momentum)) return NULL;
577 // Init exchange container
578 Int_t s(AliTRDpidUtil::Pdg2Pid(pdg));
579 AliTRDpidInfo *pid = new AliTRDpidInfo(s);
581 //(const_cast<AliTRDrecoParam*>(AliTRDinfoGen::Reconstructor()->GetRecoParam()))->SetPIDNeuralNetwork(kTRUE);
583 // Float_t sumdEdx(0.);
584 Int_t iBin = FindBin(s, momentum);
585 AliTRDseedV1 *tracklet = NULL;
586 for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
587 tracklet = cTrack.GetTracklet(ily);
588 if(!tracklet) continue;
589 //tracklet -> CookdEdx(AliTRDpidUtil::kNNslices);
591 // fill exchange container
592 pid->PushBack(tracklet->GetPlane(),
593 AliTRDpidUtil::GetMomentumBin(tracklet->GetMomentum()), tracklet->GetdEdx());
596 for(Int_t i = AliTRDpidUtil::kNNslices; i--;) sumdEdx += tracklet->GetdEdx()[i];
597 sumdEdx /= AliTRDCalPIDNN::kMLPscale;*/
598 hdEdx -> Fill(iBin, tracklet->GetdQdl());
606 //_______________________________________________________
607 TH1 *AliTRDcheckPID::PlotdEdxSlice(const AliTRDtrackV1 *track)
610 // Plot the probabilities for electrons using 2-dim LQ
613 if(track) fkTrack = track;
615 AliDebug(2, "No Track defined.");
619 if(!CheckTrackQuality(fkTrack)) return NULL;
622 if(!(hdEdxSlice = dynamic_cast<TH2F*>(fContainer->At(kdEdxSlice)))){
623 AliWarning("No Histogram defined.");
627 AliTRDtrackV1 cTrack(*fkTrack);
628 cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
630 Float_t momentum = 0.;
632 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
633 pdg = fkMC->GetPDG();
635 //AliWarning("No MC info available!");
636 momentum = cTrack.GetMomentum(0);
637 pdg = CalcPDG(&cTrack);
639 if(!IsInRange(momentum)) return NULL;
641 // (const_cast<AliTRDrecoParam*>(AliTRDinfoGen::Reconstructor()->GetRecoParam()))->SetPIDNeuralNetwork(kFALSE);
642 Int_t iMomBin(-1);/* = fMomentumAxis->FindBin(momentum);*/
643 Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
644 Double32_t *dedxIter = fkESD->GetSliceIter(), *momIter = &dedxIter[AliTRDpidUtil::kNNslices*AliTRDgeometry::kNlayer];
647 while(islice<fkESD->GetNSlices()){
648 printf(" slice[%2d] = %f\n", islice, dedxIter[islice]);
651 for(Int_t ily(0); ily < AliTRDgeometry::kNlayer; ily++){
652 iMomBin = fMomentumAxis->FindBin(momIter[ily]);
653 for(Int_t islice(0); islice < AliTRDpidUtil::kNNslices; islice++){
654 hdEdxSlice -> Fill(species * fMomentumAxis->GetNbins() * AliTRDpidUtil::kNNslices + (iMomBin-1) * AliTRDpidUtil::kNNslices + islice, dedxIter[AliTRDpidUtil::kNNslices*ily+islice]);
658 // AliTRDseedV1 *tracklet = NULL;
659 // for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
660 // tracklet = cTrack.GetTracklet(iChamb);
661 // if(!tracklet) continue;
662 // // tracklet -> CookdEdx(AliTRDpidUtil::kNNslices);
663 // fdEdx = const_cast<Float_t *>(tracklet->GetdEdx());
664 // for(Int_t iSlice = 0; iSlice < AliTRDpidUtil::kNNslices; iSlice++){
665 // hdEdxSlice -> Fill(species * fMomentumAxis->GetNbins() * AliTRDpidUtil::kNNslices + (iMomBin-1) * AliTRDpidUtil::kNNslices + iSlice, fdEdx[iSlice]);
673 //_______________________________________________________
674 TH1 *AliTRDcheckPID::PlotPH(const AliTRDtrackV1 *track)
677 // Plot the probabilities for electrons using 2-dim LQ
680 if(track) fkTrack = track;
682 AliDebug(2, "No Track defined.");
686 if(!CheckTrackQuality(fkTrack)) return NULL;
688 TObjArray *arr = NULL;
689 TProfile2D *hPHX, *hPHT;
690 if(!(arr = dynamic_cast<TObjArray *>(fContainer->At(kPH)))){
691 AliWarning("No Histogram defined.");
694 hPHT = (TProfile2D*)arr->At(0);
695 hPHX = (TProfile2D*)arr->At(1);
698 Float_t momentum = 0.;
700 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
701 pdg = fkMC->GetPDG();
703 //AliWarning("No MC info available!");
704 AliTRDtrackV1 cTrack(*fkTrack);
705 cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
706 momentum = cTrack.GetMomentum(0);
707 pdg = CalcPDG(&cTrack);
709 if(!IsInRange(momentum)) return NULL;;
711 AliTRDseedV1 *tracklet = NULL;
712 AliTRDcluster *cluster = NULL;
713 Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
714 Int_t iBin = FindBin(species, momentum);
715 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
716 tracklet = fkTrack->GetTracklet(iChamb);
717 if(!tracklet) continue;
718 Float_t x0 = tracklet->GetX0();
719 for(Int_t ic = 0; ic < AliTRDseedV1::kNclusters; ic++){
720 if(!(cluster = tracklet->GetClusters(ic))) continue;
721 hPHT -> Fill(iBin, cluster->GetLocalTimeBin(), TMath::Abs(cluster->GetQ()));
722 if(ic<AliTRDseedV1::kNtb) hPHX -> Fill(iBin, x0 - cluster->GetX(), tracklet->GetdQdl(ic));
729 //_______________________________________________________
730 TH1 *AliTRDcheckPID::PlotNClus(const AliTRDtrackV1 *track)
733 // Plot the probabilities for electrons using 2-dim LQ
736 if(track) fkTrack = track;
738 AliDebug(2, "No Track defined.");
742 if(!CheckTrackQuality(fkTrack)) return NULL;
745 if(!(hNClus = dynamic_cast<TH2F *>(fContainer->At(kNClus)))){
746 AliWarning("No Histogram defined.");
752 Float_t momentum = 0.;
754 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
755 pdg = fkMC->GetPDG();
757 //AliWarning("No MC info available!");
758 AliTRDtrackV1 cTrack(*fkTrack);
759 cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
760 momentum = cTrack.GetMomentum(0);
761 pdg = CalcPDG(&cTrack);
763 if(!IsInRange(momentum)) return NULL;
765 Int_t species = AliTRDpidUtil::AliTRDpidUtil::Pdg2Pid(pdg);
766 Int_t iBin = FindBin(species, momentum);
767 AliTRDseedV1 *tracklet = NULL;
768 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
769 tracklet = fkTrack->GetTracklet(iChamb);
770 if(!tracklet) continue;
771 hNClus -> Fill(iBin, tracklet->GetN());
777 //_______________________________________________________
778 TH1 *AliTRDcheckPID::PlotNTracklets(const AliTRDtrackV1 *track)
781 // Plot the probabilities for electrons using 2-dim LQ
784 if(track) fkTrack = track;
786 AliDebug(2, "No Track defined.");
791 if(!(hTracklets = dynamic_cast<TH2F *>(fContainer->At(kNTracklets)))){
792 AliWarning("No Histogram defined.");
796 AliTRDtrackV1 cTrack(*fkTrack);
797 cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
799 Float_t momentum = 0.;
801 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
802 pdg = fkMC->GetPDG();
804 //AliWarning("No MC info available!");
805 momentum = cTrack.GetMomentum();
806 pdg = CalcPDG(&cTrack);
808 Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
809 if(!IsInRange(momentum)) return NULL;
811 Int_t iBin = FindBin(species, momentum);
812 hTracklets -> Fill(iBin, cTrack.GetNumberOfTracklets());
816 //_______________________________________________________
817 TH1 *AliTRDcheckPID::PlotMom(const AliTRDtrackV1 *track)
820 // Plot the probabilities for electrons using 2-dim LQ
823 if(track) fkTrack = track;
825 AliDebug(2, "No Track defined.");
829 if(!CheckTrackQuality(fkTrack)) return NULL;
832 if(!(hMom = dynamic_cast<TH1F *>(fContainer->At(kMomentum)))){
833 AliWarning("No Histogram defined.");
839 Float_t momentum = 0.;
841 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
842 //pdg = fkMC->GetPDG();
844 //AliWarning("No MC info available!");
845 AliTRDtrackV1 cTrack(*fkTrack);
846 cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
847 momentum = cTrack.GetMomentum(0);
848 //pdg = CalcPDG(&cTrack);
850 if(IsInRange(momentum)) hMom -> Fill(momentum);
855 //_______________________________________________________
856 TH1 *AliTRDcheckPID::PlotMomBin(const AliTRDtrackV1 *track)
859 // Plot the probabilities for electrons using 2-dim LQ
862 if(track) fkTrack = track;
864 AliDebug(2, "No Track defined.");
868 if(!CheckTrackQuality(fkTrack)) return NULL;
871 if(!(hMomBin = dynamic_cast<TH1F *>(fContainer->At(kMomentumBin)))){
872 AliWarning("No Histogram defined.");
878 Float_t momentum = 0.;
881 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
882 //pdg = fkMC->GetPDG();
884 //AliWarning("No MC info available!");
885 AliTRDtrackV1 cTrack(*fkTrack);
886 cTrack.SetReconstructor(AliTRDinfoGen::Reconstructor());
887 momentum = cTrack.GetMomentum(0);
889 if(IsInRange(momentum)) hMomBin -> Fill(fMomentumAxis->FindBin(momentum));
893 //_______________________________________________________
894 TH1 *AliTRDcheckPID::PlotV0(const AliTRDtrackV1 *track)
897 // Plot the V0 performance against MC
900 if(track) fkTrack = track;
902 AliDebug(2, "No Track defined.");
905 if(!fkESD->HasV0()) return NULL;
907 AliDebug(1, "No MC defined.");
911 AliWarning("No output container defined.");
914 AliDebug(2, Form("TRACK[%d] species[%s][%d]\n", fkESD->GetId(), fkMC->GetPID()>=0?AliPID::ParticleShortName(fkMC->GetPID()):"none", fkMC->GetPDG()));
917 if(!(h = dynamic_cast<TH1F*>(fContainer->At(kV0)))) return NULL;
918 Int_t sgn(0), n(0); AliTRDv0Info *v0(NULL);
919 for(Int_t iv0(fV0s->GetEntriesFast()); iv0--;){
920 if(!(v0=(AliTRDv0Info*)fV0s->At(iv0))) continue;
921 if(!(sgn = v0->HasTrack(fkESD->GetId()))) continue;
922 //for(Int_t is=AliPID::kSPECIES; is--;) v0->GetPID(is, track);
931 //________________________________________________________
932 Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig)
934 // Steering function to retrieve performance plots
936 Bool_t kFIRST = kTRUE;
937 TGraphErrors *g = NULL;
939 TObjArray *arr = NULL;
940 TH1 *h1 = NULL, *h=NULL;
942 TList *content = NULL;
945 /* gPad->Divide(2, 1, 1.e-5, 1.e-5);
946 TList *l=gPad->GetListOfPrimitives();
947 TVirtualPad *pad = ((TVirtualPad*)l->At(0));pad->cd();
948 pad->SetMargin(0.1, 0.01, 0.1, 0.01);
950 TLegend *legEff = new TLegend(.64, .84, .98, .98);
951 legEff->SetBorderSize(1);legEff->SetTextSize(0.03255879);
952 legEff->SetFillColor(0);
953 h=new TH1S("hEff", "", 1, .5, 11.);
954 h->SetLineColor(1);h->SetLineWidth(1);
956 ax->SetTitle("p [GeV/c]");
957 ax->SetRangeUser(.5, 11.);
958 ax->SetMoreLogLabels();
960 ax->SetTitle("#epsilon_{#pi} [%]");
962 ax->SetRangeUser(1.e-2, 10.);
964 content = (TList *)fGraph->FindObject(Form("Eff_%s", AliTRDCalPID::GetPartName(AliPID::kPion)));
965 if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kLQ))) break;
966 if(!g->GetN()) break;
967 legEff->SetHeader("PID Method [PION]");
968 g->Draw("pc"); legEff->AddEntry(g, "LQ 2D", "pl");
969 if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kNN))) break;
970 g->Draw("pc"); legEff->AddEntry(g, "NN", "pl");
971 if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kESD))) break;
972 g->Draw("p"); legEff->AddEntry(g, "ESD", "pl");
980 pad = ((TVirtualPad*)l->At(1));pad->cd();
981 pad->SetMargin(0.1, 0.01, 0.1, 0.01);
982 h=new TH1S("hThr", "", 1, .5, 11.);
983 h->SetLineColor(1);h->SetLineWidth(1);
985 ax->SetTitle("p [GeV/c]");
986 ax->SetMoreLogLabels();
988 ax->SetTitle("Threshold [%]");
989 ax->SetRangeUser(5.e-2, 1.);
991 content = (TList *)fGraph->FindObject("Thres");
992 if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kLQ))) break;
993 if(!g->GetN()) break;
995 if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kNN))) break;
997 if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kESD))) break;
1004 case kEfficiencyKa:{
1005 gPad->Divide(2, 1, 1.e-5, 1.e-5);
1006 TList *l=gPad->GetListOfPrimitives();
1007 TVirtualPad *pad = ((TVirtualPad*)l->At(0));pad->cd();
1008 pad->SetMargin(0.1, 0.01, 0.1, 0.01);
1010 TLegend *legEff = new TLegend(.64, .84, .98, .98);
1011 legEff->SetBorderSize(1);legEff->SetTextSize(0.03255879);
1012 legEff->SetFillColor(0);
1013 h = (TH1S*)gROOT->FindObject("hEff");
1014 h=(TH1S*)h->Clone("hEff_K");
1015 h->SetYTitle("#epsilon_{K} [%]");
1016 h->GetYaxis()->SetRangeUser(1.e-2, 1.e2);
1018 content = (TList *)fGraph->FindObject(Form("Eff_%s", AliTRDCalPID::GetPartName(AliPID::kKaon)));
1019 if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kLQ))) break;
1020 if(!g->GetN()) break;
1021 legEff->SetHeader("PID Method [KAON]");
1022 g->Draw("pc"); legEff->AddEntry(g, "LQ 2D", "pl");
1023 if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kNN))) break;
1024 g->Draw("pc"); legEff->AddEntry(g, "NN", "pl");
1025 if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kESD))) break;
1026 g->Draw("p"); legEff->AddEntry(g, "ESD", "pl");
1033 TLegend *legEff2 = new TLegend(.64, .84, .98, .98);
1034 legEff2->SetBorderSize(1);legEff2->SetTextSize(0.03255879);
1035 legEff2->SetFillColor(0);
1036 pad = ((TVirtualPad*)l->At(1));pad->cd();
1037 pad->SetMargin(0.1, 0.01, 0.1, 0.01);
1038 h=(TH1S*)h->Clone("hEff_p");
1039 h->SetYTitle("#epsilon_{p} [%]");
1040 h->GetYaxis()->SetRangeUser(1.e-2, 1.e2);
1042 content = (TList *)fGraph->FindObject(Form("Eff_%s", AliTRDCalPID::GetPartName(AliPID::kProton)));
1043 if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kLQ))) break;
1044 if(!g->GetN()) break;
1045 legEff2->SetHeader("PID Method [PROTON]");
1046 g->Draw("pc"); legEff2->AddEntry(g, "LQ 2D", "pl");
1047 if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kNN))) break;
1048 g->Draw("pc"); legEff2->AddEntry(g, "NN", "pl");
1049 if(! (g = (TGraphErrors*)content->At(AliTRDpidUtil::kESD))) break;
1050 g->Draw("p"); legEff2->AddEntry(g, "ESD", "pl");
1059 // save 2.0 GeV projection as reference
1060 TLegend *legdEdx = new TLegend(.7, .7, .98, .98);
1061 legdEdx->SetBorderSize(1);
1063 if(!(h2 = (TH2F*)(fContainer->At(kdEdx)))) break;
1064 legdEdx->SetHeader("Particle Species");
1065 gPad->SetMargin(0.1, 0.01, 0.1, 0.01);
1066 for(Int_t is = AliPID::kSPECIES-1; is>=0; is--){
1067 Int_t bin = FindBin(is, 2.);
1068 h1 = h2->ProjectionY(Form("px%d", is), bin, bin);
1069 if(!h1->GetEntries()) continue;
1070 h1->Scale(1./h1->Integral());
1071 h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
1073 h1->GetXaxis()->SetTitle("dE/dx (a.u.)");
1074 h1->GetYaxis()->SetTitle("<Entries>");
1076 h = (TH1F*)h1->DrawClone(kFIRST ? "c" : "samec");
1077 legdEdx->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "l");
1091 /* gPad->Divide(2, 1, 1.e-5, 1.e-5);
1092 TList *l=gPad->GetListOfPrimitives();
1094 // save 2.0 GeV projection as reference
1095 TLegend *legPH = new TLegend(0.6865278,0.6882035,0.9655359,0.9688384,NULL,"brNDC");
1096 legPH->SetBorderSize(1);legPH->SetFillColor(0);
1097 legPH->SetHeader("Particle Species");
1098 if(!(arr = (TObjArray*)(fContainer->At(kPH)))) break;
1099 if(!(h2 = (TProfile2D*)(arr->At(0)))) break;
1101 /* TVirtualPad *pad = ((TVirtualPad*)l->At(0));pad->cd();
1102 pad->SetMargin(0.1, 0.01, 0.1, 0.01);*/
1104 for(Int_t is=0; is<AliPID::kSPECIES; is++){
1105 Int_t bin = FindBin(is, 2.);
1106 h1 = h2->ProjectionY(Form("pyt%d", is), bin, bin);
1107 if(!h1->GetEntries()) continue;
1108 h1->SetMarkerStyle(24);
1109 h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
1110 h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
1112 h1->GetXaxis()->SetTitle("t_{drift} [100*ns]");
1113 h1->GetYaxis()->SetTitle("<dQ/dt> @ 2GeV/c [a.u.]");
1114 h1->GetYaxis()->CenterTitle();h1->GetYaxis()->SetTitleOffset(1.2);
1116 h = (TH1F*)h1->DrawClone(kFIRST ? "c" : "samec");
1117 legPH->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "pl");
1121 // pad = ((TVirtualPad*)l->At(1));pad->cd();
1122 // pad->SetMargin(0.1, 0.01, 0.1, 0.01);
1123 // if(!(h2 = (TProfile2D*)(arr->At(1)))) break;
1125 // for(Int_t is=0; is<AliPID::kSPECIES; is++){
1126 // Int_t bin = FindBin(is, 2.);
1127 // h1 = h2->ProjectionY(Form("pyx%d", is), bin, bin);
1128 // if(!h1->GetEntries()) continue;
1129 // h1->SetMarkerStyle(24);
1130 // h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
1131 // h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
1133 // h1->GetXaxis()->SetTitle("x_{drift} [cm]");
1134 // h1->GetYaxis()->SetTitle("<dQ/dl> [a.u./cm]");
1136 // h1->DrawClone(kFIRST ? "c" : "samec");
1149 // save 2.0 GeV projection as reference
1150 // TLegend *legNClus = new TLegend(.13, .7, .4, .98);
1151 // legNClus->SetBorderSize(1);
1152 // legNClus->SetFillColor(0);
1155 if(!(h2 = (TH2F*)(fContainer->At(kNClus)))) break;
1156 // legNClus->SetHeader("Particle Species");
1157 for(Int_t is=0; is<AliPID::kSPECIES; is++){
1158 Int_t bin = FindBin(is, 2.);
1159 h1 = h2->ProjectionY(Form("pyNClus%d", is), bin, bin);
1160 if(!h1->GetEntries()) continue;
1161 h1->Scale(100./h1->Integral());
1162 //h1->SetMarkerStyle(24);
1163 //h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
1164 h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
1166 h1->GetXaxis()->SetTitle("N^{cl}/tracklet @ 2GeV/c");
1167 h1->GetYaxis()->SetTitle("Probability [%]");
1168 h1->GetYaxis()->CenterTitle();h1->GetYaxis()->SetTitleOffset(1.2);
1169 h = (TH1F*)h1->DrawClone("c");
1171 h->GetXaxis()->SetRangeUser(0., 35.);
1173 } else /*h = (TH1F*)*/h1->DrawClone("samec");
1175 // legNClus->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "l");
1178 // legNClus->Draw();
1189 /* TLegend *legNClus = new TLegend(.4, .7, .68, .98);
1190 legNClus->SetBorderSize(1);*/
1192 if(!(h2 = (TH2F*)(fContainer->At(kNTracklets)))) break;
1193 // legNClus->SetHeader("Particle Species");
1194 for(Int_t is=0; is<AliPID::kSPECIES; is++){
1195 Int_t bin = FindBin(is, 2.);
1196 h1 = h2->ProjectionY(Form("pyNTracklets%d", is), bin, bin);
1197 if(!h1->GetEntries()) continue;
1198 h1->Scale(100./h1->Integral());
1199 //h1->SetMarkerStyle(24);
1200 //h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
1201 h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
1203 h1->GetXaxis()->SetTitle("N^{trklt}/track @ 2GeV/c");
1204 h1->GetXaxis()->SetRangeUser(1.,6.);
1205 h1->GetYaxis()->SetTitle("Probability [%]");
1206 h1->GetYaxis()->CenterTitle();h1->GetYaxis()->SetTitleOffset(1.2);
1207 h1->GetYaxis()->SetRangeUser(0.,40.);
1209 (TH1F*)h1->DrawClone(kFIRST ? "c" : "samec");
1210 //legNClus->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "l");
1222 AliInfo(Form("Reference plot [%d] missing result", ifig));
1226 //________________________________________________________________________
1227 void AliTRDcheckPID::MakeSummary()
1229 // Put all representative pictures here
1230 if(!fGraph || !fContainer){
1231 AliError("Missing results");
1234 TVirtualPad *p(NULL); TCanvas *cOut(NULL);
1236 const Int_t /*nx(2048),*/ ny(1536);
1237 cOut = new TCanvas(GetName(), "TRD PID", ny, ny);
1238 cOut->Divide(2,2, 1.e-5, 1.e-5);
1239 p=cOut->cd(1); p->SetRightMargin(0.01882353);p->SetTopMargin(0.01785714);
1241 p=cOut->cd(2); p->SetRightMargin(0.01882353);p->SetTopMargin(0.01785714);
1242 GetRefFigure(kNClus);
1243 p=cOut->cd(3); p->SetRightMargin(0.01882353);p->SetTopMargin(0.01785714);
1245 p=cOut->cd(4); p->SetRightMargin(0.01882353);p->SetTopMargin(0.01785714);
1246 GetRefFigure(kNTracklets);
1247 cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1250 //________________________________________________________________________
1251 Bool_t AliTRDcheckPID::PostProcess()
1253 // Draw result to the screen
1254 // Called once at the end of the query
1257 Printf("ERROR: list not available");
1261 TObjArray *eff(NULL);
1263 fGraph = new TObjArray(2*AliPID::kSPECIES);
1266 if(!(eff = dynamic_cast<TObjArray *>(fContainer->At(kEfficiency)))){
1267 AliError("Efficiency container for Electrons missing.");
1270 EvaluateEfficiency(eff, fGraph, AliPID::kPion, 0.9);
1271 EvaluateEfficiency(eff, fGraph, AliPID::kKaon, 0.9);
1272 EvaluateEfficiency(eff, fGraph, AliPID::kProton, 0.9);
1278 //________________________________________________________________________
1279 void AliTRDcheckPID::EvaluateEfficiency(const TObjArray * const histoContainer, TObjArray *results, Int_t species, Float_t electronEfficiency){
1280 // Process PID information for pion efficiency
1282 fUtil->SetElectronEfficiency(electronEfficiency);
1285 const Int_t kNmethodsPID=Int_t(sizeof(fgMethod)/sizeof(Char_t*));
1286 Color_t colors[kNmethodsPID] = {kBlue, kGreen+2, kRed};
1287 Int_t markerStyle[kNmethodsPID] = {7, 7, 24};
1288 // efficiency graphs
1289 TGraphErrors *g(NULL);
1290 TObjArray *eff = new TObjArray(kNmethodsPID); eff->SetOwner(); eff->SetName(Form("Eff_%s", AliTRDCalPID::GetPartName(species)));
1291 results->AddAt(eff, species);
1292 for(Int_t iMethod = 0; iMethod < kNmethodsPID; iMethod++){
1293 eff->AddAt(g = new TGraphErrors(), iMethod);
1294 g->SetName(Form("%s", fgMethod[iMethod]));
1295 g->SetLineColor(colors[iMethod]);
1296 g->SetMarkerColor(colors[iMethod]);
1297 g->SetMarkerStyle(markerStyle[iMethod]);
1300 // Threshold graphs if not already
1301 TObjArray *thres(NULL);
1302 if(!(results->At(AliPID::kSPECIES))){
1303 thres = new TObjArray(kNmethodsPID); thres->SetOwner();
1304 thres->SetName("Thres");
1305 results->AddAt(thres, AliPID::kSPECIES);
1306 for(Int_t iMethod = 0; iMethod < kNmethodsPID; iMethod++){
1307 thres->AddAt(g = new TGraphErrors(), iMethod);
1308 g->SetName(Form("%s", fgMethod[iMethod]));
1309 g->SetLineColor(colors[iMethod]);
1310 g->SetMarkerColor(colors[iMethod]);
1311 g->SetMarkerStyle(markerStyle[iMethod]);
1315 TH2F *hPtr[kNmethodsPID]={
1316 (TH2F*)histoContainer->At(AliTRDpidUtil::kLQ),
1317 (TH2F*)histoContainer->At(AliTRDpidUtil::kNN),
1318 (TH2F*)histoContainer->At(AliTRDpidUtil::kESD)
1320 for(Int_t iMom = 0; iMom < fMomentumAxis->GetNbins(); iMom++){
1321 Float_t mom(fMomentumAxis->GetBinCenter(iMom+1));
1323 Int_t binEl(fMomentumAxis->GetNbins() * AliPID::kElectron + iMom + 1),
1324 binXX(fMomentumAxis->GetNbins() * species + iMom + 1);
1326 for(Int_t iMethod = 2; iMethod < kNmethodsPID; iMethod++){
1327 // Calculate the Species Efficiency at electronEfficiency% electron efficiency for each Method
1329 TH1D *histo1 = hPtr[iMethod] -> ProjectionY(Form("%s_el", fgMethod[iMethod]), binEl, binEl);
1330 TH1D *histo2 = hPtr[iMethod] -> ProjectionY(Form("%s_%s", fgMethod[iMethod], AliTRDCalPID::GetPartName(species)), binXX, binXX);
1332 if(!fUtil->CalculatePionEffi(histo1, histo2)) continue;
1334 g=(TGraphErrors*)eff->At(iMethod);
1335 g->SetPoint(iMom, mom, 1.e2*fUtil->GetPionEfficiency());
1336 g->SetPointError(iMom, 0., 1.e2*fUtil->GetError());
1337 AliDebug(2, Form("%s Efficiency for %s[p=%6.2fGeV/c] = %f +/- %f", fgMethod[iMethod], AliTRDCalPID::GetPartName(species), mom, fUtil->GetPionEfficiency(), fUtil->GetError()));
1339 if(!thres) continue;
1340 g=(TGraphErrors*)thres->At(iMethod);
1341 g->SetPoint(iMom, mom, fUtil->GetThreshold());
1342 g->SetPointError(iMom, 0., 0.);