source: trunk/MagicSoft/Mars/melectronics/MAvalanchePhotoDiode.cc@ 9565

Last change on this file since 9565 was 9565, checked in by tbretz, 13 years ago
*** empty log message ***
File size: 11.1 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#include "MAvalanchePhotoDiode.h"
43
44#include <TRandom.h>
45
46#include "MMath.h"
47
48ClassImp(APD);
49
50using namespace std;
51
52/*
53class MyProfile : public TProfile2D
54{
55public:
56 void AddBinEntry(Int_t cell) { fBinEntries.fArray[cell]++; }
57};
58*/
59
60// --------------------------------------------------------------------------
61//
62// Default Constructor.
63//
64// n is the number od cells in x or y. The APD is assumed to
65// be square.
66// prob is the crosstalk probability, i.e., the probability that a
67// photon which produced an avalanche will create another
68// photon in a neighboring cell
69// dt is the deadtime, i.e., the time in which the APD cell will show
70// no response to a photon after a hit
71// rt is the recovering tims, i.e. the exponential (e^(-dt/rt))
72// with which the cell is recovering after being dead
73//
74// prob, dt and ar can be set to 0 to switch the effect off.
75// 0 is also the dfeault for all three.
76//
77APD::APD(Int_t n, Float_t prob, Float_t dt, Float_t rt)
78 : fHist("APD", "", n, 0.5, n+0.5, n, 0.5, n+0.5),
79 fCrosstalkProb(prob), fDeadTime(dt), fRecoveryTime(rt), fTime(-1)
80{
81 fHist.SetDirectory(0);
82}
83
84// --------------------------------------------------------------------------
85//
86// This is the time the cell needs after a hit until less than threshold
87// (0.001 == one permille) of the signal is lost.
88//
89// If all cells of a G-APD have fired simultaneously and they are fired
90// once more after the relaxation time (with the defaultthreshold of 0.001)
91// the chip will show a signal which is one permille too small.
92//
93Float_t APD::GetRelaxationTime(Float_t threshold) const
94{
95 return fDeadTime-TMath::Log(threshold)*fRecoveryTime;
96}
97
98// --------------------------------------------------------------------------
99//
100// This is the recursive implementation of a hit. If a photon hits a cell
101// at x and y (must be a valid cell!) at time t, at first we check if the
102// cell is still dead. If it is not dead we calculate the signal height
103// from the recovery time. Now we check with the crosstalk probability
104// whether another photon is created. If another photon is created we
105// calculate randomly which of the four neighbor cells are hit.
106// If the cell is outside the APD the photon is ignored. As many
107// new photons are created until our random number is below the crosstak-
108// probability.
109//
110// The total height of the signal (in units of photons) is returned.
111// Note, that this can be a fractional number.
112//
113// This function looks a bit fancy accessing the histogram and works around
114// a few histogram functions. This is a speed optimization which works
115// around a lot of sanity checks which are obsolete in our case.
116//
117// The default time is 0.
118//
119Float_t APD::HitCellImp(Int_t x, Int_t y, Float_t t)
120{
121 // if (x<1 || x>fHist.GetNbinsX() ||
122 // y<1 || y>fHist.GetNbinsY())
123 // return 0;
124
125 // Number of the x/y cell in the one dimensional array
126 // const Int_t cell = fHist.GetBin(x, y);
127 const Int_t cell = x + (fHist.GetNbinsX()+2)*y;
128
129 // Getting a reference to the float is the fastes way to
130 // access the bin-contents in fArray
131 Float_t &cont = fHist.GetArray()[cell];
132
133 // Calculate the time since the last breakdown
134 // const Double_t dt = t-fHist.GetBinContent(x, y)-fDeadTime; //
135 const Float_t dt = t-cont-fDeadTime;
136
137 // Photons within the dead time are just ignored
138 if (/*hx.GetBinContent(x,y)>0 &&*/ dt<=0)
139 return 0;
140
141 // The signal height (in units of one photon) produced after dead time
142 // depends on the recovery of the cell - described by an exponential.
143 const Float_t weight = fRecoveryTime<=0 ? 1. : 1-TMath::Exp(-dt/fRecoveryTime);
144
145 // The probability that the cell emits a photon causing crosstalk
146 // scales as the signal height.
147 const Float_t prob = weight*fCrosstalkProb;
148
149 // Set the contents to the time of the last breakdown (now)
150 cont = t; // fHist.SetBinContent(x, y, t)
151
152 // Counter for the numbers of produced photons
153 Float_t n = weight;
154
155 /*
156 // Check if a photon in a neighboring cell is produced (crosstalk)
157 while (gRandom->Rndm()<fCrosstalkProb)
158 {
159 // Get a random neighbor which is hit.
160 switch (gRandom->Integer(4))
161 {
162 case 0: x++; if (x>fHist.GetNbinsX()) continue; break;
163 case 1: x--; if (x<1) continue; break;
164 case 2: y++; if (y>fHist.GetNbinsY()) continue; break;
165 case 3: y--; if (y<1) continue; break;
166 }
167
168 n += HitCellImp(x, y, t);
169 }
170 */
171
172
173 //for (int i=0; i<1; i++)
174 while (1)
175 {
176 const Double_t rndm = gRandom->Rndm();
177 if (rndm>=prob/*fCrosstalkProb*/)
178 break;
179
180 // We can re-use the random number because it is uniformely
181 // distributed. This saves cpu power
182 const Int_t dir = TMath::FloorNint(4*rndm/prob/*fCrosstalkProb*/);
183
184 // Get a random neighbor which is hit.
185 switch (dir)
186 {
187 case 0: if (x<fHist.GetNbinsX()) n += HitCellImp(x+1, y, t); break;
188 case 1: if (x>1) n += HitCellImp(x-1, y, t); break;
189 case 2: if (y<fHist.GetNbinsY()) n += HitCellImp(x, y+1, t); break;
190 case 3: if (y>1) n += HitCellImp(x, y-1, t); break;
191 }
192
193 // In the unlikely case the calculated direction is out-of-range,
194 // i.e. <0 or >3, we would just try to fill the same cell again which
195 }
196
197 return n;
198}
199
200// --------------------------------------------------------------------------
201//
202// Check if x and y is a valid cell. If not return 0, otherwise
203// HitCelImp(x, y, t)
204//
205// The default time is 0.
206//
207Float_t APD::HitCell(Int_t x, Int_t y, Float_t t)
208{
209 if (x<1 || x>fHist.GetNbinsX() ||
210 y<1 || y>fHist.GetNbinsY())
211 return 0;
212
213 return HitCellImp(x, y, t);
214}
215
216// --------------------------------------------------------------------------
217//
218// Determine randomly (uniformly) a cell which was hit. Return
219// HitCellImp for this cell and the given time.
220//
221// The default time is 0.
222//
223// If you want t w.r.t. fTime use HitRandomCellRelative istead.
224//
225Float_t APD::HitRandomCell(Float_t t)
226{
227 const UInt_t nx = fHist.GetNbinsX();
228 const UInt_t ny = fHist.GetNbinsY();
229
230 const UInt_t idx = gRandom->Integer(nx*ny);
231
232 const UInt_t x = idx%nx;
233 const UInt_t y = idx/nx;
234
235 return HitCellImp(x+1, y+1, t);
236}
237
238// --------------------------------------------------------------------------
239//
240// Sets all cells with a contents whihc is well before the time t such that
241// the chip is "virgin". Therefore all cells are set to a time which
242// is twice the deadtime before the given time and 1000 times the recovery
243// time.
244//
245// If deadtime and recovery time are 0 then t-1 is set.
246//
247// Sets fTime to t
248//
249// The default time is 0.
250//
251void APD::FillEmpty(Float_t t)
252{
253 const Int_t n = (fHist.GetNbinsX()+2)*(fHist.GetNbinsY()+2);
254
255 const Double_t tm = fDeadTime<=0 && fRecoveryTime<=0 ? t-1 : t-2*fDeadTime-1000*fRecoveryTime;
256
257 for (int i=0; i<n; i++)
258 fHist.GetArray()[i] = tm;
259
260 fHist.SetEntries(1);
261
262 fTime = t;
263}
264
265// --------------------------------------------------------------------------
266//
267// First call FillEmpty for the given time t. Then fill each cell by
268// by calling HitCellImp with time t-gRandom->Exp(n/rate) with n being
269// the total number of cells.
270//
271// Sets fTime to t
272//
273// The default time is 0.
274//
275void APD::FillRandom(Float_t rate, Float_t t)
276{
277 FillEmpty(t);
278
279 const Int_t nx = fHist.GetNbinsX();
280 const Int_t ny = fHist.GetNbinsY();
281
282 const Double_t f = (nx*ny)/rate;
283
284 // FIXME: This is not perfect, is it? What about the dead time?
285
286 for (int x=1; x<=nx; x++)
287 for (int y=1; y<=ny; y++)
288 HitCellImp(x, y, t-MMath::RndmExp(f));
289
290 fTime = t;
291}
292
293// --------------------------------------------------------------------------
294//
295// Evolve the chip from fTime to fTime+dt if it with a given frequency
296// freq. Returns the total signal "recorded".
297//
298// fTime is set to the fTime+dt
299//
300// If you want to evolve over a default relaxation time (relax the chip
301// from a signal) use Relax instead.
302//
303Float_t APD::Evolve(Double_t freq, Double_t dt)
304{
305 const Double_t avglen = 1./freq;
306
307 const Double_t end = fTime+dt;
308
309 Float_t hits = 0;
310
311 Double_t time = fTime;
312 while (1)
313 {
314 time += MMath::RndmExp(avglen);
315 if (time>end)
316 break;
317
318 hits += HitRandomCell(time);
319 }
320
321 fTime = end;
322
323 return hits;
324}
325
326// --------------------------------------------------------------------------
327//
328// Retunrs the number of cells which have a time t<=fDeadTime, i.e. which are
329// dead.
330// The default time is 0.
331//
332Int_t APD::CountDeadCells(Float_t t) const
333{
334 const Int_t nx = fHist.GetNbinsX();
335 const Int_t ny = fHist.GetNbinsY();
336
337 Int_t n=0;
338 for (int x=1; x<=nx; x++)
339 for (int y=1; y<=ny; y++)
340 if ((t-fHist.GetBinContent(x, y))<=fDeadTime)
341 n++;
342
343 return n;
344}
345
346// --------------------------------------------------------------------------
347//
348// Returs the number of cells which have a time t<=fDeadTime+fRecoveryTime.
349// The default time is 0.
350//
351Int_t APD::CountRecoveringCells(Float_t t) const
352{
353 const Int_t nx = fHist.GetNbinsX();
354 const Int_t ny = fHist.GetNbinsY();
355
356 Int_t n=0;
357 for (int x=1; x<=nx; x++)
358 for (int y=1; y<=ny; y++)
359 {
360 Float_t dt = t-fHist.GetBinContent(x, y);
361 if (dt>fDeadTime && dt<=fDeadTime+fRecoveryTime)
362 n++;
363 }
364 return n;
365}
Note: See TracBrowser for help on using the repository browser.