source: branches/MarsMoreSimulationTruth/msimcamera/MSimCamera.cc@ 20067

Last change on this file since 20067 was 18333, checked in by Jens Buss, 9 years ago
Reintegrated private branch for the simulation of a gapd time jitter to the trunk
File size: 19.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@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 // Dominik Neise and Sebastian Mueller on fix time offsets:
126 // We obtain the fix temporal offsets for the FACT camera pixels out of
127 // a text file. The textfile must be mentioned in the ceres.rc file.
128 // There are no default offsets on purporse. The filename must be specified
129 // in ceres.rc and the file must be parsed without errors and it must
130 // provide exactly 1440 floating point numbers.
131 fFixTimeOffsetsBetweenPixelsInNs =
132 (MMatrix*)pList->FindObject("MFixTimeOffset");
133 if (!fFixTimeOffsetsBetweenPixelsInNs)
134 {
135 // the key value pair providing the text file is not present in the
136 // environment env.
137 *fLog << err << "In Source: "<< __FILE__ <<" in line: "<< __LINE__;
138 *fLog << " in function: "<< __func__ <<"\n";
139 *fLog << "MFixTimeOffset not found... aborting." << endl;
140 return kFALSE;
141
142 }
143 else if ( fFixTimeOffsetsBetweenPixelsInNs->fM.size() != 1440 )
144 {
145 // The number of time offsets must match the number of pixels in the
146 // FACT camera.
147 *fLog << err << "In Source: "<< __FILE__ <<" in line: "<< __LINE__;
148 *fLog << " in function: "<< __func__ <<"\n";
149 *fLog << "MFixTimeOffset has the wrong dimension! ";
150 *fLog << "There should be "<< 1440 <<" time offsets ";
151 *fLog << "(one for each pixel in FACT) but there are: ";
152 *fLog << fFixTimeOffsetsBetweenPixelsInNs->fM.size() << "! ";
153 *fLog << "... aborting." << endl;
154 return kFALSE;
155 }
156 // Check all entries for inf and nan. Those are not accepted here.
157 for( int row_index=0; row_index<fFixTimeOffsetsBetweenPixelsInNs->fM.size(); row_index++){
158 std::vector<double> row = fFixTimeOffsetsBetweenPixelsInNs->fM.at(row_index);
159 for( int col_index=0; col_index<row.size(); col_index++){
160 double specific_delay = row.at(col_index);
161 if( std::isnan(specific_delay) || std::isinf(specific_delay) ){
162 *fLog << err << "In Source: "<< __FILE__ <<" in line: ";
163 *fLog << __LINE__;
164 *fLog << " in function: "<< __func__ <<"\n";
165 *fLog << "There is a non normal specific_delay in the fix temporal ";
166 *fLog << "pixel offsets. This is that at least one specific_delay is ";
167 *fLog << "NaN or Inf. This here is >"<< specific_delay;
168 *fLog << "<... aborting." << endl;
169 return kFALSE;
170 }
171 }
172
173 }
174 // -------------------------------------------------------------------
175/*
176 fPulsePos = (MParameterD*)pList->FindObject("IntendedPulsePos", "MParameterD");
177 if (!fPulsePos)
178 {
179 *fLog << err << "IntendedPulsePos [MParameterD] not found... aborting." << endl;
180 return kFALSE;
181 }
182 */
183 fResidualTimeSpread = (MParameterD*)pList->FindObject("ResidualTimeSpread");
184 if (!fResidualTimeSpread)
185 {
186 *fLog << err << "ResidualTimeSpread [MParameterD] not found... aborting." << endl;
187 return kFALSE;
188 }
189
190 // Get GapdTimeJitter from parameter list
191 fGapdTimeJitter = (MParameterD*)pList->FindObject("GapdTimeJitter");
192 if (!fGapdTimeJitter)
193 {
194 *fLog << err << "GapdTimeJitter [MParameterD] not found... aborting." << endl;
195 return kFALSE;
196 }
197
198 // Create it here to make sure that MGeomApply will set the correct size
199 fElectronicNoise = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", "ElectronicNoise");
200 if (!fElectronicNoise)
201 return kFALSE;
202
203 fGain = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", "Gain");
204 if (!fGain)
205 return kFALSE;
206
207 fAccidentalPhotons = (MPedestalCam*)pList->FindObject("AccidentalPhotonRates","MPedestalCam");
208 if(!fAccidentalPhotons)
209 {
210 *fLog << err << "AccidentalPhotonRates [MPedestalCam] not found... aborting." << endl;
211 return kFALSE;
212 }
213
214 fCrosstalkCoeffParam = (MParameterD*)pList->FindCreateObj("MParameterD","CrosstalkCoeffParam");
215 if (!fCrosstalkCoeffParam)
216 {
217 *fLog << err << "CrosstalkCoeffParam [MParameterD] not found... aborting." << endl;
218 return kFALSE;
219 }
220
221 fTruePhotons = (MTruePhotonsPerPixelCont*)pList->FindCreateObj("MTruePhotonsPerPixelCont");
222 if (!fTruePhotons)
223 {
224 *fLog << err << "MTruePhotonsPerPixelCont not found... aborting." << endl;
225 return kFALSE;
226 }
227
228 MParSpline *pulse = (MParSpline*)pList->FindObject("PulseShape", "MParSpline");
229 if (!pulse)
230 {
231 *fLog << err << "PulseShape [MParSpline] not found... aborting." << endl;
232 return kFALSE;
233 }
234
235// if (fRunHeader->GetFreqSampling()!=1000)
236// {
237// *fLog << err << "ERROR - Sampling frequencies others than 1GHz are not yet supported." << endl;
238// *fLog << warn << "FIXME - SCALE MPulsShape WITH THE SAMPLING FREQUENCY." << endl;
239// return kFALSE;
240// }
241
242 fSpline = pulse->GetSpline();
243 if (!fSpline)
244 {
245 *fLog << err << "No spline initialized." << endl;
246 return kFALSE;
247 }
248
249 // ---------------- Information output ----------------------
250
251 if (fBaselineGain)
252 *fLog << inf << "Gain is also applied to the electronic noise." << endl;
253
254 return kTRUE;
255}
256
257// --------------------------------------------------------------------------
258//
259// FIXME: For now this is a workaround to set a baseline and the
260// electronic (guassian noise)
261//
262Bool_t MSimCamera::ReInit(MParList *plist)
263{
264 for (int i=0; i<fElectronicNoise->GetSize(); i++)
265 {
266 MPedestalPix &ped = (*fElectronicNoise)[i];
267 ped.SetPedestal(fDefaultOffset);
268 if (fDefaultNoise>0)
269 ped.SetPedestalRms(fDefaultNoise);
270
271 ped.SetPedestalABoffset(0);
272 ped.SetNumEvents(0);
273
274
275 MPedestalPix &gain = (*fGain)[i];
276 if (fDefaultGain>0)
277 gain.SetPedestal(fDefaultGain);
278
279 gain.SetPedestalRms(0);
280 gain.SetPedestalABoffset(0);
281 gain.SetNumEvents(0);
282 }
283
284 return kTRUE;
285}
286
287// --------------------------------------------------------------------------
288//
289// fStat->GetMaxIndex must return the maximum index possible
290// (equiv. number of pixels) not just the maximum index stored!
291//
292Int_t MSimCamera::Process()
293{
294 // Calculate start time, end time and corresponding number of samples
295 const Double_t freq = fRunHeader->GetFreqSampling()/1000.;
296
297 // FIXME: Should we use a higher sampling here?
298
299 const Double_t start = fStat->GetTimeFirst()*freq;
300 const Double_t end = fStat->GetTimeLast() *freq;
301
302 const UInt_t nlen = TMath::CeilNint(end-start);
303
304 // Get number of pixels/channels
305 const UInt_t npix = fStat->GetMaxIndex()+1;
306
307 if (npix>(UInt_t)fElectronicNoise->GetSize())
308 {
309 *fLog << err << "ERROR - More indices (" << npix << ") ";
310 *fLog << "assigned than existing in camera (";
311 *fLog << fElectronicNoise->GetSize() << ")!" << endl;
312 return kERROR;
313 }
314
315 const Double_t pl = fSpline->GetXmin()*freq;
316 const Double_t pr = fSpline->GetXmax()*freq;
317
318 // Init the arrays and set the range which will contain valid data
319 fCamera->Init(npix, nlen);
320 fCamera->SetValidRange(TMath::FloorNint(pr), TMath::CeilNint(nlen+pl));
321
322 Double_t timeoffset[npix];
323
324
325 // Add electronic noise to empty channels
326 for (UInt_t i=0; i<npix; i++)
327 {
328
329 // Get the ResidualTimeSpread Parameter
330 const Double_t residualTimeSpread = fResidualTimeSpread->GetVal();
331
332 // Jens Buss on residual time spread:
333 // randomly draw an additional time offset to be added to the arrivaltime
334 // from a gaussian normal distribution with a given standard deviation
335 timeoffset[i] = gRandom->Gaus(0.0, residualTimeSpread);
336 const MPedestalPix &pix = (*fElectronicNoise)[i];
337
338 const Double_t val = pix.GetPedestal();
339 const Double_t rms = pix.GetPedestalRms();
340
341 // FTemme: Implementation of AC-coupling:
342 // to calculate the value of the accoupling per slice I use the
343 // following equation:
344 // accouplingPerSlice = accidentalPhotonRate * (1 + crossTalkProb)
345 // * areaOfOnePulse / samplingRate;
346 // Therefore I need the following variables
347 // Double_t accidentalPhotonRate; // [MHz]
348 // Float_t crossTalkProb; // [1]
349 // Double_t areaOfOnePulse; // [ADC-Counts * s]
350 // Double_t samplingRate; // [slices * MHz]
351
352 // The accidental photon rate is stored in GHz, so we have to multiply
353 // with 1E3 to get MHz:
354 const MPedestalPix &accPhoPix = (*fAccidentalPhotons)[i];
355
356 const Double_t accidentalPhotonRate = accPhoPix.GetPedestal() * 1e3; //[MHz]
357
358 Double_t currentAccidentalPhotonRate = accidentalPhotonRate;
359 if (fACTimeConstant!=0)
360 {
361 const Double_t accidentalPhotons = fACTimeConstant * accidentalPhotonRate;
362 const Double_t sigmaAccidentalPhotons = TMath::Sqrt(accidentalPhotons);
363
364 const Double_t gaus = gRandom->Gaus(accidentalPhotons,sigmaAccidentalPhotons);
365
366 currentAccidentalPhotonRate = gaus / fACTimeConstant;
367 }
368
369 // Get the CrosstalkCoefficient Parameter
370 const Double_t crossTalkProb = fCrosstalkCoeffParam->GetVal();
371
372 // To get the area of one Pulse, I only need to calculate the Integral
373 // of the Pulse Shape, which is stored in fSpline. Because the spline is
374 // normalized to a maximal amplitude of 1.0, I had to multiply it with
375 // the Default gain [ADC-Counts * s]
376 const Double_t areaOfOnePulse = fSpline->Integral() * fDefaultGain;
377
378 // The sampling rate I get from the RunHeader:
379 const Double_t samplingRate = fRunHeader->GetFreqSampling(); // [slices * MHz]
380
381 const Double_t accouplingPerSlice = currentAccidentalPhotonRate
382 * (1 + crossTalkProb + fACFudgeFactor)
383 * areaOfOnePulse / samplingRate;
384
385 // The accoupling is substracted from the timeline by decreasing the
386 // mean of the gaussian noise which is added
387
388 if (!fBaselineGain)
389 {
390 (*fCamera)[i].AddGaussianNoise(rms, val - accouplingPerSlice);
391 continue;
392 }
393 // Sorry, the name "pedestal" is misleading here
394 // FIXME: Simulate gain fluctuations
395 const Double_t gain = (*fGain)[i].GetPedestal();
396
397 // FIXME: We might add the base line here already.
398 // FIXME: How stable is the offset?
399 // FIXME: Should we write a container AppliedGain for MSImTrigger?
400
401 (*fCamera)[i].AddGaussianNoise(rms*gain, (val - accouplingPerSlice)*gain);
402 }
403
404 // FIXME: Simulate correlations with neighboring pixels
405
406 const Int_t num = fEvt->GetNumPhotons();
407
408 // A random shift, uniformely distributed within one slice, to make sure that
409 // the first photon is not always aligned identically with a sample edge.
410 // FIXME: Make it switchable
411 const Float_t rndm = gRandom->Uniform();
412
413 // FIXME: Shell we add a random shift of [0,1] samples per channel?
414 // Or maybe per channel and run?
415
416 Double_t tot = 0;
417
418 for (int i=0 ; i<1440 ; i++)
419 {
420 (*fTruePhotons->cherenkov_photons_weight)[i] = 0;
421 (*fTruePhotons->cherenkov_photons_number)[i] = 0;
422 (*fTruePhotons->cherenkov_arrival_time_mean)[i] = 0;
423 (*fTruePhotons->cherenkov_arrival_time_variance)[i] = 0;
424 (*fTruePhotons->muon_cherenkov_photons_weight)[i] = 0;
425 (*fTruePhotons->muon_cherenkov_photons_number)[i] = 0;
426 (*fTruePhotons->cherenkov_arrival_time_min)[i] = 10000;
427 (*fTruePhotons->cherenkov_arrival_time_max)[i] = 0;
428 (*fTruePhotons->noise_photons_weight)[i] = 0;
429 }
430
431 //--------------------------------------------------------------------------
432
433 // Get the ResidualTimeSpread Parameter
434 const Double_t gapdTimeJitter = fGapdTimeJitter->GetVal();
435
436 // Simulate pulses
437 for (Int_t i=0; i<num; i++)
438 {
439 const MPhotonData &ph = (*fEvt)[i];
440
441 const UInt_t idx = ph.GetTag();
442 Double_t t = (ph.GetTime()-fStat->GetTimeFirst())*freq+rndm;// - fSpline->GetXmin();
443
444 // Sebastian Mueller and Dominik Neise on fix time offsets:
445 // We add a fix temporal offset to the relative arrival time of the
446 // individual pixel. The offsets are stored in the
447 // fFixTimeOffsetsBetweenPixelsInNs -> fM matrix. We identify the first
448 // column to hold the offsets in ns.
449 t = t + freq*fFixTimeOffsetsBetweenPixelsInNs->fM[idx][0];
450
451 // Jens Buss on residual time spread:
452 // add random time offset to the arrivaltimes
453 t = t + timeoffset[idx];
454
455 // FIXME: Time jitter?
456 // Jens Buss on GapdTimeJitter
457 // add also a time offset to arrival times of single photons
458 // TODO: change to ns, use: fRunHeader->GetFreqSampling()
459 Double_t timeJitter = gRandom->Gaus(0.0, gapdTimeJitter);
460 t = t + timeJitter;
461
462 // FIXME: Add additional routing here?
463 // FIMXE: How stable is the gain?
464
465 if (ph.GetPrimary()!=MMcEvt::kNightSky && ph.GetPrimary()!=MMcEvt::kArtificial)
466 {
467 tot += ph.GetWeight();
468 (*fTruePhotons->cherenkov_photons_weight)[idx] += ph.GetWeight();
469 (*fTruePhotons->cherenkov_photons_number)[idx] += 1;
470
471 (*fTruePhotons->cherenkov_arrival_time_mean)[idx] += t;
472 (*fTruePhotons->cherenkov_arrival_time_variance)[idx] += t*t;
473
474 if (ph.GetPrimary()==MMcEvt::kMUON)
475 {
476 (*fTruePhotons->muon_cherenkov_photons_weight)[idx] += ph.GetWeight();
477 (*fTruePhotons->muon_cherenkov_photons_number)[idx] += 1;
478 }
479
480 // find min
481 if (t < (*fTruePhotons->cherenkov_arrival_time_min)[idx] )
482 {
483 (*fTruePhotons->cherenkov_arrival_time_min)[idx] = t;
484 }
485 // find max
486 if (t > (*fTruePhotons->cherenkov_arrival_time_max)[idx] )
487 {
488 (*fTruePhotons->cherenkov_arrival_time_max)[idx] = t;
489 }
490 }
491 else
492 {
493 (*fTruePhotons->noise_photons_weight)[idx] += ph.GetWeight();
494 }
495
496 // Sorry, the name "pedestal" is misleading here
497 // FIXME: Simulate gain fluctuations
498 const Double_t gain = (*fGain)[idx].GetPedestal();
499
500 // === FIXME === FIXME === FIXME === Frequency!!!!
501 (*fCamera)[idx].AddPulse(*fSpline, t, ph.GetWeight()*gain);
502 }
503
504 for (unsigned int i=0 ; i < 1440 ; i++)
505 {
506 float number = (*fTruePhotons->cherenkov_photons_number)[i];
507 (*fTruePhotons->cherenkov_arrival_time_mean)[i] /= number;
508 float mean = (*fTruePhotons->cherenkov_arrival_time_mean)[i];
509 float sum_tt = (*fTruePhotons->cherenkov_arrival_time_variance)[i];
510 (*fTruePhotons->cherenkov_arrival_time_variance)[i] = (sum_tt / number - mean*mean) /(number - 1);
511 }
512
513 fMcEvt->SetPhotElfromShower(TMath::Nint(tot));
514
515 return kTRUE;
516}
517
518// --------------------------------------------------------------------------
519//
520// BaselineGain: Off
521//
522Int_t MSimCamera::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
523{
524 Bool_t rc = kFALSE;
525 if (IsEnvDefined(env, prefix, "BaselineGain", print))
526 {
527 rc = kTRUE;
528 fBaselineGain = GetEnvValue(env, prefix, "BaselineGain", fBaselineGain);
529 }
530
531 if (IsEnvDefined(env, prefix, "DefaultOffset", print))
532 {
533 rc = kTRUE;
534 fDefaultOffset = GetEnvValue(env, prefix, "DefaultOffset", fDefaultOffset);
535 }
536 if (IsEnvDefined(env, prefix, "DefaultNoise", print))
537 {
538 rc = kTRUE;
539 fDefaultNoise = GetEnvValue(env, prefix, "DefaultNoise", fDefaultNoise);
540 }
541 if (IsEnvDefined(env, prefix, "DefaultGain", print))
542 {
543 rc = kTRUE;
544 fDefaultGain = GetEnvValue(env, prefix, "DefaultGain", fDefaultGain);
545 }
546 if (IsEnvDefined(env, prefix, "ACFudgeFactor", print))
547 {
548 rc = kTRUE;
549 fACFudgeFactor = GetEnvValue(env, prefix, "ACFudgeFactor", fACFudgeFactor);
550 }
551 if (IsEnvDefined(env, prefix, "ACTimeConstant", print))
552 {
553 rc = kTRUE;
554 fACTimeConstant = GetEnvValue(env, prefix, "ACTimeConstant", fACTimeConstant);
555 }
556
557 return rc;
558}
Note: See TracBrowser for help on using the repository browser.