00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00029
00030
00031
00032
00033
00035
00036 #include "G4AdjointPosOnPhysVolGenerator.hh"
00037 #include "G4VSolid.hh"
00038 #include "G4VoxelLimits.hh"
00039 #include "G4AffineTransform.hh"
00040 #include "Randomize.hh"
00041 #include "G4VPhysicalVolume.hh"
00042 #include "G4PhysicalVolumeStore.hh"
00043 #include "G4LogicalVolumeStore.hh"
00044
00045 G4AdjointPosOnPhysVolGenerator* G4AdjointPosOnPhysVolGenerator::theInstance = 0;
00046
00048
00049 G4AdjointPosOnPhysVolGenerator* G4AdjointPosOnPhysVolGenerator::GetInstance()
00050 {
00051 if(theInstance == 0) {
00052 static G4AdjointPosOnPhysVolGenerator manager;
00053 theInstance = &manager;
00054 }
00055 return theInstance;
00056 }
00057
00059
00060 G4AdjointPosOnPhysVolGenerator::~G4AdjointPosOnPhysVolGenerator()
00061 {
00062 }
00063
00065
00066 G4AdjointPosOnPhysVolGenerator::G4AdjointPosOnPhysVolGenerator()
00067 : theSolid(0), thePhysicalVolume(0),
00068 UseSphere(true), ModelOfSurfaceSource("OnSolid"),
00069 AreaOfExtSurfaceOfThePhysicalVolume(0.), CosThDirComparedToNormal(0.)
00070 {
00071 }
00072
00074
00075 G4VPhysicalVolume* G4AdjointPosOnPhysVolGenerator::DefinePhysicalVolume(const G4String& aName)
00076 {
00077 thePhysicalVolume = 0;
00078 theSolid =0;
00079 G4PhysicalVolumeStore* thePhysVolStore =G4PhysicalVolumeStore::GetInstance();
00080 for ( unsigned int i=0; i< thePhysVolStore->size();i++){
00081 G4String vol_name =(*thePhysVolStore)[i]->GetName();
00082 if (vol_name == ""){
00083 vol_name = (*thePhysVolStore)[i]->GetLogicalVolume()->GetName();
00084 }
00085 if (vol_name == aName){
00086 thePhysicalVolume = (*thePhysVolStore)[i];
00087 }
00088 }
00089 if (thePhysicalVolume){
00090 theSolid = thePhysicalVolume->GetLogicalVolume()->GetSolid();
00091 ComputeTransformationFromPhysVolToWorld();
00092
00093
00094 }
00095 else {
00096 G4cout<<"The physical volume with name "<<aName<<" does not exist!!"<<std::endl;
00097 G4cout<<"Before generating a source on an external surface of a volume you should select another physical volume"<<std::endl;
00098 }
00099 return thePhysicalVolume;
00100 }
00102
00103 void G4AdjointPosOnPhysVolGenerator::DefinePhysicalVolume1(const G4String& aName)
00104 {
00105 thePhysicalVolume = DefinePhysicalVolume(aName);
00106 }
00108
00109 G4double G4AdjointPosOnPhysVolGenerator::ComputeAreaOfExtSurface()
00110 {
00111 return ComputeAreaOfExtSurface(theSolid);
00112 }
00114
00115 G4double G4AdjointPosOnPhysVolGenerator::ComputeAreaOfExtSurface(G4int NStats)
00116 {
00117 return ComputeAreaOfExtSurface(theSolid,NStats);
00118 }
00120
00121 G4double G4AdjointPosOnPhysVolGenerator::ComputeAreaOfExtSurface(G4double eps)
00122 {
00123 return ComputeAreaOfExtSurface(theSolid,eps);
00124 }
00126
00127 G4double G4AdjointPosOnPhysVolGenerator::ComputeAreaOfExtSurface(G4VSolid* aSolid)
00128 {
00129 return ComputeAreaOfExtSurface(aSolid,1.e-3);
00130 }
00132
00133 G4double G4AdjointPosOnPhysVolGenerator::ComputeAreaOfExtSurface(G4VSolid* aSolid,G4int NStats)
00134 {
00135 if (ModelOfSurfaceSource == "OnSolid" ){
00136 if (UseSphere){
00137 return ComputeAreaOfExtSurfaceStartingFromSphere(aSolid,NStats);
00138 }
00139 else {
00140 return ComputeAreaOfExtSurfaceStartingFromBox(aSolid,NStats);
00141 }
00142 }
00143 else {
00144 G4ThreeVector p,dir;
00145 if (ModelOfSurfaceSource == "ExternalSphere" ) return GenerateAPositionOnASphereBoundary(aSolid, p,dir);
00146 return GenerateAPositionOnABoxBoundary(aSolid, p,dir);
00147 }
00148 }
00150
00151 G4double G4AdjointPosOnPhysVolGenerator::ComputeAreaOfExtSurface(G4VSolid* aSolid,G4double eps)
00152 {
00153 G4int Nstats = G4int(1./(eps*eps));
00154 return ComputeAreaOfExtSurface(aSolid,Nstats);
00155 }
00157 void G4AdjointPosOnPhysVolGenerator::GenerateAPositionOnTheExtSurfaceOfASolid(G4VSolid* aSolid,G4ThreeVector& p, G4ThreeVector& direction)
00158 {
00159 if (ModelOfSurfaceSource == "OnSolid" ){
00160 GenerateAPositionOnASolidBoundary(aSolid, p,direction);
00161 return;
00162 }
00163 if (ModelOfSurfaceSource == "ExternalSphere" ) {
00164 GenerateAPositionOnASphereBoundary(aSolid, p, direction);
00165 return;
00166 }
00167 GenerateAPositionOnABoxBoundary(aSolid, p, direction);
00168 return;
00169 }
00171 void G4AdjointPosOnPhysVolGenerator::GenerateAPositionOnTheExtSurfaceOfTheSolid(G4ThreeVector& p, G4ThreeVector& direction)
00172 {
00173 GenerateAPositionOnTheExtSurfaceOfASolid(theSolid,p,direction);
00174 }
00176
00177 G4double G4AdjointPosOnPhysVolGenerator::ComputeAreaOfExtSurfaceStartingFromBox(G4VSolid* aSolid,G4int Nstat)
00178 {
00179 G4double area=1.;
00180 G4int i=0;
00181 G4int j=0;
00182 while (i<Nstat){
00183 G4ThreeVector p, direction;
00184 area = GenerateAPositionOnABoxBoundary( aSolid,p, direction);
00185 G4double dist_to_in = aSolid->DistanceToIn(p,direction);
00186 if (dist_to_in<kInfinity/2.) i++;
00187 j++;
00188 }
00189 area=area*double(i)/double(j);
00190 return area;
00191 }
00193
00194 G4double G4AdjointPosOnPhysVolGenerator::ComputeAreaOfExtSurfaceStartingFromSphere(G4VSolid* aSolid,G4int Nstat)
00195 {
00196 G4double area=1.;
00197 G4int i=0;
00198 G4int j=0;
00199 while (i<Nstat){
00200 G4ThreeVector p, direction;
00201 area = GenerateAPositionOnASphereBoundary( aSolid,p, direction);
00202 G4double dist_to_in = aSolid->DistanceToIn(p,direction);
00203 if (dist_to_in<kInfinity/2.) i++;
00204 j++;
00205 }
00206 area=area*double(i)/double(j);
00207
00208 return area;
00209 }
00211
00212 void G4AdjointPosOnPhysVolGenerator::GenerateAPositionOnASolidBoundary(G4VSolid* aSolid,G4ThreeVector& p, G4ThreeVector& direction)
00213 {
00214 G4bool find_pos =false;
00215 while (!find_pos){
00216 if (UseSphere) GenerateAPositionOnASphereBoundary( aSolid,p, direction);
00217 else GenerateAPositionOnABoxBoundary( aSolid,p, direction);
00218 G4double dist_to_in = aSolid->DistanceToIn(p,direction);
00219 if (dist_to_in<kInfinity/2.) {
00220 find_pos =true;
00221 G4ThreeVector p1=p+ 0.99999*direction*dist_to_in;
00222 G4ThreeVector norm =aSolid->SurfaceNormal(p1);
00223 p+= 0.999999*direction*dist_to_in;
00224 CosThDirComparedToNormal=direction.dot(-norm);
00225 }
00226 }
00227 }
00229
00230 G4double G4AdjointPosOnPhysVolGenerator::GenerateAPositionOnASphereBoundary(G4VSolid* aSolid,G4ThreeVector& p, G4ThreeVector& direction)
00231 {
00232 G4double minX,maxX,minY,maxY,minZ,maxZ;
00233
00234
00235
00236 G4VoxelLimits limit;
00237 G4AffineTransform origin;
00238
00239
00240
00241 aSolid->CalculateExtent(kXAxis,limit,origin,minX,maxX);
00242 aSolid->CalculateExtent(kYAxis,limit,origin,minY,maxY);
00243 aSolid->CalculateExtent(kZAxis,limit,origin,minZ,maxZ);
00244
00245 G4ThreeVector center = G4ThreeVector((minX+maxX)/2.,(minY+maxY)/2.,(minZ+maxZ)/2.);
00246
00247 G4double dX=(maxX-minX)/2.;
00248 G4double dY=(maxY-minY)/2.;
00249 G4double dZ=(maxZ-minZ)/2.;
00250 G4double scale=1.01;
00251 G4double r=scale*std::sqrt(dX*dX+dY*dY+dZ*dZ);
00252
00253 G4double cos_th2 = G4UniformRand();
00254 G4double theta = std::acos(std::sqrt(cos_th2));
00255 G4double phi=G4UniformRand()*3.1415926*2;
00256 direction.setRThetaPhi(1.,theta,phi);
00257 direction=-direction;
00258 G4double cos_th = (1.-2.*G4UniformRand());
00259 theta = std::acos(cos_th);
00260 if (G4UniformRand() <0.5) theta=3.1415926-theta;
00261 phi=G4UniformRand()*3.1415926*2;
00262 p.setRThetaPhi(r,theta,phi);
00263 p+=center;
00264 direction.rotateY(theta);
00265 direction.rotateZ(phi);
00266 return 4.*3.1415926*r*r;;
00267 }
00269
00270 G4double G4AdjointPosOnPhysVolGenerator::GenerateAPositionOnABoxBoundary(G4VSolid* aSolid,G4ThreeVector& p, G4ThreeVector& direction)
00271 {
00272
00273 G4double ran_var,px,py,pz,minX,maxX,minY,maxY,minZ,maxZ;
00274
00275
00276
00277 G4VoxelLimits limit;
00278 G4AffineTransform origin;
00279
00280
00281
00282 aSolid->CalculateExtent(kXAxis,limit,origin,minX,maxX);
00283 aSolid->CalculateExtent(kYAxis,limit,origin,minY,maxY);
00284 aSolid->CalculateExtent(kZAxis,limit,origin,minZ,maxZ);
00285
00286 G4double scale=.1;
00287 minX-=scale*std::abs(minX);
00288 minY-=scale*std::abs(minY);
00289 minZ-=scale*std::abs(minZ);
00290 maxX+=scale*std::abs(maxX);
00291 maxY+=scale*std::abs(maxY);
00292 maxZ+=scale*std::abs(maxZ);
00293
00294 G4double dX=(maxX-minX);
00295 G4double dY=(maxY-minY);
00296 G4double dZ=(maxZ-minZ);
00297
00298 G4double XY_prob=2.*dX*dY;
00299 G4double YZ_prob=2.*dY*dZ;
00300 G4double ZX_prob=2.*dZ*dX;
00301 G4double area=XY_prob+YZ_prob+ZX_prob;
00302 XY_prob/=area;
00303 YZ_prob/=area;
00304 ZX_prob/=area;
00305
00306 ran_var=G4UniformRand();
00307 G4double cos_th2 = G4UniformRand();
00308 G4double sth = std::sqrt(1.-cos_th2);
00309 G4double cth = std::sqrt(cos_th2);
00310 G4double phi=G4UniformRand()*3.1415926*2;
00311 G4double dirX = sth*std::cos(phi);
00312 G4double dirY = sth*std::sin(phi);
00313 G4double dirZ = cth;
00314 if (ran_var <=XY_prob){
00315 G4double ran_var1=ran_var/XY_prob;
00316 G4double ranX=ran_var1;
00317 if (ran_var1<=0.5){
00318 pz=minZ;
00319 direction=G4ThreeVector(dirX,dirY,dirZ);
00320 ranX=ran_var1*2.;
00321 }
00322 else{
00323 pz=maxZ;
00324 direction=-G4ThreeVector(dirX,dirY,dirZ);
00325 ranX=(ran_var1-0.5)*2.;
00326 }
00327 G4double ranY=G4UniformRand();
00328 px=minX+(maxX-minX)*ranX;
00329 py=minY+(maxY-minY)*ranY;
00330 }
00331 else if (ran_var <=(XY_prob+YZ_prob)){
00332 G4double ran_var1=(ran_var-XY_prob)/YZ_prob;
00333 G4double ranY=ran_var1;
00334 if (ran_var1<=0.5){
00335 px=minX;
00336 direction=G4ThreeVector(dirZ,dirX,dirY);
00337 ranY=ran_var1*2.;
00338 }
00339 else{
00340 px=maxX;
00341 direction=-G4ThreeVector(dirZ,dirX,dirY);
00342 ranY=(ran_var1-0.5)*2.;
00343 }
00344 G4double ranZ=G4UniformRand();
00345 py=minY+(maxY-minY)*ranY;
00346 pz=minZ+(maxZ-minZ)*ranZ;
00347
00348 }
00349 else{
00350 G4double ran_var1=(ran_var-XY_prob-YZ_prob)/ZX_prob;
00351 G4double ranZ=ran_var1;
00352 if (ran_var1<=0.5){
00353 py=minY;
00354 direction=G4ThreeVector(dirY,dirZ,dirX);
00355 ranZ=ran_var1*2.;
00356 }
00357 else{
00358 py=maxY;
00359 direction=-G4ThreeVector(dirY,dirZ,dirX);
00360 ranZ=(ran_var1-0.5)*2.;
00361 }
00362 G4double ranX=G4UniformRand();
00363 px=minX+(maxX-minX)*ranX;
00364 pz=minZ+(maxZ-minZ)*ranZ;
00365 }
00366
00367 p=G4ThreeVector(px,py,pz);
00368 return area;
00369 }
00371
00372 void G4AdjointPosOnPhysVolGenerator::GenerateAPositionOnTheExtSurfaceOfThePhysicalVolume(G4ThreeVector& p, G4ThreeVector& direction)
00373 {
00374 if (!thePhysicalVolume) {
00375 G4cout<<"Before generating a source on an external surface of volume you should select a physical volume"<<std::endl;
00376 return;
00377 };
00378 GenerateAPositionOnTheExtSurfaceOfTheSolid(p,direction);
00379 p = theTransformationFromPhysVolToWorld.TransformPoint(p);
00380 direction = theTransformationFromPhysVolToWorld.TransformAxis(direction);
00381 }
00383
00384 void G4AdjointPosOnPhysVolGenerator::GenerateAPositionOnTheExtSurfaceOfThePhysicalVolume(G4ThreeVector& p, G4ThreeVector& direction,
00385 G4double& costh_to_normal)
00386 {
00387 GenerateAPositionOnTheExtSurfaceOfThePhysicalVolume(p, direction);
00388 costh_to_normal = CosThDirComparedToNormal;
00389 }
00391
00392 void G4AdjointPosOnPhysVolGenerator::ComputeTransformationFromPhysVolToWorld()
00393 {
00394 G4VPhysicalVolume* daughter =thePhysicalVolume;
00395 G4LogicalVolume* mother = thePhysicalVolume->GetMotherLogical();
00396 theTransformationFromPhysVolToWorld = G4AffineTransform();
00397 G4PhysicalVolumeStore* thePhysVolStore =G4PhysicalVolumeStore::GetInstance();
00398 while (mother){
00399 theTransformationFromPhysVolToWorld *=
00400 G4AffineTransform(daughter->GetFrameRotation(),daughter->GetObjectTranslation());
00401 for ( unsigned int i=0; i< thePhysVolStore->size();i++){
00402 if ((*thePhysVolStore)[i]->GetLogicalVolume() == mother){
00403 daughter = (*thePhysVolStore)[i];
00404 mother =daughter->GetMotherLogical();
00405 break;
00406 };
00407 }
00408 }
00409 }
00410
