| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4AdjointPosOnPhysVolGenerator class implem
27 //
28 // Author: L. Desorgher, SpaceIT GmbH - 01.06.
29 // Contract: ESA contract 21435/08/NL/AT
30 // Customer: ESA/ESTEC
31 // -------------------------------------------
32
33 #include "G4AdjointPosOnPhysVolGenerator.hh"
34 #include "G4VSolid.hh"
35 #include "G4VoxelLimits.hh"
36 #include "G4AffineTransform.hh"
37 #include "Randomize.hh"
38 #include "G4VPhysicalVolume.hh"
39 #include "G4PhysicalVolumeStore.hh"
40 #include "G4LogicalVolumeStore.hh"
41
42 G4ThreadLocal G4AdjointPosOnPhysVolGenerator*
43 G4AdjointPosOnPhysVolGenerator::theInstance =
44
45 // -------------------------------------------
46 //
47 G4AdjointPosOnPhysVolGenerator* G4AdjointPosOn
48 {
49 if(theInstance == nullptr)
50 {
51 theInstance = new G4AdjointPosOnPhysVolGen
52 }
53 return theInstance;
54 }
55
56 // -------------------------------------------
57 //
58 G4VPhysicalVolume*
59 G4AdjointPosOnPhysVolGenerator::DefinePhysical
60 {
61 thePhysicalVolume = nullptr;
62 theSolid = nullptr;
63 G4PhysicalVolumeStore* thePhysVolStore = G4P
64 for ( unsigned int i=0; i< thePhysVolStore->
65 {
66 G4String vol_name =(*thePhysVolStore)[i]->
67 if (vol_name.empty())
68 {
69 vol_name = (*thePhysVolStore)[i]->GetLog
70 }
71 if (vol_name == aName)
72 {
73 thePhysicalVolume = (*thePhysVolStore)[i
74 }
75 }
76 if (thePhysicalVolume != nullptr)
77 {
78 theSolid = thePhysicalVolume->GetLogicalVo
79 ComputeTransformationFromPhysVolToWorld();
80 }
81 else
82 {
83 G4cout << "The physical volume with name "
84 << " does not exist!!" << G4endl;
85 G4cout << "Before generating a source on a
86 << "of a volume you should select a
87 << G4endl;
88 }
89 return thePhysicalVolume;
90 }
91
92 // -------------------------------------------
93 //
94 void
95 G4AdjointPosOnPhysVolGenerator::DefinePhysical
96 {
97 thePhysicalVolume = DefinePhysicalVolume(aNa
98 }
99
100 // -------------------------------------------
101 //
102 G4double G4AdjointPosOnPhysVolGenerator::Compu
103 {
104 return ComputeAreaOfExtSurface(theSolid);
105 }
106
107 // -------------------------------------------
108 //
109 G4double G4AdjointPosOnPhysVolGenerator::Compu
110 {
111 return ComputeAreaOfExtSurface(theSolid,NSta
112 }
113
114 // -------------------------------------------
115 //
116 G4double G4AdjointPosOnPhysVolGenerator::Compu
117 {
118 return ComputeAreaOfExtSurface(theSolid,eps)
119 }
120
121 // -------------------------------------------
122 //
123 G4double
124 G4AdjointPosOnPhysVolGenerator::ComputeAreaOfE
125 {
126 return ComputeAreaOfExtSurface(aSolid,1.e-3)
127 }
128
129 // -------------------------------------------
130 //
131 G4double
132 G4AdjointPosOnPhysVolGenerator::ComputeAreaOfE
133
134 {
135 if (ModelOfSurfaceSource == "OnSolid")
136 {
137 if (UseSphere)
138 {
139 return ComputeAreaOfExtSurfaceStartingFr
140 }
141
142 return ComputeAreaOfExtSurfaceStartingFrom
143 }
144
145 G4ThreeVector p, dir;
146 if (ModelOfSurfaceSource == "ExternalSphere"
147 {
148 return GenerateAPositionOnASphereBoundary(
149 }
150
151 return GenerateAPositionOnABoxBoundary(aSoli
152 }
153
154 // -------------------------------------------
155 //
156 G4double
157 G4AdjointPosOnPhysVolGenerator::ComputeAreaOfE
158
159 {
160 auto Nstats = G4int(1./(eps*eps));
161 return ComputeAreaOfExtSurface(aSolid,Nstats
162 }
163
164 // -------------------------------------------
165 //
166 void G4AdjointPosOnPhysVolGenerator::
167 GenerateAPositionOnTheExtSurfaceOfASolid(G4VSo
168 G4Thr
169 {
170 if (ModelOfSurfaceSource == "OnSolid")
171 {
172 GenerateAPositionOnASolidBoundary(aSolid,
173 return;
174 }
175 if (ModelOfSurfaceSource == "ExternalSphere"
176 {
177 GenerateAPositionOnASphereBoundary(aSolid,
178 return;
179 }
180 GenerateAPositionOnABoxBoundary(aSolid, p, d
181 return;
182 }
183
184 // -------------------------------------------
185 //
186 void G4AdjointPosOnPhysVolGenerator::
187 GenerateAPositionOnTheExtSurfaceOfTheSolid(G4T
188 G4T
189 {
190 GenerateAPositionOnTheExtSurfaceOfASolid(the
191 }
192
193 // -------------------------------------------
194 //
195 G4double G4AdjointPosOnPhysVolGenerator::
196 ComputeAreaOfExtSurfaceStartingFromBox(G4VSoli
197 {
198 if ( Nstat <= 0 ) { return 0.; }
199 G4double area=1.;
200 G4int i=0, j=0;
201 while (i<Nstat)
202 {
203 G4ThreeVector p, direction;
204 area = GenerateAPositionOnABoxBoundary( aS
205 G4double dist_to_in = aSolid->DistanceToIn
206 if (dist_to_in<kInfinity/2.) { ++i; }
207 ++j;
208 }
209 area=area*G4double(i)/G4double(j);
210 return area;
211 }
212
213 // -------------------------------------------
214 //
215 G4double G4AdjointPosOnPhysVolGenerator::
216 ComputeAreaOfExtSurfaceStartingFromSphere(G4VS
217 {
218 if ( Nstat <= 0 ) { return 0.; }
219 G4double area=1.;
220 G4int i=0, j=0;
221 while (i<Nstat)
222 {
223 G4ThreeVector p, direction;
224 area = GenerateAPositionOnASphereBoundary(
225 G4double dist_to_in = aSolid->DistanceToIn
226 if (dist_to_in<kInfinity/2.) { ++i; }
227 ++j;
228 }
229 area=area*G4double(i)/G4double(j);
230 return area;
231 }
232
233 // -------------------------------------------
234 //
235 void G4AdjointPosOnPhysVolGenerator::
236 GenerateAPositionOnASolidBoundary(G4VSolid* aS
237 G4ThreeVecto
238 {
239 G4bool find_pos = false;
240 while (!find_pos)
241 {
242 if (UseSphere)
243 {
244 GenerateAPositionOnASphereBoundary( aSol
245 }
246 else
247 {
248 GenerateAPositionOnABoxBoundary( aSolid,
249 }
250 G4double dist_to_in = aSolid->DistanceToIn
251 if (dist_to_in<kInfinity/2.)
252 {
253 find_pos = true;
254 p += 0.999999*direction*dist_to_in;
255 }
256 }
257 }
258
259 // -------------------------------------------
260 //
261 G4double G4AdjointPosOnPhysVolGenerator::
262 GenerateAPositionOnASphereBoundary(G4VSolid* a
263 G4ThreeVect
264 {
265 G4double minX,maxX,minY,maxY,minZ,maxZ;
266
267 // values needed for CalculateExtent signatu
268
269 G4VoxelLimits limit; // Unlim
270 G4AffineTransform origin;
271
272 // min max extents of pSolid along X,Y,Z
273
274 aSolid->CalculateExtent(kXAxis,limit,origin,
275 aSolid->CalculateExtent(kYAxis,limit,origin,
276 aSolid->CalculateExtent(kZAxis,limit,origin,
277
278 G4ThreeVector center = G4ThreeVector((minX+m
279 (minY+m
280 (minZ+m
281 G4double dX=(maxX-minX)/2.;
282 G4double dY=(maxY-minY)/2.;
283 G4double dZ=(maxZ-minZ)/2.;
284 G4double scale=1.01;
285 G4double r=scale*std::sqrt(dX*dX+dY*dY+dZ*dZ
286
287 G4double cos_th2 = G4UniformRand();
288 G4double theta = std::acos(std::sqrt(cos_th2
289 G4double phi=G4UniformRand()*CLHEP::twopi;
290 direction.setRThetaPhi(1.,theta,phi);
291 direction=-direction;
292 G4double cos_th = (1.-2.*G4UniformRand());
293 theta = std::acos(cos_th);
294 if (G4UniformRand() < 0.5) { theta=CLHEP::p
295 phi=G4UniformRand()*CLHEP::twopi;
296 p.setRThetaPhi(r,theta,phi);
297 p+=center;
298 direction.rotateY(theta);
299 direction.rotateZ(phi);
300 return 4.*CLHEP::pi*r*r;;
301 }
302
303 // -------------------------------------------
304 //
305 G4double G4AdjointPosOnPhysVolGenerator::
306 GenerateAPositionOnABoxBoundary(G4VSolid* aSol
307 G4ThreeVector&
308 {
309
310 G4double ran_var,px,py,pz,minX,maxX,minY,max
311
312 // values needed for CalculateExtent signatu
313
314 G4VoxelLimits limit; // Unlim
315 G4AffineTransform origin;
316
317 // min max extents of pSolid along X,Y,Z
318
319 aSolid->CalculateExtent(kXAxis,limit,origin,
320 aSolid->CalculateExtent(kYAxis,limit,origin,
321 aSolid->CalculateExtent(kZAxis,limit,origin,
322
323 G4double scale=.1;
324 minX-=scale*std::abs(minX);
325 minY-=scale*std::abs(minY);
326 minZ-=scale*std::abs(minZ);
327 maxX+=scale*std::abs(maxX);
328 maxY+=scale*std::abs(maxY);
329 maxZ+=scale*std::abs(maxZ);
330
331 G4double dX=(maxX-minX);
332 G4double dY=(maxY-minY);
333 G4double dZ=(maxZ-minZ);
334
335 G4double XY_prob=2.*dX*dY;
336 G4double YZ_prob=2.*dY*dZ;
337 G4double ZX_prob=2.*dZ*dX;
338 G4double area=XY_prob+YZ_prob+ZX_prob;
339 XY_prob/=area;
340 YZ_prob/=area;
341 ZX_prob/=area;
342
343 ran_var=G4UniformRand();
344 G4double cos_th2 = G4UniformRand();
345 G4double sth = std::sqrt(1.-cos_th2);
346 G4double cth = std::sqrt(cos_th2);
347 G4double phi = G4UniformRand()*CLHEP::twopi;
348 G4double dirX = sth*std::cos(phi);
349 G4double dirY = sth*std::sin(phi);
350 G4double dirZ = cth;
351 if (ran_var <=XY_prob) // on the XY faces
352 {
353 G4double ran_var1=ran_var/XY_prob;
354 G4double ranX=ran_var1;
355 if (ran_var1<=0.5)
356 {
357 pz=minZ;
358 direction=G4ThreeVector(dirX,dirY,dirZ);
359 ranX=ran_var1*2.;
360 }
361 else
362 {
363 pz=maxZ;
364 direction=-G4ThreeVector(dirX,dirY,dirZ)
365 ranX=(ran_var1-0.5)*2.;
366 }
367 G4double ranY=G4UniformRand();
368 px=minX+(maxX-minX)*ranX;
369 py=minY+(maxY-minY)*ranY;
370 }
371 else if (ran_var <=(XY_prob+YZ_prob)) // on
372 {
373 G4double ran_var1=(ran_var-XY_prob)/YZ_pro
374 G4double ranY=ran_var1;
375 if (ran_var1<=0.5)
376 {
377 px=minX;
378 direction=G4ThreeVector(dirZ,dirX,dirY);
379 ranY=ran_var1*2.;
380 }
381 else
382 {
383 px=maxX;
384 direction=-G4ThreeVector(dirZ,dirX,dirY)
385 ranY=(ran_var1-0.5)*2.;
386 }
387 G4double ranZ=G4UniformRand();
388 py=minY+(maxY-minY)*ranY;
389 pz=minZ+(maxZ-minZ)*ranZ;
390 }
391 else // on the ZX faces
392 {
393 G4double ran_var1=(ran_var-XY_prob-YZ_prob
394 G4double ranZ=ran_var1;
395 if (ran_var1<=0.5)
396 {
397 py=minY;
398 direction=G4ThreeVector(dirY,dirZ,dirX);
399 ranZ=ran_var1*2.;
400 }
401 else
402 {
403 py=maxY;
404 direction=-G4ThreeVector(dirY,dirZ,dirX)
405 ranZ=(ran_var1-0.5)*2.;
406 }
407 G4double ranX=G4UniformRand();
408 px=minX+(maxX-minX)*ranX;
409 pz=minZ+(maxZ-minZ)*ranZ;
410 }
411
412 p=G4ThreeVector(px,py,pz);
413 return area;
414 }
415
416 // -------------------------------------------
417 //
418 void G4AdjointPosOnPhysVolGenerator::
419 GenerateAPositionOnTheExtSurfaceOfThePhysicalV
420
421 {
422 if (thePhysicalVolume == nullptr)
423 {
424 G4cout << "Before generating a source on a
425 << "of volume you should select a p
426 return;
427 }
428 GenerateAPositionOnTheExtSurfaceOfTheSolid(p
429 p = theTransformationFromPhysVolToWorld.Tran
430 direction = theTransformationFromPhysVolToWo
431 }
432
433 // -------------------------------------------
434 //
435 void G4AdjointPosOnPhysVolGenerator::
436 GenerateAPositionOnTheExtSurfaceOfThePhysicalV
437
438
439 {
440 GenerateAPositionOnTheExtSurfaceOfThePhysica
441 costh_to_normal = CosThDirComparedToNormal;
442 }
443
444 // -------------------------------------------
445 //
446 void G4AdjointPosOnPhysVolGenerator::ComputeTr
447 {
448 G4VPhysicalVolume* daughter = thePhysicalVol
449 G4LogicalVolume* mother = thePhysicalVolume-
450 theTransformationFromPhysVolToWorld = G4Affi
451 G4PhysicalVolumeStore* thePhysVolStore = G4P
452 while (mother != nullptr)
453 {
454 theTransformationFromPhysVolToWorld *=
455 G4AffineTransform(daughter->GetFrameRota
456 daughter->GetObjectTra
457 for ( unsigned int i=0; i<thePhysVolStore-
458 {
459 if ((*thePhysVolStore)[i]->GetLogicalVol
460 {
461 daughter = (*thePhysVolStore)[i];
462 mother = daughter->GetMotherLogical();
463 break;
464 }
465 }
466 }
467 }
468