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

Last change on this file since 1347 was 1337, checked in by garnier, 14 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.