13 #include "G4VPhysicalVolume.hh" 16 #include "Randomize.hh" 17 #include <CLHEP/Units/SystemOfUnits.h> 29 emPDG = theParticleTable->FindParticle(parName =
"e-")->GetPDGEncoding();
30 epPDG = theParticleTable->FindParticle(parName =
"e+")->GetPDGEncoding();
31 gammaPDG = theParticleTable->FindParticle(parName =
"gamma")->GetPDGEncoding();
32 pi0PDG = theParticleTable->FindParticle(parName =
"pi0")->GetPDGEncoding();
33 etaPDG = theParticleTable->FindParticle(parName =
"eta")->GetPDGEncoding();
34 nuePDG = theParticleTable->FindParticle(parName =
"nu_e")->GetPDGEncoding();
35 numuPDG = theParticleTable->FindParticle(parName =
"nu_mu")->GetPDGEncoding();
36 nutauPDG = theParticleTable->FindParticle(parName =
"nu_tau")->GetPDGEncoding();
37 anuePDG = theParticleTable->FindParticle(parName =
"anti_nu_e")->GetPDGEncoding();
38 anumuPDG = theParticleTable->FindParticle(parName =
"anti_nu_mu")->GetPDGEncoding();
39 anutauPDG = theParticleTable->FindParticle(parName =
"anti_nu_tau")->GetPDGEncoding();
40 geantinoPDG = theParticleTable->FindParticle(parName =
"geantino")->GetPDGEncoding();
49 G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
50 G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
51 G4Track*
track = aStep->GetTrack();
53 const G4DynamicParticle* aParticle =
track->GetDynamicParticle();
54 const G4ThreeVector& momDir = aParticle->GetMomentumDirection();
55 double energy = preStepPoint->GetKineticEnergy();
56 G4ThreeVector hitPoint = preStepPoint->GetPosition();
57 const G4ThreeVector& hitPointOrig = preStepPoint->GetPosition();
58 G4int parCode =
track->GetDefinition()->GetPDGEncoding();
68 G4ThreeVector posLocal;
69 double tSlice = (postStepPoint->GetGlobalTime()) / CLHEP::nanosecond;
76 double xxlocal, yylocal, zzlocal;
79 side = (hitPointOrig.z() > 0.) ?
true :
false;
85 for (
int i = 0;
i <
npe;
i++) {
94 zz = (hitPointOrig.z() > 0.) ? zzlocal +
Z0 : zzlocal -
Z0;
96 posLocal = G4ThreeVector(xxlocal, yylocal, zzlocal);
107 zz = (hitPointOrig.z() > 0.) ? zzlocal +
Z0 : zzlocal -
Z0;
109 posLocal = G4ThreeVector(xxlocal, yylocal, zzlocal);
114 oneHit.
depth = channel;
115 oneHit.
time = tSlice;
119 hitPoint.setX(hitPointOrig.x() -
X0);
120 hitPoint.setY(hitPointOrig.y() -
Y0);
121 double setZ = (hitPointOrig.z() > 0.) ? hitPointOrig.z() -
Z0 : fabs(hitPointOrig.z()) -
Z0;
135 hits.push_back(oneHit);
137 edm::LogVerbatim(
"ZdcShower") <<
"\nZdcShowerLibrary:Generated Hit " << nHit <<
" orig hit pos " << hitPointOrig
138 <<
" orig hit pos local coord" << hitPoint <<
" new position " 139 << (
hits[nHit].position) <<
" Channel " << (
hits[nHit].
depth) <<
" side " <<
side 140 <<
" Time " << (
hits[nHit].time) <<
" DetectorID " << (
hits[nHit].detID)
141 <<
" Had Energy " << (
hits[nHit].DeHad) <<
" EM Energy " << (
hits[nHit].DeEM)
149 const G4ThreeVector& momDir,
159 edm::LogVerbatim(
"ZdcShower") <<
"\n ZdcShowerLibrary::getEnergyFromLibrary input/output variables:" 160 <<
" phi: " << momDir.phi() / CLHEP::deg <<
" theta: " << momDir.theta() / CLHEP::deg
161 <<
" xin : " << hitPoint.x() <<
" yin : " << hitPoint.y() <<
" zin : " << hitPoint.z()
162 <<
" track en: " <<
energy <<
"(GeV)" 163 <<
" section: " <<
section <<
" side: " <<
side <<
" channel: " << channel
164 <<
" partID: " << parCode;
179 float xin = hitPoint.x();
180 float yin = hitPoint.y();
183 if (
section == 1 && iparCode != 0) {
192 if (
section == 2 && iparCode != 0) {
202 if (
section == 1 && iparCode == 0) {
212 yin = yin / CLHEP::cm;
213 xin = xin / CLHEP::cm;
216 eav = ((((((-0.0002 * xin - 2.0e-13) * xin + 0.0022) * xin + 1.0e-11) * xin - 0.0217) * xin - 3.0e-10) * xin +
218 (((0.0001 * yin + 0.0056) * yin + 0.0508) * yin + 1.0) * 300.0 *
pow((
energy / 300.0), 0.99);
219 esig = ((((((0.0005 * xin - 1.0e-12) * xin - 0.0052) * xin + 5.0e-11) * xin + 0.032) * xin - 2.0e-10) * xin + 1.0) *
220 (((0.0006 * yin + 0.0071) * yin - 0.031) * yin + 1.0) * 30.0 *
pow((
energy / 300.0), 0.54);
223 eav = ((((((-0.0002 * xin - 2.0e-13) * xin + 0.0022) * xin + 1.0e-11) * xin - 0.0217) * xin - 3.0e-10) * xin +
225 (((0.0001 * yin + 0.0056) * yin + 0.0508) * yin + 1.0) * 300.0 *
pow((
energy / 300.0), 1.12);
226 esig = ((((((0.0005 * xin - 1.0e-12) * xin - 0.0052) * xin + 5.0e-11) * xin + 0.032) * xin - 2.0e-10) * xin + 1.0) *
227 (((0.0006 * yin + 0.0071) * yin - 0.031) * yin + 1.0) * 54.0 *
pow((
energy / 300.0), 0.93);
231 if (eav < 0. || edis < 0.) {
232 edm::LogVerbatim(
"ZdcShower") <<
" Negative everage energy from parametrization \n" 233 <<
" xin: " << xin <<
"(cm)" 234 <<
" yin: " << yin <<
"(cm)" 235 <<
" track en: " <<
energy <<
"(GeV)" 236 <<
" eaverage: " << eav / CLHEP::GeV <<
" (GeV)" 237 <<
" esigma: " << esig / CLHEP::GeV <<
" (GeV)" 238 <<
" edist: " << edis <<
" (GeV)" 239 <<
" dE hit: " <<
nphotons / CLHEP::GeV <<
" (GeV)";
244 eav = eav * CLHEP::GeV;
245 esig = esig * CLHEP::GeV;
251 <<
" eaverage: " << eav / CLHEP::GeV <<
" (GeV)" 252 <<
" esigma: " << esig / CLHEP::GeV <<
" (GeV)" 253 <<
" edist: " << edis <<
" (GeV)" 254 <<
" dE hit: " <<
nphotons / CLHEP::GeV <<
" (GeV)";
263 efluct = eav + esig * CLHEP::RandGaussQ::shoot();
265 efluct = eav + esig * CLHEP::RandLandau::shoot();
Log< level::Info, true > LogVerbatim
ALPAKA_FN_ACC int side(int ieta, int iphi)
ZdcShowerLibrary(const std::string &name, edm::ParameterSet const &p)
T getUntrackedParameter(std::string const &, T const &) const
static const double theXChannelBoundaries[]
int getEnergyFromLibrary(const G4ThreeVector &posHit, const G4ThreeVector &momDir, double energy, G4int parCode, HcalZDCDetId::Section section, bool side, int channel)
static const double theZHadChannelBoundaries[]
std::vector< ZdcShowerLibrary::Hit > hits
int photonFluctuation(double eav, double esig, double edis)
std::vector< Hit > & getHits(const G4Step *aStep, bool &ok)
int encodePartID(G4int parCode)
static const double theZSectionBoundaries[]
void initRun(G4ParticleTable *theParticleTable)
Power< A, B >::type pow(const A &a, const B &b)