diff --git a/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx b/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx index 8b11a3f..014f927 100644 --- a/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx +++ b/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx @@ -1,913 +1,924 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #include "Smear_SVDUnfold_Propagation_Osc.h" #include "SmearceptanceUtils.h" #include "OscWeightEngine.h" #include "HistogramInputHandler.h" void Smear_SVDUnfold_Propagation_Osc::AddNDInputs(nuiskey &samplekey) { NDSample nds; // Plot Setup ------------------------------------------------------- // Check that we have the relevant near detector histograms specified. if (!NDSamples.size()) { // If this is the first ND sample, take from the // sample input InputHandlerBase *InputBase = GetInput(); if (InputBase->GetType() != kHISTO) { THROW( "Smear_SVDUnfold_Propagation_Osc expects a Histogram input that " "contains the ND observed spectrum."); } HistoInputHandler *HInput = dynamic_cast(InputBase); if (!HInput) { THROW( "Smear_SVDUnfold_Propagation_Osc expects a Histogram input that " "contains the ND observed spectrum."); } if (HInput->NHistograms() != 2) { THROW( "Input expected to contain 2 histograms. " "HISTO:input.root[NDObs_TH1D,NDSmear_TH2D]"); } nds.NDDataHist = dynamic_cast(HInput->GetHistogram(0)); nds.NDToSpectrumSmearingMatrix = dynamic_cast(HInput->GetHistogram(1)); if (!nds.NDDataHist) { THROW("Expected a valid TH1D input for the ND observed spectrum."); } if (!nds.NDToSpectrumSmearingMatrix) { THROW("Expected a valid TH2D input for the ND observed smearing."); } } else { std::vector NDObsInputs = PlotUtils::GetTH1sFromRootFile(samplekey.GetS("ObsInput")); if (NDObsInputs.size() < 2) { THROW( "Near detector sample must contain the observed ERec spectrum and " "the " "ND ETrue/ERec smearing matrix. e.g. " "ObsInput=\"input.root[NDObs_species,NDSmearing_species]\""); } nds.NDDataHist = dynamic_cast(NDObsInputs[0]); if (!nds.NDDataHist) { ERROR(FTL, "First histogram from ObsInput attribute was not a TH1D containing " "the near detector observed ERec spectrum (" << samplekey.GetS("ObsInput") << ")."); THROW( "Near detector sample must contain the observed ERec spectrum and " "the " "ND ETrue/ERec smearing matrix. e.g. " "ObsInput=\"input.root[FDObs_species,FDSmearing_species]\""); } nds.NDToSpectrumSmearingMatrix = dynamic_cast(NDObsInputs[1]); if (!nds.NDToSpectrumSmearingMatrix) { ERROR( FTL, "Second histogram from ObsInput attribute was not a TH2D containing " "the near detector ETrue/ERec smearing matrix (" << samplekey.GetS("ObsInput") << ")."); THROW( "Near detector sample must contain the observed ERec spectrum and " "the " "ND ETrue/ERec smearing matrix. e.g. " "ObsInput=\"input.root[FDObs_species,FDSmearing_species]\""); } } nds.NDDataHist->Scale(ScalePOT); if (UseRateErrors) { for (Int_t bi_it = 1; bi_it < nds.NDDataHist->GetXaxis()->GetNbins() + 1; ++bi_it) { nds.NDDataHist->SetBinError(bi_it, sqrt(nds.NDDataHist->GetBinContent(bi_it))); } } nds.TruncateStart = 0; if (samplekey.Has("TruncateStart")) { nds.TruncateStart = samplekey.GetI("TruncateStart"); } nds.TruncateUpTo = 0; if (samplekey.Has("TruncateUpTo")) { nds.TruncateUpTo = samplekey.GetI("TruncateUpTo"); QLOG(SAM, "\tAllowed to truncate unfolding matrix by up to " << nds.TruncateUpTo << " singular values to limit negative ENu spectra."); } nds.NuPDG = 14; if (samplekey.Has("NuPDG")) { nds.NuPDG = samplekey.GetI("NuPDG"); } NDSamples.push_back(nds); } void Smear_SVDUnfold_Propagation_Osc::SetupNDInputs() { for (size_t nd_it = 0; nd_it < NDSamples.size(); ++nd_it) { NDSample &nds = NDSamples[nd_it]; TMatrixD NDToSpectrumResponseMatrix_l = SmearceptanceUtils::GetMatrix( SmearceptanceUtils::SVDGetInverse(nds.NDToSpectrumSmearingMatrix)); nds.NDToSpectrumResponseMatrix.ResizeTo(NDToSpectrumResponseMatrix_l); nds.NDToSpectrumResponseMatrix = NDToSpectrumResponseMatrix_l; if (nds.TruncateStart != 0) { TMatrixD NDToSpectrumResponseMatrix_l = SmearceptanceUtils::GetMatrix(SmearceptanceUtils::SVDGetInverse( nds.NDToSpectrumSmearingMatrix, nds.TruncateStart)); nds.NDToSpectrumResponseMatrix.ResizeTo(NDToSpectrumResponseMatrix_l); nds.NDToSpectrumResponseMatrix = NDToSpectrumResponseMatrix_l; } if (nds.TruncateStart >= nds.TruncateUpTo) { nds.TruncateUpTo = nds.TruncateStart + 1; } UnfoldToNDETrueSpectrum(nd_it); } } void Smear_SVDUnfold_Propagation_Osc::ReadExtraConfig(nuiskey &samplekey) { UseRateErrors = false; if (samplekey.Has("SetErrorsFromRate")) { UseRateErrors = samplekey.GetI("SetErrorsFromRate"); } NDetectorInfo.first = 0xdeadbeef; if (samplekey.Has("DetectorVolume") && samplekey.Has("DetectorDensity")) { NDetectorInfo.first = samplekey.GetD("DetectorVolume"); NDetectorInfo.second = samplekey.GetD("DetectorDensity"); double TargetMass_kg = NDetectorInfo.first * NDetectorInfo.second; QLOG(SAM, "\tND sample detector mass : "); QLOG(SAM, "\t\tTarget volume : " << NDetectorInfo.first); QLOG(SAM, "\t\tTarget density : " << NDetectorInfo.second); QLOG(SAM, "\t\tTarget mass : " << TargetMass_kg << " kg"); } ScalePOT = 1; if (samplekey.Has("ScalePOT")) { ScalePOT = samplekey.GetD("ScalePOT"); } } void Smear_SVDUnfold_Propagation_Osc::AddFDTarget(nuiskey &nk) { FDSample fds; fds.FitRegion_Min = 0xdeadbeef; if (nk.Has("FitRegion_Min")) { fds.FitRegion_Min = nk.GetD("FitRegion_Min"); QLOG(SAM, "FD Sample [" << FDSamples.size() << "] imposes FitRegion E_nu > " << fds.FitRegion_Min); } if ((FitRegion_Min == 0xdeadbeef) || FitRegion_Min > fds.FitRegion_Min) { FitRegion_Min = fds.FitRegion_Min; } nk.Print(); fds.FitRegion_Max = 0xdeadbeef; if (nk.Has("FitRegion_Max")) { fds.FitRegion_Max = nk.GetD("FitRegion_Max"); QLOG(SAM, "FD Sample [" << FDSamples.size() << "] imposes FitRegion E_nu < " << fds.FitRegion_Max); } if ((FitRegion_Max == 0xdeadbeef) || FitRegion_Max < fds.FitRegion_Max) { FitRegion_Max = fds.FitRegion_Max; } fds.OscillateToPDG = 0; if (nk.Has("OscillateToPDG")) { fds.OscillateToPDG = nk.GetD("OscillateToPDG"); } std::vector FDNDRatioElements = nk.GetListOfChildNodes("FDNDRatio"); for (size_t fdnd_it = 0; fdnd_it < FDNDRatioElements.size(); ++fdnd_it) { nuiskey &fnr = FDNDRatioElements[fdnd_it]; if (fnr.Has("FromPDG") && fnr.Has("DivergenceFactor")) { fds.FDNDRatios[fnr.GetI("FromPDG")] = fnr.GetD("DivergenceFactor"); QLOG(SAM, "FDND DivergenceFactor for far detector sample index: " << FDSamples.size() << " for PDG: " << fnr.GetI("FromPDG") << " -> " << fds.OscillateToPDG << " = " << fnr.GetD("DivergenceFactor")); } else { THROW( "Far detector sample contained FDNDRatio element, but couldn't find " "both FromPDG and Factor attributes."); } } fds.FDNDMassRatio = 1; if (NDetectorInfo.first != 0xdeadbeef) { if ((!nk.Has("DetectorVolume")) || (!nk.Has("DetectorDensity"))) { THROW( "Near detector sample has specified volume but FD doesn't. This is " "needed to scale the predicted event rate by the mass ratio."); } fds.DetectorInfo.first = nk.GetD("DetectorVolume"); fds.DetectorInfo.second = nk.GetD("DetectorDensity"); double TargetMass_kg = fds.DetectorInfo.first * fds.DetectorInfo.second; fds.FDNDMassRatio = TargetMass_kg / (NDetectorInfo.first * NDetectorInfo.second); QLOG(SAM, "\tFD[" << FDSamples.size() << "] Event rate prediction : "); QLOG(SAM, "\t\tTarget volume : " << fds.DetectorInfo.first); QLOG(SAM, "\t\tTarget density : " << fds.DetectorInfo.second); QLOG(SAM, "\t\tFD/ND mass : " << fds.FDNDMassRatio); } if (!nk.Has("ObsInput")) { THROW("Far detector sample must specify at least ObsInput."); } std::vector FDObsInputs = PlotUtils::GetTH1sFromRootFile(nk.GetS("ObsInput")); if (FDObsInputs.size() < 2) { THROW( "Far detector sample must contain the observed ERec spectrum and the " "FD ETrue/ERec smearing matrix. " "ObsInput=\"input.root[FDObs_species,FDSmearing_species]\""); } fds.FDDataHist = NULL; for (size_t hist_it = 0; hist_it < FDObsInputs.size() - 1; ++hist_it) { if (!dynamic_cast(FDObsInputs[hist_it])) { ERROR(FTL, "Input spectrum index " << hist_it << " from ObsInput attribute was not a TH1D containing " "a far detector observed ERec spectrum (" << nk.GetS("ObsInput") << ")."); THROW( "Far detector sample must contain the observed ERec spectrum and the " "FD ETrue/ERec smearing matrix. " "ObsInput=\"input.root[FDObs_species,(FDObs_species2),FDSmearing_" "species]\""); } FDObsInputs[hist_it]->Scale(ScalePOT); if (!fds.FDDataHist) { fds.FDDataHist = dynamic_cast(FDObsInputs[hist_it]); } else { fds.FDDataHist->Add(dynamic_cast(FDObsInputs[hist_it])); } QLOG(SAM, "Added " << (FDObsInputs.size() - 1) << " far detector component spectra to form Observed " "spectra for sample index " << FDSamples.size() << "."); } fds.SpectrumToFDSmearingMatrix_TH2 = dynamic_cast(FDObsInputs.back()); if (!fds.SpectrumToFDSmearingMatrix_TH2) { ERROR(FTL, "last histogram from ObsInput attribute was not a TH2D containing " "the far detector ETrue/ERec smearing matrix (" << nk.GetS("ObsInput") << ")."); THROW( "Far detector sample must contain the observed ERec spectrum and the " "FD ETrue/ERec smearing matrix. " "ObsInput=\"input.root[FDObs_species,FDSmearing_species]\""); } TMatrixD SpectrumToFDSmearingMatrix_l = SmearceptanceUtils::GetMatrix(fds.SpectrumToFDSmearingMatrix_TH2); fds.SpectrumToFDSmearingMatrix.ResizeTo(SpectrumToFDSmearingMatrix_l); fds.SpectrumToFDSmearingMatrix = SpectrumToFDSmearingMatrix_l; FDSamples.push_back(fds); } void Smear_SVDUnfold_Propagation_Osc::FinaliseFDSamples() { std::stringstream ss(""); for (size_t fds_it = 0; fds_it < FDSamples.size(); ++fds_it) { FDSample &fds = FDSamples[fds_it]; // Set up FD histograms. // ============================== for (size_t nd_it = 0; nd_it < NDSamples.size(); ++nd_it) { NDSample &nds = NDSamples[nd_it]; TH1D *sampleHist = static_cast(nds.ND_Unfolded_Spectrum_Hist->Clone()); sampleHist->Reset(); ss.str(""); ss << "FD_Propagated_Spectrum_Hist_" << fds_it << "_NDSample_" << nd_it; sampleHist->SetName(ss.str().c_str()); fds.FD_Propagated_Spectrum_Hist_NDSamples.push_back(sampleHist); } fds.FD_Propagated_Spectrum_Hist = static_cast( fds.FD_Propagated_Spectrum_Hist_NDSamples[0]->Clone()); fds.FD_Propagated_Spectrum_Hist->Reset(); ss.str(""); ss << "FD_Propagated_Spectrum_Hist_" << fds_it; fds.FD_Propagated_Spectrum_Hist->SetName(ss.str().c_str()); for (size_t nd_it = 0; nd_it < NDSamples.size(); ++nd_it) { NDSample &nds = NDSamples[nd_it]; TH1D *sampleHist = static_cast(fds.FD_Propagated_Spectrum_Hist->Clone()); sampleHist->Reset(); ss.str(""); ss << "NDFD_Corrected_Spectrum_Hist_" << fds_it << "_NDSample_" << nd_it; sampleHist->SetName(ss.str().c_str()); fds.NDFD_Corrected_Spectrum_Hist_NDSamples.push_back(sampleHist); } fds.NDFD_Corrected_Spectrum_Hist = static_cast(fds.FD_Propagated_Spectrum_Hist->Clone()); fds.NDFD_Corrected_Spectrum_Hist->Reset(); ss.str(""); ss << "NDFD_Corrected_Spectrum_Hist_" << fds_it; fds.NDFD_Corrected_Spectrum_Hist->SetName(ss.str().c_str()); fds.FD_Smeared_Spectrum_Hist = new TH1D(*fds.SpectrumToFDSmearingMatrix_TH2->ProjectionX()); fds.FD_Smeared_Spectrum_Hist->SetDirectory(NULL); fds.FD_Smeared_Spectrum_Hist->Reset(); ss.str(""); ss << "FD_Smeared_Spectrum_Hist_" << fds_it; fds.FD_Smeared_Spectrum_Hist->SetName(ss.str().c_str()); for (size_t nd_it = 0; nd_it < NDSamples.size(); ++nd_it) { NDSample &nds = NDSamples[nd_it]; if (!fds.FDNDRatios.count(nds.NuPDG)) { ERROR(WRN, "Have an ND sample that provides PDG:" << nds.NuPDG << " neutrinos but far detector sample index " << fds_it << " doesn't have an NDFD ratio for this sample. " "Setting to 0 (contribution from sample ignored.)"); fds.FDNDRatios[nds.NuPDG] = 0; } } } } Int_t Smear_SVDUnfold_Propagation_Osc::GetFDSampleNAnalysisBins(size_t fds_it) { if (fds_it >= FDSamples.size()) { THROW("Requested FD sample index " << fds_it << " but only initialised " << FDSamples.size()); } FDSample &fds = FDSamples[fds_it]; Int_t NAnalysisBins = 0; for (Int_t bi_it = 1; bi_it < fds.FDDataHist->GetXaxis()->GetNbins() + 1; ++bi_it) { if ((fds.FitRegion_Min != 0xdeadbeef) && (fds.FDDataHist->GetXaxis()->GetBinUpEdge(bi_it) <= fds.FitRegion_Min)) { continue; } if ((fds.FitRegion_Max != 0xdeadbeef) && (fds.FDDataHist->GetXaxis()->GetBinLowEdge(bi_it) > fds.FitRegion_Max)) { continue; } NAnalysisBins++; } return NAnalysisBins; } void Smear_SVDUnfold_Propagation_Osc::SetupChi2Hists() { fDataHist = new TH1D("SmearSVDUnfold", "", NFDAnalysisBins, 0, NFDAnalysisBins); fDataHist->SetNameTitle((fSettings.GetName() + "_data").c_str(), (fSettings.GetFullTitles()).c_str()); Int_t CurrAnalysisBin = 1; for (size_t fds_it = 0; fds_it < FDSamples.size(); ++fds_it) { FDSample &fds = FDSamples[fds_it]; for (Int_t bi_it = 1; bi_it < fds.FDDataHist->GetXaxis()->GetNbins() + 1; ++bi_it) { if ((fds.FitRegion_Min != 0xdeadbeef) && (fds.FDDataHist->GetXaxis()->GetBinUpEdge(bi_it) <= fds.FitRegion_Min)) { continue; } if ((fds.FitRegion_Max != 0xdeadbeef) && (fds.FDDataHist->GetXaxis()->GetBinLowEdge(bi_it) > fds.FitRegion_Max)) { continue; } fDataHist->SetBinContent(CurrAnalysisBin, fds.FDDataHist->GetBinContent(bi_it)); if (UseRateErrors) { fDataHist->SetBinError(CurrAnalysisBin, sqrt(fds.FDDataHist->GetBinContent(bi_it))); } else { fDataHist->SetBinError(CurrAnalysisBin, fds.FDDataHist->GetBinError(bi_it)); } CurrAnalysisBin++; } } fMCHist = static_cast(fDataHist->Clone()); fMCHist->SetNameTitle((fSettings.GetName() + "_MC").c_str(), (fSettings.GetFullTitles()).c_str()); fMCHist->Reset(); } void Smear_SVDUnfold_Propagation_Osc::UpdateChi2Hists() { Int_t CurrAnalysisBin = 1; for (size_t fds_it = 0; fds_it < FDSamples.size(); ++fds_it) { FDSample &fds = FDSamples[fds_it]; for (Int_t bi_it = 1; bi_it < fds.FD_Smeared_Spectrum_Hist->GetXaxis()->GetNbins() + 1; ++bi_it) { if ((fds.FitRegion_Min != 0xdeadbeef) && (fds.FD_Smeared_Spectrum_Hist->GetXaxis()->GetBinUpEdge(bi_it) <= fds.FitRegion_Min)) { continue; } if ((fds.FitRegion_Max != 0xdeadbeef) && (fds.FD_Smeared_Spectrum_Hist->GetXaxis()->GetBinLowEdge(bi_it) > fds.FitRegion_Max)) { continue; } fMCHist->SetBinContent( CurrAnalysisBin, fds.FD_Smeared_Spectrum_Hist->GetBinContent(bi_it)); fMCHist->SetBinError(CurrAnalysisBin, fds.FD_Smeared_Spectrum_Hist->GetBinError(bi_it)); CurrAnalysisBin++; } } } //******************************************************************** Smear_SVDUnfold_Propagation_Osc::Smear_SVDUnfold_Propagation_Osc( nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "Simple measurement class for doing fake data oscillation studies.\n"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetTitle("Osc Studies"); fSettings.SetDescription(descrip); fSettings.SetXTitle("XXX"); fSettings.SetYTitle("Number of events"); fSettings.SetEnuRange(0.0, 50); fSettings.SetAllowedTypes("EVT/SHAPE/DIAG", "EVT/SHAPE/DIAG"); fSettings.DefineAllowedTargets("*"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = GetEventHistogram()->Integral("width") / (fNEvents + 0.); fMCHist = NULL; fDataHist = NULL; ReadExtraConfig(samplekey); AddNDInputs(samplekey); std::vector NDTargets = samplekey.GetListOfChildNodes("NDObs"); for (size_t nd_it = 0; nd_it < NDTargets.size(); ++nd_it) { AddNDInputs(NDTargets[nd_it]); } std::vector FDTargets = samplekey.GetListOfChildNodes("FDObs"); NFDAnalysisBins = 0; if (!FDTargets.size()) { // If no child elements, assume that everything is // contained on the sample element AddFDTarget(samplekey); } else { for (size_t fd_it = 0; fd_it < FDTargets.size(); ++fd_it) { AddFDTarget(FDTargets[fd_it]); NFDAnalysisBins += GetFDSampleNAnalysisBins(fd_it); } } if ((FitRegion_Min != 0xdeadbeef) || (FitRegion_Max != 0xdeadbeef)) { QLOG(SAM, "When unfolding, interested region limited to:"); if (FitRegion_Min != 0xdeadbeef) { QLOG(SAM, "\tE_nu > " << FitRegion_Min); } if (FitRegion_Max != 0xdeadbeef) { QLOG(SAM, "\tE_nu < " << FitRegion_Max); } } SetupNDInputs(); FinaliseFDSamples(); QLOG(SAM, "Set up " << FDSamples.size() << " samples for oscillation analysis with " << NFDAnalysisBins << " analysis bins."); SetupChi2Hists(); SetCovarFromDiagonal(); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; void Smear_SVDUnfold_Propagation_Osc::FillEventVariables(FitEvent *event){}; bool Smear_SVDUnfold_Propagation_Osc::isSignal(FitEvent *event) { return false; } void DumpUnfoldDebugInfo(Smear_SVDUnfold_Propagation_Osc::NDSample &nds, size_t nd_it, size_t truncations) { TDirectory *ogDir = gDirectory; std::stringstream ss; ss.str(""); ss << "DEBUG_FailedInvert_NDSample_" << nd_it; Config::Get().out->mkdir(ss.str().c_str()); Config::Get().out->cd(ss.str().c_str()); ss.str(""); ss << "ND_Smearing_Matrix_" << nd_it; nds.NDToSpectrumSmearingMatrix->Write(ss.str().c_str(), TObject::kOverwrite); ss.str(""); ss << "ND_Inverse_Smearing_Matrix_" << nd_it; SmearceptanceUtils::SVDGetInverse(nds.NDToSpectrumSmearingMatrix, 0) ->Write(ss.str().c_str(), TObject::kOverwrite); ss.str(""); ss << "ND_Inverse_Smearing_Matrix_" << nd_it << "_Trunc_" << truncations; SmearceptanceUtils::SVDGetInverse(nds.NDToSpectrumSmearingMatrix, truncations) ->Write(ss.str().c_str(), TObject::kOverwrite); ss.str(""); ss << "ND_Obs_" << nd_it; nds.NDDataHist->Write(ss.str().c_str(), TObject::kOverwrite); TMatrixD NDToSpectrumResponseMatrix_notrunc = SmearceptanceUtils::GetMatrix( SmearceptanceUtils::SVDGetInverse(nds.NDToSpectrumSmearingMatrix, 0)); nds.NDToSpectrumResponseMatrix.ResizeTo(NDToSpectrumResponseMatrix_notrunc); nds.NDToSpectrumResponseMatrix = NDToSpectrumResponseMatrix_notrunc; SmearceptanceUtils::PushTH1ThroughMatrixWithErrors( nds.NDDataHist, nds.ND_Unfolded_Spectrum_Hist, nds.NDToSpectrumResponseMatrix, 1000, false); ss.str(""); ss << "ND_Unfolded_" << nd_it; nds.ND_Unfolded_Spectrum_Hist->Write(ss.str().c_str(), TObject::kOverwrite); TMatrixD NDToSpectrumResponseMatrix_trunc = SmearceptanceUtils::GetMatrix(SmearceptanceUtils::SVDGetInverse( nds.NDToSpectrumSmearingMatrix, truncations)); nds.NDToSpectrumResponseMatrix.ResizeTo(NDToSpectrumResponseMatrix_trunc); nds.NDToSpectrumResponseMatrix = NDToSpectrumResponseMatrix_trunc; SmearceptanceUtils::PushTH1ThroughMatrixWithErrors( nds.NDDataHist, nds.ND_Unfolded_Spectrum_Hist, nds.NDToSpectrumResponseMatrix, 1000, false); ss.str(""); ss << "ND_Unfolded_" << nd_it << "_Trunc_" << truncations; nds.ND_Unfolded_Spectrum_Hist->Write(ss.str().c_str(), TObject::kOverwrite); if (ogDir) { ogDir->cd(); } else { Config::Get().out->cd(); } } void Smear_SVDUnfold_Propagation_Osc::UnfoldToNDETrueSpectrum(size_t nd_it) { if (nd_it >= NDSamples.size()) { THROW("Attempting to unfold ND sample index " << nd_it << " but only have " << NDSamples.size() << " ND samples."); } NDSample &nds = NDSamples[nd_it]; nds.ND_Unfolded_Spectrum_Hist = new TH1D(*nds.NDToSpectrumSmearingMatrix->ProjectionY()); nds.ND_Unfolded_Spectrum_Hist->SetDirectory(NULL); nds.ND_Unfolded_Spectrum_Hist->Clear(); std::stringstream ss(""); ss << "ND_Unfolded_Spectrum_Hist_" << nd_it; nds.ND_Unfolded_Spectrum_Hist->SetName(ss.str().c_str()); bool HasNegValue = false; int truncations = nds.TruncateStart; do { if (truncations >= nds.TruncateUpTo) { DumpUnfoldDebugInfo(nds, nd_it, truncations); THROW("Unfolded enu spectrum had negative values even after " << truncations << " SVD singular value truncations."); } // Unfold ND ERec -> Enu spectrum // ------------------------------ SmearceptanceUtils::PushTH1ThroughMatrixWithErrors( nds.NDDataHist, nds.ND_Unfolded_Spectrum_Hist, nds.NDToSpectrumResponseMatrix, 1000, false); HasNegValue = false; for (Int_t bi_it = 1; bi_it < nds.ND_Unfolded_Spectrum_Hist->GetXaxis()->GetNbins() + 1; ++bi_it) { if ((FitRegion_Min != 0xdeadbeef) && (nds.ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinUpEdge(bi_it) <= FitRegion_Min)) { continue; } if ((FitRegion_Max != 0xdeadbeef) && (nds.ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinLowEdge(bi_it) > FitRegion_Max)) { continue; } if (nds.ND_Unfolded_Spectrum_Hist->GetBinContent(bi_it) < 0) { HasNegValue = true; QLOG(SAM, "After " << truncations << " truncations, bin " << (bi_it - 1) << " [" << nds.ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinLowEdge( bi_it) << " -- " << nds.ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinUpEdge( bi_it) << " ] has value: " << nds.ND_Unfolded_Spectrum_Hist->GetBinContent(bi_it)); break; } } if (HasNegValue) { TMatrixD NDToSpectrumResponseMatrix_l = SmearceptanceUtils::GetMatrix(SmearceptanceUtils::SVDGetInverse( nds.NDToSpectrumSmearingMatrix, truncations)); nds.NDToSpectrumResponseMatrix.ResizeTo(NDToSpectrumResponseMatrix_l); nds.NDToSpectrumResponseMatrix = NDToSpectrumResponseMatrix_l; } truncations++; } while (HasNegValue); } void Smear_SVDUnfold_Propagation_Osc::PropagateFDSample(size_t fds_it) { if (fds_it >= FDSamples.size()) { THROW("Requested FD sample index " << fds_it << " but only initialised " << FDSamples.size()); } FDSample &fds = FDSamples[fds_it]; // Apply Oscillations // ------------------------------ FitWeight *fw = FitBase::GetRW(); OscWeightEngine *oscWE = dynamic_cast(fw->GetRWEngine(kOSCILLATION)); if (!oscWE) { THROW( "Couldn't load oscillation weight engine for sample: " "Smear_SVDUnfold_Propagation_Osc."); } for (Int_t bi_it = 1; bi_it < fds.NDFD_Corrected_Spectrum_Hist->GetXaxis()->GetNbins() + 1; ++bi_it) { double content_osc = 0; double error_osc = 0; double content_fdnd = 0; double error_fdnd = 0; // Oscillate each ND sample to FD neutrino for (size_t nd_it = 0; nd_it < NDSamples.size(); ++nd_it) { NDSample &nds = NDSamples[nd_it]; double oscWeight = oscWE->CalcWeight( nds.ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinCenter(bi_it), nds.NuPDG, fds.OscillateToPDG); double sample_content_osc = nds.ND_Unfolded_Spectrum_Hist->GetBinContent(bi_it) * oscWeight; double sample_error_osc = nds.ND_Unfolded_Spectrum_Hist->GetBinError(bi_it) * oscWeight; fds.FD_Propagated_Spectrum_Hist_NDSamples[nd_it]->SetBinContent( bi_it, sample_content_osc); fds.FD_Propagated_Spectrum_Hist_NDSamples[nd_it]->SetBinError( bi_it, sample_error_osc); double sample_content_fdnd = sample_content_osc * fds.FDNDRatios[nds.NuPDG] * fds.FDNDMassRatio; double sample_error_fdnd = sample_error_osc * fds.FDNDRatios[nds.NuPDG] * fds.FDNDMassRatio; fds.NDFD_Corrected_Spectrum_Hist_NDSamples[nd_it]->SetBinContent( bi_it, sample_content_fdnd); fds.NDFD_Corrected_Spectrum_Hist_NDSamples[nd_it]->SetBinError( bi_it, sample_error_fdnd); content_osc += sample_content_osc; error_osc += sample_error_osc * sample_error_osc; content_fdnd += sample_content_fdnd; error_fdnd += sample_error_fdnd * sample_error_fdnd; } fds.FD_Propagated_Spectrum_Hist->SetBinContent(bi_it, content_osc); fds.FD_Propagated_Spectrum_Hist->SetBinError(bi_it, sqrt(error_osc)); fds.NDFD_Corrected_Spectrum_Hist->SetBinContent(bi_it, content_fdnd); fds.NDFD_Corrected_Spectrum_Hist->SetBinError(bi_it, sqrt(error_fdnd)); } // Forward fold Spectrum -> ERec FD // ------------------------------ SmearceptanceUtils::PushTH1ThroughMatrixWithErrors( fds.NDFD_Corrected_Spectrum_Hist, fds.FD_Smeared_Spectrum_Hist, fds.SpectrumToFDSmearingMatrix, 1000, true); } void Smear_SVDUnfold_Propagation_Osc::ConvertEventRates(void) { + + /// Get topology weights + /// for each ND sample + /// /// Build NDSample::NDToSpectrumSmearingMatrix from sum of topology smearing histos scaled by topology weights + /// + /// SetupNDInputs(); + /// + /// for each FD sample + /// /// Build FDSample::SpectrumToFDSmearingMatrix_TH2 as above. + /// /// Convert to TMatrixD, FDSample::SpectrumToFDSmearingMatrix + for (size_t fds_it = 0; fds_it < FDSamples.size(); ++fds_it) { PropagateFDSample(fds_it); } UpdateChi2Hists(); } void Smear_SVDUnfold_Propagation_Osc::Write(std::string drawOpt) { TDirectory *ogDir = gDirectory; ConvertEventRates(); FitWeight *fw = FitBase::GetRW(); OscWeightEngine *oscWE = dynamic_cast(fw->GetRWEngine(kOSCILLATION)); if (!oscWE) { THROW( "Couldn't load oscillation weight engine for sample: " "Smear_SVDUnfold_Propagation_Osc."); } oscWE->Print(); // Write ND samples //---------------- for (size_t nd_it = 0; nd_it < NDSamples.size(); ++nd_it) { NDSample &nds = NDSamples[nd_it]; std::stringstream ss(""); ss << "Smear_SVDUnfold_Propagation_Osc_NDSample_" << nd_it; if (ogDir) { ogDir->mkdir(ss.str().c_str()); ogDir->cd(ss.str().c_str()); } else { Config::Get().out->mkdir(ss.str().c_str()); Config::Get().out->cd(ss.str().c_str()); } nds.NDToSpectrumSmearingMatrix->Write("SmearingMatrix_ND", TObject::kOverwrite); nds.NDDataHist->Write("Obs_ND", TObject::kOverwrite); nds.ND_Unfolded_Spectrum_Hist->Write( nds.ND_Unfolded_Spectrum_Hist->GetName(), TObject::kOverwrite); if (ogDir) { ogDir->cd(); } else { Config::Get().out->cd(); } } // Write FD samples //---------------- for (size_t fds_it = 0; fds_it < FDSamples.size(); ++fds_it) { FDSample &fds = FDSamples[fds_it]; std::stringstream ss(""); ss << "Smear_SVDUnfold_Propagation_Osc_FDSample_" << fds_it; if (ogDir) { ogDir->mkdir(ss.str().c_str()); ogDir->cd(ss.str().c_str()); } else { Config::Get().out->mkdir(ss.str().c_str()); Config::Get().out->cd(ss.str().c_str()); } // Calc oscillation probability // ------------------------------ FitWeight *fw = FitBase::GetRW(); OscWeightEngine *oscWE = dynamic_cast(fw->GetRWEngine(kOSCILLATION)); if (!oscWE) { THROW( "Couldn't load oscillation weight engine for sample: " "Smear_SVDUnfold_Propagation_Osc."); } for (size_t nd_it = 0; nd_it < NDSamples.size(); ++nd_it) { NDSample &nds = NDSamples[nd_it]; ss.str(""); ss << "NDSample_" << nd_it << "_contribution"; fds.FD_Propagated_Spectrum_Hist_NDSamples[nd_it]->Write( ss.str().c_str(), TObject::kOverwrite); ss.str(""); ss << "NDSample_" << nd_it << "_FDNDCorrected_contribution"; fds.NDFD_Corrected_Spectrum_Hist_NDSamples[nd_it]->Write( ss.str().c_str(), TObject::kOverwrite); ss.str(""); ss << "NDSample_" << nd_it << "_OscillationProb"; TGraph POsc; POsc.Set(1E4 - 1); double min = nds.ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinLowEdge(1); double step = (nds.ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinUpEdge( nds.ND_Unfolded_Spectrum_Hist->GetXaxis()->GetNbins()) - nds.ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinLowEdge(1)) / double(1E4); for (size_t i = 1; i < 1E4; ++i) { double enu = min + i * step; double ow = oscWE->CalcWeight(enu, nds.NuPDG, fds.OscillateToPDG); if (ow != ow) { std::cout << "Bad osc weight for ENu: " << enu << std::endl; } POsc.SetPoint(i - 1, enu, ow); } POsc.Write(ss.str().c_str(), TObject::kOverwrite); } fds.FDDataHist->Write("Obs_FD", TObject::kOverwrite); fds.FD_Propagated_Spectrum_Hist->Write( fds.FD_Propagated_Spectrum_Hist->GetName(), TObject::kOverwrite); fds.NDFD_Corrected_Spectrum_Hist->Write( fds.NDFD_Corrected_Spectrum_Hist->GetName(), TObject::kOverwrite); fds.SpectrumToFDSmearingMatrix_TH2->Write("SmearingMatrix_FD", TObject::kOverwrite); fds.FD_Smeared_Spectrum_Hist->Write(fds.FD_Smeared_Spectrum_Hist->GetName(), TObject::kOverwrite); if (ogDir) { ogDir->cd(); } else { Config::Get().out->cd(); } } fMCHist->Write("Pred_FD", TObject::kOverwrite); fDataHist->Write("Obs_FD", TObject::kOverwrite); Measurement1D::Write(drawOpt); } diff --git a/src/Smearceptance/ISmearcepter.cxx b/src/Smearceptance/ISmearcepter.cxx index b22d81a..7f12bf8 100644 --- a/src/Smearceptance/ISmearcepter.cxx +++ b/src/Smearceptance/ISmearcepter.cxx @@ -1,11 +1,11 @@ #include "ISmearcepter.h" void ISmearcepter::Setup(nuiskey& nk) { InstanceName = nk.GetS("Name"); ElementName = nk.GetElementName(); QLOG(SAM, "Setting up smearcepter (Type:" << ElementName << ", InstanceName: " - << InstanceName); + << InstanceName << ")."); SpecifcSetup(nk); } diff --git a/src/Smearceptance/Smearcepterton.cxx b/src/Smearceptance/Smearcepterton.cxx index 5db2e91..88a2c7d 100644 --- a/src/Smearceptance/Smearcepterton.cxx +++ b/src/Smearceptance/Smearcepterton.cxx @@ -1,338 +1,342 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #include "Smearcepterton.h" #include "EfficiencyApplicator.h" #include "GaussianSmearer.h" #include "MetaSimpleSmearcepter.h" #include "ThresholdAccepter.h" #include "TrackedMomentumMatrixSmearer.h" #include "VisECoalescer.h" #include #ifdef __USE_DYNSAMPLES__ #include "TRegexp.h" #include // linux #include DynamicSmearceptorFactory::DynamicSmearceptorFactory() : NSmearceptors(0), NManifests(0) { LoadPlugins(); QLOG(FIT, "Loaded " << NSmearceptors << " from " << NManifests << " shared object libraries."); } DynamicSmearceptorFactory* DynamicSmearceptorFactory::glblDSF = NULL; DynamicSmearceptorFactory::PluginManifest::~PluginManifest() { for (size_t i_it = 0; i_it < Instances.size(); ++i_it) { (*(DSF_DestroySmearceptor))(Instances[i_it]); } } std::string EnsureTrailingSlash(std::string const& inp) { if (!inp.length()) { return "/"; } if (inp[inp.length() - 1] == '/') { return inp; } return inp + "/"; } void DynamicSmearceptorFactory::LoadPlugins() { std::vector SearchDirectories; if (Config::HasPar("dynamic_smearceptor.path")) { SearchDirectories = GeneralUtils::ParseToStr( Config::GetParS("dynamic_smearceptor.path"), ":"); } char const* envPath = getenv("NUISANCE_DS_PATH"); if (envPath) { std::vector envPaths = GeneralUtils::ParseToStr(envPath, ":"); for (size_t ep_it = 0; ep_it < envPaths.size(); ++ep_it) { SearchDirectories.push_back(envPaths[ep_it]); } } if (!SearchDirectories.size()) { char const* pwdPath = getenv("PWD"); if (pwdPath) { SearchDirectories.push_back(pwdPath); } } for (size_t sp_it = 0; sp_it < SearchDirectories.size(); ++sp_it) { std::string dirpath = EnsureTrailingSlash(SearchDirectories[sp_it]); QLOG(FIT, "Searching for dynamic smearceptor manifests in: " << dirpath); Ssiz_t len = 0; DIR* dir; struct dirent* ent; dir = opendir(dirpath.c_str()); if (dir != NULL) { TRegexp matchExp("*.so", true); while ((ent = readdir(dir)) != NULL) { if (matchExp.Index(TString(ent->d_name), &len) != Ssiz_t(-1)) { QLOG(FIT, "\tFound shared object: " << ent->d_name << " checking for relevant methods..."); void* dlobj = dlopen((dirpath + ent->d_name).c_str(), RTLD_NOW | RTLD_GLOBAL); char const* dlerr_cstr = dlerror(); std::string dlerr; if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { ERROR(WRN, "\tDL Load Error: " << dlerr); continue; } PluginManifest plgManif; plgManif.dllib = dlobj; plgManif.soloc = (dirpath + ent->d_name); plgManif.DSF_NSmearceptors = reinterpret_cast( dlsym(dlobj, "DSF_NSmearceptors")); dlerr = ""; dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { ERROR(WRN, "\tFailed to load symbol \"DSF_NSmearceptors\" from " << (dirpath + ent->d_name) << ": " << dlerr); dlclose(dlobj); continue; } plgManif.DSF_GetSmearceptorName = reinterpret_cast( dlsym(dlobj, "DSF_GetSmearceptorName")); dlerr = ""; dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { ERROR(WRN, "\tFailed to load symbol \"DSF_GetSmearceptorName\" from " << (dirpath + ent->d_name) << ": " << dlerr); dlclose(dlobj); continue; } plgManif.DSF_GetSmearceptor = reinterpret_cast( dlsym(dlobj, "DSF_GetSmearceptor")); dlerr = ""; dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { ERROR(WRN, "\tFailed to load symbol \"DSF_GetSmearceptor\" from " << (dirpath + ent->d_name) << ": " << dlerr); dlclose(dlobj); continue; } plgManif.DSF_DestroySmearceptor = reinterpret_cast( dlsym(dlobj, "DSF_DestroySmearceptor")); dlerr = ""; dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { ERROR(WRN, "Failed to load symbol \"DSF_DestroySmearceptor\" from " << (dirpath + ent->d_name) << ": " << dlerr); dlclose(dlobj); continue; } plgManif.NSmearceptors = (*(plgManif.DSF_NSmearceptors))(); QLOG(FIT, "\tSuccessfully loaded dynamic smearceptor manifest: " << plgManif.soloc << ". Contains " << plgManif.NSmearceptors << " smearceptors."); for (size_t smp_it = 0; smp_it < plgManif.NSmearceptors; ++smp_it) { char const* smp_name = (*(plgManif.DSF_GetSmearceptorName))(smp_it); if (!smp_name) { THROW("Could not load smearceptor " << smp_it << " / " << plgManif.NSmearceptors << " from " << plgManif.soloc); } if (Smearceptors.count(smp_name)) { ERROR(WRN, "Already loaded a smearceptor named: \"" << smp_name << "\". cannot load duplciates. This " "smearceptor will be skipped."); continue; } plgManif.SmearceptorsProvided.push_back(smp_name); Smearceptors[smp_name] = std::make_pair(plgManif.soloc, smp_it); QLOG(FIT, "\t\t" << smp_name); } if (plgManif.SmearceptorsProvided.size()) { Manifests[plgManif.soloc] = plgManif; NSmearceptors += plgManif.SmearceptorsProvided.size(); NManifests++; } else { dlclose(dlobj); } } } closedir(dir); } else { ERROR(WRN, "Tried to open non-existant directory."); } } } DynamicSmearceptorFactory& DynamicSmearceptorFactory::Get() { if (!glblDSF) { glblDSF = new DynamicSmearceptorFactory(); } return *glblDSF; } void DynamicSmearceptorFactory::Print() { std::map > ManifestSmearceptors; for (std::map >::iterator smp_it = Smearceptors.begin(); smp_it != Smearceptors.end(); ++smp_it) { if (!ManifestSmearceptors.count(smp_it->second.first)) { ManifestSmearceptors[smp_it->second.first] = std::vector(); } ManifestSmearceptors[smp_it->second.first].push_back(smp_it->first); } QLOG(FIT, "Dynamic smearceptor manifest: "); for (std::map >::iterator m_it = ManifestSmearceptors.begin(); m_it != ManifestSmearceptors.end(); ++m_it) { QLOG(FIT, "\tLibrary " << m_it->first << " contains: "); for (size_t s_it = 0; s_it < m_it->second.size(); ++s_it) { QLOG(FIT, "\t\t" << m_it->second[s_it]); } } } bool DynamicSmearceptorFactory::HasSmearceptor(std::string const& name) { return Smearceptors.count(name); } bool DynamicSmearceptorFactory::HasSmearceptor(nuiskey& smearceptorkey) { return HasSmearceptor(smearceptorkey.GetElementName()); } ISmearcepter* DynamicSmearceptorFactory::CreateSmearceptor( nuiskey& smearceptorkey) { if (!HasSmearceptor(smearceptorkey)) { ERROR(WRN, "Asked to load unknown smearceptor: \"" << smearceptorkey.GetElementName() << "\"."); return NULL; } std::pair smearceptor = Smearceptors[smearceptorkey.GetElementName()]; QLOG(SAM, "\tLoading smearceptor " << smearceptor.second << " from " << smearceptor.first); - ISmearcepter* smear = (*(Manifests[smearceptor.first].DSF_GetSmearceptor))( + ISmearcepter* smear = (*(Manifests[smearceptor.first].DSF_GetSmearceptor))( smearceptor.second, &smearceptorkey); return smear; } DynamicSmearceptorFactory::~DynamicSmearceptorFactory() { Manifests.clear(); } #endif Smearcepterton* Smearcepterton::_inst = NULL; Smearcepterton& Smearcepterton::Get() { if (!_inst) { _inst = new Smearcepterton(); } return *_inst; } Smearcepterton::Smearcepterton() { InitialiserSmearcepters(); } void Smearcepterton::InitialiserSmearcepters() { // hard coded list of tag name -> smearcepter factories, add here to add your // own. std::map factories; factories["ThresholdAccepter"] = &BuildSmearcepter; factories["EfficiencyApplicator"] = &BuildSmearcepter; factories["GaussianSmearer"] = &BuildSmearcepter; factories["TrackedMomentumMatrixSmearer"] = &BuildSmearcepter; factories["VisECoalescer"] = &BuildSmearcepter; factories["MetaSimpleSmearcepter"] = &BuildSmearcepter; std::vector smearcepterBlocks = Config::QueryKeys("smearcepters"); for (size_t smearB_it = 0; smearB_it < smearcepterBlocks.size(); ++smearB_it) { std::vector smearcepters = smearcepterBlocks[smearB_it].GetListOfChildNodes(); for (size_t smear_it = 0; smear_it < smearcepters.size(); ++smear_it) { std::string const& smearType = smearcepters[smear_it].GetElementName(); ISmearcepter* smearer = NULL; if (DynamicSmearceptorFactory::Get().HasSmearceptor(smearType)) { smearer = DynamicSmearceptorFactory::Get().CreateSmearceptor( smearcepters[smear_it]); } else { if (!factories.count(smearType)) { ERROR(WRN, "No known smearer accepts elements named: \"" << smearType << "\""); continue; } smearer = factories[smearType](smearcepters[smear_it]); } - if(!smearer){ + if (!smearer) { THROW("Failed to load smearceptor."); } + if (!smearer->GetName().length()) { + THROW("Smearcepter type " << smearer->GetElementName() + << " had no instance name."); + } Smearcepters[smearer->GetName()] = smearer; QLOG(FIT, "Configured smearer named: " << smearer->GetName() << " of type: " << smearer->GetElementName()); } } }