I had once made a half sphere volume to simulate the process of the gamma ray travel through the axis. But I found the ray penetrate the object as if it doesn't exit. I don't know why, and my pal also met the problem, and it still happen to us today. After that occured, I and my pal changed the object from sphere to box and found the phenomenon disapeared.
Can you please provide the details of you geometry setup ? Did you try using one of the latest releases of Geant4 ? Fixes have introduced to G4Sphere since release 6.0 ... Waiting for your reply.
HI, all Yes, I try the latest version, which is Geant4.6.0 then, and the problem seems still exit. And the attachments are my files. The first file is my scfPrimaryGeneratorAction.cc #include "scfPrimaryGeneratorAction.hh" #include "scfDetectorConstruction.hh" #include "G4Event.hh" #include "G4ParticleGun.hh" #include "G4ParticleTable.hh" #include "G4ParticleDefinition.hh" #include "Randomize.hh" #include "G4RunManager.hh" #include "G4ios.hh" //#include "fstream.h" #include "globals.hh" #include "math.h" extern G4double PrimaryBeamEDeposit; extern long IntegrateEnergySpectrum[250]; extern G4double AirCoef[37][2]; extern G4int runNm; scfPrimaryGeneratorAction::scfPrimaryGeneratorAction( scfDetectorConstruction* scfDC) :scfDetector(scfDC),rndmFlag("on") { G4int n_particle = 1; particleGun = new G4ParticleGun(n_particle); G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable(); G4String particleName; G4ParticleDefinition* particle = particleTable->FindParticle(particleName="gamma"); particleGun->SetParticleDefinition(particle); } scfPrimaryGeneratorAction::~scfPrimaryGeneratorAction() { delete particleGun; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void scfPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) {//G4cout<<"enter primaryGun"<<G4endl; //this function is called at the begining of event //rndmFlag=="on" G4int ID_energy=G4UniformRand()*IntegrateEnergySpectrum[249]; G4int Ibot=0; G4int Imid=124; G4int Itop=249; do{ if(ID_energy<IntegrateEnergySpectrum[Imid-1]) { Itop=Imid; Imid=(Itop+Ibot)/2; } else if(ID_energy>IntegrateEnergySpectrum[Imid-1]) { Ibot=Imid; Imid=(Ibot+Itop)/2; } else { Ibot=Imid-1; Itop=Imid; } } while(Ibot!=Imid&&Imid!=Itop); G4double ParticleEnergy=(Ibot+G4UniformRand())*10*keV; // //G4cout <<"ParticleEnergy="<<ParticleEnergy/keV<<"keV"<<G4endl; // //G4cout <<"ParticleEnergy="<<ParticleEnergy/MeV<<"MeV"<<G4endl; particleGun->SetParticleEnergy(ParticleEnergy); G4double theta=(runNm+1.)*deg/rad; G4double Momentx,Momenty,Momentz; Momentx=cos(theta); Momenty=0.; Momentz=-sin(theta); particleGun->SetParticleMomentumDirection(G4ThreeVector(Momentx,Momenty, Momentz)); particleGun->SetParticlePosition(G4ThreeVector(-Momentx*cm,Momenty*cm,- Momentz*cm)); particleGun->GeneratePrimaryVertex(anEvent); PrimaryBeamEDeposit+=ParticleEnergy*CalculateAirCoef(ParticleEnergy); } G4double scfPrimaryGeneratorAction::CalculateAirCoef(G4double GammaEnergy) { //G4cout <<"aircof[36][1]="<<AirCoef[36][1]<<"cm^2/g"<<G4endl;//for debug only G4double Factor[4]; G4double AirCoefficient=0.; int PosInTheTable=0; int EnergyRatio=1; for( int i=0;i<4;i++) { Factor[i]=1.; } EnergyRatio=GammaEnergy/AirCoef[0][0]; PosInTheTable=log(double(EnergyRatio))/log(10); PosInTheTable=PosInTheTable*8+1+EnergyRatio/pow(10,PosInTheTable); if(PosInTheTable<2) {AirCoefficient=AirCoef[0][1]; //cout<<"energy too lower"<<endl; } else if(PosInTheTable>=35) { for (int J=33;J<=36;J++) {for (int K=33;K<=36;K++) {if (J!=K) Factor[J-33]=Factor[J-33]*(GammaEnergy-AirCoef[K][0])/ (AirCoef[J][0]-AirCoef[K][0]); } AirCoefficient=AirCoefficient+Factor[J-33]*AirCoef[J][1]; } } else{ int ILOW=PosInTheTable-4; int IMED=PosInTheTable; int IUPPR=PosInTheTable+1; while(ILOW!=IMED&&IMED!=IUPPR) {if (GammaEnergy<AirCoef[IMED][0]) {IUPPR=IMED; IMED=(IUPPR+ILOW)/2; } else if(GammaEnergy>AirCoef[IMED][0]) {ILOW=IMED; IMED=(ILOW+IUPPR)/2; } else ILOW=IMED; } PosInTheTable=ILOW; if (fabs((GammaEnergy-AirCoef[PosInTheTable][0])/AirCoef [PosInTheTable][0])<0.01) AirCoefficient=AirCoef[PosInTheTable][1]; else if(fabs((GammaEnergy-AirCoef[PosInTheTable+1][0])/AirCoef [PosInTheTable+1][0])<0.01) AirCoefficient=AirCoef[PosInTheTable+1][1]; else {for (int J=PosInTheTable-1;J<=PosInTheTable+2;J++) {for (int K=PosInTheTable-1;K<=PosInTheTable+2;K++) if (J!=K) Factor[J-PosInTheTable+1]=Factor[J-PosInTheTable+1]* (GammaEnergy-AirCoef[K][0])/ (AirCoef[J][0]-AirCoef[K][0]); AirCoefficient=AirCoefficient+Factor[J-PosInTheTable+1] *AirCoef[J][1]; } } } <<GammaEnergy<<"MeV"<<"AirCoefficient"<<AirCoefficient<<endl;//for debug only! return(AirCoefficient); } The second one is my scfDetectorConstruction.cc #include "scfDetectorConstruction.hh" #include "G4Material.hh" #include "G4Sphere.hh" #include "G4LogicalVolume.hh" #include "G4PVPlacement.hh" #include "G4PVReplica.hh" #include "G4FieldManager.hh" #include "G4UniformMagField.hh" #include "G4TransportationManager.hh" #include "G4IdentityTrajectoryFilter.hh" #include "G4PropagatorInField.hh" #include "G4SDManager.hh" #include "G4RunManager.hh" #include "G4PhysicalVolumeStore.hh" #include "G4LogicalVolumeStore.hh" #include "G4SolidStore.hh" #include "G4VisAttributes.hh" #include "G4Colour.hh" #include "G4ios.hh" scfDetectorConstruction::scfDetectorConstruction():MaterialSteel(0),MaterialAir (0),MaterialVacuum(0) { // create commands for interactive definition of the calorimeter // detectorMessenger = new scfDetectorMessenger(this); DefineMaterials(); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... scfDetectorConstruction::~scfDetectorConstruction() { // delete detectorMessenger; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void scfDetectorConstruction::DefineMaterials() { //This function illustrates the possible ways to define materials G4String name, symbol; //a=mass of a mole; G4double a, z, density; //z=mean number of protons; //G4int iz; //iz=number of protons in an isotope; G4int ncomponents; G4double temperature, pressure; // // define Elements // a = 1.01*g/mole; G4Element* H = new G4Element(name="Hydrogen", symbol="H", z=1., a); a = 12.01*g/mole; G4Element* C = new G4Element(name="Carbon", symbol="C", z=6., a); a = 14.01*g/mole; G4Element* N = new G4Element(name="Nitrogen", symbol="N", z=7., a); a = 16.00*g/mole; G4Element* O = new G4Element(name="Oxygen", symbol="O", z=8., a); a = 55.85*g/mole; G4Element* Fe = new G4Element(name="Iron", symbol="Fe", z=26., a); // mixture by fractional mass density = 1.290*mg/cm3; G4Material* Air = new G4Material(name="Air", density, ncomponents=2); Air->AddElement(N,0.7); Air->AddElement(O,0.3); //Steel density = 7.83*g/cm3; G4Material* Steel = new G4Material(name="Steel",density, ncomponents=2); Steel->AddElement(Fe, 0.99); Steel->AddElement(C, 0.01); //Vacuum density = universe_mean_density; pressure = 3.e-18*pascal; temperature = 2.73*kelvin; G4Material* Vacuum = new G4Material(name="Galactic", z=1., a=1.01*g/mole, density,kStateGas,temperature,pressure); //default materials used in the plant MaterialSteel=Steel; MaterialAir=Air; MaterialVacuum=Vacuum; } G4VPhysicalVolume* scfDetectorConstruction::Construct() { G4PhysicalVolumeStore::GetInstance()->Clean(); G4LogicalVolumeStore::GetInstance()->Clean(); G4SolidStore::GetInstance()->Clean(); //XY is horizontal, Z is the height. G4Sphere* scfPlant = new G4Sphere("Plant", //its name 0.,1.8*m,0.,2*M_PI,0.,M_PI); //its size scfPlant_log = new G4LogicalVolume(scfPlant, //its solid MaterialVacuum, //its material "Plant"); //its name scfPlant_phy = new G4PVPlacement(0, //no rotation G4ThreeVector(), //at (0,0,0) scfPlant_log, //its logical volume "Plant", //its name 0, //its mother logic volume false, //no boolean operation 0); //copy number G4Sphere* scatVolume = new G4Sphere("scaVolume", //its name 0.,1.8*m,0.,2*M_PI,M_PI/2.,M_PI); //its size scatVolume_log = new G4LogicalVolume(scatVolume, //its solid MaterialSteel, //its material "scaVolume"); //its name scatVolume_phy = new G4PVPlacement(0, //no rotation G4ThreeVector(), //at (0,0,0) scatVolume_log, //its logical volume "scaVolume", //its name scfPlant_log, //its mother logic volume false, //no boolean operation 0); //copy number G4Sphere* testVolume = new G4Sphere("tstVolume", //its name 1.0*m,1.01*m,0.,2*M_PI,0.,M_PI/2.); //its size testVolume_log = new G4LogicalVolume(testVolume, //its solid MaterialAir, //its material "tstVolume"); //its name testVolume_phy = new G4PVPlacement(0, //no rotation G4ThreeVector(), //at (0,0,0) testVolume_log, //its logical volume "tstVolume", //its name scfPlant_log, //its mother logic volume false, //no boolean operation 0); //copy number // Visualization attributes scfPlant_log->SetVisAttributes (G4VisAttributes::Invisible); return scfPlant_phy; } You are highly appreciate for answering
One can check that for half sphere distance to in for p(0,0,0) and v(0,0,1) returns the correct value Rmin. Therefore the problem is not related to the sphere. Other reasons for strange behaviout could be: 1. zero gap betwenn world volume and steel, just increase a bit world radius; 2. Too small density of vacuum. Are you sure your vacuum is ~10^-18 Pascal? In any case a tar file with application (as simple as possible) will help to resolve the source of volume skip.
