Fixes for reading zero-suppressed data. These should be propagated to
[u/mrichter/AliRoot.git] / FMD / scripts / PoissonHit.C
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  */
31 class PoissonHit : public Poisson
32 {
33 protected:
34   TH2D*  fHits; // Histogram 
35   TH2D*  fDiff; // Histogram 
36 public:
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 //