#include <EsPmtEffectPulseTool.h>
Inheritance diagram for EsPmtEffectPulseTool:
Public Types | |
SUCCESS | |
NO_INTERFACE | |
VERSMISMATCH | |
LAST_ERROR | |
SUCCESS | |
NO_INTERFACE | |
VERSMISMATCH | |
LAST_ERROR | |
enum | Status |
Public Member Functions | |
EsPmtEffectPulseTool (const std::string &type, const std::string &name, const IInterface *parent) | |
virtual | ~EsPmtEffectPulseTool () |
virtual StatusCode | generatePulses (DayaBay::SimHitCollection *, DayaBay::ElecPulseCollection *) |
This is the extension. Sub classes must provide it. | |
virtual StatusCode | initialize () |
virtual StatusCode | finalize () |
INTupleSvc * | ntupleSvc () const |
INTupleSvc * | evtColSvc () const |
IDataProviderSvc * | detSvc () const |
IDataProviderSvc * | evtSvc () const |
IIncidentSvc * | incSvc () const |
IChronoStatSvc * | chronoSvc () const |
IHistogramSvc * | histoSvc () const |
IAlgContextSvc * | contextSvc () const |
DataObject * | put (IDataProviderSvc *svc, DataObject *object, const std::string &address, const bool useRootInTES=true) const |
DataObject * | put (DataObject *object, const std::string &address, const bool useRootInTES=true) const |
Gaudi::Utils::GetData< TYPE >::return_type | get (IDataProviderSvc *svc, const std::string &location, const bool useRootInTES=true) const |
Gaudi::Utils::GetData< TYPE >::return_type | get (const std::string &location, const bool useRootInTES=true) const |
TYPE * | getDet (IDataProviderSvc *svc, const std::string &location) const |
TYPE * | getDet (const std::string &location) const |
bool | exist (IDataProviderSvc *svc, const std::string &location, const bool useRootInTES=true) const |
bool | exist (const std::string &location, const bool useRootInTES=true) const |
bool | existDet (IDataProviderSvc *svc, const std::string &location) const |
bool | existDet (const std::string &location) const |
TYPE * | getOrCreate (IDataProviderSvc *svc, const std::string &location, const bool useRootInTES=true) const |
TYPE * | getOrCreate (const std::string &location, const bool useRootInTES=true) const |
TOOL * | tool (const std::string &type, const std::string &name, const IInterface *parent=0, bool create=true) const |
TOOL * | tool (const std::string &type, const IInterface *parent=0, bool create=true) const |
SERVICE * | svc (const std::string &name, const bool create=true) const |
IUpdateManagerSvc * | updMgrSvc () const |
IDataProviderSvc * | fastContainersSvc () const |
StatusCode | Error (const std::string &msg, const StatusCode st=StatusCode::FAILURE, const size_t mx=10) const |
StatusCode | Warning (const std::string &msg, const StatusCode st=StatusCode::FAILURE, const size_t mx=10) const |
StatusCode | Print (const std::string &msg, const StatusCode st=StatusCode::SUCCESS, const MSG::Level lev=MSG::INFO) const |
StatusCode | Assert (const bool ok, const std::string &message="", const StatusCode sc=StatusCode(StatusCode::FAILURE, true)) const |
StatusCode | Assert (const bool ok, const char *message, const StatusCode sc=StatusCode(StatusCode::FAILURE, true)) const |
StatusCode | Exception (const std::string &msg, const GaudiException &exc, const StatusCode sc=StatusCode(StatusCode::FAILURE, true)) const |
StatusCode | Exception (const std::string &msg, const std::exception &exc, const StatusCode sc=StatusCode(StatusCode::FAILURE, true)) const |
StatusCode | Exception (const std::string &msg="no message", const StatusCode sc=StatusCode(StatusCode::FAILURE, true)) const |
MsgStream & | msgStream (const MSG::Level level) const |
MsgStream & | always () const |
MsgStream & | fatal () const |
MsgStream & | err () const |
MsgStream & | error () const |
MsgStream & | warning () const |
MsgStream & | info () const |
MsgStream & | debug () const |
MsgStream & | verbose () const |
MsgStream & | msg () const |
const Statistics & | counters () const |
StatEntity & | counter (const std::string &tag) const |
MSG::Level | msgLevel () const |
bool | msgLevel (const MSG::Level level) const |
void | resetMsgStream () const |
bool | typePrint () const |
bool | propsPrint () const |
bool | statPrint () const |
bool | errorsPrint () const |
long | printStat (const MSG::Level level=MSG::ALWAYS) const |
long | printErrors (const MSG::Level level=MSG::ALWAYS) const |
long | printProps (const MSG::Level level=MSG::ALWAYS) const |
void | registerCondition (const std::string &condition, StatusCode(CallerClass::*mf)()=NULL) |
void | registerCondition (const std::string &condition, CondType *&condPtrDest, StatusCode(CallerClass::*mf)()=NULL) |
void | registerCondition (char *condition, StatusCode(CallerClass::*mf)()=NULL) |
void | registerCondition (TargetClass *condition, StatusCode(CallerClass::*mf)()=NULL) |
StatusCode | runUpdate () |
TransientFastContainer< T > * | getFastContainer (const std::string &location, typename TransientFastContainer< T >::size_type initial=0) |
StatusCode | release (const IInterface *interface) const |
virtual unsigned long | release () |
const std::string & | context () const |
const std::string & | rootInTES () const |
double | globalTimeOffset () const |
virtual StatusCode | queryInterface (const InterfaceID &riid, void **ppvUnknown) |
virtual unsigned long | addRef () |
virtual const std::string & | name () const |
virtual const std::string & | type () const |
virtual const IInterface * | parent () const |
virtual StatusCode | configure () |
virtual StatusCode | start () |
virtual StatusCode | stop () |
virtual StatusCode | terminate () |
virtual StatusCode | reinitialize () |
virtual StatusCode | restart () |
virtual Gaudi::StateMachine::State | FSMState () const |
virtual Gaudi::StateMachine::State | targetFSMState () const |
virtual StatusCode | sysInitialize () |
virtual StatusCode | sysStart () |
virtual StatusCode | sysStop () |
virtual StatusCode | sysFinalize () |
virtual StatusCode | sysReinitialize () |
virtual StatusCode | sysRestart () |
virtual StatusCode | setProperty (const Property &p) |
virtual StatusCode | setProperty (const std::string &s) |
virtual StatusCode | setProperty (const std::string &n, const std::string &v) |
StatusCode | setProperty (const std::string &name, const TYPE &value) |
virtual StatusCode | getProperty (Property *p) const |
virtual const Property & | getProperty (const std::string &name) const |
virtual StatusCode | getProperty (const std::string &n, std::string &v) const |
virtual const std::vector< Property * > & | getProperties () const |
PropertyMgr * | getPropertyMgr () |
ISvcLocator * | serviceLocator () const |
ISvcLocator * | svcLoc () const |
IMessageSvc * | msgSvc () const |
IToolSvc * | toolSvc () const |
StatusCode | setProperties () |
StatusCode | service (const std::string &name, T *&svc, bool createIf=true) const |
StatusCode | service (const std::string &type, const std::string &name, T *&svc) const |
void | declInterface (const InterfaceID &, void *) |
Property * | declareProperty (const std::string &name, T &property, const std::string &doc="none") const |
Property * | declareRemoteProperty (const std::string &name, IProperty *rsvc, const std::string &rname="") const |
IAuditorSvc * | auditorSvc () const |
IMonitorSvc * | monitorSvc () const |
void | declareInfo (const std::string &name, const T &var, const std::string &desc) const |
void | declareInfo (const std::string &name, const std::string &format, const void *var, int size, const std::string &desc) const |
virtual const std::string & | type () const =0 |
virtual const IInterface * | parent () const =0 |
virtual StatusCode | configure ()=0 |
virtual StatusCode | start ()=0 |
virtual StatusCode | stop ()=0 |
virtual StatusCode | terminate ()=0 |
virtual StatusCode | reinitialize ()=0 |
virtual StatusCode | restart ()=0 |
virtual Gaudi::StateMachine::State | FSMState () const =0 |
virtual StatusCode | sysInitialize ()=0 |
virtual StatusCode | sysStart ()=0 |
virtual StatusCode | sysStop ()=0 |
virtual StatusCode | sysFinalize ()=0 |
virtual StatusCode | sysReinitialize ()=0 |
virtual StatusCode | sysRestart ()=0 |
virtual unsigned long | refCount () const =0 |
virtual const std::string & | name () const =0 |
virtual StatusCode | queryInterface (const InterfaceID &riid, void **ppvInterface)=0 |
virtual unsigned long | addRef ()=0 |
virtual unsigned long | release ()=0 |
Static Public Member Functions | |
static const InterfaceID & | interfaceID () |
static const InterfaceID & | interfaceID () |
static const InterfaceID & | interfaceID () |
static const InterfaceID & | interfaceID () |
Retrieve interface ID. | |
Public Attributes | |
SUCCESS | |
NO_INTERFACE | |
VERSMISMATCH | |
LAST_ERROR | |
Protected Types | |
typedef std::map< std::string, StatEntity > | Statistics |
typedef std::map< std::string, unsigned int > | Counter |
typedef std::vector< IAlgTool * > | AlgTools |
typedef std::pair< IInterface *, std::string > | ServiceEntry |
typedef std::vector< ServiceEntry > | Services |
Protected Member Functions | |
StatusCode | releaseTool (const IAlgTool *tool) const |
StatusCode | releaseSvc (const IInterface *svc) const |
int | outputLevel () const |
virtual unsigned long | refCount () const |
IntegerProperty & | outputLevelProperty () |
void | initOutputLevel (Property &prop) |
Static Protected Attributes | |
static const bool | IgnoreRootInTES |
static const bool | UseRootInTES |
Private Member Functions | |
DayaBay::ElecPulse * | makePrePulse (DayaBay::ElecPulse *simpulse, const DayaBay::PmtSimData &pmtData) |
DayaBay::ElecPulse * | makeAfterPulse (DayaBay::ElecPulse *simpulse, const DayaBay::PmtSimData &pmtData) |
DayaBay::ElecPulse * | makeDarkPulse (const DayaBay::ElecChannelId &channelId, const DayaBay::PmtSimData &pmtData) |
float | PulseAmp (float weight, float gain, float gainFwhm) |
double | ConvertPdfRand (double rand, std::vector< double > pdf, std::vector< double > edges) |
double | ConvertPdfRand01 (double rand, std::vector< double > pdf, std::vector< double > edges) |
int | PoissonRand (double mean) |
double | NumAfterPulse (const int numPmtHit) |
void | getAfterPulseAmpPdf (std::vector< double > &pdf) |
void | getAfterPulseAmpEdges (std::vector< double > &edges) |
Private Attributes | |
std::string | m_cableSvcName |
std::string | m_simDataSvcName |
std::vector< double > | m_prePulsePdf |
std::vector< double > | m_prePulseEdges |
std::vector< double > | m_afterPulsePdf |
std::vector< double > | m_afterPulseEdges |
std::vector< double > | m_afterPulseTime |
bool | m_enablePrePulse |
bool | m_enableAfterPulse |
std::string | m_afterPulseAmpMode |
std::string | m_darkPulseAmpMode |
std::vector< double > | m_afterPulseAmpPdf |
std::vector< double > | m_afterPulseAmpEdges |
bool | m_enableNonlinearAfterpulse |
double | m_linearAfterpulseThreshold |
bool | m_enableVeryLongTimeSuppression |
double | m_veryLongTime |
double | m_expWeight |
double | m_speExpDecay |
double | m_speCutoff |
ICableSvc * | m_cableSvc |
ISimDataSvc * | m_simDataSvc |
Rndm::Numbers | m_uni |
Rndm::Numbers | m_gauss |
Rndm::Numbers | m_exp |
Rndm::Numbers | m_randPrePulseTime |
Rndm::Numbers | m_randAfterPulseTime |
TimeStamp | m_simTimeEarliest |
TimeStamp | m_simTimeLatest |
double | m_timeInterval |
Definition at line 33 of file EsPmtEffectPulseTool.h.
EsPmtEffectPulseTool::EsPmtEffectPulseTool | ( | const std::string & | type, | |
const std::string & | name, | |||
const IInterface * | parent | |||
) |
Definition at line 45 of file EsPmtEffectPulseTool.cc.
00048 : GaudiTool(type,name,parent) 00049 , m_cableSvc(0) 00050 , m_simDataSvc(0) 00051 , m_simTimeEarliest(0) 00052 , m_simTimeLatest(0) 00053 { 00054 // Initialization 00055 m_prePulsePdf.clear(); m_prePulsePdf.push_back(1); 00056 m_prePulseEdges.clear(); m_prePulseEdges.push_back(-50.0*ns); m_prePulseEdges.push_back(0.0*ns); 00057 // m_afterPulsePdf.clear(); m_afterPulsePdf.push_back(1); 00058 // m_afterPulseEdges.clear(); m_afterPulseEdges.push_back(0.0*ns); m_afterPulseEdges.push_back(10.0e3*ns); 00059 m_afterPulsePdf.clear(); 00060 m_afterPulseEdges.clear(); 00061 00062 for(int ii=0; ii < NUM_BINS_DIST; ii++) { 00063 m_afterPulsePdf.push_back(afterPusleTimingDist[ii]); 00064 m_afterPulseEdges.push_back(ii*100.+500.); 00065 } 00066 00067 debug() << "time(ns) Afterpulse PDF (integrated) " <<endreq; 00068 for(int numLine=0; numLine < NUM_BINS_DIST; numLine++) 00069 debug() << m_afterPulseEdges[numLine] << ", " << m_afterPulsePdf[numLine] << endreq; 00070 // std::cout << m_afterPulseEdges[numLine] << ", " << m_afterPulsePdf[numLine] << std::endl; 00071 00072 declareInterface< IEsPulseTool >(this) ; 00073 declareProperty("CableSvcName",m_cableSvcName="StaticCableSvc", 00074 "Name of service to map between detector, hardware, and electronic IDs"); 00075 declareProperty("SimDataSvcName",m_simDataSvcName="StaticSimDataSvc", 00076 "Name of service to provide pmt properties for simulation"); 00077 declareProperty("EnablePrePulse",m_enablePrePulse=true,"Switch on/off prepulse simulation"); 00078 declareProperty("EnableAfterPulse",m_enableAfterPulse=true,"Switch on/off afterpulse simulation"); 00079 declareProperty("PrePulsePdf",m_prePulsePdf,"User-defined PDF for timing distribution of pre-pulse"); 00080 declareProperty("PrePulsePdfEdges",m_prePulseEdges,"Bin edges of PrePulsePdf"); 00081 declareProperty("AfterPulseAmpMode",m_afterPulseAmpMode="SinglePE","Mode for afterpulse amplitude distribution"); 00082 declareProperty("DarkPulseAmpMode",m_darkPulseAmpMode="SinglePE","Mode for dark pulse amplitude distribution"); 00083 declareProperty("AfterPulsePdf",m_afterPulsePdf,"User-defined PDF for timing distribution of after-pulse"); 00084 declareProperty("AfterPulsePdfEdges",m_afterPulseEdges,"Bin edges of AfterPulsePdf"); 00085 declareProperty("AfterPulseTimeInterval", m_timeInterval=20*ns, 00086 "PMT hit counting window size"); 00087 declareProperty("EnableNonlinearAfterpulse", m_enableNonlinearAfterpulse=true, 00088 "Enable nonlinear suppression of afterpulsing."); 00089 declareProperty("LinearAfterpulseThreshold", m_linearAfterpulseThreshold=400, 00090 "Upper limit of linear afterpulsing in number of PE."); 00091 00092 00093 declareProperty("EnableVeryLongTimeSuppression",m_enableVeryLongTimeSuppression=false, 00094 "Enable suppression of hit times in the far future"); 00095 declareProperty("VeryLongTime",m_veryLongTime=3e7*s, 00096 "Definition of very long time for hit time suppression"); 00097 declareProperty("ExpWeight", m_expWeight=0,"Weight of the exponential contribution to the spe response function"); 00098 declareProperty("SpeExpDecay", m_speExpDecay=0.5,"Decay time of the exponential contribution to the spe response function"); 00099 }
EsPmtEffectPulseTool::~EsPmtEffectPulseTool | ( | ) | [virtual] |
StatusCode EsPmtEffectPulseTool::generatePulses | ( | DayaBay::SimHitCollection * | , | |
DayaBay::ElecPulseCollection * | ||||
) | [virtual] |
This is the extension. Sub classes must provide it.
Implements IEsPulseTool.
Definition at line 190 of file EsPmtEffectPulseTool.cc.
00192 { 00193 debug() << "running generatePulses()" << endreq; 00194 00195 // Define simulation time window 00196 m_simTimeEarliest = pulses->header()->header()->earliest(); 00197 m_simTimeLatest = pulses->header()->header()->latest(); 00198 TimeStamp headerTime = pulses->header()->header()->timeStamp(); 00199 00200 // Maintain a list of primary pulses for each PMT 00201 std::map<DayaBay::DetectorSensor, 00202 std::vector<DayaBay::ElecPulse*> > pulsesByPmt; 00203 std::map<DayaBay::DetectorSensor, 00204 std::vector<DayaBay::ElecPulse*> >::iterator pulsePmtIter, 00205 pulsePmtEnd; 00206 00207 const DayaBay::SimHitCollection::hit_container& htcon = hits->collection(); 00208 DayaBay::SimHitCollection::hit_container::const_iterator hcIter, 00209 hcDone = htcon.end(); 00210 DayaBay::ElecPulseCollection::PulseContainer& pulsecon = pulses->pulses(); 00211 00212 // Context, ServiceMode for this data 00213 Context context(pulses->detector().site(), SimFlag::kMC, 00214 pulses->header()->header()->timeStamp(), 00215 pulses->detector().detectorId()); 00216 int task = 0; 00217 ServiceMode svcMode(context, task); 00218 00219 debug() << "Processing " << htcon.size() << " simulated hits from detector " 00220 << hits->detector() << endreq; 00221 00222 00223 for (hcIter=htcon.begin(); hcIter != hcDone; ++hcIter) { 00224 00225 DayaBay::DetectorSensor sensDetId((*hcIter)->sensDetId()); 00226 00227 if (sensDetId.isAD()) { 00228 DayaBay::AdPmtSensor pmtid(sensDetId.fullPackedData()); 00229 if (pmtid.bogus()) { 00230 warning() << "Got bogus AD PMT: " << pmtid << endreq; 00231 } 00232 } 00233 else { 00234 DayaBay::PoolPmtSensor pmtid(sensDetId.fullPackedData()); 00235 if (pmtid.bogus()) { 00236 warning() << "Got bogus Pool PMT: " << pmtid << endreq; 00237 } 00238 } 00239 00240 //temporarily ignore very long time hit to avoid program break 00241 //if((*hcIter)->hitTime() > 1e7) { 00242 if((*hcIter)->hitTime() > m_veryLongTime ) { 00243 warning() << "Skipping very long hit time of " << (*hcIter)->hitTime() 00244 << " in PMT " << sensDetId << endreq; 00245 continue; 00246 } 00247 00248 verbose() << "Requesting pmtData for " << sensDetId << endreq; 00249 const DayaBay::PmtSimData* pmtData = 00250 m_simDataSvc->pmtSimData(sensDetId, svcMode); 00251 if(!pmtData){ 00252 warning() << "No Simulation input properties for PMT: " << sensDetId 00253 << " Ignoring..." << endreq; 00254 continue; 00255 } 00256 // Ignore SimHit due to efficiency? 00257 if (m_uni() > pmtData->m_efficiency) continue; 00258 00259 00260 const DayaBay::ElecChannelId& elecChannelId = m_cableSvc->elecChannelId( 00261 sensDetId, 00262 svcMode); 00263 if(elecChannelId.fullPackedData() == 0){ 00264 warning() << "PMT: " << sensDetId 00265 << " is not connected in cable map. Ignoring..." << endreq; 00266 continue; 00267 } 00268 DayaBay::ElecPulse* simpulse = new DayaBay::ElecPmtPulse(); 00269 00270 // Primary vertex time (absolute, in seconds) 00271 TimeStamp vertexTime = (*hcIter)->hc()->header()->header()->timeStamp(); 00272 // Need to subtract ElecHeader time (absolute, in seconds) to get 00273 // pulse time wrt ElecHeader 00274 // Set pulse time, include transit time (nanoseconds) 00275 double transitTime = m_gauss() * pmtData->m_timeSpread 00276 + pmtData->m_timeOffset; 00277 double headerOffset = double(vertexTime - headerTime)*Gaudi::Units::second; 00278 double pulseHitTime = headerOffset + (*hcIter)->hitTime() + transitTime; 00279 00280 simpulse->setTime(pulseHitTime); 00281 00282 simpulse->setAmplitude(PulseAmp((*hcIter)->weight(),pmtData->m_gain, 00283 pmtData->m_sigmaGain)); // Include SPE effects, relative gain 00284 simpulse->setAncestor((*hcIter)); 00285 simpulse->setType(PulseType::kPmtHit); 00286 simpulse->setChannelId(elecChannelId); 00287 00288 pulsecon.push_back(simpulse); 00289 00290 // PMT pulse counting 00291 if(m_enableAfterPulse){ 00292 if( m_enableNonlinearAfterpulse ){ 00293 // Add pulse to the list for this PMT 00294 pulsesByPmt[sensDetId].push_back(simpulse); 00295 }else{ 00296 // Immediatly add linear afterpulsing 00297 if (m_uni() < pmtData->m_afterPulseProb) 00298 pulsecon.push_back(makeAfterPulse(simpulse, *pmtData)); 00299 } 00300 } 00301 // Include pre-pulse? 00302 if (m_uni() < pmtData->m_prePulseProb && m_enablePrePulse){ 00303 pulsecon.push_back(makePrePulse(simpulse, *pmtData)); 00304 } 00305 } 00306 00307 // Add afterpulses in nonlinear mode 00308 if(m_enableAfterPulse){ 00309 if( m_enableNonlinearAfterpulse){ 00310 pulsePmtEnd = pulsesByPmt.end(); 00311 for(pulsePmtIter = pulsesByPmt.begin();pulsePmtIter != pulsePmtEnd;pulsePmtIter++){ 00312 const DayaBay::DetectorSensor& sensDetId = pulsePmtIter->first; 00313 std::vector<DayaBay::ElecPulse*>& pulseList = pulsePmtIter->second; 00314 // Get PMT properties 00315 const DayaBay::PmtSimData* pmtData = 00316 m_simDataSvc->pmtSimData(sensDetId, svcMode); 00317 if(!pmtData){ 00318 error() << "Nonlinear: No Simulation input properties for PMT: " 00319 << sensDetId << endreq; 00320 return StatusCode::FAILURE; 00321 } 00322 00323 // Prepare a handle for memory buffer 00324 int* pulseCount = 0; 00325 00326 // Process hits to make afterpulses 00327 std::vector<DayaBay::ElecPulse*>::iterator pulseIter, 00328 pulseEnd = pulseList.end(); 00329 if(pulseList.size() > m_linearAfterpulseThreshold){ 00330 // Add Nonlinear afterpulsing 00331 // Step 1: assemble pulse counts in a specified time window 00332 int numTimeSlots = double(m_simTimeLatest - headerTime) 00333 *Gaudi::Units::second/m_timeInterval + 1; 00334 // allocate memory only if needed 00335 if(!pulseCount) pulseCount = new int[numTimeSlots]; 00336 // Clear memory buffer 00337 for(int i=0; i<numTimeSlots; i++) pulseCount[i] = 0; 00338 // Loop once to create pulse counts 00339 for(pulseIter=pulseList.begin(); pulseIter != pulseEnd; pulseIter++){ 00340 int hitTimeSlot = int((*pulseIter)->time()/m_timeInterval); 00341 int pulseNo = int((*pulseIter)->amplitude()+0.5); 00342 pulseCount[hitTimeSlot] += pulseNo; 00343 } 00344 // Loop again to create afterpulses 00345 for(pulseIter=pulseList.begin(); pulseIter != pulseEnd; pulseIter++){ 00346 int hitTimeSlot = int((*pulseIter)->time()/m_timeInterval); 00347 double nonlinearProb = pmtData->m_afterPulseProb/0.0164 00348 *NumAfterPulse(pulseCount[hitTimeSlot]) / pulseCount[hitTimeSlot]; 00349 if (m_uni() < nonlinearProb) 00350 pulsecon.push_back(makeAfterPulse(*pulseIter, *pmtData)); 00351 } 00352 }else{ 00353 // Add Linear afterpulsing 00354 for(pulseIter=pulseList.begin(); pulseIter != pulseEnd; pulseIter++){ 00355 if (m_uni() < pmtData->m_afterPulseProb) 00356 pulsecon.push_back(makeAfterPulse(*pulseIter, *pmtData)); 00357 } 00358 } 00359 // Free memory buffer if it was used 00360 if(pulseCount) delete [] pulseCount; 00361 } 00362 } 00363 } 00364 // Process dark pulses 00365 // Get list of all detector sensors 00366 debug() << "Adding dark hits" << endreq; 00367 00368 const std::vector<DayaBay::DetectorSensor>& detSensors = 00369 m_cableSvc->sensors( svcMode ); 00370 00371 for (size_t idet = 0; idet < detSensors.size(); ++idet) { 00372 DayaBay::DetectorSensor detSens = detSensors[idet]; 00373 if (detSens.isAD() ) { 00374 DayaBay::AdPmtSensor pmtid(detSens.fullPackedData()); 00375 if (pmtid.bogus()) { 00376 warning() << "bogus AD PMT id: " << pmtid << endreq; 00377 } 00378 } 00379 else { 00380 DayaBay::PoolPmtSensor pmtid(detSens.fullPackedData()); 00381 if (pmtid.bogus()) { 00382 warning() << "bogus Pool PMT id: " << pmtid << endreq; 00383 } 00384 } 00385 } 00386 00387 00388 for (unsigned int idet=0; idet < detSensors.size(); ++idet) { 00389 // Skip this detector sensor if not in detector of ElecPulseCollection/SimHitCollection 00390 DayaBay::DetectorSensor detSens = detSensors[idet]; 00391 00392 if (detSens.detName() != pulses->detector().detName()) continue; 00393 00394 const DayaBay::ElecChannelId& channelId = m_cableSvc->elecChannelId( 00395 detSens, 00396 svcMode); 00397 00398 const DayaBay::PmtSimData* pmtData = 00399 m_simDataSvc->pmtSimData(detSens, svcMode); 00400 if(!pmtData){ 00401 error() << "No Simulation input properties for PMT #" << idet 00402 << ": id = " << detSens 00403 << " svcMode = " << svcMode 00404 << endreq; 00405 return StatusCode::FAILURE; 00406 } 00407 // Generate Poisson-distributed number around mean number of dark hits in simulation time window 00408 int Ndark = PoissonRand(pmtData->m_darkRate 00409 *double(m_simTimeLatest-m_simTimeEarliest) 00410 *Gaudi::Units::second); 00411 verbose() << "Adding " << Ndark << " dark hits on pmt " 00412 << detSens << " with dark rate " << pmtData->m_darkRate 00413 << " and dt " 00414 << double(m_simTimeLatest-m_simTimeEarliest)*Gaudi::Units::second 00415 << endreq; 00416 for (int dummy = 0; dummy < Ndark; ++dummy) { 00417 pulsecon.push_back(makeDarkPulse(channelId, *pmtData)); 00418 } 00419 } 00420 00421 debug() << "Adding " << pulsecon.size() << " pulses to detector " 00422 << pulses->detector() << endreq; 00423 return StatusCode::SUCCESS; 00424 }
StatusCode EsPmtEffectPulseTool::initialize | ( | ) | [virtual] |
Reimplemented from GaudiTool.
Definition at line 103 of file EsPmtEffectPulseTool.cc.
00104 { 00105 // Get Cable Service 00106 m_cableSvc = svc<ICableSvc>(m_cableSvcName,true); 00107 00108 // Get PMT simulation input data service 00109 m_simDataSvc = svc<ISimDataSvc>(m_simDataSvcName,true); 00110 00111 // Random number service 00112 IRndmGenSvc *rgs = 0; 00113 if (service("RndmGenSvc",rgs,true).isFailure()) { 00114 fatal() << "Failed to get random service" << endreq; 00115 return StatusCode::FAILURE; 00116 } 00117 00118 StatusCode sc; 00119 sc = m_uni.initialize(rgs, Rndm::Flat(0,1)); 00120 if (sc.isFailure()) { 00121 fatal() << "Failed to initialize uniform random numbers" << endreq; 00122 return StatusCode::FAILURE; 00123 } 00124 sc = m_gauss.initialize(rgs, Rndm::Gauss(0,1)); 00125 if (sc.isFailure()) { 00126 fatal() << "Failed to initialize gaussian random numbers" << endreq; 00127 return StatusCode::FAILURE; 00128 } 00129 sc = m_exp.initialize(rgs, Rndm::Exponential(m_speExpDecay)); 00130 if (sc.isFailure()) { 00131 fatal() << "Failed to initialize exponential random numbers" << endreq; 00132 return StatusCode::FAILURE; 00133 } 00134 if(m_enablePrePulse){ 00135 info() << "Prepulse simulation enabled" << endreq; 00136 sc = m_randPrePulseTime.initialize(rgs, Rndm::DefinedPdf(m_prePulsePdf,0)); 00137 if (sc.isFailure()) { 00138 fatal() << "Failed to initialize random numbers for timing distribution" << endreq; 00139 return StatusCode::FAILURE; 00140 } 00141 } 00142 else info() << "Prepulse simulation disabled" << endreq; 00143 00144 if(m_enableAfterPulse){ 00145 info() << "Afterpulse simulation enabled" << endreq; 00146 sc = m_randAfterPulseTime.initialize(rgs, Rndm::DefinedPdf(m_afterPulsePdf,0)); 00147 if (sc.isFailure()) { 00148 fatal() << "Failed to initialize random numbers for timing distribution" << endreq; 00149 return StatusCode::FAILURE; 00150 } 00151 } 00152 else info() << "Afterpulse simulation disabled" << endreq; 00153 00154 // Check user input for user-defined PDFs 00155 if(m_enablePrePulse){ 00156 if (m_prePulsePdf.size()+1 != m_prePulseEdges.size()) { 00157 // || m_afterPulsePdf.size()+1 != m_afterPulseEdges.size() ) { 00158 fatal() << "There should be 1 more number of bin edges than bins when defining pdf" << endreq; 00159 return StatusCode::FAILURE; 00160 } 00161 } 00162 00163 // DWL: Check user input for user-defined PMT hit counting time interval 00164 if (m_timeInterval <= 0) { 00165 fatal() << "PMT hit couting time interval should be great than zero !" << endreq; 00166 return StatusCode::FAILURE; 00167 } 00168 if(m_enableAfterPulse){ 00169 if ((m_afterPulseAmpMode != "SinglePE") && (m_afterPulseAmpMode != "PDF")) { 00170 fatal() << "Not a recognized afterpulse amplitude distribution mode" << endreq; 00171 return StatusCode::FAILURE; 00172 } 00173 if (m_afterPulseAmpMode == "SinglePE") debug() << "Assume single p.e. amplitude for afterpulses" << endreq; 00174 if (m_afterPulseAmpMode == "PDF" || m_darkPulseAmpMode == "PDF") { 00175 debug() << "Set up amplitude PDF for afterpulses" << endreq; 00176 m_afterPulseAmpPdf.clear(); 00177 m_afterPulseAmpEdges.clear(); 00178 getAfterPulseAmpPdf(m_afterPulseAmpPdf); 00179 getAfterPulseAmpEdges(m_afterPulseAmpEdges); 00180 } 00181 } 00182 return StatusCode::SUCCESS; 00183 }
StatusCode EsPmtEffectPulseTool::finalize | ( | ) | [virtual] |
Reimplemented from GaudiTool.
Definition at line 185 of file EsPmtEffectPulseTool.cc.
00186 { 00187 return StatusCode::SUCCESS; 00188 }
DayaBay::ElecPulse * EsPmtEffectPulseTool::makePrePulse | ( | DayaBay::ElecPulse * | simpulse, | |
const DayaBay::PmtSimData & | pmtData | |||
) | [private] |
Definition at line 426 of file EsPmtEffectPulseTool.cc.
00428 { 00429 DayaBay::ElecPulse* prepulse = new DayaBay::ElecPmtPulse(); 00430 00431 // Time offset from main pulse based on time PDF of pre-pulses 00432 double current_rand = m_randPrePulseTime(); 00433 double prePulseTime = ConvertPdfRand(current_rand,m_prePulsePdf,m_prePulseEdges); 00434 00435 prepulse->setTime(simpulse->time() + prePulseTime); 00436 prepulse->setAmplitude(PulseAmp(simpulse->amplitude(), pmtData.m_gain, 00437 pmtData.m_sigmaGain)); 00438 prepulse->setAncestor(simpulse->ancestor()); 00439 prepulse->setType(PulseType::kPrePulse); 00440 prepulse->setChannelId(simpulse->channelId()); 00441 00442 return prepulse; 00443 }
DayaBay::ElecPulse * EsPmtEffectPulseTool::makeAfterPulse | ( | DayaBay::ElecPulse * | simpulse, | |
const DayaBay::PmtSimData & | pmtData | |||
) | [private] |
Definition at line 445 of file EsPmtEffectPulseTool.cc.
00447 { 00448 DayaBay::ElecPulse* afterpulse = new DayaBay::ElecPmtPulse(); 00449 00450 // Time offset from main pulse based on time PDF of after-pulses 00451 // double current_rand = m_randAfterPulseTime(); 00452 double current_rand = m_uni(); 00453 double afterPulseTime = ConvertPdfRand01(current_rand,m_afterPulsePdf,m_afterPulseEdges); 00454 // std::cout << "current_rand = " << current_rand << ", time = " << afterPulseTime << std::endl; 00455 00456 afterpulse->setTime(simpulse->time() + afterPulseTime); 00457 if(m_afterPulseAmpMode == "SinglePE") 00458 afterpulse->setAmplitude(PulseAmp(simpulse->amplitude(), pmtData.m_gain, pmtData.m_sigmaGain)); 00459 if(m_afterPulseAmpMode == "PDF"){ 00460 current_rand = m_uni(); 00461 double amplitude = ConvertPdfRand01(current_rand, m_afterPulseAmpPdf, m_afterPulseAmpEdges); 00462 afterpulse->setAmplitude(PulseAmp(amplitude, pmtData.m_gain, pmtData.m_sigmaGain)); 00463 } 00464 afterpulse->setAncestor(simpulse->ancestor()); 00465 afterpulse->setType(PulseType::kAfterPulse); 00466 afterpulse->setChannelId(simpulse->channelId()); 00467 00468 return afterpulse; 00469 }
DayaBay::ElecPulse * EsPmtEffectPulseTool::makeDarkPulse | ( | const DayaBay::ElecChannelId & | channelId, | |
const DayaBay::PmtSimData & | pmtData | |||
) | [private] |
Definition at line 471 of file EsPmtEffectPulseTool.cc.
00473 { 00474 DayaBay::ElecPulse* darkpulse = new DayaBay::ElecPmtPulse(); 00475 double darkPulseTime 00476 = m_uni()*double(m_simTimeLatest-m_simTimeEarliest)*Gaudi::Units::second; 00477 00478 darkpulse->setTime(darkPulseTime); 00479 if(m_darkPulseAmpMode == "SinglePE") darkpulse->setAmplitude(PulseAmp(1.0, pmtData.m_gain, pmtData.m_sigmaGain)); 00480 if(m_darkPulseAmpMode == "PDF"){ 00481 double current_rand = m_uni(); 00482 double amplitude = ConvertPdfRand01(current_rand, m_afterPulseAmpPdf, m_afterPulseAmpEdges); 00483 darkpulse->setAmplitude(PulseAmp(amplitude, pmtData.m_gain, pmtData.m_sigmaGain)); 00484 } 00485 darkpulse->setAncestor(0); 00486 darkpulse->setType(PulseType::kDarkPulse); 00487 darkpulse->setChannelId(channelId); 00488 00489 return darkpulse; 00490 }
float EsPmtEffectPulseTool::PulseAmp | ( | float | weight, | |
float | gain, | |||
float | gainFwhm | |||
) | [private] |
Definition at line 492 of file EsPmtEffectPulseTool.cc.
00492 { 00493 00494 // Include relative gain 00495 float randW = m_uni(); 00496 if (randW > m_expWeight ){ 00497 return m_gauss() * sigmaGain + gain; 00498 } 00499 else { 00500 return (m_exp() + m_speCutoff) * gain; 00501 } 00502 }
double EsPmtEffectPulseTool::ConvertPdfRand | ( | double | rand, | |
std::vector< double > | pdf, | |||
std::vector< double > | edges | |||
) | [private] |
Definition at line 505 of file EsPmtEffectPulseTool.cc.
00505 { 00506 // Defined PDF returns random number in [0,1] distributed according to user-defined histogram. 00507 // It assumes even bin sizes, so accomodate uneven bin sizes for generality. 00508 int current_bin; 00509 int Nbins = pdf.size(); 00510 00511 for (int bin=0;bin<Nbins;bin++) { // find which bin rand is in 00512 if (rand >= (double)(bin/Nbins) && rand < (double)((bin+1)/Nbins)) current_bin = bin; 00513 else current_bin = Nbins-1; // else rand=1, so current_bin is last bin (Nbins-1) 00514 } 00515 return edges[current_bin] + (rand*Nbins - current_bin)*(edges[current_bin+1]-edges[current_bin]); 00516 }
double EsPmtEffectPulseTool::ConvertPdfRand01 | ( | double | rand, | |
std::vector< double > | pdf, | |||
std::vector< double > | edges | |||
) | [private] |
Definition at line 518 of file EsPmtEffectPulseTool.cc.
00518 { 00519 // Defined PDF returns random number in [0,1] distributed according to user-defined histogram. 00520 // It assumes even bin sizes, so accomodate uneven bin sizes for generality. 00521 int current_bin; 00522 int Nbins = pdf.size(); 00523 00524 for(int bin=0; bin<Nbins; bin++) { 00525 if(rand >= pdf[bin] && rand < pdf[bin+1]) { 00526 current_bin = bin; 00527 break; 00528 } 00529 else 00530 current_bin = Nbins-1; 00531 } 00532 00533 return edges[current_bin] + (rand-pdf[current_bin])*(edges[current_bin+1]-edges[current_bin]) 00534 /(pdf[current_bin+1]-pdf[current_bin]); 00535 }
int EsPmtEffectPulseTool::PoissonRand | ( | double | mean | ) | [private] |
Definition at line 537 of file EsPmtEffectPulseTool.cc.
00537 { 00538 // Using source code from ROOT's TRandom::Poisson 00539 // Note: ROOT uses different algorithms depending on the mean, but mean is small 00540 // for our purposes, so use algorithm for mean<25 00541 00542 int n; 00543 if (mean <= 0) return 0; 00544 00545 double expmean = exp(-mean); 00546 double pir = 1; 00547 n = -1; 00548 while(1) { 00549 n++; 00550 pir *= m_uni(); 00551 if (pir <= expmean) break; 00552 } 00553 return n; 00554 }
double EsPmtEffectPulseTool::NumAfterPulse | ( | const int | numPmtHit | ) | [private] |
Definition at line 557 of file EsPmtEffectPulseTool.cc.
00557 { 00558 // Get number of after-pulses per PMT pulse based on the number of PMT hits 00559 00560 double cnt; 00561 00562 if (numPmtHit <= 400) 00563 cnt = 0.016 * numPmtHit; 00564 else if (numPmtHit <= 2500) 00565 cnt = -1.307853 + 0.02666278 * numPmtHit - 0.2001321e-4 * pow(numPmtHit, 2) 00566 +0.7330343e-8 * pow(numPmtHit, 3) - 0.102098e-11 * pow(numPmtHit, 4); 00567 else 00568 cnt = 15.; 00569 00570 return cnt; 00571 }
void EsPmtEffectPulseTool::getAfterPulseAmpPdf | ( | std::vector< double > & | ) | [private] |
Definition at line 589 of file EsPmtEffectPulseTool.cc.
00589 { 00590 pdf.push_back(0);pdf.push_back(0.0219574);pdf.push_back(0.0931247);pdf.push_back(0.179757); 00591 pdf.push_back(0.264803);pdf.push_back(0.342568);pdf.push_back(0.411712);pdf.push_back(0.472498); 00592 pdf.push_back(0.525729);pdf.push_back(0.572335);pdf.push_back(0.613205);pdf.push_back(0.649137); 00593 pdf.push_back(0.702162);pdf.push_back(0.745131);pdf.push_back(0.780295);pdf.push_back(0.809341); 00594 pdf.push_back(0.842677);pdf.push_back(0.868769);pdf.push_back(0.889482);pdf.push_back(0.906131); 00595 pdf.push_back(0.930777);pdf.push_back(0.947671);pdf.push_back(0.95962);pdf.push_back(0.968294); 00596 pdf.push_back(0.974731);pdf.push_back(0.983341);pdf.push_back(0.988564);pdf.push_back(0.99189); 00597 pdf.push_back(0.994094);pdf.push_back(0.995601);pdf.push_back(0.996661);pdf.push_back(0.997423); 00598 pdf.push_back(0.997983);pdf.push_back(0.998401);pdf.push_back(0.998717);pdf.push_back(0.999151); 00599 pdf.push_back(0.999418);pdf.push_back(0.999591);pdf.push_back(0.999705);pdf.push_back(0.999783); 00600 pdf.push_back(0.999838);pdf.push_back(0.999876);pdf.push_back(0.999905);pdf.push_back(0.999926); 00601 pdf.push_back(0.999941);pdf.push_back(0.999953);pdf.push_back(0.999962);pdf.push_back(0.999969); 00602 pdf.push_back(0.999975);pdf.push_back(0.999979);pdf.push_back(0.999983);pdf.push_back(1.0); 00603 }
void EsPmtEffectPulseTool::getAfterPulseAmpEdges | ( | std::vector< double > & | edges | ) | [private] |
Definition at line 573 of file EsPmtEffectPulseTool.cc.
00573 { 00574 edges.push_back(0.5);edges.push_back(0.7);edges.push_back(0.9);edges.push_back(1.1); 00575 edges.push_back(1.3);edges.push_back(1.5);edges.push_back(1.7);edges.push_back(1.9); 00576 edges.push_back(2.1);edges.push_back(2.3);edges.push_back(2.5);edges.push_back(2.7); 00577 edges.push_back(3.05);edges.push_back(3.4);edges.push_back(3.75);edges.push_back(4.1); 00578 edges.push_back(4.6);edges.push_back(5.1);edges.push_back(5.6);edges.push_back(6.1); 00579 edges.push_back(7.1);edges.push_back(8.1);edges.push_back(9.1);edges.push_back(10.1); 00580 edges.push_back(11.1);edges.push_back(13.1);edges.push_back(15.1);edges.push_back(17.1); 00581 edges.push_back(19.1);edges.push_back(21.1);edges.push_back(23.1);edges.push_back(25.1); 00582 edges.push_back(27.1);edges.push_back(29.1);edges.push_back(31.1);edges.push_back(35.1); 00583 edges.push_back(39.1);edges.push_back(43.1);edges.push_back(47.1);edges.push_back(51.1); 00584 edges.push_back(55.1);edges.push_back(59.1);edges.push_back(63.1);edges.push_back(67.1); 00585 edges.push_back(71.1);edges.push_back(75.1);edges.push_back(79.1);edges.push_back(83.1); 00586 edges.push_back(87.1);edges.push_back(91.1);edges.push_back(95.1);edges.push_back(100); 00587 }
const InterfaceID & IEsPulseTool::interfaceID | ( | ) | [static, inherited] |
Retrieve interface ID.
Reimplemented from IAlgTool.
Definition at line 8 of file IEsPulseTool.cc.
00009 { 00010 return IID_IEsPulseTool; 00011 }
std::string EsPmtEffectPulseTool::m_cableSvcName [private] |
Definition at line 51 of file EsPmtEffectPulseTool.h.
std::string EsPmtEffectPulseTool::m_simDataSvcName [private] |
Definition at line 55 of file EsPmtEffectPulseTool.h.
std::vector<double> EsPmtEffectPulseTool::m_prePulsePdf [private] |
Definition at line 60 of file EsPmtEffectPulseTool.h.
std::vector<double> EsPmtEffectPulseTool::m_prePulseEdges [private] |
Definition at line 61 of file EsPmtEffectPulseTool.h.
std::vector<double> EsPmtEffectPulseTool::m_afterPulsePdf [private] |
Definition at line 62 of file EsPmtEffectPulseTool.h.
std::vector<double> EsPmtEffectPulseTool::m_afterPulseEdges [private] |
Definition at line 63 of file EsPmtEffectPulseTool.h.
std::vector<double> EsPmtEffectPulseTool::m_afterPulseTime [private] |
Definition at line 64 of file EsPmtEffectPulseTool.h.
bool EsPmtEffectPulseTool::m_enablePrePulse [private] |
Definition at line 66 of file EsPmtEffectPulseTool.h.
bool EsPmtEffectPulseTool::m_enableAfterPulse [private] |
Definition at line 67 of file EsPmtEffectPulseTool.h.
std::string EsPmtEffectPulseTool::m_afterPulseAmpMode [private] |
Definition at line 72 of file EsPmtEffectPulseTool.h.
std::string EsPmtEffectPulseTool::m_darkPulseAmpMode [private] |
Definition at line 73 of file EsPmtEffectPulseTool.h.
std::vector<double> EsPmtEffectPulseTool::m_afterPulseAmpPdf [private] |
Definition at line 76 of file EsPmtEffectPulseTool.h.
std::vector<double> EsPmtEffectPulseTool::m_afterPulseAmpEdges [private] |
Definition at line 77 of file EsPmtEffectPulseTool.h.
bool EsPmtEffectPulseTool::m_enableNonlinearAfterpulse [private] |
Definition at line 81 of file EsPmtEffectPulseTool.h.
double EsPmtEffectPulseTool::m_linearAfterpulseThreshold [private] |
Definition at line 84 of file EsPmtEffectPulseTool.h.
bool EsPmtEffectPulseTool::m_enableVeryLongTimeSuppression [private] |
Definition at line 88 of file EsPmtEffectPulseTool.h.
double EsPmtEffectPulseTool::m_veryLongTime [private] |
Definition at line 91 of file EsPmtEffectPulseTool.h.
double EsPmtEffectPulseTool::m_expWeight [private] |
Definition at line 94 of file EsPmtEffectPulseTool.h.
double EsPmtEffectPulseTool::m_speExpDecay [private] |
Definition at line 96 of file EsPmtEffectPulseTool.h.
double EsPmtEffectPulseTool::m_speCutoff [private] |
Definition at line 97 of file EsPmtEffectPulseTool.h.
ICableSvc* EsPmtEffectPulseTool::m_cableSvc [private] |
Definition at line 100 of file EsPmtEffectPulseTool.h.
ISimDataSvc* EsPmtEffectPulseTool::m_simDataSvc [private] |
Definition at line 102 of file EsPmtEffectPulseTool.h.
Rndm::Numbers EsPmtEffectPulseTool::m_uni [private] |
Definition at line 105 of file EsPmtEffectPulseTool.h.
Rndm::Numbers EsPmtEffectPulseTool::m_gauss [private] |
Definition at line 105 of file EsPmtEffectPulseTool.h.
Rndm::Numbers EsPmtEffectPulseTool::m_exp [private] |
Definition at line 105 of file EsPmtEffectPulseTool.h.
Definition at line 106 of file EsPmtEffectPulseTool.h.
Definition at line 106 of file EsPmtEffectPulseTool.h.
Definition at line 109 of file EsPmtEffectPulseTool.h.
Definition at line 110 of file EsPmtEffectPulseTool.h.
double EsPmtEffectPulseTool::m_timeInterval [private] |
Definition at line 113 of file EsPmtEffectPulseTool.h.