G4LambertianRand(G4ThreeVector* normal) function falls into the infinite loop when its argument is zero vector. G4LambertianRand is called in the G4OpBounadryProcess::DoReflection(). It uses theGlobalNormal as argument without checking if it valid. Code to get normal=(0., 0., 0.): G4Tubs* blk = new G4Tubs("blk",2.249740e+03, 2.252740e+03, 4.990000e02, 0.000000e+00, 2.609116e-01); G4ThreeVector localpoint(-1.59250955522304866463e+03, 1.60239571861107833683e+03, 1.55565646370305239543e+03); G4ThreeVector dir(-2.97231637897997502673e-01, -9.46902667985754731284e-01, -1.22591560859357129321e-01); bool valid; G4ThreeVector norm; double dist=blk->DistanceToOut(localpoint, dir, true, &valid, &norm); printf("dist=%20f, norm(%20f, %20f, %20f), valid=%i\n", dist,norm.x(),norm.y(),norm.z(),valid);
G4OpBoundary::PostStepDoIt now throws an G4Exception EventMustBeAborted when an invalid surface normal is returnd. The underlying problem has to do with G4Tubs and I will reassign.
A geometry defined that way with hard-coded numbers is ill by definition, it can lead to undefined behavior due to precision truncation. Please, provide a case reproducing the problem with precise parameters.
Dear Maxim, DistanceToOut(localpoint, dir, true, &valid, &norm) can return SurfaceNormal of (0,0,0), it happends for all non-convex surfaces, like the Inner Radius Surface for G4Tubs. In this case a boolean flag for validnormal set to false and SurfaceNormal() is called for SurfaceNormal calculations. To investigate the problem more information about your geometrical setup is needed, like the placement of the G4Tubs in the mother volume, the mother volume information(solid,placement,...). Best Regards, Tatiana Nikitina
Dear Gabriele, numbers actually are not hard-coded. The geometry is described via xml in Gaudi framework. It's rather complicated, so I've tried to extract only the volume which causes this infinite loop. The chosen precision is enough to reproduce zero-normal. Dear Tatiana, I didn't mean that there are errors in geometry. The geometry information is added only to show the zero-normal. As I've said it is not handy to extract whole information about this volume. As for the SurfaceNormal() --- I didn't notice it. I'll look and if more information about geometry is still needed i'll send it. regards, Maxim Gonchar