/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

Public Member Functions | Static Public Member Functions | Private Member Functions | Private Attributes
MuonProphet Class Reference

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

#include <MuonProphet.h>

Inheritance diagram for MuonProphet:
Inheritance graph
[legend]
Collaboration diagram for MuonProphet:
Collaboration graph
[legend]

List of all members.

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)
 $$$$$

Static Public Member Functions

static const InterfaceID & interfaceID ()

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::GenParticlem_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
RectangularBoxm_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
 MineralOil.
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
RpcPlanem_rpcPlane
bool m_crossRpc
HepGeom::Point3D< double > m_rpcIntersect
MpTriggerStat m_poolTriggerStatus
MpTriggerStat m_rpcTriggerStatus
TF1 * m_nEnergyPDF
 In principle the following two distributions are muon energy dependent.
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.

Spallation production related part.

Place a background into its right space and time position, including altering its momentum.

All geometry related code

Zhe Wang, 12 Oct. 2009 at BNL

Zhe Wang, 25 Nov. 2009 at BNL

Zhe Wang, 15 Oct. 2009 at BNL

All Trigger related code. Now only two statuses are implemented; kTriggered and kNeedSim. to do: A more sophisticated logic is needed. kIneffi, kFarAway ...

Zhe Wang, 2 Oct. 2009 at BNL

Determine the fate of cosmogenic muons: kTriggered, kMissed, kUnknown. For the case of kUnknown, a full simulation will be invoked for it, while for kTriggered and kMissed, call a fast muon simulation.

Zhe Wang, 2 Oct. 2009 at BNL

Determine the fate of cosmogenic muons: kUnknown, kTriggered, kMissed and kNeedSim. For the case of kNeedSim, a full simulation will be invoked for it, while for kTriggered and kMissed, call a fast muon simulation. kUnknown is reserved for errors, like no muon track is found.

Zhe Wang, 2 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 
)

Trigger related configuration

Definition at line 18 of file MuonProphet.cc.

  : GaudiTool(type,name,parent),
    m_muon(0),
    m_topGeoInfo(0),
    m_topLVolume(0)
{
  declareInterface<IHepMCEventMutator>(this);
  //
  declareProperty("Active",            m_active=true, "Use false to turn off this module");
  //
  declareProperty("GenTools",          m_genToolNames,"Tools to generate HepMC::GenEvents");
  declareProperty("GenYields",         m_genYields,   "Yield for each generator");
  declareProperty("GenYieldMeasuredAt",m_genYieldMeasuredAt,"The energy at which the yield is measured");
  declareProperty("GenLifetimes",      m_genLifetimes,"Lifetime of each generator");
  //
  declareProperty("ActNeutron",        m_actNeutron=true,"Turn on or off neutron production");
  declareProperty("NeutronYields",     m_neutronYield=-1,"Neutron yield. See MuonProphet.h for description");
  //
  declareProperty("Site",              m_siteName="DayaBay", "Site");

  declareProperty("TrkLengthInWaterThres", m_trkLengthInWaterThres, 
                  "If a track length in water is longer than this, it has high trigger rate.");

  declareProperty("WaterPoolTriggerEff",   m_waterPoolTriggerEff=-1,
                  "Trigger efficiency in high trigger rate region, not in use!");

  declareProperty("TopDetectorElementName",
                  m_topDetectorElementName="/dd/Structure/DayaBay",
                  "Name of the DetectorElement that holds all others");
}
MuonProphet::~MuonProphet ( ) [virtual]

Definition at line 54 of file MuonProphet.cc.

{}

Member Function Documentation

StatusCode MuonProphet::initialize ( ) [virtual]

Initialize geometry

Definition at line 58 of file MuonProphet.cc.

{
  info() << "Initializing MuonProphet" << endreq;
  
  size_t s=m_genToolNames.size();

  if(s != m_genYields.size() ||
     s != m_genLifetimes.size() ||
     s != m_genYieldMeasuredAt.size() )  {
    error()<<"Imcomplete parameter set"<<endreq;
    return StatusCode::FAILURE;
  }

  // print all user configuration
  info()<<"MuonProphet is active: "<<m_active<<endreq;
  m_site=Site::FromString(m_siteName.c_str());
  info()<<"MuonProphet is set to site: "<<m_siteName<<" "<<m_site<<endreq;
  info()<<"Name  ||  Yield  ||  Yield measured at  ||  Lifetime"<<endreq;
  for (size_t idx=0; idx<s; ++idx) {
    info() << m_genToolNames[idx]<< " "
           << m_genYields[idx]/(CLHEP::cm2/CLHEP::gram) <<"cm2/g "
           << m_genYieldMeasuredAt[idx]/CLHEP::GeV <<"GeV "
           << m_genLifetimes[idx]/CLHEP::second<<"s" <<endreq;
  }

  // get all spallation background gentool
  string gtn;
  for (size_t idx=0; idx<s; ++idx) {
    gtn = m_genToolNames[idx];
    try {
      m_genTools.push_back(tool<IHepMCEventMutator>(gtn));
    }
    catch(const GaudiException& exg) {
      fatal() << "Failed to get generator: \"" << gtn << "\"" << endreq;
      return StatusCode::FAILURE;
    }
    info () << "Added spallation background gentool " << gtn << endreq;
  }

  if(geometryCalculationInit().isFailure()) return StatusCode::FAILURE;
  // Spallation neutron related distribution initialization
  if(spallationNeutronsInit().isFailure()) return StatusCode::FAILURE;


  return StatusCode::SUCCESS;
}
StatusCode MuonProphet::finalize ( ) [virtual]

