genie::rew::GReWeightINukeParams::Fates Class Reference

#include <GReWeightINukeParams.h>

Public Member Functions
	Fates (HadronType_t hadtype=kRwINukeUndefined)
	~Fates ()
double	ScaleFactor (GSyst_t s, const TLorentzVector &p4) const
	see next More...
double	ScaleFactor (GSyst_t s, double KE=-1.) const
	fate fraction scale factor = 1 + twk_dial * fractional_err More...
bool	IsIncluded (GSyst_t s) const
	is included? More...
bool	IsCushionTerm (GSyst_t s) const
	is it a cushion term? More...
bool	IsTweaked (GSyst_t s) const
	is included & tweaked to non-def value? More...
bool	IsTweaked (void) const
	is any param tweaked More...
void	Reset (void)
void	Reconfigure (void)
double	ChisqPenalty (void) const
void	SetTwkDial (GSyst_t s, double val)

Private Member Functions
bool	IsHandled (GSyst_t s) const
void	AddCushionTerms (void)
double	ActualTwkDial (GSyst_t s, double KE=-1.) const
	actual tweaking dial for input systematic at input kinetic energy More...

Private Attributes
HadronType_t	fHadType
map< GSyst_t, double >	fSystValuesUser
	List of systematics included & values set by the user. More...
map< GSyst_t, double >	fSystValuesActual
	List of systematics included & values actually used (user values limited to physical range) More...
map< GSyst_t, bool >	fIsCushion
	cushion term flag More...

Detailed Description

Definition at line 72 of file GReWeightINukeParams.h.

Constructor & Destructor Documentation

GReWeightINukeParams::Fates::Fates

(

GReWeightINukeParams::HadronType_t

ht = kRwINukeUndefined

)

Definition at line 146 of file GReWeightINukeParams.cxx.

{
  if (ht != kRwINukePion && ht != kRwINukeNucl) {
    LOG("ReW", pFATAL) 
      << "Mean-free path reweighting code handles only pions or nucleons";
    exit(1);
  }

  fHadType = ht;

  this->Reset();
}

GReWeightINukeParams::Fates::~Fates

(

void

)

Definition at line 159 of file GReWeightINukeParams.cxx.

{
  this->Reset();
}

Member Function Documentation

double GReWeightINukeParams::Fates::ActualTwkDial ( GSyst_t  s, double  KE = -1.  ) const

private

actual tweaking dial for input systematic at input kinetic energy

Definition at line 222 of file GReWeightINukeParams.cxx.

