HCalSD.cc

Go to the documentation of this file.

// File: HCalSD.cc
// Description: Sensitive Detector class for calorimeters

#include "SimG4CMS/Calo/interface/HCalSD.h"
#include "SimG4CMS/Calo/interface/HcalTestNumberingScheme.h"
#include "SimG4CMS/Calo/interface/HFFibreFiducial.h"
#include "SimG4Core/Notification/interface/TrackInformation.h"
#include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.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/DDValue.h"

#include "Geometry/Records/interface/HcalSimNumberingRecord.h"
#include "FWCore/Framework/interface/ESTransientHandle.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "CondFormats/DataRecord/interface/HBHEDarkeningRecord.h"

#include "G4LogicalVolumeStore.hh"
#include "G4LogicalVolume.hh"
#include "G4Step.hh"
#include "G4Track.hh"

#include "G4SystemOfUnits.hh"
#include "G4PhysicalConstants.hh"
#include "Randomize.hh"

#include <iostream>
#include <fstream>
#include <iomanip>
#include <sstream>

//#define EDM_ML_DEBUG
//#define plotDebug

#ifdef plotDebug
#include <TH1F.h>
#endif

HCalSD::HCalSD(const std::string& name, const DDCompactView & cpv,
               const SensitiveDetectorCatalog & clg,
               edm::ParameterSet const & p, const SimTrackManager* manager) : 
  CaloSD(name, cpv, clg, p, manager,
         (float)(p.getParameter<edm::ParameterSet>("HCalSD").getParameter<double>("TimeSliceUnit")),
         p.getParameter<edm::ParameterSet>("HCalSD").getParameter<bool>("IgnoreTrackID")), 
  hcalConstants(nullptr), m_HBDarkening(nullptr), m_HEDarkening(nullptr), isHF(false),
  weight_(1.0), depth_(1) {

  numberingFromDDD.reset(nullptr); 
  numberingScheme.reset(nullptr); 
  showerLibrary.reset(nullptr); 
  hfshower.reset(nullptr); 
  showerParam.reset(nullptr);
  showerPMT.reset(nullptr); 
  showerBundle.reset(nullptr);
  m_HFDarkening.reset(nullptr);
  m_HcalTestNS.reset(nullptr);

  //static SimpleConfigurable<double> bk1(0.013, "HCalSD:BirkC1");
  //static SimpleConfigurable<double> bk2(0.0568,"HCalSD:BirkC2");
  //static SimpleConfigurable<double> bk3(1.75,  "HCalSD:BirkC3");
  // Values from NIM 80 (1970) 239-244: as implemented in Geant3

  edm::ParameterSet m_HC = p.getParameter<edm::ParameterSet>("HCalSD");
  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");
  useShowerLibrary = m_HC.getParameter<bool>("UseShowerLibrary");
  useParam         = m_HC.getParameter<bool>("UseParametrize");
  testNumber       = m_HC.getParameter<bool>("TestNumberingScheme");
  neutralDensity   = m_HC.getParameter<bool>("doNeutralDensityFilter");
  usePMTHit        = m_HC.getParameter<bool>("UsePMTHits");
  betaThr          = m_HC.getParameter<double>("BetaThreshold");
  eminHitHB        = m_HC.getParameter<double>("EminHitHB")*MeV;
  eminHitHE        = m_HC.getParameter<double>("EminHitHE")*MeV;
  eminHitHO        = m_HC.getParameter<double>("EminHitHO")*MeV;
  eminHitHF        = m_HC.getParameter<double>("EminHitHF")*MeV;
  useFibreBundle   = m_HC.getParameter<bool>("UseFibreBundleHits");
  deliveredLumi    = m_HC.getParameter<double>("DelivLuminosity");
  agingFlagHB      = m_HC.getParameter<bool>("HBDarkening");
  agingFlagHE      = m_HC.getParameter<bool>("HEDarkening");
  bool agingFlagHF = m_HC.getParameter<bool>("HFDarkening");
  useHF            = m_HC.getUntrackedParameter<bool>("UseHF",true);
  bool forTBH2     = m_HC.getUntrackedParameter<bool>("ForTBH2",false);
  useLayerWt       = m_HC.getUntrackedParameter<bool>("UseLayerWt",false);
  std::string file = m_HC.getUntrackedParameter<std::string>("WtFile","None");
  testNS_          = m_HC.getUntrackedParameter<bool>("TestNS",false);
  edm::ParameterSet m_HF  = p.getParameter<edm::ParameterSet>("HFShower");
  applyFidCut             = m_HF.getParameter<bool>("ApplyFiducialCut");

#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HcalSim") << "***************************************************" 
                  << "\n"
                  << "* Constructing a HCalSD  with name " << name << "\n"
                  << "\n"
                  << "***************************************************";
#endif
  edm::LogVerbatim("HcalSim") << "HCalSD:: Use of HF code is set to " << useHF
                          << "\nUse of shower parametrization set to "
                          << useParam << "\nUse of shower library is set to " 
                          << useShowerLibrary << "\nUse PMT Hit is set to "
                          << usePMTHit << " with beta Threshold "<< betaThr
                          << "\nUSe of FibreBundle Hit set to "<<useFibreBundle
                          << "\n         Use of Birks law is set to      " 
                          << useBirk << "  with three constants kB = "
                          << birk1 << ", C1 = " << birk2 << ", C2 = " << birk3;
  edm::LogVerbatim("HcalSim") << "HCalSD:: Suppression Flag " << suppressHeavy
                          << " protons below " << kmaxProton << " MeV,"
                          << " neutrons below " << kmaxNeutron << " MeV and"
                          << " ions below " << kmaxIon << " MeV\n"
                          << "         Threshold for storing hits in HB: "
                          << eminHitHB << " HE: " << eminHitHE << " HO: "
                          << eminHitHO << " HF: " << eminHitHF << "\n"
                          << "Delivered luminosity for Darkening " 
                          << deliveredLumi << " Flag (HE) " << agingFlagHE
                          << " Flag (HB) " << agingFlagHB
                          << " Flag (HF) " << agingFlagHF << "\n"
                          << "Application of Fiducial Cut " << applyFidCut
                          << "Flag for test number|neutral density filter "
                          << testNumber << " " << neutralDensity;

  HcalNumberingScheme* scheme;
  if (testNumber || forTBH2) {
    scheme = dynamic_cast<HcalNumberingScheme*>(new HcalTestNumberingScheme(forTBH2));
  } else { 
    scheme = new HcalNumberingScheme();
  }
  setNumberingScheme(scheme);

  // always call getFromLibrary() method to identify HF region
  setParameterized(true);

  const G4LogicalVolumeStore * lvs = G4LogicalVolumeStore::GetInstance();
  //  std::vector<G4LogicalVolume *>::const_iterator lvcite;
  const G4LogicalVolume* lv;
  std::string attribute, value;

  if (useHF) {
    if (useParam) {
      showerParam.reset(new HFShowerParam(name, cpv, p));
    }  else {
      if (useShowerLibrary) { showerLibrary.reset(new HFShowerLibrary(name, cpv, p)); }
      hfshower.reset(new HFShower(name, cpv, p, 0));
    }

    // HF volume names
    attribute = "Volume";
    value     = "HF";
    DDSpecificsMatchesValueFilter filter0{DDValue(attribute,value,0)};
    DDFilteredView fv0(cpv,filter0);
    hfNames = getNames(fv0);
    fv0.firstChild();
    DDsvalues_type sv0(fv0.mergedSpecifics());
    std::vector<double> temp =  getDDDArray("Levels",sv0);

    unsigned int nhf = hfNames.size();
    std::stringstream ss;
    ss << "HCalSD: Names to be tested for " << attribute 
       << " = " << value << " has " << nhf << " elements";
    for (unsigned int i=0; i < nhf; ++i) {
      G4String namv = hfNames[i];
      lv            = nullptr;
      for(auto lvol : *lvs) {
        if(lvol->GetName() == namv) {
          lv = lvol;
          break;
        }
      }
      hfLV.push_back(lv);
      int level = static_cast<int>(temp[i]);
      hfLevels.push_back(level);
      ss << "\n        HF[" << i << "] = " << namv
         << " LV " << hfLV[i] << " at level " << hfLevels[i];
    }
    edm::LogVerbatim("HcalSim") << ss.str();
  
    // HF Fibre volume names
    fillLogVolumeVector(attribute, "HFFibre", cpv, fibreLV, fibreNames);
    std::vector<G4String> tempNames;
    fillLogVolumeVector(attribute, "HFPMT", cpv, pmtLV, tempNames);
    fillLogVolumeVector(attribute, "HFFibreBundleStraight", cpv, fibre1LV, tempNames);
    fillLogVolumeVector(attribute, "HFFibreBundleConical", cpv, fibre2LV, tempNames);
  }

  //Material list for HB/HE/HO sensitive detectors
  const G4MaterialTable * matTab = G4Material::GetMaterialTable();
  std::vector<G4Material*>::const_iterator matite;
  attribute = "OnlyForHcalSimNumbering"; 
  DDSpecificsHasNamedValueFilter filter2{attribute};
  DDFilteredView fv2(cpv,filter2);
  bool dodet = fv2.firstChild();
  DDsvalues_type sv(fv2.mergedSpecifics());

  while (dodet) {
    const DDLogicalPart & log = fv2.logicalPart();
    G4String namx = log.name().name();
    if (!isItHF(namx) && !isItFibre(namx)) {
      bool notIn = true;
      for (unsigned int i=0; i<matNames.size(); ++i) {
        if (!strcmp(matNames[i].c_str(),log.material().name().name().c_str())){
          notIn = false;
          break;
        }
      }
      if (notIn) {
        namx = log.material().name().name();
        matNames.push_back(namx);
        const G4Material* mat = nullptr;
        for (matite = matTab->begin(); matite != matTab->end(); ++matite) {
          if ((*matite)->GetName() == namx) {
            mat = (*matite);
            break;
          }
        }
        materials.push_back(mat);
      }
    }
    dodet = fv2.next();
  }

  std::stringstream sss;
  for (unsigned int i=0; i<matNames.size(); ++i) {
    if(i/10*10 == i) { sss << "\n"; }
    sss << "  "  << matNames[i];
  }
  edm::LogVerbatim("HcalSim") << "HCalSD: Material names for " << attribute 
                          << " = " << name << ":" << sss.str();

  if (useLayerWt) { readWeightFromFile(file); }

  for (int i=0;  i<9; ++i) { hit_[i] = time_[i]= dist_[i] = nullptr; }
  hzvem = hzvhad = nullptr;

  if (agingFlagHF) {
    m_HFDarkening.reset(new HFDarkening(m_HC.getParameter<edm::ParameterSet>("HFDarkeningParameterBlock")));
  }
#ifdef plotDebug
  edm::Service<TFileService> tfile;

  if ( tfile.isAvailable() ) {
    static const char * const labels[] = {"HB", "HE", "HO", "HF Absorber", "HF PMT",
                                          "HF Absorber Long", "HF Absorber Short",
                                          "HF PMT Long", "HF PMT Short"};
    TFileDirectory hcDir = tfile->mkdir("ProfileFromHCalSD");
    char name[20], title[60];
    for (int i=0; i<9; ++i) {
      sprintf (title, "Hit energy in %s", labels[i]);
      sprintf (name, "HCalSDHit%d", i);
      hit_[i] = hcDir.make<TH1F>(name, title, 2000, 0., 2000.);
      sprintf (title, "Energy (MeV)");
      hit_[i]->GetXaxis()->SetTitle(title);
      hit_[i]->GetYaxis()->SetTitle("Hits");
      sprintf (title, "Time of the hit in %s", labels[i]);
      sprintf (name, "HCalSDTime%d", i);
      time_[i] = hcDir.make<TH1F>(name, title, 2000, 0., 2000.);
      sprintf (title, "Time (ns)");
      time_[i]->GetXaxis()->SetTitle(title);
      time_[i]->GetYaxis()->SetTitle("Hits");
      sprintf (title, "Longitudinal profile in %s", labels[i]);
      sprintf (name, "HCalSDDist%d", i);
      dist_[i] = hcDir.make<TH1F>(name, title, 2000, 0., 2000.);
      sprintf (title, "Distance (mm)");
      dist_[i]->GetXaxis()->SetTitle(title);
      dist_[i]->GetYaxis()->SetTitle("Hits");
    }
    if (useHF && (!useParam)) {
      hzvem  = hcDir.make<TH1F>("hzvem", "Longitudinal Profile (EM Part)",330,0.0,1650.0);
      hzvem->GetXaxis()->SetTitle("Longitudinal Profile (EM Part)");
      hzvhad = hcDir.make<TH1F>("hzvhad","Longitudinal Profile (Had Part)",330,0.0,1650.0);
      hzvhad->GetXaxis()->SetTitle("Longitudinal Profile (Hadronic Part)");
    }
  }
#endif
}

