source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/detector/electronics/src/Photomultiplier.cc @ 114

// implementation for class Photomultiplier
// M. Pallavicini
//
#include "Photomultiplier.hh"
#include "PmtGeometry.hh"
#include "FrontEndChip.hh"
#include "MacroCell.hh"
#include <unistd.h>
#include <math.h>
#include <iostream>
#include "Config.hh"
#include "EsafRandom.hh"
#include "EEvent.hh"
#include "Etypes.hh"
#include "EusoMapping.hh"
#include "EDetectorPhotonAdder.hh"
#include "EusoElectronics.hh"
#include "EusoDetector.hh"
#include "EConst.hh"

ClassImp(Photomultiplier)

Double_t Photomultiplier::fgQuantum = 0;           // quantum efficiency
Double_t Photomultiplier::fgGain = 0;              // nominal gain
Double_t Photomultiplier::fgGainSigma = 0;         // gain spread
Double_t Photomultiplier::fgWidth = 0;             // time width
Double_t Photomultiplier::fgDarkNoiseRate = 0;     // dark noise rate
Double_t Photomultiplier::fgGtuLength = 0;         // gtu length in ns

//______________________________________________________________________________
Photomultiplier::Photomultiplier(Int_t id, PmtGeometry* g ): fId(id),
    fFrontEnd(NULL), fGeometry(g), fCell(NULL), fEC(NULL) {
    //
    // Constructor
    //

    Geometry()->SetPmt( this );


    SetState(kPmtIdle);
    SetEmpty( true );

    if ( fgQuantum == 0 ) {
        fgQuantum = Conf()->GetNum("Photomultiplier.PmtQuantum");
        fgGain = Conf()->GetNum("Photomultiplier.PmtGain");
        fgGainSigma = Conf()->GetNum("Photomultiplier.PmtGainSigma");
        fgWidth = Conf()->GetNum("Photomultiplier.PmtTimeWidth");
        fgDarkNoiseRate = Conf()->GetNum("Photomultiplier.fDarkNoiseRate")/sou::microsecond;
        fgGtuLength = Config::Get()->GetCF("Electronics","MacroCell")->GetNum("MacroCell.fGtuTimeLength");
    }
}


//______________________________________________________________________________
Photomultiplier::~Photomultiplier() {
    //
    // Destructor
    //

    SafeDelete(fGeometry);
    ClearMemory();
}

//______________________________________________________________________________
Bool_t Photomultiplier::ClearMemory() {
    //
    // Physically release the memory allocated by the arrays of this object
    //

    // The frontend holds a list of pointers to the vectors of pmt signals.
    // They are going to be deleted and the pointers to be invalid. The
    // frontend must be reseted
    FrontEnd()->Reset();

    map<Int_t, vector<PmtSignal*>* >::const_iterator itVSig;

    for(itVSig=fPmtHits.begin(); itVSig != fPmtHits.end(); itVSig++)  {
        delete itVSig->second;
    }

    fPmtHits.clear();



    return kTRUE;

}

//______________________________________________________________________________
void Photomultiplier::Reset() {
    // Reset pmt and get ready for next event list are deleted
    // PmtSignal objects are deleted by Front End


    map<Int_t, vector<PmtSignal*>* >::const_iterator itVSig;
    vector<PmtSignal*>::const_iterator itSig;

    for(itVSig=fPmtHits.begin(); itVSig != fPmtHits.end(); itVSig++)  {
        if( itVSig->second ) {
            vector<PmtSignal*>* ptrVSig = itVSig->second;
            for(itSig = ptrVSig->begin(); itSig != ptrVSig->end(); itSig++)
                delete *itSig ;

            ptrVSig->clear();

        }
    }
//    for ( Int_t ch=0; ch < Geometry()->NumPads(); ch++ ) {
//        if ( fPmtHits[ch] ) {
//            fPmtHits[ch]->clear();
//        }
//    }
    SetState(kPmtIdle);
    fStartTime = HUGE;
    fEndTime = -HUGE;
    SetEmpty( true );
}


//______________________________________________________________________________
void Photomultiplier::ResetClass() {
    //
    // Reset static members of the class
    //

    fgQuantum = 0;
    fgGain = 0;
    fgGainSigma = 0;
    fgWidth = 0;
    fgGtuLength = 0;
}

//______________________________________________________________________________
Int_t Photomultiplier::FEChannel(Int_t n) {
    //
    // Returns front end channel number corresponding to pmt channel n
    //

    if ( isValid(n) )
        return fFeMap[n];
    else
        return -1;
}

//______________________________________________________________________________
ChannelUniqueId Photomultiplier::GetUniqueId(Int_t ch) const {
    //
    // Returns channel unique id if the channel belongs to this pmt
    //

    if ( isValid(ch))
        return Geometry()->GetUniqueId(ch);
    else
        MsgForm(EsafMsg::Panic, "Photomultiplier: %d out of range",ch);
        return -1;
}

//______________________________________________________________________________
void Photomultiplier::SetFrontEnd( FrontEndChip* fe, Int_t row, Int_t col ) {
    //
    // Re-assign a new front end chip to this pmt
    //

    if ( fe == NULL ) {
        Msg(EsafMsg::Warning) << "NULL front end pointer in Photomultiplier::SetFrontEnd()" << MsgDispatch;
        return;
    }
    if ( fFrontEnd ) {
        Msg(EsafMsg::Warning)  << "Front End Chip " << fFrontEnd->Id()
            << " will be destroyed in Photomultiplier::SetFrontEnd()" << MsgDispatch;
        delete fFrontEnd;
    }
    fFrontEnd = fe;
    BuildFrontEndMap(row, col);
}

