source: trunk/Mars/msimcamera/MSimAPD.cc@ 17248

Last change on this file since 17248 was 17202, checked in by dneise, 11 years ago
in commit 17200 & 17201 I committed accidentally a lot of changes, that were totally unnecessary. So I reverted back to 17199. In this revision the svn property svn:ignore is still set, such that all the Cint.cc and Cint.h files as well as .d files are ignored. In addition I changed melectronics/MAvalanchePhotoDiode.h: I gave the Afterpulse class a ClassDef macro... and I added the link pragma into LinkDef.h file. So now one can easier write test macros for the APD class, without the need to compile the macro. Sorry for the hubbub.
File size: 13.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@astro.uni-wuerzburg.de>
19!
20! Copyright: CheObs Software Development, 2000-2009
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MSimAPD
28//
29// This tasks simulates the individual APDs. Before starting the APD is
30// initialized randomly according to the photon rate hitting the APD. Such
31// it is assumed that the initial condition of the APD is similar to the real
32// one. In this context it is assumed that the events are independent, so
33// that the APD is always in the same condition.
34//
35// For every photon and event the behaviour of the APD is simulated. The
36// output is set as weight to the MPhotonData containers.
37//
38// Remark:
39// - The photon rate used to initialize the APD must match the one used
40// to "fill" the random photons. (FIXME: This should be stored somewhere)
41//
42// Input Containers:
43// fNameGeomCam [MGeomCam]
44// MPhotonEvent
45// MPhotonStatistics
46//
47// Output Containers:
48// MPhotonEvent
49//
50//////////////////////////////////////////////////////////////////////////////
51#include "MSimAPD.h"
52
53#include <TH2.h>
54#include <TRandom.h>
55
56#include "MLog.h"
57#include "MLogManip.h"
58
59#include "MMath.h"
60#include "MParList.h"
61
62#include "MGeomCam.h"
63
64#include "MPhotonEvent.h"
65#include "MPhotonData.h"
66
67#include "MPedestalCam.h"
68#include "MPedestalPix.h"
69
70#include "MAvalanchePhotoDiode.h"
71
72ClassImp(MSimAPD);
73
74using namespace std;
75
76// --------------------------------------------------------------------------
77//
78// Default Constructor.
79//
80MSimAPD::MSimAPD(const char* name, const char *title)
81 : fGeom(0), fEvt(0), fStat(0), fType(1),
82 fNumCells(60), fCrosstalkCoeff(0), fDeadTime(3),
83 fRecoveryTime(8.75), fAfterpulseProb1(0.11), fAfterpulseProb2(0.14)
84
85{
86 fName = name ? name : "MSimAPD";
87 fTitle = title ? title : " Task to simulate the detection behaviour of APDs";
88}
89
90// --------------------------------------------------------------------------
91//
92// Get the necessary parameter containers
93//
94Int_t MSimAPD::PreProcess(MParList *pList)
95{
96 if (fNameGeomCam.IsNull())
97 {
98 *fLog << inf << "No geometry container... skipping." << endl;
99 return kSKIP;
100 }
101
102 fGeom = (MGeomCam*)pList->FindObject(fNameGeomCam, "MGeomCam");
103 if (!fGeom)
104 {
105 *fLog << inf << fNameGeomCam << " [MGeomCam] not found..." << endl;
106
107 fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
108 if (!fGeom)
109 {
110 *fLog << err << "MGeomCam not found... aborting." << endl;
111 return kFALSE;
112 }
113 }
114
115 fStat = (MPhotonStatistics*)pList->FindObject("MPhotonStatistics");
116 if (!fStat)
117 {
118 *fLog << err << "MPhotonStatistics not found... aborting." << endl;
119 return kFALSE;
120 }
121
122 fEvt = (MPhotonEvent*)pList->FindObject("MPhotonEvent");
123 if (!fEvt)
124 {
125 *fLog << err << "MPhotonEvent not found... aborting." << endl;
126 return kFALSE;
127 }
128
129 fRates = (MPedestalCam*)pList->FindObject("AccidentalPhotonRates", "MPedestalCam");
130 if (!fRates)
131 {
132 *fLog << inf;
133 *fLog << "AccidentalPhotonRates [MPedestalCam] not found..." << endl;
134 *fLog << " using " << fFreq << " as default for all G-APDs." << endl;
135 }
136
137 return kTRUE;
138}
139
140// --------------------------------------------------------------------------
141//
142// Initialize as many APDs as we have pixels in the fGeomCam
143//
144Bool_t MSimAPD::ReInit(MParList *plist)
145{
146 if (UInt_t(fAPDs.GetEntriesFast())==fGeom->GetNumPixels())
147 return kTRUE;
148
149 fAPDs.Delete();
150
151 // FIXME:
152 // * initialize an empty APD and read the APD setup from a file to
153 // allow different APDs.
154 // * Make the arguments a data member of MSimAPD
155
156 Int_t ncells = 0;
157 Float_t crosstalk = 0;
158 Float_t deadtime = 0;
159 Float_t recovery = 0;
160 Float_t afterprob1 = 0;
161 Float_t afterprob2 = 0;
162
163 switch (fType)
164 {
165 case 0:
166 ncells = fNumCells;
167 crosstalk = fCrosstalkCoeff;
168 deadtime = fDeadTime;
169 recovery = fRecoveryTime;
170 afterprob1 = fAfterpulseProb1;
171 afterprob2 = fAfterpulseProb2;
172
173 gLog << inf << "Using custom G-APD parameters:" << endl;
174 gLog << " - Num cells: " << fNumCells << endl;
175 gLog << " - Crosstalk coeff.: " << fCrosstalkCoeff << endl;
176 gLog << " - Dead time: " << fDeadTime << " ns" << endl;
177 gLog << " - Recovery time: " << fRecoveryTime << " ns" << endl;
178 gLog << " - Afterpulse 1 prob.: " << fAfterpulseProb1 << endl;
179 gLog << " - Afterpulse 2 prob.: " << fAfterpulseProb2 << endl;
180 break;
181
182 case 1:
183 ncells = 30;
184 crosstalk = 0.2;
185 deadtime = 3;
186 recovery = 8.75*4;
187 afterprob1 = 0;
188 afterprob2 = 0;
189 break;
190
191 case 2:
192 ncells = 60;
193 crosstalk = 0.2;
194 deadtime = 3;
195 recovery = 8.75;
196 afterprob1 = 0;
197 afterprob2 = 0;
198 break;
199
200 case 3:
201 ncells = 60;
202 crosstalk = 0.15;
203 deadtime = 3;
204 recovery = 8.75;
205 afterprob1 = 0;
206 afterprob2 = 0;
207 break;
208
209 case 4:
210 ncells = 60;
211 crosstalk = 0.15;
212 deadtime = 3;
213 recovery = 8.75;
214 afterprob1 = 0.14;
215 afterprob2 = 0.11;
216 break;
217
218 default:
219 *fLog << err << "ERROR - APD type " << fType << " undefined." << endl;
220 return kFALSE;
221 }
222
223 for (UInt_t i=0; i<fGeom->GetNumPixels(); i++)
224 {
225 APD *apd = new APD(ncells, crosstalk, deadtime, recovery);
226 apd->SetAfterpulseProb(afterprob1, afterprob2);
227
228 fAPDs.Add(apd);
229 }
230
231 return kTRUE;
232}
233
234// --------------------------------------------------------------------------
235//
236// Process all photons through the corresponding APD and set the output
237// (weight) accordingly.
238//
239Int_t MSimAPD::Process()
240{
241 // Make all APDs look neutral for the first hit by a photon according to the
242 // average hit rate
243 const UInt_t npix = fAPDs.GetEntriesFast();
244
245 // Check if we can safely proceed (this can fail if we either haven't been
246 // ReInit'ed or the max index in MPhotonStatistics is wrong)
247 if ((Int_t)npix<fStat->GetMaxIndex())
248 {
249 *fLog << err << "ERROR - MSimAPD::Process: Only " << npix << " APDs initialized. At least " << fStat->GetMaxIndex() << " needed... abort." << endl;
250 return kERROR;
251 }
252
253 // To ensure that the photons are always sorted when calling
254 // HitRandomCellRelative the array is sorted here. If it is sorted
255 // already nothing will be done since the status is stored.
256 // FIXME: Check that this is true and check that it is really necessary
257 fEvt->Sort();
258
259 // This tries to initialize dead and relaxing cells properly. If
260 // the APD has not been initialized before the chip is randomsly
261 // filled, otherwise a time window of the default relaxing time
262 // is simulated, so that the previous influence is less than a permille.
263 for (UInt_t idx=0; idx<npix; idx++)
264 {
265 const Double_t freq = fRates ? (*fRates)[idx].GetPedestal() : fFreq;
266
267 // Init creates an empty G-APD, i.e. without external pulses
268 // but the correct history for afterpulses and relaxation.
269 // If it was already initialized once it evolves the G-APD
270 // for a time until the effect of relaxation and afterpulses
271 // is below 0.1%. The also creates the possible afterpulses
272 // of the future and deletes later afterpulses from the list.
273 // After the the time stamp fTime is set to 0.
274 static_cast<APD*>(fAPDs.UncheckedAt(idx))->Init(freq);
275 }
276
277 // Get number of photons
278 const Int_t num = fEvt->GetNumPhotons();
279
280 // Loop over all photons
281 for (Int_t i=0; i<num; i++)
282 {
283 // Get i-th photon
284 MPhotonData &ph = (*fEvt)[i];
285
286 // Get arrival time of photon wrt to left edge of window and its index
287 const Double_t t = ph.GetTime()-fStat->GetTimeFirst();
288 const Int_t idx = ph.GetTag();
289 if (idx<0)
290 {
291 *fLog << err << "ERROR - MSimAPD: Invalid index -1." << endl;
292 return kERROR;
293 }
294
295 if (ph.GetWeight()!=1)
296 {
297 *fLog << err << "ERROR - MSimAPD: Weight of " << i << "-th photon not 1, but " << ph.GetWeight() << endl;
298 ph.Print();
299 return kERROR;
300 }
301
302 // Simulate hitting the APD at a time t after T0 (APD::fTime).
303 // Crosstalk is taken into account and the resulting signal height
304 // in effective "number of photons" is returned. Afterpulses until
305 // this time "hit" the G-APD and newly created afterpulses
306 // are stored in the list of afterpulses
307 const Double_t hits = static_cast<APD*>(fAPDs.UncheckedAt(idx))->HitRandomCellRelative(t);
308
309 // Set the weight to the input
310 ph.SetWeight(hits);
311 }
312
313 // Now we have to shift the evolved time of all APDs to the end of our
314 // simulated time.
315 for (UInt_t idx=0; idx<npix; idx++)
316 {
317 APD *a = static_cast<APD*>(fAPDs.UncheckedAt(idx));
318
319 const Double_t end = fStat->GetTimeLast()-fStat->GetTimeFirst();
320
321 // This moves T0 (APD::fTime) at the right edge of our time-
322 // window. For the next event this means that afterpulses of past
323 // noise and afterpulses will be available already.
324 // FIXME: Note, that this might mean that a cosmics-pulse
325 // might increase the noise above the normal level.
326 a->IncreaseTime(end);
327
328 // Get the afterpulses and add them to the signal
329 TIter Next(&a->GetListOfAfterpulses());
330 Afterpulse *ap = 0;
331 while ((ap=static_cast<Afterpulse*>(Next())))
332 {
333 // Skip afterpulses later than that which have been
334 // already produced
335 if (ap->GetTime()>=end)
336 continue;
337
338 // Add a new photon
339 // FIXME: SLOW!
340 MPhotonData &ph = fEvt->Add();
341
342 // Set source to Artificial (noise), the amplitude produced by the
343 // afterpulse (includes crosstalk), its arrival time
344 // and its amplitude, as well as the index of the channel it
345 // corresponds to.
346 ph.SetPrimary(MMcEvtBasic::kArtificial);
347 ph.SetWeight(ap->GetAmplitude());
348 ph.SetTime(ap->GetTime()+fStat->GetTimeFirst());
349 ph.SetTag(idx);
350 }
351
352 // It seems to make sense to delete the previous afterpulses now
353 // but this is not necessary. We have to loop over them in any
354 // case. So we omit the additional loop for deletion but instead
355 // do the deletion in the next loop at the end of Init()
356 // If needed this could be used to keep the total memory
357 // consumption slighly lower.
358 }
359
360 // Now the newly added afterpulses have to be sorted into the array correctly
361 fEvt->Sort();
362
363 return kTRUE;
364}
365
366// --------------------------------------------------------------------------
367//
368// NameGeomCam
369// Type: 1
370//
371// Example for a custom made G-APD:
372//
373// Type: 0
374// NumCells: 60
375// CrosstalkCoefficient: 0.1 // [prob]
376// DeadTime: 3 // [ns]
377// RecoveryTime: 8.75 // [ns]
378// AfterpulseProb1: 0.14 // [prob]
379// AfterpulseProb2: 0.11 // [prob]
380//
381Int_t MSimAPD::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
382{
383 Bool_t rc = kFALSE;
384 if (IsEnvDefined(env, prefix, "NameGeomCam", print))
385 {
386 rc = kTRUE;
387 fNameGeomCam = GetEnvValue(env, prefix, "NameGeomCam", fNameGeomCam);
388 }
389
390 if (IsEnvDefined(env, prefix, "Type", print))
391 {
392 rc = kTRUE;
393 fType = GetEnvValue(env, prefix, "Type", fType);
394 }
395
396 if (fType!=0)
397 return rc;
398
399 if (IsEnvDefined(env, prefix, "NumCells", print))
400 {
401 rc = kTRUE;
402 fNumCells = GetEnvValue(env, prefix, "NumCells", fNumCells);
403 }
404 if (IsEnvDefined(env, prefix, "CrosstalkCoefficient", print))
405 {
406 rc = kTRUE;
407 fCrosstalkCoeff = GetEnvValue(env, prefix, "CrosstalkCoefficient", fCrosstalkCoeff);
408 }
409 if (IsEnvDefined(env, prefix, "DeadTime", print))
410 {
411 rc = kTRUE;
412 fDeadTime = GetEnvValue(env, prefix, "DeadTime", fDeadTime);
413 }
414 if (IsEnvDefined(env, prefix, "RecoveryTime", print))
415 {
416 rc = kTRUE;
417 fRecoveryTime = GetEnvValue(env, prefix, "RecoveryTime", fRecoveryTime);
418 }
419 if (IsEnvDefined(env, prefix, "AfterpulseProb1", print))
420 {
421 rc = kTRUE;
422 fAfterpulseProb1 = GetEnvValue(env, prefix, "AfterpulseProb1", fAfterpulseProb1);
423 }
424 if (IsEnvDefined(env, prefix, "AfterpulseProb2", print))
425 {
426 rc = kTRUE;
427 fAfterpulseProb2 = GetEnvValue(env, prefix, "AfterpulseProb2", fAfterpulseProb2);
428 }
429
430 return rc;
431}
Note: See TracBrowser for help on using the repository browser.