Definition at line 109 of file MuonProphet.cc.

{
  return StatusCode::SUCCESS;
}
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.

{
  // Not active. Don't bother to run it.
  if(!m_active) return StatusCode::SUCCESS;
  debug()<<"Processing a GenEvent"<<endreq;

  printout(event);

  m_crossRpc = false;
  m_crossWater = false;

  // Predict the fate of muon and change it.
  MpMuonFate fate = predictMuonFate(event);
  debug()<<fate<<endreq;

  if (fate.getCode() == MpMuonFate::kNeedSim) {
    // do nothing, return
    // to do:
    // But if geant4 can't produce these background, I still need to do it here.
    // For example, Li9 and He8.
    // Then it will be consistent within whole simulation scheme.
    // Right now no kNeedSim decision has been made.
  } else if (fate.getCode() == MpMuonFate::kUnknown) {
    // do nothing, return. Let geant4 to handle it.
  } else { 
    int ntrk=0;
    size_t s=m_genTools.size();

    for (size_t idx=0; idx<s; ++idx) {
      
      ntrk=spallationProducts(event, idx);
      
      if(ntrk<0) {
        error()<<"failed to run one spallation background production"<<endreq;
        return StatusCode::FAILURE;
      }
    }
    
    // spallation neutron
    // Since it will be triggered, won't be bothered with gamma, 
    // because it can only cause prompt background and it is buried in muon signal.
    if( m_actNeutron ) {
      ntrk=spallationNeutrons(event);    
      if(ntrk<0) {
        error()<<"failed to run spallation neutron background production"<<endreq;
        return StatusCode::FAILURE;
      }
    }
  }

  HepMC::GenParticle* particle = new HepMC::GenParticle(HepMC::FourVector(0.01*CLHEP::eV,
                                                                          0,
                                                                          0,
                                                                          0.01*CLHEP::eV),
                                                        22,
                                                        1/*=status*/);
  (*(event.vertices_begin()))->add_particle_out(particle);
  
  printout(event);

  return StatusCode::SUCCESS;
}
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.

{
  
  HepMC::GenEvent::particle_iterator p;
  HepMC::GenEvent::particle_iterator particleEnd=event.particles_end();
  
  int nMuon=0;
  
  m_muon=0;
  for ( p = event.particles_begin(); p!=particleEnd; ++p) {

    // muon pdg id is 13 or -13
    if( (*p)->pdg_id() == 13 || (*p)->pdg_id() == -13 ) {
      // Get the Muon pointer
      m_muon=*p;
      // In case of two or more than two muon tracks, don't work.
      ++nMuon;      
      if(nMuon>=2) {
        return MpMuonFate::kUnknown;
      }
    }
  }

  if(!m_muon) {
    info()<<"No muon track found"<<endreq;
    return MpMuonFate::kUnknown;
  }

  // calculate all needed geometry parameters
  geometryCalculation(m_muon);

  // get the trigger status of RPC and waterpool.  
  //   Note: No MpTriggerStat::kNeedSim will be issued.
  m_rpcTriggerStatus = triggeredByRpc(m_muon);
  //   Note: No MpTriggerStat::kNeedSim will be issued.
  m_poolTriggerStatus = triggeredByWaterPool(m_muon);

  debug()<<"RPC trigger status:  "<<m_rpcTriggerStatus<<endreq;
  debug()<<"Pool trigger status: "<<m_poolTriggerStatus<<endreq;

  
  MpMuonFate fate(m_rpcTriggerStatus.getCode(), m_poolTriggerStatus.getCode());
  m_muon->set_status(fate.getCode());

  return fate;

  return MpMuonFate::kUnknown;
}
StatusCode MuonProphet::geometryCalculationInit ( ) [private]

try lhcb DetDesc. This will give a realistic geometry.

RPC hit

OWS

IWS

MineralOil

rock

SmartDataPtr is not good. It will dispear when get into another function The traditional way to retrive a data is better. Stable. MixGas

OwsWater

IwsWater

Definition at line 17 of file MpGeometry.cc.

{
  /*
  double width_near = 10.0*CLHEP::meter ; // width of near pool                          
  double width_far  = 16.0*CLHEP::meter ; // width of far pool                           
  double length     = 16.0*CLHEP::meter ; // length of near and far pool                 
  double depth      = 10.0*CLHEP::meter ; // depth of near and far pool 
  double PoolDeadThickness = 84*CLHEP::mm;

  // OWS box is constructed.
  if(m_site == Site::kDayaBay || m_site == Site::kLingAo) {
    
    m_waterPool = new RectangularBox(length/2, -length/2,
                                     width_near/2, -width_near/2,
                                     depth/2, -depth/2+PoolDeadThickness);
  } else if ( m_site == Site::kFar ) {
    m_waterPool = new RectangularBox(length/2, -length/2,
                                     width_far/2, -width_far/2,
                                     depth/2, -depth/2+PoolDeadThickness);
  }
  
  info()<<"Water pool geometry "<<*m_waterPool<<endreq;
  info()<<"TrkLengthInWaterThres= "<<m_trkLengthInWaterThres/CLHEP::cm<<"cm"<<endreq;


  double height_RPC = 0.50*CLHEP::meter ; // RPCs are nominally 50cm above  surface of shield volume (includes 20cm of air above water)
  double thick_RPC  = (0.056+0.030)*CLHEP::meter ; // 1.4cm/layer*4layers + 1cm spacing*3spaces between layers
  double air_thick  =  0.2*CLHEP::meter ; // 20cm of air above water
  double extend_RPC = 1*CLHEP::meter;
  if(m_site == Site::kDayaBay || m_site == Site::kLingAo) {

    m_rpcPlane = new RpcPlane(length/2+extend_RPC, -length/2-extend_RPC,
                              width_near/2+extend_RPC, -width_near/2-extend_RPC,
                              depth/2+air_thick+height_RPC+thick_RPC/2);
  } else if ( m_site == Site::kFar ) {

    m_rpcPlane = new RpcPlane(length/2+extend_RPC, -length/2-extend_RPC,
                              width_far/2+extend_RPC, -width_far/2-extend_RPC,
                              depth/2+air_thick+height_RPC+thick_RPC/2);
  }
  info()<<"RPC geometry "<<*m_rpcPlane<<endreq;
  */

  // 1. DetectorDataSvc
  IDataProviderSvc* detSvc = 0;
  StatusCode sc = service("DetectorDataSvc",detSvc,true);
  if (sc.isFailure()) return sc;
 
  // 2. Top Detector Element
  debug()<<"DetectorElement input: "<<m_topDetectorElementName<<endreq;
  SmartDataPtr<IDetectorElement> topDE(detSvc,m_topDetectorElementName);
  if (!topDE) return StatusCode::FAILURE;
  debug()<<"DetectorElement output: "<<topDE->name()<<endreq;
 
  // 3. GeometryInfo
  m_topGeoInfo = topDE->geometry();
  if(!m_topGeoInfo) return StatusCode::FAILURE;
  debug()<<"Top Geometry "<<m_topGeoInfo->lvolumeName()<<endreq;

  // 4. get top LVolume
  const ILVolume* cILV = m_topGeoInfo->lvolume();
  m_topLVolume = const_cast<ILVolume*>(cILV);
  if(!m_topLVolume) return StatusCode::FAILURE;
  info()<<"Got top logic volume: "<<m_topLVolume->name()<<endreq;

  // 5. get the pointer to all charactoristic materials
  m_mixGas = 0;    
  m_owsWater = 0;  
  m_iwsWater = 0;  
  m_mineralOil=0;  
  m_rock = 0;      

  DataObject* pObj;

  sc = detSvc->retrieveObject("/dd/Materials/MixGas", pObj);
  if(sc.isFailure()) return sc;
  m_mixGas = dynamic_cast<Material *> (pObj);
  //SmartDataPtr<Material> m_mixGas(detSvc,"/dd/Materials/MixGas");
  if(!m_mixGas) return StatusCode::FAILURE;
  debug()<<"Got charactoristic material for RPC: "<<m_mixGas->name()<<endreq;



  sc = detSvc->retrieveObject("/dd/Materials/OwsWater", pObj);
  if(sc.isFailure()) return sc;
  m_owsWater = dynamic_cast<Material *> (pObj);
  //SmartDataPtr<Material> m_mixGas(detSvc,"/dd/Materials/OwsWater");
  if(!m_owsWater) return StatusCode::FAILURE;
  debug()<<"Got charactoristic material for OWS: "<<m_owsWater->name()<<endreq;

  sc = detSvc->retrieveObject("/dd/Materials/IwsWater", pObj);
  if(sc.isFailure()) return sc;
  m_iwsWater = dynamic_cast<Material *> (pObj);
  //SmartDataPtr<Material> m_mixGas(detSvc,"/dd/Materials/IwsWater");
  if(!m_iwsWater) return StatusCode::FAILURE;
  debug()<<"Got charactoristic material for IWS: "<<m_iwsWater->name()<<endreq;

  // Rock
  sc = detSvc->retrieveObject("/dd/Materials/Rock", pObj);
  if(sc.isFailure()) return sc;
  m_rock = dynamic_cast<Material *> (pObj);
  //SmartDataPtr<Material> m_rock(detSvc,"/dd/Materials/Rock");
  if(!m_rock) return StatusCode::FAILURE;
  debug()<<"Got charactoristic material for rock: "<<m_rock->name()<<endreq;

  // MineralOil
  sc = detSvc->retrieveObject("/dd/Materials/MineralOil", pObj);
  if(sc.isFailure()) return sc;
  m_mineralOil = dynamic_cast<Material *> (pObj);
  //SmartDataPtr<Material> m_rock(detSvc,"/dd/Materials/MineralOil");
  if(!m_mineralOil) return StatusCode::FAILURE;
  debug()<<"Got charactoristic material for AD: "<<m_mineralOil->name()<<endreq;

  //  debug()<<*m_mixGas<<endreq;
  //  debug()<<int(m_mixGas)<<" "<<int(m_owsWater)<<" "<<int(m_rock)<<endreq;

  vector<string> ADNames;
  
  if(m_site == Site::kDayaBay) 
    {
      ADNames.push_back("/dd/Structure/AD/db-oil1");
      ADNames.push_back("/dd/Structure/AD/db-oil2");
    } 
  if(m_site == Site::kLingAo)
    {
      ADNames.push_back("/dd/Structure/AD/la-oil1");
      ADNames.push_back("/dd/Structure/AD/la-oil2");
    }
  if(m_site == Site::kFar)  
    {
      ADNames.push_back("/dd/Structure/AD/far-oil1");
      ADNames.push_back("/dd/Structure/AD/far-oil2");
      ADNames.push_back("/dd/Structure/AD/far-oil3");
      ADNames.push_back("/dd/Structure/AD/far-oil4");
    }

  IGeometryInfo* pGI;
  for(unsigned int idx = 0; idx < ADNames.size(); ++idx)
    {
      // Get Detector Element
      SmartDataPtr<IDetectorElement> pDE(detSvc,ADNames[idx]);
      if (!pDE) return StatusCode::FAILURE;
      debug()<<"Got DetectorElement: "<<pDE->name()<<endreq;
 
      // Then geometryInfo
      pGI = pDE->geometry();
      if(!pGI) return StatusCode::FAILURE;
      debug()<<"Got GeometryInfo: "<<pGI->lvolumeName()<<endreq;
      
      m_ADs.push_back(pGI);
    }

  return StatusCode::SUCCESS;

}
StatusCode MuonProphet::geometryCalculation ( HepMC::GenParticle pMuon) [private]

