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

Last change on this file since 19623 was 19623, checked in by tbretz, 12 months ago
Sanity check for a bad programming style.
File size: 19.8 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(size_t row_index=0; row_index<fFixTimeOffsetsBetweenPixelsInNs->fM.size(); row_index++){
158        const vector<double> row = fFixTimeOffsetsBetweenPixelsInNs->fM.at(row_index);
159        for(size_t col_index=0; col_index<row.size(); col_index++){
160            const double specific_delay = row.at(col_index);
161            if (!TMath::Finite(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        return kFALSE;
217
218    fTruePhotons = (MTruePhotonsPerPixelCont*)pList->FindCreateObj("MTruePhotonsPerPixelCont");
219    if (!fTruePhotons)
220        return kFALSE;
221
222    MParSpline *pulse = (MParSpline*)pList->FindObject("PulseShape", "MParSpline");
223    if (!pulse)
224    {
225        *fLog << err << "PulseShape [MParSpline] not found... aborting." << endl;
226        return kFALSE;
227    }
228
229//    if (fRunHeader->GetFreqSampling()!=1000)
230//    {
231//        *fLog << err  << "ERROR - Sampling frequencies others than 1GHz are not yet supported." << endl;
232//        *fLog << warn << "FIXME - SCALE MPulsShape WITH THE SAMPLING FREQUENCY." << endl;
233//        return kFALSE;
234//    }
235
236    fSpline = pulse->GetSpline();
237    if (!fSpline)
238    {
239        *fLog << err << "No spline initialized." << endl;
240        return kFALSE;
241    }
242
243    // ---------------- Information output ----------------------
244
245    if (fBaselineGain)
246        *fLog << inf << "Gain is also applied to the electronic noise." << endl;
247
248    // To get the area of one Pulse, I only need to calculate the Integral
249    // of the Pulse Shape, which is stored in fSpline. Because the spline is
250    // normalized to a maximal amplitude of 1.0, I had to multiply it with
251    // the Default gain [ADC-Counts * s]
252    fAreaOfOnePulse = fSpline->Integral() * fDefaultGain;
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        // The sampling rate I get from the RunHeader:
373        const Double_t samplingRate = fRunHeader->GetFreqSampling(); // [slices * MHz]
374
375        const Double_t accouplingPerSlice = currentAccidentalPhotonRate
376            * (1 + crossTalkProb + fACFudgeFactor)
377            * fAreaOfOnePulse / samplingRate;
378
379        // The accoupling is substracted from the timeline by decreasing the
380        // mean of the gaussian noise which is added
381
382        if (!fBaselineGain)
383        {
384            (*fCamera)[i].AddGaussianNoise(rms, val - accouplingPerSlice);
385            continue;
386        }
387        // Sorry, the name "pedestal" is misleading here
388        // FIXME: Simulate gain fluctuations
389        const Double_t gain = (*fGain)[i].GetPedestal();
390
391        // FIXME: We might add the base line here already.
392        // FIXME: How stable is the offset?
393        // FIXME: Should we write a container AppliedGain for MSImTrigger?
394
395        (*fCamera)[i].AddGaussianNoise(rms*gain, (val - accouplingPerSlice)*gain);
396    }
397
398    // FIXME: Simulate correlations with neighboring pixels
399
400    const Int_t num = fEvt->GetNumPhotons();
401
402    // A random shift, uniformely distributed within one slice, to make sure that
403    // the first photon is not always aligned identically with a sample edge.
404    // FIXME: Make it switchable
405    const Float_t rndm = gRandom->Uniform();
406
407    // FIXME: Shell we add a random shift of [0,1] samples per channel?
408    //        Or maybe per channel and run?
409
410    Double_t tot = 0;
411
412    // Sanity check for bad programming style
413    if (npix==1440)
414    {
415        for (int i=0 ; i<1440 ; i++)
416        {
417            (*fTruePhotons->cherenkov_photons_weight)[i] = 0;
418            (*fTruePhotons->cherenkov_photons_number)[i] = 0;
419            (*fTruePhotons->cherenkov_arrival_time_mean)[i] = 0;
420            (*fTruePhotons->cherenkov_arrival_time_variance)[i] = 0;
421            (*fTruePhotons->muon_cherenkov_photons_weight)[i] = 0;
422            (*fTruePhotons->muon_cherenkov_photons_number)[i] = 0;
423            (*fTruePhotons->cherenkov_arrival_time_min)[i] = 10000;
424            (*fTruePhotons->cherenkov_arrival_time_max)[i] = 0;
425            (*fTruePhotons->noise_photons_weight)[i] = 0;
426        }
427    }
428
429    //--------------------------------------------------------------------------
430
431    // Get the ResidualTimeSpread Parameter
432    const Double_t gapdTimeJitter = fGapdTimeJitter->GetVal();
433   
434    // Simulate pulses
435    for (Int_t i=0; i<num; i++)
436    {
437        const MPhotonData &ph = (*fEvt)[i];
438
439        const UInt_t   idx = ph.GetTag();
440        Double_t t   = (ph.GetTime()-fStat->GetTimeFirst())*freq+rndm;// - fSpline->GetXmin();
441
442        // Sebastian Mueller and Dominik Neise on fix time offsets:
443        // We add a fix temporal offset to the relative arrival time of the
444        // individual pixel. The offsets are stored in the
445        // fFixTimeOffsetsBetweenPixelsInNs -> fM matrix. We identify the first
446        // column to hold the offsets in ns.
447        t += freq*fFixTimeOffsetsBetweenPixelsInNs->fM[idx][0];
448
449        // Jens Buss on residual time spread:
450        // add random time offset to the arrivaltimes
451        t += timeoffset[idx];
452
453        // FIXME: Time jitter?
454        // Jens Buss on GapdTimeJitter
455        // add also a time offset to arrival times of single photons
456        // TODO: change to ns, use: fRunHeader->GetFreqSampling()
457        t += gRandom->Gaus(0.0, gapdTimeJitter);
458
459        // FIXME: Add additional routing here?
460        // FIMXE: How stable is the gain?
461
462        if (ph.GetPrimary()!=MMcEvt::kNightSky && ph.GetPrimary()!=MMcEvt::kArtificial)
463        {
464            tot += ph.GetWeight();
465
466            // Sanity check for bad programming style
467            if (npix==1440)
468            {
469                (*fTruePhotons->cherenkov_photons_weight)[idx] += ph.GetWeight();
470                (*fTruePhotons->cherenkov_photons_number)[idx] += 1;
471
472                (*fTruePhotons->cherenkov_arrival_time_mean)[idx] += t;
473                (*fTruePhotons->cherenkov_arrival_time_variance)[idx] += t*t;
474
475                if (ph.GetPrimary()==MMcEvt::kMUON)
476                {
477                    (*fTruePhotons->muon_cherenkov_photons_weight)[idx] += ph.GetWeight();
478                    (*fTruePhotons->muon_cherenkov_photons_number)[idx] += 1;
479                }
480
481                // find min
482                if (t < (*fTruePhotons->cherenkov_arrival_time_min)[idx] )
483                {
484                    (*fTruePhotons->cherenkov_arrival_time_min)[idx] = t;
485                }
486                // find max
487                if (t > (*fTruePhotons->cherenkov_arrival_time_max)[idx] )
488                {
489                    (*fTruePhotons->cherenkov_arrival_time_max)[idx] = t;
490                }
491            }
492        }
493        else
494        {
495            // Sanity check for bad programming style
496            if (npix==1440)
497            {
498                (*fTruePhotons->noise_photons_weight)[idx] += ph.GetWeight();
499            }
500        }
501
502        // Sorry, the name "pedestal" is misleading here
503        // FIXME: Simulate gain fluctuations
504        const Double_t gain = (*fGain)[idx].GetPedestal();
505
506        // === FIXME === FIXME === FIXME === Frequency!!!!
507        (*fCamera)[idx].AddPulse(*fSpline, t, ph.GetWeight()*gain);
508    }
509
510    // Sanity check for bad programming style
511    if (npix==1440)
512    {
513        for (unsigned int i=0 ; i < 1440 ; i++)
514        {
515            float number = (*fTruePhotons->cherenkov_photons_number)[i];
516            (*fTruePhotons->cherenkov_arrival_time_mean)[i] /= number;
517            float mean = (*fTruePhotons->cherenkov_arrival_time_mean)[i];
518            float sum_tt = (*fTruePhotons->cherenkov_arrival_time_variance)[i];
519            (*fTruePhotons->cherenkov_arrival_time_variance)[i] = (sum_tt / number - mean*mean) /(number - 1);
520        }
521    }
522
523    fMcEvt->SetPhotElfromShower(TMath::Nint(tot));
524
525    return kTRUE;
526}
527
528// --------------------------------------------------------------------------
529//
530// BaselineGain: Off
531//
532Int_t MSimCamera::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
533{
534    Bool_t rc = kFALSE;
535    if (IsEnvDefined(env, prefix, "BaselineGain", print))
536    {
537        rc = kTRUE;
538        fBaselineGain = GetEnvValue(env, prefix, "BaselineGain", fBaselineGain);
539    }
540
541    if (IsEnvDefined(env, prefix, "DefaultOffset", print))
542    {
543        rc = kTRUE;
544        fDefaultOffset = GetEnvValue(env, prefix, "DefaultOffset", fDefaultOffset);
545    }
546    if (IsEnvDefined(env, prefix, "DefaultNoise", print))
547    {
548        rc = kTRUE;
549        fDefaultNoise = GetEnvValue(env, prefix, "DefaultNoise", fDefaultNoise);
550    }
551    if (IsEnvDefined(env, prefix, "DefaultGain", print))
552    {
553        rc = kTRUE;
554        fDefaultGain = GetEnvValue(env, prefix, "DefaultGain", fDefaultGain);
555    }
556    if (IsEnvDefined(env, prefix, "ACFudgeFactor", print))
557    {
558        rc = kTRUE;
559        fACFudgeFactor = GetEnvValue(env, prefix, "ACFudgeFactor", fACFudgeFactor);
560    }
561    if (IsEnvDefined(env, prefix, "ACTimeConstant", print))
562    {
563        rc = kTRUE;
564        fACTimeConstant = GetEnvValue(env, prefix, "ACTimeConstant", fACTimeConstant);
565    }
566
567    return rc;
568}
Note: See TracBrowser for help on using the repository browser.