Additiona TOF data members (S.Arcelli)
[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>
19#include <iostream>
20
9b48326f 21/** @class Poisson
22 @brief Make a poisson reconstruction
23 @code
24 Root> .L Compile.C
25 Root> Compile("Poisson.C")
26 Root> Poisson c
27 Root> c.Run();
28 @endcode
29 @ingroup FMD_script
30 */
bf000c32 31class Poisson : public AliFMDInput
32{
33private:
34 TH2D* fEmpty; // Histogram
35 TH2D* fTotal; // Histogram
36 TH2D* fMult; // Histogram
37 TFile* fFile; // File
38 Int_t fEv; // Event number
39public:
40 Poisson(Double_t threshold=.3,
41 Int_t nEta=120, Float_t minEta=-6, Float_t maxEta=6,
42 Int_t nPhi=4, Float_t minPhi=0, Float_t maxPhi=2*TMath::Pi())
43 : fFile(0), fEv(0)
44 {
45 AddLoad(kESD);
46
47 fEmpty = new TH2D("empty", "# of empty strips", nEta, minEta, maxEta,
48 nPhi, minPhi, maxPhi);
49 fTotal = new TH2D("total", "Total # of strips", nEta, minEta, maxEta,
50 nPhi, minPhi, maxPhi);
51 fMult = new TH2D("mult", "Multiplicity", nEta, minEta, maxEta,
52 nPhi, minPhi, maxPhi);
53 fEmpty->SetXTitle("#eta");
54 fEmpty->SetYTitle("#phi");
55 fEmpty->SetZTitle("N");
56 fTotal->SetXTitle("#eta");
57 fTotal->SetYTitle("#phi");
58 fTotal->SetZTitle("N");
59 fMult->SetXTitle("#eta");
60 fMult->SetYTitle("#phi");
61 fMult->SetZTitle("<M_{ch}>");
62
63 }
64 Bool_t Init()
65 {
66 if (!AliFMDInput::Begin(event)) return kFALSE;
67 fFile = TFile::Open("poisson.root");
68 if (!fFile) return kFALSE;
69 return kTRUE;
70 }
9b48326f 71 /** Begining of event
72 @param event Event number
73 @return @c false on error */
bf000c32 74 Bool_t Begin(Int_t event)
75 {
76 if (!AliFMDInput::Begin(event)) return kFALSE;
77 fEv = event;
78 fEmpty->Clear();
79 fTotal->Clear();
80 fMult->Clear();
81 return kTRUE;
82 }
83 Bool_t ProcessESD(AliESDFMDHit* esd)
84 {
85 for (UShort_t det = 1; det <= esd->MaxDetector(); det++) {
86 for (UShort_t rng = 0; rng < esd->MaxRing(); rng++) {
87 Char_t ring = (rng == 0 ? 'I' : 'O');
88 // Not covered channels
89 for (UShort_t sec = 0; sec < esd->MaxSector(); sec++) {
90 for (UShort_t str = 0; str < esd->MaxStrip(); str++) {
91 Float_t mult = esd->Multiplicity(det, ring, sec, str);
92 Float_t eta = esd->Eta(det, ring, sec, str);
93 // Dead channels, or not covered.
94 if (mult >= AliESDFMD::kInvalidMult) continue;
95 if (esd >= AliESDFMD::kInvalidEta) continue;
96 Float_t phi;
97 switch (ring) {
98 case 'I': phi = (sec + .5) * 2 * TMath::Pi() / 20; break;
99 case 'O': phi = (sec + .5) * 2 * TMath::Pi() / 40; break;
100 }
101 fTotal->Fill(eta, phi);
102 if (mult < threshold) fEmpty->Fill(eta, phi);
103 } // Loop over strips
104 } // Loop over sectors
105 } // Loop over rings
106 } // Loop over detectors
107 }
9b48326f 108 /** For each bin, reconstruct the charge particle multiplicity as
109 @f[
110 m = - N_{total} \log\left(\frac{N_{empty}}{N_{total}}\right)
111 @f]
112 where @f$ N_{total}@f$ is the total number of strips in the bin,
113 and @f$ N_{empty}@f$ is the number of strips in the bin that did
114 not fire.
115 @return @c true */
bf000c32 116 Bool_t End()
117 {
118 for (Int_t etaBin = 1; etaBin <= fEmpty->GetNbinsX(); etaBin++) {
119 for (Int_t phiBin = 1; phiBin <= fEmpty->GetNbinsY(); phiBin++) {
120 Double_t empty = fEmpty->GetBinContent(etaBin, phiBin);
121 Double_t total = fTotal->GetBinContent(etaBin, phiBin);
122 Double_t lambda = (empty > 0 ? - TMath::Log(empty / nTotal) : 1);
123 Double_t mult = lambda * nTotal;
124 fMult->SetBinContent(etaBin, phiBin, mult);
125 }
126 }
127 fFile->cd();
128 fMult->Write(Form("mult%03d", fEv));
129 return AliFMDInput::End();
130 }
131
132 Bool_t Finish()
133 {
134 fFile->Write();
135 fFile->Close();
136 fFile = 0;
137 return AliFMDInput::Finish();
138 }
139 ClassDef(Poisson,0);
140};
141
142//____________________________________________________________________
143//
144// EOF
145//