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

/* ======================================================================== *\
!
! *
! * This file is part of CheObs, the Modular Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appears in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Thomas Bretz,  1/2009 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: CheObs Software Development, 2000-2009
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//  APD
//
// All times in this class are relative times. Therefor the unit for the
// time is not intrinsically fixed. In fact the dead-time and recovery-
// time given in the constructor must have the same units. This is what
// defines the unit of the times given in the function and the unit of
// rates given.
// For example, if recovery and dead time are given in nanoseconds the
// all times must be in nanoseconds and rates are given per nanosecond,
// i.e. GHz.
//
//  Hamamatsu 30x30 cells:  APD(30, 0.2, 3, 35)
//  Hamamatsu 60x60 cells:  APD(60, 0.2, 3, 8.75)
//
//////////////////////////////////////////////////////////////////////////////
#include "MAvalanchePhotoDiode.h"

#include <TRandom.h>

#include "MMath.h"

ClassImp(APD);

using namespace std;

/*
class MyProfile : public TProfile2D
{
public:
    void AddBinEntry(Int_t cell) { fBinEntries.fArray[cell]++; }
};
*/

// --------------------------------------------------------------------------
//
//  Default Constructor.
//
//    n     is the number od cells in x or y. The APD is assumed to
//          be square.
//    prob  is the crosstalk probability, i.e., the probability that a
//          photon which produced an avalanche will create another
//          photon in a neighboring cell
//    dt    is the deadtime, i.e., the time in which the APD cell will show
//          no response to a photon after a hit
//    rt    is the recovering tims, i.e. the exponential (e^(-dt/rt))
//          with which the cell is recovering after being dead
//
//    prob, dt and ar can be set to 0 to switch the effect off.
//    0 is also the dfeault for all three.
//
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)
{
    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.
//
Float_t APD::GetRelaxationTime(Float_t threshold) const
{
    return fDeadTime-TMath::Log(threshold)*fRecoveryTime;
}

// --------------------------------------------------------------------------
//
// This is the recursive implementation of a hit. If a photon hits a cell
// at x and y (must be a valid cell!) at time t, at first we check if the
// cell is still dead. If it is not dead we calculate the signal height
// from the recovery time. Now we check with the crosstalk probability
// whether another photon is created. If another photon is created we
// calculate randomly which of the four neighbor cells are hit.
// If the cell is outside the APD the photon is ignored. As many
// new photons are created until our random number is below the crosstak-
// probability.
//
// The total height of the signal (in units of photons) is returned.
// Note, that this can be a fractional number.
//
// This function looks a bit fancy accessing the histogram and works around
// a few histogram functions. This is a speed optimization which works
// around a lot of sanity checks which are obsolete in our case.
//
// The default time is 0.
//
Float_t APD::HitCellImp(Int_t x, Int_t y, Float_t t)
{
    //        if (x<1 || x>fHist.GetNbinsX() ||
    //            y<1 || y>fHist.GetNbinsY())
    //            return 0;

    // const Int_t cell = fHist.GetBin(x, y);
    const Int_t cell = x + (fHist.GetNbinsX()+2)*y;

    // This is the fastes way to access the bin-contents in fArray
    Float_t &cont = fHist.GetArray()[cell];

    // const Double_t dt = t-fHist.GetBinContent(x, y)-fDeadTime; //
    const Float_t dt = t-cont-fDeadTime;

    // Photons within the dead time are just ignored
    if (/*hx.GetBinContent(x,y)>0 &&*/ dt<=0)
        return 0;

    // Signal height (in units of one photon) produced after dead time
    const Float_t weight = fRecoveryTime<=0 ? 1 : 1.-exp(-dt/fRecoveryTime);

    cont = t; // fHist.SetBinContent(x, y, t)

    // Counter for the numbers of produced photons
    Float_t n = weight;

    /*
     // Check if a photon in a neighboring cell is produced (crosstalk)
     while (gRandom->Rndm()<fCrosstalkProb)
     {
        // Get a random neighbor which is hit.
        switch (gRandom->Integer(4))
        {
        case 0: x++; if (x>fHist.GetNbinsX()) continue; break;
        case 1: x--; if (x<1) continue;                 break;
        case 2: y++; if (y>fHist.GetNbinsY()) continue; break;
        case 3: y--; if (y<1) continue;                 break;
        }

        n += HitCellImp(x, y, t);
     }
     */

    //for (int i=0; i<1; i++)
    while (1)
    {
        const Double_t rndm = gRandom->Rndm();
        if (rndm>=fCrosstalkProb)
            break;

        // We can re-use the random number becuase it is uniformely
        // distributed. This saves cpu power
        const Int_t dir = TMath::FloorNint(4*rndm/fCrosstalkProb);

        // Get a random neighbor which is hit.
        switch (dir)
        {
        case 0: if (x<fHist.GetNbinsX()) n += HitCellImp(x+1, y, t); break;
        case 1: if (x>1)                 n += HitCellImp(x-1, y, t); break;
        case 2: if (y<fHist.GetNbinsY()) n += HitCellImp(x, y+1, t); break;
        case 3: if (y>1)                 n += HitCellImp(x, y-1, t); break;
        }

        // In the unlikely case the calculated direction is out-of-range,
        // i.e. <0 or >3, we would just try to fill the same cell again which
    }

    return n;
}

// --------------------------------------------------------------------------
//
// Check if x and y is a valid cell. If not return 0, otherwise
// HitCelImp(x, y, t)
//
// The default time is 0.
//
Float_t APD::HitCell(Int_t x, Int_t y, Float_t t)
{
    if (x<1 || x>fHist.GetNbinsX() ||
        y<1 || y>fHist.GetNbinsY())
        return 0;

    return HitCellImp(x, y, t);
}

// --------------------------------------------------------------------------
//
// Determine randomly (uniformly) a cell which was hit. Return
// HitCellImp for this cell and the given time.
//
// The default time is 0.
//
// If you want t w.r.t. fTime use HitRandomCellRelative istead.
//
Float_t APD::HitRandomCell(Float_t t)
{
    const UInt_t nx  = fHist.GetNbinsX();
    const UInt_t ny  = fHist.GetNbinsY();

    const UInt_t idx = gRandom->Integer(nx*ny);

    const UInt_t x   = idx%nx;
    const UInt_t y   = idx/nx;

    return HitCellImp(x+1, y+1, t);
}

// --------------------------------------------------------------------------
//
// Sets all cells with a contents whihc 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.
//
// If deadtime and recovery time are 0 then t-1 is set.
//
// Sets fTime to t
//
// The default time is 0.
//
void APD::FillEmpty(Float_t t)
{
    const Int_t n = (fHist.GetNbinsX()+2)*(fHist.GetNbinsY()+2);

    const Double_t tm = fDeadTime<=0 && fRecoveryTime<=0 ? t-1 : t-2*fDeadTime-1000*fRecoveryTime;

    for (int i=0; i<n; i++)
        fHist.GetArray()[i] = tm;

    fHist.SetEntries(1);

    fTime = t;
}

// --------------------------------------------------------------------------
//
// 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.
//
// Sets fTime to t
//
// The default time is 0.
//
void APD::FillRandom(Float_t rate, Float_t t)
{
    FillEmpty(t);

    const Int_t nx = fHist.GetNbinsX();
    const Int_t ny = fHist.GetNbinsY();

    const Double_t f = (nx*ny)/rate;

    // FIXME: This is not perfect, is it? What about the dead time?

    for (int x=1; x<=nx; x++)
        for (int y=1; y<=ny; y++)
            HitCellImp(x, y, t-MMath::RndmExp(f));

    fTime = t;
}

// --------------------------------------------------------------------------
//
// Evolve the chip from fTime to fTime+dt if it with a given frequency
// freq. Returns the total signal "recorded".
//
// fTime is set to the fTime+dt
//
// If you want to evolve over a default relaxation time (relax the chip
// from a signal) use Relax instead.
//
Float_t APD::Evolve(Double_t freq, Double_t dt)
{
    const Double_t avglen = 1./freq;

    const Double_t end  = fTime+dt;

    Float_t hits = 0;

    Double_t time = fTime;
    while (1)
    {
        time += MMath::RndmExp(avglen);
        if (time>end)
            break;

        hits += HitRandomCell(time);
    }

    fTime = end;

    return hits;
}

// --------------------------------------------------------------------------
//
// Retunrs the number of cells which have a time t<=fDeadTime, i.e. which are
// dead.
// The default time is 0.
//
Int_t APD::CountDeadCells(Float_t t) const
{
    const Int_t nx = fHist.GetNbinsX();
    const Int_t ny = fHist.GetNbinsY();

    Int_t n=0;
    for (int x=1; x<=nx; x++)
        for (int y=1; y<=ny; y++)
            if ((t-fHist.GetBinContent(x, y))<=fDeadTime)
                n++;

    return n;
}

// --------------------------------------------------------------------------
//
// Returs the number of cells which have a time t<=fDeadTime+fRecoveryTime.
// The default time is 0.
//
Int_t APD::CountRecoveringCells(Float_t t) const
{
    const Int_t nx = fHist.GetNbinsX();
    const Int_t ny = fHist.GetNbinsY();

    Int_t n=0;
    for (int x=1; x<=nx; x++)
        for (int y=1; y<=ny; y++)
        {
            Float_t dt = t-fHist.GetBinContent(x, y);
            if (dt>fDeadTime && dt<=fDeadTime+fRecoveryTime)
                n++;
        }
    return n;
}