source: trunk/source/tracking/src/G4AdjointCrossSurfChecker.cc@ 1357

Last change on this file since 1357 was 1337, checked in by garnier, 15 years ago

tag geant4.9.4 beta 1 + modifs locales
File size: 15.3 KB
Line 
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,v 1.2 2009/11/18 18:04:11 gcosmo Exp $
27// GEANT4 tag $Name: geant4-09-04-beta-01 $
28//
29/////////////////////////////////////////////////////////////////////////////
30// Class Name: G4AdjointCrossSurfChecker
31// Author: L. Desorgher
32// Organisation: SpaceIT GmbH
33// Contract: ESA contract 21435/08/NL/AT
34// Customer: ESA/ESTEC
35/////////////////////////////////////////////////////////////////////////////
36
37#include"G4AdjointCrossSurfChecker.hh"
38#include"G4Step.hh"
39#include"G4StepPoint.hh"
40#include"G4PhysicalVolumeStore.hh"
41#include"G4VSolid.hh"
42#include"G4AffineTransform.hh"
43
44//////////////////////////////////////////////////////////////////////////////
45//
46G4AdjointCrossSurfChecker* G4AdjointCrossSurfChecker::instance = 0;
47
48//////////////////////////////////////////////////////////////////////////////
49//
50G4AdjointCrossSurfChecker::G4AdjointCrossSurfChecker()
51{;
52}
53///////////////////////////////////////////////////////////////////////////////
54//
55G4AdjointCrossSurfChecker::~G4AdjointCrossSurfChecker()
56{
57 delete instance;
58}
59////////////////////////////////////////////////////////////////////////////////
60//
61G4AdjointCrossSurfChecker* G4AdjointCrossSurfChecker::GetInstance()
62{
63 if (!instance) instance = new G4AdjointCrossSurfChecker();
64 return instance;
65}
66/////////////////////////////////////////////////////////////////////////////////
67//
68G4bool 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}
104/////////////////////////////////////////////////////////////////////////////////
105//
106G4bool 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}
132/////////////////////////////////////////////////////////////////////////////////
133//
134G4bool 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}
161/////////////////////////////////////////////////////////////////////////////////
162//
163G4bool 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}
172/////////////////////////////////////////////////////////////////////////////////
173//
174G4bool 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}
199/////////////////////////////////////////////////////////////////////////////////
200//
201//
202G4bool 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}
212/////////////////////////////////////////////////////////////////////////////////
213//
214G4bool 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
244/////////////////////////////////////////////////////////////////////////////////
245//
246G4bool 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}
269/////////////////////////////////////////////////////////////////////////////////
270//
271G4bool 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 G4PhysicalVolumeStore* thePhysVolStore =G4PhysicalVolumeStore::GetInstance();
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}
314/////////////////////////////////////////////////////////////////////////////////
315//
316G4bool 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}
356/////////////////////////////////////////////////////////////////////////////////
357//
358G4bool 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}
382/////////////////////////////////////////////////////////////////////////////////
383//
384void G4AdjointCrossSurfChecker::ClearListOfSelectedSurface()
385{
386 ListOfSurfaceName.clear();
387 ListOfSurfaceType.clear();
388 ListOfSphereRadius.clear();
389 ListOfSphereCenter.clear();
390 ListOfVol1Name.clear();
391 ListOfVol2Name.clear();
392}
393/////////////////////////////////////////////////////////////////////////////////
394//
395G4int 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
Note: See TracBrowser for help on using the repository browser.