Correcly handle centrality bins (Diego)
[u/mrichter/AliRoot.git] / PWGPP / comparison / AliTOFComparisonTask.cxx
CommitLineData
18693982 1//-------------------------------------------------------------------------
2//
3// This is the PROOF-enabled version of TOF/AliTOFComparison.C macro.
4// Origin: Andrei.Zalite@cern.ch
5//
6//-------------------------------------------------------------------------
7
f2247c19 8#include "TChain.h"
9#include "TTree.h"
10#include "TH1F.h"
11#include "TH2F.h"
12#include "TList.h"
13#include "TClonesArray.h"
14#include "TCanvas.h"
15#include "TStyle.h"
16#include "TLine.h"
17#include "TText.h"
18#include "TFile.h"
f2247c19 19
20#include "AliLog.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"
27
28#include "AliTOFComparisonTask.h"
29
30ClassImp(AliTOFComparisonTask)
31
18693982 32extern TStyle *gStyle;
33
f2247c19 34AliTOFComparisonTask::AliTOFComparisonTask()
35 : AliAnalysisTaskSE("AliTOFComparisonTask"),
36 fListOfHistos(0),
18693982 37 fGood(0),
38 fFound(0),
39 fFake(0),
40 fGoodPhi(0),
41 fFoundPhi(0),
42 fGoodl(0),
43 fFakel(0),
44 fFoundl(0)
f2247c19 45{
46 // Default constructor
47 AliInfo("Defaulut Constructor AliTOFComparisonTask");
48 // Define input and output slots here
49 // Input slot #0 works with a TChain
50 DefineInput(0, TChain::Class());
51 // Output slot #1 TList
52 DefineOutput(1, TList::Class());
53}
54
55AliTOFComparisonTask::AliTOFComparisonTask(const char* name)
56 : AliAnalysisTaskSE(name),
57 fListOfHistos(0),
18693982 58 fGood(0),
59 fFound(0),
60 fFake(0),
61 fGoodPhi(0),
62 fFoundPhi(0),
63 fGoodl(0),
64 fFakel(0),
65 fFoundl(0)
f2247c19 66{
67 // Constructor
68 AliInfo("Constructor AliTOFComparisonTask");
69 // Define input and output slots here
70 // Input slot #0 works with a TChain
71 DefineInput(0, TChain::Class());
72 // Output slot #1 TList
73 DefineOutput(1, TList::Class());
74}
75
76
77
78void AliTOFComparisonTask::UserCreateOutputObjects()
79{
80 // Create histograms
81 // Called once
82 AliInfo("AliTOFComparisonTask::UserCreateOutputObjects");
83 // Create output container
84 fListOfHistos = new TList();
85
18693982 86 fGood = new TH1F("fGood", "Pt for good tracks", 34, 0.2, 7.0);
87 fFound = new TH1F("fFound", "Pt for found tracks", 34, 0.2, 7.0);
88 fFake = new TH1F("fFake", "Pt for fake tracks", 34, 0.2, 7.0);
89 fGoodPhi = new TH1F("fGoodPhi", "Phi for Good tracks", 90, 0., 2.*TMath::Pi());
90 fFoundPhi = new TH1F("fFoundPhi", "Phi for Found tracks", 90, 0., 2.*TMath::Pi());
91 fGoodl = new TH1F("fGoodl", "Good tracks", 30, -1., 1.);
92 fFakel = new TH1F("fFakel", "Mismatched tracks", 30, -1., 1.);
93 fFoundl = new TH1F("fFoundl", "Matched tracks", 30, -1., 1.);
f2247c19 94
18693982 95 fListOfHistos->Add(fGood);
96 fListOfHistos->Add(fFound);
97 fListOfHistos->Add(fFake);
98 fListOfHistos->Add(fGoodPhi);
99 fListOfHistos->Add(fFoundPhi);
100 fListOfHistos->Add(fGoodl);
101 fListOfHistos->Add(fFakel);
102 fListOfHistos->Add(fFoundl);
f2247c19 103}
104
105
106void AliTOFComparisonTask::UserExec(Option_t *)
107{
108 // Main loop
109 // Called for each event
110
111 // MC information
112 AliMCEvent* mcEvent = MCEvent();
113 if (!mcEvent) {
114 Printf("ERROR: Could not retrieve MC event");
115 return;
116 }
117
118 Int_t nt = 0;
119
120 TClonesArray dummy("AliMCComparisonTrack",1000), *refs=&dummy;
121
122 // Loop over all MC tracks and select "good" tracks
123 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
84933864 124 AliMCParticle* track = (AliMCParticle*)mcEvent->GetTrack(iTracks);
f2247c19 125 if (!track) {
126 Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
127 continue;
128 }
129
130 // Track selection
131 if (track->Pt() < 0.2) continue;
132 if (track->Pt() > 7.0) continue;
133 if (TMath::Abs(track->Pz()/track->Pt()) > 0.999) continue;
134
135 Double_t vx = track->Xv(), vy = track->Yv(), vz = track->Zv();
136 if (TMath::Sqrt(vx*vx+vy*vy) > 3.5) continue;
137 if (TMath::Abs(vz) > 50.) continue;
138
139 if (TMath::Abs(vx) > 0.1) continue;
140 if (TMath::Abs(vy) > 0.1) continue;
141
142
143 // Loop over Track References
18693982 144 Bool_t labelTPC = kFALSE;
145 Bool_t labelTOF = kFALSE;
f2247c19 146 AliTrackReference* trackRef = 0;
18693982 147 UInt_t iTSLayerMap = 0;
f2247c19 148 for (Int_t iTrackRef = 0; iTrackRef < track->GetNumberOfTrackReferences(); iTrackRef++) {
149 trackRef = track->GetTrackReference(iTrackRef);
150 if(trackRef) {
151 Int_t detectorId = trackRef->DetectorId();
152 if (detectorId == AliTrackReference::kITS) {
153 Float_t Radius = trackRef->R();
154 if (Radius > 2.5 && Radius < 5.5)
18693982 155 iTSLayerMap |= 1;
f2247c19 156 else if (Radius > 5.5 && Radius < 8.5)
18693982 157 iTSLayerMap |= 1 << 1;
f2247c19 158 else if (Radius > 13.5 && Radius < 16.5)
18693982 159 iTSLayerMap |= 1 << 2;
f2247c19 160 else if (Radius > 22. && Radius < 26.)
18693982 161 iTSLayerMap |= 1 << 3;
f2247c19 162 else if (Radius > 36. && Radius < 41.)
18693982 163 iTSLayerMap |= 1 << 4;
f2247c19 164 else if (Radius > 42. && Radius < 46.)
18693982 165 iTSLayerMap |= 1 << 5;
f2247c19 166 else {
167 Printf("Wrong radius %f ", Radius);
168 return;
169 }
170
171 }
172 if (detectorId == AliTrackReference::kTPC) {
18693982 173 labelTPC = kTRUE;
f2247c19 174 }
175 if (detectorId == AliTrackReference::kTOF) {
18693982 176 labelTOF = kTRUE;
f2247c19 177 }
178 }
179 } // track references loop
180
181 // Skip tracks that passed not all ITS layers
18693982 182 if (iTSLayerMap != 0x3F) continue;
f2247c19 183
184 // "Good" tracks
18693982 185 if (labelTPC && labelTOF) {
f2247c19 186 AliMCComparisonTrack* ref = new((*refs)[nt]) AliMCComparisonTrack();
187 ref->SetLabel(iTracks);
188 TParticle* particle = track->Particle();
189 Int_t pdg = particle->GetPdgCode();
190 ref->SetPDG(pdg);
191 ref->SetPz(track->Pz());
192 Float_t Pt = track->Pt();
193 ref->SetPt(Pt);
18693982 194 fGood->Fill(Pt);
f2247c19 195 Float_t phig = track->Phi();
196 ref->SetPhi(phig);
18693982 197 fGoodPhi->Fill(phig);
f2247c19 198 Double_t tgl = track->Pz()/Pt;
18693982 199 fGoodl->Fill(tgl);
f2247c19 200 nt++;
201 }
202 } // track loop
203
204 // ESD information
205 AliVEvent* event = InputEvent();
206 if (!event) {
207 Printf("ERROR: Could not retrieve event");
208 return;
209 }
210
211 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
212
213 Bool_t iFound;
214 Int_t nfound = 0;
215 Int_t nfake = 0;
216 Int_t nlost = 0;
217
18693982 218 Int_t mcGoods = refs->GetEntriesFast();
f2247c19 219
220 // Loop over all "good" MC tracks
18693982 221 for (Int_t k = 0; k < mcGoods; k++) {
f2247c19 222 AliMCComparisonTrack* ref = (AliMCComparisonTrack*)refs->UncheckedAt(k);
223 if (!ref) continue;
224 Int_t MCLabel = ref->GetLabel();
225 Float_t ptg = ref->GetPt();
226 Float_t Phig = ref->GetPhi();
227 iFound = kFALSE;
228 Int_t iTrack;
229 AliESDtrack* track = 0;
230 for (iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
231 track = esd->GetTrack(iTrack);
232 if (!track) {
233 Printf("ERROR: Could not receive track %d", iTrack);
234 continue;
235 }
236
237 Int_t Label = track->GetLabel();
238 if (MCLabel != TMath::Abs(Label)) continue;
239
240 if (track->GetTOFsignal() <= 0.) continue;
241
242 Double_t tgl = ref->GetPz()/ref->GetPt();
243 Int_t TOFLabel[3];
244 track->GetTOFLabel(TOFLabel);
245 if (TOFLabel[0] != MCLabel)
246 if (TOFLabel[1] != MCLabel)
247 if (TOFLabel[2] != MCLabel) {
18693982 248 fFake->Fill(ptg);
249 fFakel->Fill(tgl);
f2247c19 250 nfake++;
251 break;
252 }
18693982 253 fFound->Fill(ptg);
254 fFoundl->Fill(tgl);
255 fFoundPhi->Fill(Phig);
f2247c19 256 nfound++;
257 break;
258 }
259 if(iTrack == esd->GetNumberOfTracks()) {
260 Printf("Not matched: MCLabel %d ", MCLabel);
261 nlost++;
262 if (track) {
263 Printf(" kITSout %d kTPCout %d kTRDout %d kTIME %d",
264 (track->GetStatus()&AliESDtrack::kITSout),
265 (track->GetStatus()&AliESDtrack::kTPCout),
266 (track->GetStatus()&AliESDtrack::kTRDout),
267 (track->GetStatus()&AliESDtrack::kTIME));
268 }
269 else
270 Printf(" No ESD track !");
271 }
272 }
273
274 Printf(" Results: " );
275 Printf(" Good %d Found %d Fake %d Lost %d", nt, nfound, nfake, nlost);
276
277 refs->Clear();
278
279 // Post output data.
280 PostData(1, fListOfHistos);
281
282}
283
284
285void AliTOFComparisonTask::Terminate(Option_t *)
286{
287 // Draw result to the screen
288 // Called once at the end of the query
289 fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
290 if (!fListOfHistos) {
291 Printf("ERROR: fListOfHistos not available");
292 return;
293 }
294
18693982 295 fGood = dynamic_cast<TH1F*>(fListOfHistos->At(0));
296 fFound = dynamic_cast<TH1F*>(fListOfHistos->At(1));
297 fFake = dynamic_cast<TH1F*>(fListOfHistos->At(2));
298 fGoodPhi = dynamic_cast<TH1F*>(fListOfHistos->At(3));
299 fFoundPhi = dynamic_cast<TH1F*>(fListOfHistos->At(4));
300 fGoodl = dynamic_cast<TH1F*>(fListOfHistos->At(5));
301 fFakel = dynamic_cast<TH1F*>(fListOfHistos->At(6));
302 fFoundl = dynamic_cast<TH1F*>(fListOfHistos->At(7));
f2247c19 303
304 gStyle->SetOptStat(111110);
305 gStyle->SetOptFit(1);
306
307 TCanvas* c1 = new TCanvas("c1","",0,0,600,900);
308
309 TPad* p1 = new TPad("p1", "", 0., 0.5, 1., 1.); p1->Draw();
310 p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
311
18693982 312 fFound->Sumw2(); fGood->Sumw2(); fFake->Sumw2();
f2247c19 313 TH1F* hgp = new TH1F("hgp", " ", 34, 0.2, 7.0);
314 TH1F* hfp = new TH1F("hfp", "Probability of mismatching", 34, 0.2, 7.0);
18693982 315 hgp->Divide(fFound, fGood, 1., 1., "b");
316 hfp->Divide(fFake, fGood, 1., 1., "b");
f2247c19 317 hgp->SetLineColor(4); hgp->SetLineWidth(2);
318 hgp->SetMaximum(1.4);
319 hgp->SetYTitle("Tracking efficiency");
320 hgp->SetXTitle("Pt (GeV/c)");
321 hfp->SetFillColor(1); hfp->SetFillStyle(3013); hfp->SetLineWidth(2);
322 hfp->SetLineColor(2);
323
324 hgp->Draw();
325 hfp->Draw("histsame");
326 TLine *line1 = new TLine(0.2,1.0,7.0,1.0); line1->SetLineStyle(4);
327 line1->Draw("same");
328 TLine *line2 = new TLine(0.2,0.9,7.0,0.9); line2->SetLineStyle(4);
329 line2->Draw("same");
330 c1->cd();
331
332 TPad* p2 = new TPad("p2", "", 0., 0., 1., 0.5); p2->Draw();
333 p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10);
334
18693982 335 fFoundl->Sumw2(); fGoodl->Sumw2(); fFakel->Sumw2();
f2247c19 336 TH1F* hgl = new TH1F("hgl", "", 30, -1., 1.);
337 TH1F* hfl = new TH1F("hfl", "Probability of mismatching", 30, -1., 1.);
18693982 338 hgl->Divide(fFoundl, fGoodl, 1., 1., "b");
339 hfl->Divide(fFakel, fGoodl, 1., 1., "b");
f2247c19 340 hgl->SetLineColor(4); hgl->SetLineWidth(2);
341 hgl->SetMaximum(1.4);
342 hgl->SetYTitle("Tracking efficiency");
343 hgl->SetXTitle("Tan(lambda)");
344 hfl->SetFillColor(1); hfl->SetFillStyle(3013); hfl->SetLineWidth(2);
345 hfl->SetLineColor(2);
346
347 hgl->Draw();
348 hfl->Draw("histsame");
349 TLine *line3 = new TLine(-1,1.0,1,1.0); line3->SetLineStyle(4);
350 line3->Draw("same");
351 TLine *line4 = new TLine(-1,0.9,1,0.9); line4->SetLineStyle(4);
352 line4->Draw("same");
353
354 c1->Update();
355
356 TCanvas* c2 = new TCanvas("c2", "", 10, 10, 510, 510);
357 c2->SetFillColor(42); c2->SetFrameFillColor(10);
18693982 358 fFoundPhi->Sumw2(); fGoodPhi->Sumw2();
f2247c19 359 TH1F* hgphi = new TH1F("hgphi", "Efficiency for good tracks (Phi)",
360 90, 0., 2.*TMath::Pi());
18693982 361 hgphi->Divide(fFoundPhi, fGoodPhi, 1., 1., "B");
f2247c19 362 hgphi->SetYTitle("Tracking efficiency");
363 hgphi->SetXTitle("Phi (rad)");
364 hgphi->Draw();
365
366 TFile fc("AliTOFComparison.root","RECREATE");
367 c1->Write();
368 c2->Write();
369 fc.Close();
370}