Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

R8900M36Photomultiplier - source file

// $Id: R8900M36Photomultiplier.cc,v 1.19 2005/10/02 14:17:55 thea Exp $

// A.Thea created Mar, 26 2003

//

// S.Ogawa & M.Bertaina

// 2004-11-26

// LastUpdate: 2005-01-20


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

//_____________________________________________________________________________

//  

//   R8900-M36 Weakly focused Photomultiplier

//   ========================================

//

//   This class implements the R8900 36 pixels weakly focused photomultiplier.

//   In weakly focused photomultipliers border pixel behaviour can be sensibly

//   differet from the centra one. A static map, with the average pixel

//   detector efficiency has been introduced to simulate the border effect.

//   Efficiency Wavelenght dependency has been also introduced. 

//   

//   Config file parameters

//   ======================

//

//   fgUsePmtAdd [string]: if yes Add() method calls Photomultiplier::Add() and

//   uses average values for all pixels.

//


#include "R8900M36Photomultiplier.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(R8900M36Photomultiplier)

R8900M36Photomultiplier::EStatus R8900M36Photomultiplier::fgStatus = kUndefined;

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

//______________________________________________________________________________

 R8900M36Photomultiplier::R8900M36Photomultiplier( Int_t id, PmtGeometry* g ) : 
    Photomultiplier(id, g ) {
    //

    // Constructor

    // 


    // initialize static data members

    if ( fgStatus == kUndefined )  {

        if ( Conf()->GetBool("R8900M36Photomultiplier.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("R8900M36Photomultiplier.fCollectionEfficiency.FileName");
            NumbersFileParser pce(s,4);
            vector<Double_t> dummy=pce.GetCol(3);

            if(dummy.size()==0)
                Msg(EsafMsg::Panic) << "R8900M36Photomultipiler: "+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. Datas load..." << MsgDispatch;

            // load Quantum Efficiency from file

            fgQuantumEff = new Interpolate(cfgdir+Conf()->GetStr("R8900M36Photomultiplier.fQuantumEfficiency.FileName"));
            fgQuantumEff->SetXUnit(nm);
            Msg(EsafMsg::Debug) << "Quantum Efficiency Datas load..." << MsgDispatch;

            // load Angular Dependence from file

            fgAngularDependence = new Interpolate(cfgdir+Conf()->GetStr("R8900M36Photomultiplier.fAngularDependence.FileName"));
            fgAngularDependence->SetXUnit(deg);
            Msg(EsafMsg::Debug) << "Angular Dependece Datas load..." << MsgDispatch;

            // load Cross Talk form file

            // Cross Talk File Format:

            // (ch.),x,y,(data1),(data2),...(data36)[CR] for each line

            fgCrossTalkTable= new vector<Double_t>[36];
            s = cfgdir + Conf()->GetStr("R8900M36Photomultiplier.fCrossTalk.FileName");
            NumbersFileParser *dummy2 = new NumbersFileParser(s,39);
            for(size_t j(0);j<36;j++){
                // Get row No."j"

                vector<Double_t> dummy2_line=dummy2->GetRow(j);
                // Get data from colum No.3 to No.39 

                for(size_t k(3);k<dummy2_line.size();k++){
                    // Load fgCrossTalkTable

                    fgCrossTalkTable[j].push_back(dummy2_line[k]);
                }
            }

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


}

//______________________________________________________________________________

 R8900M36Photomultiplier::~R8900M36Photomultiplier() {
    // Destructor

}


//______________________________________________________________________________

 Bool_t R8900M36Photomultiplier::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;

    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_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; 
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005