122 double NCherPhot = 0.;
125 G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
126 G4VPhysicalVolume* currentPV = preStepPoint->GetPhysicalVolume();
127 const G4String& nameVolume = currentPV->GetName();
129 const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
130 const G4ThreeVector& hit_mom = preStepPoint->GetMomentumDirection();
131 G4double stepL = aStep->GetStepLength()/cm;
132 G4double
beta = preStepPoint->GetBeta();
133 G4double
charge = preStepPoint->GetCharge();
136 G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
137 G4VPhysicalVolume* postPV = postStepPoint->GetPhysicalVolume();
138 const G4String& postnameVolume = postPV->GetName();
141 G4Track* theTrack = aStep->GetTrack();
142 G4String
particleType = theTrack->GetDefinition()->GetParticleName();
143 const G4ThreeVector& vert_mom = theTrack->GetVertexMomentumDirection();
144 G4ThreeVector localPoint = theTrack->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
147 float costheta = vert_mom.z()/
sqrt(vert_mom.x()*vert_mom.x()+
148 vert_mom.y()*vert_mom.y()+
149 vert_mom.z()*vert_mom.z());
153 if (vert_mom.x() != 0) phi = std::atan2(vert_mom.y(),vert_mom.x());
154 if (phi < 0.) phi += twopi;
157 double stepE = aStep->GetTotalEnergyDeposit();
159 <<
"ZdcSD:: getEnergyDeposit: \n" 160 <<
" preStepPoint: " << nameVolume <<
"," << stepL <<
"," << stepE
161 <<
"," << beta <<
"," << charge <<
"\n" 162 <<
" postStepPoint: " << postnameVolume <<
"," << costheta <<
"," 163 << theta <<
"," << eta <<
"," << phi <<
"," << particleType <<
" id= " 164 << theTrack->GetTrackID() <<
" Etot(GeV)= " << theTrack->GetTotalEnergy()/
GeV;
166 const double bThreshold = 0.67;
167 if ((beta > bThreshold) && (charge != 0) && (nameVolume ==
"ZDC_EMFiber" || nameVolume ==
"ZDC_HadFiber")) {
168 LogDebug(
"ForwardSim") <<
"ZdcSD:: getEnergyDeposit: pass ";
170 const float nMedium = 1.4925;
174 const float photEnSpectrDE = 1.24;
180 const float effPMTandTransport = 0.15;
183 const float thFullRefl = 23.;
184 float thFullReflRad = thFullRefl*
pi/180.;
192 float costh = hit_mom.z()/
sqrt(hit_mom.x()*hit_mom.x()+
193 hit_mom.y()*hit_mom.y()+
194 hit_mom.z()*hit_mom.z());
197 if (th < 0.) th += twopi;
200 float costhcher =1./(nMedium*
beta);
204 float DelFibPart =
std::abs(th - thFibDirRad);
217 if (DelFibPart > (thFullReflRad + thcher) ) {
218 variant = 0.; d_qz = 0.;
221 if ((th + thcher) < (thFibDirRad+thFullReflRad) && (th - thcher) > (thFibDirRad-thFullReflRad) ) {
222 variant = 1.; d_qz = 1.;
225 if ((thFibDirRad + thFullReflRad) < (th + thcher) && (thFibDirRad - thFullReflRad) > (th - thcher) ) {
226 variant = 2.; d_qz = 0.;
231 float arg_arcos = 0.;
232 float tan_arcos = 2.*a*
d;
233 if (tan_arcos != 0.) arg_arcos =(r*r-a*a-d*
d)/tan_arcos;
239 d_qz = th_arcos/twopi;
246 double meanNCherPhot = 0.;
247 int poissNCherPhot = 0;
249 meanNCherPhot = 370.*charge*charge*( 1. - 1./(nMedium*nMedium*beta*
beta) ) * photEnSpectrDE * stepL;
251 poissNCherPhot =
std::max((
int)G4Poisson(meanNCherPhot),0);
252 NCherPhot = poissNCherPhot * effPMTandTransport * d_qz;
256 <<
"ZdcSD:: getEnergyDeposit: gED: " 275 <<
"," << meanNCherPhot
276 <<
"," << poissNCherPhot
294 if (beta <= bThreshold)
296 <<
"ZdcSD:: getEnergyDeposit: fail beta=" <<
beta;
299 <<
"ZdcSD:: getEnergyDeposit: fail charge=0";
300 if ( !(nameVolume ==
"ZDC_EMFiber" || nameVolume ==
"ZDC_HadFiber") )
302 <<
"ZdcSD:: getEnergyDeposit: fail nv=" << nameVolume;
Geom::Theta< T > theta() const
Tan< T >::type tan(const T &t)
Abs< T >::type abs(const T &t)