diff --git a/src/event/Particle.hxx b/src/event/Particle.hxx index 04fcfe8..f71d2e9 100644 --- a/src/event/Particle.hxx +++ b/src/event/Particle.hxx @@ -1,78 +1,89 @@ // Copyright 2018 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 . *******************************************************************************/ -#ifndef EVENT_PARTICLE_HXX_SEEN -#define EVENT_PARTICLE_HXX_SEEN +#pragma once #include "event/types.hxx" #include "TLorentzVector.h" namespace nuis { namespace event { class Particle { public: NEW_NUIS_EXCEPT(invalid_particle); #define STATUS_LIST \ X(kNuclearLeaving, 0) \ X(kPrimaryInitialState, 1) \ X(kPrimaryFinalState, 2) \ X(kIntermediate, 3) \ X(kUnknown, 4) \ X(kBlocked, 5) \ X(kNParticleStatus, 6) #define X(A, B) A = B, enum class Status_t { STATUS_LIST }; #undef X Particle(); Particle(Particle const &); PDG_t pdg; TLorentzVector P4; bool operator!() const { return (pdg == std::numeric_limits::max()); } double E() const { return P4.E(); } double KE() const { return P4.E() - P4.M(); } double P() const { return P4.Vect().Mag(); } TVector3 P3() const { return P4.Vect(); } + TVector3 Dir() const { return P3().Unit(); } double M() const { return P4.M(); } double Theta() const { return P4.Vect().Theta(); } double CosTheta() const { return P4.Vect().CosTheta(); } }; } // namespace event } // namespace nuis #define X(A, B) \ case nuis::event::Particle::Status_t::A: { \ return os << #A; \ } inline std::ostream &operator<<(std::ostream &os, nuis::event::Particle::Status_t te) { switch (te) { STATUS_LIST } return os; } #undef X #undef STATUS_LIST -#endif +inline bool operator==(nuis::event::Particle const &l, + nuis::event::Particle const &r) { + if (l.pdg != r.pdg) { + return false; + } + for (size_t i = 0; i < 4; ++i) { + if (l.P4[i] != r.P4[i]) { + return false; + } + } + return true; +} diff --git a/src/samples/MCTools/CMakeLists.txt b/src/samples/MCTools/CMakeLists.txt index 236fdfa..6d1455b 100644 --- a/src/samples/MCTools/CMakeLists.txt +++ b/src/samples/MCTools/CMakeLists.txt @@ -1,16 +1,16 @@ if(USE_NUWRO) LIST(APPEND MC_TOOL_IMPL NuisToNuWro.cxx) endif(USE_NUWRO) -LIST(APPEND MC_TOOL_IMPL VerboseEventSummary.cxx) +LIST(APPEND MC_TOOL_IMPL VerboseEventSummary.cxx NuisToCAFTruth) add_library(MCTools SHARED ${MC_TOOL_IMPL}) target_link_libraries(MCTools ${SAMPLES_MCTOOLS_LIBS}) if(USE_NUWRO) target_compile_options(MCTools PUBLIC ${NUWRO_CXX_FLAGS}) include_directories(${NUWRO_INCLUDE_DIRS}) target_link_libraries(MCTools -Wl,--no-as-needed ${NUWRO_IMPORTED_TARGETS}) endif() install(TARGETS MCTools DESTINATION plugins) diff --git a/src/samples/MCTools/NuisToCAFTruth.cxx b/src/samples/MCTools/NuisToCAFTruth.cxx new file mode 100644 index 0000000..5ba6935 --- /dev/null +++ b/src/samples/MCTools/NuisToCAFTruth.cxx @@ -0,0 +1,253 @@ +// Copyright 2018 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 "samples/SimpleMCStudy.hxx" + +#include "utility/EventTopologyUtility.hxx" +#include "utility/FullEventUtility.hxx" +#include "utility/KinematicUtility.hxx" +#include "utility/PDGCodeUtility.hxx" +#include "utility/ROOTUtility.hxx" + +#include "persistency/ROOTOutput.hxx" + +using namespace nuis::event; +using namespace nuis::utility; + +class NuisToCAFTruth : public SimpleMCStudy { + + TreeFile outputTree; + + // int mode; + + int isFD; + int isFHC; + int isCC; + + int nuPDG; + int nuPDGunosc; + int LepPDG; + + double NuMomX; + double NuMomY; + double NuMomZ; + double Ev; + double LepMomX; + double LepMomY; + double LepMomZ; + double LepE; + double LepNuAngle; + + double Q2; + double W; + double X; + double Y; + + int nP; + int nN; + int nipip; + int nipim; + int nipi0; + int nikp; + int nikm; + int nik0; + int niem; + int niother; + int nNucleus; + int nUNKNOWN; + + double eP; + double eN; + double ePip; + double ePim; + double ePi0; + double eOther; + + double xsecWeight; + + bool isOscSwap; + +public: + NuisToCAFTruth() { ReadGlobalConfigDefaults(); } + + void Initialize(fhicl::ParameterSet const &instance_sample_configuration) { + + outputTree = AddNewTTreeToFile(nuis::persistency::GetOutputFile().get(), + "cafTruthTree"); + + isFD = instance_sample_configuration.get("is_fd", false); + isFHC = instance_sample_configuration.get("is_numode", true); + isOscSwap = instance_sample_configuration.get("is_nueswap", false); + + // Perform any per-sample configuration in the base class + SimpleMCStudy::Initialize(instance_sample_configuration); + + // outputTree->Branch("mode", &mode, "mode/I"); + outputTree->Branch("isFD", &isFD, "isFD/I"); + outputTree->Branch("isFHC", &isFHC, "isFHC/I"); + outputTree->Branch("isCC", &isCC, "isCC/I"); + outputTree->Branch("nuPDG", &nuPDG, "nuPDG/I"); + outputTree->Branch("nuPDGunosc", &nuPDGunosc, "nuPDGunosc/I"); + outputTree->Branch("LepPDG", &LepPDG, "LepPDG/I"); + outputTree->Branch("NuMomX", &NuMomX, "NuMomX/D"); + outputTree->Branch("NuMomY", &NuMomY, "NuMomY/D"); + outputTree->Branch("NuMomZ", &NuMomZ, "NuMomZ/D"); + outputTree->Branch("Ev", &Ev, "Ev/D"); + outputTree->Branch("LepMomX", &LepMomX, "LepMomX/D"); + outputTree->Branch("LepMomY", &LepMomY, "LepMomY/D"); + outputTree->Branch("LepMomZ", &LepMomZ, "LepMomZ/D"); + outputTree->Branch("LepE", &LepE, "LepE/D"); + outputTree->Branch("LepNuAngle", &LepNuAngle, "LepNuAngle/D"); + outputTree->Branch("Q2", &Q2, "Q2/D"); + outputTree->Branch("W", &W, "W/D"); + outputTree->Branch("X", &X, "X/D"); + outputTree->Branch("Y", &Y, "Y/D"); + outputTree->Branch("nP", &nP, "nP/I"); + outputTree->Branch("nN", &nN, "nN/I"); + outputTree->Branch("nipip", &nipip, "nipip/I"); + outputTree->Branch("nipim", &nipim, "nipim/I"); + outputTree->Branch("nipi0", &nipi0, "nipi0/I"); + outputTree->Branch("nikp", &nikp, "nikp/I"); + outputTree->Branch("nikm", &nikm, "nikm/I"); + outputTree->Branch("nik0", &nik0, "nik0/I"); + outputTree->Branch("niem", &niem, "niem/I"); + outputTree->Branch("niother", &niother, "niother/I"); + outputTree->Branch("nNucleus", &nNucleus, "nNucleus/I"); + outputTree->Branch("nUNKNOWN", &nUNKNOWN, "nUNKNOWN/I"); + outputTree->Branch("eP", &eP, "eP/D"); + outputTree->Branch("eN", &eN, "eN/D"); + outputTree->Branch("ePip", &ePip, "ePip/D"); + outputTree->Branch("ePim", &ePim, "ePim/D"); + outputTree->Branch("ePi0", &ePi0, "ePi0/D"); + outputTree->Branch("eOther", &eOther, "eOther/D"); + outputTree->Branch("xsecWeight", &xsecWeight, "xsecWeight/D"); + + //! Here you do whatever you want to do on a per-event basis. + //! This method will be run for every event given by the input handler, + //! the weight allows differential xsecs to be plotted easily and also + //! contains any requested reweightable parameter variations. + ProcessEventFunction = [&](nuis::event::FullEvent const &ev, + double weight) -> void { + xsecWeight = weight; + + nP = 0; + nN = 0; + nipip = 0; + nipim = 0; + nipi0 = 0; + nikp = 0; + nikm = 0; + nik0 = 0; + niem = 0; + niother = 0; + nNucleus = 0; + nUNKNOWN = 0; + eP = 0; + eN = 0; + ePip = 0; + ePim = 0; + ePi0 = 0; + eOther = 0; + + Particle ISLep = GetHMISNeutralLepton(ev); + Particle FSLep = GetHMFSLepton(ev); + + isCC = IsChargedLepton(FSLep.pdg); + + nuPDG = ISLep.pdg; + nuPDGunosc = isOscSwap + ? (ISLep.pdg < 0 ? pdgcodes::kNuMuBar : pdgcodes::kNuMu) + : ISLep.pdg; + LepPDG = FSLep.pdg; + + NuMomX = ISLep.P4.X(); + NuMomY = ISLep.P4.Y(); + NuMomZ = ISLep.P4.Z(); + Ev = ISLep.P4.E(); + + LepMomX = FSLep.P4.X(); + LepMomY = FSLep.P4.Y(); + LepMomZ = FSLep.P4.Z(); + LepE = FSLep.P4.E(); + + LepNuAngle = FSLep.Dir().Dot(ISLep.Dir()); + + TLorentzVector FourMomTransfer = GetEnergyMomentumTransfer(ev); + + Q2 = -FourMomTransfer.Mag2(); + W = sqrt(-Q2 + + 2 * pdgmasses::kNucleonAverageMass_MeV * FourMomTransfer.E() + + pdgmasses::kNucleonAverageMass_MeV * + pdgmasses::kNucleonAverageMass_MeV); + X = Q2 / (2 * pdgmasses::kNucleonAverageMass_MeV * FourMomTransfer.E()); + Y = 1 - LepE / Ev; + + for (Particle const &p : GetNuclearLeavingParticles(ev)) { + if (p == FSLep) { + continue; + } + + if (IsProton(p.pdg)) { + nP++; + eP += p.KE(); + } else if (IsNeutron(p.pdg)) { + nN++; + eN += p.KE(); + } else if (IsPositivePion(p.pdg)) { + nipip++; + ePip += p.KE(); + } else if (IsNegativePion(p.pdg)) { + nipim++; + ePim += p.KE(); + } else if (IsNeutralPion(p.pdg)) { + nipi0++; + ePi0 += p.KE(); + } else if (IsPositiveKaon(p.pdg)) { + nikp++; + eOther += p.KE(); + } else if (IsNegativeKaon(p.pdg)) { + nikm++; + eOther += p.KE(); + } else if (IsNeutralKaon(p.pdg)) { + nik0++; + eOther += p.KE(); + } else if (IsGamma(p.pdg)) { + niem++; + eOther += p.KE(); + } else if (IsNuclearPDG(p.pdg)) { + nNucleus++; + } else { + niother++; + eOther += p.KE(); + } + } + + outputTree->Fill(); + }; + } + + std::string Name() { return "NuisToCAFTruth"; } + + void Write() {} +}; + +//! These declarations allow your class to be loaded dynamically by NUISANCE +DECLARE_PLUGIN(IEventProcessor, NuisToCAFTruth); diff --git a/src/utility/FullEventUtility.cxx b/src/utility/FullEventUtility.cxx index 1406c6b..7046b2e 100644 --- a/src/utility/FullEventUtility.cxx +++ b/src/utility/FullEventUtility.cxx @@ -1,169 +1,172 @@ #include "utility/FullEventUtility.hxx" #include "event/FullEvent.hxx" #include "utility/PDGCodeUtility.hxx" using namespace nuis::event; namespace nuis { namespace utility { event::Particle GetHMParticle(std::vector const &parts) { if (parts.size()) { return event::Particle(); } return *std::max_element( parts.begin(), parts.end(), [](event::Particle const &l, event::Particle const &r) { return l.P() < r.P(); }); } std::vector GetParticles(FullEvent const &ev, std::vector const &pdgs, Particle::Status_t status, bool include_matching_pdg) { std::vector selected_particles; for (auto const &parts : ev.ParticleStack) { if (parts.status != status) { continue; } for (Particle const &part : parts.particles) { bool matched_pdg = false; for (PDG_t pdg : pdgs) { matched_pdg = matched_pdg || (part.pdg == pdg); } bool keep = ((include_matching_pdg && matched_pdg) || (!include_matching_pdg && !matched_pdg)); if (!keep) { continue; } selected_particles.push_back(part); } } return selected_particles; } std::vector const &GetISParticles(FullEvent const &ev) { return ev .ParticleStack[static_cast( Particle::Status_t::kPrimaryInitialState)] .particles; } std::vector const &GetPrimaryFSParticles(FullEvent const &ev) { return ev .ParticleStack[static_cast( Particle::Status_t::kPrimaryFinalState)] .particles; } std::vector const &GetNuclearLeavingParticles(FullEvent const &ev) { return ev .ParticleStack[static_cast(Particle::Status_t::kNuclearLeaving)] .particles; } Particle GetHMParticle(FullEvent const &ev, std::vector const &pdgs, Particle::Status_t status, bool include_matching_pdg) { Particle HMParticle; for (auto const &parts : ev.ParticleStack) { if (parts.status != status) { continue; } for (Particle const &part : parts.particles) { bool matched_pdg = false; for (PDG_t pdg : pdgs) { matched_pdg = matched_pdg || (part.pdg == pdg); } bool keep = ((include_matching_pdg && matched_pdg) || (!include_matching_pdg && !matched_pdg)); if (!keep) { continue; } if (part.P() > HMParticle.P()) { HMParticle = part; } } } return HMParticle; } event::Particle GetHMISParticle(event::FullEvent const &ev, std::vector const &pdgs) { return GetHMParticle(ev, pdgs, Particle::Status_t::kPrimaryInitialState); } event::Particle GetHMFSParticle(event::FullEvent const &ev, std::vector const &pdgs) { return GetHMParticle(ev, pdgs); } std::vector GetFSChargedLeptons(FullEvent const &ev) { return GetParticles(ev, pdgcodes::ChargedLeptons); } std::vector GetFSNeutralLeptons(FullEvent const &ev) { return GetParticles(ev, pdgcodes::NeutralLeptons); } std::vector GetISNeutralLeptons(FullEvent const &ev) { return GetParticles(ev, pdgcodes::NeutralLeptons, Particle::Status_t::kPrimaryInitialState); } std::vector GetFSChargedPions(FullEvent const &ev) { return GetParticles(ev, pdgcodes::ChargedPions); } std::vector GetFSNeutralPions(FullEvent const &ev) { return GetParticles(ev, pdgcodes::NeutralPions); } std::vector GetFSPions(FullEvent const &ev) { return GetParticles(ev, pdgcodes::Pions); } std::vector GetFSProtons(FullEvent const &ev) { return GetParticles(ev, pdgcodes::Protons); } std::vector GetFSNeutrons(FullEvent const &ev) { return GetParticles(ev, pdgcodes::Neutron); } std::vector GetFSNucleons(FullEvent const &ev) { return GetParticles(ev, pdgcodes::Nucleons); } std::vector GetFSOthers(FullEvent const &ev) { return GetParticles(ev, pdgcodes::CommonParticles, Particle::Status_t::kNuclearLeaving, false); } Particle GetHMFSChargedLepton(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::ChargedLeptons); } Particle GetHMFSNeutralLepton(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::NeutralLeptons); } +Particle GetHMFSLepton(event::FullEvent const &ev) { + return GetHMParticle(ev, pdgcodes::Leptons); +} Particle GetHMISNeutralLepton(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::NeutralLeptons, Particle::Status_t::kPrimaryInitialState); } Particle GetHMFSChargedPion(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::ChargedPions); } Particle GetHMFSNeutralPion(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::NeutralPions); } Particle GetHMFSPion(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::Pions); } Particle GetHMFSProton(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::Protons); } Particle GetHMFSNeutron(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::Neutron); } Particle GetHMFSNucleon(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::Nucleons); } Particle GetHMFSOther(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::CommonParticles, Particle::Status_t::kNuclearLeaving, false); } } // namespace utility } // namespace nuis diff --git a/src/utility/FullEventUtility.hxx b/src/utility/FullEventUtility.hxx index aa25bea..22c7ac9 100644 --- a/src/utility/FullEventUtility.hxx +++ b/src/utility/FullEventUtility.hxx @@ -1,84 +1,85 @@ // Copyright 2018 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 . *******************************************************************************/ #pragma once #include "event/types.hxx" #include "event/Particle.hxx" #include namespace nuis { namespace event { class FullEvent; } // namespace event } // namespace nuis namespace nuis { namespace utility { event::Particle GetHMParticle(std::vector); std::vector GetParticles(event::FullEvent const &, std::vector const &, event::Particle::Status_t status = event::Particle::Status_t::kNuclearLeaving, bool include_matching_pdg = true); std::vector const &GetISParticles(event::FullEvent const &); std::vector const & GetPrimaryFSParticles(event::FullEvent const &); std::vector const & GetNuclearLeavingParticles(event::FullEvent const &); event::Particle GetHMParticle(event::FullEvent const &, std::vector const &, event::Particle::Status_t status = event::Particle::Status_t::kNuclearLeaving, bool include_matching_pdg = true); event::Particle GetHMISParticle(event::FullEvent const &, std::vector const &); event::Particle GetHMFSParticle(event::FullEvent const &, std::vector const &); std::vector GetFSChargedLeptons(event::FullEvent const &); std::vector GetFSNeutralLeptons(event::FullEvent const &); std::vector GetISNeutralLeptons(event::FullEvent const &); std::vector GetFSChargedPions(event::FullEvent const &); std::vector GetFSNeutralPions(event::FullEvent const &); std::vector GetFSPions(event::FullEvent const &); std::vector GetFSProtons(event::FullEvent const &); std::vector GetFSNeutrons(event::FullEvent const &); std::vector GetFSNucleons(event::FullEvent const &); std::vector GetFSOthers(event::FullEvent const &); event::Particle GetHMFSChargedLepton(event::FullEvent const &); event::Particle GetHMFSNeutralLepton(event::FullEvent const &); +event::Particle GetHMFSLepton(event::FullEvent const &); event::Particle GetHMISNeutralLepton(event::FullEvent const &); event::Particle GetHMFSChargedPion(event::FullEvent const &); event::Particle GetHMFSNeutralPion(event::FullEvent const &); event::Particle GetHMFSPion(event::FullEvent const &); event::Particle GetHMFSProton(event::FullEvent const &); event::Particle GetHMFSNeutron(event::FullEvent const &); event::Particle GetHMFSNucleon(event::FullEvent const &); event::Particle GetHMFSOther(event::FullEvent const &); } // namespace utility } // namespace nuis diff --git a/src/utility/KinematicUtility.cxx b/src/utility/KinematicUtility.cxx index cfdbd69..be2a305 100644 --- a/src/utility/KinematicUtility.cxx +++ b/src/utility/KinematicUtility.cxx @@ -1,228 +1,228 @@ #include "utility/KinematicUtility.hxx" #include "event/FullEvent.hxx" #include "event/Particle.hxx" #include "utility/FullEventUtility.hxx" #include "utility/PDGCodeUtility.hxx" using namespace nuis::event; namespace nuis { namespace utility { double GetNeutrinoEQERec(FullEvent const &fev, double SeparationEnergy_MeV) { Particle const &neutrino = GetHMISNeutralLepton(fev); if (!neutrino) { throw Particle::invalid_particle() << "[ERROR]: In GetNeutrinoEQERec, expected to be able to get IS " "neutral lepton, but found none:" << "\n" << fev.to_string(); } Particle const &charged_lepton = GetHMFSChargedLepton(fev); if (!charged_lepton) { throw Particle::invalid_particle() << "[ERROR]: In GetNeutrinoEQERec, expected to be able to get FS " "charged lepton, but found none:" << "\n" << fev.to_string(); } double mass_ISN = (IsMatter(neutrino.pdg) ? (pdgmasses::kNeutronMass_MeV - SeparationEnergy_MeV) : (pdgmasses::kProtonMass_MeV - SeparationEnergy_MeV)) / 1000.0; double mass_FSN = (IsMatter(neutrino.pdg) ? pdgmasses::kProtonMass_MeV : pdgmasses::kNeutronMass_MeV) / 1000.0; double el_GeV = charged_lepton.E() / 1000.0; double pl_GeV = charged_lepton.P() / 1000.0; // momentum of lepton double ml_GeV = charged_lepton.M() / 1000.0; // lepton mass double Theta_nu_mu = neutrino.P3().Angle(charged_lepton.P3()); return (2.0 * mass_ISN * el_GeV - ml_GeV * ml_GeV + mass_FSN * mass_FSN - mass_ISN * mass_ISN) / (2.0 * (mass_ISN - el_GeV + pl_GeV * cos(Theta_nu_mu))); } double GetNeutrinoQ2QERec(FullEvent const &fev, double SeparationEnergy_MeV) { Particle const &neutrino = GetHMISNeutralLepton(fev); if (!neutrino) { throw Particle::invalid_particle() << "[ERROR]: In GetNeutrinoQ2QERec, expected to be able to get IS " "neutral lepton, but found none:" << "\n" << fev.to_string(); } Particle const &charged_lepton = GetHMFSChargedLepton(fev); if (!charged_lepton) { throw Particle::invalid_particle() << "[ERROR]: In GetNeutrinoQ2QERec, expected to be able to get FS " "charged lepton, but found none:" << "\n" << fev.to_string(); } double el_GeV = charged_lepton.E() / 1000.0; double pl_GeV = charged_lepton.P() / 1000.0; // momentum of lepton double ml_GeV = charged_lepton.M() / 1000.0; // lepton mass double Theta_nu_mu = neutrino.P3().Angle(charged_lepton.P3()); return -ml_GeV * ml_GeV + 2.0 * GetNeutrinoEQERec(fev, SeparationEnergy_MeV) * (el_GeV - pl_GeV * cos(Theta_nu_mu)); } TVector3 GetVectorInTPlane(const TVector3 &inp, const TVector3 &planarNormal) { TVector3 pnUnit = planarNormal.Unit(); double inpProjectPN = inp.Dot(pnUnit); TVector3 InPlanarInput = inp - (inpProjectPN * pnUnit); return InPlanarInput; } TVector3 GetUnitVectorInTPlane(const TVector3 &inp, const TVector3 &planarNormal) { return GetVectorInTPlane(inp, planarNormal).Unit(); } double GetDeltaPhiT(TVector3 const &V_lepton, TVector3 const &V_other, TVector3 const &Normal, bool PiMinus) { TVector3 V_lepton_T = GetUnitVectorInTPlane(V_lepton, Normal); TVector3 V_other_T = GetUnitVectorInTPlane(V_other, Normal); return PiMinus ? acos(V_lepton_T.Dot(V_other_T)) : (M_PI - acos(V_lepton_T.Dot(V_other_T))); } TVector3 GetDeltaPT(TVector3 const &V_lepton, TVector3 const &V_other, TVector3 const &Normal) { TVector3 V_lepton_T = GetVectorInTPlane(V_lepton, Normal); TVector3 V_other_T = GetVectorInTPlane(V_other, Normal); return V_lepton_T + V_other_T; } double GetDeltaAlphaT(TVector3 const &V_lepton, TVector3 const &V_other, TVector3 const &Normal, bool PiMinus) { TVector3 DeltaPT = GetDeltaPT(V_lepton, V_other, Normal); return GetDeltaPhiT(V_lepton, DeltaPT, Normal, PiMinus); } TVector3 GetDeltaPT_CC0PiN(event::FullEvent const &fev, event::PDG_t islep_pdg, event::PDG_t fslep_pdg, event::PDG_t fsnuc_pdg) { if (fslep_pdg == nuis::utility::pdgcodes::kDefault) { // Assume neutrino CC fslep_pdg = islep_pdg > 0 ? islep_pdg - 1 : islep_pdg + 1; } if (fsnuc_pdg == nuis::utility::pdgcodes::kDefault) { fsnuc_pdg = islep_pdg > 0 ? nuis::utility::pdgcodes::kProton : nuis::utility::pdgcodes::kNeutron; } Particle ISSLep = GetHMISParticle(fev, {islep_pdg}); if (!ISSLep) { return {-std::numeric_limits::max(), -std::numeric_limits::max(), -std::numeric_limits::max()}; } Particle FSLep = GetHMFSParticle(fev, {fslep_pdg}); if (!FSLep) { return {-std::numeric_limits::max(), -std::numeric_limits::max(), -std::numeric_limits::max()}; } Particle FSNuc = GetHMFSParticle(fev, {fsnuc_pdg}); if (!FSNuc) { return {-std::numeric_limits::max(), -std::numeric_limits::max(), -std::numeric_limits::max()}; } return GetDeltaPT(FSLep.P3(), FSNuc.P3(), ISSLep.P3()); } double GetDeltaPhiT_CC0PiN(event::FullEvent const &fev, event::PDG_t islep_pdg, event::PDG_t fslep_pdg, event::PDG_t fsnuc_pdg) { if (fslep_pdg == nuis::utility::pdgcodes::kDefault) { // Assume neutrino CC fslep_pdg = islep_pdg > 0 ? islep_pdg - 1 : islep_pdg + 1; } if (fsnuc_pdg == nuis::utility::pdgcodes::kDefault) { fsnuc_pdg = islep_pdg > 0 ? nuis::utility::pdgcodes::kProton : nuis::utility::pdgcodes::kNeutron; } Particle ISSLep = GetHMISParticle(fev, {islep_pdg}); if (!ISSLep) { return -std::numeric_limits::max(); } Particle FSLep = GetHMFSParticle(fev, {fslep_pdg}); if (!FSLep) { return -std::numeric_limits::max(); } Particle FSNuc = GetHMFSParticle(fev, {fsnuc_pdg}); if (!FSNuc) { return -std::numeric_limits::max(); } return GetDeltaPhiT(FSLep.P3(), FSNuc.P3(), ISSLep.P3()); } double GetDeltaAlphaT_CC0PiN(event::FullEvent const &fev, event::PDG_t islep_pdg, event::PDG_t fslep_pdg, event::PDG_t fsnuc_pdg) { if (fslep_pdg == nuis::utility::pdgcodes::kDefault) { // Assume neutrino CC fslep_pdg = islep_pdg > 0 ? islep_pdg - 1 : islep_pdg + 1; } if (fsnuc_pdg == nuis::utility::pdgcodes::kDefault) { fsnuc_pdg = islep_pdg > 0 ? nuis::utility::pdgcodes::kProton : nuis::utility::pdgcodes::kNeutron; } Particle ISSLep = GetHMISParticle(fev, {islep_pdg}); if (!ISSLep) { return -std::numeric_limits::max(); } Particle FSLep = GetHMFSParticle(fev, {fslep_pdg}); if (!FSLep) { return -std::numeric_limits::max(); } Particle FSNuc = GetHMFSParticle(fev, {fsnuc_pdg}); if (!FSNuc) { return -std::numeric_limits::max(); } return GetDeltaAlphaT(FSLep.P3(), FSNuc.P3(), ISSLep.P3()); } TLorentzVector GetEnergyMomentumTransfer(event::FullEvent const &fev) { Particle const &neutrino = GetHMISNeutralLepton(fev); if (!neutrino) { throw Particle::invalid_particle() << "[ERROR]: In GetEnergyMomentumTransfer, expected to be able to get " "IS neutral lepton, but found none: \n" << fev.to_string(); } - Particle const &charged_lepton = GetHMFSChargedLepton(fev); - if (!charged_lepton) { + Particle const &fs_lepton = GetHMFSLepton(fev); + if (!fs_lepton) { throw Particle::invalid_particle() << "[ERROR]: In GetEnergyMomentumTransfer, expected to be able to get " - "FS charged lepton, but found none: \n" + "FS lepton, but found none: \n" << fev.to_string(); } - return (neutrino.P4 - charged_lepton.P4); + return (neutrino.P4 - fs_lepton.P4); } } // namespace utility } // namespace nuis diff --git a/src/utility/PDGCodeUtility.cxx b/src/utility/PDGCodeUtility.cxx index f777df9..d2d7e4b 100644 --- a/src/utility/PDGCodeUtility.cxx +++ b/src/utility/PDGCodeUtility.cxx @@ -1,118 +1,128 @@ #include "utility/PDGCodeUtility.hxx" #include #include using namespace nuis::event; using namespace nuis::utility::pdgcodes; using namespace nuis::utility::pdgmasses; namespace nuis { namespace utility { static std::map const PDGMasses{ {kNuMu, kNuMuMass_MeV}, {kNuMuBar, kNuMuBarMass_MeV}, {kMu, kMuMass_MeV}, {kMuPlus, kMuPlusMass_MeV}, {kNue, kNueMass_MeV}, {kNueBar, kNueBarMass_MeV}, {kPiPlus, kPiPlusMass_MeV}, {kPiMinus, kPiMinusMass_MeV}, {kPi0, kPi0Mass_MeV}, {kProton, kNeutronMass_MeV}, {kNeutron, kProtonMass_MeV}}; double GetPDGMass(PDG_t pdg) { if (PDGMasses.find(pdg) == PDGMasses.end()) { throw unhandled_pdg_code() << "[ERROR]: Unknown mass for particle with PDG code: " << pdg << ", please add it to src/utility/PDGCodeUtility.cxx:PDGMasses"; } return PDGMasses.at(pdg); } NEW_NUIS_EXCEPT(invalid_MatterType); bool IsInPDGList(PDG_t pdg, std::vector const &MatterList, std::vector const &AntiMatterList, MatterType type) { switch (type) { case kMatter: { return std::count(MatterList.begin(), MatterList.end(), pdg); } case kMatterAntimatter: { return std::count(MatterList.begin(), MatterList.end(), pdg) || std::count(AntiMatterList.begin(), AntiMatterList.end(), pdg); } case kAntimatter: { return std::count(AntiMatterList.begin(), AntiMatterList.end(), pdg); } default: { throw invalid_MatterType(); } } } bool IsNeutralLepton(PDG_t pdg, MatterType type) { return IsInPDGList(pdg, NeutralLeptons_matter, NeutralLeptons_antimatter, type); } bool IsChargedLepton(PDG_t pdg, MatterType type) { return IsInPDGList(pdg, ChargedLeptons_matter, ChargedLeptons_antimatter, type); } bool IsProton(PDG_t pdg, MatterType type) { return IsInPDGList(pdg, Proton_matter, Proton_antimatter, type); } bool IsNeutron(PDG_t pdg) { return pdg == Neutron[0]; } bool IsChargedPion(PDG_t pdg) { return std::count(ChargedPions.begin(), ChargedPions.end(), pdg); } -bool IsNeutralPion(PDG_t pdg) { - return std::count(NeutralPions.begin(), NeutralPions.end(), pdg); -} +bool IsNeutralPion(PDG_t pdg) { return (pdg == kPi0); } +bool IsPositivePion(PDG_t pdg) { return (pdg == kPiPlus); } +bool IsNegativePion(PDG_t pdg) { return (pdg == kPiMinus); } bool IsPion(PDG_t pdg) { return std::count(Pions.begin(), Pions.end(), pdg); } +bool IsChargedKaon(PDG_t pdg) { + return std::count(ChargedKaons.begin(), ChargedKaons.end(), pdg); +} +bool IsNeutralKaon(PDG_t pdg) { + return std::count(NeutralKaons.begin(), NeutralKaons.end(), pdg); +} +bool IsPositiveKaon(PDG_t pdg) { return (pdg == kKPlus); } +bool IsNegativeKaon(PDG_t pdg) { return (pdg == kKMinus); } +bool IsKaon(PDG_t pdg) { return std::count(Kaons.begin(), Kaons.end(), pdg); } +bool IsGamma(PDG_t pdg) { return (pdg == kGamma); } bool IsOther(PDG_t pdg) { return !std::count(CommonParticles.begin(), CommonParticles.end(), pdg); } bool IsMatter(PDG_t pdg) { // Special cases switch (pdg) { case kPiPlus: { return true; } case kPiMinus: { return true; } case kPi0: { return true; } case kNeutron: { return true; } } return (pdg > 0); } bool IsAntiMatter(PDG_t pdg) { // Special cases switch (pdg) { case kPiPlus: { return true; } case kPiMinus: { return true; } case kPi0: { return true; } case kNeutron: { return true; } } return (pdg < 0); } bool IsNuclearPDG(event::PDG_t pdg) { return (pdg > 1000000000); } PDG_t MakeNuclearPDG(size_t A, size_t Z) { return 1000 * Z + 10 * A + 1000000000; } size_t GetA(PDG_t pdg) { return ((pdg / 10) % 1000); } size_t GetZ(PDG_t pdg) { return ((pdg / 1000) % 1000); } } // namespace utility } // namespace nuis diff --git a/src/utility/PDGCodeUtility.hxx b/src/utility/PDGCodeUtility.hxx index 3cf290c..6cf5c7f 100644 --- a/src/utility/PDGCodeUtility.hxx +++ b/src/utility/PDGCodeUtility.hxx @@ -1,133 +1,157 @@ // Copyright 2018 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 . *******************************************************************************/ #pragma once #include "event/types.hxx" #include "exception/exception.hxx" #include namespace nuis { namespace utility { namespace pdgcodes { NEW_NUIS_EXCEPT(unhandled_pdg_code); enum MatterType { kMatter = 1, kMatterAntimatter = 0, kAntimatter = -1 }; static event::PDG_t const kDefault = 0; static event::PDG_t const kNuMu = 14; static event::PDG_t const kNuMuBar = -14; static event::PDG_t const kMu = 13; static event::PDG_t const kMuPlus = -13; static event::PDG_t const kElectron = 11; static event::PDG_t const kPositron = -11; static event::PDG_t const kNue = 12; static event::PDG_t const kNueBar = 12; static event::PDG_t const kPiPlus = 211; static event::PDG_t const kPiMinus = -211; static event::PDG_t const kPi0 = 111; +static event::PDG_t const kK0Long = 130; +static event::PDG_t const kK0Short = 130; +static event::PDG_t const kK0 = 311; +static event::PDG_t const kKPlus = 321; +static event::PDG_t const kKMinus = -321; + static event::PDG_t const kProton = 2212; static event::PDG_t const kNeutron = 2112; +static event::PDG_t const kGamma = 22; + static std::vector const ChargedLeptons{kElectron, kMu, 15, kPositron, kMuPlus, -15}; static std::vector const ChargedLeptons_matter{kElectron, kMu, 15}; static std::vector const ChargedLeptons_antimatter{kPositron, kMuPlus, -15}; static std::vector const NeutralLeptons{kNue, kNuMu, 16, kNueBar, kNuMuBar, -16}; static std::vector const NeutralLeptons_matter{kNue, kNuMu, 16}; static std::vector const NeutralLeptons_antimatter{kNueBar, kNuMuBar, -16}; +static std::vector const Leptons{ + kElectron, kMu, 15, kPositron, kMuPlus, -15, + kNue, kNuMu, 16, kNueBar, kNuMuBar, -16}; static std::vector const ChargedPions{kPiPlus, kPiMinus}; static std::vector const NeutralPions{kPi0}; static std::vector const Pions{kPiPlus, kPiMinus, kPi0}; +static std::vector const ChargedKaons{kKPlus, kKMinus}; +static std::vector const NeutralKaons{kK0, kK0Long, kK0Short}; +static std::vector const Kaons{kKPlus, kKMinus, kK0, kK0Long, + kK0Short}; static std::vector const Protons{kProton, -kProton}; static std::vector const Proton_matter{kProton}; static std::vector const Proton_antimatter{-kProton}; static std::vector const Neutron{kNeutron}; static std::vector const Nucleons{kProton, kNeutron, -kProton}; static std::vector const Nucleons_matter{kProton, kNeutron}; static std::vector const Nucleons_antimatter{-kProton, kNeutron}; static std::vector const CommonParticles{ kElectron, kMu, 15, kPositron, kMuPlus, -15, kNue, kNuMu, 16, kNueBar, kNuMuBar, -16, kPiPlus, kPiMinus, kPi0, kProton, kNeutron}; } // namespace pdgcodes namespace pdgmasses { static double const kNuMuMass_MeV = 0; static double const kNuMuBarMass_MeV = 0; static double const kMuMass_MeV = 105.65; static double const kMuPlusMass_MeV = 105.65; static double const kNueMass_MeV = 0; static double const kNueBarMass_MeV = 0; static double const kPiPlusMass_MeV = 139.57; static double const kPiMinusMass_MeV = 139.57; static double const kPi0Mass_MeV = 134.97; static double const kNeutronMass_MeV = 939.56; static double const kProtonMass_MeV = 938.27; +static double const kNucleonAverageMass_MeV = 938.27; } // namespace pdgmasses double GetPDGMass(event::PDG_t); bool IsInPDGList(event::PDG_t pdg, std::vector const &MatterList, std::vector const &AntiMatterList, pdgcodes::MatterType type = pdgcodes::kMatterAntimatter); bool IsNeutralLepton(event::PDG_t pdg, pdgcodes::MatterType type = pdgcodes::kMatterAntimatter); bool IsChargedLepton(event::PDG_t pdg, pdgcodes::MatterType type = pdgcodes::kMatterAntimatter); bool IsProton(event::PDG_t pdg, pdgcodes::MatterType type = pdgcodes::kMatterAntimatter); bool IsNeutron(event::PDG_t pdg); bool IsChargedPion(event::PDG_t pdg); +bool IsPositivePion(event::PDG_t pdg); +bool IsNegativePion(event::PDG_t pdg); bool IsNeutralPion(event::PDG_t pdg); bool IsPion(event::PDG_t pdg); +bool IsChargedKaon(event::PDG_t pdg); +bool IsPositiveKaon(event::PDG_t pdg); +bool IsNegativeKaon(event::PDG_t pdg); +bool IsNeutralKaon(event::PDG_t pdg); +bool IsKaon(event::PDG_t pdg); +bool IsGamma(event::PDG_t pdg); bool IsOther(event::PDG_t pdg); bool IsMatter(event::PDG_t pdg); bool IsAntiMatter(event::PDG_t pdg); bool IsNuclearPDG(event::PDG_t pdg); event::PDG_t MakeNuclearPDG(size_t A, size_t Z); size_t GetA(event::PDG_t); size_t GetZ(event::PDG_t); } // namespace utility } // namespace nuis diff --git a/src/utility/ROOTUtility.hxx b/src/utility/ROOTUtility.hxx index bc62623..60551e8 100644 --- a/src/utility/ROOTUtility.hxx +++ b/src/utility/ROOTUtility.hxx @@ -1,142 +1,155 @@ // Copyright 2018 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 . *******************************************************************************/ #ifndef UTILITY_ROOTUTILITY_HXX_SEEN #define UTILITY_ROOTUTILITY_HXX_SEEN #include "exception/exception.hxx" #include "TFile.h" #include "TTree.h" #include #include #include #include namespace nuis { namespace utility { NEW_NUIS_EXCEPT(failed_to_open_TFile); NEW_NUIS_EXCEPT(failed_to_get_TTree); inline void CloseTFile(TFile *f = nullptr) { if (f) { std::cout << "[INFO]: Shutting TFile: " << f->GetName() << ", " << f << std::endl; f->Close(); } delete f; } using TFile_ptr = std::unique_ptr; inline TFile_ptr make_unique_TFile(TFile *f) { return TFile_ptr(f, &CloseTFile); } inline TFile_ptr CheckOpenTFile(std::string const &fname, char const *opts = "") { TDirectory *ogDir = gDirectory; TFile *inpF = new TFile(fname.c_str(), opts); if (!inpF || !inpF->IsOpen()) { throw failed_to_open_TFile() << "[ERROR]: Couldn't open input file: " << std::quoted(fname); } if (ogDir) { ogDir->cd(); } return make_unique_TFile(inpF); } struct TreeFile { TFile_ptr file; TTree *tree; bool file_owned; TreeFile() : file(make_unique_TFile(nullptr)), tree(nullptr), file_owned(false) {} TreeFile(TreeFile const &other) = delete; TreeFile &operator=(TreeFile &&other) { file = std::move(other.file); tree = other.tree; file_owned = other.file_owned; other.file = nullptr; other.tree = nullptr; other.file_owned = false; return *this; } TreeFile(TreeFile &&other) : file(std::move(other.file)), tree(other.tree), file_owned(other.file_owned) { other.file = nullptr; other.tree = nullptr; other.file_owned = false; } + + TTree *operator->() { return tree; } + ~TreeFile() { if (!file_owned) { file.release(); } } }; inline TreeFile CheckGetTTree(TFile *file, std::string const &tname) { TreeFile tf; tf.file = make_unique_TFile(file); tf.tree = dynamic_cast(tf.file->Get(tname.c_str())); tf.file_owned = false; if (!tf.tree) { throw failed_to_get_TTree() << "[ERROR]: Failed to get TTree named: " << std::quoted(tname) << " from TFile: " << std::quoted(file->GetName()); } std::cout << "[INFO]: Opened TFile: " << file->GetName() << ", " << file << std::endl; return tf; } inline TreeFile CheckGetTTree(std::string const &fname, std::string const &tname, char const *opts = "") { TreeFile tf = CheckGetTTree(CheckOpenTFile(fname, opts).release(), tname); tf.file_owned = true; return tf; } inline TreeFile MakeNewTTree(std::string const &fname, std::string const &tname, char const *opts = "") { TreeFile tf; tf.file = CheckOpenTFile(fname, opts); tf.tree = new TTree(tname.c_str(), ""); tf.tree->SetDirectory(tf.file.get()); tf.file_owned = true; return tf; } +inline TreeFile AddNewTTreeToFile(TFile *f, std::string const &tname) { + TreeFile tf; + // This is pretty dumb... + tf.file = make_unique_TFile(f); + tf.tree = new TTree(tname.c_str(), ""); + tf.tree->SetDirectory(tf.file.get()); + tf.file_owned = false; + return tf; +} + inline std::string SanitizeROOTObjectName(std::string name) { return name; } } // namespace utility } // namespace nuis #endif