Old fast simplified geometry. However not precise.

look for intersections with RPC

look for intersections with water pool

For any positive track length in water, crossWater = true

Regulate intersections. Only retain information with rock, water and oil. all volumes within oil are called oil iwsWater and owsWater are combined.

If nothing added, so create a new one.

combine the volume with the same material

Definition at line 184 of file MpGeometry.cc.

{
  /*
  HepGeom::Point3D<double>  aPoint(pMuon->production_vertex()->position().x(),
                                   pMuon->production_vertex()->position().y(),
                                   pMuon->production_vertex()->position().z());
  HepGeom::Point3D<double> aSlope(pMuon->momentum().px(),
                                  pMuon->momentum().py(),
                                  pMuon->momentum().pz());
 
  aSlope = aSlope.unit();

  m_crossWater = m_waterPool->intersect(aPoint,aSlope,m_crossWaterA,m_crossWaterB);
  
  if(m_crossWater) {
    debug() << "Intersect with water OWS at " << m_crossWaterA <<" and "<< m_crossWaterB <<endreq;
    
    m_trkLengthInWater = m_crossWaterA.distance(m_crossWaterB);

    debug() << "Track length in water is "<<m_trkLengthInWater/CLHEP::cm<<"cm"<<endreq;
  }

  
  m_crossRpc = m_rpcPlane->intersect(aPoint,aSlope,m_rpcIntersect);
  
  if(m_crossRpc) {
    debug()<<"Intersect with RPC at "<<m_rpcIntersect <<endreq;
  }
  */

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

  m_point.SetXYZ(pMuon->production_vertex()->position().x(),
                 pMuon->production_vertex()->position().y(),
                 pMuon->production_vertex()->position().z());
  m_vec.SetXYZ(pMuon->momentum().px(),
               pMuon->momentum().py(),
               pMuon->momentum().pz());
  m_unit = m_vec.Unit();

  //Gaudi::XYZPoint point(0,0,0);
  debug()<<"Point  "<<m_point<<endreq;
  string path = m_topGeoInfo->belongsToPath(m_point,-1);
  debug()<<"Belong to path "<<path<<endreq;
  
  //Gaudi::XYZPoint point(0,0,0);
  //Gaudi::XYZVector vec(1,0,0);

  unsigned int N = m_topLVolume->intersectLine(m_point, 
                                               m_unit, 
                                               m_intersections, 
                                               0,100*CLHEP::meter,
                                               -1);
  
  debug()<<"Found "<<N<<" intersections"<<endreq;
  unsigned int idx;
  for (idx=0; idx<N; ++idx) {
    debug()<<"intersection: "<<m_intersections[idx].first<<" "      
          <<m_intersections[idx].second->name()<<"  "
          <<m_intersections[idx].second->density()<<endreq;
    //<<" "<<int(m_intersections[idx].second)<<endreq;
  }

  // Stop muon. Stopping power = 2 [MeVcm2/g]
  // Begin with an approximation
  double ke = pMuon->momentum().e() - pMuon->momentum().m();
  double sp = 2 * CLHEP::MeV * CLHEP::cm2 / CLHEP::g;
  double trkLength = ke/sp;
  //debug()<< CLHEP::MeV << "  " << CLHEP::cm2 << "  " << CLHEP::g <<endreq;
  debug()<<"ke= "<<ke<<"  sp="<<sp<<"  trkLength="<<trkLength<<endreq;
  
  for (idx=0; idx<N; ++idx) {
    double Lpred = m_intersections[idx].first.second - m_intersections[idx].first.first;
    debug()<<"Length in this volume "<< Lpred*m_intersections[idx].second->density()<<endreq;
    if( Lpred*m_intersections[idx].second->density() < trkLength )  {
      trkLength -= Lpred*m_intersections[idx].second->density();
    } else { // end of track
      m_intersections[idx].first.second = 
        m_intersections[idx].first.first + 
        trkLength / m_intersections[idx].second->density();
      break;
    }
  }
  
  if( idx<N-1 )  {  // erase [ first, last ) 
    m_intersections.erase( m_intersections.begin()+idx+1, m_intersections.end() );
  }
  
  N = m_intersections.size();
  debug()<<"Left "<<N<<" intersections"<<endreq;
  for (idx=0; idx<N; ++idx) {
    debug()<<"intersection: "<<m_intersections[idx].first<<" "
          <<m_intersections[idx].second->name()<<"  "
          <<m_intersections[idx].second->density()<<endreq;
    //<<" "<<int(m_intersections[idx].second)<<endreq;
  }



  for (idx=0; idx<N; ++idx) {
    if(m_intersections[idx].second == m_mixGas) {
      m_crossRpc = true;
      break;
    }
  }
  
  bool first = true;
  m_tickWaterMin = 0;
  m_tickWaterMax = 0;

  for (idx=0; idx<N; ++idx) {
    if(m_intersections[idx].second == m_owsWater ||
       m_intersections[idx].second == m_iwsWater ) {
      // first time find OwsWater
      if(first) {
        m_tickWaterMin = m_intersections[idx].first.first;
        m_tickWaterMax = m_intersections[idx].first.second;
        first = false;
      } else {
        if(m_intersections[idx].first.first < m_tickWaterMin ) {
          m_tickWaterMin = m_intersections[idx].first.first;
        }
        if(m_intersections[idx].first.second > m_tickWaterMax ) {
          m_tickWaterMax = m_intersections[idx].first.second;
        }
      }
    }
  }

   
  m_trkLengthInWater = m_tickWaterMax - m_tickWaterMin;
  debug()<<"Track length in water "<< m_trkLengthInWater/CLHEP::m <<"m"<<endreq;

  if( m_trkLengthInWater>0 ) m_crossWater = true;

  m_crucialSegments.clear();
  
  for (idx=0; idx<N; ++idx) {
    if(m_intersections[idx].second == m_rock ||
       m_intersections[idx].second == m_owsWater ||
       m_intersections[idx].second == m_iwsWater ||
       m_intersections[idx].second == m_mineralOil ) 
      {
        if(m_crucialSegments.empty())    
          {
            ILVolume::Intersection * seg = new ILVolume::Intersection(m_intersections[idx]);
            m_crucialSegments.push_back(*seg);
          } else {   
            // check whether last one has the same material as current volume
            // iwsWater and owsWater are considered as the same
            if((m_crucialSegments.back().second == m_intersections[idx].second) ||
               (m_intersections[idx].second == m_iwsWater && m_crucialSegments.back().second == m_owsWater) ||
               (m_intersections[idx].second == m_owsWater && m_crucialSegments.back().second == m_iwsWater) )
              { // If yes, combine them. Assume LS and oil are same
                m_crucialSegments.back().first.second = m_intersections[idx].first.second;
              } else {
                ILVolume::Intersection * seg = new ILVolume::Intersection(m_intersections[idx]);
                m_crucialSegments.push_back(*seg);
              }
          }
      }
  }

  unsigned int NJiaodian = m_crucialSegments.size();
  for (idx=0; idx<NJiaodian; ++idx) {
    debug()<<"intersection: "<<m_crucialSegments[idx].first<<" "
           <<m_crucialSegments[idx].second->name()<<endreq;
    //           <<" "<<int(m_crucialSegments[idx].second)<<endreq;
  }

  return StatusCode::SUCCESS;
  
}
MpTriggerStat MuonProphet::triggeredByRpc ( HepMC::GenParticle pMuon) [private]

