diff --git a/cmake/BuildDynamicSample.in b/cmake/BuildDynamicSample.in
index 2838e8c..6ed170f 100644
--- a/cmake/BuildDynamicSample.in
+++ b/cmake/BuildDynamicSample.in
@@ -1,87 +1,87 @@
# 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 .
################################################################################
#!/bin/bash
if [ ! "${1}" ] || [ ! -e ${1} ] || [ ! "${2}" ]; then
echo "[USAGE]: ${0} input.cxx outputLibName.so [classname]"
exit 1
fi
if [ ! "${3}" ]; then
CN=$(grep "class .*" $1 | sed "s/^class \([0-9a-zA-Z]\+\).*$/\1/g")
else
CN=${3}
fi
if [ ! "${CN}" ]; then
echo "[ERROR]: Couldn't find class name -- Expected to find a line like: \"class XXXX : public Measurement1D\" in \"$1\". You can also forcibly specify your classes name by passing a third argument to this script."
exit 1
fi
if [ ! -e compile.tmp ]; then
mkdir compile.tmp
fi
cat $1 > compile.tmp/$1
echo -e "static char const * SampleNames[] = {\"${CN}\"};\n"\
"static int const NSamples = 1;\n"\
"\n"\
"extern \"C\" {\n"\
"int DSF_NSamples() { return NSamples; }\n"\
"char const* DSF_GetSampleName(int i) {\n"\
" if (i < NSamples) {\n"\
" return SampleNames[i];\n"\
" }\n"\
" return 0;\n"\
"}\n"\
"MeasurementBase* DSF_GetSample(int i, void* samplekey) {\n"\
" nuiskey* sk = reinterpret_cast(samplekey);\n"\
" if (!sk) {\n"\
" return 0;\n"\
" }\n"\
"\n"\
" if (sk->GetS(\"name\") != DSF_GetSampleName(i)) {\n"\
" std::cout\n"\
" << \"[ERROR]: When instantiating dynamic sample. Samplekey named: \"\n"\
" << sk->GetS(\"name\") << \", but requested sample named: \"\n"\
" << DSF_GetSampleName(i)\n"\
" << \". It is possible that the nuiskey object is lost in translation. \"\n"\
" \"Was NUISANCE and this dynamic sample manifest built with the same \"\n"\
" \"environment and compiler?\"\n"\
" << std::endl;\n"\
" }\n"\
"\n"\
" if (i == 0) {\n"\
" return new ${CN}(*sk);\n"\
" }\n"\
" return 0;\n"\
"}\n"\
"void DSF_DestroySample(MeasurementBase* mb) { delete mb; }\n"\
"}" >> compile.tmp/$1
-echo "g++ compile.tmp/$1 -shared -o $2 -fPIC @CMAKE_CXX_FLAGS@ -I. -I@ALL_INCLUDES_STR@ -L@CMAKE_INSTALL_PREFIX@/lib -l@ALL_MODULETARGETS_STR@ @CMAKE_LINK_FLAGS@ @CMAKE_DEPENDLIB_FLAGS@"
+echo "g++ ${DYNSAMPLE_EXTRA_BUILD_FLAGS} compile.tmp/$1 -shared -o $2 -fPIC @CMAKE_CXX_FLAGS@ -I. -I@ALL_INCLUDES_STR@ -L@CMAKE_INSTALL_PREFIX@/lib -l@ALL_MODULETARGETS_STR@ @CMAKE_LINK_FLAGS@ @CMAKE_DEPENDLIB_FLAGS@"
-if ! g++ compile.tmp/$1 -shared -o $2 -fPIC @CMAKE_CXX_FLAGS@ -I. -I@ALL_INCLUDES_STR@ -L@CMAKE_INSTALL_PREFIX@/lib -l@ALL_MODULETARGETS_STR@ @CMAKE_LINK_FLAGS@ @CMAKE_DEPENDLIB_FLAGS@; then
+if ! g++ ${DYNSAMPLE_EXTRA_BUILD_FLAGS} compile.tmp/$1 -shared -o $2 -fPIC @CMAKE_CXX_FLAGS@ -I. -I@ALL_INCLUDES_STR@ -L@CMAKE_INSTALL_PREFIX@/lib -l@ALL_MODULETARGETS_STR@ @CMAKE_LINK_FLAGS@ @CMAKE_DEPENDLIB_FLAGS@; then
echo "[ERROR]: Failed to compile $1. Generated code can be found in ./compile.tmp/$1"
else
rm -r compile.tmp
echo "Successfully build: $2."
fi
diff --git a/cmake/Prob3++Setup.cmake b/cmake/Prob3++Setup.cmake
index e82e0ee..03de820 100644
--- a/cmake/Prob3++Setup.cmake
+++ b/cmake/Prob3++Setup.cmake
@@ -1,42 +1,42 @@
# 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 .
################################################################################
if(USE_PROB3PP)
LIST(APPEND EXTRA_CXX_FLAGS -D__PROB3PP_ENABLED__)
ExternalProject_Add(prob3pp
PREFIX "${CMAKE_BINARY_DIR}/Ext"
URL "http://webhome.phy.duke.edu/~raw22/public/Prob3++/Prob3++.20121225.tar.gz"
CONFIGURE_COMMAND ""
BUILD_IN_SOURCE 1
UPDATE_COMMAND ""
- BUILD_COMMAND make
+ BUILD_COMMAND CXXFLAGS=-fPIC CFLAGS=-fPIC make
INSTALL_COMMAND ""
)
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${CMAKE_BINARY_DIR}/Ext/src/prob3pp)
LIST(APPEND EXTRA_LINK_DIRS ${CMAKE_BINARY_DIR}/Ext/src/prob3pp)
LIST(APPEND EXTRA_LIBS ThreeProb_2.10)
cmessage(STATUS "Using Prob3++ 2.10")
endif()
diff --git a/src/FitBase/MeasurementBase.cxx b/src/FitBase/MeasurementBase.cxx
index facfa13..9c8c621 100644
--- a/src/FitBase/MeasurementBase.cxx
+++ b/src/FitBase/MeasurementBase.cxx
@@ -1,599 +1,597 @@
// 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 "MeasurementBase.h"
/*
Constructor/Destructors
*/
//********************************************************************
// 2nd Level Constructor (Inherits From MeasurementBase.h)
MeasurementBase::MeasurementBase(void) {
//********************************************************************
fScaleFactor = 1.0;
fMCFilled = false;
fNoData = false;
fInput = NULL;
NSignal = 0;
// Set the default values
// After-wards this gets set in SetupMeasurement
EnuMin = 0.;
EnuMax = 1.E5;
fMeasurementSpeciesType = kSingleSpeciesMeasurement;
fEventVariables = NULL;
fIsJoint = false;
fNPOT = 0xdeadbeef;
fFluxIntegralOverride = 0xdeadbeef;
fTargetVolume = 0xdeadbeef;
fTargetMaterialDensity = 0xdeadbeef;
fEvtRateScaleFactor = 0xdeadbeef;
};
void MeasurementBase::FinaliseMeasurement() {
// Used to setup default data hists, covars, etc.
}
//********************************************************************
// 2nd Level Destructor (Inherits From MeasurementBase.h)
MeasurementBase::~MeasurementBase(){
//********************************************************************
};
//********************************************************************
double MeasurementBase::TotalIntegratedFlux(std::string intOpt, double low,
double high) {
//********************************************************************
// Set Energy Limits
if (low == -9999.9) low = this->EnuMin;
if (high == -9999.9) high = this->EnuMax;
return GetInput()->TotalIntegratedFlux(low, high, intOpt);
};
//********************************************************************
double MeasurementBase::PredictedEventRate(std::string intOpt, double low,
double high) {
//********************************************************************
// Set Energy Limits
if (low == -9999.9) low = this->EnuMin;
if (high == -9999.9) high = this->EnuMax;
return GetInput()->PredictedEventRate(low, high, intOpt) * 1E-38;
};
//********************************************************************
void MeasurementBase::SetupInputs(std::string inputfile) {
//********************************************************************
// Add this infile to the global manager
if (FitPar::Config().GetParB("EventManager")) {
fInput = FitBase::AddInput(fName, inputfile);
} else {
std::vector file_descriptor =
GeneralUtils::ParseToStr(inputfile, ":");
if (file_descriptor.size() != 2) {
ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfile
<< "\". expected \"FILETYPE:file.root\"" << std::endl;
throw;
}
InputUtils::InputType inpType =
InputUtils::ParseInputType(file_descriptor[0]);
fInput = InputUtils::CreateInputHandler(fName, inpType, file_descriptor[1]);
}
fNEvents = fInput->GetNEvents();
// Expect INPUTTYPE:FileLocation(s)
std::vector file_descriptor =
GeneralUtils::ParseToStr(inputfile, ":");
if (file_descriptor.size() != 2) {
ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfile
<< "\". expected \"FILETYPE:file.root\"" << std::endl;
throw;
}
fInputType = InputUtils::ParseInputType(file_descriptor[0]);
fInputFileName = file_descriptor[1];
if (EnuMin == 0 && EnuMax == 1.E5) {
EnuMin = fInput->GetFluxHistogram()->GetBinLowEdge(1);
EnuMax = fInput->GetFluxHistogram()->GetBinLowEdge(
fInput->GetFluxHistogram()->GetNbinsX() + 1);
}
fFluxHist = fInput->GetFluxHistogram();
fEventHist = fInput->GetEventHistogram();
}
//***********************************************
int MeasurementBase::GetInputID() {
//***********************************************
return FitBase::GetInputID(fInputFileName);
}
//***********************************************
SampleSettings MeasurementBase::LoadSampleSettings(nuiskey samplekey) {
//***********************************************
SampleSettings setting = SampleSettings(samplekey);
fName = setting.GetS("name");
// Used as an initial setup function incase we need to do anything here.
LOG(SAM) << "Loading Sample : " << setting.GetName() << std::endl;
fEvtRateScaleFactor = 0xdeadbeef;
if (!fIsJoint) {
SetupInputs(setting.GetS("input"));
fNPOT = samplekey.Has("NPOT") ? samplekey.GetD("NPOT") : 1;
fFluxIntegralOverride = samplekey.Has("FluxIntegralOverride")
? samplekey.GetD("FluxIntegralOverride")
: 0xdeadbeef;
fTargetVolume = samplekey.Has("TargetVolume")
? samplekey.GetD("TargetVolume")
: 0xdeadbeef;
fTargetMaterialDensity = samplekey.Has("TargetMaterialDensity")
? samplekey.GetD("TargetMaterialDensity")
: 0xdeadbeef;
if ((fTargetVolume != 0xdeadbeef) &&
(fTargetMaterialDensity != 0xdeadbeef)) {
double TargetMass_kg = fTargetVolume * fTargetMaterialDensity;
double NNucleons = TargetMass_kg / PhysConst::mass_nucleon_kg;
double NNeutrinos =
((fFluxIntegralOverride == 0xdeadbeef) ? TotalIntegratedFlux()
: fFluxIntegralOverride) *
fNPOT;
fEvtRateScaleFactor = NNeutrinos * NNucleons;
QLOG(SAM, "\tEvent rate prediction : ");
QLOG(SAM, "\t\tTarget volume : " << fTargetVolume << " m^3");
QLOG(SAM, "\t\tTarget density : " << fTargetMaterialDensity << " kg/m^3");
QLOG(SAM, "\t\tTarget mass : " << TargetMass_kg << " kg");
QLOG(SAM, "\t\tNTarget Nucleons : " << NNucleons);
if ((fNPOT != 1)) {
QLOG(SAM, "\t\tTotal POT : " << fNPOT);
}
QLOG(SAM, "\t\tNNeutrinos : " << NNeutrinos
<< ((fNPOT != 1) ? " /cm^2" : " /POT /cm^2"));
QLOG(SAM, "\t\tXSec -> EvtRate scale factor : " << fEvtRateScaleFactor);
}
}
return setting;
}
//***********************************************
SampleSettings MeasurementBase::LoadSampleSettings(std::string name,
std::string input,
std::string type) {
//***********************************************
nuiskey samplekey = Config::CreateKey("sample");
samplekey.SetS("name", name);
samplekey.SetS("input", input);
samplekey.SetS("type", type);
return LoadSampleSettings(samplekey);
}
void MeasurementBase::FinaliseSampleSettings() {
EnuMin = fSettings.GetD("enu_min");
EnuMax = fSettings.GetD("enu_max");
}
//***********************************************
void MeasurementBase::Reconfigure() {
//***********************************************
LOG(REC) << " Reconfiguring sample " << fName << std::endl;
// Reset Histograms
ResetExtraHistograms();
AutoResetExtraTH1();
this->ResetAll();
// FitEvent* cust_event = fInput->GetEventPointer();
int fNEvents = fInput->GetNEvents();
int countwidth = (fNEvents / 5);
// MAIN EVENT LOOP
FitEvent* cust_event = fInput->FirstNuisanceEvent();
int i = 0;
int npassed = 0;
while (cust_event) {
cust_event->RWWeight = fRW->CalcWeight(cust_event);
cust_event->Weight = cust_event->RWWeight * cust_event->InputWeight;
Weight = cust_event->Weight;
// Initialize
fXVar = -999.9;
fYVar = -999.9;
fZVar = -999.9;
Signal = false;
Mode = cust_event->Mode;
// Extract Measurement Variables
this->FillEventVariables(cust_event);
Signal = this->isSignal(cust_event);
if (Signal) npassed++;
GetBox()->SetX(fXVar);
GetBox()->SetY(fYVar);
GetBox()->SetZ(fZVar);
GetBox()->SetMode(Mode);
// GetBox()->fSignal = Signal;
// Fill Histogram Values
GetBox()->FillBoxFromEvent(cust_event);
// this->FillExtraHistograms(GetBox(), Weight);
this->FillHistogramsFromBox(GetBox(), Weight);
// Print Out
if (LOG_LEVEL(REC) && countwidth > 0 && !(i % countwidth)) {
std::stringstream ss("");
ss.unsetf(std::ios_base::fixed);
ss << std::setw(7) << std::right << i << "/" << fNEvents << " events ("
<< std::setw(2) << double(i) / double(fNEvents) * 100. << std::left
<< std::setw(5) << "%) "
<< "[S,X,Y,Z,M,W] = [" << std::fixed << std::setprecision(2)
<< std::right << Signal << ", " << std::setw(5) << fXVar << ", "
<< std::setw(5) << fYVar << ", " << std::setw(5) << fYVar << ", "
<< std::setw(3) << (int)Mode << ", " << std::setw(5) << Weight << "] "
<< std::endl;
LOG(SAM) << ss.str();
}
// iterate
cust_event = fInput->NextNuisanceEvent();
i++;
}
LOG(SAM) << npassed << "/" << fNEvents << " passed selection " << std::endl;
if (npassed == 0) {
LOG(SAM) << "WARNING: NO EVENTS PASSED SELECTION!" << std::endl;
}
LOG(REC) << std::setw(10) << std::right << NSignal << "/" << fNEvents
<< " events passed selection + binning after reweight" << std::endl;
// Finalise Histograms
fMCFilled = true;
this->ConvertEventRates();
}
void MeasurementBase::FillHistogramsFromBox(MeasurementVariableBox* var,
double weight) {
fXVar = var->GetX();
fYVar = var->GetY();
fZVar = var->GetZ();
- // Signal = var->fSignal;
- // Mode = var->fMode;
Weight = weight;
fEventVariables = var;
FillHistograms();
FillExtraHistograms(var, weight);
}
void MeasurementBase::FillHistograms(double weight) {
Weight = weight * GetBox()->GetSampleWeight();
FillHistograms();
FillExtraHistograms(GetBox(), Weight);
}
MeasurementVariableBox* MeasurementBase::FillVariableBox(FitEvent* event) {
GetBox()->Reset();
Mode = event->Mode;
Weight = 1.0; // event->Weight;
this->FillEventVariables(event);
Signal = this->isSignal(event);
GetBox()->FillBoxFromEvent(event);
GetBox()->SetX(fXVar);
GetBox()->SetY(fYVar);
GetBox()->SetZ(fZVar);
GetBox()->SetMode(event->Mode);
GetBox()->SetSampleWeight(Weight);
// GetBox()->fSignal = Signal;
return GetBox();
}
MeasurementVariableBox* MeasurementBase::GetBox() {
if (!fEventVariables) fEventVariables = CreateBox();
return fEventVariables;
}
//***********************************************
void MeasurementBase::ReconfigureFast() {
//***********************************************
this->Reconfigure();
}
//***********************************************
void MeasurementBase::ConvertEventRates() {
//***********************************************
AutoScaleExtraTH1();
ScaleExtraHistograms(GetBox());
this->ScaleEvents();
double normval = fRW->GetSampleNorm(this->fName);
if (normval < 0.01 or normval > 10.0) {
ERR(WRN)
<< "Norm Value inside MeasurementBase::ConvertEventRates() looks off!"
<< std::endl;
ERR(WRN) << "It could have become out of sync with the minimizer norm list."
<< std::endl;
ERR(WRN) << "Setting it to 1.0" << std::endl;
normval = 1.0;
}
AutoNormExtraTH1(normval);
NormExtraHistograms(GetBox(), normval);
this->ApplyNormScale(normval);
}
//***********************************************
InputHandlerBase* MeasurementBase::GetInput() {
//***********************************************
if (!fInput) {
ERR(FTL) << "MeasurementBase::fInput not set. Please submit your command "
"line options and input cardfile with a bug report to: "
"nuisance@projects.hepforge.org"
<< std::endl;
throw;
}
return fInput;
};
//***********************************************
void MeasurementBase::Renormalise() {
//***********************************************
// Called when the fitter has changed a measurements normalisation but not any
// reweight dials
// Means we don't have to call the time consuming reconfigure when this
// happens.
double norm = fRW->GetDialValue(this->fName + "_norm");
if ((this->fCurrentNorm == 0.0 and norm != 0.0) or not fMCFilled) {
this->ReconfigureFast();
return;
}
if (this->fCurrentNorm == norm) return;
this->ApplyNormScale(1.0 / this->fCurrentNorm);
this->ApplyNormScale(norm);
return;
};
//***********************************************
void MeasurementBase::SetSignal(bool sig) {
//***********************************************
Signal = sig;
}
//***********************************************
void MeasurementBase::SetSignal(FitEvent* evt) {
//***********************************************
Signal = this->isSignal(evt);
}
//***********************************************
void MeasurementBase::SetWeight(double wght) {
//***********************************************
Weight = wght;
}
//***********************************************
void MeasurementBase::SetMode(int md) {
//***********************************************
Mode = md;
}
//***********************************************
std::vector MeasurementBase::GetFluxList() {
//***********************************************
return GetInput()->GetFluxList();
}
//***********************************************
std::vector MeasurementBase::GetEventRateList() {
//***********************************************
return GetInput()->GetEventList();
}
//***********************************************
std::vector MeasurementBase::GetXSecList() {
//***********************************************
return GetInput()->GetXSecList();
}
void MeasurementBase::ProcessExtraHistograms(int cmd,
MeasurementVariableBox* vars,
double weight) {
// This should be overriden if we have extra histograms!!!
// Add a flag to tell user this...
return;
}
void MeasurementBase::FillExtraHistograms(MeasurementVariableBox* vars,
double weight) {
ProcessExtraHistograms(kCMD_Fill, vars, weight);
}
void MeasurementBase::ScaleExtraHistograms(MeasurementVariableBox* vars) {
ProcessExtraHistograms(kCMD_Scale, vars, 1.0);
}
void MeasurementBase::ResetExtraHistograms() {
ProcessExtraHistograms(kCMD_Reset, NULL, 1.0);
}
void MeasurementBase::NormExtraHistograms(MeasurementVariableBox* vars,
double norm) {
ProcessExtraHistograms(kCMD_Norm, vars, norm);
}
void MeasurementBase::WriteExtraHistograms() {
ProcessExtraHistograms(kCMD_Write, NULL, 1.00);
}
void MeasurementBase::SetAutoProcessTH1(TH1* hist, int c1, int c2, int c3,
int c4, int c5) {
FakeStack* fake = new FakeStack(hist);
SetAutoProcessTH1(fake, c1, c2, c3, c4,
c5); // Need to add a destroy command!
}
void MeasurementBase::SetAutoProcess(TH1* hist, int c1, int c2, int c3, int c4,
int c5) {
FakeStack* fake = new FakeStack(hist);
SetAutoProcessTH1(fake, c1, c2, c3, c4,
c5); // Need to add a destroy command!
}
void MeasurementBase::SetAutoProcess(TGraph* g, int c1, int c2, int c3, int c4,
int c5) {
FakeStack* fake = new FakeStack(g);
SetAutoProcessTH1(fake, c1, c2, c3, c4,
c5); // Need to add a destroy command!
}
void MeasurementBase::SetAutoProcess(TF1* f, int c1, int c2, int c3, int c4,
int c5) {
FakeStack* fake = new FakeStack(f);
SetAutoProcessTH1(fake, c1, c2, c3, c4,
c5); // Need to add a destroy command!
}
void MeasurementBase::SetAutoProcess(StackBase* hist, int c1, int c2, int c3,
int c4, int c5) {
SetAutoProcessTH1(hist, c1, c2, c3, c4, c5);
}
void MeasurementBase::SetAutoProcessTH1(StackBase* hist, int c1, int c2, int c3,
int c4, int c5) {
// Set Defaults
// int ncommands = kCMD_extraplotflags;
bool autoflags[5];
autoflags[0] = false;
autoflags[1] = false;
autoflags[2] = false;
autoflags[3] = false;
autoflags[4] = false;
int givenflags[5];
givenflags[0] = c1;
givenflags[1] = c2;
givenflags[2] = c3;
givenflags[3] = c4;
givenflags[4] = c5;
fExtraTH1s[hist] = std::vector(5, 0);
// Setup a default one.
if (c1 == -1 && c2 == -1 && c3 == -1 && c4 == -1 && c5 == -1) {
fExtraTH1s[hist][kCMD_Reset] = 1;
fExtraTH1s[hist][kCMD_Scale] = 1;
fExtraTH1s[hist][kCMD_Norm] = 1;
fExtraTH1s[hist][kCMD_Write] = 1;
}
for (int i = 0; i < 5; i++) {
switch (givenflags[i]) {
// Skip over...
case -1:
break;
case kCMD_Reset:
case kCMD_Scale:
case kCMD_Norm:
case kCMD_Write:
fExtraTH1s[hist][givenflags[i]] = 1;
break;
case kCMD_Fill:
ERR(FTL) << "Can't auto fill yet!" << std::endl;
autoflags[givenflags[i]] = 1;
break;
default:
break;
}
}
// LOG(SAM) << "AutoProcessing " << hist->GetName() << std::endl;
};
void MeasurementBase::AutoFillExtraTH1() {
ERR(FTL) << "Can't auto fill yet! it's too inefficent!" << std::endl;
return;
}
void MeasurementBase::AutoResetExtraTH1() {
for (std::map >::iterator iter =
fExtraTH1s.begin();
iter != fExtraTH1s.end(); iter++) {
if (!((*iter).second)[kCMD_Reset]) continue;
(*iter).first->Reset();
}
};
void MeasurementBase::AutoScaleExtraTH1() {
for (std::map >::iterator iter =
fExtraTH1s.begin();
iter != fExtraTH1s.end(); iter++) {
if (!((*iter).second)[kCMD_Scale]) continue;
if (fIsNoWidth) {
(*iter).first->Scale(fScaleFactor);
} else {
(*iter).first->Scale(fScaleFactor, "width");
}
}
};
void MeasurementBase::AutoNormExtraTH1(double norm) {
double sfactor = 0.0;
if (norm != 0.0) sfactor = 1.0 / norm;
for (std::map >::iterator iter =
fExtraTH1s.begin();
iter != fExtraTH1s.end(); iter++) {
if (!((*iter).second)[kCMD_Norm]) continue;
(*iter).first->Scale(sfactor);
}
};
void MeasurementBase::AutoWriteExtraTH1() {
for (std::map >::iterator iter =
fExtraTH1s.begin();
iter != fExtraTH1s.end(); iter++) {
if (!(((*iter).second)[kCMD_Write])) continue;
(*iter).first->Write();
}
};
diff --git a/src/InputHandler/FitEvent.cxx b/src/InputHandler/FitEvent.cxx
index 13a4e14..61751f8 100644
--- a/src/InputHandler/FitEvent.cxx
+++ b/src/InputHandler/FitEvent.cxx
@@ -1,459 +1,429 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is pddrt 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 "FitEvent.h"
#include
#include "TObjArray.h"
FitEvent::FitEvent() {
-
fGenInfo = NULL;
kRemoveFSIParticles = true;
kRemoveUndefParticles = true;
AllocateParticleStack(400);
-
};
-
void FitEvent::AddGeneratorInfo(GeneratorInfoBase* gen) {
fGenInfo = gen;
gen->AllocateParticleStack(kMaxParticles);
}
void FitEvent::AllocateParticleStack(int stacksize) {
LOG(DEB) << "Allocating particle stack of size: " << stacksize << std::endl;
kMaxParticles = stacksize;
- fParticleList = new FitParticle*[kMaxParticles];
+ fParticleList = new FitParticle*[kMaxParticles];
- fParticleMom = new double*[kMaxParticles];
+ fParticleMom = new double*[kMaxParticles];
fParticleState = new UInt_t[kMaxParticles];
- fParticlePDG = new int[kMaxParticles];
+ fParticlePDG = new int[kMaxParticles];
- fOrigParticleMom = new double*[kMaxParticles];
+ fOrigParticleMom = new double*[kMaxParticles];
fOrigParticleState = new UInt_t[kMaxParticles];
- fOrigParticlePDG = new int[kMaxParticles];
+ fOrigParticlePDG = new int[kMaxParticles];
for (size_t i = 0; i < kMaxParticles; i++) {
- fParticleList[i] = NULL;
- fParticleMom[i] = new double[4];
+ fParticleList[i] = NULL;
+ fParticleMom[i] = new double[4];
fOrigParticleMom[i] = new double[4];
}
if (fGenInfo) fGenInfo->AllocateParticleStack(kMaxParticles);
-
}
void FitEvent::ExpandParticleStack(int stacksize) {
DeallocateParticleStack();
AllocateParticleStack(stacksize);
}
void FitEvent::DeallocateParticleStack() {
-
for (size_t i = 0; i < kMaxParticles; i++) {
if (fParticleList[i]) delete fParticleList[i];
delete fParticleMom[i];
delete fOrigParticleMom[i];
}
delete fParticleMom;
delete fOrigParticleMom;
delete fParticleList;
delete fParticleState;
delete fParticlePDG;
delete fOrigParticleState;
delete fOrigParticlePDG;
- if (fGenInfo) fGenInfo -> DeallocateParticleStack();
+ if (fGenInfo) fGenInfo->DeallocateParticleStack();
kMaxParticles = 0;
-
}
void FitEvent::ClearFitParticles() {
for (size_t i = 0; i < kMaxParticles; i++) {
fParticleList[i] = NULL;
}
}
void FitEvent::FreeFitParticles() {
for (size_t i = 0; i < kMaxParticles; i++) {
FitParticle* fp = fParticleList[i];
if (fp) delete fp;
fParticleList[i] = NULL;
}
}
void FitEvent::ResetParticleList() {
for (unsigned int i = 0; i < kMaxParticles; i++) {
FitParticle* fp = fParticleList[i];
if (fp) delete fp;
fParticleList[i] = NULL;
}
}
void FitEvent::HardReset() {
for (unsigned int i = 0; i < kMaxParticles; i++) {
fParticleList[i] = NULL;
}
}
void FitEvent::ResetEvent() {
-
- fMode = 9999;
Mode = 9999;
fEventNo = -1;
fTotCrs = -1.0;
fTargetA = -1;
fTargetZ = -1;
fTargetH = -1;
fBound = false;
fNParticles = 0;
if (fGenInfo) fGenInfo->Reset();
for (unsigned int i = 0; i < kMaxParticles; i++) {
if (fParticleList[i]) delete fParticleList[i];
fParticleList[i] = NULL;
continue;
fParticlePDG[i] = 0;
fParticleState[i] = kUndefinedState;
fParticleMom[i][0] = 0.0;
fParticleMom[i][1] = 0.0;
fParticleMom[i][2] = 0.0;
fParticleMom[i][3] = 0.0;
fOrigParticlePDG[i] = 0;
fOrigParticleState[i] = kUndefinedState;
fOrigParticleMom[i][0] = 0.0;
fOrigParticleMom[i][1] = 0.0;
fOrigParticleMom[i][2] = 0.0;
fOrigParticleMom[i][3] = 0.0;
}
-
}
void FitEvent::OrderStack() {
-
// Copy current stack
int npart = fNParticles;
for (int i = 0; i < npart; i++) {
- fOrigParticlePDG[i] = fParticlePDG[i];
- fOrigParticleState[i] = fParticleState[i];
+ fOrigParticlePDG[i] = fParticlePDG[i];
+ fOrigParticleState[i] = fParticleState[i];
fOrigParticleMom[i][0] = fParticleMom[i][0];
fOrigParticleMom[i][1] = fParticleMom[i][1];
fOrigParticleMom[i][2] = fParticleMom[i][2];
fOrigParticleMom[i][3] = fParticleMom[i][3];
}
// Now run loops for each particle
fNParticles = 0;
int stateorder[6] = {kInitialState, kFinalState, kFSIState,
- kNuclearInitial, kNuclearRemnant, kUndefinedState
- };
+ kNuclearInitial, kNuclearRemnant, kUndefinedState};
for (int s = 0; s < 6; s++) {
for (int i = 0; i < npart; i++) {
if ((UInt_t)fOrigParticleState[i] != (UInt_t)stateorder[s]) continue;
- fParticlePDG[fNParticles] = fOrigParticlePDG[i];
- fParticleState[fNParticles] = fOrigParticleState[i];
+ fParticlePDG[fNParticles] = fOrigParticlePDG[i];
+ fParticleState[fNParticles] = fOrigParticleState[i];
fParticleMom[fNParticles][0] = fOrigParticleMom[i][0];
fParticleMom[fNParticles][1] = fOrigParticleMom[i][1];
fParticleMom[fNParticles][2] = fOrigParticleMom[i][2];
fParticleMom[fNParticles][3] = fOrigParticleMom[i][3];
fNParticles++;
}
}
if (LOG_LEVEL(DEB)) {
LOG(DEB) << "Ordered stack" << std::endl;
for (int i = 0; i < fNParticles; i++) {
LOG(DEB) << "Particle " << i << ". " << fParticlePDG[i] << " "
<< fParticleMom[i][0] << " " << fParticleMom[i][1] << " "
<< fParticleMom[i][2] << " " << fParticleMom[i][3] << " "
<< fParticleState[i] << std::endl;
}
}
if (fNParticles != npart) {
ERR(FTL) << "Dropped some particles when ordering the stack!" << std::endl;
}
return;
}
-
void FitEvent::Print() {
-
if (LOG_LEVEL(FIT)) {
LOG(FIT) << "FITEvent print" << std::endl;
- LOG(FIT) << "Mode: " << fMode << std::endl;
+ LOG(FIT) << "Mode: " << Mode << std::endl;
LOG(FIT) << "Particles: " << fNParticles << std::endl;
LOG(FIT) << " -> Particle Stack " << std::endl;
for (int i = 0; i < fNParticles; i++) {
LOG(FIT) << " -> -> " << i << ". " << fParticlePDG[i] << " "
<< fParticleState[i] << " "
<< " Mom(" << fParticleMom[i][0] << ", " << fParticleMom[i][1]
- << ", " << fParticleMom[i][2] << ", " << fParticleMom[i][3] << ")."
- << std::endl;
+ << ", " << fParticleMom[i][2] << ", " << fParticleMom[i][3]
+ << ")." << std::endl;
}
}
return;
}
/* Read/Write own event class */
void FitEvent::SetBranchAddress(TChain* tn) {
-
- tn->SetBranchAddress("Mode", &fMode);
- tn->SetBranchAddress("Mode", &Mode);
+ tn->SetBranchAddress("Mode", &Mode);
tn->SetBranchAddress("EventNo", &fEventNo);
- tn->SetBranchAddress("TotCrs", &fTotCrs);
+ tn->SetBranchAddress("TotCrs", &fTotCrs);
tn->SetBranchAddress("TargetA", &fTargetA);
tn->SetBranchAddress("TargetH", &fTargetH);
- tn->SetBranchAddress("Bound", &fBound);
+ tn->SetBranchAddress("Bound", &fBound);
- tn->SetBranchAddress("RWWeight", &SavedRWWeight);
+ tn->SetBranchAddress("RWWeight", &SavedRWWeight);
tn->SetBranchAddress("InputWeight", &InputWeight);
-
-
-
- // This has to be setup by the handler now :(
- // tn->SetBranchAddress("NParticles", &fNParticles);
- // tn->SetBranchAddress("ParticleState", &fParticleState);
- // tn->SetBranchAddress("ParticlePDG", &fParticlePDG);
- // tn->SetBranchAddress("ParticleMom", &fParticleMom);
-
}
void FitEvent::AddBranchesToTree(TTree* tn) {
-
tn->Branch("Mode", &Mode, "Mode/I");
tn->Branch("EventNo", &fEventNo, "EventNo/i");
- tn->Branch("TotCrs", &fTotCrs, "TotCrs/D");
+ tn->Branch("TotCrs", &fTotCrs, "TotCrs/D");
tn->Branch("TargetA", &fTargetA, "TargetA/I");
tn->Branch("TargetH", &fTargetH, "TargetH/I");
- tn->Branch("Bound", &fBound, "Bound/O");
+ tn->Branch("Bound", &fBound, "Bound/O");
- tn->Branch("RWWeight", &RWWeight, "RWWeight/D");
+ tn->Branch("RWWeight", &RWWeight, "RWWeight/D");
tn->Branch("InputWeight", &InputWeight, "InputWeight/D");
- tn->Branch("NParticles", &fNParticles, "NParticles/I");
- tn->Branch("ParticleState", fOrigParticleState, "ParticleState[NParticles]/i");
- tn->Branch("ParticlePDG", fOrigParticlePDG, "ParticlePDG[NParticles]/I");
- tn->Branch("ParticleMom", fOrigParticleMom, "ParticleMom[NParticles][4]/D");
-
+ tn->Branch("NParticles", &fNParticles, "NParticles/I");
+ tn->Branch("ParticleState", fOrigParticleState,
+ "ParticleState[NParticles]/i");
+ tn->Branch("ParticlePDG", fOrigParticlePDG, "ParticlePDG[NParticles]/I");
+ tn->Branch("ParticleMom", fOrigParticleMom, "ParticleMom[NParticles][4]/D");
}
-
// ------- EVENT ACCESS FUNCTION --------- //
-TLorentzVector FitEvent::GetParticleP4 (int index) const {
+TLorentzVector FitEvent::GetParticleP4(int index) const {
if (index == -1 or index >= fNParticles) return TLorentzVector();
- return TLorentzVector( fParticleMom[index][0],
- fParticleMom[index][1],
- fParticleMom[index][2],
- fParticleMom[index][3] );
+ return TLorentzVector(fParticleMom[index][0], fParticleMom[index][1],
+ fParticleMom[index][2], fParticleMom[index][3]);
}
-TVector3 FitEvent::GetParticleP3 (int index) const {
+TVector3 FitEvent::GetParticleP3(int index) const {
if (index == -1 or index >= fNParticles) return TVector3();
- return TVector3( fParticleMom[index][0],
- fParticleMom[index][1],
- fParticleMom[index][2] );
+ return TVector3(fParticleMom[index][0], fParticleMom[index][1],
+ fParticleMom[index][2]);
}
-double FitEvent::GetParticleMom (int index) const {
+double FitEvent::GetParticleMom(int index) const {
if (index == -1 or index >= fNParticles) return 0.0;
return sqrt(fParticleMom[index][0] * fParticleMom[index][0] +
fParticleMom[index][1] * fParticleMom[index][1] +
fParticleMom[index][2] * fParticleMom[index][2]);
}
-double FitEvent::GetParticleMom2 (int index) const {
+double FitEvent::GetParticleMom2(int index) const {
if (index == -1 or index >= fNParticles) return 0.0;
return fabs((fParticleMom[index][0] * fParticleMom[index][0] +
fParticleMom[index][1] * fParticleMom[index][1] +
fParticleMom[index][2] * fParticleMom[index][2]));
}
-double FitEvent::GetParticleE (int index) const {
+double FitEvent::GetParticleE(int index) const {
if (index == -1 or index >= fNParticles) return 0.0;
return fParticleMom[index][3];
}
-int FitEvent::GetParticleState (int index) const {
+int FitEvent::GetParticleState(int index) const {
if (index == -1 or index >= fNParticles) return kUndefinedState;
return (fParticleState[index]);
}
-int FitEvent::GetParticlePDG (int index) const {
+int FitEvent::GetParticlePDG(int index) const {
if (index == -1 or index >= fNParticles) return 0;
return (fParticlePDG[index]);
}
-FitParticle* FitEvent::GetParticle (int const i) {
-
+FitParticle* FitEvent::GetParticle(int const i) {
// Check Valid Index
- if (i == -1){
+ if (i == -1) {
return NULL;
}
// Check Valid
if (i > fNParticles) {
ERR(FTL) << "Requesting particle beyond stack!" << std::endl
<< "i = " << i << " N = " << fNParticles << std::endl
- << "Mode = " << fMode << std::endl;
+ << "Mode = " << Mode << std::endl;
throw;
}
if (!fParticleList[i]) {
/*
std::cout << "Creating particle with values i " << i << " ";
- std::cout << fParticleMom[i][0] << " " << fParticleMom[i][1] << " " << fParticleMom[i][2] << " " << fParticleMom[i][3] << " ";
+ std::cout << fParticleMom[i][0] << " " << fParticleMom[i][1] << " " <<
+ fParticleMom[i][2] << " " << fParticleMom[i][3] << " ";
std::cout << fParticlePDG[i] << " " << fParticleState[i] << std::endl;
*/
fParticleList[i] = new FitParticle(fParticleMom[i][0], fParticleMom[i][1],
fParticleMom[i][2], fParticleMom[i][3],
fParticlePDG[i], fParticleState[i]);
} else {
/*
std::cout << "Filling particle with values i " << i << " ";
- std::cout << fParticleMom[i][0] << " " << fParticleMom[i][1] << " " << fParticleMom[i][2] << " " << fParticleMom[i][3] << " ";
+ std::cout << fParticleMom[i][0] << " " << fParticleMom[i][1] << " " <<
+ fParticleMom[i][2] << " " << fParticleMom[i][3] << " ";
std::cout << fParticlePDG[i] << " "<< fParticleState[i] <SetValues(fParticleMom[i][0], fParticleMom[i][1],
fParticleMom[i][2], fParticleMom[i][3],
fParticlePDG[i], fParticleState[i]);
-
}
return fParticleList[i];
}
bool FitEvent::HasParticle(int const pdg, int const state) const {
bool found = false;
for (int i = 0; i < fNParticles; i++) {
if (state != -1 && fParticleState[i] != (uint)state) continue;
if (fParticlePDG[i] == pdg) found = true;
}
return found;
}
int FitEvent::NumParticle(int const pdg, int const state) const {
int nfound = 0;
for (int i = 0; i < fNParticles; i++) {
if (state != -1 and fParticleState[i] != (uint)state) continue;
if (pdg == 0 or fParticlePDG[i] == pdg) nfound += 1;
}
return nfound;
}
-std::vector FitEvent::GetAllParticleIndices (int const pdg, int const state) const {
+std::vector FitEvent::GetAllParticleIndices(int const pdg,
+ int const state) const {
std::vector indexlist;
for (int i = 0; i < fNParticles; i++) {
if (state != -1 and fParticleState[i] != (uint)state) continue;
if (pdg == 0 or fParticlePDG[i] == pdg) {
indexlist.push_back(i);
}
}
return indexlist;
}
-std::vector FitEvent::GetAllParticle(int const pdg, int const state) {
+std::vector FitEvent::GetAllParticle(int const pdg,
+ int const state) {
std::vector indexlist = GetAllParticleIndices(pdg, state);
std::vector plist;
for (std::vector::iterator iter = indexlist.begin();
iter != indexlist.end(); iter++) {
- plist.push_back( GetParticle((*iter)) );
+ plist.push_back(GetParticle((*iter)));
}
return plist;
}
-int FitEvent::GetHMParticleIndex (int const pdg, int const state) const {
+int FitEvent::GetHMParticleIndex(int const pdg, int const state) const {
double maxmom2 = -9999999.9;
int maxind = -1;
for (int i = 0; i < fNParticles; i++) {
if (state != -1 and fParticleState[i] != (uint)state) continue;
if (pdg == 0 or fParticlePDG[i] == pdg) {
-
double newmom2 = GetParticleMom2(i);
if (newmom2 > maxmom2) {
maxind = i;
maxmom2 = newmom2;
}
}
}
return maxind;
}
-int FitEvent::GetBeamNeutrinoIndex (void) const {
- for (int i = 0; i < fNParticles; i++){
+int FitEvent::GetBeamNeutrinoIndex(void) const {
+ for (int i = 0; i < fNParticles; i++) {
if (fParticleState[i] != kInitialState) continue;
int pdg = abs(fParticlePDG[i]);
- if (pdg == 12 or pdg == 14 or pdg == 16){
+ if (pdg == 12 or pdg == 14 or pdg == 16) {
return i;
}
}
return 0;
}
-
-int FitEvent::GetBeamElectronIndex (void) const {
- return GetHMISParticleIndex( 11 );
+int FitEvent::GetBeamElectronIndex(void) const {
+ return GetHMISParticleIndex(11);
}
-int FitEvent::GetBeamPionIndex (void) const {
- return GetHMISParticleIndex( PhysConst::pdg_pions );
+int FitEvent::GetBeamPionIndex(void) const {
+ return GetHMISParticleIndex(PhysConst::pdg_pions);
}
int FitEvent::NumFSMesons() {
int nMesons = 0;
for (int i = 0; i < fNParticles; i++) {
if (fParticleState[i] != kFinalState) continue;
if (abs(fParticlePDG[i]) >= 111 && abs(fParticlePDG[i]) <= 557)
nMesons += 1;
}
return nMesons;
}
-int FitEvent::NumFSLeptons(void) const{
-
+int FitEvent::NumFSLeptons(void) const {
int nLeptons = 0;
for (int i = 0; i < fNParticles; i++) {
if (fParticleState[i] != kFinalState) continue;
if (abs(fParticlePDG[i]) == 11 || abs(fParticlePDG[i]) == 13 ||
- abs(fParticlePDG[i]) == 15)
+ abs(fParticlePDG[i]) == 15)
nLeptons += 1;
}
return nLeptons;
}
diff --git a/src/InputHandler/FitEvent.h b/src/InputHandler/FitEvent.h
index 07653b3..ddc5f99 100644
--- a/src/InputHandler/FitEvent.h
+++ b/src/InputHandler/FitEvent.h
@@ -1,640 +1,636 @@
// 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 .
*******************************************************************************/
#ifndef FITEVENT2_H_SEEN
#define FITEVENT2_H_SEEN
/*!
* \addtogroup InputHandler
* @{
*/
#include
#include
#include
#include "FitParticle.h"
#include "TLorentzVector.h"
#include "TSpline.h"
#include "BaseFitEvt.h"
#include "FitLogger.h"
#include "TArrayD.h"
#include "TTree.h"
#include "TChain.h"
#include "PhysConst.h"
/// Common container for event particles
class FitEvent : public BaseFitEvt {
public:
FitEvent();
~FitEvent() {};
void FreeFitParticles();
void ClearFitParticles();
void ResetEvent(void);
void ResetParticleList(void);
void HardReset();
void OrderStack();
void SetBranchAddress(TChain* tn);
void AddBranchesToTree(TTree* tn);
void Print();
void PrintChris();
void DeallocateParticleStack();
void AllocateParticleStack(int stacksize);
void ExpandParticleStack(int stacksize);
void AddGeneratorInfo(GeneratorInfoBase* gen);
// ---- HELPER/ACCESS FUNCTIONS ---- //
/// Return True Interaction ID
- inline int GetMode (void) const { return fMode; };
+ inline int GetMode (void) const { return Mode; };
/// Return Target Atomic Number
inline int GetTargetA (void) const { return fTargetA; };
/// Return Target Nuclear Charge
inline int GetTargetZ (void) const { return fTargetZ; };
/// Get Event Total Cross-section
inline int GetTotCrs (void) const { return fTotCrs; };
/// Is Event Charged Current?
- inline bool IsCC() const { return (abs(fMode) <= 30); };
+ inline bool IsCC() const { return (abs(Mode) <= 30); };
/// Is Event Neutral Current?
- inline bool IsNC() const { return (abs(fMode) > 30); };
+ inline bool IsNC() const { return (abs(Mode) > 30); };
/// Return Particle 4-momentum for given index in particle stack
TLorentzVector GetParticleP4 (int index) const;
/// Return Particle 3-momentum for given index in particle stack
TVector3 GetParticleP3 (int index) const;
/// Return Particle absolute momentum for given index in particle stack
double GetParticleMom (int index) const;
/// Return Particle absolute momentum-squared for given index in particle stack
double GetParticleMom2 (int index) const;
/// Return Particle energy for given index in particle stack
double GetParticleE (int index) const;
/// Return Particle State for given index in particle stack
int GetParticleState (int index) const;
/// Return Particle PDG for given index in particle stack
int GetParticlePDG (int index) const;
/// Allows the removal of KE up to total KE.
inline void RemoveKE(int index, double KE){
FitParticle *fp = GetParticle(index);
double mass = fp->M();
double oKE = fp->KE();
double nE = mass + (oKE - KE);
if(nE < mass){ // Can't take more KE than it has
nE = mass;
}
double n3Mom = sqrt(nE*nE - mass*mass);
TVector3 np3 = fp->P3().Unit()*n3Mom;
fParticleMom[index][0] = np3[0];
fParticleMom[index][1] = np3[1];
fParticleMom[index][2] = np3[2];
fParticleMom[index][3] = nE;
}
/// Allows the removal of KE up to total KE.
inline void GiveKE(int index, double KE){
RemoveKE(index,-KE);
}
/// Return Particle for given index in particle stack
FitParticle* GetParticle(int const index);
/// Get Total Number of Particles in stack
inline uint NParticles (void) const { return fNParticles; };
/// Check if event contains a particle given a pdg and state.
/// If no state is passed all states are considered.
bool HasParticle (int const pdg = 0, int const state = -1) const ;
template
inline bool HasParticle(int const (&pdgs)[N], int const state = -1) const {
for (size_t i = 0; i < N; i++) {
if (HasParticle( pdgs[i], state )) {
return false;
}
}
return false;
}
/// Get total number of particles given a pdg and state.
/// If no state is passed all states are considered.
int NumParticle (int const pdg = 0, int const state = -1) const;
template
inline int NumParticle(int const (&pdgs)[N], int const state = -1) const {
int ncount = 0;
for (size_t i = 0; i < N; i++) {
ncount += NumParticle( pdgs[i], state );
}
return ncount;
}
/// Return a vector of particle indices that can be used with 'GetParticle'
/// functions given a particle pdg and state.
/// If no state is passed all states are considered.
std::vector GetAllParticleIndices (int const pdg = -1, int const state = -1) const;
template
inline std::vector GetAllParticleIndices(int const (&pdgs)[N], const int state = -1) const {
std::vector plist;
for (size_t i = 0; i < N; i++) {
std::vector plisttemp = GetAllParticleIndices(pdgs[i], state);
plist.insert( plist.end(), plisttemp.begin(), plisttemp.end() );
}
return plist;
}
/// Return a vector of FitParticles given a particle pdg and state.
/// This is memory intensive and slow than GetAllParticleIndices,
/// but is slightly easier to use.
std::vector GetAllParticle (int const pdg = -1, int const state = -1);
template
inline std::vector GetAllParticle(int const (&pdgs)[N], int const state = -1) {
std::vector plist;
for (size_t i = 0; i < N; i++) {
std::vector plisttemp = GetAllParticle(pdgs[i], state);
plist.insert( plist.end(), plisttemp.begin(), plisttemp.end() );
}
return plist;
}
inline std::vector GetAllNuElectronIndices (void) { return GetAllParticleIndices(12); };
inline std::vector GetAllNuMuonIndices (void) { return GetAllParticleIndices(14); };
inline std::vector GetAllNuTauIndices (void) { return GetAllParticleIndices(16); };
inline std::vector GetAllElectronIndices (void) { return GetAllParticleIndices(11); };
inline std::vector GetAllMuonIndices (void) { return GetAllParticleIndices(13); };
inline std::vector GetAllTauIndices (void) { return GetAllParticleIndices(15); };
inline std::vector GetAllProtonIndices (void) { return GetAllParticleIndices(2212); };
inline std::vector GetAllNeutronIndices (void) { return GetAllParticleIndices(2112); };
inline std::vector GetAllPiZeroIndices (void) { return GetAllParticleIndices(111); };
inline std::vector GetAllPiPlusIndices (void) { return GetAllParticleIndices(211); };
inline std::vector GetAllPiMinusIndices (void) { return GetAllParticleIndices(-211); };
inline std::vector GetAllPhotonIndices (void) { return GetAllParticleIndices(22); };
inline std::vector GetAllNuElectron (void) { return GetAllParticle(12); };
inline std::vector GetAllNuMuon (void) { return GetAllParticle(14); };
inline std::vector GetAllNuTau (void) { return GetAllParticle(16); };
inline std::vector GetAllElectron (void) { return GetAllParticle(11); };
inline std::vector GetAllMuon (void) { return GetAllParticle(13); };
inline std::vector GetAllTau (void) { return GetAllParticle(15); };
inline std::vector GetAllProton (void) { return GetAllParticle(2212); };
inline std::vector GetAllNeutron (void) { return GetAllParticle(2112); };
inline std::vector GetAllPiZero (void) { return GetAllParticle(111); };
inline std::vector GetAllPiPlus (void) { return GetAllParticle(211); };
inline std::vector GetAllPiMinus (void) { return GetAllParticle(-211); };
inline std::vector GetAllPhoton (void) { return GetAllParticle(22); };
// --- Highest Momentum Search Functions --- //
/// Returns the Index of the highest momentum particle given a pdg and state.
/// If no state is given all states are considered, but that will just return the
/// momentum of the beam in most cases so is not advised.
int GetHMParticleIndex (int const pdg = 0, int const state = -1) const;
template
inline int GetHMParticleIndex (int const (&pdgs)[N], int const state = -1) const {
double mom = -999.9;
int rtnindex = -1;
for (size_t i = 0; i < N; ++i) {
// Use ParticleMom as doesn't require Particle Mem alloc
int pindex = GetHMParticleIndex(pdgs[i], state);
if (pindex != -1){
double momnew = GetParticleMom2(pindex);
if (momnew > mom) {
rtnindex = pindex;
mom = momnew;
}
}
}
return rtnindex;
};
/// Returns the highest momentum particle given a pdg and state.
/// If no state is given all states are considered, but that will just return the
/// momentum of the beam in most cases so is not advised.
inline FitParticle* GetHMParticle(int const pdg = 0, int const state = -1) {
return GetParticle( GetHMParticleIndex(pdg, state) );
}
template
inline FitParticle* GetHMParticle(int const (&pdgs)[N], int const state) {
return GetParticle(GetHMParticleIndex(pdgs, state));
};
// ---- Initial State Helpers --- //
/// Checks the event has a particle of a given pdg in the initial state.
inline bool HasISParticle(int const pdg) const {
return HasParticle(pdg, kInitialState);
};
template
inline bool HasISParticle(int const (&pdgs)[N]) const {
return HasParticle(pdgs, kInitialState);
};
/// Returns the number of particles with a given pdg in the initial state.
inline int NumISParticle(int const pdg = 0) const {
return NumParticle(pdg, kInitialState);
};
template
inline int NumISParticle(int const (&pdgs)[N]) const {
return NumParticle(pdgs, kInitialState);
};
/// Returns a list of indices for all particles with a given pdg
/// in the initial state. These can be used with the 'GetParticle' functions.
inline std::vector GetAllISParticleIndices(int const pdg = -1) const {
return GetAllParticleIndices(pdg, kInitialState);
};
template
inline std::vector GetAllISParticleIndices(int const (&pdgs)[N]) const {
return GetAllParticleIndices(pdgs, kInitialState);
};
/// Returns a list of particles with a given pdg in the initial state.
/// This function is more memory intensive and slower than the Indices function.
inline std::vector GetAllISParticle(int const pdg = -1) {
return GetAllParticle(pdg, kInitialState);
};
template
inline std::vector GetAllISParticle(int const (&pdgs)[N]) {
return GetAllParticle(pdgs, kInitialState);
};
/// Returns the highest momentum particle with a given pdg in the initial state.
inline FitParticle* GetHMISParticle(int const pdg) {
return GetHMParticle(pdg, kInitialState);
};
template
inline FitParticle* GetHMISParticle(int const (&pdgs)[N]) {
return GetHMParticle(pdgs, kInitialState);
};
/// Returns the highest momentum particle index with a given pdg in the initial state.
inline int GetHMISParticleIndex(int const pdg) const {
return GetHMParticleIndex(pdg, kInitialState);
};
template
inline int GetHMISParticleIndex(int const (&pdgs)[N]) const {
return GetHMParticleIndex(pdgs, kInitialState);
};
inline bool HasISNuElectron (void) const { return HasISParticle(12); };
inline bool HasISNuMuon (void) const { return HasISParticle(14); };
inline bool HasISNuTau (void) const { return HasISParticle(16); };
inline bool HasISElectron (void) const { return HasISParticle(11); };
inline bool HasISMuon (void) const { return HasISParticle(13); };
inline bool HasISTau (void) const { return HasISParticle(15); };
inline bool HasISProton (void) const { return HasISParticle(2212); };
inline bool HasISNeutron (void) const { return HasISParticle(2112); };
inline bool HasISPiZero (void) const { return HasISParticle(111); };
inline bool HasISPiPlus (void) const { return HasISParticle(211); };
inline bool HasISPiMinus (void) const { return HasISParticle(-211); };
inline bool HasISPhoton (void) const { return HasISParticle(22); };
inline bool HasISLeptons (void) const { return HasISParticle(PhysConst::pdg_leptons); };
inline bool HasISPions (void) const { return HasISParticle(PhysConst::pdg_pions); };
inline bool HasISChargePions (void) const { return HasISParticle(PhysConst::pdg_charged_pions); };
inline int NumISNuElectron (void) const { return NumISParticle(12); };
inline int NumISNuMuon (void) const { return NumISParticle(14); };
inline int NumISNuTau (void) const { return NumISParticle(16); };
inline int NumISElectron (void) const { return NumISParticle(11); };
inline int NumISMuon (void) const { return NumISParticle(13); };
inline int NumISTau (void) const { return NumISParticle(15); };
inline int NumISProton (void) const { return NumISParticle(2212); };
inline int NumISNeutron (void) const { return NumISParticle(2112); };
inline int NumISPiZero (void) const { return NumISParticle(111); };
inline int NumISPiPlus (void) const { return NumISParticle(211); };
inline int NumISPiMinus (void) const { return NumISParticle(-211); };
inline int NumISPhoton (void) const { return NumISParticle(22); };
inline int NumISLeptons (void) const { return NumISParticle(PhysConst::pdg_leptons); };
inline int NumISPions (void) const { return NumISParticle(PhysConst::pdg_pions); };
inline int NumISChargePions (void) const { return NumISParticle(PhysConst::pdg_charged_pions); };
inline std::vector GetAllISNuElectronIndices (void) const { return GetAllISParticleIndices(12); };
inline std::vector GetAllISNuMuonIndices (void) const { return GetAllISParticleIndices(14); };
inline std::vector GetAllISNuTauIndices (void) const { return GetAllISParticleIndices(16); };
inline std::vector GetAllISElectronIndices (void) const { return GetAllISParticleIndices(11); };
inline std::vector GetAllISMuonIndices (void) const { return GetAllISParticleIndices(13); };
inline std::vector GetAllISTauIndices (void) const { return GetAllISParticleIndices(15); };
inline std::vector GetAllISProtonIndices (void) const { return GetAllISParticleIndices(2212); };
inline std::vector GetAllISNeutronIndices (void) const { return GetAllISParticleIndices(2112); };
inline std::vector GetAllISPiZeroIndices (void) const { return GetAllISParticleIndices(111); };
inline std::vector GetAllISPiPlusIndices (void) const { return GetAllISParticleIndices(211); };
inline std::vector GetAllISPiMinusIndices (void) const { return GetAllISParticleIndices(-211); };
inline std::vector GetAllISPhotonIndices (void) const { return GetAllISParticleIndices(22); };
inline std::vector GetAllISLeptonsIndices (void) const { return GetAllISParticleIndices(PhysConst::pdg_leptons); };
inline std::vector GetAllISPionsIndices (void) const { return GetAllISParticleIndices(PhysConst::pdg_pions); };
inline std::vector GetAllISChargePionsIndices(void) const { return GetAllISParticleIndices(PhysConst::pdg_charged_pions); };
inline std::vector GetAllISNuElectron (void) { return GetAllISParticle(12); };
inline std::vector GetAllISNuMuon (void) { return GetAllISParticle(14); };
inline std::vector GetAllISNuTau (void) { return GetAllISParticle(16); };
inline std::vector GetAllISElectron (void) { return GetAllISParticle(11); };
inline std::vector GetAllISMuon (void) { return GetAllISParticle(13); };
inline std::vector GetAllISTau (void) { return GetAllISParticle(15); };
inline std::vector GetAllISProton (void) { return GetAllISParticle(2212); };
inline std::vector GetAllISNeutron (void) { return GetAllISParticle(2112); };
inline std::vector GetAllISPiZero (void) { return GetAllISParticle(111); };
inline std::vector GetAllISPiPlus (void) { return GetAllISParticle(211); };
inline std::vector GetAllISPiMinus (void) { return GetAllISParticle(-211); };
inline std::vector GetAllISPhoton (void) { return GetAllISParticle(22); };
inline std::vector GetAllISLeptons (void) { return GetAllISParticle(PhysConst::pdg_leptons); };
inline std::vector GetAllISPions (void) { return GetAllISParticle(PhysConst::pdg_pions); };
inline std::vector GetAllISChargePions(void) { return GetAllISParticle(PhysConst::pdg_charged_pions); };
inline FitParticle* GetHMISNuElectron (void) { return GetHMISParticle(12); };
inline FitParticle* GetHMISNuMuon (void) { return GetHMISParticle(14); };
inline FitParticle* GetHMISNuTau (void) { return GetHMISParticle(16); };
inline FitParticle* GetHMISElectron (void) { return GetHMISParticle(11); };
inline FitParticle* GetHMISMuon (void) { return GetHMISParticle(13); };
inline FitParticle* GetHMISTau (void) { return GetHMISParticle(15); };
inline FitParticle* GetHMISAnyLeptons (void) { return GetHMISParticle(PhysConst::pdg_all_leptons); };
inline FitParticle* GetHMISProton (void) { return GetHMISParticle(2212); };
inline FitParticle* GetHMISNeutron (void) { return GetHMISParticle(2112); };
inline FitParticle* GetHMISPiZero (void) { return GetHMISParticle(111); };
inline FitParticle* GetHMISPiPlus (void) { return GetHMISParticle(211); };
inline FitParticle* GetHMISPiMinus (void) { return GetHMISParticle(-211); };
inline FitParticle* GetHMISPhoton (void) { return GetHMISParticle(22); };
inline FitParticle* GetHMISLepton (void) { return GetHMISParticle(PhysConst::pdg_leptons); };
inline FitParticle* GetHMISPions (void) { return GetHMISParticle(PhysConst::pdg_pions); };
inline FitParticle* GetHMISChargePions(void) { return GetHMISParticle(PhysConst::pdg_charged_pions); };
inline int GetHMISNuElectronIndex (void) { return GetHMISParticleIndex(12); };
inline int GetHMISNuMuonIndex (void) { return GetHMISParticleIndex(14); };
inline int GetHMISNuTauIndex (void) { return GetHMISParticleIndex(16); };
inline int GetHMISElectronIndex (void) { return GetHMISParticleIndex(11); };
inline int GetHMISMuonIndex (void) { return GetHMISParticleIndex(13); };
inline int GetHMISTauIndex (void) { return GetHMISParticleIndex(15); };
inline int GetHMISProtonIndex (void) { return GetHMISParticleIndex(2212); };
inline int GetHMISNeutronIndex (void) { return GetHMISParticleIndex(2112); };
inline int GetHMISPiZeroIndex (void) { return GetHMISParticleIndex(111); };
inline int GetHMISPiPlusIndex (void) { return GetHMISParticleIndex(211); };
inline int GetHMISPiMinusIndex (void) { return GetHMISParticleIndex(-211); };
inline int GetHMISPhotonIndex (void) { return GetHMISParticleIndex(22); };
inline int GetHMISLeptonsIndex (void) { return GetHMISParticleIndex(PhysConst::pdg_leptons); };
inline int GetHMISPionsIndex (void) { return GetHMISParticleIndex(PhysConst::pdg_pions); };
inline int GetHMISChargePionsIndex(void) { return GetHMISParticleIndex(PhysConst::pdg_charged_pions); };
// ---- Final State Helpers --- //
inline bool HasFSParticle(int const pdg) const {
return HasParticle(pdg, kFinalState);
};
template
inline bool HasFSParticle(int const (&pdgs)[N]) const {
return HasParticle(pdgs, kFinalState);
};
inline int NumFSParticle(int const pdg = 0) const {
return NumParticle(pdg, kFinalState);
};
template
inline int NumFSParticle(int const (&pdgs)[N]) const {
return NumParticle(pdgs, kFinalState);
};
inline std::vector GetAllFSParticleIndices(int const pdg = -1) const {
return GetAllParticleIndices(pdg, kFinalState);
};
template
inline std::vector GetAllFSParticleIndices(int const (&pdgs)[N]) const {
return GetAllParticleIndices(pdgs, kFinalState);
};
inline std::vector GetAllFSParticle(int const pdg = -1) {
return GetAllParticle(pdg, kFinalState);
};
template
inline std::vector GetAllFSParticle(int const (&pdgs)[N]) {
return GetAllParticle(pdgs, kFinalState);
};
inline FitParticle* GetHMFSParticle(int const pdg) {
return GetHMParticle(pdg, kFinalState);
};
template
inline FitParticle* GetHMFSParticle(int const (&pdgs)[N]) {
return GetHMParticle(pdgs, kFinalState);
};
inline int GetHMFSParticleIndex(int const pdg) const {
return GetHMParticleIndex(pdg, kFinalState);
};
template
inline int GetHMFSParticleIndex(int const (&pdgs)[N]) const {
return GetHMParticleIndex(pdgs, kFinalState);
};
inline bool HasFSNuElectron (void) const { return HasFSParticle(12); };
inline bool HasFSNuMuon (void) const { return HasFSParticle(14); };
inline bool HasFSNuTau (void) const { return HasFSParticle(16); };
inline bool HasFSElectron (void) const { return HasFSParticle(11); };
inline bool HasFSMuon (void) const { return HasFSParticle(13); };
inline bool HasFSTau (void) const { return HasFSParticle(15); };
inline bool HasFSProton (void) const { return HasFSParticle(2212); };
inline bool HasFSNeutron (void) const { return HasFSParticle(2112); };
inline bool HasFSPiZero (void) const { return HasFSParticle(111); };
inline bool HasFSPiPlus (void) const { return HasFSParticle(211); };
inline bool HasFSPiMinus (void) const { return HasFSParticle(-211); };
inline bool HasFSPhoton (void) const { return HasFSParticle(22); };
inline bool HasFSLeptons (void) const { return HasFSParticle(PhysConst::pdg_leptons); };
inline bool HasFSPions (void) const { return HasFSParticle(PhysConst::pdg_pions); };
inline bool HasFSChargePions (void) const { return HasFSParticle(PhysConst::pdg_charged_pions); };
inline int NumFSNuElectron (void) const { return NumFSParticle(12); };
inline int NumFSNuMuon (void) const { return NumFSParticle(14); };
inline int NumFSNuTau (void) const { return NumFSParticle(16); };
inline int NumFSElectron (void) const { return NumFSParticle(11); };
inline int NumFSMuon (void) const { return NumFSParticle(13); };
inline int NumFSTau (void) const { return NumFSParticle(15); };
inline int NumFSProton (void) const { return NumFSParticle(2212); };
inline int NumFSNeutron (void) const { return NumFSParticle(2112); };
inline int NumFSPiZero (void) const { return NumFSParticle(111); };
inline int NumFSPiPlus (void) const { return NumFSParticle(211); };
inline int NumFSPiMinus (void) const { return NumFSParticle(-211); };
inline int NumFSPhoton (void) const { return NumFSParticle(22); };
int NumFSLeptons (void) const; // { return NumFSParticle(PhysConst::pdg_leptons); };
inline int NumFSPions (void) const { return NumFSParticle(PhysConst::pdg_pions); };
inline int NumFSChargePions (void) const { return NumFSParticle(PhysConst::pdg_charged_pions); };
inline std::vector GetAllFSNuElectronIndices (void) const { return GetAllFSParticleIndices(12); };
inline std::vector GetAllFSNuMuonIndices (void) const { return GetAllFSParticleIndices(14); };
inline std::vector GetAllFSNuTauIndices (void) const { return GetAllFSParticleIndices(16); };
inline std::vector GetAllFSElectronIndices (void) const { return GetAllFSParticleIndices(11); };
inline std::vector GetAllFSMuonIndices (void) const { return GetAllFSParticleIndices(13); };
inline std::vector GetAllFSTauIndices (void) const { return GetAllFSParticleIndices(15); };
inline std::vector GetAllFSProtonIndices (void) const { return GetAllFSParticleIndices(2212); };
inline std::vector GetAllFSNeutronIndices (void) const { return GetAllFSParticleIndices(2112); };
inline std::vector GetAllFSPiZeroIndices (void) const { return GetAllFSParticleIndices(111); };
inline std::vector GetAllFSPiPlusIndices (void) const { return GetAllFSParticleIndices(211); };
inline std::vector GetAllFSPiMinusIndices (void) const { return GetAllFSParticleIndices(-211); };
inline std::vector GetAllFSPhotonIndices (void) const { return GetAllFSParticleIndices(22); };
inline std::vector GetAllFSLeptonsIndices (void) const { return GetAllFSParticleIndices(PhysConst::pdg_leptons); };
inline std::vector GetAllFSPionsIndices (void) const { return GetAllFSParticleIndices(PhysConst::pdg_pions); };
inline std::vector GetAllFSChargePionsIndices(void) const { return GetAllFSParticleIndices(PhysConst::pdg_charged_pions); };
inline std::vector GetAllFSNuElectron (void) { return GetAllFSParticle(12); };
inline std::vector GetAllFSNuMuon (void) { return GetAllFSParticle(14); };
inline std::vector GetAllFSNuTau (void) { return GetAllFSParticle(16); };
inline std::vector GetAllFSElectron (void) { return GetAllFSParticle(11); };
inline std::vector GetAllFSMuon (void) { return GetAllFSParticle(13); };
inline std::vector GetAllFSTau (void) { return GetAllFSParticle(15); };
inline std::vector GetAllFSProton (void) { return GetAllFSParticle(2212); };
inline std::vector GetAllFSNeutron (void) { return GetAllFSParticle(2112); };
inline std::vector GetAllFSPiZero (void) { return GetAllFSParticle(111); };
inline std::vector GetAllFSPiPlus (void) { return GetAllFSParticle(211); };
inline std::vector GetAllFSPiMinus (void) { return GetAllFSParticle(-211); };
inline std::vector GetAllFSPhoton (void) { return GetAllFSParticle(22); };
inline std::vector GetAllFSLeptons (void) { return GetAllFSParticle(PhysConst::pdg_leptons); };
inline std::vector GetAllFSPions (void) { return GetAllFSParticle(PhysConst::pdg_pions); };
inline std::vector GetAllFSChargePions (void) { return GetAllFSParticle(PhysConst::pdg_charged_pions); };
inline FitParticle* GetHMFSNuElectron (void) { return GetHMFSParticle(12); };
inline FitParticle* GetHMFSNuMuon (void) { return GetHMFSParticle(14); };
inline FitParticle* GetHMFSNuTau (void) { return GetHMFSParticle(16); };
inline FitParticle* GetHMFSElectron (void) { return GetHMFSParticle(11); };
inline FitParticle* GetHMFSMuon (void) { return GetHMFSParticle(13); };
inline FitParticle* GetHMFSTau (void) { return GetHMFSParticle(15); };
inline FitParticle* GetHMFSAnyLeptons (void) { return GetHMFSParticle(PhysConst::pdg_all_leptons); };
inline FitParticle* GetHMFSProton (void) { return GetHMFSParticle(2212); };
inline FitParticle* GetHMFSNeutron (void) { return GetHMFSParticle(2112); };
inline FitParticle* GetHMFSPiZero (void) { return GetHMFSParticle(111); };
inline FitParticle* GetHMFSPiPlus (void) { return GetHMFSParticle(211); };
inline FitParticle* GetHMFSPiMinus (void) { return GetHMFSParticle(-211); };
inline FitParticle* GetHMFSPhoton (void) { return GetHMFSParticle(22); };
inline FitParticle* GetHMFSLeptons (void) { return GetHMFSParticle(PhysConst::pdg_leptons); };
inline FitParticle* GetHMFSAnyLepton (void) { return GetHMFSParticle(PhysConst::pdg_all_leptons); };
inline FitParticle* GetHMFSPions (void) { return GetHMFSParticle(PhysConst::pdg_pions); };
inline FitParticle* GetHMFSChargePions(void) { return GetHMFSParticle(PhysConst::pdg_charged_pions); };
inline int GetHMFSNuElectronIndex (void) const { return GetHMFSParticleIndex(12); };
inline int GetHMFSNuMuonIndex (void) const { return GetHMFSParticleIndex(14); };
inline int GetHMFSNuTauIndex (void) const { return GetHMFSParticleIndex(16); };
inline int GetHMFSElectronIndex (void) const { return GetHMFSParticleIndex(11); };
inline int GetHMFSMuonIndex (void) const { return GetHMFSParticleIndex(13); };
inline int GetHMFSTauIndex (void) const { return GetHMFSParticleIndex(15); };
inline int GetHMFSProtonIndex (void) const { return GetHMFSParticleIndex(2212); };
inline int GetHMFSNeutronIndex (void) const { return GetHMFSParticleIndex(2112); };
inline int GetHMFSPiZeroIndex (void) const { return GetHMFSParticleIndex(111); };
inline int GetHMFSPiPlusIndex (void) const { return GetHMFSParticleIndex(211); };
inline int GetHMFSPiMinusIndex (void) const { return GetHMFSParticleIndex(-211); };
inline int GetHMFSPhotonIndex (void) const { return GetHMFSParticleIndex(22); };
inline int GetHMFSLeptonsIndex (void) const { return GetHMFSParticleIndex(PhysConst::pdg_leptons); };
inline int GetHMFSAnyLeptonIndex (void) const { return GetHMFSParticleIndex(PhysConst::pdg_all_leptons); };
inline int GetHMFSPionsIndex (void) const { return GetHMFSParticleIndex(PhysConst::pdg_pions); };
inline int GetHMFSChargePionsIndex(void) const { return GetHMFSParticleIndex(PhysConst::pdg_charged_pions); };
// ---- NEUTRINO INCOMING Related Functions
int GetBeamNeutrinoIndex (void) const;
inline TLorentzVector GetBeamNeutrinoP4 (void) const { return GetParticleP4(GetBeamNeutrinoIndex()); };
inline TVector3 GetBeamNeutrinoP3 (void) const { return GetParticleP3(GetBeamNeutrinoIndex()); };
inline double GetBeamNeutrinoMom (void) const { return GetParticleMom(GetBeamNeutrinoIndex()); };
inline double GetBeamNeutrinoMom2 (void) const { return GetParticleMom2(GetBeamNeutrinoIndex()); };
inline double GetBeamNeutrinoE (void) const { return GetParticleE(GetBeamNeutrinoIndex()); };
inline double Enu (void) const { return GetBeamNeutrinoE(); };
inline int GetBeamNeutrinoPDG (void) const { return GetParticlePDG(GetBeamNeutrinoIndex()); };
inline int PDGnu (void) const { return GetBeamNeutrinoPDG(); };
inline int GetNeutrinoInPos (void) const { return GetBeamNeutrinoIndex(); };
inline FitParticle* GetBeamNeutrino (void) { return GetParticle(GetBeamNeutrinoIndex()); };
inline FitParticle* GetNeutrinoIn (void) { return GetParticle(GetBeamNeutrinoIndex()); };
// ---- Electron INCOMING Related Functions
int GetBeamElectronIndex (void) const;
inline TLorentzVector GetBeamElectronP4 (void) const { return GetParticleP4(GetBeamElectronIndex()); };
inline TVector3 GetBeamElectronP3 (void) const { return GetParticleP3(GetBeamElectronIndex()); };
inline double GetBeamElectronMom (void) const { return GetParticleMom(GetBeamElectronIndex()); };
inline double GetBeamElectronMom2 (void) const { return GetParticleMom2(GetBeamElectronIndex()); };
inline double GetBeamElectronE (void) const { return GetParticleE(GetBeamElectronIndex()); };
inline FitParticle* GetBeamElectron (void) { return GetParticle(GetBeamElectronIndex()); };
// ---- Pion INCOMING Functions
int GetBeamPionIndex (void) const;
inline TLorentzVector GetBeamPionP4 (void) const { return GetParticleP4(GetBeamPionIndex()); };
inline TVector3 GetBeamPionP3 (void) const { return GetParticleP3(GetBeamPionIndex()); };
inline double GetBeamPionMom (void) const { return GetParticleMom(GetBeamPionIndex()); };
inline double GetBeamPionMom2 (void) const { return GetParticleMom2(GetBeamPionIndex()); };
inline double GetBeamPionE (void) const { return GetParticleE(GetBeamPionIndex()); };
inline FitParticle* GetBeamPion (void) { return GetParticle(GetBeamPionIndex()); };
/// Legacy Functions
inline FitParticle* PartInfo(uint i) { return GetParticle(i); };
inline UInt_t Npart (void) const { return NPart(); };
inline UInt_t NPart (void) const { return fNParticles; };
// Other Functions
int NumFSMesons();
inline int GetBeamPartPos(void) const { return 0; };
FitParticle* GetBeamPart(void);
int GetLeptonOutPos(void) const;
FitParticle* GetLeptonOut(void);
// Event Information
- double weight; // need for silly reason
- int Mode; // Public access needed
-
- int fMode;
UInt_t fEventNo;
double fTotCrs;
int fTargetA;
int fTargetZ;
int fTargetH;
bool fBound;
int fDistance;
int fTargetPDG;
// Reduced Particle Stack
UInt_t kMaxParticles;
int fNParticles;
double** fParticleMom;
UInt_t* fParticleState;
int* fParticlePDG;
FitParticle** fParticleList;
double** fOrigParticleMom;
UInt_t* fOrigParticleState;
int* fOrigParticlePDG;
double* fNEUT_ParticleStatusCode;
double* fNEUT_ParticleAliveCode;
GeneratorInfoBase* fGenInfo;
// Config Options for this class
bool kRemoveFSIParticles;
bool kRemoveUndefParticles;
};
/*! @} */
#endif
diff --git a/src/InputHandler/GENIEInputHandler.cxx b/src/InputHandler/GENIEInputHandler.cxx
index fa85737..7ee3454 100644
--- a/src/InputHandler/GENIEInputHandler.cxx
+++ b/src/InputHandler/GENIEInputHandler.cxx
@@ -1,451 +1,450 @@
// 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 .
*******************************************************************************/
#ifdef __GENIE_ENABLED__
#include "GENIEInputHandler.h"
#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(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((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;
// 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");
// 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(
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]);
}
// 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());
}
// Registor all our file inputs
SetupJointInputs();
// Assign to tree
fEventType = kGENIE;
fGenieNtpl = NULL;
fGENIETree->SetBranchAddress("gmcrec", &fGenieNtpl);
fGENIETree->GetEntry(0);
// 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;
// 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(fGenieNtpl->event);
if (!fGenieGHep) return NULL;
TObjArrayIter iter(fGenieGHep);
genie::GHepParticle* p;
while ((p = (dynamic_cast((iter).Next())))) {
if (!p) {
continue;
}
// Get Status
- int state = GetGENIEParticleStatus(p, fNUISANCEEvent->fMode);
+ 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(fGenieNtpl->event);
if (!fGenieGHep) return;
// Convert GENIE Reaction Code
- fNUISANCEEvent->fMode = ConvertGENIEReactionCode(fGenieGHep);
+ fNUISANCEEvent->Mode = ConvertGENIEReactionCode(fGenieGHep);
// Set Event Info
- fNUISANCEEvent->Mode = fNUISANCEEvent->fMode;
fNUISANCEEvent->fEventNo = 0.0;
fNUISANCEEvent->fTotCrs = fGenieGHep->XSec();
fNUISANCEEvent->fTargetA = 0.0;
fNUISANCEEvent->fTargetZ = 0.0;
fNUISANCEEvent->fTargetH = 0;
fNUISANCEEvent->fBound = 0.0;
fNUISANCEEvent->InputWeight =
1.0; //(1E+38 / genie::units::cm2) * fGenieGHep->XSec();
// 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((iter).Next())))) {
if (!p) continue;
// Get Status
- int state = GetGENIEParticleStatus(p, fNUISANCEEvent->fMode);
+ 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();
// 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/InputHandler/GIBUUInputHandler.cxx b/src/InputHandler/GIBUUInputHandler.cxx
index 880a563..a8c961b 100644
--- a/src/InputHandler/GIBUUInputHandler.cxx
+++ b/src/InputHandler/GIBUUInputHandler.cxx
@@ -1,301 +1,300 @@
#ifdef __GiBUU_ENABLED__
#include "GIBUUInputHandler.h"
#include "InputUtils.h"
GIBUUGeneratorInfo::~GIBUUGeneratorInfo() { DeallocateParticleStack(); }
void GIBUUGeneratorInfo::AddBranchesToTree(TTree* tn) {
// tn->Branch("NEUTParticleN", fNEUTParticleN, "NEUTParticleN/I");
// tn->Branch("NEUTParticleStatusCode", fNEUTParticleStatusCode,
// "NEUTParticleStatusCode[NEUTParticleN]/I");
// tn->Branch("NEUTParticleAliveCode", fNEUTParticleAliveCode,
// "NEUTParticleAliveCode[NEUTParticleN]/I");
}
void GIBUUGeneratorInfo::SetBranchesFromTree(TTree* tn) {
// tn->SetBranchAddress("NEUTParticleN", &fNEUTParticleN );
// tn->SetBranchAddress("NEUTParticleStatusCode", &fNEUTParticleStatusCode );
// tn->SetBranchAddress("NEUTParticleAliveCode", &fNEUTParticleAliveCode );
}
void GIBUUGeneratorInfo::AllocateParticleStack(int stacksize) {
// fNEUTParticleN = 0;
// fNEUTParticleStatusCode = new int[stacksize];
// fNEUTParticleStatusCode = new int[stacksize];
}
void GIBUUGeneratorInfo::DeallocateParticleStack() {
// delete fNEUTParticleStatusCode;
// delete fNEUTParticleAliveCode;
}
void GIBUUGeneratorInfo::FillGeneratorInfo(GiBUUStdHepReader* nevent) {
Reset();
// for (int i = 0; i < nevent->Npart(); i++) {
// fNEUTParticleStatusCode[i] = nevent->PartInfo(i)->fStatus;
// fNEUTParticleAliveCode[i] = nevent->PartInfo(i)->fIsAlive;
// fNEUTParticleN++;
// }
}
void GIBUUGeneratorInfo::Reset() {
// for (int i = 0; i < fNEUTParticleN; i++) {
// fNEUTParticleStatusCode[i] = -1;
// fNEUTParticleAliveCode[i] = 9;
// }
// fNEUTParticleN = 0;
}
GIBUUInputHandler::GIBUUInputHandler(std::string const& handle,
std::string const& rawinputs) {
LOG(SAM) << "Creating GiBUUInputHandler : " << handle << std::endl;
// Run a joint input handling
fName = handle;
fEventType = kGiBUU;
fGIBUUTree = new TChain("giRooTracker");
// 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
LOG(SAM) << "Opening event file " << inputs[inp_it] << std::endl;
TFile* inp_file = new TFile(inputs[inp_it].c_str(), "READ");
if ((!inp_file) || (!inp_file->IsOpen())) {
THROW("GiBUU file !IsOpen() at : '"
<< inputs[inp_it] << "'" << std::endl
<< "Check that your file paths are correct and the file exists!");
}
int NFluxes = bool(dynamic_cast(inp_file->Get("numu_flux"))) +
bool(dynamic_cast(inp_file->Get("numub_flux"))) +
bool(dynamic_cast(inp_file->Get("nue_flux"))) +
bool(dynamic_cast(inp_file->Get("nueb_flux")));
if (NFluxes != 1) {
THROW("Found " << NFluxes << " input fluxes in " << inputs[inp_it]
<< ". The NUISANCE GiBUU interface expects to be "
"passed multiple species vectors as separate "
"input files like: "
"\"GiBUU:(MINERVA_FHC_numu_evts.root,MINERVA_FHC_"
"numubar_evts.root,[...])\"");
}
// Get Flux/Event hist
TH1D* fluxhist = dynamic_cast(inp_file->Get("flux"));
TH1D* eventhist = dynamic_cast(inp_file->Get("evt"));
if (!fluxhist || !eventhist) {
ERROR(FTL, "Input File Contents: " << inputs[inp_it]);
inp_file->ls();
THROW(
"GiBUU FILE doesn't contain flux/xsec info. You may have to "
"regenerate your MC!");
}
// Get N Events
TTree* giRooTracker = dynamic_cast(inp_file->Get("giRooTracker"));
if (!giRooTracker) {
ERROR(FTL,
"giRooTracker Tree not located in NEUT file: " << inputs[inp_it]);
THROW("Check your inputs, they may need to be completely regenerated!");
throw;
}
int nevents = giRooTracker->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
fGIBUUTree->AddFile(inputs[inp_it].c_str());
}
// Registor all our file inputs
SetupJointInputs();
// Create Fit Event
fNUISANCEEvent = new FitEvent();
fGiReader = new GiBUUStdHepReader();
fGiReader->SetBranchAddresses(fGIBUUTree);
fNUISANCEEvent->HardReset();
};
FitEvent* GIBUUInputHandler::GetNuisanceEvent(const UInt_t entry,
const bool lightweight) {
// Check out of bounds
if (entry >= (UInt_t)fNEvents) return NULL;
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fGIBUUTree->GetEntry(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
#ifdef __PROB3PP_ENABLED__
else {
for (int i = 0; i < fGiReader->StdHepN; i++) {
int state = GetGIBUUParticleStatus(fGiReader->StdHepStatus[i],
fGiReader->StdHepPdg[i]);
if (state != kInitialState) {
continue;
}
if (std::count(PhysConst::pdg_neutrinos, PhysConst::pdg_neutrinos + 4,
fGiReader->StdHepPdg[i])) {
fNUISANCEEvent->probe_E = fGiReader->StdHepP4[i][3] * 1.E3;
fNUISANCEEvent->probe_pdg = fGiReader->StdHepPdg[i];
break;
}
}
}
#endif
fNUISANCEEvent->InputWeight *= GetInputWeight(entry);
return fNUISANCEEvent;
}
int GetGIBUUParticleStatus(int status, int pdg) {
int state = kUndefinedState;
switch (status) {
case 0: // Incoming
case 11: // Struck nucleon
state = kInitialState;
break;
case 1: // Good Final State
state = kFinalState;
break;
default: // Other
break;
}
// Set Nuclear States Flag
if (pdg > 1000000) {
if (state == kInitialState)
state = kNuclearInitial;
else if (state == kFinalState)
state = kNuclearRemnant;
else
state = kUndefinedState;
}
return state;
}
void GIBUUInputHandler::CalcNUISANCEKinematics() {
// Reset all variables
fNUISANCEEvent->ResetEvent();
FitEvent* evt = fNUISANCEEvent;
- evt->fMode = fGiReader->GiBUU2NeutCode;
- evt->Mode = evt->fMode;
+ evt->Mode = fGiReader->GiBUU2NeutCode;
evt->fEventNo = 0.0;
evt->fTotCrs = 0;
evt->fTargetA = 0.0; // Change to get these from nuclear remnant.
evt->fTargetZ = 0.0;
evt->fTargetH = 0;
evt->fBound = 0.0;
// Extra GiBUU Input Weight
evt->InputWeight = fGiReader->EvtWght;
// Check Stack N
int npart = fGiReader->StdHepN;
int kmax = evt->kMaxParticles;
if ((UInt_t)npart > (UInt_t)kmax) {
ERROR(WRN, "GiBUU has too many particles. Expanding Stack.");
fNUISANCEEvent->ExpandParticleStack(npart);
}
// Create Stack
evt->fNParticles = 0;
for (int i = 0; i < npart; i++) {
// State
int state = GetGIBUUParticleStatus(fGiReader->StdHepStatus[i],
fGiReader->StdHepPdg[i]);
int curpart = evt->fNParticles;
// Set State
evt->fParticleState[evt->fNParticles] = state;
// Mom
evt->fParticleMom[curpart][0] = fGiReader->StdHepP4[i][0] * 1.E3;
evt->fParticleMom[curpart][1] = fGiReader->StdHepP4[i][1] * 1.E3;
evt->fParticleMom[curpart][2] = fGiReader->StdHepP4[i][2] * 1.E3;
evt->fParticleMom[curpart][3] = fGiReader->StdHepP4[i][3] * 1.E3;
// PDG
evt->fParticlePDG[curpart] = fGiReader->StdHepPdg[i];
// Add to total particles
evt->fNParticles++;
}
// 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 GIBUUInputHandler::Print() {}
void GIBUUInputHandler::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) {
THROW("Can only handle joint inputs when config MAXEVENTS = -1!");
}
for (size_t i = 0; i < jointeventinputs.size(); i++) {
double scale = double(fNEvents) / fEventHist->Integral("width");
scale *= jointfluxinputs.at(i)->Integral("width");
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;
}
}
#endif
diff --git a/src/InputHandler/NEUTInputHandler.cxx b/src/InputHandler/NEUTInputHandler.cxx
index d00defd..189026d 100644
--- a/src/InputHandler/NEUTInputHandler.cxx
+++ b/src/InputHandler/NEUTInputHandler.cxx
@@ -1,460 +1,459 @@
#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) {
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 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()) {
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) {
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) {
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);
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 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);
// 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;
}
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();
FitParticle* ISNeutralLepton =
fNUISANCEEvent->GetHMISParticle(PhysConst::pdg_neutrinos);
if (ISNeutralLepton) {
fNUISANCEEvent->probe_E = ISNeutralLepton->E();
fNUISANCEEvent->probe_pdg = ISNeutralLepton->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();
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/InputHandler/NUANCEInputHandler.cxx b/src/InputHandler/NUANCEInputHandler.cxx
index ad98dcd..7814f4a 100644
--- a/src/InputHandler/NUANCEInputHandler.cxx
+++ b/src/InputHandler/NUANCEInputHandler.cxx
@@ -1,837 +1,836 @@
#ifdef __NUANCE_ENABLED__
#include "NUANCEInputHandler.h"
#include "InputUtils.h"
NUANCEGeneratorInfo::~NUANCEGeneratorInfo() {
DeallocateParticleStack();
}
void NUANCEGeneratorInfo::AddBranchesToTree(TTree * tn) {
// tn->Branch("NEUTParticleN", fNEUTParticleN, "NEUTParticleN/I");
// tn->Branch("NEUTParticleStatusCode", fNEUTParticleStatusCode, "NEUTParticleStatusCode[NEUTParticleN]/I");
// tn->Branch("NEUTParticleAliveCode", fNEUTParticleAliveCode, "NEUTParticleAliveCode[NEUTParticleN]/I");
}
void NUANCEGeneratorInfo::SetBranchesFromTree(TTree* tn) {
// tn->SetBranchAddress("NEUTParticleN", &fNEUTParticleN );
// tn->SetBranchAddress("NEUTParticleStatusCode", &fNEUTParticleStatusCode );
// tn->SetBranchAddress("NEUTParticleAliveCode", &fNEUTParticleAliveCode );
}
void NUANCEGeneratorInfo::AllocateParticleStack(int stacksize) {
// fNEUTParticleN = 0;
// fNEUTParticleStatusCode = new int[stacksize];
// fNEUTParticleStatusCode = new int[stacksize];
}
void NUANCEGeneratorInfo::DeallocateParticleStack() {
// delete fNEUTParticleStatusCode;
// delete fNEUTParticleAliveCode;
}
void NUANCEGeneratorInfo::FillGeneratorInfo(NuanceEvent* nevent) {
Reset();
// for (int i = 0; i < nevent->Npart(); i++) {
// fNEUTParticleStatusCode[i] = nevent->PartInfo(i)->fStatus;
// fNEUTParticleAliveCode[i] = nevent->PartInfo(i)->fIsAlive;
// fNEUTParticleN++;
// }
}
void NUANCEGeneratorInfo::Reset() {
// for (int i = 0; i < fNEUTParticleN; i++) {
// fNEUTParticleStatusCode[i] = -1;
// fNEUTParticleAliveCode[i] = 9;
// }
// fNEUTParticleN = 0;
}
NUANCEInputHandler::NUANCEInputHandler(std::string const& handle, std::string const& rawinputs) {
LOG(SAM) << "Creating NUANCEInputHandler : " << handle << std::endl;
// Run a joint input handling
fName = handle;
fSaveExtra = FitPar::Config().GetParB("SaveExtraNUANCE");
fCacheSize = FitPar::Config().GetParI("CacheSize");
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
// Parse Inputs
std::vector inputs = InputUtils::ParseInputFileList(rawinputs);
if (inputs.size() > 1) {
ERR(FTL) << "NUANCE is not currently setup to handle joint inputs sorry!" << std::endl
<< "If you know how to correctly normalise the events for this"
<< " please let us know!" << std::endl;
}
// Read in NUANCE Tree
fNUANCETree = new TChain("h3");
fNUANCETree->AddFile(rawinputs.c_str());
// Get entries and fNuwroEvent
int nevents = fNUANCETree->GetEntries();
double EnuMin = 0.0;
double EnuMax = 1000.0;
TH1D* fluxhist = new TH1D((fName + "_FLUX").c_str(), (fName + "_FLUX").c_str(),
100, EnuMin, EnuMax);
for (int i = 0; i < fluxhist->GetNbinsX(); i++) {
fluxhist->SetBinContent(i + 1, 1.0);
}
fluxhist->Scale(1.0 / fluxhist->Integral());
TH1D* eventhist = new TH1D((fName + "_EVT").c_str(), (fName + "_EVT").c_str(), 100,
EnuMin, EnuMax);
for (int i = 0; i < fluxhist->GetNbinsX(); i++) {
eventhist->SetBinContent(i + 1, 1.0);
}
eventhist->Scale(1.0 / eventhist->Integral());
RegisterJointInput( rawinputs, nevents, fluxhist, eventhist );
SetupJointInputs();
// Setup Reader
fNuanceEvent = new NuanceEvent();
fNuanceEvent->SetBranchAddresses(fNUANCETree);
fNUANCETree->GetEntry(0);
// Setup Event in FitEvent
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->SetNuanceEvent(fNuanceEvent);
// Setup extra if needed
if (fSaveExtra){
ERR(FTL) << "NO SAVEExtra Implemented for NUANCE YET!" << std::endl;
throw;
// fNuanceInfo = new NUANCEGeneratorInfo();
// fNUISANCEEvent->AddGeneratorInfo(fNuanceInfo);
}
};
NUANCEInputHandler::~NUANCEInputHandler() {
if (fNuanceEvent) delete fNuanceEvent;
if (fNUANCETree) delete fNUANCETree;
// if (fNuanceInfo) delete fNuanceInfo;
}
void NUANCEInputHandler::CreateCache() {
if (fCacheSize > 0) {
fNUANCETree->SetCacheEntryRange(0, fNEvents);
fNUANCETree->AddBranchToCache("h3", 1);
fNUANCETree->SetCacheSize(fCacheSize);
}
}
void NUANCEInputHandler::RemoveCache() {
fNUANCETree->SetCacheEntryRange(0, fNEvents);
fNUANCETree->AddBranchToCache("h3", 0);
fNUANCETree->SetCacheSize(0);
}
FitEvent* NUANCEInputHandler::GetNuisanceEvent(const UInt_t entry, const bool lightweight) {
// Check out of bounds
if (entry >= (UInt_t)fNEvents) return NULL;
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fNUANCETree->GetEntry(entry);
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
return fNUISANCEEvent;
}
void NUANCEInputHandler::CalcNUISANCEKinematics() {
// Reset all variables
fNUISANCEEvent->ResetEvent();
// Get shortened pointer
FitEvent* evt = fNUISANCEEvent;
// Fill Global
- evt->fMode = ConvertNuanceMode(fNuanceEvent);
- evt->Mode = evt->fMode;
+ evt->Mode = ConvertNuanceMode(fNuanceEvent);
evt->fEventNo = 0.0;
evt->fTotCrs = 1.0;
evt->fTargetA = 0.0;
evt->fTargetZ = 0.0;
evt->fTargetH = 0;
evt->fBound = 0.0;
// Fill particle Stack
evt->fNParticles = 0;
// Check Particle Stack
UInt_t npart = 2 + fNuanceEvent->n_leptons + fNuanceEvent->n_hadrons;
UInt_t kmax = evt->kMaxParticles;
if (npart > kmax) {
ERR(FTL) << "NUANCE has too many particles" << std::endl;
ERR(FTL) << "npart=" << npart << " kMax=" << kmax << std::endl;
throw;
}
// Fill Neutrino
evt->fParticleState[0] = kInitialState;
evt->fParticleMom[0][0] = fNuanceEvent->p_neutrino[0];
evt->fParticleMom[0][1] = fNuanceEvent->p_neutrino[1];
evt->fParticleMom[0][2] = fNuanceEvent->p_neutrino[2];
evt->fParticleMom[0][3] = fNuanceEvent->p_neutrino[3];
evt->fParticlePDG[0] = fNuanceEvent->neutrino;
// Fill Target Nucleon
evt->fParticleState[1] = kInitialState;
evt->fParticleMom[1][0] = fNuanceEvent->p_targ[0];
evt->fParticleMom[1][1] = fNuanceEvent->p_targ[1];
evt->fParticleMom[1][2] = fNuanceEvent->p_targ[2];
evt->fParticleMom[1][3] = fNuanceEvent->p_targ[3];
evt->fParticlePDG[1] = fNuanceEvent->target;
evt->fNParticles = 2;
// Fill Outgoing Leptons
for (int i = 0; i < fNuanceEvent->n_leptons; i++) {
evt->fParticleState[evt->fNParticles] = kFinalState;
evt->fParticleMom[evt->fNParticles][0] = fNuanceEvent->p_lepton[i][0];
evt->fParticleMom[evt->fNParticles][1] = fNuanceEvent->p_lepton[i][1];
evt->fParticleMom[evt->fNParticles][2] = fNuanceEvent->p_lepton[i][2];
evt->fParticleMom[evt->fNParticles][3] = fNuanceEvent->p_lepton[i][3];
evt->fParticlePDG[evt->fNParticles] = fNuanceEvent->lepton[i];
evt->fNParticles++;
}
// Fill Outgoing Hadrons
for (int i = 0; i < fNuanceEvent->n_hadrons; i++) {
evt->fParticleState[evt->fNParticles] = kFinalState;
evt->fParticleMom[evt->fNParticles][0] = fNuanceEvent->p_hadron[i][0];
evt->fParticleMom[evt->fNParticles][1] = fNuanceEvent->p_hadron[i][1];
evt->fParticleMom[evt->fNParticles][2] = fNuanceEvent->p_hadron[i][2];
evt->fParticleMom[evt->fNParticles][3] = fNuanceEvent->p_hadron[i][3];
evt->fParticlePDG[evt->fNParticles] = fNuanceEvent->hadron[i];
evt->fNParticles++;
}
// Save Extra info
if (fSaveExtra) {
// fNuanceInfo->FillGeneratorInfo(fNuanceEvent);
}
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent-> OrderStack();
return;
}
void NUANCEInputHandler::Print() {}
int NUANCEInputHandler::ConvertNuanceMode(NuanceEvent * evt) {
int ch = evt->channel;
int sg = 1;
if (evt->neutrino < 0) sg = -1;
switch (ch) {
// 1 NUANCE CCQE -> NEUT CCQE 1
case 1: return sg * 1;
// 2 NUANCE NCEL -> NEUT NCEL 51,52 -> Set from whether target is p or n
case 2:
if (evt->target == 2212) return sg * 51;
else return sg * 52;
// 3 NUANCE CCPIP -> NEUT CCPIP 11
case 3: return sg * 11;
// 4 NUANCE CCPI0 -> NEUT CCPI0 = 12
case 4: return sg * 12;
// 5 NUANCE CCPIPn -> NEUT CCPIPn 13
case 5: return sg * 13;
// 6 NUANCE NCpPI0 -> NEUT NCpPI0 32
case 6: return sg * 32;
// 7 NUANCE NCpPI+ -> NEUT NCpPI+ 34
case 7: return sg * 34;
// 8 NUANCE NCnPI0 -> NEUT NCnPI0 31
case 8: return sg * 31;
// 9 NUANCE NCnPIM -> NEUT NCnPIM 33
case 9: return sg * 33;
// 10 NUANCE CCPIP -> NEUT CCPIP -11
case 10: return sg * 11;
// 11 NUANCE CCPI0 -> NEUT CCPI0 -12
case 11: return sg * 12;
// 12 NUANCE CCPIPn -> NEUT CCPIPn 13
case 12: return sg * 13;
// 13 NUANCE NCpPI0 -> NEUT NCnPI0 -32
case 13: return sg * 32;
// 14 NUANCE NCpPI+ -> NEUT NCpPI+ -34
case 14: return sg * 34;
// 15 NUANCE NCnPI0 -> NEUT NCnPI0 -31
case 15: return sg * 31;
// 16 NUANCE NCnPIM -> NEUT NCnPIM -33
case 16: return sg * 33;
// 17 NUANCE -> NEUT 21 CC MULTIPI
case 17: return sg * 21;
// 18 NUANCE -> NEUT 21 CC MULTIPI
case 18: return sg * 21;
// 19 NUANCE -> NEUT 21 CC MULTIPI
case 19: return sg * 21;
// 20 NUANCE -> NEUT 21 CC MULTIPI
case 20: return sg * 21;
// 21 NUANCE -> NEUT 21 CC MULTIPI
case 21: return sg * 21;
// 22 NUANCE -> NEUT 41 NC MULTIPI
case 22: return sg * 41;
// 23 NUANCE -> NEUT 41 NC MULTIPI
case 23: return sg * 41;
// 24 NUANCE -> NEUT 41 NC MULTIPI
case 24: return sg * 41;
// 25 NUANCE -> NEUT 41 NC MULTIPI
case 25: return sg * 41;
// 26 NUANCE -> NEUT 41 NC MULTIPI
case 26: return sg * 41;
// 27 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV)
case 27: return sg * 41;
// 28 NUANCE -> NEUT -21 CC (1.3 < W < 2 GeV)
case 28: return sg * 21;
// 29 NUANCE -> NEUT -21 CC (1.3 < W < 2 GeV)
case 29: return sg * 21;
// 30 NUANCE -> NEUT -21 CC (1.3 < W < 2 GeV)
case 30: return sg * 21;
// 31 NUANCE -> NEUT -21 CC (1.3 < W < 2 GeV)
case 31: return sg * 21;
// 32 NUANCE -> NEUT -21 CC (1.3 < W < 2 GeV)
case 32: return sg * 21;
// 33 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV)
case 33: return sg * 41;
// 34 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV)
case 34: return sg * 41;
// 35 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV)
case 35: return sg * 41;
// 36 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV)
case 36: return sg * 41;
// 37 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV)
case 37: return sg * 41;
// 38 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV)
case 38: return sg * 41;
// 39 NUANCE -> NEUT 22
case 39: return sg * 22;
// 40 NUANCE -> NEUT 22
case 40: return sg * 22;
// 41 NUANCE -> NEUT 22
case 41: return sg * 22;
// 42 NUANCE -> NEUT 43
case 42: return sg * 43;
// 43 NUANCE -> NEUT 43
case 43: return sg * 43;
// 44 NUANCE -> NUET 42
case 44: return sg * 42;
// 45 NUANCE -> NEUT -42
case 45: return sg * 42;
// 46 NUANCE -> NEUT -22
case 46: return sg * 22;
// 47 NUANCE -> NEUT -22
case 47: return sg * 22;
// 48 NUANCE -> NEUT -22
case 48: return sg * 22;
// 49 NUANCE -> NEUT -43
case 49: return sg * 43;
// 50 NUANCE -> NEUT -43
case 50: return sg * 43;
// 51 NUANCE -> NEUT -42
case 51: return sg * 42;
// 52 NUANCE -> NEUT -42
case 52: return sg * 42;
// 53 NUANCE -> NEUT 23 CC 1K
case 53: return sg * 23;
// 54 NUANCE -> NEUT 23 CC 1K
case 54: return sg * 23;
// 55 NUANCE -> NEUT 23 CC 1K
case 55: return sg * 23;
// 56 NUANCE -> NEUT 45 NC 1K
case 56: return sg * 45;
// 57 NUANCE -> NEUT 44 NC 1K
case 57: return sg * 44;
// 58 NUANCE -> NEUT 44 NC 1K
case 58: return sg * 44;
// 59 NUANCE -> NEUT 44 NC 1K
case 59: return sg * 44;
// 60 NUANCE -> NEUT -23 CC 1K
case 60: return sg * 23;
// 61 NUANCE -> NEUT -23 CC 1K
case 61: return sg * 23;
// 62 NUANCE -> NEUT -23 CC 1K
case 62: return sg * 23;
// 63 NUANCE -> NEUT -23 CC 1K
case 63: return sg * 23;
// 64 NUANCE -> NEUT -44 NC 1K
case 64: return sg * 44;
// 65 NUANCE -> NEUT -44 NC 1K
case 65: return sg * 44;
// 66 NUANCE -> NEUT -45 NC 1K
case 66: return sg * 45;
// 67 NUANCE -> NEUT 22 CC1eta
case 67: return sg * 22;
// 68 NUANCE -> NEUT 43 NC p eta
case 68: return sg * 43;
// 69 NUANCE -> NEUT 43 NC n eta
case 69: return sg * 43;
// 70 NUANCE -> NEUT -22 CC1eta
case 70: return sg * 22;
// 71 NUANCE -> NEUT -43 NC p eta
case 71: return sg * 43;
// 72 NUANCE -> NEUT 42 NC n eta
case 72: return sg * 42;
// 73 NUANCE -> NEUT 21 CC Multi Pi
case 73: return sg * 21;
// 74 NUANCE -> NEUT 41 NC Multi Pi
case 74: return sg * 41;
// 75 NUANCE -> NEUT 41 NC Multi Pi
case 75: return sg * 41;
// 76 NUANCE -> NEUT -21 CC Multi Pi
case 76: return sg * 21;
// 77 NUANCE -> NEUT -41 NC Multi Pi
case 77: return sg * 41;
// 78 NUANCE -> NEUT -41 NC Multi Pi
case 78: return sg * 41;
// 79 NUANCE -> NEUT 21 CC Multi Pi
case 79: return sg * 21;
// 80 NUANCE -> NEUT 21 CC Multi Pi
case 80: return sg * 21;
// 81 NUANCE -> NEUT 41 NC Multi Pi
case 81: return sg * 41;
// 82 NUANCE -> NEUT 41 NC Multi Pi
case 82: return sg * 41;
// 83 NUANCE -> NEUT 41 NC Multi Pi
case 83: return sg * 41;
// 84 NUANCE -> NEUT 41 NC Multi Pi
case 84: return sg * 84;
// 85 NUANCE -> NEUT -21 CC Multi Pi
case 85: return sg * 21;
// 86 NUANCE -> NEUT -21 CC Multi Pi
case 86: return sg * 21;
// 87 NUANCE -> NEUT -41 CC Multi Pi
case 87: return sg * 41;
// 88 NUANCE -> NEUT -41
case 88: return sg * 41;
// 89 NUANCE -> NEUT -41
case 89: return sg * 41;
// 90 NUANCE -> NEUT -41
case 90: return sg * 41;
// 91 NUANCE -> NEUT 26 CC DIS
case 91: return sg * 26;
// 92 NUANCE -> NEUT 46 NC DIS
case 92: return sg * 46;
// 93 NUANCE -> NEUT 17 1#gamma from #Delta
case 93: return sg * 17;
// 94 NUANCE -> NEUT 39 1#gamma from #Delta
case 94: return sg * 39;
// 95 -> UNKOWN NEUT MODE
case 95: return sg * 0;
// 96 NUANCE -> NEUT 36 NC COH
case 96: return sg * 36;
// 97 NUANCE -> NEUT 16
case 97: return sg * 16;
// 98 -> UNKNOWN NEUT MODE
case 98: return sg * 0;
// 99 -> UNKNOWN NEUT MODE
case 99: return sg * 0;
default:
ERR(FTL) << "Unknown Nuance Channel ID = " << ch << std::endl;
throw ("Exiting.");
return 0;
}
return 0;
}
/*
// Notes copied from NuanceChannels.pdf
1 NUANCE CCQE -> NEUT CCQE 1
CC, numu n --> mu- p
Cabibbo-allowed quasi-elastic scattering from nucleons
2 NUANCE NCEL -> NEUT NCEL 51,52 -> Set from whether target is p or n
NC, numu N --> num N, (N=n,p)
(quasi)-elastic scattering from nucleons
3 NUANCE CCPIP -> NEUT CCPIP 11
CC, numu p --> mu- p pi+
resonant single pion production
4 NUANCE CCPI0 -> NEUT CCPI0 = 12
CC, numu n --> mu- p pi0
resonant single pion production
5 NUANCE CCPIPn -> NEUT CCPIPn 13
CC, numu n --> mu- n pi+
resonant single pion production
6 NUANCE NCpPI0 -> NEUT NCpPI0 32
NC, numu p --> numu p pi0
resonant single pion production
7 NUANCE NCpPI+ -> NEUT NCpPI+ 34
NC, numu p --> numu n pi+
resonant single pion production
8 NUANCE NCnPI0 -> NEUT NCnPI0 31
NC, numu n --> numu n pi0
resonant single pion production
9 NUANCE NCnPIM -> NEUT NCnPIM 33
NC, numu n --> numu p pi-
resonant single pion production
10 NUANCE CCPIP -> NEUT CCPIP -11
CC, numubar p --> mu- p pi+
resonant single pion production
11 NUANCE CCPI0 -> NEUT CCPI0 -12
CC, numubar n --> mu- p pi0
resonant single pion production
12 NUANCE CCPIPn -> NEUT CCPIPn -13
CC, numubar n --> mu- n pi+
resonant single pion production
13 NUANCE NCpPI0 -> NEUT NCnPI0 -32
NC, numubar p --> numubar p pi0
resonant single pion production
14 NUANCE NCpPI+ -> NEUT NCpPI+ -34
NC, numubar p --> numubar n pi+
resonant single pion production
15 NUANCE NCnPI0 -> NEUT NCnPI0 -31
NC, numubar n --> numubar n pi0
resonant single pion production
16 NUANCE NCnPIM -> NEUT NCnPIM -33
NC, numubar n --> numubar p pi-
resonant single pion production
17 NUANCE -> NEUT 21 CC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow l^{-} N' multi-#pi
CC, numu p --> mu- Delta+ pi+
resonant processes involving more than a single pion
18 NUANCE -> NEUT 21 CC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow l^{-} N' multi-#pi
CC, numu p --> mu- Delta++ pi0
resonant processes involving more than a single pion
19 NUANCE -> NEUT 21 CC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow l^{-} N' multi-#pi
CC, numu n --> mu- Delta+ pi0
resonant processes involving more than a single pion
20 NUANCE -> NEUT 21 CC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow l^{-} N' multi-#pi
CC, numu n --> mu- Delta0 pi+
resonant processes involving more than a single pion
21 NUANCE -> NEUT 21 CC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow l^{-} N' multi-#pi
CC, numu n --> mu- Delta++ pi-
resonant processes involving more than a single pion
22 NUANCE -> NEUT 41 "NC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow #nu_{l} N multi-#pi"
NC, numu p+ --> numu Delta+ pi0
resonant processes involving more than a single pion
23 NUANCE -> NEUT 41 "NC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow #nu_{l} N multi-#pi"
NC,numu p --> numu Delta0 pi+
resonant processes involving more than a single pion
24 NUANCE -> NEUT 41 "NC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow #nu_{l} N multi-#pi"
NC, numu p --> numu Delta++ pi-
resonant processes involving more than a single pion
25 NUANCE -> NEUT 41 "NC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow #nu_{l} N multi-#pi"
NC, numu n --> numu Delta+ pi-
resonant processes involving more than a single pion
26 NUANCE -> NEUT 41 "NC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow #nu_{l} N multi-#pi"
NC, numu n --> numu Delta0 pi0
resonant processes involving more than a single pion
27 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow #nu_{l} N multi-#pi"
NC, numubar n --> numubar Delta- pi+
resonant processes involving more than a single pion
28 NUANCE -> NEUT -21 CC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow l^{-} N' multi-#pi
CC, numubar p --> mu- Delta+ pi+
resonant processes involving more than a single pion
29 UANCE -> NEUT -21 CC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow l^{-} N' multi-#pi
CC, numubar p --> mu- Delta++ pi0
resonant processes involving more than a single pion
30 NUANCE -> NEUT -21 CC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow l^{-} N' multi-#pi
CC, numubar n --> mu- Delta+ pi0
resonant processes involving more than a single pion
31 NUANCE -> NEUT -21 CC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow l^{-} N' multi-#pi
CC, numubar n --> mu- Delta0 pi+
resonant processes involving more than a single pion
32 NUANCE -> NEUT -21 CC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow l^{-} N' multi-#pi
CC, numubar n --> mu- Delta++ pi-
resonant processes involving more than a single pion
33 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow #nu_{l} N multi-#pi"
NC, numubar p+ --> numubar Delta+ pi0
resonant processes involving more than a single pion
34 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow #nu_{l} N multi-#pi"
NC,numubar p --> numubar Delta0 pi+
resonant processes involving more than a single pion
35 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow #nu_{l} N multi-#pi"
NC, numubar p --> numubar Delta++ pi-
resonant processes involving more than a single pion
36 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow #nu_{l} N multi-#pi"
NC, numubar n --> numubar Delta+ pi-
resonant processes involving more than a single pion
37 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow #nu_{l} N multi-#pi"
NC, numubar n --> numubar Delta0 pi0
resonant processes involving more than a single pion
38 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow #nu_{l} N multi-#pi"
NC, numubar n --> numubar Delta- pi+
resonant processes involving more than a single pion
// RHO Production lumped in with eta production
22 CCeta
43 NCeta on p
42 NCeta on n
39 NUANCE -> NEUT 22
CC, numu p --> mu- p rho+(770)
resonant processes involving more than a single pion
40 NUANCE -> NEUT 22
CC, numu n --> mu- p rho0(770)
resonant processes involving more than a single pion
41 NUANCE -> NEUT 22
CC, numu n --> mu- n rho+(770)
resonant processes involving more than a single pion
42 NUANCE -> NEUT 43
NC, numu p --> numu p rho0(770)
resonant processes involving more than a single pion
43 NUANCE -> NEUT 43
NC, numu p --> numu n rho+(770)
resonant processes involving more than a single pion
44 NUANCE -> NUET 42
NC, numu n --> numu n rho0(770)
resonant processes involving more than a single pion
45 NUANCE -> NEUT -42
NC, numubar n --> numubar p rho-(770)
resonant processes involving more than a single pion
46 NUANCE -> NEUT -22
CC, numubar p --> mu- p rho+(770)
resonant processes involving more than a single pion
47 NUANCE -> NEUT -22
CC, numubar n --> mu- p rho0(770)
resonant processes involving more than a single pion
48 NUANCE -> NEUT -22
CC, numubar n --> mu- n rho+(770)
resonant processes involving more than a single pion
49 NUANCE -> NEUT -43
NC, numubar p --> numubar p rho0(770)
resonant processes involving more than a single pion
50 NUANCE -> NEUT -43
NC, numubar p --> numubar n rho+(770)
resonant processes involving more than a single pion
51 NUANCE -> NEUT -42
NC, numubar n --> numubar n rho0(770)
resonant processes involving more than a single pion
52 NUANCE -> NEUT -42
NC, numubar n --> numubar p rho-(770)
resonant processes involving more than a single pion
53 NUANCE -> NEUT 23 CC 1K: #nu_{l} n #rightarrow l^{-} #Lambda K^{+}
CC, numu p --> mu- Sigma+ K+
resonant processes involving more than a single pion
54 NUANCE -> NEUT 23 CC 1K: #nu_{l} n #rightarrow l^{-} #Lambda K^{+}
CC, numu n --> mu- Sigma0 K+
resonant processes involving more than a single pion
55 NUANCE -> NEUT 23 CC 1K: #nu_{l} n #rightarrow l^{-} #Lambda K^{+}
CC, numu n --> mu- Sigma+ K0
resonant processes involving more than a single pion
56 NUANCE -> NEUT 45 NC 1K: #nu_{l} n #rightarrow #nu_{l} #Lambda K^{+}
NC, numu p --> numu Sigma0 K+
resonant processes involving more than a single pion
57 NUANCE -> NEUT 44 NC 1K: #nu_{l} n #rightarrow #nu_{l} #Lambda K^{0}
NC, numu p --> numu Sigma+ K0
resonant processes involving more than a single pion
58 NUANCE -> NEUT 44 NC 1K: #nu_{l} n #rightarrow #nu_{l} #Lambda K^{0}
NC, numu n --> numu Sigma0 K0
resonant processes involving more than a single pion
59 NUANCE -> NEUT 45 NC 1K: #nu_{l} n #rightarrow #nu_{l} #Lambda K^{+}
NC, numu n --> numu Sigma- K+
resonant processes involving more than a single pion
60 NUANCE -> NEUT -23 CC 1K: #nu_{l} n #rightarrow l^{-} #Lambda K^{+}
CC, numubar p --> mu- Sigma+ K+
resonant processes involving more than a single pion
61 NUANCE -> NEUT -23 CC 1K: #nu_{l} n #rightarrow l^{-} #Lambda K^{+}
CC, numubar n --> mu- Sigma0 K+
resonant processes involving more than a single pion
62 NUANCE -> NEUT -23 CC 1K: #nu_{l} n #rightarrow l^{-} #Lambda K^{+}
CC, numubar n --> mu- Sigma+ K0
resonant processes involving more than a single pion
63 NUANCE -> NEUT -45 NC 1K: #nu_{l} n #rightarrow #nu_{l} #Lambda K^{+}
NC, numubar p --> numubar Sigma0 K+
resonant processes involving more than a single pion
64 NUANCE -> NEUT -44 NC 1K: #nu_{l} n #rightarrow #nu_{l} #Lambda K^{0}
NC, numubar p --> numubar Sigma+ K0
resonant processes involving more than a single pion
65 NUANCE -> NEUT -44 NC 1K: #nu_{l} n #rightarrow #nu_{l} #Lambda K^{0}
NC, numubar n --> numubar Sigma0 K0
resonant processes involving more than a single pion
66 NUANCE -> NEUT -45 NC 1K: #nu_{l} n #rightarrow #nu_{l} #Lambda K^{+}
NC, numubar n --> numubar Sigma- K+
resonant processes involving more than a single pion
67 NUANCE -> NEUT 22
ModeStack[22]->SetTitle("CC1#eta^{0} on n");
CC, numu n --> mu- p eta
resonant processes involving more than a single pion
68 NUANCE -> NEUT 43
NC, numu p --> numu p eta
resonant processes involving more than a single pion
69 NUANCE -> NEUT 42
NC, numu n --> numu n eta
resonant processes involving more than a single pion
70 NUANCE -> NEUT -22
ModeStack[22]->SetTitle("CC1#eta^{0} on n");
CC, numubar n --> mu- p eta
resonant processes involving more than a single pion
71 NUANCE -> NEUT -43
ModeStack[43]->SetTitle("NC1#eta^{0} on p");
NC, numubar p --> numubar p eta
resonant processes involving more than a single pion
72 NUANCE -> NEUT -42
ModeStack[42]->SetTitle("NC1#eta^{0} on n");
NC, numubar n --> numubar n eta
resonant processes involving more than a single pion
73 NUANCE -> NEUT 21
ModeStack[21]->SetTitle("Multi #pi (1.3 < W < 2.0)");
CC, numu n --> mu- K+ Lambda
resonant processes involving more than a single pion
74 NUANCE -> NEUT 41
ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)");
NC, numu p --> numu K+ Lambda
resonant processes involving more than a single pion
75 NUANCE -> NEUT 41
ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)");
NC, numu n --> numu K0 Lambda
resonant processes involving more than a single pion
76 NUANCE -> NEUT -21
ModeStack[21]->SetTitle("Multi #pi (1.3 < W < 2.0)");
CC, numubar n --> mu- K+ Lambda
resonant processes involving more than a single pion
77 NUANCE -> NEUT -41
ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)");
NC, numubar p --> numubar K+ Lambda
resonant processes involving more than a single pion
78 NUANCE -> NEUT -41
ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)");
NC, numubar n --> numubar K0 Lambda
resonant processes involving more than a single pion
CC Multipi ModeStack[21]->SetTitle("Multi #pi (1.3 < W < 2.0)");
NC Multipi ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)");
79 NUANCE -> NEUT 21
CC, numu n --> mu- p pi+ pi-
two pion production
80 NUANCE -> NEUT 21
CC, numu n --> mu- p pi0 pi0
two pion production
81 NUANCE -> NEUT 41
NC, numu p --> numu p pi+ pi-
two pion production
82 NUANCE -> NEUT 41
NC, numu p --> numu p pi0 pi0
two pion production
83 NUANCE -> NEUT 41
NC, numu n --> numu n pi+ pi-
two pion production
84 NUANCE -> NEUT 41
NC, numu n --> numu n pi0 pi0
two pion production
85 NUANCE -> NEUT -21
CC, numubar n --> mu- p pi+ pi-
two pion production
86 NUANCE -> NEUT -21
CC, numubar n --> mu- p pi0 pi0
two pion production
87 NUANCE -> NEUT -41
ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)");
NC, numubar p --> numubar p pi+ pi-
two pion production
88 NUANCE -> NEUT -41
ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)");
NC, numubar p --> numubar p pi0 pi0
two pion production
89 NUANCE -> NEUT -41
ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)");
NC, numubar n --> numubar n pi+ pi-
two pion production
90 NUANCE -> NEUT -41
ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)");
NC, numubar n --> numubar n pi0 pi0
two pion production
91 NUANCE -> NEUT 26
ModeStack[26]->SetTitle("DIS (W > 2.0)");
CC, numu N --> mu- X (where N=n,p)
deep inelastic scattering (nu or nubar)
92 NUANCE -> NEUT 46
ModeStack[46]->SetTitle("DIS (W > 2.0)");
NC, numu N --> numu X (where N=n,p)
deep inelastic scattering (nu or nubar)
93 NUANCE -> NEUT 17 1#gamma from #Delta: #nu_{l} n #rightarrow l^{-} p #gamma
CC, numu n --> mu- p gamma
Delta radiative decay, Delta --> N gamma (only in NUANCE versions v3 and and higher)
94 NUANCE -> NEUT 39 1#gamma from #Delta: #nu_{l} p #rightarrow #nu_{l} p #gamma
neutModeID[15] = 38; neutModeName[15] = "ncngam"; neutModeTitle[15] = "1#gamma from #Delta: #nu_{l} n #rightarrow #nu_{l} n #gamma";
neutModeID[16] = 39; neutModeName[16] = "ncpgam"; neutModeTitle[16] = "1#gamma from #Delta: #nu_{l} p #rightarrow #nu_{l} p #gamma";
NC, numu N --> numu N gamma
Delta radiative decay, Delta --> N gamma (only in NUANCE versions v3 and and higher)
95 -> UNKOWN NEUT MODE
CC, numubar p --> mu+ Lambda, numubar n -- > mu+ Sigma-, numubar p --> mu+ Sigma0
Cabibbo-suppressed QE hyperon production from nucleons
96 NUANCE -> NEUT 36
neutModeID[14] = 36; neutModeName[14] = "nccoh"; neutModeTitle[14] = "NC coherent-#pi: #nu_{l} ^{16}O #rightarrow #nu_{l} ^{16}O #pi^{0}";
NC, numu A --> numu pi0 A
coherent or diffractive pi0 production
97 NUANCE -> NEUT 16
neutModeID[4] = 16; neutModeName[4] = "cccoh"; neutModeTitle[4] = "CC coherent-#pi: #nu_{l} ^{16}O #rightarrow l^{-} ^{16}O #pi^{+}";
CC, numu A --> mu- pi+ A (or numubar A -->
coherent or diffractive pi0 production
98 -> UNKNOWN NEUT MODE
NC, numu e- --> numu e- (or numubar e- -->
neutrino + electron elastic scattering
99 -> UNKNOWN NEUT MODE
CC, numu e- --> mu- nue
neutrino + electron inverse muon decay
NEUT Modes:
// CC Modes
neutModeID[0] = 1; neutModeName[0] = "ccqe"; neutModeTitle[0] = "CCQE: #nu_{l} n #rightarrow l^{-} p";
neutModeID[1] = 11; neutModeName[1] = "ccppip"; neutModeTitle[1] = "CC 1#pi: #nu_{l} p #rightarrow l^{-} p #pi^{+}";
neutModeID[2] = 12; neutModeName[2] = "ccppi0"; neutModeTitle[2] = "CC 1#pi: #nu_{l} n #rightarrow l^{-} p #pi^{0}";
neutModeID[3] = 13; neutModeName[3] = "ccnpip"; neutModeTitle[3] = "CC 1#pi: #nu_{l} n #rightarrow l^{-} n #pi^{+}";
neutModeID[4] = 16; neutModeName[4] = "cccoh"; neutModeTitle[4] = "CC coherent-#pi: #nu_{l} ^{16}O #rightarrow l^{-} ^{16}O #pi^{+}";
neutModeID[5] = 17; neutModeName[5] = "ccgam"; neutModeTitle[5] = "1#gamma from #Delta: #nu_{l} n #rightarrow l^{-} p #gamma";
neutModeID[6] = 21; neutModeName[6] = "ccmpi"; neutModeTitle[6] = "CC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow l^{-} N' multi-#pi";
neutModeID[7] = 22; neutModeName[7] = "cceta"; neutModeTitle[7] = "CC 1#eta: #nu_{l} n #rightarrow l^{-} p #eta";
neutModeID[8] = 23; neutModeName[8] = "cck"; neutModeTitle[8] = "CC 1K: #nu_{l} n #rightarrow l^{-} #Lambda K^{+}";
neutModeID[9] = 26; neutModeName[9] = "ccdis"; neutModeTitle[9] = "CC DIS (2 GeV < W): #nu_{l} N #rightarrow l^{-} N' mesons";
neutModeID[10] = 31; neutModeName[10] = "ncnpi0"; neutModeTitle[10] = "NC 1#pi: #nu_{l} n #rightarrow #nu_{l} n #pi^{0}";
neutModeID[11] = 32; neutModeName[11] = "ncppi0"; neutModeTitle[11] = "NC 1#pi: #nu_{l} p #rightarrow #nu_{l} p #pi^{0}";
neutModeID[12] = 33; neutModeName[12] = "ncppim"; neutModeTitle[12] = "NC 1#pi: #nu_{l} n #rightarrow #nu_{l} p #pi^{-}";
neutModeID[13] = 34; neutModeName[13] = "ncnpip"; neutModeTitle[13] = "NC 1#pi: #nu_{l} p #rightarrow #nu_{l} n #pi^{+}";
neutModeID[14] = 36; neutModeName[14] = "nccoh"; neutModeTitle[14] = "NC coherent-#pi: #nu_{l} ^{16}O #rightarrow #nu_{l} ^{16}O #pi^{0}";
neutModeID[15] = 38; neutModeName[15] = "ncngam"; neutModeTitle[15] = "1#gamma from #Delta: #nu_{l} n #rightarrow #nu_{l} n #gamma";
neutModeID[16] = 39; neutModeName[16] = "ncpgam"; neutModeTitle[16] = "1#gamma from #Delta: #nu_{l} p #rightarrow #nu_{l} p #gamma";
neutModeID[17] = 41; neutModeName[17] = "ncmpi"; neutModeTitle[17] = "NC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow #nu_{l} N multi-#pi";
neutModeID[18] = 42; neutModeName[18] = "ncneta"; neutModeTitle[18] = "NC 1#eta: #nu_{l} n #rightarrow #nu_{l} n #eta";
neutModeID[19] = 43; neutModeName[19] = "ncpeta"; neutModeTitle[19] = "NC 1#eta: #nu_{l} p #rightarrow #nu_{l} p #eta";
neutModeID[20] = 44; neutModeName[20] = "nck0"; neutModeTitle[20] = "NC 1K: #nu_{l} n #rightarrow #nu_{l} #Lambda K^{0}";
neutModeID[21] = 45; neutModeName[21] = "nckp"; neutModeTitle[21] = "NC 1K: #nu_{l} n #rightarrow #nu_{l} #Lambda K^{+}";
neutModeID[22] = 46; neutModeName[22] = "ncdis"; neutModeTitle[22] = "NC DIS (2 GeV < W): #nu_{l} N #rightarrow #nu_{l} N' mesons";
neutModeID[23] = 51; neutModeName[23] = "ncqep"; neutModeTitle[23] = "NC elastic: #nu_{l} p #rightarrow #nu_{l} p";
neutModeID[24] = 52; neutModeName[24] = "ncqen"; neutModeTitle[24] = "NC elastic: #nu_{l} n #rightarrow #nu_{l} n";
*/
#endif
diff --git a/src/InputHandler/NuWroInputHandler.cxx b/src/InputHandler/NuWroInputHandler.cxx
index 62c6b00..f4c62ee 100644
--- a/src/InputHandler/NuWroInputHandler.cxx
+++ b/src/InputHandler/NuWroInputHandler.cxx
@@ -1,477 +1,473 @@
#ifdef __NUWRO_ENABLED__
#include "NuWroInputHandler.h"
#include "InputUtils.h"
NuWroGeneratorInfo::~NuWroGeneratorInfo() { delete fNuWroParticlePDGs; }
void NuWroGeneratorInfo::AddBranchesToTree(TTree* tn) {
tn->Branch("NuWroParticlePDGs", &fNuWroParticlePDGs, "NuWroParticlePDGs/I");
}
void NuWroGeneratorInfo::SetBranchesFromTree(TTree* tn) {
tn->SetBranchAddress("NuWroParticlePDGs", &fNuWroParticlePDGs);
}
void NuWroGeneratorInfo::AllocateParticleStack(int stacksize) {
fNuWroParticlePDGs = new int[stacksize];
}
void NuWroGeneratorInfo::DeallocateParticleStack() {
delete fNuWroParticlePDGs;
}
void NuWroGeneratorInfo::FillGeneratorInfo(event* e) { Reset(); }
void NuWroGeneratorInfo::Reset() {
for (int i = 0; i < kMaxParticles; i++) {
fNuWroParticlePDGs[i] = 0;
}
}
int event1_nof(event* e, int pdg) {
int c = 0;
for (size_t i = 0; i < e->out.size(); i++)
if (e->out[i].pdg == pdg) c++;
return c;
}
NuWroInputHandler::NuWroInputHandler(std::string const& handle,
std::string const& rawinputs) {
LOG(SAM) << "Creating NuWroInputHandler : " << handle << std::endl;
// Run a joint input handling
fName = handle;
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
fSaveExtra = false; // FitPar::Config().GetParB("NuWroSaveExtra");
// Setup the TChain
fNuWroTree = new TChain("treeout");
// Loop over all inputs and grab flux, eventhist, and nevents
std::vector 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) << "nuwro File IsZombie() at " << inputs[inp_it] << std::endl;
throw;
}
// Get Flux/Event hist
TH1D* fluxhist = (TH1D*)inp_file->Get(
(PlotUtils::GetObjectWithName(inp_file, "FluxHist")).c_str());
TH1D* eventhist = (TH1D*)inp_file->Get(
(PlotUtils::GetObjectWithName(inp_file, "EvtHist")).c_str());
if (!fluxhist or !eventhist) {
ERR(FTL) << "nuwro FILE doesn't contain flux/xsec info" << std::endl;
if (FitPar::Config().GetParB("regennuwro")) {
ERR(FTL) << "Regen NuWro has not been added yet. Email the developers!"
<< std::endl;
// ProcessNuWroInputFlux(inputs[inp_it]);
throw;
} else {
ERR(FTL) << "If you would like NUISANCE to generate these for you "
<< "please set parameter regennuwro=1 and re-run."
<< std::endl;
throw;
}
}
// Get N Events
TTree* nuwrotree = (TTree*)inp_file->Get("treeout");
if (!nuwrotree) {
ERR(FTL) << "treeout not located in nuwro file! " << inputs[inp_it]
<< std::endl;
throw;
}
int nevents = nuwrotree->GetEntries();
// Register input to form flux/event rate hists
RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist);
// Add to TChain
fNuWroTree->Add(inputs[inp_it].c_str());
}
// Registor all our file inputs
SetupJointInputs();
// Setup Events
fNuWroEvent = NULL;
fNuWroTree->SetBranchAddress("e", &fNuWroEvent);
fNuWroTree->GetEntry(0);
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->fType = kNUWRO;
fNUISANCEEvent->fNuwroEvent = fNuWroEvent;
fNUISANCEEvent->HardReset();
if (fSaveExtra) {
fNuWroInfo = new NuWroGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fNuWroInfo);
}
};
NuWroInputHandler::~NuWroInputHandler() {
if (fNuWroTree) delete fNuWroTree;
}
void NuWroInputHandler::CreateCache() {
// fNuWroTree->SetCacheEntryRange(0, fNEvents);
// fNuWroTree->AddBranchToCache("*", 1);
// fNuWroTree->SetCacheSize(fCacheSize);
}
void NuWroInputHandler::RemoveCache() {
// fNuWroTree->SetCacheEntryRange(0, fNEvents);
// fNuWroTree->AddBranchToCache("*", 0);
// fNuWroTree->SetCacheSize(0);
}
void NuWroInputHandler::ProcessNuWroInputFlux(const std::string file) {}
FitEvent* NuWroInputHandler::GetNuisanceEvent(const UInt_t 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;
fNuWroTree->GetEntry(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
#ifdef __PROB3PP_ENABLED__
for (size_t i = 0; i < fNUISANCEEvent->fNuwroEvent->in.size(); i++) {
if (std::count(PhysConst::pdg_neutrinos, PhysConst::pdg_neutrinos + 4,
fNUISANCEEvent->fNuwroEvent->in[i].pdg)) {
fNUISANCEEvent->probe_E = fNUISANCEEvent->fNuwroEvent->in[i].t;
fNUISANCEEvent->probe_pdg = fNUISANCEEvent->fNuwroEvent->in[i].pdg;
break;
}
}
#endif
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
#ifdef __USE_NUWRO_SRW_EVENTS__
if (!rwEvs.size()) {
fNuwroParams = fNuWroEvent->par;
}
if (entry >= rwEvs.size()) {
rwEvs.push_back(BaseFitEvt());
rwEvs.back().fNuwroSRWEvent = SRW::SRWEvent(*fNuWroEvent);
rwEvs.back().fNuwroEvent = NULL;
rwEvs.back().fNuwroParams = &fNuwroParams;
rwEvs.back().probe_E = rwEvs.back().fNuwroSRWEvent.NeutrinoEnergy;
rwEvs.back().probe_pdg = rwEvs.back().fNuwroSRWEvent.NeutrinoPDG;
}
#endif
return fNUISANCEEvent;
}
int NuWroInputHandler::ConvertNuwroMode(event* e) {
Int_t proton_pdg, neutron_pdg, pion_pdg, pion_plus_pdg, pion_minus_pdg,
lambda_pdg, eta_pdg, kaon_pdg, kaon_plus_pdg;
proton_pdg = 2212;
eta_pdg = 221;
neutron_pdg = 2112;
pion_pdg = 111;
pion_plus_pdg = 211;
pion_minus_pdg = -211;
// O_16_pdg = 100069; // oznacznie z Neuta
lambda_pdg = 3122;
kaon_pdg = 311;
kaon_plus_pdg = 321;
if (e->flag.qel) // kwiazielastyczne oddziaływanie
{
if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem
{
if (e->flag.cc)
return -1;
else {
if (event1_nof(e, proton_pdg))
return -51;
else if (event1_nof(e, neutron_pdg))
return -52; // sprawdzam dodatkowo ?
}
} else // oddziaływanie z neutrinem
{
if (e->flag.cc)
return 1;
else {
if (event1_nof(e, proton_pdg))
return 51;
else if (event1_nof(e, neutron_pdg))
return 52;
}
}
}
if (e->flag.mec) {
if (e->flag.anty)
return -2;
else
return 2;
}
if (e->flag.res) // rezonansowa produkcja: pojedynczy pion, pojed.eta, kaon,
// multipiony
{
Int_t liczba_pionow, liczba_kaonow;
liczba_pionow = event1_nof(e, pion_pdg) + event1_nof(e, pion_plus_pdg) +
event1_nof(e, pion_minus_pdg);
liczba_kaonow = event1_nof(e, kaon_pdg) + event1_nof(e, kaon_pdg);
if (liczba_pionow > 1 || liczba_pionow == 0) // multipiony
{
if (e->flag.anty) {
if (e->flag.cc)
return -21;
else
return -41;
} else {
if (e->flag.cc)
return 21;
else
return 41;
}
}
if (liczba_pionow == 1) {
if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem
{
if (e->flag.cc) {
if (event1_nof(e, neutron_pdg) && event1_nof(e, pion_minus_pdg))
return -11;
if (event1_nof(e, neutron_pdg) && event1_nof(e, pion_pdg)) return -12;
if (event1_nof(e, proton_pdg) && event1_nof(e, pion_minus_pdg))
return -13;
} else {
if (event1_nof(e, proton_pdg)) {
if (event1_nof(e, pion_minus_pdg))
return -33;
else if (event1_nof(e, pion_pdg))
return -32;
} else if (event1_nof(e, neutron_pdg)) {
if (event1_nof(e, pion_plus_pdg))
return -34;
else if (event1_nof(e, pion_pdg))
return -31;
}
}
} else // oddziaływanie z neutrinem
{
if (e->flag.cc) {
if (event1_nof(e, proton_pdg) && event1_nof(e, pion_plus_pdg))
return 11;
if (event1_nof(e, proton_pdg) && event1_nof(e, pion_pdg)) return 12;
if (event1_nof(e, neutron_pdg) && event1_nof(e, pion_plus_pdg))
return 13;
} else {
if (event1_nof(e, proton_pdg)) {
if (event1_nof(e, pion_minus_pdg))
return 33;
else if (event1_nof(e, pion_pdg))
return 32;
} else if (event1_nof(e, neutron_pdg)) {
if (event1_nof(e, pion_plus_pdg))
return 34;
else if (event1_nof(e, pion_pdg))
return 31;
}
}
}
}
if (event1_nof(e, eta_pdg)) // produkcja rezonansowa ety
{
if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem
{
if (e->flag.cc)
return -22;
else {
if (event1_nof(e, neutron_pdg))
return -42;
else if (event1_nof(e, proton_pdg))
return -43; // sprawdzam dodatkowo ?
}
} else // oddziaływanie z neutrinem
{
if (e->flag.cc)
return 22;
else {
if (event1_nof(e, neutron_pdg))
return 42;
else if (event1_nof(e, proton_pdg))
return 43;
}
}
}
if (event1_nof(e, lambda_pdg) == 1 &&
liczba_kaonow == 1) // produkcja rezonansowa kaonu
{
if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem
{
if (e->flag.cc && event1_nof(e, kaon_pdg))
return -23;
else {
if (event1_nof(e, kaon_pdg))
return -44;
else if (event1_nof(e, kaon_plus_pdg))
return -45;
}
} else // oddziaływanie z neutrinem
{
if (e->flag.cc && event1_nof(e, kaon_plus_pdg))
return 23;
else {
if (event1_nof(e, kaon_pdg))
return 44;
else if (event1_nof(e, kaon_plus_pdg))
return 45;
}
}
}
}
if (e->flag.coh) // koherentne oddziaływanie tylko na O(16)
{
Int_t _target;
_target = e->par.nucleus_p + e->par.nucleus_n; // liczba masowa O(16)
if (_target == 16) {
if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem
{
if (e->flag.cc && event1_nof(e, pion_minus_pdg))
return -16;
else if (event1_nof(e, pion_pdg))
return -36;
} else // oddziaływanie z neutrinem
{
if (e->flag.cc && event1_nof(e, pion_plus_pdg))
return 16;
else if (event1_nof(e, pion_pdg))
return 36;
}
}
}
// gleboko nieelastyczne rozpraszanie
if (e->flag.dis) {
if (e->flag.anty) {
if (e->flag.cc)
return -26;
else
return -46;
} else {
if (e->flag.cc)
return 26;
else
return 46;
}
}
return 9999;
}
void NuWroInputHandler::CalcNUISANCEKinematics() {
// std::cout << "NuWro Event Address " << fNuWroEvent << std::endl;
// Reset all variables
fNUISANCEEvent->ResetEvent();
FitEvent* evt = fNUISANCEEvent;
// Sort Event Info
- evt->fMode = ConvertNuwroMode(fNuWroEvent);
+ evt->Mode = ConvertNuwroMode(fNuWroEvent);
- if (abs(evt->fMode) > 60) {
- evt->fMode = 0;
-
- // Remove failed mode converts
- // return;
+ if (abs(evt->Mode) > 60) {
+ evt->Mode = 0;
}
- evt->Mode = evt->fMode;
evt->fEventNo = 0.0;
evt->fTotCrs = 0.0;
evt->fTargetA = fNuWroEvent->par.nucleus_p + fNuWroEvent->par.nucleus_n;
evt->fTargetZ = fNuWroEvent->par.nucleus_p;
evt->fTargetH = 0;
evt->fBound = (evt->fTargetA) == 1;
// Check Particle Stack
UInt_t npart_in = fNuWroEvent->in.size();
UInt_t npart_out = fNuWroEvent->out.size();
UInt_t npart_post = fNuWroEvent->post.size();
UInt_t npart = npart_in + npart_out + npart_post;
UInt_t kmax = evt->kMaxParticles;
if (npart > kmax) {
ERR(WRN) << "NUWRO has too many particles. Expanding stack." << std::endl;
fNUISANCEEvent->ExpandParticleStack(npart);
}
// Sort Particles
evt->fNParticles = 0;
std::vector::iterator p_iter;
// Initial State
for (p_iter = fNuWroEvent->in.begin(); p_iter != fNuWroEvent->in.end();
p_iter++) {
AddNuWroParticle(fNUISANCEEvent, (*p_iter), kInitialState);
}
// FSI State
// for (size_t i = 0; i < npart_in; i++ ) {
// AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFSIState);
// }
// Final State
for (p_iter = fNuWroEvent->post.begin(); p_iter != fNuWroEvent->post.end();
p_iter++) {
AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFinalState);
}
// Fill Generator Info
if (fSaveExtra) fNuWroInfo->FillGeneratorInfo(fNuWroEvent);
// 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 NuWroInputHandler::AddNuWroParticle(FitEvent* evt, particle& p,
int state) {
// Add Mom
evt->fParticleMom[evt->fNParticles][0] = static_cast(p).x;
evt->fParticleMom[evt->fNParticles][1] = static_cast(p).y;
evt->fParticleMom[evt->fNParticles][2] = static_cast(p).z;
evt->fParticleMom[evt->fNParticles][3] = static_cast(p).t;
// Status/PDG
evt->fParticleState[evt->fNParticles] = state;
evt->fParticlePDG[evt->fNParticles] = p.pdg;
// Add to particle count
evt->fNParticles++;
}
void NuWroInputHandler::Print() {}
#endif
diff --git a/src/InputHandler/SigmaQ0HistogramInputHandler.cxx b/src/InputHandler/SigmaQ0HistogramInputHandler.cxx
index 3db265f..dc6c2ab 100644
--- a/src/InputHandler/SigmaQ0HistogramInputHandler.cxx
+++ b/src/InputHandler/SigmaQ0HistogramInputHandler.cxx
@@ -1,265 +1,264 @@
#include "SigmaQ0HistogramInputHandler.h"
SigmaQ0HistogramInputHandler::SigmaQ0HistogramInputHandler(std::string const& handle, std::string const& rawinputs) {
LOG(SAM) << "Creating SigmaQ0HistogramInputHandler : " << handle << std::endl;
// Run a joint input handling
fName = handle;
// Assign to hist event format
fEventType = kSIGMAQ0HIST;
// Parse our input file
// SIGMAQ0HIST:path.txt,ENERGY=0.25,THETA=0.34,Q0COL=1,XSCOL=2,SCL=10
std::vector parsedinputs = GeneralUtils::ParseToStr(rawinputs, ",");
fFilePath = parsedinputs[0];
// Setup Defaults incase none given
fEnergy = -1.0;
fTheta = -1.0;
fQ0Column = 0;
fSigmaColumn = 1;
fBeamPDG = -1; // Assume electron by default
fDelim = " ";
fScaling = 1.E38;
// Now Loop through and parse possible inputs
for (uint i = 1; i < parsedinputs.size(); i++) {
std::vector parsedspec = GeneralUtils::ParseToStr(parsedinputs[i], "=");
if (parsedspec.size() < 2) {
THROW("NO VALUE GIVEN TO SPECIFIER : " << parsedinputs[0]);
}
std::string spec = parsedspec[0];
std::string value = parsedspec[1];
// Read Energy
if (!spec.compare("ENERGY") or
!spec.compare("E")) {
fEnergy = GeneralUtils::StrToDbl(value);
// Read Theta
} else if (!spec.compare("THETA") or
!spec.compare("T")) {
fTheta = GeneralUtils::StrToDbl(value);
// Set Q0 Column
} else if (!spec.compare("Q0") or
!spec.compare("Q0COL")) {
fQ0Column = GeneralUtils::StrToInt(value);
// Set XSec Column
} else if (!spec.compare("XSCOL") or
!spec.compare("XS") or
!spec.compare("X") or
!spec.compare("SIG") or
!spec.compare("SIGMA")) {
fSigmaColumn = GeneralUtils::StrToInt(value);
} else if (!spec.compare("BEAM") or
!spec.compare("P") or
!spec.compare("PDG")) {
fBeamPDG = GeneralUtils::StrToInt(value);
} else if (!spec.compare("SCL") or
!spec.compare("SC") or
!spec.compare("S")) {
fScaling = GeneralUtils::StrToDbl(value);
// Throw
} else {
THROW("Unknown argument given to SigmaQ0 InputHandler!" );
}
}
// Check we got everything
if (fEnergy == -1 or fEnergy < 0) {
THROW("Energy for SigmaQ0 Handler either invalid or not given! Specify with ',E=ENERGY' in the input spec.");
}
// Default to Electron
if (fBeamPDG == -1) {
ERR(WRN) << "No Beam PDG Given in SigmaQ0 Handler so assuming electron. " <<
"If input is not electron scattering please specify using : PDG=PDGCODE" << std::endl;
}
// Convert Theta
fTheta = fTheta * M_PI / 180.0;
LOG(FIT) << "Set Theta to = " << fTheta << std::endl;
// Have to create a dummy flux and event rate histogram :(
fFluxHist = new TH1D("fluxhist", "fluxhist", 100.0, 0.0, 2.0 * fEnergy);
fFluxHist->SetBinContent( fFluxHist->FindBin(fEnergy), 1.0 );
fEventHist = new TH1D("eventhist", "eventhist", 100.0, 0.0, 2.0 * fEnergy);
fEventHist->SetBinContent( fEventHist->FindBin(fEnergy), 1.0 );
fEnergy = fEnergy * 1.E3;
LOG(FIT) << "Set E energy" << std::endl;
// Now parse the lines in our input file.
fApplyInterpolation = FitPar::Config().GetParB("InterpolateSigmaQ0Histogram");
double interpolation_res = FitPar::Config().GetParD("InterpolateSigmaQ0HistogramRes");
fMaxValue = 0.0;
fMaxX = 0.0;
fMinX = 1.E10;
// Create a TGraph of Points
TGraph* gr = new TGraph();
std::vector inputlines = GeneralUtils::ParseFileToStr(fFilePath, "\n");
for (uint i = 0; i < inputlines.size(); i++) {
std::vector splitline = GeneralUtils::ParseToDbl(inputlines[i], fDelim.c_str());
double q0 = splitline[fQ0Column];
double sig = splitline[fSigmaColumn];
gr->SetPoint(gr->GetN(), q0, sig);
if (sig > fMaxValue) fMaxValue = sig;
if (q0 > fMaxX) fMaxX = q0;
if (q0 < fMinX) fMinX = q0;
}
fInputGraph = new TGraph();
for (int i = 0; i < gr->GetN(); i++){
fInputGraph->SetPoint(fInputGraph->GetN(), gr->GetX()[i], gr->GetY()[i]);
fNEvents++;
if (fApplyInterpolation){
// If not last point create extra events
if (i != gr->GetN()-1){
double xlow = gr->GetX()[i];
double xhigh = gr->GetX()[i+1];
for (int j = 1; j < int(interpolation_res); j++){
double x = xlow + double(j) * (xhigh - xlow) / interpolation_res;
fInputGraph->SetPoint(fInputGraph->GetN(),x, gr->Eval(x));
fNEvents++;
}
}
}
}
delete gr;
fFluxHist->Scale(1.0 / fFluxHist->Integral("width") );
fEventHist->Scale( fScaling * double(fNEvents) / fEventHist->Integral() );
if (fApplyInterpolation){
fEventHist->Scale(1.0 / interpolation_res);
}
fUseAcceptReject = FitPar::Config().GetParB("InterpolateSigmaQ0HistogramThrow");
if (fUseAcceptReject){
std::cout << "USING ACCEPT REJECT" << std::endl;
fEventHist->Scale( fMaxValue / double(fNEvents) );
fNEvents = FitPar::Config().GetParI("InterpolateSigmaQ0HistogramNTHROWS");
std::cout << "NEvents = " << fNEvents << std::endl;
sleep(1);
}
fNUISANCEEvent = new FitEvent();
};
SigmaQ0HistogramInputHandler::~SigmaQ0HistogramInputHandler() {
};
FitEvent* SigmaQ0HistogramInputHandler::GetNuisanceEvent(const UInt_t entry, const bool lightweight) {
// Catch too large entries
if (entry >= (UInt_t)fNEvents) return NULL;
// Evaluate Graph to Create an Event
if (!fUseAcceptReject){
double q0 = fInputGraph->GetX()[entry];
double sig = fInputGraph->GetY()[entry];
FillNuisanceEvent(q0, sig);
} else {
double q0 = ThrowQ0();
FillNuisanceEvent(q0, 1.0);
}
// Return event pointer
return fNUISANCEEvent;
}
double SigmaQ0HistogramInputHandler::ThrowQ0(){
// Use input graph to Throw Q0.
int count = 0;
while (count < 1E7){
double x = fRandom.Uniform(fMinX,fMaxX);
if (fRandom.Uniform(0.0,1.0) <= fInputGraph->Eval(x) / fMaxValue){ return x; };
std::cout << "THROW " << count << " : " << x << std::endl;
}
return 0.0;
}
void SigmaQ0HistogramInputHandler::FillNuisanceEvent(double q0, double sig) {
// Reset all variables
fNUISANCEEvent->ResetEvent();
// Fill Globals
- fNUISANCEEvent->fMode = 1; // Assume CCQE for now...
fNUISANCEEvent->Mode = 1;
fNUISANCEEvent->fEventNo = 0;
fNUISANCEEvent->fTargetA = 0; // Should the User Specify these?
fNUISANCEEvent->fTargetZ = 0; // Should the User Specify these?
fNUISANCEEvent->fTargetH = 0; // Should the User Specify these?
fNUISANCEEvent->fBound = 1; // Should the User Specify these?
fNUISANCEEvent->InputWeight = sig;
// Add incoming beam particle along Z with energy E and outgoing with theta
if (fBeamPDG == 11) {
// Initial Beam
fNUISANCEEvent->fParticleState[0] = kInitialState;
fNUISANCEEvent->fParticlePDG[0] = 11;
double mass = 0.511;
// Get Momentum of Electron moving along Z
fNUISANCEEvent->fParticleMom[0][0] = 0.0;
fNUISANCEEvent->fParticleMom[0][1] = 0.0;
LOG(FIT) << "Setting Energy = " << fEnergy << std::endl;
fNUISANCEEvent->fParticleMom[0][2] = sqrt(fEnergy * fEnergy - mass * mass);
fNUISANCEEvent->fParticleMom[0][3] = fEnergy;
// Outgoing particle
fNUISANCEEvent->fParticleState[1] = kFinalState;
fNUISANCEEvent->fParticlePDG[1] = 11;
// Get Momentum of Electron outgoing minus q0
double oute = fEnergy - q0;
double outp = sqrt(oute * oute - mass * mass);
fNUISANCEEvent->fParticleMom[1][0] = 0.0;
fNUISANCEEvent->fParticleMom[1][1] = sin(fTheta) * outp;
fNUISANCEEvent->fParticleMom[1][2] = cos(fTheta) * outp;
fNUISANCEEvent->fParticleMom[1][3] = oute;
}
// Update Particles
fNUISANCEEvent->fNParticles = 2;
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent-> OrderStack();
// Check Q0
// FitParticle* ein = fNUISANCEEvent->PartInfo(0);
// FitParticle* eout = fNUISANCEEvent->GetHMFSParticle(11);
// double newq0 = fabs(ein->fP.E() - eout->fP.E()) / 1000.0;
// double E = ein->fP.E() / 1000.0;
// double theta = ein->fP.Vect().Angle(eout->fP.Vect()) * 180. / M_PI;
// LOG(FIT) << "Extracted event from line: theirs-" << q0/1.E3 << " ours-" << newq0 << " E-" << fEnergy << " T-" << theta << " X-" << sig << std::endl;
return;
}
diff --git a/src/Reweight/FitWeight.cxx b/src/Reweight/FitWeight.cxx
index 56939c4..3559a3c 100644
--- a/src/Reweight/FitWeight.cxx
+++ b/src/Reweight/FitWeight.cxx
@@ -1,273 +1,274 @@
#include "FitWeight.h"
#include "GENIEWeightEngine.h"
#include "LikelihoodWeightEngine.h"
#include "ModeNormEngine.h"
#include "NEUTWeightEngine.h"
#include "NIWGWeightEngine.h"
#include "NUISANCEWeightEngine.h"
#include "NuWroWeightEngine.h"
#include "OscWeightEngine.h"
#include "SampleNormEngine.h"
#include "SplineWeightEngine.h"
#include "T2KWeightEngine.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;
case kOSCILLATION:
fAllRW[type] = new OscWeightEngine();
break;
case kMODENORM:
fAllRW[type] = new ModeNormEngine();
+ break;
default:
THROW("CANNOT ADD RW Engine for unknown dial type: " << type);
break;
}
}
WeightEngineBase* FitWeight::GetRWEngine(int type) {
if (HasRWEngine(type)) {
return fAllRW[type];
}
THROW("CANNOT get RW Engine for dial type: " << type);
}
bool FitWeight::HasRWEngine(int type) {
switch (type) {
case kNEUT:
case kNUWRO:
case kGENIE:
case kNORM:
case kLIKEWEIGHT:
case kT2K:
case kCUSTOM:
case kSPLINEPARAMETER:
case kNIWG:
case kOSCILLATION: {
return fAllRW.count(type);
}
default: { THROW("CANNOT get RW Engine for dial type: " << type); }
}
}
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);
if (nuisenum == -1) {
THROW("NUISENUM == " << nuisenum << " for " << name);
}
// 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::iterator iter = fAllRW.begin();
iter != fAllRW.end(); iter++) {
(*iter).second->Reconfigure(silent);
}
}
void FitWeight::SetDialValue(std::string name, double val) {
// Add extra check, if name not found look for one with name in it.
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 = Reweight::GetDialType(nuisenum);
if (fAllRW.find(dialtype) == fAllRW.end()) {
THROW("Cannot find RW Engine for dialtype = "
<< dialtype << ", " << Reweight::RemoveDialType(nuisenum));
}
// 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) {
// Add extra check, if name not found look for one with name in it.
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::iterator iter = fAllRW.begin();
iter != fAllRW.end(); iter++) {
double w = (*iter).second->CalcWeight(evt);
rwweight *= w;
}
return rwweight;
}
void FitWeight::UpdateWeightEngine(const double* x) {
size_t count = 0;
for (std::vector::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::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 + "_norm") != fAllEnums.end()) {
if (fAllValues.find(fAllEnums[name + "_norm"]) != fAllValues.end()) {
return fAllValues[fAllEnums[name + "_norm"]];
} 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;
}
}
diff --git a/src/Reweight/WeightUtils.cxx b/src/Reweight/WeightUtils.cxx
index b6787f3..a17eb22 100644
--- a/src/Reweight/WeightUtils.cxx
+++ b/src/Reweight/WeightUtils.cxx
@@ -1,614 +1,615 @@
#include "WeightUtils.h"
#include "FitLogger.h"
#ifdef __T2KREW_ENABLED__
#include "T2KGenieReWeight.h"
#include "T2KNIWGReWeight.h"
#include "T2KNIWGUtils.h"
#include "T2KNeutReWeight.h"
#include "T2KNeutUtils.h"
#include "T2KReWeight.h"
using namespace t2krew;
#endif
#ifdef __NIWG_ENABLED__
#include "NIWGReWeight.h"
#include "NIWGReWeight1piAngle.h"
#include "NIWGReWeight2010a.h"
#include "NIWGReWeight2012a.h"
#include "NIWGReWeight2014a.h"
#include "NIWGReWeightDeltaMass.h"
#include "NIWGReWeightEffectiveRPA.h"
#include "NIWGReWeightHadronMultSwitch.h"
#include "NIWGReWeightMEC.h"
#include "NIWGReWeightPiMult.h"
#include "NIWGReWeightProtonFSIbug.h"
#include "NIWGReWeightRPA.h"
#include "NIWGReWeightSpectralFunc.h"
#include "NIWGReWeightSplineEnu.h"
#include "NIWGSyst.h"
#include "NIWGSystUncertainty.h"
#endif
#ifdef __NEUT_ENABLED__
#include "NReWeight.h"
#include "NReWeightCasc.h"
#include "NReWeightNuXSecCCQE.h"
#include "NReWeightNuXSecCCRES.h"
#include "NReWeightNuXSecCOH.h"
#include "NReWeightNuXSecDIS.h"
#include "NReWeightNuXSecNC.h"
#include "NReWeightNuXSecNCEL.h"
#include "NReWeightNuXSecNCRES.h"
#include "NReWeightNuXSecRES.h"
#include "NReWeightNuclPiless.h"
#include "NSyst.h"
#include "NSystUncertainty.h"
#include "neutpart.h"
#include "neutvect.h"
#endif
#ifdef __NUWRO_ENABLED__
#include "event1.h"
#endif
#ifdef __NUWRO_REWEIGHT_ENABLED__
#include "NuwroReWeight.h"
#include "NuwroReWeight_FlagNorm.h"
#include "NuwroReWeight_QEL.h"
#include "NuwroReWeight_SPP.h"
#include "NuwroSyst.h"
#include "NuwroSystUncertainty.h"
#endif
#ifdef __GENIE_ENABLED__
#include "EVGCore/EventRecord.h"
#include "EVGCore/EventRecord.h"
#include "GHEP/GHepRecord.h"
#include "GSyst.h"
#include "GSystUncertainty.h"
#include "Ntuple/NtpMCEventRecord.h"
#include "ReWeight/GReWeight.h"
#include "ReWeight/GReWeightAGKY.h"
#include "ReWeight/GReWeightDISNuclMod.h"
#include "ReWeight/GReWeightFGM.h"
#include "ReWeight/GReWeightFZone.h"
#include "ReWeight/GReWeightINuke.h"
#include "ReWeight/GReWeightNonResonanceBkg.h"
#include "ReWeight/GReWeightNuXSecCCQE.h"
#include "ReWeight/GReWeightNuXSecCCQEvec.h"
#include "ReWeight/GReWeightNuXSecCCRES.h"
#include "ReWeight/GReWeightNuXSecCOH.h"
#include "ReWeight/GReWeightNuXSecDIS.h"
#include "ReWeight/GReWeightNuXSecNC.h"
#include "ReWeight/GReWeightNuXSecNCEL.h"
#include "ReWeight/GReWeightNuXSecNCRES.h"
#include "ReWeight/GReWeightResonanceDecay.h"
using namespace genie;
using namespace genie::rew;
#endif
#include "GlobalDialList.h"
#include "ModeNormEngine.h"
#include "NUISANCESyst.h"
#include "OscWeightEngine.h"
//********************************************************************
TF1 FitBase::GetRWConvFunction(std::string const &type,
std::string const &name) {
//********************************************************************
std::string dialfunc = "x";
std::string parType = type;
double low = -10000.0;
double high = 10000.0;
if (parType.find("parameter") == std::string::npos) parType += "_parameter";
std::string line;
ifstream card(
(GeneralUtils::GetTopLevelDir() + "/parameters/dial_conversion.card")
.c_str(),
ifstream::in);
while (std::getline(card >> std::ws, line, '\n')) {
std::vector inputlist = GeneralUtils::ParseToStr(line, " ");
// Check the line length
if (inputlist.size() < 4) continue;
// Check whether this is a comment
if (inputlist[0].c_str()[0] == '#') continue;
// Check whether this is the correct parameter type
if (inputlist[0].compare(parType) != 0) continue;
// Check the parameter name
if (inputlist[1].compare(name) != 0) continue;
// inputlist[2] should be the units... ignore for now
dialfunc = inputlist[3];
// High and low are optional, check whether they exist
if (inputlist.size() > 4) low = GeneralUtils::StrToDbl(inputlist[4]);
if (inputlist.size() > 5) high = GeneralUtils::StrToDbl(inputlist[5]);
}
TF1 convfunc = TF1((name + "_convfunc").c_str(), dialfunc.c_str(), low, high);
return convfunc;
}
//********************************************************************
std::string FitBase::GetRWUnits(std::string const &type,
std::string const &name) {
//********************************************************************
std::string unit = "sig.";
std::string parType = type;
if (parType.find("parameter") == std::string::npos) {
parType += "_parameter";
}
std::string line;
std::ifstream card(
(GeneralUtils::GetTopLevelDir() + "/parameters/dial_conversion.card")
.c_str(),
ifstream::in);
while (std::getline(card >> std::ws, line, '\n')) {
std::vector inputlist = GeneralUtils::ParseToStr(line, " ");
// Check the line length
if (inputlist.size() < 3) continue;
// Check whether this is a comment
if (inputlist[0].c_str()[0] == '#') continue;
// Check whether this is the correct parameter type
if (inputlist[0].compare(parType) != 0) continue;
// Check the parameter name
if (inputlist[1].compare(name) != 0) continue;
unit = inputlist[2];
break;
}
return unit;
}
//********************************************************************
double FitBase::RWAbsToSigma(std::string const &type, std::string const &name,
double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.GetX(val);
if (fabs(conv_val) < 1E-10) conv_val = 0.0;
+
+ std::cout << "AbsToSigma(" << name << ") = " << val << " -> " << conv_val
+ << std::endl;
return conv_val;
}
//********************************************************************
double FitBase::RWSigmaToAbs(std::string const &type, std::string const &name,
double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.Eval(val);
return conv_val;
}
//********************************************************************
double FitBase::RWFracToSigma(std::string const &type, std::string const &name,
double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.GetX((val * f1.Eval(0.0)));
if (fabs(conv_val) < 1E-10) conv_val = 0.0;
return conv_val;
}
//********************************************************************
double FitBase::RWSigmaToFrac(std::string const &type, std::string const &name,
double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.Eval(val) / f1.Eval(0.0);
return conv_val;
}
int FitBase::ConvDialType(std::string const &type) {
if (!type.compare("neut_parameter"))
return kNEUT;
else if (!type.compare("niwg_parameter"))
return kNIWG;
else if (!type.compare("nuwro_parameter"))
return kNUWRO;
else if (!type.compare("t2k_parameter"))
return kT2K;
else if (!type.compare("genie_parameter"))
return kGENIE;
else if (!type.compare("custom_parameter"))
return kCUSTOM;
else if (!type.compare("norm_parameter"))
return kNORM;
else if (!type.compare("likeweight_parameter"))
return kLIKEWEIGHT;
else if (!type.compare("spline_parameter"))
return kSPLINEPARAMETER;
else if (!type.compare("osc_parameter"))
return kOSCILLATION;
else if (!type.compare("modenorm_parameter"))
return kMODENORM;
else
return kUNKNOWN;
}
std::string FitBase::ConvDialType(int type) {
switch (type) {
case kNEUT: {
return "neut_parameter";
}
case kNIWG: {
return "niwg_parameter";
}
case kNUWRO: {
return "nuwro_parameter";
}
case kT2K: {
return "t2k_parameter";
}
case kGENIE: {
return "genie_parameter";
}
case kNORM: {
return "norm_parameter";
}
case kCUSTOM: {
return "custom_parameter";
}
case kLIKEWEIGHT: {
return "likeweight_parameter";
}
case kSPLINEPARAMETER: {
return "spline_parameter";
}
case kOSCILLATION: {
return "osc_parameter";
}
case kMODENORM: {
return "modenorm_parameter";
}
default:
return "unknown_parameter";
}
}
int FitBase::GetDialEnum(std::string const &type, std::string const &name) {
return FitBase::GetDialEnum(FitBase::ConvDialType(type), name);
}
int FitBase::GetDialEnum(int type, std::string const &name) {
int offset = type * 1000;
int this_enum = Reweight::kNoDialFound; // Not Found
std::cout << "Getting dial enum " << type << " " << name << std::endl;
// Select Types
switch (type) {
// NEUT DIAL TYPE
case kNEUT: {
#ifdef __NEUT_ENABLED__
int neut_enum = (int)neut::rew::NSyst::FromString(name);
if (neut_enum != 0) {
this_enum = neut_enum + offset;
}
#else
this_enum = Reweight::kNoTypeFound; // Not enabled
#endif
break;
}
// NIWG DIAL TYPE
case kNIWG: {
#ifdef __NIWG_ENABLED__
int niwg_enum = (int)niwg::rew::NIWGSyst::FromString(name);
if (niwg_enum != 0) {
this_enum = niwg_enum + offset;
}
#else
this_enum = Reweight::kNoTypeFound;
#endif
break;
}
// NUWRO DIAL TYPE
case kNUWRO: {
#ifdef __NUWRO_REWEIGHT_ENABLED__
int nuwro_enum = (int)nuwro::rew::NuwroSyst::FromString(name);
if (nuwro_enum > 0) {
this_enum = nuwro_enum + offset;
}
#else
this_enum = Reweight::kNoTypeFound;
#endif
}
// GENIE DIAL TYPE
case kGENIE: {
#ifdef __GENIE_ENABLED__
int genie_enum = (int)genie::rew::GSyst::FromString(name);
if (genie_enum > 0) {
this_enum = genie_enum + offset;
}
#else
this_enum = Reweight::kNoTypeFound;
#endif
break;
}
case kCUSTOM: {
int custom_enum = 0; // PLACEHOLDER
this_enum = custom_enum + offset;
break;
}
// T2K DIAL TYPE
case kT2K: {
#ifdef __T2KREW_ENABLED__
int t2k_enum = (int)t2krew::T2KSyst::FromString(name);
if (t2k_enum > 0) {
this_enum = t2k_enum + offset;
}
#else
this_enum = Reweight::kNoTypeFound;
#endif
break;
}
case kNORM: {
if (gNormEnums.find(name) == gNormEnums.end()) {
gNormEnums[name] = gNormEnums.size() + 1 + offset;
}
this_enum = gNormEnums[name];
break;
}
case kLIKEWEIGHT: {
if (gLikeWeightEnums.find(name) == gLikeWeightEnums.end()) {
gLikeWeightEnums[name] = gLikeWeightEnums.size() + 1 + offset;
}
this_enum = gLikeWeightEnums[name];
break;
}
case kSPLINEPARAMETER: {
if (gSplineParameterEnums.find(name) == gSplineParameterEnums.end()) {
gSplineParameterEnums[name] = gSplineParameterEnums.size() + 1 + offset;
}
this_enum = gSplineParameterEnums[name];
break;
}
case kOSCILLATION: {
#ifdef __PROB3PP_ENABLED__
int oscEnum = OscWeightEngine::SystEnumFromString(name);
if (oscEnum != 0) {
this_enum = oscEnum + offset;
}
#else
this_enum = Reweight::kNoTypeFound; // Not enabled
#endif
}
case kMODENORM: {
size_t us_pos = name.find_first_of('_');
std::string numstr = name.substr(us_pos + 1);
int mode_num = std::atoi(numstr.c_str());
LOG(FTL) << "Getting mode num " << mode_num << std::endl;
if (!mode_num) {
ERR(FTL) << "Attempting to parse dial name: \"" << name
<< "\" as a mode norm dial but failed." << std::endl;
throw;
}
this_enum = 60 + mode_num + offset;
break;
}
}
// If Not Enabled
if (this_enum == Reweight::kNoTypeFound) {
ERR(FTL) << "RW Engine not supported for " << FitBase::ConvDialType(type)
<< std::endl;
ERR(FTL) << "Check dial " << name << std::endl;
}
// If Not Found
if (this_enum == Reweight::kNoDialFound) {
ERR(FTL) << "Dial " << name << " not found." << std::endl;
}
return this_enum;
}
int Reweight::ConvDialType(std::string const &type) {
return FitBase::ConvDialType(type);
}
std::string Reweight::ConvDialType(int type) {
return FitBase::ConvDialType(type);
}
int Reweight::GetDialType(int type) {
int t = (type / 1000);
return t > kMODENORM ? Reweight::kNoDialFound : t;
}
-int Reweight::RemoveDialType(int type){
- return (type%1000);
-}
+int Reweight::RemoveDialType(int type) { return (type % 1000); }
int Reweight::NEUTEnumFromName(std::string const &name) {
#ifdef __NEUT_ENABLED__
int neutenum = (int)neut::rew::NSyst::FromString(name);
return (neutenum > 0) ? neutenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::NIWGEnumFromName(std::string const &name) {
#ifdef __NIWG_ENABLED__
int niwgenum = (int)niwg::rew::NIWGSyst::FromString(name);
return (niwgenum != 0) ? niwgenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::NUWROEnumFromName(std::string const &name) {
#ifdef __NUWRO_REWEIGHT_ENABLED__
int nuwroenum = (int)nuwro::rew::NuwroSyst::FromString(name);
return (nuwroenum > 0) ? nuwroenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::GENIEEnumFromName(std::string const &name) {
#ifdef __GENIE_ENABLED__
int genieenum = (int)genie::rew::GSyst::FromString(name);
return (genieenum > 0) ? genieenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::T2KEnumFromName(std::string const &name) {
#ifdef __T2KREW_ENABLED__
int t2kenum = (int)t2krew::T2KSyst::FromString(name);
return (t2kenum > 0) ? t2kenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::OscillationEnumFromName(std::string const &name) {
#ifdef __PROB3PP_ENABLED__
int oscEnum = OscWeightEngine::SystEnumFromString(name);
return (oscEnum > 0) ? oscEnum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::NUISANCEEnumFromName(std::string const &name, int type) {
int nuisenum = Reweight::DialList().EnumFromNameAndType(name, type);
return nuisenum;
}
int Reweight::CustomEnumFromName(std::string const &name) {
int custenum = Reweight::ConvertNUISANCEDial(name);
return custenum;
}
int Reweight::ConvDial(std::string const &name, std::string const &type,
bool exceptions) {
return Reweight::ConvDial(name, Reweight::ConvDialType(type), exceptions);
}
int Reweight::ConvDial(std::string const &fullname, int type, bool exceptions) {
std::string name =
GeneralUtils::ParseToStr(fullname, ",")[0]; // Only use first dial given
// Produce offset seperating each type.
int offset = type * 1000;
int genenum = Reweight::kNoDialFound;
switch (type) {
case kNEUT:
genenum = NEUTEnumFromName(name);
break;
case kNIWG:
genenum = NIWGEnumFromName(name);
break;
case kNUWRO:
genenum = NUWROEnumFromName(name);
break;
case kGENIE:
genenum = GENIEEnumFromName(name);
break;
case kT2K:
genenum = T2KEnumFromName(name);
break;
case kCUSTOM:
genenum = CustomEnumFromName(name);
break;
case kNORM:
case kLIKEWEIGHT:
case kSPLINEPARAMETER:
case kNEWSPLINE:
genenum = NUISANCEEnumFromName(name, type);
break;
case kOSCILLATION:
genenum = OscillationEnumFromName(name);
break;
case kMODENORM:
genenum = ModeNormEngine::SystEnumFromString(name);
break;
default:
genenum = Reweight::kNoTypeFound;
break;
}
// Throw if required.
if (exceptions) {
// If Not Enabled
if (genenum == Reweight::kGeneratorNotBuilt) {
ERR(FTL) << "RW Engine not supported for " << FitBase::ConvDialType(type)
<< std::endl;
ERR(FTL) << "Check dial " << name << std::endl;
throw;
}
// If no type enabled
if (genenum == Reweight::kNoTypeFound) {
ERR(FTL) << "Type mismatch inside ConvDialEnum" << std::endl;
throw;
}
// If Not Found
if (genenum == Reweight::kNoDialFound) {
ERR(FTL) << "Dial " << name << " not found." << std::endl;
throw;
}
}
// Add offset if no issue
int nuisenum = genenum;
if ((genenum != Reweight::kGeneratorNotBuilt) &&
(genenum != Reweight::kNoTypeFound) &&
(genenum != Reweight::kNoDialFound)) {
nuisenum += offset;
}
// Now register dial
Reweight::DialList().RegisterDialEnum(name, type, nuisenum);
return nuisenum;
}
std::string Reweight::ConvDial(int nuisenum) {
for (size_t i = 0; i < Reweight::DialList().fAllDialEnums.size(); i++) {
if (Reweight::DialList().fAllDialEnums[i] == nuisenum) {
return Reweight::DialList().fAllDialNames[i];
}
}
LOG(FIT) << "Cannot find dial with enum = " << nuisenum << std::endl;
return "";
}
diff --git a/src/Tests/SignalDefTests.cxx b/src/Tests/SignalDefTests.cxx
index 70d762b..551f66b 100644
--- a/src/Tests/SignalDefTests.cxx
+++ b/src/Tests/SignalDefTests.cxx
@@ -1,451 +1,447 @@
#include
#include
#include "TLorentzVector.h"
#include "TRandom3.h"
#include "FitEvent.h"
#include "SignalDef.h"
struct ConstructibleFitEvent : public FitEvent {
ConstructibleFitEvent() : FitEvent() { fNParticles = 0; }
void AddPart(double Mom[4], size_t State, int PDG) {
fParticleMom[fNParticles][0] = Mom[0];
fParticleMom[fNParticles][1] = Mom[1];
fParticleMom[fNParticles][2] = Mom[2];
fParticleMom[fNParticles][3] = Mom[3];
fParticleState[fNParticles] = State;
fParticlePDG[fNParticles] = PDG;
fNParticles++;
}
- void SetMode(int mode) {
- fMode = mode;
- Mode = mode;
- }
+ void SetMode(int mode) { Mode = mode; }
std::string ToString() {
std::stringstream ss("");
- ss << "Mode: " << fMode << std::endl;
+ ss << "Mode: " << Mode << std::endl;
ss << "Particles: " << fNParticles << std::endl;
ss << " -> Particle Stack " << std::endl;
for (int i = 0; i < fNParticles; i++) {
ss << " -> -> " << i << ". " << fParticlePDG[i] << " "
<< fParticleState[i] << " "
<< " Mom(" << fParticleMom[i][0] << ", " << fParticleMom[i][1] << ", "
<< fParticleMom[i][2] << ", " << fParticleMom[i][3] << ")."
<< std::endl;
}
return ss.str();
}
};
template
ConstructibleFitEvent MakePDGStackEvent(int (&ISpdgs)[N], int (&FSpdgs)[M],
int Mode = 1) {
ConstructibleFitEvent fe;
double MomHolder[4] = {0, 0, 1E3, 1E3};
for (size_t p_it = 0; p_it < N; ++p_it) {
fe.AddPart(MomHolder, kInitialState, ISpdgs[p_it]);
}
TRandom3 rnd;
TLorentzVector rnd4M;
TVector3 rn3M;
for (size_t p_it = 0; p_it < M; ++p_it) {
/// Could do better and actually get the correct masses...
rn3M.SetMagThetaPhi(fabs(rnd.Gaus(200, 50)), rnd.Uniform(M_PI),
2 * rnd.Uniform(M_PI));
rnd4M.SetVectM(rn3M, 105);
MomHolder[0] = rnd4M[0];
MomHolder[1] = rnd4M[1];
MomHolder[2] = rnd4M[2];
MomHolder[3] = rnd4M[3];
fe.AddPart(MomHolder, kFinalState, FSpdgs[p_it]);
}
fe.SetMode(Mode);
fe.OrderStack();
return fe;
}
int main(int argc, char const *argv[]) {
-
- bool FailOnFail = (argc>1);
+ bool FailOnFail = (argc > 1);
LOG_VERB(SAM);
LOG(FIT) << "* Running SignalDef Tests" << std::endl;
LOG(FIT) << "***************************************************"
<< std::endl;
int IS[] = {14};
int FS_CC0pi_1[] = {13, 2112, 2212, 2112, 2212};
ConstructibleFitEvent fe_CC0pi_1 = MakePDGStackEvent(IS, FS_CC0pi_1);
int FS_CC0pi_2[] = {13, 2112, 2212, 2112, 2212};
ConstructibleFitEvent fe_CC0pi_2 = MakePDGStackEvent(IS, FS_CC0pi_2, 2);
int FS_CC0pi_3[] = {-13, 2112, 2212, 2112, 2212};
ConstructibleFitEvent fe_CC0pi_3 = MakePDGStackEvent(IS, FS_CC0pi_3);
int FS_CC0pi_4[] = {13, 2112, 2212};
ConstructibleFitEvent fe_CC0pi_4 = MakePDGStackEvent(IS, FS_CC0pi_4, 12);
int FS_CC1pip_1[] = {13, 2212, 2112, 211};
ConstructibleFitEvent fe_CC1pip_1 = MakePDGStackEvent(IS, FS_CC1pip_1, 2);
int FS_CC1pim_1[] = {13, -211, 2212, 2112};
ConstructibleFitEvent fe_CC1pim_1 = MakePDGStackEvent(IS, FS_CC1pim_1, 11);
int FS_CC1pi0_1[] = {13, 2212, 111, 2112};
ConstructibleFitEvent fe_CC1pi0_1 = MakePDGStackEvent(IS, FS_CC1pi0_1, 12);
int FS_CC1pi0_2[] = {11, 2212, 111, 2112};
ConstructibleFitEvent fe_CC1pi0_2 = MakePDGStackEvent(IS, FS_CC1pi0_2, 12);
int FS_CCNpi_1[] = {13, 2212, 111, 2112, 211};
ConstructibleFitEvent fe_CCNpi_1 = MakePDGStackEvent(IS, FS_CCNpi_1, 21);
int FS_CCNpi_2[] = {13, -211, 211, 2112, 211};
ConstructibleFitEvent fe_CCNpi_2 = MakePDGStackEvent(IS, FS_CCNpi_2, 21);
int FS_CCNpi_3[] = {13, 2212, 111, 211, -211};
ConstructibleFitEvent fe_CCNpi_3 = MakePDGStackEvent(IS, FS_CCNpi_3, 26);
int FS_CCNpi_4[] = {13, 2212, 111, 111, 111};
ConstructibleFitEvent fe_CCNpi_4 = MakePDGStackEvent(IS, FS_CCNpi_4, 26);
int FS_CCCOH_1[] = {13, 211};
ConstructibleFitEvent fe_CCCOH_1 = MakePDGStackEvent(IS, FS_CCCOH_1, 16);
int FS_NCel_1[] = {14, 2112};
ConstructibleFitEvent fe_NCel_1 = MakePDGStackEvent(IS, FS_NCel_1, 52);
int FS_NCel_2[] = {12, 2212};
ConstructibleFitEvent fe_NCel_2 = MakePDGStackEvent(IS, FS_NCel_2, 51);
int FS_NC1pi_1[] = {14, 2112, -211};
ConstructibleFitEvent fe_NC1pi_1 = MakePDGStackEvent(IS, FS_NC1pi_1, 31);
int FS_NCNpi_1[] = {14, 2212, 211};
ConstructibleFitEvent fe_NCNpi_1 = MakePDGStackEvent(IS, FS_NCNpi_1, 32);
LOG(FIT) << "* Testing: SignalDef::isCCINC" << std::endl;
std::map isCCINC_PassExpectations;
isCCINC_PassExpectations[&fe_CC0pi_1] = true; // numu CC0pi
isCCINC_PassExpectations[&fe_CC0pi_2] = true; // numu CC0pi (2p2h)
isCCINC_PassExpectations[&fe_CC0pi_3] = false; // numub CC0pi
isCCINC_PassExpectations[&fe_CC0pi_4] = true; // numu CC0pi (RES)
isCCINC_PassExpectations[&fe_CC1pip_1] = true; // numu CC1pip (2p2h)
isCCINC_PassExpectations[&fe_CC1pim_1] = true; // numu CC1pim
isCCINC_PassExpectations[&fe_CC1pi0_1] = true; // numu CC1pi0
isCCINC_PassExpectations[&fe_CC1pi0_2] = false; // nue CC1pi0
isCCINC_PassExpectations[&fe_CCNpi_1] = true; // numu CC multi pi
isCCINC_PassExpectations[&fe_CCNpi_2] = true; // numu CC multi pi
isCCINC_PassExpectations[&fe_CCNpi_3] = true; // numu CC multi pi
isCCINC_PassExpectations[&fe_CCNpi_4] = true; // numu CC multi pi
isCCINC_PassExpectations[&fe_CCCOH_1] = true; // numu CC COH pi
isCCINC_PassExpectations[&fe_NCel_1] = false; // numu NCEl
isCCINC_PassExpectations[&fe_NCel_2] = false; // nue NCEl
isCCINC_PassExpectations[&fe_NC1pi_1] = false; // numu NC1pi
isCCINC_PassExpectations[&fe_NCNpi_1] = false; // numu NC multi pi
size_t ctr = 0;
for (std::map::iterator
fe_it = isCCINC_PassExpectations.begin();
fe_it != isCCINC_PassExpectations.end(); ++fe_it, ++ctr) {
bool res = SignalDef::isCCINC(fe_it->first, 14);
if (res != fe_it->second) {
ERR(FTL) << "Event: (" << ctr << ")\n"
<< fe_it->first->ToString() << std::endl;
ERR(FTL) << (res ? "passed" : "failed")
<< " SignalDef::isCCINC unexpectedly." << std::endl;
} else {
LOG(SAM) << "Event: (" << ctr << ") " << (res ? "passed" : "failed")
<< " as expected." << std::endl;
}
if (FailOnFail) {
assert(res == fe_it->second);
}
}
LOG(FIT) << "* Testing: SignalDef::isNCINC" << std::endl;
std::map isNCINC_PassExpectations;
isNCINC_PassExpectations[&fe_CC0pi_1] = false; // numu CC0pi
isNCINC_PassExpectations[&fe_CC0pi_2] = false; // numu CC0pi (2p2h)
isNCINC_PassExpectations[&fe_CC0pi_3] = false; // numub CC0pi
isNCINC_PassExpectations[&fe_CC0pi_4] = false; // numu CC0pi (RES)
isNCINC_PassExpectations[&fe_CC1pip_1] = false; // numu CC1pip (2p2h)
isNCINC_PassExpectations[&fe_CC1pim_1] = false; // numu CC1pim
isNCINC_PassExpectations[&fe_CC1pi0_1] = false; // numu CC1pi0
isNCINC_PassExpectations[&fe_CC1pi0_2] = false; // nue CC1pi0
isNCINC_PassExpectations[&fe_CCNpi_1] = false; // numu CC multi pi
isNCINC_PassExpectations[&fe_CCNpi_2] = false; // numu CC multi pi
isNCINC_PassExpectations[&fe_CCNpi_3] = false; // numu CC multi pi
isNCINC_PassExpectations[&fe_CCNpi_4] = false; // numu CC multi pi
isNCINC_PassExpectations[&fe_CCCOH_1] = false; // numu CC COH pi
isNCINC_PassExpectations[&fe_NCel_1] = true; // numu NCEl
isNCINC_PassExpectations[&fe_NCel_2] = false; // nue NCEl
isNCINC_PassExpectations[&fe_NC1pi_1] = true; // numu NC1pi
isNCINC_PassExpectations[&fe_NCNpi_1] = true; // numu NC multi pi
ctr = 0;
for (std::map::iterator
fe_it = isNCINC_PassExpectations.begin();
fe_it != isNCINC_PassExpectations.end(); ++fe_it, ++ctr) {
bool res = SignalDef::isNCINC(fe_it->first, 14);
if (res != fe_it->second) {
ERR(FTL) << "Event: (" << ctr << ")\n"
<< fe_it->first->ToString() << std::endl;
ERR(FTL) << (res ? "passed" : "failed")
<< " SignalDef::isNCINC unexpectedly." << std::endl;
} else {
LOG(SAM) << "Event: (" << ctr << ") " << (res ? "passed" : "failed")
<< " as expected." << std::endl;
}
if (FailOnFail) {
assert(res == fe_it->second);
}
}
LOG(FIT) << "* Testing: SignalDef::isCC0pi" << std::endl;
std::map isCC0pi_PassExpectations;
isCC0pi_PassExpectations[&fe_CC0pi_1] = true; // numu CC0pi
isCC0pi_PassExpectations[&fe_CC0pi_2] = true; // numu CC0pi (2p2h)
isCC0pi_PassExpectations[&fe_CC0pi_3] = false; // numub CC0pi
isCC0pi_PassExpectations[&fe_CC0pi_4] = true; // numu CC0pi (RES)
isCC0pi_PassExpectations[&fe_CC1pip_1] = false; // numu CC1pip (2p2h)
isCC0pi_PassExpectations[&fe_CC1pim_1] = false; // numu CC1pim
isCC0pi_PassExpectations[&fe_CC1pi0_1] = false; // numu CC1pi0
isCC0pi_PassExpectations[&fe_CC1pi0_2] = false; // nue CC1pi0
isCC0pi_PassExpectations[&fe_CCNpi_1] = false; // numu CC multi pi
isCC0pi_PassExpectations[&fe_CCNpi_2] = false; // numu CC multi pi
isCC0pi_PassExpectations[&fe_CCNpi_3] = false; // numu CC multi pi
isCC0pi_PassExpectations[&fe_CCNpi_4] = false; // numu CC multi pi
isCC0pi_PassExpectations[&fe_CCCOH_1] = false; // numu CC COH pi
isCC0pi_PassExpectations[&fe_NCel_1] = false; // numu NCEl
isCC0pi_PassExpectations[&fe_NCel_2] = false; // nue NCEl
isCC0pi_PassExpectations[&fe_NC1pi_1] = false; // numu NC1pi
isCC0pi_PassExpectations[&fe_NCNpi_1] = false; // numu NC multi pi
ctr = 0;
for (std::map::iterator
fe_it = isCC0pi_PassExpectations.begin();
fe_it != isCC0pi_PassExpectations.end(); ++fe_it, ++ctr) {
bool res = SignalDef::isCC0pi(fe_it->first, 14);
if (res != fe_it->second) {
ERR(FTL) << "Event: (" << ctr << ")\n"
<< fe_it->first->ToString() << " " << std::endl;
ERR(FTL) << (res ? "passed" : "failed")
<< " SignalDef::isCC0pi unexpectedly." << std::endl;
} else {
LOG(SAM) << "Event: (" << ctr << ") " << (res ? "passed" : "failed")
<< " as expected." << std::endl;
}
if (FailOnFail) {
assert(res == fe_it->second);
}
}
LOG(FIT) << "* Testing: SignalDef::isCCQELike" << std::endl;
std::map isCCQELike_PassExpectations;
isCCQELike_PassExpectations[&fe_CC0pi_1] = true; // numu CC0pi
isCCQELike_PassExpectations[&fe_CC0pi_2] = true; // numu CC0pi (2p2h)
isCCQELike_PassExpectations[&fe_CC0pi_3] = false; // numub CC0pi
isCCQELike_PassExpectations[&fe_CC0pi_4] = false; // numu CC0pi (RES)
- isCCQELike_PassExpectations[&fe_CC1pip_1] = true; // numu CC1pip (2p2h)
+ isCCQELike_PassExpectations[&fe_CC1pip_1] = true; // numu CC1pip (2p2h)
isCCQELike_PassExpectations[&fe_CC1pim_1] = false; // numu CC1pim
isCCQELike_PassExpectations[&fe_CC1pi0_1] = false; // numu CC1pi0
isCCQELike_PassExpectations[&fe_CC1pi0_2] = false; // nue CC1pi0
isCCQELike_PassExpectations[&fe_CCNpi_1] = false; // numu CC multi pi
isCCQELike_PassExpectations[&fe_CCNpi_2] = false; // numu CC multi pi
isCCQELike_PassExpectations[&fe_CCNpi_3] = false; // numu CC multi pi
isCCQELike_PassExpectations[&fe_CCNpi_4] = false; // numu CC multi pi
isCCQELike_PassExpectations[&fe_CCCOH_1] = false; // numu CC COH pi
isCCQELike_PassExpectations[&fe_NCel_1] = false; // numu NCEl
isCCQELike_PassExpectations[&fe_NCel_2] = false; // nue NCEl
isCCQELike_PassExpectations[&fe_NC1pi_1] = false; // numu NC1pi
isCCQELike_PassExpectations[&fe_NCNpi_1] = false; // numu NC multi pi
ctr = 0;
for (std::map::iterator
fe_it = isCCQELike_PassExpectations.begin();
fe_it != isCCQELike_PassExpectations.end(); ++fe_it, ++ctr) {
bool res = SignalDef::isCCQELike(fe_it->first, 14);
if (res != fe_it->second) {
ERR(FTL) << "Event: (" << ctr << ")\n"
<< fe_it->first->ToString() << std::endl;
ERR(FTL) << (res ? "passed" : "failed")
<< " SignalDef::isCCQELike unexpectedly." << std::endl;
} else {
LOG(SAM) << "Event: (" << ctr << ") " << (res ? "passed" : "failed")
<< " as expected." << std::endl;
}
if (FailOnFail) {
assert(res == fe_it->second);
}
}
LOG(FIT) << "* Testing: SignalDef::isCCQE" << std::endl;
std::map isCCQE_PassExpectations;
isCCQE_PassExpectations[&fe_CC0pi_1] = true; // numu CC0pi
isCCQE_PassExpectations[&fe_CC0pi_2] = false; // numu CC0pi (2p2h)
isCCQE_PassExpectations[&fe_CC0pi_3] = false; // numub CC0pi
isCCQE_PassExpectations[&fe_CC0pi_4] = false; // numu CC0pi (RES)
isCCQE_PassExpectations[&fe_CC1pip_1] = false; // numu CC1pip (2p2h)
isCCQE_PassExpectations[&fe_CC1pim_1] = false; // numu CC1pim
isCCQE_PassExpectations[&fe_CC1pi0_1] = false; // numu CC1pi0
isCCQE_PassExpectations[&fe_CC1pi0_2] = false; // nue CC1pi0
isCCQE_PassExpectations[&fe_CCNpi_1] = false; // numu CC multi pi
isCCQE_PassExpectations[&fe_CCNpi_2] = false; // numu CC multi pi
isCCQE_PassExpectations[&fe_CCNpi_3] = false; // numu CC multi pi
isCCQE_PassExpectations[&fe_CCNpi_4] = false; // numu CC multi pi
isCCQE_PassExpectations[&fe_CCCOH_1] = false; // numu CC COH pi
isCCQE_PassExpectations[&fe_NCel_1] = false; // numu NCEl
isCCQE_PassExpectations[&fe_NCel_2] = false; // nue NCEl
isCCQE_PassExpectations[&fe_NC1pi_1] = false; // numu NC1pi
isCCQE_PassExpectations[&fe_NCNpi_1] = false; // numu NC multi pi
ctr = 0;
for (std::map::iterator
fe_it = isCCQE_PassExpectations.begin();
fe_it != isCCQE_PassExpectations.end(); ++fe_it, ++ctr) {
bool res = SignalDef::isCCQE(fe_it->first, 14);
if (res != fe_it->second) {
ERR(FTL) << "Event: (" << ctr << ")\n"
<< fe_it->first->ToString() << std::endl;
ERR(FTL) << (res ? "passed" : "failed")
<< " SignalDef::isCCQE unexpectedly." << std::endl;
} else {
LOG(SAM) << "Event: (" << ctr << ") " << (res ? "passed" : "failed")
<< " as expected." << std::endl;
}
if (FailOnFail) {
assert(res == fe_it->second);
}
}
LOG(FIT) << "* Testing: SignalDef::isCCCOH" << std::endl;
std::map isCCCOH_PassExpectations;
isCCCOH_PassExpectations[&fe_CC0pi_1] = false; // numu CC0pi
isCCCOH_PassExpectations[&fe_CC0pi_2] = false; // numu CC0pi (2p2h)
isCCCOH_PassExpectations[&fe_CC0pi_3] = false; // numub CC0pi
isCCCOH_PassExpectations[&fe_CC0pi_4] = false; // numu CC0pi (RES)
isCCCOH_PassExpectations[&fe_CC1pip_1] = false; // numu CC1pip (2p2h)
isCCCOH_PassExpectations[&fe_CC1pim_1] = false; // numu CC1pim
isCCCOH_PassExpectations[&fe_CC1pi0_1] = false; // numu CC1pi0
isCCCOH_PassExpectations[&fe_CC1pi0_2] = false; // nue CC1pi0
isCCCOH_PassExpectations[&fe_CCNpi_1] = false; // numu CC multi pi
isCCCOH_PassExpectations[&fe_CCNpi_2] = false; // numu CC multi pi
isCCCOH_PassExpectations[&fe_CCNpi_3] = false; // numu CC multi pi
isCCCOH_PassExpectations[&fe_CCNpi_4] = false; // numu CC multi pi
isCCCOH_PassExpectations[&fe_CCCOH_1] = true; // numu CC COH pi
isCCCOH_PassExpectations[&fe_NCel_1] = false; // numu NCEl
isCCCOH_PassExpectations[&fe_NCel_2] = false; // nue NCEl
isCCCOH_PassExpectations[&fe_NC1pi_1] = false; // numu NC1pi
isCCCOH_PassExpectations[&fe_NCNpi_1] = false; // numu NC multi pi
ctr = 0;
for (std::map::iterator
fe_it = isCCCOH_PassExpectations.begin();
fe_it != isCCCOH_PassExpectations.end(); ++fe_it, ++ctr) {
bool res = SignalDef::isCCCOH(fe_it->first, 14, 211);
if (res != fe_it->second) {
ERR(FTL) << "Event: (" << ctr << ")\n"
<< fe_it->first->ToString() << std::endl;
ERR(FTL) << (res ? "passed" : "failed")
<< " SignalDef::isCCCOH unexpectedly." << std::endl;
} else {
LOG(SAM) << "Event: (" << ctr << ") " << (res ? "passed" : "failed")
<< " as expected." << std::endl;
}
if (FailOnFail) {
assert(res == fe_it->second);
}
}
LOG(FIT) << "* Testing: SignalDef::isCC1pi" << std::endl;
std::map isCC1pi_PassExpectations;
isCC1pi_PassExpectations[&fe_CC0pi_1] = false; // numu CC0pi
isCC1pi_PassExpectations[&fe_CC0pi_2] = false; // numu CC0pi (2p2h)
isCC1pi_PassExpectations[&fe_CC0pi_3] = false; // numub CC0pi
isCC1pi_PassExpectations[&fe_CC0pi_4] = false; // numu CC0pi (RES)
isCC1pi_PassExpectations[&fe_CC1pip_1] = true; // numu CC1pip (2p2h)
isCC1pi_PassExpectations[&fe_CC1pim_1] = false; // numu CC1pim
isCC1pi_PassExpectations[&fe_CC1pi0_1] = false; // numu CC1pi0
isCC1pi_PassExpectations[&fe_CC1pi0_2] = false; // nue CC1pi0
isCC1pi_PassExpectations[&fe_CCNpi_1] = false; // numu CC multi pi
isCC1pi_PassExpectations[&fe_CCNpi_2] = false; // numu CC multi pi
isCC1pi_PassExpectations[&fe_CCNpi_3] = false; // numu CC multi pi
isCC1pi_PassExpectations[&fe_CCNpi_4] = false; // numu CC multi pi
- isCC1pi_PassExpectations[&fe_CCCOH_1] = true; // numu CC COH pi
+ isCC1pi_PassExpectations[&fe_CCCOH_1] = true; // numu CC COH pi
isCC1pi_PassExpectations[&fe_NCel_1] = false; // numu NCEl
isCC1pi_PassExpectations[&fe_NCel_2] = false; // nue NCEl
isCC1pi_PassExpectations[&fe_NC1pi_1] = false; // numu NC1pi
isCC1pi_PassExpectations[&fe_NCNpi_1] = false; // numu NC multi pi
ctr = 0;
for (std::map::iterator
fe_it = isCC1pi_PassExpectations.begin();
fe_it != isCC1pi_PassExpectations.end(); ++fe_it, ++ctr) {
bool res = SignalDef::isCC1pi(fe_it->first, 14, 211);
if (res != fe_it->second) {
ERR(FTL) << "Event: (" << ctr << ")\n"
<< fe_it->first->ToString() << std::endl;
ERR(FTL) << (res ? "passed" : "failed")
<< " SignalDef::isCC1pi unexpectedly." << std::endl;
} else {
LOG(SAM) << "Event: (" << ctr << ") " << (res ? "passed" : "failed")
<< " as expected." << std::endl;
}
if (FailOnFail) {
assert(res == fe_it->second);
}
}
LOG(FIT) << "* Testing: SignalDef::isNC1pi" << std::endl;
std::map isNC1pi_PassExpectations;
isNC1pi_PassExpectations[&fe_CC0pi_1] = false; // numu CC0pi
isNC1pi_PassExpectations[&fe_CC0pi_2] = false; // numu CC0pi (2p2h)
isNC1pi_PassExpectations[&fe_CC0pi_3] = false; // numub CC0pi
isCCINC_PassExpectations[&fe_CC0pi_4] = false; // numu CC0pi (RES)
isNC1pi_PassExpectations[&fe_CC1pip_1] = false; // numu CC1pip (2p2h)
isNC1pi_PassExpectations[&fe_CC1pim_1] = false; // numu CC1pim
isNC1pi_PassExpectations[&fe_CC1pi0_1] = false; // numu CC1pi0
isNC1pi_PassExpectations[&fe_CC1pi0_2] = false; // nue CC1pi0
isNC1pi_PassExpectations[&fe_CCNpi_1] = false; // numu CC multi pi
isNC1pi_PassExpectations[&fe_CCNpi_2] = false; // numu CC multi pi
isNC1pi_PassExpectations[&fe_CCNpi_3] = false; // numu CC multi pi
isNC1pi_PassExpectations[&fe_CCNpi_4] = false; // numu CC multi pi
isCCINC_PassExpectations[&fe_CCCOH_1] = false; // numu CC COH pi
isNC1pi_PassExpectations[&fe_NCel_1] = false; // numu NCEl
isNC1pi_PassExpectations[&fe_NCel_2] = false; // nue NCEl
isCCINC_PassExpectations[&fe_NC1pi_1] = true; // numu NC1pi
isCCINC_PassExpectations[&fe_NCNpi_1] = false; // numu NC multi pi
ctr = 0;
for (std::map::iterator
fe_it = isNC1pi_PassExpectations.begin();
fe_it != isNC1pi_PassExpectations.end(); ++fe_it, ++ctr) {
bool res = SignalDef::isNC1pi(fe_it->first, 14, -211);
if (res != fe_it->second) {
ERR(FTL) << "Event: (" << ctr << ")\n"
<< fe_it->first->ToString() << std::endl;
ERR(FTL) << (res ? "passed" : "failed")
<< " SignalDef::isNC1pi unexpectedly." << std::endl;
} else {
LOG(SAM) << "Event: (" << ctr << ") " << (res ? "passed" : "failed")
<< " as expected." << std::endl;
}
if (FailOnFail) {
assert(res == fe_it->second);
}
}
// SignalDef::isCCWithFS(&fe,14);
}