source: trunk/Mars/msimcamera/MSimCamera.cc@ 17663

Last change on this file since 17663 was 17663, checked in by ftemme, 10 years ago
Added the creation of a MTruePhotonsPerPixelCont in MSimCamera and the filling of the arrays of the container. Also added the container to the fitsoutput
File size: 14.5 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@phys.ethz.ch>
19!
20! Copyright: CheObs Software Development, 2000-2013
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MSimCamera
28//
29// This task initializes the analog channels with analog noise and simulated
30// the analog pulses from the photon signal.
31//
32// Input Containers:
33// MPhotonEvent
34// MPhotonStatistics
35// MRawRunHeader
36//
37// Output Containers:
38// MAnalogChannels
39//
40//////////////////////////////////////////////////////////////////////////////
41#include "MSimCamera.h"
42
43#include <TF1.h>
44#include <TRandom.h> // Needed for TRandom
45
46#include "MLog.h"
47#include "MLogManip.h"
48
49#include "MTruePhotonsPerPixelCont.h"
50
51#include "MSpline3.h"
52#include "MParSpline.h"
53
54#include "MParList.h"
55
56#include "MPhotonEvent.h"
57#include "MPhotonData.h"
58
59#include "MPedestalCam.h"
60#include "MPedestalPix.h"
61
62#include "MAnalogSignal.h"
63#include "MAnalogChannels.h"
64
65#include "MParameters.h"
66
67#include "MMcEvt.hxx" // To be replaced by a CheObs class
68#include "MRawRunHeader.h"
69
70ClassImp(MSimCamera);
71
72using namespace std;
73
74// --------------------------------------------------------------------------
75//
76// Default Constructor.
77//
78MSimCamera::MSimCamera(const char* name, const char *title)
79 : fEvt(0), fStat(0), fRunHeader(0), fElectronicNoise(0), fGain(0),
80 fCamera(0), fMcEvt(0),fCrosstalkCoeffParam(0), fSpline(0), fBaselineGain(kFALSE),
81 fDefaultOffset(-1), fDefaultNoise(-1), fDefaultGain(-1), fACFudgeFactor(0),
82 fACTimeConstant(0)
83
84{
85 fName = name ? name : "MSimCamera";
86 fTitle = title ? title : "Task to simulate the electronic noise and to convert photons into pulses";
87}
88
89// --------------------------------------------------------------------------
90//
91// Search for the necessayr parameter containers.
92// Setup spline for pulse shape.
93//
94Int_t MSimCamera::PreProcess(MParList *pList)
95{
96 fMcEvt = (MMcEvt*)pList->FindCreateObj("MMcEvt");
97 if (!fMcEvt)
98 return kFALSE;
99
100 fCamera = (MAnalogChannels*)pList->FindCreateObj("MAnalogChannels");
101 if (!fCamera)
102 return kFALSE;
103
104 fEvt = (MPhotonEvent*)pList->FindObject("MPhotonEvent");
105 if (!fEvt)
106 {
107 *fLog << err << "MPhotonEvent not found... aborting." << endl;
108 return kFALSE;
109 }
110
111 fStat = (MPhotonStatistics*)pList->FindObject("MPhotonStatistics");
112 if (!fStat)
113 {
114 *fLog << err << "MPhotonStatistics not found... aborting." << endl;
115 return kFALSE;
116 }
117
118 fRunHeader = (MRawRunHeader *)pList->FindObject("MRawRunHeader");
119 if (!fRunHeader)
120 {
121 *fLog << err << "MRawRunHeader not found... aborting." << endl;
122 return kFALSE;
123 }
124/*
125 fPulsePos = (MParameterD*)pList->FindObject("IntendedPulsePos", "MParameterD");
126 if (!fPulsePos)
127 {
128 *fLog << err << "IntendedPulsePos [MParameterD] not found... aborting." << endl;
129 return kFALSE;
130 }
131 */
132
133 // Create it here to make sure that MGeomApply will set the correct size
134 fElectronicNoise = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", "ElectronicNoise");
135 if (!fElectronicNoise)
136 return kFALSE;
137
138 fGain = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", "Gain");
139 if (!fGain)
140 return kFALSE;
141
142 fAccidentalPhotons = (MPedestalCam*)pList->FindObject("AccidentalPhotonRates","MPedestalCam");
143 if(!fAccidentalPhotons)
144 {
145 *fLog << err << "AccidentalPhotonRates [MPedestalCam] not found... aborting." << endl;
146 return kFALSE;
147 }
148
149 fCrosstalkCoeffParam = (MParameterD*)pList->FindCreateObj("MParameterD","CrosstalkCoeffParam");
150 if (!fCrosstalkCoeffParam)
151 {
152 *fLog << err << "CrosstalkCoeffParam [MParameterD] not found... aborting." << endl;
153 return kFALSE;
154 }
155
156 fTruePhotons = (MTruePhotonsPerPixelCont*)pList->FindCreateObj("MTruePhotonsPerPixelCont");
157 if (!fTruePhotons)
158 {
159 *fLog << err << "MTruePhotonsPerPixelCont not found... aborting." << endl;
160 return kFALSE;
161 }
162
163 MParSpline *pulse = (MParSpline*)pList->FindObject("PulseShape", "MParSpline");
164 if (!pulse)
165 {
166 *fLog << err << "PulseShape [MParSpline] not found... aborting." << endl;
167 return kFALSE;
168 }
169
170// if (fRunHeader->GetFreqSampling()!=1000)
171// {
172// *fLog << err << "ERROR - Sampling frequencies others than 1GHz are not yet supported." << endl;
173// *fLog << warn << "FIXME - SCALE MPulsShape WITH THE SAMPLING FREQUENCY." << endl;
174// return kFALSE;
175// }
176
177 fSpline = pulse->GetSpline();
178 if (!fSpline)
179 {
180 *fLog << err << "No spline initialized." << endl;
181 return kFALSE;
182 }
183
184 // ---------------- Information output ----------------------
185
186 if (fBaselineGain)
187 *fLog << inf << "Gain is also applied to the electronic noise." << endl;
188
189 return kTRUE;
190}
191
192// --------------------------------------------------------------------------
193//
194// FIXME: For now this is a workaround to set a baseline and the
195// electronic (guassian noise)
196//
197Bool_t MSimCamera::ReInit(MParList *plist)
198{
199 for (int i=0; i<fElectronicNoise->GetSize(); i++)
200 {
201 MPedestalPix &ped = (*fElectronicNoise)[i];
202 ped.SetPedestal(fDefaultOffset);
203 if (fDefaultNoise>0)
204 ped.SetPedestalRms(fDefaultNoise);
205
206 ped.SetPedestalABoffset(0);
207 ped.SetNumEvents(0);
208
209
210 MPedestalPix &gain = (*fGain)[i];
211 if (fDefaultGain>0)
212 gain.SetPedestal(fDefaultGain);
213
214 gain.SetPedestalRms(0);
215 gain.SetPedestalABoffset(0);
216 gain.SetNumEvents(0);
217 }
218
219 return kTRUE;
220}
221
222// --------------------------------------------------------------------------
223//
224// fStat->GetMaxIndex must return the maximum index possible
225// (equiv. number of pixels) not just the maximum index stored!
226//
227Int_t MSimCamera::Process()
228{
229 // Calculate start time, end time and corresponding number of samples
230 const Double_t freq = fRunHeader->GetFreqSampling()/1000.;
231
232 // FIXME: Should we use a higher sampling here?
233
234 const Double_t start = fStat->GetTimeFirst()*freq;
235 const Double_t end = fStat->GetTimeLast() *freq;
236
237 const UInt_t nlen = TMath::CeilNint(end-start);
238
239 // Get number of pixels/channels
240 const UInt_t npix = fStat->GetMaxIndex()+1;
241
242 if (npix>(UInt_t)fElectronicNoise->GetSize())
243 {
244 *fLog << err << "ERROR - More indices (" << npix << ") ";
245 *fLog << "assigned than existing in camera (";
246 *fLog << fElectronicNoise->GetSize() << ")!" << endl;
247 return kERROR;
248 }
249
250 const Double_t pl = fSpline->GetXmin()*freq;
251 const Double_t pr = fSpline->GetXmax()*freq;
252
253 // Init the arrays and set the range which will contain valid data
254 fCamera->Init(npix, nlen);
255 fCamera->SetValidRange(TMath::FloorNint(pr), TMath::CeilNint(nlen+pl));
256
257 // Add electronic noise to empty channels
258 for (UInt_t i=0; i<npix; i++)
259 {
260 const MPedestalPix &pix = (*fElectronicNoise)[i];
261
262 const Double_t val = pix.GetPedestal();
263 const Double_t rms = pix.GetPedestalRms();
264
265 // FTemme: Implementation of AC-coupling:
266 // to calculate the value of the accoupling per slice I use the
267 // following equation:
268 // accouplingPerSlice = accidentalPhotonRate * (1 + crossTalkProb)
269 // * areaOfOnePulse / samplingRate;
270 // Therefore I need the following variables
271 // Double_t accidentalPhotonRate; // [MHz]
272 // Float_t crossTalkProb; // [1]
273 // Double_t areaOfOnePulse; // [ADC-Counts * s]
274 // Double_t samplingRate; // [slices * MHz]
275
276 // The accidental photon rate is stored in GHz, so we have to multiply
277 // with 1E3 to get MHz:
278 const MPedestalPix &accPhoPix = (*fAccidentalPhotons)[i];
279
280 const Double_t accidentalPhotonRate = accPhoPix.GetPedestal() * 1e3; //[MHz]
281
282 Double_t currentAccidentalPhotonRate = accidentalPhotonRate;
283 if (fACTimeConstant!=0)
284 {
285 const Double_t accidentalPhotons = fACTimeConstant * accidentalPhotonRate;
286 const Double_t sigmaAccidentalPhotons = TMath::Sqrt(accidentalPhotons);
287
288 const Double_t gaus = gRandom->Gaus(accidentalPhotons,sigmaAccidentalPhotons);
289
290 currentAccidentalPhotonRate = gaus / fACTimeConstant;
291 }
292
293 // Get the CrosstalkCoefficient Parameter
294 const Double_t crossTalkProb = fCrosstalkCoeffParam->GetVal();
295
296 // To get the area of one Pulse, I only need to calculate the Integral
297 // of the Pulse Shape, which is stored in fSpline. Because the spline is
298 // normalized to a maximal amplitude of 1.0, I had to multiply it with
299 // the Default gain [ADC-Counts * s]
300 const Double_t areaOfOnePulse = fSpline->Integral() * fDefaultGain;
301
302 // The sampling rate I get from the RunHeader:
303 const Double_t samplingRate = fRunHeader->GetFreqSampling(); // [slices * MHz]
304
305 const Double_t accouplingPerSlice = currentAccidentalPhotonRate
306 * (1 + crossTalkProb + fACFudgeFactor)
307 * areaOfOnePulse / samplingRate;
308
309 // The accoupling is substracted from the timeline by decreasing the
310 // mean of the gaussian noise which is added
311
312 if (!fBaselineGain)
313 {
314 (*fCamera)[i].AddGaussianNoise(rms, val - accouplingPerSlice);
315 continue;
316 }
317 // Sorry, the name "pedestal" is misleading here
318 // FIXME: Simulate gain fluctuations
319 const Double_t gain = (*fGain)[i].GetPedestal();
320
321 // FIXME: We might add the base line here already.
322 // FIXME: How stable is the offset?
323 // FIXME: Should we write a container AppliedGain for MSImTrigger?
324
325 (*fCamera)[i].AddGaussianNoise(rms*gain, (val - accouplingPerSlice)*gain);
326 }
327
328 // FIXME: Simulate correlations with neighboring pixels
329
330 const Int_t num = fEvt->GetNumPhotons();
331
332 // A random shift, uniformely distributed within one slice, to make sure that
333 // the first photon is not always aligned identically with a sample edge.
334 // FIXME: Make it switchable
335 const Float_t rndm = gRandom->Uniform();
336
337 // FIXME: Shell we add a random shift of [0,1] samples per channel?
338 // Or maybe per channel and run?
339
340 Double_t tot = 0;
341
342 // Simulate pulses
343 for (Int_t i=0; i<num; i++)
344 {
345 const MPhotonData &ph = (*fEvt)[i];
346
347 const UInt_t idx = ph.GetTag();
348 const Double_t t = (ph.GetTime()-fStat->GetTimeFirst())*freq+rndm;// - fSpline->GetXmin();
349
350 // FIXME: Time jitter?
351 // FIXME: Add additional routing here?
352 // FIMXE: How stable is the gain?
353
354 if (ph.GetPrimary()!=MMcEvt::kNightSky && ph.GetPrimary()!=MMcEvt::kArtificial)
355 {
356 tot += ph.GetWeight();
357 (*fTruePhotons->cherenkov_photons_weight)[idx] += ph.GetWeight();
358 (*fTruePhotons->cherenkov_photons_number)[idx] += 1;
359
360 (*fTruePhotons->cherenkov_arrival_time_mean)[idx] += t;
361 (*fTruePhotons->cherenkov_arrival_time_variance)[idx] += t*t;
362
363 if (ph.GetPrimary()==MMcEvt::kMUON)
364 {
365 (*fTruePhotons->muon_cherenkov_photons_weight)[idx] += ph.GetWeight();
366 (*fTruePhotons->muon_cherenkov_photons_number)[idx] += 1;
367 }
368
369 // find min
370 if (t < (*fTruePhotons->cherenkov_arrival_time_min)[idx] )
371 {
372 (*fTruePhotons->cherenkov_arrival_time_min)[idx] = t;
373 }
374 // find max
375 if (t > (*fTruePhotons->cherenkov_arrival_time_max)[idx] )
376 {
377 (*fTruePhotons->cherenkov_arrival_time_max)[idx] = t;
378 }
379 }
380 else
381 {
382 (*fTruePhotons->noise_photons_weight)[idx] += ph.GetWeight();
383 }
384
385 // Sorry, the name "pedestal" is misleading here
386 // FIXME: Simulate gain fluctuations
387 const Double_t gain = (*fGain)[idx].GetPedestal();
388
389 // === FIXME === FIXME === FIXME === Frequency!!!!
390 (*fCamera)[idx].AddPulse(*fSpline, t, ph.GetWeight()*gain);
391 }
392
393 for (unsigned int i=0 ; i < 1440 ; i++)
394 {
395 float number = (*fTruePhotons->cherenkov_photons_number)[i];
396 (*fTruePhotons->cherenkov_arrival_time_mean)[i] /= number;
397 float mean = (*fTruePhotons->cherenkov_arrival_time_mean)[i];
398 float sum_tt = (*fTruePhotons->cherenkov_arrival_time_variance)[i];
399 (*fTruePhotons->cherenkov_arrival_time_variance)[i] = (sum_tt / number - mean*mean) /(number - 1);
400 }
401
402 fMcEvt->SetPhotElfromShower(TMath::Nint(tot));
403
404 return kTRUE;
405}
406
407// --------------------------------------------------------------------------
408//
409// BaselineGain: Off
410//
411Int_t MSimCamera::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
412{
413 Bool_t rc = kFALSE;
414 if (IsEnvDefined(env, prefix, "BaselineGain", print))
415 {
416 rc = kTRUE;
417 fBaselineGain = GetEnvValue(env, prefix, "BaselineGain", fBaselineGain);
418 }
419
420 if (IsEnvDefined(env, prefix, "DefaultOffset", print))
421 {
422 rc = kTRUE;
423 fDefaultOffset = GetEnvValue(env, prefix, "DefaultOffset", fDefaultOffset);
424 }
425 if (IsEnvDefined(env, prefix, "DefaultNoise", print))
426 {
427 rc = kTRUE;
428 fDefaultNoise = GetEnvValue(env, prefix, "DefaultNoise", fDefaultNoise);
429 }
430 if (IsEnvDefined(env, prefix, "DefaultGain", print))
431 {
432 rc = kTRUE;
433 fDefaultGain = GetEnvValue(env, prefix, "DefaultGain", fDefaultGain);
434 }
435 if (IsEnvDefined(env, prefix, "fACFudgeFactor", print))
436 {
437 rc = kTRUE;
438 fACFudgeFactor = GetEnvValue(env, prefix, "fACFudgeFactor", fACFudgeFactor);
439 }
440 if (IsEnvDefined(env, prefix, "fACTimeConstant", print))
441 {
442 rc = kTRUE;
443 fACTimeConstant = GetEnvValue(env, prefix, "fACTimeConstant", fACTimeConstant);
444 }
445
446 return rc;
447}
Note: See TracBrowser for help on using the repository browser.