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