89 double tempZCalo = 26;
90 double hfcriticalEnergy = 0.021;
92 std::vector<HFGflash::Hit>
hit;
95 G4StepPoint * preStepPoint = aStep->GetPreStepPoint();
97 G4Track * track = aStep->GetTrack();
99 const G4DynamicParticle *aParticle = track->GetDynamicParticle();
100 G4ThreeVector momDir = aParticle->GetMomentumDirection();
102 G4ThreeVector hitPoint = preStepPoint->GetPosition();
103 G4String partType = track->GetDefinition()->GetParticleName();
109 const G4double energyCutoff = 1;
110 const G4int maxNumberOfSpots = 10000000;
112 G4ThreeVector showerStartingPosition = track->GetPosition()/cm;
113 G4ThreeVector showerMomentum = preStepPoint->GetMomentum()/GeV;
116 G4double logEinc =
std::log((preStepPoint->GetTotalEnergy())/GeV);
118 G4double
y = ((preStepPoint->GetTotalEnergy())/GeV) / hfcriticalEnergy;
122 G4double nSpots = 93.0 *
std::log(tempZCalo) * ((preStepPoint->GetTotalEnergy())/GeV);
123 if(preStepPoint->GetTotalEnergy()/GeV < 1.6) nSpots = 140.4 *
std::log(tempZCalo) *
std::pow(((preStepPoint->GetTotalEnergy())/GeV),0.876);
127 double charge = track->GetStep()->GetPreStepPoint()->GetCharge();
131 G4double pathLength = pathLength0;
135 G4double fluctuatedTmax =
std::log(logY - 0.7157);
136 G4double fluctuatedAlpha=
std::log(0.7996 +(0.4581 + 1.8628/tempZCalo)*logY);
138 G4double sigmaTmax = 1.0/( -1.4 + 1.26 * logY);
139 G4double sigmaAlpha = 1.0/( -0.58 + 0.86 * logY);
140 G4double
rho = 0.705 - 0.023 * logY;
141 G4double sqrtPL =
std::sqrt((1.0+rho)/2.0);
142 G4double sqrtLE =
std::sqrt((1.0-rho)/2.0);
144 G4double norm1 = G4RandGauss::shoot();
145 G4double norm2 = G4RandGauss::shoot();
146 G4double tempTmax = fluctuatedTmax + sigmaTmax*(sqrtPL*norm1 + sqrtLE*norm2);
147 G4double tempAlpha = fluctuatedAlpha + sigmaAlpha*(sqrtPL*norm1 - sqrtLE*norm2);
154 if (!alpha)
return hit;
155 if (!beta)
return hit;
156 if (alpha < 0.00001)
return hit;
157 if (beta < 0.00001)
return hit;
160 G4double averageTmax = logY-0.858;
161 G4double averageAlpha = 0.21+(0.492+2.38/tempZCalo)*logY;
162 G4double spotTmax = averageTmax * (0.698 + .00212*tempZCalo);
163 G4double spotAlpha= averageAlpha * (0.639 + .00334*tempZCalo);
164 G4double spotBeta = (spotAlpha-1.0)/spotTmax;
166 if (!spotAlpha)
return hit;
167 if (!spotBeta)
return hit;
168 if (spotAlpha < 0.00001)
return hit;
169 if (spotBeta < 0.00001)
return hit;
172 LogDebug(
"HFShower") <<
"Incoming energy = " << ((preStepPoint->GetTotalEnergy())/GeV) <<
" Position (rho,z) = (" << showerStartingPosition.rho() <<
", " << showerStartingPosition.z() <<
")";
175 em_incE->Fill(((preStepPoint->GetTotalEnergy())/GeV));
176 em_ssp_rho->Fill(showerStartingPosition.rho());
181 G4double z1=0.0251+0.00319*logEinc;
182 G4double z2=0.1162-0.000381*tempZCalo;
184 G4double k1=0.659 - 0.00309 * tempZCalo;
187 G4double k4=0.3585+ 0.0421*logEinc;
189 G4double
p1=2.623 -0.00094*tempZCalo;
190 G4double
p2=0.401 +0.00187*tempZCalo;
191 G4double
p3=1.313 -0.0686*logEinc;
195 G4double e25Scale = 1.03551;
196 z1 *= 9.76972e-01 - 3.85026e-01 * std::tanh(1.82790
e+00*
std::log(((preStepPoint->GetTotalEnergy())/GeV)) - 3.66237
e+00);
199 G4double stepLengthLeft = 10000;
201 G4double zInX0 = 0.0;
202 G4double deltaZInX0 = 0.0;
203 G4double deltaZ = 0.0;
204 G4double stepLengthLeftInX0 = 0.0;
206 const G4double divisionStepInX0 = 0.1;
207 G4double
energy = ((preStepPoint->GetTotalEnergy())/GeV);
209 Genfun::IncompleteGamma gammaDist;
211 G4double energyInGamma = 0.0;
212 G4double preEnergyInGamma = 0.0;
213 G4double sigmaInGamma = 0.;
214 G4double preSigmaInGamma = 0.0;
217 G4double deltaEnergy =0.0 ;
218 G4int spotCounter = 0;
221 G4double deltaStep = 0.0;
224 G4double timeGlobal = track->GetStep()->GetPreStepPoint()->GetGlobalTime();
228 theGflashNavigator->SetWorldVolume(G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume());
233 LogDebug(
"HFShower") <<
" Energy = " << energy <<
" Step Length Left = " << stepLengthLeft;
235 while(energy > 0.0 && stepLengthLeft > 0.0) {
238 if ( stepLengthLeftInX0 < divisionStepInX0 ) {
239 deltaZInX0 = stepLengthLeftInX0;
241 stepLengthLeft = 0.0;
243 deltaZInX0 = divisionStepInX0;
245 stepLengthLeft -= deltaZ;
251 LogDebug(
"HFShower") <<
" zInX0 = " << zInX0 <<
" spotBeta*zInX0 = " << spotBeta*zInX0;
253 if ((!zInX0) || (!spotBeta*zInX0) || (zInX0 < 0.01) ||
254 (spotBeta*zInX0 < 0.00001) || (!zInX0*beta) || (zInX0*beta < 0.00001))
257 G4int nSpotsInStep = 0;
260 LogDebug(
"HFShower") <<
" Energy - Energy Cut off = " << energy - energyCutoff;
263 if ( energy > energyCutoff ) {
264 preEnergyInGamma = energyInGamma;
265 gammaDist.a().setValue(alpha);
267 energyInGamma = gammaDist(beta*zInX0);
268 G4double energyInDeltaZ = energyInGamma - preEnergyInGamma;
269 deltaEnergy =
std::min(energy,((preStepPoint->GetTotalEnergy())/GeV)*energyInDeltaZ);
271 preSigmaInGamma = sigmaInGamma;
272 gammaDist.a().setValue(spotAlpha);
273 sigmaInGamma = gammaDist(spotBeta*zInX0);
274 nSpotsInStep =
std::max(1,
int(nSpots * (sigmaInGamma - preSigmaInGamma)));
277 preSigmaInGamma = sigmaInGamma;
278 nSpotsInStep =
std::max(1,
int(nSpots * (1.0 - preSigmaInGamma)));
281 if ( deltaEnergy > energy || (energy-deltaEnergy) < energyCutoff ) deltaEnergy =
energy;
283 energy -= deltaEnergy;
285 if ( spotCounter+nSpotsInStep > maxNumberOfSpots ) {
286 nSpotsInStep = maxNumberOfSpots - spotCounter;
287 if (nSpotsInStep < 1) nSpotsInStep = 1;
292 deltaStep += 0.5*deltaZ;
293 pathLength += deltaStep;
294 deltaStep = 0.5*deltaZ;
298 G4double tScale = tmax *alpha/(alpha-1.0) * (
std::exp(fluctuatedAlpha)-1.0)/
std::exp(fluctuatedAlpha);
299 G4double
tau =
std::min(10.0,(zInX0 - 0.5*deltaZInX0)/tScale);
300 G4double rCore = z1 + z2 *
tau;
302 G4double p23 = (p2 -
tau)/p3;
309 G4double emSpotEnergy = deltaEnergy / nSpotsInStep * e25Scale * GeV;
313 LogDebug(
"HFShower") <<
" nSpotsInStep = " << nSpotsInStep;
318 for (G4int ispot = 0 ; ispot < nSpotsInStep ; ispot++) {
320 G4double u1 = G4UniformRand();
321 G4double u2 = G4UniformRand();
322 G4double rInRM = 0.0;
324 if (u1 < probabilityWeight) {
325 rInRM = rCore *
std::sqrt( u2/(1.0-u2) );
327 rInRM = rTail *
std::sqrt( u2/(1.0-u2) );
333 G4double azimuthalAngle = twopi*G4UniformRand();
337 G4double incrementPath = (deltaZ/nSpotsInStep)*(ispot+0.5 - 0.5*nSpotsInStep);
355 timeGlobal += 0.0001*nanosecond;
362 LogDebug(
"HFShower") <<
"zInX0_spot,emSpotEnergy/GeV =" << zInX0_spot <<
" , " << emSpotEnergy/GeV <<
"emSpotEnergy/GeV =" << emSpotEnergy/GeV;
365 if ((!zInX0_spot) || (zInX0_spot < 0))
continue;
366 if ((!emSpotEnergy/GeV) || (emSpotEnergy < 0))
continue;
367 if ((!rShower/Gflash::rMoliere[
jCalorimeter]) || (rShower/Gflash::rMoliere[jCalorimeter] < 0))
continue;
373 em_long->Fill(SpotPosition0.z()/cm,emSpotEnergy/GeV);
374 em_lateral->Fill(rShower/Gflash::rMoliere[jCalorimeter],emSpotEnergy/GeV);
375 em_2d->Fill(SpotPosition0.z()/cm,rShower/Gflash::rMoliere[
jCalorimeter],emSpotEnergy/GeV);
379 if(SpotPosition0 == 0)
continue;
381 double energyratio = emSpotEnergy/(preStepPoint->GetTotalEnergy()/(nSpots*e25Scale));
383 if (emSpotEnergy/GeV < 0.0001)
continue;
384 if (energyratio > 80)
continue;
387 if(SpotPosition0.z() > 0) zshift = 18;
388 if(SpotPosition0.z() < 0) zshift = -18;
391 G4ThreeVector gfshift(0,0,zshift*(
pow(100,0.1)/
pow(preStepPoint->GetTotalEnergy()/GeV,0.1)));
393 G4ThreeVector SpotPosition = gfshift + SpotPosition0;
395 double LengthWeight = std::fabs(
std::pow(SpotPosition0.z()/11370,1));
396 if (G4UniformRand()> 0.0021 * energyratio * LengthWeight)
continue;
399 oneHit.
time = timeGlobal;
400 oneHit.
edep = emSpotEnergy/GeV;
401 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]
const T & max(const T &a, const T &b)
Cos< T >::type cos(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)