LEMSummarizer_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file LEMSummarizer_module.cc
3 // \brief Transform slices into the minimal description %LEM needs to PID them
4 // \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
8 #include "LEM/LEMInput.h"
9 #include "LEM/Util.h"
10 
11 // NOvASoft includes
12 #include "Calibrator/Calibrator.h"
13 #include "RecoBase/CellHit.h"
14 #include "RecoBase/RecoHit.h"
15 #include "RecoBase/Cluster.h"
16 #include "RecoBase/FilterList.h"
17 #include "RecoBase/Vertex.h"
18 #include "Utilities/AssociationUtil.h"
19 
20 // Framework includes
26 #include "fhiclcpp/ParameterSet.h"
27 
28 #include <cmath>
29 
30 namespace lem
31 {
32  /// Transform slices into the minimal description %LEM needs to PID them
34  {
35  public:
36  explicit LEMSummarizer(const fhicl::ParameterSet& pset);
38 
39  virtual void reconfigure(const fhicl::ParameterSet& pset);
40  virtual void produce(art::Event& evt);
41 
42  protected:
43  LEMInput SliceToLEMInput(const art::Ptr<rb::Cluster>& slice, bool isMC,
44  int vtxPlane, int vtxCell, int vtxCellOther,
45  bool reverse) const;
46 
47  double fGeVToPECorr;
48 
51 
52  /// Labels of filter lists to obey
53  std::vector<std::string> fFilterLabels;
54 
56  };
57 
58  //......................................................................
60  {
61  reconfigure(pset);
62 
63  produces<std::vector<LEMInput>>();
64  produces<art::Assns<LEMInput, rb::Cluster>>();
65 
66  // Also produce variants where the z direction is reverse, if vertexing
67  // module provided a reversed vertex.
68  produces<std::vector<LEMInput>>("reverse");
69  produces<art::Assns<LEMInput, rb::Cluster>>("reverse");
70  }
71 
72  //......................................................................
74  {
75  }
76 
77  //......................................................................
79  {
80  fRawGeVForCalE = pset.get<bool>("RawGeVForCalE");
81 
82  fGeVToPECorr = pset.get<double>("GeVToPECorr");
83 
84  fSlicerLabel = pset.get<std::string>("SlicerLabel");
85  fVertexLabel = pset.get<std::string>("VertexLabel");
86  fFilterLabels = pset.get< std::vector<std::string> >("FilterLabels");
87  }
88 
89  //......................................................................
91  {
92  std::unique_ptr<std::vector<LEMInput>> sumcol(new std::vector<LEMInput>);
93  std::unique_ptr<art::Assns<LEMInput, rb::Cluster> > assns(new art::Assns<LEMInput, rb::Cluster>);
94 
95  std::unique_ptr<std::vector<LEMInput>> sumcol_rev(new std::vector<LEMInput>);
96  std::unique_ptr<art::Assns<LEMInput, rb::Cluster> > assns_rev(new art::Assns<LEMInput, rb::Cluster>);
97 
99  evt.getByLabel(fSlicerLabel, slices);
100 
102  const bool isMC = !evt.isRealData();
103 
104  art::FindManyP<rb::Vertex> fmv(slices, evt, fVertexLabel);
105  art::FindManyP<rb::Vertex> fmv_rev(slices, evt,
106  art::InputTag(fVertexLabel, "reverse"));
107 
108  const int sliceMax = slices->size();
109  for(int sliceIdx = 0; sliceIdx < sliceMax; ++sliceIdx){
110  if(rb::IsFiltered(evt, slices, sliceIdx, fFilterLabels))
111  continue;
112 
113  const art::Ptr<rb::Cluster> slice(slices, sliceIdx);
114 
115  if(slice->IsNoise()) continue;
116 
117  for(bool reverse: {false, true}){
118 
119  // Find the associated vertex and transform it into plane/cell
120  std::vector<art::Ptr<rb::Vertex> > vtxs;
121  if(reverse){
122  // It's OK if the vertexing module didn't produce a reversed vertex
123  if(!fmv_rev.isValid()) continue;
124  vtxs = fmv_rev.at(sliceIdx);
125  if(vtxs.empty()) continue;
126  }
127  else{
128  vtxs = fmv.at(sliceIdx);
129  }
130 
131  int vtxPlane, vtxCell, vtxCellOther;
132  if(vtxs.empty()){
133  DefaultVertex(*slice, vtxPlane, vtxCell, vtxCellOther);
134  }
135  else{
136  assert(vtxs.size() == 1);
137  VertexToPlaneAndCell(vtxs[0]->GetXYZ(), *slice,
138  vtxPlane, vtxCell, vtxCellOther, reverse);
139  }
140 
141  const LEMInput sum = SliceToLEMInput(slice, isMC,
142  vtxPlane, vtxCell, vtxCellOther,
143  reverse);
144 
145  if(reverse){
146  sumcol_rev->push_back(sum);
147  util::CreateAssn(*this, evt, *sumcol_rev, slice, *assns_rev);
148  }
149  else{
150  sumcol->push_back(sum);
151  util::CreateAssn(*this, evt, *sumcol, slice, *assns);
152  }
153 
154  } // end for reverse
155  } // end for sliceIdx
156 
157  evt.put(std::move(sumcol));
158  evt.put(std::move(assns));
159 
160  evt.put(std::move(sumcol_rev), "reverse");
161  evt.put(std::move(assns_rev), "reverse");
162  }
163 
164  //......................................................................
166  SliceToLEMInput(const art::Ptr<rb::Cluster>& slice, bool isMC,
167  int vtxPlane, int vtxCell, int vtxCellOther,
168  bool reverse) const
169  {
170  LEMInput sum;
171 
173 
174  // TODO: merge hits that fall on already hit cells. This apparently
175  // happens. Does it cause the self-energy to be wrong?
176  for(unsigned int cellIdx = 0; cellIdx < slice->NCell(); ++cellIdx){
177  const art::Ptr<rb::CellHit>& chit = slice->Cell(cellIdx);
178 
179  const int planeDelta = (reverse
180  ? vtxPlane - chit->Plane()
181  : chit->Plane() - vtxPlane);
182  const int cellDelta = ((abs(planeDelta)%2) ?
183  chit->Cell()-vtxCellOther :
184  chit->Cell()-vtxCell);
185 
186  // Align to standard position and clamp into allowed range
187  int plane = std::min(255, std::max(0, planeDelta + lem::kVertexPlane));
188  int cell = std::min(255, std::max(0, cellDelta + lem::kVertexCell ));
189 
190  const rb::RecoHit rhit = slice->RecoHit(chit);
191  if(rhit.IsCalibrated()){
192  // LEM libraries were based on PECorrs and generated for the FA.
193  // We need to get the GeV of each hit, and scale the GeV to the
194  // corresponding FA-era PECorr.
195  double gev = rhit.GeV();
196  double pecorr = gev*fGeVToPECorr;
197  assert(!std::isnan(pecorr));
198  assert(!std::isinf(pecorr));
199  sum.hits.push_back(LiteHit(plane, cell, pecorr, 0));
200  }
201  } // end for cellIdx
202 
203  if(!fRawGeVForCalE){
204  sum.totalGeV = slice->CalorimetricEnergy();
205  assert(!std::isnan(sum.totalGeV));
206  assert(!std::isinf(sum.totalGeV));
207  }
208  else{
209  sum.totalGeV = 0;
210  // Mean position in each view
211  const double wx = (slice->NXCell() > 0) ? slice->W(slice->XCell(0).get()) : 0;
212  const double wy = (slice->NYCell() > 0) ? slice->W(slice->YCell(0).get()) : 0;
213 
214  for(unsigned int hitIdx = 0; hitIdx < slice->NCell(); ++hitIdx){
215  // This is slow because we have to keep recalculating the mean
216  // positions. The manual version below is much better.
217  // const rb::RecoHit rhit = slice->RecoHit(hitIdx);
218 
219  const art::Ptr<rb::CellHit>& chit = slice->Cell(hitIdx);
220 
221  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, chit->View() == geo::kX ? wx : wy));
222 
223  if(!rhit.IsCalibrated()){
224  LOG_DEBUG("LEMSummarizer") << "Not calibrated?! "
225  << chit->Plane() << " " << chit->Cell()
226  << std::endl;
227  continue;
228  }
229  sum.totalGeV += rhit.GeV();
230  }
231  }
232 
233  return sum;
234  }
235 
236 
238 
239 } //namespace lem
240 ////////////////////////////////////////////////////////////////////////
T max(const caf::Proxy< T > &a, T b)
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
int isinf(const stan::math::var &a)
Definition: std_isinf.hpp:16
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.
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
Definition: Cluster.cxx:121
LEMSummarizer(const fhicl::ParameterSet &pset)
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
Definition: Cluster.cxx:157
void VertexToPlaneAndCell(const TVector3 vtx, const rb::Cluster &slice, int &plane, int &cell, int &cellOtherView, bool reverse)
Definition: Util.cxx:60
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
Transform slices into the minimal description LEM needs to PID them.
Vertical planes which measure X.
Definition: PlaneGeo.h:28
const int kVertexPlane
Definition: EventSummary.h:22
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
void abs(TH1 *hist)
bool isRealData() const
Definition: Event.h:83
std::vector< std::string > fFilterLabels
Labels of filter lists to obey.
DEFINE_ART_MODULE(TestTMapFile)
virtual void reconfigure(const fhicl::ParameterSet &pset)
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
PID
Definition: FillPIDs.h:14
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
double totalGeV
Definition: LEMInput.h:15
unsigned short Cell() const
Definition: CellHit.h:40
LEMInput SliceToLEMInput(const art::Ptr< rb::Cluster > &slice, bool isMC, int vtxPlane, int vtxCell, int vtxCellOther, bool reverse) const
void DefaultVertex(const rb::Cluster &slice, int &plane, int &cell, int &cellOtherView, bool reverse)
Definition: Util.cxx:20
double CalorimetricEnergy(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple estimate of neutrino energy.
Definition: Cluster.cxx:439
art::Ptr< rb::CellHit > YCell(unsigned int yIdx) const
Get the ith cell in the y-view.
Definition: Cluster.cxx:165
const int kVertexCell
Definition: EventSummary.h:23
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
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
Vertex location in position and time.
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
virtual void produce(art::Event &evt)
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
float GeV() const
Definition: RecoHit.cxx:69
T const * get() const
Definition: Ptr.h:321
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
std::vector< LiteHit > hits
Definition: LEMInput.h:14
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
void geom(int which=0)
Definition: geom.C:163
assert(nhit_max >=nhit_nbins)
T min(const caf::Proxy< T > &a, T b)
Simple representation of event for LEM use.
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
Compressed hit info, basic component of LEM events.
Definition: LiteHit.h:18
Double_t sum
Definition: plot.C:31
enum BeamMode string