sim::MCRecoEdep Class Reference

#include <MCRecoEdep.h>

Public Member Functions
	MCRecoEdep (fhicl::ParameterSet const &pset)
	Default constructor with fhicl parameters. More...
void	MakeMCEdep (const std::vector< sim::SimChannel > &schArray)
void	MakeMCEdep (const std::vector< sim::SimEnergyDeposit > &sedArray)
void	MakeMCEdep (const std::vector< sim::SimEnergyDepositLite > &sedArray)
bool	ExistTrack (const unsigned int track_id) const
int	TrackToEdepIndex (unsigned int track_id) const
	Converts a track ID to MCEdep array index. Returns -1 if no corresponding array found . More...
const std::vector< sim::MCEdep > &	GetEdepArrayAt (size_t edep_index) const
	Returns a vector of MCEdep object at the given index. More...
const std::map< unsigned int, size_t >	TrackIndexMap () const
	Returns a map of track id <-> MCEdep vector index. More...
void	Clear ()

Protected Member Functions
std::vector< sim::MCEdep > &	__GetEdepArray__ (unsigned int track_id)

Protected Attributes
bool	_debug_mode
bool	_save_mchit
std::map< unsigned int, size_t >	_track_index
std::vector< std::vector < sim::MCEdep > >	_mc_edeps

Detailed Description

Definition at line 95 of file MCRecoEdep.h.

Constructor & Destructor Documentation

sim::MCRecoEdep::MCRecoEdep

(

fhicl::ParameterSet const &

pset

)

Default constructor with fhicl parameters.

Definition at line 42 of file MCRecoEdep.cxx.

  {
    _debug_mode = pset.get<bool>("DebugMode");
    _save_mchit = pset.get<bool>("SaveMCHit");
  }

Member Function Documentation

std::vector< sim::MCEdep > & sim::MCRecoEdep::__GetEdepArray__ ( unsigned int  track_id )

protected

Definition at line 56 of file MCRecoEdep.cxx.

  {
    if(ExistTrack(track_id)) return _mc_edeps.at((*_track_index.find(track_id)).second);
    _track_index.insert(std::pair<unsigned int,size_t>(track_id,_mc_edeps.size()));
    _mc_edeps.push_back(std::vector<sim::MCEdep>());
    return (*(_mc_edeps.rbegin()));
  }

void sim::MCRecoEdep::Clear ( )

inline

Definition at line 126 of file MCRecoEdep.h.

                 {
      _mc_edeps.clear();
      _track_index.clear();
      std::vector<std::vector<sim::MCEdep>>().swap(_mc_edeps);
      std::map<unsigned int,size_t>().swap(_track_index);
  }

bool sim::MCRecoEdep::ExistTrack ( const unsigned int  track_id ) const

inline

Definition at line 109 of file MCRecoEdep.h.

    { return (_track_index.find(track_id) != _track_index.end()); }

const std::vector< sim::MCEdep > & sim::MCRecoEdep::GetEdepArrayAt

(

size_t

edep_index

)

const

Returns a vector of MCEdep object at the given index.

Definition at line 49 of file MCRecoEdep.cxx.

  {
    if(edep_index >= _mc_edeps.size())
      throw cet::exception(__FUNCTION__) << Form("Track ID %zu not found!",edep_index);
    return _mc_edeps.at(edep_index);
  }

void sim::MCRecoEdep::MakeMCEdep

(

const std::vector< sim::SimChannel > &

schArray

)

