NCNNKerasVal_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: NCCosRej
3 // Module Type: analyzer
4 // \file NCNNKerasVal_module.cc
5 // \brief A module to for Validation
6 // \author Shaokai Yang - yangs9uc@gmail.com
7 ////////////////////////////////////////////////////////////////////////
8 // system includes
9 #include <memory>
10 #include <string>
11 #include <vector>
12 
13 // art includes
25 #include "cetlib_except/exception.h"
26 #include "cetlib/filesystem.h"
27 #include "cetlib/search_path.h"
28 #include "fhiclcpp/ParameterSet.h"
30 
31 // novasoft includes
32 #include "CosRej/NueCosRej.h"
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 "ShowerLID/ParticleIDAlg.h"
39 #include "ShowerLID/ShowerLID.h"
41 #include "RecoBase/Cluster.h"
42 #include "RecoBase/CellHit.h"
43 #include "RecoBase/Prong.h"
44 #include "RecoBase/RecoHit.h"
45 #include "RecoBase/Shower.h"
46 #include "RecoBase/Track.h"
47 #include "RecoBase/Vertex.h"
49 #include "Utilities/AssociationUtil.h"
50 
51 // ROOT includes
52 #include "TFile.h"
53 #include "TH1.h"
54 #include "TMath.h"
55 #include "TSystem.h"
56 #include "TTree.h"
57 #include "TVector3.h"
58 
59 
60 namespace ncid
61 {
63  {
64  public:
65  explicit NCNNKerasVal(fhicl::ParameterSet const& p);
66  virtual ~NCNNKerasVal();
67  void analyze(const art::Event& evt);
68  void beginJob();
69  void beginRun(art::Run const& run);
70  void reconfigure(fhicl::ParameterSet const& p);
71 
72  private:
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  float abdtValue;
107  float gbdtValue;
108  //slid::ParticleIDAlg* fParticleIDAlg;
110  };
111 
113  const art::Ptr<rb::Prong> b);
114 }
115 
116 //////////////////////////////////////////////////////////////////////////
117 namespace ncid
118 {
119  //////////////////////////////////////////////////////////////////////////
121  EDAnalyzer(p),
122  fMCTruthModuleLabel(p.get< std::string>("MCTruthModuleLabel")),
123  //fParticleIDAlg(0),
124  fParticleIDAlgPSet(p.get< fhicl::ParameterSet >("ParticleIDAlgPSet"))
125  {
126  this->reconfigure(p);
127  }
128 
129  //////////////////////////////////////////////////////////////////////////
131  {
132  }
133 
134  //////////////////////////////////////////////////////////////////////////
136  {
137  fMCTruthModuleLabel = p.get< std::string >("MCTruthModuleLabel");
138  fSliceLabel = p.get< std::string >("SliceLabel");
139  fCellHitLabel = p.get< std::string >("CellHitLabel");
140  fProngLabel = p.get< std::string >("ProngLabel");
141  fProngInstance = p.get< std::string >("ProngInstance");
142  fCVNLabel = p.get< std::string >("CVNLabel");
143  fShowerLabel = p.get< std::string >("ShowerLabel");
144  fShowerLIDLabel = p.get< std::string >("ShowerLIDLabel");
145  fPIDLibPath = p.get< std::string >("PIDLibPath");
146  fNCKerasFile = p.get< std::string >("NCKerasFile");
147  fCosRejLabel = p.get< std::string >("CosRejLabel");
148  fParticleIDAlgPSet = p.get< fhicl::ParameterSet >("ParticleIDAlgPSet");
149  }
150 
151  //////////////////////////////////////////////////////////////////////////
153  {
155  fTree = tfs->make<TTree>("Testtree", "Testtree");
156 
157  fTree->Branch("cosmicid", &cosmicid, "cosmicid/F");
158  fTree->Branch("shwGap", &shwGap, "shwGap/F");
159  fTree->Branch("shwnhit", &shwnhit, "shwnhit/F");
160  fTree->Branch("shwnhitx", &shwnhitx, "shwnhitx/F");
161  fTree->Branch("shwnhity", &shwnhity, "shwnhity/F");
162  fTree->Branch("shwwwidth", &shwwwidth, "shwwwidth/F");
163  fTree->Branch("shwcalE", &shwcalE, "shwcalE/F");
164  fTree->Branch("partptp", &partptp, "partptp/F");
165  fTree->Branch("shwdirY", &shwdirY, "shwdirY/F");
166  fTree->Branch("shwlen", &shwlen, "shwlen/F");
167  fTree->Branch("shwxminusy", &shwxminusy, "shwxminusy/F");
168  fTree->Branch("shwxplusy", &shwxplusy, "shwxplusy/F");
169  fTree->Branch("shwxovery", &shwxovery, "shwxovery/F");
170  fTree->Branch("nmiphit", &nmiphit, "nmiphit/F");
171  fTree->Branch("nshwlid", &nshwlid, "nshwlid/F");
172  fTree->Branch("kerasValue", &kerasValue, "kerasValue/F");
173  fTree->Branch("abdtValue", &abdtValue, "abdtValue/F");
174  fTree->Branch("gbdtValue", &gbdtValue, "gbdtValue/F");
175  }
176 
177  //////////////////////////////////////////////////////////////////////////
179  {
181  // Make sure we can find the Neural Networks weight file
182  if(!cet::file_exists(pidlibpath+fNCKerasFile)){
183  throw cet::exception("NCKeras")
184  << "Couldn't find file " << pidlibpath+fNCKerasFile <<std::endl;
185  }
186  keras_5nn = new keras::KerasModel(pidlibpath+fNCKerasFile, false);
187  }
188 
189  //////////////////////////////////////////////////////////////////////////
191  {
193  bool isMCTruthListEmpty = false;
195  try {
196  evt.getByLabel(fMCTruthModuleLabel, mclist);
197  if(mclist->empty()) isMCTruthListEmpty = true;
198  }
199  catch(...) { isMCTruthListEmpty = true; }
200  if(isMCTruthListEmpty){
201  mf::LogError("MCTruth Error") << "Error retrieving MCTruth list." << std::endl;
202  return;
203  }
204 
205  // Get the slices in the event
207  evt.getByLabel( fSliceLabel, sliceHandle);
208 
209  std::vector<art::Ptr<rb::Cluster> > slices;
210  art::fill_ptr_vector(slices, sliceHandle);
211 
213  evt.getByLabel(fCellHitLabel, htHandle);
214 
216  for(unsigned int m = 0; m < htHandle->size(); ++m){
217  allHits.push_back(art::Ptr<rb::CellHit>(htHandle,m));
218  }
219  // Get things associated with slices
220  art::FindManyP<rb::Prong> prongAssn(sliceHandle, evt, art::InputTag(fProngLabel, fProngInstance));
221  art::FindMany<slid::ShowerLID> showerLidAssn(sliceHandle, evt, fShowerLIDLabel);
223  art::FindManyP<rb::Vertex> vertexAssn(sliceHandle, evt, fVertexLabel);
224  const int nslices = slices.size();
225  cosmicid = -5.; // CVN values
226  shwGap = -5.; // LIDShower values
227  shwnhit = -5.;
228  shwnhitx = -5.;
229  shwnhity = -5.;
230  shwwwidth = -5.;
231  shwcalE = -5.;
232  partptp = -5.; //nckeras.ProngTransMom();
233  shwdirY = -5.;
234  shwlen = -5.;
235  kerasValue = -5.;
236 
237  // Loop over slices
238  for(int isli = 0; isli < nslices; ++isli){
239  art::Ptr<rb::Cluster> thisSlice = slices.at(isli);
240  if(thisSlice->IsNoise()) continue;
241  if(!prongAssn.isValid()) continue;
242  //Get True neutrino info
243  std::vector<cheat::NeutrinoEffPur> SlicePurr = bt->SliceToMCTruth(slices[isli]->AllCells(), slices[isli]->AllCells());
244 
245  if(SlicePurr.size() == 0) continue;
246 
247  const art::Ptr<simb::MCTruth> &SliceTruth = SlicePurr[0].neutrinoInt;
248  if(!SliceTruth->NeutrinoSet()) continue;
249 
250  ncid::NCKeras nckeras;
251  std::vector< art::Ptr<rb::Prong> > prongs = prongAssn.at(isli);
252  if(prongs.empty()) continue;
253 
254  art::FindManyP<cvn::Result> fmCVN(slices, evt, fCVNLabel);
255  // CVN Validation
256  if(!fmCVN.isValid()) continue;
257 
258  std::vector< art::Ptr<cosrej::NueCosRej> > cosrejs = fmcr.at(isli);
259  if(cosrejs.empty()) continue;
260 
261  std::vector< art::Ptr<cvn::Result> > cvns = fmCVN.at(isli);
262  // showerLidAssn Validation
263  if(!showerLidAssn.isValid()) continue;
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  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 
276  art::FindOneP<slid::ShowerLID> foShLID(showers, evt,fShowerLIDLabel);
277  // foShLID validation
278  if(!foShLID.isValid()) continue;
279 
280  cet::maybe_ref< art::Ptr<slid::ShowerLID> const > roShLID(foShLID.at(0));
281  art::Ptr<slid::ShowerLID> showerLID = roShLID.ref();
282  /////////////OK, Cuts///////////////
283  /// Event Quality Cuts
284  if(cosrejs[0]->HitsPerPlane() >= 8) continue;
285  if(showers.size() == 0) continue;
286  if(showerLID->Gap() >= 100.) continue;
287  if(thisSlice->MostContiguousPlanes(geo::kXorY) <= 2) continue;
288  /// NC Official Containment Cuts
289  if(cosrejs[0]->StartDistToFront() < 10) continue;
290  if(cosrejs[0]->StopDistToFront() < 10) continue;
291  if(cosrejs[0]->StartDistToBack() < 10) continue;
292  if(cosrejs[0]->StopDistToBack() < 10) continue;
293  if(cosrejs[0]->StartDistToTop() < 10) continue;
294  if(cosrejs[0]->StopDistToTop() < 10) continue;
295  if(cosrejs[0]->StartDistToWest() < 10) continue;
296  if(cosrejs[0]->StopDistToWest() < 10) continue;
297  if(cosrejs[0]->StartDistToEast() < 10) continue;
298  if(cosrejs[0]->StopDistToEast() < 10) continue;
299 
300  int nMipHitt = 0;
301  for(size_t hitIdx = 0; hitIdx < slices[isli]->NCell(); ++hitIdx){
302  const art::Ptr<rb::CellHit>& chptr = slices[isli]->Cell(hitIdx);
303  const rb::RecoHit Reco = slices[isli]->RecoHit(chptr);
304  if(!Reco.IsCalibrated()) continue;
305  if((Reco.PECorr() > 100.0) && (Reco.PECorr() < 245.0)) nMipHitt++;
306  } // end of hits loop
307 
308  if(!cvns.empty()) cosmicid = cvns[0]->fOutput[14];
309  shwGap = showerLID->Gap();
310  shwnhit = shower->NCell();
311  shwwwidth = showerLID->Radius();
312  shwcalE = showerLID->ShowerEnergy();
313  shwnhitx = shower->NXCell();
314  shwdirY = (float)shower->Dir().Y();
315  shwnhity = shower->NYCell();
316  shwlen = shower->TotalLength();
320  nmiphit = nMipHitt;
321  partptp = cosrejs[0]->ParticleShowerTransMom();
322 
323  // Keras model
324  keras::DataChunk *sample = new keras::DataChunkFlat();
325  float myfloats[] = {partptp,
326  cosmicid,
327  nmiphit,
328  shwcalE,
329  shwnhit,
330  nshwlid,
331  shwnhitx,
332  shwnhity,
333  shwxminusy,
334  shwxplusy,
335  shwxovery,
336  shwdirY,
337  shwlen,
338  shwwwidth,
339  shwGap};
340 
341  std::vector<float> v (myfloats, myfloats + sizeof(myfloats) / sizeof(float) );
342  sample->set_data(v);
343  kerasValue = keras_5nn->compute_output(sample)[0];
344  nckeras.SetCosPIDKeras(kerasValue);
345  delete sample;
346 
347  fTree->Fill();
348  } // end loop over slices
349  } // end of analyzer
350 } // end namespace ncid
351 
Definition: FillPIDs.h:17
back track the reconstruction to the simulation
NCNNKerasVal(fhicl::ParameterSet const &p)
X or Y views.
Definition: PlaneGeo.h:30
const char * p
Definition: xmltok.h:285
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Definition: EnvExpand.cxx:8
void SetCosPIDKeras(double input)
Definition: NCKeras.h:28
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
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
DEFINE_ART_MODULE(TestTMapFile)
fhicl::ParameterSet fParticleIDAlgPSet
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
Definition: Run.h:31
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
void beginRun(art::Run const &run)
unsigned short Cell() const
Definition: CellHit.h:40
Result for CVN.
keras::KerasModel * keras_5nn
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
bool CompareByCalEnergy(const art::Ptr< rb::Prong > a, const art::Ptr< rb::Prong > b)
const double a
T get(std::string const &key) const
Definition: ParameterSet.h:231
void reconfigure(fhicl::ParameterSet const &p)
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
Vertex location in position and time.
Definition: run.py:1
::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
int MostContiguousPlanes(geo::View_t view) const
Longest run of adjacent planes with hits.
Definition: Cluster.cxx:635
const hit & b
Definition: hits.cxx:21
std::string fMCTruthModuleLabel
bool file_exists(std::string const &qualified_filename)
bool NeutrinoSet() const
Definition: MCTruth.h:77
float PECorr() const
Definition: RecoHit.cxx:47
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
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)
void analyze(const art::Event &evt)