source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/simulation/detector/electronics/src/R8900M64Photomultiplier.cc @ 117

/*****************************************************************************
 * ESAF: Euso Simulation and Analysis Framework                              *
 *                                                                           *
 *  Id: R8900M64Photomultiplier                                               *
 *  Package: Electronics                                                     *
 *  Coordinator: Marco.Pallavicini                                           *
 *                                                                           *
 *****************************************************************************/

//_____________________________________________________________________________
//
//   R8900-M64 Weakly focused Photomultiplier
//   ========================================
//  
//   THIS IS TEMPORARY CLASS FOR ICRC2011 (M64 MAPMT)
//
//   ======================
//
//   fgUsePmtAdd [string]: if yes Add() method calls Photomultiplier::Add() and
//   uses average values for all pixels.
//

#include "R8900M64Photomultiplier.hh"
#include "EusoElectronics.hh"
#include "EusoDetector.hh"
#include "EsafRandom.hh"
#include "EConst.hh"
#include "EDetectorPhotonAdder.hh"
#include "NumbersFileParser.hh"
#include "Photon.hh"

using namespace sou;

ClassImp(R8900M64Photomultiplier)

R8900M64Photomultiplier::EStatus R8900M64Photomultiplier::fgStatus = kUndefined;

vector<Double_t> R8900M64Photomultiplier::fgColEff;
Interpolate* R8900M64Photomultiplier::fgQuantumEff           = NULL;
Interpolate* R8900M64Photomultiplier::fgAngularDependence    = NULL;
vector<Double_t>* R8900M64Photomultiplier::fgCrossTalkTable = NULL;

//______________________________________________________________________________
R8900M64Photomultiplier::R8900M64Photomultiplier( Int_t id, PmtGeometry* g ) :
    Photomultiplier(id, g ) {
    //
    // Constructor
    //
    // initialize static data members
    if ( fgStatus == kUndefined )  {

        if ( Conf()->GetBool("R8900M64Photomultiplier.fgUsePmtAdd") )
            fgStatus = kPmtAdd;
        else {
            fgStatus = kR8900Add;

            //--------------------  begin  Satoi & Mario 's changes  -----------------

            string cfgdir = Conf()->GetCfgDir()+"/"+Conf()->GetType()+"/"+Conf()->GetName()+"/";

            // load Collection Efficiency from file
            string s = cfgdir + Conf()->GetStr("R8900M64Photomultiplier.fCollectionEfficiency.FileName");
            NumbersFileParser pce(s,4);
            vector<Double_t> dummy=pce.GetCol(3);

            if(dummy.size()==0)
                Msg(EsafMsg::Panic) << "R8900M64Photomultipiler: "+s+"is empty!" << MsgDispatch;

            for(size_t i(0);i<dummy.size();i++)
                fgColEff.push_back(dummy[i]);
            Msg(EsafMsg::Debug) << fgColEff.size() << " C.E. data loading..." << MsgDispatch;

            // load Quantum Efficiency from file
            fgQuantumEff = new Interpolate(cfgdir+Conf()->GetStr("R8900M64Photomultiplier.fQuantumEfficiency.FileName"));
            fgQuantumEff->SetXUnit(nm);
            Msg(EsafMsg::Debug) << "Quantum Efficiency data loading..." << MsgDispatch;

            // load Angular Dependence from file
            fgAngularDependence = new Interpolate(cfgdir+Conf()->GetStr("R8900M64Photomultiplier.fAngularDependence.FileName"));
            fgAngularDependence->SetXUnit(deg);
            Msg(EsafMsg::Debug) << "Angular Dependece data loading..." << MsgDispatch;

            // load Cross Talk form file
            // Cross Talk File Format:
            // (ch.),x,y,(data1),(data2),...(data64)[CR] for each line
            fgCrossTalkTable= new vector<Double_t>[64];
            s = cfgdir + Conf()->GetStr("R8900M64Photomultiplier.fCrossTalk.FileName");
            NumbersFileParser *dummy2 = new NumbersFileParser(s,67);
            for(size_t j(0);j<64;j++){
                // Get row No."j"
                vector<Double_t> dummy2_line=dummy2->GetRow(j);
                // Get data from colum No.3 to No.67
                for(size_t k(3);k<dummy2_line.size();k++){
                    // Load fgCrossTalkTable
                    fgCrossTalkTable[j].push_back(dummy2_line[k]);
                }
            }

            Msg(EsafMsg::Debug) << "Cross Talk data loading..." << MsgDispatch;
        }
    }
    //--------------------  end  Satoi & Mario 's changes  -----------------

}