void HCalSD::fillLogVolumeVector(const std::string& attribute, const std::string& value, 
                                 const DDCompactView& cpv,
                                 std::vector<const G4LogicalVolume*>& lvvec,
                                 std::vector<G4String>& lvnames) {

  const G4LogicalVolumeStore * lvs = G4LogicalVolumeStore::GetInstance();
  const G4LogicalVolume* lv;

  DDSpecificsMatchesValueFilter filter{DDValue(attribute,value,0)};
  DDFilteredView fv(cpv,filter);
  lvnames = getNames(fv);

  unsigned int nvol = lvnames.size();
  std::stringstream ss;
  ss << "HCalSD: " << nvol << " names to be tested for " << attribute << " <" << value << ">:";
  for (unsigned int i=0; i<nvol; ++i) {
    G4String namv = lvnames[i];
    lv = nullptr;
    for (auto lvol : *lvs) {
      if (lvol->GetName() == namv) {
        lv = lvol;
        break;
      }
    }
    lvvec.push_back(lv);
    if(i/10*10 == i) { ss << "\n"; }
    ss << "  " << namv;
  }
  edm::LogVerbatim("HcalSim") << ss.str();
}

bool HCalSD::getFromLibrary(const G4Step * aStep) {

  auto const track = aStep->GetTrack();
  depth_ = (aStep->GetPreStepPoint()->GetTouchable()->GetReplicaNumber(0))%10;
  weight_= 1.0;
  bool kill(false);
  isHF = isItHF(aStep);
  if(isHF) { 
    if (m_HFDarkening) {
      G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
      const double invcm = 1./CLHEP::cm;
      double r = hitPoint.perp()*invcm;
      double z = std::abs(hitPoint.z())*invcm;
      double dose_acquired = 0.;
      if (z>=HFDarkening::lowZLimit && z <= HFDarkening::upperZLimit) {
        unsigned int hfZLayer = (unsigned int)((z - HFDarkening::lowZLimit)/5);
        if (hfZLayer >= HFDarkening::upperZLimit) hfZLayer = (HFDarkening::upperZLimit-1);
        float normalized_lumi = m_HFDarkening->int_lumi(deliveredLumi);
        for (int i = hfZLayer; i != HFDarkening::numberOfZLayers; ++i) {
          dose_acquired = m_HFDarkening->dose(i,r);
          weight_ *= m_HFDarkening->degradation(normalized_lumi*dose_acquired);
        }
      }
#ifdef EDM_ML_DEBUG
      edm::LogVerbatim("HcalSim") << "HCalSD::getFromLibrary: HFLumiDarkening at "
                  << "r= " << r << ", z= " << z << " Dose= "
                  << dose_acquired << " weight= " << weight_;
#endif
    }

    if (useParam) {
      getFromParam(aStep, kill);
#ifdef EDM_ML_DEBUG
      G4String nameVolume = lv->GetName();
      edm::LogVerbatim("HcalSim") << "HCalSD: " << getNumberOfHits()
                  << " hits from parametrization in "
                  << nameVolume << " for Track " 
                  << track->GetTrackID() <<" ("
                  << track->GetDefinition()->GetParticleName() 
                  <<")";
#endif
    } else if (useShowerLibrary && !G4TrackToParticleID::isMuon(track)) {
      if(G4TrackToParticleID::isGammaElectronPositron(track) || 
     G4TrackToParticleID::isStableHadronIon(track) ) {
#ifdef EDM_ML_DEBUG
    auto nameVolume = 
      aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
    edm::LogVerbatim("HcalSim") << "HCalSD: Starts shower library from " 
                << nameVolume << " for Track " 
                << track->GetTrackID() << " ("
                << track->GetDefinition()->GetParticleName()
                << ")";

#endif
    getFromHFLibrary(aStep, kill);
      }
    }
  }
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HcalSim") << "HCalSD::getFromLibrary ID= " 
                          << track->GetTrackID() << " ("
                          << track->GetDefinition()->GetParticleName() 
                          << ") kill= " << kill << " weight= " << weight_ 
                          << " depth= " << depth_ << " isHF: " << isHF;
#endif
  return kill;
} 

