evgb::GiBUUHelper Class Reference

#include "GiBUUHelper.h"

Inheritance diagram for evgb::GiBUUHelper:

Public Member Functions
	GiBUUHelper (fhicl::ParameterSet const &pset, TGeoManager *rootGeom, std::string const &rootFile, double const &detectorMass)
	~GiBUUHelper ()
genie::EventRecord *	GetGiBUUEventRecord ()
void	Initialize ()
bool	Stop ()
bool	Sample (simb::MCTruth &truth, simb::MCFlux &flux, simb::GTruth &gtruth)
double	TotalHistFlux ()
double	TotalExposure () const
double	SpillExposure () const
std::string	FluxType () const
std::string	DetectorLocation () const
std::vector< TH1D * >	FluxHistograms () const
double	TotalMass () const
genie::EventRecord *	GetGenieEventRecord ()
TRandom3 *	GetHelperRandom ()
genie::GFluxI *	GetFluxDriver (bool base=true)
std::string	GetTuneName () const
std::string	GetEventGeneratorList () const

Private Attributes
genie::EventRecord *	fGiBUUEventRecord

Detailed Description

Definition at line 20 of file GiBUUHelper.h.

Constructor & Destructor Documentation

evgb::GiBUUHelper::GiBUUHelper ( fhicl::ParameterSet const &  pset, TGeoManager *  rootGeom, std::string const &  rootFile, double const &  detectorMass  )

explicit

Definition at line 54 of file GiBUUHelper.cxx.

    : GENIEHelper(pset, geoManager, rootFile, detectorMass),
      fGiBUUEventRecord(0)
  {
  }

evgb::GiBUUHelper::~GiBUUHelper

(

)

Definition at line 64 of file GiBUUHelper.cxx.

References fGiBUUEventRecord.

  {
    // clean up owned genie object (other genie obj are ref ptrs)
    delete fGiBUUEventRecord;

  }

Member Function Documentation

std::string evgb::GENIEHelper::DetectorLocation ( ) const

inlineinherited

Definition at line 77 of file GENIEHelper.h.

{ return fDetLocation;    }

std::vector<TH1D*> evgb::GENIEHelper::FluxHistograms ( ) const

inlineinherited

Definition at line 81 of file GENIEHelper.h.

{ return fFluxHistograms; }

std::string evgb::GENIEHelper::FluxType ( ) const

inlineinherited

Definition at line 76 of file GENIEHelper.h.

{ return fFluxType;       }

std::string evgb::GENIEHelper::GetEventGeneratorList ( ) const

inlineinherited

Definition at line 96 of file GENIEHelper.h.

{ return fEventGeneratorList; }

genie::GFluxI* evgb::GENIEHelper::GetFluxDriver ( bool  base = true )

inlineinherited

Definition at line 91 of file GENIEHelper.h.

Referenced by evgen::TestGENIEHelper::produce(), and evgen::GENIEGen::produce().

      { return ( (base) ? fFluxD : fFluxD2GMCJD ); }

genie::EventRecord* evgb::GENIEHelper::GetGenieEventRecord ( )

inlineinherited

Definition at line 84 of file GENIEHelper.h.

Referenced by evgen::TestGENIEHelper::produce().

{ return fGenieEventRecord; }

genie::EventRecord* evgb::GiBUUHelper::GetGiBUUEventRecord ( )

inline

Definition at line 30 of file GiBUUHelper.h.

References fGiBUUEventRecord.

{ return fGiBUUEventRecord; }

TRandom3* evgb::GENIEHelper::GetHelperRandom ( )

inlineinherited

Definition at line 87 of file GENIEHelper.h.

{ return fHelperRandom; }

std::string evgb::GENIEHelper::GetTuneName ( ) const

inlineinherited

Definition at line 95 of file GENIEHelper.h.

{ return fTuneName; }

void evgb::GENIEHelper::Initialize ( )

inherited

Definition at line 606 of file GENIEHelper.cxx.

