All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CRTTruthMatchAnalysis_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: CRTTruthMatchAnalysis
3 // Plugin Type: analyzer (art v3_05_00)
4 // File: CRTTruthMatchAnalysis_module.cc
5 //
6 // Generated at Fri Apr 24 14:45:31 2020 by Christopher Hilgenberg using cetskelgen
7 // from cetlib version v3_10_00.
8 ////////////////////////////////////////////////////////////////////////
9 
10 //Framework includes
11 #include "art/Framework/Core/EDAnalyzer.h"
12 #include "art/Framework/Core/ModuleMacros.h"
13 #include "art/Framework/Principal/Event.h"
14 #include "art/Framework/Principal/Handle.h"
15 #include "art/Framework/Principal/Run.h"
16 #include "art/Framework/Principal/SubRun.h"
17 #include "canvas/Utilities/InputTag.h"
18 #include "fhiclcpp/ParameterSet.h"
19 #include "messagefacility/MessageLogger/MessageLogger.h"
20 #include "art_root_io/TFileService.h"
21 
22 //LarSoft includes
23 #include "nusimdata/SimulationBase/MCParticle.h"
24 #include "lardataobj/Simulation/AuxDetSimChannel.h"
25 
26 //icaruscode includes
27 #include "sbnobj/ICARUS/CRT/CRTData.hh"
28 #include "sbnobj/Common/CRT/CRTHit.hh"
29 #include "icaruscode/CRT/CRTUtils/CRTBackTracker.h"
30 #include "icaruscode/CRT/CRTUtils/CRTCommonUtils.h"
31 
32 //C++ includes
33 #include <vector>
34 #include <map>
35 #include <string>
36 #include <utility>
37 
38 //ROOT includes
39 #include <TH1F.h>
40 #include <TH2F.h>
41 #include <TGraph.h>
42 #include <TVector3.h>
43 
44 namespace icarus{
45  namespace crt {
46  class CRTTruthMatchAnalysis;
47  }
48 }
49 
50 using namespace icarus::crt;
51 using std::vector;
52 using std::map;
53 
54 class icarus::crt::CRTTruthMatchAnalysis : public art::EDAnalyzer {
55 
56  public:
57  using CRTHit = sbn::crt::CRTHit;
58 
59  explicit CRTTruthMatchAnalysis(fhicl::ParameterSet const& p);
60  // The compiler-generated destructor is fine for non-base
61  // classes without bare pointers or other resource use.
62 
63  // Plugins should not be copied or assigned.
64  CRTTruthMatchAnalysis(CRTTruthMatchAnalysis const&) = delete;
65  CRTTruthMatchAnalysis(CRTTruthMatchAnalysis&&) = delete;
66  CRTTruthMatchAnalysis& operator=(CRTTruthMatchAnalysis const&) = delete;
67  CRTTruthMatchAnalysis& operator=(CRTTruthMatchAnalysis&&) = delete;
68 
69  void beginJob() override;
70  void analyze(art::Event const& e) override;
71 
72  private:
73 
74  float GetDeltaR(art::Ptr<CRTHit> h1, art::Ptr<CRTHit> h2);
75  int GetDeltaT(art::Ptr<CRTHit> h1, art::Ptr<CRTHit> h2);
76  TVector3 HitToPosVec(art::Ptr<CRTHit> hit);
77 
78  art::InputTag fSimulationLabel;
79  art::InputTag fCRTTrueHitLabel;
80  art::InputTag fCRTDataLabel;
81  art::InputTag fCRTSimHitLabel;
82 
83  CRTBackTracker bt;
84  CRTCommonUtils* fCrtutils;
85 
86  const float fPosmax;
87  const float fPosbinning;
88  const float fTimemax;
89  const float fTbinning;
90 
91  const bool fVerbose;
92 
93  TH1F* fNTrueHitMu; //number of CRTTrueHits per muon track
94  TH1F* fNDataMu; //number of CRTData per muon track
95  TH1F* fNSimHitMu; //number of CRTSimHits per muon track
96  TH1F* fNTrueSimDiffMu; //difference in number of true - sim hits
97 
98  TH1F* fNTrueHitMu_top;
99  TH1F* fNTrueHitMu_side;
100  TH1F* fNTrueHitMu_bot;
101 // TH1F* fNDataMu_top;
102 // TH1F* fNDataMu_side;
103 // TH1F* fNDataMu_bot;
104  TH1F* fNSimHitMu_top;
105  TH1F* fNSimHitMu_side;
106  TH1F* fNSimHitMu_bot;
107 
108  TH1F* fNRegTrue;
109  TH1F* fNRegSim;
110 
111  TH1F* fFebsTrue;
112  TH1F* fFebsSim;
113 
114  TH1F *fRres_top, *fXres_top, *fYres_top, *fZres_top, *fTres_top;
115  TH1F *fRres_rim, *fXres_rim, *fYres_rim, *fZres_rim, *fTres_rim;
116  TH1F *fRres_lat, *fXres_lat, *fYres_lat, *fZres_lat, *fTres_lat;
117  TH1F *fRres_nor, *fXres_nor, *fYres_nor, *fZres_nor, *fTres_nor;
118  TH1F *fRres_sth, *fXres_sth, *fYres_sth, *fZres_sth, *fTres_sth;
119  TH1F *fRres_bot, *fXres_bot, *fYres_bot, *fZres_bot, *fTres_bot;
120 
121  //size_t fNViews;
122  //vector<TH2F*> fChargeViews;
123  //vector<TGraph*> fTrueHitPoints;
124  //vector<TGraph*> fSimHitPoints;
125 
126 };//class definition
127 
128 //constructor
129 CRTTruthMatchAnalysis::CRTTruthMatchAnalysis(fhicl::ParameterSet const& p)
130  : EDAnalyzer{p},
131  fSimulationLabel(p.get<art::InputTag>("SimulationLabel","largeant")),
132  fCRTTrueHitLabel(p.get<art::InputTag>("CRTTrueHitLabel","crttruehit")),
133  fCRTDataLabel(p.get<art::InputTag>("CRTDataLabel","crtdaq")),
134  fCRTSimHitLabel(p.get<art::InputTag>("CRTSimHitLabel","crthit")),
135  bt(p.get<fhicl::ParameterSet>("CRTBackTrack")),
136  fCrtutils(new CRTCommonUtils()),
137  fPosmax(p.get<float>("PosMax",1000.)),
138  fPosbinning(p.get<float>("PosBinning",1.)),
139  fTimemax(p.get<float>("TimeMax",30.)),
140  fTbinning(p.get<float>("Tbinning",1.)),
141  fVerbose(p.get<bool>("Verbose",false))
142 {}
143 
144 void CRTTruthMatchAnalysis::beginJob()
145 {
146  // Access ART's TFileService, which will handle creating and writing
147  // histograms and n-tuples for us.
148  art::ServiceHandle<art::TFileService> tfs;
149 
150  fNTrueHitMu = tfs->make<TH1F>("hNTrueHitMu", "No. True Hits per #mu", 10,0,10);
151  fNDataMu = tfs->make<TH1F>("hNDataMu", "No. FEB Triggers per #mu", 10,0,10);
152  fNSimHitMu = tfs->make<TH1F>("hNSimHitMu", "No. Sim Hits per #mu", 10,0,10);
153  fNTrueSimDiffMu = tfs->make<TH1F>("hNTrueSimDiffMu", "No. True-Sim Hits per #mu", 14,-6,8);
154 
155  fNTrueHitMu_top = tfs->make<TH1F>("hNTrueHitMu_top", "No. True Hits per #mu", 10,0,10);
156  fNTrueHitMu_side = tfs->make<TH1F>("hNTrueHitMu_side","No. True Hits per #mu", 10,0,10);
157  fNTrueHitMu_bot = tfs->make<TH1F>("hNTrueHitMu_bot", "No. True Hits per #mu", 10,0,10);
158  fNSimHitMu_top = tfs->make<TH1F>("hNSimHitMu_top", "No. Sim Hits per #mu", 10,0,10);
159  fNSimHitMu_side = tfs->make<TH1F>("hNSimHitMu_side", "No. Sim Hits per #mu", 10,0,10);
160  fNSimHitMu_bot = tfs->make<TH1F>("hNSimHitMu_bot", "No. Sim Hits per #mu", 10,0,10);
161 
162  fNRegTrue = tfs->make<TH1F>("hNRegTrue", "No. True Hits per Region", 10,0,10);
163  fNRegSim = tfs->make<TH1F>("hNRegSim", "No. Sim Hits per Region", 10,0,10);
164 
165  fFebsTrue = tfs->make<TH1F>("hFebsTrue", "mac5s in true hits", 300,0,300);
166  fFebsSim = tfs->make<TH1F>("hFebsSim", "mac5s in sim hits", 300,0,300);
167 
168  const float absmax = fPosmax+0.5*fPosbinning;
169  const float absmin = -1*absmax;
170  const int nbinsabs = 2*absmax/fPosbinning;
171 
172  const float rmax = fPosmax;
173  const int nbinsr = rmax/fPosbinning;
174 
175  const float tmax = fTimemax+0.5*fTbinning;
176  const float tmin = -1*tmax;
177  const int nbinst = 2*tmax/fTbinning;
178 
179  art::TFileDirectory resdir = tfs->mkdir("resolution");
180  fRres_top = resdir.make<TH1F>("hRres_top", "CRTHit #sigma_{r}: Top", nbinsr,0,rmax);
181  fRres_rim = resdir.make<TH1F>("hRres_rim", "CRTHit #sigma_{r}: Rim", nbinsr,0,rmax);
182  fRres_bot = resdir.make<TH1F>("hRres_bot", "CRTHit #sigma_{r}: Bottom", nbinsr,0,rmax);
183  fRres_lat = resdir.make<TH1F>("hRres_lat", "CRTHit #sigma_{r}: East/West", nbinsr,0,rmax);
184  fRres_nor = resdir.make<TH1F>("hRres_nor", "CRTHit #sigma_{r}: North", nbinsr,0,rmax);
185  fRres_sth = resdir.make<TH1F>("hRres_sth", "CRTHit #sigma_{r}: South", nbinsr,0,rmax);
186 
187  fXres_top = resdir.make<TH1F>("hXres_top", "CRTHit #sigma_{x}: Top", nbinsabs,absmin,absmax);
188  fXres_rim = resdir.make<TH1F>("hXres_rim", "CRTHit #sigma_{x}: Rim", nbinsabs,absmin,absmax);
189  fXres_bot = resdir.make<TH1F>("hXres_bot", "CRTHit #sigma_{x}: Bottom", nbinsabs,absmin,absmax);
190  fXres_lat = resdir.make<TH1F>("hXres_lat", "CRTHit #sigma_{x}: East/West", nbinsabs,absmin,absmax);
191  fXres_nor = resdir.make<TH1F>("hXres_nor", "CRTHit #sigma_{x}: North", nbinsabs,absmin,absmax);
192  fXres_sth = resdir.make<TH1F>("hXres_sth", "CRTHit #sigma_{x}: South", nbinsabs,absmin,absmax);
193 
194  fYres_top = resdir.make<TH1F>("hYres_top", "CRTHit #sigma_{y}: Top", nbinsabs,absmin,absmax);
195  fYres_rim = resdir.make<TH1F>("hYres_rim", "CRTHit #sigma_{y}: Rim", nbinsabs,absmin,absmax);
196  fYres_bot = resdir.make<TH1F>("hYres_bot", "CRTHit #sigma_{y}: Bottom", nbinsabs,absmin,absmax);
197  fYres_lat = resdir.make<TH1F>("hYres_lat", "CRTHit #sigma_{y}: East/West", nbinsabs,absmin,absmax);
198  fYres_nor = resdir.make<TH1F>("hYres_nor", "CRTHit #sigma_{y}: North", nbinsabs,absmin,absmax);
199  fYres_sth = resdir.make<TH1F>("hYres_sth", "CRTHit #sigma_{y}: South", nbinsabs,absmin,absmax);
200 
201  fZres_top = resdir.make<TH1F>("hZres_top", "CRTHit #sigma_{z}: Top", nbinsabs,absmin,absmax);
202  fZres_rim = resdir.make<TH1F>("hZres_rim", "CRTHit #sigma_{z}: Rim", nbinsabs,absmin,absmax);
203  fZres_bot = resdir.make<TH1F>("hZres_bot", "CRTHit #sigma_{z}: Bottom", nbinsabs,absmin,absmax);
204  fZres_lat = resdir.make<TH1F>("hZres_lat", "CRTHit #sigma_{z}: East/West", nbinsabs,absmin,absmax);
205  fZres_nor = resdir.make<TH1F>("hZres_nor", "CRTHit #sigma_{z}: North", nbinsabs,absmin,absmax);
206  fZres_sth = resdir.make<TH1F>("hZres_sth", "CRTHit #sigma_{z}: South", nbinsabs,absmin,absmax);
207 
208  fTres_top = resdir.make<TH1F>("hTres_top", "CRTHit #sigma_{t}: Top", nbinst,tmin,tmax);
209  fTres_rim = resdir.make<TH1F>("hTres_rim", "CRTHit #sigma_{t}: Rim", nbinst,tmin,tmax);
210  fTres_bot = resdir.make<TH1F>("hTres_bot", "CRTHit #sigma_{t}: Bottom", nbinst,tmin,tmax);
211  fTres_lat = resdir.make<TH1F>("hTres_lat", "CRTHit #sigma_{t}: East/West", nbinst,tmin,tmax);
212  fTres_nor = resdir.make<TH1F>("hTres_nor", "CRTHit #sigma_{t}: North", nbinst,tmin,tmax);
213  fTres_sth = resdir.make<TH1F>("hTres_sth", "CRTHit #sigma_{t}: South", nbinst,tmin,tmax);
214 
215  /*art::TFileDirectory viewdir = tfs->mkdir("chargeview");
216 
217  fNViews = 0;
218  for(int i=0; i<100; i++) {
219  std::string gnametrue = "gview_true"+std::to_string(i);
220  std::string gnamesim = "gview_sim"+std::to_string(i);
221  std::string hname = "hview"+std::to_string(i);
222  std::string gtitletrue = "Charge View True Point "+std::to_string(i);
223  std::string gtitlesim = "Charge View Sim Point "+std::to_string(i);
224  std::string htitle = "Charge View " + std::to_string(i);
225  fChargeViews.push_back(viewdir.make<TH2F>(hname.c_str(), htitle.c_str(),16,0,16,16,0,16));
226  fChargeViews.back()->GetXaxis()->SetTitle("X-layer channel");
227  fChargeViews.back()->GetYaxis()->SetTitle("Y-layer channel");
228  fTrueHitPoints.push_back(viewdir.makeAndRegister<TGraph>(gnametrue.c_str(),gtitletrue.c_str(),1));
229  fSimHitPoints.push_back(viewdir.makeAndRegister<TGraph>(gnamesim.c_str(),gtitlesim.c_str(),1));
230  fTrueHitPoints.back()->SetMarkerStyle(47);
231  fTrueHitPoints.back()->SetMarkerSize(10);
232  fTrueHitPoints.back()->SetMarkerColor(kRed);
233  fSimHitPoints.back()->SetMarkerStyle(8);
234  fSimHitPoints.back()->SetMarkerSize(10);
235  fSimHitPoints.back()->SetMarkerColor(kMagenta);
236  }*/
237 
238  //style options
239  fNTrueHitMu->SetLineWidth(2);
240  fNTrueHitMu_top->SetLineWidth(2);
241  fNTrueHitMu_side->SetLineWidth(2);
242  fNTrueHitMu_bot->SetLineWidth(2);
243  fNSimHitMu->SetLineWidth(2);
244  fNSimHitMu_top->SetLineWidth(2);
245  fNSimHitMu_side->SetLineWidth(2);
246  fNSimHitMu_bot->SetLineWidth(2);
247 
248  fNTrueHitMu->SetLineColor(kBlack);
249  fNTrueHitMu_top->SetLineColor(kBlue);
250  fNTrueHitMu_side->SetLineColor(kRed);
251  fNTrueHitMu_bot->SetLineColor(kGreen);
252  fNSimHitMu->SetLineColor(kBlack);
253  fNSimHitMu_top->SetLineColor(kBlue);
254  fNSimHitMu_side->SetLineColor(kRed);
255  fNSimHitMu_bot->SetLineColor(kGreen);
256 
257  fRres_top->SetLineWidth(2);
258  fRres_rim->SetLineWidth(2);
259  fRres_bot->SetLineWidth(2);
260  fRres_lat->SetLineWidth(2);
261  fRres_nor->SetLineWidth(2);
262  fRres_sth->SetLineWidth(2);
263 
264  fXres_top->SetLineWidth(2);
265  fXres_rim->SetLineWidth(2);
266  fXres_bot->SetLineWidth(2);
267  fXres_lat->SetLineWidth(2);
268  fXres_nor->SetLineWidth(2);
269  fXres_sth->SetLineWidth(2);
270 
271  fYres_top->SetLineWidth(2);
272  fYres_rim->SetLineWidth(2);
273  fYres_bot->SetLineWidth(2);
274  fYres_lat->SetLineWidth(2);
275  fYres_nor->SetLineWidth(2);
276  fYres_sth->SetLineWidth(2);
277 
278  fZres_top->SetLineWidth(2);
279  fZres_rim->SetLineWidth(2);
280  fZres_bot->SetLineWidth(2);
281  fZres_lat->SetLineWidth(2);
282  fZres_nor->SetLineWidth(2);
283  fZres_sth->SetLineWidth(2);
284 
285  fTres_top->SetLineWidth(2);
286  fTres_rim->SetLineWidth(2);
287  fTres_bot->SetLineWidth(2);
288  fTres_lat->SetLineWidth(2);
289  fTres_nor->SetLineWidth(2);
290  fTres_sth->SetLineWidth(2);
291 
292  fRres_top->GetXaxis()->SetTitle("#||{#vec{r_{true}} - #vec{r_{reco}}} [cm]");
293  fRres_rim->GetXaxis()->SetTitle("#||{#vec{r_{true}} - #vec{r_{reco}}} [cm]");
294  fRres_bot->GetXaxis()->SetTitle("#||{#vec{r_{true}} - #vec{r_{reco}}} [cm]");
295  fRres_lat->GetXaxis()->SetTitle("#||{#vec{r_{true}} - #vec{r_{reco}}} [cm]");
296  fRres_nor->GetXaxis()->SetTitle("#||{#vec{r_{true}} - #vec{r_{reco}}} [cm]");
297  fRres_sth->GetXaxis()->SetTitle("#||{#vec{r_{true}} - #vec{r_{reco}}} [cm]");
298 
299  fXres_top->GetXaxis()->SetTitle("#||{x_{true} - x_{reco}} [cm]");
300  fXres_rim->GetXaxis()->SetTitle("#||{x_{true} - x_{reco}} [cm]");
301  fXres_bot->GetXaxis()->SetTitle("#||{x_{true} - x_{reco}} [cm]");
302  fXres_lat->GetXaxis()->SetTitle("#||{x_{true} - x_{reco}} [cm]");
303  fXres_nor->GetXaxis()->SetTitle("#||{x_{true} - x_{reco}} [cm]");
304  fXres_sth->GetXaxis()->SetTitle("#||{x_{true} - x_{reco}} [cm]");
305 
306  fYres_top->GetXaxis()->SetTitle("#||{y_{true} - y_{reco}} [cm]");
307  fYres_rim->GetXaxis()->SetTitle("#||{y_{true} - y_{reco}} [cm]");
308  fYres_bot->GetXaxis()->SetTitle("#||{y_{true} - y_{reco}} [cm]");
309  fYres_lat->GetXaxis()->SetTitle("#||{y_{true} - y_{reco}} [cm]");
310  fYres_nor->GetXaxis()->SetTitle("#||{y_{true} - y_{reco}} [cm]");
311  fYres_sth->GetXaxis()->SetTitle("#||{y_{true} - y_{reco}} [cm]");
312 
313  fZres_top->GetXaxis()->SetTitle("#||{z_{true} - z_{reco}} [cm]");
314  fZres_rim->GetXaxis()->SetTitle("#||{z_{true} - z_{reco}} [cm]");
315  fZres_bot->GetXaxis()->SetTitle("#||{z_{true} - z_{reco}} [cm]");
316  fZres_lat->GetXaxis()->SetTitle("#||{z_{true} - z_{reco}} [cm]");
317  fZres_nor->GetXaxis()->SetTitle("#||{z_{true} - z_{reco}} [cm]");
318  fZres_sth->GetXaxis()->SetTitle("#||{z_{true} - z_{reco}} [cm]");
319 
320  fTres_top->GetXaxis()->SetTitle("#||{t_{true} - t_{reco}} [ns]");
321  fTres_rim->GetXaxis()->SetTitle("#||{t_{true} - t_{reco}} [ns]");
322  fTres_bot->GetXaxis()->SetTitle("#||{t_{true} - t_{reco}} [ns]");
323  fTres_lat->GetXaxis()->SetTitle("#||{t_{true} - t_{reco}} [ns]");
324  fTres_nor->GetXaxis()->SetTitle("#||{t_{true} - t_{reco}} [ns]");
325  fTres_sth->GetXaxis()->SetTitle("#||{t_{true} - t_{reco}} [ns]");
326 
327 
328 }
329 
330 //main body
331 void CRTTruthMatchAnalysis::analyze(art::Event const& e)
332 {
333 
334  size_t nmu=0, nmu_dep=0, ntrue=0, ndata=0, nsim=0;
335  size_t nmiss_true = 0, nmiss_data = 0, nmiss_sim = 0;
336 
337  //Intitialize CRTBackTracker
338  //bt.Initialize(e);
339 
340  //MCParticles
341  art::Handle< vector<simb::MCParticle> > mcHandle;
342  map< int, const simb::MCParticle*> particleMap;
343  if (!e.getByLabel(fSimulationLabel, mcHandle))
344  {
345  // If we have no MCParticles at all in an event, then we're in
346  // big trouble. Throw an exception.
347  throw cet::exception("CRTTruthMatchAnalysis")
348  << " No simb::MCParticle objects in this event - "
349  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
350  }
351 
352  //fill particleMap (trackID->MCParticle object)
353  for(auto const& particle : (*mcHandle)) {
354  particleMap[particle.TrackId()] = &particle;
355  }
356  //mf::LogPrint("CRTTruthMatchAnalysis")
357  // << "found " << particleMap.size() << " MCParticles";
358 
359  //AuxDetSimChannels
360  art::Handle< vector<sim::AuxDetSimChannel> > adscHandle;
361  vector< art::Ptr<sim::AuxDetSimChannel> > adscList;
362  if( e.getByLabel(fSimulationLabel,adscHandle) )
363  art::fill_ptr_vector(adscList,adscHandle);
364  else
365  throw cet::exception("CRTTruthMatchAnalysis")
366  << " No sim::AuxDetSimChannel objects in this event - "
367  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
368 
369  vector<int> muIds; //list of muon trackIds passing through CRT
370  for(auto const& adsc : adscList) {
371  for(auto const& ide : adsc->AuxDetIDEs()) {
372  if(particleMap.find(ide.trackID)==particleMap.end())
373  continue;
374  if(abs(particleMap[ide.trackID]->PdgCode())!=13) //muons only
375  continue;
376  muIds.push_back(ide.trackID);
377  nmu_dep++;
378  }
379  }
380  //remove repeat IDs
381  std::sort(muIds.begin(),muIds.end());
382  muIds.erase(std::unique(muIds.begin(),muIds.end()),muIds.end());
383  nmu = muIds.size();
384  //mf::LogInfo("CRTTruthMatchAnalysis")
385  // << "found " << nmu << " muons passing through CRT"
386  // << " with " << nmu_dep << " associated energy deposits ("
387  // << 1.0*nmu_dep/nmu << ")";
388 
389  //CRTTrueHits
390  art::Handle< vector<CRTHit> > trueHitHandle;
391  vector< art::Ptr<CRTHit> > trueHitList;
392  map<int,vector<art::Ptr<CRTHit>>> muToTrueHits;
393 
394  if( e.getByLabel(fCRTTrueHitLabel,trueHitHandle) ) {
395  art::fill_ptr_vector(trueHitList,trueHitHandle);
396 
397  for(auto const& hit : trueHitList) {
398  for(const int id: bt.AllTrueIds(e,*hit)) {
399  if(particleMap.find(id)==particleMap.end())
400  continue;
401  if(abs(particleMap[id]->PdgCode())!=13) //muons only
402  continue;
403  muToTrueHits[id].push_back(hit);
404  }//for trackIDs
405  }//for CRTTrueHits
406 
407  ntrue = muToTrueHits.size();
408  }//if CRTTrueHits
409  else
410  mf::LogWarning("CRTTruthMatchAnalysis") << "no CRTTrueHits found";
411 
412  //CRTData
413  art::Handle< vector<CRTData> > dataHandle;
414  vector< art::Ptr<CRTData> > dataList;
415  map<int,vector<art::Ptr<CRTData>>> muToData;
416  if( e.getByLabel(fCRTDataLabel,dataHandle)) {
417  art::fill_ptr_vector(dataList,dataHandle);
418  //mf::LogPrint("CRTTruthMatchAnalysis")
419  // << "found " << dataList.size() << " CRTData";
420 
421  for(auto const& data : dataList) {
422  for(const int id: bt.AllTrueIds(e,*data)) {
423  if(particleMap.find(id)==particleMap.end())
424  continue;
425  if(abs(particleMap[id]->PdgCode())!=13) //muons only
426  continue;
427  muToData[id].push_back(data);
428  ndata++;
429  }//for trackIDs
430  }//for data
431 
432  for(const int trk : muIds) {
433  if(muToData.find(trk)==muToData.end()){
434  nmiss_data++;
435  fNDataMu->Fill(0);
436  }
437  else
438  fNDataMu->Fill(muToData[trk].size());
439  }
440 
441  }//if data
442  else
443  mf::LogWarning("CRTTruthMatchAnalysis") << "no CRTData found";
444 
445  //CRTSimHits
446  art::Handle< vector<CRTHit> > simHitHandle;
447  vector< art::Ptr<CRTHit> > simHitList;
448  map<int,vector<art::Ptr<CRTHit>>> muToSimHits;
449  if( e.getByLabel(fCRTSimHitLabel,simHitHandle) ) {
450  art::fill_ptr_vector(simHitList,simHitHandle);
451 
452  for(auto const& hit : simHitList) {
453  for(const int id: bt.AllTrueIds(e,*hit)) {
454  if(particleMap.find(id)==particleMap.end())
455  continue;
456  if(abs(particleMap[id]->PdgCode())!=13) //muons only
457  continue;
458  muToSimHits[id].push_back(hit);
459  nsim++;
460  }//for trackIDs
461  }//for CRTSimHits
462 
463  if(fVerbose)
464  mf::LogPrint("CRTTruthMatchAnalysis")
465  << " found " << nsim << " muon-associated CRTSimHits ("
466  << 100.0*nsim/simHitList.size() << " %)";
467  }//if CRTSimHits
468  else
469  mf::LogWarning("CRTTruthMatchAnalysis") << "no CRTSimHits found";
470 
471  //====tallying matches==========
472  // loop over muon IDs
473  size_t nmatch=0, nmiss=0, nmiss_track=0, nmiss_feb=0;
474  size_t nmatch_top=0;
475  for(const int trk : muIds) {
476  //counters
477  int ntruehit=0, nsimhit=0;
478  int ntrue_top=0, ntrue_side=0, ntrue_bot=0;
479  int nsim_top=0, nsim_side=0, nsim_bot=0;
480  map<string,size_t> nregtrue, nregsim;
481 
482  bool simhasmu = true;
483  //True Hits
484  if(muToTrueHits.find(trk)!=muToTrueHits.end()){
485  ntruehit = muToTrueHits[trk].size();
486  for(auto const& hit : muToTrueHits[trk]) {
487  //mac addresses
488  for(auto const& id : hit->feb_id)
489  fFebsTrue->Fill((int)id);
490  //region
491  string truereg = hit->tagger;
492  if(nregtrue.find(truereg)==nregtrue.end())
493  nregtrue[truereg] = 1;
494  else
495  nregtrue[truereg]++;
496  //type
497  char truetype = fCrtutils->GetRegTypeFromRegName(truereg);
498  if(truetype=='c') ntrue_top++;
499  if(truetype=='m') ntrue_side++;
500  if(truetype=='d') ntrue_bot++;
501 
502  //Sim Hits
503  if(muToSimHits.find(trk)!=muToSimHits.end()) {
504  for(auto const& simhit : muToSimHits[trk]) {
505  std::string simreg = simhit->tagger; //region
506  if(simreg != truereg) continue; //make same region
507  if(truetype == 'c' || truetype=='d'){
508  if(simhit->feb_id[0] == hit->feb_id[0]){ //same mac5
509 
510  if(simreg=="Top") {
511  fRres_top->Fill(GetDeltaR(hit,simhit));
512  fXres_top->Fill(hit->x_pos - simhit->x_pos);
513  fYres_top->Fill(hit->y_pos - simhit->y_pos);
514  fZres_top->Fill(hit->z_pos - simhit->z_pos);
515  fTres_top->Fill(GetDeltaT(hit,simhit));
516  nmatch++;
517  /*if(nmatch_top<10&&fNViews<100) {
518  //need way to map global back to local coords
519 
520  TVector3 truepos =
521  fCrtutils->WorldToModuleCoords(HitToPosVec(hit),fCrtutils->MacToAuxDetID(hit->feb_id[0],0));
522  TVector3 recopos =
523  fCrtutils->WorldToModuleCoords(HitToPosVec(simhit),fCrtutils->MacToAuxDetID(hit->feb_id[0],0));
524  fTrueHitPoints[fNViews]->SetPoint(0,truepos.X(),truepos.Z());
525  fSimHitPoints[fNViews]->SetPoint(0,recopos.X(),recopos.Z());
526  float pesx[16], pesy[16];
527  for(size_t i=0; i<16; i++){pesx[i]=0.; pesy[i]=0.;}
528  if(simhit->pesmap.find(simhit->feb_id[0])==simhit->pesmap.end())
529  std::cout << "FEB not found in PEs map!" << std::endl;
530  for(auto const& chanpe : (*((simhit->pesmap).find(simhit->feb_id[0]))).second) {
531  if(chanpe.first<16) pesx[chanpe.first] = chanpe.second;
532  else pesy[chanpe.first-16] = chanpe.second;
533  }
534  for(size_t i=0; i<16; i++){
535  for(size_t j=0; j<16; j++){
536  fChargeViews[fNViews]->SetBinContent(i,j,pesx[i]+pesy[j]);
537  }
538  }
539  fNViews++;
540  }*/
541  nmatch_top++;
542  }
543  if(simreg.find("Rim")!=string::npos){
544  fRres_rim->Fill(GetDeltaR(hit,simhit));
545  fXres_rim->Fill(hit->x_pos - simhit->x_pos);
546  fYres_rim->Fill(hit->y_pos - simhit->y_pos);
547  fZres_rim->Fill(hit->z_pos - simhit->z_pos);
548  fTres_rim->Fill(GetDeltaT(hit,simhit));
549  nmatch++;
550  }
551  if(simreg=="Bottom") {
552  fRres_bot->Fill(GetDeltaR(hit,simhit));
553  fXres_bot->Fill(hit->x_pos - simhit->x_pos);
554  fYres_bot->Fill(hit->y_pos - simhit->y_pos);
555  fZres_bot->Fill(hit->z_pos - simhit->z_pos);
556  fTres_bot->Fill(GetDeltaT(hit,simhit));
557  nmatch++;
558  }
559  }
560  else {
561  nmiss++;
562  nmiss_feb++;
563  }//no feb match
564  }//c or d type
565  else {//m type
566  bool match = false;
567  for(auto const& truefeb : hit->feb_id){
568  for(auto const& simfeb : simhit->feb_id){
569  if(truefeb==simfeb) {
570  match = true;
571  break;
572  }
573  }
574  if(match) break;
575  }
576  if(match){
577  if(simreg=="North") {
578  fRres_nor->Fill(GetDeltaR(hit,simhit));
579  fXres_nor->Fill(hit->x_pos - simhit->x_pos);
580  fYres_nor->Fill(hit->y_pos - simhit->y_pos);
581  fZres_nor->Fill(hit->z_pos - simhit->z_pos);
582  fTres_nor->Fill(GetDeltaT(hit,simhit));
583  }
584  else if(simreg=="South") {
585  fRres_sth->Fill(GetDeltaR(hit,simhit));
586  fXres_sth->Fill(hit->x_pos - simhit->x_pos);
587  fYres_sth->Fill(hit->y_pos - simhit->y_pos);
588  fZres_sth->Fill(hit->z_pos - simhit->z_pos);
589  fTres_sth->Fill(GetDeltaT(hit,simhit));
590  }
591  else {
592  fRres_lat->Fill(GetDeltaR(hit,simhit));
593  fXres_lat->Fill(hit->x_pos - simhit->x_pos);
594  fYres_lat->Fill(hit->y_pos - simhit->y_pos);
595  fZres_lat->Fill(hit->z_pos - simhit->z_pos);
596  fTres_lat->Fill(GetDeltaT(hit,simhit));
597  }
598  }
599  else nmiss++;
600  }//else m type
601  }//matching sim hits
602  } //if muID trk found
603  else if(simhasmu){
604  simhasmu = false;
605  nmiss_track++;
606  nmiss++;
607  }
608  }//loop over true hits
609  }//if mu trk found in true map
610  else
611  nmiss_true++;
612 
613  //Sim Hits
614  if(muToSimHits.find(trk)!=muToSimHits.end()) {
615  nsimhit = muToSimHits[trk].size();
616  for(auto const& hit : muToSimHits[trk]) {
617  for(auto const& id : hit->feb_id)
618  fFebsSim->Fill((int)id);
619 
620  if(nregsim.find(hit->tagger)==nregsim.end())
621  nregsim[hit->tagger] = 1;
622  else
623  nregsim[hit->tagger]++;
624 
625  char type = fCrtutils->GetRegTypeFromRegName(hit->tagger);
626  if(type=='c') nsim_top++;
627  if(type=='m') nsim_side++;
628  if(type=='d') nsim_bot++;
629  }
630  }
631  else
632  nmiss_sim++;
633 
634  //Fill histos
635  fNTrueHitMu->Fill(ntruehit);
636  fNSimHitMu->Fill(nsimhit);
637  fNTrueSimDiffMu->Fill(ntruehit-nsimhit);
638 
639  if(ntrue_top>0) fNTrueHitMu_top->Fill(ntrue_top);
640  if(ntrue_side>0) fNTrueHitMu_side->Fill(ntrue_side);
641  if(ntrue_bot>0) fNTrueHitMu_bot->Fill(ntrue_bot);
642  if(nsim_top>0) fNSimHitMu_top->Fill(nsim_top);
643  if(nsim_side>0) fNSimHitMu_side->Fill(nsim_side);
644  if(nsim_bot>0) fNSimHitMu_bot->Fill(nsim_bot);
645 
646  for(auto const& nreg : nregtrue) fNRegTrue->Fill(nreg.second);
647  for(auto const& nreg : nregsim) fNRegSim->Fill(nreg.second);
648 
649  }//for muon IDs
650 
651  if(fVerbose) {
652  std::cout << "found " << '\n'
653  << " * nmu: " << nmu << '\n'
654  << " * ntrue hit: " << trueHitList.size() << " with " << ntrue << " ass'd with muons" << '\n'
655  << " * nsim hit: " << simHitList.size() << " with " << nsim << " ass'd with muons" << '\n'
656  << std::endl;
657 
658  std::cout << "matched true/sim hits:" << '\n'
659  << " * matched: " << nmatch << '\n'
660  << " * missed: " << nmiss << '\n'
661  << " - SimHits missing muon track: " << nmiss_track << '\n'
662  << " - SimHits, c or d missing FEB match: " << nmiss_feb << '\n'
663  << std::endl;
664 
665  std::cout << "missed muons (eff.):" << '\n'
666  << " * trueHit: " << nmiss_true << " (" << 1.0*(nmu-nmiss_true)/nmu << ")" << '\n'
667  << " * data: " << nmiss_data << " (" << 1.0*(nmu-nmiss_data)/nmu << ")" << '\n'
668  << " * simHit: " << nmiss_sim << " (" << 1.0*(nmu-nmiss_sim)/nmu << ")" << '\n'
669  << std::endl;
670  }
671 
672 }
673 
674 float CRTTruthMatchAnalysis::GetDeltaR(art::Ptr<CRTHit> h1, art::Ptr<CRTHit> h2){
675  float dr;
676  double dx = h1->x_pos - h2->x_pos;
677  double dy = h1->y_pos - h2->y_pos;
678  double dz = h1->z_pos - h2->z_pos;
679 
680  dr = sqrt(dx*dx + dy*dy + dz*dz);
681  return dr;
682 }
683 
684 
685 int CRTTruthMatchAnalysis::GetDeltaT(art::Ptr<CRTHit> h1, art::Ptr<CRTHit> h2){
686  int t1 = h1->ts0_ns;
687  int t2 = h2->ts0_ns;
688  return t1-t2;
689 }
690 
691 TVector3 CRTTruthMatchAnalysis::HitToPosVec(art::Ptr<CRTHit> hit){
692 
693  TVector3 vec(hit->x_pos,hit->y_pos,hit->z_pos);
694  return vec;
695 }
696 
697 DEFINE_ART_MODULE(CRTTruthMatchAnalysis)
analyze
process_name opflashCryo1 flashfilter analyze
Definition: run_flashfilter_data.fcl:35
icarus::crt::CRTTruthMatchAnalysis::GetDeltaR
float GetDeltaR(art::Ptr< CRTHit > h1, art::Ptr< CRTHit > h2)
Definition: CRTTruthMatchAnalysis_module.cc:674
icarus::crt::CRTTruthMatchAnalysis::CRTTruthMatchAnalysis
CRTTruthMatchAnalysis(fhicl::ParameterSet const &p)
Definition: CRTTruthMatchAnalysis_module.cc:129
p
pdgs p
Definition: selectors.fcl:22
icarus::ns::util::size
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
hit
process_name hit
Definition: cheaterreco.fcl:51
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
AuxDetSimChannel.h
object containing MC truth information necessary for making RawDigits and doing back tracking ...
icarus::crt::CRTTruthMatchAnalysis::HitToPosVec
TVector3 HitToPosVec(art::Ptr< CRTHit > hit)
Definition: CRTTruthMatchAnalysis_module.cc:691
e
do i e
Definition: updateVersionList.sh:68
icarus::crt::CRTBackTracker::AllTrueIds
std::vector< int > AllTrueIds(const art::Event &event, const CRTData &data)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTBackTracker.cc:208
icarus::crt::CRTTruthMatchAnalysis::analyze
void analyze(art::Event const &e) override
Definition: CRTTruthMatchAnalysis_module.cc:331
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
crt
process_name crt
Definition: run_genericCRT.fcl:82
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
icarus::crt::CRTTruthMatchAnalysis::GetDeltaT
int GetDeltaT(art::Ptr< CRTHit > h1, art::Ptr< CRTHit > h2)
Definition: CRTTruthMatchAnalysis_module.cc:685
Generated by   doxygen 1.8.5