source: trunk/Mars/melectronics/MAvalanchePhotoDiode.cc@ 17682

Last change on this file since 17682 was 17323, checked in by tbretz, 11 years ago
According to a simply toy MC this reproduces the behaviour of the sensor best and reproduces very well the measured dark count spectra.
File size: 19.0 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// APD
28//
29// All times in this class are relative times. Therefor the unit for the
30// time is not intrinsically fixed. In fact the dead-time and recovery-
31// time given in the constructor must have the same units. This is what
32// defines the unit of the times given in the function and the unit of
33// rates given.
34// For example, if recovery and dead time are given in nanoseconds the
35// all times must be in nanoseconds and rates are given per nanosecond,
36// i.e. GHz.
37//
38// Hamamatsu 30x30 cells: APD(30, 0.2, 3, 35)
39// Hamamatsu 60x60 cells: APD(60, 0.2, 3, 8.75)
40//
41//
42// The implementation of afterpulsing is based on
43// A.Du, F.Retiere
44// After-pulsing and cross-talk in multi-pixel photon counters
45// NIM A, Volume 596, Issue 3, p. 396-401
46//
47//
48// Example:
49//
50// APD apd(ncells, crosstalk, deadtime, recovery);
51// apd.SetAfterpulseProb(0.14, 0.11);
52//
53// while (1)
54// {
55// // Make the chip "empty" from the influence of external photons
56// // It also sets fTime to 0.
57// apd.Init(freq); // freq of external noise (eg. nsb)
58//
59// // Now call this function for each external photon you have. The
60// // times are relative to the the time you get by apd.GetTime()
61// // set automatically after the call to apd.Init().
62// for (int i=0; i<nphot; i++)
63// apd.HitRandomCellRelative(dt);
64//
65// [...]
66//
67// // Process and produce afterpulses until a time, wrt to GetTime(),
68// // given by the end of your simulated window. Note that this
69// // does not produce noise. This must have been properly simulated
70// // up to this time already.
71// apd.IncreaseTime(dtend);
72//
73// // Now you can excess the afterpulses by
74// TIter Next(&a->GetListOfAfterpulses());
75// Afterpulse *ap = 0;
76// while ((ap=static_cast<Afterpulse*>(Next())))
77// {
78// if (apd.GetTime()>=dtend)
79// continue;
80//
81// cout << "Amplitude: " << ap->GetAmplitude() << endl;
82// cout << "Arrival Time: " << ap->GetTime() << endl;
83// }
84// }
85//
86//
87//////////////////////////////////////////////////////////////////////////////
88#include "MAvalanchePhotoDiode.h"
89
90#include <TRandom.h>
91
92#include "MMath.h"
93
94#include "MLog.h"
95#include "MLogManip.h"
96
97ClassImp(APD);
98
99using namespace std;
100
101/*
102class MyProfile : public TProfile2D
103{
104public:
105 void AddBinEntry(Int_t cell) { fBinEntries.fArray[cell]++; }
106};
107*/
108
109// --------------------------------------------------------------------------
110//
111// Default Constructor.
112//
113// n is the number od cells in x or y. The APD is assumed to
114// be square.
115// prob is the crosstalk probability, i.e., the probability that a
116// photon which produced an avalanche will create another
117// photon in a neighboring cell
118// dt is the deadtime, i.e., the time in which the APD cell will show
119// no response to a photon after a hit
120// rt is the recovering tims, i.e. the exponential (e^(-dt/rt))
121// with which the cell is recovering after being dead
122//
123// prob, dt and ar can be set to 0 to switch the effect off.
124// 0 is also the dfeault for all three.
125//
126APD::APD(Int_t n, Float_t prob, Float_t dt, Float_t rt)
127 : fHist("APD", "", n, 0.5, n+0.5, n, 0.5, n+0.5),
128 fCrosstalkProb(prob), fDeadTime(dt), fRecoveryTime(rt),
129 fTime(-1)
130{
131 fHist.SetDirectory(0);
132
133 fAfterpulses.SetOwner();
134
135 fAfterpulseProb[0] = 0;
136 fAfterpulseProb[1] = 0;
137
138 fAfterpulseTau[0] = 15;
139 fAfterpulseTau[1] = 85;
140}
141
142// --------------------------------------------------------------------------
143//
144// This is the time a chips needs after an external signal to relax to
145// a "virgin" state, i.e. without no influence of the external pulse
146// above the given threshold.
147//
148// It takes into account the dead time of the cell, the relaxation time
149// and the two afterpulse distributions. However, in most cases the
150// afterpulse distribution will dominate (except they are switched off by
151// a zero probability).
152//
153// FIXME: Maybe the calculation of the relaxation time could be optimized?
154//
155Float_t APD::GetRelaxationTime(Float_t threshold) const
156{
157 // Calculate time until the probability of one of these
158 // events falls below the threshold probability.
159 const Double_t rt0 = - TMath::Log(threshold)*fRecoveryTime;
160 const Double_t rt1 = fAfterpulseProb[0]>0 ? -TMath::Log(threshold/fAfterpulseProb[0])*fAfterpulseTau[0] : 0;
161 const Double_t rt2 = fAfterpulseProb[1]>0 ? -TMath::Log(threshold/fAfterpulseProb[1])*fAfterpulseTau[1] : 0;
162
163 // Probability not between t and inf, but between t and t+dt
164 // -tau * log ( p / ( 1 - exp(- dt/tau) ) ) = t
165
166 return fDeadTime + TMath::Max(rt0, TMath::Max(rt1, rt2));
167}
168
169// --------------------------------------------------------------------------
170//
171// This is the recursive implementation of a hit. If a photon hits a cell
172// at x and y (must be a valid cell!) at time t, at first we check if the
173// cell is still dead. If it is not dead we calculate the signal height
174// from the recovery time. Now we check with the crosstalk probability
175// whether another photon is created. If another photon is created we
176// calculate randomly which of the four neighbor cells are hit.
177// If the cell is outside the APD the photon is ignored. As many
178// new photons are created until our random number is below the crosstak-
179// probability.
180//
181// For each photon the possible afterpulses of two distributions are
182// created and added to the list of afterpulses. This is done by calling
183// GenerateAfterpulse for the two afterpulse-distributions.
184//
185// The total height of the signal (in units of photons) is returned.
186// Note, that this can be a fractional number.
187//
188// This function looks a bit fancy accessing the histogram and works around
189// a few histogram functions. This is a speed optimization which works
190// around a lot of sanity checks which are obsolete in our case.
191//
192// The default time is 0.
193//
194Float_t APD::HitCellImp(Int_t x, Int_t y, Float_t t)
195{
196 // if (x<1 || x>fHist.GetNbinsX() ||
197 // y<1 || y>fHist.GetNbinsY())
198 // return 0;
199#ifdef DEBUG
200 cout << "Hit: " << t << endl;
201#endif
202
203 // Number of the x/y cell in the one dimensional array
204 // const Int_t cell = fHist.GetBin(x, y);
205 const Int_t cell = x + (fHist.GetNbinsX()+2)*y;
206
207 // Getting a reference to the float is the fastes way to
208 // access the bin-contents in fArray
209 Float_t &cont = fHist.GetArray()[cell];
210
211 // Calculate the time since the last breakdown
212 // const Double_t dt = t-fHist.GetBinContent(x, y)-fDeadTime; //
213 const Float_t dt = t-cont-fDeadTime;
214
215 // Photons within the dead time are just ignored
216 if (/*hx.GetBinContent(x,y)>0 &&*/ dt<=0)
217 {
218#ifdef DEBUG
219 cout << "Dead: " << t << " " << cont << " " << dt << endl;
220#endif
221 return 0;
222 }
223 // The signal height (in units of one photon) produced after dead time
224 // depends on the recovery of the cell - described by an exponential.
225 const Float_t weight = fRecoveryTime<=0 ? 1. : 1-TMath::Exp(-dt/fRecoveryTime);
226
227 // Now we know the charge in the cell and we can generate
228 // the afterpulses with both time-constants
229 GenerateAfterpulse(cell, 0, weight, t);
230 GenerateAfterpulse(cell, 1, weight, t);
231
232 // The probability that the cell emits a photon causing crosstalk
233 // scales as the signal height.
234 const Float_t prob = weight*fCrosstalkProb;
235
236 // Set the contents to the time of the last breakdown (now)
237 cont = t; // fHist.SetBinContent(x, y, t)
238
239 // Counter for the numbers of produced photons
240 Float_t n = weight;
241
242 // Get random number of emitted and possible converted crosstalk photons
243 const UInt_t rndm = gRandom->Poisson(prob);
244
245 for (int i=0; i<rndm; i++)
246 {
247 // Get a random neighbor which is hit.
248 switch (gRandom->Integer(4))
249 {
250 case 0: if (x<fHist.GetNbinsX()) n += HitCellImp(x+1, y, t); break;
251 case 1: if (x>1) n += HitCellImp(x-1, y, t); break;
252 case 2: if (y<fHist.GetNbinsY()) n += HitCellImp(x, y+1, t); break;
253 case 3: if (y>1) n += HitCellImp(x, y-1, t); break;
254 }
255 }
256
257 return n;
258}
259
260// --------------------------------------------------------------------------
261//
262// Check if x and y is a valid cell. If not return 0, otherwise
263// HitCelImp(x, y, t). HitCellImp generates Crosstalk and Afterpulses.
264//
265// The default time is 0.
266//
267Float_t APD::HitCell(Int_t x, Int_t y, Float_t t)
268{
269 if (x<1 || x>fHist.GetNbinsX() ||
270 y<1 || y>fHist.GetNbinsY())
271 return 0;
272
273 return HitCellImp(x, y, t);
274}
275
276// --------------------------------------------------------------------------
277//
278// Determine randomly (uniformly) a cell which was hit. Return
279// HitCellImp for this cell and the given time. HitCellImp
280// generates Crosstalk and Afterpulses
281//
282// The default time is 0.
283//
284// If you want t w.r.t. fTime use HitRandomCellRelative istead.
285//
286Float_t APD::HitRandomCell(Float_t t)
287{
288 const UInt_t nx = fHist.GetNbinsX();
289 const UInt_t ny = fHist.GetNbinsY();
290
291 const UInt_t idx = gRandom->Integer(nx*ny);
292
293 const UInt_t x = idx%nx;
294 const UInt_t y = idx/nx;
295
296 return HitCellImp(x+1, y+1, t);
297}
298
299// --------------------------------------------------------------------------
300//
301// Sets all cells with a contents which is well before the time t such that
302// the chip is "virgin". Therefore all cells are set to a time which
303// is twice the deadtime before the given time and 1000 times the recovery
304// time.
305//
306// The afterpulse list is deleted.
307//
308// If deadtime and recovery time are 0 then t-1 is set.
309//
310// Sets fTime to t
311//
312// The default time is 0.
313//
314void APD::FillEmpty(Float_t t)
315{
316 const Int_t n = (fHist.GetNbinsX()+2)*(fHist.GetNbinsY()+2);
317
318 const Double_t tm = fDeadTime<=0 && fRecoveryTime<=0 ? t-1 : t-2*fDeadTime-1000*fRecoveryTime;
319
320 for (int i=0; i<n; i++)
321 fHist.GetArray()[i] = tm;
322
323 fHist.SetEntries(1);
324
325 fAfterpulses.Delete();
326
327 fTime = t;
328}
329
330// --------------------------------------------------------------------------
331//
332// First call FillEmpty for the given time t. Then fill each cell by
333// by calling HitCellImp with time t-gRandom->Exp(n/rate) with n being
334// the total number of cells. This the time at which the cell was last hit.
335//
336// Sets fTime to t
337//
338// If the argument t is omitted it defaults to 0.
339//
340// Since after calling this function the chip should reflect the
341// status at the new time fTime=t, all afterpulses are processed
342// until this time. However, the produced random pulses might have produced
343// new new afterpulses.
344//
345// All afterpulses before the new timestamp are deleted.
346//
347// WARNING: Note that this might not correctly reproduce afterpulses
348// produced by earlier pulese.
349//
350void APD::FillRandom(Float_t rate, Float_t t)
351{
352 FillEmpty(t);
353
354 const Int_t nx = fHist.GetNbinsX();
355 const Int_t ny = fHist.GetNbinsY();
356
357 const Double_t f = (nx*ny)/rate;
358
359 // FIXME: Dead time is not taken into account,
360 // possible earlier afterpulses are not produced.
361
362 for (int x=1; x<=nx; x++)
363 for (int y=1; y<=ny; y++)
364 HitCellImp(x, y, t-MMath::RndmExp(f));
365
366 // Deleting of the afterpulses before fHist.GetMinimum() won't
367 // speed things because we have to loop over them once in any case
368
369 ProcessAfterpulses(fHist.GetMinimum(), t);
370 DeleteAfterpulses(t);
371
372 fTime = t;
373}
374
375// --------------------------------------------------------------------------
376//
377// Shift all times including fTime to dt (ie. add -dt to all times)
378// This allows to set a user-defined T0 or shift T0 to fTime=0.
379//
380// However, T0<0 is not allowed (dt cannot be greater than fTime)
381//
382void APD::ShiftTime(Double_t dt)
383{
384 if (dt>fTime)
385 {
386 gLog << warn << "APD::ShiftTime: WARNING - requested time shift results in fTime<0... ignored." << endl;
387 return;
388 }
389
390 // If reset was requested shift all times by end backwards
391 // so that fTime is now 0
392 const Int_t n = (fHist.GetNbinsX()+2)*(fHist.GetNbinsY()+2);
393 for (int i=0; i<n; i++)
394 fHist.GetArray()[i] -= dt;
395
396 fTime -= dt;
397}
398
399// --------------------------------------------------------------------------
400//
401// Evolve the chip from fTime to fTime+dt if it with a given frequency
402// freq. Returns the total signal "recorded".
403//
404// fTime is set to the fTime+dt.
405//
406// If you want to evolve over a default relaxation time (relax the chip
407// from a signal) use Relax instead.
408//
409// Since after calling this function the chip should reflect the
410// status at the new time fTime=fTime+dt, all afterpulses are processed
411// until this time. However, evolving the chip until this time might
412// have produced new afterpulses.
413//
414// All afterpulses before the new timestamp are deleted.
415//
416Float_t APD::Evolve(Double_t freq, Double_t dt)
417{
418 const Double_t end = fTime+dt;
419
420 Float_t hits = 0;
421
422 if (freq>0)
423 {
424 const Double_t avglen = 1./freq;
425
426 Double_t time = fTime;
427 while (1)
428 {
429 const Double_t deltat = MMath::RndmExp(avglen);
430 if (time+deltat>end)
431 break;
432
433 time += deltat;
434
435 hits += HitRandomCell(time);
436 }
437 }
438
439 // Deleting of the afterpulses before fTime won't speed things
440 // because we have to loop over them once in any case
441
442 ProcessAfterpulses(fTime, dt);
443 DeleteAfterpulses(end);
444
445 fTime = end;
446
447 return hits;
448}
449
450// --------------------------------------------------------------------------
451//
452// Retunrs the number of cells which have a time t<=fDeadTime, i.e. which are
453// dead.
454// The default time is 0.
455//
456// Note that if you want to get a correct measure of teh number of dead cells
457// at the time t, this function will only produce a valid count if the
458// afterpulses have been processed up to this time.
459//
460Int_t APD::CountDeadCells(Float_t t) const
461{
462 const Int_t nx = fHist.GetNbinsX();
463 const Int_t ny = fHist.GetNbinsY();
464
465 Int_t n=0;
466 for (int x=1; x<=nx; x++)
467 for (int y=1; y<=ny; y++)
468 if ((t-fHist.GetBinContent(x, y))<=fDeadTime)
469 n++;
470
471 return n;
472}
473
474// --------------------------------------------------------------------------
475//
476// Returs the number of cells which have a time t<=fDeadTime+fRecoveryTime.
477// The default time is 0.
478//
479// Note that if you want to get a correct measure of teh number of dead cells
480// at the time t, this function will only produce a valid count if the
481// afterpulses have been processed up to this time.
482//
483Int_t APD::CountRecoveringCells(Float_t t) const
484{
485 const Int_t nx = fHist.GetNbinsX();
486 const Int_t ny = fHist.GetNbinsY();
487
488 Int_t n=0;
489 for (int x=1; x<=nx; x++)
490 for (int y=1; y<=ny; y++)
491 {
492 Float_t dt = t-fHist.GetBinContent(x, y);
493 if (dt>fDeadTime && dt<=fDeadTime+fRecoveryTime)
494 n++;
495 }
496 return n;
497}
498
499// --------------------------------------------------------------------------
500//
501// Generate an afterpulse originating from the given cell and a pulse with
502// charge. The afterpulse distribution to use is specified by
503// the index. The "current" time to which the afterpulse delay refers must
504// be given by t.
505//
506// A generated Afterpulse is added to the list of afterpulses
507//
508void APD::GenerateAfterpulse(UInt_t cell, Int_t idx, Double_t charge, Double_t t)
509{
510 // The cell had a single avalanche with signal height weight.
511 // This cell now can produce an afterpulse photon/avalanche.
512 const Double_t p = gRandom->Uniform();
513
514 // It's probability scales with the charge of the pulse
515 if (p>charge*fAfterpulseProb[idx])
516 return;
517
518 // Afterpulses come with a well defined time-constant
519 // after the normal pulse
520 const Double_t dt = MMath::RndmExp(fAfterpulseTau[idx]);
521
522 fAfterpulses.Add(new Afterpulse(cell, t+dt));
523
524#ifdef DEBUG
525 cout << "Add : " << t << " + " << dt << " = " << t+dt << endl;
526#endif
527}
528
529// --------------------------------------------------------------------------
530//
531// Process afterpulses between time and time+dt. All afterpulses in the list
532// before t=time are ignored. All afterpulses between t=time and
533// t=time+dt are processed through HitCellImp. Afterpulses after and
534// equal t=time+dt are skipped.
535//
536// Since the afterpulse list is a sorted list newly generated afterpulses
537// are always inserted into the list behind the current entry. Consequently,
538// afterpulses generated by afterpulses will also be processed correctly.
539//
540// Afterpulses with zero amplitude are deleted from the list. All other after
541// pulses remain in the list for later evaluation.
542//
543void APD::ProcessAfterpulses(Float_t time, Float_t dt)
544{
545#ifdef DEBUG
546 cout << "Process afterpulses from " << time << " to " << time+dt << endl;
547#endif
548
549 const Float_t end = time+dt;
550
551 TObjLink *lnk = fAfterpulses.FirstLink();
552 while (lnk)
553 {
554 Afterpulse &ap = *static_cast<Afterpulse*>(lnk->GetObject());
555
556 // Skip afterpulses which have been processed already
557 // or which we do not have to process anymore
558 if (ap.GetTime()<time || ap.GetAmplitude()>0)
559 {
560 lnk = lnk->Next();
561 continue;
562 }
563
564 // No afterpulses left in correct time window
565 if (ap.GetTime()>=end)
566 break;
567
568 // Process afterpulse through HitCellImp
569 const Float_t ampl = ap.Process(*this);
570
571 // Step to the next entry already, the current one might get deleted
572 lnk = lnk->Next();
573
574 if (ampl!=0)
575 continue;
576
577#ifdef DEBUG
578 cout << "Del : " << ap.GetTime() << " (" << ampl << ")" << endl;
579#endif
580
581 // The afterpulse "took place" within the dead time of the
582 // pixel ==> No afterpulse, no crosstalk.
583 delete fAfterpulses.Remove(&ap);
584 }
585}
586
587// --------------------------------------------------------------------------
588//
589// Delete all afterpulses before and equal to t=time from the list
590//
591void APD::DeleteAfterpulses(Float_t time)
592{
593 TIter Next(&fAfterpulses);
594
595 Afterpulse *ap = 0;
596
597 while ((ap = static_cast<Afterpulse*>(Next())))
598 {
599 if (ap->GetTime()>=time)
600 break;
601
602 delete fAfterpulses.Remove(ap);
603 }
604}
Note: See TracBrowser for help on using the repository browser.