HcalTB06BeamSD.cc

Go to the documentation of this file.

// File: HcalTB06BeamSD.cc
// Description: Sensitive Detector class for beam counters in TB06 setup

#include "SimG4Core/Notification/interface/TrackInformation.h"
#include "SimDataFormats/HcalTestBeam/interface/HcalTestBeamNumbering.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/DDLogicalPart.h"
#include "DetectorDescription/Core/interface/DDMaterial.h"
#include "DetectorDescription/Core/interface/DDSplit.h"
#include "DetectorDescription/Core/interface/DDValue.h"
#include "FWCore/Framework/interface/ESTransientHandle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "HcalTB06BeamSD.h"

#include "G4Step.hh"
#include "G4Track.hh"
#include "G4Material.hh"
#include "CLHEP/Units/GlobalSystemOfUnits.h"

//#define EDM_ML_DEBUG

HcalTB06BeamSD::HcalTB06BeamSD(const std::string& name,
                               const edm::EventSetup& es,
                               const SensitiveDetectorCatalog& clg,
                               edm::ParameterSet const& p,
                               const SimTrackManager* manager)
    : CaloSD(name, es, clg, p, manager) {
  // Values from NIM 80 (1970) 239-244: as implemented in Geant3
  edm::ParameterSet m_HC = p.getParameter<edm::ParameterSet>("HcalTB06BeamSD");
  useBirk = m_HC.getParameter<bool>("UseBirkLaw");
  birk1 = m_HC.getParameter<double>("BirkC1") * (g / (MeV * cm2));
  birk2 = m_HC.getParameter<double>("BirkC2");
  birk3 = m_HC.getParameter<double>("BirkC3");

  edm::LogInfo("HcalTB06BeamSD") << "HcalTB06BeamSD:: Use of Birks law is set "
                                 << "to " << useBirk << "  with three "
                                 << "constants kB = " << birk1 << ", C1 = " << birk2 << ", C2 = " << birk3;

  std::string attribute, value;

  edm::ESTransientHandle<DDCompactView> cpv;
  es.get<IdealGeometryRecord>().get(cpv);

  // Wire Chamber volume names
  attribute = "Volume";
  value = "WireChamber";
  DDSpecificsMatchesValueFilter filter1{DDValue(attribute, value, 0)};
  DDFilteredView fv1(*cpv, filter1);
  wcNames = getNames(fv1);
  edm::LogInfo("HcalTB06BeamSD") << "HcalTB06BeamSD:: Names to be tested for " << attribute << " = " << value << ": "
                                 << wcNames.size() << " paths";
  for (unsigned int i = 0; i < wcNames.size(); i++)
    edm::LogInfo("HcalTB06BeamSD") << "HcalTB06BeamSD:: (" << i << ") " << wcNames[i];

  //Material list for scintillator detector
  attribute = "ReadOutName";
  DDSpecificsMatchesValueFilter filter2{DDValue(attribute, name, 0)};
  DDFilteredView fv2(*cpv, filter2);
  bool dodet = fv2.firstChild();

  std::vector<G4String> matNames;
  std::vector<int> nocc;
  while (dodet) {
    const DDLogicalPart& log = fv2.logicalPart();
    matName = log.material().name().name();
    bool notIn = true;
    for (unsigned int i = 0; i < matNames.size(); i++) {
      if (matName == matNames[i]) {
        notIn = false;
        nocc[i]++;
      }
    }
    if (notIn) {
      matNames.push_back(matName);
      nocc.push_back(0);
    }
    dodet = fv2.next();
  }
  if (!matNames.empty()) {
    matName = matNames[0];
    int occ = nocc[0];
    for (unsigned int i = 0; i < matNames.size(); i++) {
      if (nocc[i] > occ) {
        occ = nocc[i];
        matName = matNames[i];
      }
    }
  } else {
    matName = "Not Found";
  }
  edm::LogInfo("HcalTB06BeamSD") << "HcalTB06BeamSD: Material name for " << attribute << " = " << name << ":"
                                 << matName;
}

HcalTB06BeamSD::~HcalTB06BeamSD() {}

double HcalTB06BeamSD::getEnergyDeposit(const G4Step* aStep) {
  double destep = aStep->GetTotalEnergyDeposit();
  double weight = 1;
  if (useBirk && matName == aStep->GetPreStepPoint()->GetMaterial()->GetName())
    weight *= getAttenuation(aStep, birk1, birk2, birk3);
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HcalTB06BeamSD") << "HcalTB06BeamSD: Detector "
                                     << aStep->GetPreStepPoint()->GetTouchable()->GetVolume()->GetName() << " weight "
                                     << weight;
#endif
  return weight * destep;
}

uint32_t HcalTB06BeamSD::setDetUnitId(const G4Step* aStep) {
  const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
  const G4VTouchable* touch = preStepPoint->GetTouchable();
  G4String name = preStepPoint->GetPhysicalVolume()->GetName();

  int det = 1;
  int lay = 0, x = 0, y = 0;
  if (!isItWireChamber(name)) {
    lay = (touch->GetReplicaNumber(0));
  } else {
    det = 2;
    lay = (touch->GetReplicaNumber(1));
    G4ThreeVector localPoint = setToLocal(preStepPoint->GetPosition(), touch);
    x = (int)(localPoint.x() / (0.2 * mm));
    y = (int)(localPoint.y() / (0.2 * mm));
  }

  return HcalTestBeamNumbering::packIndex(det, lay, x, y);
}

std::vector<G4String> HcalTB06BeamSD::getNames(DDFilteredView& fv) {
  std::vector<G4String> tmp;
  bool dodet = fv.firstChild();
  while (dodet) {
    const DDLogicalPart& log = fv.logicalPart();
    G4String name = log.name().name();
    bool ok = true;
    for (unsigned int i = 0; i < tmp.size(); i++)
      if (name == tmp[i])
        ok = false;
    if (ok)
      tmp.push_back(name);
    dodet = fv.next();
  }
  return tmp;
}

bool HcalTB06BeamSD::isItWireChamber(const G4String& name) {
  std::vector<G4String>::const_iterator it = wcNames.begin();
  for (; it != wcNames.end(); it++)
    if (name == *it)
      return true;
  return false;
}