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.uniwuerzburg.de>


19  !


20  ! Copyright: CheObs Software Development, 20002009


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 deadtime 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  // Afterpulsing and crosstalk in multipixel photon counters


45  // NIM A, Volume 596, Issue 3, p. 396401


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 


97  ClassImp(APD);


98 


99  using namespace std;


100 


101  /*


102  class MyProfile : public TProfile2D


103  {


104  public:


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  //


126  APD::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  //


155  Float_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 afterpulsedistributions.


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  //


194  Float_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 bincontents in fArray


209  Float_t &cont = fHist.GetArray()[cell];


210 


211  // Calculate the time since the last breakdown


212  // const Double_t dt = tfHist.GetBinContent(x, y)fDeadTime; //


213  const Float_t dt = tcontfDeadTime;


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. : 1TMath::Exp(dt/fRecoveryTime);


226 


227  // Now we know the charge in the cell and we can generate


228  // the afterpulses with both timeconstants


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 reuse 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(x1, 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, y1, t); break;


278  }


279 


280  // In the unlikely case the calculated direction is outofrange,


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  //


294  Float_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  //


313  Float_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 t1 is set.


336  //


337  // Sets fTime to t


338  //


339  // The default time is 0.


340  //


341  void 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 ? t1 : t2*fDeadTime1000*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 tgRandom>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  //


377  void 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, tMMath::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 userdefined T0 or shift T0 to fTime=0.


406  //


407  // However, T0<0 is not allowed (dt cannot be greater than fTime)


408  //


409  void 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  //


443  Float_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  //


487  Int_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 ((tfHist.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  //


510  Int_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 = tfHist.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  //


535  void 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 timeconstant


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  //


570  void 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  //


618  void 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  }

