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); }