92 double tempZCalo = 26;
93 double hfcriticalEnergy = 0.021;
95 std::vector<HFGflash::Hit>
hit;
98 G4StepPoint * preStepPoint = aStep->GetPreStepPoint();
99 G4Track *
track = aStep->GetTrack();
100 G4String partType = track->GetDefinition()->GetParticleName();
105 const G4double energyCutoff = 1;
106 const G4int maxNumberOfSpots = 10000000;
108 G4ThreeVector showerStartingPosition = track->GetPosition()/cm;
109 G4ThreeVector showerMomentum = preStepPoint->GetMomentum()/
GeV;
112 G4double logEinc =
std::log((preStepPoint->GetTotalEnergy())/
GeV);
114 G4double
y = ((preStepPoint->GetTotalEnergy())/
GeV) / hfcriticalEnergy;
118 G4double nSpots = 93.0 *
std::log(tempZCalo) * ((preStepPoint->GetTotalEnergy())/
GeV);
119 if(preStepPoint->GetTotalEnergy()/
GeV < 1.6) nSpots = 140.4 *
std::log(tempZCalo) *
std::pow(((preStepPoint->GetTotalEnergy())/
GeV),0.876);
123 double charge = track->GetStep()->GetPreStepPoint()->GetCharge();
127 G4double pathLength = pathLength0;
131 G4double fluctuatedTmax =
std::log(logY - 0.7157);
132 G4double fluctuatedAlpha=
std::log(0.7996 +(0.4581 + 1.8628/tempZCalo)*logY);
134 G4double sigmaTmax = 1.0/( -1.4 + 1.26 *
logY);
135 G4double sigmaAlpha = 1.0/( -0.58 + 0.86 *
logY);
136 G4double
rho = 0.705 - 0.023 *
logY;
137 G4double sqrtPL =
std::sqrt((1.0+rho)/2.0);
138 G4double sqrtLE =
std::sqrt((1.0-rho)/2.0);
140 G4double norm1 = G4RandGauss::shoot();
141 G4double norm2 = G4RandGauss::shoot();
142 G4double tempTmax = fluctuatedTmax + sigmaTmax*(sqrtPL*norm1 + sqrtLE*norm2);
143 G4double tempAlpha = fluctuatedAlpha + sigmaAlpha*(sqrtPL*norm1 - sqrtLE*norm2);
150 if (!alpha)
return hit;
151 if (!beta)
return hit;
152 if (alpha < 0.00001)
return hit;
153 if (beta < 0.00001)
return hit;
156 G4double averageTmax = logY-0.858;
157 G4double averageAlpha = 0.21+(0.492+2.38/tempZCalo)*logY;
158 G4double spotTmax = averageTmax * (0.698 + .00212*tempZCalo);
159 G4double spotAlpha= averageAlpha * (0.639 + .00334*tempZCalo);
160 G4double spotBeta = (spotAlpha-1.0)/spotTmax;
162 if (!spotAlpha)
return hit;
163 if (!spotBeta)
return hit;
164 if (spotAlpha < 0.00001)
return hit;
165 if (spotBeta < 0.00001)
return hit;
168 LogDebug(
"HFShower") <<
"Incoming energy = " << ((preStepPoint->GetTotalEnergy())/
GeV) <<
" Position (rho,z) = (" << showerStartingPosition.rho() <<
", " << showerStartingPosition.z() <<
")";
171 em_incE->Fill(((preStepPoint->GetTotalEnergy())/
GeV));
172 em_ssp_rho->Fill(showerStartingPosition.rho());
177 G4double z1=0.0251+0.00319*logEinc;
178 G4double z2=0.1162-0.000381*tempZCalo;
180 G4double k1=0.659 - 0.00309 * tempZCalo;
183 G4double k4=0.3585+ 0.0421*logEinc;
185 G4double
p1=2.623 -0.00094*tempZCalo;
186 G4double
p2=0.401 +0.00187*tempZCalo;
187 G4double
p3=1.313 -0.0686*logEinc;
191 G4double e25Scale = 1.03551;
192 z1 *= 9.76972e-01 - 3.85026e-01 * std::tanh(1.82790
e+00*
std::log(((preStepPoint->GetTotalEnergy())/
GeV)) - 3.66237
e+00);
195 G4double stepLengthLeft = 10000;
197 G4double zInX0 = 0.0;
198 G4double deltaZInX0 = 0.0;
199 G4double deltaZ = 0.0;
200 G4double stepLengthLeftInX0 = 0.0;
202 const G4double divisionStepInX0 = 0.1;
203 G4double energy = ((preStepPoint->GetTotalEnergy())/
GeV);
205 Genfun::IncompleteGamma gammaDist;
207 G4double energyInGamma = 0.0;
208 G4double preEnergyInGamma = 0.0;
209 G4double sigmaInGamma = 0.;
210 G4double preSigmaInGamma = 0.0;
213 G4double deltaEnergy =0.0 ;
214 G4int spotCounter = 0;
217 G4double deltaStep = 0.0;
220 G4double timeGlobal = track->GetStep()->GetPreStepPoint()->GetGlobalTime();
224 theGflashNavigator->SetWorldVolume(G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume());
229 LogDebug(
"HFShower") <<
" Energy = " << energy <<
" Step Length Left = " << stepLengthLeft;
231 while(energy > 0.0 && stepLengthLeft > 0.0) {
234 if ( stepLengthLeftInX0 < divisionStepInX0 ) {
235 deltaZInX0 = stepLengthLeftInX0;
237 stepLengthLeft = 0.0;
239 deltaZInX0 = divisionStepInX0;
241 stepLengthLeft -= deltaZ;
247 LogDebug(
"HFShower") <<
" zInX0 = " << zInX0 <<
" spotBeta*zInX0 = " << spotBeta*zInX0;
249 if ((!zInX0) || (! (spotBeta*zInX0 != 0) ) || (zInX0 < 0.01) ||
250 (spotBeta*zInX0 < 0.00001) || (! (zInX0*beta != 0) ) || (zInX0*beta < 0.00001))
253 G4int nSpotsInStep = 0;
256 LogDebug(
"HFShower") <<
" Energy - Energy Cut off = " << energy - energyCutoff;
259 if ( energy > energyCutoff ) {
260 preEnergyInGamma = energyInGamma;
261 gammaDist.a().setValue(alpha);
263 energyInGamma = gammaDist(beta*zInX0);
264 G4double energyInDeltaZ = energyInGamma - preEnergyInGamma;
265 deltaEnergy =
std::min(energy,((preStepPoint->GetTotalEnergy())/
GeV)*energyInDeltaZ);
267 preSigmaInGamma = sigmaInGamma;
268 gammaDist.a().setValue(spotAlpha);
269 sigmaInGamma = gammaDist(spotBeta*zInX0);
270 nSpotsInStep =
std::max(1,
int(nSpots * (sigmaInGamma - preSigmaInGamma)));
272 deltaEnergy = energy;
273 preSigmaInGamma = sigmaInGamma;
274 nSpotsInStep =
std::max(1,
int(nSpots * (1.0 - preSigmaInGamma)));
277 if ( deltaEnergy > energy || (energy-deltaEnergy) < energyCutoff ) deltaEnergy = energy;
279 energy -= deltaEnergy;
281 if ( spotCounter+nSpotsInStep > maxNumberOfSpots ) {
282 nSpotsInStep = maxNumberOfSpots - spotCounter;
283 if (nSpotsInStep < 1) nSpotsInStep = 1;
288 deltaStep += 0.5*deltaZ;
289 pathLength += deltaStep;
290 deltaStep = 0.5*deltaZ;
294 G4double tScale = tmax *alpha/(alpha-1.0) * (
std::exp(fluctuatedAlpha)-1.0)/
std::exp(fluctuatedAlpha);
295 G4double
tau =
std::min(10.0,(zInX0 - 0.5*deltaZInX0)/tScale);
296 G4double rCore = z1 + z2 *
tau;
298 G4double p23 = (p2 -
tau)/p3;
305 G4double emSpotEnergy = deltaEnergy / nSpotsInStep * e25Scale *
GeV;
309 LogDebug(
"HFShower") <<
" nSpotsInStep = " << nSpotsInStep;
314 for (G4int ispot = 0 ; ispot < nSpotsInStep ; ispot++) {
316 G4double
u1 = G4UniformRand();
317 G4double
u2 = G4UniformRand();
318 G4double rInRM = 0.0;
320 if (u1 < probabilityWeight) {
321 rInRM = rCore *
std::sqrt( u2/(1.0-u2) );
323 rInRM = rTail *
std::sqrt( u2/(1.0-u2) );
329 G4double azimuthalAngle = twopi*G4UniformRand();
333 G4double incrementPath = (deltaZ/nSpotsInStep)*(ispot+0.5 - 0.5*nSpotsInStep);
352 timeGlobal += 0.0001*nanosecond;
359 LogDebug(
"HFShower") <<
"zInX0_spot,emSpotEnergy/GeV =" << zInX0_spot <<
" , " << emSpotEnergy/GeV <<
"emSpotEnergy/GeV =" << emSpotEnergy/
GeV;
362 if ((!zInX0_spot) || (zInX0_spot < 0))
continue;
363 if ((!emSpotEnergy/GeV) || (emSpotEnergy < 0))
continue;
364 if ((!rShower/Gflash::rMoliere[
jCalorimeter]) || (rShower/Gflash::rMoliere[jCalorimeter] < 0))
continue;
370 em_long->Fill(SpotPosition0.z()/cm,emSpotEnergy/
GeV);
371 em_lateral->Fill(rShower/Gflash::rMoliere[jCalorimeter],emSpotEnergy/GeV);
379 double energyratio = emSpotEnergy/(preStepPoint->GetTotalEnergy()/(nSpots*e25Scale));
381 if (emSpotEnergy/GeV < 0.0001)
continue;
382 if (energyratio > 80)
continue;
385 if(SpotPosition0.z() > 0) zshift = 18;
386 if(SpotPosition0.z() < 0) zshift = -18;
389 G4ThreeVector gfshift(0,0,zshift*(
pow(100,0.1)/
pow(preStepPoint->GetTotalEnergy()/
GeV,0.1)));
391 G4ThreeVector SpotPosition = gfshift + SpotPosition0;
393 double LengthWeight = std::fabs(
std::pow(SpotPosition0.z()/11370,1));
394 if (G4UniformRand()> 0.0021 * energyratio * LengthWeight)
continue;
397 oneHit.
time = timeGlobal;
398 oneHit.
edep = emSpotEnergy/
GeV;
399 hit.push_back(oneHit);
G4Navigator * theGflashNavigator
double getPathLengthAtZ(double z) const
Gflash3Vector getCrossUnitVector()
Sin< T >::type sin(const T &t)
Gflash::CalorimeterNumber jCalorimeter
const double rMoliere[kNumberCalorimeter]
Cos< T >::type cos(const T &t)
Abs< T >::type abs(const T &t)
Gflash3Vector getOrthogonalUnitVector()
void getGflashTrajectoryPoint(GflashTrajectoryPoint &point, double s) const
static const double tmax[3]
const double radLength[kNumberCalorimeter]
Gflash3Vector & getPosition()
GflashTrajectory * theHelix
Power< A, B >::type pow(const A &a, const B &b)
void initializeTrajectory(const HepGeom::Vector3D< double > &, const HepGeom::Point3D< double > &, double q, double Field)