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