nnet::RawWaveformClnSigDump Class Reference

Inheritance diagram for nnet::RawWaveformClnSigDump:

Public Member Functions
	RawWaveformClnSigDump (fhicl::ParameterSet const &p)
	RawWaveformClnSigDump (RawWaveformClnSigDump const &)=delete
	RawWaveformClnSigDump (RawWaveformClnSigDump &&)=delete
RawWaveformClnSigDump &	operator= (RawWaveformClnSigDump const &)=delete
RawWaveformClnSigDump &	operator= (RawWaveformClnSigDump &&)=delete
void	analyze (art::Event const &e) override
void	beginJob () override
void	endJob () override
Public Member Functions inherited from art::EDAnalyzer
	EDAnalyzer (fhicl::ParameterSet const &pset)
template<typename Config >
	EDAnalyzer (Table< Config > const &config)
std::string	workerType () const
Public Member Functions inherited from art::detail::Analyzer
virtual	~Analyzer () noexcept
	Analyzer (fhicl::ParameterSet const &pset)
template<typename Config >
	Analyzer (Table< Config > const &config)
void	doBeginJob (SharedResources const &resources)
void	doEndJob ()
void	doRespondToOpenInputFile (FileBlock const &fb)
void	doRespondToCloseInputFile (FileBlock const &fb)
void	doRespondToOpenOutputFiles (FileBlock const &fb)
void	doRespondToCloseOutputFiles (FileBlock const &fb)
bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)
bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
Public Member Functions inherited from art::Observer
	~Observer () noexcept
	Observer (Observer const &)=delete
	Observer (Observer &&)=delete
Observer &	operator= (Observer const &)=delete
Observer &	operator= (Observer &&)=delete
void	registerProducts (ProductDescriptions &, ModuleDescription const &)
void	fillDescriptions (ModuleDescription const &)
fhicl::ParameterSetID	selectorConfig () const
Public Member Functions inherited from art::ModuleBase
virtual	~ModuleBase () noexcept
	ModuleBase ()
ModuleDescription const &	moduleDescription () const
void	setModuleDescription (ModuleDescription const &)
std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const
void	sortConsumables (std::string const &current_process_name)
template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)
template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

Private Attributes
std::string	fDumpWaveformsFileName
std::string	fDumpCleanSignalFileName
std::string	fSimulationProducerLabel
	producer that tracked simulated part. through detector More...
std::string	fSimChannelLabel
	module that made simchannels More...
std::string	fDigitModuleLabel
	module that made digits More...
std::string	fCleanSignalDigitModuleLabel
	module that made the signal-only digits More...
bool	fUseFullWaveform
unsigned int	fShortWaveformSize
int	fEstIndFWForOffset
int	fEstColFWForOffset
std::string	fSelectGenLabel
std::string	fSelectProcID
int	fSelectPDGCode
std::string	fPlaneToDump
double	fMinParticleEnergyGeV
double	fMinEnergyDepositedMeV
int	fMinNumberOfElectrons
int	fMaxNumberOfElectrons
int	fMinPureSignalADCs
bool	fSaveSignal
int	fMaxNoiseChannelsPerEvent
std::string	fCollectionPlaneLabel
art::ServiceHandle< geo::Geometry >	fgeom
art::ServiceHandle< cheat::ParticleInventoryService >	PIS
CLHEP::RandFlat	fRandFlat
c2numpy_writer	npywriter
c2numpy_writer	npywriter2

Additional Inherited Members
Public Types inherited from art::EDAnalyzer
using	WorkerType = WorkerT< EDAnalyzer >
using	ModuleType = EDAnalyzer
Protected Member Functions inherited from art::Observer
std::string const &	processName () const
bool	wantAllEvents () const noexcept
bool	wantEvent (ScheduleID id, Event const &e) const
Handle< TriggerResults >	getTriggerResults (Event const &e) const
	Observer (fhicl::ParameterSet const &config)
	Observer (std::vector< std::string > const &select_paths, std::vector< std::string > const &reject_paths, fhicl::ParameterSet const &config)
Protected Member Functions inherited from art::ModuleBase
ConsumesCollector &	consumesCollector ()
template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)
template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)
template<typename T , BranchType = InEvent>
void	consumesMany ()
template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)
template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)
template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

Detailed Description

Definition at line 85 of file RawWaveformClnSigDump_module.cc.

Constructor & Destructor Documentation

nnet::RawWaveformClnSigDump::RawWaveformClnSigDump ( fhicl::ParameterSet const &  p )

explicit

