FillLSTME_module.cc
Go to the documentation of this file.
10 
11 #include "Geometry/Geometry.h"
12 #include "NovaDAQConventions/DAQConventions.h"
13 #include "RecoBase/Cluster.h"
14 #include "RecoBase/Vertex.h"
15 #include "RecoBase/FilterList.h"
16 #include "SummaryData/SpillData.h"
17 #include "Utilities/AssociationUtil.h"
19 
20 
21 #include <memory>
22 
29 
30 namespace LSTME
31 {
32 
33 class FillLSTME : public art::EDProducer
34 {
35 
36 public:
37  explicit FillLSTME(const fhicl::ParameterSet &pset);
38 
39  FillLSTME(FillLSTME const &) = delete;
40  FillLSTME(FillLSTME &&) = delete;
41  FillLSTME & operator = (FillLSTME const &) = delete;
42  FillLSTME & operator = (FillLSTME &&) = delete;
43 
44  void produce(art::Event &evt) override;
45 
46  void reconfigure(const fhicl::ParameterSet &pset);
47 
48  void fillEnergies(
49  art::Event &evt,
50  std::unique_ptr<std::vector<LSTMEnergy>> &sliceEnergies,
51  std::unique_ptr<art::Assns<LSTMEnergy, rb::Cluster>> &sliceEnergyAssoc
52  );
53 
54 private:
56 
59 
64 
66  std::vector<std::string> fPreselectionLabels;
67 
69 
70  bool getIsRHC(const art::Event &evt);
71  Model& selectModel(const art::Event &evt);
72 
74  Model &model,
75  const art::Event &evt,
76  const art::Handle<std::vector<rb::Cluster>> &slices_h,
77  size_t sliceIdx,
78  const std::vector<art::Ptr<rb::Vertex>> &elastics,
79  const std::vector<bool> &sliceIsMC
80  );
81 
82 };
83 
84 
86  : modelFDFHC(util::EnvExpansion(pset.get<std::string>("modelFDFHC"))),
87  modelFDRHC(util::EnvExpansion(pset.get<std::string>("modelFDRHC"))),
88  modelNDFHC(util::EnvExpansion(pset.get<std::string>("modelNDFHC"))),
89  modelNDRHC(util::EnvExpansion(pset.get<std::string>("modelNDRHC"))),
90  fObeyPreselection (pset.get<bool> ("ObeyPreselection" )),
91  fPreselectionLabels(pset.get<std::vector<std::string>> ("PreselectionLabels")),
92  useOppositeHornCurrentNetwork(pset.get<bool>("useOppositeHornCurrentNetwork"))
93 {
94  reconfigure(pset);
95 
96  produces<std::vector<LSTMEnergy>>();
97  produces<art::Assns<LSTMEnergy, rb::Cluster>>();
98 }
99 
101 {
102  beamLabel = pset.get<std::string>("beamLabel");
103  generatorLabel = pset.get<std::string>("generatorLabel");
104 
105  config.parse(pset);
106 }
107 
109 {
110  /* TODO: ask if this should be cached per file/run/job/etc ? */
111 
112  /*
113  * NOTE: Won't work for CRY files.
114  * Apparently need to change config file for them and set
115  * generatorLabel = beamLabel
116  */
118  if (! evt.isRealData()) {
119  evt.getByLabel(generatorLabel, spillPot);
120  }
121  else {
122  evt.getByLabel(beamLabel, spillPot);
123  }
124 
125  if (spillPot.failedToGet())
126  {
127  mf::LogError("FillLSTME") <<
128  "Spill Data not found, aborting without horn current information";
129  abort();
130  }
131 
132  if(useOppositeHornCurrentNetwork) return !spillPot->isRHC;
133 
134  // NB - the logic here will cause 0HC to use the FHC network
135  return spillPot->isRHC;
136 }
137 
139 {
141  novadaq::cnv::DetId detId = geom->DetId();
142 
143  bool isRHC = getIsRHC(evt);
144 
145  if (detId == novadaq::cnv::kFARDET)
146  {
147  if (isRHC) {
148  return modelFDRHC;
149  }
150  else {
151  return modelFDFHC;
152  }
153  }
154 
155  if (detId == novadaq::cnv::kNEARDET)
156  {
157  if (isRHC) {
158  return modelNDRHC;
159  }
160  else {
161  return modelNDFHC;
162  }
163  }
164 
165  mf::LogError("FillLSTME") << "Unsupported detector type, aborting";
166  abort();
167 }
168 
170  Model &model,
171  const art::Event &evt,
172  const art::Handle<std::vector<rb::Cluster>> &slices_h,
173  size_t sliceIdx,
174  const std::vector<art::Ptr<rb::Vertex>> &elastics,
175  const std::vector<bool> &sliceIsMC
176 )
177 {
178  VarDictBuilder builder(
179  evt, slices_h, sliceIdx, elastics, config
180  );
181  VarDict varDict = builder.build();
182 
183  /* TODO: this must be removed once geometry is fixed */
185  TrackLengthCorrection trkLenCorr(
186  sliceIsMC[sliceIdx], (geom->DetId() == novadaq::cnv::kFARDET)
187  );
188  trkLenCorr.shift(varDict);
189 
190  return model.predict(varDict);
191 }
192 
194  art::Event &evt,
195  std::unique_ptr<std::vector<LSTMEnergy>> &sliceEnergies,
196  std::unique_ptr<art::Assns<LSTMEnergy, rb::Cluster>> &sliceEnergyAssoc
197 )
198 {
200  if (! evt.getByLabel(config.sliceLabel, slices_h)) {
201  mf::LogError("FillLSTME") << "Failed to get slice info" << std::endl;
202  return;
203  }
204 
205  art::FindManyP<rb::Vertex> elasticVxt_fmp(
206  slices_h, evt, config.elasticArmsLabel
207  );
208  if (! elasticVxt_fmp.isValid()) {
209  mf::LogError("FillLSTME") << "No elastics found" << std::endl;
210  return;
211  }
212 
213  Model &model = selectModel(evt);
214  std::vector<bool> sliceIsMC = getSliceIsMC(*slices_h);
215 
216  for (size_t sliceIdx = 0; sliceIdx < slices_h->size(); ++sliceIdx)
217  {
218  const art::Ptr<rb::Cluster> slicePtr(slices_h, sliceIdx);
219 
220  if (slicePtr->IsNoise()) {
221  continue;
222  }
223 
224  if(fObeyPreselection && rb::IsFiltered(evt, slices_h, sliceIdx, fPreselectionLabels)) continue;
225 
226  auto elastics = elasticVxt_fmp.at(sliceIdx);
227 
229  model, evt, slices_h, sliceIdx, elastics, sliceIsMC
230  );
231 
232  sliceEnergies->emplace_back(std::move(energy));
233 
235  *this, evt,
236  *(sliceEnergies.get()), slicePtr, *(sliceEnergyAssoc.get())
237  );
238  }
239 }
240 
242 {
243  std::unique_ptr<std::vector<LSTMEnergy>> sliceEnergies(
244  new std::vector<LSTMEnergy>
245  );
246  std::unique_ptr<art::Assns<LSTMEnergy, rb::Cluster>> sliceEnergyAssoc(
248  );
249 
250  fillEnergies(evt, sliceEnergies, sliceEnergyAssoc);
251 
252  evt.put(std::move(sliceEnergies));
253  evt.put(std::move(sliceEnergyAssoc));
254 }
255 
257 
258 }
bool isRHC
is the beam in antineutrino mode, aka RHC
Definition: SpillData.h:28
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
bool useOppositeHornCurrentNetwork
Filter events based on their run/event numbers.
FillLSTME(const fhicl::ParameterSet &pset)
void fillEnergies(art::Event &evt, std::unique_ptr< std::vector< LSTMEnergy >> &sliceEnergies, std::unique_ptr< art::Assns< LSTMEnergy, rb::Cluster >> &sliceEnergyAssoc)
void reconfigure(const fhicl::ParameterSet &pset)
void produce(art::Event &evt) override
std::vector< std::string > fPreselectionLabels
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Definition: EnvExpand.cxx:8
Definition: config.py:1
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
bool isRealData() const
Definition: Event.h:83
DEFINE_ART_MODULE(TestTMapFile)
LSTMEnergy predictEnergy(Model &model, const art::Event &evt, const art::Handle< std::vector< rb::Cluster >> &slices_h, size_t sliceIdx, const std::vector< art::Ptr< rb::Vertex >> &elastics, const std::vector< bool > &sliceIsMC)
LSTMEnergy predict(const VarDict &varDict)
Definition: Model.cxx:27
VarDict build() const
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
Far Detector at Ash River, MN.
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
int evt
double energy
Definition: plottest35.C:25
Near Detector in the NuMI cavern.
void shift(VarDict &varDict) const
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
Definition: FilterList.h:96
Model & selectModel(const art::Event &evt)
Vertex location in position and time.
std::string generatorLabel
FillLSTME & operator=(FillLSTME const &)=delete
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void geom(int which=0)
Definition: geom.C:163
Definition: VarDict.h:7
std::vector< bool > getSliceIsMC(const std::vector< rb::Cluster > &slices)
Definition: utils.cxx:6
std::string beamLabel
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
const XML_Char XML_Content * model
Definition: expat.h:151
Definition: fwd.h:28
Encapsulate the geometry of one entire detector (near, far, ndos)
bool getIsRHC(const art::Event &evt)
bool failedToGet() const
Definition: Handle.h:196
enum BeamMode string