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
21 #include "art_root_io/TFileDirectory.h"
22 #include "art_root_io/TFileService.h"
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  EDProducer(p),
123  fMCTruthModuleLabel(p.get<std::string>("MCTruthModuleLabel")),
124  //fParticleIDAlg(0),
125  fParticleIDAlgPSet(p.get<fhicl::ParameterSet>("ParticleIDAlgPSet"))
126  {
127  produces< std::vector<ncid::NCKeras> > ();
128  produces< art::Assns<ncid::NCKeras, rb::Cluster> > ();
129  this->reconfigure(p);
130  }
131 
132  //////////////////////////////////////////////////////////////////////////
134  {
135  }
136 
137  //////////////////////////////////////////////////////////////////////////
139  {
141  fTree = tfs->make<TTree>("Testtree", "Testtree");
142 
143  fTree->Branch("cosmicid", &cosmicid, "cosmicid/F");
144  fTree->Branch("shwGap", &shwGap, "shwGap/F");
145  fTree->Branch("shwnhit", &shwnhit, "shwnhit/F");
146  fTree->Branch("shwnhitx", &shwnhitx, "shwnhitx/F");
147  fTree->Branch("shwnhity", &shwnhity, "shwnhity/F");
148  fTree->Branch("shwwwidth", &shwwwidth, "shwwwidth/F");
149  fTree->Branch("shwcalE", &shwcalE, "shwcalE/F");
150  fTree->Branch("partptp", &partptp, "partptp/F");
151  fTree->Branch("shwdirY", &shwdirY, "shwdirY/F");
152  fTree->Branch("shwlen", &shwlen, "shwlen/F");
153  fTree->Branch("shwxminusy", &shwxminusy, "shwxminusy/F");
154  fTree->Branch("shwxplusy", &shwxplusy, "shwxplusy/F");
155  fTree->Branch("shwxovery", &shwxovery, "shwxovery/F");
156  fTree->Branch("nmiphit", &nmiphit, "nmiphit/F");
157  fTree->Branch("nshwlid", &nshwlid, "nshwlid/F");
158  fTree->Branch("kerasValue", &kerasValue, "kerasValue/F");
159  }
160 
161  //////////////////////////////////////////////////////////////////////////
163  {
165  // Make sure we can find the Neural Networks weight file
166  if(!cet::file_exists(pidlibpath + fNCKerasFile)){
167  throw cet::exception("NCKeras")
168  << "Couldn't find file " << pidlibpath + fNCKerasFile <<std::endl;
169  }
170  keras_5nn = new keras::KerasModel(pidlibpath + fNCKerasFile, false);
171  }
172 
173  //////////////////////////////////////////////////////////////////////////
175  {
176  fMCTruthModuleLabel = p.get< std::string >("MCTruthModuleLabel");
177  fSliceLabel = p.get< std::string >("SliceLabel");
178  fProngLabel = p.get< std::string >("ProngLabel");
179  fProngInstance = p.get< std::string >("ProngInstance");
180  fCVNLabel = p.get< std::string >("CVNLabel");
181  fShowerLabel = p.get< std::string >("ShowerLabel");
182  fShowerLIDLabel = p.get< std::string >("ShowerLIDLabel");
183  fPIDLibPath = p.get< std::string >("PIDLibPath");
184  fNCKerasFile = p.get< std::string >("NCKerasFile");
185  fParticleIDAlgPSet = p.get< fhicl::ParameterSet >("ParticleIDAlgPSet");
186  }
187 
188  //////////////////////////////////////////////////////////////////////////
190  {
191 
192  bool isMCTruthListEmpty = false;
194  try {
195  evt.getByLabel(fMCTruthModuleLabel, mclist);
196  if(mclist->empty()) { isMCTruthListEmpty = true; }
197  }
198  catch(...) { isMCTruthListEmpty = true; }
199  if(isMCTruthListEmpty){
200  mf::LogError("MCTruth Error") << "Error retrieving MCTruth list." << std::endl;
201  return;
202  }
203 
206 
207  std::unique_ptr< std::vector<ncid::NCKeras> >
208  ncCosKeras(new std::vector<ncid::NCKeras>);
209  std::unique_ptr< art::Assns<ncid::NCKeras, rb::Cluster> >
211 
212  // Get the slices in the event
214  evt.getByLabel(fSliceLabel, sliceHandle);
215  std::vector<art::Ptr<rb::Cluster> > slices;
216  art::fill_ptr_vector(slices, sliceHandle);
217 
218  // Get things associated with slices
219  art::FindManyP<rb::Prong> prongAssn(sliceHandle, evt, art::InputTag(fProngLabel, fProngInstance));
220  art::FindMany<slid::ShowerLID> showerLidAssn(sliceHandle, evt, fShowerLIDLabel);
221 
222  const int nslices = slices.size();
223  for(int isli = 0; isli < nslices; ++isli){
224 
225  std::vector<cheat::NeutrinoEffPur> SlicePurr = bt->SliceToMCTruth(slices[isli]->AllCells(), slices[isli]->AllCells());
226 
227  if(SlicePurr.size() == 0) continue;
228 
229  const art::Ptr<simb::MCTruth> &SliceTruth = SlicePurr[0].neutrinoInt;
230  if(!SliceTruth->NeutrinoSet()) continue;
231  const simb::MCNeutrino& nu = SliceTruth->GetNeutrino();
232  art::Ptr<rb::Cluster> thisSlice = slices.at(isli);
233 
234  if(thisSlice->IsNoise()) continue;
235  if(!prongAssn.isValid()) continue;
236 
237  int intCurrent = nu.CCNC();
238 
239  if(intCurrent != simb::kNC) continue; // select only CC
240 
241  ncid::NCKeras nckeras;
242  std::vector< art::Ptr<rb::Prong> > prongs = prongAssn.at(isli);
243  if(prongs.empty()) continue;
244 
245  nckeras.SetProngTransMom(TransMomFraction(prongs));
246 
247  art::FindManyP<cvn::Result> fmCVN(slices, evt, fCVNLabel);
248  if(!fmCVN.isValid()) continue; // CVN Validation
249  std::vector< art::Ptr<cvn::Result> > cvns = fmCVN.at(isli);
250  cosmicid = -5.; // CVN values
251  if(!cvns.empty()) cosmicid = cvns[0]->fOutput[14];
252  shwGap = -5.; // LIDShower values
253  shwnhit = -5.;
254  shwnhitx = -5.;
255  shwnhity = -5.;
256  shwwwidth = -5.;
257  shwcalE = -5.;
258  partptp = nckeras.ProngTransMom();
259  shwdirY = -5.;
260  shwlen = -5.;
261  kerasValue = -5.;
262 
263  if(!showerLidAssn.isValid()) continue; // showerLidAssn Validation
264 
265  std::vector<const slid::ShowerLID*> shwlids = showerLidAssn.at(isli);
266  if(shwlids.empty()) continue;
267 
268  sort(shwlids.begin(), shwlids.end(), slid::CompareByEnergy);
269  art::FindOneP<rb::Shower> foSh(shwlids, evt, fShowerLIDLabel);
270  nshwlid = shwlids.size();
271 
272  cet::maybe_ref< art::Ptr<rb::Shower> const > roSh(foSh.at(0));
273  art::Ptr<rb::Shower> shower = roSh.ref();
274  std::vector< art::Ptr<rb::Shower> > showers;
275  showers.push_back(shower);
276  art::FindOneP<slid::ShowerLID> foShLID(showers, evt,fShowerLIDLabel);
277  if(!foShLID.isValid()) continue; // foShLID validation
278 
279  cet::maybe_ref< art::Ptr<slid::ShowerLID> const > roShLID(foShLID.at(0));
280  art::Ptr<slid::ShowerLID> showerLID = roShLID.ref();
281  shwGap = showerLID->Gap();
282  shwnhit = shower->NCell();
283  shwwwidth = showerLID->Radius();
284  shwcalE = showerLID->ShowerEnergy();
285  shwnhitx = shower->NXCell();
286  shwdirY = (float)shower->Dir().Y();
287  shwnhity = shower->NYCell();
288  shwlen = shower->TotalLength();
292  nmiphit = thisSlice->NCell();
293 
294  // Keras model
295  keras::DataChunk *sample = new keras::DataChunkFlat();
296  float myfloats[] = {partptp,
297  cosmicid,
298  nmiphit,
299  shwcalE,
300  shwnhit,
301  nshwlid,
302  shwnhitx,
303  shwnhity,
304  shwxminusy,
305  shwxplusy,
306  shwxovery,
307  shwdirY,
308  shwlen,
309  shwwwidth,
310  shwGap};
311 
312  std::vector<float> v (myfloats, myfloats + sizeof(myfloats) / sizeof(float));
313  sample->set_data(v);
314  kerasValue = keras_5nn->compute_output(sample)[0];
315 
316  nckeras.SetCosPIDKeras(kerasValue);
317 
318  delete sample;
319  ncCosKeras->push_back(nckeras);
320 
321  util::CreateAssn(evt, *ncCosKeras, thisSlice ,*sliceAssn);
322  fTree->Fill();
323  } // end loop over slices
324  evt.put(std::move(sliceAssn));
325  evt.put(std::move(ncCosKeras));
326  } // end of producer
327 
328  ///////////////////////////////////////////////////////////////////////
329  double NCNNKeras::TransMomFraction(const std::vector< art::Ptr<rb::Prong> >& prongs) const
330  {
331  if(prongs.empty()) return 0.;
332 
333  TVector3 ret;
334  for(const art::Ptr<rb::Prong>& prong: prongs){
335  if(!prong->Is3D()) continue;
336  const double w = prong->CalorimetricEnergy();
337  ret += w*prong->Dir();
338  }
339 
340  ret = ret.Unit();
341 
343  const TVector3 beamDir = geom->NuMIBeamDirection();
344  const double beamProj = ret.Dot(beamDir);
345 
346  return (ret-beamProj*beamDir).Mag();
347  }
348 
350  const art::Ptr<rb::Prong> b)
351  {
352  return a->CalorimetricEnergy() > b->CalorimetricEnergy();
353  }
354 
355 } // end namespace ncid
356 
358 ////////////////////////////////////////////////////////////////////////
Definition: FillPIDs.h:18
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:77
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
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
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:21
double TransMomFraction(const std::vector< art::Ptr< rb::Prong > > &prongs) const
Return transverse momentum fraction of the event. Calculation based on reconstructed prongs...
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
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
int evt
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)
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:78
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:291
void beginRun(art::Run &evt) override
Float_t e
Definition: plot.C:35
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
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:29
Encapsulate the geometry of one entire detector (near, far, ndos)
enum BeamMode string