00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 #include "SimpleGeometrySvc.hh"
00016
00017 #include "myglobals.hh"
00018 #include "G4LogicalVolume.hh"
00019 #include "G4VSensitiveDetector.hh"
00020 #include "G4Box.hh"
00021 #include "G4EllipticalTube.hh"
00022 #include "G4Tubs.hh"
00023 #include "G4Torus.hh"
00024 #include "G4Hype.hh"
00025 #include "G4TwistedTubs.hh"
00026 #include "G4Cons.hh"
00027 #include "G4Polycone.hh"
00028 #include "G4ExtrudedSolid.hh"
00029 #include "G4Sphere.hh"
00030 #include "G4IntersectionSolid.hh"
00031 #include "G4SubtractionSolid.hh"
00032 #include "G4UnionSolid.hh"
00033 #include "G4SolidStore.hh"
00034 #include "G4PVPlacement.hh"
00035 #include "G4RotationMatrix.hh"
00036 #include "G4VisAttributes.hh"
00037 #include "G4UserLimits.hh"
00038
00039 #include <iostream>
00040 #include <vector>
00041
00042 #include "MyVGeometrySvc.hh"
00043 #include "SimpleGeometryParameter.hh"
00044 #include "MyDetectorManager.hh"
00045
00046 SimpleGeometrySvc::SimpleGeometrySvc(G4String name, G4String opt )
00047 : MyVGeometrySvc(name, "S")
00048 {
00049 if ( opt == "" ){
00050 SimpleGeometryParameter *pointer = new SimpleGeometryParameter(name);
00051
00052 set_GeometryParameter(pointer);
00053 }
00054 }
00055
00056 SimpleGeometrySvc::~SimpleGeometrySvc(){
00057 printf("~SimpleGeometrySvc");
00058 }
00059
00060
00061 void SimpleGeometrySvc::set_GeometryParameter( SimpleGeometryParameter* val ){
00062 m_GeometryParameter = val;
00063 MyVGeometrySvc::set_GeometryParameter(m_GeometryParameter);
00064 }
00065
00066
00067
00068 G4VPhysicalVolume* SimpleGeometrySvc::SetGeometry(){
00069 m_GeometryParameter->Dump();
00070 ConstructVolumes();
00071 G4VPhysicalVolume *current = PlaceVolumes();
00072 return current;
00073 }
00074
00075
00076
00077 void SimpleGeometrySvc::ConstructVolumes(){
00078 MyVGeometrySvc::ConstructVolumes();
00079 bool vis;
00080 double r, g, b, t;
00081 std::stringstream buffer;
00082 G4Material* pttoMaterial;
00083
00084 int nVol = m_GeometryParameter->get_VolNo();
00085 for ( int i_Vol = 0; i_Vol < nVol; i_Vol++ ){
00086 G4String name, mat_name, SDName;
00087 G4String SolidType;
00088 G4int SolidIndex;
00089 G4int SRepNo;
00090 G4int RepNo;
00091 G4double max_step_length = 0;
00092 SolidIndex = m_GeometryParameter->get_SolidIndex(i_Vol);
00093 SolidType = m_GeometryParameter->get_SolidType(i_Vol);
00094 name = m_GeometryParameter->get_name(i_Vol);
00095 SRepNo = m_GeometryParameter->get_SRepNo(i_Vol);
00096 RepNo = m_GeometryParameter->get_RepNo(i_Vol);
00097 for ( int i = 0; i < RepNo; i++ ){
00098 G4VSolid* sol_Vol = 0;
00099 G4String iname = name;
00100 if (RepNo>1){
00101 buffer.str("");
00102 buffer.clear();
00103 buffer<<name<<"_"<<i;
00104 iname = buffer.str();
00105 }
00106 if ( SolidType == "Box" ){
00107 G4double halfX, halfY, halfZ;
00108 halfX = m_GeometryParameter->get_Box_X(SolidIndex,i)/2;
00109 halfY = m_GeometryParameter->get_Box_Y(SolidIndex,i)/2;
00110 halfZ = m_GeometryParameter->get_Box_Z(SolidIndex,i)/2;
00111 sol_Vol=new G4Box(iname,halfX,halfY,halfZ);
00112 max_step_length = 0.25*std::min(halfX, std::min(halfY, halfZ));
00113 }
00114 else if ( SolidType == "EllipticalTube" ){
00115 G4double halfX, halfY, halfZ;
00116 halfX = m_GeometryParameter->get_EllipticalTube_X(SolidIndex,i)/2;
00117 halfY = m_GeometryParameter->get_EllipticalTube_Y(SolidIndex,i)/2;
00118 halfZ = m_GeometryParameter->get_EllipticalTube_Z(SolidIndex,i)/2;
00119 sol_Vol=new G4EllipticalTube(iname,halfX,halfY,halfZ);
00120 max_step_length = 0.25*std::min(halfX, std::min(halfY, halfZ));
00121 }
00122 else if ( SolidType == "Tubs" ){
00123 G4double RMax, RMin, halfLength, StartAng, SpanAng;
00124 RMax = m_GeometryParameter->get_Tubs_RMax(SolidIndex,i);
00125 RMin = m_GeometryParameter->get_Tubs_RMin(SolidIndex,i);
00126 halfLength = m_GeometryParameter->get_Tubs_Length(SolidIndex,i)/2;
00127 StartAng = m_GeometryParameter->get_Tubs_StartAng(SolidIndex,i);
00128 SpanAng = m_GeometryParameter->get_Tubs_SpanAng(SolidIndex,i);
00129 sol_Vol=new G4Tubs(iname,RMin,RMax,halfLength,StartAng,SpanAng);
00130 max_step_length = 0.25*RMin;
00131 }
00132 else if ( SolidType == "Torus" ){
00133 G4double RMax, RMin, Rtor, StartAng, SpanAng;
00134 RMax = m_GeometryParameter->get_Torus_RMax(SolidIndex,i);
00135 RMin = m_GeometryParameter->get_Torus_RMin(SolidIndex,i);
00136 Rtor = m_GeometryParameter->get_Torus_Rtor(SolidIndex,i);
00137 StartAng = m_GeometryParameter->get_Torus_StartAng(SolidIndex,i);
00138 SpanAng = m_GeometryParameter->get_Torus_SpanAng(SolidIndex,i);
00139 sol_Vol=new G4Torus(iname,RMin,RMax,Rtor,StartAng,SpanAng);
00140 max_step_length = 0.25*RMin;
00141 }
00142 else if ( SolidType == "Sphere" ){
00143 G4double RMax, RMin, StartPhi, SpanPhi, StartTheta, SpanTheta;
00144 RMax = m_GeometryParameter->get_Sphere_RMax(SolidIndex,i);
00145 RMin = m_GeometryParameter->get_Sphere_RMin(SolidIndex,i);
00146 StartPhi = m_GeometryParameter->get_Sphere_StartPhi(SolidIndex,i);
00147 SpanPhi = m_GeometryParameter->get_Sphere_SpanPhi(SolidIndex,i);
00148 StartTheta = m_GeometryParameter->get_Sphere_StartTheta(SolidIndex,i);
00149 SpanTheta = m_GeometryParameter->get_Sphere_SpanTheta(SolidIndex,i);
00150 sol_Vol=new G4Sphere(iname,RMin,RMax,StartPhi,SpanPhi,StartTheta,SpanTheta);
00151 max_step_length = 0.25*RMin;
00152 }
00153 else if ( SolidType == "Hype" ){
00154 G4double innerRadius, outerRadius, halfLength, innerStereo, outerStereo;
00155 innerRadius = m_GeometryParameter->get_Hype_innerRadius(SolidIndex,i);
00156 outerRadius = m_GeometryParameter->get_Hype_outerRadius(SolidIndex,i);
00157 halfLength = m_GeometryParameter->get_Hype_Length(SolidIndex,i)/2;
00158 innerStereo = m_GeometryParameter->get_Hype_innerStereo(SolidIndex,i);
00159 outerStereo = m_GeometryParameter->get_Hype_outerStereo(SolidIndex,i);
00160 sol_Vol=new G4Hype(iname,innerRadius,outerRadius,innerStereo,outerStereo,halfLength);
00161 }
00162 else if ( SolidType == "TwistedTubs" ){
00163 G4double endinnerrad, endouterrad, halfLength, twistedangle, dphi;
00164 twistedangle = m_GeometryParameter->get_TwistedTubs_twistedangle(SolidIndex,i);
00165 endinnerrad = m_GeometryParameter->get_TwistedTubs_endinnerrad(SolidIndex,i);
00166 endouterrad = m_GeometryParameter->get_TwistedTubs_endouterrad(SolidIndex,i);
00167 halfLength = m_GeometryParameter->get_TwistedTubs_Length(SolidIndex,i)/2;
00168 dphi = m_GeometryParameter->get_TwistedTubs_dphi(SolidIndex,i);
00169 sol_Vol=new G4TwistedTubs(iname,twistedangle,endinnerrad,endouterrad,halfLength,dphi);
00170 }
00171 else if ( SolidType == "Cons" ){
00172 G4double RMax1, RMin1, RMax2, RMin2, halfLength, StartAng, SpanAng;
00173 RMax1 = m_GeometryParameter->get_Cons_RMax1(SolidIndex,i);
00174 RMin1 = m_GeometryParameter->get_Cons_RMin1(SolidIndex,i);
00175 RMax2 = m_GeometryParameter->get_Cons_RMax2(SolidIndex,i);
00176 RMin2 = m_GeometryParameter->get_Cons_RMin2(SolidIndex,i);
00177 halfLength = m_GeometryParameter->get_Cons_Length(SolidIndex,i)/2;
00178 StartAng = m_GeometryParameter->get_Cons_StartAng(SolidIndex,i);
00179 SpanAng = m_GeometryParameter->get_Cons_SpanAng(SolidIndex,i);
00180 sol_Vol=new G4Cons(iname,RMin1,RMax1,RMin2,RMax2,halfLength,StartAng,SpanAng);
00181 max_step_length = 0.25*std::min(RMin2, RMin2);
00182
00183 }
00184 else if ( SolidType == "Polycone" ){
00185 G4double StartAng, SpanAng;
00186 G4int numZ;
00187 numZ = m_GeometryParameter->get_Polycone_numZ(SolidIndex,i);
00188 G4double *RMax = new G4double[numZ];
00189 G4double *RMin = new G4double[numZ];
00190 G4double *Z = new G4double[numZ];
00191 for ( int j = 0; j < numZ; j++ ){
00192 RMax[j] = m_GeometryParameter->get_Polycone_RMax(SolidIndex,j,i);
00193 RMin[j] = m_GeometryParameter->get_Polycone_RMin(SolidIndex,j,i);
00194 Z[j] = m_GeometryParameter->get_Polycone_Z(SolidIndex,j,i);
00195 std::cout<<j<<": RMax = "<<RMax[j]/mm<<"mm, RMin = "<<RMin[j]/mm<<"mm, "<<Z[j]/mm<<"mm"<<std::endl;
00196 }
00197 StartAng = m_GeometryParameter->get_Polycone_StartAng(SolidIndex,i);
00198 SpanAng = m_GeometryParameter->get_Polycone_SpanAng(SolidIndex,i);
00199 sol_Vol=new G4Polycone(iname,StartAng,SpanAng,numZ,Z,RMin,RMax);
00200 }
00201 else if ( SolidType == "ExtrudedSolid" ){
00202 G4int numZ, numP;
00203 numZ = m_GeometryParameter->get_ExtrudedSolid_numZ(SolidIndex,i);
00204 numP = m_GeometryParameter->get_ExtrudedSolid_numP(SolidIndex,i);
00205 std::vector<G4TwoVector> polygon;
00206 std::vector<G4ExtrudedSolid::ZSection> zsections;
00207 for ( int j = 0; j < numZ; j++ ){
00208 zsections.push_back(G4ExtrudedSolid::ZSection(m_GeometryParameter->get_ExtrudedSolid_Z(SolidIndex,j,i),
00209 G4TwoVector(m_GeometryParameter->get_ExtrudedSolid_X0(SolidIndex,j,i),
00210 m_GeometryParameter->get_ExtrudedSolid_Y0(SolidIndex,j,i)),
00211 m_GeometryParameter->get_ExtrudedSolid_Scale(SolidIndex,j,i)));
00212 }
00213 for ( int j = 0; j < numP; j++ ){
00214 polygon.push_back(G4TwoVector(m_GeometryParameter->get_ExtrudedSolid_X(SolidIndex,j,i),
00215 m_GeometryParameter->get_ExtrudedSolid_Y(SolidIndex,j,i)));
00216 }
00217 sol_Vol=new G4ExtrudedSolid(iname,polygon,zsections);
00218 }
00219 else if ( SolidType == "BooleanSolid" ){
00220 G4double Ephi = m_GeometryParameter->get_BooleanSolid_Ephi(SolidIndex,i);
00221 G4double Etheta = m_GeometryParameter->get_BooleanSolid_Etheta(SolidIndex,i);
00222 G4double Epsi = m_GeometryParameter->get_BooleanSolid_Epsi(SolidIndex,i);
00223 G4double PosX = m_GeometryParameter->get_BooleanSolid_PosX(SolidIndex,i);
00224 G4double PosY = m_GeometryParameter->get_BooleanSolid_PosY(SolidIndex,i);
00225 G4double PosZ = m_GeometryParameter->get_BooleanSolid_PosZ(SolidIndex,i);
00226 G4RotationMatrix* rot =new G4RotationMatrix(Ephi,Etheta,Epsi);
00227 G4ThreeVector trans(PosX ,PosY ,PosZ);
00228 G4String type = m_GeometryParameter->get_BooleanSolid_type(SolidIndex);
00229 G4String sol_name1 = m_GeometryParameter->get_BooleanSolid_sol1(SolidIndex);
00230 G4String sol_name2 = m_GeometryParameter->get_BooleanSolid_sol2(SolidIndex);
00231 G4String isol_name1 = sol_name1;
00232 G4String isol_name2 = sol_name2;
00233 if (RepNo>1){
00234 buffer.str("");
00235 buffer.clear();
00236 buffer<<sol_name1<<"_"<<i;
00237 isol_name1 = buffer.str();
00238 buffer.str("");
00239 buffer.clear();
00240 buffer<<sol_name2<<"_"<<i;
00241 isol_name2 = buffer.str();
00242 }
00243 G4SolidStore *pSolidStore = G4SolidStore::GetInstance();
00244 G4VSolid *sol1 = pSolidStore->GetSolid(isol_name1);
00245 G4VSolid *sol2 = pSolidStore->GetSolid(isol_name2);
00246 if (!sol1||!sol2){
00247 std::cout<<"ERROR: in SimpleGeometrySvc::ConstructVolumes(), cannot find solids for BooleanSolid \""<<iname<<"\" !!!"<<std::endl;
00248 continue;
00249 }
00250 if ( type == "plus" ){
00251 sol_Vol = new G4UnionSolid(iname,sol1,sol2,rot,trans);
00252 }
00253 else if ( type == "minus" ){
00254 sol_Vol = new G4SubtractionSolid(iname,sol1,sol2,rot,trans);
00255 }
00256 else if ( type == "times" ){
00257 sol_Vol = new G4IntersectionSolid(iname,sol1,sol2,rot,trans);
00258 }
00259 else {
00260 std::cout<<"ERROR: in SimpleGeometrySvc::ConstructVolumes(), cannot recognize BooleanSolid type \""<<type<<"\""<<"for \""<<iname<<"\"!!!"<<std::endl;
00261 continue;
00262 }
00263 }
00264 else {
00265 std::cout<<"SolidType "<<SolidType<<" is not supported yet!"<<std::endl;
00266 G4Exception("SimpleGeometrySvc::ConstructVolumes()",
00267 "InvalidSetup", FatalException,
00268 "unknown SolidType");
00269 }
00270 if (m_GeometryParameter->get_SolidBoolean(i_Vol)) continue;
00271 SDName = m_GeometryParameter->get_SDName(i_Vol);
00272 mat_name = m_GeometryParameter->get_material(i_Vol);
00273 pttoMaterial = G4Material::GetMaterial(mat_name);
00274
00275
00276
00277
00278
00279 if (!pttoMaterial){
00280 std::cout<<"Material "<<mat_name<<" is not defined!"<<std::endl;
00281 G4Exception("SimpleGeometrySvc::ConstructVolumes()",
00282 "InvalidSetup", FatalException,
00283 "unknown material name");
00284 }
00285 G4LogicalVolume* log_Vol;
00286 G4UserLimits* user_limit;
00287 if (max_step_length != 0) {
00288 user_limit = new G4UserLimits(max_step_length);
00289
00290 }
00291 else {
00292 user_limit = new G4UserLimits();
00293 }
00294 log_Vol = new G4LogicalVolume(sol_Vol, pttoMaterial, iname,0,0,user_limit);
00295 G4VSensitiveDetector* aSD;
00296 aSD = MyDetectorManager::GetMyDetectorManager()->GetSD(iname, SDName, const_cast<SimpleGeometrySvc*>(this) );
00297 log_Vol->SetSensitiveDetector(aSD);
00298
00299 vis = m_GeometryParameter->get_vis(i_Vol);
00300 if (!vis){
00301 log_Vol->SetVisAttributes(G4VisAttributes::Invisible);
00302 }
00303 else{
00304 r = m_GeometryParameter->get_r(i_Vol);
00305 g = m_GeometryParameter->get_g(i_Vol);
00306 b = m_GeometryParameter->get_b(i_Vol);
00307 t = m_GeometryParameter->get_t(i_Vol);
00308 G4VisAttributes* visAttributes = new G4VisAttributes;
00309 visAttributes->SetVisibility(true);
00310 visAttributes->SetColour(r,g,b,t);
00311 log_Vol->SetVisAttributes(visAttributes);
00312 }
00313 }
00314 }
00315 }
00316
00317
00318 G4VPhysicalVolume* SimpleGeometrySvc::PlaceVolumes(){
00319 std::stringstream buffer;
00320 G4VPhysicalVolume* world_pvol = 0;
00321 bool checkOverlap = m_GeometryParameter->get_checkoverlap();
00322 int nVol = m_GeometryParameter->get_VolNo();
00323 for ( int i_Vol = 0; i_Vol < nVol; i_Vol++ ){
00324 if (m_GeometryParameter->get_SolidBoolean(i_Vol)) continue;
00325 G4double PosX, PosY, PosZ;
00326 G4String name, motVolName;
00327 G4int SRepNo;
00328 G4int RepNo;
00329 name = m_GeometryParameter->get_name(i_Vol);
00330
00331 motVolName = m_GeometryParameter->get_MotherName(i_Vol);
00332 SRepNo = m_GeometryParameter->get_SRepNo(i_Vol);
00333 RepNo = m_GeometryParameter->get_RepNo(i_Vol);
00334 G4LogicalVolume* log_container;
00335 int motRepNo = 1;
00336 if ( name == "World" ) motRepNo = 1;
00337 else{
00338 SimpleGeometryParameter* motPara = MyDetectorManager::GetMyDetectorManager()->GetParaFromVolume(motVolName);
00339 if (motPara)
00340 motRepNo = motPara->get_RepNo(motPara->get_VolIndex(motVolName));
00341 else{
00342 std::cout<<"We can't find parameter pointer!!"<<std::endl;
00343 }
00344 }
00345 for ( int im = 0; im < motRepNo; im++ ){
00346 if ( name == "World" ) log_container = 0;
00347 else {
00348 G4String imotVolName = motVolName;
00349 if (motRepNo>1){
00350 buffer.str("");
00351 buffer.clear();
00352 buffer<<motVolName<<"_"<<im;
00353 imotVolName = buffer.str();
00354 }
00355 log_container = get_logicalVolume(imotVolName);
00356 }
00357 for ( int i = 0; i < RepNo; i++ ){
00358 G4String iname = name;
00359 if (RepNo>1){
00360 buffer.str("");
00361 buffer.clear();
00362 buffer<<name<<"_"<<i;
00363 iname = buffer.str();
00364 }
00365 G4LogicalVolume* log_Vol = get_logicalVolume(iname);
00366 PosX = m_GeometryParameter->get_PosX(i_Vol,i);
00367 PosY = m_GeometryParameter->get_PosY(i_Vol,i);
00368 PosZ = m_GeometryParameter->get_PosZ(i_Vol,i);
00369 G4ThreeVector pos(PosX ,PosY ,PosZ);
00370 G4double Ephi = m_GeometryParameter->get_Ephi(i_Vol,i);
00371 G4double Etheta = m_GeometryParameter->get_Etheta(i_Vol,i);
00372 G4double Epsi = m_GeometryParameter->get_Epsi(i_Vol,i);
00373 G4RotationMatrix* rotateMatrix=new G4RotationMatrix(Ephi,Etheta,Epsi);
00374 if ( name == "World" ) rotateMatrix = 0;
00375 G4VPhysicalVolume* phy_Vol =
00376 new G4PVPlacement(rotateMatrix,
00377 pos,
00378 log_Vol,
00379 name,
00380 log_container,
00381 false,
00382 i+SRepNo,
00383 checkOverlap);
00384 if ( name == "World" ) world_pvol = phy_Vol;
00385 }
00386 }
00387 }
00388
00389 G4VPhysicalVolume *former = MyVGeometrySvc::PlaceVolumes();
00390 if (!world_pvol) world_pvol=former;
00391 return world_pvol;
00392 }
