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