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 "G4AdjointCrossSurfChecker.hh"
00037 #include "G4PhysicalConstants.hh"
00038 #include "G4SystemOfUnits.hh"
00039 #include "G4Step.hh"
00040 #include "G4StepPoint.hh"
00041 #include "G4PhysicalVolumeStore.hh"
00042 #include "G4VSolid.hh"
00043 #include "G4AffineTransform.hh"
00044
00046
00047 G4AdjointCrossSurfChecker* G4AdjointCrossSurfChecker::instance = 0;
00048
00050
00051 G4AdjointCrossSurfChecker::G4AdjointCrossSurfChecker()
00052 {;
00053 }
00055
00056 G4AdjointCrossSurfChecker::~G4AdjointCrossSurfChecker()
00057 {
00058 delete instance;
00059 }
00061
00062 G4AdjointCrossSurfChecker* G4AdjointCrossSurfChecker::GetInstance()
00063 {
00064 if (!instance) instance = new G4AdjointCrossSurfChecker();
00065 return instance;
00066 }
00068
00069 G4bool G4AdjointCrossSurfChecker::CrossingASphere(const G4Step* aStep,G4double sphere_radius, G4ThreeVector sphere_center,G4ThreeVector& crossing_pos, G4double& cos_th , G4bool& GoingIn)
00070 {
00071 G4ThreeVector pos1= aStep->GetPreStepPoint()->GetPosition() - sphere_center;
00072 G4ThreeVector pos2= aStep->GetPostStepPoint()->GetPosition() - sphere_center;
00073 G4double r1= pos1.mag();
00074 G4double r2= pos2.mag();
00075 G4bool did_cross =false;
00076
00077 if (r1<=sphere_radius && r2>sphere_radius){
00078 did_cross=true;
00079 GoingIn=false;
00080 }
00081 else if (r2<=sphere_radius && r1>sphere_radius){
00082 did_cross=true;
00083 GoingIn=true;
00084 }
00085
00086 if (did_cross) {
00087
00088 G4ThreeVector dr=pos2-pos1;
00089 G4double r12 = r1*r1;
00090 G4double rdr = dr.mag();
00091 G4double a,b,c,d;
00092 a = rdr*rdr;
00093 b = 2.*pos1.dot(dr);
00094 c = r12-sphere_radius*sphere_radius;
00095 d=std::sqrt(b*b-4.*a*c);
00096 G4double l= (-b+d)/2./a;
00097 if (l > 1.) l=(-b-d)/2./a;
00098 crossing_pos=pos1+l*dr;
00099 cos_th = std::abs(dr.cosTheta(crossing_pos));
00100
00101
00102 }
00103 return did_cross;
00104 }
00106
00107 G4bool G4AdjointCrossSurfChecker::GoingInOrOutOfaVolume(const G4Step* aStep,const G4String& volume_name, G4double& , G4bool& GoingIn)
00108 {
00109 G4bool step_at_boundary = (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
00110 G4bool did_cross =false;
00111 if (step_at_boundary){
00112 const G4VTouchable* postStepTouchable = aStep->GetPostStepPoint()->GetTouchable();
00113 const G4VTouchable* preStepTouchable = aStep->GetPreStepPoint()->GetTouchable();
00114 if (preStepTouchable && postStepTouchable && postStepTouchable->GetVolume() && preStepTouchable->GetVolume()){
00115 G4String post_vol_name = postStepTouchable->GetVolume()->GetName();
00116 G4String pre_vol_name = preStepTouchable->GetVolume()->GetName();
00117
00118 if (post_vol_name == volume_name ){
00119 GoingIn=true;
00120 did_cross=true;
00121 }
00122 else if (pre_vol_name == volume_name){
00123 GoingIn=false;
00124 did_cross=true;
00125
00126 }
00127
00128 }
00129 }
00130 return did_cross;
00131
00132 }
00134
00135 G4bool G4AdjointCrossSurfChecker::GoingInOrOutOfaVolumeByExtSurface(const G4Step* aStep,const G4String& volume_name, const G4String& mother_logical_vol_name, G4double& , G4bool& GoingIn)
00136 {
00137 G4bool step_at_boundary = (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
00138 G4bool did_cross =false;
00139 if (step_at_boundary){
00140 const G4VTouchable* postStepTouchable = aStep->GetPostStepPoint()->GetTouchable();
00141 const G4VTouchable* preStepTouchable = aStep->GetPreStepPoint()->GetTouchable();
00142 if (preStepTouchable && postStepTouchable && postStepTouchable->GetVolume() && preStepTouchable->GetVolume()){
00143 G4String post_vol_name = postStepTouchable->GetVolume()->GetName();
00144 G4String post_log_vol_name = postStepTouchable->GetVolume()->GetLogicalVolume()->GetName();
00145 G4String pre_vol_name = preStepTouchable->GetVolume()->GetName();
00146 G4String pre_log_vol_name = preStepTouchable->GetVolume()->GetLogicalVolume()->GetName();
00147 if (post_vol_name == volume_name && pre_log_vol_name == mother_logical_vol_name){
00148 GoingIn=true;
00149 did_cross=true;
00150 }
00151 else if (pre_vol_name == volume_name && post_log_vol_name == mother_logical_vol_name ){
00152 GoingIn=false;
00153 did_cross=true;
00154
00155 }
00156
00157 }
00158 }
00159 return did_cross;
00160
00161 }
00163
00164 G4bool G4AdjointCrossSurfChecker::CrossingAGivenRegisteredSurface(const G4Step* aStep,const G4String& surface_name,G4ThreeVector& crossing_pos, G4double& cos_to_surface, G4bool& GoingIn)
00165 {
00166 G4int ind = FindRegisteredSurface(surface_name);
00167 G4bool did_cross = false;
00168 if (ind >=0){
00169 did_cross = CrossingAGivenRegisteredSurface(aStep, ind, crossing_pos,cos_to_surface, GoingIn);
00170 }
00171 return did_cross;
00172 }
00174
00175 G4bool G4AdjointCrossSurfChecker::CrossingAGivenRegisteredSurface(const G4Step* aStep, int ind,G4ThreeVector& crossing_pos, G4double& cos_to_surface, G4bool& GoingIn)
00176 {
00177 G4String surf_type = ListOfSurfaceType[ind];
00178 G4double radius = ListOfSphereRadius[ind];
00179 G4ThreeVector center = ListOfSphereCenter[ind];
00180 G4String vol1 = ListOfVol1Name[ind];
00181 G4String vol2 = ListOfVol2Name[ind];
00182
00183 G4bool did_cross = false;
00184 if (surf_type == "Sphere"){
00185 did_cross = CrossingASphere(aStep, radius, center,crossing_pos, cos_to_surface, GoingIn);
00186 }
00187 else if (surf_type == "ExternalSurfaceOfAVolume"){
00188
00189 did_cross = GoingInOrOutOfaVolumeByExtSurface(aStep, vol1, vol2, cos_to_surface, GoingIn);
00190 crossing_pos= aStep->GetPostStepPoint()->GetPosition();
00191
00192 }
00193 else if (surf_type == "BoundaryBetweenTwoVolumes"){
00194 did_cross = CrossingAnInterfaceBetweenTwoVolumes(aStep, vol1, vol2,crossing_pos, cos_to_surface, GoingIn);
00195 }
00196 return did_cross;
00197
00198
00199 }
00201
00202
00203 G4bool G4AdjointCrossSurfChecker::CrossingOneOfTheRegisteredSurface(const G4Step* aStep,G4String& surface_name,G4ThreeVector& crossing_pos, G4double& cos_to_surface, G4bool& GoingIn)
00204 {
00205 for (size_t i=0;i <ListOfSurfaceName.size();i++){
00206 if (CrossingAGivenRegisteredSurface(aStep, int(i),crossing_pos, cos_to_surface, GoingIn)){
00207 surface_name = ListOfSurfaceName[i];
00208 return true;
00209 }
00210 }
00211 return false;
00212 }
00214
00215 G4bool G4AdjointCrossSurfChecker::CrossingAnInterfaceBetweenTwoVolumes(const G4Step* aStep,const G4String& vol1_name,const G4String& vol2_name,G4ThreeVector& , G4double& , G4bool& GoingIn)
00216 {
00217 G4bool step_at_boundary = (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
00218 G4bool did_cross =false;
00219 if (step_at_boundary){
00220 const G4VTouchable* postStepTouchable = aStep->GetPostStepPoint()->GetTouchable();
00221 const G4VTouchable* preStepTouchable = aStep->GetPreStepPoint()->GetTouchable();
00222 if (preStepTouchable && postStepTouchable){
00223
00224 G4String post_vol_name = postStepTouchable->GetVolume()->GetName();
00225 if (post_vol_name =="") post_vol_name = postStepTouchable->GetVolume()->GetLogicalVolume()->GetName();
00226 G4String pre_vol_name = preStepTouchable->GetVolume()->GetName();
00227 if (pre_vol_name =="") pre_vol_name = preStepTouchable->GetVolume()->GetLogicalVolume()->GetName();
00228
00229
00230 if ( pre_vol_name == vol1_name && post_vol_name == vol2_name){
00231 GoingIn=true;
00232 did_cross=true;
00233 }
00234 else if (pre_vol_name == vol2_name && post_vol_name == vol1_name){
00235 GoingIn=false;
00236 did_cross=true;
00237 }
00238
00239 }
00240 }
00241 return did_cross;
00242
00243 }
00244
00246
00247 G4bool G4AdjointCrossSurfChecker::AddaSphericalSurface(const G4String& SurfaceName, G4double radius, G4ThreeVector pos, G4double& Area)
00248 {
00249 G4int ind = FindRegisteredSurface(SurfaceName);
00250 Area= 4.*pi*radius*radius;
00251 if (ind>=0) {
00252 ListOfSurfaceType[ind]="Sphere";
00253 ListOfSphereRadius[ind]=radius;
00254 ListOfSphereCenter[ind]=pos;
00255 ListOfVol1Name[ind]="";
00256 ListOfVol2Name[ind]="";
00257 AreaOfSurface[ind]=Area;
00258 }
00259 else {
00260 ListOfSurfaceName.push_back(SurfaceName);
00261 ListOfSurfaceType.push_back("Sphere");
00262 ListOfSphereRadius.push_back(radius);
00263 ListOfSphereCenter.push_back(pos);
00264 ListOfVol1Name.push_back("");
00265 ListOfVol2Name.push_back("");
00266 AreaOfSurface.push_back(Area);
00267 }
00268 return true;
00269 }
00271
00272 G4bool G4AdjointCrossSurfChecker::AddaSphericalSurfaceWithCenterAtTheCenterOfAVolume(const G4String& SurfaceName, G4double radius, const G4String& volume_name, G4ThreeVector & center, G4double& area)
00273 {
00274
00275 G4VPhysicalVolume* thePhysicalVolume = 0;
00276 G4PhysicalVolumeStore* thePhysVolStore =G4PhysicalVolumeStore::GetInstance();
00277 for ( unsigned int i=0; i< thePhysVolStore->size();i++){
00278 if ((*thePhysVolStore)[i]->GetName() == volume_name){
00279 thePhysicalVolume = (*thePhysVolStore)[i];
00280 };
00281
00282 }
00283 if (thePhysicalVolume){
00284 G4VPhysicalVolume* daughter =thePhysicalVolume;
00285 G4LogicalVolume* mother = thePhysicalVolume->GetMotherLogical();
00286 G4AffineTransform theTransformationFromPhysVolToWorld = G4AffineTransform();
00287 while (mother){
00288 theTransformationFromPhysVolToWorld *=
00289 G4AffineTransform(daughter->GetFrameRotation(),daughter->GetObjectTranslation());
00290
00291
00292
00293
00294 for ( unsigned int i=0; i< thePhysVolStore->size();i++){
00295 if ((*thePhysVolStore)[i]->GetLogicalVolume() == mother){
00296 daughter = (*thePhysVolStore)[i];
00297 mother =daughter->GetMotherLogical();
00298 break;
00299 };
00300 }
00301
00302 }
00303 center = theTransformationFromPhysVolToWorld.NetTranslation();
00304 G4cout<<"Center of the spherical surface is at the position: "<<center/cm<<" cm"<<std::endl;
00305
00306 }
00307 else {
00308 G4cout<<"The physical volume with name "<<volume_name<<" does not exist!!"<<std::endl;
00309 return false;
00310 }
00311 return AddaSphericalSurface(SurfaceName, radius, center, area);
00312
00313 }
00315
00316 G4bool G4AdjointCrossSurfChecker::AddanExtSurfaceOfAvolume(const G4String& SurfaceName, const G4String& volume_name, G4double& Area)
00317 {
00318 G4int ind = FindRegisteredSurface(SurfaceName);
00319
00320 G4VPhysicalVolume* thePhysicalVolume = 0;
00321 G4PhysicalVolumeStore* thePhysVolStore =G4PhysicalVolumeStore::GetInstance();
00322 for ( unsigned int i=0; i< thePhysVolStore->size();i++){
00323 if ((*thePhysVolStore)[i]->GetName() == volume_name){
00324 thePhysicalVolume = (*thePhysVolStore)[i];
00325 };
00326
00327 }
00328 if (!thePhysicalVolume){
00329 G4cout<<"The physical volume with name "<<volume_name<<" does not exist!!"<<std::endl;
00330 return false;
00331 }
00332 Area = thePhysicalVolume->GetLogicalVolume()->GetSolid()->GetSurfaceArea();
00333 G4String mother_vol_name = "";
00334 G4LogicalVolume* theMother = thePhysicalVolume->GetMotherLogical();
00335
00336 if (theMother) mother_vol_name= theMother->GetName();
00337 if (ind>=0) {
00338 ListOfSurfaceType[ind]="ExternalSurfaceOfAVolume";
00339 ListOfSphereRadius[ind]=0.;
00340 ListOfSphereCenter[ind]=G4ThreeVector(0.,0.,0.);
00341 ListOfVol1Name[ind]=volume_name;
00342 ListOfVol2Name[ind]=mother_vol_name;
00343 AreaOfSurface[ind]=Area;
00344 }
00345 else {
00346 ListOfSurfaceName.push_back(SurfaceName);
00347 ListOfSurfaceType.push_back("ExternalSurfaceOfAVolume");
00348 ListOfSphereRadius.push_back(0.);
00349 ListOfSphereCenter.push_back(G4ThreeVector(0.,0.,0.));
00350 ListOfVol1Name.push_back(volume_name);
00351 ListOfVol2Name.push_back(mother_vol_name);
00352 AreaOfSurface.push_back(Area);
00353 }
00354 return true;
00355 }
00357
00358 G4bool G4AdjointCrossSurfChecker::AddanInterfaceBetweenTwoVolumes(const G4String& SurfaceName, const G4String& volume_name1, const G4String& volume_name2,G4double& Area)
00359 {
00360 G4int ind = FindRegisteredSurface(SurfaceName);
00361 Area=-1.;
00362 if (ind>=0) {
00363 ListOfSurfaceType[ind]="BoundaryBetweenTwoVolumes";
00364 ListOfSphereRadius[ind]=0.;
00365 ListOfSphereCenter[ind]=G4ThreeVector(0.,0.,0.);
00366 ListOfVol1Name[ind]=volume_name1;
00367 ListOfVol2Name[ind]=volume_name2;
00368 AreaOfSurface[ind]=Area;
00369
00370 }
00371 else {
00372 ListOfSurfaceName.push_back(SurfaceName);
00373 ListOfSurfaceType.push_back("BoundaryBetweenTwoVolumes");
00374 ListOfSphereRadius.push_back(0.);
00375 ListOfSphereCenter.push_back(G4ThreeVector(0.,0.,0.));
00376 ListOfVol1Name.push_back(volume_name1);
00377 ListOfVol2Name.push_back(volume_name2);
00378 AreaOfSurface.push_back(Area);
00379 }
00380 return true;
00381 }
00383
00384 void G4AdjointCrossSurfChecker::ClearListOfSelectedSurface()
00385 {
00386 ListOfSurfaceName.clear();
00387 ListOfSurfaceType.clear();
00388 ListOfSphereRadius.clear();
00389 ListOfSphereCenter.clear();
00390 ListOfVol1Name.clear();
00391 ListOfVol2Name.clear();
00392 }
00394
00395 G4int G4AdjointCrossSurfChecker::FindRegisteredSurface(const G4String& name)
00396 {
00397 G4int ind=-1;
00398 for (size_t i = 0; i<ListOfSurfaceName.size();i++){
00399 if (name == ListOfSurfaceName[i]) {
00400 ind = int (i);
00401 return ind;
00402 }
00403 }
00404 return ind;
00405 }
00406
