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

Last change on this file since 5457 was 5435, checked in by wittek, 20 years ago
*** empty log message ***
File size: 15.2 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// - the 2D-histogram Theta vs. Phi //
35// //
36//////////////////////////////////////////////////////////////////////////////
37
38#include "MHSigmaTheta.h"
39
40#include <TCanvas.h>
41
42#include "MTime.h"
43#include "MPointingPos.h"
44
45#include "MBinning.h"
46#include "MParList.h"
47
48#include "MGeomCam.h"
49#include "MGeomPix.h"
50#include "MBadPixelsCam.h"
51
52#include "MPedPhotCam.h"
53#include "MPedPhotPix.h"
54
55#include "MCerPhotEvt.h"
56#include "MCerPhotPix.h"
57
58#include "MLog.h"
59#include "MLogManip.h"
60
61ClassImp(MHSigmaTheta);
62
63using namespace std;
64
65// --------------------------------------------------------------------------
66//
67// Constructor.
68//
69MHSigmaTheta::MHSigmaTheta(const char *name, const char *title)
70{
71 fName = name ? name : "MHSigmaTheta";
72 fTitle = title ? title : "2D histogram sigmabar vs. Theta";
73
74 fThetaPhi = new TH2D;
75 fThetaPhi->SetDirectory(NULL);
76 fThetaPhi->UseCurrentStyle();
77 fThetaPhi->SetName("2D-ThetaPhi");
78 fThetaPhi->SetTitle("2D: \\Theta vs. \\Phi");
79 fThetaPhi->SetXTitle("\\phi [\\circ]");
80 fThetaPhi->SetYTitle("\\Theta [\\circ]");
81 fThetaPhi->SetTitleOffset(1.2,"Y");
82
83 fSigmaTheta = new TH2D;
84 fSigmaTheta->SetDirectory(NULL);
85 fSigmaTheta->UseCurrentStyle();
86 fSigmaTheta->SetName("2D-ThetaSigmabar(Inner)");
87 fSigmaTheta->SetTitle("2D: \\bar{\\sigma}, \\Theta");
88 fSigmaTheta->SetXTitle("\\Theta [\\circ]");
89 fSigmaTheta->SetYTitle("Sigmabar(Inner)");
90 fSigmaTheta->SetTitleOffset(1.2,"Y");
91
92 fSigmaThetaOuter = new TH2D;
93 fSigmaThetaOuter->SetDirectory(NULL);
94 fSigmaThetaOuter->UseCurrentStyle();
95 fSigmaThetaOuter->SetName("2D-ThetaSigmabar(Outer)");
96 fSigmaThetaOuter->SetTitle("2D: \\bar{\\sigma}, \\Theta");
97 fSigmaThetaOuter->SetXTitle("\\Theta [\\circ]");
98 fSigmaThetaOuter->SetYTitle("Sigmabar(Outer)");
99 fSigmaThetaOuter->SetTitleOffset(1.2,"Y");
100
101 fSigmaPixTheta = new TH3D;
102 fSigmaPixTheta->SetDirectory(NULL);
103 fSigmaPixTheta->UseCurrentStyle();
104 fSigmaPixTheta->SetName("3D-ThetaPixSigma");
105 fSigmaPixTheta->SetTitle("3D: \\Theta, Pixel Id, \\sigma");
106 fSigmaPixTheta->SetXTitle("\\Theta [\\circ]");
107 fSigmaPixTheta->SetYTitle("Pixel Id");
108 fSigmaPixTheta->SetZTitle("Sigma");
109
110 fDiffPixTheta = new TH3D;
111 fDiffPixTheta->SetDirectory(NULL);
112 fDiffPixTheta->UseCurrentStyle();
113 fDiffPixTheta->SetName("3D-ThetaPixDiff");
114 //fDiffPixTheta->SetTitle("3D: \\Theta, Pixel Id, {\\sigma}^{2}-\\bar{\\sigma}^{2}");
115 fDiffPixTheta->SetTitle("3D: \\Theta, Pixel Id, (Sigma2-Sigmabar2)/Area");
116 fDiffPixTheta->SetXTitle("\\Theta [\\circ]");
117 fDiffPixTheta->SetYTitle("Pixel Id");
118 fDiffPixTheta->SetZTitle("(Sigma2 - Sigmabar2)/Area");
119
120 // Define default binning
121 fBinsPhi = new MBinning;
122 fBinsPhi->SetEdges( 20, 0.0, 360.0);
123
124 Double_t fThetaLo = 0.0;
125 Double_t fThetaHi = 90.0;
126 fBinsTheta = new MBinning;
127 fBinsTheta->SetEdgesCos( 10, fThetaLo, fThetaHi); // theta
128 //fBinsTheta->SetEdges( 1, fThetaLo, fThetaHi); // theta
129
130 fBinsSigma = new MBinning;
131 fBinsSigma->SetEdges( 100, 0.0, 120.0); // sigma
132
133 fBinsSigmabarIn = new MBinning;
134 fBinsSigmabarOut = new MBinning;;
135 fBinsSigmabarIn->SetEdges( 100, 0.0, 25.0); // sigmabar (inner)
136 fBinsSigmabarOut->SetEdges(100, 0.0, 60.0); // sigmabar (outer)
137
138 fBinsPix = new MBinning;
139 fBinsPix->SetEdges(578, -0.5, 577.5);
140
141 fBinsDiff = new MBinning;
142 fBinsDiff->SetEdges( 100, -500.0, 1500.0); // (sigma2-sigmabar2)/area
143
144 SetBinning(fThetaPhi, fBinsPhi, fBinsTheta);
145 SetBinning(fSigmaTheta, fBinsTheta, fBinsSigmabarIn);
146 SetBinning(fSigmaThetaOuter, fBinsTheta, fBinsSigmabarOut);
147 SetBinning(fSigmaPixTheta, fBinsTheta, fBinsPix, fBinsSigma);
148 SetBinning(fDiffPixTheta, fBinsTheta, fBinsPix, fBinsDiff);
149
150 //--------------------------------------------
151 fNamePedPhotCam = "MPedPhotCamFromData";
152}
153
154
155// --------------------------------------------------------------------------
156//
157// Destructor.
158//
159MHSigmaTheta::~MHSigmaTheta()
160{
161 delete fThetaPhi;
162 delete fSigmaTheta;
163 delete fSigmaThetaOuter;
164 delete fSigmaPixTheta;
165 delete fDiffPixTheta;
166
167 delete fBinsPhi;
168 delete fBinsTheta;
169 delete fBinsSigma;
170 delete fBinsSigmabarIn;
171 delete fBinsSigmabarOut;
172 delete fBinsPix;
173 delete fBinsDiff;
174}
175
176// --------------------------------------------------------------------------
177//
178// Set the binnings and prepare the filling of the histogram
179//
180Bool_t MHSigmaTheta::SetupFill(const MParList *plist)
181{
182 fCam = (MGeomCam*)plist->FindObject("MGeomCam");
183 if (!fCam)
184 {
185 *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl;
186 return kFALSE;
187 }
188
189 fPointPos = (MPointingPos*)plist->FindObject("MPointingPos");
190 if (!fPointPos)
191 {
192 *fLog << err << "MPointingPos not found... aborting." << endl;
193 return kFALSE;
194 }
195
196 fPed = (MPedPhotCam*)plist->FindObject(AddSerialNumber(fNamePedPhotCam), "MPedPhotCam");
197 if (!fPed)
198 {
199 *fLog << err << AddSerialNumber(fNamePedPhotCam)
200 << "[MPedPhotCam] not found... aborting." << endl;
201 return kFALSE;
202 }
203 fPed->InitSize(fCam->GetNumPixels());
204
205
206 fBad = (MBadPixelsCam*)plist->FindObject("MBadPixelsCam");
207 if (!fBad)
208 {
209 *fLog << err << "MBadPixelsCam not found... continue. " << endl;
210 }
211
212
213 fEvt = (MCerPhotEvt*)plist->FindObject("MCerPhotEvt");
214 if (!fEvt)
215 {
216 *fLog << err << "MCerPhotEvt not found... aborting." << endl;
217 return kFALSE;
218 }
219
220
221 //---------------------------------------------------------------
222 *fLog << inf << "Name of MPedPhotCam container : "
223 << fNamePedPhotCam << endl;
224
225
226 //---------------------------------------------------------------
227
228 // Get Phi Binning
229 MBinning* binsphi = (MBinning*)plist->FindObject("BinningPhi", "MBinning");
230 if (!binsphi)
231 {
232 *fLog << warn << "Object 'BinningPhi' [MBinning] not found... use default binning." << endl;
233 binsphi = fBinsPhi;
234 }
235
236 // Get Theta Binning
237 MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta", "MBinning");
238 if (!binstheta)
239 {
240 *fLog << warn << "Object 'BinningTheta' [MBinning] not found... use default binning." << endl;
241 binstheta = fBinsTheta;
242 }
243
244 // Get Sigma binning
245 MBinning* binssig = (MBinning*)plist->FindObject("BinningSigma", "MBinning");
246 if (!binssig)
247 {
248 *fLog << warn << "Object 'BinningSigma' [MBinning] not found... use default binning." << endl;
249 binssig = fBinsSigma;
250 }
251
252 // Get SigmabarIn binning
253 MBinning* binssigmain = (MBinning*)plist->FindObject("BinningSigmabarIn", "MBinning");
254 if (!binssigmain)
255 {
256 *fLog << warn << "Object 'BinningSigmabarIn' [MBinning] not found... use default binning." << endl;
257 binssigmain = fBinsSigmabarIn;
258 }
259
260 // Get SigmabarOut binning
261 MBinning* binssigmaout = (MBinning*)plist->FindObject("BinningSigmabarOut", "MBinning");
262 if (!binssigmaout)
263 {
264 *fLog << warn << "Object 'BinningSigmabarOut' [MBinning] not found... use default binning." << endl;
265 binssigmaout = fBinsSigmabarOut;
266 }
267
268 // Get binning for (sigma^2-sigmabar^2)
269 MBinning* binsdiff = (MBinning*)plist->FindObject("BinningDiffsigma2", "MBinning");
270 if (!binsdiff)
271 {
272 *fLog << warn << "Object 'BinningDiffsigma2' [MBinning] not found... use default binning." << endl;
273 binsdiff = fBinsDiff;
274 }
275
276
277 //---------------------------------------------------------------
278
279 // Get binning for pixel number
280 const UInt_t npix1 = fPed->GetSize()+1;
281 //*fLog << "MHSigmaTheata::SetupFill(); npix1 = " << npix1 << endl;
282 MBinning binspix("BinningPixel");
283 binspix.SetEdges(npix1, -0.5, npix1-0.5);
284
285 // Set binnings in histograms
286 SetBinning(fThetaPhi, binsphi, binstheta);
287 SetBinning(fSigmaTheta, binstheta, binssigmain);
288 SetBinning(fSigmaThetaOuter, binstheta, binssigmaout);
289 SetBinning(fSigmaPixTheta, binstheta, &binspix, binssig);
290 SetBinning(fDiffPixTheta, binstheta, &binspix, binsdiff);
291
292 *fLog << "MHSigmaTheta::SetupFill(); binnings were set" << endl;
293
294 return kTRUE;
295}
296
297// --------------------------------------------------------------------------
298//
299// Fill the histograms
300//
301// ignore pixels if they are unused or blind
302//
303Bool_t MHSigmaTheta::Fill(const MParContainer *par, const Stat_t w)
304{
305 Double_t theta = fPointPos->GetZd();
306 Double_t phi = fPointPos->GetAz();
307
308 fPed->ReCalc(*fCam, fBad);
309
310 Double_t mysig = (fPed->GetArea(0)).GetRms();
311 Double_t mysigouter = (fPed->GetArea(1)).GetRms();
312
313 //*fLog << "theta, mysig, mysigouter = " << theta << ", " << mysig
314 // << ", " << mysigouter << endl;
315
316 fThetaPhi->Fill(phi, theta);
317 fSigmaTheta->Fill(theta, mysig);
318 fSigmaThetaOuter->Fill(theta, mysigouter);
319
320 MCerPhotPix *cerpix = 0;
321 TIter Next(*fEvt);
322 while ( (cerpix=(MCerPhotPix*)Next()) )
323 {
324 const Int_t id = cerpix->GetPixId();
325
326 if (!cerpix->IsPixelUsed())
327 {
328 //*fLog << all << "MHSigmaTheta::Fill; unused pixel found, id = "
329 // << id << endl;
330 continue;
331 }
332
333 const MPedPhotPix &pix = (*fPed)[id];
334
335 // ratio is the area of pixel 0
336 // divided by the area of the current pixel
337 const Double_t ratio = fCam->GetPixRatio(id);
338 const Double_t sigma = pix.GetRms();
339
340 fSigmaPixTheta->Fill(theta, (Double_t)id, sigma);
341
342 Double_t diff;
343 const Byte_t aidx = (*fCam)[id].GetAidx();
344 if (aidx == 0)
345 {
346 // inner pixel
347 diff = (sigma*sigma - mysig*mysig) * ratio;
348 }
349 else
350 {
351 // outer pixel
352 diff = (sigma*sigma - mysigouter*mysigouter) * ratio;
353 }
354
355 fDiffPixTheta->Fill(theta, (Double_t)id, diff);
356 }
357
358 return kTRUE;
359}
360
361// --------------------------------------------------------------------------
362//
363// Update the projections and (if possible) set log scales before painting
364//
365void MHSigmaTheta::Paint(Option_t *opt)
366{
367 TVirtualPad *padsave = gPad;
368
369 TH1D* h;
370
371 padsave->cd(1);
372 if ((h = (TH1D*)gPad->FindObject("ProjX-Theta")))
373 {
374 ProjectionX(*h, *fSigmaTheta);
375 if (h->GetEntries()!=0)
376 gPad->SetLogy();
377 }
378
379 padsave->cd(4);
380 if ((h = (TH1D*)gPad->FindObject("ProjY-sigma")))
381 ProjectionY(*h, *fSigmaTheta);
382
383 gPad = padsave;
384}
385
386// --------------------------------------------------------------------------
387//
388// Draw the histogram
389//
390void MHSigmaTheta::Draw(Option_t *opt)
391{
392 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
393 pad->SetBorderMode(0);
394 AppendPad("");
395
396 pad->Divide(3, 3);
397
398 // draw the 2D histogram Sigmabar versus Theta
399 TH1 *h;
400
401 pad->cd(1);
402 gPad->SetBorderMode(0);
403 h = fSigmaTheta->ProjectionX("ProjX-Theta", -1, 9999, "E");
404 h->SetDirectory(NULL);
405 h->UseCurrentStyle();
406 h->SetTitle("Distribution of \\Theta");
407 h->SetXTitle("\\Theta [\\circ]");
408 h->SetYTitle("No.of events");
409 h->SetTitleOffset(1.2,"Y");
410 h->Draw(opt);
411 h->SetBit(kCanDelete);
412
413 pad->cd(2);
414 gPad->SetBorderMode(0);
415 h = fDiffPixTheta->Project3D("zx");
416 h->SetDirectory(NULL);
417 h->UseCurrentStyle();
418 //h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. \\Theta (all pixels)");
419 h->SetTitle("(Sigma2-Sigmabar2)/Area vs. \\Theta (all pixels)");
420 h->SetXTitle("\\Theta [\\circ]");
421 h->SetYTitle("(Sigma2 - Sigmabar2) / Area");
422 h->SetTitleOffset(1.2,"Y");
423 h->Draw("box");
424 h->SetBit(kCanDelete);
425
426 pad->cd(3);
427 gPad->SetBorderMode(0);
428 h = fSigmaPixTheta->Project3D("zx");
429 h->SetDirectory(NULL);
430 h->UseCurrentStyle();
431 h->SetTitle("\\sigma_{ped} vs. \\Theta (all pixels)");
432 h->SetXTitle("\\Theta [\\circ]");
433 h->SetYTitle("Sigma");
434 h->SetTitleOffset(1.2,"Y");
435 h->Draw("box");
436 h->SetBit(kCanDelete);
437
438 //pad->cd(7);
439 //gPad->SetBorderMode(0);
440 //h = fSigmaTheta->ProjectionY("ProjY-sigma", -1, 9999, "E");
441 //h->SetDirectory(NULL);
442 //h->UseCurrentStyle();
443 //h->SetTitle("Distribution of \\bar{\\sigma}_{ped}");
444 //h->SetXTitle("\\bar{\\sigma}_{ped}");
445 //h->SetYTitle("No.of events");
446 //h->SetTitleOffset(1.2,"Y");
447 //h->Draw(opt);
448 //h->SetBit(kCanDelete);
449
450 pad->cd(5);
451 gPad->SetBorderMode(0);
452 h = fDiffPixTheta->Project3D("zy");
453 h->SetDirectory(NULL);
454 h->UseCurrentStyle();
455 //h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. pixel Id (all \\Theta)");
456 h->SetTitle("(Sigma2-Sigmabar2)/Area vs. pixel Id (all \\Theta)");
457 h->SetXTitle("Pixel Id");
458 h->SetYTitle("(Sigma2 - sigmabar2) / Area");
459 h->SetTitleOffset(1.2,"Y");
460 h->Draw("box");
461 h->SetBit(kCanDelete);
462
463 pad->cd(6);
464 gPad->SetBorderMode(0);
465 h = fSigmaPixTheta->Project3D("zy");
466 h->SetDirectory(NULL);
467 h->UseCurrentStyle();
468 h->SetTitle("\\sigma_{ped} vs. pixel Id (all \\Theta)");
469 h->SetXTitle("Pixel Id");
470 h->SetYTitle("Sigma");
471 h->SetTitleOffset(1.2,"Y");
472 h->Draw("box");
473 h->SetBit(kCanDelete);
474
475 pad->cd(4);
476 gPad->SetBorderMode(0);
477 fSigmaTheta->Draw();
478
479 pad->cd(7);
480 gPad->SetBorderMode(0);
481 fSigmaThetaOuter->Draw();
482
483 pad->cd(8);
484 gPad->SetBorderMode(0);
485 fThetaPhi->Draw();
486
487 pad->cd(9);
488 gPad->SetBorderMode(0);
489 h = fThetaPhi->ProjectionX("ProjX-Phi", -1, 9999, "E");
490 h->SetDirectory(NULL);
491 h->UseCurrentStyle();
492 h->SetTitle("Distribution of \\Phi");
493 h->SetXTitle("\\phi [\\circ]");
494 h->SetYTitle("No. of events");
495 h->SetTitleOffset(1.2,"Y");
496 h->Draw(opt);
497 h->SetBit(kCanDelete);
498}
499
500
501
502
503
504
505
506
507
508
509
510
511
512
Note: See TracBrowser for help on using the repository browser.