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;
147 G4double
y =
energy / hfcriticalEnergy;
155 double charge =
track->GetStep()->GetPreStepPoint()->GetCharge();
158 G4double pathLength0 =
theHelix->getPathLengthAtZ(showerStartingPosition.getZ());
159 G4double pathLength = pathLength0;
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;
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);
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)
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;
240 G4double stepLengthLeft = 10000;
241 G4double zInX0 = 0.0;
242 G4double deltaZInX0 = 0.0;
244 G4double stepLengthLeftInX0 = 0.0;
246 const G4double divisionStepInX0 = 0.1;
248 Genfun::IncompleteGamma gammaDist;
250 G4double energyInGamma = 0.0;
251 G4double preEnergyInGamma = 0.0;
252 G4double sigmaInGamma = 0.;
253 G4double preSigmaInGamma = 0.0;
256 G4double deltaEnergy = 0.0;
257 G4int spotCounter = 0;
260 G4double deltaStep = 0.0;
263 G4double timeGlobal = preStepPoint->GetGlobalTime();
266 G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume());
271 edm::LogVerbatim(
"HFShower") <<
"HFGflash: Energy = " <<
energy <<
" Step Length Left = " << stepLengthLeft;
273 while (
energy > 0.0 && stepLengthLeft > 0.0) {
276 if (stepLengthLeftInX0 < divisionStepInX0) {
277 deltaZInX0 = stepLengthLeftInX0;
279 stepLengthLeft = 0.0;
281 deltaZInX0 = divisionStepInX0;
289 edm::LogVerbatim(
"HFShower") <<
"HFGflash: zInX0 = " << zInX0 <<
" spotBeta*zInX0 = " << spotBeta * zInX0;
291 if ((!zInX0) || (!(spotBeta * zInX0 != 0)) || (zInX0 < 0.01) || (spotBeta * zInX0 < 0.00001) ||
292 (!(zInX0 *
beta != 0)) || (zInX0 *
beta < 0.00001))
295 G4int nSpotsInStep = 0;
301 if (
energy > energyCutoff) {
302 preEnergyInGamma = energyInGamma;
303 gammaDist.a().setValue(
alpha);
305 energyInGamma = gammaDist(
beta * zInX0);
306 G4double energyInDeltaZ = energyInGamma - preEnergyInGamma;
309 preSigmaInGamma = sigmaInGamma;
310 gammaDist.a().setValue(spotAlpha);
311 sigmaInGamma = gammaDist(spotBeta * zInX0);
312 nSpotsInStep =
std::max(1,
int(nSpots * (sigmaInGamma - preSigmaInGamma)));
315 preSigmaInGamma = sigmaInGamma;
316 nSpotsInStep =
std::max(1,
int(nSpots * (1.0 - preSigmaInGamma)));
319 if (deltaEnergy >
energy || (
energy - deltaEnergy) < energyCutoff)
324 if (spotCounter + nSpotsInStep > maxNumberOfSpots) {
325 nSpotsInStep = maxNumberOfSpots - spotCounter;
326 if (nSpotsInStep < 1)
331 deltaStep += 0.5 *
deltaZ;
332 pathLength += deltaStep;
337 G4double
tau =
std::min(10.0, (zInX0 - 0.5 * deltaZInX0) / tScale);
338 G4double rCore = z1 +
z2 *
tau;
340 G4double p23 = (
p2 -
tau) /
p3;
346 G4double emSpotEnergy = deltaEnergy / (nSpotsInStep * e25Scale * GeV);
349 edm::LogVerbatim(
"HFShower") <<
"HFGflash: nSpotsInStep = " << nSpotsInStep;
352 for (G4int ispot = 0; ispot < nSpotsInStep; ispot++) {
354 G4double
u1 = G4UniformRand();
355 G4double
u2 = G4UniformRand();
356 G4double rInRM = 0.0;
358 if (
u1 < probabilityWeight) {
367 G4double azimuthalAngle = twopi * G4UniformRand();
371 G4double incrementPath = (
deltaZ / nSpotsInStep) * (ispot + 0.5 - 0.5 * nSpotsInStep);
375 theHelix->getGflashTrajectoryPoint(trajectoryPoint, pathLength + incrementPath);
377 G4ThreeVector SpotPosition0 = trajectoryPoint.
getPosition() +
390 timeGlobal += 0.0001 * nanosecond;
397 edm::LogVerbatim(
"HFShower") <<
"HFGflash: zInX0_spot,emSpotEnergy/GeV =" << zInX0_spot <<
" , " 398 << emSpotEnergy / GeV <<
"emSpotEnergy/GeV =" << emSpotEnergy / GeV;
403 if (emSpotEnergy <= 0)
412 em_long->Fill(SpotPosition0.z() / cm, emSpotEnergy *
invgev);
421 double energyratio = emSpotEnergy / (preStepPoint->GetTotalEnergy() / (nSpots * e25Scale));
423 if (emSpotEnergy / GeV < 0.0001)
425 if (energyratio > 80)
429 if (SpotPosition0.z() > 0)
431 if (SpotPosition0.z() < 0)
434 G4ThreeVector gfshift(0, 0, zshift * (
pow(100, 0.1) /
pow(
energy, 0.1)));
436 G4ThreeVector SpotPosition = gfshift + SpotPosition0;
438 double LengthWeight = std::fabs(
std::pow(SpotPosition0.z() / 11370, 1));
439 if (G4UniformRand() > 0.0021 * energyratio * LengthWeight)
443 oneHit.
time = timeGlobal;
445 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)