From a1704165bc91d1c2dcfde14efaa123bebc5ad269 Mon Sep 17 00:00:00 2001 From: Luke Pickering Date: Fri, 18 Oct 2024 17:42:41 +0100 Subject: [PATCH] adds mysillyvar example demo of user-definable projection operators - a very simple proof of concept this is not integrated at all and meant to display the user API and not how it would be implemented, comments welcome on the concept --- .../Samples/SamplePDFDune_FHC_numuselec.yaml | 10 +- samplePDFDUNE/samplePDFDUNEBeamFD.cpp | 219 +++++++++++------- samplePDFDUNE/samplePDFDUNEBeamFD.h | 18 +- src/EventRates.cpp | 9 +- 4 files changed, 163 insertions(+), 93 deletions(-) diff --git a/configs/Samples/SamplePDFDune_FHC_numuselec.yaml b/configs/Samples/SamplePDFDune_FHC_numuselec.yaml index a2df937..01552de 100644 --- a/configs/Samples/SamplePDFDune_FHC_numuselec.yaml +++ b/configs/Samples/SamplePDFDune_FHC_numuselec.yaml @@ -26,12 +26,12 @@ Binning: - VarStr: "ELepRec" Uniform: [25,0,5] Title: "E_{lep.}^{Rec} [GeV]" - - VarStr: "EHadRec" + - VarStr: "mysillyvar" Uniform: [25,0,5] - Title: "E_{had.}^{Vis} [GeV]" - - VarStr: "RecoNeutrinoEnergy" - Uniform: [10,0.5,5.5] - Title: "E_{#nu}^{Rec} [GeV]" + Title: "dumbdumb" + # - VarStr: "RecoNeutrinoEnergy" + # Uniform: [10,0.5,5.5] + # Title: "E_{#nu}^{Rec} [GeV]" # if you want to use the generic binning with a single # axis, useful for calculated VarStr that cannot be referenced ForceGeneric: True diff --git a/samplePDFDUNE/samplePDFDUNEBeamFD.cpp b/samplePDFDUNE/samplePDFDUNEBeamFD.cpp index 6f79c89..77f4842 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamFD.cpp +++ b/samplePDFDUNE/samplePDFDUNEBeamFD.cpp @@ -553,11 +553,12 @@ TH1D* samplePDFDUNEBeamFD::get1DVarHist(KinematicTypes Var1,std::vector< std::ve return _h1DVar; } -double const& samplePDFDUNEBeamFD::ReturnKinematicParameterByReference(int KinematicParameter, int iSample, int iEvent) { +double const &samplePDFDUNEBeamFD::ReturnKinematicParameterByReference( + int KinematicParameter, int iSample, int iEvent) { - switch(KinematicParameter){ + switch (KinematicParameter) { case kTrueNeutrinoEnergy: - return dunemcSamples[iSample].rw_etru[iEvent]; + return dunemcSamples[iSample].rw_etru[iEvent]; case kRecoNeutrinoEnergy: return dunemcSamples[iSample].rw_erec_shifted[iEvent]; case kTrueXPos: @@ -580,98 +581,128 @@ double const& samplePDFDUNEBeamFD::ReturnKinematicParameterByReference(int Kinem case kq3: return dunemcSamples[iSample].true_q3[iEvent]; default: - std::stringstream ss; - ss << "[ERROR]: " << __FILE__ << ":" << __LINE__ - << " ReturnKinematicParameterByReference Did not recognise " - "Kinematic Parameter type:" - << KinematicParameter - << ". Is it possibly only available via ReturnKinematicParameter " - "(not by reference?) if you need it here, give it storage in " - "dunemc_base and move it to here." - ; - throw std::runtime_error(ss.str()); - } - } - - double samplePDFDUNEBeamFD::ReturnKinematicParameter(int KinematicParameter, - int iSample, - int iEvent) { - switch (KinematicParameter) { - case kERecQE: { - constexpr double V = 0; // 0 binding energy for now - constexpr double mn = 939.565; // neutron mass - constexpr double mp = 938.272; // proton mass - - double mN_eff = mn - V; - double mN_oth = mp; - - if (dunemcSamples[iSample].rw_nuPDGunosc[iEvent] < - 0) { // if anti-neutrino, swap target/out masses - mN_eff = mp - V; - mN_oth = mn; - } - double el = dunemcSamples[iSample].rw_erec_lep[iEvent]; + if (KinematicParameter >= kNDefaultProjections) { + MACH3LOG_ERROR("ReturnKinematicParameterByReference passed " + "KinematicParameter: {}, which appears to refer to user " + "projection: {}, but user projections cannot be evaluated " + "by reference.", + KinematicParameter, + KinematicParameter - kNDefaultProjections); + throw MaCh3Exception(__FILE__, __LINE__); + } - // this is funky, but don't be scared, it defines an annonymous function - // in place that grabs the lepton mass in MeV when given the neutrino PDG - // and whether the interaction was CC or NC and then immediately calls it. - // It's basically a generalisation of the ternary operator. - double ml = - [](int nupdg, bool isCC) { - switch (std::abs(nupdg)) { - case 12: { - return isCC ? 0.511 : 0; - } - case 14: { - return isCC ? 105.66 : 0; - } - case 16: { - return isCC ? 1777.0 : 0; - } - } - }(dunemcSamples[iSample].rw_nuPDGunosc[iEvent], - dunemcSamples[iSample].rw_isCC[iEvent]); + MACH3LOG_ERROR( + "ReturnKinematicParameterByReference Did not recognise " + "Kinematic Parameter type: {}. Is it possibly only available via " + "ReturnKinematicParameter " + "(not by reference?) if you need it here, give it storage in " + "dunemc_base and move it to here.", + KinematicParameter); + throw MaCh3Exception(__FILE__, __LINE__); + } +} - double pl = std::sqrt(el*el - ml*ml); // momentum of lepton +double samplePDFDUNEBeamFD::ReturnKinematicParameter(int KinematicParameter, + int iSample, int iEvent) { - double rEnu = - (2 * mN_eff * el - ml * ml + mN_oth * mN_oth - mN_eff * mN_eff) / - (2 * (mN_eff - el + - pl * std::cos(dunemcSamples[iSample].rw_theta[iEvent]))); + if (KinematicParameter >= kNDefaultProjections) { - return rEnu; + if ((KinematicParameter - kNDefaultProjections) > user_projections.size()) { + MACH3LOG_ERROR( + "Passed KinematicParameter: {}, which appears to refer to user " + "projection: {}, but we only have {} user projections defined.", + KinematicParameter, KinematicParameter - kNDefaultProjections, + user_projections.size()); + throw MaCh3Exception(__FILE__, __LINE__); } - case kEHadRec: { - return dunemcSamples[iSample].rw_eRecoP[iEvent] + - dunemcSamples[iSample].rw_eRecoPip[iEvent] + - dunemcSamples[iSample].rw_eRecoPim[iEvent] + - dunemcSamples[iSample].rw_eRecoPi0[iEvent] + - dunemcSamples[iSample].rw_eRecoN[iEvent]; - } - default: { - return ReturnKinematicParameterByReference(KinematicParameter, iSample, - iEvent); - } - } + return user_projections[KinematicParameter - kNDefaultProjections].proj( + dunemcSamples[iSample], iEvent); + } + + switch (KinematicParameter) { + case kERecQE: { + constexpr double V = 0; // 0 binding energy for now + constexpr double mn = 939.565; // neutron mass + constexpr double mp = 938.272; // proton mass + + double mN_eff = mn - V; + double mN_oth = mp; + + if (dunemcSamples[iSample].rw_nuPDGunosc[iEvent] < + 0) { // if anti-neutrino, swap target/out masses + mN_eff = mp - V; + mN_oth = mn; + } + + double el = dunemcSamples[iSample].rw_erec_lep[iEvent]; + + // this is funky, but don't be scared, it defines an annonymous function + // in place that grabs the lepton mass in MeV when given the neutrino PDG + // and whether the interaction was CC or NC and then immediately calls it. + // It's basically a generalisation of the ternary operator. + double ml = + [](int nupdg, bool isCC) { + switch (std::abs(nupdg)) { + case 12: { + return isCC ? 0.511 : 0; + } + case 14: { + return isCC ? 105.66 : 0; + } + case 16: { + return isCC ? 1777.0 : 0; + } + } + }(dunemcSamples[iSample].rw_nuPDGunosc[iEvent], + dunemcSamples[iSample].rw_isCC[iEvent]); + + double pl = std::sqrt(el * el - ml * ml); // momentum of lepton + + double rEnu = + (2 * mN_eff * el - ml * ml + mN_oth * mN_oth - mN_eff * mN_eff) / + (2 * (mN_eff - el + + pl * std::cos(dunemcSamples[iSample].rw_theta[iEvent]))); + + return rEnu; } + case kEHadRec: { + + return dunemcSamples[iSample].rw_eRecoP[iEvent] + + dunemcSamples[iSample].rw_eRecoPip[iEvent] + + dunemcSamples[iSample].rw_eRecoPim[iEvent] + + dunemcSamples[iSample].rw_eRecoPi0[iEvent] + + dunemcSamples[iSample].rw_eRecoN[iEvent]; + } + default: { + return ReturnKinematicParameterByReference(KinematicParameter, iSample, + iEvent); + } + } +} int samplePDFDUNEBeamFD::ReturnKinematicParameterFromString(std::string KinematicParameterStr){ - if (KinematicParameterStr.find("TrueNeutrinoEnergy") != std::string::npos) {return kTrueNeutrinoEnergy;} - if (KinematicParameterStr.find("RecoNeutrinoEnergy") != std::string::npos) {return kRecoNeutrinoEnergy;} - if (KinematicParameterStr.find("TrueXPos") != std::string::npos) {return kTrueXPos;} - if (KinematicParameterStr.find("TrueYPos") != std::string::npos) {return kTrueYPos;} - if (KinematicParameterStr.find("TrueZPos") != std::string::npos) {return kTrueZPos;} - if (KinematicParameterStr.find("CVNNumu") != std::string::npos) {return kCVNNumu;} - if (KinematicParameterStr.find("CVNNue") != std::string::npos) {return kCVNNue;} - if (KinematicParameterStr.find("M3Mode") != std::string::npos) {return kM3Mode;} - if (KinematicParameterStr.find("global_bin_number") != std::string::npos) {return kGlobalBinNumber;} - if (KinematicParameterStr.find("q0") != std::string::npos) {return kq0;} - if (KinematicParameterStr.find("q3") != std::string::npos) {return kq3;} - if (KinematicParameterStr.find("ERecQE") != std::string::npos) {return kERecQE;} - if (KinematicParameterStr.find("ELepRec") != std::string::npos) {return kELepRec;} - if (KinematicParameterStr.find("EHadRec") != std::string::npos) {return kEHadRec;} + if (KinematicParameterStr == "TrueNeutrinoEnergy") {return kTrueNeutrinoEnergy;} + if (KinematicParameterStr == "RecoNeutrinoEnergy") {return kRecoNeutrinoEnergy;} + if (KinematicParameterStr == "TrueXPos") {return kTrueXPos;} + if (KinematicParameterStr == "TrueYPos") {return kTrueYPos;} + if (KinematicParameterStr == "TrueZPos") {return kTrueZPos;} + if (KinematicParameterStr == "CVNNumu") {return kCVNNumu;} + if (KinematicParameterStr == "CVNNue") {return kCVNNue;} + if (KinematicParameterStr == "M3Mode") {return kM3Mode;} + if (KinematicParameterStr == "global_bin_number") {return kGlobalBinNumber;} + if (KinematicParameterStr == "q0") {return kq0;} + if (KinematicParameterStr == "q3") {return kq3;} + if (KinematicParameterStr == "ERecQE") {return kERecQE;} + if (KinematicParameterStr == "ELepRec") {return kELepRec;} + if (KinematicParameterStr == "EHadRec") {return kEHadRec;} + + for(size_t up_it = 0; up_it < user_projections.size(); ++ up_it){ + if(KinematicParameterStr == user_projections[up_it].name){ + return kNDefaultProjections + up_it; + } + } std::stringstream ss; ss << "[ERROR]: " << __FILE__ << ":" << __LINE__ @@ -842,3 +873,21 @@ std::vector samplePDFDUNEBeamFD::ReturnKinematicParameterBinning(int Kin return binningVector; } + +int samplePDFDUNEBeamFD::AddProjection( + std::string const &name, + std::function proj) { + for (auto const &up : user_projections) { + if (name == up.name) { + MACH3LOG_ERROR("Attempting to overwrite existing UserProjection: {}, " + "please pick a new name", + name); + throw MaCh3Exception(__FILE__, __LINE__); + } + } + user_projections.push_back(UserProjection{name, proj}); + return kNDefaultProjections + (user_projections.size() - 1); +} + +std::vector samplePDFDUNEBeamFD::user_projections; + diff --git a/samplePDFDUNE/samplePDFDUNEBeamFD.h b/samplePDFDUNE/samplePDFDUNEBeamFD.h index 0f86875..eadcf5d 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamFD.h +++ b/samplePDFDUNE/samplePDFDUNEBeamFD.h @@ -26,7 +26,7 @@ class samplePDFDUNEBeamFD : virtual public samplePDFFDBase ~samplePDFDUNEBeamFD(); enum KinematicTypes { - kTrueNeutrinoEnergy, + kTrueNeutrinoEnergy = 0, kRecoNeutrinoEnergy, kTrueXPos, kTrueYPos, @@ -40,14 +40,28 @@ class samplePDFDUNEBeamFD : virtual public samplePDFFDBase kq3, kERecQE, kELepRec, - kEHadRec + kEHadRec, + kNDefaultProjections }; //More robust getters to make plots in different variables, mode, osc channel, systematic weighting and with bin range TH1D* get1DVarHist(KinematicTypes Var1, int fModeToFill=-1, int fSampleToFill=-1, int WeightStyle=0, TAxis* Axis=0); TH1D* get1DVarHist(KinematicTypes Var1, std::vector< std::vector > Selection, int WeightStyle=0, TAxis* Axis=0); + + struct UserProjection { + std::string name; + std::function proj; + }; + + // add a user projection, capture the return value if you don't want to look the + // KinematiceParameter value up by string + static int AddProjection(std::string const &, std::function); + protected: + + static std::vector user_projections; + void Init(); int setupExperimentMC(int iSample); void setupFDMC(int iSample); diff --git a/src/EventRates.cpp b/src/EventRates.cpp index 42f03d0..3bb5061 100644 --- a/src/EventRates.cpp +++ b/src/EventRates.cpp @@ -15,6 +15,7 @@ #include "samplePDF/GenericBinningTools.h" #include "samplePDFDUNE/MaCh3DUNEFactory.h" +#include "samplePDFDUNE/samplePDFDUNEBeamFD.h" void Write1DHistogramsToFile(std::string OutFileName, std::vector Histograms) { @@ -66,6 +67,13 @@ int main(int argc, char *argv[]) { // #################################################################################### // Create samplePDFFD objects + samplePDFDUNEBeamFD::AddProjection( + "mysillyvar", [](dunemc_base const &dunemcSample, int iEvent) { + return (dunemcSample.rw_erec_shifted[iEvent] - + dunemcSample.rw_erec_lep[iEvent]) / + (dunemcSample.true_q0[iEvent]); + }); + std::vector DUNEPdfs; MakeMaCh3DuneInstance(fitMan.get(), DUNEPdfs, xsec, osc); @@ -114,7 +122,6 @@ int main(int argc, char *argv[]) { gc1->Print("GenericBinTest.pdf]"); - std::string OutFileName = GetFromManager( fitMan->raw()["General"]["OutputFile"], "EventRatesOutput.root");