diff --git a/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx b/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx index 014f927..4d56759 100644 --- a/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx +++ b/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx @@ -1,924 +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]; + //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/MCStudies/Smearceptance_Tester.cxx b/src/MCStudies/Smearceptance_Tester.cxx index e19e54b..2904770 100644 --- a/src/MCStudies/Smearceptance_Tester.cxx +++ b/src/MCStudies/Smearceptance_Tester.cxx @@ -1,881 +1,881 @@ // 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 "Smearceptance_Tester.h" #include "SmearceptanceUtils.h" #include "Smearcepterton.h" //#define DEBUG_SMEARTESTER //******************************************************************** /// @brief Class to perform smearceptance MC Studies on a custom measurement Smearceptance_Tester::Smearceptance_Tester(nuiskey samplekey) { //******************************************************************** samplekey.Print(); // Sample overview --------------------------------------------------- std::string descrip = "Simple measurement class for producing an event summary tree of smeared " "events.\n"; if (Config::HasPar("NPOT")) { samplekey.SetS("NPOT", Config::GetParS("NPOT")); } if (Config::HasPar("FluxIntegralOverride")) { samplekey.SetS("FluxIntegralOverride", Config::GetParS("FluxIntegralOverride")); } if (Config::HasPar("TargetVolume")) { samplekey.SetS("TargetVolume", Config::GetParS("TargetVolume")); } if (Config::HasPar("TargetMaterialDensity")) { samplekey.SetS("TargetMaterialDensity", Config::GetParS("TargetMaterialDensity")); } OutputSummaryTree = true; if (Config::HasPar("smear.OutputSummaryTree")) { OutputSummaryTree = Config::GetParI("smear.OutputSummaryTree"); } // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetTitle("Smearceptance Studies"); fSettings.SetDescription(descrip); fSettings.SetXTitle("XXX"); fSettings.SetYTitle("Number of events"); fSettings.SetEnuRange(0.0, 1E5); fSettings.SetAllowedTypes("EVT/SHAPE/DIAG", "EVT/SHAPE/DIAG"); fSettings.DefineAllowedTargets("*"); fSettings.DefineAllowedSpecies("*"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) / TotalIntegratedFlux(); // Measurement Details std::vector splitName = GeneralUtils::ParseToStr(fName, "_"); - size_t firstUS = fName.find_first_of("_"); + //size_t firstUS = fName.find_first_of("_"); std::string smearceptorName = samplekey.GetS("smearceptor"); QLOG(SAM, "Using smearceptor: " << smearceptorName << " (parsed from: " << fName << ")."); fDataHist = new TH1D(("empty_data"), ("empty-data"), 1, 0, 1); SetupDefaultHist(); fFullCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist); covar = StatUtils::GetInvert(fFullCovar); eventVariables = NULL; QLOG(SAM, "Smearceptance Flux Scaling Factor = " << fScaleFactor); if (fScaleFactor <= 0.0) { ERROR(WRN, "SCALE FACTOR TOO LOW "); sleep(20); } // Setup our TTrees AddEventVariablesToTree(); smearceptor = &Smearcepterton::Get().GetSmearcepter(smearceptorName); Int_t RecNBins = 20, TrueNBins = 20; double RecBinL = 0xdeadbeef, TrueBinL = 0, RecBinH = 10, TrueBinH = 10; if (Config::HasPar("smear.reconstructed.binning")) { std::vector args = GeneralUtils::ParseToStr( Config::GetParS("smear.reconstructed.binning"), ","); RecNBins = GeneralUtils::StrToInt(args[0]); RecBinL = GeneralUtils::StrToDbl(args[1]); RecBinH = GeneralUtils::StrToDbl(args[2]); TrueNBins = RecNBins; TrueBinL = RecBinL; TrueBinH = RecBinH; } if (Config::HasPar("smear.true.binning")) { std::vector args = GeneralUtils::ParseToStr(Config::GetParS("smear.true.binning"), ","); TrueNBins = GeneralUtils::StrToInt(args[0]); TrueBinL = GeneralUtils::StrToDbl(args[1]); TrueBinH = GeneralUtils::StrToDbl(args[2]); } SVDTruncation = 0; if (Config::HasPar("smear.true.binning")) { SVDTruncation = Config::GetParI("smear.SVD.truncation"); QLOG(SAM, "Applying SVD truncation of: " << SVDTruncation) } ETrueDistrib = NULL; ETrueDistrib_noweight = NULL; ERecDistrib = NULL; RecoSmear = NULL; if (RecBinL != 0xdeadbeef) { QLOG(SAM, "Using binning True: " << TrueNBins << ", [" << TrueBinL << " -- " << TrueBinH << "], Rec: " << RecNBins << ", [" << RecBinL << " -- " << RecBinH << "]"); ETrueDistrib = new TH1D("ELep_rate", ";True E_{#nu};Count", TrueNBins, TrueBinL, TrueBinH); ETrueDistrib_noweight = new TH1D("ELep_rate_noweight", ";True E_{#nu};Count", TrueNBins, TrueBinL, TrueBinH); ERecDistrib = new TH1D("ELepRec_rate", ";Rec E_{#nu};Count", RecNBins, RecBinL, RecBinH); ETrueDistrib->Sumw2(); ERecDistrib->Sumw2(); RecoSmear = new TH2D("ELepHadVis_Recon", ";True E_{#nu};Recon. E_{#nu}", RecNBins, RecBinL, RecBinH, TrueNBins, TrueBinL, TrueBinH); RecoSmear->Sumw2(); } // Final setup --------------------------------------------------- FinaliseMeasurement(); } void Smearceptance_Tester::AddEventVariablesToTree() { if (OutputSummaryTree) { // Setup the TTree to save everything if (!eventVariables) { Config::Get().out->cd(); eventVariables = new TTree((fName + "_VARS").c_str(), (fName + "_VARS").c_str()); } LOG(SAM) << "Adding Event Variables" << std::endl; eventVariables->Branch("Omega_true", &Omega_true, "Omega_true/F"); eventVariables->Branch("Q2_true", &Q2_true, "Q2_true/F"); eventVariables->Branch("Mode_true", &Mode_true, "Mode_true/I"); eventVariables->Branch("EISLep_true", &EISLep_true, "EISLep_true/F"); eventVariables->Branch("HMFS_clep_true", &HMFS_clep_true); eventVariables->Branch("HMFS_pip_true", &HMFS_pip_true); eventVariables->Branch("HMFS_pim_true", &HMFS_pim_true); eventVariables->Branch("HMFS_cpi_true", &HMFS_cpi_true); eventVariables->Branch("HMFS_pi0_true", &HMFS_pi0_true); eventVariables->Branch("HMFS_cK_true", &HMFS_cK_true); eventVariables->Branch("HMFS_K0_true", &HMFS_K0_true); eventVariables->Branch("HMFS_p_true", &HMFS_p_true); eventVariables->Branch("KEFSHad_cpip_true", &KEFSHad_cpip_true, "KEFSHad_cpip_true/F"); eventVariables->Branch("KEFSHad_cpim_true", &KEFSHad_cpim_true, "KEFSHad_cpim_true/F"); eventVariables->Branch("KEFSHad_cpi_true", &KEFSHad_cpi_true, "KEFSHad_cpi_true/F"); eventVariables->Branch("TEFSHad_pi0_true", &TEFSHad_pi0_true, "TEFSHad_pi0_true/F"); eventVariables->Branch("KEFSHad_cK_true", &KEFSHad_cK_true, "KEFSHad_cK_true/F"); eventVariables->Branch("KEFSHad_K0_true", &KEFSHad_K0_true, "KEFSHad_K0_true/F"); eventVariables->Branch("KEFSHad_p_true", &KEFSHad_p_true, "KEFSHad_p_true/F"); eventVariables->Branch("KEFSHad_n_true", &KEFSHad_n_true, "KEFSHad_n_true/F"); eventVariables->Branch("EFSHad_true", &EFSHad_true, "EFSHad_true/F"); eventVariables->Branch("EFSChargedEMHad_true", &EFSChargedEMHad_true, "EFSChargedEMHad_true/F"); eventVariables->Branch("EFSLep_true", &EFSLep_true, "EFSLep_true/F"); eventVariables->Branch("EFSgamma_true", &EFSgamma_true, "EFSgamma_true/F"); eventVariables->Branch("PDGISLep_true", &PDGISLep_true, "PDGISLep_true/I"); eventVariables->Branch("PDGFSLep_true", &PDGFSLep_true, "PDGFSLep_true/I"); eventVariables->Branch("Nprotons_true", &Nprotons_true, "Nprotons_true/I"); eventVariables->Branch("Nneutrons_true", &Nneutrons_true, "Nneutrons_true/I"); eventVariables->Branch("Ncpiplus_true", &Ncpiplus_true, "Ncpiplus_true/I"); eventVariables->Branch("Ncpiminus_true", &Ncpiminus_true, "Ncpiminus_true/I"); eventVariables->Branch("Ncpi_true", &Ncpi_true, "Ncpi_true/I"); eventVariables->Branch("Npi0_true", &Npi0_true, "Npi0_true/I"); eventVariables->Branch("NcK_true", &NcK_true, "NcK_true/I"); eventVariables->Branch("NK0_true", &NK0_true, "NK0_true/I"); eventVariables->Branch("HMFS_clep_rec", &HMFS_clep_rec); eventVariables->Branch("HMFS_pip_rec", &HMFS_pip_rec); eventVariables->Branch("HMFS_pim_rec", &HMFS_pim_rec); eventVariables->Branch("HMFS_cpi_rec", &HMFS_cpi_rec); eventVariables->Branch("HMFS_pi0_rec", &HMFS_pi0_rec); eventVariables->Branch("HMFS_cK_rec", &HMFS_cK_rec); eventVariables->Branch("HMFS_K0_rec", &HMFS_K0_rec); eventVariables->Branch("HMFS_p_rec", &HMFS_p_rec); eventVariables->Branch("KEFSHad_cpip_rec", &KEFSHad_cpip_rec, "KEFSHad_cpip_rec/F"); eventVariables->Branch("KEFSHad_cpim_rec", &KEFSHad_cpim_rec, "KEFSHad_cpim_rec/F"); eventVariables->Branch("KEFSHad_cpi_rec", &KEFSHad_cpi_rec, "KEFSHad_cpi_rec/F"); eventVariables->Branch("TEFSHad_pi0_rec", &TEFSHad_pi0_rec, "TEFSHad_pi0_rec/F"); eventVariables->Branch("KEFSHad_cK_rec", &KEFSHad_cK_rec, "KEFSHad_cK_rec/F"); eventVariables->Branch("KEFSHad_K0_rec", &KEFSHad_K0_rec, "KEFSHad_K0_rec/F"); eventVariables->Branch("KEFSHad_p_rec", &KEFSHad_p_rec, "KEFSHad_p_rec/F"); eventVariables->Branch("KEFSHad_n_rec", &KEFSHad_n_rec, "KEFSHad_n_rec/F"); eventVariables->Branch("EFSHad_rec", &EFSHad_rec, "EFSHad_rec/F"); eventVariables->Branch("EFSLep_rec", &EFSLep_rec, "EFSLep_rec/F"); eventVariables->Branch("EFSVis_cpip", &EFSVis_cpip, "EFSVis_cpip/F"); eventVariables->Branch("EFSVis_cpim", &EFSVis_cpim, "EFSVis_cpim/F"); eventVariables->Branch("EFSVis_cpi", &EFSVis_cpi, "EFSVis_cpi/F"); eventVariables->Branch("EFSVis_pi0", &EFSVis_pi0, "EFSVis_pi0/F"); eventVariables->Branch("EFSVis_cK", &EFSVis_cK, "EFSVis_cK/F"); eventVariables->Branch("EFSVis_K0", &EFSVis_K0, "EFSVis_K0/F"); eventVariables->Branch("EFSVis_p", &EFSVis_p, "EFSVis_p/F"); eventVariables->Branch("EFSVis_n", &EFSVis_n, "EFSVis_n/F"); eventVariables->Branch("EFSVis_gamma", &EFSVis_gamma, "EFSVis_gamma/F"); eventVariables->Branch("EFSVis_other", &EFSVis_other, "EFSVis_other/F"); eventVariables->Branch("EFSVis", &EFSVis, "EFSVis/F"); eventVariables->Branch("FSCLep_seen", &FSCLep_seen, "FSCLep_seen/I"); eventVariables->Branch("Nprotons_seen", &Nprotons_seen, "Nprotons_seen/I"); eventVariables->Branch("Nneutrons_seen", &Nneutrons_seen, "Nneutrons_seen/I"); eventVariables->Branch("Ncpip_seen", &Ncpip_seen, "Ncpip_seen/I"); eventVariables->Branch("Ncpim_seen", &Ncpim_seen, "Ncpim_seen/I"); eventVariables->Branch("Ncpi_seen", &Ncpi_seen, "Ncpi_seen/I"); eventVariables->Branch("Npi0_seen", &Npi0_seen, "Npi0_seen/I"); eventVariables->Branch("NcK_seen", &NcK_seen, "NcK_seen/I"); eventVariables->Branch("NK0_seen", &NK0_seen, "NK0_seen/I"); eventVariables->Branch("Nothers_seen", &Nothers_seen, "Nothers_seen/I"); eventVariables->Branch("EISLep_QE_rec", &EISLep_QE_rec, "EISLep_QE_rec/F"); eventVariables->Branch("EISLep_LepHad_rec", &EISLep_LepHad_rec, "EISLep_LepHad_rec/F"); eventVariables->Branch("EISLep_LepHadVis_rec", &EISLep_LepHadVis_rec, "EISLep_LepHadVis_rec/F"); eventVariables->Branch("Nprotons_contributed", &Nprotons_contributed, "Nprotons_contributed/I"); eventVariables->Branch("Nneutrons_contributed", &Nneutrons_contributed, "Nneutrons_contributed/I"); eventVariables->Branch("Ncpip_contributed", &Ncpip_contributed, "Ncpip_contributed/I"); eventVariables->Branch("Ncpim_contributed", &Ncpim_contributed, "Ncpim_contributed/I"); eventVariables->Branch("Ncpi_contributed", &Ncpi_contributed, "Ncpi_contributed/I"); eventVariables->Branch("Npi0_contributed", &Npi0_contributed, "Npi0_contributed/I"); eventVariables->Branch("NcK_contributed", &NcK_contributed, "NcK_contributed/I"); eventVariables->Branch("NK0_contributed", &NK0_contributed, "NK0_contributed/I"); eventVariables->Branch("Ngamma_contributed", &Ngamma_contributed, "Ngamma_contributed/I"); eventVariables->Branch("Nothers_contibuted", &Nothers_contibuted, "Nothers_contibuted/I"); eventVariables->Branch("Weight", &Weight, "Weight/F"); eventVariables->Branch("RWWeight", &RWWeight, "RWWeight/F"); eventVariables->Branch("InputWeight", &InputWeight, "InputWeight/F"); eventVariables->Branch("FluxWeight", &FluxWeight, "FluxWeight/F"); eventVariables->Branch("EffWeight", &EffWeight, "EffWeight/F"); xsecScaling = fScaleFactor; eventVariables->Branch("xsecScaling", &xsecScaling, "xsecScaling/F"); eventVariables->Branch("flagCCINC_true", &flagCCINC_true, "flagCCINC_true/O"); eventVariables->Branch("flagCC0K_true", &flagCC0K_true, "flagCC0K_true/O"); eventVariables->Branch("flagCC0Pi_true", &flagCC0Pi_true, "flagCC0Pi_true/O"); eventVariables->Branch("flagCC1Pi_true", &flagCC1Pi_true, "flagCC1Pi_true/O"); eventVariables->Branch("flagCCINC_rec", &flagCCINC_rec, "flagCCINC_rec/O"); eventVariables->Branch("flagCC0K_rec", &flagCC0K_rec, "flagCC0K_rec/O"); eventVariables->Branch("flagCC0Pi_rec", &flagCC0Pi_rec, "flagCC0Pi_rec/O"); eventVariables->Branch("flagCC1Pi_rec", &flagCC1Pi_rec, "flagCC1Pi_rec/O"); } PredEvtRateWeight = 1; if (fEvtRateScaleFactor != 0xdeadbeef) { if (OutputSummaryTree) { eventVariables->Branch("PredEvtRateWeight", &PredEvtRateWeight, "PredEvtRateWeight/F"); } PredEvtRateWeight = fScaleFactor * fEvtRateScaleFactor; } } template int CountNPdgsSeen(RecoInfo ri, int const (&pdgs)[N]) { int sum = 0; for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) { sum += std::count(ri.RecObjClass.begin(), ri.RecObjClass.end(), pdgs[pdg_it]); } return sum; } template int CountNNotPdgsSeen(RecoInfo ri, int const (&pdgs)[N]) { int sum = 0; for (size_t p_it = 0; p_it < ri.RecObjClass.size(); ++p_it) { if (!std::count(pdgs, pdgs + N, ri.RecObjClass[p_it])) { sum++; } } return sum; } template int CountNPdgsContributed(RecoInfo ri, int const (&pdgs)[N]) { int sum = 0; for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) { sum += std::count(ri.TrueContribPDGs.begin(), ri.TrueContribPDGs.end(), pdgs[pdg_it]); } return sum; } template int CountNNotPdgsContributed(RecoInfo ri, int const (&pdgs)[N]) { int sum = 0; for (size_t p_it = 0; p_it < ri.TrueContribPDGs.size(); ++p_it) { if (!std::count(pdgs, pdgs + N, ri.TrueContribPDGs[p_it])) { sum++; } } return sum; } TLorentzVector GetHMFSRecParticles(RecoInfo ri, int pdg) { TLorentzVector mom(0, 0, 0, 0); for (size_t p_it = 0; p_it < ri.RecObjMom.size(); ++p_it) { if ((ri.RecObjClass[p_it] == pdg) && (mom.Mag() < ri.RecObjMom[p_it].Mag())) { mom.SetXYZM(ri.RecObjMom[p_it].X(), ri.RecObjMom[p_it].Y(), ri.RecObjMom[p_it].Z(), PhysConst::GetMass(ri.RecObjClass[p_it]) * 1.0E3); } } return mom; } template double SumKE_RecoInfo(RecoInfo ri, int const (&pdgs)[N], double mass) { double sum = 0; for (size_t p_it = 0; p_it < ri.RecObjMom.size(); ++p_it) { if (!std::count(pdgs, pdgs + N, ri.RecObjClass[p_it])) { // If we don't care about this // particle type. continue; } sum += sqrt(ri.RecObjMom[p_it].Mag2() + mass * mass) - mass; } return sum; } template double SumTE_RecoInfo(RecoInfo ri, int const (&pdgs)[N], double mass) { double sum = 0; for (size_t p_it = 0; p_it < ri.RecObjMom.size(); ++p_it) { if (!std::count(pdgs, pdgs + N, ri.RecObjClass[p_it])) { // If we don't care about this // particle type. continue; } sum += sqrt(ri.RecObjMom[p_it].Mag2() + mass * mass); } return sum; } template double SumVisE_RecoInfo(RecoInfo ri, int const (&pdgs)[N]) { double sum = 0; for (size_t p_it = 0; p_it < ri.RecVisibleEnergy.size(); ++p_it) { if (!std::count(pdgs, pdgs + N, ri.TrueContribPDGs[p_it])) { // If we don't care about this // particle type. continue; } sum += ri.RecVisibleEnergy[p_it]; } return sum; } template double SumVisE_RecoInfo_NotPdgs(RecoInfo ri, int const (&pdgs)[N]) { double sum = 0; for (size_t p_it = 0; p_it < ri.RecVisibleEnergy.size(); ++p_it) { if (std::count(pdgs, pdgs + N, ri.TrueContribPDGs[p_it])) { // If we know about this // particle type. continue; } sum += ri.RecVisibleEnergy[p_it]; } return sum; } //******************************************************************** void Smearceptance_Tester::FillEventVariables(FitEvent *event) { //******************************************************************** static int const cpipPDG[] = {211}; static int const cpimPDG[] = {-211}; static int const pi0PDG[] = {111}; static int const cKPDG[] = {321, -321}; static int const K0PDG[] = {311, 310, 130}; static int const ProtonPDG[] = {2212}; static int const NeutronPDG[] = {2112}; static int const GammaPDG[] = {22}; static int const CLeptonPDGs[] = {11, 13, 15, -11, -13, -15}; static int const ExplicitPDGs[] = {211, -211, 111, 321, -321, 311, 310, 130, 2212, 2112, 22, 11, 13, 15, 12, 14, 16, -11, -13, -15, -12, -14, -16}; RecoInfo *ri = smearceptor->Smearcept(event); //** START Pions HMFS_clep_true = TLorentzVector(0, 0, 0, 0); HMFS_clep_rec = TLorentzVector(0, 0, 0, 0); FitParticle *fsCLep = event->GetHMFSParticle(CLeptonPDGs); if (fsCLep) { HMFS_clep_true = fsCLep->P4(); HMFS_clep_rec = GetHMFSRecParticles(*ri, fsCLep->PDG()); } //** END Charged leptons //** START Pions HMFS_pip_true = TLorentzVector(0, 0, 0, 0); HMFS_pip_rec = TLorentzVector(0, 0, 0, 0); FitParticle *fsPip = event->GetHMFSPiPlus(); if (fsPip) { HMFS_pip_true = fsPip->P4(); HMFS_pip_rec = GetHMFSRecParticles(*ri, fsPip->PDG()); } HMFS_pim_true = TLorentzVector(0, 0, 0, 0); HMFS_pim_rec = TLorentzVector(0, 0, 0, 0); FitParticle *fsPim = event->GetHMFSPiMinus(); if (fsPim) { HMFS_pim_true = fsPim->P4(); HMFS_pim_rec = GetHMFSRecParticles(*ri, fsPim->PDG()); } HMFS_cpi_true = TLorentzVector(0, 0, 0, 0); HMFS_cpi_rec = TLorentzVector(0, 0, 0, 0); if (fsPip || fsPim) { if (!fsPip) { HMFS_cpi_true = HMFS_pim_true; HMFS_cpi_rec = HMFS_pim_rec; } else if (!fsPim) { HMFS_cpi_true = HMFS_pip_true; HMFS_cpi_rec = HMFS_pip_rec; } else { HMFS_cpi_true = (fsPip->p2() > fsPim->p2()) ? HMFS_pip_true : HMFS_pim_true; HMFS_cpi_rec = (fsPip->p2() > fsPim->p2()) ? HMFS_pip_rec : HMFS_pim_rec; } } HMFS_pi0_true = TLorentzVector(0, 0, 0, 0); HMFS_pi0_rec = TLorentzVector(0, 0, 0, 0); FitParticle *fsPi0 = event->GetHMFSPiZero(); if (fsPi0) { HMFS_pi0_true = fsPi0->P4(); HMFS_pi0_rec = GetHMFSRecParticles(*ri, fsPi0->PDG()); } //** END Pions //** START Kaons HMFS_cK_true = TLorentzVector(0, 0, 0, 0); HMFS_cK_rec = TLorentzVector(0, 0, 0, 0); FitParticle *fscK = event->GetHMFSParticle(cKPDG); if (fscK) { HMFS_cK_true = fscK->P4(); HMFS_cK_rec = GetHMFSRecParticles(*ri, fscK->PDG()); } HMFS_K0_true = TLorentzVector(0, 0, 0, 0); HMFS_K0_rec = TLorentzVector(0, 0, 0, 0); FitParticle *fsK0 = event->GetHMFSParticle(K0PDG); if (fsK0) { HMFS_K0_true = fsK0->P4(); HMFS_K0_rec = GetHMFSRecParticles(*ri, fsK0->PDG()); } //** END Kaons //** START Nucleons HMFS_p_true = TLorentzVector(0, 0, 0, 0); HMFS_p_rec = TLorentzVector(0, 0, 0, 0); FitParticle *fsP = event->GetHMFSProton(); if (fsP) { HMFS_p_true = fsP->P4(); HMFS_p_rec = GetHMFSRecParticles(*ri, fsP->PDG()); } TLorentzVector FourMomentumTransfer = (event->GetHMISAnyLeptons()->P4() - event->GetHMFSAnyLeptons()->P4()); Omega_true = FourMomentumTransfer.E(); Q2_true = -1 * FourMomentumTransfer.Mag2(); Mode_true = event->Mode; EISLep_true = event->GetHMISAnyLeptons()->E(); KEFSHad_cpip_true = FitUtils::SumTE_PartVect(event->GetAllFSPiPlus()); KEFSHad_cpim_true = FitUtils::SumTE_PartVect(event->GetAllFSPiMinus()); KEFSHad_cpi_true = KEFSHad_cpip_true + KEFSHad_cpim_true; TEFSHad_pi0_true = FitUtils::SumTE_PartVect(event->GetAllFSPiZero()); KEFSHad_cK_true = FitUtils::SumTE_PartVect(event->GetAllFSParticle(cKPDG)); KEFSHad_K0_true = FitUtils::SumTE_PartVect(event->GetAllFSParticle(K0PDG)); KEFSHad_p_true = FitUtils::SumKE_PartVect(event->GetAllFSProton()); KEFSHad_n_true = FitUtils::SumKE_PartVect(event->GetAllFSNeutron()); EFSHad_true = KEFSHad_cpi_true + TEFSHad_pi0_true + KEFSHad_p_true + KEFSHad_n_true; EFSChargedEMHad_true = KEFSHad_cpi_true + TEFSHad_pi0_true + KEFSHad_p_true + KEFSHad_cK_true + KEFSHad_K0_true; EFSLep_true = event->GetHMFSAnyLeptons()->E(); EFSgamma_true = FitUtils::SumTE_PartVect(event->GetAllFSPhoton()); PDGISLep_true = event->GetHMISAnyLeptons()->PDG(); PDGFSLep_true = event->GetHMFSAnyLeptons()->PDG(); Nprotons_true = event->GetAllFSProton().size(); Nneutrons_true = event->GetAllFSNeutron().size(); Ncpiplus_true = event->GetAllFSPiPlus().size(); Ncpiminus_true = event->GetAllFSPiMinus().size(); Ncpi_true = Ncpiplus_true + Ncpiminus_true; Npi0_true = event->GetAllFSPiZero().size(); NcK_true = event->GetAllFSParticle(cKPDG).size(); NK0_true = event->GetAllFSParticle(K0PDG).size(); KEFSHad_cpip_rec = SumKE_RecoInfo(*ri, cpipPDG, PhysConst::mass_cpi * PhysConst::mass_MeV); KEFSHad_cpim_rec = SumKE_RecoInfo(*ri, cpimPDG, PhysConst::mass_cpi * PhysConst::mass_MeV); KEFSHad_cpi_rec = KEFSHad_cpip_rec + KEFSHad_cpim_rec; TEFSHad_pi0_rec = SumTE_RecoInfo(*ri, pi0PDG, PhysConst::mass_pi0 * PhysConst::mass_MeV); KEFSHad_cK_rec = SumKE_RecoInfo(*ri, cKPDG, PhysConst::mass_cK * PhysConst::mass_MeV); KEFSHad_K0_rec = SumKE_RecoInfo(*ri, K0PDG, PhysConst::mass_K0 * PhysConst::mass_MeV); KEFSHad_p_rec = SumKE_RecoInfo(*ri, ProtonPDG, PhysConst::mass_proton * PhysConst::mass_MeV); KEFSHad_n_rec = SumKE_RecoInfo(*ri, NeutronPDG, PhysConst::mass_neutron * PhysConst::mass_MeV); EFSHad_rec = KEFSHad_cpi_rec + TEFSHad_pi0_rec + KEFSHad_p_rec + KEFSHad_n_rec; TLorentzVector FSLepMom_rec(0, 0, 0, 0); if (event->GetHMFSAnyLeptons()) { FSLepMom_rec = GetHMFSRecParticles(*ri, event->GetHMFSAnyLeptons()->PDG()); EFSLep_rec = FSLepMom_rec.E(); } else { EFSLep_rec = 0; } EFSVis_cpip = SumVisE_RecoInfo(*ri, cpipPDG); EFSVis_cpim = SumVisE_RecoInfo(*ri, cpimPDG); EFSVis_cpi = EFSVis_cpip + EFSVis_cpim; EFSVis_pi0 = SumVisE_RecoInfo(*ri, pi0PDG); EFSVis_cK = SumVisE_RecoInfo(*ri, cKPDG); EFSVis_K0 = SumVisE_RecoInfo(*ri, K0PDG); EFSVis_p = SumVisE_RecoInfo(*ri, ProtonPDG); EFSVis_n = SumVisE_RecoInfo(*ri, NeutronPDG); EFSVis_gamma = SumVisE_RecoInfo(*ri, GammaPDG); EFSVis_other = SumVisE_RecoInfo_NotPdgs(*ri, ExplicitPDGs); EFSVis = EFSVis_cpi + EFSVis_pi0 + EFSVis_p + EFSVis_n + EFSVis_gamma + EFSVis_cK + EFSVis_K0; FSCLep_seen = CountNPdgsSeen(*ri, CLeptonPDGs); Nprotons_seen = CountNPdgsSeen(*ri, ProtonPDG); Nneutrons_seen = CountNPdgsSeen(*ri, NeutronPDG); Ncpip_seen = CountNPdgsSeen(*ri, cpipPDG); Ncpim_seen = CountNPdgsSeen(*ri, cpimPDG); Ncpi_seen = Ncpip_seen + Ncpim_seen; Npi0_seen = CountNPdgsSeen(*ri, pi0PDG); NcK_seen = CountNPdgsSeen(*ri, cKPDG); NK0_seen = CountNPdgsSeen(*ri, K0PDG); Nothers_seen = CountNNotPdgsSeen(*ri, ExplicitPDGs); if (FSCLep_seen && (FSLepMom_rec.Mag() > 1E-8)) { EISLep_QE_rec = FitUtils::EnuQErec(FSLepMom_rec.Mag() / 1000.0, FSLepMom_rec.CosTheta(), 34, PDGFSLep_true > 0) * 1000.0; } else { EISLep_QE_rec = 0; } EISLep_LepHad_rec = EFSLep_rec + EFSHad_rec; EISLep_LepHadVis_rec = EFSLep_rec + EFSHad_rec + EFSVis; Nprotons_contributed = CountNPdgsContributed(*ri, ProtonPDG); Nneutrons_contributed = CountNPdgsContributed(*ri, NeutronPDG); Ncpip_contributed = CountNPdgsContributed(*ri, cpipPDG); Ncpim_contributed = CountNPdgsContributed(*ri, cpimPDG); Ncpi_contributed = Ncpip_contributed + Ncpim_contributed; Npi0_contributed = CountNPdgsContributed(*ri, pi0PDG); NcK_contributed = CountNPdgsContributed(*ri, cKPDG); NK0_contributed = CountNPdgsContributed(*ri, K0PDG); Ngamma_contributed = CountNPdgsContributed(*ri, GammaPDG); Nothers_contibuted = CountNNotPdgsContributed(*ri, ExplicitPDGs); Weight = event->RWWeight * event->InputWeight; RWWeight = event->RWWeight; InputWeight = event->InputWeight; FluxWeight = GetFluxHistogram()->GetBinContent( GetFluxHistogram()->FindBin(EISLep_true)) / GetFluxHistogram()->Integral(); EffWeight = ri->Weight; flagCCINC_true = PDGFSLep_true & 1; flagCC0K_true = (NcK_true + NK0_true) == 0; flagCC0Pi_true = flagCCINC_true && flagCC0K_true && ((Ncpi_true + Npi0_true) == 0); flagCC1Pi_true = flagCCINC_true && flagCC0K_true && ((Ncpi_true + Npi0_true) == 1); flagCCINC_rec = FSCLep_seen && flagCCINC_true; flagCC0K_rec = (NcK_seen + NK0_seen) == 0; flagCC0Pi_rec = flagCCINC_rec && flagCC0K_rec && ((Ncpi_seen + Npi0_seen) == 0); flagCC1Pi_rec = flagCCINC_rec && flagCC0K_rec && ((Ncpi_seen + Npi0_seen) == 1); if (OutputSummaryTree) { // Fill the eventVariables Tree eventVariables->Fill(); } if (RecoSmear) { RecoSmear->Fill(EISLep_true / 1000.0, flagCCINC_rec ? EISLep_LepHadVis_rec / 1000.0 : -1, Weight); ETrueDistrib_noweight->Fill(EISLep_true / 1000.0, flagCCINC_true ? Weight : 0); ETrueDistrib->Fill(EISLep_true / 1000.0, flagCCINC_true ? Weight * PredEvtRateWeight : 0); ERecDistrib->Fill(EISLep_LepHadVis_rec / 1000.0, flagCCINC_rec ? Weight * PredEvtRateWeight : 0); } }; //******************************************************************** void Smearceptance_Tester::Write(std::string drawOpt) { //******************************************************************** if (OutputSummaryTree) { // First save the TTree eventVariables->Write(); } // Save Flux and Event Histograms too GetInput()->GetFluxHistogram()->Write(); GetInput()->GetEventHistogram()->Write(); if (!RecoSmear) { return; } TH2D *SmearMatrix_ev = static_cast(RecoSmear->Clone("ELepHadVis_Smear_ev")); for (Int_t trueAxis_it = 1; trueAxis_it < RecoSmear->GetXaxis()->GetNbins() + 1; ++trueAxis_it) { double NEISLep = ETrueDistrib_noweight->GetBinContent(trueAxis_it); for (Int_t recoAxis_it = 1; recoAxis_it < RecoSmear->GetYaxis()->GetNbins() + 1; ++recoAxis_it) { if (NEISLep > std::numeric_limits::epsilon()) { SmearMatrix_ev->SetBinContent( trueAxis_it, recoAxis_it, SmearMatrix_ev->GetBinContent(trueAxis_it, recoAxis_it) / NEISLep); } } } ETrueDistrib_noweight->Write(); ETrueDistrib->Write(); ERecDistrib->Write(); RecoSmear->Write(); SmearMatrix_ev->Write(); TH2D *ResponseMatrix_ev = SmearceptanceUtils::SVDGetInverse(SmearMatrix_ev, SVDTruncation); ResponseMatrix_ev->SetName("ResponseMatrix_ev"); ResponseMatrix_ev->Write(); #ifdef DEBUG_SMEARTESTER TMatrixD SmearMatrix_ev_md = SmearceptanceUtils::GetMatrix(SmearMatrix_ev); TH1D *SmearedEvt = static_cast(ERecDistrib->Clone()); SmearedEvt->SetNameTitle("SmearedEvt", ";Rec E_{#nu}; count"); SmearceptanceUtils::PushTH1ThroughMatrixWithErrors( ETrueDistrib, SmearedEvt, SmearMatrix_ev_md, 5000, false); SmearedEvt->Write(); SmearedEvt->Scale(1, "width"); SmearedEvt->SetName("SmearedEvt_bw"); SmearedEvt->Write(); #endif #ifdef __PROB3PP_ENABLED__ FitWeight *fw = FitBase::GetRW(); if (fw->HasRWEngine(kOSCILLATION)) { OscWeightEngine *oscWE = dynamic_cast(fw->GetRWEngine(kOSCILLATION)); TGraph POsc; POsc.Set(1E4 - 1); double min = ETrueDistrib->GetXaxis()->GetBinLowEdge(1); double step = (ETrueDistrib->GetXaxis()->GetBinUpEdge( ETrueDistrib->GetXaxis()->GetNbins()) - ETrueDistrib->GetXaxis()->GetBinLowEdge(1)) / double(1E4); for (size_t i = 1; i < 1E4; ++i) { double enu = min + i * step; double ow = oscWE->CalcWeight(enu, 14); if (ow != ow) { std::cout << "Bad osc weight for ENu: " << enu << std::endl; } POsc.SetPoint(i - 1, enu, ow); } POsc.Write("POsc", TObject::kOverwrite); } #endif TMatrixD ResponseMatrix_evt_md = SmearceptanceUtils::GetMatrix(ResponseMatrix_ev); TH1D *Unfolded_enu_obs = static_cast(ETrueDistrib->Clone()); Unfolded_enu_obs->SetNameTitle("UnfoldedENu_evt", ";True E_{#nu};count"); SmearceptanceUtils::PushTH1ThroughMatrixWithErrors( ERecDistrib, Unfolded_enu_obs, ResponseMatrix_evt_md, 5000, false); Unfolded_enu_obs->Write(); Unfolded_enu_obs->Scale(1, "width"); Unfolded_enu_obs->SetName("UnfoldedENu_evt_bw"); Unfolded_enu_obs->Write(); ETrueDistrib->Scale(1, "width"); ETrueDistrib->SetName("ELep_rate_bw"); ETrueDistrib->Write(); ERecDistrib->Scale(1, "width"); ERecDistrib->SetName("ELepRec_rate_bw"); ERecDistrib->Write(); } // ------------------------------------------------------------------- // Purely MC Plot // Following functions are just overrides to handle this // ------------------------------------------------------------------- //******************************************************************** /// Everything is classed as signal... bool Smearceptance_Tester::isSignal(FitEvent *event) { //******************************************************************** (void)event; return true; }; //******************************************************************** void Smearceptance_Tester::ScaleEvents() { //******************************************************************** // Saving everything to a TTree so no scaling required return; } //******************************************************************** void Smearceptance_Tester::ApplyNormScale(float norm) { //******************************************************************** // Saving everything to a TTree so no scaling required fCurrentNorm = norm; return; } //******************************************************************** void Smearceptance_Tester::FillHistograms() { //******************************************************************** // No Histograms need filling........ return; } //******************************************************************** void Smearceptance_Tester::ResetAll() { //******************************************************************** if (OutputSummaryTree) { eventVariables->Reset(); } return; } //******************************************************************** float Smearceptance_Tester::GetChi2() { //******************************************************************** // No Likelihood to test, purely MC return 0.0; }