Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-05-10 02:21:17 

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 namespace {
0023   HFShowerLibrary::Params paramsFrom(edm::ParameterSet const& hfShower,
0024                                      edm::ParameterSet const& hfShowerLibrary,
0025                                      double iDeltaPhi) {
0026     HFShowerLibrary::Params params;
0027 
0028     params.dphi_ = iDeltaPhi;
0029 
0030     params.probMax_ = hfShower.getParameter<double>("ProbMax");
0031     params.equalizeTimeShift_ = hfShower.getParameter<bool>("EqualizeTimeShift");
0032 
0033     params.backProb_ = hfShowerLibrary.getParameter<double>("BackProbability");
0034     params.verbose_ = hfShowerLibrary.getUntrackedParameter<bool>("Verbosity", false);
0035     params.applyFidCut_ = hfShowerLibrary.getParameter<bool>("ApplyFiducialCut");
0036 
0037     return params;
0038   }
0039 
0040   HFShowerLibrary::FileParams fileParamsFrom(edm::ParameterSet const& hfShowerLibrary) {
0041     HFShowerLibrary::FileParams params;
0042 
0043     edm::FileInPath fp = hfShowerLibrary.getParameter<edm::FileInPath>("FileName");
0044     params.fileName_ = fp.fullPath();
0045     if (params.fileName_.find('.') == 0)
0046       params.fileName_.erase(0, 2);
0047 
0048     std::string emName = hfShowerLibrary.getParameter<std::string>("TreeEMID");
0049     std::string hadName = hfShowerLibrary.getParameter<std::string>("TreeHadID");
0050     std::string branchEvInfo = hfShowerLibrary.getUntrackedParameter<std::string>(
0051         "BranchEvt", "HFShowerLibraryEventInfos_hfshowerlib_HFShowerLibraryEventInfo");
0052     std::string branchPre =
0053         hfShowerLibrary.getUntrackedParameter<std::string>("BranchPre", "HFShowerPhotons_hfshowerlib_");
0054     std::string branchPost = hfShowerLibrary.getUntrackedParameter<std::string>("BranchPost", "_R.obj");
0055     params.fileVersion_ = hfShowerLibrary.getParameter<int>("FileVersion");
0056 
0057     params.emBranchName_ = branchPre + emName + branchPost;
0058     params.hadBranchName_ = branchPre + hadName + branchPost;
0059 
0060     if (not branchEvInfo.empty()) {
0061       params.branchEvInfo_ = branchEvInfo + branchPost;
0062     }
0063 
0064     params.cacheBranches_ = hfShowerLibrary.getUntrackedParameter<bool>("cacheBranches", false);
0065     return params;
0066   }
0067 }  // namespace
0068 
0069 HFShowerLibrary::HFShowerLibrary(const HcalDDDSimConstants* hcons,
0070                                  const HcalSimulationParameters* hps,
0071                                  edm::ParameterSet const& hfShower,
0072                                  edm::ParameterSet const& hfShowerLibrary)
0073     : HFShowerLibrary(paramsFrom(hfShower, hfShowerLibrary, hcons->getPhiTableHF().front()),
0074                       fileParamsFrom(hfShowerLibrary),
0075                       HFFibre::Params(hfShower.getParameter<double>("CFibre"), hcons, hps)) {}
0076 
0077 HFShowerLibrary::HFShowerLibrary(const std::string& name,
0078                                  const HcalDDDSimConstants* hcons,
0079                                  const HcalSimulationParameters* hps,
0080                                  edm::ParameterSet const& p)
0081     : HFShowerLibrary(
0082           hcons,
0083           hps,
0084           p.getParameter<edm::ParameterSet>("HFShower").getParameter<edm::ParameterSet>("HFShowerBlock"),
0085           p.getParameter<edm::ParameterSet>("HFShowerLibrary").getParameter<edm::ParameterSet>("HFLibraryFileBlock")) {}
0086 
0087 HFShowerLibrary::HFShowerLibrary(const Params& iParams, const FileParams& iFileParams, HFFibre::Params iFibreParams)
0088     : fibre_(iFibreParams),
0089       hf_(),
0090       emBranch_(),
0091       hadBranch_(),
0092       verbose_(iParams.verbose_),
0093       applyFidCut_(iParams.applyFidCut_),
0094       equalizeTimeShift_(iParams.equalizeTimeShift_),
0095       probMax_(iParams.probMax_),
0096       backProb_(iParams.backProb_),
0097       dphi_(iParams.dphi_),
0098       rMin_(iFibreParams.rTableHF_.front()),
0099       rMax_(iFibreParams.rTableHF_.back()),
0100       gpar_(iFibreParams.gParHF_) {
0101   std::string pTreeName = iFileParams.fileName_;
0102 
0103   const char* nTree = pTreeName.c_str();
0104   {
0105     //It is not safe to open a Tfile in one thread and close in another without adding the following:
0106     TDirectory::TContext context;
0107     hf_ = std::unique_ptr<TFile>(TFile::Open(nTree));
0108   }
0109   if (!hf_->IsOpen()) {
0110     edm::LogError("HFShower") << "HFShowerLibrary: opening " << nTree << " failed";
0111     throw cms::Exception("Unknown", "HFShowerLibrary") << "Opening of " << pTreeName << " fails\n";
0112   } else {
0113     edm::LogVerbatim("HFShower") << "HFShowerLibrary: opening " << nTree << " successfully";
0114   }
0115 
0116   auto fileFormat = FileFormat::kOld;
0117   const int fileVersion = iFileParams.fileVersion_;
0118 
0119   auto newForm = iFileParams.branchEvInfo_.empty();
0120   TTree* event(nullptr);
0121   if (newForm) {
0122     fileFormat = FileFormat::kNew;
0123     event = (TTree*)hf_->Get("HFSimHits");
0124   } else {
0125     event = (TTree*)hf_->Get("Events");
0126   }
0127   VersionInfo versionInfo;
0128   if (event) {
0129     TBranch* evtInfo(nullptr);
0130     if (!newForm) {
0131       std::string info = iFileParams.branchEvInfo_;
0132       evtInfo = event->GetBranch(info.c_str());
0133     }
0134     if (evtInfo || newForm) {
0135       versionInfo = loadEventInfo(evtInfo, fileVersion);
0136     } else {
0137       edm::LogError("HFShower") << "HFShowerLibrary: HFShowerLibrayEventInfo"
0138                                 << " Branch does not exist in Event";
0139       throw cms::Exception("Unknown", "HFShowerLibrary") << "Event information absent\n";
0140     }
0141   } else {
0142     edm::LogError("HFShower") << "HFShowerLibrary: Events Tree does not "
0143                               << "exist";
0144     throw cms::Exception("Unknown", "HFShowerLibrary") << "Events tree absent\n";
0145   }
0146 
0147   edm::LogVerbatim("HFShower").log([&](auto& logger) {
0148     logger << "HFShowerLibrary: Library " << versionInfo.libVers_ << " ListVersion " << versionInfo.listVersion_
0149            << " File version " << fileVersion << " Events Total " << totEvents_ << " and " << evtPerBin_
0150            << " per bin\n";
0151     logger << "HFShowerLibrary: Energies (GeV) with " << nMomBin_ << " bins\n";
0152     for (int i = 0; i < nMomBin_; ++i) {
0153       if (i / 10 * 10 == i && i > 0) {
0154         logger << "\n";
0155       }
0156       logger << "  " << pmom_[i] / CLHEP::GeV;
0157     }
0158   });
0159 
0160   std::string nameBr = iFileParams.emBranchName_;
0161   auto emBranch = event->GetBranch(nameBr.c_str());
0162   if (verbose_)
0163     emBranch->Print();
0164   nameBr = iFileParams.hadBranchName_;
0165   auto hadBranch = event->GetBranch(nameBr.c_str());
0166   if (verbose_)
0167     hadBranch->Print();
0168 
0169   if (emBranch->GetClassName() == std::string("vector<float>")) {
0170     assert(fileFormat == FileFormat::kNew);
0171     fileFormat = FileFormat::kNewV3;
0172   }
0173   emBranch_ = BranchReader(emBranch, fileFormat, 0, iFileParams.cacheBranches_ ? totEvents_ : 0);
0174   size_t offset = 0;
0175   if ((fileFormat == FileFormat::kNewV3 && fileVersion < 3) || (fileFormat == FileFormat::kNew && fileVersion < 2)) {
0176     //NOTE: for this format, the hadBranch is all empty up to
0177     // totEvents_ (which is more like 1/2*GenEntries())
0178     offset = totEvents_;
0179   }
0180   hadBranch_ = BranchReader(hadBranch, fileFormat, offset, iFileParams.cacheBranches_ ? totEvents_ : 0);
0181 
0182   edm::LogVerbatim("HFShower") << " HFShowerLibrary:Branch " << iFileParams.emBranchName_ << " has "
0183                                << emBranch->GetEntries() << " entries and Branch " << iFileParams.hadBranchName_
0184                                << " has " << hadBranch->GetEntries()
0185                                << " entries\n HFShowerLibrary::No packing information - Assume x, y, z are not in "
0186                                   "packed form\n Maximum probability cut off "
0187                                << probMax_ << "  Back propagation of light probability " << backProb_
0188                                << " Flag for equalizing Time Shift for different eta " << equalizeTimeShift_;
0189 
0190   edm::LogVerbatim("HFShower") << "HFShowerLibrary: rMIN " << rMin_ / CLHEP::cm << " cm and rMax " << rMax_ / CLHEP::cm
0191                                << " (Half) Phi Width of wedge " << dphi_ / CLHEP::deg;
0192   if (iFileParams.cacheBranches_) {
0193     hf_.reset();
0194   }
0195 }
0196 
0197 HFShowerLibrary::~HFShowerLibrary() {
0198   if (hf_)
0199     hf_->Close();
0200 }
0201 
0202 std::vector<HFShowerLibrary::Hit> HFShowerLibrary::getHits(const G4Step* aStep,
0203                                                            bool& isKilled,
0204                                                            double weight,
0205                                                            bool onlyLong) {
0206   auto const preStepPoint = aStep->GetPreStepPoint();
0207   auto const postStepPoint = aStep->GetPostStepPoint();
0208   auto const track = aStep->GetTrack();
0209   // Get Z-direction
0210   auto const aParticle = track->GetDynamicParticle();
0211   const G4ThreeVector& momDir = aParticle->GetMomentumDirection();
0212   const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
0213   int parCode = track->GetDefinition()->GetPDGEncoding();
0214 
0215   // VI: for ions use internally pdg code of alpha in order to keep
0216   // consistency with previous simulation
0217   if (track->GetDefinition()->IsGeneralIon()) {
0218     parCode = 1000020040;
0219   }
0220 
0221 #ifdef EDM_ML_DEBUG
0222   G4String partType = track->GetDefinition()->GetParticleName();
0223   const G4ThreeVector localPos = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
0224   double zoff = localPos.z() + 0.5 * gpar_[1];
0225 
0226   edm::LogVerbatim("HFShower") << "HFShowerLibrary::getHits " << partType << " of energy "
0227                                << track->GetKineticEnergy() / CLHEP::GeV << " GeV weight= " << weight
0228                                << " onlyLong: " << onlyLong << "  dir.orts " << momDir.x() << ", " << momDir.y() << ", "
0229                                << momDir.z() << "  Pos x,y,z = " << hitPoint.x() << "," << hitPoint.y() << ","
0230                                << hitPoint.z() << " (" << zoff << ")   sphi,cphi,stheta,ctheta  = " << sin(momDir.phi())
0231                                << "," << cos(momDir.phi()) << ", " << sin(momDir.theta()) << "," << cos(momDir.theta());
0232 #endif
0233 
0234   double tSlice = (postStepPoint->GetGlobalTime()) / CLHEP::nanosecond;
0235 
0236   // use kinetic energy for protons and ions
0237   double pin = (track->GetDefinition()->GetBaryonNumber() > 0) ? preStepPoint->GetKineticEnergy()
0238                                                                : preStepPoint->GetTotalEnergy();
0239 
0240   return fillHits(hitPoint, momDir, parCode, pin, isKilled, weight, tSlice, onlyLong);
0241 }
0242 
0243 std::vector<HFShowerLibrary::Hit> HFShowerLibrary::fillHits(const G4ThreeVector& hitPoint,
0244                                                             const G4ThreeVector& momDir,
0245                                                             int parCode,
0246                                                             double pin,
0247                                                             bool& ok,
0248                                                             double weight,
0249                                                             double tSlice,
0250                                                             bool onlyLong) {
0251   std::vector<HFShowerLibrary::Hit> hit;
0252   ok = false;
0253   bool isEM = G4TrackToParticleID::isGammaElectronPositron(parCode);
0254   // shower is built only for gamma, e+- and stable hadrons
0255   if (!isEM && !G4TrackToParticleID::isStableHadron(parCode)) {
0256     return hit;
0257   }
0258   ok = true;
0259 
0260   // remove low-energy component
0261   const double threshold = 50 * CLHEP::MeV;
0262   if (pin < threshold) {
0263     return hit;
0264   }
0265 
0266   double pz = momDir.z();
0267   double zint = hitPoint.z();
0268 
0269   // if particle moves from interaction point or "backwards (halo)
0270   bool backward = (pz * zint < 0.) ? true : false;
0271 
0272   double sphi = sin(momDir.phi());
0273   double cphi = cos(momDir.phi());
0274   double ctheta = cos(momDir.theta());
0275   double stheta = sin(momDir.theta());
0276 
0277   HFShowerPhotonCollection pe;
0278   if (isEM) {
0279     if (pin < pmom_[nMomBin_ - 1]) {
0280       pe = interpolate(0, pin);
0281     } else {
0282       pe = extrapolate(0, pin);
0283     }
0284   } else {
0285     if (pin < pmom_[nMomBin_ - 1]) {
0286       pe = interpolate(1, pin);
0287     } else {
0288       pe = extrapolate(1, pin);
0289     }
0290   }
0291 
0292   std::size_t nHit = 0;
0293   HFShowerLibrary::Hit oneHit;
0294 #ifdef EDM_ML_DEBUG
0295   int i = 0;
0296 #endif
0297   for (auto const& photon : pe) {
0298     double zv = std::abs(photon.z());  // abs local z
0299 #ifdef EDM_ML_DEBUG
0300     edm::LogVerbatim("HFShower") << "HFShowerLibrary: Hit " << i++ << " " << photon << " zv " << zv;
0301 #endif
0302     if (zv <= gpar_[1] && photon.lambda() > 0 && (photon.z() >= 0 || (zv > gpar_[0] && (!onlyLong)))) {
0303       int depth = 1;
0304       if (onlyLong) {
0305       } else if (!backward) {  // fully valid only for "front" particles
0306         if (photon.z() < 0)
0307           depth = 2;                 // with "front"-simulated shower lib.
0308       } else {                       // for "backward" particles - almost equal
0309         double r = G4UniformRand();  // share between L and S fibers
0310         if (r > 0.5)
0311           depth = 2;
0312       }
0313 
0314       // Updated coordinate transformation from local
0315       //  back to global using two Euler angles: phi and theta
0316       double pex = photon.x();
0317       double pey = photon.y();
0318 
0319       double xx = pex * ctheta * cphi - pey * sphi + zv * stheta * cphi;
0320       double yy = pex * ctheta * sphi + pey * cphi + zv * stheta * sphi;
0321       double zz = -pex * stheta + zv * ctheta;
0322 
0323       G4ThreeVector pos = hitPoint + G4ThreeVector(xx, yy, zz);
0324       zv = std::abs(pos.z()) - gpar_[4] - 0.5 * gpar_[1];
0325       G4ThreeVector lpos = G4ThreeVector(pos.x(), pos.y(), zv);
0326 
0327       zv = fibre_.zShift(lpos, depth, 0);  // distance to PMT !
0328 
0329       double r = pos.perp();
0330       double p = fibre_.attLength(photon.lambda());
0331       double fi = pos.phi();
0332       if (fi < 0)
0333         fi += CLHEP::twopi;
0334       int isect = int(fi / dphi_) + 1;
0335       isect = (isect + 1) / 2;
0336       double dfi = ((isect * 2 - 1) * dphi_ - fi);
0337       if (dfi < 0)
0338         dfi = -dfi;
0339       double dfir = r * sin(dfi);
0340 #ifdef EDM_ML_DEBUG
0341       edm::LogVerbatim("HFShower") << "HFShowerLibrary: Position shift " << xx << ", " << yy << ", " << zz << ": "
0342                                    << pos << " R " << r << " Phi " << fi << " Section " << isect << " R*Dfi " << dfir
0343                                    << " Dist " << zv;
0344 #endif
0345       zz = std::abs(pos.z());
0346       double r1 = G4UniformRand();
0347       double r2 = G4UniformRand();
0348       double r3 = backward ? G4UniformRand() : -9999.;
0349       if (!applyFidCut_)
0350         dfir += gpar_[5];
0351 
0352 #ifdef EDM_ML_DEBUG
0353       edm::LogVerbatim("HFShower") << "HFShowerLibrary: rLimits " << rInside(r) << " attenuation " << r1 << ":"
0354                                    << exp(-p * zv) << " r2 " << r2 << " r3 " << r3 << " rDfi " << gpar_[5] << " zz "
0355                                    << zz << " zLim " << gpar_[4] << ":" << gpar_[4] + gpar_[1] << "\n"
0356                                    << "  rInside(r) :" << rInside(r) << "  r1 <= exp(-p*zv) :" << (r1 <= exp(-p * zv))
0357                                    << "  r2 <= probMax :" << (r2 <= probMax_ * weight)
0358                                    << "  r3 <= backProb :" << (r3 <= backProb_)
0359                                    << "  dfir > gpar[5] :" << (dfir > gpar_[5])
0360                                    << "  zz >= gpar[4] :" << (zz >= gpar_[4])
0361                                    << "  zz <= gpar[4]+gpar[1] :" << (zz <= gpar_[4] + gpar_[1]);
0362 #endif
0363       if (rInside(r) && r1 <= exp(-p * zv) && r2 <= probMax_ * weight && dfir > gpar_[5] && zz >= gpar_[4] &&
0364           zz <= gpar_[4] + gpar_[1] && r3 <= backProb_ && (depth != 2 || zz >= gpar_[4] + gpar_[0])) {
0365         double tdiff = (equalizeTimeShift_) ? (fibre_.tShift(lpos, depth, -1)) : (fibre_.tShift(lpos, depth, 1));
0366         oneHit.position = pos;
0367         oneHit.depth = depth;
0368         oneHit.time = (tSlice + (photon.t()) + tdiff);
0369         hit.push_back(oneHit);
0370 
0371 #ifdef EDM_ML_DEBUG
0372         edm::LogVerbatim("HFShower") << "HFShowerLibrary: Final Hit " << nHit << " position " << (hit[nHit].position)
0373                                      << " Depth " << (hit[nHit].depth) << " Time " << tSlice << ":" << photon.t() << ":"
0374                                      << tdiff << ":" << (hit[nHit].time);
0375 #endif
0376         ++nHit;
0377       }
0378 #ifdef EDM_ML_DEBUG
0379       else
0380         edm::LogVerbatim("HFShower") << "HFShowerLibrary: REJECTED !!!";
0381 #endif
0382       if (onlyLong && zz >= gpar_[4] + gpar_[0] && zz <= gpar_[4] + gpar_[1]) {
0383         r1 = G4UniformRand();
0384         r2 = G4UniformRand();
0385         if (rInside(r) && r1 <= exp(-p * zv) && r2 <= probMax_ && dfir > gpar_[5]) {
0386           double tdiff = (equalizeTimeShift_) ? (fibre_.tShift(lpos, 2, -1)) : (fibre_.tShift(lpos, 2, 1));
0387           oneHit.position = pos;
0388           oneHit.depth = 2;
0389           oneHit.time = (tSlice + (photon.t()) + tdiff);
0390           hit.push_back(oneHit);
0391 #ifdef EDM_ML_DEBUG
0392           edm::LogVerbatim("HFShower") << "HFShowerLibrary: Final Hit " << nHit << " position " << (hit[nHit].position)
0393                                        << " Depth " << (hit[nHit].depth) << " Time " << tSlice << ":" << photon.t()
0394                                        << ":" << tdiff << ":" << (hit[nHit].time);
0395 #endif
0396           ++nHit;
0397         }
0398       }
0399     }
0400   }
0401 
0402 #ifdef EDM_ML_DEBUG
0403   edm::LogVerbatim("HFShower") << "HFShowerLibrary: Total Hits " << nHit << " out of " << pe.size() << " PE";
0404 #endif
0405   if (nHit > pe.size() && !onlyLong) {
0406     edm::LogWarning("HFShower") << "HFShowerLibrary: Hit buffer " << pe.size() << " smaller than " << nHit << " Hits";
0407   }
0408   return hit;
0409 }
0410 
0411 bool HFShowerLibrary::rInside(double r) const { return (r >= rMin_ && r <= rMax_); }
0412 
0413 HFShowerLibrary::BranchCache::BranchCache(HFShowerLibrary::BranchReader& iReader, size_t maxRecordsToCache) {
0414   auto nRecords = iReader.numberOfRecords();
0415   if (nRecords > maxRecordsToCache) {
0416     nRecords = maxRecordsToCache;
0417   }
0418   offsets_.reserve(nRecords + 1);
0419 
0420   //first pass is to figure out how much space will be needed
0421   size_t nPhotons = 0;
0422   for (size_t r = 0; r < nRecords; ++r) {
0423     auto shower = iReader.getRecord(r + 1);
0424     nPhotons += shower.size();
0425   }
0426 
0427   photons_.reserve(nPhotons);
0428 
0429   size_t offset = 0;
0430   for (size_t r = 0; r < nRecords; ++r) {
0431     offsets_.emplace_back(offset);
0432     auto shower = iReader.getRecord(r + 1);
0433     offset += shower.size();
0434     std::copy(shower.begin(), shower.end(), std::back_inserter(photons_));
0435   }
0436   offsets_.emplace_back(offset);
0437   photons_.shrink_to_fit();
0438 }
0439 
0440 void HFShowerLibrary::BranchReader::doCaching(size_t maxRecordsToCache) {
0441   cache_ = makeCache(*this, maxRecordsToCache, branch_->GetDirectory()->GetFile()->GetName(), branch_->GetName());
0442 }
0443 
0444 std::shared_ptr<HFShowerLibrary::BranchCache const> HFShowerLibrary::BranchReader::makeCache(
0445     BranchReader& iReader, size_t maxRecordsToCache, std::string const& iFileName, std::string const& iBranchName) {
0446   //This allows sharing of the same cached data across the different modules (e.g. the per stream instances of OscarMTProducer
0447   CMS_SA_ALLOW static std::mutex s_mutex;
0448   std::lock_guard<std::mutex> guard(s_mutex);
0449   CMS_SA_ALLOW static std::map<std::pair<std::string, std::string>, std::weak_ptr<BranchCache>> s_map;
0450   auto v = s_map[{iFileName, iBranchName}].lock();
0451   if (v) {
0452     return v;
0453   }
0454   v = std::make_shared<BranchCache>(iReader, maxRecordsToCache);
0455   s_map[{iFileName, iBranchName}] = v;
0456   return v;
0457 }
0458 
0459 HFShowerPhotonCollection HFShowerLibrary::BranchCache::getRecord(int iRecord) const {
0460   assert(iRecord > 0);
0461   assert(static_cast<size_t>(iRecord + 1) < offsets_.size());
0462 
0463   auto start = offsets_[iRecord - 1];
0464   auto end = offsets_[iRecord];
0465 
0466   return HFShowerPhotonCollection(photons_.begin() + start, photons_.begin() + end);
0467 }
0468 
0469 HFShowerPhotonCollection HFShowerLibrary::BranchReader::getRecordOldForm(TBranch* iBranch, int iEntry) {
0470   HFShowerPhotonCollection photo;
0471   iBranch->SetAddress(&photo);
0472   iBranch->GetEntry(iEntry);
0473   return photo;
0474 }
0475 HFShowerPhotonCollection HFShowerLibrary::BranchReader::getRecordNewForm(TBranch* iBranch, int iEntry) {
0476   HFShowerPhotonCollection photo;
0477 
0478   auto temp = std::make_unique<HFShowerPhotonCollection>();
0479   iBranch->SetAddress(&temp);
0480   iBranch->GetEntry(iEntry);
0481   photo = std::move(*temp);
0482 
0483   return photo;
0484 }
0485 HFShowerPhotonCollection HFShowerLibrary::BranchReader::getRecordNewFormV3(TBranch* iBranch, int iEntry) {
0486   HFShowerPhotonCollection photo;
0487 
0488   std::vector<float> t;
0489   std::vector<float>* tp = &t;
0490   iBranch->SetAddress(&tp);
0491   iBranch->GetEntry(iEntry);
0492   unsigned int tSize = t.size() / 5;
0493   photo.reserve(tSize);
0494   for (unsigned int i = 0; i < tSize; i++) {
0495     photo.emplace_back(t[i], t[1 * tSize + i], t[2 * tSize + i], t[3 * tSize + i], t[4 * tSize + i]);
0496   }
0497 
0498   return photo;
0499 }
0500 
0501 HFShowerPhotonCollection HFShowerLibrary::BranchReader::getRecord(int record) const {
0502   if (cache_) {
0503     return cache_->getRecord(record);
0504   }
0505   int nrc = record - 1;
0506   HFShowerPhotonCollection photo;
0507 
0508   switch (format_) {
0509     case FileFormat::kNew: {
0510       photo = getRecordNewForm(branch_, nrc + offset_);
0511       break;
0512     }
0513     case FileFormat::kNewV3: {
0514       photo = getRecordNewFormV3(branch_, nrc + offset_);
0515       break;
0516     }
0517     case FileFormat::kOld: {
0518       photo = getRecordOldForm(branch_, nrc);
0519       break;
0520     }
0521   }
0522   return photo;
0523 }
0524 
0525 size_t HFShowerLibrary::BranchReader::numberOfRecords() const { return branch_->GetEntries() - offset_; }
0526 
0527 HFShowerPhotonCollection HFShowerLibrary::getRecord(int type, int record) const {
0528   HFShowerPhotonCollection photo;
0529   if (type > 0) {
0530     photo = hadBranch_.getRecord(record);
0531   } else {
0532     photo = emBranch_.getRecord(record);
0533   }
0534 #ifdef EDM_ML_DEBUG
0535   int nPhoton = photo.size();
0536   edm::LogVerbatim("HFShower") << "HFShowerLibrary::getRecord: Record " << record << " of type " << type << " with "
0537                                << nPhoton << " photons";
0538   for (int j = 0; j < nPhoton; j++)
0539     edm::LogVerbatim("HFShower") << "Photon " << j << " " << photo[j];
0540 #endif
0541   return photo;
0542 }
0543 
0544 HFShowerLibrary::VersionInfo HFShowerLibrary::loadEventInfo(TBranch* branch, int fileVersion) {
0545   VersionInfo versionInfo;
0546   if (branch) {
0547     std::vector<HFShowerLibraryEventInfo> eventInfoCollection;
0548     branch->SetAddress(&eventInfoCollection);
0549     branch->GetEntry(0);
0550     edm::LogVerbatim("HFShower") << "HFShowerLibrary::loadEventInfo loads EventInfo Collection of size "
0551                                  << eventInfoCollection.size() << " records";
0552 
0553     totEvents_ = eventInfoCollection[0].totalEvents();
0554     nMomBin_ = eventInfoCollection[0].numberOfBins();
0555     evtPerBin_ = eventInfoCollection[0].eventsPerBin();
0556     versionInfo.libVers_ = eventInfoCollection[0].showerLibraryVersion();
0557     versionInfo.listVersion_ = eventInfoCollection[0].physListVersion();
0558     pmom_ = eventInfoCollection[0].energyBins();
0559   } else {
0560     edm::LogVerbatim("HFShower") << "HFShowerLibrary::loadEventInfo loads EventInfo from hardwired"
0561                                  << " numbers";
0562 
0563     nMomBin_ = 16;
0564     evtPerBin_ = (fileVersion == 0) ? 5000 : 10000;
0565     totEvents_ = nMomBin_ * evtPerBin_;
0566     versionInfo.libVers_ = (fileVersion == 0) ? 1.1 : 1.2;
0567     versionInfo.listVersion_ = 3.6;
0568     pmom_ = {2, 3, 5, 7, 10, 15, 20, 30, 50, 75, 100, 150, 250, 350, 500, 1000};
0569   }
0570   for (int i = 0; i < nMomBin_; i++)
0571     pmom_[i] *= CLHEP::GeV;
0572   return versionInfo;
0573 }
0574 
0575 HFShowerPhotonCollection HFShowerLibrary::interpolate(int type, double pin) {
0576 #ifdef EDM_ML_DEBUG
0577   edm::LogVerbatim("HFShower") << "HFShowerLibrary:: Interpolate for Energy " << pin / CLHEP::GeV << " GeV with "
0578                                << nMomBin_ << " momentum bins and " << evtPerBin_ << " entries/bin -- total "
0579                                << totEvents_;
0580 #endif
0581   int irc[2] = {0, 0};
0582   double w = 0.;
0583   double r = G4UniformRand();
0584 
0585   if (pin < pmom_[0]) {
0586     w = pin / pmom_[0];
0587     irc[1] = int(evtPerBin_ * r) + 1;
0588     irc[0] = 0;
0589   } else {
0590     for (int j = 0; j < nMomBin_ - 1; j++) {
0591       if (pin >= pmom_[j] && pin < pmom_[j + 1]) {
0592         w = (pin - pmom_[j]) / (pmom_[j + 1] - pmom_[j]);
0593         if (j == nMomBin_ - 2) {
0594           irc[1] = int(evtPerBin_ * 0.5 * r);
0595         } else {
0596           irc[1] = int(evtPerBin_ * r);
0597         }
0598         irc[1] += (j + 1) * evtPerBin_ + 1;
0599         r = G4UniformRand();
0600         irc[0] = int(evtPerBin_ * r) + 1 + j * evtPerBin_;
0601         if (irc[0] < 0) {
0602           edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[0] = " << irc[0] << " now set to 0";
0603           irc[0] = 0;
0604         } else if (irc[0] > totEvents_) {
0605           edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[0] = " << irc[0] << " now set to "
0606                                       << totEvents_;
0607           irc[0] = totEvents_;
0608         }
0609       }
0610     }
0611   }
0612   if (irc[1] < 1) {
0613     edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[1] = " << irc[1] << " now set to 1";
0614     irc[1] = 1;
0615   } else if (irc[1] > totEvents_) {
0616     edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[1] = " << irc[1] << " now set to " << totEvents_;
0617     irc[1] = totEvents_;
0618   }
0619 
0620 #ifdef EDM_ML_DEBUG
0621   edm::LogVerbatim("HFShower") << "HFShowerLibrary:: Select records " << irc[0] << " and " << irc[1] << " with weights "
0622                                << 1 - w << " and " << w;
0623 #endif
0624   HFShowerPhotonCollection pe;
0625 
0626   std::size_t npold = 0;
0627   for (int ir = 0; ir < 2; ir++) {
0628     if (irc[ir] > 0) {
0629       auto photons = getRecord(type, irc[ir]);
0630       int nPhoton = photons.size();
0631       npold += nPhoton;
0632       pe.reserve(pe.size() + nPhoton);
0633       for (auto const& photon : photons) {
0634         r = G4UniformRand();
0635         if ((ir == 0 && r > w) || (ir > 0 && r < w)) {
0636           storePhoton(photon, pe);
0637         }
0638       }
0639     }
0640   }
0641 
0642   if ((pe.size() > npold || (npold == 0 && irc[0] > 0)) && !(pe.empty() && npold == 0))
0643     edm::LogWarning("HFShower") << "HFShowerLibrary: Interpolation Warning =="
0644                                 << " records " << irc[0] << " and " << irc[1] << " gives a buffer of " << npold
0645                                 << " photons and fills " << pe.size() << " *****";
0646 #ifdef EDM_ML_DEBUG
0647   else
0648     edm::LogVerbatim("HFShower") << "HFShowerLibrary: Interpolation == records " << irc[0] << " and " << irc[1]
0649                                  << " gives a buffer of " << npold << " photons and fills " << pe.size() << " PE";
0650   for (std::size_t j = 0; j < pe.size(); j++)
0651     edm::LogVerbatim("HFShower") << "Photon " << j << " " << pe[j];
0652 #endif
0653   return pe;
0654 }
0655 
0656 HFShowerPhotonCollection HFShowerLibrary::extrapolate(int type, double pin) {
0657   int nrec = int(pin / pmom_[nMomBin_ - 1]);
0658   double w = (pin - pmom_[nMomBin_ - 1] * nrec) / pmom_[nMomBin_ - 1];
0659   nrec++;
0660 #ifdef EDM_ML_DEBUG
0661   edm::LogVerbatim("HFShower") << "HFShowerLibrary:: Extrapolate for Energy " << pin / CLHEP::GeV << " GeV with "
0662                                << nMomBin_ << " momentum bins and " << evtPerBin_ << " entries/bin -- "
0663                                << "total " << totEvents_ << " using " << nrec << " records";
0664 #endif
0665   std::vector<int> irc(nrec);
0666 
0667   for (int ir = 0; ir < nrec; ir++) {
0668     double r = G4UniformRand();
0669     irc[ir] = int(evtPerBin_ * 0.5 * r) + (nMomBin_ - 1) * evtPerBin_ + 1;
0670     if (irc[ir] < 1) {
0671       edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[" << ir << "] = " << irc[ir] << " now set to 1";
0672       irc[ir] = 1;
0673     } else if (irc[ir] > totEvents_) {
0674       edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[" << ir << "] = " << irc[ir] << " now set to "
0675                                   << totEvents_;
0676       irc[ir] = totEvents_;
0677 #ifdef EDM_ML_DEBUG
0678     } else {
0679       edm::LogVerbatim("HFShower") << "HFShowerLibrary::Extrapolation use irc[" << ir << "] = " << irc[ir];
0680 #endif
0681     }
0682   }
0683 
0684   HFShowerPhotonCollection pe;
0685   std::size_t npold = 0;
0686   for (int ir = 0; ir < nrec; ir++) {
0687     if (irc[ir] > 0) {
0688       auto const photons = getRecord(type, irc[ir]);
0689       int nPhoton = photons.size();
0690       npold += nPhoton;
0691       pe.reserve(pe.size() + nPhoton);
0692       for (auto const& photon : photons) {
0693         double r = G4UniformRand();
0694         if (ir != nrec - 1 || r < w) {
0695           storePhoton(photon, pe);
0696         }
0697       }
0698 #ifdef EDM_ML_DEBUG
0699       edm::LogVerbatim("HFShower") << "HFShowerLibrary: Record [" << ir << "] = " << irc[ir] << " npold = " << npold;
0700 #endif
0701     }
0702   }
0703 #ifdef EDM_ML_DEBUG
0704   edm::LogVerbatim("HFShower") << "HFShowerLibrary:: uses " << npold << " photons";
0705 #endif
0706 
0707   if (pe.size() > npold || npold == 0)
0708     edm::LogWarning("HFShower") << "HFShowerLibrary: Extrapolation Warning == " << nrec << " records " << irc[0] << ", "
0709                                 << irc[1] << ", ... gives a buffer of " << npold << " photons and fills " << pe.size()
0710                                 << " *****";
0711 #ifdef EDM_ML_DEBUG
0712   else
0713     edm::LogVerbatim("HFShower") << "HFShowerLibrary: Extrapolation == " << nrec << " records " << irc[0] << ", "
0714                                  << irc[1] << ", ... gives a buffer of " << npold << " photons and fills " << pe.size()
0715                                  << " PE";
0716   for (std::size_t j = 0; j < pe.size(); j++)
0717     edm::LogVerbatim("HFShower") << "Photon " << j << " " << pe[j];
0718 #endif
0719   return pe;
0720 }
0721 
0722 void HFShowerLibrary::storePhoton(HFShowerPhoton const& iPhoton, HFShowerPhotonCollection& iPhotons) const {
0723   iPhotons.push_back(iPhoton);
0724 #ifdef EDM_ML_DEBUG
0725   edm::LogVerbatim("HFShower") << "HFShowerLibrary: storePhoton " << iPhoton << " npe " << iPhotons.size() << " "
0726                                << iPhotons.back();
0727 #endif
0728 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!