source: trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc@ 4675

Last change on this file since 4675 was 4584, checked in by wittek, 21 years ago
*** empty log message ***
File size: 11.6 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#include "MSigmabar.h"
47
48#include "MGeomCam.h"
49#include "MBlindPixels.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
118// --------------------------------------------------------------------------
119//
120// Set the binnings and prepare the filling of the histogram
121//
122Bool_t MHSigmaTheta::SetupFill(const MParList *plist)
123{
124 fCam = (MGeomCam*)plist->FindObject("MGeomCam");
125 if (!fCam)
126 {
127 *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl;
128 return kFALSE;
129 }
130
131 fPointPos = (MPointingPos*)plist->FindObject("MPointingPos");
132 if (!fPointPos)
133 *fLog << warn << "MPointingPos not found... aborting." << endl;
134
135
136 fPed = (MPedPhotCam*)plist->FindObject("MPedPhotCam");
137 if (!fPed)
138 {
139 *fLog << err << "MPedPhotCam not found... aborting." << endl;
140 return kFALSE;
141 }
142 fPed->InitSize(fCam->GetNumPixels());
143
144
145 fBlindPix = (MBlindPixels*)plist->FindObject("MBlindPixels");
146 if (!fBlindPix)
147 {
148 *fLog << err << "MBlindPixels not found... continue. " << endl;
149 }
150
151
152 fEvt = (MCerPhotEvt*)plist->FindObject("MCerPhotEvt");
153 if (!fEvt)
154 {
155 *fLog << err << "MCerPhotEvt not found... aborting." << endl;
156 return kFALSE;
157 }
158
159 fSigmabar = (MSigmabar*)plist->FindObject("MSigmabar");
160 if (!fSigmabar)
161 {
162 *fLog << err << "MSigmabar not found... aborting." << endl;
163 return kFALSE;
164 }
165
166 // Get Theta Binning
167 MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta", "MBinning");
168 if (!binstheta)
169 {
170 *fLog << warn << "Object 'BinningTheta' [MBinning] not found... no binning applied." << endl;
171 return kTRUE;
172 }
173
174 // Get Sigmabar binning
175 MBinning* binssigma = (MBinning*)plist->FindObject("BinningSigmabar", "MBinning");
176 if (!binssigma)
177 *fLog << warn << "Object 'BinningSigmabar' [MBinning] not found... no binning applied." << endl;
178
179 // Get binning for (sigma^2-sigmabar^2)
180 MBinning* binsdiff = (MBinning*)plist->FindObject("BinningDiffsigma2", "MBinning");
181 if (!binsdiff)
182 *fLog << warn << "Object 'BinningDiffsigma2' [MBinning] not found... no binning applied." << endl;
183
184 //FIXME: Missing: Apply new binning to only one axis!
185
186 // Get binning for pixel number
187 const UInt_t npix1 = fPed->GetSize()+1;
188 *fLog << "npix1 = " << npix1 << endl;
189 MBinning binspix("BinningPixel");
190 binspix.SetEdges(npix1, -0.5, npix1-0.5);
191
192 // Set binnings in histograms
193 if (binssigma)
194 {
195 SetBinning(&fSigmaTheta, binstheta, binssigma);
196 SetBinning(&fSigmaThetaOuter, binstheta, binssigma);
197 SetBinning(&fSigmaPixTheta, binstheta, &binspix, binssigma);
198 }
199
200 if (binsdiff)
201 SetBinning(&fDiffPixTheta, binstheta, &binspix, binsdiff);
202
203 return kTRUE;
204}
205
206// --------------------------------------------------------------------------
207//
208// Fill the histograms
209//
210// ignore pixels if they are unused or blind
211//
212Bool_t MHSigmaTheta::Fill(const MParContainer *par, const Stat_t w)
213{
214 Double_t theta = fPointPos->GetZd();
215 fSigmabar->Calc(*fCam, *fPed, *fEvt);
216 Double_t mysig = fSigmabar->GetSigmabarInner();
217 Double_t mysigouter = fSigmabar->GetSigmabarOuter();
218
219 //*fLog << "theta, mysig, mysigouter = " << theta << ", " << mysig
220 // << ", " << mysigouter << endl;
221
222 fSigmaTheta.Fill(theta, mysig);
223 fSigmaThetaOuter.Fill(theta, mysigouter);
224
225 const UInt_t npix = fEvt->GetNumPixels();
226
227 for (UInt_t i=0; i<npix; i++)
228 {
229 MCerPhotPix &cerpix = (*fEvt)[i];
230 const Int_t id = cerpix.GetPixId();
231
232 if (!cerpix.IsPixelUsed())
233 {
234 //*fLog << all << "MHSigmaTheta::Fill; unused pixel found, id = "
235 // << id << endl;
236 continue;
237 }
238
239 const MPedPhotPix &pix = (*fPed)[id];
240
241 if ( fBlindPix != NULL && fBlindPix->IsBlind(id) )
242 {
243 // this should never occur, because blind pixels should have
244 // been set unused by MBlindPixelsCalc2::UnMap()
245 //*fLog << all << "MHSigmaTheta::Fill; blind pixel found which is used, id = "
246 // << id << "... go to next pixel." << endl;
247 continue;
248 }
249
250 // ratio is the area of pixel 0
251 // divided by the area of the current pixel
252 const Double_t ratio = fCam->GetPixRatio(id);
253 const Double_t sigma = pix.GetRms();
254
255 fSigmaPixTheta.Fill(theta, (Double_t)id, sigma*sqrt(ratio));
256
257 Double_t diff;
258 if (ratio > 0.5)
259 {
260 // inner pixel
261 diff = sigma*sigma*ratio - mysig*mysig;
262 }
263 else
264 {
265 // outer pixel
266 diff = sigma*sigma*ratio - mysigouter*mysigouter;
267 }
268
269 fDiffPixTheta.Fill(theta, (Double_t)id, diff);
270 }
271
272 return kTRUE;
273}
274
275// --------------------------------------------------------------------------
276//
277// Update the projections and (if possible) set log scales before painting
278//
279void MHSigmaTheta::Paint(Option_t *opt)
280{
281 TVirtualPad *padsave = gPad;
282
283 TH1D* h;
284
285 padsave->cd(1);
286 if ((h = (TH1D*)gPad->FindObject("ProjX-Theta")))
287 {
288 ProjectionX(*h, fSigmaTheta);
289 if (h->GetEntries()!=0)
290 gPad->SetLogy();
291 }
292
293 padsave->cd(4);
294 if ((h = (TH1D*)gPad->FindObject("ProjY-sigma")))
295 ProjectionY(*h, fSigmaTheta);
296
297 gPad = padsave;
298}
299
300// --------------------------------------------------------------------------
301//
302// Draw the histogram
303//
304void MHSigmaTheta::Draw(Option_t *opt)
305{
306 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
307 pad->SetBorderMode(0);
308
309 AppendPad("");
310
311 pad->Divide(3, 3);
312
313 // draw the 2D histogram Sigmabar versus Theta
314 TH1 *h;
315
316 pad->cd(1);
317 gPad->SetBorderMode(0);
318 h = fSigmaTheta.ProjectionX("ProjX-Theta", -1, 9999, "E");
319 h->SetDirectory(NULL);
320 h->SetTitle("Distribution of \\Theta");
321 h->SetXTitle("\\Theta [\\circ]");
322 h->SetYTitle("No.of events");
323 h->Draw(opt);
324 h->SetBit(kCanDelete);
325
326 pad->cd(2);
327 gPad->SetBorderMode(0);
328 h = fDiffPixTheta.Project3D("zx");
329 h->SetDirectory(NULL);
330 h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. \\Theta (all pixels)");
331 h->SetXTitle("\\Theta [\\circ]");
332 h->SetYTitle("(Sigma2 - Sigmabar2) / Area");
333 h->Draw("box");
334 h->SetBit(kCanDelete);
335
336 pad->cd(3);
337 gPad->SetBorderMode(0);
338 h = fSigmaPixTheta.Project3D("zx");
339 h->SetDirectory(NULL);
340 h->SetTitle("\\sigma_{ped} vs. \\Theta (all pixels)");
341 h->SetXTitle("\\Theta [\\circ]");
342 h->SetYTitle("Sigma");
343 h->Draw("box");
344 h->SetBit(kCanDelete);
345
346 //pad->cd(7);
347 //gPad->SetBorderMode(0);
348 //h = fSigmaTheta.ProjectionY("ProjY-sigma", -1, 9999, "E");
349 //h->SetDirectory(NULL);
350 //h->SetTitle("Distribution of \\bar{\\sigma}_{ped}");
351 //h->SetXTitle("\\bar{\\sigma}_{ped}");
352 //h->SetYTitle("No.of events");
353 //h->Draw(opt);
354 //h->SetBit(kCanDelete);
355
356 pad->cd(5);
357 gPad->SetBorderMode(0);
358 h = fDiffPixTheta.Project3D("zy");
359 h->SetDirectory(NULL);
360 h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. pixel Id (all \\Theta)");
361 h->SetXTitle("Pixel Id");
362 h->SetYTitle("Sigma2 - sigmabar2) / Area");
363 h->Draw("box");
364 h->SetBit(kCanDelete);
365
366 pad->cd(6);
367 gPad->SetBorderMode(0);
368 h = fSigmaPixTheta.Project3D("zy");
369 h->SetDirectory(NULL);
370 h->SetTitle("\\sigma_{ped} vs. pixel Id (all \\Theta)");
371 h->SetXTitle("Pixel Id");
372 h->SetYTitle("Sigma");
373 h->Draw("box");
374 h->SetBit(kCanDelete);
375
376 pad->cd(4);
377 fSigmaTheta.Draw(opt);
378
379 pad->cd(7);
380 fSigmaThetaOuter.Draw(opt);
381
382 //pad->cd(8);
383 //fDiffPixTheta.Draw(opt);
384
385 //pad->cd(9);
386 //fSigmaPixTheta.Draw(opt);
387}
388
389
390
391
392
393
394
395
396
397
398
399
400
401
Note: See TracBrowser for help on using the repository browser.