HcalTB02SD.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:     HcalTestBeam
// Class  :     HcalTB02SD
//
// Implementation:
//     Sensitive Detector class for Hcal Test Beam 2002 detectors
//
// Original Author:
//         Created:  Sun 21 10:14:34 CEST 2006
//

// system include files

// user include files
#include "FWCore/Framework/interface/ESTransientHandle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "Geometry/Records/interface/IdealGeometryRecord.h"
#include "DetectorDescription/Core/interface/DDCompactView.h"
#include "DetectorDescription/Core/interface/DDFilter.h"
#include "DetectorDescription/Core/interface/DDFilteredView.h"
#include "DetectorDescription/Core/interface/DDSolid.h"
#include "DetectorDescription/Core/interface/DDSplit.h"
#include "DetectorDescription/Core/interface/DDValue.h"
#include "HcalTB02SD.h"
#include "HcalTB02HcalNumberingScheme.h"
#include "HcalTB02XtalNumberingScheme.h"

#include "G4Step.hh"
#include "G4Track.hh"
#include "G4VProcess.hh"

#include "G4SystemOfUnits.hh"

//
// constructors and destructor
//

HcalTB02SD::HcalTB02SD(const std::string& name,
                       const edm::EventSetup& es,
                       const SensitiveDetectorCatalog& clg,
                       edm::ParameterSet const& p,
                       const SimTrackManager* manager)
    : CaloSD(name, es, clg, p, manager), numberingScheme(nullptr) {
  edm::ParameterSet m_SD = p.getParameter<edm::ParameterSet>("HcalTB02SD");
  useBirk = m_SD.getUntrackedParameter<bool>("UseBirkLaw", false);
  birk1 = m_SD.getUntrackedParameter<double>("BirkC1", 0.013) * (g / (MeV * cm2));
  birk2 = m_SD.getUntrackedParameter<double>("BirkC2", 0.0568);
  birk3 = m_SD.getUntrackedParameter<double>("BirkC3", 1.75);
  useWeight = true;

  HcalTB02NumberingScheme* scheme = nullptr;
  if (name == "EcalHitsEB") {
    scheme = dynamic_cast<HcalTB02NumberingScheme*>(new HcalTB02XtalNumberingScheme());
    useBirk = false;
  } else if (name == "HcalHits") {
    scheme = dynamic_cast<HcalTB02NumberingScheme*>(new HcalTB02HcalNumberingScheme());
    useWeight = false;
  } else {
    edm::LogWarning("HcalTBSim") << "HcalTB02SD: ReadoutName " << name << " not supported\n";
  }

  if (scheme)
    setNumberingScheme(scheme);
  LogDebug("HcalTBSim") << "***************************************************"
                        << "\n"
                        << "*                                                 *"
                        << "\n"
                        << "* Constructing a HcalTB02SD  with name " << GetName() << "\n"
                        << "*                                                 *"
                        << "\n"
                        << "***************************************************";
  edm::LogInfo("HcalTBSim") << "HcalTB02SD:: Use of Birks law is set to      " << useBirk
                            << "        with three constants kB = " << birk1 << ", C1 = " << birk2
                            << ", C2 = " << birk3;

  if (useWeight)
    initMap(name, es);
}

HcalTB02SD::~HcalTB02SD() {
  if (numberingScheme)
    delete numberingScheme;
}

//
// member functions
//

double HcalTB02SD::getEnergyDeposit(const G4Step* aStep) {
  auto const preStepPoint = aStep->GetPreStepPoint();
  auto const& nameVolume = preStepPoint->GetPhysicalVolume()->GetName();

  // take into account light collection curve for crystals
  double weight = 1.;
  if (useWeight)
    weight *= curve_LY(nameVolume, preStepPoint);
  if (useBirk)
    weight *= getAttenuation(aStep, birk1, birk2, birk3);
  double edep = aStep->GetTotalEnergyDeposit() * weight;
  LogDebug("HcalTBSim") << "HcalTB02SD:: " << nameVolume << " Light Collection Efficiency " << weight
                        << " Weighted Energy Deposit " << edep / MeV << " MeV";
  return edep;
}

uint32_t HcalTB02SD::setDetUnitId(const G4Step* aStep) {
  return (numberingScheme == nullptr ? 0 : (uint32_t)(numberingScheme->getUnitID(aStep)));
}

void HcalTB02SD::setNumberingScheme(HcalTB02NumberingScheme* scheme) {
  if (scheme != nullptr) {
    edm::LogInfo("HcalTBSim") << "HcalTB02SD: updates numbering scheme for " << GetName();
    if (numberingScheme)
      delete numberingScheme;
    numberingScheme = scheme;
  }
}

void HcalTB02SD::initMap(const std::string& sd, const edm::EventSetup& es) {
  edm::ESTransientHandle<DDCompactView> cpv;
  es.get<IdealGeometryRecord>().get(cpv);

  G4String attribute = "ReadOutName";
  DDSpecificsMatchesValueFilter filter{DDValue(attribute, sd, 0)};
  DDFilteredView fv(*cpv, filter);
  fv.firstChild();

  bool dodet = true;
  while (dodet) {
    const DDSolid& sol = fv.logicalPart().solid();
    const std::vector<double>& paras = sol.parameters();
    G4String name = sol.name().name();
    LogDebug("HcalTBSim") << "HcalTB02SD::initMap (for " << sd << "): Solid " << name << " Shape " << sol.shape()
                          << " Parameter 0 = " << paras[0];
    if (sol.shape() == DDSolidShape::ddtrap) {
      double dz = 2 * paras[0];
      lengthMap.insert(std::pair<G4String, double>(name, dz));
    }
    dodet = fv.next();
  }
  LogDebug("HcalTBSim") << "HcalTB02SD: Length Table for " << attribute << " = " << sd << ":";
  std::map<G4String, double>::const_iterator it = lengthMap.begin();
  int i = 0;
  for (; it != lengthMap.end(); it++, i++) {
    LogDebug("HcalTBSim") << " " << i << " " << it->first << " L = " << it->second;
  }
}

double HcalTB02SD::curve_LY(const G4String& nameVolume, const G4StepPoint* stepPoint) {
  double weight = 1.;
  G4ThreeVector localPoint = setToLocal(stepPoint->GetPosition(), stepPoint->GetTouchable());
  double crlength = crystalLength(nameVolume);
  double dapd = 0.5 * crlength - localPoint.z();
  if (dapd >= -0.1 || dapd <= crlength + 0.1) {
    if (dapd <= 100.)
      weight = 1.05 - dapd * 0.0005;
  } else {
    edm::LogWarning("HcalTBSim") << "HcalTB02SD: light coll curve : wrong "
                                 << "distance to APD " << dapd << " crlength = " << crlength
                                 << " crystal name = " << nameVolume << " z of localPoint = " << localPoint.z()
                                 << " take weight = " << weight;
  }
  LogDebug("HcalTBSim") << "HcalTB02SD, light coll curve : " << dapd << " crlength = " << crlength
                        << " crystal name = " << nameVolume << " z of localPoint = " << localPoint.z()
                        << " take weight = " << weight;
  return weight;
}

double HcalTB02SD::crystalLength(const G4String& name) {
  double length = 230.;
  std::map<G4String, double>::const_iterator it = lengthMap.find(name);
  if (it != lengthMap.end())
    length = it->second;
  return length;
}