diff --git a/CMakeLists.txt b/CMakeLists.txt index bbbe20d40..0df1801dc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -77,6 +77,7 @@ find_package( CLHEP REQUIRED ) find_package( ROOT REQUIRED ) find_package( Geant4 REQUIRED ) find_package( Boost COMPONENTS system filesystem REQUIRED ) +find_package( geant4reweight REQUIRED ) # macros for dictionary and simple_plugin include(ArtDictionary) diff --git a/fcl/caf/cafmaker_common_defs.fcl b/fcl/caf/cafmaker_common_defs.fcl index 66fcdf7ab..b4d45abbf 100644 --- a/fcl/caf/cafmaker_common_defs.fcl +++ b/fcl/caf/cafmaker_common_defs.fcl @@ -1,3 +1,4 @@ +#include "eventweight_geant4_sbn.fcl" #include "eventweight_genie_sbn.fcl" #include "eventweight_flux_sbn.fcl" @@ -15,6 +16,7 @@ cafmaker_common_producers: { rns: { module_type: "RandomNumberSaver" } genieweight: @local::sbn_eventweight_genie fluxweight: @local::sbn_eventweight_flux + geant4weight: @local::sbn_eventweight_geant4 } END_PROLOG diff --git a/sbncode/BeamSpillInfoRetriever/EXTRetriever/EXTRetriever_module.cc b/sbncode/BeamSpillInfoRetriever/BNBEXTRetriever/BNBEXTRetriever_module.cc similarity index 79% rename from sbncode/BeamSpillInfoRetriever/EXTRetriever/EXTRetriever_module.cc rename to sbncode/BeamSpillInfoRetriever/BNBEXTRetriever/BNBEXTRetriever_module.cc index 77d0661e1..9f14e29d3 100644 --- a/sbncode/BeamSpillInfoRetriever/EXTRetriever/EXTRetriever_module.cc +++ b/sbncode/BeamSpillInfoRetriever/BNBEXTRetriever/BNBEXTRetriever_module.cc @@ -1,7 +1,7 @@ //////////////////////////////////////////////////////////////////////// -// Class: EXTRetriever +// Class: BNBEXTRetriever // Plugin Type: producer -// File: EXTRetriever_module.cc +// File: BNBEXTRetriever_module.cc // // Created by hand Thurs June 24th 2021 by J. Zennamo (FNAL) // @@ -22,7 +22,7 @@ #include "larcorealg/Geometry/Exceptions.h" #include "artdaq-core/Data/Fragment.hh" -#include "sbndaq-artdaq-core/Overlays/ICARUS/ICARUSTriggerUDPFragment.hh" +#include "sbndaq-artdaq-core/Overlays/ICARUS/ICARUSTriggerV3Fragment.hh" #include "lardataalg/DetectorInfo/DetectorPropertiesStandard.h" #include "lardata/DetectorInfoServices/DetectorPropertiesService.h" @@ -36,10 +36,10 @@ #include namespace sbn { - class EXTRetriever; + class BNBEXTRetriever; } -class sbn::EXTRetriever : public art::EDProducer { +class sbn::BNBEXTRetriever : public art::EDProducer { public: struct Config { @@ -57,15 +57,15 @@ class sbn::EXTRetriever : public art::EDProducer { using Parameters = art::EDProducer::Table; - explicit EXTRetriever(Parameters const& params); + explicit BNBEXTRetriever(Parameters const& params); // The compiler-generated destructor is fine for non-base // classes without bare pointers or other resource use. // Plugins should not be copied or assigned. - EXTRetriever(EXTRetriever const&) = delete; - EXTRetriever(EXTRetriever&&) = delete; - EXTRetriever& operator=(EXTRetriever const&) = delete; - EXTRetriever& operator=(EXTRetriever&&) = delete; + BNBEXTRetriever(BNBEXTRetriever const&) = delete; + BNBEXTRetriever(BNBEXTRetriever&&) = delete; + BNBEXTRetriever& operator=(BNBEXTRetriever const&) = delete; + BNBEXTRetriever& operator=(BNBEXTRetriever&&) = delete; // Required functions. void produce(art::Event& e) override; @@ -82,7 +82,7 @@ class sbn::EXTRetriever : public art::EDProducer { }; -sbn::EXTRetriever::EXTRetriever(Parameters const& params) +sbn::BNBEXTRetriever::BNBEXTRetriever(Parameters const& params) : EDProducer{params}, raw_data_label_(params().RawDataLabel()) { @@ -91,7 +91,7 @@ sbn::EXTRetriever::EXTRetriever(Parameters const& params) TotalEXTCounts = 0; } -void sbn::EXTRetriever::produce(art::Event& e) +void sbn::BNBEXTRetriever::produce(art::Event& e) { //Here we read in the artdaq Fragments and extract three pieces of information: @@ -100,19 +100,22 @@ void sbn::EXTRetriever::produce(art::Event& e) // 3. the number of beam spills since the previously triggered event, number_of_gates_since_previous_event int gate_type = 0; - auto const & raw_data = e.getProduct< std::vector >({ raw_data_label_, "ICARUSTriggerUDP" }); + auto const & raw_data = e.getProduct< std::vector >({ raw_data_label_, "ICARUSTriggerV3" }); unsigned int number_of_gates_since_previous_event = 0; for(auto raw_datum : raw_data){ - icarus::ICARUSTriggerUDPFragment frag(raw_datum); + icarus::ICARUSTriggerV3Fragment frag(raw_datum); std::string data = frag.GetDataString(); char *buffer = const_cast(data.c_str()); - icarus::ICARUSTriggerInfo datastream_info = icarus::parse_ICARUSTriggerString(buffer); + icarus::ICARUSTriggerInfo datastream_info = icarus::parse_ICARUSTriggerV3String(buffer); gate_type = datastream_info.gate_type; - number_of_gates_since_previous_event = frag.getDeltaGatesBNB(); + number_of_gates_since_previous_event = frag.getDeltaGatesBNBOffMaj(); + std::cout << "BNB OFF MAJ : " << frag.getDeltaGatesBNBOffMaj() << std::endl; + std::cout << "NuMI OFF MAJ : " << frag.getDeltaGatesNuMIOffMaj() << std::endl; + } //We only want to process EXT gates, i.e. type 3 @@ -135,13 +138,13 @@ void sbn::EXTRetriever::produce(art::Event& e) } //end loop over events -void sbn::EXTRetriever::beginSubRun(art::SubRun& sr) +void sbn::BNBEXTRetriever::beginSubRun(art::SubRun& sr) { return; } //____________________________________________________________________________ -void sbn::EXTRetriever::endSubRun(art::SubRun& sr) +void sbn::BNBEXTRetriever::endSubRun(art::SubRun& sr) { // We will add all of the EXTCountInfo data-products to the // art::SubRun so it persists @@ -152,4 +155,4 @@ void sbn::EXTRetriever::endSubRun(art::SubRun& sr) return; } -DEFINE_ART_MODULE(sbn::EXTRetriever) +DEFINE_ART_MODULE(sbn::BNBEXTRetriever) diff --git a/sbncode/BeamSpillInfoRetriever/EXTRetriever/CMakeLists.txt b/sbncode/BeamSpillInfoRetriever/BNBEXTRetriever/CMakeLists.txt similarity index 92% rename from sbncode/BeamSpillInfoRetriever/EXTRetriever/CMakeLists.txt rename to sbncode/BeamSpillInfoRetriever/BNBEXTRetriever/CMakeLists.txt index 316419a52..1c49f8646 100644 --- a/sbncode/BeamSpillInfoRetriever/EXTRetriever/CMakeLists.txt +++ b/sbncode/BeamSpillInfoRetriever/BNBEXTRetriever/CMakeLists.txt @@ -1,5 +1,5 @@ -cet_build_plugin(EXTRetriever art::module +cet_build_plugin(BNBEXTRetriever art::module LIBRARIES art::Persistency_Common art::Utilities canvas::canvas diff --git a/sbncode/BeamSpillInfoRetriever/BNBRetriever/BNBRetriever_module.cc b/sbncode/BeamSpillInfoRetriever/BNBRetriever/BNBRetriever_module.cc index a907341d5..db124cf36 100644 --- a/sbncode/BeamSpillInfoRetriever/BNBRetriever/BNBRetriever_module.cc +++ b/sbncode/BeamSpillInfoRetriever/BNBRetriever/BNBRetriever_module.cc @@ -20,7 +20,7 @@ #include "larcorealg/CoreUtils/counter.h" #include "artdaq-core/Data/Fragment.hh" -#include "sbndaq-artdaq-core/Overlays/ICARUS/ICARUSTriggerV2Fragment.hh" +#include "sbndaq-artdaq-core/Overlays/ICARUS/ICARUSTriggerV3Fragment.hh" #include "sbnobj/Common/POTAccounting/BNBSpillInfo.h" @@ -245,23 +245,23 @@ sbn::BNBRetriever::TriggerInfo_t sbn::BNBRetriever::extractTriggerInfo(art::Even // 2. the time of the previously triggered event, t_previous_event (NOTE: Events are non-sequential!) // 3. the number of beam spills since the previously triggered event, number_of_gates_since_previous_event - auto const & raw_data = e.getProduct< std::vector >({ raw_data_label, "ICARUSTriggerV2" }); + auto const & raw_data = e.getProduct< std::vector >({ raw_data_label, "ICARUSTriggerV3" }); TriggerInfo_t triggerInfo; for(auto raw_datum : raw_data){ uint64_t artdaq_ts = raw_datum.timestamp(); - icarus::ICARUSTriggerV2Fragment frag(raw_datum); + icarus::ICARUSTriggerV3Fragment frag(raw_datum); std::string data = frag.GetDataString(); char *buffer = const_cast(data.c_str()); - icarus::ICARUSTriggerInfo datastream_info = icarus::parse_ICARUSTriggerV2String(buffer); + icarus::ICARUSTriggerInfo datastream_info = icarus::parse_ICARUSTriggerV3String(buffer); triggerInfo.gate_type = datastream_info.gate_type; - triggerInfo.number_of_gates_since_previous_event = frag.getDeltaGatesBNB(); + triggerInfo.number_of_gates_since_previous_event = frag.getDeltaGatesBNBMaj(); triggerInfo.t_current_event = static_cast(artdaq_ts)/(1000000000.0); //check this offset... if(triggerInfo.gate_type == 1) - triggerInfo.t_previous_event = (static_cast(frag.getLastTimestampBNB()))/(1e9); + triggerInfo.t_previous_event = (static_cast(frag.getLastTimestampBNBMaj()))/(1e9); else triggerInfo.t_previous_event = (static_cast(frag.getLastTimestampOther()))/(1000000000.0); diff --git a/sbncode/BeamSpillInfoRetriever/CMakeLists.txt b/sbncode/BeamSpillInfoRetriever/CMakeLists.txt index 1a6fecf80..0caadf2a5 100644 --- a/sbncode/BeamSpillInfoRetriever/CMakeLists.txt +++ b/sbncode/BeamSpillInfoRetriever/CMakeLists.txt @@ -1,6 +1,7 @@ add_subdirectory(BNBRetriever) add_subdirectory(NuMIRetriever) -add_subdirectory(EXTRetriever) +add_subdirectory(BNBEXTRetriever) +add_subdirectory(NuMIEXTRetriever) add_subdirectory(job) install_headers() diff --git a/sbncode/BeamSpillInfoRetriever/NuMIEXTRetriever/CMakeLists.txt b/sbncode/BeamSpillInfoRetriever/NuMIEXTRetriever/CMakeLists.txt new file mode 100644 index 000000000..f11657b00 --- /dev/null +++ b/sbncode/BeamSpillInfoRetriever/NuMIEXTRetriever/CMakeLists.txt @@ -0,0 +1,20 @@ + +cet_build_plugin(NuMIEXTRetriever art::module + LIBRARIES + art::Persistency_Common + art::Utilities canvas::canvas + cetlib::cetlib cetlib_except::cetlib_except + messagefacility::MF_MessageLogger + sbnobj::Common_POTAccounting + sbndaq_artdaq_core::sbndaq-artdaq-core_Overlays_Common + sbndaq_artdaq_core::sbndaq-artdaq-core_Overlays + sbndaq_artdaq_core::sbndaq-artdaq-core_Overlays_ICARUS + artdaq_core::artdaq-core_Utilities + lardata::Utilities + larcore::Geometry_AuxDetGeometry_service +) + +install_headers() +install_fhicl() +install_source() + diff --git a/sbncode/BeamSpillInfoRetriever/NuMIEXTRetriever/NuMIEXTRetriever_module.cc b/sbncode/BeamSpillInfoRetriever/NuMIEXTRetriever/NuMIEXTRetriever_module.cc new file mode 100644 index 000000000..7fdc974e2 --- /dev/null +++ b/sbncode/BeamSpillInfoRetriever/NuMIEXTRetriever/NuMIEXTRetriever_module.cc @@ -0,0 +1,155 @@ +//////////////////////////////////////////////////////////////////////// +// Class: NuMIEXTRetriever +// Plugin Type: producer +// File: NuMIEXTRetriever_module.cc +// +// Created by hand Thurs June 24th 2021 by J. Zennamo (FNAL) +// +//////////////////////////////////////////////////////////////////////// + +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Principal/SubRun.h" +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/types/Atom.h" +#include "messagefacility/MessageLogger/MessageLogger.h" +#include "lardata/Utilities/AssociationUtil.h" +#include "lardataobj/Utilities/sparse_vector.h" +#include "larcoreobj/SimpleTypesAndConstants/RawTypes.h" +#include "larcorealg/Geometry/Exceptions.h" + +#include "artdaq-core/Data/Fragment.hh" +#include "sbndaq-artdaq-core/Overlays/ICARUS/ICARUSTriggerV3Fragment.hh" + +#include "lardataalg/DetectorInfo/DetectorPropertiesStandard.h" +#include "lardata/DetectorInfoServices/DetectorPropertiesService.h" +#include "sbnobj/Common/POTAccounting/EXTCountInfo.h" + +#include +#include +#include +#include +#include +#include + +namespace sbn { + class NuMIEXTRetriever; +} + +class sbn::NuMIEXTRetriever : public art::EDProducer { +public: + + struct Config { + + using Name = fhicl::Name; + using Comment = fhicl::Comment; + + fhicl::Atom RawDataLabel { + Name{ "raw_data_label" }, + Comment{ "art data product instance name for trigger information (product label is 'daq')" } + }; + + }; // Config + + using Parameters = art::EDProducer::Table; + + + explicit NuMIEXTRetriever(Parameters const& params); + // The compiler-generated destructor is fine for non-base + // classes without bare pointers or other resource use. + + // Plugins should not be copied or assigned. + NuMIEXTRetriever(NuMIEXTRetriever const&) = delete; + NuMIEXTRetriever(NuMIEXTRetriever&&) = delete; + NuMIEXTRetriever& operator=(NuMIEXTRetriever const&) = delete; + NuMIEXTRetriever& operator=(NuMIEXTRetriever&&) = delete; + + // Required functions. + void produce(art::Event& e) override; + void beginSubRun(art::SubRun& sr) override; + void endSubRun(art::SubRun& sr) override; + +private: + std::vector< sbn::EXTCountInfo > fOutExtInfos; + + // input labels + std::string raw_data_label_; + int TotalEXTCounts; + +}; + + +sbn::NuMIEXTRetriever::NuMIEXTRetriever(Parameters const& params) + : EDProducer{params}, + raw_data_label_(params().RawDataLabel()) +{ + + produces< std::vector< sbn::EXTCountInfo >, art::InSubRun >(); + TotalEXTCounts = 0; +} + +void sbn::NuMIEXTRetriever::produce(art::Event& e) +{ + + //Here we read in the artdaq Fragments and extract three pieces of information: + // 1. The time of the current event, t_current_event + // 2. the time of the previously triggered event, t_previous_event (NOTE: Events are non-sequential!) + // 3. the number of beam spills since the previously triggered event, number_of_gates_since_previous_event + + //int gate_type = 0; + auto const & raw_data = e.getProduct< std::vector >({ raw_data_label_, "ICARUSTriggerV3" }); + + unsigned int number_of_gates_since_previous_event = 0; + + for(auto raw_datum : raw_data){ + + icarus::ICARUSTriggerV3Fragment frag(raw_datum); + std::string data = frag.GetDataString(); + char *buffer = const_cast(data.c_str()); + icarus::ICARUSTriggerInfo datastream_info = icarus::parse_ICARUSTriggerV3String(buffer); + //gate_type = datastream_info.gate_type; + number_of_gates_since_previous_event = frag.getDeltaGatesNuMIOffMaj(); + + + } + + //We only want to process EXT gates, i.e. type 3 + + // Keep track of the number of beam gates the DAQ thinks + // are in this file + TotalEXTCounts += number_of_gates_since_previous_event; + + //Store everything in our data-product + sbn::EXTCountInfo extInfo; + extInfo.gates_since_last_trigger = number_of_gates_since_previous_event; + + fOutExtInfos.push_back(extInfo); + // We do not write these to the art::Events because + // we can filter events but want to keep all the POT + // information, so we'll write it to the SubRun + + +} //end loop over events + + +void sbn::NuMIEXTRetriever::beginSubRun(art::SubRun& sr) +{ + return; +} + +//____________________________________________________________________________ +void sbn::NuMIEXTRetriever::endSubRun(art::SubRun& sr) +{ + // We will add all of the EXTCountInfo data-products to the + // art::SubRun so it persists + auto p = std::make_unique< std::vector< sbn::EXTCountInfo > >(fOutExtInfos); + + sr.put(std::move(p)); + + return; +} + +DEFINE_ART_MODULE(sbn::NuMIEXTRetriever) diff --git a/sbncode/BeamSpillInfoRetriever/NuMIRetriever/NuMIRetriever_module.cc b/sbncode/BeamSpillInfoRetriever/NuMIRetriever/NuMIRetriever_module.cc index 84e4c6d03..a0eb1e48a 100644 --- a/sbncode/BeamSpillInfoRetriever/NuMIRetriever/NuMIRetriever_module.cc +++ b/sbncode/BeamSpillInfoRetriever/NuMIRetriever/NuMIRetriever_module.cc @@ -19,7 +19,7 @@ #include "larcorealg/Geometry/Exceptions.h" #include "artdaq-core/Data/Fragment.hh" -#include "sbndaq-artdaq-core/Overlays/ICARUS/ICARUSTriggerV2Fragment.hh" +#include "sbndaq-artdaq-core/Overlays/ICARUS/ICARUSTriggerV3Fragment.hh" #include "lardataalg/DetectorInfo/DetectorPropertiesStandard.h" #include "lardata/DetectorInfoServices/DetectorPropertiesService.h" @@ -104,9 +104,9 @@ void sbn::NuMIRetriever::produce(art::Event &e) // TODO: long-term goal -- can we fix this? if (e.event() == 1) return; - int gate_type = 0; + // int gate_type = 0; art::Handle< std::vector > raw_data_ptr; - e.getByLabel(raw_data_label_, "ICARUSTriggerV2", raw_data_ptr); + e.getByLabel(raw_data_label_, "ICARUSTriggerV3", raw_data_ptr); auto const & raw_data = (*raw_data_ptr); double t_current_event = 0; @@ -116,145 +116,143 @@ void sbn::NuMIRetriever::produce(art::Event &e) for(auto raw_datum : raw_data){ uint64_t artdaq_ts = raw_datum.timestamp(); - icarus::ICARUSTriggerV2Fragment frag(raw_datum); + icarus::ICARUSTriggerV3Fragment frag(raw_datum); std::string data = frag.GetDataString(); char *buffer = const_cast(data.c_str()); - icarus::ICARUSTriggerInfo datastream_info = icarus::parse_ICARUSTriggerV2String(buffer); - gate_type = datastream_info.gate_type; - number_of_gates_since_previous_event = frag.getDeltaGatesNuMI(); + icarus::ICARUSTriggerInfo datastream_info = icarus::parse_ICARUSTriggerV3String(buffer); + //gate_type = datastream_info.gate_type; + number_of_gates_since_previous_event = frag.getDeltaGatesNuMIMaj(); t_current_event = static_cast(artdaq_ts)/(1000000000.); //check this offset... - if(gate_type == 2) - t_previous_event = (static_cast(frag.getLastTimestampNuMI()))/(1000000000.); - else - t_previous_event = (static_cast(frag.getLastTimestampOther()))/(1000000000.); + + t_previous_event = (static_cast(frag.getLastTimestampNuMIMaj()))/(1000000000.); + } std::cout << std::setprecision(19) << "Previous : " << t_previous_event << ", Current : " << t_current_event << std::endl; //We only want to process NuMI gates, i.e. type 2 - if(gate_type == 2) - { - // Keep track of the number of beam gates the DAQ thinks - // are in this job - TotalBeamSpills += number_of_gates_since_previous_event; + + // Keep track of the number of beam gates the DAQ thinks + // are in this job + TotalBeamSpills += number_of_gates_since_previous_event; + + // Fill up the BFP cache with times starting at the previous event + // + // If the difference in time between events is bigger than the fcl-provided TimeWindow, then + // we won't be able to fill up the cache to be big enough. In that case, resize the window as + // necessary. + double this_window_size = t_current_event - t_previous_event + 2*fTimePad; + if (this_window_size > fTimeWindow) { + std::cout << "Resizing time window from: " << fTimeWindow << " to: " << this_window_size << std::endl; + fTimeWindow = this_window_size; + MakeBFP(); + } + + // DO NOT CHANGE THESE LINES OR THEIR ORDER + // If you really think you need to, please reach out to grayputnam uchicago.edu + bfp->FillCache(t_current_event + fTimePad); + bfp->FillCache(t_previous_event - fTimePad); + std::vector times_temps = bfp->GetTimeList(fDeviceUsedForTiming); + + int spill_count = 0; + // Iterating through each of the beamline times + for (size_t i = 0; i < times_temps.size(); i++) { - // Fill up the BFP cache with times starting at the previous event - // - // If the difference in time between events is bigger than the fcl-provided TimeWindow, then - // we won't be able to fill up the cache to be big enough. In that case, resize the window as - // necessary. - double this_window_size = t_current_event - t_previous_event + 2*fTimePad; - if (this_window_size > fTimeWindow) { - std::cout << "Resizing time window from: " << fTimeWindow << " to: " << this_window_size << std::endl; - fTimeWindow = this_window_size; - MakeBFP(); + // Only continue if these times are matched to our DAQ time + // plus or minus some time padding, currently using 3.3 ms + // which is half the Booster Rep Rate + if(e.event() != 1){//We already addressed the "first event" above + if(times_temps[i] > t_current_event){continue;} + if(times_temps[i] <= t_previous_event){continue;} } - - // DO NOT CHANGE THESE LINES OR THEIR ORDER - // If you really think you need to, please reach out to grayputnam uchicago.edu - bfp->FillCache(t_current_event + fTimePad); - bfp->FillCache(t_previous_event - fTimePad); - std::vector times_temps = bfp->GetTimeList(fDeviceUsedForTiming); - - int spill_count = 0; - // Iterating through each of the beamline times - for (size_t i = 0; i < times_temps.size(); i++) { - - // Only continue if these times are matched to our DAQ time - // plus or minus some time padding, currently using 3.3 ms - // which is half the Booster Rep Rate - if(e.event() != 1){//We already addressed the "first event" above - if(times_temps[i] > t_current_event){continue;} - if(times_temps[i] <= t_previous_event){continue;} - } - - //count found spills - spill_count++; - - //initialize all devices found in NuMISpillInfo.h in sbnobj - double HRNDIR = -1.; - double NSLINA = -1.; - double NSLINB = -1.; - double NSLINC = -1.; - double NSLIND = -1.; - double TOR101 = -1.; - double TORTGT = -1.; - double TR101D = -1.; - double TRTGTD = -1.; - std::vector< double > HP121; - std::vector< double > VP121; - std::vector< double > HPTGT; - std::vector< double > VPTGT; - std::vector< double > HITGT; - std::vector< double > VITGT; - std::vector< double > MTGTDS; - double TRTGTD_time = -1.; - std::cout << "Grabbing IFBeam info!" << std::endl; - try{bfp->GetNamedData(times_temps[i], "E:TRTGTD@",&TRTGTD,&TRTGTD_time);}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " << "got exception: " << we.what() << "\n";} - try{bfp->GetNamedData(times_temps[i], "E:TR101D",&TR101D);}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " << "got exception: " << we.what() << "\n";} - try{bfp->GetNamedData(times_temps[i], "E:HRNDIR",&HRNDIR);}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " << "got exception: " << we.what() << "\n";} - try{bfp->GetNamedData(times_temps[i], "E:NSLINA",&NSLINA);}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " << "got exception: " << we.what() << "\n";} - try{bfp->GetNamedData(times_temps[i], "E:NSLINB",&NSLINB);}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " << "got exception: " << we.what() << "\n";} - try{bfp->GetNamedData(times_temps[i], "E:NSLINC",&NSLINC);}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " << "got exception: " << we.what() << "\n";} - try{bfp->GetNamedData(times_temps[i], "E:NSLIND",&NSLIND);}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " << "got exception: " << we.what() << "\n";} - try{bfp->GetNamedData(times_temps[i], "E:TOR101",&TOR101);}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " << "got exception: " << we.what() << "\n";} - try{bfp->GetNamedData(times_temps[i], "E:TORTGT",&TORTGT);}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " << "got exception: " << we.what() << "\n";} - // BPM Positions and Intensities - they each have 7 elements - // First is an average value of a few batches (often 2,3,4) - // used for auto-tuning, so we should disregard it - - try{HP121 = bfp->GetNamedVector(times_temps[i], "E:HP121[]");}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " << "got exception: " << we.what() << "\n";} - try{VP121 = bfp->GetNamedVector(times_temps[i], "E:VP121[]");}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " <<"got exception: " << we.what() << "\n";} - try{HPTGT = bfp->GetNamedVector(times_temps[i], "E:HPTGT[]");}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " <<"got exception: " << we.what() << "\n";} - try{VPTGT = bfp->GetNamedVector(times_temps[i], "E:VPTGT[]");}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " <<"got exception: " << we.what() << "\n";} - try{HITGT = bfp->GetNamedVector(times_temps[i], "E:HITGT[]");}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " <<"got exception: " << we.what() << "\n";} - try{VITGT = bfp->GetNamedVector(times_temps[i], "E:VITGT[]");}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " <<"got exception: " << we.what() << "\n";} - try{MTGTDS = bfp->GetNamedVector(times_temps[i], "E:MTGTDS[]");}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " <<"got exception: " << we.what() << "\n";} - - std::cout << "Finished getting IFBeam info" << std::endl; - std::cout << "BFP Time: " << times_temps[i] << " TOROID Time: " << TRTGTD_time << " TOROID COUNT: " << TRTGTD << std::endl; - unsigned long int time_closest_int = (int) TRTGTD_time; - double time_closest_ns = (TRTGTD_time - time_closest_int)*1000000000; - - sbn::NuMISpillInfo NuMIbeamInfo; - NuMIbeamInfo.TORTGT = TORTGT*1e12; //include factor of 1e12 protons in POT calculation - NuMIbeamInfo.TOR101 = TOR101*1e12; //include factor of 1e12 protons in POT calculation - NuMIbeamInfo.TRTGTD = TRTGTD*1e12; //include factor of 1e12 protons in POT calculation - NuMIbeamInfo.TR101D = TR101D*1e12; //include factor of 1e12 protons in POT calculation - NuMIbeamInfo.HRNDIR = HRNDIR; - NuMIbeamInfo.NSLINA = NSLINA; - NuMIbeamInfo.NSLINB = NSLINB; - NuMIbeamInfo.NSLINC = NSLINC; - NuMIbeamInfo.NSLIND = NSLIND; - - NuMIbeamInfo.HP121 = HP121; - NuMIbeamInfo.VP121 = VP121; - NuMIbeamInfo.HPTGT = HPTGT; - NuMIbeamInfo.VPTGT = VPTGT; - NuMIbeamInfo.HITGT = HITGT; - NuMIbeamInfo.VITGT = VITGT; - NuMIbeamInfo.MTGTDS = MTGTDS; - - NuMIbeamInfo.time = times_temps[i]; - NuMIbeamInfo.event = e.event(); - NuMIbeamInfo.spill_time_s = time_closest_int; - NuMIbeamInfo.spill_time_ns = time_closest_ns; + + //count found spills + spill_count++; + + //initialize all devices found in NuMISpillInfo.h in sbnobj + double HRNDIR = -1.; + double NSLINA = -1.; + double NSLINB = -1.; + double NSLINC = -1.; + double NSLIND = -1.; + double TOR101 = -1.; + double TORTGT = -1.; + double TR101D = -1.; + double TRTGTD = -1.; + std::vector< double > HP121; + std::vector< double > VP121; + std::vector< double > HPTGT; + std::vector< double > VPTGT; + std::vector< double > HITGT; + std::vector< double > VITGT; + std::vector< double > MTGTDS; + double TRTGTD_time = -1.; + std::cout << "Grabbing IFBeam info!" << std::endl; + try{bfp->GetNamedData(times_temps[i], "E:TRTGTD@",&TRTGTD,&TRTGTD_time);}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " << "got exception: " << we.what() << "\n";} + try{bfp->GetNamedData(times_temps[i], "E:TR101D",&TR101D);}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " << "got exception: " << we.what() << "\n";} + try{bfp->GetNamedData(times_temps[i], "E:HRNDIR",&HRNDIR);}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " << "got exception: " << we.what() << "\n";} + try{bfp->GetNamedData(times_temps[i], "E:NSLINA",&NSLINA);}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " << "got exception: " << we.what() << "\n";} + try{bfp->GetNamedData(times_temps[i], "E:NSLINB",&NSLINB);}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " << "got exception: " << we.what() << "\n";} + try{bfp->GetNamedData(times_temps[i], "E:NSLINC",&NSLINC);}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " << "got exception: " << we.what() << "\n";} + try{bfp->GetNamedData(times_temps[i], "E:NSLIND",&NSLIND);}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " << "got exception: " << we.what() << "\n";} + try{bfp->GetNamedData(times_temps[i], "E:TOR101",&TOR101);}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " << "got exception: " << we.what() << "\n";} + try{bfp->GetNamedData(times_temps[i], "E:TORTGT",&TORTGT);}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " << "got exception: " << we.what() << "\n";} + // BPM Positions and Intensities - they each have 7 elements + // First is an average value of a few batches (often 2,3,4) + // used for auto-tuning, so we should disregard it + + try{HP121 = bfp->GetNamedVector(times_temps[i], "E:HP121[]");}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " << "got exception: " << we.what() << "\n";} + try{VP121 = bfp->GetNamedVector(times_temps[i], "E:VP121[]");}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " <<"got exception: " << we.what() << "\n";} + try{HPTGT = bfp->GetNamedVector(times_temps[i], "E:HPTGT[]");}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " <<"got exception: " << we.what() << "\n";} + try{VPTGT = bfp->GetNamedVector(times_temps[i], "E:VPTGT[]");}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " <<"got exception: " << we.what() << "\n";} + try{HITGT = bfp->GetNamedVector(times_temps[i], "E:HITGT[]");}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " <<"got exception: " << we.what() << "\n";} + try{VITGT = bfp->GetNamedVector(times_temps[i], "E:VITGT[]");}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " <<"got exception: " << we.what() << "\n";} + try{MTGTDS = bfp->GetNamedVector(times_temps[i], "E:MTGTDS[]");}catch (WebAPIException &we) {std::cout << "At time : " << times_temps[i] << " " <<"got exception: " << we.what() << "\n";} + + std::cout << "Finished getting IFBeam info" << std::endl; + std::cout << "BFP Time: " << times_temps[i] << " TOROID Time: " << TRTGTD_time << " TOROID COUNT: " << TRTGTD << std::endl; + unsigned long int time_closest_int = (int) TRTGTD_time; + double time_closest_ns = (TRTGTD_time - time_closest_int)*1000000000; + + sbn::NuMISpillInfo NuMIbeamInfo; + NuMIbeamInfo.TORTGT = TORTGT*1e12; //include factor of 1e12 protons in POT calculation + NuMIbeamInfo.TOR101 = TOR101*1e12; //include factor of 1e12 protons in POT calculation + NuMIbeamInfo.TRTGTD = TRTGTD*1e12; //include factor of 1e12 protons in POT calculation + NuMIbeamInfo.TR101D = TR101D*1e12; //include factor of 1e12 protons in POT calculation + NuMIbeamInfo.HRNDIR = HRNDIR; + NuMIbeamInfo.NSLINA = NSLINA; + NuMIbeamInfo.NSLINB = NSLINB; + NuMIbeamInfo.NSLINC = NSLINC; + NuMIbeamInfo.NSLIND = NSLIND; + + NuMIbeamInfo.HP121 = HP121; + NuMIbeamInfo.VP121 = VP121; + NuMIbeamInfo.HPTGT = HPTGT; + NuMIbeamInfo.VPTGT = VPTGT; + NuMIbeamInfo.HITGT = HITGT; + NuMIbeamInfo.VITGT = VITGT; + NuMIbeamInfo.MTGTDS = MTGTDS; + + NuMIbeamInfo.time = times_temps[i]; + NuMIbeamInfo.event = e.event(); + NuMIbeamInfo.spill_time_s = time_closest_int; + NuMIbeamInfo.spill_time_ns = time_closest_ns; // Save the Number of DAQ Gates in the first saved spill - if (spill_count == 1) { - NuMIbeamInfo.daq_gates = number_of_gates_since_previous_event; - } - else { - NuMIbeamInfo.daq_gates = 0; + if (spill_count == 1) { + NuMIbeamInfo.daq_gates = number_of_gates_since_previous_event; } - - fOutbeamInfos.push_back(NuMIbeamInfo); + else { + NuMIbeamInfo.daq_gates = 0; } - if(spill_count > number_of_gates_since_previous_event) - std::cout << "Event Spills : " << spill_count << ", DAQ Spills : " << number_of_gates_since_previous_event << " \t \t ::: WRONG!"<< std::endl; - else - std::cout << "Event Spills : " << spill_count << ", DAQ Spills : " << number_of_gates_since_previous_event << std::endl; - } + + fOutbeamInfos.push_back(NuMIbeamInfo); + } + if(spill_count > number_of_gates_since_previous_event) + std::cout << "Event Spills : " << spill_count << ", DAQ Spills : " << number_of_gates_since_previous_event << " \t \t ::: WRONG!"<< std::endl; + else + std::cout << "Event Spills : " << spill_count << ", DAQ Spills : " << number_of_gates_since_previous_event << std::endl; + } void sbn::NuMIRetriever::beginSubRun(art::SubRun& sr) diff --git a/sbncode/BeamSpillInfoRetriever/job/bnbextcountinfo.fcl b/sbncode/BeamSpillInfoRetriever/job/bnbextcountinfo.fcl new file mode 100644 index 000000000..e8191e2dc --- /dev/null +++ b/sbncode/BeamSpillInfoRetriever/job/bnbextcountinfo.fcl @@ -0,0 +1,10 @@ + +BEGIN_PROLOG + +bnbextcountinfo: { + + module_type: "BNBEXTRetriever" + raw_data_label: "daq" + +} +END_PROLOG diff --git a/sbncode/BeamSpillInfoRetriever/job/extcountinfo.fcl b/sbncode/BeamSpillInfoRetriever/job/extcountinfo.fcl deleted file mode 100644 index 1ca870996..000000000 --- a/sbncode/BeamSpillInfoRetriever/job/extcountinfo.fcl +++ /dev/null @@ -1,10 +0,0 @@ - -BEGIN_PROLOG - -extcountinfo: { - - module_type: "EXTRetriever" - raw_data_label: "daq" - -} -END_PROLOG diff --git a/sbncode/BeamSpillInfoRetriever/job/numiextcountinfo.fcl b/sbncode/BeamSpillInfoRetriever/job/numiextcountinfo.fcl new file mode 100644 index 000000000..f9a9583f5 --- /dev/null +++ b/sbncode/BeamSpillInfoRetriever/job/numiextcountinfo.fcl @@ -0,0 +1,10 @@ + +BEGIN_PROLOG + +numiextcountinfo: { + + module_type: "NuMIEXTRetriever" + raw_data_label: "daq" + +} +END_PROLOG diff --git a/sbncode/BeamSpillInfoRetriever/job/run_extinfo_sbn.fcl b/sbncode/BeamSpillInfoRetriever/job/run_bnbextinfo_sbn.fcl similarity index 53% rename from sbncode/BeamSpillInfoRetriever/job/run_extinfo_sbn.fcl rename to sbncode/BeamSpillInfoRetriever/job/run_bnbextinfo_sbn.fcl index e55cb29fe..fbee5871f 100644 --- a/sbncode/BeamSpillInfoRetriever/job/run_extinfo_sbn.fcl +++ b/sbncode/BeamSpillInfoRetriever/job/run_bnbextinfo_sbn.fcl @@ -1,6 +1,6 @@ -#include "extcountinfo.fcl" +#include "bnbextcountinfo.fcl" -process_name: EXTInfoGen +process_name: BNBEXTInfoGen services:{ IFBeam:{} @@ -13,10 +13,10 @@ source: { physics: { producers: { - extinfo: @local::extcountinfo + bnbextinfo: @local::bnbextcountinfo } - simulate: [extinfo ] + simulate: [bnbextinfo ] stream1: [ out1 ] } @@ -24,7 +24,7 @@ physics: { outputs: { out1: { module_type: RootOutput - fileName: "%ifb_%tc_extinfo.root" + fileName: "%ifb_%tc_bnbextinfo.root" dataTier: "raw" compressionLevel: 1 } diff --git a/sbncode/BeamSpillInfoRetriever/job/run_numiextinfo_sbn.fcl b/sbncode/BeamSpillInfoRetriever/job/run_numiextinfo_sbn.fcl new file mode 100644 index 000000000..e41261cda --- /dev/null +++ b/sbncode/BeamSpillInfoRetriever/job/run_numiextinfo_sbn.fcl @@ -0,0 +1,32 @@ +#include "numiextcountinfo.fcl" + +process_name: NuMIEXTInfoGen + +services:{ + IFBeam:{} +} + + +source: { + +} + +physics: { + producers: { + numiextinfo: @local::numiextcountinfo + } + + simulate: [numiextinfo ] + stream1: [ out1 ] + +} + +outputs: { + out1: { + module_type: RootOutput + fileName: "%ifb_%tc_numiextinfo.root" + dataTier: "raw" + compressionLevel: 1 + } +} + diff --git a/sbncode/CAFMaker/CAFMakerParams.h b/sbncode/CAFMaker/CAFMakerParams.h index ad75a102d..42f1c9d81 100644 --- a/sbncode/CAFMaker/CAFMakerParams.h +++ b/sbncode/CAFMaker/CAFMakerParams.h @@ -95,6 +95,18 @@ namespace caf "numiinfo" }; + Atom OffbeamBNBCountDataLabel { + Name("OffbeamBNBCountDataLabel"), + Comment("Label of BNB EXT module"), + "bnbextinfo" + }; + + Atom OffbeamNuMICountDataLabel { + Name("OffbeamNuMICountDataLabel"), + Comment("Label of NuMI EXT module"), + "numiextinfo" + }; + Atom G4Label { Name("G4Label"), Comment("Label of G4 module."), diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index 46a4eaea9..0f67659bb 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/sbncode/CAFMaker/CAFMaker_module.cc @@ -106,6 +106,7 @@ #include "sbnobj/Common/Reco/ScatterClosestApproach.h" #include "sbnobj/Common/Reco/StoppingChi2Fit.h" #include "sbnobj/Common/POTAccounting/BNBSpillInfo.h" +#include "sbnobj/Common/POTAccounting/EXTCountInfo.h" #include "sbnobj/Common/POTAccounting/NuMISpillInfo.h" #include "sbnobj/Common/Trigger/ExtraTriggerInfo.h" #include "sbnobj/Common/Reco/CRUMBSResult.h" @@ -184,6 +185,8 @@ class CAFMaker : public art::EDProducer { int fFileNumber; double fTotalPOT; double fSubRunPOT; + unsigned int fOffbeamBNBGates; + unsigned int fOffbeamNuMIGates; double fTotalSinglePOT; double fTotalEvents; double fBlindEvents; @@ -726,25 +729,59 @@ void CAFMaker::beginSubRun(art::SubRun& sr) { fBNBInfo.clear(); fNuMIInfo.clear(); fSubRunPOT = 0; + fOffbeamBNBGates = 0; + fOffbeamNuMIGates = 0; + + auto bnb_spill = sr.getHandle>(fParams.BNBPOTDataLabel()); + auto numi_spill = sr.getHandle>(fParams.NuMIPOTDataLabel()); + auto bnb_offbeam_spill = sr.getHandle>(fParams.OffbeamBNBCountDataLabel()); + auto numi_offbeam_spill = sr.getHandle>(fParams.OffbeamNuMICountDataLabel()); + + if(bool(bnb_spill) + bool(numi_spill) + bool(bnb_offbeam_spill) + bool(numi_offbeam_spill) > 1) { + std::cout << "Expected at most one of " << fParams.BNBPOTDataLabel() << ", " + << fParams.NuMIPOTDataLabel() << ", " << fParams.OffbeamBNBCountDataLabel() << ", and " + << fParams.OffbeamNuMICountDataLabel() << ". Found "; + if(bnb_spill) std::cout << fParams.BNBPOTDataLabel() << " "; + if(numi_spill) std::cout << fParams.NuMIPOTDataLabel() << " "; + if(bnb_offbeam_spill) std::cout << fParams.OffbeamBNBCountDataLabel() << " "; + if(numi_offbeam_spill) std::cout << fParams.OffbeamNuMICountDataLabel(); + std::cout << std::endl; + abort(); + } - if(auto bnb_spill = sr.getHandle>(fParams.BNBPOTDataLabel())){ + if(bnb_spill){ FillExposure(*bnb_spill, fBNBInfo, fSubRunPOT); fTotalPOT += fSubRunPOT; } - else if (auto numi_spill = sr.getHandle>(fParams.NuMIPOTDataLabel())) { + else if (numi_spill) { FillExposureNuMI(*numi_spill, fNuMIInfo, fSubRunPOT); fTotalPOT += fSubRunPOT; } + else if (bnb_offbeam_spill){ + for(const auto& spill: *bnb_offbeam_spill) { + fOffbeamBNBGates += spill.gates_since_last_trigger; + } + } + else if (numi_offbeam_spill){ + for(const auto& spill: *numi_offbeam_spill) { + fOffbeamNuMIGates += spill.gates_since_last_trigger; + } + } else if(auto pot_handle = sr.getHandle(fParams.GenLabel())){ fSubRunPOT = pot_handle->totgoodpot; fTotalPOT += fSubRunPOT; } else{ - if(!fParams.BNBPOTDataLabel().empty() || !fParams.GenLabel().empty() || !fParams.NuMIPOTDataLabel().empty()){ + if(!fParams.BNBPOTDataLabel().empty() || !fParams.GenLabel().empty() || !fParams.NuMIPOTDataLabel().empty() || + !fParams.OffbeamBNBCountDataLabel().empty() || !fParams.OffbeamNuMICountDataLabel().empty()){ std::cout << "Found neither BNB data POT info under '" << fParams.BNBPOTDataLabel() - << "' not NuMIdata POT info under '" + << "' nor NuMIdata POT info under '" << fParams.NuMIPOTDataLabel() + << "' nor BNB EXT Count info under '" + << fParams.OffbeamBNBCountDataLabel() + << "' nor NuMI EXT Count info under '" + << fParams.OffbeamNuMICountDataLabel() << "' nor MC POT info under '" << fParams.GenLabel() << "'" << std::endl; @@ -1735,9 +1772,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); } @@ -1888,6 +1926,8 @@ void CAFMaker::produce(art::Event& evt) noexcept { rec.hdr.bnbinfo = fBNBInfo; rec.hdr.nnumiinfo = fNuMIInfo.size(); rec.hdr.numiinfo = fNuMIInfo; + rec.hdr.noffbeambnb = fOffbeamBNBGates; + rec.hdr.noffbeamnumi = fOffbeamNuMIGates; rec.hdr.pot = fSubRunPOT; } 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) { 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; } 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 diff --git a/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx b/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx index 8ae583973..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" @@ -15,7 +18,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 +32,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 +57,21 @@ 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 { + // 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); @@ -69,6 +83,8 @@ void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj auto nextValid=traj.FirstValidPoint(); breakpoints.push_back(nextValid); auto pos0 = traj.LocationAtPoint(nextValid); + bool pos0good = isInVolume(fiducialVolumes, pos0) && !isInVolume(excludeVolumes, pos0); + breakpointsgood.push_back(pos0good); nextValid = traj.NextValidPoint(nextValid+1); int npoints = 0; while (nextValid!=recob::TrackTrajectory::InvalidIndex) { @@ -77,7 +93,9 @@ void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj pos0=pos1; npoints++; if (thislen>=thisSegLen) { + 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); else segradlengths.push_back(-999.); cumseglens.push_back(cumseglens.back()+thislen); @@ -88,7 +106,10 @@ void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj } //then add last segment if (thislen>0.) { + auto endpointpos = traj.LocationAtPoint(nextValid); + 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); cumseglens.push_back(cumseglens.back()+thislen); } @@ -202,7 +223,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 +322,73 @@ double TrajectoryMCSFitter::GetE(const double initial_E, const double length_tra } return current_E; } + +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)); + } + 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 07b63c443..6d3b48032 100644 --- a/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h +++ b/sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h @@ -3,11 +3,15 @@ #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" #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 { /** @@ -87,10 +91,18 @@ namespace trkf::sbn { Comment("Angular resolution parameter used in modified Highland formula. Unit is mrad."), 3.0 }; + 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){ + 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; @@ -101,9 +113,11 @@ namespace trkf::sbn { pMax_ = pMax; pStep_ = pStep; angResol_ = angResol; + 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()) {} + : 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); } @@ -117,7 +131,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 +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; + std::vector setFiducialVolumes() const; + std::vector setExcludeVolumes() const; + bool isInVolume(const std::vector &volumes, const geo::Point_t &point) const; // private: int pIdHyp_; @@ -158,6 +175,8 @@ namespace trkf::sbn { double pMax_; double pStep_; double angResol_; + std::vector fiducialVolumeInsets_; + std::vector excludeVolumes_; }; } diff --git a/sbncode/LArRecoProducer/mcsproducer.fcl b/sbncode/LArRecoProducer/mcsproducer.fcl index e41a6d00b..5db98243c 100644 --- a/sbncode/LArRecoProducer/mcsproducer.fcl +++ b/sbncode/LArRecoProducer/mcsproducer.fcl @@ -1,7 +1,10 @@ BEGIN_PROLOG mcs_sbn: { module_type: MCSFitAllPID - MCS: {} + MCS: { + fiducialVolumeInsets: [] + excludeVolumes: [] + } TrackLabel: pandoraTrack } diff --git a/sbncode/LArRecoProducer/mcsproducer_icarus.fcl b/sbncode/LArRecoProducer/mcsproducer_icarus.fcl new file mode 100644 index 000000000..4926fdfc7 --- /dev/null +++ b/sbncode/LArRecoProducer/mcsproducer_icarus.fcl @@ -0,0 +1,13 @@ +BEGIN_PROLOG +mcs_sbn: { + module_type: MCSFitAllPID + MCS: { + eLossMode: 2 + angResol: 10.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 +} + +END_PROLOG diff --git a/sbncode/SBNEventWeight/App/CMakeLists.txt b/sbncode/SBNEventWeight/App/CMakeLists.txt index 2989659b7..335a19557 100644 --- a/sbncode/SBNEventWeight/App/CMakeLists.txt +++ b/sbncode/SBNEventWeight/App/CMakeLists.txt @@ -1,8 +1,11 @@ cet_build_plugin(SBNEventWeight art::module LIBRARIES - sbnobj::Common_SBNEventWeight sbncode_SBNEventWeight_Base + sbnobj::Common_SBNEventWeight + sbncode_SBNEventWeight_Base sbncode_SBNEventWeight_Calculators_CrossSection - sbncode_SBNEventWeight_Calculators_BNBFlux nusimdata::SimulationBase + sbncode_SBNEventWeight_Calculators_BNBFlux + sbncode_SBNEventWeight_Calculators_Geant4 + nusimdata::SimulationBase ROOT::Geom) cet_build_plugin(SystToolsEventWeight art::module diff --git a/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc b/sbncode/SBNEventWeight/App/SystToolsEventWeight_module.cc index 4793171b9..6e836054c 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,22 +139,32 @@ 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; + 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" << "[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 " @@ -159,8 +172,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,25 +204,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 (fDebugMode) { - std::cout << " sph.prettyName = " << sph.prettyName << std::endl; + 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; + } + } + } - std::string rwmode = "multisigma"; - if(sph.isRandomlyThrown) rwmode = "multisim"; + for( auto &sph : smd ){ - if (fDebugMode) { - std::cout << " rwmode = " << rwmode << std::endl; + // responsless + if(sph.isResponselessParam){ + if (fDebugMode) { + std::cout << "Responsless dial found: " << sph.prettyName << ", thus skipping" << std::endl; + } + continue; } - 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()); + 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{ + 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 { paramVars_dep.begin(), paramVars_dep.end() }; + fParameterSet.AddParameter(sph_dep.prettyName, std::move(widths_dep)); + + } + + + + } + else{ + + rwmode = sph.isRandomlyThrown ? "multisim" : "multisigma"; + + std::vector widths { paramVars.begin(), paramVars.end() }; + + fParameterSet.AddParameter(sph.prettyName, std::move(widths)); + + } + + if(rwmode==""){ + throw cet::exception{ "SystToolsEventWeight" } + << "rwmode not set for " << sph.prettyName << "\n"; + } + + if (fDebugMode) { + std::cout << " sph.prettyName = " << sph.prettyName << ", rwmode = " << rwmode << " is added to the header" << std::endl; + } + + fParameterSet.Configure(sp->GetFullyQualifiedName()+"_"+sph.prettyName, rwmode, paramVars.size()); fParameterSet.FillKnobValues(); p->push_back( std::move(fParameterSet) ); 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/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..568ecc7f6 --- /dev/null +++ b/sbncode/SBNEventWeight/Calculators/Geant4/CMakeLists.txt @@ -0,0 +1,24 @@ +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 + LIBRARIES + sbncode_SBNEventWeight_Base + nusimdata::SimulationBase + nurandom::RandomUtils_NuRandomService_service + art::Framework_Principal + art::Framework_Services_Registry + art_root_io::TFileService_service + larcore::Geometry_Geometry_service + cetlib_except::cetlib_except + ReweightBaseLib + PropBaseLib + ROOT::Tree +) + +install_headers() +install_fhicl() +install_source() + diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx new file mode 100644 index 000000000..a0cead5d2 --- /dev/null +++ b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx @@ -0,0 +1,492 @@ +/** + * \class evwgh::Geant4WeightCalc + * \brief Updated hadron reinteraction event reweighting, using geant4reweight + * \author K. Duffy , 2019/10 + * + * Reweight events based on hadron reinteraction probabilities. + */ + +#include +#include +#include "TDirectory.h" +#include "TFile.h" +#include "TH1D.h" +#include "TTree.h" +#include "Geant4/G4LossTableManager.hh" +#include "Geant4/G4ParticleTable.hh" +#include "Geant4/G4ParticleDefinition.hh" +#include "Geant4/G4Material.hh" +#include "Geant4/G4MaterialCutsCouple.hh" +#include "art_root_io/TFileService.h" +#include "art_root_io/TFileDirectory.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Persistency/Provenance/ModuleContext.h" +#include "art/Framework/Principal/Event.h" +#include "canvas/Persistency/Common/FindManyP.h" +#include "canvas/Persistency/Common/Ptr.h" +#include "nusimdata/SimulationBase/MCParticle.h" +#include "nusimdata/SimulationBase/MCTruth.h" +#include "CLHEP/Random/RandGaussQ.h" +#include "CLHEP/Units/PhysicalConstants.h" +#include "CLHEP/Units/SystemOfUnits.h" +#include "sbncode/SBNEventWeight/Base/WeightCalcCreator.h" +#include "sbncode/SBNEventWeight/Base/WeightCalc.h" +#include "larcore/Geometry/Geometry.h" +#include "geant4reweight/src/ReweightBase/G4ReweighterFactory.hh" +#include "geant4reweight/src/ReweightBase/G4Reweighter.hh" +#include "geant4reweight/src/ReweightBase/G4ReweightTraj.hh" +#include "geant4reweight/src/ReweightBase/G4ReweightStep.hh" +#include "geant4reweight/src/PropBase/G4ReweightParameterMaker.hh" + +// local include +#include "BetheBlochForG4ReweightValid.h" + +namespace sbn { + +namespace evwgh { + +class Geant4WeightCalc : public WeightCalc { +public: + Geant4WeightCalc() {} + + void Configure(fhicl::ParameterSet const& p, CLHEP::HepRandomEngine& engine); + + std::vector GetWeight(art::Event& e, size_t inu); + +private: + std::string fMCParticleProducer; //!< Label for the MCParticle producer + std::string fMCTruthProducer; //!< Label for the MCTruth producer + CLHEP::RandGaussQ* fGaussRandom; //!< Random number generator + // std::map fParticles; //!< Particles to reweight + unsigned fNsims; //!< Number of multisims + int fPdg; //!< PDG value for particles that a given weight calculator should apply to. Note that for now this module can only handle weights for one particle species at a time. + // float fXSUncertainty; //!< Flat cross section uncertainty + G4ReweighterFactory RWFactory; //!< Base class to handle all Geant4Reweighters (right now "all" means pi+, pi-, p) + G4Reweighter *theReweighter; //!< Geant4Reweighter -- this is what provides the weights + G4ReweightParameterMaker *ParMaker; + std::vector> UniverseVals; //!< Vector of maps relating parameter name to value (defines parameter values that will be evaluated in universes). Each map should have one entry per parameter we are considering + + art::ServiceHandle < geo::Geometry > fGeometryService; + + bool fMakeOutputTrees; ///!< Fcl parameter to decide whether to save output tree (useful for validations but not necessary when running in production) + TTree *fOutTree_MCTruth; //!< Output tree for quick validations: on entry per MCTruth object + TTree *fOutTree_Particle; //!< Output tree for quick validations: one entry per neutrino-induced pi+, pi-, or proton + int event_num; //!< Variables for both output trees + int run_num; //!< Variables for both output trees + int subrun_num; //!< Variables for both output trees + double p_track_length; //!< Variables for by-particle output tree + int p_PDG; //!< Variables for by-particle output tree + std::string p_final_proc; //!< Variables for by-particle output tree + double p_init_momentum; //!< Variables for by-particle output tree + double p_final_momentum; //!< Variables for by-particle output tree + std::vector< double > p_energies_el; //!< Variables for by-particle output tree + std::vector< int > p_sliceInts_el; //!< Variables for by-particle output tree + std::vector< double > p_energies_inel; //!< Variables for by-particle output tree + std::vector< int > p_sliceInts_inel; //!< Variables for by-particle output tree + int p_nElasticScatters; //!< Variables for by-particle output tree + std::vector p_inel_weight; //!< Variables for by-particle output tree + std::vector p_elastic_weight; //!< Variables for by-particle output + + bool fDebug; //!< print debug statements + + DECLARE_WEIGHTCALC(Geant4WeightCalc) +}; + + +void Geant4WeightCalc::Configure(fhicl::ParameterSet const& p, + CLHEP::HepRandomEngine& engine) +{ + std::cout << "Using Geant4WeightCalc for reinteraction weights" << std::endl; + + // Get configuration for this function + fhicl::ParameterSet const& pset = p.get(GetName()); + // std::cout << pset.GetName() << std::endl; + fMCParticleProducer = pset.get("MCParticleProducer", "largeant"); + fMCTruthProducer = pset.get("MCTruthProducer", "generator"); + fMakeOutputTrees = pset.get< bool >( "makeoutputtree", false ); + std::string mode = pset.get("mode"); + std::string FracsFileName = pset.get< std::string >( "fracsfile" ); + std::string XSecFileName = pset.get< std::string >( "xsecfile" ); + std::vector< fhicl::ParameterSet > FitParSets = pset.get< std::vector< fhicl::ParameterSet > >("parameters"); + fNsims = pset.get ("number_of_multisims", 0); + fPdg = pset.get ("pdg_to_reweight"); + fDebug = pset.get ("debug",false); + + // Prepare random generator + fGaussRandom = new CLHEP::RandGaussQ(engine); + + // Get input files + TFile FracsFile( FracsFileName.c_str(), "OPEN" ); + TFile XSecFile( XSecFileName.c_str(), "OPEN" ); + + // Configure G4Reweighter + bool totalonly = false; + if (fPdg==2212) totalonly = true; + ParMaker = new G4ReweightParameterMaker( FitParSets, totalonly ); + theReweighter = RWFactory.BuildReweighter(fPdg, &XSecFile, &FracsFile, ParMaker->GetFSHists(), ParMaker->GetElasticHist() ); + + // Make output trees to save things for quick and easy validation + art::ServiceHandle tfs; + + if (fMakeOutputTrees){ + fOutTree_Particle = tfs->make(Form("%s_%i","ByParticleValidTree",fPdg),""); + fOutTree_Particle->Branch("event_num",&event_num); + fOutTree_Particle->Branch("run_num",&run_num); + fOutTree_Particle->Branch("subrun_num",&subrun_num); + fOutTree_Particle->Branch("UniverseVals", &UniverseVals); + fOutTree_Particle->Branch("pdg_to_reweight",&fPdg); + fOutTree_Particle->Branch("inelastic_weight",&p_inel_weight); + fOutTree_Particle->Branch("elastic_weight",&p_elastic_weight); + fOutTree_Particle->Branch("track_length",&p_track_length); + fOutTree_Particle->Branch("PDG",&p_PDG); + fOutTree_Particle->Branch("final_proc",&p_final_proc); + fOutTree_Particle->Branch("init_momentum",&p_init_momentum); + fOutTree_Particle->Branch("final_momentum",&p_final_momentum); + fOutTree_Particle->Branch("energies_el",&p_energies_el); + fOutTree_Particle->Branch("sliceInts_el",&p_sliceInts_el); + fOutTree_Particle->Branch("energies_inel",&p_energies_inel); + fOutTree_Particle->Branch("sliceInts_inel",&p_sliceInts_inel); + fOutTree_Particle->Branch("nElasticScatters",&p_nElasticScatters); + + fOutTree_MCTruth = tfs->make(Form("%s_%i","ByMCTruthValidTree",fPdg),""); + fOutTree_MCTruth->Branch("event_num",&event_num); + fOutTree_MCTruth->Branch("run_num",&run_num); + fOutTree_MCTruth->Branch("subrun_num",&subrun_num); + fOutTree_MCTruth->Branch("UniverseVals", &UniverseVals); + fOutTree_MCTruth->Branch("pdg_to_reweight",&fPdg); + } + + + // Read input parameter sets and set up universes + size_t n_parsets = FitParSets.size(); + std::vector FitParNames; + std::vector FitParNominals; + std::vector FitParSigmas; + std::map theNominals; + + fParameterSet.Configure(GetFullName(), mode, fNsims); + + for (size_t i_parset=0; i_parset("Name"); + double theNominal = theSet.get("Nominal",1.); + double theSigma = theSet.get("Sigma",0.); + + FitParNames.push_back(theName); + FitParNominals.push_back(theNominal); + FitParSigmas.push_back(theSigma); + + theNominals[theName] = theNominal; + + fParameterSet.AddParameter(theName, theSigma); + } + + if (mode=="pm1sigma"){ + // pm1sigma mode: 0 = +1sigma, 1 = -1sigma of a single parameter. All other parameters at nominal + for (size_t i_parset=0; i_parset tmp_vals_p1sigma(theNominals); + std::map tmp_vals_m1sigma(theNominals); + // Now reset the +1sigma and -1sigma values for this parameter set only + tmp_vals_p1sigma[FitParNames.at(i_parset)] = FitParNominals.at(i_parset)+FitParSigmas.at(i_parset); + tmp_vals_m1sigma[FitParNames.at(i_parset)] = FitParNominals.at(i_parset)-FitParSigmas.at(i_parset); + + if (fDebug){ + std::cout << "Universe " << i_parset*2 << ": " << FitParNames.at(i_parset) << " = " << FitParNominals.at(i_parset)+FitParSigmas.at(i_parset) << std::endl; + std::cout << "Universe " << i_parset*2+1 << ": " << FitParNames.at(i_parset) << " = " << FitParNominals.at(i_parset)-FitParSigmas.at(i_parset) << std::endl; + } + + // Finally, add these universes into the vector + UniverseVals.push_back(tmp_vals_p1sigma); + UniverseVals.push_back(tmp_vals_m1sigma); + } // end loop over parsets (i) + } // pm1sigma + else if (mode=="multisim"){ + // multisim mode: parameter values sample within the given uncertainty for all parameters simultaneously + // Loop over universes j + for (unsigned j=0; j tmp_vals; + for (size_t i_parset=0; i_parsetfire(0.0,1.0); + tmp_vals[FitParNames.at(i_parset)] = FitParNominals.at(i_parset)+(FitParSigmas.at(i_parset)*r); + } // loop over parameters (i_parset) + // Now save this universe + UniverseVals.push_back(tmp_vals); + } // loop over Nsims (j) + } // multisim + else{ + // Anything else mode: Set parameters to user-defined nominal value + UniverseVals.push_back(theNominals); + } // any other mode + + fNsims = UniverseVals.size(); + if (fDebug) std::cout << "Running mode: " << mode <<". Nsims = " << fNsims << std::endl; +} + + +std::vector Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) { + + // Get event/run/subrun numbers for output + run_num = e.run(); + subrun_num = e.subRun(); + event_num = e.id().event(); + + // Get MCParticles for each MCTruth in this event + art::Handle > truthHandle; + e.getByLabel(fMCTruthProducer, truthHandle); + const art::FindManyP truthParticles(truthHandle, e, fMCParticleProducer); + assert(truthParticles.isValid()); + + // Initialize the vector of event weights + std::vector weight; + + // Initialize weight vector for this MCTruth + weight.clear(); + weight.resize(fNsims, 1.0); + + // Loop over MCParticles in the event + auto const& mcparticles = truthParticles.at(itruth); + + for (size_t i=0; i trajpoint_X; + std::vector trajpoint_Y; + std::vector trajpoint_Z; + std::vector trajpoint_PX; + std::vector trajpoint_PY; + std::vector trajpoint_PZ; + std::vector elastic_indices; + + //Get the list of processes from the true trajectory + const std::vector< std::pair< size_t, unsigned char > > & processes = p.Trajectory().TrajectoryProcesses(); + std::map< size_t, std::string > process_map; + for( auto it = processes.begin(); it != processes.end(); ++it ){ + process_map[ it->first ] = p.Trajectory().KeyToProcess( it->second ); + } + + for( size_t i = 0; i < p.NumberTrajectoryPoints(); ++i ){ + double X = p.Position(i).X(); + double Y = p.Position(i).Y(); + double Z = p.Position(i).Z(); + geo::Point_t testpoint1 { X, Y, Z }; + const TGeoMaterial* testmaterial1 = fGeometryService->Material( testpoint1 ); + //For now, just going to reweight the points within the LAr of the TPC + // TODO check if this is right + if ( !strcmp( testmaterial1->GetName(), "LAr" ) ){ + trajpoint_X.push_back( X ); + trajpoint_Y.push_back( Y ); + trajpoint_Z.push_back( Z ); + + trajpoint_PX.push_back( p.Px(i) ); + trajpoint_PY.push_back( p.Py(i) ); + trajpoint_PZ.push_back( p.Pz(i) ); + + auto itProc = process_map.find(i); + if( itProc != process_map.end() && itProc->second == "hadElastic" ){ + //Push back the index relative to the start of the reweightable steps + elastic_indices.push_back( trajpoint_X.size() - 1 ); + // if (fDebug) std::cout << "Elastic index: " << trajpoint_X.size() - 1 << std::endl; + } + + } + + } // end loop over trajectory points + + // Now find daughters of the MCP + std::vector daughter_PDGs; + std::vector daughter_IDs; + for( int i_mcp = 0; i_mcp < p.NumberDaughters(); i_mcp++ ){ + int daughterID = p.Daughter(i_mcp); + for (auto test_mcp : mcparticles){ + if (test_mcp->TrackId() == daughterID){ + int pid = test_mcp->PdgCode(); + daughter_PDGs.push_back(pid); + daughter_IDs.push_back( test_mcp->TrackId() ); + break; + } + } + } // end loop over daughters + + // --- Now that we have all the information about the track we need, here comes the reweighting part! --- // + + //Make a G4ReweightTraj -- This is the reweightable object + G4ReweightTraj theTraj(mcpID, p_PDG, 0, event_num, std::make_pair(0,0)); + + //Create its set of G4ReweightSteps and add them to the Traj (note: this needs to be done once per MCParticle but will be valid for all weight calculations) + std::vector< G4ReweightStep* > allSteps; + + size_t nSteps = trajpoint_PX.size(); + + if( nSteps < 2 ) continue; + + p_nElasticScatters = elastic_indices.size(); + for( size_t istep = 1; istep < nSteps; ++istep ){ + + // if( istep == trajpoint_PX.size() - 1 && std::find( elastic_indices.begin(), elastic_indices.end(), j ) != elastic_indices.end() ) + // std::cout << "Warning: last step an elastic process" << std::endl; + + std::string proc = "default"; + if( istep == trajpoint_PX.size() - 1 ) + proc = EndProcess; + else if( std::find( elastic_indices.begin(), elastic_indices.end(), istep ) != elastic_indices.end() ){ + proc = "hadElastic"; + } + + + double deltaX = ( trajpoint_X.at(istep) - trajpoint_X.at(istep-1) ); + double deltaY = ( trajpoint_Y.at(istep) - trajpoint_Y.at(istep-1) ); + double deltaZ = ( trajpoint_Z.at(istep) - trajpoint_Z.at(istep-1) ); + + double len = sqrt( + std::pow( deltaX, 2 ) + + std::pow( deltaY, 2 ) + + std::pow( deltaZ, 2 ) + ); + + double preStepP[3] = { + trajpoint_PX.at(istep-1)*1.e3, + trajpoint_PY.at(istep-1)*1.e3, + trajpoint_PZ.at(istep-1)*1.e3 + }; + + double postStepP[3] = { + trajpoint_PX.at(istep)*1.e3, + trajpoint_PY.at(istep)*1.e3, + trajpoint_PZ.at(istep)*1.e3 + }; + + if( istep == 1 ){ + theTraj.SetEnergy( sqrt( preStepP[0]*preStepP[0] + preStepP[1]*preStepP[1] + preStepP[2]*preStepP[2] + mass*mass ) ); + } + + G4ReweightStep * theStep = new G4ReweightStep( mcpID, p_PDG, 0, event_num, preStepP, postStepP, len, proc ); + theStep->SetDeltaX( deltaX ); + theStep->SetDeltaY( deltaY ); + theStep->SetDeltaZ( deltaZ ); + + theTraj.AddStep( theStep ); + + for( size_t k = 0; k < daughter_PDGs.size(); ++k ){ + theTraj.AddChild( new G4ReweightTraj(daughter_IDs[k], daughter_PDGs[k], mcpID, event_num, std::make_pair(0,0) ) ); + } + } // end loop over nSteps (istep) + p_track_length = theTraj.GetTotalLength(); + + p_init_momentum = sqrt( theTraj.GetEnergy()*theTraj.GetEnergy() - mass*mass ); + p_final_momentum = sqrt( + std::pow( theTraj.GetStep( theTraj.GetNSteps() - 1 )->GetPreStepPx(), 2 ) + + std::pow( theTraj.GetStep( theTraj.GetNSteps() - 1 )->GetPreStepPy(), 2 ) + + std::pow( theTraj.GetStep( theTraj.GetNSteps() - 1 )->GetPreStepPz(), 2 ) + ); + + std::vector< std::pair< double, int > > thin_slice_inelastic = ThinSliceBetheBloch( &theTraj, .5, mass , false); + std::vector< std::pair< double, int > > thin_slice_elastic = ThinSliceBetheBloch( &theTraj, .5, mass , true); + + p_energies_inel.clear(); + p_sliceInts_inel.clear(); + for( size_t islice = 0; islice < thin_slice_inelastic.size(); ++islice ){ + p_energies_inel.push_back( thin_slice_inelastic[islice].first ); + p_sliceInts_inel.push_back( thin_slice_inelastic[islice].second ); + } + + p_energies_el.clear(); + p_sliceInts_el.clear(); + for( size_t islice = 0; islice < thin_slice_elastic.size(); ++islice ){ + p_energies_el.push_back( thin_slice_elastic[islice].first ); + p_sliceInts_el.push_back( thin_slice_elastic[islice].second ); + } + + // Loop through universes (j) + for (size_t j=0; jSetNewVals(UniverseVals.at(j)); + theReweighter->SetNewHists(ParMaker->GetFSHists()); + theReweighter->SetNewElasticHists(ParMaker->GetElasticHist()); + //Get the weight from the G4ReweightTraj + w = theReweighter->GetWeight( &theTraj ); + // Total weight is the product of track weights in the event + weight[j] *= std::max((float)0.0, w); + + // Do the same for elastic weight (should be 1 unless set to non-nominal ) + el_w = theReweighter->GetElasticWeight( &theTraj ); + weight[j] *= std::max((float)0.0,el_w); + + // just for the output tree + p_inel_weight[j] = w; + p_elastic_weight[j] = el_w; + + if (fDebug){ + std::cout << " Universe " << j << ": "; + // std::cout << UniverseVals.at(j) << std::endl; + std::cout << " w = " << w << ", el_w = " << el_w << std::endl; + } + + + } // loop through universes (j) + + if (fDebug){ + std::cout << "PDG = " << p_PDG << std::endl; + std::cout << "inel weights by particle: "; + for (unsigned int j=0; jFill(); + } // loop over mcparticles (i) + if (fMakeOutputTrees) fOutTree_MCTruth->Fill(); + +return weight; + +} + +REGISTER_WEIGHTCALC(Geant4WeightCalc) + +} // namespace evwgh + +} // namespace sbn diff --git a/sbncode/SBNEventWeight/jobs/CMakeLists.txt b/sbncode/SBNEventWeight/jobs/CMakeLists.txt index 441c7dbde..7897943c9 100644 --- a/sbncode/SBNEventWeight/jobs/CMakeLists.txt +++ b/sbncode/SBNEventWeight/jobs/CMakeLists.txt @@ -4,4 +4,5 @@ install_source() add_subdirectory(genie) add_subdirectory(flux) +add_subdirectory(geant4) diff --git a/sbncode/SBNEventWeight/jobs/geant4/CMakeLists.txt b/sbncode/SBNEventWeight/jobs/geant4/CMakeLists.txt new file mode 100644 index 000000000..589bc695a --- /dev/null +++ b/sbncode/SBNEventWeight/jobs/geant4/CMakeLists.txt @@ -0,0 +1,6 @@ +install_fhicl() + +FILE(GLOB fcl_files *.fcl) + +install_source( EXTRAS ${fcl_files} ) + diff --git a/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl b/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl new file mode 100644 index 000000000..e45065880 --- /dev/null +++ b/sbncode/SBNEventWeight/jobs/geant4/eventweight_geant4_sbn.fcl @@ -0,0 +1,103 @@ +########################################################## +## Hadron reinteraction uncertainties +## +## References: +## +## * K. Duffy, Pion Secondary Interaction Systematics (from GEANTReweight), DocDB 25379 +## * J. Calcutt, GEANTReweight documentation DocDB 25084, repository https://cdcvs.fnal.gov/redmine/projects/geant4reweight/repository +## +## From MicroBooNE, ubsim UBOONE_SUITE_v08_00_00_69 +## (uBooNE author: Kirsty Duffy (kduffy@fnal.gov)) +## +########################################################## + +#include "piplus_reweight_parameters.fcl" +#include "piminus_reweight_parameters.fcl" +#include "proton_reweight_parameters.fcl" + +BEGIN_PROLOG + +sbn_eventweight_geant4: { + module_type: "SBNEventWeight" + AllowMissingTruth: true + min_weight: 0.0 + max_weight: 100 + + weight_functions_reint: [ reinteractions_piplus, reinteractions_piminus, reinteractions_proton ] + + reinteractions_piplus: { + type: Geant4 + random_seed: 58 + parameters: @local::PiPlusParameters + mode: multisim + number_of_multisims: 1000 + fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_piplus.root" + xsecfile: "$SBNDATA_DIR/systematics/reint/g4_cross_section_piplus.root" + makeoutputtree: false + pdg_to_reweight: 211 + debug: false + } + + reinteractions_piminus: { + type: Geant4 + random_seed: 59 + parameters: @local::PiMinusParameters + mode: multisim + number_of_multisims: 1000 + fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_piminus.root" + xsecfile: "$SBNDATA_DIR/systematics/reint/g4_cross_section_piminus.root" + makeoutputtree: false + pdg_to_reweight: -211 + debug: false + } + + reinteractions_proton: { + type: Geant4 + random_seed: 60 + parameters: @local::ProtonParameters + mode: multisim + number_of_multisims: 1000 + fracsfile: "$SBNDATA_DIR/systematics/reint/g4_fracs_proton.root" + xsecfile: "$SBNDATA_DIR/systematics/reint/g4_cross_section_proton.root" + makeoutputtree: false + pdg_to_reweight: 2212 + debug: false + } +} + +### +# Reweighting parameters should be defined as +# +# TheParameters: [ +# { +# Cut: "reac" +# Name: "fReacLow" +# Range: [10., 200.] +# Nominal: 1.0 +# Sigma: 0.3 +# }, +# {...} +# ] +# +# - Range defines the energy range over which that parameter has an effect (MeV) +# - Nominal is not used in "multisim" mode. In "pm1sigma" mode it defines the +# nominal (around which you calculate +/- 1 sigma variations). In any other +# mode, the parameter is set to the value under Nominal (so it's not really a +# "nominal" value in that case, just the set value) +# - Sigma defines the 1-sigma range for multisims and unisims. Currently the +# code can only accept a single value for sigma, giving symmetric +# uncertainty bounds: nominal-sigma and nominal+sigma +# - Cut must be one of the following: +# "reac" <- total inelastic scattering cross section +# "abs" <- absorption cross section (pi+ and pi- only at the moment -- Oct 2019) +# "cex" <- charge exchange cross section (pi+ and pi- only at the moment -- Oct 2019) +# "dcex" <- double charge exchange cross section (pi+ and pi- only at the moment -- Oct 2019) +# "prod" <- pion production cross section (pi+ and pi- only at the moment -- Oct 2019) +# "inel" <- quasi-elastic inelastic scattering cross section (pi+ and pi- only at the moment -- Oct 2019) +# "elast" <- total elastic scattering cross section +# - Name can be anything you like (but should uniquely identify that parameter/ +# variation) +### + +END_PROLOG + diff --git a/sbncode/SBNEventWeight/jobs/geant4/piminus_reweight_parameters.fcl b/sbncode/SBNEventWeight/jobs/geant4/piminus_reweight_parameters.fcl new file mode 100644 index 000000000..b64756415 --- /dev/null +++ b/sbncode/SBNEventWeight/jobs/geant4/piminus_reweight_parameters.fcl @@ -0,0 +1,66 @@ +BEGIN_PROLOG + +# From MicroBooNE, ubsim UBOONE_SUITE_v08_00_00_69 + +PiMinusParameters: [ + { + Name: "fPiMinusReacLow" + Cut: "reac" + Range: [10., 200.] + Nominal: 1 + Sigma: 0.2 + } , + { + Name: "fPiMinusReacHigh" + Cut: "reac" + Range: [700., 2005.] + Nominal: 1.0 + Sigma: 0.2 + } , + + { + Name: "fPiMinusAbs" + Cut: "abs" + # Range: [200.0, 700.00] + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + }, + + { + Name: "fPiMinusCex" + Cut: "cex" + # Range: [200.0, 700.00] + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + }, + + { + Name: "fPiMinusDCex" + Cut: "dcex" + # Range: [200.0, 700.00] + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + }, + + { + Name: "fPiMinusPiProd" + Cut: "prod" + # Range: [200.0, 700.00] + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + }, + + { + Name: "fPiMinusElast" + Cut: "elast" + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + } +] + +END_PROLOG diff --git a/sbncode/SBNEventWeight/jobs/geant4/piplus_reweight_parameters.fcl b/sbncode/SBNEventWeight/jobs/geant4/piplus_reweight_parameters.fcl new file mode 100644 index 000000000..c9cfc3e23 --- /dev/null +++ b/sbncode/SBNEventWeight/jobs/geant4/piplus_reweight_parameters.fcl @@ -0,0 +1,66 @@ +BEGIN_PROLOG + +# From MicroBooNE, ubsim UBOONE_SUITE_v08_00_00_69 + +PiPlusParameters: [ + { + Name: "fPiPlusReacLow" + Cut: "reac" + Range: [10., 200.] + Nominal: 1 + Sigma: 0.2 + } , + { + Name: "fPiPlusReacHigh" + Cut: "reac" + Range: [700., 2005.] + Nominal: 1.0 + Sigma: 0.2 + } , + + { + Name: "fPiPlusAbs" + Cut: "abs" + # Range: [200.0, 700.00] + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + }, + + { + Name: "fPiPlusCex" + Cut: "cex" + # Range: [200.0, 700.00] + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + }, + + { + Name: "fPiPlusDCex" + Cut: "dcex" + # Range: [200.0, 700.00] + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + }, + + { + Name: "fPiPlusPiProd" + Cut: "prod" + # Range: [200.0, 700.00] + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + }, + + { + Name: "fPiPlusElast" + Cut: "elast" + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + } +] + +END_PROLOG diff --git a/sbncode/SBNEventWeight/jobs/geant4/proton_reweight_parameters.fcl b/sbncode/SBNEventWeight/jobs/geant4/proton_reweight_parameters.fcl new file mode 100644 index 000000000..13113db75 --- /dev/null +++ b/sbncode/SBNEventWeight/jobs/geant4/proton_reweight_parameters.fcl @@ -0,0 +1,23 @@ +BEGIN_PROLOG + +# From MicroBooNE, ubsim UBOONE_SUITE_v08_00_00_69 + +ProtonParameters: [ + { + Name: "fProtonReac" + Cut: "reac" + Range: [10., 2005.] + Nominal: 1 + Sigma: 0.2 + }, + + { + Name: "fProtonElast" + Cut: "elast" + Range: [10.0, 2005.00] + Nominal: 1.0 + Sigma: 0.2 + } +] + +END_PROLOG diff --git a/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 diff --git a/ups/product_deps b/ups/product_deps index 8d2b04768..e9aae535b 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 @@ -257,10 +257,11 @@ larcv2 v2_1_0 - larsoft v09_72_00 - sbnanaobj v09_20_06 - sbndaq_artdaq_core v1_06_00of0 - -sbndata v01_04 - +sbndata v01_05 - sbnobj v09_16_00 - -systematicstools v01_02_00 - -nusystematics v01_02_10 - +systematicstools v01_02_01 - +nusystematics v01_02_11 - +geant4reweight v01_00_03e - 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 ####################################