ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4AdjointCrossSurfChecker.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4AdjointCrossSurfChecker.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
28 // Class Name: G4AdjointCrossSurfChecker
29 // Author: L. Desorgher
30 // Organisation: SpaceIT GmbH
31 // Contract: ESA contract 21435/08/NL/AT
32 // Customer: ESA/ESTEC
34 
35 #include "G4AdjointCrossSurfChecker.hh"
36 #include "G4PhysicalConstants.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4Step.hh"
39 #include "G4StepPoint.hh"
40 #include "G4PhysicalVolumeStore.hh"
41 #include "G4VSolid.hh"
42 #include "G4AffineTransform.hh"
43 
45 //
46 G4ThreadLocal G4AdjointCrossSurfChecker* G4AdjointCrossSurfChecker::instance = 0;
47 
49 //
50 G4AdjointCrossSurfChecker::G4AdjointCrossSurfChecker()
51 {;
52 }
54 //
55 G4AdjointCrossSurfChecker::~G4AdjointCrossSurfChecker()
56 {
57  delete instance;
58 }
60 //
61 G4AdjointCrossSurfChecker* G4AdjointCrossSurfChecker::GetInstance()
62 {
63  if (!instance) instance = new G4AdjointCrossSurfChecker();
64  return instance;
65 }
67 //
68 G4bool G4AdjointCrossSurfChecker::CrossingASphere(const G4Step* aStep,G4double sphere_radius, G4ThreeVector sphere_center,G4ThreeVector& crossing_pos, G4double& cos_th , G4bool& GoingIn)
69 {
70  G4ThreeVector pos1= aStep->GetPreStepPoint()->GetPosition() - sphere_center;
71  G4ThreeVector pos2= aStep->GetPostStepPoint()->GetPosition() - sphere_center;
72  G4double r1= pos1.mag();
73  G4double r2= pos2.mag();
74  G4bool did_cross =false;
75 
76  if (r1<=sphere_radius && r2>sphere_radius){
77  did_cross=true;
78  GoingIn=false;
79  }
80  else if (r2<=sphere_radius && r1>sphere_radius){
81  did_cross=true;
82  GoingIn=true;
83  }
84 
85  if (did_cross) {
86 
87  G4ThreeVector dr=pos2-pos1;
88  G4double r12 = r1*r1;
89  G4double rdr = dr.mag();
90  G4double a,b,c,d;
91  a = rdr*rdr;
92  b = 2.*pos1.dot(dr);
93  c = r12-sphere_radius*sphere_radius;
94  d=std::sqrt(b*b-4.*a*c);
95  G4double l= (-b+d)/2./a;
96  if (l > 1.) l=(-b-d)/2./a;
97  crossing_pos=pos1+l*dr;
98  cos_th = std::abs(dr.cosTheta(crossing_pos));
99 
100 
101  }
102  return did_cross;
103 }
105 //
106 G4bool G4AdjointCrossSurfChecker::GoingInOrOutOfaVolume(const G4Step* aStep,const G4String& volume_name, G4double& , G4bool& GoingIn) //from external surface
107 {
108  G4bool step_at_boundary = (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
109  G4bool did_cross =false;
110  if (step_at_boundary){
111  const G4VTouchable* postStepTouchable = aStep->GetPostStepPoint()->GetTouchable();
112  const G4VTouchable* preStepTouchable = aStep->GetPreStepPoint()->GetTouchable();
113  if (preStepTouchable && postStepTouchable && postStepTouchable->GetVolume() && preStepTouchable->GetVolume()){
114  G4String post_vol_name = postStepTouchable->GetVolume()->GetName();
115  G4String pre_vol_name = preStepTouchable->GetVolume()->GetName();
116 
117  if (post_vol_name == volume_name ){
118  GoingIn=true;
119  did_cross=true;
120  }
121  else if (pre_vol_name == volume_name){
122  GoingIn=false;
123  did_cross=true;
124 
125  }
126 
127  }
128  }
129  return did_cross;
130  //still need to compute the cosine of the direction
131 }
133 //
134 G4bool G4AdjointCrossSurfChecker::GoingInOrOutOfaVolumeByExtSurface(const G4Step* aStep,const G4String& volume_name, const G4String& mother_logical_vol_name, G4double& , G4bool& GoingIn) //from external surface
135 {
136  G4bool step_at_boundary = (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
137  G4bool did_cross =false;
138  if (step_at_boundary){
139  const G4VTouchable* postStepTouchable = aStep->GetPostStepPoint()->GetTouchable();
140  const G4VTouchable* preStepTouchable = aStep->GetPreStepPoint()->GetTouchable();
141  if (preStepTouchable && postStepTouchable && postStepTouchable->GetVolume() && preStepTouchable->GetVolume()){
142  G4String post_vol_name = postStepTouchable->GetVolume()->GetName();
143  G4String post_log_vol_name = postStepTouchable->GetVolume()->GetLogicalVolume()->GetName();
144  G4String pre_vol_name = preStepTouchable->GetVolume()->GetName();
145  G4String pre_log_vol_name = preStepTouchable->GetVolume()->GetLogicalVolume()->GetName();
146  if (post_vol_name == volume_name && pre_log_vol_name == mother_logical_vol_name){
147  GoingIn=true;
148  did_cross=true;
149  }
150  else if (pre_vol_name == volume_name && post_log_vol_name == mother_logical_vol_name ){
151  GoingIn=false;
152  did_cross=true;
153 
154  }
155 
156  }
157  }
158  return did_cross;
159  //still need to compute the cosine of the direction
160 }
162 //
163 G4bool G4AdjointCrossSurfChecker::CrossingAGivenRegisteredSurface(const G4Step* aStep,const G4String& surface_name,G4ThreeVector& crossing_pos, G4double& cos_to_surface, G4bool& GoingIn)
164 {
165  G4int ind = FindRegisteredSurface(surface_name);
166  G4bool did_cross = false;
167  if (ind >=0){
168  did_cross = CrossingAGivenRegisteredSurface(aStep, ind, crossing_pos,cos_to_surface, GoingIn);
169  }
170  return did_cross;
171 }
173 //
174 G4bool G4AdjointCrossSurfChecker::CrossingAGivenRegisteredSurface(const G4Step* aStep, int ind,G4ThreeVector& crossing_pos, G4double& cos_to_surface, G4bool& GoingIn)
175 {
176  G4String surf_type = ListOfSurfaceType[ind];
177  G4double radius = ListOfSphereRadius[ind];
178  G4ThreeVector center = ListOfSphereCenter[ind];
179  G4String vol1 = ListOfVol1Name[ind];
180  G4String vol2 = ListOfVol2Name[ind];
181 
182  G4bool did_cross = false;
183  if (surf_type == "Sphere"){
184  did_cross = CrossingASphere(aStep, radius, center,crossing_pos, cos_to_surface, GoingIn);
185  }
186  else if (surf_type == "ExternalSurfaceOfAVolume"){
187 
188  did_cross = GoingInOrOutOfaVolumeByExtSurface(aStep, vol1, vol2, cos_to_surface, GoingIn);
189  crossing_pos= aStep->GetPostStepPoint()->GetPosition();
190 
191  }
192  else if (surf_type == "BoundaryBetweenTwoVolumes"){
193  did_cross = CrossingAnInterfaceBetweenTwoVolumes(aStep, vol1, vol2,crossing_pos, cos_to_surface, GoingIn);
194  }
195  return did_cross;
196 
197 
198 }
200 //
201 //
202 G4bool G4AdjointCrossSurfChecker::CrossingOneOfTheRegisteredSurface(const G4Step* aStep,G4String& surface_name,G4ThreeVector& crossing_pos, G4double& cos_to_surface, G4bool& GoingIn)
203 {
204  for (size_t i=0;i <ListOfSurfaceName.size();i++){
205  if (CrossingAGivenRegisteredSurface(aStep, int(i),crossing_pos, cos_to_surface, GoingIn)){
206  surface_name = ListOfSurfaceName[i];
207  return true;
208  }
209  }
210  return false;
211 }
213 //
214 G4bool G4AdjointCrossSurfChecker::CrossingAnInterfaceBetweenTwoVolumes(const G4Step* aStep,const G4String& vol1_name,const G4String& vol2_name,G4ThreeVector& , G4double& , G4bool& GoingIn)
215 {
216  G4bool step_at_boundary = (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
217  G4bool did_cross =false;
218  if (step_at_boundary){
219  const G4VTouchable* postStepTouchable = aStep->GetPostStepPoint()->GetTouchable();
220  const G4VTouchable* preStepTouchable = aStep->GetPreStepPoint()->GetTouchable();
221  if (preStepTouchable && postStepTouchable){
222 
223  G4String post_vol_name = postStepTouchable->GetVolume()->GetName();
224  if (post_vol_name =="") post_vol_name = postStepTouchable->GetVolume()->GetLogicalVolume()->GetName();
225  G4String pre_vol_name = preStepTouchable->GetVolume()->GetName();
226  if (pre_vol_name =="") pre_vol_name = preStepTouchable->GetVolume()->GetLogicalVolume()->GetName();
227 
228 
229  if ( pre_vol_name == vol1_name && post_vol_name == vol2_name){
230  GoingIn=true;
231  did_cross=true;
232  }
233  else if (pre_vol_name == vol2_name && post_vol_name == vol1_name){
234  GoingIn=false;
235  did_cross=true;
236  }
237 
238  }
239  }
240  return did_cross;
241  //still need to compute the cosine of the direction
242 }
243 
245 //
246 G4bool G4AdjointCrossSurfChecker::AddaSphericalSurface(const G4String& SurfaceName, G4double radius, G4ThreeVector pos, G4double& Area)
247 {
248  G4int ind = FindRegisteredSurface(SurfaceName);
249  Area= 4.*pi*radius*radius;
250  if (ind>=0) {
251  ListOfSurfaceType[ind]="Sphere";
252  ListOfSphereRadius[ind]=radius;
253  ListOfSphereCenter[ind]=pos;
254  ListOfVol1Name[ind]="";
255  ListOfVol2Name[ind]="";
256  AreaOfSurface[ind]=Area;
257  }
258  else {
259  ListOfSurfaceName.push_back(SurfaceName);
260  ListOfSurfaceType.push_back("Sphere");
261  ListOfSphereRadius.push_back(radius);
262  ListOfSphereCenter.push_back(pos);
263  ListOfVol1Name.push_back("");
264  ListOfVol2Name.push_back("");
265  AreaOfSurface.push_back(Area);
266  }
267  return true;
268 }
270 //
271 G4bool G4AdjointCrossSurfChecker::AddaSphericalSurfaceWithCenterAtTheCenterOfAVolume(const G4String& SurfaceName, G4double radius, const G4String& volume_name, G4ThreeVector & center, G4double& area)
272 {
273 
274  G4VPhysicalVolume* thePhysicalVolume = 0;
275  G4PhysicalVolumeStore* thePhysVolStore =G4PhysicalVolumeStore::GetInstance();
276  for ( unsigned int i=0; i< thePhysVolStore->size();i++){
277  if ((*thePhysVolStore)[i]->GetName() == volume_name){
278  thePhysicalVolume = (*thePhysVolStore)[i];
279  };
280 
281  }
282  if (thePhysicalVolume){
283  G4VPhysicalVolume* daughter =thePhysicalVolume;
284  G4LogicalVolume* mother = thePhysicalVolume->GetMotherLogical();
285  G4AffineTransform theTransformationFromPhysVolToWorld = G4AffineTransform();
286  while (mother){
287  theTransformationFromPhysVolToWorld *=
288  G4AffineTransform(daughter->GetFrameRotation(),daughter->GetObjectTranslation());
289  /*G4cout<<"Mother "<<mother->GetName()<<std::endl;
290  G4cout<<"Daughter "<<daughter->GetName()<<std::endl;
291  G4cout<<daughter->GetObjectTranslation()<<std::endl;
292  G4cout<<theTransformationFromPhysVolToWorld.NetTranslation()<<std::endl;*/
293  for ( unsigned int i=0; i< thePhysVolStore->size();i++){
294  if ((*thePhysVolStore)[i]->GetLogicalVolume() == mother){
295  daughter = (*thePhysVolStore)[i];
296  mother =daughter->GetMotherLogical();
297  break;
298  };
299  }
300 
301  }
302  center = theTransformationFromPhysVolToWorld.NetTranslation();
303  G4cout<<"Center of the spherical surface is at the position: "<<center/cm<<" cm"<<std::endl;
304 
305  }
306  else {
307  G4cout<<"The physical volume with name "<<volume_name<<" does not exist!!"<<std::endl;
308  return false;
309  }
310  return AddaSphericalSurface(SurfaceName, radius, center, area);
311 
312 }
314 //
315 G4bool G4AdjointCrossSurfChecker::AddanExtSurfaceOfAvolume(const G4String& SurfaceName, const G4String& volume_name, G4double& Area)
316 {
317  G4int ind = FindRegisteredSurface(SurfaceName);
318 
319  G4VPhysicalVolume* thePhysicalVolume = 0;
320  G4PhysicalVolumeStore* thePhysVolStore =G4PhysicalVolumeStore::GetInstance();
321  for ( unsigned int i=0; i< thePhysVolStore->size();i++){
322  if ((*thePhysVolStore)[i]->GetName() == volume_name){
323  thePhysicalVolume = (*thePhysVolStore)[i];
324  };
325 
326  }
327  if (!thePhysicalVolume){
328  G4cout<<"The physical volume with name "<<volume_name<<" does not exist!!"<<std::endl;
329  return false;
330  }
331  Area = thePhysicalVolume->GetLogicalVolume()->GetSolid()->GetSurfaceArea();
332  G4String mother_vol_name = "";
333  G4LogicalVolume* theMother = thePhysicalVolume->GetMotherLogical();
334 
335  if (theMother) mother_vol_name= theMother->GetName();
336  if (ind>=0) {
337  ListOfSurfaceType[ind]="ExternalSurfaceOfAVolume";
338  ListOfSphereRadius[ind]=0.;
339  ListOfSphereCenter[ind]=G4ThreeVector(0.,0.,0.);
340  ListOfVol1Name[ind]=volume_name;
341  ListOfVol2Name[ind]=mother_vol_name;
342  AreaOfSurface[ind]=Area;
343  }
344  else {
345  ListOfSurfaceName.push_back(SurfaceName);
346  ListOfSurfaceType.push_back("ExternalSurfaceOfAVolume");
347  ListOfSphereRadius.push_back(0.);
348  ListOfSphereCenter.push_back(G4ThreeVector(0.,0.,0.));
349  ListOfVol1Name.push_back(volume_name);
350  ListOfVol2Name.push_back(mother_vol_name);
351  AreaOfSurface.push_back(Area);
352  }
353  return true;
354 }
356 //
357 G4bool G4AdjointCrossSurfChecker::AddanInterfaceBetweenTwoVolumes(const G4String& SurfaceName, const G4String& volume_name1, const G4String& volume_name2,G4double& Area)
358 {
359  G4int ind = FindRegisteredSurface(SurfaceName);
360  Area=-1.; //the way to compute the surface is not known yet
361  if (ind>=0) {
362  ListOfSurfaceType[ind]="BoundaryBetweenTwoVolumes";
363  ListOfSphereRadius[ind]=0.;
364  ListOfSphereCenter[ind]=G4ThreeVector(0.,0.,0.);
365  ListOfVol1Name[ind]=volume_name1;
366  ListOfVol2Name[ind]=volume_name2;
367  AreaOfSurface[ind]=Area;
368 
369  }
370  else {
371  ListOfSurfaceName.push_back(SurfaceName);
372  ListOfSurfaceType.push_back("BoundaryBetweenTwoVolumes");
373  ListOfSphereRadius.push_back(0.);
374  ListOfSphereCenter.push_back(G4ThreeVector(0.,0.,0.));
375  ListOfVol1Name.push_back(volume_name1);
376  ListOfVol2Name.push_back(volume_name2);
377  AreaOfSurface.push_back(Area);
378  }
379  return true;
380 }
382 //
383 void G4AdjointCrossSurfChecker::ClearListOfSelectedSurface()
384 {
385  ListOfSurfaceName.clear();
386  ListOfSurfaceType.clear();
387  ListOfSphereRadius.clear();
388  ListOfSphereCenter.clear();
389  ListOfVol1Name.clear();
390  ListOfVol2Name.clear();
391 }
393 //
394 G4int G4AdjointCrossSurfChecker::FindRegisteredSurface(const G4String& name)
395 {
396  G4int ind=-1;
397  for (size_t i = 0; i<ListOfSurfaceName.size();i++){
398  if (name == ListOfSurfaceName[i]) {
399  ind = int (i);
400  return ind;
401  }
402  }
403  return ind;
404 }
405