This class is responsible to compute the MEC scaling factor given Q0, Q3. The scaling is done as a function of the hadronic invariant mass. More...

#include <MECScaleVsW.h>

Inheritance diagram for genie::MECScaleVsW:

[legend]

Collaboration diagram for genie::MECScaleVsW:

[legend]

Public Member Functions
	MECScaleVsW ()

	MECScaleVsW (string config)

virtual	~MECScaleVsW ()

virtual double	GetScaling (const Interaction &) const override

Public Member Functions inherited from genie::XSecScaleI
virtual	~XSecScaleI ()

Public Member Functions inherited from genie::Algorithm
virtual	~Algorithm ()

virtual void	FindConfig (void)

virtual const Registry &	GetConfig (void) const

Registry *	GetOwnedConfig (void)

virtual const AlgId &	Id (void) const
	Get algorithm ID. More...

virtual AlgStatus_t	GetStatus (void) const
	Get algorithm status. More...

virtual bool	AllowReconfig (void) const

virtual AlgCmp_t	Compare (const Algorithm *alg) const
	Compare with input algorithm. More...

virtual void	SetId (const AlgId &id)
	Set algorithm ID. More...

virtual void	SetId (string name, string config)

const Algorithm *	SubAlg (const RgKey &registry_key) const

void	AdoptConfig (void)

void	AdoptSubstructure (void)

virtual void	Print (ostream &stream) const
	Print algorithm info. More...

Protected Member Functions
virtual void	LoadConfig (void) override

virtual double	GetScaling (const double Q0, const double Q3) const

weight_type_map	GetMapWithLimits (const double Q0, const double Q3) const

virtual double	ScaleFunction (const double W, const weight_type_pair min, const weight_type_pair max) const

Protected Member Functions inherited from genie::XSecScaleI
	XSecScaleI (string name, string config="Default")

void	Configure (const Registry &config) override

virtual void	Configure (string config) override

Protected Member Functions inherited from genie::Algorithm
	Algorithm ()

	Algorithm (string name)

	Algorithm (string name, string config)

void	Initialize (void)

void	DeleteConfig (void)

void	DeleteSubstructure (void)

Registry *	ExtractLocalConfig (const Registry &in) const

Registry *	ExtractLowerConfig (const Registry &in, const string &alg_key) const
	Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key. More...

template<class T >
bool	GetParam (const RgKey &name, T &p, bool is_top_call=true) const

template<class T >
bool	GetParamDef (const RgKey &name, T &p, const T &def) const

