MuonProphet.cc

Go to the documentation of this file.

#include "MuonProphet.h"
// 
MuonProphet::MuonProphet(const std::string& type,
                         const std::string& name,
                         const IInterface* parent)
  : 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()
{}

  
StatusCode MuonProphet::initialize()
{
  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()
{
  return StatusCode::SUCCESS;
}
  
// HepMCEventMutator interface
StatusCode MuonProphet::mutate(HepMC::GenEvent& event)
{
  // 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;
}
  

// Predict the fate of a muon, alter the status of the muon track
// This will affect all the following simulation behaviour, FAST or FULL?
MpMuonFate MuonProphet::predictMuonFate(HepMC::GenEvent& event)
{
  
  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::printout(HepMC::GenEvent& event)
{
  // 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;
}