| Summary: | infitite loop in G4SubtractionSolid::DistanceToIn | ||
|---|---|---|---|
| Product: | Geant4 | Reporter: | andrea.di.simone |
| Component: | geometry/solids | Assignee: | tatiana.nikitina |
| Status: | RESOLVED WORKSFORME | ||
| Severity: | major | CC: | andrea.di.simone |
| Priority: | P3 | ||
| Version: | 9.3 | ||
| Hardware: | PC | ||
| OS: | Linux | ||
| Attachments: |
the TRD and the two tubs, as UNION
the TRD and the two tubs, as SUBTRACTION the gdml file for the problematic volume First fix to the problem |
||
Created attachment 125 [details]
the TRD and the two tubs, as UNION
This shows the TRD and the two tubs we subtract. To give you an idea the two tubs are added, rather than subtracted
Created attachment 126 [details]
the TRD and the two tubs, as SUBTRACTION
The final solid, after both the drills.
The orange point is the problematic one. It is on the surface of the first drill, but definitely outside the 2nd drill.
I attached a couple of images, to better clarify our geometry. Comments are attache. One shows the TRD and the two tubs, in an UNION rather than a subtraction, just to visualize the solids. I also added the real subtraction. The orange point is the problematic one. It is on the surcafe of the first drill, but definitely away from the second one. Hi Andrea, indeed such protection exists already in the latest release 9.4.p02 (I would actually suggest you to try out your setup on it, to verify if you still experience the problem!). The case you report deserves anyhow verification, as, assuming the construct is correctly done, it may indicate a mis-behavior of the solid. Should be also verified that the solids being subtracted have NOT shared surfaces (or vertices) subject to the subtraction, as this may lead to problems of this kind (although, from the pictures you provide this doesn't seem the case...). Gabriele Dear Andrea, Could you provide geometry file to test? If you have any additional information, like position, direction of particle when it enters in infinite loop, it will help to understand the problem. Tatiana Hi Tatiana, Gabriele, thanks for your feedback. I'll try 9.4.p02 and let you know, but it will take some time since I'll have to find a problematic event on that release (I'm assuming the simulation results will change from 9.3 to 9.4). I have just seen the new protection in 9.4.p02, and it looks great, thanks. My only concern with it is that if the problem turns out to be solid-dependent (as you seem to suggest), the policy you chose (i.e. to abort the whole event) would dramatically bias the results of the simulation: events with a high track density on that solid will most likely be aborted. Since this code is used for accelerator background estimation, the effect of the bias is actually visible. The same happens to us, btw, with the "StuckTrack" exception: we also have several stuck tracks, and we had to implement a workaround (a custom exception handler) to avoid the bias due to the aborted events. Of course this is unrelated to this bug, but please just keep in mind that navigation problems need a better solution than just aborting the events. I am attaching a gdml file reproducing the problematic volume. The point where the problem occurs (the orange dot in the images attached) is (24.981741,-27.40732,104.24314)*mm in the local ref system of the solid. I don't have the direction at hand, but I'll find it and post it here in the next hours. Thanks for your help, Andrea. Created attachment 127 [details]
the gdml file for the problematic volume
These are p,v of the problem: (24.981741,-27.40732,104.24314) (0.73468058,-0.43257898,0.52260871) in mm, in the local ref frame of the subtraction. Hi, I tried the latest release, and I can still find the infinite loop (though at a different event). The bright side is that I will be able to test any patch you may provide. However, I think the fix is not working the way you wanted it to. You are not breaking the infinite loop, hence ou are just throwing an infinite number of exceptions. Information about the event being killed is not propagated down to the solid and the loop keeps on without stopping. So to summarize: I still have the infinite loop (in the same solid, btw). The exception helps in the sense that it makes the problem more evident, but the loop is not stopped (this deserves to be treated like a separate bug, probably: if you want I can open a separate bug report). Cheers, Andrea. Since I assume you are going to fix again that loop, I have a couple more suggestions: - please add to the exception text the (p,v) values - please notice that if solidB is a G4DisplacedSolid, its name is always "placedB", which makes it useless for debugging. I would recommend something similar to (dynamic_cast<G4DisplacedSolid*> (fPtrSolidB))->GetConstituentMovedSolid()->GetName() Thanks for your help, Andrea. Created attachment 128 [details]
First fix to the problem
This a first fix. DistanceToIn(p,v) will return 0, then G4Navigator will take care of track, pushing if needed.
Dear Andrea, Thank you for comments and for trails. Here is a first fix, based on your comments. In case of infinite loop DistanceToIn(p,v) will replay 0 and the navigation will take care of the track. I hope it can help as first fix, but I continue further investigations. I have included provaded geometry in Geant4 tests. I don't have infinite loop in the given point and direction, but by moving it around, I hope to detect this infinite loop. Best Regards, Tatiana Hi Tatiana, thanks for tyhe fix. I'll try it now. Concerning the solid problem, maybe you could give me some hints on where to look at. I have all setup to reproduce the problem, and I am stuck anyway until this is solved, so I'm happy to help debugging. Could you confirm that that call to fPtrSolidB->DistanceToOut(p,v) when p is already kOutside of solidB is really what you mean? In this case I'll have a look at G4Tubs::DistanceToOut and see why it returns zero, since that is clearly wrong. Thanks for your help, Andrea. Hi Andrea, Thank you for helping. I'm looking in the algorithm of G4SubtractionSolid::DistanceToIn(p,v). It seems to be strange asking for DistanceToOut(p,v) from point OutsideB. And it seems to be normal that Solid(G4Tubs) returns DistanceToOut=0 for point outside. I continue investigate. On the first look the geometry description is correct. There aren't subtractions of solids with shared surfaces. Could you check if it really the case? Tatiana Hi Tatiana, yes the consistency of the solids is the first thing we checked, and there are no shared surfaces (not that we could spot). Concerning G4SubtractionSolid::DistanceToIn, there are two things I am not sure of (I stressed them in the comments to the code snippet): - the distToOut call for a point that is kOutside: as you said it sounds weird. if the result in such situation must be zero by design (as you seem to imply), then the loop is very likely to become infinite (though this would not be happening always) - the kOutside returned for a point on the surface of the 1st drill. Is this the expected answer for such "boudary" conditions? If you are sure that both of these are correct, then I am afraid we'll relly need to dig into G4Tubs. Cheers, Andrea.
Hi andrea,
According the code in G4SubtractionSolid::Inside(p)
...
if(( positionA == kInside && positionB == kSurface) ||
( positionB == kOutside && positionA == kSurface) ||
( positionA == kSurface && positionB == kSurface &&
( fPtrSolidA->SurfaceNormal(p) -
fPtrSolidB->SurfaceNormal(p) ).mag2() >
1000*G4GeometryTolerance::GetInstance()->GetRadialTolerance() ) )
{
return kSurface;
}
else
{
return kOutside;
}
...
If point is on Surface of SolidA, bur outside SolidB :
positionA=kSurface, positionB=kOutside => solid(A-B) has to return kSurface.
If it is returning kOutside, may be there is a problem with G4Tubs::Inside(),
as you said.
Could you add a printsout to G4Tubs::Inside(p)
to verify what answer it is giving and coordinates of p?
Thank you for your help.
Tatiana
Hi Andrea, Thank you again for all help. Finally I could reproduce the infinite loop in my tests. It seems to be a problem in G4Tubs::Inside(p) indeed. I'm preparing the fix. Tatiana Hi Tatiana, thanks, that's indeed great news. I'll be waiting for your fix. Cheers, Andrea. Can I ask what the plans are to release this fix? Of course I could apply it privately to my code, but I woul rather prefer it to be included in an official release. Any chance it will make it into 9.5? Or maybe a third patch of 9.4? Cheers, Andrea. Hi Andrea, I'm working on the fix for G4Tubs,it will be ready in next few days. Then it will be included together with the fix for G4SubtractionSolid in following patches and releases. Cheers, Tatiana Hi Tatiana, thanks the info. In the meantime, would you mind explaining what the bug is and what its consequences are? Apart from the stuck jobs we are observing all sorts of weird behavior in our simulation, and of course G4Tubs::Inside is something we use quite a lot... Thanks for your help, Andrea. Hi Andrea, I am still investigating a potential problem in G4Tubs. But in the same time, I tried to change the way how the geometry setup is created. When you subtract the second Tube, it has crossing surfaces with the previous Tube( long Z) and this kind of Subtractions can be a potential problem in tracking. Instead of Subtraction of the first tube from Trd and then of the second, Gabriele proposed first make the union of two tubes and then make the subtraction of this union from Trd. In this case crossing surfaces became a common surface, it is safer way to do Boolean operations with solids. Indeed, I tested it in your geometry setup, it works, no infinite loops are detected. Cheers, Tatiana Dear Andrea, By further investigations I found that the problem is comming from given direction v, it is near to unit vector, but is not a unit vector. May be it come from truncation when the values were printed. But when I normalise this vector v, no more erros are reported. Best Regards, Tatiana Nikitina Hi Tatiana, thanks for your tests. I just don't understand what I should be doing from my side. Are you going to add an extra normalization on the G4 side? Cheers, Andrea. I also notice that this has been classified ad "resolved". Can you please point me to the relevant patch, so that I can test it? Andrea. Dear Andrea, As it was mentioned before, the proposal is to change the way how the boolean operations are done. In your case first G4Tubs is subtracted from G4Trd and then G4Tubs is subtracted from the first subtraction. In this case you subtract solids with shared surfaces, it can lead to problems in the navigation. The new proposal is following : first make a Union of two Tubs and then subtract this Union from G4Trd. In this way you avoid subtraction of shared surfaces. I tested it, it works without entering in infinite loop. For the given direction, could be it possible to check if you enter a unit vector? If it is the case, it has to be verified where precision is lost. Cheers, Tatiana |
Hi, I was investigating some cases of stuck jobs, and I traced them back to G4SubtractionSolid::DistanceToIn. Our situation (with some level of simplification) is the following: we have a TRD from which we subtract two tubs (actually, it is the same solid, subtracted at different offsets). so the steps are: 1) 1st drill: subtract tube from TRD 2) 2nd drill: subtract tube from 1st drill The problem occurs when in the "final" solid, i.e. the 2nd drill, I ask what DistToIn is from a point on the surface of the 1st drill. I paste here the code from latest release in LXR, and some annotations with my understanding of the problem: G4double G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p, const G4ThreeVector& v ) const { G4double dist = 0.0,disTmp = 0.0 ; // we are processing the 2nd drill. hence solidA is 1st drill and solidB is a tubs // the point is kOutside of solidB, so this if is skipped and we enter the else a few lines from here if ( fPtrSolidB->Inside(p) != kOutside ) { dist = fPtrSolidB->DistanceToOut(p,v) ; // ,calcNorm,validNorm,n) ; if( fPtrSolidA->Inside(p+dist*v) != kInside ) { do { disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ; if(disTmp == kInfinity) { return kInfinity ; } dist += disTmp ; if( Inside(p+dist*v) == kOutside ) { disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ; dist += disTmp ; } } while( Inside(p+dist*v) == kOutside ) ; } } else // p outside A, start in A { // ok we enter here, and we descend one level. we ask distToIn wrt to the 1st drill. Since the point is exactly on the surface of the first tube we subtracted, this returns 0 dist = fPtrSolidA->DistanceToIn(p,v) ; if( dist == kInfinity ) // past A, hence past A\B { return kInfinity ; } else { // dist=0, p+dist*v=p, and Inside(p) returns kOutside, hence the loop below starts // One possible problem here: if the point is on the surface of the drill, is kOutside the expected answer? while( Inside(p+dist*v) == kOutside ) // pushing loop { disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ; dist += disTmp ; // Note: this is asking solidB->DistanceToOut(p) (which returns 0), but if we are here then we are already kOutside of solidB. From visual inspection of the problematic point, it actually is outside of solidB, but not at zero distance from it. if( Inside(p+dist*v) == kOutside ) { // Inside(p) is still kOutside, so we enter here as well disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ; // solidA is the TRD with one hole, and the point is on the surface of the hole. // fPtrSolidA->DistanceToIn(p,v) returns 0 if(disTmp == kInfinity) // past A, hence past A\B { return kInfinity ; } dist += distTmp ; // since distTmp is zero, dist stays zero. this leads to an infinite loop } } } } return dist ; } Do you have any idea of what may be going on? We checked the geometry of the volume (which I could provide in form of a gdml file), and it looks perfectly fine. On the other hand, I am not sure I understand the logic of the DistToOut call from a tube we are already kOutside of. In any case, I think having this kind of loops without any protection against infinite iterations is a major flaw. I would like to ask you (at the very least) to throw an exception and return 0 (or pushing by a few tolerances, or whatever solution you think is more sensible). Of course, I would be happy to test any patch or try any workaround you suggest. Thanks for your help, Andrea.