new one 3 Dec. 2009. Never return NeedSim

Definition at line 25 of file MpTrigger.cc.

{
  /* old one
  if(m_crossRpc) {
    return MpTriggerStat::kTriggered;
  } else {
    return MpTriggerStat::kNeedSim;
  }
  */
  
  if(m_crossRpc) {
    return MpTriggerStat::kTriggered;
  } else {
    return MpTriggerStat::kFarAway;
  }
  //
}
MpTriggerStat MuonProphet::triggeredByWaterPool ( HepMC::GenParticle pMuon) [private]

new one 3 Dec. 2009. Never return NeedSim

Definition at line 45 of file MpTrigger.cc.

{
  /* old one
  if(m_crossWater) {
    if(m_trkLengthInWater > m_trkLengthInWaterThres) {
      return MpTriggerStat::kTriggered;
    }
  }

  return MpTriggerStat::kNeedSim;
  */

  if(m_crossWater) {
    if(m_trkLengthInWater > m_trkLengthInWaterThres) {    // 100% triggered
      return MpTriggerStat::kTriggered;
    } else {
      if(m_rnd.Rndm()<0.5) {   // 50% - 50%  triggered or inefficent
        return MpTriggerStat::kTriggered;
      } else {
        return MpTriggerStat::kIneffi;
      }
    } 
  } else {
    if(m_rnd.Rndm()<0.05) {   // 5% triggered  95% faraway
      return MpTriggerStat::kTriggered;
    } else {
      return MpTriggerStat::kFarAway;
    }
  }
}
int MuonProphet::spallationProducts ( HepMC::GenEvent event,
int  idx 
) [private]

$$$$$ The idx, index, is for parameter seaching in m_gen...

For a given muon track, add spallation products, return number of tracks added.

vector series. Return number of tracks attached, negative for error

The idx, index, is for parameter seaching in m_gen... vector series.

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.

