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)

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 50 of file GiBUUHelper.cxx.

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

evgb::GiBUUHelper::~GiBUUHelper

(

)

Definition at line 60 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 72 of file GENIEHelper.h.

{ return fDetLocation;    }

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

inlineinherited

Definition at line 76 of file GENIEHelper.h.

{ return fFluxHistograms; }

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

inlineinherited

Definition at line 71 of file GENIEHelper.h.

{ return fFluxType;       }

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

inlineinherited

Definition at line 86 of file GENIEHelper.h.

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

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

inlineinherited

Definition at line 79 of file GENIEHelper.h.

{ 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 82 of file GENIEHelper.h.

{ return fHelperRandom; }

void evgb::GENIEHelper::Initialize ( )

inherited

Definition at line 601 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::fXSecMassPOT, evgb::GENIEHelper::InitializeFluxDriver(), and evgb::GENIEHelper::InitializeGeometry().

Referenced by evgen::GENIEGen::beginJob().

  {
    fDriver = new genie::GMCJDriver(); // needs to be before ConfigGeomScan
#ifndef GENIE_USE_ENVVAR
    // this configuration happened via $GEVGL beore R-2_8_0
    fDriver->SetEventGeneratorList(fEventGeneratorList);
#endif

    // 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();  // trigger GeomDriver::ComputeMaxPathLengths()
    fDriver->UseSplines();
    fDriver->ForceSingleProbScale();

    if ( fFluxType.compare("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.
    // This should not be necessary (GENIE should do it automagically)
    // but the current version (as of 3665) doesn't.

    // alas, at this time there is no "unified" interface for Clear()
    // as there is for GetTotalExposure (replacing UsedPOTs())
    double preUsedFluxPOTs = 0;
    bool   wasCleared = true;
    bool   doprintpre = false;
    if( fFluxType.compare("numi") == 0 ) {
      genie::flux::GNuMIFlux* gnumiflux =
        dynamic_cast<genie::flux::GNuMIFlux *>(fFluxD);
      preUsedFluxPOTs = gnumiflux->UsedPOTs();
      if ( preUsedFluxPOTs > 0 ) {
        doprintpre = true;
        gnumiflux->Clear("CycleHistory");
        if ( gnumiflux->UsedPOTs() != 0 ) wasCleared = false;
      }
    } else if ( fFluxType.compare("simple") == 0 ) {
      genie::flux::GSimpleNtpFlux* gsimpleflux =
        dynamic_cast<genie::flux::GSimpleNtpFlux *>(fFluxD);
      preUsedFluxPOTs = gsimpleflux->UsedPOTs();
      if ( preUsedFluxPOTs > 0 ) {
        doprintpre = true;
        gsimpleflux->Clear("CycleHistory");
        if ( gsimpleflux->UsedPOTs() != 0 ) wasCleared = false;
      }
#if __GENIE_RELEASE_CODE__ >= GRELCODE(2,11,0)
    } else if ( fFluxType.compare("dk2nu") == 0 ) {
      genie::flux::GDk2NuFlux* dk2nuflux =
        dynamic_cast<genie::flux::GDk2NuFlux *>(fFluxD);
      preUsedFluxPOTs = dk2nuflux->UsedPOTs();
      if ( preUsedFluxPOTs > 0 ) {
        doprintpre = true;
        dk2nuflux->Clear("CycleHistory");
        if ( dk2nuflux->UsedPOTs() != 0 ) wasCleared = false;
      }
#endif
    }
    if ( doprintpre ) {
      double probscale = fDriver->GlobProbScale();
      mf::LogInfo("GENIEHelper")
        << "Pre-Event Generation: "
        << " FluxDriver base " << preUsedFluxPOTs
        << " / GMCJDriver GlobProbScale " << probscale
        << " = used POTS " << preUsedFluxPOTs/TMath::Max(probscale,1.0e-100)
        << " "
        << (wasCleared?"successfully":"failed to") << " cleared count for "
        << fFluxType;
    }
    return;
  }

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

inherited

Definition at line 1486 of file GENIEHelper.cxx.

References bin, simb::MCParticle::E(), 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, 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::fUseHelperRndGen4GENIE, evgb::GENIEHelper::fWorldVolume, simb::MCTruth::GetNeutrino(), simb::kHistPlusFocus, simb::kNtuple, evgb::kNue, evgb::kNueBar, evgb::kNuMu, evgb::kNuMuBar, evgb::kNuTau, evgb::kNuTauBar, simb::kSimple_Flux, simb::MCNeutrino::Nu(), evgb::GENIEHelper::PackGTruth(), evgb::GENIEHelper::PackMCTruth(), evgb::GENIEHelper::PackNuMIFlux(), evgb::GENIEHelper::PackSimpleFlux(), simb::MCFlux::SetFluxGen(), and evgb::EvtTimeShiftI::TimeOffset().

Referenced by 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 (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


#if __GENIE_RELEASE_CODE__ >= GRELCODE(2,11,0)
    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);
#else
    // pack the flux information no support for dk2nu
    if ( fFluxType.compare("numi") == 0){
      fSpillExposure = (dynamic_cast<genie::flux::GNuMIFlux *>(fFluxD)->UsedPOTs()/fDriver->GlobProbScale() - fTotalExposure);
      flux.fFluxType = simb::kNtuple;
      PackNuMIFlux(flux);
    }
    else if ( fFluxType.compare("simple")==0 ) {
      // pack the flux information
      fSpillExposure = (dynamic_cast<genie::flux::GSimpleNtpFlux *>(fFluxD)->UsedPOTs()/fDriver->GlobProbScale() - fTotalExposure);
      flux.fFluxType = simb::kSimple_Flux;
      PackSimpleFlux(flux);
    }
#endif

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

#if __GENIE_RELEASE_CODE__ >= GRELCODE(2,11,0)
    // 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 spill time (ns) to shift the vertex times by:
    double spilltime = fGlobalTimeOffset;
    if ( ! fTimeShifter ) {
      spilltime += fHelperRandom->Uniform()*fRandomTimeOffset;
    } else {
      spilltime += fTimeShifter->TimeOffset();
    }

    evgb::FillMCTruth(fGenieEventRecord, spilltime, truth);
    evgb::FillGTruth(fGenieEventRecord, gtruth);
#else
    // fill the MC truth information as we have a good interaction
    PackMCTruth(fGenieEventRecord,truth);
    // fill the Generator (genie) truth information
    PackGTruth(fGenieEventRecord, gtruth);
#endif

    // check to see if we are using flux ntuples but want to
    // make n events per spill
    if ( fEventsPerSpill > 0 &&
         ( fFluxType.compare("numi")    == 0 ||
           fFluxType.compare("simple")  == 0 ||
           fFluxType.compare("dk2nu")   == 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.compare("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.compare("mono") == 0 || fFluxType.compare("function") == 0){
      ++fSpillEvents;
    }

    else if(fFluxType.compare("atmo_FLUKA") == 0 || fFluxType.compare("atmo_BARTOL") == 0 ||
            fFluxType.compare("atmo_HAKKM") == 0 || fFluxType.compare("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 ) {
      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 70 of file GENIEHelper.h.

{ return fSpillExposure;  }

bool evgb::GENIEHelper::Stop ( )

inherited

Definition at line 1431 of file GENIEHelper.cxx.

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

Referenced by 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 ( fFluxType.compare("atmo_FLUKA") == 0 || fFluxType.compare("atmo_BARTOL") == 0 ||
         fFluxType.compare("atmo_HONDA") == 0 || fFluxType.compare("atmo_HAKKM")  == 0    ) {
      if ( (fEventsPerSpill > 0) && (fSpillEvents < fEventsPerSpill) ) {
        return false;
      }
    }

    else if ( fEventsPerSpill > 0 ) {
      if ( fSpillEvents < fEventsPerSpill )
        return false;
    } else {
      if ( ( fFluxType.compare("numi")   == 0 ||
             fFluxType.compare("simple") == 0 ||
             fFluxType.compare("dk2nu")  == 0    ) &&
           fSpillExposure < fPOTPerSpill ) return false;
      else if ( fFluxType.compare("histogram") == 0 ){
        if ( fSpillEvents < fHistEventsPerSpill ) return false;
        else fSpillExposure = fPOTPerSpill;
      }
    }

    // made it to here, means need to reset the counters

    if(fFluxType.compare("atmo_FLUKA") == 0 || fFluxType.compare("atmo_BARTOL") == 0 ||
       fFluxType.compare("atmo_HONDA") == 0 || fFluxType.compare("atmo_HAKKM") == 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
      fTotalExposure = (1e4 * (dynamic_cast<genie::flux::GAtmoFlux *>(fFluxD)->NFluxNeutrinos())) / (TMath::Pi() * fAtmoRt*fAtmoRt);

      LOG_DEBUG("GENIEHelper") << "===> Atmo EXPOSURE = " << fTotalExposure << " seconds";
    }

    else{
      fTotalExposure += fSpillExposure;
    }

    fSpillEvents   = 0;
    fSpillExposure = 0.;
    fHistEventsPerSpill = fHelperRandom->Poisson(fXSecMassPOT*fTotalHistFlux);
    return true;
  }

double evgb::GENIEHelper::TotalExposure ( ) const

inlineinherited

Definition at line 66 of file GENIEHelper.h.

Referenced by evgen::GENIEGen::endSubRun().

{ return fTotalExposure;  }

double evgb::GENIEHelper::TotalHistFlux ( )

inherited

Definition at line 586 of file GENIEHelper.cxx.

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

  {
    if ( fFluxType.compare("mono")     == 0 ||
         fFluxType.compare("function") == 0 ||
         fFluxType.compare("numi")     == 0 ||
         fFluxType.compare("simple")   == 0 ||
         fFluxType.compare("dk2nu")    == 0    ) {
      // shouldn't be asking for this for any of the above
      return -999.;
    }

    return fTotalHistFlux;
  }

double evgb::GENIEHelper::TotalMass ( ) const

inlineinherited

Definition at line 77 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:

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