Index: /trunk/Mars/Changelog =================================================================== --- /trunk/Mars/Changelog (revision 10092) +++ /trunk/Mars/Changelog (revision 10093) @@ -18,4 +18,24 @@ -*-*- END OF LINE -*-*- + + 2011/01/07 Thomas Bretz + + * msimcamera/MSimAPD.[h,cc]: + - added setting of the afterpulse probabilities + - handle afterpulses in processing + + * melectronics/MAvalanchePhotoDiode.[h,cc]: + - restructured header file + - added treatment of afterpulses in HitCellImp + - added time constants and afterpulse probability + data memebers + - added TSortedList to hold the afterpulses + - added afterpulses to calculation of relaxation time + - added call to Relax to Init + - added treatment of afterpulses to FillEmpty, FillRandom and + Evolve + - added member functions to process the afterpulses + + 2011/01/06 Thomas Bretz Index: /trunk/Mars/melectronics/MAvalanchePhotoDiode.cc =================================================================== --- /trunk/Mars/melectronics/MAvalanchePhotoDiode.cc (revision 10092) +++ /trunk/Mars/melectronics/MAvalanchePhotoDiode.cc (revision 10093) @@ -39,4 +39,50 @@ // Hamamatsu 60x60 cells: APD(60, 0.2, 3, 8.75) // +// +// The implementation of afterpulsing is based on +// A.Du, F.Retiere +// After-pulsing and cross-talk in multi-pixel photon counters +// NIM A, Volume 596, Issue 3, p. 396-401 +// +// +// Example: +// +// APD apd(ncells, crosstalk, deadtime, recovery); +// apd.SetAfterpulseProb(0.14, 0.11); +// +// while (1) +// { +// // Make the chip "empty" from the influence of external photons +// // It also sets fTime to 0. +// apd.Init(freq); // freq of external noise (eg. nsb) +// +// // Now call this function for each external photon you have. The +// // times are relative to the the time you get by apd.GetTime() +// // set automatically after the call to apd.Init(). +// for (int i=0; iGetListOfAfterpulses()); +// Afterpulse *ap = 0; +// while ((ap=static_cast(Next()))) +// { +// if (apd.GetTime()>=dtend) +// continue; +// +// cout << "Amplitude: " << ap->GetAmplitude() << endl; +// cout << "Arrival Time: " << ap->GetTime() << endl; +// } +// } +// +// ////////////////////////////////////////////////////////////////////////////// #include "MAvalanchePhotoDiode.h" @@ -45,4 +91,7 @@ #include "MMath.h" + +#include "MLog.h" +#include "MLogManip.h" ClassImp(APD); @@ -77,21 +126,43 @@ APD::APD(Int_t n, Float_t prob, Float_t dt, Float_t rt) : fHist("APD", "", n, 0.5, n+0.5, n, 0.5, n+0.5), - fCrosstalkProb(prob), fDeadTime(dt), fRecoveryTime(rt), fTime(-1) + fCrosstalkProb(prob), fDeadTime(dt), fRecoveryTime(rt), + fTime(-1) { fHist.SetDirectory(0); -} - -// -------------------------------------------------------------------------- -// -// This is the time the cell needs after a hit until less than threshold -// (0.001 == one permille) of the signal is lost. -// -// If all cells of a G-APD have fired simultaneously and they are fired -// once more after the relaxation time (with the defaultthreshold of 0.001) -// the chip will show a signal which is one permille too small. + + fAfterpulses.SetOwner(); + + fAfterpulseProb[0] = 0; + fAfterpulseProb[1] = 0; + + fAfterpulseTau[0] = 15; + fAfterpulseTau[1] = 85; +} + +// -------------------------------------------------------------------------- +// +// This is the time a chips needs after an external signal to relax to +// a "virgin" state, i.e. without no influence of the external pulse +// above the given threshold. +// +// It takes into account the dead time of the cell, the relaxation time +// and the two afterpulse distributions. However, in most cases the +// afterpulse distribution will dominate (except they are switched off by +// a zero probability). +// +// FIXME: Maybe the calculation of the relaxation time could be optimized? // Float_t APD::GetRelaxationTime(Float_t threshold) const { - return fDeadTime-TMath::Log(threshold)*fRecoveryTime; + // Calculate time until the probability of one of these + // events falls below the threshold probability. + const Double_t rt0 = - TMath::Log(threshold)*fRecoveryTime; + const Double_t rt1 = fAfterpulseProb[0]>0 ? -TMath::Log(threshold/fAfterpulseProb[0])*fAfterpulseTau[0] : 0; + const Double_t rt2 = fAfterpulseProb[0]>0 ? -TMath::Log(threshold/fAfterpulseProb[1])*fAfterpulseTau[1] : 0; + + // Probability not between t and inf, but between t and t+dt + // -tau * log ( p / ( 1 - exp(- dt/tau) ) ) = t + + return fDeadTime + TMath::Max(rt0, TMath::Max(rt1, rt2)); } @@ -108,4 +179,8 @@ // probability. // +// For each photon the possible afterpulses of two distributions are +// created and added to the list of afterpulses. This is done by calling +// GenerateAfterpulse for the two afterpulse-distributions. +// // The total height of the signal (in units of photons) is returned. // Note, that this can be a fractional number. @@ -122,4 +197,7 @@ // y<1 || y>fHist.GetNbinsY()) // return 0; +#ifdef DEBUG + cout << "Hit: " << t << endl; +#endif // Number of the x/y cell in the one dimensional array @@ -137,9 +215,18 @@ // Photons within the dead time are just ignored if (/*hx.GetBinContent(x,y)>0 &&*/ dt<=0) + { +#ifdef DEBUG + cout << "Dead: " << t << " " << cont << " " << dt << endl; +#endif return 0; - + } // The signal height (in units of one photon) produced after dead time // depends on the recovery of the cell - described by an exponential. const Float_t weight = fRecoveryTime<=0 ? 1. : 1-TMath::Exp(-dt/fRecoveryTime); + + // Now we know the charge in the cell and we can generate + // the afterpulses with both time-constants + GenerateAfterpulse(cell, 0, weight, t); + GenerateAfterpulse(cell, 1, weight, t); // The probability that the cell emits a photon causing crosstalk @@ -201,5 +288,5 @@ // // Check if x and y is a valid cell. If not return 0, otherwise -// HitCelImp(x, y, t) +// HitCelImp(x, y, t). HitCellImp generates Crosstalk and Afterpulses. // // The default time is 0. @@ -217,5 +304,6 @@ // // Determine randomly (uniformly) a cell which was hit. Return -// HitCellImp for this cell and the given time. +// HitCellImp for this cell and the given time. HitCellImp +// generates Crosstalk and Afterpulses // // The default time is 0. @@ -238,9 +326,11 @@ // -------------------------------------------------------------------------- // -// Sets all cells with a contents whihc is well before the time t such that +// Sets all cells with a contents which is well before the time t such that // the chip is "virgin". Therefore all cells are set to a time which // is twice the deadtime before the given time and 1000 times the recovery // time. // +// The afterpulse list is deleted. +// // If deadtime and recovery time are 0 then t-1 is set. // @@ -260,4 +350,6 @@ fHist.SetEntries(1); + fAfterpulses.Delete(); + fTime = t; } @@ -267,9 +359,19 @@ // First call FillEmpty for the given time t. Then fill each cell by // by calling HitCellImp with time t-gRandom->Exp(n/rate) with n being -// the total number of cells. +// the total number of cells. This the time at which the cell was last hit. // // Sets fTime to t // -// The default time is 0. +// If the argument t is omitted it defaults to 0. +// +// Since after calling this function the chip should reflect the +// status at the new time fTime=t, all afterpulses are processed +// until this time. However, the produced random pulses might have produced +// new new afterpulses. +// +// All afterpulses before the new timestamp are deleted. +// +// WARNING: Note that this might not correctly reproduce afterpulses +// produced by earlier pulese. // void APD::FillRandom(Float_t rate, Float_t t) @@ -282,5 +384,6 @@ const Double_t f = (nx*ny)/rate; - // FIXME: This is not perfect, is it? What about the dead time? + // FIXME: Dead time is not taken into account, + // possible earlier afterpulses are not produced. for (int x=1; x<=nx; x++) @@ -288,5 +391,35 @@ HitCellImp(x, y, t-MMath::RndmExp(f)); + // Deleting of the afterpulses before fHist.GetMinimum() won't + // speed things because we have to loop over them once in any case + + ProcessAfterpulses(fHist.GetMinimum(), t); + DeleteAfterpulses(t); + fTime = t; +} + +// -------------------------------------------------------------------------- +// +// Shift all times including fTime to dt (ie. add -dt to all times) +// This allows to set a user-defined T0 or shift T0 to fTime=0. +// +// However, T0<0 is not allowed (dt cannot be greater than fTime) +// +void APD::ShiftTime(Double_t dt) +{ + if (dt>fTime) + { + gLog << warn << "APD::ShiftTime: WARNING - requested time shift results in fTime<0... ignored." << endl; + return; + } + + // If reset was requested shift all times by end backwards + // so that fTime is now 0 + const Int_t n = (fHist.GetNbinsX()+2)*(fHist.GetNbinsY()+2); + for (int i=0; i0) { - time += MMath::RndmExp(avglen); - if (time>end) - break; - - hits += HitRandomCell(time); + const Double_t avglen = 1./freq; + + Double_t time = fTime; + while (1) + { + const Double_t deltat = MMath::RndmExp(avglen); + if (time+deltat>end) + break; + + time += deltat; + + hits += HitRandomCell(time); + } } + + // Deleting of the afterpulses before fTime won't speed things + // because we have to loop over them once in any case + + ProcessAfterpulses(fTime, dt); + DeleteAfterpulses(end); fTime = end; @@ -329,4 +480,8 @@ // dead. // The default time is 0. +// +// Note that if you want to get a correct measure of teh number of dead cells +// at the time t, this function will only produce a valid count if the +// afterpulses have been processed up to this time. // Int_t APD::CountDeadCells(Float_t t) const @@ -348,4 +503,8 @@ // Returs the number of cells which have a time t<=fDeadTime+fRecoveryTime. // The default time is 0. +// +// Note that if you want to get a correct measure of teh number of dead cells +// at the time t, this function will only produce a valid count if the +// afterpulses have been processed up to this time. // Int_t APD::CountRecoveringCells(Float_t t) const @@ -364,2 +523,109 @@ return n; } + +// -------------------------------------------------------------------------- +// +// Generate an afterpulse originating from the given cell and a pulse with +// charge. The afterpulse distribution to use is specified by +// the index. The "current" time to which the afterpulse delay refers must +// be given by t. +// +// A generated Afterpulse is added to the list of afterpulses +// +void APD::GenerateAfterpulse(UInt_t cell, Int_t idx, Double_t charge, Double_t t) +{ + // The cell had a single avalanche with signal height weight. + // This cell now can produce an afterpulse photon/avalanche. + const Double_t p = gRandom->Uniform(); + + // It's probability scales with the charge of the pulse + if (p>charge*fAfterpulseProb[idx]) + return; + + // Afterpulses come with a well defined time-constant + // after the normal pulse + const Double_t dt = MMath::RndmExp(fAfterpulseTau[idx]); + + fAfterpulses.Add(new Afterpulse(cell, t+dt)); + +#ifdef DEBUG + cout << "Add : " << t << " + " << dt << " = " << t+dt << endl; +#endif +} + +// -------------------------------------------------------------------------- +// +// Process afterpulses between time and time+dt. All afterpulses in the list +// before t=time are ignored. All afterpulses between t=time and +// t=time+dt are processed through HitCellImp. Afterpulses after and +// equal t=time+dt are skipped. +// +// Since the afterpulse list is a sorted list newly generated afterpulses +// are always inserted into the list behind the current entry. Consequently, +// afterpulses generated by afterpulses will also be processed correctly. +// +// Afterpulses with zero amplitude are deleted from the list. All other after +// pulses remain in the list for later evaluation. +// +void APD::ProcessAfterpulses(Float_t time, Float_t dt) +{ +#ifdef DEBUG + cout << "Process afterpulses from " << time << " to " << time+dt << endl; +#endif + + const Float_t end = time+dt; + + TObjLink *lnk = fAfterpulses.FirstLink(); + while (lnk) + { + Afterpulse &ap = *static_cast(lnk->GetObject()); + + // Skip afterpulses which have been processed already + // or which we do not have to process anymore + if (ap.GetTime()