evgen::NeutronOsc Class Reference

Inheritance diagram for evgen::NeutronOsc:

Public Member Functions
	NeutronOsc (fhicl::ParameterSet const &p)
	NeutronOsc (NeutronOsc const &)=delete
	NeutronOsc (NeutronOsc &&)=delete
NeutronOsc &	operator= (NeutronOsc const &)=delete
NeutronOsc &	operator= (NeutronOsc &&)=delete
void	produce (art::Event &e) override
void	beginRun (art::Run &run) override
Public Member Functions inherited from art::EDProducer
	EDProducer (fhicl::ParameterSet const &pset)
template<typename Config >
	EDProducer (Table< Config > const &config)
std::string	workerType () const
Public Member Functions inherited from art::detail::Producer
virtual	~Producer () noexcept
	Producer (fhicl::ParameterSet const &)
	Producer (Producer const &)=delete
	Producer (Producer &&)=delete
Producer &	operator= (Producer const &)=delete
Producer &	operator= (Producer &&)=delete
void	doBeginJob (SharedResources const &resources)
void	doEndJob ()
void	doRespondToOpenInputFile (FileBlock const &fb)
void	doRespondToCloseInputFile (FileBlock const &fb)
void	doRespondToOpenOutputFiles (FileBlock const &fb)
void	doRespondToCloseOutputFiles (FileBlock const &fb)
bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)
bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
Public Member Functions inherited from art::Modifier
	~Modifier () noexcept
	Modifier ()
	Modifier (Modifier const &)=delete
	Modifier (Modifier &&)=delete
Modifier &	operator= (Modifier const &)=delete
Modifier &	operator= (Modifier &&)=delete
Public Member Functions inherited from art::ModuleBase
virtual	~ModuleBase () noexcept
	ModuleBase ()
ModuleDescription const &	moduleDescription () const
void	setModuleDescription (ModuleDescription const &)
std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const
void	sortConsumables (std::string const &current_process_name)
template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)
template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

Private Member Functions
int	SelectAnnihilationMode (int pdg_code)

Private Attributes
const genie::EventRecordVisitorI *	mcgen
genie::NNBarOscMode_t	gOptDecayMode = genie::kNONull
CLHEP::RandFlat	flatDist

Additional Inherited Members
Public Types inherited from art::EDProducer
using	ModuleType = EDProducer
using	WorkerType = WorkerT< EDProducer >
Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using	Table = Modifier::Table< UserConfig, KeysToIgnore >
Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using	Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >
Static Public Member Functions inherited from art::EDProducer
static void	commitEvent (EventPrincipal &ep, Event &e)
Protected Member Functions inherited from art::ModuleBase
ConsumesCollector &	consumesCollector ()
template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)
template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)
template<typename T , BranchType = InEvent>
void	consumesMany ()
template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)
template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)
template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

Detailed Description

Definition at line 71 of file NeutronOsc_module.cc.

Constructor & Destructor Documentation

evgen::NeutronOsc::NeutronOsc ( fhicl::ParameterSet const &  p )

explicit

Definition at line 102 of file NeutronOsc_module.cc.

  : art::EDProducer{p}
  // create a default random engine; obtain the random seed from NuRandomService,
  // unless overridden in configuration with key "Seed"
  , flatDist{art::ServiceHandle<rndm::NuRandomService>{}->createEngine(*this, p, "Seed")}
{
  genie::PDGLibrary::Instance(); //Ensure Messenger is started first in GENIE.

  string sname   = "genie::EventGenerator";
  // GENIE v2 // string sconfig = "NeutronOsc";
  string sconfig = "NNBarOsc";
  // GENIE v3 needs a tune (even if irrelevant)
  evgb::SetEventGeneratorListAndTune("Default","Default");

  genie::AlgFactory * algf = genie::AlgFactory::Instance();
  mcgen =
    dynamic_cast<const genie::EventRecordVisitorI *> (algf->GetAlgorithm(sname,sconfig));
  if(!mcgen) {
    throw cet::exception("NeutronOsc") << "Couldn't instantiate the neutron-antineutron oscillation generator";
  }
  int fDecayMode = p.get<int>("DecayMode");
  gOptDecayMode = (genie::NNBarOscMode_t) fDecayMode;

  produces< std::vector<simb::MCTruth> >();
  produces< sumdata::RunData, art::InRun >();

  unsigned int seed = art::ServiceHandle<rndm::NuRandomService>()->getSeed();
  genie::utils::app_init::RandGen(seed);
}

evgen::NeutronOsc::NeutronOsc ( NeutronOsc const &  )

delete

evgen::NeutronOsc::NeutronOsc ( NeutronOsc &&  )

delete

Member Function Documentation

void evgen::NeutronOsc::beginRun ( art::Run &  run )

overridevirtual

Reimplemented from art::EDProducer.

Definition at line 216 of file NeutronOsc_module.cc.

{
  art::ServiceHandle<geo::Geometry const> geo;
  run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()));
}

NeutronOsc& evgen::NeutronOsc::operator= ( NeutronOsc const &  )

delete

NeutronOsc& evgen::NeutronOsc::operator= ( NeutronOsc &&  )

delete

void evgen::NeutronOsc::produce ( art::Event &  e )

overridevirtual

Implements art::EDProducer.

Definition at line 132 of file NeutronOsc_module.cc.

