RegCNNEventDump_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file RegCNNEventDump_module.cc
3 // \brief Analyzer module for creating RegCNN PixelMap objects modified from CVNEventDump_module.cc
4 // \author Ilsoo Seong - iseong@uci.edu
5 //
6 // Modifications to implement 3D maps
7 // - Wenjie Wu - wenjieww@uci.edu
8 ////////////////////////////////////////////////////////////////////////
9 
10 // C/C++ includes
11 #include <iostream>
12 #include <sstream>
13 
14 // ROOT includes
15 #include "TTree.h"
16 #include "TH2F.h"
17 
18 // Framework includes
23 #include "art_root_io/TFileDirectory.h"
24 #include "art_root_io/TFileService.h"
26 #include "fhiclcpp/ParameterSet.h"
29 //#include "art/Framework/Core/FindManyP.h"
30 
31 // LArSoft includes
37 
43 
46 
47 namespace cnn {
49  public:
50 
51  explicit RegCNNEventDump(fhicl::ParameterSet const& pset);
53 
54  void analyze(const art::Event& evt) override;
55  void reconfigure(const fhicl::ParameterSet& pset);
56  void beginJob() override;
57  void endJob() override;
58 
59  private:
60 
68 
69  //art::ServiceHandle<cheat::BackTracker> fBT;
70 
71  /// Function to extract TH2 from PixelMap and write to TFile
72  void WriteMapTH2(const art::Event& evt, int slice, const RegPixelMap& pm);
73  void WriteMapTH3(const art::Event& evt, int slice, const RegPixelMap3D& pm);
74  void WriteCroppedMapTH3(const art::Event& evt, int slice, const RegPixelMap3D& pm);
75 
76  };
77 
78  //.......................................................................
80  : EDAnalyzer(pset)
81  {
82  this->reconfigure(pset);
83  }
84 
85  //......................................................................
87  { }
88 
89  //......................................................................
91  {
92  fPixelMapInput = pset.get<std::string>("PixelMapInput");
93  fPixelMap3DInput = pset.get<std::string>("PixelMap3DInput");
94  fGenieGenModuleLabel = pset.get<std::string>("GenieGenModuleLabel");
95  fWriteMapTH2 = pset.get<bool> ("WriteMapTH2");
96  fWriteMapTH3 = pset.get<bool> ("WriteMapTH3");
97  fWriteCroppedMapTH3 = pset.get<bool> ("WriteCroppedMapTH3");
98  fApplyFidVol = pset.get<bool> ("ApplyFidVol");
99 
100  }
101 
102  //......................................................................
104  {
105 
106  }
107 
108  //......................................................................
110  {
111 
112  }
113 
114  //......................................................................
116  {
117 
118  // Get the pixel maps
119  std::vector< art::Ptr< cnn::RegPixelMap > > pixelmaplist;
121  auto pixelmapListHandle = evt.getHandle< std::vector< cnn::RegPixelMap > >(itag1);
122  if (pixelmapListHandle)
123  art::fill_ptr_vector(pixelmaplist, pixelmapListHandle);
124 
125  unsigned short nhits = pixelmaplist.size();
126 
127  // Get the 3D pixel maps
128  std::vector<art::Ptr<cnn::RegPixelMap3D> > pm3Dlist;
130  auto pm3DListHandle = evt.getHandle<std::vector<cnn::RegPixelMap3D> >(itag2);
131  if (pm3DListHandle)
132  art::fill_ptr_vector(pm3Dlist, pm3DListHandle);
133 
134  unsigned short npm3D = pm3Dlist.size();
135  std::cout<<"npm3D: "<<npm3D<<std::endl;
136 
137  // If we have no hits then return
138  if (npm3D < 1 && nhits < 1) return;
139 
140  // * monte carlo
141  std::vector<art::Ptr<simb::MCTruth> > mclist;
142  auto mctruthListHandle = evt.getHandle< std::vector<simb::MCTruth> >(fGenieGenModuleLabel);
143  if (mctruthListHandle)
144  art::fill_ptr_vector(mclist, mctruthListHandle);
145 
146  //unsigned short nmclist= mclist.size();
147 
148  //std::cout<<"mctruth: "<<nmclist<<std::endl;
149 
150  art::Ptr<simb::MCTruth> truth = mclist[0];
151  simb::MCNeutrino truthN=truth->GetNeutrino();
152  std::cout<<"nuPDG : "<<truthN.Nu().PdgCode()<<", ";
153  std::cout<<"lepPDG : "<<truthN.Lepton().PdgCode()<<std::endl;
154  //truth = mclist[0];
155 
156  //float nuEnergy = 0;
157  //float lepEnergy = 0;
158  ////if(truth.NeutrinoSet()){
159  //nuEnergy = truthN.Nu().E();
160  //lepEnergy = truthN.Lepton().E();
161  ////}
162 
163  // If outside the fiducial volume don't waste any time filling other variables
164  if(fApplyFidVol){
165  // Get the interaction vertex from the end point of the neutrino. This is
166  // because the start point of the lepton doesn't make sense for taus as they
167  // are decayed by the generator and not GEANT
168  TVector3 vtx = truthN.Nu().EndPosition().Vect();
169  bool isFid = (fabs(vtx.X())<310 && fabs(vtx.Y())<550 && vtx.Z()>50 && vtx.Z()<1244);
170  if(!isFid) return;
171  }
172 
173  // Make a plot of the pixel map if required
174  if (fWriteMapTH2) WriteMapTH2(evt, 0, *pixelmaplist[0]);
175  if (fWriteMapTH3) WriteMapTH3(evt, 0, *pm3Dlist[0]);
176  if (fWriteCroppedMapTH3) WriteCroppedMapTH3(evt, 0, *pm3Dlist[0]);
177 
178  }
179 
180  //----------------------------------------------------------------------
181 
182 
183 
184  void RegCNNEventDump::WriteMapTH2(const art::Event& evt, int slice, const RegPixelMap& pm)
185  {
186  std::stringstream name;
187  name << "RegPixelMap_r" << evt.run() << "_s" << evt.subRun()<< "_e" << evt.event() << "_sl" << slice;
188  std::stringstream nameL;
189  nameL << "PixelTruthMap_r" << evt.run() << "_s" << evt.subRun()<< "_e" << evt.event() << "_sl" << slice;
190  std::stringstream nameX;
191  nameX << "RegPixelMap_X_r" << evt.run() << "_s" << evt.subRun()<< "_e" << evt.event() << "_sl" << slice;
192  std::stringstream nameY;
193  nameY << "RegPixelMap_Y_r" << evt.run() << "_s" << evt.subRun()<< "_e" << evt.event() << "_sl" << slice;
194  std::stringstream nameZ;
195  nameZ << "RegPixelMap_Z_r" << evt.run() << "_s" << evt.subRun()<< "_e" << evt.event() << "_sl" << slice;
196  TH2F* hist = pm.ToTH2();
197  TH2F* histL = pm.ToLabTH2();
198  TH2F* histX = pm.SingleViewToTH2(0);
199  TH2F* histY = pm.SingleViewToTH2(1);
200  TH2F* histZ = pm.SingleViewToTH2(2);
201  hist->SetName(name.str().c_str());
202  histL->SetName(nameL.str().c_str());
203  histX->SetName(nameX.str().c_str());
204  histY->SetName(nameY.str().c_str());
205  histZ->SetName(nameZ.str().c_str());
206 
208 
209  TH2F* histWrite = tfs->make<TH2F>(*hist);
210  histWrite->Write();
211  TH2F* histWriteL = tfs->make<TH2F>(*histL);
212  histWriteL->GetZaxis()->SetRangeUser(0,10);
213  histWriteL->Write();
214  TH2F* histWriteX = tfs->make<TH2F>(*histX);
215  histWriteX->Write();
216  TH2F* histWriteY = tfs->make<TH2F>(*histY);
217  histWriteY->Write();
218  TH2F* histWriteZ = tfs->make<TH2F>(*histZ);
219  histWriteZ->Write();
220 
221  delete hist;
222  delete histWrite;
223  delete histL;
224  delete histWriteL;
225  delete histX;
226  delete histWriteX;
227  delete histY;
228  delete histWriteY;
229  delete histZ;
230  delete histWriteZ;
231 
232  }
233 
234  void RegCNNEventDump::WriteMapTH3(const art::Event& evt, int slice, const RegPixelMap3D& pm) {
235  std::cout<<"Write Map to TH3"<<std::endl;
236  std::stringstream name;
237  name << "RegPixelMap3D_r" << evt.run() << "_s" << evt.subRun() << "_e" << evt.event() << "_sl" << slice;
238  TH3F* hist = pm.ToTH3();
239  hist->SetName(name.str().c_str());
240 
242  TH3F* histWrite = tfs->make<TH3F>(*hist);
243  histWrite->Write();
244 
245  delete hist;
246  delete histWrite;
247  }
248 
250  std::cout<<"Write Cropped Map to TH3"<<std::endl;
251  std::stringstream name;
252  name << "RegCroppedPixelMap3D_r" << evt.run() << "_s" << evt.subRun() << "_e" << evt.event() << "_sl" << slice;
253  TH3F* hist = pm.ToCroppedTH3();
254  hist->SetName(name.str().c_str());
255 
257  TH3F* histWrite = tfs->make<TH3F>(*hist);
258  histWrite->Write();
259 
260  delete hist;
261  delete histWrite;
262  }
263 
265 } // end namespace cnn
266 ////////////////////////////////////////////////////////////////////////
267 
268 
static QCString name
Definition: declinfo.cpp:673
TH3F * ToCroppedTH3() const
int PdgCode() const
Definition: MCParticle.h:212
EventNumber_t event() const
Definition: DataViewImpl.cc:85
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:77
TH3F * ToTH3() const
void WriteMapTH3(const art::Event &evt, int slice, const RegPixelMap3D &pm)
const TLorentzVector & EndPosition() const
Definition: MCParticle.h:225
void WriteCroppedMapTH3(const art::Event &evt, int slice, const RegPixelMap3D &pm)
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
std::string string
Definition: nybbler.cc:12
RegPixelMap, basic input to CNN neural net.
Definition: RegPixelMap.h:22
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
void reconfigure(const fhicl::ParameterSet &pset)
Particle class.
art framework interface to geometry description
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
TH2F * ToLabTH2() const
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
TH2F * SingleViewToTH2(const unsigned int &view) const
T get(std::string const &key) const
Definition: ParameterSet.h:271
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:78
RegPixelMap3D for RegCNN modified from PixelMap.h.
RunNumber_t run() const
Definition: DataViewImpl.cc:71
void analyze(const art::Event &evt) override
Defines an enumeration for cellhit classification.
Declaration of signal hit object.
Access the description of detector geometry.
RegPixelMap3D, input to 3D CNN neural net.
Definition: RegPixelMap3D.h:22
TCEvent evt
Definition: DataStructs.cxx:7
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
RegPixelMap for RegCNN modified from PixelMap.h.
TH2F * ToTH2() const
Return the pixel map as a 2D histogram for visualization.
Event generator information.
Definition: MCNeutrino.h:18
RegCNNEventDump(fhicl::ParameterSet const &pset)
void WriteMapTH2(const art::Event &evt, int slice, const RegPixelMap &pm)
Function to extract TH2 from PixelMap and write to TFile.
QTextStream & endl(QTextStream &s)