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