double HCalSD::getEnergyDeposit(const G4Step* aStep) {
  double destep(0.0);
  auto const lv = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
  auto const theTrack = aStep->GetTrack();

  if(isHF) {
    if (useShowerLibrary && G4TrackToParticleID::isMuon(theTrack) && isItFibre(lv)) {
#ifdef EDM_ML_DEBUG
      edm::LogVerbatim("HcalSim") << "HCalSD: Hit at Fibre in LV " << lv->GetName()
                              << " for track " 
                  << aStep->GetTrack()->GetTrackID() <<" ("
                              << aStep->GetTrack()->GetDefinition()->GetParticleName()
                  << ")";
#endif
      hitForFibre(aStep);
    }
    return destep;
  }

  if (isItPMT(lv)) {
    if(usePMTHit && showerPMT) { getHitPMT(aStep); }
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HcalSim") << "HCalSD: Hit from PMT parametrization in LV " 
                            <<  lv->GetName() << " for Track " 
                            << aStep->GetTrack()->GetTrackID() << " ("
                            << aStep->GetTrack()->GetDefinition()->GetParticleName()
                << ")";
#endif
    return destep;

  } else if (isItStraightBundle(lv)) {
    if(useFibreBundle && showerBundle) { getHitFibreBundle(aStep, false); }
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HcalSim") << "HCalSD: Hit from straight FibreBundle in LV: "
                << lv->GetName() << " for track " 
                << aStep->GetTrack()->GetTrackID() << " ("
                << aStep->GetTrack()->GetDefinition()->GetParticleName()
                << ")";
#endif
    return destep;

  } else if(isItConicalBundle(lv)) {
    if(useFibreBundle && showerBundle) { getHitFibreBundle(aStep, true); }
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HcalSim") << "HCalSD: Hit from conical FibreBundle PV: "
                << lv->GetName() << " for track " 
                << aStep->GetTrack()->GetTrackID() << " ("
                << aStep->GetTrack()->GetDefinition()->GetParticleName()
                << ")";
#endif
    return destep;
  }

  // normal hit
  destep = aStep->GetTotalEnergyDeposit();

  const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
  uint32_t detid = setDetUnitId(aStep);
  int det(0), ieta(0), phi(0), z(0), lay, depth(-1);
  if (testNumber) {
    HcalTestNumbering::unpackHcalIndex(detid,det,z,depth,ieta,phi,lay);
    if (z==0) { z = -1; }
  } else {
    HcalDetId hcid(detid);
    det  = hcid.subdetId();
    ieta = hcid.ietaAbs();
    phi  = hcid.iphi();
    z    = hcid.zside();
  }
  lay   = (touch->GetReplicaNumber(0)/10)%100 + 1;
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HcalSim") << "HCalSD: det: " << det << " ieta: "<< ieta 
                          << " iphi: " << phi << " zside " << z << "  lay: " 
                          << lay-2;
#endif 
  if (depth_==0 && (det==1 || det==2) && ((!testNumber) || neutralDensity))
    weight_ = hcalConstants->getLayer0Wt(det,phi,z);
  if (useLayerWt) {
    G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
    weight_ = layerWeight(det+2, hitPoint, depth_, lay);
  }

  if (m_HBDarkening && det == 1) {
    double dweight = m_HBDarkening->degradation(deliveredLumi,ieta,lay);
    weight_ *= dweight;
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HcalSim") << "HCalSD: HB Lumi: " << deliveredLumi
                            << " coefficient = " << dweight << " Weight= "
                << weight_;
#endif  
  }

  if (m_HEDarkening && det == 2) {
    double dweight = m_HEDarkening->degradation(deliveredLumi,ieta,lay);
    weight_ *= dweight;
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HcalSim") << "HCalSD: HB Lumi: " << deliveredLumi
                            << " coefficient = " << dweight << " Weight= "
                << weight_;
#endif  
  }

  if (suppressHeavy) {
    TrackInformation * trkInfo = (TrackInformation *)(theTrack->GetUserInformation());
    if (trkInfo) {
      int pdg = theTrack->GetDefinition()->GetPDGEncoding();
      if (!(trkInfo->isPrimary())) {         // Only secondary particles
        double ke = theTrack->GetKineticEnergy();
        if ( pdg/1000000000 == 1 && (pdg/10000)%100 > 0 &&
             (pdg/10)%100 > 0    && ke <kmaxIon   ) weight_ = 0;
        if ((pdg == 2212) && (ke < kmaxProton))     weight_ = 0;
        if ((pdg == 2112) && (ke < kmaxNeutron))    weight_ = 0;
      }
    }
  }
  if (useBirk) {
    const G4Material* mat = aStep->GetPreStepPoint()->GetMaterial();
    if (isItScintillator(mat)) weight_ *= getAttenuation(aStep, birk1, birk2, birk3);
  }
  double wt1 = getResponseWt(theTrack);
  double wt2 = theTrack->GetWeight();
  double edep = weight_*wt1*destep;
  if (wt2 > 0.0) { edep *= wt2; }
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HcalSim") 
    << "HCalSD: edep= " << edep << " Det: " << det+2 << " depth= " << depth_
    << " weight= " << weight_ << " wt1= " << wt1 << " wt2= " << wt2; 
#endif
  return edep;
}

