BTLBarDeviceSim.cc

Go to the documentation of this file.

#include "SimFastTiming/FastTimingCommon/interface/BTLBarDeviceSim.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/ForwardDetId/interface/MTDDetId.h"
#include "DataFormats/ForwardDetId/interface/BTLDetId.h"

#include "FWCore/Framework/interface/ESHandle.h"
#include "Geometry/CommonDetUnit/interface/GeomDetType.h"
#include "Geometry/MTDGeometryBuilder/interface/ProxyMTDTopology.h"
#include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h"
#include "Geometry/MTDNumberingBuilder/interface/GeometricTimingDet.h"

#include "CLHEP/Random/RandGaussQ.h"

BTLBarDeviceSim::BTLBarDeviceSim(const edm::ParameterSet& pset) : 
  geom_(nullptr),
  topo_(nullptr),
  bxTime_(pset.getParameter<double>("bxTime") ),
  LightYield_(pset.getParameter<double>("LightYield")),
  LightCollEff_(pset.getParameter<double>("LightCollectionEff")),
  LightCollSlopeR_(pset.getParameter<double>("LightCollectionSlopeR")),
  LightCollSlopeL_(pset.getParameter<double>("LightCollectionSlopeL")),
  PDE_(pset.getParameter<double>("PhotonDetectionEff")) { }

void BTLBarDeviceSim::getEventSetup(const edm::EventSetup& evs) {
 
  edm::ESHandle<MTDGeometry> geom;
  evs.get<MTDDigiGeometryRecord>().get(geom);
  geom_ = geom.product();

  edm::ESHandle<MTDTopology> mtdTopo;
  evs.get<MTDTopologyRcd>().get(mtdTopo);
  topo_ = mtdTopo.product();

}

void BTLBarDeviceSim::getHitsResponse(const std::vector<std::tuple<int,uint32_t,float> > &hitRefs, 
                      const edm::Handle<edm::PSimHitContainer> &hits,
                      mtd_digitizer::MTDSimHitDataAccumulator *simHitAccumulator,
                      CLHEP::HepRandomEngine *hre){

  //loop over sorted simHits
  for (auto const& hitRef: hitRefs) {

    const int hitidx   = std::get<0>(hitRef);
    const uint32_t id  = std::get<1>(hitRef);
    const MTDDetId detId(id);
    const PSimHit &hit = hits->at( hitidx );     
    
    // --- Safety check on the detector ID
    if ( detId.det()!=DetId::Forward || detId.mtdSubDetector()!=1 ) continue;

    if(id==0) continue; // to be ignored at RECO level                                                              

    BTLDetId btlid(detId);
    const int boundRef = btlid.modulesPerType(static_cast<BTLDetId::CrysLayout>(topo_->getMTDTopologyMode()));
    DetId geoId = BTLDetId(btlid.mtdSide(),btlid.mtdRR(),btlid.module()+boundRef*(btlid.modType()-1),0,1);
    const MTDGeomDet* thedet = geom_->idToDet(geoId);

    if( thedet == nullptr ) {
      throw cms::Exception("BTLBarDeviceSim") << "GeographicalID: " << std::hex
                          << geoId.rawId()
                          << " (" << detId.rawId()<< ") is invalid!" << std::dec
                          << std::endl;
    }
    const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(thedet->topology());
    const RectangularMTDTopology& topo = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());    
    // calculate the simhit row and column                                  
    const auto& pentry = hit.entryPoint();
    Local3DPoint simscaled(0.1*pentry.x(),0.1*pentry.y(),0.1*pentry.z()); // mm -> cm here is the switch
    // translate from crystal-local coordinates to module-local coordinates to get the row and column
    simscaled = topo.pixelToModuleLocalPoint(simscaled,btlid.row(topo.nrows()),btlid.column(topo.nrows()));
    const auto& thepixel = topo.pixel(simscaled);
    uint8_t row(thepixel.first), col(thepixel.second);

    if( btlid.row(topo.nrows()) != row || btlid.column(topo.nrows()) != col ) {
      edm::LogWarning("BTLBarDeviceSim")
    << "BTLDetId (row,column): (" << btlid.row(topo.nrows()) << ',' << btlid.column(topo.nrows()) <<") is not equal to "
    << "topology (row,column): (" << uint32_t(row) << ',' << uint32_t(col) <<"), overriding to detid";
      row = btlid.row(topo.nrows());
      col = btlid.column(topo.nrows()); 
    }


    // --- Store the detector element ID as a key of the MTDSimHitDataAccumulator map
    auto simHitIt = simHitAccumulator->emplace(mtd_digitizer::MTDCellId(id,row,col),
                           mtd_digitizer::MTDCellInfo()).first;

    // --- Get the simHit energy and convert it from MeV to photo-electrons
    float Npe = 1000.*hit.energyLoss()*LightYield_*LightCollEff_*PDE_;

    // --- Get the simHit time of arrival
    float toa = std::get<2>(hitRef);

    if ( toa > bxTime_ || toa < 0 ) //just consider BX==0
      continue;
 
    // --- Accumulate the energy of simHits in the same crystal for the BX==0
    // this is to simulate the charge integration in a 25 ns window
    (simHitIt->second).hit_info[0][0] += Npe;
    (simHitIt->second).hit_info[0][1] += Npe;
 
    double distR = 0.5*topo.pitch().second - 0.1*hit.localPosition().y();
    double distL = 0.5*topo.pitch().second + 0.1*hit.localPosition().y();
    
    // This is for the layout with bars along phi
    if ( topo_->getMTDTopologyMode() == (int) BTLDetId::CrysLayout::bar ||
     topo_->getMTDTopologyMode() == (int) BTLDetId::CrysLayout::barphiflat 
     )
      {
    distR = 0.5*topo.pitch().first - 0.1*hit.localPosition().x();
    distL = 0.5*topo.pitch().first + 0.1*hit.localPosition().x();
      }
    
    double tR = toa + LightCollSlopeR_*distR;
    double tL = toa + LightCollSlopeL_*distL;
   
    // --- Store the time of the first SimHit
    if ( (simHitIt->second).hit_info[1][0] == 0 
     || tR < (simHitIt->second).hit_info[1][0] )
      (simHitIt->second).hit_info[1][0] = tR;
    
    if ( (simHitIt->second).hit_info[1][1] == 0 
     || tL < (simHitIt->second).hit_info[1][1] )
      (simHitIt->second).hit_info[1][1] = tL;

  } // hitRef loop

}