LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
RawWaveformDump_module.cc
Go to the documentation of this file.
1 //
3 // Module to dump raw data for ROI training
4 //
5 // mwang@fnal.gov
6 // tjyang@fnal.gov
7 // wwu@fnal.gov
8 //
9 // This version:
10 // a) uses RawDigit or WireProducer
11 // b) saves full waveform or short waveform
12 // - for short waveform, picks signal
13 // associated with largest energy
14 // deposit & saves only tdcmin/tdcmaxes
15 // that are within window extents
16 //
18 
19 // Framework includes
30 #include "cetlib_except/exception.h"
31 #include "fhiclcpp/ParameterSet.h"
33 
34 // LArSoft libraries
37 #include "larcoreobj/SimpleTypesAndConstants/RawTypes.h" // raw::ChannelID_t
41 #include "lardataobj/RawData/raw.h"
45 #include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h"
46 #include "larevt/CalibrationDBI/Interface/ChannelStatusService.h"
51 
52 #include "CLHEP/Random/RandFlat.h"
53 #include "c2numpy.h"
54 
55 #include <algorithm>
56 #include <chrono>
57 #include <cstdlib>
58 #include <map>
59 #include <random>
60 #include <set>
61 #include <string>
62 #include <utility>
63 #include <vector>
64 
65 using std::cout;
66 using std::endl;
67 using std::ofstream;
68 using std::string;
69 
70 struct WireSigInfo {
71  int pdgcode;
72  std::string genlab;
73  std::string procid;
74  unsigned int tdcmin;
75  unsigned int tdcmax;
76  int numel;
77  double edep;
78 };
79 
80 namespace nnet {
81  class RawWaveformDump;
82 }
83 
85 
86 public:
87  explicit RawWaveformDump(fhicl::ParameterSet const& p);
88 
89  // Plugins should not be copied or assigned.
90  RawWaveformDump(RawWaveformDump const&) = delete;
91  RawWaveformDump(RawWaveformDump&&) = delete;
92  RawWaveformDump& operator=(RawWaveformDump const&) = delete;
93  RawWaveformDump& operator=(RawWaveformDump&&) = delete;
94 
95  // Required functions.
96  void analyze(art::Event const& e) override;
97 
98  //void reconfigure(fhicl::ParameterSet const & p);
99 
100  void beginJob() override;
101  void endJob() override;
102 
103 private:
105 
107  std::string fSimChannelLabel;
108  std::string fDigitModuleLabel;
109  std::string fWireProducerLabel;
111  unsigned int fShortWaveformSize;
112 
113  std::string fSelectGenLabel;
114  std::string fSelectProcID;
116  std::string fPlaneToDump;
124  geo::WireReadoutGeom const* fWireReadoutGeom{
127 
128  CLHEP::RandFlat fRandFlat;
129 
131 };
132 
133 //-----------------------------------------------------------------------
134 struct genFinder {
135 private:
136  typedef std::pair<int, std::string> track_id_to_string;
137  std::vector<track_id_to_string> track_id_map;
138  std::set<std::string> generator_names;
139  bool isSorted = false;
140 
141 public:
142  void sort_now()
143  {
144  std::sort(this->track_id_map.begin(),
145  this->track_id_map.end(),
146  [](const auto& a, const auto& b) { return (a.first < b.first); });
147  isSorted = true;
148  }
149  void add(const int& track_id, const std::string& gname)
150  {
151  this->track_id_map.push_back(std::make_pair(track_id, gname));
152  generator_names.emplace(gname);
153  isSorted = false;
154  }
155  bool has_gen(std::string gname) { return static_cast<bool>(generator_names.count(gname)); };
156  std::string get_gen(int tid)
157  {
158  if (!isSorted) { this->sort_now(); }
159  return std::lower_bound(track_id_map.begin(),
160  track_id_map.end(),
161  tid,
162  [](const auto& a, const auto& b) { return (a.first < b); })
163  ->second;
164  };
165 };
166 
167 // Create the random number generator
168 namespace {
169  std::string const instanceName = "RawWaveformDump";
170 }
171 
172 //-----------------------------------------------------------------------
174  : EDAnalyzer{p}
175  , fDumpWaveformsFileName(p.get<std::string>("DumpWaveformsFileName", "dumpwaveforms"))
176  , fSimulationProducerLabel(p.get<std::string>("SimulationProducerLabel", "larg4Main"))
177  , fSimChannelLabel(p.get<std::string>("SimChannelLabel", "elecDrift"))
178  , fDigitModuleLabel(p.get<std::string>("DigitModuleLabel", "simWire"))
179  , fWireProducerLabel(p.get<std::string>("WireProducerLabel"))
180  , fUseFullWaveform(p.get<bool>("UseFullWaveform", true))
181  , fShortWaveformSize(p.get<unsigned int>("ShortWaveformSize"))
182  , fSelectGenLabel(p.get<std::string>("SelectGenLabel", "ANY"))
183  , fSelectProcID(p.get<std::string>("SelectProcID", "ANY"))
184  , fSelectPDGCode(p.get<int>("SelectPDGCode", 0))
185  , fPlaneToDump(p.get<std::string>("PlaneToDump"))
186  , fMinParticleEnergyGeV(p.get<double>("MinParticleEnergyGeV", 0.))
187  , fMinEnergyDepositedMeV(p.get<double>("MinEnergyDepositedMeV", 0.))
188  , fMinNumberOfElectrons(p.get<int>("MinNumberOfElectrons", 600))
189  , fMaxNumberOfElectrons(p.get<int>("MaxNumberOfElectrons", -1))
190  , fSaveSignal(p.get<bool>("SaveSignal", true))
191  , fMaxNoiseChannelsPerEvent(p.get<int>("MaxNoiseChannelsPerEvent"))
192  , fCollectionPlaneLabel(p.get<std::string>("CollectionPlaneLabel"))
195  -> declareEngine(instanceName, p, "SeedForRawWaveformDump"),
196  "HepJamesRandom",
197  instanceName)}
198 {
199  if (std::getenv("CLUSTER") && std::getenv("PROCESS")) {
201  string(std::getenv("CLUSTER")) + "-" + string(std::getenv("PROCESS")) + "-";
202  }
203 
204  if (fDigitModuleLabel.empty() && fWireProducerLabel.empty()) {
205  throw cet::exception("RawWaveformDump")
206  << "Both DigitModuleLabel and WireProducerLabel are empty";
207  }
208 
209  if ((!fDigitModuleLabel.empty()) && (!fWireProducerLabel.empty())) {
210  throw cet::exception("RawWaveformDump")
211  << "Only one of DigitModuleLabel and WireProducerLabel should be set";
212  }
213 }
214 
215 //-----------------------------------------------------------------------
217 {
218  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob();
219 
225 
226  for (unsigned int i = 0; i < 5; i++) {
227  std::ostringstream name;
228 
229  name.str("");
230  name << "tid" << i;
231  c2numpy_addcolumn(&npywriter, name.str().c_str(), C2NUMPY_INT32);
232 
233  name.str("");
234  name << "pdg" << i;
235  c2numpy_addcolumn(&npywriter, name.str().c_str(), C2NUMPY_INT32);
236 
237  name.str("");
238  name << "gen" << i;
239  c2numpy_addcolumn(&npywriter, name.str().c_str(), (c2numpy_type)((int)C2NUMPY_STRING + 6));
240 
241  name.str("");
242  name << "pid" << i;
243  c2numpy_addcolumn(&npywriter, name.str().c_str(), (c2numpy_type)((int)C2NUMPY_STRING + 7));
244 
245  name.str("");
246  name << "edp" << i;
247  c2numpy_addcolumn(&npywriter, name.str().c_str(), C2NUMPY_FLOAT32);
248 
249  name.str("");
250  name << "nel" << i;
251  c2numpy_addcolumn(&npywriter, name.str().c_str(), C2NUMPY_UINT32);
252 
253  name.str("");
254  name << "sti" << i;
255  c2numpy_addcolumn(&npywriter, name.str().c_str(), C2NUMPY_UINT16);
256 
257  name.str("");
258  name << "stf" << i;
259  c2numpy_addcolumn(&npywriter, name.str().c_str(), C2NUMPY_UINT16);
260  }
261 
262  for (unsigned int i = 0;
263  i < (fUseFullWaveform ? detProp.ReadOutWindowSize() : fShortWaveformSize);
264  i++) {
265  std::ostringstream name;
266  name << "tck_" << i;
267  c2numpy_addcolumn(&npywriter, name.str().c_str(), C2NUMPY_INT16);
268  }
269 }
270 
271 //-----------------------------------------------------------------------
273 {
275 }
276 
277 //-----------------------------------------------------------------------
279 {
280  cout << "Event "
281  << " " << evt.id().run() << " " << evt.id().subRun() << " " << evt.id().event() << endl;
282 
283  std::unique_ptr<genFinder> gf(new genFinder());
284 
285  // ... Read in the digit List object(s).
287  std::vector<art::Ptr<raw::RawDigit>> rawdigitlist;
288  if (evt.getByLabel(fDigitModuleLabel, digitVecHandle)) {
289  art::fill_ptr_vector(rawdigitlist, digitVecHandle);
290  }
291 
292  // ... Read in the wire List object(s).
293  art::Handle<std::vector<recob::Wire>> wireListHandle;
294  std::vector<art::Ptr<recob::Wire>> wirelist;
295  if (evt.getByLabel(fWireProducerLabel, wireListHandle)) {
296  art::fill_ptr_vector(wirelist, wireListHandle);
297  }
298 
299  if (rawdigitlist.empty() && wirelist.empty()) return;
300  if (rawdigitlist.size() && wirelist.size()) return;
301 
302  // channel status
303  lariov::ChannelStatusProvider const& channelStatus =
305 
306  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
307  auto const detProp =
309 
310  // ... Use the handle to get a particular (0th) element of collection.
311  unsigned int dataSize;
312  if (rawdigitlist.size()) {
313  art::Ptr<raw::RawDigit> digitVec0(digitVecHandle, 0);
314  dataSize = digitVec0->Samples(); //size of raw data vectors
315  }
316  else {
317  dataSize = (wirelist[0]->Signal()).size();
318  }
319  if (dataSize != detProp.ReadOutWindowSize()) {
320  std::cout << "!!!!! Bad dataSize: " << dataSize << std::endl;
321  return;
322  }
323 
324  // ... Build a map from channel number -> rawdigitVec
325  std::map<raw::ChannelID_t, art::Ptr<raw::RawDigit>> rawdigitMap;
326  raw::ChannelID_t chnum = raw::InvalidChannelID; // channel number
327  if (rawdigitlist.size()) {
328  for (size_t rdIter = 0; rdIter < digitVecHandle->size(); ++rdIter) {
329  art::Ptr<raw::RawDigit> digitVec(digitVecHandle, rdIter);
330  chnum = digitVec->Channel();
331  if (chnum == raw::InvalidChannelID) continue;
332  rawdigitMap[chnum] = digitVec;
333  }
334  }
335  // ... Build a map from channel number -> wire
336  std::map<raw::ChannelID_t, art::Ptr<recob::Wire>> wireMap;
337  if (wirelist.size()) {
338  for (size_t ich = 0; ich < wirelist.size(); ++ich) {
339  art::Ptr<recob::Wire> wire = wirelist[ich];
340  chnum = wire->Channel();
341  if (chnum == raw::InvalidChannelID) continue;
342  wireMap[chnum] = wire;
343  }
344  }
345 
346  // ... Read in MC particle list
348  if (!evt.getByLabel(fSimulationProducerLabel, particleHandle)) {
349  throw cet::exception("AnalysisExample")
350  << " No simb::MCParticle objects in this event - "
351  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
352  }
353 
354  // ... Read in sim channel list
355  auto simChannelHandle = evt.getValidHandle<std::vector<sim::SimChannel>>(fSimChannelLabel);
356 
357  if (!simChannelHandle->size()) return;
358 
359  // ... Create a map of track IDs to generator labels
360  auto mcHandles = evt.getMany<std::vector<simb::MCTruth>>();
361  std::vector<std::pair<int, std::string>> track_id_to_label;
362 
363  for (auto const& mcHandle : mcHandles) {
364  const std::string& sModuleLabel = mcHandle.provenance()->moduleLabel();
366  std::vector<art::Ptr<simb::MCParticle>> mcParts = findMCParts.at(0);
367  for (const art::Ptr<simb::MCParticle> ptr : mcParts) {
368  int track_id = ptr->TrackId();
369  gf->add(track_id, sModuleLabel);
370  }
371  }
372 
373  std::string dummystr6 = "none ";
374  std::string dummystr7 = "none ";
375 
376  if (fSaveSignal) {
377  // .. create a channel number to trackid-wire signal info map
378  std::map<raw::ChannelID_t, std::map<int, WireSigInfo>> Ch2TrkWSInfoMap;
379 
380  // .. create a track ID to vector of channel numbers (in w/c this track deposited energy) map
381  std::map<int, std::vector<raw::ChannelID_t>> Trk2ChVecMap;
382 
383  // ... Loop over simChannels
384  for (auto const& channel : (*simChannelHandle)) {
385 
386  // .. get simChannel channel number
387  const raw::ChannelID_t ch1 = channel.Channel();
388  if (ch1 == raw::InvalidChannelID) continue;
389  if (geo::PlaneGeo::ViewName(fWireReadoutGeom->View(ch1)) != fPlaneToDump[0]) continue;
390 
391  bool selectThisChannel = false;
392 
393  // .. create a track ID to wire signal info map
394  std::map<int, WireSigInfo> Trk2WSInfoMap;
395 
396  // ... Loop over all ticks with ionization energy deposited
397  auto const& timeSlices = channel.TDCIDEMap();
398  for (auto const& timeSlice : timeSlices) {
399 
400  auto const& energyDeposits = timeSlice.second;
401  auto const tpctime = timeSlice.first;
402  unsigned int tdctick = static_cast<unsigned int>(clockData.TPCTDC2Tick(double(tpctime)));
403  if (tdctick < 0 || tdctick > (dataSize - 1)) continue;
404 
405  // ... Loop over all energy depositions in this tick
406  for (auto const& energyDeposit : energyDeposits) {
407 
408  if (!energyDeposit.trackID) continue;
409  int trkid = energyDeposit.trackID;
410  simb::MCParticle particle = PIS->TrackIdToMotherParticle(trkid);
411 
412  // .. ignore this energy deposition if incident particle energy below some threshold
413  if (particle.E() < fMinParticleEnergyGeV) continue;
414 
415  int eve_id = PIS->TrackIdToEveTrackId(trkid);
416  if (!eve_id) continue;
417  std::string genlab = gf->get_gen(eve_id);
418 
419  if (Trk2WSInfoMap.find(trkid) == Trk2WSInfoMap.end()) {
420  WireSigInfo wsinf;
421  wsinf.pdgcode = particle.PdgCode();
422  wsinf.genlab = genlab;
423  wsinf.procid = particle.Process();
424  wsinf.tdcmin = dataSize - 1;
425  wsinf.tdcmax = 0;
426  wsinf.edep = 0.;
427  wsinf.numel = 0;
428  Trk2WSInfoMap.insert(std::pair<int, WireSigInfo>(trkid, wsinf));
429  }
430  if (tdctick < Trk2WSInfoMap.at(trkid).tdcmin) Trk2WSInfoMap.at(trkid).tdcmin = tdctick;
431  if (tdctick > Trk2WSInfoMap.at(trkid).tdcmax) Trk2WSInfoMap.at(trkid).tdcmax = tdctick;
432  Trk2WSInfoMap.at(trkid).edep += energyDeposit.energy;
433  Trk2WSInfoMap.at(trkid).numel += energyDeposit.numElectrons;
434  }
435  }
436 
437  if (!Trk2WSInfoMap.empty()) {
438  for (std::pair<int, WireSigInfo> itmap : Trk2WSInfoMap) {
439  if (fSelectGenLabel != "ANY") {
440  if (itmap.second.genlab != fSelectGenLabel) continue;
441  }
442  if (fSelectProcID != "ANY") {
443  if (itmap.second.procid != fSelectProcID) continue;
444  }
445  if (fSelectPDGCode != 0) {
446  if (itmap.second.pdgcode != fSelectPDGCode) continue;
447  }
448  itmap.second.genlab.resize(6, ' ');
449  itmap.second.procid.resize(7, ' ');
450  if (itmap.second.numel >= fMinNumberOfElectrons &&
451  itmap.second.edep >= fMinEnergyDepositedMeV) {
452  if (fMaxNumberOfElectrons >= 0 && itmap.second.numel >= fMaxNumberOfElectrons) {
453  continue;
454  }
455  else {
456  int trkid = itmap.first;
457  if (Trk2ChVecMap.find(trkid) == Trk2ChVecMap.end()) {
458  std::vector<raw::ChannelID_t> chvec;
459  Trk2ChVecMap.insert(std::pair<int, std::vector<raw::ChannelID_t>>(trkid, chvec));
460  }
461  Trk2ChVecMap.at(trkid).push_back(ch1);
462  selectThisChannel = true;
463  }
464  }
465  } // loop over Trk2WSinfoMap
466  if (selectThisChannel) {
467  Ch2TrkWSInfoMap.insert(
468  std::pair<raw::ChannelID_t, std::map<int, WireSigInfo>>(ch1, Trk2WSInfoMap));
469  }
470  } // if Trk2WSInfoMap not empty
471 
472  } // loop over SimChannels
473 
474  std::set<raw::ChannelID_t> selected_channels;
475 
476  // ... Now write out the signal waveforms for each track
477  if (!Trk2ChVecMap.empty()) {
478  for (auto const& ittrk : Trk2ChVecMap) {
479  int i = fRandFlat.fireInt(
480  ittrk.second.size()); // randomly select one channel with a signal from this particle
481  chnum = ittrk.second[i];
482 
483  if (not selected_channels.insert(chnum).second) { continue; }
484 
485  std::map<raw::ChannelID_t, std::map<int, WireSigInfo>>::iterator itchn;
486  itchn = Ch2TrkWSInfoMap.find(chnum);
487  if (itchn != Ch2TrkWSInfoMap.end()) {
488 
489  std::vector<short> adcvec(dataSize); // vector to hold zero-padded full waveform
490 
491  if (rawdigitlist.size()) {
492  auto search = rawdigitMap.find(chnum);
493  if (search == rawdigitMap.end()) continue;
494  art::Ptr<raw::RawDigit> rawdig = (*search).second;
495  std::vector<short> rawadc(dataSize); // vector to hold uncompressed adc values later
496  raw::Uncompress(rawdig->ADCs(), rawadc, rawdig->GetPedestal(), rawdig->Compression());
497  for (size_t j = 0; j < rawadc.size(); ++j) {
498  adcvec[j] = rawadc[j] - rawdig->GetPedestal();
499  }
500  }
501  else if (wirelist.size()) {
502  auto search = wireMap.find(chnum);
503  if (search == wireMap.end()) continue;
504  art::Ptr<recob::Wire> wire = (*search).second;
505  const auto& signal = wire->Signal();
506  for (size_t j = 0; j < adcvec.size(); ++j) {
507  adcvec[j] = signal[j];
508  }
509  }
510 
511  // .. write out info for each peak
512  // a full waveform has at least one peak; the output will save up to 5 peaks (if there is
513  // only 1 peak, will fill the other 4 with 0);
514  // for fShortWaveformSize: only use the first peak's start_tick
515 
516  if (fUseFullWaveform) {
517 
518  c2numpy_uint32(&npywriter, evt.id().event());
519  c2numpy_uint32(&npywriter, chnum);
522  c2numpy_uint16(&npywriter, itchn->second.size()); // size of Trk2WSInfoMap, or #peaks
523  unsigned int icnt = 0;
524  for (auto const& it : itchn->second) {
525  c2numpy_int32(&npywriter, it.first); // trackid
526  c2numpy_int32(&npywriter, it.second.pdgcode); // pdgcode
527  c2numpy_string(&npywriter, it.second.genlab.c_str()); // genlab
528  c2numpy_string(&npywriter, it.second.procid.c_str()); // procid
529  c2numpy_float32(&npywriter, it.second.edep); // edepo
530  c2numpy_uint32(&npywriter, it.second.numel); // numelec
531 
532  c2numpy_uint16(&npywriter, it.second.tdcmin); // stck1
533  c2numpy_uint16(&npywriter, it.second.tdcmax); // stc2
534 
535  icnt++;
536  if (icnt == 5) break;
537  }
538 
539  // .. pad with 0's if number of peaks less than 5
540  for (unsigned int i = icnt; i < 5; ++i) {
543  c2numpy_string(&npywriter, dummystr6.c_str());
544  c2numpy_string(&npywriter, dummystr7.c_str());
549  }
550 
551  for (unsigned int itck = 0; itck < dataSize; ++itck) {
552  c2numpy_int16(&npywriter, adcvec[itck]);
553  }
554  }
555  else {
556 
557  // .. first loop to find largest signal
558  double EDep = 0.;
559  unsigned int TDCMin, TDCMax;
560  bool foundmaxsig = false;
561  for (auto& it : itchn->second) {
562  if (it.second.edep > EDep && it.second.numel > 0) {
563  EDep = it.second.edep;
564  TDCMin = it.second.tdcmin;
565  TDCMax = it.second.tdcmax;
566  foundmaxsig = true;
567  }
568  }
569  if (foundmaxsig) {
570  int sigtdc1, sigtdc2, sighwid, sigfwid, sigtdcm;
572  sigtdc1 = TDCMin - 14 / 2;
573  sigtdc2 = TDCMax + 3 * 14 / 2;
574  }
575  else {
576  sigtdc1 = TDCMin - 32 / 2;
577  sigtdc2 = TDCMax + 32 / 2;
578  }
579  sigfwid = sigtdc2 - sigtdc1;
580  sighwid = sigfwid / 2;
581  sigtdcm = sigtdc1 + sighwid;
582 
583  int start_tick = -1;
584  int end_tick = -1;
585  // .. set window edges to contain the largest signal
586  if (sigfwid < (int)fShortWaveformSize) {
587  // --> case 1: signal range fits within window
588  int dt = fShortWaveformSize - sigfwid;
589  start_tick = sigtdc1 - dt * fRandFlat.fire(0, 1);
590  }
591  else {
592  // --> case 2: signal range larger than window
593  int mrgn = fShortWaveformSize / 20;
594  int dt = fShortWaveformSize - 2 * mrgn;
595  start_tick = sigtdcm - mrgn - dt * fRandFlat.fire(0, 1);
596  }
597  if (start_tick < 0) start_tick = 0;
598  end_tick = start_tick + fShortWaveformSize - 1;
599  if (end_tick > int(dataSize - 1)) {
600  end_tick = dataSize - 1;
601  start_tick = end_tick - fShortWaveformSize + 1;
602  }
603 
604  c2numpy_uint32(&npywriter, evt.id().event());
605  c2numpy_uint32(&npywriter, chnum);
608 
609  // .. second loop to select only signals that are within the window
610 
611  int it_trk[5], it_pdg[5], it_nel[5];
612  unsigned int stck_1[5], stck_2[5];
613  std::string it_glb[5], it_prc[5];
614  double it_edp[5];
615 
616  unsigned int icnt = 0;
617 
618  for (auto& it : itchn->second) {
619  if ((it.second.tdcmin >= (unsigned int)start_tick &&
620  it.second.tdcmin < (unsigned int)end_tick) ||
621  (it.second.tdcmax > (unsigned int)start_tick &&
622  it.second.tdcmax <= (unsigned int)end_tick)) {
623 
624  it_trk[icnt] = it.first;
625  it_pdg[icnt] = it.second.pdgcode;
626  it_glb[icnt] = it.second.genlab;
627  it_prc[icnt] = it.second.procid;
628  it_edp[icnt] = it.second.edep;
629  it_nel[icnt] = it.second.numel;
630 
631  unsigned int mintdc = it.second.tdcmin;
632  unsigned int maxtdc = it.second.tdcmax;
633  if (mintdc < (unsigned int)start_tick) mintdc = start_tick;
634  if (maxtdc > (unsigned int)end_tick) maxtdc = end_tick;
635 
636  stck_1[icnt] = mintdc - start_tick;
637  stck_2[icnt] = maxtdc - start_tick;
638 
639  icnt++;
640  if (icnt == 5) break;
641  }
642  }
643 
644  c2numpy_uint16(&npywriter, icnt); // number of peaks
645 
646  for (unsigned int i = 0; i < icnt; ++i) {
647  c2numpy_int32(&npywriter, it_trk[i]); // trackid
648  c2numpy_int32(&npywriter, it_pdg[i]); // pdgcode
649  c2numpy_string(&npywriter, it_glb[i].c_str()); // genlab
650  c2numpy_string(&npywriter, it_prc[i].c_str()); // procid
651  c2numpy_float32(&npywriter, it_edp[i]); // edepo
652  c2numpy_uint32(&npywriter, it_nel[i]); // numelec
653  c2numpy_uint16(&npywriter, stck_1[i]); // stck1
654  c2numpy_uint16(&npywriter, stck_2[i]); // stck2
655  }
656 
657  // .. pad with 0's if number of peaks less than 5
658  for (unsigned int i = icnt; i < 5; ++i) {
661  c2numpy_string(&npywriter, dummystr6.c_str());
662  c2numpy_string(&npywriter, dummystr7.c_str());
667  }
668 
669  for (unsigned int itck = start_tick; itck < (start_tick + fShortWaveformSize);
670  ++itck) {
671  c2numpy_int16(&npywriter, adcvec[itck]);
672  }
673 
674  } // foundmaxsig
675  }
676  }
677  }
678  }
679  }
680  else {
681  //save noise
682  int noisechancount = 0;
683  std::map<raw::ChannelID_t, bool> signalMap;
684  for (auto const& channel : (*simChannelHandle)) {
685  signalMap[channel.Channel()] = true;
686  }
687  // .. create a vector for shuffling the wire channel indices
688  auto seed = std::chrono::high_resolution_clock::now().time_since_epoch().count();
689  size_t nchan = (rawdigitlist.empty() ? wirelist.size() : rawdigitlist.size());
690  std::vector<size_t> randigitmap;
691  for (size_t i = 0; i < nchan; ++i)
692  randigitmap.push_back(i);
693  std::shuffle(randigitmap.begin(), randigitmap.end(), std::mt19937(seed));
694 
695  for (size_t rdIter = 0; rdIter < (rawdigitlist.empty() ? wirelist.size() : rawdigitlist.size());
696  ++rdIter) {
697 
698  if (noisechancount == fMaxNoiseChannelsPerEvent) break;
699 
700  std::vector<short> adcvec(dataSize); // vector to wire adc values
701  if (rawdigitlist.size()) {
702  size_t ranIdx = randigitmap[rdIter];
703  art::Ptr<raw::RawDigit> digitVec(digitVecHandle, ranIdx);
704  if (signalMap[digitVec->Channel()]) continue;
705 
706  std::vector<short> rawadc(dataSize); // vector to hold uncompressed adc values later
708  continue;
709  raw::Uncompress(digitVec->ADCs(), rawadc, digitVec->GetPedestal(), digitVec->Compression());
710  for (size_t j = 0; j < rawadc.size(); ++j) {
711  adcvec[j] = rawadc[j] - digitVec->GetPedestal();
712  }
713  c2numpy_uint32(&npywriter, evt.id().event());
714  c2numpy_uint32(&npywriter, digitVec->Channel());
717  }
718  else if (wirelist.size()) {
719  size_t ranIdx = randigitmap[rdIter];
720  art::Ptr<recob::Wire> wire = wirelist[ranIdx];
721  if (signalMap[wire->Channel()]) continue;
722  if (channelStatus.IsBad(wire->Channel())) continue;
724  continue;
725  const auto& signal = wire->Signal();
726  for (size_t j = 0; j < adcvec.size(); ++j) {
727  adcvec[j] = signal[j];
728  }
729  c2numpy_uint32(&npywriter, evt.id().event());
730  c2numpy_uint32(&npywriter, wire->Channel());
733  }
734 
735  c2numpy_uint16(&npywriter, 0); //number of peaks
736  for (unsigned int i = 0; i < 5; ++i) {
739  c2numpy_string(&npywriter, dummystr6.c_str());
740  c2numpy_string(&npywriter, dummystr7.c_str());
745  }
746 
747  if (fUseFullWaveform) {
748  for (unsigned int itck = 0; itck < dataSize; ++itck) {
749  c2numpy_int16(&npywriter, adcvec[itck]);
750  }
751  }
752  else {
753  int start_tick = int((dataSize - fShortWaveformSize) * fRandFlat.fire(0, 1));
754  for (unsigned int itck = start_tick; itck < (start_tick + fShortWaveformSize); ++itck) {
755  c2numpy_int16(&npywriter, adcvec[itck]);
756  }
757  }
758 
759  ++noisechancount;
760  }
761  std::cout << "Total number of noise channels " << noisechancount << std::endl;
762  }
763 }
double E(const int i=0) const
Definition: MCParticle.h:234
float GetPedestal() const
Definition: RawDigit.h:221
int c2numpy_init(c2numpy_writer *writer, const std::string outputFilePrefix, int32_t numRowsPerFile)
Definition: c2numpy.h:140
std::string fDigitModuleLabel
module that made digits
base_engine_t & createEngine(seed_t seed)
bool has_gen(std::string gname)
int c2numpy_int32(c2numpy_writer *writer, int32_t data)
Definition: c2numpy.h:298
const ADCvector_t & ADCs() const
Reference to the compressed ADC count vector.
Definition: RawDigit.h:209
int PdgCode() const
Definition: MCParticle.h:213
ULong64_t Samples() const
Number of samples in the uncompressed ADC data.
Definition: RawDigit.h:217
c2numpy_type
Definition: c2numpy.h:29
static std::string ViewName(View_t view)
Returns the name of the specified view.
Definition: PlaneGeo.cxx:584
void add(const int &track_id, const std::string &gname)
RawWaveformDump(fhicl::ParameterSet const &p)
Declaration of signal hit object.
std::string fSimChannelLabel
module that made simchannels
ChannelID_t Channel() const
DAQ channel this raw data was read from.
Definition: RawDigit.h:213
Definition of basic raw digits.
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
std::string Process() const
Definition: MCParticle.h:216
Particle class.
void analyze(art::Event const &e) override
RunNumber_t run() const
Definition: EventID.h:98
int c2numpy_close(c2numpy_writer *writer)
Definition: c2numpy.h:425
std::pair< int, std::string > track_id_to_string
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
geo::WireReadoutGeom const * fWireReadoutGeom
constexpr ChannelID_t InvalidChannelID
ID of an invalid channel.
Definition: RawTypes.h:31
int c2numpy_addcolumn(c2numpy_writer *writer, const std::string name, c2numpy_type type)
Definition: c2numpy.h:157
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
void beginJob()
Definition: Breakpoints.cc:14
long seed
Definition: chem4.cc:67
raw::ChannelID_t Channel() const
Returns the ID of the channel (or InvalidChannelID)
Definition: Wire.h:223
Interface for a class providing readout channel mapping to geometry.
Collect all the RawData header files together.
Definition: EmTrack.h:40
int c2numpy_uint32(c2numpy_writer *writer, uint32_t data)
Definition: c2numpy.h:334
An art service to assist in the distribution of guaranteed unique seeds to all engines within an art ...
std::vector< float > Signal() const
Return a zero-padded full length vector filled with RoI signal.
Definition: Wire.cxx:30
int c2numpy_float32(c2numpy_writer *writer, float data)
Definition: c2numpy.h:369
raw::Compress_t Compression() const
Compression algorithm used to store the ADC counts.
Definition: RawDigit.h:229
simb::MCParticle TrackIdToMotherParticle(int const id) const
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Encapsulate the construction of a single detector plane .
View_t View(raw::ChannelID_t const channel) const
Returns the view (wire orientation) on the specified TPC channel.
int c2numpy_uint16(c2numpy_writer *writer, uint16_t data)
Definition: c2numpy.h:325
object containing MC truth information necessary for making RawDigits and doing back tracking ...
int c2numpy_int16(c2numpy_writer *writer, int16_t data)
Definition: c2numpy.h:289
EventNumber_t event() const
Definition: EventID.h:116
Declaration of basic channel signal object.
art::ServiceHandle< cheat::ParticleInventoryService > PIS
TCEvent evt
Definition: DataStructs.cxx:8
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
void Uncompress(const std::vector< short > &adc, std::vector< short > &uncompressed, raw::Compress_t compress)
Uncompresses a raw data buffer.
Definition: raw.cxx:744
Float_t e
Definition: plot.C:35
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:82
std::string fSimulationProducerLabel
producer that tracked simulated part. through detector
std::string get_gen(int tid)
SubRunNumber_t subRun() const
Definition: EventID.h:110
int c2numpy_string(c2numpy_writer *writer, const char *data)
Definition: c2numpy.h:411
EventID id() const
Definition: Event.cc:23
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::vector< Handle< PROD > > getMany(SelectorBase const &selector=MatchAllSelector{}) const