[1316] | 1 | // |
---|
| 2 | // ******************************************************************** |
---|
| 3 | // * License and Disclaimer * |
---|
| 4 | // * * |
---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
| 7 | // * conditions of the Geant4 Software License, included in the file * |
---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
| 9 | // * include a list of copyright holders. * |
---|
| 10 | // * * |
---|
| 11 | // * Neither the authors of this software system, nor their employing * |
---|
| 12 | // * institutes,nor the agencies providing financial support for this * |
---|
| 13 | // * work make any representation or warranty, express or implied, * |
---|
| 14 | // * regarding this software system or assume any liability for its * |
---|
| 15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
| 16 | // * for the full disclaimer and the limitation of liability. * |
---|
| 17 | // * * |
---|
| 18 | // * This code implementation is the result of the scientific and * |
---|
| 19 | // * technical work of the GEANT4 collaboration. * |
---|
| 20 | // * By using, copying, modifying or distributing the software (or * |
---|
| 21 | // * any work based on the software) you agree to acknowledge its * |
---|
| 22 | // * use in resulting scientific publications, and indicate your * |
---|
| 23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
| 24 | // ******************************************************************** |
---|
| 25 | // |
---|
| 26 | // |
---|
| 27 | // $Id: testG4PVDivision.cc,v 1.5 2009/05/14 14:19:32 ivana Exp $ |
---|
| 28 | // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ |
---|
| 29 | // |
---|
| 30 | // test for G4PVDivision classes |
---|
| 31 | // |
---|
| 32 | // 26.05.03 - P.Arce Initial version |
---|
| 33 | // ******************************************************************** |
---|
| 34 | |
---|
| 35 | #include <assert.h> |
---|
| 36 | #include <fstream> |
---|
| 37 | #include <vector> |
---|
| 38 | |
---|
| 39 | #include "G4ios.hh" |
---|
| 40 | #include "globals.hh" |
---|
| 41 | |
---|
| 42 | #include "G4Navigator.hh" |
---|
| 43 | |
---|
| 44 | #include "G4LogicalVolume.hh" |
---|
| 45 | #include "G4VPhysicalVolume.hh" |
---|
| 46 | #include "G4PVPlacement.hh" |
---|
| 47 | #include "G4PVParameterised.hh" |
---|
| 48 | #include "G4VPVParameterisation.hh" |
---|
| 49 | #include "G4Box.hh" |
---|
| 50 | #include "G4Trd.hh" |
---|
| 51 | #include "G4Trap.hh" |
---|
| 52 | #include "G4Tubs.hh" |
---|
| 53 | #include "G4Cons.hh" |
---|
| 54 | #include "G4Para.hh" |
---|
| 55 | #include "G4Polycone.hh" |
---|
| 56 | #include "G4Polyhedra.hh" |
---|
| 57 | |
---|
| 58 | #include "G4GeometryManager.hh" |
---|
| 59 | |
---|
| 60 | #include "G4RotationMatrix.hh" |
---|
| 61 | #include "G4ThreeVector.hh" |
---|
| 62 | |
---|
| 63 | #include "Randomize.hh" |
---|
| 64 | #include "G4PVReplica.hh" |
---|
| 65 | |
---|
| 66 | #include "G4PVDivision.hh" |
---|
| 67 | |
---|
| 68 | G4String theSolidTypeStr; |
---|
| 69 | G4String thePVTypeStr; |
---|
| 70 | G4double theWorldDim = 1*m; |
---|
| 71 | G4int numberOfPoints = 1000; |
---|
| 72 | G4int theNReplicas; |
---|
| 73 | G4double theWidth; |
---|
| 74 | G4double theOffset; |
---|
| 75 | G4VSolid* theParentSolid; |
---|
| 76 | std::vector<G4LogicalVolume*> theParentLogs; |
---|
| 77 | std::vector<G4VPhysicalVolume*> theParentPhyss; |
---|
| 78 | std::vector<G4VSolid*> theChildSolids; |
---|
| 79 | std::vector<G4LogicalVolume*> theChildLogs; |
---|
| 80 | std::vector<EAxis> theAxis; |
---|
| 81 | std::vector<G4double> theWidths; |
---|
| 82 | std::vector<G4double> theExtraPars; |
---|
| 83 | G4int theDivType; |
---|
| 84 | |
---|
| 85 | enum SolidType{g4box, g4trd, g4tube, g4tubs, g4cone, g4cons, g4polycone, g4polyhedra }; |
---|
| 86 | enum PVType{pvDivision, pvReplica, pvPlacement }; |
---|
| 87 | |
---|
| 88 | //-------------------------------------------------------------------------- |
---|
| 89 | void initialize(); |
---|
| 90 | void calculateParentSolid( SolidType solType ); |
---|
| 91 | void calculateChildSolids( SolidType solType, G4VSolid* pSolid ); |
---|
| 92 | void calculateAxis( SolidType solType ); |
---|
| 93 | void buildOutputName( SolidType& soltype, PVType pvtype ); |
---|
| 94 | void starttest( const std::vector<G4String>& vsarg ); |
---|
| 95 | G4int checkNumberOfArguments( const G4String& st, G4int narg ); |
---|
| 96 | PVType getPVType( const G4String& pvt ); |
---|
| 97 | SolidType getSolidType( const G4String& st ); |
---|
| 98 | void generateRandomPoints(); |
---|
| 99 | void generateScanPointsForBox(); |
---|
| 100 | void generateScanPointsForTube(); |
---|
| 101 | void generateScanPointsForTrd(); |
---|
| 102 | void generateScanPointsForPolycone(); |
---|
| 103 | G4VPhysicalVolume* BuildGeometry( SolidType solType, PVType pvType ); |
---|
| 104 | G4bool testG4Navigator1(G4VPhysicalVolume *pTopNode); |
---|
| 105 | G4bool testG4Navigator2(G4VPhysicalVolume *pTopNode); |
---|
| 106 | // World geometry is a box 1 X 1 X 3. |
---|
| 107 | // User select the type of solid |
---|
| 108 | // Inside it three solids of the chosen type are placed along Z |
---|
| 109 | // Each of these solids is divided along a different axis |
---|
| 110 | // |
---|
| 111 | // Then a set of points is generated and it is tested in which division copy they are and which is the DistanceToOut in the three directions (X,Y,Z) |
---|
| 112 | |
---|
| 113 | //-------------------------------------------------------------------------- |
---|
| 114 | int main( G4int argc, char** argv ) |
---|
| 115 | { |
---|
| 116 | // first argument is type of divisioning (repli/divis), second is type of solid |
---|
| 117 | std::vector<G4String> vsarg; |
---|
| 118 | for( G4int jj = 0; jj < argc; jj ++ ) { |
---|
| 119 | vsarg.push_back( G4String(argv[jj] ) ); |
---|
| 120 | } |
---|
| 121 | |
---|
| 122 | starttest( vsarg ); |
---|
| 123 | } |
---|
| 124 | |
---|
| 125 | //-------------------------------------------------------------------------- |
---|
| 126 | void starttest( const std::vector<G4String>& vsarg ) |
---|
| 127 | { |
---|
| 128 | G4int narg = vsarg.size(); |
---|
| 129 | if( narg == 1 ) { |
---|
| 130 | thePVTypeStr = "divis"; |
---|
| 131 | theSolidTypeStr = "box"; |
---|
| 132 | } else if( narg == 2 ){ |
---|
| 133 | // wrong number |
---|
| 134 | checkNumberOfArguments( " ", narg ); |
---|
| 135 | } else { |
---|
| 136 | G4int divTypeSet = checkNumberOfArguments( vsarg[2], narg ); |
---|
| 137 | if( divTypeSet ) { |
---|
| 138 | theDivType = atoi( vsarg[divTypeSet] ); |
---|
| 139 | } else { |
---|
| 140 | theDivType = 0; |
---|
| 141 | } |
---|
| 142 | thePVTypeStr = G4String(vsarg[1]); |
---|
| 143 | theSolidTypeStr = G4String(vsarg[2]); |
---|
| 144 | } |
---|
| 145 | if( narg > 3 ){ |
---|
| 146 | for( G4int ii = 3; ii < narg; ii++ ){ |
---|
| 147 | theExtraPars.push_back( atof( vsarg[ii] ) ); |
---|
| 148 | } |
---|
| 149 | } |
---|
| 150 | PVType pvtype = getPVType( thePVTypeStr ); |
---|
| 151 | SolidType soltype = getSolidType( theSolidTypeStr ); |
---|
| 152 | |
---|
| 153 | // PVType pvtype = pvDivision; |
---|
| 154 | initialize(); |
---|
| 155 | if( soltype == g4tube || soltype == g4tubs ) { |
---|
| 156 | // generateRandomPoints(); |
---|
| 157 | generateScanPointsForTube(); |
---|
| 158 | }else if( soltype == g4cone || soltype == g4cons ) { |
---|
| 159 | // generateRandomPoints(); |
---|
| 160 | generateScanPointsForTube(); |
---|
| 161 | } else if( soltype == g4box ) { |
---|
| 162 | // generateRandomPoints(); |
---|
| 163 | generateScanPointsForBox(); |
---|
| 164 | } else if( soltype == g4trd ) { |
---|
| 165 | // generateRandomPoints(); |
---|
| 166 | generateScanPointsForTrd(); |
---|
| 167 | } else if( soltype == g4polycone ) { |
---|
| 168 | // generateRandomPoints(); |
---|
| 169 | generateScanPointsForPolycone(); |
---|
| 170 | } else { |
---|
| 171 | generateRandomPoints(); |
---|
| 172 | } |
---|
| 173 | |
---|
| 174 | G4VPhysicalVolume *myTopNode; |
---|
| 175 | myTopNode=BuildGeometry( soltype, pvtype ); // Build the geometry |
---|
| 176 | G4GeometryManager::GetInstance()->CloseGeometry(false); |
---|
| 177 | testG4Navigator1(myTopNode); |
---|
| 178 | testG4Navigator2(myTopNode); |
---|
| 179 | // Repeat tests but with full voxels |
---|
| 180 | G4GeometryManager::GetInstance()->OpenGeometry(); |
---|
| 181 | G4GeometryManager::GetInstance()->CloseGeometry(true); |
---|
| 182 | testG4Navigator1(myTopNode); |
---|
| 183 | testG4Navigator2(myTopNode); |
---|
| 184 | |
---|
| 185 | G4cout << " theParentSolid " << G4endl; |
---|
| 186 | G4cout << *theParentSolid << G4endl; |
---|
| 187 | for( size_t ii = 0; ii < theChildSolids.size(); ii++) { |
---|
| 188 | G4cout << " theChildSolid after tracking " << " "<< ii << G4endl; |
---|
| 189 | G4cout << *theChildSolids[ii] << G4endl; |
---|
| 190 | } |
---|
| 191 | |
---|
| 192 | G4GeometryManager::GetInstance()->OpenGeometry(); |
---|
| 193 | } |
---|
| 194 | |
---|
| 195 | //-------------------------------------------------------------------------- |
---|
| 196 | void generateRandomPoints() |
---|
| 197 | { |
---|
| 198 | std::ofstream fout("points.lis"); |
---|
| 199 | G4double posX, posY, posZ; |
---|
| 200 | |
---|
| 201 | for( G4int ii = 0; ii < numberOfPoints; ii++ ) { |
---|
| 202 | posX = CLHEP::RandFlat::shoot(-theWorldDim, theWorldDim); |
---|
| 203 | posY = CLHEP::RandFlat::shoot(-theWorldDim, theWorldDim); |
---|
| 204 | posZ = CLHEP::RandFlat::shoot(-3*theWorldDim, 3*theWorldDim); |
---|
| 205 | fout << posX << " " << posY << " " << posZ << G4endl; |
---|
| 206 | } |
---|
| 207 | } |
---|
| 208 | |
---|
| 209 | //-------------------------------------------------------------------------- |
---|
| 210 | void generateScanPointsForBox() |
---|
| 211 | { |
---|
| 212 | std::ofstream fout("points.lis"); |
---|
| 213 | G4int ii; |
---|
| 214 | |
---|
| 215 | G4int nPointsPerReplica = 2; |
---|
| 216 | G4int numberOfPoints = theNReplicas*nPointsPerReplica; |
---|
| 217 | // For division along X |
---|
| 218 | G4ThreeVector centre(0.,0.,-2*theWorldDim); |
---|
| 219 | for( ii = 0; ii < numberOfPoints; ii++ ) { |
---|
| 220 | // any Z, any Y |
---|
| 221 | G4ThreeVector pR( 0., theWorldDim/100., theWorldDim/100. ); |
---|
| 222 | G4double X = -theWorldDim + (ii+0.001) * 2*theWorldDim/numberOfPoints; |
---|
| 223 | pR.setX( X ); |
---|
| 224 | pR += centre; |
---|
| 225 | fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl; |
---|
| 226 | } |
---|
| 227 | |
---|
| 228 | // For division along Y |
---|
| 229 | centre = G4ThreeVector(0.,0.,0.); |
---|
| 230 | for( ii = 0; ii < numberOfPoints; ii++ ) { |
---|
| 231 | // any X, any Z |
---|
| 232 | G4ThreeVector pR( theWorldDim/100., 0., theWorldDim/100. ); |
---|
| 233 | G4double Y = -theWorldDim + (ii+0.001) * 2*theWorldDim/numberOfPoints; |
---|
| 234 | pR.setY( Y ); |
---|
| 235 | pR += centre; |
---|
| 236 | fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl; |
---|
| 237 | } |
---|
| 238 | |
---|
| 239 | // For division along Z |
---|
| 240 | centre = G4ThreeVector(0.,0.,2*theWorldDim); |
---|
| 241 | for( ii = 0; ii < numberOfPoints; ii++ ) { |
---|
| 242 | // any X, any Y |
---|
| 243 | G4ThreeVector pR( theWorldDim/100., 0., theWorldDim/100. ); |
---|
| 244 | G4double Z = -theWorldDim + (ii+0.001) * 2*theWorldDim/numberOfPoints; |
---|
| 245 | pR.setZ( Z ); |
---|
| 246 | pR += centre; |
---|
| 247 | fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl; |
---|
| 248 | } |
---|
| 249 | } |
---|
| 250 | |
---|
| 251 | //-------------------------------------------------------------------------- |
---|
| 252 | void generateScanPointsForTrd() |
---|
| 253 | { |
---|
| 254 | std::ofstream fout("points.lis"); |
---|
| 255 | G4int ii; |
---|
| 256 | |
---|
| 257 | G4int nPointsPerReplica = 2; |
---|
| 258 | G4int numberOfPoints = theNReplicas*nPointsPerReplica; |
---|
| 259 | // For division along X |
---|
| 260 | G4ThreeVector centre(0.,0.,-2*theWorldDim); |
---|
| 261 | for( ii = 0; ii < numberOfPoints; ii++ ) { |
---|
| 262 | // any Z, any Y |
---|
| 263 | G4ThreeVector pR( 0., theWorldDim/100., theWorldDim/100. ); |
---|
| 264 | G4double X = -theWorldDim*0.75 + (ii+0.001) * 1.5*theWorldDim/numberOfPoints; |
---|
| 265 | pR.setX( X ); |
---|
| 266 | pR += centre; |
---|
| 267 | fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl; |
---|
| 268 | } |
---|
| 269 | |
---|
| 270 | // For division along Y |
---|
| 271 | centre = G4ThreeVector(0.,0.,0.); |
---|
| 272 | for( ii = 0; ii < numberOfPoints; ii++ ) { |
---|
| 273 | // any X, any Z |
---|
| 274 | G4ThreeVector pR( theWorldDim/100., 0., theWorldDim/100. ); |
---|
| 275 | G4double Y = -theWorldDim*0.75 + (ii+0.001) * 1.5*theWorldDim/numberOfPoints; |
---|
| 276 | pR.setY( Y ); |
---|
| 277 | pR += centre; |
---|
| 278 | fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl; |
---|
| 279 | } |
---|
| 280 | |
---|
| 281 | // For division along Z |
---|
| 282 | centre = G4ThreeVector(0.,0.,2*theWorldDim); |
---|
| 283 | for( ii = 0; ii < numberOfPoints; ii++ ) { |
---|
| 284 | // any X, any Y |
---|
| 285 | G4ThreeVector pR( theWorldDim/100., 0., theWorldDim/100. ); |
---|
| 286 | G4double Z = -theWorldDim + (ii+0.001) * 2*theWorldDim/numberOfPoints; |
---|
| 287 | pR.setZ( Z ); |
---|
| 288 | pR += centre; |
---|
| 289 | fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl; |
---|
| 290 | } |
---|
| 291 | } |
---|
| 292 | |
---|
| 293 | //-------------------------------------------------------------------------- |
---|
| 294 | void generateScanPointsForTube() |
---|
| 295 | { |
---|
| 296 | std::ofstream fout("points.lis"); |
---|
| 297 | G4int ii; |
---|
| 298 | |
---|
| 299 | G4int nPointsPerReplica = 2; |
---|
| 300 | G4int numberOfPoints = theNReplicas*nPointsPerReplica; |
---|
| 301 | G4cout << " numberOfPoints " << numberOfPoints << G4endl; |
---|
| 302 | // For division along R |
---|
| 303 | G4ThreeVector centre(0.,0.,-2*theWorldDim); |
---|
| 304 | for( ii = 0; ii < numberOfPoints; ii++ ) { |
---|
| 305 | // any Z, phi close initial phi |
---|
| 306 | G4double phi = 1.*deg; |
---|
| 307 | //t? if( theExtraPars.size() > 0 ) phi = (theExtraPars[0]+1.)*deg; |
---|
| 308 | G4ThreeVector pR( std::cos(phi), std::sin(phi), theWorldDim/100. ); |
---|
| 309 | G4double rho = (ii+0.001) * theWorldDim/numberOfPoints; |
---|
| 310 | pR.setRho( rho ); |
---|
| 311 | pR += centre; |
---|
| 312 | fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl; |
---|
| 313 | } |
---|
| 314 | |
---|
| 315 | // For division along phi |
---|
| 316 | centre = G4ThreeVector(0.,0.,0.); |
---|
| 317 | for( ii = 0; ii < numberOfPoints; ii++ ) { |
---|
| 318 | G4double phi = (ii+0.001) * 360*deg/numberOfPoints; |
---|
| 319 | // any rho (=1), any Z |
---|
| 320 | G4ThreeVector pR( std::cos(phi), std::sin(phi), theWorldDim/100. ); |
---|
| 321 | pR += centre; |
---|
| 322 | fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl; |
---|
| 323 | } |
---|
| 324 | |
---|
| 325 | // For division along Z |
---|
| 326 | centre = G4ThreeVector(0.,0.,2*theWorldDim); |
---|
| 327 | for( ii = 0; ii < numberOfPoints; ii++ ) { |
---|
| 328 | //any rho (=1), phi close initial phi |
---|
| 329 | G4double phi = 1.*deg; |
---|
| 330 | //t? if( theExtraPars.size() > 0 ) phi = (theExtraPars[0]+1.)*deg; |
---|
| 331 | G4ThreeVector pR( std::cos(phi), std::sin(phi),0.); |
---|
| 332 | G4double Z = -theWorldDim + (ii+0.001) * 2*theWorldDim / numberOfPoints; |
---|
| 333 | pR.setZ( Z ); |
---|
| 334 | pR += centre; |
---|
| 335 | fout << pR.x() << " " << pR.y() << " " << pR.z() << G4endl; |
---|
| 336 | } |
---|
| 337 | } |
---|
| 338 | |
---|
| 339 | //-------------------------------------------------------------------------- |
---|
| 340 | void generateScanPointsForPolycone() |
---|
| 341 | { |
---|
| 342 | generateScanPointsForTube(); |
---|
| 343 | } |
---|
| 344 | |
---|
| 345 | //-------------------------------------------------------------------------- |
---|
| 346 | // Build simple geometry: |
---|
| 347 | // world is |
---|
| 348 | G4VPhysicalVolume* BuildGeometry( SolidType solType, PVType pvType ) |
---|
| 349 | { |
---|
| 350 | G4int ii; |
---|
| 351 | // The world volume |
---|
| 352 | // |
---|
| 353 | G4Box *worldSolid= new G4Box ("Big Cube", theWorldDim, theWorldDim, 3*theWorldDim); |
---|
| 354 | |
---|
| 355 | G4LogicalVolume *worldLog=new G4LogicalVolume(worldSolid,0, |
---|
| 356 | "World",0,0,0); |
---|
| 357 | // Logical with no material,field, |
---|
| 358 | // sensitive detector or user limits |
---|
| 359 | |
---|
| 360 | G4PVPlacement *worldPhys=new G4PVPlacement(0,G4ThreeVector(0,0,0), |
---|
| 361 | "World",worldLog, |
---|
| 362 | 0,false,0); |
---|
| 363 | // Note: no parent pointer set |
---|
| 364 | |
---|
| 365 | // build three boxes |
---|
| 366 | // A set of boxes |
---|
| 367 | calculateParentSolid( solType ); |
---|
| 368 | |
---|
| 369 | //parent logical volumes do not depend on SolidType, PVType |
---|
| 370 | theParentLogs.push_back( new G4LogicalVolume(theParentSolid,0, "Parent1",0,0,0) ); |
---|
| 371 | theParentLogs.push_back( new G4LogicalVolume(theParentSolid,0, "Parent2",0,0,0) ); |
---|
| 372 | theParentLogs.push_back( new G4LogicalVolume(theParentSolid,0, "Parent3",0,0,0) ); |
---|
| 373 | |
---|
| 374 | //parent physical volumes positions do not depend on SolidType, PVType |
---|
| 375 | for( ii = 0; ii < 3; ii++ ) { |
---|
| 376 | theParentPhyss.push_back( new G4PVPlacement( 0, G4ThreeVector(0.,0.,(ii-1)*2*theWorldDim), theParentLogs[ii] , "parent", worldLog, 0, 0 ) ); |
---|
| 377 | } |
---|
| 378 | |
---|
| 379 | // children |
---|
| 380 | calculateChildSolids( solType, theParentSolid ); |
---|
| 381 | calculateAxis( solType ); |
---|
| 382 | for( ii = 0; ii < 3; ii++ ) { |
---|
| 383 | theChildLogs.push_back( new G4LogicalVolume(theChildSolids[ii], 0, "child",0,0,0) ); |
---|
| 384 | } |
---|
| 385 | |
---|
| 386 | if( pvType == pvDivision ) { |
---|
| 387 | for( ii = 0; ii < 3; ii++ ) { |
---|
| 388 | if( theDivType == 0 ) { |
---|
| 389 | new G4PVDivision("child", theChildLogs[ii], theParentLogs[ii], |
---|
| 390 | theAxis[ii], |
---|
| 391 | theNReplicas, |
---|
| 392 | theWidths[ii], |
---|
| 393 | theOffset ); |
---|
| 394 | G4cout << "division NDIVandWIDTH axis " |
---|
| 395 | << theNReplicas << " " << theWidths[ii]<< " " << theOffset << " " << theAxis[ii] << G4endl; |
---|
| 396 | } else if( theDivType == 1 ) { |
---|
| 397 | new G4PVDivision("child", theChildLogs[ii], theParentLogs[ii], |
---|
| 398 | theAxis[ii], |
---|
| 399 | theWidths[ii], |
---|
| 400 | theOffset ); |
---|
| 401 | G4cout << "division WIDTH " << theWidths[ii]<< " " << theOffset << G4endl; |
---|
| 402 | }else if( theDivType == 2 ) { |
---|
| 403 | new G4PVDivision("child", theChildLogs[ii], theParentLogs[ii], |
---|
| 404 | theAxis[ii], |
---|
| 405 | theNReplicas, |
---|
| 406 | theOffset ); |
---|
| 407 | G4cout << "division NDIV " << theNReplicas << " " << theOffset << G4endl; |
---|
| 408 | } |
---|
| 409 | } |
---|
| 410 | } else if( pvType == pvReplica ) { |
---|
| 411 | |
---|
| 412 | for( ii = 0; ii < 3; ii++ ) { |
---|
| 413 | new G4PVReplica("child", theChildLogs[ii], theParentLogs[ii], |
---|
| 414 | theAxis[ii], |
---|
| 415 | theNReplicas, |
---|
| 416 | theWidths[ii], |
---|
| 417 | theOffset ); |
---|
| 418 | G4cout << "replica " << theNReplicas << " " << theWidths[ii]<< " " << theOffset << G4endl; |
---|
| 419 | } |
---|
| 420 | } |
---|
| 421 | return worldPhys; |
---|
| 422 | } |
---|
| 423 | |
---|
| 424 | //-------------------------------------------------------------------------- |
---|
| 425 | // |
---|
| 426 | // Test LocateGlobalPointAndSetup |
---|
| 427 | // |
---|
| 428 | G4bool testG4Navigator1(G4VPhysicalVolume *pTopNode) |
---|
| 429 | { |
---|
| 430 | G4Navigator myNav; |
---|
| 431 | G4VPhysicalVolume *located; |
---|
| 432 | myNav.SetWorldVolume(pTopNode); |
---|
| 433 | |
---|
| 434 | G4String foutname = "points." + thePVTypeStr + "." + theSolidTypeStr + ".out"; |
---|
| 435 | std::ofstream fout(foutname); |
---|
| 436 | std::ifstream fin("points.lis"); |
---|
| 437 | G4double posX, posY, posZ; |
---|
| 438 | for( G4int ii = 0; ii < numberOfPoints; ii++ ) { |
---|
| 439 | fin >> posX >> posY >> posZ; |
---|
| 440 | if( fin.eof() ) break; |
---|
| 441 | located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(posX,posY,posZ),0,false); |
---|
| 442 | G4ThreeVector pos(posX,posY,posZ); |
---|
| 443 | fout << ii << pos << " " << located->GetName() << " " << located->GetCopyNo(); |
---|
| 444 | if( theSolidTypeStr == "tubs" || theSolidTypeStr == "tube" ){ |
---|
| 445 | // fout << " " << pos.phi()/deg << G4endl; |
---|
| 446 | fout << G4endl; |
---|
| 447 | } else { |
---|
| 448 | fout << G4endl; |
---|
| 449 | } |
---|
| 450 | } |
---|
| 451 | |
---|
| 452 | return true; |
---|
| 453 | } |
---|
| 454 | |
---|
| 455 | //-------------------------------------------------------------------------- |
---|
| 456 | // |
---|
| 457 | // Test Stepping |
---|
| 458 | // |
---|
| 459 | G4bool testG4Navigator2(G4VPhysicalVolume *pTopNode) |
---|
| 460 | { |
---|
| 461 | G4Navigator myNav; |
---|
| 462 | G4VPhysicalVolume *located; |
---|
| 463 | G4double Step,physStep,safety; |
---|
| 464 | G4ThreeVector* Hat = new G4ThreeVector[3]; |
---|
| 465 | Hat[0] = G4ThreeVector(1,0,0); |
---|
| 466 | Hat[1] = G4ThreeVector(0,1,0); |
---|
| 467 | Hat[2] = G4ThreeVector(0,0,1); |
---|
| 468 | |
---|
| 469 | myNav.SetWorldVolume(pTopNode); |
---|
| 470 | |
---|
| 471 | G4String foutname = "step." + thePVTypeStr + "." + theSolidTypeStr + ".out"; |
---|
| 472 | std::ofstream fout(foutname); |
---|
| 473 | std::ifstream fin("points.lis"); |
---|
| 474 | G4double posX, posY, posZ; |
---|
| 475 | for( G4int ii = 0; ii < numberOfPoints; ii++ ) { |
---|
| 476 | if( fin.eof() ) break; |
---|
| 477 | fin >> posX >> posY >> posZ; |
---|
| 478 | located=myNav.LocateGlobalPointAndSetup(G4ThreeVector(posX,posY,posZ)); |
---|
| 479 | physStep=kInfinity; |
---|
| 480 | for( G4int jj = 0; jj < 3; jj++ ) { |
---|
| 481 | Step=myNav.ComputeStep(G4ThreeVector(posX,posY,posZ),Hat[jj],physStep,safety); |
---|
| 482 | fout << Step << " "; |
---|
| 483 | } |
---|
| 484 | fout << G4endl; |
---|
| 485 | } |
---|
| 486 | |
---|
| 487 | return true; |
---|
| 488 | } |
---|
| 489 | |
---|
| 490 | //-------------------------------------------------------------------------- |
---|
| 491 | void initialize() |
---|
| 492 | { |
---|
| 493 | theNReplicas = 10; |
---|
| 494 | theOffset = 0.; |
---|
| 495 | } |
---|
| 496 | |
---|
| 497 | //-------------------------------------------------------------------------- |
---|
| 498 | void calculateParentSolid( SolidType solType ) |
---|
| 499 | { |
---|
| 500 | if( solType == g4box ) { |
---|
| 501 | theParentSolid = new G4Box("parent", theWorldDim, theWorldDim, theWorldDim); |
---|
| 502 | }else if( solType == g4trd ) { |
---|
| 503 | theParentSolid = new G4Trd("parent", theWorldDim/2., theWorldDim, theWorldDim/2., theWorldDim, theWorldDim); |
---|
| 504 | }else if( solType == g4tube ) { |
---|
| 505 | theParentSolid = new G4Tubs("parent", 0., theWorldDim, theWorldDim, 0., 360.*deg); |
---|
| 506 | }else if( solType == g4tubs ) { |
---|
| 507 | theParentSolid = new G4Tubs("parent", 0., theWorldDim, theWorldDim, theExtraPars[0]*deg, theExtraPars[1]*deg); |
---|
| 508 | }else if( solType == g4cone ) { |
---|
| 509 | theParentSolid = new G4Cons("parent", 1.E-6*mm, theWorldDim, 1.E-6*mm, theWorldDim/2., theWorldDim, 0., 360.*deg); |
---|
| 510 | }else if( solType == g4cons ) { |
---|
| 511 | theParentSolid = new G4Cons("parent", 1.E-6*mm, theWorldDim/2., 1.E-6*mm, theWorldDim, theWorldDim, theExtraPars[0]*deg, theExtraPars[1]*deg); |
---|
| 512 | }else if( solType == g4polycone ) { |
---|
| 513 | |
---|
| 514 | G4double* zPlane = new G4double[4]; |
---|
| 515 | zPlane[0] = -theWorldDim; zPlane[1] = -theWorldDim/2.; zPlane[2] = theWorldDim/2; zPlane[3] = theWorldDim; |
---|
| 516 | G4double* rInner = new G4double[4]; |
---|
| 517 | rInner[0] = theWorldDim/10.; rInner[1] = 0.; rInner[2] = 0.; rInner[3] = theWorldDim/10.; |
---|
| 518 | G4double* rOuter = new G4double[4]; |
---|
| 519 | rOuter[0] = theWorldDim; rOuter[1] = theWorldDim*0.9; rOuter[2] = theWorldDim*0.9; rOuter[3] = theWorldDim; |
---|
| 520 | |
---|
| 521 | theParentSolid = new G4Polycone("parent", 0., 360.*deg, 4, zPlane, rInner, rOuter); |
---|
| 522 | }else if( solType == g4polyhedra ) { |
---|
| 523 | |
---|
| 524 | G4double* zPlane = new G4double[4]; |
---|
| 525 | zPlane[0] = -theWorldDim; zPlane[1] = -theWorldDim/2.; zPlane[2] = theWorldDim/2; zPlane[3] = theWorldDim; |
---|
| 526 | G4double* rInner = new G4double[4]; |
---|
| 527 | rInner[0] = theWorldDim/10.; rInner[1] = 0.; rInner[2] = 0.; rInner[3] = theWorldDim/10.; |
---|
| 528 | G4double* rOuter = new G4double[4]; |
---|
| 529 | rOuter[0] = theWorldDim; rOuter[1] = theWorldDim*0.9; rOuter[2] = theWorldDim*0.9; rOuter[3] = theWorldDim; |
---|
| 530 | |
---|
| 531 | theParentSolid = new G4Polyhedra("parent", 0., 360.*deg, theNReplicas, 4, zPlane, rInner, rOuter); |
---|
| 532 | } |
---|
| 533 | } |
---|
| 534 | |
---|
| 535 | //-------------------------------------------------------------------------- |
---|
| 536 | void calculateChildSolids( SolidType solType, G4VSolid* pSolid ) |
---|
| 537 | { |
---|
| 538 | if( solType == g4box ) { |
---|
| 539 | theWidths.push_back( 2*theWorldDim / theNReplicas ); |
---|
| 540 | theWidths.push_back( 2*theWorldDim / theNReplicas ); |
---|
| 541 | theWidths.push_back( 2*theWorldDim / theNReplicas ); |
---|
| 542 | |
---|
| 543 | theChildSolids.push_back( new G4Box("child", theWorldDim/theNReplicas, theWorldDim, theWorldDim) ); |
---|
| 544 | theChildSolids.push_back( new G4Box("child", theWorldDim, theWorldDim/theNReplicas, theWorldDim) ); |
---|
| 545 | theChildSolids.push_back( new G4Box("child", theWorldDim, theWorldDim, theWorldDim/theNReplicas) ); |
---|
| 546 | |
---|
| 547 | }else if( solType == g4trd ) { |
---|
| 548 | theWidths.push_back( 1.5*theWorldDim / theNReplicas ); |
---|
| 549 | theWidths.push_back( 1.5*theWorldDim / theNReplicas ); |
---|
| 550 | theWidths.push_back( 2*theWorldDim / theNReplicas ); |
---|
| 551 | |
---|
| 552 | theChildSolids.push_back( new G4Trap("child", theWorldDim/2./theNReplicas, theWorldDim/theNReplicas, theWorldDim, theWorldDim, theWorldDim) ); // solid dimensions will be recalculated |
---|
| 553 | theChildSolids.push_back( new G4Trap("child", theWorldDim, theWorldDim, theWorldDim/2./theNReplicas, theWorldDim/theNReplicas, theWorldDim) ); // solid dimensions will be recalculated |
---|
| 554 | theChildSolids.push_back( new G4Trd("child", theWorldDim, theWorldDim, theWorldDim, theWorldDim, theWorldDim/theNReplicas) ); |
---|
| 555 | |
---|
| 556 | }else if( solType == g4tube || solType == g4tubs ) { |
---|
| 557 | G4Tubs* pTubs = reinterpret_cast<G4Tubs*>( pSolid ); |
---|
| 558 | |
---|
| 559 | theWidths.push_back( (pTubs->GetInnerRadius()+pTubs->GetOuterRadius()) / theNReplicas ); |
---|
| 560 | theWidths.push_back( pTubs->GetDeltaPhiAngle() / theNReplicas ); |
---|
| 561 | theWidths.push_back( 2*pTubs->GetZHalfLength() / theNReplicas ); |
---|
| 562 | |
---|
| 563 | theChildSolids.push_back( new G4Tubs("child", pTubs->GetInnerRadius(), pTubs->GetInnerRadius()+theWidths[0], pTubs->GetZHalfLength(), pTubs->GetStartPhiAngle(), pTubs->GetDeltaPhiAngle() )); |
---|
| 564 | theChildSolids.push_back( new G4Tubs("child", pTubs->GetInnerRadius(), pTubs->GetOuterRadius(), pTubs->GetZHalfLength(), pTubs->GetStartPhiAngle(), theWidths[1] )); |
---|
| 565 | theChildSolids.push_back( new G4Tubs("child", pTubs->GetInnerRadius(), pTubs->GetOuterRadius(), theWidths[2]/2., pTubs->GetStartPhiAngle(), pTubs->GetDeltaPhiAngle() )); |
---|
| 566 | }else if( solType == g4cone || solType == g4cons ) { |
---|
| 567 | G4Cons* pCons = reinterpret_cast<G4Cons*>( pSolid ); |
---|
| 568 | |
---|
| 569 | theWidths.push_back( (pCons->GetInnerRadiusMinusZ()+pCons->GetOuterRadiusMinusZ()) / theNReplicas ); |
---|
| 570 | theWidths.push_back( pCons->GetDeltaPhiAngle() / theNReplicas ); |
---|
| 571 | theWidths.push_back( 2*pCons->GetZHalfLength() / theNReplicas ); |
---|
| 572 | |
---|
| 573 | G4double fwidthPlus = 0.; |
---|
| 574 | if( pCons->GetInnerRadiusMinusZ() != 0. ) fwidthPlus = theWidths[0] * pCons->GetInnerRadiusPlusZ() / pCons->GetInnerRadiusMinusZ(); |
---|
| 575 | theChildSolids.push_back( new G4Cons("child", pCons->GetInnerRadiusMinusZ(), pCons->GetInnerRadiusMinusZ()+theWidths[0], pCons->GetInnerRadiusPlusZ(), pCons->GetInnerRadiusPlusZ()+fwidthPlus, pCons->GetZHalfLength(), pCons->GetStartPhiAngle(), pCons->GetDeltaPhiAngle() )); |
---|
| 576 | theChildSolids.push_back( new G4Cons("child", pCons->GetInnerRadiusMinusZ(), pCons->GetOuterRadiusMinusZ(), pCons->GetInnerRadiusPlusZ(), pCons->GetOuterRadiusPlusZ(), pCons->GetZHalfLength(), pCons->GetStartPhiAngle(), theWidths[1] )); |
---|
| 577 | theChildSolids.push_back( new G4Cons("child", pCons->GetInnerRadiusMinusZ(), pCons->GetOuterRadiusMinusZ(), pCons->GetInnerRadiusPlusZ(), pCons->GetOuterRadiusPlusZ(), theWidths[2]/2., pCons->GetStartPhiAngle(), pCons->GetDeltaPhiAngle() )); |
---|
| 578 | |
---|
| 579 | }else if( solType == g4polycone ) { |
---|
| 580 | G4Polycone* pPCone = reinterpret_cast<G4Polycone*>( pSolid ); |
---|
| 581 | G4PolyconeHistorical* origparam = pPCone->GetOriginalParameters() ; |
---|
| 582 | G4int nz = origparam->Num_z_planes; |
---|
| 583 | theWidths.push_back( (origparam->Rmax[0] - origparam->Rmin[0] ) / theNReplicas ); |
---|
| 584 | theWidths.push_back( (pPCone->GetEndPhi() - pPCone->GetStartPhi() ) / theNReplicas ); |
---|
| 585 | theWidths.push_back( (origparam->Z_values[1] - origparam->Z_values[0] ) / theNReplicas ); |
---|
| 586 | |
---|
| 587 | G4double* zPlane = new G4double[4]; |
---|
| 588 | zPlane[0] = -theWorldDim; zPlane[1] = -theWorldDim/2.; zPlane[2] = theWorldDim/2; zPlane[3] = theWorldDim; |
---|
| 589 | G4double* rInner = new G4double[4]; |
---|
| 590 | rInner[0] = theWorldDim/10.; rInner[1] = 0.; rInner[2] = 0.; rInner[3] = theWorldDim/10.; |
---|
| 591 | G4double* rOuter = new G4double[4]; |
---|
| 592 | rOuter[0] = theWorldDim; rOuter[1] = theWorldDim*0.9; rOuter[2] = theWorldDim*0.9; rOuter[3] = theWorldDim; |
---|
| 593 | |
---|
| 594 | theChildSolids.push_back( new G4Polycone("child", 0., 360.*deg, nz, zPlane, rInner, rOuter) ); |
---|
| 595 | theChildSolids.push_back( new G4Polycone("child", 0., 360.*deg, nz, zPlane, rInner, rOuter) ); |
---|
| 596 | theChildSolids.push_back( new G4Polycone("child", 0., 360.*deg, nz, zPlane, rInner, rOuter) ); |
---|
| 597 | |
---|
| 598 | }else if( solType == g4polyhedra ) { |
---|
| 599 | G4Polyhedra* pPhedra = reinterpret_cast<G4Polyhedra*>( pSolid ); |
---|
| 600 | G4PolyhedraHistorical* origparam = pPhedra->GetOriginalParameters() ; |
---|
| 601 | G4int ns = origparam->numSide; |
---|
| 602 | G4int nz = origparam->Num_z_planes; |
---|
| 603 | theWidths.push_back( (origparam->Rmax[0] - origparam->Rmin[0] ) / theNReplicas ); |
---|
| 604 | theWidths.push_back( (pPhedra->GetEndPhi() - pPhedra->GetStartPhi() ) / theNReplicas ); |
---|
| 605 | theWidths.push_back( (origparam->Z_values[1] - origparam->Z_values[0] ) / theNReplicas ); |
---|
| 606 | |
---|
| 607 | G4double* zPlane = new G4double[4]; |
---|
| 608 | zPlane[0] = -theWorldDim; zPlane[1] = -theWorldDim/2.; zPlane[2] = theWorldDim/2; zPlane[3] = theWorldDim; |
---|
| 609 | G4double* rInner = new G4double[4]; |
---|
| 610 | rInner[0] = theWorldDim/10.; rInner[1] = 0.; rInner[2] = 0.; rInner[3] = theWorldDim/10.; |
---|
| 611 | G4double* rOuter = new G4double[4]; |
---|
| 612 | rOuter[0] = theWorldDim; rOuter[1] = theWorldDim*0.9; rOuter[2] = theWorldDim*0.9; rOuter[3] = theWorldDim; |
---|
| 613 | |
---|
| 614 | theChildSolids.push_back( new G4Polyhedra("child", 0., 360.*deg, ns, nz, zPlane, rInner, rOuter) ); |
---|
| 615 | theChildSolids.push_back( new G4Polyhedra("child", 0., 360.*deg, ns, nz, zPlane, rInner, rOuter) ); |
---|
| 616 | theChildSolids.push_back( new G4Polyhedra("child", 0., 360.*deg, ns, nz, zPlane, rInner, rOuter) ); |
---|
| 617 | |
---|
| 618 | } |
---|
| 619 | for( size_t ii = 0; ii < theChildSolids.size(); ii++) { |
---|
| 620 | G4cout << theChildSolids[0] << "child solid built 0 " << *theChildSolids[0] << G4endl; |
---|
| 621 | G4cout << " theChildSolid built " << ii << G4endl; |
---|
| 622 | G4cout << *theChildSolids[ii] << G4endl; |
---|
| 623 | } |
---|
| 624 | } |
---|
| 625 | |
---|
| 626 | //-------------------------------------------------------------------------- |
---|
| 627 | void calculateAxis( SolidType solType ) |
---|
| 628 | { |
---|
| 629 | if( solType == g4box ) { |
---|
| 630 | theAxis.push_back( kXAxis ); |
---|
| 631 | theAxis.push_back( kYAxis ); |
---|
| 632 | theAxis.push_back( kZAxis ); |
---|
| 633 | }else if( solType == g4trd ) { |
---|
| 634 | theAxis.push_back( kXAxis ); |
---|
| 635 | theAxis.push_back( kYAxis ); |
---|
| 636 | theAxis.push_back( kZAxis ); |
---|
| 637 | }else if( solType == g4tube || solType == g4tubs ) { |
---|
| 638 | theAxis.push_back( kRho ); |
---|
| 639 | theAxis.push_back( kPhi ); |
---|
| 640 | theAxis.push_back( kZAxis ); |
---|
| 641 | }else if( solType == g4cone || solType == g4cons || solType == g4polycone || solType == g4polyhedra ) { |
---|
| 642 | theAxis.push_back( kRho ); |
---|
| 643 | theAxis.push_back( kPhi ); |
---|
| 644 | theAxis.push_back( kZAxis ); |
---|
| 645 | } |
---|
| 646 | } |
---|
| 647 | |
---|
| 648 | //-------------------------------------------------------------------------- |
---|
| 649 | SolidType getSolidType( const G4String& st ) |
---|
| 650 | { |
---|
| 651 | SolidType stype = g4box; |
---|
| 652 | |
---|
| 653 | if( st == "box" ) { |
---|
| 654 | stype = g4box; |
---|
| 655 | }else if( st == "trd" ){ |
---|
| 656 | stype = g4trd; |
---|
| 657 | }else if( st == "tube" ){ |
---|
| 658 | stype = g4tube; |
---|
| 659 | }else if( st == "tubs" ){ |
---|
| 660 | stype = g4tubs; |
---|
| 661 | }else if( st == "cone" ){ |
---|
| 662 | stype = g4cone; |
---|
| 663 | }else if( st == "cons" ){ |
---|
| 664 | stype = g4cons; |
---|
| 665 | }else if( st == "polycone" ){ |
---|
| 666 | stype = g4polycone; |
---|
| 667 | }else if( st == "polyhedra" ){ |
---|
| 668 | stype = g4polyhedra; |
---|
| 669 | } else { |
---|
| 670 | G4Exception(" Input solid type can only be 'box', 'trd', 'tube', 'tubs', 'cons', 'polycone', 'polyhedra' " ); |
---|
| 671 | } |
---|
| 672 | return stype; |
---|
| 673 | } |
---|
| 674 | |
---|
| 675 | //-------------------------------------------------------------------------- |
---|
| 676 | G4int checkNumberOfArguments( const G4String& st, G4int narg ) |
---|
| 677 | { |
---|
| 678 | G4bool nok = false; |
---|
| 679 | if( narg == 1 ) return 0; |
---|
| 680 | G4int nArgExcess; |
---|
| 681 | if( st == "box" || st == "trd" || st == "tube" || st == "cone" ){ |
---|
| 682 | nArgExcess = narg - 3; |
---|
| 683 | }else if( st == "tubs" || st == "cons" ){ |
---|
| 684 | nArgExcess = narg - 5; |
---|
| 685 | }else if( st == "polycone" ){ |
---|
| 686 | nArgExcess = narg - 3; |
---|
| 687 | }else if( st == "polyhedra" ){ |
---|
| 688 | nArgExcess = narg - 3; |
---|
| 689 | } else { |
---|
| 690 | nArgExcess = 10; |
---|
| 691 | } |
---|
| 692 | |
---|
| 693 | G4int divTypeSet = 0; |
---|
| 694 | if( nArgExcess == 0 ) { |
---|
| 695 | nok = true; |
---|
| 696 | } else if( nArgExcess == 1 ) { |
---|
| 697 | divTypeSet = narg - 1; //last argument sets the division type |
---|
| 698 | nok = true; |
---|
| 699 | } |
---|
| 700 | |
---|
| 701 | if( !nok ){ |
---|
| 702 | G4cerr << " Number Of arguments " << narg << " for solid type " << st << G4endl; |
---|
| 703 | G4Exception(" Number of arguments incorrect for input solid type, please check method checkNumberOfArguments " ); |
---|
| 704 | } |
---|
| 705 | |
---|
| 706 | return divTypeSet; |
---|
| 707 | } |
---|
| 708 | |
---|
| 709 | //-------------------------------------------------------------------------- |
---|
| 710 | PVType getPVType( const G4String& pvt ) |
---|
| 711 | { |
---|
| 712 | G4cout << pvt << G4endl; |
---|
| 713 | PVType vtype = pvPlacement; |
---|
| 714 | |
---|
| 715 | if( pvt == "divis"){ |
---|
| 716 | vtype = pvDivision; |
---|
| 717 | }else if( pvt == "repli") { |
---|
| 718 | vtype = pvReplica; |
---|
| 719 | }else if( pvt == "place") { |
---|
| 720 | vtype = pvPlacement; |
---|
| 721 | } else { |
---|
| 722 | G4Exception(" Input PV type can only be 'divis', 'repli' or 'place' " ); |
---|
| 723 | } |
---|
| 724 | return vtype; |
---|
| 725 | } |
---|