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