References evgb::GENIEHelper::ConfigGeomScan(), e, evgb::GENIEHelper::fDetectorMass, evgb::GENIEHelper::fDriver, evgb::GENIEHelper::fEventGeneratorList, evgb::GENIEHelper::fEventsPerSpill, evgb::GENIEHelper::fFluxD, evgb::GENIEHelper::fFluxD2GMCJD, evgb::GENIEHelper::fFluxType, evgb::GENIEHelper::fGeomD, evgb::GENIEHelper::fHelperRandom, evgb::GENIEHelper::fHistEventsPerSpill, evgb::GENIEHelper::fPOTPerSpill, evgb::GENIEHelper::fSpillEvents, evgb::GENIEHelper::fSpillExposure, evgb::GENIEHelper::fSurroundingMass, evgb::GENIEHelper::fTotalExposure, evgb::GENIEHelper::fTotalHistFlux, evgb::GENIEHelper::fTuneName, evgb::GENIEHelper::fXSecMassPOT, evgb::GENIEHelper::InitializeFluxDriver(), evgb::GENIEHelper::InitializeGeometry(), evgb::GENIEHelper::ReadXSecTable(), and evgb::SetEventGeneratorListAndTune().

Referenced by evgen::GENIEGen::beginJob(), and evgen::TestGENIEHelper::beginRun().

  {
#ifdef GENIE_PRE_R3
  // trigger early initialization of PDG database & GENIE message service
  // just to get it out of the way and not intermixed with other output
  genie::PDGLibrary::Instance();
#else
  // get the GENIE banner out of the way
  // no longer can use genie::PDGLibrary::Instance() to do this
  // because that must happen, in some cases in v3_02_xx, after the tune
  // is determined
  // banner is triggered by first use of GENIE Messenger
  // avoid using GENIE macros (possible conflict with mf macros)
  // LOG("GENIE",pInfo) << "Trigger GENIE banner";
  (*genie::Messenger::Instance())("GENIE") << log4cpp::Priority::INFO
                                           << "Trigger GENIE banner";
#endif


    mf::LogInfo("GENIEHelper") << "GENIE EventGeneratorList using \""
                               << fEventGeneratorList << "\"";

    evgb::SetEventGeneratorListAndTune(fEventGeneratorList,fTuneName);

#ifdef GENIE_PRE_R3
    // no tunes ... so can report current setting if different from fcl
#else
    genie::RunOpt* grunopt = genie::RunOpt::Instance();

    // RunOpt's SetEventGeneratorList() wasn't introduced until R-3
    // so above call could not have modified it in RunOpt
    std::string currentEvtGenListName = grunopt->EventGeneratorList();
    if ( currentEvtGenListName != fEventGeneratorList ) {
      mf::LogInfo("GENIEHelper") << " EventGeneratorList name changed from \""
                                 << fEventGeneratorList << "\" to \""
                                 << currentEvtGenListName << "\""
                                 << " by evgb::SetEventGeneratorListAndTune";
      fEventGeneratorList = currentEvtGenListName;
    }

    std::string currentTuneName = grunopt->Tune()->Name();
    if ( currentTuneName != fTuneName ) {
      mf::LogInfo("GENIEHelper") << " TuneName name changed from \""
                                 << fTuneName << "\" to \""
                                 << currentTuneName << "\""
                                 << " by evgb::SetEventGeneratorListAndTune";
      fTuneName = currentTuneName;
    }
#endif


    fDriver = new genie::GMCJDriver(); // this needs to be before ConfigGeomScan
    // let the driver know which list
    // (for R-3 this is in _addition_ to setting it in RunOpt)
    // pre-R-3 one could not set it in RunOpt directly (ie. without parsing
    //    command line args), thus don't try to fetch it from there.
    // post-R-3 one *-could-* use
    //   fDriver->SetEventGeneratorList(RunOpt::Instance()->EventGeneratorList());
    fDriver->SetEventGeneratorList(fEventGeneratorList);

    // Figure out which cross section file to use
    // post R-2_8_0 this actually triggers reading the file
    ReadXSecTable();

    // initialize the Geometry and Flux drivers
    InitializeGeometry();
    InitializeFluxDriver();

    fDriver->UseFluxDriver(fFluxD2GMCJD);
    fDriver->UseGeomAnalyzer(fGeomD);

    // must come after creation of Geom, Flux and GMCJDriver
    ConfigGeomScan();      // could trigger fDriver->UseMaxPathLengths(*xmlfile*)

    fDriver->Configure();  // could trigger GeomDriver::ComputeMaxPathLengths()
    fDriver->UseSplines();
    fDriver->ForceSingleProbScale();

    if ( fFluxType.find("histogram") == 0 && fEventsPerSpill < 0.01 ) {
      // fluxes are assumed to be given in units of neutrinos/cm^2/1e20POT/energy
      // integral over all fluxes removes energy dependence
      // histograms should have bin width that reflects the value of the /energy bit
      // ie if /energy = /50MeV then the bin width should be 50 MeV

      // determine product of pot/spill, mass, and cross section
      // events = flux * pot * 10^-38 cm^2 (xsec) * (mass detector (in kg) / nucleon mass (in kg))
      fXSecMassPOT  = 1.e-38*1.e-20;
      fXSecMassPOT *= fPOTPerSpill*(fDetectorMass+fSurroundingMass)/(1.67262158e-27);

      mf::LogInfo("GENIEHelper") << "Number of events per spill will be based on poisson mean of "
                                 << fXSecMassPOT*fTotalHistFlux;

      fHistEventsPerSpill = fHelperRandom->Poisson(fXSecMassPOT*fTotalHistFlux);
    }

    // set the pot/event counters to zero
    fSpillEvents   = 0;
    fSpillExposure = 0.;
    fTotalExposure = 0.;

    // If the flux driver knows how to keep track of exposure (time,pots)
    // reset it now as some might have been used in determining
    // the geometry maxpathlength or internally scanning for weights.
    // Should be happen for all cases since R-2_8_0 .. but no harm doing it ourselves

    fFluxD->Clear("CycleHistory");

    return;
  }

