From 19032c9b5354d295dceaa1ad284e20c89b7996ca Mon Sep 17 00:00:00 2001 From: Steven Gardiner Date: Tue, 25 Jul 2023 15:54:01 -0500 Subject: [PATCH 01/16] cherry-pick d7281db, fix merge conflict --- CMakeLists.txt | 1 + .../SBNEventWeight/Calculators/CMakeLists.txt | 2 +- .../Geant4/BetheBlochForG4ReweightValid.h | 145 +++++ .../Calculators/Geant4/CMakeLists.txt | 22 + .../Calculators/Geant4/Geant4WeightCalc.cxx | 524 ++++++++++++++++++ ups/product_deps | 11 +- 6 files changed, 699 insertions(+), 6 deletions(-) create mode 100644 sbncode/SBNEventWeight/Calculators/Geant4/BetheBlochForG4ReweightValid.h create mode 100644 sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt create mode 100644 sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx diff --git a/CMakeLists.txt b/CMakeLists.txt index fc36b1322..d5b3927ad 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -77,6 +77,7 @@ find_package( CLHEP REQUIRED ) find_package( ROOT REQUIRED ) find_package( Geant4 REQUIRED ) find_package( Boost COMPONENTS system filesystem REQUIRED ) +find_package( geant4reweight REQUIRED ) # macros for dictionary and simple_plugin include(ArtDictionary) diff --git a/sbncode/SBNEventWeight/Calculators/CMakeLists.txt b/sbncode/SBNEventWeight/Calculators/CMakeLists.txt index e2080685f..383d41a3a 100644 --- a/sbncode/SBNEventWeight/Calculators/CMakeLists.txt +++ b/sbncode/SBNEventWeight/Calculators/CMakeLists.txt @@ -1,4 +1,4 @@ add_subdirectory(CrossSections) add_subdirectory(BNBFlux) -#add_subdirectory(Geant4) +add_subdirectory(Geant4) diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/BetheBlochForG4ReweightValid.h b/sbncode/SBNEventWeight/Calculators/Geant4/BetheBlochForG4ReweightValid.h new file mode 100644 index 000000000..58cfdc99e --- /dev/null +++ b/sbncode/SBNEventWeight/Calculators/Geant4/BetheBlochForG4ReweightValid.h @@ -0,0 +1,145 @@ + + +double BetheBloch(double energy, double mass){ + + //Need to make this configurable? Or delete... + double K = .307075; + double rho = 1.390; + double Z = 18; + double A = 40; + double I = 188E-6; + double me = .511; + //Need to make sure this is total energy, not KE + double gamma = energy/mass; + double beta = sqrt( 1. - (1. / (gamma*gamma)) ); + double Tmax = 2 * me * beta*beta * gamma*gamma; + + double first = K * (Z/A) * rho / (beta*beta); + double second = .5 * log(Tmax*Tmax/(I*I)) - beta*beta; + + double dEdX = first*second; + return dEdX; +} + +std::vector< std::pair > ThinSliceBetheBloch(G4ReweightTraj * theTraj, double res, double mass, bool isElastic){ + + std::vector< std::pair > result; + + //First slice position +// double sliceEdge = res; +// double lastPos = 0.; +// double nextPos = 0.; +// double px,py,pz; + int interactInSlice = 0; + + //Get total distance traveled in z +// double totalDeltaZ = 0.; + double disp = 0.; +// double oldDisp = 0.; +// int crossedSlices = 0; + + int currentSlice = 0; + int oldSlice = 0; + + double sliceEnergy = theTraj->GetEnergy(); + size_t nSteps = theTraj->GetNSteps(); + for(size_t is = 0; is < nSteps; ++is){ + auto theStep = theTraj->GetStep(is); + + disp += theStep->GetStepLength(); + currentSlice = floor(disp/res); + + std::string theProc = theStep->GetStepChosenProc(); + + //Check to see if in a new slice and it's not the end + if( oldSlice != currentSlice && is < nSteps - 1){ + + + //Save Interaction info of the prev slice + //and reset + result.push_back( std::make_pair(sliceEnergy, interactInSlice) ); + interactInSlice = 0; + + //Update the energy + sliceEnergy = sliceEnergy - res*BetheBloch(sliceEnergy, mass); + if( sliceEnergy - mass < 0.){ + //std::cout << "Warning! Negative energy " << sliceEnergy - mass << std::endl; + //std::cout << "Crossed " << oldSlice - currentSlice << std::endl; + sliceEnergy = 0.0001; + } + //If it's more than 1 slice, add in non-interacting slices + for(int ic = 1; ic < abs( oldSlice - currentSlice ); ++ic){ + //std::cout << ic << std::endl; + + result.push_back( std::make_pair(sliceEnergy, 0) ); + + //Update the energy again + sliceEnergy = sliceEnergy - res*BetheBloch(sliceEnergy, mass); + if( sliceEnergy - mass < 0.){ + //std::cout << "Warning! Negative energy " << sliceEnergy - mass << std::endl; + //std::cout << "Crossed " << oldSlice - currentSlice << std::endl; + sliceEnergy = 0.0001; + } + } + + if ((!isElastic && theProc.find(std::string("Inelastic")) != std::string::npos) || (isElastic && theProc.find(std::string("hadElastic")) != std::string::npos)) { + // std::cout << "found! " << theProc << '\n'; + interactInSlice = 1; + } + } + //It's crossed a slice and it's the last step. Save both info + else if( oldSlice != currentSlice && is == nSteps - 1 ){ + result.push_back( std::make_pair(sliceEnergy, interactInSlice) ); + interactInSlice = 0; + + //Update the energy + sliceEnergy = sliceEnergy - res*BetheBloch(sliceEnergy, mass); + if( sliceEnergy - mass < 0.){ + //std::cout << "Warning! Negative energy " << sliceEnergy - mass << std::endl; + //std::cout << "Crossed " << oldSlice - currentSlice << std::endl; + sliceEnergy = 0.0001; + } + //If it's more than 1 slice, add in non-interacting slices + for(int ic = 1; ic < abs( oldSlice - currentSlice ); ++ic){ + //std::cout << ic << std::endl; + + result.push_back( std::make_pair(sliceEnergy, 0) ); + + //Update the energy again + sliceEnergy = sliceEnergy - res*BetheBloch(sliceEnergy, mass); + if( sliceEnergy - mass < 0.){ + //std::cout << "Warning! Negative energy " << sliceEnergy - mass << std::endl; + //std::cout << "Crossed " << oldSlice - currentSlice << std::endl; + sliceEnergy = 0.0001; + } + } + + //Save the last slice + if ((!isElastic && theProc.find(std::string("Inelastic")) != std::string::npos) || (isElastic && theProc.find(std::string("hadElastic")) != std::string::npos)) { + // std::cout << "found! " << theProc << '\n'; + interactInSlice = 1; + } + result.push_back( std::make_pair(sliceEnergy, interactInSlice) ); + } + //It's the end, so just save this last info + else if( oldSlice == currentSlice && is == nSteps - 1 ){ + if ((!isElastic && theProc.find(std::string("Inelastic")) != std::string::npos) || (isElastic && theProc.find(std::string("hadElastic")) != std::string::npos)) { + // std::cout << "found! " << theProc << '\n'; + interactInSlice = 1; + } + result.push_back( std::make_pair(sliceEnergy, interactInSlice) ); + } + //Same slice, not the end. Check for interactions + else{ + if ((!isElastic && theProc.find(std::string("Inelastic")) != std::string::npos) || (isElastic && theProc.find(std::string("hadElastic")) != std::string::npos)) { + // std::cout << "found! " << theProc << '\n'; + interactInSlice = 1; + } + } + + //Update oldslice + oldSlice = currentSlice; + } + + return result; +} diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt b/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt new file mode 100644 index 000000000..e1289208d --- /dev/null +++ b/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt @@ -0,0 +1,22 @@ +include_directories ( $ENV{GEANT4_FQ_DIR}/include ) +include_directories ( $ENV{GEANT4REWEIGHT_INC} ) + +art_make_library( + LIBRARY_NAME sbncode_SBNEventWeight_Calculators_Geant4 + LIBRARIES + sbncode_SBNEventWeight_Base + nusimdata::SimulationBase + nurandom::RandomUtils_NuRandomService_service + art::Framework_Principal + art::Framework_Services_Registry + art_root_io::TFileService_service + larcore::Geometry_Geometry_service + cetlib_except::cetlib_except + ReweightBaseLib + PropBaseLib +) + +install_headers() +install_fhicl() +install_source() + diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx new file mode 100644 index 000000000..f4ba26774 --- /dev/null +++ b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx @@ -0,0 +1,524 @@ +/** + * \class evwgh::Geant4WeightCalc + * \brief Updated hadron reinteraction event reweighting, using geant4reweight + * \author K. Duffy , 2019/10 + * + * Reweight events based on hadron reinteraction probabilities. + */ + +#include +#include +#include "TDirectory.h" +#include "TFile.h" +#include "TH1D.h" +#include "TTree.h" +#include "Geant4/G4LossTableManager.hh" +#include "Geant4/G4ParticleTable.hh" +#include "Geant4/G4ParticleDefinition.hh" +#include "Geant4/G4Material.hh" +#include "Geant4/G4MaterialCutsCouple.hh" +#include "art_root_io/TFileService.h" +#include "art_root_io/TFileDirectory.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Persistency/Provenance/ModuleContext.h" +#include "art/Framework/Principal/Event.h" +#include "canvas/Persistency/Common/FindManyP.h" +#include "canvas/Persistency/Common/Ptr.h" +#include "nusimdata/SimulationBase/MCParticle.h" +#include "nusimdata/SimulationBase/MCTruth.h" +#include "CLHEP/Random/RandGaussQ.h" +#include "CLHEP/Units/PhysicalConstants.h" +#include "CLHEP/Units/SystemOfUnits.h" +#include "sbncode/SBNEventWeight/Base/WeightCalcCreator.h" +#include "sbncode/SBNEventWeight/Base/WeightCalc.h" +#include "larcore/Geometry/Geometry.h" +#include "geant4reweight/ReweightBase/G4ReweighterFactory.hh" +#include "geant4reweight/ReweightBase/G4Reweighter.hh" +#include "geant4reweight/ReweightBase/G4ReweightTraj.hh" +#include "geant4reweight/ReweightBase/G4ReweightStep.hh" +#include "geant4reweight/PropBase/G4ReweightParameterMaker.hh" + +// local include +#include "BetheBlochForG4ReweightValid.h" + +namespace sbn { + +namespace evwgh { + +class Geant4WeightCalc : public WeightCalc { +public: + Geant4WeightCalc() {} + + void Configure(fhicl::ParameterSet const& p, CLHEP::HepRandomEngine& engine); + + std::vector > GetWeight(art::Event& e); + +private: + std::string fMCParticleProducer; //!< Label for the MCParticle producer + std::string fMCTruthProducer; //!< Label for the MCTruth producer + CLHEP::RandGaussQ* fGaussRandom; //!< Random number generator + // std::map 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 + G4ReweighterFactory RWFactory; //!< Base class to handle all Geant4Reweighters (right now "all" means pi+, pi-, p) + G4Reweighter *theReweighter; //!< Geant4Reweighter -- this is what provides the weights + G4ReweightParameterMaker *ParMaker; + std::vector> 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; + + bool fMakeOutputTrees; ///!< Fcl parameter to decide whether to save output tree (useful for validations but not necessary when running in production) + TTree *fOutTree_MCTruth; //!< Output tree for quick validations: on entry per MCTruth object + TTree *fOutTree_Particle; //!< Output tree for quick validations: one entry per neutrino-induced pi+, pi-, or proton + int event_num; //!< Variables for both output trees + int run_num; //!< Variables for both output trees + int subrun_num; //!< Variables for both output trees + double p_track_length; //!< Variables for by-particle output tree + int p_PDG; //!< Variables for by-particle output tree + std::string p_final_proc; //!< Variables for by-particle output tree + double p_init_momentum; //!< Variables for by-particle output tree + double p_final_momentum; //!< Variables for by-particle output tree + std::vector< double > p_energies_el; //!< Variables for by-particle output tree + std::vector< int > p_sliceInts_el; //!< Variables for by-particle output tree + std::vector< double > p_energies_inel; //!< Variables for by-particle output tree + std::vector< int > p_sliceInts_inel; //!< Variables for by-particle output tree + int p_nElasticScatters; //!< Variables for by-particle output tree + std::vector p_inel_weight; //!< Variables for by-particle output tree + std::vector p_elastic_weight; //!< Variables for by-particle output + std::vector > e_inel_weight; //!< Variables for by-event output tree + std::vector > e_elastic_weight; //!< Variables for by-event output tree + + bool fDebug; //!< print debug statements + + DECLARE_WEIGHTCALC(Geant4WeightCalc) +}; + + +void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, + CLHEP::HepRandomEngine& engine) +{ + std::cout << "Using Geant4WeightCalc for reinteraction weights" << std::endl; + + // Get configuration for this function + fhicl::ParameterSet const& pset = p.get(GetName()); + // std::cout << pset.GetName() << std::endl; + fMCParticleProducer = pset.get("MCParticleProducer", "largeant"); + fMCTruthProducer = pset.get("MCTruthProducer", "generator"); + fMakeOutputTrees = pset.get< bool >( "makeoutputtree", false ); + std::string mode = pset.get("mode"); + std::string FracsFileName = pset.get< std::string >( "fracsfile" ); + std::string XSecFileName = pset.get< std::string >( "xsecfile" ); + std::vector< fhicl::ParameterSet > FitParSets = pset.get< std::vector< fhicl::ParameterSet > >("parameters"); + fNsims = pset.get ("number_of_multisims", 0); + fPdg = pset.get ("pdg_to_reweight"); + fDebug = pset.get ("debug",false); + + // Prepare random generator + fGaussRandom = new CLHEP::RandGaussQ(engine); + + // Get input files + TFile FracsFile( FracsFileName.c_str(), "OPEN" ); + TFile XSecFile( XSecFileName.c_str(), "OPEN" ); + + // Configure G4Reweighter + bool totalonly = false; + if (fPdg==2212) totalonly = true; + ParMaker = new G4ReweightParameterMaker( FitParSets, totalonly ); + theReweighter = RWFactory.BuildReweighter(fPdg, &XSecFile, &FracsFile, ParMaker->GetFSHists(), ParMaker->GetElasticHist() ); + + // Make output trees to save things for quick and easy validation + art::ServiceHandle tfs; + + if (fMakeOutputTrees){ + fOutTree_Particle = tfs->make(Form("%s_%i","ByParticleValidTree",fPdg),""); + fOutTree_Particle->Branch("event_num",&event_num); + fOutTree_Particle->Branch("run_num",&run_num); + fOutTree_Particle->Branch("subrun_num",&subrun_num); + fOutTree_Particle->Branch("UniverseVals", &UniverseVals); + fOutTree_Particle->Branch("pdg_to_reweight",&fPdg); + fOutTree_Particle->Branch("inelastic_weight",&p_inel_weight); + fOutTree_Particle->Branch("elastic_weight",&p_elastic_weight); + fOutTree_Particle->Branch("track_length",&p_track_length); + fOutTree_Particle->Branch("PDG",&p_PDG); + fOutTree_Particle->Branch("final_proc",&p_final_proc); + fOutTree_Particle->Branch("init_momentum",&p_init_momentum); + fOutTree_Particle->Branch("final_momentum",&p_final_momentum); + fOutTree_Particle->Branch("energies_el",&p_energies_el); + fOutTree_Particle->Branch("sliceInts_el",&p_sliceInts_el); + fOutTree_Particle->Branch("energies_inel",&p_energies_inel); + fOutTree_Particle->Branch("sliceInts_inel",&p_sliceInts_inel); + fOutTree_Particle->Branch("nElasticScatters",&p_nElasticScatters); + + fOutTree_MCTruth = tfs->make(Form("%s_%i","ByMCTruthValidTree",fPdg),""); + fOutTree_MCTruth->Branch("event_num",&event_num); + fOutTree_MCTruth->Branch("run_num",&run_num); + fOutTree_MCTruth->Branch("subrun_num",&subrun_num); + fOutTree_MCTruth->Branch("UniverseVals", &UniverseVals); + fOutTree_MCTruth->Branch("pdg_to_reweight",&fPdg); + fOutTree_MCTruth->Branch("inelastic_weight",&e_inel_weight); + fOutTree_MCTruth->Branch("elastic_weight",&e_elastic_weight); + + } + + + // Read input parameter sets and set up universes + size_t n_parsets = FitParSets.size(); + std::vector FitParNames; + std::vector FitParNominals; + std::vector FitParSigmas; + std::map theNominals; + + for (size_t i_parset=0; i_parset("Name"); + double theNominal = theSet.get("Nominal",1.); + double theSigma = theSet.get("Sigma",0.); + + FitParNames.push_back(theName); + FitParNominals.push_back(theNominal); + FitParSigmas.push_back(theSigma); + + theNominals[theName] = theNominal; + } + + if (mode=="pm1sigma"){ + // pm1sigma mode: 0 = +1sigma, 1 = -1sigma of a single parameter. All other parameters at nominal + for (size_t i_parset=0; i_parset tmp_vals_p1sigma(theNominals); + std::map tmp_vals_m1sigma(theNominals); + // Now reset the +1sigma and -1sigma values for this parameter set only + tmp_vals_p1sigma[FitParNames.at(i_parset)] = FitParNominals.at(i_parset)+FitParSigmas.at(i_parset); + tmp_vals_m1sigma[FitParNames.at(i_parset)] = FitParNominals.at(i_parset)-FitParSigmas.at(i_parset); + + if (fDebug){ + std::cout << "Universe " << i_parset*2 << ": " << FitParNames.at(i_parset) << " = " << FitParNominals.at(i_parset)+FitParSigmas.at(i_parset) << std::endl; + std::cout << "Universe " << i_parset*2+1 << ": " << FitParNames.at(i_parset) << " = " << FitParNominals.at(i_parset)-FitParSigmas.at(i_parset) << std::endl; + } + + // Finally, add these universes into the vector + UniverseVals.push_back(tmp_vals_p1sigma); + UniverseVals.push_back(tmp_vals_m1sigma); + } // end loop over parsets (i) + } // pm1sigma + else if (mode=="multisim"){ + // multisim mode: parameter values sample within the given uncertainty for all parameters simultaneously + // Loop over universes j + for (unsigned j=0; j tmp_vals; + for (size_t i_parset=0; i_parsetfire(0.0,1.0); + tmp_vals[FitParNames.at(i_parset)] = FitParNominals.at(i_parset)+(FitParSigmas.at(i_parset)*r); + } // loop over parameters (i_parset) + // Now save this universe + UniverseVals.push_back(tmp_vals); + } // loop over Nsims (j) + } // multisim + else{ + // Anything else mode: Set parameters to user-defined nominal value + UniverseVals.push_back(theNominals); + } // any other mode + + fNsims = UniverseVals.size(); + if (fDebug) std::cout << "Running mode: " << mode <<". Nsims = " << fNsims << std::endl; +} + + +std::vector > +Geant4WeightCalc::GetWeight(art::Event& e) { + + // Get event/run/subrun numbers for output + run_num = e.run(); + subrun_num = e.subRun(); + event_num = e.id().event(); + + // Get MCParticles for each MCTruth in this event + art::Handle > truthHandle; + e.getByLabel(fMCTruthProducer, truthHandle); + const art::FindManyP truthParticles(truthHandle, e, fMCParticleProducer); + assert(truthParticles.isValid()); + + // Initialize the vector of event weights + std::vector > weight(truthHandle->size()); + // These two are just for saving to the output tree for fast validation + e_inel_weight.clear(); + e_inel_weight.resize(truthHandle->size()); + e_elastic_weight.clear(); + e_elastic_weight.resize(truthHandle->size()); + + // Loop over sets of MCTruth-associated particles + for (size_t itruth=0; itruth trajpoint_X; + std::vector trajpoint_Y; + std::vector trajpoint_Z; + std::vector trajpoint_PX; + std::vector trajpoint_PY; + std::vector trajpoint_PZ; + std::vector elastic_indices; + + //Get the list of processes from the true trajectory + const std::vector< std::pair< size_t, unsigned char > > & processes = p.Trajectory().TrajectoryProcesses(); + std::map< size_t, std::string > process_map; + for( auto it = processes.begin(); it != processes.end(); ++it ){ + process_map[ it->first ] = p.Trajectory().KeyToProcess( it->second ); + } + + for( size_t i = 0; i < p.NumberTrajectoryPoints(); ++i ){ + double X = p.Position(i).X(); + double Y = p.Position(i).Y(); + double Z = p.Position(i).Z(); + geo::Point_t testpoint1 { X, Y, Z }; + const TGeoMaterial* testmaterial1 = fGeometryService->Material( testpoint1 ); + //For now, just going to reweight the points within the LAr of the TPC + // TODO check if this is right + if ( !strcmp( testmaterial1->GetName(), "LAr" ) ){ + trajpoint_X.push_back( X ); + trajpoint_Y.push_back( Y ); + trajpoint_Z.push_back( Z ); + + trajpoint_PX.push_back( p.Px(i) ); + trajpoint_PY.push_back( p.Py(i) ); + trajpoint_PZ.push_back( p.Pz(i) ); + + auto itProc = process_map.find(i); + if( itProc != process_map.end() && itProc->second == "hadElastic" ){ + //Push back the index relative to the start of the reweightable steps + elastic_indices.push_back( trajpoint_X.size() - 1 ); + // if (fDebug) std::cout << "Elastic index: " << trajpoint_X.size() - 1 << std::endl; + } + + } + + } // end loop over trajectory points + + // Now find daughters of the MCP + std::vector daughter_PDGs; + std::vector daughter_IDs; + for( int i_mcp = 0; i_mcp < p.NumberDaughters(); i_mcp++ ){ + int daughterID = p.Daughter(i_mcp); + for (auto test_mcp : mcparticles){ + if (test_mcp->TrackId() == daughterID){ + int pid = test_mcp->PdgCode(); + daughter_PDGs.push_back(pid); + daughter_IDs.push_back( test_mcp->TrackId() ); + break; + } + } + } // end loop over daughters + + // --- Now that we have all the information about the track we need, here comes the reweighting part! --- // + + //Make a G4ReweightTraj -- This is the reweightable object + G4ReweightTraj theTraj(mcpID, p_PDG, 0, event_num, std::make_pair(0,0)); + + //Create its set of G4ReweightSteps and add them to the Traj (note: this needs to be done once per MCParticle but will be valid for all weight calculations) + std::vector< G4ReweightStep* > allSteps; + + size_t nSteps = trajpoint_PX.size(); + + if( nSteps < 2 ) continue; + + p_nElasticScatters = elastic_indices.size(); + for( size_t istep = 1; istep < nSteps; ++istep ){ + + // if( istep == trajpoint_PX.size() - 1 && std::find( elastic_indices.begin(), elastic_indices.end(), j ) != elastic_indices.end() ) + // std::cout << "Warning: last step an elastic process" << std::endl; + + std::string proc = "default"; + if( istep == trajpoint_PX.size() - 1 ) + proc = EndProcess; + else if( std::find( elastic_indices.begin(), elastic_indices.end(), istep ) != elastic_indices.end() ){ + proc = "hadElastic"; + } + + + double deltaX = ( trajpoint_X.at(istep) - trajpoint_X.at(istep-1) ); + double deltaY = ( trajpoint_Y.at(istep) - trajpoint_Y.at(istep-1) ); + double deltaZ = ( trajpoint_Z.at(istep) - trajpoint_Z.at(istep-1) ); + + double len = sqrt( + std::pow( deltaX, 2 ) + + std::pow( deltaY, 2 ) + + std::pow( deltaZ, 2 ) + ); + + double preStepP[3] = { + trajpoint_PX.at(istep-1)*1.e3, + trajpoint_PY.at(istep-1)*1.e3, + trajpoint_PZ.at(istep-1)*1.e3 + }; + + double postStepP[3] = { + trajpoint_PX.at(istep)*1.e3, + trajpoint_PY.at(istep)*1.e3, + trajpoint_PZ.at(istep)*1.e3 + }; + + if( istep == 1 ){ + theTraj.SetEnergy( sqrt( preStepP[0]*preStepP[0] + preStepP[1]*preStepP[1] + preStepP[2]*preStepP[2] + mass*mass ) ); + } + + G4ReweightStep * theStep = new G4ReweightStep( mcpID, p_PDG, 0, event_num, preStepP, postStepP, len, proc ); + theStep->SetDeltaX( deltaX ); + theStep->SetDeltaY( deltaY ); + theStep->SetDeltaZ( deltaZ ); + + theTraj.AddStep( theStep ); + + for( size_t k = 0; k < daughter_PDGs.size(); ++k ){ + theTraj.AddChild( new G4ReweightTraj(daughter_IDs[k], daughter_PDGs[k], mcpID, event_num, std::make_pair(0,0) ) ); + } + } // end loop over nSteps (istep) + p_track_length = theTraj.GetTotalLength(); + + p_init_momentum = sqrt( theTraj.GetEnergy()*theTraj.GetEnergy() - mass*mass ); + p_final_momentum = sqrt( + std::pow( theTraj.GetStep( theTraj.GetNSteps() - 1 )->GetPreStepPx(), 2 ) + + std::pow( theTraj.GetStep( theTraj.GetNSteps() - 1 )->GetPreStepPy(), 2 ) + + std::pow( theTraj.GetStep( theTraj.GetNSteps() - 1 )->GetPreStepPz(), 2 ) + ); + + std::vector< std::pair< double, int > > thin_slice_inelastic = ThinSliceBetheBloch( &theTraj, .5, mass , false); + std::vector< std::pair< double, int > > thin_slice_elastic = ThinSliceBetheBloch( &theTraj, .5, mass , true); + + p_energies_inel.clear(); + p_sliceInts_inel.clear(); + for( size_t islice = 0; islice < thin_slice_inelastic.size(); ++islice ){ + p_energies_inel.push_back( thin_slice_inelastic[islice].first ); + p_sliceInts_inel.push_back( thin_slice_inelastic[islice].second ); + } + + p_energies_el.clear(); + p_sliceInts_el.clear(); + for( size_t islice = 0; islice < thin_slice_elastic.size(); ++islice ){ + p_energies_el.push_back( thin_slice_elastic[islice].first ); + p_sliceInts_el.push_back( thin_slice_elastic[islice].second ); + } + + // Loop through universes (j) + for (size_t j=0; jSetNewVals(UniverseVals.at(j)); + theReweighter->SetNewHists(ParMaker->GetFSHists()); + theReweighter->SetNewElasticHists(ParMaker->GetElasticHist()); + + //Get the weight from the G4ReweightTraj + w = theReweighter->GetWeight( &theTraj ); + // Total weight is the product of track weights in the event + weight[itruth][j] *= std::max((float)0.0, w); + + // Do the same for elastic weight (should be 1 unless set to non-nominal ) + el_w = theReweighter->GetElasticWeight( &theTraj ); + weight[itruth][j] *= std::max((float)0.0,el_w); + + // just for the output tree + p_inel_weight[j] = w; + p_elastic_weight[j] = el_w; + e_inel_weight[itruth][j] *= std::max((float)0.0,w); + e_elastic_weight[itruth][j] *= std::max((float)0.0,el_w); + + if (fDebug){ + std::cout << " Universe " << j << ": "; + // std::cout << UniverseVals.at(j) << std::endl; + std::cout << " w = " << w << ", el_w = " << el_w << std::endl; + } + + + } // loop through universes (j) + + if (fDebug){ + std::cout << "PDG = " << p_PDG << std::endl; + std::cout << "inel weights by particle: "; + for (unsigned int j=0; jFill(); + } // loop over mcparticles (i) + if (fMakeOutputTrees) fOutTree_MCTruth->Fill(); + } // loop over sets of MCtruth-associated particles (itruth) + +return weight; + +} + +REGISTER_WEIGHTCALC(Geant4WeightCalc) + +} // namespace evwgh + +} // namespace sbn diff --git a/ups/product_deps b/ups/product_deps index 9e0fd280d..9586bb584 100644 --- a/ups/product_deps +++ b/ups/product_deps @@ -261,6 +261,7 @@ sbndata v01_07 - sbnobj v09_19_04 - systematicstools v01_04_04 - nusystematics v1_05_01 - +geant4reweight v01_20_08 - cetmodules v3_24_01 - only_for_build end_product_list #################################### @@ -317,11 +318,11 @@ end_product_list # #################################### -qualifier larsoft sbnobj sbnanaobj sbndaq_artdaq_core genie_xsec sbndata larcv2 systematicstools nusystematics notes -c14:debug c14:debug c14:debug c14:debug c14:s131:debug AR2320i00000:e1000:k250 -nq- c14:debug:p3915 c14:debug c14:debug -nq- -c14:prof c14:prof c14:prof c14:prof c14:s131:prof AR2320i00000:e1000:k250 -nq- c14:p3915:prof c14:prof c14:prof -nq- -e26:debug e26:debug e26:debug e26:debug e26:s131:debug AR2320i00000:e1000:k250 -nq- debug:e26:p3915 e26:debug e26:debug -nq- -e26:prof e26:prof e26:prof e26:prof e26:s131:prof AR2320i00000:e1000:k250 -nq- e26:p3915:prof e26:prof e26:prof -nq- +qualifier larsoft sbnobj sbnanaobj sbndaq_artdaq_core genie_xsec sbndata larcv2 systematicstools nusystematics geant4reweight notes +c14:debug c14:debug c14:debug c14:debug c14:s131:debug AR2320i00000:e1000:k250 -nq- c14:debug:p3915 c14:debug c14:debug c14:debug:s131 -nq- +c14:prof c14:prof c14:prof c14:prof c14:s131:prof AR2320i00000:e1000:k250 -nq- c14:p3915:prof c14:prof c14:prof c14:prof:s131 -nq- +e26:debug e26:debug e26:debug e26:debug e26:s131:debug AR2320i00000:e1000:k250 -nq- debug:e26:p3915 e26:debug e26:debug e26:debug:s131 -nq- +e26:prof e26:prof e26:prof e26:prof e26:s131:prof AR2320i00000:e1000:k250 -nq- e26:p3915:prof e26:prof e26:prof e26:prof:s131 -nq- end_qualifier_list #################################### From 5ede08a2adf1cb827865607c45b64be348fa02d5 Mon Sep 17 00:00:00 2001 From: Steven Gardiner Date: Tue, 25 Jul 2023 17:25:52 -0500 Subject: [PATCH 02/16] Some tweaks to adjust for larsim EventWeight (used by uB) vs. SBNEventWeight differences. Also make some adjustments to get it to build against geant4reweight v01_00_03e -q e20:s120a:prof --- .../Calculators/Geant4/CMakeLists.txt | 1 + .../Calculators/Geant4/Geant4WeightCalc.cxx | 463 ++++++++---------- ups/product_deps | 2 +- 3 files changed, 216 insertions(+), 250 deletions(-) diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt b/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt index e1289208d..f05f8d274 100644 --- a/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt +++ b/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt @@ -14,6 +14,7 @@ art_make_library( cetlib_except::cetlib_except ReweightBaseLib PropBaseLib + ROOT::Tree ) install_headers() diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx index f4ba26774..8b4736aa0 100644 --- a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx +++ b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx @@ -52,7 +52,7 @@ class Geant4WeightCalc : public WeightCalc { void Configure(fhicl::ParameterSet const& p, CLHEP::HepRandomEngine& engine); - std::vector > GetWeight(art::Event& e); + std::vector GetWeight(art::Event& e, size_t inu); private: std::string fMCParticleProducer; //!< Label for the MCParticle producer @@ -87,17 +87,18 @@ class Geant4WeightCalc : public WeightCalc { int p_nElasticScatters; //!< Variables for by-particle output tree std::vector p_inel_weight; //!< Variables for by-particle output tree std::vector p_elastic_weight; //!< Variables for by-particle output - std::vector > e_inel_weight; //!< Variables for by-event output tree - std::vector > e_elastic_weight; //!< Variables for by-event output tree bool fDebug; //!< print debug statements DECLARE_WEIGHTCALC(Geant4WeightCalc) }; +} // namespace evwgh + +} // namespace sbn -void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, - CLHEP::HepRandomEngine& engine) +void sbn::evwgh::Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, + CLHEP::HepRandomEngine& engine) { std::cout << "Using Geant4WeightCalc for reinteraction weights" << std::endl; @@ -157,9 +158,6 @@ void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, fOutTree_MCTruth->Branch("subrun_num",&subrun_num); fOutTree_MCTruth->Branch("UniverseVals", &UniverseVals); fOutTree_MCTruth->Branch("pdg_to_reweight",&fPdg); - fOutTree_MCTruth->Branch("inelastic_weight",&e_inel_weight); - fOutTree_MCTruth->Branch("elastic_weight",&e_elastic_weight); - } @@ -227,8 +225,8 @@ void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, } -std::vector > -Geant4WeightCalc::GetWeight(art::Event& e) { +std::vector +sbn::evwgh::Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { // Get event/run/subrun numbers for output run_num = e.run(); @@ -242,283 +240,250 @@ Geant4WeightCalc::GetWeight(art::Event& e) { assert(truthParticles.isValid()); // Initialize the vector of event weights - std::vector > weight(truthHandle->size()); - // These two are just for saving to the output tree for fast validation - e_inel_weight.clear(); - e_inel_weight.resize(truthHandle->size()); - e_elastic_weight.clear(); - e_elastic_weight.resize(truthHandle->size()); + std::vector weight; - // Loop over sets of MCTruth-associated particles - for (size_t itruth=0; itruth trajpoint_X; - std::vector trajpoint_Y; - std::vector trajpoint_Z; - std::vector trajpoint_PX; - std::vector trajpoint_PY; - std::vector trajpoint_PZ; - std::vector elastic_indices; - - //Get the list of processes from the true trajectory - const std::vector< std::pair< size_t, unsigned char > > & processes = p.Trajectory().TrajectoryProcesses(); - std::map< size_t, std::string > process_map; - for( auto it = processes.begin(); it != processes.end(); ++it ){ - process_map[ it->first ] = p.Trajectory().KeyToProcess( it->second ); - } - - for( size_t i = 0; i < p.NumberTrajectoryPoints(); ++i ){ - double X = p.Position(i).X(); - double Y = p.Position(i).Y(); - double Z = p.Position(i).Z(); - geo::Point_t testpoint1 { X, Y, Z }; - const TGeoMaterial* testmaterial1 = fGeometryService->Material( testpoint1 ); - //For now, just going to reweight the points within the LAr of the TPC - // TODO check if this is right - if ( !strcmp( testmaterial1->GetName(), "LAr" ) ){ - trajpoint_X.push_back( X ); - trajpoint_Y.push_back( Y ); - trajpoint_Z.push_back( Z ); - - trajpoint_PX.push_back( p.Px(i) ); - trajpoint_PY.push_back( p.Py(i) ); - trajpoint_PZ.push_back( p.Pz(i) ); - - auto itProc = process_map.find(i); - if( itProc != process_map.end() && itProc->second == "hadElastic" ){ - //Push back the index relative to the start of the reweightable steps - elastic_indices.push_back( trajpoint_X.size() - 1 ); - // if (fDebug) std::cout << "Elastic index: " << trajpoint_X.size() - 1 << std::endl; - } + p_inel_weight.clear(); + p_inel_weight.resize(fNsims,1.0); + p_elastic_weight.clear(); + p_elastic_weight.resize(fNsims,1.0); + p_track_length=-9999; + p_PDG=-9999; + p_final_proc="dummy"; + p_init_momentum=-9999; + p_final_momentum=-9999; + p_energies_el.clear(); + p_sliceInts_el.clear(); + p_energies_inel.clear(); + p_sliceInts_inel.clear(); + p_nElasticScatters=-9999; + + + const simb::MCParticle& p = *mcparticles.at(i); + p_PDG = p.PdgCode(); + int mcpID = p.TrackId(); + std::string EndProcess = p.EndProcess(); + + double mass = 0.; + if( TMath::Abs(p_PDG) == 211 ) mass = 139.57; + else if( p_PDG == 2212 ) mass = 938.28; + + // We only want to record weights for one type of particle (defined by fPDG from the fcl file), so skip other particles + if (p_PDG == fPdg){ + // Get GEANT trajectory points: weighting will depend on position and momentum at each trajectory point so calculate those + std::vector trajpoint_X; + std::vector trajpoint_Y; + std::vector trajpoint_Z; + std::vector trajpoint_PX; + std::vector trajpoint_PY; + std::vector trajpoint_PZ; + std::vector elastic_indices; + + //Get the list of processes from the true trajectory + const std::vector< std::pair< size_t, unsigned char > > & processes = p.Trajectory().TrajectoryProcesses(); + std::map< size_t, std::string > process_map; + for( auto it = processes.begin(); it != processes.end(); ++it ){ + process_map[ it->first ] = p.Trajectory().KeyToProcess( it->second ); + } + for( size_t i = 0; i < p.NumberTrajectoryPoints(); ++i ){ + double X = p.Position(i).X(); + double Y = p.Position(i).Y(); + double Z = p.Position(i).Z(); + geo::Point_t testpoint1 { X, Y, Z }; + const TGeoMaterial* testmaterial1 = fGeometryService->Material( testpoint1 ); + //For now, just going to reweight the points within the LAr of the TPC + // TODO check if this is right + if ( !strcmp( testmaterial1->GetName(), "LAr" ) ){ + trajpoint_X.push_back( X ); + trajpoint_Y.push_back( Y ); + trajpoint_Z.push_back( Z ); + + trajpoint_PX.push_back( p.Px(i) ); + trajpoint_PY.push_back( p.Py(i) ); + trajpoint_PZ.push_back( p.Pz(i) ); + + auto itProc = process_map.find(i); + if( itProc != process_map.end() && itProc->second == "hadElastic" ){ + //Push back the index relative to the start of the reweightable steps + elastic_indices.push_back( trajpoint_X.size() - 1 ); + // if (fDebug) std::cout << "Elastic index: " << trajpoint_X.size() - 1 << std::endl; } - } // end loop over trajectory points - - // Now find daughters of the MCP - std::vector daughter_PDGs; - std::vector daughter_IDs; - for( int i_mcp = 0; i_mcp < p.NumberDaughters(); i_mcp++ ){ - int daughterID = p.Daughter(i_mcp); - for (auto test_mcp : mcparticles){ - if (test_mcp->TrackId() == daughterID){ - int pid = test_mcp->PdgCode(); - daughter_PDGs.push_back(pid); - daughter_IDs.push_back( test_mcp->TrackId() ); - break; - } + } + + } // end loop over trajectory points + + // Now find daughters of the MCP + std::vector daughter_PDGs; + std::vector daughter_IDs; + for( int i_mcp = 0; i_mcp < p.NumberDaughters(); i_mcp++ ){ + int daughterID = p.Daughter(i_mcp); + for (auto test_mcp : mcparticles){ + if (test_mcp->TrackId() == daughterID){ + int pid = test_mcp->PdgCode(); + daughter_PDGs.push_back(pid); + daughter_IDs.push_back( test_mcp->TrackId() ); + break; } - } // end loop over daughters + } + } // end loop over daughters - // --- Now that we have all the information about the track we need, here comes the reweighting part! --- // + // --- Now that we have all the information about the track we need, here comes the reweighting part! --- // - //Make a G4ReweightTraj -- This is the reweightable object - G4ReweightTraj theTraj(mcpID, p_PDG, 0, event_num, std::make_pair(0,0)); + //Make a G4ReweightTraj -- This is the reweightable object + G4ReweightTraj theTraj(mcpID, p_PDG, 0, event_num, std::make_pair(0,0)); - //Create its set of G4ReweightSteps and add them to the Traj (note: this needs to be done once per MCParticle but will be valid for all weight calculations) - std::vector< G4ReweightStep* > allSteps; + //Create its set of G4ReweightSteps and add them to the Traj (note: this needs to be done once per MCParticle but will be valid for all weight calculations) + std::vector< G4ReweightStep* > allSteps; - size_t nSteps = trajpoint_PX.size(); + size_t nSteps = trajpoint_PX.size(); - if( nSteps < 2 ) continue; + if( nSteps < 2 ) continue; - p_nElasticScatters = elastic_indices.size(); - for( size_t istep = 1; istep < nSteps; ++istep ){ + p_nElasticScatters = elastic_indices.size(); + for( size_t istep = 1; istep < nSteps; ++istep ){ - // if( istep == trajpoint_PX.size() - 1 && std::find( elastic_indices.begin(), elastic_indices.end(), j ) != elastic_indices.end() ) - // std::cout << "Warning: last step an elastic process" << std::endl; + // if( istep == trajpoint_PX.size() - 1 && std::find( elastic_indices.begin(), elastic_indices.end(), j ) != elastic_indices.end() ) + // std::cout << "Warning: last step an elastic process" << std::endl; - std::string proc = "default"; - if( istep == trajpoint_PX.size() - 1 ) - proc = EndProcess; - else if( std::find( elastic_indices.begin(), elastic_indices.end(), istep ) != elastic_indices.end() ){ - proc = "hadElastic"; - } + std::string proc = "default"; + if( istep == trajpoint_PX.size() - 1 ) + proc = EndProcess; + else if( std::find( elastic_indices.begin(), elastic_indices.end(), istep ) != elastic_indices.end() ){ + proc = "hadElastic"; + } - double deltaX = ( trajpoint_X.at(istep) - trajpoint_X.at(istep-1) ); - double deltaY = ( trajpoint_Y.at(istep) - trajpoint_Y.at(istep-1) ); - double deltaZ = ( trajpoint_Z.at(istep) - trajpoint_Z.at(istep-1) ); + double deltaX = ( trajpoint_X.at(istep) - trajpoint_X.at(istep-1) ); + double deltaY = ( trajpoint_Y.at(istep) - trajpoint_Y.at(istep-1) ); + double deltaZ = ( trajpoint_Z.at(istep) - trajpoint_Z.at(istep-1) ); - double len = sqrt( - std::pow( deltaX, 2 ) + - std::pow( deltaY, 2 ) + - std::pow( deltaZ, 2 ) - ); + double len = sqrt( + std::pow( deltaX, 2 ) + + std::pow( deltaY, 2 ) + + std::pow( deltaZ, 2 ) + ); - double preStepP[3] = { - trajpoint_PX.at(istep-1)*1.e3, - trajpoint_PY.at(istep-1)*1.e3, - trajpoint_PZ.at(istep-1)*1.e3 - }; + double preStepP[3] = { + trajpoint_PX.at(istep-1)*1.e3, + trajpoint_PY.at(istep-1)*1.e3, + trajpoint_PZ.at(istep-1)*1.e3 + }; - double postStepP[3] = { - trajpoint_PX.at(istep)*1.e3, - trajpoint_PY.at(istep)*1.e3, - trajpoint_PZ.at(istep)*1.e3 - }; + double postStepP[3] = { + trajpoint_PX.at(istep)*1.e3, + trajpoint_PY.at(istep)*1.e3, + trajpoint_PZ.at(istep)*1.e3 + }; - if( istep == 1 ){ - theTraj.SetEnergy( sqrt( preStepP[0]*preStepP[0] + preStepP[1]*preStepP[1] + preStepP[2]*preStepP[2] + mass*mass ) ); - } + if( istep == 1 ){ + theTraj.SetEnergy( sqrt( preStepP[0]*preStepP[0] + preStepP[1]*preStepP[1] + preStepP[2]*preStepP[2] + mass*mass ) ); + } - G4ReweightStep * theStep = new G4ReweightStep( mcpID, p_PDG, 0, event_num, preStepP, postStepP, len, proc ); - theStep->SetDeltaX( deltaX ); - theStep->SetDeltaY( deltaY ); - theStep->SetDeltaZ( deltaZ ); + G4ReweightStep * theStep = new G4ReweightStep( mcpID, p_PDG, 0, event_num, preStepP, postStepP, len, proc ); + theStep->SetDeltaX( deltaX ); + theStep->SetDeltaY( deltaY ); + theStep->SetDeltaZ( deltaZ ); - theTraj.AddStep( theStep ); + theTraj.AddStep( theStep ); - for( size_t k = 0; k < daughter_PDGs.size(); ++k ){ - theTraj.AddChild( new G4ReweightTraj(daughter_IDs[k], daughter_PDGs[k], mcpID, event_num, std::make_pair(0,0) ) ); - } - } // end loop over nSteps (istep) - p_track_length = theTraj.GetTotalLength(); - - p_init_momentum = sqrt( theTraj.GetEnergy()*theTraj.GetEnergy() - mass*mass ); - p_final_momentum = sqrt( - std::pow( theTraj.GetStep( theTraj.GetNSteps() - 1 )->GetPreStepPx(), 2 ) + - std::pow( theTraj.GetStep( theTraj.GetNSteps() - 1 )->GetPreStepPy(), 2 ) + - std::pow( theTraj.GetStep( theTraj.GetNSteps() - 1 )->GetPreStepPz(), 2 ) - ); + for( size_t k = 0; k < daughter_PDGs.size(); ++k ){ + theTraj.AddChild( new G4ReweightTraj(daughter_IDs[k], daughter_PDGs[k], mcpID, event_num, std::make_pair(0,0) ) ); + } + } // end loop over nSteps (istep) + p_track_length = theTraj.GetTotalLength(); - std::vector< std::pair< double, int > > thin_slice_inelastic = ThinSliceBetheBloch( &theTraj, .5, mass , false); - std::vector< std::pair< double, int > > thin_slice_elastic = ThinSliceBetheBloch( &theTraj, .5, mass , true); + p_init_momentum = sqrt( theTraj.GetEnergy()*theTraj.GetEnergy() - mass*mass ); + p_final_momentum = sqrt( + std::pow( theTraj.GetStep( theTraj.GetNSteps() - 1 )->GetPreStepPx(), 2 ) + + std::pow( theTraj.GetStep( theTraj.GetNSteps() - 1 )->GetPreStepPy(), 2 ) + + std::pow( theTraj.GetStep( theTraj.GetNSteps() - 1 )->GetPreStepPz(), 2 ) + ); - p_energies_inel.clear(); - p_sliceInts_inel.clear(); - for( size_t islice = 0; islice < thin_slice_inelastic.size(); ++islice ){ - p_energies_inel.push_back( thin_slice_inelastic[islice].first ); - p_sliceInts_inel.push_back( thin_slice_inelastic[islice].second ); - } + std::vector< std::pair< double, int > > thin_slice_inelastic = ThinSliceBetheBloch( &theTraj, .5, mass , false); + std::vector< std::pair< double, int > > thin_slice_elastic = ThinSliceBetheBloch( &theTraj, .5, mass , true); - p_energies_el.clear(); - p_sliceInts_el.clear(); - for( size_t islice = 0; islice < thin_slice_elastic.size(); ++islice ){ - p_energies_el.push_back( thin_slice_elastic[islice].first ); - p_sliceInts_el.push_back( thin_slice_elastic[islice].second ); - } + p_energies_inel.clear(); + p_sliceInts_inel.clear(); + for( size_t islice = 0; islice < thin_slice_inelastic.size(); ++islice ){ + p_energies_inel.push_back( thin_slice_inelastic[islice].first ); + p_sliceInts_inel.push_back( thin_slice_inelastic[islice].second ); + } - // Loop through universes (j) - for (size_t j=0; jSetNewVals(UniverseVals.at(j)); - theReweighter->SetNewHists(ParMaker->GetFSHists()); - theReweighter->SetNewElasticHists(ParMaker->GetElasticHist()); - - //Get the weight from the G4ReweightTraj - w = theReweighter->GetWeight( &theTraj ); - // Total weight is the product of track weights in the event - weight[itruth][j] *= std::max((float)0.0, w); - - // Do the same for elastic weight (should be 1 unless set to non-nominal ) - el_w = theReweighter->GetElasticWeight( &theTraj ); - weight[itruth][j] *= std::max((float)0.0,el_w); - - // just for the output tree - p_inel_weight[j] = w; - p_elastic_weight[j] = el_w; - e_inel_weight[itruth][j] *= std::max((float)0.0,w); - e_elastic_weight[itruth][j] *= std::max((float)0.0,el_w); - - if (fDebug){ - std::cout << " Universe " << j << ": "; - // std::cout << UniverseVals.at(j) << std::endl; - std::cout << " w = " << w << ", el_w = " << el_w << std::endl; - } + p_energies_el.clear(); + p_sliceInts_el.clear(); + for( size_t islice = 0; islice < thin_slice_elastic.size(); ++islice ){ + p_energies_el.push_back( thin_slice_elastic[islice].first ); + p_sliceInts_el.push_back( thin_slice_elastic[islice].second ); + } + + // Loop through universes (j) + for (size_t j=0; jSetNewVals(UniverseVals.at(j)); + theReweighter->SetNewHists(ParMaker->GetFSHists()); + theReweighter->SetNewElasticHists(ParMaker->GetElasticHist()); + //Get the weight from the G4ReweightTraj + w = theReweighter->GetWeight( &theTraj ); + // Total weight is the product of track weights in the event + weight[j] *= std::max((float)0.0, w); - } // loop through universes (j) + // Do the same for elastic weight (should be 1 unless set to non-nominal ) + el_w = theReweighter->GetElasticWeight( &theTraj ); + weight[j] *= std::max((float)0.0,el_w); + + // just for the output tree + p_inel_weight[j] = w; + p_elastic_weight[j] = el_w; if (fDebug){ - std::cout << "PDG = " << p_PDG << std::endl; - std::cout << "inel weights by particle: "; - for (unsigned int j=0; jFill(); + } // loop over mcparticles (i) + if (fMakeOutputTrees) fOutTree_MCTruth->Fill(); -} // namespace evwgh +return weight; -} // namespace sbn +} + +REGISTER_WEIGHTCALC(sbn::evwgh::Geant4WeightCalc) diff --git a/ups/product_deps b/ups/product_deps index 9586bb584..b54043517 100644 --- a/ups/product_deps +++ b/ups/product_deps @@ -248,7 +248,7 @@ libdir fq_dir lib # cetmodules), an entry for cetmodules is required. # # * It is an error for more than one non-( == "-default-") -# entry to match for a given product. + entry to match for a given product. # #################################### product version qual flags From c30d8d153982f1476ec18e2183a3bbfa620ee64f Mon Sep 17 00:00:00 2001 From: Andrew Mastbaum Date: Wed, 26 Jul 2023 15:59:01 -0500 Subject: [PATCH 03/16] geant4reweight interface updates + generic sbn fcls --- sbncode/SBNEventWeight/App/CMakeLists.txt | 7 +- .../Calculators/Geant4/Geant4WeightCalc.cxx | 14 +-- sbncode/SBNEventWeight/jobs/CMakeLists.txt | 1 + .../SBNEventWeight/jobs/geant4/CMakeLists.txt | 6 ++ .../jobs/geant4/eventweight_g4rwt_reint.fcl | 99 +++++++++++++++++++ .../geant4/piminus_reweight_parameters.fcl | 66 +++++++++++++ .../geant4/piplus_reweight_parameters.fcl | 66 +++++++++++++ .../geant4/proton_reweight_parameters.fcl | 23 +++++ 8 files changed, 273 insertions(+), 9 deletions(-) create mode 100644 sbncode/SBNEventWeight/jobs/geant4/CMakeLists.txt create mode 100644 sbncode/SBNEventWeight/jobs/geant4/eventweight_g4rwt_reint.fcl create mode 100644 sbncode/SBNEventWeight/jobs/geant4/piminus_reweight_parameters.fcl create mode 100644 sbncode/SBNEventWeight/jobs/geant4/piplus_reweight_parameters.fcl create mode 100644 sbncode/SBNEventWeight/jobs/geant4/proton_reweight_parameters.fcl diff --git a/sbncode/SBNEventWeight/App/CMakeLists.txt b/sbncode/SBNEventWeight/App/CMakeLists.txt index 2989659b7..335a19557 100644 --- a/sbncode/SBNEventWeight/App/CMakeLists.txt +++ b/sbncode/SBNEventWeight/App/CMakeLists.txt @@ -1,8 +1,11 @@ cet_build_plugin(SBNEventWeight art::module LIBRARIES - sbnobj::Common_SBNEventWeight sbncode_SBNEventWeight_Base + sbnobj::Common_SBNEventWeight + sbncode_SBNEventWeight_Base sbncode_SBNEventWeight_Calculators_CrossSection - sbncode_SBNEventWeight_Calculators_BNBFlux nusimdata::SimulationBase + sbncode_SBNEventWeight_Calculators_BNBFlux + sbncode_SBNEventWeight_Calculators_Geant4 + nusimdata::SimulationBase ROOT::Geom) cet_build_plugin(SystToolsEventWeight art::module diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx index 8b4736aa0..0ed593b04 100644 --- a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx +++ b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx @@ -93,11 +93,8 @@ class Geant4WeightCalc : public WeightCalc { DECLARE_WEIGHTCALC(Geant4WeightCalc) }; -} // namespace evwgh - -} // namespace sbn -void sbn::evwgh::Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, +void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, CLHEP::HepRandomEngine& engine) { std::cout << "Using Geant4WeightCalc for reinteraction weights" << std::endl; @@ -225,8 +222,7 @@ void sbn::evwgh::Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, } -std::vector -sbn::evwgh::Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { +std::vector Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { // Get event/run/subrun numbers for output run_num = e.run(); @@ -486,4 +482,8 @@ return weight; } -REGISTER_WEIGHTCALC(sbn::evwgh::Geant4WeightCalc) +REGISTER_WEIGHTCALC(Geant4WeightCalc) + +} // namespace evwgh + +} // namespace sbn diff --git a/sbncode/SBNEventWeight/jobs/CMakeLists.txt b/sbncode/SBNEventWeight/jobs/CMakeLists.txt index e36b8b0cc..b5aae6724 100644 --- a/sbncode/SBNEventWeight/jobs/CMakeLists.txt +++ b/sbncode/SBNEventWeight/jobs/CMakeLists.txt @@ -5,4 +5,5 @@ install_source() add_subdirectory(SystTools) add_subdirectory(genie) add_subdirectory(flux) +add_subdirectory(geant4) diff --git a/sbncode/SBNEventWeight/jobs/geant4/CMakeLists.txt b/sbncode/SBNEventWeight/jobs/geant4/CMakeLists.txt new file mode 100644 index 000000000..589bc695a --- /dev/null +++ b/sbncode/SBNEventWeight/jobs/geant4/CMakeLists.txt @@ -0,0 +1,6 @@ +install_fhicl() + +FILE(GLOB fcl_files *.fcl) + +install_source( EXTRAS ${fcl_files} ) + diff --git a/sbncode/SBNEventWeight/jobs/geant4/eventweight_g4rwt_reint.fcl b/sbncode/SBNEventWeight/jobs/geant4/eventweight_g4rwt_reint.fcl new file mode 100644 index 000000000..d08990c2c --- /dev/null +++ b/sbncode/SBNEventWeight/jobs/geant4/eventweight_g4rwt_reint.fcl @@ -0,0 +1,99 @@ +########################################################## +## Hadron reinteraction uncertainties +## +## References: +## +## * K. Duffy, Pion Secondary Interaction Systematics (from GEANTReweight), DocDB 25379 +## * J. Calcutt, GEANTReweight documentation DocDB 25084, repository https://cdcvs.fnal.gov/redmine/projects/geant4reweight/repository +## +## From MicroBooNE, ubsim UBOONE_SUITE_v08_00_00_69 +## (uBooNE author: Kirsty Duffy (kduffy@fnal.gov)) +## +########################################################## + +#include "piplus_reweight_parameters.fcl" +#include "piminus_reweight_parameters.fcl" +#include "proton_reweight_parameters.fcl" + +geant4weight_sbn: { + module_type: "SBNEventWeight" + AllowMissingTruth: true + min_weight: 0.0 + max_weight: 100 + + weight_functions_reint: [ reinteractions_piplus, reinteractions_piminus, reinteractions_proton ] + + reinteractions_piplus: { + type: Geant4 + 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 + pdg_to_reweight: 211 + debug: false + } + + reinteractions_piminus: { + type: Geant4 + 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 + pdg_to_reweight: -211 + debug: false + } + + reinteractions_proton: { + type: Geant4 + 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 + pdg_to_reweight: 2212 + debug: false + } +} + +### +# Reweighting parameters should be defined as +# +# TheParameters: [ +# { +# Cut: "reac" +# Name: "fReacLow" +# Range: [10., 200.] +# Nominal: 1.0 +# Sigma: 0.3 +# }, +# {...} +# ] +# +# - Range defines the energy range over which that parameter has an effect (MeV) +# - Nominal is not used in "multisim" mode. In "pm1sigma" mode it defines the +# nominal (around which you calculate +/- 1 sigma variations). In any other +# mode, the parameter is set to the value under Nominal (so it's not really a +# "nominal" value in that case, just the set value) +# - Sigma defines the 1-sigma range for multisims and unisims. Currently the +# code can only accept a single value for sigma, giving symmetric +# uncertainty bounds: nominal-sigma and nominal+sigma +# - Cut must be one of the following: +# "reac" <- total inelastic scattering cross section +# "abs" <- absorption cross section (pi+ and pi- only at the moment -- Oct 2019) +# "cex" <- charge exchange cross section (pi+ and pi- only at the moment -- Oct 2019) +# "dcex" <- double charge exchange cross section (pi+ and pi- only at the moment -- Oct 2019) +# "prod" <- pion production cross section (pi+ and pi- only at the moment -- Oct 2019) +# "inel" <- quasi-elastic inelastic scattering cross section (pi+ and pi- only at the moment -- Oct 2019) +# "elast" <- total elastic scattering cross section +# - Name can be anything you like (but should uniquely identify that parameter/ +# variation) +### + diff --git a/sbncode/SBNEventWeight/jobs/geant4/piminus_reweight_parameters.fcl b/sbncode/SBNEventWeight/jobs/geant4/piminus_reweight_parameters.fcl new file mode 100644 index 000000000..b64756415 --- /dev/null +++ b/sbncode/SBNEventWeight/jobs/geant4/piminus_reweight_parameters.fcl @@ -0,0 +1,66 @@ +BEGIN_PROLOG + +# From MicroBooNE, ubsim UBOONE_SUITE_v08_00_00_69 + +PiMinusParameters: [ + { + Name: "fPiMinusReacLow" + Cut: "reac" + Range: [10., 200.] + Nominal: 1 + Sigma: 0.2 + } , + { + Name: "fPiMinusReacHigh" + Cut: "reac" + Range: [700., 2005.] + Nominal: 1.0 + Sigma: 0.2 + } , + + { + Name: "fPiMinusAbs" + Cut: "abs" + # Range: [200.0, 700.00] + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + }, + + { + Name: "fPiMinusCex" + Cut: "cex" + # Range: [200.0, 700.00] + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + }, + + { + Name: "fPiMinusDCex" + Cut: "dcex" + # Range: [200.0, 700.00] + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + }, + + { + Name: "fPiMinusPiProd" + Cut: "prod" + # Range: [200.0, 700.00] + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + }, + + { + Name: "fPiMinusElast" + Cut: "elast" + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + } +] + +END_PROLOG diff --git a/sbncode/SBNEventWeight/jobs/geant4/piplus_reweight_parameters.fcl b/sbncode/SBNEventWeight/jobs/geant4/piplus_reweight_parameters.fcl new file mode 100644 index 000000000..c9cfc3e23 --- /dev/null +++ b/sbncode/SBNEventWeight/jobs/geant4/piplus_reweight_parameters.fcl @@ -0,0 +1,66 @@ +BEGIN_PROLOG + +# From MicroBooNE, ubsim UBOONE_SUITE_v08_00_00_69 + +PiPlusParameters: [ + { + Name: "fPiPlusReacLow" + Cut: "reac" + Range: [10., 200.] + Nominal: 1 + Sigma: 0.2 + } , + { + Name: "fPiPlusReacHigh" + Cut: "reac" + Range: [700., 2005.] + Nominal: 1.0 + Sigma: 0.2 + } , + + { + Name: "fPiPlusAbs" + Cut: "abs" + # Range: [200.0, 700.00] + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + }, + + { + Name: "fPiPlusCex" + Cut: "cex" + # Range: [200.0, 700.00] + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + }, + + { + Name: "fPiPlusDCex" + Cut: "dcex" + # Range: [200.0, 700.00] + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + }, + + { + Name: "fPiPlusPiProd" + Cut: "prod" + # Range: [200.0, 700.00] + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + }, + + { + Name: "fPiPlusElast" + Cut: "elast" + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + } +] + +END_PROLOG diff --git a/sbncode/SBNEventWeight/jobs/geant4/proton_reweight_parameters.fcl b/sbncode/SBNEventWeight/jobs/geant4/proton_reweight_parameters.fcl new file mode 100644 index 000000000..13113db75 --- /dev/null +++ b/sbncode/SBNEventWeight/jobs/geant4/proton_reweight_parameters.fcl @@ -0,0 +1,23 @@ +BEGIN_PROLOG + +# From MicroBooNE, ubsim UBOONE_SUITE_v08_00_00_69 + +ProtonParameters: [ + { + Name: "fProtonReac" + Cut: "reac" + Range: [10., 2005.] + Nominal: 1 + Sigma: 0.2 + }, + + { + Name: "fProtonElast" + Cut: "elast" + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + } +] + +END_PROLOG From 82c781cbdeb59c339d4b0d615d74e0c9581f1cfc Mon Sep 17 00:00:00 2001 From: Andrew Mastbaum Date: Thu, 27 Jul 2023 13:40:59 -0500 Subject: [PATCH 04/16] g4reweight bugfix and temporary parameter storage * Placeholder for sampled parameter storage, for CAF compatibility. * Renamed fcl for consistency --- fcl/caf/cafmaker_common_defs.fcl | 2 ++ .../SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx | 7 +++++-- ...ntweight_g4rwt_reint.fcl => eventweight_geant4_sbn.fcl} | 6 +++++- 3 files changed, 12 insertions(+), 3 deletions(-) rename sbncode/SBNEventWeight/jobs/geant4/{eventweight_g4rwt_reint.fcl => eventweight_geant4_sbn.fcl} (98%) diff --git a/fcl/caf/cafmaker_common_defs.fcl b/fcl/caf/cafmaker_common_defs.fcl index 66fcdf7ab..b4d45abbf 100644 --- a/fcl/caf/cafmaker_common_defs.fcl +++ b/fcl/caf/cafmaker_common_defs.fcl @@ -1,3 +1,4 @@ +#include "eventweight_geant4_sbn.fcl" #include "eventweight_genie_sbn.fcl" #include "eventweight_flux_sbn.fcl" @@ -15,6 +16,7 @@ cafmaker_common_producers: { rns: { module_type: "RandomNumberSaver" } genieweight: @local::sbn_eventweight_genie fluxweight: @local::sbn_eventweight_flux + geant4weight: @local::sbn_eventweight_geant4 } END_PROLOG diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx index 0ed593b04..de2dda819 100644 --- a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx +++ b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx @@ -165,6 +165,8 @@ void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, std::vector FitParSigmas; std::map theNominals; + fParameterSet.Configure(GetFullName(), mode, fNsims); + for (size_t i_parset=0; i_parset("Name"); @@ -176,6 +178,8 @@ void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, FitParSigmas.push_back(theSigma); theNominals[theName] = theNominal; + + fParameterSet.AddParameter(theName, theSigma); } if (mode=="pm1sigma"){ @@ -240,7 +244,7 @@ std::vector Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { // Initialize weight vector for this MCTruth weight.clear(); - weight.resize(1.0); + weight.resize(fNsims, 1.0); // Loop over MCParticles in the event auto const& mcparticles = truthParticles.at(itruth); @@ -431,7 +435,6 @@ std::vector Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { ParMaker->SetNewVals(UniverseVals.at(j)); theReweighter->SetNewHists(ParMaker->GetFSHists()); theReweighter->SetNewElasticHists(ParMaker->GetElasticHist()); - //Get the weight from the G4ReweightTraj w = theReweighter->GetWeight( &theTraj ); // Total weight is the product of track weights in the event diff --git a/sbncode/SBNEventWeight/jobs/geant4/eventweight_g4rwt_reint.fcl b/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl similarity index 98% rename from sbncode/SBNEventWeight/jobs/geant4/eventweight_g4rwt_reint.fcl rename to sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl index d08990c2c..e45065880 100644 --- a/sbncode/SBNEventWeight/jobs/geant4/eventweight_g4rwt_reint.fcl +++ b/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl @@ -15,7 +15,9 @@ #include "piminus_reweight_parameters.fcl" #include "proton_reweight_parameters.fcl" -geant4weight_sbn: { +BEGIN_PROLOG + +sbn_eventweight_geant4: { module_type: "SBNEventWeight" AllowMissingTruth: true min_weight: 0.0 @@ -97,3 +99,5 @@ geant4weight_sbn: { # variation) ### +END_PROLOG + From 47810c44af32287061351cf050ca067b3dbf803d Mon Sep 17 00:00:00 2001 From: gputnam Date: Sat, 21 Oct 2023 20:50:57 -0500 Subject: [PATCH 05/16] Include Geant4 libraries for buiding --- sbncode/SBNEventWeight/CMakeLists.txt | 3 +++ sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt | 1 + 2 files changed, 4 insertions(+) diff --git a/sbncode/SBNEventWeight/CMakeLists.txt b/sbncode/SBNEventWeight/CMakeLists.txt index de59387f6..9f4018c28 100644 --- a/sbncode/SBNEventWeight/CMakeLists.txt +++ b/sbncode/SBNEventWeight/CMakeLists.txt @@ -1,3 +1,6 @@ +include_directories ( $ENV{GEANT4REWEIGHT_INC} ) +link_directories ( $ENV{GEANT4REWEIGHT_LIB} ) + add_subdirectory(Base) add_subdirectory(Calculators) add_subdirectory(App) diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt b/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt index f05f8d274..568ecc7f6 100644 --- a/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt +++ b/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt @@ -1,5 +1,6 @@ include_directories ( $ENV{GEANT4_FQ_DIR}/include ) include_directories ( $ENV{GEANT4REWEIGHT_INC} ) +link_directories ( $ENV{GEANT4REWEIGHT_LIB} ) art_make_library( LIBRARY_NAME sbncode_SBNEventWeight_Calculators_Geant4 From 0509def6c3c4ca29537c21364e96c7cfc82d03eb Mon Sep 17 00:00:00 2001 From: Jack Smedley Date: Tue, 23 Jul 2024 19:40:50 -0500 Subject: [PATCH 06/16] Added G4ReweightManager and collapsed elastic and inelastic weights into one, successfully compiled --- .../Calculators/Geant4/Geant4WeightCalc.cxx | 49 +++++++++++++------ 1 file changed, 34 insertions(+), 15 deletions(-) diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx index de2dda819..401c08896 100644 --- a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx +++ b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx @@ -33,6 +33,7 @@ #include "sbncode/SBNEventWeight/Base/WeightCalcCreator.h" #include "sbncode/SBNEventWeight/Base/WeightCalc.h" #include "larcore/Geometry/Geometry.h" +#include "geant4reweight/ReweightBase/G4ReweightManager.hh" #include "geant4reweight/ReweightBase/G4ReweighterFactory.hh" #include "geant4reweight/ReweightBase/G4Reweighter.hh" #include "geant4reweight/ReweightBase/G4ReweightTraj.hh" @@ -62,7 +63,8 @@ class Geant4WeightCalc : public WeightCalc { 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 - G4ReweighterFactory RWFactory; //!< Base class to handle all Geant4Reweighters (right now "all" means pi+, pi-, p) + 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; std::vector> 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 @@ -85,8 +87,12 @@ class Geant4WeightCalc : public WeightCalc { std::vector< double > p_energies_inel; //!< Variables for by-particle output tree std::vector< int > p_sliceInts_inel; //!< Variables for by-particle output tree int p_nElasticScatters; //!< Variables for by-particle output tree - std::vector p_inel_weight; //!< Variables for by-particle output tree - std::vector p_elastic_weight; //!< Variables for by-particle output + std::vector p_weight; //!< Variables for by-particle output tree + + + // **IF** I understand correctly, then after commit 11077f2 (el_weight*inel_weight) is returned by GetWeight() + //std::vector p_inel_weight; //!< Variables for by-particle output tree + //std::vector p_elastic_weight; //!< Variables for by-particle output bool fDebug; //!< print debug statements @@ -121,10 +127,9 @@ void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, TFile XSecFile( XSecFileName.c_str(), "OPEN" ); // Configure G4Reweighter - bool totalonly = false; - if (fPdg==2212) totalonly = true; - ParMaker = new G4ReweightParameterMaker( FitParSets, totalonly ); - theReweighter = RWFactory.BuildReweighter(fPdg, &XSecFile, &FracsFile, ParMaker->GetFSHists(), ParMaker->GetElasticHist() ); + 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() ); // Make output trees to save things for quick and easy validation art::ServiceHandle tfs; @@ -136,8 +141,10 @@ void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, fOutTree_Particle->Branch("subrun_num",&subrun_num); fOutTree_Particle->Branch("UniverseVals", &UniverseVals); fOutTree_Particle->Branch("pdg_to_reweight",&fPdg); - fOutTree_Particle->Branch("inelastic_weight",&p_inel_weight); - fOutTree_Particle->Branch("elastic_weight",&p_elastic_weight); + + fOutTree_Particle->Branch("weight",&p_weight); + //fOutTree_Particle->Branch("inelastic_weight",&p_inel_weight); + //fOutTree_Particle->Branch("elastic_weight",&p_elastic_weight); fOutTree_Particle->Branch("track_length",&p_track_length); fOutTree_Particle->Branch("PDG",&p_PDG); fOutTree_Particle->Branch("final_proc",&p_final_proc); @@ -252,10 +259,12 @@ std::vector Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { for (size_t i=0; i Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { // Loop through universes (j) for (size_t j=0; jSetNewVals(UniverseVals.at(j)); @@ -440,6 +449,10 @@ std::vector Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { // Total weight is the product of track weights in the event weight[j] *= std::max((float)0.0, w); + // just for the output tree + p_weight[j] = w; + +/* // Do the same for elastic weight (should be 1 unless set to non-nominal ) el_w = theReweighter->GetElasticWeight( &theTraj ); weight[j] *= std::max((float)0.0,el_w); @@ -453,12 +466,13 @@ std::vector Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { // std::cout << UniverseVals.at(j) << std::endl; std::cout << " w = " << w << ", el_w = " << el_w << std::endl; } - +*/ } // loop through universes (j) if (fDebug){ std::cout << "PDG = " << p_PDG << std::endl; +/* std::cout << "inel weights by particle: "; for (unsigned int j=0; j Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { for (unsigned int j=0; j Date: Thu, 25 Jul 2024 20:16:24 -0500 Subject: [PATCH 07/16] Add material input; This build ~works~ as long as you only run one particle type at time --- .../Calculators/Geant4/Geant4WeightCalc.cxx | 26 ++++++++++--------- .../jobs/geant4/eventweight_geant4_sbn.fcl | 20 +++++++++----- .../jobs/geant4/material_parameters.fcl | 18 +++++++++++++ 3 files changed, 45 insertions(+), 19 deletions(-) create mode 100644 sbncode/SBNEventWeight/jobs/geant4/material_parameters.fcl diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx index 401c08896..dd200f14c 100644 --- a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx +++ b/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 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> 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 ("number_of_multisims", 0); fPdg = pset.get ("pdg_to_reweight"); + fMaterial = pset.get ("material"); fDebug = pset.get ("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 tfs; @@ -441,11 +443,11 @@ std::vector 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 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 diff --git a/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl b/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl index e45065880..997b40157 100644 --- a/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl +++ b/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 # diff --git a/sbncode/SBNEventWeight/jobs/geant4/material_parameters.fcl b/sbncode/SBNEventWeight/jobs/geant4/material_parameters.fcl new file mode 100644 index 000000000..ed735621f --- /dev/null +++ b/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 From 5a02db64b0da60edc07af46ade92db2f82b9471a Mon Sep 17 00:00:00 2001 From: Jack Smedley Date: Fri, 26 Jul 2024 19:37:13 -0500 Subject: [PATCH 08/16] Iitialize G4ReweightManager as static so only one G4RunManager is started --- .../SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx | 5 ++--- .../SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl | 2 +- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx index dd200f14c..ebb5ad44d 100644 --- a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx +++ b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx @@ -64,7 +64,6 @@ class Geant4WeightCalc : public WeightCalc { 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 *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; @@ -130,8 +129,8 @@ void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, // Configure G4Reweighter 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() ); + static G4ReweightManager* RWManager = new G4ReweightManager( {fMaterial} ); //Constructor asks for a vector, but SBN only cares about one material + fReweighter = fRWFactory.BuildReweighter(fPdg, &FracsFile, fParMaker->GetFSHists(), fMaterial, RWManager, fParMaker->GetElasticHist() ); // Make output trees to save things for quick and easy validation art::ServiceHandle tfs; diff --git a/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl b/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl index 997b40157..aa565bf9e 100644 --- a/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl +++ b/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl @@ -24,7 +24,7 @@ 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 From afc823dc080b4e2eea40f995bfc72d39f629bebb Mon Sep 17 00:00:00 2001 From: Jack Smedley Date: Thu, 1 Aug 2024 15:23:14 -0500 Subject: [PATCH 09/16] Add neutron and kaon parameters, remove dependence on input XS file --- .../Calculators/Geant4/Geant4WeightCalc.cxx | 4 +- .../jobs/geant4/eventweight_geant4_sbn.fcl | 49 +++++++++++++++++-- .../geant4/kminus_reweight_parameters.fcl | 21 ++++++++ .../jobs/geant4/kplus_reweight_parameters.fcl | 21 ++++++++ .../geant4/neutron_reweight_parameters.fcl | 21 ++++++++ 5 files changed, 109 insertions(+), 7 deletions(-) create mode 100644 sbncode/SBNEventWeight/jobs/geant4/kminus_reweight_parameters.fcl create mode 100644 sbncode/SBNEventWeight/jobs/geant4/kplus_reweight_parameters.fcl create mode 100644 sbncode/SBNEventWeight/jobs/geant4/neutron_reweight_parameters.fcl diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx index ebb5ad44d..16d684f72 100644 --- a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx +++ b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx @@ -113,7 +113,6 @@ void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, fMakeOutputTrees = pset.get< bool >( "makeoutputtree", false ); std::string mode = pset.get("mode"); std::string FracsFileName = pset.get< std::string >( "fracsfile" ); - std::string XSecFileName = pset.get< std::string >( "xsecfile" ); std::vector< fhicl::ParameterSet > FitParSets = pset.get< std::vector< fhicl::ParameterSet > >("parameters"); fNsims = pset.get ("number_of_multisims", 0); fPdg = pset.get ("pdg_to_reweight"); @@ -123,9 +122,8 @@ void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, // Prepare random generator fGaussRandom = new CLHEP::RandGaussQ(engine); - // Get input files + // Get input file TFile FracsFile( FracsFileName.c_str(), "OPEN" ); - TFile XSecFile( XSecFileName.c_str(), "OPEN" ); // Configure G4Reweighter fParMaker = new G4ReweightParameterMaker( FitParSets, true , fPdg ); //TODO:Do we want check_overlap? Maybe a fcl switch? diff --git a/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl b/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl index aa565bf9e..6711f56d5 100644 --- a/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl +++ b/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl @@ -14,6 +14,10 @@ #include "piplus_reweight_parameters.fcl" #include "piminus_reweight_parameters.fcl" #include "proton_reweight_parameters.fcl" +#include "neutron_reweight_parameters.fcl" +#include "kplus_reweight_parameters.fcl" +#include "kminus_reweight_parameters.fcl" + #include "material_parameters.fcl" BEGIN_PROLOG @@ -24,7 +28,7 @@ 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_neutron , reinteractions_kplus , reinteractions_kminus ] reinteractions_piplus: { type: Geant4 @@ -34,7 +38,6 @@ sbn_eventweight_geant4: { 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: true pdg_to_reweight: 211 debug: true @@ -48,7 +51,6 @@ sbn_eventweight_geant4: { 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: true pdg_to_reweight: -211 debug: true @@ -62,11 +64,50 @@ sbn_eventweight_geant4: { 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: true pdg_to_reweight: 2212 debug: true } + + reinteractions_neutron: { + type: Geant4 + material: @local::LAr + random_seed: 60 + parameters: @local::NeutronParameters + mode: multisim + number_of_multisims: 1000 + fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_proton.root" + makeoutputtree: true + pdg_to_reweight: 2112 + debug: true + } + + reinteractions_kplus: { + type: Geant4 + material: @local::LAr + random_seed: 60 + parameters: @local::KPlusParameters + mode: multisim + number_of_multisims: 1000 + fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_proton.root" + makeoutputtree: true + pdg_to_reweight: 321 + debug: true + } + + reinteractions_kminus: { + type: Geant4 + material: @local::LAr + random_seed: 60 + parameters: @local::KMinusParameters + mode: multisim + number_of_multisims: 1000 + fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_proton.root" + makeoutputtree: true + pdg_to_reweight: -321 + debug: true + } + } # reinteractions_neutron diff --git a/sbncode/SBNEventWeight/jobs/geant4/kminus_reweight_parameters.fcl b/sbncode/SBNEventWeight/jobs/geant4/kminus_reweight_parameters.fcl new file mode 100644 index 000000000..3713f0183 --- /dev/null +++ b/sbncode/SBNEventWeight/jobs/geant4/kminus_reweight_parameters.fcl @@ -0,0 +1,21 @@ +BEGIN_PROLOG + +KMinusParameters: [ + { + Name: "fKMinusReac" + Cut: "reac" + Range: [10., 2005.] + Nominal: 1 + Sigma: 0.2 + }, + + { + Name: "fKMinusElast" + Cut: "elast" + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + } +] + +END_PROLOG diff --git a/sbncode/SBNEventWeight/jobs/geant4/kplus_reweight_parameters.fcl b/sbncode/SBNEventWeight/jobs/geant4/kplus_reweight_parameters.fcl new file mode 100644 index 000000000..b426eb6bd --- /dev/null +++ b/sbncode/SBNEventWeight/jobs/geant4/kplus_reweight_parameters.fcl @@ -0,0 +1,21 @@ +BEGIN_PROLOG + +KPlusParameters: [ + { + Name: "fKPlusReac" + Cut: "reac" + Range: [10., 2005.] + Nominal: 1 + Sigma: 0.2 + }, + + { + Name: "fKPlusElast" + Cut: "elast" + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + } +] + +END_PROLOG diff --git a/sbncode/SBNEventWeight/jobs/geant4/neutron_reweight_parameters.fcl b/sbncode/SBNEventWeight/jobs/geant4/neutron_reweight_parameters.fcl new file mode 100644 index 000000000..a5c462d84 --- /dev/null +++ b/sbncode/SBNEventWeight/jobs/geant4/neutron_reweight_parameters.fcl @@ -0,0 +1,21 @@ +BEGIN_PROLOG + +NeutronParameters: [ + { + Name: "fNeutronReac" + Cut: "reac" + Range: [10., 2005.] + Nominal: 1 + Sigma: 0.2 + }, + + { + Name: "fNeutronElast" + Cut: "elast" + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + } +] + +END_PROLOG From 6ab4aff6f7fc718d13f558532d3e5467afe42be6 Mon Sep 17 00:00:00 2001 From: Patrick Green Date: Thu, 24 Jul 2025 07:21:57 -0500 Subject: [PATCH 10/16] fixing dependencies --- ups/product_deps | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/ups/product_deps b/ups/product_deps index 9dfd1799f..467fe7bf5 100644 --- a/ups/product_deps +++ b/ups/product_deps @@ -248,19 +248,19 @@ libdir fq_dir lib # cetmodules), an entry for cetmodules is required. # # * It is an error for more than one non-( == "-default-") - entry to match for a given product. +# entry to match for a given product. # #################################### product version qual flags genie_xsec v3_04_00 - larcv2 v2_2_6 - -larsoft v10_06_00 - -sbnalg v10_06_00_01 - +larsoft v10_06_00_01 - +sbnalg v10_06_00_03 - sbndaq_artdaq_core v1_10_06 - sbndata v01_07 - systematicstools v01_04_04 - nusystematics v1_05_07 - -geant4reweight v01_20_08 - +geant4reweight v01_20_11 - cetmodules v3_24_01 - only_for_build end_product_list #################################### @@ -316,12 +316,11 @@ end_product_list # case it is optional. # #################################### - -qualifier larsoft sbnalg sbndaq_artdaq_core genie_xsec sbndata larcv2 systematicstools nusystematics notes -c14:debug c14:debug c14:debug c14:s131:debug AR2320i00000:e1000:k250 -nq- c14:debug:p3915 c14:debug c14:debug -nq- -c14:prof c14:prof c14:prof c14:s131:prof AR2320i00000:e1000:k250 -nq- c14:p3915:prof c14:prof c14:prof -nq- -e26:debug e26:debug e26:debug e26:s131:debug AR2320i00000:e1000:k250 -nq- debug:e26:p3915 e26:debug e26:debug -nq- -e26:prof e26:prof e26:prof e26:s131:prof AR2320i00000:e1000:k250 -nq- e26:p3915:prof e26:prof e26:prof -nq- +qualifier larsoft sbnalg sbndaq_artdaq_core genie_xsec sbndata larcv2 systematicstools nusystematics geant4reweight notes +c14:debug c14:debug c14:debug c14:s131:debug AR2320i00000:e1000:k250 -nq- c14:debug:p3915 c14:debug c14:debug c14:debug:s131 -nq- +c14:prof c14:prof c14:prof c14:s131:prof AR2320i00000:e1000:k250 -nq- c14:p3915:prof c14:prof c14:prof c14:prof:s131 -nq- +e26:debug e26:debug e26:debug e26:s131:debug AR2320i00000:e1000:k250 -nq- debug:e26:p3915 e26:debug e26:debug e26:debug:s131 -nq- +e26:prof e26:prof e26:prof e26:s131:prof AR2320i00000:e1000:k250 -nq- e26:p3915:prof e26:prof e26:prof e26:prof:s131 -nq- end_qualifier_list #################################### From e506332515b8bcb452ee53d1d1fce0e4aa38abc6 Mon Sep 17 00:00:00 2001 From: Jack Smedley Date: Tue, 22 Jul 2025 14:06:10 -0500 Subject: [PATCH 11/16] Testing the switch from multisim to multisigma --- .../Calculators/Geant4/Geant4WeightCalc.cxx | 19 ++++++++++++++++++- .../jobs/geant4/eventweight_geant4_sbn.fcl | 14 ++++++++++---- .../geant4/neutron_reweight_parameters.fcl | 9 +++++++++ .../geant4/proton_reweight_parameters.fcl | 8 ++++++++ 4 files changed, 45 insertions(+), 5 deletions(-) diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx index 16d684f72..2a551eab6 100644 --- a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx +++ b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx @@ -61,7 +61,8 @@ class Geant4WeightCalc : public WeightCalc { CLHEP::RandGaussQ* fGaussRandom; //!< Random number generator fhicl::ParameterSet fMaterial; //!< Detector material, i.e. LAr // std::map fParticles; //!< Particles to reweight - unsigned fNsims; //!< Number of multisims + unsigned fNsims; //!< Number of universes (multisim mode) + unsigned fNsigmas; //!< Number of sigmas (multisigma mode) 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 G4ReweighterFactory fRWFactory; //!< Base class to handle all Geant4Reweighters @@ -115,6 +116,7 @@ void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, std::string FracsFileName = pset.get< std::string >( "fracsfile" ); std::vector< fhicl::ParameterSet > FitParSets = pset.get< std::vector< fhicl::ParameterSet > >("parameters"); fNsims = pset.get ("number_of_multisims", 0); + fNsigmas = pset.get ("number_of_sigmas", 0); fPdg = pset.get ("pdg_to_reweight"); fMaterial = pset.get ("material"); fDebug = pset.get ("debug",false); @@ -171,6 +173,7 @@ void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, std::vector FitParSigmas; std::map theNominals; + std::cout << "Configuring parameter: " << GetFullName() << " with mode: " << mode << std::endl; fParameterSet.Configure(GetFullName(), mode, fNsims); for (size_t i_parset=0; i_parset tmp_vals(theNominals); + tmp_vals[FitParNames.at(i_parset)] = FitParNominals.at(i_parset) + i_sigma*FitParSigmas.at(i_parset); + UniverseVals.push_back(tmp_vals); + } // ened loop over sigmas + } // end loop over parsets (i) + } // multisigma + else{ // Anything else mode: Set parameters to user-defined nominal value UniverseVals.push_back(theNominals); diff --git a/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl b/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl index 6711f56d5..f4eeab56f 100644 --- a/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl +++ b/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl @@ -61,8 +61,10 @@ sbn_eventweight_geant4: { material: @local::LAr random_seed: 60 parameters: @local::ProtonParameters - mode: multisim - number_of_multisims: 1000 + mode: multisigma + number_of_sigmas: 3 + #mode: multisim + #number_of_multisims: 1000 fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_proton.root" makeoutputtree: true pdg_to_reweight: 2212 @@ -74,8 +76,10 @@ sbn_eventweight_geant4: { material: @local::LAr random_seed: 60 parameters: @local::NeutronParameters - mode: multisim - number_of_multisims: 1000 + mode: multisigma + number_of_sigmas: 3 + #mode: multisim + #number_of_multisims: 1000 fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_proton.root" makeoutputtree: true pdg_to_reweight: 2112 @@ -87,6 +91,7 @@ sbn_eventweight_geant4: { material: @local::LAr random_seed: 60 parameters: @local::KPlusParameters + #mode: pm1sigma mode: multisim number_of_multisims: 1000 fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_proton.root" @@ -100,6 +105,7 @@ sbn_eventweight_geant4: { material: @local::LAr random_seed: 60 parameters: @local::KMinusParameters + #mode: pm1sigma mode: multisim number_of_multisims: 1000 fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_proton.root" diff --git a/sbncode/SBNEventWeight/jobs/geant4/neutron_reweight_parameters.fcl b/sbncode/SBNEventWeight/jobs/geant4/neutron_reweight_parameters.fcl index a5c462d84..b0a3e782a 100644 --- a/sbncode/SBNEventWeight/jobs/geant4/neutron_reweight_parameters.fcl +++ b/sbncode/SBNEventWeight/jobs/geant4/neutron_reweight_parameters.fcl @@ -15,7 +15,16 @@ NeutronParameters: [ Range: [10.0, 2005.00] Nominal: 1.0 Sigma: 0.2 + }, + + { + Name: "fNeutronTotal" + Cut: "total" + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 } + ] END_PROLOG diff --git a/sbncode/SBNEventWeight/jobs/geant4/proton_reweight_parameters.fcl b/sbncode/SBNEventWeight/jobs/geant4/proton_reweight_parameters.fcl index 13113db75..8bdc9435b 100644 --- a/sbncode/SBNEventWeight/jobs/geant4/proton_reweight_parameters.fcl +++ b/sbncode/SBNEventWeight/jobs/geant4/proton_reweight_parameters.fcl @@ -17,6 +17,14 @@ ProtonParameters: [ Range: [10.0, 2005.00] Nominal: 1.0 Sigma: 0.2 + }, + + { + Name: "fProtonTotal" + Cut: "total" + Range: [10., 2005.] + Nominal: 1 + Sigma: 0.2 } ] From c486b2943a98dfa8ab3961835deea252aaa5267c Mon Sep 17 00:00:00 2001 From: Patrick Green Date: Wed, 30 Jul 2025 09:36:46 -0500 Subject: [PATCH 12/16] clean up configurations to only keep cases that have implementations in geant4reweight --- .../jobs/geant4/eventweight_geant4_sbn.fcl | 38 +++++++++---------- .../geant4/kminus_reweight_parameters.fcl | 16 +++----- .../jobs/geant4/kplus_reweight_parameters.fcl | 16 +++----- .../geant4/neutron_reweight_parameters.fcl | 21 ++-------- .../geant4/piminus_reweight_parameters.fcl | 11 +----- .../geant4/piplus_reweight_parameters.fcl | 13 ++----- .../geant4/proton_reweight_parameters.fcl | 22 ++--------- 7 files changed, 39 insertions(+), 98 deletions(-) diff --git a/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl b/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl index f4eeab56f..af11d7ea5 100644 --- a/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl +++ b/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl @@ -61,11 +61,11 @@ sbn_eventweight_geant4: { material: @local::LAr random_seed: 60 parameters: @local::ProtonParameters - mode: multisigma - number_of_sigmas: 3 - #mode: multisim - #number_of_multisims: 1000 - fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_proton.root" + #mode: multisigma + #number_of_sigmas: 3 + mode: multisim + number_of_multisims: 1000 + fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_proton.root" # dummy file, cross section pulled directly from Geant4 makeoutputtree: true pdg_to_reweight: 2212 debug: true @@ -74,13 +74,13 @@ sbn_eventweight_geant4: { reinteractions_neutron: { type: Geant4 material: @local::LAr - random_seed: 60 + random_seed: 61 parameters: @local::NeutronParameters - mode: multisigma - number_of_sigmas: 3 - #mode: multisim - #number_of_multisims: 1000 - fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_proton.root" + #mode: multisigma + #number_of_sigmas: 3 + mode: multisim + number_of_multisims: 1000 + fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_proton.root" # dummy file, cross section pulled directly from Geant4 makeoutputtree: true pdg_to_reweight: 2112 debug: true @@ -89,12 +89,12 @@ sbn_eventweight_geant4: { reinteractions_kplus: { type: Geant4 material: @local::LAr - random_seed: 60 + random_seed: 62 parameters: @local::KPlusParameters #mode: pm1sigma mode: multisim number_of_multisims: 1000 - fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_proton.root" + fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_proton.root" # dummy file, cross section pulled directly from Geant4 makeoutputtree: true pdg_to_reweight: 321 debug: true @@ -103,12 +103,12 @@ sbn_eventweight_geant4: { reinteractions_kminus: { type: Geant4 material: @local::LAr - random_seed: 60 + random_seed: 63 parameters: @local::KMinusParameters #mode: pm1sigma mode: multisim number_of_multisims: 1000 - fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_proton.root" + fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_proton.root" # dummy file, cross section pulled directly from Geant4 makeoutputtree: true pdg_to_reweight: -321 debug: true @@ -116,8 +116,6 @@ sbn_eventweight_geant4: { } -# reinteractions_neutron - ### # Reweighting parameters should be defined as # @@ -139,15 +137,15 @@ sbn_eventweight_geant4: { # "nominal" value in that case, just the set value) # - Sigma defines the 1-sigma range for multisims and unisims. Currently the # code can only accept a single value for sigma, giving symmetric -# uncertainty bounds: nominal-sigma and nominal+sigma +# uncertainty bounds: nominal-sigma and nominal+sigma +# [Only in use for fitting purposes, instead range defined through input files - Jul 2025] # - Cut must be one of the following: -# "reac" <- total inelastic scattering cross section +# "total" <- total scattering cross section, pulled from Geant4 # "abs" <- absorption cross section (pi+ and pi- only at the moment -- Oct 2019) # "cex" <- charge exchange cross section (pi+ and pi- only at the moment -- Oct 2019) # "dcex" <- double charge exchange cross section (pi+ and pi- only at the moment -- Oct 2019) # "prod" <- pion production cross section (pi+ and pi- only at the moment -- Oct 2019) # "inel" <- quasi-elastic inelastic scattering cross section (pi+ and pi- only at the moment -- Oct 2019) -# "elast" <- total elastic scattering cross section # - Name can be anything you like (but should uniquely identify that parameter/ # variation) ### diff --git a/sbncode/SBNEventWeight/jobs/geant4/kminus_reweight_parameters.fcl b/sbncode/SBNEventWeight/jobs/geant4/kminus_reweight_parameters.fcl index 3713f0183..ac330892e 100644 --- a/sbncode/SBNEventWeight/jobs/geant4/kminus_reweight_parameters.fcl +++ b/sbncode/SBNEventWeight/jobs/geant4/kminus_reweight_parameters.fcl @@ -1,18 +1,12 @@ BEGIN_PROLOG -KMinusParameters: [ - { - Name: "fKMinusReac" - Cut: "reac" - Range: [10., 2005.] - Nominal: 1 - Sigma: 0.2 - }, +## KMinus available process keys: "total" +KMinusParameters: [ { - Name: "fKMinusElast" - Cut: "elast" - Range: [10.0, 2005.00] + Name: "fKMinusTotal" + Cut: "total" + Range: [10.0, 2005.0] Nominal: 1.0 Sigma: 0.2 } diff --git a/sbncode/SBNEventWeight/jobs/geant4/kplus_reweight_parameters.fcl b/sbncode/SBNEventWeight/jobs/geant4/kplus_reweight_parameters.fcl index b426eb6bd..08ec5e9f2 100644 --- a/sbncode/SBNEventWeight/jobs/geant4/kplus_reweight_parameters.fcl +++ b/sbncode/SBNEventWeight/jobs/geant4/kplus_reweight_parameters.fcl @@ -1,18 +1,12 @@ BEGIN_PROLOG -KPlusParameters: [ - { - Name: "fKPlusReac" - Cut: "reac" - Range: [10., 2005.] - Nominal: 1 - Sigma: 0.2 - }, +## KPlus available process keys: "total" +KPlusParameters: [ { - Name: "fKPlusElast" - Cut: "elast" - Range: [10.0, 2005.00] + Name: "fKPlusTotal" + Cut: "total" + Range: [10.0, 2005.0] Nominal: 1.0 Sigma: 0.2 } diff --git a/sbncode/SBNEventWeight/jobs/geant4/neutron_reweight_parameters.fcl b/sbncode/SBNEventWeight/jobs/geant4/neutron_reweight_parameters.fcl index b0a3e782a..f33d5e030 100644 --- a/sbncode/SBNEventWeight/jobs/geant4/neutron_reweight_parameters.fcl +++ b/sbncode/SBNEventWeight/jobs/geant4/neutron_reweight_parameters.fcl @@ -1,30 +1,15 @@ BEGIN_PROLOG -NeutronParameters: [ - { - Name: "fNeutronReac" - Cut: "reac" - Range: [10., 2005.] - Nominal: 1 - Sigma: 0.2 - }, - - { - Name: "fNeutronElast" - Cut: "elast" - Range: [10.0, 2005.00] - Nominal: 1.0 - Sigma: 0.2 - }, +## Neutron available process keys: "total" +NeutronParameters: [ { Name: "fNeutronTotal" Cut: "total" - Range: [10.0, 2005.00] + Range: [10.0, 2005.0] Nominal: 1.0 Sigma: 0.2 } - ] END_PROLOG diff --git a/sbncode/SBNEventWeight/jobs/geant4/piminus_reweight_parameters.fcl b/sbncode/SBNEventWeight/jobs/geant4/piminus_reweight_parameters.fcl index b64756415..3bf335877 100644 --- a/sbncode/SBNEventWeight/jobs/geant4/piminus_reweight_parameters.fcl +++ b/sbncode/SBNEventWeight/jobs/geant4/piminus_reweight_parameters.fcl @@ -1,6 +1,6 @@ BEGIN_PROLOG -# From MicroBooNE, ubsim UBOONE_SUITE_v08_00_00_69 +## PiMinus available process keys: "abs", "cex", "prod", "dcex", "inel" PiMinusParameters: [ { @@ -52,15 +52,8 @@ PiMinusParameters: [ Range: [10.0, 2005.00] Nominal: 1.0 Sigma: 0.2 - }, - - { - Name: "fPiMinusElast" - Cut: "elast" - Range: [10.0, 2005.00] - Nominal: 1.0 - Sigma: 0.2 } + ] END_PROLOG diff --git a/sbncode/SBNEventWeight/jobs/geant4/piplus_reweight_parameters.fcl b/sbncode/SBNEventWeight/jobs/geant4/piplus_reweight_parameters.fcl index c9cfc3e23..9524a28fd 100644 --- a/sbncode/SBNEventWeight/jobs/geant4/piplus_reweight_parameters.fcl +++ b/sbncode/SBNEventWeight/jobs/geant4/piplus_reweight_parameters.fcl @@ -1,13 +1,13 @@ BEGIN_PROLOG -# From MicroBooNE, ubsim UBOONE_SUITE_v08_00_00_69 +## PiPlus available process keys: "abs", "cex", "prod", "dcex", "inel" PiPlusParameters: [ { Name: "fPiPlusReacLow" Cut: "reac" Range: [10., 200.] - Nominal: 1 + Nominal: 1.0 Sigma: 0.2 } , { @@ -52,15 +52,8 @@ PiPlusParameters: [ Range: [10.0, 2005.00] Nominal: 1.0 Sigma: 0.2 - }, - - { - Name: "fPiPlusElast" - Cut: "elast" - Range: [10.0, 2005.00] - Nominal: 1.0 - Sigma: 0.2 } + ] END_PROLOG diff --git a/sbncode/SBNEventWeight/jobs/geant4/proton_reweight_parameters.fcl b/sbncode/SBNEventWeight/jobs/geant4/proton_reweight_parameters.fcl index 8bdc9435b..3dd017d93 100644 --- a/sbncode/SBNEventWeight/jobs/geant4/proton_reweight_parameters.fcl +++ b/sbncode/SBNEventWeight/jobs/geant4/proton_reweight_parameters.fcl @@ -1,29 +1,13 @@ BEGIN_PROLOG -# From MicroBooNE, ubsim UBOONE_SUITE_v08_00_00_69 +## Proton available process keys: "total" ProtonParameters: [ - { - Name: "fProtonReac" - Cut: "reac" - Range: [10., 2005.] - Nominal: 1 - Sigma: 0.2 - }, - - { - Name: "fProtonElast" - Cut: "elast" - Range: [10.0, 2005.00] - Nominal: 1.0 - Sigma: 0.2 - }, - { Name: "fProtonTotal" Cut: "total" - Range: [10., 2005.] - Nominal: 1 + Range: [10.0, 2005.0] + Nominal: 1.0 Sigma: 0.2 } ] From 70bf4195bbcc5ea482ddca83fc78dff24c84d08b Mon Sep 17 00:00:00 2001 From: Patrick Green Date: Wed, 30 Jul 2025 10:35:11 -0500 Subject: [PATCH 13/16] adding missing particle masses --- .../SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx index 2a551eab6..56afd66a8 100644 --- a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx +++ b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx @@ -300,7 +300,11 @@ std::vector Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { double mass = 0.; if( TMath::Abs(p_PDG) == 211 ) mass = 139.57; + else if ( p_PDG == 111 ) mass = 134.9768; else if( p_PDG == 2212 ) mass = 938.28; + else if ( p_PDG == 2112 ) mass = 939.565; + else if ( TMath::Abs(p_PDG) == 321) mass = 493.677; + else if ( p_PDG == 311) mass = 497.61; // We only want to record weights for one type of particle (defined by fPDG from the fcl file), so skip other particles if (p_PDG == fPdg){ From 081b5ba89ffea13de4d2a2055b7027b736e92d16 Mon Sep 17 00:00:00 2001 From: Patrick Green Date: Mon, 4 Aug 2025 15:12:38 -0500 Subject: [PATCH 14/16] adding neutron threshold and cleaning up configurations --- .../Calculators/Geant4/Geant4WeightCalc.cxx | 5 +++-- .../jobs/geant4/eventweight_geant4_sbn.fcl | 1 - .../jobs/geant4/neutron_reweight_parameters.fcl | 2 +- .../jobs/geant4/piminus_reweight_parameters.fcl | 12 ++++-------- .../jobs/geant4/piplus_reweight_parameters.fcl | 12 ++++-------- 5 files changed, 12 insertions(+), 20 deletions(-) diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx index 56afd66a8..1889e9c90 100644 --- a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx +++ b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx @@ -186,6 +186,8 @@ void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, FitParNominals.push_back(theNominal); FitParSigmas.push_back(theSigma); + if (fDebug) std::cout << "Name: " << theName << ", Nominal: " << theNominal << ", Sigma: " << theSigma << std::endl; + theNominals[theName] = theNominal; fParameterSet.AddParameter(theName, theSigma); @@ -514,8 +516,7 @@ std::vector Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { } std::cout << std::endl; } - - } // if ( ( TMath::Abs(p_PDG) == 211 || p_PDG == 2212 ) ) + } if (fMakeOutputTrees) fOutTree_Particle->Fill(); } // loop over mcparticles (i) if (fMakeOutputTrees) fOutTree_MCTruth->Fill(); diff --git a/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl b/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl index af11d7ea5..081964d3f 100644 --- a/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl +++ b/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl @@ -138,7 +138,6 @@ sbn_eventweight_geant4: { # - Sigma defines the 1-sigma range for multisims and unisims. Currently the # code can only accept a single value for sigma, giving symmetric # uncertainty bounds: nominal-sigma and nominal+sigma -# [Only in use for fitting purposes, instead range defined through input files - Jul 2025] # - Cut must be one of the following: # "total" <- total scattering cross section, pulled from Geant4 # "abs" <- absorption cross section (pi+ and pi- only at the moment -- Oct 2019) diff --git a/sbncode/SBNEventWeight/jobs/geant4/neutron_reweight_parameters.fcl b/sbncode/SBNEventWeight/jobs/geant4/neutron_reweight_parameters.fcl index f33d5e030..65743db98 100644 --- a/sbncode/SBNEventWeight/jobs/geant4/neutron_reweight_parameters.fcl +++ b/sbncode/SBNEventWeight/jobs/geant4/neutron_reweight_parameters.fcl @@ -6,7 +6,7 @@ NeutronParameters: [ { Name: "fNeutronTotal" Cut: "total" - Range: [10.0, 2005.0] + Range: [100.0, 2005.0] Nominal: 1.0 Sigma: 0.2 } diff --git a/sbncode/SBNEventWeight/jobs/geant4/piminus_reweight_parameters.fcl b/sbncode/SBNEventWeight/jobs/geant4/piminus_reweight_parameters.fcl index 3bf335877..fb27aeed4 100644 --- a/sbncode/SBNEventWeight/jobs/geant4/piminus_reweight_parameters.fcl +++ b/sbncode/SBNEventWeight/jobs/geant4/piminus_reweight_parameters.fcl @@ -4,15 +4,15 @@ BEGIN_PROLOG PiMinusParameters: [ { - Name: "fPiMinusReacLow" - Cut: "reac" + Name: "fPiMinusInelLow" + Cut: "inel" Range: [10., 200.] Nominal: 1 Sigma: 0.2 } , { - Name: "fPiMinusReacHigh" - Cut: "reac" + Name: "fPiMinusInelHigh" + Cut: "inel" Range: [700., 2005.] Nominal: 1.0 Sigma: 0.2 @@ -21,7 +21,6 @@ PiMinusParameters: [ { Name: "fPiMinusAbs" Cut: "abs" - # Range: [200.0, 700.00] Range: [10.0, 2005.00] Nominal: 1.0 Sigma: 0.2 @@ -30,7 +29,6 @@ PiMinusParameters: [ { Name: "fPiMinusCex" Cut: "cex" - # Range: [200.0, 700.00] Range: [10.0, 2005.00] Nominal: 1.0 Sigma: 0.2 @@ -39,7 +37,6 @@ PiMinusParameters: [ { Name: "fPiMinusDCex" Cut: "dcex" - # Range: [200.0, 700.00] Range: [10.0, 2005.00] Nominal: 1.0 Sigma: 0.2 @@ -48,7 +45,6 @@ PiMinusParameters: [ { Name: "fPiMinusPiProd" Cut: "prod" - # Range: [200.0, 700.00] Range: [10.0, 2005.00] Nominal: 1.0 Sigma: 0.2 diff --git a/sbncode/SBNEventWeight/jobs/geant4/piplus_reweight_parameters.fcl b/sbncode/SBNEventWeight/jobs/geant4/piplus_reweight_parameters.fcl index 9524a28fd..4aeaeb96d 100644 --- a/sbncode/SBNEventWeight/jobs/geant4/piplus_reweight_parameters.fcl +++ b/sbncode/SBNEventWeight/jobs/geant4/piplus_reweight_parameters.fcl @@ -4,15 +4,15 @@ BEGIN_PROLOG PiPlusParameters: [ { - Name: "fPiPlusReacLow" - Cut: "reac" + Name: "fPiPlusInelLow" + Cut: "inel" Range: [10., 200.] Nominal: 1.0 Sigma: 0.2 } , { - Name: "fPiPlusReacHigh" - Cut: "reac" + Name: "fPiPlusInelHigh" + Cut: "inel" Range: [700., 2005.] Nominal: 1.0 Sigma: 0.2 @@ -21,7 +21,6 @@ PiPlusParameters: [ { Name: "fPiPlusAbs" Cut: "abs" - # Range: [200.0, 700.00] Range: [10.0, 2005.00] Nominal: 1.0 Sigma: 0.2 @@ -30,7 +29,6 @@ PiPlusParameters: [ { Name: "fPiPlusCex" Cut: "cex" - # Range: [200.0, 700.00] Range: [10.0, 2005.00] Nominal: 1.0 Sigma: 0.2 @@ -39,7 +37,6 @@ PiPlusParameters: [ { Name: "fPiPlusDCex" Cut: "dcex" - # Range: [200.0, 700.00] Range: [10.0, 2005.00] Nominal: 1.0 Sigma: 0.2 @@ -48,7 +45,6 @@ PiPlusParameters: [ { Name: "fPiPlusPiProd" Cut: "prod" - # Range: [200.0, 700.00] Range: [10.0, 2005.00] Nominal: 1.0 Sigma: 0.2 From 83479a0eca91984f1e89faccd7cf71bd15a90b6e Mon Sep 17 00:00:00 2001 From: Patrick Green Date: Mon, 4 Aug 2025 15:35:22 -0500 Subject: [PATCH 15/16] disable debug mode by default --- .../jobs/geant4/eventweight_geant4_sbn.fcl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl b/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl index 081964d3f..2431af89d 100644 --- a/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl +++ b/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl @@ -40,7 +40,7 @@ sbn_eventweight_geant4: { fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_piplus.root" makeoutputtree: true pdg_to_reweight: 211 - debug: true + debug: false } reinteractions_piminus: { @@ -53,7 +53,7 @@ sbn_eventweight_geant4: { fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_piminus.root" makeoutputtree: true pdg_to_reweight: -211 - debug: true + debug: false } reinteractions_proton: { @@ -68,7 +68,7 @@ sbn_eventweight_geant4: { fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_proton.root" # dummy file, cross section pulled directly from Geant4 makeoutputtree: true pdg_to_reweight: 2212 - debug: true + debug: false } reinteractions_neutron: { @@ -83,7 +83,7 @@ sbn_eventweight_geant4: { fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_proton.root" # dummy file, cross section pulled directly from Geant4 makeoutputtree: true pdg_to_reweight: 2112 - debug: true + debug: false } reinteractions_kplus: { @@ -97,7 +97,7 @@ sbn_eventweight_geant4: { fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_proton.root" # dummy file, cross section pulled directly from Geant4 makeoutputtree: true pdg_to_reweight: 321 - debug: true + debug: false } reinteractions_kminus: { @@ -111,7 +111,7 @@ sbn_eventweight_geant4: { fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_proton.root" # dummy file, cross section pulled directly from Geant4 makeoutputtree: true pdg_to_reweight: -321 - debug: true + debug: false } } From 6051a113bdf70cc2e82900fdeb25d39c1af75b44 Mon Sep 17 00:00:00 2001 From: Patrick Green Date: Sun, 17 Aug 2025 22:24:16 -0500 Subject: [PATCH 16/16] cleaning up code --- .../Calculators/Geant4/Geant4WeightCalc.cxx | 72 +------------------ 1 file changed, 3 insertions(+), 69 deletions(-) diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx index 1889e9c90..ec9d8cc6e 100644 --- a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx +++ b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx @@ -60,11 +60,9 @@ class Geant4WeightCalc : public WeightCalc { std::string fMCTruthProducer; //!< Label for the MCTruth producer CLHEP::RandGaussQ* fGaussRandom; //!< Random number generator fhicl::ParameterSet fMaterial; //!< Detector material, i.e. LAr - // std::map fParticles; //!< Particles to reweight unsigned fNsims; //!< Number of universes (multisim mode) unsigned fNsigmas; //!< Number of sigmas (multisigma mode) 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 G4ReweighterFactory fRWFactory; //!< Base class to handle all Geant4Reweighters G4Reweighter *fReweighter; //!< Geant4Reweighter -- this is what provides the weights G4ReweightParameterMaker *fParMaker; @@ -90,10 +88,7 @@ class Geant4WeightCalc : public WeightCalc { int p_nElasticScatters; //!< Variables for by-particle output tree std::vector p_weight; //!< Variables for by-particle output tree - // **IF** I understand correctly, then after commit 11077f2 (el_weight*inel_weight) is returned by GetWeight() - //std::vector p_inel_weight; //!< Variables for by-particle output tree - //std::vector p_elastic_weight; //!< Variables for by-particle output bool fDebug; //!< print debug statements @@ -108,7 +103,6 @@ void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, // Get configuration for this function fhicl::ParameterSet const& pset = p.get(GetName()); - // std::cout << pset.GetName() << std::endl; fMCParticleProducer = pset.get("MCParticleProducer", "largeant"); fMCTruthProducer = pset.get("MCTruthProducer", "generator"); fMakeOutputTrees = pset.get< bool >( "makeoutputtree", false ); @@ -144,8 +138,6 @@ void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, fOutTree_Particle->Branch("pdg_to_reweight",&fPdg); fOutTree_Particle->Branch("weight",&p_weight); - //fOutTree_Particle->Branch("inelastic_weight",&p_inel_weight); - //fOutTree_Particle->Branch("elastic_weight",&p_elastic_weight); fOutTree_Particle->Branch("track_length",&p_track_length); fOutTree_Particle->Branch("PDG",&p_PDG); fOutTree_Particle->Branch("final_proc",&p_final_proc); @@ -165,7 +157,6 @@ void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, fOutTree_MCTruth->Branch("pdg_to_reweight",&fPdg); } - // Read input parameter sets and set up universes size_t n_parsets = FitParSets.size(); std::vector FitParNames; @@ -193,29 +184,7 @@ void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, fParameterSet.AddParameter(theName, theSigma); } - //DEPRECATED - if (mode=="pm1sigma"){ - // pm1sigma mode: 0 = +1sigma, 1 = -1sigma of a single parameter. All other parameters at nominal - for (size_t i_parset=0; i_parset tmp_vals_p1sigma(theNominals); - std::map tmp_vals_m1sigma(theNominals); - // Now reset the +1sigma and -1sigma values for this parameter set only - tmp_vals_p1sigma[FitParNames.at(i_parset)] = FitParNominals.at(i_parset)+FitParSigmas.at(i_parset); - tmp_vals_m1sigma[FitParNames.at(i_parset)] = FitParNominals.at(i_parset)-FitParSigmas.at(i_parset); - - if (fDebug){ - std::cout << "Universe " << i_parset*2 << ": " << FitParNames.at(i_parset) << " = " << FitParNominals.at(i_parset)+FitParSigmas.at(i_parset) << std::endl; - std::cout << "Universe " << i_parset*2+1 << ": " << FitParNames.at(i_parset) << " = " << FitParNominals.at(i_parset)-FitParSigmas.at(i_parset) << std::endl; - } - - // Finally, add these universes into the vector - UniverseVals.push_back(tmp_vals_p1sigma); - UniverseVals.push_back(tmp_vals_m1sigma); - } // end loop over parsets (i) - } // pm1sigma - - else if (mode=="multisim"){ + if (mode=="multisim"){ // multisim mode: parameter values sample within the given uncertainty for all parameters simultaneously // Loop over universes j for (unsigned j=0; j Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { // Reset things to be saved in the output tree for fast validation p_weight.clear(); p_weight.resize(fNsims,1.0); - //p_inel_weight.clear(); - //p_inel_weight.resize(fNsims,1.0); - //p_elastic_weight.clear(); - //p_elastic_weight.resize(fNsims,1.0); p_track_length=-9999; p_PDG=-9999; p_final_proc="dummy"; @@ -384,9 +349,6 @@ std::vector Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { p_nElasticScatters = elastic_indices.size(); for( size_t istep = 1; istep < nSteps; ++istep ){ - // if( istep == trajpoint_PX.size() - 1 && std::find( elastic_indices.begin(), elastic_indices.end(), j ) != elastic_indices.end() ) - // std::cout << "Warning: last step an elastic process" << std::endl; - std::string proc = "default"; if( istep == trajpoint_PX.size() - 1 ) proc = EndProcess; @@ -394,7 +356,6 @@ std::vector Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { proc = "hadElastic"; } - double deltaX = ( trajpoint_X.at(istep) - trajpoint_X.at(istep-1) ); double deltaY = ( trajpoint_Y.at(istep) - trajpoint_Y.at(istep-1) ); double deltaZ = ( trajpoint_Z.at(istep) - trajpoint_Z.at(istep-1) ); @@ -460,7 +421,7 @@ std::vector Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { // Loop through universes (j) for (size_t j=0; jSetNewVals(UniverseVals.at(j)); @@ -474,37 +435,10 @@ std::vector Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { // just for the output tree p_weight[j] = w; -/* - // Do the same for elastic weight (should be 1 unless set to non-nominal ) - el_w = fReweighter->GetElasticWeight( &theTraj ); - weight[j] *= std::max((float)0.0,el_w); - - // just for the output tree - p_inel_weight[j] = w; - p_elastic_weight[j] = el_w; - - if (fDebug){ - std::cout << " Universe " << j << ": "; - // std::cout << UniverseVals.at(j) << std::endl; - std::cout << " w = " << w << ", el_w = " << el_w << std::endl; - } -*/ - } // loop through universes (j) if (fDebug){ std::cout << "PDG = " << p_PDG << std::endl; -/* - std::cout << "inel weights by particle: "; - for (unsigned int j=0; j