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