{
  if(m_muon == 0) {
    error()<<"No muon track found"<<endreq;
    // Since no muon track, it should not reach this point.
    return -1;
  }
  
  HepMC::ThreeVector bkgPoint;
  bool bkgGenerated = false;
  bkgGenerated = generateOneBkg(idx, 
                                bkgPoint);

  if( !bkgGenerated ) return 0;

  HepMC::GenEvent * bkgEvt=new HepMC::GenEvent;
  if(m_genTools[idx]->mutate(*bkgEvt).isFailure()) {
    error()<<"Failed to run background gentool "<<m_genToolNames[idx]<<endreq;
    return -1;
  }
  
  // debug
#if 0
  {
    HepMC::GenEvent::particle_iterator p;
    HepMC::GenEvent::particle_iterator particleEnd = bkgEvt->particles_end();
    for ( p = bkgEvt->particles_begin(); p!=particleEnd; ++p) {
      debug()<<*(*p)<<endreq;
    }
    HepMC::GenEvent::vertex_iterator v;
    HepMC::GenEvent::vertex_iterator vertexEnd = bkgEvt->vertices_end();
    for ( v = bkgEvt->vertices_begin(); v!=vertexEnd; ++v) {
      debug()<<*(*v)<<endreq;
    }
  }
#endif

  // change the vertex position first
  setPosition(*bkgEvt,bkgPoint);

  // change the time of the vertex
  // Get muon vertex time 
  double t0 = m_muon->production_vertex()->position().t();
  // to do: A delta should be added to t0 to account for flight time of muon
  //        and flight time of muon's daughter particle which reached the 
  //        production vertex of the background.
  
  setTime(*bkgEvt,t0,m_genLifetimes[idx]);

  // drag vertices out from background event 
  HepMC::GenVertex *bkgVtx;
  while( !bkgEvt->vertices_empty() ){

    bkgVtx = *(bkgEvt->vertices_begin());

    // add to muon event
    // This will automatically transfer the ownership from the old event to new event.
    // So the original event can be safely deleted now.
    event.add_vertex(bkgVtx);
  }

  // delete the bkgEvt
  delete bkgEvt;

  return 0;
}
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 generated it will return the position which should be in one AD.

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.

{
  if( !m_crossWater ) {
    return false;
  }

  double eMuon = m_muon->momentum().e();

  // Mineral oil and LS densities are assumed to be the same.
  double density = m_mineralOil->density();
  double yield = m_genYields[idx] ;
  double alpha = 0.77;

  double yieldAtThisEnergy = yield * pow( (eMuon / m_genYieldMeasuredAt[idx]), alpha);

  double ProbThres = yieldAtThisEnergy * density * m_trkLengthInWater;

  debug()<<m_genToolNames[idx]<<" eMuon= "<<eMuon/CLHEP::GeV<<"GeV "
         <<" trkLengthInWater= "<<m_trkLengthInWater/CLHEP::meter<<"m "
         <<"ProbThres= "<<ProbThres<<endreq;

  if( ProbThres > 1 || ProbThres < 0 ) {
    debug()<<"Yield too high or too low. Invalid probability "<<ProbThres<<endreq;
  }

  if( m_rnd.Rndm() > ProbThres) {
    return false;
  }

  
  // The position of background in the cylinder coordinate system of muon track
  // Muon direction is z.
  double r,phi,z;

  z = m_rnd.Rndm() * (m_tickWaterMax - m_tickWaterMin) + m_tickWaterMin;
  r = m_nLateralPDF->GetRandom() * CLHEP::cm ;
  phi = m_rnd.Rndm() * 2 * 3.14159265358979323846;


  // try point in coordinate system of intersect 1st
  HepGeom::Point3D<double> tryPoint;

  tryPoint.setX( r*cos(phi) );
  tryPoint.setY( r*sin(phi) );
  tryPoint.setZ( z );

  /* This is a muon track shown below.
   *
   *            intersect 1st (tick1)  intersect 2nd (tick2)
   * Primary vertex -----1--------C----------2------>
   *                        z     | theta
   *                              |  
   *                              | r
   *                              |
   *                              |
   *                           tryPoint
   *
   */

  // direction of muon
  HepGeom::Point3D<double> muonDir(m_unit.X(),
                                   m_unit.Y(),
                                   m_unit.Z());
 
  // need a y axis which is perpendicular to muonSlope
  HepGeom::Point3D<double> y( m_unit.X(),
                              m_unit.Y(),
                              0);
  y.rotateZ(3.14159265358979323846/2);
  if(y.r()<=0) y.setX(1);

  // tryPoint in ows coordinate system
  HepGeom::Point3D<double> tryPointGlobal = muonDir * tryPoint.r();

  // rotation
  tryPointGlobal.rotate(tryPoint.theta(), y);
  tryPointGlobal.rotate(tryPoint.phi(),   muonDir);

  // transform                                                                                                              
  tryPointGlobal.setX(tryPointGlobal.x() + m_point.X());
  tryPointGlobal.setY(tryPointGlobal.y() + m_point.Y());
  tryPointGlobal.setZ(tryPointGlobal.z() + m_point.Z());

  debug()<<tryPoint<<endreq;
  debug()<<tryPointGlobal<<endreq;

  Gaudi::XYZPoint p(tryPointGlobal.x(),
                    tryPointGlobal.y(),
                    tryPointGlobal.z());
  
  bool inside = false;
  for (unsigned int idx = 0; idx<m_ADs.size(); ++idx) 
    {
      if(m_ADs[idx]->isInside(p)) inside = true;
    }
  
  if(!inside) return false;
  
  // done
  bkgPoint.setX(tryPointGlobal.x());
  bkgPoint.setY(tryPointGlobal.y());
  bkgPoint.setZ(tryPointGlobal.z());

  return true;
}
StatusCode MuonProphet::spallationNeutronsInit ( ) [private]

Neutron production Neutron production will run on its own set of parameters' value.

Spallation neutron This package runs on all kinds of distributions in ../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.

