Fixes for reading zero-suppressed data. These should be propagated to
[u/mrichter/AliRoot.git] / FMD / scripts / PoissonHit.C
CommitLineData
c2fc1258 1//____________________________________________________________________
2//
3// $Id$
4//
5// Script that contains a class to draw hits, using the
6// AliFMDInputHits class in the util library.
7//
8// It draws the energy loss versus the p/(mq^2). It can be overlayed
9// with the Bethe-Bloc curve to show how the simulation behaves
10// relative to the expected.
11//
12// Use the script `Compile.C' to compile this class using ACLic.
13//
14#include "Poisson.C"
15#include <TMath.h>
16#include <TCanvas.h>
17#include <AliFMDHit.h>
18#include <AliFMDGeometry.h>
19
20/** @class PoissonHit
21 @brief Make a poisson reconstruction and compare to simulated hits
22 @code
23 Root> .L Compile.C
24 Root> Compile("Poisson.C")
25 Root> Compile("PoissonHit.C")
26 Root> PoissonHit c
27 Root> c.Run();
28 @endcode
29 @ingroup FMD_script
30 */
31class PoissonHit : public Poisson
32{
33protected:
34 TH2D* fHits; // Histogram
35 TH2D* fDiff; // Histogram
36public:
37 /** Constructor
38 @param threshold Threshold
39 @param nEta # of @f$ \eta@f$ bins
40 @param minEta minimum @f$ \eta@f$
41 @param maxEta maximum @f$ \eta@f$
42 @param nPhi # of @f$ \eta@f$ bins
43 @param minPhi minimum @f$ \varphi@f$
44 @param maxPhi maximum @f$ \varphi@f$ */
45 PoissonHit(Double_t threshold=.3,
46 Int_t nEta=120, Float_t minEta=-6, Float_t maxEta=6,
47 Int_t nPhi=4, Float_t minPhi=0, Float_t maxPhi=2*TMath::Pi())
48 : Poisson(threshold, nEta, minEta, maxEta, nPhi, minPhi, maxPhi)
49 {
50 AddLoad(kHits);
51 AddLoad(kGeometry);
52 fHits = new TH2D(*fEmpty);
53 fHits->SetName("hits");
54 fHits->SetTitle("# of hits");
55 fDiff = new TH2D(*fEmpty);
56 fDiff->SetName("diff");
57 fDiff->SetTitle("Difference between poisson and hits");
58 fHits->SetXTitle("#eta");
59 fHits->SetYTitle("#phi");
60 fHits->SetZTitle("N");
61 fDiff->SetXTitle("#eta");
62 fDiff->SetYTitle("#phi");
63 fDiff->SetZTitle("#frac{N_{hit}-N_{poisson}}{N_{hit}}");
64 }
65 /** Initialize the analyser. Opens the output file.
66 @return @c true on success. */
67 virtual Bool_t Init()
68 {
69 if (!Poisson::Init()) return kFALSE;
70 AliFMDGeometry::Instance()->Init();
71 AliFMDGeometry::Instance()->InitTransformations();
72 return kTRUE;
73 }
74 /** Get the @f$ \eta@f$ and @f$\varphi@f$ corresponding to the
75 spatial coordinates @f$ \mathbf{v} = (x,y,z)@f$
76 @param x X coordinate
77 @param y Y coordinate
78 @param z Z coordinate
79 @param eta Psuedo rapidity @f$ \eta@f$
80 @param phi Azimuthal angle @f$\varphi@f$
81 */
82 void PhysicalCoordinates(Double_t x, Double_t y, Double_t z,
83 Double_t& eta, Double_t& phi)
84 {
85 Double_t r, theta;
86 phi = TMath::ATan2(y, x);
87 r = TMath::Sqrt(y * y + x * x);
88 theta = TMath::ATan2(r, z);
89 eta = -TMath::Log(TMath::Tan(theta / 2));
90 if (phi < 0) phi += 2 * TMath::Pi();
91 }
92 /** Process one hit. Increment bin corresponding to strip.
93 @param hit Hit.
94 @return @c true on success */
95 virtual Bool_t ProcessHit(AliFMDHit* hit, TParticle*)
96 {
97 Double_t x, y, z;
98#if 0
99 AliFMDGeometry* geom = AliFMDGeometry::Instance();
100 geom->Detector2XYZ(hit->Detector(),hit->Ring(),hit->Sector(),
101 hit->Strip(),x,y,z);
102#else
103 x = hit->X();
104 y = hit->Y();
105 z = hit->Z();
106#endif
107 Double_t eta, phi;
108 PhysicalCoordinates(x, y, z, eta, phi);
109 fHits->Fill(eta, phi);
110 return kTRUE;
111 }
112 /** Begining of event
113 @param event Event number
114 @return @c false on error */
115 virtual Bool_t Begin(Int_t event)
116 {
117 if (!Poisson::Begin(event)) return kFALSE;
118 fHits->Clear();
119 fDiff->Clear();
120 return kTRUE;
121 }
122 /** Let the poisson code do it's job, and then compare to the
123 numberr of hits.
124 @return @c true */
125 virtual Bool_t End()
126 {
127 if (!Poisson::End()) return kFALSE;
128 fDiff->Add(fMult,fHits,-1.,1.);
129 fDiff->Divide(fHits);
130 if (!gROOT->IsBatch()) {
131 gStyle->SetPalette(1);
132 TCanvas* c1 = new TCanvas("hits", "Hit multiplicity");
133 c1->SetFillColor(0);
134 fHits->Draw("colz");
135 TCanvas* c2 = new TCanvas("diff", "Difference between Hit and poisson");
136 c2->SetFillColor(0);
137 fDiff->Draw("colz");
138 TCanvas* c3 = new TCanvas("empty", "# of Empty strips");
139 c3->SetFillColor(0);
140 fEmpty->Draw("colz");
141 TCanvas* c4 = new TCanvas("total", "Total # of strips");
142 c4->SetFillColor(0);
143 fTotal->Draw("colz");
144 }
145 return kTRUE;
146 }
147
148 ClassDef(PoissonHit,0);
149};
150
151//____________________________________________________________________
152//
153// EOF
154//