I have used G4Tubs and G4Torus (the faulty code is the same) between 0 and Pi. It happened that I was sending some geantino (RayTracer) almost parallel to the plane X=0. When the solid is within a UnionSolid (See example), GEANT4 found correctly the intersection with the surface, but when it tried to compute the DistanceToOut, the part taking care of the cutting plane messed up, probably because of Round Off errors; the distanceToOut returned is 0. I propose the following change in G4Tubs.cc:1417 ->1555 by something like that if (!fPhiFullTube) { G4ThreeVector exit = p+sr*v; G4double exitPhi = std::atan2(exit.y(),exit.x()); if (exitPhi < 0) exitPhi + twopi; G4double srSE[2]; G4double phi[2] = {fSPhi, fSPhi + fDPhi}; if (exitPhi < fSPhi - halfAngTolerance && exitPhi> fSPhi + fDPhi + halfAngTolerance){ G4double phi; for (int i=0 ;i<2 ;++i){ if (std::cos(phi) == 0.) srSE[i] = (v.x() != 0.) ? - p.x()/v.x() : kInfinity; else { if (std::sin(phi) == 0.) srSE[i] = (v.y() != 0.) ? - p.x()/v.y() : kInfinity; else srSE[i] = (p.x()*std::tan(phi) - p.y())/(v.y() - v.x()*std::tan(phi)); } if (srSE[i] < 0.) srSE[i] = kInfinity; } if ( srSE[1] < sr) sr = srSE[1]; if ( srSE[2] < sr) sr = srSE[2]; } } Now I have to admit I did not completly understand the code, and I may have forgotten some particular case (I certainly hope not). The code in G4Torus is the same, and maybe also in other solids using the limitation (StartPhi, DeltaPhi) Regards
Can you please provide a test case where the faulty behavior can be reproduced ? Also can you please specify the version of Geant4 you're using and the platform (system, compiler) ?
Created attachment 52 [details] Modified ExampleN02 Since the problem seems to happen only when the volumes are in a combination of UnionSolid, I recreated the same volume I initially have problem with. The particle are geantino in the correct position/direction to trigger the problem.
Sorry, I thought I had attached the example already. Compiled with gcc 4.3.2 and gcc4.1.2 on Scientific Linux, tested on GEANT 4.9.0.p02 and 4.9.2.p02. The symptoms are different on OSX (Leopard) with g++4.0 : the particle is just stacked. To my understanding, the fact that the solid is in a UnionSolid makes the matter worse : The program loops infinitely in UnionSolid::DistanceToOut. Regards.
Dear Pierre, Thank you for sending the example and for your comments. By looking in the example, I found that you use a Boolean operations for solids with sharing surfaces(phi planes). Most probably, this causing the problem with tracking. Boolean operations in Geant4 have to be used with solids with overlaps. If the solids have sharing surfaces, the tracking can be confused, but also the visualisation can be wrong. In your case, you can add small overlaps( 1 micron), in order to avoid sharing surfaces. Please, let us know, if this helps. Best Regards, Tatiana Nikitina
Thanks for the quick answer. Unfortunately, having overlapping volume does not seem to solve the problem. I keep looking around to see if I can find some reason to that.
Dear Pierre, Thank you for the message and for the test. After your message, I tried to simplify provided example in order to check the replay of the Solids. If I use only one Rotated LargeTorus and I shoot the particle along the phi axe, the problem is appears(bouncing between Target and World). I think, there is a problem in the Replay of G4Torus for phi section, I will continue to investigate. I have also tested the G4Tubs in the same configurations. The Replay of G4Tubs is correct. Best Regards, Tatiana Nikitina
Dear Tatiana, I am sorry, but I am not sure about what the "Replay" is. I have to correct my previous answer: it seems that the overlapping helps but does not solve the problem. In any case I thought that the tolerance are making to solids put beside each other slightly overlapping... Is it wrong ? Anyway, in G4Torus::DistanceToIn, something is bothering me: Lines 965 and 997 there are if ( sphi < 0 ) { sphi = 0 ; } I do not understand that test. It would be more logical to me if instead of having if (sphi < snxt) we had : if (0 < sphi && nsphi < snxt) which would reflect the fact that we are close to the limit and we have a possibility to cross the phi plane. For sure that modification makes my problem disappear... Cheers, Pierre
Hello again (the last one today I hope) In the case SPhi<0, the RayTracer shows only the part of the Torus which has Phi>0. I think it comes from SolveNumericJT at the line 287 if (theta < 0) { theta += twopi; } In the case SPhi < 0, the following condition cannot be true if ( (theta - fSPhi >= - kAngTolerance*0.5) && (theta - (fSPhi + fDPhi) <= kAngTolerance*0.5) ) Since the constructor is taking care of the case where fSPhi < -Pi, then I guess it would be enough to have if (fSPhi >=0 && theta < 0) { theta += twopi; } Cheers, Pierre
Once again me (I forgot to commit that message before) in G4Torus::DistanceToOut(p,v) if I understand correctly, there can be an intersection between the track and the phi-planes only if the "pDistX" (X=E or S) and "compX" have the same sign. If it is the case. I believe that the lines 1356 and 1424 are incorrect and sphi should receive kInfinity instead of "0". I also wonder about the case managed line 1283 if ( (compS>=0) && (compE>=0) ) { sphi=kInfinity; } else { sphi=0; } I have a feeling that it is wrong but I cannot be sure. Cheers, Pierre
Dear Pierre, "Replay" of the solid is what replaying DistanceToIn(), DistanceToOut() or SurfaceNormal() for given points, directions. I case of the G4Torus it seems to be a problem with DistanceToIn() or DistanceToOut() when the point is On phi surface. It can be a problem with use of kAngTolerance, specialy for Torus of (0,pi), beacause 'angle' becomming a plane. Thank you for the comments about the code, I will have a look. Yes, overlapping is what you did. Best Regards, Tatiana
*** Problem 1086 has been marked as a duplicate of this problem. ***
Dear Pierre, Thank very much you for usefull comments and suggestions. You are right(comment#8) there is a problem with theta. Your proposal is correct, only the kAngTolerance has to be added due to numerical accupacy. Also case with Sphi<-pi has to be added with theta=-twopi; This fix will be available in the next december release. The other problem is the wrong 'Replay' (DistanceToIn=DistanceToOut=0) for the geantinos on the Phi Surface and going almost parallel to PhiSection. These errors are due to the incorrect use of kAngTolerance. It has to be carrefully added and tested. Best Regards, Tatiana Nikitina
Hi Tatiana, I am sorry to come back with it, but several versions of GEANT later, I still have bugs in G4Torus and G4Tubs when I have rays grazing the Phi-planes. I think I localized the problem, at least for G4Torus. In G4Torus::DistanceToOut( const G4ThreeVector& p, const G4ThreeVector& v, const G4bool calcNorm, G4bool *validNorm, G4ThreeVector *n ) const when we are testing the Phi-planes, we are comparing pDist (shortest distance between the point and the plane) (E/S) to the tolerance (lines 1232, 1259, 1237 and 1395 of G4Torus.cc). This makes little sense, considering that we just calculated sphi which is the distance to the phi-plane in the direction of the momentum. On the other hand, I agree that if sphi < Tolerance/2, then it can be approximate to zero. DistanceToIn seems to follow that policy, even if sphi is calculated only if Dist>Tolerance/2. Which means that if we are very close to the phi-plane, the phi-plane is not taken into account for DistToIn... Is that logic? I have the feeling that we have a similar problem with G4Tubs, although since the code is rather different I am not yet sure. By the way, is it wise to have so many different codes for the solids of revolution? The phi-plane calculations for instance should be very similar for the yub, sphere, torus... Well, I guess you ruled that out for a good reason. regards, Pierre
Dear Pierre, Thank you for your comments. The problem in Phi section was not solved yet. I agree with you that it will be good make code in Tubs and Torus similar, where it is possible. I was working on it. When it is ready, I will let you know. Best Regards, Tatiana Nikitina
Dear Pierre, Thank you again for your remarques. You was right to pointing out the code in DistanceToOut(p,v,...) with ... if(pDistS>-kCarTolerance*0.5){sphi = 0;} ... This is causing a problem. This condition (pDistS>-kCarTolerance*0.5) means only that the point p is situated almost on Phi plane. But the distance to out is depending on direction v, it is NOT always zero. I'm testing a fix for this problem. It will be available in the next release. If you are interested, I can send you the fix already. Best Regards, Tatiana Nikitina
The fix is now available in the patched release 9.3.p02.
