#include <MuonProphet.h>
Inheritance diagram for MuonProphet:
Public Types | |
SUCCESS | |
NO_INTERFACE | |
VERSMISMATCH | |
LAST_ERROR | |
SUCCESS | |
NO_INTERFACE | |
VERSMISMATCH | |
LAST_ERROR | |
enum | Status |
Public Member Functions | |
MuonProphet (const std::string &type, const std::string &name, const IInterface *parent) | |
virtual | ~MuonProphet () |
virtual StatusCode | initialize () |
virtual StatusCode | finalize () |
virtual StatusCode | mutate (HepMC::GenEvent &event) |
$$$$$ | |
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 () |
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 | |
MpMuonFate | predictMuonFate (HepMC::GenEvent &event) |
$$$$$ | |
StatusCode | geometryCalculationInit () |
StatusCode | geometryCalculation (HepMC::GenParticle *pMuon) |
MpTriggerStat | triggeredByRpc (HepMC::GenParticle *pMuon) |
MpTriggerStat | triggeredByWaterPool (HepMC::GenParticle *pMuon) |
int | spallationProducts (HepMC::GenEvent &event, int idx) |
$$$$$ The idx, index, is for parameter seaching in m_gen. | |
bool | generateOneBkg (int idx, HepMC::ThreeVector &bkgPoint) |
This will determine whether one type of background will be generated. | |
StatusCode | spallationNeutronsInit () |
Neutron production Neutron production will run on its own set of parameters' value. | |
int | spallationNeutrons (HepMC::GenEvent &event) |
Return number of tracks attached, negative for error. | |
int | generateNeutrons (HepMC::GenEvent *bkgEvt, ISolid::Tick tickBegin, ISolid::Tick tickEnd, double multiplicity) |
The code for generate neutrons in one segment. | |
StatusCode | printout (HepMC::GenEvent &event) |
void | setPosition (HepMC::GenEvent &event, HepMC::ThreeVector &position) |
.......................................... | |
void | setTime (HepMC::GenEvent &event, double t0, double lifetime) |
Private Attributes | |
bool | m_active |
HepMC::GenParticle * | m_muon |
vector< std::string > | m_genToolNames |
GtGenerator --------------------------- Tool to generate spallation background (non-neutron). | |
vector< int > | m_genId |
vector< IHepMCEventMutator * > | m_genTools |
vector< double > | m_genYields |
vector< double > | m_genYieldMeasuredAt |
vector< double > | m_genLifetimes |
bool | m_actNeutron |
double | m_neutronYield |
string | m_siteName |
Site::Site_t | m_site |
TRandom3 | m_rnd |
random number | |
RectangularBox * | m_waterPool |
A simplified water pool geometry. | |
string | m_topDetectorElementName |
Try the interfaces from lhcb DecDesc top detector name. | |
IGeometryInfo * | m_topGeoInfo |
top geoInfo | |
ILVolume * | m_topLVolume |
top logical volume | |
vector< IGeometryInfo * > | m_ADs |
hold the IGeometryInfo for each AD | |
ILVolume::Intersections | m_intersections |
A container to hold all intersections and materials. | |
ILVolume::Intersections | m_crucialSegments |
A container to hold only interesting intersections with thin material removed. | |
Material * | m_mixGas |
Some important materials to identify some volume. | |
Material * | m_owsWater |
RPC hit. | |
Material * | m_iwsWater |
OWS. | |
Material * | m_rock |
IWS. | |
Material * | m_mineralOil |
rock | |
Gaudi::XYZPoint | m_point |
muon track is expressed as x(t)= m_point + m_unit * t | |
Gaudi::XYZVector | m_vec |
Gaudi::XYZVector | m_unit |
bool | m_crossWater |
double | m_trkLengthInWater |
ISolid::Tick | m_tickWaterMin |
ISolid::Tick | m_tickWaterMax |
double | m_trkLengthInWaterThres |
double | m_waterPoolTriggerEff |
RpcPlane * | m_rpcPlane |
bool | m_crossRpc |
HepGeom::Point3D< double > | m_rpcIntersect |
MpTriggerStat | m_poolTriggerStatus |
MpTriggerStat | m_rpcTriggerStatus |
TF1 * | m_nEnergyPDF |
Spallation neutron energy distribution. | |
TF1 * | m_nAnglePDF |
Spallation neutron angular distribution. | |
TF1 * | m_nLateralPDF |
Spallation background lateral profile universal to neutron and other isotopes. |
All geometry related code
Zhe Wang, 12 Oct. 2009 at BNL
Definition at line 47 of file MuonProphet.h.
MuonProphet::MuonProphet | ( | const std::string & | type, | |
const std::string & | name, | |||
const IInterface * | parent | |||
) |
Definition at line 18 of file MuonProphet.cc.
00021 : GaudiTool(type,name,parent), 00022 m_muon(0), 00023 m_topGeoInfo(0), 00024 m_topLVolume(0) 00025 { 00026 declareInterface<IHepMCEventMutator>(this); 00027 // 00028 declareProperty("Active", m_active=true, "Use false to turn off this module"); 00029 // 00030 declareProperty("GenTools", m_genToolNames,"Tools to generate HepMC::GenEvents"); 00031 declareProperty("GenYields", m_genYields, "Yield for each generator"); 00032 declareProperty("GenYieldMeasuredAt",m_genYieldMeasuredAt,"The energy at which the yield is measured"); 00033 declareProperty("GenLifetimes", m_genLifetimes,"Lifetime of each generator"); 00034 // 00035 declareProperty("ActNeutron", m_actNeutron=true,"Turn on or off neutron production"); 00036 declareProperty("NeutronYields", m_neutronYield=-1,"Neutron yield. See MuonProphet.h for description"); 00037 // 00038 declareProperty("Site", m_siteName="DayaBay", "Site"); 00039 00041 declareProperty("TrkLengthInWaterThres", m_trkLengthInWaterThres, 00042 "If a track length in water is longer than this, it has high trigger rate."); 00043 00044 declareProperty("WaterPoolTriggerEff", m_waterPoolTriggerEff=-1, 00045 "Trigger efficiency in high trigger rate region, not in use!"); 00046 00048 declareProperty("TopDetectorElementName", 00049 m_topDetectorElementName="/dd/Structure/DayaBay", 00050 "Name of the DetectorElement that holds all others"); 00051 }
MuonProphet::~MuonProphet | ( | ) | [virtual] |
StatusCode MuonProphet::initialize | ( | ) | [virtual] |
Reimplemented from GaudiTool.
Definition at line 58 of file MuonProphet.cc.
00059 { 00060 info() << "Initializing MuonProphet" << endreq; 00061 00062 size_t s=m_genToolNames.size(); 00063 00064 if(s != m_genYields.size() || 00065 s != m_genLifetimes.size() || 00066 s != m_genYieldMeasuredAt.size() ) { 00067 error()<<"Imcomplete parameter set"<<endreq; 00068 return StatusCode::FAILURE; 00069 } 00070 00071 // print all user configuration 00072 info()<<"MuonProphet is active: "<<m_active<<endreq; 00073 m_site=Site::FromString(m_siteName.c_str()); 00074 info()<<"MuonProphet is set to site: "<<m_siteName<<" "<<m_site<<endreq; 00075 info()<<"Name || Yield || Yield measured at || Lifetime"<<endreq; 00076 for (size_t idx=0; idx<s; ++idx) { 00077 info() << m_genToolNames[idx]<< " " 00078 << m_genYields[idx]/(CLHEP::cm2/CLHEP::gram) <<"cm2/g " 00079 << m_genYieldMeasuredAt[idx]/CLHEP::GeV <<"GeV " 00080 << m_genLifetimes[idx]/CLHEP::second<<"s" <<endreq; 00081 } 00082 00083 // get all spallation background gentool 00084 string gtn; 00085 for (size_t idx=0; idx<s; ++idx) { 00086 gtn = m_genToolNames[idx]; 00087 try { 00088 m_genTools.push_back(tool<IHepMCEventMutator>(gtn)); 00089 } 00090 catch(const GaudiException& exg) { 00091 fatal() << "Failed to get generator: \"" << gtn << "\"" << endreq; 00092 return StatusCode::FAILURE; 00093 } 00094 info () << "Added spallation background gentool " << gtn << endreq; 00095 } 00096 00099 if(geometryCalculationInit().isFailure()) return StatusCode::FAILURE; 00101 // Spallation neutron related distribution initialization 00102 if(spallationNeutronsInit().isFailure()) return StatusCode::FAILURE; 00103 00104 00105 return StatusCode::SUCCESS; 00106 }
StatusCode MuonProphet::finalize | ( | ) | [virtual] |
Reimplemented from GaudiTool.
Definition at line 109 of file MuonProphet.cc.
00110 { 00111 return StatusCode::SUCCESS; 00112 }
StatusCode MuonProphet::mutate | ( | HepMC::GenEvent & | event | ) | [virtual] |
$$$$$
For all the rest of situations Attach some spallation background to muon track For all non-neutron products
At the end I add a tiny optical photon which won't affect any physics result. But it will make all the following program happy.
attach it to bkgEvt
Implements IHepMCEventMutator.
Definition at line 115 of file MuonProphet.cc.
00116 { 00117 // Not active. Don't bother to run it. 00118 if(!m_active) return StatusCode::SUCCESS; 00119 debug()<<"Processing a GenEvent"<<endreq; 00120 00121 printout(event); 00122 00123 m_crossRpc = false; 00124 m_crossWater = false; 00125 00126 // Predict the fate of muon and change it. 00127 MpMuonFate fate = predictMuonFate(event); 00128 debug()<<fate<<endreq; 00129 00130 if (fate.getCode() == MpMuonFate::kNeedSim) { 00131 // do nothing, return 00132 // to do: 00133 // But if geant4 can't produce these background, I still need to do it here. 00134 // For example, Li9 and He8. 00135 // Then it will be consistent within whole simulation scheme. 00136 // Right now no kNeedSim decision has been made. 00137 } else if (fate.getCode() == MpMuonFate::kUnknown) { 00138 // do nothing, return. Let geant4 to handle it. 00139 } else { 00143 int ntrk=0; 00144 size_t s=m_genTools.size(); 00145 00146 for (size_t idx=0; idx<s; ++idx) { 00147 00148 ntrk=spallationProducts(event, idx); 00149 00150 if(ntrk<0) { 00151 error()<<"failed to run one spallation background production"<<endreq; 00152 return StatusCode::FAILURE; 00153 } 00154 } 00155 00156 // spallation neutron 00157 // Since it will be triggered, won't be bothered with gamma, 00158 // because it can only cause prompt background and it is buried in muon signal. 00159 if( m_actNeutron ) { 00160 ntrk=spallationNeutrons(event); 00161 if(ntrk<0) { 00162 error()<<"failed to run spallation neutron background production"<<endreq; 00163 return StatusCode::FAILURE; 00164 } 00165 } 00166 } 00167 00170 HepMC::GenParticle* particle = new HepMC::GenParticle(HepMC::FourVector(0.01*CLHEP::eV, 00171 0, 00172 0, 00173 0.01*CLHEP::eV), 00174 22, 00175 1/*=status*/); 00177 (*(event.vertices_begin()))->add_particle_out(particle); 00178 00180 printout(event); 00181 00182 return StatusCode::SUCCESS; 00183 }
MpMuonFate MuonProphet::predictMuonFate | ( | HepMC::GenEvent & | event | ) | [private] |
$$$$$
Find the muon track here. All the following methods can use this pointer.
Definition at line 188 of file MuonProphet.cc.
00189 { 00190 00191 HepMC::GenEvent::particle_iterator p; 00192 HepMC::GenEvent::particle_iterator particleEnd=event.particles_end(); 00193 00194 int nMuon=0; 00195 00196 m_muon=0; 00199 for ( p = event.particles_begin(); p!=particleEnd; ++p) { 00200 00201 // muon pdg id is 13 or -13 00202 if( (*p)->pdg_id() == 13 || (*p)->pdg_id() == -13 ) { 00203 // Get the Muon pointer 00204 m_muon=*p; 00205 // In case of two or more than two muon tracks, don't work. 00206 ++nMuon; 00207 if(nMuon>=2) { 00208 return MpMuonFate::kUnknown; 00209 } 00210 } 00211 } 00212 00213 if(!m_muon) { 00214 info()<<"No muon track found"<<endreq; 00215 return MpMuonFate::kUnknown; 00216 } 00217 00218 // calculate all needed geometry parameters 00219 geometryCalculation(m_muon); 00220 00221 // get the trigger status of RPC and waterpool. 00222 // Note: No MpTriggerStat::kNeedSim will be issued. 00223 m_rpcTriggerStatus = triggeredByRpc(m_muon); 00224 // Note: No MpTriggerStat::kNeedSim will be issued. 00225 m_poolTriggerStatus = triggeredByWaterPool(m_muon); 00226 00227 debug()<<"RPC trigger status: "<<m_rpcTriggerStatus<<endreq; 00228 debug()<<"Pool trigger status: "<<m_poolTriggerStatus<<endreq; 00229 00230 00231 MpMuonFate fate(m_rpcTriggerStatus.getCode(), m_poolTriggerStatus.getCode()); 00232 m_muon->set_status(fate.getCode()); 00233 00234 return fate; 00235 00236 return MpMuonFate::kUnknown; 00237 }
StatusCode MuonProphet::geometryCalculationInit | ( | ) | [private] |
Definition at line 17 of file MpGeometry.cc.
00018 { 00019 /* 00021 double width_near = 10.0*CLHEP::meter ; // width of near pool 00022 double width_far = 16.0*CLHEP::meter ; // width of far pool 00023 double length = 16.0*CLHEP::meter ; // length of near and far pool 00024 double depth = 10.0*CLHEP::meter ; // depth of near and far pool 00025 double PoolDeadThickness = 84*CLHEP::mm; 00026 00027 // OWS box is constructed. 00028 if(m_site == Site::kDayaBay || m_site == Site::kLingAo) { 00029 00030 m_waterPool = new RectangularBox(length/2, -length/2, 00031 width_near/2, -width_near/2, 00032 depth/2, -depth/2+PoolDeadThickness); 00033 } else if ( m_site == Site::kFar ) { 00034 m_waterPool = new RectangularBox(length/2, -length/2, 00035 width_far/2, -width_far/2, 00036 depth/2, -depth/2+PoolDeadThickness); 00037 } 00038 00039 info()<<"Water pool geometry "<<*m_waterPool<<endreq; 00040 info()<<"TrkLengthInWaterThres= "<<m_trkLengthInWaterThres/CLHEP::cm<<"cm"<<endreq; 00041 00042 00044 double height_RPC = 0.50*CLHEP::meter ; // RPCs are nominally 50cm above surface of shield volume (includes 20cm of air above water) 00045 double thick_RPC = (0.056+0.030)*CLHEP::meter ; // 1.4cm/layer*4layers + 1cm spacing*3spaces between layers 00046 double air_thick = 0.2*CLHEP::meter ; // 20cm of air above water 00047 double extend_RPC = 1*CLHEP::meter; 00048 if(m_site == Site::kDayaBay || m_site == Site::kLingAo) { 00049 00050 m_rpcPlane = new RpcPlane(length/2+extend_RPC, -length/2-extend_RPC, 00051 width_near/2+extend_RPC, -width_near/2-extend_RPC, 00052 depth/2+air_thick+height_RPC+thick_RPC/2); 00053 } else if ( m_site == Site::kFar ) { 00054 00055 m_rpcPlane = new RpcPlane(length/2+extend_RPC, -length/2-extend_RPC, 00056 width_far/2+extend_RPC, -width_far/2-extend_RPC, 00057 depth/2+air_thick+height_RPC+thick_RPC/2); 00058 } 00059 info()<<"RPC geometry "<<*m_rpcPlane<<endreq; 00060 */ 00061 00063 // 1. DetectorDataSvc 00064 IDataProviderSvc* detSvc = 0; 00065 StatusCode sc = service("DetectorDataSvc",detSvc,true); 00066 if (sc.isFailure()) return sc; 00067 00068 // 2. Top Detector Element 00069 debug()<<"DetectorElement input: "<<m_topDetectorElementName<<endreq; 00070 SmartDataPtr<IDetectorElement> topDE(detSvc,m_topDetectorElementName); 00071 if (!topDE) return StatusCode::FAILURE; 00072 debug()<<"DetectorElement output: "<<topDE->name()<<endreq; 00073 00074 // 3. GeometryInfo 00075 m_topGeoInfo = topDE->geometry(); 00076 if(!m_topGeoInfo) return StatusCode::FAILURE; 00077 debug()<<"Top Geometry "<<m_topGeoInfo->lvolumeName()<<endreq; 00078 00079 // 4. get top LVolume 00080 const ILVolume* cILV = m_topGeoInfo->lvolume(); 00081 m_topLVolume = const_cast<ILVolume*>(cILV); 00082 if(!m_topLVolume) return StatusCode::FAILURE; 00083 info()<<"Got top logic volume: "<<m_topLVolume->name()<<endreq; 00084 00085 // 5. get the pointer to all charactoristic materials 00086 m_mixGas = 0; 00087 m_owsWater = 0; 00088 m_iwsWater = 0; 00089 m_mineralOil=0; 00090 m_rock = 0; 00091 00092 DataObject* pObj; 00093 00097 sc = detSvc->retrieveObject("/dd/Materials/MixGas", pObj); 00098 if(sc.isFailure()) return sc; 00099 m_mixGas = dynamic_cast<Material *> (pObj); 00100 //SmartDataPtr<Material> m_mixGas(detSvc,"/dd/Materials/MixGas"); 00101 if(!m_mixGas) return StatusCode::FAILURE; 00102 debug()<<"Got charactoristic material for RPC: "<<m_mixGas->name()<<endreq; 00103 00104 00105 00107 sc = detSvc->retrieveObject("/dd/Materials/OwsWater", pObj); 00108 if(sc.isFailure()) return sc; 00109 m_owsWater = dynamic_cast<Material *> (pObj); 00110 //SmartDataPtr<Material> m_mixGas(detSvc,"/dd/Materials/OwsWater"); 00111 if(!m_owsWater) return StatusCode::FAILURE; 00112 debug()<<"Got charactoristic material for OWS: "<<m_owsWater->name()<<endreq; 00113 00115 sc = detSvc->retrieveObject("/dd/Materials/IwsWater", pObj); 00116 if(sc.isFailure()) return sc; 00117 m_iwsWater = dynamic_cast<Material *> (pObj); 00118 //SmartDataPtr<Material> m_mixGas(detSvc,"/dd/Materials/IwsWater"); 00119 if(!m_iwsWater) return StatusCode::FAILURE; 00120 debug()<<"Got charactoristic material for IWS: "<<m_iwsWater->name()<<endreq; 00121 00122 // Rock 00123 sc = detSvc->retrieveObject("/dd/Materials/Rock", pObj); 00124 if(sc.isFailure()) return sc; 00125 m_rock = dynamic_cast<Material *> (pObj); 00126 //SmartDataPtr<Material> m_rock(detSvc,"/dd/Materials/Rock"); 00127 if(!m_rock) return StatusCode::FAILURE; 00128 debug()<<"Got charactoristic material for rock: "<<m_rock->name()<<endreq; 00129 00130 // MineralOil 00131 sc = detSvc->retrieveObject("/dd/Materials/MineralOil", pObj); 00132 if(sc.isFailure()) return sc; 00133 m_mineralOil = dynamic_cast<Material *> (pObj); 00134 //SmartDataPtr<Material> m_rock(detSvc,"/dd/Materials/MineralOil"); 00135 if(!m_mineralOil) return StatusCode::FAILURE; 00136 debug()<<"Got charactoristic material for AD: "<<m_mineralOil->name()<<endreq; 00137 00138 // debug()<<*m_mixGas<<endreq; 00139 // debug()<<int(m_mixGas)<<" "<<int(m_owsWater)<<" "<<int(m_rock)<<endreq; 00140 00142 vector<string> ADNames; 00143 00144 if(m_site == Site::kDayaBay) 00145 { 00146 ADNames.push_back("/dd/Structure/AD/db-oil1"); 00147 ADNames.push_back("/dd/Structure/AD/db-oil2"); 00148 } 00149 if(m_site == Site::kLingAo) 00150 { 00151 ADNames.push_back("/dd/Structure/AD/la-oil1"); 00152 ADNames.push_back("/dd/Structure/AD/la-oil2"); 00153 } 00154 if(m_site == Site::kFar) 00155 { 00156 ADNames.push_back("/dd/Structure/AD/far-oil1"); 00157 ADNames.push_back("/dd/Structure/AD/far-oil2"); 00158 ADNames.push_back("/dd/Structure/AD/far-oil3"); 00159 ADNames.push_back("/dd/Structure/AD/far-oil4"); 00160 } 00161 00162 IGeometryInfo* pGI; 00163 for(unsigned int idx = 0; idx < ADNames.size(); ++idx) 00164 { 00165 // Get Detector Element 00166 SmartDataPtr<IDetectorElement> pDE(detSvc,ADNames[idx]); 00167 if (!pDE) return StatusCode::FAILURE; 00168 debug()<<"Got DetectorElement: "<<pDE->name()<<endreq; 00169 00170 // Then geometryInfo 00171 pGI = pDE->geometry(); 00172 if(!pGI) return StatusCode::FAILURE; 00173 debug()<<"Got GeometryInfo: "<<pGI->lvolumeName()<<endreq; 00174 00175 m_ADs.push_back(pGI); 00176 } 00177 00178 return StatusCode::SUCCESS; 00179 00180 }
StatusCode MuonProphet::geometryCalculation | ( | HepMC::GenParticle * | pMuon | ) | [private] |
Definition at line 184 of file MpGeometry.cc.
00185 { 00187 /* 00188 HepGeom::Point3D<double> aPoint(pMuon->production_vertex()->position().x(), 00189 pMuon->production_vertex()->position().y(), 00190 pMuon->production_vertex()->position().z()); 00191 HepGeom::Point3D<double> aSlope(pMuon->momentum().px(), 00192 pMuon->momentum().py(), 00193 pMuon->momentum().pz()); 00194 00195 aSlope = aSlope.unit(); 00196 00197 m_crossWater = m_waterPool->intersect(aPoint,aSlope,m_crossWaterA,m_crossWaterB); 00198 00199 if(m_crossWater) { 00200 debug() << "Intersect with water OWS at " << m_crossWaterA <<" and "<< m_crossWaterB <<endreq; 00201 00202 m_trkLengthInWater = m_crossWaterA.distance(m_crossWaterB); 00203 00204 debug() << "Track length in water is "<<m_trkLengthInWater/CLHEP::cm<<"cm"<<endreq; 00205 } 00206 00207 00208 m_crossRpc = m_rpcPlane->intersect(aPoint,aSlope,m_rpcIntersect); 00209 00210 if(m_crossRpc) { 00211 debug()<<"Intersect with RPC at "<<m_rpcIntersect <<endreq; 00212 } 00213 */ 00214 00216 // 00217 // muon track is expressed as x(t)=m_point+m_unit*t 00218 00219 m_point.SetXYZ(pMuon->production_vertex()->position().x(), 00220 pMuon->production_vertex()->position().y(), 00221 pMuon->production_vertex()->position().z()); 00222 m_vec.SetXYZ(pMuon->momentum().px(), 00223 pMuon->momentum().py(), 00224 pMuon->momentum().pz()); 00225 m_unit = m_vec.Unit(); 00226 00228 //Gaudi::XYZPoint point(0,0,0); 00229 debug()<<"Point "<<m_point<<endreq; 00230 string path = m_topGeoInfo->belongsToPath(m_point,-1); 00231 debug()<<"Belong to path "<<path<<endreq; 00232 00233 //Gaudi::XYZPoint point(0,0,0); 00234 //Gaudi::XYZVector vec(1,0,0); 00235 00236 unsigned int N = m_topLVolume->intersectLine(m_point, 00237 m_unit, 00238 m_intersections, 00239 0,100*CLHEP::meter, 00240 -1); 00241 00242 debug()<<"Found "<<N<<" intersections"<<endreq; 00243 unsigned int idx; 00244 for (idx=0; idx<N; ++idx) { 00245 debug()<<"intersection: "<<m_intersections[idx].first<<" " 00246 <<m_intersections[idx].second->name()<<" " 00247 <<m_intersections[idx].second->density()<<endreq; 00248 //<<" "<<int(m_intersections[idx].second)<<endreq; 00249 } 00250 00251 // Stop muon. Stopping power = 2 [MeVcm2/g] 00252 // Begin with an approximation 00253 double ke = pMuon->momentum().e() - pMuon->momentum().m(); 00254 double sp = 2 * CLHEP::MeV * CLHEP::cm2 / CLHEP::g; 00255 double trkLength = ke/sp; 00256 //debug()<< CLHEP::MeV << " " << CLHEP::cm2 << " " << CLHEP::g <<endreq; 00257 debug()<<"ke= "<<ke<<" sp="<<sp<<" trkLength="<<trkLength<<endreq; 00258 00259 for (idx=0; idx<N; ++idx) { 00260 double Lpred = m_intersections[idx].first.second - m_intersections[idx].first.first; 00261 debug()<<"Length in this volume "<< Lpred*m_intersections[idx].second->density()<<endreq; 00262 if( Lpred*m_intersections[idx].second->density() < trkLength ) { 00263 trkLength -= Lpred*m_intersections[idx].second->density(); 00264 } else { // end of track 00265 m_intersections[idx].first.second = 00266 m_intersections[idx].first.first + 00267 trkLength / m_intersections[idx].second->density(); 00268 break; 00269 } 00270 } 00271 00272 if( idx<N-1 ) { // erase [ first, last ) 00273 m_intersections.erase( m_intersections.begin()+idx+1, m_intersections.end() ); 00274 } 00275 00276 N = m_intersections.size(); 00277 debug()<<"Left "<<N<<" intersections"<<endreq; 00278 for (idx=0; idx<N; ++idx) { 00279 debug()<<"intersection: "<<m_intersections[idx].first<<" " 00280 <<m_intersections[idx].second->name()<<" " 00281 <<m_intersections[idx].second->density()<<endreq; 00282 //<<" "<<int(m_intersections[idx].second)<<endreq; 00283 } 00284 00285 00286 00288 for (idx=0; idx<N; ++idx) { 00289 if(m_intersections[idx].second == m_mixGas) { 00290 m_crossRpc = true; 00291 break; 00292 } 00293 } 00294 00296 bool first = true; 00297 m_tickWaterMin = 0; 00298 m_tickWaterMax = 0; 00299 00300 for (idx=0; idx<N; ++idx) { 00301 if(m_intersections[idx].second == m_owsWater || 00302 m_intersections[idx].second == m_iwsWater ) { 00303 // first time find OwsWater 00304 if(first) { 00305 m_tickWaterMin = m_intersections[idx].first.first; 00306 m_tickWaterMax = m_intersections[idx].first.second; 00307 first = false; 00308 } else { 00309 if(m_intersections[idx].first.first < m_tickWaterMin ) { 00310 m_tickWaterMin = m_intersections[idx].first.first; 00311 } 00312 if(m_intersections[idx].first.second > m_tickWaterMax ) { 00313 m_tickWaterMax = m_intersections[idx].first.second; 00314 } 00315 } 00316 } 00317 } 00318 00319 00320 m_trkLengthInWater = m_tickWaterMax - m_tickWaterMin; 00321 debug()<<"Track length in water "<< m_trkLengthInWater/CLHEP::m <<"m"<<endreq; 00322 00324 if( m_trkLengthInWater>0 ) m_crossWater = true; 00326 00330 m_crucialSegments.clear(); 00331 00332 for (idx=0; idx<N; ++idx) { 00333 if(m_intersections[idx].second == m_rock || 00334 m_intersections[idx].second == m_owsWater || 00335 m_intersections[idx].second == m_iwsWater || 00336 m_intersections[idx].second == m_mineralOil ) 00337 { 00338 if(m_crucialSegments.empty()) 00339 { 00340 ILVolume::Intersection * seg = new ILVolume::Intersection(m_intersections[idx]); 00341 m_crucialSegments.push_back(*seg); 00342 } else { 00343 // check whether last one has the same material as current volume 00344 // iwsWater and owsWater are considered as the same 00345 if((m_crucialSegments.back().second == m_intersections[idx].second) || 00346 (m_intersections[idx].second == m_iwsWater && m_crucialSegments.back().second == m_owsWater) || 00347 (m_intersections[idx].second == m_owsWater && m_crucialSegments.back().second == m_iwsWater) ) 00348 { // If yes, combine them. Assume LS and oil are same 00349 m_crucialSegments.back().first.second = m_intersections[idx].first.second; 00350 } else { 00351 ILVolume::Intersection * seg = new ILVolume::Intersection(m_intersections[idx]); 00352 m_crucialSegments.push_back(*seg); 00353 } 00354 } 00355 } 00356 } 00357 00358 unsigned int NJiaodian = m_crucialSegments.size(); 00359 for (idx=0; idx<NJiaodian; ++idx) { 00360 debug()<<"intersection: "<<m_crucialSegments[idx].first<<" " 00361 <<m_crucialSegments[idx].second->name()<<endreq; 00362 // <<" "<<int(m_crucialSegments[idx].second)<<endreq; 00363 } 00364 00365 return StatusCode::SUCCESS; 00366 00367 }
MpTriggerStat MuonProphet::triggeredByRpc | ( | HepMC::GenParticle * | pMuon | ) | [private] |
Definition at line 25 of file MpTrigger.cc.
00026 { 00027 /* old one 00028 if(m_crossRpc) { 00029 return MpTriggerStat::kTriggered; 00030 } else { 00031 return MpTriggerStat::kNeedSim; 00032 } 00033 */ 00034 00036 if(m_crossRpc) { 00037 return MpTriggerStat::kTriggered; 00038 } else { 00039 return MpTriggerStat::kFarAway; 00040 } 00041 // 00042 }
MpTriggerStat MuonProphet::triggeredByWaterPool | ( | HepMC::GenParticle * | pMuon | ) | [private] |
Definition at line 45 of file MpTrigger.cc.
00046 { 00047 /* old one 00048 if(m_crossWater) { 00049 if(m_trkLengthInWater > m_trkLengthInWaterThres) { 00050 return MpTriggerStat::kTriggered; 00051 } 00052 } 00053 00054 return MpTriggerStat::kNeedSim; 00055 */ 00056 00058 if(m_crossWater) { 00059 if(m_trkLengthInWater > m_trkLengthInWaterThres) { // 100% triggered 00060 return MpTriggerStat::kTriggered; 00061 } else { 00062 if(m_rnd.Rndm()<0.5) { // 50% - 50% triggered or inefficent 00063 return MpTriggerStat::kTriggered; 00064 } else { 00065 return MpTriggerStat::kIneffi; 00066 } 00067 } 00068 } else { 00069 if(m_rnd.Rndm()<0.05) { // 5% triggered 95% faraway 00070 return MpTriggerStat::kTriggered; 00071 } else { 00072 return MpTriggerStat::kFarAway; 00073 } 00074 } 00076 }
int MuonProphet::spallationProducts | ( | HepMC::GenEvent & | event, | |
int | idx | |||
) | [private] |
$$$$$ The idx, index, is for parameter seaching in m_gen.
.. vector series. Return number of tracks attached, negative for error
Let's see if a background should be generated first. Should it be produced? Where is it?
Nothing generated, go back.
Second, start to work. Produce a spallation background
Fourth set its position and time
Fifth, give all vertices and particles to muon event
Definition at line 15 of file MpSpallation.cc.
00016 { 00017 if(m_muon == 0) { 00018 error()<<"No muon track found"<<endreq; 00019 // Since no muon track, it should not reach this point. 00020 return -1; 00021 } 00022 00025 HepMC::ThreeVector bkgPoint; 00026 bool bkgGenerated = false; 00027 bkgGenerated = generateOneBkg(idx, 00028 bkgPoint); 00029 00031 if( !bkgGenerated ) return 0; 00032 00034 HepMC::GenEvent * bkgEvt=new HepMC::GenEvent; 00035 if(m_genTools[idx]->mutate(*bkgEvt).isFailure()) { 00036 error()<<"Failed to run background gentool "<<m_genToolNames[idx]<<endreq; 00037 return -1; 00038 } 00039 00040 // debug 00041 #if 0 00042 { 00043 HepMC::GenEvent::particle_iterator p; 00044 HepMC::GenEvent::particle_iterator particleEnd = bkgEvt->particles_end(); 00045 for ( p = bkgEvt->particles_begin(); p!=particleEnd; ++p) { 00046 debug()<<*(*p)<<endreq; 00047 } 00048 HepMC::GenEvent::vertex_iterator v; 00049 HepMC::GenEvent::vertex_iterator vertexEnd = bkgEvt->vertices_end(); 00050 for ( v = bkgEvt->vertices_begin(); v!=vertexEnd; ++v) { 00051 debug()<<*(*v)<<endreq; 00052 } 00053 } 00054 #endif 00055 00057 // change the vertex position first 00058 setPosition(*bkgEvt,bkgPoint); 00059 00060 // change the time of the vertex 00061 // Get muon vertex time 00062 double t0 = m_muon->production_vertex()->position().t(); 00063 // to do: A delta should be added to t0 to account for flight time of muon 00064 // and flight time of muon's daughter particle which reached the 00065 // production vertex of the background. 00066 00067 setTime(*bkgEvt,t0,m_genLifetimes[idx]); 00068 00070 // drag vertices out from background event 00071 HepMC::GenVertex *bkgVtx; 00072 while( !bkgEvt->vertices_empty() ){ 00073 00074 bkgVtx = *(bkgEvt->vertices_begin()); 00075 00076 // add to muon event 00077 // This will automatically transfer the ownership from the old event to new event. 00078 // So the original event can be safely deleted now. 00079 event.add_vertex(bkgVtx); 00080 } 00081 00082 // delete the bkgEvt 00083 delete bkgEvt; 00084 00085 return 0; 00086 }
bool MuonProphet::generateOneBkg | ( | int | idx, | |
HepMC::ThreeVector & | bkgPoint | |||
) | [private] |
This will determine whether one type of background will be generated.
If generated it will return the position which should be in one AD. For non-neutron background only
if a track has never touch water, that means the track is far away from AD. Don't bother with this situation. Those heavy nuclei won't fly that long.
get yield
get probability
less than the predicted probability
Let's generate something now. If it failed get into an AD, it still will fail. Let's build a muon shower lateral profile function now.
First of all, get the values of all these parameters
Transfer and rotate it to global coordinate system according to muon's oritation
The last thing is to see whether it is in an AD or not. gaudi use its point and vector, not compatible with others. CLHEP
Definition at line 90 of file MpSpallation.cc.
00092 { 00095 if( !m_crossWater ) { 00096 return false; 00097 } 00098 00100 double eMuon = m_muon->momentum().e(); 00101 00102 // Mineral oil and LS densities are assumed to be the same. 00103 double density = m_mineralOil->density(); 00104 double yield = m_genYields[idx] ; 00105 double alpha = 0.77; 00106 00108 double yieldAtThisEnergy = yield * pow( (eMuon / m_genYieldMeasuredAt[idx]), alpha); 00109 00111 double ProbThres = yieldAtThisEnergy * density * m_trkLengthInWater; 00112 00113 debug()<<m_genToolNames[idx]<<" eMuon= "<<eMuon/CLHEP::GeV<<"GeV " 00114 <<" trkLengthInWater= "<<m_trkLengthInWater/CLHEP::meter<<"m " 00115 <<"ProbThres= "<<ProbThres<<endreq; 00116 00117 if( ProbThres > 1 || ProbThres < 0 ) { 00118 debug()<<"Yield too high or too low. Invalid probability "<<ProbThres<<endreq; 00119 } 00120 00121 if( m_rnd.Rndm() > ProbThres) { 00123 return false; 00124 } 00125 00129 00130 // The position of background in the cylinder coordinate system of muon track 00131 // Muon direction is z. 00132 double r,phi,z; 00133 00135 z = m_rnd.Rndm() * (m_tickWaterMax - m_tickWaterMin) + m_tickWaterMin; 00136 r = m_nLateralPDF->GetRandom() * CLHEP::cm ; 00137 phi = m_rnd.Rndm() * 2 * 3.14159265358979323846; 00138 00139 00140 // try point in coordinate system of intersect 1st 00141 HepGeom::Point3D<double> tryPoint; 00142 00143 tryPoint.setX( r*cos(phi) ); 00144 tryPoint.setY( r*sin(phi) ); 00145 tryPoint.setZ( z ); 00146 00148 /* This is a muon track shown below. 00149 * 00150 * intersect 1st (tick1) intersect 2nd (tick2) 00151 * Primary vertex -----1--------C----------2------> 00152 * z | theta 00153 * | 00154 * | r 00155 * | 00156 * | 00157 * tryPoint 00158 * 00159 */ 00160 00161 // direction of muon 00162 HepGeom::Point3D<double> muonDir(m_unit.X(), 00163 m_unit.Y(), 00164 m_unit.Z()); 00165 00166 // need a y axis which is perpendicular to muonSlope 00167 HepGeom::Point3D<double> y( m_unit.X(), 00168 m_unit.Y(), 00169 0); 00170 y.rotateZ(3.14159265358979323846/2); 00171 if(y.r()<=0) y.setX(1); 00172 00173 // tryPoint in ows coordinate system 00174 HepGeom::Point3D<double> tryPointGlobal = muonDir * tryPoint.r(); 00175 00176 // rotation 00177 tryPointGlobal.rotate(tryPoint.theta(), y); 00178 tryPointGlobal.rotate(tryPoint.phi(), muonDir); 00179 00180 // transform 00181 tryPointGlobal.setX(tryPointGlobal.x() + m_point.X()); 00182 tryPointGlobal.setY(tryPointGlobal.y() + m_point.Y()); 00183 tryPointGlobal.setZ(tryPointGlobal.z() + m_point.Z()); 00184 00185 debug()<<tryPoint<<endreq; 00186 debug()<<tryPointGlobal<<endreq; 00187 00190 Gaudi::XYZPoint p(tryPointGlobal.x(), 00191 tryPointGlobal.y(), 00192 tryPointGlobal.z()); 00193 00194 bool inside = false; 00195 for (unsigned int idx = 0; idx<m_ADs.size(); ++idx) 00196 { 00197 if(m_ADs[idx]->isInside(p)) inside = true; 00198 } 00199 00200 if(!inside) return false; 00201 00202 // done 00203 bkgPoint.setX(tryPointGlobal.x()); 00204 bkgPoint.setY(tryPointGlobal.y()); 00205 bkgPoint.setZ(tryPointGlobal.z()); 00206 00207 return true; 00208 }
StatusCode MuonProphet::spallationNeutronsInit | ( | ) | [private] |
Neutron production Neutron production will run on its own set of parameters' value.
./functions
Zhe Wang 24 Nov, 2009
In principle the following two distributions are muon energy dependent. They will vary as the muon's energy changes. However this is simplified in simulation. Use the mean muon energy instead.
Spallation neutron energy distribution Energy range: 0.001 - 5 GeV 0 is a pole for this function Number of parameter: One
This must be big, else some energy will be missing.
Wierd root problem
Spallation neutron angular distribution cos(thera) range: -1 to 1 Number of parameter: One
Get the lateral profile PDF
Definition at line 12 of file MpNeutron.cc.
00013 { 00014 m_nEnergyPDF=0; 00015 m_nAnglePDF=0; 00016 m_nLateralPDF=0; 00017 00022 00023 // Avarage muon energy is from the estimation of Mengyun Guan. (Ph.D. thesis) 00024 double avaMuonE; 00025 if(m_site == Site::kDayaBay) { 00026 avaMuonE = 55.3*CLHEP::GeV; 00027 } else if (m_site == Site::kLingAo) { 00028 avaMuonE = 61.4*CLHEP::GeV; 00029 } else if (m_site == Site::kFar ) { 00030 avaMuonE = 140.3*CLHEP::GeV; 00031 } 00032 00037 m_nEnergyPDF = new TF1("m_nEnergyPDF",NeutronEnergyPDF,0.001,4,1); 00038 m_nEnergyPDF->SetNpx(200); 00039 00040 m_nEnergyPDF->SetParameter(0,avaMuonE/CLHEP::GeV); 00041 if(!m_nEnergyPDF) { 00042 error()<<"Failed to get NeutronEnergyPDF"<<endreq; 00043 return StatusCode::FAILURE; 00044 } 00045 00049 m_nAnglePDF = new TF1("m_nAnglePDF",NeutronAngularDis,-1,1,1); 00050 m_nAnglePDF->SetParameter(0,avaMuonE/CLHEP::GeV); 00051 if(!m_nAnglePDF) { 00052 error()<<"Failed to get NeutronAngularDis"<<endreq; 00053 return StatusCode::FAILURE; 00054 } 00055 00056 00058 m_nLateralPDF = new TF1("m_nLateralPDF",LateralProfile,0,100,1); 00059 m_nLateralPDF->SetParameter(0,0); 00060 if(!m_nLateralPDF) { 00061 error()<<"Failed to get LateralProfile"<<endreq; 00062 return StatusCode::FAILURE; 00063 } 00064 00065 return StatusCode::SUCCESS; 00066 }
int MuonProphet::spallationNeutrons | ( | HepMC::GenEvent & | event | ) | [private] |
Return number of tracks attached, negative for error.
Build a background event, then later add its vertex(es) to the orginal event
Probability = yield * track length * density / multiplicity
generate some neutrons
Move all vertices and particles to muon event
Definition at line 70 of file MpNeutron.cc.
00071 { 00072 debug()<<"Spallation Neutron production"<<endreq; 00073 00074 //-------------------------------------------------------- 00075 unsigned int NSeg = m_crucialSegments.size(); 00076 00077 double yield; 00078 if(m_neutronYield>=0) { 00079 yield = m_neutronYield; 00080 } else { 00081 yield = NeutronYield( m_muon->momentum().e()/CLHEP::GeV )*(CLHEP::cm2/CLHEP::gram); 00082 } 00083 00084 double tracklength; 00085 double multiplicity; 00086 00087 int nNewTrack = 0; 00088 00090 HepMC::GenEvent * bkgEvt=new HepMC::GenEvent; 00091 00092 for (unsigned int idx=0; idx<NSeg; ++idx) 00093 { 00094 double prob; 00095 00101 00102 tracklength = m_crucialSegments[idx].first.second - m_crucialSegments[idx].first.first; 00103 multiplicity = NeutronMulti( m_muon->momentum().e()/CLHEP::GeV, 00104 tracklength ); 00105 00106 prob = yield 00107 * tracklength 00108 * m_crucialSegments[idx].second->density() 00109 / multiplicity; 00110 00111 debug()<<"Yield= "<<NeutronYield( m_muon->momentum().e()/CLHEP::GeV )<<"cm2/g " 00112 <<"TrackLength= "<< tracklength/CLHEP::m <<"m " 00113 <<"Density= "<<m_crucialSegments[idx].second->density()/(CLHEP::g/CLHEP::cm3)<<"g/cm3 " 00114 <<"Multiplicity= "<< multiplicity <<endreq; 00115 debug()<<"prob= "<<prob<<endreq; 00116 00117 00118 00119 if( m_rnd.Rndm() < prob ) 00120 { 00122 nNewTrack += generateNeutrons(bkgEvt, 00123 m_crucialSegments[idx].first.first, // tick one 00124 m_crucialSegments[idx].first.second, // tick two 00125 multiplicity); 00126 } 00127 } 00128 00130 // drag vertices out from background event 00131 00132 HepMC::GenVertex *bkgVtx; 00133 while( !bkgEvt->vertices_empty() ) 00134 { 00135 bkgVtx = *(bkgEvt->vertices_begin()); 00136 00137 // add to muon event 00138 // This will automatically transfer the ownership from the old event to new event. 00139 // So the original event can be safely deleted now. 00140 event.add_vertex(bkgVtx); 00141 } 00142 00143 // delete the bkgEvt 00144 delete bkgEvt; 00145 00146 return nNewTrack; 00147 00148 }
int MuonProphet::generateNeutrons | ( | HepMC::GenEvent * | bkgEvt, | |
ISolid::Tick | tickBegin, | |||
ISolid::Tick | tickEnd, | |||
double | multiplicity | |||
) | [private] |
The code for generate neutrons in one segment.
The position of background in the local cylinder coordinate system of muon track Muon direction is in z direction. Muon vertex, m_point, is the origin.
(z,r,phi) is the vertex position of neutron
eng (ke, momentum), Neutron energy, (phiN,costh) gives the neutron direction
First of all, get the values of all these parameters
to do: build correlation with neutron energy
Second, put a neutron to its place
Neutron vertex : local
Neutron momentum vector : local
transfer to muon's global coordination system
transfer vertex first
rotate the neutron momentum direction
Create a neutron vertex and track
attach it to bkgEvt
Definition at line 155 of file MpNeutron.cc.
00159 { 00160 int nNeutron = m_rnd.Poisson(multiplicity); 00161 00162 debug()<<"nNeutron= "<<nNeutron<<endreq; 00163 00166 00168 double z; 00169 double r; 00170 double phi; 00171 00173 double ke; 00174 double eng; 00175 double momentum; 00176 double massN = 939.565346*CLHEP::MeV; 00177 double phiN; 00178 double costh; 00179 double sinth; 00180 00181 for(int idx =0; idx<nNeutron; ++idx) 00182 { 00184 z = m_rnd.Rndm() * (tickEnd - tickBegin) + tickBegin; 00185 r = m_nLateralPDF->GetRandom() * CLHEP::cm ; 00186 phi = m_rnd.Rndm() * 2 * 3.14159265358979323846; 00187 00188 ke = m_nEnergyPDF->GetRandom() * CLHEP::GeV ; 00189 eng = ke + massN; 00190 momentum = sqrt (eng*eng - massN*massN); 00191 phiN = m_rnd.Rndm() * 2 * 3.14159265358979323846; 00193 costh = m_nAnglePDF->GetRandom(); 00194 sinth = sqrt(1-costh*costh); 00195 00197 00198 debug()<<"z "<<z/CLHEP::m<<"m"<<endreq; 00199 debug()<<"r "<<r/CLHEP::cm<<"cm"<<endreq; 00200 debug()<<"ke "<<ke/CLHEP::GeV<<"GeV"<<endreq; 00201 debug()<<"costh "<<costh<<endreq; 00202 00204 HepGeom::Point3D<double> genPoint; 00205 00206 genPoint.setX( r*cos(phi) ); 00207 genPoint.setY( r*sin(phi) ); 00208 genPoint.setZ( z ); 00209 00211 HepGeom::Point3D<double> genVec; 00212 00213 genVec.setX ( momentum*sinth*cos(phiN) ); 00214 genVec.setY ( momentum*sinth*sin(phiN) ); 00215 genVec.setZ ( momentum*costh ); 00216 00218 HepGeom::Point3D<double> muonDir(m_unit.X(), 00219 m_unit.Y(), 00220 m_unit.Z()); 00221 00222 // need a y axis which is perpendicular to muonSlope 00223 HepGeom::Point3D<double> y( m_unit.X(), 00224 m_unit.Y(), 00225 0); 00226 y.rotateZ(3.14159265358979323846/2); 00227 if(y.r()<=0) y.setX(1); 00228 00230 HepGeom::Point3D<double> genPointGlobal = muonDir * genPoint.r(); 00231 debug()<<m_unit<<endreq; 00232 debug()<<muonDir<<endreq; 00233 debug()<<genPointGlobal<<endreq; 00234 00235 // rotation 00236 genPointGlobal.rotate(genPoint.theta(), y); 00237 genPointGlobal.rotate(genPoint.phi(), muonDir); 00238 00239 // transform 00240 genPointGlobal.setX(genPointGlobal.x() + m_point.X()); 00241 genPointGlobal.setY(genPointGlobal.y() + m_point.Y()); 00242 genPointGlobal.setZ(genPointGlobal.z() + m_point.Z()); 00243 00245 HepGeom::Point3D<double> genVectorGlobal = muonDir * genVec.r(); 00246 00247 genVectorGlobal.rotate(genVec.theta(), y); 00248 genVectorGlobal.rotate(genVec.phi(), muonDir); 00249 00250 00252 double t0 = m_muon->production_vertex()->position().t(); 00253 debug()<<"t0= "<<t0<<endreq; 00254 t0 += tickEnd / (3e8 * CLHEP::m/CLHEP::s); 00255 debug()<<"dt= "<<tickEnd / (3e8 * CLHEP::m/CLHEP::s)<<endreq; 00256 debug()<<"t0= "<<t0<<endreq; 00257 HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(genPointGlobal.x(), 00258 genPointGlobal.y(), 00259 genPointGlobal.z(), 00260 t0)); 00261 00262 HepMC::GenParticle* particle = new HepMC::GenParticle(HepMC::FourVector(genVectorGlobal.x(), 00263 genVectorGlobal.y(), 00264 genVectorGlobal.z(), 00265 eng), 00266 2112, 00267 1/*=status*/); 00269 vertex->add_particle_out(particle); 00270 bkgEvt->add_vertex(vertex); 00271 } 00272 00273 return nNeutron; 00274 }
StatusCode MuonProphet::printout | ( | HepMC::GenEvent & | event | ) | [private] |
Definition at line 240 of file MuonProphet.cc.
00241 { 00242 // every vertex 00243 HepMC::GenEvent::vertex_iterator vtx, done = event.vertices_end(); 00244 for (vtx = event.vertices_begin(); vtx != done; ++vtx) { 00245 debug()<<*(*vtx)<<endreq; 00246 // every particle 00247 HepMC::GenVertex::particles_in_const_iterator p; 00248 HepMC::GenVertex::particles_in_const_iterator inEnd=(*vtx)->particles_in_const_end(); 00249 for ( p = (*vtx)->particles_in_const_begin(); p!=inEnd; ++p) { 00250 debug()<<*(*p)<<endreq; 00251 } 00252 HepMC::GenVertex::particles_out_const_iterator outEnd=(*vtx)->particles_out_const_end(); 00253 for ( p = (*vtx)->particles_out_const_begin(); p!=outEnd; ++p) { 00254 debug()<<*(*p)<<endreq; 00255 } 00256 } 00257 return StatusCode::SUCCESS; 00258 }
void MuonProphet::setPosition | ( | HepMC::GenEvent & | event, | |
HepMC::ThreeVector & | position | |||
) | [private] |
..........................................
Definition at line 13 of file MpPlace.cc.
00014 { 00015 HepMC::FourVector position(aPoint.x(), 00016 aPoint.y(), 00017 aPoint.z(), 00018 0); 00019 00020 HepMC::GenEvent::vertex_iterator vtx, done = event.vertices_end(); 00021 for (vtx = event.vertices_begin(); vtx != done; ++vtx) { 00022 00023 // still use the original vertex time 00024 position.setT( (*vtx)->position().t() ); 00025 00026 (*vtx)->set_position(position); 00027 } 00028 00029 }
void MuonProphet::setTime | ( | HepMC::GenEvent & | event, | |
double | t0, | |||
double | lifetime | |||
) | [private] |
Definition at line 33 of file MpPlace.cc.
00034 { 00035 // Use negtive lifetime to skip the exponential random lifetime 00036 00037 // 00038 // debug()<<"t0= "<<t0<<endreq; 00039 // 00040 double dt = 0; 00041 00042 if(lifetime<0) { 00043 dt=0; 00044 } else { 00045 double u = m_rnd.Rndm(); 00046 // Exponential Distribution: t(n) = t(n-1) * exp(-t/tL) 00047 dt = ((-1.0 * log(u)) * lifetime); 00048 } 00049 00051 dt = dt / CLHEP::nanosecond; 00052 00053 HepMC::GenEvent::vertex_iterator vtx, done = event.vertices_end(); 00054 for (vtx = event.vertices_begin(); vtx != done; ++vtx) { 00055 00056 HepMC::FourVector position = (*vtx)->position(); 00057 00058 // relative to muon + mother particle's lifetime + internal delay 00059 double vertex_time = t0+dt+position.t(); 00060 00061 position.setT(vertex_time); 00062 (*vtx)->set_position(position); 00063 } 00064 }
bool MuonProphet::m_active [private] |
Definition at line 114 of file MuonProphet.h.
HepMC::GenParticle* MuonProphet::m_muon [private] |
Definition at line 117 of file MuonProphet.h.
vector<std::string> MuonProphet::m_genToolNames [private] |
GtGenerator --------------------------- Tool to generate spallation background (non-neutron).
Definition at line 122 of file MuonProphet.h.
vector<int> MuonProphet::m_genId [private] |
Definition at line 124 of file MuonProphet.h.
vector<IHepMCEventMutator*> MuonProphet::m_genTools [private] |
Definition at line 126 of file MuonProphet.h.
vector<double> MuonProphet::m_genYields [private] |
Definition at line 142 of file MuonProphet.h.
vector<double> MuonProphet::m_genYieldMeasuredAt [private] |
Definition at line 144 of file MuonProphet.h.
vector<double> MuonProphet::m_genLifetimes [private] |
Definition at line 146 of file MuonProphet.h.
bool MuonProphet::m_actNeutron [private] |
Definition at line 150 of file MuonProphet.h.
double MuonProphet::m_neutronYield [private] |
Definition at line 158 of file MuonProphet.h.
string MuonProphet::m_siteName [private] |
Definition at line 163 of file MuonProphet.h.
Site::Site_t MuonProphet::m_site [private] |
Definition at line 164 of file MuonProphet.h.
TRandom3 MuonProphet::m_rnd [private] |
RectangularBox* MuonProphet::m_waterPool [private] |
string MuonProphet::m_topDetectorElementName [private] |
Try the interfaces from lhcb DecDesc top detector name.
Definition at line 174 of file MuonProphet.h.
IGeometryInfo* MuonProphet::m_topGeoInfo [private] |
ILVolume* MuonProphet::m_topLVolume [private] |
vector<IGeometryInfo *> MuonProphet::m_ADs [private] |
A container to hold only interesting intersections with thin material removed.
Only Rock, OwsWater, MineralOil are left
Definition at line 186 of file MuonProphet.h.
Material* MuonProphet::m_mixGas [private] |
Material* MuonProphet::m_owsWater [private] |
Material* MuonProphet::m_iwsWater [private] |
Material* MuonProphet::m_rock [private] |
Material* MuonProphet::m_mineralOil [private] |
Gaudi::XYZPoint MuonProphet::m_point [private] |
Gaudi::XYZVector MuonProphet::m_vec [private] |
Definition at line 202 of file MuonProphet.h.
Gaudi::XYZVector MuonProphet::m_unit [private] |
Definition at line 203 of file MuonProphet.h.
bool MuonProphet::m_crossWater [private] |
Definition at line 207 of file MuonProphet.h.
double MuonProphet::m_trkLengthInWater [private] |
Definition at line 215 of file MuonProphet.h.
ISolid::Tick MuonProphet::m_tickWaterMin [private] |
Definition at line 218 of file MuonProphet.h.
ISolid::Tick MuonProphet::m_tickWaterMax [private] |
Definition at line 219 of file MuonProphet.h.
double MuonProphet::m_trkLengthInWaterThres [private] |
Definition at line 223 of file MuonProphet.h.
double MuonProphet::m_waterPoolTriggerEff [private] |
Definition at line 232 of file MuonProphet.h.
RpcPlane* MuonProphet::m_rpcPlane [private] |
Definition at line 235 of file MuonProphet.h.
bool MuonProphet::m_crossRpc [private] |
Definition at line 237 of file MuonProphet.h.
HepGeom::Point3D<double> MuonProphet::m_rpcIntersect [private] |
Definition at line 238 of file MuonProphet.h.
Definition at line 241 of file MuonProphet.h.
MpTriggerStat MuonProphet::m_rpcTriggerStatus [private] |
Definition at line 242 of file MuonProphet.h.
TF1* MuonProphet::m_nEnergyPDF [private] |
Spallation neutron energy distribution.
They will vary as the muon's energy changes. However this is simplified in simulation. Use the mean muon energy instead.
Definition at line 252 of file MuonProphet.h.
TF1* MuonProphet::m_nAnglePDF [private] |
TF1* MuonProphet::m_nLateralPDF [private] |
Spallation background lateral profile universal to neutron and other isotopes.
Definition at line 257 of file MuonProphet.h.