9 #include "TProfile2D.h"
11 #include "TGraphErrors.h"
14 #include <TClonesArray.h>
15 #include <TObjArray.h>
18 #include "AliESDEvent.h"
19 #include "AliESDInputHandler.h"
20 #include "AliTrackReference.h"
22 #include "AliAnalysisTask.h"
24 #include "AliTRDtrackerV1.h"
25 #include "AliTRDtrackV1.h"
26 #include "AliTRDcluster.h"
27 #include "AliTRDReconstructor.h"
28 #include "AliCDBManager.h"
29 #include "AliTRDpidUtil.h"
31 #include "AliTRDpidChecker.h"
32 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
34 // calculate pion efficiency at 90% electron efficiency for 11 momentum bins
35 // this task should be used with simulated data only
37 ClassImp(AliTRDpidChecker)
39 //________________________________________________________________________
40 AliTRDpidChecker::AliTRDpidChecker()
41 :AliTRDrecoTask("PID", "PID Checker")
48 // Default constructor
51 fReconstructor = new AliTRDReconstructor();
52 fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
54 fUtil = new AliTRDpidUtil();
60 //________________________________________________________________________
61 AliTRDpidChecker::~AliTRDpidChecker()
63 if(fGraph){fGraph->Delete(); delete fGraph;}
64 if(fReconstructor) delete fReconstructor;
65 if(fUtil) delete fUtil;
69 //________________________________________________________________________
70 void AliTRDpidChecker::CreateOutputObjects()
75 OpenFile(0, "RECREATE");
76 fContainer = Histos();
80 //_______________________________________________________
81 TObjArray * AliTRDpidChecker::Histos(){
84 // Create QA histograms
86 if(fContainer) return fContainer;
88 Int_t xBins = AliPID::kSPECIES*AliTRDCalPID::kNMom;
89 fContainer = new TObjArray(); fContainer->Expand(7);
91 const Float_t epsilon = 1./(2*(AliTRDpidUtil::kBins-1)); // get nice histos with bin center at 0 and 1
93 // histos of the electron probability of all 5 particle species and 11 momenta for the 2-dim LQ method
94 fEfficiency = new TObjArray(); fEfficiency->Expand(3);
95 fContainer->AddAt(fEfficiency, kEfficiency);
98 if(!(h = (TH2F*)gROOT->FindObject("PID_LQ"))){
99 h = new TH2F("PID_LQ", "", xBins, -0.5, xBins - 0.5,
100 AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon);
102 fEfficiency->AddAt(h, kLQ);
104 // histos of the electron probability of all 5 particle species and 11 momenta for the neural network method
105 if(!(h = (TH2F*)gROOT->FindObject("PID_NN"))){
106 h = new TH2F("PID_NN", "",
107 xBins, -0.5, xBins - 0.5,
108 AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon);
110 fEfficiency->AddAt(h, kNN);
112 // histos of the electron probability of all 5 particle species and 11 momenta for the ESD output
113 if(!(h = (TH2F*)gROOT->FindObject("PID_ESD"))){
114 h = new TH2F("PID_ESD", "",
115 xBins, -0.5, xBins - 0.5,
116 AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon);
118 fEfficiency->AddAt(h, kESD);
120 // histos of the dE/dx distribution for all 5 particle species and 11 momenta
121 if(!(h = (TH2F*)gROOT->FindObject("dEdx"))){
122 h = new TH2F("dEdx", "",
123 xBins, -0.5, xBins - 0.5,
126 fContainer->AddAt(h, kdEdx);
128 // histos of the dE/dx slices for all 5 particle species and 11 momenta
129 if(!(h = (TH2F*)gROOT->FindObject("dEdxSlice"))){
130 h = new TH2F("dEdxSlice", "",
131 xBins*AliTRDReconstructor::kLQslices, -0.5, xBins*AliTRDReconstructor::kLQslices - 0.5,
134 fContainer->AddAt(h, kdEdxSlice);
136 // histos of the pulse height distribution for all 5 particle species and 11 momenta
137 if(!(h = (TH2F*)gROOT->FindObject("PH"))){
138 h = new TProfile2D("PH", "",
139 xBins, -0.5, xBins - 0.5,
140 AliTRDtrackerV1::GetNTimeBins(), -0.5, AliTRDtrackerV1::GetNTimeBins() - 0.5);
142 fContainer->AddAt(h, kPH);
144 // histos of the number of clusters distribution for all 5 particle species and 11 momenta
145 if(!(h = (TH2F*)gROOT->FindObject("NClus"))){
146 h = new TH2F("NClus", "",
147 xBins, -0.5, xBins - 0.5,
148 AliTRDtrackerV1::GetNTimeBins(), -0.5, AliTRDtrackerV1::GetNTimeBins() - 0.5);
150 fContainer->AddAt(h, kNClus);
153 // momentum distributions - absolute and in momentum bins
154 if(!(h = (TH1F*)gROOT->FindObject("hMom"))){
155 h = new TH1F("hMom", "momentum distribution", 100, 0., 12.);
157 fContainer->AddAt(h, kMomentum);
159 if(!(h = (TH1F*)gROOT->FindObject("hMomBin"))){
160 h = new TH1F("hMomBin", "momentum distribution in momentum bins", AliTRDCalPID::kNMom, 0.5, 11.5);
162 fContainer->AddAt(h, kMomentumBin);
169 //________________________________________________________________________
170 Bool_t AliTRDpidChecker::CheckTrackQuality(const AliTRDtrackV1* track)
173 // Check if the track is ok for PID
176 if(track->GetNumberOfTracklets() == AliTRDgeometry::kNlayer) return 1;
184 //________________________________________________________________________
185 Int_t AliTRDpidChecker::CalcPDG(AliTRDtrackV1* track)
188 track -> SetReconstructor(fReconstructor);
190 fReconstructor -> SetOption("nn");
193 if(track -> GetPID(AliPID::kElectron) > track -> GetPID(AliPID::kMuon) + track -> GetPID(AliPID::kPion) + track -> GetPID(AliPID::kKaon) + track -> GetPID(AliPID::kProton)){
196 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)){
199 else if(track -> GetPID(AliPID::kKaon) > track -> GetPID(AliPID::kMuon) && track -> GetPID(AliPID::kKaon) > track -> GetPID(AliPID::kPion)){
202 else if(track -> GetPID(AliPID::kMuon) > track -> GetPID(AliPID::kPion)){
211 //_______________________________________________________
212 TH1 *AliTRDpidChecker::PlotLQ(const AliTRDtrackV1 *track)
215 // Plot the probabilities for electrons using 2-dim LQ
218 if(track) fTrack = track;
220 AliWarning("No Track defined.");
224 if(!CheckTrackQuality(fTrack)) return 0x0;
226 if(!(fEfficiency = dynamic_cast<TObjArray *>(fContainer->At(kEfficiency)))){
227 AliWarning("No Histogram defined.");
231 if(!(hPIDLQ = dynamic_cast<TH2F *>(fEfficiency->At(kLQ)))){
232 AliWarning("No Histogram defined.");
236 AliTRDtrackV1 cTrack(*fTrack);
237 cTrack.SetReconstructor(fReconstructor);
240 Float_t momentum = 0.;
243 if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
246 //AliWarning("No MC info available!");
247 momentum = cTrack.GetMomentum(0);
248 pdg = CalcPDG(&cTrack);
250 if(momentum < 0.4) return 0x0;;
251 if(momentum > 12.) return 0x0;;
253 fReconstructor -> SetOption("!nn");
255 Int_t iMomBin = fUtil->GetMomentumBin(momentum);
260 hPIDLQ -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, cTrack.GetPID(AliPID::kElectron));
264 hPIDLQ -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, cTrack .GetPID(AliPID::kElectron));
268 hPIDLQ -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, cTrack .GetPID(AliPID::kElectron));
272 hPIDLQ -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, cTrack .GetPID(AliPID::kElectron));
276 hPIDLQ -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, cTrack.GetPID(AliPID::kElectron));
284 //_______________________________________________________
285 TH1 *AliTRDpidChecker::PlotNN(const AliTRDtrackV1 *track)
288 // Plot the probabilities for electrons using 2-dim LQ
291 if(track) fTrack = track;
293 AliWarning("No Track defined.");
297 if(!CheckTrackQuality(fTrack)) return 0x0;
299 if(!(fEfficiency = dynamic_cast<TObjArray *>(fContainer->At(kEfficiency)))){
300 AliWarning("No Histogram defined.");
304 if(!(hPIDNN = dynamic_cast<TH2F *>(fEfficiency->At(kNN)))){
305 AliWarning("No Histogram defined.");
310 AliTRDtrackV1 cTrack(*fTrack);
311 cTrack.SetReconstructor(fReconstructor);
314 Float_t momentum = 0.;
316 if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
319 //AliWarning("No MC info available!");
320 momentum = cTrack.GetMomentum(0);
321 pdg = CalcPDG(&cTrack);
323 if(momentum < 0.4) return 0x0;;
324 if(momentum > 12.) return 0x0;;
326 fReconstructor -> SetOption("nn");
328 Int_t iMomBin = fUtil -> GetMomentumBin(momentum);
333 hPIDNN -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, cTrack.GetPID(AliPID::kElectron));
337 hPIDNN -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, cTrack.GetPID(AliPID::kElectron));
341 hPIDNN -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, cTrack.GetPID(AliPID::kElectron));
345 hPIDNN -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, cTrack.GetPID(AliPID::kElectron));
349 hPIDNN -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, cTrack.GetPID(AliPID::kElectron));
356 //_______________________________________________________
357 TH1 *AliTRDpidChecker::PlotESD(const AliTRDtrackV1 *track)
360 // Plot the probabilities for electrons using 2-dim LQ
364 AliWarning("No ESD info available.");
368 if(track) fTrack = track;
370 AliWarning("No Track defined.");
374 if(!CheckTrackQuality(fTrack)) return 0x0;
376 if(!(fEfficiency = dynamic_cast<TObjArray *>(fContainer->At(kEfficiency)))){
377 AliWarning("No Histogram defined.");
381 if(!(hPIDESD = dynamic_cast<TH2F *>(fEfficiency->At(kESD)))){
382 AliWarning("No Histogram defined.");
388 Float_t momentum = 0.;
390 if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
393 //AliWarning("No MC info available!");
394 AliTRDtrackV1 cTrack(*fTrack);
395 cTrack.SetReconstructor(fReconstructor);
396 momentum = cTrack.GetMomentum(0);
397 pdg = CalcPDG(&cTrack);
399 if(momentum < 0.4) return 0x0;;
400 if(momentum > 12.) return 0x0;;
403 Int_t iMomBin = fUtil->GetMomentumBin(momentum);
406 // Double32_t pidESD[AliPID::kSPECIES];
407 const Double32_t *pidESD = fESD->GetResponseIter();
412 hPIDESD -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, pidESD[0]);
416 hPIDESD -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, pidESD[0]);
420 hPIDESD -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, pidESD[0]);
424 hPIDESD -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, pidESD[0]);
428 hPIDESD -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, pidESD[0]);
435 //_______________________________________________________
436 TH1 *AliTRDpidChecker::PlotdEdx(const AliTRDtrackV1 *track)
439 // Plot the probabilities for electrons using 2-dim LQ
442 if(track) fTrack = track;
444 AliWarning("No Track defined.");
448 if(!CheckTrackQuality(fTrack)) return 0x0;
451 if(!(hdEdx = dynamic_cast<TH2F *>(fContainer->At(kdEdx)))){
452 AliWarning("No Histogram defined.");
458 Float_t momentum = 0.;
460 if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
463 //AliWarning("No MC info available!");
464 AliTRDtrackV1 cTrack(*fTrack);
465 cTrack.SetReconstructor(fReconstructor);
466 momentum = cTrack.GetMomentum(0);
467 pdg = CalcPDG(&cTrack);
469 if(momentum < 0.4) return 0x0;;
470 if(momentum > 12.) return 0x0;;
472 Int_t iMomBin = fUtil -> GetMomentumBin(momentum);
476 Float_t SumdEdx[AliTRDgeometry::kNlayer];
477 AliTRDseedV1 *TRDtracklet[AliTRDgeometry::kNlayer];
478 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++) TRDtracklet[iChamb] = 0x0;
479 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++) TRDtracklet[iChamb] = fTrack -> GetTracklet(iChamb);
482 Float_t dEdxSlice[AliTRDgeometry::kNlayer][AliTRDReconstructor::kLQslices];
484 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
485 SumdEdx[iChamb] = 0.;
486 fdEdx = TRDtracklet[iChamb] -> GetdEdx();
487 SumdEdx[iChamb] += fdEdx[0] + fdEdx[1] + fdEdx[2];
488 for(Int_t iSlice = 0; iSlice < AliTRDReconstructor::kLQslices; iSlice++){
489 dEdxSlice[iChamb][iSlice] = fdEdx[iSlice];
496 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++)
497 hdEdx -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
501 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++)
502 hdEdx -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
506 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++)
507 hdEdx -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
511 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++)
512 hdEdx -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
516 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++)
517 hdEdx -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
525 //_______________________________________________________
526 TH1 *AliTRDpidChecker::PlotdEdxSlice(const AliTRDtrackV1 *track)
529 // Plot the probabilities for electrons using 2-dim LQ
532 if(track) fTrack = track;
534 AliWarning("No Track defined.");
538 if(!CheckTrackQuality(fTrack)) return 0x0;
541 if(!(hdEdxSlice = dynamic_cast<TH2F *>(fContainer->At(kdEdxSlice)))){
542 AliWarning("No Histogram defined.");
547 Float_t momentum = 0.;
549 if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
552 //AliWarning("No MC info available!");
553 AliTRDtrackV1 cTrack(*fTrack);
554 cTrack.SetReconstructor(fReconstructor);
555 momentum = cTrack.GetMomentum(0);
556 pdg = CalcPDG(&cTrack);
558 if(momentum < 0.4) return 0x0;
559 if(momentum > 12.) return 0x0;;
561 Int_t iMomBin = fUtil -> GetMomentumBin(momentum);
565 AliTRDseedV1 *TRDtracklet[AliTRDgeometry::kNlayer];
566 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++) TRDtracklet[iChamb] = 0x0;
567 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++) TRDtracklet[iChamb] = fTrack -> GetTracklet(iChamb);
570 Float_t dEdxSlice[AliTRDgeometry::kNlayer][AliTRDReconstructor::kLQslices];
572 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
573 fdEdx = TRDtracklet[iChamb] -> GetdEdx();
574 for(Int_t iSlice = 0; iSlice < AliTRDReconstructor::kLQslices; iSlice++){
575 dEdxSlice[iChamb][iSlice] = fdEdx[iSlice];
582 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
583 for(Int_t iSlice = 0; iSlice < AliTRDReconstructor::kLQslices; iSlice++){
584 hdEdxSlice -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom * AliTRDReconstructor::kLQslices + iMomBin * AliTRDReconstructor::kLQslices + iSlice, dEdxSlice[iChamb][iSlice]);
590 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
591 for(Int_t iSlice = 0; iSlice < AliTRDReconstructor::kLQslices; iSlice++){
592 hdEdxSlice -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom * AliTRDReconstructor::kLQslices + iMomBin * AliTRDReconstructor::kLQslices + iSlice,
593 dEdxSlice[iChamb][iSlice]);
599 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++)
600 for(Int_t iSlice = 0; iSlice < AliTRDReconstructor::kLQslices; iSlice++)
601 hdEdxSlice -> Fill(AliPID::kPion * AliTRDCalPID::kNMom * AliTRDReconstructor::kLQslices + iMomBin * AliTRDReconstructor::kLQslices + iSlice,
602 dEdxSlice[iChamb][iSlice]);
606 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++)
607 for(Int_t iSlice = 0; iSlice < AliTRDReconstructor::kLQslices; iSlice++)
608 hdEdxSlice -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom * AliTRDReconstructor::kLQslices + iMomBin * AliTRDReconstructor::kLQslices + iSlice,
609 dEdxSlice[iChamb][iSlice]);
613 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++)
614 for(Int_t iSlice = 0; iSlice < AliTRDReconstructor::kLQslices; iSlice++)
615 hdEdxSlice -> Fill(AliPID::kProton * AliTRDCalPID::kNMom * AliTRDReconstructor::kLQslices + iMomBin * AliTRDReconstructor::kLQslices + iSlice,
616 dEdxSlice[iChamb][iSlice]);
625 //_______________________________________________________
626 TH1 *AliTRDpidChecker::PlotPH(const AliTRDtrackV1 *track)
629 // Plot the probabilities for electrons using 2-dim LQ
632 if(track) fTrack = track;
634 AliWarning("No Track defined.");
638 if(!CheckTrackQuality(fTrack)) return 0x0;
641 if(!(hPH = dynamic_cast<TProfile2D *>(fContainer->At(kPH)))){
642 AliWarning("No Histogram defined.");
647 Float_t momentum = 0.;
649 if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
652 //AliWarning("No MC info available!");
653 AliTRDtrackV1 cTrack(*fTrack);
654 cTrack.SetReconstructor(fReconstructor);
655 momentum = cTrack.GetMomentum(0);
656 pdg = CalcPDG(&cTrack);
658 if(momentum < 0.4) return 0x0;;
659 if(momentum > 12.) return 0x0;;
661 Int_t iMomBin = fUtil -> GetMomentumBin(momentum);
663 AliTRDseedV1 *TRDtracklet[AliTRDgeometry::kNlayer];
664 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++) TRDtracklet[iChamb] = 0x0;
665 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++) TRDtracklet[iChamb] = fTrack -> GetTracklet(iChamb);
667 AliTRDcluster *TRDcluster = 0x0;
672 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
673 for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
674 if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
676 hPH -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
682 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
683 for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
684 if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
686 hPH -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
692 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
693 for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
694 if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
696 hPH -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
702 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
703 for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
704 if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
706 hPH -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
712 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
713 for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
714 if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
716 hPH -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
726 //_______________________________________________________
727 TH1 *AliTRDpidChecker::PlotNClus(const AliTRDtrackV1 *track)
730 // Plot the probabilities for electrons using 2-dim LQ
733 if(track) fTrack = track;
735 AliWarning("No Track defined.");
739 if(!CheckTrackQuality(fTrack)) return 0x0;
742 if(!(hNClus = dynamic_cast<TH2F *>(fContainer->At(kNClus)))){
743 AliWarning("No Histogram defined.");
749 Float_t momentum = 0.;
751 if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
754 //AliWarning("No MC info available!");
755 AliTRDtrackV1 cTrack(*fTrack);
756 cTrack.SetReconstructor(fReconstructor);
757 momentum = cTrack.GetMomentum(0);
758 pdg = CalcPDG(&cTrack);
760 if(momentum < 0.4) return 0x0;;
761 if(momentum > 12.) return 0x0;;
763 Int_t iMomBin = fUtil -> GetMomentumBin(momentum);
766 Int_t iNClus[AliTRDgeometry::kNlayer];
767 memset(iNClus, 0, sizeof(Int_t) * AliTRDgeometry::kNlayer);
769 AliTRDseedV1 *TRDtracklet[AliTRDgeometry::kNlayer];
770 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++) TRDtracklet[iChamb] = 0x0;
771 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
772 TRDtracklet[iChamb] = fTrack -> GetTracklet(iChamb);
773 iNClus[iChamb] = TRDtracklet[iChamb] -> GetN();
779 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++)
780 hNClus -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, iNClus[iChamb]);
784 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++)
785 hNClus -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, iNClus[iChamb]);
789 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++)
790 hNClus -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, iNClus[iChamb]);
794 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++)
795 hNClus -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, iNClus[iChamb]);
799 for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++)
800 hNClus -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, iNClus[iChamb]);
808 //_______________________________________________________
809 TH1 *AliTRDpidChecker::PlotMom(const AliTRDtrackV1 *track)
812 // Plot the probabilities for electrons using 2-dim LQ
815 if(track) fTrack = track;
817 AliWarning("No Track defined.");
821 if(!CheckTrackQuality(fTrack)) return 0x0;
824 if(!(hMom = dynamic_cast<TH1F *>(fContainer->At(kMomentum)))){
825 AliWarning("No Histogram defined.");
831 Float_t momentum = 0.;
833 if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
836 //AliWarning("No MC info available!");
837 AliTRDtrackV1 cTrack(*fTrack);
838 cTrack.SetReconstructor(fReconstructor);
839 momentum = cTrack.GetMomentum(0);
840 pdg = CalcPDG(&cTrack);
842 if(momentum < 0.4) return 0x0;
843 if(momentum > 12.) return 0x0;;
845 hMom -> Fill(momentum);
850 //_______________________________________________________
851 TH1 *AliTRDpidChecker::PlotMomBin(const AliTRDtrackV1 *track)
854 // Plot the probabilities for electrons using 2-dim LQ
857 if(track) fTrack = track;
859 AliWarning("No Track defined.");
863 if(!CheckTrackQuality(fTrack)) return 0x0;
866 if(!(hMomBin = dynamic_cast<TH1F *>(fContainer->At(kMomentumBin)))){
867 AliWarning("No Histogram defined.");
873 Float_t momentum = 0.;
876 if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
879 //AliWarning("No MC info available!");
880 AliTRDtrackV1 cTrack(*fTrack);
881 cTrack.SetReconstructor(fReconstructor);
882 momentum = cTrack.GetMomentum(0);
883 pdg = CalcPDG(&cTrack);
885 if(momentum < 0.4) return 0x0;
886 if(momentum > 12.) return 0x0;;
888 Int_t iMomBin = fUtil -> GetMomentumBin(momentum);
889 hMomBin -> Fill(iMomBin);
894 //________________________________________________________
895 Bool_t AliTRDpidChecker::GetRefFigure(Int_t ifig)
897 Bool_t FIRST = kTRUE;
898 TLegend *leg = new TLegend(.7, .7, .98, .98);
899 leg->SetBorderSize(1);
900 TGraphErrors *g = 0x0;
902 TH1 *h1 = 0x0, *h=0x0;
906 if(!(g = (TGraphErrors*)fGraph->At(kLQ))) break;
907 if(!g->GetN()) break;
908 leg->SetHeader("PID Method");
910 ax = g->GetHistogram()->GetXaxis();
911 ax->SetTitle("p [GeV/c]");
912 ax->SetRangeUser(.6, 10.5);
913 ax->SetMoreLogLabels();
914 ax = g->GetHistogram()->GetYaxis();
915 ax->SetTitle("#epsilon_{#pi} [%]");
916 ax->SetRangeUser(1.e-3, 1.e-1);
917 leg->AddEntry(g, "2D LQ", "pl");
918 if(! (g = (TGraphErrors*)fGraph->At(kNN))) break;
920 leg->AddEntry(g, "NN", "pl");
921 if(! (g = (TGraphErrors*)fGraph->At(kESD))) break;
923 leg->AddEntry(g, "ESD", "pl");
931 // save 2.0 GeV projection as reference
933 if(!(h2 = (TH2F*)(fContainer->At(kdEdx)))) break;
934 leg->SetHeader("Particle Species");
935 for(Int_t is = AliPID::kSPECIES-1; is>=0; is--){
936 Int_t bin = is*AliTRDCalPID::kNMom+4;
937 h1 = h2->ProjectionY("px", bin, bin);
938 if(!h1->GetEntries()) continue;
939 h1->Scale(1./h1->Integral());
940 h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
941 h = (TH1F*)h1->DrawClone(FIRST ? "c" : "samec");
942 leg->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "l");
955 // save 2.0 GeV projection as reference
957 if(!(h2 = (TH2F*)(fContainer->At(kPH)))) break;;
958 leg->SetHeader("Particle Species");
959 for(Int_t is=0; is<AliPID::kSPECIES; is++){
960 Int_t bin = is*AliTRDCalPID::kNMom+4;
961 h1 = h2->ProjectionY("py", bin, bin);
962 if(!h1->GetEntries()) continue;
963 h1->SetMarkerStyle(24);
964 h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
965 h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
967 h1->GetXaxis()->SetTitle("tb[1/100 ns^{-1}]");
968 h1->GetYaxis()->SetTitle("<PH> [a.u.]");
970 h = (TH1F*)h1->DrawClone(FIRST ? "c" : "samec");
971 leg->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "pl");
982 // save 2.0 GeV projection as reference
984 if(!(h2 = (TH2F*)(fContainer->At(kNClus)))) break;
985 leg->SetHeader("Particle Species");
986 for(Int_t is=0; is<AliPID::kSPECIES; is++){
987 Int_t bin = is*AliTRDCalPID::kNMom+4;
988 h1 = h2->ProjectionY("py", bin, bin);
989 if(!h1->GetEntries()) continue;
990 //h1->SetMarkerStyle(24);
991 //h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
992 h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
993 if(FIRST) h1->GetXaxis()->SetTitle("N^{cl}/tracklet");
994 h = (TH1F*)h1->DrawClone(FIRST ? "c" : "samec");
995 leg->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "l");
1009 if(!(g = (TGraphErrors*)fGraph->At(kLQ+3))) break;
1010 if(!g->GetN()) break;
1011 leg->SetHeader("PID Method");
1013 ax = g->GetHistogram()->GetXaxis();
1014 ax->SetTitle("p [GeV/c]");
1015 ax->SetRangeUser(.6, 10.5);
1016 ax->SetMoreLogLabels();
1017 ax = g->GetHistogram()->GetYaxis();
1018 ax->SetTitle("threshold");
1019 ax->SetRangeUser(5.e-2, 1.);
1020 leg->AddEntry(g, "2D LQ", "pl");
1021 if(!(g = (TGraphErrors*)fGraph->At(kNN+3))) break;
1023 leg->AddEntry(g, "NN", "pl");
1024 if(!(g = (TGraphErrors*)fGraph->At(kESD+3))) break;
1026 leg->AddEntry(g, "ESD", "pl");
1033 AliInfo(Form("Reference plot [%d] missing result", ifig));
1037 //________________________________________________________________________
1038 Bool_t AliTRDpidChecker::PostProcess()
1040 // Draw result to the screen
1041 // Called once at the end of the query
1044 Printf("ERROR: list not available");
1047 if(!(fEfficiency = dynamic_cast<TObjArray *>(fContainer->At(kEfficiency)))){
1048 AliError("Efficiency container missing.");
1052 TGraphErrors *g = 0x0;
1054 fGraph = new TObjArray(6);
1057 fGraph->AddAt(g = new TGraphErrors(), kLQ);
1058 g->SetLineColor(kBlue);
1059 g->SetMarkerColor(kBlue);
1060 g->SetMarkerStyle(7);
1062 fGraph->AddAt(g = new TGraphErrors(), kNN);
1063 g->SetLineColor(kGreen);
1064 g->SetMarkerColor(kGreen);
1065 g->SetMarkerStyle(7);
1067 fGraph -> AddAt(g = new TGraphErrors(), kESD);
1068 g->SetLineColor(kRed);
1069 g->SetMarkerColor(kRed);
1070 g->SetMarkerStyle(24);
1072 fGraph->AddAt(g = new TGraphErrors(), 3+kLQ);
1073 g->SetLineColor(kBlue);
1074 g->SetMarkerColor(kBlue);
1075 g->SetMarkerStyle(7);
1077 fGraph->AddAt(g = new TGraphErrors(), 3+kNN);
1078 g->SetLineColor(kGreen);
1079 g->SetMarkerColor(kGreen);
1080 g->SetMarkerStyle(7);
1082 fGraph -> AddAt(g = new TGraphErrors(), 3+kESD);
1083 g->SetLineColor(kRed);
1084 g->SetMarkerColor(kRed);
1085 g->SetMarkerStyle(24);
1089 TH1D *Histo1=0x0, *Histo2=0x0;
1091 // calculate the pion efficiencies and the errors for 90% electron efficiency (2-dim LQ)
1092 TH2F *hPIDLQ = (TH2F*)fEfficiency->At(kLQ);
1093 for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
1094 mom = AliTRDCalPID::GetMomentum(iMom);
1096 Histo1 = hPIDLQ -> ProjectionY("LQ_ele",AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1,AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1);
1097 Histo2 = hPIDLQ -> ProjectionY("LQ_pio",AliTRDCalPID::kNMom*AliPID::kPion+iMom+1,AliTRDCalPID::kNMom*AliPID::kPion+iMom+1);
1099 if(!fUtil->CalculatePionEffi(Histo1, Histo2)) continue;
1101 g = (TGraphErrors*)fGraph->At(kLQ);
1102 g->SetPoint(iMom, mom, fUtil->GetPionEfficiency());
1103 g->SetPointError(iMom, 0., fUtil->GetError());
1104 g = (TGraphErrors*)fGraph->At(3 + kLQ);
1105 g->SetPoint(iMom, mom, fUtil->GetThreshold());
1106 g->SetPointError(iMom, 0., 0.);
1108 if(fDebugLevel>=2) Printf("Pion Efficiency for 2-dim LQ is : %f +/- %f\n\n", fUtil->GetPionEfficiency(), fUtil->GetError());
1112 // calculate the pion efficiencies and the errors for 90% electron efficiency (NN)
1113 TH2F *hPIDNN = (TH2F*)fEfficiency->At(kNN);
1114 for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
1115 mom = AliTRDCalPID::GetMomentum(iMom);
1117 Histo1 = hPIDNN -> ProjectionY("NN_ele",AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1,AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1);
1118 Histo2 = hPIDNN -> ProjectionY("NN_pio",AliTRDCalPID::kNMom*AliPID::kPion+iMom+1,AliTRDCalPID::kNMom*AliPID::kPion+iMom+1);
1120 if(!fUtil -> CalculatePionEffi(Histo1, Histo2)) continue;
1122 g = (TGraphErrors*)fGraph->At(kNN);
1123 g->SetPoint(iMom, mom, fUtil->GetPionEfficiency());
1124 g->SetPointError(iMom, 0., fUtil->GetError());
1125 g = (TGraphErrors*)fGraph->At(3+kNN);
1126 g->SetPoint(iMom, mom, fUtil->GetThreshold());
1127 g->SetPointError(iMom, 0., 0.);
1129 if(fDebugLevel>=2) Printf("Pion Efficiency for NN is : %f +/- %f\n\n", fUtil->GetPionEfficiency(), fUtil->GetError());
1133 // calculate the pion efficiencies and the errors for 90% electron efficiency (ESD)
1134 TH2F *hPIDESD = (TH2F*)fEfficiency->At(kESD);
1135 for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
1136 mom = AliTRDCalPID::GetMomentum(iMom);
1138 Histo1 = hPIDESD -> ProjectionY("NN_ele",AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1,AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1);
1139 Histo2 = hPIDESD -> ProjectionY("NN_pio",AliTRDCalPID::kNMom*AliPID::kPion+iMom+1,AliTRDCalPID::kNMom*AliPID::kPion+iMom+1);
1141 if(!fUtil->CalculatePionEffi(Histo1, Histo2)) continue;
1143 g = (TGraphErrors*)fGraph->At(kESD);
1144 g->SetPoint(iMom, mom, fUtil->GetPionEfficiency());
1145 g->SetPointError(iMom, 0., fUtil->GetError());
1146 g = (TGraphErrors*)fGraph->At(3+kESD);
1147 g->SetPoint(iMom, mom, fUtil->GetThreshold());
1148 g->SetPointError(iMom, 0., 0.);
1150 if(fDebugLevel>=2) Printf("Pion Efficiency for ESD is : %f +/- %f\n\n", fUtil->GetPionEfficiency(), fUtil->GetError());
1159 //________________________________________________________________________
1160 void AliTRDpidChecker::Terminate(Option_t *)
1162 // Draw result to the screen
1163 // Called once at the end of the query
1165 fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
1167 Printf("ERROR: list not available");