Definition at line 64 of file MCRecoEdep.cxx.

  {
    _mc_edeps.clear();
    _track_index.clear();

    art::ServiceHandle<geo::Geometry const> geom;
//    const detinfo::DetectorProperties* detp = lar::providerFrom<detinfo::DetectorPropertiesService>();

    // Key map to identify a unique particle energy deposition point
    std::map<std::pair<UniquePosition, unsigned int>, int> hit_index_m;

    auto pindex = details::createPlaneIndexMap();

    if(_debug_mode) std::cout<<"Processing "<<schArray.size()<<" channels..."<<std::endl;
    // Loop over channels
    for(size_t i=0; i<schArray.size(); ++i) {

      // Get data to loop over
      auto const& sch = schArray[i];
      const auto &sch_map(sch.TDCIDEMap());
      // Channel
      UInt_t ch = sch.Channel();
      // Loop over ticks
      for(auto tdc_iter = sch_map.begin(); tdc_iter!=sch_map.end(); ++tdc_iter) {
        // for c2: hit_time is unused
        //unsigned short hit_time = (*tdc_iter).first;
        // Loop over IDEs
        for(auto const &ide : (*tdc_iter).second) {

          int track_id = ide.trackID;
          if(track_id < 0) track_id = track_id * (-1);
          unsigned int real_track_id = track_id;

          UniquePosition pos(ide.x, ide.y, ide.z);

          int hit_index = -1;
          auto key = std::make_pair(pos, real_track_id);
          auto hit_index_track_iter = hit_index_m.find(key);
          if(hit_index_track_iter == hit_index_m.end()) {
            int new_hit_index = this->__GetEdepArray__(real_track_id).size();
            hit_index_m[key]= new_hit_index;
          }
          else {
            hit_index = (*hit_index_track_iter).second;
          }
          auto const pid = geom->ChannelToWire(ch)[0].planeID();
          auto const channel_id = pindex[pid];
          double charge = ide.numElectrons;
          if(hit_index < 0) {
            // This particle energy deposition is never recorded so far. Create a new Edep
            //float charge = ide.numElectrons * detp->ElectronsToADC();
            this->__GetEdepArray__(real_track_id).emplace_back(pos, pid, pindex.size(), ide.energy, charge, channel_id);
          } else {
            // Append charge to the relevant edep (@ hit_index)
            //float charge = ide.numElectrons * detp->ElectronsToADC();
            MCEdep &edep = this->__GetEdepArray__(real_track_id).at(hit_index);
            edep.deps[channel_id].charge += charge;
            edep.deps[channel_id].energy += ide.energy;
          }
        } // end looping over ides in this tick
      } // end looping over ticks in this channel
    }// end looping over channels

    if(_debug_mode) {
      std::cout<< Form("  Collected %zu particles' energy depositions...",_mc_edeps.size()) << std::endl;
      // for c2: disable the entire loop instead of just the print statement
      //for(auto const& track_id_index : _track_index ) {
        //auto track_id   = track_id_index.first;
        //auto edep_index = track_id_index.second;
        //      std::cout<< Form("    Track ID: %d ... %zu Edep!", track_id, edep_index) << std::endl;
      //}
      std::cout<<std::endl;
    }
  }

void sim::MCRecoEdep::MakeMCEdep

(

const std::vector< sim::SimEnergyDeposit > &

sedArray

)