bool evgb::GENIEHelper::Sample ( simb::MCTruth &  truth, simb::MCFlux &  flux, simb::GTruth &  gtruth  )

inherited

Definition at line 1662 of file GENIEHelper.cxx.

References bin, simb::MCParticle::E(), evgb::GENIEHelper::fAddGenieVtxTime, evgb::GENIEHelper::fAtmoSpectralIndex, evgb::GENIEHelper::fDebugFlags, simb::MCFlux::fdk2gen, evgb::GENIEHelper::fDriver, evgb::GENIEHelper::fEventsPerSpill, evgb::GENIEHelper::fFluxD, evgb::GENIEHelper::fFluxD2GMCJD, evgb::GENIEHelper::fFluxHistograms, simb::MCFlux::fFluxType, evgb::GENIEHelper::fFluxType, evgb::GENIEHelper::fForceApplyFlxWgt, simb::MCFlux::fgen2vtx, evgb::GENIEHelper::fGenFlavors, evgb::GENIEHelper::fGenieEventRecord, simb::MCFlux::fgenx, simb::MCFlux::fgeny, simb::MCFlux::fgenz, evgb::GENIEHelper::fGeoManager, evgb::GENIEHelper::fGHepPrintLevel, evgb::GENIEHelper::fGlobalTimeOffset, evgb::GENIEHelper::fHelperRandom, evgb::FillGTruth(), evgb::FillMCFlux(), evgb::FillMCTruth(), evgb::GENIEHelper::fRandomTimeOffset, evgb::GENIEHelper::fSpillEvents, evgb::GENIEHelper::fSpillExposure, evgb::GENIEHelper::fTimeShifter, evgb::GENIEHelper::fTopVolume, evgb::GENIEHelper::fTotalExposure, evgb::GENIEHelper::fTuneName, evgb::GENIEHelper::fUseBlenderDist, evgb::GENIEHelper::fUseHelperRndGen4GENIE, evgb::GENIEHelper::fWorldVolume, simb::MCTruth::GetNeutrino(), simb::kHistPlusFocus, evgb::kNue, evgb::kNueBar, evgb::kNuMu, evgb::kNuMuBar, evgb::kNuTau, evgb::kNuTauBar, simb::MCNeutrino::Nu(), simb::MCFlux::SetFluxGen(), evgb::EvtTimeShiftI::TimeOffset(), and util::to_string().

