source: branches/MarsResidualTimeSpread/msimcamera/MSimCamera.cc@ 18038

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