{
  // Implementation of required member function here.
  genie::EventRecord * event = new genie::EventRecord;
  int target = 1000180400;  //Only use argon target
  int decay  = SelectAnnihilationMode(target);
  genie::Interaction * interaction = genie::Interaction::NOsc(target,decay);
  event->AttachSummary(interaction);

  // Simulate decay
  mcgen->ProcessEventRecord(event);

//  genie::Interaction *inter = event->Summary();
//  const genie::InitialState &initState = inter->InitState();
//  std::cout<<"initState = "<<initState.AsString()<<std::endl;
//  const genie::ProcessInfo &procInfo = inter->ProcInfo();
//  std::cout<<"procInfo = "<<procInfo.AsString()<<std::endl;
  MF_LOG_DEBUG("NeutronOsc")
    << "Generated event: " << *event;

  std::unique_ptr< std::vector<simb::MCTruth> > truthcol(new std::vector<simb::MCTruth>);
  simb::MCTruth truth;

  art::ServiceHandle<geo::Geometry const> geo;

  // Find boundary of active volume
  double minx = 1e9;
  double maxx = -1e9;
  double miny = 1e9;
  double maxy = -1e9;
  double minz = 1e9;
  double maxz = -1e9;
  for (size_t i = 0; i<geo->NTPC(); ++i){
    const geo::TPCGeo &tpc = geo->TPC(i);
    if (minx>tpc.MinX()) minx = tpc.MinX();
    if (maxx<tpc.MaxX()) maxx = tpc.MaxX();
    if (miny>tpc.MinY()) miny = tpc.MinY();
    if (maxy<tpc.MaxY()) maxy = tpc.MaxY();
    if (minz>tpc.MinZ()) minz = tpc.MinZ();
    if (maxz<tpc.MaxZ()) maxz = tpc.MaxZ();
  }

  // Assign vertice position
  double X0 = flatDist.fire( minx, maxx );
  double Y0 = flatDist.fire( miny, maxy );
  double Z0 = flatDist.fire( minz, maxz );

  TIter partitr(event);
  genie::GHepParticle *part = 0;
  // GHepParticles return units of GeV/c for p.  the V_i are all in fermis
  // and are relative to the center of the struck nucleus.
  // add the vertex X/Y/Z to the V_i for status codes 0 and 1
  int trackid = 0;
  std::string primary("primary");

  while( (part = dynamic_cast<genie::GHepParticle *>(partitr.Next())) ){

    simb::MCParticle tpart(trackid,
                           part->Pdg(),
                           primary,
                           part->FirstMother(),
                           part->Mass(),
                           part->Status());

    TLorentzVector pos(X0, Y0, Z0, 0);
    TLorentzVector mom(part->Px(), part->Py(), part->Pz(), part->E());
    tpart.AddTrajectoryPoint(pos,mom);
    if(part->PolzIsSet()) {
      TVector3 polz;
      part->GetPolarization(polz);
      tpart.SetPolarization(polz);
    }
    truth.Add(tpart);

    ++trackid;
  }// end loop to convert GHepParticles to MCParticles
  truth.SetOrigin(simb::kUnknown);
  truthcol->push_back(truth);
  //FillHistograms(truth);
  e.put(std::move(truthcol));

  delete event;
}

int evgen::NeutronOsc::SelectAnnihilationMode ( int  pdg_code )

private

Definition at line 223 of file NeutronOsc_module.cc.

{
  // if the mode is set to 'random' (the default), pick one at random!
  if ((int)gOptDecayMode == 0) {
    int mode;

    std::string pdg_string = std::to_string(static_cast<long long>(pdg_code));
    if (pdg_string.size() != 10) {
      std::cout << "Expecting PDG code to be a 10-digit integer; instead, it's the following: " << pdg_string << std::endl;
      exit(1);
    }

    // count number of protons & neutrons
    int n_nucleons = std::stoi(pdg_string.substr(6,3)) - 1;
    int n_protons  = std::stoi(pdg_string.substr(3,3));

    // factor proton / neutron ratio into branching ratios
    double proton_frac  = ((double)n_protons) / ((double)n_nucleons);
    double neutron_frac = 1 - proton_frac;

    // set branching ratios, taken from bubble chamber data
    const int n_modes = 16;
    double br [n_modes] = { 0.010, 0.080, 0.100, 0.220,
                            0.360, 0.160, 0.070, 0.020,
                            0.015, 0.065, 0.110, 0.280,
                            0.070, 0.240, 0.100, 0.100 };

    for (int i = 0; i < n_modes; i++) {
      if (i < 7)
        br[i] *= proton_frac;
      else
        br[i] *= neutron_frac;
    }

    // randomly generate a number between 1 and 0
    double p = flatDist.fire();

    // loop through all modes, figure out which one our random number corresponds to
    double threshold = 0;
    for (int i = 0; i < n_modes; i++) {
      threshold += br[i];
      if (p < threshold) {
        // once we've found our mode, return it!
        mode = i + 1;
        return mode;
      }
    }

    // error message, in case the random number selection fails
    std::cout << "Random selection of final state failed!" << std::endl;
    exit(1);
  }

  // if specific annihilation mode specified, just use that
  else {
    int mode = (int) gOptDecayMode;
    return mode;
  }
}

Member Data Documentation

CLHEP::RandFlat evgen::NeutronOsc::flatDist

private

Definition at line 97 of file NeutronOsc_module.cc.

genie::NNBarOscMode_t evgen::NeutronOsc::gOptDecayMode = genie::kNONull

private

Definition at line 96 of file NeutronOsc_module.cc.

const genie::EventRecordVisitorI* evgen::NeutronOsc::mcgen

private

Definition at line 95 of file NeutronOsc_module.cc.

---

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

• larsim/larsim/EventGenerator/GENIE/NeutronOsc_module.cc