All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
sbndcode/sbndcode/CRT/CRTAna/CRTT0MatchingAna_module.cc
Go to the documentation of this file.
1 
2 ////////////////////////////////////////////////////////////////////////
3 // Class: CRTT0MatchingAnaAna
4 // Module Type: analyzer
5 // File: CRTT0MatchingAnaAna_module.cc
6 //
7 // Tom Brooks (tbrooks@fnal.gov)
8 ////////////////////////////////////////////////////////////////////////
9 
10 // sbndcode includes
18 
19 // LArSoft includes
26 #include "nusimdata/SimulationBase/MCParticle.h"
27 #include "larcorealg/CoreUtils/NumericUtils.h" // util::absDiff()
28 
29 // Framework includes
30 #include "art/Framework/Core/EDAnalyzer.h"
31 #include "art/Framework/Principal/Event.h"
32 #include "art/Framework/Principal/Handle.h"
33 #include "art/Framework/Services/Registry/ServiceHandle.h"
34 #include "art_root_io/TFileService.h"
35 #include "art/Framework/Core/ModuleMacros.h"
36 #include "canvas/Persistency/Common/FindManyP.h"
37 #include "canvas/Utilities/Exception.h"
40 
41 // Utility libraries
42 #include "messagefacility/MessageLogger/MessageLogger.h"
43 #include "fhiclcpp/ParameterSet.h"
44 #include "fhiclcpp/types/Table.h"
45 #include "fhiclcpp/types/Atom.h"
46 #include "cetlib/pow.h" // cet::sum_of_squares()
47 
48 // ROOT includes. Note: To look up the properties of the ROOT classes,
49 // use the ROOT web site; e.g.,
50 // <https://root.cern.ch/doc/master/annotated.html>
51 #include "TH1.h"
52 #include "TH2.h"
53 #include "TVector3.h"
54 
55 // C++ includes
56 #include <map>
57 #include <vector>
58 #include <string>
59 
60 namespace sbnd {
61 
62  class CRTT0MatchingAna : public art::EDAnalyzer {
63  public:
64 
65  // Describes configuration parameters of the module
66  struct Config {
67  using Name = fhicl::Name;
68  using Comment = fhicl::Comment;
69 
70  // One Atom for each parameter
71  fhicl::Atom<art::InputTag> SimModuleLabel {
72  Name("SimModuleLabel"),
73  Comment("tag of detector simulation data product")
74  };
75 
76  fhicl::Atom<art::InputTag> CRTHitLabel {
77  Name("CRTHitLabel"),
78  Comment("tag of CRT hit producer data product")
79  };
80 
81  fhicl::Atom<art::InputTag> TPCTrackLabel {
82  Name("TPCTrackLabel"),
83  Comment("tag of tpc track producer data product")
84  };
85 
86  fhicl::Atom<bool> Verbose {
87  Name("Verbose"),
88  Comment("Print information about what's going on")
89  };
90 
91  fhicl::Table<CRTT0MatchAlg::Config> CRTT0Alg {
92  Name("CRTT0Alg"),
93  };
94 
95  fhicl::Table<CRTBackTracker::Config> CrtBackTrack {
96  Name("CrtBackTrack"),
97  };
98 
99  }; // Config
100 
101  using Parameters = art::EDAnalyzer::Table<Config>;
102 
103  // Constructor: configures module
104  explicit CRTT0MatchingAna(Parameters const& config);
105 
106  // Called once, at start of the job
107  virtual void beginJob() override;
108 
109  // Called once per event
110  virtual void analyze(const art::Event& event) override;
111 
112  // Called once, at end of the job
113  virtual void endJob() override;
114 
115  // Calculate the distance from the track crossing point to CRT overlap coordinates
116  double DistToCrtHit(TVector3 trackPos, sbn::crt::CRTHit crtHit);
117 
118  private:
119 
120  // fcl file parameters
121  art::InputTag fSimModuleLabel; ///< name of detsim producer
122  art::InputTag fCRTHitLabel; ///< name of CRT producer
123  art::InputTag fTPCTrackLabel; ///< name of CRT producer
124  bool fVerbose; ///< print information about what's going on
125 
127 
130 
132 
133  // Histograms
134  std::map<std::string, TH1D*> hDCA;
135  std::map<std::string, TH1D*> hMatchDCA;
136  std::map<std::string, TH1D*> hNoMatchDCA;
137 
138  std::map<std::string, TH1D*> hDoL;
139  std::map<std::string, TH1D*> hMatchDoL;
140  std::map<std::string, TH1D*> hNoMatchDoL;
141 
142  std::map<std::string, TH1D*> hT0;
143  std::map<std::string, TH1D*> hMatchT0;
144  std::map<std::string, TH1D*> hNoMatchT0;
145 
146  std::map<std::string, TH1D*> hEffDCATotal;
147  std::map<std::string, TH1D*> hEffDCAReco;
148  std::map<std::string, TH1D*> hEffDoLTotal;
149  std::map<std::string, TH1D*> hEffDoLReco;
150  std::map<std::string, TH1D*> hEffLengthTotal;
151  std::map<std::string, TH1D*> hEffLengthReco;
152 
153  std::map<std::string, TH1D*> hPurityDCATotal;
154  std::map<std::string, TH1D*> hPurityDCAReco;
155  std::map<std::string, TH1D*> hPurityDoLTotal;
156  std::map<std::string, TH1D*> hPurityDoLReco;
157  std::map<std::string, TH1D*> hPurityLengthTotal;
158  std::map<std::string, TH1D*> hPurityLengthReco;
159 
160  }; // class CRTT0MatchingAna
161 
162 
163  // Constructor
165  : EDAnalyzer(config)
166  , fSimModuleLabel (config().SimModuleLabel())
167  , fCRTHitLabel (config().CRTHitLabel())
168  , fTPCTrackLabel (config().TPCTrackLabel())
169  , fVerbose (config().Verbose())
170  , t0Alg (config().CRTT0Alg())
171  , fCrtBackTrack (config().CrtBackTrack())
172  {
173 
174  } //CRTT0MatchingAna()
175 
176 
178  {
179 
180  // Access tfileservice to handle creating and writing histograms
181  art::ServiceHandle<art::TFileService> tfs;
182  for(size_t i = 0; i < fCrtGeo.NumTaggers() + 1; i++){
183  std::string tagger = "All";
184  if(i < fCrtGeo.NumTaggers()){
185  tagger = fCrtGeo.GetTagger(i).name;
186  }
187  hDCA[tagger] = tfs->make<TH1D>(Form("DCA_%s", tagger.c_str()), "", 50, 0, 100);
188  hMatchDCA[tagger] = tfs->make<TH1D>(Form("MatchDCA_%s", tagger.c_str()), "", 50, 0, 100);
189  hNoMatchDCA[tagger] = tfs->make<TH1D>(Form("NoMatchDCA_%s", tagger.c_str()), "", 50, 0, 100);
190 
191  hDoL[tagger] = tfs->make<TH1D>(Form("DoL_%s", tagger.c_str()), "", 100, 0, 0.25);
192  hMatchDoL[tagger] = tfs->make<TH1D>(Form("MatchDoL_%s", tagger.c_str()), "", 100, 0, 0.25);
193  hNoMatchDoL[tagger] = tfs->make<TH1D>(Form("NoMatchDoL_%s", tagger.c_str()), "", 100, 0, 0.25);
194 
195  hT0[tagger] = tfs->make<TH1D>(Form("T0_%s", tagger.c_str()), "", 600, -3000, 3000);
196  hMatchT0[tagger] = tfs->make<TH1D>(Form("MatchT0_%s", tagger.c_str()), "", 600, -3000, 3000);
197  hNoMatchT0[tagger] = tfs->make<TH1D>(Form("NoMatchT0_%s", tagger.c_str()), "", 600, -3000, 3000);
198 
199  hEffDCATotal[tagger] = tfs->make<TH1D>(Form("EffDCATotal_%s", tagger.c_str()), "", 50, 0, 100);
200  hEffDCAReco[tagger] = tfs->make<TH1D>(Form("EffDCAReco_%s", tagger.c_str()), "", 50, 0, 100);
201  hEffDoLTotal[tagger] = tfs->make<TH1D>(Form("EffDoLTotal_%s", tagger.c_str()), "", 100, 0, 0.25);
202  hEffDoLReco[tagger] = tfs->make<TH1D>(Form("EffDoLReco_%s", tagger.c_str()), "", 100, 0, 0.25);
203  hEffLengthTotal[tagger] = tfs->make<TH1D>(Form("EffLengthTotal_%s", tagger.c_str()), "", 20, 0, 600);
204  hEffLengthReco[tagger] = tfs->make<TH1D>(Form("EffLengthReco_%s", tagger.c_str()), "", 20, 0, 600);
205 
206  hPurityDCATotal[tagger] = tfs->make<TH1D>(Form("PurityDCATotal_%s", tagger.c_str()), "", 50, 0, 100);
207  hPurityDCAReco[tagger] = tfs->make<TH1D>(Form("PurityDCAReco_%s", tagger.c_str()), "", 50, 0, 100);
208  hPurityDoLTotal[tagger] = tfs->make<TH1D>(Form("PurityDoLTotal_%s", tagger.c_str()), "", 100, 0, 0.25);
209  hPurityDoLReco[tagger] = tfs->make<TH1D>(Form("PurityDoLReco_%s", tagger.c_str()), "", 100, 0, 0.25);
210  hPurityLengthTotal[tagger] = tfs->make<TH1D>(Form("PurityLengthTotal_%s", tagger.c_str()), "", 20, 0, 600);
211  hPurityLengthReco[tagger] = tfs->make<TH1D>(Form("PurityLengthReco_%s", tagger.c_str()), "", 20, 0, 600);
212  }
213 
214  // Initial output
215  if(fVerbose) std::cout<<"----------------- CRT T0 Matching Ana Module -------------------"<<std::endl;
216 
217  } // CRTT0MatchingAna::beginJob()
218 
219 
220  void CRTT0MatchingAna::analyze(const art::Event& event)
221  {
222 
223  // Fetch basic event info
224  if(fVerbose){
225  std::cout<<"============================================"<<std::endl
226  <<"Run = "<<event.run()<<", SubRun = "<<event.subRun()<<", Event = "<<event.id().event()<<std::endl
227  <<"============================================"<<std::endl;
228  }
229 
230  //----------------------------------------------------------------------------------------------------------
231  // GETTING PRODUCTS
232  //----------------------------------------------------------------------------------------------------------
233 
234  // Get g4 particles
235  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
236  auto particleHandle = event.getValidHandle<std::vector<simb::MCParticle>>(fSimModuleLabel);
237 
238  // Get CRT hits from the event
239  art::Handle< std::vector<sbn::crt::CRTHit>> crtHitHandle;
240  std::vector<art::Ptr<sbn::crt::CRTHit> > crtHitList;
241  if (event.getByLabel(fCRTHitLabel, crtHitHandle))
242  art::fill_ptr_vector(crtHitList, crtHitHandle);
243 
244  // Get reconstructed tracks from the event
245  auto tpcTrackHandle = event.getValidHandle<std::vector<recob::Track>>(fTPCTrackLabel);
246  art::FindManyP<recob::Hit> findManyHits(tpcTrackHandle, event, fTPCTrackLabel);
247 
248  fCrtBackTrack.Initialize(event);
249 
250  //----------------------------------------------------------------------------------------------------------
251  // TRUTH MATCHING
252  //----------------------------------------------------------------------------------------------------------
253 
254  std::map<int, simb::MCParticle> particles;
255  // Loop over the true particles
256  for (auto const& particle: (*particleHandle)){
257 
258  // Make map with ID
259  int partID = particle.TrackId();
260  particles[partID] = particle;
261 
262  }
263 
264  std::map<int, std::vector<std::string>> crtTaggerMap;
265  std::vector<sbn::crt::CRTHit> crtHits;
266  int hit_i = 0;
267  double minHitTime = 99999;
268  double maxHitTime = -99999;
269 
270  for(auto const& hit : (*crtHitHandle)){
271  double hitTime = (double)(int)hit.ts1_ns * 1e-3;
272  if(hitTime < minHitTime) minHitTime = hitTime;
273  if(hitTime > maxHitTime) maxHitTime = hitTime;
274 
275  crtHits.push_back(hit);
276  int trueId = fCrtBackTrack.TrueIdFromHitId(event, hit_i);
277  hit_i++;
278  if(trueId == -99999) continue;
279  if(crtTaggerMap.find(trueId) == crtTaggerMap.end()){
280  crtTaggerMap[trueId].push_back(hit.tagger);
281  }
282  else if(std::find(crtTaggerMap[trueId].begin(), crtTaggerMap[trueId].end(), hit.tagger) == crtTaggerMap[trueId].end()){
283  crtTaggerMap[trueId].push_back(hit.tagger);
284  }
285  }
286 
287  // std::cout << " New event" << std::endl;
288  //----------------------------------------------------------------------------------------------------------
289  // DISTANCE OF CLOSEST APPROACH ANALYSIS
290  //----------------------------------------------------------------------------------------------------------
291 
292  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
293  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(event, clockData);
294 
295  // Loop over reconstructed tracks
296  for (auto const& tpcTrack : (*tpcTrackHandle)){
297  // Get the associated hits
298  std::vector<art::Ptr<recob::Hit>> hits = findManyHits.at(tpcTrack.ID());
299  int trackTrueID = RecoUtils::TrueParticleIDFromTotalRecoHits(clockData, hits, false);
300  if(particles.find(trackTrueID) == particles.end()) continue;
301  // Only consider primary muons
302  if(!(std::abs(particles[trackTrueID].PdgCode()) == 13 && particles[trackTrueID].Mother() == 0)) continue;
303 
304  // Only consider particles inside the reco hit time window
305  double trueTime = particles[trackTrueID].T() * 1e-3;
306  if(trueTime < minHitTime || trueTime > maxHitTime) continue;
307 
308  // std::cout << "new track " << trueTime << std::endl;
309  // Calculate t0 from CRT Hit matching
310  matchCand closest = t0Alg.GetClosestCRTHit(detProp, tpcTrack, crtHits, event);
311  // std::cout << "closest match " << closest.t0 << std::endl;
312  double sin_angle = -99999;
313  if(closest.dca != -99999){
314  hDCA[closest.thishit.tagger]->Fill(closest.dca);
315  hDCA["All"]->Fill(closest.dca);
316  sin_angle = closest.dca/closest.extrapLen;
317  hDoL[closest.thishit.tagger]->Fill(sin_angle);
318  hDoL["All"]->Fill(sin_angle);
319  hT0[closest.thishit.tagger]->Fill(closest.t0);
320  hT0["All"]->Fill(closest.t0);
321 
322 
323  // Is hit matched to that track
324  int hitTrueID = fCrtBackTrack.TrueIdFromTotalEnergy(event, closest.thishit);
325  if(hitTrueID == trackTrueID && hitTrueID != -99999){
326  hMatchDCA[closest.thishit.tagger]->Fill(closest.dca);
327  hMatchDCA["All"]->Fill(closest.dca);
328  hMatchDoL[closest.thishit.tagger]->Fill(sin_angle);
329  hMatchDoL["All"]->Fill(sin_angle);
330  hMatchT0[closest.thishit.tagger]->Fill(closest.t0);
331  hMatchT0["All"]->Fill(closest.t0);
332  }
333  else{
334  hNoMatchDCA[closest.thishit.tagger]->Fill(closest.dca);
335  hNoMatchDCA["All"]->Fill(closest.dca);
336  hNoMatchDoL[closest.thishit.tagger]->Fill(sin_angle);
337  hNoMatchDoL["All"]->Fill(sin_angle);
338  hNoMatchT0[closest.thishit.tagger]->Fill(closest.t0);
339  hNoMatchT0["All"]->Fill(closest.t0);
340  }
341  }
342 
343  int nbins = hEffDCATotal.begin()->second->GetNbinsX();
344  for(int i = 0; i < nbins; i++){
345  double DCAcut = hEffDCATotal.begin()->second->GetBinCenter(i);
346 
347  // Fill total efficiency histogram with each cut if track matches any hits
348  if(crtTaggerMap.find(trackTrueID) != crtTaggerMap.end()){
349  for(auto const& tagger : crtTaggerMap[trackTrueID]){
350  hEffDCATotal[tagger]->Fill(DCAcut);
351 
352  // If closest hit is below limit and track matches any hits then fill efficiency
353  if(closest.dca < DCAcut && closest.dca != -99999){
354  hEffDCAReco[tagger]->Fill(DCAcut);
355  }
356  }
357  // Fill total efficiency histograms
358  hEffDCATotal["All"]->Fill(DCAcut);
359  if(closest.dca < DCAcut && closest.dca != -99999){
360  hEffDCAReco["All"]->Fill(DCAcut);
361  }
362  }
363 
364  // Fill total purity histogram with each cut if closest hit is below limit
365  if(closest.dca < DCAcut && closest.dca != -99999){
366  hPurityDCATotal[closest.thishit.tagger]->Fill(DCAcut);
367  hPurityDCATotal["All"]->Fill(DCAcut);
368 
369  // If closest hit is below limit and matched time is correct then fill purity
370  double hitTime = closest.thishit.ts1_ns * 1e-3;
371  if(particles.find(trackTrueID) != particles.end()){
372  if(std::abs(hitTime - trueTime) < 2.){
373  hPurityDCAReco[closest.thishit.tagger]->Fill(DCAcut);
374  hPurityDCAReco["All"]->Fill(DCAcut);
375  }
376  }
377 
378  }
379  }
380 
381  nbins = hEffDoLTotal.begin()->second->GetNbinsX();
382 
383  for(int i = 0; i < nbins; i++){
384  double DCAcut = hEffDoLTotal.begin()->second->GetBinCenter(i);
385 
386  // Fill total efficiency histogram with each cut if track matches any hits
387  if(crtTaggerMap.find(trackTrueID) != crtTaggerMap.end()){
388  for(auto const& tagger : crtTaggerMap[trackTrueID]){
389  hEffDoLTotal[tagger]->Fill(DCAcut);
390 
391  // If closest hit is below limit and track matches any hits then fill efficiency
392  if(sin_angle < DCAcut && closest.dca != -99999){
393  hEffDoLReco[tagger]->Fill(DCAcut);
394  }
395  }
396  // Fill total efficiency histograms
397  hEffDoLTotal["All"]->Fill(DCAcut);
398  if(sin_angle < DCAcut && closest.dca != -99999){
399  hEffDoLReco["All"]->Fill(DCAcut);
400  }
401  }
402 
403  // Fill total purity histogram with each cut if closest hit is below limit
404  if(sin_angle < DCAcut && closest.dca != -99999){
405  hPurityDoLTotal[closest.thishit.tagger]->Fill(DCAcut);
406  hPurityDoLTotal["All"]->Fill(DCAcut);
407 
408  // If closest hit is below limit and matched time is correct then fill purity
409  double hitTime = closest.thishit.ts1_ns * 1e-3;
410  if(particles.find(trackTrueID) != particles.end()){
411  if(std::abs(hitTime - trueTime) < 2.){
412  hPurityDoLReco[closest.thishit.tagger]->Fill(DCAcut);
413  hPurityDoLReco["All"]->Fill(DCAcut);
414  }
415  }
416 
417  }
418  }
419 
420  double fixedCut = 30.;
421 
422  // Fill total efficiency histogram with each cut if track matches any hits
423  if(crtTaggerMap.find(trackTrueID) != crtTaggerMap.end()){
424  for(auto const& tagger : crtTaggerMap[trackTrueID]){
425  hEffLengthTotal[tagger]->Fill(tpcTrack.Length());
426 
427  // If closest hit is below limit and track matches any hits then fill efficiency
428  if(closest.dca < fixedCut && closest.dca >=0 ){
429  hEffLengthReco[tagger]->Fill(tpcTrack.Length());
430  }
431  }
432  // Fill total efficiency histograms
433  hEffLengthTotal["All"]->Fill(tpcTrack.Length());
434  if(closest.dca < fixedCut && closest.dca >=0){
435  hEffLengthReco["All"]->Fill(tpcTrack.Length());
436  }
437  }
438 
439  // Fill total purity histogram with each cut if closest hit is below limit
440  if(closest.dca < fixedCut && closest.dca >= 0){
441  hPurityLengthTotal[closest.thishit.tagger]->Fill(tpcTrack.Length());
442  hPurityLengthTotal["All"]->Fill(tpcTrack.Length());
443 
444  // If closest hit is below limit and matched time is correct then fill purity
445  double hitTime = closest.thishit.ts1_ns * 1e-3;
446  if(particles.find(trackTrueID) != particles.end()){
447  double trueTime = particles[trackTrueID].T() * 1e-3;
448  if(std::abs(hitTime - trueTime) < 2.){
449  hPurityLengthReco[closest.thishit.tagger]->Fill(tpcTrack.Length());
450  hPurityLengthReco["All"]->Fill(tpcTrack.Length());
451  }
452  }
453 
454  }
455 
456  }
457 
458 
459  } // CRTT0MatchingAna::analyze()
460 
461 
463 
464 
465  } // CRTT0MatchingAna::endJob()
466 
467 
468  DEFINE_ART_MODULE(CRTT0MatchingAna)
469 } // namespace sbnd
BEGIN_PROLOG Verbose
matchCand GetClosestCRTHit(detinfo::DetectorPropertiesData const &detProp, recob::Track tpcTrack, std::pair< double, double > t0MinMax, std::vector< sbn::crt::CRTHit > crtHits, int driftDirection)
Functions to help with numbers.
double ts1_ns
Timestamp T1 ([signal time w.r.t. Trigger time]), in UTC absolute time scale in nanoseconds from the ...
Definition: CRTHit.hh:34
CRTTaggerGeo GetTagger(std::string taggerName) const
Declaration of signal hit object.
virtual void analyze(const art::Event &event) override
process_name hit
Definition: cheaterreco.fcl:51
int TrueIdFromTotalEnergy(const art::Event &event, const sbnd::crt::CRTData &data)
int TrueIdFromHitId(const art::Event &event, int hit_i)
T abs(T value)
double DistToCrtHit(TVector3 trackPos, sbn::crt::CRTHit crtHit)
auto end(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:585
BEGIN_PROLOG vertical distance to the surface Name
int TrueParticleIDFromTotalRecoHits(detinfo::DetectorClocksData const &clockData, const std::vector< art::Ptr< recob::Hit > > &hits, bool rollup_unsaved_ids=1)
BEGIN_PROLOG don t mess with this TPCTrackLabel
Definition of data types for geometry description.
Provides recob::Track data product.
auto begin(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:573
bool fVerbose
print information about what&#39;s going on
object containing MC truth information necessary for making RawDigits and doing back tracking ...
do i e
stream1 can override from command line with o or output services user sbnd
art::ServiceHandle< art::TFileService > tfs
std::string tagger
Name of the CRT wall (in the form of strings).
Definition: CRTHit.hh:45
BEGIN_PROLOG could also be cout
auto const detProp