Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/event/src/G4AdjointPosOnPhysVolGenerator.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /event/src/G4AdjointPosOnPhysVolGenerator.cc (Version 11.3.0) and /event/src/G4AdjointPosOnPhysVolGenerator.cc (Version 11.1.3)


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