127 double tempZCalo = 26;
128 double hfcriticalEnergy = 0.021;
130 std::vector<HFGflash::Hit>
hit;
133 auto const preStepPoint = aStep->GetPreStepPoint();
134 auto const track = aStep->GetTrack();
139 const G4double energyCutoff = 1;
140 const G4int maxNumberOfSpots = 10000000;
142 G4ThreeVector showerStartingPosition =
track->GetPosition() / cm;
143 G4ThreeVector showerMomentum = preStepPoint->GetMomentum() / GeV;
146 const G4double
invgev = 1.0 / CLHEP::GeV;
147 G4double
energy = preStepPoint->GetTotalEnergy() *
invgev;
150 G4double
y =
energy / hfcriticalEnergy;
158 double charge =
track->GetStep()->GetPreStepPoint()->GetCharge();
161 G4double pathLength0 =
theHelix->getPathLengthAtZ(showerStartingPosition.getZ());
162 G4double pathLength = pathLength0;
167 G4double fluctuatedAlpha =
std::log(0.7996 + (0.4581 + 1.8628 / tempZCalo) *
logY);
169 G4double sigmaTmax = 1.0 / (-1.4 + 1.26 *
logY);
170 G4double sigmaAlpha = 1.0 / (-0.58 + 0.86 *
logY);
171 G4double
rho = 0.705 - 0.023 *
logY;
175 G4double norm1 = G4RandGauss::shoot();
176 G4double norm2 = G4RandGauss::shoot();
177 G4double tempTmax = fluctuatedTmax + sigmaTmax * (sqrtPL * norm1 + sqrtLE * norm2);
178 G4double tempAlpha = fluctuatedAlpha + sigmaAlpha * (sqrtPL * norm1 - sqrtLE * norm2);
195 G4double averageTmax =
logY - 0.858;
196 G4double averageAlpha = 0.21 + (0.492 + 2.38 / tempZCalo) *
logY;
197 G4double spotTmax = averageTmax * (0.698 + .00212 * tempZCalo);
198 G4double spotAlpha = averageAlpha * (0.639 + .00334 * tempZCalo);
199 G4double spotBeta = (spotAlpha - 1.0) / spotTmax;
205 if (spotAlpha < 0.00001)
207 if (spotBeta < 0.00001)
212 << showerStartingPosition.rho() <<
", " << showerStartingPosition.z() <<
")";
216 em_ssp_rho->Fill(showerStartingPosition.rho());
221 G4double z1 = 0.0251 + 0.00319 * logEinc;
222 G4double
z2 = 0.1162 - 0.000381 * tempZCalo;
224 G4double k1 = 0.659 - 0.00309 * tempZCalo;
227 G4double k4 = 0.3585 + 0.0421 * logEinc;
229 G4double
p1 = 2.623 - 0.00094 * tempZCalo;
230 G4double
p2 = 0.401 + 0.00187 * tempZCalo;
231 G4double
p3 = 1.313 - 0.0686 * logEinc;
235 const G4double e25Scale = 1.03551;
243 G4double stepLengthLeft = 10000;
244 G4double zInX0 = 0.0;
245 G4double deltaZInX0 = 0.0;
247 G4double stepLengthLeftInX0 = 0.0;
249 const G4double divisionStepInX0 = 0.1;
251 Genfun::IncompleteGamma gammaDist;
253 G4double energyInGamma = 0.0;
254 G4double preEnergyInGamma = 0.0;
255 G4double sigmaInGamma = 0.;
256 G4double preSigmaInGamma = 0.0;
259 G4double deltaEnergy = 0.0;
260 G4int spotCounter = 0;
263 G4double deltaStep = 0.0;
266 G4double timeGlobal = preStepPoint->GetGlobalTime();
269 G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume());
274 edm::LogVerbatim(
"HFShower") <<
"HFGflash: Energy = " <<
energy <<
" Step Length Left = " << stepLengthLeft;
276 while (
energy > 0.0 && stepLengthLeft > 0.0) {
279 if (stepLengthLeftInX0 < divisionStepInX0) {
280 deltaZInX0 = stepLengthLeftInX0;
282 stepLengthLeft = 0.0;
284 deltaZInX0 = divisionStepInX0;
292 edm::LogVerbatim(
"HFShower") <<
"HFGflash: zInX0 = " << zInX0 <<
" spotBeta*zInX0 = " << spotBeta * zInX0;
294 if ((!zInX0) || (!(spotBeta * zInX0 != 0)) || (zInX0 < 0.01) || (spotBeta * zInX0 < 0.00001) ||
295 (!(zInX0 *
beta != 0)) || (zInX0 *
beta < 0.00001))
298 G4int nSpotsInStep = 0;
304 if (
energy > energyCutoff) {
305 preEnergyInGamma = energyInGamma;
306 gammaDist.a().setValue(
alpha);
308 energyInGamma = gammaDist(
beta * zInX0);
309 G4double energyInDeltaZ = energyInGamma - preEnergyInGamma;
312 preSigmaInGamma = sigmaInGamma;
313 gammaDist.a().setValue(spotAlpha);
314 sigmaInGamma = gammaDist(spotBeta * zInX0);
315 nSpotsInStep =
std::max(1,
int(nSpots * (sigmaInGamma - preSigmaInGamma)));
318 preSigmaInGamma = sigmaInGamma;
319 nSpotsInStep =
std::max(1,
int(nSpots * (1.0 - preSigmaInGamma)));
322 if (deltaEnergy >
energy || (
energy - deltaEnergy) < energyCutoff)
327 if (spotCounter + nSpotsInStep > maxNumberOfSpots) {
328 nSpotsInStep = maxNumberOfSpots - spotCounter;
329 if (nSpotsInStep < 1)
334 deltaStep += 0.5 *
deltaZ;
335 pathLength += deltaStep;
340 G4double
tau =
std::min(10.0, (zInX0 - 0.5 * deltaZInX0) / tScale);
341 G4double rCore = z1 +
z2 *
tau;
343 G4double p23 = (
p2 -
tau) /
p3;
349 G4double emSpotEnergy = deltaEnergy / (nSpotsInStep * e25Scale * GeV);
352 edm::LogVerbatim(
"HFShower") <<
"HFGflash: nSpotsInStep = " << nSpotsInStep;
355 for (G4int ispot = 0; ispot < nSpotsInStep; ispot++) {
357 G4double
u1 = G4UniformRand();
358 G4double
u2 = G4UniformRand();
359 G4double rInRM = 0.0;
361 if (
u1 < probabilityWeight) {
370 G4double azimuthalAngle = twopi * G4UniformRand();
374 G4double incrementPath = (
deltaZ / nSpotsInStep) * (ispot + 0.5 - 0.5 * nSpotsInStep);
378 theHelix->getGflashTrajectoryPoint(trajectoryPoint, pathLength + incrementPath);
380 G4ThreeVector SpotPosition0 = trajectoryPoint.
getPosition() +
393 timeGlobal += 0.0001 * nanosecond;
400 edm::LogVerbatim(
"HFShower") <<
"HFGflash: zInX0_spot,emSpotEnergy/GeV =" << zInX0_spot <<
" , " 401 << emSpotEnergy / GeV <<
"emSpotEnergy/GeV =" << emSpotEnergy / GeV;
406 if (emSpotEnergy <= 0)
415 em_long->Fill(SpotPosition0.z() / cm, emSpotEnergy *
invgev);
424 double energyratio = emSpotEnergy / (preStepPoint->GetTotalEnergy() / (nSpots * e25Scale));
426 if (emSpotEnergy / GeV < 0.0001)
428 if (energyratio > 80)
432 if (SpotPosition0.z() > 0)
434 if (SpotPosition0.z() < 0)
437 G4ThreeVector gfshift(0, 0, zshift * (
pow(100, 0.1) /
pow(
energy, 0.1)));
439 G4ThreeVector SpotPosition = gfshift + SpotPosition0;
441 double LengthWeight = std::fabs(
std::pow(SpotPosition0.z() / 11370, 1));
442 if (G4UniformRand() > 0.0021 * energyratio * LengthWeight)
446 oneHit.
time = timeGlobal;
448 hit.push_back(oneHit);
Log< level::Info, true > LogVerbatim
void tanh(data_T data[CONFIG_T::n_in], res_T res[CONFIG_T::n_in])
Gflash3Vector getCrossUnitVector()
Sin< T >::type sin(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)