Page MenuHomeHEPForge
Files F7878304

No OneTemporary
Actions

View Options

diff --git a/app/PrepareGENIE.cxx b/app/PrepareGENIE.cxx
index fb48a61..68c1477 100644
--- a/app/PrepareGENIE.cxx
+++ b/app/PrepareGENIE.cxx
@@ -1,356 +1,348 @@
#include <stdio.h>
#include <stdlib.h>
#include "TFile.h"
#include "TH1D.h"
#include "TTree.h"
#include "PlotUtils.h"
#include "FitLogger.h"
#ifdef __GENIE_ENABLED__
#include "Conventions/Units.h"
#endif
bool gFlagMerge = false;
std::string gInputFiles = "";
std::string gOutputFile = "";
std::string gFluxFile = "";
std::string gTarget = "";
void PrintOptions();
void ParseOptions(int argc, char* argv[]);
void RunGENIEMerger(std::string inputs, 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 (gFlagMerge) RunGENIEMerger(gInputFiles, gOutputFile);
else RunGENIEPrepare(gInputFiles, gFluxFile, gTarget, gOutputFile);
}
void RunGENIEMerger(std::string inputs, std::string output) {
};
void RunGENIEPrepare(std::string input, std::string flux, std::string target, std::string output) {
LOG(FIT) << "Running GENIE Prepare" << 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);
// Get Flux Hist
std::vector<std::string> fluxvect = GeneralUtils::ParseToStr(flux, ",");
TH1D* fluxhist = NULL;
if (fluxvect.size() > 1) {
TFile* fluxfile = new TFile(fluxvect[0].c_str(), "READ");
if (!fluxfile->IsZombie()) {
fluxhist = (TH1D*) fluxfile->Get(fluxvect[1].c_str());
fluxhist->SetDirectory(0);
} else {
// Function with EnuRange
if (fluxvect.size() == 3){
ERR(FTL) << "FUNCTION WITH ENU RANGE NOT SUPPORTED SORRY!" << std::endl;
throw;
// Single Enu
}
}
} else if (fluxvect.size() == 1){
double E = GeneralUtils::StrToDbl(fluxvect[0]);
fluxhist = new TH1D("fluxhist","fluxhist",1, E-0.00001, E+0.00001);
fluxhist->SetBinContent(1, 1.0);
-
-
// if (fluxvect[0] == '1.0') {
// fluxhist = new TH1D("fluxhist", "fluxhist", 40, GeneralUtils::StrToDbl(fluxvect[1]), GeneralUtils::StrToDbl(fluxvect[2]));
// for (int i = 0; i < fluxhist->GetNbinsX(); i++) {
// fluxhist->SetBinContent(i + 1, 1.0);
// }
// } else {
// // TF1 f1 = TF1("1.0", GeneralUtils::StrToDbl(fluxvect[1]), GeneralUtils::StrToDbl(fluxvect[2]), 1);
// // std::cout << "Created TF1 with '" << fluxvect[0].c_str()<<"' " << GeneralUtils::StrToDbl(fluxvect[1]) << " " << GeneralUtils::StrToDbl(fluxvect[2]) << std::endl;
// fluxhist = new TH1D("fluxhist", "fluxhist", 40, GeneralUtils::StrToDbl(fluxvect[1]), GeneralUtils::StrToDbl(fluxvect[2]));
// for (int i = 0; i < fluxhist->GetNbinsX(); i++) {
// fluxhist->SetBinContent(i + 1, 1.0); //f1.Eval(fluxhist->GetXaxis()->GetBinCenter(i+1)));
// // std::cout << "Filling Flux Hist " << f1.Eval(fluxhist->GetXaxis()->GetBinCenter(i+1)) << std::endl;
// }
// sleep(10);
} else {
LOG(FTL) << "NO FLUX SPECIFIED" << std::endl;
throw;
}
// 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());
- std::cout << "Electron Energy = " << neu->E() << std::endl;
// 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;
// Once event loop is done we can start saving stuff into the file
// bool savesplines = FitPar::Config().GetParB("save_genie_splines"); // Currently not implemented
TFile* outputfile = new TFile(input.c_str(), "UPDATE");
outputfile->cd();
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();
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]);
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 << " : " << xsec->Integral("width") << std::endl;
+ 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*)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;
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 << " [ -t target ] : Target that GHepRecords were generated with. Comma seperated list. E.g. for CH2 target=1000060120,1000010010,1000010010" << std::endl;
/*
std::cout << "Merger Mode [activate with -m] : Takes the list of input files assuming 'Prepare Mode' has already been ran on them and merges them "
<< "into a single file with associated Friend Tree to help with conserving ratios of events (e.g. adding nue and nueb beams together into a single file" << std::endl;
std::cout << "Following optoins are required for Merger Mode:" << std::endl;
std::cout << " [ -i inputfile1.root,inputfile2.root ] : Comma Seperated list of files to be merged. " << std::endl;
std::cout << " [ -o outputfile.root ] : Output file for merger." << 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 (!std::strcmp(argv[i], "-m")) { gFlagMerge = 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 {
ERR(FTL) << "ERROR: unknown command line option given! - '"
<< argv[i] << " " << argv[i + 1] << "'" << std::endl;
PrintOptions();
break;
}
}
}
/*
if (gOutputFile == "" && !flagopt){
ERR(FTL) << "No output file specificed!" << std::endl;
flagopt = true;
}
*/
if (gInputFiles == "" && !flagopt) {
ERR(FTL) << "No input file(s) specified!" << std::endl;
flagopt = true;
}
if (!gFlagMerge && gFluxFile == "" && !flagopt) {
ERR(FTL) << "No flux input specified for Prepare Mode" << std::endl;
flagopt = true;
}
if (!gFlagMerge && 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/app/gEvGen_NUISANCE.cxx b/app/gEvGen_NUISANCE.cxx
index b23d07f..ba402dc 100644
--- a/app/gEvGen_NUISANCE.cxx
+++ b/app/gEvGen_NUISANCE.cxx
@@ -1,919 +1,960 @@
+#ifdef __GENIE_ENABLED__
//____________________________________________________________________________
/*!
\program gevgen
\brief A simple 'generic' GENIE v+A event generation driver (gevgen).
It handles:
a) event generation for a fixed init state (v+A) at fixed energy, or
b) event generation for simple fluxes (specified either via some
functional form, tabular text file or a ROOT histogram) and for
simple 'geometries' (a target mix with its corresponding weights)
See the GENIE manual for other apps handling experiment-specific
event generation cases using the outputs of detailed neutrino flux
simulations and realistic detector geometry descriptions.
Syntax :
gevgen [-h]
[-r run#]
-n nev
-e energy (or energy range)
-p neutrino_pdg
-t target_pdg
[-f flux_description]
[-w]
[--seed random_number_seed]
[--cross-sections xml_file]
[--event-generator-list list_name]
[--tune genie_tune]
[--message-thresholds xml_file]
[--unphysical-event-mask mask]
[--event-record-print-level level]
[--mc-job-status-refresh-rate rate]
[--cache-file root_file]
Options :
[] Denotes an optional argument.
-h
Prints-out help on using gevgen and exits.
-n
Specifies the number of events to generate.
-r
Specifies the MC run number.
-e
Specifies the neutrino energy.
If what follows the -e option is a comma separated pair of values
it will be interpreted as an energy range for the flux specified
via the -f option (see below).
-p
Specifies the neutrino PDG code.
-t
Specifies the target PDG code (pdg format: 10LZZZAAAI) _or_ a target
mix (pdg codes with corresponding weights) typed as a comma-separated
list of pdg codes with the corresponding weight fractions in brackets,
eg code1[fraction1],code2[fraction2],...
For example, to use a target mix of 95% O16 and 5% H type:
`-t 1000080160[0.95],1000010010[0.05]'.
-f
Specifies the neutrino flux spectrum.
It can be any of:
-- A function:
eg ` -f x*x+4*exp(-x)'
-- A vector file:
The vector file should contain 2 columns corresponding to
energy,flux (see $GENIE/data/flux/ for few examples).
-- A 1-D ROOT histogram (TH1D):
The general syntax is `-f /full/path/file.root,object_name'
-w
Forces generation of weighted events.
This option is relevant only if a neutrino flux is specified.
Note that 'weighted' refers to the selection of the primary
flux neutrino + target that were forced to interact. A weighting
scheme for the generated kinematics of individual processes can
still be in effect if enabled..
** Only use that option if you understand what it means **
--seed
Random number seed.
--cross-sections
Name (incl. full path) of an XML file with pre-computed
cross-section values used for constructing splines.
--event-generator-list
List of event generators to load in event generation drivers.
[default: "Default"].
--tune
Specifies a GENIE comprehensive neutrino interaction model tune.
[default: "Default"].
--message-thresholds
Allows users to customize the message stream thresholds.
The thresholds are specified using an XML file.
See $GENIE/config/Messenger.xml for the XML schema.
--unphysical-event-mask
Allows users to specify a 16-bit mask to allow certain types of
unphysical events to be written in the output file.
[default: all unphysical events are rejected]
--event-record-print-level
Allows users to set the level of information shown when the event
record is printed in the screen. See GHepRecord::Print().
--mc-job-status-refresh-rate
Allows users to customize the refresh rate of the status file.
--cache-file
Allows users to specify a cache file so that the cache can be
re-used in subsequent MC jobs.
*** See the User Manual for more details and examples. ***
\author Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
University of Liverpool & STFC Rutherford Appleton Lab
\created October 05, 2004
\cpright Copyright (c) 2003-2016, GENIE Neutrino MC Generator Collaboration
For the full text of the license visit http://copyright.genie-mc.org
or see $GENIE/LICENSE
*/
//____________________________________________________________________________
#include <cstdlib>
#include <cassert>
#include <sstream>
#include <string>
#include <vector>
#include <map>
#if defined(HAVE_FENV_H) && defined(HAVE_FEENABLEEXCEPT)
#include <fenv.h> // for `feenableexcept`
#endif
#include <TFile.h>
#include <TTree.h>
#include <TSystem.h>
#include <TVector3.h>
#include <TH1.h>
#include <TF1.h>
#include "Conventions/XmlParserStatus.h"
#include "Conventions/GBuild.h"
#include "Conventions/Controls.h"
+#include "Conventions/Units.h"
+
#include "EVGCore/EventRecord.h"
#include "EVGDrivers/GFluxI.h"
#include "EVGDrivers/GEVGDriver.h"
//#include "EVGDrivers/GMCJDriver.h"
#include "GNUISANCEMCJDriver.h"
#include "EVGDrivers/GMCJMonitor.h"
#include "Interaction/Interaction.h"
#include "Messenger/Messenger.h"
#include "Ntuple/NtpWriter.h"
#include "Ntuple/NtpMCFormat.h"
#include "Numerical/RandomGen.h"
#include "Numerical/Spline.h"
#include "PDG/PDGCodes.h"
#include "PDG/PDGUtils.h"
#include "Utils/AppInit.h"
#include "Utils/RunOpt.h"
#include "Utils/XSecSplineList.h"
#include "Utils/StringUtils.h"
#include "Utils/PrintUtils.h"
#include "Utils/SystemUtils.h"
#include "Utils/CmdLnArgParser.h"
//#include "FitBase/GNUISANCEFlux.h"
#include "GNUISANCEFlux.h"
//#ifdef __GENIE_FLUX_DRIVERS_ENABLED__
//#ifdef __GENIE_GEOM_DRIVERS_ENABLED__
#define __CAN_GENERATE_EVENTS_USING_A_FLUX_OR_TGTMIX__
#include "GNUISANCEFlux.h"
#include "FluxDrivers/GMonoEnergeticFlux.h"
#include "Geo/PointGeomAnalyzer.h"
//#endif
//#endif
using std::string;
using std::vector;
using std::map;
using std::ostringstream;
using namespace genie;
using namespace genie::controls;
using namespace genie::flux;
void GetCommandLineArgs (int argc, char ** argv);
void Initialize (void);
void PrintSyntax (void);
#ifdef __CAN_GENERATE_EVENTS_USING_A_FLUX_OR_TGTMIX__
void GenerateEventsUsingFluxOrTgtMix();
GeomAnalyzerI * GeomDriver (void);
GFluxI * FluxDriver (void);
GFluxI * MonoEnergeticFluxDriver (void);
GFluxI * TH1FluxDriver (void);
string ConvertTargetIDs (string);
string ConvertFluxIDs (string);
void ListTargetIDs(void);
void ListFluxIDs(void);
#endif
void GenerateEventsAtFixedInitState (void);
//Default options (override them using the command line arguments):
int kDefOptNevents = 0; // n-events to generate
NtpMCFormat_t kDefOptNtpFormat = kNFGHEP; // ntuple format
Long_t kDefOptRunNu = 0; // default run number
//User-specified options:
int gOptNevents; // n-events to generate
double gOptNuEnergy; // neutrino E, or min neutrino energy in spectrum
double gOptNuEnergyRange;// energy range in input spectrum
int gOptNuPdgCode; // neutrino PDG code
map<int,double> gOptTgtMix; // target mix (each with its relative weight)
Long_t gOptRunNu; // run number
string gOptFlux; //
bool gOptWeighted; //
bool gOptUsingFluxOrTgtMix = false;
long int gOptRanSeed; // random number seed
string gOptInpXSecFile; // cross-section splines
vector<int> gOptNuPDGs; // Beam PDGS
vector<TH1D*> gOptNuFluxs; // Beam Fluxes
vector<int> gOptTgtPDGs; // Target PDGS
int gOptNumberNucleons; // Total number of target nucleons
//____________________________________________________________________________
int main(int argc, char ** argv)
{
GetCommandLineArgs(argc,argv);
Initialize();
// throw on NaNs and Infs...
#if defined(HAVE_FENV_H) && defined(HAVE_FEENABLEEXCEPT)
feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
#endif
//
// Generate neutrino events
//
if(gOptUsingFluxOrTgtMix) {
#ifdef __CAN_GENERATE_EVENTS_USING_A_FLUX_OR_TGTMIX__
GenerateEventsUsingFluxOrTgtMix();
#else
LOG("gevgen", pERROR)
<< "\n To be able to generate neutrino events from a flux and/or a target mix"
<< "\n you need to add the following config options at your GENIE installation:"
<< "\n --enable-flux-drivers --enable-geom-drivers \n" ;
#endif
} else {
GenerateEventsAtFixedInitState();
}
return 0;
}
//____________________________________________________________________________
void Initialize()
{
// Initialization of random number generators, cross-section table,
// messenger thresholds, cache file
utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
utils::app_init::CacheFile(RunOpt::Instance()->CacheFile());
utils::app_init::RandGen(gOptRanSeed);
utils::app_init::XSecTable(gOptInpXSecFile, false);
// Set GHEP print level
GHepRecord::SetPrintLevel(RunOpt::Instance()->EventRecordPrintLevel());
}
//____________________________________________________________________________
void GenerateEventsAtFixedInitState(void)
{
int neutrino = gOptNuPdgCode;
int target = gOptTgtMix.begin()->first;
double Ev = gOptNuEnergy;
TLorentzVector nu_p4(0.,0.,Ev,Ev); // px,py,pz,E (GeV)
// Create init state
InitialState init_state(target, neutrino);
// Create/config event generation driver
GEVGDriver evg_driver;
evg_driver.SetEventGeneratorList(RunOpt::Instance()->EventGeneratorList());
evg_driver.SetUnphysEventMask(*RunOpt::Instance()->UnphysEventMask());
evg_driver.Configure(init_state);
// Initialize an Ntuple Writer
NtpWriter ntpw(kDefOptNtpFormat, gOptRunNu);
ntpw.Initialize();
// Create an MC Job Monitor
GMCJMonitor mcjmonitor(gOptRunNu);
mcjmonitor.SetRefreshRate(RunOpt::Instance()->MCJobStatusRefreshRate());
LOG("gevgen", pNOTICE)
<< "\n ** Will generate " << gOptNevents << " events for \n"
<< init_state << " at Ev = " << Ev << " GeV";
// Generate events / print the GHEP record / add it to the ntuple
int ievent = 0;
while (ievent < gOptNevents) {
LOG("gevgen", pNOTICE)
<< " *** Generating event............ " << ievent;
// generate a single event
EventRecord * event = evg_driver.GenerateEvent(nu_p4);
if(!event) {
LOG("gevgen", pNOTICE)
<< "Last attempt failed. Re-trying....";
continue;
}
LOG("gevgen", pNOTICE)
<< "Generated Event GHEP Record: " << *event;
// add event at the output ntuple, refresh the mc job monitor & clean up
ntpw.AddEventRecord(ievent, event);
mcjmonitor.Update(ievent,event);
ievent++;
delete event;
}
// Save the generated MC events
ntpw.Save();
}
//____________________________________________________________________________
#ifdef __CAN_GENERATE_EVENTS_USING_A_FLUX_OR_TGTMIX__
//............................................................................
void GenerateEventsUsingFluxOrTgtMix(void)
{
// Get flux and geom drivers
GFluxI * flux_driver = FluxDriver();
GeomAnalyzerI * geom_driver = GeomDriver();
// Create the monte carlo job driver
GNUISANCEMCJDriver * mcj_driver = new GNUISANCEMCJDriver;
mcj_driver->SetEventGeneratorList(RunOpt::Instance()->EventGeneratorList());
mcj_driver->SetUnphysEventMask(*RunOpt::Instance()->UnphysEventMask());
mcj_driver->UseFluxDriver(flux_driver);
mcj_driver->UseGeomAnalyzer(geom_driver);
mcj_driver->Configure();
mcj_driver->UseSplines();
if(!gOptWeighted)
mcj_driver->ForceSingleProbScale();
// Initialize an Ntuple Writer to save GHEP records into a TTree
NtpWriter ntpw(kDefOptNtpFormat, gOptRunNu);
ntpw.Initialize();
// Now Create flux histograms
TH1D* totalflux = static_cast<GNUISANCEFlux*>(flux_driver)->GetTotalSpectrum();
vector<TH1D*> fluxinputs = static_cast<GNUISANCEFlux*>(flux_driver)->GetSpectrum();
totalflux->SetNameTitle("nuisance_flux","NUISANCE Total Flux");
+ TH1D* totalevent = (TH1D*) totalflux->Clone();
+ totalevent->SetNameTitle("nuisance_events", "NUISANCE Total Events");
+ // totalevent->Reset();
+
TH1D* totalxsec = (TH1D*)totalflux->Clone();
totalxsec->SetNameTitle("nuisance_xsec", "NUISANCE Total XSec");
totalxsec->Reset();
- TH1D* totalevent = (TH1D*) totalflux->Clone();
- totalevent->SetNameTitle("nuisance_events", "NUISANCE Total Events");
- totalevent->Reset();
-
// Save the cross-section histograms
for (UInt_t i = 0; i < gOptNuPDGs.size(); i++){
for (UInt_t j = 0; j < gOptTgtPDGs.size(); j++){
int npdg = gOptNuPDGs[i];
TH1D* fluxhist = gOptNuFluxs[i];
int tpdg = gOptTgtPDGs[j];
- double tfrac = gOptTgtMix[tpdg];
+ // Determine the total number of these targets by using the tgt ratios backwards
+ double mixfrac = gOptTgtMix[tpdg];
+
+ int tnucl = tpdg % 10000 / 10; // Total Single Target Nucleons
+ double enucl = double(mixfrac * double(gOptNumberNucleons)); // Expected nucleons from weight fractions
+
+ double tfrac = double(enucl) / double(tnucl); // Scaling factor.
+
// Get Initial state driver
InitialState init_state(tpdg, npdg);
GEVGPool* gpool = mcj_driver->GetConfigPool();
GEVGDriver * evgdriver = gpool->FindDriver(init_state);
const Spline * totxsecspl = evgdriver->XSecSumSpline();
// make new empty xsec hist
TH1D* xsechist = (TH1D*) totalflux->Clone();
xsechist->SetNameTitle((init_state.AsString()+"_xsec").c_str(),init_state.AsString().c_str());
xsechist->Reset();
TH1D* eventhist = (TH1D*) totalflux->Clone();
eventhist->SetNameTitle((init_state.AsString()+"_events").c_str(),init_state.AsString().c_str());
eventhist->Reset();
+ std::cout << "Total Frac for " << init_state.AsString() << " = " << tfrac << std::endl;
+ std::cout << "TNUCL = " << tnucl << std::endl;
+ std::cout << "ENUCL = " << enucl << std::endl;
+ std::cout << "MIXFRAC = " << mixfrac << std::endl;
+ std::cout << "TPDG = " << tpdg << std::endl;
+
// Fill xsec hist
for (int k = 0; k < totalflux->GetNbinsX(); k++){
- double E = fluxhist->GetXaxis()->GetBinCenter(k+1);
- double xsec = totxsecspl->Evaluate( E );
- xsechist->SetBinContent(k+1, xsec);
- eventhist->SetBinContent(k+1, xsec * fluxhist->GetBinContent(k+1));
+ double avgxsec = 0.0;
+ int res = 100;
+ for (int l = 0; l < res; l++){
+ double E = fluxhist->GetXaxis()->GetBinLowEdge(k+1) + (double(l)*fluxhist->GetXaxis()->GetBinWidth(k+1) / double(res));
+ double xsec = totxsecspl->Evaluate( E ) / (1E-38 * genie::units::cm2);
+ avgxsec += xsec;
+ }
+ avgxsec /= double(res);
+ xsechist->SetBinContent(k+1, avgxsec);
+ eventhist->SetBinContent(k+1, avgxsec * fluxhist->GetBinContent(k+1));
}
// Scale totals by target fractions
xsechist->Scale(tfrac);
eventhist->Scale(tfrac);
// Add to totals
totalxsec->Add(xsechist);
- totalevent->Add(eventhist);
+ //totalevent->Add(eventhist);
+ std::cout << "Total XSec Hist for " << init_state.AsString() << " = " << xsechist->Integral("width") << std::endl;
+
// Scale by nnucleons
- xsechist->Scale(1.0 / double(gOptNumberNucleons));
- eventhist->Scale(1.0 / double(gOptNumberNucleons));
+ // xsechist->Scale(1.0 / double(gOptNumberNucleons));
+ // eventhist->Scale(1.0 / double(gOptNumberNucleons));
// Write our new xsec histto file
fluxhist->SetNameTitle((init_state.AsString()+"_flux").c_str(),init_state.AsString().c_str());
fluxhist->Write();
xsechist->Write();
eventhist->Write();
delete xsechist;
delete eventhist;
+ // sleep(10);
}
}
- totalxsec->Scale(double(gOptNumberNucleons));
- totalevent->Scale(double(gOptNumberNucleons));
+ // Determine NUISANCE Style Histograms
+ totalevent->SetNameTitle("nuisance_events", "NUISANCE Total Events");
+ totalevent->Multiply(totalxsec);
+ totalevent->Scale(1.0 / double(gOptNumberNucleons));
+ std::cout << "GOPTNumberNucleons = " << gOptNumberNucleons << std::endl;
+ std::cout << "Inclusive XSec Per Nucleon = " << totalevent->Integral("width") * 1E-38 / totalflux->Integral("width") << std::endl;
+ std::cout << "XSec Hist Integral = " << totalxsec->Integral("width") << std::endl;
+ // sleep(20);
+
+ // totalxsec->Scale(1.0/double(gOptNumberNucleons));
+ // totalevent->Scale(1.0/double(gOptNumberNucleons));
// Create an MC Job Monitor
GMCJMonitor mcjmonitor(gOptRunNu);
mcjmonitor.SetRefreshRate(RunOpt::Instance()->MCJobStatusRefreshRate());
// Generate events / print the GHEP record / add it to the ntuple
int ievent = 0;
while ( ievent < gOptNevents) {
LOG("gevgen", pNOTICE) << " *** Generating event............ " << ievent;
// generate a single event for neutrinos coming from the specified flux
EventRecord * event = mcj_driver->GenerateEvent();
LOG("gevgen", pNOTICE) << "Generated Event GHEP Record: " << *event;
// add event at the output ntuple, refresh the mc job monitor & clean-up
ntpw.AddEventRecord(ievent, event);
mcjmonitor.Update(ievent,event);
ievent++;
delete event;
}
// Save the generated MC events
ntpw.Save();
delete flux_driver;
delete geom_driver;
delete mcj_driver;;
}
//____________________________________________________________________________
GeomAnalyzerI * GeomDriver(void)
{
// create a trivial point geometry with the specified target or target mix
GeomAnalyzerI * geom_driver = new geometry::PointGeomAnalyzer(gOptTgtMix);
return geom_driver;
}
//____________________________________________________________________________
GFluxI * FluxDriver(void)
{
// create & configure one of the generic flux drivers
//
GFluxI * flux_driver = 0;
if(gOptNuEnergyRange<0) flux_driver = MonoEnergeticFluxDriver();
else flux_driver = TH1FluxDriver();
return flux_driver;
}
//____________________________________________________________________________
GFluxI * MonoEnergeticFluxDriver(void)
{
//
//
flux::GMonoEnergeticFlux * flux =
new flux::GMonoEnergeticFlux(gOptNuEnergy, gOptNuPdgCode);
GFluxI * flux_driver = dynamic_cast<GFluxI *>(flux);
return flux_driver;
}
//____________________________________________________________________________
GFluxI * TH1FluxDriver(void)
{
// GeVGEN NUISANCE Can only Read flux files from ROOT
flux::GNUISANCEFlux * flux = new flux::GNUISANCEFlux;
// Initial beam directions
TVector3 bdir (0,0,1);
TVector3 bspot(0,0,0);
flux->SetNuDirection (bdir);
flux->SetBeamSpot (bspot);
flux->SetTransverseRadius (-1);
// Flux inputs
int flux_entries = 100000;
double emin = gOptNuEnergy;
double emax = gOptNuEnergy+gOptNuEnergyRange;
double de = gOptNuEnergyRange;
// check whether the input flux is a file or a functional form
bool input_is_root_file = gOptFlux.find(".root") != string::npos &&
gOptFlux.find(",") != string::npos;
if (!input_is_root_file){
std::cout << "ONLY ROOT INPUTS ALLOWED" << std::endl;
throw;
}
//
// ** extract specified flux histogram from the input root file
//
vector<string> fv = utils::str::Split(gOptFlux,",");
assert(fv.size()>=2);
assert( !gSystem->AccessPathName(fv[0].c_str()) );
// First entry is the root file
LOG("gevgen", pNOTICE) << "Getting input flux from root file: " << fv[0];
TFile * flux_file = new TFile(fv[0].c_str(),"read");
// Later entries are flux inputs
for (UInt_t i = 1; i < fv.size(); i++){
string fluxid = fv[i];
// Parse out the beam pdg and flux name
string flux_with_pdg = fluxid;
string::size_type open_bracket = flux_with_pdg.find("[");
string::size_type close_bracket = flux_with_pdg.find("]");
string::size_type ibeg = 0;
string::size_type iend = open_bracket;
string::size_type jbeg = open_bracket+1;
string::size_type jend = close_bracket-1;
string fluxname = flux_with_pdg.substr(ibeg,iend);
int fluxpdg = atoi(flux_with_pdg.substr(jbeg,jend).c_str());
// Get histogram
LOG("gevgen", pNOTICE) << "Flux name: " << fluxname << " pdg = " << fluxpdg;
TH1D * hst = (TH1D *)flux_file->Get(fluxname.c_str());
assert(hst);
LOG("gevgen", pNOTICE) << hst->GetEntries();
// Copy in the flux histogram from the root file and remove bins outside the emin,emax range
TH1D* spectrum = (TH1D*)hst->Clone();
spectrum->SetNameTitle("spectrum","neutrino_flux");
spectrum->SetDirectory(0);
- for(int ibin = 1; ibin <= hst->GetNbinsX(); ibin++) {
- if(hst->GetBinLowEdge(ibin) + hst->GetBinWidth(ibin) > emax ||
- hst->GetBinLowEdge(ibin) < emin) {
- spectrum->SetBinContent(ibin, 0);
- }
- }
+ // for(int ibin = 1; ibin <= hst->GetNbinsX(); ibin++) {
+ // if(spectrum->GetBinCenter(ibin) > emax ||
+ // spectrum->GetBinCenter(ibin) < emin){
+ // spectrum->SetBinContent(ibin, 0);
+ // }
+ //}
LOG("gevgen", pNOTICE) << spectrum->GetEntries();
// Add energy spectrum to the flux driver.
flux->AddEnergySpectrum( fluxpdg, spectrum );
gOptNuPDGs.push_back(fluxpdg);
gOptNuFluxs.push_back(spectrum);
}
// flux->AddAllFluxes();
// Close inputs
flux_file->Close();
delete flux_file;
// Return flux driver.
GFluxI * flux_driver = dynamic_cast<GFluxI *>(flux);
return flux_driver;
}
///____________________________________________________________________________
void ListTargetIDs(){
// Keep in sync with ConvertTargetIDs
LOG("gevgen", pNOTICE) << "Possible Target IDs: \n"
<< "\n H : " << ConvertTargetIDs("H")
<< "\n C : " << ConvertTargetIDs("C")
<< "\n CH : " << ConvertTargetIDs("CH")
<< "\n CH2 : " << ConvertTargetIDs("CH2")
- << "\n H2O : " << ConvertTargetIDs("H2O");
-
+ << "\n H2O : " << ConvertTargetIDs("H2O")
+ << "\n Fe : " << ConvertTargetIDs("Fe")
+ << "\n Pb : " << ConvertTargetIDs("Pb");
+
}
//____________________________________________________________________________
string ConvertTargetIDs(string id){
if (!id.compare("H")) return "1000010010";
else if (!id.compare("C")) return "1000060120";
else if (!id.compare("CH")) return "13,1000060120[0.9231],1000010010[0.0769]";
else if (!id.compare("CH2")) return "14,1000060120[0.8571],1000010010[0.1429]";
else if (!id.compare("H2O")) return "18,1000080160[0.8888],1000010010[0.1111]";
+ else if (!id.compare("Fe")) return "1000260560";
+ else if (!id.compare("Pb")) return "1000822070";
else return "";
};
///____________________________________________________________________________
void ListFluxIDs(){
// Keep in sync with ConvertTargetIDs
LOG("gevgen", pNOTICE) << "Possible Flux IDs: \n"
<< "\n MINERvA_fhc_numu : " << ConvertTargetIDs("MINERvA_fhc_numu")
<< "\n MINERvA_fhc_numunumubar : " << ConvertTargetIDs("MINERvA_fhc_numunumubar")
<< "\n MINERvA_fhc_nue : " << ConvertTargetIDs("MINERvA_fhc_nue")
<< "\n MINERvA_fhc_nuenuebar : " << ConvertTargetIDs("MINERvA_fhc_nuenuebar")
<< "\n MINERvA_fhc_all : " << ConvertTargetIDs("MINERvA_fhc_all")
<< "\n MINERvA_rhc_numubar : " << ConvertTargetIDs("MINERvA_rhc_numubar")
<< "\n MINERvA_rhc_numubarnumu : " << ConvertTargetIDs("MINERvA_rhc_numubarnumu")
<< "\n MINERvA_rhc_nuebar : " << ConvertTargetIDs("MINERvA_rhc_nuebar")
<< "\n MINERvA_rhc_nuebarnue : " << ConvertTargetIDs("MINERvA_rhc_nuebarnue")
<< "\n MINERvA_rhc_all : " << ConvertTargetIDs("MINERvA_rhc_all");
}
//____________________________________________________________________________
string ConvertFluxIDs(string id){
char * const var = getenv("EXT_FIT");
if (!var) {
std::cout << "Cannot find top level directory! Set the EXT_FIT environmental variable" << std::endl;
exit(-1);
}
string topnuisancedir = string(var);
string fluxfolder = topnuisancedir + "/data/flux/";
- if (!id.compare("MINERvA_fhc_numu")) return "minerva_flux.root,numu_fhc[14]";
- else if (!id.compare("MINERvA_fhc_numunumubar")) return "minerva_flux.root,numu_fhc[14],numubar_fhc[-14]";
- else if (!id.compare("MINERvA_fhc_numu")) return "minerva_flux.root,nue_fhc[12]";
- else if (!id.compare("MINERvA_fhc_nuenuebar")) return "minerva_flux.root,nue_fhc[12],nuebar_fhc[-12]";
- else if (!id.compare("MINERvA_fhc_all")) return "minerva_flux.root,numu_fhc[14],numubar_fhc[-14],nue_fhc[12],nuebar_fhc[-12]";
-
- else if (!id.compare("MINERvA_rhc_numubar")) return "minerva_flux.root,numubar_rhc[-14]";
- else if (!id.compare("MINERvA_rhc_numubarnumu")) return "minerva_flux.root,numubar_rhc[-14],numu_rhc[14]";
- else if (!id.compare("MINERvA_rhc_nuebar")) return "minerva_flux.root,nuebar_rhc[-12]";
- else if (!id.compare("MINERvA_rhc_nuebarnue")) return "minerva_flux.root,nuebar_rhc[-12],nue_rhc[12]";
- else if (!id.compare("MINERvA_rhc_all")) return "minerva_flux.root,numu_rhc[14],numubar_rhc[-14],nue_rhc[12],nuebar_rhc[-12]";
-
-
+ string inputs = "";
+
+ if (!id.compare("MINERvA_fhc_numu")) inputs="minerva_flux.root,numu_fhc[14]";
+ else if (!id.compare("MINERvA_fhc_numunumubar")) inputs="minerva_flux.root,numu_fhc[14],numubar_fhc[-14]";
+ else if (!id.compare("MINERvA_fhc_numu")) inputs="minerva_flux.root,nue_fhc[12]";
+ else if (!id.compare("MINERvA_fhc_nuenuebar")) inputs="minerva_flux.root,nue_fhc[12],nuebar_fhc[-12]";
+ else if (!id.compare("MINERvA_fhc_all")) inputs="minerva_flux.root,numu_fhc[14],numubar_fhc[-14],nue_fhc[12],nuebar_fhc[-12]";
+
+ else if (!id.compare("MINERvA_rhc_numubar")) inputs="minerva_flux.root,numubar_rhc[-14]";
+ else if (!id.compare("MINERvA_rhc_numubarnumu")) inputs="minerva_flux.root,numubar_rhc[-14],numu_rhc[14]";
+ else if (!id.compare("MINERvA_rhc_nuebar")) inputs="minerva_flux.root,nuebar_rhc[-12]";
+ else if (!id.compare("MINERvA_rhc_nuebarnue")) inputs="minerva_flux.root,nuebar_rhc[-12],nue_rhc[12]";
+ else if (!id.compare("MINERvA_rhc_all")) inputs="minerva_flux.root,numu_rhc[14],numubar_rhc[-14],nue_rhc[12],nuebar_rhc[-12]";
else return "";
+ return fluxfolder + inputs;
+
};
//............................................................................
#endif
//____________________________________________________________________________
void GetCommandLineArgs(int argc, char ** argv)
{
LOG("gevgen", pINFO) << "Parsing command line arguments";
// Common run options. Set defaults and read.
RunOpt::Instance()->EnableBareXSecPreCalc(true);
RunOpt::Instance()->ReadFromCommandLine(argc,argv);
// Parse run options for this app
CmdLnArgParser parser(argc,argv);
// help?
bool help = parser.OptionExists('h');
if(help) {
PrintSyntax();
exit(0);
}
// number of events
if( parser.OptionExists('n') ) {
LOG("gevgen", pINFO) << "Reading number of events to generate";
gOptNevents = parser.ArgAsInt('n');
} else {
LOG("gevgen", pINFO)
<< "Unspecified number of events to generate - Using default";
gOptNevents = kDefOptNevents;
}
// run number
if( parser.OptionExists('r') ) {
LOG("gevgen", pINFO) << "Reading MC run number";
gOptRunNu = parser.ArgAsLong('r');
} else {
LOG("gevgen", pINFO) << "Unspecified run number - Using default";
gOptRunNu = kDefOptRunNu;
}
// flux functional form
bool using_flux = false;
if( parser.OptionExists('f') ) {
LOG("gevgen", pINFO) << "Reading flux function";
gOptFlux = parser.ArgAsString('f');
// Convert for known strings.
string fluxid = ConvertFluxIDs(gOptFlux);
if (!fluxid.empty()) gOptFlux = fluxid;
using_flux = true;
}
if(parser.OptionExists('s')) {
LOG("gevgen", pWARN)
<< "-s option no longer available. Please read the revised code documentation";
gAbortingInErr = true;
exit(1);
}
// generate weighted events option (only relevant if using a flux)
gOptWeighted = parser.OptionExists('w');
// neutrino energy
if( parser.OptionExists('e') ) {
LOG("gevgen", pINFO) << "Reading neutrino energy";
string nue = parser.ArgAsString('e');
// is it just a value or a range (comma separated set of values)
if(nue.find(",") != string::npos) {
// split the comma separated list
vector<string> nurange = utils::str::Split(nue, ",");
assert(nurange.size() == 2);
double emin = atof(nurange[0].c_str());
double emax = atof(nurange[1].c_str());
assert(emax>emin && emin>=0);
gOptNuEnergy = emin;
gOptNuEnergyRange = emax-emin;
if(!using_flux) {
LOG("gevgen", pWARN)
<< "No flux was specified but an energy range was input!";
LOG("gevgen", pWARN)
<< "Events will be generated at fixed E = " << gOptNuEnergy << " GeV";
gOptNuEnergyRange = -1;
}
} else {
gOptNuEnergy = atof(nue.c_str());
gOptNuEnergyRange = -1;
}
} else {
LOG("gevgen", pFATAL) << "Unspecified neutrino energy - Exiting";
PrintSyntax();
exit(1);
}
// neutrino PDG code
if( parser.OptionExists('p') && !parser.OptionExists('f')) {
LOG("gevgen", pINFO) << "Reading neutrino PDG code";
gOptNuPdgCode = parser.ArgAsInt('p');
} else if (!parser.OptionExists('p') && !parser.OptionExists('f')){
LOG("gevgen", pFATAL) << "Unspecified neutrino PDG code or Flux Inputs - Exiting";
PrintSyntax();
exit(1);
}
// target mix (their PDG codes with their corresponding weights)
bool using_tgtmix = false;
if( parser.OptionExists('t') ) {
LOG("gevgen", pINFO) << "Reading target mix";
string stgtmix = parser.ArgAsString('t');
gOptTgtMix.clear();
// Check for ID Strings
string tgtids = ConvertTargetIDs(stgtmix);
if (!tgtids.empty()){
stgtmix = tgtids;
}
// Parse Targets
vector<string> tgtmix = utils::str::Split(stgtmix,",");
if(tgtmix.size()==1) {
int pdg = atoi(tgtmix[0].c_str());
double wgt = 1.0;
gOptTgtPDGs.push_back(pdg);
gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
// For single target number of free nucleons
// automatic
- gOptNumberNucleons = pdg % 1000/ 10;
+ gOptNumberNucleons = pdg % 10000/ 10;
} else {
using_tgtmix = true;
vector<string>::const_iterator tgtmix_iter = tgtmix.begin();
// For multiple targets N nucleons must be specified first.
gOptNumberNucleons = atoi((*tgtmix_iter).c_str());
tgtmix_iter++;
// Now remainder is the full target list
for( ; tgtmix_iter != tgtmix.end(); ++tgtmix_iter) {
string tgt_with_wgt = *tgtmix_iter;
string::size_type open_bracket = tgt_with_wgt.find("[");
string::size_type close_bracket = tgt_with_wgt.find("]");
string::size_type ibeg = 0;
string::size_type iend = open_bracket;
string::size_type jbeg = open_bracket+1;
string::size_type jend = close_bracket-1;
int pdg = atoi(tgt_with_wgt.substr(ibeg,iend).c_str());
double wgt = atof(tgt_with_wgt.substr(jbeg,jend).c_str());
LOG("Main", pNOTICE)
<< "Adding to target mix: pdg = " << pdg << ", wgt = " << wgt;
gOptTgtPDGs.push_back(pdg);
gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
}//tgtmix_iter
}//>1
} else {
LOG("gevgen", pFATAL) << "Unspecified target PDG code - Exiting";
PrintSyntax();
exit(1);
}
gOptUsingFluxOrTgtMix = using_flux || using_tgtmix;
// random number seed
if( parser.OptionExists("seed") ) {
LOG("gevgen", pINFO) << "Reading random number seed";
gOptRanSeed = parser.ArgAsLong("seed");
} else {
LOG("gevgen", pINFO) << "Unspecified random number seed - Using default";
gOptRanSeed = -1;
}
// input cross-section file
if( parser.OptionExists("cross-sections") ) {
LOG("gevgen", pINFO) << "Reading cross-section file";
gOptInpXSecFile = parser.ArgAsString("cross-sections");
} else {
LOG("gevgen", pINFO) << "Unspecified cross-section file";
gOptInpXSecFile = "";
}
//
// print-out the command line options
//
LOG("gevgen", pNOTICE)
<< "\n"
<< utils::print::PrintFramedMesg("gevgen job configuration");
LOG("gevgen", pNOTICE)
<< "MC Run Number: " << gOptRunNu;
if(gOptRanSeed != -1) {
LOG("gevgen", pNOTICE)
<< "Random number seed: " << gOptRanSeed;
} else {
LOG("gevgen", pNOTICE)
<< "Random number seed was not set, using default";
}
LOG("gevgen", pNOTICE)
<< "Number of events requested: " << gOptNevents;
if(gOptInpXSecFile.size() > 0) {
LOG("gevgen", pNOTICE)
<< "Using cross-section splines read from: " << gOptInpXSecFile;
} else {
LOG("gevgen", pNOTICE)
<< "No input cross-section spline file";
}
LOG("gevgen", pNOTICE)
<< "Flux: " << gOptFlux;
LOG("gevgen", pNOTICE)
<< "Generate weighted events? " << gOptWeighted;
if(gOptNuEnergyRange>0) {
LOG("gevgen", pNOTICE)
<< "Neutrino energy: ["
<< gOptNuEnergy << ", " << gOptNuEnergy+gOptNuEnergyRange << "]";
} else {
LOG("gevgen", pNOTICE)
<< "Neutrino energy: " << gOptNuEnergy;
}
LOG("gevgen", pNOTICE)
<< "Neutrino code (PDG): " << gOptNuPdgCode;
LOG("gevgen", pNOTICE)
<< "Target code (PDG) & weight fraction (in case of multiple targets): ";
map<int,double>::const_iterator iter;
for(iter = gOptTgtMix.begin(); iter != gOptTgtMix.end(); ++iter) {
int tgtpdgc = iter->first;
double wgt = iter->second;
LOG("gevgen", pNOTICE)
<< " >> " << tgtpdgc << " (weight fraction = " << wgt << ")";
}
LOG("gevgen", pNOTICE) << "\n";
LOG("gevgen", pNOTICE) << *RunOpt::Instance();
}
//____________________________________________________________________________
void PrintSyntax(void)
{
LOG("gevgen", pNOTICE)
<< "\n\n" << "Syntax:" << "\n"
<< "\n gevgen [-h]"
<< "\n [-r run#]"
<< "\n -n nev"
<< "\n -e energy (or energy range) "
<< "\n -p neutrino_pdg"
<< "\n -t target_pdg "
<< "\n [-f flux_description]"
<< "\n [-w]"
<< "\n [--seed random_number_seed]"
<< "\n [--cross-sections xml_file]"
<< "\n [--event-generator-list list_name]"
<< "\n [--message-thresholds xml_file]"
<< "\n [--unphysical-event-mask mask]"
<< "\n [--event-record-print-level level]"
<< "\n [--mc-job-status-refresh-rate rate]"
<< "\n [--cache-file root_file]"
<< "\n\n" ;
ListTargetIDs();
ListFluxIDs();
}
//____________________________________________________________________________
+#endif
diff --git a/cmake/setup.sh.in b/cmake/setup.sh.in
index f658bc2..e4a28cf 100644
--- a/cmake/setup.sh.in
+++ b/cmake/setup.sh.in
@@ -1,120 +1,120 @@
# 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/>.
################################################################################
#!/bin/sh
if ! [[ ":$PATH:" == *":@CMAKE_INSTALL_PREFIX@/bin:"* ]]; then
export PATH=@CMAKE_INSTALL_PREFIX@/bin:$PATH
fi
if ! [[ ":$LD_LIBRARY_PATH:" == *":@CMAKE_INSTALL_PREFIX@/lib:"* ]]; then
export LD_LIBRARY_PATH=@CMAKE_INSTALL_PREFIX@/lib:$LD_LIBRARY_PATH
fi
if [[ ! "${ROOTSYS}" ]]; then
echo "[INFO]: Sourcing ROOT from: @CMAKE_ROOTSYS@"
source "@CMAKE_ROOTSYS@/bin/thisroot.sh"
fi
if [[ "@USE_NEUT@" != "0" ]]; then
echo "[INFO]: Adding NEUT library paths to the environment."
export NEUT_ROOT=@NEUT_ROOT@
export CERN=@CERN@
export CERN_LEVEL=@CERN_LEVEL@
if ! [[ ":$LD_LIBRARY_PATH:" == *":${NEUT_ROOT}/lib/Linux_pc:"* ]]; then
export LD_LIBRARY_PATH=${NEUT_ROOT}/lib/Linux_pc:$LD_LIBRARY_PATH
fi
if ! [[ ":$LD_LIBRARY_PATH:" == *":${NEUT_ROOT}/src/reweight:"* ]]; then
export LD_LIBRARY_PATH=${NEUT_ROOT}/src/reweight:$LD_LIBRARY_PATH
fi
fi
if [[ "@USE_NuWro@" != "0" ]]; then
echo "[INFO]: Adding NuWro library paths to the environment."
export NUWRO="@NUWRO@"
if ! [[ ":$LD_LIBRARY_PATH:" == *":@NUWRO@/build/@CMAKE_SYSTEM_NAME@/lib:"* ]]; then
export LD_LIBRARY_PATH=@NUWRO@/build/@CMAKE_SYSTEM_NAME@/lib:$LD_LIBRARY_PATH
fi
fi
if [[ "@NEED_PYTHIA6@" != "0" ]]; then
echo "[INFO]: Adding PYTHIA6 library paths to the environment."
export PYTHIA6="@PYTHIA6@"
if ! [[ ":$LD_LIBRARY_PATH:" == *":@PYTHIA6@:"* ]]; then
export LD_LIBRARY_PATH=@PYTHIA6@:$LD_LIBRARY_PATH
fi
fi
if [[ "@USE_GENIE@" != "0" ]]; then
echo "[INFO]: Adding GENIE paths to the environment."
export GENIE="@GENIE@"
export LHAPDF_LIB="@LHAPDF_LIB@"
export LHAPDF_INC="@LHAPDF_INC@"
export LIBXML2_LIB="@LIBXML2_LIB@"
export LIBXML2_INC="@LIBXML2_INC@"
export LOG4CPP_LIB="@LOG4CPP_LIB@"
export LOG4CPP_INC="@LOG4CPP_INC@"
export LHAPATH="@LHAPATH@"
if ! [[ ":$LD_LIBRARY_PATH:" == *":@GENIE@/lib:"* ]]; then
export LD_LIBRARY_PATH=@GENIE@/lib:$LD_LIBRARY_PATH
fi
if ! [[ ":$LD_LIBRARY_PATH:" == *":@LHAPDF_LIB@:"* ]]; then
export LD_LIBRARY_PATH=@LHAPDF_LIB@:$LD_LIBRARY_PATH
fi
if ! [[ ":$LD_LIBRARY_PATH:" == *":@LIBXML2_LIB@:"* ]]; then
export LD_LIBRARY_PATH=@LIBXML2_LIB@:$LD_LIBRARY_PATH
fi
if ! [[ ":$LD_LIBRARY_PATH:" == *":@LOG4CPP_LIB@:"* ]]; then
export LD_LIBRARY_PATH=@LOG4CPP_LIB@:$LD_LIBRARY_PATH
fi
fi
if [[ "@USE_NIWG@" != "0" ]]; then
echo "[INFO]: Adding NIWG paths to the environment."
export NIWG=@NIWG@
export NIWGREWEIGHT_INPUTS=@NIWG@/inputs
if ! [[ ":$LD_LIBRARY_PATH:" == *":@NIWG@:"* ]]; then
export LD_LIBRARY_PATH=${NIWG}:${LD_LIBRARY_PATH}
fi
fi
if [[ "@USE_T2K@" != "0" ]]; then
echo "[INFO]: Adding T2K paths to the environment."
export T2KREWEIGHT=@T2KREWEIGHT@
if ! [[ ":$LD_LIBRARY_PATH:" == *":@T2KREWEIGHT@/lib:"* ]]; then
export LD_LIBRARY_PATH=${T2KREWEIGHT}/lib:${LD_LIBRARY_PATH}
fi
fi
-if [[ "@USE_GiBUU@" != "0" ]]; then
+if [[ "@BUILD_GiBUU@" != "0" ]]; then
echo "[INFO]: Sourcing GiBUU tools."
source @CMAKE_BINARY_DIR@/GiBUUTools/src/GiBUUTools-build/Linux/setup.sh
fi
export EXT_FIT="@CMAKE_SOURCE_DIR@"
diff --git a/src/FitBase/CMakeLists.txt b/src/FitBase/CMakeLists.txt
index 7d1a1eb..cd7b95d 100644
--- a/src/FitBase/CMakeLists.txt
+++ b/src/FitBase/CMakeLists.txt
@@ -1,105 +1,109 @@
# 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/>.
################################################################################
set(IMPLFILES
ParamPull.cxx
BaseFitEvt.cxx
FitParticle.cxx
Measurement1D.cxx
EventManager.cxx
Measurement2D.cxx
FitEvent.cxx
JointMeas1D.cxx
MeasurementBase.cxx
GeneratorUtils.cxx
StdHepEvt.cxx
TemplateMeas1D.cxx
InputUtils.cxx
NEUTInputHandler.cxx
GENIEInputHandler.cxx
NuWroInputHandler.cxx
GIBUUInputHandler.cxx
NUANCEInputHandler.cxx
NuanceEvent.cxx
SampleSettings.cxx
FitEventInputHandler.cxx
SplineInputHandler.cxx
MeasurementVariableBox.cxx
MeasurementVariableBox2D.cxx
MeasurementVariableBox1D.cxx
CustomVariableBoxes.h
GeneratorInfoBase.h
InputTypes.h
GNUISANCEFlux.cxx
GNUISANCEMCJDriver.cxx
)
set(HEADERFILES
ParamPull.h
NuanceEvent.h
BaseFitEvt.h FitEvent.h JointMeas1D.h Measurement2D.h
EventManager.h FitParticle.h MeasurementBase.h Measurement1D.h GeneratorUtils.h StdHepEvt.h TemplateMeas1D.h
InputUtils.h
NEUTInputHandler.h
GENIEInputHandler.h
NuWroInputHandler.h
InputHandler2.h
GIBUUInputHandler.h
NUANCEInputHandler.h
SampleSettings.h
FitEventInputHandler.h
SplineInputHandler.h
MeasurementVariableBox.h
MeasurementVariableBox2D.h
MeasurementVariableBox1D.h
CustomVariableBoxes.h
GeneratorInfoBase.h
GNUISANCEFlux.h
GNUISANCEMCJDriver.h
)
set(LIBNAME FitBase)
if(CMAKE_BUILD_TYPE MATCHES DEBUG)
add_library(${LIBNAME} STATIC ${IMPLFILES})
else(CMAKE_BUILD_TYPE MATCHES RELEASE)
add_library(${LIBNAME} SHARED ${IMPLFILES})
endif()
include_directories(${RWENGINE_INCLUDE_DIRECTORIES})
-include_directories(${GENIE_INCLUDES}/Apps)
-include_directories(${GENIE_INCLUDES}/FluxDrivers)
-include_directories(${GENIE_INCLUDES}/EVGDrivers)
+
+if (USE_GENIE)
+ include_directories(${GENIE_INCLUDES}/Apps)
+ include_directories(${GENIE_INCLUDES}/FluxDrivers)
+ include_directories(${GENIE_INCLUDES}/EVGDrivers)
+endif()
+
include_directories(${CMAKE_SOURCE_DIR}/src/Reweight)
include_directories(${CMAKE_SOURCE_DIR}/src/Splines)
include_directories(${CMAKE_SOURCE_DIR}/src/Utils)
include_directories(${CMAKE_CURRENT_SOURCE_DIR})
set_target_properties(${LIBNAME} PROPERTIES VERSION
"${ExtFit_VERSION_MAJOR}.${ExtFit_VERSION_MINOR}.${ExtFit_VERSION_REVISION}")
set_target_properties(${LIBNAME} PROPERTIES LINK_FLAGS ${ROOT_LD_FLAGS})
if(DEFINED PROJECTWIDE_EXTRA_DEPENDENCIES)
add_dependencies(${LIBNAME} ${PROJECTWIDE_EXTRA_DEPENDENCIES})
endif()
install(TARGETS ${LIBNAME} DESTINATION lib)
#Can uncomment this to install the headers... but is it really neccessary?
install(FILES ${HEADERFILES} DESTINATION include)
set(MODULETargets ${MODULETargets} ${LIBNAME} PARENT_SCOPE)
diff --git a/src/FitBase/GNUISANCEFlux.cxx b/src/FitBase/GNUISANCEFlux.cxx
index 4a0103b..ac65722 100644
--- a/src/FitBase/GNUISANCEFlux.cxx
+++ b/src/FitBase/GNUISANCEFlux.cxx
@@ -1,302 +1,304 @@
+#ifdef __GENIE_ENABLED__
//____________________________________________________________________________
/*
Copyright (c) 2003-2016, GENIE Neutrino MC Generator Collaboration
For the full text of the license visit http://copyright.genie-mc.org
or see $GENIE/LICENSE
Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
University of Liverpool & STFC Rutherford Appleton Lab - July 04, 2005
For the class documentation see the corresponding header file.
Important revisions after version 2.0.0 :
@ Feb 22, 2011 - JD
Implemented dummy versions of the new GFluxI::Clear, GFluxI::Index and
GFluxI::GenerateWeighted methods needed for pre-generation of flux
interaction probabilities in GMCJDriver.
*/
//____________________________________________________________________________
#include <cassert>
#include <algorithm>
#include <TH1D.h>
#include <TF1.h>
#include <TVector3.h>
#include "Conventions/Constants.h"
#include "GNUISANCEFlux.h"
#include "Messenger/Messenger.h"
#include "Numerical/RandomGen.h"
#include "PDG/PDGCodeList.h"
#include "Utils/PrintUtils.h"
#include "FluxDrivers/GFluxDriverFactory.h"
FLUXDRIVERREG4(genie,flux,GNUISANCEFlux,genie::flux::GNUISANCEFlux)
using namespace genie;
using namespace genie::constants;
using namespace genie::flux;
//____________________________________________________________________________
GNUISANCEFlux::GNUISANCEFlux()
{
this->Initialize();
}
//___________________________________________________________________________
GNUISANCEFlux::~GNUISANCEFlux()
{
this->CleanUp();
}
//___________________________________________________________________________
bool GNUISANCEFlux::GenerateNext(void)
{
//-- Reset previously generated neutrino code / 4-p / 4-x
this->ResetSelection();
//-- Generate an energy from the 'combined' spectrum histogram
// and compute the momentum vector
if (!fTotSpectrum) AddAllFluxes();
std::cout << "fTotSpectrum = " << fTotSpectrum << std::endl;
std::cout << "TotalSpectrum Integral = " << fTotSpectrum->Integral() << std::endl;
double Ev = (double) fTotSpectrum->GetRandom();
TVector3 p3(*fDirVec); // momentum along the neutrino direction
p3.SetMag(Ev); // with |p|=Ev
fgP4.SetPxPyPzE(p3.Px(), p3.Py(), p3.Pz(), Ev);
//-- Select a neutrino species from the flux fractions at the
// selected energy
fgPdgC = (*fPdgCList)[this->SelectNeutrino(Ev)];
//-- Compute neutrino 4-x
if(fRt <= 0) {
fgX4.SetXYZT(0.,0.,0.,0.);
}
else {
// Create a vector (vec) that points to a random position at a disk
// of radius Rt passing through the origin, perpendicular to the
// input direction.
TVector3 vec0(*fDirVec); // vector along the input direction
TVector3 vec = vec0.Orthogonal(); // orthogonal vector
double psi = this->GeneratePhi(); // rndm angle [0,2pi]
double Rt = this->GenerateRt(); // rndm R [0,Rtransverse]
vec.Rotate(psi,vec0); // rotate around original vector
vec.SetMag(Rt); // set new norm
// Set the neutrino position as beam_spot + vec
double x = fBeamSpot->X() + vec.X();
double y = fBeamSpot->Y() + vec.Y();
double z = fBeamSpot->Z() + vec.Z();
fgX4.SetXYZT(x,y,z,0.);
}
LOG("Flux", pINFO) << "Generated neutrino pdg-code: " << fgPdgC;
LOG("Flux", pINFO)
<< "Generated neutrino p4: " << utils::print::P4AsShortString(&fgP4);
LOG("Flux", pINFO)
<< "Generated neutrino x4: " << utils::print::X4AsString(&fgX4);
return true;
}
//___________________________________________________________________________
void GNUISANCEFlux::Clear(Option_t * opt)
{
// Dummy clear method needed to conform to GFluxI interface
//
LOG("Flux", pERROR) <<
"No Clear(Option_t * opt) method implemented for opt: "<< opt;
}
//___________________________________________________________________________
void GNUISANCEFlux::GenerateWeighted(bool gen_weighted)
{
// Dummy implementation needed to conform to GFluxI interface
//
LOG("Flux", pERROR) <<
"No GenerateWeighted(bool gen_weighted) method implemented for " <<
"gen_weighted: " << gen_weighted;
}
//___________________________________________________________________________
void GNUISANCEFlux::Initialize(void)
{
LOG("Flux", pNOTICE) << "Initializing GNUISANCEFlux driver";
fMaxEv = 0;
fPdgCList = new PDGCodeList;
fTotSpectrum = 0;
fDirVec = 0;
fBeamSpot = 0;
fRt =-1;
fRtDep = 0;
this->ResetSelection();
this->SetRtDependence("x");
//eg, other example: this->SetRtDependence("pow(x,2)");
}
//___________________________________________________________________________
void GNUISANCEFlux::ResetSelection(void)
{
// initializing running neutrino pdg-code, 4-position, 4-momentum
fgPdgC = 0;
fgP4.SetPxPyPzE (0.,0.,0.,0.);
fgX4.SetXYZT (0.,0.,0.,0.);
}
//___________________________________________________________________________
void GNUISANCEFlux::CleanUp(void)
{
LOG("Flux", pNOTICE) << "Cleaning up...";
if (fDirVec ) delete fDirVec;
if (fBeamSpot ) delete fBeamSpot;
if (fPdgCList ) delete fPdgCList;
// if (fTotSpectrum) delete fTotSpectrum;
if (fRtDep ) delete fRtDep;
unsigned int nspectra = fSpectrum.size();
for(unsigned int i = 0; i < nspectra; i++) {
TH1D * spectrum = fSpectrum[i];
delete spectrum;
spectrum = 0;
}
}
//___________________________________________________________________________
void GNUISANCEFlux::SetNuDirection(const TVector3 & direction)
{
if(fDirVec) delete fDirVec;
fDirVec = new TVector3(direction);
}
//___________________________________________________________________________
void GNUISANCEFlux::SetBeamSpot(const TVector3 & spot)
{
if(fBeamSpot) delete fBeamSpot;
fBeamSpot = new TVector3(spot);
}
//___________________________________________________________________________
void GNUISANCEFlux::SetTransverseRadius(double Rt)
{
LOG ("Flux", pNOTICE) << "Setting R[transverse] = " << Rt;
fRt = Rt;
if(fRtDep) fRtDep->SetRange(0,Rt);
}
//___________________________________________________________________________
void GNUISANCEFlux::AddEnergySpectrum(int nu_pdgc, TH1D * spectrum)
{
LOG("Flux", pNOTICE) << "Adding flux spectrum for pdg = " << nu_pdgc;
fPdgCList->push_back(nu_pdgc);
bool accepted = (count(fPdgCList->begin(),fPdgCList->end(),nu_pdgc) == 1);
if(!accepted) {
LOG ("Flux", pWARN)
<< "The pdg-code isn't recognized and the spectrum was ignored";
} else {
fSpectrum.push_back(spectrum);
int nb = spectrum->GetNbinsX();
Axis_t max = spectrum->GetBinLowEdge(nb)+spectrum->GetBinWidth(nb);
fMaxEv = TMath::Max(fMaxEv, (double)max);
LOG("Flux", pNOTICE)
<< "Updating maximum energy of flux particles to: " << fMaxEv;
// this->AddAllFluxes(); // update combined flux
}
}
//___________________________________________________________________________
void GNUISANCEFlux::SetRtDependence(string rdep)
{
// Set the (functional form of) Rt dependence as string, eg "x*x+sin(x)"
// You do not need to set this method. The default behaviour is to generate
// flux neutrinos uniformly over the area of the cylinder's cross section.
if(fRtDep) delete fRtDep;
fRtDep = new TF1("rdep", rdep.c_str(), 0,fRt);
}
//___________________________________________________________________________
void GNUISANCEFlux::AddAllFluxes(void)
{
LOG("Flux", pNOTICE) << "Computing combined flux";
//if(fTotSpectrum) delete fTotSpectrum;
vector<TH1D *>::const_iterator spectrum_iter;
unsigned int inu=0;
for(spectrum_iter = fSpectrum.begin();
spectrum_iter != fSpectrum.end(); ++spectrum_iter) {
TH1D * spectrum = *spectrum_iter;
if(inu==0) { fTotSpectrum = new TH1D(*spectrum); }
else { fTotSpectrum->Add(spectrum); }
LOG("Flux",pNOTICE) << "Added spectrum ";
inu++;
}
LOG("Flux", pNOTICE) << "Flux Pointer = " << fTotSpectrum;
LOG("Flux", pNOTICE) << "Total flux integral = " << fTotSpectrum->Integral();
}
//___________________________________________________________________________
int GNUISANCEFlux::SelectNeutrino(double Ev)
{
const unsigned int n = fPdgCList->size();
double fraction[n];
vector<TH1D *>::const_iterator spectrum_iter;
unsigned int inu=0;
for(spectrum_iter = fSpectrum.begin();
spectrum_iter != fSpectrum.end(); ++spectrum_iter) {
TH1D * spectrum = *spectrum_iter;
fraction[inu++] = spectrum->GetBinContent(spectrum->FindBin(Ev));
}
double sum = 0;
for(inu = 0; inu < n; inu++) {
sum += fraction[inu];
fraction[inu] = sum;
LOG("Flux", pDEBUG) << "SUM-FRACTION(0->" << inu <<") = " << sum;
}
RandomGen * rnd = RandomGen::Instance();
double R = sum * rnd->RndFlux().Rndm();
LOG("Flux", pDEBUG) << "R e [0,SUM] = " << R;
for(inu = 0; inu < n; inu++) {if ( R < fraction[inu] ) return inu;}
LOG("Flux", pERROR) << "Could not select a neutrino species";
assert(false);
return -1;
}
//___________________________________________________________________________
double GNUISANCEFlux::GeneratePhi(void) const
{
RandomGen * rnd = RandomGen::Instance();
double phi = 2.*kPi * rnd->RndFlux().Rndm(); // [0,2pi]
return phi;
}
//___________________________________________________________________________
double GNUISANCEFlux::GenerateRt(void) const
{
double Rt = fRtDep->GetRandom(); // rndm R [0,Rtransverse]
return Rt;
}
//___________________________________________________________________________
TH1D* GNUISANCEFlux::GetTotalSpectrum(void){
// if (fTotSpectrum) delete fTotSpectrum;
this->AddAllFluxes();
return fTotSpectrum;
}
+#endif
diff --git a/src/FitBase/GNUISANCEFlux.h b/src/FitBase/GNUISANCEFlux.h
index 0b6c3f7..ac34b2b 100644
--- a/src/FitBase/GNUISANCEFlux.h
+++ b/src/FitBase/GNUISANCEFlux.h
@@ -1,101 +1,103 @@
+#ifdef __GENIE_ENABLED__
//____________________________________________________________________________
/*!
\class genie::flux::GNUISANCEFlux
\brief A generic GENIE flux driver.
Generates a 'cylindrical' neutrino beam along the input direction,
with the input transverse radius and centered at the input position.
The energies are generated from the input energy spectrum (TH1D).
Multiple neutrino species can be generated (you will need to supply
an energy spectrum for each).
\author Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
University of Liverpool & STFC Rutherford Appleton Lab
\created July 4, 2005
\cpright Copyright (c) 2003-2016, GENIE Neutrino MC Generator Collaboration
For the full text of the license visit http://copyright.genie-mc.org
or see $GENIE/LICENSE
*/
//____________________________________________________________________________
#ifndef _G_NUISANCE_FLUX_H_
#define _G_NUISANCE_FLUX_H_
#include <string>
#include <vector>
#include <TLorentzVector.h>
#include "EVGDrivers/GFluxI.h"
class TH1D;
class TF1;
class TVector3;
using std::string;
using std::vector;
namespace genie {
namespace flux {
class GNUISANCEFlux: public GFluxI {
public :
GNUISANCEFlux();
~GNUISANCEFlux();
// methods specific to this flux object
void SetNuDirection (const TVector3 & direction);
void SetBeamSpot (const TVector3 & spot);
void SetTransverseRadius (double Rt);
void AddEnergySpectrum (int nu_pdgc, TH1D * spectrum);
void SetRtDependence (string rdep);
// methods implementing the GENIE GFluxI interface
const PDGCodeList & FluxParticles (void) { return *fPdgCList; }
double MaxEnergy (void) { return fMaxEv; }
bool GenerateNext (void);
int PdgCode (void) { return fgPdgC; }
double Weight (void) { return 1.0; }
const TLorentzVector & Momentum (void) { return fgP4; }
const TLorentzVector & Position (void) { return fgX4; }
bool End (void) { return false; }
long int Index (void) { return -1; }
void Clear (Option_t * opt);
void GenerateWeighted (bool gen_weighted);
TH1D* GetTotalSpectrum();
inline vector<TH1D*> GetSpectrum(){ return fSpectrum; };
void AddAllFluxes (void);
private:
// private methods
void Initialize (void);
void CleanUp (void);
void ResetSelection (void);
int SelectNeutrino (double Ev);
double GeneratePhi (void) const;
double GenerateRt (void) const;
// private data members
double fMaxEv; ///< maximum energy
PDGCodeList * fPdgCList; ///< list of neutrino pdg-codes
int fgPdgC; ///< running generated nu pdg-code
TLorentzVector fgP4; ///< running generated nu 4-momentum
TLorentzVector fgX4; ///< running generated nu 4-position
vector<TH1D *> fSpectrum; ///< flux = f(Ev), 1/neutrino species
TH1D * fTotSpectrum; ///< combined flux = f(Ev)
TVector3 * fDirVec; ///< neutrino direction
TVector3 * fBeamSpot; ///< beam spot position
double fRt; ///< transverse size of neutrino beam
TF1 * fRtDep; ///< transverse radius dependence
};
} // flux namespace
} // genie namespace
#endif // _G_TH1_CYLICDRICAL_FLUX_H_
+#endif
diff --git a/src/FitBase/GNUISANCEMCJDriver.cxx b/src/FitBase/GNUISANCEMCJDriver.cxx
index 3b667b8..fe3a1fe 100644
--- a/src/FitBase/GNUISANCEMCJDriver.cxx
+++ b/src/FitBase/GNUISANCEMCJDriver.cxx
@@ -1,1370 +1,1372 @@
+#ifdef __GENIE_ENABLED__
//____________________________________________________________________________
/*
Copyright (c) 2003-2016, GENIE Neutrino MC Generator Collaboration
For the full text of the license visit http://copyright.genie-mc.org
or see $GENIE/LICENSE
Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
University of Liverpool & STFC Rutherford Appleton Lab
For the class documentation see the corresponding header file.
Important revisions after version 2.0.0 :
@ Feb 08, 2008 - CA
Modified the global probability scale to be the maximum amongst the maximum
interaction probabilities for each neutrino (rather than the sum of maximum
probabilities). The modified probability scale still gives unbiased event
generation & reduces the 'no-interaction' probability.
@ Feb 14, 2008 - CA
Significant speed improvements - Most of the rejected flux neutrinos, are
rejected before having them propagated through the detector geometry
(flux neutrinos are pre-selected using the maximum path-lengths / for most
realistic fluxes -high energy tail, low energy peak- most of the selection
inefficiency is caused not because the path-lengths are not close the max
possible ones, but because the energy is not close to the max possible one).
Also some speed improvement was gained by properly using the total cross
section splines (before, han repeatedly summing-up the
The driver does not assert that _each_ flux neutrino generation & geometry
navigation will be succesfull - In the rare event that this may happen, it
prints an err mesg and tries again. In next revision I will limit the
number of successive trials something may go wrong to prevent the driver
from hanging in truly problematic cases.
The driver was adapted to handle flux drivers that -at some point- may stop
generating more flux neutrinos (eg because they read flux neutrinos by
looping over a beam simulation ntuple and they reached its last entry).
Code was appropriately restructured and some methods have been renamed.
@ Feb 29, 2008 - CA
Modified the InteractionProbability() to calculate absolute interaction
probabilities. Added NFluxNeutrinos() and GlobProbScale() to get the
number of neutrinos thrown by the flux driver towards the geometry and
the global interaction probability scale so as to be able to calculate
event sample normalization factors.
@ Jan 15, 2009 - CA
Stopped GMCJDriver from initializing the unphysical event mask so as not
to overwrite the values that each GEVGDriver obtains from the environment.
@ Jan 16, 2009 - CA
Added methods to return pointers to the flux and geometry drivers.
@ Mar 11, 2009 - CA
In GenerateEvent1Try() handle failure to generate kinematics. Added sanity
check on the no interaction probability.
@ Mar 04, 2010 - CA
Remove unused FilterUnphysical(TBits) method. Now set exclusively via the
GUNPHYSMASK env.var.
@ Apr 15, 2010 - CA
Fix unit error in ComputeEventProbability() - Reported by Corey Reed.
The probability stored at the output event was wrong but this doesn't
affect any of the existing applications as this number wasn't actually
used anywhere.
@ Dec 07, 2010 - CA
Don't use a fixed bin size in ComputeProbScales() as this was causing
errors for low energy applications. Addresses a problem reported by
Joachim Kopp.
@ Feb 22, 2011 - JD
Added a number of new methods to allow pre-calculation of exact flux
interaction probabilities for a given set of flux neutrinos from the
flux driver. See the comments for the new LoadFluxProbabilities,
SaveFluxProbabilities, PreCalcFluxProbabilities and PreSelectEvents
methods for details. Using these methods mean that there is no need
to generate maximum path lengths as instead use the exact interaction
probabilities to pre-select. This can result in very significant speed
increases (between factor of 5 and ~300) for event generation over complex
detector geometries and with realistic flux drivers. See
src/support/t2k/EvGen/gT2KEvGen.cxx for an example of how to use.
@ Mar, 7, 2011 - JD
Store sum totals of the flux interaction probabilities for various neutrino
type in a map relating pdg code to total interaction probability. Also add
public getter method so that this can be used in applications to work out
expected event rates. See gT2KEvGen.cxx for an example of how to do this.
Also save the PDG code for each entry in the flux interaction probabilities
tree.
@ Mar, 11, 2011 - JD
Set the directory of fFluxIntTree to the output file fFluxIntProbFile if
saving it later. This is so that it is incrementally saved and fixes bug
where getting std::bad_alloc when trying to Write large trees
fFluxIntProbFile.
@ Jan 31, 2013 - CA
Added SetEventGeneratorList(string listname). $GEVGL var no longer in use.
@ Feb 01, 2013 - CA
The GUNPHYSMASK env. var is no longer used. Added SetUnphysEventMask(const
TBits &). Input is propagated accordingly.
@ Feb 06, 2013 - CA
Fix small problem introduced with recent changes.
In PopulateEventGenDriverPool() calls to GEVGDriver::SetEventGeneratorList()
and GEVGDriver::Configure() were reversed. Problem reported by W.Huelsnitz.
@ July 15, 2014 - HG
Incorporated code provided by Jason Koskinen - IceCube
Modified ComputeProbScales to evalulate the cross sections at both the high
and low edges of the energy bin when calculating the max interaction
probability.
*/
//____________________________________________________________________________
#include <cassert>
#include <TVector3.h>
#include <TSystem.h>
#include <TStopwatch.h>
#include "Algorithm/AlgConfigPool.h"
#include "Conventions/GBuild.h"
#include "Conventions/Constants.h"
#include "Conventions/Units.h"
#include "Conventions/Controls.h"
#include "EVGCore/EventRecord.h"
#include "EVGDrivers/GMCJDriver.h"
#include "GNUISANCEMCJDriver.h"
#include "EVGDrivers/GEVGDriver.h"
#include "EVGDrivers/GEVGPool.h"
#include "EVGDrivers/GFluxI.h"
#include "EVGDrivers/GeomAnalyzerI.h"
#include "GHEP/GHepFlags.h"
#include "GHEP/GHepParticle.h"
#include "Interaction/InitialState.h"
#include "Messenger/Messenger.h"
#include "Numerical/RandomGen.h"
#include "Numerical/Spline.h"
#include "PDG/PDGUtils.h"
#include "Utils/PrintUtils.h"
#include "Utils/XSecSplineList.h"
#include "Conventions/Constants.h"
using namespace genie;
using namespace genie::constants;
//____________________________________________________________________________
GNUISANCEMCJDriver::GNUISANCEMCJDriver()
{
this->InitJob();
}
//___________________________________________________________________________
GNUISANCEMCJDriver::~GNUISANCEMCJDriver()
{
if(fUnphysEventMask) delete fUnphysEventMask;
if (fGPool) delete fGPool;
map<int,TH1D*>::iterator pmax_iter = fPmax.begin();
for( ; pmax_iter != fPmax.end(); ++pmax_iter) {
TH1D * pmax = pmax_iter->second;
if(pmax) {
delete pmax; pmax = 0;
}
}
fPmax.clear();
if(fFluxIntTree) delete fFluxIntTree;
if(fFluxIntProbFile) delete fFluxIntProbFile;
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::SetEventGeneratorList(string listname)
{
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Setting event generator list: " << listname;
fEventGenList = listname;
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::SetUnphysEventMask(const TBits & mask)
{
*fUnphysEventMask = mask;
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Setting unphysical event mask (bits: " << GHepFlags::NFlags() - 1
<< " -> 0) : " << *fUnphysEventMask;
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::UseFluxDriver(GFluxI * flux_driver)
{
fFluxDriver = flux_driver;
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::UseGeomAnalyzer(GeomAnalyzerI * geom_analyzer)
{
fGeomAnalyzer = geom_analyzer;
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::UseSplines(bool useLogE)
{
fUseSplines = true;
fUseLogE = useLogE;
}
//___________________________________________________________________________
bool GNUISANCEMCJDriver::UseMaxPathLengths(string xml_filename)
{
// If you supply the maximum path lengths for all materials in your geometry
// (eg for ROOT/GEANT geometries they can be computed running GENIE's gmxpl
// application, see $GENIE/src/stdapp/gMaxPathLengths.cxx ) you can speed up
// the driver init phase by quite a bit (especially for complex geometries).
fMaxPlXmlFilename = xml_filename;
bool is_accessible = !(gSystem->AccessPathName(fMaxPlXmlFilename.c_str()));
if ( is_accessible ) fUseExtMaxPl = true;
else {
fUseExtMaxPl = false;
LOG("GNUISANCEMCJDriver", pWARN)
<< "UseMaxPathLengths could not find file: \"" << xml_filename << "\"";
}
return fUseExtMaxPl;
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::KeepOnThrowingFluxNeutrinos(bool keep_on)
{
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Keep on throwing flux neutrinos till one interacts? : "
<< utils::print::BoolAsYNString(keep_on);
fKeepThrowingFluxNu = keep_on;
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::ForceSingleProbScale()
{
// Use a single probability scale. That generates unweighted events.
// (Note that generating unweighted event kinematics is a different thing)
//
fGenerateUnweighted = true;
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "GNUISANCEMCJDriver will generate un-weighted events. "
<< "Note: That does not force unweighted event kinematics!";
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::PreSelectEvents(bool preselect)
{
// Set whether to pre-select events based on a max-path lengths file. This
// should be turned off if using pre-generated interaction probabilities
// calculated from a given flux file.
fPreSelect = preselect;
}
//___________________________________________________________________________
bool GNUISANCEMCJDriver::PreCalcFluxProbabilities(void)
{
// Loop over complete set of flux entries satisfying input config options
// (such as neutrino type) and save the interaction probability in a tree
// relating flux index (entry number in input flux tree) to interaction
// probability. If a pre-generated flux interaction probability tree has
// already been loaded then just returns true. Also save tree to a TFile
// for use in later jobs if flag is set
//
bool success = true;
bool save_to_file = fFluxIntProbFile == 0 && fFluxIntFileName.size()>0;
// Clear map storing sum(fBrFluxWeight*fBrFluxIntProb) for each neutrino pdg
fSumFluxIntProbs.clear();
// check if already loaded flux interaction probs using LoadFluxProbTree
if(fFluxIntTree){
LOG("GNUISANCEMCJDriver", pNOTICE) <<
"Skipping pre-generation of flux interaction probabilities - "<<
"using pre-generated file";
success = true;
}
// otherwise create them on the fly now
else {
if(save_to_file){
fFluxIntProbFile = new TFile(fFluxIntFileName.c_str(), "CREATE");
if(fFluxIntProbFile->IsZombie()){
LOG("GNUISANCEMCJDriver", pFATAL) << "Cannot overwrite an existing file. Exiting!";
exit(1);
}
}
// Create the tree to store flux probs
fFluxIntTree = new TTree(fFluxIntTreeName.c_str(),
"Tree storing pre-calculated flux interaction probs");
fFluxIntTree->Branch("FluxIndex", &fBrFluxIndex, "FluxIndex/I");
fFluxIntTree->Branch("FluxIntProb", &fBrFluxIntProb, "FluxIntProb/D");
fFluxIntTree->Branch("FluxEnu", &fBrFluxEnu, "FluxEnu/D");
fFluxIntTree->Branch("FluxWeight", &fBrFluxWeight, "FluxWeight/D");
fFluxIntTree->Branch("FluxPDG", &fBrFluxPDG, "FluxPDG/I");
// Associate to file otherwise get std::bad_alloc when writing large trees
if(save_to_file) fFluxIntTree->SetDirectory(fFluxIntProbFile);
fFluxDriver->GenerateWeighted(true);
fGlobPmax = 1.0; // Force ComputeInteractionProbabilities to return absolute value
// Loop over flux entries and calculate interaction probabilities
TStopwatch stopwatch;
stopwatch.Start();
long int first_index = -1;
bool first_loop = true;
// loop until at end of flux ntuple
while(fFluxDriver->End() == false){
// get the next flux neutrino
bool gotnext = fFluxDriver->GenerateNext();
if(!gotnext){
LOG("GNUISANCEMCJDriver", pWARN) << "*** Couldn't generate next flux ray! ";
continue;
}
// stop if completed a full cycle (this check is necessary as fluxdriver
// may be set to loop over more than one cycle before reaching end)
bool already_been_here = first_loop ? false : first_index == fFluxDriver->Index();
if(already_been_here) break;
// compute the path lengths for current flux neutrino
if(this->ComputePathLengths() == false){ success = false; break;}
// compute and store the interaction probability
double psum = this->ComputeInteractionProbabilities(false /*Based on actual PLs*/);
assert(psum+controls::kASmallNum > 0.);
fBrFluxIntProb = psum;
fBrFluxIndex = fFluxDriver->Index();
fBrFluxEnu = fFluxDriver->Momentum().E();
fBrFluxWeight = fFluxDriver->Weight();
fBrFluxPDG = fFluxDriver->PdgCode();
fFluxIntTree->Fill();
// store the first index so know when have cycled exactly once
if(first_loop){
first_index = fFluxDriver->Index();
first_loop = false;
}
} // flux loop
stopwatch.Stop();
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Finished pre-calculating flux interaction probabilities. "
<< "Total CPU time to process "<< fFluxIntTree->GetEntries()
<< " entries: "<< stopwatch.CpuTime();
// reset the flux driver so can be used at next stage. N.B. This
// should also reset flux driver to throw de-weighted flux neutrinos
fFluxDriver->Clear("CycleHistory");
}
// If successfully calculated/loaded interaction probabilities then set global
// probability scale and, if requested, save tree to output file
if(success){
fGlobPmax = 0.0;
double safety_factor = 1.01;
for(int i = 0; i< fFluxIntTree->GetEntries(); i++){
fFluxIntTree->GetEntry(i);
// Check have non-negative probabilities
assert(fBrFluxIntProb+controls::kASmallNum > 0.0);
assert(fBrFluxWeight+controls::kASmallNum > 0.0);
// Update the global maximum
fGlobPmax = TMath::Max(fGlobPmax, fBrFluxIntProb*safety_factor);
// Update the sum of fBrFluxIntProb*fBrFluxWeight for different species
if(fSumFluxIntProbs.find(fBrFluxPDG) == fSumFluxIntProbs.end()){
fSumFluxIntProbs[fBrFluxPDG] = 0.0;
}
fSumFluxIntProbs[fBrFluxPDG] += fBrFluxIntProb * fBrFluxWeight;
}
LOG("GNUISANCEMCJDriver", pNOTICE) <<
"Updated global probability scale to fGlobPmax = "<< fGlobPmax;
if(save_to_file){
LOG("GNUISANCEMCJDriver", pNOTICE) <<
"Saving pre-generated interaction probabilities to file: "<<
fFluxIntProbFile->GetName();
fFluxIntProbFile->cd();
fFluxIntTree->Write();
}
// Also build index for use later
if(fFluxIntTree->BuildIndex("FluxIndex") != fFluxIntTree->GetEntries()){
LOG("GNUISANCEMCJDriver", pFATAL) <<
"Cannot build index using branch \"FluxIndex\" for flux prob tree!";
exit(1);
}
// Now that have pre-generated flux probabilities need to trun off event
// preselection as this is only advantages when using max path lengths
this->PreSelectEvents(false);
LOG("GNUISANCEMCJDriver", pNOTICE) << "Successfully generated/loaded pre-calculate flux interaction probabilities";
}
// Otherwise clean up
else if(fFluxIntTree){
delete fFluxIntTree;
fFluxIntTree = 0;
}
// Return whether have successfully pre-calculated flux interaction probabilities
return success;
}
//___________________________________________________________________________
bool GNUISANCEMCJDriver::LoadFluxProbabilities(string filename)
{
// Load a pre-generated set of flux interaction probabilities from an external
// file. This is recommended when using large flux files (>100k entries) as
// for these the time to calculate the interaction probabilities can exceed
// ~20 minutes. After loading the input tree we call PreCalcFluxProbabilities
// to check that has successfully loaded
//
if(fFluxIntProbFile){
LOG("GNUISANCEMCJDriver", pWARN)
<< "Can't load flux interaction prob file as one is already loaded";
return false;
}
fFluxIntProbFile = new TFile(filename.c_str(), "OPEN");
if(fFluxIntProbFile){
fFluxIntTree = dynamic_cast<TTree*>(fFluxIntProbFile->Get(fFluxIntTreeName.c_str()));
if(fFluxIntTree){
bool set_addresses =
fFluxIntTree->SetBranchAddress("FluxIntProb", &fBrFluxIntProb) >= 0 &&
fFluxIntTree->SetBranchAddress("FluxIndex", &fBrFluxIndex) >= 0 &&
fFluxIntTree->SetBranchAddress("FluxPDG", &fBrFluxPDG) >= 0 &&
fFluxIntTree->SetBranchAddress("FluxWeight", &fBrFluxWeight) >= 0 &&
fFluxIntTree->SetBranchAddress("FluxEnu", &fBrFluxEnu) >= 0;
if(set_addresses){
// Finally check that can use them
if(this->PreCalcFluxProbabilities()) {
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Successfully loaded pre-generated flux interaction probabilities";
return true;
}
}
// If cannot load then delete tree
LOG("GNUISANCEMCJDriver", pERROR) <<
"Cannot find expected branches in input flux probability tree!";
delete fFluxIntTree; fFluxIntTree = 0;
}
else LOG("GNUISANCEMCJDriver", pERROR)
<< "Cannot find tree: "<< fFluxIntTreeName.c_str();
}
LOG("GNUISANCEMCJDriver", pWARN)
<< "Unable to load flux interaction probabilities file";
return false;
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::SaveFluxProbabilities(string outfilename)
{
// Configue the flux driver to save the calculated flux interaction
// probabilities to the specified output file name for use in later jobs. See
// the LoadFluxProbTree method for how they are fed into a later job.
//
fFluxIntFileName = outfilename;
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::Configure(bool calc_prob_scales)
{
LOG("GNUISANCEMCJDriver", pNOTICE)
<< utils::print::PrintFramedMesg("Configuring GNUISANCEMCJDriver");
// Get the list of neutrino types from the input flux driver and the list
// of target materials from the input geometry driver
this->GetParticleLists();
// Ask the input GFluxI for the max. neutrino energy (to compute Pmax)
this->GetMaxFluxEnergy();
// Create all possible initial states and for each one initialize,
// configure & store an GEVGDriver event generation driver object.
// Once an 'initial state' has been selected from the input flux / geom,
// the responsibility for generating the neutrino interaction will be
// delegated to one of these drivers.
this->PopulateEventGenDriverPool();
// If the user wants to use cross section splines in order to speed things
// up, then coordinate spline creation from all GEVGDriver objects pushed
// into GEVGPool. This will create all xsec splines needed for all (enabled)
// processes that can be simulated involving the particles in the input flux
// and geometry.
// Spline creation will be skipped for every spline that has been pre-loaded
// into the the XSecSplineList.
// Once more it is noted that computing cross section splines is a huge
// overhead. The user is encouraged to generate them in advance and load
// them into the XSecSplineList
this->BootstrapXSecSplines();
// Create cross section splines describing the total interaction xsec
// for a given initial state (Create them by summing all xsec splines
// for each possible initial state)
this->BootstrapXSecSplineSummation();
if(calc_prob_scales){
// Ask the input geometry driver to compute the max. path length for each
// material in the list of target materials (or load a precomputed list)
this->GetMaxPathLengthList();
// Compute the max. interaction probability to scale all interaction
// probabilities to be computed by this driver
this->ComputeProbScales();
}
LOG("GNUISANCEMCJDriver", pNOTICE) << "Finished configuring GNUISANCEMCJDriver\n\n";
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::InitJob(void)
{
fEventGenList = "Default"; // <-- set of event generators to be loaded by this driver
fUnphysEventMask = new TBits(GHepFlags::NFlags()); //<-- unphysical event mask
//fUnphysEventMask->ResetAllBits(true);
for(unsigned int i = 0; i < GHepFlags::NFlags(); i++) {
fUnphysEventMask->SetBitNumber(i, true);
}
fFluxDriver = 0; // <-- flux driver
fGeomAnalyzer = 0; // <-- geometry driver
fGPool = 0; // <-- pool of GEVGDriver event generation drivers
fEmax = 0; // <-- maximum neutrino energy
fMaxPlXmlFilename = ""; // <-- XML file with external path lengths
fUseExtMaxPl = false;
fUseSplines = false;
fNFluxNeutrinos = 0; // <-- number of flux neutrinos thrown so far
fGlobPmax = 0; // <-- maximum interaction probability (global prob scale)
fPmax.clear(); // <-- maximum interaction probability per neutrino & per energy bin
fGenerateUnweighted = false; // <-- default opt to generate weighted events
fPreSelect = true; // <-- default to use pre-selection based on maximum path lengths
fSelTgtPdg = 0;
fCurEvt = 0;
fCurVtx.SetXYZT(0.,0.,0.,0.);
fFluxIntProbFile = 0;
fFluxIntTreeName = "gFlxIntProb";
fFluxIntFileName = "";
fFluxIntTree = 0;
fBrFluxIntProb = -1.;
fBrFluxIndex = -1;
fBrFluxEnu = -1.;
fBrFluxWeight = -1.;
fBrFluxPDG = 0;
fSumFluxIntProbs.clear();
// Throw as many flux neutrinos as necessary till one has interacted
// so that GenerateEvent() never returns NULL (except when in error)
this->KeepOnThrowingFluxNeutrinos(true);
// Allow the selected GEVGDriver to go into recursive mode and regenerate
// an interaction that turns out to be unphysical.
//TBits unphysmask(GHepFlags::NFlags());
//unphysmask.ResetAllBits(false);
//this->FilterUnphysical(unphysmask);
// Force early initialization of singleton objects that are typically
// would be initialized at their first use later on.
// This is purely cosmetic and I do it to send the banner and some prolific
// initialization printout at the top.
assert( Messenger::Instance() );
assert( AlgConfigPool::Instance() );
// Autoload splines (from the XML file pointed at the $GSPLOAD env. var.,
// if the env. var. has been set);
XSecSplineList * xspl = XSecSplineList::Instance();
xspl->AutoLoad();
// Clear the target and neutrino lists
fNuList.clear();
fTgtList.clear();
// Clear the maximum path length list
fMaxPathLengths.clear();
fCurPathLengths.clear();
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::GetParticleLists(void)
{
// Get the list of flux neutrinos from the flux driver
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Asking the flux driver for its list of neutrinos";
fNuList = fFluxDriver->FluxParticles();
LOG("GNUISANCEMCJDriver", pNOTICE) << "Flux particles: " << fNuList;
// Get the list of target materials from the geometry driver
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Asking the geometry driver for its list of targets";
fTgtList = fGeomAnalyzer->ListOfTargetNuclei();
LOG("GNUISANCEMCJDriver", pNOTICE) << "Target materials: " << fTgtList;
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::GetMaxPathLengthList(void)
{
if(fUseExtMaxPl) {
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Loading external max path-length list for input geometry from "
<< fMaxPlXmlFilename;
fMaxPathLengths.LoadFromXml(fMaxPlXmlFilename);
} else {
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Querying the geometry driver to compute the max path-length list";
fMaxPathLengths = fGeomAnalyzer->ComputeMaxPathLengths();
}
// Print maximum path lengths & neutrino energy
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Maximum path length list: " << fMaxPathLengths;
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::GetMaxFluxEnergy(void)
{
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Querying the flux driver for the maximum energy of flux neutrinos";
fEmax = fFluxDriver->MaxEnergy();
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Maximum flux neutrino energy = " << fEmax << " GeV";
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::PopulateEventGenDriverPool(void)
{
LOG("GNUISANCEMCJDriver", pDEBUG)
<< "Creating GEVGPool & adding a GEVGDriver object per init-state";
if (fGPool) delete fGPool;
fGPool = new GEVGPool;
PDGCodeList::const_iterator nuiter;
PDGCodeList::const_iterator tgtiter;
for(nuiter = fNuList.begin(); nuiter != fNuList.end(); ++nuiter) {
for(tgtiter = fTgtList.begin(); tgtiter != fTgtList.end(); ++tgtiter) {
int target_pdgc = *tgtiter;
int neutrino_pdgc = *nuiter;
InitialState init_state(target_pdgc, neutrino_pdgc);
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "\n\n ---- Creating a GEVGDriver object configured for init-state: "
<< init_state.AsString() << " ----\n\n";
GEVGDriver * evgdriver = new GEVGDriver;
evgdriver->SetEventGeneratorList(fEventGenList); // specify list of generators
evgdriver->Configure(init_state);
evgdriver->UseSplines(); // check if all splines needed are loaded
LOG("GNUISANCEMCJDriver", pDEBUG) << "Adding new GEVGDriver object to GEVGPool";
fGPool->insert( GEVGPool::value_type(init_state.AsString(), evgdriver) );
} // targets
} // neutrinos
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "All necessary GEVGDriver object were pushed into GEVGPool\n";
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::BootstrapXSecSplines(void)
{
// Bootstrap cross section spline generation by the event generation drivers
// that handle each initial state.
if(!fUseSplines) return;
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Asking event generation drivers to compute all needed xsec splines";
PDGCodeList::const_iterator nuiter;
PDGCodeList::const_iterator tgtiter;
for(nuiter = fNuList.begin(); nuiter != fNuList.end(); ++nuiter){
for(tgtiter = fTgtList.begin(); tgtiter != fTgtList.end(); ++tgtiter) {
int target_pdgc = *tgtiter;
int neutrino_pdgc = *nuiter;
InitialState init_state(target_pdgc, neutrino_pdgc);
LOG("GNUISANCEMCJDriver", pINFO)
<< "Computing all splines needed for init-state: "
<< init_state.AsString();
GEVGDriver * evgdriver = fGPool->FindDriver(init_state);
evgdriver->CreateSplines(-1,-1,fUseLogE);
} // targets
} // neutrinos
LOG("GNUISANCEMCJDriver", pINFO) << "Finished creating cross section splines\n";
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::BootstrapXSecSplineSummation(void)
{
// Sum-up the cross section splines for all the interaction that can be
// simulated for each initial state
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Summing-up splines to get total cross section for each init state";
GEVGPool::iterator diter;
for(diter = fGPool->begin(); diter != fGPool->end(); ++diter) {
string init_state = diter->first;
GEVGDriver * evgdriver = diter->second;
assert(evgdriver);
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "**** Summing xsec splines for init-state = " << init_state;
Range1D_t rE = evgdriver->ValidEnergyRange();
if (fEmax>rE.max || fEmax<rE.min)
LOG("GNUISANCEMCJDriver",pFATAL)
<< " rE (validEnergyRange) [" << rE.min << "," << rE.max << "] "
<< " fEmax " << fEmax;
assert(fEmax<rE.max && fEmax>rE.min);
// decide the energy range for the sum spline - extend the spline a little
// bit above the maximum beam energy (but below the maximum valid energy)
// to avoid the evaluation of the cubic spline around the viscinity of
// knots with zero y values (although the GENIE Spline object handles it)
double dE = fEmax/10.;
double min = rE.min;
double max = (fEmax+dE < rE.max) ? fEmax+dE : rE.max;
evgdriver->CreateXSecSumSpline(100,min,max,true);
}
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Finished summing all interaction xsec splines per initial state";
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::ComputeProbScales(void)
{
// Computing interaction probability scales.
// To minimize the numbers or trials before choosing a neutrino+target init
// state for generating an event (note: each trial means selecting a flux
// neutrino, navigating it through the detector to compute path lengths,
// computing interaction probabilities for each material and so on...)
// a set of probability scales can be used for different neutrino species
// and at different energy bins.
// A global probability scale is also being constructed for keeping the correct
// proportions between differect flux neutrino species or flux neutrinos of
// different energies.
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Computing the max. interaction probability (probability scale)";
// clean up global probability scale and maximum probabilties per neutrino
// type & energy bin
{
fGlobPmax = 0;
map<int,TH1D*>::iterator pmax_iter = fPmax.begin();
for( ; pmax_iter != fPmax.end(); ++pmax_iter) {
TH1D * pmax = pmax_iter->second;
if(pmax) {
delete pmax; pmax = 0;
}
}
fPmax.clear();
}
// for maximum interaction probability vs E /for given geometry/ I will
// be using 300 bins up to the maximum energy for the input flux
// double de = fEmax/300.;//djk june 5, 2013
double de = fEmax/300.;//djk june 5, 2013
double emin = 0.0;
double emax = fEmax + de;
int n = 1 + (int) ((emax-emin)/de);
PDGCodeList::const_iterator nuiter;
PDGCodeList::const_iterator tgtiter;
// loop over all neutrino types generated by the flux driver
for(nuiter = fNuList.begin(); nuiter != fNuList.end(); ++nuiter) {
int neutrino_pdgc = *nuiter;
TH1D * pmax_hst = new TH1D("pmax_hst",
"max interaction probability vs E | geom",n,emin,emax);
pmax_hst->SetDirectory(0);
// loop over energy bins
for(int ie = 1; ie <= pmax_hst->GetNbinsX(); ie++) {
double EvLow = pmax_hst->GetBinCenter(ie) - 0.5*pmax_hst->GetBinWidth(ie);
double EvHigh = pmax_hst->GetBinCenter(ie) + 0.5*pmax_hst->GetBinWidth(ie);
//double Ev = pmax_hst->GetBinCenter(ie);
// loop over targets in input geometry, form initial state and compute
// the sum of maximum interaction probabilities at the current energy bin
//
for(tgtiter = fTgtList.begin(); tgtiter != fTgtList.end(); ++tgtiter) {
int target_pdgc = *tgtiter;
InitialState init_state(target_pdgc, neutrino_pdgc);
LOG("GNUISANCEMCJDriver", pDEBUG)
<< "Computing Pmax for init-state: " << init_state.AsString() << " E from " << EvLow << "-" << EvHigh;
// get the appropriate driver
GEVGDriver * evgdriver = fGPool->FindDriver(init_state);
// get xsec sum over all modelled processes for given neutrino+target)
double sxsecLow = evgdriver->XSecSumSpline()->Evaluate(EvLow);
double sxsecHigh = evgdriver->XSecSumSpline()->Evaluate(EvHigh);
// get max{path-length x density}
double plmax = fMaxPathLengths.PathLength(target_pdgc);
// compute/store the max interaction probabiity (for given energy)
int A = pdg::IonPdgCodeToA(target_pdgc);
double pmaxLow = this->InteractionProbability(sxsecLow, plmax, A);
double pmaxHigh = this->InteractionProbability(sxsecHigh, plmax, A);
double pmax = pmaxHigh;
if ( pmaxLow > pmaxHigh){
pmax = pmaxLow;
LOG("GNUISANCEMCJDriver", pWARN)
<< "Lower energy neutrinos have a higher probability of interacting than those at higher energy."
<< " pmaxLow(E=" << EvLow << ")=" << pmaxLow << " and " << " pmaxHigh(E=" << EvHigh << ")=" << pmaxHigh;
}
pmax_hst->SetBinContent(ie, pmax_hst->GetBinContent(ie) + pmax);
LOG("GNUISANCEMCJDriver", pDEBUG)
<< "Pmax[" << init_state.AsString() << ", Ev from " << EvLow << "-" << EvHigh << "] = " << pmax;
} // targets
pmax_hst->SetBinContent(ie, 1.2 * pmax_hst->GetBinContent(ie));
LOG("GNUISANCEMCJDriver", pINFO)
<< "Pmax[nu=" << neutrino_pdgc << ", Ev from " << EvLow << "-" << EvHigh << "] = "
<< pmax_hst->GetBinContent(ie);
} // E
fPmax.insert(map<int,TH1D*>::value_type(neutrino_pdgc,pmax_hst));
} // nu
// Compute global probability scale
// Sum Probabilities {
// all neutrinos, all targets, @ max path length, @ max energy}
//
{
for(nuiter = fNuList.begin(); nuiter != fNuList.end(); ++nuiter) {
int neutrino_pdgc = *nuiter;
map<int,TH1D*>::const_iterator pmax_iter = fPmax.find(neutrino_pdgc);
assert(pmax_iter != fPmax.end());
TH1D * pmax_hst = pmax_iter->second;
assert(pmax_hst);
// double pmax = pmax_hst->GetBinContent(pmax_hst->FindBin(fEmax));
double pmax = pmax_hst->GetMaximum();
assert(pmax>0);
// fGlobPmax += pmax;
fGlobPmax = TMath::Max(pmax, fGlobPmax); // ?;
}
LOG("GNUISANCEMCJDriver", pNOTICE) << "*** Probability scale = " << fGlobPmax;
}
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::InitEventGeneration(void)
{
fCurPathLengths.clear();
fCurEvt = 0;
fSelTgtPdg = 0;
fCurVtx.SetXYZT(0.,0.,0.,0.);
}
//___________________________________________________________________________
EventRecord * GNUISANCEMCJDriver::GenerateEvent(void)
{
LOG("GNUISANCEMCJDriver", pNOTICE) << "Generating next event...";
this->InitEventGeneration();
while(1) {
bool flux_end = fFluxDriver->End();
if(flux_end) {
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "No more neutrinos can be thrown by the flux driver";
return 0;
}
EventRecord * event = this->GenerateEvent1Try();
if(event) return event;
if(fKeepThrowingFluxNu) {
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Flux neutrino didn't interact - Trying the next one...";
continue;
}
break;
} // (w(1)
LOG("GNUISANCEMCJDriver", pINFO) << "Returning NULL event!";
return 0;
}
//___________________________________________________________________________
EventRecord * GNUISANCEMCJDriver::GenerateEvent1Try(void)
{
// attempt generating a neutrino interaction by firing a single flux neutrino
//
RandomGen * rnd = RandomGen::Instance();
double Pno=0, Psum=0;
double R = rnd->RndEvg().Rndm();
LOG("GNUISANCEMCJDriver", pDEBUG) << "Rndm [0,1] = " << R;
// Generate a neutrino using the input GFluxI & get current pdgc/p4/x4
bool flux_ok = this->GenerateFluxNeutrino();
if(!flux_ok) {
LOG("GNUISANCEMCJDriver", pERROR)
<< "** Rejecting current flux neutrino (flux driver err)";
return 0;
}
// Compute the interaction probabilities assuming max. path lengths
// and decide whether the neutrino would interact --
// Many flux neutrinos should be rejected here, drastically reducing
// the number of neutrinos that I need to propagate through the
// actual detector geometry (this is skipped when using
// pre-calculated flux interaction probabilities)
if(fPreSelect) {
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Computing interaction probabilities for max. path lengths";
Psum = this->ComputeInteractionProbabilities(true /* <- max PL*/);
Pno = 1-Psum;
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "The no-interaction probability (max. path lengths) is: "
<< 100*Pno << " %";
if(Pno<0.) {
LOG("GNUISANCEMCJDriver", pFATAL)
<< "Negative no-interaction probability! (P = " << 100*Pno << " %)"
<< " Particle E=" << fFluxDriver->Momentum().E() << " type=" << fFluxDriver->PdgCode() << "Psum=" << Psum;
gAbortingInErr=true;
exit(1);
}
if(R>=1-Pno) {
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "** Rejecting current flux neutrino";
return 0;
}
} // preselect
bool pl_ok = false;
// If possible use pre-generated flux neutrino interaction probabilities
if(fFluxIntTree){
Psum = this->PreGenFluxInteractionProbability();
}
// Else compute them in the usual manner
else {
// Compute (pathLength x density x weight fraction) for all materials
// in the input geometry, for the neutrino generated by the flux driver
pl_ok = this->ComputePathLengths();
if(!pl_ok) {
LOG("GNUISANCEMCJDriver", pERROR)
<< "** Rejecting current flux neutrino (err computing path-lengths)";
return 0;
}
if(fCurPathLengths.AreAllZero()) {
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "** Rejecting current flux neutrino (misses generation volume)";
return 0;
}
Psum = this->ComputeInteractionProbabilities(false /* <- actual PL */);
}
if(TMath::Abs(Psum) < controls::kASmallNum){
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "** Rejecting current flux neutrino (has null interaction probability)";
return 0;
}
// Now decide whether the current neutrino interacts
Pno = 1-Psum;
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "The actual 'no interaction' probability is: " << 100*Pno << " %";
if(Pno<0.) {
LOG("GNUISANCEMCJDriver", pFATAL)
<< "Negative no interactin probability! (P = " << 100*Pno << " %)";
// print info about what caused the problem
int nupdg = fFluxDriver -> PdgCode ();
const TLorentzVector & nup4 = fFluxDriver -> Momentum ();
const TLorentzVector & nux4 = fFluxDriver -> Position ();
LOG("GNUISANCEMCJDriver", pWARN)
<< "\n [-] Problematic neutrino: "
<< "\n |----o PDG-code : " << nupdg
<< "\n |----o 4-momentum : " << utils::print::P4AsString(&nup4)
<< "\n |----o 4-position : " << utils::print::X4AsString(&nux4)
<< "\n Emax : " << fEmax;
LOG("GNUISANCEMCJDriver", pWARN)
<< "\n Problematic path lengths:" << fCurPathLengths;
LOG("GNUISANCEMCJDriver", pWARN)
<< "\n Maximum path lengths:" << fMaxPathLengths;
exit(1);
}
if(R>=1-Pno) {
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "** Rejecting current flux neutrino";
return 0;
}
//
// The flux neutrino interacts!
//
// Calculate path lengths for first time and check potential mismatch if
// used pre-generated flux interaction probabilities
if(fFluxIntTree){
pl_ok = this->ComputePathLengths();
if(!pl_ok) {
LOG("GNUISANCEMCJDriver", pFATAL) << "** Cannot calculate path lenths!";
exit(1);
}
double Psum_curr = this->ComputeInteractionProbabilities(false /* <- actual PL */);
bool mismatch = TMath::Abs(Psum-Psum_curr) > controls::kASmallNum;
if(mismatch){
LOG("GNUISANCEMCJDriver", pFATAL) <<
"** Mismatch between pre-calculated and current interaction "<<
"probabilities!";
exit(1);
}
}
// Select a target material
fSelTgtPdg = this->SelectTargetMaterial(R);
if(fSelTgtPdg==0) {
LOG("GNUISANCEMCJDriver", pERROR)
<< "** Rejecting current flux neutrino (failed to select tgt!)";
return 0;
}
// Ask the GEVGDriver object to select and generate an interaction and
// its kinematics for the selected initial state & neutrino 4-momentum
this->GenerateEventKinematics();
if(!fCurEvt) {
LOG("GNUISANCEMCJDriver", pWARN)
<< "** Couldn't generate kinematics for selected interaction";
return 0;
}
// Generate an 'interaction position' in the selected material (in the
// detector coord system), along the direction of nup4 & set it
this->GenerateVertexPosition();
// Set the event probability (probability for this event to happen given
// the detector setup & the selected flux neutrino)
// Note for users:
// The above probability is stored at GHepRecord::Probability()
// For normalization purposes make sure that you take into account the
// GHepRecord::Weight() -if event generation is weighted-, and
// GFluxI::Weight() -if beam simulation is weighted-.
this->ComputeEventProbability();
return fCurEvt;
}
//___________________________________________________________________________
bool GNUISANCEMCJDriver::GenerateFluxNeutrino(void)
{
// Ask the neutrino flux driver to generate a flux neutrino and make sure
// that things look ok...
//
LOG("GNUISANCEMCJDriver", pNOTICE) << "Generating a flux neutrino";
bool ok = fFluxDriver->GenerateNext();
if(!ok) {
LOG("GNUISANCEMCJDriver", pERROR)
<< "*** The flux driver couldn't generate a flux neutrino!!";
return false;
}
fNFluxNeutrinos++;
int nupdg = fFluxDriver -> PdgCode ();
const TLorentzVector & nup4 = fFluxDriver -> Momentum ();
const TLorentzVector & nux4 = fFluxDriver -> Position ();
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "\n [-] Generated flux neutrino: "
<< "\n |----o PDG-code : " << nupdg
<< "\n |----o 4-momentum : " << utils::print::P4AsString(&nup4)
<< "\n |----o 4-position : " << utils::print::X4AsString(&nux4);
if(nup4.Energy() > fEmax) {
LOG("GNUISANCEMCJDriver", pFATAL)
<< "\n *** Flux driver error ***"
<< "\n Generated flux v with E = " << nup4.Energy() << " GeV"
<< "\n Max v energy (declared by flux driver) = " << fEmax << " GeV"
<< "\n My interaction probability scaling is invalidated!!";
return false;
}
if(!fNuList.ExistsInPDGCodeList(nupdg)) {
LOG("GNUISANCEMCJDriver", pFATAL)
<< "\n *** Flux driver error ***"
<< "\n Generated flux v with pdg = " << nupdg
<< "\n It does not belong to the declared list of flux neutrinos"
<< "\n I was not configured to handle this!!";
return false;
}
return true;
}
//___________________________________________________________________________
bool GNUISANCEMCJDriver::ComputePathLengths(void)
{
// Ask the geometry driver to compute (pathLength x density x weight frac.)
// for all detector materials for the neutrino generated by the flux driver
// and make sure that things look ok...
fCurPathLengths.clear();
const TLorentzVector & nup4 = fFluxDriver -> Momentum ();
const TLorentzVector & nux4 = fFluxDriver -> Position ();
fCurPathLengths = fGeomAnalyzer->ComputePathLengths(nux4, nup4);
LOG("GNUISANCEMCJDriver", pNOTICE) << fCurPathLengths;
if(fCurPathLengths.size() == 0) {
LOG("GNUISANCEMCJDriver", pFATAL)
<< "\n *** Geometry driver error ***"
<< "\n Got an empty PathLengthList - No material found in geometry?";
return false;
}
if(fCurPathLengths.AreAllZero()) {
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "current flux v doesn't cross any geometry material...";
}
return true;
}
//___________________________________________________________________________
double GNUISANCEMCJDriver::ComputeInteractionProbabilities(bool use_max_path_length)
{
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Computing relative interaction probabilities for each material";
// current flux neutrino code & 4-p
int nupdg = fFluxDriver->PdgCode();
const TLorentzVector & nup4 = fFluxDriver->Momentum();
fCurCumulProbMap.clear();
const PathLengthList & path_length_list =
(use_max_path_length) ? fMaxPathLengths : fCurPathLengths;
double probsum=0;
PathLengthList::const_iterator pliter;
for(pliter = path_length_list.begin();
pliter != path_length_list.end(); ++pliter) {
int mpdg = pliter->first; // material PDG code
double pl = pliter->second; // density x path-length
int A = pdg::IonPdgCodeToA(mpdg);
double xsec = 0.; // sum of xsecs for all modelled processes for given init state
double prob = 0.; // interaction probability
double probn = 0.; // normalized interaction probability
// find the GEVGDriver object that is handling the current init state
InitialState init_state(mpdg, nupdg);
GEVGDriver * evgdriver = fGPool->FindDriver(init_state);
if(!evgdriver) {
LOG("GNUISANCEMCJDriver", pFATAL)
<< "\n * The MC Job driver isn't properly configured!"
<< "\n * No event generation driver could be found for init state: "
<< init_state.AsString();
exit(1);
}
// compute the interaction xsec and probability (if path-length>0)
if(pl>0.) {
const Spline * totxsecspl = evgdriver->XSecSumSpline();
if(!totxsecspl) {
LOG("GNUISANCEMCJDriver", pFATAL)
<< "\n * The MC Job driver isn't properly configured!"
<< "\n * Couldn't retrieve total cross section spline for init state: "
<< init_state.AsString();
exit(1);
} else {
xsec = totxsecspl->Evaluate( nup4.Energy() );
}
prob = this->InteractionProbability(xsec,pl,A);
LOG("GNUISANCEMCJDriver", pDEBUG)
<< " (xsec, pl, A)=(" << xsec << "," << pl << "," << A << ")";
// scale the interaction probability to the maximum one so as not
// to have to throw few billions of flux neutrinos before getting
// an interaction...
double pmax = 0;
if(fGenerateUnweighted) pmax = fGlobPmax;
else {
map<int,TH1D*>::const_iterator pmax_iter = fPmax.find(nupdg);
assert(pmax_iter != fPmax.end());
TH1D * pmax_hst = pmax_iter->second;
assert(pmax_hst);
int ie = pmax_hst->FindBin(nup4.Energy());
pmax = pmax_hst->GetBinContent(ie);
}
assert(pmax>0);
LOG("GNUISANCEMCJDriver", pDEBUG)
<< "Pmax=" << pmax;
probn = prob/pmax;
}
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "tgt: " << mpdg << " -> TotXSec = "
<< xsec/units::cm2 << " cm^2, Norm.Prob = " << 100*probn << "%";
#endif
probsum += probn;
fCurCumulProbMap.insert(map<int,double>::value_type(mpdg,probsum));
}
return probsum;
}
//___________________________________________________________________________
int GNUISANCEMCJDriver::SelectTargetMaterial(double R)
{
// Pick a target material using the pre-computed interaction probabilities
// for a flux neutrino that has already been determined that interacts
LOG("GNUISANCEMCJDriver", pNOTICE) << "Selecting target material";
int tgtpdg = 0;
map<int,double>::const_iterator probiter = fCurCumulProbMap.begin();
for( ; probiter != fCurCumulProbMap.end(); ++probiter) {
double prob = probiter->second;
if(R<prob) {
tgtpdg = probiter->first;
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Selected target material = " << tgtpdg;
return tgtpdg;
}
}
LOG("GNUISANCEMCJDriver", pERROR)
<< "Could not select target material for an interacting neutrino";
return 0;
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::GenerateEventKinematics(void)
{
int nupdg = fFluxDriver->PdgCode();
const TLorentzVector & nup4 = fFluxDriver->Momentum();
// Find the GEVGDriver object that generates interactions for the
// given initial state (neutrino + target)
InitialState init_state(fSelTgtPdg, nupdg);
GEVGDriver * evgdriver = fGPool->FindDriver(init_state);
if(!evgdriver) {
LOG("GNUISANCEMCJDriver", pFATAL)
<< "No GEVGDriver object for init state: " << init_state.AsString();
exit(1);
}
// propagate current unphysical event mask
evgdriver->SetUnphysEventMask(*fUnphysEventMask);
// Ask the GEVGDriver object to select and generate an interaction for
// the selected initial state & neutrino 4-momentum
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Asking the selected GEVGDriver object to generate an event";
fCurEvt = evgdriver->GenerateEvent(nup4);
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::GenerateVertexPosition(void)
{
// Generate an 'interaction position' in the selected material, along
// the direction of nup4
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "Asking the geometry analyzer to generate a vertex";
const TLorentzVector & p4 = fFluxDriver->Momentum ();
const TLorentzVector & x4 = fFluxDriver->Position ();
const TVector3 & vtx = fGeomAnalyzer->GenerateVertex(x4, p4, fSelTgtPdg);
TVector3 origin(x4.X(), x4.Y(), x4.Z());
origin-=vtx; // computes vector dr = origin - vtx
double dL = origin.Mag();
double c = kLightSpeed /(units::meter/units::second);
double dt = dL/c;
LOG("GNUISANCEMCJDriver", pNOTICE)
<< "|vtx - origin|: dL = " << dL << " m, dt = " << dt << " sec";
fCurVtx.SetXYZT(vtx.x(), vtx.y(), vtx.z(), x4.T() + dt);
fCurEvt->SetVertex(fCurVtx);
}
//___________________________________________________________________________
void GNUISANCEMCJDriver::ComputeEventProbability(void)
{
// Compute event probability for the given flux neutrino & detector geometry
// get interaction cross section
double xsec = fCurEvt->XSec();
// get path length in detector along v direction for specified target material
PathLengthList::const_iterator pliter = fCurPathLengths.find(fSelTgtPdg);
double path_length = pliter->second;
// get target material mass number
int A = pdg::IonPdgCodeToA(fSelTgtPdg);
// calculate interaction probability
double P = this->InteractionProbability(xsec, path_length, A);
//
// get weight for selected event
//
GHepParticle * nu = fCurEvt->Probe();
int nu_pdg = nu->Pdg();
double Ev = nu->P4()->Energy();
double weight = 1.0;
if(!fGenerateUnweighted) {
map<int,TH1D*>::const_iterator pmax_iter = fPmax.find(nu_pdg);
assert(pmax_iter != fPmax.end());
TH1D * pmax_hst = pmax_iter->second;
assert(pmax_hst);
double pmax = pmax_hst->GetBinContent(pmax_hst->FindBin(Ev));
assert(pmax>0);
weight = pmax/fGlobPmax;
}
// set probability & update weight
fCurEvt->SetProbability(P);
fCurEvt->SetWeight(weight * fCurEvt->Weight());
}
//___________________________________________________________________________
double GNUISANCEMCJDriver::InteractionProbability(double xsec, double pL, int A)
{
// P = Na (Avogadro number, atoms/mole) *
// 1/A (1/mass number, mole/gr) *
// xsec (total interaction cross section, cm^2) *
// pL (density-weighted path-length, gr/cm^2)
//
xsec = xsec / units::cm2;
pL = pL * ((units::kilogram/units::m2)/(units::gram/units::cm2));
return kNA*(xsec*pL)/A;
}
//___________________________________________________________________________
double GNUISANCEMCJDriver::PreGenFluxInteractionProbability()
{
// Return the pre-computed interaction probability for the current flux
// neutrino index (entry number in flux file). Exit if not possible as
// using meaningless interaction probability leads to incorrect physics
//
if(!fFluxIntTree){
LOG("GNUISANCEMCJDriver", pERROR) <<
"Cannot get pre-computed flux interaction probability as no tree!";
exit(1);
}
assert(fFluxDriver->Index() >= 0); // Check trying to find meaningfull index
// Check if can find relevant entry and no mismatch in energies -->
// using correct pre-gen interaction prob file
bool found_entry = fFluxIntTree->GetEntryWithIndex(fFluxDriver->Index()) > 0;
bool enu_match = false;
if(found_entry){
double rel_err = fBrFluxEnu-fFluxDriver->Momentum().E();
if(fBrFluxEnu > controls::kASmallNum) rel_err /= fBrFluxEnu;
enu_match = TMath::Abs(rel_err)<controls::kASmallNum;
if(enu_match == false){
LOG("GNUISANCEMCJDriver", pERROR) <<
"Mismatch between: Enu_curr = "<< fFluxDriver->Momentum().E() <<
", Enu_pre_gen = "<< fBrFluxEnu;
}
}
else {
LOG("GNUISANCEMCJDriver", pERROR) << "Cannot find flux entry in interaction prob tree!";
}
// Exit if not successful
bool success = found_entry && enu_match;
if(!success){
LOG("GNUISANCEMCJDriver", pFATAL) <<
"Cannot find pre-generated interaction probability! Check you "<<
"are using the correct pre-generated interaction prob file " <<
"generated using current flux input file with same input " <<
"config (same geom TopVol, neutrino species list)";
exit(1);
}
assert(fGlobPmax+controls::kASmallNum>0.0);
return fBrFluxIntProb/fGlobPmax;
}
//___________________________________________________________________________
+#endif
diff --git a/src/FitBase/GNUISANCEMCJDriver.h b/src/FitBase/GNUISANCEMCJDriver.h
index 71ee4d0..600bbc5 100644
--- a/src/FitBase/GNUISANCEMCJDriver.h
+++ b/src/FitBase/GNUISANCEMCJDriver.h
@@ -1,143 +1,145 @@
//____________________________________________________________________________
/*!
\class genie::GMCJDriver
\brief A GENIE `MC Job Driver'. Can be used for setting up complicated event
generation cases involving detailed flux descriptions and detector
geometry descriptions.
\author Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
University of Liverpool & STFC Rutherford Appleton Lab
\created May 25, 2005
\cpright Copyright (c) 2003-2016, GENIE Neutrino MC Generator Collaboration
For the full text of the license visit http://copyright.genie-mc.org
or see $GENIE/LICENSE
*/
//____________________________________________________________________________
#ifndef _GENIE_NUISANCE_MC_JOB_DRIVER_H_
#define _GENIE_NUISANCE_MC_JOB_DRIVER_H_
+#ifdef __GENIE_ENABLED__
#include <string>
#include <map>
#include <TH1D.h>
#include <TLorentzVector.h>
#include <TFile.h>
#include <TTree.h>
#include <TBits.h>
#include "EVGDrivers/PathLengthList.h"
#include "PDG/PDGCodeList.h"
#include "EVGDrivers/GEVGPool.h"
using std::string;
using std::map;
namespace genie {
class EventRecord;
class GFluxI;
class GeomAnalyzerI;
class GENIE;
//class GEVGPool;
class GNUISANCEMCJDriver {
public :
GNUISANCEMCJDriver();
~GNUISANCEMCJDriver();
// configure MC job
void SetEventGeneratorList (string listname);
void SetUnphysEventMask (const TBits & mask);
void UseFluxDriver (GFluxI * flux);
void UseGeomAnalyzer (GeomAnalyzerI * geom);
void UseSplines (bool useLogE = true);
bool UseMaxPathLengths (string xml_filename);
void KeepOnThrowingFluxNeutrinos (bool keep_on);
void ForceSingleProbScale (void);
void PreSelectEvents (bool preselect = true);
bool PreCalcFluxProbabilities (void);
bool LoadFluxProbabilities (string filename);
void SaveFluxProbabilities (string outfilename);
void Configure (bool calc_prob_scales = true);
GEVGPool* GetConfigPool(void){return this->fGPool;};
// generate single neutrino event for input flux & geometry
EventRecord * GenerateEvent (void);
// info needed for computing the generated sample normalization
double GlobProbScale (void) const { return fGlobPmax; }
long int NFluxNeutrinos (void) const { return (long int) fNFluxNeutrinos; }
map<int, double> SumFluxIntProbs(void) const { return fSumFluxIntProbs; }
// input flux and geometry drivers
const GFluxI & FluxDriver (void) const { return *fFluxDriver; }
const GeomAnalyzerI & GeomAnalyzer (void) const { return *fGeomAnalyzer; }
GFluxI * FluxDriverPtr (void) const { return fFluxDriver; }
GeomAnalyzerI * GeomAnalyzerPtr (void) const { return fGeomAnalyzer; }
// private methods:
void InitJob (void);
void InitEventGeneration (void);
void GetParticleLists (void);
void GetMaxPathLengthList (void);
void GetMaxFluxEnergy (void);
void PopulateEventGenDriverPool (void);
void BootstrapXSecSplines (void);
void BootstrapXSecSplineSummation (void);
void ComputeProbScales (void);
EventRecord * GenerateEvent1Try (void);
bool GenerateFluxNeutrino (void);
bool ComputePathLengths (void);
double ComputeInteractionProbabilities (bool use_max_path_length);
int SelectTargetMaterial (double R);
void GenerateEventKinematics (void);
void GenerateVertexPosition (void);
void ComputeEventProbability (void);
double InteractionProbability (double xsec, double pl, int A);
double PreGenFluxInteractionProbability(void);
// private data members:
GEVGPool * fGPool; ///< A pool of GEVGDrivers properly configured event generation drivers / one per init state
GFluxI * fFluxDriver; ///< [input] neutrino flux driver
GeomAnalyzerI * fGeomAnalyzer; ///< [input] detector geometry analyzer
double fEmax; ///< [declared by the flux driver] maximum neutrino energy
PDGCodeList fNuList; ///< [declared by the flux driver] list of neutrino codes
PDGCodeList fTgtList; ///< [declared by the geom driver] list of target codes
PathLengthList fMaxPathLengths; ///< [declared by the geom driver] maximum path length list
PathLengthList fCurPathLengths; ///< [current] path length list for current flux neutrino
TLorentzVector fCurVtx; ///< [current] interaction vertex
EventRecord * fCurEvt; ///< [current] generated event
int fSelTgtPdg; ///< [current] selected target material PDG code
map<int,double> fCurCumulProbMap; ///< [current] cummulative interaction probabilities
double fNFluxNeutrinos; ///< [current] number of flux nuetrinos fired by the flux driver so far
map<int,TH1D*> fPmax; ///< [computed at init] interaction probability scale /neutrino /energy for given geometry
double fGlobPmax; ///< [computed at init] global interaction probability scale for given flux & geometry
string fEventGenList; ///< [config] list of event generators loaded by this driver (what used to be the $GEVGL setting)
TBits * fUnphysEventMask; ///< [config] controls whether unphysical events are returned (what used to be the $GUNPHYSMASK setting)
string fMaxPlXmlFilename; ///< [config] input file with max density-weighted path lengths for all materials
bool fUseExtMaxPl; ///< [config] using external max path length estimate?
bool fUseSplines; ///< [config] compute all needed & not-loaded splines at init
bool fUseLogE; ///< [config] build splines = f(logE) (rather than f(E)) ?
bool fKeepThrowingFluxNu; ///< [config] keep firing flux neutrinos till one of them interacts
bool fGenerateUnweighted; ///< [config] force single probability scale?
bool fPreSelect; ///< [config] set whether to pre-select events using max interaction paths
TFile* fFluxIntProbFile; ///< [input] pre-generated flux interaction probability file
TTree* fFluxIntTree; ///< [computed-or-loaded] pre-computed flux interaction probabilities (expected tree name is "gFlxIntProbs")
double fBrFluxIntProb; ///< flux interaction probability (set to branch:"FluxIntProb")
int fBrFluxIndex; ///< corresponding entry in flux input tree (set to address of branch:"FluxEntry")
double fBrFluxEnu; ///< corresponding flux P4 (set to address of branch:"FluxP4")
double fBrFluxWeight; ///< corresponding flux weight (set to address of branch: "FluxWeight")
int fBrFluxPDG; ///< corresponding flux pdg code (set to address of branch: "FluxPDG")
string fFluxIntFileName; ///< whether to save pre-generated flux tree for use in later jobs
string fFluxIntTreeName; ///< name for tree holding flux probabilities
map<int, double> fSumFluxIntProbs; ///< map where the key is flux pdg code and the value is sum of fBrFluxWeight * fBrFluxIntProb for all these flux neutrinos
};
} // genie namespace
#endif // _GENIE_MC_JOB_DRIVER_H_
+#endif
diff --git a/src/MINERvA/MINERvA_CCQE_XSec_1DQ2_antinu.cxx b/src/MINERvA/MINERvA_CCQE_XSec_1DQ2_antinu.cxx
index f5a7f84..7adf843 100644
--- a/src/MINERvA/MINERvA_CCQE_XSec_1DQ2_antinu.cxx
+++ b/src/MINERvA/MINERvA_CCQE_XSec_1DQ2_antinu.cxx
@@ -1,146 +1,146 @@
// 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 "MINERvA_SignalDef.h"
#include "MINERvA_CCQE_XSec_1DQ2_antinu.h"
//********************************************************************
MINERvA_CCQE_XSec_1DQ2_antinu::MINERvA_CCQE_XSec_1DQ2_antinu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MINERvA_CCQE_XSec_1DQ2_antinu sample. \n" \
- "Target: CH \n" \
- "Flux: MINERvA Forward Horn Current nue + nuebar \n" \
- "Signal: Any event with 1 electron, any nucleons, and no other FS particles \n";
-
+ "Target: CH \n" \
+ "Flux: MINERvA Forward Horn Current Numubar \n" \
+ "Signal: True CCQE/2p2h defined at the vertex level \n";
+
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("Q^{2}_{QE} (GeV^{2})");
fSettings.SetYTitle("d#sigma/dQ_{QE}^{2} (cm^{2}/GeV^{2})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(1.5, 10.0);
fSettings.DefineAllowedTargets("C,H");
isFluxFix = !fSettings.Found("name", "_oldflux");
fullphasespace = !fSettings.Found("name", "_20deg");
// CCQELike plot information
fSettings.SetTitle("MINERvA_CCQE_XSec_1DQ2_antinu");
std::string basedir = FitPar::GetDataBase() + "/MINERvA/CCQE/";
std::string datafilename = "";
std::string covarfilename = "";
// Full Phase Space
if (fullphasespace) {
if (isFluxFix) {
if (fIsShape) {
ERR(WRN) << "SHAPE likelihood comparison not available for MINERvA "
<< "datasets with fixed flux information. NUISANCE will scale MC to match "
<< "data normalization but full covariance will be used. " << std::endl;
}
datafilename = "Q2QE_numubar_data_fluxfix.txt";
covarfilename = "Q2QE_numubar_covar_fluxfix.txt";
} else {
if (fIsShape) {
datafilename = "Q2QE_numubar_data_SHAPE-extracted.txt";
covarfilename = "Q2QE_numubar_covar_SHAPE-extracted.txt";
} else {
datafilename = "Q2QE_numubar_data.txt";
covarfilename = "Q2QE_numubar_covar.txt";
}
}
// Restricted Phase Space
} else {
if (isFluxFix) {
if (fIsShape) {
ERR(WRN) << "SHAPE likelihood comparison not available for MINERvA "
<< "datasets with fixed flux information. NUISANCE will scale MC to match "
<< "data normalization but full covariance will be used. " << std::endl;
}
datafilename = "20deg_Q2QE_numubar_data_fluxfix.txt";
covarfilename = "20deg_Q2QE_numubar_covar_fluxfix.txt";
} else {
if (fIsShape) {
datafilename = "20deg_Q2QE_numubar_data_SHAPE-extracted.txt";
covarfilename = "20deg_Q2QE_numubar_covar_SHAPE-extracted.txt";
} else {
datafilename = "20deg_Q2QE_numubar_data.txt";
covarfilename = "20deg_Q2QE_numubar_covar.txt";
}
}
}
fSettings.SetDataInput( basedir + datafilename );
fSettings.SetCovarInput( basedir + covarfilename );
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) * 13. / 7. / TotalIntegratedFlux();
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
if (!isFluxFix or !fullphasespace) {
SetCorrelationFromTextFile( fSettings.GetCovarInput() );
ScaleCovar(1E76);
} else {
SetCovarFromTextFile( fSettings.GetCovarInput() );
}
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
void MINERvA_CCQE_XSec_1DQ2_antinu::FillEventVariables(FitEvent *event) {
//********************************************************************
if (event->NumFSParticle(-13) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pmu = event->GetHMFSParticle(-13)->fP;
double ThetaMu = Pnu.Vect().Angle(Pmu.Vect());
double q2qe = FitUtils::Q2QErec(Pmu, cos(ThetaMu), 30., true);
fXVar = q2qe;
return;
}
//********************************************************************
bool MINERvA_CCQE_XSec_1DQ2_antinu::isSignal(FitEvent *event) {
//*******************************************************************
return SignalDef::isCCQEnumubar_MINERvA(event, EnuMin, EnuMax, fullphasespace);
}
diff --git a/src/MINERvA/MINERvA_CCQE_XSec_1DQ2_nu.cxx b/src/MINERvA/MINERvA_CCQE_XSec_1DQ2_nu.cxx
index e4d3716..1dbf931 100644
--- a/src/MINERvA/MINERvA_CCQE_XSec_1DQ2_nu.cxx
+++ b/src/MINERvA/MINERvA_CCQE_XSec_1DQ2_nu.cxx
@@ -1,148 +1,148 @@
// 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 "MINERvA_SignalDef.h"
#include "MINERvA_CCQE_XSec_1DQ2_nu.h"
//********************************************************************
MINERvA_CCQE_XSec_1DQ2_nu::MINERvA_CCQE_XSec_1DQ2_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MINERvA_CCQE_XSec_1DQ2_nu sample. \n" \
"Target: CH \n" \
- "Flux: MINERvA Forward Horn Current nue + nuebar \n" \
- "Signal: Any event with 1 electron, any nucleons, and no other FS particles \n";
+ "Flux: MINERvA Forward Horn Current Numu \n" \
+ "Signal: True CCQE/2p2h defined at the vertex level \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("Q^{2}_{QE} (GeV^{2})");
fSettings.SetYTitle("d#sigma/dQ_{QE}^{2} (cm^{2}/GeV^{2})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(1.5, 10.0);
fSettings.DefineAllowedTargets("C,H");
isFluxFix = !fSettings.Found("name", "_oldflux");
fullphasespace = !fSettings.Found("name", "_20deg");
// CCQELike plot information
fSettings.SetTitle("MINERvA_CCQE_XSec_1DQ2_nu");
std::string basedir = FitPar::GetDataBase() + "/MINERvA/CCQE/";
std::string datafilename = "";
std::string covarfilename = "";
// Full Phase Space
if (fullphasespace) {
if (isFluxFix) {
if (fIsShape) {
ERR(WRN) << "SHAPE likelihood comparison not available for MINERvA "
<< "datasets with fixed flux information. NUISANCE will scale MC to match "
<< "data normalization but full covariance will be used. " << std::endl;
}
datafilename = "Q2QE_numu_data_fluxfix.txt";
covarfilename = "Q2QE_numu_covar_fluxfix.txt";
} else {
if (fIsShape) {
datafilename = "Q2QE_numu_data_SHAPE-extracted.txt";
covarfilename = "Q2QE_numu_covar_SHAPE-extracted.txt";
} else {
datafilename = "Q2QE_numu_data.txt";
covarfilename = "Q2QE_numu_covar.txt";
}
}
// Restricted Phase Space
} else {
if (isFluxFix) {
if (fIsShape) {
ERR(WRN) << "SHAPE likelihood comparison not available for MINERvA "
<< "datasets with fixed flux information. NUISANCE will scale MC to match "
<< "data normalization but full covariance will be used. " << std::endl;
}
datafilename = "20deg_Q2QE_numu_data_fluxfix.txt";
covarfilename = "20deg_Q2QE_numu_covar_fluxfix.txt";
} else {
if (fIsShape) {
datafilename = "20deg_Q2QE_numu_data_SHAPE-extracted.txt";
covarfilename = "20deg_Q2QE_numu_covar_SHAPE-extracted.txt";
} else {
datafilename = "20deg_Q2QE_numu_data.txt";
covarfilename = "20deg_Q2QE_numu_covar.txt";
}
}
}
fSettings.SetDataInput( basedir + datafilename );
fSettings.SetCovarInput( basedir + covarfilename );
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 * 13.0 / 6.0 / (fNEvents + 0.)) / TotalIntegratedFlux();
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
if (!isFluxFix or !fullphasespace){
SetCorrelationFromTextFile( fSettings.GetCovarInput() );
ScaleCovar(1E76);
} else {
SetCovarFromTextFile( fSettings.GetCovarInput() );
}
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
void MINERvA_CCQE_XSec_1DQ2_nu::FillEventVariables(FitEvent *event) {
//********************************************************************
if (event->NumFSParticle(13) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double ThetaMu = Pnu.Vect().Angle(Pmu.Vect());
double q2qe = FitUtils::Q2QErec(Pmu, cos(ThetaMu), 34., true);
// Set binning variable
fXVar = q2qe;
return;
}
//********************************************************************
bool MINERvA_CCQE_XSec_1DQ2_nu::isSignal(FitEvent *event) {
//*******************************************************************
return SignalDef::isCCQEnumu_MINERvA(event, EnuMin, EnuMax, fullphasespace);
}
diff --git a/src/MINERvA/SAMPLEREADME b/src/MINERvA/SAMPLEREADME
deleted file mode 100644
index fd8c0f0..0000000
--- a/src/MINERvA/SAMPLEREADME
+++ /dev/null
@@ -1,9 +0,0 @@
-./MINERvA_CC0pi_XSec_1DQ2_nu_proton.cxx
-./MINERvA_CC0pi_XSec_1DEe_nue.cxx
-./MINERvA_CC0pi_XSec_1DFill_nue.cxx
-./MINERvA_CC0pi_XSec_1DQ2_nue.cxx
-./MINERvA_CC0pi_XSec_1DThetae_nue.cxx
-./MINERvA_CCQE_XSec_1DQ2_antinu.cxx
-./MINERvA_CCQE_XSec_1DQ2_joint.cxx
-./MINERvA_CCQE_XSec_1DQ2_nu.cxx
-./MINERvA_CCinc_XSec_2DEavq3_nu.cxx

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 5:42 PM (1 d, 17 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805435
Default Alt Text
(145 KB)

Event Timeline