CVNEventDump_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file CVNEventDump_module.cc
3 // \brief Analyzer module for creating CVN PixelMap objects
4 // \author Alexander Radovic - a.radovic@gmail.com
5 // Leigh Whitehead - leigh.howard.whitehead@cern.ch
6 // - Added in truth based fiducial volume cuts
7 ////////////////////////////////////////////////////////////////////////
8 
9 // C/C++ includes
10 #include <iostream>
11 #include <sstream>
12 
13 // ROOT includes
14 #include "TTree.h"
15 #include "TH2F.h"
16 
17 // Framework includes
22 #include "art_root_io/TFileDirectory.h"
23 #include "art_root_io/TFileService.h"
25 #include "fhiclcpp/ParameterSet.h"
28 
29 // LArSoft includes
33 
40 
41 
42 
43 namespace cvn {
44  class CVNEventDump : public art::EDAnalyzer {
45  public:
46 
47  explicit CVNEventDump(fhicl::ParameterSet const& pset);
48  ~CVNEventDump();
49 
50  void analyze(const art::Event& evt) override;
51  void reconfigure(const fhicl::ParameterSet& pset);
52  void beginJob() override;
53  void endJob() override;
54 
55  private:
56 
58 
69  unsigned int fTopologyHits; // Number of hits for a track to be considered detectable
70  // for topology definitions.
71 
73  TTree* fTrainTree;
74 
75  /// Function to extract TH2 from PixelMap and write to TFile
76  void WriteMapTH2(const art::Event& evt, int slice, const PixelMap& pm);
77 
78  };
79 
80  //.......................................................................
82  : EDAnalyzer(pset),
83  fMVAAlg(pset.get<fhicl::ParameterSet>("MVAAlg"))
84  {
85  this->reconfigure(pset);
86  }
87 
88  //......................................................................
90  { }
91 
92  //......................................................................
94  {
95  fPixelMapInput = pset.get<std::string>("PixelMapInput");
96  fGenieGenModuleLabel = pset.get<std::string>("GenieGenModuleLabel");
97  fWriteMapTH2 = pset.get<bool> ("WriteMapTH2");
98  fApplyFidVol = pset.get<bool>("ApplyFidVol");
99  fEnergyNueLabel = pset.get<std::string> ("EnergyNueLabel");
100  fEnergyNumuLabel = pset.get<std::string> ("EnergyNumuLabel");
101  fEnergyNutauLabel = pset.get<std::string> ("EnergyNutauLabel");
102  fGetEnergyOutput = pset.get<bool> ("GetEnergyOutput");
103  fGetEventWeight = pset.get<bool> ("GetEventWeight");
104  fUseTopology = pset.get<bool>("UseTopology");
105  fTopologyHits = pset.get<unsigned int>("TopologyHitsCut");
106  }
107 
108  //......................................................................
110  {
111 
112 
114 
115  fTrainTree = tfs->make<TTree>("CVNTrainTree", "Training records");
116  fTrainTree->Branch("train", "cvn::TrainingData", &fTrain);
117 
118 
119  }
120 
121  //......................................................................
123  {
124 
125  }
126 
127  //......................................................................
129  {
130 
131  // Get the pixel maps
132  std::vector< art::Ptr< cvn::PixelMap > > pixelmaplist;
134  auto pixelmapListHandle = evt.getHandle< std::vector< cvn::PixelMap > >(itag1);
135  if (pixelmapListHandle)
136  art::fill_ptr_vector(pixelmaplist, pixelmapListHandle);
137 
138  // If we have no pixel map then return
139  if(pixelmaplist.size() == 0) return;
140 
142 
143  // * monte carlo
144  std::vector<art::Ptr<simb::MCTruth> > mclist;
145  auto mctruthListHandle = evt.getHandle< std::vector<simb::MCTruth> >(fGenieGenModuleLabel);
146  if (mctruthListHandle)
147  art::fill_ptr_vector(mclist, mctruthListHandle);
148 
149  //unsigned short nmclist= mclist.size();
150 
151  //std::cout<<"mctruth: "<<nmclist<<std::endl;
152 
153  art::Ptr<simb::MCTruth> truth = mclist[0];
154  simb::MCNeutrino truthN=truth->GetNeutrino();
155  //truth = mclist[0];
156 
158 
159  interaction = labels.GetInteractionType(truthN);
160  if(fUseTopology){
161  labels.GetTopology(truth,fTopologyHits);
162 // labels.PrintTopology();
163  }
164  float nuEnergy = 0;
165  float lepEnergy = 0;
166 // if(truth.NeutrinoSet()){
167  nuEnergy = truthN.Nu().E();
168  lepEnergy = truthN.Lepton().E();
169 // }
170 
171  // If outside the fiducial volume don't waste any time filling other variables
172  if(fApplyFidVol){
173  // Get the interaction vertex from the end point of the neutrino. This is
174  // because the start point of the lepton doesn't make sense for taus as they
175  // are decayed by the generator and not GEANT
176  TVector3 vtx = truthN.Nu().EndPosition().Vect();
177  bool isFid = (fabs(vtx.X())<310 && fabs(vtx.Y())<550 && vtx.Z()>50 && vtx.Z()<1244);
178  if(!isFid) return;
179  }
180 
181  float recoNueEnergy = 0.;
182  float recoNumuEnergy = 0.;
183  float recoNutauEnergy = 0.;
184  // Should we use the EnergyReco_module reconstructed energies?
185  if(fGetEnergyOutput){
186  // Get the nue info
187  if(fEnergyNueLabel != ""){
188  auto energyRecoNueHandle = evt.getHandle<dune::EnergyRecoOutput>(fEnergyNueLabel);
189  recoNueEnergy = energyRecoNueHandle->fNuLorentzVector.E();
190  }
191  // And the numu
192  if(fEnergyNumuLabel != ""){
193  auto energyRecoNumuHandle = evt.getHandle<dune::EnergyRecoOutput>(fEnergyNumuLabel);
194  recoNumuEnergy = energyRecoNumuHandle->fNuLorentzVector.E();
195  }
196  // And the nutau
197  if(fEnergyNutauLabel != ""){
198  auto energyRecoNutauHandle = evt.getHandle<dune::EnergyRecoOutput>(fEnergyNutauLabel);
199  recoNutauEnergy = energyRecoNutauHandle->fNuLorentzVector.E();
200  }
201  }
202 
203  // If we don't want to get the event weight then leave it as 1.0.
204  double eventWeight = 1.;
205  if(fGetEventWeight){
206  double mvaResult = 0.; // We don't care about this, but need it for the call
207  fMVAAlg.Run(evt,mvaResult,eventWeight);
208  }
209 
210  // Create the training data and add it to the tree
211  TrainingData train(interaction, nuEnergy, lepEnergy, recoNueEnergy,
212  recoNumuEnergy, recoNutauEnergy, eventWeight, *pixelmaplist[0]);
213  // Set the topology information
214  int topPDG = labels.GetPDG();
215  int nprot = labels.GetNProtons();
216  int npion = labels.GetNPions();
217  int npi0 = labels.GetNPizeros();
218  int nneut = labels.GetNNeutrons();
219  int toptype = labels.GetTopologyType();
220  int toptypealt = labels.GetTopologyTypeAlt();
221  if(fUseTopology){
222  train.SetTopologyInformation(topPDG, nprot, npion, npi0, nneut, toptype, toptypealt);
223  }
224  fTrain = &train;
225  fTrainTree->Fill();
226 
227  // Make a plot of the pixel map if required
228  if (fWriteMapTH2) WriteMapTH2(evt, 0, train.fPMap);
229 
230  }
231 
232  //----------------------------------------------------------------------
233 
234 
235 
236  void CVNEventDump::WriteMapTH2(const art::Event& evt, int slice, const PixelMap& pm)
237  {
238  std::stringstream name;
239  name << "PixelMap_r" << evt.run() << "_s" << evt.subRun()<< "_e" << evt.event() << "_sl" << slice;
240  std::stringstream nameL;
241  nameL << "PixelTruthMap_r" << evt.run() << "_s" << evt.subRun()<< "_e" << evt.event() << "_sl" << slice;
242  std::stringstream nameX;
243  nameX << "PixelMap_X_r" << evt.run() << "_s" << evt.subRun()<< "_e" << evt.event() << "_sl" << slice;
244  std::stringstream nameY;
245  nameY << "PixelMap_Y_r" << evt.run() << "_s" << evt.subRun()<< "_e" << evt.event() << "_sl" << slice;
246  std::stringstream nameZ;
247  nameZ << "PixelMap_Z_r" << evt.run() << "_s" << evt.subRun()<< "_e" << evt.event() << "_sl" << slice;
248  TH2F* hist = pm.ToTH2();
249  TH2F* histL = pm.ToLabTH2();
250  TH2F* histX = pm.SingleViewToTH2(0);
251  TH2F* histY = pm.SingleViewToTH2(1);
252  TH2F* histZ = pm.SingleViewToTH2(2);
253  hist->SetName(name.str().c_str());
254  histL->SetName(nameL.str().c_str());
255  histX->SetName(nameX.str().c_str());
256  histY->SetName(nameY.str().c_str());
257  histZ->SetName(nameZ.str().c_str());
258 
260 
261  TH2F* histWrite = tfs->make<TH2F>(*hist);
262  histWrite->Write();
263  TH2F* histWriteL = tfs->make<TH2F>(*histL);
264  histWriteL->GetZaxis()->SetRangeUser(0,10);
265  histWriteL->Write();
266  TH2F* histWriteX = tfs->make<TH2F>(*histX);
267  histWriteX->Write();
268  TH2F* histWriteY = tfs->make<TH2F>(*histY);
269  histWriteY->Write();
270  TH2F* histWriteZ = tfs->make<TH2F>(*histZ);
271  histWriteZ->Write();
272 
273  delete hist;
274  delete histWrite;
275  delete histL;
276  delete histWriteL;
277  delete histX;
278  delete histWriteX;
279  delete histY;
280  delete histWriteY;
281  delete histZ;
282  delete histWriteZ;
283 
284  }
285 
287 } // end namespace cvn
288 ////////////////////////////////////////////////////////////////////////
289 
290 
291 
292 
293 
294 
295 
static QCString name
Definition: declinfo.cpp:673
double E(const int i=0) const
Definition: MCParticle.h:233
EventNumber_t event() const
Definition: DataViewImpl.cc:85
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:77
const TLorentzVector & EndPosition() const
Definition: MCParticle.h:225
unsigned short GetNPizeros()
Definition: AssignLabels.h:36
unsigned int fTopologyHits
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
std::string string
Definition: nybbler.cc:12
TH2F * SingleViewToTH2(const unsigned int &view) const
Definition: PixelMap.cxx:206
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
void WriteMapTH2(const art::Event &evt, int slice, const PixelMap &pm)
Function to extract TH2 from PixelMap and write to TFile.
unsigned short GetNNeutrons()
Definition: AssignLabels.h:37
enum cvn::Interaction InteractionType
unsigned short GetTopologyType()
TH2F * ToLabTH2() const
Definition: PixelMap.cxx:186
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
PixelMap for CVN.
Utility class for truth labels.
Particle class.
InteractionType GetInteractionType(simb::MCNeutrino &truth)
unsigned short GetNProtons()
Definition: AssignLabels.h:34
void Run(const art::Event &evt, std::vector< double > &result, double &wgt)
TH2F * ToTH2() const
Return the pixel map as a 2D histogram for visualization.
Definition: PixelMap.cxx:166
unsigned short GetTopologyTypeAlt()
PixelMap fPMap
PixelMap for the event.
Definition: TrainingData.h:59
def train(model, train_files, valid_files, maskpatterns, epochs, batchsize, info)
Definition: train.py:16
void GetTopology(const art::Ptr< simb::MCTruth > truth, unsigned int nTopologyHits)
void beginJob() override
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
T get(std::string const &key) const
Definition: ParameterSet.h:271
std::string fEnergyNumuLabel
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:78
RunNumber_t run() const
Definition: DataViewImpl.cc:71
dunemva::MVAAlg fMVAAlg
CVNEventDump(fhicl::ParameterSet const &pset)
void analyze(const art::Event &evt) override
std::string fGenieGenModuleLabel
void reconfigure(const fhicl::ParameterSet &pset)
Something else. Tau? Hopefully we don&#39;t use this.
The TrainingData objects contains a PixelMap and the output class type, and any other bit that goes i...
Definition: TrainingData.h:20
PixelMap, basic input to CVN neural net.
Definition: PixelMap.h:22
std::string fEnergyNutauLabel
void endJob() override
TCEvent evt
Definition: DataStructs.cxx:7
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
unsigned short GetNPions()
Definition: AssignLabels.h:35
Event generator information.
Definition: MCNeutrino.h:18
void SetTopologyInformation(int pdg, int nproton, int npion, int npizero, int nneutron, int toptype, int toptypealt)