FTLDigitizer.h

Go to the documentation of this file.

#ifndef FastTimingSimProducers_FastTimingCommon_FTLDigitizer_h
#define FastTimingSimProducers_FastTimingCommon_FTLDigitizer_h

#include "SimFastTiming/FastTimingCommon/interface/FTLDigitizerBase.h"

#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"

#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/ESWatcher.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "Geometry/Records/interface/IdealGeometryRecord.h"
#include "Geometry/HGCalCommonData/interface/FastTimeDDDConstants.h"

#include "SimGeneral/MixingModule/interface/PileUpEventPrincipal.h"

#include <vector>
#include <unordered_map>
#include <unordered_set>
#include <memory>
#include <tuple>



namespace std {
  template<>
  struct hash<DetId> {
    std::size_t operator()(const DetId& detid) const {
      return hash<uint32_t>()(detid.rawId());
    }
  };
}

namespace ftl_digitizer {
  
  namespace FTLHelpers {
    // index , det id, time
    typedef std::tuple<int,uint32_t,float> FTLCaloHitTuple_t;
    
    bool orderByDetIdThenTime(const FTLCaloHitTuple_t &a, const FTLCaloHitTuple_t &b)
    {
      unsigned int detId_a(std::get<1>(a)), detId_b(std::get<1>(b));
      
      if(detId_a<detId_b) return true;
      if(detId_a>detId_b) return false;
      
      double time_a(std::get<2>(a)), time_b(std::get<2>(b));
      if(time_a<time_b) return true;
      
      return false;
    }
  }

  template<class SensorPhysics, class ElectronicsSim>
  class FTLDigitizer : public FTLDigitizerBase
  {
  public:
    
  FTLDigitizer(const edm::ParameterSet& config, 
           edm::ConsumesCollector& iC,
           edm::stream::EDProducerBase& parent) :
    FTLDigitizerBase(config,iC,parent),
    deviceSim_( config.getParameterSet("DeviceSimulation") ),
    electronicsSim_( config.getParameterSet("ElectronicsSimulation") ),        
    maxSimHitsAccTime_( config.getParameter< uint32_t >("maxSimHitsAccTime") ),
    bxTime_( config.getParameter< double >("bxTime") ),         
    tofDelay_( config.getParameter< double >("tofDelay") ) { }
    
    ~FTLDigitizer() override { }
    
    void accumulate(edm::Event const& e, edm::EventSetup const& c, CLHEP::HepRandomEngine* hre) override;
    void accumulate(PileUpEventPrincipal const& e, edm::EventSetup const& c, CLHEP::HepRandomEngine* hre) override;
    void accumulate(edm::Handle<edm::PSimHitContainer> const &hits, int bxCrossing, CLHEP::HepRandomEngine* hre) override;
    
    void initializeEvent(edm::Event const& e, edm::EventSetup const& c) override;
    void finalizeEvent(edm::Event& e, edm::EventSetup const& c, CLHEP::HepRandomEngine* hre) override;
    
    void beginRun(const edm::EventSetup & es) override; 
    void endRun() override {}
    
  private :
    
    void resetSimHitDataAccumulator() {
      FTLSimHitDataAccumulator().swap(simHitAccumulator_);
    }
    
    // implementations
    SensorPhysics deviceSim_; // processes a given simhit into an entry in a FTLSimHitDataAccumulator
    ElectronicsSim electronicsSim_; // processes a FTLSimHitDataAccumulator into a FTLDigiCollection
        
    //handle sim hits
    const int maxSimHitsAccTime_;
    const double bxTime_;
    FTLSimHitDataAccumulator simHitAccumulator_;  
        
    //delay to apply after evaluating time of arrival at the sensitive detector
    const float tofDelay_;

    //geometries
    std::unordered_set<DetId> validIds_;
    edm::ESWatcher<IdealGeometryRecord> idealGeomWatcher_;
    edm::ESHandle<FastTimeDDDConstants> dddFTL_;   
  };

  template<class SensorPhysics, class ElectronicsSim>
  void FTLDigitizer<SensorPhysics,ElectronicsSim>::accumulate(edm::Event const& e, 
                                  edm::EventSetup const& c, 
                                  CLHEP::HepRandomEngine* hre) {
    edm::Handle<edm::PSimHitContainer> simHits;
    e.getByLabel(inputSimHits_, simHits);
    accumulate(simHits,0,hre);
  }

  template<class SensorPhysics, class ElectronicsSim>
  void FTLDigitizer<SensorPhysics,ElectronicsSim>::accumulate(PileUpEventPrincipal const& e, 
                                  edm::EventSetup const& c, 
                                  CLHEP::HepRandomEngine* hre){
    edm::Handle<edm::PSimHitContainer> simHits;
    e.getByLabel(inputSimHits_, simHits);
    accumulate(simHits,e.bunchCrossing(),hre);
  }

