CheckBackTracking_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief CheckBackTracking module
3 ///
4 /// \author brebel@fnal.gov
5 ///
6 ////////////////////////////////////////////////////////////////////////
7 
8 #include <string>
9 #include <vector>
10 
11 // Framework includes
16 #include "fhiclcpp/ParameterSet.h"
21 
22 // NOvASoft includes
23 #include "Geometry/Geometry.h"
24 #include "MCCheater/BackTracker.h"
25 #include "RecoBase/CellHit.h"
26 #include "Simulation/Particle.h"
29 #include "Simulation/FLSHit.h"
30 
31 namespace cheat {
33  public:
34  explicit CheckBackTracking(fhicl::ParameterSet const& pset);
35  virtual ~CheckBackTracking();
36 
37  void analyze(art::Event const& evt);
38 
39  void reconfigure(fhicl::ParameterSet const& pset);
40 
41  private:
42 
46  void CheckPurityAndEfficiency(std::vector< art::Ptr<rb::CellHit> > const& hits);
47 
49  std::string fHitModuleLabel; ///< label for module creating rb::CellHit objects
50 
51  };
52 }
53 
54 namespace cheat{
55 
56  //--------------------------------------------------------------------
58  : EDAnalyzer(pset)
59  {
60  this->reconfigure(pset);
61  }
62 
63  //--------------------------------------------------------------------
65  {
66  }
67 
68  //--------------------------------------------------------------------
70  {
71  fHitModuleLabel = pset.get<std::string>("HitModuleLabel", "calhit");
72  }
73 
74  //--------------------------------------------------------------------
76  {
77 
78  // grab the hits that have been reconstructed
80  evt.getByLabel(fHitModuleLabel, hitcol);
81 
82  // make a vector of them - we aren't writing anything out to a file
83  // so no need for a art::PtrVector here
84  std::vector< art::Ptr<rb::CellHit> > hits;
85  art::fill_ptr_vector(hits, hitcol);
86 
87  // loop over the hits and figure out which particle contributed to each one
88  std::vector< art::Ptr<rb::CellHit> >::iterator itr = hits.begin();
89 
90  while( itr != hits.end() ){
91 
92  // print the information for this hit
93  mf::LogVerbatim("CheckBackTracking") << *(*itr).get();
94 
95  // check is this is a noise hit
96  if (fBT->IsNoise((*itr))){
97  mf::LogVerbatim("CheckBackTracking") << "This hit is a noise hit."
98  << " No truth info will be printed.";
99  itr++;
100  continue;
101  }
102 
103  this->CheckFLSHitConsistency(*itr);
104 
105  TVector3 xyz = fBT->HitToXYZ((*itr));
106  mf::LogVerbatim("CheckBackTracking") << "hit weighted mean position is ("
107  << xyz[0] << "," << xyz[1] << "," << xyz[2] << ")";
108 
109  std::vector<cheat::TrackIDE> trackides = fBT->HitToTrackIDE((*itr));
110  for(size_t t = 0; t < trackides.size(); ++t){
111  mf::LogVerbatim("CheckBackTracking") << "track id: " << trackides[t].trackID
112  << " contributed " << trackides[t].energy << "/"
113  << trackides[t].energyFrac
114  << " to the current hit";
115  }
116 
117  const sim::Particle* part = fBT->HitToParticle(*itr);
118  if(part){
119  mf::LogVerbatim("CheckBackTracking") << "largest contributing particle is: " << *part;
120  }
121  else{
122  mf::LogVerbatim("CheckBackTracking") << "no particle corresponding to CellHit: " << *itr;
123  }
124 
125  itr++;
126  }// end loop over hits
127 
128  this->CheckTrackIdToParticle();
129  this->CheckPurityAndEfficiency(hits);
130 
131  return;
132 
133  } // end analyze
134 
135  //----------------------------------------------------------------------------
137  {
138  unsigned int plane = hit->Plane();
139  unsigned int cell = hit->Cell();
140 
141  const std::vector<sim::FLSHit> flshits = fBT->CellToFLSHit(plane, cell);
142  const std::vector<sim::FLSHit> flshits2 = fBT->HitToFLSHit(hit);
143 
144  for (size_t f = 0; f < flshits2.size(); ++f) {
145  if ( std::find(flshits.begin(), flshits.end(), flshits2[f]) == flshits.end() ) {
146  throw cet::exception("CheckBackTracking")
147  << "FLSHits from Cell do not match FLSHits from Hit\n"
148  << __FILE__ << ":" << __LINE__ << "\n";
149  }
150  }
151 
152  return;
153  }
154 
155  //----------------------------------------------------------------------------
157  {
158  // loop over the eveID values and calculate the purity and efficiency for each
159  std::set<int>::iterator setitr = fBT->GetTrackIDList().begin();
160  while( setitr != fBT->GetTrackIDList().end() ){
161 
162  sim::ParticleHistory phist(&(fBT->ParticleNavigator()), *setitr);
163 
164  const sim::Particle* pfromID = fBT->TrackIDToParticle(*setitr);
165  const sim::Particle* pfromNav = fBT->ParticleNavigator().find(*setitr)->second;
166  if ( pfromID != pfromNav ) {
167  throw cet::exception("CheckBackTracking")
168  << "Same track ID returns different pointers: "
169  << *pfromID << "\n"
170  << *pfromNav << "\n"
171  << __FILE__ << ":" << __LINE__ << "\n";
172  }
173 
174  // attempt to recover the MCTruth object for this particle
176  mf::LogVerbatim("CheckBackTracking") << "MCTruth for " << *setitr << " is "
177  << *(mct.get());
178 
179  // check that you get the same art::Ptr<simb::MCTruth> from both accessors
180  if ( mct.id() != fBT->ParticleToMCTruth(pfromID).id() ) {
181  throw cet::exception("CheckBackTrack")
182  << "Track ID and Particle return different MCTruth:\n"
183  << *(mct.get()) << "\n"
184  << *(fBT->ParticleToMCTruth(pfromID).get()) << "\n"
185  << __FILE__ << ":" << __LINE__ << "\n";
186  }
187 
188  // get the list of particles for this MCTruth
189  std::vector<const sim::Particle*> parts = fBT->MCTruthToParticles(mct);
190  if ( std::find(parts.begin(), parts.end(), pfromID) == parts.end() ) {
191  throw cet::exception("CheckBackTracking")
192  << "Cannot map particle to MCTruth and back again\n"
193  << __FILE__ << ":" << __LINE__ << "\n";
194  }
195 
196  // check the mother particle
197  pfromNav = fBT->ParticleNavigator().find(pfromID->Mother())->second;
198  pfromID = fBT->TrackIDToMotherParticle(*setitr);
199  if ( pfromID != pfromNav ) {
200  throw cet::exception("CheckBackTracking")
201  << "Same track ID returns different mothers: "
202  << *pfromID << "\n"
203  << *pfromNav << "\n\t Parentage is:\n\t"
204  << phist << "\n"
205  << __FILE__ << ":" << __LINE__ << "\n";
206  }
207 
208  }
209 
210  return;
211  }
212 
213  //----------------------------------------------------------------------------
215  {
216  art::PtrVector<rb::CellHit> hitPtrVec;
217  for(size_t h = 0; h < hits.size(); ++h) hitPtrVec.push_back(hits[h]);
218 
219  // loop over the track ID values and calculate the purity and efficiency for each
220  std::set<int>::iterator setitr = fBT->GetTrackIDList().begin();
221  while( setitr != fBT->GetTrackIDList().end() ){
222 
223  std::set<int> id;
224  id.insert(*setitr);
225  mf::LogVerbatim ("CheckBackTracking") << "track ID: " << *setitr
226  << " purity: "
227  << fBT->HitCollectionPurity(id, hitPtrVec)
228  << " efficiency: "
229  << fBT->HitCollectionEfficiency(id, hitPtrVec, hitPtrVec, geo::kXorY);
230 
231 
232  setitr++;
233  }// end loop over track ids
234 
235  return;
236  }
237 
238 } // end namespace
239 
240 
241 namespace cheat{
242 
244 
245 }
::xsd::cxx::tree::id< char, ncname > id
Definition: Database.h:165
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
back track the reconstruction to the simulation
std::string fHitModuleLabel
label for module creating rb::CellHit objects
const sim::Particle * TrackIDToMotherParticle(int const &id) const
std::vector< TrackIDE > HitToTrackIDE(const rb::CellHit &hit, bool useBirksE=false) const
Convenience function. HitsToTrackIDE but for a single hit.
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
X or Y views.
Definition: PlaneGeo.h:30
unsigned short Plane() const
Definition: CellHit.h:39
void reconfigure(fhicl::ParameterSet const &pset)
int Mother() const
Definition: MCParticle.h:212
CheckBackTracking(fhicl::ParameterSet const &pset)
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
DEFINE_ART_MODULE(TestTMapFile)
const sim::Particle * HitToParticle(art::Ptr< rb::CellHit > const &hit, bool quiet=false) const
Returns the sim::Particle object contributing the most light to the specified rb::CellHit.
bool IsNoise(const art::Ptr< rb::CellHit > &hit) const
Is this hit not associated with any particles?
double HitCollectionEfficiency(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, const std::vector< const rb::CellHit * > &allhits, const geo::View_t &view, std::map< int, double > *effMap=0, bool energyEff=false, double *desiredEnergy=0, double *totalEnergy=0, int *desiredHits=0, int *totalHits=0) const
Returns the fraction of all energy in an event from a specific set of Geant4 track IDs that are repre...
unsigned short Cell() const
Definition: CellHit.h:40
TString part[npart]
Definition: Style.C:32
void hits()
Definition: readHits.C:15
if(dump)
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
std::vector< sim::FLSHit > HitToFLSHit(const rb::CellHit &hit) const
All the FLSHits that contributed to this hit, sorted from most to least light.
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
const art::Ptr< simb::MCTruth > & ParticleToMCTruth(const sim::Particle *p) const
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
void CheckPurityAndEfficiency(std::vector< art::Ptr< rb::CellHit > > const &hits)
ProductID id() const
Definition: Ptr.h:349
double HitCollectionPurity(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, std::map< int, double > *purMap=0, std::map< int, int > *parents=0, bool energyPur=false) const
Returns the fraction of hits in a collection that come from the specified Geant4 track ids...
const sim::Particle * TrackIDToParticle(int const &id) const
Returns a pointer to the sim::Particle object corresponding to the given TrackID. ...
art::ServiceHandle< cheat::BackTracker > fBT
BackTracker handle.
code to link reconstructed objects back to the MC truth information
Definition: FillTruth.h:17
const art::Ptr< simb::MCTruth > & TrackIDToMCTruth(int const &id) const
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
void CheckFLSHitConsistency(art::Ptr< rb::CellHit > const &hit)
Definition: structs.h:12
const std::set< int > GetTrackIDList() const
Get all G4 track ids present in the event.
Definition: BackTracker.h:750
void CheckParticleConsistency(art::Ptr< rb::CellHit > const &hit)
const std::vector< sim::FLSHit > CellToFLSHit(unsigned int const &plane, unsigned int const &cell) const
Returns the FLSHits contributing the signal in the specified cell. WARNING: Use with extreme caution!...
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
TVector3 HitToXYZ(art::Ptr< rb::CellHit > const &hit, bool useBirksE=false) const
Returns the XYZ position of the energy deposition for a given hit.
void analyze(art::Event const &evt)
iterator find(const key_type &key)
Encapsulate the geometry of one entire detector (near, far, ndos)
std::vector< const sim::Particle * > MCTruthToParticles(art::Ptr< simb::MCTruth > const &mct) const
enum BeamMode string