Changeset 10093


Ignore:
Timestamp:
01/07/11 15:48:35 (12 years ago)
Author:
tbretz
Message:
Added afterpulse treatment to MAvalanchePhotoDiode and MSimAPD
Location:
trunk/Mars
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/Mars/Changelog

    r10092 r10093  
    1818
    1919                                                 -*-*- END OF LINE -*-*-
     20
     21 2011/01/07 Thomas Bretz
     22
     23   * msimcamera/MSimAPD.[h,cc]:
     24     - added setting of the afterpulse probabilities
     25     - handle afterpulses in processing
     26
     27   * melectronics/MAvalanchePhotoDiode.[h,cc]:
     28     - restructured header file
     29     - added treatment of afterpulses in HitCellImp
     30     - added time constants and afterpulse probability
     31       data memebers
     32     - added TSortedList to hold the afterpulses
     33     - added afterpulses to calculation of relaxation time
     34     - added call to Relax to Init
     35     - added treatment of afterpulses to FillEmpty, FillRandom and
     36       Evolve
     37     - added member functions to process the afterpulses
     38
     39
    2040
    2141 2011/01/06 Thomas Bretz
  • trunk/Mars/melectronics/MAvalanchePhotoDiode.cc

    r9565 r10093  
    3939//  Hamamatsu 60x60 cells:  APD(60, 0.2, 3, 8.75)
    4040//
     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//
    4187//////////////////////////////////////////////////////////////////////////////
    4288#include "MAvalanchePhotoDiode.h"
     
    4591
    4692#include "MMath.h"
     93
     94#include "MLog.h"
     95#include "MLogManip.h"
    4796
    4897ClassImp(APD);
     
    77126APD::APD(Int_t n, Float_t prob, Float_t dt, Float_t rt)
    78127    : fHist("APD", "", n, 0.5, n+0.5, n, 0.5, n+0.5),
    79     fCrosstalkProb(prob), fDeadTime(dt), fRecoveryTime(rt), fTime(-1)
     128    fCrosstalkProb(prob), fDeadTime(dt), fRecoveryTime(rt),
     129    fTime(-1)
    80130{
    81131    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.
     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?
    92154//
    93155Float_t APD::GetRelaxationTime(Float_t threshold) const
    94156{
    95     return fDeadTime-TMath::Log(threshold)*fRecoveryTime;
     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));
    96167}
    97168
     
    108179// probability.
    109180//
     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//
    110185// The total height of the signal (in units of photons) is returned.
    111186// Note, that this can be a fractional number.
     
    122197    //            y<1 || y>fHist.GetNbinsY())
    123198    //            return 0;
     199#ifdef DEBUG
     200    cout << "Hit: " << t << endl;
     201#endif
    124202
    125203    // Number of the x/y cell in the one dimensional array
     
    137215    // Photons within the dead time are just ignored
    138216    if (/*hx.GetBinContent(x,y)>0 &&*/ dt<=0)
     217    {
     218#ifdef DEBUG
     219        cout << "Dead: " << t << " " << cont << " " << dt << endl;
     220#endif
    139221        return 0;
    140 
     222    }
    141223    // The signal height (in units of one photon) produced after dead time
    142224    // depends on the recovery of the cell - described by an exponential.
    143225    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);
    144231
    145232    // The probability that the cell emits a photon causing crosstalk
     
    201288//
    202289// Check if x and y is a valid cell. If not return 0, otherwise
    203 // HitCelImp(x, y, t)
     290// HitCelImp(x, y, t). HitCellImp generates Crosstalk and Afterpulses.
    204291//
    205292// The default time is 0.
     
    217304//
    218305// Determine randomly (uniformly) a cell which was hit. Return
    219 // HitCellImp for this cell and the given time.
     306// HitCellImp for this cell and the given time. HitCellImp
     307// generates Crosstalk and Afterpulses
    220308//
    221309// The default time is 0.
     
    238326// --------------------------------------------------------------------------
    239327//
    240 // Sets all cells with a contents whihc is well before the time t such that
     328// Sets all cells with a contents which is well before the time t such that
    241329// the chip is "virgin". Therefore all cells are set to a time which
    242330// is twice the deadtime before the given time and 1000 times the recovery
    243331// time.
    244332//
     333// The afterpulse list is deleted.
     334//
    245335// If deadtime and recovery time are 0 then t-1 is set.
    246336//
     
    260350    fHist.SetEntries(1);
    261351
     352    fAfterpulses.Delete();
     353
    262354    fTime = t;
    263355}
     
    267359// First call FillEmpty for the given time t. Then fill each cell by
    268360// by calling HitCellImp with time t-gRandom->Exp(n/rate) with n being
    269 // the total number of cells.
     361// the total number of cells. This the time at which the cell was last hit.
    270362//
    271363// Sets fTime to t
    272364//
    273 // The default time is 0.
     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.
    274376//
    275377void APD::FillRandom(Float_t rate, Float_t t)
     
    282384    const Double_t f = (nx*ny)/rate;
    283385
    284     // FIXME: This is not perfect, is it? What about the dead time?
     386    // FIXME: Dead time is not taken into account,
     387    //        possible earlier afterpulses are not produced.
    285388
    286389    for (int x=1; x<=nx; x++)
     
    288391            HitCellImp(x, y, t-MMath::RndmExp(f));
    289392
     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
    290399    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;
    291424}
    292425
     
    296429// freq. Returns the total signal "recorded".
    297430//
    298 // fTime is set to the fTime+dt
     431// fTime is set to the fTime+dt.
    299432//
    300433// If you want to evolve over a default relaxation time (relax the chip
    301434// from a signal) use Relax instead.
    302435//
     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//
    303443Float_t APD::Evolve(Double_t freq, Double_t dt)
    304444{
    305     const Double_t avglen = 1./freq;
    306 
    307     const Double_t end  = fTime+dt;
     445    const Double_t end = fTime+dt;
    308446
    309447    Float_t hits = 0;
    310 
    311     Double_t time = fTime;
    312     while (1)
     448 
     449    if (freq>0)
    313450    {
    314         time += MMath::RndmExp(avglen);
    315         if (time>end)
    316             break;
    317 
    318         hits += HitRandomCell(time);
     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        }
    319464    }
     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);
    320471
    321472    fTime = end;
     
    329480// dead.
    330481// 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.
    331486//
    332487Int_t APD::CountDeadCells(Float_t t) const
     
    348503// Returs the number of cells which have a time t<=fDeadTime+fRecoveryTime.
    349504// 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.
    350509//
    351510Int_t APD::CountRecoveringCells(Float_t t) const
     
    364523    return n;
    365524}
     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}
  • trunk/Mars/msimcamera/MSimAPD.cc

    r9518 r10093  
    3737//
    3838// Remark:
    39 //   - The photons MUST be sorted increasing in time.
    4039//   - The photon rate used to initialize the APD must match the one used
    4140//     to "fill" the random photons. (FIXME: This should be stored somewhere)
     
    152151    //   * Make the arguments a data member of MSimAPD
    153152
    154     Int_t   ncells    = 0;
    155     Float_t crosstalk = 0;
    156     Float_t deadtime  = 0;
    157     Float_t recovery  = 0;
     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;
    158159
    159160    switch (fType)
    160161    {
    161162    case 1:
    162         ncells    = 30;
    163         crosstalk = 0.2;
    164         deadtime  = 3;
    165         recovery  = 8.75*4;
     163        ncells     = 30;
     164        crosstalk  = 0.2;
     165        deadtime   = 3;
     166        recovery   = 8.75*4;
     167        afterprob1 = 0;
     168        afterprob2 = 0;
    166169        break;
    167170
    168171    case 2:
    169         ncells    = 60;
    170         crosstalk = 0.2;
    171         deadtime  = 3;
    172         recovery  = 8.75;
     172        ncells     = 60;
     173        crosstalk  = 0.2;
     174        deadtime   = 3;
     175        recovery   = 8.75;
     176        afterprob1 = 0;
     177        afterprob2 = 0;
    173178        break;
    174179
    175180    case 3:
    176         ncells    = 60;
    177         crosstalk = 0.15;
    178         deadtime  = 3;
    179         recovery  = 8.75;
     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;
    180196        break;
    181197
     
    186202
    187203    for (UInt_t i=0; i<fGeom->GetNumPixels(); i++)
    188         fAPDs.Add(new APD(ncells, crosstalk, deadtime, recovery));
     204    {
     205        APD *apd = new APD(ncells, crosstalk, deadtime, recovery);
     206        apd->SetAfterpulseProb(afterprob1, afterprob2);
     207
     208        fAPDs.Add(apd);
     209    }
    189210
    190211    return kTRUE;
     
    198219Int_t MSimAPD::Process()
    199220{
    200     //const Double_t rate = 40e9;
    201     // FIXME: Where do we get this number from??
    202     // const Double_t rate = 0.04;
    203 
    204221    // Make all APDs look neutral for the first hit by a photon according to the
    205222    // average hit rate
     
    213230        return kERROR;
    214231    }
    215 /*
    216     for (UInt_t idx=0; idx<npix; idx++)
    217     {
    218         const Double_t freq = fRates ? (*fRates)[idx].GetPedestal() : fFreq;
    219         static_cast<APD*>(fAPDs.UncheckedAt(idx))->FillRandom(freq, fStat->GetTimeFirst());
    220     }
    221 */
     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();
    222238
    223239    // This tries to initialize dead and relaxing cells properly. If
     
    228244    {
    229245        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.
    230254        static_cast<APD*>(fAPDs.UncheckedAt(idx))->Init(freq);
    231255    }
     
    240264        MPhotonData &ph = (*fEvt)[i];
    241265
    242         // Get arrival time of photon and idx
     266        // Get arrival time of photon wrt to left edge of window and its index
    243267        const Double_t t = ph.GetTime()-fStat->GetTimeFirst();
    244268        const Int_t  idx = ph.GetTag();
     
    255279            return kERROR;
    256280        }
    257         // Simulate hitting the APD (the signal height in effective
    258         // "number of photons" is returned)
     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
    259287        const Double_t hits = static_cast<APD*>(fAPDs.UncheckedAt(idx))->HitRandomCellRelative(t);
    260 
    261         // FIXME: Make a proper simulation of the excess noise!!!
    262         //const Double_t signal = gRandom->Gaus(hits, 0.2*TMath::Sqrt(hits));
    263288
    264289        // Set the weight to the input
     
    269294    // simulated time.
    270295    for (UInt_t idx=0; idx<npix; idx++)
    271         static_cast<APD*>(fAPDs.UncheckedAt(idx))->IncreaseTime(fStat->GetTimeLast());
     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();
    272342
    273343    return kTRUE;
Note: See TracChangeset for help on using the changeset viewer.