RegCVNTF_module.cc
Go to the documentation of this file.
1 #include<iostream>
2 
10 #include "fhiclcpp/ParameterSet.h"
12 
13 #include "RecoBase/Cluster.h"
14 #include "RecoBase/FilterList.h"
15 #include "RecoBase/Track.h"
16 #include "SummaryData/SpillData.h"
17 #include "Utilities/AssociationUtil.h"
18 
19 
22 #include "CVN/func/PixelMap.h"
23 #include "CVN/func/RegResult.h"
24 #include "CVN/func/Result.h"
25 
28 #include <memory>
29 
30 
32 
33 namespace regcvntf{
34  class RegCVNTF : public art::EDProducer{
35  public:
36 
37  explicit RegCVNTF(fhicl::ParameterSet const &pset);
38  virtual ~RegCVNTF();
39 
40  void produce(art::Event& evt);
41 
42  std::vector< tensorflow::Tensor > vector_to_tensor_nue(std::vector<unsigned char>, unsigned int ncells, unsigned int nplanes);
43  tensorflow::Tensor vector_to_tensor_numu(std::vector<unsigned char>);
44  bool IsRHC(const art::Event &evt);
45 
46  std::pair<RegModel*, RegModel*> GetModelsNue(const art::Event &evt);
47  RegModel* GetModelsNumu(const art::Event &evt);
48  bool PassCVNCut(const cvn::Result& result);
49 
50  protected:
55 
60 
65  double fCVNCut;
66 
74 
75  std::pair<RegModel*, RegModel*> fTFFHC;
76  std::pair<RegModel*, RegModel*> fTFRHC;
77 
80 
81  };
82 }
83 
84 namespace regcvntf{
85 
87  fSliceLabel (pset.get<std::string>("SliceLabel")),
88  fPixelMapInput (pset.get<std::string>("PixelMapInput")),
89  fGeneratorLabel (pset.get<std::string>("GeneratorLabel")),
90  fNuMILabel (pset.get<std::string>("NuMILabel")),
91  fProngInput (pset.get<bool> ("ProngInput")),
92  fProngModLabel (pset.get<std::string> ("ProngModLabel")),
93  fProng3DLabel (pset.get<std::string> ("Prong3DLabel")),
94  fTrack3DLabel (pset.get<std::string> ("Track3DLabel")),
95  fCVNLabel1 (pset.get<std::string> ("CVNLabel1")),
96  fCVNLabel2 (pset.get<std::string> ("CVNLabel2")),
97  fCVNLabel3 (pset.get<std::string> ("CVNLabel3")),
98  fApplyLooseCVN (pset.get<bool> ("ApplyLooseCVN")),
99  fCVNCut (pset.get<double> ("CVNCut")),
100  fLibPath (pset.get<std::string> ("LibPath")),
101  fTFProtoBufEvtE (pset.get<std::string> ("TFProtoBufEvtE")),
102  fTFProtoBufShE (pset.get<std::string> ("TFProtoBufShE")),
103  fTFProtoBufRHCEvtE (pset.get<std::string> ("TFProtoBufRHCEvtE")),
104  fTFProtoBufRHCShE (pset.get<std::string> ("TFProtoBufRHCShE")),
105  fTFProtoBufHadE (pset.get<std::string> ("TFProtoBufHadE")),
106  fTFProtoBufRHCHadE (pset.get<std::string> ("TFProtoBufRHCHadE"))
107  {
113 
116 
117  fTFFHC = std::make_pair(new RegModel(fTFProtoBufEvtE),
118  new RegModel(fTFProtoBufShE));
119  fTFRHC = std::make_pair(new RegModel(fTFProtoBufRHCEvtE),
120  new RegModel(fTFProtoBufRHCShE));
121  fTFFHCHad = new RegModel(fTFProtoBufHadE);
122  fTFRHCHad = new RegModel(fTFProtoBufRHCHadE);
123 
124  produces< std::vector<cvn::RegNuResult> >();
125  produces< std::vector<cvn::RegProngResult> >();
126  produces< art::Assns<cvn::RegNuResult, rb::Cluster> >();
127  produces< art::Assns<cvn::RegProngResult, rb::Prong> >();
128  produces< std::vector<cvn::RegHadronResult> >();
129  produces< art::Assns<cvn::RegHadronResult, rb::Cluster> >();
130  }
131 
133  {
134  if(fTFFHC.first) delete fTFFHC.first;
135  if(fTFFHC.second) delete fTFFHC.second;
136  if(fTFRHC.first) delete fTFRHC.first;
137  if(fTFRHC.second) delete fTFRHC.second;
138  }
139 
141  {
143  if (!evt.isRealData())
144  evt.getByLabel(fGeneratorLabel, spillPot);
145  else
146  evt.getByLabel(fNuMILabel, spillPot);
147 
148  if (spillPot.failedToGet())
149  {
150  mf::LogError("RegCVNTF") <<
151  "Spill Data not found, aborting without horn current information";
152  abort();
153  }
154 
155  return spillPot->isRHC;
156  }
157 
158  std::pair<RegModel*, RegModel*> RegCVNTF::GetModelsNue(const art::Event &evt)
159  {
160  if (IsRHC(evt))
161  return fTFRHC;
162  else
163  return fTFFHC;
164  }
165 
167  {
168  if (IsRHC(evt))
169  return fTFRHCHad;
170  else
171  return fTFFHCHad;
172  }
173 
175  {
176  double nueid = result.fOutput[1];
177  // auto output = result.fOutput;
178  // double nueid = 0;
179  // for(int idx = 0; idx < caf::kNumCVNFinalStates; ++idx){
180  // const caf::CVNFinalState state = caf::cvnStates[idx];
181  // const double prob = result.fOutput[idx];
182  //
183  // if(state.nuPdg == 12) nueid += prob;
184  // }
185  // nueid += result.fOutput[caf::kCVN_nE_Other];
186 
187  return nueid > fCVNCut;
188  }
189 
190  std::vector< tensorflow::Tensor > RegCVNTF::vector_to_tensor_nue(std::vector<unsigned char> pm, unsigned int nplanes, unsigned int ncells)
191  {
192 
193  std::size_t const half_size = pm.size() / 2;
194  std::vector<unsigned char> pm_x(pm.begin(), pm.begin() + half_size);
195  std::vector<unsigned char> pm_y(pm.begin() + half_size, pm.end());
196 
197  long long int samples = 1, rows = nplanes, cols = ncells;
198 
199  std::vector< tensorflow::Tensor > _x;
200  for (unsigned int ii = 0; ii < 2; ++ii){
201  tensorflow::Tensor _xtemp(tensorflow::DT_FLOAT, tensorflow::TensorShape({ samples, rows, cols, 1 }));
202  _x.push_back(_xtemp);
203  }
204 
205  for (long long int s = 0; s < samples; ++s) {
206  for (long long int r = 0; r < rows; ++r) {
207  for (long long int c = 0; c < cols; ++c) {
208  unsigned int element = c + cols * r;
209  _x[0].tensor<float, 4>()(s, r, c, 0) = pm_x[element];
210  _x[1].tensor<float, 4>()(s, r, c, 0) = pm_y[element];
211  }
212  }
213  }
214  return _x;
215 
216  }
217 
218  tensorflow::Tensor RegCVNTF::vector_to_tensor_numu(std::vector<unsigned char> pm) {
219  const unsigned int vectorSize = pm.size();
220 
221  // Initialize the tensors
222  tensorflow::Tensor tensor(tensorflow::DT_FLOAT, {1, vectorSize});
223  auto rel = tensor.tensor<float,2>();
224 
225  // Loop over each element
226  for(unsigned int i = 0; i < vectorSize; ++i) rel(0, i) = pm[i];
227 
228  return tensor;
229  }
230 
232  {
233  std::unique_ptr< std::vector<cvn::RegNuResult> >
234  resultEvtCol(new std::vector<cvn::RegNuResult>);
235  std::unique_ptr< std::vector<cvn::RegProngResult> >
236  resultShCol(new std::vector<cvn::RegProngResult>);
237  std::unique_ptr< art::Assns<cvn::RegNuResult, rb::Cluster> >
239  std::unique_ptr< art::Assns<cvn::RegProngResult, rb::Prong> >
241  std::unique_ptr< std::vector<cvn::RegHadronResult> >
242  resultHadCol(new std::vector<cvn::RegHadronResult>);
243  std::unique_ptr< art::Assns<cvn::RegHadronResult, rb::Cluster> >
245 
246 
247  std::pair<RegModel*, RegModel*> models = GetModelsNue(evt);
248  RegModel* hadModel = GetModelsNumu(evt);
249 
250  // Get slices
252  evt.getByLabel(fSliceLabel, slicecol);
253  art::PtrVector<rb::Cluster> slicelist;
254  for(unsigned int i = 0; i < slicecol->size(); ++i){
255  slicelist.push_back(art::Ptr<rb::Cluster>(slicecol, i));
256  }
257 
258  // Get pixel maps
259  art::FindManyP<cvn::PixelMap> fmPixelMap(slicecol, evt, fPixelMapInput);
260  // Get cvn nue results
261  art::FindManyP<cvn::Result> fmCVN1(slicecol, evt, fCVNLabel1);
262  art::FindManyP<cvn::Result> fmCVN2(slicecol, evt, fCVNLabel2);
263  art::FindManyP<cvn::Result> fmCVN3(slicecol, evt, fCVNLabel3);
264 
265 
266  //loop over slices
267  for(size_t iClust = 0; iClust < slicelist.size(); ++iClust) {
268  if(!fmPixelMap.isValid()) continue;
269  if(slicelist[iClust]->IsNoise()) continue;
270 
271  const std::vector<art::Ptr<cvn::PixelMap> > pixelMaps = fmPixelMap.at(iClust);
272  if(pixelMaps.empty()) continue;
273 
274  if(fApplyLooseCVN && fmCVN1.isValid() && fmCVN2.isValid()){
275  cvn::Result cvnNue1;
276  cvn::Result cvnNue2;
277  cvn::Result cvnNue3;
278  const std::vector<art::Ptr<cvn::Result>> results1 = fmCVN1.at(iClust);
279  const std::vector<art::Ptr<cvn::Result>> results2 = fmCVN2.at(iClust);
280  const std::vector<art::Ptr<cvn::Result>> results3 = fmCVN3.at(iClust);
281  if ( !results1.empty() && !results2.empty() && !results3.empty() ){
282  cvnNue1 = *results1[0];
283  cvnNue2 = *results2[0];
284  cvnNue3 = *results3[0];
285  if (!PassCVNCut(cvnNue1) && !PassCVNCut(cvnNue2) &&!PassCVNCut(cvnNue3)) continue;
286  }
287  }
288 
289  std::vector<unsigned char> evtpm = (*pixelMaps[0]).PixelMapToVector(true);
290  std::vector<tensorflow::Tensor> tensorEvtE = vector_to_tensor_nue(evtpm, (*pixelMaps[0]).NPlanePerView(), (*pixelMaps[0]).NCell());
291 
292  std::vector<tensorflow::Tensor> resultEvtE = models.first->Predict({{"input_1",tensorEvtE[0]},{"input_2",tensorEvtE[1]}}, {"output_node0"});
293  auto tfoutputEvtE = resultEvtE[0].tensor<float,2>();
294 
295  resultEvtCol->emplace_back(tfoutputEvtE(0,0));
296  util::CreateAssn(*this, evt, *(resultEvtCol.get()),
297  slicelist[iClust], *(assoc.get()), UINT_MAX);
298 
299 
301  art::FindManyP<rb::Track> fmTrack3D(slicecol, evt, fTrack3DLabel);
302  std::vector<art::Ptr<rb::Prong>> prongs3D;
303 
304  if( fProngInput && fmProng3D.isValid())
305  prongs3D = fmProng3D.at(iClust);
306  if( !fProngInput && fmTrack3D.isValid() ){
307  std::vector<art::Ptr<rb::Track>> tracks3D;
308  tracks3D.clear();
309  tracks3D = fmTrack3D.at(iClust);
310  for (unsigned int ip=0; ip<tracks3D.size(); ++ip){
311  prongs3D.push_back(tracks3D[ip]);
312  }
313  }
314 
315  art::FindManyP<cvn::PixelMap> fmPixelMap3D(prongs3D, evt, fPixelMapInput);
316  for( unsigned int iProng = 0; iProng < prongs3D.size(); ++iProng ){
317  if(!fmPixelMap3D.isValid()) continue;
318 
319  const std::vector< art::Ptr<cvn::PixelMap> > pixelMaps3D = fmPixelMap3D.at(iProng);
320  if(pixelMaps3D.empty()) continue;
321 
322  std::vector<unsigned char> prongpm = (*pixelMaps3D[0]).PixelMapToVector(true);
323  std::vector<tensorflow::Tensor> tensorShE = vector_to_tensor_nue(prongpm, (*pixelMaps3D[0]).NPlanePerView(), (*pixelMaps3D[0]).NCell());
324 
325  std::vector<tensorflow::Tensor> resultShE = models.second->Predict({{"input_1",tensorShE[0]},{"input_2",tensorShE[1]}}, {"output_node0"});
326  auto tfoutputShE = resultShE[0].tensor<float,2>();
327 
328  resultShCol->emplace_back(tfoutputShE(0,0));
329  util::CreateAssn(*this, evt, *(resultShCol.get()),
330  prongs3D[iProng], *(assocp.get()), UINT_MAX);
331 
332  } // prong
333  } // nue slices
334  for(size_t iClust = 0; iClust < slicelist.size(); ++iClust) {
335  if(!fmPixelMap.isValid()) continue;
336  if(slicelist[iClust]->IsNoise()) continue;
337 
338  const std::vector<art::Ptr<cvn::PixelMap> > pixelMaps = fmPixelMap.at(iClust);
339  if(pixelMaps.empty()) continue;
340 
341  std::vector<unsigned char> evtpm = (*pixelMaps[0]).PixelMapToVector(true);
342  tensorflow::Tensor tensorEvtE = vector_to_tensor_numu(evtpm);
343 
344  std::vector<tensorflow::Tensor> resultHadE = hadModel->Predict({{"input",tensorEvtE}}, {"output_node0"});
345  auto tfoutputHadE = resultHadE[0].tensor<float,2>();
346  resultHadCol->emplace_back(tfoutputHadE(0,0));
347  util::CreateAssn(*this, evt, *(resultHadCol.get()),
348  slicelist[iClust], *(assocHad.get()), UINT_MAX);
349 
350  } // slices
351 
352  evt.put(std::move(resultEvtCol));
353  evt.put(std::move(resultShCol));
354  evt.put(std::move(assoc));
355  evt.put(std::move(assocp));
356  evt.put(std::move(resultHadCol));
357  evt.put(std::move(assocHad));
358  } // produce
359 
360 }
361 
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.
std::string fTFProtoBufHadE
bool PassCVNCut(const cvn::Result &result)
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Definition: EnvExpand.cxx:8
Results for Regression CVN.
std::string fTFProtoBufEvtE
void produce(art::Event &evt)
std::vector< tensorflow::Tensor > vector_to_tensor_nue(std::vector< unsigned char >, unsigned int ncells, unsigned int nplanes)
Definition: models.py:1
TString ip
Definition: loadincs.C:5
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
bool isRealData() const
Definition: Event.h:83
DEFINE_ART_MODULE(TestTMapFile)
PixelMap for CVN.
std::vector< Tensor > Predict(std::vector< std::pair< std::string, Tensor >> inputs, std::vector< std::string > outputLabels)
Definition: TFHandler.cxx:64
std::vector< float > fOutput
Vector of outputs from neural net.
Definition: Result.h:30
std::string fGeneratorLabel
std::pair< RegModel *, RegModel * > fTFFHC
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
const XML_Char * s
Definition: expat.h:262
std::string fCVNLabel1
std::pair< RegModel *, RegModel * > fTFRHC
Result for CVN.
const int cols[3]
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
static Var nueid(const std::shared_ptr< CAFAnaModel > &model)
Definition: SliceLIDVar.h:83
bool IsRHC(const art::Event &evt)
RegModel * GetModelsNumu(const art::Event &evt)
std::string fTFProtoBufRHCEvtE
string rel
Definition: shutoffs.py:11
std::string fCVNLabel3
size_type size() const
Definition: PtrVector.h:308
int nplanes
Definition: geom.C:145
std::string fSliceLabel
std::string fTFProtoBufRHCShE
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Result, basic output of CVN neural net.
Definition: Result.h:15
tensorflow::TFHandler RegModel
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
std::string fCVNLabel2
std::string fProng3DLabel
TRandom3 r(0)
std::string fProngModLabel
int ncells
Definition: geom.C:124
std::string fNuMILabel
RegCVNTF(fhicl::ParameterSet const &pset)
std::pair< RegModel *, RegModel * > GetModelsNue(const art::Event &evt)
std::string fTFProtoBufShE
Wrapper for Tensorflow which handles construction and prediction.
Definition: TFHandler.h:19
std::string fTFProtoBufRHCHadE
tensorflow::Tensor vector_to_tensor_numu(std::vector< unsigned char >)
bool failedToGet() const
Definition: Handle.h:196
std::string fPixelMapInput
std::string fTrack3DLabel