A lot of changes after detector review:
[u/mrichter/AliRoot.git] / FMD / scripts / Poisson.C
CommitLineData
bf000c32 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>
c2fc1258 19#include <TCanvas.h>
20#include <TStyle.h>
21#include <TROOT.h>
22#include <TFile.h>
bf000c32 23#include <iostream>
24
9b48326f 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 */
bf000c32 35class Poisson : public AliFMDInput
36{
c2fc1258 37protected:
bf000c32 38 TH2D* fEmpty; // Histogram
39 TH2D* fTotal; // Histogram
40 TH2D* fMult; // Histogram
41 TFile* fFile; // File
42 Int_t fEv; // Event number
c2fc1258 43 Double_t fThreshold;
bf000c32 44public:
c2fc1258 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$ */
bf000c32 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())
c2fc1258 56 : fFile(0), fEv(0), fThreshold(threshold)
bf000c32 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 }
c2fc1258 77 /** Initialize the analyser. Opens the output file.
78 @return @c true on success. */
79 virtual Bool_t Init()
bf000c32 80 {
c2fc1258 81 if (!AliFMDInput::Init()) return kFALSE;
82 fFile = TFile::Open("poisson.root", "RECREATE");
bf000c32 83 if (!fFile) return kFALSE;
84 return kTRUE;
85 }
9b48326f 86 /** Begining of event
87 @param event Event number
88 @return @c false on error */
c2fc1258 89 virtual Bool_t Begin(Int_t event)
bf000c32 90 {
91 if (!AliFMDInput::Begin(event)) return kFALSE;
92 fEv = event;
93 fEmpty->Clear();
94 fTotal->Clear();
95 fMult->Clear();
96 return kTRUE;
97 }
c2fc1258 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)
bf000c32 104 {
c2fc1258 105 for (UShort_t det = 1; det <= 3; det++) {
106 for (UShort_t rng = 0; rng < 2; rng++) {
bf000c32 107 Char_t ring = (rng == 0 ? 'I' : 'O');
108 // Not covered channels
c2fc1258 109 for (UShort_t sec = 0; sec < 40; sec++) {
110 for (UShort_t str = 0; str < 512; str++) {
bf000c32 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;
c2fc1258 115 if (eta >= AliESDFMD::kInvalidEta) continue;
bf000c32 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);
c2fc1258 122 if (mult < fThreshold) fEmpty->Fill(eta, phi);
bf000c32 123 } // Loop over strips
124 } // Loop over sectors
125 } // Loop over rings
126 } // Loop over detectors
c2fc1258 127 return kTRUE;
bf000c32 128 }
9b48326f 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 */
c2fc1258 137 virtual Bool_t End()
bf000c32 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);
c2fc1258 143 Double_t lambda = (empty > 0 ? - TMath::Log(empty / total) : 1);
144 Double_t mult = lambda * total;
bf000c32 145 fMult->SetBinContent(etaBin, phiBin, mult);
146 }
147 }
148 fFile->cd();
149 fMult->Write(Form("mult%03d", fEv));
c2fc1258 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 }
bf000c32 156 return AliFMDInput::End();
157 }
158
c2fc1258 159 /** At end of run. Write and close output file @c poisson.root
160 @return @c true on success */
161 virtual Bool_t Finish()
bf000c32 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//