Page MenuHomeHEPForge
Files F8310140

No OneTemporary
Actions

View Options

diff --git a/src/Reweight/WeightUtils.cxx b/src/Reweight/WeightUtils.cxx
index a17eb22..1a1bd77 100644
--- a/src/Reweight/WeightUtils.cxx
+++ b/src/Reweight/WeightUtils.cxx
@@ -1,615 +1,614 @@
#include "WeightUtils.h"
#include "FitLogger.h"
#ifdef __T2KREW_ENABLED__
#include "T2KGenieReWeight.h"
#include "T2KNIWGReWeight.h"
#include "T2KNIWGUtils.h"
#include "T2KNeutReWeight.h"
#include "T2KNeutUtils.h"
#include "T2KReWeight.h"
using namespace t2krew;
#endif
#ifdef __NIWG_ENABLED__
#include "NIWGReWeight.h"
#include "NIWGReWeight1piAngle.h"
#include "NIWGReWeight2010a.h"
#include "NIWGReWeight2012a.h"
#include "NIWGReWeight2014a.h"
#include "NIWGReWeightDeltaMass.h"
#include "NIWGReWeightEffectiveRPA.h"
#include "NIWGReWeightHadronMultSwitch.h"
#include "NIWGReWeightMEC.h"
#include "NIWGReWeightPiMult.h"
#include "NIWGReWeightProtonFSIbug.h"
#include "NIWGReWeightRPA.h"
#include "NIWGReWeightSpectralFunc.h"
#include "NIWGReWeightSplineEnu.h"
#include "NIWGSyst.h"
#include "NIWGSystUncertainty.h"
#endif
#ifdef __NEUT_ENABLED__
#include "NReWeight.h"
#include "NReWeightCasc.h"
#include "NReWeightNuXSecCCQE.h"
#include "NReWeightNuXSecCCRES.h"
#include "NReWeightNuXSecCOH.h"
#include "NReWeightNuXSecDIS.h"
#include "NReWeightNuXSecNC.h"
#include "NReWeightNuXSecNCEL.h"
#include "NReWeightNuXSecNCRES.h"
#include "NReWeightNuXSecRES.h"
#include "NReWeightNuclPiless.h"
#include "NSyst.h"
#include "NSystUncertainty.h"
#include "neutpart.h"
#include "neutvect.h"
#endif
#ifdef __NUWRO_ENABLED__
#include "event1.h"
#endif
#ifdef __NUWRO_REWEIGHT_ENABLED__
#include "NuwroReWeight.h"
#include "NuwroReWeight_FlagNorm.h"
#include "NuwroReWeight_QEL.h"
#include "NuwroReWeight_SPP.h"
#include "NuwroSyst.h"
#include "NuwroSystUncertainty.h"
#endif
#ifdef __GENIE_ENABLED__
#include "EVGCore/EventRecord.h"
#include "EVGCore/EventRecord.h"
#include "GHEP/GHepRecord.h"
#include "GSyst.h"
#include "GSystUncertainty.h"
#include "Ntuple/NtpMCEventRecord.h"
#include "ReWeight/GReWeight.h"
#include "ReWeight/GReWeightAGKY.h"
#include "ReWeight/GReWeightDISNuclMod.h"
#include "ReWeight/GReWeightFGM.h"
#include "ReWeight/GReWeightFZone.h"
#include "ReWeight/GReWeightINuke.h"
#include "ReWeight/GReWeightNonResonanceBkg.h"
#include "ReWeight/GReWeightNuXSecCCQE.h"
#include "ReWeight/GReWeightNuXSecCCQEvec.h"
#include "ReWeight/GReWeightNuXSecCCRES.h"
#include "ReWeight/GReWeightNuXSecCOH.h"
#include "ReWeight/GReWeightNuXSecDIS.h"
#include "ReWeight/GReWeightNuXSecNC.h"
#include "ReWeight/GReWeightNuXSecNCEL.h"
#include "ReWeight/GReWeightNuXSecNCRES.h"
#include "ReWeight/GReWeightResonanceDecay.h"
using namespace genie;
using namespace genie::rew;
#endif
#include "GlobalDialList.h"
#include "ModeNormEngine.h"
#include "NUISANCESyst.h"
#include "OscWeightEngine.h"
//********************************************************************
TF1 FitBase::GetRWConvFunction(std::string const &type,
std::string const &name) {
//********************************************************************
std::string dialfunc = "x";
std::string parType = type;
double low = -10000.0;
double high = 10000.0;
if (parType.find("parameter") == std::string::npos) parType += "_parameter";
std::string line;
- ifstream card(
+ std::ifstream card(
(GeneralUtils::GetTopLevelDir() + "/parameters/dial_conversion.card")
.c_str(),
- ifstream::in);
+ std::ifstream::in);
while (std::getline(card >> std::ws, line, '\n')) {
std::vector<std::string> inputlist = GeneralUtils::ParseToStr(line, " ");
// Check the line length
if (inputlist.size() < 4) continue;
// Check whether this is a comment
if (inputlist[0].c_str()[0] == '#') continue;
// Check whether this is the correct parameter type
if (inputlist[0].compare(parType) != 0) continue;
// Check the parameter name
if (inputlist[1].compare(name) != 0) continue;
// inputlist[2] should be the units... ignore for now
dialfunc = inputlist[3];
// High and low are optional, check whether they exist
if (inputlist.size() > 4) low = GeneralUtils::StrToDbl(inputlist[4]);
if (inputlist.size() > 5) high = GeneralUtils::StrToDbl(inputlist[5]);
}
TF1 convfunc = TF1((name + "_convfunc").c_str(), dialfunc.c_str(), low, high);
return convfunc;
}
//********************************************************************
std::string FitBase::GetRWUnits(std::string const &type,
std::string const &name) {
//********************************************************************
std::string unit = "sig.";
std::string parType = type;
if (parType.find("parameter") == std::string::npos) {
parType += "_parameter";
}
std::string line;
std::ifstream card(
(GeneralUtils::GetTopLevelDir() + "/parameters/dial_conversion.card")
.c_str(),
- ifstream::in);
+ std::ifstream::in);
while (std::getline(card >> std::ws, line, '\n')) {
std::vector<std::string> inputlist = GeneralUtils::ParseToStr(line, " ");
// Check the line length
if (inputlist.size() < 3) continue;
// Check whether this is a comment
if (inputlist[0].c_str()[0] == '#') continue;
// Check whether this is the correct parameter type
if (inputlist[0].compare(parType) != 0) continue;
// Check the parameter name
if (inputlist[1].compare(name) != 0) continue;
unit = inputlist[2];
break;
}
return unit;
}
//********************************************************************
double FitBase::RWAbsToSigma(std::string const &type, std::string const &name,
double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.GetX(val);
if (fabs(conv_val) < 1E-10) conv_val = 0.0;
- std::cout << "AbsToSigma(" << name << ") = " << val << " -> " << conv_val
- << std::endl;
+ QLOG(FIT, "AbsToSigma(" << name << ") = " << val << " -> " << conv_val);
return conv_val;
}
//********************************************************************
double FitBase::RWSigmaToAbs(std::string const &type, std::string const &name,
double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.Eval(val);
return conv_val;
}
//********************************************************************
double FitBase::RWFracToSigma(std::string const &type, std::string const &name,
double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.GetX((val * f1.Eval(0.0)));
if (fabs(conv_val) < 1E-10) conv_val = 0.0;
return conv_val;
}
//********************************************************************
double FitBase::RWSigmaToFrac(std::string const &type, std::string const &name,
double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.Eval(val) / f1.Eval(0.0);
return conv_val;
}
int FitBase::ConvDialType(std::string const &type) {
if (!type.compare("neut_parameter"))
return kNEUT;
else if (!type.compare("niwg_parameter"))
return kNIWG;
else if (!type.compare("nuwro_parameter"))
return kNUWRO;
else if (!type.compare("t2k_parameter"))
return kT2K;
else if (!type.compare("genie_parameter"))
return kGENIE;
else if (!type.compare("custom_parameter"))
return kCUSTOM;
else if (!type.compare("norm_parameter"))
return kNORM;
else if (!type.compare("likeweight_parameter"))
return kLIKEWEIGHT;
else if (!type.compare("spline_parameter"))
return kSPLINEPARAMETER;
else if (!type.compare("osc_parameter"))
return kOSCILLATION;
else if (!type.compare("modenorm_parameter"))
return kMODENORM;
else
return kUNKNOWN;
}
std::string FitBase::ConvDialType(int type) {
switch (type) {
case kNEUT: {
return "neut_parameter";
}
case kNIWG: {
return "niwg_parameter";
}
case kNUWRO: {
return "nuwro_parameter";
}
case kT2K: {
return "t2k_parameter";
}
case kGENIE: {
return "genie_parameter";
}
case kNORM: {
return "norm_parameter";
}
case kCUSTOM: {
return "custom_parameter";
}
case kLIKEWEIGHT: {
return "likeweight_parameter";
}
case kSPLINEPARAMETER: {
return "spline_parameter";
}
case kOSCILLATION: {
return "osc_parameter";
}
case kMODENORM: {
return "modenorm_parameter";
}
default:
return "unknown_parameter";
}
}
int FitBase::GetDialEnum(std::string const &type, std::string const &name) {
return FitBase::GetDialEnum(FitBase::ConvDialType(type), name);
}
int FitBase::GetDialEnum(int type, std::string const &name) {
int offset = type * 1000;
int this_enum = Reweight::kNoDialFound; // Not Found
- std::cout << "Getting dial enum " << type << " " << name << std::endl;
+ QLOG(DEB, "Getting dial enum " << type << " " << name);
// Select Types
switch (type) {
// NEUT DIAL TYPE
case kNEUT: {
#ifdef __NEUT_ENABLED__
int neut_enum = (int)neut::rew::NSyst::FromString(name);
if (neut_enum != 0) {
this_enum = neut_enum + offset;
}
#else
this_enum = Reweight::kNoTypeFound; // Not enabled
#endif
break;
}
// NIWG DIAL TYPE
case kNIWG: {
#ifdef __NIWG_ENABLED__
int niwg_enum = (int)niwg::rew::NIWGSyst::FromString(name);
if (niwg_enum != 0) {
this_enum = niwg_enum + offset;
}
#else
this_enum = Reweight::kNoTypeFound;
#endif
break;
}
// NUWRO DIAL TYPE
case kNUWRO: {
#ifdef __NUWRO_REWEIGHT_ENABLED__
int nuwro_enum = (int)nuwro::rew::NuwroSyst::FromString(name);
if (nuwro_enum > 0) {
this_enum = nuwro_enum + offset;
}
#else
this_enum = Reweight::kNoTypeFound;
#endif
}
// GENIE DIAL TYPE
case kGENIE: {
#ifdef __GENIE_ENABLED__
int genie_enum = (int)genie::rew::GSyst::FromString(name);
if (genie_enum > 0) {
this_enum = genie_enum + offset;
}
#else
this_enum = Reweight::kNoTypeFound;
#endif
break;
}
case kCUSTOM: {
int custom_enum = 0; // PLACEHOLDER
this_enum = custom_enum + offset;
break;
}
// T2K DIAL TYPE
case kT2K: {
#ifdef __T2KREW_ENABLED__
int t2k_enum = (int)t2krew::T2KSyst::FromString(name);
if (t2k_enum > 0) {
this_enum = t2k_enum + offset;
}
#else
this_enum = Reweight::kNoTypeFound;
#endif
break;
}
case kNORM: {
if (gNormEnums.find(name) == gNormEnums.end()) {
gNormEnums[name] = gNormEnums.size() + 1 + offset;
}
this_enum = gNormEnums[name];
break;
}
case kLIKEWEIGHT: {
if (gLikeWeightEnums.find(name) == gLikeWeightEnums.end()) {
gLikeWeightEnums[name] = gLikeWeightEnums.size() + 1 + offset;
}
this_enum = gLikeWeightEnums[name];
break;
}
case kSPLINEPARAMETER: {
if (gSplineParameterEnums.find(name) == gSplineParameterEnums.end()) {
gSplineParameterEnums[name] = gSplineParameterEnums.size() + 1 + offset;
}
this_enum = gSplineParameterEnums[name];
break;
}
case kOSCILLATION: {
#ifdef __PROB3PP_ENABLED__
int oscEnum = OscWeightEngine::SystEnumFromString(name);
if (oscEnum != 0) {
this_enum = oscEnum + offset;
}
#else
this_enum = Reweight::kNoTypeFound; // Not enabled
#endif
}
case kMODENORM: {
size_t us_pos = name.find_first_of('_');
std::string numstr = name.substr(us_pos + 1);
int mode_num = std::atoi(numstr.c_str());
LOG(FTL) << "Getting mode num " << mode_num << std::endl;
if (!mode_num) {
ERR(FTL) << "Attempting to parse dial name: \"" << name
<< "\" as a mode norm dial but failed." << std::endl;
throw;
}
this_enum = 60 + mode_num + offset;
break;
}
}
// If Not Enabled
if (this_enum == Reweight::kNoTypeFound) {
ERR(FTL) << "RW Engine not supported for " << FitBase::ConvDialType(type)
<< std::endl;
ERR(FTL) << "Check dial " << name << std::endl;
}
// If Not Found
if (this_enum == Reweight::kNoDialFound) {
ERR(FTL) << "Dial " << name << " not found." << std::endl;
}
return this_enum;
}
int Reweight::ConvDialType(std::string const &type) {
return FitBase::ConvDialType(type);
}
std::string Reweight::ConvDialType(int type) {
return FitBase::ConvDialType(type);
}
int Reweight::GetDialType(int type) {
int t = (type / 1000);
return t > kMODENORM ? Reweight::kNoDialFound : t;
}
int Reweight::RemoveDialType(int type) { return (type % 1000); }
int Reweight::NEUTEnumFromName(std::string const &name) {
#ifdef __NEUT_ENABLED__
int neutenum = (int)neut::rew::NSyst::FromString(name);
return (neutenum > 0) ? neutenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::NIWGEnumFromName(std::string const &name) {
#ifdef __NIWG_ENABLED__
int niwgenum = (int)niwg::rew::NIWGSyst::FromString(name);
return (niwgenum != 0) ? niwgenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::NUWROEnumFromName(std::string const &name) {
#ifdef __NUWRO_REWEIGHT_ENABLED__
int nuwroenum = (int)nuwro::rew::NuwroSyst::FromString(name);
return (nuwroenum > 0) ? nuwroenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::GENIEEnumFromName(std::string const &name) {
#ifdef __GENIE_ENABLED__
int genieenum = (int)genie::rew::GSyst::FromString(name);
return (genieenum > 0) ? genieenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::T2KEnumFromName(std::string const &name) {
#ifdef __T2KREW_ENABLED__
int t2kenum = (int)t2krew::T2KSyst::FromString(name);
return (t2kenum > 0) ? t2kenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::OscillationEnumFromName(std::string const &name) {
#ifdef __PROB3PP_ENABLED__
int oscEnum = OscWeightEngine::SystEnumFromString(name);
return (oscEnum > 0) ? oscEnum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::NUISANCEEnumFromName(std::string const &name, int type) {
int nuisenum = Reweight::DialList().EnumFromNameAndType(name, type);
return nuisenum;
}
int Reweight::CustomEnumFromName(std::string const &name) {
int custenum = Reweight::ConvertNUISANCEDial(name);
return custenum;
}
int Reweight::ConvDial(std::string const &name, std::string const &type,
bool exceptions) {
return Reweight::ConvDial(name, Reweight::ConvDialType(type), exceptions);
}
int Reweight::ConvDial(std::string const &fullname, int type, bool exceptions) {
std::string name =
GeneralUtils::ParseToStr(fullname, ",")[0]; // Only use first dial given
// Produce offset seperating each type.
int offset = type * 1000;
int genenum = Reweight::kNoDialFound;
switch (type) {
case kNEUT:
genenum = NEUTEnumFromName(name);
break;
case kNIWG:
genenum = NIWGEnumFromName(name);
break;
case kNUWRO:
genenum = NUWROEnumFromName(name);
break;
case kGENIE:
genenum = GENIEEnumFromName(name);
break;
case kT2K:
genenum = T2KEnumFromName(name);
break;
case kCUSTOM:
genenum = CustomEnumFromName(name);
break;
case kNORM:
case kLIKEWEIGHT:
case kSPLINEPARAMETER:
case kNEWSPLINE:
genenum = NUISANCEEnumFromName(name, type);
break;
case kOSCILLATION:
genenum = OscillationEnumFromName(name);
break;
case kMODENORM:
genenum = ModeNormEngine::SystEnumFromString(name);
break;
default:
genenum = Reweight::kNoTypeFound;
break;
}
// Throw if required.
if (exceptions) {
// If Not Enabled
if (genenum == Reweight::kGeneratorNotBuilt) {
ERR(FTL) << "RW Engine not supported for " << FitBase::ConvDialType(type)
<< std::endl;
ERR(FTL) << "Check dial " << name << std::endl;
throw;
}
// If no type enabled
if (genenum == Reweight::kNoTypeFound) {
ERR(FTL) << "Type mismatch inside ConvDialEnum" << std::endl;
throw;
}
// If Not Found
if (genenum == Reweight::kNoDialFound) {
ERR(FTL) << "Dial " << name << " not found." << std::endl;
throw;
}
}
// Add offset if no issue
int nuisenum = genenum;
if ((genenum != Reweight::kGeneratorNotBuilt) &&
(genenum != Reweight::kNoTypeFound) &&
(genenum != Reweight::kNoDialFound)) {
nuisenum += offset;
}
// Now register dial
Reweight::DialList().RegisterDialEnum(name, type, nuisenum);
return nuisenum;
}
std::string Reweight::ConvDial(int nuisenum) {
for (size_t i = 0; i < Reweight::DialList().fAllDialEnums.size(); i++) {
if (Reweight::DialList().fAllDialEnums[i] == nuisenum) {
return Reweight::DialList().fAllDialNames[i];
}
}
LOG(FIT) << "Cannot find dial with enum = " << nuisenum << std::endl;
return "";
}
diff --git a/src/Reweight/weightRPA.h b/src/Reweight/weightRPA.h
index 6f969fd..c44acdd 100644
--- a/src/Reweight/weightRPA.h
+++ b/src/Reweight/weightRPA.h
@@ -1,414 +1,409 @@
#ifndef weightRPA_h
#define weightRPA_h
#include <fstream> //ifstream
#include <iostream> //cout
#include <TString.h>
#include <TH2D.h>
#include <TH1D.h>
#include <TFile.h>
#include "math.h"
#include "assert.h"
//For Compatibility with ROOT compiler
//uncomment the following:
//#define ROOT
/*!
* Code example to extract the RPA effect central weight
* and its uncertainties from the prepared files.
* Heidi Schellman (Oregon State) and Rik Gran (Minnesota Duluth)
* for use in MINERvA experiment analysis
* must compile with the ROOT libraries
* g++ `root-config --glibs --cflags` -O3 weightRPAtest.cxx -o weightRPA
* The underlying model is from the IFIC Valencia group
* see (public) minerva docdb:12688 for the full physics discussion
*/
// NOTE UNITS ARE GEV in the calculation
// make sure you convert MeV to GeV before calling these functions
// Class for getting RPA paramaters inputs from a given file
// Includes methods that return all five RPA weights at once
// (is the most cpu efficient way to get them)
// Or return just the central value
// (skipping the uncertainties code completely if only cv wanted)
// Or return each CV and uncertainty one at a time
// (is the least cpu efficient, repeats some calculations 5 times)
class weightRPA {
public:
//Constructor: Read in params from a filename
weightRPA(const TString f) { read(f); } //Read in params from file
TString filename;
TFile* fRPAratio;
TH2D *hRPArelratio;
TH2D *hRPAnonrelratio;
TH1D *hQ2relratio;
TH1D *hQ2nonrelratio;
TArrayD *TADrpapolyrel;
Double_t *rpapolyrel ;
TArrayD *TADrpapolynonrel;
Double_t *rpapolynonrel;
static const int CENTRAL=0;
static const int LOWQ2 = 1;
static const int HIGHQ2 = 2;
// true to take from histogram, false use parameterization
static const bool Q2histORparam = true;
// MINERvA holds kinematics in MeV, but all these functions require GeV
// So make sure you pass them in GeV.
inline double getWeight(const double mc_q0, const double mc_q3, double * weights); //in GeV
inline double getWeight(const double mc_q0, const double mc_q3); //in GeV
inline double getWeight(const double mc_q0, const double mc_q3, int type, int sign); //in GeV
inline double getWeight(const double Q2); //in GeV^2
inline double getWeightLowQ2(const double mc_q0, const double mc_q3, const int sign);
inline double getWeightHighQ2(const double mc_q0, const double mc_q3, const int sign);
inline double getWeightQ2(const double mc_Q2, const bool relORnonrel=true);
//Initializer
inline void read(const TString f);
// q0 and q3 in GeV, type = 1 for low Q2, 2 for high Q2, 0 for central
//double getWeightInternal(const double mc_q0, const double mc_q3,int type, int sign);
private:
inline double getWeightInternal(const double mc_Q2);
inline double getWeightInternal(const double mc_q0, const double mc_q3, const int type, const int sign);
inline double getWeightInternal(const double mc_q0, const double mc_q3, double *weights=0);
inline double getWeightQ2parameterization(const double mc_Q2, const bool relORnonrel);
inline double getWeightQ2fromhistogram(const double mc_Q2, const bool relORnonrel);
};
void weightRPA::read(const TString f)
//Read in the params doubles from a file
//argument: valid filename
{
fRPAratio = TFile::Open(f,"READONLY");
if (fRPAratio){
hRPArelratio = (TH2D*)fRPAratio->Get("hrelratio");
hRPAnonrelratio = (TH2D*)fRPAratio->Get("hnonrelratio");
hQ2relratio = (TH1D*)fRPAratio->Get("hQ2relratio");
hQ2nonrelratio = (TH1D*)fRPAratio->Get("hQ2nonrelratio");
TADrpapolyrel = (TArrayD*)fRPAratio->Get("rpapolyrel");
rpapolyrel = TADrpapolyrel->GetArray();
TADrpapolynonrel = (TArrayD*)fRPAratio->Get("rpapolynonrel");
rpapolynonrel = TADrpapolynonrel->GetArray();
hRPArelratio->Print();
std::cout << "have read in ratios from file " << f <<std::endl;
-
- }
- else{
+ } else {
//File could not be read
std::cout << "File could not be read" << std::endl;
-
}
}
double weightRPA::getWeightInternal(const double mc_q0, const double mc_q3, double *weights){
assert(hRPArelratio);
if (weights != 0){
for (int i = 0; i < 5; i++){
weights[i] = 1.0;
}
}
// the construction here is that the histogram bins
// line up in mev-sized steps e.g. from 0.018 to 0.019
// and the value stored is the value at bin-center.
// double gevmev = 0.001 ;
const Double_t gevlimit = 3.; // upper limit for 2d
Int_t rpamevlimit = gevlimit*1000.;
//Double_t Q2gev = mc_Q2*gevmev*gevmev;
Double_t Q2gev = mc_q3*mc_q3-mc_q0*mc_q0;
Int_t q3bin = Int_t(mc_q3*1000.);
Int_t q0bin = Int_t(mc_q0*1000.);
if(mc_q0 >= gevlimit) q0bin = rpamevlimit - 1;
if(mc_q3 >= gevlimit) q3bin = rpamevlimit - 1;
// Nieves does not calculate anything below binding energy.
// I don't know what GENIE does, but lets be soft about this.
// Two things lurking here at once.
// One, we are taking out a 10 MeV offset between GENIE and Valencia.
// Second, we are protecting against being asked for a weight that is too low in q0.
// It actually shouldn't happen for real GENIE events,
// but this protection does something that doesn't suck, just in case.
// you would see the artifact in a plot for sure, but better than writing 1.0.
Int_t q0offsetValenciaGENIE = 10;
if(mc_q0 < 0.018) q0bin = 18+q0offsetValenciaGENIE;
Double_t thisrwtemp = hRPArelratio->GetBinContent(q3bin,q0bin-q0offsetValenciaGENIE);
// now trap bogus entries. Not sure why they happen, but set to 1.0 not 0.0
if(thisrwtemp <= 0.001)thisrwtemp = 1.0;
// events in genie but not in valencia should get a weight
// related to a similar q0 from the bulk distribution.
if(mc_q0 < 0.15 && thisrwtemp > 0.9){
thisrwtemp = hRPArelratio->GetBinContent(q3bin+150, q0bin-q0offsetValenciaGENIE);
}
//Double_t *mypoly;
//mypoly = rpapoly;
if(Q2gev >= 9.0){
thisrwtemp = 1.0;
}
else if(Q2gev > 3.0) {
// hiding option, use the Q2 parameterization everywhere
// } else if(Q2gev > 3.0 || rwRPAQ2) {
// turn rwRPAQ2 all the way on to override the 2D histogram
// illustrates the old-style Q2 suppression many folks still use.
thisrwtemp = getWeightQ2(Q2gev,true);
// double powerQ2 = 1.0;
//thisrwtemp = 0.0;
//for(int ii=0; ii<10; ii++){
// thisrwtemp += rpapoly[ii]*powerQ2;
// powerQ2 *= Q2gev;
//}
//std::cout << "test temp " << thisrwtemp << " " << rpamypoly[2] << std::endl;
}
if(!(thisrwtemp >= 0.001 && thisrwtemp <= 2.0))thisrwtemp = 1.0;
// hiding option, turn off the enhancement.
//if(rwoffSRC && thisrwtemp > 1.0)thisrwtemp = 1.0;
if (0 == weights) return thisrwtemp;
// if this was called without passing an array,
// the user didn't want us to calculate the +/- 1-sigma bounds
// so the above line returned before even trying.
weights[0] = thisrwtemp;
//if (type == 0) return thisrwtemp;
//if (type == 1) {
// Construct the error bands on the low Q2 suppression.
// Double_t thisrwSupP1 = 1.0;
// Double_t thisrwSupM1 = 1.0;
if( thisrwtemp < 1.0){
// make the suppression stronger or weaker to muon capture uncertainty
// rwRPAonesig is either +1 or -1, which is 0.25 (25%).
// possible to be re-written to produce 2 and 3 sigma.
weights[1] = thisrwtemp + 1.0 * (0.25)*(1.0 - thisrwtemp);
weights[2] = thisrwtemp - 1.0 * (0.25)*(1.0 - thisrwtemp);
}
else{
weights[1] = thisrwtemp;
weights[2] = thisrwtemp;
}
//std::cout << "check " << thisrwtemp << " " << weights[1] << " " << weights[2] << std::endl;
// Construct the rest of the error bands on the low Q2 suppression.
// this involves getting the weight from the non-relativistic ratio
//if (type == 2){
Double_t thisrwEnhP1 = 1.0;
Double_t thisrwEnhM1 = 1.0;
// make enhancement stronger or weaker to Federico Sanchez uncertainty
// this does NOT mean two sigma, its overloading the option.
Double_t thisrwextreme = hRPAnonrelratio->GetBinContent(q3bin,q0bin-q0offsetValenciaGENIE);
// now trap bogus entries. Not sure why they happen, but set to 1.0 not 0.0
if(thisrwextreme <= 0.001)thisrwextreme = 1.0;
if(mc_q0 < 0.15 && thisrwextreme > 0.9){
thisrwextreme = hRPAnonrelratio->GetBinContent(q3bin+150, q0bin-q0offsetValenciaGENIE);
}
//std::cout << "ext " << thisrwextreme << " " << thisrwtemp << std::endl;
// get the same for the Q2 dependent thing,
// but from the nonrelativistic polynomial
if(Q2gev >= 9.0){
thisrwextreme = 1.0;
}
else if(Q2gev > 3.0 ) {
thisrwextreme = getWeightQ2(Q2gev,false);
//double powerQ2 = 1.0;
//thisrwextreme = 0.0;
//for(int ii=0; ii<10; ii++){
// thisrwextreme += rpapolynonrel[ii]*powerQ2;
// powerQ2 *= Q2gev;
//}
//std::cout << "test extreme " << thisrwextreme << " " << mypolynonrel[2] << std::endl;
}
if(!(thisrwextreme >= 0.001 && thisrwextreme <= 2.0))thisrwextreme = 1.0;
//std::cout << "test extreme " << Q2gev << " " << thisrwextreme << " " << thisrwtemp << std::endl;
Double_t RelToNonRel = 0.6;
// based on some distance between central value and extreme
thisrwEnhP1 = thisrwtemp + RelToNonRel * (thisrwextreme-thisrwtemp);
Double_t thisrwEnhP1max = thisrwextreme;
if(Q2gev < 0.9)thisrwEnhP1 += 1.5*(0.9 - Q2gev)*(thisrwEnhP1max - thisrwEnhP1);
// sanity check, don't let the upper error bound go above the nonrel limit.
if(thisrwEnhP1 > thisrwEnhP1max)thisrwEnhP1 = thisrwEnhP1max;
// don't let positive error bound be closer than 3% above the central value
// will happen at very high Q2 and very close to Q2 = 0
if(thisrwEnhP1 < thisrwtemp + 0.03)thisrwEnhP1 = thisrwtemp + 0.03;
thisrwEnhM1 = thisrwtemp - RelToNonRel * (thisrwextreme-thisrwtemp);
// don't let negative error bound be closer than 3% below the central value
if(thisrwEnhM1 > thisrwtemp - 0.03)thisrwEnhM1 = thisrwtemp - 0.03;
// even still, don't let the lower error bound go below 1.0 at high-ish Q2
if(Q2gev > 1.0 && thisrwEnhM1 < 1.0)thisrwEnhM1 = 1.0;
// whew. so now return the main weight
// and return the array of all five weights in some array
// thisrwtemp, thisrwSupP1, thisrwSupM1, thisrwEnhP1, thisrwEnhM1
//if (sign == 1) return thisrwEnhP1;
//if (sign == -1) return thisrwEnhM1;
weights[3] = thisrwEnhP1;
weights[4] = thisrwEnhM1;
// still return the central value
return thisrwtemp;
}
double weightRPA::getWeight(const double mc_q0, const double mc_q3){
return getWeightInternal(mc_q0, mc_q3);
}
double weightRPA::getWeight(const double mc_q0, const double mc_q3, double *weights){
return getWeightInternal(mc_q0, mc_q3, weights);
}
double weightRPA::getWeight(const double mc_q0, const double mc_q3, int type, int sign){
return getWeightInternal(mc_q0, mc_q3, type, sign);
}
double weightRPA::getWeight(const double mc_Q2){
return getWeightQ2(mc_Q2);
}
double weightRPA::getWeightInternal(const double mc_q0, const double mc_q3, int type, int sign){
double weights[5] = {1., 1., 1., 1., 1.};
double cv = getWeightInternal(mc_q0, mc_q3, weights);
if(type==0)return cv;
else if(type==weightRPA::LOWQ2 && sign == 1)return weights[1];
else if(type==weightRPA::LOWQ2 && sign == -1)return weights[2];
else if(type==weightRPA::HIGHQ2 && sign == 1)return weights[3];
else if(type==weightRPA::HIGHQ2 && sign == -1)return weights[4];
//else {
// // should never happen? Bork ?
// return cv; //?
//}
return cv;
}
double weightRPA::getWeightQ2(const double mc_Q2, const bool relORnonrel){
if(mc_Q2 < 0.0)return 1.0; // this is Q2 actually, not sure hw
if(mc_Q2 > 9.0)return 1.0;
// this function needs to know two options.
// does user want rel (cv) or nonrel
// does user want to use the histogram or parameterization
if(Q2histORparam)return getWeightQ2fromhistogram(mc_Q2, relORnonrel);
else return getWeightQ2parameterization(mc_Q2, relORnonrel);
}
double weightRPA::getWeightQ2parameterization(const double mc_Q2, const bool relORnonrel){
if(mc_Q2 < 0.0)return 1.0;
if(mc_Q2 > 9.0)return 1.0;
// this one returns just the polynomial Q2 version
// for special tests. Poor answer for baseline MINERvA QE events.
// No uncertainty assigned to this usecase at this time.
//double gevmev = 0.001; // minerva sends in MeV.
double Q2gev = mc_Q2;
double powerQ2 = 1.0;
double thisrwtemp = 0.0;
thisrwtemp = 0.0;
for(int ii=0; ii<10; ii++){
if(relORnonrel)thisrwtemp += rpapolyrel[ii]*powerQ2;
else thisrwtemp += rpapolynonrel[ii]*powerQ2;
powerQ2 *= Q2gev;
}
return thisrwtemp;
}
double weightRPA::getWeightQ2fromhistogram(const double mc_Q2, const bool relORnonrel){
if(mc_Q2 < 0.0)return 1.0;
if(mc_Q2 > 9.0) return 1.0;
if(relORnonrel)return hQ2relratio->GetBinContent( hQ2relratio->FindBin(mc_Q2) );
else return hQ2nonrelratio->GetBinContent( hQ2nonrelratio->FindBin(mc_Q2) );
// interpolation might be overkill for such a finely binned histogram, 0.01%
// but the extra cpu cycles maybe small.
// save it here for some future use.
//if(relORnonrel)return hQ2relratio->Interpolate(mc_Q2);
//else return hQ2nonrelratio->Interpolate(mc_Q2);
}
double weightRPA::getWeightLowQ2(const double mc_q0, const double mc_q3, int sign){
return getWeightInternal(mc_q0,mc_q3,weightRPA::LOWQ2,sign);
}
double weightRPA::getWeightHighQ2(const double mc_q0, const double mc_q3, int sign){
return getWeightInternal(mc_q0,mc_q3,weightRPA::HIGHQ2,sign);
}
#endif
-
-

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 5:28 PM (9 h, 37 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4016321
Default Alt Text
(32 KB)

Event Timeline