// $Id: R8900M36Photomultiplier.cc,v 1.18 2005/04/26 15:38:53 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"
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;
}