source: trunk/Mars/msim/MHPhotonEvent.cc@ 12105

Last change on this file since 12105 was 10167, checked in by tbretz, 14 years ago
Treat kArtificial in addition to kNightSky.
File size: 12.6 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of CheObs, the Modular 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 appears 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, 1/2009 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: CheObs Software Development, 2000-2009
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MHPhotonEvent
28//
29// This is a histogram class to visualize the contents of a MPhotonEvent.
30//
31// Histograms for the distribution in x/y (z is just ignored), the
32// arrival direction, a profile of the arrival time vs position and
33// a spectrum is filles and displayed.
34//
35// There are several types of histograms (in the sense of binnings)
36// for several purposes:
37//
38// Type 1:
39// The maximum in x and y is determined from MCorsikaRunHeader
40// (not yet implemented. Fixed to 25000, 100 bins)
41//
42// Type 2:
43// The maximum in x and y is determined from MReflector->GetMaxR();
44// (100 bins)
45//
46// Type 3:
47// The maximum in x and y is determined from MGeomCam->GetMaxR();
48// (roughly 10 bins per pixel)
49//
50// Type 4:
51// Two times the traversal size of pixel #0 ( 4*MGeomCam[0].GetT() )
52// (101 bins, units cm)
53//
54// Type 5:
55// As type 3 but in cm instead of mm
56//
57// The binning is optimized using MH::FindGoodLimits.
58//
59// Fill expects a MPhotonEvent (the second argumnet in MFillH).
60//
61//
62// Class Version 2:
63// ----------------
64// + TH1D fHistH;
65//
66// Class Version 3:
67// ----------------
68// + replaced TH1D by TH1F
69// + replaced TH2D by TH2F
70//
71/////////////////////////////////////////////////////////////////////////////
72#include "MHPhotonEvent.h"
73
74#include <TCanvas.h>
75
76#include "MLog.h"
77#include "MLogManip.h"
78
79#include "MBinning.h"
80#include "MParList.h"
81
82#include "MGeomCam.h"
83#include "MGeomPix.h"
84
85#include "MPhotonEvent.h"
86#include "MPhotonData.h"
87
88#include "MCorsikaRunHeader.h"
89
90#include "MReflector.h"
91
92ClassImp(MHPhotonEvent);
93
94using namespace std;
95
96void MHPhotonEvent::Init(const char *name, const char *title)
97{
98 fName = name ? name : "MHPhotonEvent";
99 fTitle = title ? title : "Histogram to display the information of MPhotonEvents";
100
101 fHistXY.SetName("Position");
102 fHistXY.SetTitle("Histogram of position x/y");
103 fHistXY.SetXTitle("X [cm]");
104 fHistXY.SetYTitle("Y [cm]");
105 fHistXY.SetZTitle("Counts");
106 fHistXY.SetDirectory(NULL);
107 fHistXY.Sumw2();
108
109 fHistUV.SetName("Direction");
110 fHistUV.SetTitle("Histogram of arrival direction CosU/CosV");
111 fHistUV.SetXTitle("CosU");
112 fHistUV.SetYTitle("CosV");
113 fHistUV.SetZTitle("Counts");
114 fHistUV.SetDirectory(NULL);
115 fHistUV.Sumw2();
116
117 fHistT.SetName("Time");
118 fHistT.SetTitle("Time profile in x/y");
119 fHistT.SetXTitle("X [cm]");
120 fHistT.SetYTitle("Y [cm]");
121 fHistT.SetZTitle("T [ns]");
122 fHistT.SetDirectory(NULL);
123
124 fHistWL.SetName("Spectrum");
125 fHistWL.SetTitle("Wavelength distribution");
126 fHistWL.SetXTitle("\\lambda [nm]");
127 fHistWL.SetYTitle("Counts");
128 fHistWL.SetDirectory(NULL);
129
130 fHistH.SetName("Height");
131 fHistH.SetTitle("Production Height");
132 fHistH.SetXTitle("h [km]");
133 fHistH.SetYTitle("Counts");
134 fHistH.SetDirectory(NULL);
135}
136
137// --------------------------------------------------------------------------
138//
139// Default constructor. Creates and initializes the histograms.
140//
141MHPhotonEvent::MHPhotonEvent(Double_t max, const char *name, const char *title)
142 : fHistT("", "", 1, 0, 1, 1, 0, 1), fType(-1), fPermanentReset(kFALSE)
143{
144 // pre-initialization of the profile is necessary to get fZmin and fZmax set
145
146 Init(name, title);
147
148 const MBinning binsd(51, -max, max);
149 const MBinning binsa(51, -1, 1);
150
151 SetBinning(fHistXY, binsd, binsd);
152 SetBinning(fHistUV, binsa, binsa);
153 SetBinning(fHistT, binsd, binsd);
154}
155
156// --------------------------------------------------------------------------
157//
158// Creates and initializes the histograms.
159//
160MHPhotonEvent::MHPhotonEvent(Int_t type, const char *name, const char *title)
161 : fHistT("", "", 1, 0, 1, 1, 0, 1), fType(type), fPermanentReset(kFALSE)
162{
163 // pre-initialization of the profile is necessary to get fZmin and fZmax set
164
165 Init(name, title);
166
167 const MBinning bins(51, -1, 1);
168
169 SetBinning(fHistUV, bins, bins);
170}
171
172// --------------------------------------------------------------------------
173//
174// Find good limits for a binning num x [-max;max]
175// and apply it to fHistXY and fHistT.
176//
177void MHPhotonEvent::SetBinningXY(Int_t num, Double_t max)
178{
179 Double_t min = -max;
180
181 MH::FindGoodLimits(num, num, min, max, kFALSE);
182 MH::FindGoodLimits(num, num, min, max, kFALSE);
183
184 const MBinning binsd(num+1, min, max);
185
186 SetBinning(fHistXY, binsd, binsd);
187 SetBinning(fHistT, binsd, binsd);
188}
189
190// --------------------------------------------------------------------------
191//
192// Search for MRflEvtData, MRflEvtHeader and MGeomCam
193//
194Bool_t MHPhotonEvent::SetupFill(const MParList *pList)
195{
196 Double_t xmax = -1;
197 Int_t num = 100;
198
199 const Int_t f = fPermanentReset ? 2 : 1;
200 MBinning(100/f, 0, 25).Apply(fHistH);
201 MBinning(70, 275, 625).Apply(fHistWL);
202
203 switch (fType)
204 {
205 case -1:
206 return kTRUE;
207 // case0: Take a value defined by the user
208 case 1:
209 xmax = 25000;
210 break;
211 case 2:
212 {
213 MReflector *r = (MReflector*)pList->FindObject("Reflector", "MReflector");
214 if (!r)
215 {
216 *fLog << err << "Reflector [MReflector] not found... aborting." << endl;
217 return kFALSE;
218 }
219 xmax = r->GetMaxR();
220 break;
221 }
222 case 3: // The maximum radius (mm)
223 fHistXY.SetXTitle("X [mm]");
224 fHistXY.SetYTitle("Y [mm]");
225 fHistT.SetXTitle("X [mm]");
226 fHistT.SetYTitle("Y [mm]");
227 // *Fallthrough*
228
229 case 4: // Two times the pixel-0 traversal size
230 case 5: // The maximum radius (cm) -- ~ 10 bins / pix
231 case 6: // The maximum radius (cm) -- ~ 5 bins / pix
232 case 7: // The maximum radius (cm) -- ~ 3 bins / pix
233 case 8: // The maximum radius (cm) -- ~ 2 bins / pix
234 {
235 MGeomCam *c = (MGeomCam*)pList->FindObject("MGeomCam");
236 if (!c)
237 {
238 *fLog << err << "MGeomCam not found... aborting." << endl;
239 return kFALSE;
240 }
241 // Type 3: Define ~10 bins per pixel
242 xmax = fType!=4 ? c->GetMaxRadius() : 4*(*c)[0].GetT()/10;
243
244 switch (fType)
245 {
246 case 3:
247 case 5:
248 // ~10 bins per pixel
249 num = TMath::Nint(10*c->GetMaxRadius()/(*c)[0].GetT());
250 break;
251 case 6:
252 // ~5 bins per pixel
253 num = TMath::Nint(5*c->GetMaxRadius()/(*c)[0].GetT());
254 break;
255 case 7:
256 // ~3 bins per pixel
257 num = TMath::Nint(3*c->GetMaxRadius()/(*c)[0].GetT());
258 break;
259 case 8:
260 // ~2 bins per pixel
261 num = TMath::Nint(2*c->GetMaxRadius()/(*c)[0].GetT());
262 break;
263 case 4:
264 num = 101;
265 break;
266 }
267
268 if (fType==5 || fType==6 || fType==7)
269 xmax /= 10;
270
271 break;
272 }
273 default:
274 *fLog << err << "No valid binning (Type=" << fType << ") given... aborting." << endl;
275 return kFALSE;
276 }
277
278 SetBinningXY(num, xmax);
279
280 return kTRUE;
281}
282
283// --------------------------------------------------------------------------
284//
285// ReInit overwrites the binning of the Wavelength histogram
286// by the wavlenegth band from MCorsikaRunHeader.
287// The bin-width is 5nm (fPermanentReset==kTRUE) or 10nm (kFALSE).
288//
289Bool_t MHPhotonEvent::ReInit(MParList *pList)
290{
291 MCorsikaRunHeader *h = (MCorsikaRunHeader*)pList->FindObject("MCorsikaRunHeader");
292 if (!h)
293 {
294 *fLog << err << "MCorsikaRunHeader not found... aborting." << endl;
295 return kFALSE;
296 }
297
298 /*
299 // What is the size of the light pool onm the ground?
300 if (fType==1)
301 SetBinningXY(100, h->GetImpactMax());
302 */
303
304 const Int_t f = fPermanentReset ? 10 : 2;
305
306 Double_t xmin = h->GetWavelengthMin()-20;
307 Double_t xmax = h->GetWavelengthMax()+20;
308 Int_t num = TMath::CeilNint((xmax-xmin)/f);
309
310 MH::FindGoodLimits(num, num, xmin, xmax, kTRUE);
311
312 MBinning(TMath::Abs(num), xmin-.5, xmax-.5).Apply(fHistWL);
313
314 return kTRUE;
315}
316
317// --------------------------------------------------------------------------
318//
319// Fill contents of MPhotonEvent into histograms
320//
321Int_t MHPhotonEvent::Fill(const MParContainer *par, const Stat_t weight)
322{
323 const MPhotonEvent *evt = dynamic_cast<const MPhotonEvent*>(par);
324 if (!evt)
325 {
326 *fLog << err << dbginf << "No MPhotonEvent found..." << endl;
327 return kERROR;
328 }
329
330 // Check if we want to use this class as a single event display
331 if (fPermanentReset && evt->GetNumPhotons()>0)
332 Clear();
333
334 // Get number of photons
335 const Int_t num = evt->GetNumPhotons();
336
337 // Set minimum to maximum possible value
338 Double_t min = FLT_MAX;
339
340 // Loop over all photons and determine the time of the first photon
341 // FIXME: Should we get it from some statistics container?
342 for (Int_t idx=0; idx<num; idx++)
343 min = TMath::Min(min, (*evt)[idx].GetTime());
344
345 // Now fill all histograms
346 for (Int_t idx=0; idx<num; idx++)
347 {
348 const MPhotonData &ph = (*evt)[idx];
349
350 if (ph.GetPrimary()==MMcEvtBasic::kNightSky || ph.GetPrimary()==MMcEvtBasic::kArtificial)
351 continue;
352
353 const Double_t x = ph.GetPosX();
354 const Double_t y = ph.GetPosY();
355 const Double_t u = ph.GetCosU();
356 const Double_t v = ph.GetCosV();
357 const Double_t t = ph.GetTime()-min;
358 const Double_t w = ph.GetWavelength();
359 const Double_t h = ph.GetProductionHeight()/100000;
360
361 //TVector3 dir = dat->GetDir3();
362 //dir *= -1./dir.Z();
363
364 fHistXY.Fill(x, y);
365 fHistUV.Fill(u, v);
366 fHistT.Fill(x, y, t);
367 fHistWL.Fill(w);
368 fHistH.Fill(h);
369 }
370
371 return kTRUE;
372}
373
374// --------------------------------------------------------------------------
375//
376/*
377Bool_t MHPhotonEvent::Finalize()
378{
379 const TAxis &axey = *fHistRad.GetYaxis();
380 for (int y=1; y<=fHistRad.GetNbinsY(); y++)
381 {
382 const Double_t lo = axey.GetBinLowEdge(y);
383 const Double_t hi = axey.GetBinLowEdge(y+1);
384
385 const Double_t A = (hi*hi-lo*lo)*TMath::Pi()*TMath::Pi();
386
387 for (int x=1; x<=fHistRad.GetNbinsX(); x++)
388 fHistRad.SetBinContent(x, y, fHistRad.GetBinContent(x, y)/A);
389 }
390 return kTRUE;
391}
392*/
393
394// --------------------------------------------------------------------------
395//
396// Make sure that the TProfile2D doesn't fix it's displayed minimum at 0.
397// Set the pretty-palette.
398//
399void MHPhotonEvent::Paint(Option_t *opt)
400{
401 SetPalette("pretty");
402
403 fHistT.SetMinimum();
404 fHistT.SetMinimum(fHistT.GetMinimum(0));
405}
406
407// --------------------------------------------------------------------------
408//
409void MHPhotonEvent::Draw(Option_t *o)
410{
411 TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
412 pad->SetBorderMode(0);
413
414 AppendPad();
415
416 pad->Divide(3,2);
417
418 pad->cd(1);
419 gPad->SetBorderMode(0);
420 gPad->SetGrid();
421 fHistXY.Draw("colz");
422
423 pad->cd(2);
424 gPad->SetBorderMode(0);
425 gPad->SetGrid();
426 fHistT.Draw("colz");
427
428 pad->cd(4);
429 gPad->SetBorderMode(0);
430 gPad->SetGrid();
431 fHistUV.Draw("colz");
432
433 pad->cd(5);
434 gPad->SetBorderMode(0);
435 gPad->SetGrid();
436 fHistWL.Draw();
437
438 pad->cd(6);
439 gPad->SetBorderMode(0);
440 gPad->SetGrid();
441 fHistH.Draw();
442}
443
444// --------------------------------------------------------------------------
445//
446// PermanentReset: Off
447//
448Int_t MHPhotonEvent::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
449{
450 Bool_t rc = kFALSE;
451 if (IsEnvDefined(env, prefix, "PermanentReset", print))
452 {
453 rc = kTRUE;
454 fPermanentReset = GetEnvValue(env, prefix, "PermanentReset", fPermanentReset);
455 }
456
457 return rc;
458}
Note: See TracBrowser for help on using the repository browser.