Merge pull request #208 from igor144/smith_moniz_update

Updates of Smith Moniz model; Spline, CacheBranchFx and KineGeneratorWithCache classes; adding new Normalization channel
GENIE-MC · Jul 7, 2023 · 8aba654 · 8aba654
2 parents 70f72d1 + ee3305e
commit 8aba654
Show file tree

Hide file tree

Showing 37 changed files with 1,521 additions and 711 deletions.
diff --git a/config/EventGenerator.xml b/config/EventGenerator.xml
@@ -830,5 +830,15 @@ XSecModel                   alg     Yes  Cross section model used at the thread
   </param_set>
 
 
+  <param_set name="NORM">
+     <param type="string" name="VldContext"> </param>
+     <param type="int"    name="NModules">   3                                                   </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::NormGenerator/Default                        </param>
+     <param type="alg"    name="ILstGen">    genie::NormInteractionListGenerator/Default         </param>
+     <param type="alg"    name="XSecModel@genie::EventGenerator/NORM">  genie::NormXSec/Default  </param>
+  </param_set>
+
 </alg_conf>
 
diff --git a/config/EventGeneratorListAssembler.xml b/config/EventGeneratorListAssembler.xml
@@ -562,5 +562,16 @@ Generator-%d                alg     No
      <param type="alg" name="Generator-0">  genie::EventGenerator/PhotonCOH     </param>
   </param_set>
 
+  <param_set name="NORM"> 
+     <param type="int" name="NGenerators">    1                               </param>
+  <param type="alg" name="Generator-0">    genie::EventGenerator/NORM         </param>
+  </param_set>
+
+  <param_set name="CCQE+NORM"> 
+     <param type="int" name="NGenerators">   2                                     </param>
+     <param type="alg" name="Generator-0">   genie::EventGenerator/QEL-CC          </param>
+     <param type="alg" name="Generator-1">   genie::EventGenerator/NORM            </param>
+  </param_set>
+
 </alg_conf>
 
diff --git a/config/NormGenerator.xml b/config/NormGenerator.xml
@@ -0,0 +1,20 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+
+<!--
+Configuration for the NormCOHGenerator EventRecordVisitorI
+
+Algorithm Configurable Parameters:
+.......................................................................................................................
+Name                     Type    Opt   Comment                                        Default
+.......................................................................................................................
+
+-->
+
+<alg_conf>
+
+  <param_set name="Default"> 
+
+  </param_set>
+
+</alg_conf>
+
diff --git a/config/NormInteractionListGenerator.xml b/config/NormInteractionListGenerator.xml
@@ -0,0 +1,21 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+
+<alg_conf>
+
+  <!--
+      Configuration for the COHInteractionListGenerator InteractionListGeneratorI
+
+Configurable Parameters:
+.......................................................................................................
+Name             Type     Optional   Comment               Default
+.......................................................................................................
+  -->
+
+  <param_set name="Default"> 
+
+  </param_set>
+
+
+
+</alg_conf>
+
diff --git a/config/NormXSec.xml b/config/NormXSec.xml
@@ -0,0 +1,20 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+
+<!--
+Configuration for the NormCOHXSec EventRecordVisitorI
+
+Algorithm Configurable Parameters:
+.......................................................................................................................
+Name                     Type    Opt   Comment                                        Default
+NormScale                double  Yes   Scale factor                                   1.0
+.......................................................................................................................
+
+-->
+
+<alg_conf>
+
+  <param_set name="Default"> 
+  </param_set>
+
+</alg_conf>
+
diff --git a/config/QELEventGeneratorSM.xml b/config/QELEventGeneratorSM.xml
@@ -4,28 +4,28 @@
 Configuration for the QELEventGeneratorSM EventRecordVisitorI
 
 Configurable Parameters:
