92 double tempZCalo = 26;
93 double hfcriticalEnergy = 0.021;
95 std::vector<HFGflash::Hit>
hit;
98 G4StepPoint * preStepPoint = aStep->GetPreStepPoint();
100 G4Track * track = aStep->GetTrack();
102 const G4DynamicParticle *aParticle = track->GetDynamicParticle();
103 G4ThreeVector momDir = aParticle->GetMomentumDirection();
105 G4ThreeVector hitPoint = preStepPoint->GetPosition();
106 G4String partType = track->GetDefinition()->GetParticleName();
112 const G4double energyCutoff = 1;
113 const G4int maxNumberOfSpots = 10000000;
115 G4ThreeVector showerStartingPosition = track->GetPosition()/cm;
116 G4ThreeVector showerMomentum = preStepPoint->GetMomentum()/
GeV;
119 G4double logEinc =
std::log((preStepPoint->GetTotalEnergy())/
GeV);
121 G4double
y = ((preStepPoint->GetTotalEnergy())/
GeV) / hfcriticalEnergy;
125 G4double nSpots = 93.0 *
std::log(tempZCalo) * ((preStepPoint->GetTotalEnergy())/
GeV);
126 if(preStepPoint->GetTotalEnergy()/
GeV < 1.6) nSpots = 140.4 *
std::log(tempZCalo) *
std::pow(((preStepPoint->GetTotalEnergy())/
GeV),0.876);
130 double charge = track->GetStep()->GetPreStepPoint()->GetCharge();
134 G4double pathLength = pathLength0;
138 G4double fluctuatedTmax =
std::log(logY - 0.7157);
139 G4double fluctuatedAlpha=
std::log(0.7996 +(0.4581 + 1.8628/tempZCalo)*logY);
141 G4double sigmaTmax = 1.0/( -1.4 + 1.26 * logY);
142 G4double sigmaAlpha = 1.0/( -0.58 + 0.86 * logY);
143 G4double
rho = 0.705 - 0.023 * logY;
144 G4double sqrtPL =
std::sqrt((1.0+rho)/2.0);
145 G4double sqrtLE =
std::sqrt((1.0-rho)/2.0);
147 G4double norm1 = G4RandGauss::shoot();
148 G4double norm2 = G4RandGauss::shoot();
149 G4double tempTmax = fluctuatedTmax + sigmaTmax*(sqrtPL*norm1 + sqrtLE*norm2);
150 G4double tempAlpha = fluctuatedAlpha + sigmaAlpha*(sqrtPL*norm1 - sqrtLE*norm2);
157 if (!alpha)
return hit;
158 if (!beta)
return hit;
159 if (alpha < 0.00001)
return hit;
160 if (beta < 0.00001)
return hit;
163 G4double averageTmax = logY-0.858;
164 G4double averageAlpha = 0.21+(0.492+2.38/tempZCalo)*logY;
165 G4double spotTmax = averageTmax * (0.698 + .00212*tempZCalo);
166 G4double spotAlpha= averageAlpha * (0.639 + .00334*tempZCalo);
167 G4double spotBeta = (spotAlpha-1.0)/spotTmax;
169 if (!spotAlpha)
return hit;
170 if (!spotBeta)
return hit;
171 if (spotAlpha < 0.00001)
return hit;
172 if (spotBeta < 0.00001)
return hit;
175 LogDebug(
"HFShower") <<
"Incoming energy = " << ((preStepPoint->GetTotalEnergy())/
GeV) <<
" Position (rho,z) = (" << showerStartingPosition.rho() <<
", " << showerStartingPosition.z() <<
")";
178 em_incE->Fill(((preStepPoint->GetTotalEnergy())/
GeV));
179 em_ssp_rho->Fill(showerStartingPosition.rho());
184 G4double z1=0.0251+0.00319*logEinc;
185 G4double z2=0.1162-0.000381*tempZCalo;
187 G4double k1=0.659 - 0.00309 * tempZCalo;
190 G4double k4=0.3585+ 0.0421*logEinc;
192 G4double
p1=2.623 -0.00094*tempZCalo;
193 G4double
p2=0.401 +0.00187*tempZCalo;
194 G4double
p3=1.313 -0.0686*logEinc;
198 G4double e25Scale = 1.03551;
199 z1 *= 9.76972e-01 - 3.85026e-01 * std::tanh(1.82790
e+00*
std::log(((preStepPoint->GetTotalEnergy())/
GeV)) - 3.66237
e+00);
202 G4double stepLengthLeft = 10000;
204 G4double zInX0 = 0.0;
205 G4double deltaZInX0 = 0.0;
206 G4double deltaZ = 0.0;
207 G4double stepLengthLeftInX0 = 0.0;
209 const G4double divisionStepInX0 = 0.1;
210 G4double
energy = ((preStepPoint->GetTotalEnergy())/
GeV);
212 Genfun::IncompleteGamma gammaDist;
214 G4double energyInGamma = 0.0;
215 G4double preEnergyInGamma = 0.0;
216 G4double sigmaInGamma = 0.;
217 G4double preSigmaInGamma = 0.0;
220 G4double deltaEnergy =0.0 ;
221 G4int spotCounter = 0;
224 G4double deltaStep = 0.0;
227 G4double timeGlobal = track->GetStep()->GetPreStepPoint()->GetGlobalTime();
231 theGflashNavigator->SetWorldVolume(G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume());
236 LogDebug(
"HFShower") <<
" Energy = " << energy <<
" Step Length Left = " << stepLengthLeft;
238 while(energy > 0.0 && stepLengthLeft > 0.0) {
241 if ( stepLengthLeftInX0 < divisionStepInX0 ) {
242 deltaZInX0 = stepLengthLeftInX0;
244 stepLengthLeft = 0.0;
246 deltaZInX0 = divisionStepInX0;
248 stepLengthLeft -= deltaZ;
254 LogDebug(
"HFShower") <<
" zInX0 = " << zInX0 <<
" spotBeta*zInX0 = " << spotBeta*zInX0;
256 if ((!zInX0) || (!spotBeta*zInX0) || (zInX0 < 0.01) ||
257 (spotBeta*zInX0 < 0.00001) || (!zInX0*beta) || (zInX0*beta < 0.00001))
260 G4int nSpotsInStep = 0;
263 LogDebug(
"HFShower") <<
" Energy - Energy Cut off = " << energy - energyCutoff;
266 if ( energy > energyCutoff ) {
267 preEnergyInGamma = energyInGamma;
268 gammaDist.a().setValue(alpha);
270 energyInGamma = gammaDist(beta*zInX0);
271 G4double energyInDeltaZ = energyInGamma - preEnergyInGamma;
272 deltaEnergy =
std::min(energy,((preStepPoint->GetTotalEnergy())/
GeV)*energyInDeltaZ);
274 preSigmaInGamma = sigmaInGamma;
275 gammaDist.a().setValue(spotAlpha);
276 sigmaInGamma = gammaDist(spotBeta*zInX0);
277 nSpotsInStep =
std::max(1,
int(nSpots * (sigmaInGamma - preSigmaInGamma)));
280 preSigmaInGamma = sigmaInGamma;
281 nSpotsInStep =
std::max(1,
int(nSpots * (1.0 - preSigmaInGamma)));
284 if ( deltaEnergy > energy || (energy-deltaEnergy) < energyCutoff ) deltaEnergy =
energy;
286 energy -= deltaEnergy;
288 if ( spotCounter+nSpotsInStep > maxNumberOfSpots ) {
289 nSpotsInStep = maxNumberOfSpots - spotCounter;
290 if (nSpotsInStep < 1) nSpotsInStep = 1;
295 deltaStep += 0.5*deltaZ;
296 pathLength += deltaStep;
297 deltaStep = 0.5*deltaZ;
301 G4double tScale = tmax *alpha/(alpha-1.0) * (
std::exp(fluctuatedAlpha)-1.0)/
std::exp(fluctuatedAlpha);
302 G4double
tau =
std::min(10.0,(zInX0 - 0.5*deltaZInX0)/tScale);
303 G4double rCore = z1 + z2 *
tau;
305 G4double p23 = (p2 -
tau)/p3;
312 G4double emSpotEnergy = deltaEnergy / nSpotsInStep * e25Scale *
GeV;
316 LogDebug(
"HFShower") <<
" nSpotsInStep = " << nSpotsInStep;
321 for (G4int ispot = 0 ; ispot < nSpotsInStep ; ispot++) {
323 G4double u1 = G4UniformRand();
324 G4double u2 = G4UniformRand();
325 G4double rInRM = 0.0;
327 if (u1 < probabilityWeight) {
328 rInRM = rCore *
std::sqrt( u2/(1.0-u2) );
330 rInRM = rTail *
std::sqrt( u2/(1.0-u2) );
336 G4double azimuthalAngle = twopi*G4UniformRand();
340 G4double incrementPath = (deltaZ/nSpotsInStep)*(ispot+0.5 - 0.5*nSpotsInStep);
359 timeGlobal += 0.0001*nanosecond;
366 LogDebug(
"HFShower") <<
"zInX0_spot,emSpotEnergy/GeV =" << zInX0_spot <<
" , " << emSpotEnergy/GeV <<
"emSpotEnergy/GeV =" << emSpotEnergy/
GeV;
369 if ((!zInX0_spot) || (zInX0_spot < 0))
continue;
370 if ((!emSpotEnergy/GeV) || (emSpotEnergy < 0))
continue;
371 if ((!rShower/Gflash::rMoliere[
jCalorimeter]) || (rShower/Gflash::rMoliere[jCalorimeter] < 0))
continue;
377 em_long->Fill(SpotPosition0.z()/cm,emSpotEnergy/
GeV);
378 em_lateral->Fill(rShower/Gflash::rMoliere[jCalorimeter],emSpotEnergy/GeV);
386 double energyratio = emSpotEnergy/(preStepPoint->GetTotalEnergy()/(nSpots*e25Scale));
388 if (emSpotEnergy/GeV < 0.0001)
continue;
389 if (energyratio > 80)
continue;
392 if(SpotPosition0.z() > 0) zshift = 18;
393 if(SpotPosition0.z() < 0) zshift = -18;
396 G4ThreeVector gfshift(0,0,zshift*(
pow(100,0.1)/
pow(preStepPoint->GetTotalEnergy()/
GeV,0.1)));
398 G4ThreeVector SpotPosition = gfshift + SpotPosition0;
400 double LengthWeight = std::fabs(
std::pow(SpotPosition0.z()/11370,1));
401 if (G4UniformRand()> 0.0021 * energyratio * LengthWeight)
continue;
404 oneHit.
time = timeGlobal;
405 oneHit.
edep = emSpotEnergy/
GeV;
406 hit.push_back(oneHit);
G4Navigator * theGflashNavigator
static std::vector< std::string > checklist log
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)