Bug fix for reading raw data - thanks Frederic YERMIA, Diego, and Gines MARTINEZ
[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>
2b893216 30#include <TF1.h>
8f6ee336 31
bf000c32 32//____________________________________________________________________
9b48326f 33/** @class DrawHitsRecs
34 @brief Draw hit energy loss versus rec point mult
35 @code
36 Root> .L Compile.C
37 Root> Compile("DrawHitsRecs.C")
38 Root> DrawHitsRecs c
39 Root> c.Run();
40 @endcode
41 @ingroup FMD_script
42 */
15b17c89 43class DrawHitsRecs : public AliFMDInput
8f6ee336 44{
45private:
15b17c89 46 TH2D* fHitEvsAdc;
47 TH2D* fHitEvsRecM; // Histogram
48 TH2D* fHitEvsRecE; // Histogram
49 TH1D* fDiffE; // Histogram
50 TH2D* fHitsVsRecM; // Histogram
f95a63c4 51 TH2D* fDiffM; // Histogram
2b893216 52 TH1* fHitEloss;
53 TH1* fRecEloss;
15b17c89 54 AliFMDEdepMap fMap;
8f6ee336 55 AliFMDFloatMap fEta;
56 AliFMDFloatMap fPhi;
57 AliFMDFloatMap fMult;
58 Bool_t fPrimary;
59public:
bf000c32 60 //__________________________________________________________________
8f6ee336 61 TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max)
62 {
63 TArrayF bins(n+1);
64 Float_t dp = n / TMath::Log10(max / min);
65 Float_t pmin = TMath::Log10(min);
66 bins[0] = min;
67 for (Int_t i = 1; i < n+1; i++) {
68 Float_t p = pmin + i / dp;
69 bins[i] = TMath::Power(10, p);
70 }
71 return bins;
72 }
bf000c32 73 //__________________________________________________________________
8f6ee336 74 DrawHitsRecs(Bool_t primary=kFALSE,
75 Int_t n=900, Double_t emin=1e-3, Double_t emax=10,
76 Int_t m=21, Double_t mmin=-0.5, Double_t mmax=20.5)
77 {
78 fPrimary = primary;
79 AddLoad(kRecPoints);
8ec606c2 80 AddLoad(kHits);
15b17c89 81 AddLoad(kDigits);
8f6ee336 82 if (fPrimary) AddLoad(kKinematics);
83 TArrayF eloss(MakeLogScale(n, emin, emax));
84 TArrayF mults(m+1);
85 mults[0] = mmin;
86 for (Int_t i = 1; i < m+1; i++) mults[i] = mults[i-1] + (mmax-mmin)/m;
15b17c89 87
88 fHitEvsAdc = new TH2D("hitEvsAdc", "#Delta E_{sim} vs. ADC",
89 n, emin, emax, 1025, -.5, 1024.5);
90 fHitEvsAdc->SetXTitle("#Delta E_{sim} [MeV]");
91 fHitEvsAdc->SetYTitle("ADC");
92
93 fHitEvsRecM = new TH2D("hitEvsRecM", "#Delta E_{sim} vs. M_{rec}",
8f6ee336 94 eloss.fN-1, eloss.fArray, mults.fN-1, mults.fArray);
15b17c89 95 fHitEvsRecM->SetXTitle("#Delta E_{sim} [MeV]");
96 fHitEvsRecM->SetYTitle("M_{rec}");
97
98 fHitEvsRecE = new TH2D("hitEvsRecE", "#Delta E_{sim} vs. #Delta E_{rec}",
99 n, emin, emax, n, emin, emax);
100 fHitEvsRecE->SetXTitle("#Delta E_{sim} [MeV]");
101 fHitEvsRecE->SetYTitle("#Delta E_{rec} [MeV]");
8f6ee336 102
15b17c89 103
104 fDiffE = new TH1D("diffE",
105 "#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}",
106 1100, -1, 1.1);
107 fDiffE->SetXTitle("#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}");
108
8f6ee336 109 Double_t omin = -.5;
110 Double_t omax = 7.5;
111 Int_t o = 8;
15b17c89 112 fHitsVsRecM = new TH2D("hitsVsStrip", "# of Hits vs. M_{rec}",
8f6ee336 113 o, omin, omax, m, mmin, mmax);
15b17c89 114 fHitsVsRecM->SetXTitle("# of Hits");
115 fHitsVsRecM->SetYTitle("M_{rec}");
f95a63c4 116
117 fDiffM = new TH2D("diffM", "M_{sim} - M_{rec}",
118 41, -20.5, 20.5, 70, 1.5, 5);
119 // 36, -TMath::Pi(),TMath::Pi());
120 fDiffM->SetXTitle("M_{sim} - M_{rec}");
121 fDiffM->SetYTitle("|#eta|");
122 // fDiffM->SetYTitle("Detector");
2b893216 123
124 fHitEloss = new TH1D("hitEloss", "#frac{#Delta E_{sim}}{#Delta x} (MeV/cm)",
125 100, 0, 10);
126 fHitEloss->SetFillColor(2);
127 fHitEloss->SetFillStyle(3001);
f95a63c4 128
2b893216 129 fRecEloss = new TH1D("recEloss", "#frac{#Delta E_{rec}}{#Delta x} (MeV/cm)",
130 100, 0, 10);
131 fRecEloss->SetFillColor(4);
132 fRecEloss->SetFillStyle(3001);
8f6ee336 133 }
bf000c32 134 //__________________________________________________________________
9b48326f 135 /** Begining of event
136 @param ev Event number
137 @return @c false on error */
8f6ee336 138 Bool_t Begin(Int_t ev)
139 {
140 fMap.Reset();
15b17c89 141 return AliFMDInput::Begin(ev);
8f6ee336 142 }
bf000c32 143 //__________________________________________________________________
8f6ee336 144 Bool_t ProcessHit(AliFMDHit* hit, TParticle*)
145 {
146 // Cache the energy loss
147 if (!hit) return kFALSE;
148 UShort_t det = hit->Detector();
149 Char_t rng = hit->Ring();
150 UShort_t sec = hit->Sector();
151 UShort_t str = hit->Strip();
152 if (str > 511) {
153 AliWarning(Form("Bad strip number %d in hit", str));
154 return kTRUE;
155 }
156 if (fPrimary) {
157 TParticle* kine = fStack->Particle(hit->Track());
158 if (!kine) return kTRUE;
159 if (!kine->IsPrimary()) return kTRUE;
160 }
161
2b893216 162 if (hit->Edep()/hit->Length() > 0.1) fHitEloss->Fill(hit->Edep() / hit->Length());
8f6ee336 163 fMap(det, rng, sec, str).fEdep += hit->Edep();
164 fMap(det, rng, sec, str).fN++;
165 return kTRUE;
166 }
15b17c89 167
168 //__________________________________________________________________
169 Bool_t ProcessDigit(AliFMDDigit* digit)
170 {
171 if (!digit) return kTRUE;
172
173 UShort_t det = digit->Detector();
174 Char_t rng = digit->Ring();
175 UShort_t sec = digit->Sector();
176 UShort_t str = digit->Strip();
177 if (str > 511) {
178 AliWarning(Form("Bad strip number %d in digit", str));
179 return kTRUE;
180 }
181 Double_t edep = fMap(det, rng, sec, str).fEdep;
182 if (edep > 0) fHitEvsAdc->Fill(edep, digit->Counts());
183
184 return kTRUE;
185 }
186
bf000c32 187 //__________________________________________________________________
188 Bool_t ProcessRecPoint(AliFMDRecPoint* single)
8f6ee336 189 {
15b17c89 190 if (!single) return kTRUE;
bf000c32 191 UShort_t det = single->Detector();
192 Char_t rng = single->Ring();
193 UShort_t sec = single->Sector();
194 UShort_t str = single->Strip();
195 if (str > 511) {
196 AliWarning(Form("Bad strip number %d in single", str));
15b17c89 197 return kTRUE;
198 }
199 Double_t edep = fMap(det, rng, sec, str).fEdep;
200 Int_t nhit = fMap(det, rng, sec, str).fN;
201 if (edep > 0) {
202 fHitEvsRecM->Fill(edep, single->Particles());
203 fHitEvsRecE->Fill(edep, single->Edep());
204 fDiffE->Fill((single->Edep() - edep) / edep);
8f6ee336 205 }
15b17c89 206 if (nhit > 0) fHitsVsRecM->Fill(nhit, single->Particles());
f95a63c4 207 fDiffM->Fill(nhit - single->Particles(), TMath::Abs(single->Eta()));
2b893216 208 if (single->Edep()/.03 > 0.1) fRecEloss->Fill(single->Edep() / 0.0300);
8f6ee336 209 return kTRUE;
210 }
bf000c32 211 //__________________________________________________________________
8f6ee336 212 Bool_t Finish()
213 {
214 gStyle->SetPalette(1);
215 gStyle->SetOptTitle(0);
216 gStyle->SetCanvasColor(0);
217 gStyle->SetCanvasBorderSize(0);
218 gStyle->SetPadColor(0);
219 gStyle->SetPadBorderSize(0);
f95a63c4 220 TCanvas* c = 0;
8f6ee336 221
f95a63c4 222 c = new TCanvas("c0", fHitEvsAdc->GetTitle());
15b17c89 223 fHitEvsAdc->SetStats(kFALSE);
224 fHitEvsAdc->Draw("COLZ");
225
f95a63c4 226 c = new TCanvas("c1", fHitEvsRecM->GetTitle());
15b17c89 227 fHitEvsRecM->SetStats(kFALSE);
228 fHitEvsRecM->Draw("COLZ");
229
f95a63c4 230 c = new TCanvas("c2", fHitEvsRecE->GetTitle());
15b17c89 231 fHitEvsRecE->SetStats(kFALSE);
232 fHitEvsRecE->Draw("COLZ");
233
f95a63c4 234 c = new TCanvas("c3", fDiffE->GetTitle());
235 c->SetLogz();
15b17c89 236 fDiffE->Draw();
8f6ee336 237
f95a63c4 238 c = new TCanvas("c4", fHitsVsRecM->GetTitle());
239 c->SetLogz();
15b17c89 240 fHitsVsRecM->SetStats(kFALSE);
241 fHitsVsRecM->Draw("COLZ");
8f6ee336 242
f95a63c4 243 c = new TCanvas("c5", fDiffM->GetTitle());
244 fDiffM->SetFillColor(2);
245 fDiffM->SetFillStyle(3001);
246 c->SetLogz();
247 fDiffM->Draw("colz");
248
2b893216 249 c = new TCanvas("c6", "Hit Eloss, Reco Eloss");
250 fRecEloss->Scale(1./fRecEloss->GetMaximum());
251 fRecEloss->Draw();
252 fRecEloss->Fit("landau", "", "SAME", 2, 4);
253 TF1* recResp = new TF1(*fRecEloss->GetFunction("landau"));
254 fHitEloss->Scale(1./fHitEloss->GetMaximum());
255 fHitEloss->Draw("same");
256 fHitEloss->Fit("landau", "", "SAME", 2, 10);
257 TF1* hitResp = new TF1(*fHitEloss->GetFunction("landau"));
258 std::cout << "Hit MPV,width: " << hitResp->GetParameter(1) << ","
259 << hitResp->GetParameter(2) << "\n"
260 << "Rec MPV,width: " << recResp->GetParameter(1) << ","
261 << recResp->GetParameter(2) << std::endl;
262 c->SetLogy();
f95a63c4 263
8f6ee336 264 return kTRUE;
265 }
266
267 ClassDef(DrawHitsRecs,0);
268};
269
270//____________________________________________________________________
271//
272// EOF
273//