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

Last change on this file since 10093 was 10093, checked in by tbretz, 12 years ago
Added afterpulse treatment to MAvalanchePhotoDiode and MSimAPD
File size: 19.9 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[0]>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 /*
243 // Check if a photon in a neighboring cell is produced (crosstalk)
244 while (gRandom->Rndm()<fCrosstalkProb)
245 {
246 // Get a random neighbor which is hit.
247 switch (gRandom->Integer(4))
248 {
249 case 0: x++; if (x>fHist.GetNbinsX()) continue; break;
250 case 1: x--; if (x<1) continue; break;
251 case 2: y++; if (y>fHist.GetNbinsY()) continue; break;
252 case 3: y--; if (y<1) continue; break;
253 }
254
255 n += HitCellImp(x, y, t);
256 }
257 */
258
259
260 //for (int i=0; i<1; i++)
261 while (1)
262 {
263 const Double_t rndm = gRandom->Rndm();
264 if (rndm>=prob/*fCrosstalkProb*/)
265 break;
266
267 // We can re-use the random number because it is uniformely
268 // distributed. This saves cpu power
269 const Int_t dir = TMath::FloorNint(4*rndm/prob/*fCrosstalkProb*/);
270
271 // Get a random neighbor which is hit.
272 switch (dir)
273 {
274 case 0: if (x<fHist.GetNbinsX()) n += HitCellImp(x+1, y, t); break;
275 case 1: if (x>1) n += HitCellImp(x-1, y, t); break;
276 case 2: if (y<fHist.GetNbinsY()) n += HitCellImp(x, y+1, t); break;
277 case 3: if (y>1) n += HitCellImp(x, y-1, t); break;
278 }
279
280 // In the unlikely case the calculated direction is out-of-range,
281 // i.e. <0 or >3, we would just try to fill the same cell again which
282 }
283
284 return n;
285}
286
287// --------------------------------------------------------------------------
288//
289// Check if x and y is a valid cell. If not return 0, otherwise
290// HitCelImp(x, y, t). HitCellImp generates Crosstalk and Afterpulses.
291//
292// The default time is 0.
293//
294Float_t APD::HitCell(Int_t x, Int_t y, Float_t t)
295{
296 if (x<1 || x>fHist.GetNbinsX() ||
297 y<1 || y>fHist.GetNbinsY())
298 return 0;
299
300 return HitCellImp(x, y, t);
301}
302
303// --------------------------------------------------------------------------
304//
305// Determine randomly (uniformly) a cell which was hit. Return
306// HitCellImp for this cell and the given time. HitCellImp
307// generates Crosstalk and Afterpulses
308//
309// The default time is 0.
310//
311// If you want t w.r.t. fTime use HitRandomCellRelative istead.
312//
313Float_t APD::HitRandomCell(Float_t t)
314{
315 const UInt_t nx = fHist.GetNbinsX();
316 const UInt_t ny = fHist.GetNbinsY();
317
318 const UInt_t idx = gRandom->Integer(nx*ny);
319
320 const UInt_t x = idx%nx;
321 const UInt_t y = idx/nx;
322
323 return HitCellImp(x+1, y+1, t);
324}
325
326// --------------------------------------------------------------------------
327//
328// Sets all cells with a contents which is well before the time t such that
329// the chip is "virgin". Therefore all cells are set to a time which
330// is twice the deadtime before the given time and 1000 times the recovery
331// time.
332//
333// The afterpulse list is deleted.
334//
335// If deadtime and recovery time are 0 then t-1 is set.
336//
337// Sets fTime to t
338//
339// The default time is 0.
340//
341void APD::FillEmpty(Float_t t)
342{
343 const Int_t n = (fHist.GetNbinsX()+2)*(fHist.GetNbinsY()+2);
344
345 const Double_t tm = fDeadTime<=0 && fRecoveryTime<=0 ? t-1 : t-2*fDeadTime-1000*fRecoveryTime;
346
347 for (int i=0; i<n; i++)
348 fHist.GetArray()[i] = tm;
349
350 fHist.SetEntries(1);
351
352 fAfterpulses.Delete();
353
354 fTime = t;
355}
356
357// --------------------------------------------------------------------------
358//
359// First call FillEmpty for the given time t. Then fill each cell by
360// by calling HitCellImp with time t-gRandom->Exp(n/rate) with n being
361// the total number of cells. This the time at which the cell was last hit.
362//
363// Sets fTime to t
364//
365// If the argument t is omitted it defaults to 0.
366//
367// Since after calling this function the chip should reflect the
368// status at the new time fTime=t, all afterpulses are processed
369// until this time. However, the produced random pulses might have produced
370// new new afterpulses.
371//
372// All afterpulses before the new timestamp are deleted.
373//
374// WARNING: Note that this might not correctly reproduce afterpulses
375// produced by earlier pulese.
376//
377void APD::FillRandom(Float_t rate, Float_t t)
378{
379 FillEmpty(t);
380
381 const Int_t nx = fHist.GetNbinsX();
382 const Int_t ny = fHist.GetNbinsY();
383
384 const Double_t f = (nx*ny)/rate;
385
386 // FIXME: Dead time is not taken into account,
387 // possible earlier afterpulses are not produced.
388
389 for (int x=1; x<=nx; x++)
390 for (int y=1; y<=ny; y++)
391 HitCellImp(x, y, t-MMath::RndmExp(f));
392
393 // Deleting of the afterpulses before fHist.GetMinimum() won't
394 // speed things because we have to loop over them once in any case
395
396 ProcessAfterpulses(fHist.GetMinimum(), t);
397 DeleteAfterpulses(t);
398
399 fTime = t;
400}
401
402// --------------------------------------------------------------------------
403//
404// Shift all times including fTime to dt (ie. add -dt to all times)
405// This allows to set a user-defined T0 or shift T0 to fTime=0.
406//
407// However, T0<0 is not allowed (dt cannot be greater than fTime)
408//
409void APD::ShiftTime(Double_t dt)
410{
411 if (dt>fTime)
412 {
413 gLog << warn << "APD::ShiftTime: WARNING - requested time shift results in fTime<0... ignored." << endl;
414 return;
415 }
416
417 // If reset was requested shift all times by end backwards
418 // so that fTime is now 0
419 const Int_t n = (fHist.GetNbinsX()+2)*(fHist.GetNbinsY()+2);
420 for (int i=0; i<n; i++)
421 fHist.GetArray()[i] -= dt;
422
423 fTime -= dt;
424}
425
426// --------------------------------------------------------------------------
427//
428// Evolve the chip from fTime to fTime+dt if it with a given frequency
429// freq. Returns the total signal "recorded".
430//
431// fTime is set to the fTime+dt.
432//
433// If you want to evolve over a default relaxation time (relax the chip
434// from a signal) use Relax instead.
435//
436// Since after calling this function the chip should reflect the
437// status at the new time fTime=fTime+dt, all afterpulses are processed
438// until this time. However, evolving the chip until this time might
439// have produced new afterpulses.
440//
441// All afterpulses before the new timestamp are deleted.
442//
443Float_t APD::Evolve(Double_t freq, Double_t dt)
444{
445 const Double_t end = fTime+dt;
446
447 Float_t hits = 0;
448
449 if (freq>0)
450 {
451 const Double_t avglen = 1./freq;
452
453 Double_t time = fTime;
454 while (1)
455 {
456 const Double_t deltat = MMath::RndmExp(avglen);
457 if (time+deltat>end)
458 break;
459
460 time += deltat;
461
462 hits += HitRandomCell(time);
463 }
464 }
465
466 // Deleting of the afterpulses before fTime won't speed things
467 // because we have to loop over them once in any case
468
469 ProcessAfterpulses(fTime, dt);
470 DeleteAfterpulses(end);
471
472 fTime = end;
473
474 return hits;
475}
476
477// --------------------------------------------------------------------------
478//
479// Retunrs the number of cells which have a time t<=fDeadTime, i.e. which are
480// dead.
481// The default time is 0.
482//
483// Note that if you want to get a correct measure of teh number of dead cells
484// at the time t, this function will only produce a valid count if the
485// afterpulses have been processed up to this time.
486//
487Int_t APD::CountDeadCells(Float_t t) const
488{
489 const Int_t nx = fHist.GetNbinsX();
490 const Int_t ny = fHist.GetNbinsY();
491
492 Int_t n=0;
493 for (int x=1; x<=nx; x++)
494 for (int y=1; y<=ny; y++)
495 if ((t-fHist.GetBinContent(x, y))<=fDeadTime)
496 n++;
497
498 return n;
499}
500
501// --------------------------------------------------------------------------
502//
503// Returs the number of cells which have a time t<=fDeadTime+fRecoveryTime.
504// The default time is 0.
505//
506// Note that if you want to get a correct measure of teh number of dead cells
507// at the time t, this function will only produce a valid count if the
508// afterpulses have been processed up to this time.
509//
510Int_t APD::CountRecoveringCells(Float_t t) const
511{
512 const Int_t nx = fHist.GetNbinsX();
513 const Int_t ny = fHist.GetNbinsY();
514
515 Int_t n=0;
516 for (int x=1; x<=nx; x++)
517 for (int y=1; y<=ny; y++)
518 {
519 Float_t dt = t-fHist.GetBinContent(x, y);
520 if (dt>fDeadTime && dt<=fDeadTime+fRecoveryTime)
521 n++;
522 }
523 return n;
524}
525
526// --------------------------------------------------------------------------
527//
528// Generate an afterpulse originating from the given cell and a pulse with
529// charge. The afterpulse distribution to use is specified by
530// the index. The "current" time to which the afterpulse delay refers must
531// be given by t.
532//
533// A generated Afterpulse is added to the list of afterpulses
534//
535void APD::GenerateAfterpulse(UInt_t cell, Int_t idx, Double_t charge, Double_t t)
536{
537 // The cell had a single avalanche with signal height weight.
538 // This cell now can produce an afterpulse photon/avalanche.
539 const Double_t p = gRandom->Uniform();
540
541 // It's probability scales with the charge of the pulse
542 if (p>charge*fAfterpulseProb[idx])
543 return;
544
545 // Afterpulses come with a well defined time-constant
546 // after the normal pulse
547 const Double_t dt = MMath::RndmExp(fAfterpulseTau[idx]);
548
549 fAfterpulses.Add(new Afterpulse(cell, t+dt));
550
551#ifdef DEBUG
552 cout << "Add : " << t << " + " << dt << " = " << t+dt << endl;
553#endif
554}
555
556// --------------------------------------------------------------------------
557//
558// Process afterpulses between time and time+dt. All afterpulses in the list
559// before t=time are ignored. All afterpulses between t=time and
560// t=time+dt are processed through HitCellImp. Afterpulses after and
561// equal t=time+dt are skipped.
562//
563// Since the afterpulse list is a sorted list newly generated afterpulses
564// are always inserted into the list behind the current entry. Consequently,
565// afterpulses generated by afterpulses will also be processed correctly.
566//
567// Afterpulses with zero amplitude are deleted from the list. All other after
568// pulses remain in the list for later evaluation.
569//
570void APD::ProcessAfterpulses(Float_t time, Float_t dt)
571{
572#ifdef DEBUG
573 cout << "Process afterpulses from " << time << " to " << time+dt << endl;
574#endif
575
576 const Float_t end = time+dt;
577
578 TObjLink *lnk = fAfterpulses.FirstLink();
579 while (lnk)
580 {
581 Afterpulse &ap = *static_cast<Afterpulse*>(lnk->GetObject());
582
583 // Skip afterpulses which have been processed already
584 // or which we do not have to process anymore
585 if (ap.GetTime()<time || ap.GetAmplitude()>0)
586 {
587 lnk = lnk->Next();
588 continue;
589 }
590
591 // No afterpulses left in correct time window
592 if (ap.GetTime()>=end)
593 break;
594
595 // Process afterpulse through HitCellImp
596 const Float_t ampl = ap.Process(*this);
597
598 // Step to the next entry already, the current one might get deleted
599 lnk = lnk->Next();
600
601 if (ampl!=0)
602 continue;
603
604#ifdef DEBUG
605 cout << "Del : " << ap.GetTime() << " (" << ampl << ")" << endl;
606#endif
607
608 // The afterpulse "took place" within the dead time of the
609 // pixel ==> No afterpulse, no crosstalk.
610 delete fAfterpulses.Remove(&ap);
611 }
612}
613
614// --------------------------------------------------------------------------
615//
616// Delete all afterpulses before and equal to t=time from the list
617//
618void APD::DeleteAfterpulses(Float_t time)
619{
620 TIter Next(&fAfterpulses);
621
622 Afterpulse *ap = 0;
623
624 while ((ap = static_cast<Afterpulse*>(Next())))
625 {
626 if (ap->GetTime()>=time)
627 break;
628
629 delete fAfterpulses.Remove(ap);
630 }
631}
Note: See TracBrowser for help on using the repository browser.