HFShowerLibrary.cc

Go to the documentation of this file.

// File: HFShowerLibrary.cc
// Description: Shower library for Very forward hadron calorimeter
 
#include "SimG4CMS/Calo/interface/HFShowerLibrary.h"
#include "SimDataFormats/CaloHit/interface/HFShowerLibraryEventInfo.h"
#include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
 
#include "FWCore/Utilities/interface/Exception.h"
 
#include "G4VPhysicalVolume.hh"
#include "G4NavigationHistory.hh"
#include "G4Step.hh"
#include "G4Track.hh"
#include "G4ParticleTable.hh"
#include "Randomize.hh"
#include "CLHEP/Units/SystemOfUnits.h"
#include "CLHEP/Units/PhysicalConstants.h"
 
//#define EDM_ML_DEBUG
 
HFShowerLibrary::HFShowerLibrary(const std::string& name,
                                 const HcalDDDSimConstants* hcons,
                                 const HcalSimulationParameters* hps,
                                 edm::ParameterSet const& p)
    : hcalConstant_(hcons), hf(nullptr), emBranch(nullptr), hadBranch(nullptr), npe(0) {
  edm::ParameterSet m_HF =
      (p.getParameter<edm::ParameterSet>("HFShower")).getParameter<edm::ParameterSet>("HFShowerBlock");
  probMax = m_HF.getParameter<double>("ProbMax");
  equalizeTimeShift_ = m_HF.getParameter<bool>("EqualizeTimeShift");
 
  edm::ParameterSet m_HS =
      (p.getParameter<edm::ParameterSet>("HFShowerLibrary")).getParameter<edm::ParameterSet>("HFLibraryFileBlock");
  edm::FileInPath fp = m_HS.getParameter<edm::FileInPath>("FileName");
  std::string pTreeName = fp.fullPath();
  backProb = m_HS.getParameter<double>("BackProbability");
  std::string emName = m_HS.getParameter<std::string>("TreeEMID");
  std::string hadName = m_HS.getParameter<std::string>("TreeHadID");
  std::string branchEvInfo = m_HS.getUntrackedParameter<std::string>(
      "BranchEvt", "HFShowerLibraryEventInfos_hfshowerlib_HFShowerLibraryEventInfo");
  std::string branchPre = m_HS.getUntrackedParameter<std::string>("BranchPre", "HFShowerPhotons_hfshowerlib_");
  std::string branchPost = m_HS.getUntrackedParameter<std::string>("BranchPost", "_R.obj");
  verbose = m_HS.getUntrackedParameter<bool>("Verbosity", false);
  applyFidCut = m_HS.getParameter<bool>("ApplyFiducialCut");
  fileVersion_ = m_HS.getParameter<int>("FileVersion");
 
  if (pTreeName.find('.') == 0)
    pTreeName.erase(0, 2);
  const char* nTree = pTreeName.c_str();
  hf = TFile::Open(nTree);
 
  if (!hf->IsOpen()) {
    edm::LogError("HFShower") << "HFShowerLibrary: opening " << nTree << " failed";
    throw cms::Exception("Unknown", "HFShowerLibrary") << "Opening of " << pTreeName << " fails\n";
  } else {
    edm::LogVerbatim("HFShower") << "HFShowerLibrary: opening " << nTree << " successfully";
  }
 
  newForm = (branchEvInfo.empty());
  TTree* event(nullptr);
  if (newForm)
    event = (TTree*)hf->Get("HFSimHits");
  else
    event = (TTree*)hf->Get("Events");
  if (event) {
    TBranch* evtInfo(nullptr);
    if (!newForm) {
      std::string info = branchEvInfo + branchPost;
      evtInfo = event->GetBranch(info.c_str());
    }
    if (evtInfo || newForm) {
      loadEventInfo(evtInfo);
    } else {
      edm::LogError("HFShower") << "HFShowerLibrary: HFShowerLibrayEventInfo"
                                << " Branch does not exist in Event";
      throw cms::Exception("Unknown", "HFShowerLibrary") << "Event information absent\n";
    }
  } else {
    edm::LogError("HFShower") << "HFShowerLibrary: Events Tree does not "
                              << "exist";
    throw cms::Exception("Unknown", "HFShowerLibrary") << "Events tree absent\n";
  }
 
  std::stringstream ss;
  ss << "HFShowerLibrary: Library " << libVers << " ListVersion " << listVersion << " File version " << fileVersion_
     << " Events Total " << totEvents << " and " << evtPerBin << " per bin\n";
  ss << "HFShowerLibrary: Energies (GeV) with " << nMomBin << " bins\n";
  for (int i = 0; i < nMomBin; ++i) {
    if (i / 10 * 10 == i && i > 0) {
      ss << "\n";
    }
    ss << "  " << pmom[i] / CLHEP::GeV;
  }
  edm::LogVerbatim("HFShower") << ss.str();
 
  std::string nameBr = branchPre + emName + branchPost;
  emBranch = event->GetBranch(nameBr.c_str());
  if (verbose)
    emBranch->Print();
  nameBr = branchPre + hadName + branchPost;
  hadBranch = event->GetBranch(nameBr.c_str());
  if (verbose)
    hadBranch->Print();
 
  v3version = false;
  if (emBranch->GetClassName() == std::string("vector<float>")) {
    v3version = true;
  }
 
  edm::LogVerbatim("HFShower") << " HFShowerLibrary:Branch " << emName << " has " << emBranch->GetEntries()
                               << " entries and Branch " << hadName << " has " << hadBranch->GetEntries()
                               << " entries\n HFShowerLibrary::No packing information - Assume x, y, z are not in "
                                  "packed form\n Maximum probability cut off "
                               << probMax << "  Back propagation of light probability " << backProb
                               << " Flag for equalizing Time Shift for different eta " << equalizeTimeShift_;
 
  fibre_ = std::make_unique<HFFibre>(name, hcalConstant_, hps, p);
  photo = new HFShowerPhotonCollection;
 
  //Radius (minimum and maximum)
  std::vector<double> rTable = hcalConstant_->getRTableHF();
  rMin = rTable[0];
  rMax = rTable[rTable.size() - 1];
 
  //Delta phi
  std::vector<double> phibin = hcalConstant_->getPhiTableHF();
  dphi = phibin[0];
 
  edm::LogVerbatim("HFShower") << "HFShowerLibrary: rMIN " << rMin / CLHEP::cm << " cm and rMax " << rMax / CLHEP::cm
                               << " (Half) Phi Width of wedge " << dphi / CLHEP::deg;
 
  //Special Geometry parameters
  gpar = hcalConstant_->getGparHF();
}
 
HFShowerLibrary::~HFShowerLibrary() {
  if (hf)
    hf->Close();
  delete photo;
}
 
std::vector<HFShowerLibrary::Hit> HFShowerLibrary::getHits(const G4Step* aStep,
                                                           bool& isKilled,
                                                           double weight,
                                                           bool onlyLong) {
  auto const preStepPoint = aStep->GetPreStepPoint();
  auto const postStepPoint = aStep->GetPostStepPoint();
  auto const track = aStep->GetTrack();
  // Get Z-direction
  auto const aParticle = track->GetDynamicParticle();
  const G4ThreeVector& momDir = aParticle->GetMomentumDirection();
  const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
  int parCode = track->GetDefinition()->GetPDGEncoding();
 
  // VI: for ions use internally pdg code of alpha in order to keep
  // consistency with previous simulation
  if (track->GetDefinition()->IsGeneralIon()) {
    parCode = 1000020040;
  }
 
#ifdef EDM_ML_DEBUG
  G4String partType = track->GetDefinition()->GetParticleName();
  const G4ThreeVector localPos = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
  double zoff = localPos.z() + 0.5 * gpar[1];
 
  edm::LogVerbatim("HFShower") << "HFShowerLibrary::getHits " << partType << " of energy "
                               << track->GetKineticEnergy() / CLHEP::GeV << " GeV weight= " << weight
                               << " onlyLong: " << onlyLong << "  dir.orts " << momDir.x() << ", " << momDir.y() << ", "
                               << momDir.z() << "  Pos x,y,z = " << hitPoint.x() << "," << hitPoint.y() << ","
                               << hitPoint.z() << " (" << zoff << ")   sphi,cphi,stheta,ctheta  = " << sin(momDir.phi())
                               << "," << cos(momDir.phi()) << ", " << sin(momDir.theta()) << "," << cos(momDir.theta());
#endif
 
  double tSlice = (postStepPoint->GetGlobalTime()) / CLHEP::nanosecond;
 
  // use kinetic energy for protons and ions
  double pin = (track->GetDefinition()->GetBaryonNumber() > 0) ? preStepPoint->GetKineticEnergy()
                                                               : preStepPoint->GetTotalEnergy();
 
  return fillHits(hitPoint, momDir, parCode, pin, isKilled, weight, tSlice, onlyLong);
}
 
std::vector<HFShowerLibrary::Hit> HFShowerLibrary::fillHits(const G4ThreeVector& hitPoint,
                                                            const G4ThreeVector& momDir,
                                                            int parCode,
                                                            double pin,
                                                            bool& ok,
                                                            double weight,
                                                            double tSlice,
                                                            bool onlyLong) {
  std::vector<HFShowerLibrary::Hit> hit;
  ok = false;
  bool isEM = G4TrackToParticleID::isGammaElectronPositron(parCode);
  // shower is built only for gamma, e+- and stable hadrons
  if (!isEM && !G4TrackToParticleID::isStableHadron(parCode)) {
    return hit;
  }
  ok = true;
 
  // remove low-energy component
  const double threshold = 50 * MeV;
  if (pin < threshold) {
    return hit;
  }
 
  double pz = momDir.z();
  double zint = hitPoint.z();
 
  // if particle moves from interaction point or "backwards (halo)
  bool backward = (pz * zint < 0.) ? true : false;
 
  double sphi = sin(momDir.phi());
  double cphi = cos(momDir.phi());
  double ctheta = cos(momDir.theta());
  double stheta = sin(momDir.theta());
 
  if (isEM) {
    if (pin < pmom[nMomBin - 1]) {
      interpolate(0, pin);
    } else {
      extrapolate(0, pin);
    }
  } else {
    if (pin < pmom[nMomBin - 1]) {
      interpolate(1, pin);
    } else {
      extrapolate(1, pin);
    }
  }
 
  int nHit = 0;
  HFShowerLibrary::Hit oneHit;
  for (int i = 0; i < npe; ++i) {
    double zv = std::abs(pe[i].z());  // abs local z
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HFShower") << "HFShowerLibrary: Hit " << i << " " << pe[i] << " zv " << zv;
#endif
    if (zv <= gpar[1] && pe[i].lambda() > 0 && (pe[i].z() >= 0 || (zv > gpar[0] && (!onlyLong)))) {
      int depth = 1;
      if (onlyLong) {
      } else if (!backward) {  // fully valid only for "front" particles
        if (pe[i].z() < 0)
          depth = 2;                 // with "front"-simulated shower lib.
      } else {                       // for "backward" particles - almost equal
        double r = G4UniformRand();  // share between L and S fibers
        if (r > 0.5)
          depth = 2;
      }
 
      // Updated coordinate transformation from local
      //  back to global using two Euler angles: phi and theta
      double pex = pe[i].x();
      double pey = pe[i].y();
 
      double xx = pex * ctheta * cphi - pey * sphi + zv * stheta * cphi;
      double yy = pex * ctheta * sphi + pey * cphi + zv * stheta * sphi;
      double zz = -pex * stheta + zv * ctheta;
 
      G4ThreeVector pos = hitPoint + G4ThreeVector(xx, yy, zz);
      zv = std::abs(pos.z()) - gpar[4] - 0.5 * gpar[1];
      G4ThreeVector lpos = G4ThreeVector(pos.x(), pos.y(), zv);
 
      zv = fibre_->zShift(lpos, depth, 0);  // distance to PMT !
 
      double r = pos.perp();
      double p = fibre_->attLength(pe[i].lambda());
      double fi = pos.phi();
      if (fi < 0)
        fi += CLHEP::twopi;
      int isect = int(fi / dphi) + 1;
      isect = (isect + 1) / 2;
      double dfi = ((isect * 2 - 1) * dphi - fi);
      if (dfi < 0)
        dfi = -dfi;
      double dfir = r * sin(dfi);
#ifdef EDM_ML_DEBUG
      edm::LogVerbatim("HFShower") << "HFShowerLibrary: Position shift " << xx << ", " << yy << ", " << zz << ": "
                                   << pos << " R " << r << " Phi " << fi << " Section " << isect << " R*Dfi " << dfir
                                   << " Dist " << zv;
#endif
      zz = std::abs(pos.z());
      double r1 = G4UniformRand();
      double r2 = G4UniformRand();
      double r3 = backward ? G4UniformRand() : -9999.;
      if (!applyFidCut)
        dfir += gpar[5];
 
#ifdef EDM_ML_DEBUG
      edm::LogVerbatim("HFShower") << "HFShowerLibrary: rLimits " << rInside(r) << " attenuation " << r1 << ":"
                                   << exp(-p * zv) << " r2 " << r2 << " r3 " << r3 << " rDfi " << gpar[5] << " zz "
                                   << zz << " zLim " << gpar[4] << ":" << gpar[4] + gpar[1] << "\n"
                                   << "  rInside(r) :" << rInside(r) << "  r1 <= exp(-p*zv) :" << (r1 <= exp(-p * zv))
                                   << "  r2 <= probMax :" << (r2 <= probMax * weight)
                                   << "  r3 <= backProb :" << (r3 <= backProb)
                                   << "  dfir > gpar[5] :" << (dfir > gpar[5]) << "  zz >= gpar[4] :" << (zz >= gpar[4])
                                   << "  zz <= gpar[4]+gpar[1] :" << (zz <= gpar[4] + gpar[1]);
#endif
      if (rInside(r) && r1 <= exp(-p * zv) && r2 <= probMax * weight && dfir > gpar[5] && zz >= gpar[4] &&
          zz <= gpar[4] + gpar[1] && r3 <= backProb && (depth != 2 || zz >= gpar[4] + gpar[0])) {
        double tdiff = (equalizeTimeShift_) ? (fibre_->tShift(lpos, depth, -1)) : (fibre_->tShift(lpos, depth, 1));
        oneHit.position = pos;
        oneHit.depth = depth;
        oneHit.time = (tSlice + (pe[i].t()) + tdiff);
        hit.push_back(oneHit);
 
#ifdef EDM_ML_DEBUG
        edm::LogVerbatim("HFShower") << "HFShowerLibrary: Final Hit " << nHit << " position " << (hit[nHit].position)
                                     << " Depth " << (hit[nHit].depth) << " Time " << tSlice << ":" << pe[i].t() << ":"
                                     << tdiff << ":" << (hit[nHit].time);
#endif
        ++nHit;
      }
#ifdef EDM_ML_DEBUG
      else
        edm::LogVerbatim("HFShower") << "HFShowerLibrary: REJECTED !!!";
#endif
      if (onlyLong && zz >= gpar[4] + gpar[0] && zz <= gpar[4] + gpar[1]) {
        r1 = G4UniformRand();
        r2 = G4UniformRand();
        if (rInside(r) && r1 <= exp(-p * zv) && r2 <= probMax && dfir > gpar[5]) {
          double tdiff = (equalizeTimeShift_) ? (fibre_->tShift(lpos, 2, -1)) : (fibre_->tShift(lpos, 2, 1));
          oneHit.position = pos;
          oneHit.depth = 2;
          oneHit.time = (tSlice + (pe[i].t()) + tdiff);
          hit.push_back(oneHit);
#ifdef EDM_ML_DEBUG
          edm::LogVerbatim("HFShower") << "HFShowerLibrary: Final Hit " << nHit << " position " << (hit[nHit].position)
                                       << " Depth " << (hit[nHit].depth) << " Time " << tSlice << ":" << pe[i].t()
                                       << ":" << tdiff << ":" << (hit[nHit].time);
#endif
          ++nHit;
        }
      }
    }
  }
 
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HFShower") << "HFShowerLibrary: Total Hits " << nHit << " out of " << npe << " PE";
#endif
  if (nHit > npe && !onlyLong) {
    edm::LogWarning("HFShower") << "HFShowerLibrary: Hit buffer " << npe << " smaller than " << nHit << " Hits";
  }
  return hit;
}
 
