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 10.6.p3)


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