Use VMC id's rather than TGeo id's
[u/mrichter/AliRoot.git] / FMD / scripts / DrawHitsRecs.C
CommitLineData
bf000c32 1//____________________________________________________________________
8f6ee336 2//
3// $Id$
4//
5// Script that contains a class to draw eloss from hits, versus ADC
6// counts from digits, using the 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 <TH2D.h>
15#include <AliFMDHit.h>
8f6ee336 16#include <AliFMDDigit.h>
15b17c89 17#include <AliFMDRecPoint.h>
8f6ee336 18#include <AliFMDInput.h>
19#include <AliFMDEdepMap.h>
20#include <AliFMDFloatMap.h>
21#include <iostream>
22#include <TStyle.h>
23#include <TArrayF.h>
24#include <TCanvas.h>
25#include <TParticle.h>
26#include <TClonesArray.h>
27#include <TTree.h>
28#include <AliStack.h>
29#include <AliLog.h>
30
bf000c32 31//____________________________________________________________________
9b48326f 32/** @class DrawHitsRecs
33 @brief Draw hit energy loss versus rec point mult
34 @code
35 Root> .L Compile.C
36 Root> Compile("DrawHitsRecs.C")
37 Root> DrawHitsRecs c
38 Root> c.Run();
39 @endcode
40 @ingroup FMD_script
41 */
15b17c89 42class DrawHitsRecs : public AliFMDInput
8f6ee336 43{
44private:
15b17c89 45 TH2D* fHitEvsAdc;
46 TH2D* fHitEvsRecM; // Histogram
47 TH2D* fHitEvsRecE; // Histogram
48 TH1D* fDiffE; // Histogram
49 TH2D* fHitsVsRecM; // Histogram
f95a63c4 50 TH2D* fDiffM; // Histogram
15b17c89 51 AliFMDEdepMap fMap;
8f6ee336 52 AliFMDFloatMap fEta;
53 AliFMDFloatMap fPhi;
54 AliFMDFloatMap fMult;
55 Bool_t fPrimary;
56public:
bf000c32 57 //__________________________________________________________________
8f6ee336 58 TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max)
59 {
60 TArrayF bins(n+1);
61 Float_t dp = n / TMath::Log10(max / min);
62 Float_t pmin = TMath::Log10(min);
63 bins[0] = min;
64 for (Int_t i = 1; i < n+1; i++) {
65 Float_t p = pmin + i / dp;
66 bins[i] = TMath::Power(10, p);
67 }
68 return bins;
69 }
bf000c32 70 //__________________________________________________________________
8f6ee336 71 DrawHitsRecs(Bool_t primary=kFALSE,
72 Int_t n=900, Double_t emin=1e-3, Double_t emax=10,
73 Int_t m=21, Double_t mmin=-0.5, Double_t mmax=20.5)
74 {
75 fPrimary = primary;
76 AddLoad(kRecPoints);
8ec606c2 77 AddLoad(kHits);
15b17c89 78 AddLoad(kDigits);
8f6ee336 79 if (fPrimary) AddLoad(kKinematics);
80 TArrayF eloss(MakeLogScale(n, emin, emax));
81 TArrayF mults(m+1);
82 mults[0] = mmin;
83 for (Int_t i = 1; i < m+1; i++) mults[i] = mults[i-1] + (mmax-mmin)/m;
15b17c89 84
85 fHitEvsAdc = new TH2D("hitEvsAdc", "#Delta E_{sim} vs. ADC",
86 n, emin, emax, 1025, -.5, 1024.5);
87 fHitEvsAdc->SetXTitle("#Delta E_{sim} [MeV]");
88 fHitEvsAdc->SetYTitle("ADC");
89
90 fHitEvsRecM = new TH2D("hitEvsRecM", "#Delta E_{sim} vs. M_{rec}",
8f6ee336 91 eloss.fN-1, eloss.fArray, mults.fN-1, mults.fArray);
15b17c89 92 fHitEvsRecM->SetXTitle("#Delta E_{sim} [MeV]");
93 fHitEvsRecM->SetYTitle("M_{rec}");
94
95 fHitEvsRecE = new TH2D("hitEvsRecE", "#Delta E_{sim} vs. #Delta E_{rec}",
96 n, emin, emax, n, emin, emax);
97 fHitEvsRecE->SetXTitle("#Delta E_{sim} [MeV]");
98 fHitEvsRecE->SetYTitle("#Delta E_{rec} [MeV]");
8f6ee336 99
15b17c89 100
101 fDiffE = new TH1D("diffE",
102 "#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}",
103 1100, -1, 1.1);
104 fDiffE->SetXTitle("#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}");
105
8f6ee336 106 Double_t omin = -.5;
107 Double_t omax = 7.5;
108 Int_t o = 8;
15b17c89 109 fHitsVsRecM = new TH2D("hitsVsStrip", "# of Hits vs. M_{rec}",
8f6ee336 110 o, omin, omax, m, mmin, mmax);
15b17c89 111 fHitsVsRecM->SetXTitle("# of Hits");
112 fHitsVsRecM->SetYTitle("M_{rec}");
f95a63c4 113
114 fDiffM = new TH2D("diffM", "M_{sim} - M_{rec}",
115 41, -20.5, 20.5, 70, 1.5, 5);
116 // 36, -TMath::Pi(),TMath::Pi());
117 fDiffM->SetXTitle("M_{sim} - M_{rec}");
118 fDiffM->SetYTitle("|#eta|");
119 // fDiffM->SetYTitle("Detector");
120
8f6ee336 121 }
bf000c32 122 //__________________________________________________________________
9b48326f 123 /** Begining of event
124 @param ev Event number
125 @return @c false on error */
8f6ee336 126 Bool_t Begin(Int_t ev)
127 {
128 fMap.Reset();
15b17c89 129 return AliFMDInput::Begin(ev);
8f6ee336 130 }
bf000c32 131 //__________________________________________________________________
8f6ee336 132 Bool_t ProcessHit(AliFMDHit* hit, TParticle*)
133 {
134 // Cache the energy loss
135 if (!hit) return kFALSE;
136 UShort_t det = hit->Detector();
137 Char_t rng = hit->Ring();
138 UShort_t sec = hit->Sector();
139 UShort_t str = hit->Strip();
140 if (str > 511) {
141 AliWarning(Form("Bad strip number %d in hit", str));
142 return kTRUE;
143 }
144 if (fPrimary) {
145 TParticle* kine = fStack->Particle(hit->Track());
146 if (!kine) return kTRUE;
147 if (!kine->IsPrimary()) return kTRUE;
148 }
149
150 fMap(det, rng, sec, str).fEdep += hit->Edep();
151 fMap(det, rng, sec, str).fN++;
152 return kTRUE;
153 }
15b17c89 154
155 //__________________________________________________________________
156 Bool_t ProcessDigit(AliFMDDigit* digit)
157 {
158 if (!digit) return kTRUE;
159
160 UShort_t det = digit->Detector();
161 Char_t rng = digit->Ring();
162 UShort_t sec = digit->Sector();
163 UShort_t str = digit->Strip();
164 if (str > 511) {
165 AliWarning(Form("Bad strip number %d in digit", str));
166 return kTRUE;
167 }
168 Double_t edep = fMap(det, rng, sec, str).fEdep;
169 if (edep > 0) fHitEvsAdc->Fill(edep, digit->Counts());
170
171 return kTRUE;
172 }
173
bf000c32 174 //__________________________________________________________________
175 Bool_t ProcessRecPoint(AliFMDRecPoint* single)
8f6ee336 176 {
15b17c89 177 if (!single) return kTRUE;
bf000c32 178 UShort_t det = single->Detector();
179 Char_t rng = single->Ring();
180 UShort_t sec = single->Sector();
181 UShort_t str = single->Strip();
182 if (str > 511) {
183 AliWarning(Form("Bad strip number %d in single", str));
15b17c89 184 return kTRUE;
185 }
186 Double_t edep = fMap(det, rng, sec, str).fEdep;
187 Int_t nhit = fMap(det, rng, sec, str).fN;
188 if (edep > 0) {
189 fHitEvsRecM->Fill(edep, single->Particles());
190 fHitEvsRecE->Fill(edep, single->Edep());
191 fDiffE->Fill((single->Edep() - edep) / edep);
8f6ee336 192 }
15b17c89 193 if (nhit > 0) fHitsVsRecM->Fill(nhit, single->Particles());
f95a63c4 194 fDiffM->Fill(nhit - single->Particles(), TMath::Abs(single->Eta()));
8f6ee336 195 return kTRUE;
196 }
bf000c32 197 //__________________________________________________________________
8f6ee336 198 Bool_t Finish()
199 {
200 gStyle->SetPalette(1);
201 gStyle->SetOptTitle(0);
202 gStyle->SetCanvasColor(0);
203 gStyle->SetCanvasBorderSize(0);
204 gStyle->SetPadColor(0);
205 gStyle->SetPadBorderSize(0);
f95a63c4 206 TCanvas* c = 0;
8f6ee336 207
f95a63c4 208 c = new TCanvas("c0", fHitEvsAdc->GetTitle());
15b17c89 209 fHitEvsAdc->SetStats(kFALSE);
210 fHitEvsAdc->Draw("COLZ");
211
f95a63c4 212 c = new TCanvas("c1", fHitEvsRecM->GetTitle());
15b17c89 213 fHitEvsRecM->SetStats(kFALSE);
214 fHitEvsRecM->Draw("COLZ");
215
f95a63c4 216 c = new TCanvas("c2", fHitEvsRecE->GetTitle());
15b17c89 217 fHitEvsRecE->SetStats(kFALSE);
218 fHitEvsRecE->Draw("COLZ");
219
f95a63c4 220 c = new TCanvas("c3", fDiffE->GetTitle());
221 c->SetLogz();
15b17c89 222 fDiffE->Draw();
8f6ee336 223
f95a63c4 224 c = new TCanvas("c4", fHitsVsRecM->GetTitle());
225 c->SetLogz();
15b17c89 226 fHitsVsRecM->SetStats(kFALSE);
227 fHitsVsRecM->Draw("COLZ");
8f6ee336 228
f95a63c4 229 c = new TCanvas("c5", fDiffM->GetTitle());
230 fDiffM->SetFillColor(2);
231 fDiffM->SetFillStyle(3001);
232 c->SetLogz();
233 fDiffM->Draw("colz");
234
235
8f6ee336 236 return kTRUE;
237 }
238
239 ClassDef(DrawHitsRecs,0);
240};
241
242//____________________________________________________________________
243//
244// EOF
245//