Definition at line 139 of file MCRecoEdep.cxx.

  {
    _mc_edeps.clear();
    _track_index.clear();

    art::ServiceHandle<geo::Geometry const> geom;
//    const detinfo::DetectorProperties* detp = lar::providerFrom<detinfo::DetectorPropertiesService>();

    // Key map to identify a unique particle energy deposition point
    std::map<std::pair<UniquePosition, unsigned int>, int> hit_index_m;

    auto pindex = details::createPlaneIndexMap();

    if(_debug_mode) std::cout<<"Processing "<<sedArray.size()<<" energy deposits..."<<std::endl;
    // Loop over channels
    for(size_t i=0; i<sedArray.size(); ++i) {

      // Get data to loop over
      auto const& sed = sedArray[i];

      // David Caratelli:
      // much of the code below is taken from the module:
      // https://cdcvs.fnal.gov/redmine/projects/larsim/repository/revisions/develop/entry/larsim/ElectronDrift/SimDriftElectrons_module.cc

      // given this SimEnergyDeposit, find the TPC channel information
      // "xyz" is the position of the energy deposit in world
      // coordinates. Note that the units of distance in
      // sim::SimEnergyDeposit are supposed to be cm.
      auto const mp = sed.MidPoint();
      double const xyz[3] = { mp.X(), mp.Y(), mp.Z() };
      // From the position in world coordinates, determine the
      // cryostat and tpc. If somehow the step is outside a tpc
      // (e.g., cosmic rays in rock) just move on to the next one.
      unsigned int cryostat = 0;
      try {
        geom->PositionToCryostat(xyz, cryostat);
      }
      catch(cet::exception &e){
        mf::LogWarning("SimDriftElectrons") << "step "// << energyDeposit << "\n"
                                            << "cannot be found in a cryostat\n"
                                            << e;
        continue;
      }
      unsigned int tpc = 0;
      try {
        geom->PositionToTPC(xyz, tpc, cryostat);
      }
      catch(cet::exception &e){
        mf::LogWarning("SimDriftElectrons") << "step "// << energyDeposit << "\n"
                                            << "cannot be found in a TPC\n"
                                            << e;
        continue;
      }
      const geo::TPCGeo& tpcGeo = geom->TPC(tpc, cryostat);

      //Define charge drift direction: driftcoordinate (x, y or z) and driftsign (positive or negative). Also define coordinates perpendicular to drift direction.
      // unused int driftcoordinate = std::abs(tpcGeo.DetectDriftDirection())-1;  //x:0, y:1, z:2

      // make a collection of electrons for each plane
      for(size_t p = 0; p < tpcGeo.Nplanes(); ++p){

        // grab the nearest channel to the fDriftClusterPos position
        // David Caratelli, comment begin:
        // NOTE: the below code works only when the drift coordinate is indeed in x (i.e. 0th coordinate)
        // see code linked above for a much more careful treatment of the coordinate system
        // David Caratelli, comment end.
        raw::ChannelID_t ch;
        try {
          ch = geom->NearestChannel(xyz, p, tpc, cryostat);
        }
        catch(cet::exception &e){
          mf::LogWarning("SimDriftElectrons") << "step "// << energyDeposit << "\n"
                                              << "nearest wire not in TPC\n"
                                              << e;
          continue;
        }

        int track_id = sed.TrackID();

        if(track_id < 0) track_id = track_id * (-1);
        unsigned int real_track_id = track_id;

        UniquePosition pos(xyz[0], xyz[1], xyz[2]);

        int hit_index = -1;
        auto key = std::make_pair(pos, real_track_id);
        auto hit_index_track_iter = hit_index_m.find(key);
        if(hit_index_track_iter == hit_index_m.end()) {
          int new_hit_index = this->__GetEdepArray__(real_track_id).size();
          hit_index_m[key]= new_hit_index;
        }
        else {
          hit_index = (*hit_index_track_iter).second;
        }
        auto const pid = geom->ChannelToWire(ch)[0].planeID();
        auto const channel_id = pindex[pid];
        double charge = sed.NumElectrons();
        if(hit_index < 0) {
          // This particle energy deposition is never recorded so far. Create a new Edep
          //float charge = ide.numElectrons * detp->ElectronsToADC();
          this->__GetEdepArray__(real_track_id).emplace_back(pos, pid, pindex.size(), sed.Energy(), charge, channel_id);
        } else {
          // Append charge to the relevant edep (@ hit_index)
          //float charge = ide.numElectrons * detp->ElectronsToADC();
          MCEdep &edep = this->__GetEdepArray__(real_track_id).at(hit_index);
          edep.deps[channel_id].charge += charge;
          edep.deps[channel_id].energy += sed.Energy();
        }
      } // end looping over planes
    }// end looping over SimEnergyDeposits

    if(_debug_mode) {
      std::cout<< Form("  Collected %zu particles' energy depositions...",_mc_edeps.size()) << std::endl;
      // for c2: disable the entire loop instead of just the print statement
      //for(auto const& track_id_index : _track_index ) {
      //auto track_id   = track_id_index.first;
      //auto edep_index = track_id_index.second;
      //        std::cout<< Form("    Track ID: %d ... %zu Edep!", track_id, edep_index) << std::endl;
      //}
      std::cout<<std::endl;
    }

    return;
  }

void sim::MCRecoEdep::MakeMCEdep

(

const std::vector< sim::SimEnergyDepositLite > &

sedArray

)

Definition at line 264 of file MCRecoEdep.cxx.

                                                                                {
    // Create a substitue array of sim::SimEnergyDeposit to avoid duplicating code...
    // Note that this makes use of the explicit conversion operator
    // defined in SimEnergyDepositLite class. Information will be partial.
    // Most notably for MakeMCEdep, charge (numElectrons) will be 0.
    std::vector<sim::SimEnergyDeposit> new_sedArray(sedArray.begin(), sedArray.end());
    MakeMCEdep(new_sedArray);
  }

const std::map<unsigned int,size_t> sim::MCRecoEdep::TrackIndexMap ( ) const

inline

Returns a map of track id <-> MCEdep vector index.

Definition at line 123 of file MCRecoEdep.h.

    { return _track_index; }

int sim::MCRecoEdep::TrackToEdepIndex ( unsigned int  track_id ) const

inline

Converts a track ID to MCEdep array index. Returns -1 if no corresponding array found .

Definition at line 113 of file MCRecoEdep.h.

    {
      auto iter = _track_index.find(track_id);
      return (iter == _track_index.end() ? -1 : (int)((*iter).second));
    }

Member Data Documentation

bool sim::MCRecoEdep::_debug_mode

protected

Definition at line 136 of file MCRecoEdep.h.

std::vector<std::vector<sim::MCEdep> > sim::MCRecoEdep::_mc_edeps

protected

Definition at line 139 of file MCRecoEdep.h.

bool sim::MCRecoEdep::_save_mchit

protected

Definition at line 137 of file MCRecoEdep.h.

std::map<unsigned int,size_t> sim::MCRecoEdep::_track_index

protected

Definition at line 138 of file MCRecoEdep.h.

---

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

• MCRecoEdep.h
• MCRecoEdep.cxx