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

Last change on this file since 13417 was 10093, checked in by tbretz, 14 years ago
Added afterpulse treatment to MAvalanchePhotoDiode and MSimAPD
File size: 11.3 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{
83 fName = name ? name : "MSimAPD";
84 fTitle = title ? title : " Task to simulate the detection behaviour of APDs";
85}
86
87// --------------------------------------------------------------------------
88//
89// Get the necessary parameter containers
90//
91Int_t MSimAPD::PreProcess(MParList *pList)
92{
93 if (fNameGeomCam.IsNull())
94 {
95 *fLog << inf << "No geometry container... skipping." << endl;
96 return kSKIP;
97 }
98
99 fGeom = (MGeomCam*)pList->FindObject(fNameGeomCam, "MGeomCam");
100 if (!fGeom)
101 {
102 *fLog << inf << fNameGeomCam << " [MGeomCam] not found..." << endl;
103
104 fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
105 if (!fGeom)
106 {
107 *fLog << err << "MGeomCam not found... aborting." << endl;
108 return kFALSE;
109 }
110 }
111
112 fStat = (MPhotonStatistics*)pList->FindObject("MPhotonStatistics");
113 if (!fStat)
114 {
115 *fLog << err << "MPhotonStatistics not found... aborting." << endl;
116 return kFALSE;
117 }
118
119 fEvt = (MPhotonEvent*)pList->FindObject("MPhotonEvent");
120 if (!fEvt)
121 {
122 *fLog << err << "MPhotonEvent not found... aborting." << endl;
123 return kFALSE;
124 }
125
126 fRates = (MPedestalCam*)pList->FindObject("AccidentalPhotonRates", "MPedestalCam");
127 if (!fRates)
128 {
129 *fLog << inf;
130 *fLog << "AccidentalPhotonRates [MPedestalCam] not found..." << endl;
131 *fLog << " using " << fFreq << " as default for all G-APDs." << endl;
132 }
133
134 return kTRUE;
135}
136
137// --------------------------------------------------------------------------
138//
139// Initialize as many APDs as we have pixels in the fGeomCam
140//
141Bool_t MSimAPD::ReInit(MParList *plist)
142{
143 if (UInt_t(fAPDs.GetEntriesFast())==fGeom->GetNumPixels())
144 return kTRUE;
145
146 fAPDs.Delete();
147
148 // FIXME:
149 // * initialize an empty APD and read the APD setup from a file to
150 // allow different APDs.
151 // * Make the arguments a data member of MSimAPD
152
153 Int_t ncells = 0;
154 Float_t crosstalk = 0;
155 Float_t deadtime = 0;
156 Float_t recovery = 0;
157 Float_t afterprob1 = 0;
158 Float_t afterprob2 = 0;
159
160 switch (fType)
161 {
162 case 1:
163 ncells = 30;
164 crosstalk = 0.2;
165 deadtime = 3;
166 recovery = 8.75*4;
167 afterprob1 = 0;
168 afterprob2 = 0;
169 break;
170
171 case 2:
172 ncells = 60;
173 crosstalk = 0.2;
174 deadtime = 3;
175 recovery = 8.75;
176 afterprob1 = 0;
177 afterprob2 = 0;
178 break;
179
180 case 3:
181 ncells = 60;
182 crosstalk = 0.15;
183 deadtime = 3;
184 recovery = 8.75;
185 afterprob1 = 0;
186 afterprob2 = 0;
187 break;
188
189 case 4:
190 ncells = 60;
191 crosstalk = 0.15;
192 deadtime = 3;
193 recovery = 8.75;
194 afterprob1 = 0.14;
195 afterprob2 = 0.11;
196 break;
197
198 default:
199 *fLog << err << "ERROR - APD type " << fType << " undefined." << endl;
200 return kFALSE;
201 }
202
203 for (UInt_t i=0; i<fGeom->GetNumPixels(); i++)
204 {
205 APD *apd = new APD(ncells, crosstalk, deadtime, recovery);
206 apd->SetAfterpulseProb(afterprob1, afterprob2);
207
208 fAPDs.Add(apd);
209 }
210
211 return kTRUE;
212}
213
214// --------------------------------------------------------------------------
215//
216// Process all photons through the corresponding APD and set the output
217// (weight) accordingly.
218//
219Int_t MSimAPD::Process()
220{
221 // Make all APDs look neutral for the first hit by a photon according to the
222 // average hit rate
223 const UInt_t npix = fAPDs.GetEntriesFast();
224
225 // Check if we can safely proceed (this can fail if we either haven't been
226 // ReInit'ed or the max index in MPhotonStatistics is wrong)
227 if ((Int_t)npix<fStat->GetMaxIndex())
228 {
229 *fLog << err << "ERROR - MSimAPD::Process: Only " << npix << " APDs initialized. At least " << fStat->GetMaxIndex() << " needed... abort." << endl;
230 return kERROR;
231 }
232
233 // To ensure that the photons are always sorted when calling
234 // HitRandomCellRelative the array is sorted here. If it is sorted
235 // already nothing will be done since the status is stored.
236 // FIXME: Check that this is true and check that it is really necessary
237 fEvt->Sort();
238
239 // This tries to initialize dead and relaxing cells properly. If
240 // the APD has not been initialized before the chip is randomsly
241 // filled, otherwise a time window of the default relaxing time
242 // is simulated, so that the previous influence is less than a permille.
243 for (UInt_t idx=0; idx<npix; idx++)
244 {
245 const Double_t freq = fRates ? (*fRates)[idx].GetPedestal() : fFreq;
246
247 // Init creates an empty G-APD, i.e. without external pulses
248 // but the correct history for afterpulses and relaxation.
249 // If it was already initialized once it evolves the G-APD
250 // for a time until the effect of relaxation and afterpulses
251 // is below 0.1%. The also creates the possible afterpulses
252 // of the future and deletes later afterpulses from the list.
253 // After the the time stamp fTime is set to 0.
254 static_cast<APD*>(fAPDs.UncheckedAt(idx))->Init(freq);
255 }
256
257 // Get number of photons
258 const Int_t num = fEvt->GetNumPhotons();
259
260 // Loop over all photons
261 for (Int_t i=0; i<num; i++)
262 {
263 // Get i-th photon
264 MPhotonData &ph = (*fEvt)[i];
265
266 // Get arrival time of photon wrt to left edge of window and its index
267 const Double_t t = ph.GetTime()-fStat->GetTimeFirst();
268 const Int_t idx = ph.GetTag();
269 if (idx<0)
270 {
271 *fLog << err << "ERROR - MSimAPD: Invalid index -1." << endl;
272 return kERROR;
273 }
274
275 if (ph.GetWeight()!=1)
276 {
277 *fLog << err << "ERROR - MSimAPD: Weight of " << i << "-th photon not 1, but " << ph.GetWeight() << endl;
278 ph.Print();
279 return kERROR;
280 }
281
282 // Simulate hitting the APD at a time t after T0 (APD::fTime).
283 // Crosstalk is taken into account and the resulting signal height
284 // in effective "number of photons" is returned. Afterpulses until
285 // this time "hit" the G-APD and newly created afterpulses
286 // are stored in the list of afterpulses
287 const Double_t hits = static_cast<APD*>(fAPDs.UncheckedAt(idx))->HitRandomCellRelative(t);
288
289 // Set the weight to the input
290 ph.SetWeight(hits);
291 }
292
293 // Now we have to shift the evolved time of all APDs to the end of our
294 // simulated time.
295 for (UInt_t idx=0; idx<npix; idx++)
296 {
297 APD *a = static_cast<APD*>(fAPDs.UncheckedAt(idx));
298
299 const Double_t end = fStat->GetTimeLast()-fStat->GetTimeFirst();
300
301 // This moves T0 (APD::fTime) at the right edge of our time-
302 // window. For the next event this means that afterpulses of past
303 // noise and afterpulses will be available already.
304 // FIXME: Note, that this might mean that a cosmics-pulse
305 // might increase the noise above the normal level.
306 a->IncreaseTime(end);
307
308 // Get the afterpulses and add them to the signal
309 TIter Next(&a->GetListOfAfterpulses());
310 Afterpulse *ap = 0;
311 while ((ap=static_cast<Afterpulse*>(Next())))
312 {
313 // Skip afterpulses later than that which have been
314 // already produced
315 if (ap->GetTime()>=end)
316 continue;
317
318 // Add a new photon
319 // FIXME: SLOW!
320 MPhotonData &ph = fEvt->Add();
321
322 // Set source to Artificial (noise), the amplitude produced by the
323 // afterpulse (includes crosstalk), its arrival time
324 // and its amplitude, as well as the index of the channel it
325 // corresponds to.
326 ph.SetPrimary(MMcEvtBasic::kArtificial);
327 ph.SetWeight(ap->GetAmplitude());
328 ph.SetTime(ap->GetTime()+fStat->GetTimeFirst());
329 ph.SetTag(idx);
330 }
331
332 // It seems to make sense to delete the previous afterpulses now
333 // but this is not necessary. We have to loop over them in any
334 // case. So we omit the additional loop for deletion but instead
335 // do the deletion in the next loop at the end of Init()
336 // If needed this could be used to keep the total memory
337 // consumption slighly lower.
338 }
339
340 // Now the newly added afterpulses have to be sorted into the array correctly
341 fEvt->Sort();
342
343 return kTRUE;
344}
345
346// --------------------------------------------------------------------------
347//
348// NameGeomCam
349// Type: 1
350//
351Int_t MSimAPD::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
352{
353 Bool_t rc = kFALSE;
354 if (IsEnvDefined(env, prefix, "NameGeomCam", print))
355 {
356 rc = kTRUE;
357 fNameGeomCam = GetEnvValue(env, prefix, "NameGeomCam", fNameGeomCam);
358 }
359
360 if (IsEnvDefined(env, prefix, "Type", print))
361 {
362 rc = kTRUE;
363 fType = GetEnvValue(env, prefix, "Type", fType);
364 }
365
366 return rc;
367}
Note: See TracBrowser for help on using the repository browser.