-.............................................................................................................
-Name                      Type    Optional   Comment                                                    Default
-.............................................................................................................
-MaxXSec-SafetyFactor      double   Yes       Safety factor for the maximum                              1.2
-                                             differential cross section
-MaxDiffv-SafetyFactor     double   Yes       Safety factor for the maximum                              1.2
-                                             vmax(Q2)-vmin(Q2)
-MaxXSec-DiffTolerance     double   Yes       Max allowed 200*(xsec-xsecmax)/(xsec+xsecmax)              999999.0
-                                             if xsec>xsecmax
-UniformOverPhaseSpace     bool     Yes       Kinematics uniformly over allowed phase space              false
-                                             wgt = (phase_space_volume)*(diff_xsec)/(xsec)
-IsNucleonInNucleus        bool     Yes       Generate struck nucleon in nucleus                         true
-Threshold-Q2              double   Yes       Q2-threshold for seeking the second maximum                2.00
-Cache-MinEnergy           double   Yes       min E for which maxxsec is cached -                        1.00
-                                             forcing explicit calc
+................................................................................................................................
+Name                     Type        Optional   Comment                                                                  Default
+................................................................................................................................
+MaxXSec-SafetyFactor     double      Yes        Safety factor: ComputeMaxXSec -> ComputeMaxXSec * MaxXSec-SafetyFactor   1.0
+Safety-Factors           vec-double  Yes        Safety factors: MaxXSec[nkey] -> MaxXSec[nkey] * Safety-Factors[nkey]    all 1.0
+Interpolator-Types       vec-string  Yes        Type of interpolator for each key in a branch                            all ""
+MaxXSec-DiffTolerance    double      Yes        Max allowed 200*(xsec-xsecmax)/(xsec+xsecmax)                            999999.0
+                                                if xsec>xsecmax                                                          
+UniformOverPhaseSpace    bool        Yes        Kinematics uniformly over allowed phase space                            false
+                                                wgt = (phase_space_volume)*(diff_xsec)/(xsec)                            
+IsNucleonInNucleus       bool        Yes        Generate struck nucleon in nucleus                                       true
+Threshold-Q2             double      Yes        Q2-threshold for seeking the second maximum                              2.00
+Cache-MinEnergy          double      Yes        min E for which maxxsec is cached -                                      1.00
+                                                forcing explicit calc
 -->
 
 
 <alg_conf>
 
   <param_set name="Default"> 
-
+    <param type="vec-double" name="Safety-Factors"        delim=";">  1.06; 1.06; 1.06  </param>
+    <param type="vec-string" name="Interpolator-Types"    delim=";">  AKIMA; AKIMA; LINEAR  </param>
   </param_set>
 
 

diff --git a/config/master_config.xml b/config/master_config.xml
@@ -264,4 +264,8 @@
    <config alg="genie::evtlib::EventLibraryInterface">          EventLibraryInterface.xml          </config>
    <config alg="genie::evtlib::EvtLibInteractionListGenerator"> EvtLibInteractionListGenerator.xml </config>
    <config alg="genie::evtlib::EvtLibPXSec">                    EvtLibPXSec.xml                    </config>
+
+   <config alg="genie::NormGenerator">                       NormGenerator.xml                  </config>
+   <config alg="genie::NormXSec">                            NormXSec.xml                       </config>
+   <config alg="genie::NormInteractionListGenerator">        NormInteractionListGenerator.xml   </config>
 </genie_config>
diff --git a/src/Apps/gNtpConv.cxx b/src/Apps/gNtpConv.cxx
@@ -324,6 +324,7 @@ void ConvertToGST(void)
   bool   brIsMec       = false;  // Is MEC?
   bool   brIsDfr       = false;  // Is Diffractive?
   bool   brIsImd       = false;  // Is IMD?
+  bool   brIsNrm       = false;  // Is Norm?
   bool   brIsSingleK   = false;  // Is single kaon?  
   bool   brIsImdAnh    = false;  // Is IMD annihilation?
   bool   brIsNuEL      = false;  // Is ve elastic?
