All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CRTPMTMatchingAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: CRTPMTMatchingAna
3 // Plugin Type: analyzer (Unknown Unknown)
4 // File: CRTPMTMatchingAna_module.cc
5 //
6 // Generated at Thu Feb 17 13:31:00 2022 by Biswaranjan Behera using cetskelgen
7 // from version .
8 ////////////////////////////////////////////////////////////////////////
9 
10 #include "art/Framework/Core/EDAnalyzer.h"
11 #include "art/Framework/Core/ModuleMacros.h"
12 #include "art/Framework/Principal/Event.h"
13 #include "art/Framework/Principal/Handle.h"
14 #include "art/Framework/Principal/Run.h"
15 #include "art/Framework/Principal/SubRun.h"
16 #include "canvas/Utilities/InputTag.h"
17 #include "fhiclcpp/ParameterSet.h"
18 #include "messagefacility/MessageLogger/MessageLogger.h"
19 #include "art_root_io/TFileService.h"
20 #include "art/Framework/Services/Registry/ServiceHandle.h"
21 #include "canvas/Persistency/Common/FindManyP.h"
22 
23 //LArSoft includes
28 //#include "larsim/MCCheater/PhotonBackTrackerService.h"
29 
30 //Data product includes
31 #include "nusimdata/SimulationBase/MCGeneratorInfo.h"
32 #include "nusimdata/SimulationBase/MCTruth.h"
33 #include "nusimdata/SimulationBase/MCParticle.h"
41 
42 //C++ includes
43 #include <vector>
44 #include <map>
45 
46 //ROOT includes
47 #include "TTree.h"
48 #include "TVector3.h"
49 
50 using std::vector;
51 using std::map;
52 
53 
54 namespace icarus {
55  namespace crt {
56  class CRTPMTMatchingAna;
57  }
58 }
59 
60 using namespace icarus::crt;
61 
62 class icarus::crt::CRTPMTMatchingAna : public art::EDAnalyzer {
63 public:
64 
66 
67  explicit CRTPMTMatchingAna(fhicl::ParameterSet const& p);
68  // The compiler-generated destructor is fine for non-base
69  // classes without bare pointers or other resource use.
70 
71  // Plugins should not be copied or assigned.
72  CRTPMTMatchingAna(CRTPMTMatchingAna const&) = delete;
74  CRTPMTMatchingAna& operator=(CRTPMTMatchingAna const&) = delete;
75  CRTPMTMatchingAna& operator=(CRTPMTMatchingAna&&) = delete;
76 
77  // Required functions.
78  void analyze(art::Event const& e) override;
79 
80 private:
81 
82  // Declare member data here.
83 
84  bool HitCompare(const art::Ptr<CRTHit>& h1, const art::Ptr<CRTHit>& h2);
85  void ClearVecs();
86 
87  art::InputTag fOpHitModuleLabel;
88  art::InputTag fOpFlashModuleLabel0;
89  art::InputTag fOpFlashModuleLabel1;
90  art::InputTag fOpFlashModuleLabel2;
91  art::InputTag fOpFlashModuleLabel3;
92  art::InputTag fCrtHitModuleLabel;
93  art::InputTag fTriggerLabel;
94  //tart::InputTag fCrtTrackModuleLabel;
95 
96  int fEvent; ///< number of the event being processed
97  int fRun; ///< number of the run being processed
98  int fSubRun; ///< number of the sub-run being processed
99 
100  //add trigger data product vars
101  unsigned int m_gate_type;
102  std::string m_gate_name;
106  uint64_t m_gate_crt_diff;
107 
108  double fCoinWindow;
109  double fOpDelay;
110  double fCrtDelay;
117 
118  // CRTBackTracker* bt;
120 
121  map<int,art::InputTag> fFlashLabels;
122 
123  TTree* fMatchTree;
124 
125  //matchTree vars
126  int fNCrt; //number of CRT hits per event
127  vector<vector<double>> fCrtXYZT; //CRT hit x,y,z,t [cm/ns]
128  vector<vector<double>> fCrtXYZErr; //CRT hit x,y,z,t error [cm/ns]
129  vector<int> fCrtRegion; //CRT hit region code
130  vector<double> fCrtPE; //total PE's for CRT hit
131  vector<double> fTofHit; //CRT - PMT TOF using OpHit
132  vector<double> fTofFlash; //CRT - PMT TOF using OpFlash
133  vector<double> fTofFlashHit;
134  vector<double> fDistHit; //CRT - PMT distance [cm]
135  vector<double> fDistFlash; //CRT - flash barycenter distance [cm]
136  vector<double> fDistFlashHit;
137  vector<double> fTofPeHit; //total PE for matched OpHit
138  vector<double> fTofPeFlash; //total PE for matched OpFlash
139  vector<double> fTofPeFlashHit;
140 
141  vector<vector<double>> fTofXYZTHit; //position/time [cm/ns] for matched OpHit
142  vector<vector<double>> fTofXYZTFlash; //position/time [cm/ns] for matched OpFlash
143  vector<vector<double>> fTofXYZTFlashHit;
144  vector<int> fTofTpcHit; //TPC for matched OpHit
145  vector<int> fTofTpcFlash; //TPC for matched OpFlash
146  vector<bool> fMatchHit; //was CRT hit matched to OpHit?
147  vector<bool> fMatchFlash; //was CRT hit matched to OpFlash?
148  vector<bool> fTrackFilt; //excluded by track filter
149 
150  geo::GeometryCore const* fGeometryService; ///< pointer to Geometry provider
151 };
152 
153 
155  :EDAnalyzer{p} // ,
156  ,fOpHitModuleLabel(p.get<art::InputTag>("OpHitModuleLabel","ophit"))
157  ,fOpFlashModuleLabel0(p.get<art::InputTag>("OpFlashModuleLabel0","opflashTPC0"))
158  ,fOpFlashModuleLabel1(p.get<art::InputTag>("OpFlashModuleLabel1","opflashTPC1"))
159  // ,fOpFlashModuleLabel2(p.get<art::InputTag>("OpFlashModuleLabel2","opflashTPC2"))
160  //,fOpFlashModuleLabel3(p.get<art::InputTag>("OpFlashModuleLabel3","opflashTPC3"))
161  ,fCrtHitModuleLabel(p.get<art::InputTag>("CrtHitModuleLabel","crthit"))
162  ,fTriggerLabel(p.get<art::InputTag>("TriggerLabel","daqTrigger"))
163  // ,fCrtTrackModuleLabel(p.get<art::InputTag>("CrtTrackModuleLabel","crttrack"))
164  ,fCoinWindow(p.get<double>("CoincidenceWindow",60.0))
165  ,fOpDelay(p.get<double>("OpDelay",55.1))
166  ,fCrtDelay(p.get<double>("CrtDelay",1.6e6))
167  ,fFlashPeThresh(p.get<int>("FlashPeThresh",9000))
168  ,fHitPeThresh(p.get<int>("HitPeThresh",700))
169  ,fFlashVelocity(p.get<double>("FlashVelocityThresh",-40.))
170  ,fFlashZOffset(p.get<double>("FlashZOffset",0.))
171  ,fHitVelocityMax(p.get<double>("HitVelocityMax",20.))
172  ,fHitVelocityMin(p.get<double>("HitVelocityMin",1.))
173  // bt(new CRTBackTracker(p.get<fhicl::ParameterSet>("CRTBackTrack"))),
174  ,crtutil(new CRTCommonUtils())
175  // More initializers here.
176 {
177  // Call appropriate consumes<>() for any products to be retrieved by this module.
178  fFlashLabels[0] = fOpFlashModuleLabel0;
179  fFlashLabels[1] = fOpFlashModuleLabel1;
180  fFlashLabels[2] = fOpFlashModuleLabel2;
181  fFlashLabels[3] = fOpFlashModuleLabel3;
182 
183  // Get a pointer to the geometry service provider.
184  fGeometryService = lar::providerFrom<geo::Geometry>();
185 
186  art::ServiceHandle<art::TFileService> tfs;
187 
188  fMatchTree = tfs->make<TTree>("matchTree","CRTHit - OpHit/Flash matching analysis");
189 
190  fMatchTree->Branch("event", &fEvent, "event/I");
191  fMatchTree->Branch("run", &fRun, "run/I");
192  fMatchTree->Branch("subrun", &fSubRun, "subrun/I");
193  fMatchTree->Branch("ncrt", &fNCrt, "nCrtHit/I");
194  fMatchTree->Branch("crtXYZT", &fCrtXYZT);
195  fMatchTree->Branch("crtXYZErr", &fCrtXYZErr);
196  fMatchTree->Branch("crtPE", &fCrtPE);
197  fMatchTree->Branch("crtRegion", &fCrtRegion);
198  fMatchTree->Branch("tofHit", &fTofHit);
199  fMatchTree->Branch("tofFlash", &fTofFlash);
200  fMatchTree->Branch("tofFlashHit", &fTofFlashHit);
201  fMatchTree->Branch("peHit", &fTofPeHit);
202  fMatchTree->Branch("peFlash", &fTofPeFlash);
203  fMatchTree->Branch("peFlashHit", &fTofPeFlashHit);
204  fMatchTree->Branch("xyztHit", &fTofXYZTHit);
205  fMatchTree->Branch("xyztFlash", &fTofXYZTFlash);
206  fMatchTree->Branch("xyztFlashHit", &fTofXYZTFlashHit);
207  fMatchTree->Branch("distHit", &fDistHit);
208  fMatchTree->Branch("distFlash", &fDistFlash);
209  fMatchTree->Branch("distFlashHit", &fDistFlashHit);
210  fMatchTree->Branch("tpcHit", &fTofTpcHit);
211  fMatchTree->Branch("tpcFlash", &fTofTpcFlash);
212  fMatchTree->Branch("gate_type", &m_gate_type, "gate_type/b");
213  fMatchTree->Branch("gate_name", &m_gate_name);
214  fMatchTree->Branch("trigger_timestamp", &m_trigger_timestamp, "trigger_timestamp/l");
215  fMatchTree->Branch("gate_start_timestamp", &m_gate_start_timestamp, "gate_start_timestamp/l");
216  fMatchTree->Branch("trigger_gate_diff", &m_trigger_gate_diff, "trigger_gate_diff/l");
217  fMatchTree->Branch("gate_crt_diff",&m_gate_crt_diff, "gate_crt_diff/l");
218 }
219 
221 {
222  // Implementation of required member function here.
223  /*
224  geo::CryostatGeo const& cryo0 = fGeometryService->Cryostat(0);
225  geo::CryostatGeo const& cryo1 = fGeometryService->Cryostat(1);
226  geo::TPCGeo const& tpc00 = cryo0.TPC(0);
227  geo::TPCGeo const& tpc01 = cryo0.TPC(1);
228  geo::TPCGeo const& tpc10 = cryo1.TPC(0);
229  geo::TPCGeo const& tpc11 = cryo1.TPC(1);
230  */
231  MF_LOG_DEBUG("CRTPMTMatching: ") << "beginning analyis" << '\n';
232 
233  // Start by fetching some basic event information for our n-tuple.
234  fEvent = e.id().event();
235  fRun = e.run();
236  fSubRun = e.subRun();
237 
238  ClearVecs();
239 
240  //add trigger info
241  if( !fTriggerLabel.empty() ) {
242 
243  art::Handle<sbn::ExtraTriggerInfo> trigger_handle;
244  e.getByLabel( fTriggerLabel, trigger_handle );
245  if( trigger_handle.isValid() ) {
246  sbn::triggerSource bit = trigger_handle->sourceType;
247  m_gate_type = (unsigned int)bit;
248  m_gate_name = bitName(bit);
249  m_trigger_timestamp = trigger_handle->triggerTimestamp;
250  m_gate_start_timestamp = trigger_handle->beamGateTimestamp;
251  m_trigger_gate_diff = trigger_handle->triggerTimestamp - trigger_handle->beamGateTimestamp;
252 
253  }
254  else{
255  mf::LogError("CRTPMTMatching:") << "No raw::Trigger associated to label: " << fTriggerLabel.label() << "\n" ;
256  }
257  }
258  else {
259  std::cout << "Trigger Data product " << fTriggerLabel.label() << " not found!\n" ;
260  }
261 
262 
263  //OpHits
264  art::Handle< std::vector<recob::OpHit> > opHitListHandle;
265  std::vector< art::Ptr<recob::OpHit> > opHitList;
266  if( e.getByLabel(fOpHitModuleLabel,opHitListHandle) )
267  art::fill_ptr_vector(opHitList, opHitListHandle);
268 
269 
270  //OpFlash
271  map<int, art::Handle< std::vector<recob::OpFlash> > > flashHandles;
272  std::map<int,std::vector< art::Ptr<recob::OpFlash> >> opFlashLists;
273  // fNFlash = 0;
274  for(int i=0; i<2; i++) {
275  if( e.getByLabel(fFlashLabels[i],flashHandles[i]) )
276  art::fill_ptr_vector(opFlashLists[i], flashHandles[i]);
277  }
278 
279 
280  //CRTHits
281  art::Handle< std::vector<CRTHit> > crtHitListHandle;
282  std::vector< art::Ptr<CRTHit> > crtHitList;
283  if( e.getByLabel(fCrtHitModuleLabel,crtHitListHandle))
284  art::fill_ptr_vector(crtHitList, crtHitListHandle);
285 
286  fNCrt = crtHitList.size();
287 
288  for(auto const& crt : crtHitList){
289  vector<double> xyzt, xyzerr;
290  TVector3 rcrt(crt->x_pos,crt->y_pos,crt->z_pos);
291  m_gate_crt_diff = m_gate_start_timestamp - crt->ts0_ns;
292  xyzt.push_back(rcrt.X());
293  xyzt.push_back(rcrt.Y());
294  xyzt.push_back(rcrt.Z());
295 
296  double tcrt = double(m_gate_start_timestamp - crt->ts0_ns)/1e3;
297  tcrt = -tcrt+1e6;
298 
299  // double tcrt = (int32_t)crt->ts0_ns;// - fCrtDelay;
300  //uint64_t tcrt = crt->ts0_ns;
301 
302  // std::cout <<"T0: " << crt->ts0_ns << " , tcrt : " << tcrt << std::endl;
303  xyzt.push_back(tcrt);
304  fCrtXYZT.push_back(xyzt);
305 
306  xyzerr.push_back(crt->x_err);
307  xyzerr.push_back(crt->y_err);
308  xyzerr.push_back(crt->z_err);
309  fCrtXYZErr.push_back(xyzerr);
310 
311  fCrtPE.push_back(crt->peshit);
312  fCrtRegion.push_back(crtutil->AuxDetRegionNameToNum(crt->tagger));
313 
314  // -- flash match --
315  int matchtpc = -1;
316  double tdiff = DBL_MAX, rdiff=DBL_MAX, peflash=DBL_MAX;
317  bool matched=false;
318  xyzt.clear();
319  double flashHitT = DBL_MAX, flashHitPE=DBL_MAX, flashHitDiff=DBL_MAX;
320  vector<double> flashHitxyzt;
321 
322  for(auto const& flashList : opFlashLists) {
323 
324  art::FindManyP<recob::OpHit> findManyHits(flashHandles[flashList.first], e, fFlashLabels[flashList.first]);
325 
326  for(size_t iflash=0; iflash<flashList.second.size(); iflash++) {
327 
328  auto const& flash = flashList.second[iflash];
329  if(flash->TotalPE()<fFlashPeThresh){
330  continue;
331  }
332 
333  double tflash = flash->Time();//*1e3;//-fOpDelay;
334 
335  TVector3 rflash(0,flash->YCenter(),flash->ZCenter());
336  TVector3 vdiff = rcrt-rflash;
337  //std::cout << "flash time: "<<flash->Time() << " absdiff : "<< abs(tcrt-tflash)<< std::endl;
338  if(abs(tcrt-tflash)<abs(tdiff)) {
339  peflash = flash->TotalPE();
340  tdiff = tcrt-tflash;
341  rdiff = vdiff.Mag();
342  xyzt.clear();
343  xyzt.push_back(rflash.X());
344  xyzt.push_back(rflash.Y());
345  xyzt.push_back(rflash.Z());
346  xyzt.push_back(tflash);
347  matched = true;
348  matchtpc = flashList.first;
349 
350  vector<art::Ptr<recob::OpHit>> hits = findManyHits.at(iflash);
351  for(auto const& hit : hits) {
352  double tPmt = hit->PeakTime();//*1.e3;//-fOpDelay;
353  if( tPmt < flashHitT) {
354  flashHitT = tPmt;
355  flashHitPE = hit->PE();
356  //FlashHit position/time
357  double pos[3];
358  fGeometryService->OpDetGeoFromOpChannel(hit->OpChannel()).GetCenter(pos);
359  flashHitxyzt.clear();
360  for(int i=0; i<3; i++) flashHitxyzt.push_back(pos[i]);
361  flashHitxyzt.push_back(flashHitT);
362 
363  //FlashHit distance
364  TVector3 rflashHit(pos[0],pos[1],pos[2]);
365  TVector3 vdiffHit = rcrt-rflashHit;
366  flashHitDiff = vdiffHit.Mag();
367  }
368  }//loop over flash hits
369  }//if minimum tdiff
370  }//for OpFlash in this flash list
371  }//for flash lists
372  if(!matched) {
373  peflash = DBL_MAX;
374  for(int i=0; i<4; i++) xyzt.push_back(DBL_MAX);
375  }
376 
377  fMatchFlash.push_back(matched);
378  fTofFlash.push_back(tdiff);
379  fTofPeFlash.push_back(peflash);
380  fTofXYZTFlash.push_back(xyzt);
381  fDistFlash.push_back(rdiff);
382  fTofTpcFlash.push_back(matchtpc);
383  fTofFlashHit.push_back(tcrt-flashHitT);
384  fTofPeFlashHit.push_back(flashHitPE);
385  fTofXYZTFlashHit.push_back(flashHitxyzt);
386  fDistFlashHit.push_back(flashHitDiff);
387 
388  // -- match OpHits to CRTHits --
389  tdiff = DBL_MAX;
390  matched = false;
391  peflash = DBL_MAX;
392  rdiff = DBL_MAX;
393  xyzt.clear();
394 
395  double pemax = 0.;
396  for(auto const& hit : opHitList) {
397  double thit = hit->PeakTime();//*1e3-fOpDelay;
398  if(hit->PE()<fHitPeThresh){
399  continue;
400  }
401 
402 
403  //if(abs(tcrt-thit)<abs(tdiff)) {
404  if(abs(tcrt-thit)<fCoinWindow && hit->PE()>pemax) {
405  pemax = hit->PE();
406 
407  //hitXYZT
408  double pos[3];
409  fGeometryService->OpDetGeoFromOpChannel(hit->OpChannel()).GetCenter(pos);
410 
411  //distHit
412  TVector3 rhit (pos[0],pos[1],pos[2]);
413  TVector3 vdiff = rcrt-rhit;
414  rdiff = vdiff.Mag();
415  peflash = hit->PE();
416  tdiff = tcrt-thit;
417  xyzt.clear();
418  for(int i=0; i<3; i++) xyzt.push_back(pos[i]);
419  xyzt.push_back(thit);
420 
421  matched = true;
422  std::cout << "thit: "<<thit << " , tcrt: " << tcrt << " , tdiff " << tdiff << std::endl;
423  }//if min tof
424  }//for OpHits
425  if(!matched) {
426  tdiff = DBL_MAX;
427  peflash = DBL_MAX;
428  rdiff = DBL_MAX;
429  xyzt.clear();
430  for(int i=0; i<4; i++) xyzt.push_back(DBL_MAX);
431  }
432  fMatchHit.push_back(matched);
433  fTofHit.push_back(tdiff);
434  fTofPeHit.push_back(peflash);
435  fTofXYZTHit.push_back(xyzt);
436  fDistHit.push_back(rdiff);
437  fTofTpcHit.push_back(matchtpc);
438  }//for CRTHits
439 
440  fMatchTree->Fill();
441 }
442 
444 {
445  //matchTree
446  fCrtXYZT.clear();
447  fCrtXYZErr.clear();
448  fCrtRegion.clear();
449  fCrtPE.clear();
450  fTofHit.clear();
451  fTofFlash.clear();
452  fTofFlashHit.clear();
453  fDistHit.clear();
454  fDistFlash.clear();
455  fDistFlashHit.clear();
456  fTofPeHit.clear();
457  fTofPeFlash.clear();
458  fTofPeFlashHit.clear();
459  fTofXYZTHit.clear();
460  fTofXYZTFlash.clear();
461  fTofXYZTFlashHit.clear();
462  fTofTpcHit.clear();
463  fTofTpcFlash.clear();
464  fMatchHit.clear();
465  fMatchFlash.clear();
466 }
467 
468 bool icarus::crt::CRTPMTMatchingAna::HitCompare(const art::Ptr<CRTHit>& hit1, const art::Ptr<CRTHit>& hit2) {
469 
470  if(hit1->ts1_ns != hit2->ts1_ns) return false;
471  if(hit1->plane != hit2->plane) return false;
472  if(hit1->x_pos != hit2->x_pos) return false;
473  if(hit1->y_pos != hit2->y_pos) return false;
474  if(hit1->z_pos != hit2->z_pos) return false;
475  if(hit1->x_err != hit2->x_err) return false;
476  if(hit1->y_err != hit2->y_err) return false;
477  if(hit1->z_err != hit2->z_err) return false;
478  if(hit1->tagger != hit2->tagger) return false;
479 
480  return true;
481 }
482 
483 DEFINE_ART_MODULE(icarus::crt::CRTPMTMatchingAna)
vector< vector< double > > fTofXYZTFlashHit
process_name opflashCryo1 flashfilter analyze
int fEvent
number of the event being processed
Utilities related to art service access.
pdgs p
Definition: selectors.fcl:22
map< int, art::InputTag > fFlashLabels
vector< vector< double > > fTofXYZTHit
process_name hit
Definition: cheaterreco.fcl:51
std::string bitName(triggerSource bit)
Returns a mnemonic short name of the beam type.
Definition: BeamBits.h:267
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
Access the description of detector geometry.
T abs(T value)
vector< vector< double > > fTofXYZTFlash
Description of geometry of one entire detector.
bool HitCompare(const art::Ptr< CRTHit > &h1, const art::Ptr< CRTHit > &h2)
vector< vector< double > > fCrtXYZErr
int fSubRun
number of the sub-run being processed
CRTPMTMatchingAna(fhicl::ParameterSet const &p)
int fRun
number of the run being processed
do i e
static double tdiff(const art::Timestamp &ts1, const art::Timestamp &ts2)
triggerSource
Type of beam or beam gate or other trigger source.
Definition: BeamBits.h:97
Data product holding additional trigger information.
art::ServiceHandle< art::TFileService > tfs
geo::GeometryCore const * fGeometryService
pointer to Geometry provider
process_name crt
void analyze(art::Event const &e) override
art framework interface to geometry description
BEGIN_PROLOG could also be cout