{
  if(! this->IsIncluded(syst)) {
    LOG("ReW", pWARN) 
      << "Systematic " << GSyst::AsString(syst) << " not included";
    return 0.;
  }
  if(KE < 0.) {
    LOG("ReW", pERROR) << "Negative kinetic energy! (KE = " << KE << ")";
    return 0.;
  }

  GSystUncertainty * uncert = GSystUncertainty::Instance();
  map<GSyst_t, double>::const_iterator iter;  

  //
  // initialize actual systematic parameter values to the ones set by the user
  //
  fSystValuesActual.clear();
  fSystValuesActual.insert(fSystValuesUser.begin(), fSystValuesUser.end());

  //
  // Check the tweaking dials set by the user.
  // If a systematic variation is too extreme so that the scale factor for 
  // the corresponding fate fraction becomes negative, automatically re-adjust
  // the systematic variation so that the scale factor is exactly 0
  //

  for(iter = fSystValuesActual.begin(); iter != fSystValuesActual.end(); ++iter)
  {
     GSyst_t curr_syst     = iter->first;
     double  curr_twk_dial = iter->second;

     bool curr_is_cushion = this->IsCushionTerm(curr_syst);
     if(curr_is_cushion) continue;

     double fractional_frac_err = uncert->OneSigmaErr(curr_syst); // fractional 1 sigma error
     double frac_scale = 1. + curr_twk_dial * fractional_frac_err;

     if(frac_scale < 0) {
       double curr_twk_dial_min = -1/fractional_frac_err;
       fSystValuesActual[curr_syst] = curr_twk_dial_min;
       LOG("ReW", pNOTICE) 
         << "Too large systematic variation for non-cushion systematic " << GSyst::AsString(curr_syst) 
         << " makes the corresponding fate scaling factor negative! Taking corrective action...";
       LOG("ReW", pINFO) 
         << "Tweak dial automatically re-adjusted from " << curr_twk_dial
         << " to " << curr_twk_dial_min;
     }
  }

  //
  // Loop over all non-cushion terms and calculate the change in their total fraction
  //

  double sum_nocushion_fate_fraction_nom = 0.;
  double sum_nocushion_fate_fraction_twk = 0.;

  for(iter = fSystValuesActual.begin(); iter != fSystValuesActual.end(); ++iter)
  {
     GSyst_t curr_syst     = iter->first;
     double  curr_twk_dial = iter->second;

     bool curr_is_cushion = this->IsCushionTerm(curr_syst);
     if(curr_is_cushion) continue;

     double fractional_frac_err = uncert->OneSigmaErr(curr_syst); // fractional 1 sigma error

     double frac_scale = 1. + curr_twk_dial * fractional_frac_err;
     double curr_frac  = genie::utils::rew::FateFraction(curr_syst, KE, frac_scale);
     double nom_frac   = genie::utils::rew::FateFraction(curr_syst, KE, 1.);

     sum_nocushion_fate_fraction_nom += nom_frac; 
     sum_nocushion_fate_fraction_twk += curr_frac;

     LOG("ReW", pDEBUG) 
         << "Non-cushion systematic " << GSyst::AsString(curr_syst) 
         << " (1 sigma frac. err = " << 100*fractional_frac_err
         << "%, tweak dial = " << curr_twk_dial 
         << "): fate fraction (KE = " << KE << " GeV) = "
         << nom_frac  << " (nominal), "
         << curr_frac << " (tweaked)";
  }

  LOG("ReW", pDEBUG) 
    << "Sum of non-cushion fate fractions = "
    << sum_nocushion_fate_fraction_nom << " (nominal), "
    << sum_nocushion_fate_fraction_twk << " (tweaked)";

  //
  // If the systematic variation specified by the user was so extreme that the
  // tweaked total fate fraction of non-cushion terms is above 1, automatically re-scale the
  // tweak dials so that the sum of their fate fraction is 1.
  // In that case, all cushion terms will be set to 0.
  //
  
  if(sum_nocushion_fate_fraction_twk > 1.) {
 
    LOG("ReW", pNOTICE) 
       << "Current sum of non-cushion fate fractions = "
       << sum_nocushion_fate_fraction_twk << " ( > 1 ). "
       << "Can not maintain unitarity with current set of systematic parameter values. "
       << "Taking corrective action...";

    for(iter = fSystValuesActual.begin(); iter != fSystValuesActual.end(); ++iter)
    {
       GSyst_t curr_syst     = iter->first;
       double  curr_twk_dial = iter->second;

       bool curr_is_cushion = this->IsCushionTerm(curr_syst);

       double fractional_frac_err = uncert->OneSigmaErr(curr_syst); // fractional 1 sigma error
  
       double frac_scale = 1. + curr_twk_dial * fractional_frac_err;
       double curr_frac  = genie::utils::rew::FateFraction(curr_syst, KE, frac_scale);

       double curr_frac_new     = (curr_is_cushion) ? 0 : curr_frac * (1./sum_nocushion_fate_fraction_twk);
       double frac_scale_new    = genie::utils::rew::WhichFateFractionScaleFactor(curr_syst, KE, curr_frac_new);
       double curr_twk_dial_new = (frac_scale_new != 0) ? (frac_scale_new - 1.)/fractional_frac_err : 0;

       fSystValuesActual[curr_syst] = curr_twk_dial_new;

       LOG("ReW", pINFO) 
         << ((curr_is_cushion) ? " Cushion " : " Non-cushion ")
         << "systematic " << GSyst::AsString(curr_syst) 
         << ": fate fraction re-adjusted from " << curr_frac << " to " << curr_frac_new
         << " (tweak dial re-adjusted from " << curr_twk_dial
         << " to " << curr_twk_dial_new << ")";
    }

    // update sum of fate fraction (set at limit)
    sum_nocushion_fate_fraction_twk = 1.;
  }

  // 
  // In the normal case where the sum of all non-cushion terms was less than 1,
  // leave them as they are. Adjust all all non-cushion terms accordingly so as to respect unitarity.
  // There are many ways to adjust the cushion terms, if more than one such term exists.
  // Chose to tweak them so that they all change the same amount, in units of the corresponding 1 sigma error.
  //

  else {

      double delta_fate_fraction = 
          sum_nocushion_fate_fraction_twk - sum_nocushion_fate_fraction_nom;

      LOG("ReW", pDEBUG) 
        << "Sum of non-cushion fate fractions changed by " << delta_fate_fraction
        << " - Change to be absorbed by the selected cushion terms...";

      double sum = 0;
      for(iter = fSystValuesActual.begin(); iter != fSystValuesActual.end(); ++iter)
      {
         GSyst_t curr_syst = iter->first;

         bool curr_is_cushion = this->IsCushionTerm(curr_syst);
         if(!curr_is_cushion) continue;

         double fractional_frac_err = uncert->OneSigmaErr(curr_syst); // fractional 1 sigma error
         double nom_frac = genie::utils::rew::FateFraction(curr_syst, KE, 1.);
         sum += (nom_frac * fractional_frac_err);
      }

      double twk_dial = -1.*delta_fate_fraction / sum;

      for(iter = fSystValuesActual.begin(); iter != fSystValuesActual.end(); ++iter)
      {
         GSyst_t curr_syst  = iter->first;
         bool curr_is_cushion = this->IsCushionTerm(curr_syst);
         if(!curr_is_cushion) continue;

         fSystValuesActual[curr_syst] = twk_dial;
      }

      LOG("ReW", pDEBUG) 
        << "To absorb change in the sum of non-cushion fate fractions, "
        << "the tweaking dial for all cushion terms needs to be set to "
        << twk_dial;

  } //check sum_nocushion_fate_fraction_twk < 1.

  //
  // Confirm unitarity
  //

  LOG("ReW", pINFO) << "Confirming unitarity for current fate fractions";

  double sum_fate_fraction_all = 0;
  for(iter = fSystValuesActual.begin(); iter != fSystValuesActual.end(); ++iter)
  {
     GSyst_t curr_syst     = iter->first;
     double  curr_twk_dial = iter->second;

     bool curr_is_cushion = this->IsCushionTerm(curr_syst);

     double fractional_frac_err = uncert->OneSigmaErr(curr_syst); // fractional 1 sigma error

     double frac_scale = 1. + curr_twk_dial * fractional_frac_err;
     double curr_frac  = genie::utils::rew::FateFraction(curr_syst, KE, frac_scale);

     LOG("ReW", pINFO) 
         << "Systematic " << GSyst::AsString(curr_syst) 
         << " (1 sigma frac. err = " << 100*fractional_frac_err
         << "%, tweak dial = " << curr_twk_dial 
         << "): Current fate fraction (KE = " << KE << " GeV) = " << curr_frac
         << ((curr_is_cushion) ? " ** cushion term ** " : "");

     sum_fate_fraction_all += curr_frac;
  }

  LOG("ReW", pINFO) << "Current sum of all fate fractions = " << sum_fate_fraction_all;

  if(TMath::Abs(sum_fate_fraction_all-1) > 0.01) {   
      LOG("ReW", pWARN) << "Unitarity violation level exceeded 1 part in 100.";
      LOG("ReW", pWARN) << "Current sum of all fate fractions = " << sum_fate_fraction_all;
  }

  map<GSyst_t, double>::const_iterator dial_iter = fSystValuesActual.find(syst);
  return dial_iter->second;
}

