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/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/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/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..f05f8d274 --- /dev/null +++ b/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt @@ -0,0 +1,23 @@ +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 + ROOT::Tree +) + +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..a0cead5d2 --- /dev/null +++ b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx @@ -0,0 +1,492 @@ +/** + * \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/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" + +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, size_t inu); + +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 + + 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); + } + + + // 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; + + fParameterSet.Configure(GetFullName(), mode, fNsims); + + 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; + + fParameterSet.AddParameter(theName, theSigma); + } + + 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, size_t itruth ) { + + // 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; + + // Initialize weight vector for this MCTruth + weight.clear(); + weight.resize(fNsims, 1.0); + + // Loop over MCParticles in the event + auto const& mcparticles = truthParticles.at(itruth); + + for (size_t i=0; i 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[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[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; jFill(); + } // loop over mcparticles (i) + if (fMakeOutputTrees) fOutTree_MCTruth->Fill(); + +return weight; + +} + +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_geant4_sbn.fcl b/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl new file mode 100644 index 000000000..e45065880 --- /dev/null +++ b/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl @@ -0,0 +1,103 @@ +########################################################## +## 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" + +BEGIN_PROLOG + +sbn_eventweight_geant4: { + 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) +### + +END_PROLOG + 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 diff --git a/ups/product_deps b/ups/product_deps index ec3b03e8b..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,6 +261,7 @@ sbndata v01_04 - sbnobj v09_17_04 - systematicstools v01_03_01 - nusystematics v01_03_09 - +geant4reweight v01_00_03e - 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: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 ####################################