66{
67 if(verboseLevel>5)
68 {
69 G4cout << "AFP_TDSensitiveDetector::ProcessHits" << G4endl;
70 }
71
72 bool bRes=false;
73
74 int nTrackID=-1;
75 int nParticleEncoding=-1;
76 float fKineticEnergy=0.0;
77 float fEnergyDeposit=0.0;
78 float fWaveLength=0.0;
79 float fPreStepX=0.0;
80 float fPreStepY=0.0;
81 float fPreStepZ=0.0;
82 float fPostStepX=0.0;
83 float fPostStepY=0.0;
84 float fPostStepZ=0.0;
85 float fGlobalTime=0.0;
86 int nStationID=-1;
87 int nDetectorID=-1;
88 int nQuarticID=-1;
89
90
91
92
93 G4Track* pTrack = pStep->GetTrack();
94 G4ParticleDefinition* pParticleDefinition = pTrack->GetDefinition();
95 G4StepPoint* pPreStepPoint = pStep->GetPreStepPoint();
96 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
97 G4ThreeVector PreStepPointPos = pPreStepPoint->GetPosition();
98 G4ThreeVector PostStepPointPos = pPostStepPoint->GetPosition();
99
100 nTrackID=pTrack->GetTrackID();
101 fKineticEnergy = pPreStepPoint->GetKineticEnergy();
102 fEnergyDeposit = pStep->GetTotalEnergyDeposit();
103
104 fPreStepX = PreStepPointPos.x();
105 fPreStepY = PreStepPointPos.y();
106 fPreStepZ = PreStepPointPos.z();
107 fPostStepX = PostStepPointPos.x();
108 fPostStepY = PostStepPointPos.y();
109 fPostStepZ = PostStepPointPos.z();
110 nParticleEncoding = pParticleDefinition->GetPDGEncoding();
111 fGlobalTime = pStep->GetPreStepPoint()->GetGlobalTime()/CLHEP::picosecond;
112
113
114 G4TouchableHandle touch1 = pPreStepPoint->GetTouchableHandle();
115 G4VPhysicalVolume* volume = touch1->GetVolume();
116 G4String VolumeName = volume->GetName();
117
118
119
120
121
122 if(verboseLevel>5)
123 {
124 G4cout << "hit volume name is " << VolumeName << G4endl;
125
126 G4cout << "global, x_pre: " << fPreStepX << ", y_pre: " << fPreStepY << ", z_pre: " << fPreStepZ << G4endl;
127 G4cout << "global, x_post: " << fPostStepX << ", y_post: " << fPostStepY << ", z_post: " << fPostStepZ << G4endl;
128 }
129
130 char* ppv1, *ppv2;
131 char szbuff[32];
132 memset(&szbuff[0],0,sizeof(szbuff));
133 strncpy(szbuff,VolumeName.data(),sizeof(szbuff));
134 szbuff[sizeof(szbuff)-1] = '\0';
135 ppv1=strchr(szbuff,'[');
136 ppv2=strchr(szbuff,']');
137 if(!ppv2 || !ppv1){
138 G4cout << "ERROR: Invalid format of volume name " << VolumeName << G4endl;
139 return false;
140 }
141 else *ppv2='\0';
142
143 nStationID=10*(szbuff[3]-0x30)+(szbuff[4]-0x30);
144 nDetectorID=
atoi(ppv1+1);
145
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186#if G4VERSION_NUMBER < 1100
187 if ( (VolumeName.contains("TDQuarticBarVacBorder")) && pParticleDefinition->GetPDGCharge() !=0 )
188#else
189 if ( (G4StrUtil::contains(VolumeName, "TDQuarticBarVacBorder")) && pParticleDefinition->GetPDGCharge() !=0 )
190#endif
191 {
192 nQuarticID=szbuff[7]-0x30;
193
194
195
196
197
198
199
200
201
202 }
203
205#if G4VERSION_NUMBER < 1100
206 if ( (bRes=VolumeName.contains("TDQuarticBar[")) )
207#else
208 if ( (bRes=G4StrUtil::contains(VolumeName, "TDQuarticBar[")) )
209#endif
210 {
211 nQuarticID=szbuff[7]-0x30;
212
213
218
219
220 const G4TouchableHistory* myTouch = static_cast<const G4TouchableHistory*>(pPreStepPoint->GetTouchable());
221
222
223
224 const G4AffineTransform transformation = myTouch->GetHistory()->GetTopTransform();
225 G4ThreeVector PreStepPointPos2 = transformation.TransformPoint(PreStepPointPos);
226 G4ThreeVector PostStepPointPos2 = transformation.TransformPoint(PostStepPointPos);
227
228 G4String shape( myTouch->GetSolid()->GetEntityType() );
229
230 G4ThreeVector normpX( 1., 0., 0.);
231 G4ThreeVector normnX(-1., 0., 0.);
232 G4ThreeVector normpY( 0., 1., 0.);
233 G4ThreeVector normnY( 0.,-1., 0.);
234 G4ThreeVector normpZ( 0., 0., 1.);
235 G4ThreeVector normnZ( 0., 0.,-1.);
236
237
238
239
240
241
242
243
244
245
246
247
248 G4double PreProtonX = PreStepPointPos2.x();
249 G4double PreProtonY = PreStepPointPos2.y();
250 G4double PreProtonZ = PreStepPointPos2.z();
251
252 G4double PostProtonX = PostStepPointPos2.x();
253 G4double PostProtonY = PostStepPointPos2.y();
254 G4double PostProtonZ = PostStepPointPos2.z();
255
256 G4ThreeVector
p0 = pStep->GetDeltaPosition().unit();
257
258 G4Material*
mat = pStep->GetTrack()->GetMaterial();
259 G4MaterialPropertiesTable *matPropTable =
mat->GetMaterialPropertiesTable();
260
261 G4MaterialPropertyVector* Rind = matPropTable->GetProperty("RINDEX");
262
263 G4MaterialPropertyVector* Alen = matPropTable->GetProperty("ABSLENGTH");
264
265 const G4double
charge = pParticleDefinition->GetPDGCharge();
266 const G4double
beta = (pPreStepPoint->GetBeta() + pPostStepPoint->GetBeta()) / 2.;
267 G4double BetaInverse = 1. /
beta;
268
269
270
271#if G4VERSION_NUMBER < 1100
272 G4double Pmin = Rind->GetMinLowEdgeEnergy();
273
274 G4double Pmax = Rind->GetMaxLowEdgeEnergy();
275#else
276 G4double Pmin = Rind->GetMinEnergy();
277 G4double Pmax = Rind->GetMaxEnergy();
278#endif
279 G4double
dp = Pmax - Pmin;
280
281 G4double maxCosTheta = BetaInverse / Rind->GetMinValue();
282 G4double maxSin2Theta = (1.0 - maxCosTheta) * (1.0 + maxCosTheta);
283
284
285 G4double meanRI = .5*(Rind->GetMinValue() + Rind->GetMaxValue());
286
287
288 G4double MeanNumberOfPhotons = 370.*(
charge/CLHEP::eplus)*(
charge/CLHEP::eplus)* (1.0 - 1.0/(
beta * meanRI *
beta * meanRI)) / (CLHEP::cm*CLHEP::eV);
289 if (MeanNumberOfPhotons <= 0.0) return 1;
290
291 G4double step_length = pStep->GetStepLength();
292
293 MeanNumberOfPhotons = MeanNumberOfPhotons * step_length *
dp;
294 G4int NumPhotons = (G4int) G4Poisson( MeanNumberOfPhotons );
295 if (NumPhotons <= 0) return 1;
296
297
298 G4int NumPhotonsCuts=0;
299 for (G4int
I = 0;
I < NumPhotons;
I++) {
300
302 G4double sampledEnergy, sampledRI;
303 G4double cosTheta, sin2Theta;
304
305
306 do {
307 rand = G4UniformRand();
308 sampledEnergy = Pmin +
rand *
dp;
309
310 sampledRI = Rind->Value(sampledEnergy);
311 cosTheta = BetaInverse / sampledRI;
312
313 sin2Theta = (1.0 - cosTheta)*(1.0 + cosTheta);
314 rand = G4UniformRand();
315
316 } while (rand * maxSin2Theta > sin2Theta);
317
318
319 rand = G4UniformRand();
321 G4double sinPhi =
sin(
phi);
322 G4double cosPhi =
cos(
phi);
323
324
325
326
327 G4double sinTheta = sqrt(sin2Theta);
328 G4double
px = sinTheta*cosPhi;
329 G4double
py = sinTheta*sinPhi;
330 G4double
pz = cosTheta;
331 G4ParticleMomentum photonMomentum(px, py, pz);
332 photonMomentum.rotateUz(p0);
333
334 G4double PX = photonMomentum.getX();
335 G4double PY = photonMomentum.getY();
336 G4double PZ = photonMomentum.getZ();
337
338
339
340 G4double PYp = PY/sqrt(PX*PX+PY*PY+PZ*PZ);
341 G4double PZp = PZ/sqrt(PX*PX+PY*PY+PZ*PZ);
342
343 G4double PYt = PY/sqrt(PY*PY+PZ*PZ);
344 G4double PZt = PZ/sqrt(PY*PY+PZ*PZ);
345
346
347 cosPhi = (PYp*PYt + PZp*PZt);
348 if (nStationID == 0) cosTheta = ( -PYt*
sin(48.*CLHEP::deg) + PZt*
cos(48.*CLHEP::deg) );
349 else cosTheta = ( -PYt*
sin(48.*CLHEP::deg) - PZt*
cos(48.*CLHEP::deg) );
350
351
352 G4double cosThetaC = sqrt(1.-1./sampledRI/sampledRI);
353 if (sqrt(1.-cosPhi*cosPhi)>cosThetaC) continue;
354 if (sqrt(1.-cosTheta*cosTheta)*cosPhi>cosThetaC) continue;
355
356
357 rand = G4UniformRand();
358 G4double PhotonX = PreProtonX + (PostProtonX-PreProtonX)*rand;
359 G4double PhotonY = PreProtonY + (PostProtonY-PreProtonY)*rand;
360 G4double PhotonZ = PreProtonZ + (PostProtonZ-PreProtonZ)*rand;
361 G4ThreeVector PhotonPos(PhotonX,PhotonY,PhotonZ);
362 G4double PhotonR = sqrt( (PreProtonX-PhotonX)*(PreProtonX-PhotonX) + (PreProtonY-PhotonY)*(PreProtonY-PhotonY) + (PreProtonZ-PhotonZ)*(PreProtonZ-PhotonZ) );
363
364 G4double Y0;
365
366 if (cosTheta>=0)
367 {
368 Y0 = static_cast<G4ReflectedSolid *>(myTouch->GetSolid())->DistanceToOut(PhotonPos, normnY);
369 Y0 = Y0/cosTheta/cosPhi;
370 }
371 else
372 {
373 continue;
374
375
376 }
377
378
379
380 float Pabs = 1. -
exp( - Y0/Alen->Value(sampledEnergy) );
381 rand = G4UniformRand();
382 if (Pabs>rand) continue;
383
384
385 rand = G4UniformRand();
388
389 NumPhotonsCuts++;
390
391 float fGlobalTime2 = fGlobalTime;
392 fGlobalTime2 += ( PhotonR * BetaInverse / CLHEP::c_light )/CLHEP::picosecond;
393
394
395 float EdndE;
396
397 if (sampledEnergy > (Pmin+.5*dp)) EdndE = (sampledRI - Rind->Value(sampledEnergy-0.0001*CLHEP::eV))/0.0001*sampledEnergy/CLHEP::eV;
398
399 else EdndE = (Rind->Value(sampledEnergy+0.0001*CLHEP::eV) - sampledRI)/0.0001*sampledEnergy/CLHEP::eV;
400 fGlobalTime2 += ( (sampledRI + EdndE)* Y0 * BetaInverse / CLHEP::c_light )/CLHEP::picosecond;
401
402 if (verboseLevel>5)
403 {
404 G4cout << "FastCher EdndE: " << EdndE << G4endl;
405 }
406 fWaveLength = 2.*
M_PI*CLHEP::hbarc/sampledEnergy/(CLHEP::MeV*CLHEP::nm);
407
408
413
414 int nSensitiveElementID=-1;
415 if(nQuarticID==0) { nSensitiveElementID=1; }
416 else if(nQuarticID==1) { nSensitiveElementID=3; }
417
418 m_HitColl->Emplace(
m_nHitID,nTrackID,nParticleEncoding,fKineticEnergy,fEnergyDeposit,
419 fWaveLength,PhotonX,PhotonY,PhotonZ,(PhotonX+PX),(PhotonY+PY),(PhotonZ+PZ),
420 fGlobalTime2,nStationID,nDetectorID,nSensitiveElementID);
422
424 else if(nStationID==0 && nQuarticID==1)
m_nNOfTDSimHits[1][nDetectorID]++;
425 else if(nStationID==3 && nQuarticID==0)
m_nNOfTDSimHits[2][nDetectorID]++;
426 else if(nStationID==3 && nQuarticID==1)
m_nNOfTDSimHits[3][nDetectorID]++;
427 }
428 if(verboseLevel>5)
429 {
430 G4cout << "FastCher number of photons: " << NumPhotonsCuts << G4endl;
431 }
432 }
434 return true;
435}
Scalar phi() const
phi method
double charge(const T &p)
static constexpr int TDMaxCnt
static constexpr double TDMaxQEff
int atoi(std::string_view str)
Helper functions to unpack numbers decoded in string into integers and doubles The strings are requir...