uint32_t HCalSD::setDetUnitId(const G4Step * aStep) { 

  auto const prePoint  = aStep->GetPreStepPoint(); 
  auto const touch     = prePoint->GetTouchable();
  const G4ThreeVector& hitPoint = prePoint->GetPosition();

  int depth = (touch->GetReplicaNumber(0))%10 + 1;
  int lay   = (touch->GetReplicaNumber(0)/10)%100 + 1;
  int det   = (touch->GetReplicaNumber(1))/1000;

  return setDetUnitId (det, hitPoint, depth, lay);
}

void HCalSD::setNumberingScheme(HcalNumberingScheme * scheme) {
  if (scheme != nullptr) {
    edm::LogVerbatim("HcalSim") << "HCalSD: updates numbering scheme for "
                << GetName();
    numberingScheme.reset(scheme);
  }
}

void HCalSD::update(const BeginOfJob * job) {

  const edm::EventSetup* es = (*job)();
  edm::ESHandle<HcalDDDSimConstants>    hdc;
  es->get<HcalSimNumberingRecord>().get(hdc);
  if (hdc.isValid()) {
    hcalConstants = hdc.product();
  } else {
    edm::LogError("HcalSim") << "HCalSD : Cannot find HcalDDDSimConstant";
    throw cms::Exception("Unknown", "HCalSD") << "Cannot find HcalDDDSimConstant" << "\n";
  }

  numberingFromDDD.reset(new HcalNumberingFromDDD(hcalConstants));

  //Special Geometry parameters
  gpar = hcalConstants->getGparHF();
  std::stringstream sss;
  for (unsigned int ig=0; ig<gpar.size(); ig++) {
    sss << "\n         gpar[" << ig << "] = " << gpar[ig]/cm << " cm";
  }
  edm::LogVerbatim("HcalSim") << "Maximum depth for HF " 
              << hcalConstants->getMaxDepth(2)
                          << gpar.size()<< " gpar (cm)" << sss.str();
  //Test Hcal Numbering Scheme
  if (testNS_) m_HcalTestNS.reset(new HcalTestNS(es));

  if (agingFlagHB) {
    edm::ESHandle<HBHEDarkening> hdark;
    es->get<HBHEDarkeningRecord>().get("HB",hdark);
    m_HBDarkening = &*hdark;
  }
  if (agingFlagHE) {
    edm::ESHandle<HBHEDarkening> hdark;
    es->get<HBHEDarkeningRecord>().get("HE",hdark);
    m_HEDarkening = &*hdark;
  }
}

