121 double NCherPhot = 0.;
124 G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
125 G4VPhysicalVolume* currentPV = preStepPoint->GetPhysicalVolume();
128 const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
129 const G4ThreeVector& hit_mom = preStepPoint->GetMomentumDirection();
130 G4double stepL = aStep->GetStepLength() / cm;
131 G4double
beta = preStepPoint->GetBeta();
132 G4double
charge = preStepPoint->GetCharge();
135 G4Track* theTrack = aStep->GetTrack();
136 G4String
particleType = theTrack->GetDefinition()->GetParticleName();
137 G4ThreeVector localPoint = theTrack->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
140 const G4ThreeVector& vert_mom = theTrack->GetVertexMomentumDirection();
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();
157 G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
158 G4VPhysicalVolume* postPV = postStepPoint->GetPhysicalVolume();
161 <<
" preStepPoint: " << nameVolume <<
"," << stepL <<
"," << stepE <<
"," << beta
162 <<
"," << charge <<
"\n"
163 <<
" postStepPoint: " << postnameVolume <<
"," << costheta <<
"," << theta <<
","
164 << eta <<
"," << phi <<
"," << particleType <<
" id= " << theTrack->GetTrackID()
165 <<
" Etot(GeV)= " << theTrack->GetTotalEnergy() / GeV;
167 const double bThreshold = 0.67;
168 if ((beta > bThreshold) && (charge != 0) && (nameVolume ==
"ZDC_EMFiber" || nameVolume ==
"ZDC_HadFiber")) {
172 const float nMedium = 1.4925;
176 const float photEnSpectrDE = 1.24;
182 const float effPMTandTransport = 0.15;
185 const float thFullRefl = 23.;
186 float thFullReflRad = thFullRefl *
pi / 180.;
194 float costh = hit_mom.z() /
sqrt(hit_mom.x() * hit_mom.x() + hit_mom.y() * hit_mom.y() + hit_mom.z() * hit_mom.z());
201 float costhcher = 1. / (nMedium *
beta);
205 float DelFibPart =
std::abs(th - thFibDirRad);
219 if (DelFibPart > (thFullReflRad + thcher)) {
226 if ((th + thcher) < (thFibDirRad + thFullReflRad) && (th - thcher) > (thFibDirRad - thFullReflRad)) {
233 if ((thFibDirRad + thFullReflRad) < (th + thcher) && (thFibDirRad - thFullReflRad) > (th - thcher)) {
243 float arg_arcos = 0.;
244 float tan_arcos = 2. * a *
d;
246 arg_arcos = (r * r - a * a - d *
d) / tan_arcos;
249 d_qz = th_arcos / twopi;
252 edm::LogVerbatim(
"ForwardSim") <<
" d_qz: " << r <<
"," << a <<
"," << d <<
" " << tan_arcos <<
" "
262 double meanNCherPhot = 0.;
263 int poissNCherPhot = 0;
265 meanNCherPhot = 370. * charge * charge * (1. - 1. / (nMedium * nMedium * beta *
beta)) * photEnSpectrDE * stepL;
267 poissNCherPhot =
std::max((
int)G4Poisson(meanNCherPhot), 0);
268 NCherPhot = poissNCherPhot * effPMTandTransport * d_qz;
272 edm::LogVerbatim(
"ForwardSim") <<
"ZdcSD:: getEnergyDeposit: gED: " << stepE <<
"," << costh <<
"," << th <<
","
273 << costhcher <<
"," << thcher <<
"," << DelFibPart <<
"," << d <<
"," << a <<
","
274 << r <<
"," << hitPoint <<
"," << hit_mom <<
"," << vert_mom <<
"," << localPoint
275 <<
"," << charge <<
"," << beta <<
"," << stepL <<
"," << d_qz <<
"," << variant
276 <<
"," << meanNCherPhot <<
"," << poissNCherPhot <<
"," << NCherPhot;
294 if (beta <= bThreshold)
297 edm::LogVerbatim(
"ForwardSim") <<
"ZdcSD:: getEnergyDeposit: fail charge=0";
298 if (!(nameVolume ==
"ZDC_EMFiber" || nameVolume ==
"ZDC_HadFiber"))
299 edm::LogVerbatim(
"ForwardSim") <<
"ZdcSD:: getEnergyDeposit: fail nv=" << nameVolume;
Log< level::Info, true > LogVerbatim
static std::vector< std::string > checklist log
Geom::Theta< T > theta() const
Tan< T >::type tan(const T &t)
Abs< T >::type abs(const T &t)
std::string getName(const G4String &)