{
  m_nEnergyPDF=0;
  m_nAnglePDF=0;
  m_nLateralPDF=0;


  // Avarage muon energy is from the estimation of Mengyun Guan. (Ph.D. thesis)
  double avaMuonE;
  if(m_site == Site::kDayaBay) {
    avaMuonE = 55.3*CLHEP::GeV;
  } else if (m_site == Site::kLingAo) {
    avaMuonE = 61.4*CLHEP::GeV;
  } else if (m_site == Site::kFar ) {
    avaMuonE = 140.3*CLHEP::GeV;
  }
  
  m_nEnergyPDF = new TF1("m_nEnergyPDF",NeutronEnergyPDF,0.001,4,1);
  m_nEnergyPDF->SetNpx(200);  

  m_nEnergyPDF->SetParameter(0,avaMuonE/CLHEP::GeV);
  if(!m_nEnergyPDF) {
    error()<<"Failed to get NeutronEnergyPDF"<<endreq;
    return StatusCode::FAILURE;
  }

  m_nAnglePDF = new TF1("m_nAnglePDF",NeutronAngularDis,-1,1,1);
  m_nAnglePDF->SetParameter(0,avaMuonE/CLHEP::GeV);
  if(!m_nAnglePDF) {
    error()<<"Failed to get NeutronAngularDis"<<endreq;
    return StatusCode::FAILURE;
  }

  
  m_nLateralPDF = new TF1("m_nLateralPDF",LateralProfile,0,100,1);
  m_nLateralPDF->SetParameter(0,0);
  if(!m_nLateralPDF) {
    error()<<"Failed to get LateralProfile"<<endreq;
    return StatusCode::FAILURE;
  }

  return StatusCode::SUCCESS;
}
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.

{
  debug()<<"Spallation Neutron production"<<endreq;

  //--------------------------------------------------------
  unsigned int NSeg = m_crucialSegments.size();
  
  double yield;
  if(m_neutronYield>=0) {
    yield = m_neutronYield;
  } else {
    yield = NeutronYield( m_muon->momentum().e()/CLHEP::GeV )*(CLHEP::cm2/CLHEP::gram);
  }

  double tracklength;
  double multiplicity;

  int nNewTrack = 0;

  HepMC::GenEvent * bkgEvt=new HepMC::GenEvent;

  for (unsigned int idx=0; idx<NSeg; ++idx) 
    {
      double prob;

      
      tracklength = m_crucialSegments[idx].first.second - m_crucialSegments[idx].first.first;
      multiplicity = NeutronMulti( m_muon->momentum().e()/CLHEP::GeV,
                                   tracklength );
      
      prob = yield
        * tracklength
        * m_crucialSegments[idx].second->density()
        / multiplicity;
      
      debug()<<"Yield= "<<NeutronYield( m_muon->momentum().e()/CLHEP::GeV )<<"cm2/g "
             <<"TrackLength= "<< tracklength/CLHEP::m <<"m "
             <<"Density= "<<m_crucialSegments[idx].second->density()/(CLHEP::g/CLHEP::cm3)<<"g/cm3 "
             <<"Multiplicity= "<< multiplicity <<endreq;
      debug()<<"prob= "<<prob<<endreq;

      

      if( m_rnd.Rndm() < prob ) 
        {
          nNewTrack += generateNeutrons(bkgEvt,
                                        m_crucialSegments[idx].first.first,   // tick one
                                        m_crucialSegments[idx].first.second,  // tick two
                                        multiplicity);
        }
    }  
  
  // drag vertices out from background event

  HepMC::GenVertex *bkgVtx;
  while( !bkgEvt->vertices_empty() )
    {
    bkgVtx = *(bkgEvt->vertices_begin());
    
    // add to muon event
    // This will automatically transfer the ownership from the old event to new event.
    // So the original event can be safely deleted now.
    event.add_vertex(bkgVtx);
    }
 
  // delete the bkgEvt
  delete bkgEvt;

  return nNewTrack;
  
}
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.

