1 //-------------------------------------------------------------------------
3 // This is the PROOF-enabled version of TRD/Macros/AliTRDComparisonV2.C macro.
4 // Origin: Andrei.Zalite@cern.ch
6 //-------------------------------------------------------------------------
13 #include "TClonesArray.h"
21 #include "AliVEvent.h"
22 #include "AliESDEvent.h"
23 #include "AliMCEvent.h"
24 #include "AliESDtrack.h"
25 #include "AliTrackReference.h"
26 #include "AliMCComparisonTrack.h"
28 #include "AliTRDComparisonTask.h"
30 ClassImp(AliTRDComparisonTask)
32 extern TStyle *gStyle;
34 AliTRDComparisonTask::AliTRDComparisonTask()
35 : AliAnalysisTaskSE("AliTRDComparisonTask"),
51 // Default constructor
52 AliInfo("Default constructor AliTRDComparisonTask");
53 // Define input and output slots here
54 // Input slot #0 works with a TChain
55 DefineInput(0, TChain::Class());
56 // Output slot #1 TList
57 DefineOutput(1, TList::Class());
60 AliTRDComparisonTask::AliTRDComparisonTask(const char* name)
61 : AliAnalysisTaskSE(name),
78 AliInfo("Constructor AliTRDComparisonTask");
79 // Define input and output slots here
80 // Input slot #0 works with a TChain
81 DefineInput(0, TChain::Class());
82 // Output slot #1 TList
83 DefineOutput(1, TList::Class());
88 void AliTRDComparisonTask::UserCreateOutputObjects()
92 AliInfo("AliTRDComparisonTask::UserCreateOutputObjects");
93 // Create output container
94 fListOfHistos = new TList();
96 fGood = new TH1F("fGood", "Pt for good tracks", 34, 0.2, 7.0);
97 fFound = new TH1F("fFound", "Pt for found tracks", 34, 0.2, 7.0);
98 fFake = new TH1F("fFake", "Pt for fake tracks", 34, 0.2, 7.0);
99 fP = new TH1F("fP", "PHI resolution", 50, -20., 20.);
100 fL = new TH1F("fL", "LAMBDA resolution", 50, -20., 20.);
101 fPt = new TH1F("fPt", "Relative Pt resolution", 30, -10., 10.);
102 fHmpt = new TH1F("fHmpt", "Y and Z resolution", 30, -30., 30.);
103 fE = new TH1F("fE", "dE/dX for pions with 0.4<p<0.5 GeV/c", 50, 0., 1000.);
104 fEp = new TH2F("fEp", "dE/dX vs momentum", 50, 0., 2., 50, 0., 2000.);
105 fGoodPhi = new TH1F("fGoodPhi", "Phi for Good tracks", 90, 0., 2.*TMath::Pi());
106 fFoundPhi = new TH1F("fFoundPhi", "Phi for Found tracks", 90, 0., 2.*TMath::Pi());
107 fZ = new TH1F("fZ", "Z resolution", 30, -30., 30.);
108 fC = new TH1F("fC", "Number of the assigned clusters", 25, 110., 135.);
110 fListOfHistos->Add(fGood);
111 fListOfHistos->Add(fFound);
112 fListOfHistos->Add(fFake);
113 fListOfHistos->Add(fP);
114 fListOfHistos->Add(fL);
115 fListOfHistos->Add(fPt);
116 fListOfHistos->Add(fHmpt);
117 fListOfHistos->Add(fE);
118 fListOfHistos->Add(fEp);
119 fListOfHistos->Add(fGoodPhi);
120 fListOfHistos->Add(fFoundPhi);
121 fListOfHistos->Add(fZ);
122 fListOfHistos->Add(fC);
126 void AliTRDComparisonTask::UserExec(Option_t *)
130 // Called for each event
133 AliMCEvent* mcEvent = MCEvent();
135 Printf("ERROR: Could not retrieve MC event");
141 TClonesArray dummy("AliMCComparisonTrack",1000), *refs=&dummy;
143 // Loop over all MC tracks and select "good" tracks
144 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
145 AliMCParticle* track = (AliMCParticle*)mcEvent->GetTrack(iTracks);
147 Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
152 if (track->Pt() < 0.2) continue;
153 if (track->Pt() > 7.0) continue;
154 if (TMath::Abs(track->Pz()/track->Pt()) > 0.999) continue;
156 Double_t vx = track->Xv(), vy = track->Yv(), vz = track->Zv();
157 if (TMath::Sqrt(vx*vx+vy*vy) > 3.5) continue;
158 if (TMath::Abs(vz) > 50.) continue;
160 // Check TRD information
161 Int_t nRefs = track->GetNumberOfTrackReferences();
162 if (nRefs <= 1) continue;
164 AliTrackReference* trackRef0 = 0;
165 for (Int_t iTrackRef = 0; iTrackRef < nRefs; iTrackRef++) {
166 trackRef0 = track->GetTrackReference(iTrackRef);
168 Int_t detectorId = trackRef0->DetectorId();
169 if (detectorId == AliTrackReference::kTRD) {
175 if (!trackRef0) continue;
177 if (trackRef0->LocalX() > 300.) continue;
179 AliTrackReference* trackRef = 0;
180 for (Int_t iTrackRef = nRefs - 1; iTrackRef >= 0; --iTrackRef) {
181 trackRef = track->GetTrackReference(iTrackRef);
182 if (!trackRef) continue;
183 if (trackRef->LocalX() > 363. &&
184 trackRef->DetectorId() == AliTrackReference::kTRD) break;
187 if (!trackRef) continue;
189 if (TMath::Abs(trackRef0->Alpha() - trackRef->Alpha()) > 1e-5) continue;
191 // Check TPC information
192 Bool_t labelTPC = kFALSE;
193 for (Int_t iTrackRef = 0; iTrackRef < nRefs; iTrackRef++) {
194 trackRef = track->GetTrackReference(iTrackRef);
196 Int_t detectorId = trackRef->DetectorId();
197 if (detectorId == AliTrackReference::kTPC) {
202 } // track references loop
206 AliMCComparisonTrack* ref = new((*refs)[nt]) AliMCComparisonTrack();
207 ref->SetLabel(iTracks);
208 TParticle* particle = track->Particle();
209 Int_t pdg = particle->GetPdgCode();
211 ref->SetPz(trackRef->Pz());
212 Float_t pt = trackRef->Pt();
215 Float_t phig = trackRef->Phi();
217 fGoodPhi->Fill(phig);
218 ref->SetLocalX(trackRef->LocalX());
219 ref->SetLocalY(trackRef->LocalY());
220 ref->SetZ(trackRef->Z());
227 AliVEvent* event = InputEvent();
229 Printf("ERROR: Could not retrieve event");
233 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
240 Int_t mcGoods = refs->GetEntriesFast();
242 // Loop over all "good" MC tracks
243 for (Int_t k = 0; k < mcGoods; k++) {
244 AliMCComparisonTrack* ref = (AliMCComparisonTrack*)refs->UncheckedAt(k);
246 Int_t mcLabel = ref->GetLabel();
247 Float_t ptg = ref->GetPt();
248 Float_t phiG = ref->GetPhi();
250 for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
251 AliESDtrack* track = esd->GetTrack(iTrack);
253 Printf("ERROR: Could not receive track %d", iTrack);
257 if (! track->IsOn(AliESDtrack::kTRDout)) continue;
258 const AliExternalTrackParam * outorig = track->GetOuterParam();
259 if (!outorig) continue;
261 Int_t lable = track->GetTRDLabel();
263 if (mcLabel == TMath::Abs(lable)) {
264 if (mcLabel == lable) {
267 fFoundPhi->Fill(phiG);
275 AliExternalTrackParam out(*outorig);
277 Double_t bz = esd->GetMagneticField();
278 Double_t xg = ref->GetLocalX();
279 out.PropagateTo(xg,bz);
281 Float_t phi = TMath::ASin(out.GetSnp()) + out.GetAlpha();
282 if (phi < 0) phi += 2*TMath::Pi();
283 if (phi >= 2*TMath::Pi()) phi -= 2*TMath::Pi();
284 Float_t phig = ref->GetPhi();
285 fP->Fill((phi - phig)*1000.);
287 Float_t lam = TMath::ATan(out.GetTgl());
288 Float_t lamg = TMath::ATan2(ref->GetPz(), ptg);
289 fL->Fill((lam - lamg)*1000.);
291 fC->Fill(track->GetTRDclusters(0));
293 Float_t pt1 = out.OneOverPt();
294 fPt->Fill((pt1 - 1/ptg)/(1/ptg)*100.);
296 Float_t y = out.GetY();
297 Float_t yg = ref->GetLocalY();
298 fHmpt->Fill(10*(y - yg));
300 Float_t z = out.GetZ();
301 Float_t zg = ref->GetZ();
302 fZ->Fill(10.*(z - zg));
304 Float_t mom = 1./(pt1*TMath::Cos(lam));
305 Float_t dedx = track->GetTRDsignal();
306 Printf (" dedx = %f ", dedx);
307 fEp->Fill(mom, dedx, 1.);
309 Int_t pdg = ref->GetPDG();
310 if (TMath::Abs(pdg)==211) //pions
311 if (mom>0.4 && mom<0.5) fE->Fill(dedx,1.);
321 Printf(" Results: " );
322 Printf(" Good %d Found %d Fake %d Lost %d ", nt, nfound, nfake, nlost);
327 PostData(1, fListOfHistos);
331 void AliTRDComparisonTask::Terminate(Option_t *) {
332 // Draw result to the screen
333 // Called once at the end of the query
334 fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
335 if (!fListOfHistos) {
336 Printf("ERROR: fListOfHistos not available");
340 fGood = dynamic_cast<TH1F*>(fListOfHistos->At(0));
341 fFound = dynamic_cast<TH1F*>(fListOfHistos->At(1));
342 fFake = dynamic_cast<TH1F*>(fListOfHistos->At(2));
343 fP = dynamic_cast<TH1F*>(fListOfHistos->At(3));
344 fL = dynamic_cast<TH1F*>(fListOfHistos->At(4));
345 fPt = dynamic_cast<TH1F*>(fListOfHistos->At(5));
346 fHmpt = dynamic_cast<TH1F*>(fListOfHistos->At(6));
347 fE = dynamic_cast<TH1F*>(fListOfHistos->At(7));
348 fEp = dynamic_cast<TH2F*>(fListOfHistos->At(8));
349 fGoodPhi = dynamic_cast<TH1F*>(fListOfHistos->At(9));
350 fFoundPhi = dynamic_cast<TH1F*>(fListOfHistos->At(10));
351 fZ = dynamic_cast<TH1F*>(fListOfHistos->At(11));
352 fC = dynamic_cast<TH1F*>(fListOfHistos->At(12));
354 gStyle->SetOptStat(111110);
355 gStyle->SetOptFit(1);
357 TCanvas* c1 = new TCanvas("c1", "", 0, 0, 700, 850);
361 TPad* p1 = new TPad("p1", "", 0., 0.3, 0.5, 0.6); p1->Draw();
362 p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
363 fP->SetFillColor(4); fP->SetXTitle("(mrad)");
364 if (fP->GetEntries() < minc) fP->Draw(); else fP->Fit("gaus"); c1->cd();
366 TPad* p2 = new TPad("p2", "", 0.5, 0.3, 1., 0.6); p2->Draw();
367 p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10);
368 fL->SetFillColor(4); fL->SetXTitle("(mrad)");
369 if (fL->GetEntries() < minc) fL->Draw(); else fL->Fit("gaus"); c1->cd();
371 TPad* p3 = new TPad("p3", "", 0., 0., 0.5, 0.3); p3->Draw();
372 p3->cd(); p3->SetFillColor(42); p3->SetFrameFillColor(10);
373 fPt->SetFillColor(2); fPt->SetXTitle("(%)");
374 if (fPt->GetEntries() < minc) fPt->Draw(); else fPt->Fit("gaus"); c1->cd();
376 TPad* p4 = new TPad("p4", "", 0.5, 0., 1., 0.3); p4->Draw();
377 p4->cd(); p4->SetFillColor(42); p4->SetFrameFillColor(10);
378 fHmpt->SetFillColor(6); fHmpt->SetXTitle("(mm)");
379 if (fHmpt->GetEntries() < minc) fHmpt->Draw(); else fHmpt->Fit("gaus");
380 fZ->Draw("same"); c1->cd();
382 TPad* p5 = new TPad("p5", "", 0., 0.6, 1., 1.); p5->Draw(); p5->cd();
383 p5->SetFillColor(41); p5->SetFrameFillColor(10);
384 fFound->Sumw2(); fGood->Sumw2(); fFake->Sumw2();
385 TH1F* hg = new TH1F("hg", "Efficiency for good tracks", 34, 0.2, 7.0);
386 TH1F* hf = new TH1F("hf", "Efficiency for fake tracks", 34, 0.2, 7.0);
387 hg->Divide(fFound,fGood,1.,1.,"B");
388 hf->Divide(fFake,fGood,1.,1.,"B");
390 hg->SetYTitle("Tracking efficiency");
391 hg->SetXTitle("Pt (GeV/c)");
394 TLine* line1 = new TLine(0.2, 1.0, 7.0, 1.0); line1->SetLineStyle(4);
396 TLine* line2 = new TLine(0.2, 0.9, 7.0, 0.9); line2->SetLineStyle(4);
400 hf->SetFillStyle(3013);
403 hf->Draw("histsame");
404 TText* text = new TText(0.461176, 0.248448, "Fake tracks");
405 text->SetTextSize(0.05);
407 text = new TText(0.453919, 1.11408, "Good tracks");
408 text->SetTextSize(0.05);
411 TCanvas* c2 = new TCanvas("c2", "", 320, 32, 530, 590);
413 TPad* p6 = new TPad("p6", "", 0., 0., 1., .5); p6->Draw();
414 p6->cd(); p6->SetFillColor(42); p6->SetFrameFillColor(10);
415 fE->SetFillColor(2); fE->SetFillStyle(3005);
416 fE->SetXTitle("Arbitrary Units");
417 if (fE->GetEntries() < minc) fE->Draw(); else fE->Fit("gaus"); c2->cd();
419 TPad* p7 = new TPad("p7", "", 0., 0.5, 1., 1.); p7->Draw();
420 p7->cd(); p7->SetFillColor(42); p7->SetFrameFillColor(10);
421 fEp->SetMarkerStyle(8); fEp->SetMarkerSize(0.4);
422 fEp->SetFillColor(2); fEp->SetFillStyle(3005);
423 fEp->SetXTitle("p (GeV/c)"); fEp->SetYTitle("dE/dX (Arb. Units)");
424 fEp->Draw(); c1->cd();
426 TCanvas* c3 = new TCanvas("c3", "", 10, 10, 510, 510);
427 c3->SetFillColor(42); c3->SetFrameFillColor(10);
428 fFoundPhi->Sumw2(); fGoodPhi->Sumw2();
429 TH1F* hgphi = new TH1F("hgphi", "Efficiency for good tracks (Phi)",
430 90, 0., 2.*TMath::Pi());
431 hgphi->Divide(fFoundPhi, fGoodPhi, 1., 1., "B");
432 hgphi->SetYTitle("Tracking efficiency");
433 hgphi->SetXTitle("Phi (rad)");
436 TFile fc("AliTRDComparison.root","RECREATE");