bool HFShowerLibrary::rInside(double r) { return (r >= rMin && r <= rMax); }
 
void HFShowerLibrary::getRecord(int type, int record) {
  int nrc = record - 1;
  photon.clear();
  photo->clear();
  if (type > 0) {
    if (newForm) {
      if (!v3version) {
        hadBranch->SetAddress(&photo);
        hadBranch->GetEntry(nrc + totEvents);
      } else {
        std::vector<float> t;
        std::vector<float>* tp = &t;
        hadBranch->SetAddress(&tp);
        hadBranch->GetEntry(nrc + totEvents);
        unsigned int tSize = t.size() / 5;
        photo->reserve(tSize);
        for (unsigned int i = 0; i < tSize; i++) {
          photo->push_back(
              HFShowerPhoton(t[i], t[1 * tSize + i], t[2 * tSize + i], t[3 * tSize + i], t[4 * tSize + i]));
        }
      }
    } else {
      hadBranch->SetAddress(&photon);
      hadBranch->GetEntry(nrc);
    }
  } else {
    if (newForm) {
      if (!v3version) {
        emBranch->SetAddress(&photo);
        emBranch->GetEntry(nrc);
      } else {
        std::vector<float> t;
        std::vector<float>* tp = &t;
        emBranch->SetAddress(&tp);
        emBranch->GetEntry(nrc);
        unsigned int tSize = t.size() / 5;
        photo->reserve(tSize);
        for (unsigned int i = 0; i < tSize; i++) {
          photo->push_back(
              HFShowerPhoton(t[i], t[1 * tSize + i], t[2 * tSize + i], t[3 * tSize + i], t[4 * tSize + i]));
        }
      }
    } else {
      emBranch->SetAddress(&photon);
      emBranch->GetEntry(nrc);
    }
  }
#ifdef EDM_ML_DEBUG
  int nPhoton = (newForm) ? photo->size() : photon.size();
  edm::LogVerbatim("HFShower") << "HFShowerLibrary::getRecord: Record " << record << " of type " << type << " with "
                               << nPhoton << " photons";
  for (int j = 0; j < nPhoton; j++)
    if (newForm)
      edm::LogVerbatim("HFShower") << "Photon " << j << " " << photo->at(j);
    else
      edm::LogVerbatim("HFShower") << "Photon " << j << " " << photon[j];
#endif
}
 
void HFShowerLibrary::loadEventInfo(TBranch* branch) {
  if (branch) {
    std::vector<HFShowerLibraryEventInfo> eventInfoCollection;
    branch->SetAddress(&eventInfoCollection);
    branch->GetEntry(0);
    edm::LogVerbatim("HFShower") << "HFShowerLibrary::loadEventInfo loads EventInfo Collection of size "
                                 << eventInfoCollection.size() << " records";
 
    totEvents = eventInfoCollection[0].totalEvents();
    nMomBin = eventInfoCollection[0].numberOfBins();
    evtPerBin = eventInfoCollection[0].eventsPerBin();
    libVers = eventInfoCollection[0].showerLibraryVersion();
    listVersion = eventInfoCollection[0].physListVersion();
    pmom = eventInfoCollection[0].energyBins();
  } else {
    edm::LogVerbatim("HFShower") << "HFShowerLibrary::loadEventInfo loads EventInfo from hardwired"
                                 << " numbers";
 
    nMomBin = 16;
    evtPerBin = (fileVersion_ == 0) ? 5000 : 10000;
    totEvents = nMomBin * evtPerBin;
    libVers = (fileVersion_ == 0) ? 1.1 : 1.2;
    listVersion = 3.6;
    pmom = {2, 3, 5, 7, 10, 15, 20, 30, 50, 75, 100, 150, 250, 350, 500, 1000};
  }
  for (int i = 0; i < nMomBin; i++)
    pmom[i] *= CLHEP::GeV;
}
 
void HFShowerLibrary::interpolate(int type, double pin) {
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HFShower") << "HFShowerLibrary:: Interpolate for Energy " << pin / CLHEP::GeV << " GeV with "
                               << nMomBin << " momentum bins and " << evtPerBin << " entries/bin -- total "
                               << totEvents;
#endif
  int irc[2] = {0, 0};
  double w = 0.;
  double r = G4UniformRand();
 
  if (pin < pmom[0]) {
    w = pin / pmom[0];
    irc[1] = int(evtPerBin * r) + 1;
    irc[0] = 0;
  } else {
    for (int j = 0; j < nMomBin - 1; j++) {
      if (pin >= pmom[j] && pin < pmom[j + 1]) {
        w = (pin - pmom[j]) / (pmom[j + 1] - pmom[j]);
        if (j == nMomBin - 2) {
          irc[1] = int(evtPerBin * 0.5 * r);
        } else {
          irc[1] = int(evtPerBin * r);
        }
        irc[1] += (j + 1) * evtPerBin + 1;
        r = G4UniformRand();
        irc[0] = int(evtPerBin * r) + 1 + j * evtPerBin;
        if (irc[0] < 0) {
          edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[0] = " << irc[0] << " now set to 0";
          irc[0] = 0;
        } else if (irc[0] > totEvents) {
          edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[0] = " << irc[0] << " now set to " << totEvents;
          irc[0] = totEvents;
        }
      }
    }
  }
  if (irc[1] < 1) {
    edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[1] = " << irc[1] << " now set to 1";
    irc[1] = 1;
  } else if (irc[1] > totEvents) {
    edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[1] = " << irc[1] << " now set to " << totEvents;
    irc[1] = totEvents;
  }
 
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HFShower") << "HFShowerLibrary:: Select records " << irc[0] << " and " << irc[1] << " with weights "
                               << 1 - w << " and " << w;
#endif
  pe.clear();
  npe = 0;
  int npold = 0;
  for (int ir = 0; ir < 2; ir++) {
    if (irc[ir] > 0) {
      getRecord(type, irc[ir]);
      int nPhoton = (newForm) ? photo->size() : photon.size();
      npold += nPhoton;
      for (int j = 0; j < nPhoton; j++) {
        r = G4UniformRand();
        if ((ir == 0 && r > w) || (ir > 0 && r < w)) {
          storePhoton(j);
        }
      }
    }
  }
 
  if ((npe > npold || (npold == 0 && irc[0] > 0)) && !(npe == 0 && npold == 0))
    edm::LogWarning("HFShower") << "HFShowerLibrary: Interpolation Warning =="
                                << " records " << irc[0] << " and " << irc[1] << " gives a buffer of " << npold
                                << " photons and fills " << npe << " *****";
#ifdef EDM_ML_DEBUG
  else
    edm::LogVerbatim("HFShower") << "HFShowerLibrary: Interpolation == records " << irc[0] << " and " << irc[1]
                                 << " gives a buffer of " << npold << " photons and fills " << npe << " PE";
  for (int j = 0; j < npe; j++)
    edm::LogVerbatim("HFShower") << "Photon " << j << " " << pe[j];
#endif
}
 
void HFShowerLibrary::extrapolate(int type, double pin) {
  int nrec = int(pin / pmom[nMomBin - 1]);
  double w = (pin - pmom[nMomBin - 1] * nrec) / pmom[nMomBin - 1];
  nrec++;
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HFShower") << "HFShowerLibrary:: Extrapolate for Energy " << pin / CLHEP::GeV << " GeV with "
                               << nMomBin << " momentum bins and " << evtPerBin << " entries/bin -- "
                               << "total " << totEvents << " using " << nrec << " records";
#endif
  std::vector<int> irc(nrec);
 
  for (int ir = 0; ir < nrec; ir++) {
    double r = G4UniformRand();
    irc[ir] = int(evtPerBin * 0.5 * r) + (nMomBin - 1) * evtPerBin + 1;
    if (irc[ir] < 1) {
      edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[" << ir << "] = " << irc[ir] << " now set to 1";
      irc[ir] = 1;
    } else if (irc[ir] > totEvents) {
      edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[" << ir << "] = " << irc[ir] << " now set to "
                                  << totEvents;
      irc[ir] = totEvents;
#ifdef EDM_ML_DEBUG
    } else {
      edm::LogVerbatim("HFShower") << "HFShowerLibrary::Extrapolation use irc[" << ir << "] = " << irc[ir];
#endif
    }
  }
 
  pe.clear();
  npe = 0;
  int npold = 0;
  for (int ir = 0; ir < nrec; ir++) {
    if (irc[ir] > 0) {
      getRecord(type, irc[ir]);
      int nPhoton = (newForm) ? photo->size() : photon.size();
      npold += nPhoton;
      for (int j = 0; j < nPhoton; j++) {
        double r = G4UniformRand();
        if (ir != nrec - 1 || r < w) {
          storePhoton(j);
        }
      }
#ifdef EDM_ML_DEBUG
      edm::LogVerbatim("HFShower") << "HFShowerLibrary: Record [" << ir << "] = " << irc[ir] << " npold = " << npold;
#endif
    }
  }
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HFShower") << "HFShowerLibrary:: uses " << npold << " photons";
#endif
 
  if (npe > npold || npold == 0)
    edm::LogWarning("HFShower") << "HFShowerLibrary: Extrapolation Warning == " << nrec << " records " << irc[0] << ", "
                                << irc[1] << ", ... gives a buffer of " << npold << " photons and fills " << npe
                                << " *****";
#ifdef EDM_ML_DEBUG
  else
    edm::LogVerbatim("HFShower") << "HFShowerLibrary: Extrapolation == " << nrec << " records " << irc[0] << ", "
                                 << irc[1] << ", ... gives a buffer of " << npold << " photons and fills " << npe
                                 << " PE";
  for (int j = 0; j < npe; j++)
    edm::LogVerbatim("HFShower") << "Photon " << j << " " << pe[j];
#endif
}
 
void HFShowerLibrary::storePhoton(int j) {
  if (newForm)
    pe.push_back(photo->at(j));
  else
    pe.push_back(photon[j]);
  npe++;
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HFShower") << "HFShowerLibrary: storePhoton " << j << " npe " << npe << " " << pe[npe - 1];
#endif
}