Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4AdjointCrossSurfChecker.cc
Go to the documentation of this file.
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 // $Id: G4AdjointCrossSurfChecker.cc 67009 2013-01-29 16:00:21Z gcosmo $
27 //
29 // Class Name: G4AdjointCrossSurfChecker
30 // Author: L. Desorgher
31 // Organisation: SpaceIT GmbH
32 // Contract: ESA contract 21435/08/NL/AT
33 // Customer: ESA/ESTEC
35 
36 #include "G4AdjointCrossSurfChecker.hh"
37 #include "G4PhysicalConstants.hh"
38 #include "G4SystemOfUnits.hh"
39 #include "G4Step.hh"
40 #include "G4StepPoint.hh"
41 #include "G4PhysicalVolumeStore.hh"
42 #include "G4VSolid.hh"
43 #include "G4AffineTransform.hh"
44 
46 //
47 G4AdjointCrossSurfChecker* G4AdjointCrossSurfChecker::instance = 0;
48 
50 //
51 G4AdjointCrossSurfChecker::G4AdjointCrossSurfChecker()
52 {;
53 }
55 //
56 G4AdjointCrossSurfChecker::~G4AdjointCrossSurfChecker()
57 {
58  delete instance;
59 }
61 //
62 G4AdjointCrossSurfChecker* G4AdjointCrossSurfChecker::GetInstance()
63 {
64  if (!instance) instance = new G4AdjointCrossSurfChecker();
65  return instance;
66 }
68 //
69 G4bool G4AdjointCrossSurfChecker::CrossingASphere(const G4Step* aStep,G4double sphere_radius, G4ThreeVector sphere_center,G4ThreeVector& crossing_pos, G4double& cos_th , G4bool& GoingIn)
70 {
71  G4ThreeVector pos1= aStep->GetPreStepPoint()->GetPosition() - sphere_center;
72  G4ThreeVector pos2= aStep->GetPostStepPoint()->GetPosition() - sphere_center;
73  G4double r1= pos1.mag();
74  G4double r2= pos2.mag();
75  G4bool did_cross =false;
76 
77  if (r1<=sphere_radius && r2>sphere_radius){
78  did_cross=true;
79  GoingIn=false;
80  }
81  else if (r2<=sphere_radius && r1>sphere_radius){
82  did_cross=true;
83  GoingIn=true;
84  }
85 
86  if (did_cross) {
87 
88  G4ThreeVector dr=pos2-pos1;
89  G4double r12 = r1*r1;
90  G4double rdr = dr.mag();
91  G4double a,b,c,d;
92  a = rdr*rdr;
93  b = 2.*pos1.dot(dr);
94  c = r12-sphere_radius*sphere_radius;
95  d=std::sqrt(b*b-4.*a*c);
96  G4double l= (-b+d)/2./a;
97  if (l > 1.) l=(-b-d)/2./a;
98  crossing_pos=pos1+l*dr;
99  cos_th = std::abs(dr.cosTheta(crossing_pos));
100 
101 
102  }
103  return did_cross;
104 }
106 //
107 G4bool G4AdjointCrossSurfChecker::GoingInOrOutOfaVolume(const G4Step* aStep,const G4String& volume_name, G4double& , G4bool& GoingIn) //from external surface
108 {
109  G4bool step_at_boundary = (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
110  G4bool did_cross =false;
111  if (step_at_boundary){
112  const G4VTouchable* postStepTouchable = aStep->GetPostStepPoint()->GetTouchable();
113  const G4VTouchable* preStepTouchable = aStep->GetPreStepPoint()->GetTouchable();
114  if (preStepTouchable && postStepTouchable && postStepTouchable->GetVolume() && preStepTouchable->GetVolume()){
115  G4String post_vol_name = postStepTouchable->GetVolume()->GetName();
116  G4String pre_vol_name = preStepTouchable->GetVolume()->GetName();
117 
118  if (post_vol_name == volume_name ){
119  GoingIn=true;
120  did_cross=true;
121  }
122  else if (pre_vol_name == volume_name){
123  GoingIn=false;
124  did_cross=true;
125 
126  }
127 
128  }
129  }
130  return did_cross;
131  //still need to compute the cosine of the direction
132 }
134 //
135 G4bool G4AdjointCrossSurfChecker::GoingInOrOutOfaVolumeByExtSurface(const G4Step* aStep,const G4String& volume_name, const G4String& mother_logical_vol_name, G4double& , G4bool& GoingIn) //from external surface
136 {
137  G4bool step_at_boundary = (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
138  G4bool did_cross =false;
139  if (step_at_boundary){
140  const G4VTouchable* postStepTouchable = aStep->GetPostStepPoint()->GetTouchable();
141  const G4VTouchable* preStepTouchable = aStep->GetPreStepPoint()->GetTouchable();
142  if (preStepTouchable && postStepTouchable && postStepTouchable->GetVolume() && preStepTouchable->GetVolume()){
143  G4String post_vol_name = postStepTouchable->GetVolume()->GetName();
144  G4String post_log_vol_name = postStepTouchable->GetVolume()->GetLogicalVolume()->GetName();
145  G4String pre_vol_name = preStepTouchable->GetVolume()->GetName();
146  G4String pre_log_vol_name = preStepTouchable->GetVolume()->GetLogicalVolume()->GetName();
147  if (post_vol_name == volume_name && pre_log_vol_name == mother_logical_vol_name){
148  GoingIn=true;
149  did_cross=true;
150  }
151  else if (pre_vol_name == volume_name && post_log_vol_name == mother_logical_vol_name ){
152  GoingIn=false;
153  did_cross=true;
154 
155  }
156 
157  }
158  }
159  return did_cross;
160  //still need to compute the cosine of the direction
161 }
163 //
164 G4bool G4AdjointCrossSurfChecker::CrossingAGivenRegisteredSurface(const G4Step* aStep,const G4String& surface_name,G4ThreeVector& crossing_pos, G4double& cos_to_surface, G4bool& GoingIn)
165 {
166  G4int ind = FindRegisteredSurface(surface_name);
167  G4bool did_cross = false;
168  if (ind >=0){
169  did_cross = CrossingAGivenRegisteredSurface(aStep, ind, crossing_pos,cos_to_surface, GoingIn);
170  }
171  return did_cross;
172 }
174 //
175 G4bool G4AdjointCrossSurfChecker::CrossingAGivenRegisteredSurface(const G4Step* aStep, int ind,G4ThreeVector& crossing_pos, G4double& cos_to_surface, G4bool& GoingIn)
176 {
177  G4String surf_type = ListOfSurfaceType[ind];
178  G4double radius = ListOfSphereRadius[ind];
179  G4ThreeVector center = ListOfSphereCenter[ind];
180  G4String vol1 = ListOfVol1Name[ind];
181  G4String vol2 = ListOfVol2Name[ind];
182 
183  G4bool did_cross = false;
184  if (surf_type == "Sphere"){
185  did_cross = CrossingASphere(aStep, radius, center,crossing_pos, cos_to_surface, GoingIn);
186  }
187  else if (surf_type == "ExternalSurfaceOfAVolume"){
188 
189  did_cross = GoingInOrOutOfaVolumeByExtSurface(aStep, vol1, vol2, cos_to_surface, GoingIn);
190  crossing_pos= aStep->GetPostStepPoint()->GetPosition();
191 
192  }
193  else if (surf_type == "BoundaryBetweenTwoVolumes"){
194  did_cross = CrossingAnInterfaceBetweenTwoVolumes(aStep, vol1, vol2,crossing_pos, cos_to_surface, GoingIn);
195  }
196  return did_cross;
197 
198 
199 }
201 //
202 //
203 G4bool G4AdjointCrossSurfChecker::CrossingOneOfTheRegisteredSurface(const G4Step* aStep,G4String& surface_name,G4ThreeVector& crossing_pos, G4double& cos_to_surface, G4bool& GoingIn)
204 {
205  for (size_t i=0;i <ListOfSurfaceName.size();i++){
206  if (CrossingAGivenRegisteredSurface(aStep, int(i),crossing_pos, cos_to_surface, GoingIn)){
207  surface_name = ListOfSurfaceName[i];
208  return true;
209  }
210  }
211  return false;
212 }
214 //
215 G4bool G4AdjointCrossSurfChecker::CrossingAnInterfaceBetweenTwoVolumes(const G4Step* aStep,const G4String& vol1_name,const G4String& vol2_name,G4ThreeVector& , G4double& , G4bool& GoingIn)
216 {
217  G4bool step_at_boundary = (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
218  G4bool did_cross =false;
219  if (step_at_boundary){
220  const G4VTouchable* postStepTouchable = aStep->GetPostStepPoint()->GetTouchable();
221  const G4VTouchable* preStepTouchable = aStep->GetPreStepPoint()->GetTouchable();
222  if (preStepTouchable && postStepTouchable){
223 
224  G4String post_vol_name = postStepTouchable->GetVolume()->GetName();
225  if (post_vol_name =="") post_vol_name = postStepTouchable->GetVolume()->GetLogicalVolume()->GetName();
226  G4String pre_vol_name = preStepTouchable->GetVolume()->GetName();
227  if (pre_vol_name =="") pre_vol_name = preStepTouchable->GetVolume()->GetLogicalVolume()->GetName();
228 
229 
230  if ( pre_vol_name == vol1_name && post_vol_name == vol2_name){
231  GoingIn=true;
232  did_cross=true;
233  }
234  else if (pre_vol_name == vol2_name && post_vol_name == vol1_name){
235  GoingIn=false;
236  did_cross=true;
237  }
238 
239  }
240  }
241  return did_cross;
242  //still need to compute the cosine of the direction
243 }
244 
246 //
247 G4bool G4AdjointCrossSurfChecker::AddaSphericalSurface(const G4String& SurfaceName, G4double radius, G4ThreeVector pos, G4double& Area)
248 {
249  G4int ind = FindRegisteredSurface(SurfaceName);
250  Area= 4.*pi*radius*radius;
251  if (ind>=0) {
252  ListOfSurfaceType[ind]="Sphere";
253  ListOfSphereRadius[ind]=radius;
254  ListOfSphereCenter[ind]=pos;
255  ListOfVol1Name[ind]="";
256  ListOfVol2Name[ind]="";
257  AreaOfSurface[ind]=Area;
258  }
259  else {
260  ListOfSurfaceName.push_back(SurfaceName);
261  ListOfSurfaceType.push_back("Sphere");
262  ListOfSphereRadius.push_back(radius);
263  ListOfSphereCenter.push_back(pos);
264  ListOfVol1Name.push_back("");
265  ListOfVol2Name.push_back("");
266  AreaOfSurface.push_back(Area);
267  }
268  return true;
269 }
271 //
272 G4bool G4AdjointCrossSurfChecker::AddaSphericalSurfaceWithCenterAtTheCenterOfAVolume(const G4String& SurfaceName, G4double radius, const G4String& volume_name, G4ThreeVector & center, G4double& area)
273 {
274 
275  G4VPhysicalVolume* thePhysicalVolume = 0;
276  G4PhysicalVolumeStore* thePhysVolStore =G4PhysicalVolumeStore::GetInstance();
277  for ( unsigned int i=0; i< thePhysVolStore->size();i++){
278  if ((*thePhysVolStore)[i]->GetName() == volume_name){
279  thePhysicalVolume = (*thePhysVolStore)[i];
280  };
281 
282  }
283  if (thePhysicalVolume){
284  G4VPhysicalVolume* daughter =thePhysicalVolume;
285  G4LogicalVolume* mother = thePhysicalVolume->GetMotherLogical();
286  G4AffineTransform theTransformationFromPhysVolToWorld = G4AffineTransform();
287  while (mother){
288  theTransformationFromPhysVolToWorld *=
289  G4AffineTransform(daughter->GetFrameRotation(),daughter->GetObjectTranslation());
290  /*G4cout<<"Mother "<<mother->GetName()<<std::endl;
291  G4cout<<"Daughter "<<daughter->GetName()<<std::endl;
292  G4cout<<daughter->GetObjectTranslation()<<std::endl;
293  G4cout<<theTransformationFromPhysVolToWorld.NetTranslation()<<std::endl;*/
294  for ( unsigned int i=0; i< thePhysVolStore->size();i++){
295  if ((*thePhysVolStore)[i]->GetLogicalVolume() == mother){
296  daughter = (*thePhysVolStore)[i];
297  mother =daughter->GetMotherLogical();
298  break;
299  };
300  }
301 
302  }
303  center = theTransformationFromPhysVolToWorld.NetTranslation();
304  G4cout<<"Center of the spherical surface is at the position: "<<center/cm<<" cm"<<std::endl;
305 
306  }
307  else {
308  G4cout<<"The physical volume with name "<<volume_name<<" does not exist!!"<<std::endl;
309  return false;
310  }
311  return AddaSphericalSurface(SurfaceName, radius, center, area);
312 
313 }
315 //
316 G4bool G4AdjointCrossSurfChecker::AddanExtSurfaceOfAvolume(const G4String& SurfaceName, const G4String& volume_name, G4double& Area)
317 {
318  G4int ind = FindRegisteredSurface(SurfaceName);
319 
320  G4VPhysicalVolume* thePhysicalVolume = 0;
321  G4PhysicalVolumeStore* thePhysVolStore =G4PhysicalVolumeStore::GetInstance();
322  for ( unsigned int i=0; i< thePhysVolStore->size();i++){
323  if ((*thePhysVolStore)[i]->GetName() == volume_name){
324  thePhysicalVolume = (*thePhysVolStore)[i];
325  };
326 
327  }
328  if (!thePhysicalVolume){
329  G4cout<<"The physical volume with name "<<volume_name<<" does not exist!!"<<std::endl;
330  return false;
331  }
332  Area = thePhysicalVolume->GetLogicalVolume()->GetSolid()->GetSurfaceArea();
333  G4String mother_vol_name = "";
334  G4LogicalVolume* theMother = thePhysicalVolume->GetMotherLogical();
335 
336  if (theMother) mother_vol_name= theMother->GetName();
337  if (ind>=0) {
338  ListOfSurfaceType[ind]="ExternalSurfaceOfAVolume";
339  ListOfSphereRadius[ind]=0.;
340  ListOfSphereCenter[ind]=G4ThreeVector(0.,0.,0.);
341  ListOfVol1Name[ind]=volume_name;
342  ListOfVol2Name[ind]=mother_vol_name;
343  AreaOfSurface[ind]=Area;
344  }
345  else {
346  ListOfSurfaceName.push_back(SurfaceName);
347  ListOfSurfaceType.push_back("ExternalSurfaceOfAVolume");
348  ListOfSphereRadius.push_back(0.);
349  ListOfSphereCenter.push_back(G4ThreeVector(0.,0.,0.));
350  ListOfVol1Name.push_back(volume_name);
351  ListOfVol2Name.push_back(mother_vol_name);
352  AreaOfSurface.push_back(Area);
353  }
354  return true;
355 }
357 //
358 G4bool G4AdjointCrossSurfChecker::AddanInterfaceBetweenTwoVolumes(const G4String& SurfaceName, const G4String& volume_name1, const G4String& volume_name2,G4double& Area)
359 {
360  G4int ind = FindRegisteredSurface(SurfaceName);
361  Area=-1.; //the way to compute the surface is not known yet
362  if (ind>=0) {
363  ListOfSurfaceType[ind]="BoundaryBetweenTwoVolumes";
364  ListOfSphereRadius[ind]=0.;
365  ListOfSphereCenter[ind]=G4ThreeVector(0.,0.,0.);
366  ListOfVol1Name[ind]=volume_name1;
367  ListOfVol2Name[ind]=volume_name2;
368  AreaOfSurface[ind]=Area;
369 
370  }
371  else {
372  ListOfSurfaceName.push_back(SurfaceName);
373  ListOfSurfaceType.push_back("BoundaryBetweenTwoVolumes");
374  ListOfSphereRadius.push_back(0.);
375  ListOfSphereCenter.push_back(G4ThreeVector(0.,0.,0.));
376  ListOfVol1Name.push_back(volume_name1);
377  ListOfVol2Name.push_back(volume_name2);
378  AreaOfSurface.push_back(Area);
379  }
380  return true;
381 }
383 //
384 void G4AdjointCrossSurfChecker::ClearListOfSelectedSurface()
385 {
386  ListOfSurfaceName.clear();
387  ListOfSurfaceType.clear();
388  ListOfSphereRadius.clear();
389  ListOfSphereCenter.clear();
390  ListOfVol1Name.clear();
391  ListOfVol2Name.clear();
392 }
394 //
395 G4int G4AdjointCrossSurfChecker::FindRegisteredSurface(const G4String& name)
396 {
397  G4int ind=-1;
398  for (size_t i = 0; i<ListOfSurfaceName.size();i++){
399  if (name == ListOfSurfaceName[i]) {
400  ind = int (i);
401  return ind;
402  }
403  }
404  return ind;
405 }
406 
Generated on Sat May 25 2013 14:34:12 for Geant4 by   doxygen 1.8.4