evgen::GENIEHelper Class Reference

#include <GENIEHelper.h>

Inheritance diagram for evgen::GENIEHelper:

List of all members.

Public Member Functions
	GENIEHelper ()
	GENIEHelper (std::string config, std::string version)
	~GENIEHelper ()
void	Update (const cfg::Config &c)
void	Initialize ()
bool	Sample (sim::MCTruth &truth)
Private Member Functions
void	InitializeGeometry ()
void	InitializeFluxDriver ()
void	PackNuMIFlux (sim::MCTruth &truth)
void	PackMCTruth (genie::EventRecord *record, sim::MCTruth &truth)
Private Attributes
genie::GeomAnalyzerI *	fGeomD
genie::GFluxI *	fFluxD
genie::GMCJDriver *	fDriver
rawdata::det_id_	fDetector
std::string	fFluxType
	ipnd, nd, fd
std::string	fFluxFile
	histogram or ntuple
std::string	fBeamName
	name of file containing histograms or ntuples, or txt
std::vector< TH1D * >	fFluxHistograms
	name of the beam we are simulating
int	fEvents
	histograms for each nu species
double	fTotalGen
	number of events to generate
double	fTotalRequested
	total number of events/pot generated to this point
double	fPOT
	total number of events/pot requested to be generated
double	fFluxNormalization
	total POT to generate
double	fMonoEnergy
	POT per file.
std::vector< int >	fGenFlavors
	energy of monoenergetic neutrinos
std::vector< string >	fEnvironment
	pdg codes for flavors to generate
std::map< int, int >	fFlavorMap
	environmental variables and settings used by genie

---

Constructor & Destructor Documentation

GENIEHelper::GENIEHelper (   )

Definition at line 54 of file GENIEHelper.cxx. References fFlavorMap, and cfg::Observer::SetWatch(). : fTotalGen(0.) { fFlavorMap[12] = kNue; fFlavorMap[-12] = kNueBar; fFlavorMap[14] = kNuMu; fFlavorMap[-14] = kNuMuBar; fFlavorMap[16] = kNuTau; fFlavorMap[-16] = kNuTauBar; this->SetWatch("GENIE","default"); }

GENIEHelper::GENIEHelper (  std::string  config, std::string  version )

Definition at line 69 of file GENIEHelper.cxx. References fFlavorMap, and cfg::Observer::SetWatch(). { fFlavorMap[12] = kNue; fFlavorMap[-12] = kNueBar; fFlavorMap[14] = kNuMu; fFlavorMap[-14] = kNuMuBar; fFlavorMap[16] = kNuTau; fFlavorMap[-16] = kNuTauBar; this->SetWatch(config.c_str(), version.c_str()); }

GENIEHelper::~GENIEHelper (   )

Definition at line 82 of file GENIEHelper.cxx. { }

---

Member Function Documentation

void GENIEHelper::Initialize (   )

