diff --git a/src/InputHandler/InputHandler.cxx b/src/InputHandler/InputHandler.cxx
index 304b445..9101762 100644
--- a/src/InputHandler/InputHandler.cxx
+++ b/src/InputHandler/InputHandler.cxx
@@ -1,311 +1,313 @@
// 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 .
*******************************************************************************/
#include "InputHandler.h"
#include "InputUtils.h"
InputHandlerBase::InputHandlerBase() {
fName = "";
fFluxHist = NULL;
fEventHist = NULL;
fNEvents = 0;
fNUISANCEEvent = NULL;
fBaseEvent = NULL;
kRemoveUndefParticles = FitPar::Config().GetParB("RemoveUndefParticles");
kRemoveFSIParticles = FitPar::Config().GetParB("RemoveFSIParticles");
kRemoveNuclearParticles = FitPar::Config().GetParB("RemoveNuclearParticles");
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
fTTreePerformance = NULL;
fSkip = 0;
if (FitPar::Config().HasConfig("NSKIPEVENTS")) {
fSkip = FitPar::Config().GetParI("NSKIPEVENTS");
+ std::cout << "Skipping " << fSkip << " events when reading input trees."
+ << std::endl;
}
};
InputHandlerBase::~InputHandlerBase() {
if (fFluxHist)
delete fFluxHist;
if (fEventHist)
delete fEventHist;
// if (fXSecHist) delete fXSecHist;
// if (fNUISANCEEvent) delete fNUISANCEEvent;
jointfluxinputs.clear();
jointeventinputs.clear();
jointindexlow.clear();
jointindexhigh.clear();
jointindexallowed.clear();
jointindexscale.clear();
// if (fTTreePerformance) {
// fTTreePerformance->SaveAs(("ttreeperfstats_" + fName +
// ".root").c_str());
// }
}
void InputHandlerBase::Print(){};
TH1D *InputHandlerBase::GetXSecHistogram(void) {
fXSecHist = (TH1D *)fFluxHist->Clone();
fXSecHist->Divide(fEventHist);
return fXSecHist;
};
double InputHandlerBase::PredictedEventRate(double low, double high,
std::string intOpt) {
Int_t minBin = fEventHist->GetXaxis()->FindFixBin(low);
Int_t maxBin = fEventHist->GetXaxis()->FindFixBin(high);
if ((fEventHist->IsBinOverflow(minBin) && (low != -9999.9))) {
minBin = 1;
}
if ((fEventHist->IsBinOverflow(maxBin) && (high != -9999.9))) {
maxBin = fEventHist->GetXaxis()->GetNbins() + 1;
}
// If we are within a single bin
if (minBin == maxBin) {
// Get the contained fraction of the single bin's width
return ((high - low) / fEventHist->GetXaxis()->GetBinWidth(minBin)) *
fEventHist->Integral(minBin, minBin, intOpt.c_str());
}
double lowBinUpEdge = fEventHist->GetXaxis()->GetBinUpEdge(minBin);
double highBinLowEdge = fEventHist->GetXaxis()->GetBinLowEdge(maxBin);
double lowBinfracIntegral =
((lowBinUpEdge - low) / fEventHist->GetXaxis()->GetBinWidth(minBin)) *
fEventHist->Integral(minBin, minBin, intOpt.c_str());
double highBinfracIntegral =
((high - highBinLowEdge) / fEventHist->GetXaxis()->GetBinWidth(maxBin)) *
fEventHist->Integral(maxBin, maxBin, intOpt.c_str());
// If they are neighbouring bins
if ((minBin + 1) == maxBin) {
std::cout << "Get lowfrac + highfrac" << std::endl;
// Get the contained fraction of the two bin's width
return lowBinfracIntegral + highBinfracIntegral;
}
double ContainedIntegral =
fEventHist->Integral(minBin + 1, maxBin - 1, intOpt.c_str());
// If there are filled bins between them
return lowBinfracIntegral + highBinfracIntegral + ContainedIntegral;
};
double InputHandlerBase::TotalIntegratedFlux(double low, double high,
std::string intOpt) {
Int_t minBin = fFluxHist->GetXaxis()->FindFixBin(low);
Int_t maxBin = fFluxHist->GetXaxis()->FindFixBin(high);
if ((fFluxHist->IsBinOverflow(minBin) && (low != -9999.9))) {
minBin = 1;
}
if ((fFluxHist->IsBinOverflow(maxBin) && (high != -9999.9))) {
maxBin = fFluxHist->GetXaxis()->GetNbins();
high = fFluxHist->GetXaxis()->GetBinLowEdge(maxBin + 1);
}
// If we are within a single bin
if (minBin == maxBin) {
// Get the contained fraction of the single bin's width
return ((high - low) / fFluxHist->GetXaxis()->GetBinWidth(minBin)) *
fFluxHist->Integral(minBin, minBin, intOpt.c_str());
}
double lowBinUpEdge = fFluxHist->GetXaxis()->GetBinUpEdge(minBin);
double highBinLowEdge = fFluxHist->GetXaxis()->GetBinLowEdge(maxBin);
double lowBinfracIntegral =
((lowBinUpEdge - low) / fFluxHist->GetXaxis()->GetBinWidth(minBin)) *
fFluxHist->Integral(minBin, minBin, intOpt.c_str());
double highBinfracIntegral =
((high - highBinLowEdge) / fFluxHist->GetXaxis()->GetBinWidth(maxBin)) *
fFluxHist->Integral(maxBin, maxBin, intOpt.c_str());
// If they are neighbouring bins
if ((minBin + 1) == maxBin) {
std::cout << "Get lowfrac + highfrac" << std::endl;
// Get the contained fraction of the two bin's width
return lowBinfracIntegral + highBinfracIntegral;
}
double ContainedIntegral =
fFluxHist->Integral(minBin + 1, maxBin - 1, intOpt.c_str());
// If there are filled bins between them
return lowBinfracIntegral + highBinfracIntegral + ContainedIntegral;
}
std::vector InputHandlerBase::GetFluxList(void) {
return std::vector(1, fFluxHist);
};
std::vector InputHandlerBase::GetEventList(void) {
return std::vector(1, fEventHist);
};
std::vector InputHandlerBase::GetXSecList(void) {
return std::vector(1, GetXSecHistogram());
};
FitEvent *InputHandlerBase::FirstNuisanceEvent() {
fCurrentIndex = 0;
return GetNuisanceEvent(fCurrentIndex);
};
FitEvent *InputHandlerBase::NextNuisanceEvent() {
fCurrentIndex++;
if ((fMaxEvents != -1) && (fCurrentIndex > fMaxEvents)) {
return NULL;
}
return GetNuisanceEvent(fCurrentIndex);
};
BaseFitEvt *InputHandlerBase::FirstBaseEvent() {
fCurrentIndex = 0;
return GetBaseEvent(fCurrentIndex);
};
BaseFitEvt *InputHandlerBase::NextBaseEvent() {
fCurrentIndex++;
if (jointinput and fMaxEvents != -1) {
while (fCurrentIndex < jointindexlow[jointindexswitch] ||
fCurrentIndex >= jointindexhigh[jointindexswitch]) {
jointindexswitch++;
// Loop Around
if (jointindexswitch == jointindexlow.size()) {
jointindexswitch = 0;
}
}
if (fCurrentIndex >
jointindexlow[jointindexswitch] + jointindexallowed[jointindexswitch]) {
fCurrentIndex = jointindexlow[jointindexswitch];
}
}
return GetBaseEvent(fCurrentIndex);
};
void InputHandlerBase::RegisterJointInput(std::string input, int n, TH1D *f,
TH1D *e) {
if (jointfluxinputs.size() == 0) {
jointindexswitch = 0;
fNEvents = 0;
}
// Push into individual input vectors
jointfluxinputs.push_back((TH1D *)f->Clone());
jointeventinputs.push_back((TH1D *)e->Clone());
jointindexlow.push_back(fNEvents);
jointindexhigh.push_back(fNEvents + n);
fNEvents += n;
// Add to the total flux/event hist
if (!fFluxHist)
fFluxHist = (TH1D *)f->Clone();
else
fFluxHist->Add(f);
if (!fEventHist)
fEventHist = (TH1D *)e->Clone();
else
fEventHist->Add(e);
}
void InputHandlerBase::SetupJointInputs() {
if (jointeventinputs.size() <= 1) {
jointinput = false;
} else if (jointeventinputs.size() > 1) {
jointinput = true;
jointindexswitch = 0;
}
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
if (fMaxEvents != -1 and jointeventinputs.size() > 1) {
NUIS_ABORT("Can only handle joint inputs when config MAXEVENTS = -1!");
}
if (jointeventinputs.size() > 1) {
NUIS_ERR(
WRN,
"GiBUU sample contains multiple inputs. This will only work for "
"samples that expect multi-species inputs. If this sample does, you "
"can ignore this warning.");
}
for (size_t i = 0; i < jointeventinputs.size(); i++) {
double scale = double(fNEvents) / fEventHist->Integral("width");
scale *= jointeventinputs.at(i)->Integral("width");
scale /= double(jointindexhigh[i] - jointindexlow[i]);
jointindexscale.push_back(scale);
}
fEventHist->SetNameTitle((fName + "_EVT").c_str(), (fName + "_EVT").c_str());
fFluxHist->SetNameTitle((fName + "_FLUX").c_str(), (fName + "_FLUX").c_str());
// Setup Max Events
if (fMaxEvents > 1 && fMaxEvents < fNEvents) {
if (LOG_LEVEL(SAM)) {
std::cout << "\t\t|-> Read Max Entries : " << fMaxEvents << std::endl;
}
fNEvents = fMaxEvents;
}
// Print out Status
if (LOG_LEVEL(SAM)) {
std::cout << "\t\t|-> Total Entries : " << fNEvents << std::endl
<< "\t\t|-> Event Integral : "
<< fEventHist->Integral("width") * 1.E-38 << " events/nucleon"
<< std::endl
<< "\t\t|-> Flux Integral : " << fFluxHist->Integral("width")
<< " /cm2" << std::endl
<< "\t\t|-> Event/Flux : "
<< fEventHist->Integral("width") * 1.E-38 /
fFluxHist->Integral("width")
<< " cm2/nucleon" << std::endl;
}
}
BaseFitEvt *InputHandlerBase::GetBaseEvent(const UInt_t entry) {
// Do some light processing: don't calculate the kinematics
return static_cast(GetNuisanceEvent(entry, true));
}
double InputHandlerBase::GetInputWeight(int entry) {
if (!jointinput)
return 1.0;
// Find Switch Scale
while (entry < jointindexlow[jointindexswitch] ||
entry >= jointindexhigh[jointindexswitch]) {
jointindexswitch++;
// Loop Around
if (jointindexswitch >= jointindexlow.size()) {
jointindexswitch = 0;
}
}
return jointindexscale[jointindexswitch];
};
diff --git a/src/InputHandler/NEUTInputHandler.cxx b/src/InputHandler/NEUTInputHandler.cxx
index 56b41af..7a55754 100644
--- a/src/InputHandler/NEUTInputHandler.cxx
+++ b/src/InputHandler/NEUTInputHandler.cxx
@@ -1,525 +1,520 @@
#ifdef __NEUT_ENABLED__
#include "NEUTInputHandler.h"
#include "InputUtils.h"
NEUTGeneratorInfo::~NEUTGeneratorInfo() { DeallocateParticleStack(); }
void NEUTGeneratorInfo::AddBranchesToTree(TTree *tn) {
tn->Branch("NEUTParticleN", fNEUTParticleN, "NEUTParticleN/I");
tn->Branch("NEUTParticleStatusCode", fNEUTParticleStatusCode,
"NEUTParticleStatusCode[NEUTParticleN]/I");
tn->Branch("NEUTParticleAliveCode", fNEUTParticleAliveCode,
"NEUTParticleAliveCode[NEUTParticleN]/I");
}
void NEUTGeneratorInfo::SetBranchesFromTree(TTree *tn) {
tn->SetBranchAddress("NEUTParticleN", &fNEUTParticleN);
tn->SetBranchAddress("NEUTParticleStatusCode", &fNEUTParticleStatusCode);
tn->SetBranchAddress("NEUTParticleAliveCode", &fNEUTParticleAliveCode);
}
void NEUTGeneratorInfo::AllocateParticleStack(int stacksize) {
fNEUTParticleN = 0;
fNEUTParticleStatusCode = new int[stacksize];
fNEUTParticleStatusCode = new int[stacksize];
}
void NEUTGeneratorInfo::DeallocateParticleStack() {
delete fNEUTParticleStatusCode;
delete fNEUTParticleAliveCode;
}
void NEUTGeneratorInfo::FillGeneratorInfo(NeutVect *nevent) {
Reset();
for (int i = 0; i < nevent->Npart(); i++) {
fNEUTParticleStatusCode[i] = nevent->PartInfo(i)->fStatus;
fNEUTParticleAliveCode[i] = nevent->PartInfo(i)->fIsAlive;
fNEUTParticleN++;
}
}
void NEUTGeneratorInfo::Reset() {
for (int i = 0; i < fNEUTParticleN; i++) {
fNEUTParticleStatusCode[i] = -1;
fNEUTParticleAliveCode[i] = 9;
}
fNEUTParticleN = 0;
}
NEUTInputHandler::NEUTInputHandler(std::string const &handle,
std::string const &rawinputs) {
NUIS_LOG(SAM, "Creating NEUTInputHandler : " << handle);
// Run a joint input handling
fName = handle;
// Setup the TChain
fNEUTTree = new TChain("neuttree");
fSaveExtra = FitPar::Config().GetParB("SaveExtraNEUT");
fCacheSize = FitPar::Config().GetParI("CacheSize");
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
// Loop over all inputs and grab flux, eventhist, and nevents
std::vector 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(
"NEUT 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(
(PlotUtils::GetObjectWithName(inp_file, "flux")).c_str());
TH1D *eventhist = (TH1D *)inp_file->Get(
(PlotUtils::GetObjectWithName(inp_file, "evt")).c_str());
if (!fluxhist or !eventhist) {
NUIS_ERR(FTL, "Input File Contents: " << inputs[inp_it]);
inp_file->ls();
NUIS_ABORT("NEUT FILE doesn't contain flux/xsec info. You may have to "
"regenerate your MC!");
}
// Get N Events
TTree *neuttree = (TTree *)inp_file->Get("neuttree");
if (!neuttree) {
NUIS_ERR(FTL, "neuttree not located in NEUT file: " << inputs[inp_it]);
NUIS_ABORT(
"Check your inputs, they may need to be completely regenerated!");
throw;
}
int nevents = neuttree->GetEntries();
if (nevents <= 0) {
NUIS_ABORT("Trying to a TTree with "
<< nevents << " to TChain from : " << inputs[inp_it]);
}
// Register input to form flux/event rate hists
RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist);
// Add To TChain
fNEUTTree->AddFile(inputs[inp_it].c_str());
}
// Registor all our file inputs
SetupJointInputs();
// Assign to tree
fEventType = kNEUT;
fNeutVect = NULL;
- fNEUTTree->SetBranchStatus("*", false);
- fNEUTTree->SetBranchStatus("vectorbranch", true);
- fNEUTTree->SetBranchAddress("vectorbranch", &fNeutVect);
- fNEUTTree->GetBranch("vectorbranch")->SetAutoDelete(true);
-
fNEUTTree->SetBranchAddress("vectorbranch", &fNeutVect);
fNEUTTree->GetEntry(0);
// Create Fit Event
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->SetNeutVect(fNeutVect);
if (fSaveExtra) {
fNeutInfo = new NEUTGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fNeutInfo);
}
fNUISANCEEvent->HardReset();
};
NEUTInputHandler::~NEUTInputHandler(){
// if (fNEUTTree) delete fNEUTTree;
// if (fNeutVect) delete fNeutVect;
// if (fNeutInfo) delete fNeutInfo;
};
void NEUTInputHandler::CreateCache() {
if (fCacheSize > 0) {
// fNEUTTree->SetCacheEntryRange(0, fNEvents);
fNEUTTree->AddBranchToCache("vectorbranch", 1);
fNEUTTree->SetCacheSize(fCacheSize);
}
}
void NEUTInputHandler::RemoveCache() {
// fNEUTTree->SetCacheEntryRange(0, fNEvents);
fNEUTTree->AddBranchToCache("vectorbranch", 0);
fNEUTTree->SetCacheSize(0);
}
FitEvent *NEUTInputHandler::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;
fNEUTTree->GetEntry(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
#ifdef __PROB3PP_ENABLED__
else {
UInt_t npart = fNeutVect->Npart();
for (size_t i = 0; i < npart; i++) {
NeutPart *part = fNUISANCEEvent->fNeutVect->PartInfo(i);
if ((part->fIsAlive == false) && (part->fStatus == -1) &&
std::count(PhysConst::pdg_neutrinos, PhysConst::pdg_neutrinos + 4,
part->fPID)) {
fNUISANCEEvent->probe_E = part->fP.T();
fNUISANCEEvent->probe_pdg = part->fPID;
break;
} else {
continue;
}
}
}
#endif
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
// Return event pointer
return fNUISANCEEvent;
}
// From NEUT neutclass/neutpart.h
// Bool_t fIsAlive; // Particle should be tracked or not
// ( in the detector simulator )
//
// Int_t fStatus; // Status flag of this particle
// -2: Non existing particle
// -1: Initial state particle
// 0: Normal
// 1: Decayed to the other particle
// 2: Escaped from the detector
// 3: Absorped
// 4: Charge exchanged
// 5: Pauli blocked
// 6: N/A
// 7: Produced child particles
// 8: Inelastically scattered
//
int NEUTInputHandler::GetNeutParticleStatus(NeutPart *part) {
// State
int state = kUndefinedState;
// Remove Pauli blocked events, probably just single pion events
if (part->fStatus == 5) {
state = kFSIState;
// fStatus == -1 means initial state
} else if (part->fIsAlive == false && part->fStatus == -1) {
state = kInitialState;
// NEUT has a bit of a strange convention for fIsAlive and fStatus
// combinations
// for NC and neutrino particle isAlive true/false and status 2 means
// final state particle
// for other particles in NC status 2 means it's an FSI particle
// for CC it means it was an FSI particle
} else if (part->fStatus == 2) {
// NC case is a little strange... The outgoing neutrino might be alive or
// not alive. Remaining particles with status 2 are FSI particles that
// reinteracted
if (abs(fNeutVect->Mode) > 30 &&
(abs(part->fPID) == 16 || abs(part->fPID) == 14 ||
abs(part->fPID) == 12)) {
state = kFinalState;
// The usual CC case
} else if (part->fIsAlive == true) {
state = kFSIState;
}
} else if (part->fIsAlive == true && part->fStatus == 2 &&
(abs(part->fPID) == 16 || abs(part->fPID) == 14 ||
abs(part->fPID) == 12)) {
state = kFinalState;
} else if (part->fIsAlive == true && part->fStatus == 0) {
state = kFinalState;
} else if (!part->fIsAlive &&
(part->fStatus == 1 || part->fStatus == 3 || part->fStatus == 4 ||
part->fStatus == 7 || part->fStatus == 8)) {
state = kFSIState;
// There's one hyper weird case where fStatus = -3. This apparently
// corresponds to a nucleon being ejected via pion FSI when there is "data
// available"
} else if (!part->fIsAlive && (part->fStatus == -3)) {
state = kUndefinedState;
// NC neutrino outgoing
} else if (!part->fIsAlive && part->fStatus == 0 &&
(abs(part->fPID) == 16 || abs(part->fPID) == 14 ||
abs(part->fPID) == 12)) {
state = kFinalState;
// Warn if we still find alive particles without classifying them
} else if (part->fIsAlive == true) {
NUIS_ABORT("Undefined NEUT state "
<< " Alive: " << part->fIsAlive << " Status: " << part->fStatus
<< " PDG: " << part->fPID);
// Warn if we find dead particles that we haven't classified
} else {
NUIS_ABORT("Undefined NEUT state "
<< " Alive: " << part->fIsAlive << " Status: " << part->fStatus
<< " PDG: " << part->fPID);
}
return state;
}
void NEUTInputHandler::CalcNUISANCEKinematics() {
// Reset all variables
fNUISANCEEvent->ResetEvent();
// Fill Globals
fNUISANCEEvent->Mode = fNeutVect->Mode;
fNUISANCEEvent->fEventNo = fNeutVect->EventNo;
fNUISANCEEvent->fTargetA = fNeutVect->TargetA;
fNUISANCEEvent->fTargetZ = fNeutVect->TargetZ;
fNUISANCEEvent->fTargetH = fNeutVect->TargetH;
fNUISANCEEvent->fBound = bool(fNeutVect->Ibound);
if (fNUISANCEEvent->fBound) {
fNUISANCEEvent->fTargetPDG = TargetUtils::GetTargetPDGFromZA(
fNUISANCEEvent->fTargetZ, fNUISANCEEvent->fTargetA);
} else {
fNUISANCEEvent->fTargetPDG = 1000010010;
}
// Check Particle Stack
UInt_t npart = fNeutVect->Npart();
UInt_t kmax = fNUISANCEEvent->kMaxParticles;
if (npart > kmax) {
NUIS_ERR(WRN, "NEUT has too many particles. Expanding stack.");
fNUISANCEEvent->ExpandParticleStack(npart);
}
int nprimary = fNeutVect->Nprimary();
// Fill Particle Stack
for (size_t i = 0; i < npart; i++) {
// Get Current Count
int curpart = fNUISANCEEvent->fNParticles;
// Get NEUT Particle
NeutPart *part = fNeutVect->PartInfo(i);
// State
int state = GetNeutParticleStatus(part);
// Remove Undefined
if (kRemoveUndefParticles && state == kUndefinedState)
continue;
// Remove FSI
if (kRemoveFSIParticles && state == kFSIState)
continue;
// Remove Nuclear
if (kRemoveNuclearParticles &&
(state == kNuclearInitial || state == kNuclearRemnant))
continue;
// State
fNUISANCEEvent->fParticleState[curpart] = state;
// Is the paricle associated with the primary vertex?
bool primary = false;
// NEUT events are just popped onto the stack as primary, then continues to
// be non-primary
if (i < nprimary)
primary = true;
fNUISANCEEvent->fPrimaryVertex[curpart] = primary;
// Mom
fNUISANCEEvent->fParticleMom[curpart][0] = part->fP.X();
fNUISANCEEvent->fParticleMom[curpart][1] = part->fP.Y();
fNUISANCEEvent->fParticleMom[curpart][2] = part->fP.Z();
fNUISANCEEvent->fParticleMom[curpart][3] = part->fP.T();
// PDG
fNUISANCEEvent->fParticlePDG[curpart] = part->fPID;
// Add up particle count
fNUISANCEEvent->fNParticles++;
}
// Save Extra Generator Info
if (fSaveExtra) {
fNeutInfo->FillGeneratorInfo(fNeutVect);
}
// 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 NEUTUtils::FillNeutCommons(NeutVect *nvect) {
// WARNING: This has only been implemented for a neuttree and not GENIE
// This should be kept in sync with T2KNIWGUtils::GetNIWGEvent(TTree)
// NEUT version info. Can't get it to compile properly with this yet
// neutversion_.corev = nvect->COREVer;
// neutversion_.nucev = nvect->NUCEVer;
// neutversion_.nuccv = nvect->NUCCVer;
// Documentation: See nework.h
nework_.modene = nvect->Mode;
nework_.numne = nvect->Npart();
#ifdef NEUT_COMMON_QEAV
nemdls_.mdlqeaf = nvect->QEAVForm;
#else
nemdls_.mdlqeaf = nvect->QEVForm;
#endif
nemdls_.mdlqe = nvect->QEModel;
nemdls_.mdlspi = nvect->SPIModel;
nemdls_.mdldis = nvect->DISModel;
nemdls_.mdlcoh = nvect->COHModel;
neutcoh_.necohepi = nvect->COHModel;
nemdls_.xmaqe = nvect->QEMA;
nemdls_.xmvqe = nvect->QEMV;
nemdls_.kapp = nvect->KAPPA;
// nemdls_.sccfv = SCCFVdef;
// nemdls_.sccfa = SCCFAdef;
// nemdls_.fpqe = FPQEdef;
nemdls_.xmaspi = nvect->SPIMA;
nemdls_.xmvspi = nvect->SPIMV;
nemdls_.xmares = nvect->RESMA;
nemdls_.xmvres = nvect->RESMV;
neut1pi_.xmanffres = nvect->SPIMA;
neut1pi_.xmvnffres = nvect->SPIMV;
neut1pi_.xmarsres = nvect->RESMA;
neut1pi_.xmvrsres = nvect->RESMV;
neut1pi_.neiff = nvect->SPIForm;
neut1pi_.nenrtype = nvect->SPINRType;
neut1pi_.rneca5i = nvect->SPICA5I;
neut1pi_.rnebgscl = nvect->SPIBGScale;
nemdls_.xmacoh = nvect->COHMA;
nemdls_.rad0nu = nvect->COHR0;
// nemdls_.fa1coh = nvect->COHA1err;
// nemdls_.fb1coh = nvect->COHb1err;
// neutdis_.nepdf = NEPDFdef;
// neutdis_.nebodek = NEBODEKdef;
neutcard_.nefrmflg = nvect->FrmFlg;
neutcard_.nepauflg = nvect->PauFlg;
neutcard_.nenefo16 = nvect->NefO16;
neutcard_.nemodflg = nvect->ModFlg;
// neutcard_.nenefmodl = 1;
// neutcard_.nenefmodh = 1;
// neutcard_.nenefkinh = 1;
// neutpiabs_.neabspiemit = 1;
nenupr_.iformlen = nvect->FormLen;
neutpiless_.ipilessdcy = nvect->IPilessDcy;
neutpiless_.rpilessdcy = nvect->RPilessDcy;
neutpiless_.ipilessdcy = nvect->IPilessDcy;
neutpiless_.rpilessdcy = nvect->RPilessDcy;
neffpr_.fefqe = nvect->NuceffFactorPIQE;
neffpr_.fefqeh = nvect->NuceffFactorPIQEH;
neffpr_.fefinel = nvect->NuceffFactorPIInel;
neffpr_.fefabs = nvect->NuceffFactorPIAbs;
neffpr_.fefcx = nvect->NuceffFactorPICX;
neffpr_.fefcxh = nvect->NuceffFactorPICXH;
neffpr_.fefcoh = nvect->NuceffFactorPICoh;
neffpr_.fefqehf = nvect->NuceffFactorPIQEHKin;
neffpr_.fefcxhf = nvect->NuceffFactorPICXKin;
neffpr_.fefcohf = nvect->NuceffFactorPIQELKin;
for (int i = 0; i < nework_.numne; i++) {
nework_.ipne[i] = nvect->PartInfo(i)->fPID;
nework_.pne[i][0] =
(float)nvect->PartInfo(i)->fP.X() / 1000; // VC(NE)WORK in M(G)eV
nework_.pne[i][1] =
(float)nvect->PartInfo(i)->fP.Y() / 1000; // VC(NE)WORK in M(G)eV
nework_.pne[i][2] =
(float)nvect->PartInfo(i)->fP.Z() / 1000; // VC(NE)WORK in M(G)eV
}
// fsihist.h
// neutroot fills a dummy object for events with no FSI to prevent memory leak
// when
// reading the TTree, so check for it here
if ((int)nvect->NfsiVert() ==
1) { // An event with FSI must have at least two vertices
// if (nvect->NfsiPart()!=1 || nvect->Fsiprob!=-1)
// ERR(WRN) << "T2KNeutUtils::fill_neut_commons(TTree) NfsiPart!=1 or
// Fsiprob!=-1 when NfsiVert==1" << std::endl;
fsihist_.nvert = 0;
fsihist_.nvcvert = 0;
fsihist_.fsiprob = 1;
} else { // Real FSI event
fsihist_.nvert = (int)nvect->NfsiVert();
for (int ivert = 0; ivert < fsihist_.nvert; ivert++) {
fsihist_.iflgvert[ivert] = nvect->FsiVertInfo(ivert)->fVertID;
fsihist_.posvert[ivert][0] = (float)nvect->FsiVertInfo(ivert)->fPos.X();
fsihist_.posvert[ivert][1] = (float)nvect->FsiVertInfo(ivert)->fPos.Y();
fsihist_.posvert[ivert][2] = (float)nvect->FsiVertInfo(ivert)->fPos.Z();
}
fsihist_.nvcvert = nvect->NfsiPart();
for (int ip = 0; ip < fsihist_.nvcvert; ip++) {
fsihist_.abspvert[ip] = (float)nvect->FsiPartInfo(ip)->fMomLab;
fsihist_.abstpvert[ip] = (float)nvect->FsiPartInfo(ip)->fMomNuc;
fsihist_.ipvert[ip] = nvect->FsiPartInfo(ip)->fPID;
fsihist_.iverti[ip] = nvect->FsiPartInfo(ip)->fVertStart;
fsihist_.ivertf[ip] = nvect->FsiPartInfo(ip)->fVertEnd;
fsihist_.dirvert[ip][0] = (float)nvect->FsiPartInfo(ip)->fDir.X();
fsihist_.dirvert[ip][1] = (float)nvect->FsiPartInfo(ip)->fDir.Y();
fsihist_.dirvert[ip][2] = (float)nvect->FsiPartInfo(ip)->fDir.Z();
}
fsihist_.fsiprob = nvect->Fsiprob;
}
neutcrscom_.crsx = nvect->Crsx;
neutcrscom_.crsy = nvect->Crsy;
neutcrscom_.crsz = nvect->Crsz;
neutcrscom_.crsphi = nvect->Crsphi;
neutcrscom_.crsq2 = nvect->Crsq2;
neuttarget_.numbndn = nvect->TargetA - nvect->TargetZ;
neuttarget_.numbndp = nvect->TargetZ;
neuttarget_.numfrep = nvect->TargetH;
neuttarget_.numatom = nvect->TargetA;
posinnuc_.ibound = nvect->Ibound;
// put empty nucleon FSI history (since it is not saved in the NeutVect
// format)
// Comment out as NEUT does not have the necessary proton FSI information yet
// nucleonfsihist_.nfnvert = 0;
// nucleonfsihist_.nfnstep = 0;
}
#endif