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.1.p1)


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