Referenced by evgen::TestGENIEHelper::produce(), and evgen::GENIEGen::produce().

  {
    // set the top volume for the geometry
    fGeoManager->SetTopVolume(fGeoManager->FindVolumeFast(fTopVolume.c_str()));

    if ( fGenieEventRecord ) delete fGenieEventRecord;

    // ART Framework plays games with gRandom, undo that if requested
    TRandom* old_gRandom = gRandom;
    if (fUseHelperRndGen4GENIE) gRandom = fHelperRandom;

    fGenieEventRecord = fDriver->GenerateEvent();

    if (fForceApplyFlxWgt) {
      //const issue:  double flxweight = fDriver->FluxDriver().Weight();
      //genie::GFluxI& flx = const_cast<genie::GFluxI&>(fDriver->FluxDriver());
      //double flxweight = flx.Weight();
      // use the flux driver handle we already have, rather than fetch const one
      //   from the GMCJDriver fDriver
      double flxweight = fFluxD->Weight();
      double curweight = fGenieEventRecord->Weight();
      fGenieEventRecord->SetWeight(flxweight*curweight);
    }

    if (fUseHelperRndGen4GENIE) gRandom = old_gRandom;

    // now check if we produced a viable event record
    bool viableInteraction = true;
    if ( ! fGenieEventRecord ) viableInteraction = false;

    // update the spill total information, then check to see
    // if we got an event record that was valid

    genie::flux::GFluxExposureI* fexposure =
      dynamic_cast<genie::flux::GFluxExposureI*>(fFluxD);
    if ( fexposure ) {
      fSpillExposure =
        (fexposure->GetTotalExposure()/fDriver->GlobProbScale()) - fTotalExposure;
    }
    // use GENIE2ART code to fill MCFlux
    evgb::FillMCFlux(fFluxD,flux);

    // if no interaction generated return false
    if ( ! viableInteraction ) return false;

    // fill the MCTruth & GTruth information as we have a good interaction
    // these two objects are enough to reconstruct the GENIE record
    // use the new external functions in GENIE2ART

    // choose a time within the spill (ns) to shift the vertex times by:
    double timeoffset = 0;
    if ( ! fTimeShifter ) {
      timeoffset = fHelperRandom->Uniform()*fRandomTimeOffset;
    } else {
      timeoffset = fTimeShifter->TimeOffset();
    }
    // mf::LogInfo("GENIEHelper") << "TimeShifter adding " << timeoffset;
    double spilltime  = fGlobalTimeOffset + timeoffset;

    std::unordered_map<std::string, std::string> genConfig;

    if(fFluxType.find("PowerSpectrum") != std::string::npos){
      genConfig.emplace("SpectralIndex", std::to_string(fAtmoSpectralIndex));
    }

    evgb::FillMCTruth(fGenieEventRecord, spilltime, truth,
                      __GENIE_RELEASE__, fTuneName, fAddGenieVtxTime );
    evgb::FillGTruth(fGenieEventRecord, gtruth);

    // check to see if we are using flux ntuples but want to
    // make n events per spill
    if ( ( fEventsPerSpill > 0 ) &&
         ( fFluxType.find("tree_") == 0 ) ) ++fSpillEvents;

    // now check if using either histogram or mono fluxes, using
    // either n events per spill or basing events on POT per spill for the
    // histogram case
    if ( fFluxType.find("histogram") == 0 ) {
      // set the flag in the parent object that says the
      // fluxes came from histograms and fill related values
      flux.fFluxType = simb::kHistPlusFocus;

      // save the fluxes - fluxes were added to the vector in the same
      // order that the flavors appear in fGenFlavors
      int ctr = 0;
      int bin = fFluxHistograms[0]->FindBin(truth.GetNeutrino().Nu().E());
      std::vector<double> fluxes(6, 0.);
      for(std::vector<int>::iterator i = fGenFlavors.begin(); i != fGenFlavors.end(); i++){
        if(*i ==  12) fluxes[kNue]      = fFluxHistograms[ctr]->GetBinContent(bin);
        if(*i == -12) fluxes[kNueBar]   = fFluxHistograms[ctr]->GetBinContent(bin);
        if(*i ==  14) fluxes[kNuMu]     = fFluxHistograms[ctr]->GetBinContent(bin);
        if(*i == -14) fluxes[kNuMuBar]  = fFluxHistograms[ctr]->GetBinContent(bin);
        if(*i ==  16) fluxes[kNuTau]    = fFluxHistograms[ctr]->GetBinContent(bin);
        if(*i == -16) fluxes[kNuTauBar] = fFluxHistograms[ctr]->GetBinContent(bin);
        ++ctr;
      }

      // get the flux for each neutrino flavor of this energy
      flux.SetFluxGen(fluxes[kNue],   fluxes[kNueBar],
                      fluxes[kNuMu],  fluxes[kNuMuBar],
                      fluxes[kNuTau], fluxes[kNuTauBar]);

      ++fSpillEvents;
    }
    else if ( fFluxType.find("mono")     == 0 ||
              fFluxType.find("function") == 0    ){
      ++fSpillEvents;
    }
    else if ( fFluxType.find("atmo_FLUKA")  == 0 ||
              fFluxType.find("atmo_BARTOL") == 0 ||
              fFluxType.find("atmo_PowerSpectrum") == 0 ||
              fFluxType.find("atmo_BGLRS")  == 0 ||
              fFluxType.find("atmo_HAKKM")  == 0 ||
              fFluxType.find("atmo_HONDA")  == 0    ) {
      if ( fEventsPerSpill > 0 ) ++fSpillEvents;
      flux.fFluxType = simb::kHistPlusFocus;
    }


    // fill these after the Pack[NuMI|Simple]Flux because those
    // will Reset() the values at the start
    TLorentzVector *vertex = fGenieEventRecord->Vertex();
    TLorentzVector nuray_pos = fFluxD->Position();
    TVector3 ray2vtx = nuray_pos.Vect() - vertex->Vect();
    flux.fgenx    = nuray_pos.X();
    flux.fgeny    = nuray_pos.Y();
    flux.fgenz    = nuray_pos.Z();
    flux.fgen2vtx = ray2vtx.Mag();

    genie::flux::GFluxBlender* blender =
      dynamic_cast<genie::flux::GFluxBlender*>(fFluxD2GMCJD);
    if ( blender ) {
      if ( fUseBlenderDist ) flux.fdk2gen = blender->TravelDist();
      // / if mixing flavors print the state of the blender
      if ( fDebugFlags & 0x02 ) blender->PrintState();
    }

    if ( fDebugFlags & 0x04 ) {
      mf::LogInfo("GENIEHelper") << "vertex loc " << vertex->X() << ","
                                 << vertex->Y() << "," << vertex->Z() << std::endl
                                 << " flux ray start " << nuray_pos.X() << ","
                                 << nuray_pos.Y() << "," << nuray_pos.Z() << std::endl
                                 << " ray2vtx = " << flux.fgen2vtx
                                 << " dk2ray = " << flux.fdk2gen;
    }
    if ( fGHepPrintLevel >= 0 ) {
      std::cout << *fGenieEventRecord;
    }

    // set the top volume of the geometry back to the world volume
    fGeoManager->SetTopVolume(fGeoManager->FindVolumeFast(fWorldVolume.c_str()));

    return true;
  }