//______________________________________________________________________________
void Photomultiplier::BuildFrontEndMap( Int_t row, Int_t col) {
    // Connection between PMT channels and Front-End channels
    // this is non trivial where a single front end chip reads more than one PMT
    // the other possibility (pmt channels exceeding fe channels) is not implemented

    Int_t r,c;

    for ( Int_t i=0; i<Geometry()->NumPads(); i++) {
        r = row+(i / Geometry()->Rows());
        c = col+(i % Geometry()->Rows());

        if ( fFrontEnd->ChanRowCol(r,c) < FrontEnd()->Channels() ) {

            fFeMap[i] = fFrontEnd->ChanRowCol(r,c);

        } else {
            Msg(EsafMsg::Warning) << "FE Channels: " << FrontEnd()->Channels() << MsgDispatch;
            Msg(EsafMsg::Warning) << "Pmt Channels: " << Geometry()->NumPads() << MsgDispatch;
            Msg(EsafMsg::Warning) << "Current channel row, col: " << r <<","<< c << MsgDispatch;
            Msg(EsafMsg::Warning) << "Probable mismatch between PMT and FE types"<< MsgDispatch;
            Msg(EsafMsg::Panic) << "Bad PMT-Front End Mapping" << MsgDispatch;
        }
    }
}
//
//______________________________________________________________________________
void Photomultiplier::AddTest(Double_t tm, Int_t ch) {
    // Add method for testing purpose only
    // private; only ElecTestDetTransManager can call it

    TRandom* rndm = EsafRandom::Get();

    // compute charge
    Double_t delta = fgGainSigma * rndm->Gaus();
    Double_t charge = fgGain + delta;
    if ( charge < 0. ) charge = 0.;
    charge *= EConst::ElectronCharge();   // coulomb

    // create the list if this does not exist
    if ( fPmtHits[ch] == NULL ) {
        fPmtHits[ch] = new vector<PmtSignal*>;
    }

    // create the PmtSignal and add it to the list
    PmtSignal* sig = new PmtSignal(tm,charge,fgWidth,GetUniqueId(ch),ch+1);
    fPmtHits[ch]->push_back( sig );

    SetState(kPmtFilling);

    // keep memory of the time interval of this event
    if ( sig->Time() > fEndTime )
        fEndTime = sig->Time();
    if ( sig->Time() < fStartTime )
        fStartTime = sig->Time();

    SetEmpty( false );
}

//______________________________________________________________________________
Bool_t Photomultiplier::Add(Photon& ph) {
    //
    // add a photon hit to be simulated
    //
    // check if this photon hits this pmt
    if (!Geometry()->IsInside( ph )) {
        TVector3 r = (ph.pos-Geometry()->Position()) ;
        Double_t x = r.Dot(Geometry()->GetX());
        Double_t y = r.Dot(Geometry()->GetY());
        Double_t z = r.Dot(Geometry()->GetZ());
        Printf("Wrong ph. DIFF=(%.3e, %.3e, %.3e)  PROJ=(%.3e, %.3e, %.3e)",
                r[0], r[1], r[2], x, y, z);
        return false;
    }

    // get the channel number
    Int_t ch = Geometry()->Pad(ph);

    // negative if Photon is lost for geometrical reasons (dead spaces)
    if ( ch < 0 ) return false;

    // saving pixel id on the photon
    ph.pixelUid = GetUniqueId(ch);

    // if required, do association mapping between original
    // theta phi in field of view and this channel
    //FIXME: EusoMapping::Get()->Associate(this,ch,ph);

    // if signal simulation is disabled, stop here
    if ( !(GetEusoElectronics()->GetSimulationStatus()) ) return true;

    TRandom* rndm = EsafRandom::Get();

    // handle quantum efficiency
    Double_t shot = rndm->Rndm();

    if ( shot > fgQuantum )
        return false;

    // compute charge
    Double_t delta = fgGainSigma * rndm->Gaus();
    Double_t charge = fgGain + delta;
    charge *= EConst::ElectronCharge();

    // create the list if this does not exist
    if ( fPmtHits[ch] == NULL ) {
        fPmtHits[ch] = new vector<PmtSignal*>;
    }

    // create the PmtSignal and add it to the list
    PmtSignal* sig = new PmtSignal(ph.time,charge,fgWidth,GetUniqueId(ch),ch+1);
    fPmtHits[ch]->push_back( sig );

    SetState(kPmtFilling);

    ph.madeSignal = true;

    // add this level of information to photon history in root file
    if ( EEvent::GetCurrent() ) {
        EDetectorPhotonAdder a(&ph,sig,false);
        EEvent::GetCurrent()->Fill(a);
    }

    // keep memory of the time interval of this event
    if ( sig->Time() > fEndTime )
        fEndTime = sig->Time();
    if ( sig->Time() < fStartTime )
        fStartTime = sig->Time();

    SetEmpty( false );

    return true;
}

//______________________________________________________________________________
void Photomultiplier::Simulate() {
    // Add PmtHits to the Front-End chip
    // the Photomultipliers job ends here. The Front End simulation
    // and MacroCell simulation is assumed to be handled by some one else.

    // do not do it twice or if empty
    if ( Status() != kPmtFilling)
        return;

    // loop on all channels
    for( Int_t ch=0; ch<Geometry()->NumPads(); ch++) {
        if ( fPmtHits[ch] )
            FrontEnd()->Add( fPmtHits[ch]  , fFeMap[ch] );
    }

    // set the status to the right value
    SetState(kPmtDone);
}