ProngCheater_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief ProngCheater module
3 ///
4 /// \author brebel@fnal.gov
5 ///
6 ////////////////////////////////////////////////////////////////////////
7 
8 #include <string>
9 #include <vector>
10 
11 // NOvASoft includes
12 #include "MCCheater/BackTracker.h"
13 #include "RecoBase/CellHit.h"
14 #include "RecoBase/Prong.h"
16 #include "Utilities/AssociationUtil.h"
17 
18 // Framework includes
20 #include "fhiclcpp/ParameterSet.h"
29 
30 namespace cheat {
31  class ProngCheater : public art::EDProducer {
32  public:
33  explicit ProngCheater(fhicl::ParameterSet const& pset);
34  virtual ~ProngCheater();
35 
36  void produce(art::Event& evt);
37 
38  void reconfigure(fhicl::ParameterSet const& pset);
39 
40  private:
41 
42  std::string fCheatedClusterLabel; ///< label for module creating recob::Cluster objects
43  std::string fG4ModuleLabel; ///< label for module running G4 and making particles, etc
44 
45  };
46 }
47 
48 namespace cheat{
49 
50  //--------------------------------------------------------------------
52  {
53  this->reconfigure(pset);
54 
55  produces< std::vector<rb::Prong> >();
56  produces< art::Assns<rb::Prong, rb::CellHit> >();
57  }
58 
59  //--------------------------------------------------------------------
61  {
62  }
63 
64  //--------------------------------------------------------------------
66  {
67  fCheatedClusterLabel = pset.get< std::string >("CheatedClusterLabel", "cluster" );
68  fG4ModuleLabel = pset.get< std::string >("G4ModuleLabel", "geantgen");
69 
70  return;
71  }
72 
73  //--------------------------------------------------------------------
75  {
77  const sim::ParticleNavigator& pnav = bt->ParticleNavigator();
78 
79  // grab the clusters that have been reconstructed
81  evt.getByLabel(fCheatedClusterLabel, clustercol);
82 
83  // make a vector of them - we aren't writing anything out to a file
84  // so no need for a art::PtrVector here
85  std::vector< art::Ptr<rb::Cluster> > clusters;
86  art::fill_ptr_vector(clusters, clustercol);
87 
88  // loop over the clusters and figure out which particle contributed to each one
89  std::vector< art::Ptr<rb::Cluster> >::iterator itr = clusters.begin();
90 
91  // make a map of vectors of art::Ptrs keyed by trackID values
92  std::map< int, std::vector< art::Ptr<rb::Cluster> > > tidClusterMap;
93  std::map< int, std::vector< art::Ptr<rb::Cluster> > >::iterator clusterMapItr = tidClusterMap.begin();
94 
95  // loop over all clusters and fill in the map
96  while( itr != clusters.end() ){
97 
98  // in the ClusterCheater module we set the cluster ID to be
99  // the particle track ID*1000 + view*10. The
100  // floor function on the cluster ID / 1000 will give us
101  // the track ID
102  int tID = floor((*itr)->ID()/1000.);
103 
104  tidClusterMap[tID].push_back((*itr));
105 
106  itr++;
107  }// end loop over clusters
108 
109 
110  // loop over the map and make prongs
111  std::unique_ptr< std::vector<rb::Prong> > prongcol(new std::vector<rb::Prong>);
112  std::unique_ptr< art::Assns<rb::Prong, rb::CellHit> > assn(new art::Assns<rb::Prong, rb::CellHit>);
113 
114  for(clusterMapItr = tidClusterMap.begin(); clusterMapItr != tidClusterMap.end(); clusterMapItr++){
115 
116  // is this prong electro-magnetic in nature or
117  // hadronic/muonic? EM --> shower, everything else is a track
118  if( abs(pnav[(*clusterMapItr).first]->PdgCode()) == 11 ||
119  abs(pnav[(*clusterMapItr).first]->PdgCode()) == 22 ||
120  abs(pnav[(*clusterMapItr).first]->PdgCode()) == 111 ){
121 
122  LOG_DEBUG("ProngCheater") << "prong of " << (*clusterMapItr).first
123  << " is a EM shower with pdg code "
124  << pnav[(*clusterMapItr).first]->PdgCode();
125 
126  // separate out the hits for each particle into the different views
127  std::vector< art::Ptr<rb::Cluster> > tidClusters( (*clusterMapItr).second );
128 
129  // get the hits from the different clusters
131  for(size_t c = 0; c < tidClusters.size(); ++c){
132  for(size_t h = 0; h < tidClusters[c]->NCell(); ++h)
133  ptrvs.push_back(tidClusters[c]->Cell(h));
134  }
135 
136  // get the initial position and direction of the particle
137  TVector3 pos(pnav[(*clusterMapItr).first]->Vx(),
138  pnav[(*clusterMapItr).first]->Vy(),
139  pnav[(*clusterMapItr).first]->Vz());
140 
141  TVector3 dir(pnav[(*clusterMapItr).first]->Px()/pnav[(*clusterMapItr).first]->P(),
142  pnav[(*clusterMapItr).first]->Py()/pnav[(*clusterMapItr).first]->P(),
143  pnav[(*clusterMapItr).first]->Pz()/pnav[(*clusterMapItr).first]->P());
144 
145  // add a prong to the collection. Make the prong
146  // ID be the same as the track ID for the eve particle
147  prongcol->push_back(rb::Prong(ptrvs, pos, dir, (*clusterMapItr).first));
148 
149  // get the hits associated with each cluster and associate those with the track
150  art::FindManyP<rb::CellHit> fmh(tidClusters, evt, fCheatedClusterLabel);
151  for(size_t c = 0; c < tidClusters.size(); ++c){
152  std::vector< art::Ptr<rb::CellHit> > hits = fmh.at(c);
153  util::CreateAssn(*this, evt, *prongcol, hits, *assn);
154  }
155 
156  LOG_DEBUG("ProngCheater") << "adding prong: \n"
157  << prongcol->back()
158  << "\nto collection.";
159 
160  }//end if this is a track
161 
162  } // end loop over the map
163 
164  evt.put(std::move(prongcol));
165  evt.put(std::move(assn));
166 
167  return;
168 
169  } // end produce
170 
171 } // end namespace
172 
173 namespace cheat{
174 
176 
177 }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
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.
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
void abs(TH1 *hist)
void produce(art::Event &evt)
DEFINE_ART_MODULE(TestTMapFile)
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
void hits()
Definition: readHits.C:15
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
T get(std::string const &key) const
Definition: ParameterSet.h:231
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::string fCheatedClusterLabel
label for module creating recob::Cluster objects
A Cluster with defined start position and direction.
Definition: Prong.h:19
code to link reconstructed objects back to the MC truth information
Definition: FillTruth.h:16
void reconfigure(fhicl::ParameterSet const &pset)
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
ProngCheater(fhicl::ParameterSet const &pset)
std::string fG4ModuleLabel
label for module running G4 and making particles, etc
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464