NCNNKeras_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: NCNNKeras
3 // Module Type: producer
4 // \file NCNNKeras_module.cc
5 // \brief A module to write out the NC Cosmic Rejection NN (Keras)
6 //
7 // \author Shaokai Yang - yangs9uc@gmail.com
8 ////////////////////////////////////////////////////////////////////////
9 // system includes
10 #include <memory>
11 #include <string>
12 #include <vector>
13 
14 // art includes
26 #include "cetlib_except/exception.h"
27 #include "cetlib/filesystem.h"
28 #include "cetlib/search_path.h"
29 #include "fhiclcpp/ParameterSet.h"
31 
32 // novasoft includes
33 #include "CVN/func/Result.h"
34 #include "Geometry/Geometry.h"
35 #include "MCCheater/BackTracker.h"
36 #include "NCID/func/KerasModel.h"
37 #include "NCID/NCKeras.h"
38 #include "RecoBase/Prong.h"
39 #include "RecoBase/Cluster.h"
40 #include "RecoBase/Shower.h"
41 #include "RecoBase/Track.h"
42 #include "RecoBase/Vertex.h"
44 //#include "ShowerLID/ParticleIDAlg.h"
45 #include "ShowerLID/ShowerLID.h"
47 #include "Utilities/AssociationUtil.h"
48 
49 // ROOT includes
50 #include "TFile.h"
51 #include "TH1D.h"
52 #include "TMath.h"
53 #include "TMatrixD.h"
54 #include "TSystem.h"
55 #include "TTree.h"
56 #include "TVector3.h"
57 #include "TVectorD.h"
58 
59 namespace ncid
60 {
61  class NCNNKeras : public art::EDProducer
62  {
63  public:
64 
65  explicit NCNNKeras(fhicl::ParameterSet const & p);
66  ~NCNNKeras();
67 
68  void beginJob() override;
69  void beginRun(art::Run &evt) override;
70  void reconfigure(const fhicl::ParameterSet& pset);
71  void produce(art::Event & e) override;
72 
73  private:
74 
86 
88  TTree* fTree;
89 
90  float cosmicid;
91  float shwGap;
92  float shwnhit;
93  float shwnhitx;
94  float shwnhity;
95  float shwwwidth;
96  float shwcalE;
97  float partptp;
98  float shwdirY;
99  float shwlen;
100  float shwxminusy;
101  float shwxplusy;
102  float shwxovery;
103  float nmiphit;
104  float nshwlid;
105  float kerasValue;
106 
107  /// Particle ID alorithm (loglikelihoods and dE/dx)
108  //slid::ParticleIDAlg* fParticleIDAlg;
109  /// FCL parameters for particle ID alorithm (loglikelihoods and dE/dx)
111  /// Return transverse momentum fraction of the event. Calculation based on reconstructed prongs.
112  double TransMomFraction(const std::vector< art::Ptr<rb::Prong> >& prongs) const;
113  };
114 
116  const art::Ptr<rb::Prong> b);
117 }
118 //////////////////////////////////////////////////////////////////////////
119 namespace ncid
120 {
122  fMCTruthModuleLabel(p.get<std::string>("MCTruthModuleLabel")),
123  //fParticleIDAlg(0),
124  fParticleIDAlgPSet(p.get<fhicl::ParameterSet>("ParticleIDAlgPSet"))
125  {
126  produces< std::vector<ncid::NCKeras> > ();
127  produces< art::Assns<ncid::NCKeras, rb::Cluster> > ();
128  this->reconfigure(p);
129  }
130 
131  //////////////////////////////////////////////////////////////////////////
133  {
134  }
135 
136  //////////////////////////////////////////////////////////////////////////
138  {
140  fTree = tfs->make<TTree>("Testtree", "Testtree");
141 
142  fTree->Branch("cosmicid", &cosmicid, "cosmicid/F");
143  fTree->Branch("shwGap", &shwGap, "shwGap/F");
144  fTree->Branch("shwnhit", &shwnhit, "shwnhit/F");
145  fTree->Branch("shwnhitx", &shwnhitx, "shwnhitx/F");
146  fTree->Branch("shwnhity", &shwnhity, "shwnhity/F");
147  fTree->Branch("shwwwidth", &shwwwidth, "shwwwidth/F");
148  fTree->Branch("shwcalE", &shwcalE, "shwcalE/F");
149  fTree->Branch("partptp", &partptp, "partptp/F");
150  fTree->Branch("shwdirY", &shwdirY, "shwdirY/F");
151  fTree->Branch("shwlen", &shwlen, "shwlen/F");
152  fTree->Branch("shwxminusy", &shwxminusy, "shwxminusy/F");
153  fTree->Branch("shwxplusy", &shwxplusy, "shwxplusy/F");
154  fTree->Branch("shwxovery", &shwxovery, "shwxovery/F");
155  fTree->Branch("nmiphit", &nmiphit, "nmiphit/F");
156  fTree->Branch("nshwlid", &nshwlid, "nshwlid/F");
157  fTree->Branch("kerasValue", &kerasValue, "kerasValue/F");
158  }
159 
160  //////////////////////////////////////////////////////////////////////////
162  {
164  // Make sure we can find the Neural Networks weight file
165  if(!cet::file_exists(pidlibpath + fNCKerasFile)){
166  throw cet::exception("NCKeras")
167  << "Couldn't find file " << pidlibpath + fNCKerasFile <<std::endl;
168  }
169  keras_5nn = new keras::KerasModel(pidlibpath + fNCKerasFile, false);
170  }
171 
172  //////////////////////////////////////////////////////////////////////////
174  {
175  fMCTruthModuleLabel = p.get< std::string >("MCTruthModuleLabel");
176  fSliceLabel = p.get< std::string >("SliceLabel");
177  fProngLabel = p.get< std::string >("ProngLabel");
178  fProngInstance = p.get< std::string >("ProngInstance");
179  fCVNLabel = p.get< std::string >("CVNLabel");
180  fShowerLabel = p.get< std::string >("ShowerLabel");
181  fShowerLIDLabel = p.get< std::string >("ShowerLIDLabel");
182  fPIDLibPath = p.get< std::string >("PIDLibPath");
183  fNCKerasFile = p.get< std::string >("NCKerasFile");
184  fParticleIDAlgPSet = p.get< fhicl::ParameterSet >("ParticleIDAlgPSet");
185  }
186 
187  //////////////////////////////////////////////////////////////////////////
189  {
190 
191  bool isMCTruthListEmpty = false;
193  try {
194  evt.getByLabel(fMCTruthModuleLabel, mclist);
195  if(mclist->empty()) { isMCTruthListEmpty = true; }
196  }
197  catch(...) { isMCTruthListEmpty = true; }
198  if(isMCTruthListEmpty){
199  mf::LogError("MCTruth Error") << "Error retrieving MCTruth list." << std::endl;
200  return;
201  }
202 
205 
206  std::unique_ptr< std::vector<ncid::NCKeras> >
207  ncCosKeras(new std::vector<ncid::NCKeras>);
208  std::unique_ptr< art::Assns<ncid::NCKeras, rb::Cluster> >
210 
211  // Get the slices in the event
213  evt.getByLabel(fSliceLabel, sliceHandle);
214  std::vector<art::Ptr<rb::Cluster> > slices;
215  art::fill_ptr_vector(slices, sliceHandle);
216 
217  // Get things associated with slices
218  art::FindManyP<rb::Prong> prongAssn(sliceHandle, evt, art::InputTag(fProngLabel, fProngInstance));
219  art::FindMany<slid::ShowerLID> showerLidAssn(sliceHandle, evt, fShowerLIDLabel);
220 
221  const int nslices = slices.size();
222  for(int isli = 0; isli < nslices; ++isli){
223 
224  std::vector<cheat::NeutrinoEffPur> SlicePurr = bt->SliceToMCTruth(slices[isli]->AllCells(), slices[isli]->AllCells());
225 
226  if(SlicePurr.size() == 0) continue;
227 
228  const art::Ptr<simb::MCTruth> &SliceTruth = SlicePurr[0].neutrinoInt;
229  if(!SliceTruth->NeutrinoSet()) continue;
230  const simb::MCNeutrino& nu = SliceTruth->GetNeutrino();
231  art::Ptr<rb::Cluster> thisSlice = slices.at(isli);
232 
233  if(thisSlice->IsNoise()) continue;
234  if(!prongAssn.isValid()) continue;
235 
236  int intCurrent = nu.CCNC();
237 
238  if(intCurrent != simb::kNC) continue; // select only CC
239 
240  ncid::NCKeras nckeras;
241  std::vector< art::Ptr<rb::Prong> > prongs = prongAssn.at(isli);
242  if(prongs.empty()) continue;
243 
244  nckeras.SetProngTransMom(TransMomFraction(prongs));
245 
246  art::FindManyP<cvn::Result> fmCVN(slices, evt, fCVNLabel);
247  if(!fmCVN.isValid()) continue; // CVN Validation
248  std::vector< art::Ptr<cvn::Result> > cvns = fmCVN.at(isli);
249  cosmicid = -5.; // CVN values
250  if(!cvns.empty()) cosmicid = cvns[0]->fOutput[14];
251  shwGap = -5.; // LIDShower values
252  shwnhit = -5.;
253  shwnhitx = -5.;
254  shwnhity = -5.;
255  shwwwidth = -5.;
256  shwcalE = -5.;
257  partptp = nckeras.ProngTransMom();
258  shwdirY = -5.;
259  shwlen = -5.;
260  kerasValue = -5.;
261 
262  if(!showerLidAssn.isValid()) continue; // showerLidAssn Validation
263 
264  std::vector<const slid::ShowerLID*> shwlids = showerLidAssn.at(isli);
265  if(shwlids.empty()) continue;
266 
267  sort(shwlids.begin(), shwlids.end(), slid::CompareByEnergy);
268  art::FindOneP<rb::Shower> foSh(shwlids, evt, fShowerLIDLabel);
269  nshwlid = shwlids.size();
270 
271  cet::maybe_ref< art::Ptr<rb::Shower> const > roSh(foSh.at(0));
272  art::Ptr<rb::Shower> shower = roSh.ref();
273  std::vector< art::Ptr<rb::Shower> > showers;
274  showers.push_back(shower);
275  art::FindOneP<slid::ShowerLID> foShLID(showers, evt,fShowerLIDLabel);
276  if(!foShLID.isValid()) continue; // foShLID validation
277 
278  cet::maybe_ref< art::Ptr<slid::ShowerLID> const > roShLID(foShLID.at(0));
279  art::Ptr<slid::ShowerLID> showerLID = roShLID.ref();
280  shwGap = showerLID->Gap();
281  shwnhit = shower->NCell();
282  shwwwidth = showerLID->Radius();
283  shwcalE = showerLID->ShowerEnergy();
284  shwnhitx = shower->NXCell();
285  shwdirY = (float)shower->Dir().Y();
286  shwnhity = shower->NYCell();
287  shwlen = shower->TotalLength();
291  nmiphit = thisSlice->NCell();
292 
293  // Keras model
294  keras::DataChunk *sample = new keras::DataChunkFlat();
295  float myfloats[] = {partptp,
296  cosmicid,
297  nmiphit,
298  shwcalE,
299  shwnhit,
300  nshwlid,
301  shwnhitx,
302  shwnhity,
303  shwxminusy,
304  shwxplusy,
305  shwxovery,
306  shwdirY,
307  shwlen,
308  shwwwidth,
309  shwGap};
310 
311  std::vector<float> v (myfloats, myfloats + sizeof(myfloats) / sizeof(float));
312  sample->set_data(v);
313  kerasValue = keras_5nn->compute_output(sample)[0];
314 
315  nckeras.SetCosPIDKeras(kerasValue);
316 
317  delete sample;
318  ncCosKeras->push_back(nckeras);
319 
320  util::CreateAssn(*this, evt, *ncCosKeras, thisSlice ,*sliceAssn);
321  fTree->Fill();
322  } // end loop over slices
323  evt.put(std::move(sliceAssn));
324  evt.put(std::move(ncCosKeras));
325  } // end of producer
326 
327  ///////////////////////////////////////////////////////////////////////
328  double NCNNKeras::TransMomFraction(const std::vector< art::Ptr<rb::Prong> >& prongs) const
329  {
330  if(prongs.empty()) return 0.;
331 
332  TVector3 ret;
333  for(const art::Ptr<rb::Prong>& prong: prongs){
334  if(!prong->Is3D()) continue;
335  const double w = prong->CalorimetricEnergy();
336  ret += w*prong->Dir();
337  }
338 
339  ret = ret.Unit();
340 
342  const TVector3 beamDir = geom->NuMIBeamDirection();
343  const double beamProj = ret.Dot(beamDir);
344 
345  return (ret-beamProj*beamDir).Mag();
346  }
347 
349  const art::Ptr<rb::Prong> b)
350  {
351  return a->CalorimetricEnergy() > b->CalorimetricEnergy();
352  }
353 
354 } // end namespace ncid
355 
357 ////////////////////////////////////////////////////////////////////////
Definition: FillPIDs.h:17
back track the reconstruction to the simulation
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.
fhicl::ParameterSet fParticleIDAlgPSet
Particle ID alorithm (loglikelihoods and dE/dx)
int CCNC() const
Definition: MCNeutrino.h:148
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
std::string fNCKerasFile
std::string fPIDLibPath
std::string fCosmicTrackLabel
void beginJob() override
std::string fMCTruthModuleLabel
std::string fProngLabel
const char * p
Definition: xmltok.h:285
std::string fVertexLabel
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Definition: EnvExpand.cxx:8
void SetCosPIDKeras(double input)
Definition: NCKeras.h:28
std::string fCVNLabel
virtual void set_data(std::vector< std::vector< std::vector< float > > > const &)
Definition: KerasModel.h:38
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void produce(art::Event &e) override
std::string fShowerLabel
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
double ProngTransMom() const
Transverse component of the energy-weighted average of prong directions.
Definition: NCKeras.h:20
DEFINE_ART_MODULE(TestTMapFile)
std::vector< float > compute_output(keras::DataChunk *dc)
Definition: KerasModel.cxx:354
bool CompareByEnergy(const slid::ShowerLID *a, const slid::ShowerLID *b)
Definition: ShowerLID.cxx:51
std::string fProngInstance
std::string fShowerLIDLabel
Definition: Run.h:31
double TransMomFraction(const std::vector< art::Ptr< rb::Prong > > &prongs) const
Return transverse momentum fraction of the event. Calculation based on reconstructed prongs...
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
Result for CVN.
double CalorimetricEnergy(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple estimate of neutrino energy.
Definition: Cluster.cxx:439
bool CompareByCalEnergy(const art::Ptr< rb::Prong > a, const art::Ptr< rb::Prong > b)
const double a
NCNNKeras(fhicl::ParameterSet const &p)
T get(std::string const &key) const
Definition: ParameterSet.h:231
Vertex location in position and time.
Definition: run.py:1
keras::KerasModel * keras_5nn
TVector3 NuMIBeamDirection() const
Direction of neutrinos from the NuMI beam (unit vector)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void reconfigure(const fhicl::ParameterSet &pset)
void geom(int which=0)
Definition: geom.C:163
const hit & b
Definition: hits.cxx:21
void SetProngTransMom(double input)
Definition: NCKeras.h:25
std::string fSliceLabel
bool file_exists(std::string const &qualified_filename)
bool NeutrinoSet() const
Definition: MCTruth.h:77
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
void beginRun(art::Run &evt) override
Float_t e
Definition: plot.C:35
Event generator information.
Definition: MCNeutrino.h:18
Float_t w
Definition: plot.C:20
std::vector< NeutrinoEffPur > SliceToMCTruth(const std::vector< const rb::CellHit * > &sliceHits, const std::vector< const rb::CellHit * > &allHits, bool sortPur=false) const
Given a collection of hits (often a slice), returns vector of structures of MCTruth, efficiency, and purity of that neutrino interaction relative to the slice. Efficiency is defined as FLS energy from neutrino interaction in slice / total FLS energy from neutrino interaction in event. This vector is sorted from the highest efficiency interaction to lowest. This function returns all MCTruth, including those without neutrino interactions, i.e. cosmics.
Definition: fwd.h:28
Encapsulate the geometry of one entire detector (near, far, ndos)