104 #include <TVector3.h> 106 #include <TStopwatch.h> 130 using namespace genie;
141 if(fUnphysEventMask)
delete fUnphysEventMask;
142 if (fGPool)
delete fGPool;
144 map<int,TH1D*>::iterator pmax_iter = fPmax.begin();
145 for( ; pmax_iter != fPmax.end(); ++pmax_iter) {
146 TH1D *
pmax = pmax_iter->second;
148 delete pmax; pmax = 0;
153 if(fFluxIntTree)
delete fFluxIntTree;
154 if(fFluxIntProbFile)
delete fFluxIntProbFile;
160 <<
"Setting event generator list: " << listname;
162 fEventGenList = listname;
167 *fUnphysEventMask =
mask;
171 <<
" -> 0) : " << *fUnphysEventMask;
176 fFluxDriver = flux_driver;
181 fGeomAnalyzer = geom_analyzer;
197 fMaxPlXmlFilename = xml_filename;
199 bool is_accessible = !(gSystem->AccessPathName(fMaxPlXmlFilename.c_str()));
201 if ( is_accessible ) fUseExtMaxPl =
true;
203 fUseExtMaxPl =
false;
205 <<
"UseMaxPathLengths could not find file: \"" << xml_filename <<
"\"";
214 <<
"Keep on throwing flux neutrinos till one interacts? : " 216 fKeepThrowingFluxNu = keep_on;
224 fGenerateUnweighted =
true;
227 <<
"GMCJDriver will generate un-weighted events. " 228 <<
"Note: That does not force unweighted event kinematics!";
236 fPreSelect = preselect;
250 bool save_to_file = fFluxIntProbFile == 0 && fFluxIntFileName.size()>0;
253 fSumFluxIntProbs.clear();
258 "Skipping pre-generation of flux interaction probabilities - "<<
259 "using pre-generated file";
266 fFluxIntProbFile =
new TFile(fFluxIntFileName.c_str(),
"CREATE");
267 if(fFluxIntProbFile->IsZombie()){
268 LOG(
"GMCJDriver",
pFATAL) <<
"Cannot overwrite an existing file. Exiting!";
274 fFluxIntTree =
new TTree(fFluxIntTreeName.c_str(),
275 "Tree storing pre-calculated flux interaction probs");
276 fFluxIntTree->Branch(
"FluxIndex", &fBrFluxIndex,
"FluxIndex/I");
277 fFluxIntTree->Branch(
"FluxIntProb", &fBrFluxIntProb,
"FluxIntProb/D");
278 fFluxIntTree->Branch(
"FluxEnu", &fBrFluxEnu,
"FluxEnu/D");
279 fFluxIntTree->Branch(
"FluxWeight", &fBrFluxWeight,
"FluxWeight/D");
280 fFluxIntTree->Branch(
"FluxPDG", &fBrFluxPDG,
"FluxPDG/I");
282 if(save_to_file) fFluxIntTree->SetDirectory(fFluxIntProbFile);
284 fFluxDriver->GenerateWeighted(
true);
289 TStopwatch stopwatch;
291 long int first_index = -1;
292 bool first_loop =
true;
294 while(fFluxDriver->End() ==
false){
297 bool gotnext = fFluxDriver->GenerateNext();
299 LOG(
"GMCJDriver",
pWARN) <<
"*** Couldn't generate next flux ray! ";
305 bool already_been_here = first_loop ?
false : first_index == fFluxDriver->Index();
306 if(already_been_here)
break;
309 if(this->ComputePathLengths() ==
false){ success =
false;
break;}
312 double psum = this->ComputeInteractionProbabilities(
false );
314 fBrFluxIntProb = psum;
315 fBrFluxIndex = fFluxDriver->Index();
316 fBrFluxEnu = fFluxDriver->Momentum().E();
317 fBrFluxWeight = fFluxDriver->Weight();
318 fBrFluxPDG = fFluxDriver->PdgCode();
319 fFluxIntTree->Fill();
323 first_index = fFluxDriver->Index();
329 <<
"Finished pre-calculating flux interaction probabilities. " 330 <<
"Total CPU time to process "<< fFluxIntTree->GetEntries()
331 <<
" entries: "<< stopwatch.CpuTime();
335 fFluxDriver->Clear(
"CycleHistory");
342 double safety_factor = 1.01;
343 for(
int i = 0;
i< fFluxIntTree->GetEntries();
i++){
344 fFluxIntTree->GetEntry(
i);
349 fGlobPmax = TMath::Max(fGlobPmax, fBrFluxIntProb*safety_factor);
351 if(fSumFluxIntProbs.find(fBrFluxPDG) == fSumFluxIntProbs.end()){
352 fSumFluxIntProbs[fBrFluxPDG] = 0.0;
354 fSumFluxIntProbs[fBrFluxPDG] += fBrFluxIntProb * fBrFluxWeight;
357 "Updated global probability scale to fGlobPmax = "<< fGlobPmax;
361 "Saving pre-generated interaction probabilities to file: "<<
362 fFluxIntProbFile->GetName();
363 fFluxIntProbFile->cd();
364 fFluxIntTree->Write();
368 if(fFluxIntTree->BuildIndex(
"FluxIndex") != fFluxIntTree->GetEntries()){
370 "Cannot build index using branch \"FluxIndex\" for flux prob tree!";
376 this->PreSelectEvents(
false);
378 LOG(
"GMCJDriver",
pNOTICE) <<
"Successfully generated/loaded pre-calculate flux interaction probabilities";
381 else if(fFluxIntTree){
398 if(fFluxIntProbFile){
400 <<
"Can't load flux interaction prob file as one is already loaded";
404 fFluxIntProbFile =
new TFile(filename.c_str(),
"OPEN");
406 if(fFluxIntProbFile){
407 fFluxIntTree =
dynamic_cast<TTree*
>(fFluxIntProbFile->Get(fFluxIntTreeName.c_str()));
410 fFluxIntTree->SetBranchAddress(
"FluxIntProb", &fBrFluxIntProb) >= 0 &&
411 fFluxIntTree->SetBranchAddress(
"FluxIndex", &fBrFluxIndex) >= 0 &&
412 fFluxIntTree->SetBranchAddress(
"FluxPDG", &fBrFluxPDG) >= 0 &&
413 fFluxIntTree->SetBranchAddress(
"FluxWeight", &fBrFluxWeight) >= 0 &&
414 fFluxIntTree->SetBranchAddress(
"FluxEnu", &fBrFluxEnu) >= 0;
417 if(this->PreCalcFluxProbabilities()) {
419 <<
"Successfully loaded pre-generated flux interaction probabilities";
425 "Cannot find expected branches in input flux probability tree!";
426 delete fFluxIntTree; fFluxIntTree = 0;
429 <<
"Cannot find tree: "<< fFluxIntTreeName.c_str();
433 <<
"Unable to load flux interaction probabilities file";
453 this->GetParticleLists();
456 this->GetMaxFluxEnergy();
463 this->PopulateEventGenDriverPool();
475 this->BootstrapXSecSplines();
480 this->BootstrapXSecSplineSummation();
482 if(calc_prob_scales){
485 this->GetMaxPathLengthList();
489 this->ComputeProbScales();
491 LOG(
"GMCJDriver",
pNOTICE) <<
"Finished configuring GMCJDriver\n\n";
496 fEventGenList =
"Default";
501 fUnphysEventMask->SetBitNumber(
i,
true);
508 fMaxPlXmlFilename =
"";
509 fUseExtMaxPl =
false;
516 fGenerateUnweighted =
false;
521 fCurVtx.SetXYZT(0.,0.,0.,0.);
523 fFluxIntProbFile = 0;
524 fFluxIntTreeName =
"gFlxIntProb";
525 fFluxIntFileName =
"";
527 fBrFluxIntProb = -1.;
532 fSumFluxIntProbs.clear();
536 this->KeepOnThrowingFluxNeutrinos(
true);
556 fMaxPathLengths.clear();
557 fCurPathLengths.clear();
564 <<
"Asking the flux driver for its list of neutrinos";
565 fNuList = fFluxDriver->FluxParticles();
567 LOG(
"GMCJDriver",
pNOTICE) <<
"Flux particles: " << fNuList;
571 <<
"Asking the geometry driver for its list of targets";
572 fTgtList = fGeomAnalyzer->ListOfTargetNuclei();
574 LOG(
"GMCJDriver",
pNOTICE) <<
"Target materials: " << fTgtList;
581 <<
"Loading external max path-length list for input geometry from " 582 << fMaxPlXmlFilename;
583 fMaxPathLengths.LoadFromXml(fMaxPlXmlFilename);
587 <<
"Querying the geometry driver to compute the max path-length list";
588 fMaxPathLengths = fGeomAnalyzer->ComputeMaxPathLengths();
592 <<
"Maximum path length list: " << fMaxPathLengths;
598 <<
"Querying the flux driver for the maximum energy of flux neutrinos";
599 fEmax = fFluxDriver->MaxEnergy();
602 <<
"Maximum flux neutrino energy = " << fEmax <<
" GeV";
608 <<
"Creating GEVGPool & adding a GEVGDriver object per init-state";
610 if (fGPool)
delete fGPool;
613 PDGCodeList::const_iterator nuiter;
614 PDGCodeList::const_iterator tgtiter;
616 for(nuiter = fNuList.begin(); nuiter != fNuList.end(); ++nuiter) {
617 for(tgtiter = fTgtList.begin(); tgtiter != fTgtList.end(); ++tgtiter) {
619 int target_pdgc = *tgtiter;
620 int neutrino_pdgc = *nuiter;
625 <<
"\n\n ---- Creating a GEVGDriver object configured for init-state: " 626 << init_state.
AsString() <<
" ----\n\n";
633 LOG(
"GMCJDriver",
pDEBUG) <<
"Adding new GEVGDriver object to GEVGPool";
634 fGPool->insert( GEVGPool::value_type(init_state.
AsString(), evgdriver) );
639 <<
"All necessary GEVGDriver object were pushed into GEVGPool\n";
647 if(!fUseSplines)
return;
650 <<
"Asking event generation drivers to compute all needed xsec splines";
652 PDGCodeList::const_iterator nuiter;
653 PDGCodeList::const_iterator tgtiter;
654 for(nuiter = fNuList.begin(); nuiter != fNuList.end(); ++nuiter){
655 for(tgtiter = fTgtList.begin(); tgtiter != fTgtList.end(); ++tgtiter) {
656 int target_pdgc = *tgtiter;
657 int neutrino_pdgc = *nuiter;
660 <<
"Computing all splines needed for init-state: " 662 GEVGDriver * evgdriver = fGPool->FindDriver(init_state);
666 LOG(
"GMCJDriver",
pINFO) <<
"Finished creating cross section splines\n";
675 <<
"Summing-up splines to get total cross section for each init state";
677 GEVGPool::iterator diter;
678 for(diter = fGPool->begin(); diter != fGPool->end(); ++diter) {
679 string init_state = diter->first;
683 <<
"**** Summing xsec splines for init-state = " << init_state;
686 if (fEmax>rE.
max || fEmax<rE.
min)
688 <<
" rE (validEnergyRange) [" << rE.
min <<
"," << rE.
max <<
"] " 689 <<
" fEmax " << fEmax;
696 double dE = fEmax/10.;
698 double max = (fEmax+dE < rE.
max) ? fEmax+dE : rE.
max;
702 <<
"Finished summing all interaction xsec splines per initial state";
719 <<
"Computing the max. interaction probability (probability scale)";
725 map<int,TH1D*>::iterator pmax_iter = fPmax.begin();
726 for( ; pmax_iter != fPmax.end(); ++pmax_iter) {
727 TH1D *
pmax = pmax_iter->second;
729 delete pmax; pmax = 0;
738 double de = fEmax/300.;
740 double emax = fEmax + de;
741 int n = 1 + (
int) ((emax-emin)/de);
743 PDGCodeList::const_iterator nuiter;
744 PDGCodeList::const_iterator tgtiter;
747 for(nuiter = fNuList.begin(); nuiter != fNuList.end(); ++nuiter) {
748 int neutrino_pdgc = *nuiter;
749 TH1D * pmax_hst =
new TH1D(
"pmax_hst",
750 "max interaction probability vs E | geom",n,emin,emax);
751 pmax_hst->SetDirectory(0);
754 for(
int ie = 1; ie <= pmax_hst->GetNbinsX(); ie++) {
755 double EvLow = pmax_hst->GetBinCenter(ie) - 0.5*pmax_hst->GetBinWidth(ie);
756 double EvHigh = pmax_hst->GetBinCenter(ie) + 0.5*pmax_hst->GetBinWidth(ie);
762 for(tgtiter = fTgtList.begin(); tgtiter != fTgtList.end(); ++tgtiter) {
763 int target_pdgc = *tgtiter;
768 <<
"Computing Pmax for init-state: " << init_state.
AsString() <<
" E from " << EvLow <<
"-" << EvHigh;
771 GEVGDriver * evgdriver = fGPool->FindDriver(init_state);
778 double plmax = fMaxPathLengths.PathLength(target_pdgc);
782 double pmaxLow = this->InteractionProbability(sxsecLow, plmax, A);
783 double pmaxHigh = this->InteractionProbability(sxsecHigh, plmax, A);
785 double pmax = pmaxHigh;
786 if ( pmaxLow > pmaxHigh){
789 <<
"Lower energy neutrinos have a higher probability of interacting than those at higher energy." 790 <<
" pmaxLow(E=" << EvLow <<
")=" << pmaxLow <<
" and " <<
" pmaxHigh(E=" << EvHigh <<
")=" << pmaxHigh;
793 pmax_hst->SetBinContent(ie, pmax_hst->GetBinContent(ie) +
pmax);
796 <<
"Pmax[" << init_state.
AsString() <<
", Ev from " << EvLow <<
"-" << EvHigh <<
"] = " <<
pmax;
799 pmax_hst->SetBinContent(ie, 1.2 * pmax_hst->GetBinContent(ie));
802 <<
"Pmax[nu=" << neutrino_pdgc <<
", Ev from " << EvLow <<
"-" << EvHigh <<
"] = " 803 << pmax_hst->GetBinContent(ie);
806 fPmax.insert(map<int,TH1D*>::value_type(neutrino_pdgc,pmax_hst));
814 for(nuiter = fNuList.begin(); nuiter != fNuList.end(); ++nuiter) {
815 int neutrino_pdgc = *nuiter;
816 map<int,TH1D*>::const_iterator pmax_iter = fPmax.find(neutrino_pdgc);
817 assert(pmax_iter != fPmax.end());
818 TH1D * pmax_hst = pmax_iter->second;
821 double pmax = pmax_hst->GetMaximum();
824 fGlobPmax = TMath::Max(pmax, fGlobPmax);
826 LOG(
"GMCJDriver",
pNOTICE) <<
"*** Probability scale = " << fGlobPmax;
832 fCurPathLengths.clear();
835 fCurVtx.SetXYZT(0.,0.,0.,0.);
840 LOG(
"GMCJDriver",
pNOTICE) <<
"Generating next event...";
842 this->InitEventGeneration();
845 bool flux_end = fFluxDriver->End();
848 <<
"No more neutrinos can be thrown by the flux driver";
853 if(event)
return event;
855 if(fKeepThrowingFluxNu) {
857 <<
"Flux neutrino didn't interact - Trying the next one...";
863 LOG(
"GMCJDriver",
pINFO) <<
"Returning NULL event!";
873 double Pno=0, Psum=0;
874 double R = rnd->
RndEvg().Rndm();
875 LOG(
"GMCJDriver",
pDEBUG) <<
"Rndm [0,1] = " <<
R;
878 bool flux_ok = this->GenerateFluxNeutrino();
881 <<
"** Rejecting current flux neutrino (flux driver err)";
893 <<
"Computing interaction probabilities for max. path lengths";
895 Psum = this->ComputeInteractionProbabilities(
true );
898 <<
"The no-interaction probability (max. path lengths) is: " 902 <<
"Negative no-interaction probability! (P = " << 100*Pno <<
" %)" 903 <<
" Particle E=" << fFluxDriver->Momentum().E() <<
" type=" << fFluxDriver->PdgCode() <<
"Psum=" << Psum;
909 <<
"** Rejecting current flux neutrino";
919 Psum = this->PreGenFluxInteractionProbability();
925 pl_ok = this->ComputePathLengths();
928 <<
"** Rejecting current flux neutrino (err computing path-lengths)";
931 if(fCurPathLengths.AreAllZero()) {
933 <<
"** Rejecting current flux neutrino (misses generation volume)";
936 Psum = this->ComputeInteractionProbabilities(
false );
942 <<
"** Rejecting current flux neutrino (has null interaction probability)";
949 <<
"The actual 'no interaction' probability is: " << 100*Pno <<
" %";
952 <<
"Negative no interactin probability! (P = " << 100*Pno <<
" %)";
955 int nupdg = fFluxDriver ->
PdgCode ();
956 const TLorentzVector & nup4 = fFluxDriver -> Momentum ();
957 const TLorentzVector & nux4 = fFluxDriver -> Position ();
960 <<
"\n [-] Problematic neutrino: " 961 <<
"\n |----o PDG-code : " << nupdg
964 <<
"\n Emax : " << fEmax;
967 <<
"\n Problematic path lengths:" << fCurPathLengths;
970 <<
"\n Maximum path lengths:" << fMaxPathLengths;
976 <<
"** Rejecting current flux neutrino";
987 pl_ok = this->ComputePathLengths();
989 LOG(
"GMCJDriver",
pFATAL) <<
"** Cannot calculate path lenths!";
992 double Psum_curr = this->ComputeInteractionProbabilities(
false );
996 "** Mismatch between pre-calculated and current interaction "<<
1003 fSelTgtPdg = this->SelectTargetMaterial(R);
1006 <<
"** Rejecting current flux neutrino (failed to select tgt!)";
1012 this->GenerateEventKinematics();
1015 <<
"** Couldn't generate kinematics for selected interaction";
1021 this->GenerateVertexPosition();
1030 this->ComputeEventProbability();
1040 LOG(
"GMCJDriver",
pNOTICE) <<
"Generating a flux neutrino";
1042 bool ok = fFluxDriver->GenerateNext();
1045 <<
"*** The flux driver couldn't generate a flux neutrino!!";
1050 int nupdg = fFluxDriver ->
PdgCode ();
1051 const TLorentzVector & nup4 = fFluxDriver -> Momentum ();
1052 const TLorentzVector & nux4 = fFluxDriver -> Position ();
1055 <<
"\n [-] Generated flux neutrino: " 1056 <<
"\n |----o PDG-code : " << nupdg
1060 if(nup4.Energy() > fEmax) {
1062 <<
"\n *** Flux driver error ***" 1063 <<
"\n Generated flux v with E = " << nup4.Energy() <<
" GeV" 1064 <<
"\n Max v energy (declared by flux driver) = " << fEmax <<
" GeV" 1065 <<
"\n My interaction probability scaling is invalidated!!";
1068 if(!fNuList.ExistsInPDGCodeList(nupdg)) {
1070 <<
"\n *** Flux driver error ***" 1071 <<
"\n Generated flux v with pdg = " << nupdg
1072 <<
"\n It does not belong to the declared list of flux neutrinos" 1073 <<
"\n I was not configured to handle this!!";
1085 fCurPathLengths.clear();
1087 const TLorentzVector & nup4 = fFluxDriver -> Momentum ();
1088 const TLorentzVector & nux4 = fFluxDriver -> Position ();
1090 fCurPathLengths = fGeomAnalyzer->ComputePathLengths(nux4, nup4);
1092 LOG(
"GMCJDriver",
pNOTICE) << fCurPathLengths;
1094 if(fCurPathLengths.size() == 0) {
1096 <<
"\n *** Geometry driver error ***" 1097 <<
"\n Got an empty PathLengthList - No material found in geometry?";
1101 if(fCurPathLengths.AreAllZero()) {
1103 <<
"current flux v doesn't cross any geometry material...";
1111 <<
"Computing relative interaction probabilities for each material";
1114 int nupdg = fFluxDriver->PdgCode();
1115 const TLorentzVector & nup4 = fFluxDriver->Momentum();
1117 fCurCumulProbMap.clear();
1120 (use_max_path_length) ? fMaxPathLengths : fCurPathLengths;
1123 PathLengthList::const_iterator pliter;
1125 for(pliter = path_length_list.begin();
1126 pliter != path_length_list.end(); ++pliter) {
1127 int mpdg = pliter->first;
1128 double pl = pliter->second;
1136 GEVGDriver * evgdriver = fGPool->FindDriver(init_state);
1139 <<
"\n * The MC Job driver isn't properly configured!" 1140 <<
"\n * No event generation driver could be found for init state: " 1149 <<
"\n * The MC Job driver isn't properly configured!" 1150 <<
"\n * Couldn't retrieve total cross section spline for init state: " 1154 xsec = totxsecspl->
Evaluate( nup4.Energy() );
1156 prob = this->InteractionProbability(xsec,pl,A);
1158 <<
" (xsec, pl, A)=(" << xsec <<
"," << pl <<
"," << A <<
")";
1164 if(fGenerateUnweighted) pmax = fGlobPmax;
1166 map<int,TH1D*>::const_iterator pmax_iter = fPmax.find(nupdg);
1167 assert(pmax_iter != fPmax.end());
1168 TH1D * pmax_hst = pmax_iter->second;
1170 int ie = pmax_hst->FindBin(nup4.Energy());
1171 pmax = pmax_hst->GetBinContent(ie);
1178 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 1180 <<
"tgt: " << mpdg <<
" -> TotXSec = " 1181 << xsec/
units::cm2 <<
" cm^2, Norm.Prob = " << 100*probn <<
"%";
1185 fCurCumulProbMap.insert(map<int,double>::value_type(mpdg,probsum));
1195 LOG(
"GMCJDriver",
pNOTICE) <<
"Selecting target material";
1197 map<int,double>::const_iterator probiter = fCurCumulProbMap.begin();
1198 for( ; probiter != fCurCumulProbMap.end(); ++probiter) {
1199 double prob = probiter->second;
1201 tgtpdg = probiter->first;
1203 <<
"Selected target material = " << tgtpdg;
1208 <<
"Could not select target material for an interacting neutrino";
1214 int nupdg = fFluxDriver->PdgCode();
1215 const TLorentzVector & nup4 = fFluxDriver->Momentum();
1220 GEVGDriver * evgdriver = fGPool->FindDriver(init_state);
1223 <<
"No GEVGDriver object for init state: " << init_state.
AsString();
1233 <<
"Asking the selected GEVGDriver object to generate an event";
1242 <<
"Asking the geometry analyzer to generate a vertex";
1244 const TLorentzVector &
p4 = fFluxDriver->Momentum ();
1245 const TLorentzVector & x4 = fFluxDriver->Position ();
1247 const TVector3 & vtx = fGeomAnalyzer->GenerateVertex(x4, p4, fSelTgtPdg);
1249 TVector3 origin(x4.X(), x4.Y(), x4.Z());
1252 double dL = origin.Mag();
1257 <<
"|vtx - origin|: dL = " << dL <<
" m, dt = " << dt <<
" sec";
1259 fCurVtx.SetXYZT(vtx.x(), vtx.y(), vtx.z(), x4.T() +
dt);
1261 fCurEvt->SetVertex(fCurVtx);
1269 double xsec = fCurEvt->XSec();
1272 PathLengthList::const_iterator pliter = fCurPathLengths.find(fSelTgtPdg);
1273 double path_length = pliter->second;
1279 double P = this->InteractionProbability(xsec, path_length, A);
1286 int nu_pdg = nu->
Pdg();
1287 double Ev = nu->
P4()->Energy();
1290 if(!fGenerateUnweighted) {
1291 map<int,TH1D*>::const_iterator pmax_iter = fPmax.find(nu_pdg);
1292 assert(pmax_iter != fPmax.end());
1293 TH1D * pmax_hst = pmax_iter->second;
1295 double pmax = pmax_hst->GetBinContent(pmax_hst->FindBin(Ev));
1297 weight = pmax/fGlobPmax;
1301 fCurEvt->SetProbability(P);
1302 fCurEvt->SetWeight(weight * fCurEvt->Weight());
1315 return kNA*(xsec*pL)/A;
1326 "Cannot get pre-computed flux interaction probability as no tree!";
1330 assert(fFluxDriver->Index() >= 0);
1334 bool found_entry = fFluxIntTree->GetEntryWithIndex(fFluxDriver->Index()) > 0;
1335 bool enu_match =
false;
1337 double rel_err = fBrFluxEnu-fFluxDriver->Momentum().E();
1340 if(enu_match ==
false){
1342 "Mismatch between: Enu_curr = "<< fFluxDriver->Momentum().E() <<
1343 ", Enu_pre_gen = "<< fBrFluxEnu;
1347 LOG(
"GMCJDriver",
pERROR) <<
"Cannot find flux entry in interaction prob tree!";
1351 bool success = found_entry && enu_match;
1354 "Cannot find pre-generated interaction probability! Check you "<<
1355 "are using the correct pre-generated interaction prob file " <<
1356 "generated using current flux input file with same input " <<
1357 "config (same geom TopVol, neutrino species list)";
1361 return fBrFluxIntProb/fGlobPmax;
static constexpr double dL
static const double second
void CreateXSecSumSpline(int nk, double Emin, double Emax, bool inlogE=true)
void SetUnphysEventMask(const TBits &mask)
void PopulateEventGenDriverPool(void)
bool PreCalcFluxProbabilities(void)
THE MAIN GENIE PROJECT NAMESPACE
void GetMaxPathLengthList(void)
static RandomGen * Instance()
Access instance.
void SetEventGeneratorList(string listname)
const TLorentzVector * P4(void) const
void InitEventGeneration(void)
A simple [min,max] interval for doubles.
string BoolAsYNString(bool b)
A numeric analysis tool class for interpolating 1-D functions.
int IonPdgCodeToA(int pdgc)
bool GenerateFluxNeutrino(void)
bool ComputePathLengths(void)
string P4AsString(const TLorentzVector *p)
const Spline * XSecSumSpline(void) const
double Evaluate(double x) const
A singleton holding random number generator classes. All random number generation in GENIE should tak...
void ComputeEventProbability(void)
Range1D_t ValidEnergyRange(void) const
void UseFluxDriver(GFluxI *flux)
#define P(a, b, c, d, e, x)
Object to be filled with the neutrino path-length, for all detector geometry materials, when starting from a position x and travelling along the direction of the neutrino 4-momentum.
void KeepOnThrowingFluxNeutrinos(bool keep_on)
void GetMaxFluxEnergy(void)
string outfilename
knobs that need extra care
static unsigned int NFlags(void)
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
void SetEventGeneratorList(string listname)
void ForceSingleProbScale(void)
double PreGenFluxInteractionProbability(void)
static Messenger * Instance(void)
bool UseMaxPathLengths(string xml_filename)
void Configure(bool calc_prob_scales=true)
static const double meter
string AsString(void) const
static const double kilogram
static const double kASmallNum
TRandom3 & RndEvg(void) const
rnd number generator used by the event generation drivers
GENIE Event Generation Driver. A minimalist user interface object for generating neutrino interaction...
void BootstrapXSecSplines(void)
double ComputeInteractionProbabilities(bool use_max_path_length)
static const double kLightSpeed
EventRecord * GenerateEvent(void)
static float min(const float a, const float b, const float c)
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
double InteractionProbability(double xsec, double pl, int A)
void PreSelectEvents(bool preselect=true)
void GenerateEventKinematics(void)
int SelectTargetMaterial(double R)
void GenerateVertexPosition(void)
void ComputeProbScales(void)
void Configure(int nu_pdgc, int Z, int A)
bool LoadFluxProbabilities(string filename)
EventRecord * GenerateEvent1Try(void)
assert(nhit_max >=nhit_nbins)
string PrintFramedMesg(string mesg, unsigned int nl=1, const char f='*')
void CreateSplines(int nknots=-1, double emax=-1, bool inLogE=true)
void BootstrapXSecSplineSummation(void)
Defines the GENIE Geometry Analyzer Interface.
string X4AsString(const TLorentzVector *x)
void SetUnphysEventMask(const TBits &mask)
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
STDHEP-like event record entry that can fit a particle or a nucleus.
void UseGeomAnalyzer(GeomAnalyzerI *geom)
void UseSplines(bool useLogE=true)
void SaveFluxProbabilities(string outfilename)
void GetParticleLists(void)
static AlgConfigPool * Instance()
Initial State information.
GENIE Interface for user-defined flux classes.
EventRecord * GenerateEvent(const TLorentzVector &nu4p)
A pool of GEVGDriver objects with an initial state key.