CMS 3D CMS Logo

List of all members | Classes | Public Member Functions | Protected Member Functions | Private Attributes
HFShowerLibrary Class Reference

#include <HFShowerLibrary.h>

Classes

struct  Hit
 

Public Member Functions

std::vector< HitfillHits (const G4ThreeVector &p, const G4ThreeVector &v, int parCode, double parEnergy, bool &ok, double weight, double time, bool onlyLong=false)
 
std::vector< HitgetHits (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
 
bool equalizeTimeShift_
 
int evtPerBin
 
std::unique_ptr< HFFibrefibre_
 
int fileVersion_
 
std::vector< double > gpar
 
TBranch * hadBranch
 
const HcalDDDSimConstantshcalConstant_
 
TFile * hf
 
float libVers
 
float listVersion
 
bool newForm
 
int nMomBin
 
int npe
 
HFShowerPhotonCollection pe
 
std::unique_ptr< HFShowerPhotonCollectionphoto
 
HFShowerPhotonCollection photon
 
std::vector< double > pmom
 
double probMax
 
double rMax
 
double rMin
 
int totEvents
 
bool v3version
 
bool verbose
 

Detailed Description

Definition at line 28 of file HFShowerLibrary.h.

Constructor & Destructor Documentation

◆ HFShowerLibrary()

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, equalizeTimeShift_, edmPickEvents::event, evtPerBin, Exception, fibre_, fileVersion_, personalPlayback::fp, HcalDDDSimConstants::getGparHF(), edm::ParameterSet::getParameter(), HcalDDDSimConstants::getPhiTableHF(), HcalDDDSimConstants::getRTableHF(), 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, v3version, and verbose.

27  : hcalConstant_(hcons), hf(nullptr), emBranch(nullptr), hadBranch(nullptr), npe(0) {
28  edm::ParameterSet m_HF =
29  (p.getParameter<edm::ParameterSet>("HFShower")).getParameter<edm::ParameterSet>("HFShowerBlock");
30  probMax = m_HF.getParameter<double>("ProbMax");
31  equalizeTimeShift_ = m_HF.getParameter<bool>("EqualizeTimeShift");
32 
33  edm::ParameterSet m_HS =
34  (p.getParameter<edm::ParameterSet>("HFShowerLibrary")).getParameter<edm::ParameterSet>("HFLibraryFileBlock");
35  edm::FileInPath fp = m_HS.getParameter<edm::FileInPath>("FileName");
36  std::string pTreeName = fp.fullPath();
37  backProb = m_HS.getParameter<double>("BackProbability");
38  std::string emName = m_HS.getParameter<std::string>("TreeEMID");
39  std::string hadName = m_HS.getParameter<std::string>("TreeHadID");
40  std::string branchEvInfo = m_HS.getUntrackedParameter<std::string>(
41  "BranchEvt", "HFShowerLibraryEventInfos_hfshowerlib_HFShowerLibraryEventInfo");
42  std::string branchPre = m_HS.getUntrackedParameter<std::string>("BranchPre", "HFShowerPhotons_hfshowerlib_");
43  std::string branchPost = m_HS.getUntrackedParameter<std::string>("BranchPost", "_R.obj");
44  verbose = m_HS.getUntrackedParameter<bool>("Verbosity", false);
45  applyFidCut = m_HS.getParameter<bool>("ApplyFiducialCut");
46  fileVersion_ = m_HS.getParameter<int>("FileVersion");
47 
48  if (pTreeName.find('.') == 0)
49  pTreeName.erase(0, 2);
50  const char* nTree = pTreeName.c_str();
51  hf = TFile::Open(nTree);
52 
53  if (!hf->IsOpen()) {
54  edm::LogError("HFShower") << "HFShowerLibrary: opening " << nTree << " failed";
55  throw cms::Exception("Unknown", "HFShowerLibrary") << "Opening of " << pTreeName << " fails\n";
56  } else {
57  edm::LogVerbatim("HFShower") << "HFShowerLibrary: opening " << nTree << " successfully";
58  }
59 
60  newForm = (branchEvInfo.empty());
61  TTree* event(nullptr);
62  if (newForm)
63  event = (TTree*)hf->Get("HFSimHits");
64  else
65  event = (TTree*)hf->Get("Events");
66  if (event) {
67  TBranch* evtInfo(nullptr);
68  if (!newForm) {
69  std::string info = branchEvInfo + branchPost;
70  evtInfo = event->GetBranch(info.c_str());
71  }
72  if (evtInfo || newForm) {
73  loadEventInfo(evtInfo);
74  } else {
75  edm::LogError("HFShower") << "HFShowerLibrary: HFShowerLibrayEventInfo"
76  << " Branch does not exist in Event";
77  throw cms::Exception("Unknown", "HFShowerLibrary") << "Event information absent\n";
78  }
79  } else {
80  edm::LogError("HFShower") << "HFShowerLibrary: Events Tree does not "
81  << "exist";
82  throw cms::Exception("Unknown", "HFShowerLibrary") << "Events tree absent\n";
83  }
84 
85  std::stringstream ss;
86  ss << "HFShowerLibrary: Library " << libVers << " ListVersion " << listVersion << " File version " << fileVersion_
87  << " Events Total " << totEvents << " and " << evtPerBin << " per bin\n";
88  ss << "HFShowerLibrary: Energies (GeV) with " << nMomBin << " bins\n";
89  for (int i = 0; i < nMomBin; ++i) {
90  if (i / 10 * 10 == i && i > 0) {
91  ss << "\n";
92  }
93  ss << " " << pmom[i] / CLHEP::GeV;
94  }
95  edm::LogVerbatim("HFShower") << ss.str();
96 
97  std::string nameBr = branchPre + emName + branchPost;
98  emBranch = event->GetBranch(nameBr.c_str());
99  if (verbose)
100  emBranch->Print();
101  nameBr = branchPre + hadName + branchPost;
102  hadBranch = event->GetBranch(nameBr.c_str());
103  if (verbose)
104  hadBranch->Print();
105 
106  v3version = false;
107  if (emBranch->GetClassName() == std::string("vector<float>")) {
108  v3version = true;
109  }
110 
111  edm::LogVerbatim("HFShower") << " HFShowerLibrary:Branch " << emName << " has " << emBranch->GetEntries()
112  << " entries and Branch " << hadName << " has " << hadBranch->GetEntries()
113  << " entries\n HFShowerLibrary::No packing information - Assume x, y, z are not in "
114  "packed form\n Maximum probability cut off "
115  << probMax << " Back propagation of light probability " << backProb
116  << " Flag for equalizing Time Shift for different eta " << equalizeTimeShift_;
117 
118  fibre_ = std::make_unique<HFFibre>(name, hcalConstant_, hps, p);
119  photo = std::make_unique<HFShowerPhotonCollection>();
120 
121  //Radius (minimum and maximum)
122  std::vector<double> rTable = hcalConstant_->getRTableHF();
123  rMin = rTable[0];
124  rMax = rTable[rTable.size() - 1];
125 
126  //Delta phi
127  std::vector<double> phibin = hcalConstant_->getPhiTableHF();
128  dphi = phibin[0];
129 
130  edm::LogVerbatim("HFShower") << "HFShowerLibrary: rMIN " << rMin / CLHEP::cm << " cm and rMax " << rMax / CLHEP::cm
131  << " (Half) Phi Width of wedge " << dphi / CLHEP::deg;
132 
133  //Special Geometry parameters
135 }
Log< level::Info, true > LogVerbatim
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
static const TGPicture * info(bool iBackgroundIsBlack)
const HcalDDDSimConstants * hcalConstant_
std::unique_ptr< HFFibre > fibre_
Log< level::Error, false > LogError
const std::vector< double > & getGparHF() const
const std::vector< double > & getRTableHF() const
void loadEventInfo(TBranch *)
std::vector< double > gpar
const std::vector< double > & getPhiTableHF() const
TBranch * emBranch
std::vector< double > pmom
std::unique_ptr< HFShowerPhotonCollection > photo
Definition: event.py:1
TBranch * hadBranch

◆ ~HFShowerLibrary()

HFShowerLibrary::~HFShowerLibrary ( )

Definition at line 137 of file HFShowerLibrary.cc.

References hf.

137  {
138  if (hf)
139  hf->Close();
140 }

Member Function Documentation

◆ extrapolate()

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

Definition at line 515 of file HFShowerLibrary.cc.

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

Referenced by fillHits().

515  {
516  int nrec = int(pin / pmom[nMomBin - 1]);
517  double w = (pin - pmom[nMomBin - 1] * nrec) / pmom[nMomBin - 1];
518  nrec++;
519 #ifdef EDM_ML_DEBUG
520  edm::LogVerbatim("HFShower") << "HFShowerLibrary:: Extrapolate for Energy " << pin / CLHEP::GeV << " GeV with "
521  << nMomBin << " momentum bins and " << evtPerBin << " entries/bin -- "
522  << "total " << totEvents << " using " << nrec << " records";
523 #endif
524  std::vector<int> irc(nrec);
525 
526  for (int ir = 0; ir < nrec; ir++) {
527  double r = G4UniformRand();
528  irc[ir] = int(evtPerBin * 0.5 * r) + (nMomBin - 1) * evtPerBin + 1;
529  if (irc[ir] < 1) {
530  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[" << ir << "] = " << irc[ir] << " now set to 1";
531  irc[ir] = 1;
532  } else if (irc[ir] > totEvents) {
533  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[" << ir << "] = " << irc[ir] << " now set to "
534  << totEvents;
535  irc[ir] = totEvents;
536 #ifdef EDM_ML_DEBUG
537  } else {
538  edm::LogVerbatim("HFShower") << "HFShowerLibrary::Extrapolation use irc[" << ir << "] = " << irc[ir];
539 #endif
540  }
541  }
542 
543  pe.clear();
544  npe = 0;
545  int npold = 0;
546  for (int ir = 0; ir < nrec; ir++) {
547  if (irc[ir] > 0) {
548  getRecord(type, irc[ir]);
549  int nPhoton = (newForm) ? photo->size() : photon.size();
550  npold += nPhoton;
551  for (int j = 0; j < nPhoton; j++) {
552  double r = G4UniformRand();
553  if (ir != nrec - 1 || r < w) {
554  storePhoton(j);
555  }
556  }
557 #ifdef EDM_ML_DEBUG
558  edm::LogVerbatim("HFShower") << "HFShowerLibrary: Record [" << ir << "] = " << irc[ir] << " npold = " << npold;
559 #endif
560  }
561  }
562 #ifdef EDM_ML_DEBUG
563  edm::LogVerbatim("HFShower") << "HFShowerLibrary:: uses " << npold << " photons";
564 #endif
565 
566  if (npe > npold || npold == 0)
567  edm::LogWarning("HFShower") << "HFShowerLibrary: Extrapolation Warning == " << nrec << " records " << irc[0] << ", "
568  << irc[1] << ", ... gives a buffer of " << npold << " photons and fills " << npe
569  << " *****";
570 #ifdef EDM_ML_DEBUG
571  else
572  edm::LogVerbatim("HFShower") << "HFShowerLibrary: Extrapolation == " << nrec << " records " << irc[0] << ", "
573  << irc[1] << ", ... gives a buffer of " << npold << " photons and fills " << npe
574  << " PE";
575  for (int j = 0; j < npe; j++)
576  edm::LogVerbatim("HFShower") << "Photon " << j << " " << pe[j];
577 #endif
578 }
Log< level::Info, true > LogVerbatim
void getRecord(int, int)
T w() const
HFShowerPhotonCollection pe
std::vector< double > pmom
std::unique_ptr< HFShowerPhotonCollection > photo
void storePhoton(int j)
Log< level::Warning, false > LogWarning
HFShowerPhotonCollection photon

◆ fillHits()

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 183 of file HFShowerLibrary.cc.

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

Referenced by getHits().

190  {
191  std::vector<HFShowerLibrary::Hit> hit;
192  ok = false;
194  // shower is built only for gamma, e+- and stable hadrons
195  if (!isEM && !G4TrackToParticleID::isStableHadron(parCode)) {
196  return hit;
197  }
198  ok = true;
199 
200  // remove low-energy component
201  const double threshold = 50 * MeV;
202  if (pin < threshold) {
203  return hit;
204  }
205 
206  double pz = momDir.z();
207  double zint = hitPoint.z();
208 
209  // if particle moves from interaction point or "backwards (halo)
210  bool backward = (pz * zint < 0.) ? true : false;
211 
212  double sphi = sin(momDir.phi());
213  double cphi = cos(momDir.phi());
214  double ctheta = cos(momDir.theta());
215  double stheta = sin(momDir.theta());
216 
217  if (isEM) {
218  if (pin < pmom[nMomBin - 1]) {
219  interpolate(0, pin);
220  } else {
221  extrapolate(0, pin);
222  }
223  } else {
224  if (pin < pmom[nMomBin - 1]) {
225  interpolate(1, pin);
226  } else {
227  extrapolate(1, pin);
228  }
229  }
230 
231  int nHit = 0;
232  HFShowerLibrary::Hit oneHit;
233  for (int i = 0; i < npe; ++i) {
234  double zv = std::abs(pe[i].z()); // abs local z
235 #ifdef EDM_ML_DEBUG
236  edm::LogVerbatim("HFShower") << "HFShowerLibrary: Hit " << i << " " << pe[i] << " zv " << zv;
237 #endif
238  if (zv <= gpar[1] && pe[i].lambda() > 0 && (pe[i].z() >= 0 || (zv > gpar[0] && (!onlyLong)))) {
239  int depth = 1;
240  if (onlyLong) {
241  } else if (!backward) { // fully valid only for "front" particles
242  if (pe[i].z() < 0)
243  depth = 2; // with "front"-simulated shower lib.
244  } else { // for "backward" particles - almost equal
245  double r = G4UniformRand(); // share between L and S fibers
246  if (r > 0.5)
247  depth = 2;
248  }
249 
250  // Updated coordinate transformation from local
251  // back to global using two Euler angles: phi and theta
252  double pex = pe[i].x();
253  double pey = pe[i].y();
254 
255  double xx = pex * ctheta * cphi - pey * sphi + zv * stheta * cphi;
256  double yy = pex * ctheta * sphi + pey * cphi + zv * stheta * sphi;
257  double zz = -pex * stheta + zv * ctheta;
258 
259  G4ThreeVector pos = hitPoint + G4ThreeVector(xx, yy, zz);
260  zv = std::abs(pos.z()) - gpar[4] - 0.5 * gpar[1];
261  G4ThreeVector lpos = G4ThreeVector(pos.x(), pos.y(), zv);
262 
263  zv = fibre_->zShift(lpos, depth, 0); // distance to PMT !
264 
265  double r = pos.perp();
266  double p = fibre_->attLength(pe[i].lambda());
267  double fi = pos.phi();
268  if (fi < 0)
269  fi += CLHEP::twopi;
270  int isect = int(fi / dphi) + 1;
271  isect = (isect + 1) / 2;
272  double dfi = ((isect * 2 - 1) * dphi - fi);
273  if (dfi < 0)
274  dfi = -dfi;
275  double dfir = r * sin(dfi);
276 #ifdef EDM_ML_DEBUG
277  edm::LogVerbatim("HFShower") << "HFShowerLibrary: Position shift " << xx << ", " << yy << ", " << zz << ": "
278  << pos << " R " << r << " Phi " << fi << " Section " << isect << " R*Dfi " << dfir
279  << " Dist " << zv;
280 #endif
281  zz = std::abs(pos.z());
282  double r1 = G4UniformRand();
283  double r2 = G4UniformRand();
284  double r3 = backward ? G4UniformRand() : -9999.;
285  if (!applyFidCut)
286  dfir += gpar[5];
287 
288 #ifdef EDM_ML_DEBUG
289  edm::LogVerbatim("HFShower") << "HFShowerLibrary: rLimits " << rInside(r) << " attenuation " << r1 << ":"
290  << exp(-p * zv) << " r2 " << r2 << " r3 " << r3 << " rDfi " << gpar[5] << " zz "
291  << zz << " zLim " << gpar[4] << ":" << gpar[4] + gpar[1] << "\n"
292  << " rInside(r) :" << rInside(r) << " r1 <= exp(-p*zv) :" << (r1 <= exp(-p * zv))
293  << " r2 <= probMax :" << (r2 <= probMax * weight)
294  << " r3 <= backProb :" << (r3 <= backProb)
295  << " dfir > gpar[5] :" << (dfir > gpar[5]) << " zz >= gpar[4] :" << (zz >= gpar[4])
296  << " zz <= gpar[4]+gpar[1] :" << (zz <= gpar[4] + gpar[1]);
297 #endif
298  if (rInside(r) && r1 <= exp(-p * zv) && r2 <= probMax * weight && dfir > gpar[5] && zz >= gpar[4] &&
299  zz <= gpar[4] + gpar[1] && r3 <= backProb && (depth != 2 || zz >= gpar[4] + gpar[0])) {
300  double tdiff = (equalizeTimeShift_) ? (fibre_->tShift(lpos, depth, -1)) : (fibre_->tShift(lpos, depth, 1));
301  oneHit.position = pos;
302  oneHit.depth = depth;
303  oneHit.time = (tSlice + (pe[i].t()) + tdiff);
304  hit.push_back(oneHit);
305 
306 #ifdef EDM_ML_DEBUG
307  edm::LogVerbatim("HFShower") << "HFShowerLibrary: Final Hit " << nHit << " position " << (hit[nHit].position)
308  << " Depth " << (hit[nHit].depth) << " Time " << tSlice << ":" << pe[i].t() << ":"
309  << tdiff << ":" << (hit[nHit].time);
310 #endif
311  ++nHit;
312  }
313 #ifdef EDM_ML_DEBUG
314  else
315  edm::LogVerbatim("HFShower") << "HFShowerLibrary: REJECTED !!!";
316 #endif
317  if (onlyLong && zz >= gpar[4] + gpar[0] && zz <= gpar[4] + gpar[1]) {
318  r1 = G4UniformRand();
319  r2 = G4UniformRand();
320  if (rInside(r) && r1 <= exp(-p * zv) && r2 <= probMax && dfir > gpar[5]) {
321  double tdiff = (equalizeTimeShift_) ? (fibre_->tShift(lpos, 2, -1)) : (fibre_->tShift(lpos, 2, 1));
322  oneHit.position = pos;
323  oneHit.depth = 2;
324  oneHit.time = (tSlice + (pe[i].t()) + tdiff);
325  hit.push_back(oneHit);
326 #ifdef EDM_ML_DEBUG
327  edm::LogVerbatim("HFShower") << "HFShowerLibrary: Final Hit " << nHit << " position " << (hit[nHit].position)
328  << " Depth " << (hit[nHit].depth) << " Time " << tSlice << ":" << pe[i].t()
329  << ":" << tdiff << ":" << (hit[nHit].time);
330 #endif
331  ++nHit;
332  }
333  }
334  }
335  }
336 
337 #ifdef EDM_ML_DEBUG
338  edm::LogVerbatim("HFShower") << "HFShowerLibrary: Total Hits " << nHit << " out of " << npe << " PE";
339 #endif
340  if (nHit > npe && !onlyLong) {
341  edm::LogWarning("HFShower") << "HFShowerLibrary: Hit buffer " << npe << " smaller than " << nHit << " Hits";
342  }
343  return hit;
344 }
Log< level::Info, true > LogVerbatim
std::unique_ptr< HFFibre > fibre_
float *__restrict__ zv
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
HFShowerPhotonCollection pe
Definition: weight.py:1
static bool isStableHadron(int pdgCode)
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< double > gpar
std::vector< double > pmom
G4ThreeVector position
void extrapolate(int, double)
static bool isGammaElectronPositron(int pdgCode)
Log< level::Warning, false > LogWarning
bool rInside(double r)
void interpolate(int, double)

◆ getHits()

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

Definition at line 142 of file HFShowerLibrary.cc.

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

145  {
146  auto const preStepPoint = aStep->GetPreStepPoint();
147  auto const postStepPoint = aStep->GetPostStepPoint();
148  auto const track = aStep->GetTrack();
149  // Get Z-direction
150  auto const aParticle = track->GetDynamicParticle();
151  const G4ThreeVector& momDir = aParticle->GetMomentumDirection();
152  const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
153  int parCode = track->GetDefinition()->GetPDGEncoding();
154 
155  // VI: for ions use internally pdg code of alpha in order to keep
156  // consistency with previous simulation
157  if (track->GetDefinition()->IsGeneralIon()) {
158  parCode = 1000020040;
159  }
160 
161 #ifdef EDM_ML_DEBUG
162  G4String partType = track->GetDefinition()->GetParticleName();
163  const G4ThreeVector localPos = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
164  double zoff = localPos.z() + 0.5 * gpar[1];
165 
166  edm::LogVerbatim("HFShower") << "HFShowerLibrary::getHits " << partType << " of energy "
167  << track->GetKineticEnergy() / CLHEP::GeV << " GeV weight= " << weight
168  << " onlyLong: " << onlyLong << " dir.orts " << momDir.x() << ", " << momDir.y() << ", "
169  << momDir.z() << " Pos x,y,z = " << hitPoint.x() << "," << hitPoint.y() << ","
170  << hitPoint.z() << " (" << zoff << ") sphi,cphi,stheta,ctheta = " << sin(momDir.phi())
171  << "," << cos(momDir.phi()) << ", " << sin(momDir.theta()) << "," << cos(momDir.theta());
172 #endif
173 
174  double tSlice = (postStepPoint->GetGlobalTime()) / CLHEP::nanosecond;
175 
176  // use kinetic energy for protons and ions
177  double pin = (track->GetDefinition()->GetBaryonNumber() > 0) ? preStepPoint->GetKineticEnergy()
178  : preStepPoint->GetTotalEnergy();
179 
180  return fillHits(hitPoint, momDir, parCode, pin, isKilled, weight, tSlice, onlyLong);
181 }
Log< level::Info, true > LogVerbatim
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Definition: weight.py:1
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
std::vector< double > gpar
std::vector< Hit > fillHits(const G4ThreeVector &p, const G4ThreeVector &v, int parCode, double parEnergy, bool &ok, double weight, double time, bool onlyLong=false)

◆ getRecord()

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

Definition at line 348 of file HFShowerLibrary.cc.

References emBranch, fileVersion_, hadBranch, mps_fire::i, dqmiolumiharvest::j, newForm, l1ctLayer1_cff::nPhoton, photo, photon, position, AlCaHarvesting_cff::record, submitPVValidationJobs::t, totEvents, cmsswSequenceInfo::tp, and v3version.

Referenced by extrapolate(), and interpolate().

348  {
349  int nrc = record - 1;
350  photon.clear();
351  photo->clear();
352  if (type > 0) {
353  if (newForm) {
354  if (!v3version) {
355  hadBranch->SetAddress(&photo);
356  int position = (fileVersion_ >= 2) ? nrc : (nrc + totEvents);
357  hadBranch->GetEntry(position);
358  } else {
359  std::vector<float> t;
360  std::vector<float>* tp = &t;
361  hadBranch->SetAddress(&tp);
362  hadBranch->GetEntry(nrc + totEvents);
363  unsigned int tSize = t.size() / 5;
364  photo->reserve(tSize);
365  for (unsigned int i = 0; i < tSize; i++) {
366  photo->push_back(
367  HFShowerPhoton(t[i], t[1 * tSize + i], t[2 * tSize + i], t[3 * tSize + i], t[4 * tSize + i]));
368  }
369  }
370  } else {
371  hadBranch->SetAddress(&photon);
372  hadBranch->GetEntry(nrc);
373  }
374  } else {
375  if (newForm) {
376  if (!v3version) {
377  emBranch->SetAddress(&photo);
378  emBranch->GetEntry(nrc);
379  } else {
380  std::vector<float> t;
381  std::vector<float>* tp = &t;
382  emBranch->SetAddress(&tp);
383  emBranch->GetEntry(nrc);
384  unsigned int tSize = t.size() / 5;
385  photo->reserve(tSize);
386  for (unsigned int i = 0; i < tSize; i++) {
387  photo->push_back(
388  HFShowerPhoton(t[i], t[1 * tSize + i], t[2 * tSize + i], t[3 * tSize + i], t[4 * tSize + i]));
389  }
390  }
391  } else {
392  emBranch->SetAddress(&photon);
393  emBranch->GetEntry(nrc);
394  }
395  }
396 #ifdef EDM_ML_DEBUG
397  int nPhoton = (newForm) ? photo->size() : photon.size();
398  edm::LogVerbatim("HFShower") << "HFShowerLibrary::getRecord: Record " << record << " of type " << type << " with "
399  << nPhoton << " photons";
400  for (int j = 0; j < nPhoton; j++)
401  if (newForm)
402  edm::LogVerbatim("HFShower") << "Photon " << j << " " << photo->at(j);
403  else
404  edm::LogVerbatim("HFShower") << "Photon " << j << " " << photon[j];
405 #endif
406 }
Log< level::Info, true > LogVerbatim
TBranch * emBranch
std::unique_ptr< HFShowerPhotonCollection > photo
static int position[264][3]
Definition: ReadPGInfo.cc:289
HFShowerPhotonCollection photon
TBranch * hadBranch

◆ interpolate()

void HFShowerLibrary::interpolate ( int  type,
double  pin 
)
protected

Definition at line 437 of file HFShowerLibrary.cc.

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

Referenced by fillHits().

437  {
438 #ifdef EDM_ML_DEBUG
439  edm::LogVerbatim("HFShower") << "HFShowerLibrary:: Interpolate for Energy " << pin / CLHEP::GeV << " GeV with "
440  << nMomBin << " momentum bins and " << evtPerBin << " entries/bin -- total "
441  << totEvents;
442 #endif
443  int irc[2] = {0, 0};
444  double w = 0.;
445  double r = G4UniformRand();
446 
447  if (pin < pmom[0]) {
448  w = pin / pmom[0];
449  irc[1] = int(evtPerBin * r) + 1;
450  irc[0] = 0;
451  } else {
452  for (int j = 0; j < nMomBin - 1; j++) {
453  if (pin >= pmom[j] && pin < pmom[j + 1]) {
454  w = (pin - pmom[j]) / (pmom[j + 1] - pmom[j]);
455  if (j == nMomBin - 2) {
456  irc[1] = int(evtPerBin * 0.5 * r);
457  } else {
458  irc[1] = int(evtPerBin * r);
459  }
460  irc[1] += (j + 1) * evtPerBin + 1;
461  r = G4UniformRand();
462  irc[0] = int(evtPerBin * r) + 1 + j * evtPerBin;
463  if (irc[0] < 0) {
464  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[0] = " << irc[0] << " now set to 0";
465  irc[0] = 0;
466  } else if (irc[0] > totEvents) {
467  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[0] = " << irc[0] << " now set to " << totEvents;
468  irc[0] = totEvents;
469  }
470  }
471  }
472  }
473  if (irc[1] < 1) {
474  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[1] = " << irc[1] << " now set to 1";
475  irc[1] = 1;
476  } else if (irc[1] > totEvents) {
477  edm::LogWarning("HFShower") << "HFShowerLibrary:: Illegal irc[1] = " << irc[1] << " now set to " << totEvents;
478  irc[1] = totEvents;
479  }
480 
481 #ifdef EDM_ML_DEBUG
482  edm::LogVerbatim("HFShower") << "HFShowerLibrary:: Select records " << irc[0] << " and " << irc[1] << " with weights "
483  << 1 - w << " and " << w;
484 #endif
485  pe.clear();
486  npe = 0;
487  int npold = 0;
488  for (int ir = 0; ir < 2; ir++) {
489  if (irc[ir] > 0) {
490  getRecord(type, irc[ir]);
491  int nPhoton = (newForm) ? photo->size() : photon.size();
492  npold += nPhoton;
493  for (int j = 0; j < nPhoton; j++) {
494  r = G4UniformRand();
495  if ((ir == 0 && r > w) || (ir > 0 && r < w)) {
496  storePhoton(j);
497  }
498  }
499  }
500  }
501 
502  if ((npe > npold || (npold == 0 && irc[0] > 0)) && !(npe == 0 && npold == 0))
503  edm::LogWarning("HFShower") << "HFShowerLibrary: Interpolation Warning =="
504  << " records " << irc[0] << " and " << irc[1] << " gives a buffer of " << npold
505  << " photons and fills " << npe << " *****";
506 #ifdef EDM_ML_DEBUG
507  else
508  edm::LogVerbatim("HFShower") << "HFShowerLibrary: Interpolation == records " << irc[0] << " and " << irc[1]
509  << " gives a buffer of " << npold << " photons and fills " << npe << " PE";
510  for (int j = 0; j < npe; j++)
511  edm::LogVerbatim("HFShower") << "Photon " << j << " " << pe[j];
512 #endif
513 }
Log< level::Info, true > LogVerbatim
void getRecord(int, int)
T w() const
HFShowerPhotonCollection pe
std::vector< double > pmom
std::unique_ptr< HFShowerPhotonCollection > photo
void storePhoton(int j)
Log< level::Warning, false > LogWarning
HFShowerPhotonCollection photon

◆ loadEventInfo()

void HFShowerLibrary::loadEventInfo ( TBranch *  branch)
protected

Definition at line 408 of file HFShowerLibrary.cc.

References MicroEventContent_cff::branch, evtPerBin, fileVersion_, mps_fire::i, libVers, listVersion, nMomBin, pmom, and totEvents.

Referenced by HFShowerLibrary().

408  {
409  if (branch) {
410  std::vector<HFShowerLibraryEventInfo> eventInfoCollection;
411  branch->SetAddress(&eventInfoCollection);
412  branch->GetEntry(0);
413  edm::LogVerbatim("HFShower") << "HFShowerLibrary::loadEventInfo loads EventInfo Collection of size "
414  << eventInfoCollection.size() << " records";
415 
416  totEvents = eventInfoCollection[0].totalEvents();
417  nMomBin = eventInfoCollection[0].numberOfBins();
418  evtPerBin = eventInfoCollection[0].eventsPerBin();
419  libVers = eventInfoCollection[0].showerLibraryVersion();
420  listVersion = eventInfoCollection[0].physListVersion();
421  pmom = eventInfoCollection[0].energyBins();
422  } else {
423  edm::LogVerbatim("HFShower") << "HFShowerLibrary::loadEventInfo loads EventInfo from hardwired"
424  << " numbers";
425 
426  nMomBin = 16;
427  evtPerBin = (fileVersion_ == 0) ? 5000 : 10000;
429  libVers = (fileVersion_ == 0) ? 1.1 : 1.2;
430  listVersion = 3.6;
431  pmom = {2, 3, 5, 7, 10, 15, 20, 30, 50, 75, 100, 150, 250, 350, 500, 1000};
432  }
433  for (int i = 0; i < nMomBin; i++)
434  pmom[i] *= CLHEP::GeV;
435 }
Log< level::Info, true > LogVerbatim
std::vector< double > pmom

◆ rInside()

bool HFShowerLibrary::rInside ( double  r)
protected

Definition at line 346 of file HFShowerLibrary.cc.

References alignCSCRings::r, rMax, and rMin.

Referenced by fillHits().

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

◆ storePhoton()

void HFShowerLibrary::storePhoton ( int  j)
protected

Definition at line 580 of file HFShowerLibrary.cc.

References dqmiolumiharvest::j, newForm, npe, pe, photo, and photon.

Referenced by extrapolate(), and interpolate().

580  {
581  if (newForm)
582  pe.push_back(photo->at(j));
583  else
584  pe.push_back(photon[j]);
585  npe++;
586 #ifdef EDM_ML_DEBUG
587  edm::LogVerbatim("HFShower") << "HFShowerLibrary: storePhoton " << j << " npe " << npe << " " << pe[npe - 1];
588 #endif
589 }
Log< level::Info, true > LogVerbatim
HFShowerPhotonCollection pe
std::unique_ptr< HFShowerPhotonCollection > photo
HFShowerPhotonCollection photon

Member Data Documentation

◆ applyFidCut

bool HFShowerLibrary::applyFidCut
private

Definition at line 69 of file HFShowerLibrary.h.

Referenced by fillHits(), and HFShowerLibrary().

◆ backProb

double HFShowerLibrary::backProb
private

Definition at line 76 of file HFShowerLibrary.h.

Referenced by fillHits(), and HFShowerLibrary().

◆ dphi

double HFShowerLibrary::dphi
private

Definition at line 77 of file HFShowerLibrary.h.

Referenced by fillHits(), and HFShowerLibrary().

◆ emBranch

TBranch* HFShowerLibrary::emBranch
private

Definition at line 67 of file HFShowerLibrary.h.

Referenced by getRecord(), and HFShowerLibrary().

◆ equalizeTimeShift_

bool HFShowerLibrary::equalizeTimeShift_
private

Definition at line 75 of file HFShowerLibrary.h.

Referenced by fillHits(), and HFShowerLibrary().

◆ evtPerBin

int HFShowerLibrary::evtPerBin
private

Definition at line 70 of file HFShowerLibrary.h.

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

◆ fibre_

std::unique_ptr<HFFibre> HFShowerLibrary::fibre_
private

Definition at line 65 of file HFShowerLibrary.h.

Referenced by fillHits(), and HFShowerLibrary().

◆ fileVersion_

int HFShowerLibrary::fileVersion_
private

Definition at line 74 of file HFShowerLibrary.h.

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

◆ gpar

std::vector<double> HFShowerLibrary::gpar
private

Definition at line 78 of file HFShowerLibrary.h.

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

◆ hadBranch

TBranch * HFShowerLibrary::hadBranch
private

Definition at line 67 of file HFShowerLibrary.h.

Referenced by getRecord(), and HFShowerLibrary().

◆ hcalConstant_

const HcalDDDSimConstants* HFShowerLibrary::hcalConstant_
private

Definition at line 64 of file HFShowerLibrary.h.

Referenced by HFShowerLibrary().

◆ hf

TFile* HFShowerLibrary::hf
private

Definition at line 66 of file HFShowerLibrary.h.

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

◆ libVers

float HFShowerLibrary::libVers
private

Definition at line 71 of file HFShowerLibrary.h.

Referenced by HFShowerLibrary(), and loadEventInfo().

◆ listVersion

float HFShowerLibrary::listVersion
private

Definition at line 71 of file HFShowerLibrary.h.

Referenced by HFShowerLibrary(), and loadEventInfo().

◆ newForm

bool HFShowerLibrary::newForm
private

Definition at line 69 of file HFShowerLibrary.h.

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

◆ nMomBin

int HFShowerLibrary::nMomBin
private

Definition at line 70 of file HFShowerLibrary.h.

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

◆ npe

int HFShowerLibrary::npe
private

Definition at line 80 of file HFShowerLibrary.h.

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

◆ pe

HFShowerPhotonCollection HFShowerLibrary::pe
private

Definition at line 81 of file HFShowerLibrary.h.

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

◆ photo

std::unique_ptr<HFShowerPhotonCollection> HFShowerLibrary::photo
private

Definition at line 82 of file HFShowerLibrary.h.

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

◆ photon

HFShowerPhotonCollection HFShowerLibrary::photon
private

Definition at line 83 of file HFShowerLibrary.h.

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

◆ pmom

std::vector<double> HFShowerLibrary::pmom
private

Definition at line 72 of file HFShowerLibrary.h.

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

◆ probMax

double HFShowerLibrary::probMax
private

Definition at line 76 of file HFShowerLibrary.h.

Referenced by fillHits(), and HFShowerLibrary().

◆ rMax

double HFShowerLibrary::rMax
private

Definition at line 77 of file HFShowerLibrary.h.

Referenced by HFShowerLibrary(), and rInside().

◆ rMin

double HFShowerLibrary::rMin
private

Definition at line 77 of file HFShowerLibrary.h.

Referenced by HFShowerLibrary(), and rInside().

◆ totEvents

int HFShowerLibrary::totEvents
private

Definition at line 70 of file HFShowerLibrary.h.

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

◆ v3version

bool HFShowerLibrary::v3version
private

Definition at line 69 of file HFShowerLibrary.h.

Referenced by getRecord(), and HFShowerLibrary().

◆ verbose

bool HFShowerLibrary::verbose
private

Definition at line 69 of file HFShowerLibrary.h.

Referenced by HFShowerLibrary().