I am running a modified version of the electromagnetic extended example TestEm2 whose geometry consists of a tapered cylinder containing liquid xenon with 2 0.5 mm partitions of polystyrene at 2 and 5 radiation lengths. The model runs with this geometry. However, if change even one of the 0.5 mm partitions to be a G4BREPSolidPolyhedra, the executable will hang when I try to run it. I'm not sure if it is a bug or simply a programming error on my part. G4VPhysicalVolume* Em2DetectorConstruction::ConstructVolumes() { G4double Radl = myMaterial->GetRadlen(); G4double dL = dLradl*Radl, dR = dRradl*Radl; EcalLength = dL; EcalSRadius = dR; EcalERadius = EndRadius(EcalLength,EcalSRadius); //EndRadius does a simple calculation G4double Thickness = 0.5*mm; G4double factor = 2*mm; // // Ecal // solidEcal = new G4Cons("Ecal",0.,2*EcalSRadius+Thickness, 0.,2*EcalERadius+Thickness,0.5*EcalLength,0.,360*deg); logicEcal = new G4LogicalVolume(solidEcal,myMaterial,"Ecal",0,0,0); physiEcal = new G4PVPlacement(0,G4ThreeVector(0.,0.,0.), "Ecal",logicEcal,0,false,0); // // Slice1 // //solidSlice1 = new G4Cons("Slice1",0.,EndRadius(2*Radl,EcalSRadius),0., // EndRadius((2*Radl+Thickness),EcalSRadius),0.5*Thickness,0.,360*deg); //works fine with solidSlice1 defined as above, however if defined like this: G4double S1RMAXVec[2]={1.1547*EndRadius(2*Radl,EcalSRadius), 1.1547*EndRadius((2*Radl+Thickness),EcalSRadius)}; G4double S1RMINVec[2]={0.,0.}; G4double S1Z_Values[2]={0.,Thickness}; solidSlice1 = new G4BREPSolidPolyhedra("Slice1",0.,2.*pi,6.,2.,0.,S1Z_Values,S1RMINVec,S1RMAXVec); logicSlice1 = new G4LogicalVolume(solidSlice1,Plastic,"Slice1",0,0,0); physiSlice1 = new G4PVPlacement(0,G4ThreeVector(0.,0.,(-7)*Radl+0.25*mm), "Slice1",logicSlice1,physiEcal,false,0); program hangs. It generates the geometry fine, and will start to generate an event, but then it never finishes. (I was curious and made RunAction cout dEstep, so a run does indeed start, I'm just not sure what happens after that): =============================================================== ===== ### Run 0 start. --------- Ranecu engine status --------- Initial seed (index) = 0 Current couple of seeds = 9876, 54321 ---------------------------------------- Start Run processing. ---> Begin of Event: 0 dEstep 1.77329 dEstep 0.632948 dEstep 0.0107371 dEstep 1.37221 dEstep 1.41794 dEstep 0.0259452 dEstep 0.689808 dEstep 3.40421 dEstep 0.213846 dEstep 0.0254319 dEstep 0.16502 dEstep 0.220984 dEstep 0.894355 dEstep 0.479399 dEstep 0.0446573 dEstep 0.0640842 dEstep 1.34006 dEstep 1.46371 dEstep 0.42053 dEstep 0.553527 dEstep 0.574093 dEstep 2.81689 dEstep 3.84388 dEstep 0.0874613 The exact same thing happens if I change any of the other volumes to be a G4BREPSolidPolyhedra (of tapered hexagonal geometry), i.e. it hangs at the same location. Any help would be greatly appreciated! Thanks :)
Can you please specify the exact values (with units) for: Radl, dLradl, EcalERadius ? Did you also verify that your final geometry setup does not contain overlapping volumes (making the physical dimensions of a solid depend on the material radiation-length is somehow... original and risky in this respect) ?
Yes, the same thing occurs when that slice is the only volume, besides Ecal (the mother volume). I'm using only one material, so the radiation length is actually a constant. Radl 2.80839 cm dRradl 1 dLradl 18 EcalERadius 4.58297 cm EcalSRadius 2.80839 cm EcalLength 50.551 cm
I still need some missing information in order to reproduce your exact geometry and setup: - Definition of the 'Plastic' material. - EndRadius(2*Radl,EcalSRadius) - EndRadius((2*Radl+Thickness),EcalSRadius) - Possibly also the particle with position/direction you're shooting. Is the rest (physics-list) the same as for TestEm2 ?
EndRadius(2*Radl,EcalSRadius) 3.00557 cm EndRadius((2*Radl+Thickness),EcalSRadius) 3.00732 cm //Plastic density = 1.032*g/cm3; G4Material* PStyrene = new G4Material(name="PolyStyrene", density, ncomponents=2); PStyrene->AddElement(C, natoms=8); PStyrene->AddElement(H, natoms=8); density = 1.032*g/cm3; G4Material* Sci = new G4Material(name="Scintillator", density, ncomponents=2); Sci->AddElement(C, natoms=9); Sci->AddElement(H, natoms=10); const G4int NENTRIES = 20; G4double TPB_PP[NENTRIES] = { 2.0*eV, 2.1*eV, 2.2*eV, 2.3*eV, 2.4*eV, 2.5*eV, 2.6*eV, 2.7*eV, 2.8*eV, 2.9*eV, 7.0*eV, 7.03*eV, 7.05*eV, 7.07*eV, 7.08*eV, 7.09*eV, 7.10*eV, 7.11*eV, 7.12*eV, 7.14*eV }; G4double TPB_RIND[NENTRIES] = { 1.59, 1.59, 1.59, 1.59, 1.59, 1.59, 1.59, 1.59, 1.59, 1.57, 1.57, 1.57, 1.57, 1.57, 1.57, 1.57, 1.57, 1.57, 1.54, 1.54 }; G4double TPB_ABSL[NENTRIES] = { 35.*cm, 35.*cm, 35.*cm, 35*cm, 35.*cm, 35.*cm, 35.*cm, 35.*cm, 35*cm, 35.*cm, 35.*cm, 35.*cm, 35.*cm, 35*cm, 35.*cm, 35.*cm, 35.*cm, 35.*cm, 35*cm, 35.*cm }; G4double TPB_EMISSPEC[NENTRIES] = { 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3 }; G4MaterialPropertiesTable *TPB_mt = new G4MaterialPropertiesTable(); //TPB_mt->AddProperty("RINDEX", TPB_PP, TPB_RIND, NENTRIES); TPB_mt->AddProperty("WLSABSLENGTH", TPB_PP, TPB_ABSL, NENTRIES); TPB_mt->AddProperty("WLSCOMPONENT", TPB_PP, TPB_EMISSPEC, NENTRIES); Sci->SetMaterialPropertiesTable(TPB_mt); Plastic = Sci; The physics list is the same except for the addition of optical photons and their processes: #include "G4Cerenkov.hh" #include "G4OpAbsorption.hh" #include "G4OpRayleigh.hh" #include "G4OpBoundaryProcess.hh" #include "G4Scintillation.hh" void LXePhysListEmStandard::ConstructOp() { // G4WLSProcess* theWLSProcess = new G4OpWLS("OpWLS"); G4Cerenkov* theCerenkovProcess = new G4Cerenkov("Cerenkov"); G4Scintillation* theScintillationProcess = new G4Scintillation("Scintillation"); G4OpAbsorption* theAbsorptionProcess = new G4OpAbsorption(); G4OpRayleigh* theRayleighScatteringProcess = new G4OpRayleigh(); G4OpBoundaryProcess* theBoundaryProcess = new G4OpBoundaryProcess(); // theCerenkovProcess->DumpPhysicsTable(); // theScintillationProcess->DumpPhysicsTable(); // theAbsorptionProcess->DumpPhysicsTable(); // theRayleighScatteringProcess->DumpPhysicsTable(); theCerenkovProcess->SetVerboseLevel(1); theScintillationProcess->SetVerboseLevel(1); theAbsorptionProcess->SetVerboseLevel(1); theRayleighScatteringProcess->SetVerboseLevel(1); theBoundaryProcess->SetVerboseLevel(1); G4int MaxNumPhotons = 200; theCerenkovProcess->SetTrackSecondariesFirst(true); theCerenkovProcess->SetMaxNumPhotonsPerStep(MaxNumPhotons); theScintillationProcess->SetScintillationYieldFactor(1.); theScintillationProcess->SetTrackSecondariesFirst(true); G4OpticalSurfaceModel themodel = unified; theBoundaryProcess->SetModel(themodel); theParticleIterator->reset(); while( (*theParticleIterator)() ){ G4ParticleDefinition* particle = theParticleIterator->value(); G4ProcessManager* pmanager = particle->GetProcessManager(); G4String particleName = particle->GetParticleName(); if (theCerenkovProcess->IsApplicable(*particle)) { pmanager->AddContinuousProcess(theCerenkovProcess); } if (theScintillationProcess->IsApplicable(*particle)) { pmanager->AddProcess(theScintillationProcess); pmanager->SetProcessOrderingToLast(theScintillationProcess, idxAtRest); pmanager->SetProcessOrderingToLast(theScintillationProcess, idxPostStep); } if (particleName == "opticalphoton") { G4cout << " AddDiscreteProcess to OpticalPhoton " << G4endl; pmanager->AddDiscreteProcess(theAbsorptionProcess); pmanager->AddDiscreteProcess(theRayleighScatteringProcess); pmanager->AddDiscreteProcess(theBoundaryProcess); } } // set ordering for AtRestDoIt } the particle/direction is just default: G4ParticleDefinition* particle = G4ParticleTable::GetParticleTable()->FindParticle("e-"); particleGun->SetParticleDefinition(particle); particleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.)); particleGun->SetParticleEnergy(500.*MeV);
OK, the situation is more clear now. The problem you're observing is eventually due only partially to the solid you're using. The loop you get eventually happens in the optical processes involved, and this problem has already been fixed in more recent releases. I would suggest you to upgrade to release 6.2 where several important fixes have been introduced in the optical processes (release 5.2 is rather outdated in this respect...). You can also find a patched version of G4BREPPolyhedra.cc here: http://pcitapiww.cern.ch/geant4/source/source/geant4/source/geometry/solids/BREPS/src/G4BREPSolidPolyhedra.cc which should appear soon in a patch to release 6.2. You can also alternatively use the G4Polyhedra CSG-like version instead. The new module for optical processes can also be found here: http://pcitapiww.cern.ch/geant4/source/source/geant4/source/processes/optical .