Addes script to compare Naiive, Poisson to Hits, Primaries
[u/mrichter/AliRoot.git] / FMD / scripts / CompareMethods.C
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
26 Double_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
41 void 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 }