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): Thomas Bretz, 12/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
|
---|
19 | !
|
---|
20 | ! Copyright: MAGIC Software Development, 2000-2003
|
---|
21 | !
|
---|
22 | !
|
---|
23 | \* ======================================================================== */
|
---|
24 |
|
---|
25 | /////////////////////////////////////////////////////////////////////////////
|
---|
26 | //
|
---|
27 | // MHCamEvent
|
---|
28 | //
|
---|
29 | // A histogram to store the sum of camera events. This can be photons,
|
---|
30 | // currents or enything else derived from MCamEvent
|
---|
31 | //
|
---|
32 | // To plot the variance instead of the rms use:
|
---|
33 | // MHCamEvent::SetBit(MHCamera::kVariance);
|
---|
34 | //
|
---|
35 | // Axis titles
|
---|
36 | // ===========
|
---|
37 | //
|
---|
38 | // 1) If no other title is given 'a.u.' is used.
|
---|
39 | // 2) If the title of MHCamEvent is different from the default,
|
---|
40 | // it is used as histogram title. You can use this to set the
|
---|
41 | // axis title, too. For more information see TH1::SetTitle, eg.
|
---|
42 | // SetTitle("MyHist;;y[cm];Counts");
|
---|
43 | // Make sure to set an empty x-axis title.
|
---|
44 | //
|
---|
45 | //
|
---|
46 | // For example:
|
---|
47 | // MHCamEvent myhist("Titele;;y [cm]");
|
---|
48 | //
|
---|
49 | //
|
---|
50 | /////////////////////////////////////////////////////////////////////////////
|
---|
51 | #include "MHCamEvent.h"
|
---|
52 |
|
---|
53 | #include <TCanvas.h>
|
---|
54 | #include <TPaveStats.h>
|
---|
55 |
|
---|
56 | #include "MLog.h"
|
---|
57 | #include "MLogManip.h"
|
---|
58 |
|
---|
59 | #include "MParList.h"
|
---|
60 | #include "MCamEvent.h"
|
---|
61 | #include "MHCamera.h"
|
---|
62 |
|
---|
63 | #include "MGeomCam.h"
|
---|
64 |
|
---|
65 | ClassImp(MHCamEvent);
|
---|
66 |
|
---|
67 | using namespace std;
|
---|
68 |
|
---|
69 | const TString MHCamEvent::gsDefName = "MHCamEvent";
|
---|
70 | const TString MHCamEvent::gsDefTitle = "Average of MCamEvents";
|
---|
71 |
|
---|
72 | // --------------------------------------------------------------------------
|
---|
73 | //
|
---|
74 | // Initialize the name and title of the task
|
---|
75 | //
|
---|
76 | void MHCamEvent::Init(const char *name, const char *title)
|
---|
77 | {
|
---|
78 | //
|
---|
79 | // set the name and title of this object
|
---|
80 | //
|
---|
81 | fName = name ? name : gsDefName.Data();
|
---|
82 | fTitle = title ? title : gsDefTitle.Data();
|
---|
83 | }
|
---|
84 |
|
---|
85 | // --------------------------------------------------------------------------
|
---|
86 | //
|
---|
87 | // Initialize the name and title of the task. Set fType to 0
|
---|
88 | //
|
---|
89 | MHCamEvent::MHCamEvent(const char *name, const char *title)
|
---|
90 | : fSum(NULL), fEvt(NULL), fType(0)
|
---|
91 | {
|
---|
92 | Init(name, title);
|
---|
93 | }
|
---|
94 |
|
---|
95 | // --------------------------------------------------------------------------
|
---|
96 | //
|
---|
97 | // Initialize the name and title of the task. Set fType to type
|
---|
98 | //
|
---|
99 | MHCamEvent::MHCamEvent(Int_t type, const char *name, const char *title)
|
---|
100 | : fSum(NULL), fEvt(NULL), fType(type)
|
---|
101 | {
|
---|
102 | Init(name, title);
|
---|
103 | }
|
---|
104 |
|
---|
105 | // --------------------------------------------------------------------------
|
---|
106 | //
|
---|
107 | // Delete the corresponding camera display if available
|
---|
108 | //
|
---|
109 | MHCamEvent::~MHCamEvent()
|
---|
110 | {
|
---|
111 | if (fSum)
|
---|
112 | delete fSum;
|
---|
113 | }
|
---|
114 |
|
---|
115 | // --------------------------------------------------------------------------
|
---|
116 | //
|
---|
117 | // Get the event (MCamEvent) the histogram might be filled with. If
|
---|
118 | // it is not given, it is assumed, that it is filled with the argument
|
---|
119 | // of the Fill function.
|
---|
120 | // Looks for the camera geometry MGeomCam and resets the sum histogram.
|
---|
121 | //
|
---|
122 | Bool_t MHCamEvent::SetupFill(const MParList *plist)
|
---|
123 | {
|
---|
124 | fEvt = (MCamEvent*)plist->FindObject(fNameEvt, "MCamEvent");
|
---|
125 | if (!fEvt)
|
---|
126 | {
|
---|
127 | if (!fNameEvt.IsNull())
|
---|
128 | {
|
---|
129 | *fLog << err << GetDescriptor() << ": No " << fNameEvt <<" [MCamEvent] available..." << endl;
|
---|
130 | return kFALSE;
|
---|
131 | }
|
---|
132 | *fLog << warn << GetDescriptor() << ": No MCamEvent available..." << endl;
|
---|
133 | }
|
---|
134 |
|
---|
135 | MGeomCam *cam = (MGeomCam*)plist->FindObject("MGeomCam");
|
---|
136 | if (!cam)
|
---|
137 | {
|
---|
138 | *fLog << err << GetDescriptor() << ": No MGeomCam found." << endl;
|
---|
139 | return kFALSE;
|
---|
140 | }
|
---|
141 |
|
---|
142 | if (fSum)
|
---|
143 | delete (fSum);
|
---|
144 |
|
---|
145 | const TString name = fNameEvt.IsNull() ? fName : fNameEvt;
|
---|
146 |
|
---|
147 | fSum = new MHCamera(*cam, name+";avg");
|
---|
148 | if (fTitle!=gsDefTitle)
|
---|
149 | fSum->SetTitle(fTitle);
|
---|
150 | if (!fTitle.Contains(";"))
|
---|
151 | fSum->SetYTitle("a.u.");
|
---|
152 | fSum->SetBit(MHCamera::kProfile);
|
---|
153 | if (TestBit(MHCamera::kVariance))
|
---|
154 | fSum->SetBit(MHCamera::kVariance);
|
---|
155 |
|
---|
156 | return kTRUE;
|
---|
157 | }
|
---|
158 |
|
---|
159 | // --------------------------------------------------------------------------
|
---|
160 | //
|
---|
161 | // Fill the histograms with data from a MCamEvent-Container.
|
---|
162 | //
|
---|
163 | Bool_t MHCamEvent::Fill(const MParContainer *par, const Stat_t w)
|
---|
164 | {
|
---|
165 | const MCamEvent *evt = par ? dynamic_cast<const MCamEvent*>(par) : fEvt;
|
---|
166 | if (!evt)
|
---|
167 | {
|
---|
168 | *fLog << err << dbginf << "Got no MCamEvent as argument of Fill()..." << endl;
|
---|
169 | return kFALSE;
|
---|
170 | }
|
---|
171 | fSum->AddCamContent(*evt, fType);
|
---|
172 | return kTRUE;
|
---|
173 | }
|
---|
174 |
|
---|
175 | // --------------------------------------------------------------------------
|
---|
176 | //
|
---|
177 | // Take the mean of the sum histogram and print all pixel indices
|
---|
178 | // which are above sum+s*rms
|
---|
179 | //
|
---|
180 | void MHCamEvent::PrintOutliers(Float_t s) const
|
---|
181 | {
|
---|
182 | const Double_t mean = fSum->GetMean();
|
---|
183 | const Double_t rms = fSum->GetRMS();
|
---|
184 |
|
---|
185 | *fLog << all << underline << GetDescriptor() << ": Mean=" << mean << ", Rms=" << rms << endl;
|
---|
186 |
|
---|
187 | for (UInt_t i=0; i<fSum->GetNumPixels(); i++)
|
---|
188 | {
|
---|
189 | if (!fSum->IsUsed(i))
|
---|
190 | continue;
|
---|
191 |
|
---|
192 | if ((*fSum)[i+1]>mean+s*rms)
|
---|
193 | *fLog << "Contents of Pixel-Index #" << i << ": " << (*fSum)[i+1] << " > " << s << "*rms" << endl;
|
---|
194 | }
|
---|
195 | }
|
---|
196 |
|
---|
197 | // --------------------------------------------------------------------------
|
---|
198 | //
|
---|
199 | // Return fSum for "sum" and fRms for "rms"
|
---|
200 | //
|
---|
201 | TH1 *MHCamEvent::GetHistByName(const TString name)
|
---|
202 | {
|
---|
203 | return fSum;
|
---|
204 | }
|
---|
205 |
|
---|
206 | void MHCamEvent::Paint(Option_t *)
|
---|
207 | {
|
---|
208 | TVirtualPad *pad = gPad;
|
---|
209 |
|
---|
210 | pad->cd(2);
|
---|
211 | if (gPad->FindObject(Form("Proj_%p", this)))
|
---|
212 | {
|
---|
213 | TH1 *h=fSum->Projection(Form("Proj_%p", this));
|
---|
214 | if (h->GetMaximum()>0)
|
---|
215 | gPad->SetLogy();
|
---|
216 | }
|
---|
217 |
|
---|
218 | pad->cd(5);
|
---|
219 | if (gPad->FindObject(Form("ProfR_%p", this)))
|
---|
220 | fSum->RadialProfile(Form("ProfR_%p", this));
|
---|
221 |
|
---|
222 | pad->cd(6);
|
---|
223 | if (gPad->FindObject(Form("ProfA_%p", this)))
|
---|
224 | fSum->AzimuthProfile(Form("ProfA_%p", this));
|
---|
225 |
|
---|
226 | pad->cd(4);
|
---|
227 | gPad->cd(1);
|
---|
228 | MHCamera *cam = (MHCamera*)gPad->FindObject(Form("Err_%p", this));
|
---|
229 | if (cam)
|
---|
230 | cam->SetCamContent(*fSum, 1);
|
---|
231 | }
|
---|
232 |
|
---|
233 | void MHCamEvent::Draw(Option_t *)
|
---|
234 | {
|
---|
235 | TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
|
---|
236 | pad->SetBorderMode(0);
|
---|
237 |
|
---|
238 | AppendPad();
|
---|
239 |
|
---|
240 | TString name = Form("%s_5", pad->GetName());
|
---|
241 | TPad *p = new TPad(name,name,6./8,0.25,0.99,0.5,pad->GetFillColor(),0,0);
|
---|
242 | p->SetNumber(5);
|
---|
243 | p->Draw();
|
---|
244 |
|
---|
245 | name = Form("%s_6", pad->GetName());
|
---|
246 | p = new TPad(name,name,6./8,0.01,0.99,0.25,pad->GetFillColor(),0,0);
|
---|
247 | p->SetNumber(6);
|
---|
248 | p->Draw();
|
---|
249 |
|
---|
250 | pad->Divide(2,2);
|
---|
251 |
|
---|
252 | pad->cd(1);
|
---|
253 | gPad->SetBorderMode(0);
|
---|
254 | gPad->SetPad(0.01, 0.5, 0.66, 0.99);
|
---|
255 | fSum->Draw("EPhist");
|
---|
256 |
|
---|
257 | pad->cd(2);
|
---|
258 | gPad->SetBorderMode(0);
|
---|
259 | gPad->SetPad(0.66, 0.5, 0.99, 0.99);
|
---|
260 | TH1 *h = fSum->Projection(Form("Proj_%p", this), 50);
|
---|
261 | h->SetTitle("Projection");
|
---|
262 | h->SetBit(kCanDelete);
|
---|
263 | h->Draw();
|
---|
264 |
|
---|
265 | pad->cd(3);
|
---|
266 | gPad->SetPad(0.01, 0.01, 3./8, 0.5);
|
---|
267 | gPad->SetBorderMode(0);
|
---|
268 | fSum->Draw();
|
---|
269 |
|
---|
270 | pad->cd(4);
|
---|
271 | gPad->SetPad(3./8, 0.01, 6./8, 0.5);
|
---|
272 | gPad->SetBorderMode(0);
|
---|
273 |
|
---|
274 | MHCamera *cam = new MHCamera(*fSum->GetGeometry());
|
---|
275 | cam->SetName(Form("Err_%p", this));
|
---|
276 | cam->SetTitle(TestBit(MHCamera::kVariance)?"Variance":"Root Mean Squared (rms)");
|
---|
277 | cam->SetYTitle(fSum->GetYaxis()->GetTitle());
|
---|
278 | cam->SetCamContent(*fSum, 1);
|
---|
279 | cam->SetBit(kCanDelete);
|
---|
280 | cam->Draw();
|
---|
281 |
|
---|
282 | pad->cd(5);
|
---|
283 | h = (TH1*)fSum->RadialProfile(Form("ProfR_%p", this), 20);
|
---|
284 | h->SetTitle("Radial Profile");
|
---|
285 | h->SetBit(kCanDelete|TH1::kNoStats);
|
---|
286 | h->Draw();
|
---|
287 |
|
---|
288 | pad->cd(6);
|
---|
289 | h = (TH1*)fSum->AzimuthProfile(Form("ProfA_%p", this), 30);
|
---|
290 | h->SetTitle("Azimuth Profile");
|
---|
291 | h->SetBit(kCanDelete|TH1::kNoStats);
|
---|
292 | h->Draw();
|
---|
293 | }
|
---|