SliceLIDBuilder_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief Make LL based NueCC ID with Tensorflow implementation of LSTM
3 ///
4 /// \author Andrew Vold voldx034@umn.edu
5 ////////////////////////////////////////////////////////////////////////
6 
7 #include "TMath.h"
8 
9 // NOvA includes
10 #include "RecoBase/Cluster.h"
11 #include "RecoBase/PID.h"
12 #include "RecoBase/FilterList.h"
13 #include "RecoBase/RecoHit.h"
14 #include "Utilities/AssociationUtil.h"
16 #include "ShowerLID/SliceLID.h"
17 #include "ShowerLID/ShowerLID.h"
18 #include "ShowerLID/NuEEnergyAlg.h"
20 #include "Geometry/Geometry.h"
21 
22 // ART includes
34 #include "fhiclcpp/ParameterSet.h"
35 
36 namespace slid {
37 
39  public:
40  explicit SliceLIDBuilder(const fhicl::ParameterSet& pset);
42 
43  virtual void produce(art::Event& evt);
44  virtual void reconfigure(const fhicl::ParameterSet& pset);
45  virtual void beginRun(art::Run& run);
46 
47  protected:
51  std::vector<std::string> fFilterLabels;
52 
55 
60  };
62 } // namespace
63 
64 namespace slid {
65 
67  : fTF(0),fNuEEnergyAlg(0),
68  fTFPSet(pset.get< fhicl::ParameterSet >("TensorflowHandler")),
69  fNuEEnergyAlgPSet(pset.get< fhicl::ParameterSet >("NuEEnergyAlgPSet"))
70  {
71  reconfigure(pset);
72  // Define output structures produced.
73  produces<std::vector<slid::SliceLID> >();
74  produces< art::Assns<slid::SliceLID, rb::Cluster> >();
75  }
76 
77  //......................................................................
78 
80  {
81  if(fTF) delete fTF;
82  if (fNuEEnergyAlg) delete fNuEEnergyAlg;
83  }
84 
85  //......................................................................
86 
88  {
89  fTFPSet = pset.get< fhicl::ParameterSet >("TensorflowHandler");
90  fNuEEnergyAlgPSet = pset.get< fhicl::ParameterSet >("NuEEnergyAlgPSet");
91  fSlicerLabel = pset.get< std::string >("SlicerLabel");
92  fShowerLabel = pset.get< std::string >("ShowerLabel");
93  fShowerLIDLabel = pset.get< std::string >("ShowerLIDLabel");
94  fObeyPreselection = pset.get<bool>("ObeyPreselection");
95  fSkipNoiseSlices = pset.get<bool>("SkipNoiseSlices");
96  fFilterLabels = pset.get< std::vector<std::string> >("FilterLabels");
97  }
98 
100  if(fTF) delete fTF;
101  if (fNuEEnergyAlg) delete fNuEEnergyAlg;
104  }
105  //......................................................................
106 
108  {
109 
110  LOG_DEBUG("SliceLIDBuilderDEBUG") << "In SliceLIDBuilder::produce" << '\n';
111 
112  std::unique_ptr< std::vector<slid::SliceLID> > slicelidcol(new std::vector<slid::SliceLID>);
113  std::unique_ptr< art::Assns<slid::SliceLID, rb::Cluster> > assnSliceLidSlice(new art::Assns<slid::SliceLID, rb::Cluster>);
114 
116  evt.getByLabel(fSlicerLabel, slicecol);
117  art::PtrVector<rb::Cluster> slicelist;
119  for (unsigned int i = 0; i < slicecol->size(); ++i) {
120  art::Ptr<rb::Cluster> slice(slicecol, i);
121  slicelist.push_back(slice);
122  }
123 
124  // loop over all of the slices
125  for (unsigned int iSlice = 0; iSlice < slicelist.size(); ++iSlice) {
126 
127  if(fSkipNoiseSlices && slicelist[iSlice]->IsNoise()) continue;
128  if(fObeyPreselection && rb::IsFiltered(evt, slicecol, iSlice, fFilterLabels)) continue;
129 
130  std::vector<art::Ptr<slid::ShowerLID> > slidcol;
131  slidcol.clear();
132  art::FindManyP<slid::ShowerLID> fmslid(slicecol, evt, fShowerLIDLabel);
133  if(fmslid.isValid()){
134  slidcol = fmslid.at(iSlice);
135  }
136 
137  art::FindOneP<rb::Shower> fos(slidcol, evt, fShowerLIDLabel);
138  std::sort(slidcol.begin(),slidcol.end(),CompareByE);
139  std::vector<art::Ptr<rb::Shower>> showercol;
140 
141  if(slidcol.size()==0) continue;
142 
143  for(unsigned int i = 0; i < slidcol.size(); i++){
144  cet::maybe_ref<art::Ptr<rb::Shower> const> rshw(fos.at(i));
145  art::Ptr<rb::Shower> shw = rshw.ref();
146  showercol.push_back(shw);
147  }
148 
149  //Form feature map for evaluation
150  float sliceE = slicelist[iSlice]->TotalGeV();
151  std::map<int,std::vector<float>> featuremap;
152  std::vector<float> shwvec;
153 
154  for(int i = 0; i < 7; i++){
155  shwvec.clear();
156  if(i < (int)showercol.size()){
157  art::Ptr<rb::Shower> shw = showercol[i];
158  float Ex = 0;
159  float Ey = 0;
160  for(int j = 0; j < (int)shw->NCell(); j++){
162  rb::RecoHit rhit = shw->RecoHit(hit);
163  if(!rhit.IsCalibrated()) continue;
164  if(hit->View()==geo::kX) Ex += rhit.GeV();
165  else if(hit->View()==geo::kY) Ey += rhit.GeV();
166  }
167  art::Ptr<slid::ShowerLID> slid = slidcol[i];
168  float E = shw->TotalGeV();
169  float len = shw->TotalLength();
170  float Eview = Ex/Ey;
171  double theta = shw->Dir().Angle(geom->NuMIBeamDirection());
172  float cos = TMath::Cos(theta);
173  float gap = slid->Gap();
174  float miss = shw->NMissingPlanes(geo::kXorY);
175  float emulll = slid->EMuLLL();
176  float epi0lll = slid->EPi0LLL();
177  float eplll = slid->EPLLL();
178  float enlll = slid->ENLLL();
179  float epilll = slid->EPiLLL();
180  float eglll = slid->EGLLL();
181  float emullt = slid->EMuLLT();
182  float epi0llt = slid->EPi0LLT();
183  float epllt = slid->EPLLT();
184  float enllt = slid->ENLLT();
185  float epillt = slid->EPiLLT();
186  float egllt = slid->EGLLT();
187  shwvec.push_back(E);
188  shwvec.push_back(len);
189  shwvec.push_back(Eview);
190  shwvec.push_back(cos);
191  shwvec.push_back(gap);
192  shwvec.push_back(miss);
193  shwvec.push_back(emulll);
194  shwvec.push_back(epi0lll);
195  shwvec.push_back(eplll);
196  shwvec.push_back(enlll);
197  shwvec.push_back(epilll);
198  shwvec.push_back(eglll);
199  shwvec.push_back(emullt);
200  shwvec.push_back(epi0llt);
201  shwvec.push_back(epllt);
202  shwvec.push_back(enllt);
203  shwvec.push_back(epillt);
204  shwvec.push_back(egllt);
205  }
206  else{
207  for(int j = 0; j < 19; j++){shwvec.push_back(0.);}
208  }
209  featuremap[i] = shwvec;
210  }
211 
212  float score = fTF->Predict(sliceE,featuremap,false);
213  slid::SliceLID slicelid;
214  slicelid.SetValue(score);
215  slicelidcol->push_back(slicelid);
216  util::CreateAssn(*this,evt,*(slicelidcol.get()),slicelist[iSlice],*(assnSliceLidSlice.get()));
217  //create new associations for the most recent set of showers
218 
219  } //end loop over slices
220 
221  evt.put(std::move(slicelidcol));
222  evt.put(std::move(assnSliceLidSlice));
223 
224  }// end of producer
225 
226  ////////////////////////////////////////////////////////////////////////
228  {
229  return a->ShowerEnergy() > b->ShowerEnergy();
230  }
231 
232 
234 }
fhicl::ParameterSet fNuEEnergyAlgPSet
const XML_Char int len
Definition: expat.h:262
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
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.
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
float ENLLT() const
Definition: ShowerLID.h:114
float EGLLL() const
Definition: ShowerLID.h:105
pdg code and pid value
float EPiLLT() const
Definition: ShowerLID.h:116
X or Y views.
Definition: PlaneGeo.h:30
geo::View_t View() const
Definition: CellHit.h:41
Vertical planes which measure X.
Definition: PlaneGeo.h:28
virtual void beginRun(art::Run &run)
bool CompareByE(const slid::ShowerLID &a, const slid::ShowerLID &b)
float EPLLL() const
Definition: ShowerLID.h:111
std::vector< std::string > fFilterLabels
DEFINE_ART_MODULE(TestTMapFile)
fhicl::ParameterSet fTFPSet
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
virtual void reconfigure(const fhicl::ParameterSet &pset)
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
float EPi0LLT() const
Definition: ShowerLID.h:110
virtual double TotalLength() const
Length (cm) of a shower.
Definition: Shower.cxx:43
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
TensorflowHandler for ShowerLID.
SliceLIDBuilder(const fhicl::ParameterSet &pset)
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
Float_t E
Definition: plot.C:20
const double a
T get(std::string const &key) const
Definition: ParameterSet.h:231
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
float EMuLLL() const
Definition: ShowerLID.h:107
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
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
float EPi0LLL() const
Definition: ShowerLID.h:109
Definition: run.py:1
size_type size() const
Definition: PtrVector.h:308
float EPiLLL() const
Definition: ShowerLID.h:115
TVector3 NuMIBeamDirection() const
Direction of neutrinos from the NuMI beam (unit vector)
int NMissingPlanes(geo::View_t view) const
Total number of missing planes in cluster.
Definition: Cluster.cxx:693
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
virtual void produce(art::Event &evt)
float GeV() const
Definition: RecoHit.cxx:69
float ShowerEnergy() const
Definition: ShowerLID.h:157
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Calculates deposited and corrected energy of the electron shower and of electron flavoured neutrino...
Definition: event.h:1
Wrapper for Tensorflow which handles construction and prediction.
void geom(int which=0)
Definition: geom.C:163
const hit & b
Definition: hits.cxx:21
float Predict(float sliceE, std::map< int, std::vector< float >> featuremap, bool fRHC)
float EMuLLT() const
Definition: ShowerLID.h:108
float ENLLL() const
Definition: ShowerLID.h:113
T cos(T number)
Definition: d0nt_math.hpp:78
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
Definition: Cluster.cxx:378
Build slid::LID objects to store electron ID, if asked for, otherwise, calculate LID info and make av...
Definition: FillPIDs.h:12
tensorflow::TensorflowHandler * fTF
float Gap() const
Definition: ShowerLID.h:117
void SetValue(float in)
Definition: SliceLID.h:24
float EGLLT() const
Definition: ShowerLID.h:106
void clear()
Definition: PtrVector.h:537
Definition: fwd.h:28
Encapsulate the geometry of one entire detector (near, far, ndos)
float EPLLT() const
Definition: ShowerLID.h:112