No OneTemporary
Actions

Size

37 KB

Subscribers

None

View Options

	diff --git a/src/InputHandler/NEUTInputHandler.cxx b/src/InputHandler/NEUTInputHandler.cxx
	index 7df81d1..9154fe6 100644
	--- a/src/InputHandler/NEUTInputHandler.cxx
	+++ b/src/InputHandler/NEUTInputHandler.cxx
	@@ -1,503 +1,507 @@
	#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<std::string> inputs = InputUtils::ParseInputFileList(rawinputs);
	for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) {
	// Open File for histogram access
	TFile* inp_file = new TFile(inputs[inp_it].c_str(), "READ");
	if (!inp_file or inp_file->IsZombie()) {
	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;
	}

	// From NEUT neutclass/neutpart.h
	// Bool_t fIsAlive; // Particle should be tracked or not
	// ( in the detector simulator )
	//
	// Int_t fStatus; // Status flag of this particle
	// -2: Non existing particle
	// -1: Initial state particle
	// 0: Normal
	// 1: Decayed to the other particle
	// 2: Escaped from the detector
	// 3: Absorped
	// 4: Charge exchanged
	// 5: Pauli blocked
	// 6: N/A
	// 7: Produced child particles
	// 8: Inelastically scattered
	//
	int NEUTInputHandler::GetNeutParticleStatus(NeutPart* part) {
	// State
	int state = kUndefinedState;

	// Remove Pauli blocked events, probably just single pion events
	if (part->fStatus == 5) {
	state = kFSIState;

	// fStatus == -1 means initial state
	} else if (part->fIsAlive == false && part->fStatus == -1) {
	state = kInitialState;

	// NEUT has a bit of a strange convention for fIsAlive and fStatus
	// combinations
	// for NC and neutrino particle isAlive true/false and status 2 means
	// final state particle
	// for other particles in NC status 2 means it's an FSI particle
	// for CC it means it was an FSI particle
	} else if (part->fStatus == 2) {
	// NC case is a little strange... The outgoing neutrino might be alive or
	// not alive. Remaining particles with status 2 are FSI particles that
	// reinteracted
	if (abs(fNeutVect->Mode) > 30 &&
	(abs(part->fPID) == 16 \|\| abs(part->fPID) == 14 \|\| abs(part->fPID) == 12)) {
	state = kFinalState;
	// The usual CC case
	} else if (part->fIsAlive == true) {
	state = kFSIState;
	}

	} else if (part->fIsAlive == true && part->fStatus == 2 &&
	(abs(part->fPID) == 16 \|\| abs(part->fPID) == 14 \|\| abs(part->fPID) == 12)) {
	state = kFinalState;

	} else if (part->fIsAlive == true && part->fStatus == 0) {
	state = kFinalState;

	} else if (!part->fIsAlive && (part->fStatus == 1 \|\| part->fStatus == 3 \|\| part->fStatus == 4 \|\| part->fStatus == 7 \|\| part->fStatus == 8)) {
	state = kFSIState;

	// There's one hyper weird case where fStatus = -3. This apparently corresponds to a nucleon being ejected via pion FSI when there is "data available"
	} else if (!part->fIsAlive && (part->fStatus == -3)) {
	state = kUndefinedState;
	// NC neutrino outgoing
	} else if (!part->fIsAlive && part->fStatus == 0 && (abs(part->fPID) == 16 \|\| abs(part->fPID) == 14 \|\| abs(part->fPID) == 12)) {
	state = kFinalState;

	// Warn if we still find alive particles without classifying them
	} else if (part->fIsAlive == true) {
	ERR(WRN) << "Undefined NEUT state "
	<< " Alive: " << part->fIsAlive << " Status: " << part->fStatus
	<< " PDG: " << part->fPID << std::endl;
	throw;
	// Warn if we find dead particles that we haven't classified
	} else {
	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->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;

	+ // Is the paricle associated with the primary vertex?
	+ fNUISANCEEvent->fPrimaryVertex[curpart] = fNeutVect->ParentIdx(i) == 0;
	+ std::cout << fNeutVect->ParentIdx(i) << std::endl;
	+
	// 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();

	#ifdef NEUT_COMMON_QEAV
	nemdls_.mdlqeaf = nvect->QEAVForm;
	#else
	nemdls_.mdlqeaf = nvect->QEVForm;
	#endif
	nemdls_.mdlqe = nvect->QEModel;
	nemdls_.mdlspi = nvect->SPIModel;
	nemdls_.mdldis = nvect->DISModel;
	nemdls_.mdlcoh = nvect->COHModel;
	neutcoh_.necohepi = nvect->COHModel;

	nemdls_.xmaqe = nvect->QEMA;
	nemdls_.xmvqe = nvect->QEMV;
	nemdls_.kapp = nvect->KAPPA;

	// nemdls_.sccfv = SCCFVdef;
	// nemdls_.sccfa = SCCFAdef;
	// nemdls_.fpqe = FPQEdef;

	nemdls_.xmaspi = nvect->SPIMA;
	nemdls_.xmvspi = nvect->SPIMV;
	nemdls_.xmares = nvect->RESMA;
	nemdls_.xmvres = nvect->RESMV;

	neut1pi_.xmanffres = nvect->SPIMA;
	neut1pi_.xmvnffres = nvect->SPIMV;
	neut1pi_.xmarsres = nvect->RESMA;
	neut1pi_.xmvrsres = nvect->RESMV;
	neut1pi_.neiff = nvect->SPIForm;
	neut1pi_.nenrtype = nvect->SPINRType;
	neut1pi_.rneca5i = nvect->SPICA5I;
	neut1pi_.rnebgscl = nvect->SPIBGScale;

	nemdls_.xmacoh = nvect->COHMA;
	nemdls_.rad0nu = nvect->COHR0;
	// nemdls_.fa1coh = nvect->COHA1err;
	// nemdls_.fb1coh = nvect->COHb1err;

	// neutdis_.nepdf = NEPDFdef;
	// neutdis_.nebodek = NEBODEKdef;

	neutcard_.nefrmflg = nvect->FrmFlg;
	neutcard_.nepauflg = nvect->PauFlg;
	neutcard_.nenefo16 = nvect->NefO16;
	neutcard_.nemodflg = nvect->ModFlg;
	// neutcard_.nenefmodl = 1;
	// neutcard_.nenefmodh = 1;
	// neutcard_.nenefkinh = 1;
	// neutpiabs_.neabspiemit = 1;

	nenupr_.iformlen = nvect->FormLen;

	neutpiless_.ipilessdcy = nvect->IPilessDcy;
	neutpiless_.rpilessdcy = nvect->RPilessDcy;

	neutpiless_.ipilessdcy = nvect->IPilessDcy;
	neutpiless_.rpilessdcy = nvect->RPilessDcy;

	neffpr_.fefqe = nvect->NuceffFactorPIQE;
	neffpr_.fefqeh = nvect->NuceffFactorPIQEH;
	neffpr_.fefinel = nvect->NuceffFactorPIInel;
	neffpr_.fefabs = nvect->NuceffFactorPIAbs;
	neffpr_.fefcx = nvect->NuceffFactorPICX;
	neffpr_.fefcxh = nvect->NuceffFactorPICXH;

	neffpr_.fefcoh = nvect->NuceffFactorPICoh;
	neffpr_.fefqehf = nvect->NuceffFactorPIQEHKin;
	neffpr_.fefcxhf = nvect->NuceffFactorPICXKin;
	neffpr_.fefcohf = nvect->NuceffFactorPIQELKin;

	for (int i = 0; i < nework_.numne; i++) {
	nework_.ipne[i] = nvect->PartInfo(i)->fPID;
	nework_.pne[i][0] =
	(float)nvect->PartInfo(i)->fP.X() / 1000; // VC(NE)WORK in M(G)eV
	nework_.pne[i][1] =
	(float)nvect->PartInfo(i)->fP.Y() / 1000; // VC(NE)WORK in M(G)eV
	nework_.pne[i][2] =
	(float)nvect->PartInfo(i)->fP.Z() / 1000; // VC(NE)WORK in M(G)eV
	}
	// fsihist.h

	// neutroot fills a dummy object for events with no FSI to prevent memory leak
	// when
	// reading the TTree, so check for it here

	if ((int)nvect->NfsiVert() ==
	1) { // An event with FSI must have at least two vertices
	// if (nvect->NfsiPart()!=1 \|\| nvect->Fsiprob!=-1)
	// ERR(WRN) << "T2KNeutUtils::fill_neut_commons(TTree) NfsiPart!=1 or
	// Fsiprob!=-1 when NfsiVert==1" << std::endl;

	fsihist_.nvert = 0;
	fsihist_.nvcvert = 0;
	fsihist_.fsiprob = 1;
	} else { // Real FSI event
	fsihist_.nvert = (int)nvect->NfsiVert();
	for (int ivert = 0; ivert < fsihist_.nvert; ivert++) {
	fsihist_.iflgvert[ivert] = nvect->FsiVertInfo(ivert)->fVertID;
	fsihist_.posvert[ivert][0] = (float)nvect->FsiVertInfo(ivert)->fPos.X();
	fsihist_.posvert[ivert][1] = (float)nvect->FsiVertInfo(ivert)->fPos.Y();
	fsihist_.posvert[ivert][2] = (float)nvect->FsiVertInfo(ivert)->fPos.Z();
	}

	fsihist_.nvcvert = nvect->NfsiPart();
	for (int ip = 0; ip < fsihist_.nvcvert; ip++) {
	fsihist_.abspvert[ip] = (float)nvect->FsiPartInfo(ip)->fMomLab;
	fsihist_.abstpvert[ip] = (float)nvect->FsiPartInfo(ip)->fMomNuc;
	fsihist_.ipvert[ip] = nvect->FsiPartInfo(ip)->fPID;
	fsihist_.iverti[ip] = nvect->FsiPartInfo(ip)->fVertStart;
	fsihist_.ivertf[ip] = nvect->FsiPartInfo(ip)->fVertEnd;
	fsihist_.dirvert[ip][0] = (float)nvect->FsiPartInfo(ip)->fDir.X();
	fsihist_.dirvert[ip][1] = (float)nvect->FsiPartInfo(ip)->fDir.Y();
	fsihist_.dirvert[ip][2] = (float)nvect->FsiPartInfo(ip)->fDir.Z();
	}
	fsihist_.fsiprob = nvect->Fsiprob;
	}

	neutcrscom_.crsx = nvect->Crsx;
	neutcrscom_.crsy = nvect->Crsy;
	neutcrscom_.crsz = nvect->Crsz;
	neutcrscom_.crsphi = nvect->Crsphi;
	neutcrscom_.crsq2 = nvect->Crsq2;

	neuttarget_.numbndn = nvect->TargetA - nvect->TargetZ;
	neuttarget_.numbndp = nvect->TargetZ;
	neuttarget_.numfrep = nvect->TargetH;
	neuttarget_.numatom = nvect->TargetA;
	posinnuc_.ibound = nvect->Ibound;

	// put empty nucleon FSI history (since it is not saved in the NeutVect
	// format)
	// Comment out as NEUT does not have the necessary proton FSI information yet
	// nucleonfsihist_.nfnvert = 0;
	// nucleonfsihist_.nfnstep = 0;
	}

	#endif
	diff --git a/src/InputHandler/NuWroInputHandler.cxx b/src/InputHandler/NuWroInputHandler.cxx
	index 83dc5f2..1853d95 100644
	--- a/src/InputHandler/NuWroInputHandler.cxx
	+++ b/src/InputHandler/NuWroInputHandler.cxx
	@@ -1,497 +1,507 @@
	#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<std::string> inputs = InputUtils::ParseInputFileList(rawinputs);
	for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) {
	// Open File for histogram access
	TFile* inp_file = new TFile(inputs[inp_it].c_str(), "READ");
	if (!inp_file or inp_file->IsZombie()) {
	ERR(FTL) << "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().fType = kNUWRO;
	rwEvs.back().Mode = fNUISANCEEvent->Mode;
	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;
	}

	fNUISANCEEvent->fNuwroSRWEvent = SRW::SRWEvent(*fNuWroEvent);
	fNUISANCEEvent->fNuwroParams = &fNuwroParams;
	fNUISANCEEvent->probe_E = fNUISANCEEvent->fNuwroSRWEvent.NeutrinoEnergy;
	fNUISANCEEvent->probe_pdg = fNUISANCEEvent->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->Mode = ConvertNuwroMode(fNuWroEvent);

	if (abs(evt->Mode) > 60) {
	evt->Mode = 0;
	}

	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->fTargetPDG = TargetUtils::GetTargetPDGFromZA(evt->fTargetZ, evt->fTargetA);
	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<particle>::iterator p_iter;

	- // Initial State
	- for (p_iter = fNuWroEvent->in.begin(); p_iter != fNuWroEvent->in.end();
	- p_iter++) {
	- AddNuWroParticle(fNUISANCEEvent, (*p_iter), kInitialState);
	+ // Get the Initial State
	+ for (p_iter = fNuWroEvent->in.begin(); p_iter != fNuWroEvent->in.end(); p_iter++) {
	+ AddNuWroParticle(fNUISANCEEvent, (*p_iter), kInitialState, true);
	}

	- for (p_iter = fNuWroEvent->all.begin(); p_iter != fNuWroEvent->all.end(); p_iter++) {
	- std::cout << (*p_iter)->pdg << std::endl;
	- std::cout << (*p_iter)->ks << std::endl;
	- std::cout << (*p_iter)->orgig << std::endl;
	- std::cout << (*p_iter)->id << std::endl;
	- std::cout << (*p_iter)->mother << std::endl;
	- std::cout << (*p_iter)->endproc << std::endl;
	- std::cout << (*p_iter)->his_fermi << std::endl;
	- std::cout << (*p_iter)->primary << std::endl;
	- std::cout << (*p_iter)->E() << std::endl;
	- }
	-
	- // NuWro saves the incoming, pre-FSI and post-FSI particles
	- // So we fill the incoming and outgoing particles, but no the FSI particles
	- // This is because sometimes the pre-FSI particles are the same as the post-FSI particles, so we would double count our particle stack
	- /*
	- for (p_iter = fNuWroEvent->out.begin(); p_iter != fNuWroEvent->out.end();
	- p_iter++) {
	- AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFSIState);
	+ // Try to find the FSI state particles
	+ // Loop over the primary vertex particles
	+ // If they match the post-FSI they haven't undergone FSI.
	+ // If they don't match post-FSI they have undergone FSI.
	+ for (p_iter = fNuWroEvent->out.begin(); p_iter != fNuWroEvent->out.end(); p_iter++) {
	+ // Get the particle
	+ particle p = (*p_iter);
	+ // Check against all the post particles, match them
	+ std::vector<particle>::iterator p2_iter;
	+ bool match = false;
	+ for (p2_iter = fNuWroEvent->post.begin(); p2_iter != fNuWroEvent->post.end(); p2_iter++) {
	+ particle p2 = (*p2_iter);
	+ // Check energy and pdg
	+ // A very small cascade which changes the energy by 1E-5 MeV should be matched
	+ match = (fabs(p2.E()-p.E()) < 1E-5 && p2.pdg == p.pdg);
	+ // If we match p to p2 break the loop
	+ if (match) break;
	+ }
	+ // If we've looped through the whole particle stack of post-FSI and haven't found a match it's a primary particle that has been FSIed
	+ if (!match) AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFSIState, true);
	}
	- */

	- // Final State
	- for (p_iter = fNuWroEvent->post.begin(); p_iter != fNuWroEvent->post.end();
	- p_iter++) {
	- AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFinalState);
	+ // Loop over the final state particles
	+ for (p_iter = fNuWroEvent->post.begin(); p_iter != fNuWroEvent->post.end(); p_iter++) {
	+ particle p = (*p_iter);
	+ // To find if it's primary or not we have to loop through the primary ones and match, just like above
	+ bool match = false;
	+ std::vector<particle>::iterator p2_iter;
	+ for (p2_iter = fNuWroEvent->out.begin(); p2_iter != fNuWroEvent->out.end(); p2_iter++) {
	+ particle p2 = (*p2_iter);
	+ match = (fabs(p2.E()-p.E()) < 1E-5 && p2.pdg == p.pdg);
	+ if (match) break;
	+ }
	+ AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFinalState, match);
	}

	// 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) {
	+ int state, bool primary = false) {
	// Add Mom
	evt->fParticleMom[evt->fNParticles][0] = static_cast<vect&>(p).x;
	evt->fParticleMom[evt->fNParticles][1] = static_cast<vect&>(p).y;
	evt->fParticleMom[evt->fNParticles][2] = static_cast<vect&>(p).z;
	evt->fParticleMom[evt->fNParticles][3] = static_cast<vect&>(p).t;

	+ // For NuWro a particle that we've given a FSI state is a pre-FSI particle
	+ // An initial state particle is also a primary vertex praticle
	+ evt->fPrimaryVertex[evt->fNParticles] = primary;
	+
	// 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/NuWroInputHandler.h b/src/InputHandler/NuWroInputHandler.h
	index 5620a0f..43d7443 100644
	--- a/src/InputHandler/NuWroInputHandler.h
	+++ b/src/InputHandler/NuWroInputHandler.h
	@@ -1,95 +1,95 @@
	#ifndef NuWroINPUTHANDLER_H
	#define NuWroINPUTHANDLER_H
	#ifdef __NUWRO_ENABLED__
	#include "GeneratorUtils.h"
	#include "InputHandler.h"
	#include "PlotUtils.h"

	/// NuWro Generator Container to save extra particle status codes.
	class NuWroGeneratorInfo : public GeneratorInfoBase {
	public:
	NuWroGeneratorInfo(){};
	virtual ~NuWroGeneratorInfo();

	/// Assigns information to branches
	void AddBranchesToTree(TTree* tn);

	/// Setup reading information from branches
	void SetBranchesFromTree(TTree* tn);

	/// Allocate any dynamic arrays for a new particle stack size
	void AllocateParticleStack(int stacksize);

	/// Clear any dynamic arrays
	void DeallocateParticleStack();

	/// Read extra genie information from the event
	void FillGeneratorInfo(event* e);

	/// Reset extra information to default/empty values
	void Reset();

	int kMaxParticles; ///< Number of particles in stack
	int* fNuWroParticlePDGs; ///< NuWro Particle PDGs (example)
	};

	/// Main NuWro Input Reader. Requires events have flux and xsec TH1Ds saved into
	/// them.
	class NuWroInputHandler : public InputHandlerBase {
	#ifdef __USE_NUWRO_SRW_EVENTS__
	params fNuwroParams;
	std::vector<BaseFitEvt> rwEvs;
	#endif

	public:
	/// Constructor. Can handle single and joint inputs.
	NuWroInputHandler(std::string const& handle, std::string const& rawinputs);
	~NuWroInputHandler();

	/// Create a TTree Cache to speed up file read
	void CreateCache();

	/// Remove TTree Cache to save memory
	void RemoveCache();

	/// Returns filled NUISANCEEvent for given entry.
	FitEvent* GetNuisanceEvent(const UInt_t entry,
	const bool lightweight = false);

	#ifdef __USE_NUWRO_SRW_EVENTS__
	// Returns filled BaseFitEvent for a given entry;
	BaseFitEvt* GetBaseEvent(const UInt_t entry) {
	if (rwEvs.size() <= entry) {
	THROW("Tried to get cached BaseFitEv[" << entry << "], but only have "
	<< rwEvs.size()
	<< " in the cache.");
	}
	return &rwEvs[entry];
	}
	#endif

	/// Fills fNUISANCEEvent from fNuWroEvent
	void CalcNUISANCEKinematics();

	/// (LEGACY) Automatically creates nuwro flux/event histograms that
	/// nuisance needs to normalise events.
	void ProcessNuWroInputFlux(const std::string file);

	/// Calculates a True Interaction code for NuWro events
	int ConvertNuwroMode(event* e);

	/// Adds a new particle to NUISANCE stack for given NuWro particle
	- void AddNuWroParticle(FitEvent* evt, particle& p, int state);
	+ void AddNuWroParticle(FitEvent* evt, particle& p, int state, bool primary);

	event* fNuWroEvent; ///< Pointer to NuWro Format Events

	/// Print Event Information
	void Print();

	TChain* fNuWroTree; ///< TTree for reading NuWro event vectors
	bool fSaveExtra; ///< Save Extra NuWro info into Nuisance Event
	NuWroGeneratorInfo* fNuWroInfo; ///< Extra NuWro Generator Info
	};
	/! @} /
	#endif
	#endif

File Metadata

Mime Type: text/x-diff
Expires: Wed, May 14, 11:20 AM (17 h, 14 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 5093849
Default Alt Text: (37 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions