Problem 648 - Problem with G4BREPSolidPolyhedra
Summary: Problem with G4BREPSolidPolyhedra
Status: RESOLVED WORKSFORME
Alias: None
Product: Geant4
Classification: Unclassified
Component: geometry/solids (show other problems)
Version: 5.2
Hardware: HP Linux
: P2 normal
Assignee: Gabriele Cosmo
URL:
Depends on:
Blocks:
 
Reported: 2004-07-14 17:53 CEST by lisas
Modified: 2004-07-19 04:01 CEST (History)
0 users

See Also:


Attachments

Note You need to log in before you can comment on or make changes to this problem.
Description lisas 2004-07-14 17:53:06 CEST
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 :)
Comment 1 Gabriele Cosmo 2004-07-15 02:23:59 CEST
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) ?
Comment 2 lisas 2004-07-15 12:00:59 CEST
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
Comment 3 Gabriele Cosmo 2004-07-16 07:57:59 CEST
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 ?
Comment 4 lisas 2004-07-16 12:27:59 CEST
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);
Comment 5 Gabriele Cosmo 2004-07-19 04:01:59 CEST
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 .