diff --git a/covariance/covarianceXsec.cpp b/covariance/covarianceXsec.cpp index 602d925e..f8353185 100644 --- a/covariance/covarianceXsec.cpp +++ b/covariance/covarianceXsec.cpp @@ -1,4 +1,5 @@ #include "covariance/covarianceXsec.h" +#include "samplePDF/Structs.h" // ******************************************** // ETA - YAML constructor @@ -114,10 +115,29 @@ const std::vector covarianceXsec::GetSplineParsNamesFromDetID(const for (auto &pair : _fSystToGlobalSystIndexMap[SystType::kSpline]) { auto &SplineIndex = pair.first; auto &SystIndex = pair.second; - if ((GetParDetID(SystIndex) & DetID )){ + if (AppliesToDetID(SystIndex, DetID)) { //If parameter applies to required DetID returnVec.push_back(_fSplineNames.at(SplineIndex)); } + + } + return returnVec; +} + +const std::vector covarianceXsec::GetSplineInterpolationFromDetID(const int DetID) { + std::vector returnVec; + for (auto &pair : _fSystToGlobalSystIndexMap[SystType::kSpline]) { + auto &SplineIndex = pair.first; + auto &SystIndex = pair.second; + + if (AppliesToDetID(SystIndex, DetID)) { //If parameter applies to required DetID + returnVec.push_back(SplineParams.at(SplineIndex)._SplineInterpolationType); + } + + // if ((GetParDetID(SystIndex) & DetID )){ + // returnVec.push_back(SplineParams.at(SplineIndex)._SplineInterpolationType); + // } } + return returnVec; } diff --git a/covariance/covarianceXsec.h b/covariance/covarianceXsec.h index 1dd61aff..b6fe0ea9 100644 --- a/covariance/covarianceXsec.h +++ b/covariance/covarianceXsec.h @@ -32,6 +32,8 @@ class covarianceXsec : public covarianceBase { /// @brief Get interpolation type for a given parameter /// @param i spline parameter index, not confuse with global index inline SplineInterpolation GetParSplineInterpolation(const int i) {return SplineParams.at(i)._SplineInterpolationType;} + /// @brief Get the interpolation types for splines affecting a particular DetID + const std::vector GetSplineInterpolationFromDetID(int DetID); /// @brief Get the name of the spline associated with the spline at index i /// @param i spline parameter index, not to be confused with global index std::string GetParSplineName(const int i) {return _fSplineNames[i];} diff --git a/mcmc/MaCh3Factory.h b/mcmc/MaCh3Factory.h index 160c6f8c..b6db2378 100644 --- a/mcmc/MaCh3Factory.h +++ b/mcmc/MaCh3Factory.h @@ -124,9 +124,7 @@ std::vector MaCh3SamplePDFFactory(const std::vector& S for (size_t i = 0; i < SampleConfig.size(); ++i) { // Instantiate the sample using the specified class type - SampleType* Sample = new SampleType(SampleConfig[i], xsec); - Sample->SetXsecCov(xsec); - if (osc != nullptr) Sample->SetOscCov(osc); + SampleType* Sample = new SampleType(SampleConfig[i], xsec, osc); Sample->reweight(); // Obtain sample name and create a TString version for histogram naming diff --git a/python/samplePDF.cpp b/python/samplePDF.cpp index 1d154acb..04994d89 100644 --- a/python/samplePDF.cpp +++ b/python/samplePDF.cpp @@ -274,19 +274,6 @@ void initSamplePDF(py::module &m){ py::arg("mc_version"), py::arg("xsec_cov") ) - - .def( - "set_xsec_cov", - &samplePDFFDBase::SetXsecCov, - "Set the cross section covariance matrix object." - ) - - .def( - "set_osc_cov", - &samplePDFFDBase::SetOscCov, - "Set the oscillation parameter covariance matrix object." - ) - ; /* Not sure if this will be needed in future versions of MaCh3 so leaving commented for now diff --git a/samplePDF/FarDetectorCoreInfoStruct.h b/samplePDF/FarDetectorCoreInfoStruct.h index 10b36005..2a059799 100644 --- a/samplePDF/FarDetectorCoreInfoStruct.h +++ b/samplePDF/FarDetectorCoreInfoStruct.h @@ -10,16 +10,14 @@ struct FarDetectorCoreInfo { ~FarDetectorCoreInfo(){ delete [] isNC; } - int nutype; // 2 = numu/signue | -2 = numub | 1 = nue | -1 = nueb - int oscnutype; - int nupdg; - int nupdgUnosc; bool signal; // true if signue int nEvents; // how many MC events are there double ChannelIndex; std::string flavourName; std::vector Target; // target the interaction was on + std::vector nupdg; + std::vector nupdgUnosc; int SampleDetID; diff --git a/samplePDF/Structs.cpp b/samplePDF/Structs.cpp index ca0ac272..ee3cf688 100644 --- a/samplePDF/Structs.cpp +++ b/samplePDF/Structs.cpp @@ -1,7 +1,13 @@ #include "samplePDF/Structs.h" +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wfloat-conversion" +#include "Constants/OscillatorConstants.h" +#pragma GCC diagnostic pop + #include "TList.h" #include "TObjArray.h" +#include "manager/MaCh3Logger.h" namespace MaCh3Utils { @@ -71,6 +77,31 @@ namespace MaCh3Utils { return 0; } + int PDGToNuOscillatorFlavour(int NuPdg){ + int NuOscillatorFlavour = _BAD_INT_; + switch(std::abs(NuPdg)){ + case NuPDG::kNue: + NuOscillatorFlavour = NuOscillator::kElectron; + break; + case NuPDG::kNumu: + NuOscillatorFlavour = NuOscillator::kMuon; + break; + case NuPDG::kNutau: + NuOscillatorFlavour = NuOscillator::kTau; + break; + default: + MACH3LOG_ERROR("Unknown Nuetrino PDG {}, cannot convert to NuOscillator type", NuPdg); + break; + } + + //This is very cheeky but if the PDG is negative then multiply the PDG by -1 + // This is consistent with the treatment that NuOscillator expects as enums only + // exist for the generic matter flavour and not the anti-matter version + if(NuPdg < 0){NuOscillatorFlavour *= -1;} + + return NuOscillatorFlavour; + } + //DB Anything added here must be of the form 2^X, where X is an integer // diff --git a/samplePDF/Structs.h b/samplePDF/Structs.h index 9fd181c9..e200927e 100644 --- a/samplePDF/Structs.h +++ b/samplePDF/Structs.h @@ -270,6 +270,13 @@ namespace MaCh3Utils { /// @brief Return mass for given PDG double GetMassFromPDG(int PDG); // *************************** + + // *************************** + /// @brief Convert from PDG flavour to NuOscillator type + /// beware that in the case of anti-neutrinos the NuOscillator + /// type simply gets multiplied by -1 + int PDGToNuOscillatorFlavour(int PDG); + // *************************** extern std::unordered_map KnownDetIDsMap; extern int nKnownDetIDs; diff --git a/samplePDF/samplePDFFDBase.cpp b/samplePDF/samplePDFFDBase.cpp index fa9653aa..dd7e4570 100644 --- a/samplePDF/samplePDFFDBase.cpp +++ b/samplePDF/samplePDFFDBase.cpp @@ -1,4 +1,5 @@ #include "samplePDFFDBase.h" +#include "samplePDF/Structs.h" #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wfloat-conversion" @@ -7,23 +8,29 @@ #pragma GCC diagnostic pop #include +#include -samplePDFFDBase::samplePDFFDBase(std::string ConfigFileName, covarianceXsec* xsec_cov) : samplePDFBase() +samplePDFFDBase::samplePDFFDBase(std::string ConfigFileName, covarianceXsec* xsec_cov, covarianceOsc* osc_cov) : samplePDFBase() { MACH3LOG_INFO("-------------------------------------------------------------------"); - MACH3LOG_INFO("Ceating SamplePDFFDBase object"); + MACH3LOG_INFO("Creating SamplePDFFDBase object"); - //ETA - safety feature so you can't pass a nullptr xsec_cov - if(xsec_cov == nullptr){ - MACH3LOG_ERROR("You've passed me a nullptr xsec covariance matrix... I need this to setup splines!"); + //ETA - safety feature so you can't pass a NULL xsec_cov + if(!xsec_cov){ + MACH3LOG_ERROR("You've passed me a nullptr to a covarianceXsec... I need this to setup splines!"); throw MaCh3Exception(__FILE__, __LINE__); } - SetXsecCov(xsec_cov); + XsecCov = xsec_cov; + + if(!osc_cov){ + MACH3LOG_WARN("You have passed a nullptr to a covarianceOsc, this means I will not calculate oscillation weights"); + } + OscCov = osc_cov; samplePDFFD_array = nullptr; samplePDFFD_data = nullptr; - SampleManager = new manager(ConfigFileName.c_str()); + SampleManager = std::unique_ptr(new manager(ConfigFileName.c_str())); } samplePDFFDBase::~samplePDFFDBase() @@ -61,7 +68,7 @@ void samplePDFFDBase::ReadSampleConfig() throw MaCh3Exception(__FILE__, __LINE__); } nSamples = SampleManager->raw()["NSubSamples"].as(); - + if (!CheckNodeExists(SampleManager->raw(), "DetID")) { MACH3LOG_ERROR("ID not defined in {}, please add this!", SampleManager->GetFileName()); throw MaCh3Exception(__FILE__, __LINE__); @@ -134,8 +141,8 @@ void samplePDFFDBase::ReadSampleConfig() mc_files.push_back(mtupleprefix+osc_channel["mtuplefile"].as()+mtuplesuffix); spline_files.push_back(splineprefix+osc_channel["splinefile"].as()+splinesuffix); sample_vecno.push_back(osc_channel["samplevecno"].as()); - sample_nutype.push_back(PDGToProbs(static_cast(osc_channel["nutype"].as()))); - sample_oscnutype.push_back(PDGToProbs(static_cast(osc_channel["oscnutype"].as()))); + sample_nupdgunosc.push_back(static_cast(osc_channel["nutype"].as())); + sample_nupdg.push_back(static_cast(osc_channel["oscnutype"].as())); sample_signal.push_back(osc_channel["signal"].as()); } @@ -184,13 +191,15 @@ void samplePDFFDBase::Initialise() { MACH3LOG_INFO("Total number of events is: {}", TotalMCEvents); MACH3LOG_INFO("Setting up NuOscillator.."); - SetupNuOscillator(); + SetupNuOscillator(); MACH3LOG_INFO("Setting up Sample Binning.."); SetupSampleBinning(); MACH3LOG_INFO("Setting up Splines.."); SetupSplines(); MACH3LOG_INFO("Setting up Normalisation Pointers.."); SetupNormParameters(); + MACH3LOG_INFO("Setting up Functional Pointers.."); + SetupFunctionalParameters(); MACH3LOG_INFO("Setting up Weight Pointers.."); SetupWeightPointers(); @@ -347,15 +356,19 @@ void samplePDFFDBase::reweight() // Reweight function - Depending on Osc Calcula //KS: Reset the histograms before reweight ResetHistograms(); - std::vector OscVec(OscCov->GetNumParams()); - for (int iPar=0;iParGetNumParams();iPar++) { - #pragma GCC diagnostic push - #pragma GCC diagnostic ignored "-Wuseless-cast" - OscVec[iPar] = M3::float_t(OscCov->getParProp(iPar)); - #pragma GCC diagnostic pop - } - for (int iSample=0;iSampleCalculateProbabilities(OscVec); + //You only need to do these things if OscCov has been initialised + //if not then you're not considering oscillations + if (OscCov) { + std::vector OscVec(OscCov->GetNumParams()); + for (int iPar=0;iParGetNumParams();iPar++) { + #pragma GCC diagnostic push + #pragma GCC diagnostic ignored "-Wuseless-cast" + OscVec[iPar] = M3::float_t(OscCov->getParProp(iPar)); + #pragma GCC diagnostic pop + } + for (int iSample=0;iSampleCalculateProbabilities(OscVec); + } } fillArray(); @@ -406,6 +419,8 @@ void samplePDFFDBase::fillArray() { continue; } + std::cout << "Event passed selection, here we go!!" << std::endl; + double splineweight = 1.0; double normweight = 1.0; double funcweight = 1.0; @@ -556,109 +571,111 @@ void samplePDFFDBase::fillArray_MP() } for (unsigned int iSample=0;iSample= MCSamples[iSample].rw_lower_xbinedge[iEvent]) { - XBinToFill = MCSamples[iSample].NomXBin[iEvent]; - } - //DB - Second, check to see if the event is outside of the binning range and skip event if it is - else if (XVar < XBinEdges[0] || XVar >= XBinEdges[nXBins]) { - continue; - } - //DB - Thirdly, check the adjacent bins first as Eb+CC+EScale shifts aren't likely to move an Erec more than 1bin width - //Shifted down one bin from the event bin at nominal - else if (XVar < MCSamples[iSample].rw_lower_xbinedge[iEvent] && XVar >= MCSamples[iSample].rw_lower_lower_xbinedge[iEvent]) { - XBinToFill = MCSamples[iSample].NomXBin[iEvent]-1; - } - //Shifted up one bin from the event bin at nominal - else if (XVar < MCSamples[iSample].rw_upper_upper_xbinedge[iEvent] && XVar >= MCSamples[iSample].rw_upper_xbinedge[iEvent]) { - XBinToFill = MCSamples[iSample].NomXBin[iEvent]+1; - } - //DB - If we end up in this loop, the event has been shifted outside of its nominal bin, but is still within the allowed binning range - else { - for(unsigned int iBin=0;iBin<(XBinEdges.size()-1);iBin++) { - if (XVar >= XBinEdges[iBin] && XVar < XBinEdges[iBin+1]) { - XBinToFill = iBin; - } - } - } - - //ETA - we can probably remove this final if check on the -1? - //Maybe we can add an overflow bin to the array and assign any events to this bin? - //Might save us an extra if call? - //DB Fill relevant part of thread array - if (XBinToFill != -1 && YBinToFill != -1) { + //DB Catch negative norm weights and skip any event with a negative event. Previously we would set weight to zero and continue but that is inefficient + if (normweight <= 0.){ + MCSamples[iSample].xsec_w[iEvent] = 0.; + continue; + } + + funcweight = CalcXsecWeightFunc(iSample,iEvent); + //DB Catch negative func weights and skip any event with a negative event. Previously we would set weight to zero and continue but that is inefficient + if (funcweight <= 0.){ + MCSamples[iSample].xsec_w[iEvent] = 0.; + continue; + } + + MCSamples[iSample].xsec_w[iEvent] = splineweight*normweight*funcweight; + + totalweight = GetEventWeight(iSample, iEvent); + + //DB Catch negative weights and skip any event with a negative event + if (totalweight <= 0.){ + //std::cout << "NEGATIVE WEIGHT!!" << std::endl; + MCSamples[iSample].xsec_w[iEvent] = 0.; + continue; + } + + //DB Switch on BinningOpt to allow different binning options to be implemented + //The alternative would be to have inheritance based on BinningOpt + double XVar = (*(MCSamples[iSample].x_var[iEvent])); + + //DB Commented out by default but if we ever want to consider shifts in theta this will be needed + //double YVar = MCSamples[iSample].rw_theta[iEvent]; + //ETA - this would actually be with (*(MCSamples[iSample].y_var[iEvent])) and done extremely + //similarly to XVar now + + //DB Find the relevant bin in the PDF for each event + int XBinToFill = -1; + int YBinToFill = MCSamples[iSample].NomYBin[iEvent]; + + //DB Check to see if momentum shift has moved bins + //DB - First, check to see if the event is still in the nominal bin + if (XVar < MCSamples[iSample].rw_upper_xbinedge[iEvent] && XVar >= MCSamples[iSample].rw_lower_xbinedge[iEvent]) { + XBinToFill = MCSamples[iSample].NomXBin[iEvent]; + } + //DB - Second, check to see if the event is outside of the binning range and skip event if it is + else if (XVar < XBinEdges[0] || XVar >= XBinEdges[nXBins]) { + std::cout << "XVAR BEYOND BIN EDGES!!" << std::endl; + continue; + } + //DB - Thirdly, check the adjacent bins first as Eb+CC+EScale shifts aren't likely to move an Erec more than 1bin width + //Shifted down one bin from the event bin at nominal + else if (XVar < MCSamples[iSample].rw_lower_xbinedge[iEvent] && XVar >= MCSamples[iSample].rw_lower_lower_xbinedge[iEvent]) { + XBinToFill = MCSamples[iSample].NomXBin[iEvent]-1; + } + //Shifted up one bin from the event bin at nominal + else if (XVar < MCSamples[iSample].rw_upper_upper_xbinedge[iEvent] && XVar >= MCSamples[iSample].rw_upper_xbinedge[iEvent]) { + XBinToFill = MCSamples[iSample].NomXBin[iEvent]+1; + } + //DB - If we end up in this loop, the event has been shifted outside of its nominal bin, but is still within the allowed binning range + else { + for(unsigned int iBin=0;iBin<(XBinEdges.size()-1);iBin++) { + if (XVar >= XBinEdges[iBin] && XVar < XBinEdges[iBin+1]) { + XBinToFill = iBin; + } + } + } + + //ETA - we can probably remove this final if check on the -1? + //Maybe we can add an overflow bin to the array and assign any events to this bin? + //Might save us an extra if call? + //DB Fill relevant part of thread array + if (XBinToFill != -1 && YBinToFill != -1) { samplePDFFD_array_private[YBinToFill][XBinToFill] += totalweight; samplePDFFD_array_private_w2[YBinToFill][XBinToFill] += totalweight*totalweight; - } + } } } @@ -732,31 +749,19 @@ double samplePDFFDBase::CalcXsecWeightNorm(const int iSample, const int iEvent) return xsecw; } -void samplePDFFDBase::SetXsecCov(covarianceXsec *xsec){ - XsecCov = xsec; - - // Get the map between the normalisation parameters index, their name, what mode they should apply to, and what target - //This information is used later in CalcXsecNormsBins to decide if a parameter applies to an event - - //DB Now get this information using the DetID from the config - xsec_norms = XsecCov->GetNormParsFromDetID(SampleDetID); - nFuncParams = XsecCov->GetNumParamsFromDetID(SampleDetID, SystType::kFunc); - funcParsNames = XsecCov->GetParsNamesFromDetID(SampleDetID, SystType::kFunc); - funcParsIndex = XsecCov->GetParsIndexFromDetID(SampleDetID, SystType::kFunc); - - return; -} - void samplePDFFDBase::SetupNormParameters(){ + + xsec_norms = XsecCov->GetNormParsFromDetID(SampleDetID); + std::cout << "FOUND " << xsec_norms.size() << " norm parameters that affect this sample" << std::endl; if(!XsecCov){ - MACH3LOG_ERROR("XsecCov is not setup!"); - throw MaCh3Exception(__FILE__ , __LINE__ ); + MACH3LOG_ERROR("XsecCov is not setup!"); + throw MaCh3Exception(__FILE__ , __LINE__ ); } // Assign xsec norm bins in MCSamples tree - for (int iSample = 0; iSample < int(MCSamples.size()); ++iSample) { - CalcXsecNormsBins(iSample); + for (unsigned int iSample = 0; iSample < MCSamples.size(); ++iSample) { + CalcXsecNormsBins(iSample); } //DB @@ -790,96 +795,96 @@ void samplePDFFDBase::CalcXsecNormsBins(int iSample){ std::vector< int > XsecBins = {}; if (XsecCov) { for (std::vector::iterator it = xsec_norms.begin(); it != xsec_norms.end(); ++it) { - // Skip oscillated NC events - // Not strictly needed, but these events don't get included in oscillated predictions, so - // no need to waste our time calculating and storing information about xsec parameters - // that will never be used. - if (fdobj->isNC[iEvent] && fdobj->signal) {continue;} //DB Abstract check on MaCh3Modes to determine which apply to neutral current - - //Now check that the target of an interaction matches with the normalisation parameters - bool TargetMatch=false; - //If no target specified then apply to all modes - if ((*it).targets.size()==0) { - TargetMatch=true; - } else { - for (unsigned iTarget=0;iTarget<(*it).targets.size();iTarget++) { - if ((*it).targets.at(iTarget)== *(fdobj->Target[iEvent])) { - TargetMatch=true; - } - } - } - if (!TargetMatch) {continue;} - - //Now check that the neutrino flavour in an interaction matches with the normalisation parameters - bool FlavourMatch=false; - //If no mode specified then apply to all modes - if ((*it).pdgs.size()==0) { - FlavourMatch=true; - } else { - for (unsigned iPDG=0;iPDG<(*it).pdgs.size();iPDG++) { - if ((*it).pdgs.at(iPDG)== fdobj->nupdg) { - FlavourMatch=true; - } - } - } - if (!FlavourMatch){continue;} - - //Now check that the unoscillated neutrino flavour in an interaction matches with the normalisation parameters - bool FlavourUnoscMatch=false; - //If no mode specified then apply to all modes - if ((*it).preoscpdgs.size()==0) { - FlavourUnoscMatch=true; - } else { - for (unsigned iPDG=0;iPDG<(*it).preoscpdgs.size();iPDG++) { - if ((*it).preoscpdgs.at(iPDG) == fdobj->nupdgUnosc) { - FlavourUnoscMatch=true; - } - } - } - if (!FlavourUnoscMatch){continue;} - - //Now check that the mode of an interaction matches with the normalisation parameters - bool ModeMatch=false; - //If no mode specified then apply to all modes - if ((*it).modes.size()==0) { - ModeMatch=true; - } else { - for (unsigned imode=0;imode<(*it).modes.size();imode++) { - if ((*it).modes.at(imode)== *(fdobj->mode[iEvent])) { - ModeMatch=true; - } - } - } - if (!ModeMatch) {continue;} - - //Now check whether the norm has kinematic bounds - //i.e. does it only apply to events in a particular kinematic region? - bool IsSelected = true; - if ((*it).hasKinBounds) { - for (unsigned int iKinematicParameter = 0 ; iKinematicParameter < (*it).KinematicVarStr.size() ; ++iKinematicParameter ) { - if (ReturnKinematicParameter((*it).KinematicVarStr[iKinematicParameter], iSample, iEvent) <= (*it).Selection[iKinematicParameter][0]) { - IsSelected = false; - continue; - } - else if (ReturnKinematicParameter((*it).KinematicVarStr[iKinematicParameter], iSample, iEvent) > (*it).Selection[iKinematicParameter][1]) { - IsSelected = false; - continue; - } - } - } - //Need to then break the event loop - if(!IsSelected){ - continue; - } - - // Now set 'index bin' for each normalisation parameter - // All normalisations are just 1 bin for 2015, so bin = index (where index is just the bin for that normalisation) - int bin = (*it).index; - - //If syst on applies to a particular detector - if ((XsecCov->GetParDetID(bin) & SampleDetID)==SampleDetID) { - XsecBins.push_back(bin); - } + // Skip oscillated NC events + // Not strictly needed, but these events don't get included in oscillated predictions, so + // no need to waste our time calculating and storing information about xsec parameters + // that will never be used. + if (fdobj->isNC[iEvent] && fdobj->signal) {continue;} //DB Abstract check on MaCh3Modes to determine which apply to neutral current + + //Now check that the target of an interaction matches with the normalisation parameters + bool TargetMatch=false; + //If no target specified then apply to all modes + if ((*it).targets.size()==0) { + TargetMatch=true; + } else { + for (unsigned iTarget=0;iTarget<(*it).targets.size();iTarget++) { + if ((*it).targets.at(iTarget)== *(fdobj->Target[iEvent])) { + TargetMatch=true; + } + } + } + if (!TargetMatch) {continue;} + + //Now check that the neutrino flavour in an interaction matches with the normalisation parameters + bool FlavourMatch=false; + //If no mode specified then apply to all modes + if ((*it).pdgs.size()==0) { + FlavourMatch=true; + } else { + for (unsigned iPDG=0;iPDG<(*it).pdgs.size();iPDG++) { + if ((*it).pdgs.at(iPDG)== (*fdobj->nupdg[iEvent])) { + FlavourMatch=true; + } + } + } + if (!FlavourMatch){continue;} + + //Now check that the unoscillated neutrino flavour in an interaction matches with the normalisation parameters + bool FlavourUnoscMatch=false; + //If no mode specified then apply to all modes + if ((*it).preoscpdgs.size()==0) { + FlavourUnoscMatch=true; + } else { + for (unsigned iPDG=0;iPDG<(*it).preoscpdgs.size();iPDG++) { + if ((*it).preoscpdgs.at(iPDG) == (*fdobj->nupdgUnosc[iEvent])) { + FlavourUnoscMatch=true; + } + } + } + if (!FlavourUnoscMatch){continue;} + + //Now check that the mode of an interaction matches with the normalisation parameters + bool ModeMatch=false; + //If no mode specified then apply to all modes + if ((*it).modes.size()==0) { + ModeMatch=true; + } else { + for (unsigned imode=0;imode<(*it).modes.size();imode++) { + if ((*it).modes.at(imode)== *(fdobj->mode[iEvent])) { + ModeMatch=true; + } + } + } + if (!ModeMatch) {continue;} + + //Now check whether the norm has kinematic bounds + //i.e. does it only apply to events in a particular kinematic region? + bool IsSelected = true; + if ((*it).hasKinBounds) { + for (unsigned int iKinematicParameter = 0 ; iKinematicParameter < (*it).KinematicVarStr.size() ; ++iKinematicParameter ) { + if (ReturnKinematicParameter((*it).KinematicVarStr[iKinematicParameter], iSample, iEvent) <= (*it).Selection[iKinematicParameter][0]) { + IsSelected = false; + continue; + } + else if (ReturnKinematicParameter((*it).KinematicVarStr[iKinematicParameter], iSample, iEvent) > (*it).Selection[iKinematicParameter][1]) { + IsSelected = false; + continue; + } + } + } + //Need to then break the event loop + if(!IsSelected){ + continue; + } + + // Now set 'index bin' for each normalisation parameter + // All normalisations are just 1 bin for 2015, so bin = index (where index is just the bin for that normalisation) + int bin = (*it).index; + + //If syst on applies to a particular detector + if ((XsecCov->GetParDetID(bin) & SampleDetID)==SampleDetID) { + XsecBins.push_back(bin); + } } // end iteration over xsec_norms } // end if (xsecCov) fdobj->xsec_norms_bins[iEvent]=XsecBins; @@ -1320,44 +1325,44 @@ void samplePDFFDBase::addData(TH2D* Data) { } void samplePDFFDBase::SetupNuOscillator() { - OscillatorFactory* OscillFactory = new OscillatorFactory(); - + OscillatorFactory* OscillFactory = new OscillatorFactory(); + NuOscProbCalcers = std::vector(int(MCSamples.size())); for (size_t iSample=0;iSampleCreateOscillator(NuOscillatorConfigFile); - - if (!NuOscProbCalcers[iSample]->EvalPointsSetInConstructor()) { - std::vector EnergyArray; - for (int iEvent=0;iEventCreateOscillator(NuOscillatorConfigFile); + + if (!NuOscProbCalcers[iSample]->EvalPointsSetInConstructor()) { + std::vector EnergyArray; + for (int iEvent=0;iEventSetEnergyArrayInCalcer(EnergyArray); + + //============================================================================ + //DB Atmospheric only part + if (MCSamples[iSample].rw_truecz.size() > 0 && int(MCSamples[iSample].rw_truecz.size()) == MCSamples[iSample].nEvents) { //Can only happen if truecz has been initialised within the experiment specific code + std::vector CosineZArray; + for (int iEvent=0;iEventSetCosineZArrayInCalcer(CosineZArray); + } + //============================================================================ } - std::sort(EnergyArray.begin(),EnergyArray.end()); - - NuOscProbCalcers[iSample]->SetEnergyArrayInCalcer(EnergyArray); - - //============================================================================ - //DB Atmospheric only part - if (MCSamples[iSample].rw_truecz.size() > 0 && int(MCSamples[iSample].rw_truecz.size()) == MCSamples[iSample].nEvents) { //Can only happen if truecz has been initialised within the experiment specific code - std::vector CosineZArray; - for (int iEvent=0;iEventSetCosineZArrayInCalcer(CosineZArray); - } - //============================================================================ - - } - - NuOscProbCalcers[iSample]->Setup(); + + NuOscProbCalcers[iSample]->Setup(); + } for (int iEvent=0;iEvent 0) { //Can only happen if truecz has been initialised within the experiment specific code - //Atmospherics - MCSamples[iSample].osc_w_pointer[iEvent] = NuOscProbCalcers[iSample]->ReturnWeightPointer(InitFlav,FinalFlav,FLOAT_T(*(MCSamples[iSample].rw_etru[iEvent])),FLOAT_T(*(MCSamples[iSample].rw_truecz[iEvent]))); - } else { - //Beam - MCSamples[iSample].osc_w_pointer[iEvent] = NuOscProbCalcers[iSample]->ReturnWeightPointer(InitFlav,FinalFlav,FLOAT_T(*(MCSamples[iSample].rw_etru[iEvent]))); - } + int InitFlav = _BAD_INT_; + int FinalFlav = _BAD_INT_; + + InitFlav = MaCh3Utils::PDGToNuOscillatorFlavour((*MCSamples[iSample].nupdgUnosc[iEvent])); + FinalFlav = MaCh3Utils::PDGToNuOscillatorFlavour((*MCSamples[iSample].nupdg[iEvent])); + + if (InitFlav == _BAD_INT_ || FinalFlav == _BAD_INT_) { + MACH3LOG_ERROR("Something has gone wrong in the mapping between MCSamples[iSample].nutype and the enum used within NuOscillator"); + MACH3LOG_ERROR("MCSamples[iSample].nupdgUnosc: {}", (*MCSamples[iSample].nupdgUnosc[iEvent])); + MACH3LOG_ERROR("InitFlav: {}", InitFlav); + MACH3LOG_ERROR("MCSamples[iSample].nupdg: {}", (*MCSamples[iSample].nupdg[iEvent])); + MACH3LOG_ERROR("FinalFlav: {}", FinalFlav); + throw MaCh3Exception(__FILE__, __LINE__); + } + + if(OscCov){ + if (MCSamples[iSample].rw_truecz.size() > 0) { //Can only happen if truecz has been initialised within the experiment specific code + //Atmospherics + MCSamples[iSample].osc_w_pointer[iEvent] = NuOscProbCalcers[iSample]->ReturnWeightPointer(InitFlav,FinalFlav,FLOAT_T(*(MCSamples[iSample].rw_etru[iEvent])),FLOAT_T(*(MCSamples[iSample].rw_truecz[iEvent]))); + } else { + //Beam + MCSamples[iSample].osc_w_pointer[iEvent] = NuOscProbCalcers[iSample]->ReturnWeightPointer(InitFlav,FinalFlav,FLOAT_T(*(MCSamples[iSample].rw_etru[iEvent]))); + } + } } // end if NC } // end loop over events }// end loop over channels @@ -1435,7 +1424,9 @@ double samplePDFFDBase::GetEventWeight(int iSample, int iEntry) { #endif for (int iParam=0;iParamnEvents = nEvents_; - fdobj->nutype = -9; - fdobj->oscnutype = -9; fdobj->signal = false; fdobj->Unity = 1.; fdobj->Unity_Int = 1.; + int nEvents = fdobj->nEvents; fdobj->x_var.resize(nEvents, &fdobj->Unity); @@ -1542,6 +1532,8 @@ void samplePDFFDBase::InitialiseSingleFDMCObject(int iSample, int nEvents_) { fdobj->xsec_norm_pointers.resize(nEvents); fdobj->xsec_norms_bins.resize(nEvents); fdobj->xsec_w.resize(nEvents, 1.0); + fdobj->nupdg.resize(nEvents); + fdobj->nupdgUnosc.resize(nEvents); fdobj->isNC = new bool[nEvents]; fdobj->nxsec_spline_pointers.resize(nEvents); fdobj->xsec_spline_pointers.resize(nEvents); @@ -1586,6 +1578,8 @@ void samplePDFFDBase::InitialiseSplineObject() { MACH3LOG_INFO("Setup Far Detector splines"); fillSplineBins(); + + SplineHandler->cleanUpMemory(); } TH1* samplePDFFDBase::get1DVarHist(std::string ProjectionVar_Str, std::vector< std::vector > SelectionVec, int WeightStyle, TAxis* Axis) { diff --git a/samplePDF/samplePDFFDBase.h b/samplePDF/samplePDFFDBase.h index b11a5881..6bbfac1f 100644 --- a/samplePDF/samplePDFFDBase.h +++ b/samplePDF/samplePDFFDBase.h @@ -1,12 +1,5 @@ #pragma once -//C++ includes -#include - -//ROOT includes -#include "THStack.h" -#include "TLegend.h" - //MaCh3 includes #include "splines/splineFDBase.h" @@ -26,7 +19,7 @@ class samplePDFFDBase : public samplePDFBase //######################################### Functions ######################################### samplePDFFDBase(){}; - samplePDFFDBase(std::string mc_version, covarianceXsec* xsec_cov); + samplePDFFDBase(std::string mc_version, covarianceXsec* xsec_cov, covarianceOsc* osc_cov = nullptr); virtual ~samplePDFFDBase(); int GetNDim(){return nDimensions;} //DB Function to differentiate 1D or 2D binning @@ -47,12 +40,6 @@ class samplePDFFDBase : public samplePDFBase void reweight(); double GetEventWeight(int iSample, int iEntry); - - - void SetXsecCov(covarianceXsec* xsec_cov); - - /// @brief setup the Oscillation covariance object to get values to calculate probailities from - void SetOscCov(covarianceOsc* osc_cov){OscCov = osc_cov;}; /// @brief including Dan's magic NuOscillator void SetupNuOscillator(); @@ -131,6 +118,8 @@ class samplePDFFDBase : public samplePDFBase //=============================================================================== /// @brief ETA - a function to setup and pass values to functional parameters where you need to pass a value to some custom reweight calc or engine + virtual void SetupFunctionalParameters(){}; + /// @brief Update the functional parameter values to the latest propsed values. Needs to be called before every new reweight so is called in fillArray virtual void PrepFunctionalParameters(){}; /// @brief ETA - generic function applying shifts virtual void applyShifts(int iSample, int iEvent){(void) iSample; (void) iEvent;}; @@ -154,10 +143,6 @@ class samplePDFFDBase : public samplePDFBase virtual const double* GetPointerToKinematicParameter(std::string KinematicParamter, int iSample, int iEvent) = 0; virtual const double* GetPointerToKinematicParameter(double KinematicVariable, int iSample, int iEvent) = 0; - // Function to setup Functional and shift parameters. This isn't idea but - // do this in your experiment specific code for now as we don't have a - // generic treatment for func and shift pars yet - virtual void SetupFuncParameters(){return;}; void SetupNormParameters(); //=============================================================================== @@ -207,15 +192,15 @@ class samplePDFFDBase : public samplePDFBase //DB Covariance Objects //ETA - All experiments will need an xsec, det and osc cov //these should be added to samplePDFBase to be honest - covarianceXsec *XsecCov; - covarianceOsc *OscCov; + covarianceXsec *XsecCov = nullptr; + covarianceOsc *OscCov = nullptr; //=============================================================================== /// @brief Keep track of the dimensions of the sample binning - int nDimensions; + int nDimensions = _BAD_INT_; /// @brief A unique ID for each sample based on powers of two for quick binary operator comparisons - int SampleDetID; + int SampleDetID = _BAD_INT_; /// holds "TrueNeutrinoEnergy" and the strings used for the sample binning. std::vector SplineBinnedVars; @@ -224,10 +209,6 @@ class samplePDFFDBase : public samplePDFBase /// @brief Information to store for normalisation pars std::vector xsec_norms; - /// @brief the total number of function parameters found in the xsec model - int nFuncParams; - std::vector funcParsNames; - std::vector funcParsIndex; //=========================================================================== //DB Vectors to store which kinematic cuts we apply @@ -235,7 +216,7 @@ class samplePDFFDBase : public samplePDFBase //What gets used in IsEventSelected, which gets set equal to user input plus //all the vectors in StoreSelection /// @brief the Number of selections in the - int NSelections; + int NSelections = _BAD_INT_; /// @brief What gets pulled from config options, these are constant after loading in /// this is of length 3: 0th index is the value, 1st is lower bound, 2nd is upper bound @@ -249,7 +230,7 @@ class samplePDFFDBase : public samplePDFBase std::vector< std::vector > Selection; //=========================================================================== - manager* SampleManager; + std::unique_ptr SampleManager; void InitialiseSingleFDMCObject(int iSample, int nEvents); void InitialiseSplineObject(); @@ -262,7 +243,7 @@ class samplePDFFDBase : public samplePDFBase std::vector mc_files; std::vector spline_files; std::vector sample_vecno; - std::vector sample_oscnutype; - std::vector sample_nutype; + std::vector sample_nupdg; + std::vector sample_nupdgunosc; std::vector sample_signal; }; diff --git a/splines/splineFDBase.cpp b/splines/splineFDBase.cpp index 6f93c575..d6d3eb75 100644 --- a/splines/splineFDBase.cpp +++ b/splines/splineFDBase.cpp @@ -1,4 +1,6 @@ #include "splineFDBase.h" +#include +#include "samplePDF/Structs.h" #pragma GCC diagnostic ignored "-Wuseless-cast" #pragma GCC diagnostic ignored "-Wfloat-conversion" @@ -8,7 +10,7 @@ splineFDBase::splineFDBase(covarianceXsec *xsec_) : SplineBase() { //**************************************** if (!xsec_) { - MACH3LOG_ERROR("Trying to create splineSKBase with null covariance object"); + MACH3LOG_ERROR("Trying to create splineFDBase with uninitialised covariance object"); throw MaCh3Exception(__FILE__, __LINE__); } xsec = xsec_; @@ -42,6 +44,10 @@ void splineFDBase::cleanUpMemory() { GlobalSystIndex.shrink_to_fit(); UniqueSystNames.clear(); UniqueSystNames.shrink_to_fit(); + //Really make sure all the memory is cleared + for(auto Spline : splinevec_Monolith){ + if(Spline){delete Spline;} + } splinevec_Monolith.clear(); splinevec_Monolith.shrink_to_fit(); if(isflatarray != nullptr) delete isflatarray; @@ -65,6 +71,10 @@ bool splineFDBase::AddSample(std::string SampleName, int DetID, std::vector SplineInterpolation_Sample = xsec->GetSplineInterpolationFromDetID(DetID); + // Keep track of this for all samples + SplineInterpolationTypes.push_back(SplineInterpolation_Sample); + //std::vector SplineParsIndex_Sample_temp = xsec->GetSplineParsIndexFromDetID(DetID); std::vector SplineFileParPrefixNames_Sample = xsec->GetSplineParsNamesFromDetID(DetID); @@ -74,6 +84,11 @@ bool splineFDBase::AddSample(std::string SampleName, int DetID, std::vector> SplineModeVecs_Sample = StripDuplicatedModes(xsec->GetSplineModeVecFromDetID(DetID)); MACH3LOG_INFO("SplineModeVecs_Sample is of size {}", SplineModeVecs_Sample.size()); SplineModeVecs.push_back(SplineModeVecs_Sample); + // for(auto SplineMode : SplineModeVecs_Sample){ + // std::cout << SplineMode << ","; + // } + std::cout << "/n" << std::endl; + MACH3LOG_INFO("SplineModeVecs is of size {}", SplineModeVecs.size()); int nOscChan = int(OscChanFileNames.size()); @@ -163,7 +178,7 @@ void splineFDBase::TransferToMonolith() } int splineKnots; - if(splinevec_Monolith[splineindex]!=NULL){ + if(splinevec_Monolith[splineindex]){ isflatarray[splineindex]=false; splineKnots=splinevec_Monolith[splineindex]->GetNp(); @@ -182,8 +197,8 @@ void splineFDBase::TransferToMonolith() manycoeff_arr[(iCoeff+i)*4+j]=tmpManyCoeffArr[i*4+j]; } } - delete tmpXCoeffArr; - delete tmpManyCoeffArr; + delete[] tmpXCoeffArr; + delete[] tmpManyCoeffArr; } else { isflatarray[splineindex]=true; @@ -227,8 +242,7 @@ void splineFDBase::FindSplineSegment() #ifdef MULTITHREAD #pragma omp parallel //for schedule(dynamic) #endif - for (int iSyst = 0; iSyst < nUniqueSysts; iSyst++) - { + for (int iSyst = 0; iSyst < nUniqueSysts; iSyst++) { int nPoints = UniqueSystNKnots[iSyst]; std::vector xArray = UniqueSystXPts[iSyst]; @@ -238,23 +252,24 @@ void splineFDBase::FindSplineSegment() M3::float_t xvar = M3::float_t(xsec->getParProp(GlobalIndex)); xVarArray[iSyst]=xvar; - + M3::int_t segment = 0; - M3::int_t kHigh = M3::int_t(nPoints - 1); + M3::int_t kHigh = M3::int_t(nPoints - 1); //KS: We expect new segment is very close to previous const M3::int_t PreviousSegment = M3::int_t(UniqueSystCurrSegment[iSyst]); //KS: It is quite probable the new segment is same as in previous step so try to avoid binary search - if( xArray[PreviousSegment+1] > xvar && xvar >= xArray[PreviousSegment] ){segment = PreviousSegment;} + if( xArray[PreviousSegment+1] > xvar && xvar >= xArray[PreviousSegment] ) { + segment = PreviousSegment; + } else if (xvar <= xArray[0]) { // If the variation is below the lowest saved spline point - else if (xvar <= xArray[0]) { - segment = 0; - // If the variation is above the highest saved spline point - } else if (xvar >= xArray[nPoints-1]) { - //CW: Yes, the -2 is indeed correct, see TSpline.cxx:814 and //see: https://savannah.cern.ch/bugs/?71651 - segment = kHigh; - //KS: It is quite probable the new segment is same as in previous step so try to avoid binary search - } else { + segment = 0; + // If the variation is above the highest saved spline point + } else if (xvar >= xArray[nPoints-1]) { + //CW: Yes, the -2 is indeed correct, see TSpline.cxx:814 and //see: https://savannah.cern.ch/bugs/?71651 + segment = kHigh; + //KS: It is quite probable the new segment is same as in previous step so try to avoid binary search + } else { // The top point we've got M3::int_t kHalf = 0; // While there is still a difference in the points (we haven't yet found the segment) @@ -274,23 +289,23 @@ void splineFDBase::FindSplineSegment() if (segment >= nPoints-1 && nPoints > 1){segment = M3::int_t(nPoints-2);} UniqueSystCurrSegment[iSyst] = segment; - -//#ifdef DEBUG -// if (SplineInfoArray[i].xPts[segment] > xvar && segment != 0) { -// std::cerr << "Found a segment which is _ABOVE_ the variation!" << std::endl; -// std::cerr << "IT SHOULD ALWAYS BE BELOW! (except when segment 0)" << std::endl; -// std::cerr << "Spline: "<< i << std::endl; -// -// std::cerr << "Found segment = " << segment << std::endl; -// std::cerr << "Doing variation = " << xvar << std::endl; -// std::cerr << "x in spline = " << SplineInfoArray[i].xPts[segment] << std::endl; -// for (_M3::int_t_ j = 0; j < SplineInfoArray[j].nPts; ++j) { -// std::cerr << " " << j << " = " << SplineInfoArray[i].xPts[j] << std::endl; -// } -// std::cerr << __FILE__ << ":" << __LINE__ << std::endl; -// throw; -// } -//#endif + + //#ifdef DEBUG + // if (SplineInfoArray[i].xPts[segment] > xvar && segment != 0) { + // std::cerr << "Found a segment which is _ABOVE_ the variation!" << std::endl; + // std::cerr << "IT SHOULD ALWAYS BE BELOW! (except when segment 0)" << std::endl; + // std::cerr << "Spline: "<< i << std::endl; + // + // std::cerr << "Found segment = " << segment << std::endl; + // std::cerr << "Doing variation = " << xvar << std::endl; + // std::cerr << "x in spline = " << SplineInfoArray[i].xPts[segment] << std::endl; + // for (_M3::int_t_ j = 0; j < SplineInfoArray[j].nPts; ++j) { + // std::cerr << " " << j << " = " << SplineInfoArray[i].xPts[j] << std::endl; + // } + // std::cerr << __FILE__ << ":" << __LINE__ << std::endl; + // throw; + // } + //#endif } //end loop over params } @@ -394,6 +409,12 @@ std::vector splineFDBase::FindSplineBinning(std::string FileName, std:: std::vector ReturnVec; int iSample=getSampleIndex(SampleName); + //Try declaring these outside of TFile so they aren't owned by File + int nDummyBins = 1; + double DummyEdges[2]; + DummyEdges[0] = -1e15; + DummyEdges[1] = 1e15; + TAxis* DummyAxis = new TAxis(nDummyBins, DummyEdges); TH2F* Hist2D = nullptr; TH3F* Hist3D = nullptr; @@ -406,16 +427,18 @@ std::vector splineFDBase::FindSplineBinning(std::string FileName, std:: } MACH3LOG_INFO("Finding binning for:"); - spdlog::info("{}", FileName); + MACH3LOG_INFO("{}", FileName); bool isHist2D = false; bool isHist3D = false; - TObject *Obj = File->Get("dev_tmp_0_0"); + std::string TemplateName = "dev_tmp_0_0"; + TObject *Obj = File->Get(TemplateName.c_str()); //If you can't find dev_tmp_0_0 then this will cause a problem if (!Obj) { - Obj = File->Get("dev_tmp.0.0"); + TemplateName = "dev_tmp.0.0"; + Obj = File->Get(TemplateName.c_str()); if (!Obj) { MACH3LOG_ERROR("Error: could not find dev_tmp_0_0 in spline file. Spline binning cannot be set!"); @@ -447,9 +470,10 @@ std::vector splineFDBase::FindSplineBinning(std::string FileName, std:: if (Dimensions[iSample] != 2) { MACH3LOG_ERROR("Trying to load a 2D spline template when nDim={}", Dimensions[iSample]); - throw MaCh3Exception(__FILE__, __LINE__); } - Hist2D = File->Get("dev_tmp_0_0"); + throw MaCh3Exception(__FILE__, __LINE__); + //Hist2D = std::unique_ptr(File->Get("dev_tmp_0_0")); + Hist2D = File->Get(TemplateName.c_str()); } if (isHist3D) @@ -461,47 +485,36 @@ std::vector splineFDBase::FindSplineBinning(std::string FileName, std:: MACH3LOG_ERROR("Trying to load a 3D spline template when nDim={}", Dimensions[iSample]); throw MaCh3Exception(__FILE__ , __LINE__ ); } - } - - double DummyEdges[2] = {-1e-15, 1e15}; - auto DummyAxis = std::unique_ptr(new TAxis(1, DummyEdges)); - - if (Dimensions[iSample] == 2) - { - if(isHist2D){ - ReturnVec.push_back(static_cast(Hist2D->GetXaxis()->Clone())); - ReturnVec.push_back(static_cast(Hist2D->GetYaxis()->Clone())); - ReturnVec.push_back(static_cast(DummyAxis->Clone())); - } - else if(isHist3D){ - ReturnVec.push_back(static_cast(Hist3D->GetXaxis()->Clone())); - ReturnVec.push_back(static_cast(Hist3D->GetYaxis()->Clone())); - ReturnVec.push_back(static_cast(DummyAxis->Clone())); - } - } - else if (Dimensions[iSample] == 3) - { - ReturnVec.push_back(static_cast(Hist3D->GetXaxis()->Clone())); - ReturnVec.push_back(static_cast(Hist3D->GetYaxis()->Clone())); - ReturnVec.push_back(static_cast(Hist3D->GetZaxis()->Clone())); - } - else - { + Hist3D = File->Get(TemplateName.c_str()); + } + + if (Dimensions[iSample] == 2) { + if(isHist2D){ + ReturnVec.push_back(static_cast(Hist2D->GetXaxis()->Clone())); + ReturnVec.push_back(static_cast(Hist2D->GetYaxis()->Clone())); + ReturnVec.push_back(static_cast(DummyAxis->Clone())); + } else if (isHist3D) { + ReturnVec.push_back(static_cast(Hist3D->GetXaxis()->Clone())); + ReturnVec.push_back(static_cast(Hist3D->GetYaxis()->Clone())); + ReturnVec.push_back(static_cast(DummyAxis->Clone())); + } + } else if (Dimensions[iSample] == 3) { + ReturnVec.push_back(static_cast(Hist3D->GetXaxis()->Clone())); + ReturnVec.push_back(static_cast(Hist3D->GetYaxis()->Clone())); + ReturnVec.push_back(static_cast(Hist3D->GetZaxis()->Clone())); + } else { MACH3LOG_ERROR("Number of dimensions not valid! Given: {}", Dimensions[iSample]); throw MaCh3Exception(__FILE__ , __LINE__ ); } - for (unsigned int iAxis = 0; iAxis < ReturnVec.size(); ++iAxis) - { + for (unsigned int iAxis = 0; iAxis < ReturnVec.size(); ++iAxis) { PrintBinning(ReturnVec[iAxis]); } MACH3LOG_INFO("Left PrintBinning now tidying up"); - - //Make sure these actually get deleted - delete Hist2D; - delete Hist3D; + delete DummyAxis; File->Close(); + return ReturnVec; } @@ -534,7 +547,7 @@ int splineFDBase::CountNumberOfLoadedSplines(bool NonFlat, int Verbosity) if (isValidSplineIndex(SampleName, iOscChan, iSyst, iMode, iVar1, iVar2, iVar3)) { int splineindex = indexvec[iSample][iOscChan][iSyst][iMode][iVar1][iVar2][iVar3]; - if (splinevec_Monolith[splineindex] != NULL) + if (splinevec_Monolith[splineindex]) { SampleCounter_NonFlat += 1; } @@ -613,7 +626,7 @@ void splineFDBase::PrepForReweight() for (unsigned int iVar3 = 0; iVar3 < indexvec[iSample][iOscChan][iSyst][iMode][iVar1][iVar2].size(); iVar3++) { // Loop over third dimension int splineindex=indexvec[iSample][iOscChan][iSyst][iMode][iVar1][iVar2][iVar3]; - if (splinevec_Monolith[splineindex] != NULL) + if (splinevec_Monolith[splineindex]) { UniqueSystSplines.push_back(splinevec_Monolith[splineindex]); UniqueSystIndices.push_back(GlobalSystIndex[iSample][iSyst]); @@ -700,7 +713,7 @@ void splineFDBase::PrepForReweight() //Isn't this just doing what CountNumberOfLoadedSplines() does? int nCombinations_FlatSplines = 0; int nCombinations_All = 0; - // DB Now actually loop over splines to determine which are all NULL + // DB Now actually loop over splines to determine which are all null i.e. flat for (unsigned int iSample = 0; iSample < indexvec.size(); iSample++) { // Loop over systematics for (unsigned int iOscChan = 0; iOscChan < indexvec[iSample].size(); iOscChan++) @@ -717,7 +730,7 @@ void splineFDBase::PrepForReweight() for (unsigned int iVar3 = 0; iVar3 < indexvec[iSample][iOscChan][iSyst][iMode][iVar1][iVar2].size(); iVar3++) { // Loop over third dimension int splineindex=indexvec[iSample][iOscChan][iSyst][iMode][iVar1][iVar2][iVar3]; - if (splinevec_Monolith[splineindex] != NULL) + if (splinevec_Monolith[splineindex]) { nCombinations_All += 1; } else{ diff --git a/splines/splineFDBase.h b/splines/splineFDBase.h index a49a5ca0..fcecaa36 100644 --- a/splines/splineFDBase.h +++ b/splines/splineFDBase.h @@ -4,6 +4,7 @@ #include "TH3F.h" //MaCh3 +#include "samplePDF/Structs.h" #include "splines/SplineBase.h" /// @brief Bin-by-bin class calculating response for spline parameters. @@ -91,6 +92,9 @@ class splineFDBase : public SplineBase { /// to evaluate splines at. Each internal vector will be of size of the number of spline /// systematics which affect that sample. std::vector< std::vector > GlobalSystIndex; + /// @brief spline interpolation types for each sample. These vectors are from + /// a call to GetSplineInterpolationFromDetID() + std::vector< std::vector > SplineInterpolationTypes; int nUniqueSysts; std::vector UniqueSystNames;