124 double tempZCalo = 26;
125 double hfcriticalEnergy = 0.021;
127 std::vector<HFGflash::Hit>
hit;
130 auto const preStepPoint = aStep->GetPreStepPoint();
131 auto const track = aStep->GetTrack();
136 const G4double energyCutoff = 1;
137 const G4int maxNumberOfSpots = 10000000;
139 G4ThreeVector showerStartingPosition =
track->GetPosition() / cm;
140 G4ThreeVector showerMomentum = preStepPoint->GetMomentum() / GeV;
143 const G4double
invgev = 1.0 / CLHEP::GeV;
144 G4double
energy = preStepPoint->GetTotalEnergy() *
invgev;
145 G4double logEinc =
std::log(energy);
147 G4double
y = energy / hfcriticalEnergy;
155 double charge =
track->GetStep()->GetPreStepPoint()->GetCharge();
156 theHelix->initializeTrajectory(showerMomentum, showerStartingPosition, charge,
theBField);
158 G4double pathLength0 =
theHelix->getPathLengthAtZ(showerStartingPosition.getZ());
159 G4double pathLength = pathLength0;
163 G4double fluctuatedTmax =
std::log(logY - 0.7157);
164 G4double fluctuatedAlpha =
std::log(0.7996 + (0.4581 + 1.8628 / tempZCalo) * logY);
166 G4double sigmaTmax = 1.0 / (-1.4 + 1.26 * logY);
167 G4double sigmaAlpha = 1.0 / (-0.58 + 0.86 * logY);
168 G4double
rho = 0.705 - 0.023 * logY;
169 G4double sqrtPL =
std::sqrt((1.0 + rho) / 2.0);
170 G4double sqrtLE =
std::sqrt((1.0 - rho) / 2.0);
172 G4double norm1 = G4RandGauss::shoot();
173 G4double norm2 = G4RandGauss::shoot();
174 G4double tempTmax = fluctuatedTmax + sigmaTmax * (sqrtPL * norm1 + sqrtLE * norm2);
175 G4double tempAlpha = fluctuatedAlpha + sigmaAlpha * (sqrtPL * norm1 - sqrtLE * norm2);
180 G4double
beta = (alpha - 1.0) /
tmax;
192 G4double averageTmax = logY - 0.858;
193 G4double averageAlpha = 0.21 + (0.492 + 2.38 / tempZCalo) * logY;
194 G4double spotTmax = averageTmax * (0.698 + .00212 * tempZCalo);
195 G4double spotAlpha = averageAlpha * (0.639 + .00334 * tempZCalo);
196 G4double spotBeta = (spotAlpha - 1.0) / spotTmax;
202 if (spotAlpha < 0.00001)
204 if (spotBeta < 0.00001)
208 edm::LogVerbatim(
"HFShower") <<
"HFGflash: Incoming energy = " << energy <<
" Position (rho,z) = ("
209 << showerStartingPosition.rho() <<
", " << showerStartingPosition.z() <<
")";
213 em_ssp_rho->Fill(showerStartingPosition.rho());
218 G4double z1 = 0.0251 + 0.00319 * logEinc;
219 G4double z2 = 0.1162 - 0.000381 * tempZCalo;
221 G4double k1 = 0.659 - 0.00309 * tempZCalo;
224 G4double k4 = 0.3585 + 0.0421 * logEinc;
226 G4double
p1 = 2.623 - 0.00094 * tempZCalo;
227 G4double
p2 = 0.401 + 0.00187 * tempZCalo;
228 G4double p3 = 1.313 - 0.0686 * logEinc;
232 const G4double e25Scale = 1.03551;
233 z1 *= 9.76972e-01 - 3.85026e-01 * std::tanh(1.82790
e+00 *
std::log(energy) - 3.66237
e+00);
236 G4double stepLengthLeft = 10000;
238 G4double zInX0 = 0.0;
239 G4double deltaZInX0 = 0.0;
240 G4double deltaZ = 0.0;
241 G4double stepLengthLeftInX0 = 0.0;
243 const G4double divisionStepInX0 = 0.1;
245 Genfun::IncompleteGamma gammaDist;
247 G4double energyInGamma = 0.0;
248 G4double preEnergyInGamma = 0.0;
249 G4double sigmaInGamma = 0.;
250 G4double preSigmaInGamma = 0.0;
253 G4double deltaEnergy = 0.0;
254 G4int spotCounter = 0;
257 G4double deltaStep = 0.0;
260 G4double timeGlobal = preStepPoint->GetGlobalTime();
263 G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume());
268 edm::LogVerbatim(
"HFShower") <<
"HFGflash: Energy = " << energy <<
" Step Length Left = " << stepLengthLeft;
270 while (energy > 0.0 && stepLengthLeft > 0.0) {
273 if (stepLengthLeftInX0 < divisionStepInX0) {
274 deltaZInX0 = stepLengthLeftInX0;
276 stepLengthLeft = 0.0;
278 deltaZInX0 = divisionStepInX0;
280 stepLengthLeft -= deltaZ;
286 edm::LogVerbatim(
"HFShower") <<
"HFGflash: zInX0 = " << zInX0 <<
" spotBeta*zInX0 = " << spotBeta * zInX0;
288 if ((!zInX0) || (!(spotBeta * zInX0 != 0)) || (zInX0 < 0.01) || (spotBeta * zInX0 < 0.00001) ||
289 (!(zInX0 * beta != 0)) || (zInX0 * beta < 0.00001))
292 G4int nSpotsInStep = 0;
295 edm::LogVerbatim(
"HFShower") <<
"HFGflash: Energy - Energy Cut off = " << energy - energyCutoff;
298 if (energy > energyCutoff) {
299 preEnergyInGamma = energyInGamma;
300 gammaDist.a().setValue(alpha);
302 energyInGamma = gammaDist(beta * zInX0);
303 G4double energyInDeltaZ = energyInGamma - preEnergyInGamma;
304 deltaEnergy =
std::min(energy, energy * energyInDeltaZ);
306 preSigmaInGamma = sigmaInGamma;
307 gammaDist.a().setValue(spotAlpha);
308 sigmaInGamma = gammaDist(spotBeta * zInX0);
309 nSpotsInStep =
std::max(1,
int(nSpots * (sigmaInGamma - preSigmaInGamma)));
312 preSigmaInGamma = sigmaInGamma;
313 nSpotsInStep =
std::max(1,
int(nSpots * (1.0 - preSigmaInGamma)));
316 if (deltaEnergy > energy || (energy - deltaEnergy) < energyCutoff)
319 energy -= deltaEnergy;
321 if (spotCounter + nSpotsInStep > maxNumberOfSpots) {
322 nSpotsInStep = maxNumberOfSpots - spotCounter;
323 if (nSpotsInStep < 1)
328 deltaStep += 0.5 * deltaZ;
329 pathLength += deltaStep;
330 deltaStep = 0.5 * deltaZ;
333 G4double tScale = tmax * alpha / (alpha - 1.0) * (1.0 -
std::exp(-fluctuatedAlpha));
334 G4double
tau =
std::min(10.0, (zInX0 - 0.5 * deltaZInX0) / tScale);
335 G4double rCore = z1 + z2 *
tau;
336 G4double rTail = k1 * (
std::exp(k3 * (tau - k2)) +
std::exp(k4 * (tau - k2)));
337 G4double p23 = (p2 -
tau) / p3;
343 G4double emSpotEnergy = deltaEnergy / (nSpotsInStep * e25Scale * GeV);
346 edm::LogVerbatim(
"HFShower") <<
"HFGflash: nSpotsInStep = " << nSpotsInStep;
349 for (G4int ispot = 0; ispot < nSpotsInStep; ispot++) {
351 G4double
u1 = G4UniformRand();
352 G4double
u2 = G4UniformRand();
353 G4double rInRM = 0.0;
355 if (u1 < probabilityWeight) {
356 rInRM = rCore *
std::sqrt(u2 / (1.0 - u2));
358 rInRM = rTail *
std::sqrt(u2 / (1.0 - u2));
364 G4double azimuthalAngle = twopi * G4UniformRand();
368 G4double incrementPath = (deltaZ / nSpotsInStep) * (ispot + 0.5 - 0.5 * nSpotsInStep);
372 theHelix->getGflashTrajectoryPoint(trajectoryPoint, pathLength + incrementPath);
374 G4ThreeVector SpotPosition0 = trajectoryPoint.
getPosition() +
387 timeGlobal += 0.0001 * nanosecond;
394 edm::LogVerbatim(
"HFShower") <<
"HFGflash: zInX0_spot,emSpotEnergy/GeV =" << zInX0_spot <<
" , "
395 << emSpotEnergy / GeV <<
"emSpotEnergy/GeV =" << emSpotEnergy / GeV;
400 if (emSpotEnergy <= 0)
409 em_long->Fill(SpotPosition0.z() / cm, emSpotEnergy *
invgev);
418 double energyratio = emSpotEnergy / (preStepPoint->GetTotalEnergy() / (nSpots * e25Scale));
420 if (emSpotEnergy / GeV < 0.0001)
422 if (energyratio > 80)
426 if (SpotPosition0.z() > 0)
428 if (SpotPosition0.z() < 0)
431 G4ThreeVector gfshift(0, 0, zshift * (
pow(100, 0.1) /
pow(energy, 0.1)));
433 G4ThreeVector SpotPosition = gfshift + SpotPosition0;
435 double LengthWeight = std::fabs(
std::pow(SpotPosition0.z() / 11370, 1));
436 if (G4UniformRand() > 0.0021 * energyratio * LengthWeight)
440 oneHit.
time = timeGlobal;
442 hit.push_back(oneHit);
Log< level::Info, true > LogVerbatim
static std::vector< std::string > checklist log
Gflash3Vector getCrossUnitVector()
Sin< T >::type sin(const T &t)
Exp< T >::type exp(const T &t)
Gflash::CalorimeterNumber jCalorimeter
const double rMoliere[kNumberCalorimeter]
static const float invgev
std::unique_ptr< GflashTrajectory > theHelix
Cos< T >::type cos(const T &t)
Abs< T >::type abs(const T &t)
Gflash3Vector getOrthogonalUnitVector()
static const double tmax[3]
const double radLength[kNumberCalorimeter]
Gflash3Vector & getPosition()
std::unique_ptr< G4Navigator > theGflashNavigator
Power< A, B >::type pow(const A &a, const B &b)