source: trunk/Mars/msimcamera/MSimRandomPhotons.cc@ 10092

Last change on this file since 10092 was 10092, checked in by tbretz, 12 years ago
Fixed a mistake in the log-output of MSimRandomPhotons.
File size: 18.9 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// MSimRandomPhotons
28//
29// Simulate poissonian photons. Since the distribution of the arrival time
30// differences of these photons is an exonential we can simulate them
31// using exponentially distributed time differences between two consecutive
32// photons.
33//
34// FIXME: We should add the wavelength distribution.
35//
36// The artificial night sky background rate is calculated as follows:
37//
38// * The photon detection efficiency vs. wavelength of the detector is obtained
39// from "PhotonDetectionEfficiency" of type "MParSpline"
40//
41// * The angular acceptance of the light collectors is obtained
42// from "ConesAngularAcceptance" of type "MParSpline"
43//
44// * The spectral acceptance of the light collectors is obtained
45// from "ConesTransmission" of type "MParSpline"
46//
47// * The reflectivity of the mirrors vs wavelength is obtained
48// from "MirrorReflectivity" of type "MParSpline"
49//
50// The rate is then calculated as
51//
52// R = R0 * Ai * f
53//
54// R0 is the night sky background rate as given in Eckart's paper (divided
55// by the wavelength window). Ai the area of the cones acceptance window,
56// f is given as:
57//
58// f = nm * Min(Ar, sr*d^2)
59//
60// with
61//
62// nm being the integral of the product of the mirror reflectivity, the cone
63// transmission and the photon detection efficiency.
64//
65// d the distance of the focal plane to the mirror
66//
67// Ar is the total reflective area of the reflector
68//
69// sr is the effective solid angle corresponding to the integral of the
70// cones angular acceptance
71//
72// Alternatively, the night-sky background rate can be calculated from
73// a spectrum as given in Fig. 1 (but versus Nanometers) in
74//
75// Chris R. Benn & Sara L. Ellison La Palma Night-Sky Brightness
76//
77// After proper conversion of the units, the rate of the pixel 0
78// is then calculated by
79//
80// rate = f * nsb
81//
82// With nsb
83//
84// nsb = Integral(nsb spectrum * combines efficiencies)
85//
86// and f can be either
87//
88// Eff. angular acceptance Cones (e.g. 20deg) * Cone-Area (mm^2)
89// f = sr * A0
90//
91// or
92//
93// Mirror-Area * Field of view of cones (deg^2)
94// f = Ar * A0;
95//
96//
97// Input Containers:
98// fNameGeomCam [MGeomCam]
99// MPhotonEvent
100// MPhotonStatistics
101// MCorsikaEvtHeader
102// [MCorsikaRunHeader]
103//
104// Output Containers:
105// MPhotonEvent
106// AccidentalPhotonRate [MPedestalCam]
107//
108//////////////////////////////////////////////////////////////////////////////
109#include "MSimRandomPhotons.h"
110
111#include <TRandom.h>
112
113#include "MMath.h" // RndmExp
114
115#include "MLog.h"
116#include "MLogManip.h"
117
118#include "MParList.h"
119
120#include "MGeomCam.h"
121#include "MGeom.h"
122
123#include "MPhotonEvent.h"
124#include "MPhotonData.h"
125
126#include "MPedestalCam.h"
127#include "MPedestalPix.h"
128
129#include "MCorsikaRunHeader.h"
130
131#include "MSpline3.h"
132#include "MParSpline.h"
133#include "MReflector.h"
134
135ClassImp(MSimRandomPhotons);
136
137using namespace std;
138
139// --------------------------------------------------------------------------
140//
141// Default Constructor.
142//
143MSimRandomPhotons::MSimRandomPhotons(const char* name, const char *title)
144 : fGeom(0), fEvt(0), fStat(0), /*fEvtHeader(0),*/ fRunHeader(0),
145 fRates(0), fSimulateWavelength(kFALSE), fNameGeomCam("MGeomCam"),
146 fFileNameNSB("resmc/night-sky-la-palma.txt")
147{
148 fName = name ? name : "MSimRandomPhotons";
149 fTitle = title ? title : "Simulate possonian photons (like NSB or dark current)";
150}
151
152// --------------------------------------------------------------------------
153//
154// Check for the necessary containers
155//
156Int_t MSimRandomPhotons::PreProcess(MParList *pList)
157{
158 fGeom = (MGeomCam*)pList->FindObject(fNameGeomCam, "MGeomCam");
159 if (!fGeom)
160 {
161 *fLog << inf << fNameGeomCam << " [MGeomCam] not found..." << endl;
162
163 fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
164 if (!fGeom)
165 {
166 *fLog << err << "MGeomCam not found... aborting." << endl;
167 return kFALSE;
168 }
169 }
170
171 fEvt = (MPhotonEvent*)pList->FindObject("MPhotonEvent");
172 if (!fEvt)
173 {
174 *fLog << err << "MPhotonEvent not found... aborting." << endl;
175 return kFALSE;
176 }
177
178 fStat = (MPhotonStatistics*)pList->FindObject("MPhotonStatistics");
179 if (!fStat)
180 {
181 *fLog << err << "MPhotonStatistics not found... aborting." << endl;
182 return kFALSE;
183 }
184
185 fRates = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", "AccidentalPhotonRates");
186 if (!fRates)
187 return kFALSE;
188
189 /*
190 fEvtHeader = (MCorsikaEvtHeader*)pList->FindObject("MCorsikaEvtHeader");
191 if (!fEvtHeader)
192 {
193 *fLog << err << "MCorsikaEvtHeader not found... aborting." << endl;
194 return kFALSE;
195 }*/
196
197 fRunHeader = (MCorsikaRunHeader*)pList->FindObject("MCorsikaRunHeader");
198 if (fSimulateWavelength && !fRunHeader)
199 {
200 *fLog << err << "MCorsikaRunHeader not found... aborting." << endl;
201 return kFALSE;
202 }
203
204 MReflector *r = (MReflector*)pList->FindObject("Reflector", "MReflector");
205 if (!r)
206 {
207 *fLog << err << "Reflector [MReflector] not found... aborting." << endl;
208 return kFALSE;
209 }
210
211 const MParSpline *s1 = (MParSpline*)pList->FindObject("PhotonDetectionEfficiency", "MParSpline");
212 const MParSpline *s2 = (MParSpline*)pList->FindObject("ConesTransmission", "MParSpline");
213 const MParSpline *s3 = (MParSpline*)pList->FindObject("MirrorReflectivity", "MParSpline");
214 const MParSpline *s4 = (MParSpline*)pList->FindObject("ConesAngularAcceptance", "MParSpline");
215
216 // Ensure that all efficiencies are at least defined in the raneg of the
217 // photon detection efficiency. We assume that this is the limiting factor
218 // and has to be zero at both ends.
219 if (s2->GetXmin()>s1->GetXmin())
220 {
221 *fLog << err << "ERROR - ConeTransmission range must be defined down to " << s1->GetXmin() << "nm (PhotonDetectionEffciency)." << endl;
222 return kFALSE;
223 }
224 if (s2->GetXmax()<s1->GetXmax())
225 {
226 *fLog << err << "ERROR - ConeTransmission range must be defined up to " << s1->GetXmax() << "nm (PhotonDetectionEffciency)." << endl;
227 return kFALSE;
228 }
229 if (s3->GetXmin()>s1->GetXmin())
230 {
231 *fLog << err << "ERROR - MirrorReflectivity range must be defined down to " << s1->GetXmin() << "nm (PhotonDetectionEffciency)." << endl;
232 return kFALSE;
233 }
234 if (s3->GetXmax()<s1->GetXmax())
235 {
236 *fLog << err << "ERROR - MirrorReflectivity range must be defined up to " << s1->GetXmax() << "nm (PhotonDetectionEffciency)." << endl;
237 return kFALSE;
238 }
239
240 // If the simulated wavelength range exists and is smaller, reduce the
241 // range to it. Later it is checked that at both edges the transmission
242 // is 0. This must be true in both cases: The simulated wavelength range
243 // exceed the PDE or the PDE range exceeds the simulated waveband.
244 const Float_t wmin = fRunHeader && fRunHeader->GetWavelengthMin()>s1->GetXmin() ? fRunHeader->GetWavelengthMin() : s1->GetXmin();
245 const Float_t wmax = fRunHeader && fRunHeader->GetWavelengthMax()<s1->GetXmax() ? fRunHeader->GetWavelengthMax() : s1->GetXmax();
246
247 const Int_t min = TMath::FloorNint(wmin);
248 const Int_t max = TMath::CeilNint(wmax);
249
250 // Multiply all relevant efficiencies to get the total transmission
251 MParSpline eff;
252 eff.SetFunction("1", max-min, min, max);
253
254 eff.Multiply(*s1->GetSpline());
255 eff.Multiply(*s2->GetSpline());
256 eff.Multiply(*s3->GetSpline());
257
258 // Effectively transmitted wavelength band in the simulated range
259 const Double_t nm = eff.GetSpline()->Integral();
260
261 // Angular acceptance of the cones
262 const Double_t sr = s4 && s4->GetSpline() ? s4->GetSpline()->IntegralSolidAngle() : 1;
263
264 {
265 const Double_t d2 = fGeom->GetCameraDist()*fGeom->GetCameraDist();
266 const Double_t conv = fGeom->GetConvMm2Deg()*TMath::DegToRad();
267 const Double_t f1 = TMath::Min(r->GetA()/1e4, sr*d2) * conv*conv;
268
269 // Rate in GHz / mm^2
270 fScale = fFreqNSB * nm * f1; // [GHz/mm^2] efficiency * m^2 *rad^2 *mm^2
271
272 const Double_t freq0 = fScale*(*fGeom)[0].GetA()*1000;
273
274 *fLog << inf << "Resulting Freq. in " << fNameGeomCam << "[0]: " << Form("%.2f", freq0) << "MHz" << endl;
275
276 // FIXME: Scale with the number of pixels
277 if (freq0>1000)
278 {
279 *fLog << err << "ERROR - Frequency exceeds 1GHz, this might leed to too much memory consumption." << endl;
280 return kFALSE;
281 }
282 }
283
284 if (fFileNameNSB.IsNull())
285 return kTRUE;
286
287 // const MMcRunHeader *mcrunheader = (MMcRunHeader*)pList->FindObject("MMcRunHeader");
288 // Set NumPheFromDNSB
289
290 // # Number of photons from the diffuse NSB (nphe / ns 0.1*0.1 deg^2 239 m^2) and
291 // nsb_mean 0.20
292 // Magic pixel: 0.00885361 deg
293 // dnsbpix = 0.2*50/15
294 // ampl = MMcFadcHeader->GetAmplitud()
295 // sqrt(pedrms*pedrms + dnsbpix*ampl*ampl/ratio)
296
297 // Conversion of the y-axis
298 // ------------------------
299 // Double_t ff = 1; // myJy / arcsec^2 per nm
300 // ff *= 1e-6; // Jy / arcsec^2 per nm
301 // ff *= 3600*3600; // Jy / deg^2
302 // ff *= 1./TMath::DegToRad()/TMath::DegToRad(); // Jy/sr = 1e-26J/s/m^2/Hz/sr
303 // ff *= 1e-26; // J/s/m^2/Hz/sr per nm
304
305 const Double_t arcsec2rad = TMath::DegToRad()/3600.;
306 const Double_t f = 1e-32 / (arcsec2rad*arcsec2rad);
307
308 // Read night sky background flux from file
309 MParSpline flux;
310 if (!flux.ReadFile(fFileNameNSB))
311 return kFALSE;
312
313 if (flux.GetXmin()>wmin)
314 {
315 *fLog << err << "ERROR - NSB flux from " << fFileNameNSB << " must be defined down to " << wmin << "nm." << endl;
316 return kFALSE;
317 }
318 if (flux.GetXmax()<wmax)
319 {
320 *fLog << err << "ERROR - NSB flux from " << fFileNameNSB << " must be defined up to " << wmax << "nm." << endl;
321 return kFALSE;
322 }
323
324 MParSpline nsb;
325
326 // Normalization to our units,
327 // conversion from energy flux to photon flux
328 nsb.SetFunction(Form("%.12e/(x*TMath::H())", f), max-min, min, max);
329
330 // multiply night sky background flux with normalization
331 nsb.Multiply(*flux.GetSpline());
332
333 // Multiply with the total transmission
334 nsb.Multiply(*eff.GetSpline());
335
336 // Check if the photon flux is zero at both ends of the NSB
337 if (eff.GetSpline()->Eval(min)>1e-5)
338 {
339 *fLog << warn << "WARNING - Total transmission efficiency at ";
340 *fLog << min << "nm is not zero, but " << eff.GetSpline()->Eval(min) << "... abort." << endl;
341 }
342 if (eff.GetSpline()->Eval(max)>1e-5)
343 {
344 *fLog << warn << "WARNING - Total transmission efficiency at ";
345 *fLog << max << "nm is not zero, but " << eff.GetSpline()->Eval(max) << "... abort." << endl;
346 }
347
348 // Check if the photon flux is zero at both ends of the simulated region
349 if (eff.GetSpline()->Eval(wmin)>1e-5)
350 {
351 *fLog << err << "ERROR - Total transmission efficiency at ";
352 *fLog << wmin << "nm is not zero... abort." << endl;
353 *fLog << " PhotonDetectionEfficency: " << s1->GetSpline()->Eval(wmin) << endl;
354 *fLog << " ConeTransmission: " << s2->GetSpline()->Eval(wmin) << endl;
355 *fLog << " MirrorReflectivity: " << s3->GetSpline()->Eval(wmin) << endl;
356 *fLog << " TotalEfficiency: " << eff.GetSpline()->Eval(wmin) << endl;
357 return kFALSE;
358 }
359 if (eff.GetSpline()->Eval(wmax)>1e-5)
360 {
361 *fLog << err << "ERROR - Total transmission efficiency at ";
362 *fLog << wmax << "nm is not zero... abort." << endl;
363 *fLog << " PhotonDetectionEfficency: " << s1->GetSpline()->Eval(wmax) << endl;
364 *fLog << " ConeTransmission: " << s2->GetSpline()->Eval(wmax) << endl;
365 *fLog << " MirrorReflectivity: " << s3->GetSpline()->Eval(wmax) << endl;
366 *fLog << " TotalEfficiency: " << eff.GetSpline()->Eval(wmax) << endl;
367 return kFALSE;
368 }
369
370 // Conversion from m to radians
371 const Double_t conv = fGeom->GetConvMm2Deg()*TMath::DegToRad()*1e3;
372
373 // Angular acceptance of the cones
374 //const Double_t sr = s5.GetSpline()->IntegralSolidAngle(); // sr
375 // Absolute reflector area
376 const Double_t Ar = r->GetA()/1e4; // m^2
377 // Size of the cone's entrance window
378 const Double_t A0 = (*fGeom)[0].GetA()*1e-6; // m^2
379
380 // Rate * m^2 * Solid Angle
381 // -------------------------
382
383 // Angular acceptance Cones (e.g. 20deg) * Cone-Area
384 const Double_t f1 = A0 * sr; // m^2 sr
385
386 // Mirror-Area * Field of view of cones (e.g. 0.1deg)
387 const Double_t f2 = Ar * A0*conv*conv; // m^2 sr
388
389 // FIXME: Calculate the reflectivity of the bottom by replacing
390 // MirrorReflectivity by bottom reflectivity and reflect
391 // and use it to reflect the difference between f1 and f2
392 // if any.
393
394 // Total NSB rate in MHz per m^2 and sr
395 const Double_t rate = nsb.GetSpline()->Integral() * 1e-6;
396
397 *fLog << inf;
398
399 // Resulting rates as if Razmick's constant had been used
400 // *fLog << 1.75e6/(600-300) * f1 * eff.GetSpline()->Integral() << " MHz" << endl;
401 // *fLog << 1.75e6/(600-300) * f2 * eff.GetSpline()->Integral() << " MHz" << endl;
402
403 *fLog << "Conversion factor Fnu: " << f << endl;
404 *fLog << "Total reflective area: " << Form("%.2f", Ar) << " m" << UTF8::kSquare << endl;
405 *fLog << "Acceptance area of cone 0: " << Form("%.2f", A0*1e6) << " mm" << UTF8::kSquare << " = ";
406 *fLog << A0*conv*conv << " sr" << endl;
407 *fLog << "Cones angular acceptance: " << sr << " sr" << endl;
408 *fLog << "ConeArea*ConeSolidAngle (f1): " << f1 << " m^2 sr" << endl;
409 *fLog << "MirrorArea*ConeSkyAngle (f2): " << f2 << " m^2 sr" << endl;
410 *fLog << "Effective. transmission: " << Form("%.1f", nm) << " nm" << endl;
411 *fLog << "NSB freq. in " << fNameGeomCam << "[0] (f1): " << Form("%.2f", rate * f1) << " MHz" << endl;
412 *fLog << "NSB freq. in " << fNameGeomCam << "[0] (f2): " << Form("%.2f", rate * f2) << " MHz" << endl;
413 *fLog << "Using f2." << endl;
414
415 // Scale the rate per mm^2 and to GHz
416 fScale = rate * f2 / (*fGeom)[0].GetA() / 1000;
417
418 // FIXME: Scale with the number of pixels
419 if (rate*f2>1000)
420 {
421 *fLog << err << "ERROR - Frequency exceeds 1GHz, this might leed to too much memory consumption." << endl;
422 return kFALSE;
423 }
424
425 return kTRUE;
426}
427
428Bool_t MSimRandomPhotons::ReInit(MParList *pList)
429{
430 // Overwrite the default set by MGeomApply
431 fRates->Init(*fGeom);
432 return kTRUE;
433}
434
435// --------------------------------------------------------------------------
436//
437// Check for the necessary containers
438//
439Int_t MSimRandomPhotons::Process()
440{
441 // Get array from event container
442 // const Int_t num = fEvt->GetNumPhotons();
443 //
444 // Do not produce pure pedestal events!
445 // if (num==0)
446 // return kTRUE;
447
448 // Get array from event container
449 // FIXME: Use statistics container instead
450 const UInt_t npix = fGeom->GetNumPixels();
451
452 // This is the possible window in which the triggered digitization
453 // may take place.
454 const Double_t start = fStat->GetTimeFirst();
455 const Double_t end = fStat->GetTimeLast();
456
457 // Loop over all pixels
458 for (UInt_t idx=0; idx<npix; idx++)
459 {
460 // Scale the rate with the pixel size.
461 const Double_t rate = fFreqFixed + fScale*(*fGeom)[idx].GetA();
462
463 (*fRates)[idx].SetPedestal(rate);
464
465 // Calculate the average distance between two consequtive photons
466 const Double_t avglen = 1./rate;
467
468 // Start producing photons at time "start"
469 Double_t t = start;
470 while (1)
471 {
472 // Get a random time for the photon.
473 // The differences are exponentially distributed.
474 t += MMath::RndmExp(avglen);
475
476 // Check if we reached the end of the useful time window
477 if (t>end)
478 break;
479
480 // Add a new photon
481 // FIXME: SLOW!
482 MPhotonData &ph = fEvt->Add();
483
484 // Set source to NightSky, time to t and tag to pixel index
485 ph.SetPrimary(MMcEvtBasic::kNightSky);
486 ph.SetWeight();
487 ph.SetTime(t);
488 ph.SetTag(idx);
489
490 // fProductionHeight, fPosX, fPosY, fCosU, fCosV (irrelevant) FIXME: Reset?
491
492 if (fSimulateWavelength)
493 {
494 const Float_t wmin = fRunHeader->GetWavelengthMin();
495 const Float_t wmax = fRunHeader->GetWavelengthMax();
496
497 ph.SetWavelength(TMath::Nint(gRandom->Uniform(wmin, wmax)));
498 }
499 }
500 }
501
502 // Re-sort the photons by time!
503 fEvt->Sort(kTRUE);
504
505 // Update maximum index
506 fStat->SetMaxIndex(npix-1);
507
508 // Shrink
509 return kTRUE;
510}
511
512// --------------------------------------------------------------------------
513//
514// Read the parameters from the resource file.
515//
516// FrequencyFixed: 0.040
517// FrequencyNSB: 5.8
518//
519// The fixed frequency is given in units fitting the units of the time.
520// Usually the time is given in nanoseconds thus, e.g., 0.040 means 40MHz.
521//
522// The FrequencyNSB is scaled by the area of the pixel in cm^2. Therefore
523// 0.040 would mean 40MHz/cm^2
524//
525Int_t MSimRandomPhotons::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
526{
527 Bool_t rc = kFALSE;
528 if (IsEnvDefined(env, prefix, "FrequencyFixed", print))
529 {
530 rc = kTRUE;
531 fFreqFixed = GetEnvValue(env, prefix, "FrequencyFixed", fFreqFixed);
532 }
533
534 if (IsEnvDefined(env, prefix, "FrequencyNSB", print))
535 {
536 rc = kTRUE;
537 fFreqNSB = GetEnvValue(env, prefix, "FrequencyNSB", fFreqNSB);
538 }
539
540 if (IsEnvDefined(env, prefix, "FileNameNSB", print))
541 {
542 rc = kTRUE;
543 fFileNameNSB = GetEnvValue(env, prefix, "FileNameNSB", fFileNameNSB);
544 }
545
546 if (IsEnvDefined(env, prefix, "SimulateCherenkovSpectrum", print))
547 {
548 rc = kTRUE;
549 fSimulateWavelength = GetEnvValue(env, prefix, "SimulateCherenkovSpectrum", fSimulateWavelength);
550 }
551
552 return rc;
553}
Note: See TracBrowser for help on using the repository browser.