Definition at line 185 of file RawWaveformClnSigDump_module.cc.

  : EDAnalyzer{p}
  , fDumpWaveformsFileName(p.get<std::string>("DumpWaveformsFileName", "dumpwaveforms"))
  , fDumpCleanSignalFileName(p.get<std::string>("CleanSignalFileName", "dumpcleansignal"))
  , fSimulationProducerLabel(p.get<std::string>("SimulationProducerLabel", "larg4Main"))
  , fSimChannelLabel(p.get<std::string>("SimChannelLabel", "elecDrift"))
  , fDigitModuleLabel(p.get<std::string>("DigitModuleLabel", "simWire"))
  , fCleanSignalDigitModuleLabel(p.get<std::string>("CleanSignalDigitModuleLabel", "simWire:signal"))
  , fUseFullWaveform(p.get<bool>("UseFullWaveform", true))
  , fShortWaveformSize(p.get<unsigned int>("ShortWaveformSize"))
  , fEstIndFWForOffset(p.get<int>("EstIndFWForOffset"))
  , fEstColFWForOffset(p.get<int>("EstIndFWForOffset"))
  , fSelectGenLabel(p.get<std::string>("SelectGenLabel", "ANY"))
  , fSelectProcID(p.get<std::string>("SelectProcID", "ANY"))
  , fSelectPDGCode(p.get<int>("SelectPDGCode", 0))
  , fPlaneToDump(p.get<std::string>("PlaneToDump"))
  , fMinParticleEnergyGeV(p.get<double>("MinParticleEnergyGeV", 0.))
  , fMinEnergyDepositedMeV(p.get<double>("MinEnergyDepositedMeV", 0.))
  , fMinNumberOfElectrons(p.get<int>("MinNumberOfElectrons", 600))
  , fMaxNumberOfElectrons(p.get<int>("MaxNumberOfElectrons", -1))
  , fMinPureSignalADCs(p.get<int>("MinPureSignalADCs", 0))
  , fSaveSignal(p.get<bool>("SaveSignal", true))
  , fMaxNoiseChannelsPerEvent(p.get<int>("MaxNoiseChannelsPerEvent"))
  , fCollectionPlaneLabel(p.get<std::string>("CollectionPlaneLabel"))
  , fRandFlat{createEngine(art::ServiceHandle<rndm::NuRandomService>{}->declareEngine(
                           instanceName, p, "SeedForRawWaveformDump"),
                           "HepJamesRandom", instanceName)}
{
  if (std::getenv("CLUSTER") && std::getenv("PROCESS")) {
    fDumpWaveformsFileName += string(std::getenv("CLUSTER")) + "-" + string(std::getenv("PROCESS")) + "-";
    fDumpCleanSignalFileName += string(std::getenv("CLUSTER")) + "-" + string(std::getenv("PROCESS")) + "-";
  }

  if (fDigitModuleLabel.empty() && fCleanSignalDigitModuleLabel.empty()) {
    throw cet::exception("RawWaveformClnSigDump")
      << "Both DigitModuleLabel and CleanSignalModuleLabel are empty";
  }
}

nnet::RawWaveformClnSigDump::RawWaveformClnSigDump ( RawWaveformClnSigDump const &  )

delete

nnet::RawWaveformClnSigDump::RawWaveformClnSigDump ( RawWaveformClnSigDump &&  )

delete

Member Function Documentation

void nnet::RawWaveformClnSigDump::analyze ( art::Event const &  e )

overridevirtual

Implements art::EDAnalyzer.

Definition at line 310 of file RawWaveformClnSigDump_module.cc.

{
  cout << "Event "
       << " " << evt.id().run() << " " << evt.id().subRun() << " " << evt.id().event() << endl;

  std::unique_ptr<genFinder> gf(new genFinder());

  // ... Read in the digit List object(s).
  art::Handle<std::vector<raw::RawDigit>> digitVecHandle;
  std::vector<art::Ptr<raw::RawDigit>> rawdigitlist;
  if (evt.getByLabel(fDigitModuleLabel, digitVecHandle)) {
    //std::cout << " !!!! RawWaveformClnSigDump: fDigitModuleLabel -> " << fDigitModuleLabel << std::endl; 
    art::fill_ptr_vector(rawdigitlist, digitVecHandle);
  }

  // ... Read in the signal-only digit List object(s).
  art::Handle<std::vector<raw::RawDigit>> digitVecHandle2;
  std::vector<art::Ptr<raw::RawDigit>> rawdigitlist2;
  if (evt.getByLabel(fCleanSignalDigitModuleLabel, digitVecHandle2)) {
    //std::cout << " !!!! RawWaveformClnSigDump: fCleanSignalDigitModuleLabel -> " << fCleanSignalDigitModuleLabel << std::endl; 
    art::fill_ptr_vector(rawdigitlist2, digitVecHandle2);
  }

  if (rawdigitlist.empty() && rawdigitlist2.empty()) return;

  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
  auto const detProp =
    art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);

  // ... Use the handle to get a particular (0th) element of collection.
  unsigned int dataSize;
  art::Ptr<raw::RawDigit> digitVec0(digitVecHandle, 0);
  dataSize = digitVec0->Samples(); //size of raw data vectors
  if (dataSize != detProp.ReadOutWindowSize()) {
    throw cet::exception("RawWaveformClnSigDump") << "Bad dataSize: " << dataSize;
  }
  art::Ptr<raw::RawDigit> digitVec20(digitVecHandle2, 0);
  unsigned int dataSize2 = digitVec20->Samples();
  if (dataSize != dataSize2) {
    throw cet::exception("RawWaveformClnSigDump")
      << "RawDigits from the 2 data products have different dataSizes: " << dataSize << "not eq to" << dataSize2;
  }

  // ... Build a map from channel number -> rawdigitVec
  std::map<raw::ChannelID_t, art::Ptr<raw::RawDigit>> rawdigitMap;
  raw::ChannelID_t chnum = raw::InvalidChannelID; // channel number
  if (rawdigitlist.size()) {
    for (size_t rdIter = 0; rdIter < digitVecHandle->size(); ++rdIter) {
      art::Ptr<raw::RawDigit> digitVec(digitVecHandle, rdIter);
      chnum = digitVec->Channel();
      if (chnum == raw::InvalidChannelID) continue;
      rawdigitMap[chnum] = digitVec;
    }
  }
  std::map<raw::ChannelID_t, art::Ptr<raw::RawDigit>> rawdigitMap2;
  raw::ChannelID_t chnum2 = raw::InvalidChannelID; // channel number
  if (rawdigitlist2.size()) {
    for (size_t rdIter = 0; rdIter < digitVecHandle2->size(); ++rdIter) {
      art::Ptr<raw::RawDigit> digitVec2(digitVecHandle2, rdIter);
      chnum2 = digitVec2->Channel();
      if (chnum2 == raw::InvalidChannelID) continue;
      rawdigitMap2[chnum2] = digitVec2;
    }
  }

  // ... Read in MC particle list
  art::Handle<std::vector<simb::MCParticle>> particleHandle;
  if (!evt.getByLabel(fSimulationProducerLabel, particleHandle)) {
    throw cet::exception("AnalysisExample")
      << " No simb::MCParticle objects in this event - "
      << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
  }

  // ... Read in sim channel list
  auto simChannelHandle =
    evt.getValidHandle<std::vector<sim::SimChannel>>(fSimChannelLabel);

  if (!simChannelHandle->size()) return;

  // ... Create a map of track IDs to generator labels
  //Get a list of generator names.
  //std::vector<art::Handle<std::vector<simb::MCTruth>>> mcHandles;
  //evt.getManyByType(mcHandles);
  auto mcHandles = evt.getMany<std::vector<simb::MCTruth>>();
  std::vector<std::pair<int, std::string>> track_id_to_label;

  for (auto const& mcHandle : mcHandles) {
    const std::string& sModuleLabel = mcHandle.provenance()->moduleLabel();
    art::FindManyP<simb::MCParticle> findMCParts(mcHandle, evt, fSimulationProducerLabel);
    std::vector<art::Ptr<simb::MCParticle>> mcParts = findMCParts.at(0);
    for (const art::Ptr<simb::MCParticle> ptr : mcParts) {
      int track_id = ptr->TrackId();
      gf->add(track_id, sModuleLabel);
    }
  }

  std::string dummystr6 = "none  ";
  std::string dummystr7 = "none   ";

  if (fSaveSignal) {
    // .. create a channel number to trackid-wire signal info map
    std::map<raw::ChannelID_t, std::map<int, WireSigInfo>> Ch2TrkWSInfoMap;

    // .. create a track ID to vector of channel numbers (in w/c this track deposited energy) map
    std::map<int, std::vector<raw::ChannelID_t>> Trk2ChVecMap;

    // ... Loop over simChannels
    for (auto const& channel : (*simChannelHandle)) {

      // .. get simChannel channel number
      const raw::ChannelID_t ch1 = channel.Channel();
      if (ch1 == raw::InvalidChannelID) continue;
      if (geo::PlaneGeo::ViewName(fgeom->View(ch1)) != fPlaneToDump[0]) continue;

      bool selectThisChannel = false;

      // .. create a track ID to wire signal info map
      std::map<int, WireSigInfo> Trk2WSInfoMap;

      // ... Loop over all ticks with ionization energy deposited
      auto const& timeSlices = channel.TDCIDEMap();
      for (auto const& timeSlice : timeSlices) {

        auto const& energyDeposits = timeSlice.second;
        auto const tpctime = timeSlice.first;
        unsigned int tdctick = static_cast<unsigned int>(clockData.TPCTDC2Tick(double(tpctime)));
        if (tdctick < 0 || tdctick > (dataSize - 1)) continue;

        // ... Loop over all energy depositions in this tick
        for (auto const& energyDeposit : energyDeposits) {

          if (!energyDeposit.trackID) continue;
          int trkid = energyDeposit.trackID;
          simb::MCParticle particle = PIS->TrackIdToMotherParticle(trkid);
          //std::cout << energyDeposit.trackID << " " << trkid << " " << particle.TrackId() << std::endl;

          // .. ignore this energy deposition if incident particle energy below some threshold
          if (particle.E() < fMinParticleEnergyGeV) continue;

          int eve_id = PIS->TrackIdToEveTrackId(trkid);
          if (!eve_id) continue;
          std::string genlab = gf->get_gen(eve_id);

          if (Trk2WSInfoMap.find(trkid) == Trk2WSInfoMap.end()) {
            WireSigInfo wsinf;
            wsinf.pdgcode = particle.PdgCode();
            wsinf.genlab = genlab;
            wsinf.procid = particle.Process();
            wsinf.tdcmin = dataSize - 1;
            wsinf.tdcmax = 0;
            wsinf.tdcpeak = -1;
            wsinf.adcpeak = 0;
            wsinf.edep = 0.;
            wsinf.numel = 0;
            Trk2WSInfoMap.insert(std::pair<int, WireSigInfo>(trkid, wsinf));
          }
          if (tdctick < Trk2WSInfoMap.at(trkid).tdcmin) Trk2WSInfoMap.at(trkid).tdcmin = tdctick;
          if (tdctick > Trk2WSInfoMap.at(trkid).tdcmax) Trk2WSInfoMap.at(trkid).tdcmax = tdctick;
          Trk2WSInfoMap.at(trkid).edep += energyDeposit.energy;
          Trk2WSInfoMap.at(trkid).numel += energyDeposit.numElectrons;
        }
      } // loop over timeSlices

      auto search2 = rawdigitMap2.find(ch1);
      if (search2 == rawdigitMap2.end()) continue;
      art::Ptr<raw::RawDigit> rawdig2 = (*search2).second;
      std::vector<short> rawadc(dataSize);
      raw::Uncompress(rawdig2->ADCs(), rawadc, rawdig2->GetPedestal(), rawdig2->Compression());
      std::vector<short> adcvec2(dataSize);
      for (size_t j = 0; j < rawadc.size(); ++j) {
        adcvec2[j] = rawadc[j] - rawdig2->GetPedestal();
      }

      if (!Trk2WSInfoMap.empty()) {
        for (std::pair<int, WireSigInfo> itmap : Trk2WSInfoMap) {
          // find the peak adc value in the signal-only raw digits within the range tdcmin->tdcmax
          int pkadc = 0;
          int pktdc = -1;
          for (size_t i = itmap.second.tdcmin; i <= itmap.second.tdcmax; i++) {
            if (abs(adcvec2[i]) > abs(pkadc)) {
              pkadc = adcvec2[i];
              pktdc = i;
            }
          }
          Trk2WSInfoMap.at(itmap.first).tdcpeak = pktdc;
          Trk2WSInfoMap.at(itmap.first).adcpeak = pkadc;

          if (fSelectGenLabel != "ANY") {
            if (itmap.second.genlab != fSelectGenLabel) continue;
          }
          if (fSelectProcID != "ANY") {
            if (itmap.second.procid != fSelectProcID) continue;
          }
          if (fSelectPDGCode != 0) {
            if (itmap.second.pdgcode != fSelectPDGCode) continue;
          }
          itmap.second.genlab.resize(6, ' ');
          itmap.second.procid.resize(7, ' ');
          if (itmap.second.numel >= fMinNumberOfElectrons &&
              itmap.second.edep >= fMinEnergyDepositedMeV && abs(pkadc) >= fMinPureSignalADCs) {
            if (fMaxNumberOfElectrons >= 0 && itmap.second.numel >= fMaxNumberOfElectrons) {
              continue;
            }
            else {
              int trkid = itmap.first;
              if (Trk2ChVecMap.find(trkid) == Trk2ChVecMap.end()) {
                std::vector<raw::ChannelID_t> chvec;
                Trk2ChVecMap.insert(std::pair<int, std::vector<raw::ChannelID_t>>(trkid, chvec));
              }
              Trk2ChVecMap.at(trkid).push_back(ch1);
              selectThisChannel = true;
            }
          }
        } // loop over Trk2WSinfoMap
        if (selectThisChannel) {
          Ch2TrkWSInfoMap.insert(
            std::pair<raw::ChannelID_t, std::map<int, WireSigInfo>>(ch1, Trk2WSInfoMap));
        }
      } // if Trk2WSInfoMap not empty

    } // loop over SimChannels

    std::set<raw::ChannelID_t> selected_channels;

    // ... Now write out the signal waveforms for each track
    if (!Trk2ChVecMap.empty()) {
      for (auto const& ittrk : Trk2ChVecMap) {
        int i = fRandFlat.fireInt(ittrk.second.size()); // randomly select one channel with a signal from this particle
        chnum = ittrk.second[i];

        if (not selected_channels.insert(chnum).second) {
          continue;
        }

        std::map<raw::ChannelID_t, std::map<int, WireSigInfo>>::iterator itchn;
        itchn = Ch2TrkWSInfoMap.find(chnum);
        if (itchn != Ch2TrkWSInfoMap.end()) {

          std::vector<short> adcvec(dataSize); // vector to hold zero-padded full waveform

          auto search = rawdigitMap.find(chnum);
          if (search == rawdigitMap.end()) continue;
          art::Ptr<raw::RawDigit> rawdig = (*search).second;
          std::vector<short> rawadc(dataSize); // vector to hold uncompressed adc values later
          raw::Uncompress(rawdig->ADCs(), rawadc, rawdig->GetPedestal(), rawdig->Compression());
          for (size_t j = 0; j < rawadc.size(); ++j) {
            adcvec[j] = rawadc[j] - rawdig->GetPedestal();
          }

          std::vector<short> adcvec2(dataSize); // vector to hold zero-padded full signal-only waveform

          auto search2 = rawdigitMap2.find(chnum);
          if (search2 == rawdigitMap2.end()) continue;
          art::Ptr<raw::RawDigit> rawdig2 = (*search2).second;
          raw::Uncompress(rawdig2->ADCs(), rawadc, rawdig2->GetPedestal(), rawdig2->Compression());
          for (size_t j = 0; j < rawadc.size(); ++j) {
            adcvec2[j] = rawadc[j] - rawdig2->GetPedestal();
          }

          // .. write out info for each peak
          //    a full waveform has at least one peak; the output will save up to 5 peaks (if there is
          //    only 1 peak, will fill the other 4 with 0);
          //    for fShortWaveformSize: only use the first peak's start_tick

          if (fUseFullWaveform) {

            c2numpy_uint32(&npywriter, evt.id().event());
            c2numpy_uint32(&npywriter, chnum);
            c2numpy_string(&npywriter, geo::PlaneGeo::ViewName(fgeom->View(chnum)).c_str());
            c2numpy_uint16(&npywriter, itchn->second.size()); // size of Trk2WSInfoMap, or #peaks
            unsigned int icnt = 0;
            for (auto & it : itchn->second) {
              c2numpy_int32(&npywriter, it.first);                  // trackid
              c2numpy_int32(&npywriter, it.second.pdgcode);         // pdgcode
              c2numpy_string(&npywriter, it.second.genlab.c_str()); // genlab
              c2numpy_string(&npywriter, it.second.procid.c_str()); // procid
              c2numpy_float32(&npywriter, it.second.edep);          // edepo
              c2numpy_uint32(&npywriter, it.second.numel);          // numelec

              c2numpy_uint16(&npywriter, it.second.tdcmin);         // stck1
              c2numpy_uint16(&npywriter, it.second.tdcmax);         // stck2
              c2numpy_int32(&npywriter, it.second.tdcpeak);         // pktdc
              c2numpy_int32(&npywriter, it.second.adcpeak);         // pkadc

              icnt++;
              if (icnt == 5) break;
            }

            // .. pad with 0's if number of peaks less than 5
            for (unsigned int i = icnt; i < 5; ++i) {
              c2numpy_int32(&npywriter, 0);
              c2numpy_int32(&npywriter, 0);
              c2numpy_string(&npywriter, dummystr6.c_str());
              c2numpy_string(&npywriter, dummystr7.c_str());
              c2numpy_float32(&npywriter, 0.);
              c2numpy_uint32(&npywriter, 0);
              c2numpy_uint16(&npywriter, 0);
              c2numpy_uint16(&npywriter, 0);
              c2numpy_int32(&npywriter, 0);
              c2numpy_int32(&npywriter, 0);
            }

            for (unsigned int itck = 0; itck < dataSize; ++itck) {
              c2numpy_int16(&npywriter, adcvec[itck]);
            }
            for (unsigned int itck = 0; itck < dataSize; ++itck) {
              c2numpy_int16(&npywriter2, adcvec2[itck]);
            }

          } else {

            // .. first loop to find largest signal
            double EDep = 0.;
            unsigned int TDCMin, TDCMax;
            bool foundmaxsig = false;
            for (auto & it : itchn->second) {
              if (it.second.edep > EDep && it.second.adcpeak != 0 && it.second.numel > 0){ 
                EDep = it.second.edep;
                TDCMin = it.second.tdcmin;
                TDCMax = it.second.tdcmax;
                foundmaxsig = true;
              }
            }
            if (foundmaxsig) {
              int sigtdc1, sigtdc2, sighwid, sigfwid, sigtdcm;
              if (fPlaneToDump!=fCollectionPlaneLabel){
                sigtdc1 = TDCMin - fEstIndFWForOffset/2;
                sigtdc2 = TDCMax + 3*fEstIndFWForOffset/2;
              } else {
                sigtdc1 = TDCMin - fEstColFWForOffset/2;
                sigtdc2 = TDCMax + fEstColFWForOffset/2;
              }
              sigfwid=sigtdc2 - sigtdc1;
              sighwid=sigfwid/2;
              sigtdcm=sigtdc1+sighwid;

              int start_tick = -1;
              int end_tick = -1;
              // .. set window edges to contain the largest signal
              if (sigfwid < (int)fShortWaveformSize) {
                // --> case 1: signal range fits within window
                int dt = fShortWaveformSize - sigfwid;
                start_tick = sigtdc1 - dt * fRandFlat.fire(0,1);
              }
              else {
                // --> case 2: signal range larger than window
                int mrgn = fShortWaveformSize/20;
                int dt = fShortWaveformSize - 2*mrgn;
                start_tick = sigtdcm - mrgn - dt * fRandFlat.fire(0,1);
              }
              if (start_tick < 0) start_tick = 0;
              end_tick = start_tick + fShortWaveformSize - 1;
              if (end_tick > int (dataSize - 1)) {
                end_tick = dataSize - 1;
                start_tick = end_tick - fShortWaveformSize + 1;
              }

              c2numpy_uint32(&npywriter, evt.id().event());
              c2numpy_uint32(&npywriter, chnum);
              c2numpy_string(&npywriter, geo::PlaneGeo::ViewName(fgeom->View(chnum)).c_str());

              // .. second loop to select only signals that are within the window

              int it_trk[5],it_pdg[5],it_nel[5],pk_tdc[5],pk_adc[5];
              unsigned int stck_1[5],stck_2[5];
              std::string it_glb[5], it_prc[5];
              double it_edp[5];

              unsigned int icnt = 0;

              for (auto & it : itchn->second) {
                if (abs(it.second.adcpeak) < fMinPureSignalADCs) continue;
                if (( it.second.tdcmin >= (unsigned int)start_tick && it.second.tdcmin <  (unsigned int)end_tick) ||
                    ( it.second.tdcmax >  (unsigned int)start_tick && it.second.tdcmax <= (unsigned int)end_tick)) {

                  it_trk[icnt] = it.first;       
                  it_pdg[icnt] = it.second.pdgcode;       
                  it_glb[icnt] = it.second.genlab;          
                  it_prc[icnt] = it.second.procid;          
                  it_edp[icnt] = it.second.edep;    
                  it_nel[icnt] = it.second.numel;   

                  unsigned int mintdc = it.second.tdcmin;
                  unsigned int maxtdc = it.second.tdcmax;
                  if (mintdc < (unsigned int)start_tick)mintdc = start_tick; 
                  if (maxtdc > (unsigned int)end_tick)maxtdc = end_tick; 

                  stck_1[icnt] = mintdc - start_tick;       
                  stck_2[icnt] = maxtdc - start_tick;       
                  pk_tdc[icnt] = it.second.tdcpeak - start_tick;       
                  pk_adc[icnt] = it.second.adcpeak;       

                  icnt++;
                  if (icnt == 5) break;
                }
              }

              c2numpy_uint16(&npywriter, icnt); // number of peaks

              for (unsigned int i = 0; i < icnt; ++i) {
                  c2numpy_int32(&npywriter,  it_trk[i]);         // trackid
                  c2numpy_int32(&npywriter,  it_pdg[i]);         // pdgcode
                  c2numpy_string(&npywriter, it_glb[i].c_str()); // genlab
                  c2numpy_string(&npywriter, it_prc[i].c_str()); // procid
                  c2numpy_float32(&npywriter,it_edp[i]);         // edepo
                  c2numpy_uint32(&npywriter, it_nel[i]);         // numelec
                  c2numpy_uint16(&npywriter, stck_1[i]); // stck1
                  c2numpy_uint16(&npywriter, stck_2[i]); // stck2
                  c2numpy_int32(&npywriter,  pk_tdc[i]); // pktdc
                  c2numpy_int32(&npywriter,  pk_adc[i]); // pkadc
              }

              // .. pad with 0's if number of peaks less than 5
              for (unsigned int i = icnt; i < 5; ++i) {
                c2numpy_int32(&npywriter, 0);
                c2numpy_int32(&npywriter, 0);
                c2numpy_string(&npywriter, dummystr6.c_str());
                c2numpy_string(&npywriter, dummystr7.c_str());
                c2numpy_float32(&npywriter, 0.);
                c2numpy_uint32(&npywriter, 0);
                c2numpy_uint16(&npywriter, 0);
                c2numpy_uint16(&npywriter, 0);
                c2numpy_int32(&npywriter, 0);
                c2numpy_int32(&npywriter, 0);
              }

              for (unsigned int itck = start_tick; itck < (start_tick + fShortWaveformSize); ++itck) {
                c2numpy_int16(&npywriter, adcvec[itck]);
              }
              for (unsigned int itck = start_tick; itck < (start_tick + fShortWaveformSize); ++itck) {
                c2numpy_int16(&npywriter2, adcvec2[itck]);
              }

            } // foundmaxsig
          }
        }
      }
    }
  }
  else {
    //save noise
    int noisechancount = 0;
    std::map<raw::ChannelID_t, bool> signalMap;
    for (auto const& channel : (*simChannelHandle)) {
      signalMap[channel.Channel()] = true;
    }
    // .. create a vector for shuffling the wire channel indices
    auto seed = std::chrono::high_resolution_clock::now().time_since_epoch().count();
    std::vector<size_t> randigitmap;
    for (size_t i=0; i<rawdigitlist.size(); ++i) randigitmap.push_back(i);
    std::shuffle ( randigitmap.begin(), randigitmap.end(), std::mt19937(seed) );

    for (size_t rdIter = 0; rdIter < rawdigitlist.size(); ++rdIter) {

      if (noisechancount==fMaxNoiseChannelsPerEvent)break;

      std::vector<float> adcvec(dataSize); // vector to wire adc values

      size_t ranIdx=randigitmap[rdIter];
      art::Ptr<raw::RawDigit> digitVec(digitVecHandle, ranIdx);
      if (signalMap[digitVec->Channel()]) continue;

      std::vector<short> rawadc(dataSize); // vector to hold uncompressed adc values later
      if (geo::PlaneGeo::ViewName(fgeom->View(digitVec->Channel())) != fPlaneToDump[0]) continue;
      raw::Uncompress(digitVec->ADCs(), rawadc, digitVec->GetPedestal(), digitVec->Compression());
      for (size_t j = 0; j < rawadc.size(); ++j) {
        adcvec[j] = rawadc[j] - digitVec->GetPedestal();
      }
      c2numpy_uint32(&npywriter, evt.id().event());
      c2numpy_uint32(&npywriter, digitVec->Channel());
      c2numpy_string(&npywriter,
                     geo::PlaneGeo::ViewName(fgeom->View(digitVec->Channel())).c_str());

      c2numpy_uint16(&npywriter, 0); //number of peaks
      for (unsigned int i = 0; i < 5; ++i) {
        c2numpy_int32(&npywriter, 0);
        c2numpy_int32(&npywriter, 0);
        c2numpy_string(&npywriter, dummystr6.c_str());
        c2numpy_string(&npywriter, dummystr7.c_str());
        c2numpy_float32(&npywriter, 0.);
        c2numpy_uint32(&npywriter, 0);
        c2numpy_uint16(&npywriter, 0);
        c2numpy_uint16(&npywriter, 0);
        c2numpy_int32(&npywriter, 0);
        c2numpy_int32(&npywriter, 0);
      }

      if (fUseFullWaveform) {
        for (unsigned int itck = 0; itck < dataSize; ++itck) {
          c2numpy_int16(&npywriter, short(adcvec[itck]));
        }
      }
      else {
        int start_tick = int((dataSize - fShortWaveformSize) * fRandFlat.fire(0, 1));
        for (unsigned int itck = start_tick; itck < (start_tick + fShortWaveformSize); ++itck) {
          c2numpy_int16(&npywriter, short(adcvec[itck]));
        }
      }

      ++noisechancount;
    }
    std::cout << "Total number of noise channels " << noisechancount << std::endl;
  }
}

void nnet::RawWaveformClnSigDump::beginJob ( )

overridevirtual

Reimplemented from art::EDAnalyzer.

Definition at line 226 of file RawWaveformClnSigDump_module.cc.

{
  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob();

  c2numpy_init(&npywriter, fDumpWaveformsFileName, 50000);
  c2numpy_addcolumn(&npywriter, "evt", C2NUMPY_UINT32);
  c2numpy_addcolumn(&npywriter, "chan", C2NUMPY_UINT32);
  c2numpy_addcolumn(&npywriter, "view", (c2numpy_type)((int)C2NUMPY_STRING + 1));
  c2numpy_addcolumn(&npywriter, "ntrk", C2NUMPY_UINT16);

  for (unsigned int i = 0; i < 5; i++) {
    std::ostringstream name;

    name.str("");
    name << "tid" << i;
    c2numpy_addcolumn(&npywriter, name.str().c_str(), C2NUMPY_INT32);

    name.str("");
    name << "pdg" << i;
    c2numpy_addcolumn(&npywriter, name.str().c_str(), C2NUMPY_INT32);

    name.str("");
    name << "gen" << i;
    c2numpy_addcolumn(&npywriter, name.str().c_str(), (c2numpy_type)((int)C2NUMPY_STRING + 6));

    name.str("");
    name << "pid" << i;
    c2numpy_addcolumn(&npywriter, name.str().c_str(), (c2numpy_type)((int)C2NUMPY_STRING + 7));

    name.str("");
    name << "edp" << i;
    c2numpy_addcolumn(&npywriter, name.str().c_str(), C2NUMPY_FLOAT32);

    name.str("");
    name << "nel" << i;
    c2numpy_addcolumn(&npywriter, name.str().c_str(), C2NUMPY_UINT32);

    name.str("");
    name << "sti" << i;
    c2numpy_addcolumn(&npywriter, name.str().c_str(), C2NUMPY_UINT16);

    name.str("");
    name << "stf" << i;
    c2numpy_addcolumn(&npywriter, name.str().c_str(), C2NUMPY_UINT16);

    name.str("");
    name << "stp" << i;
    c2numpy_addcolumn(&npywriter, name.str().c_str(), C2NUMPY_INT32);

    name.str("");
    name << "adc" << i;
    c2numpy_addcolumn(&npywriter, name.str().c_str(), C2NUMPY_INT32);
  }

  for (unsigned int i = 0;
       i < (fUseFullWaveform ? detProp.ReadOutWindowSize() : fShortWaveformSize);
       i++) {
    std::ostringstream name;
    name << "tck_" << i;
    c2numpy_addcolumn(&npywriter, name.str().c_str(), C2NUMPY_INT16);
  }

  // ... this is for storing the clean signal (no noise) waveform
  c2numpy_init(&npywriter2, fDumpCleanSignalFileName, 50000);

  for (unsigned int i = 0;
       i < (fUseFullWaveform ? detProp.ReadOutWindowSize() : fShortWaveformSize);
       i++) {
    std::ostringstream name;
    name << "tck_" << i;
    c2numpy_addcolumn(&npywriter2, name.str().c_str(), C2NUMPY_INT16);
  }
}

void nnet::RawWaveformClnSigDump::endJob ( )

overridevirtual

Reimplemented from art::EDAnalyzer.

Definition at line 302 of file RawWaveformClnSigDump_module.cc.

{
  c2numpy_close(&npywriter);
  c2numpy_close(&npywriter2);
}

RawWaveformClnSigDump& nnet::RawWaveformClnSigDump::operator= ( RawWaveformClnSigDump const &  )

delete

RawWaveformClnSigDump& nnet::RawWaveformClnSigDump::operator= ( RawWaveformClnSigDump &&  )

delete

Member Data Documentation

std::string nnet::RawWaveformClnSigDump::fCleanSignalDigitModuleLabel

private

module that made the signal-only digits

Definition at line 111 of file RawWaveformClnSigDump_module.cc.

std::string nnet::RawWaveformClnSigDump::fCollectionPlaneLabel

private

Definition at line 128 of file RawWaveformClnSigDump_module.cc.

std::string nnet::RawWaveformClnSigDump::fDigitModuleLabel

private

module that made digits

Definition at line 110 of file RawWaveformClnSigDump_module.cc.

std::string nnet::RawWaveformClnSigDump::fDumpCleanSignalFileName

private

Definition at line 106 of file RawWaveformClnSigDump_module.cc.