void HCalSD::initRun() {
  if (showerLibrary.get()) showerLibrary.get()->initRun(nullptr, hcalConstants);
  if (showerParam.get())   showerParam.get()->initRun(hcalConstants);
  if (hfshower.get())      hfshower.get()->initRun(hcalConstants);
  if (showerPMT.get())     showerPMT.get()->initRun(hcalConstants);
  if (showerBundle.get())  showerBundle.get()->initRun(hcalConstants);
}

bool HCalSD::filterHit(CaloG4Hit* aHit, double time) {
  double threshold=0;
  DetId theId(aHit->getUnitID());
  switch (theId.subdetId()) {
  case HcalBarrel:
    threshold = eminHitHB; break;
  case HcalEndcap:
    threshold = eminHitHE; break;
  case HcalOuter:
    threshold = eminHitHO; break;
  case HcalForward:
    threshold = eminHitHF; break;
  default:
    break;
  }
  return ((time <= tmaxHit) && (aHit->getEnergyDeposit() > threshold));
}

uint32_t HCalSD::setDetUnitId (int det, const G4ThreeVector& pos, int depth, int lay=1) { 
  uint32_t id = 0;
  if (numberingFromDDD.get()) {
    //get the ID's as eta, phi, depth, ... indices
    HcalNumberingFromDDD::HcalID tmp = numberingFromDDD.get()->unitID(det, pos, depth, lay);
    id = setDetUnitId(tmp);
  }
  return id;
}

uint32_t HCalSD::setDetUnitId (HcalNumberingFromDDD::HcalID& tmp) { 
  modifyDepth(tmp);
  uint32_t id = (numberingScheme.get()) ? numberingScheme.get()->getUnitID(tmp) : 0;
  if ((!testNumber) && m_HcalTestNS.get()) {
    bool ok = m_HcalTestNS.get()->compare(tmp,id);
    if (!ok)
      edm::LogWarning("HcalSim") << "Det ID from HCalSD " << HcalDetId(id) 
                                 << " " << std::hex << id << std::dec 
                                 << " does not match one from relabller for " 
                                 << tmp.subdet << ":" << tmp.etaR << ":"
                                 << tmp.phi << ":" << tmp.phis << ":" 
                                 << tmp.depth << ":" << tmp.lay << std::endl;
  }
  return id;
}

std::vector<double> HCalSD::getDDDArray(const std::string & str,
                                        const DDsvalues_type & sv) {
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HcalSim") << "HCalSD:getDDDArray called for " << str;
#endif
  DDValue value(str);
  if (DDfetch(&sv,value)) {
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HcalSim") << value;
#endif
    const std::vector<double> & fvec = value.doubles();
    int nval = fvec.size();
    if (nval < 1) {
      edm::LogError("HcalSim") << "HCalSD : # of " << str << " bins " << nval
          << " < 2 ==> illegal";
      throw cms::Exception("Unknown", "HCalSD") << "nval < 2 for array " << str << "\n";
    }
    
    return fvec;
  } else {
    edm::LogError("HcalSim") << "HCalSD :  cannot get array " << str;
    throw cms::Exception("Unknown", "HCalSD") << "cannot get array " << str << "\n";
  }
}

std::vector<G4String> HCalSD::getNames(DDFilteredView& fv) {

  std::vector<G4String> tmp;
  bool dodet = fv.firstChild();
  while (dodet) {
    const DDLogicalPart & log = fv.logicalPart();
    bool ok = true;

    for (unsigned int i=0; i<tmp.size(); ++i) {
      if (!strcmp(tmp[i].c_str(), log.name().name().c_str())) {
        ok = false;
        break;
      }
    }
    if (ok) tmp.push_back(log.name().name());
    dodet = fv.next();
  }
  return tmp;
}

