Electron velocity distribution implementation

config/BohrElectronVelocity.xml

@@ -0,0 +1,17 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+
+<!--
+Configuration for the Bohr Electron Velocity
+
+Configurable Parameters:
+.............................................................................................................
+Name                          Type    Optional   Comment                                           Default
+.............................................................................................................
+
+-->
+
+<alg_conf>
+  <param_set name="Default">
+  </param_set>
+
+</alg_conf>

config/EventGenerator.xml

@@ -455,12 +455,13 @@ XSecModel                   alg     Yes  Cross section model used at the thread
 
   <param_set name="NUE-EL">
      <param type="string" name="VldContext"> </param>
-     <param type="int"    name="NModules">   5                                          </param>
+     <param type="int"    name="NModules">   6                                          </param>
      <param type="alg"    name="Module-0">   genie::InitialStateAppender/Default        </param>
      <param type="alg"    name="Module-1">   genie::VertexGenerator/Default             </param>
-     <param type="alg"    name="Module-2">   genie::NuEKinematicsGenerator/Default      </param>
-     <param type="alg"    name="Module-3">   genie::NuEPrimaryLeptonGenerator/Default   </param>
-     <param type="alg"    name="Module-4">   genie::NuETargetRemnantGenerator/Default   </param>
+     <param type="alg"    name="Module-2">   genie::BohrElectronVelocity/Default            </param>
+     <param type="alg"    name="Module-3">   genie::NuEKinematicsGenerator/Default      </param>
+     <param type="alg"    name="Module-4">   genie::NuEPrimaryLeptonGenerator/Default   </param>
+     <param type="alg"    name="Module-5">   genie::NuETargetRemnantGenerator/Default   </param>
      <param type="alg"    name="ILstGen">    genie::NuEInteractionListGenerator/NUE-EL  </param>
   </param_set>
 

config/NuElectronPXSec.xml

@@ -7,15 +7,20 @@ Configuration for the NuElectronPXSec xsec algorithm.
 
 Configurable Parameters:
 ....................................................................................
-Name           Type     Optional   Comment                 Default
-WeinbergAngle  double   No                                 CommonParam[WeakInt]
+Name                    Type     Optional   Comment                  Default
+N-Integration-Samples   int      No                                  2000
+WeinbergAngle           double   No                                  CommonParam[WeakInt]
+XSec-Integrator         alg      No                                  genie::NuElectronXSec/Default
+Electron-Velocity       alg      Yes    Integrate over vel dist.     genie::ElectronVelocity/Default
 ....................................................................................
 -->
 
   <param_set name="Default"> 
     <param type="string"  name="CommonParam"> WeakInt  </param>
-
-     <param type="alg"  name="XSec-Integrator"> genie::NuElectronXSec/Default  </param>
+    <param type="int"    name ="N-Integration-Samples" >              2000  </param>
+    <param type="double"    name ="NuE-Integration-Error-Tolerance" >     0.01  </param>
+    <param type="alg"  name="XSec-Integrator"> genie::NuElectronXSec/Default  </param>
+    <param type="alg"  name="Electron-Velocity"> genie::BohrElectronVelocity/Default  </param>
   </param_set>
 
 </alg_conf>

config/StaticElectronVelocity.xml

@@ -0,0 +1,17 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+
+<!--
+Configuration for the Static Electron Velocity
+
+Configurable Parameters:
+.............................................................................................................
+Name                          Type    Optional   Comment                                           Default
+.............................................................................................................
+
+-->
+
+<alg_conf>
+  <param_set name="Default">
+  </param_set>
+
+</alg_conf>

config/master_config.xml

@@ -5,6 +5,8 @@
    <!-- ****** CONFIGURATION FOR EVENT GENERATION MODULES ****** -->
    <config alg="genie::EventGenerator">                  EventGenerator.xml                  </config>
    <config alg="genie::FermiMover">                      FermiMover.xml                      </config>
+   <config alg="genie::BohrElectronVelocity">            BohrElectronVelocity.xml            </config>
+   <config alg="genie::StaticElectronVelocity">          StaticElectronVelocity.xml          </config>
    <config alg="genie::HadronTransporter">               HadronTransporter.xml               </config>
    <config alg="genie::HAIntranuke">                     HAIntranuke.xml                     </config>
    <config alg="genie::HAIntranuke2018">                 HAIntranuke2018.xml                 </config>

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,21 @@ 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
+             //assert(fTgt->HitEleIsSet());
+             if (!fTgt->HitEleIsSet()){
+              return nullptr;
+             }
+             double inv_mass = fTgt->HitEleP4Ptr()->Mag();
+             TLorentzVector * p4;
+             p4->SetVectM(TVector3(),inv_mass);
+             return p4;
+             break; //Do we also need this break?
+       }
        default:
              LOG("Interaction", pERROR) << "Uknown reference frame";
   }
@@ -374,6 +392,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";

src/Framework/Interaction/Target.cxx

@@ -5,6 +5,9 @@
 
  Costas Andreopoulos <constantinos.andreopoulos \at cern.ch>
  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

src/Framework/Interaction/Target.h

@@ -11,6 +11,9 @@
 \author   Costas Andreopoulos <constantinos.andreopoulos \at cern.ch>
           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-2022, 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) //Advance the index?
 };
 
 }      // genie namespace

src/Physics/NuElectron/EventGen/NuEKinematicsGenerator.cxx

@@ -5,6 +5,10 @@
 
  Costas Andreopoulos <constantinos.andreopoulos \at cern.ch>
  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;
 }
 //___________________________________________________________________________

src/Physics/NuElectron/EventGen/NuEPrimaryLeptonGenerator.cxx

@@ -5,6 +5,9 @@
 
  Costas Andreopoulos <constantinos.andreopoulos \at cern.ch>
  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;
+}
+//___________________________________________________________________________