double evgb::GENIEHelper::SpillExposure ( ) const

inlineinherited

Definition at line 75 of file GENIEHelper.h.

{ return fSpillExposure;  }

bool evgb::GENIEHelper::Stop ( )

inherited

Definition at line 1599 of file GENIEHelper.cxx.

References evgb::GENIEHelper::fAtmoRt, evgb::GENIEHelper::fDriver, evgb::GENIEHelper::fEventsPerSpill, evgb::GENIEHelper::fFluxD, evgb::GENIEHelper::fFluxType, evgb::GENIEHelper::fGenFlavors, evgb::GENIEHelper::fHelperRandom, evgb::GENIEHelper::fHistEventsPerSpill, evgb::GENIEHelper::fPOTPerSpill, evgb::GENIEHelper::fSpillEvents, evgb::GENIEHelper::fSpillExposure, evgb::GENIEHelper::fTotalExposure, evgb::GENIEHelper::fTotalHistFlux, and evgb::GENIEHelper::fXSecMassPOT.

Referenced by evgen::TestGENIEHelper::produce(), and evgen::GENIEGen::produce().

  {
    //   std::cout << "in GENIEHelper::Stop(), fEventsPerSpill = " << fEventsPerSpill
    //      << " fPOTPerSpill = " << fPOTPerSpill << " fSpillExposure = " << fSpillExposure
    //      << " fSpillEvents = " << fSpillEvents
    //      << " fHistEventsPerSpill = " << fHistEventsPerSpill << std::endl;

    // determine if we should keep throwing neutrinos for
    // this spill or move on

    if ( fEventsPerSpill > 0 ) {
      if ( fSpillEvents < fEventsPerSpill ) return false;
    } else {
      // exposure (POT) based
      if ( fFluxType.find("tree_") == 0 ) {
        if ( fSpillExposure < fPOTPerSpill ) {
          return false;
        }
      }
      else if ( fFluxType.find("histogram") == 0 ) {
        if ( fSpillEvents < fHistEventsPerSpill ) return false;
        else fSpillExposure = fPOTPerSpill;
      }
    }

    if ( fFluxType.find("atmo_")  == 0 ) {
      // the exposure for atmo is in SECONDS. In order to get seconds,
      // it needs to be normalized by 1e4 to take into account the units
      // discrepency between AtmoFluxDriver(/m2) and Generate(/cm2)
      // and it need to be normalized by the generation surface area since
      // it's not taken into accoutn in the flux driver

      long int nNeutrinos;

      if(fFluxType.find("PowerSpectrum") != std::string::npos){
        nNeutrinos = dynamic_cast<genie::flux::GPowerSpectrumAtmoFlux *>(fFluxD)->NFluxNeutrinos();
        fTotalExposure = nNeutrinos/fGenFlavors.size()/fDriver->GlobProbScale();
        mf::LogInfo("GENIEHelper")
        << "===> Atmo Pscale*Ngen/Nflavours = " << fTotalExposure;
      }
      else{
        nNeutrinos = dynamic_cast<genie::flux::GAtmoFlux *>(fFluxD)->NFluxNeutrinos();
        fTotalExposure = nNeutrinos * 1.0e4 / (TMath::Pi() * fAtmoRt*fAtmoRt);
        mf::LogInfo("GENIEHelper")
        << "===> Atmo EXPOSURE = " << fTotalExposure << " seconds.";
      }

      

      

    } else {
      fTotalExposure += fSpillExposure;
    }

    // made it to here, means need to reset the counters
    fSpillEvents   = 0;
    fSpillExposure = 0.;
    fHistEventsPerSpill = fHelperRandom->Poisson(fXSecMassPOT*fTotalHistFlux);
    return true;
  }

