1 //____________________________________________________________________
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.
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.
12 // Use the script `Compile.C' to compile this class using ACLic.
15 #include <AliFMDHit.h>
16 #include <AliFMDDigit.h>
17 #include <AliFMDRecPoint.h>
18 #include <AliFMDInput.h>
19 #include <AliFMDEdepMap.h>
20 #include <AliFMDFloatMap.h>
25 #include <TParticle.h>
26 #include <TClonesArray.h>
32 //____________________________________________________________________
33 /** @class DrawHitsRecs
34 @brief Draw hit energy loss versus rec point mult
37 Root> Compile("DrawHitsRecs.C")
43 class DrawHitsRecs : public AliFMDInput
47 TH2D* fHitEvsRecM; // Histogram
48 TH2D* fHitEvsRecE; // Histogram
49 TH1D* fDiffE; // Histogram
50 TH2D* fHitsVsRecM; // Histogram
51 TH2D* fDiffM; // Histogram
60 //__________________________________________________________________
61 DrawHitsRecs(Bool_t primary=kFALSE,
62 Int_t n=900, Double_t emin=1e-3, Double_t emax=10,
63 Int_t m=21, Double_t mmin=-0.5, Double_t mmax=20.5)
69 if (fPrimary) AddLoad(kKinematics);
70 TArrayF eloss(MakeLogScale(n, emin, emax));
73 for (Int_t i = 1; i < m+1; i++) mults[i] = mults[i-1] + (mmax-mmin)/m;
75 fHitEvsAdc = new TH2D("hitEvsAdc", "#Delta E_{sim} vs. ADC",
76 n, emin, emax, 1025, -.5, 1024.5);
77 fHitEvsAdc->SetXTitle("#Delta E_{sim} [MeV]");
78 fHitEvsAdc->SetYTitle("ADC");
81 fHitEvsRecM = new TH2D("hitEvsRecM", "#Delta E_{sim} vs. M_{rec}",
82 eloss.fN-1, eloss.fArray, mults.fN-1, mults.fArray);
83 fHitEvsRecM->SetXTitle("#Delta E_{sim} [MeV]");
84 fHitEvsRecM->SetYTitle("M_{rec}");
87 fHitEvsRecE = new TH2D("hitEvsRecE", "#Delta E_{sim} vs. #Delta E_{rec}",
88 n, emin, emax, n, emin, emax);
89 fHitEvsRecE->SetXTitle("#Delta E_{sim} [MeV]");
90 fHitEvsRecE->SetYTitle("#Delta E_{rec} [MeV]");
94 fDiffE = new TH1D("diffE",
95 "#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}",
97 fDiffE->SetXTitle("#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}");
100 Double_t omin = mmin; // -.5;
101 Double_t omax = mmax; // 7.5;
103 fHitsVsRecM = new TH2D("hitsVsStrip", "# of Hits vs. M_{rec}",
104 o, omin, omax, m, mmin, mmax);
105 fHitsVsRecM->SetXTitle("# of Hits");
106 fHitsVsRecM->SetYTitle("M_{rec}");
107 fHitsVsRecM->Sumw2();
109 fDiffM = new TH2D("diffM", "M_{sim} - M_{rec}",
110 41, -20.5, 20.5, 70, 1.5, 5);
111 // 36, -TMath::Pi(),TMath::Pi());
112 fDiffM->SetXTitle("M_{sim} - M_{rec}");
113 fDiffM->SetYTitle("|#eta|");
115 // fDiffM->SetYTitle("Detector");
117 fHitEloss = new TH1D("hitEloss","#frac{#Delta E_{sim}}{#Delta x} (MeV/cm)",
119 fHitEloss->SetFillColor(2);
120 fHitEloss->SetFillStyle(3001);
123 fRecEloss = new TH1D("recEloss","#frac{#Delta E_{rec}}{#Delta x} (MeV/cm)",
125 fRecEloss->SetFillColor(4);
126 fRecEloss->SetFillStyle(3001);
129 //__________________________________________________________________
130 /** Begining of event
131 @param ev Event number
132 @return @c false on error */
133 Bool_t Begin(Int_t ev)
136 return AliFMDInput::Begin(ev);
138 //__________________________________________________________________
139 Bool_t ProcessHit(AliFMDHit* hit, TParticle*)
141 // Cache the energy loss
142 if (!hit) return kFALSE;
143 UShort_t det = hit->Detector();
144 Char_t rng = hit->Ring();
145 UShort_t sec = hit->Sector();
146 UShort_t str = hit->Strip();
148 AliWarning(Form("Bad strip number %d in hit", str));
152 TParticle* kine = fStack->Particle(hit->Track());
153 if (!kine) return kTRUE;
154 if (!kine->IsPrimary()) return kTRUE;
157 if (hit->Edep()/hit->Length() > 0.1)
158 fHitEloss->Fill(hit->Edep() / hit->Length());
159 fMap(det, rng, sec, str).fEdep += hit->Edep();
160 fMap(det, rng, sec, str).fN++;
164 //__________________________________________________________________
165 Bool_t ProcessDigit(AliFMDDigit* digit)
167 if (!digit) return kTRUE;
169 UShort_t det = digit->Detector();
170 Char_t rng = digit->Ring();
171 UShort_t sec = digit->Sector();
172 UShort_t str = digit->Strip();
174 AliWarning(Form("Bad strip number %d in digit", str));
177 Double_t edep = fMap(det, rng, sec, str).fEdep;
178 if (edep > 0) fHitEvsAdc->Fill(edep, digit->Counts());
183 //__________________________________________________________________
184 Bool_t ProcessRecPoint(AliFMDRecPoint* single)
186 if (!single) return kTRUE;
187 UShort_t det = single->Detector();
188 Char_t rng = single->Ring();
189 UShort_t sec = single->Sector();
190 UShort_t str = single->Strip();
192 AliWarning(Form("Bad strip number %d in single", str));
195 Double_t edep = fMap(det, rng, sec, str).fEdep;
196 Int_t nhit = fMap(det, rng, sec, str).fN;
198 fHitEvsRecM->Fill(edep, single->Particles());
199 fHitEvsRecE->Fill(edep, single->Edep());
200 fDiffE->Fill((single->Edep() - edep) / edep);
202 if (nhit > 0) fHitsVsRecM->Fill(nhit, single->Particles());
203 fDiffM->Fill(nhit - single->Particles(), TMath::Abs(single->Eta()));
204 if (single->Edep()/.03 > 0.1) fRecEloss->Fill(single->Edep() / 0.0300);
207 //__________________________________________________________________
210 gStyle->SetPalette(1);
211 gStyle->SetOptTitle(0);
212 gStyle->SetCanvasColor(0);
213 gStyle->SetCanvasBorderSize(0);
214 gStyle->SetPadColor(0);
215 gStyle->SetPadBorderSize(0);
218 c = new TCanvas("c0", fHitEvsAdc->GetTitle());
219 fHitEvsAdc->SetStats(kFALSE);
220 fHitEvsAdc->Draw("COLZ");
222 c = new TCanvas("c1", fHitEvsRecM->GetTitle());
223 fHitEvsRecM->SetStats(kFALSE);
224 fHitEvsRecM->Draw("COLZ");
226 c = new TCanvas("c2", fHitEvsRecE->GetTitle());
227 fHitEvsRecE->SetStats(kFALSE);
228 fHitEvsRecE->Draw("COLZ");
230 c = new TCanvas("c3", fDiffE->GetTitle());
234 c = new TCanvas("c4", fHitsVsRecM->GetTitle());
236 fHitsVsRecM->SetStats(kFALSE);
237 fHitsVsRecM->Draw("COLZ");
239 c = new TCanvas("c5", fDiffM->GetTitle());
240 fDiffM->SetFillColor(2);
241 fDiffM->SetFillStyle(3001);
243 fDiffM->Draw("colz");
245 c = new TCanvas("c6", "Hit Eloss, Reco Eloss");
246 fRecEloss->Scale(1./fRecEloss->GetEntries());
247 fRecEloss->Draw("hist e");
248 fRecEloss->Fit("landau", "", "SAME", 2, 4);
249 TF1* recResp = new TF1(*fRecEloss->GetFunction("landau"));
250 fHitEloss->Scale(1./fHitEloss->GetEntries());
251 fHitEloss->Draw("same hist e");
252 fHitEloss->Fit("landau", "", "SAME", 2, 10);
253 TF1* hitResp = new TF1(*fHitEloss->GetFunction("landau"));
254 std::cout << "From hits: MPV="
255 << hitResp->GetParameter(1) << "+/-"
256 << hitResp->GetParError(1) << " Width="
257 << hitResp->GetParameter(2) << "+/-"
258 << hitResp->GetParError(2) << "\nFrom recs: MPV="
259 << recResp->GetParameter(1) << "+/-"
260 << recResp->GetParError(1) << " Width="
261 << recResp->GetParameter(2) << "+/-"
262 << recResp->GetParError(2) << "\nRatio: MPV(hit/rec)="
263 << hitResp->GetParameter(1) / recResp->GetParameter(1)
270 ClassDef(DrawHitsRecs,0);
273 //____________________________________________________________________