Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4AdjointPosOnPhysVolGenerator.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$
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 "G4AdjointPosOnPhysVolGenerator.hh"
37 #include "G4VSolid.hh"
38 #include "G4VoxelLimits.hh"
39 #include "G4AffineTransform.hh"
40 #include "Randomize.hh"
41 #include "G4VPhysicalVolume.hh"
42 #include "G4PhysicalVolumeStore.hh"
43 #include "G4LogicalVolumeStore.hh"
44 
45 G4AdjointPosOnPhysVolGenerator* G4AdjointPosOnPhysVolGenerator::theInstance = 0;
46 
48 //
49 G4AdjointPosOnPhysVolGenerator* G4AdjointPosOnPhysVolGenerator::GetInstance()
50 {
51  if(theInstance == 0) {
52  static G4AdjointPosOnPhysVolGenerator manager;
53  theInstance = &manager;
54  }
55  return theInstance;
56 }
57 
59 //
60 G4AdjointPosOnPhysVolGenerator::~G4AdjointPosOnPhysVolGenerator()
61 {
62 }
63 
65 //
66 G4AdjointPosOnPhysVolGenerator::G4AdjointPosOnPhysVolGenerator()
67  : theSolid(0), thePhysicalVolume(0),
68  UseSphere(true), ModelOfSurfaceSource("OnSolid"),
69  AreaOfExtSurfaceOfThePhysicalVolume(0.), CosThDirComparedToNormal(0.)
70 {
71 }
72 
74 //
75 G4VPhysicalVolume* G4AdjointPosOnPhysVolGenerator::DefinePhysicalVolume(const G4String& aName)
76 {
77  thePhysicalVolume = 0;
78  theSolid =0;
79  G4PhysicalVolumeStore* thePhysVolStore =G4PhysicalVolumeStore::GetInstance();
80  for ( unsigned int i=0; i< thePhysVolStore->size();i++){
81  G4String vol_name =(*thePhysVolStore)[i]->GetName();
82  if (vol_name == ""){
83  vol_name = (*thePhysVolStore)[i]->GetLogicalVolume()->GetName();
84  }
85  if (vol_name == aName){
86  thePhysicalVolume = (*thePhysVolStore)[i];
87  }
88  }
89  if (thePhysicalVolume){
90  theSolid = thePhysicalVolume->GetLogicalVolume()->GetSolid();
91  ComputeTransformationFromPhysVolToWorld();
92  /*AreaOfExtSurfaceOfThePhysicalVolume=ComputeAreaOfExtSurface(1.e-3);
93  G4cout<<"Monte Carlo Estimate of the area of the external surface :"<<AreaOfExtSurfaceOfThePhysicalVolume/m/m<<" m2"<<std::endl;*/
94  }
95  else {
96  G4cout<<"The physical volume with name "<<aName<<" does not exist!!"<<std::endl;
97  G4cout<<"Before generating a source on an external surface of a volume you should select another physical volume"<<std::endl;
98  }
99  return thePhysicalVolume;
100 }
102 //
103 void G4AdjointPosOnPhysVolGenerator::DefinePhysicalVolume1(const G4String& aName)
104 {
105  thePhysicalVolume = DefinePhysicalVolume(aName);
106 }
108 //
109 G4double G4AdjointPosOnPhysVolGenerator::ComputeAreaOfExtSurface()
110 {
111  return ComputeAreaOfExtSurface(theSolid);
112 }
114 //
115 G4double G4AdjointPosOnPhysVolGenerator::ComputeAreaOfExtSurface(G4int NStats)
116 {
117  return ComputeAreaOfExtSurface(theSolid,NStats);
118 }
120 //
121 G4double G4AdjointPosOnPhysVolGenerator::ComputeAreaOfExtSurface(G4double eps)
122 {
123  return ComputeAreaOfExtSurface(theSolid,eps);
124 }
126 //
127 G4double G4AdjointPosOnPhysVolGenerator::ComputeAreaOfExtSurface(G4VSolid* aSolid)
128 {
129  return ComputeAreaOfExtSurface(aSolid,1.e-3);
130 }
132 //
133 G4double G4AdjointPosOnPhysVolGenerator::ComputeAreaOfExtSurface(G4VSolid* aSolid,G4int NStats)
134 {
135  if (ModelOfSurfaceSource == "OnSolid" ){
136  if (UseSphere){
137  return ComputeAreaOfExtSurfaceStartingFromSphere(aSolid,NStats);
138  }
139  else {
140  return ComputeAreaOfExtSurfaceStartingFromBox(aSolid,NStats);
141  }
142  }
143  else {
144  G4ThreeVector p,dir;
145  if (ModelOfSurfaceSource == "ExternalSphere" ) return GenerateAPositionOnASphereBoundary(aSolid, p,dir);
146  return GenerateAPositionOnABoxBoundary(aSolid, p,dir);
147  }
148 }
150 //
151 G4double G4AdjointPosOnPhysVolGenerator::ComputeAreaOfExtSurface(G4VSolid* aSolid,G4double eps)
152 {
153  G4int Nstats = G4int(1./(eps*eps));
154  return ComputeAreaOfExtSurface(aSolid,Nstats);
155 }
157 void G4AdjointPosOnPhysVolGenerator::GenerateAPositionOnTheExtSurfaceOfASolid(G4VSolid* aSolid,G4ThreeVector& p, G4ThreeVector& direction)
158 {
159  if (ModelOfSurfaceSource == "OnSolid" ){
160  GenerateAPositionOnASolidBoundary(aSolid, p,direction);
161  return;
162  }
163  if (ModelOfSurfaceSource == "ExternalSphere" ) {
164  GenerateAPositionOnASphereBoundary(aSolid, p, direction);
165  return;
166  }
167  GenerateAPositionOnABoxBoundary(aSolid, p, direction);
168  return;
169 }
171 void G4AdjointPosOnPhysVolGenerator::GenerateAPositionOnTheExtSurfaceOfTheSolid(G4ThreeVector& p, G4ThreeVector& direction)
172 {
173  GenerateAPositionOnTheExtSurfaceOfASolid(theSolid,p,direction);
174 }
176 //
177 G4double G4AdjointPosOnPhysVolGenerator::ComputeAreaOfExtSurfaceStartingFromBox(G4VSolid* aSolid,G4int Nstat)
178 {
179  G4double area=1.;
180  G4int i=0;
181  G4int j=0;
182  while (i<Nstat){
183  G4ThreeVector p, direction;
184  area = GenerateAPositionOnABoxBoundary( aSolid,p, direction);
185  G4double dist_to_in = aSolid->DistanceToIn(p,direction);
186  if (dist_to_in<kInfinity/2.) i++;
187  j++;
188  }
189  area=area*double(i)/double(j);
190  return area;
191 }
193 //
194 G4double G4AdjointPosOnPhysVolGenerator::ComputeAreaOfExtSurfaceStartingFromSphere(G4VSolid* aSolid,G4int Nstat)
195 {
196  G4double area=1.;
197  G4int i=0;
198  G4int j=0;
199  while (i<Nstat){
200  G4ThreeVector p, direction;
201  area = GenerateAPositionOnASphereBoundary( aSolid,p, direction);
202  G4double dist_to_in = aSolid->DistanceToIn(p,direction);
203  if (dist_to_in<kInfinity/2.) i++;
204  j++;
205  }
206  area=area*double(i)/double(j);
207 
208  return area;
209 }
211 //
212 void G4AdjointPosOnPhysVolGenerator::GenerateAPositionOnASolidBoundary(G4VSolid* aSolid,G4ThreeVector& p, G4ThreeVector& direction)
213 {
214  G4bool find_pos =false;
215  while (!find_pos){
216  if (UseSphere) GenerateAPositionOnASphereBoundary( aSolid,p, direction);
217  else GenerateAPositionOnABoxBoundary( aSolid,p, direction);
218  G4double dist_to_in = aSolid->DistanceToIn(p,direction);
219  if (dist_to_in<kInfinity/2.) {
220  find_pos =true;
221  G4ThreeVector p1=p+ 0.99999*direction*dist_to_in;
222  G4ThreeVector norm =aSolid->SurfaceNormal(p1);
223  p+= 0.999999*direction*dist_to_in;
224  CosThDirComparedToNormal=direction.dot(-norm);
225  }
226  }
227 }
229 //
230 G4double G4AdjointPosOnPhysVolGenerator::GenerateAPositionOnASphereBoundary(G4VSolid* aSolid,G4ThreeVector& p, G4ThreeVector& direction)
231 {
232  G4double minX,maxX,minY,maxY,minZ,maxZ;
233 
234  // values needed for CalculateExtent signature
235 
236  G4VoxelLimits limit; // Unlimited
237  G4AffineTransform origin;
238 
239  // min max extents of pSolid along X,Y,Z
240 
241  aSolid->CalculateExtent(kXAxis,limit,origin,minX,maxX);
242  aSolid->CalculateExtent(kYAxis,limit,origin,minY,maxY);
243  aSolid->CalculateExtent(kZAxis,limit,origin,minZ,maxZ);
244 
245  G4ThreeVector center = G4ThreeVector((minX+maxX)/2.,(minY+maxY)/2.,(minZ+maxZ)/2.);
246 
247  G4double dX=(maxX-minX)/2.;
248  G4double dY=(maxY-minY)/2.;
249  G4double dZ=(maxZ-minZ)/2.;
250  G4double scale=1.01;
251  G4double r=scale*std::sqrt(dX*dX+dY*dY+dZ*dZ);
252 
253  G4double cos_th2 = G4UniformRand();
254  G4double theta = std::acos(std::sqrt(cos_th2));
255  G4double phi=G4UniformRand()*3.1415926*2;
256  direction.setRThetaPhi(1.,theta,phi);
257  direction=-direction;
258  G4double cos_th = (1.-2.*G4UniformRand());
259  theta = std::acos(cos_th);
260  if (G4UniformRand() <0.5) theta=3.1415926-theta;
261  phi=G4UniformRand()*3.1415926*2;
262  p.setRThetaPhi(r,theta,phi);
263  p+=center;
264  direction.rotateY(theta);
265  direction.rotateZ(phi);
266  return 4.*3.1415926*r*r;;
267 }
269 //
270 G4double G4AdjointPosOnPhysVolGenerator::GenerateAPositionOnABoxBoundary(G4VSolid* aSolid,G4ThreeVector& p, G4ThreeVector& direction)
271 {
272 
273  G4double ran_var,px,py,pz,minX,maxX,minY,maxY,minZ,maxZ;
274 
275  // values needed for CalculateExtent signature
276 
277  G4VoxelLimits limit; // Unlimited
278  G4AffineTransform origin;
279 
280  // min max extents of pSolid along X,Y,Z
281 
282  aSolid->CalculateExtent(kXAxis,limit,origin,minX,maxX);
283  aSolid->CalculateExtent(kYAxis,limit,origin,minY,maxY);
284  aSolid->CalculateExtent(kZAxis,limit,origin,minZ,maxZ);
285 
286  G4double scale=.1;
287  minX-=scale*std::abs(minX);
288  minY-=scale*std::abs(minY);
289  minZ-=scale*std::abs(minZ);
290  maxX+=scale*std::abs(maxX);
291  maxY+=scale*std::abs(maxY);
292  maxZ+=scale*std::abs(maxZ);
293 
294  G4double dX=(maxX-minX);
295  G4double dY=(maxY-minY);
296  G4double dZ=(maxZ-minZ);
297 
298  G4double XY_prob=2.*dX*dY;
299  G4double YZ_prob=2.*dY*dZ;
300  G4double ZX_prob=2.*dZ*dX;
301  G4double area=XY_prob+YZ_prob+ZX_prob;
302  XY_prob/=area;
303  YZ_prob/=area;
304  ZX_prob/=area;
305 
306  ran_var=G4UniformRand();
307  G4double cos_th2 = G4UniformRand();
308  G4double sth = std::sqrt(1.-cos_th2);
309  G4double cth = std::sqrt(cos_th2);
310  G4double phi=G4UniformRand()*3.1415926*2;
311  G4double dirX = sth*std::cos(phi);
312  G4double dirY = sth*std::sin(phi);
313  G4double dirZ = cth;
314  if (ran_var <=XY_prob){ //on the XY faces
315  G4double ran_var1=ran_var/XY_prob;
316  G4double ranX=ran_var1;
317  if (ran_var1<=0.5){
318  pz=minZ;
319  direction=G4ThreeVector(dirX,dirY,dirZ);
320  ranX=ran_var1*2.;
321  }
322  else{
323  pz=maxZ;
324  direction=-G4ThreeVector(dirX,dirY,dirZ);
325  ranX=(ran_var1-0.5)*2.;
326  }
327  G4double ranY=G4UniformRand();
328  px=minX+(maxX-minX)*ranX;
329  py=minY+(maxY-minY)*ranY;
330  }
331  else if (ran_var <=(XY_prob+YZ_prob)){ //on the YZ faces
332  G4double ran_var1=(ran_var-XY_prob)/YZ_prob;
333  G4double ranY=ran_var1;
334  if (ran_var1<=0.5){
335  px=minX;
336  direction=G4ThreeVector(dirZ,dirX,dirY);
337  ranY=ran_var1*2.;
338  }
339  else{
340  px=maxX;
341  direction=-G4ThreeVector(dirZ,dirX,dirY);
342  ranY=(ran_var1-0.5)*2.;
343  }
344  G4double ranZ=G4UniformRand();
345  py=minY+(maxY-minY)*ranY;
346  pz=minZ+(maxZ-minZ)*ranZ;
347 
348  }
349  else{ //on the ZX faces
350  G4double ran_var1=(ran_var-XY_prob-YZ_prob)/ZX_prob;
351  G4double ranZ=ran_var1;
352  if (ran_var1<=0.5){
353  py=minY;
354  direction=G4ThreeVector(dirY,dirZ,dirX);
355  ranZ=ran_var1*2.;
356  }
357  else{
358  py=maxY;
359  direction=-G4ThreeVector(dirY,dirZ,dirX);
360  ranZ=(ran_var1-0.5)*2.;
361  }
362  G4double ranX=G4UniformRand();
363  px=minX+(maxX-minX)*ranX;
364  pz=minZ+(maxZ-minZ)*ranZ;
365  }
366 
367  p=G4ThreeVector(px,py,pz);
368  return area;
369 }
371 //
372 void G4AdjointPosOnPhysVolGenerator::GenerateAPositionOnTheExtSurfaceOfThePhysicalVolume(G4ThreeVector& p, G4ThreeVector& direction)
373 {
374  if (!thePhysicalVolume) {
375  G4cout<<"Before generating a source on an external surface of volume you should select a physical volume"<<std::endl;
376  return;
377  };
378  GenerateAPositionOnTheExtSurfaceOfTheSolid(p,direction);
379  p = theTransformationFromPhysVolToWorld.TransformPoint(p);
380  direction = theTransformationFromPhysVolToWorld.TransformAxis(direction);
381 }
383 //
384 void G4AdjointPosOnPhysVolGenerator::GenerateAPositionOnTheExtSurfaceOfThePhysicalVolume(G4ThreeVector& p, G4ThreeVector& direction,
385  G4double& costh_to_normal)
386 {
387  GenerateAPositionOnTheExtSurfaceOfThePhysicalVolume(p, direction);
388  costh_to_normal = CosThDirComparedToNormal;
389 }
391 //
392 void G4AdjointPosOnPhysVolGenerator::ComputeTransformationFromPhysVolToWorld()
393 {
394  G4VPhysicalVolume* daughter =thePhysicalVolume;
395  G4LogicalVolume* mother = thePhysicalVolume->GetMotherLogical();
396  theTransformationFromPhysVolToWorld = G4AffineTransform();
397  G4PhysicalVolumeStore* thePhysVolStore =G4PhysicalVolumeStore::GetInstance();
398  while (mother){
399  theTransformationFromPhysVolToWorld *=
400  G4AffineTransform(daughter->GetFrameRotation(),daughter->GetObjectTranslation());
401  for ( unsigned int i=0; i< thePhysVolStore->size();i++){
402  if ((*thePhysVolStore)[i]->GetLogicalVolume() == mother){
403  daughter = (*thePhysVolStore)[i];
404  mother =daughter->GetMotherLogical();
405  break;
406  };
407  }
408  }
409 }
410 
Generated on Sat May 25 2013 14:33:05 for Geant4 by   doxygen 1.8.4