diff --git a/config/master_config.xml b/config/master_config.xml index 15460c2891..a4e10ed273 100644 --- a/config/master_config.xml +++ b/config/master_config.xml @@ -275,4 +275,7 @@ NormGenerator.xml NormXSec.xml NormInteractionListGenerator.xml + NucleusGenerator.xml + NucleusGenINCL.xml + NucleusGenTraditional.xml diff --git a/src/Physics/NuclearState/INCLNucleus.cxx b/src/Physics/NuclearState/INCLNucleus.cxx new file mode 100644 index 0000000000..4f668b5db1 --- /dev/null +++ b/src/Physics/NuclearState/INCLNucleus.cxx @@ -0,0 +1,465 @@ +//____________________________________________________________________________ +/*! + + \class genie::INCLNucleus + + \brief INCLXX nuclear model. Implements the NuclearModelI + interface. + + \ref + + \author Liang Liu, (liangliu@fnal.gov) + + \created Oct. 2024 + + \cpright Copyright (c) 2003-2024, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + +*/ +//____________________________________________________________________________ + +#include "Framework/Conventions/GBuild.h" +#ifdef __GENIE_INCL_ENABLED__ + +#include +#include + +#include +#include +#include + +#include "Framework/Messenger/Messenger.h" +#include "Physics/NuclearState/INCLNucleus.h" +#include "Framework/Numerical/Spline.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/GHEP/GHepRecord.h" +#include "Framework/Numerical/RandomGen.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Framework/ParticleData/PDGUtils.h" + + +#include "G4INCLCascade.hh" +#include "G4INCLRandom.hh" +#include "G4INCLStandardPropagationModel.hh" +#include "G4INCLParticleTable.hh" +#include "G4INCLParticle.hh" +#include "G4INCLNuclearMassTable.hh" +#include "G4INCLGlobalInfo.hh" +#include "G4INCLNucleus.hh" + +#include "G4INCLPauliBlocking.hh" + +#include "G4INCLCrossSections.hh" + +#include "G4INCLPhaseSpaceGenerator.hh" + +#include "G4INCLLogger.hh" +#include "G4INCLGlobals.hh" +#include "G4INCLNuclearDensityFactory.hh" + +#include "G4INCLINuclearPotential.hh" + +#include "G4INCLCoulombDistortion.hh" + +#include "G4INCLClustering.hh" + +#include "G4INCLIntersection.hh" + +#include "G4INCLBinaryCollisionAvatar.hh" + +#include "G4INCLCascadeAction.hh" +#include "G4INCLAvatarDumpAction.hh" + +#include +#include +#include + +#include "G4INCLPbarAtrestEntryChannel.hh" + + +#include "G4INCLGeant4Compat.hh" +#include "G4INCLCascade.hh" + +#include "G4INCLClustering.hh" +#include "G4INCLParticle.hh" +#include "G4INCLIAvatar.hh" +#include "G4INCLIPropagationModel.hh" +#include "G4INCLNucleus.hh" +#include "G4INCLStandardPropagationModel.hh" +#include "G4INCLRandom.hh" +#include "G4INCLRanecu.hh" +#include "G4INCLFinalState.hh" +#include "G4INCLParticleTable.hh" +#include "G4INCLKinematicsUtils.hh" +#include "g4inclpauliblocking.hh" +#include "g4inclphasespacegenerator.hh" +#include "g4inclcoulombdistortion.hh" +#include "g4inclbinarycollisionavatar.hh" + +// signal handler (for Linux and GCC) +#include "G4INCLSignalHandling.hh" + +// For I/O +#include "IWriter.hh" +#include "ASCIIWriter.hh" +#include "ProtobufWriter.hh" +#include "INCLTree.hh" +#include "ROOTWriter.hh" +#include "HDF5Writer.hh" + +// For configuration +#include "G4INCLConfig.hh" + +// For logging +#include "G4INCLLogger.hh" + +// Generic de-excitation interface +#include "G4INCLIDeExcitation.hh" + +// ABLA v3p de-excitation +#ifdef INCL_DEEXCITATION_ABLAXX +#include "G4INCLAblaInterface.hh" +#endif + +// ABLACXX de-excitation +#ifdef INCL_DEEXCITATION_ABLACXX +#include "G4INCLAblaxxInterface.hh" +#endif + +// ABLA07 de-excitation +#ifdef INCL_DEEXCITATION_ABLA07 +#include "G4INCLAbla07Interface.hh" +#endif + +// SMM de-excitation +#ifdef INCL_DEEXCITATION_SMM +#include "G4INCLSMMInterface.hh" +#endif + +// GEMINIXX de-excitation +#ifdef INCL_DEEXCITATION_GEMINIXX +#include "G4INCLGEMINIXXInterface.hh" +#endif + + + +// INCL++ +#include "G4INCLConfig.hh" +#include "G4INCLCascade.hh" +#include "G4INCLConfigEnums.hh" +#include "G4INCLParticle.hh" +// signal handler (for Linux and GCC) +#include "G4INCLSignalHandling.hh" + +// Generic de-excitation interface +#include "G4INCLIDeExcitation.hh" + +// ABLA v3p de-excitation +#ifdef INCL_DEEXCITATION_ABLAXX +#include "G4INCLAblaInterface.hh" +#endif + +// ABLA07 de-excitation +#ifdef INCL_DEEXCITATION_ABLA07 +#include "G4INCLAbla07Interface.hh" +#endif + +// SMM de-excitation +#ifdef INCL_DEEXCITATION_SMM +#include "G4INCLSMMInterface.hh" +#endif + +// GEMINIXX de-excitation +#ifdef INCL_DEEXCITATION_GEMINIXX +#include "G4INCLGEMINIXXInterface.hh" +#endif + + + +using std::cout; +using std::endl; + +using namespace genie; + +//____________________________________________________________________________ +INCLNucleus * INCLNucleus::fInstance = 0; + +//____________________________________________________________________________ +INCLNucleus::INCLNucleus():propagationModel_(0) +{ + // this->Load(); + fInstance = 0; + nucleon_index_ = -1; + nucleus_ = nullptr; + theConfig_ = nullptr; +} +//____________________________________________________________________________ +INCLNucleus::~INCLNucleus() +{ + // if(!gAbortingInErr) { + // cout << "INCLNucleus singleton dtor: Deleting inputs... " << endl; + // } + // delete fNuclSupprD2; +} +//____________________________________________________________________________ +INCLNucleus * INCLNucleus::Instance() +{ + if(fInstance == 0) { + LOG("NuclData", pINFO) << "INCLNucleus late initialization"; + // static INCLNucleus::Cleaner cleaner; + // cleaner.DummyMethodAndSilentCompiler(); + fInstance = new INCLNucleus; + fInstance->theConfig_ = new G4INCL::Config(); + fInstance->nucleus_ = nullptr; + fInstance->init(); + } + return fInstance; +} + + +void INCLNucleus::init(){ + LOG("NuclData", pINFO) << "init()"; + theConfig_->init(); + theConfig_->setINCLXXDataFilePath("/root/inclxx/inclxx-v6.33.1-e5857a1/data"); // FIXME:: using config to set path + // initialize INCL model + G4INCL::Random::initialize(theConfig_); + // Select the Pauli and CDPP blocking algorithms + G4INCL::Pauli::initialize(theConfig_); + // Set the phase-space generator + G4INCL::PhaseSpaceGenerator::initialize(theConfig_); + // Select the Coulomb-distortion algorithm: + G4INCL::CoulombDistortion::initialize(theConfig_); + // Select the clustering algorithm: + G4INCL::Clustering::initialize(theConfig_); + // Initialize the INCL particle table: + G4INCL::ParticleTable::initialize(theConfig_); + // Initialize the value of cutNN in BinaryCollisionAvatar + G4INCL::BinaryCollisionAvatar::setCutNN(theConfig_->getCutNN()); + // Initialize the value of strange cross section bias + G4INCL::BinaryCollisionAvatar::setBias(theConfig_->getBias()); + // Set the cross-section set + G4INCL::CrossSections::initialize(theConfig_); + //theConfig_->setLocalEnergyBBType(G4INCL::NeverLocalEnergy); + //theConfig_->setLocalEnergyPiType(G4INCL::NeverLocalEnergy); + + + // Propagation model is responsible for finding avatars and + // transporting the particles. In principle this step is "hidden" + // behind an abstract interface and the rest of the system does not + // care how the transportation and avatar finding is done. This + // should allow us to "easily" experiment with different avatar + // finding schemes and even to support things like curved + // trajectories in the future. + propagationModel_ = new G4INCL::StandardPropagationModel(theConfig_->getLocalEnergyBBType(),theConfig_->getLocalEnergyPiType(),theConfig_->getHadronizationTime()); + if(theConfig_->getCascadeActionType() == G4INCL::AvatarDumpActionType) + cascadeAction_ = new G4INCL::AvatarDumpAction(); + else + cascadeAction_ = new G4INCL::CascadeAction(); +// cascadeAction_->beforeRunAction(theConfig_); + + + +} + + +void INCLNucleus::initialize(const GHepRecord * evrec){ + + LOG("NuclData", pINFO) << "initialize()"; + + // initialize according process Event in INCL + + GHepParticle * nucleus = evrec->TargetNucleus(); + G4INCL::ParticleSpecies targetSpecies = G4INCL::ParticleSpecies(nucleus->Name()); + theConfig_->setTargetA(targetSpecies.theA); + theConfig_->setTargetZ(targetSpecies.theZ); + theConfig_->setTargetS(targetSpecies.theS); + LOG("NuclData", pINFO) << "initialize()"; + // define Nucleus and initialize it + + if(nucleus_){ + if(!nucleus_->getStore()->getParticles().empty()){ + if(nucleus_->getA() == evrec->TargetNucleus()->A() && nucleus_->getZ() == evrec->TargetNucleus()->Z()){ +// LOG("NuclData", pINFO) << nucleus_->print(); + return; + } + else{ + nucleus_->deleteParticles(); + nucleus_->getStore()->clear(); + delete nucleus_; + } + } + } + + // ReInitialize the bias vector + G4INCL::Particle::INCLBiasVector.clear(); + //Particle::INCLBiasVector.Clear(); + G4INCL::Particle::nextBiasedCollisionID = 0; + + // Set the target and the projectile + // implement prepare reaction + + // Reset the forced-transparent flag + // forceTransparent = false; FIXME + // + // Initialise the maximum universe radius + // INCL initialize universe radius according to particle species, + // kenetic energy, and nucleus type. void INCL::initUniverseRadius(ParticleSpecies const &p, const double kineticEnergy, const int A, const int Z) + // FIXME: I'm not sure if we need the interaction distance for neutrino scattering + double rMax = 0.0; + const double pMaximumRadius = G4INCL::ParticleTable::getMaximumNuclearRadius(G4INCL::Proton, targetSpecies.theA, targetSpecies.theZ); + const double nMaximumRadius = G4INCL::ParticleTable::getMaximumNuclearRadius(G4INCL::Neutron, targetSpecies.theA, targetSpecies.theZ); + const double maximumRadius = std::max(pMaximumRadius, nMaximumRadius); + rMax = std::max(maximumRadius, rMax); + maxUniverseRadius_ = rMax; + + // FIXME the last two parameters need to be configed + // theConfig_, G4INCL::ParticleTable::getMaximumNuclearRadius(G4INCL::Proton, targetSpecies.theA, targetSpecies.theZ) + // G4INCL::NType + nucleus_ = new G4INCL::Nucleus(targetSpecies.theA, targetSpecies.theZ, targetSpecies.theS, + theConfig_, maxUniverseRadius_, G4INCL::Def); + nucleus_->getStore()->getBook().reset(); + nucleus_->initializeParticles(); + propagationModel_->setNucleus(nucleus_); + + // initialize max interaction distance + // FIXME: in INCL, composite has non-zero max interaction distance. + maxInteractionDistance_ = 0; + + // set the min Remnant size to be 4 + // the min remnant is alpha particle + // FIXME: it is only works for nuclei with large A + minRemnantSize_ = 4; + + // cascade action is not related to simulation + // it is just output the casade to file FIXME + //cascadeAction_->beforeCascadeAction(propagationModel_); + // + // INCL need to decide whether the cascade can be ran or not + // For genie, we need to run casecade for every events + // const bool canRunCascade = preCascade(projectileSpecies, kineticEnergy); + + LOG("INCLNucleus", pDEBUG) << nucleus_->getStore()->getParticles().at(2)->getPotentialEnergy() ; + LOG("INCLNucleus", pDEBUG) << nucleus_->getStore()->getParticles().at(2)->getEnergy() - nucleus_->getStore()->getParticles().at(2)->getPotentialEnergy() ; + LOG("INCLNucleus", pDEBUG) << nucleus_->getStore()->getParticles().at(2)->getMomentum().print() ; + LOG("INCLNucleus", pDEBUG) << nucleus_->getStore()->getParticles().at(2)->getPosition().print() ; + + GHepParticle * nucleon = evrec->HitNucleon(); + LOG("INCLNucleus", pDEBUG) << "hit nucleon pdg : " << nucleon->Pdg(); + + RandomGen * rnd = RandomGen::Instance(); + if(pdg::IsProton(nucleon->Pdg())) + nucleon_index_ = rnd->RndGen().Integer(targetSpecies.theZ); + else if(pdg::IsNeutron(nucleon->Pdg())) + nucleon_index_ = rnd->RndGen().Integer(targetSpecies.theA - targetSpecies.theZ) + targetSpecies.theZ; + else + exit(1); + hitNucleon_ = nucleus_->getStore()->getParticles().at(nucleon_index_); +} + +void INCLNucleus::reset(const GHepRecord * evrec){ + // nucleus must exsit! + if(!nucleus_) LOG("INCLNucleus", pFATAL) << "nucleus doesn't exsit!"; + // can't reset a nucleus with different type + if(!(nucleus_->getA() == evrec->TargetNucleus()->A() && nucleus_->getZ() == evrec->TargetNucleus()->Z())) + LOG("INCLNucleus", pFATAL) << "you are try to reset a nucleus with different type!"; + // reset the nucleus + if(nucleus_){ + nucleus_->deleteParticles(); + nucleus_->getStore()->clear(); + nucleus_->initializeParticles(); + nucleus_->getStore()->getBook().reset(); + + GHepParticle * nucleon = evrec->HitNucleon(); + LOG("INCLNucleus", pDEBUG) << "hit nucleon pdg : " << nucleon->Pdg(); + RandomGen * rnd = RandomGen::Instance(); + if(pdg::IsProton(nucleon->Pdg())) + nucleon_index_ = rnd->RndGen().Integer(evrec->TargetNucleus()->Z()); + else if(pdg::IsNeutron(nucleon->Pdg())) + nucleon_index_ = rnd->RndGen().Integer(evrec->TargetNucleus()->A() - evrec->TargetNucleus()->Z()) + evrec->TargetNucleus()->Z(); + else + exit(1); + hitNucleon_ = nucleus_->getStore()->getParticles().at(nucleon_index_); + + + delete propagationModel_; + propagationModel_ = new G4INCL::StandardPropagationModel(theConfig_->getLocalEnergyBBType(),theConfig_->getLocalEnergyPiType(),theConfig_->getHadronizationTime()); + + } +} +TVector3 INCLNucleus::getHitNucleonPosition(){ + if(nucleus_ && nucleon_index_ != -1){ + TVector3 v3_(999999.,999999.,999999.); + v3_.SetXYZ(nucleus_->getStore()->getParticles().at(nucleon_index_)->getPosition().getX(), + nucleus_->getStore()->getParticles().at(nucleon_index_)->getPosition().getY(), + nucleus_->getStore()->getParticles().at(nucleon_index_)->getPosition().getZ()); + return v3_; + } + else + exit(1); +} +TVector3 INCLNucleus::getHitNucleonMomentum(){ + if(nucleus_ && nucleon_index_ != -1){ + TVector3 v3_(999999.,999999.,999999.); + v3_.SetXYZ(nucleus_->getStore()->getParticles().at(nucleon_index_)->getMomentum().getX(), + nucleus_->getStore()->getParticles().at(nucleon_index_)->getMomentum().getY(), + nucleus_->getStore()->getParticles().at(nucleon_index_)->getMomentum().getZ()); + return v3_; + } + else + exit(1); + +} +double INCLNucleus::getHitNucleonEnergy(){ + if(nucleus_ && nucleon_index_ != -1){ + return nucleus_->getStore()->getParticles().at(nucleon_index_)->getEnergy(); + } + else + exit(1); +} + +double INCLNucleus::getHitNucleonMass(){ + if(nucleus_ && nucleon_index_ != -1){ + return nucleus_->getStore()->getParticles().at(nucleon_index_)->getMass(); + } + else + exit(1); +} + +double INCLNucleus::getMass(){ + if(nucleus_ && nucleon_index_ != -1){ + return nucleus_->getMass(); + } + else + exit(1); +} + +G4INCL::Nucleus * INCLNucleus::getNuclues(){ + if(nucleus_ && nucleon_index_ != -1) + return nucleus_; + else + exit(1); +} + +G4INCL::Particle * INCLNucleus::getHitParticle(){ + return hitNucleon_; +} + +G4INCL::StandardPropagationModel * INCLNucleus::getPropagationModel(){ + return propagationModel_; +} + +double INCLNucleus::getRemovalEnergy(){ + if(nucleus_ && nucleon_index_ != -1){ + double removal_energy = 0; + double nucleon_mass = nucleus_->getStore()->getParticles().at(nucleon_index_)->getRealMass(); + double mag = nucleus_->getStore()->getParticles().at(nucleon_index_)->getMomentum().mag(); + removal_energy = TMath::Sqrt(mag*mag + nucleon_mass*nucleon_mass) - nucleus_->getStore()->getParticles().at(nucleon_index_)->getEnergy(); + return removal_energy; + } + else + exit(1); +} + +#endif // __GENIE_INCL_ENABLED__ + diff --git a/src/Physics/NuclearState/INCLNucleus.h b/src/Physics/NuclearState/INCLNucleus.h new file mode 100644 index 0000000000..251cdcecd5 --- /dev/null +++ b/src/Physics/NuclearState/INCLNucleus.h @@ -0,0 +1,98 @@ +//____________________________________________________________________________ +/*! + + \class genie::INCLNucleus + + \brief INCLXX nuclear model. Implements the NuclearModelI + interface. + + \ref + + \author Liang Liu, (liangliu@fnal.gov) + + \created Oct. 2024 + + \cpright Copyright (c) 2003-2024, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + +*/ +//____________________________________________________________________________ + + +#include "Framework/Conventions/GBuild.h" +#ifdef __GENIE_INCL_ENABLED__ + +#ifndef _INCL_NUCLEUS_H_ +#define _INCL_NUCLEUS_H_ + +// ROOT +#include "TVector3.h" + + +// INCL++ +// For configuration +#include "G4INCLConfig.hh" +#include "G4INCLNucleus.hh" +// GENIE +#include "Framework/GHEP/GHepRecord.h" +#include "G4INCLParticle.hh" +#include "G4INCLNucleus.hh" +#include "G4INCLIPropagationModel.hh" +#include "G4INCLStandardPropagationModel.hh" +#include "G4INCLCascadeAction.hh" +#include "G4INCLEventInfo.hh" +#include "G4INCLGlobalInfo.hh" +#include "G4INCLLogger.hh" +#include "G4INCLConfig.hh" +#include "G4INCLRootFinder.hh" + + + +namespace genie { + + class INCLNucleus { + + public: + static INCLNucleus * Instance (void); + void initialize(const GHepRecord * evrec); + void reset(const GHepRecord * evrec); + + TVector3 getHitNucleonPosition(); + TVector3 getHitNucleonMomentum(); + double getHitNucleonEnergy(); + double getHitNucleonMass(); + double getMass(); + double getRemovalEnergy(); + G4INCL::Nucleus * getNuclues(); + G4INCL::StandardPropagationModel * getPropagationModel(); + G4INCL::Particle *getHitParticle(); + G4INCL::Config *getConfig(){return theConfig_;} + + private: + INCLNucleus(); + ~INCLNucleus(); + void init(); + + static INCLNucleus *fInstance; + + TVector3 v3_; // position of initial nucleon + TVector3 p3_; // fermi momentum of initial nucleon + double energy_; // off-shell energy of initial nucleon + G4INCL::Config *theConfig_; + G4INCL::Nucleus *nucleus_; + G4INCL::StandardPropagationModel *propagationModel_; + G4INCL::CascadeAction *cascadeAction_; + G4INCL::Particle *hitNucleon_; + + int nucleon_index_; + + double maxUniverseRadius_; + double maxInteractionDistance_; + double minRemnantSize_; + + + }; + +} // genie namespace +#endif // _INCL_NUCLEUS_H_ +#endif // __GENIE_INCL_ENABLED__ diff --git a/src/Physics/NuclearState/LinkDef.h b/src/Physics/NuclearState/LinkDef.h index c9207f1ca8..269a14c93b 100644 --- a/src/Physics/NuclearState/LinkDef.h +++ b/src/Physics/NuclearState/LinkDef.h @@ -21,6 +21,9 @@ #pragma link C++ class genie::FermiMomentumTablePool; #pragma link C++ class genie::EffectiveSF; #pragma link C++ class genie::FermiMover; +#pragma link C++ class genie::NucleusGenerator; +#pragma link C++ class genie::NucleusGenINCL; +#pragma link C++ class genie::NucleusGenTraditional; #pragma link C++ class genie::PauliBlocker; #pragma link C++ class genie::SRCNuclearRecoil; #pragma link C++ class genie::SecondNucleonEmissionI; diff --git a/src/Physics/NuclearState/Makefile b/src/Physics/NuclearState/Makefile index f1e69f1034..4464854b2a 100644 --- a/src/Physics/NuclearState/Makefile +++ b/src/Physics/NuclearState/Makefile @@ -25,4 +25,19 @@ install : install-inc install-lib # include $(GENIE)/src/make/Make.std-package-targets +# Only need these here ... not in every package + +ifeq ($(strip $(GOPT_ENABLE_INCL)),YES) + # extra flags, include paths, and libraries to link to + CXXFLAGS += $(INCL_DFLAGS) + CPP_INCLUDES += $(INCL_INCLUDES) + ROOT_DICT_GEN_INCLUDES += $(INCL_INCLUDES) $(INCL_DFLAGS) + EXTRA_EXT_LIBS += $(INCL_LIBRARIES) +else + $(info $(PACKAGE) not built against INCL++) +endif + + + + FORCE: diff --git a/src/Physics/NuclearState/NucleusGenINCL.cxx b/src/Physics/NuclearState/NucleusGenINCL.cxx new file mode 100644 index 0000000000..fa67ab7cf1 --- /dev/null +++ b/src/Physics/NuclearState/NucleusGenINCL.cxx @@ -0,0 +1,336 @@ +//____________________________________________________________________________ +/*! + +\class genie::NucleusGenINCL + +\brief It visits the event record & computes a Fermi motion momentum for + initial state nucleons bound in nuclei. + Is a concrete implementation of the EventRecordVisitorI interface. + +\author Liang Liu + Fermi National Accelerator Laboratory + +\created October 17, 2024 + +\cpright Copyright (c) 2003-2024, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + +*/ +//____________________________________________________________________________ + +#include "Framework/Conventions/GBuild.h" +#ifdef __GENIE_INCL_ENABLED__ +#include + +#include +#include +#include +#include + +#include "Framework/Algorithm/AlgFactory.h" +#include "Framework/Algorithm/AlgConfigPool.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/Conventions/Units.h" +#include "Physics/NuclearState/NucleusGenINCL.h" + +#include "Physics/NuclearState/NuclearModel.h" +#include "Physics/NuclearState/NuclearModelI.h" +#include "Framework/EventGen/EVGThreadException.h" +#include "Framework/GHEP/GHepRecord.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/GHEP/GHepStatus.h" +#include "Framework/GHEP/GHepFlags.h" +#include "Framework/Interaction/Interaction.h" +#include "Framework/Messenger/Messenger.h" +#include "Physics/NuclearState/FermiMomentumTablePool.h" +#include "Physics/NuclearState/FermiMomentumTable.h" +#include "Framework/Numerical/RandomGen.h" +#include "Framework/ParticleData/PDGLibrary.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Framework/Utils/KineUtils.h" +#include "Physics/NuclearState/NuclearUtils.h" +#include "Physics/NuclearState/INCLNucleus.h" + +using namespace genie; +using namespace genie::constants; + +//___________________________________________________________________________ +NucleusGenINCL::NucleusGenINCL() : +EventRecordVisitorI("genie::NucleusGenINCL") +{ + +} +//___________________________________________________________________________ +NucleusGenINCL::NucleusGenINCL(string config) : +EventRecordVisitorI("genie::NucleusGenINCL", config) +{ + +} +//___________________________________________________________________________ +NucleusGenINCL::~NucleusGenINCL() +{ + +} + +//___________________________________________________________________________ +void NucleusGenINCL::ProcessEventRecord(GHepRecord * evrec) const +{ + // skip if not a nuclear target + if(! evrec->Summary()->InitState().Tgt().IsNucleus()) return; + + // skip if no hit nucleon is set + if(! evrec->HitNucleon()) return; + + LOG("NucleusGenINCL", pINFO) << "Adding final state nucleus"; + INCLNucleus *incl_nucleus = INCLNucleus::Instance(); + LOG("NucleusGenINCL", pINFO) << "Adding final state nucleus"; + incl_nucleus->initialize(evrec); + incl_nucleus->reset(evrec); + incl_nucleus->initialize(evrec); + LOG("NucleusGenINCL", pINFO) << "Adding final state nucleus"; + + // give hit nucleon a vertex + this->setInitialStateVertex(evrec); + LOG("NucleusGenINCL", pINFO) << "Adding final state nucleus"; + // give hit nucleon a Fermi momentum + this->setInitialStateMomentum(evrec); + LOG("NucleusGenINCL", pINFO) << "Adding final state nucleus"; + + // handle the addition of the recoil nucleon + // TODO: INCL has it own SRC model +// if ( fSecondEmitter ) fSecondEmitter -> ProcessEventRecord( evrec ) ; + + // add a recoiled nucleus remnant + this->setTargetNucleusRemnant(evrec); + LOG("NucleusGenINCL", pINFO) << "Adding final state nucleus"; +} + +//___________________________________________________________________________ +// using INCL model to get the position and momentum of +// Hit nucleon +void NucleusGenINCL::setInitialStateVertex(GHepRecord * evrec) const{ + + INCLNucleus *incl_nucleus = INCLNucleus::Instance(); + +// generate a vtx and set it to all GHEP physical particles + Interaction * interaction = evrec->Summary(); + GHepParticle * nucltgt = evrec->TargetNucleus(); + TVector3 vtx(9999999.,999999.,999999.); + if(!nucltgt){ + vtx.SetXYZ(0.,0.,0.); + }else{ + double A = nucltgt->A(); + // vtx = GenerateVertex(interaction,A); + + const ProcessInfo & proc_info = interaction->ProcInfo(); + bool is_coh = proc_info.IsCoherentProduction() || proc_info.IsCoherentElastic(); + bool is_ve = proc_info.IsInverseMuDecay() || + proc_info.IsIMDAnnihilation() || + proc_info.IsNuElectronElastic() || + proc_info.IsGlashowResonance() || + proc_info.IsPhotonResonance() || + proc_info.IsPhotonCoherent(); + + + if(is_coh||is_ve) { + // ** For COH or ve- set a vertex positon on the nuclear boundary + // + // LOG("Vtx", pINFO) << "Setting vertex on the nuclear boundary"; + // double phi = 2*kPi * rnd->RndFsi().Rndm(); + // double cosphi = TMath::Cos(phi); + // double sinphi = TMath::Sin(phi); + // double costheta = -1 + 2 * rnd->RndFsi().Rndm(); + // double sintheta = TMath::Sqrt(1-costheta*costheta); + // vtx.SetX(R*sintheta*cosphi); + // vtx.SetY(R*sintheta*sinphi); + // vtx.SetZ(R*costheta); + } + else { + vtx = incl_nucleus->getHitNucleonPosition(); + } + + } + + LOG("NucleusGenINCL", pINFO) << "Position"; + vtx.Print(); + // Copy the vertex info to the particles already in the event record + // + TObjArrayIter piter(evrec); + GHepParticle * p = 0; + while( (p = (GHepParticle *) piter.Next()) ) + { + if(pdg::IsPseudoParticle(p->Pdg())) continue; + if(pdg::IsIon (p->Pdg())) continue; + + LOG("NucleusGenINCL", pINFO) << "Setting vertex position for: " << p->Name(); + p->SetPosition(vtx.x(), vtx.y(), vtx.z(), 0.); + } + + + +} + +void NucleusGenINCL::setInitialStateMomentum(GHepRecord * evrec) const{ + Interaction * interaction = evrec -> Summary(); + InitialState * init_state = interaction -> InitStatePtr(); + Target * tgt = init_state -> TgtPtr(); + + // do nothing for non-nuclear targets + if(!tgt->IsNucleus()) return; + + TLorentzVector * p4 = tgt->HitNucP4Ptr(); + + // do nothing if the struct nucleon 4-momentum was set (eg as part of the + // initial state selection) + if(p4->Px()>0 || p4->Py()>0 || p4->Pz()>0) return; + + // access the hit nucleon and target nucleus at the GHEP record + GHepParticle * nucleon = evrec->HitNucleon(); + GHepParticle * nucleus = evrec->TargetNucleus(); + assert(nucleon); + assert(nucleus); + + // initialize INCL nucleus model + // INCL nucleus model sample all nucleons with r-p correlation + INCLNucleus *incl_nucleus = INCLNucleus::Instance(); + G4INCL::Nucleus *incl_nuc = incl_nucleus->getNuclues(); + TLorentzVector p4tgt; + p4tgt.SetPx(incl_nuc->getMomentum().getX() / 1000.); + p4tgt.SetPy(incl_nuc->getMomentum().getY() / 1000. ); + p4tgt.SetPz(incl_nuc->getMomentum().getZ() / 1000. ); + p4tgt.SetE(incl_nuc->getEnergy() / 1000.); + init_state->SetTgtP4(p4tgt); + nucleus->SetMomentum(p4tgt); + + + // get a random nucleon with respect to the isospin of evrec->HitNucleon(); + // the removal energy maybe not necessary + TVector3 p3 = incl_nucleus->getHitNucleonMomentum(); + double hit_nucleon_energy = incl_nucleus->getHitNucleonEnergy(); + double w = incl_nucleus->getRemovalEnergy(); + //-- update the struck nucleon 4p at the interaction summary and at + // the GHEP record + p4->SetPx(p3.Px()/1000.); + p4->SetPy(p3.Py()/1000.); + p4->SetPz(p3.Pz()/1000.); + p4->SetE (hit_nucleon_energy/1000.); + LOG("NucleusGenINCL", pINFO) << "Momentum"; + p4->Print(); + + nucleon->SetMomentum(*p4); // update GHEP value + nucleon->SetRemovalEnergy(w); // FIXME this may be not necessary + + + // Sometimes, for interactions near threshold, Fermi momentum might bring + // the neutrino energy in the nucleon rest frame below threshold (for the + // selected interaction). In this case mark the event as unphysical and + // abort the current thread. + const KPhaseSpace & kps = interaction->PhaseSpace(); + if(!kps.IsAboveThreshold()) { + LOG("NucleusGenINCL", pNOTICE) + << "Event below threshold after generating Fermi momentum"; + + double Ethr = kps.Threshold(); + double Ev = init_state->ProbeE(kRfHitNucRest); + LOG("NucleusGenINCL", pNOTICE) + << "Ev (@ nucleon rest frame) = " << Ev << ", Ethr = " << Ethr; + + evrec->EventFlags()->SetBitNumber(kBelowThrNRF, true); + genie::exceptions::EVGThreadException exception; + exception.SetReason("E < Ethr after generating nucleon Fermi momentum"); + exception.SwitchOnFastForward(); + throw exception; + } + +} + +void NucleusGenINCL::setTargetNucleusRemnant(GHepRecord * evrec)const{ +// add the remnant nuclear target at the GHEP record + + LOG("NucleusGenINCL", pINFO) << "Adding final state nucleus"; + + double Px = 0; + double Py = 0; + double Pz = 0; + double E = 0; + + GHepParticle * nucleus = evrec->TargetNucleus(); + int A = nucleus->A(); + int Z = nucleus->Z(); + + int fd = nucleus->FirstDaughter(); + int ld = nucleus->LastDaughter(); + + INCLNucleus *incl_nucleus = INCLNucleus::Instance(); + + G4INCL::Nucleus *incl_nuc = incl_nucleus->getNuclues(); + + + for(int id = fd; id <= ld; id++) { + + // compute A,Z for final state nucleus & get its PDG code and its mass + GHepParticle * particle = evrec->Particle(id); + assert(particle); + int pdgc = particle->Pdg(); + bool is_p = pdg::IsProton (pdgc); + bool is_n = pdg::IsNeutron(pdgc); + + if (is_p) Z--; + if (is_p || is_n) A--; + + Px += particle->Px(); + Py += particle->Py(); + Pz += particle->Pz(); + E += particle->E(); + + }//daughters + + TParticlePDG * remn = 0; + int ipdgc = pdg::IonPdgCode(A, Z); + remn = PDGLibrary::Instance()->Find(ipdgc); + if(!remn) { + LOG("HadronicVtx", pFATAL) + << "No particle with [A = " << A << ", Z = " << Z + << ", pdgc = " << ipdgc << "] in PDGLibrary!"; + assert(remn); + } + + double Mi = nucleus->Mass(); + Px = incl_nuc->getMomentum().getX()/1000. - Px; + Py = incl_nuc->getMomentum().getY()/1000. - Py; + Pz = incl_nuc->getMomentum().getZ()/1000. - Pz; + E = incl_nuc->getEnergy()/1000. - E; + + // Add the nucleus to the event record + LOG("FermiMover", pINFO) + << "Adding nucleus [A = " << A << ", Z = " << Z + << ", pdgc = " << ipdgc << "]"; + + int imom = evrec->TargetNucleusPosition(); + evrec->AddParticle( + ipdgc,kIStStableFinalState, imom,-1,-1,-1, Px,Py,Pz,E, 0,0,0,0); +} + + +//___________________________________________________________________________ +void NucleusGenINCL::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void NucleusGenINCL::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void NucleusGenINCL::LoadConfig(void) +{ + +} +//____________________________________________________________________________ + + +#endif // end __GENIE_INCL_ENABLED__ diff --git a/src/Physics/NuclearState/NucleusGenINCL.h b/src/Physics/NuclearState/NucleusGenINCL.h new file mode 100644 index 0000000000..a09c093bbc --- /dev/null +++ b/src/Physics/NuclearState/NucleusGenINCL.h @@ -0,0 +1,77 @@ +//____________________________________________________________________________ +/*! + +\class genie::NucleusGenINCL + +\brief It visits the event record & computes a Fermi motion momentum for + initial state nucleons bound in nuclei. + Is a concrete implementation of the EventRecordVisitorI interface. + +\author Liang Liu + Fermi National Accelerator Laboratory + +\created October 17, 2024 + +\cpright Copyright (c) 2003-2024, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + +*/ +//____________________________________________________________________________ + +#include "Framework/Conventions/GBuild.h" +#ifdef __GENIE_INCL_ENABLED__ + +#ifndef _NUCLEUS_GEN_INCL_H_ +#define _NUCLEUS_GEN_INCL_H_ + +#include "Framework/EventGen/EventRecordVisitorI.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Physics/NuclearState/FermiMomentumTable.h" +#include "Framework/Interaction/Target.h" +#include "Physics/NuclearState/SRCNuclearRecoil.h" +#include "Physics/NuclearState/SecondNucleonEmissionI.h" + +namespace genie { + + +class NucleusGenINCL : public EventRecordVisitorI { + +public : + NucleusGenINCL(); + NucleusGenINCL(string config); + ~NucleusGenINCL(); + + //-- implement the EventRecordVisitorI interface + void ProcessEventRecord(GHepRecord * event_rec) const; + + //-- overload the Algorithm::Configure() methods to load private data + // members from configuration options + void Configure(const Registry & config); + void Configure(string config); + +private: + + void setInitialStateVertex (GHepRecord * evrec) const; ///< give hit nucleon a position + + void setInitialStateMomentum (GHepRecord * evrec) const; ///< give hit nucleon a momentum + + void setTargetNucleusRemnant (GHepRecord * evrec) const; ///< add a recoiled nucleus remnant + + // TODO; check the energy-momentum conservation + + void LoadConfig (void); + + // bool fKeepNuclOnMassShell; ///< keep hit bound nucleon on the mass shell? + // bool fMomDepErmv; ///< use momentum dependent calculation of Ermv + // const NuclearModelI * fNuclModel; ///< nuclear model + + // const SecondNucleonEmissionI * fSecondEmitter ; + +}; + +} // genie namespace + + + +#endif // _NUCLEUS_GEN_INCL_H_ +#endif // end __GENIE_INCL_ENABLED__ diff --git a/src/Physics/NuclearState/NucleusGenTraditional.cxx b/src/Physics/NuclearState/NucleusGenTraditional.cxx new file mode 100644 index 0000000000..c7f412f2b9 --- /dev/null +++ b/src/Physics/NuclearState/NucleusGenTraditional.cxx @@ -0,0 +1,111 @@ +//____________________________________________________________________________ +/*! + +\class genie::NucleusGenTraditional + +\brief It visits the event record & computes a Fermi motion momentum for + initial state nucleons bound in nuclei. + Is a concrete implementation of the EventRecordVisitorI interface. + +\author Liang Liu + Fermi National Accelerator Laboratory + +\created October 17, 2024 + +\cpright Copyright (c) 2003-2024, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + +*/ +//____________________________________________________________________________ + +#include + +#include +#include +#include +#include + +#include "Framework/Algorithm/AlgFactory.h" +#include "Framework/Algorithm/AlgConfigPool.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/Conventions/Units.h" +#include "Physics/NuclearState/NucleusGenTraditional.h" + +#include "Physics/NuclearState/NuclearModel.h" +#include "Physics/NuclearState/NuclearModelI.h" +#include "Framework/EventGen/EVGThreadException.h" +#include "Framework/GHEP/GHepRecord.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/GHEP/GHepStatus.h" +#include "Framework/GHEP/GHepFlags.h" +#include "Framework/Interaction/Interaction.h" +#include "Framework/Messenger/Messenger.h" +#include "Physics/NuclearState/FermiMomentumTablePool.h" +#include "Physics/NuclearState/FermiMomentumTable.h" +#include "Framework/Numerical/RandomGen.h" +#include "Framework/ParticleData/PDGLibrary.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Framework/Utils/KineUtils.h" +#include "Physics/NuclearState/NuclearUtils.h" +#include "Physics/NuclearState/INCLNucleus.h" + +using namespace genie; +using namespace genie::constants; + +//___________________________________________________________________________ +NucleusGenTraditional::NucleusGenTraditional() : +EventRecordVisitorI("genie::NucleusGenTraditional") +{ + +} +//___________________________________________________________________________ +NucleusGenTraditional::NucleusGenTraditional(string config) : +EventRecordVisitorI("genie::NucleusGenTraditional", config) +{ + +} +//___________________________________________________________________________ +NucleusGenTraditional::~NucleusGenTraditional() +{ + +} + +//___________________________________________________________________________ +void NucleusGenTraditional::ProcessEventRecord(GHepRecord * evrec) const +{ + // skip if not a nuclear target + if(! evrec->Summary()->InitState().Tgt().IsNucleus()) return; + + LOG("NucleusGenTraditional", pINFO) << "Adding final state nucleus"; + + fVertexGenerator->ProcessEventRecord(evrec); + fFermiMover->ProcessEventRecord(evrec); + +} + +//___________________________________________________________________________ +void NucleusGenTraditional::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void NucleusGenTraditional::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void NucleusGenTraditional::LoadConfig(void) +{ + + + fFermiMover = nullptr; + fVertexGenerator = nullptr; + fFermiMover = dynamic_cast (this->SubAlg("FermiMover")); + fVertexGenerator = dynamic_cast (this->SubAlg("VertexGenerator")); + +} +//____________________________________________________________________________ + diff --git a/src/Physics/NuclearState/NucleusGenTraditional.h b/src/Physics/NuclearState/NucleusGenTraditional.h new file mode 100644 index 0000000000..257c3eed9f --- /dev/null +++ b/src/Physics/NuclearState/NucleusGenTraditional.h @@ -0,0 +1,57 @@ +//____________________________________________________________________________ +/*! + +\class genie::NucleusGenTraditional + +\brief It visits the event record & computes a Fermi motion momentum for + initial state nucleons bound in nuclei. + Is a concrete implementation of the EventRecordVisitorI interface. + +\author Liang Liu + Fermi National Accelerator Laboratory + +\created October 17, 2024 + +\cpright Copyright (c) 2003-2024, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + +*/ +//____________________________________________________________________________ + + +#ifndef _NUCLEUS_GEN_TRADITIONAL_H_ +#define _NUCLEUS_GEN_TRADITIONAL_H_ + +#include "Framework/EventGen/EventRecordVisitorI.h" +#include "Framework/GHEP/GHepParticle.h" + +namespace genie { + + +class NucleusGenTraditional : public EventRecordVisitorI { + +public : + NucleusGenTraditional(); + NucleusGenTraditional(string config); + ~NucleusGenTraditional(); + + //-- implement the EventRecordVisitorI interface + void ProcessEventRecord(GHepRecord * event_rec) const; + + //-- overload the Algorithm::Configure() methods to load private data + // members from configuration options + void Configure(const Registry & config); + void Configure(string config); + +private: + + void LoadConfig (void); + + const EventRecordVisitorI *fFermiMover; + const EventRecordVisitorI *fVertexGenerator; + +}; + +} // genie namespace + +#endif // _NUCLEUS_GEN_TRADITIONAL_H_ diff --git a/src/Physics/NuclearState/NucleusGenerator.cxx b/src/Physics/NuclearState/NucleusGenerator.cxx new file mode 100644 index 0000000000..9c37174178 --- /dev/null +++ b/src/Physics/NuclearState/NucleusGenerator.cxx @@ -0,0 +1,87 @@ +//____________________________________________________________________________ +/*! + +\class genie::NucleusGenerator + +\brief It visits the event record & computes a Fermi motion momentum for + initial state nucleons bound in nuclei. + Is a concrete implementation of the EventRecordVisitorI interface. + +\author Liang Liu + Fermi National Accelerator Laboratory + +\created October 17, 2024 + +\cpright Copyright (c) 2003-2024, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + +*/ +//____________________________________________________________________________ + +#include + +#include +#include +#include +#include + +#include "Framework/Algorithm/AlgFactory.h" +#include "Framework/Algorithm/AlgConfigPool.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/Conventions/Units.h" +#include "Physics/NuclearState/NucleusGenerator.h" + +using namespace genie; +using namespace genie::constants; + +//___________________________________________________________________________ +NucleusGenerator::NucleusGenerator() : +EventRecordVisitorI("genie::NucleusGenerator") +{ + +} +//___________________________________________________________________________ +NucleusGenerator::NucleusGenerator(string config) : +EventRecordVisitorI("genie::NucleusGenerator", config) +{ + +} +//___________________________________________________________________________ +NucleusGenerator::~NucleusGenerator() +{ + +} + +//___________________________________________________________________________ +void NucleusGenerator::ProcessEventRecord(GHepRecord * evrec) const +{ + // skip if not a nuclear target + //if(! evrec->Summary()->InitState().Tgt().IsNucleus()) return; + LOG("NucleusGenerator", pINFO) << "Adding final state nucleus"; + fNucleusGen->ProcessEventRecord(evrec); +} + +//___________________________________________________________________________ +void NucleusGenerator::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void NucleusGenerator::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void NucleusGenerator::LoadConfig(void) +{ + RgKey nuclkey = "NucleusGenerator"; + fNucleusGen = nullptr; + //fNuclModel = dynamic_cast (this->SubAlg(nuclkey)); + // + LOG("NucleusGenerator", pINFO) << "Hello world!"; + fNucleusGen = dynamic_cast (this->SubAlg("Nuclear-Model")); +} +//____________________________________________________________________________ + diff --git a/src/Physics/NuclearState/NucleusGenerator.h b/src/Physics/NuclearState/NucleusGenerator.h new file mode 100644 index 0000000000..0c994fdc56 --- /dev/null +++ b/src/Physics/NuclearState/NucleusGenerator.h @@ -0,0 +1,56 @@ +//____________________________________________________________________________ +/*! + + \class genie::NucleusGenerator + + \brief It visits the event record & computes a Fermi motion momentum for + initial state nucleons bound in nuclei. + Is a concrete implementation of the EventRecordVisitorI interface. + + \author Liang Liu + Fermi National Accelerator Laboratory + + \created October 17, 2024 + + \cpright Copyright (c) 2003-2024, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + +*/ +//____________________________________________________________________________ + +#ifndef _NUCLEUS_GENERATOR_H_ +#define _NUCLEUS_GENERATOR_H_ + +#include "Framework/EventGen/EventRecordVisitorI.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Physics/NuclearState/FermiMomentumTable.h" +#include "Framework/Interaction/Target.h" +#include "Physics/NuclearState/SRCNuclearRecoil.h" +#include "Physics/NuclearState/SecondNucleonEmissionI.h" + +namespace genie { + + + class NucleusGenerator : public EventRecordVisitorI { + + public : + NucleusGenerator(); + NucleusGenerator(string config); + ~NucleusGenerator(); + + //-- implement the EventRecordVisitorI interface + void ProcessEventRecord(GHepRecord * event_rec) const; + + //-- overload the Algorithm::Configure() methods to load private data + // members from configuration options + void Configure(const Registry & config); + void Configure(string config); + + private: + + void LoadConfig (void); + + const EventRecordVisitorI *fNucleusGen; + }; +} // genie namespace +#endif // _NUCLEUS_GENERATOR_H_