Skip to content

Commit 41d2c59

Browse files
seokjuchungseokjuchung
committed
CC Delta Rad filter
1 parent 8442cd7 commit 41d2c59

2 files changed

Lines changed: 318 additions & 0 deletions

File tree

‎sbncode/SinglePhotonAnalysis/CCDeltaRadiative_module.cc‎

Lines changed: 212 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,212 @@
1+
////////////////////////////////////////////////////////////////////////
2+
// Class: CCDeltaRadiative
3+
// Plugin Type: filter (art v2_05_00)
4+
// File: CCDeltaRadiative_module.cc
5+
//
6+
// Selects CC events with a photon from a Delta radiative decay.
7+
// Mirrors NCDeltaRadiative_module.cc but requires CCNC == 0 (charged
8+
// current) and checks all four Delta charge states (1114, 2114, 2214, 2224).
9+
////////////////////////////////////////////////////////////////////////
10+
11+
#include "art/Framework/Core/EDFilter.h"
12+
#include "art/Framework/Core/ModuleMacros.h"
13+
#include "art/Framework/Principal/Event.h"
14+
#include "art/Framework/Principal/Handle.h"
15+
#include "art/Framework/Principal/Run.h"
16+
#include "art/Framework/Principal/SubRun.h"
17+
#include "canvas/Utilities/InputTag.h"
18+
#include "fhiclcpp/ParameterSet.h"
19+
#include "messagefacility/MessageLogger/MessageLogger.h"
20+
21+
#include "art_root_io/TFileService.h"
22+
#include "art/Framework/Services/Registry/ServiceHandle.h"
23+
#include "art_root_io/TFileDirectory.h"
24+
25+
#include <memory>
26+
27+
#include "nusimdata/SimulationBase/MCTruth.h"
28+
29+
#include "TTree.h"
30+
31+
class CCDeltaRadiative : public art::EDFilter {
32+
33+
TTree * ftree;
34+
35+
int frun;
36+
int fsubrun;
37+
int fevent;
38+
int fnu_pdg;
39+
int fccnc;
40+
int fmode;
41+
int finteraction_type;
42+
int fis_cc_delta_radiative;
43+
int fparent_status_code;
44+
int fparent_pdg;
45+
46+
public:
47+
explicit CCDeltaRadiative(fhicl::ParameterSet const & p);
48+
49+
CCDeltaRadiative(CCDeltaRadiative const &) = delete;
50+
CCDeltaRadiative(CCDeltaRadiative &&) = delete;
51+
CCDeltaRadiative & operator = (CCDeltaRadiative const &) = delete;
52+
CCDeltaRadiative & operator = (CCDeltaRadiative &&) = delete;
53+
54+
void cout_stuff(art::Event & e, bool passed);
55+
void FillTree(art::Event & e,
56+
size_t const mct_index,
57+
size_t const parent_index,
58+
bool const is_cc_delta_radiative);
59+
void Reset();
60+
bool filter(art::Event & e) override;
61+
62+
};
63+
64+
65+
CCDeltaRadiative::CCDeltaRadiative(fhicl::ParameterSet const & p) :
66+
art::EDFilter(p),
67+
ftree(nullptr) {
68+
69+
art::ServiceHandle<art::TFileService> tfs;
70+
ftree = tfs->make<TTree>("CCDeltaRadiativeFilter", "");
71+
72+
ftree->Branch("run", &frun, "run/I");
73+
ftree->Branch("subrun", &fsubrun, "subrun/I");
74+
ftree->Branch("event", &fevent, "event/I");
75+
ftree->Branch("nu_pdg", &fnu_pdg, "nu_pdg/I");
76+
ftree->Branch("ccnc", &fccnc, "ccnc/I");
77+
ftree->Branch("mode", &fmode, "mode/I");
78+
ftree->Branch("is_cc_delta_radiative", &fis_cc_delta_radiative, "is_cc_delta_radiative/I");
79+
ftree->Branch("parent_status_code", &fparent_status_code, "parent_status_code/I");
80+
ftree->Branch("parent_pdg", &fparent_pdg, "parent_pdg/I");
81+
82+
}
83+
84+
85+
void CCDeltaRadiative::cout_stuff(art::Event & e, bool passed = false) {
86+
87+
art::ValidHandle<std::vector<simb::MCTruth>> const & ev_mct =
88+
e.getValidHandle<std::vector<simb::MCTruth>>("generator");
89+
90+
std::cout << passed << "\n"
91+
<< "==========================\n";
92+
for(simb::MCTruth const & mct : *ev_mct) {
93+
std::cout << "----------------------------\n";
94+
for(int i = 0; i < mct.NParticles(); ++i) {
95+
simb::MCParticle const & mcp = mct.GetParticle(i);
96+
std::cout << mcp.TrackId() << " " << mcp.PdgCode() << " " << mcp.Mother() << " " << mcp.StatusCode() << "\n";
97+
}
98+
}
99+
100+
}
101+
102+
103+
void CCDeltaRadiative::Reset() {
104+
105+
frun = -1;
106+
fsubrun = -1;
107+
fevent = -1;
108+
fnu_pdg = 0;
109+
fccnc = -1;
110+
fmode = -2;
111+
finteraction_type = -2;
112+
fis_cc_delta_radiative = -1;
113+
fparent_status_code = -1;
114+
fparent_pdg = 0;
115+
116+
}
117+
118+
119+
void CCDeltaRadiative::FillTree(art::Event & e,
120+
size_t const mct_index,
121+
size_t const parent_index,
122+
bool const is_cc_delta_radiative) {
123+
124+
Reset();
125+
126+
frun = e.id().run();
127+
fsubrun = e.id().subRun();
128+
fevent = e.id().event();
129+
130+
art::ValidHandle<std::vector<simb::MCTruth>> const & ev_mct =
131+
e.getValidHandle<std::vector<simb::MCTruth>>("generator");
132+
simb::MCNeutrino const & mcn = ev_mct->at(mct_index).GetNeutrino();
133+
134+
fnu_pdg = mcn.Nu().PdgCode();
135+
fccnc = mcn.CCNC();
136+
fmode = mcn.Mode();
137+
finteraction_type = mcn.InteractionType();
138+
139+
fis_cc_delta_radiative = is_cc_delta_radiative;
140+
if(parent_index != SIZE_MAX) {
141+
fparent_status_code = ev_mct->at(mct_index).GetParticle(parent_index).StatusCode();
142+
fparent_pdg = ev_mct->at(mct_index).GetParticle(parent_index).PdgCode();
143+
}
144+
145+
ftree->Fill();
146+
147+
}
148+
149+
150+
// Returns true if |pdg| matches any Delta charge state (Δ-, Δ0, Δ+, Δ++).
151+
static bool isDelta(int pdg) {
152+
int apdg = std::abs(pdg);
153+
return apdg == 1114 || apdg == 2114 || apdg == 2214 || apdg == 2224;
154+
}
155+
156+
157+
bool CCDeltaRadiative::filter(art::Event & e) {
158+
159+
art::ValidHandle<std::vector<simb::MCTruth>> const & ev_mct =
160+
e.getValidHandle<std::vector<simb::MCTruth>>("generator");
161+
162+
for(size_t i = 0; i < ev_mct->size(); ++i) {
163+
164+
simb::MCTruth const & mct = ev_mct->at(i);
165+
// Require charged-current interaction (CCNC == 0).
166+
if(mct.GetNeutrino().CCNC() != 0) continue;
167+
168+
std::vector<size_t> exiting_photon_parents;
169+
for(int j = 0; j < mct.NParticles(); ++j) {
170+
simb::MCParticle const & mcp = mct.GetParticle(j);
171+
if(mcp.TrackId() != j) {
172+
std::cout << "ERROR: " << __LINE__ << " " << __PRETTY_FUNCTION__ << "\nTrackId does not match index\n";
173+
return false;
174+
}
175+
if(!(mcp.StatusCode() == 1 && mcp.PdgCode() == 22)) continue;
176+
exiting_photon_parents.push_back(mcp.Mother());
177+
}
178+
179+
std::vector<size_t> in_nucleus_photons;
180+
for(size_t const s : exiting_photon_parents) {
181+
simb::MCParticle const & mcp = mct.GetParticle(s);
182+
// Require parent vertex inside TPC.
183+
if(std::abs(mcp.Vx()) > 210 || std::abs(mcp.Vy()) > 210 || mcp.Vz() > 510 || mcp.Vz() < -1) {
184+
std::cout << "OUTSIDE TPC x y z =" << mcp.Vx() << " " << mcp.Vy() << " " << mcp.Vz() << std::endl;
185+
return false;
186+
}
187+
if(isDelta(mcp.PdgCode())) {
188+
if(ftree) FillTree(e, i, s, true);
189+
return true;
190+
}
191+
else if(mcp.PdgCode() == 22) {
192+
in_nucleus_photons.push_back(mcp.Mother());
193+
}
194+
}
195+
196+
for(size_t const s : in_nucleus_photons) {
197+
simb::MCParticle const & mcp = mct.GetParticle(s);
198+
if(isDelta(mcp.PdgCode())) {
199+
if(ftree) FillTree(e, i, s, true);
200+
return true;
201+
}
202+
}
203+
204+
}
205+
206+
if(ftree) FillTree(e, 0, SIZE_MAX, false);
207+
return false;
208+
209+
}
210+
211+
212+
DEFINE_ART_MODULE(CCDeltaRadiative)

‎sbncode/SinglePhotonAnalysis/jobs/prodgenie_nu_singleinteraction_tpc_CCDeltaRadiative_filtered_sbnd.fcl‎

Lines changed: 106 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,106 @@
1+
# Simulates GENIE neutrino interactions from the BNB beam
2+
# forcing one interaction per event, inside the TPC volume
3+
# (with 10 cm padding on each side),
4+
# selecting only CC events with photons coming out from the Delta
5+
6+
7+
#include "prodgenie_nu_singleinteraction_tpc_sbnd.fcl"
8+
9+
#------ this could be a separated file
10+
#
11+
# services
12+
#
13+
14+
#
15+
# modules
16+
#
17+
18+
19+
20+
process_name: GenieGenFiltered
21+
22+
23+
#
24+
# services
25+
#
26+
services: {
27+
28+
TFileService: { fileName: "hists_genie.root" }
29+
IFDH: {} # required by GENIEGen
30+
@table::sbnd_basic_services # from simulationservices_sbnd.fcl
31+
@table::sbnd_random_services # from simulationservices_sbnd.fcl
32+
FileCatalogMetadata: @local::sbnd_file_catalog_mc # from sam_sbnd.fcl
33+
34+
# since this is a configuration expected to be run for production,
35+
# we set up message configuration accordingly:
36+
message: @local::sbnd_message_services_prod
37+
38+
} # services
39+
40+
41+
#
42+
# input
43+
#
44+
45+
46+
#
47+
# processing
48+
#
49+
physics:
50+
{
51+
52+
producers:
53+
{
54+
generator: @local::sbnd_genie_simple
55+
rns: { module_type: "RandomNumberSaver" }
56+
}
57+
58+
filters:
59+
{
60+
CCDeltaRadFilter:
61+
{
62+
module_type: "CCDeltaRadiative"
63+
}
64+
65+
}
66+
67+
68+
#define the producer and filter modules for this path, order matters,
69+
simulate: [ rns, generator, CCDeltaRadFilter]
70+
71+
#define the output stream, there could be more than one if using filters
72+
stream1: [ out1 ]
73+
74+
#trigger_paths is a keyword and contains the paths that modify the art::event,
75+
#ie filters and producers
76+
trigger_paths: [simulate]
77+
78+
#end_paths is a keyword and contains the paths that do not modify the art::Event,
79+
#ie analyzers and output streams. these all run simultaneously
80+
end_paths: [stream1]
81+
82+
} # physics
83+
84+
85+
#
86+
# output
87+
#
88+
outputs:
89+
{
90+
out1:
91+
{
92+
module_type: RootOutput
93+
fileName: "prodgenie_bnb_nu_filtered_CCDeltaRad_sbnd_%p-%tc.root" # default file name, can override from command line with -o or --output
94+
SelectEvents: [simulate]
95+
dataTier: "generated"
96+
compressionLevel: 1
97+
}
98+
} # outputs
99+
100+
101+
#
102+
# override
103+
# THIS DOES NOT WORK, CHECK!
104+
physics.producers.generator.TopVolume: "volTPCActive"
105+
physics.producers.generator.BeamName: "booster"
106+
physics.producers.generator.EventGeneratorList: "CCRES"

0 commit comments

Comments
 (0)