SkimmingUtils.h
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////
2 // \version $ $
3 //
4 // \brief Utility object to make skimming easier
5 //
6 // \author brebel@fnal.gov
7 //
8 ////////////////////////////////////////////////////////////////////////////
9 
10 #ifndef SKIMMINGUTILS_H
11 #define SKIMMINGUTILS_H
12 
13 #include <map>
14 
28 
29 #include "Utilities/AssociationUtil.h"
33 #include "Geometry/Geometry.h"
34 #include "RecoBase/Cluster.h"
35 #include "RecoBase/Prong.h"
36 #include "RecoBase/Track.h"
37 #include "RecoBase/Shower.h"
38 #include "RecoBase/Vertex.h"
39 #include "ShowerLID/ShowerLID.h"
40 #include "MCCheater/BackTracker.h"
41 #include "MCReweight/FluxWeights.h"
45 
46 namespace skim {
47  class SkimmingUtils;
48 }
49 
50 namespace skim {
51  class SkimmingUtils {
52 
53  public:
54 
55  SkimmingUtils();
56  SkimmingUtils(art::Handle<std::vector<rb::CellHit> > const& cellHandle,
57  std::string const& instance=std::string());
58  virtual ~SkimmingUtils();
59 
60  // method that fills a map of std::vector< T const* > to an index for a slice
61  template<class T> void FillSliceMap(size_t const& s,
62  art::FindMany<T> & fmp,
63  std::map<size_t, std::vector<T const*> > & ptrMap);
64 
65  // method that fills a std::vector< T const* > for a slice
66  template<class T> void FillSliceVector(size_t const& s,
67  art::FindMany<T> & fmp,
68  std::vector<T const*> & ptrMap);
69 
70  // method that fills a std::vector< T const* > for a slice based on another type
71  // U where there is a one to one association between T and U. Use this sparingly
72  // as this method creates a FindOne which is expensive - can't seem to find a
73  // better way of doing it though
74  template<class T, class U> void FillSliceVector(art::Event const& e,
75  std::string const& label,
76  std::vector<U const*> & otherPtrs,
77  std::vector<T const*> & fillPtrs);
78 
79  template<class T, class U> void FillSliceVector(art::Event const& e,
80  art::InputTag const& tag,
81  std::vector<U const*> & otherPtrs,
82  std::vector<T const*> & fillPtrs);
83 
84 
85  // method that copies the requested data products and creates an association between the
86  // product and the slice from which it came. The data product must not contain any kind
87  // of art::Ptrs or similar external references because we are using a copy constructor
88  // and those references will not copy correctly
89  template<class T> bool CopyProductAndSliceAssn(art::EDProducer const& prod,
90  art::Event & e,
91  std::vector<art::Ptr<T> > const& objIn,
92  std::vector<T> & objOut,
93  std::vector<rb::Cluster> & sliceOut,
94  art::Assns<rb::Cluster, T> & slcCROAssn);
95 
96  bool CopySlice(art::Ptr<rb::Cluster> const& sliceIn,
97  std::vector<rb::Cluster> & sliceOut,
98  size_t & sliceNum);
99 
100  bool CopyVertex(art::EDProducer const& prod,
101  art::Event & e,
102  art::Ptr<rb::Vertex> const& vertexIn,
103  std::vector<rb::Vertex> & vertexOut,
104  std::vector<rb::Cluster> & sliceOut,
106 
107  bool CopyShowers(art::EDProducer const& prod,
108  art::Event & e,
109  std::vector< art::Ptr<rb::Shower> > const& showerIn,
110  std::vector<rb::Shower> & showerOut,
111  std::vector<rb::Cluster> & sliceOut,
113 
114  bool CopyTracks(art::EDProducer const& prod,
115  art::Event & e,
116  std::vector< art::Ptr<rb::Track> > const& trackIn,
117  std::vector<rb::Track> & trackOut,
118  std::vector<rb::Cluster> & sliceOut,
120  std::string const& slcInstance="",
121  std::string const& trkInstance="");
122 
123  bool CopyProngs(art::EDProducer const& prod,
124  art::Event & e,
125  std::vector< art::Ptr<rb::Prong> > const& prongIn,
126  std::vector<rb::Prong> & prongOut,
127  std::vector<rb::Cluster> & sliceOut,
129 
130  bool SliceMCTruthAssociation(art::EDProducer const& prod,
131  art::Event & e,
132  size_t const& slcIdx,
133  std::vector<rb::Cluster> & sliceOut,
139  std::string const& fluxWeightLabel);
140 
141  private:
142 
143  std::map<std::pair<geo::CellUniqueId, float>, art::Ptr<rb::CellHit> > fCellIdToHitPtr;
146  };
147 
148 } // end namespace
149 
150 //------------------------------------------------------------------------------
152  std::vector<rb::Cluster> & sliceOut,
153  size_t & sliceNum)
154 {
155  if(fCellIdToHitPtr.size() < 1)
156  throw cet::exception("SkimmingUtils")
157  << "CellId/TNS to art::Ptr<rb::CellHit> map is empty, no way to copy the slice";
158 
160  // make a art::PtrVector<rb::CellHit> for all the hits in the slice
161  auto allCells = sliceIn->AllCells();
162 
163  std::set<std::pair<geo::CellUniqueId, float> > cuid;
164 
166  std::vector<double> slcWeights;
167 
168  // figure out which cells belong in the slice and the weight of each
169  for(size_t ch = 0; ch < allCells.size(); ++ch){
170  auto key = std::make_pair(geom->Plane(allCells[ch]->Plane())->Cell(allCells[ch]->Cell())->Id(),
171  allCells[ch]->TNS());
172 
173  // don't duplicate cells in the output list
174  if( !cuid.insert(key).second ){
175  LOG_WARNING("SkimmingUtils")
176  << "trying to insert a duplicate cell " << (*allCells[ch]);
177  continue;
178  }
179 
180  if( fCellIdToHitPtr.count(key) > 0 ){
181  slcCells.push_back( fCellIdToHitPtr.find(key)->second );
182  slcWeights.emplace_back( sliceIn->Weight(ch) );
183  }
184  }
185 
186  // make the new slice on the output collection
187  sliceOut.emplace_back(sliceIn->View(), sliceNum);
188  sliceOut.back().Add(slcCells, slcWeights);
189 
190  return true;
191 }
192 
193 //------------------------------------------------------------------------------
194 template<class T> inline void skim::SkimmingUtils::FillSliceMap(size_t const& s,
195  art::FindMany<T> & fmp,
196  std::map<size_t, std::vector<T const*> > & ptrMap)
197 {
198  std::vector<T const*> ptrs;
199  if( fmp.isValid() ){
200  fmp.get(s, ptrs);
201  ptrMap[s] = ptrs;
202  }
203 
204  return;
205 }
206 
207 //------------------------------------------------------------------------------
208 template<class T, class U> inline void skim::SkimmingUtils::FillSliceVector(art::Event const& e,
209  std::string const& label,
210  std::vector<U const*> & otherPtrs,
211  std::vector<T const*> & fillPtrs)
212 {
213  // make the FindOne
214  art::InputTag tag(label);
215 
216  this->FillSliceVector(e, tag, otherPtrs, fillPtrs);
217 
218  return;
219 }
220 
221 //------------------------------------------------------------------------------
222 template<class T, class U> inline void skim::SkimmingUtils::FillSliceVector(art::Event const& e,
223  art::InputTag const& tag,
224  std::vector<U const*> & otherPtrs,
225  std::vector<T const*> & fillPtrs)
226 {
227 
228  // make the FindOne
229  art::FindOne<T> fo(otherPtrs, e, tag);
230 
231  //clear the input vector
232  fillPtrs.clear();
233  fillPtrs.reserve(otherPtrs.size());
234 
235  // loop over the thesePtrs and get the associated others
236  for(size_t t = 0; t < otherPtrs.size(); ++t){
237  if(fo.at(t).isValid()) fillPtrs.emplace_back(&(fo.at(t).ref()));
238  }
239 
240  fillPtrs.shrink_to_fit();
241 
242  return;
243 }
244 
245 
246 //------------------------------------------------------------------------------
247 template<class T> inline void skim::SkimmingUtils::FillSliceVector(size_t const& s,
248  art::FindMany<T> & fmp,
249  std::vector<T const*> & ptrVec)
250 {
251  if( fmp.isValid() ) fmp.get(s, ptrVec);
252 
253  return;
254 }
255 
256 //------------------------------------------------------------------------------
258  art::Event & e,
259  std::vector<art::Ptr<T> > const& objIn,
260  std::vector<T> & objOut,
261  std::vector<rb::Cluster> & sliceOut,
262  art::Assns<rb::Cluster, T> & slcObjAssn)
263 {
264 
265  for(auto in : objIn){
266  objOut.push_back( *in );
267  util::CreateAssn(prod,
268  e,
269  sliceOut,
270  objOut,
271  slcObjAssn,
272  objOut.size()-1,
273  objOut.size(),
274  UINT_MAX,
275  fInstance,
276  fInstance);
277  }
278 
279  return true;
280 }
281 
282 //------------------------------------------------------------------------------
284  art::Event & e,
285  art::Ptr<rb::Vertex> const& vertexIn,
286  std::vector<rb::Vertex> & vertexOut,
287  std::vector<rb::Cluster> & sliceOut,
289 {
290  // I am assured there is only one vertex per slice by the elastic arms module
291  vertexOut.emplace_back(vertexIn->GetXYZ(), vertexIn->GetT() );
292  util::CreateAssn(prod,
293  e,
294  sliceOut,
295  vertexOut,
296  slcVtxAssn,
297  vertexOut.size()-1,
298  vertexOut.size(),
299  UINT_MAX,
300  fInstance,
301  fInstance);
302 
303  return true;
304 }
305 
306 //------------------------------------------------------------------------------
308  art::Event & e,
309  std::vector< art::Ptr<rb::Shower> > const& showerIn,
310  std::vector<rb::Shower> & showerOut,
311  std::vector<rb::Cluster> & sliceOut,
313 {
314  if(fCellIdToHitPtr.size() < 1)
315  throw cet::exception("SkimmingUtils")
316  << "CellId to art::Ptr<rb::CellHit> map is empty, no way to copy the showers";
317 
318  // loop over all the showers
321  for(auto shw : showerIn){
322 
323  allCells = shw->AllCells();
324 
325  std::vector<double> shwCellWt(shw->NCell(), 1.);
326  art::PtrVector<rb::CellHit> shwHitPtrs;
327 
328  // loop over all the hits in the shower to get their weights for the
329  // shower and to determine which of the Ptrs of the new slice collection
330  // are needed
331  for(size_t ch = 0; ch < allCells.size(); ++ch){
332  auto key = std::make_pair(geom->Plane(allCells[ch]->Plane())->Cell(allCells[ch]->Cell())->Id(),
333  allCells[ch]->TNS());
334 
335  if( fCellIdToHitPtr.count(key) > 0 ){
336  shwHitPtrs.push_back( fCellIdToHitPtr.find(key)->second );
337  shwCellWt[ch] = shw->Weight(ch);
338  }
339  } // end loop over cells in the current shower
340 
341  // make the copy of the shower
342  showerOut.emplace_back(shwHitPtrs,
343  shwCellWt,
344  shw->TotalLength(),
345  shw->Start(),
346  shw->Dir(),
347  shw->ID() );
348 
349  util::CreateAssn(prod,
350  e,
351  sliceOut,
352  showerOut,
353  slcShwAssn,
354  showerOut.size()-1,
355  showerOut.size(),
356  UINT_MAX,
357  fInstance,
358  fInstance);
359 
360  } // end loop over showers
361 
362  return true;
363 }
364 
365 //------------------------------------------------------------------------------
367  art::Event & e,
368  std::vector< art::Ptr<rb::Track> > const& trackIn,
369  std::vector<rb::Track> & trackOut,
370  std::vector<rb::Cluster> & sliceOut,
372  std::string const& slcInstance,
373  std::string const& trkInstance)
374 {
375 
376  std::string slcInst = (slcInstance.empty()) ? fInstance : slcInstance;
377  std::string trkInst = (trkInstance.empty()) ? fInstance : trkInstance;
378 
379  if(fCellIdToHitPtr.size() < 1)
380  throw cet::exception("SkimmingUtils")
381  << "CellId to art::Ptr<rb::CellHit> map is empty, no way to copy the tracks";
382 
383  // loop over all the showers
386  for(auto trk : trackIn){
387 
388  //////
389  // if( trk->View() != geo::kXorY) continue; // Copy Only 3D Tracks
390  /////
391  allCells = trk->AllCells();
392 
393  std::vector<double> trkCellWt(trk->NCell(), 1.);
394  art::PtrVector<rb::CellHit> trkHitPtrs;
395 
396  // loop over all the hits in the shower to get their weights for the
397  // shower and to determine which of the Ptrs of the new slice collection
398  // are needed
399  for(size_t ch = 0; ch < allCells.size(); ++ch){
400  auto key = std::make_pair(geom->Plane(allCells[ch]->Plane())->Cell(allCells[ch]->Cell())->Id(),
401  allCells[ch]->TNS());
402 
403  if( fCellIdToHitPtr.count(key) > 0 ){
404  trkHitPtrs.push_back( fCellIdToHitPtr.find(key)->second );
405  trkCellWt[ch] = trk->Weight(ch);
406  }
407  } // end loop over cells in the current shower
408 
409  // make the copy of the track
410  trackOut.emplace_back( trkHitPtrs,
411  trkCellWt,
412  trk->Trajectory(),
413  trk->Start(),
414  trk->Dir(),
415  trk->ID(),
416  trk->View());
417 
418  util::CreateAssn(prod,
419  e,
420  sliceOut,
421  trackOut,
422  slcTrkAssn,
423  trackOut.size()-1,
424  trackOut.size(),
425  UINT_MAX,
426  slcInst,
427  trkInst);
428 
429  } // end loop over tracks
430 
431  return true;
432 }
433 
434 
435 //------------------------------------------------------------------------------
437  art::Event & e,
438  std::vector< art::Ptr<rb::Prong> > const& prongIn,
439  std::vector<rb::Prong> & prongOut,
440  std::vector<rb::Cluster> & sliceOut,
442 {
443  if(fCellIdToHitPtr.size() < 1)
444  throw cet::exception("SkimmingUtils")
445  << "CellId to art::Ptr<rb::CellHit> map is empty, no way to copy the prongs";
446 
447  // loop over all the prongs
450  for(auto png : prongIn){
451 
452  allCells = png->AllCells();
453 
454  std::vector<double> pngCellWt(png->NCell(), 1.);
455  art::PtrVector<rb::CellHit> pngHitPtrs;
456 
457  // loop over all the hits in the prong to get their weights for the
458  // prong and to determine which of the Ptrs of the new slice collection
459  // are needed
460  for(size_t ch = 0; ch < allCells.size(); ++ch){
461  auto key = std::make_pair(geom->Plane(allCells[ch]->Plane())->Cell(allCells[ch]->Cell())->Id(),
462  allCells[ch]->TNS());
463 
464  if( fCellIdToHitPtr.count(key) > 0 ){
465  pngHitPtrs.push_back( fCellIdToHitPtr.find(key)->second );
466  pngCellWt[ch] = png->Weight(ch);
467  }
468  } // end loop over cells in the current prong
469 
470  // make the copy of the prong
471  prongOut.emplace_back(pngHitPtrs,
472  png->Start(),
473  png->Dir(),
474  png->ID() );
475 
476  util::CreateAssn(prod,
477  e,
478  sliceOut,
479  prongOut,
480  slcPrgAssn,
481  prongOut.size()-1,
482  prongOut.size(),
483  UINT_MAX,
484  fInstance,
485  fInstance);
486 
487  } // end loop over showers
488 
489  return true;
490 }
491 
492 //------------------------------------------------------------------------------
494  art::Event & e,
495  size_t const& slcIdx,
496  std::vector<rb::Cluster> & sliceOut,
502  std::string const& fluxWeightLabel)
503 {
504 
505  // have to get the BackTracker service to be able to make the truth to slice
506  // association
508 
509  // We want to get the slice truth, and to match the CAF analysis we use the
510  // energy matching metric
511  std::vector<rb::Cluster> slice({sliceOut[slcIdx]});
512  auto sliceTruthTable = bt->SlicesToMCTruthsTable(slice);
513 
514  if(sliceTruthTable[0].size() < 1){
515  LOG_WARNING("SkimmingUtils")
516  << "Could not find MCTruth objects for this slice in backtracker, so "
517  << "return without setting MCTruth or MCFlux associations for this slice";
518  return false;
519  }
520 
521  auto nuIdx = bt->SliceToOrderedNuIds(bt->allMCTruth(),
522  sliceTruthTable,
525 
526  // the first entry in nuIdx is supposed to be the best match of this slice
527  // to a true interaction there is only one slice in sliceTruthTable because
528  // we only passed one in the vector that made it
529  auto trueInt = sliceTruthTable[0][nuIdx[0]].neutrinoInt;
530 
531  util::CreateAssn(prod,
532  e,
533  sliceOut,
534  trueInt,
535  slcTruAssn,
536  slcIdx,
537  fInstance);
538 
539  // find the MCTruth to MCFlux association and associate the flux to the slice
541  mctp.push_back(trueInt);
542  art::FindOneP<simb::MCFlux> fofx(mctp, e, generatorLabel);
543  if(fofx.isValid()) util::CreateAssn(prod,
544  e,
545  sliceOut,
546  fofx.at(0),
547  slcFlxAssn,
548  slcIdx,
549  fInstance);
550 
551  art::FindOneP<fxwgt::FluxWeights> fofw(mctp, e, fluxWeightLabel);
552 
553  if(fofw.isValid()) util::CreateAssn(prod,
554  e,
555  sliceOut,
556  fofw.at(0),
557  slcFlxWgtAssn,
558  slcIdx,
559  fInstance);
560 
561  LOG_VERBATIM("SkimmingUtils")
562  << "The FindOne for the FluxWeights is valid: "
563  << fofw.isValid()
564  << " and there are "
565  << slcFlxWgtAssn.size()
566  << " associations "
567  << fofw.at(0)->GetHadronProductionCentralValue();
568 
569  art::FindOneP<simb::GTruth> fogt(mctp, e, generatorLabel);
570  if(fogt.isValid()) util::CreateAssn(prod,
571  e,
572  sliceOut,
573  fogt.at(0),
574  slcGTruAssn,
575  slcIdx,
576  fInstance);
577 
578 
579  return true;
580 }
581 
582 
583 #endif // SkimmingUtils
bool CopyProngs(art::EDProducer const &prod, art::Event &e, std::vector< art::Ptr< rb::Prong > > const &prongIn, std::vector< rb::Prong > &prongOut, std::vector< rb::Cluster > &sliceOut, art::Assns< rb::Cluster, rb::Prong > &slcPrgAssn)
virtual geo::View_t View() const
kXorY for 3D clusters.
Definition: Cluster.h:99
size_type size() const
Definition: Assns.h:440
bool CopyShowers(art::EDProducer const &prod, art::Event &e, std::vector< art::Ptr< rb::Shower > > const &showerIn, std::vector< rb::Shower > &showerOut, std::vector< rb::Cluster > &sliceOut, art::Assns< rb::Cluster, rb::Shower > &slcShwAssn)
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 FillSliceMap(size_t const &s, art::FindMany< T > &fmp, std::map< size_t, std::vector< T const * > > &ptrMap)
Atom< std::string > generatorLabel
std::vector< NeutrinoWithIndex > allMCTruth() const
bool SliceMCTruthAssociation(art::EDProducer const &prod, art::Event &e, size_t const &slcIdx, std::vector< rb::Cluster > &sliceOut, art::Assns< rb::Cluster, simb::MCTruth > &slcTruAssn, art::Assns< rb::Cluster, simb::MCFlux > &slcFlxAssn, art::Assns< rb::Cluster, fxwgt::FluxWeights > &slcFlxWgtAssn, art::Assns< rb::Cluster, simb::GTruth > &slcGTruAssn, std::string const &generatorLabel, std::string const &fluxWeightLabel)
size_type get(size_type i, reference item, data_reference data) const
Definition: FindMany.h:469
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
bool CopyProductAndSliceAssn(art::EDProducer const &prod, art::Event &e, std::vector< art::Ptr< T > > const &objIn, std::vector< T > &objOut, std::vector< rb::Cluster > &sliceOut, art::Assns< rb::Cluster, T > &slcCROAssn)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
const PlaneGeo * Plane(unsigned int i) const
std::string fInstance
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...
bool CopySlice(art::Ptr< rb::Cluster > const &sliceIn, std::vector< rb::Cluster > &sliceOut, size_t &sliceNum)
Module to create a summary of total POT seen in a job.
Definition: Evaluator.h:27
object containing MC flux information
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
Give every cell in the geometry a unique ID number based on the TGeo path to the node.
std::vector< int > SliceToOrderedNuIds(const std::vector< cheat::NeutrinoWithIndex > &nusWithIdx, const std::vector< std::vector< cheat::NeutrinoEffPur >> &slTruthTable, std::function< double(const cheat::NeutrinoEffPur &)> slMetric, std::function< double(const cheat::NeutrinoEffPur &)> nuMetric) const
Given a vector of indexed neutrino interaction and a vector of slice truth tables, returns a vector of neutrino interaction indices ordered by best match to the corresponding slice. Here, best match is determined according to the given cheat::NeutrinoEffPur functions for the slice and the nu.
void FillSliceVector(size_t const &s, art::FindMany< T > &fmp, std::vector< T const * > &ptrMap)
const char * label
Track finder for cosmic rays.
const XML_Char * s
Definition: expat.h:262
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
TVector3 GetXYZ() const
Definition: Vertex.cxx:45
double EnergyMetric(const cheat::NeutrinoEffPur &ep)
Function for NeutrinoEffPur&#39;s nu interaction to slice energy.
Vertex location in position and time.
T prod(const std::vector< T > &v)
Definition: prod.hpp:17
size_type size() const
Definition: PtrVector.h:308
#define LOG_WARNING(category)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
ifstream in
Definition: comparison.C:7
art::ServiceHandle< geo::Geometry > fGeom
void geom(int which=0)
Definition: geom.C:163
std::map< std::pair< geo::CellUniqueId, float >, art::Ptr< rb::CellHit > > fCellIdToHitPtr
bool CopyTracks(art::EDProducer const &prod, art::Event &e, std::vector< art::Ptr< rb::Track > > const &trackIn, std::vector< rb::Track > &trackOut, std::vector< rb::Cluster > &sliceOut, art::Assns< rb::Cluster, rb::Track > &slcTrkAssn, std::string const &slcInstance="", std::string const &trkInstance="")
#define LOG_VERBATIM(category)
double GetT() const
Definition: Vertex.h:26
Float_t e
Definition: plot.C:35
bool CopyVertex(art::EDProducer const &prod, art::Event &e, art::Ptr< rb::Vertex > const &vertexIn, std::vector< rb::Vertex > &vertexOut, std::vector< rb::Cluster > &sliceOut, art::Assns< rb::Cluster, rb::Vertex > &slcVtxAssn)
Definition: fwd.h:28
double Weight(unsigned int globalIdx) const
Weight assigned to the cell.
Definition: Cluster.cxx:209
Encapsulate the geometry of one entire detector (near, far, ndos)