bool HCalSD::isItHF(const G4Step * aStep) {
  const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
  int levels = (touch->GetHistoryDepth()) + 1;
  for (unsigned int it=0; it < hfNames.size(); ++it) {
    if (levels >= hfLevels[it]) {
      const G4LogicalVolume* lv = touch->GetVolume(levels-hfLevels[it])->GetLogicalVolume();
      if (lv == hfLV[it]) return true;
    }
  }
  return false;
}

bool HCalSD::isItHF (const G4String& name) {
  for (auto nam : hfNames) if (name == nam) { return true; }
  return false;
}

bool HCalSD::isItFibre (const G4LogicalVolume* lv) {
  for (auto lvol : fibreLV) if (lv == lvol) { return true; }
  return false;
}

bool HCalSD::isItFibre (const G4String& name) {
  for (auto nam : fibreNames) if (name == nam) { return true; }
  return false;
}

bool HCalSD::isItPMT (const G4LogicalVolume* lv) {
  for (auto lvol : pmtLV) if (lv == lvol) { return true; }
  return false;
}

bool HCalSD::isItStraightBundle (const G4LogicalVolume* lv) {
  for (auto lvol : fibre1LV) if (lv == lvol) { return true; }
  return false;
}

bool HCalSD::isItConicalBundle (const G4LogicalVolume* lv) {
  for (auto lvol : fibre2LV) if (lv == lvol) { return true; }
  return false;
}

bool HCalSD::isItScintillator (const G4Material* mat) {
  for (auto amat : materials) if (amat == mat) { return true; }
  return false;
}

bool HCalSD::isItinFidVolume (const G4ThreeVector& hitPoint) {
  bool flag = true;
  if (applyFidCut) {
    int npmt = HFFibreFiducial::PMTNumber(hitPoint);
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HcalSim") << "HCalSD::isItinFidVolume:#PMT= " << npmt 
                            << " for hit point " << hitPoint;
#endif
    if (npmt <= 0) flag = false;
  }
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HcalSim") << "HCalSD::isItinFidVolume: point " << hitPoint
                            << " return flag " << flag;
#endif
  return flag;
}

void HCalSD::getFromHFLibrary (const G4Step* aStep, bool& isKilled) {

  std::vector<HFShowerLibrary::Hit> hits = showerLibrary.get()->getHits(aStep, isKilled, weight_, false);
  if(!isKilled || hits.empty()) { return; } 

  int primaryID = setTrackID(aStep);

  // Reset entry point for new primary
  resetForNewPrimary(aStep);

  auto const theTrack = aStep->GetTrack();   
  int det = 5;

  if (G4TrackToParticleID::isGammaElectronPositron(theTrack)) {
    edepositEM  = 1.*GeV;
    edepositHAD = 0.;
  } else {
    edepositEM  = 0.;
    edepositHAD = 1.*GeV;
  }
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HcalSim") << "HCalSD::getFromLibrary " <<hits.size() 
                          << " hits for " << GetName() << " of " << primaryID 
                          << " with " << theTrack->GetDefinition()->GetParticleName() 
                          << " of " << aStep->GetPreStepPoint()->GetKineticEnergy()/GeV << " GeV";
#endif
  for (unsigned int i=0; i<hits.size(); ++i) {
    G4ThreeVector hitPoint = hits[i].position;
    if (isItinFidVolume (hitPoint)) {
      int depth              = hits[i].depth;
      double time            = hits[i].time;
      unsigned int unitID    = setDetUnitId(det, hitPoint, depth);
      currentID.setID(unitID, time, primaryID, 0);
#ifdef plotDebug
      plotProfile(aStep, hitPoint, 1.0*GeV, time, depth);
      bool emType = G4TrackToParticleID::isGammaElectronPositron(parCode);
      plotHF(hitPoint,emType);
#endif
      processHit(aStep);
    }
  }
}

void HCalSD::hitForFibre (const G4Step* aStep) { // if not ParamShower

  std::vector<HFShower::Hit> hits = hfshower.get()->getHits(aStep, weight_);
  if(hits.empty()) { return; }

  auto const theTrack = aStep->GetTrack();
  int primaryID = setTrackID(aStep);
  int det   = 5;

  if (G4TrackToParticleID::isGammaElectronPositron(theTrack)) {
    edepositEM  = 1.*GeV;
    edepositHAD = 0.;
  } else {
    edepositEM  = 0.;
    edepositHAD = 1.*GeV;
  }
 
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HcalSim") << "HCalSD::hitForFibre " << hits.size() 
                          << " hits for " << GetName() << " of " << primaryID 
                          << " with " << theTrack->GetDefinition()->GetParticleName() 
                          << " of " << preStepPoint->GetKineticEnergy()/GeV 
                          << " GeV in detector type " << det;
#endif

  for (unsigned int i=0; i<hits.size(); ++i) {
    G4ThreeVector hitPoint = hits[i].position;
    if (isItinFidVolume (hitPoint)) {
      int depth              = hits[i].depth;
      double time            = hits[i].time;
      unsigned int unitID = setDetUnitId(det, hitPoint, depth);
      currentID.setID(unitID, time, primaryID, 0);
#ifdef plotDebug
      plotProfile(aStep, hitPoint, edepositEM, time, depth);
      bool emType = (edepositEM > 0.) ? true : false;
      plotHF(hitPoint,emType);
#endif
      processHit(aStep);
    }
  }
}

void HCalSD::getFromParam (const G4Step* aStep, bool& isKilled) {

  std::vector<HFShowerParam::Hit> hits = showerParam.get()->getHits(aStep, weight_, isKilled);
  if(!isKilled || hits.empty()) { return; }

  int primaryID = setTrackID(aStep);   
  int det   = 5;

#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HcalSim") << "HCalSD::getFromParam " << hits.size() << " hits for " 
                            << GetName() << " of " << primaryID << " with " 
                            <<  aStep->GetTrack()->GetDefinition()->GetParticleName()
                            << " of " << preStepPoint->GetKineticEnergy()/GeV 
                            << " GeV in detector type " << det;
#endif
  for (unsigned int i=0; i<hits.size(); ++i) {
    G4ThreeVector hitPoint = hits[i].position;
    int depth              = hits[i].depth;
    double time            = hits[i].time;
    unsigned int unitID    = setDetUnitId(det, hitPoint, depth);
    currentID.setID(unitID, time, primaryID, 0);
    edepositEM             = hits[i].edep*GeV; 
    edepositHAD            = 0.;
#ifdef plotDebug
    plotProfile(aStep, hitPoint, edepositEM, time, depth);
#endif
    processHit(aStep);
  }
}

void HCalSD::getHitPMT (const G4Step * aStep) {

  auto const preStepPoint = aStep->GetPreStepPoint();
  auto const theTrack     = aStep->GetTrack();
  double edep  = showerPMT->getHits(aStep);

  if (edep >= 0.) {
    edep *= GeV;
    double etrack    = preStepPoint->GetKineticEnergy();
    int    primaryID = 0;
    if (etrack >= energyCut) {
      primaryID    = theTrack->GetTrackID();
    } else {
      primaryID    = theTrack->GetParentID();
      if (primaryID == 0) primaryID = theTrack->GetTrackID();
    }
    // Reset entry point for new primary
    resetForNewPrimary(aStep);
    //
    int    det      = static_cast<int>(HcalForward);
    const G4ThreeVector& hitPoint = preStepPoint->GetPosition();   
    double rr       = (hitPoint.x()*hitPoint.x() + hitPoint.y()*hitPoint.y());
    double phi      = (rr == 0. ? 0. :atan2(hitPoint.y(),hitPoint.x()));
    double etaR     = showerPMT->getRadius();
    int depth       = 1;
    if (etaR < 0) {
      depth         = 2;
      etaR          =-etaR;
    }
    if (hitPoint.z() < 0) etaR =-etaR;
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HcalSim") << "HCalSD::Hit for Detector " << det << " etaR "
                            << etaR << " phi " << phi/deg << " depth " <<depth;
#endif
    double time = (aStep->GetPostStepPoint()->GetGlobalTime());
    uint32_t unitID = 0;
    if (numberingFromDDD) {
      HcalNumberingFromDDD::HcalID tmp = numberingFromDDD.get()->unitID(det,etaR,phi,
                                                                  depth,1);
      unitID = setDetUnitId(tmp);
    }
    currentID.setID(unitID, time, primaryID, 1);

    edepositHAD = aStep->GetTotalEnergyDeposit();
    edepositEM  =-edepositHAD + edep;
#ifdef plotDebug
    plotProfile(aStep, hitPoint, edep, time, depth);
#endif
#ifdef EDM_ML_DEBUG
    double beta = preStepPoint->GetBeta();
    edm::LogVerbatim("HcalSim") << "HCalSD::getHitPMT 1 hit for " << GetName() 
                        << " of " << primaryID << " with " 
                        << theTrack->GetDefinition()->GetParticleName()
                        << " of " << preStepPoint->GetKineticEnergy()/GeV 
                        << " GeV with velocity " << beta << " UnitID "
                        << std::hex << unitID << std::dec;
#endif
    processHit(aStep);
  }
}

void HCalSD::getHitFibreBundle (const G4Step* aStep, bool type) {
  auto const preStepPoint = aStep->GetPreStepPoint();
  auto const theTrack     = aStep->GetTrack();
  double edep  = showerBundle.get()->getHits(aStep, type);

  if (edep >= 0.0) {
    edep *= GeV;
    double etrack    = preStepPoint->GetKineticEnergy();
    int    primaryID = 0;
    if (etrack >= energyCut) {
      primaryID    = theTrack->GetTrackID();
    } else {
      primaryID    = theTrack->GetParentID();
      if (primaryID == 0) primaryID = theTrack->GetTrackID();
    }
    // Reset entry point for new primary
    resetForNewPrimary(aStep);
    //
    int    det      = static_cast<int>(HcalForward);
    const G4ThreeVector& hitPoint = preStepPoint->GetPosition();   
    double rr       = hitPoint.x()*hitPoint.x() + hitPoint.y()*hitPoint.y();
    double phi      = rr == 0. ? 0. :atan2(hitPoint.y(),hitPoint.x());
    double etaR     = showerBundle->getRadius();
    int depth       = 1;
    if (etaR < 0.) {
      depth         = 2;
      etaR          =-etaR;
    }
    if (hitPoint.z() < 0.) etaR =-etaR;
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HcalSim") << "HCalSD::Hit for Detector " << det 
                << " etaR " << etaR << " phi " << phi/deg
                << " depth " << depth;
#endif
    double time = (aStep->GetPostStepPoint()->GetGlobalTime());
    uint32_t unitID = 0;
    if (numberingFromDDD) {
      HcalNumberingFromDDD::HcalID tmp = numberingFromDDD.get()->unitID(det,etaR,phi,depth,1);
      unitID = setDetUnitId(tmp);
    }
    if (type) currentID.setID(unitID, time, primaryID, 3);
    else      currentID.setID(unitID, time, primaryID, 2);

    edepositHAD = aStep->GetTotalEnergyDeposit();
    edepositEM  =-edepositHAD + edep;
#ifdef plotDebug
    plotProfile(aStep, hitPoint, edep, time, depth);
#endif
#ifdef EDM_ML_DEBUG
    double beta = preStepPoint->GetBeta();
    edm::LogVerbatim("HcalSim") << "HCalSD::getHitFibreBundle 1 hit for "
                << GetName() << " of " << primaryID << " with "
                << theTrack->GetDefinition()->GetParticleName()
                << " of " << preStepPoint->GetKineticEnergy()/GeV 
                << " GeV with velocity " << beta << " UnitID "
                << std::hex << unitID << std::dec;
#endif
    processHit(aStep);
  } // non-zero energy deposit
}

