AHCalSD.cc

Go to the documentation of this file.

#include "SimG4CMS/HGCalTestBeam/interface/AHCalSD.h"
#include "SimG4CMS/HGCalTestBeam/interface/AHCalDetId.h"
#include "SimG4Core/Notification/interface/TrackInformation.h"
#include "Geometry/HGCalCommonData/interface/AHCalParameters.h"

#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"

#include "G4LogicalVolume.hh"
#include "G4LogicalVolumeStore.hh"
#include "G4ParticleTable.hh"
#include "G4Track.hh"
#include "G4VProcess.hh"

#include "G4PhysicalConstants.hh"
#include <CLHEP/Units/SystemOfUnits.h>

#include <iomanip>
#include <map>

//#define EDM_ML_DEBUG

AHCalSD::AHCalSD(const std::string& name,
                 const SensitiveDetectorCatalog& clg,
                 edm::ParameterSet const& p,
                 const SimTrackManager* manager)
    : CaloSD(name,
             clg,
             p,
             manager,
             (float)(p.getParameter<edm::ParameterSet>("AHCalSD").getParameter<double>("TimeSliceUnit")),
             p.getParameter<edm::ParameterSet>("AHCalSD").getParameter<bool>("IgnoreTrackID")) {
  edm::ParameterSet m_HC = p.getParameter<edm::ParameterSet>("AHCalSD");
  useBirk = m_HC.getParameter<bool>("UseBirkLaw");
  birk1 = m_HC.getParameter<double>("BirkC1") * (CLHEP::g / (CLHEP::MeV * CLHEP::cm2));
  birk2 = m_HC.getParameter<double>("BirkC2");
  birk3 = m_HC.getParameter<double>("BirkC3");
  eminHit = m_HC.getParameter<double>("EminHit") * CLHEP::MeV;

  edm::LogVerbatim("HcalSim") << "AHCalSD::  Use of Birks law is set to      " << useBirk
                              << "  with three constants kB = " << birk1 << ", C1 = " << birk2 << ", C2 = " << birk3
                              << "\nAHCalSD:: Threshold for storing"
                              << " hits: " << eminHit;
}

double AHCalSD::getEnergyDeposit(const G4Step* aStep) {
  double destep = aStep->GetTotalEnergyDeposit();
  double wt2 = aStep->GetTrack()->GetWeight();
  double weight = (wt2 > 0.0) ? wt2 : 1.0;
#ifdef EDM_ML_DEBUG
  double weight0 = weight;
#endif
  if (useBirk)
    weight *= getAttenuation(aStep, birk1, birk2, birk3);
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HcalSim") << "AHCalSD: weight " << weight0 << " " << weight;
#endif
  double edep = weight * destep;
  return edep;
}

uint32_t AHCalSD::setDetUnitId(const G4Step* aStep) {
  const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
  const G4VTouchable* touch = preStepPoint->GetTouchable();

  int depth = (touch->GetReplicaNumber(1));
  int incol = ((touch->GetReplicaNumber(0)) % AHCalParameters::kColumn_);
  int inrow = ((touch->GetReplicaNumber(0)) / AHCalParameters::kColumn_) % AHCalParameters::kRow_;
  int jncol = ((touch->GetReplicaNumber(0)) / AHCalParameters::kRowColumn_) % AHCalParameters::kSign_;
  int jnrow = ((touch->GetReplicaNumber(0)) / AHCalParameters::kSignRowColumn_) % AHCalParameters::kSign_;
  int col = (jncol == 0) ? incol : -incol;
  int row = (jnrow == 0) ? inrow : -inrow;
  uint32_t index = AHCalDetId(row, col, depth).rawId();
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HcalSim") << "AHCalSD: det = " << HcalOther << " depth = " << depth << " row = " << row
                              << " column = " << col << " packed index = 0x" << std::hex << index << std::dec;
  bool flag = unpackIndex(index, row, col, depth);
  edm::LogVerbatim("HcalSim") << "Results from unpacker for 0x" << std::hex << index << std::dec << " Flag " << flag
                              << " Row " << row << " Col " << col << " Depth " << depth;
#endif
  return index;
}

bool AHCalSD::unpackIndex(const uint32_t& idx, int& row, int& col, int& depth) {
  DetId gen(idx);
  HcalSubdetector subdet = (HcalSubdetector(gen.subdetId()));
  bool rcode = (gen.det() == DetId::Hcal && subdet != HcalOther);
  row = col = depth = 0;
  if (rcode) {
    row = AHCalDetId(idx).irow();
    col = AHCalDetId(idx).icol();
    depth = AHCalDetId(idx).depth();
  }
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HcalSim") << "AHCalSD: packed index = 0x" << std::hex << idx << std::dec << " Row " << row
                              << " Column " << col << " Depth " << depth << " OK " << rcode;
#endif
  return rcode;
}

bool AHCalSD::filterHit(CaloG4Hit* aHit, double time) {
  return ((time <= tmaxHit) && (aHit->getEnergyDeposit() > eminHit));
}