ShowerCheater_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief ShowerCheater module
3 ///
4 /// \author brebel@fnal.gov
5 ////////////////////////////////////////////////////////////////////////
6 
7 #include <string>
8 #include <vector>
9 
10 // NOvASoft includes
11 #include "MCCheater/BackTracker.h"
12 #include "RecoBase/Prong.h"
13 #include "RecoBase/CellHit.h"
14 #include "RecoBase/Cluster.h"
16 #include "Utilities/AssociationUtil.h"
17 
18 // Framework includes
20 #include "fhiclcpp/ParameterSet.h"
29 
30 namespace cheat {
31  class ShowerCheater : public art::EDProducer {
32  public:
33  explicit ShowerCheater(fhicl::ParameterSet const& pset);
34  virtual ~ShowerCheater();
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 track ID 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 eve particle track ID*1000 + view*10. The
100  // floor function on the cluster ID / 1000 will give us
101  // the eve 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  // loop over the map and make prongs
110  std::unique_ptr< std::vector<rb::Prong> > showercol(new std::vector<rb::Prong>);
111  std::unique_ptr< art::Assns<rb::Prong, rb::CellHit> > assn(new art::Assns<rb::Prong, rb::CellHit>);
112 
113  for(clusterMapItr = tidClusterMap.begin(); clusterMapItr != tidClusterMap.end(); clusterMapItr++){
114 
115  // is this prong electro-magnetic in nature or
116  // hadronic/muonic? EM --> shower, everything else is a track
117  if( abs(pnav[(*clusterMapItr).first]->PdgCode()) == 11 ||
118  abs(pnav[(*clusterMapItr).first]->PdgCode()) == 22 ||
119  abs(pnav[(*clusterMapItr).first]->PdgCode()) == 111 ){
120 
121  LOG_DEBUG("ShowerCheater") << "prong of " << (*clusterMapItr).first
122  << " is a shower with pdg code "
123  << pnav[(*clusterMapItr).first]->PdgCode();
124 
125  // separate out the hits for each particle into the different views
126  std::vector< art::Ptr<rb::Cluster> > tidClusters( (*clusterMapItr).second );
127 
129  for(size_t c = 0; c < tidClusters.size(); ++c){
130  for(size_t h = 0; h < tidClusters[c]->NCell(); ++h)
131  ptrvs.push_back(tidClusters[c]->Cell(h));
132  }
133 
134  // get the initial position and direction of the particle
135  TVector3 pos(pnav[(*clusterMapItr).first]->Vx(),
136  pnav[(*clusterMapItr).first]->Vy(),
137  pnav[(*clusterMapItr).first]->Vz());
138 
139  TVector3 dir(pnav[(*clusterMapItr).first]->Px()/pnav[(*clusterMapItr).first]->P(),
140  pnav[(*clusterMapItr).first]->Py()/pnav[(*clusterMapItr).first]->P(),
141  pnav[(*clusterMapItr).first]->Pz()/pnav[(*clusterMapItr).first]->P());
142 
143  // add a prong to the collection. Make the prong
144  // ID be the same as the track ID for the eve particle
145  showercol->push_back(rb::Prong(ptrvs, pos, dir, (*clusterMapItr).first));
146 
147  // get the hits associated with each cluster and associate those with the track
148  art::FindManyP<rb::CellHit> fmh(tidClusters, evt, fCheatedClusterLabel);
149  for(size_t c = 0; c < tidClusters.size(); ++c){
150  std::vector< art::Ptr<rb::CellHit> > hits = fmh.at(c);
151  util::CreateAssn(*this, evt, *showercol, hits, *assn);
152  }
153 
154  LOG_DEBUG("ShowerCheater") << "adding shower: \n"
155  << showercol->back()
156  << "\nto collection.";
157 
158  }// end if this is a shower
159  } // end loop over the map
160 
161  evt.put(std::move(showercol));
162  evt.put(std::move(assn));
163 
164  return;
165 
166  } // end produce
167 
168 } // end namespace
169 
170 
171 namespace cheat{
172 
174 
175 }
#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
std::string fCheatedClusterLabel
label for module creating recob::Cluster objects
ShowerCheater(fhicl::ParameterSet const &pset)
void abs(TH1 *hist)
DEFINE_ART_MODULE(TestTMapFile)
std::string fG4ModuleLabel
label for module running G4 and making particles, etc
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
void produce(art::Event &evt)
void reconfigure(fhicl::ParameterSet const &pset)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
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
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
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464