GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MECScaleVsW.cxx
Go to the documentation of this file.
1 //_________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Original code contributed by J.Tena and M.Roda
7  Substantial code refactorizations by the core GENIE group.
8 */
9 //_________________________________________________________________________
10 
15 
16 using namespace genie;
17 
18 //_________________________________________________________________________
20  XSecScaleI("genie::MECScaleVsW")
21 {
22 
23 }
24 //_________________________________________________________________________
25 MECScaleVsW::MECScaleVsW(string config) :
26  XSecScaleI("genie::MECScaleVsW",config)
27 {
28 
29 }
30 //_________________________________________________________________________
32 {
33 
34 }
35 //_________________________________________________________________________
36 double MECScaleVsW::GetScaling( const Interaction & interaction ) const
37 {
38  double Q0 = interaction.Kine().GetKV(kKVQ0) ;
39  double Q3 = interaction.Kine().GetKV(kKVQ3) ;
40 
41  return GetScaling( Q0, Q3 ) ;
42 }
43 //_________________________________________________________________________
44 double MECScaleVsW::GetScaling( const double Q0, const double Q3 ) const
45 {
46  // Get the vectors that include the kinematic limits of W for a given event:
47  MECScaleVsW::weight_type_map weight_map = GetMapWithLimits( Q0, Q3 ) ;
48 
49  // The Scaling is done using the "experimenter's W", which assumes a single nucleon
50  // See motivation in : https://arxiv.org/pdf/1601.02038.pdf
51  // Calculate event W:
52  static double Mn = ( PDGLibrary::Instance()->Find(kPdgProton)->Mass() + PDGLibrary::Instance()->Find(kPdgNeutron)->Mass() ) * 0.5 ; // Nucleon mass
53  double W = pow(Mn,2) + 2*Mn*Q0 - pow(Q3,2) + pow(Q0,2) ;
54  // Do not scale if W<0. This can happen while we try to get the correct kinematics.
55  // If the kinematics is not correct, W can be negative, and we scale with a nan.
56  // To avoid this we do this check.
57  if ( W < 0 ) return 1. ;
58  W = sqrt( W ) ;
59 
60  // Calculate scaling:
61  MECScaleVsW::weight_type_map::iterator it_min = weight_map.begin() ;
62  MECScaleVsW::weight_type_map::iterator it_max = std::next( weight_map.begin(), weight_map.size() -1 ) ;
63 
64  if ( W < it_min->first || W > it_max->first ) return fDefaultWeight ;
65 
66  while ( std::distance( it_min, it_max ) > 1 ) {
67  unsigned int step = std::distance( weight_map.begin(), it_min ) + std::distance( it_min, it_max ) / 2 ;
68  MECScaleVsW::weight_type_map::iterator it_middle = std::next( weight_map.begin(), step ) ;
69  if ( W < it_middle->first ) it_max = it_middle ;
70  else it_min = it_middle ;
71  }
72 
73  return ScaleFunction( W, *it_min, *it_max ) ;
74 }
75 
76 //_________________________________________________________________________
77 MECScaleVsW::weight_type_map MECScaleVsW::GetMapWithLimits( const double Q0, const double Q3 ) const {
78  // This function is responsible to add the phase space limits in the WValues vector in case they are not included
79  // in the configuration setup.
80 
81  static double Mn = ( PDGLibrary::Instance()->Find(kPdgProton)->Mass() + PDGLibrary::Instance()->Find(kPdgNeutron)->Mass() ) * 0.5 ; // Nucleon mass
82  double W_min = sqrt( pow(Mn,2) + 2*Mn*fW1_Q0Q3_limits.Eval(Q3) - pow(Q3,2) + pow(fW1_Q0Q3_limits.Eval(Q3),2) ) ;
83  double W_max = sqrt( pow(Mn,2) + 2*Mn*Q0 ) ; // Imposing Q2 = 0
84 
85  // Insert phase space limits:
87  w_map.insert( weight_type_pair( W_max, fUpperLimitWeight ) ) ;
88  w_map.insert( weight_type_pair( W_min, fLowLimitWeight ) ) ;
89 
90  return w_map ;
91 }
92 //_________________________________________________________________________
93 
94 double MECScaleVsW::ScaleFunction( const double W, const weight_type_pair min, const weight_type_pair max ) const
95 {
96  // This function is responsible to calculate the scale at a given W
97  // It interpolates the value between scale_min (W_min) and scale_max (W_max) linearly
98  return ( max.second - min.second ) * ( W - min.first ) / ( max.first - min.first ) + min.second ;
99 }
100 
101 //_________________________________________________________________________
102 
104 {
105  bool good_config = true ;
106  // Reset members
107  fWeightsMap.clear();
108 
109  if( GetConfig().Exists("MECScaleVsW-Default-Weight") ) {
110  GetParam( "MECScaleVsW-Default-Weight", fDefaultWeight ) ;
111  } else {
112  good_config = false ;
113  LOG("MECScaleVsW", pERROR) << "Default weight is not specified." ;
114  }
115 
116  std::vector<double> Weights, WValues ;
117  GetParamVect( "MECScaleVsW-Weights", Weights, false ) ;
118  GetParamVect( "MECScaleVsW-WValues", WValues, false ) ;
119 
120  if( Weights.size() != WValues.size() ) {
121  good_config = false ;
122  LOG("MECScaleVsW", pERROR) << "Entries don't match" ;
123  LOG("MECScaleVsW", pERROR) << "Weights size: " << Weights.size() ;
124  LOG("MECScaleVsW", pERROR) << "WValues size: " << WValues.size() ;
125  }
126 
127  // Store weights and WValues in map:
128  for( unsigned int i = 0 ; i<Weights.size() ; ++i ) {
129  fWeightsMap.insert( weight_type_pair( WValues[i], Weights[i] ) ) ;
130  }
131 
132  std::vector<double> limit_Q0, limit_Q3 ;
133  if( GetParamVect( "MECScaleVsW-LowerLimitQ0", limit_Q0 ) == 0 ) {
134  good_config = false ;
135  LOG("MECScaleVsW", pERROR) << "MECScaleVsW-LowerLimitQ0 is empty" ;
136  }
137 
138  if( GetParamVect( "MECScaleVsW-LowerLimitQ3", limit_Q3 ) == 0 ) {
139  good_config = false ;
140  LOG("MECScaleVsW", pERROR) << "MECScaleVsW-LowerLimitQ3 is empty" ;
141  }
142 
143  if( limit_Q0.size() != limit_Q3.size() ) {
144  good_config = false ;
145  LOG("MECScaleVsW", pERROR) << "Entries don't match" ;
146  LOG("MECScaleVsW", pERROR) << "Lower limit for Q0 size: " << limit_Q0.size() ;
147  LOG("MECScaleVsW", pERROR) << "Lower limit for Q3 size: " << limit_Q3.size() ;
148  }
149 
150  GetParamDef("MECScleVsW-LowerLimit-Weight", fLowLimitWeight, fDefaultWeight ) ;
151  GetParamDef("MECScleVsW-UpperLimit-Weight", fUpperLimitWeight, fDefaultWeight ) ;
152 
153  if( ! good_config ) {
154  LOG("MECScaleVsW", pERROR) << "Configuration has failed.";
155  exit(78) ;
156  }
157 
158  fW1_Q0Q3_limits = TSpline3("fW1_Q0Q3_limits",limit_Q3.data(),limit_Q0.data(),limit_Q3.size());
159 
160 }
161 
162 //_________________________________________________________________________
#define pERROR
Definition: Messenger.h:59
virtual void LoadConfig(void) override
This class is responsible to compute a scaling factor for the XSec.
Definition: XSecScaleI.h:25
int GetParamVect(const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
Handle to load vectors of parameters.
virtual double GetScaling(const Interaction &) const override
Definition: MECScaleVsW.cxx:36
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:246
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
Summary information for an interaction.
Definition: Interaction.h:56
std::pair< double, double > weight_type_pair
Definition: MECScaleVsW.h:30
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
weight_type_map GetMapWithLimits(const double Q0, const double Q3) const
Definition: MECScaleVsW.cxx:77
const Kinematics & Kine(void) const
Definition: Interaction.h:71
double GetKV(KineVar_t kv) const
Definition: Kinematics.cxx:323
std::map< double, double > weight_type_map
Definition: MECScaleVsW.h:29
virtual double ScaleFunction(const double W, const weight_type_pair min, const weight_type_pair max) const
Definition: MECScaleVsW.cxx:94
double fLowLimitWeight
Definition: MECScaleVsW.h:61
weight_type_map fWeightsMap
Definition: MECScaleVsW.h:57
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
const int kPdgProton
Definition: PDGCodes.h:81
double fUpperLimitWeight
Definition: MECScaleVsW.h:62
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
virtual ~MECScaleVsW()
Definition: MECScaleVsW.cxx:31
const int kPdgNeutron
Definition: PDGCodes.h:83
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
TSpline3 fW1_Q0Q3_limits
Definition: MECScaleVsW.h:59