diff --git a/CMakeLists.txt b/CMakeLists.txt index d797dc2d1..3bafc97d9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,7 +16,7 @@ cmake_minimum_required(VERSION 3.20 FATAL_ERROR) find_package(cetmodules 3.20.00 REQUIRED) -project(sbncode VERSION 10.14.02.03 LANGUAGES CXX) +project(sbncode VERSION 10.14.02.04 LANGUAGES CXX) message(STATUS "\n\n ========================== ${PROJECT_NAME} ==========================") diff --git a/sbncode/CAFMaker/CAFMakerParams.h b/sbncode/CAFMaker/CAFMakerParams.h index 0250e740f..6843b4ad9 100644 --- a/sbncode/CAFMaker/CAFMakerParams.h +++ b/sbncode/CAFMaker/CAFMakerParams.h @@ -635,6 +635,11 @@ namespace caf "cvn" }; + Atom fBlipTag { + Name("BlipTag"), + Comment("Provide a string to label the blip input"), + "blipreco" + }; }; } diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index 9f2c01590..46c1d475c 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/sbncode/CAFMaker/CAFMaker_module.cc @@ -20,6 +20,7 @@ #include "sbncode/CAFMaker/FillExposure.h" #include "sbncode/CAFMaker/FillTrigger.h" #include "sbncode/CAFMaker/Utils.h" +#include "sbncode/CAFMaker/FillBlip.h" // C/C++ includes #include @@ -1928,6 +1929,13 @@ void CAFMaker::produce(art::Event& evt) noexcept { } } + //Fill blips. art::handle for blips and then call fill blips for each one. Make a vector to hold all of them. I handle for loop in Fill blip + art::Handle> blipHandle; + std::vector srblips; + if(evt.getByLabel( fParams.fBlipTag(), blipHandle)) //fill SR blips + { + FillBlip( *blipHandle, srblips); + } // collect the TPC slices std::vector> slices; std::vector slice_tag_suffixes; @@ -2605,6 +2613,7 @@ void CAFMaker::produce(art::Event& evt) noexcept { rec.sbnd_crt_veto = srsbndcrtveto; rec.opflashes = srflashes; rec.nopflashes = srflashes.size(); + rec.blips = srblips; rec.sbnd_frames = srsbndframeshiftinfo; rec.sbnd_timings = srsbndtiminginfo; rec.soft_trig = srsbndsofttrig; diff --git a/sbncode/CAFMaker/CMakeLists.txt b/sbncode/CAFMaker/CMakeLists.txt index 5dd4f2a8f..7ffa9eb93 100644 --- a/sbncode/CAFMaker/CMakeLists.txt +++ b/sbncode/CAFMaker/CMakeLists.txt @@ -39,6 +39,7 @@ art_make_library( LIBRARY_NAME sbncode_CAFMaker sbnobj::Common_Trigger sbnobj::SBND_CRT sbnobj::SBND_Timing + sbnobj::SBND_Blip lardataalg::DetectorInfo art::Framework_Services_System_TriggerNamesService_service sbncode_Metadata_MetadataSBN_service diff --git a/sbncode/CAFMaker/FillBlip.cxx b/sbncode/CAFMaker/FillBlip.cxx new file mode 100644 index 000000000..5d71edb35 --- /dev/null +++ b/sbncode/CAFMaker/FillBlip.cxx @@ -0,0 +1,96 @@ +#include "sbncode/CAFMaker/FillBlip.h" + +namespace caf +{ + void FillBlip( const std::vector& LAr_Blips, std::vector& CAF_Blips) + { + for(blip::Blip const& CurrentBlip: LAr_Blips) + { + caf::SRBlip NewBlip; + NewBlip.ID = CurrentBlip.ID; + NewBlip.isValid = CurrentBlip.isValid; + NewBlip.cryostat = CurrentBlip.Cryostat; + NewBlip.TPC = CurrentBlip.TPC; + NewBlip.nPlanes = CurrentBlip.NPlanes; + NewBlip.maxWireSpan = CurrentBlip.MaxWireSpan; + NewBlip.timeTick = CurrentBlip.TimeTick; + NewBlip.time = CurrentBlip.Time; + NewBlip.charge = CurrentBlip.Charge; + NewBlip.energy = CurrentBlip.Energy/1000.; //convert to GeV + NewBlip.energyESTAR = CurrentBlip.EnergyESTAR/1000.; //convert to GeV + NewBlip.energyPSTAR = CurrentBlip.EnergyPSTAR/1000.; //convert to GeV + NewBlip.proxTrkDist = CurrentBlip.ProxTrkDist; + NewBlip.proxTrkID = CurrentBlip.ProxTrkID; + NewBlip.inCylinder = CurrentBlip.inCylinder; + NewBlip.position.SetXYZ(CurrentBlip.Position.X(), CurrentBlip.Position.Y(), CurrentBlip.Position.Z()); + NewBlip.sigmaYZ = CurrentBlip.SigmaYZ; + NewBlip.dX = CurrentBlip.dX; + NewBlip.dYZ = CurrentBlip.dYZ; + if(CurrentBlip.truth.ID >= 0 ) //MC Blip + { + FillMCTruthBlip( CurrentBlip.truth, NewBlip.truthBlip ); + } + for(int iPlane=0; iPlane + +namespace caf +{ + + void FillBlip( const std::vector& LAr_Blips, std::vector& CAF_Blips); + void FillMCTruthBlip( blip::TrueBlip const & TrueLAr_Blip, caf::SRTrueBlip &TrueCAF_Blip ); + void FillBlipRelatedHitCluster(blip::HitClust const & LAr_HitClust, caf::SRBlipHitClust &CAF_HitClust); +} + +#endif diff --git a/sbncode/SinglePhotonAnalysis/CCDeltaRadiative_module.cc b/sbncode/SinglePhotonAnalysis/CCDeltaRadiative_module.cc new file mode 100644 index 000000000..54ee42f59 --- /dev/null +++ b/sbncode/SinglePhotonAnalysis/CCDeltaRadiative_module.cc @@ -0,0 +1,212 @@ +//////////////////////////////////////////////////////////////////////// +// Class: CCDeltaRadiative +// Plugin Type: filter (art v2_05_00) +// File: CCDeltaRadiative_module.cc +// +// Selects CC events with a photon from a Delta radiative decay. +// Mirrors NCDeltaRadiative_module.cc but requires CCNC == 0 (charged +// current) and checks all four Delta charge states (1114, 2114, 2214, 2224). +//////////////////////////////////////////////////////////////////////// + +#include "art/Framework/Core/EDFilter.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/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" + +#include "art_root_io/TFileService.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "art_root_io/TFileDirectory.h" + +#include + +#include "nusimdata/SimulationBase/MCTruth.h" + +#include "TTree.h" + +class CCDeltaRadiative : public art::EDFilter { + + TTree * ftree; + + int frun; + int fsubrun; + int fevent; + int fnu_pdg; + int fccnc; + int fmode; + int finteraction_type; + int fis_cc_delta_radiative; + int fparent_status_code; + int fparent_pdg; + +public: + explicit CCDeltaRadiative(fhicl::ParameterSet const & p); + + CCDeltaRadiative(CCDeltaRadiative const &) = delete; + CCDeltaRadiative(CCDeltaRadiative &&) = delete; + CCDeltaRadiative & operator = (CCDeltaRadiative const &) = delete; + CCDeltaRadiative & operator = (CCDeltaRadiative &&) = delete; + + void cout_stuff(art::Event & e, bool passed); + void FillTree(art::Event & e, + size_t const mct_index, + size_t const parent_index, + bool const is_cc_delta_radiative); + void Reset(); + bool filter(art::Event & e) override; + +}; + + +CCDeltaRadiative::CCDeltaRadiative(fhicl::ParameterSet const & p) : + art::EDFilter(p), + ftree(nullptr) { + + art::ServiceHandle tfs; + ftree = tfs->make("CCDeltaRadiativeFilter", ""); + + ftree->Branch("run", &frun, "run/I"); + ftree->Branch("subrun", &fsubrun, "subrun/I"); + ftree->Branch("event", &fevent, "event/I"); + ftree->Branch("nu_pdg", &fnu_pdg, "nu_pdg/I"); + ftree->Branch("ccnc", &fccnc, "ccnc/I"); + ftree->Branch("mode", &fmode, "mode/I"); + ftree->Branch("is_cc_delta_radiative", &fis_cc_delta_radiative, "is_cc_delta_radiative/I"); + ftree->Branch("parent_status_code", &fparent_status_code, "parent_status_code/I"); + ftree->Branch("parent_pdg", &fparent_pdg, "parent_pdg/I"); + +} + + +void CCDeltaRadiative::cout_stuff(art::Event & e, bool passed = false) { + + art::ValidHandle> const & ev_mct = + e.getValidHandle>("generator"); + + std::cout << passed << "\n" + << "==========================\n"; + for(simb::MCTruth const & mct : *ev_mct) { + std::cout << "----------------------------\n"; + for(int i = 0; i < mct.NParticles(); ++i) { + simb::MCParticle const & mcp = mct.GetParticle(i); + std::cout << mcp.TrackId() << " " << mcp.PdgCode() << " " << mcp.Mother() << " " << mcp.StatusCode() << "\n"; + } + } + +} + + +void CCDeltaRadiative::Reset() { + + frun = -1; + fsubrun = -1; + fevent = -1; + fnu_pdg = 0; + fccnc = -1; + fmode = -2; + finteraction_type = -2; + fis_cc_delta_radiative = -1; + fparent_status_code = -1; + fparent_pdg = 0; + +} + + +void CCDeltaRadiative::FillTree(art::Event & e, + size_t const mct_index, + size_t const parent_index, + bool const is_cc_delta_radiative) { + + Reset(); + + frun = e.id().run(); + fsubrun = e.id().subRun(); + fevent = e.id().event(); + + art::ValidHandle> const & ev_mct = + e.getValidHandle>("generator"); + simb::MCNeutrino const & mcn = ev_mct->at(mct_index).GetNeutrino(); + + fnu_pdg = mcn.Nu().PdgCode(); + fccnc = mcn.CCNC(); + fmode = mcn.Mode(); + finteraction_type = mcn.InteractionType(); + + fis_cc_delta_radiative = is_cc_delta_radiative; + if(parent_index != SIZE_MAX) { + fparent_status_code = ev_mct->at(mct_index).GetParticle(parent_index).StatusCode(); + fparent_pdg = ev_mct->at(mct_index).GetParticle(parent_index).PdgCode(); + } + + ftree->Fill(); + +} + + +// Returns true if |pdg| matches any Delta charge state (Δ-, Δ0, Δ+, Δ++). +static bool isDelta(int pdg) { + int apdg = std::abs(pdg); + return apdg == 1114 || apdg == 2114 || apdg == 2214 || apdg == 2224; +} + + +bool CCDeltaRadiative::filter(art::Event & e) { + + art::ValidHandle> const & ev_mct = + e.getValidHandle>("generator"); + + for(size_t i = 0; i < ev_mct->size(); ++i) { + + simb::MCTruth const & mct = ev_mct->at(i); + // Require charged-current interaction (CCNC == 0). + if(mct.GetNeutrino().CCNC() != 0) continue; + + std::vector exiting_photon_parents; + for(int j = 0; j < mct.NParticles(); ++j) { + simb::MCParticle const & mcp = mct.GetParticle(j); + if(mcp.TrackId() != j) { + std::cout << "ERROR: " << __LINE__ << " " << __PRETTY_FUNCTION__ << "\nTrackId does not match index\n"; + return false; + } + if(!(mcp.StatusCode() == 1 && mcp.PdgCode() == 22)) continue; + exiting_photon_parents.push_back(mcp.Mother()); + } + + std::vector in_nucleus_photons; + for(size_t const s : exiting_photon_parents) { + simb::MCParticle const & mcp = mct.GetParticle(s); + // Require parent vertex inside TPC. + if(std::abs(mcp.Vx()) > 210 || std::abs(mcp.Vy()) > 210 || mcp.Vz() > 510 || mcp.Vz() < -1) { + std::cout << "OUTSIDE TPC x y z =" << mcp.Vx() << " " << mcp.Vy() << " " << mcp.Vz() << std::endl; + return false; + } + if(isDelta(mcp.PdgCode())) { + if(ftree) FillTree(e, i, s, true); + return true; + } + else if(mcp.PdgCode() == 22) { + in_nucleus_photons.push_back(mcp.Mother()); + } + } + + for(size_t const s : in_nucleus_photons) { + simb::MCParticle const & mcp = mct.GetParticle(s); + if(isDelta(mcp.PdgCode())) { + if(ftree) FillTree(e, i, s, true); + return true; + } + } + + } + + if(ftree) FillTree(e, 0, SIZE_MAX, false); + return false; + +} + + +DEFINE_ART_MODULE(CCDeltaRadiative) diff --git a/sbncode/SinglePhotonAnalysis/CMakeLists.txt b/sbncode/SinglePhotonAnalysis/CMakeLists.txt index 35f24c2bb..4d144d72c 100644 --- a/sbncode/SinglePhotonAnalysis/CMakeLists.txt +++ b/sbncode/SinglePhotonAnalysis/CMakeLists.txt @@ -23,6 +23,13 @@ art_make_library( BASENAME_ONLY ROOT::Core ) +cet_build_plugin(CCDeltaRadiative art::module BASENAME_ONLY LIBRARIES + sbncode_SinglePhotonAnalysis_Libraries + sbncode_SinglePhotonAnalysis_SEAview + art_root_io::TFileService_service + ROOT::Tree + ) + cet_build_plugin(NCDeltaRadiative art::module BASENAME_ONLY LIBRARIES sbncode_SinglePhotonAnalysis_Libraries sbncode_SinglePhotonAnalysis_SEAview diff --git a/sbncode/SinglePhotonAnalysis/jobs/prodgenie_nu_singleinteraction_tpc_CCDeltaRadiative_filtered_sbnd.fcl b/sbncode/SinglePhotonAnalysis/jobs/prodgenie_nu_singleinteraction_tpc_CCDeltaRadiative_filtered_sbnd.fcl new file mode 100644 index 000000000..71178e589 --- /dev/null +++ b/sbncode/SinglePhotonAnalysis/jobs/prodgenie_nu_singleinteraction_tpc_CCDeltaRadiative_filtered_sbnd.fcl @@ -0,0 +1,106 @@ +# Simulates GENIE neutrino interactions from the BNB beam +# forcing one interaction per event, inside the TPC volume +# (with 10 cm padding on each side), +# selecting only CC events with photons coming out from the Delta + + +#include "prodgenie_nu_singleinteraction_tpc_sbnd.fcl" + +#------ this could be a separated file +# +# services +# + +# +# modules +# + + + +process_name: GenieGenFiltered + + +# +# services +# +services: { + + TFileService: { fileName: "hists_genie.root" } + IFDH: {} # required by GENIEGen + @table::sbnd_basic_services # from simulationservices_sbnd.fcl + @table::sbnd_random_services # from simulationservices_sbnd.fcl + FileCatalogMetadata: @local::sbnd_file_catalog_mc # from sam_sbnd.fcl + + # since this is a configuration expected to be run for production, + # we set up message configuration accordingly: + message: @local::sbnd_message_services_prod + +} # services + + +# +# input +# + + +# +# processing +# +physics: +{ + + producers: + { + generator: @local::sbnd_genie_simple + rns: { module_type: "RandomNumberSaver" } + } + + filters: + { + CCDeltaRadFilter: + { + module_type: "CCDeltaRadiative" + } + + } + + + #define the producer and filter modules for this path, order matters, + simulate: [ rns, generator, CCDeltaRadFilter] + + #define the output stream, there could be more than one if using filters + stream1: [ out1 ] + + #trigger_paths is a keyword and contains the paths that modify the art::event, + #ie filters and producers + trigger_paths: [simulate] + + #end_paths is a keyword and contains the paths that do not modify the art::Event, + #ie analyzers and output streams. these all run simultaneously + end_paths: [stream1] + +} # physics + + +# +# output +# +outputs: +{ + out1: + { + module_type: RootOutput + fileName: "prodgenie_bnb_nu_filtered_CCDeltaRad_sbnd_%p-%tc.root" # default file name, can override from command line with -o or --output + SelectEvents: [simulate] + dataTier: "generated" + compressionLevel: 1 + } +} # outputs + + +# +# override +# THIS DOES NOT WORK, CHECK! +physics.producers.generator.TopVolume: "volTPCActive" +physics.producers.generator.BeamName: "booster" +physics.producers.generator.EventGeneratorList: "CCRES" diff --git a/ups/product_deps b/ups/product_deps index 483437ae1..07db71db3 100644 --- a/ups/product_deps +++ b/ups/product_deps @@ -255,7 +255,7 @@ product version qual flags genie_xsec v3_06_02_sbn2 - larcv2 v2_2_6 - larsoft v10_14_02_02 - -sbnalg v10_14_02_01 - +sbnalg v10_14_02_04 - sbndaq_artdaq_core v1_10_06 - sbndata v01_08 - systematicstools v01_04_04 -