0. General code clean-up, including messages, and the like.
[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
21class Poisson : public AliFMDInput
22{
23private:
24 TH2D* fEmpty; // Histogram
25 TH2D* fTotal; // Histogram
26 TH2D* fMult; // Histogram
27 TFile* fFile; // File
28 Int_t fEv; // Event number
29public:
30 Poisson(Double_t threshold=.3,
31 Int_t nEta=120, Float_t minEta=-6, Float_t maxEta=6,
32 Int_t nPhi=4, Float_t minPhi=0, Float_t maxPhi=2*TMath::Pi())
33 : fFile(0), fEv(0)
34 {
35 AddLoad(kESD);
36
37 fEmpty = new TH2D("empty", "# of empty strips", nEta, minEta, maxEta,
38 nPhi, minPhi, maxPhi);
39 fTotal = new TH2D("total", "Total # of strips", nEta, minEta, maxEta,
40 nPhi, minPhi, maxPhi);
41 fMult = new TH2D("mult", "Multiplicity", nEta, minEta, maxEta,
42 nPhi, minPhi, maxPhi);
43 fEmpty->SetXTitle("#eta");
44 fEmpty->SetYTitle("#phi");
45 fEmpty->SetZTitle("N");
46 fTotal->SetXTitle("#eta");
47 fTotal->SetYTitle("#phi");
48 fTotal->SetZTitle("N");
49 fMult->SetXTitle("#eta");
50 fMult->SetYTitle("#phi");
51 fMult->SetZTitle("<M_{ch}>");
52
53 }
54 Bool_t Init()
55 {
56 if (!AliFMDInput::Begin(event)) return kFALSE;
57 fFile = TFile::Open("poisson.root");
58 if (!fFile) return kFALSE;
59 return kTRUE;
60 }
61 Bool_t Begin(Int_t event)
62 {
63 if (!AliFMDInput::Begin(event)) return kFALSE;
64 fEv = event;
65 fEmpty->Clear();
66 fTotal->Clear();
67 fMult->Clear();
68 return kTRUE;
69 }
70 Bool_t ProcessESD(AliESDFMDHit* esd)
71 {
72 for (UShort_t det = 1; det <= esd->MaxDetector(); det++) {
73 for (UShort_t rng = 0; rng < esd->MaxRing(); rng++) {
74 Char_t ring = (rng == 0 ? 'I' : 'O');
75 // Not covered channels
76 for (UShort_t sec = 0; sec < esd->MaxSector(); sec++) {
77 for (UShort_t str = 0; str < esd->MaxStrip(); str++) {
78 Float_t mult = esd->Multiplicity(det, ring, sec, str);
79 Float_t eta = esd->Eta(det, ring, sec, str);
80 // Dead channels, or not covered.
81 if (mult >= AliESDFMD::kInvalidMult) continue;
82 if (esd >= AliESDFMD::kInvalidEta) continue;
83 Float_t phi;
84 switch (ring) {
85 case 'I': phi = (sec + .5) * 2 * TMath::Pi() / 20; break;
86 case 'O': phi = (sec + .5) * 2 * TMath::Pi() / 40; break;
87 }
88 fTotal->Fill(eta, phi);
89 if (mult < threshold) fEmpty->Fill(eta, phi);
90 } // Loop over strips
91 } // Loop over sectors
92 } // Loop over rings
93 } // Loop over detectors
94 }
95 Bool_t End()
96 {
97 for (Int_t etaBin = 1; etaBin <= fEmpty->GetNbinsX(); etaBin++) {
98 for (Int_t phiBin = 1; phiBin <= fEmpty->GetNbinsY(); phiBin++) {
99 Double_t empty = fEmpty->GetBinContent(etaBin, phiBin);
100 Double_t total = fTotal->GetBinContent(etaBin, phiBin);
101 Double_t lambda = (empty > 0 ? - TMath::Log(empty / nTotal) : 1);
102 Double_t mult = lambda * nTotal;
103 fMult->SetBinContent(etaBin, phiBin, mult);
104 }
105 }
106 fFile->cd();
107 fMult->Write(Form("mult%03d", fEv));
108 return AliFMDInput::End();
109 }
110
111 Bool_t Finish()
112 {
113 fFile->Write();
114 fFile->Close();
115 fFile = 0;
116 return AliFMDInput::Finish();
117 }
118 ClassDef(Poisson,0);
119};
120
121//____________________________________________________________________
122//
123// EOF
124//