diff --git a/parameters/config.xml b/parameters/config.xml
index 3eb548b..48492b7 100644
--- a/parameters/config.xml
+++ b/parameters/config.xml
@@ -1,240 +1,232 @@
<nuisance>
<!-- # ###################################################### -->
<!-- # NUISANCE CONFIGURATION OPTIONS -->
<!-- # This file is read in by default at runtime -->
<!-- # If you want to override on a case by case bases use -q at runtime -->
<!-- # ###################################################### -->
<!-- # MAIN Configs -->
<!-- # ###################################################### -->
<!-- # Logger goes from -->
<!-- # 1 Quiet -->
<!-- # 2 Fitter -->
<!-- # 3 Samples -->
<!-- # 4 Reconfigure Loops -->
<!-- # 5 Every Event print out (SHOUT) -->
<!-- # -1 DEBUGGING -->
<config VERBOSITY='4'/>
<!-- # ERROR goes from -->
<!-- # 0 NONE -->
<!-- # 1 FATAL -->
<!-- # 2 WARN -->
<config ERROR='2'/>
<config TRACE='0'/>
<config cores='1' />
<config spline_test_throws='50' />
<config spline_cores='1' />
<config spline_chunks='20' />
<config spline_procchunk='-1' />
<config Electron_NThetaBins='4' />
<config Electron_NEnergyBins='4' />
<config Electron_ThetaWidth='1.0' />
<config Electron_EnergyWidth='0.10' />
<!-- Do we want to remove FSI, undefined and nuclear particles from the GENIE particle stack? -->
<config RemoveFSIParticles='0' />
<config RemoveUndefParticles='0' />
<config RemoveNuclearParticles='0'/>
<config MINERvASaveExtraCCQE='0' />
<!-- # Input Configs -->
<!-- # ###################################################### -->
<!-- # Default Requirements file for the externalDataFitter Package -->
<!-- # MAX Events : -1 is no cut. Events will be scaled automatically to give good xsec predictions. -->
<config MAXEVENTS='-1'/>
<!-- Include empty stacks in the THStack -->
<config includeemptystackhists='0'/>
<!-- # Turn on/off event manager -->
<!-- # EventManager enables us to only loop number of events once for multiple projections of the same measurements -->
<!-- # e.g. MiniBooNE CC1pi+ Q2 and MiniBooNE CC1pi+ Tmu would ordinarily require 2 reconfigures, but with this enabled it requires only one -->
<config EventManager='1'/>
<!-- # Event Directories -->
<!-- # Can setup default directories and use @EVENT_DIR/path to link to it -->
<config EVENT_DIR='/data2/stowell/NIWG/'/>
<config NEUT_EVENT_DIR='/data2/stowell/NIWG/neut/fit_samples_neut5.3.3/'/>
<config GENIE_EVENT_DIR='/data2/stowell/NIWG/genie/fit_samples_R.2.10.0/'/>
<config NUWRO_EVENT_DIR='/data2/stowell/NIWG/nuwro/fit_samples/'/>
<config GIBUU_EVENT_DIR='/data/GIBUU/DIR/'/>
<config SaveNuWroExtra='0' />
<!-- # In PrepareGENIE the reconstructed splines can be saved into the file -->
<config save_genie_splines='1'/>
<!-- # In InputHandler the option to regenerate NuWro flux/xsec plots is available -->
<!-- # Going to move this to its own app soon -->
<!-- # DEVEL CONFIG OPTION, don't touch! -->
<config CacheSize='0'/>
<!-- # ReWeighting Configuration Options -->
<!-- # ###################################################### -->
<!-- # Convert Dials in output statements using dial conversion card -->
<config convert_dials='0'/>
<!-- # Vetos can be used to specify RW dials NOT to be loaded in -->
<!-- # Useful if one specific one has an issue -->
<config FitWeight_fNIWGRW_veto=''/>
<config FitWeight_fNuwroRW_veto=''/>
<config FitWeight_fNeutRW_veto=''/>
<config FitWeight_fGenieRW_veto=''/>
<!-- # Output Options -->
<!-- # ###################################################### -->
<!-- # Save Nominal prediction with all rw engines at default -->
<config savenominal='0'/>
<!-- # Save prefit with values at starting values -->
<config saveprefit='0'/>
<!-- # Here's the full list of drawing options -->
<!-- DATA/MC/EVT/FINE/RATIO/MODES/SHAPE/WGHT/WEIGHTS/FLUX/XSEC/MASK/COV/INCOV/DECOMP/CANVPDG/CANVMC/SETTINGS'/>
<!-- #config drawopts DATA/MC/EVT/FINE/RATIO/MODES/SHAPE/RESIDUAL/MATRIX/FLUX/MASK/MAP -->
<!-- #config drawopts DATA/MC -->
<config drawopts='DATA/MC/EVT/FINE/RATIO/MODES/SHAPE/FLUX/XSEC/MASK/COV/INCOV/DECOMP/CANVPDG/CANVMC/SETTINGS/PROJ/CANVSLICEMC'/>
<config nuisflat_SavePreFSI='true' />
<config InterpolateSigmaQ0Histogram='1' />
<config InterpolateSigmaQ0HistogramRes='100' />
<config InterpolateSigmaQ0HistogramThrow='1' />
<config InterpolateSigmaQ0HistogramNTHROWS='100000' />
<!-- # Save the shape scaling applied with option SHAPE into the main MC hist -->
<config saveshapescaling='0'/>
<config CorrectGENIEMECNorm='1'/>
<!-- # Set style of 1D output histograms -->
<config linecolour='1'/>
<config linestyle='1'/>
<config linewidth='1'/>
<!-- # For GenericFlux -->
<config isLiteMode='0'/>
<!-- # Statistical Options -->
<!-- # ###################################################### -->
<!-- # Add MC Statistical error to likelihoods -->
<config addmcerror='0'/>
<!-- # NUISMIN Configurations -->
<!-- # ###################################################### -->
<config MAXCALLS='1000000'/>
<config MAXITERATIONS='1000000'/>
<config TOLERANCE='0.001'/>
<!-- # Number of events required in low stats routines -->
<config LOWSTATEVENTS='25000'/>
<!-- # Error band generator configs -->
<!-- # ###################################################### -->
<!-- # For -f ErrorBands creates error bands for given measurements -->
<!-- # How many throws do we want (The higher the better precision) -->
<config error_throws='500'/>
<!-- # Are we throwing uniform or according to Gaussian? -->
<!-- # Only use uniform if wanting to study the limits of a dial. -->
<config error_uniform='0'/>
<config WriteSeperateStacks='1'/>
<!-- # Other Individual Case Configs -->
<!-- # ###################################################### -->
<!-- # Covariance throw options for fake data studies with MiniBooNE data. -->
<config thrown_covariance='FULL'/>
<config throw_mc_stat='0.0'/>
<config throw_diag_syst='0'/>
<config throw_corr_syst='0'/>
<config throw_mc_stat='0'/>
<!-- # Apply a shift to the muon momentum before calculation of Q2 -->
<config muon_momentum_shift='0.0'/>
<config muon_momentum_throw='0'/>
<!-- # MINERvA Specific Configs -->
<config MINERvA_CCinc_XSec_2DEavq3_nu_hadron_cut='0'/>
<config MINERvA_CCinc_XSec_2DEavq3_nu_useq3true='0'/>
<config Modes_split_PN_NN='0'/>
<!-- Use only signal events when reconfiguring -->
<config SignalReconfigures='false'/>
<config FullEventOnSignalReconfigure="true"/>
<!-- # SciBooNE specific -->
<config SciBarDensity='1.04'/>
<config SciBarRecoDist='12.0'/>
<config PenetratingMuonEnergy='1.4'/>
<config FlatEfficiency='1.0'/>
<config NumRangeSteps='50'/>
<config UseProton='true'/>
<config UseZackEff='false'/>
<config MINERvADensity='1.04'/>
<config MINERvARecoDist='10.0'/>
<!-- Different way of reweighting in GENIE -->
<config GENIEWeightEngine_CCRESMode="kModeMaMv"/>
<!--
<config GENIEWeightEngine_CCRESMode="kModeMa"/>
<<<<<<< HEAD
-->
<config GENIEWeightEngine_CCQEMode="kModeMa"/>
-||||||| merged common ancestors
--->
-<config GENIEWeightEngine_CCQEMode="kModeZExp"/> Using z-expansion
-=======
--->
-<!-- Normal MAQE Llewelyn-Smith -->
-<config GENIEWeightEngine_CCQEMode="kModeMa"/>
->>>>>>> feature_prefsi
<!--
<config GENIEWeightEngine_CCQEMode="kModeNormAndMaShape"/> Taking shape into account (don't really use this...)
<config GENIEWeightEngine_CCQEMode="kModeZExp"/> Using z-expansion
-->
<!-- CCQE/2p2h/1pi Gaussian enhancement method -->
<!-- Apply tilt-shift weights or normal Gaussian parameters-->
<config Gaussian_Enhancement="Normal" />
<!--
<config Gaussian_Enhancement="Tilt-Shift" />
-->
<config NToyThrows='100000' />
<!-- Use NOvA Weights or not -->
<config NOvA_Weights="false" />
<!-- Tune name, for GENIE v3+ -->
<config GENIETune="G18_02a_00_000" />
<config GENIEEventGeneratorList="Default" />
<!--
<config GENIEXSecModelCCRES="genie::ReinSehgalRESPXSec" />
<config GENIEXSecModelCOH="genie::ReinSehgalCOHPiPXSec" />
<config GENIEXSecModelCCQE="genie::LwlynSmithQELCCPXSec" />
-->
<!--
<config GENIEXSecModelCCQE="genie::NievesQELCCPXSec" />
<config GENIEXSecModelCCRES="genie::BergerSehgalRESPXSec2014" />
<config GENIEXSecModelCOH="genie::BergerSehgalCOHPiPXSec2015" />
-->
<config UseShapeCovar="0" />
<config CC0piNBINS="156" />
</nuisance>
diff --git a/src/InputHandler/GENIEInputHandler.cxx b/src/InputHandler/GENIEInputHandler.cxx
index 5d5839d..b3b5420 100644
--- a/src/InputHandler/GENIEInputHandler.cxx
+++ b/src/InputHandler/GENIEInputHandler.cxx
@@ -1,586 +1,583 @@
// 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/>.
*******************************************************************************/
#ifdef __GENIE_ENABLED__
#include "GENIEInputHandler.h"
#pragma push_macro("LOG")
#undef LOG
#pragma push_macro("ERROR")
#undef ERROR
#ifdef GENIE_PRE_R3
#include "Messenger/Messenger.h"
#else
#include "Framework/Messenger/Messenger.h"
#endif
#pragma pop_macro("LOG")
#pragma pop_macro("ERROR")
#include "InputUtils.h"
GENIEGeneratorInfo::~GENIEGeneratorInfo() { DeallocateParticleStack(); }
void GENIEGeneratorInfo::AddBranchesToTree(TTree* tn) {
tn->Branch("GenieParticlePDGs", &fGenieParticlePDGs, "GenieParticlePDGs/I");
}
void GENIEGeneratorInfo::SetBranchesFromTree(TTree* tn) {
tn->SetBranchAddress("GenieParticlePDGs", &fGenieParticlePDGs);
}
void GENIEGeneratorInfo::AllocateParticleStack(int stacksize) {
fGenieParticlePDGs = new int[stacksize];
}
void GENIEGeneratorInfo::DeallocateParticleStack() {
delete fGenieParticlePDGs;
}
void GENIEGeneratorInfo::FillGeneratorInfo(NtpMCEventRecord* ntpl) {
Reset();
// Check for GENIE Event
if (!ntpl) return;
if (!ntpl->event) return;
// Cast Event Record
GHepRecord* ghep = static_cast<GHepRecord*>(ntpl->event);
if (!ghep) return;
// Fill Particle Stack
GHepParticle* p = 0;
TObjArrayIter iter(ghep);
// Loop over all particles
int i = 0;
while ((p = (dynamic_cast<genie::GHepParticle*>((iter).Next())))) {
if (!p) continue;
// Get PDG
fGenieParticlePDGs[i] = p->Pdg();
i++;
}
}
void GENIEGeneratorInfo::Reset() {
for (int i = 0; i < kMaxParticles; i++) {
fGenieParticlePDGs[i] = 0;
}
}
GENIEInputHandler::GENIEInputHandler(std::string const& handle,
std::string const& rawinputs) {
LOG(SAM) << "Creating GENIEInputHandler : " << handle << std::endl;
genie::Messenger::Instance()->SetPriorityLevel("GHepUtils", pFATAL);
// Run a joint input handling
fName = handle;
// Setup the TChain
fGENIETree = new TChain("gtree");
fSaveExtra = FitPar::Config().GetParB("SaveExtraGenie");
fCacheSize = FitPar::Config().GetParI("CacheSize");
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
// Are we running with NOvA weights
fNOvAWeights = FitPar::Config().GetParB("NOvA_Weights");
MAQEw = 1.0;
NonResw = 1.0;
RPAQEw = 1.0;
RPARESw = 1.0;
MECw = 1.0;
DISw = 1.0;
NOVAw = 1.0;
// 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(
InputUtils::ExpandInputDirectories(inputs[inp_it]).c_str(), "READ");
if (!inp_file or inp_file->IsZombie()) {
THROW("GENIE File IsZombie() at : '"
<< inputs[inp_it] << "'" << std::endl
<< "Check that your file paths are correct and the file exists!"
<< std::endl
<< "$ ls -lh " << inputs[inp_it]);
}
// Get Flux/Event hist
TH1D* fluxhist = (TH1D*)inp_file->Get("nuisance_flux");
TH1D* eventhist = (TH1D*)inp_file->Get("nuisance_events");
if (!fluxhist or !eventhist) {
ERROR(FTL, "Input File Contents: " << inputs[inp_it]);
inp_file->ls();
THROW("GENIE FILE doesn't contain flux/xsec info."
<< std::endl
<< "Try running the app PrepareGENIE first on :" << inputs[inp_it]
<< std::endl
<< "$ PrepareGENIE -h");
}
// Get N Events
TTree* genietree = (TTree*)inp_file->Get("gtree");
if (!genietree) {
ERROR(FTL, "gtree not located in GENIE file: " << inputs[inp_it]);
THROW("Check your inputs, they may need to be completely regenerated!");
throw;
}
int nevents = genietree->GetEntries();
if (nevents <= 0) {
THROW("Trying to a TTree with "
<< nevents << " to TChain from : " << inputs[inp_it]);
}
// Check for precomputed weights
TTree *weighttree = (TTree*)inp_file->Get("nova_wgts");
if (fNOvAWeights) {
if (!weighttree) {
THROW("Did not find nova_wgts tree in file " << inputs[inp_it] << " but you specified it" << std::endl);
} else {
LOG(FIT) << "Found nova_wgts tree in file " << inputs[inp_it] << std::endl;
}
}
// Register input to form flux/event rate hists
RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist);
// Add To TChain
fGENIETree->AddFile(inputs[inp_it].c_str());
if (weighttree != NULL) fGENIETree->AddFriend(weighttree);
}
// Registor all our file inputs
SetupJointInputs();
// Assign to tree
fEventType = kGENIE;
fGenieNtpl = NULL;
fGENIETree->SetBranchAddress("gmcrec", &fGenieNtpl);
// Set up the custom weights
if (fNOvAWeights) {
fGENIETree->SetBranchAddress("MAQEwgt", &MAQEw);
fGENIETree->SetBranchAddress("nonResNormWgt", &NonResw);
fGENIETree->SetBranchAddress("RPAQEWgt", &RPAQEw);
fGENIETree->SetBranchAddress("RPARESWgt", &RPARESw);
fGENIETree->SetBranchAddress("MECWgt", &MECw);
fGENIETree->SetBranchAddress("DISWgt", &DISw);
fGENIETree->SetBranchAddress("nova2018CVWgt", &NOVAw);
}
// Libraries should be seen but not heard...
StopTalking();
fGENIETree->GetEntry(0);
StartTalking();
// Create Fit Event
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->SetGenieEvent(fGenieNtpl);
if (fSaveExtra) {
fGenieInfo = new GENIEGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fGenieInfo);
}
fNUISANCEEvent->HardReset();
};
GENIEInputHandler::~GENIEInputHandler() {
// if (fGenieGHep) delete fGenieGHep;
// if (fGenieNtpl) delete fGenieNtpl;
// if (fGENIETree) delete fGENIETree;
// if (fGenieInfo) delete fGenieInfo;
}
void GENIEInputHandler::CreateCache() {
if (fCacheSize > 0) {
// fGENIETree->SetCacheEntryRange(0, fNEvents);
fGENIETree->AddBranchToCache("*", 1);
fGENIETree->SetCacheSize(fCacheSize);
}
}
void GENIEInputHandler::RemoveCache() {
// fGENIETree->SetCacheEntryRange(0, fNEvents);
fGENIETree->AddBranchToCache("*", 0);
fGENIETree->SetCacheSize(0);
}
FitEvent* GENIEInputHandler::GetNuisanceEvent(const UInt_t entry, const bool lightweight) {
if (entry >= (UInt_t)fNEvents) return NULL;
// Clear the previous event (See Note 1 in ROOT TClonesArray documentation)
if (fGenieNtpl) {
fGenieNtpl->Clear();
}
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fGENIETree->GetEntry(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
#ifdef __PROB3PP_ENABLED__
else {
// Check for GENIE Event
if (!fGenieNtpl) return NULL;
if (!fGenieNtpl->event) return NULL;
// Cast Event Record
fGenieGHep = static_cast<GHepRecord*>(fGenieNtpl->event);
if (!fGenieGHep) return NULL;
TObjArrayIter iter(fGenieGHep);
genie::GHepParticle* p;
while ((p = (dynamic_cast<genie::GHepParticle*>((iter).Next())))) {
if (!p) {
continue;
}
// Get Status
int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode);
if (state != genie::kIStInitialState) {
continue;
}
fNUISANCEEvent->probe_E = p->E() * 1.E3;
fNUISANCEEvent->probe_pdg = p->Pdg();
break;
}
}
#endif
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
return fNUISANCEEvent;
}
int GENIEInputHandler::GetGENIEParticleStatus(genie::GHepParticle* p, int mode) {
/*
kIStUndefined = -1,
kIStInitialState = 0, / generator-level initial state /
kIStStableFinalState = 1, / generator-level final state:
particles to be tracked by detector-level MC /
kIStIntermediateState = 2,
kIStDecayedState = 3,
kIStCorrelatedNucleon = 10,
kIStNucleonTarget = 11,
kIStDISPreFragmHadronicState = 12,
kIStPreDecayResonantState = 13,
kIStHadronInTheNucleus = 14, / hadrons inside the nucleus: marked
for hadron transport modules to act on /
kIStFinalStateNuclearRemnant = 15, / low energy nuclear fragments
entering the record collectively as a 'hadronic blob' pseudo-particle /
kIStNucleonClusterTarget = 16, // for composite nucleons before
phase space decay
*/
int state = kUndefinedState;
switch (p->Status()) {
case genie::kIStNucleonTarget:
case genie::kIStInitialState:
case genie::kIStCorrelatedNucleon:
case genie::kIStNucleonClusterTarget:
state = kInitialState;
break;
case genie::kIStStableFinalState:
state = kFinalState;
break;
case genie::kIStHadronInTheNucleus:
if (abs(mode) == 2)
state = kInitialState;
else
state = kFSIState;
break;
case genie::kIStPreDecayResonantState:
case genie::kIStDISPreFragmHadronicState:
case genie::kIStIntermediateState:
state = kFSIState;
break;
case genie::kIStFinalStateNuclearRemnant:
case genie::kIStUndefined:
case genie::kIStDecayedState:
default:
break;
}
// Flag to remove nuclear part in genie
if (p->Pdg() > 1000000) {
if (state == kInitialState)
state = kNuclearInitial;
else if (state == kFinalState)
state = kNuclearRemnant;
}
return state;
}
#endif
#ifdef __GENIE_ENABLED__
int GENIEInputHandler::ConvertGENIEReactionCode(GHepRecord* gheprec) {
// Electron Scattering
if (gheprec->Summary()->ProcInfo().IsEM()) {
if (gheprec->Summary()->InitState().ProbePdg() == 11) {
if (gheprec->Summary()->ProcInfo().IsQuasiElastic())
return 1;
else if (gheprec->Summary()->ProcInfo().IsMEC())
return 2;
else if (gheprec->Summary()->ProcInfo().IsResonant())
return 13;
else if (gheprec->Summary()->ProcInfo().IsDeepInelastic())
return 26;
else {
ERROR(WRN,
"Unknown GENIE Electron Scattering Mode!"
<< std::endl
<< "ScatteringTypeId = "
<< gheprec->Summary()->ProcInfo().ScatteringTypeId() << " "
<< "InteractionTypeId = "
<< gheprec->Summary()->ProcInfo().InteractionTypeId()
<< std::endl
<< genie::ScatteringType::AsString(
gheprec->Summary()->ProcInfo().ScatteringTypeId())
<< " "
<< genie::InteractionType::AsString(
gheprec->Summary()->ProcInfo().InteractionTypeId())
<< " " << gheprec->Summary()->ProcInfo().IsMEC());
return 0;
}
}
// Weak CC
} else if (gheprec->Summary()->ProcInfo().IsWeakCC()) {
// CC MEC
if (gheprec->Summary()->ProcInfo().IsMEC()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 2;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -2;
// CC OTHER
} else {
return utils::ghep::NeutReactionCode(gheprec);
}
// Weak NC
} else if (gheprec->Summary()->ProcInfo().IsWeakNC()) {
// NC MEC
if (gheprec->Summary()->ProcInfo().IsMEC()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 32;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -32;
// NC OTHER
} else {
return utils::ghep::NeutReactionCode(gheprec);
}
}
return 0;
}
void GENIEInputHandler::CalcNUISANCEKinematics() {
// Reset all variables
fNUISANCEEvent->ResetEvent();
// Check for GENIE Event
if (!fGenieNtpl) return;
if (!fGenieNtpl->event) return;
// Cast Event Record
fGenieGHep = static_cast<GHepRecord*>(fGenieNtpl->event);
if (!fGenieGHep) return;
// Convert GENIE Reaction Code
fNUISANCEEvent->Mode = ConvertGENIEReactionCode(fGenieGHep);
// Set Event Info
fNUISANCEEvent->fEventNo = 0.0;
fNUISANCEEvent->fTotCrs = fGenieGHep->XSec();
// Have a bool storing if interaction happened on free or bound nucleon
bool IsFree = false;
// Set the TargetPDG
if (fGenieGHep->TargetNucleus() != NULL) {
fNUISANCEEvent->fTargetPDG = fGenieGHep->TargetNucleus()->Pdg();
IsFree = false;
// Sometimes GENIE scatters off free nucleons, electrons, photons
// In which TargetNucleus is NULL and we need to find the initial state particle
} else {
// Check the particle is an initial state particle
// Follows GHepRecord::TargetNucleusPosition but doesn't do check on pdg::IsIon
GHepParticle *p = fGenieGHep->Particle(1);
// Check that particle 1 actually exists
if (!p) {
ERR(FTL) << "Can't find particle 1 for GHepRecord" << std::endl;
throw;
}
// If not an ion but is an initial state particle
if (!pdg::IsIon(p->Pdg()) &&
p->Status() == kIStInitialState) {
IsFree = true;
fNUISANCEEvent->fTargetPDG = p->Pdg();
// Catch if something strange happens:
// Here particle 1 is not an initial state particle OR
// particle 1 is an ion OR
// both
} else {
if (pdg::IsIon(p->Pdg())) {
ERR(FTL) << "Particle 1 in GHepRecord stack is an ion but isn't an initial state particle" << std::endl;
throw;
} else {
ERR(FTL) << "Particle 1 in GHepRecord stack is not an ion but is an initial state particle" << std::endl;
throw;
}
}
}
// Set the A and Z and H from the target PDG
// Depends on if we scattered off a free or bound nucleon
if (!IsFree) {
fNUISANCEEvent->fTargetA = TargetUtils::GetTargetAFromPDG(fNUISANCEEvent->fTargetPDG);
fNUISANCEEvent->fTargetZ = TargetUtils::GetTargetZFromPDG(fNUISANCEEvent->fTargetPDG);
fNUISANCEEvent->fTargetH = 0;
} else {
// If free proton scattering
if (fNUISANCEEvent->fTargetPDG == 2212) {
fNUISANCEEvent->fTargetA = 1;
fNUISANCEEvent->fTargetZ = 1;
fNUISANCEEvent->fTargetH = 1;
// If free neutron scattering
} else if (fNUISANCEEvent->fTargetPDG == 2112) {
fNUISANCEEvent->fTargetA = 0;
fNUISANCEEvent->fTargetZ = 1;
fNUISANCEEvent->fTargetH = 0;
// If neither
} else {
fNUISANCEEvent->fTargetA = 0;
fNUISANCEEvent->fTargetZ = 0;
fNUISANCEEvent->fTargetH = 0;
}
}
fNUISANCEEvent->fBound = !IsFree;
fNUISANCEEvent->InputWeight = 1.0; //(1E+38 / genie::units::cm2) * fGenieGHep->XSec();
// And the custom weights
if (fNOvAWeights) {
fNUISANCEEvent->CustomWeight = NOVAw;
fNUISANCEEvent->CustomWeightArray[0] = MAQEw;
fNUISANCEEvent->CustomWeightArray[1] = NonResw;
fNUISANCEEvent->CustomWeightArray[2] = RPAQEw;
fNUISANCEEvent->CustomWeightArray[3] = RPARESw;
fNUISANCEEvent->CustomWeightArray[4] = MECw;
fNUISANCEEvent->CustomWeightArray[5] = NOVAw;
} else {
fNUISANCEEvent->CustomWeight = 1.0;
fNUISANCEEvent->CustomWeightArray[0] = 1.0;
fNUISANCEEvent->CustomWeightArray[1] = 1.0;
fNUISANCEEvent->CustomWeightArray[2] = 1.0;
fNUISANCEEvent->CustomWeightArray[3] = 1.0;
fNUISANCEEvent->CustomWeightArray[4] = 1.0;
fNUISANCEEvent->CustomWeightArray[5] = 1.0;
}
// Get N Particle Stack
unsigned int npart = fGenieGHep->GetEntries();
unsigned int kmax = fNUISANCEEvent->kMaxParticles;
if (npart > kmax) {
ERR(WRN) << "GENIE has too many particles, expanding stack." << std::endl;
fNUISANCEEvent->ExpandParticleStack(npart);
}
// Fill Particle Stack
GHepParticle* p = 0;
TObjArrayIter iter(fGenieGHep);
fNUISANCEEvent->fNParticles = 0;
// Loop over all particles
while ((p = (dynamic_cast<genie::GHepParticle*>((iter).Next())))) {
if (!p) continue;
// Get Status
int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode);
// Remove Undefined
if (kRemoveUndefParticles && state == kUndefinedState) continue;
// Remove FSI
if (kRemoveFSIParticles && state == kFSIState) continue;
if (kRemoveNuclearParticles &&
(state == kNuclearInitial || state == kNuclearRemnant))
continue;
// Fill Vectors
int curpart = fNUISANCEEvent->fNParticles;
fNUISANCEEvent->fParticleState[curpart] = state;
// Mom
fNUISANCEEvent->fParticleMom[curpart][0] = p->Px() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][1] = p->Py() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][2] = p->Pz() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][3] = p->E() * 1.E3;
// PDG
fNUISANCEEvent->fParticlePDG[curpart] = p->Pdg();
// Set if the particle was on the fundamental vertex
fNUISANCEEvent->fPrimaryVertex[curpart] = (p->FirstMother() < 2);
- //std::cout << curpart << " " << p->Pdg() << " " << p->FirstMother() << std::endl;
- //std::cout << p->LastMother() << std::endl;
- //std::cout << "**" << std::endl;
// Add to N particle count
fNUISANCEEvent->fNParticles++;
// Extra Check incase GENIE fails.
if ((UInt_t)fNUISANCEEvent->fNParticles == kmax) {
ERR(WRN) << "Number of GENIE Particles exceeds maximum!" << std::endl;
ERR(WRN) << "Extend kMax, or run without including FSI particles!"
<< std::endl;
break;
}
}
// Fill Extra Stack
if (fSaveExtra) fGenieInfo->FillGeneratorInfo(fGenieNtpl);
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent->OrderStack();
FitParticle* ISNeutralLepton =
fNUISANCEEvent->GetHMISParticle(PhysConst::pdg_neutrinos);
if (ISNeutralLepton) {
fNUISANCEEvent->probe_E = ISNeutralLepton->E();
fNUISANCEEvent->probe_pdg = ISNeutralLepton->PDG();
}
return;
}
void GENIEInputHandler::Print() {}
#endif
diff --git a/src/MCStudies/GenericFlux_Vectors.cxx b/src/MCStudies/GenericFlux_Vectors.cxx
index a8cec51..59b267a 100644
--- a/src/MCStudies/GenericFlux_Vectors.cxx
+++ b/src/MCStudies/GenericFlux_Vectors.cxx
@@ -1,383 +1,382 @@
// 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 "GenericFlux_Vectors.h"
GenericFlux_Vectors::GenericFlux_Vectors(std::string name,
std::string inputfile, FitWeight *rw,
std::string type,
std::string fakeDataFile) {
// Measurement Details
fName = name;
eventVariables = NULL;
// Define our energy range for flux calcs
EnuMin = 0.;
EnuMax = 1E10; // Arbritrarily high energy limit
if (Config::HasPar("EnuMin")) {
EnuMin = Config::GetParD("EnuMin");
}
if (Config::HasPar("EnuMax")) {
EnuMax = Config::GetParD("EnuMax");
}
SavePreFSI = Config::Get().GetParB("nuisflat_SavePreFSI");
LOG(SAM) << "Running GenericFlux_Vectors saving pre-FSI particles? " << SavePreFSI << std::endl;
// 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);
// 1. The generator is organised in SetupMeasurement so it gives the
// cross-section in "per nucleon" units.
// So some extra scaling for a specific measurement may be required. For
// Example to get a "per neutron" measurement on carbon
// which we do here, we have to multiple by the number of nucleons 12 and
// divide by the number of neutrons 6.
// N.B. MeasurementBase::PredictedEventRate includes the 1E-38 factor that is
// often included here in other classes that directly integrate the event
// histogram. This method is used here as it now respects EnuMin and EnuMax
// correctly.
this->fScaleFactor =
(this->PredictedEventRate("width", 0, EnuMax) / double(fNEvents)) /
this->TotalIntegratedFlux("width");
LOG(SAM) << " Generic Flux Scaling Factor = " << fScaleFactor
<< " [= " << (GetEventHistogram()->Integral("width") * 1E-38) << "/("
<< (fNEvents + 0.) << "*" << TotalIntegratedFlux("width") << ")]"
<< std::endl;
if (fScaleFactor <= 0.0) {
ERR(WRN) << "SCALE FACTOR TOO LOW " << std::endl;
throw;
}
// Setup our TTrees
this->AddEventVariablesToTree();
this->AddSignalFlagsToTree();
}
void GenericFlux_Vectors::AddEventVariablesToTree() {
// Setup the TTree to save everything
if (!eventVariables) {
Config::Get().out->cd();
eventVariables = new TTree((this->fName + "_VARS").c_str(),
(this->fName + "_VARS").c_str());
}
LOG(SAM) << "Adding Event Variables" << std::endl;
eventVariables->Branch("Mode", &Mode, "Mode/I");
eventVariables->Branch("cc", &cc, "cc/B");
eventVariables->Branch("PDGnu", &PDGnu, "PDGnu/I");
eventVariables->Branch("Enu_true", &Enu_true, "Enu_true/F");
eventVariables->Branch("tgt", &tgt, "tgt/I");
eventVariables->Branch("tgta", &tgta, "tgta/I");
eventVariables->Branch("tgtz", &tgtz, "tgtz/I");
eventVariables->Branch("PDGLep", &PDGLep, "PDGLep/I");
eventVariables->Branch("ELep", &ELep, "ELep/F");
eventVariables->Branch("CosLep", &CosLep, "CosLep/F");
// Basic interaction kinematics
eventVariables->Branch("Q2", &Q2, "Q2/F");
eventVariables->Branch("q0", &q0, "q0/F");
eventVariables->Branch("q3", &q3, "q3/F");
eventVariables->Branch("Enu_QE", &Enu_QE, "Enu_QE/F");
eventVariables->Branch("Q2_QE", &Q2_QE, "Q2_QE/F");
eventVariables->Branch("W_nuc_rest", &W_nuc_rest, "W_nuc_rest/F");
eventVariables->Branch("W", &W, "W/F");
eventVariables->Branch("W_genie", &W_genie, "W_genie/F");
eventVariables->Branch("x", &x, "x/F");
eventVariables->Branch("y", &y, "y/F");
eventVariables->Branch("Eav", &Eav, "Eav/F");
eventVariables->Branch("EavAlt", &EavAlt, "EavAlt/F");
eventVariables->Branch("dalphat", &dalphat, "dalphat/F");
eventVariables->Branch("dpt", &dpt, "dpt/F");
eventVariables->Branch("dphit", &dphit, "dphit/F");
eventVariables->Branch("pnreco_C", &pnreco_C, "pnreco_C/F");
// Save outgoing particle vectors
eventVariables->Branch("nfsp", &nfsp, "nfsp/I");
eventVariables->Branch("px", px, "px[nfsp]/F");
eventVariables->Branch("py", py, "py[nfsp]/F");
eventVariables->Branch("pz", pz, "pz[nfsp]/F");
eventVariables->Branch("E", E, "E[nfsp]/F");
eventVariables->Branch("pdg", pdg, "pdg[nfsp]/I");
eventVariables->Branch("status", status, "status[nfsp]/I");
eventVariables->Branch("isalive", isalive, "isalive[nfsp]/O");
eventVariables->Branch("isprimary", isprimary, "isprimary[nfsp]/O");
// Event Scaling Information
eventVariables->Branch("Weight", &Weight, "Weight/F");
eventVariables->Branch("InputWeight", &InputWeight, "InputWeight/F");
eventVariables->Branch("RWWeight", &RWWeight, "RWWeight/F");
// Should be a double because may be 1E-39 and less
eventVariables->Branch("fScaleFactor", &fScaleFactor, "fScaleFactor/D");
// The customs
eventVariables->Branch("CustomWeight", &CustomWeight, "CustomWeight/F");
eventVariables->Branch("CustomWeightArray", CustomWeightArray, "CustomWeightArray[6]/F");
return;
}
void GenericFlux_Vectors::FillEventVariables(FitEvent *event) {
ResetVariables();
// Fill Signal Variables
FillSignalFlags(event);
LOG(DEB) << "Filling signal" << std::endl;
// Now fill the information
Mode = event->Mode;
cc = (abs(event->Mode) < 30);
// Get the incoming neutrino and outgoing lepton
FitParticle *nu = event->GetNeutrinoIn();
FitParticle *lep = event->GetHMFSAnyLepton();
PDGnu = nu->fPID;
Enu_true = nu->fP.E() / 1E3;
tgt = event->fTargetPDG;
tgta = event->fTargetA;
tgtz = event->fTargetZ;
if (lep != NULL) {
PDGLep = lep->fPID;
ELep = lep->fP.E() / 1E3;
CosLep = cos(nu->fP.Vect().Angle(lep->fP.Vect()));
// Basic interaction kinematics
Q2 = -1 * (nu->fP - lep->fP).Mag2() / 1E6;
q0 = (nu->fP - lep->fP).E() / 1E3;
q3 = (nu->fP - lep->fP).Vect().Mag() / 1E3;
// These assume C12 binding from MINERvA... not ideal
Enu_QE = FitUtils::EnuQErec(lep->fP, CosLep, 34., true);
Q2_QE = FitUtils::Q2QErec(lep->fP, CosLep, 34., true);
Eav = FitUtils::GetErecoil_MINERvA_LowRecoil(event)/1.E3;
EavAlt = FitUtils::Eavailable(event)/1.E3;
// Get W_true with assumption of initial state nucleon at rest
float m_n = (float)PhysConst::mass_proton;
// Q2 assuming nucleon at rest
W_nuc_rest = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n);
// True Q2
W = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n);
x = Q2 / (2 * m_n * q0);
y = 1 - ELep / Enu_true;
dalphat = FitUtils::Get_STV_dalphat(event, PDGnu, true);
dpt = FitUtils::Get_STV_dpt(event, PDGnu, true);
dphit = FitUtils::Get_STV_dphit(event, PDGnu, true);
pnreco_C = FitUtils::Get_pn_reco_C(event, PDGnu, true);
}
// Loop over the particles and store all the final state particles in a vector
for (UInt_t i = 0; i < event->Npart(); ++i) {
bool part_alive = event->PartInfo(i)->fIsAlive &&
event->PartInfo(i)->Status() == kFinalState;
if (!SavePreFSI) {
if (!part_alive) continue;
}
- //std::cout << event->PartInfo(i)->Status() << std::endl;
partList.push_back(event->PartInfo(i));
}
// Save outgoing particle vectors
nfsp = (int)partList.size();
for (int i = 0; i < nfsp; ++i) {
px[i] = partList[i]->fP.X() / 1E3;
py[i] = partList[i]->fP.Y() / 1E3;
pz[i] = partList[i]->fP.Z() / 1E3;
E[i] = partList[i]->fP.E() / 1E3;
pdg[i] = partList[i]->fPID;
status[i] = partList[i]->Status();
isalive[i] = partList[i]->fIsAlive;
isprimary[i] = event->fPrimaryVertex[i];
}
#ifdef __GENIE_ENABLED__
if (event->fType == kGENIE) {
EventRecord * gevent = static_cast<EventRecord*>(event->genie_event->event);
const Interaction * interaction = gevent->Summary();
const Kinematics & kine = interaction->Kine();
double W_genie = kine.W();
}
#endif
// Fill event weights
Weight = event->RWWeight * event->InputWeight;
RWWeight = event->RWWeight;
InputWeight = event->InputWeight;
// And the Customs
CustomWeight = event->CustomWeight;
for (int i = 0; i < 6; ++i) {
CustomWeightArray[i] = event->CustomWeightArray[i];
}
// Fill the eventVariables Tree
eventVariables->Fill();
return;
};
//********************************************************************
void GenericFlux_Vectors::ResetVariables() {
//********************************************************************
cc = false;
// Reset all Function used to extract any variables of interest to the event
Mode = PDGnu = tgt = tgta = tgtz = PDGLep = 0;
Enu_true = ELep = CosLep = Q2 = q0 = q3 = Enu_QE = Q2_QE = W_nuc_rest = W = x = y = Eav = EavAlt = -999.9;
W_genie = -999;
// Other fun variables
// MINERvA-like ones
dalphat = dpt = dphit = pnreco_C = -999.99;
nfsp = 0;
for (int i = 0; i < kMAX; ++i){
px[i] = py[i] = pz[i] = E[i] = -999;
pdg[i] = 0;
status[i] = -999;
isalive[i] = false;
isprimary[i] = false;
}
Weight = InputWeight = RWWeight = 0.0;
CustomWeight = 0.0;
for (int i = 0; i < 6; ++i) CustomWeightArray[i] = 0.0;
partList.clear();
flagCCINC = flagNCINC = flagCCQE = flagCC0pi = flagCCQELike = flagNCEL = flagNC0pi = flagCCcoh = flagNCcoh = flagCC1pip = flagNC1pip = flagCC1pim = flagNC1pim = flagCC1pi0 = flagNC1pi0 = false;
}
//********************************************************************
void GenericFlux_Vectors::FillSignalFlags(FitEvent *event) {
//********************************************************************
// Some example flags are given from SignalDef.
// See src/Utils/SignalDef.cxx for more.
int nuPDG = event->PartInfo(0)->fPID;
// Generic signal flags
flagCCINC = SignalDef::isCCINC(event, nuPDG);
flagNCINC = SignalDef::isNCINC(event, nuPDG);
flagCCQE = SignalDef::isCCQE(event, nuPDG);
flagCCQELike = SignalDef::isCCQELike(event, nuPDG);
flagCC0pi = SignalDef::isCC0pi(event, nuPDG);
flagNCEL = SignalDef::isNCEL(event, nuPDG);
flagNC0pi = SignalDef::isNC0pi(event, nuPDG);
flagCCcoh = SignalDef::isCCCOH(event, nuPDG, 211);
flagNCcoh = SignalDef::isNCCOH(event, nuPDG, 111);
flagCC1pip = SignalDef::isCC1pi(event, nuPDG, 211);
flagNC1pip = SignalDef::isNC1pi(event, nuPDG, 211);
flagCC1pim = SignalDef::isCC1pi(event, nuPDG, -211);
flagNC1pim = SignalDef::isNC1pi(event, nuPDG, -211);
flagCC1pi0 = SignalDef::isCC1pi(event, nuPDG, 111);
flagNC1pi0 = SignalDef::isNC1pi(event, nuPDG, 111);
}
void GenericFlux_Vectors::AddSignalFlagsToTree() {
if (!eventVariables) {
Config::Get().out->cd();
eventVariables = new TTree((this->fName + "_VARS").c_str(),
(this->fName + "_VARS").c_str());
}
LOG(SAM) << "Adding signal flags" << std::endl;
// Signal Definitions from SignalDef.cxx
eventVariables->Branch("flagCCINC", &flagCCINC, "flagCCINC/O");
eventVariables->Branch("flagNCINC", &flagNCINC, "flagNCINC/O");
eventVariables->Branch("flagCCQE", &flagCCQE, "flagCCQE/O");
eventVariables->Branch("flagCC0pi", &flagCC0pi, "flagCC0pi/O");
eventVariables->Branch("flagCCQELike", &flagCCQELike, "flagCCQELike/O");
eventVariables->Branch("flagNCEL", &flagNCEL, "flagNCEL/O");
eventVariables->Branch("flagNC0pi", &flagNC0pi, "flagNC0pi/O");
eventVariables->Branch("flagCCcoh", &flagCCcoh, "flagCCcoh/O");
eventVariables->Branch("flagNCcoh", &flagNCcoh, "flagNCcoh/O");
eventVariables->Branch("flagCC1pip", &flagCC1pip, "flagCC1pip/O");
eventVariables->Branch("flagNC1pip", &flagNC1pip, "flagNC1pip/O");
eventVariables->Branch("flagCC1pim", &flagCC1pim, "flagCC1pim/O");
eventVariables->Branch("flagNC1pim", &flagNC1pim, "flagNC1pim/O");
eventVariables->Branch("flagCC1pi0", &flagCC1pi0, "flagCC1pi0/O");
eventVariables->Branch("flagNC1pi0", &flagNC1pi0, "flagNC1pi0/O");
};
void GenericFlux_Vectors::Write(std::string drawOpt) {
// First save the TTree
eventVariables->Write();
// Save Flux and Event Histograms too
GetInput()->GetFluxHistogram()->Write();
GetInput()->GetEventHistogram()->Write();
return;
}
// Override functions which aren't really necessary
bool GenericFlux_Vectors::isSignal(FitEvent *event) {
(void)event;
return true;
};
void GenericFlux_Vectors::ScaleEvents() { return; }
void GenericFlux_Vectors::ApplyNormScale(float norm) {
this->fCurrentNorm = norm;
return;
}
void GenericFlux_Vectors::FillHistograms() { return; }
void GenericFlux_Vectors::ResetAll() {
eventVariables->Reset();
return;
}
float GenericFlux_Vectors::GetChi2() { return 0.0; }