From 706c4d8895a1458164feed5a31de6924a19d4dcf Mon Sep 17 00:00:00 2001 From: Jaesung Kim Date: Fri, 2 Jun 2023 12:25:54 -0500 Subject: [PATCH 01/28] Update SystToolsEventWeight for dependent dials --- .../App/SystToolsEventWeight_module.cc | 102 +++++++++++++++--- 1 file changed, 87 insertions(+), 15 deletions(-) diff --git a/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc b/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc index 4793171b9..78fee3c44 100644 --- a/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc +++ b/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc @@ -116,19 +116,22 @@ void SystToolsEventWeight::produce(art::Event& e) { << sp->GetFullyQualifiedName() << "\n"; } - //==== syst_resp->size() is (Number of MCTruth) - //==== Each index corresponds to each of MCTruth + // syst_resp->size() is (Number of MCTruth) + // Each index corresponds to each of MCTruth int nMCTruthIndex = 0; if(fDebugMode) std::cout << "[SystToolsEventWeight::produce] syst_resp.size() (= Number of MCTruth) of this SystProvider = " << syst_resp->size() << std::endl; - //==== Looping over syst_resp is identical to looping over MCTruth + // Looping over syst_resp is identical to looping over MCTruth for(systtools::event_unit_response_t const& resp: *syst_resp) { - //==== resp.size() corresponds to number of knobs we altered; - //==== e.g., MaCCQE, MaCCRES, MvCCRE -> resp.size() = 3 + // resp.size() corresponds to number of knobs we altered; + // e.g., MaCCQE, MaCCRES, MvCCRE -> resp.size() = 3 if(fDebugMode){ std::cout << "[SystToolsEventWeight::produce] sp->GetSystMetaData().size() (expected) = " << sp->GetSystMetaData().size() << "\n" << "[SystToolsEventWeight::produce] resp.size() of this syst_resp (produced) = " << resp.size() << std::endl; } + + // Below check is not valid if there is a dependent dial + /* if(sp->GetSystMetaData().size()!=resp.size()){ throw cet::exception{ "SystToolsEventWeight" } << "sp->GetFullyQualifiedName() = " << sp->GetFullyQualifiedName() << "\n" @@ -136,12 +139,13 @@ void SystToolsEventWeight::produce(art::Event& e) { << resp.size() << " are produced. " << "Probably this particular event is not relevant to this systematic variation.\n"; } + */ sbn::evwgh::EventWeightMap mcwgh; for( auto const& r: resp){ - //==== responses.size() is the number of universes + // responses.size() is the number of universes systtools::SystParamHeader const &sph = fParamHeaderHelper.GetHeader( r.pid ); std::string prettyName = sph.prettyName; @@ -159,8 +163,8 @@ void SystToolsEventWeight::produce(art::Event& e) { << "Probably this particular event is not relevant to this systematic variation.\n"; } - //==== r.responses : std::vector - //==== std::map > EventWeightMap + // r.responses : std::vector + // std::map > EventWeightMap mcwgh[sp->GetFullyQualifiedName()+"_"+prettyName].assign(r.responses.cbegin(), r.responses.cend()); } // END loop over knobs @@ -191,24 +195,92 @@ void SystToolsEventWeight::beginRun(art::Run& run) { << "sp->GetFullyQualifiedName() = " << sp->GetFullyQualifiedName() << "\n" << "sp->GetInstanceName() = " << sp->GetInstanceName() << "\n" << "Printing each SystParamHeader of this ISystProviderTool.."; - //==== Note: typedef std::vector SystMetaData; + // Note: typedef std::vector SystMetaData; auto const& smd = sp->GetSystMetaData(); + + // make a map of responsless-response params + std::map> map_resp_to_respless; + for( auto &sph : smd ){ + if(sph.isResponselessParam){ + auto it = map_resp_to_respless.find( sph.responseParamId ); + if( it != map_resp_to_respless.end() ){ + it->second.push_back( sph.systParamId ); + } + else{ + map_resp_to_respless[sph.responseParamId] = {}; + map_resp_to_respless[sph.responseParamId].push_back( sph.systParamId ); + } + } + } + if (fDebugMode) { + for(const auto& it: map_resp_to_respless){ + const auto& sph = systtools::GetParam(smd, it.first); + std::cout << "Found a dependent dial: " << sph.prettyName << std::endl; + for(const auto& depdialid: it.second){ + const auto& sph_dep = systtools::GetParam(smd, depdialid); + std::cout << "- dep dial: " << sph_dep.prettyName << std::endl; + } + } + } + for( auto &sph : smd ){ if (fDebugMode) { std::cout << " sph.prettyName = " << sph.prettyName << std::endl; } - std::string rwmode = "multisigma"; - if(sph.isRandomlyThrown) rwmode = "multisim"; + // responsless + if(sph.isResponselessParam){ + if (fDebugMode) { + std::cout << "Responsless dial found: " << sph.prettyName << ", thus skipping" << std::endl; + } + continue; + } + + sbn::evwgh::EventWeightParameterSet fParameterSet; + std::string rwmode = ""; + + auto it = map_resp_to_respless.find( sph.systParamId ); + if(it!=map_resp_to_respless.end()){ + + for(const auto depdialid: it->second){ + const auto& sph_dep = systtools::GetParam(smd, depdialid); + + if(rwmode=="") rwmode = sph_dep.isRandomlyThrown ? "multisim" : "multisigma"; + else{ + if(rwmode!= (sph_dep.isRandomlyThrown ? "multisim" : "multisigma")){ + throw cet::exception{ "SystToolsEventWeight" } + << sph.prettyName << " depends on other dials, but the remode are different between the deps dials\n"; + } + } + + std::vector widths_dep { sph_dep.paramVariations.begin(), sph_dep.paramVariations.end() }; + fParameterSet.AddParameter(sph_dep.prettyName, std::move(widths_dep)); + + } + + + + } + else{ + + rwmode = sph.isRandomlyThrown ? "multisim" : "multisigma"; + + std::vector widths { sph.paramVariations.begin(), sph.paramVariations.end() }; + + sbn::evwgh::EventWeightParameterSet fParameterSet; + fParameterSet.AddParameter(sph.prettyName, std::move(widths)); + + } + + if(rwmode==""){ + throw cet::exception{ "SystToolsEventWeight" } + << "rwmode not set for " << sph.prettyName << "\n"; + } if (fDebugMode) { std::cout << " rwmode = " << rwmode << std::endl; } - std::vector widths { sph.paramVariations.begin(), sph.paramVariations.end() }; - - sbn::evwgh::EventWeightParameterSet fParameterSet; - fParameterSet.AddParameter(sph.prettyName, std::move(widths)); fParameterSet.Configure(sp->GetFullyQualifiedName()+"_"+sph.prettyName, rwmode, sph.paramVariations.size()); fParameterSet.FillKnobValues(); From f0c3c49448d9aeaadc0df911305103d5bc2a3556 Mon Sep 17 00:00:00 2001 From: Jaesung Kim Date: Mon, 12 Jun 2023 09:38:25 -0500 Subject: [PATCH 02/28] Duplicated declaration of fParameterSet removed --- sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc b/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc index 78fee3c44..ed3c870c0 100644 --- a/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc +++ b/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc @@ -149,6 +149,8 @@ void SystToolsEventWeight::produce(art::Event& e) { systtools::SystParamHeader const &sph = fParamHeaderHelper.GetHeader( r.pid ); std::string prettyName = sph.prettyName; + if(sph.isResponselessParam) continue; + if(fDebugMode){ std::cout << "[SystToolsEventWeight::produce] pid of this resp = " << r.pid << "\n" << "[SystToolsEventWeight::produce] prettyName of this resp = " << prettyName << "\n" @@ -267,7 +269,6 @@ void SystToolsEventWeight::beginRun(art::Run& run) { std::vector widths { sph.paramVariations.begin(), sph.paramVariations.end() }; - sbn::evwgh::EventWeightParameterSet fParameterSet; fParameterSet.AddParameter(sph.prettyName, std::move(widths)); } From 12aa7494968f529df72a2e4e9f740b07d2ee474b Mon Sep 17 00:00:00 2001 From: Jaesung Kim Date: Sat, 24 Jun 2023 22:19:53 -0500 Subject: [PATCH 03/28] Update logger for SystToolsEventWeight_module.cc --- sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc b/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc index ed3c870c0..62f03b6ce 100644 --- a/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc +++ b/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc @@ -226,9 +226,6 @@ void SystToolsEventWeight::beginRun(art::Run& run) { } for( auto &sph : smd ){ - if (fDebugMode) { - std::cout << " sph.prettyName = " << sph.prettyName << std::endl; - } // responsless if(sph.isResponselessParam){ @@ -279,7 +276,7 @@ void SystToolsEventWeight::beginRun(art::Run& run) { } if (fDebugMode) { - std::cout << " rwmode = " << rwmode << std::endl; + std::cout << " sph.prettyName = " << sph.prettyName << ", rwmode = " << rwmode << " is added to the header" << std::endl; } fParameterSet.Configure(sp->GetFullyQualifiedName()+"_"+sph.prettyName, rwmode, sph.paramVariations.size()); From 58c6946968495db9184c2d4777d64fe361436bd0 Mon Sep 17 00:00:00 2001 From: Jaesung Kim Date: Thu, 27 Jul 2023 13:24:07 -0500 Subject: [PATCH 04/28] Update SystToolsEventWeight module to treat isCorrection=true paramheader --- .../App/SystToolsEventWeight_module.cc | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc b/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc index 62f03b6ce..6e836054c 100644 --- a/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc +++ b/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc @@ -157,7 +157,14 @@ void SystToolsEventWeight::produce(art::Event& e) { << "[SystToolsEventWeight::produce] paramVariations.size() of this resp (expected) = " << sph.paramVariations.size() << "\n" << "[SystToolsEventWeight::produce] responses.size() of this resp (produced) = " << r.responses.size() << std::endl; } - if(sph.paramVariations.size()!=r.responses.size()){ + + if(sph.isCorrection && (r.responses.size()!=1)){ + throw cet::exception{ "SystToolsEventWeight" } + << "prettyName of this resp = " << prettyName << "\n" + << "We expect to have 1 universes as it is a correction, but " + << r.responses.size() << " are produced\n"; + } + if(!sph.isCorrection && sph.paramVariations.size()!=r.responses.size()){ throw cet::exception{ "SystToolsEventWeight" } << "prettyName of this resp = " << prettyName << "\n" << "We expect to have " << sph.paramVariations.size() << " universes, but " @@ -238,11 +245,14 @@ void SystToolsEventWeight::beginRun(art::Run& run) { sbn::evwgh::EventWeightParameterSet fParameterSet; std::string rwmode = ""; + std::vector paramVars = sph.isCorrection ? std::vector(1, sph.centralParamValue) : sph.paramVariations; + auto it = map_resp_to_respless.find( sph.systParamId ); if(it!=map_resp_to_respless.end()){ for(const auto depdialid: it->second){ const auto& sph_dep = systtools::GetParam(smd, depdialid); + std::vector paramVars_dep = sph_dep.isCorrection ? std::vector(1, sph_dep.centralParamValue) : sph_dep.paramVariations; if(rwmode=="") rwmode = sph_dep.isRandomlyThrown ? "multisim" : "multisigma"; else{ @@ -252,7 +262,7 @@ void SystToolsEventWeight::beginRun(art::Run& run) { } } - std::vector widths_dep { sph_dep.paramVariations.begin(), sph_dep.paramVariations.end() }; + std::vector widths_dep { paramVars_dep.begin(), paramVars_dep.end() }; fParameterSet.AddParameter(sph_dep.prettyName, std::move(widths_dep)); } @@ -264,7 +274,7 @@ void SystToolsEventWeight::beginRun(art::Run& run) { rwmode = sph.isRandomlyThrown ? "multisim" : "multisigma"; - std::vector widths { sph.paramVariations.begin(), sph.paramVariations.end() }; + std::vector widths { paramVars.begin(), paramVars.end() }; fParameterSet.AddParameter(sph.prettyName, std::move(widths)); @@ -279,7 +289,7 @@ void SystToolsEventWeight::beginRun(art::Run& run) { std::cout << " sph.prettyName = " << sph.prettyName << ", rwmode = " << rwmode << " is added to the header" << std::endl; } - fParameterSet.Configure(sp->GetFullyQualifiedName()+"_"+sph.prettyName, rwmode, sph.paramVariations.size()); + fParameterSet.Configure(sp->GetFullyQualifiedName()+"_"+sph.prettyName, rwmode, paramVars.size()); fParameterSet.FillKnobValues(); p->push_back( std::move(fParameterSet) ); From fd58a141f178b7bdf57b576ed493cbd102d16cc1 Mon Sep 17 00:00:00 2001 From: Jaesung Kim Date: Thu, 15 Feb 2024 15:20:41 -0600 Subject: [PATCH 05/28] Saving non primaries to nu.prim for G4 study --- sbncode/CAFMaker/FillTrue.cxx | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/sbncode/CAFMaker/FillTrue.cxx b/sbncode/CAFMaker/FillTrue.cxx index 3b9da6340..5d7393832 100644 --- a/sbncode/CAFMaker/FillTrue.cxx +++ b/sbncode/CAFMaker/FillTrue.cxx @@ -470,7 +470,9 @@ namespace caf { for(const caf::SRTrueParticle& part: srparticles){ // save the G4 particles that came from this interaction if(part.interaction_id == (int)i) { - if(part.start_process == caf::kG4primary) srneutrino.prim.push_back(part); + //if(part.start_process == caf::kG4primary) srneutrino.prim.push_back(part); + // Feb. 15th 2024, Jaesung Kim: Saving non-primaries for G4Reweight study + srneutrino.prim.push_back(part); // total up the deposited energy for(int p = 0; p < 3; ++p) { From 43893b9b84e91cf3f2aa6676f1ab7cded4317d8a Mon Sep 17 00:00:00 2001 From: Jaesung Kim Date: Thu, 15 Feb 2024 22:53:55 -0600 Subject: [PATCH 06/28] Update systematicstools/nusystematics versions --- ups/product_deps | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ups/product_deps b/ups/product_deps index 8d2b04768..3c0df5cca 100644 --- a/ups/product_deps +++ b/ups/product_deps @@ -259,8 +259,8 @@ sbnanaobj v09_20_06 - sbndaq_artdaq_core v1_06_00of0 - sbndata v01_04 - sbnobj v09_16_00 - -systematicstools v01_02_00 - -nusystematics v01_02_10 - +systematicstools v01_02_01 - +nusystematics v01_02_11 - cetmodules v3_20_00 - only_for_build end_product_list #################################### From 658fb3bbfe467951db12dc26e1254020d77b42c5 Mon Sep 17 00:00:00 2001 From: Steven Gardiner Date: Tue, 25 Jul 2023 15:54:01 -0500 Subject: [PATCH 07/28] 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 183b0b83b..0b969ec7f 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 7457d3fbf..98dfb2983 100644 --- a/ups/product_deps +++ b/ups/product_deps @@ -261,6 +261,7 @@ sbndata v01_04 - sbnobj v09_16_00_02 - systematicstools v01_02_00 - nusystematics v01_02_10 - +geant4reweight v01_14_04 - cetmodules v3_20_00 - 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:s117:debug AR2320i00000:e1000:k250 -nq- c7:p392:debug c7:debug c7:debug -nq- -c7:prof c7:prof c7:prof c7:prof c7:s117:prof AR2320i00000:e1000:k250 -nq- c7:p392:prof c7:prof c7:prof -nq- -e20:debug e20:debug e20:debug e20:debug e20:s117:debug AR2320i00000:e1000:k250 -nq- e20:p392:debug e20:debug e20:debug -nq- -e20:prof e20:prof e20:prof e20:prof e20:s117:prof AR2320i00000:e1000:k250 -nq- e20:p392: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:s117:debug AR2320i00000:e1000:k250 -nq- c7:p392:debug c7:debug c7:debug c7:s117:debug -nq- +c7:prof c7:prof c7:prof c7:prof c7:s117:prof AR2320i00000:e1000:k250 -nq- c7:p392:prof c7:prof c7:prof c7:s117:prof -nq- +e20:debug e20:debug e20:debug e20:debug e20:s117:debug AR2320i00000:e1000:k250 -nq- e20:p392:debug e20:debug e20:debug e20:s117:debug -nq- +e20:prof e20:prof e20:prof e20:prof e20:s117:prof AR2320i00000:e1000:k250 -nq- e20:p392:prof e20:prof e20:prof e20:s117:prof -nq- end_qualifier_list #################################### From 1f3c410ebcfff03bbddec876da215b724f5cb4dd Mon Sep 17 00:00:00 2001 From: Steven Gardiner Date: Tue, 25 Jul 2023 17:25:52 -0500 Subject: [PATCH 08/28] 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:s117:prof --- .../Calculators/Geant4/CMakeLists.txt | 1 + .../Calculators/Geant4/Geant4WeightCalc.cxx | 473 ++++++++---------- ups/product_deps | 4 +- 3 files changed, 222 insertions(+), 256 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 98dfb2983..8963bf559 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_16_00_02 - systematicstools v01_02_00 - nusystematics v01_02_10 - -geant4reweight v01_14_04 - +geant4reweight v01_00_03e - cetmodules v3_20_00 - only_for_build end_product_list #################################### From 38da5270785337518dadbe9bcbcd23d00d3966a5 Mon Sep 17 00:00:00 2001 From: Andrew Mastbaum Date: Wed, 26 Jul 2023 15:59:01 -0500 Subject: [PATCH 09/28] 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 ce8b26c71caf8eaf65c086afa73bac7ca8173f38 Mon Sep 17 00:00:00 2001 From: Andrew Mastbaum Date: Thu, 27 Jul 2023 13:40:59 -0500 Subject: [PATCH 10/28] 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 + From dc8fffddb06ddd8b50a4044600674ea9c130ba1d Mon Sep 17 00:00:00 2001 From: gputnam Date: Fri, 28 Jul 2023 16:23:21 -0500 Subject: [PATCH 11/28] Fill true trajectory points. --- .../Calibration/TrackCaloSkimmer_module.cc | 35 +++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/sbncode/Calibration/TrackCaloSkimmer_module.cc b/sbncode/Calibration/TrackCaloSkimmer_module.cc index 8cbfcaa19..b31b86c9e 100644 --- a/sbncode/Calibration/TrackCaloSkimmer_module.cc +++ b/sbncode/Calibration/TrackCaloSkimmer_module.cc @@ -831,6 +831,41 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle, } } + // Save the true trajectory + for (unsigned i_traj = 0; i_traj < particle.NumberTrajectoryPoints(); i_traj++) { + // Get trajectory point + TVector3 traj = particle.Position(i_traj).Vect(); + geo::Point_t traj_p(traj.X(), traj.Y(), traj.Z()); + + // lookup TPC + geo::TPCID tpc; // invalid by default + for (auto const &cryo: geo->Iterate()) { + for (auto const& TPC : geo->Iterate(cryo.ID())) { + if (TPC.ActiveBoundingBox().ContainsPosition(traj_p)) { + tpc = TPC.ID(); + break; + } + } + if (tpc.isValid) break; + } + + // add in space-charge-deflected position if applicable + geo::Point_t traj_p_sce = tpc.isValid ? TrajectoryToWirePosition(traj_p, tpc) : traj_p; + + sbn::Vector3D traj_v; + traj_v.x = traj_p.x(); + traj_v.y = traj_p.y(); + traj_v.z = traj_p.z(); + + sbn::Vector3D traj_v_sce; + traj_v_sce.x = traj_p_sce.x(); + traj_v_sce.y = traj_p_sce.y(); + traj_v_sce.z = traj_p_sce.z(); + + trueparticle.traj.push_back(traj_v); + trueparticle.traj_sce.push_back(traj_v_sce); + } + return trueparticle; } From 53c9ee92566d09a5efb15054aef2fcb901ea4dff Mon Sep 17 00:00:00 2001 From: gputnam Date: Fri, 28 Jul 2023 16:06:04 -0500 Subject: [PATCH 12/28] Disable cross-plane stub merging. --- sbncode/TPCReco/VertexStub/config/sbn_stub_merge_tools.fcl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/sbncode/TPCReco/VertexStub/config/sbn_stub_merge_tools.fcl b/sbncode/TPCReco/VertexStub/config/sbn_stub_merge_tools.fcl index ec4bb6a13..78e14279f 100644 --- a/sbncode/TPCReco/VertexStub/config/sbn_stub_merge_tools.fcl +++ b/sbncode/TPCReco/VertexStub/config/sbn_stub_merge_tools.fcl @@ -13,6 +13,8 @@ twoplane_stub_merge: { SaveOldStubs: false } -stub_merge: [@local::plane_stub_merge, @local::twoplane_stub_merge] +# stub_merge: [@local::plane_stub_merge, @local::twoplane_stub_merge] +# NOTE: disable multi-plane matching while other planes are not calibrated +stub_merge: [@local::plane_stub_merge] END_PROLOG From 2eebbe0ebd6da1aab6c88395192aa7ed22baced7 Mon Sep 17 00:00:00 2001 From: Guadalupe Moreno Granados Date: Tue, 19 Sep 2023 09:39:04 -0500 Subject: [PATCH 13/28] For cosmics that are crossing the cathode we added T0 information in calorimetry variables --- sbncode/CAFMaker/CAFMaker_module.cc | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index 8a44dfc2c..d45fa0a7e 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/sbncode/CAFMaker/CAFMaker_module.cc @@ -1828,9 +1828,10 @@ void CAFMaker::produce(art::Event& evt) noexcept { if (fmTrackDazzle.isValid() && fmTrackDazzle.at(iPart).size()==1) { FillTrackDazzle(fmTrackDazzle.at(iPart).front(), trk); } + std::cout<< "t0 calo: " << pfp.t0 <(), dprop, trk); } From 5dd91405dc12a66b11dfc9da22caff2b88efc5c3 Mon Sep 17 00:00:00 2001 From: Mun Jung Jung Date: Fri, 15 Sep 2023 16:24:11 -0500 Subject: [PATCH 14/28] exclude region feature --- .../LArReco/TrajectoryMCSFitter.cxx | 39 ++++++++++++++++--- .../LArReco/TrajectoryMCSFitter.h | 14 +++++-- sbncode/LArRecoProducer/mcsproducer.fcl | 4 +- .../LArRecoProducer/mcsproducer_icarus.fcl | 12 ++++++ 4 files changed, 59 insertions(+), 10 deletions(-) create mode 100644 sbncode/LArRecoProducer/mcsproducer_icarus.fcl diff --git a/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx b/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx index 8ae583973..d07dcd630 100644 --- a/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx +++ b/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx @@ -15,7 +15,8 @@ recob::MCSFitResult TrajectoryMCSFitter::fitMcs(const recob::TrackTrajectory& tr vector breakpoints; vector segradlengths; vector cumseglens; - breakTrajInSegments(traj, breakpoints, segradlengths, cumseglens); + vector breakpointsgood; + breakTrajInSegments(traj, breakpoints, segradlengths, cumseglens, breakpointsgood); // // Fit segment directions, and get 3D angles between them // @@ -28,7 +29,11 @@ recob::MCSFitResult TrajectoryMCSFitter::fitMcs(const recob::TrackTrajectory& tr if (p>0) { if (segradlengths[p]<-100. || segradlengths[p-1]<-100.) { dtheta.push_back(-999.); - } else { + } + if (!breakpointsgood[p] || !breakpointsgood[p-1]) { + dtheta.push_back(-999.); + } + else { const double cosval = pcdir0.X()*pcdir1.X()+pcdir0.Y()*pcdir1.Y()+pcdir0.Z()*pcdir1.Z(); //assert(std::abs(cosval)<=1); //units are mrad @@ -49,15 +54,19 @@ recob::MCSFitResult TrajectoryMCSFitter::fitMcs(const recob::TrackTrajectory& tr } const ScanResult fwdResult = doLikelihoodScan(dtheta, segradlengths, cumLenFwd, true, momDepConst, pid); const ScanResult bwdResult = doLikelihoodScan(dtheta, segradlengths, cumLenBwd, false, momDepConst, pid); - // + // std::cout << "fwdResult.p=" << fwdResult.p << " fwdResult.pUnc=" << fwdResult.pUnc << " fwdResult.logL=" << fwdResult.logL << std::endl; return recob::MCSFitResult(pid, fwdResult.p,fwdResult.pUnc,fwdResult.logL, bwdResult.p,bwdResult.pUnc,bwdResult.logL, segradlengths,dtheta); } -void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj, vector& breakpoints, vector& segradlengths, vector& cumseglens) const { - // +void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj, vector& breakpoints, vector& segradlengths, vector& cumseglens, vector& breakpointsgood) const { + // check that length of excludeRegions vector is in the correct format (a multiple of 6) + bool excludeRegionsSizeGood = (excludeRegions_.size()%6==0); + if (!excludeRegionsSizeGood) { + std::cout << "Error: excludeRegions vector must have length multiple of 6, not excluding any regions" << std::endl; + } const double trajlen = traj.Length(); const int nseg = std::max(minNSegs_,int(trajlen/segLen_)); const double thisSegLen = trajlen/double(nseg); @@ -69,6 +78,8 @@ void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj auto nextValid=traj.FirstValidPoint(); breakpoints.push_back(nextValid); auto pos0 = traj.LocationAtPoint(nextValid); + bool pos0good = isInGoodRegion(excludeRegionsSizeGood, pos0.X(), pos0.Y(), pos0.Z()); + breakpointsgood.push_back(pos0good); nextValid = traj.NextValidPoint(nextValid+1); int npoints = 0; while (nextValid!=recob::TrackTrajectory::InvalidIndex) { @@ -77,7 +88,9 @@ void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj pos0=pos1; npoints++; if (thislen>=thisSegLen) { + bool pos1good = isInGoodRegion(excludeRegionsSizeGood, pos1.X(), pos1.Y(), pos1.Z()); breakpoints.push_back(nextValid); + breakpointsgood.push_back(pos1good); if (npoints>=minHitsPerSegment_) segradlengths.push_back(thislen*lar_radl_inv); else segradlengths.push_back(-999.); cumseglens.push_back(cumseglens.back()+thislen); @@ -88,7 +101,10 @@ void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj } //then add last segment if (thislen>0.) { + auto endpointpos = traj.LocationAtPoint(nextValid); + bool endpointposgood = isInGoodRegion(excludeRegionsSizeGood, endpointpos.X(), endpointpos.Y(), endpointpos.Z()); breakpoints.push_back(traj.LastValidPoint()+1); + breakpointsgood.push_back(endpointposgood); segradlengths.push_back(thislen*lar_radl_inv); cumseglens.push_back(cumseglens.back()+thislen); } @@ -202,7 +218,6 @@ double TrajectoryMCSFitter::mcsLikelihood(double p, double theta0x, std::vector< double result = 0; for (int i = beg; i != end; i+=incr ) { if (dthetaij[i]<0) { - //cout << "skip segment with too few points" << endl; continue; } // @@ -302,3 +317,15 @@ double TrajectoryMCSFitter::GetE(const double initial_E, const double length_tra } return current_E; } + +bool TrajectoryMCSFitter::isInGoodRegion(const bool sizeGood, const double x, const double y, const double z) const { + if (!sizeGood) return true; // if in wrong foramt, do not exclude any regions + for (unsigned int i=0; i=excludeRegions_[6*i+0] && x<=excludeRegions_[6*i+1] && + y>=excludeRegions_[6*i+2] && y<=excludeRegions_[6*i+3] && + z>=excludeRegions_[6*i+4] && z<=excludeRegions_[6*i+5]) { + return false; + } + } + return true; +} diff --git a/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h b/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h index 07b63c443..1ea2f1dd8 100644 --- a/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h +++ b/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h @@ -3,6 +3,7 @@ #include "fhiclcpp/ParameterSet.h" #include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/Sequence.h" #include "fhiclcpp/types/Table.h" #include "canvas/Persistency/Common/Ptr.h" #include "lardataobj/RecoBase/MCSFitResult.h" @@ -87,10 +88,14 @@ namespace trkf::sbn { Comment("Angular resolution parameter used in modified Highland formula. Unit is mrad."), 3.0 }; + fhicl::Sequence excludeRegions { + Name("excludeRegions"), + Comment("regions to exclude") + }; }; using Parameters = fhicl::Table; // - TrajectoryMCSFitter(int pIdHyp, int minNSegs, double segLen, int minHitsPerSegment, int nElossSteps, int eLossMode, double pMin, double pMax, double pStep, double angResol){ + TrajectoryMCSFitter(int pIdHyp, int minNSegs, double segLen, int minHitsPerSegment, int nElossSteps, int eLossMode, double pMin, double pMax, double pStep, double angResol, std::vector excludeRegions){ pIdHyp_ = pIdHyp; minNSegs_ = minNSegs; segLen_ = segLen; @@ -101,9 +106,10 @@ namespace trkf::sbn { pMax_ = pMax; pStep_ = pStep; angResol_ = angResol; + excludeRegions_ = excludeRegions; } explicit TrajectoryMCSFitter(const Parameters & p) - : TrajectoryMCSFitter(p().pIdHypothesis(),p().minNumSegments(),p().segmentLength(),p().minHitsPerSegment(),p().nElossSteps(),p().eLossMode(),p().pMin(),p().pMax(),p().pStep(),p().angResol()) {} + : TrajectoryMCSFitter(p().pIdHypothesis(),p().minNumSegments(),p().segmentLength(),p().minHitsPerSegment(),p().nElossSteps(),p().eLossMode(),p().pMin(),p().pMax(),p().pStep(),p().angResol(),p().excludeRegions()) {} // recob::MCSFitResult fitMcs(const recob::TrackTrajectory& traj, bool momDepConst = true) const { return fitMcs(traj,pIdHyp_,momDepConst); } recob::MCSFitResult fitMcs(const recob::Track& track, bool momDepConst = true) const { return fitMcs(track,pIdHyp_,momDepConst); } @@ -117,7 +123,7 @@ namespace trkf::sbn { return fitMcs(tt,pid,momDepConst); } // - void breakTrajInSegments(const recob::TrackTrajectory& traj, std::vector& breakpoints, std::vector& segradlengths, std::vector& cumseglens) const; + void breakTrajInSegments(const recob::TrackTrajectory& traj, std::vector& breakpoints, std::vector& segradlengths, std::vector& cumseglens, std::vector& breakpointsgood) const; void linearRegression(const recob::TrackTrajectory& traj, const size_t firstPoint, const size_t lastPoint, recob::tracking::Vector_t& pcdir) const; double mcsLikelihood(double p, double theta0x, std::vector& dthetaij, std::vector& seg_nradl, std::vector& cumLen, bool fwd, bool momDepConst, int pid) const; // @@ -146,6 +152,7 @@ namespace trkf::sbn { double energyLossLandau(const double mass2,const double E2, const double x) const; // double GetE(const double initial_E, const double length_travelled, const double mass) const; + bool isInGoodRegion(const bool sizeGood, const double x, const double y, const double z) const; // private: int pIdHyp_; @@ -158,6 +165,7 @@ namespace trkf::sbn { double pMax_; double pStep_; double angResol_; + std::vector excludeRegions_; }; } diff --git a/sbncode/LArRecoProducer/mcsproducer.fcl b/sbncode/LArRecoProducer/mcsproducer.fcl index e41a6d00b..31ca1df13 100644 --- a/sbncode/LArRecoProducer/mcsproducer.fcl +++ b/sbncode/LArRecoProducer/mcsproducer.fcl @@ -1,7 +1,9 @@ BEGIN_PROLOG mcs_sbn: { module_type: MCSFitAllPID - MCS: {} + MCS: { + excludeRegions: [] + } TrackLabel: pandoraTrack } diff --git a/sbncode/LArRecoProducer/mcsproducer_icarus.fcl b/sbncode/LArRecoProducer/mcsproducer_icarus.fcl new file mode 100644 index 000000000..09e04041b --- /dev/null +++ b/sbncode/LArRecoProducer/mcsproducer_icarus.fcl @@ -0,0 +1,12 @@ +BEGIN_PROLOG +mcs_sbn: { + module_type: MCSFitAllPID + MCS: { + eLossMode: 2 + angResol: 10.0 + excludeRegions: [190.0, 240.0, -9999.0, 9999.0, -9999.0, 9999.0] + } + TrackLabel: pandoraTrack +} + +END_PROLOG From b7657b38c9fa0347a6b5126cc71e9920fe15e0e5 Mon Sep 17 00:00:00 2001 From: Mun Jung Jung Date: Mon, 18 Sep 2023 18:32:39 -0500 Subject: [PATCH 15/28] restructure, skip out of fv hits --- .../LArReco/TrajectoryMCSFitter.cxx | 95 +++++++++++++++---- .../LArReco/TrajectoryMCSFitter.h | 25 +++-- sbncode/LArRecoProducer/mcsproducer.fcl | 3 +- .../LArRecoProducer/mcsproducer_icarus.fcl | 3 +- 4 files changed, 101 insertions(+), 25 deletions(-) diff --git a/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx b/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx index d07dcd630..cd0b4c53f 100644 --- a/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx +++ b/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx @@ -1,6 +1,9 @@ #include "TrajectoryMCSFitter.h" #include "lardataobj/RecoBase/Track.h" #include "larcorealg/Geometry/geo_vectors_utils.h" +#include "larcore/Geometry/Geometry.h" +#include "larcore/CoreUtils/ServiceUtil.h" +#include "larcorealg/Geometry/GeometryCore.h" #include "TMatrixDSym.h" #include "TMatrixDSymEigen.h" @@ -62,11 +65,13 @@ recob::MCSFitResult TrajectoryMCSFitter::fitMcs(const recob::TrackTrajectory& tr } void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj, vector& breakpoints, vector& segradlengths, vector& cumseglens, vector& breakpointsgood) const { - // check that length of excludeRegions vector is in the correct format (a multiple of 6) - bool excludeRegionsSizeGood = (excludeRegions_.size()%6==0); - if (!excludeRegionsSizeGood) { - std::cout << "Error: excludeRegions vector must have length multiple of 6, not excluding any regions" << std::endl; - } + // set fiducial volume + std::vector fiducialVolumes; + fiducialVolumes = setFiducialVolumes(); + // set excluded volumes + std::vector excludeVolumes; + excludeVolumes = setExcludeVolumes(); + const double trajlen = traj.Length(); const int nseg = std::max(minNSegs_,int(trajlen/segLen_)); const double thisSegLen = trajlen/double(nseg); @@ -78,7 +83,7 @@ void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj auto nextValid=traj.FirstValidPoint(); breakpoints.push_back(nextValid); auto pos0 = traj.LocationAtPoint(nextValid); - bool pos0good = isInGoodRegion(excludeRegionsSizeGood, pos0.X(), pos0.Y(), pos0.Z()); + bool pos0good = isInVolume(fiducialVolumes, pos0) && !isInVolume(excludeVolumes, pos0); breakpointsgood.push_back(pos0good); nextValid = traj.NextValidPoint(nextValid+1); int npoints = 0; @@ -88,7 +93,7 @@ void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj pos0=pos1; npoints++; if (thislen>=thisSegLen) { - bool pos1good = isInGoodRegion(excludeRegionsSizeGood, pos1.X(), pos1.Y(), pos1.Z()); + bool pos1good = isInVolume(fiducialVolumes, pos1) && !isInVolume(excludeVolumes, pos1); breakpoints.push_back(nextValid); breakpointsgood.push_back(pos1good); if (npoints>=minHitsPerSegment_) segradlengths.push_back(thislen*lar_radl_inv); @@ -102,7 +107,7 @@ void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj //then add last segment if (thislen>0.) { auto endpointpos = traj.LocationAtPoint(nextValid); - bool endpointposgood = isInGoodRegion(excludeRegionsSizeGood, endpointpos.X(), endpointpos.Y(), endpointpos.Z()); + bool endpointposgood = isInVolume(fiducialVolumes, endpointpos) && !isInVolume(excludeVolumes, endpointpos); breakpoints.push_back(traj.LastValidPoint()+1); breakpointsgood.push_back(endpointposgood); segradlengths.push_back(thislen*lar_radl_inv); @@ -318,14 +323,72 @@ double TrajectoryMCSFitter::GetE(const double initial_E, const double length_tra return current_E; } -bool TrajectoryMCSFitter::isInGoodRegion(const bool sizeGood, const double x, const double y, const double z) const { - if (!sizeGood) return true; // if in wrong foramt, do not exclude any regions - for (unsigned int i=0; i=excludeRegions_[6*i+0] && x<=excludeRegions_[6*i+1] && - y>=excludeRegions_[6*i+2] && y<=excludeRegions_[6*i+3] && - z>=excludeRegions_[6*i+4] && z<=excludeRegions_[6*i+5]) { - return false; +std::vector TrajectoryMCSFitter::setFiducialVolumes() const { + std::vector fiducialVolumes; + std::vector> TPCVolumes; + std::vector ActiveVolumes; + + const geo::GeometryCore *geometry = lar::providerFrom(); + for (auto const &cryoID: geometry->Iterate()) { + std::vector thisTPCVolumes; + for (auto const& tpc: geometry->Iterate(cryoID)) { + thisTPCVolumes.push_back(tpc.ActiveBoundingBox()); } + TPCVolumes.push_back(std::move(thisTPCVolumes)); } - return true; + for (const std::vector &tpcs: TPCVolumes) { + double xMin = std::min_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MinX() < rhs.MinX(); })->MinX(); + double xMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxX() < rhs.MaxX(); })->MaxX(); + double yMin = std::min_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MinY() < rhs.MinY(); })->MinY(); + double yMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxY() < rhs.MaxY(); })->MaxY(); + double zMin = std::min_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MinZ() < rhs.MinZ(); })->MinZ(); + double zMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxZ() < rhs.MaxZ(); })->MaxZ(); + ActiveVolumes.emplace_back(xMin, xMax, yMin, yMax, zMin, zMax); + } + double fidInsetMinX = fiducialVolumeInsets_[0]; + double fidInsetMaxX = fiducialVolumeInsets_[1]; + double fidInsetMinY = fiducialVolumeInsets_[2]; + double fidInsetMaxY = fiducialVolumeInsets_[3]; + double fidInsetMinZ = fiducialVolumeInsets_[4]; + double fidInsetMaxZ = fiducialVolumeInsets_[5]; + + if (fiducialVolumeInsets_.size() != 6) { + std::cout << "Error: fiducialVolumeInsets vector must have length of 6, not fiducializing" << std::endl; + fidInsetMinX = 0.0; + fidInsetMaxX = 0.0; + fidInsetMinY = 0.0; + fidInsetMaxY = 0.0; + fidInsetMinZ = 0.0; + fidInsetMaxZ = 0.0; + } + for (const geo::BoxBoundedGeo &AV: ActiveVolumes) { + fiducialVolumes.emplace_back(AV.MinX() + fidInsetMinX, AV.MaxX() - fidInsetMaxX, + AV.MinY() + fidInsetMinY, AV.MaxY() - fidInsetMaxY, + AV.MinZ() + fidInsetMinZ, AV.MaxZ() - fidInsetMaxZ); + } + return fiducialVolumes; } + +std::vector TrajectoryMCSFitter::setExcludeVolumes() const { + std::vector excludeVolumes; + if (excludeVolumes_.size()%6 != 0) { + std::cout << "Error: excludeVolumes vector must have length multiple of 6, not excluding any regions" << std::endl; + excludeVolumes.emplace_back(-9999.0, -9999.0, -9999.0, -9999.0, -9999.0, -9999.0); + return excludeVolumes; + } + for (unsigned int i=0; i &volumes, const geo::Point_t &point) const { + for (const geo::BoxBoundedGeo &volume: volumes) { + if (point.X()>=volume.MinX() && point.X()<=volume.MaxX() && + point.Y()>=volume.MinY() && point.Y()<=volume.MaxY() && + point.Z()>=volume.MinZ() && point.Z()<=volume.MaxZ()) { + return true; + } + } + return false; +} \ No newline at end of file diff --git a/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h b/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h index 1ea2f1dd8..6d3b48032 100644 --- a/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h +++ b/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h @@ -9,6 +9,9 @@ #include "lardataobj/RecoBase/MCSFitResult.h" #include "lardataobj/RecoBase/Track.h" #include "lardata/RecoObjects/TrackState.h" +#include "larcore/Geometry/Geometry.h" +#include "larcore/CoreUtils/ServiceUtil.h" +#include "larcorealg/Geometry/GeometryCore.h" namespace trkf::sbn { /** @@ -88,14 +91,18 @@ namespace trkf::sbn { Comment("Angular resolution parameter used in modified Highland formula. Unit is mrad."), 3.0 }; - fhicl::Sequence excludeRegions { - Name("excludeRegions"), + fhicl::Sequence fiducialVolumeInsets { + Name("fiducialVolumeInsets"), + Comment("insets from the TPC boundaries to exclude from the fit") + }; + fhicl::Sequence excludeVolumes { + Name("excludeVolumes"), Comment("regions to exclude") }; }; using Parameters = fhicl::Table; // - TrajectoryMCSFitter(int pIdHyp, int minNSegs, double segLen, int minHitsPerSegment, int nElossSteps, int eLossMode, double pMin, double pMax, double pStep, double angResol, std::vector excludeRegions){ + TrajectoryMCSFitter(int pIdHyp, int minNSegs, double segLen, int minHitsPerSegment, int nElossSteps, int eLossMode, double pMin, double pMax, double pStep, double angResol, std::vectorfiducialVolumeInsets, std::vector excludeVolumes){ pIdHyp_ = pIdHyp; minNSegs_ = minNSegs; segLen_ = segLen; @@ -106,10 +113,11 @@ namespace trkf::sbn { pMax_ = pMax; pStep_ = pStep; angResol_ = angResol; - excludeRegions_ = excludeRegions; + fiducialVolumeInsets_ = fiducialVolumeInsets; + excludeVolumes_ = excludeVolumes; } explicit TrajectoryMCSFitter(const Parameters & p) - : TrajectoryMCSFitter(p().pIdHypothesis(),p().minNumSegments(),p().segmentLength(),p().minHitsPerSegment(),p().nElossSteps(),p().eLossMode(),p().pMin(),p().pMax(),p().pStep(),p().angResol(),p().excludeRegions()) {} + : TrajectoryMCSFitter(p().pIdHypothesis(),p().minNumSegments(),p().segmentLength(),p().minHitsPerSegment(),p().nElossSteps(),p().eLossMode(),p().pMin(),p().pMax(),p().pStep(),p().angResol(),p().fiducialVolumeInsets(),p().excludeVolumes()) {} // recob::MCSFitResult fitMcs(const recob::TrackTrajectory& traj, bool momDepConst = true) const { return fitMcs(traj,pIdHyp_,momDepConst); } recob::MCSFitResult fitMcs(const recob::Track& track, bool momDepConst = true) const { return fitMcs(track,pIdHyp_,momDepConst); } @@ -152,7 +160,9 @@ namespace trkf::sbn { double energyLossLandau(const double mass2,const double E2, const double x) const; // double GetE(const double initial_E, const double length_travelled, const double mass) const; - bool isInGoodRegion(const bool sizeGood, const double x, const double y, const double z) const; + std::vector setFiducialVolumes() const; + std::vector setExcludeVolumes() const; + bool isInVolume(const std::vector &volumes, const geo::Point_t &point) const; // private: int pIdHyp_; @@ -165,7 +175,8 @@ namespace trkf::sbn { double pMax_; double pStep_; double angResol_; - std::vector excludeRegions_; + std::vector fiducialVolumeInsets_; + std::vector excludeVolumes_; }; } diff --git a/sbncode/LArRecoProducer/mcsproducer.fcl b/sbncode/LArRecoProducer/mcsproducer.fcl index 31ca1df13..5db98243c 100644 --- a/sbncode/LArRecoProducer/mcsproducer.fcl +++ b/sbncode/LArRecoProducer/mcsproducer.fcl @@ -2,7 +2,8 @@ BEGIN_PROLOG mcs_sbn: { module_type: MCSFitAllPID MCS: { - excludeRegions: [] + fiducialVolumeInsets: [] + excludeVolumes: [] } TrackLabel: pandoraTrack } diff --git a/sbncode/LArRecoProducer/mcsproducer_icarus.fcl b/sbncode/LArRecoProducer/mcsproducer_icarus.fcl index 09e04041b..4926fdfc7 100644 --- a/sbncode/LArRecoProducer/mcsproducer_icarus.fcl +++ b/sbncode/LArRecoProducer/mcsproducer_icarus.fcl @@ -4,7 +4,8 @@ mcs_sbn: { MCS: { eLossMode: 2 angResol: 10.0 - excludeRegions: [190.0, 240.0, -9999.0, 9999.0, -9999.0, 9999.0] + fiducialVolumeInsets: [10.0, 10.0, 10.0, 10.0, 10.0, 10.0] + excludeVolumes: [190.0, 240.0, -9999.0, 9999.0, -9999.0, 9999.0] } TrackLabel: pandoraTrack } From 2dbf3b35546c56607729c13fbc2219688de43c27 Mon Sep 17 00:00:00 2001 From: gputnam Date: Fri, 28 Jul 2023 15:50:15 -0500 Subject: [PATCH 16/28] Allow time shift in trigger to IFBeam timing. Skip first event in stream in each run. --- sbncode/CAFMaker/FillTrigger.cxx | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/sbncode/CAFMaker/FillTrigger.cxx b/sbncode/CAFMaker/FillTrigger.cxx index fe68e00a1..f6bb51317 100644 --- a/sbncode/CAFMaker/FillTrigger.cxx +++ b/sbncode/CAFMaker/FillTrigger.cxx @@ -13,6 +13,15 @@ namespace caf triggerInfo.global_trigger_det_time = trig.TriggerTime(); double diff_ts = triggerInfo.global_trigger_det_time - triggerInfo.beam_gate_det_time; triggerInfo.trigger_within_gate = diff_ts; + + triggerInfo.prev_global_trigger_time = addltrig_info.previousTriggerTimestamp; + triggerInfo.source_type = (int)addltrig_info.sourceType; + triggerInfo.trigger_type = (int)addltrig_info.triggerType; + triggerInfo.trigger_id = addltrig_info.triggerID; + triggerInfo.gate_id = addltrig_info.gateID; + triggerInfo.trigger_count = addltrig_info.triggerCount; + triggerInfo.gate_count = addltrig_info.gateCount; + triggerInfo.gate_delta = (int)addltrig_info.gateCountFromPreviousTrigger; } void FillTriggerMC(double absolute_time, caf::SRTrigger& triggerInfo) { From 79920dd26bdcc65787a39a30855787bb09ba52bb Mon Sep 17 00:00:00 2001 From: gputnam Date: Sat, 21 Oct 2023 20:50:19 -0500 Subject: [PATCH 17/28] Update sbndata --- ups/product_deps | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ups/product_deps b/ups/product_deps index 8963bf559..2c9bf7711 100644 --- a/ups/product_deps +++ b/ups/product_deps @@ -257,7 +257,7 @@ larcv2 v2_1_0 - larsoft v09_72_00 - sbnanaobj v09_20_06_03 - sbndaq_artdaq_core v1_06_00of0 - -sbndata v01_04 - +sbndata v01_05 - sbnobj v09_16_00_02 - systematicstools v01_02_00 - nusystematics v01_02_10 - From 782e10c11955f611bfabf76a6f5c9babd74eeb6e Mon Sep 17 00:00:00 2001 From: gputnam Date: Sat, 21 Oct 2023 20:50:57 -0500 Subject: [PATCH 18/28] Include Geant4 libraries for buiding --- sbncode/SBNEventWeight/CMakeLists.txt | 3 +++ sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt | 1 + 2 files changed, 4 insertions(+) diff --git a/sbncode/SBNEventWeight/CMakeLists.txt b/sbncode/SBNEventWeight/CMakeLists.txt index de59387f6..9f4018c28 100644 --- a/sbncode/SBNEventWeight/CMakeLists.txt +++ b/sbncode/SBNEventWeight/CMakeLists.txt @@ -1,3 +1,6 @@ +include_directories ( $ENV{GEANT4REWEIGHT_INC} ) +link_directories ( $ENV{GEANT4REWEIGHT_LIB} ) + add_subdirectory(Base) add_subdirectory(Calculators) add_subdirectory(App) diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt b/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt index f05f8d274..568ecc7f6 100644 --- a/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt +++ b/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt @@ -1,5 +1,6 @@ include_directories ( $ENV{GEANT4_FQ_DIR}/include ) include_directories ( $ENV{GEANT4REWEIGHT_INC} ) +link_directories ( $ENV{GEANT4REWEIGHT_LIB} ) art_make_library( LIBRARY_NAME sbncode_SBNEventWeight_Calculators_Geant4 From 2591aa1cef7d6a57c9f7c3f0c5ec6dfc7eaf5edd Mon Sep 17 00:00:00 2001 From: gputnam Date: Sat, 21 Oct 2023 20:51:49 -0500 Subject: [PATCH 19/28] Fix ALP mass suppression. --- .../MeVPrtl/Tools/ALP/Meson2ALP_tool.cc | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/sbncode/EventGenerator/MeVPrtl/Tools/ALP/Meson2ALP_tool.cc b/sbncode/EventGenerator/MeVPrtl/Tools/ALP/Meson2ALP_tool.cc index fad80f2a1..2f709e026 100644 --- a/sbncode/EventGenerator/MeVPrtl/Tools/ALP/Meson2ALP_tool.cc +++ b/sbncode/EventGenerator/MeVPrtl/Tools/ALP/Meson2ALP_tool.cc @@ -76,6 +76,7 @@ class Meson2ALP : public IMeVPrtlFlux double fcG; //!< Axion coupling to gluon double fcB; //!< Axion coupling to U(1) B boson (before EW symmetry breaking) double fcW; //!< Axion coupling to SU(2) W bosons (before EW symmetry breaking) + bool fIncludeMassSuppressionFactor; // Whether to include a phenomonological suppression to mixing at high mass // branching ratios double fPi0BR; @@ -115,6 +116,7 @@ void Meson2ALP::configure(fhicl::ParameterSet const &pset) fMaxWeightPi0 = pset.get("MaxWeightPi0", 1.); fMaxWeightEta = pset.get("MaxWeightEta", 1.); fMaxWeightEtaP = pset.get("MaxWeightEtaP", 1.); + fIncludeMassSuppressionFactor = pset.get("IncludeMassSuppressionFactor", false); fPi0BR = Pi0BR(); fEtaBR = EtaBR(); @@ -136,7 +138,7 @@ double Meson2ALP::EtaBR() const { double mixing_angle = (1. / sqrt(6)) * (fcG * fpion / ffa) * (fM*fM - (4./9.)*pzero_mass*pzero_mass) / (fM*fM - eta_mass*eta_mass); double qcd_rate_f = 1.; - if (fM > eta_mass) qcd_rate_f = pow(fM/eta_mass, -1.6); + if (fIncludeMassSuppressionFactor && fM > eta_mass) qcd_rate_f = pow(fM/eta_mass, -1.6); return mixing_angle*mixing_angle*qcd_rate_f; } @@ -148,7 +150,7 @@ double Meson2ALP::EtaPBR() const { double mixing_angle = (1. / sqrt(12)) * (fcG * fpion / ffa) * (fM*fM - (16./9.)*pzero_mass*pzero_mass) / (fM*fM - etap_mass*etap_mass); double qcd_rate_f = 1.; - if (fM > etap_mass) qcd_rate_f = pow(fM/etap_mass, -1.6); + if (fIncludeMassSuppressionFactor && fM > etap_mass) qcd_rate_f = pow(fM/etap_mass, -1.6); return mixing_angle*mixing_angle*qcd_rate_f; } @@ -160,7 +162,7 @@ double Meson2ALP::Pi0BR() const { double mixing_angle = (1./6) * (fcG * fpion / ffa) * fM*fM / (fM*fM - pzero_mass*pzero_mass); double qcd_rate_f = 1.; - if (fM > pzero_mass) qcd_rate_f = pow(fM/pzero_mass, -1.6); + if (fIncludeMassSuppressionFactor && fM > pzero_mass) qcd_rate_f = pow(fM/pzero_mass, -1.6); return mixing_angle*mixing_angle*qcd_rate_f; } @@ -177,6 +179,10 @@ bool Meson2ALP::MakeFlux(const simb::MCFlux &flux, evgen::ldm::MeVPrtlFlux &alp, // Energy is same as for meson (don't worry about momentum conservation) double alp_energy = meson.mom.E(); + + // ignore if alp isn't energetic enough + if (alp_energy < fM) return false; + double alp_momentum = sqrt(alp_energy*alp_energy - fM*fM); // Momentum 4-vector From 79e67c8b77e87371873cbde853f14f71001b2ffd Mon Sep 17 00:00:00 2001 From: Bruce Howard Date: Thu, 25 Apr 2024 00:26:02 -0500 Subject: [PATCH 20/28] Revert "Include Geant4 libraries for buiding" This reverts commit 782e10c11955f611bfabf76a6f5c9babd74eeb6e. --- sbncode/SBNEventWeight/CMakeLists.txt | 3 --- sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt | 1 - 2 files changed, 4 deletions(-) diff --git a/sbncode/SBNEventWeight/CMakeLists.txt b/sbncode/SBNEventWeight/CMakeLists.txt index 9f4018c28..de59387f6 100644 --- a/sbncode/SBNEventWeight/CMakeLists.txt +++ b/sbncode/SBNEventWeight/CMakeLists.txt @@ -1,6 +1,3 @@ -include_directories ( $ENV{GEANT4REWEIGHT_INC} ) -link_directories ( $ENV{GEANT4REWEIGHT_LIB} ) - add_subdirectory(Base) add_subdirectory(Calculators) add_subdirectory(App) diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt b/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt index 568ecc7f6..f05f8d274 100644 --- a/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt +++ b/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt @@ -1,6 +1,5 @@ include_directories ( $ENV{GEANT4_FQ_DIR}/include ) include_directories ( $ENV{GEANT4REWEIGHT_INC} ) -link_directories ( $ENV{GEANT4REWEIGHT_LIB} ) art_make_library( LIBRARY_NAME sbncode_SBNEventWeight_Calculators_Geant4 From 0a96a7e416bfb6dd2d11938a7a32b91a6c771cea Mon Sep 17 00:00:00 2001 From: Bruce Howard Date: Thu, 25 Apr 2024 00:29:20 -0500 Subject: [PATCH 21/28] Revert "g4reweight bugfix and temporary parameter storage" This reverts commit ce8b26c71caf8eaf65c086afa73bac7ca8173f38. --- fcl/caf/cafmaker_common_defs.fcl | 2 -- .../SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx | 7 ++----- ...ntweight_geant4_sbn.fcl => eventweight_g4rwt_reint.fcl} | 6 +----- 3 files changed, 3 insertions(+), 12 deletions(-) rename sbncode/SBNEventWeight/jobs/geant4/{eventweight_geant4_sbn.fcl => eventweight_g4rwt_reint.fcl} (98%) diff --git a/fcl/caf/cafmaker_common_defs.fcl b/fcl/caf/cafmaker_common_defs.fcl index b4d45abbf..66fcdf7ab 100644 --- a/fcl/caf/cafmaker_common_defs.fcl +++ b/fcl/caf/cafmaker_common_defs.fcl @@ -1,4 +1,3 @@ -#include "eventweight_geant4_sbn.fcl" #include "eventweight_genie_sbn.fcl" #include "eventweight_flux_sbn.fcl" @@ -16,7 +15,6 @@ 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 a0cead5d2..a84ca7d5f 100644 --- a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx +++ b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx @@ -165,8 +165,6 @@ 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"); @@ -178,8 +176,6 @@ void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, FitParSigmas.push_back(theSigma); theNominals[theName] = theNominal; - - fParameterSet.AddParameter(theName, theSigma); } if (mode=="pm1sigma"){ @@ -244,7 +240,7 @@ std::vector Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { // Initialize weight vector for this MCTruth weight.clear(); - weight.resize(fNsims, 1.0); + weight.resize(1.0); // Loop over MCParticles in the event auto const& mcparticles = truthParticles.at(itruth); @@ -435,6 +431,7 @@ 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_geant4_sbn.fcl b/sbncode/SBNEventWeight/jobs/geant4/eventweight_g4rwt_reint.fcl similarity index 98% rename from sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl rename to sbncode/SBNEventWeight/jobs/geant4/eventweight_g4rwt_reint.fcl index e45065880..d08990c2c 100644 --- a/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl +++ b/sbncode/SBNEventWeight/jobs/geant4/eventweight_g4rwt_reint.fcl @@ -15,9 +15,7 @@ #include "piminus_reweight_parameters.fcl" #include "proton_reweight_parameters.fcl" -BEGIN_PROLOG - -sbn_eventweight_geant4: { +geant4weight_sbn: { module_type: "SBNEventWeight" AllowMissingTruth: true min_weight: 0.0 @@ -99,5 +97,3 @@ sbn_eventweight_geant4: { # variation) ### -END_PROLOG - From b4a7daa3929b9affb1a6b2e0e89492ec1701ac0f Mon Sep 17 00:00:00 2001 From: Bruce Howard Date: Thu, 25 Apr 2024 00:30:54 -0500 Subject: [PATCH 22/28] Revert "geant4reweight interface updates + generic sbn fcls" This reverts commit 38da5270785337518dadbe9bcbcd23d00d3966a5. --- 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, 9 insertions(+), 273 deletions(-) delete mode 100644 sbncode/SBNEventWeight/jobs/geant4/CMakeLists.txt delete mode 100644 sbncode/SBNEventWeight/jobs/geant4/eventweight_g4rwt_reint.fcl delete mode 100644 sbncode/SBNEventWeight/jobs/geant4/piminus_reweight_parameters.fcl delete mode 100644 sbncode/SBNEventWeight/jobs/geant4/piplus_reweight_parameters.fcl delete 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 335a19557..2989659b7 100644 --- a/sbncode/SBNEventWeight/App/CMakeLists.txt +++ b/sbncode/SBNEventWeight/App/CMakeLists.txt @@ -1,11 +1,8 @@ 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 - sbncode_SBNEventWeight_Calculators_Geant4 - nusimdata::SimulationBase + sbncode_SBNEventWeight_Calculators_BNBFlux 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 a84ca7d5f..aaf13c5f3 100644 --- a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx +++ b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx @@ -93,8 +93,11 @@ class Geant4WeightCalc : public WeightCalc { DECLARE_WEIGHTCALC(Geant4WeightCalc) }; +} // namespace evwgh + +} // namespace sbn -void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, +void sbn::evwgh::Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, CLHEP::HepRandomEngine& engine) { std::cout << "Using Geant4WeightCalc for reinteraction weights" << std::endl; @@ -222,7 +225,8 @@ void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, } -std::vector Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { +std::vector +sbn::evwgh::Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { // Get event/run/subrun numbers for output run_num = e.run(); @@ -482,8 +486,4 @@ return weight; } -REGISTER_WEIGHTCALC(Geant4WeightCalc) - -} // namespace evwgh - -} // namespace sbn +REGISTER_WEIGHTCALC(sbn::evwgh::Geant4WeightCalc) diff --git a/sbncode/SBNEventWeight/jobs/CMakeLists.txt b/sbncode/SBNEventWeight/jobs/CMakeLists.txt index 7897943c9..441c7dbde 100644 --- a/sbncode/SBNEventWeight/jobs/CMakeLists.txt +++ b/sbncode/SBNEventWeight/jobs/CMakeLists.txt @@ -4,5 +4,4 @@ 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 deleted file mode 100644 index 589bc695a..000000000 --- a/sbncode/SBNEventWeight/jobs/geant4/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -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 deleted file mode 100644 index d08990c2c..000000000 --- a/sbncode/SBNEventWeight/jobs/geant4/eventweight_g4rwt_reint.fcl +++ /dev/null @@ -1,99 +0,0 @@ -########################################################## -## 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 deleted file mode 100644 index b64756415..000000000 --- a/sbncode/SBNEventWeight/jobs/geant4/piminus_reweight_parameters.fcl +++ /dev/null @@ -1,66 +0,0 @@ -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 deleted file mode 100644 index c9cfc3e23..000000000 --- a/sbncode/SBNEventWeight/jobs/geant4/piplus_reweight_parameters.fcl +++ /dev/null @@ -1,66 +0,0 @@ -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 deleted file mode 100644 index 13113db75..000000000 --- a/sbncode/SBNEventWeight/jobs/geant4/proton_reweight_parameters.fcl +++ /dev/null @@ -1,23 +0,0 @@ -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 caf0d386adac651411f1232c09bf32975bd4ef36 Mon Sep 17 00:00:00 2001 From: Bruce Howard Date: Thu, 25 Apr 2024 00:34:14 -0500 Subject: [PATCH 23/28] Revert "Some tweaks to adjust for larsim EventWeight (used by uB) vs." This reverts commit 1f3c410ebcfff03bbddec876da215b724f5cb4dd. --- .../Calculators/Geant4/CMakeLists.txt | 1 - .../Calculators/Geant4/Geant4WeightCalc.cxx | 473 ++++++++++-------- ups/product_deps | 8 +- 3 files changed, 258 insertions(+), 224 deletions(-) diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt b/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt index f05f8d274..e1289208d 100644 --- a/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt +++ b/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt @@ -14,7 +14,6 @@ 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 aaf13c5f3..f4ba26774 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/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" +#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" @@ -52,7 +52,7 @@ class Geant4WeightCalc : public WeightCalc { void Configure(fhicl::ParameterSet const& p, CLHEP::HepRandomEngine& engine); - std::vector GetWeight(art::Event& e, size_t inu); + std::vector > GetWeight(art::Event& e); private: std::string fMCParticleProducer; //!< Label for the MCParticle producer @@ -87,18 +87,17 @@ 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 sbn::evwgh::Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, - CLHEP::HepRandomEngine& engine) +void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, + CLHEP::HepRandomEngine& engine) { std::cout << "Using Geant4WeightCalc for reinteraction weights" << std::endl; @@ -158,6 +157,9 @@ void sbn::evwgh::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); + } @@ -225,8 +227,8 @@ 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) { // Get event/run/subrun numbers for output run_num = e.run(); @@ -240,250 +242,283 @@ sbn::evwgh::Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { assert(truthParticles.isValid()); // Initialize the vector of event weights - std::vector weight; + 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()); - // Initialize weight vector for this MCTruth - weight.clear(); - weight.resize(1.0); + // 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; - } + e_inel_weight[itruth].clear(); + e_inel_weight[itruth].resize(fNsims,1.0); + e_elastic_weight[itruth].clear(); + e_elastic_weight[itruth].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 ); } - } // end loop over daughters - // --- Now that we have all the information about the track we need, here comes the reweighting part! --- // + 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; + } - //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; + } // 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 - size_t nSteps = trajpoint_PX.size(); + // --- Now that we have all the information about the track we need, here comes the reweighting part! --- // - if( nSteps < 2 ) continue; + //Make a G4ReweightTraj -- This is the reweightable object + G4ReweightTraj theTraj(mcpID, p_PDG, 0, event_num, std::make_pair(0,0)); - p_nElasticScatters = elastic_indices.size(); - for( size_t istep = 1; istep < nSteps; ++istep ){ + //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; - // 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; + size_t nSteps = trajpoint_PX.size(); - 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"; - } + if( nSteps < 2 ) continue; + p_nElasticScatters = elastic_indices.size(); + for( size_t istep = 1; istep < nSteps; ++istep ){ - 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) ); + // 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; - double len = sqrt( - std::pow( deltaX, 2 ) + - std::pow( deltaY, 2 ) + - std::pow( deltaZ, 2 ) - ); + 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 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 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) ); - if( istep == 1 ){ - theTraj.SetEnergy( sqrt( preStepP[0]*preStepP[0] + preStepP[1]*preStepP[1] + preStepP[2]*preStepP[2] + mass*mass ) ); - } + double len = sqrt( + std::pow( deltaX, 2 ) + + std::pow( deltaY, 2 ) + + std::pow( deltaZ, 2 ) + ); - G4ReweightStep * theStep = new G4ReweightStep( mcpID, p_PDG, 0, event_num, preStepP, postStepP, len, proc ); - theStep->SetDeltaX( deltaX ); - theStep->SetDeltaY( deltaY ); - theStep->SetDeltaZ( deltaZ ); + double preStepP[3] = { + trajpoint_PX.at(istep-1)*1.e3, + trajpoint_PY.at(istep-1)*1.e3, + trajpoint_PZ.at(istep-1)*1.e3 + }; - theTraj.AddStep( theStep ); + double postStepP[3] = { + trajpoint_PX.at(istep)*1.e3, + trajpoint_PY.at(istep)*1.e3, + trajpoint_PZ.at(istep)*1.e3 + }; - 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(); + if( istep == 1 ){ + theTraj.SetEnergy( sqrt( preStepP[0]*preStepP[0] + preStepP[1]*preStepP[1] + preStepP[2]*preStepP[2] + mass*mass ) ); + } - 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 ) - ); + G4ReweightStep * theStep = new G4ReweightStep( mcpID, p_PDG, 0, event_num, preStepP, postStepP, len, proc ); + theStep->SetDeltaX( deltaX ); + theStep->SetDeltaY( deltaY ); + theStep->SetDeltaZ( deltaZ ); - 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); + theTraj.AddStep( theStep ); - 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 ); - } + 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 ) + ); - 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 ); - } + 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); - // Loop through universes (j) - for (size_t j=0; jSetNewVals(UniverseVals.at(j)); - theReweighter->SetNewHists(ParMaker->GetFSHists()); - theReweighter->SetNewElasticHists(ParMaker->GetElasticHist()); + 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 ); + } - //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) + 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; + } - // 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; + } // loop through universes (j) 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(); + } // if ( ( TMath::Abs(p_PDG) == 211 || p_PDG == 2212 ) ) + if (fMakeOutputTrees) fOutTree_Particle->Fill(); + } // loop over mcparticles (i) + if (fMakeOutputTrees) fOutTree_MCTruth->Fill(); + } // loop over sets of MCtruth-associated particles (itruth) return weight; } -REGISTER_WEIGHTCALC(sbn::evwgh::Geant4WeightCalc) +REGISTER_WEIGHTCALC(Geant4WeightCalc) + +} // namespace evwgh + +} // namespace sbn diff --git a/ups/product_deps b/ups/product_deps index a23a4c052..0ba2b657f 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 @@ -259,9 +259,9 @@ sbnanaobj v09_20_06_03 - sbndaq_artdaq_core v1_06_00of0 - sbndata v01_05 - sbnobj v09_16_00_02 - -systematicstools v01_02_01 - -nusystematics v01_02_11 - -geant4reweight v01_00_03e - +systematicstools v01_02_00 - +nusystematics v01_02_10 - +geant4reweight v01_14_04 - cetmodules v3_20_00 - only_for_build end_product_list #################################### From e94da95d71e14f551ed5445b1f62597101800f50 Mon Sep 17 00:00:00 2001 From: Bruce Howard Date: Thu, 25 Apr 2024 00:34:49 -0500 Subject: [PATCH 24/28] Revert "Add the MicroBooNE weight calculator that interfaces with geant4reweight" This reverts commit 658fb3bbfe467951db12dc26e1254020d77b42c5. --- 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, 6 insertions(+), 699 deletions(-) delete mode 100644 sbncode/SBNEventWeight/Calculators/Geant4/BetheBlochForG4ReweightValid.h delete mode 100644 sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt delete mode 100644 sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx diff --git a/CMakeLists.txt b/CMakeLists.txt index 0b969ec7f..183b0b83b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -77,7 +77,6 @@ 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 383d41a3a..e2080685f 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 deleted file mode 100644 index 58cfdc99e..000000000 --- a/sbncode/SBNEventWeight/Calculators/Geant4/BetheBlochForG4ReweightValid.h +++ /dev/null @@ -1,145 +0,0 @@ - - -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 deleted file mode 100644 index e1289208d..000000000 --- a/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt +++ /dev/null @@ -1,22 +0,0 @@ -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 deleted file mode 100644 index f4ba26774..000000000 --- a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx +++ /dev/null @@ -1,524 +0,0 @@ -/** - * \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 0ba2b657f..4627e291e 100644 --- a/ups/product_deps +++ b/ups/product_deps @@ -261,7 +261,6 @@ sbndata v01_05 - sbnobj v09_16_00_02 - systematicstools v01_02_00 - nusystematics v01_02_10 - -geant4reweight v01_14_04 - cetmodules v3_20_00 - only_for_build end_product_list #################################### @@ -318,11 +317,11 @@ 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:s117:debug AR2320i00000:e1000:k250 -nq- c7:p392:debug c7:debug c7:debug c7:s117:debug -nq- -c7:prof c7:prof c7:prof c7:prof c7:s117:prof AR2320i00000:e1000:k250 -nq- c7:p392:prof c7:prof c7:prof c7:s117:prof -nq- -e20:debug e20:debug e20:debug e20:debug e20:s117:debug AR2320i00000:e1000:k250 -nq- e20:p392:debug e20:debug e20:debug e20:s117:debug -nq- -e20:prof e20:prof e20:prof e20:prof e20:s117:prof AR2320i00000:e1000:k250 -nq- e20:p392:prof e20:prof e20:prof e20:s117:prof -nq- +qualifier larsoft sbnobj sbnanaobj sbndaq_artdaq_core genie_xsec sbndata larcv2 systematicstools nusystematics notes +c7:debug c7:debug c7:debug c7:debug c7:s117:debug AR2320i00000:e1000:k250 -nq- c7:p392:debug c7:debug c7:debug -nq- +c7:prof c7:prof c7:prof c7:prof c7:s117:prof AR2320i00000:e1000:k250 -nq- c7:p392:prof c7:prof c7:prof -nq- +e20:debug e20:debug e20:debug e20:debug e20:s117:debug AR2320i00000:e1000:k250 -nq- e20:p392:debug e20:debug e20:debug -nq- +e20:prof e20:prof e20:prof e20:prof e20:s117:prof AR2320i00000:e1000:k250 -nq- e20:p392:prof e20:prof e20:prof -nq- end_qualifier_list #################################### From 2b1442db491cfd688f74d181b52dc53da6f5dccb Mon Sep 17 00:00:00 2001 From: Bruce Howard <61521590+brucehoward-physics@users.noreply.github.com> Date: Thu, 25 Apr 2024 10:17:15 -0500 Subject: [PATCH 25/28] Update product_deps I think I should have left these 1 version number higher but pushed them back in working to revert things. Putting this back where I think it should be. --- ups/product_deps | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ups/product_deps b/ups/product_deps index 4627e291e..b6cd7bd31 100644 --- a/ups/product_deps +++ b/ups/product_deps @@ -259,8 +259,8 @@ sbnanaobj v09_20_06_03 - sbndaq_artdaq_core v1_06_00of0 - sbndata v01_05 - sbnobj v09_16_00_02 - -systematicstools v01_02_00 - -nusystematics v01_02_10 - +systematicstools v01_02_01 - +nusystematics v01_02_11 - cetmodules v3_20_00 - only_for_build end_product_list #################################### From 4ac9e90633fae83603073de0ab3a2bcc585bc06a Mon Sep 17 00:00:00 2001 From: Bruce Howard Date: Thu, 25 Apr 2024 14:36:22 -0500 Subject: [PATCH 26/28] Revert "For cosmics that are crossing the cathode we added T0 information in calorimetry variables" This reverts commit 2eebbe0ebd6da1aab6c88395192aa7ed22baced7. --- sbncode/CAFMaker/CAFMaker_module.cc | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index d45fa0a7e..8a44dfc2c 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/sbncode/CAFMaker/CAFMaker_module.cc @@ -1828,10 +1828,9 @@ void CAFMaker::produce(art::Event& evt) noexcept { if (fmTrackDazzle.isValid() && fmTrackDazzle.at(iPart).size()==1) { FillTrackDazzle(fmTrackDazzle.at(iPart).front(), trk); } - std::cout<< "t0 calo: " << pfp.t0 <(), dprop, trk); } From 54dff1edd22615b40e118f73f68914790be85d5b Mon Sep 17 00:00:00 2001 From: Mun Jung Jung Date: Tue, 7 May 2024 12:42:16 -0500 Subject: [PATCH 27/28] last valid point for last segment validity --- sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx b/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx index cd0b4c53f..153189a29 100644 --- a/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx +++ b/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx @@ -106,7 +106,7 @@ void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj } //then add last segment if (thislen>0.) { - auto endpointpos = traj.LocationAtPoint(nextValid); + auto endpointpos = traj.LocationAtPoint(traj.LastValidPoint()); bool endpointposgood = isInVolume(fiducialVolumes, endpointpos) && !isInVolume(excludeVolumes, endpointpos); breakpoints.push_back(traj.LastValidPoint()+1); breakpointsgood.push_back(endpointposgood); From ad8a28f6505dbe086e3d3f66c92d88b64ec2da70 Mon Sep 17 00:00:00 2001 From: Bruce Howard Date: Wed, 8 May 2024 10:15:09 -0500 Subject: [PATCH 28/28] After talking to Moon and Gianluca, adding Gianluca's suggestion to change default behavior around fiducializing in TrajectoryMCSFitter --- .../LArReco/TrajectoryMCSFitter.cxx | 30 ++++++++++--------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx b/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx index 153189a29..caf10e17f 100644 --- a/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx +++ b/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx @@ -345,21 +345,23 @@ std::vector TrajectoryMCSFitter::setFiducialVolumes() const double zMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxZ() < rhs.MaxZ(); })->MaxZ(); ActiveVolumes.emplace_back(xMin, xMax, yMin, yMax, zMin, zMax); } - double fidInsetMinX = fiducialVolumeInsets_[0]; - double fidInsetMaxX = fiducialVolumeInsets_[1]; - double fidInsetMinY = fiducialVolumeInsets_[2]; - double fidInsetMaxY = fiducialVolumeInsets_[3]; - double fidInsetMinZ = fiducialVolumeInsets_[4]; - double fidInsetMaxZ = fiducialVolumeInsets_[5]; + double fidInsetMinX = 0; + double fidInsetMaxX = 0; + double fidInsetMinY = 0; + double fidInsetMaxY = 0; + double fidInsetMinZ = 0; + double fidInsetMaxZ = 0; - if (fiducialVolumeInsets_.size() != 6) { + if (fiducialVolumeInsets_.size() == 6) { + fidInsetMinX = fiducialVolumeInsets_[0]; + fidInsetMaxX = fiducialVolumeInsets_[1]; + fidInsetMinY = fiducialVolumeInsets_[2]; + fidInsetMaxY = fiducialVolumeInsets_[3]; + fidInsetMinZ = fiducialVolumeInsets_[4]; + fidInsetMaxZ = fiducialVolumeInsets_[5]; + } + else if (!fiducialVolumeInsets_.empty()) { std::cout << "Error: fiducialVolumeInsets vector must have length of 6, not fiducializing" << std::endl; - fidInsetMinX = 0.0; - fidInsetMaxX = 0.0; - fidInsetMinY = 0.0; - fidInsetMaxY = 0.0; - fidInsetMinZ = 0.0; - fidInsetMaxZ = 0.0; } for (const geo::BoxBoundedGeo &AV: ActiveVolumes) { fiducialVolumes.emplace_back(AV.MinX() + fidInsetMinX, AV.MaxX() - fidInsetMaxX, @@ -391,4 +393,4 @@ bool TrajectoryMCSFitter::isInVolume(const std::vector &volu } } return false; -} \ No newline at end of file +}