template<class T >
int	GetParamVect (const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
	Handle to load vectors of parameters. More...

int	GetParamVectKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const

template<class T >
int	GetParamMat (const std::string &comm_name, TMatrixT< T > &mat, bool is_top_call=true) const
	Handle to load matrix of parameters. More...

template<class T >
int	GetParamMatSym (const std::string &comm_name, TMatrixTSym< T > &mat, bool is_top_call=true) const

int	GetParamMatKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const

int	AddTopRegistry (Registry *rp, bool owns=true)
	add registry with top priority, also update ownership More...

int	AddLowRegistry (Registry *rp, bool owns=true)
	add registry with lowest priority, also update ownership More...

int	MergeTopRegistry (const Registry &r)

int	AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
	Add registries with top priority, also udated Ownerships. More...

Private Types
using	weight_type_map = std::map< double, double >

using	weight_type_pair = std::pair< double, double >

Private Attributes
double	fDefaultWeight

weight_type_map	fWeightsMap

TSpline3	fW1_Q0Q3_limits

double	fLowLimitWeight

double	fUpperLimitWeight

Additional Inherited Members
Static Public Member Functions inherited from genie::Algorithm
static string	BuildParamVectKey (const std::string &comm_name, unsigned int i)

static string	BuildParamVectSizeKey (const std::string &comm_name)

static string	BuildParamMatKey (const std::string &comm_name, unsigned int i, unsigned int j)

static string	BuildParamMatRowSizeKey (const std::string &comm_name)

static string	BuildParamMatColSizeKey (const std::string &comm_name)

Protected Attributes inherited from genie::Algorithm
bool	fAllowReconfig

bool	fOwnsSubstruc
	true if it owns its substructure (sub-algs,...) More...

AlgId	fID
	algorithm name and configuration set More...

vector< Registry * >	fConfVect

vector< bool >	fOwnerships
	ownership for every registry in fConfVect More...

AlgStatus_t	fStatus
	algorithm execution status More...

AlgMap *	fOwnedSubAlgMp
	local pool for owned sub-algs (taken out of the factory pool) More...

Detailed Description

This class is responsible to compute the MEC scaling factor given Q0, Q3. The scaling is done as a function of the hadronic invariant mass.

Author: Code contributed by J.Tena Vidal and M.Roda

Created:: June, 2020

License:: Copyright (c) 2003-2024, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Definition at line 27 of file MECScaleVsW.h.

Member Typedef Documentation

using genie::MECScaleVsW::weight_type_map = std::map<double,double>

private

Definition at line 29 of file MECScaleVsW.h.

using genie::MECScaleVsW::weight_type_pair = std::pair<double,double>

private

Definition at line 30 of file MECScaleVsW.h.

Constructor & Destructor Documentation

MECScaleVsW::MECScaleVsW ( )

Definition at line 19 of file MECScaleVsW.cxx.

                          : 
   XSecScaleI("genie::MECScaleVsW")
 {
   
 }

MECScaleVsW::MECScaleVsW ( string config )

Definition at line 25 of file MECScaleVsW.cxx.

                                       : 
   XSecScaleI("genie::MECScaleVsW",config) 
 {
   
 }

MECScaleVsW::~MECScaleVsW ( )

virtual

Definition at line 31 of file MECScaleVsW.cxx.

 {
 
 }

Member Function Documentation

MECScaleVsW::weight_type_map MECScaleVsW::GetMapWithLimits	(	const double	Q0,
		const double	Q3
	)		const

protected

Definition at line 77 of file MECScaleVsW.cxx.

References genie::PDGLibrary::Find(), fLowLimitWeight, fUpperLimitWeight, fW1_Q0Q3_limits, fWeightsMap, genie::PDGLibrary::Instance(), genie::kPdgNeutron, and genie::kPdgProton.

Referenced by GetScaling().

                                                                                                  {
   // This function is responsible to add the phase space limits in the WValues vector in case they are not included
   // in the configuration setup.
 
   static double Mn = ( PDGLibrary::Instance()->Find(kPdgProton)->Mass() + PDGLibrary::Instance()->Find(kPdgNeutron)->Mass() ) * 0.5 ;  // Nucleon mass
   double W_min = sqrt( pow(Mn,2) + 2*Mn*fW1_Q0Q3_limits.Eval(Q3) - pow(Q3,2) + pow(fW1_Q0Q3_limits.Eval(Q3),2) ) ;
   double W_max = sqrt( pow(Mn,2) + 2*Mn*Q0 ) ; // Imposing Q2 = 0 
 
   // Insert phase space limits:
   MECScaleVsW::weight_type_map w_map = fWeightsMap ; 
   w_map.insert( weight_type_pair( W_max, fUpperLimitWeight ) ) ;
   w_map.insert( weight_type_pair( W_min, fLowLimitWeight ) ) ;
 
   return w_map ; 
 } 

double MECScaleVsW::GetScaling ( const Interaction & interaction ) const

overridevirtual

Implements genie::XSecScaleI.

Definition at line 36 of file MECScaleVsW.cxx.

References genie::Kinematics::GetKV(), genie::Interaction::Kine(), genie::kKVQ0, and genie::kKVQ3.

 {
   double Q0 = interaction.Kine().GetKV(kKVQ0) ; 
   double Q3 = interaction.Kine().GetKV(kKVQ3) ; 
   
   return GetScaling( Q0, Q3 ) ; 
 }

double MECScaleVsW::GetScaling	(	const double	Q0,
		const double	Q3
	)		const

protectedvirtual

Definition at line 44 of file MECScaleVsW.cxx.

References fDefaultWeight, genie::PDGLibrary::Find(), GetMapWithLimits(), genie::PDGLibrary::Instance(), genie::kPdgNeutron, genie::kPdgProton, ScaleFunction(), and genie::utils::kinematics::W().

 {
   // Get the vectors that include the kinematic limits of W for a given event:
   MECScaleVsW::weight_type_map weight_map = GetMapWithLimits( Q0, Q3 ) ;
     
   // The Scaling is done using the "experimenter's W", which assumes a single nucleon
   // See motivation in : https://arxiv.org/pdf/1601.02038.pdf
   // Calculate event W:
   static double Mn = ( PDGLibrary::Instance()->Find(kPdgProton)->Mass() + PDGLibrary::Instance()->Find(kPdgNeutron)->Mass() ) * 0.5 ;  // Nucleon mass
   double W = pow(Mn,2) + 2*Mn*Q0 - pow(Q3,2) + pow(Q0,2) ;
   // Do not scale if W<0. This can happen while we try to get the correct kinematics.
   // If the kinematics is not correct, W can be negative, and we scale with a nan.
   // To avoid this we do this check.
   if ( W < 0 ) return 1. ; 
   W = sqrt( W ) ; 
 
   // Calculate scaling:
   MECScaleVsW::weight_type_map::iterator it_min = weight_map.begin() ; 
   MECScaleVsW::weight_type_map::iterator it_max = std::next( weight_map.begin(), weight_map.size() -1 ) ; 
 
   if ( W < it_min->first || W > it_max->first ) return fDefaultWeight ; 
 
   while ( std::distance( it_min, it_max ) > 1 ) {
     unsigned int step = std::distance( weight_map.begin(), it_min ) + std::distance( it_min, it_max ) / 2 ; 
     MECScaleVsW::weight_type_map::iterator it_middle = std::next( weight_map.begin(), step ) ; 
     if ( W < it_middle->first ) it_max = it_middle ; 
     else it_min = it_middle ; 
   }
 
   return ScaleFunction( W, *it_min, *it_max ) ; 
 } 

void MECScaleVsW::LoadConfig ( void )

overrideprotectedvirtual

Implements genie::XSecScaleI.

Definition at line 103 of file MECScaleVsW.cxx.

References fDefaultWeight, fLowLimitWeight, fUpperLimitWeight, fW1_Q0Q3_limits, fWeightsMap, genie::Algorithm::GetConfig(), genie::Algorithm::GetParam(), genie::Algorithm::GetParamDef(), genie::Algorithm::GetParamVect(), LOG, and pERROR.

 {
   bool good_config = true ; 
   // Reset members
   fWeightsMap.clear(); 
 
   if( GetConfig().Exists("MECScaleVsW-Default-Weight") ) {
     GetParam( "MECScaleVsW-Default-Weight", fDefaultWeight ) ;
   } else {
       good_config = false ; 
       LOG("MECScaleVsW", pERROR) << "Default weight is not specified." ;
   }
 
   std::vector<double> Weights, WValues ;
   GetParamVect( "MECScaleVsW-Weights", Weights, false ) ; 
   GetParamVect( "MECScaleVsW-WValues", WValues, false ) ; 
   
   if( Weights.size() != WValues.size() ) {
     good_config = false ; 
     LOG("MECScaleVsW", pERROR) << "Entries don't match" ;
     LOG("MECScaleVsW", pERROR) << "Weights size: " << Weights.size() ;
     LOG("MECScaleVsW", pERROR) << "WValues size: " << WValues.size() ;
   }
 
   // Store weights and WValues in map:
   for( unsigned int i = 0 ; i<Weights.size() ; ++i ) {
     fWeightsMap.insert( weight_type_pair( WValues[i], Weights[i] ) ) ;
   }
   
   std::vector<double> limit_Q0, limit_Q3 ;
   if( GetParamVect( "MECScaleVsW-LowerLimitQ0", limit_Q0 ) == 0 ) {
     good_config = false ; 
     LOG("MECScaleVsW", pERROR) << "MECScaleVsW-LowerLimitQ0 is empty" ;
   }
 
   if( GetParamVect( "MECScaleVsW-LowerLimitQ3", limit_Q3 ) == 0 ) {
     good_config = false ; 
     LOG("MECScaleVsW", pERROR) << "MECScaleVsW-LowerLimitQ3 is empty" ;
   }
   
   if( limit_Q0.size() != limit_Q3.size() ) {
     good_config = false ; 
     LOG("MECScaleVsW", pERROR) << "Entries don't match" ;
     LOG("MECScaleVsW", pERROR) << "Lower limit for Q0 size: " << limit_Q0.size() ;
     LOG("MECScaleVsW", pERROR) << "Lower limit for Q3 size: " << limit_Q3.size() ;
   }
 
   GetParamDef("MECScleVsW-LowerLimit-Weight", fLowLimitWeight, fDefaultWeight ) ; 
   GetParamDef("MECScleVsW-UpperLimit-Weight", fUpperLimitWeight, fDefaultWeight ) ; 
 
   if( ! good_config ) {
     LOG("MECScaleVsW", pERROR) << "Configuration has failed.";
     exit(78) ;
   }
 
   fW1_Q0Q3_limits = TSpline3("fW1_Q0Q3_limits",limit_Q3.data(),limit_Q0.data(),limit_Q3.size()); 
   
 }

double MECScaleVsW::ScaleFunction	(	const double	W,
		const weight_type_pair	min,
		const weight_type_pair	max
	)		const

protectedvirtual

Definition at line 94 of file MECScaleVsW.cxx.

Referenced by GetScaling().

 {
   // This function is responsible to calculate the scale at a given W
   // It interpolates the value between scale_min (W_min) and scale_max (W_max) linearly 
   return ( max.second - min.second ) * ( W - min.first ) / ( max.first - min.first ) + min.second ; 
 } 