source: tags/Mars-V0.8.6/mhist/MHSigmaTheta.cc

Last change on this file was 4841, checked in by wittek, 20 years ago
*** empty log message ***
File size: 11.5 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 1/2003 <mailto:wittek@mppmu.mpg.de>
19! Author(s): Thomas Bretz 4/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2003
22!
23!
24\* ======================================================================== */
25
26//////////////////////////////////////////////////////////////////////////////
27// //
28// MHSigmaTheta //
29// //
30// calculates - the 2D-histogram sigmabar vs. Theta (Inner), and //
31// - the 2D-histogram sigmabar vs. Theta (Outer) //
32// - the 3D-histogram sigma, pixel no., Theta //
33// - the 3D-histogram (sigma^2-sigmabar^2), pixel no., Theta //
34// //
35//////////////////////////////////////////////////////////////////////////////
36
37#include "MHSigmaTheta.h"
38
39#include <TCanvas.h>
40
41#include "MTime.h"
42#include "MPointingPos.h"
43
44#include "MBinning.h"
45#include "MParList.h"
46
47#include "MGeomCam.h"
48#include "MGeomPix.h"
49#include "MBadPixelsCam.h"
50
51#include "MPedPhotCam.h"
52#include "MPedPhotPix.h"
53
54#include "MCerPhotEvt.h"
55#include "MCerPhotPix.h"
56
57#include "MLog.h"
58#include "MLogManip.h"
59
60ClassImp(MHSigmaTheta);
61
62using namespace std;
63
64// --------------------------------------------------------------------------
65//
66// Default Constructor. It sets name and title of the histogram.
67//
68MHSigmaTheta::MHSigmaTheta(const char *name, const char *title)
69{
70 fName = name ? name : "MHSigmaTheta";
71 fTitle = title ? title : "2D histogram sigmabar vs. Theta";
72
73 fSigmaTheta.SetDirectory(NULL);
74 fSigmaTheta.SetName("2D-ThetaSigmabar(Inner)");
75 fSigmaTheta.SetTitle("2D: \\bar{\\sigma}, \\Theta");
76 fSigmaTheta.SetXTitle("\\Theta [\\circ]");
77 fSigmaTheta.SetYTitle("Sigmabar(Inner) / SQRT(Area)");
78
79 fSigmaThetaOuter.SetDirectory(NULL);
80 fSigmaThetaOuter.SetName("2D-ThetaSigmabar(Outer)");
81 fSigmaThetaOuter.SetTitle("2D: \\bar{\\sigma}, \\Theta");
82 fSigmaThetaOuter.SetXTitle("\\Theta [\\circ]");
83 fSigmaThetaOuter.SetYTitle("Sigmabar(Outer) / SQRT(Area)");
84
85 fSigmaPixTheta.SetDirectory(NULL);
86 fSigmaPixTheta.SetName("3D-ThetaPixSigma");
87 fSigmaPixTheta.SetTitle("3D: \\Theta, Pixel Id, \\sigma");
88 fSigmaPixTheta.SetXTitle("\\Theta [\\circ]");
89 fSigmaPixTheta.SetYTitle("Pixel Id");
90 fSigmaPixTheta.SetZTitle("Sigma");
91
92 fDiffPixTheta.SetDirectory(NULL);
93 fDiffPixTheta.SetName("3D-ThetaPixDiff");
94 fDiffPixTheta.SetTitle("3D: \\Theta, Pixel Id, {\\sigma}^{2}-\\bar{\\sigma}^{2}");
95 fDiffPixTheta.SetXTitle("\\Theta [\\circ]");
96 fDiffPixTheta.SetYTitle("Pixel Id");
97 fDiffPixTheta.SetZTitle("(Sigma2 - Sigmabar2)/Area");
98
99 // Set default binning
100 // FIXME: Maybe ist's necessary to adapt the value to the
101 // Magic default values
102 MBinning binsb;
103 MBinning binst;
104 MBinning binsd;
105 binsd.SetEdges(100, -10, 20);
106 binsb.SetEdges(100, 0, 10);
107 binst.SetEdgesCos(10, 0, 90);
108
109 MBinning binspix("BinningPixel");
110 binspix.SetEdges(578, -0.5, 577.5);
111
112 SetBinning(&fSigmaTheta, &binst, &binsb);
113 SetBinning(&fSigmaThetaOuter, &binst, &binsb);
114 SetBinning(&fSigmaPixTheta, &binst, &binspix, &binsb);
115 SetBinning(&fDiffPixTheta, &binst, &binspix, &binsd);
116
117 fNamePedPhotCam = "MPedPhotCamFromData";
118}
119
120// --------------------------------------------------------------------------
121//
122// Set the binnings and prepare the filling of the histogram
123//
124Bool_t MHSigmaTheta::SetupFill(const MParList *plist)
125{
126 fCam = (MGeomCam*)plist->FindObject("MGeomCam");
127 if (!fCam)
128 {
129 *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl;
130 return kFALSE;
131 }
132
133 fPointPos = (MPointingPos*)plist->FindObject("MPointingPos");
134 if (!fPointPos)
135 *fLog << warn << "MPointingPos not found... aborting." << endl;
136
137
138
139 fPed = (MPedPhotCam*)plist->FindObject(AddSerialNumber(fNamePedPhotCam), "MPedPhotCam");
140 if (!fPed)
141 {
142 *fLog << err << AddSerialNumber(fNamePedPhotCam)
143 << "[MPedPhotCam] not found... aborting." << endl;
144 return kFALSE;
145 }
146 //fPed->InitSize(fCam->GetNumPixels());
147
148
149 fBad = (MBadPixelsCam*)plist->FindObject("MBadPixelsCam");
150 if (!fBad)
151 {
152 *fLog << err << "MBadPixelsCam not found... continue. " << endl;
153 }
154
155
156 fEvt = (MCerPhotEvt*)plist->FindObject("MCerPhotEvt");
157 if (!fEvt)
158 {
159 *fLog << err << "MCerPhotEvt not found... aborting." << endl;
160 return kFALSE;
161 }
162
163
164 // Get Theta Binning
165 MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta", "MBinning");
166 if (!binstheta)
167 {
168 *fLog << warn << "Object 'BinningTheta' [MBinning] not found... no binning applied." << endl;
169 return kTRUE;
170 }
171
172 // Get Sigmabar binning
173 MBinning* binssigma = (MBinning*)plist->FindObject("BinningSigmabar", "MBinning");
174 if (!binssigma)
175 *fLog << warn << "Object 'BinningSigmabar' [MBinning] not found... no binning applied." << endl;
176
177 // Get binning for (sigma^2-sigmabar^2)
178 MBinning* binsdiff = (MBinning*)plist->FindObject("BinningDiffsigma2", "MBinning");
179 if (!binsdiff)
180 *fLog << warn << "Object 'BinningDiffsigma2' [MBinning] not found... no binning applied." << endl;
181
182 //FIXME: Missing: Apply new binning to only one axis!
183
184 // Get binning for pixel number
185 const UInt_t npix1 = fPed->GetSize()+1;
186 *fLog << "npix1 = " << npix1 << endl;
187 MBinning binspix("BinningPixel");
188 binspix.SetEdges(npix1, -0.5, npix1-0.5);
189
190 // Set binnings in histograms
191 if (binssigma)
192 {
193 SetBinning(&fSigmaTheta, binstheta, binssigma);
194 SetBinning(&fSigmaThetaOuter, binstheta, binssigma);
195 SetBinning(&fSigmaPixTheta, binstheta, &binspix, binssigma);
196 }
197
198 if (binsdiff)
199 SetBinning(&fDiffPixTheta, binstheta, &binspix, binsdiff);
200
201 return kTRUE;
202}
203
204// --------------------------------------------------------------------------
205//
206// Fill the histograms
207//
208// ignore pixels if they are unused or blind
209//
210Bool_t MHSigmaTheta::Fill(const MParContainer *par, const Stat_t w)
211{
212 Double_t theta = fPointPos->GetZd();
213
214 Double_t mysig = (fPed->GetArea(0)).GetRms();
215 Double_t mysigouter = (fPed->GetArea(1)).GetRms();
216
217 *fLog << "theta, mysig, mysigouter = " << theta << ", " << mysig
218 << ", " << mysigouter << endl;
219
220 fSigmaTheta.Fill(theta, mysig);
221 fSigmaThetaOuter.Fill(theta, mysigouter);
222
223 const UInt_t npix = fEvt->GetNumPixels();
224
225 for (UInt_t i=0; i<npix; i++)
226 {
227 MCerPhotPix &cerpix = (*fEvt)[i];
228 const Int_t id = cerpix.GetPixId();
229
230 if (!cerpix.IsPixelUsed())
231 {
232 //*fLog << all << "MHSigmaTheta::Fill; unused pixel found, id = "
233 // << id << endl;
234 continue;
235 }
236
237 const MPedPhotPix &pix = (*fPed)[id];
238
239 if ( fBad != NULL && ((*fBad)[id]).IsUnsuitable() )
240 {
241 // this should never occur, because bad pixels should have
242 // been set unused
243 *fLog << all << "MHSigmaTheta::Fill; bad pixel found which is used, id = "
244 << id << "... go to next pixel." << endl;
245 continue;
246 }
247
248 // ratio is the area of pixel 0
249 // divided by the area of the current pixel
250 const Double_t ratio = fCam->GetPixRatio(id);
251 const Double_t sigma = pix.GetRms();
252
253 fSigmaPixTheta.Fill(theta, (Double_t)id, sigma*sqrt(ratio));
254
255 Double_t diff;
256 const Byte_t aidx = (*fCam)[id].GetAidx();
257 if (aidx == 0)
258 {
259 // inner pixel
260 diff = sigma*sigma*ratio - mysig*mysig;
261 }
262 else
263 {
264 // outer pixel
265 diff = sigma*sigma*ratio - mysigouter*mysigouter;
266 }
267
268 fDiffPixTheta.Fill(theta, (Double_t)id, diff);
269 }
270
271 return kTRUE;
272}
273
274// --------------------------------------------------------------------------
275//
276// Update the projections and (if possible) set log scales before painting
277//
278void MHSigmaTheta::Paint(Option_t *opt)
279{
280 TVirtualPad *padsave = gPad;
281
282 TH1D* h;
283
284 padsave->cd(1);
285 if ((h = (TH1D*)gPad->FindObject("ProjX-Theta")))
286 {
287 ProjectionX(*h, fSigmaTheta);
288 if (h->GetEntries()!=0)
289 gPad->SetLogy();
290 }
291
292 padsave->cd(4);
293 if ((h = (TH1D*)gPad->FindObject("ProjY-sigma")))
294 ProjectionY(*h, fSigmaTheta);
295
296 gPad = padsave;
297}
298
299// --------------------------------------------------------------------------
300//
301// Draw the histogram
302//
303void MHSigmaTheta::Draw(Option_t *opt)
304{
305 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
306 pad->SetBorderMode(0);
307
308 AppendPad("");
309
310 pad->Divide(3, 3);
311
312 // draw the 2D histogram Sigmabar versus Theta
313 TH1 *h;
314
315 pad->cd(1);
316 gPad->SetBorderMode(0);
317 h = fSigmaTheta.ProjectionX("ProjX-Theta", -1, 9999, "E");
318 h->SetDirectory(NULL);
319 h->SetTitle("Distribution of \\Theta");
320 h->SetXTitle("\\Theta [\\circ]");
321 h->SetYTitle("No.of events");
322 h->Draw(opt);
323 h->SetBit(kCanDelete);
324
325 pad->cd(2);
326 gPad->SetBorderMode(0);
327 h = fDiffPixTheta.Project3D("zx");
328 h->SetDirectory(NULL);
329 h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. \\Theta (all pixels)");
330 h->SetXTitle("\\Theta [\\circ]");
331 h->SetYTitle("(Sigma2 - Sigmabar2) / Area");
332 h->Draw("box");
333 h->SetBit(kCanDelete);
334
335 pad->cd(3);
336 gPad->SetBorderMode(0);
337 h = fSigmaPixTheta.Project3D("zx");
338 h->SetDirectory(NULL);
339 h->SetTitle("\\sigma_{ped} vs. \\Theta (all pixels)");
340 h->SetXTitle("\\Theta [\\circ]");
341 h->SetYTitle("Sigma");
342 h->Draw("box");
343 h->SetBit(kCanDelete);
344
345 //pad->cd(7);
346 //gPad->SetBorderMode(0);
347 //h = fSigmaTheta.ProjectionY("ProjY-sigma", -1, 9999, "E");
348 //h->SetDirectory(NULL);
349 //h->SetTitle("Distribution of \\bar{\\sigma}_{ped}");
350 //h->SetXTitle("\\bar{\\sigma}_{ped}");
351 //h->SetYTitle("No.of events");
352 //h->Draw(opt);
353 //h->SetBit(kCanDelete);
354
355 pad->cd(5);
356 gPad->SetBorderMode(0);
357 h = fDiffPixTheta.Project3D("zy");
358 h->SetDirectory(NULL);
359 h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. pixel Id (all \\Theta)");
360 h->SetXTitle("Pixel Id");
361 h->SetYTitle("Sigma2 - sigmabar2) / Area");
362 h->Draw("box");
363 h->SetBit(kCanDelete);
364
365 pad->cd(6);
366 gPad->SetBorderMode(0);
367 h = fSigmaPixTheta.Project3D("zy");
368 h->SetDirectory(NULL);
369 h->SetTitle("\\sigma_{ped} vs. pixel Id (all \\Theta)");
370 h->SetXTitle("Pixel Id");
371 h->SetYTitle("Sigma");
372 h->Draw("box");
373 h->SetBit(kCanDelete);
374
375 pad->cd(4);
376 fSigmaTheta.Draw(opt);
377
378 pad->cd(7);
379 fSigmaThetaOuter.Draw(opt);
380
381 //pad->cd(8);
382 //fDiffPixTheta.Draw(opt);
383
384 //pad->cd(9);
385 //fSigmaPixTheta.Draw(opt);
386}
387
388
389
390
391
392
393
394
395
396
397
398
399
400
Note: See TracBrowser for help on using the repository browser.