In This Package:

MuonProphet Class Reference

Predict the fate of muons and the generation of spallation background. More...

#include <MuonProphet.h>

Inheritance diagram for MuonProphet:

[legend]Collaboration diagram for MuonProphet:[legend]List of all members.


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.

Detailed Description

Predict the fate of muons and the generation of spallation background.

All geometry related code

Zhe Wang, 12 Oct. 2009 at BNL

Definition at line 47 of file MuonProphet.h.

Constructor & Destructor Documentation

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]

Definition at line 54 of file MuonProphet.cc.

00055 {}

Member Function Documentation

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 }

Member Data Documentation

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]

random number

Definition at line 167 of file MuonProphet.h.

RectangularBox* MuonProphet::m_waterPool [private]

A simplified water pool geometry.

Definition at line 170 of file MuonProphet.h.

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]

top geoInfo

Definition at line 176 of file MuonProphet.h.

ILVolume* MuonProphet::m_topLVolume [private]

top logical volume

Definition at line 178 of file MuonProphet.h.

vector<IGeometryInfo *> MuonProphet::m_ADs [private]

hold the IGeometryInfo for each AD

Definition at line 180 of file MuonProphet.h.

ILVolume::Intersections MuonProphet::m_intersections [private]

A container to hold all intersections and materials.

Definition at line 183 of file MuonProphet.h.

ILVolume::Intersections MuonProphet::m_crucialSegments [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]

Some important materials to identify some volume.

Definition at line 189 of file MuonProphet.h.

Material* MuonProphet::m_owsWater [private]

RPC hit.

Definition at line 190 of file MuonProphet.h.

Material* MuonProphet::m_iwsWater [private]

OWS.

Definition at line 191 of file MuonProphet.h.

Material* MuonProphet::m_rock [private]

IWS.

Definition at line 192 of file MuonProphet.h.

Material* MuonProphet::m_mineralOil [private]

rock

Definition at line 193 of file MuonProphet.h.

Gaudi::XYZPoint MuonProphet::m_point [private]

muon track is expressed as x(t)= m_point + m_unit * t

Definition at line 201 of file MuonProphet.h.

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.

MpTriggerStat MuonProphet::m_poolTriggerStatus [private]

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]

Spallation neutron angular distribution.

Definition at line 254 of file MuonProphet.h.

TF1* MuonProphet::m_nLateralPDF [private]

Spallation background lateral profile universal to neutron and other isotopes.

Definition at line 257 of file MuonProphet.h.

The documentation for this class was generated from the following files:

/home/dayabay/installs/trunk_2011_04_11_opt/NuWa-trunk/dybgaudi/Generators/MuonProphet/src/components/MuonProphet.h
/home/dayabay/installs/trunk_2011_04_11_opt/NuWa-trunk/dybgaudi/Generators/MuonProphet/src/components/MpGeometry.cc
/home/dayabay/installs/trunk_2011_04_11_opt/NuWa-trunk/dybgaudi/Generators/MuonProphet/src/components/MpNeutron.cc
/home/dayabay/installs/trunk_2011_04_11_opt/NuWa-trunk/dybgaudi/Generators/MuonProphet/src/components/MpPlace.cc
/home/dayabay/installs/trunk_2011_04_11_opt/NuWa-trunk/dybgaudi/Generators/MuonProphet/src/components/MpSpallation.cc
/home/dayabay/installs/trunk_2011_04_11_opt/NuWa-trunk/dybgaudi/Generators/MuonProphet/src/components/MpTrigger.cc
/home/dayabay/installs/trunk_2011_04_11_opt/NuWa-trunk/dybgaudi/Generators/MuonProphet/src/components/MuonProphet.cc

Generated on Mon Apr 11 21:01:52 2011 for MuonProphet by 1.4.7

In This Package:

MuonProphet Class Reference

Public Types

Public Member Functions

Static Public Member Functions

Public Attributes

Protected Types

Protected Member Functions

Static Protected Attributes

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation