ATLAS Offline Software
Loading...
Searching...
No Matches
Public Member Functions | Static Public Attributes | Private Member Functions | Private Attributes | List of all members
AFP_SensitiveDetector Class Reference

#include <AFP_SensitiveDetector.h>

Inheritance diagram for AFP_SensitiveDetector:
Collaboration diagram for AFP_SensitiveDetector:

Public Member Functions

 AFP_SensitiveDetector (const std::string &name, const std::string &TDhitCollectionName, const std::string &SIDhitCollectionName)
 ~AFP_SensitiveDetector ()
void Initialize (G4HCofThisEvent *) override final
G4bool ProcessHits (G4Step *, G4TouchableHistory *) override final

Static Public Attributes

static constexpr double TDMaxQEff = 0.15
 Templated method to stuff a single hit into the sensitive detector class.

Private Member Functions

 FRIEND_TEST (AFP_SensitiveDetectortest, ProcessHits1)
 FRIEND_TEST (AFP_SensitiveDetectortest, ProcessHits2)
AFP_TDSimHitCollectionBuilder * getTDHitCollection () const
AFP_SIDSimHitCollectionBuilder * getSIDHitCollection () const

Private Attributes

int m_nHitID
float m_delta_pixel_x
float m_delta_pixel_y
float m_death_edge [4][10]
float m_lower_edge [4][10]
std::string m_TDHitCollectionName
std::string m_SIDHitCollectionName
AFP_TDSimHitCollectionBuilder * m_pTDSimHitCollection {}
AFP_SIDSimHitCollectionBuilder * m_pSIDSimHitCollection {}

Detailed Description

Definition at line 22 of file AFP_SensitiveDetector.h.

Constructor & Destructor Documentation

◆ AFP_SensitiveDetector()

AFP_SensitiveDetector::AFP_SensitiveDetector ( const std::string & name,
const std::string & TDhitCollectionName,
const std::string & SIDhitCollectionName )

Definition at line 33 of file AFP_SensitiveDetector.cxx.

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}
AFP_SensitiveDetector::m_delta_pixel_x
float m_delta_pixel_x
Definition AFP_SensitiveDetector.h:51
AFP_SensitiveDetector::m_TDHitCollectionName
std::string m_TDHitCollectionName
Definition AFP_SensitiveDetector.h:56
AFP_SensitiveDetector::m_nHitID
int m_nHitID
Definition AFP_SensitiveDetector.h:49
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_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
TauHitVars::j
float j(const xAOD::IParticle &, const xAOD::TrackMeasurementValidation &hit, const Eigen::Matrix3d &jab_inv)
Definition ConstituentLoaderTauHit.cxx:102
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

◆ ~AFP_SensitiveDetector()

AFP_SensitiveDetector::~AFP_SensitiveDetector ( )
inline

Definition at line 32 of file AFP_SensitiveDetector.h.

32{ /* I don't own myHitColl if all has gone well */ }

Member Function Documentation

◆ FRIEND_TEST() [1/2]

AFP_SensitiveDetector::FRIEND_TEST ( AFP_SensitiveDetectortest ,
ProcessHits1  )
private

◆ FRIEND_TEST() [2/2]

AFP_SensitiveDetector::FRIEND_TEST ( AFP_SensitiveDetectortest ,
ProcessHits2  )
private

◆ getSIDHitCollection()

AFP_SIDSimHitCollectionBuilder * AFP_SensitiveDetector::getSIDHitCollection ( ) const
private

Definition at line 749 of file AFP_SensitiveDetector.cxx.

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}
AtlasG4EventUserInfo::GetEventUserInfo
static AtlasG4EventUserInfo * GetEventUserInfo()
Definition AtlasG4EventUserInfo.h:130

◆ getTDHitCollection()

AFP_TDSimHitCollectionBuilder * AFP_SensitiveDetector::getTDHitCollection ( ) const
private

Definition at line 739 of file AFP_SensitiveDetector.cxx.

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}

◆ Initialize()

void AFP_SensitiveDetector::Initialize ( G4HCofThisEvent * )
finaloverride

Definition at line 52 of file AFP_SensitiveDetector.cxx.

53{
54 m_pTDSimHitCollection = getTDHitCollection();
55 m_pSIDSimHitCollection = getSIDHitCollection();
56}
AFP_SensitiveDetector::getTDHitCollection
AFP_TDSimHitCollectionBuilder * getTDHitCollection() const
Definition AFP_SensitiveDetector.cxx:739
AFP_SensitiveDetector::getSIDHitCollection
AFP_SIDSimHitCollectionBuilder * getSIDHitCollection() const
Definition AFP_SensitiveDetector.cxx:749
AFP_SensitiveDetector::m_pSIDSimHitCollection
AFP_SIDSimHitCollectionBuilder * m_pSIDSimHitCollection
Definition AFP_SensitiveDetector.h:60
AFP_SensitiveDetector::m_pTDSimHitCollection
AFP_TDSimHitCollectionBuilder * m_pTDSimHitCollection
Definition AFP_SensitiveDetector.h:59

◆ ProcessHits()

bool AFP_SensitiveDetector::ProcessHits ( G4Step * pStep,
G4TouchableHistory *  )
finaloverride

Definition at line 60 of file AFP_SensitiveDetector.cxx.

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}
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_SensitiveDetector::TDMaxQEff
static constexpr double TDMaxQEff
Templated method to stuff a single hit into the sensitive detector class.
Definition AFP_SensitiveDetector.h:43
CxxUtils::atoi
int atoi(std::string_view str)
Helper functions to unpack numbers decoded in string into integers and doubles The strings are requir...
Definition Control/CxxUtils/Root/StringUtils.cxx:40
DeMoScan.first
bool first
Definition DeMoScan.py:534
TRTCalib_cfilter.p0
p0
Definition TRTCalib_cfilter.py:129
conifer::pow
constexpr int pow(int x)
Definition conifer.h:27
pdg_comparison.mat
mat
Definition pdg_comparison.py:322
test_pyathena.px
px
Definition test_pyathena.py:18

Member Data Documentation

◆ m_death_edge

float AFP_SensitiveDetector::m_death_edge[4][10]
private

Definition at line 52 of file AFP_SensitiveDetector.h.

◆ m_delta_pixel_x

float AFP_SensitiveDetector::m_delta_pixel_x
private

Definition at line 51 of file AFP_SensitiveDetector.h.

◆ m_delta_pixel_y

float AFP_SensitiveDetector::m_delta_pixel_y
private

Definition at line 51 of file AFP_SensitiveDetector.h.

◆ m_lower_edge

float AFP_SensitiveDetector::m_lower_edge[4][10]
private

Definition at line 53 of file AFP_SensitiveDetector.h.

◆ m_nHitID

int AFP_SensitiveDetector::m_nHitID
private

Definition at line 49 of file AFP_SensitiveDetector.h.

◆ m_pSIDSimHitCollection

AFP_SIDSimHitCollectionBuilder* AFP_SensitiveDetector::m_pSIDSimHitCollection {}
private

Definition at line 60 of file AFP_SensitiveDetector.h.

60{};

◆ m_pTDSimHitCollection

AFP_TDSimHitCollectionBuilder* AFP_SensitiveDetector::m_pTDSimHitCollection {}
private

Definition at line 59 of file AFP_SensitiveDetector.h.

59{};

◆ m_SIDHitCollectionName

std::string AFP_SensitiveDetector::m_SIDHitCollectionName
private

Definition at line 57 of file AFP_SensitiveDetector.h.

◆ m_TDHitCollectionName

std::string AFP_SensitiveDetector::m_TDHitCollectionName
private

Definition at line 56 of file AFP_SensitiveDetector.h.

◆ TDMaxQEff

double AFP_SensitiveDetector::TDMaxQEff = 0.15
staticconstexpr

Templated method to stuff a single hit into the sensitive detector class.

This could get rather tricky, but the idea is to allow fast simulations to use the very same SD classes as the standard simulation.

Definition at line 43 of file AFP_SensitiveDetector.h.

The documentation for this class was generated from the following files:
Generated on for ATLAS Offline Software by doxygen 1.16.1