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 // 26 // 27 // 27 // 28 // by I.Hrivnacova, V.Berejnoi 13.10.99 28 // by I.Hrivnacova, V.Berejnoi 13.10.99 29 29 30 #include <assert.h> 30 #include <assert.h> 31 31 32 #include "G3Division.hh" 32 #include "G3Division.hh" 33 #include "G3VolTableEntry.hh" 33 #include "G3VolTableEntry.hh" 34 #include "G3toG4MakeSolid.hh" 34 #include "G3toG4MakeSolid.hh" 35 #include "G4Para.hh" 35 #include "G4Para.hh" 36 #include "G3Pos.hh" 36 #include "G3Pos.hh" 37 #include "G4SystemOfUnits.hh" 37 #include "G4SystemOfUnits.hh" 38 #include "G4LogicalVolume.hh" 38 #include "G4LogicalVolume.hh" 39 #include "G4VPhysicalVolume.hh" 39 #include "G4VPhysicalVolume.hh" 40 #include "G4PVPlacement.hh" 40 #include "G4PVPlacement.hh" 41 #include "G4PVReplica.hh" 41 #include "G4PVReplica.hh" 42 #ifndef G3G4_NO_REFLECTION 42 #ifndef G3G4_NO_REFLECTION 43 #include "G4ReflectionFactory.hh" 43 #include "G4ReflectionFactory.hh" 44 #endif 44 #endif 45 45 46 G3VolTableEntry* G4CreateVTE(G4String vname, G 46 G3VolTableEntry* G4CreateVTE(G4String vname, G4String shape, G4int nmed, 47 G4double Rpar[] 47 G4double Rpar[], G4int npar); 48 48 49 G3Division::G3Division(G3DivType type, G3VolTa 49 G3Division::G3Division(G3DivType type, G3VolTableEntry* vte, 50 G3VolTableEntry* mvte, G4int n 50 G3VolTableEntry* mvte, G4int nofDivisions, 51 G4int iaxis, G4int nmed, G4double c0, G4do 51 G4int iaxis, G4int nmed, G4double c0, G4double step) 52 : fType(type), 52 : fType(type), 53 fVTE(vte), 53 fVTE(vte), 54 fMVTE(mvte), 54 fMVTE(mvte), 55 fNofDivisions(nofDivisions), 55 fNofDivisions(nofDivisions), 56 fIAxis(iaxis), 56 fIAxis(iaxis), 57 fNmed(nmed), 57 fNmed(nmed), 58 fC0(c0), 58 fC0(c0), 59 fStep(step), 59 fStep(step), 60 fLowRange(0.), 60 fLowRange(0.), 61 fHighRange(0.), 61 fHighRange(0.), 62 fWidth(0.), 62 fWidth(0.), 63 fOffset(0.), 63 fOffset(0.), 64 fAxis(kXAxis) 64 fAxis(kXAxis) 65 { 65 { 66 fVTE->SetHasNegPars(true); 66 fVTE->SetHasNegPars(true); 67 } 67 } 68 68 69 G3Division::G3Division(G3VolTableEntry* vte, G 69 G3Division::G3Division(G3VolTableEntry* vte, G3VolTableEntry* mvte, 70 const G3Division& divis 70 const G3Division& division) 71 : fVTE(vte), 71 : fVTE(vte), 72 fMVTE(mvte) 72 fMVTE(mvte) 73 { 73 { 74 // only "input" parameters are copied from d 74 // only "input" parameters are copied from division 75 fType = division.fType; 75 fType = division.fType; 76 fNofDivisions = division.fNofDivisions; 76 fNofDivisions = division.fNofDivisions; 77 fIAxis = division.fIAxis; 77 fIAxis = division.fIAxis; 78 fNmed = division.fNmed; 78 fNmed = division.fNmed; 79 fC0 = division.fC0; 79 fC0 = division.fC0; 80 fStep = division.fStep; 80 fStep = division.fStep; 81 81 82 // other parameters are set as in standard c 82 // other parameters are set as in standard constructor 83 fLowRange = 0.; 83 fLowRange = 0.; 84 fHighRange = 0.; 84 fHighRange = 0.; 85 fWidth = 0.; 85 fWidth = 0.; 86 fOffset = 0.; 86 fOffset = 0.; 87 fAxis = kXAxis; 87 fAxis = kXAxis; 88 fVTE->SetHasNegPars(true); 88 fVTE->SetHasNegPars(true); 89 } 89 } 90 90 91 G3Division::~G3Division() 91 G3Division::~G3Division() 92 {} 92 {} 93 93 94 // public methods 94 // public methods 95 95 96 void G3Division::UpdateVTE() 96 void G3Division::UpdateVTE() 97 { 97 { 98 if (fVTE->HasNegPars() && !(fMVTE->HasNegPar 98 if (fVTE->HasNegPars() && !(fMVTE->HasNegPars())) { 99 99 100 // set nmed from mother 100 // set nmed from mother 101 if (fNmed == 0) fNmed = fMVTE->GetNmed(); 101 if (fNmed == 0) fNmed = fMVTE->GetNmed(); 102 fVTE->SetNmed(fNmed); 102 fVTE->SetNmed(fNmed); 103 103 104 SetRangeAndAxis(); 104 SetRangeAndAxis(); 105 105 106 // create envelope (if necessary) 106 // create envelope (if necessary) 107 // and solid 107 // and solid 108 G3VolTableEntry* envVTE = 0; 108 G3VolTableEntry* envVTE = 0; 109 if (fType == kDvn) envVTE = Dvn(); 109 if (fType == kDvn) envVTE = Dvn(); 110 else if (fType == kDvn2) envVTE = Dvn2(); 110 else if (fType == kDvn2) envVTE = Dvn2(); 111 else if (fType == kDvt) envVTE = Dvt(); 111 else if (fType == kDvt) envVTE = Dvt(); 112 else if (fType == kDvt2) envVTE = Dvt2(); 112 else if (fType == kDvt2) envVTE = Dvt2(); 113 113 114 if (envVTE) { 114 if (envVTE) { 115 // reset mother <-> daughter 115 // reset mother <-> daughter 116 fMVTE->ReplaceDaughter(fVTE, envVTE); 116 fMVTE->ReplaceDaughter(fVTE, envVTE); 117 fVTE->ReplaceMother(fMVTE, envVTE); 117 fVTE->ReplaceMother(fMVTE, envVTE); 118 envVTE->AddDaughter(fVTE); 118 envVTE->AddDaughter(fVTE); 119 envVTE->AddMother(fMVTE); 119 envVTE->AddMother(fMVTE); 120 120 121 // replace mother with envelope 121 // replace mother with envelope 122 fMVTE = envVTE; 122 fMVTE = envVTE; 123 } 123 } 124 } 124 } 125 } 125 } 126 126 127 void G3Division::CreatePVReplica() 127 void G3Division::CreatePVReplica() 128 { 128 { 129 G4String name = fVTE->GetName(); 129 G4String name = fVTE->GetName(); 130 G4LogicalVolume* lv = fVTE->GetLV(); 130 G4LogicalVolume* lv = fVTE->GetLV(); 131 G4LogicalVolume* mlv = fMVTE->GetLV(); 131 G4LogicalVolume* mlv = fMVTE->GetLV(); 132 132 133 G4String shape = fMVTE->GetShape(); 133 G4String shape = fMVTE->GetShape(); 134 if (shape == "PARA") { 134 if (shape == "PARA") { 135 // The para volume cannot be replicated us 135 // The para volume cannot be replicated using G4PVReplica. 136 // (Replicating a volume along a cartesian 136 // (Replicating a volume along a cartesian axis means "slicing" it 137 // with slices -perpendicular- to that axi 137 // with slices -perpendicular- to that axis.) 138 138 139 // position the replicated elements 139 // position the replicated elements 140 for (G4int i=0; i<fNofDivisions; i++) { 140 for (G4int i=0; i<fNofDivisions; i++) { 141 G4ThreeVector position = G4ThreeVector( 141 G4ThreeVector position = G4ThreeVector(); 142 position[fIAxis-1] = fLowRange + fWidth 142 position[fIAxis-1] = fLowRange + fWidth/2. + i*fWidth; 143 if (position.y()!=0.) 143 if (position.y()!=0.) 144 position.setX(position.y()*((G4Para*) 144 position.setX(position.y()*((G4Para*)lv->GetSolid())->GetTanAlpha()); 145 145 146 #ifndef G3G4_NO_REFLECTION 146 #ifndef G3G4_NO_REFLECTION 147 G4ReflectionFactory::Instance() 147 G4ReflectionFactory::Instance() 148 ->Place(G4Translate3D(position), name 148 ->Place(G4Translate3D(position), name, lv, mlv, 0, i); 149 149 150 #else 150 #else 151 new G4PVPlacement(0, position, lv, name 151 new G4PVPlacement(0, position, lv, name, mlv, 0, i); 152 152 153 #endif 153 #endif 154 } 154 } 155 155 156 // G4PVReplica cannot be created 156 // G4PVReplica cannot be created 157 return; 157 return; 158 } 158 } 159 159 160 #ifdef G3G4DEBUG 160 #ifdef G3G4DEBUG 161 G4cout << "Create G4PVReplica name " << na 161 G4cout << "Create G4PVReplica name " << name << " logical volume name " 162 << lv->GetName() << " mother logical volm 162 << lv->GetName() << " mother logical volme name " 163 << mlv->GetName() << " axis " << fAxis << 163 << mlv->GetName() << " axis " << fAxis << " ndivisions " 164 << fNofDivisions << " width " << fWidth < 164 << fNofDivisions << " width " << fWidth << " Offset " 165 << fOffset << G4endl; 165 << fOffset << G4endl; 166 #endif 166 #endif 167 167 168 #ifndef G3G4_NO_REFLECTION 168 #ifndef G3G4_NO_REFLECTION 169 G4ReflectionFactory::Instance() 169 G4ReflectionFactory::Instance() 170 ->Replicate(name, lv, mlv, fAxis, fNofDivi 170 ->Replicate(name, lv, mlv, fAxis, fNofDivisions, fWidth, fOffset); 171 171 172 #else 172 #else 173 new G4PVReplica(name, lv, mlv, fAxis, fNofDi 173 new G4PVReplica(name, lv, mlv, fAxis, fNofDivisions, fWidth, fOffset); 174 174 175 #endif 175 #endif 176 } 176 } 177 177 178 // private methods 178 // private methods 179 179 180 void G3Division::Exception(G4String where, G4S 180 void G3Division::Exception(G4String where, G4String what) 181 { 181 { 182 G4String err_message = "G3Division::" + wher 182 G4String err_message = "G3Division::" + where + " for " 183 + what + " is not imple 183 + what + " is not implemented"; 184 G4Exception("G3Division::Exception()", "G3to 184 G4Exception("G3Division::Exception()", "G3toG40004", 185 FatalException, err_message); 185 FatalException, err_message); 186 return; 186 return; 187 } 187 } 188 188 189 void G3Division::SetRangeAndAxis() 189 void G3Division::SetRangeAndAxis() 190 // set fHighRange, fLowRange, fAxis 190 // set fHighRange, fLowRange, fAxis 191 { 191 { 192 G4String shape = fMVTE->GetShape(); 192 G4String shape = fMVTE->GetShape(); 193 G4double *Rpar = fMVTE->GetRpar(); 193 G4double *Rpar = fMVTE->GetRpar(); 194 194 195 switch (fIAxis) { 195 switch (fIAxis) { 196 case 1: fAxis = kXAxis; 196 case 1: fAxis = kXAxis; 197 break; 197 break; 198 case 2: fAxis = kYAxis; 198 case 2: fAxis = kYAxis; 199 break; 199 break; 200 case 3: fAxis = kZAxis; 200 case 3: fAxis = kZAxis; 201 break; 201 break; 202 default: G4Exception("G3Division::SetRan 202 default: G4Exception("G3Division::SetRangeAndAxis()", "G3toG40005", 203 FatalException, "W 203 FatalException, "Wrong axis index definition!"); 204 } 204 } 205 205 206 if ( shape == "BOX" ) { 206 if ( shape == "BOX" ) { 207 fHighRange = Rpar[fIAxis-1]*cm; 207 fHighRange = Rpar[fIAxis-1]*cm; 208 fLowRange = -fHighRange; 208 fLowRange = -fHighRange; 209 } 209 } 210 else if ( shape == "TRD1" ) { 210 else if ( shape == "TRD1" ) { 211 if (fIAxis == 1){ 211 if (fIAxis == 1){ 212 fHighRange = std::max(Rpar[0]*cm, Rpar 212 fHighRange = std::max(Rpar[0]*cm, Rpar[1]*cm); 213 } 213 } 214 else if( fIAxis == 2) { 214 else if( fIAxis == 2) { 215 fHighRange = Rpar[2]*cm; 215 fHighRange = Rpar[2]*cm; 216 } 216 } 217 else if( fIAxis == 3) { 217 else if( fIAxis == 3) { 218 fHighRange = Rpar[3]*cm; 218 fHighRange = Rpar[3]*cm; 219 } 219 } 220 fLowRange = - fHighRange; 220 fLowRange = - fHighRange; 221 } 221 } 222 else if ( shape == "TRD2" ) { 222 else if ( shape == "TRD2" ) { 223 if (fIAxis == 1){ 223 if (fIAxis == 1){ 224 fHighRange = std::max(Rpar[0]*cm, Rpar 224 fHighRange = std::max(Rpar[0]*cm, Rpar[1]*cm); 225 } 225 } 226 else if( fIAxis == 2) { 226 else if( fIAxis == 2) { 227 fHighRange = std::max(Rpar[2]*cm, Rpar 227 fHighRange = std::max(Rpar[2]*cm, Rpar[3]*cm); 228 } 228 } 229 else if( fIAxis == 3) { 229 else if( fIAxis == 3) { 230 fHighRange = Rpar[4]*cm; 230 fHighRange = Rpar[4]*cm; 231 } 231 } 232 } 232 } 233 else if ( shape == "TRAP" ) { 233 else if ( shape == "TRAP" ) { 234 if ( fIAxis == 3 ) fHighRange = Rpar[0]* 234 if ( fIAxis == 3 ) fHighRange = Rpar[0]*cm; 235 else fHighRange = 0.; 235 else fHighRange = 0.; 236 fLowRange = -fHighRange; 236 fLowRange = -fHighRange; 237 } 237 } 238 else if ( shape == "TUBE" ) { 238 else if ( shape == "TUBE" ) { 239 if (fIAxis == 1){ 239 if (fIAxis == 1){ 240 fHighRange = Rpar[1]*cm; 240 fHighRange = Rpar[1]*cm; 241 fLowRange = Rpar[0]*cm; 241 fLowRange = Rpar[0]*cm; 242 fAxis = kRho; 242 fAxis = kRho; 243 } 243 } 244 else if( fIAxis == 2) { 244 else if( fIAxis == 2) { 245 fHighRange = 360.*deg; 245 fHighRange = 360.*deg; 246 fLowRange = 0.; 246 fLowRange = 0.; 247 fAxis = kPhi; 247 fAxis = kPhi; 248 } 248 } 249 else if( fIAxis == 3) { 249 else if( fIAxis == 3) { 250 fHighRange = Rpar[2]*cm; 250 fHighRange = Rpar[2]*cm; 251 fLowRange = -fHighRange; 251 fLowRange = -fHighRange; 252 } 252 } 253 } 253 } 254 else if ( shape == "TUBS" ) { 254 else if ( shape == "TUBS" ) { 255 if (fIAxis == 1){ 255 if (fIAxis == 1){ 256 fHighRange = Rpar[1]*cm; 256 fHighRange = Rpar[1]*cm; 257 fLowRange = Rpar[0]*cm; 257 fLowRange = Rpar[0]*cm; 258 fAxis = kRho; 258 fAxis = kRho; 259 } 259 } 260 else if( fIAxis == 2) { 260 else if( fIAxis == 2) { 261 261 262 fLowRange = Rpar[3]*deg; 262 fLowRange = Rpar[3]*deg; 263 fHighRange = Rpar[4]*deg - fLowRange; 263 fHighRange = Rpar[4]*deg - fLowRange; 264 if ( Rpar[4]*deg <= fLowRange )fHighRan 264 if ( Rpar[4]*deg <= fLowRange )fHighRange = fHighRange + 360.*deg; 265 fHighRange = fHighRange + fLowRange; 265 fHighRange = fHighRange + fLowRange; 266 fAxis = kPhi; 266 fAxis = kPhi; 267 } 267 } 268 else if( fIAxis == 3) { 268 else if( fIAxis == 3) { 269 fHighRange = Rpar[2]*cm; 269 fHighRange = Rpar[2]*cm; 270 fLowRange = -fHighRange; 270 fLowRange = -fHighRange; 271 } 271 } 272 } 272 } 273 else if ( shape == "CONE" ) { 273 else if ( shape == "CONE" ) { 274 if (fIAxis == 1){ 274 if (fIAxis == 1){ 275 fHighRange = std::max(Rpar[2]*cm,Rpar[ 275 fHighRange = std::max(Rpar[2]*cm,Rpar[4]*cm); 276 fLowRange = std::max(Rpar[1]*cm,Rpar[3 276 fLowRange = std::max(Rpar[1]*cm,Rpar[3]*cm); 277 fAxis = kRho; 277 fAxis = kRho; 278 } 278 } 279 else if( fIAxis == 2) { 279 else if( fIAxis == 2) { 280 280 281 fLowRange = 0.; 281 fLowRange = 0.; 282 fHighRange = 360.*deg; 282 fHighRange = 360.*deg; 283 fAxis = kPhi; 283 fAxis = kPhi; 284 } 284 } 285 else if( fIAxis == 3) { 285 else if( fIAxis == 3) { 286 fHighRange = Rpar[0]*cm; 286 fHighRange = Rpar[0]*cm; 287 fLowRange = -fHighRange; 287 fLowRange = -fHighRange; 288 } 288 } 289 } 289 } 290 else if ( shape == "CONS" ) { 290 else if ( shape == "CONS" ) { 291 if (fIAxis == 1){ 291 if (fIAxis == 1){ 292 fHighRange = std::max(Rpar[2]*cm,Rpar[ 292 fHighRange = std::max(Rpar[2]*cm,Rpar[4]*cm); 293 fLowRange = std::max(Rpar[1]*cm,Rpar[3 293 fLowRange = std::max(Rpar[1]*cm,Rpar[3]*cm); 294 fAxis = kRho; 294 fAxis = kRho; 295 } 295 } 296 else if( fIAxis == 2) { 296 else if( fIAxis == 2) { 297 297 298 fLowRange = Rpar[5]*deg; 298 fLowRange = Rpar[5]*deg; 299 fHighRange = Rpar[6]*deg - fLowRange; 299 fHighRange = Rpar[6]*deg - fLowRange; 300 if ( Rpar[6]*deg <= fLowRange )fHighRan 300 if ( Rpar[6]*deg <= fLowRange )fHighRange = fHighRange + 360.*deg; 301 fHighRange = fHighRange + fLowRange; 301 fHighRange = fHighRange + fLowRange; 302 fAxis = kPhi; 302 fAxis = kPhi; 303 } 303 } 304 else if( fIAxis == 3) { 304 else if( fIAxis == 3) { 305 fHighRange = Rpar[2]*cm; 305 fHighRange = Rpar[2]*cm; 306 fLowRange = -fHighRange; 306 fLowRange = -fHighRange; 307 } 307 } 308 } 308 } 309 else if ( shape == "SPHE" ) { 309 else if ( shape == "SPHE" ) { 310 if (fIAxis == 1){ 310 if (fIAxis == 1){ 311 fHighRange = Rpar[1]*cm; 311 fHighRange = Rpar[1]*cm; 312 fLowRange = Rpar[0]*cm; 312 fLowRange = Rpar[0]*cm; 313 fAxis = kRho; 313 fAxis = kRho; 314 } 314 } 315 else if( fIAxis == 2) { 315 else if( fIAxis == 2) { 316 fLowRange = std::min(Rpar[2]*deg,Rpar[3 316 fLowRange = std::min(Rpar[2]*deg,Rpar[3]*deg); 317 fHighRange = std::max(Rpar[2]*deg,Rpar[ 317 fHighRange = std::max(Rpar[2]*deg,Rpar[3]*deg); 318 fAxis = kPhi; 318 fAxis = kPhi; 319 } 319 } 320 else if( fIAxis == 3) { 320 else if( fIAxis == 3) { 321 fLowRange = std::min(Rpar[4]*deg,Rpar[5 321 fLowRange = std::min(Rpar[4]*deg,Rpar[5]*deg); 322 fHighRange = std::max(Rpar[4]*deg,Rpar[ 322 fHighRange = std::max(Rpar[4]*deg,Rpar[5]*deg); 323 fAxis = kPhi; // ?????? 323 fAxis = kPhi; // ?????? 324 } 324 } 325 } 325 } 326 else if ( shape == "PARA" ) { 326 else if ( shape == "PARA" ) { 327 fHighRange = Rpar[fIAxis-1]*cm; 327 fHighRange = Rpar[fIAxis-1]*cm; 328 fLowRange = -fHighRange; 328 fLowRange = -fHighRange; 329 } 329 } 330 else if ( shape == "PGON" ) { 330 else if ( shape == "PGON" ) { 331 G4int i; 331 G4int i; 332 G4int nz = G4int(Rpar[3]); 332 G4int nz = G4int(Rpar[3]); 333 333 334 G4double pPhi1 = Rpar[0]*deg; 334 G4double pPhi1 = Rpar[0]*deg; 335 G4double dPhi = Rpar[1]*deg; 335 G4double dPhi = Rpar[1]*deg; 336 336 337 G4double *DzArray = new G4double[nz]; 337 G4double *DzArray = new G4double[nz]; 338 G4double *Rmax = new G4double[nz]; 338 G4double *Rmax = new G4double[nz]; 339 G4double *Rmin = new G4double[nz]; 339 G4double *Rmin = new G4double[nz]; 340 G4double rangehi[3], rangelo[3]; 340 G4double rangehi[3], rangelo[3]; 341 rangehi[0] = -kInfinity ; 341 rangehi[0] = -kInfinity ; 342 rangelo[0] = kInfinity ; 342 rangelo[0] = kInfinity ; 343 rangehi[2] = -kInfinity ; 343 rangehi[2] = -kInfinity ; 344 rangelo[2] = kInfinity ; 344 rangelo[2] = kInfinity ; 345 345 346 for(i=0; i<nz; i++) 346 for(i=0; i<nz; i++) 347 { 347 { 348 G4int i4=3*i+4; 348 G4int i4=3*i+4; 349 G4int i5=i4+1; 349 G4int i5=i4+1; 350 G4int i6=i4+2; 350 G4int i6=i4+2; 351 351 352 DzArray[i] = Rpar[i4]*cm; 352 DzArray[i] = Rpar[i4]*cm; 353 Rmin[i] = Rpar[i5]*cm; 353 Rmin[i] = Rpar[i5]*cm; 354 Rmax[i] = Rpar[i6]*cm; 354 Rmax[i] = Rpar[i6]*cm; 355 rangelo[0] = std::min(rangelo[0], 355 rangelo[0] = std::min(rangelo[0], Rmin[i]); 356 rangehi[0] = std::max(rangehi[0], 356 rangehi[0] = std::max(rangehi[0], Rmax[i]); 357 rangelo[2] = std::min(rangelo[2], 357 rangelo[2] = std::min(rangelo[2], DzArray[i]); 358 rangehi[2] = std::max(rangehi[2], 358 rangehi[2] = std::max(rangehi[2], DzArray[i]); 359 } 359 } 360 for (i=0;i<nz;i++){ 360 for (i=0;i<nz;i++){ 361 assert(Rmin[i]>=0 && Rmax[i]>=Rmin 361 assert(Rmin[i]>=0 && Rmax[i]>=Rmin[i]); 362 } 362 } 363 rangehi[1] = pPhi1 + dPhi; 363 rangehi[1] = pPhi1 + dPhi; 364 rangelo[1] = pPhi1; 364 rangelo[1] = pPhi1; 365 fHighRange = rangehi[fIAxis-1]; 365 fHighRange = rangehi[fIAxis-1]; 366 fLowRange = rangelo[fIAxis-1]; 366 fLowRange = rangelo[fIAxis-1]; 367 if (fIAxis == 1)fAxis = kRho; 367 if (fIAxis == 1)fAxis = kRho; 368 else if (fIAxis == 2)fAxis = kPhi; 368 else if (fIAxis == 2)fAxis = kPhi; 369 else if (fIAxis == 3)fAxis = kZAxis; 369 else if (fIAxis == 3)fAxis = kZAxis; 370 370 371 delete [] DzArray; 371 delete [] DzArray; 372 delete [] Rmin; 372 delete [] Rmin; 373 delete [] Rmax; 373 delete [] Rmax; 374 374 375 } 375 } 376 else if ( shape == "PCON" ) { 376 else if ( shape == "PCON" ) { 377 377 378 G4int i; 378 G4int i; 379 G4double pPhi1 = Rpar[0]*deg; 379 G4double pPhi1 = Rpar[0]*deg; 380 G4double dPhi = Rpar[1]*deg; 380 G4double dPhi = Rpar[1]*deg; 381 G4int nz = G4int(Rpar[2]); 381 G4int nz = G4int(Rpar[2]); 382 382 383 G4double *DzArray = new G4double[nz]; 383 G4double *DzArray = new G4double[nz]; 384 G4double *Rmax = new G4double[nz]; 384 G4double *Rmax = new G4double[nz]; 385 G4double *Rmin = new G4double[nz]; 385 G4double *Rmin = new G4double[nz]; 386 G4double rangehi[3],rangelo[3]; 386 G4double rangehi[3],rangelo[3]; 387 387 388 rangehi[0] = -kInfinity ; 388 rangehi[0] = -kInfinity ; 389 rangelo[0] = kInfinity ; 389 rangelo[0] = kInfinity ; 390 rangehi[2] = -kInfinity ; 390 rangehi[2] = -kInfinity ; 391 rangelo[2] = kInfinity ; 391 rangelo[2] = kInfinity ; 392 392 393 for(i=0; i<nz; i++){ 393 for(i=0; i<nz; i++){ 394 G4int i4=3*i+3; 394 G4int i4=3*i+3; 395 G4int i5=i4+1; 395 G4int i5=i4+1; 396 G4int i6=i4+2; 396 G4int i6=i4+2; 397 397 398 DzArray[i] = Rpar[i4]*cm; 398 DzArray[i] = Rpar[i4]*cm; 399 Rmin[i] = Rpar[i5]*cm; 399 Rmin[i] = Rpar[i5]*cm; 400 Rmax[i] = Rpar[i6]*cm; 400 Rmax[i] = Rpar[i6]*cm; 401 rangelo[0] = std::min(rangelo[0], 401 rangelo[0] = std::min(rangelo[0], Rmin[i]); 402 rangehi[0] = std::max(rangehi[0], 402 rangehi[0] = std::max(rangehi[0], Rmax[i]); 403 rangelo[2] = std::min(rangelo[2], 403 rangelo[2] = std::min(rangelo[2], DzArray[i]); 404 rangehi[2] = std::max(rangehi[2], 404 rangehi[2] = std::max(rangehi[2], DzArray[i]); 405 } 405 } 406 for (i=0;i<nz;i++){ 406 for (i=0;i<nz;i++){ 407 assert(Rmin[i]>=0 && Rmax[i]>=Rmin 407 assert(Rmin[i]>=0 && Rmax[i]>=Rmin[i]); 408 } 408 } 409 rangehi[1] = pPhi1 + dPhi; 409 rangehi[1] = pPhi1 + dPhi; 410 rangelo[1] = pPhi1; 410 rangelo[1] = pPhi1; 411 fHighRange = rangehi[fIAxis-1]; 411 fHighRange = rangehi[fIAxis-1]; 412 fLowRange = rangelo[fIAxis-1]; 412 fLowRange = rangelo[fIAxis-1]; 413 if (fIAxis == 1)fAxis = kRho; 413 if (fIAxis == 1)fAxis = kRho; 414 else if (fIAxis == 2)fAxis = kPhi; 414 else if (fIAxis == 2)fAxis = kPhi; 415 else if (fIAxis == 3)fAxis = kZAxis; 415 else if (fIAxis == 3)fAxis = kZAxis; 416 416 417 417 418 delete [] DzArray; 418 delete [] DzArray; 419 delete [] Rmin; 419 delete [] Rmin; 420 delete [] Rmax; 420 delete [] Rmax; 421 } 421 } 422 else if ( shape == "ELTU" || shape == "HY 422 else if ( shape == "ELTU" || shape == "HYPE" || shape == "GTRA" || 423 shape == "CTUB") { 423 shape == "CTUB") { 424 Exception("SetRangeAndAxis", shape); 424 Exception("SetRangeAndAxis", shape); 425 } 425 } 426 else { 426 else { 427 Exception("SetRangeAndAxis", "Unknown s 427 Exception("SetRangeAndAxis", "Unknown shape" + shape); 428 } 428 } 429 429 430 // verbose 430 // verbose 431 #ifdef G3G4DEBUG 431 #ifdef G3G4DEBUG 432 G4cout << "Shape " << shape << " SetRang 432 G4cout << "Shape " << shape << " SetRangeAndAxis: " 433 << fLowRange << " " << fHighRange << " 433 << fLowRange << " " << fHighRange << " " << fAxis << G4endl; 434 #endif 434 #endif 435 } 435 } 436 436 437 G3VolTableEntry* G3Division::CreateEnvelope(G4 437 G3VolTableEntry* G3Division::CreateEnvelope(G4String shape, G4double hi, 438 G4double lo, G4 438 G4double lo, G4double par[], G4int npar) 439 // create new VTE with G3Pos corresponding to 439 // create new VTE with G3Pos corresponding to the 440 // envelope of divided volume 440 // envelope of divided volume 441 { 441 { 442 // verbose 442 // verbose 443 // G4cout << " G3Division::CreateEnvelope 443 // G4cout << " G3Division::CreateEnvelope " << "fIAaxis= " << fIAxis 444 // << " hi= " << hi 444 // << " hi= " << hi 445 // << " lo= " << lo 445 // << " lo= " << lo 446 // << G4endl; 446 // << G4endl; 447 447 448 G4double *Rpar = new G4double[npar+2]; 448 G4double *Rpar = new G4double[npar+2]; 449 for (G4int i=0; i<npar; ++i){ Rpar[i] = pa 449 for (G4int i=0; i<npar; ++i){ Rpar[i] = par[i];} 450 G4double pos[3] = {0.,0.,0.}; 450 G4double pos[3] = {0.,0.,0.}; 451 451 452 if ( shape == "BOX" ) { 452 if ( shape == "BOX" ) { 453 Rpar[fIAxis-1] = (hi - lo)/2./cm; 453 Rpar[fIAxis-1] = (hi - lo)/2./cm; 454 pos [fIAxis-1] = (hi + lo)/2.; 454 pos [fIAxis-1] = (hi + lo)/2.; 455 } 455 } 456 else if ( shape == "TRD1" ) { 456 else if ( shape == "TRD1" ) { 457 if ( fIAxis == 1 || fIAxis == 2 ) { 457 if ( fIAxis == 1 || fIAxis == 2 ) { 458 Exception("CreateEnvelope","TRD1-x,y") 458 Exception("CreateEnvelope","TRD1-x,y"); 459 } 459 } 460 else if ( fIAxis == 3 ) { 460 else if ( fIAxis == 3 ) { 461 // x = x1 + (c-z1)(x2 -x1)/(z2-z1) 461 // x = x1 + (c-z1)(x2 -x1)/(z2-z1) 462 G4double tn, x1, z1; 462 G4double tn, x1, z1; 463 tn = (Rpar[1] - Rpar[0])/(2.* Rpar[3]) 463 tn = (Rpar[1] - Rpar[0])/(2.* Rpar[3]); 464 x1 = Rpar[0]; z1 = -Rpar[3]; 464 x1 = Rpar[0]; z1 = -Rpar[3]; 465 Rpar[0] = x1 + tn * (lo/cm - z1); 465 Rpar[0] = x1 + tn * (lo/cm - z1); 466 Rpar[1] = x1 + tn * (hi/cm - z1); 466 Rpar[1] = x1 + tn * (hi/cm - z1); 467 Rpar[3] = (hi - lo)/2./cm; 467 Rpar[3] = (hi - lo)/2./cm; 468 pos[2] = (hi + lo)/2.; 468 pos[2] = (hi + lo)/2.; 469 } 469 } 470 } 470 } 471 else if ( shape == "TRD2" ) { 471 else if ( shape == "TRD2" ) { 472 if ( fIAxis == 1 || fIAxis == 2) { 472 if ( fIAxis == 1 || fIAxis == 2) { 473 Exception("CreateEnvelope","TRD2-x,y") 473 Exception("CreateEnvelope","TRD2-x,y"); 474 } 474 } 475 else if ( fIAxis == 3 ) { 475 else if ( fIAxis == 3 ) { 476 // x = x1 + (c-z1)(x2 -x1)/(z2-z1) 476 // x = x1 + (c-z1)(x2 -x1)/(z2-z1) 477 // y = y1 + (c-z1)(y2 -y1)/(z2-z1) 477 // y = y1 + (c-z1)(y2 -y1)/(z2-z1) 478 G4double tn1, tn2, x1, y1, z1; 478 G4double tn1, tn2, x1, y1, z1; 479 tn1 = (Rpar[1] - Rpar[0])/(2.* Rpar[4] 479 tn1 = (Rpar[1] - Rpar[0])/(2.* Rpar[4]); 480 tn2 = (Rpar[3] - Rpar[2])/(2.* Rpar[4] 480 tn2 = (Rpar[3] - Rpar[2])/(2.* Rpar[4]); 481 x1 = Rpar[0]; y1 = Rpar[2]; z1 = -Rpar 481 x1 = Rpar[0]; y1 = Rpar[2]; z1 = -Rpar[3]; 482 Rpar[0] = x1 + tn1 * (lo/cm - z1); 482 Rpar[0] = x1 + tn1 * (lo/cm - z1); 483 Rpar[1] = x1 + tn1 * (hi/cm - z1); 483 Rpar[1] = x1 + tn1 * (hi/cm - z1); 484 Rpar[2] = y1 + tn2 * (lo/cm - z1); 484 Rpar[2] = y1 + tn2 * (lo/cm - z1); 485 Rpar[3] = y1 + tn2 * (hi/cm - z1); 485 Rpar[3] = y1 + tn2 * (hi/cm - z1); 486 Rpar[4] = (hi - lo)/2./cm; 486 Rpar[4] = (hi - lo)/2./cm; 487 pos[2] = (hi + lo)/2.; 487 pos[2] = (hi + lo)/2.; 488 } 488 } 489 } 489 } 490 else if ( shape == "TRAP" ) { 490 else if ( shape == "TRAP" ) { 491 Exception("CreateEnvelope","TRAP-x,y,z") 491 Exception("CreateEnvelope","TRAP-x,y,z"); 492 } 492 } 493 else if ( shape == "TUBE" ) { 493 else if ( shape == "TUBE" ) { 494 if ( fIAxis == 1 ) { 494 if ( fIAxis == 1 ) { 495 Rpar[0] = lo/cm; 495 Rpar[0] = lo/cm; 496 Rpar[1] = hi/cm; 496 Rpar[1] = hi/cm; 497 } 497 } 498 else if ( fIAxis == 2 ) { 498 else if ( fIAxis == 2 ) { 499 Rpar[3] = lo/deg; 499 Rpar[3] = lo/deg; 500 Rpar[4] = hi/deg; 500 Rpar[4] = hi/deg; 501 npar = npar + 2; 501 npar = npar + 2; 502 shape = "TUBS"; 502 shape = "TUBS"; 503 } 503 } 504 else if ( fIAxis == 3 ) { 504 else if ( fIAxis == 3 ) { 505 Rpar[2] = (hi - lo)/2./cm; 505 Rpar[2] = (hi - lo)/2./cm; 506 pos [2] = (hi + lo)/2.; 506 pos [2] = (hi + lo)/2.; 507 } 507 } 508 } 508 } 509 else if ( shape == "TUBS" ) { 509 else if ( shape == "TUBS" ) { 510 if ( fIAxis == 1 ) { 510 if ( fIAxis == 1 ) { 511 Rpar[0] = lo/cm; 511 Rpar[0] = lo/cm; 512 Rpar[1] = hi/cm; 512 Rpar[1] = hi/cm; 513 } 513 } 514 else if ( fIAxis == 2 ) { 514 else if ( fIAxis == 2 ) { 515 Rpar[3] = lo/deg; 515 Rpar[3] = lo/deg; 516 Rpar[4] = hi/deg; 516 Rpar[4] = hi/deg; 517 } 517 } 518 else if ( fIAxis == 3 ) { 518 else if ( fIAxis == 3 ) { 519 Rpar[2] = (hi - lo)/2./cm; 519 Rpar[2] = (hi - lo)/2./cm; 520 pos [2] = (hi + lo)/2.; 520 pos [2] = (hi + lo)/2.; 521 } 521 } 522 } 522 } 523 else if ( shape == "CONE" ) { 523 else if ( shape == "CONE" ) { 524 if ( fIAxis == 1) { 524 if ( fIAxis == 1) { 525 Exception("CreateEnvelope","CONE-x,z") 525 Exception("CreateEnvelope","CONE-x,z"); 526 } 526 } 527 else if ( fIAxis == 2 ) { 527 else if ( fIAxis == 2 ) { 528 Rpar[5] = lo/deg; 528 Rpar[5] = lo/deg; 529 Rpar[6] = hi/deg; 529 Rpar[6] = hi/deg; 530 npar = npar + 2; 530 npar = npar + 2; 531 shape = "CONS"; 531 shape = "CONS"; 532 } 532 } 533 else if ( fIAxis == 3 ) { 533 else if ( fIAxis == 3 ) { 534 G4double tn1, tn2, rmin, rmax, z1; 534 G4double tn1, tn2, rmin, rmax, z1; 535 tn1 = (Rpar[3] - Rpar[1])/(2.* Rpar[0] 535 tn1 = (Rpar[3] - Rpar[1])/(2.* Rpar[0]); 536 tn2 = (Rpar[4] - Rpar[2])/(2.* Rpar[0] 536 tn2 = (Rpar[4] - Rpar[2])/(2.* Rpar[0]); 537 rmin = Rpar[1]; rmax = Rpar[2]; z1 = - 537 rmin = Rpar[1]; rmax = Rpar[2]; z1 = -Rpar[0]; 538 Rpar[1] = rmin + tn1 * (lo/cm - z1); 538 Rpar[1] = rmin + tn1 * (lo/cm - z1); 539 Rpar[3] = rmin + tn1 * (hi/cm - z1); 539 Rpar[3] = rmin + tn1 * (hi/cm - z1); 540 Rpar[2] = rmax + tn2 * (lo/cm - z1); 540 Rpar[2] = rmax + tn2 * (lo/cm - z1); 541 Rpar[4] = rmax + tn2 * (hi/cm - z1); 541 Rpar[4] = rmax + tn2 * (hi/cm - z1); 542 Rpar[0] = (hi - lo)/2./cm; 542 Rpar[0] = (hi - lo)/2./cm; 543 pos[2] = (hi + lo)/2.; 543 pos[2] = (hi + lo)/2.; 544 } 544 } 545 } 545 } 546 else if ( shape == "CONS" ) { 546 else if ( shape == "CONS" ) { 547 if ( fIAxis == 1 ) { 547 if ( fIAxis == 1 ) { 548 Exception("CreateEnvelope","CONS-x"); 548 Exception("CreateEnvelope","CONS-x"); 549 } 549 } 550 else if ( fIAxis == 2 ) { 550 else if ( fIAxis == 2 ) { 551 Rpar[5] = lo/deg; 551 Rpar[5] = lo/deg; 552 Rpar[6] = hi/deg; 552 Rpar[6] = hi/deg; 553 } 553 } 554 else if ( fIAxis == 3 ) { 554 else if ( fIAxis == 3 ) { 555 G4double tn1, tn2, rmin, rmax, z1; 555 G4double tn1, tn2, rmin, rmax, z1; 556 tn1 = (Rpar[3] - Rpar[1])/(2.* Rpar[0] 556 tn1 = (Rpar[3] - Rpar[1])/(2.* Rpar[0]); 557 tn2 = (Rpar[4] - Rpar[2])/(2.* Rpar[0] 557 tn2 = (Rpar[4] - Rpar[2])/(2.* Rpar[0]); 558 rmin = Rpar[1]; rmax = Rpar[2]; z1 = - 558 rmin = Rpar[1]; rmax = Rpar[2]; z1 = -Rpar[0]; 559 Rpar[1] = rmin + tn1 * (lo/cm - z1); 559 Rpar[1] = rmin + tn1 * (lo/cm - z1); 560 Rpar[3] = rmin + tn1 * (hi/cm - z1); 560 Rpar[3] = rmin + tn1 * (hi/cm - z1); 561 Rpar[2] = rmax + tn2 * (lo/cm - z1); 561 Rpar[2] = rmax + tn2 * (lo/cm - z1); 562 Rpar[4] = rmax + tn2 * (hi/cm - z1); 562 Rpar[4] = rmax + tn2 * (hi/cm - z1); 563 Rpar[0] = (hi - lo)/2./cm; 563 Rpar[0] = (hi - lo)/2./cm; 564 pos[2] = (hi + lo)/2.; 564 pos[2] = (hi + lo)/2.; 565 } 565 } 566 } 566 } 567 else if ( shape == "SPHE" ) { 567 else if ( shape == "SPHE" ) { 568 Exception("CreateEnvelope","SPHE-x,y,z") 568 Exception("CreateEnvelope","SPHE-x,y,z"); 569 } 569 } 570 else if ( shape == "PARA" ) { 570 else if ( shape == "PARA" ) { 571 Exception("CreateEnvelope","PARA-x,y,z") 571 Exception("CreateEnvelope","PARA-x,y,z"); 572 } 572 } 573 else if ( shape == "PGON" ) { 573 else if ( shape == "PGON" ) { 574 if ( fIAxis == 2) { 574 if ( fIAxis == 2) { 575 Rpar[0] = lo/deg; 575 Rpar[0] = lo/deg; 576 Rpar[1] = hi/deg; 576 Rpar[1] = hi/deg; 577 // rotm = ??? 577 // rotm = ??? 578 } 578 } 579 else { 579 else { 580 Exception("CreateEnvelope","PGON-x,z") 580 Exception("CreateEnvelope","PGON-x,z"); 581 } 581 } 582 } 582 } 583 else if ( shape == "PCON" ) { 583 else if ( shape == "PCON" ) { 584 if ( fIAxis == 2) { 584 if ( fIAxis == 2) { 585 Rpar[0] = lo/deg; 585 Rpar[0] = lo/deg; 586 Rpar[1] = hi/deg; 586 Rpar[1] = hi/deg; 587 // rotm = ??? 587 // rotm = ??? 588 } 588 } 589 else { 589 else { 590 Exception("CreateEnvelope","PCON-x,z") 590 Exception("CreateEnvelope","PCON-x,z"); 591 } 591 } 592 } 592 } 593 else { 593 else { 594 Exception("CreateEnvelope", "Unknown sh 594 Exception("CreateEnvelope", "Unknown shape" + shape); 595 } 595 } 596 596 597 // create new VTE corresponding to envelop 597 // create new VTE corresponding to envelope 598 G4String envName = fVTE->GetName() + "_ENV 598 G4String envName = fVTE->GetName() + "_ENV"; 599 G3VolTableEntry* envVTE 599 G3VolTableEntry* envVTE 600 = G4CreateVTE(envName, shape, fNmed, Rpa 600 = G4CreateVTE(envName, shape, fNmed, Rpar, npar); 601 601 602 // create a G3Pos object and add it to env 602 // create a G3Pos object and add it to envVTE 603 G4String motherName = fMVTE->GetMasterClon 603 G4String motherName = fMVTE->GetMasterClone()->GetName(); 604 G4ThreeVector* offset = new G4ThreeVector( 604 G4ThreeVector* offset = new G4ThreeVector(pos[0],pos[1],pos[2]); 605 G4String only = "ONLY"; 605 G4String only = "ONLY"; 606 G3Pos* aG3Pos = new G3Pos(motherName, 1, o 606 G3Pos* aG3Pos = new G3Pos(motherName, 1, offset, 0, only); 607 envVTE->AddG3Pos(aG3Pos); 607 envVTE->AddG3Pos(aG3Pos); 608 608 609 delete [] Rpar; 609 delete [] Rpar; 610 610 611 return envVTE; 611 return envVTE; 612 } 612 } 613 613 614 void G3Division::CreateSolid(G4String shape, G 614 void G3Division::CreateSolid(G4String shape, G4double par[], G4int npar) 615 // create the solid corresponding to divided v 615 // create the solid corresponding to divided volume 616 // and set the fOffset for replica 616 // and set the fOffset for replica 617 { 617 { 618 G4double *Rpar = new G4double[npar+2]; 618 G4double *Rpar = new G4double[npar+2]; 619 for (G4int i=0; i<npar; ++i){ Rpar[i] = pa 619 for (G4int i=0; i<npar; ++i){ Rpar[i] = par[i];} 620 620 621 // verbose 621 // verbose 622 // G4cout << "G3Division::CreateSolid volu 622 // G4cout << "G3Division::CreateSolid volume before: " 623 // << fVTE->GetName() << " " << sha 623 // << fVTE->GetName() << " " << shape << G4endl; 624 // G4cout << " npar,Rpar: " << npar; 624 // G4cout << " npar,Rpar: " << npar; 625 // for (G4int ii = 0; ii < npar; ++ii) G4c 625 // for (G4int ii = 0; ii < npar; ++ii) G4cout << " " << Rpar[ii]; 626 // G4cout << G4endl; 626 // G4cout << G4endl; 627 627 628 if ( shape == "BOX" ) { 628 if ( shape == "BOX" ) { 629 if ( fIAxis == 1 ) Rpar[0] = fWidth 629 if ( fIAxis == 1 ) Rpar[0] = fWidth/2./cm; 630 else if ( fIAxis == 2 ) Rpar[1] = fWidth 630 else if ( fIAxis == 2 ) Rpar[1] = fWidth/2./cm; 631 else if ( fIAxis == 3 ) Rpar[2] = fWidth 631 else if ( fIAxis == 3 ) Rpar[2] = fWidth/2./cm; 632 } 632 } 633 else if ( shape == "TRD1" ) { 633 else if ( shape == "TRD1" ) { 634 if ( fIAxis == 1 || fIAxis == 2 ) { 634 if ( fIAxis == 1 || fIAxis == 2 ) { 635 Exception("CreateSolid", "TRD1-x,y"); 635 Exception("CreateSolid", "TRD1-x,y"); 636 } 636 } 637 else if ( fIAxis == 3 ) { 637 else if ( fIAxis == 3 ) { 638 Rpar[3] = fWidth/2./cm; 638 Rpar[3] = fWidth/2./cm; 639 } 639 } 640 } 640 } 641 else if ( shape == "TRD2" ) { 641 else if ( shape == "TRD2" ) { 642 if ( fIAxis == 1 || fIAxis == 2 ) { 642 if ( fIAxis == 1 || fIAxis == 2 ) { 643 Exception("CreateSolid", "TRD2-x,y"); 643 Exception("CreateSolid", "TRD2-x,y"); 644 } 644 } 645 else if ( fIAxis == 3 ) { 645 else if ( fIAxis == 3 ) { 646 Rpar[4] = fWidth/2./cm; 646 Rpar[4] = fWidth/2./cm; 647 } 647 } 648 } 648 } 649 else if ( shape == "TRAP" ) { 649 else if ( shape == "TRAP" ) { 650 if ( fIAxis == 1 || fIAxis == 2) { 650 if ( fIAxis == 1 || fIAxis == 2) { 651 Exception("CreateSolid", "TRAP-x,y"); 651 Exception("CreateSolid", "TRAP-x,y"); 652 } 652 } 653 else if ( fIAxis == 3 ) { 653 else if ( fIAxis == 3 ) { 654 Rpar[0] = fWidth/2./cm; 654 Rpar[0] = fWidth/2./cm; 655 } 655 } 656 } 656 } 657 else if ( shape == "TUBE" ) { 657 else if ( shape == "TUBE" ) { 658 if ( fIAxis == 1 ) { 658 if ( fIAxis == 1 ) { 659 Rpar[1] = Rpar[0] + fWidth/cm; 659 Rpar[1] = Rpar[0] + fWidth/cm; 660 fOffset = Rpar[0]*cm; 660 fOffset = Rpar[0]*cm; 661 } 661 } 662 else if ( fIAxis == 2 ) { 662 else if ( fIAxis == 2 ) { 663 Rpar[3] = 0.; 663 Rpar[3] = 0.; 664 Rpar[4] = fWidth/deg; 664 Rpar[4] = fWidth/deg; 665 shape = "TUBS"; 665 shape = "TUBS"; 666 npar = npar + 2; 666 npar = npar + 2; 667 } 667 } 668 else if ( fIAxis == 3 ) { 668 else if ( fIAxis == 3 ) { 669 Rpar[2] = fWidth/2./cm; 669 Rpar[2] = fWidth/2./cm; 670 } 670 } 671 } 671 } 672 else if ( shape == "TUBS" ) { 672 else if ( shape == "TUBS" ) { 673 if ( fIAxis == 1 ) { 673 if ( fIAxis == 1 ) { 674 Rpar[1] = Rpar[0] + fWidth/cm; 674 Rpar[1] = Rpar[0] + fWidth/cm; 675 fOffset = Rpar[0]*cm; 675 fOffset = Rpar[0]*cm; 676 } 676 } 677 else if ( fIAxis == 2 ) { 677 else if ( fIAxis == 2 ) { 678 fOffset = Rpar[3]*deg; 678 fOffset = Rpar[3]*deg; 679 Rpar[3] = 0.; 679 Rpar[3] = 0.; 680 Rpar[4] = fWidth/deg; 680 Rpar[4] = fWidth/deg; 681 } 681 } 682 else if ( fIAxis == 3 ) { 682 else if ( fIAxis == 3 ) { 683 Rpar[2] = fWidth/2./cm; 683 Rpar[2] = fWidth/2./cm; 684 } 684 } 685 } 685 } 686 else if ( shape == "CONE" ) { 686 else if ( shape == "CONE" ) { 687 if ( fIAxis == 1 ) { 687 if ( fIAxis == 1 ) { 688 Exception("CreateSolid", "CONE-x"); 688 Exception("CreateSolid", "CONE-x"); 689 } 689 } 690 else if ( fIAxis == 2 ) { 690 else if ( fIAxis == 2 ) { 691 Rpar[5] = 0.; 691 Rpar[5] = 0.; 692 Rpar[6] = fWidth/deg; 692 Rpar[6] = fWidth/deg; 693 shape = "CONS"; 693 shape = "CONS"; 694 npar = npar + 2; 694 npar = npar + 2; 695 } 695 } 696 else if ( fIAxis == 3 ) { 696 else if ( fIAxis == 3 ) { 697 Rpar[0] = fWidth/2./cm; 697 Rpar[0] = fWidth/2./cm; 698 } 698 } 699 } 699 } 700 else if ( shape == "CONS" ) { 700 else if ( shape == "CONS" ) { 701 if ( fIAxis == 1 ) { 701 if ( fIAxis == 1 ) { 702 Exception("CreateSolid", "CONS-x"); 702 Exception("CreateSolid", "CONS-x"); 703 } 703 } 704 else if ( fIAxis == 2 ) { 704 else if ( fIAxis == 2 ) { 705 fOffset = Rpar[5]*deg; 705 fOffset = Rpar[5]*deg; 706 Rpar[5] = 0.; 706 Rpar[5] = 0.; 707 Rpar[6] = fWidth/deg; 707 Rpar[6] = fWidth/deg; 708 } 708 } 709 else if ( fIAxis == 3 ) { 709 else if ( fIAxis == 3 ) { 710 Rpar[0] = fWidth/2./cm; 710 Rpar[0] = fWidth/2./cm; 711 } 711 } 712 } 712 } 713 else if (shape == "PARA") { 713 else if (shape == "PARA") { 714 if ( fIAxis == 1 ) { 714 if ( fIAxis == 1 ) { 715 Rpar[0] = fWidth/2./cm; 715 Rpar[0] = fWidth/2./cm; 716 } 716 } 717 else if ( Rpar[4] == 0. && Rpar[5] == 0. 717 else if ( Rpar[4] == 0. && Rpar[5] == 0. ) { 718 // only special case for axis 2,3 is 718 // only special case for axis 2,3 is supported 719 if ( fIAxis == 2 ) { 719 if ( fIAxis == 2 ) { 720 Rpar[1] = fWidth/2./cm; 720 Rpar[1] = fWidth/2./cm; 721 } 721 } 722 else if ( fIAxis == 3) { 722 else if ( fIAxis == 3) { 723 Rpar[2] = fWidth/2./cm; 723 Rpar[2] = fWidth/2./cm; 724 } 724 } 725 } 725 } 726 else 726 else 727 Exception("CreateSolid", shape); 727 Exception("CreateSolid", shape); 728 } 728 } 729 else if (shape == "SPHE") { 729 else if (shape == "SPHE") { 730 Exception("CreateSolid", shape); 730 Exception("CreateSolid", shape); 731 } 731 } 732 else if ( shape == "PGON" ) { 732 else if ( shape == "PGON" ) { 733 if ( fIAxis == 2 ) { 733 if ( fIAxis == 2 ) { 734 fOffset = Rpar[0]*deg; 734 fOffset = Rpar[0]*deg; 735 Rpar[0] = 0.; 735 Rpar[0] = 0.; 736 Rpar[1] = fWidth/deg; 736 Rpar[1] = fWidth/deg; 737 Rpar[2] = 1.; 737 Rpar[2] = 1.; 738 } 738 } 739 else 739 else 740 Exception("CreateSolid", shape); 740 Exception("CreateSolid", shape); 741 } 741 } 742 else if ( shape == "PCON" ) { 742 else if ( shape == "PCON" ) { 743 if ( fIAxis == 2 ) { 743 if ( fIAxis == 2 ) { 744 fOffset = Rpar[0]*deg; 744 fOffset = Rpar[0]*deg; 745 Rpar[0] = 0.; 745 Rpar[0] = 0.; 746 Rpar[1] = fWidth/deg; 746 Rpar[1] = fWidth/deg; 747 } 747 } 748 else { 748 else { 749 Exception("CreateSolid", shape); 749 Exception("CreateSolid", shape); 750 } 750 } 751 } 751 } 752 else { 752 else { 753 Exception("CreateSolid", "Unknown shape 753 Exception("CreateSolid", "Unknown shape" + shape); 754 } 754 } 755 755 756 // create solid and set it to fVTE 756 // create solid and set it to fVTE 757 G4bool hasNegPars; 757 G4bool hasNegPars; 758 G4bool deferred; 758 G4bool deferred; 759 G4bool okAxis[3]; 759 G4bool okAxis[3]; 760 G4VSolid* solid 760 G4VSolid* solid 761 = G3toG4MakeSolid(fVTE->GetName(), shape, 761 = G3toG4MakeSolid(fVTE->GetName(), shape, Rpar, npar, hasNegPars, deferred, okAxis); 762 762 763 if (hasNegPars) { 763 if (hasNegPars) { 764 G4String err_message = "CreateSolid VTE 764 G4String err_message = "CreateSolid VTE " + fVTE->GetName() 765 + " has negative p 765 + " has negative parameters."; 766 G4Exception("G3Division::CreateSolid()" 766 G4Exception("G3Division::CreateSolid()", "G3toG40006", 767 FatalException, err_message 767 FatalException, err_message); 768 return; 768 return; 769 } 769 } 770 770 771 // update vte 771 // update vte 772 fVTE->SetSolid(solid); 772 fVTE->SetSolid(solid); 773 fVTE->SetNRpar(npar, Rpar); 773 fVTE->SetNRpar(npar, Rpar); 774 fVTE->SetHasNegPars(hasNegPars); 774 fVTE->SetHasNegPars(hasNegPars); 775 775 776 // verbose 776 // verbose 777 // G4cout << "G3Division::CreateSolid volu 777 // G4cout << "G3Division::CreateSolid volume after: " 778 // << fVTE->GetName() << " " << sha 778 // << fVTE->GetName() << " " << shape << G4endl; 779 // G4cout << " npar,Rpar: " << npar; 779 // G4cout << " npar,Rpar: " << npar; 780 // for (G4int iii = 0; iii < npar; ++iii) 780 // for (G4int iii = 0; iii < npar; ++iii) G4cout << " " << Rpar[iii]; 781 // G4cout << G4endl; 781 // G4cout << G4endl; 782 delete [] Rpar; 782 delete [] Rpar; 783 } 783 } 784 784 785 785 786 G3VolTableEntry* G3Division::Dvn() 786 G3VolTableEntry* G3Division::Dvn() 787 { 787 { 788 // no envelope need to be created 788 // no envelope need to be created 789 789 790 // get parameters from mother 790 // get parameters from mother 791 G4String shape = fMVTE->GetShape(); 791 G4String shape = fMVTE->GetShape(); 792 G4double* Rpar = fMVTE->GetRpar(); 792 G4double* Rpar = fMVTE->GetRpar(); 793 G4int npar = fMVTE->GetNpar(); 793 G4int npar = fMVTE->GetNpar(); 794 794 795 // set width for replica and create solid 795 // set width for replica and create solid 796 fWidth = (fHighRange - fLowRange)/fNofDivisi 796 fWidth = (fHighRange - fLowRange)/fNofDivisions; 797 CreateSolid(shape, Rpar, npar); 797 CreateSolid(shape, Rpar, npar); 798 798 799 return 0; 799 return 0; 800 } 800 } 801 801 802 G3VolTableEntry* G3Division::Dvn2() 802 G3VolTableEntry* G3Division::Dvn2() 803 { 803 { 804 // to be defined as const of this class 804 // to be defined as const of this class 805 G4double Rmin = 0.0001*cm; 805 G4double Rmin = 0.0001*cm; 806 806 807 G4String shape = fMVTE->GetShape(); 807 G4String shape = fMVTE->GetShape(); 808 G4double* Rpar = fMVTE->GetRpar(); 808 G4double* Rpar = fMVTE->GetRpar(); 809 G4int npar = fMVTE->GetNpar(); 809 G4int npar = fMVTE->GetNpar(); 810 810 811 G4double c0 = fC0; 811 G4double c0 = fC0; 812 if (fAxis == kPhi) c0 = c0*deg; 812 if (fAxis == kPhi) c0 = c0*deg; 813 else c0 = c0*cm; 813 else c0 = c0*cm; 814 814 815 // create envelope (if needed) 815 // create envelope (if needed) 816 G3VolTableEntry* envVTE = 0; 816 G3VolTableEntry* envVTE = 0; 817 if( std::abs(c0 - fLowRange) > Rmin) { 817 if( std::abs(c0 - fLowRange) > Rmin) { 818 envVTE = CreateEnvelope(shape, fHighRange, 818 envVTE = CreateEnvelope(shape, fHighRange, c0, Rpar, npar); 819 Rpar = envVTE->GetRpar(); 819 Rpar = envVTE->GetRpar(); 820 npar = envVTE->GetNpar(); 820 npar = envVTE->GetNpar(); 821 } 821 } 822 822 823 // set width for replica and create solid 823 // set width for replica and create solid 824 fWidth = (fHighRange - c0)/fNofDivisions; 824 fWidth = (fHighRange - c0)/fNofDivisions; 825 CreateSolid(shape, Rpar, npar); 825 CreateSolid(shape, Rpar, npar); 826 826 827 return envVTE; 827 return envVTE; 828 } 828 } 829 829 830 G3VolTableEntry* G3Division::Dvt() 830 G3VolTableEntry* G3Division::Dvt() 831 { 831 { 832 // to be defined as const of this class 832 // to be defined as const of this class 833 G4double Rmin = 0.0001*cm; 833 G4double Rmin = 0.0001*cm; 834 834 835 // get parameters from mother 835 // get parameters from mother 836 G4String shape = fMVTE->GetShape(); 836 G4String shape = fMVTE->GetShape(); 837 G4double* Rpar = fMVTE->GetRpar(); 837 G4double* Rpar = fMVTE->GetRpar(); 838 G4int npar = fMVTE->GetNpar(); 838 G4int npar = fMVTE->GetNpar(); 839 839 840 // calculate the number of divisions 840 // calculate the number of divisions 841 G4int ndvmx = fNofDivisions; 841 G4int ndvmx = fNofDivisions; 842 G4double step = fStep; 842 G4double step = fStep; 843 843 844 if (fAxis == kPhi) step = step*deg; 844 if (fAxis == kPhi) step = step*deg; 845 else step = step*cm; 845 else step = step*cm; 846 846 847 G4int ndiv = G4int((fHighRange - fLowRange + 847 G4int ndiv = G4int((fHighRange - fLowRange + Rmin)/step); 848 // to be added warning 848 // to be added warning 849 if (ndvmx > 255) ndvmx = 255; 849 if (ndvmx > 255) ndvmx = 255; 850 if (ndiv > ndvmx && ndvmx > 0 ) ndiv = ndvmx 850 if (ndiv > ndvmx && ndvmx > 0 ) ndiv = ndvmx; 851 851 852 // create envVTE (if needed) 852 // create envVTE (if needed) 853 G3VolTableEntry* envVTE = 0; 853 G3VolTableEntry* envVTE = 0; 854 G4double delta = std::abs((fHighRange - fLow 854 G4double delta = std::abs((fHighRange - fLowRange) - ndiv*step); 855 if (delta > Rmin) { 855 if (delta > Rmin) { 856 envVTE 856 envVTE 857 = CreateEnvelope(shape, fHighRange-delt 857 = CreateEnvelope(shape, fHighRange-delta/2., fLowRange+delta/2., 858 Rpar, npar); 858 Rpar, npar); 859 Rpar = envVTE->GetRpar(); 859 Rpar = envVTE->GetRpar(); 860 npar = envVTE->GetNpar(); 860 npar = envVTE->GetNpar(); 861 } 861 } 862 862 863 // set width for replica and create solid 863 // set width for replica and create solid 864 fWidth = step; 864 fWidth = step; 865 fNofDivisions = ndiv; 865 fNofDivisions = ndiv; 866 CreateSolid(shape, Rpar, npar); 866 CreateSolid(shape, Rpar, npar); 867 867 868 return envVTE; 868 return envVTE; 869 } 869 } 870 870 871 G3VolTableEntry* G3Division::Dvt2() 871 G3VolTableEntry* G3Division::Dvt2() 872 { 872 { 873 // to be defined as const of this class 873 // to be defined as const of this class 874 G4double Rmin = 0.0001*cm; 874 G4double Rmin = 0.0001*cm; 875 875 876 // get parameters from mother 876 // get parameters from mother 877 G4String shape = fMVTE->GetShape(); 877 G4String shape = fMVTE->GetShape(); 878 G4double* Rpar = fMVTE->GetRpar(); 878 G4double* Rpar = fMVTE->GetRpar(); 879 G4int npar = fMVTE->GetNpar(); 879 G4int npar = fMVTE->GetNpar(); 880 880 881 // calculate the number of divisions 881 // calculate the number of divisions 882 G4int ndvmx = fNofDivisions; 882 G4int ndvmx = fNofDivisions; 883 G4double step = fStep; 883 G4double step = fStep; 884 G4double c0 = fC0; 884 G4double c0 = fC0; 885 885 886 if(fAxis == kPhi){ 886 if(fAxis == kPhi){ 887 step = step*deg; 887 step = step*deg; 888 c0 = c0*deg; 888 c0 = c0*deg; 889 } 889 } 890 else { 890 else { 891 step = step*cm; 891 step = step*cm; 892 c0 = c0*cm; 892 c0 = c0*cm; 893 } 893 } 894 894 895 G4int ndiv = G4int((fHighRange - c0 + Rmin)/ 895 G4int ndiv = G4int((fHighRange - c0 + Rmin)/step); 896 // to be added warning 896 // to be added warning 897 if (ndvmx > 255) ndvmx = 255; 897 if (ndvmx > 255) ndvmx = 255; 898 if (ndiv > ndvmx && ndvmx > 0 ) ndiv = ndvmx 898 if (ndiv > ndvmx && ndvmx > 0 ) ndiv = ndvmx; 899 899 900 // create envelope (if needed) 900 // create envelope (if needed) 901 G3VolTableEntry* envVTE = 0; 901 G3VolTableEntry* envVTE = 0; 902 G4double delta = std::abs((fHighRange - c0) 902 G4double delta = std::abs((fHighRange - c0) - ndiv*step); 903 if (std::abs(c0 - fLowRange) > Rmin) { 903 if (std::abs(c0 - fLowRange) > Rmin) { 904 envVTE 904 envVTE 905 = CreateEnvelope(shape, fHighRange-delta 905 = CreateEnvelope(shape, fHighRange-delta/2., c0+delta/2., Rpar, npar); 906 Rpar = envVTE->GetRpar(); 906 Rpar = envVTE->GetRpar(); 907 npar = envVTE->GetNpar(); 907 npar = envVTE->GetNpar(); 908 } 908 } 909 909 910 // set with for replica and create solid 910 // set with for replica and create solid 911 fWidth = step; 911 fWidth = step; 912 fNofDivisions = ndiv; 912 fNofDivisions = ndiv; 913 CreateSolid(shape, Rpar, npar); 913 CreateSolid(shape, Rpar, npar); 914 914 915 return envVTE; 915 return envVTE; 916 } 916 } 917 917