initialize the Geometry and Flux drivers turn on following line to speed up driver initialization Definition at line 193 of file GENIEHelper.cxx. References fDriver, fFluxD, fGeomD, InitializeFluxDriver(), and InitializeGeometry(). Referenced by main(). { InitializeGeometry(); InitializeFluxDriver(); fDriver = new genie::GMCJDriver(); fDriver->UseFluxDriver(fFluxD); fDriver->UseGeomAnalyzer(fGeomD); //driver->UseMaxPathLengths(***supply some xml file name***); fDriver->Configure(); fDriver->UseSplines(); fDriver->ForceSingleProbScale(); return; }

void GENIEHelper::InitializeFluxDriver (   )  [private]

z=0 is the front face of the tpc set the number of cycles to run +++++++++this is stupid - to really set it i have to get a value from the MCJDriver and i am not even sure what i have below is approximately correct. end if using ntuple flux files get the geometry of the detector set the beam direction and spot based on the detector beamR, bspot need to be fixed up for FD and ND, the 0.05 angle is an estimate multiply by 1e-2 to put the values in meters, the default for genie making sure we are in genie units now add the different fluxes end loop to add flux histograms to driver weight each species equally in the generation Definition at line 231 of file GENIEHelper.cxx. References geo::Geometry::DetHalfHeight(), geo::Geometry::DetLength(), fDetector, fFlavorMap, fFluxD, fFluxFile, fFluxHistograms, fFluxNormalization, fFluxType, fGenFlavors, fMonoEnergy, fPOT, and geo::Geometry::Instance(). Referenced by Initialize(). { if(fFluxType.compare("ntuple") == 0){ genie::flux::GNuMIFlux* numiFlux = new genie::flux::GNuMIFlux(); numiFlux->LoadBeamSimData(fFluxFile); numiFlux->SetFilePOT(fFluxNormalization); numiFlux->SetUpstreamZ(-1000.); numiFlux->SetNumOfCycles(int(fPOT/fFluxNormalization)); fFluxD = dynamic_cast<genie::GFluxI *>(numiFlux); } else if(fFluxType.compare("histogram") == 0){ geo::Geometry& geo = geo::Geometry::Instance(0, fDetector); TVector3 bdir(0.015,0.104,sqrt(1. - 0.104*0.104 - 0.015*0.015)); double beamR = 0.5*(2.*geo.DetHalfHeight()*cos(0.11) + geo.DetLength()*sin(0.11)); double epsilon = beamR - geo.DetLength()*sin(0.11); double cx = 0.; double cy = epsilon*cos(0.11)-geo.DetHalfHeight(); double cz = -epsilon*sin(0.11); if(fDetector == rawdata::kNear){ bdir.SetXYZ(0.014, -2.7*3.1415/180., sqrt(1. - 0.014*0.014 - 2.7*2.7*3.1415*3.1415/(180.*180.))); beamR = 0.5*(2.*geo.DetHalfHeight()*cos(0.05) + geo.DetLength()*sin(0.05)); epsilon = beamR - geo.DetLength()*sin(0.05); cx = 0.; cy = geo.DetHalfHeight()-epsilon*cos(0.05); cz = -epsilon*sin(0.05); } else if(fDetector == rawdata::kFar){ bdir.SetXYZ(0.014, 2.7*3.1415/180., sqrt(1. - 0.014*0.014 - 2.7*2.7*3.1415*3.1415/(180.*180.))); beamR = 0.5*(2.*geo.DetHalfHeight()*cos(0.05) + geo.DetLength()*sin(0.05)); epsilon = beamR - geo.DetLength()*sin(0.05); cx = 0.; cy = epsilon*cos(0.05)-geo.DetHalfHeight(); cz = -epsilon*sin(0.05); } TVector3 bspot(cx*1.e-2, cy*1.e-2, cz*1.e-2); genie::flux::GCylindTH1Flux* histFlux = new genie::flux::GCylindTH1Flux(); histFlux->SetNuDirection(bdir); histFlux->SetBeamSpot(bspot); histFlux->SetTransverseRadius(beamR*1.e-2); std::cout << "beam r = " << beamR << " centered at (" << cx << "," << cy << "," << cz << ")" << std::endl; for(unsigned int i = 0; i < fGenFlavors.size(); ++i){ if(fFluxHistograms[fFlavorMap[fGenFlavors[i]]]) histFlux->AddEnergySpectrum(fGenFlavors[i], fFluxHistograms[fFlavorMap[fGenFlavors[i]]]); } fFluxD = dynamic_cast<genie::GFluxI *>(histFlux); }//end if using a histogram else if(fFluxType.compare("mono") == 0){ double weight = 1./(1.*fGenFlavors.size()); //make a map of pdg to weight codes std::map<int, double> pdgwmap; for(unsigned int i = 0; i < fGenFlavors.size(); ++i) pdgwmap[fGenFlavors[i]] = weight; genie::flux::GMonoEnergeticFlux *monoflux = new genie::flux::GMonoEnergeticFlux(fMonoEnergy, pdgwmap); fFluxD = dynamic_cast<genie::GFluxI *>(monoflux); }//end if using monoenergetic beam return; }

void GENIEHelper::InitializeGeometry (   )  [private]

get the geometry of the detector the detector geometry uses cgs units. casting to the GENIE geometry driver interface Definition at line 214 of file GENIEHelper.cxx. References fDetector, fGeomD, geo::Geometry::Instance(), and geo::Geometry::ROOTGeoManager(). Referenced by Initialize(). { geo::Geometry& geo = geo::Geometry::Instance(0, fDetector); genie::geometry::ROOTGeomAnalyzer *rgeom = new genie::geometry::ROOTGeomAnalyzer(geo.ROOTGeoManager()); rgeom->SetLengthUnits(genie::units::centimeter); rgeom->SetDensityUnits(genie::units::gram_centimeter3); fGeomD = dynamic_cast<genie::GeomAnalyzerI *>(rgeom); return; }

void GENIEHelper::PackMCTruth (  genie::EventRecord *  record, sim::MCTruth &  truth )  [private]

check the vertex z position against the length of the detector if the neutrino interacts beyond the end of the detector, we dont need to save it in the file get the Interaction object from the record - this is the object that talks to the event information objects and is in m get the different components making up the interaction set whether it is NC or CC set the interaction type add the particles from the interaction GHepParticles return units of GeV/c for p. the V_i are all in fermis and are relative to the center of the struck nucleus. add the vertex X/Y/Z to the V_i for status codes 0 and 1 set the vertex location for the neutrino, nucleus and everything that is to be tracked. vertex returns values in meters. end loop to convert GHepParticles to TParticles Definition at line 364 of file GENIEHelper.cxx. References sim::MCTruth::Add(), geo::Geometry::DetLength(), fDetector, geo::Geometry::Instance(), sim::MCTruth::SetCCNC(), sim::MCTruth::SetMode(), and sim::MCTruth::SetNeutrino(). Referenced by Sample(). { TLorentzVector *vertex = record->Vertex(); //std::cout << "vertex loc " << vertex->X() << "," << vertex->Y() << "," << vertex->Z() << std::endl; geo::Geometry& geo = geo::Geometry::Instance(0, fDetector); if(vertex->Z() > geo.DetLength()) return; genie::Interaction *inter = record->Summary(); const genie::InitialState &initState = inter->InitState(); const genie::ProcessInfo &procInfo = inter->ProcInfo(); const genie::Kinematics &kine = inter->Kine(); const genie::XclsTag &exclTag = inter->ExclTag(); const genie::KPhaseSpace &phaseSpace = inter->PhaseSpace(); truth.SetNeutrino(initState.ProbePdg(), initState.GetProbeP4()->Energy()); if(procInfo.IsWeakCC()) truth.SetCCNC(sim::kCC); else if(procInfo.IsWeakNC()) truth.SetCCNC(sim::kNC); if(procInfo.IsQuasiElastic()) truth.SetMode(sim::kQE); else if(procInfo.IsDeepInelastic()) truth.SetMode(sim::kDIS); else if(procInfo.IsResonant()) truth.SetMode(sim::kRes); else if(procInfo.IsCoherent()) truth.SetMode(sim::kCoh); TIter partitr(record); genie::GHepParticle *part = 0; while( (part = dynamic_cast<genie::GHepParticle *>(partitr.Next())) ){ TParticle tpart(part->Pdg(), part->Status(), part->FirstMother(), part->LastMother(), part->FirstDaughter(), part->LastDaughter(), part->Px(), part->Py(), part->Pz(), part->E(), part->Vx(), part->Vy(), part->Vz(), part->Vt()); if(part->Status() == 0 || part->Status() == 1) tpart.SetProductionVertex(100.*(part->Vx()*1.e-15 + vertex->X()), 100.*(part->Vy()*1.e-15 + vertex->Y()), 100.*(part->Vz()*1.e-15 + vertex->Z()), part->Vt()); truth.Add(tpart); } return; }

void GENIEHelper::PackNuMIFlux (  sim::MCTruth &  truth  )  [private]

cast the fFluxD pointer to be of the right type Definition at line 353 of file GENIEHelper.cxx. Referenced by Sample(). { genie::flux::GNuMIFlux *gnf = dynamic_cast<genie::flux::GNuMIFlux *>(fFluxD); const genie::flux::GNuMIFluxPassThroughInfo& nflux = gnf->PassThroughInfo(); return; }

bool GENIEHelper::Sample (  sim::MCTruth &  truth  )

pack the flux information make sure you havent gone over the limit set the flag in the parent object that says the fluxes came from histograms and fill related values get the flux for each neutrino flavor of this energy make sure you havent gone over the limit Definition at line 315 of file GENIEHelper.cxx. References sim::MCTruth::Enu(), fDriver, fFluxHistograms, fFluxType, fTotalGen, kNue, kNueBar, kNuMu, kNuMuBar, kNuTau, kNuTauBar, PackMCTruth(), PackNuMIFlux(), and sim::MCTruth::SetFluxGen(). Referenced by main(). { genie::EventRecord* record = fDriver->GenerateEvent(); if(!record) return false; PackMCTruth(record, truth); if(fFluxType.compare("ntuple") == 0) PackNuMIFlux(truth); else if(fFluxType.compare("histogram") == 0){ if(fTotalGen > fTotalRequested) return false; fTotalGen += 1.; int bin = fFluxHistograms[kNue]->FindBin(truth.Enu()); truth.SetFluxGen(fFluxHistograms[kNue]->GetBinContent(bin), fFluxHistograms[kNueBar]->GetBinContent(bin), fFluxHistograms[kNuMu]->GetBinContent(bin), fFluxHistograms[kNuMuBar]->GetBinContent(bin), fFluxHistograms[kNuTau]->GetBinContent(bin), fFluxHistograms[kNuTauBar]->GetBinContent(bin)); } else if(fFluxType.compare("mono") == 0){ if(fTotalGen > fTotalRequested) return false; fTotalGen += 1.; } return true; }

void GENIEHelper::Update (  const cfg::Config &  c  )  [virtual]

get a random number generator set the environment, the vector should come in pairs of variable name, then value figure out which detector you are using make the histograms get the detector mass and figure out how many events correspond to the requested POT if it is non-zero events = flux * pot * 10^-38 cm^2 (xsec) * mass detector (in kg) / carbon mass (in kg) integrate the flux from all species you are using end loop to add flux histograms to driver pick the number of events to generate out of a poisson distribution if doing the ipnd, account for the rock under the detector volume is now in cm^2, density of rock is 2.7 g/cm^3, convert to kg end if figuring out how many events given POT corresponds to Implements cfg::Observer. Definition at line 87 of file GENIEHelper.cxx. References geo::Geometry::DetLength(), fBeamName, fDetector, fEnvironment, fEvents, fFlavorMap, fFluxFile, fFluxHistograms, fFluxNormalization, fFluxType, fGenFlavors, fMonoEnergy, fPOT, fTotalRequested, util::RandomHandler::GetGenerator(), geo::Geometry::Instance(), and geo::Geometry::TotalMass(). { std::string geom; c("FluxNormalization").Get(fFluxNormalization); c("Detector").Get(geom); c("FluxType").Get(fFluxType); c("FluxFile").Get(fFluxFile); c("BeamName").Get(fBeamName); c("Events").Get(fEvents); c("POT").Get(fPOT); c("MonoEnergy").Get(fMonoEnergy); c("GenFlavors").Get(fGenFlavors); c("Environment").Get(fEnvironment); util::TRandomNOVA *rand = util::RandomHandler::GetGenerator(); TString junk = ""; junk += rand->Integer(100000000); std::string seed(junk); fEnvironment.push_back("GSEED"); fEnvironment.push_back(seed); for(unsigned int i = 0; i < fEnvironment.size(); i += 2){ gSystem->Setenv(fEnvironment[i].c_str(), fEnvironment[i+1].c_str()); std::cout << "setting GENIE environment " << fEnvironment[i] << " to " << fEnvironment[i+1] << std::endl; } fDetector = rawdata::kIPND; if(geom.compare("near") == 0) fDetector = rawdata::kNear; else if(geom.compare("far") == 0) fDetector = rawdata::kFar; if(fEvents > 0 && fPOT > 0.){ std::cerr << "cant specify both # of POT and # of events - using pot" << std::endl; } fTotalRequested = 1.*fEvents; if(fFluxType.compare("histogram") == 0){ TDirectory *savedir = gDirectory; fFluxHistograms.resize(6); TFile tf(fFluxFile.c_str()); tf.ls(); fFluxHistograms[kNue] = dynamic_cast<TH1D *>(tf.Get("nue")); fFluxHistograms[kNueBar] = dynamic_cast<TH1D *>(tf.Get("nuebar")); fFluxHistograms[kNuMu] = dynamic_cast<TH1D *>(tf.Get("numu")); fFluxHistograms[kNuMuBar] = dynamic_cast<TH1D *>(tf.Get("numubar")); fFluxHistograms[kNuTau] = dynamic_cast<TH1D *>(tf.Get("nutau")); fFluxHistograms[kNuTauBar] = dynamic_cast<TH1D *>(tf.Get("nutaubar")); for(unsigned int i = 0; i < fFluxHistograms.size(); ++i) if(fFluxHistograms[i]) fFluxHistograms[i]->SetDirectory(savedir); if(fPOT > 0.){ geo::Geometry& geo = geo::Geometry::Instance(0, fDetector); double mass = geo.TotalMass(); double flux = 0.; for(unsigned int i = 0; i < fGenFlavors.size(); ++i){ if(!fFluxHistograms[fFlavorMap[fGenFlavors[i]]]) continue; else flux += fFluxHistograms[fFlavorMap[fGenFlavors[i]]]->Integral(); } fTotalRequested = rand->Poisson(flux*fPOT*mass*5.02e-13*1.e-20); if(fDetector == rawdata::kIPND){ double length = geo.DetLength(); double vol = 3.141592654*(0.33*pow(length, 3.)*pow(sin(0.11)*cos(0.11), 2.) + pow(length, 3.)*pow(sin(0.11), 2.)*(pow(cos(0.11), 2.) - 1.)); double massrock = 2.7*vol*1.e-3; fTotalRequested = rand->Poisson(flux*fPOT*(mass+massrock)*5.02e-13*1.e-20); }//end if IPND } }//end if getting fluxes from histograms std::cout << "Generating events for " << fFluxType.c_str() << " from " << fFluxFile.c_str() << " with an exposure of " << fTotalRequested << "(events)/" << fPOT << "(POT) in detector " << geom.c_str() << std::endl; std::cout << "Generating the following flavors " << std::endl; for(unsigned int i = 0; i < fGenFlavors.size(); ++i) std::cout << " " << fGenFlavors[i] << std::endl; if(fFluxType.compare("mono")==0) std::cout << "Generating monoenergetic (" << fMonoEnergy << " GeV) neutrinos" << std::endl; return; }

---

Member Data Documentation

std::string evgen::GENIEHelper::fBeamName [private]

name of file containing histograms or ntuples, or txt Definition at line 57 of file GENIEHelper.h. Referenced by Update().

rawdata::det_id_ evgen::GENIEHelper::fDetector [private]

Definition at line 54 of file GENIEHelper.h. Referenced by InitializeFluxDriver(), InitializeGeometry(), PackMCTruth(), and Update().

genie::GMCJDriver* evgen::GENIEHelper::fDriver [private]

Definition at line 52 of file GENIEHelper.h. Referenced by Initialize(), and Sample().

std::vector<string> evgen::GENIEHelper::fEnvironment [private]

pdg codes for flavors to generate Definition at line 67 of file GENIEHelper.h. Referenced by Update().

int evgen::GENIEHelper::fEvents [private]

histograms for each nu species Definition at line 60 of file GENIEHelper.h. Referenced by Update().

std::map<int, int> evgen::GENIEHelper::fFlavorMap [private]

environmental variables and settings used by genie Definition at line 68 of file GENIEHelper.h. Referenced by GENIEHelper(), InitializeFluxDriver(), and Update().

genie::GFluxI* evgen::GENIEHelper::fFluxD [private]

Definition at line 51 of file GENIEHelper.h. Referenced by Initialize(), and InitializeFluxDriver().

std::string evgen::GENIEHelper::fFluxFile [private]

histogram or ntuple Definition at line 56 of file GENIEHelper.h. Referenced by InitializeFluxDriver(), and Update().

std::vector<TH1D *> evgen::GENIEHelper::fFluxHistograms [private]

name of the beam we are simulating Definition at line 58 of file GENIEHelper.h. Referenced by InitializeFluxDriver(), Sample(), and Update().

double evgen::GENIEHelper::fFluxNormalization [private]

total POT to generate Definition at line 64 of file GENIEHelper.h. Referenced by InitializeFluxDriver(), and Update().

std::string evgen::GENIEHelper::fFluxType [private]

ipnd, nd, fd Definition at line 55 of file GENIEHelper.h. Referenced by InitializeFluxDriver(), Sample(), and Update().

std::vector<int> evgen::GENIEHelper::fGenFlavors [private]

energy of monoenergetic neutrinos Definition at line 66 of file GENIEHelper.h. Referenced by InitializeFluxDriver(), and Update().

genie::GeomAnalyzerI* evgen::GENIEHelper::fGeomD [private]

Definition at line 50 of file GENIEHelper.h. Referenced by Initialize(), and InitializeGeometry().

double evgen::GENIEHelper::fMonoEnergy [private]

POT per file. Definition at line 65 of file GENIEHelper.h. Referenced by InitializeFluxDriver(), and Update().

double evgen::GENIEHelper::fPOT [private]

total number of events/pot requested to be generated Definition at line 63 of file GENIEHelper.h. Referenced by InitializeFluxDriver(), and Update().

double evgen::GENIEHelper::fTotalGen [private]

number of events to generate Definition at line 61 of file GENIEHelper.h. Referenced by Sample().

double evgen::GENIEHelper::fTotalRequested [private]

total number of events/pot generated to this point Definition at line 62 of file GENIEHelper.h. Referenced by Update().

---

The documentation for this class was generated from the following files:

• GENIEHelper.h
• GENIEHelper.cxx

---