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