12 #include "G4VPhysicalVolume.hh"
13 #include "G4NavigationHistory.hh"
16 #include "G4ParticleTable.hh"
17 #include "Randomize.hh"
18 #include "CLHEP/Units/SystemOfUnits.h"
19 #include "CLHEP/Units/PhysicalConstants.h"
27 : hcalConstant_(hcons),
hf(nullptr), emBranch(nullptr), hadBranch(nullptr), npe(0) {
36 backProb = m_HS.getParameter<
double>(
"BackProbability");
40 "BranchEvt",
"HFShowerLibraryEventInfos_hfshowerlib_HFShowerLibraryEventInfo");
41 std::string branchPre = m_HS.getUntrackedParameter<
std::string>(
"BranchPre",
"HFShowerPhotons_hfshowerlib_");
43 verbose = m_HS.getUntrackedParameter<
bool>(
"Verbosity",
false);
44 applyFidCut = m_HS.getParameter<
bool>(
"ApplyFiducialCut");
46 if (pTreeName.find(
'.') == 0)
47 pTreeName.erase(0, 2);
48 const char* nTree = pTreeName.c_str();
49 hf = TFile::Open(nTree);
52 edm::LogError(
"HFShower") <<
"HFShowerLibrary: opening " << nTree <<
" failed";
53 throw cms::Exception(
"Unknown",
"HFShowerLibrary") <<
"Opening of " << pTreeName <<
" fails\n";
55 edm::LogVerbatim(
"HFShower") <<
"HFShowerLibrary: opening " << nTree <<
" successfully";
58 newForm = (branchEvInfo.empty());
59 TTree*
event(
nullptr);
61 event = (TTree*)
hf->Get(
"HFSimHits");
63 event = (TTree*)
hf->Get(
"Events");
65 TBranch* evtInfo(
nullptr);
68 evtInfo =
event->GetBranch(
info.c_str());
73 edm::LogError(
"HFShower") <<
"HFShowerLibrary: HFShowerLibrayEventInfo"
74 <<
" Branch does not exist in Event";
75 throw cms::Exception(
"Unknown",
"HFShowerLibrary") <<
"Event information absent\n";
78 edm::LogError(
"HFShower") <<
"HFShowerLibrary: Events Tree does not "
80 throw cms::Exception(
"Unknown",
"HFShowerLibrary") <<
"Events tree absent\n";
87 ss <<
"HFShowerLibrary: Energies (GeV) with " <<
nMomBin <<
" bins\n";
89 if (
i / 10 * 10 ==
i &&
i > 0) {
96 std::string nameBr = branchPre + emName + branchPost;
97 emBranch =
event->GetBranch(nameBr.c_str());
100 nameBr = branchPre + hadName + branchPost;
101 hadBranch =
event->GetBranch(nameBr.c_str());
112 <<
" entries and Branch " << hadName <<
" has " <<
hadBranch->GetEntries() <<
" entries"
113 <<
"\n HFShowerLibrary::No packing information -"
114 <<
" Assume x, y, z are not in packed form"
115 <<
"\n Maximum probability cut off " <<
probMax <<
" Back propagation of light prob. "
124 rMax = rTable[rTable.size() - 1];
131 <<
" (Half) Phi Width of wedge " <<
dphi / CLHEP::deg;
147 auto const preStepPoint = aStep->GetPreStepPoint();
148 auto const postStepPoint = aStep->GetPostStepPoint();
149 auto const track = aStep->GetTrack();
151 auto const aParticle =
track->GetDynamicParticle();
152 const G4ThreeVector& momDir = aParticle->GetMomentumDirection();
153 const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
154 int parCode =
track->GetDefinition()->GetPDGEncoding();
158 if (
track->GetDefinition()->IsGeneralIon()) {
159 parCode = 1000020040;
163 G4String partType =
track->GetDefinition()->GetParticleName();
164 const G4ThreeVector localPos = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
165 double zoff = localPos.z() + 0.5 *
gpar[1];
167 edm::LogVerbatim(
"HFShower") <<
"HFShowerLibrary::getHits " << partType <<
" of energy "
169 <<
" onlyLong: " << onlyLong <<
" dir.orts " << momDir.x() <<
", " << momDir.y() <<
", "
170 << momDir.z() <<
" Pos x,y,z = " << hitPoint.x() <<
"," << hitPoint.y() <<
","
171 << hitPoint.z() <<
" (" << zoff <<
") sphi,cphi,stheta,ctheta = " <<
sin(momDir.phi())
172 <<
"," <<
cos(momDir.phi()) <<
", " <<
sin(momDir.theta()) <<
"," <<
cos(momDir.theta());
175 double tSlice = (postStepPoint->GetGlobalTime()) / CLHEP::nanosecond;
178 double pin = (
track->GetDefinition()->GetBaryonNumber() > 0) ? preStepPoint->GetKineticEnergy()
179 : preStepPoint->GetTotalEnergy();
181 return fillHits(hitPoint, momDir, parCode, pin, isKilled,
weight, tSlice, onlyLong);
185 const G4ThreeVector& momDir,
192 std::vector<HFShowerLibrary::Hit>
hit;
207 double pz = momDir.
z();
208 double zint = hitPoint.z();
211 bool backward = (pz * zint < 0.) ?
true :
false;
213 double sphi =
sin(momDir.phi());
214 double cphi =
cos(momDir.phi());
215 double ctheta =
cos(momDir.theta());
216 double stheta =
sin(momDir.theta());
234 for (
int i = 0;
i <
npe; ++
i) {
239 if (
zv <=
gpar[1] &&
pe[
i].lambda() > 0 && (
pe[
i].
z() >= 0 || (
zv >
gpar[0] && (!onlyLong)))) {
242 }
else if (!backward) {
246 double r = G4UniformRand();
253 double pex =
pe[
i].x();
254 double pey =
pe[
i].y();
256 double xx = pex * ctheta * cphi - pey * sphi +
zv * stheta * cphi;
257 double yy = pex * ctheta * sphi + pey * cphi +
zv * stheta * sphi;
258 double zz = -pex * stheta +
zv * ctheta;
260 G4ThreeVector
pos = hitPoint + G4ThreeVector(
xx,
yy,
zz);
262 G4ThreeVector lpos = G4ThreeVector(
pos.x(),
pos.y(),
zv);
266 double r =
pos.perp();
268 double fi =
pos.phi();
271 int isect =
int(fi /
dphi) + 1;
272 isect = (isect + 1) / 2;
273 double dfi = ((isect * 2 - 1) *
dphi - fi);
276 double dfir =
r *
sin(dfi);
278 edm::LogVerbatim(
"HFShower") <<
"HFShowerLibrary: Position shift " <<
xx <<
", " <<
yy <<
", " <<
zz <<
": "
279 <<
pos <<
" R " <<
r <<
" Phi " << fi <<
" Section " << isect <<
" R*Dfi " << dfir
283 double r1 = G4UniformRand();
284 double r2 = G4UniformRand();
285 double r3 = backward ? G4UniformRand() : -9999.;
291 <<
exp(-
p *
zv) <<
" r2 " <<
r2 <<
" r3 " << r3 <<
" rDfi " <<
gpar[5] <<
" zz "
293 <<
" rInside(r) :" <<
rInside(
r) <<
" r1 <= exp(-p*zv) :" << (
r1 <=
exp(-
p *
zv))
295 <<
" r3 <= backProb :" << (r3 <=
backProb)
296 <<
" dfir > gpar[5] :" << (dfir >
gpar[5]) <<
" zz >= gpar[4] :" << (
zz >=
gpar[4])
297 <<
" zz <= gpar[4]+gpar[1] :" << (
zz <=
gpar[4] +
gpar[1]);
304 hit.push_back(oneHit);
306 edm::LogVerbatim(
"HFShower") <<
"HFShowerLibrary: Final Hit " << nHit <<
" position " << (
hit[nHit].position)
307 <<
" Depth " << (
hit[nHit].
depth) <<
" Time " << tSlice <<
":" <<
pe[
i].t() <<
":"
317 r1 = G4UniformRand();
318 r2 = G4UniformRand();
322 oneHit.
time = (tSlice + (
pe[
i].t()) + (
fibre_->tShift(lpos, 2, 1)));
323 hit.push_back(oneHit);
325 edm::LogVerbatim(
"HFShower") <<
"HFShowerLibrary: Final Hit " << nHit <<
" position " << (
hit[nHit].position)
326 <<
" Depth " << (
hit[nHit].
depth) <<
" Time " << tSlice <<
":" <<
pe[
i].t()
327 <<
":" <<
fibre_->tShift(lpos, 2, 1) <<
":" << (
hit[nHit].time);
336 edm::LogVerbatim(
"HFShower") <<
"HFShowerLibrary: Total Hits " << nHit <<
" out of " <<
npe <<
" PE";
338 if (nHit >
npe && !onlyLong) {
339 edm::LogWarning(
"HFShower") <<
"HFShowerLibrary: Hit buffer " <<
npe <<
" smaller than " << nHit <<
" Hits";
356 std::vector<float>
t;
357 std::vector<float>*
tp = &
t;
360 unsigned int tSize =
t.size() / 5;
361 photo->reserve(tSize);
362 for (
unsigned int i = 0;
i < tSize;
i++) {
377 std::vector<float>
t;
378 std::vector<float>*
tp = &
t;
381 unsigned int tSize =
t.size() / 5;
382 photo->reserve(tSize);
383 for (
unsigned int i = 0;
i < tSize;
i++) {
396 << nPhoton <<
" photons";
397 for (
int j = 0;
j < nPhoton;
j++)
407 std::vector<HFShowerLibraryEventInfo> eventInfoCollection;
408 branch->SetAddress(&eventInfoCollection);
411 edm::LogVerbatim(
"HFShower") <<
"HFShowerLibrary::loadEventInfo loads EventInfo Collection of size "
412 << eventInfoCollection.size() <<
" records";
414 totEvents = eventInfoCollection[0].totalEvents();
415 nMomBin = eventInfoCollection[0].numberOfBins();
416 evtPerBin = eventInfoCollection[0].eventsPerBin();
417 libVers = eventInfoCollection[0].showerLibraryVersion();
418 listVersion = eventInfoCollection[0].physListVersion();
419 pmom = eventInfoCollection[0].energyBins();
422 edm::LogVerbatim(
"HFShower") <<
"HFShowerLibrary::loadEventInfo loads EventInfo from hardwired"
430 pmom = {2, 3, 5, 7, 10, 15, 20, 30, 50, 75, 100, 150, 250, 350, 500, 1000};
443 double r = G4UniformRand();
462 edm::LogWarning(
"HFShower") <<
"HFShowerLibrary:: Illegal irc[0] = " << irc[0] <<
" now set to 0";
472 edm::LogWarning(
"HFShower") <<
"HFShowerLibrary:: Illegal irc[1] = " << irc[1] <<
" now set to 1";
480 edm::LogVerbatim(
"HFShower") <<
"HFShowerLibrary:: Select records " << irc[0] <<
" and " << irc[1] <<
" with weights "
481 << 1 -
w <<
" and " <<
w;
486 for (
int ir = 0; ir < 2; ir++) {
491 for (
int j = 0;
j < nPhoton;
j++) {
493 if ((ir == 0 &&
r >
w) || (ir > 0 &&
r <
w)) {
500 if ((
npe > npold || (npold == 0 && irc[0] > 0)) && !(
npe == 0 && npold == 0))
501 edm::LogWarning(
"HFShower") <<
"HFShowerLibrary: Interpolation Warning =="
502 <<
" records " << irc[0] <<
" and " << irc[1] <<
" gives a buffer of " << npold
503 <<
" photons and fills " <<
npe <<
" *****";
506 edm::LogVerbatim(
"HFShower") <<
"HFShowerLibrary: Interpolation == records " << irc[0] <<
" and " << irc[1]
507 <<
" gives a buffer of " << npold <<
" photons and fills " <<
npe <<
" PE";
508 for (
int j = 0;
j <
npe;
j++)
519 <<
" momentum bins and " <<
evtPerBin <<
" entries/bin -- "
520 <<
"total " <<
totEvents <<
" using " << nrec <<
" records";
522 std::vector<int> irc(nrec);
524 for (
int ir = 0; ir < nrec; ir++) {
525 double r = G4UniformRand();
528 edm::LogWarning(
"HFShower") <<
"HFShowerLibrary:: Illegal irc[" << ir <<
"] = " << irc[ir] <<
" now set to 1";
531 edm::LogWarning(
"HFShower") <<
"HFShowerLibrary:: Illegal irc[" << ir <<
"] = " << irc[ir] <<
" now set to "
536 edm::LogVerbatim(
"HFShower") <<
"HFShowerLibrary::Extrapolation use irc[" << ir <<
"] = " << irc[ir];
544 for (
int ir = 0; ir < nrec; ir++) {
549 for (
int j = 0;
j < nPhoton;
j++) {
550 double r = G4UniformRand();
551 if (ir != nrec - 1 ||
r <
w) {
556 edm::LogVerbatim(
"HFShower") <<
"HFShowerLibrary: Record [" << ir <<
"] = " << irc[ir] <<
" npold = " << npold;
561 edm::LogVerbatim(
"HFShower") <<
"HFShowerLibrary:: uses " << npold <<
" photons";
564 if (
npe > npold || npold == 0)
565 edm::LogWarning(
"HFShower") <<
"HFShowerLibrary: Extrapolation Warning == " << nrec <<
" records " << irc[0] <<
", "
566 << irc[1] <<
", ... gives a buffer of " << npold <<
" photons and fills " <<
npe
570 edm::LogVerbatim(
"HFShower") <<
"HFShowerLibrary: Extrapolation == " << nrec <<
" records " << irc[0] <<
", "
571 << irc[1] <<
", ... gives a buffer of " << npold <<
" photons and fills " <<
npe
573 for (
int j = 0;
j <
npe;
j++)