HFShowerLibrary Class Reference

#include <HFShowerLibrary.h>

Classes
struct	Hit

Public Member Functions
std::vector< Hit >	fillHits (const G4ThreeVector &p, const G4ThreeVector &v, int parCode, double parEnergy, bool &ok, double weight, double time, bool onlyLong=false)
std::vector< Hit >	getHits (const G4Step *aStep, bool &ok, double weight, bool onlyLong=false)
	HFShowerLibrary (const std::string &name, const HcalDDDSimConstants *hcons, const HcalSimulationParameters *hps, edm::ParameterSet const &p)
	~HFShowerLibrary ()

Protected Member Functions
void	extrapolate (int, double)
void	getRecord (int, int)
void	interpolate (int, double)
void	loadEventInfo (TBranch *)
bool	rInside (double r)
void	storePhoton (int j)

Private Attributes
bool	applyFidCut
double	backProb
double	dphi
TBranch *	emBranch
int	evtPerBin
std::unique_ptr< HFFibre >	fibre_
std::vector< double >	gpar
TBranch *	hadBranch
const HcalDDDSimConstants *	hcalConstant_
TFile *	hf
float	libVers
float	listVersion
bool	newForm
int	nMomBin
int	npe
HFShowerPhotonCollection	pe
HFShowerPhotonCollection *	photo
HFShowerPhotonCollection	photon
std::vector< double >	pmom
double	probMax
double	rMax
double	rMin
int	totEvents
bool	v3version
bool	verbose

Detailed Description

Definition at line 29 of file HFShowerLibrary.h.

Constructor & Destructor Documentation

HFShowerLibrary::HFShowerLibrary	(	const std::string &	name,
		const HcalDDDSimConstants *	hcons,
		const HcalSimulationParameters *	hps,
		edm::ParameterSet const &	p
	)

Definition at line 23 of file HFShowerLibrary.cc.

References applyFidCut, backProb, dphi, emBranch, event(), evtPerBin, Exception, fibre_, personalPlayback::fp, edm::ParameterSet::getParameter(), GeV, gpar, hadBranch, hcalConstant_, hf, mps_fire::i, info(), libVers, listVersion, loadEventInfo(), Skims_PA_cff::name, newForm, nMomBin, AlCaHLTBitMon_ParallelJobs::p, photo, pmom, probMax, rMax, rMin, contentValuesCheck::ss, AlCaHLTBitMon_QueryRunRegistry::string, totEvents, and v3version.

    : hcalConstant_(hcons), hf(nullptr), emBranch(nullptr), hadBranch(nullptr), npe(0) {
  edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFShower");
  probMax = m_HF.getParameter<double>("ProbMax");

  edm::ParameterSet m_HS = p.getParameter<edm::ParameterSet>("HFShowerLibrary");
  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");

  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";
  }

#ifdef EDM_ML_DEBUG
  std::stringstream ss;
  ss << "HFShowerLibrary: Library " << libVers << " ListVersion " << listVersion << " 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();
#endif
  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;
  }

#ifdef EDM_ML_DEBUG
  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 prob. "
                               << backProb;
#endif
  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];
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HFShower") << "HFShowerLibrary: rMIN " << rMin / CLHEP::cm << " cm and rMax " << rMax / CLHEP::cm
                               << " (Half) Phi Width of wedge " << dphi / CLHEP::deg;
#endif
  //Special Geometry parameters
  gpar = hcalConstant_->getGparHF();
}

HFShowerLibrary::~HFShowerLibrary

(

)

Definition at line 135 of file HFShowerLibrary.cc.

References hf, and photo.

                                  {
  if (hf)
    hf->Close();
  delete photo;
}

Member Function Documentation

void HFShowerLibrary::extrapolate ( int  type, double  pin  )

protected

Definition at line 511 of file HFShowerLibrary.cc.

References evtPerBin, getRecord(), createfilelist::int, dqmiolumiharvest::j, newForm, nMomBin, npe, pe, photo, photon, pmom, alignCSCRings::r, storePhoton(), totEvents, and w.

Referenced by fillHits().

                                                      {
  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 << " 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
}

std::vector< HFShowerLibrary::Hit > HFShowerLibrary::fillHits	(	const G4ThreeVector &	p,
		const G4ThreeVector &	v,
		int	parCode,
		double	parEnergy,
		bool &	ok,
		double	weight,
		double	time,
		bool	onlyLong = false
	)

Definition at line 182 of file HFShowerLibrary.cc.

