All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ToFProducer_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: ToFProducer
3 // Plugin Type: producer (Unknown Unknown)
4 // File: ToFProducer_module.cc
5 //
6 // Generated at Mon Jan 10 18:44:13 2022 by Varuna Crishan Meddage using cetskelgen
7 // from version .
8 ////////////////////////////////////////////////////////////////////////
9 
10 #include "art/Framework/Core/ModuleMacros.h"
11 #include "art/Framework/Core/EDProducer.h"
12 #include "art/Framework/Principal/Event.h"
13 #include "art/Framework/Principal/SubRun.h"
14 #include "art/Framework/Principal/Handle.h"
15 #include "art/Framework/Principal/View.h"
16 #include "canvas/Persistency/Common/Ptr.h"
17 #include "canvas/Persistency/Common/PtrVector.h"
18 #include "art/Framework/Services/Registry/ServiceHandle.h"
19 #include "art_root_io/TFileService.h"
20 #include "art_root_io/TFileDirectory.h"
21 #include "canvas/Persistency/Common/FindMany.h"
22 #include "canvas/Persistency/Common/FindManyP.h"
23 #include "canvas/Persistency/Common/FindOneP.h"
24 #include "fhiclcpp/ParameterSet.h"
25 #include "messagefacility/MessageLogger/MessageLogger.h"
26 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
27 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
28 #include "lardataalg/DetectorInfo/DetectorPropertiesData.h"
29 
30 #include "larcore/Geometry/Geometry.h"
31 #include "nusimdata/SimulationBase/GTruth.h"
32 #include "nusimdata/SimulationBase/MCTruth.h"
33 #include "nusimdata/SimulationBase/MCFlux.h"
34 #include "lardataobj/Simulation/SimChannel.h"
35 #include "lardataobj/Simulation/AuxDetSimChannel.h"
36 #include "lardataobj/AnalysisBase/Calorimetry.h"
37 #include "lardataobj/AnalysisBase/ParticleID.h"
38 #include "lardataobj/RawData/RawDigit.h"
39 #include "lardataobj/RawData/BeamInfo.h"
40 #include "larcoreobj/SummaryData/POTSummary.h"
41 #include "larsim/MCCheater/BackTrackerService.h"
42 #include "larsim/MCCheater/ParticleInventoryService.h"
43 #include "lardataobj/RecoBase/Track.h"
44 #include "lardataobj/RecoBase/Slice.h"
45 #include "lardataobj/RecoBase/Shower.h"
46 #include "lardataobj/RecoBase/PFParticle.h"
47 #include "lardataobj/RecoBase/Cluster.h"
48 #include "lardataobj/RecoBase/Hit.h"
49 #include "lardataobj/RecoBase/EndPoint2D.h"
50 #include "lardataobj/RecoBase/Vertex.h"
51 #include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
52 #include "lardataobj/AnalysisBase/CosmicTag.h"
53 #include "lardataobj/AnalysisBase/FlashMatch.h"
54 #include "lardataobj/RecoBase/MCSFitResult.h"
55 #include "larreco/RecoAlg/TrackMomentumCalculator.h"
56 #include "larreco/RecoAlg/TrajectoryMCSFitter.h"
57 #include "lardataobj/RecoBase/MCSFitResult.h"
58 #include "larreco/RecoAlg/PMAlg/PmaTrack3D.h"
59 #include "larpandora/LArPandoraInterface/LArPandoraHelper.h"
60 #include "lardataobj/RecoBase/SpacePoint.h"
61 #include "lardataobj/RecoBase/OpFlash.h"
62 #include "lardataobj/AnalysisBase/T0.h"
63 
64 #include "sbndcode/RecoUtils/RecoUtils.h"
65 #include "lardataobj/MCBase/MCShower.h"
66 
67 // sbndcode includes
68 #include "sbnobj/Common/CRT/CRTHit.hh"
69 #include "sbnobj/Common/CRT/CRTTrack.hh"
70 #include "sbndcode/CRT/CRTUtils/CRTT0MatchAlg.h"
71 #include "sbndcode/CRT/CRTUtils/CRTTrackMatchAlg.h"
72 #include "sbndcode/CRT/CRTUtils/CRTCommonUtils.h"
73 #include "sbndcode/CRT/CRTUtils/CRTBackTracker.h"
74 
75 #include "lardataobj/RecoBase/OpHit.h"
76 
77 #include "sbnobj/SBND/ToF/ToF.hh"
78 
79 #include "TFile.h"
80 #include "TTree.h"
81 #include "TDirectory.h"
82 #include "TH1.h"
83 #include "TH2.h"
84 #include "TF1.h"
85 #include "TProfile.h"
86 #include "TROOT.h"
87 #include "TStyle.h"
88 #include "TMath.h"
89 #include "TGraphErrors.h"
90 #include "TMinuit.h"
91 #include "TString.h"
92 #include "TTimeStamp.h"
93 #include "TVectorD.h"
94 #include "TCanvas.h"
95 #include "TFrame.h"
96 #include "TLine.h"
97 #include "TAxis.h"
98 #include "TTimeStamp.h"
99 
100 #include <vector>
101 #include <fstream>
102 #include "TPaveStats.h"
103 #include <iostream>
104 #include <string>
105 #include "math.h"
106 #include "stdio.h"
107 #include <iterator>
108 #include <map>
109 #include <memory>
110 
111 using namespace std;
112 
113 namespace sbnd{
114 
115 class ToFProducer : public art::EDProducer {
116 public:
117  using CRTHit = sbn::crt::CRTHit;
118  explicit ToFProducer(fhicl::ParameterSet const& pset);
119 
120  ToFProducer(ToFProducer const&) = delete;
121  ToFProducer(ToFProducer&&) = delete;
122  ToFProducer& operator=(ToFProducer const&) = delete;
123  ToFProducer& operator=(ToFProducer&&) = delete;
124 
125  void beginJob() override;
126  void produce(art::Event& evt) override;
127  bool HitCompare(const art::Ptr<CRTHit>& h1, const art::Ptr<CRTHit>& h2);
128 
129 
130 private:
131  art::InputTag fOpHitModuleLabel;
132  art::InputTag fOpFlashModuleLabel0;
133  art::InputTag fOpFlashModuleLabel1;
134  art::InputTag fCrtHitModuleLabel;
135  art::InputTag fCrtTrackModuleLabel;
136 
137  double fCoinWindow;
138  double fOpDelay;
139  double fCRThitThresh;
140  double fFlashPeThresh;
141  double fHitPeThresh;
142  double fBeamLow;
143  double fBeamUp;
144  bool fLFlash;
145  bool fLFlash_hit;
146  bool fCFlash;
147  bool fCFlash_hit;
148  bool fLhit;
149  bool fChit;
150 
151  CRTBackTracker* bt;
152  map<int,art::InputTag> fFlashLabels;
153  geo::GeometryCore const* fGeometryService;
154 };
155 
156 ToFProducer::ToFProducer(fhicl::ParameterSet const& pset):
157 EDProducer{pset},
158 fOpHitModuleLabel(pset.get<art::InputTag>("OpHitModuleLabel")),
159 fOpFlashModuleLabel0(pset.get<art::InputTag>("OpFlashModuleLabel0")),
160 fOpFlashModuleLabel1(pset.get<art::InputTag>("OpFlashModuleLabel1")),
161 fCrtHitModuleLabel(pset.get<art::InputTag>("CrtHitModuleLabel")),
162 fCrtTrackModuleLabel(pset.get<art::InputTag>("CrtTrackModuleLabel")),
163 fCoinWindow(pset.get<double>("CoincidenceWindow")),
164 fOpDelay(pset.get<double>("OpDelay")),
165 fCRThitThresh(pset.get<double>("CRThitThresh")),
166 fFlashPeThresh(pset.get<double>("FlashPeThresh")),
167 fHitPeThresh(pset.get<double>("HitPeThresh")),
168 fBeamLow(pset.get<double>("BeamLow")),
169 fBeamUp(pset.get<double>("BeamUp")),
170 fLFlash(pset.get<bool>("LFlash")),
171 fLFlash_hit(pset.get<bool>("LFlash_hit")),
172 fCFlash(pset.get<bool>("CFlash")),
173 fCFlash_hit(pset.get<bool>("CFlash_hit")),
174 fLhit(pset.get<bool>("Lhit")),
175 fChit(pset.get<bool>("Chit")),
176 bt(new CRTBackTracker(pset.get<fhicl::ParameterSet>("CRTBackTrack")))
177 {
178  fFlashLabels[0] = fOpFlashModuleLabel0;
179  fFlashLabels[1] = fOpFlashModuleLabel1;
180  produces< std::vector<sbnd::ToF::ToF> >();
181 }
182 
183 void ToFProducer::beginJob(){
184  std::cout<<"job begin..."<<std::endl;
185  fGeometryService = lar::providerFrom<geo::Geometry>();
186 }
187 
188 void ToFProducer::produce(art::Event& evt)
189 {
190  std::unique_ptr<vector<sbnd::ToF::ToF>> ToF_vec(new vector<sbnd::ToF::ToF>);
191 
192  art::Handle< std::vector<sbn::crt::CRTTrack> > crtTrackListHandle;
193  std::vector< art::Ptr<sbn::crt::CRTTrack> > crtTrackList;
194  if( evt.getByLabel(fCrtTrackModuleLabel,crtTrackListHandle))
195  art::fill_ptr_vector(crtTrackList, crtTrackListHandle);
196 
197  art::FindManyP<sbn::crt::CRTHit> findManyHits(crtTrackListHandle, evt, fCrtTrackModuleLabel);
198  std::vector<std::vector<art::Ptr<CRTHit>>> trackhits;
199 
200  //================================================================
201 
202  for(size_t itrk=0; itrk<crtTrackList.size(); itrk++){
203  std::vector<art::Ptr<CRTHit>> trkhits = findManyHits.at(itrk);
204  std::sort(trkhits.begin(),trkhits.end(),
205  [](const art::Ptr<CRTHit>& a, const art::Ptr<CRTHit>& b)->bool
206  {
207  return a->ts1_ns < b->ts1_ns;
208  });
209  trackhits.push_back(trkhits);
210  } // Crt hits coming from are ordered on ascending order by looking into ts1_ns variable
211 
212  //==================================================================
213 
214  art::Handle< std::vector<sbn::crt::CRTHit> > crtHitListHandle;
215  std::vector< art::Ptr<sbn::crt::CRTHit> > crtHitList;
216  if( evt.getByLabel(fCrtHitModuleLabel,crtHitListHandle))
217  art::fill_ptr_vector(crtHitList, crtHitListHandle);
218 
219  map<int, std::vector<art::Ptr<CRTHit>> > tof_crt_hits;
220  map<int, std::vector<art::Ptr<recob::OpHit>> > tof_op_hits;
221  map<int, std::vector<art::Ptr<recob::OpFlash>> > tof_op_flashes;
222  map<int, std::vector<int>> tof_op_tpc;
223 
224  for(auto const& crt : crtHitList){
225  if(!(crt->ts1_ns >= fBeamLow && crt->ts1_ns<= fBeamUp)) continue;
226  if(crt->peshit < fCRThitThresh) continue;
227 
228  bool frm_trk=false;
229  int index=0;
230 
231  for(auto const& trkhits: trackhits){
232  for(size_t ihit=0; ihit<trkhits.size(); ihit++){
233  if(HitCompare(trkhits[ihit],crt)){
234  frm_trk=true;
235  break;
236  }
237  }
238  if(frm_trk) break;
239  index++;
240  }
241 
242  // ================================== Calculatin ToF values using Largest optical hit method =========================
243 
244  if(fLhit){
245  double pehit_max=0;
246  bool found_tof = false;
247  int ophit_index = -1;
248 
249  art::Handle< std::vector<recob::OpHit> > opHitListHandle;
250  std::vector< art::Ptr<recob::OpHit> > opHitList;
251  if( evt.getByLabel(fOpHitModuleLabel,opHitListHandle) )
252  art::fill_ptr_vector(opHitList, opHitListHandle);
253 
254  for(auto const& hit : opHitList){
255  if(hit->PE()<fHitPeThresh) continue;
256  double thit = hit->PeakTime()*1e3-fOpDelay;
257 
258  if(abs(crt->ts1_ns-thit)<fCoinWindow && hit->PE()>pehit_max){
259  pehit_max = hit->PE();
260  ophit_index = hit.key();
261  found_tof = true;
262  }
263  } // loop over optical hit list
264 
265  if(found_tof){
266  if(frm_trk){
267  tof_crt_hits[index].push_back(crt);
268  tof_op_hits[index].push_back(opHitList[ophit_index]);
269  } // crt hit is coming from crt track
270 
271  else{
272  sbnd::ToF::ToF new_tof;
273  new_tof.tof = crt->ts1_ns - (opHitList[ophit_index]->PeakTime()*1e3-fOpDelay);
274  new_tof.frm_hit = true;
275  new_tof.crt_time = crt->ts1_ns;
276  new_tof.pmt_time = opHitList[ophit_index]->PeakTime()*1e3-fOpDelay;
277  new_tof.crt_tagger = crt->tagger;
278  new_tof.crt_hit_id = crt.key();
279  new_tof.pmt_hit_id = ophit_index;
280  ToF_vec->push_back(new_tof);
281  } // lonely crt hit
282  } // found a tof match
283 
284  } // use Lhit method
285 
286  // =====================================================================================================================
287 
288  // ================================== Calculatin ToF values using Closest optical hit method ===========================
289 
290  if(fChit){
291  double ophit_minTOF = DBL_MAX;
292  bool found_tof = false;
293  int ophit_index = -1;
294 
295  art::Handle< std::vector<recob::OpHit> > opHitListHandle;
296  std::vector< art::Ptr<recob::OpHit> > opHitList;
297  if( evt.getByLabel(fOpHitModuleLabel,opHitListHandle) )
298  art::fill_ptr_vector(opHitList, opHitListHandle);
299 
300  for(auto const& hit : opHitList){
301  if(hit->PE()<fHitPeThresh) continue;
302  double thit = hit->PeakTime()*1e3-fOpDelay;
303 
304  if(abs(crt->ts1_ns-thit)<fCoinWindow && abs(crt->ts1_ns-thit)<ophit_minTOF){
305  ophit_minTOF = abs(crt->ts1_ns-thit);
306  ophit_index = hit.key();
307  found_tof = true;
308  }
309  } // loop over optical hit list
310 
311  if(found_tof){
312  if(frm_trk){
313  tof_crt_hits[index].push_back(crt);
314  tof_op_hits[index].push_back(opHitList[ophit_index]);
315  } // crt hit is coming from crt track
316 
317  else{
318  sbnd::ToF::ToF new_tof;
319  new_tof.tof = crt->ts1_ns - (opHitList[ophit_index]->PeakTime()*1e3-fOpDelay);
320  new_tof.frm_hit = true;
321  new_tof.crt_time = crt->ts1_ns;
322  new_tof.pmt_time = opHitList[ophit_index]->PeakTime()*1e3-fOpDelay;
323  new_tof.crt_tagger = crt->tagger;
324  new_tof.crt_hit_id = crt.key();
325  new_tof.pmt_hit_id = ophit_index;
326  ToF_vec->push_back(new_tof);
327  } // lonely crt hit
328  } // found a tof match
329 
330  } // use Chit method
331 
332  // =====================================================================================================================
333 
334  //==================================== Calculation ToF values using Largest optical flash ==============================
335 
336  if(fLFlash){
337  double peflash_max=0;
338  bool found_tof = false;
339  int opflash_index = -1;
340  int flash_tpc = -1;
341 
342  std::map<int, art::Handle< std::vector<recob::OpFlash> > > flashHandles;
343  std::map<int,std::vector< art::Ptr<recob::OpFlash> >> opFlashLists;
344 
345  for(int i=0; i<2; i++) {
346  if( evt.getByLabel(fFlashLabels[i],flashHandles[i]) )
347  art::fill_ptr_vector(opFlashLists[i], flashHandles[i]);
348  }
349 
350  for(auto const& flashList : opFlashLists){
351  for(size_t iflash=0; iflash<flashList.second.size(); iflash++){
352  auto const& flash = flashList.second[iflash];
353  if(flash->TotalPE()<fFlashPeThresh) continue;
354  double tflash = flash->Time()*1e3-fOpDelay;
355  if(abs(crt->ts1_ns-tflash)<fCoinWindow && flash->TotalPE()>peflash_max){
356  peflash_max=flash->TotalPE();
357  opflash_index = flash.key();
358  found_tof = true;
359  flash_tpc = flashList.first;
360  } // with in conincidence window and getting the largest flash
361  } // loop over flash list
362  } // loop over flash list map
363 
364  if(found_tof){
365  if(frm_trk){
366  tof_crt_hits[index].push_back(crt);
367  tof_op_flashes[index].push_back(opFlashLists[flash_tpc][opflash_index]);
368  tof_op_tpc[index].push_back(flash_tpc);
369  } // crt hit is coming from crt track
370 
371  else{
372  sbnd::ToF::ToF new_tof;
373  new_tof.tof = crt->ts1_ns - (opFlashLists[flash_tpc][opflash_index]->Time()*1e3-fOpDelay);
374  new_tof.frm_hit = true;
375  new_tof.crt_time = crt->ts1_ns;
376  new_tof.pmt_time = opFlashLists[flash_tpc][opflash_index]->Time()*1e3-fOpDelay;
377  new_tof.crt_tagger = crt->tagger;
378  new_tof.crt_hit_id = crt.key();
379  new_tof.flash_tpc_id = flash_tpc;
380  new_tof.pmt_flash_id = opflash_index;
381  ToF_vec->push_back(new_tof);
382  } // lonely crt hit
383  } // found a tof match
384 
385  } // use LFlash method
386 
387  //======================================================================================================================
388 
389  //==================================== Calculation ToF values using Closest optical flash ==============================
390 
391  if(fCFlash){
392  double flash_minTOF = DBL_MAX;
393  bool found_tof = false;
394  int opflash_index = -1;
395  int flash_tpc = -1;
396 
397  std::map<int, art::Handle< std::vector<recob::OpFlash> > > flashHandles;
398  std::map<int,std::vector< art::Ptr<recob::OpFlash> >> opFlashLists;
399 
400  for(int i=0; i<2; i++) {
401  if( evt.getByLabel(fFlashLabels[i],flashHandles[i]) )
402  art::fill_ptr_vector(opFlashLists[i], flashHandles[i]);
403  }
404 
405  for(auto const& flashList : opFlashLists){
406  for(size_t iflash=0; iflash<flashList.second.size(); iflash++){
407  auto const& flash = flashList.second[iflash];
408  if(flash->TotalPE()<fFlashPeThresh) continue;
409  double tflash = flash->Time()*1e3-fOpDelay;
410  if(abs(crt->ts1_ns-tflash)<fCoinWindow && abs(crt->ts1_ns-tflash)<flash_minTOF){
411  flash_minTOF= abs(crt->ts1_ns-tflash);
412  opflash_index = flash.key();
413  found_tof = true;
414  flash_tpc = flashList.first;
415  } // with in conincidence window and getting the largest flash
416  } // loop over flash list
417  } // loop over flash list map
418 
419  if(found_tof){
420  if(frm_trk){
421  tof_crt_hits[index].push_back(crt);
422  tof_op_flashes[index].push_back(opFlashLists[flash_tpc][opflash_index]);
423  tof_op_tpc[index].push_back(flash_tpc);
424  } // crt hit is coming from crt track
425 
426  else{
427  sbnd::ToF::ToF new_tof;
428  new_tof.tof = crt->ts1_ns - (opFlashLists[flash_tpc][opflash_index]->Time()*1e3-fOpDelay);
429  new_tof.frm_hit = true;
430  new_tof.crt_time = crt->ts1_ns;
431  new_tof.pmt_time = opFlashLists[flash_tpc][opflash_index]->Time()*1e3-fOpDelay;
432  new_tof.crt_tagger = crt->tagger;
433  new_tof.crt_hit_id = crt.key();
434  new_tof.flash_tpc_id = flash_tpc;
435  new_tof.pmt_flash_id = opflash_index;
436  ToF_vec->push_back(new_tof);
437  } // lonely crt hit
438  } // found a tof match
439 
440  } // use CFlash method
441 
442  //======================================================================================================================
443 
444  //=========================Calculation ToF values using Earliest hit of the Largest flash ==============================
445 
446  if(fLFlash_hit){
447  double peflash_max=0;
448  bool found_tof = false;
449  int opflash_index = -1;
450  int flash_tpc = -1;
451  int ophit_index = -1;
452 
453  std::map<int, art::Handle< std::vector<recob::OpFlash> > > flashHandles;
454  std::map<int,std::vector< art::Ptr<recob::OpFlash> >> opFlashLists;
455 
456  for(int i=0; i<2; i++) {
457  if( evt.getByLabel(fFlashLabels[i],flashHandles[i]) )
458  art::fill_ptr_vector(opFlashLists[i], flashHandles[i]);
459  }
460 
461  art::Handle< std::vector<recob::OpHit> > opHitListHandle;
462  std::vector< art::Ptr<recob::OpHit> > opHitList;
463  if( evt.getByLabel(fOpHitModuleLabel,opHitListHandle) )
464  art::fill_ptr_vector(opHitList, opHitListHandle);
465 
466  for(auto const& flashList : opFlashLists){
467  art::FindManyP<recob::OpHit> findManyOpHits(flashHandles[flashList.first], evt, fFlashLabels[flashList.first]);
468  for(size_t iflash=0; iflash<flashList.second.size(); iflash++){
469  auto const& flash = flashList.second[iflash];
470  if(flash->TotalPE()<fFlashPeThresh) continue;
471  double tflash = flash->AbsTime()*1e3-fOpDelay;
472  if(abs(crt->ts1_ns-tflash)<fCoinWindow && flash->TotalPE()>peflash_max){
473  peflash_max=flash->TotalPE();
474  opflash_index = flash.key();
475  found_tof = true;
476  flash_tpc = flashList.first;
477  vector<art::Ptr<recob::OpHit>> hits = findManyOpHits.at(flash.key());
478  double flashMinHitT = DBL_MAX;
479  for(auto const& hit : hits){
480  double tPmt = hit->PeakTime()*1e3-fOpDelay;
481  if(tPmt < flashMinHitT){
482  flashMinHitT = tPmt;
483  ophit_index = hit.key();
484  } // getting the earliest hit
485  } // loop over associated ophits of the flash
486  } // with in conincidence window and getting the largest flash
487  } // loop over flash list
488  } // loop over flash list map
489 
490  if(found_tof){
491  if(frm_trk){
492  tof_crt_hits[index].push_back(crt);
493  tof_op_flashes[index].push_back(opFlashLists[flash_tpc][opflash_index]);
494  tof_op_tpc[index].push_back(flash_tpc);
495  tof_op_hits[index].push_back(opHitList[ophit_index]);
496  } // crt hit is coming from crt track
497 
498  else{
499  sbnd::ToF::ToF new_tof;
500  new_tof.tof = crt->ts1_ns - (opHitList[ophit_index]->PeakTime()*1e3-fOpDelay);
501  new_tof.frm_hit = true;
502  new_tof.crt_time = crt->ts1_ns;
503  new_tof.pmt_time = opHitList[ophit_index]->PeakTime()*1e3-fOpDelay;
504  new_tof.crt_tagger = crt->tagger;
505  new_tof.crt_hit_id = crt.key();
506  new_tof.flash_tpc_id = flash_tpc;
507  new_tof.pmt_flash_id = opflash_index;
508  new_tof.pmt_hit_id = ophit_index;
509  ToF_vec->push_back(new_tof);
510  } // lonely crt hit
511  } // found a tof match
512 
513  } // use earliest optical hit of the largest flash
514 
515  //======================================================================================================================
516 
517  //=========================Calculation ToF values using Earliest hit of the Closest flash ==============================
518 
519  if(fCFlash_hit){
520  double flash_minTOF = DBL_MAX;
521  bool found_tof = false;
522  int opflash_index = -1;
523  int flash_tpc = -1;
524  int ophit_index = -1;
525 
526  std::map<int, art::Handle< std::vector<recob::OpFlash> > > flashHandles;
527  std::map<int,std::vector< art::Ptr<recob::OpFlash> >> opFlashLists;
528 
529  for(int i=0; i<2; i++) {
530  if( evt.getByLabel(fFlashLabels[i],flashHandles[i]) )
531  art::fill_ptr_vector(opFlashLists[i], flashHandles[i]);
532  }
533 
534  art::Handle< std::vector<recob::OpHit> > opHitListHandle;
535  std::vector< art::Ptr<recob::OpHit> > opHitList;
536  if( evt.getByLabel(fOpHitModuleLabel,opHitListHandle) )
537  art::fill_ptr_vector(opHitList, opHitListHandle);
538 
539  for(auto const& flashList : opFlashLists){
540  art::FindManyP<recob::OpHit> findManyOpHits(flashHandles[flashList.first], evt, fFlashLabels[flashList.first]);
541  for(size_t iflash=0; iflash<flashList.second.size(); iflash++){
542  auto const& flash = flashList.second[iflash];
543  if(flash->TotalPE()<fFlashPeThresh) continue;
544  double tflash = flash->Time()*1e3-fOpDelay;
545  if(abs(crt->ts1_ns-tflash)<fCoinWindow && abs(crt->ts1_ns-tflash)<flash_minTOF){
546  flash_minTOF= abs(crt->ts1_ns-tflash);
547  opflash_index = flash.key();
548  found_tof = true;
549  flash_tpc = flashList.first;
550  vector<art::Ptr<recob::OpHit>> hits = findManyOpHits.at(flash.key());
551  double flashMinHitT = DBL_MAX;
552  for(auto const& hit : hits){
553  double tPmt = hit->PeakTime()*1e3-fOpDelay;
554  if(tPmt < flashMinHitT){
555  flashMinHitT = tPmt;
556  ophit_index = hit.key();
557  } // getting the earliest hit
558  } // loop over associated ophits of the flash
559  } // with in conincidence window and getting the largest flash
560  } // loop over flash list
561  } // loop over flash list map
562 
563  if(found_tof){
564  if(frm_trk){
565  tof_crt_hits[index].push_back(crt);
566  tof_op_flashes[index].push_back(opFlashLists[flash_tpc][opflash_index]);
567  tof_op_tpc[index].push_back(flash_tpc);
568  tof_op_hits[index].push_back(opHitList[ophit_index]);
569  } // crt hit is coming from crt track
570 
571  else{
572  sbnd::ToF::ToF new_tof;
573  new_tof.tof = crt->ts1_ns - (opHitList[ophit_index]->PeakTime()*1e3-fOpDelay);
574  new_tof.frm_hit = true;
575  new_tof.crt_time = crt->ts1_ns;
576  new_tof.pmt_time = opHitList[ophit_index]->PeakTime()*1e3-fOpDelay;
577  new_tof.crt_tagger = crt->tagger;
578  new_tof.crt_hit_id = crt.key();
579  new_tof.flash_tpc_id = flash_tpc;
580  new_tof.pmt_flash_id = opflash_index;
581  new_tof.pmt_hit_id = ophit_index;
582  ToF_vec->push_back(new_tof);
583  } // lonely crt hit
584  } // found a tof match
585 
586  } // use earliest optical hit of the closest flash
587 
588  //======================================================================================================================
589 
590  } // loop over crt hit list
591 
592  //=================================== Calculation ToF values using Largest optical flash ===================================
593 
594  if(fLhit){
595  if(!tof_crt_hits.empty()){
596  for (auto& ele: tof_crt_hits){
597  double min_time = DBL_MAX;
598  int all_index = 0;
599  int min_index = 0;
600  for (auto const& hit: ele.second){
601  if(hit->ts1_ns < min_time){
602  min_time = hit->ts1_ns;
603  min_index = all_index;
604  }
605  all_index++;
606  }
607 
608  sbnd::ToF::ToF new_tof;
609  new_tof.tof = trackhits[ele.first].front()->ts1_ns - (tof_op_hits[ele.first][min_index]->PeakTime()*1e3-fOpDelay);
610  new_tof.crt_time = trackhits[ele.first].front()->ts1_ns;
611  new_tof.pmt_time = tof_op_hits[ele.first][min_index]->PeakTime()*1e3-fOpDelay;
612  new_tof.crt_tagger = trackhits[ele.first].front()->tagger;
613  new_tof.frm_trk = true;
614  new_tof.crt_trk_id = ele.first;
615  new_tof.pmt_hit_id = tof_op_hits[ele.first][min_index].key();
616  ToF_vec->push_back(new_tof);
617  }
618  }
619  }
620 
621  //=============================================================================================================================
622 
623  // ================================== Calculatin ToF values using Closest optical hit method ==================================
624 
625  if(fChit){
626  if(!tof_crt_hits.empty()){
627  for (auto& ele: tof_crt_hits){
628  double min_time = DBL_MAX;
629  int all_index = 0;
630  int min_index = 0;
631  for (auto const& hit: ele.second){
632  if(hit->ts1_ns < min_time){
633  min_time = hit->ts1_ns;
634  min_index = all_index;
635  }
636  all_index++;
637  }
638 
639  sbnd::ToF::ToF new_tof;
640  new_tof.tof = trackhits[ele.first].front()->ts1_ns - (tof_op_hits[ele.first][min_index]->PeakTime()*1e3-fOpDelay);
641  new_tof.crt_time = trackhits[ele.first].front()->ts1_ns;
642  new_tof.pmt_time = tof_op_hits[ele.first][min_index]->PeakTime()*1e3-fOpDelay;
643  new_tof.crt_tagger = trackhits[ele.first].front()->tagger;
644  new_tof.frm_trk = true;
645  new_tof.crt_trk_id = ele.first;
646  new_tof.pmt_hit_id = tof_op_hits[ele.first][min_index].key();
647  ToF_vec->push_back(new_tof);
648  }
649  }
650  }
651 
652  // ============================================================================================================================
653 
654  //=================================== Calculatin ToF values using Largest optical flash method ==================================
655 
656  if(fLFlash){
657  if(!tof_crt_hits.empty()){
658  for (auto& ele: tof_crt_hits){
659  double min_time = DBL_MAX;
660  int all_index = 0;
661  int min_index = 0;
662  for (auto const& hit: ele.second){
663  if(hit->ts1_ns < min_time){
664  min_time = hit->ts1_ns;
665  min_index = all_index;
666  }
667  all_index++;
668  }
669 
670  sbnd::ToF::ToF new_tof;
671  new_tof.tof = trackhits[ele.first].front()->ts1_ns - (tof_op_flashes[ele.first][min_index]->Time()*1e3-fOpDelay);
672  new_tof.crt_time = trackhits[ele.first].front()->ts1_ns;
673  new_tof.pmt_time = tof_op_flashes[ele.first][min_index]->Time()*1e3-fOpDelay;
674  new_tof.crt_tagger = trackhits[ele.first].front()->tagger;
675  new_tof.frm_trk = true;
676  new_tof.crt_trk_id = ele.first;
677  new_tof.flash_tpc_id = tof_op_tpc[ele.first][min_index];
678  new_tof.pmt_flash_id = tof_op_flashes[ele.first][min_index].key();
679  ToF_vec->push_back(new_tof);
680  }
681  }
682  }
683 
684  //=============================================================================================================================
685 
686  //=================================== Calculatin ToF values using Closest optical flash method ================================
687 
688  if(fCFlash){
689  if(!tof_crt_hits.empty()){
690  for (auto& ele: tof_crt_hits){
691  double min_time = DBL_MAX;
692  int all_index = 0;
693  int min_index = 0;
694  for (auto const& hit: ele.second){
695  if(hit->ts1_ns < min_time){
696  min_time = hit->ts1_ns;
697  min_index = all_index;
698  }
699  all_index++;
700  }
701 
702  sbnd::ToF::ToF new_tof;
703  new_tof.tof = trackhits[ele.first].front()->ts1_ns - (tof_op_flashes[ele.first][min_index]->Time()*1e3-fOpDelay);
704  new_tof.crt_time = trackhits[ele.first].front()->ts1_ns;
705  new_tof.pmt_time = tof_op_flashes[ele.first][min_index]->Time()*1e3-fOpDelay;
706  new_tof.crt_tagger = trackhits[ele.first].front()->tagger;
707  new_tof.frm_trk = true;
708  new_tof.crt_trk_id = ele.first;
709  new_tof.flash_tpc_id = tof_op_tpc[ele.first][min_index];
710  new_tof.pmt_flash_id = tof_op_flashes[ele.first][min_index].key();
711  ToF_vec->push_back(new_tof);
712  }
713  }
714  }
715 
716  //=============================================================================================================================
717 
718  //=========================Calculation ToF values using Earliest hit of the Largest flash =====================================
719 
720  if(fLFlash_hit){
721  if(!tof_crt_hits.empty()){
722  for (auto& ele: tof_crt_hits){
723  double min_time = DBL_MAX;
724  int all_index = 0;
725  int min_index = 0;
726  for (auto const& hit: ele.second){
727  if(hit->ts1_ns < min_time){
728  min_time = hit->ts1_ns;
729  min_index = all_index;
730  }
731  all_index++;
732  }
733 
734  sbnd::ToF::ToF new_tof;
735  new_tof.tof = trackhits[ele.first].front()->ts1_ns - (tof_op_hits[ele.first][min_index]->PeakTime()*1e3-fOpDelay);
736  new_tof.crt_time = trackhits[ele.first].front()->ts1_ns;
737  new_tof.pmt_time = tof_op_hits[ele.first][min_index]->PeakTime()*1e3-fOpDelay;
738  new_tof.crt_tagger = trackhits[ele.first].front()->tagger;
739  new_tof.frm_trk = true;
740  new_tof.crt_trk_id = ele.first;
741  new_tof.flash_tpc_id = tof_op_tpc[ele.first][min_index];
742  new_tof.pmt_flash_id = tof_op_flashes[ele.first][min_index].key();
743  new_tof.pmt_hit_id = tof_op_hits[ele.first][min_index].key();
744  ToF_vec->push_back(new_tof);
745  }
746  }
747  }
748 
749  //============================================================================================================================
750 
751  //=========================Calculation ToF values using Earliest hit of the Closest flash ====================================
752 
753  if(fCFlash_hit){
754  if(!tof_crt_hits.empty()){
755  for (auto& ele: tof_crt_hits){
756  double min_time = DBL_MAX;
757  int all_index = 0;
758  int min_index = 0;
759  for (auto const& hit: ele.second){
760  if(hit->ts1_ns < min_time){
761  min_time = hit->ts1_ns;
762  min_index = all_index;
763  }
764  all_index++;
765  }
766 
767  sbnd::ToF::ToF new_tof;
768  new_tof.tof = trackhits[ele.first].front()->ts1_ns - (tof_op_hits[ele.first][min_index]->PeakTime()*1e3-fOpDelay);
769  new_tof.crt_time = trackhits[ele.first].front()->ts1_ns;
770  new_tof.pmt_time = tof_op_hits[ele.first][min_index]->PeakTime()*1e3-fOpDelay;
771  new_tof.crt_tagger = trackhits[ele.first].front()->tagger;
772  new_tof.frm_trk = true;
773  new_tof.crt_trk_id = ele.first;
774  new_tof.flash_tpc_id = tof_op_tpc[ele.first][min_index];
775  new_tof.pmt_flash_id = tof_op_flashes[ele.first][min_index].key();
776  new_tof.pmt_hit_id = tof_op_hits[ele.first][min_index].key();
777  ToF_vec->push_back(new_tof);
778  }
779  }
780  }
781  //============================================================================================================================
782 
783  evt.put(std::move(ToF_vec));
784 }
785 
786 bool ToFProducer::HitCompare(const art::Ptr<CRTHit>& hit1, const art::Ptr<CRTHit>& hit2) {
787 
788  if(hit1->ts1_ns != hit2->ts1_ns) return false;
789  if(hit1->plane != hit2->plane) return false;
790  if(hit1->x_pos != hit2->x_pos) return false;
791  if(hit1->y_pos != hit2->y_pos) return false;
792  if(hit1->z_pos != hit2->z_pos) return false;
793  if(hit1->x_err != hit2->x_err) return false;
794  if(hit1->y_err != hit2->y_err) return false;
795  if(hit1->z_err != hit2->z_err) return false;
796  if(hit1->tagger != hit2->tagger) return false;
797  return true;
798 }
799 
800 DEFINE_ART_MODULE(ToFProducer)
801 }
sbnd::ToF::ToF::crt_hit_id
int crt_hit_id
Definition: ToF.hh:18
sbnd::ToFProducer::HitCompare
bool HitCompare(const art::Ptr< CRTHit > &h1, const art::Ptr< CRTHit > &h2)
Definition: ToFProducer_module.cc:786
sbnd::ToFProducer::fGeometryService
geo::GeometryCore const * fGeometryService
Definition: ToFProducer_module.cc:153
pmtana::std
double std(const std::vector< short > &wf, const double ped_mean, size_t start, size_t nsample)
Definition: UtilFunc.cxx:42
sbnd::ToF::ToF::crt_trk_id
int crt_trk_id
Definition: ToF.hh:19
PmaTrack3D.h
Implementation of the Projection Matching Algorithm.
sbnd::ToFProducer::fOpFlashModuleLabel1
art::InputTag fOpFlashModuleLabel1
Definition: ToFProducer_module.cc:133
sbnd::ToFProducer::fOpHitModuleLabel
art::InputTag fOpHitModuleLabel
Definition: ToFProducer_module.cc:131
sbnd::ToF::ToF::crt_tagger
std::string crt_tagger
Definition: ToF.hh:17
Hit.h
Declaration of signal hit object.
sbnd::ToF::ToF::crt_time
float crt_time
Definition: ToF.hh:15
sbnd::ToF::ToF::tof
float tof
Definition: ToF.hh:12
RawDigit.h
Definition of basic raw digits.
sbnd::ToFProducer::beginJob
void beginJob() override
Definition: ToFProducer_module.cc:183
hit
process_name hit
Definition: cheaterreco.fcl:51
a
process_name gaushit a
Definition: decode_signalprocess_icarus.fcl:37
sbnd::ToF::ToF::frm_trk
bool frm_trk
Definition: ToF.hh:13
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
sbnd::ToFProducer::fCrtHitModuleLabel
art::InputTag fCrtHitModuleLabel
Definition: ToFProducer_module.cc:134
AuxDetSimChannel.h
object containing MC truth information necessary for making RawDigits and doing back tracking ...
sbnd::ToFProducer::produce
void produce(art::Event &evt) override
Definition: ToFProducer_module.cc:188
sbnd::ToF::ToF::flash_tpc_id
int flash_tpc_id
Definition: ToF.hh:22
sbnd::ToFProducer::fFlashLabels
map< int, art::InputTag > fFlashLabels
Definition: ToFProducer_module.cc:152
sbnd::ToFProducer::fOpFlashModuleLabel0
art::InputTag fOpFlashModuleLabel0
Definition: ToFProducer_module.cc:132
geo::GeometryCore
Description of geometry of one entire detector.
Definition: GeometryCore.h:1476
Cluster.h
Declaration of cluster object.
geo_types.h
Definition of data types for geometry description.
Track.h
Provides recob::Track data product.
sbnd::ToFProducer::bt
CRTBackTracker * bt
Definition: ToFProducer_module.cc:151
sbnd::ToF::ToF::pmt_hit_id
int pmt_hit_id
Definition: ToF.hh:20
SimChannel.h
object containing MC truth information necessary for making RawDigits and doing back tracking ...
sbnd::ToFProducer::fCrtTrackModuleLabel
art::InputTag fCrtTrackModuleLabel
Definition: ToFProducer_module.cc:135
MCShower.h
Class def header for MCShower data container.
sbnd
stream1 can override from command line with o or output services user sbnd
Definition: prodgenie_sbnd_proj.fcl:101
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
sbnd::ToF::ToF::frm_hit
bool frm_hit
Definition: ToF.hh:14
sbnd::ToF::ToF::pmt_flash_id
int pmt_flash_id
Definition: ToF.hh:21
sbnd::ToF::ToF::pmt_time
float pmt_time
Definition: ToF.hh:16
crt
process_name crt
Definition: run_genericCRT.fcl:82
LArPandoraHelper.h
helper function for LArPandoraInterface producer module
Geometry.h
art framework interface to geometry description
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
Generated by   doxygen 1.8.5