No OneTemporary
Actions

Size

17 KB

Subscribers

None

View Options

	diff --git a/src/InputHandler/NEUTInputHandler.cxx b/src/InputHandler/NEUTInputHandler.cxx
	index 189026d..8e3d5a2 100644
	--- a/src/InputHandler/NEUTInputHandler.cxx
	+++ b/src/InputHandler/NEUTInputHandler.cxx
	@@ -1,459 +1,499 @@
	#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
	- if (part->fIsAlive == false && part->fStatus == -1) {
	+ } 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) == 14 \|\| abs(part->fPID) == 12)) {
	+ (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) == 14 \|\| abs(part->fPID) == 12)) {
	- state = kFinalState;
	+ (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;
	+ << " 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;

	// 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);
	+ 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
	+ (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
	+ (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
	+ (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

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 4:53 PM (1 d, 12 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805256
Default Alt Text: (17 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions