Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-06 04:01:07

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: HFShowerLibrary.cc
0003 // Description: Shower library for Very forward hadron calorimeter
0004 ///////////////////////////////////////////////////////////////////////////////
0005 
0006 #include "SimG4CMS/Calo/interface/HFShowerLibrary.h"
0007 #include "SimDataFormats/CaloHit/interface/HFShowerLibraryEventInfo.h"
0008 #include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
0009 
0010 #include "FWCore/Utilities/interface/Exception.h"
0011 
0012 #include "G4VPhysicalVolume.hh"
0013 #include "G4NavigationHistory.hh"
0014 #include "G4Step.hh"
0015 #include "G4Track.hh"
0016 #include "G4ParticleTable.hh"
0017 #include "Randomize.hh"
0018 #include "CLHEP/Units/SystemOfUnits.h"
0019 #include "CLHEP/Units/PhysicalConstants.h"
0020 
0021 //#define EDM_ML_DEBUG
0022 
0023 HFShowerLibrary::HFShowerLibrary(const std::string& name,
0024                                  const HcalDDDSimConstants* hcons,
0025                                  const HcalSimulationParameters* hps,
0026                                  edm::ParameterSet const& p)
0027     : hcalConstant_(hcons), hf(nullptr), emBranch(nullptr), hadBranch(nullptr), npe(0) {
0028   edm::ParameterSet m_HF =
0029       (p.getParameter<edm::ParameterSet>("HFShower")).getParameter<edm::ParameterSet>("HFShowerBlock");
0030   probMax = m_HF.getParameter<double>("ProbMax");
0031   equalizeTimeShift_ = m_HF.getParameter<bool>("EqualizeTimeShift");
0032 
0033   edm::ParameterSet m_HS =
0034       (p.getParameter<edm::ParameterSet>("HFShowerLibrary")).getParameter<edm::ParameterSet>("HFLibraryFileBlock");
0035   edm::FileInPath fp = m_HS.getParameter<edm::FileInPath>("FileName");
0036   std::string pTreeName = fp.fullPath();
0037   backProb = m_HS.getParameter<double>("BackProbability");
0038   std::string emName = m_HS.getParameter<std::string>("TreeEMID");
0039   std::string hadName = m_HS.getParameter<std::string>("TreeHadID");
0040   std::string branchEvInfo = m_HS.getUntrackedParameter<std::string>(
0041       "BranchEvt", "HFShowerLibraryEventInfos_hfshowerlib_HFShowerLibraryEventInfo");
0042   std::string branchPre = m_HS.getUntrackedParameter<std::string>("BranchPre", "HFShowerPhotons_hfshowerlib_");
0043   std::string branchPost = m_HS.getUntrackedParameter<std::string>("BranchPost", "_R.obj");
0044   verbose = m_HS.getUntrackedParameter<bool>("Verbosity", false);
0045   applyFidCut = m_HS.getParameter<bool>("ApplyFiducialCut");
0046   fileVersion_ = m_HS.getParameter<int>("FileVersion");
0047 
0048   if (pTreeName.find('.') == 0)
0049     pTreeName.erase(0, 2);
0050   const char* nTree = pTreeName.c_str();
0051   hf = TFile::Open(nTree);
0052 
0053   if (!hf->IsOpen()) {
0054     edm::LogError("HFShower") << "HFShowerLibrary: opening " << nTree << " failed";
0055     throw cms::Exception("Unknown", "HFShowerLibrary") << "Opening of " << pTreeName << " fails\n";
0056   } else {
0057     edm::LogVerbatim("HFShower") << "HFShowerLibrary: opening " << nTree << " successfully";
0058   }
0059 
0060   newForm = (branchEvInfo.empty());
0061   TTree* event(nullptr);
0062   if (newForm)
0063     event = (TTree*)hf->Get("HFSimHits");
0064   else
0065     event = (TTree*)hf->Get("Events");
0066   if (event) {
0067     TBranch* evtInfo(nullptr);
0068     if (!newForm) {
0069       std::string info = branchEvInfo + branchPost;
0070       evtInfo = event->GetBranch(info.c_str());
0071     }
0072     if (evtInfo || newForm) {
0073       loadEventInfo(evtInfo);
0074     } else {
0075       edm::LogError("HFShower") << "HFShowerLibrary: HFShowerLibrayEventInfo"
0076                                 << " Branch does not exist in Event";
0077       throw cms::Exception("Unknown", "HFShowerLibrary") << "Event information absent\n";
0078     }
0079   } else {
0080     edm::LogError("HFShower") << "HFShowerLibrary: Events Tree does not "
0081                               << "exist";
0082     throw cms::Exception("Unknown", "HFShowerLibrary") << "Events tree absent\n";
0083   }
0084 
0085   std::stringstream ss;
0086   ss << "HFShowerLibrary: Library " << libVers << " ListVersion " << listVersion << " File version " << fileVersion_
0087      << " Events Total " << totEvents << " and " << evtPerBin << " per bin\n";
0088   ss << "HFShowerLibrary: Energies (GeV) with " << nMomBin << " bins\n";
0089   for (int i = 0; i < nMomBin; ++i) {
0090     if (i / 10 * 10 == i && i > 0) {
0091       ss << "\n";
0092     }
0093     ss << "  " << pmom[i] / CLHEP::GeV;
0094   }
0095   edm::LogVerbatim("HFShower") << ss.str();
0096 
0097   std::string nameBr = branchPre + emName + branchPost;
0098   emBranch = event->GetBranch(nameBr.c_str());
0099   if (verbose)
0100     emBranch->Print();
0101   nameBr = branchPre + hadName + branchPost;
0102   hadBranch = event->GetBranch(nameBr.c_str());
0103   if (verbose)
0104     hadBranch->Print();
0105 
0106   v3version = false;
0107   if (emBranch->GetClassName() == std::string("vector<float>")) {
0108     v3version = true;
0109   }
0110 
0111   edm::LogVerbatim("HFShower") << " HFShowerLibrary:Branch " << emName << " has " << emBranch->GetEntries()
0112                                << " entries and Branch " << hadName << " has " << hadBranch->GetEntries()
0113                                << " entries\n HFShowerLibrary::No packing information - Assume x, y, z are not in "
0114                                   "packed form\n Maximum probability cut off "
0115                                << probMax << "  Back propagation of light probability " << backProb
0116                                << " Flag for equalizing Time Shift for different eta " << equalizeTimeShift_;
0117 
0118   fibre_ = std::make_unique<HFFibre>(name, hcalConstant_, hps, p);
0119   photo = std::make_unique<HFShowerPhotonCollection>();
0120 
0121   //Radius (minimum and maximum)
0122   std::vector<double> rTable = hcalConstant_->getRTableHF();
0123   rMin = rTable[0];
0124   rMax = rTable[rTable.size() - 1];
0125 
0126   //Delta phi
0127   std::vector<double> phibin = hcalConstant_->getPhiTableHF();
0128   dphi = phibin[0];
0129 
0130   edm::LogVerbatim("HFShower") << "HFShowerLibrary: rMIN " << rMin / CLHEP::cm << " cm and rMax " << rMax / CLHEP::cm
0131                                << " (Half) Phi Width of wedge " << dphi / CLHEP::deg;
0132 
0133   //Special Geometry parameters
0134   gpar = hcalConstant_->getGparHF();
0135 }
0136 
0137 HFShowerLibrary::~HFShowerLibrary() {
0138   if (hf)
0139     hf->Close();
0140 }
0141 
0142 std::vector<HFShowerLibrary::Hit> HFShowerLibrary::getHits(const G4Step* aStep,
0143                                                            bool& isKilled,
0144                                                            double weight,
0145                                                            bool onlyLong) {
0146   auto const preStepPoint = aStep->GetPreStepPoint();
0147   auto const postStepPoint = aStep->GetPostStepPoint();
0148   auto const track = aStep->GetTrack();
0149   // Get Z-direction
0150   auto const aParticle = track->GetDynamicParticle();
0151   const G4ThreeVector& momDir = aParticle->GetMomentumDirection();
0152   const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
0153   int parCode = track->GetDefinition()->GetPDGEncoding();
0154 
0155   // VI: for ions use internally pdg code of alpha in order to keep
0156   // consistency with previous simulation
0157   if (track->GetDefinition()->IsGeneralIon()) {
0158     parCode = 1000020040;
0159   }
0160 
0161 #ifdef EDM_ML_DEBUG
0162   G4String partType = track->GetDefinition()->GetParticleName();
0163   const G4ThreeVector localPos = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
0164   double zoff = localPos.z() + 0.5 * gpar[1];
0165 
0166   edm::LogVerbatim("HFShower") << "HFShowerLibrary::getHits " << partType << " of energy "
0167                                << track->GetKineticEnergy() / CLHEP::GeV << " GeV weight= " << weight
0168                                << " onlyLong: " << onlyLong << "  dir.orts " << momDir.x() << ", " << momDir.y() << ", "
0169                                << momDir.z() << "  Pos x,y,z = " << hitPoint.x() << "," << hitPoint.y() << ","
0170                                << hitPoint.z() << " (" << zoff << ")   sphi,cphi,stheta,ctheta  = " << sin(momDir.phi())
0171                                << "," << cos(momDir.phi()) << ", " << sin(momDir.theta()) << "," << cos(momDir.theta());
0172 #endif
0173 
0174   double tSlice = (postStepPoint->GetGlobalTime()) / CLHEP::nanosecond;
0175 
0176   // use kinetic energy for protons and ions
0177   double pin = (track->GetDefinition()->GetBaryonNumber() > 0) ? preStepPoint->GetKineticEnergy()
0178                                                                : preStepPoint->GetTotalEnergy();
0179 
0180   return fillHits(hitPoint, momDir, parCode, pin, isKilled, weight, tSlice, onlyLong);
0181 }
0182 
0183 std::vector<HFShowerLibrary::Hit> HFShowerLibrary::fillHits(const G4ThreeVector& hitPoint,
0184                                                             const G4ThreeVector& momDir,
0185                                                             int parCode,
0186                                                             double pin,
0187                                                             bool& ok,
0188                                                             double weight,
0189                                                             double tSlice,
0190                                                             bool onlyLong) {
0191   std::vector<HFShowerLibrary::Hit> hit;
0192   ok = false;
0193   bool isEM = G4TrackToParticleID::isGammaElectronPositron(parCode);
0194   // shower is built only for gamma, e+- and stable hadrons
0195   if (!isEM && !G4TrackToParticleID::isStableHadron(parCode)) {
0196     return hit;
0197   }
0198   ok = true;
0199 
0200   // remove low-energy component
0201   const double threshold = 50 * MeV;
0202   if (pin < threshold) {
0203     return hit;
0204   }
0205 
0206   double pz = momDir.z();
0207   double zint = hitPoint.z();
0208 
0209   // if particle moves from interaction point or "backwards (halo)
0210   bool backward = (pz * zint < 0.) ? true : false;
0211 
0212   double sphi = sin(momDir.phi());
0213   double cphi = cos(momDir.phi());
0214   double ctheta = cos(momDir.theta());
0215   double stheta = sin(momDir.theta());
0216 
0217   if (isEM) {
0218     if (pin < pmom[nMomBin - 1]) {
0219       interpolate(0, pin);
0220     } else {
0221       extrapolate(0, pin);
0222     }
0223   } else {
0224     if (pin < pmom[nMomBin - 1]) {
0225       interpolate(1, pin);
0226     } else {
0227       extrapolate(1, pin);
0228     }
0229   }
0230 
0231   int nHit = 0;
0232   HFShowerLibrary::Hit oneHit;
0233   for (int i = 0; i < npe; ++i) {
0234     double zv = std::abs(pe[i].z());  // abs local z
0235 #ifdef EDM_ML_DEBUG
0236     edm::LogVerbatim("HFShower") << "HFShowerLibrary: Hit " << i << " " << pe[i] << " zv " << zv;
0237 #endif
0238     if (zv <= gpar[1] && pe[i].lambda() > 0 && (pe[i].z() >= 0 || (zv > gpar[0] && (!onlyLong)))) {
0239       int depth = 1;
0240       if (onlyLong) {
0241       } else if (!backward) {  // fully valid only for "front" particles
0242         if (pe[i].z() < 0)
0243           depth = 2;                 // with "front"-simulated shower lib.
0244       } else {                       // for "backward" particles - almost equal
0245         double r = G4UniformRand();  // share between L and S fibers
0246         if (r > 0.5)
0247           depth = 2;
0248       }
0249 
0250       // Updated coordinate transformation from local
0251       //  back to global using two Euler angles: phi and theta
0252       double pex = pe[i].x();
0253       double pey = pe[i].y();
0254 
0255       double xx = pex * ctheta * cphi - pey * sphi + zv * stheta * cphi;
0256       double yy = pex * ctheta * sphi + pey * cphi + zv * stheta * sphi;
0257       double zz = -pex * stheta + zv * ctheta;
0258 
0259       G4ThreeVector pos = hitPoint + G4ThreeVector(xx, yy, zz);
0260       zv = std::abs(pos.z()) - gpar[4] - 0.5 * gpar[1];
0261       G4ThreeVector lpos = G4ThreeVector(pos.x(), pos.y(), zv);
0262 
0263       zv = fibre_->zShift(lpos, depth, 0);  // distance to PMT !
0264 
0265       double r = pos.perp();
0266       double p = fibre_->attLength(pe[i].lambda());
0267       double fi = pos.phi();
0268       if (fi < 0)
0269         fi += CLHEP::twopi;
0270       int isect = int(fi / dphi) + 1;
0271       isect = (isect + 1) / 2;
0272       double dfi = ((isect * 2 - 1) * dphi - fi);
0273       if (dfi < 0)
0274         dfi = -dfi;
0275       double dfir = r * sin(dfi);
0276 #ifdef EDM_ML_DEBUG
0277       edm::LogVerbatim("HFShower") << "HFShowerLibrary: Position shift " << xx << ", " << yy << ", " << zz << ": "
0278                                    << pos << " R " << r << " Phi " << fi << " Section " << isect << " R*Dfi " << dfir
0279                                    << " Dist " << zv;
0280 #endif
0281       zz = std::abs(pos.z());
0282       double r1 = G4UniformRand();
0283       double r2 = G4UniformRand();
0284       double r3 = backward ? G4UniformRand() : -9999.;
0285       if (!applyFidCut)
0286         dfir += gpar[5];
0287 
0288 #ifdef EDM_ML_DEBUG
0289       edm::LogVerbatim("HFShower") << "HFShowerLibrary: rLimits " << rInside(r) << " attenuation " << r1 << ":"
0290                                    << exp(-p * zv) << " r2 " << r2 << " r3 " << r3 << " rDfi " << gpar[5] << " zz "
0291                                    << zz << " zLim " << gpar[4] << ":" << gpar[4] + gpar[1] << "\n"
0292                                    << "  rInside(r) :" << rInside(r) << "  r1 <= exp(-p*zv) :" << (r1 <= exp(-p * zv))
0293                                    << "  r2 <= probMax :" << (r2 <= probMax * weight)
0294                                    << "  r3 <= backProb :" << (r3 <= backProb)
0295                                    << "  dfir > gpar[5] :" << (dfir > gpar[5]) << "  zz >= gpar[4] :" << (zz >= gpar[4])
0296                                    << "  zz <= gpar[4]+gpar[1] :" << (zz <= gpar[4] + gpar[1]);
0297 #endif
0298       if (rInside(r) && r1 <= exp(-p * zv) && r2 <= probMax * weight && dfir > gpar[5] && zz >= gpar[4] &&
0299           zz <= gpar[4] + gpar[1] && r3 <= backProb && (depth != 2 || zz >= gpar[4] + gpar[0])) {
0300         double tdiff = (equalizeTimeShift_) ? (fibre_->tShift(lpos, depth, -1)) : (fibre_->tShift(lpos, depth, 1));
0301         oneHit.position = pos;
0302         oneHit.depth = depth;
0303         oneHit.time = (tSlice + (pe[i].t()) + tdiff);
0304         hit.push_back(oneHit);
0305 
0306 #ifdef EDM_ML_DEBUG
0307         edm::LogVerbatim("HFShower") << "HFShowerLibrary: Final Hit " << nHit << " position " << (hit[nHit].position)
0308                                      << " Depth " << (hit[nHit].depth) << " Time " << tSlice << ":" << pe[i].t() << ":"
0309                                      << tdiff << ":" << (hit[nHit].time);
0310 #endif
0311         ++nHit;
0312       }
0313 #ifdef EDM_ML_DEBUG
0314       else
0315         edm::LogVerbatim("HFShower") << "HFShowerLibrary: REJECTED !!!";
0316 #endif
0317       if (onlyLong && zz >= gpar[4] + gpar[0] && zz <= gpar[4] + gpar[1]) {
0318         r1 = G4UniformRand();
0319         r2 = G4UniformRand();
0320         if (rInside(r) && r1 <= exp(-p * zv) && r2 <= probMax && dfir > gpar[5]) {
0321           double tdiff = (equalizeTimeShift_) ? (fibre_->tShift(lpos, 2, -1)) : (fibre_->tShift(lpos, 2, 1));
0322           oneHit.position = pos;
0323           oneHit.depth = 2;
0324           oneHit.time = (tSlice + (pe[i].t()) + tdiff);
0325           hit.push_back(oneHit);
0326 #ifdef EDM_ML_DEBUG
0327           edm::LogVerbatim("HFShower") << "HFShowerLibrary: Final Hit " << nHit << " position " << (hit[nHit].position)
0328                                        << " Depth " << (hit[nHit].depth) << " Time " << tSlice << ":" << pe[i].t()
0329                                        << ":" << tdiff << ":" << (hit[nHit].time);
0330 #endif
0331           ++nHit;
0332         }
0333       }
0334     }
0335   }
0336 
0337 #ifdef EDM_ML_DEBUG
0338   edm::LogVerbatim("HFShower") << "HFShowerLibrary: Total Hits " << nHit << " out of " << npe << " PE";
0339 #endif
0340   if (nHit > npe && !onlyLong) {
0341     edm::LogWarning("HFShower") << "HFShowerLibrary: Hit buffer " << npe << " smaller than " << nHit << " Hits";
0342   }
0343   return hit;
0344 }
0345 
0346 bool HFShowerLibrary::rInside(double r) { return (r >= rMin && r <= rMax); }
0347 
0348 void HFShowerLibrary::getRecord(int type, int record) {
0349   int nrc = record - 1;
0350   photon.clear();
0351   photo->clear();
0352   if (type > 0) {
0353     if (newForm) {
0354       if (!v3version) {
0355         hadBranch->SetAddress(&photo);
0356         int position = (fileVersion_ >= 2) ? nrc : (nrc + totEvents);
0357         hadBranch->GetEntry(position);
0358       } else {
0359         std::vector<float> t;
0360         std::vector<float>* tp = &t;
0361         hadBranch->SetAddress(&tp);
0362         hadBranch->GetEntry(nrc + totEvents);
0363         unsigned int tSize = t.size() / 5;
0364         photo->reserve(tSize);
0365         for (unsigned int i = 0; i < tSize; i++) {
0366           photo->push_back(
0367               HFShowerPhoton(t[i], t[1 * tSize + i], t[2 * tSize + i], t[3 * tSize + i], t[4 * tSize + i]));
0368         }
0369       }
0370     } else {
0371       hadBranch->SetAddress(&photon);
0372       hadBranch->GetEntry(nrc);
0373     }
0374   } else {
0375     if (newForm) {
0376       if (!v3version) {
0377         emBranch->SetAddress(&photo);
0378         emBranch->GetEntry(nrc);
0379       } else {
0380         std::vector<float> t;
0381         std::vector<float>* tp = &t;
0382         emBranch->SetAddress(&tp);
0383         emBranch->GetEntry(nrc);
0384         unsigned int tSize = t.size() / 5;
0385         photo->reserve(tSize);
0386         for (unsigned int i = 0; i < tSize; i++) {
0387           photo->push_back(
0388               HFShowerPhoton(t[i], t[1 * tSize + i], t[2 * tSize + i], t[3 * tSize + i], t[4 * tSize + i]));
0389         }
0390       }
0391     } else {
0392       emBranch->SetAddress(&photon);
0393       emBranch->GetEntry(nrc);
0394     }
0395   }
0396 #ifdef EDM_ML_DEBUG
0397   int nPhoton = (newForm) ? photo->size() : photon.size();
0398   edm::LogVerbatim("HFShower") << "HFShowerLibrary::getRecord: Record " << record << " of type " << type << " with "
0399                                << nPhoton << " photons";
0400   for (int j = 0; j < nPhoton; j++)
0401     if (newForm)
0402       edm::LogVerbatim("HFShower") << "Photon " << j << " " << photo->at(j);
0403     else
0404       edm::LogVerbatim("HFShower") << "Photon " << j << " " << photon[j];
0405 #endif
0406 }
0407 
0408 void HFShowerLibrary::loadEventInfo(TBranch* branch) {
0409   if (branch) {
0410     std::vector<HFShowerLibraryEventInfo> eventInfoCollection;
0411     branch->SetAddress(&eventInfoCollection);
0412     branch->GetEntry(0);
0413     edm::LogVerbatim("HFShower") << "HFShowerLibrary::loadEventInfo loads EventInfo Collection of size "
0414                                  << eventInfoCollection.size() << " records";
0415 
0416     totEvents = eventInfoCollection[0].totalEvents();
0417     nMomBin = eventInfoCollection[0].numberOfBins();
0418     evtPerBin = eventInfoCollection[0].eventsPerBin();
0419     libVers = eventInfoCollection[0].showerLibraryVersion();
0420     listVersion = eventInfoCollection[0].physListVersion();
0421     pmom = eventInfoCollection[0].energyBins();
0422   } else {
0423     edm::LogVerbatim("HFShower") << "HFShowerLibrary::loadEventInfo loads EventInfo from hardwired"
0424                                  << " numbers";
0425 
0426     nMomBin = 16;
0427     evtPerBin = (fileVersion_ == 0) ? 5000 : 10000;
0428     totEvents = nMomBin * evtPerBin;
0429     libVers = (fileVersion_ == 0) ? 1.1 : 1.2;
0430     listVersion = 3.6;
0431     pmom = {2, 3, 5, 7, 10, 15, 20, 30, 50, 75, 100, 150, 250, 350, 500, 1000};
0432   }
0433   for (int i = 0; i < nMomBin; i++)
0434     pmom[i] *= CLHEP::GeV;
0435 }
0436 
0437 void HFShowerLibrary::interpolate(int type, double pin) {
0438 #ifdef EDM_ML_DEBUG
0439   edm::LogVerbatim("HFShower") << "HFShowerLibrary:: Interpolate for Energy " << pin / CLHEP::GeV << " GeV with "
0440                                << nMomBin << " momentum bins and " << evtPerBin << " entries/bin -- total "
0441                                << totEvents;
0442 #endif
0443   int irc[2] = {0, 0};
0444   double w = 0.;
0445   double r = G4UniformRand();
0446 
0447   if (pin < pmom[0]) {
0448     w = pin / pmom[0];
0449     irc[1] = int(evtPerBin * r) + 1;
0450     irc[0] = 0;
0451   } else {
0452     for (int j = 0; j < nMomBin - 1; j++) {
0453       if (pin >= pmom[j] && pin < pmom[j + 1]) {
0454         w = (pin - pmom[j]) / (pmom[j + 1] - pmom[j]);
0455         if (j == nMomBin - 2) {
0456           irc[1] = int(evtPerBin * 0.5 * r);
0457         } else {
0458           irc[1] = int(evtPerBin * r);
0459         }
0460         irc[1] += (j + 1) * evtPerBin + 1;
0461         r = G4UniformRand();
0462         irc[0] = int(evtPerBin * r) + 1 + j * evtPerBin;
0463         if (irc[0] < 0) {
0464           edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[0] = " << irc[0] << " now set to 0";
0465           irc[0] = 0;
0466         } else if (irc[0] > totEvents) {
0467           edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[0] = " << irc[0] << " now set to " << totEvents;
0468           irc[0] = totEvents;
0469         }
0470       }
0471     }
0472   }
0473   if (irc[1] < 1) {
0474     edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[1] = " << irc[1] << " now set to 1";
0475     irc[1] = 1;
0476   } else if (irc[1] > totEvents) {
0477     edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[1] = " << irc[1] << " now set to " << totEvents;
0478     irc[1] = totEvents;
0479   }
0480 
0481 #ifdef EDM_ML_DEBUG
0482   edm::LogVerbatim("HFShower") << "HFShowerLibrary:: Select records " << irc[0] << " and " << irc[1] << " with weights "
0483                                << 1 - w << " and " << w;
0484 #endif
0485   pe.clear();
0486   npe = 0;
0487   int npold = 0;
0488   for (int ir = 0; ir < 2; ir++) {
0489     if (irc[ir] > 0) {
0490       getRecord(type, irc[ir]);
0491       int nPhoton = (newForm) ? photo->size() : photon.size();
0492       npold += nPhoton;
0493       for (int j = 0; j < nPhoton; j++) {
0494         r = G4UniformRand();
0495         if ((ir == 0 && r > w) || (ir > 0 && r < w)) {
0496           storePhoton(j);
0497         }
0498       }
0499     }
0500   }
0501 
0502   if ((npe > npold || (npold == 0 && irc[0] > 0)) && !(npe == 0 && npold == 0))
0503     edm::LogWarning("HFShower") << "HFShowerLibrary: Interpolation Warning =="
0504                                 << " records " << irc[0] << " and " << irc[1] << " gives a buffer of " << npold
0505                                 << " photons and fills " << npe << " *****";
0506 #ifdef EDM_ML_DEBUG
0507   else
0508     edm::LogVerbatim("HFShower") << "HFShowerLibrary: Interpolation == records " << irc[0] << " and " << irc[1]
0509                                  << " gives a buffer of " << npold << " photons and fills " << npe << " PE";
0510   for (int j = 0; j < npe; j++)
0511     edm::LogVerbatim("HFShower") << "Photon " << j << " " << pe[j];
0512 #endif
0513 }
0514 
0515 void HFShowerLibrary::extrapolate(int type, double pin) {
0516   int nrec = int(pin / pmom[nMomBin - 1]);
0517   double w = (pin - pmom[nMomBin - 1] * nrec) / pmom[nMomBin - 1];
0518   nrec++;
0519 #ifdef EDM_ML_DEBUG
0520   edm::LogVerbatim("HFShower") << "HFShowerLibrary:: Extrapolate for Energy " << pin / CLHEP::GeV << " GeV with "
0521                                << nMomBin << " momentum bins and " << evtPerBin << " entries/bin -- "
0522                                << "total " << totEvents << " using " << nrec << " records";
0523 #endif
0524   std::vector<int> irc(nrec);
0525 
0526   for (int ir = 0; ir < nrec; ir++) {
0527     double r = G4UniformRand();
0528     irc[ir] = int(evtPerBin * 0.5 * r) + (nMomBin - 1) * evtPerBin + 1;
0529     if (irc[ir] < 1) {
0530       edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[" << ir << "] = " << irc[ir] << " now set to 1";
0531       irc[ir] = 1;
0532     } else if (irc[ir] > totEvents) {
0533       edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[" << ir << "] = " << irc[ir] << " now set to "
0534                                   << totEvents;
0535       irc[ir] = totEvents;
0536 #ifdef EDM_ML_DEBUG
0537     } else {
0538       edm::LogVerbatim("HFShower") << "HFShowerLibrary::Extrapolation use irc[" << ir << "] = " << irc[ir];
0539 #endif
0540     }
0541   }
0542 
0543   pe.clear();
0544   npe = 0;
0545   int npold = 0;
0546   for (int ir = 0; ir < nrec; ir++) {
0547     if (irc[ir] > 0) {
0548       getRecord(type, irc[ir]);
0549       int nPhoton = (newForm) ? photo->size() : photon.size();
0550       npold += nPhoton;
0551       for (int j = 0; j < nPhoton; j++) {
0552         double r = G4UniformRand();
0553         if (ir != nrec - 1 || r < w) {
0554           storePhoton(j);
0555         }
0556       }
0557 #ifdef EDM_ML_DEBUG
0558       edm::LogVerbatim("HFShower") << "HFShowerLibrary: Record [" << ir << "] = " << irc[ir] << " npold = " << npold;
0559 #endif
0560     }
0561   }
0562 #ifdef EDM_ML_DEBUG
0563   edm::LogVerbatim("HFShower") << "HFShowerLibrary:: uses " << npold << " photons";
0564 #endif
0565 
0566   if (npe > npold || npold == 0)
0567     edm::LogWarning("HFShower") << "HFShowerLibrary: Extrapolation Warning == " << nrec << " records " << irc[0] << ", "
0568                                 << irc[1] << ", ... gives a buffer of " << npold << " photons and fills " << npe
0569                                 << " *****";
0570 #ifdef EDM_ML_DEBUG
0571   else
0572     edm::LogVerbatim("HFShower") << "HFShowerLibrary: Extrapolation == " << nrec << " records " << irc[0] << ", "
0573                                  << irc[1] << ", ... gives a buffer of " << npold << " photons and fills " << npe
0574                                  << " PE";
0575   for (int j = 0; j < npe; j++)
0576     edm::LogVerbatim("HFShower") << "Photon " << j << " " << pe[j];
0577 #endif
0578 }
0579 
0580 void HFShowerLibrary::storePhoton(int j) {
0581   if (newForm)
0582     pe.push_back(photo->at(j));
0583   else
0584     pe.push_back(photon[j]);
0585   npe++;
0586 #ifdef EDM_ML_DEBUG
0587   edm::LogVerbatim("HFShower") << "HFShowerLibrary: storePhoton " << j << " npe " << npe << " " << pe[npe - 1];
0588 #endif
0589 }