diff --git a/src/source.cpp b/src/source.cpp index eace931d6..b79e0fe0f 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -297,17 +297,28 @@ FileSource::FileSource(mcpl_file_t mcpl_file) const mcpl_particle_t *mcpl_particle; //extract particle from mcpl-file mcpl_particle=mcpl_read(mcpl_file); - // check if it is a neutron or a photon. otherwise skip - while ( mcpl_particle->pdgcode!=2112 && mcpl_particle->pdgcode!=22 ) { + // check if it is a neutron, photon, electron, or positron. Otherwise skip. + int pdg=mcpl_particle->pdgcode; + while ( pdg!=2112 && pdg!=22 && pdg!=11 && pdg!=-11) { mcpl_particle=mcpl_read(mcpl_file); - //should check for file exhaustion This could happen if particles are other than - //neutrons or photons + pdg=mcpl_particle->pdgcode; + //should check for file exhaustion. This could happen if particles are other than + //neutrons, photons, electrons, or positrons. } - if(mcpl_particle->pdgcode==2112) { - site_.particle=ParticleType::neutron; - } else if (mcpl_particle->pdgcode==22) { - site_.particle=ParticleType::photon; + switch(pdg){ + case 2112: + omc_particle_.particle=ParticleType::neutron; + break; + case 22: + omc_particle_.particle=ParticleType::photon; + break; + case 11: + omc_particle_.particle=ParticleType::electron; + break; + case -11: + omc_particle_.particle=ParticleType::positron; + break; } //particle is good, convert to openmc-formalism @@ -390,6 +401,100 @@ CustomSourceWrapper::~CustomSourceWrapper() #endif } +<<<<<<< HEAD +<<<<<<< HEAD +#ifdef OPENMC_MCPL +//=========================================================================== +// Read particles from an MCPL-file +//=========================================================================== +MCPLFileSource::MCPLFileSource(std::string path) +{ + // Check if source file exists + if (!file_exists(path)) { + fatal_error(fmt::format("Source file '{}' does not exist.", path)); + } + + // Read the source from a binary file instead of sampling from some + // assumed source distribution + write_message(6, "Reading mcpl source file from {}",path); + + // Open the mcpl file + mcpl_file = mcpl_open_file(path.c_str()); + + //do checks on the mcpl_file to see if particles are many enough. + // should model this on the example source shown in the docs. + n_sites=mcpl_hdr_nparticles(mcpl_file); + + read_source_bank(sites_); +} +<<<<<<< HEAD + +MCPLFileSource::~MCPLFileSource(){ + mcpl_close_file(mcpl_file); +} + +SourceSite MCPLFileSource::sample(uint64_t* seed) const +{ + size_t i_site = sites_.size() * prn(seed); + return sites_[i_site]; +} + +void MCPLFileSource::read_source_bank(vector &sites_) +{ + sites_.resize(n_sites); + for (int i=0;ipdgcode; + while ( pdg!=2112 && pdg!=22 && pdg!=11 && pdg!=-11) { + mcpl_particle=mcpl_read(mcpl_file); + pdg=mcpl_particle->pdgcode; + //should check for file exhaustion. This could happen if particles are other than + //neutrons, photons, electrons, or positrons. + } + + switch(pdg){ + case 2112: + omc_particle_.particle=ParticleType::neutron; + break; + case 22: + omc_particle_.particle=ParticleType::photon; + break; + case 11: + omc_particle_.particle=ParticleType::electron; + break; + case -11: + omc_particle_.particle=ParticleType::positron; + break; + } + + //particle is good, convert to openmc-formalism + omc_particle_.r.x=mcpl_particle->position[0]; + omc_particle_.r.y=mcpl_particle->position[1]; + omc_particle_.r.z=mcpl_particle->position[2]; + + omc_particle_.u.x=mcpl_particle->direction[0]; + omc_particle_.u.y=mcpl_particle->direction[1]; + omc_particle_.u.z=mcpl_particle->direction[2]; + + //mcpl stores kinetic energy in MeV + omc_particle_.E=mcpl_particle->ekin*1e6; + //mcpl stores time in ms + omc_particle_.time=mcpl_particle->time*1e-3; + omc_particle_.wgt=mcpl_particle->weight; + omc_particle_.delayed_group=0; + return omc_particle_; +} +#endif //OPENMC_MCPL //============================================================================== // Non-member functions