void HCalSD::readWeightFromFile(const std::string& fName) {

  std::ifstream infile;
  int entry=0;
  infile.open(fName.c_str(), std::ios::in);
  if (infile) {
    int    det, zside, etaR, phi, lay;
    double wt;
    while (infile >> det >> zside >> etaR >> phi >> lay >> wt) {
      uint32_t id = HcalTestNumbering::packHcalIndex(det,zside,1,etaR,phi,lay);
      layerWeights.insert(std::pair<uint32_t,double>(id,wt));
      ++entry;
#ifdef EDM_ML_DEBUG
      edm::LogVerbatim("HcalSim") << "HCalSD::readWeightFromFile:Entry " << entry
                              << " ID " << std::hex << id << std::dec << " ("
                              << det << "/" << zside << "/1/" << etaR << "/"
                              << phi << "/" << lay << ") Weight " << wt;
#endif
    }
    infile.close();
  }
  edm::LogVerbatim("HcalSim") << "HCalSD::readWeightFromFile: reads " << entry
                          << " weights from " << fName;
  if (entry <= 0) useLayerWt = false;
}

double HCalSD::layerWeight(int det, const G4ThreeVector& pos, int depth, int lay) { 

  double wt = 1.;
  if (numberingFromDDD) {
    //get the ID's as eta, phi, depth, ... indices
    HcalNumberingFromDDD::HcalID tmp = numberingFromDDD.get()->unitID(det, pos, 
                                                                depth, lay);
    modifyDepth(tmp);
    uint32_t id = HcalTestNumbering::packHcalIndex(tmp.subdet, tmp.zside, 1,
                                                   tmp.etaR, tmp.phis,tmp.lay);
    std::map<uint32_t,double>::const_iterator ite = layerWeights.find(id);
    if (ite != layerWeights.end()) wt = ite->second;
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HcalSim") << "HCalSD::layerWeight: ID " << std::hex << id 
                            << std::dec << " (" << tmp.subdet << "/"  
                            << tmp.zside << "/1/" << tmp.etaR << "/" 
                            << tmp.phis << "/"  << tmp.lay << ") Weight " <<wt;
#endif
  }
  return wt;
}

void HCalSD::plotProfile(const G4Step* aStep,const G4ThreeVector& global, double edep,
                         double time, int id) { 

  const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
  static const G4String modName[8] = {"HEModule", "HVQF" , "HBModule", "MBAT",
                                      "MBBT"    , "MBBTC", "MBBT_R1P", "MBBT_R1M"};
  G4ThreeVector local;
  bool found=false;
  double depth=-2000;
  int idx = 4;
  for (int n=0; n<touch->GetHistoryDepth(); ++n) {
    G4String name = touch->GetVolume(n)->GetName();
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HcalSim") << "plotProfile Depth " << n << " Name " 
                << name;
#endif
    for (unsigned int ii=0; ii<8; ++ii) {
      if (name == modName[ii]) {
        found = true;
        int dn = touch->GetHistoryDepth() - n;
        local = touch->GetHistory()->GetTransform(dn).TransformPoint(global);
        if      (ii == 0) {depth = local.z() - 4006.5; idx = 1;}
        else if (ii == 1) {depth = local.z() + 825.0;  idx = 3;}
        else if (ii == 2) {depth = local.x() - 1775.;  idx = 0;}
        else              {depth = local.y() + 15.;    idx = 2;}
        break;
      }
    }
    if (found) break;
  }
  if (!found) depth = std::abs(global.z()) - 11500;
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HcalSim") << "plotProfile Found " << found << " Global "
                  << global << " Local " << local << " depth "
                  << depth << " ID " << id << " EDEP " << edep 
                  << " Time " << time;
#endif
  if (hit_[idx]  != nullptr) hit_[idx]->Fill(edep);
  if (time_[idx] != nullptr) time_[idx]->Fill(time,edep);
  if (dist_[idx] != nullptr) dist_[idx]->Fill(depth,edep);
  int jd = 2*idx + id - 7;
  if (jd >= 0 && jd < 4) {
    jd += 5;
    if (hit_[jd]  != nullptr) hit_[jd]->Fill(edep);
    if (time_[jd] != nullptr) time_[jd]->Fill(time,edep);
    if (dist_[jd] != nullptr) dist_[jd]->Fill(depth,edep);
  }
}

void HCalSD::plotHF(const G4ThreeVector& hitPoint, bool emType) {
  double zv  = std::abs(hitPoint.z()) - gpar[4];
  if (emType) {
    if (hzvem  != nullptr) hzvem->Fill(zv);
  } else {
    if (hzvhad != nullptr) hzvhad->Fill(zv);
  }
}

void HCalSD::modifyDepth(HcalNumberingFromDDD::HcalID& id) {
  if (id.subdet == 4) {
    int ieta = (id.zside == 0) ? -id.etaR : id.etaR;
    if (hcalConstants->maxHFDepth(ieta,id.phis) > 2) {
      if (id.depth <= 2) {
        if (G4UniformRand() > 0.5) id.depth += 2;
      }
    }
  } else if ((id.subdet == 1 || id.subdet ==2) && testNumber) {
    id.depth = (depth_ == 0) ? 1 : 2;
  }
}