ATLAS Offline Software
Loading...
Searching...
No Matches
AFP_SensitiveDetector.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3*/
4
5// Class header
6#include "AFP_SensitiveDetector.h"
7
8// Athena headers
9#include "AFP_Geometry/AFP_constants.h"
10#include "HitManagement/HitCollectionMap.h"
11#include "MCTruth/AtlasG4EventUserInfo.h"
12
13// Geant4 headers
14#include "G4Version.hh"
15#include "G4TouchableHistory.hh"
16#include "G4Step.hh"
17#include "G4Track.hh"
18#include "G4ParticleDefinition.hh"
19#include "G4StepPoint.hh"
20#include "G4ThreeVector.hh"
21#include "G4Poisson.hh"
22#include "G4VSolid.hh"
23#include "G4VProcess.hh"
24#include "G4ReflectedSolid.hh"
25#include "G4Material.hh"
26#include "G4MaterialPropertyVector.hh"
27
28// STL header
29#include <sstream>
30
31//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32
33AFP_SensitiveDetector::AFP_SensitiveDetector(const std::string& name, const std::string& TDhitCollectionName, const std::string& SIDhitCollectionName)
34 : G4VSensitiveDetector( name )
35 , m_nHitID(-1)
36 , m_TDHitCollectionName(TDhitCollectionName)
37 , m_SIDHitCollectionName(SIDhitCollectionName)
38{
39 m_delta_pixel_x = AFP_CONSTANTS::SiT_Pixel_length_x;
40 m_delta_pixel_y = AFP_CONSTANTS::SiT_Pixel_length_y;
41
42 for( int i=0; i < 4; i++){
43 for( int j=0; j < 10; j++){
44 m_death_edge[i][j] = AFP_CONSTANTS::SiT_DeathEdge; //in mm, it is left edge as the movement is horizontal
45 m_lower_edge[i][j] = AFP_CONSTANTS::SiT_LowerEdge; //in mm,
46 }
47 }
48}
49
50//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
51// Initialize from G4.
57
58//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
59
60bool AFP_SensitiveDetector::ProcessHits(G4Step* pStep, G4TouchableHistory*)
61{
62 if (!m_pTDSimHitCollection) {
63 m_pTDSimHitCollection = getTDHitCollection();
64 }
65 if (!m_pSIDSimHitCollection) {
66 m_pSIDSimHitCollection = getSIDHitCollection();
67 }
68 if (!m_pTDSimHitCollection || !m_pSIDSimHitCollection) {
69 return false;
70 }
71
72 if(verboseLevel>5)
73 {
74 G4cout << "AFP_SensitiveDetector::ProcessHits" << G4endl;
75 }
76
77 bool bRes=false;
78 int nTrackID=-1;
79 int nParticleEncoding=-1;
80 float fKineticEnergy=0.0;
81 float fEnergyDeposit=0.0;
82 float fWaveLength=0.0;
83 float fPreStepX=0.0;
84 float fPreStepY=0.0;
85 float fPreStepZ=0.0;
86 float fPostStepX=0.0;
87 float fPostStepY=0.0;
88 float fPostStepZ=0.0;
89 float fGlobalTime=0.0;
90 int nStationID=-1;
91 int nDetectorID=-1;
92 int nQuarticID=-1;
93 // int nPixelRow=-1;
94 // int nPixelCol=-1;
95
96 bool bIsSIDAuxVSID=false;
97
98 // step, track and particle info
99 G4Track* pTrack = pStep->GetTrack();
100 G4ParticleDefinition* pParticleDefinition = pTrack->GetDefinition();
101 G4StepPoint* pPreStepPoint = pStep->GetPreStepPoint();
102 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
103 G4ThreeVector PreStepPointPos = pPreStepPoint->GetPosition();
104 G4ThreeVector PostStepPointPos = pPostStepPoint->GetPosition();
105
106 nTrackID=pTrack->GetTrackID();
107 fKineticEnergy = pPreStepPoint->GetKineticEnergy();
108 fEnergyDeposit = pStep->GetTotalEnergyDeposit();
109
110 fPreStepX = PreStepPointPos.x();
111 fPreStepY = PreStepPointPos.y();
112 fPreStepZ = PreStepPointPos.z();
113 fPostStepX = PostStepPointPos.x();
114 fPostStepY = PostStepPointPos.y();
115 fPostStepZ = PostStepPointPos.z();
116 nParticleEncoding = pParticleDefinition->GetPDGEncoding();
117 fGlobalTime = pStep->GetPreStepPoint()->GetGlobalTime()/CLHEP::picosecond; // time w.r.t. prestep or poststep ??
118
119 // name of physical volume
120 G4TouchableHandle touch1 = pPreStepPoint->GetTouchableHandle();
121 G4VPhysicalVolume* volume = touch1->GetVolume();
122 G4String VolumeName = volume->GetName();
123
124 if(verboseLevel>5)
125 {
126 G4cout << "hit volume name is " << VolumeName << G4endl;
127 G4cout << "particle code is " << nParticleEncoding << "kinetic energy " << fKineticEnergy << G4endl;
128 G4cout << "global, x_pre: " << fPreStepX << ", y_pre: " << fPreStepY << ", z_pre: " << fPreStepZ << G4endl;
129 G4cout << "global, x_post: " << fPostStepX << ", y_post: " << fPostStepY << ", z_post: " << fPostStepZ << G4endl;
130 }
131 //scan station and detector id
132 char* ppv1, *ppv2;
133 char szbuff[32];
134 memset(&szbuff[0],0,sizeof(szbuff));
135 strncpy(szbuff,VolumeName.data(),sizeof(szbuff));
136 szbuff[sizeof(szbuff)-1] = '\0'; // idiomatic use of strncpy...
137 ppv1=strchr(szbuff,'[');
138 ppv2=strchr(szbuff,']');
139 if(!ppv2 || !ppv1){
140 G4cout << "ERROR: Invalid format of volume name " << VolumeName << G4endl;
141 return false;
142 }
143 else *ppv2='\0';
144
145 nStationID=10*(szbuff[3]-0x30)+(szbuff[4]-0x30);
146 nDetectorID=atoi(ppv1+1);
147
148 m_nHitID++;
149
150#if G4VERSION_NUMBER < 1100
151 if ( (VolumeName.contains("TDSensor")) )
152#else
153 if ( (G4StrUtil::contains(VolumeName, "TDSensor")) )
154#endif
155 {
156 nQuarticID=szbuff[7]-0x30;
157 if ( pStep->GetPostStepPoint()->GetProcessDefinedStep() == 0 ) return 1;
158 if ( (pStep->GetPostStepPoint()->GetProcessDefinedStep())->GetProcessName()!="OpAbsorption" ) return 1;
159 fWaveLength = 2.*M_PI*CLHEP::hbarc/(CLHEP::MeV*CLHEP::nm)/fKineticEnergy;
160 if (fWaveLength > 800. || fWaveLength < 200.) return 1; // 200-800 nm cut
161
162 m_pTDSimHitCollection->Emplace(m_nHitID,nTrackID,nParticleEncoding,fKineticEnergy,fEnergyDeposit,
163 fWaveLength,fPreStepX,fPreStepY,fPreStepZ,fPostStepX,fPostStepY,
164 fPostStepZ,fGlobalTime,nStationID,nDetectorID,(1+2*nQuarticID));//Q1: 1-2, Q2: 3-4
165 m_pTDSimHitCollection->CountHit();
166 }
167
169#if G4VERSION_NUMBER < 1100
170 if ( (bRes=VolumeName.contains("TDQuarticBar[")) )
171#else
172 if ( (bRes=G4StrUtil::contains(VolumeName, "TDQuarticBar[")) )
173#endif
174 {
175 nQuarticID=szbuff[7]-0x30;
176
177 // Cut on maximum number of generated photons/bar
178 if (m_pTDSimHitCollection->HasReachedLimit(nStationID, nQuarticID, nDetectorID)) return 1;
179
180 // Get the Touchable History:
181 const G4TouchableHistory* myTouch = static_cast<const G4TouchableHistory*>(pPreStepPoint->GetTouchable());
182 // Calculate the local step position.
183 // From a G4 FAQ:
184 // http://geant4-hn.slac.stanford.edu:5090/HyperNews/public/get/geometry/17/1.html
185 const G4AffineTransform transformation = myTouch->GetHistory()->GetTopTransform();
186 G4ThreeVector PreStepPointPos2 = transformation.TransformPoint(PreStepPointPos);
187 G4ThreeVector PostStepPointPos2 = transformation.TransformPoint(PostStepPointPos);
188
189 G4String shape( myTouch->GetSolid()->GetEntityType() );
190
191 G4ThreeVector normpX( 1., 0., 0.);
192 G4ThreeVector normnX(-1., 0., 0.);
193 G4ThreeVector normpY( 0., 1., 0.);
194 G4ThreeVector normnY( 0.,-1., 0.);
195 G4ThreeVector normpZ( 0., 0., 1.);
196 G4ThreeVector normnZ( 0., 0.,-1.);
197
198 G4double PreProtonX = PreStepPointPos2.x();
199 G4double PreProtonY = PreStepPointPos2.y();
200 G4double PreProtonZ = PreStepPointPos2.z();
201
202 G4double PostProtonX = PostStepPointPos2.x();
203 G4double PostProtonY = PostStepPointPos2.y();
204 G4double PostProtonZ = PostStepPointPos2.z();
205
206 G4ThreeVector p0 = pStep->GetDeltaPosition().unit();
207
208 G4Material* mat = pStep->GetTrack()->GetMaterial();
209 G4MaterialPropertiesTable *matPropTable = mat->GetMaterialPropertiesTable();
210 // Refractive index
211 G4MaterialPropertyVector* Rind = matPropTable->GetProperty("RINDEX");
212 // Absorbtion length
213 G4MaterialPropertyVector* Alen = matPropTable->GetProperty("ABSLENGTH");
214
215 const G4double charge = pParticleDefinition->GetPDGCharge();
216 const G4double beta = (pPreStepPoint->GetBeta() + pPostStepPoint->GetBeta()) / 2.;
217 G4double BetaInverse = 1. / beta;
218
219 //G4int Rsize = Rind->Entries() - 1;
220 //G4double Pmin = Rind->GetMinPhotonEnergy(); // 800 nm
221#if G4VERSION_NUMBER < 1100
222 G4double Pmin = Rind->GetMinLowEdgeEnergy();
223 //G4double Pmax = Rind->GetMaxPhotonEnergy(); // 200 nm
224 G4double Pmax = Rind->GetMaxLowEdgeEnergy();
225#else
226 G4double Pmin = Rind->GetMinEnergy();
227 G4double Pmax = Rind->GetMaxEnergy();
228#endif
229 G4double dp = Pmax - Pmin;
230 //G4double maxCosTheta = BetaInverse / Rind->GetMinProperty();
231 G4double maxCosTheta = BetaInverse / Rind->GetMinValue();
232 G4double maxSin2Theta = (1.0 - maxCosTheta) * (1.0 + maxCosTheta);
233
234 //G4double meanRI = .5*(Rind->GetMinProperty() + Rind->GetMaxProperty());
235 G4double meanRI = .5*(Rind->GetMinValue() + Rind->GetMaxValue());
236
237 // Formula taken from G4 docu (to be changed since the integral approximation)
238 G4double MeanNumberOfPhotons = 370.*(charge/CLHEP::eplus)*(charge/CLHEP::eplus)* (1.0 - 1.0/(beta * meanRI * beta * meanRI)) / (CLHEP::cm*CLHEP::eV);
239 if (MeanNumberOfPhotons <= 0.0) return 1;
240
241 G4double step_length = pStep->GetStepLength();
242
243 MeanNumberOfPhotons = MeanNumberOfPhotons * step_length * dp;
244 G4int NumPhotons = (G4int) G4Poisson( MeanNumberOfPhotons );
245 if (NumPhotons <= 0) return 1;
246 //ATH_MSG_INFO("number of photons: " << NumPhotons);
247
248 G4int NumPhotonsCuts=0;
249 for (G4int I = 0; I < NumPhotons; I++) {
250
251 G4double rand;
252 G4double sampledEnergy, sampledRI;
253 G4double cosTheta, sin2Theta;
254
255 // Sample an energy for photon (using MC elimination method)
256 do {
257 rand = G4UniformRand();
258 sampledEnergy = Pmin + rand * dp;
259 //sampledRI = Rind->GetProperty(sampledEnergy);
260 sampledRI = Rind->Value(sampledEnergy);
261 cosTheta = BetaInverse / sampledRI;
262
263 sin2Theta = (1.0 - cosTheta)*(1.0 + cosTheta);
264 rand = G4UniformRand();
265
266 } while (rand * maxSin2Theta > sin2Theta);
267
268 // Generate random position of photon on the cone surface defined by Theta
269 rand = G4UniformRand();
270 G4double phi = 2.*M_PI*rand;
271 G4double sinPhi = sin(phi);
272 G4double cosPhi = cos(phi);
273
274 // Calculate x,y,z coordinates of photon momentum
275 // in coordinate system with primary particle direction aligned with the z-axis
276 // + Rotate momentum direction back to the global coordinate system
277 G4double sinTheta = std::sqrt(sin2Theta);
278 G4double px = sinTheta*cosPhi;
279 G4double py = sinTheta*sinPhi;
280 G4double pz = cosTheta;
281 G4ParticleMomentum photonMomentum(px, py, pz);
282 photonMomentum.rotateUz(p0);
283
284 G4double PX = photonMomentum.getX();
285 G4double PY = photonMomentum.getY();
286 G4double PZ = photonMomentum.getZ();
287
288 // calculate projections coordinates
289 //G4double PXp = PX/std::sqrt(PX*PX+PY*PY+PZ*PZ);
290 G4double PYp = PY/std::sqrt(PX*PX+PY*PY+PZ*PZ);
291 G4double PZp = PZ/std::sqrt(PX*PX+PY*PY+PZ*PZ);
292
293 G4double PYt = PY/std::sqrt(PY*PY+PZ*PZ);
294 G4double PZt = PZ/std::sqrt(PY*PY+PZ*PZ);
295
296 // Cosines (alpha, delta)
297 cosPhi = (PYp*PYt + PZp*PZt);
298 if (nStationID == 0) cosTheta = ( -PYt*sin(48.*CLHEP::deg) + PZt*cos(48.*CLHEP::deg) );
299 else cosTheta = ( -PYt*sin(48.*CLHEP::deg) - PZt*cos(48.*CLHEP::deg) );
300
301 // Total internal reflection conditions
302 G4double cosThetaC = std::sqrt(1.-1./sampledRI/sampledRI);
303 if (std::sqrt(1.-cosPhi*cosPhi)>cosThetaC) continue;
304 if (std::sqrt(1.-cosTheta*cosTheta)*cosPhi>cosThetaC) continue;
305
306 // Parametric equation of line where photons are generated
307 rand = G4UniformRand();
308 G4double PhotonX = PreProtonX + (PostProtonX-PreProtonX)*rand;
309 G4double PhotonY = PreProtonY + (PostProtonY-PreProtonY)*rand;
310 G4double PhotonZ = PreProtonZ + (PostProtonZ-PreProtonZ)*rand;
311 G4ThreeVector PhotonPos(PhotonX,PhotonY,PhotonZ);
312 G4double PhotonR = std::sqrt( (PreProtonX-PhotonX)*(PreProtonX-PhotonX) + (PreProtonY-PhotonY)*(PreProtonY-PhotonY) + (PreProtonZ-PhotonZ)*(PreProtonZ-PhotonZ) );
313
314 G4double Y0;
315 // including the scattering from the top edge of the bar (perfect absorber now)
316 if (cosTheta>=0)
317 {
318 Y0 = static_cast<G4ReflectedSolid *>(myTouch->GetSolid())->DistanceToOut(PhotonPos, normnY);
319 Y0 = Y0/cosTheta/cosPhi;
320 }
321 else
322 {
323 continue;
324 }
325
326 // absorption of photons inside the crystal
327 //float Pabs = 1. - exp( - Y0/Alen->GetProperty(sampledEnergy) );
328 float Pabs = 1. - exp( - Y0/Alen->Value(sampledEnergy) );
329 rand = G4UniformRand();
330 if (Pabs>rand) continue;
331
332 // maximum PMT efficiency cut (15%) to to avoid crashes due to the too large memory consumption
333 rand = G4UniformRand();
334 if (rand>TDMaxQEff) continue;
336
337 NumPhotonsCuts++;
338
339 float fGlobalTime2 = fGlobalTime;
340 fGlobalTime2 += ( PhotonR * BetaInverse / CLHEP::c_light )/CLHEP::picosecond;
341
342 // for group velocity of light: Edn/dE
343 float EdndE;
344
345 if (sampledEnergy > (Pmin+.5*dp)) EdndE = (sampledRI - Rind->Value(sampledEnergy-0.0001*CLHEP::eV))/0.0001*sampledEnergy/CLHEP::eV;
346 else EdndE = (Rind->Value(sampledEnergy+0.0001*CLHEP::eV) - sampledRI)/0.0001*sampledEnergy/CLHEP::eV;
347 fGlobalTime2 += ( (sampledRI + EdndE)* Y0 * BetaInverse / CLHEP::c_light )/CLHEP::picosecond;
348
349 if (verboseLevel>5)
350 {
351 G4cout << "FastCher EdndE: " << EdndE << G4endl;
352 }
353 fWaveLength = 2.*M_PI*CLHEP::hbarc/sampledEnergy/(CLHEP::MeV*CLHEP::nm);
354
355 // Cut on maximum number of generated photons/bar
356 if (m_pTDSimHitCollection->HasReachedLimit(nStationID, nQuarticID, nDetectorID)) return 1;
357
358 int nSensitiveElementID=-1;
359 if(nQuarticID==0) { nSensitiveElementID=1; }
360 else if(nQuarticID==1) { nSensitiveElementID=3; }
361
362 m_pTDSimHitCollection->Emplace(m_nHitID,nTrackID,nParticleEncoding,fKineticEnergy,fEnergyDeposit,
363 fWaveLength,PhotonX,PhotonY,PhotonZ,(PhotonX+PX),(PhotonY+PY),(PhotonZ+PZ),
364 fGlobalTime2,nStationID,nDetectorID,nSensitiveElementID);
365 m_pTDSimHitCollection->CountHit(nStationID, nQuarticID, nDetectorID);
366 }
367 if(verboseLevel>5)
368 {
369 G4cout << "FastCher number of photons: " << NumPhotonsCuts << G4endl;
370 }
371 }
373#if G4VERSION_NUMBER < 1100
374 else if (VolumeName.contains("SIDSensor") || (bIsSIDAuxVSID=VolumeName.contains("SIDVacuumSensor"))){
375#else
376 else if (G4StrUtil::contains(VolumeName, "SIDSensor") || (bIsSIDAuxVSID=G4StrUtil::contains(VolumeName, "SIDVacuumSensor"))){
377#endif
378 if(!bIsSIDAuxVSID && !(fEnergyDeposit>0.0))
379 {
380 //hit in SID sensor but with zero deposited energy (transportation)
381 }
382 else
383 {
384 if (bIsSIDAuxVSID)
385 {
386 m_pSIDSimHitCollection->Emplace(m_nHitID,nTrackID,nParticleEncoding,fKineticEnergy,fEnergyDeposit,
387 fPreStepX,fPreStepY,fPreStepZ,fPostStepX,fPostStepY,fPostStepZ,
388 fGlobalTime,nStationID,nDetectorID,bIsSIDAuxVSID,-1,-1);
389 }
390 else
391 {
392 // Cut on maximum number of SID hits/station
393 if(m_pSIDSimHitCollection->HasReachedLimit(nStationID)) return 1;
394
395 // Get the Touchable History:
396 const G4TouchableHistory* myTouch = static_cast<const G4TouchableHistory*>(pPreStepPoint->GetTouchable());
397 // Calculate the local step position.
398 // From a G4 FAQ:
399 // http://geant4-hn.slac.stanford.edu:5090/HyperNews/public/get/geometry/17/1.html
400
401 const G4AffineTransform transformation = myTouch->GetHistory()->GetTopTransform();
402 G4ThreeVector localPosition_pre = transformation.TransformPoint(PreStepPointPos);
403 G4ThreeVector localPosition_post = transformation.TransformPoint(PostStepPointPos);
404
405 G4ThreeVector normpX( 1., 0., 0.);
406 G4ThreeVector normnX(-1., 0., 0.);
407 G4ThreeVector normpY( 0., 1., 0.);
408 G4ThreeVector normnY( 0.,-1., 0.);
409 G4ThreeVector normpZ( 0., 0., 1.);
410 G4ThreeVector normnZ( 0., 0.,-1.);
411
412 G4double BarpX = static_cast<G4ReflectedSolid *>(myTouch->GetSolid())->DistanceToOut(localPosition_pre, normpX);
413 G4double BarnX = static_cast<G4ReflectedSolid *>(myTouch->GetSolid())->DistanceToOut(localPosition_pre, normnX);
414 G4double BarpY = static_cast<G4ReflectedSolid *>(myTouch->GetSolid())->DistanceToOut(localPosition_pre, normpY);
415 G4double BarnY = static_cast<G4ReflectedSolid *>(myTouch->GetSolid())->DistanceToOut(localPosition_pre, normnY);
416 G4double BarpZ = static_cast<G4ReflectedSolid *>(myTouch->GetSolid())->DistanceToOut(localPosition_pre, normpZ);
417 G4double BarnZ = static_cast<G4ReflectedSolid *>(myTouch->GetSolid())->DistanceToOut(localPosition_pre, normnZ);
418
419 G4double BarHalfX = .5 * (BarpX+BarnX);
420 G4double BarHalfY = .5 * (BarpY+BarnY);
421 G4double BarHalfZ = .5 * (BarpZ+BarnZ);
422
423 // x should be 50 mu, y should be 250 mu, z should be 250mu
424
425 G4double x_det = BarHalfX + localPosition_pre.x(); // we will remove death edge later - death_edge[nStationID][nDetectorID];
426 G4double y_det = BarHalfY + localPosition_pre.y(); // we will remove lower edge later - lower_edge[nStationID][nDetectorID];
427 G4double z_det = BarHalfZ + localPosition_pre.z();
428
429 G4double x_det_post = BarHalfX + localPosition_post.x(); //- death_edge[nStationID][nDetectorID];
430 G4double y_det_post = BarHalfY + localPosition_post.y(); //- lower_edge[nStationID][nDetectorID];
431 G4double z_det_post = BarHalfZ + localPosition_post.z();
432
433 G4double track_length_XY = std::sqrt(pow(x_det_post-x_det,2)+pow(y_det_post-y_det,2));
434
435 G4double angle_phi_global = atan2(fPostStepY-fPreStepY,fPostStepX-fPreStepX);
436 G4double angle_phi = atan2(y_det_post-y_det,x_det_post-x_det);
437
438 G4double tan_phi = (y_det_post-y_det)/(x_det_post-x_det);
439
440 if(verboseLevel>5)
441 {
442 G4cout << "AFP_SensitiveDetector::ProcessHits: local, x_det: " << x_det << ", y_det: " << y_det << ", z_det: " << z_det << G4endl;
443 G4cout << "AFP_SensitiveDetector::ProcessHits: local, x_det_post: " << x_det_post << ", y_det_post: " << y_det_post << ", z_det_post: " << z_det_post << G4endl;
444 G4cout << "AFP_SensitiveDetector::ProcessHits: angle_phi_global in -pi:pi = " << angle_phi_global << G4endl;
445 }
446 if (angle_phi_global < 0.) angle_phi_global = 2.*M_PI + angle_phi_global;
447 if(verboseLevel>5)
448 {
449 G4cout << "AFP_SensitiveDetector::ProcessHits: angle_phi_global in 0:2pi = " << angle_phi_global << G4endl;
450 G4cout << "AFP_SensitiveDetector::ProcessHits: angle_phi in -pi:pi = " << angle_phi << G4endl;
451 }
452 if (angle_phi < 0.) angle_phi = 2.*M_PI + angle_phi;
453 if(verboseLevel>5)
454 {
455 G4cout << "AFP_SensitiveDetector::ProcessHits: angle_phi in 0:2pi = " << angle_phi << G4endl;
456 }
457 signed int pre_pixel_x = (signed int)(x_det / m_delta_pixel_x);
458 signed int pre_pixel_y = (signed int)(y_det / m_delta_pixel_y);
459 signed int post_pixel_x = (signed int)(x_det_post / m_delta_pixel_x);
460 signed int post_pixel_y = (signed int)(y_det_post / m_delta_pixel_y);
461
462 signed int sign_pixels_x = 0;
463 signed int sign_pixels_y = 0;
464
465 int number_pixels_x = (int) (abs((post_pixel_x-pre_pixel_x)*1.0));
466 int number_pixels_y = (int) (abs((post_pixel_y-pre_pixel_y)*1.0));
467
468 if (number_pixels_x > 0)
469 {
470 sign_pixels_x = (post_pixel_x-pre_pixel_x)/number_pixels_x;
471 }
472 if (number_pixels_y > 0)
473 {
474 sign_pixels_y = (post_pixel_y-pre_pixel_y)/number_pixels_y;
475 }
476
477 int n_death_pixels = (int)(m_death_edge[nStationID][nDetectorID]/m_delta_pixel_x);
478 int n_lower_pixels = (int)(m_lower_edge[nStationID][nDetectorID]/m_delta_pixel_y);
479
480 if(verboseLevel>5)
481 {
482 G4cout << "AFP_SensitiveDetector::ProcessHits: pre: pixel["<< pre_pixel_x - n_death_pixels <<"]["<< pre_pixel_y - n_lower_pixels <<"] was hit" << G4endl;
483 G4cout << "AFP_SensitiveDetector::ProcessHits: post: pixel["<< post_pixel_x - n_death_pixels<<"]["<< post_pixel_y - n_lower_pixels<<"] was hit" << G4endl;
484 G4cout << "AFP_SensitiveDetector::ProcessHits: chip's length in x: " << 2.*BarHalfX << ", in y: " << 2.*BarHalfY << ", in z: " << 2.*BarHalfZ << G4endl;
485 }
486 signed int first = -1;
487
488 G4double x_next_pixel = -9999.;
489 G4double y_next_pixel = -9999.;
490
491 G4double x_border = -9999.;
492 G4double y_border = -9999.;
493
494 G4double pixel_track_length_XY = -1.;
495 G4double angle_2pixel = 10.;
496
497 //number of pixel in death and lower edge
498
499 int act_pixel_x = pre_pixel_x;
500 int act_pixel_y = pre_pixel_y;
501
502 if(verboseLevel>5)
503 {
504 G4cout << "AFP_SensitiveDetector::ProcessHits: actual pixel in x = " << act_pixel_x << ", in y = " << act_pixel_y << G4endl;
505 G4cout << "AFP_SensitiveDetector::ProcessHits: actual compensated pixel in x = " << act_pixel_x - n_death_pixels << ", in y = " << act_pixel_y - n_lower_pixels << G4endl;
506 }
507 if ((number_pixels_x == 0) && (number_pixels_y == 0))
508 {
509
510 if(verboseLevel>5)
511 {
512 G4cout << "AFP_SensitiveDetector::ProcessHits: pre and post in the same pixel " << G4endl;
513 }
514 if (( pre_pixel_y - n_lower_pixels <= 80) && (pre_pixel_x -n_death_pixels <= 336) && ( pre_pixel_y - n_lower_pixels > 0) && (pre_pixel_x - n_death_pixels > 0))
515 {
516 m_pSIDSimHitCollection->Emplace(m_nHitID,nTrackID,nParticleEncoding,fKineticEnergy,fEnergyDeposit,
517 fPreStepX,fPreStepY,fPreStepZ,fPostStepX,fPostStepY,fPostStepZ,
518 fGlobalTime,nStationID,nDetectorID,bIsSIDAuxVSID,
519 (pre_pixel_y - n_lower_pixels - 1),
520 (pre_pixel_x - n_death_pixels - 1));
521 m_pSIDSimHitCollection->CountHit(nStationID);
522 }
523 else if(verboseLevel>5)
524 {
525 G4cout << "AFP_SensitiveDetector::ProcessHits: hit outside of pixel's sensitive area " << G4endl;
526 }
527 }
528 else
529 {
530 if(verboseLevel>5)
531 {
532 G4cout << "AFP_SensitiveDetector::ProcessHits: pre and post in diferent pixels " << G4endl;
533 }
534 // still not complete logic, last step must be cut
535
536 while ( (number_pixels_x >= 0) && (number_pixels_y >= 0) )
537 {
538
539 if ((angle_phi >= 0.) && (angle_phi < M_PI_2))
540 {
541 x_next_pixel = (act_pixel_x+1)*m_delta_pixel_x; //- death_edge[nStationID][nDetectorID];
542 y_next_pixel = (act_pixel_y+1)*m_delta_pixel_y; //- death_edge[nStationID][nDetectorID];
543 angle_2pixel = atan2(y_next_pixel-y_det,x_next_pixel-x_det);
544
545 if (angle_2pixel < 0.) angle_2pixel = 2*M_PI + angle_2pixel;
546 if(verboseLevel>5) { G4cout << "AFP_SensitiveDetector::ProcessHits: angle_2pixel in 0:2pi = " << angle_2pixel << G4endl; }
547
548 if (angle_2pixel > angle_phi)
549 {
550 first = 0; // x border will be hit first
551 }
552 else
553 {
554 first = 1;
555 }
556 }
557
558 else if ((angle_phi >= M_PI_2) && (angle_phi < M_PI))
559 {
560 x_next_pixel = (act_pixel_x+0)*m_delta_pixel_x; //- death_edge[nStationID][nDetectorID];
561 y_next_pixel = (act_pixel_y+1)*m_delta_pixel_y; //- death_edge[nStationID][nDetectorID];
562 angle_2pixel = atan2(y_next_pixel-y_det,x_next_pixel-x_det);
563
564 if (angle_2pixel < 0.) angle_2pixel = 2*M_PI + angle_2pixel;
565 if(verboseLevel>5) { G4cout << "AFP_SensitiveDetector::ProcessHits: angle_2pixel in 0:2pi = " << angle_2pixel << G4endl; }
566
567 if (angle_2pixel > angle_phi)
568 {
569 first = 1; // y border will be hit first
570 }
571 else
572 {
573 first = 0;
574 }
575 }
576
577 else if ((angle_phi >= M_PI) && (angle_phi < 3.*M_PI_2))
578 {
579 x_next_pixel = (act_pixel_x+0)*m_delta_pixel_x; //- death_edge[nStationID][nDetectorID];
580 y_next_pixel = (act_pixel_y+0)*m_delta_pixel_y; //- death_edge[nStationID][nDetectorID];
581 if(verboseLevel>5) { G4cout << "AFP_SensitiveDetector::ProcessHits: next pixel corner, x = " << x_next_pixel << ", y =" << y_next_pixel << G4endl; }
582
583 angle_2pixel = atan2(y_next_pixel-y_det,x_next_pixel-x_det);
584
585 if (angle_2pixel < 0.) angle_2pixel = 2*M_PI + angle_2pixel;
586 if(verboseLevel>5) { G4cout << "AFP_SensitiveDetector::ProcessHits: angle_2pixel in 0:2pi = " << angle_2pixel << G4endl; }
587
588 if (angle_2pixel > angle_phi)
589 {
590 first = 0; // x border will be hit first
591 }
592 else
593 {
594 first = 1;
595 }
596 }
597
598 else if ((angle_phi >= 3.*M_PI_2) && (angle_phi < 2.*M_PI))
599 {
600 x_next_pixel = (act_pixel_x+1)*m_delta_pixel_x; //- death_edge[nStationID][nDetectorID];
601 y_next_pixel = (act_pixel_y+0)*m_delta_pixel_y; //- death_edge[nStationID][nDetectorID];
602 angle_2pixel = atan2(y_next_pixel-y_det,x_next_pixel-x_det);
603
604 if (angle_2pixel < 0.) angle_2pixel = 2*M_PI + angle_2pixel;
605 if(verboseLevel>5) { G4cout << "AFP_SensitiveDetector::ProcessHits: angle_2pixel in 0:2pi = " << angle_2pixel << G4endl; }
606
607 if (angle_2pixel > angle_phi)
608 {
609 first = 1; // y border will be hit first
610 }
611 else
612 {
613 first = 0;
614 }
615 }
616
617 else
618 {
619 if(verboseLevel>5) { G4cout << "AFP_SensitiveDetector::ProcessHits: something is wrong here!!! " << G4endl; }
620 }
621
622
623 if (first == -1 ) {
624 if(verboseLevel>5) { G4cout << "AFP_SensitiveDetector::ProcessHits: something is wrong here!!! " << G4endl; }
625 throw std::runtime_error("AFP_SensitiveDetector::ProcessHits: something is wrong here");
626 }
627
628 if(verboseLevel>5)
629 {
630 G4cout << "AFP_SensitiveDetector::ProcessHits: actual pixel in x = " << act_pixel_x << ", in y = " << act_pixel_y << G4endl;
631 G4cout << "AFP_SensitiveDetector::ProcessHits: actual compensated pixel in x = " << act_pixel_x - n_death_pixels << ", in y = " << act_pixel_y - n_lower_pixels << G4endl;
632 }
633
634 if (first == 0 )
635 {
636
637 if(verboseLevel>5) { G4cout << "AFP_SensitiveDetector::ProcessHits: cross is x, " << G4endl; }
638 x_border = x_next_pixel;
639
640 if ((sign_pixels_x >= 0) && (x_border > x_det_post)) x_border = x_det_post;
641 if ((sign_pixels_x < 0) && (x_border < x_det_post)) x_border = x_det_post;
642
643 y_border = tan_phi*(x_border-x_det) + y_det;
644
645 if (( act_pixel_y - n_lower_pixels <= 80) && (act_pixel_x -n_death_pixels <= 336) && ( act_pixel_y - n_lower_pixels > 0) && (act_pixel_x - n_death_pixels > 0))
646 {
647 pixel_track_length_XY = std::sqrt(pow(x_border-x_det,2)+pow(y_border-y_det,2));
648
649 if(verboseLevel>5)
650 {
651 G4cout << "AFP_SensitiveDetector::ProcessHits: overall energy = " << fEnergyDeposit << G4endl;
652 G4cout << "AFP_SensitiveDetector::ProcessHits: track XY length = " << track_length_XY << G4endl;
653 G4cout << "AFP_SensitiveDetector::ProcessHits: actual XY length = " << pixel_track_length_XY << G4endl;
654 G4cout << "AFP_SensitiveDetector::ProcessHits: deposited energy = " << fEnergyDeposit*(pixel_track_length_XY/track_length_XY) << G4endl;
655 }
656
657 // this logic has to be still checked, tracks is necessary fully in sensitive area, but logic is probably ok
658 m_pSIDSimHitCollection->Emplace(m_nHitID,nTrackID,nParticleEncoding,fKineticEnergy,
659 fEnergyDeposit*(pixel_track_length_XY/track_length_XY),
660 fPreStepX,fPreStepY,fPreStepZ,fPostStepX,fPostStepY,fPostStepZ,
661 fGlobalTime,nStationID,nDetectorID,bIsSIDAuxVSID,
662 (act_pixel_y - n_lower_pixels - 1),
663 (act_pixel_x - n_death_pixels - 1));
664
665 if(verboseLevel>5) { G4cout << "AFP_SensitiveDetector::ProcessHits:pixel["<< act_pixel_x - n_death_pixels <<"]["<< act_pixel_y - n_lower_pixels <<"] will be stored, with energy "
666 << fEnergyDeposit*(pixel_track_length_XY/track_length_XY) << G4endl; }
667
668 m_pSIDSimHitCollection->CountHit(nStationID);
669 }
670
671 x_det = x_border;
672 y_det = y_border;
673
674 x_next_pixel = x_next_pixel + sign_pixels_x*m_delta_pixel_x;
675 number_pixels_x = number_pixels_x - 1;
676
677 if(verboseLevel>5) { G4cout << "AFP_SensitiveDetector::ProcessHits: remaining number of pixels in x = " << number_pixels_x << ", in y = " << number_pixels_y << G4endl; }
678
679 act_pixel_x = act_pixel_x + sign_pixels_x;
680 }
681
682 if (first == 1 )
683 {
684
685 if(verboseLevel>5) { G4cout << "AFP_SensitiveDetector::ProcessHits: cross is y, " << G4endl; }
686 y_border = y_next_pixel;
687
688 if ((sign_pixels_y >= 0) && (y_border > y_det_post)) y_border = y_det_post;
689 if ((sign_pixels_y < 0) && (y_border < y_det_post)) y_border = y_det_post;
690
691 x_border = (y_border-y_det)/tan_phi + x_det;
692
693 if (( act_pixel_y - n_lower_pixels <= 80) && (act_pixel_x -n_death_pixels <= 336) && ( act_pixel_y - n_lower_pixels > 0) && (act_pixel_x - n_death_pixels > 0))
694 {
695 pixel_track_length_XY = std::sqrt(pow(x_border-x_det,2)+pow(y_border-y_det,2));
696
697 if(verboseLevel>5)
698 {
699 G4cout << "AFP_SensitiveDetector::ProcessHits: overall energy = " << fEnergyDeposit << G4endl;
700 G4cout << "AFP_SensitiveDetector::ProcessHits: track XY length = " << track_length_XY << G4endl;
701 G4cout << "AFP_SensitiveDetector::ProcessHits: actual XY length = " << pixel_track_length_XY << G4endl;
702 G4cout << "AFP_SensitiveDetector::ProcessHits: deposited energy = " << fEnergyDeposit*(pixel_track_length_XY/track_length_XY) << G4endl;
703 }
704
705 // this logic has to be still checked, tracks is necessary fully in sensitive area, but logic is probably ok
706 m_pSIDSimHitCollection->Emplace(m_nHitID,nTrackID,nParticleEncoding,fKineticEnergy,
707 fEnergyDeposit*(pixel_track_length_XY/track_length_XY),
708 fPreStepX,fPreStepY,fPreStepZ,fPostStepX,fPostStepY,fPostStepZ,
709 fGlobalTime,nStationID,nDetectorID,bIsSIDAuxVSID,
710 (act_pixel_y - n_lower_pixels - 1),
711 (act_pixel_x - n_death_pixels - 1));
712
713 if(verboseLevel>5) { G4cout << "AFP_SensitiveDetector::ProcessHits:pixel["<< act_pixel_x - n_death_pixels <<"]["<< act_pixel_y - n_lower_pixels <<"] will be stored, with energy "
714 << fEnergyDeposit*(pixel_track_length_XY/track_length_XY) << G4endl; }
715
716 m_pSIDSimHitCollection->CountHit(nStationID);
717 }
718
719 y_det = y_border;
720 x_det = x_border;
721
722 y_next_pixel = y_next_pixel + sign_pixels_y*m_delta_pixel_y;
723 number_pixels_y = number_pixels_y - 1;
724
725 if(verboseLevel>5) { G4cout << "AFP_SensitiveDetector::ProcessHits: remaining number of pixels in x = " << number_pixels_x << ", in y = " << number_pixels_y << G4endl; }
726
727
728 act_pixel_y = act_pixel_y + sign_pixels_y;
729 }
730 }
731 }
732 }
733 }
734 }
735
736 return true;
737}
738
739AFP_TDSimHitCollectionBuilder* AFP_SensitiveDetector::getTDHitCollection() const
740{
741 auto* eventInfo = AtlasG4EventUserInfo::GetEventUserInfo();
742 if (!eventInfo) {
743 return nullptr;
744 }
745 auto hitCollections = eventInfo->GetHitCollectionMap();
746 return hitCollections ? hitCollections->Find<AFP_TDSimHitCollectionBuilder>(m_TDHitCollectionName) : nullptr;
747}
748
749AFP_SIDSimHitCollectionBuilder* AFP_SensitiveDetector::getSIDHitCollection() const
750{
751 auto* eventInfo = AtlasG4EventUserInfo::GetEventUserInfo();
752 if (!eventInfo) {
753 return nullptr;
754 }
755 auto hitCollections = eventInfo->GetHitCollectionMap();
756 return hitCollections ? hitCollections->Find<AFP_SIDSimHitCollectionBuilder>(m_SIDHitCollectionName) : nullptr;
757}
AFP_SensitiveDetector.h
AFP_constants.h
M_PI
#define M_PI
Definition ActiveFraction.h:11
phi
Scalar phi() const
phi method
Definition AmgMatrixBasePlugin.h:67
charge
double charge(const T &p)
Definition AtlasPID.h:997
I
#define I(x, y, z)
Definition MD5.cxx:116
AFP_SIDSimHitCollectionBuilder
Definition AFP_HitCollectionBuilders.h:55
AFP_SensitiveDetector::getTDHitCollection
AFP_TDSimHitCollectionBuilder * getTDHitCollection() const
Definition AFP_SensitiveDetector.cxx:739
AFP_SensitiveDetector::m_delta_pixel_x
float m_delta_pixel_x
Definition AFP_SensitiveDetector.h:51
AFP_SensitiveDetector::getSIDHitCollection
AFP_SIDSimHitCollectionBuilder * getSIDHitCollection() const
Definition AFP_SensitiveDetector.cxx:749
AFP_SensitiveDetector::AFP_SensitiveDetector
AFP_SensitiveDetector(const std::string &name, const std::string &TDhitCollectionName, const std::string &SIDhitCollectionName)
Definition AFP_SensitiveDetector.cxx:33
AFP_SensitiveDetector::m_TDHitCollectionName
std::string m_TDHitCollectionName
Definition AFP_SensitiveDetector.h:56
AFP_SensitiveDetector::TDMaxQEff
static constexpr double TDMaxQEff
Templated method to stuff a single hit into the sensitive detector class.
Definition AFP_SensitiveDetector.h:43
AFP_SensitiveDetector::Initialize
void Initialize(G4HCofThisEvent *) override final
Definition AFP_SensitiveDetector.cxx:52
AFP_SensitiveDetector::ProcessHits
G4bool ProcessHits(G4Step *, G4TouchableHistory *) override final
Definition AFP_SensitiveDetector.cxx:60
AFP_SensitiveDetector::m_nHitID
int m_nHitID
Definition AFP_SensitiveDetector.h:49
AFP_SensitiveDetector::m_pSIDSimHitCollection
AFP_SIDSimHitCollectionBuilder * m_pSIDSimHitCollection
Definition AFP_SensitiveDetector.h:60
AFP_SensitiveDetector::m_SIDHitCollectionName
std::string m_SIDHitCollectionName
Definition AFP_SensitiveDetector.h:57
AFP_SensitiveDetector::m_delta_pixel_y
float m_delta_pixel_y
Definition AFP_SensitiveDetector.h:51
AFP_SensitiveDetector::m_pTDSimHitCollection
AFP_TDSimHitCollectionBuilder * m_pTDSimHitCollection
Definition AFP_SensitiveDetector.h:59
AFP_SensitiveDetector::m_death_edge
float m_death_edge[4][10]
Definition AFP_SensitiveDetector.h:52
AFP_SensitiveDetector::m_lower_edge
float m_lower_edge[4][10]
Definition AFP_SensitiveDetector.h:53
AFP_TDSimHitCollectionBuilder
Definition AFP_HitCollectionBuilders.h:14
AtlasG4EventUserInfo::GetEventUserInfo
static AtlasG4EventUserInfo * GetEventUserInfo()
Definition AtlasG4EventUserInfo.h:130
AFP_CONSTANTS::SiT_Pixel_length_x
static constexpr double SiT_Pixel_length_x
Definition AFP_constants.h:58
AFP_CONSTANTS::SiT_LowerEdge
static constexpr double SiT_LowerEdge
Definition AFP_constants.h:68
AFP_CONSTANTS::SiT_Pixel_length_y
static constexpr double SiT_Pixel_length_y
Definition AFP_constants.h:59
AFP_CONSTANTS::SiT_DeathEdge
static constexpr double SiT_DeathEdge
Definition AFP_constants.h:67
Generated on for ATLAS Offline Software by doxygen 1.16.1