PDSPmatch_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: PDSPmatch
3 // Plugin Type: analyzer (art v3_00_00)
4 // File: PDSPmatch_module.cc
5 //
6 // Generated at Thu Feb 14 21:47:27 2019 by Bryan Joseph Ramson using cetskelgen
7 // from cetlib version v3_04_00.
8 ////////////////////////////////////////////////////////////////////////
9 
10 #include <numeric>
11 #include <sstream>
12 #include <string.h>
13 #include <bitset>
14 #include <vector>
15 
17 
25 #include "fhiclcpp/ParameterSet.h"
27 #include "art_root_io/TFileService.h"
28 //#include "art/Framework/Services/Optional/TFileService.h"
29 
30 
31 #include "dune/Protodune/singlephase/CTB/data/pdspctb.h"
32 #include "dune/Protodune/singlephase/CRT/data/CRTTrigger.h"
33 //#include "dune/Geometry/ProtoDUNESPCRTSorter.h"
34 
36 #include "artdaq-core/Data/ContainerFragment.hh"
37 
39 #include "TTree.h"
40 #include "TH2.h"
41 #include "TLorentzVector.h"
42 #include "TVector3.h"
43 
56 #include "canvas/Persistency/Common/FindManyP.h"
57 
58 namespace pdsp {
59  class PDSPmatch;
60 }
61 
62 
63 class pdsp::PDSPmatch : public art::EDAnalyzer {
64 public:
65  explicit PDSPmatch(fhicl::ParameterSet const& p);
66  // The compiler-generated destructor is fine for non-base
67  // classes without bare pointers or other resource use.
68 
69  // Plugins should not be copied or assigned.
70  PDSPmatch(PDSPmatch const&) = delete;
71  PDSPmatch(PDSPmatch&&) = delete;
72  PDSPmatch& operator=(PDSPmatch const&) = delete;
73  PDSPmatch& operator=(PDSPmatch&&) = delete;
74 
75  // Required functions.
76  void analyze(art::Event const& e) override;
77 
78  // Selected optional functions.
79  void beginJob() override;
80 
81 private:
82  const art::InputTag fTimeLabel; //Label of timing products
83  const art::InputTag fCTBLabel; //Label of pdsctbdata products
84  const art::InputTag fCRTLabel; //Label of crt products
85  const art::InputTag fOpHitLabel; //Label of ophits after rawdecoding
86  const art::InputTag fPandoLabel; //Label for Pandora Tracks
87  const art::InputTag fPFParListLabel; // Label for PF Particle containers
88 
89  const int64_t fCRTCTBOffset; //CRT offset to global trigger in 20mus ticks
90  const uint64_t fCRTWindow; //coincidence window for CRT triggers in 20mus ticks
91  const bool fMC; //MC flag
92 
93  TTree *fTree;
94  int run;
95  int subrun;
96  int event;
97  int64_t timing_time;
98  int64_t fCTB_time;
99  uint32_t fCTBChan;
100 
101  // TH1F *CTBCRTMap[32];
102  //TH1F *CTBPairPE[16];
103  //TH1F *CosTrigger[16];
104  //TH1F *CRTPairPE[16];
105  //TH2F *PandaTrackProUp[16];
106  //TH2F *PandaTrackProDown[16];
107 
108  //TH1F *PMTrack[16];
109  //TH1F *CRTPairPE[16];
110  //TH1F *trig
111 
112  std::vector<uint32_t> fCRTChan;
113  std::vector<int64_t> fPDS_time, fPando_time, fCRT_time;
114 
115  std::vector<double> fOpChan, fPE;//, ft0, fx0, fy0, fz0;
116 
119  std::vector<double> fTrigCos_Pando;
120  std::vector<uint32_t> fTrkProUS_Pando, fTrkProDS_Pando;
121  std::vector<double> fTrkUSPixelProx_Pando, fTrkUSPixelProy_Pando;
122  std::vector<double> fTrkDSPixelProx_Pando, fTrkDSPixelProy_Pando;
123 
124  std::vector<TVector3> PandoTracks, PMTracks;
125  // Declare member data here.
126 
127  TVector3 CTBPixelCtr[32];
128 };
129 
130 
132  :
133  EDAnalyzer(p),
134  fTimeLabel(p.get<art::InputTag>("TimingLabel")),
135  fCTBLabel(p.get<art::InputTag>("CTBLabel")),
136  fCRTLabel(p.get<art::InputTag>("CRTLabel")),
137  fOpHitLabel(p.get<art::InputTag>("OpHitLabel")),
138  fPandoLabel(p.get<art::InputTag>("PandoraLabel")),
139  fPFParListLabel(p.get<art::InputTag>("PFParListLabel")),
140  fCRTCTBOffset(p.get<int64_t>("CRTCTBOffset")),
141  fCRTWindow(p.get<uint64_t>("CRTWindow")),
142  fMC(p.get<bool>("MC"))
143  //: EDAnalyzer{p} // ,
144  //More initializers here.
145 {
146  // Call appropriate consumes<>() for any products to be retrieved by this module.
147 
148  consumes<std::vector<raw::RDTimeStamp>>(fTimeLabel);
149  consumes<std::vector<raw::ctb::pdspctb>>(fCTBLabel);
150  consumes<std::vector<CRT::Trigger>>(fCRTLabel);
151  consumes<std::vector<recob::OpHit>>(fOpHitLabel);
152  consumes<std::vector<recob::Track>>(fPandoLabel);
153  consumes<std::vector<recob::PFParticle>>(fPFParListLabel);
154 }
155 
158  CTBPixelCtr[0].SetXYZ(-111.494,580.249,-972.204);
159  CTBPixelCtr[1].SetXYZ(58.5057,580.249,-972.204);
160  CTBPixelCtr[2].SetXYZ(307.506,580.249,-267.204);
161  CTBPixelCtr[3].SetXYZ(477.506,580.249,-267.204);
162  CTBPixelCtr[4].SetXYZ(477.506,410.249,-267.204);
163  CTBPixelCtr[5].SetXYZ(477.506,235.249,-267.204);
164  CTBPixelCtr[6].SetXYZ(477.506,65.2494,-267.204);
165  CTBPixelCtr[7].SetXYZ(307.506,65.2494,-267.204);
166  CTBPixelCtr[8].SetXYZ(58.5057,65.2494,-972.204);
167  CTBPixelCtr[9].SetXYZ(-111.494,65.2494,-972.204);
168  CTBPixelCtr[10].SetXYZ(-111.494,235.249,-972.204);
169  CTBPixelCtr[11].SetXYZ(-111.494,410.249,-972.204);
170  CTBPixelCtr[12].SetXYZ(58.5057,410.249,-972.204);
171  CTBPixelCtr[13].SetXYZ(307.506,410.249,-267.204);
172  CTBPixelCtr[14].SetXYZ(307.506,235.249,-267.204);
173  CTBPixelCtr[15].SetXYZ(58.5057,235.249,-972.204);
174  CTBPixelCtr[16].SetXYZ(293.506,580.249,1085.55);
175  CTBPixelCtr[17].SetXYZ(123.506,580.249,1085.55);
176  CTBPixelCtr[18].SetXYZ(-51.4993,580.249,1085.55);
177  CTBPixelCtr[19].SetXYZ(-221.494,580.249,1085.55);
178  CTBPixelCtr[20].SetXYZ(-221.494,410.249,1085.55);
179  CTBPixelCtr[21].SetXYZ(-221.494,235.249,1085.55);
180  CTBPixelCtr[22].SetXYZ(-221.494,65.2494,1085.55);
181  CTBPixelCtr[23].SetXYZ(-51.4943,65.2494,1085.55);
182  CTBPixelCtr[24].SetXYZ(123.506,65.2494,1085.55);
183  CTBPixelCtr[25].SetXYZ(293.506,65.2494,1085.55);
184  CTBPixelCtr[26].SetXYZ(293.506,235.249,1085.55);
185  CTBPixelCtr[27].SetXYZ(293.506,410.249,1085.55);
186  CTBPixelCtr[28].SetXYZ(123.506,410.249,1085.55);
187  CTBPixelCtr[29].SetXYZ(-51.4993,410.249,1085.55);
188  CTBPixelCtr[30].SetXYZ(-51.4943,235.249,1085.55);
189  CTBPixelCtr[31].SetXYZ(123.506,235.249,1085.55);
190 
191  for(int h=0;h<32;h++){
192  if(h==0 || h==1 || h==11 || h==12 || h ==10 || h==15 || h==9 || h==8) CTBPixelCtr[h].SetY(CTBPixelCtr[h].Y()-46.0); //US Jura/Beam Left
193  if(h==2 || h==3 || h==13 || h==4 || h ==14 || h==5 || h==7 || h==6) CTBPixelCtr[h].SetY(CTBPixelCtr[h].Y()-47.5); //US Saleve/Beam Right
194  if(h==16 || h==17 || h==27 || h==28 || h ==26 || h==31 || h==25 || h==24) CTBPixelCtr[h].SetY(CTBPixelCtr[h].Y()-151.0); //DS Saleve/Beam Right
195  if(h==18 || h==19 || h==29 || h==20 || h ==21 || h==30 || h==22 || h==23) CTBPixelCtr[h].SetY(CTBPixelCtr[h].Y()-150.0); //DS Jura/Beam Left
196  }
198 
199  fTree = tFileService->make<TTree>("ProtoDUNE_Evt_Match","ProtoDUNE Matched Event Tree");
200 
201  fTree->Branch("run",&run,"run/I");
202  fTree->Branch("subrun",&subrun,"subrun/I");
203  fTree->Branch("event",&event,"event/I");
204  fTree->Branch("timing_time",&timing_time,"timing_time/L");
205 
206  fTree->Branch("fCRT_time",&fCRT_time);
207  fTree->Branch("fCRTChan", &fCRTChan);
208  fTree->Branch("fCTB_time",&fCTB_time,"fCTB_time/L");
209  fTree->Branch("fCTBChan", &fCTBChan,"fCTBChan/i");
210 
211  fTree->Branch("fPDS_time",&fPDS_time);
212  fTree->Branch("fOpChan", &fOpChan);
213  fTree->Branch("fPE",&fPE);
214 
215  fTree->Branch("fPando_time",&fPando_time);
216  fTree->Branch("fTrkStartx_Pando",&fTrkStartx_Pando);
217  fTree->Branch("fTrkStarty_Pando",&fTrkStarty_Pando);
218  fTree->Branch("fTrkStartz_Pando",&fTrkStartz_Pando);
219  fTree->Branch("fTrkEndx_Pando",&fTrkEndx_Pando);
220  fTree->Branch("fTrkEndy_Pando",&fTrkEndy_Pando);
221  fTree->Branch("fTrkEndz_Pando",&fTrkEndz_Pando);
222 
223 
224 }
225 
227  run = e.run();
228  subrun = e.subRun();
229  event = e.id().event();
230 
231  fCRT_time.clear();
232  fCRTChan.clear();
233 
234  fPDS_time.clear();
235  fOpChan.clear();
236  fPE.clear();
237 
238  fPando_time.clear();
239  fTrkStartx_Pando.clear();
240  fTrkStarty_Pando.clear();
241  fTrkStartz_Pando.clear();
242  fTrkEndx_Pando.clear();
243  fTrkEndy_Pando.clear();
244  fTrkEndz_Pando.clear();
245  fTrkTheta_Pando.clear();
246  fTrkPhi_Pando.clear();
247 
248  int trigger = -1;
249  const auto timeStamps = e.getValidHandle<std::vector<raw::RDTimeStamp>>(fTimeLabel);
250  if(timeStamps.isValid() && timeStamps->size() == 1){
251  const raw::RDTimeStamp& timeStamp = timeStamps->at(0);
252  trigger = timeStamp.GetFlags();
253  timing_time = timeStamp.GetTimeStamp();
254  }
255  else {
256  mf::LogWarning("Empty Event TimeStamp") << "Invalid Event TimeStamp. Skipping. \n";
257  return;
258  }
259 
260  if (trigger !=13) return;
261 
262  const auto crtHandle = e.getValidHandle<std::vector<CRT::Trigger>>(fCRTLabel);
263  if(crtHandle->empty()){
264  mf::LogWarning("Empty CRT Fragment") << "Empty CRT fragments for this event. Skipping. \n";
265  return;
266  }
267 
268  const auto OpHitHandle = e.getValidHandle<std::vector<recob::OpHit>>(fOpHitLabel);
269  if(OpHitHandle->empty()){
270  mf::LogWarning("Empty OpHit Object") << "Empty OpHit Object. Error in retrieval. Skipping. \n";
271  return;
272  }
273 
274  const auto ctbHandle = e.getValidHandle<std::vector<raw::ctb::pdspctb>>(fCTBLabel);
275  if(ctbHandle->empty()){
276  mf::LogWarning("Empty CTB Fragment") << "Empty CTB fragment for this event. Skipping. \n";
277  return;
278  }
279 
280  if(ctbHandle->size() > 1){
281  mf::LogWarning("Multiple CTB Triggers") << "Found " << ctbHandle->size() << " CTB data products. Skipping. \n";
282  return;
283  }
284 
285  const auto& ctbStatus = ctbHandle->front();
286  const auto statuses = ctbStatus.GetChStatusAfterHLTs();
287 
288  bool pat = false;
289  for(const auto& status: statuses){
290  pat = false;
291 
292  fCRT_time.clear();
293  fCRTChan.clear();
294 
295  fPDS_time.clear();
296  fOpChan.clear();
297  fPE.clear();
298 
299  fPando_time.clear();
300  fTrkStartx_Pando.clear();
301  fTrkStarty_Pando.clear();
302  fTrkStartz_Pando.clear();
303  fTrkEndx_Pando.clear();
304  fTrkEndy_Pando.clear();
305  fTrkEndz_Pando.clear();
306 
307  if(status.crt != 0) {
308  fCTB_time = status.timestamp;
309  fCTBChan = status.crt;
310  if(!crtHandle->empty()){
311  std::vector<CRT::Trigger> inWindow(crtHandle->size());
312  const int64_t ctbetime = status.timestamp;
313  const auto newEnd = std::copy_if(crtHandle->begin(), crtHandle->end(), inWindow.begin(),[this,ctbetime](const auto& trigger){
314  return (uint64_t)(labs(ctbetime - ((int64_t)(trigger.Timestamp()) - fCRTCTBOffset))) < fCRTWindow; });
315  inWindow.resize(std::distance(inWindow.begin(), newEnd));
316  for(size_t k=0;k<inWindow.size();++k){
317  fCRT_time.push_back(inWindow[k].Timestamp());
318  fCRTChan.push_back(inWindow[k].Channel());
319  }
320  }
321  for(const auto& OpHit: *OpHitHandle){
322  if(OpHit.PE() < 10.0 || OpHit.PE() > 1000.00) continue;
323  fPDS_time.push_back((OpHit.PeakTime()/3));
324  fOpChan.push_back(OpHit.OpChannel());
325  fPE.push_back(OpHit.PE());
326  }
327 
328  std::vector<art::Ptr<recob::Track>> PandoTrk;
329  auto PandoTrkHandle = e.getHandle<std::vector<recob::Track>>(fPandoLabel);
330  if (PandoTrkHandle) art::fill_ptr_vector(PandoTrk,PandoTrkHandle);
331  else {
332  mf::LogWarning("Empty PandoTrk Fragment") << "Empty PandoTrk Vector for this event. Skipping. \n";
333  return;
334  }
335 
336  auto PFParListHandle = e.getHandle<std::vector<recob::PFParticle>>(fPFParListLabel);
337  if(!PFParListHandle){;
338  mf::LogWarning("Empty PFParticle Vector") << "Empty PFParticle Vector for this event. Skipping. \n";
339  return;
340  }
341  art::FindManyP<recob::PFParticle> PFPar(PandoTrkHandle,e,fPandoLabel);
342  art::FindManyP<anab::T0> PFT0(PFParListHandle,e,fPFParListLabel);
343 
344  for(size_t p = 0;p<PandoTrk.size();++p){
345  auto & Trk = PandoTrk[p];
346  if(!((Trk->Vertex().Z() < 40) || (Trk->End().Z() < 40))) continue;
347  if(!((Trk->Vertex().Z() > 660) || (Trk->End().Z() > 660))) continue;
348  fTrkStartx_Pando.push_back(Trk->Vertex().X());
349  fTrkStarty_Pando.push_back(Trk->Vertex().Y());
350  fTrkStartz_Pando.push_back(Trk->Vertex().Z());
351  fTrkEndx_Pando.push_back(Trk->End().X());
352  fTrkEndy_Pando.push_back(Trk->End().Y());
353  fTrkEndz_Pando.push_back(Trk->End().Z());
354  double t0temp = 0;
355  auto &PFPS = PFPar.at(Trk.key());
356  if(!PFPS.empty()){
357  auto &T0S = PFT0.at(PFPS[0].key());
358  if(!T0S.empty()){
359  t0temp = T0S[0]->Time();
360  }
361  }
362  fPando_time.push_back(t0temp);
363  }
364  if(PandoTrk.size() > 0) pat = true;
365 
366  if(pat) fTree->Fill();
367  }
368  }
369 }
PDSPmatch(fhicl::ParameterSet const &p)
std::vector< double > fPE
const art::InputTag fCTBLabel
const art::InputTag fOpHitLabel
uint16_t GetFlags() const
Definition: RDTimeStamp.h:46
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
std::vector< double > fTrkStartx_Pando
std::vector< double > fTrkDSPixelProx_Pando
const art::InputTag fCRTLabel
std::vector< int64_t > fCRT_time
TVector3 CTBPixelCtr[32]
std::vector< TVector3 > PandoTracks
const uint64_t fCRTWindow
std::vector< double > fTrkEndx_Pando
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
const art::InputTag fPandoLabel
const art::InputTag fPFParListLabel
std::vector< double > fTrkEndz_Pando
art framework interface to geometry description
void analyze(art::Event const &e) override
std::vector< int64_t > fPDS_time
std::vector< double > fTrkDSPosy_Pando
const double e
void beginJob() override
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
std::vector< double > fTrkDSPosx_Pando
std::vector< uint32_t > fCRTChan
std::vector< double > fTrkDSPixelProy_Pando
def key(type, name=None)
Definition: graph.py:13
const int64_t fCRTCTBOffset
std::vector< double > fTrkPhi_Pando
std::vector< uint32_t > fTrkProDS_Pando
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:441
p
Definition: test.py:223
std::vector< TVector3 > PMTracks
std::vector< double > fTrkUSPos1y_Pando
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:78
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
RunNumber_t run() const
Definition: DataViewImpl.cc:71
std::vector< double > fTrkUSPos2y_Pando
std::vector< double > fTrkStartz_Pando
std::vector< int64_t > fPando_time
ULong64_t GetTimeStamp() const
Definition: RDTimeStamp.h:42
Encapsulate the geometry of an optical detector.
std::vector< double > fOpChan
std::vector< double > fTrkUSPixelProx_Pando
const art::InputTag fTimeLabel
std::vector< uint32_t > fTrkProUS_Pando
std::vector< double > fTrkStarty_Pando
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Provides recob::Track data product.
EventNumber_t event() const
Definition: EventID.h:116
Access the description of detector geometry.
PDSPmatch & operator=(PDSPmatch const &)=delete
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
std::vector< double > fTrkTheta_Pando
std::vector< double > fTrkUSPixelProy_Pando
int bool
Definition: qglobal.h:345
std::vector< double > fTrkUSPos1x_Pando
EventID id() const
Definition: Event.cc:34
std::vector< double > fTrkUSPos2x_Pando
Event finding and building.
pure virtual base interface for detector clocks
std::vector< double > fTrkEndy_Pando
std::vector< double > fTrigCos_Pando