RegCVNMapper_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file Reg CVNMapper_module.cc
3 // \brief Producer module for creating Regression CVN PixelMap objects, modified from CVNMapper_module
4 // \author Jianming Bian - bianjm@uci.edu
5 ////////////////////////////////////////////////////////////////////////
6 
7 // C/C++ includes
8 #include <iostream>
9 #include <sstream>
10 
11 // ROOT includes
12 #include "TTree.h"
13 #include "TH2F.h"
14 
15 // Framework includes
22 #include "fhiclcpp/ParameterSet.h"
27 
28 // NOvASoft includes
29 #include "RecoBase/CellHit.h"
30 #include "RecoBase/Cluster.h"
31 #include "RecoBase/FilterList.h"
32 #include "RecoBase/Vertex.h"
33 #include "RecoBase/Shower.h"
34 
35 #include "Geometry/Geometry.h"
37 #include "GeometryObjects/Geo.h"
38 
39 #include "Utilities/AssociationUtil.h"
43 
44 #include "Calibrator/Calibrator.h"
45 #include "MCCheater/BackTracker.h"
46 
48 #include "CVN/func/RegPixelMap.h"
49 #include "CVN/func/TrainingData.h"
50 
51 
52 
53 namespace cvn {
54 
55  class RegCVNMapper : public art::EDProducer {
56  public:
57  explicit RegCVNMapper(fhicl::ParameterSet const& pset);
58  ~RegCVNMapper();
59 
60  void produce(art::Event& evt);
61  void beginJob();
62  void endJob();
63 
64 
65 
66  private:
67  /// Module lablel for input clusters
69 
70  /// Use Prong or Track as input
72 
73  /// Module lablel for input prongs
75 
76  /// Instance lablel for input 3D prongs
78 
79  /// Instance lablel for input 2D prongs
81 
82  /// Instance lablel for input 3D tracks
84 
85  /// Instance lablel for input 2D tracks
87 
88 
89  /// Instance lablel for Regression CVN input Shower
91 
92  /// Instance lablel for cluster pixelmaps
94 
95  /// Instance lablel for prong pixelmaps
97 
98  /// Option to save maps for prongs
99  bool fMapProngs;
100 
101  /// Check rb::IsFiltered?
103 
104  /// Filter labels
105  std::vector<std::string> fPreselectionLabels;
106 
107 
108  /// Minimum number of hits for cluster to be converted to pixel map
109  unsigned short fMinClusterHits;
110 
111  /// Width of pixel map in cells
112  unsigned short fCellWidth;
113 
114  /// Length of pixel map in planes
115  unsigned short fPlaneLength;
116 
117  /// Maximum gap in planes at front of cluster to prevent pruning of upstream
118  /// hits
119  unsigned int fMaxPlaneGap;
120 
121  /// true=GeV, false=PECorr
122  bool fUseGeV;
123 
124  /// PixelMapProducer does the work for us
126 
127  };
128 
129 
130 
131  //.......................................................................
133  fClusterToken(consumes<std::vector<rb::Cluster>>(pset.get<std::string>("ClusterLabel"))),
134  fProngInput (pset.get<bool> ("ProngInput")),
135  fProngModLabel (pset.get<std::string> ("ProngModLabel")),
136  fProng3DLabel (pset.get<std::string> ("Prong3DLabel")),
137  fProng2DLabel (pset.get<std::string> ("Prong2DLabel")),
138  fTrack3DLabel (pset.get<std::string> ("Track3DLabel")),
139  fTrack2DLabel (pset.get<std::string> ("Track2DLabel")),
140  fShowerLabel (pset.get<std::string> ("ShowerLabel")),
141  fClusterPMLabel(pset.get<std::string> ("ClusterPMLabel")),
142  fProngPMLabel (pset.get<std::string> ("ProngPMLabel")),
143  fMapProngs (pset.get<bool> ("MapProngs")),
144  fObeyPreselection(pset.get<bool> ("ObeyPreselection")),
145  fPreselectionLabels (pset.get<std::vector<std::string>>("PreselectionLabels")),
146  fMinClusterHits(pset.get<unsigned short> ("MinClusterHits")),
147  fCellWidth (pset.get<unsigned short> ("CellWidth")),
148  fPlaneLength (pset.get<unsigned short> ("PlaneLength")),
149  fMaxPlaneGap (pset.get<unsigned int> ("MaxPlaneGap")),
150  fUseGeV (pset.get<bool> ("UseGeV")),
152  {
153 
154  produces< std::vector<cvn::RegPixelMap> >(fClusterPMLabel);
155  produces< std::vector<cvn::RegPixelMap> >(fProngPMLabel);
156  produces< art::Assns<cvn::RegPixelMap, rb::Cluster> >();
157  produces< art::Assns<cvn::RegPixelMap, rb::Prong> >();
158  }
159 
160  //......................................................................
162  {
163  //======================================================================
164  // Clean up any memory allocated by your module
165  //======================================================================
166  }
167 
168  //......................................................................
170  {
171 }
172 
173  //......................................................................
175  {
176  }
177 
178  //......................................................................
180  {
181 
182 
184  evt.getByToken(fClusterToken,sliceHandle);
185 
186  if(sliceHandle->empty()){
187  mf::LogWarning ("No Slices")<<"No Slices in the input file";
188  return;
189  }
190 
192  for(unsigned int i = 0; i<sliceHandle->size();++i){
193  art::Ptr<rb::Cluster> clust(sliceHandle,i);
194  sliceCol.push_back(clust);
195  }
196 
197 
198  //Declaring containers for things to be stored in event
199  std::unique_ptr< std::vector<RegPixelMap> >
200  pmCol(new std::vector<RegPixelMap>);
201  std::unique_ptr< std::vector<RegPixelMap> >
202  prong_pmCol(new std::vector<RegPixelMap>);
203  std::unique_ptr< art::Assns<RegPixelMap, rb::Cluster> >
205  std::unique_ptr< art::Assns<RegPixelMap, rb::Prong> >
206  prong_assoc(new art::Assns<RegPixelMap,
207  rb::Prong>);
208 
210 
211  art::FindManyP<rb::Prong> fmProng3D(sliceHandle, evt,
213 
214  art::FindManyP<rb::Prong> fmProng2D(sliceHandle, evt,
216 
217  art::FindManyP<rb::Track> fmTrack3D(sliceHandle, evt, fTrack3DLabel);
218 
219  art::FindManyP<rb::Track> fmTrack2D(sliceHandle, evt, fTrack2DLabel);
220 
221  art::FindManyP<rb::Vertex> fmv(sliceHandle, evt, "elasticarmshs");
222 
223  for(unsigned int iSlice = 0; iSlice < sliceCol.size(); ++iSlice)
224  {
225  if (sliceCol[iSlice]->IsNoise() ||
226  sliceCol[iSlice]->NCell() < fMinClusterHits)
227  continue;
228  // Skip if we're obeying preselection
229  if(fObeyPreselection && rb::IsFiltered(evt, sliceHandle, iSlice, fPreselectionLabels))
230  continue;
231 // std::cout<<"CVN CreatMap, iSlice "<<iSlice<<std::endl;
232  RegPixelMap pm = fProducer.CreateMap(*sliceCol[iSlice]);
233  Boundary bound = pm.Bound();
234 
235 // pmCol->push_back(pm);
236 // util::CreateAssn(*this, evt, *(pmCol.get()),
237 // sliceCol[iSlice], *(assoc.get()), UINT_MAX,
238 // fClusterPMLabel);
239 
240  // Option to save prong level pixel maps
241  if ( fMapProngs ){
242  unsigned int vtxPlane = 0;
243  unsigned int vtxXCell = 0;
244  unsigned int vtxYCell = 0;
245 
246  if ((fmv.isValid())){
247  std::vector<art::Ptr<rb::Vertex> > v = fmv.at(iSlice);
248  for (unsigned int j=0; j<v.size(); ++j) {
249  TVector3 evtVtx(0,0,0);
250  //store vertex coordinates
251  evtVtx.SetX(v[j]->GetX());
252  evtVtx.SetY(v[j]->GetY());
253  evtVtx.SetZ(v[j]->GetZ());
254  double dmin = 9999999;
255  //find plane number closest to vertex
256  for(unsigned int pp = 0; pp < geom->NPlanes();pp++){
257  if(pp%2==0)continue;
258  Double_t cellXYZ[3]={0,0,0};
259  geom->Plane(pp)->Cell(0)->GetCenter(cellXYZ);
260  double mind = fabs(cellXYZ[2]-evtVtx[2]);
261  if(mind<=dmin){
262  dmin=mind;
263  vtxPlane = pp;
264  }
265  }
266 // if(vtxPlane%2==1)vtxPlane = vtxPlane-1;
267  double dminx = 9999999;
268  double dminy = 9999999;
269  unsigned int xp=0;
270  unsigned int yp=1;
271  if(geom->Plane(0)->View()==geo::kY){xp=1;yp=0;}
272 
273  for(unsigned int xx = 0; xx < geom->Plane(xp)->Ncells();xx++){
274  Double_t cellXYZ[3]={0,0,0};
275  geom->Plane(xp)->Cell(xx)->GetCenter(cellXYZ);
276  double mind = fabs(cellXYZ[0]-evtVtx[0]);
277  if(mind<=dminx){
278  dminx=mind;
279  vtxXCell = xx;
280  }
281  }
282 
283  for(unsigned int yy = 0; yy < geom->Plane(yp)->Ncells();yy++){
284  Double_t cellXYZ[3]={0,0,0};
285  geom->Plane(yp)->Cell(yy)->GetCenter(cellXYZ);
286  double mind = fabs(cellXYZ[1]-evtVtx[1]);
287  if(mind<=dminy){
288  dminy=mind;
289  vtxYCell = yy;
290  }
291  }
292  }
293  }
294 // std::cout<<"In PixelMapper - Evt: Run, SubRun, Event, SubEvent "<<evt.run()<<" "<<evt.subRun()<<" "<<evt.id().event()<<" "<<iSlice<<std::endl;
295  RegPixelMap pm_evt = fProducer.CreateMapGivenVertex(*sliceCol[iSlice],
296  bound, vtxPlane, vtxXCell, vtxYCell);
297 
298 // std::cout<<"Finish fProducer.CreateMapGivenVertex **************************************************************"<<std::endl;
299 // pmCol->push_back(pm_evt);
300 // util::CreateAssn(*this, evt, *(pmCol.get()),
301 // sliceCol[iSlice], *(assoc.get()), UINT_MAX,
302 // fClusterPMLabel);
303 
304 
305  std::vector<art::Ptr<rb::Prong>> prongs3D;
306  std::vector<art::Ptr<rb::Prong>> prongs2D;
307  std::vector<art::Ptr<rb::Shower>> showers3D;
308 // std::cout<<"fmProng2D.isValid(), fmProng3D.isValid() "<<fmProng2D.isValid()<<" "<<fmProng3D.isValid()<<std::endl;
309  if(fProngInput==true && fmProng2D.isValid() ){
310  prongs2D = fmProng2D.at(iSlice);
311  }
312  if(fProngInput==false && fmTrack2D.isValid() ){
313  std::vector<art::Ptr<rb::Track> > tracks2D;
314  tracks2D.clear();
315  tracks2D = fmTrack2D.at(iSlice);
316  for (unsigned int ip=0; ip<tracks2D.size(); ++ip){
317  prongs2D.push_back(tracks2D[ip]);
318  }
319  }
320 
321  if(fProngInput==true && fmProng3D.isValid() ){
322  prongs3D = fmProng3D.at(iSlice);
323  }
324  if(fProngInput==false && fmTrack3D.isValid() ){
325  std::vector<art::Ptr<rb::Track> > tracks3D;
326  tracks3D.clear();
327  tracks3D = fmTrack3D.at(iSlice);
328  for (unsigned int ip=0; ip<tracks3D.size(); ++ip){
329  prongs3D.push_back(tracks3D[ip]);
330  }
331  }
332 
333  if(!prongs3D.empty()){
334  unsigned int pngId = 0;
335  art::FindOneP<rb::Shower> fmShower3D(prongs3D, evt, fShowerLabel);
336  for(unsigned int iPng = 0; iPng < prongs3D.size(); ++iPng){
337  if(fmShower3D.isValid()){
338  art::Ptr<rb::Shower> slid = fmShower3D.at(pngId);
339  if(slid){
340  showers3D.push_back(slid);
341  }
342  }
343  pngId++;
344  }
345  }else continue;
346  if(prongs3D.size()!=showers3D.size()) continue;
347 
348 // std::cout<<"prongs3D.size(), prongs2D.size, shower3D.size() "<<prongs3D.size()<<" "<<prongs2D.size()<<" "<<showers3D.size()<<std::endl;
349 // RegPixelMap pm_evt_sh = fProducer.CreateMapGivenVertex(showers3D,
350 // bound, vtxPlane, vtxXCell, vtxYCell);
351 // pmCol->push_back(pm_evt_sh);
352  pmCol->push_back(pm_evt);
353  util::CreateAssn(*this, evt, *(pmCol.get()),
354  sliceCol[iSlice], *(assoc.get()), UINT_MAX,
356 
357  for( unsigned int iProng = 0; iProng < showers3D.size(); ++iProng ){
358 // std::cout<<"CreateMapGivenBoundary 3D Prong, NCell, NPlane "<<prongs3D[iProng]->NCell()<<" "<<prongs3D[iProng]->ExtentPlane()<<std::endl;
359 // std::cout<<"CreateMapGivenBoundary 3D Shwoer, NCell, NPlane "<<showers3D[iProng]->NCell()<<" "<<showers3D[iProng]->ExtentPlane()<<std::endl;
360 // std::cout<<"vtxPlane, vtxXCell, vtxYCell "<<vtxPlane<<" "<<vtxXCell<<" "<<vtxYCell<<std::endl;
361  RegPixelMap new_pm = fProducer.CreateMapGivenShowerVertex(*showers3D[iProng],
362  bound, vtxPlane, vtxXCell, vtxYCell);
363 
364 
365  prong_pmCol->push_back(new_pm);
366  util::CreateAssn(*this, evt, *(prong_pmCol.get()),
367  prongs3D[iProng], *(prong_assoc.get()), UINT_MAX, fProngPMLabel);
368  }// 3D prongs
369 
370  for( unsigned int iProng2 = 0; iProng2 < prongs2D.size(); ++iProng2 ){
371 
372  // make prong maps with the same boundary as the slice maps
373  RegPixelMap new_pm = fProducer.CreateMapGivenBoundary(*prongs2D[iProng2],
374  bound);
375 
376  prong_pmCol->push_back(new_pm);
377  util::CreateAssn(*this, evt, *(prong_pmCol.get()),
378  prongs2D[iProng2], *(prong_assoc.get()), UINT_MAX, fProngPMLabel);
379  }// 2D prongs
380 
381  }// MapProngs
382 
383  }
384  evt.put(std::move(pmCol), fClusterPMLabel);
385  evt.put(std::move(prong_pmCol), fProngPMLabel);
386  evt.put(std::move(assoc));
387  evt.put(std::move(prong_assoc));
388  // evt.put(std::move(prong_assoc),fProng3DLabel);
389 
390  }
391 
392  //----------------------------------------------------------------------
393 
394 
395 
397 } // end namespace cvn
398 ////////////////////////////////////////////////////////////////////////
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.
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
Double_t xx
Definition: macro.C:12
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
bool fProngInput
Use Prong or Track as input.
Producer algorithm for RegPixelMap, input to CVN neural net.
std::string fClusterPMLabel
Instance lablel for cluster pixelmaps.
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
std::vector< std::string > fPreselectionLabels
Filter labels.
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
std::string fProng3DLabel
Instance lablel for input 3D prongs.
unsigned short fPlaneLength
Length of pixel map in planes.
TString ip
Definition: loadincs.C:5
RegPixelMapProducer fProducer
PixelMapProducer does the work for us.
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
Defines an enumeration for prong classification.
Particle class.
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
std::string fProng2DLabel
Instance lablel for input 2D prongs.
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
Boundary Bound() const
Map boundary.
Definition: RegPixelMap.h:39
RegPixelMap CreateMapGivenBoundary(const rb::Cluster &cluster, const Boundary &bound)
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
RegPixelMap CreateMapGivenShowerVertex(const rb::Shower &shower, const Boundary &bound, unsigned int vtxPlane, unsigned int vtxXCell, unsigned int vtxYCell)
bool fMapProngs
Option to save maps for prongs.
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
RegPixelMapProducer for CVN.
std::string fTrack2DLabel
Instance lablel for input 2D tracks.
unsigned short fCellWidth
Width of pixel map in cells.
int evt
RegCVNMapper(fhicl::ParameterSet const &pset)
bool fUseGeV
true=GeV, false=PECorr
Collect Geo headers and supply basic geometry functions.
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
const double j
Definition: BetheBloch.cxx:29
double GetY(int dcm=1, int feb=0, int pix=0)
Definition: PlotOnMon.C:119
std::string fProngPMLabel
Instance lablel for prong pixelmaps.
Vertex location in position and time.
Perform a "2 point" Hough transform on a collection of hits.
unsigned short fMinClusterHits
Minimum number of hits for cluster to be converted to pixel map.
bool fObeyPreselection
Check rb::IsFiltered?
RegPixelMap, basic input to CVN neural net.
Definition: RegPixelMap.h:22
size_type size() const
Definition: PtrVector.h:308
RegPixelMap CreateMapGivenVertex(const rb::Cluster &cluster, const Boundary &bound, unsigned int vtxPlane, unsigned int vtxXCell, unsigned int vtxYCell)
std::string fShowerLabel
Instance lablel for Regression CVN input Shower.
A Cluster with defined start position and direction.
Definition: Prong.h:19
void produce(art::Event &evt)
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void geom(int which=0)
Definition: geom.C:163
const art::ProductToken< std::vector< rb::Cluster > > fClusterToken
Module lablel for input clusters.
std::string fTrack3DLabel
Instance lablel for input 3D tracks.
std::string fProngModLabel
Module lablel for input prongs.
RegPixelMap CreateMap(const rb::Cluster &slice)
Build slid::LID objects to store electron ID, if asked for, otherwise, calculate LID info and make av...
Definition: FillPIDs.h:13
unsigned int NPlanes() const
bool getByToken(ProductToken< PROD > const &token, Handle< PROD > &result) const
Definition: DataViewImpl.h:387
RegPixelMap for Regression modified from Dominick&#39;s PixelMap.
ProductToken< T > consumes(InputTag const &)
Encapsulate the geometry of one entire detector (near, far, ndos)
double GetX(int ndb=14, int db=1, int feb=0, int pix=0)
Definition: PlotOnMon.C:111
enum BeamMode string