source: trunk/source/event/src/G4AdjointPosOnPhysVolGenerator.cc@ 1340

Last change on this file since 1340 was 1337, checked in by garnier, 15 years ago

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