Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4AdjointPosOnPhysVolGenerator.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 // $Id: G4AdjointPosOnPhysVolGenerator.cc 89370 2015-04-08 07:50:11Z 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 
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 
46 
48 //
50 {
51  if(!theInstance)
52  {
54  }
55  return theInstance;
56 }
57 
59 //
61 {
62 }
63 
65 //
67  : theSolid(0), thePhysicalVolume(0),
68  UseSphere(true), ModelOfSurfaceSource("OnSolid"),
69  AreaOfExtSurfaceOfThePhysicalVolume(0.), CosThDirComparedToNormal(0.)
70 {
71 }
72 
74 //
76 {
78  theSolid =0;
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){
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 //
104 {
106 }
108 //
110 {
112 }
114 //
116 {
117  return ComputeAreaOfExtSurface(theSolid,NStats);
118 }
120 //
122 {
123  return ComputeAreaOfExtSurface(theSolid,eps);
124 }
126 //
128 {
129  return ComputeAreaOfExtSurface(aSolid,1.e-3);
130 }
132 //
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 {
145  if (ModelOfSurfaceSource == "ExternalSphere" ) return GenerateAPositionOnASphereBoundary(aSolid, p,dir);
146  return GenerateAPositionOnABoxBoundary(aSolid, p,dir);
147  }
148 }
150 //
152 {
153  G4int Nstats = G4int(1./(eps*eps));
154  return ComputeAreaOfExtSurface(aSolid,Nstats);
155 }
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 }
172 {
174 }
176 //
178 {
179  if ( Nstat <= 0 ) return 0.;
180  G4double area=1.;
181  G4int i=0;
182  G4int j=0;
183  while (i<Nstat){
184  G4ThreeVector p, direction;
185  area = GenerateAPositionOnABoxBoundary( aSolid,p, direction);
186  G4double dist_to_in = aSolid->DistanceToIn(p,direction);
187  if (dist_to_in<kInfinity/2.) i++;
188  j++;
189  }
190  area=area*double(i)/double(j);
191  return area;
192 }
194 //
196 {
197  if ( Nstat <= 0 ) return 0.;
198  G4double area=1.;
199  G4int i=0;
200  G4int j=0;
201  while (i<Nstat){
202  G4ThreeVector p, direction;
203  area = GenerateAPositionOnASphereBoundary( aSolid,p, direction);
204  G4double dist_to_in = aSolid->DistanceToIn(p,direction);
205  if (dist_to_in<kInfinity/2.) i++;
206  j++;
207  }
208  area=area*double(i)/double(j);
209  return area;
210 }
212 //
214 {
215  G4bool find_pos =false;
216  while (!find_pos){
217  if (UseSphere) GenerateAPositionOnASphereBoundary( aSolid,p, direction);
218  else GenerateAPositionOnABoxBoundary( aSolid,p, direction);
219  G4double dist_to_in = aSolid->DistanceToIn(p,direction);
220  if (dist_to_in<kInfinity/2.) {
221  find_pos =true;
222  p+= 0.999999*direction*dist_to_in;
223  }
224  }
225 }
227 //
229 {
230  G4double minX,maxX,minY,maxY,minZ,maxZ;
231 
232  // values needed for CalculateExtent signature
233 
234  G4VoxelLimits limit; // Unlimited
235  G4AffineTransform origin;
236 
237  // min max extents of pSolid along X,Y,Z
238 
239  aSolid->CalculateExtent(kXAxis,limit,origin,minX,maxX);
240  aSolid->CalculateExtent(kYAxis,limit,origin,minY,maxY);
241  aSolid->CalculateExtent(kZAxis,limit,origin,minZ,maxZ);
242 
243  G4ThreeVector center = G4ThreeVector((minX+maxX)/2.,(minY+maxY)/2.,(minZ+maxZ)/2.);
244 
245  G4double dX=(maxX-minX)/2.;
246  G4double dY=(maxY-minY)/2.;
247  G4double dZ=(maxZ-minZ)/2.;
248  G4double scale=1.01;
249  G4double r=scale*std::sqrt(dX*dX+dY*dY+dZ*dZ);
250 
251  G4double cos_th2 = G4UniformRand();
252  G4double theta = std::acos(std::sqrt(cos_th2));
253  G4double phi=G4UniformRand()*3.1415926*2;
254  direction.setRThetaPhi(1.,theta,phi);
255  direction=-direction;
256  G4double cos_th = (1.-2.*G4UniformRand());
257  theta = std::acos(cos_th);
258  if (G4UniformRand() <0.5) theta=3.1415926-theta;
259  phi=G4UniformRand()*3.1415926*2;
260  p.setRThetaPhi(r,theta,phi);
261  p+=center;
262  direction.rotateY(theta);
263  direction.rotateZ(phi);
264  return 4.*3.1415926*r*r;;
265 }
267 //
269 {
270 
271  G4double ran_var,px,py,pz,minX,maxX,minY,maxY,minZ,maxZ;
272 
273  // values needed for CalculateExtent signature
274 
275  G4VoxelLimits limit; // Unlimited
276  G4AffineTransform origin;
277 
278  // min max extents of pSolid along X,Y,Z
279 
280  aSolid->CalculateExtent(kXAxis,limit,origin,minX,maxX);
281  aSolid->CalculateExtent(kYAxis,limit,origin,minY,maxY);
282  aSolid->CalculateExtent(kZAxis,limit,origin,minZ,maxZ);
283 
284  G4double scale=.1;
285  minX-=scale*std::abs(minX);
286  minY-=scale*std::abs(minY);
287  minZ-=scale*std::abs(minZ);
288  maxX+=scale*std::abs(maxX);
289  maxY+=scale*std::abs(maxY);
290  maxZ+=scale*std::abs(maxZ);
291 
292  G4double dX=(maxX-minX);
293  G4double dY=(maxY-minY);
294  G4double dZ=(maxZ-minZ);
295 
296  G4double XY_prob=2.*dX*dY;
297  G4double YZ_prob=2.*dY*dZ;
298  G4double ZX_prob=2.*dZ*dX;
299  G4double area=XY_prob+YZ_prob+ZX_prob;
300  XY_prob/=area;
301  YZ_prob/=area;
302  ZX_prob/=area;
303 
304  ran_var=G4UniformRand();
305  G4double cos_th2 = G4UniformRand();
306  G4double sth = std::sqrt(1.-cos_th2);
307  G4double cth = std::sqrt(cos_th2);
308  G4double phi=G4UniformRand()*3.1415926*2;
309  G4double dirX = sth*std::cos(phi);
310  G4double dirY = sth*std::sin(phi);
311  G4double dirZ = cth;
312  if (ran_var <=XY_prob){ //on the XY faces
313  G4double ran_var1=ran_var/XY_prob;
314  G4double ranX=ran_var1;
315  if (ran_var1<=0.5){
316  pz=minZ;
317  direction=G4ThreeVector(dirX,dirY,dirZ);
318  ranX=ran_var1*2.;
319  }
320  else{
321  pz=maxZ;
322  direction=-G4ThreeVector(dirX,dirY,dirZ);
323  ranX=(ran_var1-0.5)*2.;
324  }
325  G4double ranY=G4UniformRand();
326  px=minX+(maxX-minX)*ranX;
327  py=minY+(maxY-minY)*ranY;
328  }
329  else if (ran_var <=(XY_prob+YZ_prob)){ //on the YZ faces
330  G4double ran_var1=(ran_var-XY_prob)/YZ_prob;
331  G4double ranY=ran_var1;
332  if (ran_var1<=0.5){
333  px=minX;
334  direction=G4ThreeVector(dirZ,dirX,dirY);
335  ranY=ran_var1*2.;
336  }
337  else{
338  px=maxX;
339  direction=-G4ThreeVector(dirZ,dirX,dirY);
340  ranY=(ran_var1-0.5)*2.;
341  }
342  G4double ranZ=G4UniformRand();
343  py=minY+(maxY-minY)*ranY;
344  pz=minZ+(maxZ-minZ)*ranZ;
345 
346  }
347  else{ //on the ZX faces
348  G4double ran_var1=(ran_var-XY_prob-YZ_prob)/ZX_prob;
349  G4double ranZ=ran_var1;
350  if (ran_var1<=0.5){
351  py=minY;
352  direction=G4ThreeVector(dirY,dirZ,dirX);
353  ranZ=ran_var1*2.;
354  }
355  else{
356  py=maxY;
357  direction=-G4ThreeVector(dirY,dirZ,dirX);
358  ranZ=(ran_var1-0.5)*2.;
359  }
360  G4double ranX=G4UniformRand();
361  px=minX+(maxX-minX)*ranX;
362  pz=minZ+(maxZ-minZ)*ranZ;
363  }
364 
365  p=G4ThreeVector(px,py,pz);
366  return area;
367 }
369 //
371 {
372  if (!thePhysicalVolume) {
373  G4cout<<"Before generating a source on an external surface of volume you should select a physical volume"<<std::endl;
374  return;
375  };
378  direction = theTransformationFromPhysVolToWorld.TransformAxis(direction);
379 }
381 //
383  G4double& costh_to_normal)
384 {
386  costh_to_normal = CosThDirComparedToNormal;
387 }
389 //
391 {
396  while (mother){
399  for ( unsigned int i=0; i< thePhysVolStore->size();i++){
400  if ((*thePhysVolStore)[i]->GetLogicalVolume() == mother){
401  daughter = (*thePhysVolStore)[i];
402  mother =daughter->GetMotherLogical();
403  break;
404  };
405  }
406  }
407 }
408 
const G4RotationMatrix * GetFrameRotation() const
G4VPhysicalVolume * DefinePhysicalVolume(const G4String &aName)
G4ThreeVector TransformAxis(const G4ThreeVector &axis) const
CLHEP::Hep3Vector G4ThreeVector
double minY
Definition: plot_hist.C:9
G4double ComputeAreaOfExtSurfaceStartingFromSphere(G4VSolid *aSolid, G4int NStat)
G4LogicalVolume * GetLogicalVolume() const
static const G4double kInfinity
Definition: geomdefs.hh:42
virtual G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const =0
const char * p
Definition: xmltok.h:285
double maxY
Definition: plot_hist.C:10
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const =0
#define G4ThreadLocal
Definition: tls.hh:69
static const G4int maxZ
G4double ComputeAreaOfExtSurfaceStartingFromBox(G4VSolid *aSolid, G4int NStat)
void GenerateAPositionOnTheExtSurfaceOfThePhysicalVolume(G4ThreeVector &p, G4ThreeVector &direction)
G4double GenerateAPositionOnASphereBoundary(G4VSolid *aSolid, G4ThreeVector &p, G4ThreeVector &direction)
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
Double_t scale
G4LogicalVolume * GetMotherLogical() const
static G4ThreadLocal G4AdjointPosOnPhysVolGenerator * theInstance
#define G4UniformRand()
Definition: Randomize.hh:53
void DefinePhysicalVolume1(const G4String &aName)
G4double GenerateAPositionOnABoxBoundary(G4VSolid *aSolid, G4ThreeVector &p, G4ThreeVector &direction)
void GenerateAPositionOnASolidBoundary(G4VSolid *aSolid, G4ThreeVector &p, G4ThreeVector &direction)
void GenerateAPositionOnTheExtSurfaceOfASolid(G4VSolid *aSolid, G4ThreeVector &p, G4ThreeVector &direction)
Hep3Vector & rotateZ(double)
Definition: ThreeVector.cc:110
static G4PhysicalVolumeStore * GetInstance()
int G4int
Definition: G4Types.hh:78
void setRThetaPhi(double r, double theta, double phi)
TDirectory * dir
static G4AdjointPosOnPhysVolGenerator * GetInstance()
G4GLOB_DLL std::ostream G4cout
G4VSolid * GetSolid() const
void GenerateAPositionOnTheExtSurfaceOfTheSolid(G4ThreeVector &p, G4ThreeVector &direction)
static const G4double eps
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4ThreeVector GetObjectTranslation() const
Hep3Vector & rotateY(double)
Definition: ThreeVector.cc:100