Problème d'exécution : Un cout aux pouvoirs inatendus.
Bonjour à tous,
Je me permets de créer un nouveau sujet pour vous soumettre un problème qui me laisse perplexe depuis quelques jours maintenant. Comme le titre l'indique mon problème vient du fait que mon code ne s'exécute correctement que lorsque je mets un cout de certaines variables. (voir détails plus bas)
J'ai eu beau chercher un certain temps je n'ai jamais vu ce genre d'erreur et je ne sais pas trop quoi faire la résoudre.
Je présente rapidement le contexte afin de vous expliquer un peu le but du code : j’essaye de simuler l'interaction d'une particule avec un détecteur de la forme d'un parallélépipède rectangle.
Le code qui suit a pour fonction de déterminer la distance parcourue par une particule dans un parallélépipède rectangle.
Il est loin d'être optimisé au niveau algorithmique mais avant toute chose j'aimerai comprendre mon problème pour savoir si je dois tout recoder ou non.
Le programme prend comme input :
- un vecteur de trois composantes : Celui-ci représente le vecteur directeur du mouvement de la particule.
- Un point dans l'espace : Celui-ci représente la position du centre de notre détecteur dans l'espace.
- Une série de trois longueur : celles-ci sont les dimensions de notre détecteur dans les trois directions de l'espace.
Le détecteur est supposé faire face à la source d'émission des particules (placée à l'origine du repère). Traditionnellement en physique les particules se déplacent suivant l'axe Z
On peut donc définir une face d’entrée au détecteur qui dans notre cas sera celle contenue dans un plan z = constante.
Voici la façon dont je procède (voir code ci-dessus) : je teste si le vecteur entre bien par la face d'entrée, et si c'est le cas je calcule le point d'intersection entre le vecteur directeur et le plan d'entrée.
Si ce n'est pas le cas je m'arrête là, dans le cas contraire je cherche le point de sortie. Tout d'abord dans la face arrière puis dans les autres face le cas échéant.
A partir de ces deux points j'en déduit facilement la distance parcouru par la particule dans le détecteurs.
Code:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
|
double MakeDetectionEvent(vector<double> &Vector, const vector<double> &DimensionsOfDetector, const vector<double> &CenterPositionOfDetector) {
//Getting direction vector
double Norme = sqrt(Vector.at(0) * Vector.at(0) + Vector.at(1) * Vector.at(1) + Vector.at(2) * Vector.at(2));
for (int i = 0; i<3; i++)
Vector.at(i) = Vector.at(i) / Norme;
//Test if the particle pass through the front of the Detector
bool Enter = 0;
bool Leave = 0;
vector<double> FirstIntersection;
vector<double> SecondIntersection;
double Distance = -1;
//Looking for the entrace of the particle over all detectors
FirstIntersection = CalculateIntersection(Vector,DimensionsOfDetector,CenterPositionOfDetector,4);
Enter = IsInside(DimensionsOfDetector,CenterPositionOfDetector,FirstIntersection);
if (Enter == 0)
cout << "Particle does not enter in the detector " << endl;
else {
//Looking if the Particle leaving by the back side of the detector
SecondIntersection = CalculateIntersection(Vector,DimensionsOfDetector,CenterPositionOfDetector,5);
Leave = IsInside(DimensionsOfDetector,CenterPositionOfDetector,SecondIntersection);
if (Leave == 1) {
cout << "Particle leaves the detector by the back side" << endl;
} else {
//if not look for the other sides
for (int j=0;j<4;j++) {
SecondIntersection = CalculateIntersection(Vector,DimensionsOfDetector,CenterPositionOfDetector,j);
Leave = IsInside(DimensionsOfDetector,CenterPositionOfDetector,SecondIntersection);
if (Leave == 1)
break;
else
cout << "Paricle doesn't cross the " << j << "th plane" << endl;
}
}
if (Leave == 0)
cout << "MakeDetectionEvent : Error ! Particle goes in the detector but never goes out" << endl;
else
Distance = CalculateDistance(FirstIntersection, SecondIntersection);
}
return Distance;
} |
Code:
1
2
3
4
5
6
7
8
9
10
11
|
double CalculateDistance(const std::vector<double> &Entrance, const std::vector<double> &Exit) {
double x = Entrance.at(0) - Exit.at(0);
double y = Entrance.at(1) - Exit.at(1);
double z = Entrance.at(2) - Exit.at(2);
double Distance = sqrt(x*x + y*y + z*z);
return Distance;
} |
Le calcul des points d'intersections se fait suivant l'indice du plan (0 pour gauche, 1 pour droite, 2 pour bas, 3 pour haut, 4 pour avant et 5 pour arrière ) grâce à la fonction suivante.
Code:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
|
vector<double> CalculateIntersection(const vector<double> &Vector, const vector<double> Dimensions, const vector<double> CenterPosition,const int PlaneIndex) {
vector<double> Intersection(3);
double Lambda;
switch(PlaneIndex) {
case 0 :
Intersection.at(0) = CenterPosition.at(0) - Dimensions.at(0) / 2;
Lambda = Intersection.at(0) / Vector.at(0);
Intersection.at(1) = Vector.at(1) * Lambda;
Intersection.at(2) = Vector.at(2) * Lambda;
break;
case 1 :
Intersection.at(0) = CenterPosition.at(0) + Dimensions.at(0) / 2;
Lambda = Intersection.at(0) / Vector.at(0);
Intersection.at(1) = Vector.at(1) * Lambda;
Intersection.at(2) = Vector.at(2) * Lambda;
break;
case 2 :
Intersection.at(1) = CenterPosition.at(1) - Dimensions.at(1) / 2;
Lambda = Intersection.at(1) / Vector.at(1);
Intersection.at(0) = Vector.at(0) * Lambda;
Intersection.at(2) = Vector.at(2) * Lambda;
break;
case 3 :
Intersection.at(1) = CenterPosition.at(1) + Dimensions.at(1) / 2;
Lambda = Intersection.at(1) / Vector.at(1);
Intersection.at(0) = Vector.at(0) * Lambda;
Intersection.at(2) = Vector.at(2) * Lambda;
break;
case 4 :
Intersection.at(2) = CenterPosition.at(2) - Dimensions.at(2) / 2;
Lambda = Intersection.at(2) / Vector.at(2);
Intersection.at(0) = Vector.at(0) * Lambda;
Intersection.at(1) = Vector.at(1) * Lambda;
break;
case 5 :
Intersection.at(2) = CenterPosition.at(2) + Dimensions.at(2) / 2;
Lambda = Intersection.at(2) / Vector.at(2);
Intersection.at(0) = Vector.at(0) * Lambda;
Intersection.at(1) = Vector.at(1) * Lambda;
break;
default :
cout << "CalculateIntersection : Error ! PlaneIndex is not well defined" << endl;
exit(1);
}
cout << "Intersection = " << Intersection.at(0) << " " << Intersection.at(1) << " " << Intersection.at(2) << endl;
return Intersection;
} |
Enfin pour tester si ce point d'intersection se trouve bien dans le volume voulu j'utilise la fonction suivante :
Code:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
| bool IsInside(const vector<double> &DimensionsOfDetector, const vector<double> &CenterPositionOfDetector, const std::vector<double> &Intersection) {
bool Inside = 1;
vector<double> Borders(6);
Borders.at(0) = CenterPositionOfDetector.at(0) - DimensionsOfDetector.at(0) / 2;
Borders.at(1) = CenterPositionOfDetector.at(0) + DimensionsOfDetector.at(0) / 2;
Borders.at(2) = CenterPositionOfDetector.at(1) - DimensionsOfDetector.at(1) / 2;
Borders.at(3) = CenterPositionOfDetector.at(1) + DimensionsOfDetector.at(1) / 2;
Borders.at(4) = CenterPositionOfDetector.at(2) - DimensionsOfDetector.at(2) / 2;
Borders.at(5) = CenterPositionOfDetector.at(2) + DimensionsOfDetector.at(2) / 2;
cout << "42" << endl;
// cout << "Intersection = " << Intersection.at(0) << " " << Intersection.at(1) << " " << Intersection.at(2) << endl;
if (Intersection.at(0) > Borders.at(1))
Inside = 0;
if(Intersection.at(0) < Borders.at(0))
Inside = 0;
if (Intersection.at(1) > Borders.at(3))
Inside = 0;
if(Intersection.at(1) < Borders.at(2))
Inside = 0;
if (Intersection.at(2) > Borders.at(5))
Inside = 0;
if(Intersection.at(2) < Borders.at(4))
Inside = 0;
return Inside;
} |
Et voici le main qui appelle le tout :
Code:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
| int main() {
cout.precision(12);
vector<double> CenterPositions;
CenterPositions.push_back(0);
CenterPositions.push_back(0);
CenterPositions.push_back(5.1);
vector<double> Dimensions;
Dimensions.push_back(0.20);
Dimensions.push_back(2);
Dimensions.push_back(0.20);
vector <double> VectorTest(3);
VectorTest.at(0) = 0;
VectorTest.at(1) = 0;
VectorTest.at(2) = 1;
double Distance = MakeDetectionEvent(VectorTest,Dimensions,CenterPositions);
cout << "Distance = " << Distance << endl;
return 0;
} |
Alors où ce trouve le problème dans tout ça ?
Et bien en centrant mon détecteur sur l'axe Z et en utilisant un vecteur directeur uniquement suivant l'axe des Z la valeur renvoyé par IsInside à l'entrée est false.
J'ai donc voulu faire du pas à pas en affichant les valeurs à la mains : Et la surprise le code s'exécute correctement si l'on rajoute la ligne de code 13 dans IsInside...
cout << "Intersection = " << Intersection.at(0) << " " << Intersection.at(1) << " " << Intersection.at(2) << endl;Très bien sauf que... Celle-ci n'est sensé n'avoir aucun impacte sur la suite... :aie:
Et les valeurs des intersections sont systématiquement bonnes dans les deux cas...
Pour terminer et en rajouter quant à mon incompréhension : j'ai essayé de trouver le problème avec gdb sans succès : je n'ai rien vu d'anormal. (je débute néanmoins avec gdb donc j'ai peut être raté quelque chose).
Et en utilisant valgrind j'ai eu la surprise de constater que le programme s’exécutait correctement avec ou sans le fameux commentaire ligne 13.
Du coup je commence à me demander si tout ça ne serait pas lien à des problèmes d'options de compilations ( optimisation en -O3 ?).
Au final je suis totalement perdu... J'espère que vous aurez de meilleurs idées que moi.
Merci d'avance de votre aide.
P.S. Le code fait parti d'un projet plus gros en C++ : néanmoins comme cette partie ne marchait pas je l'ai réécrite rapidement "en C". Cela explique l'utilisation des librairies std et iostream. Comme le problème n'est pas propre à l'encapsulation j'ai jugé qu'il serait plus adapté au forum C.
P.S.2 Je n'ai jamais reçu de formation en informatique / algorithmique / langage C, n'hésitez pas à me faire des remarques si vous voyez des aberrations...
Edit: J'ai oublié de préciser je compile avec g++ sur ubuntu.
j'utilise les options suivantes : g++ -O3 -Wall -g
