129 double tempZCalo = 26;
130 double hfcriticalEnergy = 0.021;
132 std::vector<HFGflash::Hit>
hit;
135 auto const preStepPoint = aStep->GetPreStepPoint();
136 auto const track = aStep->GetTrack();
141 const G4double energyCutoff = 1;
142 const G4int maxNumberOfSpots = 10000000;
144 G4ThreeVector showerStartingPosition =
track->GetPosition() / cm;
145 G4ThreeVector showerMomentum = preStepPoint->GetMomentum() /
GeV;
149 G4double
energy = preStepPoint->GetTotalEnergy() *
invgev;
152 G4double
y =
energy / hfcriticalEnergy;
160 double charge =
track->GetStep()->GetPreStepPoint()->GetCharge();
164 G4double pathLength = pathLength0;
169 G4double fluctuatedAlpha =
std::log(0.7996 + (0.4581 + 1.8628 / tempZCalo) *
logY);
171 G4double sigmaTmax = 1.0 / (-1.4 + 1.26 *
logY);
172 G4double sigmaAlpha = 1.0 / (-0.58 + 0.86 *
logY);
173 G4double
rho = 0.705 - 0.023 *
logY;
177 G4double norm1 = G4RandGauss::shoot();
178 G4double norm2 = G4RandGauss::shoot();
179 G4double tempTmax = fluctuatedTmax + sigmaTmax * (sqrtPL * norm1 + sqrtLE * norm2);
180 G4double tempAlpha = fluctuatedAlpha + sigmaAlpha * (sqrtPL * norm1 - sqrtLE * norm2);
197 G4double averageTmax =
logY - 0.858;
198 G4double averageAlpha = 0.21 + (0.492 + 2.38 / tempZCalo) *
logY;
199 G4double spotTmax = averageTmax * (0.698 + .00212 * tempZCalo);
200 G4double spotAlpha = averageAlpha * (0.639 + .00334 * tempZCalo);
201 G4double spotBeta = (spotAlpha - 1.0) / spotTmax;
207 if (spotAlpha < 0.00001)
209 if (spotBeta < 0.00001)
214 << showerStartingPosition.rho() <<
", " << showerStartingPosition.z() <<
")";
218 em_ssp_rho->Fill(showerStartingPosition.rho());
223 G4double z1 = 0.0251 + 0.00319 * logEinc;
224 G4double
z2 = 0.1162 - 0.000381 * tempZCalo;
226 G4double k1 = 0.659 - 0.00309 * tempZCalo;
229 G4double k4 = 0.3585 + 0.0421 * logEinc;
231 G4double
p1 = 2.623 - 0.00094 * tempZCalo;
232 G4double
p2 = 0.401 + 0.00187 * tempZCalo;
233 G4double
p3 = 1.313 - 0.0686 * logEinc;
237 const G4double e25Scale = 1.03551;
238 z1 *= 9.76972e-01 - 3.85026e-01 * std::tanh(1.82790
e+00 *
std::log(
energy) - 3.66237
e+00);
241 G4double stepLengthLeft = 10000;
243 G4double zInX0 = 0.0;
244 G4double deltaZInX0 = 0.0;
245 G4double deltaZ = 0.0;
246 G4double stepLengthLeftInX0 = 0.0;
248 const G4double divisionStepInX0 = 0.1;
250 Genfun::IncompleteGamma gammaDist;
252 G4double energyInGamma = 0.0;
253 G4double preEnergyInGamma = 0.0;
254 G4double sigmaInGamma = 0.;
255 G4double preSigmaInGamma = 0.0;
258 G4double deltaEnergy = 0.0;
259 G4int spotCounter = 0;
262 G4double deltaStep = 0.0;
265 G4double timeGlobal = preStepPoint->GetGlobalTime();
270 G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume());
275 edm::LogVerbatim(
"HFShower") <<
"HFGflash: Energy = " <<
energy <<
" Step Length Left = " << stepLengthLeft;
277 while (
energy > 0.0 && stepLengthLeft > 0.0) {
280 if (stepLengthLeftInX0 < divisionStepInX0) {
281 deltaZInX0 = stepLengthLeftInX0;
283 stepLengthLeft = 0.0;
285 deltaZInX0 = divisionStepInX0;
287 stepLengthLeft -= deltaZ;
293 edm::LogVerbatim(
"HFShower") <<
"HFGflash: zInX0 = " << zInX0 <<
" spotBeta*zInX0 = " << spotBeta * zInX0;
295 if ((!zInX0) || (!(spotBeta * zInX0 != 0)) || (zInX0 < 0.01) || (spotBeta * zInX0 < 0.00001) ||
296 (!(zInX0 *
beta != 0)) || (zInX0 *
beta < 0.00001))
299 G4int nSpotsInStep = 0;
305 if (
energy > energyCutoff) {
306 preEnergyInGamma = energyInGamma;
307 gammaDist.a().setValue(
alpha);
309 energyInGamma = gammaDist(
beta * zInX0);
310 G4double energyInDeltaZ = energyInGamma - preEnergyInGamma;
313 preSigmaInGamma = sigmaInGamma;
314 gammaDist.a().setValue(spotAlpha);
315 sigmaInGamma = gammaDist(spotBeta * zInX0);
316 nSpotsInStep =
std::max(1,
int(nSpots * (sigmaInGamma - preSigmaInGamma)));
319 preSigmaInGamma = sigmaInGamma;
320 nSpotsInStep =
std::max(1,
int(nSpots * (1.0 - preSigmaInGamma)));
323 if (deltaEnergy >
energy || (
energy - deltaEnergy) < energyCutoff)
328 if (spotCounter + nSpotsInStep > maxNumberOfSpots) {
329 nSpotsInStep = maxNumberOfSpots - spotCounter;
330 if (nSpotsInStep < 1)
335 deltaStep += 0.5 * deltaZ;
336 pathLength += deltaStep;
337 deltaStep = 0.5 * deltaZ;
341 G4double
tau =
std::min(10.0, (zInX0 - 0.5 * deltaZInX0) / tScale);
342 G4double rCore = z1 +
z2 *
tau;
344 G4double p23 = (
p2 -
tau) /
p3;
350 G4double emSpotEnergy = deltaEnergy / (nSpotsInStep * e25Scale *
GeV);
353 edm::LogVerbatim(
"HFShower") <<
"HFGflash: nSpotsInStep = " << nSpotsInStep;
356 for (G4int ispot = 0; ispot < nSpotsInStep; ispot++) {
358 G4double
u1 = G4UniformRand();
359 G4double
u2 = G4UniformRand();
360 G4double rInRM = 0.0;
362 if (
u1 < probabilityWeight) {
371 G4double azimuthalAngle = twopi * G4UniformRand();
375 G4double incrementPath = (deltaZ / nSpotsInStep) * (ispot + 0.5 - 0.5 * nSpotsInStep);
381 G4ThreeVector SpotPosition0 = trajectoryPoint.
getPosition() +
394 timeGlobal += 0.0001 * nanosecond;
401 edm::LogVerbatim(
"HFShower") <<
"HFGflash: zInX0_spot,emSpotEnergy/GeV =" << zInX0_spot <<
" , "
402 << emSpotEnergy /
GeV <<
"emSpotEnergy/GeV =" << emSpotEnergy /
GeV;
407 if (emSpotEnergy <= 0)
416 em_long->Fill(SpotPosition0.z() / cm, emSpotEnergy *
invgev);
425 double energyratio = emSpotEnergy / (preStepPoint->GetTotalEnergy() / (nSpots * e25Scale));
427 if (emSpotEnergy /
GeV < 0.0001)
429 if (energyratio > 80)
433 if (SpotPosition0.z() > 0)
435 if (SpotPosition0.z() < 0)
438 G4ThreeVector gfshift(0, 0, zshift * (
pow(100, 0.1) /
pow(
energy, 0.1)));
440 G4ThreeVector SpotPosition = gfshift + SpotPosition0;
442 double LengthWeight = std::fabs(
std::pow(SpotPosition0.z() / 11370, 1));
443 if (G4UniformRand() > 0.0021 * energyratio * LengthWeight)
447 oneHit.
time = timeGlobal;
449 hit.push_back(oneHit);