//______________________________________________________________________________
R8900M64Photomultiplier::~R8900M64Photomultiplier() {
    // Destructor
}


//______________________________________________________________________________
Bool_t R8900M64Photomultiplier::Add(Photon& ph) {
    //
    // add a photon hit to be simulated
    //
    // use Photomultiplier::Add, temporary disabled
    if ( fgStatus == kPmtAdd )
        return Photomultiplier::Add(ph);

    // new features after this line

    // 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());
        MsgForm(EsafMsg::Warning, "Wrong ph. DIFF=(%.3e, %.3e, %.3e)  PROJ=(%.3e, %.3e, %.3e)",
                r[0], r[1], r[2], x, y, z);
        return kFALSE;
    }

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

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

    // 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 ( !(GetEusoDetector()->GetEusoElectronics()->GetSimulationStatus()) ) return kTRUE;
    if ( !(GetEusoElectronics()->GetSimulationStatus()) ) return kTRUE;

    TRandom* rndm = EsafRandom::Get();

    //--------------------  beginning  Satoi & Mario 's changes  -----------------
    // fgQuantum is not the same as our fgQuantumEff.
    // fgQuantum    : constant.
    // fgQuantumEff : depend on wavelength

    // fgColEff depend on ch (pixel number).
    // QuantumEff depend on wavelength of photon.
    // AngularDependence depend on Theta of photon got by Detector.

    // Computing D.E.
    double ph_dir_theta = ph.dir.Theta();
    if (fGeometry){
        // calculate direction angle in relation with the geometry direction
        ph_dir_theta = (ph.dir.Angle(-fGeometry->Normal()));
    }
    else {
        Msg(EsafMsg::Panic) << "Geometry not defined." << MsgDispatch;
    }
    Double_t DetEff=fgColEff[ch]*fgQuantumEff->GetValue(ph.wl)*fgAngularDependence->GetValue(ph_dir_theta);

    // Shot for D.E.
    Double_t shot=rndm->Rndm(); // [0,1]
    // Check if photon generate electron or not
    if(shot > DetEff)
            return kFALSE;  // photon is lost

    // Shot2 for CrossTalk
    Double_t shot2=rndm->Rndm(); // [0,1]

    // Cross Talk
    //  Get Cross Talk Data
    vector<Double_t> *ChCrossTalk = &fgCrossTalkTable[ch];

    // apply CrossTalk Technique
    if (shot2 <= (*ChCrossTalk)[0] ) {
        ch=0;
    } else {
        for(size_t kch(1); kch < ChCrossTalk->size(); kch++){
            if( (shot2 >= (*ChCrossTalk)[kch-1] ) && ( shot2 <= (*ChCrossTalk)[kch] ) ){
                ch=kch;
                break;
            }
        }
    }

    // end of judge
    // Msg(EsafMsg::Debug) << "photon hit pixel no. " << ch <<  MsgDispatch;

    //-------------------------------- end Satoi & Mario 's changes --------------

    // 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 = kTRUE;

    // add this level of information to photon history in root file
    if ( EEvent::GetCurrent() ) {
        EDetectorPhotonAdder a(&ph,sig,kFALSE);
        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( kFALSE );

    return kTRUE;
}