void GReWeightINukeParams::Fates::AddCushionTerms ( void  )

private

Definition at line 545 of file GReWeightINukeParams.cxx.

{
// When this method is called for the first time, all terms not already 
// included are included as cushion terms.
// On every call, the values of the cushion terms are reset.

  LOG("ReW", pINFO) << "Adding cushion terms...";

  int ncushions = 0;

  GSyst_t syst = kNullSystematic;
  int i=0;
  while( (syst = (fHadType == kRwINukePion) ?
            GSyst::NextPionFateSystematic(i++) :
            GSyst::NextNuclFateSystematic(i++)
         ) != kNullSystematic
       )
  {
      if(this->IsIncluded(syst)) {
        bool is_cushion = this->IsCushionTerm(syst);
        LOG("ReW", pINFO) 
          << "Systematic " << GSyst::AsString(syst) 
          << " was already specified as a"
          << ((is_cushion) ? " cuhsion " : " non-cushion ")
          << "term";
        if(is_cushion) {
           fSystValuesUser[syst]  = 0.; 
           ncushions++;
        }
      }
      else {
        LOG("ReW", pINFO) 
          << "Adding systematic " << GSyst::AsString(syst) 
          << " as a cushion term.";

        fSystValuesUser[syst]  = 0.; 
        fIsCushion[syst] = true;
        ncushions++;
      }
  } // end loop over relevant fate systematics

  LOG("ReW", pINFO) << "Number of cushion terms = " << ncushions;

  if(ncushions <= 0) {
     LOG("ReW", pFATAL) 
       << "There must be at least one cushion term (" 
       << ncushions << " were set)";
     exit(1);
  }
}

