118 double NCherPhot = 0.;
121 G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
122 G4VPhysicalVolume* currentPV = preStepPoint->GetPhysicalVolume();
123 const G4String& nameVolume = currentPV->GetName();
125 const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
126 const G4ThreeVector& hit_mom = preStepPoint->GetMomentumDirection();
127 G4double stepL = aStep->GetStepLength() / cm;
128 G4double
beta = preStepPoint->GetBeta();
129 G4double
charge = preStepPoint->GetCharge();
132 G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
133 G4VPhysicalVolume* postPV = postStepPoint->GetPhysicalVolume();
134 const G4String& postnameVolume = postPV->GetName();
137 G4Track* theTrack = aStep->GetTrack();
138 G4String
particleType = theTrack->GetDefinition()->GetParticleName();
139 const G4ThreeVector& vert_mom = theTrack->GetVertexMomentumDirection();
140 G4ThreeVector localPoint = theTrack->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
144 vert_mom.z() /
sqrt(vert_mom.x() * vert_mom.x() + vert_mom.y() * vert_mom.y() + vert_mom.z() * vert_mom.z());
148 if (vert_mom.x() != 0)
149 phi = std::atan2(vert_mom.y(), vert_mom.x());
154 double stepE = aStep->GetTotalEnergyDeposit();
155 LogDebug(
"ForwardSim") <<
"ZdcSD:: getEnergyDeposit: \n" 156 <<
" preStepPoint: " << nameVolume <<
"," << stepL <<
"," << stepE <<
"," << beta <<
"," 158 <<
" postStepPoint: " << postnameVolume <<
"," << costheta <<
"," << theta <<
"," << eta
159 <<
"," << phi <<
"," << particleType <<
" id= " << theTrack->GetTrackID()
160 <<
" Etot(GeV)= " << theTrack->GetTotalEnergy() /
GeV;
162 const double bThreshold = 0.67;
163 if ((beta > bThreshold) && (charge != 0) && (nameVolume ==
"ZDC_EMFiber" || nameVolume ==
"ZDC_HadFiber")) {
164 LogDebug(
"ForwardSim") <<
"ZdcSD:: getEnergyDeposit: pass ";
166 const float nMedium = 1.4925;
170 const float photEnSpectrDE = 1.24;
176 const float effPMTandTransport = 0.15;
179 const float thFullRefl = 23.;
180 float thFullReflRad = thFullRefl *
pi / 180.;
188 float costh = hit_mom.z() /
sqrt(hit_mom.x() * hit_mom.x() + hit_mom.y() * hit_mom.y() + hit_mom.z() * hit_mom.z());
195 float costhcher = 1. / (nMedium *
beta);
199 float DelFibPart =
std::abs(th - thFibDirRad);
212 if (DelFibPart > (thFullReflRad + thcher)) {
217 if ((th + thcher) < (thFibDirRad + thFullReflRad) && (th - thcher) > (thFibDirRad - thFullReflRad)) {
222 if ((thFibDirRad + thFullReflRad) < (th + thcher) && (thFibDirRad - thFullReflRad) > (th - thcher)) {
229 float arg_arcos = 0.;
230 float tan_arcos = 2. * a *
d;
232 arg_arcos = (r * r - a * a - d *
d) / tan_arcos;
238 d_qz = th_arcos / twopi;
245 double meanNCherPhot = 0.;
246 int poissNCherPhot = 0;
248 meanNCherPhot = 370. * charge * charge * (1. - 1. / (nMedium * nMedium * beta *
beta)) * photEnSpectrDE * stepL;
250 poissNCherPhot =
std::max((
int)G4Poisson(meanNCherPhot), 0);
251 NCherPhot = poissNCherPhot * effPMTandTransport * d_qz;
254 LogDebug(
"ForwardSim") <<
"ZdcSD:: getEnergyDeposit: gED: " << stepE <<
"," << costh <<
"," << th <<
"," 255 << costhcher <<
"," << thcher <<
"," << DelFibPart <<
"," << d <<
"," << a <<
"," << r <<
"," 256 << hitPoint <<
"," << hit_mom <<
"," << vert_mom <<
"," << localPoint <<
"," << charge <<
"," 257 << beta <<
"," << stepL <<
"," << d_qz <<
"," << variant <<
"," << meanNCherPhot <<
"," 258 << poissNCherPhot <<
"," << NCherPhot;
275 if (beta <= bThreshold)
276 LogDebug(
"ForwardSim") <<
"ZdcSD:: getEnergyDeposit: fail beta=" <<
beta;
278 LogDebug(
"ForwardSim") <<
"ZdcSD:: getEnergyDeposit: fail charge=0";
279 if (!(nameVolume ==
"ZDC_EMFiber" || nameVolume ==
"ZDC_HadFiber"))
280 LogDebug(
"ForwardSim") <<
"ZdcSD:: getEnergyDeposit: fail nv=" << nameVolume;
Geom::Theta< T > theta() const
Tan< T >::type tan(const T &t)
Abs< T >::type abs(const T &t)