Add material input; This build ~works~ as long as you only run one particle type at time

Author: Jack Smedley

sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx

@@ -59,14 +59,15 @@ class Geant4WeightCalc : public WeightCalc {
   std::string fMCParticleProducer;  //!< Label for the MCParticle producer
   std::string fMCTruthProducer;  //!< Label for the MCTruth producer
   CLHEP::RandGaussQ* fGaussRandom;  //!< Random number generator
+  fhicl::ParameterSet fMaterial; //!< Detector material, i.e. LAr
   // std::map<int, ParticleDef> fParticles;  //!< Particles to reweight
   unsigned fNsims;  //!< Number of multisims
   int fPdg; //!< PDG value for particles that a given weight calculator should apply to. Note that for now this module can only handle weights for one particle species at a time.
   // float fXSUncertainty;  //!< Flat cross section uncertainty
-  G4ReweightManager *RWManager; //!< "holds the run manager and creates the detector" according to the commit message (?)
-  G4ReweighterFactory RWFactory; //!< Base class to handle all Geant4Reweighters
-  G4Reweighter *theReweighter; //!< Geant4Reweighter -- this is what provides the weights
-  G4ReweightParameterMaker *ParMaker;
+  G4ReweightManager *fRWManager; //!< "holds the run manager and creates the detector" according to the commit message (?)
+  G4ReweighterFactory fRWFactory; //!< Base class to handle all Geant4Reweighters
+  G4Reweighter *fReweighter; //!< Geant4Reweighter -- this is what provides the weights
+  G4ReweightParameterMaker *fParMaker;
   std::vector<std::map<std::string, double>> UniverseVals; //!< Vector of maps relating parameter name to value (defines parameter values that will be evaluated in universes). Each map should have one entry per parameter we are considering
 
   art::ServiceHandle < geo::Geometry > fGeometryService;
@@ -117,6 +118,7 @@ void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p,
   std::vector< fhicl::ParameterSet > FitParSets = pset.get< std::vector< fhicl::ParameterSet > >("parameters");
   fNsims = pset.get<int> ("number_of_multisims", 0);
   fPdg = pset.get<int> ("pdg_to_reweight");
+  fMaterial = pset.get<fhicl::ParameterSet> ("material");
   fDebug = pset.get<bool> ("debug",false);
 
   // Prepare random generator
@@ -127,9 +129,9 @@ void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p,
   TFile XSecFile( XSecFileName.c_str(), "OPEN" );
 
   // Configure G4Reweighter
-  ParMaker = new G4ReweightParameterMaker( FitParSets, true , fPdg ); //TODO:Do we want check_overlap? Maybe a fcl switch?
-  RWManager = new G4ReweightManager( FitParSets );
-  theReweighter = RWFactory.BuildReweighter(fPdg, &FracsFile, ParMaker->GetFSHists(), pset, RWManager, ParMaker->GetElasticHist() );
+  fParMaker = new G4ReweightParameterMaker( FitParSets, true , fPdg ); //TODO:Do we want check_overlap? Maybe a fcl switch?
+  fRWManager = new G4ReweightManager( {fMaterial} ); //Constructor asks for a vector, but SBN only cares about one material
+  fReweighter = fRWFactory.BuildReweighter(fPdg, &FracsFile, fParMaker->GetFSHists(), fMaterial, fRWManager, fParMaker->GetElasticHist() );
 
   // Make output trees to save things for quick and easy validation
   art::ServiceHandle<art::TFileService> tfs;
@@ -441,11 +443,11 @@ std::vector<float> Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) {
         float w/*, el_w*/;
 
         // I think this is the only bit that needs to change for different universes -- all the above is jut about the track, which doesn't change based on universe
-        ParMaker->SetNewVals(UniverseVals.at(j));
-        theReweighter->SetNewHists(ParMaker->GetFSHists());
-        theReweighter->SetNewElasticHists(ParMaker->GetElasticHist());
+        fParMaker->SetNewVals(UniverseVals.at(j));
+        fReweighter->SetNewHists(fParMaker->GetFSHists());
+        fReweighter->SetNewElasticHists(fParMaker->GetElasticHist());
         //Get the weight from the G4ReweightTraj
-        w = theReweighter->GetWeight( &theTraj );
+        w = fReweighter->GetWeight( &theTraj );
         // Total weight is the product of track weights in the event
         weight[j] *= std::max((float)0.0, w);
 
@@ -454,7 +456,7 @@ std::vector<float> Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) {
 
 /*
         // Do the same for elastic weight (should be 1 unless set to non-nominal )
-        el_w = theReweighter->GetElasticWeight( &theTraj );
+        el_w = fReweighter->GetElasticWeight( &theTraj );
         weight[j] *= std::max((float)0.0,el_w);
 
         // just for the output tree

sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl

@@ -14,6 +14,7 @@
 #include "piplus_reweight_parameters.fcl"
 #include "piminus_reweight_parameters.fcl"
 #include "proton_reweight_parameters.fcl"
+#include "material_parameters.fcl"
 
 BEGIN_PROLOG
 
@@ -23,48 +24,53 @@ sbn_eventweight_geant4: {
   min_weight: 0.0
   max_weight: 100
 
-  weight_functions_reint: [ reinteractions_piplus, reinteractions_piminus, reinteractions_proton ]
+  weight_functions_reint: [ reinteractions_piplus] #, reinteractions_piminus, reinteractions_proton ]
 
   reinteractions_piplus: {
     type: Geant4
+    material: @local::LAr
     random_seed: 58
     parameters: @local::PiPlusParameters
     mode: multisim
     number_of_multisims: 1000
     fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_piplus.root"
     xsecfile:  "$SBNDATA_DIR/systematics/reint/g4_cross_section_piplus.root"
-    makeoutputtree: false
+    makeoutputtree: true
     pdg_to_reweight: 211
-    debug: false
+    debug: true
   }
 
   reinteractions_piminus: {
     type: Geant4
+    material: @local::LAr
     random_seed: 59
     parameters: @local::PiMinusParameters
     mode: multisim
     number_of_multisims: 1000
     fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_piminus.root"
     xsecfile:  "$SBNDATA_DIR/systematics/reint/g4_cross_section_piminus.root"
-    makeoutputtree: false
+    makeoutputtree: true
     pdg_to_reweight: -211
-    debug: false
+    debug: true
   }
 
   reinteractions_proton: {
     type: Geant4
+    material: @local::LAr
     random_seed: 60
     parameters: @local::ProtonParameters
     mode: multisim
     number_of_multisims: 1000
     fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_proton.root"
     xsecfile:  "$SBNDATA_DIR/systematics/reint/g4_cross_section_proton.root"
-    makeoutputtree: false
+    makeoutputtree: true
     pdg_to_reweight: 2212
-    debug: false
+    debug: true
   }
 }
 
+# reinteractions_neutron
+
 ###
 # Reweighting parameters should be defined as
 #

sbncode/SBNEventWeight/jobs/geant4/material_parameters.fcl

@@ -0,0 +1,18 @@
+BEGIN_PROLOG
+
+# LAr parameters for geant4reweight, stolen out of the ND-Lar 2x2 implementation
+# A more complete solution would be to get this info out of the geometry instead...
+
+LAr: {
+  Name: "liquidArgon"
+  Density: 1.390
+  Components: [
+    {
+      Z:    18
+      Mass: 39.95
+      Fraction: 1.
+    }
+  ]
+}
+
+END_PROLOG