std::string nnet::RawWaveformClnSigDump::fDumpWaveformsFileName

private

Definition at line 105 of file RawWaveformClnSigDump_module.cc.

int nnet::RawWaveformClnSigDump::fEstColFWForOffset

private

Definition at line 115 of file RawWaveformClnSigDump_module.cc.

int nnet::RawWaveformClnSigDump::fEstIndFWForOffset

private

Definition at line 114 of file RawWaveformClnSigDump_module.cc.

art::ServiceHandle<geo::Geometry> nnet::RawWaveformClnSigDump::fgeom

private

Definition at line 129 of file RawWaveformClnSigDump_module.cc.

int nnet::RawWaveformClnSigDump::fMaxNoiseChannelsPerEvent

private

Definition at line 127 of file RawWaveformClnSigDump_module.cc.

int nnet::RawWaveformClnSigDump::fMaxNumberOfElectrons

private

Definition at line 124 of file RawWaveformClnSigDump_module.cc.

double nnet::RawWaveformClnSigDump::fMinEnergyDepositedMeV

private

Definition at line 122 of file RawWaveformClnSigDump_module.cc.

int nnet::RawWaveformClnSigDump::fMinNumberOfElectrons

private

Definition at line 123 of file RawWaveformClnSigDump_module.cc.

