No OneTemporary
Actions

Size

31 KB

Subscribers

None

View Options

	diff --git a/app/PrepareGENIE.cxx b/app/PrepareGENIE.cxx
	index 0820104..3839119 100644
	--- a/app/PrepareGENIE.cxx
	+++ b/app/PrepareGENIE.cxx
	@@ -1,611 +1,614 @@
	#include <stdio.h>
	#include <stdlib.h>
	#include "FitLogger.h"
	#include "PlotUtils.h"
	#include "TFile.h"
	#include "TH1D.h"
	#include "TTree.h"

	#ifdef __GENIE_ENABLED__
	#include "Conventions/Units.h"
	#include "GHEP/GHepParticle.h"
	#include "PDG/PDGUtils.h"
	#endif

	std::string gInputFiles = "";
	std::string gOutputFile = "";
	std::string gFluxFile = "";
	std::string gTarget = "";
	double MonoEnergy;
	bool IsMonoE = false;

	void PrintOptions();
	void ParseOptions(int argc, char* argv[]);
	void RunGENIEPrepareMono(std::string input, std::string target,
	std::string output);
	void RunGENIEPrepare(std::string input, std::string flux, std::string target,
	std::string output);

	int main(int argc, char* argv[]) {
	ParseOptions(argc, argv);
	if (IsMonoE) {
	RunGENIEPrepareMono(gInputFiles, gTarget, gOutputFile);
	} else {
	RunGENIEPrepare(gInputFiles, gFluxFile, gTarget, gOutputFile);
	}
	}

	void RunGENIEPrepareMono(std::string input, std::string target,
	std::string output) {
	+
	+ std::cout << "Running in mono" << std::endl;
	// Setup TTree
	TChain* tn = new TChain("gtree");
	tn->AddFile(input.c_str());

	int nevt = tn->GetEntries();
	NtpMCEventRecord* genientpl = NULL;
	tn->SetBranchAddress("gmcrec", &genientpl);

	TH1D* fluxhist = new TH1D("flux", "flux", 1000, 0, 10);
	fluxhist->Fill(MonoEnergy);
	fluxhist->Scale(1, "width");

	// Make Event Hist
	TH1D* eventhist = (TH1D*)fluxhist->Clone();
	eventhist->Reset();

	TH1D* xsechist = (TH1D*)eventhist->Clone();

	// Create maps
	std::map<std::string, TH1D*> modexsec;
	std::map<std::string, TH1D*> modecount;
	std::vector<std::string> genieids;
	std::vector<std::string> targetids;
	std::vector<std::string> interids;

	// Loop over all events
	for (int i = 0; i < nevt; i++) {
	tn->GetEntry(i);

	StopTalking();
	EventRecord& event = *(genientpl->event);
	GHepParticle* neu = event.Probe();
	StartTalking();

	// Get XSec From Spline
	GHepRecord genie_record = static_cast<GHepRecord>(event);
	double xsec = (genie_record.XSec() / (1E-38 * genie::units::cm2));

	// Parse Interaction String
	std::string mode = genie_record.Summary()->AsString();
	std::vector<std::string> modevec = GeneralUtils::ParseToStr(mode, ";");
	std::string targ = (modevec[0] + ";" + modevec[1]);
	std::string inter = mode;

	// Fill lists of Unique IDS
	if (std::find(targetids.begin(), targetids.end(), targ) ==
	targetids.end()) {
	targetids.push_back(targ);
	}

	if (std::find(interids.begin(), interids.end(), inter) == interids.end()) {
	interids.push_back(inter);
	}

	// Create entries Mode Maps
	if (modexsec.find(mode) == modexsec.end()) {
	genieids.push_back(mode);

	modexsec[mode] = (TH1D*)xsechist->Clone();
	modecount[mode] = (TH1D*)xsechist->Clone();
	}

	// Fill XSec Histograms
	modexsec[mode]->Fill(neu->E(), xsec);
	modecount[mode]->Fill(neu->E());

	// Fill total event hist
	eventhist->Fill(neu->E());

	// Clear Event
	genientpl->Clear();

	if (i % (nevt / 20) == 0) {
	LOG(FIT) << "Processed " << i << "/" << nevt << " GENIE events."
	<< std::endl;
	}
	}
	LOG(FIT) << "Processed all events" << std::endl;

	TFile* outputfile = new TFile(input.c_str(), "UPDATE");
	outputfile->cd();

	- LOG(FIT) << "Getting splines " << std::endl;
	+ LOG(FIT) << "Getting splines in mono" << std::endl;

	// Save each of the reconstructed splines to file
	std::map<std::string, TH1D*> modeavg;

	TDirectory* inddir = (TDirectory*)outputfile->Get("IndividualGENIESplines");
	if (!inddir)
	inddir = (TDirectory*)outputfile->mkdir("IndividualGENIESplines");
	inddir->cd();

	// Loop over GENIE ID's and get MEC count
	int MECcount = 0;
	bool MECcorrect = FitPar::Config().GetParB("CorrectGENIEMECNorm");
	for (UInt_t i = 0; i < genieids.size(); i++) {
	if (genieids[i].find("MEC") != std::string::npos) {
	MECcount++;
	}
	}
	LOG(FIT) << "Found " << MECcount << " repeated MEC instances." << std::endl;

	for (UInt_t i = 0; i < genieids.size(); i++) {
	std::string mode = genieids[i];

	modexsec[mode]->Write((mode + "_summed_xsec").c_str(), TObject::kOverwrite);
	modecount[mode]->Write((mode + "_summed_evt").c_str(), TObject::kOverwrite);

	// Form extra avg xsec map -> Reconstructed spline
	modeavg[mode] = (TH1D*)modexsec[mode]->Clone();
	modeavg[mode]->Divide(modecount[mode]);

	if (MECcorrect && (mode.find("MEC") != std::string::npos)) {
	modeavg[mode]->Scale(1.0 / double(MECcount));
	}

	modeavg[mode]->Write((mode + "_rec_spline").c_str(), TObject::kOverwrite);
	}

	TDirectory* targdir = (TDirectory*)outputfile->Get("TargetGENIESplines");
	if (!targdir) targdir = (TDirectory*)outputfile->mkdir("TargetGENIESplines");
	targdir->cd();

	LOG(FIT) << "Getting Target Splines" << std::endl;
	// For each target save a total spline
	std::map<std::string, TH1D*> targetsplines;

	for (uint i = 0; i < targetids.size(); i++) {
	LOG(FIT) << "Getting target " << i << std::endl;
	std::string targ = targetids[i];
	targetsplines[targ] = (TH1D*)xsechist->Clone();
	LOG(FIT) << "Created target spline for " << targ << std::endl;

	for (uint j = 0; j < genieids.size(); j++) {
	std::string mode = genieids[j];

	if (mode.find(targ) != std::string::npos) {
	LOG(FIT) << "Mode " << mode << " contains " << targ << " target!"
	<< std::endl;
	targetsplines[targ]->Add(modeavg[mode]);
	LOG(FIT) << "Finished with Mode " << mode << " "
	<< modeavg[mode]->Integral() << std::endl;
	}
	}

	LOG(FIT) << "Saving target spline:" << targ << std::endl;
	targetsplines[targ]->Write(("Total" + targ).c_str(), TObject::kOverwrite);
	}

	LOG(FIT) << "Getting total splines" << std::endl;
	// Now we have each of the targets we need to create a total cross-section.
	int totalnucl = 0;
	std::vector<std::string> targprs = GeneralUtils::ParseToStr(target, ",");
	TH1D* totalxsec = (TH1D*)xsechist->Clone();

	for (uint i = 0; i < targprs.size(); i++) {
	std::string targpdg = targprs[i];

	for (std::map<std::string, TH1D*>::iterator iter = targetsplines.begin();
	iter != targetsplines.end(); iter++) {
	std::string targstr = iter->first;
	TH1D* xsec = iter->second;

	if (targstr.find(targpdg) != std::string::npos) {
	LOG(FIT) << "Adding target spline " << targstr
	<< " Integral = " << xsec->Integral("width") << std::endl;
	totalxsec->Add(xsec);

	int nucl = atoi(targpdg.c_str());
	totalnucl += int((nucl % 10000) / 10);
	}
	}
	}

	outputfile->cd();
	totalxsec->Write("nuisance_Xsec", TObject::kOverwrite);
	eventhist = (TH1D*)totalxsec->Clone();
	eventhist->Multiply(fluxhist);

	eventhist->Write("nuisance_events", TObject::kOverwrite);
	fluxhist->Write("nuisance_flux", TObject::kOverwrite);

	LOG(FIT) << "Inclusive XSec Per Nucleon = "
	<< eventhist->Integral("width") * 1E-38 / fluxhist->Integral("width")
	<< std::endl;
	std::cout << "XSec Hist Integral = " << xsechist->Integral("width")
	<< std::endl;

	return;
	}

	void RunGENIEPrepare(std::string input, std::string flux, std::string target,
	std::string output) {
	LOG(FIT) << "Running GENIE Prepare" << std::endl;
	+ std::cout << "Running in prepare" << std::endl;

	// Get Flux Hist
	std::vector<std::string> fluxvect = GeneralUtils::ParseToStr(flux, ",");
	TH1* fluxhist = NULL;
	if (fluxvect.size() == 3) {
	double from = GeneralUtils::StrToDbl(fluxvect[0]);
	double to = GeneralUtils::StrToDbl(fluxvect[1]);
	double step = GeneralUtils::StrToDbl(fluxvect[2]);

	int nstep = ceil((to - from) / step);
	to = from + step * nstep;

	QLOG(FIT, "Generating flat flux histogram from "
	<< from << " to " << to << " with bins " << step
	<< " wide (NBins = " << nstep << ").");

	fluxhist =
	new TH1D("spectrum", ";E_{#nu} (GeV);Count (A.U.)", nstep, from, to);

	for (Int_t bi_it = 1; bi_it < fluxhist->GetXaxis()->GetNbins(); ++bi_it) {
	fluxhist->SetBinContent(bi_it, 1.0 / double(step * nstep));
	}
	fluxhist->SetDirectory(0);
	} else if (fluxvect.size() == 2) {
	TFile* fluxfile = new TFile(fluxvect[0].c_str(), "READ");
	if (!fluxfile->IsZombie()) {
	fluxhist = dynamic_cast<TH1*>(fluxfile->Get(fluxvect[1].c_str()));
	if (!fluxhist) {
	ERR(FTL) << "Couldn't find histogram named: \"" << fluxvect[1]
	<< "\" in file: \"" << fluxvect[0] << std::endl;
	throw;
	}
	fluxhist->SetDirectory(0);
	}
	} else if (fluxvect.size() == 1) {
	MonoEnergy = GeneralUtils::StrToDbl(fluxvect[0]);
	RunGENIEPrepareMono(input, target, output);
	return;
	} else {
	LOG(FTL) << "Bad flux specification: \"" << flux << "\"." << std::endl;
	throw;
	}

	// Setup TTree
	TChain* tn = new TChain("gtree");

	if (input.find_first_of(',') != std::string::npos) {
	std::vector<std::string> inputvect = GeneralUtils::ParseToStr(input, ",");

	for (size_t iv_it = 0; iv_it < inputvect.size(); ++iv_it) {
	tn->AddFile(inputvect[iv_it].c_str());
	QLOG(FIT, "Added input file: " << inputvect[iv_it]);
	}
	} else { // The Add form can accept wildcards.
	tn->Add(input.c_str());
	}

	int nevt = tn->GetEntries();

	if (!nevt) {
	THROW("Couldn't load any events from input specification: \""
	<< input.c_str() << "\"");
	} else {
	QLOG(FIT, "Found " << nevt << " input entries.");
	}

	NtpMCEventRecord* genientpl = NULL;
	tn->SetBranchAddress("gmcrec", &genientpl);

	// Make Event Hist
	TH1D* eventhist = (TH1D*)fluxhist->Clone();
	eventhist->Reset();

	TH1D* xsechist = (TH1D*)eventhist->Clone();

	// Create maps
	std::map<std::string, TH1D*> modexsec;
	std::map<std::string, TH1D*> modecount;
	std::vector<std::string> genieids;
	std::vector<std::string> targetids;
	std::vector<std::string> interids;

	// Loop over all events
	for (int i = 0; i < nevt; i++) {
	tn->GetEntry(i);

	StopTalking();
	EventRecord& event = *(genientpl->event);
	GHepParticle* neu = event.Probe();
	StartTalking();

	// Get XSec From Spline
	GHepRecord genie_record = static_cast<GHepRecord>(event);
	double xsec = (genie_record.XSec() / (1E-38 * genie::units::cm2));

	// Parse Interaction String
	std::string mode = genie_record.Summary()->AsString();
	std::vector<std::string> modevec = GeneralUtils::ParseToStr(mode, ";");
	std::string targ = (modevec[0] + ";" + modevec[1]);
	std::string inter = mode;

	// Fill lists of Unique IDS
	if (std::find(targetids.begin(), targetids.end(), targ) ==
	targetids.end()) {
	targetids.push_back(targ);
	}

	if (std::find(interids.begin(), interids.end(), inter) == interids.end()) {
	interids.push_back(inter);
	}

	// Create entries Mode Maps
	if (modexsec.find(mode) == modexsec.end()) {
	genieids.push_back(mode);

	modexsec[mode] = (TH1D*)xsechist->Clone();
	modecount[mode] = (TH1D*)xsechist->Clone();
	}

	// Fill XSec Histograms
	modexsec[mode]->Fill(neu->E(), xsec);
	modecount[mode]->Fill(neu->E());

	// Fill total event hist
	eventhist->Fill(neu->E());

	if (i % (nevt / 20) == 0) {
	LOG(FIT) << "Processed " << i << "/" << nevt
	<< " GENIE events (Last event: { E: " << neu->E()
	<< ", xsec: " << xsec << " }." << std::endl;
	}

	// Clear Event
	genientpl->Clear();
	}
	LOG(FIT) << "Processed all events" << std::endl;

	// Once event loop is done we can start saving stuff into the file

	TFile* outputfile;

	if (!gOutputFile.length()) {
	tn->GetEntry(0);
	outputfile = tn->GetFile();
	outputfile->cd();
	} else {
	outputfile = new TFile(gOutputFile.c_str(), "RECREATE");
	outputfile->cd();

	QLOG(FIT, "Cloning input vector to output file: " << gOutputFile);
	TTree* cloneTree = tn->CloneTree();
	cloneTree->SetDirectory(outputfile);
	cloneTree->Write();
	QLOG(FIT, "Done.");
	}

	LOG(FIT) << "Getting splines " << std::endl;

	// Save each of the reconstructed splines to file
	std::map<std::string, TH1D*> modeavg;

	TDirectory* inddir = (TDirectory*)outputfile->Get("IndividualGENIESplines");
	if (!inddir)
	inddir = (TDirectory*)outputfile->mkdir("IndividualGENIESplines");
	inddir->cd();

	// Loop over GENIE ID's and get MEC count
	int MECcount = 0;
	bool MECcorrect = FitPar::Config().GetParB("CorrectGENIEMECNorm");
	for (UInt_t i = 0; i < genieids.size(); i++) {
	if (genieids[i].find("MEC") != std::string::npos) {
	MECcount++;
	}
	}
	LOG(FIT) << "Found " << MECcount << " repeated MEC instances." << std::endl;

	for (UInt_t i = 0; i < genieids.size(); i++) {
	std::string mode = genieids[i];

	modexsec[mode]->Write((mode + "_summed_xsec").c_str(), TObject::kOverwrite);
	modecount[mode]->Write((mode + "_summed_evt").c_str(), TObject::kOverwrite);

	// Form extra avg xsec map -> Reconstructed spline
	modeavg[mode] = (TH1D*)modexsec[mode]->Clone();
	modeavg[mode]->Divide(modecount[mode]);

	if (MECcorrect && (mode.find("MEC") != std::string::npos)) {
	modeavg[mode]->Scale(1.0 / double(MECcount));
	}

	modeavg[mode]->Write((mode + "_rec_spline").c_str(), TObject::kOverwrite);
	}

	TDirectory* targdir = (TDirectory*)outputfile->Get("TargetGENIESplines");
	if (!targdir) targdir = (TDirectory*)outputfile->mkdir("TargetGENIESplines");
	targdir->cd();

	LOG(FIT) << "Getting Target Splines" << std::endl;
	// For each target save a total spline
	std::map<std::string, TH1D*> targetsplines;

	for (uint i = 0; i < targetids.size(); i++) {
	LOG(FIT) << "Getting target " << i << std::endl;
	std::string targ = targetids[i];
	targetsplines[targ] = (TH1D*)xsechist->Clone();
	LOG(FIT) << "Created target spline for " << targ << std::endl;

	for (uint j = 0; j < genieids.size(); j++) {
	std::string mode = genieids[j];

	// Look at all matching modes/targets
	if (mode.find(targ) != std::string::npos) {
	LOG(FIT) << "Mode " << mode << " contains " << targ << " target!"
	<< std::endl;
	// modeavg[mode]->Write( (mode + "_cont_" + targ).c_str() ,
	// TObject::kOverwrite);
	targetsplines[targ]->Add(modeavg[mode]);
	LOG(FIT) << "Finished with Mode " << mode << " "
	<< modeavg[mode]->Integral() << std::endl;
	}
	}

	LOG(FIT) << "Saving target spline:" << targ << std::endl;
	targetsplines[targ]->Write(("Total" + targ).c_str(), TObject::kOverwrite);
	}

	LOG(FIT) << "Getting total splines" << std::endl;
	// Now we have each of the targets we need to create a total cross-section.
	int totalnucl = 0;
	std::vector<std::string> targprs = GeneralUtils::ParseToStr(target, ",");
	TH1D* totalxsec = (TH1D*)xsechist->Clone();

	for (uint i = 0; i < targprs.size(); i++) {
	std::string targpdg = targprs[i];

	for (std::map<std::string, TH1D*>::iterator iter = targetsplines.begin();
	iter != targetsplines.end(); iter++) {
	std::string targstr = iter->first;
	TH1D* xsec = iter->second;

	if (targstr.find(targpdg) != std::string::npos) {
	LOG(FIT) << "Adding target spline " << targstr
	<< " Integral = " << xsec->Integral("width") << std::endl;
	totalxsec->Add(xsec);

	int nucl = atoi(targpdg.c_str());
	totalnucl += int((nucl % 10000) / 10);
	}
	}
	}
	LOG(FIT) << "Total XSec Integral = " << totalxsec->Integral("width")
	<< std::endl;

	outputfile->cd();
	totalxsec->Write("nuisance_xsec", TObject::kOverwrite);
	eventhist = (TH1D*)fluxhist->Clone();
	eventhist->Multiply(totalxsec);

	LOG(FIT) << "Dividing by Total Nucl = " << totalnucl << std::endl;
	eventhist->Scale(1.0 / double(totalnucl));

	eventhist->Write("nuisance_events", TObject::kOverwrite);
	fluxhist->Write("nuisance_flux", TObject::kOverwrite);

	LOG(FIT) << "Inclusive XSec Per Nucleon = "
	<< eventhist->Integral("width") * 1E-38 / fluxhist->Integral("width")
	<< std::endl;
	std::cout << "XSec Hist Integral = " << xsechist->Integral("width")
	<< std::endl;

	outputfile->Write();
	outputfile->Close();
	delete outputfile;

	return;
	};

	void PrintOptions() {
	std::cout << "PrepareGENIEEvents NUISANCE app. " << std::endl
	<< "Takes GHep Outputs and prepares events for NUISANCE."
	<< std::endl
	<< std::endl
	<< "PrepareGENIEEvents [-h,-help,--h,--help] [-i "
	"inputfile1.root,inputfile2.root,inputfile3.root,...] "
	<< "[-f flux_root_file.root,flux_hist_name] [-t "
	"target1[frac1],target2[frac2],...]"
	<< std::endl
	<< std::endl;

	std::cout << "Prepare Mode [Default] : Takes a single GHep file, "
	"reconstructs the original GENIE splines, "
	<< " and creates a duplicate file that also contains the flux, "
	"event rate, and xsec predictions that NUISANCE needs. "
	<< std::endl;
	std::cout << "Following options are required for Prepare Mode:" << std::endl;
	std::cout << " [ -i inputfile.root ] : Reads in a single GHep input file "
	"that needs the xsec calculation ran on it. "
	<< std::endl;
	std::cout << " [ -f flux_file.root,hist_name ] : Path to root file "
	"containing the flux histogram the GHep records were generated "
	"with."
	<< " A simple method is to point this to the flux histogram genie "
	"generatrs '-f /path/to/events/input-flux.root,spectrum'. "
	<< std::endl;
	std::cout << " [ -f elow,ehigh,estep ] : Energy range specification when no "
	"flux file was used."
	<< std::endl;
	std::cout << " [ -t target ] : Target that GHepRecords were generated with. "
	"Comma seperated list. E.g. for CH2 "
	"target=1000060120,1000010010,1000010010"
	<< std::endl;
	std::cout << " [ -o outputfile.root ] : File to write prepared input file to."
	<< std::endl;
	std::cout << " [ -m Mono_E_nu_GeV ] : Run in mono-energetic mode."
	<< std::endl;
	}

	void ParseOptions(int argc, char* argv[]) {
	bool flagopt = false;

	// If No Arguments print commands
	for (int i = 1; i < argc; ++i) {
	if (!std::strcmp(argv[i], "-h")) {
	flagopt = true;
	break;
	}
	if (i + 1 != argc) {
	// Cardfile
	if (!std::strcmp(argv[i], "-h")) {
	flagopt = true;
	break;
	} else if (!std::strcmp(argv[i], "-i")) {
	gInputFiles = argv[i + 1];
	++i;
	} else if (!std::strcmp(argv[i], "-o")) {
	gOutputFile = argv[i + 1];
	++i;
	} else if (!std::strcmp(argv[i], "-f")) {
	gFluxFile = argv[i + 1];
	++i;
	} else if (!std::strcmp(argv[i], "-t")) {
	gTarget = argv[i + 1];
	++i;
	} else if (!std::strcmp(argv[i], "-m")) {
	MonoEnergy = GeneralUtils::StrToDbl(argv[i + 1]);
	IsMonoE = true;
	++i;
	} else {
	ERR(FTL) << "ERROR: unknown command line option given! - '" << argv[i]
	<< " " << argv[i + 1] << "'" << std::endl;
	PrintOptions();
	break;
	}
	}
	}

	if (gInputFiles == "" && !flagopt) {
	ERR(FTL) << "No input file(s) specified!" << std::endl;
	flagopt = true;
	}

	- if (gFluxFile == "" && !flagopt) {
	+ if (gFluxFile == "" && !flagopt && !IsMonoE) {
	ERR(FTL) << "No flux input specified for Prepare Mode" << std::endl;
	flagopt = true;
	}

	if (gTarget == "" && !flagopt) {
	ERR(FTL) << "No target specified for Prepare Mode" << std::endl;
	flagopt = true;
	}

	if (argc < 1 \|\| flagopt) {
	PrintOptions();
	exit(-1);
	}

	return;
	}
	diff --git a/src/MCStudies/GenericFlux_Vectors.cxx b/src/MCStudies/GenericFlux_Vectors.cxx
	index eb75ee7..94d6d42 100644
	--- a/src/MCStudies/GenericFlux_Vectors.cxx
	+++ b/src/MCStudies/GenericFlux_Vectors.cxx
	@@ -1,237 +1,239 @@
	// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret

	/*******************************************************************************
	* This file is part of NUISANCE.
	*
	* NUISANCE is free software: you can redistribute it and/or modify
	* it under the terms of the GNU General Public License as published by
	* the Free Software Foundation, either version 3 of the License, or
	* (at your option) any later version.
	*
	* NUISANCE is distributed in the hope that it will be useful,
	* but WITHOUT ANY WARRANTY; without even the implied warranty of
	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	* GNU General Public License for more details.
	*
	* You should have received a copy of the GNU General Public License
	* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
	*******************************************************************************/

	#include "GenericFlux_Vectors.h"

	GenericFlux_Vectors::GenericFlux_Vectors(std::string name, std::string inputfile,
	FitWeight *rw, std::string type,
	std::string fakeDataFile) {
	// Measurement Details
	fName = name;
	eventVariables = NULL;

	// Define our energy range for flux calcs
	EnuMin = 0.;
	EnuMax = 100.; // Arbritrarily high energy limit

	// Set default fitter flags
	fIsDiag = true;
	fIsShape = false;
	fIsRawEvents = false;

	// This function will sort out the input files automatically and parse all the
	// inputs,flags,etc.
	// There may be complex cases where you have to do this by hand, but usually
	// this will do.
	Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile);

	eventVariables = NULL;

	// Setup fDataHist as a placeholder
	this->fDataHist = new TH1D(("empty_data"), ("empty-data"), 1, 0, 1);
	this->SetupDefaultHist();
	fFullCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist);
	covar = StatUtils::GetInvert(fFullCovar);

	// 1. The generator is organised in SetupMeasurement so it gives the
	// cross-section in "per nucleon" units.
	// So some extra scaling for a specific measurement may be required. For
	// Example to get a "per neutron" measurement on carbon
	// which we do here, we have to multiple by the number of nucleons 12 and
	// divide by the number of neutrons 6.
	- this->fScaleFactor =
	+ fScaleFactor =
	(GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) /
	this->TotalIntegratedFlux();

	LOG(SAM) << " Generic Flux Scaling Factor = " << fScaleFactor << std::endl;

	if (fScaleFactor <= 0.0) {
	ERR(WRN) << "SCALE FACTOR TOO LOW " << std::endl;
	}

	// Setup our TTrees
	this->AddEventVariablesToTree();
	}

	void GenericFlux_Vectors::AddEventVariablesToTree() {
	// Setup the TTree to save everything
	if (!eventVariables) {
	Config::Get().out->cd();
	eventVariables = new TTree((this->fName + "_VARS").c_str(),
	(this->fName + "_VARS").c_str());
	}

	LOG(SAM) << "Adding Event Variables" << std::endl;

	eventVariables->Branch("Mode", &Mode, "Mode/I" );
	eventVariables->Branch("cc", &cc, "cc/B" );
	eventVariables->Branch("PDGnu", &PDGnu, "PDGnu/I" );
	+ eventVariables->Branch("Enu_true", &Enu_true, "Enu_true/F" );
	eventVariables->Branch("tgt", &tgt, "tgt/I" );
	eventVariables->Branch("PDGLep", &PDGLep, "PDGLep/I");
	eventVariables->Branch("ELep", &ELep, "ELep/F" );
	eventVariables->Branch("CosLep", &CosLep, "CosLep/F");

	// Basic interaction kinematics
	eventVariables->Branch("Q2", &Q2, "Q2/F");
	eventVariables->Branch("q0", &q0, "q0/F");
	eventVariables->Branch("q3", &q3, "q3/F");
	eventVariables->Branch("Enu_QE", &Enu_QE, "Enu_QE/F");
	eventVariables->Branch("Q2_QE", &Q2_QE, "Q2_QE/F");
	eventVariables->Branch("W_nuc_rest", &W_nuc_rest, "W_nuc_rest/F");
	eventVariables->Branch("W", &W, "W/F");
	eventVariables->Branch("x", &x, "x/F");
	eventVariables->Branch("y", &y, "y/F");

	// Save outgoing particle vectors
	eventVariables->Branch("nfsp", &nfsp, "nfsp/I");
	eventVariables->Branch("px", px, "px[nfsp]/F");
	eventVariables->Branch("py", py, "py[nfsp]/F");
	eventVariables->Branch("pz", pz, "pz[nfsp]/F");
	eventVariables->Branch("E", E, "E[nfsp]/F");
	eventVariables->Branch("pdg", pdg, "pdg[nfsp]/I");

	// Event Scaling Information
	eventVariables->Branch("Weight", &Weight, "Weight/F");
	eventVariables->Branch("InputWeight", &InputWeight, "InputWeight/F");
	eventVariables->Branch("RWWeight", &RWWeight, "RWWeight/F");
	eventVariables->Branch("fScaleFactor", &fScaleFactor, "fScaleFactor/F");

	return;
	}


	void GenericFlux_Vectors::FillEventVariables(FitEvent *event) {

	// Reset all Function used to extract any variables of interest to the event
	- PDGnu = tgt = PDGLep = 0;
	+ Mode = PDGnu = tgt = PDGLep = 0;

	- Enu = ELep = CosLep = Q2 = q0 = q3 = Enu_QE = Q2_QE = W_nuc_rest = W = x = y =
	- -999.9;
	+ Enu_true = ELep = CosLep = Q2 = q0 = q3 = Enu_QE = Q2_QE = W_nuc_rest = W = x = y = -999.9;

	nfsp = 0;
	for (int i = 0; i < kMAX; ++i){
	px[i] = py[i] = pz[i] = E[i] = -999;
	pdg[i] = 0;
	}

	Weight = InputWeight = RWWeight = 0;

	partList.clear();

	// Now fill the information
	Mode = event->Mode;
	cc = (abs(event->Mode) < 30);

	// Get the incoming neutrino and outgoing lepton
	FitParticle *nu = event->GetNeutrinoIn();
	FitParticle *lep = event->GetHMFSAnyLepton();

	PDGnu = nu->fPID;
	- Enu = nu->fP.E()/1E3;
	+ Enu_true = nu->fP.E()/1E3;
	tgt = event->fTargetPDG;
	- PDGLep = lep->fPID;
	- ELep = lep->fP.E()/1E3;
	- CosLep = cos(nu->fP.Vect().Angle(lep->fP.Vect()));
	-
	- // Basic interaction kinematics
	- Q2 = -1*(nu->fP - lep->fP).Mag2()/1E6;
	- q0 = (nu->fP - lep->fP).E()/1E3;
	- q3 = (nu->fP - lep->fP).Vect().Mag()/1E3;
	-
	- // Get W_true with assumption of initial state nucleon at rest
	- float m_n = (float)PhysConst::mass_proton;
	- W_nuc_rest = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n);
	- W = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n);
	- x = Q2/(2 * m_n * q0);
	- y = 1 - ELep/Enu;
	-
	- // These assume C12 binding from MINERvA... not ideal
	- Q2_QE = FitUtils::Q2QErec(lep->fP, CosLep, 34., true);
	- Enu_QE = FitUtils::EnuQErec(lep->fP, CosLep, 34., true);
	+ if (lep != NULL) {
	+ PDGLep = lep->fPID;
	+ ELep = lep->fP.E()/1E3;
	+ CosLep = cos(nu->fP.Vect().Angle(lep->fP.Vect()));
	+
	+ // Basic interaction kinematics
	+ Q2 = -1*(nu->fP - lep->fP).Mag2()/1E6;
	+ q0 = (nu->fP - lep->fP).E()/1E3;
	+ q3 = (nu->fP - lep->fP).Vect().Mag()/1E3;
	+
	+ // These assume C12 binding from MINERvA... not ideal
	+ Enu_QE = FitUtils::EnuQErec(lep->fP, CosLep, 34., true);
	+ Q2_QE = FitUtils::Q2QErec(lep->fP, CosLep, 34., true);
	+
	+ // Get W_true with assumption of initial state nucleon at rest
	+ float m_n = (float)PhysConst::mass_proton;
	+ W_nuc_rest = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n);
	+ W = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n);
	+ x = Q2/(2 * m_n * q0);
	+ y = 1 - ELep/Enu_true;
	+ }

	// Loop over the particles and store all the final state particles in a vector
	for (UInt_t i = 0; i < event->Npart(); ++i) {

	bool part_alive = event->PartInfo(i)->fIsAlive and event->PartInfo(i)->Status() == kFinalState;
	if (!part_alive) continue;

	partList .push_back(event->PartInfo(i));
	}

	// Save outgoing particle vectors
	nfsp = (int)partList.size();

	for (int i = 0; i < nfsp; ++i){
	px[i] = partList[i]->fP.X()/1E3;
	py[i] = partList[i]->fP.Y()/1E3;
	pz[i] = partList[i]->fP.Z()/1E3;
	E[i] = partList[i]->fP.E()/1E3;
	pdg[i] = partList[i]->fPID;
	}

	// Fill event weights
	Weight = event->RWWeight * event->InputWeight;
	RWWeight = event->RWWeight;
	InputWeight = event->InputWeight;

	// Fill the eventVariables Tree
	eventVariables->Fill();
	return;
	};

	void GenericFlux_Vectors::Write(std::string drawOpt) {

	// First save the TTree
	eventVariables->Write();

	// Save Flux and Event Histograms too
	GetInput()->GetFluxHistogram()->Write();
	GetInput()->GetEventHistogram()->Write();

	return;
	}


	// Override functions which aren't really necessary
	bool GenericFlux_Vectors::isSignal(FitEvent *event) {
	(void)event;
	return true;
	};

	void GenericFlux_Vectors::ScaleEvents() {
	return;
	}

	void GenericFlux_Vectors::ApplyNormScale(float norm) {
	this->fCurrentNorm = norm;
	return;
	}

	void GenericFlux_Vectors::FillHistograms() {
	return;
	}

	void GenericFlux_Vectors::ResetAll() {
	eventVariables->Reset();
	return;
	}

	float GenericFlux_Vectors::GetChi2() {
	return 0.0;
	}
	diff --git a/src/MCStudies/GenericFlux_Vectors.h b/src/MCStudies/GenericFlux_Vectors.h
	index f09e7e9..8b9b7a0 100644
	--- a/src/MCStudies/GenericFlux_Vectors.h
	+++ b/src/MCStudies/GenericFlux_Vectors.h
	@@ -1,98 +1,98 @@
	// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret

	/*******************************************************************************
	* This file is part of NUISANCE.
	*
	* NUISANCE is free software: you can redistribute it and/or modify
	* it under the terms of the GNU General Public License as published by
	* the Free Software Foundation, either version 3 of the License, or
	* (at your option) any later version.
	*
	* NUISANCE is distributed in the hope that it will be useful,
	* but WITHOUT ANY WARRANTY; without even the implied warranty of
	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	* GNU General Public License for more details.
	*
	* You should have received a copy of the GNU General Public License
	* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
	*******************************************************************************/

	#ifndef GenericFlux_Vectors_H_SEEN
	#define GenericFlux_Vectors_H_SEEN
	#include "Measurement1D.h"

	class GenericFlux_Vectors : public Measurement1D {

	public:

	GenericFlux_Vectors(std::string name, std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile);
	virtual ~GenericFlux_Vectors() {};

	//! Grab info from event
	void FillEventVariables(FitEvent *event);

	//! Fill Custom Histograms
	void FillHistograms();

	//! ResetAll
	void ResetAll();

	//! Scale
	void ScaleEvents();

	//! Norm
	void ApplyNormScale(float norm);

	//! Define this samples signal
	bool isSignal(FitEvent *nvect);

	//! Write Files
	void Write(std::string drawOpt);

	//! Get Chi2
	float GetChi2();

	void AddEventVariablesToTree();

	private:

	TTree* eventVariables;
	std::vector<FitParticle*> partList;

	int Mode;
	bool cc;
	int PDGnu;
	- float Enu;
	int tgt;
	int PDGLep;
	float ELep;
	float CosLep;

	// Basic interaction kinematics
	float Q2;
	float q0;
	float q3;
	float Enu_QE;
	+ float Enu_true;
	float Q2_QE;
	float W_nuc_rest;
	float W;
	float x;
	float y;

	// Save outgoing particle vectors
	int nfsp;
	static const int kMAX = 200;
	float px[kMAX];
	float py[kMAX];
	float pz[kMAX];
	float E[kMAX];
	int pdg[kMAX];

	// Basic event info
	float Weight;
	float InputWeight;
	float RWWeight;
	float fScaleFactor;
	};

	#endif

File Metadata

Mime Type: text/x-diff
Expires: Wed, May 14, 11:07 AM (21 h, 4 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 5088632
Default Alt Text: (31 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions