runEventSelection.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include "cetlib_except/exception.h"
5 #include <ROOT/RDataFrame.hxx>
6 #include "TInterpreter.h"
7 
8 #include "SelectionDefinitions.h"
9 #include "TTree.h"
10 
11 class testing {
12  private:
13  int num;
14  public:
15  testing(int n) : num(n){}
16  int operator()() const {
17  //std::cout << num << std::endl;
18  return num;
19  }
20 };
21 
22 auto DefineMC(ROOT::RDataFrame & frame, const fhicl::ParameterSet & pset) {
23  testing testing1(1);
24  testing testing2(2);
25  auto mc = frame.Define("testing1", testing(1)/*testing1*/)
26  .Define("testing2", testing(2)/*testing2*/)
27  .Define("primary_isBeamType", isBeamType(pset.get<bool>("CheckCalo")),
28  {"reco_beam_type", "reco_beam_incidentEnergies"})
29  .Define("primary_ends_inAPA3",
30  endAPA3(pset.get<double>("EndZHigh")), {"reco_beam_endZ"})
31  .Define("shower_dists",
32  shower_dists(pset.get<double>("TrackScoreCut")),
33  {"reco_daughter_PFP_trackScore_collection",
34  "reco_daughter_allShower_startX",
35  "reco_daughter_allShower_startY",
36  "reco_daughter_allShower_startZ",
37  "reco_beam_endX", "reco_beam_endY", "reco_beam_endZ"})
38  .Define("has_shower_dist_energy",
39  has_shower_dist_energy(pset.get<double>("TrackScoreCut")),
40  {"reco_daughter_PFP_trackScore_collection",
41  "reco_daughter_allShower_startX",
42  "reco_daughter_allShower_startY",
43  "reco_daughter_allShower_startZ",
44  "reco_daughter_allShower_energy",
45  "reco_beam_endX", "reco_beam_endY", "reco_beam_endZ"})
46  .Define("reco_daughter_allTrack_truncLibo_dEdX_pos",
47  truncatedMean_pos(pset.get<double>("Limit")),
48  {"reco_daughter_allTrack_calibrated_dEdX_SCE"})
49  .Define("has_noPion_daughter",
51  pset.get<double>("TrackScoreCut"),
52  pset.get<double>("Chi2Cut"),
53  pset.get<double>("dEdXLow"),
54  pset.get<double>("dEdXMed"),
55  pset.get<double>("dEdXHigh")),
56  {"reco_daughter_PFP_trackScore_collection",
57  "reco_daughter_allTrack_ID",
58  "reco_daughter_allTrack_truncLibo_dEdX_pos",
59  "reco_daughter_allTrack_Chi2_proton",
60  "reco_daughter_allTrack_Chi2_ndof"})
61  .Define("new_interaction_topology",
62  new_interaction_topology(pset.get<double>("EndZLow"),
63  pset.get<double>("EndZHigh"),
64  pset.get<double>("Threshold")),
65  {"true_beam_PDG",
66  "true_beam_endZ", "true_beam_endProcess", "true_daughter_nPi0",
67  "true_beam_daughter_PDG", "true_beam_daughter_startP"})
68  .Define("beam_backtrack", backtrack_beam,
69  {"reco_beam_true_byHits_process", "reco_beam_true_byHits_matched",
70  "reco_beam_true_byHits_origin", "reco_beam_true_byHits_PDG"})
71  .Define("daughter_categories", categorize_daughters,
72  {"true_beam_ID", "reco_daughter_PFP_true_byHits_origin",
73  "reco_daughter_PFP_true_byHits_ID", "reco_daughter_PFP_true_byHits_PDG",
74  "reco_daughter_PFP_true_byHits_parID", "reco_daughter_PFP_true_byHits_parPDG",
75  "true_beam_daughter_ID",
76  "true_beam_grand_daughter_ID"})
77  .Define("leading_p_costheta", leading_p_costheta,
78  {"true_beam_endPx", "true_beam_endPy", "true_beam_endPz",
79  "true_beam_daughter_PDG", "true_beam_daughter_startPx",
80  "true_beam_daughter_startPy", "true_beam_daughter_startPz"})
81  .Define("leading_p_costheta2", leading_costheta(2212),
82  {"true_beam_endPx", "true_beam_endPy", "true_beam_endPz",
83  "true_beam_daughter_PDG", "true_beam_daughter_startPx",
84  "true_beam_daughter_startPy", "true_beam_daughter_startPz"})
85  .Define("leading_piplus_costheta", leading_costheta(211),
86  {"true_beam_endPx", "true_beam_endPy", "true_beam_endPz",
87  "true_beam_daughter_PDG", "true_beam_daughter_startPx",
88  "true_beam_daughter_startPy", "true_beam_daughter_startPz"})
89  .Define("leading_pi0_costheta", leading_costheta(111),
90  {"true_beam_endPx", "true_beam_endPy", "true_beam_endPz",
91  "true_beam_daughter_PDG", "true_beam_daughter_startPx",
92  "true_beam_daughter_startPy", "true_beam_daughter_startPz"})
93  .Define("daughter_PDGs_types", daughter_PDG_types,
94  {"reco_daughter_PFP_true_byHits_PDG"});
95 
96  if(pset.get<bool>("UseBI")) {
97  mc = mc.Define(
98  "passBeamCut",
99  beam_cut_BI(pset.get<double>("MCXLow"),
100  pset.get<double>("MCXHigh"),
101  pset.get<double>("MCYLow"),
102  pset.get<double>("MCYHigh"),
103  pset.get<double>("MCZLow"),
104  pset.get<double>("MCZHigh"),
105  pset.get<double>("MCCosLow")),
106  {"reco_beam_startX", "reco_beam_startY", "reco_beam_startZ",
107  "reco_beam_trackDirX", "reco_beam_trackDirY",
108  "reco_beam_trackDirZ", "beam_inst_X", "beam_inst_Y",
109  "beam_inst_dirX", "beam_inst_dirY", "beam_inst_dirZ"});
110  }
111  else {
112  mc = mc .Define(
113  "passBeamCut",
114  beam_cut_TPC(pset.get<bool>("DoAngle"),
115  pset.get<double>("XYZCut"),
116  pset.get<double>("CosCut"),
117  pset.get<double>("MCMeanX"),
118  pset.get<double>("MCMeanY"),
119  pset.get<double>("MCMeanZ"),
120  pset.get<double>("MCSigmaX"),
121  pset.get<double>("MCSigmaY"),
122  pset.get<double>("MCSigmaZ"),
123  pset.get<double>("MCMeanThetaX"),
124  pset.get<double>("MCMeanThetaY"),
125  pset.get<double>("MCMeanThetaZ")),
126  {"reco_beam_calo_startX", "reco_beam_calo_startY",
127  "reco_beam_calo_startZ", "reco_beam_calo_endX",
128  "reco_beam_calo_endY", "reco_beam_calo_endZ"});
129  }
130  mc = mc.Define("selection_ID", selection_ID,
131  {"primary_isBeamType", "primary_ends_inAPA3",
132  "has_noPion_daughter", "passBeamCut",
133  "has_shower_dist_energy"});
134  std::cout << "Filtering MC" << std::endl;
135  auto filtered = mc.Filter("true_beam_PDG == 211 || true_beam_PDG == -13");
136  return filtered;
137 }
138 
139 auto DefineData(ROOT::RDataFrame & frame, const fhicl::ParameterSet & pset) {
140  testing testing1(3);
141  testing testing2(4);
142  auto data = frame.Define("testing1", testing(1)/*testing1*/)
143  .Define("testing2", testing(2)/*testing2*/)
144  .Define("beamPID", data_beam_PID, {"beam_inst_PDG_candidates"})
145  .Define("passBeamQuality",
146  data_BI_quality(pset.get<bool>("DoNTracks")),
147  {"beam_inst_nMomenta", "beam_inst_nTracks"})
148  .Define("primary_isBeamType", isBeamType(pset.get<bool>("CheckCalo")),
149  {"reco_beam_type", "reco_beam_incidentEnergies"})
150  .Define("primary_ends_inAPA3",
151  endAPA3(pset.get<double>("EndZHigh")), {"reco_beam_endZ"})
152  .Define("shower_dists",
153  shower_dists(pset.get<double>("TrackScoreCut")),
154  {"reco_daughter_PFP_trackScore_collection",
155  "reco_daughter_allShower_startX",
156  "reco_daughter_allShower_startY",
157  "reco_daughter_allShower_startZ",
158  "reco_beam_endX", "reco_beam_endY", "reco_beam_endZ"})
159  .Define("has_shower_dist_energy",
160  has_shower_dist_energy(pset.get<double>("TrackScoreCut")),
161  {"reco_daughter_PFP_trackScore_collection",
162  "reco_daughter_allShower_startX",
163  "reco_daughter_allShower_startY",
164  "reco_daughter_allShower_startZ",
165  "reco_daughter_allShower_energy",
166  "reco_beam_endX", "reco_beam_endY", "reco_beam_endZ"})
167  .Define("reco_daughter_allTrack_truncLibo_dEdX_pos",
168  truncatedMean_pos(pset.get<double>("Limit")),
169  {"reco_daughter_allTrack_calibrated_dEdX_SCE"})
170  .Define("has_noPion_daughter",
172  pset.get<double>("TrackScoreCut"),
173  pset.get<double>("Chi2Cut"),
174  pset.get<double>("dEdXLow"),
175  pset.get<double>("dEdXMed"),
176  pset.get<double>("dEdXHigh")),
177  {"reco_daughter_PFP_trackScore_collection",
178  "reco_daughter_allTrack_ID",
179  "reco_daughter_allTrack_truncLibo_dEdX_pos",
180  "reco_daughter_allTrack_Chi2_proton",
181  "reco_daughter_allTrack_Chi2_ndof"});
182 
183  if(pset.get<bool>("UseBI")) {
184  data = data.Define(
185  "passBeamCut",
186  beam_cut_BI(pset.get<double>("DataXLow"),
187  pset.get<double>("DataXHigh"),
188  pset.get<double>("DataYLow"),
189  pset.get<double>("DataYHigh"),
190  pset.get<double>("DataZLow"),
191  pset.get<double>("DataZHigh"),
192  pset.get<double>("DataCosLow")),
193  {"reco_beam_startX", "reco_beam_startY", "reco_beam_startZ",
194  "reco_beam_trackDirX", "reco_beam_trackDirY",
195  "reco_beam_trackDirZ", "beam_inst_X", "beam_inst_Y",
196  "beam_inst_dirX", "beam_inst_dirY", "beam_inst_dirZ"});
197  }
198  else {
199  data = data .Define(
200  "passBeamCut",
201  beam_cut_TPC(pset.get<bool>("DoAngle"),
202  pset.get<double>("XYZCut"),
203  pset.get<double>("CosCut"),
204  pset.get<double>("DataMeanX"),
205  pset.get<double>("DataMeanY"),
206  pset.get<double>("DataMeanZ"),
207  pset.get<double>("DataSigmaX"),
208  pset.get<double>("DataSigmaY"),
209  pset.get<double>("DataSigmaZ"),
210  pset.get<double>("DataMeanThetaX"),
211  pset.get<double>("DataMeanThetaY"),
212  pset.get<double>("DataMeanThetaZ")),
213  {"reco_beam_calo_startX", "reco_beam_calo_startY",
214  "reco_beam_calo_startZ", "reco_beam_calo_endX",
215  "reco_beam_calo_endY", "reco_beam_calo_endZ"});
216  }
217  data = data.Define("selection_ID", selection_ID,
218  {"primary_isBeamType", "primary_ends_inAPA3",
219  "has_noPion_daughter", "passBeamCut",
220  "has_shower_dist_energy"});
221  auto filtered = data.Filter("beamPID == true");
222  return filtered;
223 }
224 
225 int main(int argc, char ** argv){
226  //gInterpreter->GenerateDictionary("vector<vector<int>>",
227  // "vector");
228 
229  std::string fcl_file;
231  std::string mc_file, data_file;
232  bool found_mc = false, found_data = false;
233  // Options to run
234  for (int iArg = 1; iArg < argc; iArg++) {
235  if (!strcasecmp(argv[iArg],"-c")) {
236  fcl_file = argv[++iArg];
237  }
238  if (!strcasecmp(argv[iArg],"-m")) {
239  mc_file = argv[++iArg];
240  found_mc = true;
241  }
242  if (!strcasecmp(argv[iArg],"-d")) {
243  data_file = argv[++iArg];
244  found_data = true;
245  }
246  if (!strcasecmp(argv[iArg],"-o")) {
247  output_file = argv[++iArg];
248  }
249  if (!strcasecmp(argv[iArg],"-h")) {
250  std::cout << "Usage: runPDSPThinSliceFit -c fclfile.fcl " <<
251  "-o outputfile.root " << std::endl;
252  return 1;
253  }
254  }
255 
256  char const* fhicl_env = getenv("FHICL_FILE_PATH");
257  std::string search_path;
258  if (fhicl_env == nullptr) {
259  std::cerr << "Expected environment variable FHICL_FILE_PATH is missing or empty: using \".\"\n";
260  search_path = ".";
261  }
262  else {
263  search_path = std::string{fhicl_env};
264  }
265 
266  cet::filepath_first_absolute_or_lookup_with_dot lookupPolicy{search_path};
267  auto const pset = fhicl::ParameterSet::make(fcl_file, lookupPolicy);
268  std::string tree_name = pset.get<std::string>("TreeName");
269 
270  if (pset.get<bool>("UseMT"))
271  ROOT::EnableImplicitMT(4);
272 
273  if (found_mc) {
274  ROOT::RDataFrame frame(tree_name, mc_file);
275  std::cout << "Calling DefineMC" << std::endl;
276  auto mc = DefineMC(frame, pset);
277  std::cout << "Snapshotting" << std::endl;
279  mc.Snapshot(tree_name, "eventSelection_mc_all.root");
281  std::cout << "Time: " <<
282  std::chrono::duration_cast<std::chrono::seconds>(time1 - time0).count() <<
283  std::endl;
284 
285  }
286 
287  if (found_data) {
288  ROOT::RDataFrame data_frame(tree_name, data_file);
289  std::cout << "Calling DefineData" << std::endl;
290  auto data = DefineData(data_frame, pset);
291  std::cout << "Snapshotting" << std::endl;
292  data.Snapshot(tree_name, "eventSelection_data_all.root");
293  data = data.Filter("passBeamQuality");
294  data.Snapshot(tree_name, "eventSelection_data_BeamQuality.root");
295  }
296 
297  return 0;
298 }
LArSoft test utilities.
std::string string
Definition: nybbler.cc:12
static ParameterSet make(intermediate_table const &tbl)
Definition: ParameterSet.cc:68
#define Define
Definition: scanner.cpp:11520
auto data_beam_PID
auto DefineData(ROOT::RDataFrame &frame, const fhicl::ParameterSet &pset)
auto leading_p_costheta
testing(int n)
second seconds
Alias for common language habits.
Definition: spacetime.h:88
std::void_t< T > n
int main(int argc, char **argv)
std::string getenv(std::string const &name)
Definition: getenv.cc:15
T get(std::string const &key) const
Definition: ParameterSet.h:271
auto backtrack_beam
int operator()() const
auto selection_ID
auto daughter_PDG_types(const std::vector< int > bt_PDGs)
auto categorize_daughters
auto DefineMC(ROOT::RDataFrame &frame, const fhicl::ParameterSet &pset)
QTextStream & endl(QTextStream &s)