@@ -438,6 +439,7 @@ void ConvertToGST(void)
   s_tree->Branch("coh",           &brIsCoh,         "coh/O"	    );
   s_tree->Branch("dfr",           &brIsDfr,         "dfr/O"	    );
   s_tree->Branch("imd",	          &brIsImd,	    "imd/O"	    );
+  s_tree->Branch("norm",          &brIsNrm,         "norm/O"	    );
   s_tree->Branch("imdanh",        &brIsImdAnh,	    "imdanh/O"	    );
   s_tree->Branch("singlek",       &brIsSingleK,     "singlek/O"     );  
   s_tree->Branch("nuel",          &brIsNuEL,        "nuel/O"	    );
@@ -657,9 +659,10 @@ void ConvertToGST(void)
     bool is_mec       = proc_info.IsMEC();
     bool is_amnugamma = proc_info.IsAMNuGamma();
     bool is_hnl       = proc_info.IsHNLDecay();
-
+    bool is_norm      = proc_info.IsNorm();
+
     if (!hitnucl && neutrino) {
-        assert(is_coh || is_imd || is_imdanh || is_nuel | is_amnugamma || is_coh_el || is_hnl);
+        assert(is_coh || is_imd || is_imdanh || is_nuel | is_amnugamma || is_coh_el || is_hnl || is_norm);
     }
 
     // Hit quark - set only for DIS events
@@ -954,6 +957,7 @@ void ConvertToGST(void)
     brIsCoh      = is_coh;  
     brIsDfr      = is_dfr;  
     brIsImd      = is_imd;
+    brIsNrm      = is_norm;
     brIsSingleK  = is_singlek;    
     brIsNuEL     = is_nuel;  
     brIsEM       = is_em;  

