From 3c7e3e8b2d756036c00eec0bf6449b1e1e3608da Mon Sep 17 00:00:00 2001 From: Steven Gardiner Date: Tue, 25 Jul 2023 15:54:01 -0500 Subject: [PATCH 1/4] Add the MicroBooNE weight calculator that interfaces with geant4reweight --- 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 e99e412ec..c4eb4abc8 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 ec3b03e8b..e7f4e6a2a 100644 --- a/ups/product_deps +++ b/ups/product_deps @@ -261,6 +261,7 @@ sbndata v01_04 - sbnobj v09_17_04 - systematicstools v01_03_01 - nusystematics v01_03_09 - +geant4reweight v01_18_01 - cetmodules v3_21_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 -c7:debug c7:debug c7:debug c7:debug c7:s120a:debug AR2320i00000:e1000:k250 -nq- c7:debug:p3913 c7:debug c7:debug -nq- -c7:prof c7:prof c7:prof c7:prof c7:s120a:prof AR2320i00000:e1000:k250 -nq- c7:p3913:prof c7:prof c7:prof -nq- -e20:debug e20:debug e20:debug e20:debug e20:s120a:debug AR2320i00000:e1000:k250 -nq- debug:e20:p3913 e20:debug e20:debug -nq- -e20:prof e20:prof e20:prof e20:prof e20:s120a:prof AR2320i00000:e1000:k250 -nq- e20:p3913:prof e20:prof e20:prof -nq- +qualifier larsoft sbnobj sbnanaobj sbndaq_artdaq_core genie_xsec sbndata larcv2 systematicstools nusystematics geant4reweight notes +c7:debug c7:debug c7:debug c7:debug c7:s120a:debug AR2320i00000:e1000:k250 -nq- c7:debug:p3913 c7:debug c7:debug c7:debug:s120a -nq- +c7:prof c7:prof c7:prof c7:prof c7:s120a:prof AR2320i00000:e1000:k250 -nq- c7:p3913:prof c7:prof c7:prof c7:prof:s120a -nq- +e20:debug e20:debug e20:debug e20:debug e20:s120a:debug AR2320i00000:e1000:k250 -nq- debug:e20:p3913 e20:debug e20:debug e20:debug:s120a -nq- +e20:prof e20:prof e20:prof e20:prof e20:s120a:prof AR2320i00000:e1000:k250 -nq- e20:p3913:prof e20:prof e20:prof e20:prof:s120a -nq- end_qualifier_list #################################### From 67ba69a49dc7aef8af6762e941848d075b0e08bf Mon Sep 17 00:00:00 2001 From: Steven Gardiner Date: Tue, 25 Jul 2023 17:25:52 -0500 Subject: [PATCH 2/4] 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 | 473 ++++++++---------- ups/product_deps | 12 +- 3 files changed, 226 insertions(+), 260 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..aaf13c5f3 100644 --- a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx +++ b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx @@ -33,11 +33,11 @@ #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" +#include "geant4reweight/src/ReweightBase/G4ReweighterFactory.hh" +#include "geant4reweight/src/ReweightBase/G4Reweighter.hh" +#include "geant4reweight/src/ReweightBase/G4ReweightTraj.hh" +#include "geant4reweight/src/ReweightBase/G4ReweightStep.hh" +#include "geant4reweight/src/PropBase/G4ReweightParameterMaker.hh" // local include #include "BetheBlochForG4ReweightValid.h" @@ -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 e7f4e6a2a..8fc8b1bdb 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 @@ -261,7 +261,7 @@ sbndata v01_04 - sbnobj v09_17_04 - systematicstools v01_03_01 - nusystematics v01_03_09 - -geant4reweight v01_18_01 - +geant4reweight v01_00_03e - cetmodules v3_21_01 - only_for_build end_product_list #################################### @@ -319,10 +319,10 @@ end_product_list #################################### qualifier larsoft sbnobj sbnanaobj sbndaq_artdaq_core genie_xsec sbndata larcv2 systematicstools nusystematics geant4reweight notes -c7:debug c7:debug c7:debug c7:debug c7:s120a:debug AR2320i00000:e1000:k250 -nq- c7:debug:p3913 c7:debug c7:debug c7:debug:s120a -nq- -c7:prof c7:prof c7:prof c7:prof c7:s120a:prof AR2320i00000:e1000:k250 -nq- c7:p3913:prof c7:prof c7:prof c7:prof:s120a -nq- -e20:debug e20:debug e20:debug e20:debug e20:s120a:debug AR2320i00000:e1000:k250 -nq- debug:e20:p3913 e20:debug e20:debug e20:debug:s120a -nq- -e20:prof e20:prof e20:prof e20:prof e20:s120a:prof AR2320i00000:e1000:k250 -nq- e20:p3913:prof e20:prof e20:prof e20:prof:s120a -nq- +c7:debug c7:debug c7:debug c7:debug c7:debug AR2320i00000:e1000:k250 -nq- c7:debug:p3913 c7:debug c7:debug c7:debug:s120a -nq- +c7:prof c7:prof c7:prof c7:prof c7:prof AR2320i00000:e1000:k250 -nq- c7:p3913:prof c7:prof c7:prof c7:prof:s120a -nq- +e20:debug e20:debug e20:debug e20:debug e20:debug:s120a AR2320i00000:e1000:k250 -nq- debug:e20:p3913 e20:debug e20:debug e20:debug:s120a -nq- +e20:prof e20:prof e20:prof e20:prof e20:prof:s120a AR2320i00000:e1000:k250 -nq- e20:p3913:prof e20:prof e20:prof e20:prof:s120a -nq- end_qualifier_list #################################### From e128b59c7b92483794b4be41743d4ae483d488ee Mon Sep 17 00:00:00 2001 From: Andrew Mastbaum Date: Wed, 26 Jul 2023 15:59:01 -0500 Subject: [PATCH 3/4] 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 aaf13c5f3..a84ca7d5f 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 441c7dbde..7897943c9 100644 --- a/sbncode/SBNEventWeight/jobs/CMakeLists.txt +++ b/sbncode/SBNEventWeight/jobs/CMakeLists.txt @@ -4,4 +4,5 @@ install_source() 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 d870e79453e50c7b70f9f1b17122cc77c70fa4da Mon Sep 17 00:00:00 2001 From: Andrew Mastbaum Date: Thu, 27 Jul 2023 13:40:59 -0500 Subject: [PATCH 4/4] 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 a84ca7d5f..a0cead5d2 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 +