ATLAS Offline Software
AFP_SensitiveDetector.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 // Class header
7 
8 // Athena headers
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 
31 AFP_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 {
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
51  }
52  }
53 }
54 
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 
58 {
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
78 void 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 
86 bool 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
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 = ((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);
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 = ((G4ReflectedSolid *)(myTouch->GetSolid()))->DistanceToOut(localPosition_pre, normpX);
440  G4double BarnX = ((G4ReflectedSolid *)(myTouch->GetSolid()))->DistanceToOut(localPosition_pre, normnX);
441  G4double BarpY = ((G4ReflectedSolid *)(myTouch->GetSolid()))->DistanceToOut(localPosition_pre, normpY);
442  G4double BarnY = ((G4ReflectedSolid *)(myTouch->GetSolid()))->DistanceToOut(localPosition_pre, normnY);
443  G4double BarpZ = ((G4ReflectedSolid *)(myTouch->GetSolid()))->DistanceToOut(localPosition_pre, normpZ);
444  G4double BarnZ = ((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));
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 
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 
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 
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++;
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 }
TileDCSDataPlotter.dp
dp
Definition: TileDCSDataPlotter.py:840
AFP_SensitiveDetector::EndOfAthenaEvent
void EndOfAthenaEvent()
Definition: AFP_SensitiveDetector.cxx:774
AFP_SensitiveDetector::ProcessHits
G4bool ProcessHits(G4Step *, G4TouchableHistory *) override final
Definition: AFP_SensitiveDetector.cxx:86
AFP_CONSTANTS::SiT_Pixel_length_x
static constexpr double SiT_Pixel_length_x
Definition: AFP_constants.h:58
test_pyathena.px
px
Definition: test_pyathena.py:18
AFP_SensitiveDetector::m_nNOfTDSimHits
int m_nNOfTDSimHits[4][32]
Definition: AFP_SensitiveDetector.h:62
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:67
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
AFP_CONSTANTS::SiT_Pixel_length_y
static constexpr double SiT_Pixel_length_y
Definition: AFP_constants.h:59
mat
GeoMaterial * mat
Definition: LArDetectorConstructionTBEC.cxx:55
AFP_constants.h
python.SystemOfUnits.MeV
int MeV
Definition: SystemOfUnits.py:154
AFP_SensitiveDetector::m_nNumberOfTDSimHits
int m_nNumberOfTDSimHits
Definition: AFP_SensitiveDetector.h:59
M_PI
#define M_PI
Definition: ActiveFraction.h:11
deg
#define deg
Definition: SbPolyhedron.cxx:17
AFP_SensitiveDetector::m_delta_pixel_y
float m_delta_pixel_y
Definition: AFP_SensitiveDetector.h:65
AFP_CONSTANTS::SiT_LowerEdge
static constexpr double SiT_LowerEdge
Definition: AFP_constants.h:68
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
drawFromPickle.exp
exp
Definition: drawFromPickle.py:36
AFP_SensitiveDetector::m_pTDSimHitCollection
SG::WriteHandle< AFP_TDSimHitCollection > m_pTDSimHitCollection
Definition: AFP_SensitiveDetector.h:70
AFP_SensitiveDetector::m_nNOfSIDSimHits
int m_nNOfSIDSimHits[4]
Definition: AFP_SensitiveDetector.h:63
AFP_SensitiveDetector::AFP_SensitiveDetector
AFP_SensitiveDetector(const std::string &name, const std::string &TDhitCollectionName, const std::string &SIDhitCollectionName)
Definition: AFP_SensitiveDetector.cxx:31
cm
const double cm
Definition: Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/FCAL_ChannelMap.cxx:25
AFP_SensitiveDetector::m_nHitID
int m_nHitID
Definition: AFP_SensitiveDetector.h:57
LArG4FSStartPointFilter.rand
rand
Definition: LArG4FSStartPointFilter.py:80
python.PhysicalConstants.hbarc
float hbarc
Definition: PhysicalConstants.py:73
AFP_SensitiveDetector::m_nEventNumber
int m_nEventNumber
Definition: AFP_SensitiveDetector.h:58
lumiFormat.i
int i
Definition: lumiFormat.py:85
AFP_SensitiveDetector::TDMaxCnt
static constexpr int TDMaxCnt
Definition: AFP_SensitiveDetector.h:53
contains
bool contains(const std::string &s, const std::string &regx)
does a string contain the substring
Definition: hcg.cxx:111
Amg::pz
@ pz
Definition: GeoPrimitives.h:40
AFP_SensitiveDetector::m_pSIDSimHitCollection
SG::WriteHandle< AFP_SIDSimHitCollection > m_pSIDSimHitCollection
Definition: AFP_SensitiveDetector.h:71
AFP_SensitiveDetector::Initialize
void Initialize(G4HCofThisEvent *) override final
Definition: AFP_SensitiveDetector.cxx:78
AFP_SensitiveDetector.h
AFP_SensitiveDetector::m_lower_edge
float m_lower_edge[4][10]
Definition: AFP_SensitiveDetector.h:67
Amg::py
@ py
Definition: GeoPrimitives.h:39
AFP_SensitiveDetector::SiDMaxCnt
static constexpr int SiDMaxCnt
Definition: AFP_SensitiveDetector.h:54
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:228
python.SystemOfUnits.eV
int eV
Definition: SystemOfUnits.py:155
AFP_CONSTANTS::SiT_DeathEdge
static constexpr double SiT_DeathEdge
Definition: AFP_constants.h:67
charge
double charge(const T &p)
Definition: AtlasPID.h:756
python.PhysicalConstants.c_light
float c_light
Definition: PhysicalConstants.py:63
AFP_SensitiveDetector::m_delta_pixel_x
float m_delta_pixel_x
Definition: AFP_SensitiveDetector.h:65
CalibCoolCompareRT.nm
nm
Definition: CalibCoolCompareRT.py:110
AFP_SensitiveDetector::m_death_edge
float m_death_edge[4][10]
Definition: AFP_SensitiveDetector.h:66
DeMoScan.first
bool first
Definition: DeMoScan.py:536
AFP_SensitiveDetector::StartOfAthenaEvent
void StartOfAthenaEvent()
Definition: AFP_SensitiveDetector.cxx:57
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:85
python.SystemOfUnits.eplus
int eplus
Definition: SystemOfUnits.py:137
AFP_SensitiveDetector::m_nNumberOfSIDSimHits
int m_nNumberOfSIDSimHits
Definition: AFP_SensitiveDetector.h:60
I
#define I(x, y, z)
Definition: MD5.cxx:116
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
MuonParameters::beta
@ beta
Definition: MuonParamDefs.h:144
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
AFP_SensitiveDetector::TDMaxQEff
static constexpr double TDMaxQEff
Templated method to stuff a single hit into the sensitive detector class.
Definition: AFP_SensitiveDetector.h:52
python.SystemOfUnits.picosecond
int picosecond
Definition: SystemOfUnits.py:123
TRTCalib_cfilter.p0
p0
Definition: TRTCalib_cfilter.py:129