{
  int nNeutron = m_rnd.Poisson(multiplicity);
  
  debug()<<"nNeutron= "<<nNeutron<<endreq;
  

  double z;
  double r;
  double phi;
  
  double ke;
  double eng;
  double momentum;
  double massN = 939.565346*CLHEP::MeV;
  double phiN;
  double costh;
  double sinth;

  for(int idx =0; idx<nNeutron; ++idx)  
    { 
      z = m_rnd.Rndm() * (tickEnd - tickBegin) + tickBegin;
      r = m_nLateralPDF->GetRandom() * CLHEP::cm ;
      phi = m_rnd.Rndm() * 2 * 3.14159265358979323846;

      ke = m_nEnergyPDF->GetRandom() * CLHEP::GeV ;      
      eng = ke + massN;
      momentum = sqrt (eng*eng - massN*massN);
      phiN = m_rnd.Rndm() * 2 * 3.14159265358979323846;
      costh = m_nAnglePDF->GetRandom();
      sinth = sqrt(1-costh*costh);


      debug()<<"z "<<z/CLHEP::m<<"m"<<endreq;
      debug()<<"r "<<r/CLHEP::cm<<"cm"<<endreq;
      debug()<<"ke "<<ke/CLHEP::GeV<<"GeV"<<endreq;
      debug()<<"costh "<<costh<<endreq;

      HepGeom::Point3D<double> genPoint;
      
      genPoint.setX( r*cos(phi) );
      genPoint.setY( r*sin(phi) );
      genPoint.setZ( z );

      HepGeom::Point3D<double> genVec;
      
      genVec.setX ( momentum*sinth*cos(phiN) );
      genVec.setY ( momentum*sinth*sin(phiN) );
      genVec.setZ ( momentum*costh );

      HepGeom::Point3D<double> muonDir(m_unit.X(),
                                       m_unit.Y(),
                                       m_unit.Z());

      // need a y axis which is perpendicular to muonSlope
      HepGeom::Point3D<double> y( m_unit.X(),
                                  m_unit.Y(),
                                  0);
      y.rotateZ(3.14159265358979323846/2);
      if(y.r()<=0) y.setX(1);
      
      HepGeom::Point3D<double> genPointGlobal = muonDir * genPoint.r();
      debug()<<m_unit<<endreq;
      debug()<<muonDir<<endreq;
      debug()<<genPointGlobal<<endreq;

      // rotation
      genPointGlobal.rotate(genPoint.theta(), y);
      genPointGlobal.rotate(genPoint.phi(),   muonDir);

      // transform
      genPointGlobal.setX(genPointGlobal.x() + m_point.X());
      genPointGlobal.setY(genPointGlobal.y() + m_point.Y());
      genPointGlobal.setZ(genPointGlobal.z() + m_point.Z());

      HepGeom::Point3D<double> genVectorGlobal = muonDir * genVec.r();
      
      genVectorGlobal.rotate(genVec.theta(),  y);
      genVectorGlobal.rotate(genVec.phi(),    muonDir);
      

      double t0 = m_muon->production_vertex()->position().t();
      debug()<<"t0= "<<t0<<endreq;
      t0 += tickEnd / (3e8 * CLHEP::m/CLHEP::s);
      debug()<<"dt= "<<tickEnd / (3e8 * CLHEP::m/CLHEP::s)<<endreq;
      debug()<<"t0= "<<t0<<endreq;
      HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(genPointGlobal.x(),
                                                                        genPointGlobal.y(),
                                                                        genPointGlobal.z(),
                                                                        t0));
      
      HepMC::GenParticle* particle = new HepMC::GenParticle(HepMC::FourVector(genVectorGlobal.x(),
                                                                              genVectorGlobal.y(),
                                                                              genVectorGlobal.z(),
                                                                              eng),
                                                                              2112,
                                                                              1/*=status*/);
      vertex->add_particle_out(particle);
      bkgEvt->add_vertex(vertex);
    } 

  return nNeutron;
}
StatusCode MuonProphet::printout ( HepMC::GenEvent event) [private]

Definition at line 240 of file MuonProphet.cc.

{
  // every vertex
  HepMC::GenEvent::vertex_iterator vtx, done = event.vertices_end();
  for (vtx = event.vertices_begin(); vtx != done; ++vtx) {
    debug()<<*(*vtx)<<endreq;
    // every particle
    HepMC::GenVertex::particles_in_const_iterator p;
    HepMC::GenVertex::particles_in_const_iterator inEnd=(*vtx)->particles_in_const_end();
    for ( p = (*vtx)->particles_in_const_begin(); p!=inEnd; ++p) {
      debug()<<*(*p)<<endreq;
    }
    HepMC::GenVertex::particles_out_const_iterator outEnd=(*vtx)->particles_out_const_end();
    for ( p = (*vtx)->particles_out_const_begin(); p!=outEnd; ++p) {
      debug()<<*(*p)<<endreq;
    }
  }
  return StatusCode::SUCCESS;
}
void MuonProphet::setPosition ( HepMC::GenEvent event,
HepMC::ThreeVector position 
) [private]

..........................................

Definition at line 13 of file MpPlace.cc.

{
  HepMC::FourVector position(aPoint.x(),
                             aPoint.y(),
                             aPoint.z(),
                             0);

  HepMC::GenEvent::vertex_iterator vtx, done = event.vertices_end();
  for (vtx = event.vertices_begin(); vtx != done; ++vtx) {

    // still use the original vertex time
    position.setT( (*vtx)->position().t() );

    (*vtx)->set_position(position);
  }

}
void MuonProphet::setTime ( HepMC::GenEvent event,
double  t0,
double  lifetime 
) [private]

The unit in geant4 is nano second. Very difficult to figure out.

Definition at line 33 of file MpPlace.cc.

{
  // Use negtive lifetime to skip the exponential random lifetime

  //
  // debug()<<"t0= "<<t0<<endreq;
  //
  double dt = 0;

  if(lifetime<0) {
    dt=0;
  } else {
    double u = m_rnd.Rndm();
    // Exponential Distribution: t(n) = t(n-1) * exp(-t/tL)
    dt = ((-1.0 * log(u)) * lifetime);
  }

  dt = dt / CLHEP::nanosecond;

  HepMC::GenEvent::vertex_iterator vtx, done = event.vertices_end();
  for (vtx = event.vertices_begin(); vtx != done; ++vtx) {

    HepMC::FourVector position = (*vtx)->position();

    // relative to muon + mother particle's lifetime + internal delay
    double vertex_time = t0+dt+position.t();

    position.setT(vertex_time);
    (*vtx)->set_position(position);
  }
}

Member Data Documentation

bool MuonProphet::m_active [private]

Definition at line 114 of file MuonProphet.h.

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.

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.

Definition at line 164 of file MuonProphet.h.

TRandom3 MuonProphet::m_rnd [private]

random number

Definition at line 167 of file MuonProphet.h.

A simplified water pool geometry.

Definition at line 170 of file MuonProphet.h.

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]

MineralOil.

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.

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.

Definition at line 223 of file MuonProphet.h.

Definition at line 232 of file MuonProphet.h.

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.

Definition at line 242 of file MuonProphet.h.

TF1* MuonProphet::m_nEnergyPDF [private]

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

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.

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:
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:21:29 for MuonProphet by doxygen 1.7.4