CVNProngEvaluatorTF_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file CVNProngEvaluatorTF_module.cc
3 // \brief Producer module creating 2D prong CVN neural net results
4 ////////////////////////////////////////////////////////////////////////
5 
6 // C/C++ includes
7 #include <iostream>
8 #include <sstream>
9 
10 // ROOT includes
11 #include "TTree.h"
12 #include "TH2F.h"
13 
14 // Framework includes
21 #include "fhiclcpp/ParameterSet.h"
26 
27 // NOvASoft includes
28 #include "RecoBase/PID.h"
29 #include "RecoBase/Cluster.h"
30 #include "RecoBase/FilterList.h"
31 #include "Utilities/AssociationUtil.h"
32 #include "SummaryData/SpillData.h"
33 
34 #include "CVN/func/PixelMap.h"
35 #include "CVN/func/Result.h"
37 #include "CVN/func/AssignLabels.h"
38 #include "CVN/func/TrainingData.h"
41 #include "CVN/func/ProngType.h"
43 
46 
47 using namespace cvn;
48 
49 namespace cvntf {
50 
52  public:
53  explicit CVNProngEvaluatorTF(fhicl::ParameterSet const& pset);
55 
56  bool filter(std::vector<unsigned char> pm, int view);
57  std::vector<tensorflow::Tensor> vector_to_tensor(std::vector<unsigned char>, std::vector<unsigned char>);
58 
59  void produce(art::Event& evt);
60  void beginJob();
61  void endJob();
62 
63 
64 
65  private:
66  /// Module lablel for input clusters
68 
69  /// Module label for input pixel maps
73 
74  /// Module lablel for input prongs
76  /// Instance lablel for input 3D prongs
77  // use 2d view
78  bool fIs2DProng;
82 
83  /// Check rb::IsFiltered?
86  bool fProngOnly;
88 
89  //
90 
91  /// Number of outputs fron neural net
92  unsigned int fNOutput;
93  bool fUseGeV;
94  std::vector<std::string> fPreselectionLabels;
98  };
99 
100 
101 
102  //.......................................................................
103  CVNProngEvaluatorTF::CVNProngEvaluatorTF(fhicl::ParameterSet const& pset):
104 
105  fClusterToken(consumes<std::vector<rb::Cluster>>(pset.get<std::string>("ClusterLabel"))),
106  fPixelMapInput (pset.get<std::string> ("PixelMapInput")),
107  fSlicePMLabel (pset.get<std::string> ("SlicePMLabel")),
108  fProngPMLabel (pset.get<std::string> ("ProngPMLabel")),
109  fProngModLabel (pset.get<std::string> ("ProngModLabel")),
110  fIs2DProng (pset.get<bool> ("Is2DProng")),
111  fProng3DLabel (pset.get<std::string> ("Prong3DLabel")),
112  fProng2DLabel (pset.get<std::string> ("Prong2DLabel")),
113  fProngLabel (pset.get<std::string> ("ProngLabel")),//(fIs2Dorong?"Prong2DLabel": "Prong3DLabel")),
114 
115  fObeyPreselection (pset.get<bool> ("ObeyPreselection")),
116  fSkipLongProngs(pset.get<bool> ("SkipLongProngs")),
117  fProngOnly (pset.get<bool> ("ProngOnly")),
118  fProngLengthCut (pset.get<double> ("ProngLengthCut")),
119  fNOutput(8),
120  fUseGeV(true),
121  fPreselectionLabels (pset.get<std::vector<std::string>> ("PreselectionLabels")),
122  fLibPath(pset.get<std::string>("LibPath")),
123  fModelName(pset.get<std::string>("ModelName")),
124  fTF(0)
125  {
128  produces< std::vector<cvn::Result> >();
129  produces< art::Assns<cvn::Result, rb::Cluster> >();
130  produces< art::Assns<cvn::Result, rb::Prong> >();
131  produces< std::vector<rb::PID> >();
132  produces< art::Assns<rb::PID, rb::Prong> >();
133  }
134 
135  //......................................................................
137  {
138  if(fTF) delete fTF;
139 
140  //======================================================================
141  // Clean up any memory allocated by your module
142  //======================================================================
143  }
144 
145  bool CVNProngEvaluatorTF::filter(std::vector<unsigned char> pm, int view)
146  {
147  float x(0), y(0), tmpx(0), tmpy(0);
148  for(int i = 0; i<8000; i++){
149  tmpx=pm[i];
150  if(tmpx>0){
151  x++;
152  }
153  }
154  for(int i = 8000; i<16000; i++){
155  tmpy = pm[i];
156  if(tmpy>0) y++;
157  }
158  if(view && y<=3) return false;
159  if((!view) && x<=3) return false;
160  else return true;
161  }
162 
163  std::vector<tensorflow::Tensor> CVNProngEvaluatorTF::vector_to_tensor(std::vector<unsigned char> prongpm, std::vector<unsigned char> slicepm)
164  {
165  tensorflow::Tensor tensorx(tensorflow::DT_FLOAT, {1,16000});
166  tensorflow::Tensor tensory(tensorflow::DT_FLOAT, {1,16000});
167 
168  auto relx = tensorx.tensor<float,2>();
169  auto rely = tensory.tensor<float,2>();
170  float x(0),y(0);
171 
172  for(int i = 0; i<8000; i++){
173  //odd-view
174  relx(0, i) = slicepm[i];
175  relx(0, i+8000) = prongpm[i];
176  x += prongpm[i];
177  }
178  for(int i = 8000; i< 16000; i++){
179  rely(0, i-8000) = slicepm[i];
180  rely(0, i) = prongpm[i];
181  y += prongpm[i];
182  }
183  return {tensorx, tensory};
184  }
185 
186  //......................................................................
188  { }
189 
190  //......................................................................
192  {
193  }
194 
195  //......................................................................
197  {
198 
199  /// Define containers for the things we're going to produce
200  std::unique_ptr< std::vector<Result> >
201  resultCol(new std::vector<Result>);
202  std::unique_ptr< std::vector<rb::PID> >
203  pidCol(new std::vector<rb::PID>);
204  std::unique_ptr< art::Assns<Result, rb::Cluster> >
206  std::unique_ptr< art::Assns<Result, rb::Prong> >
207  assocp(new art::Assns<Result, rb::Prong>);
208  std::unique_ptr< art::Assns<rb::PID, rb::Prong> >
209  assocpid(new art::Assns<rb::PID, rb::Prong>);
210 
211  // get the slices
213  evt.getByToken(fClusterToken, clusterHandle);
214  const std::vector<rb::Cluster> & clusters = *clusterHandle;
215 
216 // //get the spill data
217 // art::Handle<sumdata::SpillData> spillPot;
218 // if (!evt.isRealData()){
219 // evt.getByLabel(fGeneratorLabel, spillPot);
220 // }
221 // else{
222 // evt.getByLabel(fNuMIBeamLabel, spillPot);
223 // }
224 // // bool isRHC = false;
225 // if (spillPot.failedToGet()){
226 // mf::LogError("CVNEval") << "Spill Data not found, aborting without horn current information";
227 // abort();
228 // }
229 // else{
230 // // isRHC = spillPot->isRHC;
231 // }
232 
234  for(unsigned int i = 0; i<clusterHandle->size(); i++){
235  art::Ptr<rb::Cluster> clust(clusterHandle,i);
236  sliceCol.push_back(clust);
237  }
238 
239  art::FindManyP<PixelMap> fmSlicePM(clusterHandle, evt, fPixelMapInput);//art::InputTag(fPixelMapInput, fSlicePMLabel));
240  if(!fmSlicePM.isValid()) return;
241 
242  for(size_t iClust = 0; iClust < clusters.size(); iClust++) {
243  if(fObeyPreselection && rb::IsFiltered(evt,clusterHandle,iClust,fPreselectionLabels))
244  continue;
245  if(clusters[iClust].IsNoise()) continue;
246 
247  const std::vector<art::Ptr<PixelMap> > slicePMs = fmSlicePM.at(iClust);
248 
250  if(!fmProng.isValid()) continue;
251  //for 3d's x,y 2ds
252  std::vector<art::Ptr<rb::Prong>> prongs;
253 
254  prongs = fmProng.at(iClust);
255 
256  if(prongs.empty()) continue;
257  art::FindManyP<PixelMap> fmProngPM(prongs, evt, fPixelMapInput);
258  if( !fmProngPM.isValid()) continue;
259  // prong first
260  for( unsigned int iProng = 0; iProng < prongs.size(); iProng++ ){
261  art::Ptr<rb::Prong> prong = prongs[iProng];
262  int view = 2;
263 
264  if(fIs2DProng){
265  if(prong->NXCell()) view = 0;
266  else if(prong->NYCell()) view = 1;
267  else std::cout<<"How the hell do we get to no cells at all..."<<std::endl;
268  }
269  else{
270  if(!prong->NXCell() || !prong->NYCell()){
271  std::cout<<" 3d prong with an empty view!!!!!!"<<std::endl;
272  continue;
273  }
274  abort(); //abort here for now
275  }
276  if ( fSkipLongProngs && prong->TotalLength() > fProngLengthCut ) continue;
277  const std::vector<art::Ptr<PixelMap>> pm = fmProngPM.at(iProng);
278  if (pm.empty()) continue;
279  if (slicePMs.empty()) continue;
280  std::vector<unsigned char> prongpm = (*pm[0]).PixelMapToVector(fUseGeV);
281  std::vector<unsigned char> slicepm = (*slicePMs[0]).PixelMapToVector(fUseGeV);
282 
283  if(!filter(prongpm, view)) continue;
284  std::vector<tensorflow::Tensor> tensor = vector_to_tensor(prongpm, slicepm);
285  std::vector<tensorflow::Tensor> result = fTF->Predict({std::make_pair("input", tensor[view])}, {"output/Softmax"});
286 
287  float output[fNOutput];
288 
289  auto tfoutput = result[0].tensor<float,2>();
290 
291  std::vector<std::pair<int, double>> outputwithpdg;
292 
293  for ( unsigned int j = 0; j < fNOutput; j++ ){
294  outputwithpdg.push_back(std::make_pair(GetPDGByPType((cvn::PType)j), (double)tfoutput(0,j)));
295  pidCol->emplace_back(outputwithpdg[j].first, outputwithpdg[j].second);
296  util::CreateAssn(*this, evt, *(pidCol.get()),
297  prongs[iProng], *(assocpid.get()), UINT_MAX);
298  output[j] = (float)tfoutput(0,j);
299  }
300  const float* outputptr = output;
301  resultCol->emplace_back(outputptr, fNOutput);
302 // PixelMap tmp = *pm[0];
303 // (tmp).Print();
304 // for (unsigned int i = 0; i < fNOutput; i++) {std::cout<<*outputptr<<std::endl;outputptr++;}
305 
306  // Output object with PDGs rather than version dependednt cvn ID
307  util::CreateAssn(*this, evt, *(resultCol.get()),
308  sliceCol[iClust], *(assoc.get()), UINT_MAX);
309  util::CreateAssn(*this, evt, *(resultCol.get()),
310  prongs[iProng], *(assocp.get()), UINT_MAX);
311  }//iProng
312  }//iClust
313 
314  evt.put(std::move(pidCol));
315  evt.put(std::move(resultCol));
316  evt.put(std::move(assoc));
317  evt.put(std::move(assocp));
318  evt.put(std::move(assocpid));
319 
320  }
321 
322  //----------------------------------------------------------------------
323 
324 
325 }
327 ////////////////////////////////////////////////////////////////////////
ofstream output
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.
pdg code and pid value
bool filter(std::vector< unsigned char > pm, int view)
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Definition: EnvExpand.cxx:8
std::vector< std::string > fPreselectionLabels
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
DEFINE_ART_MODULE(TestTMapFile)
PixelMap for CVN.
Defines an enumeration for prong classification.
std::vector< Tensor > Predict(std::vector< std::pair< std::string, Tensor >> inputs, std::vector< std::string > outputLabels)
Definition: TFHandler.cxx:64
Module that kips a configurable number of events between each that it allows through. Note that this module really skips (N-1) events, it uses a simple modular division as its critera. This module will cut down the data sample to 1/N of its original size.
std::string fProngModLabel
Module lablel for input prongs.
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
const art::ProductToken< std::vector< rb::Cluster > > fClusterToken
Module lablel for input clusters.
virtual double TotalLength() const
Distance along prong to reach last cell hit.
Definition: Prong.cxx:186
Result for CVN.
void beginJob()
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
std::vector< tensorflow::Tensor > vector_to_tensor(std::vector< unsigned char >, std::vector< unsigned char >)
int evt
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
const double j
Definition: BetheBloch.cxx:29
Perform a "2 point" Hough transform on a collection of hits.
int GetPDGByPType(PType ptype)
Definition: ProngType.cxx:5
OStream cout
Definition: OStream.cxx:6
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
bool fObeyPreselection
Check rb::IsFiltered?
unsigned int fNOutput
Number of outputs fron neural net.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
bool fIs2DProng
Instance lablel for input 3D prongs.
std::string fPixelMapInput
Module label for input pixel maps.
PType
Definition: ProngType.h:18
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
bool getByToken(ProductToken< PROD > const &token, Handle< PROD > &result) const
Definition: DataViewImpl.h:387
Wrapper for Tensorflow which handles construction and prediction.
Definition: TFHandler.h:19