hitrmsrawdigits_module.cc
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////////
2 /// Class: hitrmsrawdigits ///
3 /// File: hitrmsrawdigits_module.cc ///
4 /// Description: Module for calculations with hit information based on ///
5 /// RawDigits and/or from after deconvolution. Information is ///
6 /// taken from MC samples or data files and put into TTrees ///
7 /// to be analyzed. ///
8 /// Contact Person: ehinkle@uchicago.edu ///
9 /// Written by: Ajib Paudel ///
10 /// Modified by: Elise Hinkle ///
11 /// Last Date Modified: August 20, 2021 ///
12 //////////////////////////////////////////////////////////////////////////////
13 
18 #include "art_root_io/TFileService.h"
19 #include "art_root_io/TFileDirectory.h"
21 #include "canvas/Persistency/Common/FindManyP.h"
24 #include "fhiclcpp/ParameterSet.h"
26 #include "dune/Protodune/singlephase/DataUtils/ProtoDUNEDataUtils.h"
28 
38 #include "lardataobj/RawData/raw.h"
62 
63 #include "TFile.h"
64 #include "TTree.h"
65 #include "TDirectory.h"
66 #include "TH1.h"
67 #include "TH2.h"
68 #include "TF1.h"
69 #include "TProfile.h"
70 #include "TROOT.h"
71 #include "TStyle.h"
72 #include "TMath.h"
73 #include "TGraphErrors.h"
74 #include "TMinuit.h"
75 #include "TString.h"
76 #include "TTimeStamp.h"
77 #include "TVectorD.h"
78 #include "TCanvas.h"
79 #include "TFrame.h"
80 #include "TLine.h"
81 #include "TAxis.h"
82 #include "TTimeStamp.h"
83 
84 //Root and C++ include
85 #include "TVector3.h"
86 #include "TLorentzVector.h"
87 #include <vector>
88 #include <fstream>
89 #include "TPaveStats.h"
90 #include <iostream>
91 #include <string>
92 #include "math.h"
93 #include "stdio.h"
94 #include <iterator>
95 
96 //using namespace std;
97 
98 
99 namespace protoana{
101  public:
102  explicit hitrmsrawdigits(fhicl::ParameterSet const& pset);
103  virtual ~hitrmsrawdigits();
104 
105  void beginJob();
106  void endJob();
107  void beginRun(const art::Run& run);
108  void analyze(const art::Event& evt);
109  void reset();
110 
111  private:
112  unsigned int fWaveformSize; // Full waveform size
114  TTree* fEventTree;
116 
117  //These are the tree variables I will be using
118  Int_t run;
119  Int_t subrun;
120  Int_t event;
121  Double_t evttime;
123  std::vector<float> trackthetaxz;
124  std::vector<float> trackthetayz;
125  std::vector<float> trkstartx;
126  std::vector<float> trkstarty;
127  std::vector<float> trkstartz;
128  std::vector< std::vector<float> > trkstartcosxyz;
129  std::vector< std::vector<float> > trkendcosxyz;
130  std::vector<float> trkendx;
131  std::vector<float> trkendy;
132  std::vector<float> trkendz;
133 
134  std::vector<float> trkstartx_crt2;
135  std::vector<float> trkendx_crt2;
136  std::vector<float> crtreco_x0;
137  std::vector<float> crtreco_x1;
138  std::vector<float> crtreco_y0;
139  std::vector<float> crtreco_y1;
140  std::vector<float> crtreco_z0;
141  std::vector<float> crtreco_z1;
142 
143 
144 
145  std::vector<float> trklen;
146  std::vector<int> TrkID;
147  std::vector<float> xprojectedlen;
148  std::vector<double> T0_values;
149  // std::vector<double> t0crt1;
150  std::vector<double> t0crt2;
151  std::vector<double> crt2tickoffset;
152  std::vector<int> tot_trks;
153  std::vector< std::vector<float> > hit_peakT0;
154  std::vector< std::vector<int> > hit_tpc0;
155  std::vector< std::vector<int> > hit_wire0;
156  std::vector< std::vector<float> > hit_rms0;
157  std::vector< std::vector<float> > hit_deltaT;
158  std::vector< std::vector<float> > trkhitx0;
159  std::vector< std::vector<float> > trkhity0;
160  std::vector< std::vector<float> > trkhitz0;
161  std::vector< std::vector<float> > trkdq_amp0;
162  std::vector< std::vector<float> >trkdq_int0;
163  std::vector< std::vector<float> > hit_peakT1;
164  std::vector< std::vector<int> > hit_tpc1;
165  std::vector< std::vector<int> > hit_wire1;
166  std::vector< std::vector<float> > hit_rms1;
167  std::vector< std::vector<float> > trkhitx1;
168  std::vector< std::vector<float> > trkhity1;
169  std::vector< std::vector<float> > trkhitz1;
170  std::vector< std::vector<float> > trkdq_amp1;
171  std::vector< std::vector<float> > trkdq_int1;
172  std::vector< std::vector<float> > hit_peakT2;
173  std::vector< std::vector<int> > hit_tpc2;
174  std::vector< std::vector<int> > hit_wire2;
175  std::vector< std::vector<float> > hit_rms2;
176  std::vector< std::vector<float> > hit_rmsraw2;
177  std::vector< std::vector<float> > hit_peakTraw2;
178  float hit_signal[10][5000][60];//[TrackIndex, hitindex,peakTindex]signal for each hit peaktime-30 to peaktime+30 ticks
179  //float hit_signal[10][5000][50];//[TrackIndex, hitindex,peakTindex]signal for each hit peaktime-20 to peaktime+30 ticks
180  std::vector< std::vector<float> > hit_rms_true;
181  std::vector< std::vector<float> > trkhitx2;
182  std::vector< std::vector<float> > trkhity2;
183  std::vector< std::vector<float> > trkhitz2;
184  std::vector< std::vector<float> > trkhitz_wire2;
185  std::vector< std::vector<float> > trkdq_amp2;
186  std::vector< std::vector<float> > trkdq_int2;
187  std::vector< std::vector<float> > multiplicity2;
188  std::vector< std::vector<float> > goodnessoffit2;
196  };
197 
198  //========================================================================
200  EDAnalyzer(pset),
201  fWaveformSize (pset.get<unsigned int>("WaveformSize", 6000)),
202  fDataUtils (pset.get<fhicl::ParameterSet>("DataUtils")),
203  fHitsModuleLabel (pset.get< std::string >("HitsModuleLabel","") ),
204  fTrackModuleLabel (pset.get< std::string >("TrackModuleLabel","") ),
205  fCalorimetryModuleLabel (pset.get< std::string >("CalorimetryModuleLabel","") ),
206  fSaveTrackInfo (pset.get< bool>("SaveTrackInfo",false)),
207  fSaveCaloInfo (pset.get< bool>("SaveCaloInfo",false)),
208  fRawProducerLabel (pset.get<art::InputTag>("RawProducerLabel","")),
209  fWireProducerLabel (pset.get<art::InputTag>("WireProducerLabel",""))
210 
211  {
212  if (fSaveTrackInfo == false) fSaveCaloInfo = false;
213  }
214 
215  //========================================================================
217  }
218  //========================================================================
219 
220  //========================================================================
222  std::cout<<"job begin..."<<std::endl;
224  fEventTree = tfs->make<TTree>("Event", "Event Tree from Reco");
225  fEventTree->Branch("event", &event,"event/I");
226  fEventTree->Branch("evttime",&evttime,"evttime/D");
227  fEventTree->Branch("run", &run,"run/I");
228  fEventTree->Branch("Nactivefembs",&fNactivefembs,"Nactivefembs[6]/I");
229  fEventTree->Branch("subrun", &subrun,"surbrun/I");
230  fEventTree->Branch("T0_values",&T0_values);
231  fEventTree->Branch("xprojectedlen",&xprojectedlen);
232  fEventTree->Branch("trackthetaxz",&trackthetaxz);
233  fEventTree->Branch("trackthetayz",&trackthetayz);
234  fEventTree->Branch("trkstartx",&trkstartx);
235  fEventTree->Branch("trkstarty",&trkstarty);
236  fEventTree->Branch("trkstartz",&trkstartz);
237  fEventTree->Branch("trkendx",&trkendx);
238  fEventTree->Branch("trkendy",&trkendy);
239  fEventTree->Branch("trkendz",&trkendz);
240  fEventTree->Branch("trkendx_crt2",&trkendx_crt2);
241  fEventTree->Branch("trkstartx_crt2",&trkstartx_crt2);
242  fEventTree->Branch("crtreco_x0",&crtreco_x0);
243  fEventTree->Branch("crtreco_x1",&crtreco_x1);
244  fEventTree->Branch("crtreco_y0",&crtreco_y0);
245  fEventTree->Branch("crtreco_y1",&crtreco_y1);
246  fEventTree->Branch("crtreco_z0",&crtreco_z0);
247  fEventTree->Branch("crtreco_z1",&crtreco_z1);
248  fEventTree->Branch("crt2tickoffset",&crt2tickoffset);
249  fEventTree->Branch("trklen",&trklen);
250  fEventTree->Branch("TrkID",&TrkID);
251  fEventTree->Branch("tot_trks",&tot_trks);
252  fEventTree->Branch("hit_peakT0",&hit_peakT0);
253  fEventTree->Branch("hit_tpc0",&hit_tpc0);
254  fEventTree->Branch("hit_wire0",&hit_wire0);
255  fEventTree->Branch("hit_rms0",&hit_rms0);
256  fEventTree->Branch("hit_deltaT",&hit_deltaT);
257  fEventTree->Branch("trkhitx0",&trkhitx0);
258  fEventTree->Branch("trkhity0",&trkhity0);
259  fEventTree->Branch("trkhitz0",&trkhitz0);
260  fEventTree->Branch("trkdq_int0",&trkdq_int0);
261  fEventTree->Branch("trkdq_amp0",&trkdq_amp0);
262  fEventTree->Branch("hit_peakT1",&hit_peakT1);
263  fEventTree->Branch("hit_tpc1",&hit_tpc1);
264  fEventTree->Branch("hit_wire1",&hit_wire1);
265  fEventTree->Branch("hit_rms1",&hit_rms1);
266  fEventTree->Branch("trkhitx1",&trkhitx1);
267  fEventTree->Branch("trkhity1",&trkhity1);
268  fEventTree->Branch("trkhitz1",&trkhitz1);
269  fEventTree->Branch("trkdq_int1",&trkdq_int1);
270  fEventTree->Branch("trkdq_amp1",&trkdq_amp1);
271  fEventTree->Branch("hit_peakT2",&hit_peakT2);
272  fEventTree->Branch("hit_tpc2",&hit_tpc2);
273  fEventTree->Branch("hit_wire2",&hit_wire2);
274  fEventTree->Branch("hit_rms2",&hit_rms2);
275  fEventTree->Branch("hit_rmsraw2",&hit_rmsraw2);
276  fEventTree->Branch("hit_peakTraw2",&hit_peakTraw2);
277  fEventTree->Branch("hit_signal",&hit_signal,"hit_signal[10][5000][60]/F");
278  fEventTree->Branch("hit_rms_true",&hit_rms_true);
279  fEventTree->Branch("trkhitx2",&trkhitx2);
280  fEventTree->Branch("trkhity2",&trkhity2);
281  fEventTree->Branch("trkhitz2",&trkhitz2);
282  fEventTree->Branch("trkhitz_wire2",&trkhitz_wire2);
283  fEventTree->Branch("trkdq_int2",&trkdq_int2);
284  fEventTree->Branch("trkdq_amp2",&trkdq_amp2);
285  fEventTree->Branch("trkstartcosxyz",&trkstartcosxyz);
286  fEventTree->Branch("trkendcosxyz",&trkendcosxyz);
287  fEventTree->Branch("multiplicity2",&multiplicity2);
288  fEventTree->Branch("goodnessoffit2",&goodnessoffit2);
289  // fEventTree->Branch("t0crt1",&t0crt1);
290  fEventTree->Branch("t0crt2",&t0crt2);
291  }
292 
293  //========================================================================
295 
296  }
297 
298  //========================================================================
300  mf::LogInfo("hitrmsrawdigits")<<"begin run..."<<std::endl;
301  }
302  //========================================================================
303 
304  //========================================================================
305 
306  //========================================================================
307 
308  void hitrmsrawdigits::analyze( const art::Event& evt){//analyze
309  reset();
310  std::cout<<"raw producer module label "<<fRawProducerLabel<<std::endl;
311  // art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
313  //Detector properties service
314  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService>()->DataFor(evt);
315 
316  std::vector<art::Ptr<recob::Track> > tracklist;
317  auto trackListHandle = evt.getHandle< std::vector<recob::Track> >("pandoraTrack");
318  if (trackListHandle){
319  art::fill_ptr_vector(tracklist, trackListHandle);
320  }
321  else return;
322 
323 
324  std::vector<art::Ptr<recob::PFParticle> > pfplist;
325  auto PFPListHandle = evt.getHandle< std::vector<recob::PFParticle> >("pandora");
326  if (PFPListHandle) art::fill_ptr_vector(pfplist, PFPListHandle);
327 
328  std::vector<art::Ptr<recob::Hit>> hitlist;
329  auto hitListHandle = evt.getHandle< std::vector<recob::Hit> >(fHitsModuleLabel); // to get information about the hits
330  if (hitListHandle) art::fill_ptr_vector(hitlist, hitListHandle);
331 
332  art::FindManyP<recob::Hit, recob::TrackHitMeta> fmthm(trackListHandle, evt, fTrackModuleLabel); // to associate tracks and hits
333  art::FindManyP<anab::T0> trk_t0_assn_v(PFPListHandle, evt ,"pandora");
334  art::FindManyP<recob::PFParticle> pfp_trk_assn(trackListHandle,evt,"pandoraTrack");
335  art::FindManyP<anab::T0> fmT0(trackListHandle, evt ,"pmtrack");
336 
337  //RawDigits code -- old hardcoded label: tpcrawdecoder:daq
338  std::vector<art::Ptr<raw::RawDigit> > rawlist;
339  auto rawListHandle = evt.getHandle< std::vector<raw::RawDigit> >(fRawProducerLabel);
340  if (rawListHandle) art::fill_ptr_vector(rawlist, rawListHandle);
341 
342  // old hardcoded label: wclsdatanfsp:gauss
343  std::vector<art::Ptr<recob::Wire> > wirelist;
344  auto wireListHandle = evt.getHandle< std::vector<recob::Wire> >(fWireProducerLabel);
345  if (wireListHandle) art::fill_ptr_vector(wirelist, wireListHandle);
346 
347  ///RawDigits
348  std::cout<<"rawlist size, wirelist size"<<rawlist.size()<<" "<<wirelist.size()<<std::endl;
349  std::vector<const sim::SimChannel*> fSimChannels;
350  try{
351  evt.getView("largeant", fSimChannels);
352  }catch (art::Exception const&e){
353  }
354 
355  //Get 2-CRT T0
356  art::FindManyP<anab::T0> fmt0crt2(trackListHandle, evt, "crtreco");
357  art::FindManyP<anab::CosmicTag> fmctcrt2(trackListHandle, evt, "crtreco");
358 
359  //Get 1-CRT T0
360  // art::FindManyP<anab::T0> fmt0crt1(trackListHandle, evt, "crttag");
361  // art::FindManyP<anab::CosmicTag> fmctcrt1(trackListHandle, evt, "crttag");
362 
363  run = evt.run();
364  subrun = evt.subRun();
365  event = evt.id().event();
366  art::Timestamp ts = evt.time();
367  TTimeStamp tts(ts.timeHigh(), ts.timeLow());
368  evttime=tts.AsDouble();
369 
370 
371 
372  // Get number of active fembs
373  if(!evt.isRealData()){
374  for(int k=0; k < 6; k++)
375  fNactivefembs[k] = 20;
376  }
377  else{
378  for(int k=0; k < 6; k++)
380  }
381 
382 
383  //defining the 1D temporary storage vectors
384  std::vector< float> peakT_0; std::vector< float> peakT_1; std::vector< float> peakT_2; std::vector<float> hit_peakTrawb;
385  std::vector< int > tpc_0; std::vector< int > tpc_1; std::vector<int> tpc_2;
386  std::vector< int > wire_0; std::vector< int > wire_1; std::vector<int> wire_2;
387  std::vector< float > int_0; std::vector< float > int_1; std::vector<float> int_2;
388  std::vector< float > amp_0; std::vector< float > amp_1; std::vector<float> amp_2;
389  std::vector<float> rms_0; std::vector<float> rms_1; std::vector<float> rms_2, rms_raw2;
390  std::vector<float> hitx_0; std::vector<float> hitx_1; std::vector<float> hitx_2;
391  std::vector<float> hity_0; std::vector<float> hity_1; std::vector<float> hity_2;
392  std::vector<float> hitz_0; std::vector<float> hitz_1; std::vector<float> hitz_2;
393  std::vector<float> hitz_wire2; std::vector<float> startcosxyz; std::vector<float> endcosxyz;
394  std::vector<float> dT_buffer; std::vector<float> gof, multi,truerms;
395 
396  float max_value;
397  float min_value;
398  int ntrks=0;
399  size_t NTracks = tracklist.size();
400  for(size_t i=0; i<NTracks;++i){
401  double xoffset=0.0;
402  int nhits=0;
403  //clearing the 1D temporary storage vectors
404  peakT_0.clear(); peakT_1.clear(); peakT_2.clear(); hit_peakTrawb.clear();
405  tpc_0.clear(); tpc_1.clear(); tpc_2.clear();
406  wire_0.clear(); wire_1.clear(); wire_2.clear();
407  int_0.clear(); int_1.clear() ; int_2.clear();
408  amp_0.clear(); amp_1.clear() ; amp_2.clear();
409  rms_0.clear(); rms_1.clear(); rms_2.clear();rms_raw2.clear();
410  hitx_0.clear(); hitx_1.clear(); hitx_2.clear();
411  hity_0.clear(); hity_1.clear(); hity_2.clear();
412  hitz_0.clear(); hitz_1.clear(); hitz_2.clear();
413  hitz_wire2.clear(); startcosxyz.clear();endcosxyz.clear();
414  dT_buffer.clear(); gof.clear(); multi.clear();
415  truerms.clear();
416 
417 
418  art::Ptr<recob::Track> ptrack(trackListHandle, i);
419 
420  ///this block just saves the t0 values while I include entries with no T0s as well also this saves T0 coming from pandroaTrack alg only
421  double t_zero=-999999;
422  max_value=0.0;
423  min_value=0.0;
424  /* if(fTrackModuleLabel=="pandoraTrack"){
425  std::vector<art::Ptr<recob::PFParticle>> pfps=pfp_trk_assn.at(i);
426  if(!pfps.size()) continue;
427  std::vector<art::Ptr<anab::T0>> t0s=trk_t0_assn_v.at(pfps[0].key());
428  // if(!t0s.size()) continue;
429  //auto t0 = t0s.at(0);
430  // double t_zero=t0->Time();
431  if(t0s.size()==0) continue;
432  if(t0s.size()){
433  auto t0=t0s.at(0);
434  t_zero=t0->Time();
435  }
436  }
437 
438  if(fTrackModuleLabel=="pmtrack"){
439  std::vector<art::Ptr<anab::T0>> T0s=fmT0.at(i);
440  if(T0s.size()==0)
441  continue;
442  }
443  */
444  // all_trks++;
445  const recob::Track& track = *ptrack;
446 
447 
448  double this_t0crt2=-DBL_MAX;
449  //double this_t0crt1=-DBL_MAX;
450  double this_crt2x0 = -DBL_MAX;
451  double this_crt2x1 = -DBL_MAX;
452  double this_crt2y0 = -DBL_MAX;
453  double this_crt2y1 = -DBL_MAX;
454  double this_crt2z0 = -DBL_MAX;
455  double this_crt2z1 = -DBL_MAX;
456 
457  bool test1=true;
458  // bool test2=true;
459  if(fmt0crt2.isValid()){
460  auto const& vt0crt2 = fmt0crt2.at(i);
461  if (!vt0crt2.empty()){
462  this_t0crt2 = vt0crt2[0]->Time();
463  test1=false;
464  }
465 
466  }
467  if(test1) continue;
468  if (fmctcrt2.isValid()){
469  auto const& vctcrt2 = fmctcrt2.at(i);
470  if (!vctcrt2.empty()){
471  this_crt2x0 = vctcrt2[0]->EndPoint1()[0];
472  this_crt2x1 = vctcrt2[0]->EndPoint2()[0];
473  this_crt2y0 = vctcrt2[0]->EndPoint1()[1];
474  this_crt2y1 = vctcrt2[0]->EndPoint2()[1];
475  this_crt2z0 = vctcrt2[0]->EndPoint1()[2];
476  this_crt2z1 = vctcrt2[0]->EndPoint2()[2];
477 
478  }
479  }
480 
481 
482 
483 
484  /* if (fmt0crt1.isValid()){
485  auto const& vt0crt1 = fmt0crt1.at(i);
486  if (!vt0crt1.empty()){
487  this_t0crt1 = vt0crt1[0]->Time();
488  test2=false;
489  }
490  }*/
491  // if(test1 && test2) continue;
492 
493 
494 
495 
496  auto pos = track.Vertex();
497  auto dir_start = track.VertexDirection();
498  auto dir_end = track.EndDirection();
499  auto end = track.End();
500  double theta_xz = std::atan2(dir_start.X(), dir_start.Z());
501  double theta_yz = std::atan2(dir_start.Y(), dir_start.Z());
502 
503  int planenum=999;
504  float xpos=-9999;
505  float ypos=-9999;
506  float zpos=-9999;
507  auto allHits=fmthm.at(i);
508  double ticksoffset=0;
509 
510  // if (this_t0crt2 > -DBL_MAX) ticksoffset = this_t0crt2/500.+detProp.GetXTicksOffset(allHits[0]->WireID());
511  if (this_t0crt2 > -DBL_MAX) ticksoffset = this_t0crt2/500.+detProp.GetXTicksOffset (allHits[0]->WireID().Plane, allHits[0]->WireID().TPC, allHits[0]->WireID().Cryostat);
512  // else if (this_t0crt1 > -DBL_MAX) ticksoffset = this_t0crt1/500.+detProp.GetXTicksOffset(allHits[0]->WireID());
513  xoffset = detProp.ConvertTicksToX(ticksoffset,allHits[0]->WireID());
514  std::cout<<"tickoffset , x offset "<<ticksoffset<<" "<<xoffset<<" default term "<<detProp.GetXTicksOffset (allHits[10]->WireID().Plane, allHits[10]->WireID().TPC, allHits[10]->WireID().Cryostat)<<std::endl;
515 
516 
517 
518  //hits and calorimetry loop
519  if(fmthm.isValid()){
520  auto vhit=fmthm.at(i);
521  auto vmeta=fmthm.data(i);
522  for (size_t ii = 0; ii<vhit.size(); ++ii){ //loop over all meta data hit
523  bool fBadhit = false;
524  if (vmeta[ii]->Index() == static_cast<unsigned int>(std::numeric_limits<int>::max())){
525  fBadhit = true;
526  //cout<<"fBadHit"<<fBadhit<<endl;
527  continue;
528  }
529  if (vmeta[ii]->Index()>=tracklist[i]->NumberTrajectoryPoints()){
530  throw cet::exception("Calorimetry_module.cc") << "Requested track trajectory index "<<vmeta[ii]->Index()<<" exceeds the total number of trajectory points "<<tracklist[i]->NumberTrajectoryPoints()<<" for track index "<<i<<". Something is wrong with the track reconstruction. Please contact tjyang@fnal.gov!!";
531  }
532  if (!tracklist[i]->HasValidPoint(vmeta[ii]->Index())){
533  fBadhit = true;
534  // cout<<"had valid point "<<fBadhit<<endl;
535  continue;
536  }
537 
538  auto loc = tracklist[i]->LocationAtPoint(vmeta[ii]->Index());
539  xpos=loc.X();
540  ypos=loc.Y();
541  zpos=loc.Z();
542  // cout<<"x, y, z "<<xpos<<" "<<ypos<<" "<<zpos<<endl;
543  // cout<<"BadHit"<<fBadhit<<endl;
544  if (fBadhit) continue; //HY::If BAD hit, skip this hit and go next
545  if (zpos<-100) continue; //hit not on track
546  planenum=vhit[ii]->WireID().Plane;
547  if(planenum==0){
548  peakT_0.push_back(vhit[ii]->PeakTime());
549  tpc_0.push_back(vhit[ii]->WireID().TPC);
550  wire_0.push_back(vhit[ii]->WireID().Wire);
551  int_0.push_back(vhit[ii]->Integral());
552  amp_0.push_back(vhit[ii]->PeakAmplitude());
553  rms_0.push_back(vhit[ii]->RMS());
554  hitx_0.push_back(xpos);
555  hity_0.push_back(ypos);
556  hitz_0.push_back(zpos);
557  }//planenum 0
558  if(planenum==1){
559  peakT_1.push_back(vhit[ii]->PeakTime());
560  tpc_1.push_back(vhit[ii]->WireID().TPC);
561  wire_1.push_back(vhit[ii]->WireID().Wire);
562  int_1.push_back(vhit[ii]->Integral());
563  amp_1.push_back(vhit[ii]->PeakAmplitude());
564  rms_1.push_back(vhit[ii]->RMS());
565  hitx_1.push_back(xpos);
566  hity_1.push_back(ypos);
567  hitz_1.push_back(zpos);
568  }//planenum 1
569  if(planenum==2){
570  // Remove all hits within 30tt of other hits by looping over all other hits:
571  bool removehit=false;
572  for (size_t i2=0; i2<hitlist.size(); ++i2) {
573  if (vhit[ii].key() == hitlist[i2].key()) continue;
574  if (vhit[ii]->Channel()!=hitlist[i2]->Channel()) continue;
575  if ((vhit[ii]->PeakTime()+30<hitlist[i2]->PeakTime()-30) || (vhit[ii]->PeakTime()-30>hitlist[i2]->PeakTime()+30)) continue;
576  removehit=true;
577  break;
578  }
579  if (removehit) continue;
580  // Normal procedure resumes (i.e. keep this if keeping all hits):
581  nhits++;
582  peakT_2.push_back(vhit[ii]->PeakTime()-this_t0crt2/500.0);
583  tpc_2.push_back(vhit[ii]->WireID().TPC);
584  wire_2.push_back(vhit[ii]->WireID().Wire);
585  int_2.push_back(vhit[ii]->Integral());
586  amp_2.push_back(vhit[ii]->PeakAmplitude());
587  rms_2.push_back(vhit[ii]->RMS());
588  dT_buffer.push_back(vhit[ii]->EndTick()-vhit[ii]->StartTick());
589  hitx_2.push_back(xpos);
590  hity_2.push_back(ypos);
591  hitz_2.push_back(zpos);
592  gof.push_back(vhit[ii]->GoodnessOfFit());
593  multi.push_back(vhit[ii]->Multiplicity());
594  double xyzStart[3];
595  double xyzEnd[3];
596  unsigned int wireno=vhit[ii]->WireID().Wire;
597  fGeometry->WireEndPoints(0,vhit[ii]->WireID().TPC,2,wireno, xyzStart, xyzEnd);
598  hitz_wire2.push_back(xyzStart[2]);
599  double truermsb=-1;
600 
601  //section for using tuth information
602  unsigned int t0=vhit[ii]->PeakTime()-3*(vhit[ii]->RMS());
603  if(t0<0) t0=0;
604  unsigned int t1=vhit[ii]->PeakTime()+3*(vhit[ii]->RMS());
605  if(t1>6000) t1=6000-1;
606 
607  ////RawDigits section
608  double hit_rms=-9999;
609  std::cout<<"wirelist size, rawlist size "<<wirelist.size()<<" "<<rawlist.size()<<std::endl;
610  double hit_t = -1;
611  for (unsigned int ich = 0; ich < (rawlist.empty()?wirelist.size():rawlist.size()); ++ich){
612  std::vector<float> inputsignal(fWaveformSize);
613 
614  if (!wirelist.empty() && evt.isRealData()){
615  const auto &wire = wirelist[ich];
616  if(wirelist[ich]->Channel()!=vhit[ii]->Channel()) continue;
617  //art::Ptr<recob::Wire> wire(wireHandle, wireIter);
618  const auto &signal = wire->Signal();
619  double hit_pk = -1;
620 
621  for (size_t itck = t0; itck <inputsignal.size(); ++itck){
622  if(itck>t1) continue;
623  inputsignal[itck] = signal[itck];
624  if (inputsignal[itck]>hit_pk){
625  hit_pk = inputsignal[itck];
626  hit_t = itck;
627  }
628  }//finding peak signal
629  std::cout<<"hitpeak time, hitraw peak time "<<vhit[ii]->PeakTime()<<" "<<hit_t<<std::endl;
630  int hitindx=0;
631  for(size_t it1=hit_t-30;it1<hit_t+30;it1++){
632  if(it1<1||it1>5999||it1>inputsignal.size()) continue;
633  hit_signal[ntrks][nhits-1][hitindx]=signal[it1];
634  hitindx++;
635  }//storing signal for peakT-30 to peakT+30 ticks
636  double hit_ch = 0;
637  double hit_fwhh = 0;
638  double mean_t = 0;
639  double mean_t2 = 0;
640  for (size_t itck = t0; itck < inputsignal.size(); ++itck){
641  if(itck>t1) continue;
642  if (inputsignal[itck]>=0.5*hit_pk){
643  ++hit_fwhh;
644  }
645  if (inputsignal[itck]>=0.1*hit_pk){
646  hit_ch += inputsignal[itck];
647  mean_t += itck*(inputsignal[itck]);
648  mean_t2 += itck*itck*(inputsignal[itck]);
649  }
650  }//itick loop
651  mean_t/=hit_ch;
652  mean_t2/=hit_ch;
653  hit_rms = sqrt(mean_t2-mean_t*mean_t);
654  }
655  else if (!rawlist.empty()){
656  const auto & digitVec = rawlist[ich];
657  if(rawlist[ich]->Channel()!=vhit[ii]->Channel()) continue;
658  std::vector<short> rawadc(fWaveformSize);
659  raw::Uncompress(digitVec->ADCs(), rawadc, digitVec->GetPedestal(), digitVec->Compression());
660  double hit_pk = -1;
661  // double hit_t = -1;
662  std::vector<double> signalbuffer;
663  signalbuffer.clear();
664  //calculating the median of signals
665  for(size_t it=0;it<rawadc.size();it++){
666  signalbuffer.push_back(rawadc[it]);
667  }
668  double med_signal=TMath::Median(signalbuffer.size(),&signalbuffer[0]);//median signal for a channel
669  std::cout<<"median , pedestal , rawadc size "<<med_signal<<" "<<digitVec->GetPedestal()<<" "<<rawadc.size()<<std::endl;
670  ////////////////////////////////
671 
672  for (size_t itck = t0; itck <rawadc.size(); ++itck){
673  if(itck>t1) continue;
674  // inputsignal[itck] = rawadc[itck] - digitVec->GetPedestal();
675  inputsignal[itck] = rawadc[itck] - med_signal;
676  if (inputsignal[itck]>hit_pk){
677  hit_pk = inputsignal[itck];
678  hit_t = itck;
679  }
680  }//itick loop
681  std::cout<<"hitpeak time, hitraw peak time "<<vhit[ii]->PeakTime()<<" "<<hit_t<<std::endl;
682  int hitindex=0;
683  for(size_t it1=hit_t-30;it1<hit_t+30;it1++){
684  if(it1<1|| it1>5999 || it1>rawadc.size()) continue;
685  // hit_signal[ntrks][nhits-1][hitindex]=rawadc[it1]-digitVec->GetPedestal();
686  hit_signal[ntrks][nhits-1][hitindex]=rawadc[it1]-med_signal;
687  hitindex++;
688  }
689 
690  double hit_ch = 0;
691  double hit_fwhh = 0;
692  double mean_t = 0;
693  double mean_t2 = 0;
694  for (size_t itck = t0; itck < inputsignal.size(); ++itck){
695  if(itck>t1) continue;
696  // inputsignal[itck] = rawadc[itck] - digitVec->GetPedestal();
697  inputsignal[itck] = rawadc[itck] - med_signal;
698  if (inputsignal[itck]>=0.5*hit_pk){
699  ++hit_fwhh;
700  }
701  if (inputsignal[itck]>=0.1*hit_pk){
702  hit_ch += inputsignal[itck];
703  mean_t += itck*(inputsignal[itck]);
704  mean_t2 += itck*itck*(inputsignal[itck]);
705  }
706  }//itick loop
707  mean_t/=hit_ch;
708  mean_t2/=hit_ch;
709  hit_rms = sqrt(mean_t2-mean_t*mean_t);
710  }
711  }//rawdigit channel loop
712  rms_raw2.push_back(hit_rms);
713  hit_peakTrawb.push_back(hit_t);
714  //////RawDifits section
715  if(!evt.isRealData()){
716  const sim::SimChannel* chan=0;
717  for(size_t sc=0;sc<fSimChannels.size();++sc){
718  if(fSimChannels[sc]->Channel()==vhit[ii]->Channel()) chan=fSimChannels[sc];
719  }
720  if(chan){
721  const auto &tdcidemap = chan->TDCIDEMap();
722  double hit_ch=0;
723  double mean_t=0;
724  double mean_t2=0;
725  double hit_rmsvalue=0;
726  for(auto mapitr = tdcidemap.begin(); mapitr != tdcidemap.end(); mapitr++){
727  if(!((*mapitr).first>t0 && (*mapitr).first<=t1)) continue;
728  // loop over the vector of IDE objects.
729  int hit_nelec=0;
730  const std::vector<sim::IDE> idevec = (*mapitr).second;
731  for(size_t iv = 0; iv < idevec.size(); ++iv){
732  hit_nelec+=idevec[iv].numElectrons;
733  }//iv loop
734  hit_ch+=hit_nelec;
735  // std::cout<<"hit_ch "<<hit_ch<<std::endl;
736  int j=(*mapitr).first;
737  mean_t+=double(j)*double(hit_nelec);
738  double jterm=double(j)*double(j)*double(hit_nelec);
739  mean_t2=mean_t2+jterm;
740  }//mapitr loop
741  mean_t/=hit_ch;
742  mean_t2/=hit_ch;
743  hit_rmsvalue=sqrt(mean_t2-mean_t*mean_t);
744  // std::cout<<"hit rms "<<hit_rmsvalue<<" reco rms "<<vhit[ii]->RMS()<<std::endl;
745  truermsb=hit_rmsvalue;
746  }//if(chan)
747  }//!evt.IsRealData
748  truerms.push_back(truermsb);
749 
750 
751  }//planenum 2
752  }//loop over vhit
753  }//fmthm valid
754  //hits and calorimetry loop
755  if(peakT_2.size()<10) continue;
756  max_value=*std::max_element(peakT_2.begin(),peakT_2.end());
757  min_value=*std::min_element(peakT_2.begin(),peakT_2.end());
758  // if(max_value-min_value<4300) continue;
759  std::cout<<max_value<<" "<<min_value<<std::endl;
760  ntrks++;
761  hit_peakT0.push_back(peakT_0);
762  hit_tpc0.push_back(tpc_0);
763  hit_wire0.push_back(wire_0);
764  hit_rms0.push_back(rms_0);
765  trkhitx0.push_back(hitx_0);
766  trkhity0.push_back(hity_0);
767  trkhitz0.push_back(hitz_0);
768  trkdq_amp0.push_back(amp_0);
769  trkdq_int0.push_back(int_0);
770 
771  hit_peakT1.push_back(peakT_1);
772  hit_tpc1.push_back(tpc_1);
773  hit_wire1.push_back(wire_1);
774  hit_rms1.push_back(rms_1);
775  trkhitx1.push_back(hitx_1);
776  trkhity1.push_back(hity_1);
777  trkhitz1.push_back(hitz_1);
778  trkdq_amp1.push_back(amp_1);
779  trkdq_int1.push_back(int_1);
780 
781  hit_peakT2.push_back(peakT_2);
782  hit_peakTraw2.push_back(hit_peakTrawb);
783  hit_tpc2.push_back(tpc_2);
784  hit_wire2.push_back(wire_2);
785  hit_rms2.push_back(rms_2);
786  hit_rmsraw2.push_back(rms_raw2);
787  hit_rms_true.push_back(truerms);
788  trkhitx2.push_back(hitx_2);
789  trkhity2.push_back(hity_2);
790  trkhitz2.push_back(hitz_2);
791  trkhitz_wire2.push_back(hitz_wire2);
792  trkdq_amp2.push_back(amp_2);
793  trkdq_int2.push_back(int_2);
794  hit_deltaT.push_back(dT_buffer);
795  multiplicity2.push_back(multi);
796  goodnessoffit2.push_back(gof);
797 
798 
799 
800  xprojectedlen.push_back(TMath::Abs(end.X()-pos.X()));
801  trackthetaxz.push_back(theta_xz);
802  trackthetayz.push_back(theta_yz);
803  trkstartx.push_back(pos.X());
804  trkstartx_crt2.push_back(pos.X()-xoffset);
805  trkstarty.push_back(pos.Y());
806  trkstartz.push_back(pos.Z());
807  startcosxyz.push_back(dir_start.X());
808  startcosxyz.push_back(dir_start.Y());
809  startcosxyz.push_back(dir_start.Z());
810  endcosxyz.push_back(dir_end.X());
811  endcosxyz.push_back(dir_end.Y());
812  endcosxyz.push_back(dir_end.Z());
813 
814  trkstartcosxyz.push_back(startcosxyz);
815  trkendcosxyz.push_back(endcosxyz);
816  trkendx.push_back(end.X());
817  trkendx_crt2.push_back(end.X()-xoffset);
818  crtreco_x0.push_back(this_crt2x0);
819  crtreco_x1.push_back(this_crt2x1);
820  crtreco_y0.push_back(this_crt2y0);
821  crtreco_y1.push_back(this_crt2y1);
822  crtreco_z0.push_back(this_crt2z0);
823  crtreco_z1.push_back(this_crt2z1);
824 
825  trkendy.push_back(end.Y());
826  trkendz.push_back(end.Z());
827  trklen.push_back(track.Length());
828  TrkID.push_back(track.ID());
829  T0_values.push_back(t_zero);
830  crt2tickoffset.push_back(ticksoffset);
831  t0crt2.push_back(this_t0crt2/500.0);
832  } // loop over trks...
833  tot_trks.push_back(ntrks);
834  fEventTree->Fill();
835  } // end of analyze function
836 
837  /////////////////// Defintion of reset function ///////////
839  run = -9999;
840  subrun = -9999;
841  event = -9999;
842  evttime = -9999;
843  //all_trks = -9999;
844  for(int k=0; k < 6; k++)
845  fNactivefembs[k] = -9999;
846  trackthetaxz.clear();
847  trackthetayz.clear();
848  trkstartx.clear();
849  trkstartx_crt2.clear();
850  trkendx_crt2.clear();
851  trkstarty.clear();
852  trkstartz.clear();
853  trkstartcosxyz.clear();
854  trkendcosxyz.clear();
855  trkendx.clear();
856  trkendy.clear();
857  trkendz.clear();
858  trklen.clear();
859  TrkID.clear();
860  tot_trks.clear();
861  T0_values.clear();
862  xprojectedlen.clear();
863  hit_peakT0.clear();
864  hit_tpc0.clear();
865  hit_wire0.clear();
866  trkhitx0.clear();
867  trkhity0.clear();
868  trkhitz0.clear();
869  trkdq_int0.clear();
870  trkdq_amp0.clear();
871  hit_rms0.clear();
872  hit_peakT1.clear();
873  hit_tpc1.clear();
874  hit_wire1.clear();
875  trkhitx1.clear();
876  trkhity1.clear();
877  trkhitz1.clear();
878  trkdq_int1.clear();
879  trkdq_amp1.clear();
880  hit_rms1.clear();
881  hit_peakT2.clear();
882  hit_tpc2.clear();
883  hit_wire2.clear();
884  trkhitx2.clear();
885  trkhity2.clear();
886  trkhitz2.clear();
887  trkhitz_wire2.clear();
888  trkdq_int2.clear();
889  trkdq_amp2.clear();
890  hit_rms2.clear();
891  hit_peakTraw2.clear();
892  hit_rmsraw2.clear();
893  hit_deltaT.clear();
894  multiplicity2.clear();
895  goodnessoffit2.clear();
896  hit_rms_true.clear();
897  crtreco_x0.clear();
898  crtreco_x1.clear();
899  crtreco_y0.clear();
900  crtreco_y1.clear();
901  crtreco_z0.clear();
902  crtreco_z1.clear();
903  crt2tickoffset.clear();
904 
905  // t0crt1.clear();
906  t0crt2.clear();
907  for(int i=0; i<10; i++){
908  for(int j=0; j<5000; j++){
909  for(int k=0;k<60;k++){
910  hit_signal[i][j][k]=-99999;;
911  }
912  }
913  }
914  }
915  //////////////////////// End of definition ///////////////
916 
918 }
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
code to link reconstructed objects back to the MC truth information
std::vector< std::vector< float > > trkhitz_wire2
constexpr std::uint32_t timeLow() const
Definition: Timestamp.h:29
std::vector< std::vector< float > > hit_rms1
std::vector< std::vector< int > > hit_tpc0
Energy deposited on a readout channel by simulated tracks.
Definition: SimChannel.h:140
std::vector< std::vector< int > > hit_wire1
hitrmsrawdigits(fhicl::ParameterSet const &pset)
TH3F * xpos
Definition: doAna.cpp:29
std::vector< std::vector< float > > trkdq_int1
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
std::string string
Definition: nybbler.cc:12
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
std::vector< std::vector< float > > hit_peakT0
std::vector< std::vector< float > > trkhitx2
std::vector< std::vector< float > > trkhitz2
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:34
std::vector< std::vector< float > > hit_deltaT
Class to keep data related to recob::Hit associated with recob::Track.
Vector_t VertexDirection() const
Definition: Track.h:132
std::vector< std::vector< float > > hit_peakT1
std::vector< std::vector< int > > hit_wire2
STL namespace.
void analyze(const art::Event &evt)
std::vector< std::vector< float > > hit_peakTraw2
unsigned int Index
std::vector< float > xprojectedlen
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
std::vector< std::vector< float > > trkhity1
Particle class.
object containing MC flux information
art framework interface to geometry description
Definition: Run.h:17
TH3F * ypos
Definition: doAna.cpp:30
std::vector< std::vector< float > > hit_rms_true
std::vector< std::vector< float > > trkdq_amp0
void beginRun(const art::Run &run)
std::vector< std::vector< int > > hit_wire0
std::vector< std::vector< float > > trkhitx1
bool isRealData() const
TH3F * zpos
Definition: doAna.cpp:31
const double e
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:167
Timestamp time() const
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
def key(type, name=None)
Definition: graph.py:13
std::vector< std::vector< int > > hit_tpc2
std::vector< std::vector< float > > trkdq_int2
Point_t const & Vertex() const
Definition: Track.h:124
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:78
std::vector< float > trkstartx_crt2
RunNumber_t run() const
Definition: DataViewImpl.cc:71
static int max(int a, int b)
Description of geometry of one entire detector.
Definition of data types for geometry description.
std::vector< std::vector< float > > trkendcosxyz
std::vector< std::vector< float > > multiplicity2
std::vector< std::vector< float > > trkstartcosxyz
std::vector< std::vector< float > > hit_peakT2
std::vector< std::vector< float > > trkdq_amp2
std::vector< std::vector< float > > trkdq_amp1
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
std::vector< std::vector< float > > hit_rms0
std::vector< double > crt2tickoffset
int ID() const
Definition: Track.h:198
Declaration of signal hit object.
int GetNActiveFembsForAPA(art::Event const &evt, int apa) const
Get number of active fembs in an APA.
Vector_t EndDirection() const
Definition: Track.h:133
std::vector< std::vector< float > > hit_rmsraw2
geo::GeometryCore const * fGeometry
std::vector< std::vector< float > > trkhitx0
Provides recob::Track data product.
void WireEndPoints(geo::WireID const &wireid, double *xyzStart, double *xyzEnd) const
Fills two arrays with the coordinates of the wire end points.
EventNumber_t event() const
Definition: EventID.h:116
Point_t const & End() const
Definition: Track.h:125
Declaration of basic channel signal object.
TDCIDEs_t const & TDCIDEMap() const
Returns all the deposited energy information as stored.
Definition: SimChannel.h:328
std::size_t getView(std::string const &moduleLabel, std::string const &productInstanceName, std::string const &processName, std::vector< ELEMENT const * > &result) const
Definition: DataViewImpl.h:500
std::vector< std::vector< float > > trkhity2
std::vector< std::vector< float > > hit_rms2
TCEvent evt
Definition: DataStructs.cxx:7
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< std::vector< float > > trkhitz1
void Uncompress(const std::vector< short > &adc, std::vector< short > &uncompressed, raw::Compress_t compress)
Uncompresses a raw data buffer.
Definition: raw.cxx:776
recob::tracking::Plane Plane
Definition: TrackState.h:17
int bool
Definition: qglobal.h:345
std::vector< std::vector< float > > trkhitz0
EventID id() const
Definition: Event.cc:34
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
Definition: Track.h:49
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::vector< std::vector< int > > hit_tpc1
QTextStream & endl(QTextStream &s)
Event finding and building.
std::vector< std::vector< float > > trkdq_int0
std::vector< std::vector< float > > goodnessoffit2
std::vector< std::vector< float > > trkhity0