LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
CosmicRemovalAna_module.cc
Go to the documentation of this file.
1 // Cosmic Removal Module Ana
3 //
4 // Module Designed to loop over tracks / clusters / hits that
5 // have the cosmic tag association and removeor ignore the
6 // hits and check to see what is left over
7 //
8 // Yale Workshop (Cosmics Removal Group)
9 //
11 
12 // Framework includes
18 #include "art_root_io/TFileService.h"
19 #include "fhiclcpp/ParameterSet.h"
20 
21 // LArSoft Includes
31 
32 // ROOT Includes
33 #include "TTree.h"
34 
35 #include <string>
36 #include <vector>
37 
38 namespace microboone {
39 
41 
42  public:
43  explicit CosmicRemovalAna(fhicl::ParameterSet const& pset);
44 
46  void analyze(const art::Event& evt);
47  void beginJob();
48 
50 
51  private:
52  unsigned int nCosmicTags;
53 
54  // TTree *tTagTree;
55  TTree* tEventTree;
56 
57  std::string fHitsModuleLabel;
58  std::string fMCModuleLabel;
59  std::string fMCHitsModuleLabel;
60  std::string fClusterModuleLabel;
61  std::string fTrackModuleLabel;
63  std::vector<std::string> fCosmicTagAssocLabel;
64  std::vector<float> fCosmicScoreThresholds;
65 
66  void InitEventTree(int run_number, int event_number);
67 
68  void FillMCInfo(art::Event const& e,
69  std::vector<recob::Hit> const& hitlist,
70  std::vector<hit_origin_t>& hitOrigins,
71  std::vector<sim::MCHitCollection> const& mchitCollectionVector,
72  std::map<int, const simb::MCTruth*> const& trackIDToTruthMap);
73 
74  void FillTrackInfo(size_t const& hit_iter,
75  hit_origin_t const& origin,
76  float const& charge,
77  std::vector<size_t> const& track_indices_this_hit,
78  std::vector<std::vector<const anab::CosmicTag*>> const& tags_per_cluster,
79  std::vector<bool>& hitsAccounted_per_tag,
80  std::vector<bool>& hitsAllTags);
81 
82  void FillClusterInfo(size_t const& hit_iter,
83  hit_origin_t const& origin,
84  float const& charge,
85  std::vector<size_t> const& cluster_indices_this_hit,
86  std::vector<std::vector<const anab::CosmicTag*>> const& tags_per_cluster,
87  std::vector<bool>& hitsAccounted_per_tag,
88  std::vector<bool>& hitsAllTags);
89 
90  void FillAllTagsInfo(recob::Hit const& hit, hit_origin_t const& origin);
91 
92  std::vector<float> cTaggedCharge_Cosmic;
93  std::vector<float> cTaggedCharge_NonCosmic;
94  std::vector<int> cTaggedHits_Cosmic;
95  std::vector<int> cTaggedHits_NonCosmic;
96  // std::vector<std::string> *cTagAlgorithmNames;
97 
98  /*
99  typedef struct {
100  int eventNumber;
101  int tagType;
102  float x0;
103  float x1;
104  float y0;
105  float y1;
106  float z0;
107  float z1;
108  int nHits;
109  int nGoodHits;
110  int origin;
111  float score;
112  int pdg;
113  float energy;
114  } cTagProperties_t;
115  cTagProperties_t cTagVals;
116  */
117  }; //<---End class
118 }
119 
120 typedef struct {
123 
127 
131 
135 
136  float qTrack;
139 
143 
144  float qCluster;
147 
152 
155 
156 // =====================================================
157 // fhicl::ParameterSet
158 // =====================================================
160  : EDAnalyzer(pset)
161  , fHitsModuleLabel(pset.get<std::string>("HitsModuleLabel"))
162  , fMCModuleLabel(pset.get<std::string>("MCModuleLabel"))
163  , fMCHitsModuleLabel(pset.get<std::string>("MCHitsModuleLabel"))
164  , fClusterModuleLabel(pset.get<std::string>("ClusterModuleLabel"))
165  , fTrackModuleLabel(pset.get<std::string>("TrackModuleLabel"))
166  , fHitCompareCut(pset.get<float>("HitCompareCut"))
167  , fCosmicTagAssocLabel(pset.get<std::vector<std::string>>("CosmicTagAssocLabel"))
168  , fCosmicScoreThresholds(pset.get<std::vector<float>>("CosmicScoreThresholds"))
169 {}
170 
171 // =====================================================
172 // BeginJob
173 // =====================================================
174 
176 {
177 
178  //static cEventProperties_t cEventVals;
179 
185 
187  tEventTree = (TTree*)tfs->make<TTree>("CosmicEventTree", "CosmicEventTree");
188 
189  tEventTree->Branch(
190  "event",
191  &cEventVals,
192  "runNumber/I:eventNumber/I:nHitsTotal_Unknown/I:nHitsTotal_Cosmic/I:nHitsTotal_NonCosmic/"
193  "I:qTotal_Unknown/F:qTotal_Cosmic/F:qTotal_NonCosmic/F:nHitsTrack/I:nHitsTrack_Cosmic/"
194  "I:nHitsTrack_NonCosmic/I:qTrack/F:qTrack_Cosmic/F:qTrack_NonCosmic/F:nHitsCluster/"
195  "I:nHitsCluster_Cosmic/I:nHitsCluster_NonCosmic/I:qCluster/F:qCluster_Cosmic/"
196  "F:qCluster_NonCosmic/F:TotalTaggedCharge_Cosmic/F:TotalTaggedCharge_NonCosmic/"
197  "F:TotalTaggedHits_Cosmic/I:TotalTaggedHits_NonCosmic/I");
198  tEventTree->Branch("TaggedCharge_Cosmic", &cTaggedCharge_Cosmic);
199  tEventTree->Branch("TaggedCharge_NonCosmic", &cTaggedCharge_NonCosmic);
200  tEventTree->Branch("TaggedHits_Cosmic", &cTaggedHits_Cosmic);
201  tEventTree->Branch("TaggedHits_NonCosmic", &cTaggedHits_NonCosmic);
202 }
203 
204 // =====================================================
205 // Event Loop
206 // =====================================================
208 {
209 
210  InitEventTree(evt.run(), evt.event());
211 
212  // ##################################################################
213  // ### Grabbing ALL HITS in the event to monitor the backtracking ###
214  // ##################################################################
216  evt.getByLabel(fHitsModuleLabel, hitListHandle);
217  std::vector<recob::Hit> const& hitVector(*hitListHandle);
218 
219  std::vector<hit_origin_t> hitOrigins(hitVector.size());
220 
221  //get mcHitCollection
223  evt.getByLabel(fMCHitsModuleLabel, mchitListHandle);
224  std::vector<sim::MCHitCollection> const& mchitcolVector(*mchitListHandle);
225 
226  //get mcparticles out of the event
228  evt.getByLabel(fMCModuleLabel, mcParticleHandle);
229  std::vector<simb::MCParticle> const& mcParticleVector(*mcParticleHandle);
230 
231  //get associations of mc particles to mc truth
233  evt.getByLabel(fMCModuleLabel, assnMCParticleTruthHandle);
234 
235  std::vector<const simb::MCTruth*> particle_to_truth =
236  util::GetAssociatedVectorOneP(assnMCParticleTruthHandle, mcParticleHandle);
237 
238  //make trackId to MCTruth map
239  std::map<int, const simb::MCTruth*> trackIDToTruthMap;
240  for (size_t p_iter = 0; p_iter < mcParticleVector.size(); p_iter++)
241  trackIDToTruthMap[mcParticleVector[p_iter].TrackId()] = particle_to_truth[p_iter];
242 
243  FillMCInfo(evt, hitVector, hitOrigins, mchitcolVector, trackIDToTruthMap);
244 
245  art::Handle<std::vector<recob::Track>> trackListHandle;
246  evt.getByLabel(fTrackModuleLabel, trackListHandle);
247  std::vector<recob::Track> const& trackVector(*trackListHandle);
248 
249  art::Handle<std::vector<recob::Cluster>> clusterListHandle;
250  evt.getByLabel(fClusterModuleLabel, clusterListHandle);
251  std::vector<recob::Cluster> const& clusterVector(*clusterListHandle);
252 
254  evt.getByLabel(fTrackModuleLabel, assnHitTrackHandle);
255  std::vector<std::vector<size_t>> track_indices_per_hit =
256  util::GetAssociatedVectorManyI(assnHitTrackHandle, hitListHandle);
257 
259  evt.getByLabel(fClusterModuleLabel, assnHitClusterHandle);
260  std::vector<std::vector<size_t>> cluster_indices_per_hit =
261  util::GetAssociatedVectorManyI(assnHitClusterHandle, hitListHandle);
262 
263  std::vector<art::Handle<std::vector<anab::CosmicTag>>> cosmicTagHandlesVector(
264  fCosmicTagAssocLabel.size());
265  std::vector<art::Handle<art::Assns<recob::Track, anab::CosmicTag>>> assnTrackTagHandlesVector(
266  fCosmicTagAssocLabel.size());
267  std::vector<std::vector<const anab::CosmicTag*>> tags_per_track(
268  trackVector.size(), std::vector<const anab::CosmicTag*>(fCosmicTagAssocLabel.size()));
269  std::vector<art::Handle<art::Assns<recob::Cluster, anab::CosmicTag>>> assnClusterTagHandlesVector(
270  fCosmicTagAssocLabel.size());
271  std::vector<std::vector<const anab::CosmicTag*>> tags_per_cluster(
272  clusterVector.size(), std::vector<const anab::CosmicTag*>(fCosmicTagAssocLabel.size()));
273 
274  for (size_t label_i = 0; label_i < fCosmicTagAssocLabel.size(); label_i++) {
275  try {
276  evt.getByLabel(fCosmicTagAssocLabel[label_i], cosmicTagHandlesVector[label_i]);
277  }
278  catch (...) {
279  continue;
280  }
281  try {
282  evt.getByLabel(fCosmicTagAssocLabel[label_i], assnTrackTagHandlesVector[label_i]);
283  for (auto const& pair : *assnTrackTagHandlesVector[label_i])
284  tags_per_track.at(pair.first.key())[label_i] = &(*(pair.second));
285  }
286  catch (...) {
287  }
288  try {
289  evt.getByLabel(fCosmicTagAssocLabel[label_i], assnClusterTagHandlesVector[label_i]);
290  for (auto const& pair : *assnClusterTagHandlesVector[label_i])
291  tags_per_cluster.at(pair.first.key())[label_i] = &(*(pair.second));
292  }
293  catch (...) {
294  }
295  }
296 
297  std::vector<std::vector<bool>> hitsAccounted(
298  hitVector.size(), std::vector<bool>(fCosmicTagAssocLabel.size(), false));
299  std::vector<bool> hitsAllTags(hitVector.size(), false);
300 
301  for (size_t hit_iter = 0; hit_iter < hitVector.size(); hit_iter++) {
302 
303  float charge = hitVector[hit_iter].Integral();
304  hit_origin_t origin = hitOrigins[hit_iter];
305 
306  if (track_indices_per_hit[hit_iter].size() != 0)
307  FillTrackInfo(hit_iter,
308  origin,
309  charge,
310  track_indices_per_hit[hit_iter],
311  tags_per_track,
312  hitsAccounted[hit_iter],
313  hitsAllTags);
314 
315  if (cluster_indices_per_hit[hit_iter].size() != 0)
316  FillClusterInfo(hit_iter,
317  origin,
318  charge,
319  cluster_indices_per_hit[hit_iter],
320  tags_per_cluster,
321  hitsAccounted[hit_iter],
322  hitsAllTags);
323 
324  if (hitsAllTags[hit_iter]) FillAllTagsInfo(hitVector[hit_iter], origin);
325 
326  } //end loop over all the hits
327 
328  tEventTree->Fill();
329 }
330 
331 //-------------------------------------------------------------------------------------------------------------------
332 void microboone::CosmicRemovalAna::InitEventTree(int run_number, int event_number)
333 {
334 
335  cEventVals.runNumber = run_number; //evt.run();
336  cEventVals.eventNumber = event_number; //evt.event();
337 
338  cEventVals.nHitsTotal_Unknown = 0;
339  cEventVals.nHitsTotal_Cosmic = 0;
340  cEventVals.nHitsTotal_NonCosmic = 0;
341 
342  cEventVals.qTotal_Unknown = 0;
343  cEventVals.qTotal_Cosmic = 0;
344  cEventVals.qTotal_NonCosmic = 0;
345 
346  cEventVals.nHitsTrack = 0;
347  cEventVals.nHitsTrack_Cosmic = 0;
348  cEventVals.nHitsTrack_NonCosmic = 0;
349 
350  cEventVals.qTrack = 0;
351  cEventVals.qTrack_Cosmic = 0;
352  cEventVals.qTrack_NonCosmic = 0;
353 
354  cEventVals.nHitsCluster = 0;
355  cEventVals.nHitsCluster_Cosmic = 0;
356  cEventVals.nHitsCluster_NonCosmic = 0;
357 
358  cEventVals.qCluster = 0;
359  cEventVals.qCluster_Cosmic = 0;
360  cEventVals.qCluster_NonCosmic = 0;
361 
362  cEventVals.TotalTaggedCharge_Cosmic = 0.;
363  cEventVals.TotalTaggedCharge_NonCosmic = 0.;
364  cEventVals.TotalTaggedHits_Cosmic = 0;
365  cEventVals.TotalTaggedHits_NonCosmic = 0;
366 
367  for (size_t iter = 0; iter < fCosmicTagAssocLabel.size(); iter++) {
368  cTaggedHits_Cosmic.at(iter) = 0;
369  cTaggedCharge_Cosmic.at(iter) = 0;
370  cTaggedHits_NonCosmic.at(iter) = 0;
371  cTaggedCharge_NonCosmic.at(iter) = 0;
372  }
373 }
374 
375 //-------------------------------------------------------------------------------------------------------------------
376 // take in a list of hits, and determine the origin for those hits (and fill in the tree info)
378  art::Event const& e,
379  std::vector<recob::Hit> const& hitlist,
380  std::vector<hit_origin_t>& hitOrigins,
381  std::vector<sim::MCHitCollection> const& mchitCollectionVector,
382  std::map<int, const simb::MCTruth*> const& trackIdToTruthMap)
383 {
384  auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(e);
385 
386  for (size_t itr = 0; itr < hitlist.size(); itr++) {
387 
388  recob::Hit const& this_hit = hitlist[itr];
389 
390  std::vector<int> trackIDs;
391  std::vector<double> energy;
392 
393  for (auto const& mchit : mchitCollectionVector[this_hit.Channel()]) {
394  if (std::abs(clock_data.TPCTDC2Tick(mchit.PeakTime()) - this_hit.PeakTime()) <
395  fHitCompareCut) {
396  trackIDs.push_back(mchit.PartTrackId());
397  energy.push_back(mchit.PartEnergy());
398  }
399  }
400 
401  if (trackIDs.size() == 0) {
402  hitOrigins[itr] = hit_origin_Unknown;
403  cEventVals.nHitsTotal_Unknown++;
404  cEventVals.qTotal_Unknown += this_hit.Integral();
405  continue;
406  }
407 
408  float cosmic_energy = 0;
409  float non_cosmic_energy = 0;
410 
411  for (size_t iter = 0; iter < trackIDs.size(); iter++) {
412  auto map_element = trackIdToTruthMap.find(std::abs(trackIDs[iter]));
413  if (map_element == trackIdToTruthMap.end()) continue;
414  int origin = map_element->second->Origin();
415  if (origin == simb::kBeamNeutrino)
416  non_cosmic_energy += energy[iter];
417  else
418  cosmic_energy += energy[iter];
419  }
420 
421  if (non_cosmic_energy > cosmic_energy) {
422  hitOrigins[itr] = hit_origin_NonCosmic;
423  cEventVals.nHitsTotal_NonCosmic++;
424  cEventVals.qTotal_NonCosmic += this_hit.Integral();
425  }
426  else {
427  hitOrigins[itr] = hit_origin_Cosmic;
428  cEventVals.nHitsTotal_Cosmic++;
429  cEventVals.qTotal_Cosmic += this_hit.Integral();
430  }
431  }
432 
433 } //end FillMCInfo
434 
435 //-------------------------------------------------------------------------------------------------------------------
437  size_t const& hit_iter,
438  hit_origin_t const& origin,
439  float const& charge,
440  std::vector<size_t> const& track_indices_this_hit,
441  std::vector<std::vector<const anab::CosmicTag*>> const& tags_per_track,
442  std::vector<bool>& hitsAccounted_per_tag,
443  std::vector<bool>& hitsAllTags)
444 {
445 
446  cEventVals.nHitsTrack++;
447  cEventVals.qTrack += charge;
448 
449  if (origin == hit_origin_Cosmic) {
450  cEventVals.nHitsTrack_Cosmic++;
451  cEventVals.qTrack_Cosmic += charge;
452  }
453  else if (origin == hit_origin_NonCosmic) {
454  cEventVals.nHitsTrack_NonCosmic++;
455  cEventVals.qTrack_NonCosmic += charge;
456  }
457 
458  for (unsigned int nCT = 0; nCT < fCosmicTagAssocLabel.size();
459  nCT++) { //<---This loops over the vector of cosmicTags in stored in the event
460  if (hitsAccounted_per_tag[nCT]) continue;
461 
462  for (auto const& track_index : track_indices_this_hit) {
463  if (!tags_per_track[track_index][nCT]) continue;
464  const anab::CosmicTag* currentTag(tags_per_track[track_index][nCT]);
465  if (currentTag->CosmicScore() > fCosmicScoreThresholds[nCT]) {
466 
467  hitsAccounted_per_tag[nCT] = true;
468  hitsAllTags[hit_iter] = true;
469  if (origin == hit_origin_Cosmic) {
470  cTaggedHits_Cosmic[nCT]++;
471  cTaggedCharge_Cosmic[nCT] += charge;
472  }
473  else if (origin == hit_origin_NonCosmic) {
474  cTaggedHits_NonCosmic[nCT]++;
475  cTaggedCharge_NonCosmic[nCT] += charge;
476  }
477  }
478  }
479  }
480 
481 } //end FillTrackInfo
482 
483 //-------------------------------------------------------------------------------------------------------------------
485  size_t const& hit_iter,
486  hit_origin_t const& origin,
487  float const& charge,
488  std::vector<size_t> const& cluster_indices_this_hit,
489  std::vector<std::vector<const anab::CosmicTag*>> const& tags_per_cluster,
490  std::vector<bool>& hitsAccounted_per_tag,
491  std::vector<bool>& hitsAllTags)
492 {
493 
494  cEventVals.nHitsCluster++;
495  cEventVals.qCluster += charge;
496 
497  if (origin == hit_origin_Cosmic) {
498  cEventVals.nHitsCluster_Cosmic++;
499  cEventVals.qCluster_Cosmic += charge;
500  }
501  else if (origin == hit_origin_NonCosmic) {
502  cEventVals.nHitsCluster_NonCosmic++;
503  cEventVals.qCluster_NonCosmic += charge;
504  }
505 
506  for (unsigned int nCT = 0; nCT < fCosmicTagAssocLabel.size();
507  nCT++) { //<---This loops over the vector of cosmicTags in stored in the event
508  if (hitsAccounted_per_tag[nCT]) continue;
509 
510  for (auto const& cluster_index : cluster_indices_this_hit) {
511  if (!tags_per_cluster[cluster_index][nCT]) continue;
512  const anab::CosmicTag* currentTag(tags_per_cluster[cluster_index][nCT]);
513  if (currentTag->CosmicScore() > fCosmicScoreThresholds[nCT]) {
514 
515  hitsAccounted_per_tag[nCT] = true;
516  hitsAllTags[hit_iter] = true;
517  if (origin == hit_origin_Cosmic) {
518  cTaggedHits_Cosmic[nCT]++;
519  cTaggedCharge_Cosmic[nCT] += charge;
520  }
521  else if (origin == hit_origin_NonCosmic) {
522  cTaggedHits_NonCosmic[nCT]++;
523  cTaggedCharge_NonCosmic[nCT] += charge;
524  }
525  }
526  }
527  }
528 
529 } //end FillClusterInfo
530 
531 //-------------------------------------------------------------------------------------------------------------------
533  hit_origin_t const& origin)
534 {
535  if (origin == hit_origin_Cosmic) {
536  cEventVals.TotalTaggedCharge_Cosmic += hit.Integral();
537  cEventVals.TotalTaggedHits_Cosmic++;
538  }
539  else if (origin == hit_origin_NonCosmic) {
540  cEventVals.TotalTaggedCharge_NonCosmic += hit.Integral();
541  cEventVals.TotalTaggedHits_NonCosmic++;
542  }
543 }
544 
545 namespace microboone {
546 
548 }
void analyze(const art::Event &evt)
read access to event
CosmicRemovalAna(fhicl::ParameterSet const &pset)
void FillAllTagsInfo(recob::Hit const &hit, hit_origin_t const &origin)
void FillTrackInfo(size_t const &hit_iter, hit_origin_t const &origin, float const &charge, std::vector< size_t > const &track_indices_this_hit, std::vector< std::vector< const anab::CosmicTag * >> const &tags_per_cluster, std::vector< bool > &hitsAccounted_per_tag, std::vector< bool > &hitsAllTags)
Declaration of signal hit object.
cEventProperties_t cEventVals
constexpr auto abs(T v)
Returns the absolute value of the argument.
STL namespace.
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:244
std::vector< const U * > GetAssociatedVectorOneP(art::Handle< art::Assns< T, U >> h, art::Handle< std::vector< T >> index_p)
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.cc:6
Particle class.
float & CosmicScore()
Definition: CosmicTag.h:54
std::vector< float > cTaggedCharge_NonCosmic
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
void FillMCInfo(art::Event const &e, std::vector< recob::Hit > const &hitlist, std::vector< hit_origin_t > &hitOrigins, std::vector< sim::MCHitCollection > const &mchitCollectionVector, std::map< int, const simb::MCTruth * > const &trackIDToTruthMap)
double energy
Definition: plottest35.C:25
void InitEventTree(int run_number, int event_number)
Provides recob::Track data product.
EventNumber_t event() const
Definition: Event.cc:41
Declaration of cluster object.
std::vector< float > fCosmicScoreThresholds
Detector simulation of raw signals on wires.
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:220
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Definition: StdUtils.h:120
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Utility object to perform functions of association.
std::vector< float > cTaggedCharge_Cosmic
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
TCEvent evt
Definition: DataStructs.cxx:8
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.cc:29
std::vector< std::vector< size_t > > GetAssociatedVectorManyI(art::Handle< art::Assns< T, U >> h, art::Handle< std::vector< T >> index_p)
raw::ChannelID_t Channel() const
ID of the readout channel the hit was extracted from.
Definition: Hit.h:268
constexpr Point origin()
Returns a origin position with a point of the specified type.
Definition: geo_vectors.h:229
Beam neutrinos.
Definition: MCTruth.h:23
void FillClusterInfo(size_t const &hit_iter, hit_origin_t const &origin, float const &charge, std::vector< size_t > const &cluster_indices_this_hit, std::vector< std::vector< const anab::CosmicTag * >> const &tags_per_cluster, std::vector< bool > &hitsAccounted_per_tag, std::vector< bool > &hitsAllTags)
std::vector< std::string > fCosmicTagAssocLabel