321 std::unique_ptr<genFinder> gf(
new genFinder());
325 std::vector<art::Ptr<raw::RawDigit>> rawdigitlist;
333 std::vector<art::Ptr<raw::RawDigit>> rawdigitlist2;
339 if (rawdigitlist.empty() && rawdigitlist2.empty())
return;
346 unsigned int dataSize;
348 dataSize = digitVec0->Samples();
349 if (dataSize != detProp.ReadOutWindowSize()) {
350 throw cet::exception(
"RawWaveformClnSigDump") <<
"Bad dataSize: " << dataSize;
353 unsigned int dataSize2 = digitVec20->Samples();
354 if (dataSize != dataSize2) {
356 <<
"RawDigits from the 2 data products have different dataSizes: " << dataSize <<
"not eq to" 361 std::map<raw::ChannelID_t, art::Ptr<raw::RawDigit>> rawdigitMap;
363 if (rawdigitlist.size()) {
364 for (
size_t rdIter = 0; rdIter < digitVecHandle->size(); ++rdIter) {
366 chnum = digitVec->Channel();
368 rawdigitMap[chnum] = digitVec;
371 std::map<raw::ChannelID_t, art::Ptr<raw::RawDigit>> rawdigitMap2;
373 if (rawdigitlist2.size()) {
374 for (
size_t rdIter = 0; rdIter < digitVecHandle2->size(); ++rdIter) {
376 chnum2 = digitVec2->Channel();
378 rawdigitMap2[chnum2] = digitVec2;
386 <<
" No simb::MCParticle objects in this event - " 387 <<
" Line " << __LINE__ <<
" in file " << __FILE__ << std::endl;
391 auto simChannelHandle =
evt.getValidHandle<std::vector<sim::SimChannel>>(
fSimChannelLabel);
393 if (!simChannelHandle->size())
return;
399 auto mcHandles =
evt.getMany<std::vector<simb::MCTruth>>();
400 std::vector<std::pair<int, std::string>> track_id_to_label;
402 for (
auto const& mcHandle : mcHandles) {
403 const std::string& sModuleLabel = mcHandle.provenance()->moduleLabel();
405 std::vector<art::Ptr<simb::MCParticle>> mcParts = findMCParts.at(0);
407 int track_id = ptr->TrackId();
408 gf->add(track_id, sModuleLabel);
412 std::string dummystr6 =
"none ";
413 std::string dummystr7 =
"none ";
417 std::map<raw::ChannelID_t, std::map<int, WireSigInfo>> Ch2TrkWSInfoMap;
420 std::map<int, std::vector<raw::ChannelID_t>> Trk2ChVecMap;
423 for (
auto const& channel : (*simChannelHandle)) {
430 bool selectThisChannel =
false;
433 std::map<int, WireSigInfo> Trk2WSInfoMap;
436 auto const& timeSlices = channel.TDCIDEMap();
437 for (
auto const& timeSlice : timeSlices) {
439 auto const& energyDeposits = timeSlice.second;
440 auto const tpctime = timeSlice.first;
441 unsigned int tdctick =
static_cast<unsigned int>(clockData.TPCTDC2Tick(
double(tpctime)));
442 if (tdctick < 0 || tdctick > (dataSize - 1))
continue;
445 for (
auto const& energyDeposit : energyDeposits) {
447 if (!energyDeposit.trackID)
continue;
448 int trkid = energyDeposit.trackID;
456 if (!eve_id)
continue;
457 std::string genlab = gf->get_gen(eve_id);
459 if (Trk2WSInfoMap.find(trkid) == Trk2WSInfoMap.end()) {
464 wsinf.
tdcmin = dataSize - 1;
470 Trk2WSInfoMap.insert(std::pair<int, WireSigInfo>(trkid, wsinf));
472 if (tdctick < Trk2WSInfoMap.at(trkid).tdcmin) Trk2WSInfoMap.at(trkid).tdcmin = tdctick;
473 if (tdctick > Trk2WSInfoMap.at(trkid).tdcmax) Trk2WSInfoMap.at(trkid).tdcmax = tdctick;
474 Trk2WSInfoMap.at(trkid).edep += energyDeposit.energy;
475 Trk2WSInfoMap.at(trkid).numel += energyDeposit.numElectrons;
479 auto search2 = rawdigitMap2.find(ch1);
480 if (search2 == rawdigitMap2.end())
continue;
482 std::vector<short> rawadc(dataSize);
484 std::vector<short> adcvec2(dataSize);
485 for (
size_t j = 0; j < rawadc.size(); ++j) {
489 if (!Trk2WSInfoMap.empty()) {
490 for (std::pair<int, WireSigInfo> itmap : Trk2WSInfoMap) {
494 for (
size_t i = itmap.second.tdcmin; i <= itmap.second.tdcmax; i++) {
495 if (
abs(adcvec2[i]) >
abs(pkadc)) {
500 Trk2WSInfoMap.at(itmap.first).tdcpeak = pktdc;
501 Trk2WSInfoMap.at(itmap.first).adcpeak = pkadc;
512 itmap.second.genlab.resize(6,
' ');
513 itmap.second.procid.resize(7,
' ');
520 int trkid = itmap.first;
521 if (Trk2ChVecMap.find(trkid) == Trk2ChVecMap.end()) {
522 std::vector<raw::ChannelID_t> chvec;
523 Trk2ChVecMap.insert(std::pair<
int, std::vector<raw::ChannelID_t>>(trkid, chvec));
525 Trk2ChVecMap.at(trkid).push_back(ch1);
526 selectThisChannel =
true;
530 if (selectThisChannel) {
531 Ch2TrkWSInfoMap.insert(
532 std::pair<
raw::ChannelID_t, std::map<int, WireSigInfo>>(ch1, Trk2WSInfoMap));
538 std::set<raw::ChannelID_t> selected_channels;
541 if (!Trk2ChVecMap.empty()) {
544 using trk2chvecpair = std::pair<const int, std::vector<raw::ChannelID_t>>;
545 std::vector<const trk2chvecpair*> tchpv;
546 for (
auto const& ittrk : Trk2ChVecMap) {
547 tchpv.emplace_back(&ittrk);
549 auto seed = std::chrono::high_resolution_clock::now().time_since_epoch().count();
552 int signalchancount = 0;
553 for (
auto const& itpr : tchpv) {
556 itpr->second.size());
557 chnum = itpr->second[i];
559 if (not selected_channels.insert(chnum).second) {
continue; }
561 std::map<raw::ChannelID_t, std::map<int, WireSigInfo>>
::iterator itchn;
562 itchn = Ch2TrkWSInfoMap.find(chnum);
563 if (itchn != Ch2TrkWSInfoMap.end()) {
565 std::vector<short> adcvec(dataSize);
567 auto search = rawdigitMap.find(chnum);
568 if (search == rawdigitMap.end())
continue;
570 std::vector<short> rawadc(dataSize);
572 for (
size_t j = 0; j < rawadc.size(); ++j) {
576 std::vector<short> adcvec2(
579 auto search2 = rawdigitMap2.find(chnum);
580 if (search2 == rawdigitMap2.end())
continue;
583 for (
size_t j = 0; j < rawadc.size(); ++j) {
598 unsigned int icnt = 0;
599 for (
auto& it : itchn->second) {
613 if (icnt == 5)
break;
617 for (
unsigned int i = icnt; i < 5; ++i) {
630 for (
unsigned int itck = 0; itck < dataSize; ++itck) {
633 for (
unsigned int itck = 0; itck < dataSize; ++itck) {
642 unsigned int TDCMin, TDCMax;
643 bool foundmaxsig =
false;
644 for (
auto& it : itchn->second) {
645 if (it.second.edep > EDep && it.second.adcpeak != 0 && it.second.numel > 0) {
646 EDep = it.second.edep;
647 TDCMin = it.second.tdcmin;
648 TDCMax = it.second.tdcmax;
653 int sigtdc1, sigtdc2, sighwid, sigfwid, sigtdcm;
662 sigfwid = sigtdc2 - sigtdc1;
663 sighwid = sigfwid / 2;
664 sigtdcm = sigtdc1 + sighwid;
671 int dt = fShortWaveformSize - sigfwid;
672 start_tick = sigtdc1 - dt *
fRandFlat.fire(0, 1);
676 int mrgn = fShortWaveformSize / 20;
677 int dt = fShortWaveformSize - 2 * mrgn;
678 start_tick = sigtdcm - mrgn - dt *
fRandFlat.fire(0, 1);
680 if (start_tick < 0) start_tick = 0;
681 end_tick = start_tick + fShortWaveformSize - 1;
682 if (end_tick >
int(dataSize - 1)) {
683 end_tick = dataSize - 1;
684 start_tick = end_tick - fShortWaveformSize + 1;
693 int it_trk[5], it_pdg[5], it_nel[5], pk_tdc[5], pk_adc[5];
694 unsigned int stck_1[5], stck_2[5];
695 std::string it_glb[5], it_prc[5];
698 unsigned int icnt = 0;
700 for (
auto& it : itchn->second) {
702 if ((it.second.tdcmin >= (
unsigned int)start_tick &&
703 it.second.tdcmin < (
unsigned int)end_tick) ||
704 (it.second.tdcmax > (
unsigned int)start_tick &&
705 it.second.tdcmax <= (
unsigned int)end_tick)) {
707 it_trk[icnt] = it.first;
708 it_pdg[icnt] = it.second.pdgcode;
709 it_glb[icnt] = it.second.genlab;
710 it_prc[icnt] = it.second.procid;
711 it_edp[icnt] = it.second.edep;
712 it_nel[icnt] = it.second.numel;
714 unsigned int mintdc = it.second.tdcmin;
715 unsigned int maxtdc = it.second.tdcmax;
716 if (mintdc < (
unsigned int)start_tick) mintdc = start_tick;
717 if (maxtdc > (
unsigned int)end_tick) maxtdc = end_tick;
719 stck_1[icnt] = mintdc - start_tick;
720 stck_2[icnt] = maxtdc - start_tick;
721 pk_tdc[icnt] = it.second.tdcpeak - start_tick;
722 pk_adc[icnt] = it.second.adcpeak;
725 if (icnt == 5)
break;
731 for (
unsigned int i = 0; i < icnt; ++i) {
745 for (
unsigned int i = icnt; i < 5; ++i) {
775 int noisechancount = 0;
776 std::map<raw::ChannelID_t, bool> signalMap;
777 for (
auto const& channel : (*simChannelHandle)) {
778 signalMap[channel.Channel()] =
true;
781 auto seed = std::chrono::high_resolution_clock::now().time_since_epoch().count();
782 std::vector<size_t> randigitmap;
783 for (
size_t i = 0; i < rawdigitlist.size(); ++i)
784 randigitmap.push_back(i);
785 std::shuffle(randigitmap.begin(), randigitmap.end(), std::mt19937(
seed));
787 for (
size_t rdIter = 0; rdIter < rawdigitlist.size(); ++rdIter) {
791 std::vector<float> adcvec(dataSize);
793 size_t ranIdx = randigitmap[rdIter];
795 if (signalMap[digitVec->Channel()])
continue;
797 std::vector<short> rawadc(dataSize);
799 raw::Uncompress(digitVec->ADCs(), rawadc, digitVec->GetPedestal(), digitVec->Compression());
800 for (
size_t j = 0; j < rawadc.size(); ++j) {
801 adcvec[j] = rawadc[j] - digitVec->GetPedestal();
808 for (
unsigned int i = 0; i < 5; ++i) {
822 for (
unsigned int itck = 0; itck < dataSize; ++itck) {
827 int start_tick = int((dataSize - fShortWaveformSize) *
fRandFlat.fire(0, 1));
828 for (
unsigned int itck = start_tick; itck < (start_tick +
fShortWaveformSize); ++itck) {
835 std::cout <<
"Total number of noise channels " << noisechancount << std::endl;
double E(const int i=0) const
float GetPedestal() const
int c2numpy_int32(c2numpy_writer *writer, int32_t data)
const ADCvector_t & ADCs() const
Reference to the compressed ADC count vector.
static std::string ViewName(geo::View_t view)
Returns the name of the specified view.
constexpr auto abs(T v)
Returns the absolute value of the argument.
std::string Process() const
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
int TrackIdToEveTrackId(int tid) const
constexpr ChannelID_t InvalidChannelID
ID of an invalid channel.
int c2numpy_uint32(c2numpy_writer *writer, uint32_t data)
int c2numpy_float32(c2numpy_writer *writer, float data)
raw::Compress_t Compression() const
Compression algorithm used to store the ADC counts.
simb::MCParticle TrackIdToMotherParticle(int const id) const
View_t View(PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
int c2numpy_uint16(c2numpy_writer *writer, uint16_t data)
int c2numpy_int16(c2numpy_writer *writer, int16_t data)
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
unsigned int ChannelID_t
Type representing the ID of a readout channel.
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
void Uncompress(const std::vector< short > &adc, std::vector< short > &uncompressed, raw::Compress_t compress)
Uncompresses a raw data buffer.
int c2numpy_string(c2numpy_writer *writer, const char *data)
cet::coded_exception< error, detail::translate > exception