A lot of changes after detector review:
[u/mrichter/AliRoot.git] / FMD / scripts / Poisson.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 <AliESDFMD.h>
15 #include <AliFMDInput.h>
16 #include <TH2D.h>
17 #include <TStyle.h>
18 #include <TArrayF.h>
19 #include <TCanvas.h>
20 #include <TStyle.h>
21 #include <TROOT.h>
22 #include <TFile.h>
23 #include <iostream>
24
25 /** @class Poisson
26     @brief Make a poisson reconstruction
27     @code 
28     Root> .L Compile.C
29     Root> Compile("Poisson.C")
30     Root> Poisson c
31     Root> c.Run();
32     @endcode
33     @ingroup FMD_script
34  */
35 class Poisson : public AliFMDInput
36 {
37 protected:
38   TH2D*  fEmpty; // Histogram 
39   TH2D*  fTotal; // Histogram 
40   TH2D*  fMult;  // Histogram 
41   TFile* fFile;  // File 
42   Int_t  fEv;    // Event number
43   Double_t fThreshold;
44 public:
45   /** Constructor 
46       @param threshold Threshold
47       @param nEta      # of @f$ \eta@f$ bins
48       @param minEta    minimum @f$ \eta@f$
49       @param maxEta    maximum @f$ \eta@f$
50       @param nPhi      # of @f$ \eta@f$ bins
51       @param minPhi    minimum @f$ \varphi@f$  
52       @param maxPhi    maximum @f$ \varphi@f$ */
53   Poisson(Double_t threshold=.3,
54           Int_t nEta=120, Float_t minEta=-6, Float_t maxEta=6, 
55           Int_t nPhi=4,   Float_t minPhi=0,  Float_t maxPhi=2*TMath::Pi())
56     : fFile(0), fEv(0), fThreshold(threshold)
57   { 
58     AddLoad(kESD);
59
60     fEmpty = new TH2D("empty", "# of empty strips", nEta, minEta, maxEta, 
61                       nPhi, minPhi, maxPhi);
62     fTotal = new TH2D("total", "Total # of strips", nEta, minEta, maxEta, 
63                       nPhi, minPhi, maxPhi);
64     fMult = new TH2D("mult", "Multiplicity", nEta, minEta, maxEta, 
65                       nPhi, minPhi, maxPhi);
66     fEmpty->SetXTitle("#eta");  
67     fEmpty->SetYTitle("#phi"); 
68     fEmpty->SetZTitle("N");
69     fTotal->SetXTitle("#eta");  
70     fTotal->SetYTitle("#phi"); 
71     fTotal->SetZTitle("N");
72     fMult->SetXTitle("#eta");  
73     fMult->SetYTitle("#phi"); 
74     fMult->SetZTitle("<M_{ch}>");
75
76   }
77   /** Initialize the analyser. Opens the output file. 
78       @return @c true on success. */
79   virtual Bool_t Init() 
80   {
81     if (!AliFMDInput::Init()) return kFALSE;
82     fFile = TFile::Open("poisson.root", "RECREATE");
83     if (!fFile) return kFALSE;
84     return kTRUE;
85   }
86   /** Begining of event
87       @param event Event number
88       @return @c false on error */
89   virtual Bool_t Begin(Int_t event) 
90   {
91     if (!AliFMDInput::Begin(event)) return kFALSE;
92     fEv = event;
93     fEmpty->Clear();
94     fTotal->Clear();
95     fMult->Clear();
96     return kTRUE;
97   }
98   /** Process ESD data.  For each strip, check if the
99       psuedo-multiplicity is less than the threshold.  If it is, then
100       count the strip as empty. 
101       @param esd ESD data 
102       @return @c true on success. */
103   virtual Bool_t ProcessESD(AliESDFMD* esd)
104   {
105     for (UShort_t det = 1; det <= 3; det++) {
106       for (UShort_t rng = 0; rng < 2; rng++) {
107         Char_t ring = (rng == 0 ? 'I' : 'O');
108         // Not covered channels 
109         for (UShort_t sec = 0; sec < 40; sec++) {
110           for (UShort_t str = 0; str < 512; str++) {
111             Float_t mult = esd->Multiplicity(det, ring, sec, str);
112             Float_t eta  = esd->Eta(det, ring, sec, str);
113             // Dead channels, or not covered. 
114             if (mult >= AliESDFMD::kInvalidMult) continue;
115             if (eta  >= AliESDFMD::kInvalidEta) continue;
116             Float_t phi;
117             switch (ring) { 
118             case 'I':  phi = (sec + .5) * 2 * TMath::Pi() / 20; break;
119             case 'O':  phi = (sec + .5) * 2 * TMath::Pi() / 40; break;
120             }
121             fTotal->Fill(eta, phi);
122             if (mult < fThreshold) fEmpty->Fill(eta, phi);
123           } // Loop over strips
124         } // Loop over sectors
125       } // Loop over rings
126     } // Loop over detectors
127     return kTRUE;
128   }
129   /** For each bin, reconstruct the charge particle multiplicity as 
130       @f[
131       m = - N_{total} \log\left(\frac{N_{empty}}{N_{total}}\right)
132       @f]
133       where @f$ N_{total}@f$ is the total number of strips in the bin,
134       and @f$ N_{empty}@f$ is the number of strips in the bin that did
135       not fire. 
136       @return @c true  */
137   virtual Bool_t End() 
138   {
139     for (Int_t etaBin = 1; etaBin <= fEmpty->GetNbinsX(); etaBin++) {
140       for (Int_t phiBin = 1; phiBin <= fEmpty->GetNbinsY(); phiBin++) {
141         Double_t empty  = fEmpty->GetBinContent(etaBin, phiBin);
142         Double_t total  = fTotal->GetBinContent(etaBin, phiBin);
143         Double_t lambda = (empty > 0 ? - TMath::Log(empty / total) : 1);
144         Double_t mult   = lambda * total;
145         fMult->SetBinContent(etaBin, phiBin, mult);
146       }
147     }
148     fFile->cd();
149     fMult->Write(Form("mult%03d", fEv));
150     if (!gROOT->IsBatch()) { 
151       gStyle->SetPalette(1);
152       TCanvas* c = new TCanvas("poisson", "Poisson multiplicity");
153       c->SetFillColor(0);
154       fMult->Draw("colz");
155     }
156     return AliFMDInput::End();
157   }
158   
159   /** At end of run.  Write and close output file @c poisson.root 
160       @return @c true on success */
161   virtual Bool_t Finish()
162   {
163     fFile->Write();
164     fFile->Close();
165     fFile = 0;
166     return AliFMDInput::Finish();
167   }
168   ClassDef(Poisson,0);
169 };
170
171 //____________________________________________________________________
172 //
173 // EOF
174 //