diff --git a/config/AR23_20i/ModelConfiguration.xml b/config/AR23_20i/ModelConfiguration.xml index 69e2930a1..63e847358 100644 --- a/config/AR23_20i/ModelConfiguration.xml +++ b/config/AR23_20i/ModelConfiguration.xml @@ -44,6 +44,19 @@ STFC, Rutherford Appleton Laboratory genie::SpectralFunc1d/Default --> + + genie::BohrElectronVelocity/Default + - genie::IMDXSec/Default + genie::XSecOnElectron/Default + genie::ElectronVelocityMap/Default + genie::PXSecOnElectron/Default + genie::PXSecOnElectron/Default diff --git a/config/BohrElectronVelocity.xml b/config/BohrElectronVelocity.xml new file mode 100644 index 000000000..f88208339 --- /dev/null +++ b/config/BohrElectronVelocity.xml @@ -0,0 +1,17 @@ + + + + + + + + + diff --git a/config/EventGenerator.xml b/config/EventGenerator.xml index d1a418a16..a05f75856 100644 --- a/config/EventGenerator.xml +++ b/config/EventGenerator.xml @@ -425,34 +425,37 @@ XSecModel alg Yes Cross section model used at the thread - 5 + 6 genie::InitialStateAppender/Default - genie::VertexGenerator/Default - genie::NuEKinematicsGenerator/Default - genie::NuEPrimaryLeptonGenerator/Default - genie::NuETargetRemnantGenerator/Default + genie::VertexGenerator/Default + genie::ElectronVelocityMap/Default + genie::NuEKinematicsGenerator/Default + genie::NuEPrimaryLeptonGenerator/Default + genie::NuETargetRemnantGenerator/Default genie::NuEInteractionListGenerator/IMD - 5 + 6 genie::InitialStateAppender/Default genie::VertexGenerator/Default - genie::NuEKinematicsGenerator/Default - genie::NuEPrimaryLeptonGenerator/Default - genie::NuETargetRemnantGenerator/Default + genie::ElectronVelocityMap/Default + genie::NuEKinematicsGenerator/Default + genie::NuEPrimaryLeptonGenerator/Default + genie::NuETargetRemnantGenerator/Default genie::NuEInteractionListGenerator/IMD-ANH - 5 + 6 genie::InitialStateAppender/Default genie::VertexGenerator/Default - genie::NuEKinematicsGenerator/Default - genie::NuEPrimaryLeptonGenerator/Default - genie::NuETargetRemnantGenerator/Default + genie::ElectronVelocityMap/Default + genie::NuEKinematicsGenerator/Default + genie::NuEPrimaryLeptonGenerator/Default + genie::NuETargetRemnantGenerator/Default genie::NuEInteractionListGenerator/NUE-EL diff --git a/config/G00_00a/EventGenerator.xml b/config/G00_00a/EventGenerator.xml index 2d873f1ad..911aa0130 100644 --- a/config/G00_00a/EventGenerator.xml +++ b/config/G00_00a/EventGenerator.xml @@ -322,14 +322,15 @@ XSecModel alg Yes Cross section model used at the thread genie::NuEInteractionListGenerator/IMD-ANH - + - 5 + 6 genie::InitialStateAppender/Default genie::VertexGenerator/Default - genie::NuEKinematicsGenerator/Default - genie::NuEPrimaryLeptonGenerator/Default - genie::NuETargetRemnantGenerator/Default + genie::ElectronVelocityMap/Default + genie::NuEKinematicsGenerator/Default + genie::NuEPrimaryLeptonGenerator/Default + genie::NuETargetRemnantGenerator/Default genie::NuEInteractionListGenerator/NUE-EL diff --git a/config/G00_00a/ModelConfiguration.xml b/config/G00_00a/ModelConfiguration.xml index 106bca2db..b227139fc 100644 --- a/config/G00_00a/ModelConfiguration.xml +++ b/config/G00_00a/ModelConfiguration.xml @@ -43,6 +43,17 @@ STFC, Rutherford Appleton Laboratory genie::SpectralFunc1d/Default --> + + genie::StaticElectronVelocity/Default - + + + genie::StaticElectronVelocity/Default + + genie::StaticElectronVelocity/Default + - - + + + genie::BohrElectronVelocity/Default + + genie::BohrElectronVelocity/Default + + genie::BohrElectronVelocity/Default + + + + genie::BohrElectronVelocity/Default + - + + + genie::BohrElectronVelocity/Default + + genie::BohrElectronVelocity/Default + + genie::BohrElectronVelocity/Default + + genie::BohrElectronVelocity/Default + + + genie::BohrElectronVelocity/Default + + genie::BohrElectronVelocity/Default + + genie::BohrElectronVelocity/Default + + genie::BohrElectronVelocity/Default + + genie::BohrElectronVelocity/Default + + genie::BohrElectronVelocity/Default + + genie::BohrElectronVelocity/Default + + genie::BohrElectronVelocity/Default - + + + genie::BohrElectronVelocity/Default - + + + genie::BohrElectronVelocity/Default + + + genie::BohrElectronVelocity/Default + genie::BohrElectronVelocity/Default + + - genie::IMDXSec/Default + genie::XSecOnElectron/Default + genie::ElectronVelocityMap/Default + genie::PXSecOnElectron/Default + genie::PXSecOnElectron/Default diff --git a/config/IMDXSec.xml b/config/IMDXSec.xml deleted file mode 100644 index 38eee9e52..000000000 --- a/config/IMDXSec.xml +++ /dev/null @@ -1,21 +0,0 @@ - - - - - - - - adaptive - 40000 - 0.001 - - - - diff --git a/config/Messenger.xml b/config/Messenger.xml index 52b6ff705..a1d0be951 100644 --- a/config/Messenger.xml +++ b/config/Messenger.xml @@ -109,7 +109,6 @@ WARN WARN NOTICE - NOTICE WARN WARN NOTICE diff --git a/config/Messenger_whisper.xml b/config/Messenger_whisper.xml index a8450bb2c..efd740454 100644 --- a/config/Messenger_whisper.xml +++ b/config/Messenger_whisper.xml @@ -104,7 +104,6 @@ FATAL FATAL FATAL - FATAL FATAL FATAL FATAL diff --git a/config/NuElectronPXSec.xml b/config/NuElectronPXSec.xml index 8fe2b2992..a4f191b06 100644 --- a/config/NuElectronPXSec.xml +++ b/config/NuElectronPXSec.xml @@ -7,15 +7,21 @@ Configuration for the NuElectronPXSec xsec algorithm. Configurable Parameters: .................................................................................... -Name Type Optional Comment Default -WeinbergAngle double No CommonParam[WeakInt] +Name Type Optional Comment Default +WeinbergAngle double No CommonParam[WeakInt] +XSec-Integrator alg No genie::XSecOnElectron/Default +Electron-Velocity alg No genie::ElectronVelocityMap/Default +N-Integration-Samples int No genie::PXSecOnElectron/Default +NuE-XSecRelError double No genie::PXSecOnElectron/Default .................................................................................... --> WeakInt - - genie::NuElectronXSec/Default + genie::XSecOnElectron/Default + genie::ElectronVelocityMap/Default + genie::PXSecOnElectron/Default + genie::PXSecOnElectron/Default diff --git a/config/PXSecOnElectron.xml b/config/PXSecOnElectron.xml new file mode 100644 index 000000000..d533c4aaf --- /dev/null +++ b/config/PXSecOnElectron.xml @@ -0,0 +1,27 @@ + + + + + + + + + 10000 + 0.0001 + genie::ElectronVelocityMap/Default + genie::NuElectronXSec/Default + + + + + diff --git a/config/StaticElectronVelocity.xml b/config/StaticElectronVelocity.xml new file mode 100644 index 000000000..a48332360 --- /dev/null +++ b/config/StaticElectronVelocity.xml @@ -0,0 +1,17 @@ + + + + + + + + + diff --git a/config/NuElectronXSec.xml b/config/XSecOnElectron.xml similarity index 92% rename from config/NuElectronXSec.xml rename to config/XSecOnElectron.xml index 9c30e669a..0104db344 100644 --- a/config/NuElectronXSec.xml +++ b/config/XSecOnElectron.xml @@ -3,7 +3,7 @@ EventGenerator.xml FermiMover.xml + Default.xml + BohrElectronVelocity.xml + StaticElectronVelocity.xml HadronTransporter.xml HAIntranuke.xml HAIntranuke2018.xml @@ -160,7 +163,6 @@ CEvNSXSec.xml DFRXSec.xml AlamSimoAtharVacasSKXSec.xml - IMDXSec.xml RESXSec.xml MECXSec.xml NuElectronXSec.xml @@ -198,6 +200,7 @@ KLVOxygenIBDPXSec.xml PattonCEvNSPXSec.xml NuElectronPXSec.xml + PXSecOnElectron.xml DMElectronPXSec.xml GLRESPXSec.xml HENuElPXSec.xml diff --git a/src/Framework/Interaction/InitialState.cxx b/src/Framework/Interaction/InitialState.cxx index b74b5ef75..ea96385ba 100644 --- a/src/Framework/Interaction/InitialState.cxx +++ b/src/Framework/Interaction/InitialState.cxx @@ -11,6 +11,9 @@ CMEnergy() method added by Andy Furmanski (Univ. of Manchester) and Joe Johnston (Univ of Pittsburgh) + + Modified GetProbeP4 and GetTargetP4 to boost to electron rest frame + Brinden Carlson (University of Florida) */ //____________________________________________________________________________ @@ -316,6 +319,20 @@ TLorentzVector * InitialState::GetTgtP4(RefFrame_t ref_frame) const return p4; break; } + //------------------ STRUCK ELECTRON REST FRAME: + case (kRfHitElRest) : + { + // make sure that 'struck electron' properties were set in + // the electron target object + if (!fTgt->HitEleIsSet()){ + return nullptr; + } + double inv_mass = fTgt->HitEleP4Ptr()->Mag(); + TLorentzVector * p4; + p4->SetVectM(TVector3(),inv_mass); + return p4; + break; + } default: LOG("Interaction", pERROR) << "Uknown reference frame"; } @@ -374,6 +391,23 @@ TLorentzVector * InitialState::GetProbeP4(RefFrame_t ref_frame) const break; } + //----------------- STRUCK ELECTRON REST FRAME + case (kRfHitElRest) : + { + //Ensure target is electron + if (fTgt->HitEleP4Ptr() == 0){ + return nullptr; + } + TLorentzVector * pele4 = fTgt->HitEleP4Ptr(); + + // compute velocity vector + + auto boost = pele4->BoostVector(); + TLorentzVector * p4 = new TLorentzVector(*fProbeP4); + p4->Boost(-boost); + return p4; + break; + } default: LOG("Interaction", pERROR) << "Uknown reference frame"; diff --git a/src/Framework/Interaction/Target.cxx b/src/Framework/Interaction/Target.cxx index e999998eb..bcd862218 100644 --- a/src/Framework/Interaction/Target.cxx +++ b/src/Framework/Interaction/Target.cxx @@ -5,6 +5,9 @@ Costas Andreopoulos University of Liverpool & STFC Rutherford Appleton Laboratory + + Changes required to implement the Electron Velocity module + were installed by Brinden Carlson (Univ. of Florida) */ //____________________________________________________________________________ @@ -80,7 +83,8 @@ fA(0), fTgtPDG(0), fHitNucPDG(0), fHitSeaQrk(false), -fHitNucP4(0) +fHitNucP4(0), +fHitEleP4(0) { } @@ -105,12 +109,14 @@ void Target::Init(void) fHitQrkPDG = 0; fHitSeaQrk = false; fHitNucP4 = new TLorentzVector(0,0,0,kNucleonMass); + fHitEleP4 = new TLorentzVector(0,0,0,kElectronMass); fHitNucRad = 0.; } //___________________________________________________________________________ void Target::CleanUp(void) { delete fHitNucP4; + delete fHitEleP4; } //___________________________________________________________________________ void Target::Copy(const Target & tgt) @@ -144,6 +150,10 @@ void Target::Copy(const Target & tgt) // a nucleon (p or n) or a di-nucleon cluster (p+p, p+n, n+n) this->ForceHitNucValidity(); } + if (tgt.fHitNucPDG == 0 && tgt.fHitQrkPDG == 0 && tgt.fHitSeaQrk == 0){ //No interaction with nucleus -> interaction with electron + const TLorentzVector& p4 = *(tgt.fHitEleP4); + *fHitEleP4 = *tgt.fHitEleP4 ; + } } //___________________________________________________________________________ void Target::SetId(int pdgc) @@ -192,6 +202,12 @@ void Target::SetHitNucP4(const TLorentzVector & p4) fHitNucP4 = new TLorentzVector(p4); } //___________________________________________________________________________ +void Target::SetHitEleP4(const TLorentzVector & p4) +{ + if(fHitEleP4) delete fHitEleP4; + fHitEleP4 = new TLorentzVector(p4); +} +//___________________________________________________________________________ void Target::SetHitSeaQrk(bool tf) { fHitSeaQrk = tf; @@ -254,6 +270,15 @@ TLorentzVector * Target::HitNucP4Ptr(void) const return fHitNucP4; } //___________________________________________________________________________ +TLorentzVector * Target::HitEleP4Ptr(void) const +{ + if(!fHitEleP4) { + LOG("Target", pWARN) << "Returning NULL struck electron 4-momentum"; + return 0; + } + return fHitEleP4; +} +//___________________________________________________________________________ bool Target::IsFreeNucleon(void) const { return (fA == 1 && (fZ == 0 || fZ == 1)); @@ -264,6 +289,11 @@ bool Target::IsProton(void) const return (fA == 1 && fZ == 1); } //___________________________________________________________________________ +bool Target::IsElectron(void) const +{ + return (fA == 0 && fZ == 0); //No nucleons +} +//___________________________________________________________________________ bool Target::IsNeutron(void) const { return (fA == 1 && fZ == 0); @@ -289,6 +319,15 @@ bool Target::HitNucIsSet(void) const return ok; } //___________________________________________________________________________ +bool Target::HitEleIsSet(void) const +{ + bool ok = fHitNucPDG == 0 && + fHitQrkPDG == 0 && + fHitSeaQrk == 0; //Hit no quarks + + return ok; +} +//___________________________________________________________________________ bool Target::HitQrkIsSet(void) const { return ( @@ -417,6 +456,10 @@ void Target::Print(ostream & stream) const << utils::print::BoolAsYNString(this->HitSeaQrk()) << ")"; } + if( this->HitEleIsSet() ) { + TParticlePDG * q = PDGLibrary::Instance()->Find(fHitQrkPDG); + stream << " struck electron = "; + } } //___________________________________________________________________________ bool Target::Compare(const Target & target) const diff --git a/src/Framework/Interaction/Target.h b/src/Framework/Interaction/Target.h index 1552ca145..3c53490e9 100644 --- a/src/Framework/Interaction/Target.h +++ b/src/Framework/Interaction/Target.h @@ -11,6 +11,9 @@ \author Costas Andreopoulos University of Liverpool & STFC Rutherford Appleton Laboratory + Changes required to implement the Electron Velocity module + were installed by Brinden Carlson (Univ. of Florida) + \created May 03, 2004 \cpright Copyright (c) 2003-2023, The GENIE Collaboration @@ -58,6 +61,7 @@ using TObject::Copy; void SetId (int Z, int A); void SetHitNucPdg (int pdgc); void SetHitNucP4 (const TLorentzVector & p4); + void SetHitEleP4 (const TLorentzVector & p4); void SetHitNucPosition (double r); void SetHitQrkPdg (int pdgc); void SetHitSeaQrk (bool tf); @@ -73,11 +77,13 @@ using TObject::Copy; double Charge (void) const; bool IsFreeNucleon (void) const; bool IsProton (void) const; + bool IsElectron (void) const; bool IsNeutron (void) const; bool IsNucleus (void) const; bool IsParticle (void) const; bool IsValidNucleus (void) const; bool HitNucIsSet (void) const; + bool HitEleIsSet (void) const; bool HitQrkIsSet (void) const; bool HitSeaQrk (void) const; bool IsEvenEven (void) const; @@ -90,6 +96,8 @@ using TObject::Copy; const TLorentzVector & HitNucP4 (void) const { return *this->HitNucP4Ptr(); } TLorentzVector * HitNucP4Ptr (void) const; + const TLorentzVector & HitEleP4 (void) const { return *this->HitEleP4Ptr(); } + TLorentzVector * HitEleP4Ptr (void) const; //-- Copy, reset, compare, print itself and build string code void Reset (void); @@ -121,9 +129,10 @@ using TObject::Copy; int fHitQrkPDG; ///< hit quark PDG code bool fHitSeaQrk; ///< hit quark from sea? TLorentzVector * fHitNucP4; ///< hit nucleon 4p + TLorentzVector * fHitEleP4; ///< hit electron 4p - changes index double fHitNucRad; ///< hit nucleon position -ClassDef(Target,2) +ClassDef(Target,3) }; } // genie namespace diff --git a/src/Physics/Common/ElectronVelocityMap.cxx b/src/Physics/Common/ElectronVelocityMap.cxx new file mode 100644 index 000000000..c3fc0a317 --- /dev/null +++ b/src/Physics/Common/ElectronVelocityMap.cxx @@ -0,0 +1,107 @@ +///____________________________________________________________________________ +/* + Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + \brief It provides the correct ElectronVelociy model to be used + for a given atom. + + \author Marco Roda + University of Liverpool + + \created March 28, 2023 + + For the class documentation see the corresponding header file. +*/ +//____________________________________________________________________________ + +#include "Physics/Common/ElectronVelocityMap.h" + +#include "Framework/Algorithm/AlgConfigPool.h" +#include "Framework/ParticleData/PDGUtils.h" + +using namespace genie; + +ElectronVelocityMap::ElectronVelocityMap() : + ElectronVelocity("genie::ElectronVelocityMap", "Default") { ; } +//___________________________________________________________________________ +ElectronVelocityMap::ElectronVelocityMap(string config) : +ElectronVelocity("genie::ElectronVelocityMap", config) { ; } +//___________________________________________________________________________ +void ElectronVelocityMap::Configure(string config) +{ + Algorithm::Configure(config); + + Registry * algos = AlgConfigPool::Instance() -> GlobalParameterList() ; + Registry r( "ElectronVelocitylMap", false ) ; + + // copy in local pool relevant configurations + RgIMap entries = algos -> GetItemMap(); + const std::string keyStart = "ElectronVelocity"; + for( RgIMap::const_iterator it = entries.begin(); it != entries.end(); ++it ) { + + if( it -> first.compare(0, keyStart.size(), keyStart.c_str()) == 0 ) { + r.Set( it -> first, algos -> GetAlg(it->first ) ) ; + } + + } + + ElectronVelocity::Configure(r) ; + +} +//___________________________________________________________________________ +void ElectronVelocityMap::LoadConfig() { + fDefGlobalVelocity = nullptr; + fSpecificModels.clear(); + + // load default global model (should work for all nuclei) + RgAlg dgmodel ; + GetParam( "ElectronVelocity", dgmodel ) ; + + LOG("ElcetronVelocity", pINFO) + << "Default global electron velocity model: " << dgmodel; + fDefGlobalVelocity = dynamic_cast ( this -> SubAlg( "ElectronVelocity" ) ) ; + assert(fDefGlobalVelocity); + + // We're looking for keys that match this string + const std::string keyStart = "ElectronVelocity@Pdg="; + // Looking in both of these registries + RgIMap entries = GetConfig().GetItemMap(); + + for(RgIMap::const_iterator it = entries.begin(); it != entries.end(); ++it){ + const std::string& key = it->first; + // Does it start with the right string? + if(key.compare(0, keyStart.size(), keyStart.c_str()) == 0){ + // The rest is the PDG code + const int pdg = atoi(key.c_str()+keyStart.size()); + const int Z = pdg::IonPdgCodeToZ(pdg); + //const int A = pdg::IonPdgCodeToA(pdg); + + RgAlg rgmodel = GetConfig().GetAlg(key) ; + LOG("ElectronVelocity", pNOTICE) + << "Atom =" << pdg + << " -> refined velocity model: " << rgmodel; + const ElectronVelocity * model = + dynamic_cast ( + this -> SubAlg(key) ) ; + assert(model); + fSpecificModels.insert(map::value_type(Z,model)); + } + } +} +//___________________________________________________________________________ +void ElectronVelocityMap::InitializeVelocity(Interaction & interaction) const { + + const auto & model = SelectModel(interaction.InitState().Tgt()); + + model.InitializeVelocity(interaction); + +} +//___________________________________________________________________________ +const ElectronVelocity & ElectronVelocityMap::SelectModel(const Target & t) const { + auto it = fSpecificModels.find(t.Z()); + + if ( it != fSpecificModels.end()) return *(it->second) ; + else return * fDefGlobalVelocity; +} + diff --git a/src/Physics/Common/ElectronVelocityMap.h b/src/Physics/Common/ElectronVelocityMap.h new file mode 100644 index 000000000..85064e666 --- /dev/null +++ b/src/Physics/Common/ElectronVelocityMap.h @@ -0,0 +1,59 @@ +//____________________________________________________________________________ +/*! + +\class genie::ElectronVelocityMap + +\brief This class is a hook for Electron velocity distributions and allows associating each + one of them with specific nuclei. + Is a concrete implementation of the ElectronVelocity interface. + + \author Marco Roda + University of Liverpool + + \created March 28, 2023 + +\cpright Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + +*/ +//____________________________________________________________________________ + +#ifndef _ELECTRON_VELOCITY_MAP_H_ +#define _ELECTRON_VELOCITY_MAP_H_ + +#include "Physics/NuclearState/ElectronVelocity.h" + +namespace genie { + +class ElectronVelocityMap : public ElectronVelocity { + +public : + + ElectronVelocityMap(); + ElectronVelocityMap(string config); + + ElectronVelocityMap(const ElectronVelocityMap & ) = delete; + ElectronVelocityMap(ElectronVelocityMap && ) = delete; + + + virtual ~ElectronVelocityMap() {;} + + //-- overload the ElectronVelocity::Configure() methods + // to load data from ModelConfig.xml + void Configure(string config) override final ; + + void InitializeVelocity(Interaction & interaction) const override; //Give initial velocity + +protected: + void LoadConfig () override; + const ElectronVelocity & SelectModel(const Target &) const; + +private: + const ElectronVelocity * fDefGlobalVelocity = nullptr; + map fSpecificModels ; + // the key is the Z of the atom + +}; + +} // genie namespace +#endif // _ELECTRON_VELOCITY_MAP_H_ diff --git a/src/Physics/Common/LinkDef.h b/src/Physics/Common/LinkDef.h index 8aba58b88..dc84f2e09 100644 --- a/src/Physics/Common/LinkDef.h +++ b/src/Physics/Common/LinkDef.h @@ -24,4 +24,5 @@ #pragma link C++ class genie::XSecLinearCombinations; #pragma link C++ class genie::QvalueShifter; +#pragma link C++ class genie::ElectronVelocityMap; #endif diff --git a/src/Physics/NuElectron/EventGen/NuEKinematicsGenerator.cxx b/src/Physics/NuElectron/EventGen/NuEKinematicsGenerator.cxx index cc9b4eb56..d10f32cc6 100644 --- a/src/Physics/NuElectron/EventGen/NuEKinematicsGenerator.cxx +++ b/src/Physics/NuElectron/EventGen/NuEKinematicsGenerator.cxx @@ -5,6 +5,10 @@ Costas Andreopoulos University of Liverpool & STFC Rutherford Appleton Laboratory + + Changes required to calculate energy in electron rest frame + were installed by Brinden Carlson (Univ. of Florida) + */ //____________________________________________________________________________ @@ -210,10 +214,10 @@ double NuEKinematicsGenerator::ComputeMaxXSec( double NuEKinematicsGenerator::Energy(const Interaction * interaction) const { // Override the base class Energy() method to cache the max xsec for the -// neutrino energy in the LAB rather than in the hit nucleon rest frame. +// neutrino energy in the electron rest frame. const InitialState & init_state = interaction->InitState(); - double E = init_state.ProbeE(kRfLab); + double E = init_state.ProbeE(kRfHitElRest); // return E; } //___________________________________________________________________________ diff --git a/src/Physics/NuElectron/EventGen/NuEPrimaryLeptonGenerator.cxx b/src/Physics/NuElectron/EventGen/NuEPrimaryLeptonGenerator.cxx index d7a0336f8..c05b7f533 100644 --- a/src/Physics/NuElectron/EventGen/NuEPrimaryLeptonGenerator.cxx +++ b/src/Physics/NuElectron/EventGen/NuEPrimaryLeptonGenerator.cxx @@ -5,6 +5,9 @@ Costas Andreopoulos University of Liverpool & STFC Rutherford Appleton Laboratory + + Changes required to boost from electron rest frame to lab frame + were installed by Brinden Carlson (Univ. of Florida) */ //____________________________________________________________________________ @@ -47,7 +50,14 @@ void NuEPrimaryLeptonGenerator::ProcessEventRecord(GHepRecord * evrec) const // This method generates the final state primary lepton for NuE events Interaction * interaction = evrec->Summary(); - const InitialState & init_state = interaction->InitState(); + + // Boost vector for [LAB] <-> [Electron Rest Frame] transforms + TVector3 beta = this->EleRestFrame2Lab(*evrec); // Get boost of hit + + // Neutrino 4p + TLorentzVector p4v(*evrec->Probe()->GetP4()); // v 4p @ LAB + p4v.Boost(-1.*beta); + // Get selected kinematics double y = interaction->Kine().y(true); @@ -58,7 +68,7 @@ void NuEPrimaryLeptonGenerator::ProcessEventRecord(GHepRecord * evrec) const assert(pdgc!=0); // Compute the neutrino and muon energy - double Ev = init_state.ProbeE(kRfLab); + double Ev = p4v.E(); double El = (1-y)*Ev; LOG("LeptonicVertex", pINFO) @@ -109,8 +119,9 @@ void NuEPrimaryLeptonGenerator::ProcessEventRecord(GHepRecord * evrec) const TVector3 p3l(pltx,plty,plp); p3l.RotateUz(unit_nudir); - // Lepton 4-momentum in the LAB + // Lepton 4-momentum in the Ele rest frame TLorentzVector p4l(p3l,El); + p4l.Boost(beta); //Boost back to lab // Create a GHepParticle and add it to the event record this->AddToEventRecord(evrec, pdgc, p4l); @@ -119,3 +130,17 @@ void NuEPrimaryLeptonGenerator::ProcessEventRecord(GHepRecord * evrec) const this->SetPolarization(evrec); } //___________________________________________________________________________ +TVector3 NuEPrimaryLeptonGenerator::EleRestFrame2Lab(GHepRecord & evrec) const +{ +// Velocity for an active Lorentz transform taking the final state primary +// lepton from the [electron rest frame] --> [LAB] + + Interaction * interaction = evrec.Summary(); + const InitialState & init_state = interaction->InitState(); + + const TLorentzVector & pele4 = init_state.Tgt().HitEleP4(); //[@LAB] + TVector3 beta = pele4.BoostVector(); + + return beta; +} +//___________________________________________________________________________ diff --git a/src/Physics/NuElectron/EventGen/NuEPrimaryLeptonGenerator.h b/src/Physics/NuElectron/EventGen/NuEPrimaryLeptonGenerator.h index bd0818d03..c60b05d1e 100644 --- a/src/Physics/NuElectron/EventGen/NuEPrimaryLeptonGenerator.h +++ b/src/Physics/NuElectron/EventGen/NuEPrimaryLeptonGenerator.h @@ -32,6 +32,8 @@ public : //-- implement the EventRecordVisitorI interface void ProcessEventRecord(GHepRecord * event_rec) const; + + TVector3 EleRestFrame2Lab (GHepRecord & ev) const; }; } // genie namespace diff --git a/src/Physics/NuElectron/XSection/BardinIMDRadCorPXSec.cxx b/src/Physics/NuElectron/XSection/BardinIMDRadCorPXSec.cxx index cad76e740..b2721b224 100644 --- a/src/Physics/NuElectron/XSection/BardinIMDRadCorPXSec.cxx +++ b/src/Physics/NuElectron/XSection/BardinIMDRadCorPXSec.cxx @@ -11,7 +11,6 @@ #include #include -#include "Physics/XSectionIntegration/XSecIntegratorI.h" #include "Framework/Conventions/GBuild.h" #include "Framework/Conventions/Constants.h" #include "Framework/Conventions/RefFrame.h" @@ -19,19 +18,20 @@ #include "Framework/Messenger/Messenger.h" #include "Framework/Utils/KineUtils.h" #include "Framework/Numerical/GSLUtils.h" +#include "Physics/NuElectron/XSection/PXSecOnElectron.h" using namespace genie; using namespace genie::constants; //____________________________________________________________________________ BardinIMDRadCorPXSec::BardinIMDRadCorPXSec() : -XSecAlgorithmI("genie::BardinIMDRadCorPXSec") +PXSecOnElectron::PXSecOnElectron("genie::BardinIMDRadCorPXSec","Default") { } //____________________________________________________________________________ BardinIMDRadCorPXSec::BardinIMDRadCorPXSec(string config) : -XSecAlgorithmI("genie::BardinIMDRadCorPXSec", config) +PXSecOnElectron::PXSecOnElectron("genie::BardinIMDRadCorPXSec", config) { } @@ -96,18 +96,6 @@ double BardinIMDRadCorPXSec::XSec( return xsec; } //____________________________________________________________________________ -double BardinIMDRadCorPXSec::Integral(const Interaction * interaction) const -{ - double xsec = fXSecIntegrator->Integrate(this,interaction); - return xsec; -} -//____________________________________________________________________________ -bool BardinIMDRadCorPXSec::ValidProcess(const Interaction * interaction) const -{ - if(interaction->TestBit(kISkipProcessChk)) return true; - return true; -} -//____________________________________________________________________________ double BardinIMDRadCorPXSec::Fa(double re, double r, double y) const { double y2 = y * y; @@ -253,13 +241,13 @@ double BardinIMDRadCorPXSec::C(int i, int k, double r) const //____________________________________________________________________________ void BardinIMDRadCorPXSec::Configure(const Registry & config) { - Algorithm::Configure(config); + PXSecOnElectron::Configure(config); this->LoadConfig(); } //____________________________________________________________________________ void BardinIMDRadCorPXSec::Configure(string param_set) { - Algorithm::Configure(param_set); + PXSecOnElectron::Configure(param_set); this->LoadConfig(); } //____________________________________________________________________________ @@ -268,10 +256,6 @@ void BardinIMDRadCorPXSec::LoadConfig(void) ////fIntegrator = //// dynamic_cast (this->SubAlg("Integrator")); ///// assert(fIntegrator); - - fXSecIntegrator = - dynamic_cast (this->SubAlg("XSec-Integrator")); - assert(fXSecIntegrator); } //____________________________________________________________________________ // Auxiliary scalar function for internal integration diff --git a/src/Physics/NuElectron/XSection/BardinIMDRadCorPXSec.h b/src/Physics/NuElectron/XSection/BardinIMDRadCorPXSec.h index 4f9998811..5612c9d20 100644 --- a/src/Physics/NuElectron/XSection/BardinIMDRadCorPXSec.h +++ b/src/Physics/NuElectron/XSection/BardinIMDRadCorPXSec.h @@ -19,6 +19,8 @@ \author Costas Andreopoulos University of Liverpool & STFC Rutherford Appleton Laboratory + B. Carlson modified to interface PXSecOnElectron + \created February 14, 2005 \cpright Copyright (c) 2003-2023, The GENIE Collaboration @@ -30,23 +32,19 @@ #define _BARDIN_IMD_RADIATIVE_CORRECTIONS_PARTIAL_XSEC_H_ #include -#include "Framework/EventGen/XSecAlgorithmI.h" +#include "Physics/NuElectron/XSection/PXSecOnElectron.h" namespace genie { -class XSecIntegratorI; - -class BardinIMDRadCorPXSec : public XSecAlgorithmI { +class BardinIMDRadCorPXSec : public PXSecOnElectron { public: BardinIMDRadCorPXSec(); BardinIMDRadCorPXSec(string config); virtual ~BardinIMDRadCorPXSec(); - // XSecAlgorithmI interface implementation + // PXSecOnElectron interface implementation double XSec (const Interaction * i, KinePhaseSpace_t k) const; - double Integral (const Interaction * i) const; - bool ValidProcess (const Interaction * i) const; // Override the Algorithm::Configure methods to load configuration // data to private data members @@ -67,7 +65,6 @@ class BardinIMDRadCorPXSec : public XSecAlgorithmI { // Private data members // const IntegratorI * fIntegrator; ///< num integrator for BardinIMDRadCorIntegrand - const XSecIntegratorI * fXSecIntegrator; ///< differential x-sec integrator }; } // genie namespace diff --git a/src/Physics/NuElectron/XSection/IMDAnnihilationPXSec.cxx b/src/Physics/NuElectron/XSection/IMDAnnihilationPXSec.cxx index 88845a6a2..2057743c6 100644 --- a/src/Physics/NuElectron/XSection/IMDAnnihilationPXSec.cxx +++ b/src/Physics/NuElectron/XSection/IMDAnnihilationPXSec.cxx @@ -10,13 +10,13 @@ #include "Framework/Algorithm/AlgConfigPool.h" #include "Framework/Conventions/GBuild.h" #include "Framework/Conventions/Controls.h" -#include "Physics/XSectionIntegration/XSecIntegratorI.h" #include "Framework/Conventions/Constants.h" #include "Framework/Conventions/RefFrame.h" #include "Physics/NuElectron/XSection/IMDAnnihilationPXSec.h" #include "Framework/Messenger/Messenger.h" #include "Framework/ParticleData/PDGUtils.h" #include "Framework/Utils/KineUtils.h" +#include "Physics/NuElectron/XSection/PXSecOnElectron.h" using namespace genie; using namespace genie::constants; @@ -24,13 +24,13 @@ using namespace genie::controls; //____________________________________________________________________________ IMDAnnihilationPXSec::IMDAnnihilationPXSec() : -XSecAlgorithmI("genie::IMDAnnihilationPXSec") +PXSecOnElectron::PXSecOnElectron("genie::IMDAnnihilationPXSec","Default") { } //____________________________________________________________________________ IMDAnnihilationPXSec::IMDAnnihilationPXSec(string config) : -XSecAlgorithmI("genie::IMDAnnihilationPXSec", config) +PXSecOnElectron::PXSecOnElectron("genie::IMDAnnihilationPXSec", config) { } @@ -87,12 +87,6 @@ double IMDAnnihilationPXSec::XSec( return xsec; } //____________________________________________________________________________ -double IMDAnnihilationPXSec::Integral(const Interaction * interaction) const -{ - double xsec = fXSecIntegrator->Integrate(this,interaction); - return xsec; -} -//____________________________________________________________________________ bool IMDAnnihilationPXSec::ValidProcess(const Interaction * interaction) const { if(interaction->TestBit(kISkipProcessChk)) return true; @@ -107,21 +101,18 @@ bool IMDAnnihilationPXSec::ValidKinematics(const Interaction* interaction) const //____________________________________________________________________________ void IMDAnnihilationPXSec::Configure(const Registry & config) { - Algorithm::Configure(config); + PXSecOnElectron::Configure(config); this->LoadConfig(); } //____________________________________________________________________________ void IMDAnnihilationPXSec::Configure(string config) { - Algorithm::Configure(config); + PXSecOnElectron::Configure(config); this->LoadConfig(); } //____________________________________________________________________________ void IMDAnnihilationPXSec::LoadConfig(void) { - // load XSec Integrator - fXSecIntegrator = - dynamic_cast (this->SubAlg("XSec-Integrator")); - assert(fXSecIntegrator); + } //____________________________________________________________________________ diff --git a/src/Physics/NuElectron/XSection/IMDAnnihilationPXSec.h b/src/Physics/NuElectron/XSection/IMDAnnihilationPXSec.h index b3a7f9b03..6f83e282e 100644 --- a/src/Physics/NuElectron/XSection/IMDAnnihilationPXSec.h +++ b/src/Physics/NuElectron/XSection/IMDAnnihilationPXSec.h @@ -6,13 +6,15 @@ \brief nuebar + e- -> mu- + numubar [CC] scattering differential cross section \n - Is a concrete implementation of the XSecAlgorithmI interface. \n + Is a concrete implementation of the PXSecOnElectron interface. \n \ref W.J.Marciano and Z.Parsa, Neutrino-electron scattering theory, J.Phys.G: Nucl.Part.Phys. 29 (2003) 2629-2645 \author Rosen Matev (r.matev@gmail.com) + B. Carlson implemented changes to interface PXSecOnElectron + \created October 3, 2011 \cpright Copyright (c) 2003-2023, The GENIE Collaboration @@ -23,36 +25,33 @@ #ifndef _IMD_ANNIHILATION_PXSEC_H_ #define _IMD_ANNIHILATION_PXSEC_H_ -#include "Framework/EventGen/XSecAlgorithmI.h" +#include "Physics/NuElectron/XSection/PXSecOnElectron.h" namespace genie { class IntegratorI; class XSecIntegratorI; -class IMDAnnihilationPXSec : public XSecAlgorithmI { +class IMDAnnihilationPXSec : public PXSecOnElectron { public: IMDAnnihilationPXSec(); IMDAnnihilationPXSec(string config); virtual ~IMDAnnihilationPXSec(); - //-- XSecAlgorithmI interface implementation - double XSec (const Interaction * i, KinePhaseSpace_t k) const; - double Integral (const Interaction * i) const; - bool ValidProcess (const Interaction * i) const; - bool ValidKinematics (const Interaction * i) const; + //-- PXSecOnElectron interface implementation + double XSec (const Interaction * i, KinePhaseSpace_t k) const; + bool ValidProcess (const Interaction * i) const; + bool ValidKinematics (const Interaction * i) const; //-- overload the Algorithm::Configure() methods to load private data // members from configuration options - void Configure(const Registry & config); - void Configure(string config); + void Configure(const Registry & config); + void Configure(string config); private: void LoadConfig (void); - const XSecIntegratorI * fXSecIntegrator; - }; } // genie namespace diff --git a/src/Physics/NuElectron/XSection/IMDXSec.cxx b/src/Physics/NuElectron/XSection/IMDXSec.cxx deleted file mode 100644 index b836d05bd..000000000 --- a/src/Physics/NuElectron/XSection/IMDXSec.cxx +++ /dev/null @@ -1,97 +0,0 @@ -//____________________________________________________________________________ -/* - Copyright (c) 2003-2023, The GENIE Collaboration - For the full text of the license visit http://copyright.genie-mc.org - - Costas Andreopoulos - University of Liverpool & STFC Rutherford Appleton Laboratory -*/ -//____________________________________________________________________________ - -#include -#include -#include - -#include "Framework/Conventions/GBuild.h" -#include "Framework/Conventions/Constants.h" -#include "Framework/Conventions/Units.h" -#include "Framework/Conventions/RefFrame.h" -#include "Physics/NuElectron/XSection/IMDXSec.h" -#include "Physics/XSectionIntegration/GSLXSecFunc.h" -#include "Framework/Messenger/Messenger.h" -#include "Framework/Numerical/GSLUtils.h" - -using namespace genie; -using namespace genie::constants; - -//____________________________________________________________________________ -IMDXSec::IMDXSec() : -XSecIntegratorI("genie::IMDXSec") -{ - -} -//____________________________________________________________________________ -IMDXSec::IMDXSec(string config) : -XSecIntegratorI("genie::IMDXSec", config) -{ - -} -//____________________________________________________________________________ -IMDXSec::~IMDXSec() -{ - -} -//____________________________________________________________________________ -double IMDXSec::Integrate( - const XSecAlgorithmI * model, const Interaction * in) const -{ - if(! model->ValidProcess(in) ) return 0.; - - const KPhaseSpace & kps = in->PhaseSpace(); - if(!kps.IsAboveThreshold()) { - LOG("IMDXSec", pDEBUG) << "*** below energy threshold"; - return 0; - } - Range1D_t yl = kps.Limits(kKVy); - - Interaction * interaction = new Interaction(*in); - interaction->SetBit(kISkipProcessChk); - //interaction->SetBit(kISkipKinematicChk); - - double abstol = 1; // We mostly care about relative tolerance - ROOT::Math::IBaseFunctionOneDim * func = new utils::gsl::dXSec_dy_E(model, interaction); - ROOT::Math::IntegrationOneDim::Type ig_type = utils::gsl::Integration1DimTypeFromString(fGSLIntgType); - ROOT::Math::Integrator ig(*func,ig_type,abstol,fGSLRelTol,fGSLMaxEval); - double xsec = ig.Integral(yl.min, yl.max) * (1E-38 * units::cm2); - - //LOG("IMDXSec", pDEBUG) << "XSec[IMD] (E = " << E << ") = " << xsec; - - delete interaction; - delete func; - - return xsec; -} -//____________________________________________________________________________ -void IMDXSec::Configure(const Registry & config) -{ - Algorithm::Configure(config); - this->LoadConfig(); -} -//____________________________________________________________________________ -void IMDXSec::Configure(string config) -{ - Algorithm::Configure(config); - this->LoadConfig(); -} -//____________________________________________________________________________ -void IMDXSec::LoadConfig(void) -{ - // Get GSL integration type & relative tolerance - GetParamDef( "gsl-integration-type", fGSLIntgType, string( "adaptive" ) ) ; - GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 1E-4 ) ; - int max; - GetParamDef( "gsl-max-eval", max, 100000 ) ; - fGSLMaxEval = (unsigned int) max ; - -} -//____________________________________________________________________________ diff --git a/src/Physics/NuElectron/XSection/IMDXSec.h b/src/Physics/NuElectron/XSection/IMDXSec.h deleted file mode 100644 index 8e66faa45..000000000 --- a/src/Physics/NuElectron/XSection/IMDXSec.h +++ /dev/null @@ -1,45 +0,0 @@ -//____________________________________________________________________________ -/*! - -\class genie::IMDXSec - -\brief Computes the Inverse Muon Decay cross section. - -\author Costas Andreopoulos - University of Liverpool & STFC Rutherford Appleton Laboratory - -\created Fabruary 14, 2005 - -\cpright Copyright (c) 2003-2023, The GENIE Collaboration - For the full text of the license visit http://copyright.genie-mc.org -*/ -//____________________________________________________________________________ - -#ifndef _IMD_XSEC_H_ -#define _IMD_XSEC_H_ - -#include "Physics/XSectionIntegration/XSecIntegratorI.h" - -namespace genie { - -class IMDXSec : public XSecIntegratorI { - -public: - IMDXSec(); - IMDXSec(string config); - virtual ~IMDXSec(); - - //! XSecIntegratorI interface implementation - double Integrate(const XSecAlgorithmI * model, const Interaction * i) 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); -}; - -} // genie namespace -#endif // _IMD_XSEC_H_ diff --git a/src/Physics/NuElectron/XSection/LinkDef.h b/src/Physics/NuElectron/XSection/LinkDef.h index 28f9ccc7d..c7813029e 100644 --- a/src/Physics/NuElectron/XSection/LinkDef.h +++ b/src/Physics/NuElectron/XSection/LinkDef.h @@ -6,12 +6,12 @@ #pragma link C++ namespace genie; +#pragma link C++ class genie::PXSecOnElectron; + #pragma link C++ class genie::IMDAnnihilationPXSec; #pragma link C++ class genie::BardinIMDRadCorPXSec; #pragma link C++ class genie::NuElectronPXSec; -#pragma link C++ class genie::IMDXSec; - -#pragma link C++ class genie::NuElectronXSec; +#pragma link C++ class genie::XSecOnElectron; #endif diff --git a/src/Physics/NuElectron/XSection/NuElectronPXSec.cxx b/src/Physics/NuElectron/XSection/NuElectronPXSec.cxx index b39fbd418..a1d42b02b 100644 --- a/src/Physics/NuElectron/XSection/NuElectronPXSec.cxx +++ b/src/Physics/NuElectron/XSection/NuElectronPXSec.cxx @@ -5,19 +5,28 @@ Costas Andreopoulos University of Liverpool & STFC Rutherford Appleton Laboratory + + Changes required to implement the Electron Velocity module + were installed by Brinden Carlson (Univ. of Florida) */ //____________________________________________________________________________ #include "Framework/Algorithm/AlgConfigPool.h" #include "Framework/Conventions/GBuild.h" #include "Framework/Conventions/Controls.h" -#include "Physics/XSectionIntegration/XSecIntegratorI.h" #include "Framework/Conventions/Constants.h" #include "Framework/Conventions/RefFrame.h" #include "Physics/NuElectron/XSection/NuElectronPXSec.h" #include "Framework/Messenger/Messenger.h" #include "Framework/ParticleData/PDGUtils.h" #include "Framework/Utils/KineUtils.h" +#include "Physics/NuElectron/XSection/PXSecOnElectron.h" + + +#include "TFile.h" +#include "TGraph.h" +#include +#include using namespace genie; using namespace genie::constants; @@ -25,13 +34,13 @@ using namespace genie::controls; //____________________________________________________________________________ NuElectronPXSec::NuElectronPXSec() : -XSecAlgorithmI("genie::NuElectronPXSec") +PXSecOnElectron::PXSecOnElectron("genie::NuElectronPXSec","Default") { } //____________________________________________________________________________ NuElectronPXSec::NuElectronPXSec(string config) : -XSecAlgorithmI("genie::NuElectronPXSec", config) +PXSecOnElectron::PXSecOnElectron("genie::NuElectronPXSec", config) { } @@ -52,7 +61,8 @@ double NuElectronPXSec::XSec( const Kinematics & kinematics = interaction -> Kine(); const ProcessInfo & proc_info = interaction -> ProcInfo(); - double Ev = init_state.ProbeE(kRfLab); + double Ev = init_state.ProbeE(kRfHitElRest); //Electron rest frame + double me = kElectronMass; double y = kinematics.y(); double A = kGF2*2*me*Ev/kPi; @@ -132,33 +142,15 @@ double NuElectronPXSec::XSec( return xsec; } //____________________________________________________________________________ -double NuElectronPXSec::Integral(const Interaction * interaction) const -{ - double xsec = fXSecIntegrator->Integrate(this,interaction); - return xsec; -} -//____________________________________________________________________________ -bool NuElectronPXSec::ValidProcess(const Interaction * interaction) const -{ - if(interaction->TestBit(kISkipProcessChk)) return true; - return true; -} -//____________________________________________________________________________ -bool NuElectronPXSec::ValidKinematics(const Interaction* interaction) const -{ - if(interaction->TestBit(kISkipKinematicChk)) return true; - return true; -} -//____________________________________________________________________________ void NuElectronPXSec::Configure(const Registry & config) { - Algorithm::Configure(config); + PXSecOnElectron::Configure(config); this->LoadConfig(); } //____________________________________________________________________________ void NuElectronPXSec::Configure(string config) { - Algorithm::Configure(config); + PXSecOnElectron::Configure(config); this->LoadConfig(); } //____________________________________________________________________________ @@ -169,10 +161,10 @@ void NuElectronPXSec::LoadConfig(void) GetParam( "WeinbergAngle", thw ) ; fSin28w = TMath::Power(TMath::Sin(thw), 2); fSin48w = TMath::Power(TMath::Sin(thw), 4); - - // load XSec Integrator - fXSecIntegrator = - dynamic_cast (this->SubAlg("XSec-Integrator")); - assert(fXSecIntegrator); + + if (!fElectronVelocity) { + LOG("NuElectronPXSec", pDEBUG) + << "fElectronVelocity is not initialized correctly."; + } } //____________________________________________________________________________ diff --git a/src/Physics/NuElectron/XSection/NuElectronPXSec.h b/src/Physics/NuElectron/XSection/NuElectronPXSec.h index e39775d5a..5fa628a91 100644 --- a/src/Physics/NuElectron/XSection/NuElectronPXSec.h +++ b/src/Physics/NuElectron/XSection/NuElectronPXSec.h @@ -18,6 +18,9 @@ \author Costas Andreopoulos University of Liverpool & STFC Rutherford Appleton Laboratory + Changes required to implement the Electron Velocity module + were installed by Brinden Carlson (Univ. of Florida) + \created February 10, 2006 \cpright Copyright (c) 2003-2023, The GENIE Collaboration @@ -28,14 +31,11 @@ #ifndef _NU_ELECTRON_PARTIAL_XSEC_H_ #define _NU_ELECTRON_PARTIAL_XSEC_H_ -#include "Framework/EventGen/XSecAlgorithmI.h" +#include "Physics/NuElectron/XSection/PXSecOnElectron.h" namespace genie { -class IntegratorI; -class XSecIntegratorI; - -class NuElectronPXSec : public XSecAlgorithmI { +class NuElectronPXSec : public PXSecOnElectron { public: NuElectronPXSec(); @@ -44,9 +44,6 @@ class NuElectronPXSec : public XSecAlgorithmI { //-- XSecAlgorithmI interface implementation double XSec (const Interaction * i, KinePhaseSpace_t k) const; - double Integral (const Interaction * i) const; - bool ValidProcess (const Interaction * i) const; - bool ValidKinematics (const Interaction * i) const; //-- overload the Algorithm::Configure() methods to load private data // members from configuration options @@ -56,7 +53,7 @@ class NuElectronPXSec : public XSecAlgorithmI { private: void LoadConfig (void); - const XSecIntegratorI * fXSecIntegrator; + //const ElectronVelocity * fElectronVelocity; double fSin28w; // sin^2(theta-weinberg) double fSin48w; diff --git a/src/Physics/NuElectron/XSection/PXSecOnElectron.cxx b/src/Physics/NuElectron/XSection/PXSecOnElectron.cxx new file mode 100644 index 000000000..86bf8e53c --- /dev/null +++ b/src/Physics/NuElectron/XSection/PXSecOnElectron.cxx @@ -0,0 +1,118 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Brinden Carlson bcarlson1@ufl.edu + University of Florida +*/ +//____________________________________________________________________________ + +#include "Framework/Algorithm/AlgConfigPool.h" +#include "Framework/Conventions/GBuild.h" +#include "Framework/Conventions/Controls.h" +#include "Physics/XSectionIntegration/XSecIntegratorI.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/Conventions/RefFrame.h" +#include "Physics/NuElectron/XSection/PXSecOnElectron.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Framework/Utils/KineUtils.h" +#include "Physics/NuclearState/ElectronVelocity.h" + +#include "TFile.h" +#include "TGraph.h" +#include +#include + +using namespace genie; +using namespace genie::constants; +using namespace genie::controls; + +//____________________________________________________________________________ +PXSecOnElectron::PXSecOnElectron(string name, string config) : +XSecAlgorithmI(name, config) +{ + +} +//____________________________________________________________________________ +PXSecOnElectron::PXSecOnElectron() : +XSecAlgorithmI("genie::PXSecOnElectron") +{ + +} +//____________________________________________________________________________ +PXSecOnElectron::~PXSecOnElectron() +{ + +} +//____________________________________________________________________________ +double PXSecOnElectron::Integral(const Interaction * interaction) const +{ + double xsec_sum = 0; + double xsec_sum2 = 0; + int NInt = 0; //Count number of integrations + do{ + NInt++; + Interaction in_curr(*interaction); //Copy interaction object + fElectronVelocity->InitializeVelocity(in_curr); //Modify interaction to give electron random velocity from selected distribution + double xsec = fXSecIntegrator->Integrate(this,&in_curr); // In ele at rest + //get beta - comps orthogonal to beam x,y + //scale = sqrt(1-b_t^2) + TVector3 beta = in_curr.InitState().Tgt().HitEleP4().BoostVector(); // beta + double beta_tangent = sqrt(TMath::Power(beta[0],2)+TMath::Power(beta[1],2)); //Component tangential to beam + xsec *= sqrt(1-TMath::Power(beta_tangent,2)); //Correct for lorentz factor + xsec_sum+=xsec; + xsec_sum2+=TMath::Power(xsec,2); + double xsec_mean = xsec_sum/NInt; + //var = (sum(xi^2)/N-xsec_mean^2) + //rel_err = sigma/sqrt(n)*mean + double xsec_err = sqrt((xsec_sum2/NInt-TMath::Power(xsec_mean,2))/NInt)/xsec_mean; + if (NInt > 1 && xsec_err < fErrTolerance){break;} //Break condition for dipping below set tolerance + } + while ( NInt < fNIntegration); + double xsec_avg = xsec_sum/NInt; + return xsec_avg; +} +//____________________________________________________________________________ +bool PXSecOnElectron::ValidProcess(const Interaction * interaction) const +{ + if(interaction->TestBit(kISkipProcessChk)) return true; + return true; +} +//____________________________________________________________________________ +bool PXSecOnElectron::ValidKinematics(const Interaction* interaction) const +{ + if(interaction->TestBit(kISkipKinematicChk)) return true; + return true; +} +//____________________________________________________________________________ +void PXSecOnElectron::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void PXSecOnElectron::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void PXSecOnElectron::LoadConfig(void) +{ + //Integration parameters + GetParam( "N-Integration-Samples", fNIntegration ) ; // + GetParam( "NuE-XSecRelError" , fErrTolerance ) ; // + + // load XSec Integrator + fXSecIntegrator = + dynamic_cast (this->SubAlg("XSec-Integrator")); // + fElectronVelocity = + dynamic_cast (this->SubAlg("Electron-Velocity")); // + if (!fElectronVelocity) { + LOG("PXSecOnElectron", pDEBUG) + << "fElectronVelocity is not initialized correctly."; + } +} +//____________________________________________________________________________ diff --git a/src/Physics/NuElectron/XSection/PXSecOnElectron.h b/src/Physics/NuElectron/XSection/PXSecOnElectron.h new file mode 100644 index 000000000..64613b13b --- /dev/null +++ b/src/Physics/NuElectron/XSection/PXSecOnElectron.h @@ -0,0 +1,71 @@ +//____________________________________________________________________________ +/*! + +\class genie::PXSecOnElectron + +\brief nu/nubar + e- scattering differential cross section \n + The cross section algorithm interfaces nu+e cross sections + Final states are handled by sub modules + + Is a concrete implementation of the XSecAlgorithmI interface. \n + +\ref W.J.Marciano and Z.Parsa, Neutrino-electron scattering theory, + J.Phys.G: Nucl.Part.Phys. 29 (2003) 2629-2645 +\ref Oleksandr Tomalak and Richard J. Hill, Theory of elastic neutrino-electron scattering, + Phys. Rev. D 101, 033006 – Published 21 February 2020 + +\author Brinden Carlson bcarlson1@ufl.edu + University of Florida + +\created June 1, 2023 + +\cpright Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org +*/ +//____________________________________________________________________________ + +#ifndef _NU_ON_ELECTRON_PARTIAL_XSEC_H_ +#define _NU_ON_ELECTRON_PARTIAL_XSEC_H_ + +#include "Framework/EventGen/XSecAlgorithmI.h" +#include "Physics/NuclearState/ElectronVelocity.h" + +namespace genie { + +class IntegratorI; +class XSecIntegratorI; + +class PXSecOnElectron : public XSecAlgorithmI { + +public: + PXSecOnElectron(); + PXSecOnElectron(string name,string config); + virtual ~PXSecOnElectron(); + + //-- XSecAlgorithmI interface implementation + virtual double XSec (const Interaction * i, KinePhaseSpace_t k) const = 0; + double Integral (const Interaction * i) const; + bool ValidProcess (const Interaction * i) const; + bool ValidKinematics (const Interaction * i) const; + + //-- overload the Algorithm::Configure() methods to load private data + // members from configuration options + void Configure(const Registry & config); + void Configure(string config); + +private: + + +protected: + void LoadConfig (void); + + const XSecIntegratorI * fXSecIntegrator; + const ElectronVelocity * fElectronVelocity; + + int fNIntegration; //Max number of integration samples + double fErrTolerance; //Error tolerance acceptable before returning cross section average +}; + +} // genie namespace +#endif // _NU_ON_ELECTRON_PARTIAL_XSEC_H_ + diff --git a/src/Physics/NuElectron/XSection/NuElectronXSec.cxx b/src/Physics/NuElectron/XSection/XSecOnElectron.cxx similarity index 86% rename from src/Physics/NuElectron/XSection/NuElectronXSec.cxx rename to src/Physics/NuElectron/XSection/XSecOnElectron.cxx index 5148ea85a..ec934b512 100644 --- a/src/Physics/NuElectron/XSection/NuElectronXSec.cxx +++ b/src/Physics/NuElectron/XSection/XSecOnElectron.cxx @@ -16,7 +16,7 @@ #include "Framework/Conventions/Constants.h" #include "Framework/Conventions/Units.h" #include "Framework/Conventions/RefFrame.h" -#include "Physics/NuElectron/XSection/NuElectronXSec.h" +#include "Physics/NuElectron/XSection/XSecOnElectron.h" #include "Physics/XSectionIntegration/GSLXSecFunc.h" #include "Framework/Messenger/Messenger.h" #include "Framework/Numerical/GSLUtils.h" @@ -25,24 +25,24 @@ using namespace genie; using namespace genie::constants; //____________________________________________________________________________ -NuElectronXSec::NuElectronXSec() : -XSecIntegratorI("genie::NuElectronXSec") +XSecOnElectron::XSecOnElectron() : +XSecIntegratorI("genie::XSecOnElectron") { } //____________________________________________________________________________ -NuElectronXSec::NuElectronXSec(string config) : -XSecIntegratorI("genie::NuElectronXSec", config) +XSecOnElectron::XSecOnElectron(string config) : +XSecIntegratorI("genie::XSecOnElectron", config) { } //____________________________________________________________________________ -NuElectronXSec::~NuElectronXSec() +XSecOnElectron::~XSecOnElectron() { } //____________________________________________________________________________ -double NuElectronXSec::Integrate( +double XSecOnElectron::Integrate( const XSecAlgorithmI * model, const Interaction * in) const { if(! model->ValidProcess(in) ) return 0.; @@ -71,22 +71,23 @@ double NuElectronXSec::Integrate( delete interaction; delete func; + return xsec; } //____________________________________________________________________________ -void NuElectronXSec::Configure(const Registry & config) +void XSecOnElectron::Configure(const Registry & config) { Algorithm::Configure(config); this->LoadConfig(); } //____________________________________________________________________________ -void NuElectronXSec::Configure(string config) +void XSecOnElectron::Configure(string config) { Algorithm::Configure(config); this->LoadConfig(); } //____________________________________________________________________________ -void NuElectronXSec::LoadConfig(void) +void XSecOnElectron::LoadConfig(void) { // Get GSL integration type & relative tolerance GetParamDef( "gsl-integration-type", fGSLIntgType, string( "adaptive" ) ) ; diff --git a/src/Physics/NuElectron/XSection/NuElectronXSec.h b/src/Physics/NuElectron/XSection/XSecOnElectron.h similarity index 86% rename from src/Physics/NuElectron/XSection/NuElectronXSec.h rename to src/Physics/NuElectron/XSection/XSecOnElectron.h index 156aec5e3..49396b802 100644 --- a/src/Physics/NuElectron/XSection/NuElectronXSec.h +++ b/src/Physics/NuElectron/XSection/XSecOnElectron.h @@ -1,7 +1,7 @@ //____________________________________________________________________________ /*! -\class genie::NuElectronXSec +\class genie::XSecOnElectron \brief nu/nubar + e- scattering cross section. Integrates the loaded differential cross section model. An analytical cross section @@ -19,6 +19,9 @@ \author Costas Andreopoulos University of Liverpool & STFC Rutherford Appleton Laboratory + B. Carlson change name to reflect that cross section applies + to all nu+e interactions + \created February 10, 2006 \cpright Copyright (c) 2003-2023, The GENIE Collaboration @@ -33,12 +36,12 @@ namespace genie { -class NuElectronXSec : public XSecIntegratorI { +class XSecOnElectron : public XSecIntegratorI { public: - NuElectronXSec(); - NuElectronXSec(string config); - virtual ~NuElectronXSec(); + XSecOnElectron(); + XSecOnElectron(string config); + virtual ~XSecOnElectron(); //! XSecIntegratorI interface implementation double Integrate(const XSecAlgorithmI * model, const Interaction * i) const; diff --git a/src/Physics/NuclearState/BohrElectronVelocity.cxx b/src/Physics/NuclearState/BohrElectronVelocity.cxx new file mode 100644 index 000000000..928cb83d1 --- /dev/null +++ b/src/Physics/NuclearState/BohrElectronVelocity.cxx @@ -0,0 +1,144 @@ +///____________________________________________________________________________ +/* + Copyright (c) 2003-2022, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + \brief It visits the event record & computes a Bohr Velocity for + initial state electrons bound in coloumb potential. + + \author Brinden Carlson + University of Florida & Fermilab + + \created December 5, 2022 + + For the class documentation see the corresponding header file. +*/ +//____________________________________________________________________________ + +#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/ElectronVelocity.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/BohrElectronVelocity.h" +#include "Physics/NuclearState/ElectronVelocity.h" + +#include +#include + +using namespace genie; +using namespace genie::constants; + +//___________________________________________________________________________ +BohrElectronVelocity::~BohrElectronVelocity() +{ + +} +BohrElectronVelocity::BohrElectronVelocity(const string & config) : +ElectronVelocity::ElectronVelocity("genie::BohrElectronVelocity", config) +{ + +} +BohrElectronVelocity::BohrElectronVelocity() : +ElectronVelocity::ElectronVelocity("genie::BohrElectronVelocity","Default") +{ + +} + +//___________________________________________________________________________ +void BohrElectronVelocity::InitializeVelocity(Interaction & interaction) const{ + InitialState * init_state = interaction.InitStatePtr(); + Target * tgt = init_state -> TgtPtr(); + + unsigned int fZ = tgt->Z(); //Get Z value + double nucleus_mass = tgt->Mass(); + double mu = reduced_mass(nucleus_mass); //Reduced mass + double v = random_bohr_velocity(fZ,mu); //bohr velocity + TVector3 v3 = randomize_direction_sphere(v); //Get spherically uniform random direciton + double gamma = 1/sqrt(1-v3.Mag2()); //Get boost + //Set 3 momentum + auto p3 = kElectronMass*gamma*v3; + //-- update the electron 4p + TLorentzVector * p4 = tgt->HitEleP4Ptr(); //Initialize 4 momentum pointer + p4->SetVectM( p3, kElectronMass); +} +//___________________________________________________________________________ + +//____________________________________________________________________________ +double BohrElectronVelocity::reduced_mass(double nucleus_mass) const{ + //Return reduced mass of electron + return kElectronMass*nucleus_mass/(kElectronMass+nucleus_mass); +} + +unsigned BohrElectronVelocity::random_n(unsigned int fZ) const{ + //Isotope minimum z is 118 in accordance with maximum possible in GENIE + std::array fnprobs {2,10,28,60,110,118}; //Cumulative Probability dist. + + //Get random electron orbital from atomic number Z + RandomGen * rnd = RandomGen::Instance(); //Load seed + + for (unsigned int i = 0; iRndDec().Rndm(); + unsigned int n = 0; + unsigned int sel_n = 0; + do { + sel_n = n; + } while (x > fnprobs[n++]); + + return sel_n+1; //Plus 1 to get to n +} + +double BohrElectronVelocity::bohr_velocity(unsigned int fn, unsigned int fZ, double mu) const +{ + return fZ*kAem/fn * kElectronMass/mu; +} + +double BohrElectronVelocity::random_bohr_velocity(unsigned int fZ, double mu) const{ + //Get random bohr velocity from n distribution + unsigned int fn = random_n(fZ); + return bohr_velocity(fn,fZ,mu); +} + +TVector3 BohrElectronVelocity::randomize_direction_sphere(double v) const{ + RandomGen * rnd = RandomGen::Instance(); //Load seed + TRandom3 gen = rnd->RndGen(); //Get generator + double costheta = gen.Uniform(-1,1); //Polar + double phi = gen.Uniform(0,2*kPi); //Azimuthal + double sintheta = sqrt( 1 - costheta*costheta); // which is always positive because the angle is in the 0 -> pi range! it cannot be negative + return TVector3(v*sintheta*TMath::Cos(phi), + v*sintheta*TMath::Sin(phi), + v*costheta); //Return vector with magnitude v and random direction +} + + +//___________________________________________________________________________ \ No newline at end of file diff --git a/src/Physics/NuclearState/BohrElectronVelocity.h b/src/Physics/NuclearState/BohrElectronVelocity.h new file mode 100644 index 000000000..7b94438a0 --- /dev/null +++ b/src/Physics/NuclearState/BohrElectronVelocity.h @@ -0,0 +1,45 @@ +//____________________________________________________________________________ +/*! + +\class genie::BohrElectronVelocity + +\brief It visits the event record & computes a Bohr Velocity for + initial state electrons bound in coloumb potential. + Is a concrete implementation of the ElectronVelocity interface. + +\author Brinden Carlson + University of Florida & Fermilab + +\created December 5, 2022 + +\cpright Copyright (c) 2003-2022, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + +*/ +//____________________________________________________________________________ + +// #ifndef _BOHR_ELECTRON_VELOCITY_H_ +// #define _BOHR_ELECTRON_VELOCITY_H_ + +#include "Physics/NuclearState/ElectronVelocity.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/Interaction/Target.h" +#include "Framework/Interaction/Interaction.h" + +namespace genie { + +class BohrElectronVelocity : public ElectronVelocity { +public: + void InitializeVelocity(Interaction & interaction) const final; + BohrElectronVelocity(); + BohrElectronVelocity(const string & config); + ~BohrElectronVelocity(); + +private: + double reduced_mass(double nucleus_mass) const; //Reduced mass calculation + double bohr_velocity(unsigned int fn, unsigned int fZ, double mu) const; //Bohr velocity + unsigned int random_n(unsigned int fZ) const; //Return random energy level from n_dist + double random_bohr_velocity(unsigned int fZ, double mu) const; //Generate random n, then calculate velocity from there + TVector3 randomize_direction_sphere(double v) const; +}; +} \ No newline at end of file diff --git a/src/Physics/NuclearState/ElectronVelocity.cxx b/src/Physics/NuclearState/ElectronVelocity.cxx new file mode 100644 index 000000000..f324ea2af --- /dev/null +++ b/src/Physics/NuclearState/ElectronVelocity.cxx @@ -0,0 +1,97 @@ +///____________________________________________________________________________ +/* + Copyright (c) 2003-2022, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + \brief It visits the event record & samples a velocity for + initial state electrons from a velocity distribution. + Is a concrete implementation of the EventRecordVisitorI interface. + + \author Brinden Carlson + University of Florida & Fermilab + + \created December 5, 2022 + + For the class documentation see the corresponding header file. +*/ +//____________________________________________________________________________ + +#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/ElectronVelocity.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 +#include + +using namespace genie; + +//___________________________________________________________________________ +ElectronVelocity::ElectronVelocity(const string & name, const string & config) : +EventRecordVisitorI(name, config) +{ + +} +//___________________________________________________________________________ +ElectronVelocity::~ElectronVelocity() +{ + +} +//___________________________________________________________________________ +ElectronVelocity::ElectronVelocity() +{ + +} +//___________________________________________________________________________ +void ElectronVelocity::ProcessEventRecord(GHepRecord * evrec) const +{ + // skip if not a electron target + if(!evrec->Summary()->ProcInfo().IsElectronScattering()) return; + + // give electron initial momentum + this->InitializeVelocity(*evrec->Summary()); + + //Update event record + GHepParticle * electron = evrec->HitElectron(); + electron->SetMomentum(*evrec->Summary()->InitStatePtr()->TgtPtr()->HitEleP4Ptr()); + +} +//___________________________________________________________________________ +void ElectronVelocity::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void ElectronVelocity::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} + diff --git a/src/Physics/NuclearState/ElectronVelocity.h b/src/Physics/NuclearState/ElectronVelocity.h new file mode 100644 index 000000000..cad2967a3 --- /dev/null +++ b/src/Physics/NuclearState/ElectronVelocity.h @@ -0,0 +1,59 @@ +//____________________________________________________________________________ +/*! + +\class genie::ElectronVelocity + +\brief Interface used to visit the event record & samples a velocity for + initial state electrons from a velocity distribution. + Modifies interaction object to modify electron velocity and used to + integrate over velocity distribution when calculating the cross section + Is a concrete implementation of the EventRecordVisitorI interface. + + \author Brinden Carlson + University of Florida & Fermilab + + \created December 5, 2022 + +\cpright Copyright (c) 2003-2022, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + +*/ +//____________________________________________________________________________ + +#ifndef _ELECTRON_VELOCITY_H_ +#define _ELECTRON_VELOCITY_H_ + +#include "Framework/EventGen/EventRecordVisitorI.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/Interaction/Target.h" +#include "Framework/Interaction/Interaction.h" + +namespace genie { + +class ElectronVelocity : public EventRecordVisitorI { + +public : + virtual ~ElectronVelocity(); + ElectronVelocity(); + + //-- implement the EventRecordVisitorI interface + void ProcessEventRecord(GHepRecord * event_rec) const override; + + //-- overload the Algorithm::Configure() methods to load private data + // members from configuration options + void Configure(const Registry & config) override; + void Configure(string config) override; + + //Make public to use in NuElectronPXsec + virtual void InitializeVelocity(Interaction & interaction) const = 0; //Give initial velocity + +private: + +protected: + ElectronVelocity(const string & name, const string & config); + virtual void LoadConfig (void){;} + +}; + +} // genie namespace +#endif // _FERMI_MOVER_H_ diff --git a/src/Physics/NuclearState/LinkDef.h b/src/Physics/NuclearState/LinkDef.h index c9207f1ca..b981b5c92 100644 --- a/src/Physics/NuclearState/LinkDef.h +++ b/src/Physics/NuclearState/LinkDef.h @@ -26,4 +26,7 @@ #pragma link C++ class genie::SecondNucleonEmissionI; #pragma link C++ class genie::SpectralFunction2p2h; +#pragma link C++ class genie::BohrElectronVelocity; +#pragma link C++ class genie::StaticElectronVelocity; + #endif diff --git a/src/Physics/NuclearState/StaticElectronVelocity.cxx b/src/Physics/NuclearState/StaticElectronVelocity.cxx new file mode 100644 index 000000000..b7be7bc83 --- /dev/null +++ b/src/Physics/NuclearState/StaticElectronVelocity.cxx @@ -0,0 +1,85 @@ +///____________________________________________________________________________ +/* + Copyright (c) 2003-2022, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + \brief It visits the event record & initializes a static velocity for + initial state electron. + +\author Brinden Carlson + University of Florida & Fermilab + +\created December 5, 2022 + + For the class documentation see the corresponding header file. +*/ +//____________________________________________________________________________ + +#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/ElectronVelocity.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/StaticElectronVelocity.h" +#include "Physics/NuclearState/ElectronVelocity.h" + +#include +#include + +using namespace genie; +using namespace genie::constants; +//using namespace std; + +//___________________________________________________________________________ +StaticElectronVelocity::~StaticElectronVelocity() +{ + +} +StaticElectronVelocity::StaticElectronVelocity(const string & config) : +ElectronVelocity::ElectronVelocity("genie::StaticElectronVelocity", config) +{ + +} +StaticElectronVelocity::StaticElectronVelocity() : +ElectronVelocity::ElectronVelocity("genie::StaticElectronVelocity","Default") +{ + +} + +//___________________________________________________________________________ +void StaticElectronVelocity::InitializeVelocity(Interaction & interaction) const{ + InitialState * init_state = interaction.InitStatePtr(); + Target * tgt = init_state -> TgtPtr(); + + TLorentzVector * p4 = tgt->HitEleP4Ptr(); //Initialize 4 momentum pointer + //These should be initialized like this by just in case + p4->SetVectM(TVector3(), kElectronMass); +} +//___________________________________________________________________________ + +//____________________________________________________________________________ diff --git a/src/Physics/NuclearState/StaticElectronVelocity.h b/src/Physics/NuclearState/StaticElectronVelocity.h new file mode 100644 index 000000000..d2a1e0f55 --- /dev/null +++ b/src/Physics/NuclearState/StaticElectronVelocity.h @@ -0,0 +1,36 @@ +//____________________________________________________________________________ +/*! + +\class genie::StaticElectronVelocity + +\brief It visits the event record & initializes a static velocity for + initial state electron. + +\author Brinden Carlson + University of Florida & Fermilab + +\created December 5, 2022 + +\cpright Copyright (c) 2003-2022, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + +*/ +//____________________________________________________________________________ + +#include "Physics/NuclearState/ElectronVelocity.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/Interaction/Target.h" +#include "Framework/Interaction/Interaction.h" + +namespace genie { + +class StaticElectronVelocity : public ElectronVelocity { +public: + void InitializeVelocity(Interaction & interaction) const final; + + StaticElectronVelocity(); + StaticElectronVelocity(const string & config); + ~StaticElectronVelocity(); + +}; +} \ No newline at end of file