CVNEventProngDump_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////
2 // \file CVNClassify_module.cc
3 // \brief Analyzer for creating both slice and prong by prong CVN
4 // PixelMap objects
5 // \author Fernanda Psihas - psihas@fnal.gov
6 ////////////////////////////////////////////////////////////////////////
7 
8 // C/C++ includes
9 #include <iostream>
10 #include <sstream>
11 
12 // ROOT includes
13 #include "TTree.h"
14 #include "TH2F.h"
15 #include "TVector3.h"
16 
17 // Framework includes
24 #include "fhiclcpp/ParameterSet.h"
28 
29 
30 // NOvASoft includes
31 #include "Geometry/Geometry.h"
32 #include "RecoBase/Cluster.h"
33 #include "RecoBase/FilterList.h"
34 #include "RecoBase/Prong.h"
35 #include "RecoBase/Shower.h"
36 #include "RecoBase/Vertex.h"
37 #include "MCCheater/BackTracker.h"
38 #include "ShowerLID/ShowerLID.h"
42 #include "Utilities/AssociationUtil.h"
43 
44 #include "CVN/func/AssignLabels.h"
45 #include "CVN/func/TrainingData.h"
50 #include "CVN/func/ProngType.h"
52 
55 
56 
57 namespace cvn {
59  public:
60  explicit CVNEventProngDump(fhicl::ParameterSet const& pset);
62 
63  void analyze(const art::Event& evt);
64  void reconfigure(const fhicl::ParameterSet& pset);
65  void beginJob();
66  void endJob();
67 
68 
69 
70  private:
71 
73 
74  /// Module label for input prongs
76 
77  /// Instance label for input 3D prongs
79 
80  /// Instance label for input 2D prongs
82 
83  /// Instance label for cluster pixelmaps
85 
86  /// Instance label for prong pixelmaps
88 
89  /// Instance label for vertex
91 
92  /// Instance label for showers
94 
95  /// Instance label for showers
97 
98  /// Option to save maps for prongs
99  bool fMapProngs;
100 
103 
105  std::vector<std::string> fPreselectionLabels;
106 
108 
109  TTree* fTrainTree;
110 
112 
113  /// Function to extract TH2 from PixelMap and write to TFile
114  void WriteMapTH2(const art::Event& evt, const PixelMap& pm,
115  bool isProng, unsigned int slice, unsigned int id);
116 
117  std::string fCAFLabel; ///< label for process that made standard records
118  bool fApplyCAFCuts; ///< do you want to apply CAF-level cuts?
119  int fCutType; ///< what cuts do you want to apply? (see CAFCutter.h)
120 
121  recovalid::CAFCutter fCAFCutter; ///< CAFCutter helper class
122  };
123 
124 
125 
126  //.......................................................................
128  : EDAnalyzer(pset)
129  {
130  this->reconfigure(pset);
131  }
132 
133  //......................................................................
135  {
136  }
137 
138  //......................................................................
140  {
141  fClusterLabel = pset.get<std::string>("ClusterLabel");
142  fPixelMapInput = pset.get<std::string>("PixelMapInput");
143  fProngModLabel = pset.get<std::string>("ProngModLabel");
144  fProng3DLabel = pset.get<std::string>("Prong3DLabel");
145  fProng2DLabel = pset.get<std::string>("Prong2DLabel");
146  fClusterPMLabel = pset.get<std::string>("ClusterPMLabel");
147  fProngPMLabel = pset.get<std::string>("ProngPMLabel");
148  fVertexLabel = pset.get<std::string>("VertexLabel");
149  fShowerLabel = pset.get<std::string>("ShowerLabel");
150  fShowerLIDLabel = pset.get<std::string>("ShowerLIDLabel");
151  fMapProngs = pset.get<bool> ("MapProngs");
152  fWriteMapTH2 = pset.get<bool> ("WriteMapTH2");
153  fObeyPreselection = pset.get<bool> ("ObeyPreselection");
154  fPreselectionLabels = pset.get< std::vector<std::string> >("PreselectionLabels");
155  fCAFLabel = pset.get<std::string>("CAFLabel");
156  fApplyCAFCuts = pset.get<bool> ("ApplyCAFCuts");
157  fCutType = pset.get<int> ("CutType");
158 
159  }
160 
161  //......................................................................
163  {
164 
166 
168  fTrainTree = tfs->make<TTree>("CVNTrainTree", "Training records");
169  fTrainTree->Branch("EventPMapsData", "cvn::EventLabeledPMaps", &fEventPMapsData);
170 
171  }
172 
173  //......................................................................
175  {
176  }
177 
178  //......................................................................
180  {
181 
183 
184  InteractionType interaction = kOther;
185  FinalStateType finalstate = kOtherFS;
186  FinalStateProngsType finalstateprongs = kOtherFSp;
187 
188  // get the slices
190  evt.getByLabel(fClusterLabel, clusterHandle);
191  const std::vector<rb::Cluster> & clusters = *clusterHandle;
193  for(unsigned int i=0; i<clusterHandle->size(); ++i){
194  slices.push_back(art::Ptr<rb::Cluster>(clusterHandle, i));
195  }
196 
197  // Get the Standard Records
199 
200  art::FindManyP<PixelMap> fmPixelMap(clusterHandle, evt, fPixelMapInput);
201 
202  std::vector<int> bestNuId;
203  std::vector<cheat::NeutrinoWithIndex> nus;
204  std::vector<std::vector<cheat::NeutrinoEffPur>> sEffPur;
205  nus = fBT->allMCTruth();
206  sEffPur = fBT->SlicesToMCTruthsTable(slices);
207  bestNuId =fBT->SliceToOrderedNuIdsByEnergy(nus, sEffPur);
208 
209 
210  for(size_t iClust = 0; iClust < clusters.size(); ++iClust){
211 
212  if(fObeyPreselection && rb::IsFiltered(evt, clusterHandle, iClust,
213  fPreselectionLabels)) continue;
214  // get record associated with the slice
215  std::vector< art::Ptr<caf::StandardRecord> > srs;
216  bool pass = true;
217  // Apply the CAF-level cuts
218  if(fApplyCAFCuts && fmpSR.isValid()){
219  srs = fmpSR.at(iClust);
220  if(srs.size() > 0) pass = fCAFCutter.passCuts(fCutType, srs[0].get());
221  }
222 
223  // if we don't pass the CAF cuts (when applying), skip this slice
224  if(!pass) continue;
225 
226  const rb::Cluster& cluster = clusters[iClust];
227  if(!fmPixelMap.isValid()) continue;
228 
229  const std::vector<art::Ptr<PixelMap>> pixelMaps = fmPixelMap.at(iClust);
230  if(pixelMaps.empty()) continue;
231 
233  if(bestNuId[iClust] != -1) truth = sEffPur[iClust][bestNuId[iClust]].neutrinoInt;
234  else continue;
235 
237  float nuEnergy, lepEnergy;
238  double vtxx, vtxy, vtxz;
239  interaction = SliceClassify(truth, cluster, nuEnergy, lepEnergy, vtxx, vtxy, vtxz);
240  finalstate = GetFinalStateType(truth);
241  finalstateprongs = GetFinalStateProngsType(truth);
242 
243  train = TrainingData(interaction, finalstate, finalstateprongs, nuEnergy, lepEnergy, *pixelMaps[0]);
244 
245  if(fWriteMapTH2) WriteMapTH2(evt, train.fPMap, 0, (unsigned int)iClust, 0);
246 
247 
248  if(fMapProngs){
249 
250  art::FindManyP<rb::Vertex> fmv(clusterHandle, evt, fVertexLabel);
251  art::FindManyP<rb::Prong> fmProng3D(clusterHandle, evt, art::InputTag(fProngModLabel, fProng3DLabel));
252  art::FindManyP<rb::Prong> fmProng2D(clusterHandle, evt, art::InputTag(fProngModLabel, fProng2DLabel));
253 
254  std::vector< art::Ptr<rb::Prong> > prongs3D;
255  std::vector< art::Ptr<rb::Prong> > prongs2D;
256  std::vector< art::Ptr<rb::Vertex> > vert;
257 
258  if(fmProng2D.isValid() ) prongs2D = fmProng2D.at(iClust);
259  if(fmProng3D.isValid() ) prongs3D = fmProng3D.at(iClust);
260  if(fmv.isValid() ) vert = fmv.at(iClust);
261  if(!fmProng2D.isValid() && !fmProng3D.isValid()) continue;
262 
263  art::FindManyP<PixelMap> fmPixelMap3D(prongs3D, evt, fPixelMapInput);
264  art::FindManyP<PixelMap> fmPixelMap2D(prongs2D, evt, fPixelMapInput);
265  art::FindOneP<rb::Shower> foSh(prongs3D, evt, fShowerLabel);
266 
267  TVector3 v;
268  if(!vert.empty() ) v = vert[0]->GetXYZ();
269  else v.SetXYZ(-999.,-999.,-999.);
270 
271 
272  double length;
273  TVector3 prongdir;
274  TVector3 beamdir = geom->NuMIBeamDirection();
275  float angle;
276 
277 
278  for(unsigned int iProng = 0; iProng < prongs3D.size(); ++iProng){
279 
280  length = -5.0;
281  angle = -5.0;
282 
283  if(!fmPixelMap3D.isValid()) continue;
284 
285  const std::vector< art::Ptr<PixelMap> > pixelMaps3D = fmPixelMap3D.at(iProng);
286  if(pixelMaps3D.empty()) continue;
287 
288  // Get gap to vertex from showerlid
289  double gap = -5.0;
290  if(foSh.isValid()){
291  cet::maybe_ref< art::Ptr<rb::Shower> const > roSh(foSh.at(iProng));
292  art::Ptr<rb::Shower> shower = roSh.ref();
293 
294  std::vector< art::Ptr<rb::Shower> > showers;
295  showers.push_back(shower);
296 
297  art::FindOneP<slid::ShowerLID> foShLID(showers, evt,fShowerLIDLabel);
298  if (foShLID.isValid()){
299  cet::maybe_ref< art::Ptr<slid::ShowerLID> const > roShLID(foShLID.at(0));
300  art::Ptr<slid::ShowerLID> showerLID = roShLID.ref();
301  gap = showerLID->Gap();
302  }
303  }
304 
305  // Currently not stored, for pm drawing only
306  ProngTrainingData p_train;
307  ProngType ptype3D, ptypeX, ptypeY;
308  double purity3D, purityX, purityY, recE;
309  unsigned int ncellX, ncellY;
310  bool isprimary;
311 
312  TrainingData pmtrain(interaction, finalstate, finalstateprongs,
313  nuEnergy, lepEnergy, *pixelMaps3D[0]);
314 
315  length = prongs3D[iProng]->TotalLength();
316  angle = prongdir*beamdir;
317 
318  ptype3D = ProngClassify(*prongs3D[iProng], &ptype3D, &ptypeX, &ptypeY,
319  &isprimary, &purity3D, &purityX, &purityY,
320  &recE, &ncellX, &ncellY);
321 
322  p_train = ProngTrainingData(ptype3D, ptypeX, ptypeY, isprimary,
323  purity3D, purityX, purityY, recE,
324  ncellX, ncellY, v.X(), v.Y(), v.Z(),
325  length, (double)angle, gap, *pixelMaps3D[0]);
326 
327 
328 
329  if(fWriteMapTH2) WriteMapTH2(evt, p_train.fProngPMap, 1, (unsigned int)iClust, iProng+1);
330 
333  fTrainTree->Fill();
334 
335  } // loop over 3D prongs
336 
337  for(unsigned int iProng2 = 0; iProng2 < prongs2D.size(); ++iProng2){
338 
339  if(!fmPixelMap2D.isValid()) continue;
340 
341  const std::vector< art::Ptr<PixelMap> > pixelMaps2D = fmPixelMap2D.at(iProng2);
342  if(pixelMaps2D.empty()) continue;
343 
344  ProngTrainingData p2_train;
345  ProngType ptype2D, ptypeX2D, ptypeY2D;
346  double purity2D, purityX2D, purityY2D, recE;
347  unsigned int ncellX, ncellY;
348  bool isprimary;
349 
350  ptype2D = ProngClassify(*prongs2D[iProng2], &ptype2D, &ptypeX2D, &ptypeY2D,
351  &isprimary, &purity2D, &purityX2D, &purityY2D,
352  &recE, &ncellX, &ncellY);
353 
354  p2_train = ProngTrainingData(ptype2D, ptypeX2D, ptypeY2D, isprimary,
355  purity2D, purityX2D, purityY2D, recE,
356  ncellX, ncellY, v.X(), v.Y(), v.Z(),
357  0.0, 0.0, 0.0, *pixelMaps2D[0]);
358 
360  fEventPMapsData->fProngLabeledPMs = p2_train;
361  fTrainTree->Fill();
362 
363  } // loop over 2D prongs
364  } // if(fMapProngs)
365 
366 
367  } // loop over slices
368  } // analyze
369 
370  //----------------------------------------------------------------------
371 
373  bool isProng, unsigned int slice, unsigned int id)
374  {
375 
376  TString mapid = "_r"+std::to_string(evt.run())+"_s"+std::to_string(evt.subRun())+"_e"
377  +std::to_string(evt.event())+"_sl"+std::to_string(slice)+"_p"+std::to_string(id);
378  TString data = "ProngPixelMap";
379  if (!isProng) data = "SlicePixelMap";
380 
381  TH2F* hist = pm.ToTH2(data+mapid);
382  TH2F* histX = pm.SingleViewToTH2(0, data+"ViewX"+mapid);
383  TH2F* histY = pm.SingleViewToTH2(1, data+"ViewY"+mapid);
384 
386 
387  TH2F* histWrite = tfs->make<TH2F>(*hist);
388  histWrite->Write();
389  TH2F* histWriteX = tfs->make<TH2F>(*histX);
390  histWriteX->Write();
391  TH2F* histWriteY = tfs->make<TH2F>(*histY);
392  histWriteY->Write();
393 
394  delete hist;
395  delete histWrite;
396  delete histX;
397  delete histWriteX;
398  delete histY;
399  delete histWriteY;
400 
401  }
402 
404 } // end namespace cvn
405 ////////////////////////////////////////////////////////////////////////
void reconfigure(const fhicl::ParameterSet &pset)
SubRunNumber_t subRun() const
Definition: Event.h:72
back track the reconstruction to the simulation
InteractionType SliceClassify(const art::Ptr< simb::MCTruth > truth, const rb::Cluster &slice, float &nuEnergy, float &lepEnergy, double &vtxx, double &vtxy, double &vtxz)
enum cvn::PType ProngType
Double_t angle
Definition: plot.C:86
FinalStateType GetFinalStateType(const art::Ptr< simb::MCTruth > truth)
TH2F * SingleViewToTH2(const unsigned int &view, TString name) const
Definition: PixelMap.cxx:312
The EventLabeledPMaps objects contains a PixelMap and the output class type, and any other bit that g...
std::vector< NeutrinoWithIndex > allMCTruth() const
float vtxx
A collection of associated CellHits.
Definition: Cluster.h:47
std::string fClusterPMLabel
Instance label for cluster pixelmaps.
enum cvn::Interaction InteractionType
recovalid::CAFCutter fCAFCutter
CAFCutter helper class.
DEFINE_ART_MODULE(TestTMapFile)
Defines an enumeration for prong classification.
Particle class.
std::vector< std::vector< cheat::NeutrinoEffPur > > SlicesToMCTruthsTable(const std::vector< const rb::Cluster * > &sliceList) const
Given ALL the slices in an event, including the noise slice, returns a vector of vector of structures...
std::string fCAFLabel
label for process that made standard records
Something else. Tau? Hopefully we don&#39;t use this.
PixelMap fPMap
PixelMap for the event.
Definition: TrainingData.h:69
const XML_Char const XML_Char * data
Definition: expat.h:268
std::string fProngPMLabel
Instance label for prong pixelmaps.
std::vector< int > SliceToOrderedNuIdsByEnergy(const std::vector< cheat::NeutrinoWithIndex > &nusWithIdx, const std::vector< std::vector< cheat::NeutrinoEffPur >> &slTruthTable) const
CVNEventProngDump(fhicl::ParameterSet const &pset)
TH2F * ToTH2(TString name) const
Definition: PixelMap.cxx:236
bool passCuts(int cut, const caf::StandardRecord *sr)
Definition: CAFCutter.cxx:37
length
Definition: demo0.py:21
std::string fShowerLabel
Instance label for showers.
std::string fVertexLabel
Instance label for vertex.
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
TrainingData fSliceLabeledPM
Class of the event.
Something else. Tau? Hopefully we don&#39;t use this.
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
float vtxy
FinalStateProngsType GetFinalStateProngsType(const art::Ptr< simb::MCTruth > truth)
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
Definition: FilterList.h:96
EventNumber_t event() const
Definition: Event.h:67
std::string fProng2DLabel
Instance label for input 2D prongs.
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
enum cvn::FinalState FinalStateType
Vertex location in position and time.
EventLabeledPMaps * fEventPMapsData
bool fApplyCAFCuts
do you want to apply CAF-level cuts?
The TrainingData objects contains a PixelMap and the output class type, and any other bit that goes i...
PixelMap fProngPMap
PixelMap for the event.
TVector3 NuMIBeamDirection() const
Direction of neutrinos from the NuMI beam (unit vector)
Something else. Tau? Hopefully we don&#39;t use this.
void analyze(const art::Event &evt)
float vtxz
void WriteMapTH2(const art::Event &evt, const PixelMap &pm, bool isProng, unsigned int slice, unsigned int id)
Function to extract TH2 from PixelMap and write to TFile.
T * make(ARGS...args) const
ProngTrainingData fProngLabeledPMs
PixelMap for the event.
void train()
Definition: train.C:24
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
The TrainingData objects contains a PixelMap and the output class type, and any other bit that goes i...
Definition: TrainingData.h:23
ProngType ProngClassify(const rb::Prong &prong, ProngType *pType3D, ProngType *pTypeX, ProngType *pTypeY, bool *isprimary, double *purity3D, double *purityX, double *purityY, double *recE, unsigned int *ncellX, unsigned int *ncellY)
PixelMap, basic input to CVN neural net.
Definition: PixelMap.h:23
void geom(int which=0)
Definition: geom.C:163
int fCutType
what cuts do you want to apply? (see CAFCutter.h)
std::string fShowerLIDLabel
Instance label for showers.
enum cvn::FinalStateProngs FinalStateProngsType
art::ServiceHandle< cheat::BackTracker > fBT
Helper class for Reco Validation modules.
Definition: CAFCutter.h:58
std::vector< std::string > fPreselectionLabels
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
RunNumber_t run() const
Definition: Event.h:77
std::string fProngModLabel
Module label for input prongs.
Encapsulate the geometry of one entire detector (near, far, ndos)
std::string fProng3DLabel
Instance label for input 3D prongs.
bool fMapProngs
Option to save maps for prongs.
enum BeamMode string