source: branches/removing_cpp11_features/mhistmc/MHMcEnergyMigration.cc@ 18242

Last change on this file since 18242 was 2425, checked in by moralejo, 21 years ago
*** empty log message ***
File size: 9.1 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Wolfgang Wittek 4/2002 <mailto:wittek@mppmu.mpg.de>
19! Abelardo Moralejo 5/2003 <mailto:moralejo@pd.infn.it>
20!
21! Copyright: MAGIC Software Development, 2000-2002
22!
23!
24\* ======================================================================== */
25
26//////////////////////////////////////////////////////////////////////////////
27// //
28// MHMcEnergyMigration //
29// //
30// calculates the migration matrix E-est vs. E-true //
31// for different bins in Theta //
32// //
33//////////////////////////////////////////////////////////////////////////////
34
35#include "MHMcEnergyMigration.h"
36
37#include <TCanvas.h>
38
39#include "MMcEvt.hxx"
40
41#include "MEnergyEst.h"
42#include "MBinning.h"
43#include "MHillasSrc.h"
44#include "MParList.h"
45
46#include "MLog.h"
47#include "MLogManip.h"
48#include <TProfile.h>
49
50ClassImp(MHMcEnergyMigration);
51
52using namespace std;
53
54// --------------------------------------------------------------------------
55//
56// Default Constructor. It sets name and title of the histogram.
57//
58MHMcEnergyMigration::MHMcEnergyMigration(const char *name, const char *title)
59 : fHist(), fHist2(), fHistImp()
60{
61 //
62 // set the name and title of this object
63 //
64 fName = name ? name : "MHMcEnergyMigration";
65 fTitle = title ? title : "3-D histogram E-true E-est Theta";
66
67 fHist.SetDirectory(NULL);
68 fHist.SetTitle("3D-plot E_{EST} E_{TRUE} \\Theta");
69 fHist.SetXTitle("E_{EST} [GeV]");
70 fHist.SetYTitle("E_{TRUE} [GeV]");
71 fHist.SetZTitle("\\Theta [\\circ]");
72
73 fHist2.SetDirectory(NULL);
74 fHist2.SetTitle("3D-plot \\Delta E / E vs E_{EST} E_{TRUE}");
75 fHist2.SetXTitle("E_{EST} [GeV]");
76 fHist2.SetYTitle("E_{TRUE} [GeV]");
77 fHist2.SetZTitle("\\frac{E_{EST} - E_{TRUE}}{E_{TRUE}}");
78
79 fHistImp.SetDirectory(NULL);
80 fHistImp.SetTitle("\\Delta E / E vs Impact parameter");
81 fHistImp.SetXTitle("Impact parameter (m)");
82 fHistImp.SetYTitle("\\frac{E_{EST} - E_{TRUE}}{E_{TRUE}}");
83
84}
85
86// --------------------------------------------------------------------------
87//
88// Set the binnings and prepare the filling of the histograms
89//
90Bool_t MHMcEnergyMigration::SetupFill(const MParList *plist)
91{
92 fEnergy = (MEnergyEst*)plist->FindObject("MEnergyEst");
93 if (!fEnergy)
94 {
95 *fLog << err << dbginf << "MEnergyEst not found... aborting." << endl;
96 return kFALSE;
97 }
98
99 fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt");
100 if (!fMcEvt)
101 {
102 *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
103 return kFALSE;
104 }
105
106 const MBinning* binsenergy = (MBinning*)plist->FindObject("BinningE");
107 const MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta");
108 if (!binsenergy || !binstheta)
109 {
110 *fLog << err << dbginf << "At least one MBinning not found... aborting." << endl;
111 return kFALSE;
112 }
113
114 SetBinning(&fHist, binsenergy, binsenergy, binstheta);
115 fHist.Sumw2();
116
117 const MBinning* binsde = (MBinning*)plist->FindObject("BinningDE");
118 const MBinning* binsimpact = (MBinning*)plist->FindObject("BinningImpact");
119 SetBinning(&fHistImp, binsimpact, binsde);
120
121 fHistImp.Sumw2();
122
123 SetBinning(&fHist2, binsenergy, binsenergy, binsde);
124 fHist2.Sumw2();
125
126 return kTRUE;
127}
128
129// --------------------------------------------------------------------------
130//
131// Draw the histogram
132//
133void MHMcEnergyMigration::Draw(Option_t *opt)
134{
135 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
136 pad->SetBorderMode(0);
137
138 AppendPad("");
139
140 pad->Divide(2,2);
141
142 TH1 *h;
143
144 pad->cd(1);
145 gPad->SetBorderMode(0);
146 gPad->SetLogx();
147 h = fHist.Project3D("ey_pro");
148 h->SetTitle("Distribution of E_{TRUE}");
149 h->SetXTitle("E_{TRUE} [GeV]");
150 h->SetBit(kCanDelete);
151 h->Draw(opt);
152 h->SetDirectory(NULL);
153
154 pad->cd(2);
155 gPad->SetBorderMode(0);
156 gPad->SetLogx();
157 h = fHist.Project3D("ex_pro");
158 h->SetTitle("Distribution of E_{EST}");
159 h->SetXTitle("E_{est} [GeV]");
160 h->SetBit(kCanDelete);
161 h->Draw(opt);
162 Double_t minEest = h->GetXaxis()->GetXmin();
163 h->SetDirectory(NULL);
164
165 pad->cd(3);
166 gPad->SetBorderMode(0);
167 h = fHist.Project3D("z_pro");
168 h->SetTitle("Distribution of \\Theta");
169 h->SetXTitle("\\Theta [\\circ]");
170 h->SetBit(kCanDelete);
171 h->SetLineWidth(2);
172 h->Draw(opt);
173 h->SetDirectory(NULL);
174
175 pad->cd(4);
176 gPad->SetBorderMode(0);
177 TH1D* hpx;
178 hpx = fHistImp.ProjectionX("_px", 1, fHistImp.GetNbinsY(),"e");
179 hpx->SetTitle("Distribution of Impact parameter");
180 hpx->SetXTitle("Impact parameter (m)");
181 hpx->SetBit(kCanDelete);
182 hpx->Draw(opt);
183 hpx->SetDirectory(NULL);
184 fHistImp.SetDirectory(NULL);
185
186 pad->Modified();
187 pad->Update();
188
189 TVirtualPad *pad2 = MakeDefCanvas((TString)GetName()+"2");
190 pad2->SetBorderMode(0);
191
192 AppendPad("");
193
194 pad2->Divide(2,2);
195
196 TH2D *h2;
197
198 pad2->cd(1);
199 gPad->SetBorderMode(0);
200 gPad->SetLogx();
201 gPad->SetLogy();
202 h2 = (TH2D*) fHist.Project3D("yx");
203
204 TProfile* h2pfx;
205 h2pfx = h2->ProfileX("_pfx", 1, h2->GetNbinsY(),"S");
206 h2pfx->SetXTitle("E_{TRUE} (GeV)");
207 h2pfx->SetYTitle("E_{EST} (GeV)");
208 h2pfx->SetTitle("E_{EST} vs E_{TRUE}");
209 h2pfx->SetBit(kCanDelete);
210 h2pfx->SetFillColor(41);
211 h2pfx->SetFillStyle(1001);
212 h2pfx->SetMinimum(minEest);
213 h2pfx->SetStats(kFALSE);
214 h2pfx->DrawCopy("E3");
215
216 h2->SetBit(kCanDelete);
217 h2->SetFillColor(1);
218 h2->Draw("box,same");
219 h2->SetDirectory(NULL);
220
221 h2pfx->SetLineColor(2);
222 h2pfx->SetLineWidth(2);
223 h2pfx->SetFillColor(0);
224 h2pfx->DrawCopy("C,hist,same");
225 h2pfx->SetDirectory(NULL);
226
227 pad2->cd(2);
228
229 gPad->SetBorderMode(0);
230 gPad->SetLogx();
231 gPad->SetLeftMargin(0.15);
232 h2 = (TH2D*) fHist2.Project3D("zy");
233 h2->SetBit(kCanDelete);
234 h2pfx = h2->ProfileX("_pfx", 1, h2->GetNbinsY(),"S");
235 h2pfx->SetTitle("\\Delta E / E vs E_{TRUE}");
236 h2pfx->SetXTitle("E_{TRUE} (GeV)");
237 h2pfx->SetYTitle("\\frac{E_{EST} - E_{TRUE}}{E_{TRUE}}");
238 h2pfx->SetBit(kCanDelete);
239 h2pfx->SetFillColor(41);
240 h2pfx->SetFillStyle(1001);
241 h2pfx->GetYaxis()->SetTitleOffset(1.4);
242 h2pfx->SetStats(kFALSE);
243 h2pfx->DrawCopy("E3");
244 h2->SetFillColor(1);
245 h2->Draw("same,box");
246 h2->SetDirectory(NULL);
247 h2pfx->SetLineColor(2);
248 h2pfx->SetLineWidth(2);
249 h2pfx->SetFillColor(0);
250 h2pfx->DrawCopy("C,hist,same");
251 h2pfx->SetDirectory(NULL);
252
253 pad2->cd(3);
254 gPad->SetBorderMode(0);
255 gPad->SetLogx();
256 gPad->SetLeftMargin(0.15);
257 h2 = (TH2D*) fHist2.Project3D("zx");
258 h2->SetBit(kCanDelete);
259 h2pfx = h2->ProfileX("_pfx", 1, h2->GetNbinsY(),"S");
260 h2pfx->SetTitle("\\Delta E / E vs E_{EST}");
261 h2pfx->SetXTitle("E_{EST} (GeV)");
262 h2pfx->SetYTitle("\\frac{E_{EST} - E_{TRUE}}{E_{TRUE}}");
263 h2pfx->SetBit(kCanDelete);
264 h2pfx->SetFillColor(41);
265 h2pfx->SetFillStyle(1001);
266 h2pfx->GetYaxis()->SetTitleOffset(1.4);
267 h2pfx->SetStats(kFALSE);
268 h2pfx->SetMinimum(-1.);
269 h2pfx->SetMaximum(1.);
270 h2pfx->DrawCopy("E3");
271
272 h2->SetFillColor(1);
273 h2->Draw("same,box");
274 h2->SetDirectory(NULL);
275 h2pfx->SetLineColor(2);
276 h2pfx->SetLineWidth(2);
277 h2pfx->SetFillColor(0);
278 h2pfx->DrawCopy("L,histo,same");
279 h2pfx->SetDirectory(NULL);
280
281 pad2->cd(4);
282 gPad->SetBorderMode(0);
283 h = fHist2.ProjectionZ("_pz",1,fHist2.GetNbinsX(),1,fHist2.GetNbinsY(),"e");
284 h->SetBit(kCanDelete);
285 h->Draw();
286 h->SetDirectory(NULL);
287
288 pad2->Modified();
289 pad2->Update();
290
291 fHist.SetDirectory(NULL);
292 fHist2.SetDirectory(NULL);
293
294}
295
296// --------------------------------------------------------------------------
297//
298// Fill the histogram
299//
300Bool_t MHMcEnergyMigration::Fill(const MParContainer *par, const Stat_t w)
301{
302 // get E-true from fMcEvt and E-est from fEnergy
303
304 fHist.Fill(fEnergy->GetEnergy(), fMcEvt->GetEnergy(), fMcEvt->GetTelescopeTheta()*kRad2Deg, w);
305
306 fHist2.Fill(fEnergy->GetEnergy(), fMcEvt->GetEnergy(), (fEnergy->GetEnergy()-fMcEvt->GetEnergy())/fMcEvt->GetEnergy(), w);
307
308 fHistImp.Fill(fMcEvt->GetImpact()/100., (fEnergy->GetEnergy()-fMcEvt->GetEnergy())/fMcEvt->GetEnergy(), w);
309
310 return kTRUE;
311}
Note: See TracBrowser for help on using the repository browser.