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 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 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 //ifstream #include //cout #include #include #include #include #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 <= 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 - -