double GReWeightINukeParams::Fates::ChisqPenalty

(

void

)

const

Definition at line 443 of file GReWeightINukeParams.cxx.

{
// calculate an average chisq over the energy range

  double chisq = 0.0;

  const int    kN      = 200;
  const double kKEmin  =  0;
  const double kKEmax  = 10;
  const double kKEstep = (kKEmax-kKEmin)/(kN-1);

  GSyst_t syst = kNullSystematic;
  int i=0;
  while( (syst = (fHadType == kRwINukePion) ?
            GSyst::NextPionFateSystematic(i++) :
            GSyst::NextNuclFateSystematic(i++)
         ) != kNullSystematic
       )
  {
    for(int j=0; j < kN; j++) {
       double KE = kKEmin + j * kKEstep;
       double twkdial = this->ActualTwkDial(syst,KE);
       chisq += (TMath::Power(twkdial, 2.));
    }//j
  }//i 

  chisq /= kN;

  return chisq;
}

bool GReWeightINukeParams::Fates::IsCushionTerm

(

GSyst_t

s

)

const

is it a cushion term?

Definition at line 481 of file GReWeightINukeParams.cxx.

{
  map<GSyst_t, bool>::const_iterator iter = fIsCushion.find(syst);
  if(iter != fIsCushion.end()) return iter->second;

  LOG("ReW", pERROR) 
    << "Cannot query 'is cushion' flag for systematic " 
    << GSyst::AsString(syst);

  return false;
}

bool GReWeightINukeParams::Fates::IsHandled ( GSyst_t  s ) const

private

Definition at line 527 of file GReWeightINukeParams.cxx.

{
//
//
  if(fHadType == kRwINukePion) {
    if(GSyst::IsINukePionFateSystematic(syst)) return true;
  }
  else 
  if(fHadType == kRwINukeNucl) {
    if(GSyst::IsINukeNuclFateSystematic(syst)) return true;
  }

  LOG("ReW", pWARN) 
    << "Can not handle systematic " << GSyst::AsString(syst);

  return false;
}

bool GReWeightINukeParams::Fates::IsIncluded

(

GSyst_t

s

)

const

is included?

Definition at line 474 of file GReWeightINukeParams.cxx.

{
  map<GSyst_t, double>::const_iterator iter = fSystValuesUser.find(syst);
  if(iter != fSystValuesUser.end()) return true;
  return false;
}

bool GReWeightINukeParams::Fates::IsTweaked

(

GSyst_t

s

)

const

is included & tweaked to non-def value?

Definition at line 493 of file GReWeightINukeParams.cxx.

