Updated version of the proof-enabled comparison macros: now available for ITS, TPC...
[u/mrichter/AliRoot.git] / PWG1 / comparison / AliTPCComparisonTask.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 "AliTPCComparisonTask.h"
26
27 ClassImp(AliTPCComparisonTask)
28
29
30 AliTPCComparisonTask::AliTPCComparisonTask()
31   : AliAnalysisTaskSE("AliTPCComparisonTask"),
32     fListOfHistos(0),
33     hgood(0),
34     hfound(0),
35     hfake(0),
36     hp(0),
37     hl(0),
38     hpt(0),
39     hmpt(0),
40     he(0),
41     hep(0),
42     hgoodPhi(0),
43     hfoundPhi(0)
44 {
45   // Default constructor
46   AliInfo("Default constructor AliTPCComparisonTask");
47   // Define input and output slots here
48   // Input slot #0 works with a TChain
49   DefineInput(0, TChain::Class());
50   // Output slot #1 TList
51   DefineOutput(1, TList::Class());
52 }
53
54
55 AliTPCComparisonTask::AliTPCComparisonTask(const char* name)
56   : AliAnalysisTaskSE(name),
57     fListOfHistos(0),
58     hgood(0),
59     hfound(0),
60     hfake(0),
61     hp(0),
62     hl(0),
63     hpt(0),
64     hmpt(0),
65     he(0),
66     hep(0),
67     hgoodPhi(0),
68     hfoundPhi(0)
69 {
70   // Constructor
71   AliInfo("Constructor AliTPCComparisonTask");
72   // Define input and output slots here
73   // Input slot #0 works with a TChain
74   DefineInput(0, TChain::Class());
75   // Output slot #1 TList
76   DefineOutput(1, TList::Class());
77 }
78
79
80 void AliTPCComparisonTask::UserCreateOutputObjects() 
81 {
82   // Create histograms
83   // Called once
84   AliInfo("AliTPCComparisonTask::UserCreateOutputObjects");
85   // Create output container
86   fListOfHistos = new TList();
87   
88   hgood = new TH1F("hgood", "Pt for good tracks", 34, 0.2, 7.0);
89   hfound = new TH1F("hfound", "Pt for found tracks", 34, 0.2, 7.0);
90   hfake = new TH1F("hfake",  "Pt for fake tracks", 34, 0.2, 7.0);
91   hp = new TH1F("hp", "PHI resolution", 50, -20., 20.);
92   hl = new TH1F("hl", "LAMBDA resolution", 50, -20., 20.);
93   hpt = new TH1F("hpt", "Relative Pt resolution", 30, -10., 10.);
94   hmpt = new TH1F("hmpt", "Relative Pt resolution (pt>4GeV/c)", 30, -60., 60.);
95   he = new TH1F("he", "dE/dX for pions with 0.4<p<0.5 GeV/c", 50, 0., 100.);
96   hep = new TH2F("hep", "dE/dX vs momentum", 50, 0., 2., 50, 0., 400.);
97   hgoodPhi = new TH1F("hgoodPhi", "Phi for good tracks", 90, 0., 2.*TMath::Pi());
98   hfoundPhi = new TH1F("hfoundPhi", "Phi for found tracks", 90, 0., 2.*TMath::Pi());
99   
100   fListOfHistos->Add(hgood);
101   fListOfHistos->Add(hfound);
102   fListOfHistos->Add(hfake);
103   fListOfHistos->Add(hp);
104   fListOfHistos->Add(hl);
105   fListOfHistos->Add(hpt);
106   fListOfHistos->Add(hmpt);
107   fListOfHistos->Add(he);
108   fListOfHistos->Add(hep);
109   fListOfHistos->Add(hgoodPhi);
110   fListOfHistos->Add(hfoundPhi);
111 }
112
113
114 void AliTPCComparisonTask::UserExec(Option_t *) 
115 {
116   // Main loop
117   // Called for each event
118   
119   // MC information
120   AliMCEvent* mcEvent = MCEvent();
121   if (!mcEvent) {
122     Printf("ERROR: Could not retrieve MC event");
123     return;
124   }
125   
126   Int_t nt = 0;
127   
128   TClonesArray dummy("AliMCComparisonTrack",1000), *refs=&dummy;
129   
130   // Loop over all MC tracks and select "good" tracks
131   for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
132     AliMCParticle* track = mcEvent->GetTrack(iTracks);
133     if (!track) {
134       Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
135       continue;
136     }
137     
138     // Track selection
139     if (track->Pt() < 0.2) continue; 
140     if (track->Pt() > 7.0) continue;   
141     if (TMath::Abs(track->Pz()/track->Pt()) > 0.999) continue;
142     
143     Double_t vx = track->Xv(),vy = track->Yv(),vz = track->Zv();
144     if (TMath::Sqrt(vx*vx+vy*vy) > 3.5) continue;
145     if (TMath::Abs(vz) > 50.) continue; 
146     
147     // Loop over Track References
148     Bool_t LabelTPC = kFALSE;
149     AliTrackReference* trackRef = 0;
150     for (Int_t iTrackRef = 0; iTrackRef  < track->GetNumberOfTrackReferences(); iTrackRef++) {
151       trackRef = track->GetTrackReference(iTrackRef);
152       if(trackRef) {
153         Int_t detectorId = trackRef->DetectorId(); 
154         if (detectorId == AliTrackReference::kTPC) {        
155           LabelTPC = kTRUE;
156           break;
157         }
158       }      
159     }
160     // "Good" tracks
161     if (LabelTPC) {
162       AliMCComparisonTrack* ref = new((*refs)[nt]) AliMCComparisonTrack();
163       ref->SetLabel(iTracks);
164       TParticle* particle = track->Particle();
165       Int_t pdg = particle->GetPdgCode();
166       ref->SetPDG(pdg); 
167       ref->SetPz(trackRef->Pz());
168       Float_t Pt = trackRef->Pt();
169       ref->SetPt(Pt);
170       hgood->Fill(Pt);
171       Float_t phig = trackRef->Phi();
172       ref->SetPhi(phig);
173       hgoodPhi->Fill(phig);
174       nt++;  
175     }    
176   } //track loop 
177   
178   
179   // ESD information  
180   AliVEvent* event = InputEvent();
181   if (!event) {
182     Printf("ERROR: Could not retrieve event");
183     return;
184   }
185     
186   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
187   
188   Bool_t iFound;
189   Int_t nfound = 0;
190   Int_t nfake = 0;
191   Int_t nlost = 0;
192   
193   Int_t MCgoods = refs->GetEntriesFast();
194   
195   for (Int_t k = 0; k < MCgoods; k++) {
196     AliMCComparisonTrack* ref = (AliMCComparisonTrack*)refs->UncheckedAt(k);
197     if (!ref) continue;
198     Int_t MCLabel = ref->GetLabel();
199     Float_t ptg = ref->GetPt();
200     Float_t Phig = ref->GetPhi();
201
202     iFound = kFALSE;
203     for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
204       AliESDtrack* track = esd->GetTrack(iTrack);  
205       if (!track) {
206         Printf("ERROR: Could not receive track %d", iTrack);
207         continue;
208       }
209       
210       if (! track->IsOn(AliESDtrack::kTPCrefit)) continue;     
211       
212       Int_t TPCLabel =  track->GetTPCLabel();      
213       if (MCLabel == TMath::Abs(TPCLabel)) {      
214         if (MCLabel == TPCLabel) {
215           nfound++;
216           hfound->Fill(ptg);
217           hfoundPhi->Fill(Phig);
218         } 
219         else {
220           nfake++;
221           hfake->Fill(ptg);
222         }
223         iFound = kTRUE;
224         
225       
226         Double_t pxpypz[3]; track->GetInnerPxPyPz(pxpypz);
227         Float_t phi = TMath::ATan2(pxpypz[1],pxpypz[0]);
228         if (phi < 0) phi += 2*TMath::Pi();
229         Double_t pt = TMath::Sqrt(pxpypz[0]*pxpypz[0]+pxpypz[1]*pxpypz[1]);
230         Float_t lam = TMath::ATan2(pxpypz[2],pt); 
231         Float_t pt_1 = 1/pt;
232
233         Int_t pdg = ref->GetPDG();  
234         if (TMath::Abs(pdg) == 11 && ptg>4.) {
235           //high momentum electrons
236           hmpt->Fill((pt_1 - 1/ptg)/(1/ptg)*100.);
237         } 
238         else {
239           hp->Fill((phi - Phig)*1000.);          
240           Float_t lamg = TMath::ATan2(ref->GetPz(),ptg);
241           hl->Fill((lam - lamg)*1000.);
242           hpt->Fill((pt_1 - 1/ptg)/(1/ptg)*100.);
243         }
244         
245         Float_t mom = pt/TMath::Cos(lam);
246         Float_t dedx = track->GetTPCsignal();
247         hep->Fill(mom, dedx, 1.);
248         if (TMath::Abs(pdg) == 211) { // pions
249           if (mom > 0.4 && mom < 0.5) {
250             he->Fill(dedx,1.);
251           }
252         }   
253            
254         break; 
255       }
256     }  
257     if (!iFound) {
258       nlost++;
259     } 
260   }
261      
262   Printf(" Results: " );
263   Printf(" Found %d Fake %d Lost %d ", nfound, nfake, nlost);
264   
265   refs->Clear();
266   
267   // Post output data.
268   PostData(1, fListOfHistos);
269 }      
270
271
272 void AliTPCComparisonTask::Terminate(Option_t *) 
273 {
274   // Draw result to the screen
275   // Called once at the end of the query 
276   
277   fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
278   if (!fListOfHistos) {
279     Printf("ERROR: fListOfHistos not available");
280     return;
281   }
282
283   hgood = dynamic_cast<TH1F*>(fListOfHistos->At(0));
284   hfound = dynamic_cast<TH1F*>(fListOfHistos->At(1));
285   hfake = dynamic_cast<TH1F*>(fListOfHistos->At(2));  
286   hp = dynamic_cast<TH1F*>(fListOfHistos->At(3));
287   hl = dynamic_cast<TH1F*>(fListOfHistos->At(4));
288   hpt = dynamic_cast<TH1F*>(fListOfHistos->At(5));
289   hmpt = dynamic_cast<TH1F*>(fListOfHistos->At(6));
290   he = dynamic_cast<TH1F*>(fListOfHistos->At(7));
291   hep = dynamic_cast<TH2F*>(fListOfHistos->At(8));
292   hgoodPhi = dynamic_cast<TH1F*>(fListOfHistos->At(9));
293   hfoundPhi = dynamic_cast<TH1F*>(fListOfHistos->At(10));
294
295   gStyle->SetOptStat(111110);
296   gStyle->SetOptFit(1);
297    
298   TCanvas* c1 = new TCanvas("c1", "", 0, 0, 700, 850);
299
300   Int_t minc = 33;
301   
302   TPad* p1 = new TPad("p1", "", 0., 0.3, 0.5, 0.6); p1->Draw();
303   p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
304   hp->SetFillColor(4); hp->SetXTitle("(mrad)");
305   if (hp->GetEntries() < minc) hp->Draw(); else hp->Fit("gaus"); c1->cd();
306
307   TPad* p2 = new TPad("p2", "", 0.5, 0.3, 1, 0.6); p2->Draw();
308   p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10);
309   hl->SetFillColor(4); hl->SetXTitle("(mrad)");
310   if (hl->GetEntries() < minc) hl->Draw(); else hl->Fit("gaus"); c1->cd();
311   
312   TPad* p3 = new TPad("p3", "", 0., 0., 0.5, 0.3); p3->Draw();
313   p3->cd(); p3->SetFillColor(42); p3->SetFrameFillColor(10);
314   hpt->SetFillColor(2); hpt->SetXTitle("(%)");
315   if (hpt->GetEntries() < minc) hpt->Draw(); else hpt->Fit("gaus"); c1->cd();
316   
317   TPad* p4 = new TPad("p4", "", 0.5, 0., 1., 0.3); p4->Draw();
318   p4->cd(); p4->SetFillColor(42); p4->SetFrameFillColor(10);
319   hmpt->SetFillColor(6); hmpt->SetXTitle("(%)");
320   if (hmpt->GetEntries() < minc) hmpt->Draw(); else hmpt->Fit("gaus"); c1->cd(); 
321   
322   TPad* p5 = new TPad("p5", "", 0., 0.6, 1., 1.); p5->Draw(); p5->cd();
323   p5->SetFillColor(41); p5->SetFrameFillColor(10);
324   hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
325   TH1F* hg = new TH1F("hg", "Efficiency for good tracks", 34, 0.2, 7.0);
326   TH1F* hf = new TH1F("hf", "Efficiency for fake tracks", 34, 0.2, 7.0);
327   hg->Divide(hfound, hgood, 1., 1., "B");
328   hf->Divide(hfake, hgood, 1., 1., "B");
329   hg->SetLineColor(4); hg->SetLineWidth(2);
330   hg->SetMaximum(1.4);
331   hg->SetYTitle("Tracking efficiency");
332   hg->SetXTitle("Pt (GeV/c)");
333   hg->Draw();
334   
335   TLine* line1 = new TLine(0.2, 1.0, 7.0, 1.0); line1->SetLineStyle(4);
336   line1->Draw("same");
337   TLine* line2 = new TLine(0.2, 0.9, 7.0, 0.9); line2->SetLineStyle(4);
338   line2->Draw("same");
339   
340   hf->SetFillColor(1);
341   hf->SetFillStyle(3013);
342   hf->SetLineColor(2);
343   hf->SetLineWidth(2);
344   hf->Draw("histsame");
345   TText* text = new TText(0.461176, 0.248448, "Fake tracks");
346   text->SetTextSize(0.05);
347   text->Draw();
348   text = new TText(0.453919, 1.11408, "Good tracks");
349   text->SetTextSize(0.05);
350   text->Draw();
351   
352   TCanvas* c2 = new TCanvas("c2", "", 320, 32, 530, 590);
353   
354   TPad* p6 = new TPad("p6", "", 0., 0., 1., 0.5); p6->Draw();
355   p6->cd(); p6->SetFillColor(42); p6->SetFrameFillColor(10);
356   he->SetFillColor(2); he->SetFillStyle(3005);
357   he->SetXTitle("Arbitrary Units");
358   if (he->GetEntries() < minc) he->Draw(); else he->Fit("gaus"); c2->cd();
359   
360   TPad* p7 = new TPad("p7", "", 0., 0.5, 1., 1.); p7->Draw();
361   p7->cd(); p7->SetFillColor(42); p7->SetFrameFillColor(10);
362   hep->SetFillColor(2); hep->SetFillStyle(3005);
363   hep->SetMarkerStyle(8);
364   hep->SetMarkerSize(0.4);
365   hep->SetXTitle("p (GeV/c)"); hep->SetYTitle("dE/dX (Arb. Units)");
366   hep->Draw(); c1->cd();
367   
368   TCanvas* c3 = new TCanvas("c3", "", 10, 10, 510, 510);
369   c3->SetFillColor(42); c3->SetFrameFillColor(10);
370   hfoundPhi->Sumw2(); hgoodPhi->Sumw2();
371   TH1F* hgphi = new TH1F("hgphi", "Efficiency for good tracks (Phi)", 
372                          90, 0., 2.*TMath::Pi());
373   hgphi->Divide(hfoundPhi, hgoodPhi, 1., 1., "B"); 
374   hgphi->SetYTitle("Tracking efficiency");
375   hgphi->SetXTitle("Phi (rad)");
376   hgphi->Draw();
377   
378   TFile fc("AliTPCComparison.root", "RECREATE");
379   c1->Write();
380   c2->Write();
381   c3->Write();
382   fc.Close();
383
384 }