diff --git a/src/MCStudies/Smearceptance_Tester.cxx b/src/MCStudies/Smearceptance_Tester.cxx index 7f6453e..468278e 100644 --- a/src/MCStudies/Smearceptance_Tester.cxx +++ b/src/MCStudies/Smearceptance_Tester.cxx @@ -1,645 +1,646 @@ // 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 <http://www.gnu.org/licenses/>. *******************************************************************************/ #include "Smearceptance_Tester.h" #include "SmearceptanceUtils.h" #include "Smearcepterton.h" //******************************************************************** /// @brief Class to perform smearceptance MC Studies on a custom measurement Smearceptance_Tester::Smearceptance_Tester(std::string name, std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile) { //******************************************************************** // Measurement Details std::vector<std::string> splitName = GeneralUtils::ParseToStr(name, "_"); size_t firstUS = name.find_first_of("_"); std::string smearceptorName = name.substr(firstUS + 1); fName = name.substr(0, firstUS); eventVariables = NULL; // Define our energy range for flux calcs EnuMin = 0.; EnuMax = 100.; // Arbritrarily high energy limit // Set default fitter flags fIsDiag = true; fIsShape = false; fIsRawEvents = false; // This function will sort out the input files automatically and parse all the // inputs,flags,etc. // There may be complex cases where you have to do this by hand, but usually // this will do. Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile); eventVariables = NULL; // Setup fDataHist as a placeholder this->fDataHist = new TH1D(("empty_data"), ("empty-data"), 1, 0, 1); this->SetupDefaultHist(); fFullCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist); covar = StatUtils::GetInvert(fFullCovar); this->fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) / this->TotalIntegratedFlux(); LOG(SAM) << "Smearceptance Flux Scaling Factor = " << fScaleFactor << std::endl; if (fScaleFactor <= 0.0) { ERR(WRN) << "SCALE FACTOR TOO LOW " << std::endl; sleep(20); } // Setup our TTrees this->AddEventVariablesToTree(); smearceptor = &Smearcepterton::Get().GetSmearcepter(smearceptorName); Int_t RecNBins = 20, TrueNBins = 20; double RecBinL = 0, TrueBinL = 0, RecBinH = 10, TrueBinH = 10; if (Config::Get().GetConfigNode("smear.reconstructed.binning")) { std::vector<std::string> args = GeneralUtils::ParseToStr( Config::Get().ConfS("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::Get().GetConfigNode("smear.true.binning")) { std::vector<std::string> args = GeneralUtils::ParseToStr( Config::Get().ConfS("smear.true.binning"), ","); TrueNBins = GeneralUtils::StrToInt(args[0]); TrueBinL = GeneralUtils::StrToDbl(args[1]); TrueBinH = GeneralUtils::StrToDbl(args[2]); } SVDTruncation = 0; if (Config::Get().GetConfigNode("smear.true.binning")) { SVDTruncation = Config::Get().ConfI("smear.SVD.truncation"); QLOG(SAM, "Applying SVD truncation of: " << SVDTruncation) } QLOG(SAM, "Using binning True: " << TrueNBins << ", [" << TrueBinL << " -- " << TrueBinH << "], Rec: " << RecNBins << ", [" << RecBinL << " -- " << RecBinH << "]"); ETrueDistrib = new TH1D("ELep_rate", ";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", ";Recon. E_{#nu};True E_{#nu}", RecNBins, RecBinL, RecBinH, TrueNBins, TrueBinL, TrueBinH); } void Smearceptance_Tester::AddEventVariablesToTree() { // Setup the TTree to save everything if (!eventVariables) { FitPar::Config().out->cd(); eventVariables = new TTree((this->fName + "_VARS").c_str(), (this->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("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_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("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_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_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("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("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"); eventVariables->Branch("xsecScaling", &xsecScaling, "xsecScaling/F"); eventVariables->Branch("flagCCINC_true", &flagCCINC_true, "flagCCINC_true/O"); eventVariables->Branch("flagCC0Pi_true", &flagCC0Pi_true, "flagCC0Pi_true/O"); eventVariables->Branch("flagCCINC_rec", &flagCCINC_rec, "flagCCINC_rec/O"); eventVariables->Branch("flagCC0Pi_rec", &flagCC0Pi_rec, "flagCC0Pi_rec/O"); #ifdef DEBUG_SMEARTESTER eventVariables->Branch("FSMuon_True", &FSMuon_True); eventVariables->Branch("FSMuon_Smeared", &FSMuon_Smeared); #endif } template <size_t N> int CountNPdgsSeen(RecoInfo ri, int (&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 <size_t N> int CountNNotPdgsSeen(RecoInfo ri, int (&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]) ? 0 : 1); } return sum; } template <size_t N> int CountNPdgsContributed(RecoInfo ri, int (&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 <size_t N> int CountNNotPdgsContributed(RecoInfo ri, int (&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]) ? 0 : 1); } return sum; } TVector3 GetHMFSRecParticles(RecoInfo ri, int pdg) { TVector3 mom(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 = ri.RecObjMom[p_it]; } } return mom; } template <size_t N> double SumKE_RecoInfo(RecoInfo ri, int (&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 <size_t N> double SumTE_RecoInfo(RecoInfo ri, int (&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 <size_t N> double SumVisE_RecoInfo(RecoInfo ri, int (&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 <size_t N> double SumVisE_RecoInfo_NotPdgs(RecoInfo ri, int (&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) { //******************************************************************** int cpipPDG[] = {211}; int cpimPDG[] = {-211}; int pi0PDG[] = {111}; int ProtonPDG[] = {2212}; int NeutronPDG[] = {2112}; int GammaPDG[] = {22}; int CLeptonPDGs[] = {11, 13, 15}; int ExplicitPDGs[] = {211, -211, 11, 2212, 2112, 22, 11, 13, 15, 12, 14, 16}; RecoInfo *ri = smearceptor->Smearcept(event); #ifdef DEBUG_SMEARTESTER FSMuon_True = TVector3(0, 0, 0); FSMuon_Smeared = TVector3(0, 0, 0); FitParticle *fsMu = event->GetHMFSMuon(); if (fsMu) { FSMuon_True = fsMu->P3(); FSMuon_Smeared = GetHMFSRecParticles(*ri, 13); } #endif 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_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; 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(); 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_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; TVector3 FSLepMom_rec(0, 0, 0); if (event->GetHMFSAnyLeptons()) { double massLep = event->GetHMFSAnyLeptons()->M(); FSLepMom_rec = GetHMFSRecParticles(*ri, event->GetHMFSAnyLeptons()->PDG()); EFSLep_rec = (FSLepMom_rec.Mag() > 1E-5) ? sqrt(FSLepMom_rec * FSLepMom_rec + massLep * massLep) : 0; } 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_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_cpip + EFSVis_cpim + EFSVis_cpi + EFSVis_pi0 + EFSVis_p + EFSVis_n + EFSVis_gamma; 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); 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); 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; xsecScaling = fScaleFactor; flagCCINC_true = PDGFSLep_true & 1; flagCC0Pi_true = (Ncpi_true + Npi0_true) == 0; flagCCINC_rec = FSCLep_seen && PDGFSLep_true & 1; flagCC0Pi_rec = ((Ncpi_seen + Npi0_seen) == 0) && flagCCINC_rec; // Fill the eventVariables Tree eventVariables->Fill(); RecoSmear->Fill(flagCCINC_rec ? EISLep_LepHadVis_rec / 1000.0 : -1, EISLep_true / 1000.0, Weight); ETrueDistrib->Fill(EISLep_true / 1000.0, flagCCINC_true ? Weight : 0); ERecDistrib->Fill(EISLep_LepHadVis_rec / 1000.0, flagCCINC_rec ? Weight : 0); return; }; //******************************************************************** void Smearceptance_Tester::Write(std::string drawOpt) { //******************************************************************** // First save the TTree eventVariables->Write(); // Save Flux and Event Histograms too GetInput()->GetFluxHistogram()->Write(); GetInput()->GetEventHistogram()->Write(); TH2D *SmearMatrix_ev = static_cast<TH2D *>(RecoSmear->Clone("ELepHadVis_Smear_ev")); for (Int_t trueAxis_it = 1; trueAxis_it < RecoSmear->GetYaxis()->GetNbins() + 1; ++trueAxis_it) { double NEISLep = ETrueDistrib->GetBinContent(trueAxis_it); for (Int_t recoAxis_it = 1; recoAxis_it < RecoSmear->GetXaxis()->GetNbins() + 1; ++recoAxis_it) { if (NEISLep > std::numeric_limits<double>::epsilon()) { SmearMatrix_ev->SetBinContent( recoAxis_it, trueAxis_it, SmearMatrix_ev->GetBinContent(recoAxis_it, trueAxis_it) / NEISLep); } } } ETrueDistrib->Write(); ERecDistrib->Write(); RecoSmear->Write(); SmearMatrix_ev->Write(); TH2D *ResponseMatrix_ev = SmearceptanceUtils::SVDGetInverse(SmearMatrix_ev, SVDTruncation); ResponseMatrix_ev = SmearceptanceUtils::SwapXYTH2D(ResponseMatrix_ev); ResponseMatrix_ev->SetName("ResponseMatrix_ev"); ResponseMatrix_ev->Write(); #ifdef DEBUG_SMEARTESTER TMatrixD SmearMatrix_ev_md = SmearceptanceUtils::GetMatrix( SmearceptanceUtils::SwapXYTH2D(SmearMatrix_ev)); TH1D *SmearedEvt = static_cast<TH1D *>(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 TMatrixD ResponseMatrix_evt_md = SmearceptanceUtils::GetMatrix(ResponseMatrix_ev); TH1D *Unfolded_enu_obs = static_cast<TH1D *>(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 this->fCurrentNorm = norm; return; } //******************************************************************** void Smearceptance_Tester::FillHistograms() { //******************************************************************** // No Histograms need filling........ return; } //******************************************************************** void Smearceptance_Tester::ResetAll() { //******************************************************************** eventVariables->Reset(); return; } //******************************************************************** float Smearceptance_Tester::GetChi2() { //******************************************************************** // No Likelihood to test, purely MC return 0.0; } diff --git a/src/Smearceptance/EfficiencyApplicator.cxx b/src/Smearceptance/EfficiencyApplicator.cxx index 79cc736..914b64f 100644 --- a/src/Smearceptance/EfficiencyApplicator.cxx +++ b/src/Smearceptance/EfficiencyApplicator.cxx @@ -1,305 +1,310 @@ // 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 <http://www.gnu.org/licenses/>. *******************************************************************************/ #include "EfficiencyApplicator.h" #include "TEfficiency.h" #include "TH2.h" #include "TH3.h" -// #define DEBUG_EFFAPP +//#define DEBUG_EFFAPP EfficiencyApplicator::DependVar GetVarType(std::string const &axisvar) { if (axisvar == "kMomentum") { return EfficiencyApplicator::kMomentum; } else if (axisvar == "kKE") { return EfficiencyApplicator::kKE; } else if (axisvar == "kTheta") { return EfficiencyApplicator::kTheta; } else if (axisvar == "kCosTheta") { return EfficiencyApplicator::kCosTheta; } else if (axisvar == "kPhi") { return EfficiencyApplicator::kPhi; } return EfficiencyApplicator::kNoAxis; } TH1 *GetEffHist(TFile *inputFile, std::string const &HistName) { TH1 *hist = dynamic_cast<TH1 *>(inputFile->Get(HistName.c_str())); if (hist) { return hist; } TEfficiency *effHist = dynamic_cast<TEfficiency *>(inputFile->Get(HistName.c_str())); if (effHist) { TH1D *numer = dynamic_cast<TH1D *>(effHist->GetCopyPassedHisto()); TH1D *denom = dynamic_cast<TH1D *>(effHist->GetCopyTotalHisto()); if (numer && denom) { numer->Divide(denom); denom->SetDirectory(NULL); delete denom; // Gonna be a memory leak, but I'll get over it numer->SetDirectory(NULL); return numer; } ERROR(FTL, "TEfficiency internal histograms were not TH1Ds."); } THROW("Couldn't get appropriate efficiency object named " << HistName << " from input file " << inputFile->GetName()); } /// Reads particle efficiency nodes /// /// Nodes look like: /// <EfficiencyApplicator Name="D00N_ND_LAr"> /// <EfficiencyCurve PDG="211,-211" InputFile="effs.root" /// HistName="cpion_eff_mom_ctheta" NDims="2" XAxis="kMomentum" /// YAxis="kCosTheta" Interpolate="false" /> /// </EfficiencyApplicator> void EfficiencyApplicator::SpecifcSetup(nuiskey &nk) { rand.~TRandom3(); new (&rand) TRandom3(); std::vector<nuiskey> effDescriptors = nk.GetListOfChildNodes("EfficiencyCurve"); for (size_t t_it = 0; t_it < effDescriptors.size(); ++t_it) { std::string inputFileName = effDescriptors[t_it].GetS("InputFile"); std::string HistName = effDescriptors[t_it].GetS("HistName"); TFile inputFile(inputFileName.c_str()); if (!inputFile.IsOpen()) { THROW("Couldn't open specified input root file: " << inputFileName); } TH1 *inpHist = GetEffHist(&inputFile, HistName); if (!inpHist) { THROW("Couldn't get TH1D named: " << HistName << " from input root file: " << inputFileName); } int NDims = effDescriptors[t_it].GetI("NDims"); if (NDims < 1 || NDims > 3) { THROW("Read NDims attribute as: " << NDims << ", efficiency curve can " "currently have between 1 " "and 3 dimensions."); } EfficiencyApplicator::DependVar XVar = GetVarType(effDescriptors[t_it].GetS("XAxis")); double XAxisScale = effDescriptors[t_it].Has("XAxisScaleToMeV") ? effDescriptors[t_it].GetD("XAxisScaleToMeV") : 1; EfficiencyApplicator::DependVar YVar = NDims > 1 ? GetVarType(effDescriptors[t_it].GetS("YAxis")) : EfficiencyApplicator::kNoAxis; double YAxisScale = effDescriptors[t_it].Has("YAxisScaleToMeV") ? effDescriptors[t_it].GetD("YAxisScaleToMeV") : 1; EfficiencyApplicator::DependVar ZVar = NDims > 2 ? GetVarType(effDescriptors[t_it].GetS("ZAxis")) : EfficiencyApplicator::kNoAxis; double ZAxisScale = effDescriptors[t_it].Has("ZAxisScaleToMeV") ? effDescriptors[t_it].GetD("ZAxisScaleToMeV") : 1; bool Interpolate = effDescriptors[t_it].Has("Interpolate") && effDescriptors[t_it].GetI("Interpolate"); bool ApplyAsWeight = effDescriptors[t_it].Has("Interpolate") && effDescriptors[t_it].GetI("Interpolate"); std::string pdgs_s = effDescriptors[t_it].GetS("PDG"); std::vector<int> pdgs_i = GeneralUtils::ParseToInt(pdgs_s, ","); for (size_t pdg_it = 0; pdg_it < pdgs_i.size(); ++pdg_it) { if (Efficiencies.count(pdgs_i[pdg_it])) { ERROR(WRN, "Smearceptor " << ElementName << ":" << InstanceName << " already has a efficiency for PDG: " << pdgs_i[pdg_it]); } EffMap em; em.EffCurve = static_cast<TH1 *>(inpHist->Clone()); em.EffCurve->SetDirectory(NULL); em.Interpolate = Interpolate; // em.ApplyAsWeight = ApplyAsWeight; em.NDims = NDims; em.DependVars[0] = XVar; em.DependVars[1] = YVar; em.DependVars[2] = ZVar; em.AxisScales[0] = XAxisScale; em.AxisScales[1] = YAxisScale; em.AxisScales[2] = ZAxisScale; Efficiencies[pdgs_i[pdg_it]] = em; QLOG(SAM, "Added reconstruction efficiency curve for PDG: " << pdgs_i[pdg_it]); } } } RecoInfo *EfficiencyApplicator::Smearcept(FitEvent *fe) { RecoInfo *ri = new RecoInfo(); for (size_t p_it = 0; p_it < fe->NParticles(); ++p_it) { FitParticle *fp = fe->GetParticle(p_it); #ifdef DEBUG_EFFAPP std::cout << std::endl; std::cout << "[" << p_it << "]: " << fp->PDG() << ", " << fp->Status() << ", " << fp->E() << " -- KE:" << fp->KE() << " Mom: " << fp->P3().Mag() << std::flush; #endif if (fp->Status() != kFinalState) { #ifdef DEBUG_EFFAPP std::cout << " -- Not final state." << std::flush; #endif continue; } if (!Efficiencies.count(fp->PDG())) { #ifdef DEBUG_EFFAPP std::cout << " -- Undetectable." << std::flush; #endif continue; } EffMap &em = Efficiencies[fp->PDG()]; double kineProps[3]; for (Int_t dim_it = 0; dim_it < em.NDims; ++dim_it) { switch (em.DependVars[dim_it]) { case kMomentum: { kineProps[dim_it] = fp->P3().Mag(); break; } case kKE: { kineProps[dim_it] = fp->KE(); break; } case kTheta: { kineProps[dim_it] = fp->P3().Theta(); break; } case kCosTheta: { kineProps[dim_it] = fp->P3().CosTheta(); break; } case kPhi: { kineProps[dim_it] = fp->P3().Phi(); break; } default: { THROW("Trying to find particle value for a kNoAxis."); } } kineProps[dim_it] /= em.AxisScales[dim_it]; } double effProb = 0; switch (em.NDims) { case 1: { TH1 *hist = em.EffCurve; if (em.Interpolate && (hist->GetXaxis()->GetBinCenter(1) < kineProps[0]) && (hist->GetXaxis()->GetBinCenter(hist->GetXaxis()->GetNbins()) > kineProps[0])) { effProb = hist->Interpolate(kineProps[0]); } else { Int_t xbin = hist->GetXaxis()->FindFixBin(kineProps[0]); effProb = hist->GetBinContent(xbin); } break; } case 2: { TH2 *hist = static_cast<TH2 *>(em.EffCurve); if (em.Interpolate && (hist->GetXaxis()->GetBinCenter(1) < kineProps[0]) && (hist->GetXaxis()->GetBinCenter(hist->GetXaxis()->GetNbins()) > kineProps[0]) && (hist->GetYaxis()->GetBinCenter(1) < kineProps[1]) && (hist->GetYaxis()->GetBinCenter(hist->GetYaxis()->GetNbins()) > kineProps[1])) { effProb = hist->Interpolate(kineProps[0], kineProps[1]); } else { Int_t xbin = hist->GetXaxis()->FindFixBin(kineProps[0]); Int_t ybin = hist->GetYaxis()->FindFixBin(kineProps[1]); effProb = hist->GetBinContent(xbin, ybin); #ifdef DEBUG_EFFAPP std::cout << "\t\t: XProp: " << kineProps[0] << ", YProp: " << kineProps[1] << " x/y bins: " << xbin << "/" << ybin << ". Prop ? " << effProb << std::endl; #endif } break; } case 3: { TH3 *hist = static_cast<TH3 *>(em.EffCurve); if (em.Interpolate && (hist->GetXaxis()->GetBinCenter(1) < kineProps[0]) && (hist->GetXaxis()->GetBinCenter(hist->GetXaxis()->GetNbins()) > kineProps[0]) && (hist->GetYaxis()->GetBinCenter(1) < kineProps[1]) && (hist->GetYaxis()->GetBinCenter(hist->GetYaxis()->GetNbins()) > kineProps[2]) && (hist->GetZaxis()->GetBinCenter(hist->GetZaxis()->GetNbins()) > kineProps[2])) { effProb = hist->Interpolate(kineProps[0], kineProps[1], kineProps[2]); } else { Int_t xbin = hist->GetXaxis()->FindFixBin(kineProps[0]); Int_t ybin = hist->GetYaxis()->FindFixBin(kineProps[1]); Int_t zbin = hist->GetZaxis()->FindFixBin(kineProps[2]); effProb = hist->GetBinContent(xbin, ybin, zbin); } break; } } bool accepted = (rand.Uniform() < effProb); if (accepted) { #ifdef DEBUG_EFFAPP std::cout << " -- Reconstructed with probability: " << effProb << std::flush; #endif ri->RecObjMom.push_back(fp->P3()); ri->RecObjClass.push_back(fp->PDG()); continue; } #ifdef DEBUG_EFFAPP std::cout << " -- Rejected with probability: " << effProb << std::flush; #endif } #ifdef DEBUG_EFFAPP std::cout << std::endl; #endif + +#ifdef DEBUG_EFFAPP + std::cout << "Reconstructed " << ri->RecObjMom.size() << " particles. " + << std::endl; +#endif return ri; } diff --git a/src/Smearceptance/Smearcepterton.cxx b/src/Smearceptance/Smearcepterton.cxx index 7f178b2..5382a50 100644 --- a/src/Smearceptance/Smearcepterton.cxx +++ b/src/Smearceptance/Smearcepterton.cxx @@ -1,75 +1,75 @@ // 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 <http://www.gnu.org/licenses/>. *******************************************************************************/ #include "Smearcepterton.h" #include "ThresholdAccepter.h" #include "EfficiencyApplicator.h" #include "GaussianSmearer.h" #include "TrackedMomentumMatrixSmearer.h" #include "MetaSimpleSmearcepter.h" #include <vector> 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<std::string, SmearceptionFactory_fcn> factories; factories["ThresholdAccepter"] = &BuildSmearcepter<ThresholdAccepter>; factories["EfficiencyApplicator"] = &BuildSmearcepter<EfficiencyApplicator>; factories["GaussianSmearer"] = &BuildSmearcepter<GaussianSmearer>; factories["TrackedMomentumMatrixSmearer"] = &BuildSmearcepter<TrackedMomentumMatrixSmearer>; factories["MetaSimpleSmearcepter"] = &BuildSmearcepter<MetaSimpleSmearcepter>; std::vector<nuiskey> smearcepterBlocks = Config::QueryKeys("smearcepters"); for (size_t smearB_it = 0; smearB_it < smearcepterBlocks.size(); ++smearB_it) { std::vector<nuiskey> smearcepters = smearcepterBlocks[smearB_it].GetListOfChildNodes(); - for (size_t smear_it = 0; smear_it < smearcepterBlocks.size(); ++smear_it) { + for (size_t smear_it = 0; smear_it < smearcepters.size(); ++smear_it) { std::string const& smearType = smearcepters[smear_it].GetElementName(); if (!factories.count(smearType)) { ERROR(WRN, "No known smearer accepts elements named: \"" << smearType << "\""); continue; } ISmearcepter* smearer = factories[smearType](smearcepters[smear_it]); Smearcepters[smearer->GetName()] = smearer; QLOG(FIT, "Configured smearer named: " << smearer->GetName() << " of type: " << smearer->GetElementName()); } } }