Addes script to compare Naiive, Poisson to Hits, Primaries
[u/mrichter/AliRoot.git] / FMD / scripts / CompareMethods.C
CommitLineData
6b92827e 1#ifndef __CINT__
2#include <TH1.h>
3#include <TTree.h>
4#include <TClonesArray.h>
5#include <AliFMD.h>
6#include <AliFMDHit.h>
7#include <AliFMDGeometry.h>
8#include <AliFMDMultStrip.h>
9#include <AliFMDMultRegion.h>
10#include <AliLoader.h>
11#include <AliRunLoader.h>
12#include <AliRun.h>
13#include <AliStack.h>
14#include <AliHeader.h>
15#include <AliGenEventHeader.h>
16#include <TCanvas.h>
17#include <TVector3.h>
18#include <TLegend.h>
19#include <TClassTable.h>
20#include <TParticle.h>
21#include <TLegend.h>
22#endif
23/* Script to compare the Naiive and Poisson method, to the number of
24 hits and primaries. */
25
26Double_t Strip2Eta(UShort_t detector, Char_t ring, UShort_t sector,
27 UShort_t strip, Double_t vz)
28{
29 // cout<<" Strip2Eta "<<detector<<" "<<ring<<" "<<strip;
30 AliFMDGeometry* fmdgeo = AliFMDGeometry::Instance();
31 Double_t x, y, z;
32 fmdgeo->Detector2XYZ(detector, ring, sector, strip, x, y, z);
33 Double_t realZ = z - vz;
34 Double_t r = TMath::Sqrt(x * x + y * y);
35 Double_t theta = TMath::ATan2(r, realZ);
36 Double_t eta = - TMath::Log(TMath::Tan(theta / 2));
37 return eta;
38}
39
40
41void CompareMethods()
42{
43 // Comparison of reconstruction Poisson and Naiive methods and real multiplicity
44
45 // Dynamically link some shared libs
46#ifdef __CINT__
47 if (gClassTable->GetID("AliRun") < 0) {
48 gROOT->LoadMacro("$ALICE_ROOT/macros/loadlibs.C");
49 loadlibs();
50 }
51#endif
52
53 Double_t bins1I[15]={3.61728,
54 3.7029,
55 3.8029,
56 3.9029,
57 4.0029,
58 4.1029,
59 4.2029,
60 4.3029,
61 4.4029,
62 4.5029,
63 4.6029,
64 4.7029,
65 4.8029,
66 4.9029,
67 5.0029};
68
69 Double_t bins2I[15]={2.28235,
70 2.35884,
71 2.45884,
72 2.55884,
73 2.65884,
74 2.75884,
75 2.85884,
76 2.95884,
77 3.05884,
78 3.15884,
79 3.25884,
80 3.35884,
81 3.45884,
82 3.55884,
83 3.65884};
84
85 Double_t bins3I[15]={-3.37566,
86 -3.27566,
87 -3.17566,
88 -3.07566,
89 -2.97566,
90 -2.87566,
91 -2.77566,
92 -2.67566,
93 -2.57566,
94 -2.47566,
95 -2.37566,
96 -2.27566,
97 -2.17566,
98 -2.07566,
99 -2.00644};
100
101 Double_t bins2O[7]={1.71408,
102 1.77662,
103 1.87662,
104 1.97662,
105 2.07662,
106 2.17662,
107 2.27662};
108 Double_t bins3O[7]={-2.27662,
109 -2.17662,
110 -2.07662,
111 -1.97662,
112 -1.87662,
113 -1.77662,
114 -1.71408};
115 //__________________________________________________________________
116 TH1F* hHits = new TH1F("hits", "Hits", 100, -5, 5);
117 TH1F* hPrimary = new TH1F("primary", "Primary", 100, -5, 5);
118 TH1F* hAll = new TH1F("all", "All", 100, -5, 5);
119 TH1F* hNaiive = new TH1F("naiive", "Naiive", 100, -5, 5);
120 TH1F* hPoisson = new TH1F("poisson", "Poisson", 100, -5, 5);
121 hNaiive->SetLineStyle(1); hNaiive->SetLineColor(2);
122 hNaiive->SetFillStyle(3004); hNaiive->SetFillColor(2);
123 hPoisson->SetLineStyle(1); hPoisson->SetLineColor(3);
124 hPoisson->SetFillStyle(3005); hPoisson->SetFillColor(3);
125 hHits->SetLineStyle(1); hHits->SetLineColor(1);
126 hHits->SetFillStyle(3001); hHits->SetFillColor(1);
127 hPrimary->SetLineStyle(1); hPrimary->SetLineColor(4);
128 hPrimary->SetFillStyle(3001); hPrimary->SetFillColor(4);
129 hAll->SetLineStyle(1); hAll->SetLineColor(6);
130 hAll->SetFillStyle(3001); hAll->SetFillColor(6);
131
132
133 //__________________________________________________________________
134 TH1F *h1IHits = new TH1F ("h1IHits", "Hits in FMD1I",
135 14,bins1I[0],bins1I[14]);
136 TH1F *h2IHits = new TH1F ("h2IHits", "Hits in FMD2I",
137 14,bins2I[0],bins2I[14]);
138 TH1F *h3IHits = new TH1F ("h3IHits", "Hits in FMD3I",
139 14,bins3I[0],bins3I[14]);
140 TH1F *h2OHits = new TH1F ("h2OHits", "Hits in FMD2O",
141 6,bins2O[0],bins2O[6]);
142 TH1F *h3OHits = new TH1F ("h3OHits", "Hits in FMD3O",
143 6,bins3O[0],bins3O[6]);
144 h1IHits->SetLineColor(1); // h1IHits->SetFillColor(1);
145 h1IHits->SetLineStyle(1); // h1IHits->SetFillStyle(3000 + 1);
146 h2IHits->SetLineColor(1); // h2IHits->SetFillColor(1);
147 h2IHits->SetLineStyle(1); // h2IHits->SetFillStyle(3000 + 1);
148 h2OHits->SetLineColor(1); // h2OHits->SetFillColor(1);
149 h2OHits->SetLineStyle(1); // h2OHits->SetFillStyle(3000 + 1);
150 h3IHits->SetLineColor(1); // h3IHits->SetFillColor(1);
151 h3IHits->SetLineStyle(1); // h3IHits->SetFillStyle(3000 + 1);
152 h3OHits->SetLineColor(1); // h3OHits->SetFillColor(1);
153 h3OHits->SetLineStyle(1); // h3OHits->SetFillStyle(3000 + 1);
154
155 //__________________________________________________________________
156 TH1F *h1INaiive = new TH1F ("h1INaiive", "Naiive in FMD1I",
157 14,bins1I[0],bins1I[14]);
158 TH1F *h2INaiive = new TH1F ("h2INaiive", "Naiive in FMD2I",
159 14,bins2I[0],bins2I[14]);
160 TH1F *h3INaiive = new TH1F ("h3INaiive", "Naiive in FMD3I",
161 14,bins3I[0],bins3I[14]);
162 TH1F *h2ONaiive = new TH1F ("h2ONaiive", "Naiive in FMD2O",
163 6,bins2O[0],bins2O[6]);
164 TH1F *h3ONaiive = new TH1F ("h3ONaiive", "Naiive in FMD3O",
165 6,bins3O[0],bins3O[6]);
166 h1INaiive->SetLineColor(2); // h1INaiive->SetFillColor(2);
167 h1INaiive->SetLineStyle(2); // h1INaiive->SetFillStyle(3000 + 2);
168 h2INaiive->SetLineColor(2); // h2INaiive->SetFillColor(2);
169 h2INaiive->SetLineStyle(2); // h2INaiive->SetFillStyle(3000 + 2);
170 h2ONaiive->SetLineColor(2); // h2ONaiive->SetFillColor(2);
171 h2ONaiive->SetLineStyle(2); // h2ONaiive->SetFillStyle(3000 + 2);
172 h3INaiive->SetLineColor(2); // h3INaiive->SetFillColor(2);
173 h3INaiive->SetLineStyle(2); // h3INaiive->SetFillStyle(3000 + 2);
174 h3ONaiive->SetLineColor(2); // h3ONaiive->SetFillColor(2);
175 h3ONaiive->SetLineStyle(2); // h3ONaiive->SetFillStyle(3000 + 2);
176
177 //__________________________________________________________________
178 TH1F *h1IPoisson = new TH1F ("h1IPoisson", "Poisson in FMD1I",
179 14,bins1I[0],bins1I[14]);
180 TH1F *h2IPoisson = new TH1F ("h2IPoisson", "Poisson in FMD2I",
181 14,bins2I[0],bins2I[14]);
182 TH1F *h3IPoisson = new TH1F ("h3IPoisson", "Poisson in FMD3I",
183 14,bins3I[0],bins3I[14]);
184 TH1F *h2OPoisson = new TH1F ("h2OPoisson", "Poisson in FMD2O",
185 6,bins2O[0],bins2O[6]);
186 TH1F *h3OPoisson = new TH1F ("h3OPoisson", "Poisson in FMD3O",
187 6,bins3O[0],bins3O[6]);
188 h1IPoisson->SetLineColor(3); // h1IPoisson->SetFillColor(3);
189 h1IPoisson->SetLineStyle(3); // h1IPoisson->SetFillStyle(3000 + 3);
190 h2IPoisson->SetLineColor(3); // h2IPoisson->SetFillColor(3);
191 h2IPoisson->SetLineStyle(3); // h2IPoisson->SetFillStyle(3000 + 3);
192 h2OPoisson->SetLineColor(3); // h2OPoisson->SetFillColor(3);
193 h2OPoisson->SetLineStyle(3); // h2OPoisson->SetFillStyle(3000 + 3);
194 h3IPoisson->SetLineColor(3); // h3IPoisson->SetFillColor(3);
195 h3IPoisson->SetLineStyle(3); // h3IPoisson->SetFillStyle(3000 + 3);
196 h3OPoisson->SetLineColor(3); // h3OPoisson->SetFillColor(3);
197 h3OPoisson->SetLineStyle(3); // h3OPoisson->SetFillStyle(3000 + 3);
198
199 //__________________________________________________________________
200 TH1F *h1IPrimary = new TH1F ("h1IPrimary", "Primary in FMD1I",
201 14,bins1I[0],bins1I[14]);
202 TH1F *h2IPrimary = new TH1F ("h2IPrimary", "Primary in FMD2I",
203 14,bins2I[0],bins2I[14]);
204 TH1F *h3IPrimary = new TH1F ("h3IPrimary", "Primary in FMD3I",
205 14,bins3I[0],bins3I[14]);
206 TH1F *h2OPrimary = new TH1F ("h2OPrimary", "Primary in FMD2O",
207 6,bins2O[0],bins2O[6]);
208 TH1F *h3OPrimary = new TH1F ("h3OPrimary", "Primary in FMD3O",
209 6,bins3O[0],bins3O[6]);
210 TH1F *hEdep = new TH1F("hEdep"," edep",100,0,1);
211 h1IPrimary->SetLineColor(4); // h1IPrimary->SetFillColor(4);
212 h1IPrimary->SetLineStyle(4); // h1IPrimary->SetFillStyle(3000 + 4);
213 h2IPrimary->SetLineColor(4); // h2IPrimary->SetFillColor(4);
214 h2IPrimary->SetLineStyle(4); // h2IPrimary->SetFillStyle(3000 + 4);
215 h2OPrimary->SetLineColor(4); // h2OPrimary->SetFillColor(4);
216 h2OPrimary->SetLineStyle(4); // h2OPrimary->SetFillStyle(3000 + 4);
217 h3IPrimary->SetLineColor(4); // h3IPrimary->SetFillColor(4);
218 h3IPrimary->SetLineStyle(4); // h3IPrimary->SetFillStyle(3000 + 4);
219 h3OPrimary->SetLineColor(4); // h3OPrimary->SetFillColor(4);
220 h3OPrimary->SetLineStyle(4); // h3OPrimary->SetFillStyle(3000 + 4);
221
222 //__________________________________________________________________
223 AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
224 runLoader->LoadgAlice();
225 runLoader->LoadHeader();
226 runLoader->LoadKinematics();
227 gAlice = runLoader->GetAliRun();
228 // AliFMD* fmd = static_cast<AliFMD*>(gAlice->GetDetector("FMD"));
229 AliLoader* fmdLoader = runLoader->GetLoader("FMDLoader");
230 fmdLoader->LoadHits("READ");
231 fmdLoader->LoadRecPoints("READ");
232 TArrayF vtx;
233
234 Int_t nEvents = runLoader->TreeE()->GetEntries();
235 for (Int_t event = 0; event < nEvents && event < 1 ; event++) {
236 cout << "Event # " << event << flush;
237 runLoader->GetEvent(event);
238 AliHeader* header = runLoader->GetHeader();
239 AliGenEventHeader* eventHeader = header->GenEventHeader();
240 eventHeader->PrimaryVertex(vtx);
241 Double_t vz = vtx[2];
242 cout<<" Vz= "<< vz << flush;
243
244 // Fill primaries
245 Int_t nparticles= runLoader->Stack()->GetNtrack();
246 cout << " # particles=" << nparticles << flush;
247 for(Int_t ipart=0; ipart<nparticles; ipart++) {
248 TParticle *p=runLoader->Stack()->Particle(ipart);
249 Float_t etaKine=p->Eta();
250 if (p->GetMother(0) == -1) {
251 hPrimary->Fill(etaKine);
252 h1IPrimary->Fill(etaKine);
253 h2IPrimary->Fill(etaKine);
254 h3IPrimary->Fill(etaKine);
255 h2OPrimary->Fill(etaKine);
256 h3OPrimary->Fill(etaKine);
257 }
258 hAll->Fill(etaKine);
259 }
260 cout << endl;
261
262 // Get the hits and histogram them
263 TClonesArray* hits = 0;
264 TTree* treeH = fmdLoader->TreeH();
265 TBranch* branch = treeH->GetBranch("FMD");
266 branch->SetAddress(&hits);
267 Int_t totalHits = 0;
268 Int_t nHitEntries = treeH->GetEntries();
269 cout << " # hit entries=" << nHitEntries << " " << flush;
270 for (Int_t iHitEntry = 0; iHitEntry < nHitEntries ; iHitEntry++) {
271 treeH->GetEntry(iHitEntry);
272
273 if (iHitEntry % 1000 == 0) cout << "." << flush;
274 Int_t nHits = hits->GetEntries();
275 for (Int_t ihit = 0; ihit < nHits; ihit++) {
276 AliFMDHit* hit = static_cast<AliFMDHit*>(hits->UncheckedAt(ihit));
277 UShort_t detector = hit->Detector();
278 Char_t ring = hit->Ring();
279 UShort_t sector = hit->Sector();
280 UShort_t strip = hit->Strip();
281 Float_t edep = hit->Edep();
282 if (edep >0 ) hEdep->Fill(edep);
283 Float_t eta = Strip2Eta(detector, ring, sector, strip, vz);
284 if (edep > 0.0136) {
285 totalHits++;
286 hHits->Fill(eta);
287 switch (detector) {
288 case 1: h1IHits->Fill(eta); break;
289 case 2:
290 switch (ring){
291 case 'I': h2IHits->Fill(eta); break;
292 case 'O': h2OHits->Fill(eta); break;
293 }
294 break;
295 case 3:
296 switch (ring){
297 case 'I': h3IHits->Fill(eta); break;
298 case 'O': h3OHits->Fill(eta); break;
299 }
300 break;
301 }
302 } // if (edep > 0.0136)
303 } // for (Int_t i = 0; i < nHits; i++)
304 } // for (Int_t entry = 0; entry < nEntry; entry++)
305 cout << endl;
306
307 // Get the recPoints, and histogram them
308 TClonesArray* multNaiive = 0;
309 TClonesArray* multPoisson = 0;
310 TTree* treeR = fmdLoader->TreeR();
311 TBranch* branchPoisson = treeR->GetBranch("FMDPoisson");
312 TBranch* branchNaiive = treeR->GetBranch("FMDNaiive");
313 branchPoisson->SetAddress(&multPoisson);
314 branchNaiive->SetAddress(&multNaiive);
315
316 Int_t nRecEntries = treeR->GetEntries();
317 cout << " # rec entries=" << nRecEntries << endl;
318 for (Int_t iRecEntry = 0; iRecEntry < nRecEntries; iRecEntry++) {
319 treeR->GetEntry(iRecEntry);
320
321 // Naiive reconstrution
322 Int_t nNaiives = multNaiive->GetLast();
323 cout << " # naiive=" << nNaiives << " " << flush;
324 for (Int_t inaiive = 0; inaiive < nNaiives; inaiive++) {
325 if (inaiive % 1000 == 0) cout << "." << flush;
326 AliFMDMultStrip* naiive =
327 static_cast<AliFMDMultStrip*>(multNaiive->UncheckedAt(inaiive));
328 Float_t nParticles = naiive->Particles();
329 Float_t eta = naiive->Eta();
330 Char_t ring = naiive->Ring();
331 UShort_t detector = naiive->Detector();
332 hNaiive->Fill(eta, nParticles);
333 switch (detector) {
334 case 1: h1INaiive->Fill(eta,nParticles); break;
335 case 2:
336 switch (ring){
337 case 'i':
338 case 'I': h2INaiive->Fill(eta,nParticles); break;
339 case 'o':
340 case 'O': h2ONaiive->Fill(eta,nParticles); break;
341 }
342 break;
343 case 3:
344 switch (ring){
345 case 'i':
346 case 'I': h3INaiive->Fill(eta,nParticles); break;
347 case 'o':
348 case 'O': h3ONaiive->Fill(eta,nParticles); break;
349 }
350 break;
351 }
352 } // for (inrec = 0; ...)
353 cout << endl;
354
355 // Poisson reconstruction
356 Int_t nPoissons = multPoisson->GetLast();
357 cout << " # poisson=" << nPoissons << " " << flush;
358 for (Int_t ipoisson = 0; ipoisson <= nPoissons; ipoisson++) {
359 cout << "." << flush;
360 AliFMDMultRegion* poisson =
361 static_cast<AliFMDMultRegion*>(multPoisson->UncheckedAt(ipoisson));
362
363 Float_t nParticles = poisson->Particles();
364 UShort_t detector = poisson->Detector();
365 Char_t ring = poisson->Ring();
366 Float_t eta = (poisson->MaxEta() + poisson->MinEta()) / 2;
367 // poisson->Print("EPTD");
368 // cout << " Eta returned: " << eta << endl;
369 hPoisson->Fill(eta, nParticles);
370 switch (detector) {
371 case 1: h1IPoisson->Fill(eta,nParticles); break;
372 case 2:
373 switch (ring){
374 case 'i':
375 case 'I': h2IPoisson->Fill(eta,nParticles); break;
376 case 'o':
377 case 'O': h2OPoisson->Fill(eta,nParticles); break;
378 }
379 break;
380 case 3:
381 switch (ring){
382 case 'i':
383 case 'I': h3IPoisson->Fill(eta,nParticles); break;
384 case 'o':
385 case 'O': h3OPoisson->Fill(eta,nParticles); break;
386 }
387 break;
388 }
389 } // for (ipoisson = 0; ...)
390 cout << endl;
391 }
392 }
393 cout << "All done, drawing ... " << endl;
394
395 TCanvas* nullc = new TCanvas("nullc", "Digit Data");
396 TH1* null1 = new TH1D("null1", "null", 10, 0, 1);
397 null1->Draw();
398
399 TCanvas* c2 = new TCanvas("hit", "Hit Data");
400 c2->cd();
401 hAll->Draw();
402 hPrimary->Draw("same");
403 hHits->Draw("same");
404 hNaiive->Draw("same");
405 hPoisson->Draw("same");
406 TLegend* l = new TLegend(.7, .7, 1, 1);
407 l->AddEntry(hAll, "All", "L");
408 l->AddEntry(hPrimary, "Primary", "L");
409 l->AddEntry(hHits, "Hits", "L");
410 l->AddEntry(hNaiive, "Naiive", "L");
411 l->AddEntry(hPoisson, "Poisson", "L");
412 l->Draw();
413
414
415 TCanvas* c1 = new TCanvas("perDetector", "Per detector", 1200, 800);
416 c1->Divide(3,2);
417 c1->cd(6);
418 TH1* null = new TH1D("null", "null", 10, 0, 1);
419 null->Draw("A");
420 TLegend* l2 = new TLegend(0, 0, 1, 1);
421 l2->AddEntry(h1INaiive, "Naiive", "L");
422 l2->AddEntry(h1IPoisson, "Poission", "L");
423 l2->AddEntry(h1IHits, "Hits", "L");
424 l2->AddEntry(h1IPrimary, "Primary", "L");
425 l2->Draw();
426
427 c1->cd(3);
428 h1IPrimary->SetMinimum(0); h1IPrimary->SetMaximum(300);
429 h1IPrimary->Draw();;
430 h1IHits->Draw("same");
431 h1IPoisson->Draw("same");
432 h1INaiive->Draw("same");
433
434 c1->cd(5);
435 h2OPrimary->SetMinimum(0); h2OPrimary->SetMaximum(300);
436 h2OPrimary->Draw();;
437 h2OHits->Draw("same");
438 h2OPoisson->Draw("same");
439 h2ONaiive->Draw("same");
440
441 c1->cd(2);
442 h2IPrimary->SetMinimum(0); h2IPrimary->SetMaximum(300);
443 h2IPrimary->Draw();;
444 h2IHits->Draw("same");
445 h2IPoisson->Draw("same");
446 h2INaiive->Draw("same");
447
448
449 c1->cd(4);
450 h3OPrimary->SetMinimum(0); h3OPrimary->SetMaximum(300);
451 h3OPrimary->Draw();;
452 h3OHits->Draw("same");
453 h3ONaiive->Draw("same");
454 h3OPoisson->Draw("same");
455
456 c1->cd(1);
457 h3IPrimary->SetMinimum(0); h3IPrimary->SetMaximum(300);
458 h3IPrimary->Draw();;
459 h3IHits->Draw("same");
460 h3INaiive->Draw("same");
461 h3IPoisson->Draw("same");
462
463
464
465 TFile *fileHist = new TFile("compare2000.root","RECREATE");
466 hHits->Write();
467 hPrimary->Write();
468 hNaiive->Write();
469 hPoisson->Write();
470 hAll->Write();
471
472 h1INaiive->Write();
473 h2INaiive->Write();
474 h3INaiive->Write();
475 h2ONaiive->Write();
476 h3ONaiive->Write();
477 h1IPoisson->Write();
478 h2IPoisson->Write();
479 h3IPoisson->Write();
480 h2OPoisson->Write();
481 h3OPoisson->Write();
482 h1IHits->Write();
483 h2IHits->Write();
484 h3IHits->Write();
485 h2OHits->Write();
486 h3OHits->Write();
487 h1IPrimary->Write();
488 h2IPrimary->Write();
489 h3IPrimary->Write();
490 h2OPrimary->Write();
491 h3OPrimary->Write();
492 hEdep->Write();
493 nullc->Close();
494
495}