{
  if(!this->IsHandled (syst)) return false;
  if(!this->IsIncluded(syst)) return false;

  map<GSyst_t, double>::const_iterator dial_iter = fSystValuesUser.find(syst);
  double dial = dial_iter->second;

  bool tweaked = (TMath::Abs(dial) >= controls::kASmallNum);

  LOG("ReW", pDEBUG) 
    << "Systematic " << GSyst::AsString(syst)
    << ((tweaked) ? " was " : " was not ")
    << "tweaked (tweaking dial = " << dial << ")";

  return tweaked;
}

bool GReWeightINukeParams::Fates::IsTweaked

(

void

)

const

is any param tweaked

Definition at line 511 of file GReWeightINukeParams.cxx.

{
  GSyst_t syst = kNullSystematic;
  int i=0;
  while( (syst = (fHadType == kRwINukePion) ?
            GSyst::NextPionFateSystematic(i++) :
            GSyst::NextNuclFateSystematic(i++)
         ) != kNullSystematic
       )
  {
    bool tweaked = this->IsTweaked(syst);
    if(tweaked) return true;
  }
  return false;
}

void GReWeightINukeParams::Fates::Reconfigure

(

void

)

Definition at line 171 of file GReWeightINukeParams.cxx.

{
  this->AddCushionTerms();
}

void GReWeightINukeParams::Fates::Reset

(

void

)

Definition at line 164 of file GReWeightINukeParams.cxx.

{
  fSystValuesUser.clear();
  fSystValuesActual.clear();
  fIsCushion.clear();
}

double GReWeightINukeParams::Fates::ScaleFactor	(	GSyst_t	s,
		const TLorentzVector &	p4
	)		const

see next

Definition at line 196 of file GReWeightINukeParams.cxx.

{
  double KE = p4.Energy() - p4.M(); // kinetic energy
  return this->ScaleFactor(syst, KE);
}

double GReWeightINukeParams::Fates::ScaleFactor	(	GSyst_t	s,
		double	KE = -1.
	)		const

fate fraction scale factor = 1 + twk_dial * fractional_err

Definition at line 203 of file GReWeightINukeParams.cxx.

{
  GSystUncertainty * uncert = GSystUncertainty::Instance();

  double fractional_error = uncert->OneSigmaErr(syst);
  double twk_dial         = this->ActualTwkDial(syst, KE);

  double fate_fraction_scale = 1. + twk_dial * fractional_error;

  LOG("ReW", pNOTICE) 
      << "Systematic " << GSyst::AsString(syst) 
      << " (1 sigma frac. err = " << 100*fractional_error
      << "%, tweak dial = " << twk_dial
      << "): Weight = " << fate_fraction_scale;

  return fate_fraction_scale;
}

void GReWeightINukeParams::Fates::SetTwkDial	(	GSyst_t	s,
		double	val
	)

Definition at line 176 of file GReWeightINukeParams.cxx.

{
  // check type of systematic 
  if(!this->IsHandled(syst)) return;

  // can not explicitly set params designated as cushion terms
  if(this->IsIncluded(syst)) {
    if(this->IsCushionTerm(syst)) {
       LOG("ReW", pWARN) 
         << "You may not set the value of cushion term " << GSyst::AsString(syst)
         << ". It is set automatically to maintain unitarity";
       return;
    }
  }

  // update tweaking dial
  fSystValuesUser[syst] = val;
  fIsCushion[syst]   = false;
}

Member Data Documentation

HadronType_t genie::rew::GReWeightINukeParams::Fates::fHadType

private

Definition at line 94 of file GReWeightINukeParams.h.

map<GSyst_t, bool> genie::rew::GReWeightINukeParams::Fates::fIsCushion

private

cushion term flag

Definition at line 98 of file GReWeightINukeParams.h.

map<GSyst_t, double> genie::rew::GReWeightINukeParams::Fates::fSystValuesActual

mutableprivate

List of systematics included & values actually used (user values limited to physical range)

Definition at line 97 of file GReWeightINukeParams.h.

map<GSyst_t, double> genie::rew::GReWeightINukeParams::Fates::fSystValuesUser

private

List of systematics included & values set by the user.

Definition at line 95 of file GReWeightINukeParams.h.

---

The documentation for this class was generated from the following files:

• GENIE/src/ReWeight/GReWeightINukeParams.h
• GENIE/src/ReWeight/GReWeightINukeParams.cxx