RecoCheckAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief
3 /// \author brebel@fnal.gov
4 /// \date Jul. 15 2011
5 ////////////////////////////////////////////////////////////////////////
6 
7 #include "TH1.h"
8 #include "TH2.h"
9 
15 
16 #include "MCCheater/BackTracker.h"
17 #include "RecoBase/Prong.h"
18 #include "RecoBase/Track.h"
19 #include "RecoBase/Cluster.h"
20 #include "RecoBase/CellHit.h"
21 #include "Simulation/Particle.h"
22 #include "Simulation/Simulation.h"
23 #include "Geometry/Geometry.h"
25 
26 namespace cheat {
27  class RecoCheckAna;
28 }
29 
30 class TH1D;
31 class TH2D;
32 
34 public:
35  explicit RecoCheckAna(fhicl::ParameterSet const &p);
36  virtual ~RecoCheckAna();
37 
38  virtual void analyze(art::Event const &e);
39 
40  virtual void reconfigure(fhicl::ParameterSet const & p);
41  virtual void beginJob();
42 
43 
44 
45 private:
46  struct Plots
47  {
50  effvpur(0) {}
51  std::set<int> pdgs;
52  std::set<int> allowedDaughters; ///< To not count eg Michel electrons against purity
55  TH2D* effvpur;
56  };
57 
58  void CheckRecoClusters(std::vector< art::Ptr<rb::Cluster> > const& clscol,
59  art::PtrVector<rb::CellHit> const& allhits,
60  Plots& plots);
61 
62  /// Total length in either 2D or 3D
64  geo::View_t view) const;
65 
66  std::set<int> FindVisibleProngs(const art::PtrVector<rb::CellHit>& allhits,
67  std::set<int> ids) const;
68 
69  std::string fClusterAssnLabel; ///< label for module clusters, showers, and tracks are associated with
70  std::string fShowerModuleLabel; ///< label for module making the showers
71  std::string fTrackModuleLabel; ///< label for module making the tracks
72 
73  bool fCheckShowers; ///< should we check the reconstruction of showers?
74  bool fCheckTracks; ///< should we check the reconstruction of tracks?
75  bool fLinkByMaxEff; ///< If set to true, finds mc particle with max efficiency
76  ///< and plots the corresponding efficiency and purity. If
77  ///< set to false, finds mc particle with the max purity and
78  ///< plots the corresponding efficiency and purity.
79 
80  std::vector<Plots> fPlots;
83 };
84 
85 
86 //-------------------------------------------------------------------
88  : EDAnalyzer(p)
89 {
90  this->reconfigure(p);
91 }
92 
93 //-------------------------------------------------------------------
95 {
96 }
97 
98 //-------------------------------------------------------------------
99 // We want to determine how well an event finding algorithm works by
100 // determining the following:
101 // 1. How much of the actual energy for a given interaction is
102 // found in the reconstructed event
103 // 2. How many of the particles from a given interaction are assigned
104 // to the correct reconstructed event
105 // 3. How close is the reconstructed event primary vertex to the
106 // simulated interaction primary vertex
107 // 4. What fraction of energy in the reconstructed event is from
108 // particles not associated with the simulated interaction
109 
111 {
112  // check that this is MC, stop if it isn't
113  if(e.isRealData()){
114  mf::LogWarning("RecoCheckAna") << "attempting to run MC truth check on "
115  << "real data, bail";
116  return;
117  }
118 
119 
120  // get a handle on the clusters that the clusters/tracks/showers are associated with (slices)
122  e.getByLabel(fClusterAssnLabel, hslices);
123  if(hslices.failedToGet()){
124  mf::LogWarning("RecoCheckAna") << "could not recover slices";
125  }
126 
127 
128  // get all the association finder objects
129  art::FindManyP<rb::Track> fmtrack(hslices,e,fTrackModuleLabel);
131 
132  // loop over the slices
133  for(unsigned int islc = 0; islc < hslices->size(); ++islc){
134  art::PtrVector<rb::CellHit> allhits = hslices->at(islc).AllCells();
135  if(fCheckTracks){
136  if(fmtrack.isValid()){
137  std::vector< art::Ptr<rb::Track> > trkcol = fmtrack.at(islc);
138  std::vector< art::Ptr<rb::Cluster> > tracks;
139  for(unsigned int itrk = 0; itrk < trkcol.size(); ++itrk){
140  tracks.push_back(trkcol[itrk]);
141  }
142  for(unsigned int n = 0; n < fPlots.size(); ++n){
143  CheckRecoClusters(tracks, allhits, fPlots[n]);
144  }
145  }
146  else{
147  mf::LogWarning("RecoCheckAna") << "could not recover tracks";
148  }
149  }
150 
151  if(fCheckShowers){
152  if(fmprong.isValid()){
153  std::vector< art::Ptr<rb::Prong> > shwcol = fmprong.at(islc);
154  std::vector< art::Ptr<rb::Cluster> > shws;
155  for(unsigned int ishw = 0; ishw < shwcol.size(); ++ishw){
156  shws.push_back(shwcol[ishw]);
157  }
158  this->CheckRecoClusters(shws, allhits, fShowerPlots);
159  }
160  else{
161  mf::LogWarning("RecoCheckAna") << "could not recover showers";
162  }
163  }
164 
165  }
166 
167 }
168 
169 //-------------------------------------------------------------------
171 {
172  fClusterAssnLabel = p.get< std::string >("ClusterAssnLabel" );
173  fShowerModuleLabel = p.get< std::string >("ShowerModuleLabel" );
174  fTrackModuleLabel = p.get< std::string >("TrackModuleLabel" );
175 
176  fCheckShowers = p.get< bool >("CheckShowers" );
177  fCheckTracks = p.get< bool >("CheckTracks" );
178  fLinkByMaxEff = p.get< bool >("LinkByMaxEff" );
179 }
180 
181 //-------------------------------------------------------------------
183 {
185 
186  if(fCheckTracks){
187  for(int n = 0; n < 4; ++n){
189  std::set<int> pdgs;
190  std::set<int> allowedDaughters;
191 
192  switch(n){
193  case 0:
194  name = "all";
195  break;
196  case 1:
197  name = "muon";
198  pdgs.insert(-13); pdgs.insert(+13);
199  allowedDaughters.insert(-11); allowedDaughters.insert(+11);
200  break;
201  case 2:
202  name = "proton";
203  pdgs.insert(-2212); pdgs.insert(+2212);
204  break;
205  case 3:
206  name = "piPlusMinus";
207  pdgs.insert(-211); pdgs.insert(+211);
208  break;
209  }
210 
211  Plots p;
212  p.purity = tfs->make<TH1D>((name+"Purity").c_str(), ";Purity;Tracks", 110, 0, 1.1);
213  p.efficiency = tfs->make<TH1D>((name+"Efficiency").c_str(), ";Efficiency;Tracks", 110, 0, 1.1);
214  p.deltaLength = tfs->make<TH1D>((name+"DeltaLength").c_str(), ";#DeltaL (cm);Tracks", 100, -50, +50);
215  p.purityVsLength = tfs->make<TH2D>((name+"PurityVsLength").c_str(), ";True length (cm);Purity;Tracks", 300, 0, 1500, 110, 0, 1.1);
216  p.efficiencyVsLength = tfs->make<TH2D>((name+"EfficiencyVsLength").c_str(), ";True length (cm);Efficiency;Tracks", 300, 0, 1500, 110, 0, 1.1);
217  p.deltaLengthVsLength = tfs->make<TH2D>((name+"DeltaLengthVsLength").c_str(), ";True length (cm);#DeltaL (cm);Tracks", 300, 0, 1500, 100, -50, +50);
218  p.effvpur = tfs->make<TH2D>((name+"EffvPurity").c_str(), ";Efficiency;Purity", 110, 0, 1.1, 110, 0, 1.1);
219  p.pdgs = pdgs;
220  p.allowedDaughters = allowedDaughters;
221 
222  fPlots.push_back(p);
223  }
224  }
225 
226  if(fCheckShowers){
227  fShowerPlots.purity = tfs->make<TH1D>("showerPurity", ";Purity;Showers", 110, 0., 1.1);
228  fShowerPlots.efficiency = tfs->make<TH1D>("showerEfficiency", ";Efficiency;Showers", 110, 0., 1.1);
229  }
230 
231 }
232 
233 
234 //-------------------------------------------------------------------
236  art::PtrVector<rb::CellHit> const& allhits,
237  Plots& plots)
238 {
239  // grab the set of eve IDs for this plist
241  std::set<int> trackIDs = bt->GetTrackIDList();
242 
243  // Only keep the ones associated with not completely unreasonable tracks
244  trackIDs = FindVisibleProngs(allhits, trackIDs);
245 
246  // Empty set means accept all PDGs. Otherwise, filter out the unwanted ones
247  if(!plots.pdgs.empty()){
248  std::set<int> newTrackIDs;
249  for(std::set<int>::iterator it = trackIDs.begin(); it != trackIDs.end(); ++it){
250  const int pdg = bt->TrackIDToParticle(*it)->PdgCode();
251  if(plots.pdgs.find(pdg) != plots.pdgs.end()) newTrackIDs.insert(*it);
252  } // end for it
253  trackIDs = newTrackIDs;
254  }
255 
256  std::map<int, int> parents;
257  if(!plots.allowedDaughters.empty()){
258  for(std::set<int>::iterator it = trackIDs.begin(); it != trackIDs.end(); ++it){
259  const sim::Particle* mother = bt->TrackIDToParticle(*it);
260  if(!mother) continue; // Can this happen?
261 
262  for(int i = 0; i < mother->NumberDaughters(); ++i){
263  const int id = mother->Daughter(i);
264  assert(bt->TrackIDToParticle(id));
265  if(plots.allowedDaughters.count(bt->TrackIDToParticle(id)->PdgCode())){
266  parents[id] = *it;
267  }
268  } // end for i
269  } // end for it
270  } // end if
271 
272  // Indexed by what view they're in. X, Y, or 3D
273  assert(geo::kXorY == 2); // Ensure enum still works how we think
274  std::map<int, double> bestPur[3], bestEff[3], bestLength[3];
275 
276  // Go through each reconstructed object, accumulating the best scores seen
277  // for each trackID we're interested in.
278  for(size_t c = 0; c < clscol.size(); ++c){
279  const rb::Cluster* cls = clscol[c].get();
280  const geo::View_t view = cls->View();
281 
282  // get the hits associated with this object
284 
285  // use the cheat::BackTracker to find purity and efficiency for these
286  // hits for each particle ID
287  std::map<int, double> purMap, effMap;
288  bt->HitCollectionPurity(trackIDs, hits, &purMap, &parents);
289  bt->HitCollectionEfficiency(trackIDs, hits, allhits, view, &effMap);
290 
291  for(std::set<int>::iterator it = trackIDs.begin(); it != trackIDs.end(); ++it){
292  const double p = purMap[*it];
293  const double e = effMap[*it];
294 
295  if(( fLinkByMaxEff && e > bestEff[view][*it]) ||
296  (!fLinkByMaxEff && p > bestPur[view][*it])) {
297  bestEff[view][*it] = e;
298  bestPur[view][*it] = p;
299 
300  // If we were told to fill a track length histogram then downcast
301  // should be safe.
302  if(plots.deltaLength) bestLength[view][*it] = ((rb::Track*)cls)->TotalLength();
303  }
304  } // end for it
305  } // end for c
306 
307  // For every track. If it was matched in 3D, fill the histograms.
308  // If there is no 3D match look for 2 2D matches and fill with weight 0.5
309  // If required we can add a flag to only try 3D or only try 2D.
310  for(std::set<int>::iterator it = trackIDs.begin(); it != trackIDs.end(); ++it){
311  // Try 3D first
312  for(int view = geo::kXorY; view >= geo::kX; --view){
313  // If no 3D match, continue to two 2D options
314  if(view == geo::kXorY && bestEff[view].find(*it) == bestEff[view].end()) continue;
315 
316  const double weight = (view == geo::kXorY) ? 1 : 0.5;
317 
318  if(plots.purity) plots.purity->Fill(bestPur[view][*it], weight);
319  if(plots.efficiency) plots.efficiency->Fill(bestEff[view][*it], weight);
320  if(plots.effvpur) plots.effvpur->Fill(bestEff[view][*it], bestPur[view][*it], weight);
321 
322  const double trueLength = TotalLengthInDetector(bt->TrackIDToParticle(*it), geo::View_t(view));
323 
324  if(plots.purityVsLength) plots.purityVsLength->Fill(trueLength, bestPur[view][*it], weight);
325  if(plots.efficiencyVsLength) plots.efficiencyVsLength->Fill(trueLength, bestEff[view][*it], weight);
326 
327  if(plots.deltaLength) plots.deltaLength->Fill(bestLength[view][*it]-trueLength, weight);
328  if(plots.deltaLengthVsLength) plots.deltaLengthVsLength->Fill(trueLength, bestLength[view][*it]-trueLength, weight);
329 
330  // If we did find a 3D match, don't try 2D
331  if(bestEff[geo::kXorY].find(*it) != bestEff[geo::kXorY].end()) break;
332  } // end for view
333  } // end for it
334 }
335 
336 //-------------------------------------------------------------------
338  geo::View_t view) const
339 {
340  const unsigned int N = p->NumberTrajectoryPoints();
341 
342  if(N < 2){
343  mf::LogWarning("RecoCheckAna") << "Particle has no trajectory points";
344  return 0;
345  }
346 
348 
349  double dist = 0;
350 
351  for(unsigned int n = 0; n < N-1; ++n){
352  TVector3 p0 = p->Position(n).Vect();
353  TVector3 p1 = p->Position(n+1).Vect();
354  if(p0.Z() < 0 || p1.Z() < 0 ||
355  p0.Z() > geom->DetLength() || p1.Z() > geom->DetLength() ||
356  fabs(p0.X()) > geom->DetHalfWidth() || fabs(p1.X()) > geom->DetHalfWidth() ||
357  fabs(p0.Y()) > geom->DetHalfHeight() || fabs(p1.Y()) > geom->DetHalfHeight()) continue;
358 
359  if(view == geo::kX){
360  p0.SetY(0);
361  p1.SetY(0);
362  }
363  if(view == geo::kY){
364  p0.SetX(0);
365  p1.SetX(0);
366  }
367  dist += (p1-p0).Mag();
368  }
369 
370  return dist;
371 }
372 
373 struct Counts
374 {
376  {
377  nHits[0] = nHits[1] = 0;
378  nOnlyHits[0] = nOnlyHits[1] = 0;
379  }
380 
381  int nHits[2]; // Indexed by view
382  int nOnlyHits[2]; // Hit by only one particle
383 };
384 
385 //-------------------------------------------------------------------
386 std::set<int> cheat::RecoCheckAna::
388  std::set<int> ids) const
389 {
390  std::map<int, Counts> counts;
391 
393 
394  // loop over all hits and check if this track id contributed to them
395  for(unsigned int i = 0; i < allhits.size(); ++i){
396  //get vector of track id's that contribute to this hit
397  const art::Ptr<rb::CellHit> hit = allhits[i];
398  const std::vector<TrackIDE> trackIDs = bt->HitToTrackIDE(hit);
399  // loop over all the returned track ids
400  int nTracks = 0; // number of tracks that contributed energy to the hit
401  int loneID = -1; // Which ID it was if there was only one
402  for(unsigned int j = 0; j < trackIDs.size(); ++j){
403  const int trackID = trackIDs[j].trackID;
404  if(trackID != sim::kNoiseId){
405  ++nTracks;
406  loneID = trackID;
407  }
408 
409  // check if the current id matches the input track id
410  if(ids.find(trackID) != ids.end()){
411  ++counts[trackID].nHits[hit->View()];
412  }
413  } // end for j
414  if(nTracks == 1){
415  ++counts[loneID].nOnlyHits[hit->View()];
416  }
417  } // end for i
418 
419  const int kMinHits = 2; // minimum number of hits to be considered visible
420  const int kMinOnlyHit = 1; // minimum number of hits in which it is only particle contributing energy to be visible
421 
422  std::set<int> ret;
423  for(std::set<int>::iterator it = ids.begin(); it != ids.end(); ++it){
424  const Counts c = counts[*it];
425  // Check if this track passes the visible prong requirements
426  if(c.nOnlyHits[geo::kX] >= kMinOnlyHit &&
427  c.nOnlyHits[geo::kY] >= kMinOnlyHit){
428 
429  if(c.nHits[geo::kX] >= kMinHits && c.nHits[geo::kY] >= 1)
430  ret.insert(*it);
431 
432  if(c.nHits[geo::kY] >= kMinHits && c.nHits[geo::kX] >= 1)
433  ret.insert(*it);
434  }
435  }
436  return ret;
437 }
438 
::xsd::cxx::tree::id< char, ncname > id
Definition: Database.h:165
virtual geo::View_t View() const
kXorY for 3D clusters.
Definition: Cluster.h:99
const XML_Char * name
Definition: expat.h:151
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:217
back track the reconstruction to the simulation
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
int PdgCode() const
Definition: MCParticle.h:211
std::vector< Plots > fPlots
set< int >::iterator it
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
std::vector< TrackIDE > HitToTrackIDE(const rb::CellHit &hit, bool useBirksE=false) const
Convenience function. HitsToTrackIDE but for a single hit.
static const int kNoiseId
flag for noise id
Definition: Simulation.h:13
const Var weight
X or Y views.
Definition: PlaneGeo.h:30
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
geo::View_t View() const
Definition: CellHit.h:41
const char * p
Definition: xmltok.h:285
Vertical planes which measure X.
Definition: PlaneGeo.h:28
A collection of associated CellHits.
Definition: Cluster.h:47
double DetLength() const
RecoCheckAna(fhicl::ParameterSet const &p)
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
std::set< int > allowedDaughters
To not count eg Michel electrons against purity.
bool isRealData() const
Definition: Event.h:83
DEFINE_ART_MODULE(TestTMapFile)
int NumberDaughters() const
Definition: MCParticle.h:216
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
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...
int Daughter(const int i) const
Definition: MCParticle.cxx:112
double dist
Definition: runWimpSim.h:113
virtual void reconfigure(fhicl::ParameterSet const &p)
TString part[npart]
Definition: Style.C:32
void hits()
Definition: readHits.C:15
std::set< int > FindVisibleProngs(const art::PtrVector< rb::CellHit > &allhits, std::set< int > ids) const
const std::vector< Plot > plots
T get(std::string const &key) const
Definition: ParameterSet.h:231
virtual void analyze(art::Event const &e)
const double j
Definition: BetheBloch.cxx:29
bool fCheckTracks
should we check the reconstruction of tracks?
std::string fTrackModuleLabel
label for module making the tracks
double DetHalfHeight() const
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
size_type size() const
Definition: PtrVector.h:308
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...
bool fCheckShowers
should we check the reconstruction of showers?
const sim::Particle * TrackIDToParticle(int const &id) const
Returns a pointer to the sim::Particle object corresponding to the given TrackID. ...
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
code to link reconstructed objects back to the MC truth information
Definition: FillTruth.h:16
double DetHalfWidth() const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
std::string fClusterAssnLabel
label for module clusters, showers, and tracks are associated with
Definition: structs.h:12
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void geom(int which=0)
Definition: geom.C:163
float Mag() const
const std::set< int > GetTrackIDList() const
Get all G4 track ids present in the event.
Definition: BackTracker.h:750
assert(nhit_max >=nhit_nbins)
std::string fShowerModuleLabel
label for module making the showers
Int_t trackID
Definition: plot.C:84
double TotalLengthInDetector(const sim::Particle *part, geo::View_t view) const
Total length in either 2D or 3D.
Float_t e
Definition: plot.C:35
void CheckRecoClusters(std::vector< art::Ptr< rb::Cluster > > const &clscol, art::PtrVector< rb::CellHit > const &allhits, Plots &plots)
Encapsulate the geometry of one entire detector (near, far, ndos)
bool failedToGet() const
Definition: Handle.h:196