diff --git a/parameters/config.xml b/parameters/config.xml
index 328bd1d..7dc9fb0 100644
--- a/parameters/config.xml
+++ b/parameters/config.xml
@@ -1,176 +1,176 @@
<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='5'/>
<config VERBOSITY='5'/>
<!-- # ERROR goes from -->
<!-- # 0 NONE -->
<!-- # 1 FATAL -->
<!-- # 2 WARN -->
<config ERROR='2'/>
<config TRACE='1'/>
<config cores='2' />
<config spline_test_throws='50' />
<config spline_cores='6' />
<config spline_chunks='10' />
<config spline_procchunk='-1' />
<config Electron_NThetaBins='4' />
<config Electron_NEnergyBins='4' />
<config Electron_ThetaWidth='2.0' />
<config Electron_EnergyWidth='0.10' />
<config RemoveFSIParticles='0' />
<config RemoveUndefParticles='0' />
<config RemoveNuclearParticles='0'/>
<config logging.JointFCN.cxx='4'/>
<!-- # 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 input.maxevents='-1'/>
<config MAXEVENTS='-1'/>
<config input.MAXEVENTS='-1'/>
<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 input.eventmanager='1'/>
<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 -->
<config input.regen_nuwro_plots='0'/>
<!-- # DEVEL CONFIG OPTION, don't touch! -->
-<config CacheSize='5000000'/>
+<config CacheSize='0'/>
<!-- # ReWeighting Configuration Options -->
<!-- # ###################################################### -->
<!-- # Set absolute twkdial for parameters -->
<config params.setabstwk='0'/>
<!-- # Convert Dials in output statements using dial conversion card -->
<config convert_dials='0'/>
<!-- # Make RW Calculations be quiet -->
<config params.silentweighting='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 -->
<!-- # See src/FitBase/Measurement1D::Write for more info -->
<!-- #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'/>
<!-- # 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 statutils.addmcerror='0'/>
<!-- # NUISMIN Configurations -->
<!-- # ###################################################### -->
<config minimizer.maxcalls='1000000'/>
<config minimizer.maxiterations='1000000'/>
<config minimizer.tolerance='0.001'/>
<!-- # Number of events required in low stats routines -->
<config minimizer.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='250'/>
<!-- # 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_XSec_CCinc_2DEavq3_nu.hadron_cut='0'/>
<config MINERvA_CCinc_XSec_2DEavq3_nu.useq3true='0'/>
<config Modes.split_PN_NN='0'/>
<config SignalReconfigures='0'/>
</nuisance>
diff --git a/src/InputHandler/InputHandler.cxx b/src/InputHandler/InputHandler.cxx
index 6ad52e8..c6b0b79 100644
--- a/src/InputHandler/InputHandler.cxx
+++ b/src/InputHandler/InputHandler.cxx
@@ -1,275 +1,258 @@
// 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 "InputHandler.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;
};
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 minBin = fFluxHist->GetXaxis()->FindBin(low);
int maxBin = fFluxHist->GetXaxis()->FindBin(high);
return fEventHist->Integral(minBin, maxBin + 1, intOpt.c_str());
};
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() + 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;
}
// If there are filled bins between them
return lowBinfracIntegral + highBinfracIntegral +
fFluxHist->Integral(minBin + 1, maxBin - 1, intOpt.c_str());
// return fFluxHist->Integral(minBin + 1, maxBin - 1, intOpt.c_str());
}
std::vector<TH1*> InputHandlerBase::GetFluxList(void) {
return std::vector<TH1*>(1, fFluxHist);
};
std::vector<TH1*> InputHandlerBase::GetEventList(void) {
return std::vector<TH1*>(1, fEventHist);
};
std::vector<TH1*> InputHandlerBase::GetXSecList(void) {
return std::vector<TH1*>(1, GetXSecHistogram());
};
FitEvent* InputHandlerBase::FirstNuisanceEvent() {
fCurrentIndex = 0;
return GetNuisanceEvent(fCurrentIndex);
};
FitEvent* InputHandlerBase::NextNuisanceEvent() {
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 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;
}
// 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 {
jointinput= true;
}
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
for (size_t i = 0; i < jointeventinputs.size(); i++) {
TH1D* eventhist = (TH1D*) jointeventinputs.at(i)->Clone();
double scale = double(fNEvents) / fEventHist->Integral("width");
scale *= eventhist->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) {
return static_cast<BaseFitEvt*>(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 f84c4c7..f4fea23 100644
--- a/src/InputHandler/NEUTInputHandler.cxx
+++ b/src/InputHandler/NEUTInputHandler.cxx
@@ -1,429 +1,435 @@
#ifdef __NEUT_ENABLED__
#include "NEUTInputHandler.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) {
LOG(SAM) << "Creating NEUTInputHandler : " << handle << std::endl;
// 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<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()) {
- ERR(FTL) << "NEUT File IsZombie() at " << inputs[inp_it] << std::endl;
- throw;
+ THROW( "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) {
- ERR(FTL) << "NEUT FILE doesn't contain flux/xsec info" << std::endl;
- throw;
+ ERROR(FTL, "Input File Contents: " << inputs[inp_it] );
+ inp_file->ls();
+ THROW( "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) {
- ERR(FTL) << "neuttree not located in NEUT file! " << inputs[inp_it] << std::endl;
+ ERROR(FTL, "neuttree not located in NEUT file: " << inputs[inp_it]);
+ THROW("Check your inputs, they may need to be completely regenerated!");
throw;
}
int nevents = neuttree->GetEntries();
+ if (nevents <= 0) {
+ THROW("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->SetBranchAddress("vectorbranch", &fNeutVect);
// 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 entry, const bool lightweight) {
// Catch too large entries
if (entry >= (UInt_t)fNEvents) return NULL;
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fNEUTTree->GetEntry(entry);
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
// Return event pointer
return fNUISANCEEvent;
}
int NEUTInputHandler::GetNeutParticleStatus(NeutPart * part) {
// State
int state = kUndefinedState;
// fStatus == -1 means initial state
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) == 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) == 14 || abs(part->fPID) == 12)) {
state = kFinalState;
} else if (part->fIsAlive == true && part->fStatus == 0) {
state = kFinalState;
} else if (part->fIsAlive == true) {
ERR(WRN) << "Undefined NEUT state "
<< " Alive: " << part->fIsAlive << " Status: " << part->fStatus
<< " PDG: " << part->fPID << std::endl;
throw;
}
return state;
}
void NEUTInputHandler::CalcNUISANCEKinematics() {
// Reset all variables
fNUISANCEEvent->ResetEvent();
// Fill Globals
fNUISANCEEvent->fMode = fNeutVect->Mode;
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) {
ERR(FTL) << "NEUT has too many particles. Expanding stack." << std::endl;
fNUISANCEEvent->ExpandParticleStack(npart);
throw;
}
// 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;
// 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();
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();
-
- nemdls_.mdlqeaf = nvect->QEVForm;
- 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;
+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();
+
+ nemdls_.mdlqeaf = nvect->QEVForm;
+ 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
diff --git a/src/Reweight/FitWeight.cxx b/src/Reweight/FitWeight.cxx
index b976129..0ba834a 100644
--- a/src/Reweight/FitWeight.cxx
+++ b/src/Reweight/FitWeight.cxx
@@ -1,235 +1,235 @@
#include "FitWeight.h"
void FitWeight::AddRWEngine(int type) {
switch (type) {
case kNEUT:
fAllRW[type] = new NEUTWeightEngine("neutrw");
break;
case kNUWRO:
fAllRW[type] = new NuWroWeightEngine("nuwrorw");
break;
case kGENIE:
fAllRW[type] = new GENIEWeightEngine("genierw");
break;
case kNORM:
fAllRW[type] = new SampleNormEngine("normrw");
break;
case kLIKEWEIGHT:
fAllRW[type] = new LikelihoodWeightEngine("likerw");
break;
case kT2K:
fAllRW[type] = new T2KWeightEngine("t2krw");
break;
case kCUSTOM:
fAllRW[type] = new NUISANCEWeightEngine("nuisrw");
break;
case kSPLINEPARAMETER:
fAllRW[type] = new SplineWeightEngine("splinerw");
break;
case kNIWG:
fAllRW[type] = new NIWGWeightEngine("niwgrw");
break;
}
}
void FitWeight::IncludeDial(std::string name, std::string type, double val) {
// Should register the dial here.
int typeenum = Reweight::ConvDialType(type);
IncludeDial(name, typeenum, val);
}
void FitWeight::IncludeDial(std::string name, int dialtype, double val) {
// Get the dial enum
int nuisenum = Reweight::ConvDial(name, dialtype);
// Setup RW Engine Pointer
if (fAllRW.find(dialtype) == fAllRW.end()) {
AddRWEngine(dialtype);
}
WeightEngineBase* rw = fAllRW[dialtype];
// Include the dial
rw->IncludeDial(name, val);
// Set Dial Value
if (val != -9999.9) {
rw->SetDialValue(name, val);
}
// Sort Maps
fAllEnums[name] = nuisenum;
fAllValues[nuisenum] = val;
// Sort Lists
fNameList.push_back(name);
fEnumList.push_back(nuisenum);
fValueList.push_back(val);
}
void FitWeight::Reconfigure(bool silent) {
// Reconfigure all added RW engines
for (std::map<int, WeightEngineBase*>::iterator iter = fAllRW.begin();
iter != fAllRW.end(); iter++) {
(*iter).second->Reconfigure(silent);
}
}
void FitWeight::SetDialValue(std::string name, double val) {
int nuisenum = fAllEnums[name];
SetDialValue(nuisenum, val);
}
// Allow for name aswell using GlobalList to determine sample name.
void FitWeight::SetDialValue(int nuisenum, double val) {
// Conv dial type
int dialtype = int(nuisenum - (nuisenum % 1000)) / 1000;
- if (fAllRW.find(dialtype)){
+ if (fAllRW.find(dialtype) == fAllRW.end()){
THROW("Cannot find RW Engine for dialtype = " << dialtype);
}
// Get RW Engine for this dial
fAllRW[dialtype]->SetDialValue(nuisenum, val);
fAllValues[nuisenum] = val;
// Update ValueList
for (size_t i = 0; i < fEnumList.size(); i++) {
if (fEnumList[i] == nuisenum) {
fValueList[i] = val;
}
}
}
void FitWeight::SetAllDials(const double* x, int n) {
for (size_t i = 0; i < (UInt_t)n; i++) {
int rwenum = fEnumList[i];
SetDialValue(rwenum, x[i]);
}
Reconfigure();
}
double FitWeight::GetDialValue(std::string name) {
int nuisenum = fAllEnums[name];
return GetDialValue(nuisenum);
}
double FitWeight::GetDialValue(int nuisenum) {
return fAllValues[nuisenum];
}
int FitWeight::GetDialPos(std::string name) {
int rwenum = fAllEnums[name];
return GetDialPos(rwenum);
}
int FitWeight::GetDialPos(int nuisenum) {
for (size_t i = 0; i < fEnumList.size(); i++) {
if (fEnumList[i] == nuisenum) {
return i;
}
}
ERR(FTL) << "No Dial Found! " << std::endl;
throw;
return -1;
}
bool FitWeight::DialIncluded(std::string name) {
return (fAllEnums.find(name) != fAllEnums.end());
}
bool FitWeight::DialIncluded(int rwenum) {
return (fAllValues.find(rwenum) != fAllValues.end());
}
double FitWeight::CalcWeight(BaseFitEvt* evt) {
double rwweight = 1.0;
for (std::map<int, WeightEngineBase*>::iterator iter = fAllRW.begin();
iter != fAllRW.end(); iter++) {
double w = (*iter).second->CalcWeight(evt);
// LOG(FIT) << "Iter " << (*iter).second->fCalcName << " = " << w << std::endl;
rwweight *= w;
}
return rwweight;
}
void FitWeight::UpdateWeightEngine(const double* x) {
size_t count = 0;
for (std::vector<int>::iterator iter = fEnumList.begin();
iter != fEnumList.end(); iter++) {
SetDialValue( (*iter), x[count] );
count++;
}
}
void FitWeight::GetAllDials(double* x, int n) {
for (int i = 0; i < n; i++) {
x[i] = GetDialValue( fEnumList[i] );
}
}
bool FitWeight::NeedsEventReWeight(const double* x) {
bool haschange = false;
size_t count = 0;
// Compare old to new and decide if RW needed.
for (std::vector<int>::iterator iter = fEnumList.begin();
iter != fEnumList.end(); iter++) {
int nuisenum = (*iter);
int type = (nuisenum / 1000) - (nuisenum % 1000);
// Compare old to new
double oldval = GetDialValue(nuisenum);
double newval = x[count];
if (oldval != newval) {
if (fAllRW[type]->NeedsEventReWeight()) {
haschange = true;
}
}
count++;
}
return haschange;
}
double FitWeight::GetSampleNorm(std::string name) {
if (name.empty()) return 1.0;
// Find norm dial
if (fAllEnums.find(name) != fAllEnums.end()) {
if (fAllValues.find(fAllEnums[name]) != fAllValues.end()){
return fAllValues[ fAllEnums[name] ];
} else {
return 1.0;
}
} else {
return 1.0;
}
}
void FitWeight::Print() {
LOG(REC) << "Fit Weight State: " << std::endl;
for (size_t i = 0; i < fNameList.size(); i++) {
LOG(REC) << " -> Par " << i << ". " << fNameList[i] << " " << fValueList[i] << std::endl;
}
}