References funct::abs(), applyFidCut, backProb, funct::cos(), HFShowerLibrary::Hit::depth, LEDCalibrationChannels::depth, dphi, JetChargeProducer_cfi::exp, extrapolate(), fibre_, gpar, mps_fire::i, createfilelist::int, interpolate(), G4TrackToParticleID::isGammaElectronPositron(), G4TrackToParticleID::isStableHadron(), MeV, nMomBin, npe, AlCaHLTBitMon_ParallelJobs::p, pe, pmom, HFShowerLibrary::Hit::position, probMax, alignCSCRings::r, diffTwoXMLs::r1, diffTwoXMLs::r2, rInside(), funct::sin(), MessageLogger_cff::threshold, HFShowerLibrary::Hit::time, geometryCSVtoXML::xx, geometryCSVtoXML::yy, z, and geometryCSVtoXML::zz.

Referenced by getHits().

                                                                           {
  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])) {
        oneHit.position = pos;
        oneHit.depth = depth;
        oneHit.time = (tSlice + (pe[i].t()) + (fibre_->tShift(lpos, depth, 1)));
        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() << ":"
                                     << fibre_->tShift(lpos, depth, 1) << ":" << (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]) {
          oneHit.position = pos;
          oneHit.depth = 2;
          oneHit.time = (tSlice + (pe[i].t()) + (fibre_->tShift(lpos, 2, 1)));
          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()
                                       << ":" << fibre_->tShift(lpos, 2, 1) << ":" << (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;
}

std::vector< HFShowerLibrary::Hit > HFShowerLibrary::getHits	(	const G4Step *	aStep,
		bool &	ok,
		double	weight,
		bool	onlyLong = false
	)

Definition at line 141 of file HFShowerLibrary.cc.

References funct::cos(), fillHits(), GeV, gpar, funct::sin(), and HLT_2018_cff::track.

                                                                          {
  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() / 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);
}

void HFShowerLibrary::getRecord ( int  type, int  record  )

protected

Definition at line 344 of file HFShowerLibrary.cc.

References emBranch, hadBranch, mps_fire::i, dqmiolumiharvest::j, newForm, photo, photon, OrderedSet::t, totEvents, and v3version.

Referenced by extrapolate(), and interpolate().

                                                    {
  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::interpolate ( int  type, double  pin  )

protected

Definition at line 434 of file HFShowerLibrary.cc.

References evtPerBin, getRecord(), GeV, createfilelist::int, dqmiolumiharvest::j, newForm, nMomBin, npe, pe, photo, photon, pmom, alignCSCRings::r, storePhoton(), totEvents, and w.

Referenced by fillHits().

                                                      {
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HFShower") << "HFShowerLibrary:: Interpolate for Energy " << pin / 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::loadEventInfo ( TBranch *  branch )

protected

Definition at line 403 of file HFShowerLibrary.cc.

References evtPerBin, GeV, mps_fire::i, libVers, listVersion, nMomBin, pmom, and totEvents.

Referenced by HFShowerLibrary().

                                                   {
  if (branch) {
    std::vector<HFShowerLibraryEventInfo> eventInfoCollection;
    branch->SetAddress(&eventInfoCollection);
    branch->GetEntry(0);
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HFShower") << "HFShowerLibrary::loadEventInfo loads EventInfo Collection of size "
                                 << eventInfoCollection.size() << " records";
#endif
    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 {
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HFShower") << "HFShowerLibrary::loadEventInfo loads EventInfo from hardwired"
                                 << " numbers";
#endif
    nMomBin = 16;
    evtPerBin = 5000;
    totEvents = nMomBin * evtPerBin;
    libVers = 1.1;
    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] *= GeV;
}

bool HFShowerLibrary::rInside ( double  r )

protected

Definition at line 342 of file HFShowerLibrary.cc.

References rMax, and rMin.

Referenced by fillHits().

{ return (r >= rMin && r <= rMax); }

void HFShowerLibrary::storePhoton ( int  j )

protected

Definition at line 576 of file HFShowerLibrary.cc.

References newForm, npe, pe, photo, and photon.

Referenced by extrapolate(), and interpolate().

                                       {
  if (newForm)
    pe.push_back(photo->at(j));
  else
    pe.push_back(photon[j]);
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HFShower") << "HFShowerLibrary: storePhoton " << j << " npe " << npe << " " << pe[npe];
#endif
  npe++;
}

Member Data Documentation

bool HFShowerLibrary::applyFidCut

private

Definition at line 70 of file HFShowerLibrary.h.

Referenced by fillHits(), and HFShowerLibrary().

double HFShowerLibrary::backProb

private

Definition at line 75 of file HFShowerLibrary.h.

Referenced by fillHits(), and HFShowerLibrary().

double HFShowerLibrary::dphi

private

Definition at line 76 of file HFShowerLibrary.h.

Referenced by fillHits(), and HFShowerLibrary().

TBranch* HFShowerLibrary::emBranch

private

Definition at line 68 of file HFShowerLibrary.h.

Referenced by getRecord(), and HFShowerLibrary().

int HFShowerLibrary::evtPerBin

private

Definition at line 71 of file HFShowerLibrary.h.

Referenced by extrapolate(), HFShowerLibrary(), interpolate(), and loadEventInfo().

std::unique_ptr<HFFibre> HFShowerLibrary::fibre_

private

Definition at line 66 of file HFShowerLibrary.h.

Referenced by fillHits(), and HFShowerLibrary().

std::vector<double> HFShowerLibrary::gpar

private

Definition at line 77 of file HFShowerLibrary.h.

Referenced by fillHits(), getHits(), and HFShowerLibrary().

TBranch * HFShowerLibrary::hadBranch

private

Definition at line 68 of file HFShowerLibrary.h.

Referenced by getRecord(), and HFShowerLibrary().

const HcalDDDSimConstants* HFShowerLibrary::hcalConstant_

private

Definition at line 65 of file HFShowerLibrary.h.

Referenced by HFShowerLibrary().

TFile* HFShowerLibrary::hf

private

Definition at line 67 of file HFShowerLibrary.h.

Referenced by HFShowerLibrary(), and ~HFShowerLibrary().

float HFShowerLibrary::libVers

private

Definition at line 72 of file HFShowerLibrary.h.

Referenced by HFShowerLibrary(), and loadEventInfo().

float HFShowerLibrary::listVersion

private

Definition at line 72 of file HFShowerLibrary.h.

Referenced by HFShowerLibrary(), and loadEventInfo().

bool HFShowerLibrary::newForm

private

Definition at line 70 of file HFShowerLibrary.h.

Referenced by extrapolate(), getRecord(), HFShowerLibrary(), interpolate(), and storePhoton().

int HFShowerLibrary::nMomBin

private

Definition at line 71 of file HFShowerLibrary.h.

Referenced by extrapolate(), fillHits(), HFShowerLibrary(), interpolate(), and loadEventInfo().

int HFShowerLibrary::npe

private

Definition at line 79 of file HFShowerLibrary.h.

Referenced by extrapolate(), fillHits(), interpolate(), and storePhoton().

HFShowerPhotonCollection HFShowerLibrary::pe

private

Definition at line 80 of file HFShowerLibrary.h.

Referenced by extrapolate(), fillHits(), interpolate(), and storePhoton().

HFShowerPhotonCollection* HFShowerLibrary::photo

private

Definition at line 81 of file HFShowerLibrary.h.

Referenced by extrapolate(), getRecord(), HFShowerLibrary(), interpolate(), storePhoton(), and ~HFShowerLibrary().

HFShowerPhotonCollection HFShowerLibrary::photon

private

Definition at line 82 of file HFShowerLibrary.h.

Referenced by extrapolate(), getRecord(), interpolate(), and storePhoton().

std::vector<double> HFShowerLibrary::pmom

private

Definition at line 73 of file HFShowerLibrary.h.

Referenced by extrapolate(), fillHits(), HFShowerLibrary(), interpolate(), and loadEventInfo().

double HFShowerLibrary::probMax

private

Definition at line 75 of file HFShowerLibrary.h.

Referenced by fillHits(), and HFShowerLibrary().

double HFShowerLibrary::rMax

private

Definition at line 76 of file HFShowerLibrary.h.

Referenced by HFShowerLibrary(), and rInside().

double HFShowerLibrary::rMin

private

Definition at line 76 of file HFShowerLibrary.h.

Referenced by HFShowerLibrary(), and rInside().

int HFShowerLibrary::totEvents

private

Definition at line 71 of file HFShowerLibrary.h.

Referenced by extrapolate(), getRecord(), HFShowerLibrary(), interpolate(), and loadEventInfo().

bool HFShowerLibrary::v3version

private

Definition at line 70 of file HFShowerLibrary.h.

Referenced by getRecord(), and HFShowerLibrary().

bool HFShowerLibrary::verbose

private

Definition at line 70 of file HFShowerLibrary.h.