  template<class SensorPhysics, class ElectronicsSim>
  void FTLDigitizer<SensorPhysics,ElectronicsSim>::accumulate(edm::Handle<edm::PSimHitContainer> const &hits, 
                                  int bxCrossing, 
                                  CLHEP::HepRandomEngine* hre) {
    using namespace FTLHelpers;
    //configuration to apply for the computation of time-of-flight
    bool weightToAbyEnergy(false);
    float tdcOnset(0.f);
    
    //create list of tuples (pos in container, RECO DetId, time) to be sorted first
    int nchits=(int)hits->size();  
    std::vector< FTLCaloHitTuple_t > hitRefs;
    hitRefs.reserve(nchits);
    for(int i=0; i<nchits; ++i) {
      const auto& the_hit = hits->at(i);    
      
      DetId id = ( validIds_.count(the_hit.detUnitId()) ? the_hit.detUnitId() : 0 );
      
      if (verbosity_>0) {   
    edm::LogInfo("HGCDigitizer") << " i/p " << std::hex << the_hit.detUnitId() << std::dec 
                     << " o/p " << id.rawId() << std::endl;
      }
      
      if( 0 != id.rawId() ) {      
    hitRefs.emplace_back( i, id.rawId(), the_hit.tof() );
      }
    }
    std::sort(hitRefs.begin(),hitRefs.end(),FTLHelpers::orderByDetIdThenTime);
    
    //loop over sorted hits
    nchits = hitRefs.size();
    for(int i=0; i<nchits; ++i) {
      const int hitidx   = std::get<0>(hitRefs[i]);
      const uint32_t id  = std::get<1>(hitRefs[i]);
      
      //get the data for this cell, if not available then we skip it
      
      if( !validIds_.count(id) ) continue;
      auto simHitIt = simHitAccumulator_.emplace(id,FTLCellInfo()).first;
      
      if(id==0) continue; // to be ignored at RECO level
      
      const float toa    = std::get<2>(hitRefs[i]);
      const PSimHit &hit=hits->at( hitidx );     
      const float charge = deviceSim_.getChargeForHit(hit);
      
      //distance to the center of the detector
      const float dist2center( 0.1f*hit.entryPoint().mag() );
      
      //hit time: [time()]=ns  [centerDist]=cm [refSpeed_]=cm/ns + delay by 1ns
      //accumulate in 15 buckets of 25ns (9 pre-samples, 1 in-time, 5 post-samples)
      const float tof = toa-dist2center/refSpeed_+tofDelay_ ;
      const int itime= std::floor( tof/bxTime_ ) + 9;
      
      //no need to add bx crossing - tof comes already corrected from the mixing module
      //itime += bxCrossing;
      //itime += 9;
      
      if(itime<0 || itime>14) continue;     
      
      //check if time index is ok and store energy
      if(itime >= (int)simHitIt->second.hit_info[0].size() ) continue;
      
      (simHitIt->second).hit_info[0][itime] += charge;
      float accCharge=(simHitIt->second).hit_info[0][itime];
      
      //time-of-arrival (check how to be used)
      if(weightToAbyEnergy) (simHitIt->second).hit_info[1][itime] += charge*tof;
      else if((simHitIt->second).hit_info[1][itime]==0)
    {   
      if( accCharge>tdcOnset )
        {
          //extrapolate linear using previous simhit if it concerns to the same DetId
          float fireTDC=tof;
          if(i>0)
        {
          uint32_t prev_id  = std::get<1>(hitRefs[i-1]);
          if(prev_id==id)
            {
              float prev_toa    = std::get<2>(hitRefs[i-1]);
              float prev_tof(prev_toa-dist2center/refSpeed_+tofDelay_);
              //float prev_charge = std::get<3>(hitRefs[i-1]);
              float deltaQ2TDCOnset = tdcOnset-((simHitIt->second).hit_info[0][itime]-charge);
              float deltaQ          = charge;
              float deltaT          = (tof-prev_tof);
              fireTDC               = deltaT*(deltaQ2TDCOnset/deltaQ)+prev_tof;
            }         
        }
          
          (simHitIt->second).hit_info[1][itime]=fireTDC;
        }
    }
    }
    hitRefs.clear();
  }
  
  template<class SensorPhysics, class ElectronicsSim>
  void FTLDigitizer<SensorPhysics,ElectronicsSim>::initializeEvent(edm::Event const& e, edm::EventSetup const& c) {
    deviceSim_.getEvent(e);
    electronicsSim_.getEvent(e);
  }
  
  template<class SensorPhysics, class ElectronicsSim>
  void FTLDigitizer<SensorPhysics,ElectronicsSim>::finalizeEvent(edm::Event& e, edm::EventSetup const& c, 
                                 CLHEP::HepRandomEngine* hre) {
    auto digiCollection = std::make_unique<FTLDigiCollection>();

    electronicsSim_.run(simHitAccumulator_,*digiCollection);

    e.put(std::move(digiCollection),digiCollection_);

    //release memory for next event
    resetSimHitDataAccumulator();
  }
    

  template<class SensorPhysics, class ElectronicsSim>
  void FTLDigitizer<SensorPhysics,ElectronicsSim>::beginRun(const edm::EventSetup & es) {
    if ( idealGeomWatcher_.check(es) ) {
      es.get<IdealGeometryRecord>().get(dddFTL_);
      { // force scope for the temporary nameless unordered_set
    std::unordered_set<DetId>().swap(validIds_);
      }
      // recalculate valid detids
      for( int zside = -1; zside <= 1; zside += 2 ) {     
    for( unsigned type = 1; type <= 2; ++type ) {
      for( unsigned izeta = 0; izeta < 1<<10; ++izeta ) {
        for( unsigned iphi = 0; iphi < 1<<10; ++iphi ) {
          
          if( dddFTL_->isValidXY(type, izeta, iphi) ) {
        validIds_.emplace( FastTimeDetId( type, 
                          izeta, 
                          iphi, 
                          zside ) );
          }
          
        }
      }
    }
      }      
      validIds_.reserve(validIds_.size());
    }
    deviceSim_.getEventSetup(es);
    electronicsSim_.getEventSetup(es);
  }
}


#endif