double nnet::RawWaveformClnSigDump::fMinParticleEnergyGeV

private

Definition at line 121 of file RawWaveformClnSigDump_module.cc.

int nnet::RawWaveformClnSigDump::fMinPureSignalADCs

private

Definition at line 125 of file RawWaveformClnSigDump_module.cc.

std::string nnet::RawWaveformClnSigDump::fPlaneToDump

private

Definition at line 120 of file RawWaveformClnSigDump_module.cc.

CLHEP::RandFlat nnet::RawWaveformClnSigDump::fRandFlat

private

Definition at line 132 of file RawWaveformClnSigDump_module.cc.

bool nnet::RawWaveformClnSigDump::fSaveSignal

private

Definition at line 126 of file RawWaveformClnSigDump_module.cc.

std::string nnet::RawWaveformClnSigDump::fSelectGenLabel

private

Definition at line 117 of file RawWaveformClnSigDump_module.cc.

int nnet::RawWaveformClnSigDump::fSelectPDGCode

private

Definition at line 119 of file RawWaveformClnSigDump_module.cc.

std::string nnet::RawWaveformClnSigDump::fSelectProcID

private

Definition at line 118 of file RawWaveformClnSigDump_module.cc.

unsigned int nnet::RawWaveformClnSigDump::fShortWaveformSize

private

Definition at line 113 of file RawWaveformClnSigDump_module.cc.

std::string nnet::RawWaveformClnSigDump::fSimChannelLabel

private

module that made simchannels

Definition at line 109 of file RawWaveformClnSigDump_module.cc.

std::string nnet::RawWaveformClnSigDump::fSimulationProducerLabel

private

producer that tracked simulated part. through detector

Definition at line 108 of file RawWaveformClnSigDump_module.cc.

bool nnet::RawWaveformClnSigDump::fUseFullWaveform

private

Definition at line 112 of file RawWaveformClnSigDump_module.cc.

c2numpy_writer nnet::RawWaveformClnSigDump::npywriter

private

Definition at line 134 of file RawWaveformClnSigDump_module.cc.

c2numpy_writer nnet::RawWaveformClnSigDump::npywriter2

private

Definition at line 135 of file RawWaveformClnSigDump_module.cc.

art::ServiceHandle<cheat::ParticleInventoryService> nnet::RawWaveformClnSigDump::PIS

private

Definition at line 130 of file RawWaveformClnSigDump_module.cc.

---

The documentation for this class was generated from the following file:

• larrecodnn/larrecodnn/ImagePatternAlgs/Modules/RawWaveformClnSigDump_module.cc