double evgb::GENIEHelper::TotalExposure ( ) const

inlineinherited

Definition at line 71 of file GENIEHelper.h.

Referenced by evgen::GENIEGen::beginSubRun(), evgen::TestGENIEHelper::endSubRun(), and evgen::GENIEGen::endSubRun().

{ return fTotalExposure;  }

double evgb::GENIEHelper::TotalHistFlux ( )

inherited

Definition at line 591 of file GENIEHelper.cxx.

References evgb::GENIEHelper::fFluxType, and evgb::GENIEHelper::fTotalHistFlux.

  {
    if ( fFluxType.find("mono")     == 0 ||
         fFluxType.find("function") == 0 ||
         fFluxType.find("tree_")    == 0 ||
         fFluxType.find("atmo_")    == 0    ) {
      // shouldn't be asking for this for any of the above
      // perhaps not-not "function"?
      return -999.;
    }

    return fTotalHistFlux;
  }

double evgb::GENIEHelper::TotalMass ( ) const

inlineinherited

Definition at line 82 of file GENIEHelper.h.

{ return fDetectorMass+fSurroundingMass; }

Member Data Documentation

genie::EventRecord* evgb::GiBUUHelper::fGiBUUEventRecord

private

Definition at line 34 of file GiBUUHelper.h.

Referenced by GetGiBUUEventRecord(), and ~GiBUUHelper().

---

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

• nugen/EventGeneratorBase/GiBUU/GiBUUHelper.h
• nugen/EventGeneratorBase/GiBUU/GiBUUHelper.cxx