Geant4 Cross Reference |
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