Page MenuHomeHEPForge
Files F9501143

No OneTemporary
Actions

View Options

diff --git a/src/InputHandler/NuWroInputHandler.cxx b/src/InputHandler/NuWroInputHandler.cxx
index e4446b1..1c4206b 100644
--- a/src/InputHandler/NuWroInputHandler.cxx
+++ b/src/InputHandler/NuWroInputHandler.cxx
@@ -1,465 +1,465 @@
#ifdef __NUWRO_ENABLED__
#include "NuWroInputHandler.h"
#include "InputUtils.h"
NuWroGeneratorInfo::~NuWroGeneratorInfo() { delete fNuWroParticlePDGs; }
void NuWroGeneratorInfo::AddBranchesToTree(TTree *tn) {
tn->Branch("NuWroParticlePDGs", &fNuWroParticlePDGs, "NuWroParticlePDGs/I");
}
void NuWroGeneratorInfo::SetBranchesFromTree(TTree *tn) {
tn->SetBranchAddress("NuWroParticlePDGs", &fNuWroParticlePDGs);
}
void NuWroGeneratorInfo::AllocateParticleStack(int stacksize) {
fNuWroParticlePDGs = new int[stacksize];
}
void NuWroGeneratorInfo::DeallocateParticleStack() {
delete fNuWroParticlePDGs;
}
void NuWroGeneratorInfo::FillGeneratorInfo(event *e) { Reset(); }
void NuWroGeneratorInfo::Reset() {
for (int i = 0; i < kMaxParticles; i++) {
fNuWroParticlePDGs[i] = 0;
}
}
int event1_nof(event *e, int pdg) {
int c = 0;
for (size_t i = 0; i < e->out.size(); i++)
if (e->out[i].pdg == pdg)
c++;
return c;
}
NuWroInputHandler::NuWroInputHandler(std::string const &handle,
std::string const &rawinputs) {
NUIS_LOG(SAM, "Creating NuWroInputHandler : " << handle);
// Run a joint input handling
fName = handle;
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
fSaveExtra = false; // FitPar::Config().GetParB("NuWroSaveExtra");
// Setup the TChain
fNuWroTree = new TChain("treeout");
// Loop over all inputs and grab flux, eventhist, and nevents
std::vector<std::string> inputs = InputUtils::ParseInputFileList(rawinputs);
for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) {
// Open File for histogram access
TFile *inp_file = new TFile(inputs[inp_it].c_str(), "READ");
if (!inp_file or inp_file->IsZombie()) {
NUIS_ABORT("nuwro File IsZombie() at " << inputs[inp_it]);
}
// Get Flux/Event hist
TH1D *fluxhist = (TH1D *)inp_file->Get(
(PlotUtils::GetObjectWithName(inp_file, "FluxHist")).c_str());
TH1D *eventhist = (TH1D *)inp_file->Get(
(PlotUtils::GetObjectWithName(inp_file, "EvtHist")).c_str());
if (!fluxhist or !eventhist) {
NUIS_ERR(FTL, "nuwro FILE doesn't contain flux/xsec info");
if (FitPar::Config().GetParB("regennuwro")) {
NUIS_ERR(FTL,
"Regen NuWro has not been added yet. Email the developers!");
// ProcessNuWroInputFlux(inputs[inp_it]);
throw;
} else {
NUIS_ABORT("If you would like NUISANCE to generate these for you "
<< "please set parameter regennuwro=1 and re-run.");
}
}
// Get N Events
TTree *nuwrotree = (TTree *)inp_file->Get("treeout");
if (!nuwrotree) {
NUIS_ABORT("treeout not located in nuwro file! " << inputs[inp_it]);
}
int nevents = nuwrotree->GetEntries();
// Register input to form flux/event rate hists
RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist);
// Add to TChain
fNuWroTree->Add(inputs[inp_it].c_str());
}
// Registor all our file inputs
SetupJointInputs();
// Setup Events
fNuWroEvent = NULL;
fNuWroTree->SetBranchAddress("e", &fNuWroEvent);
fNuWroTree->GetEntry(0);
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->fType = kNUWRO;
fNUISANCEEvent->fNuwroEvent = fNuWroEvent;
fNUISANCEEvent->HardReset();
if (fSaveExtra) {
fNuWroInfo = new NuWroGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fNuWroInfo);
}
};
NuWroInputHandler::~NuWroInputHandler() {
if (fNuWroTree)
delete fNuWroTree;
}
void NuWroInputHandler::CreateCache() {
// fNuWroTree->SetCacheEntryRange(0, fNEvents);
// fNuWroTree->AddBranchToCache("*", 1);
// fNuWroTree->SetCacheSize(fCacheSize);
}
void NuWroInputHandler::RemoveCache() {
// fNuWroTree->SetCacheEntryRange(0, fNEvents);
// fNuWroTree->AddBranchToCache("*", 0);
// fNuWroTree->SetCacheSize(0);
}
void NuWroInputHandler::ProcessNuWroInputFlux(const std::string file) {}
FitEvent *NuWroInputHandler::GetNuisanceEvent(const UInt_t ent,
const bool lightweight) {
UInt_t entry = ent + fSkip;
// Catch too large entries
if (entry >= (UInt_t)fNEvents)
return NULL;
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fNuWroTree->GetEntry(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
#ifdef __PROB3PP_ENABLED__
for (size_t i = 0; i < fNUISANCEEvent->fNuwroEvent->in.size(); i++) {
if (std::count(PhysConst::pdg_neutrinos, PhysConst::pdg_neutrinos + 4,
fNUISANCEEvent->fNuwroEvent->in[i].pdg)) {
fNUISANCEEvent->probe_E = fNUISANCEEvent->fNuwroEvent->in[i].t;
fNUISANCEEvent->probe_pdg = fNUISANCEEvent->fNuwroEvent->in[i].pdg;
break;
}
}
#endif
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
#ifdef __USE_NUWRO_SRW_EVENTS__
if (!rwEvs.size()) {
fNuwroParams = fNuWroEvent->par;
}
if (entry >= rwEvs.size()) {
rwEvs.push_back(BaseFitEvt());
rwEvs.back().fType = kNUWRO;
rwEvs.back().Mode = fNUISANCEEvent->Mode;
rwEvs.back().fNuwroSRWEvent = SRW::SRWEvent(*fNuWroEvent);
rwEvs.back().fNuwroEvent = NULL;
rwEvs.back().fNuwroParams = &fNuwroParams;
rwEvs.back().probe_E = rwEvs.back().fNuwroSRWEvent.NeutrinoEnergy;
rwEvs.back().probe_pdg = rwEvs.back().fNuwroSRWEvent.NeutrinoPDG;
}
fNUISANCEEvent->fNuwroSRWEvent = SRW::SRWEvent(*fNuWroEvent);
fNUISANCEEvent->fNuwroParams = &fNuwroParams;
fNUISANCEEvent->probe_E = fNUISANCEEvent->fNuwroSRWEvent.NeutrinoEnergy;
fNUISANCEEvent->probe_pdg = fNUISANCEEvent->fNuwroSRWEvent.NeutrinoPDG;
#endif
return fNUISANCEEvent;
}
int NuWroInputHandler::ConvertNuwroMode(event *e) {
Int_t proton_pdg, neutron_pdg, pion_pdg, pion_plus_pdg, pion_minus_pdg,
lambda_pdg, eta_pdg, kaon_pdg, kaon_plus_pdg;
proton_pdg = 2212;
eta_pdg = 221;
neutron_pdg = 2112;
pion_pdg = 111;
pion_plus_pdg = 211;
pion_minus_pdg = -211;
lambda_pdg = 3122;
kaon_pdg = 311;
kaon_plus_pdg = 321;
// Antineutrino modes are *-1
int nu_nubar = 1;
if (e->flag.anty) nu_nubar = -1;
// Quasielastic
if (e->flag.qel){
if (e->flag.cc)
return 1*nu_nubar;
else {
if (event1_nof(e, proton_pdg))
return 51*nu_nubar;
else if (event1_nof(e, neutron_pdg))
return 52*nu_nubar;
}
}
if (e->flag.mec) {
return 2*nu_nubar;
}
// Pion production
if (e->flag.res) {
int npions = event1_nof(e, pion_pdg) + event1_nof(e, pion_plus_pdg) +
event1_nof(e, pion_minus_pdg);
int nkaons= event1_nof(e, kaon_pdg) + event1_nof(e, kaon_plus_pdg);
// Multipion?
- if (npions > 1 || npions == 0) {
+ if ((npions+nkaons) > 1 || npions == 0) {
if (e->flag.cc)
return 21*nu_nubar;
else
return 41*nu_nubar;
}
if (npions == 1) {
if (e->flag.cc) {
// CC
if (e->flag.anty){
// Antineutrino
if (event1_nof(e, neutron_pdg) && event1_nof(e, pion_minus_pdg))
return -11;
if (event1_nof(e, neutron_pdg) && event1_nof(e, pion_pdg))
return -12;
if (event1_nof(e, proton_pdg) && event1_nof(e, pion_minus_pdg))
return -13;
} else {
// Neutrino
if (event1_nof(e, proton_pdg) && event1_nof(e, pion_plus_pdg))
return 11;
if (event1_nof(e, proton_pdg) && event1_nof(e, pion_pdg))
return 12;
if (event1_nof(e, neutron_pdg) && event1_nof(e, pion_plus_pdg))
return 13;
}
} else {
// Now NC
if (event1_nof(e, proton_pdg)) {
if (event1_nof(e, pion_minus_pdg))
return 33*nu_nubar;
else if (event1_nof(e, pion_pdg))
return 32*nu_nubar;
} else if (event1_nof(e, neutron_pdg)) {
if (event1_nof(e, pion_plus_pdg))
return 34*nu_nubar;
else if (event1_nof(e, pion_pdg))
return 31*nu_nubar;
}
}
}
// Eta production
if (event1_nof(e, eta_pdg)) {
if (e->flag.cc)
return 22*nu_nubar;
else {
if (event1_nof(e, neutron_pdg))
return 42*nu_nubar;
else if (event1_nof(e, proton_pdg))
return 43*nu_nubar;
}
}
// Resonant kaon production
if (event1_nof(e, lambda_pdg) == 1 && nkaons == 1) {
if (e->flag.cc && event1_nof(e, kaon_pdg))
return 23*nu_nubar;
else {
if (event1_nof(e, kaon_pdg))
return 44*nu_nubar;
else if (event1_nof(e, kaon_plus_pdg))
return 45*nu_nubar;
}
}
}
// Coherent
if (e->flag.coh) {
if (e->flag.cc && (event1_nof(e, pion_minus_pdg) +
event1_nof(e, pion_plus_pdg)))
return 16*nu_nubar;
else if (event1_nof(e, pion_pdg))
return 36*nu_nubar;
}
// DIS
if (e->flag.dis) {
if (e->flag.cc)
return 26*nu_nubar;
else
return 46*nu_nubar;
}
// If we got here, something is wrong, see what happened...
NUIS_ERR(WRN, "Unable to interpret NUWRO event, dumping info...");
Print();
return 9999;
}
void NuWroInputHandler::CalcNUISANCEKinematics() {
// std::cout << "NuWro Event Address " << fNuWroEvent << std::endl;
// Reset all variables
fNUISANCEEvent->ResetEvent();
FitEvent *evt = fNUISANCEEvent;
// Sort Event Info
evt->Mode = ConvertNuwroMode(fNuWroEvent);
if (abs(evt->Mode) > 60) {
evt->Mode = 0;
}
evt->fEventNo = 0.0;
evt->fTotCrs = 0.0;
evt->fTargetA = fNuWroEvent->par.nucleus_p + fNuWroEvent->par.nucleus_n;
evt->fTargetZ = fNuWroEvent->par.nucleus_p;
evt->fTargetPDG =
TargetUtils::GetTargetPDGFromZA(evt->fTargetZ, evt->fTargetA);
evt->fTargetH = 0;
evt->fBound = (evt->fTargetA != 1);
// Check Particle Stack
UInt_t npart_in = fNuWroEvent->in.size();
UInt_t npart_out = fNuWroEvent->out.size();
UInt_t npart_post = fNuWroEvent->post.size();
UInt_t npart = npart_in + npart_out + npart_post;
UInt_t kmax = evt->kMaxParticles;
if (npart > kmax) {
NUIS_ERR(WRN, "NUWRO has too many particles. Expanding stack.");
fNUISANCEEvent->ExpandParticleStack(npart);
}
evt->fNParticles = 0;
std::vector<particle>::iterator p_iter;
// Get the Initial State
for (p_iter = fNuWroEvent->in.begin(); p_iter != fNuWroEvent->in.end();
p_iter++) {
AddNuWroParticle(fNUISANCEEvent, (*p_iter), kInitialState, true);
}
// Try to find the FSI state particles
// Loop over the primary vertex particles
// If they match the post-FSI they haven't undergone FSI.
// If they don't match post-FSI they have undergone FSI.
for (p_iter = fNuWroEvent->out.begin(); p_iter != fNuWroEvent->out.end();
p_iter++) {
// Get the particle
particle p = (*p_iter);
// Check against all the post particles, match them
std::vector<particle>::iterator p2_iter;
bool match = false;
for (p2_iter = fNuWroEvent->post.begin();
p2_iter != fNuWroEvent->post.end(); p2_iter++) {
particle p2 = (*p2_iter);
// Check energy and pdg
// A very small cascade which changes the energy by 1E-5 MeV should be
// matched
match = (fabs(p2.E() - p.E()) < 1E-5 && p2.pdg == p.pdg);
// If we match p to p2 break the loop
if (match)
break;
}
// If we've looped through the whole particle stack of post-FSI and haven't
// found a match it's a primary particle that has been FSIed
if (!match)
AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFSIState, true);
}
// Loop over the final state particles
for (p_iter = fNuWroEvent->post.begin(); p_iter != fNuWroEvent->post.end();
p_iter++) {
particle p = (*p_iter);
// To find if it's primary or not we have to loop through the primary ones
// and match, just like above
bool match = false;
std::vector<particle>::iterator p2_iter;
for (p2_iter = fNuWroEvent->out.begin(); p2_iter != fNuWroEvent->out.end();
p2_iter++) {
particle p2 = (*p2_iter);
match = (fabs(p2.E() - p.E()) < 1E-5 && p2.pdg == p.pdg);
if (match)
break;
}
AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFinalState, match);
}
// Fill Generator Info
if (fSaveExtra)
fNuWroInfo->FillGeneratorInfo(fNuWroEvent);
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent->OrderStack();
FitParticle *ISAnyLepton = fNUISANCEEvent->GetHMISAnyLeptons();
if (ISAnyLepton) {
fNUISANCEEvent->probe_E = ISAnyLepton->E();
fNUISANCEEvent->probe_pdg = ISAnyLepton->PDG();
}
return;
}
void NuWroInputHandler::AddNuWroParticle(FitEvent *evt, particle &p, int state,
bool primary = false) {
// Add Mom
evt->fParticleMom[evt->fNParticles][0] = static_cast<vect &>(p).x;
evt->fParticleMom[evt->fNParticles][1] = static_cast<vect &>(p).y;
evt->fParticleMom[evt->fNParticles][2] = static_cast<vect &>(p).z;
evt->fParticleMom[evt->fNParticles][3] = static_cast<vect &>(p).t;
// For NuWro a particle that we've given a FSI state is a pre-FSI particle
// An initial state particle is also a primary vertex praticle
evt->fPrimaryVertex[evt->fNParticles] = primary;
// Status/PDG
evt->fParticleState[evt->fNParticles] = state;
evt->fParticlePDG[evt->fNParticles] = p.pdg;
// Add to particle count
evt->fNParticles++;
}
void NuWroInputHandler::Print(){
NUIS_LOG(SAM, "NuWro event information:" << std::endl
<< "\t\t|-> dyn = " << fNuWroEvent->dyn << std::endl
<< "\t\t|-> qel = " << fNuWroEvent->flag.qel << std::endl
<< "\t\t|-> res = " << fNuWroEvent->flag.res << std::endl
<< "\t\t|-> dis = " << fNuWroEvent->flag.dis << std::endl
<< "\t\t|-> coh = " << fNuWroEvent->flag.coh << std::endl
<< "\t\t|-> mec = " << fNuWroEvent->flag.mec << std::endl
<< "\t\t|-> hip = " << fNuWroEvent->flag.hip << std::endl
<< "\t\t|-> nc = " << fNuWroEvent->flag.nc << std::endl
<< "\t\t|-> cc = " << fNuWroEvent->flag.cc << std::endl
<< "\t\t|-> anty = " << fNuWroEvent->flag.anty << std::endl
<< "\t\t|-> res_delta = " << fNuWroEvent->flag.res_delta << std::endl
<< "\t\t|-> npi+ = " << event1_nof(fNuWroEvent, 211) << std::endl
<< "\t\t|-> npi+ = " << event1_nof(fNuWroEvent, 211) << std::endl
<< "\t\t|-> npi- = " << event1_nof(fNuWroEvent, -211) << std::endl
<< "\t\t|-> npi0 = " << event1_nof(fNuWroEvent, 111) << std::endl
<< "\t\t|-> neta = " << event1_nof(fNuWroEvent, 221) << std::endl
<< "\t\t|-> nK+ = " << event1_nof(fNuWroEvent, 321) << std::endl
<< "\t\t|-> nK0 = " << event1_nof(fNuWroEvent, 311) << std::endl
<< "\t\t|-> nlamda = " << event1_nof(fNuWroEvent, 3122) << std::endl
<< "\t\t|-> nproton = " << event1_nof(fNuWroEvent, 2212) << std::endl
<< "\t\t|-> nneutron = " << event1_nof(fNuWroEvent, 2112));
}
#endif

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 1:57 PM (1 h, 35 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486410
Default Alt Text
(15 KB)

Event Timeline