diff --git a/src/Framework/Conventions/KinePhaseSpace.h b/src/Framework/Conventions/KinePhaseSpace.h
@@ -24,7 +24,7 @@
 using std::string;
 
 namespace genie {
-
+// Note: please attach new phase space enum element to the end of the list .
 typedef enum EKinePhaseSpace {
   kPSNull = 0,
   kPSfE,
@@ -72,7 +72,8 @@ typedef enum EKinePhaseSpace {
   kPSlog10xlog10Q2fE,
   kPSEDNufE, // Used for Dark Neutrinos, two body final state
   kPSn1n2fE,
-  kPSn1n2n3fE
+  kPSn1n2n3fE,
+  kPSQ2vpfE
 } KinePhaseSpace_t;
 
 class KinePhaseSpace
@@ -124,6 +125,7 @@ class KinePhaseSpace
       case(kPSElOlTpifE)  : return "<{El,Omega_l,Theta_pi}|E>"; break;
       case(kPSTkTlctl)    : return "<{Tk,Tl,cos(theta_l)}|E>";  break;
       case(kPSQ2vfE)      : return "<{Q2,v}|E>"; break;
+      case(kPSQ2vpfE)     : return "<{Q2,v,p}|E>"; break;
       // TODO: update this string when the appropriate kinematic variables are known
       case(kPSQELEvGen)   : return "<QELEvGen>"; break;
       case(kPSDMELEvGen)  : return "<DMELEvGen>"; break;

diff --git a/src/Framework/EventGen/GMCJDriver.cxx b/src/Framework/EventGen/GMCJDriver.cxx
@@ -652,9 +652,12 @@ void GMCJDriver::BootstrapXSecSplineSummation(void)
     // bit above the maximum beam energy (but below the maximum valid energy)
     // to avoid the evaluation of the cubic spline around the viscinity of
     // knots with zero y values (although the GENIE Spline object handles it)
-    double dE  = fEmax/10.;
     double min = rE.min;
-    double max = (fEmax+dE < rE.max) ? fEmax+dE : rE.max;
+    double max = rE.max;
+
+    // Because of edge issue (see GENIE docdb 297) these lines are commented out
+    //    double dE  = fEmax/10.;
+    //    double max = (fEmax+dE < rE.max) ? fEmax+dE : rE.max;
 
     // in the logaritmic binning is important to have a narrow binning to
     // describe better the glashow resonance peak.

diff --git a/src/Framework/GHEP/GHepRecord.cxx b/src/Framework/GHEP/GHepRecord.cxx
@@ -1183,8 +1183,8 @@ void GHepRecord::Print(ostream & stream) const
       case ( kPSQ2fE  ) :
         stream << " dsig(Q2;E)/dQ2 =        " << setfill(' ') << setw(13) << fDiffXSec/units::cm2 << " cm^2/GeV^2 |";
         break;
-      case ( kPSQ2vfE  ) :
-        stream << " dsig(Q2,v;E)/dQ2dv =    " << setfill(' ') << setw(13) << fDiffXSec/units::cm2 << " cm^2/GeV^3 |";
+      case ( kPSQ2vpfE  ) :
+        stream << " dsig(Q2,v,p;E)/dQ2dvdp =" << setfill(' ') << setw(13) << fDiffXSec/units::cm2 << " cm^2/GeV^4 |";
         break;
       case ( kPSWQ2fE ) :
         stream << " d2sig(W,Q2;E)/dWdQ2 =   " << setfill(' ') << setw(13) << fDiffXSec/units::cm2 << " cm^2/GeV^3 |";

diff --git a/src/Framework/Interaction/KPhaseSpace.cxx b/src/Framework/Interaction/KPhaseSpace.cxx
@@ -88,6 +88,8 @@ double KPhaseSpace::Threshold(void) const
 
   if( ! pi.IsKnown() ) return 0;
 
+  if (pi.IsNorm() ) return 0;
+
   if (pi.IsSingleKaon()) {
     int kaon_pdgc = xcls.StrangeHadronPdg();
     double Mi   = tgt.HitNucP4Ptr()->M(); // initial nucleon mass

diff --git a/src/Framework/Interaction/ProcessInfo.cxx b/src/Framework/Interaction/ProcessInfo.cxx
@@ -128,6 +128,11 @@ bool ProcessInfo::IsIMDAnnihilation(void) const
   return (fScatteringType == kScIMDAnnihilation);
 }
 //____________________________________________________________________________
+bool ProcessInfo::IsNorm(void) const
+{
+  return (fScatteringType == kScNorm);
+}
+//____________________________________________________________________________
 bool ProcessInfo::IsDarkMatterElectronElastic(void) const
 {
   return (fScatteringType == kScDarkMatterElectron);

diff --git a/src/Framework/Interaction/ProcessInfo.h b/src/Framework/Interaction/ProcessInfo.h
@@ -71,6 +71,7 @@ class ProcessInfo : public TObject {
   bool IsNuElectronElastic       (void) const;
   bool IsInverseMuDecay          (void) const;
   bool IsIMDAnnihilation         (void) const;
+  bool IsNorm                    (void) const;
   bool IsDarkMatterElectronElastic (void) const;
   bool IsInverseBetaDecay        (void) const;
   bool IsGlashowResonance        (void) const;

diff --git a/src/Framework/Interaction/ScatteringType.h b/src/Framework/Interaction/ScatteringType.h
@@ -54,7 +54,8 @@ typedef enum EScatteringType {
   kScPhotonResonance,
   kScDarkMatterElastic = 101,
   kScDarkMatterDeepInelastic,
-  kScDarkMatterElectron
+  kScDarkMatterElectron,
+  kScNorm
 } ScatteringType_t;
 
 class ScatteringType
@@ -80,6 +81,7 @@ class ScatteringType
       case(kScInverseBetaDecay) :        return "IBD";       break;
       case(kScGlashowResonance) :        return "GLR";       break;
       case(kScIMDAnnihilation) :         return "IMDAnh";    break;
+      case(kScNorm) :                    return "Norm";      break;
       case(kScPhotonCoherent) :          return "PhotonCOH"; break;
       case(kScPhotonResonance) :         return "PhotonRES"; break;
       case(kScDarkMatterElastic) :       return "DMEL";      break;