diff --git a/src/samples/MCTools/NuisToCAFTruth.cxx b/src/samples/MCTools/NuisToCAFTruth.cxx
index 5ba6935..2377685 100644
--- a/src/samples/MCTools/NuisToCAFTruth.cxx
+++ b/src/samples/MCTools/NuisToCAFTruth.cxx
@@ -1,253 +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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
//********************************************************************
#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<bool>("is_fd", false);
isFHC = instance_sample_configuration.get<bool>("is_numode", true);
isOscSwap = instance_sample_configuration.get<bool>("is_nueswap", false);
// Perform any per-sample configuration in the base class
SimpleMCStudy::Initialize(instance_sample_configuration);
+ outputTree = AddNewTTreeToFile(nuis::persistency::GetOutputFile().get(),
+ "cafTruthTree", write_directory);
+
// 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/ROOTUtility.hxx b/src/utility/ROOTUtility.hxx
index 60551e8..75b0d47 100644
--- a/src/utility/ROOTUtility.hxx
+++ b/src/utility/ROOTUtility.hxx
@@ -1,155 +1,164 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifndef UTILITY_ROOTUTILITY_HXX_SEEN
#define UTILITY_ROOTUTILITY_HXX_SEEN
#include "exception/exception.hxx"
#include "TFile.h"
#include "TTree.h"
#include <iomanip>
#include <iostream>
#include <memory>
#include <string>
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<TFile, decltype(&CloseTFile)>;
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<TTree *>(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) {
+inline TreeFile AddNewTTreeToFile(TFile *f, std::string const &tname,
+ std::string const &dirname = "") {
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());
+ if (dirname.size()) {
+ TDirectory *dir = tf.file->GetDirectory(dirname.c_str());
+ if (!dir) {
+ dir = tf.file->mkdir(dirname.c_str());
+ }
+ tf.tree->SetDirectory(dir);
+ } else {
+ 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