1 | #include <stdio.h> |
---|
2 | #include <string.h> |
---|
3 | #include <stdlib.h> |
---|
4 | #include <vector> |
---|
5 | #include <cmath> |
---|
6 | |
---|
7 | |
---|
8 | static double EPSILON=1.e-10; |
---|
9 | |
---|
10 | // facteur c/ 360. pour calculer (c dphi) / (360.freq) |
---|
11 | // avec freq en Mhz et dphi en degres et résultat en cm: |
---|
12 | |
---|
13 | static double FACTEUR = 83.3333; // ameliorer la precision |
---|
14 | |
---|
15 | // contrairement a ce qu'indique la notice, dans parmdesz, les xp et yp |
---|
16 | // sont donnes en radians |
---|
17 | static double uniteAngle = 1.0e+3; // passage de radians en milliradians |
---|
18 | static double uniteLong = 1.0; // on prend l'unite de long. TRANSPORT : cm |
---|
19 | |
---|
20 | |
---|
21 | static double EMeV = 0.511; |
---|
22 | struct particle |
---|
23 | { |
---|
24 | float xx, xxp, begamx,yy,yyp,begamy,z, begamz,phi,wz; |
---|
25 | float phi0, ksi1,ksi2,ksi3; |
---|
26 | int ne,np,ngood,npart; |
---|
27 | |
---|
28 | int readFromFile(FILE* fp) |
---|
29 | { |
---|
30 | return fscanf(fp, " %e %e %e %e %e %e %e %e %e %e %d %d %d %d %e %e %e %e \n", &xx, &xxp, &begamx,&yy,&yyp,&begamy,&z, &begamz,&phi,&wz,&ne,&np,&ngood,&npart,&phi0, &ksi1,&ksi2,&ksi3); |
---|
31 | } |
---|
32 | void imprim() |
---|
33 | { |
---|
34 | printf( " %e %e %e %e %e %e %e %e %e %e %d %d %d, %d %e %e %e %e \n", xx, xxp, begamx,yy,yyp,begamy,z, begamz,phi,wz,ne,np,ngood,npart,phi0, ksi1,ksi2,ksi3); |
---|
35 | } |
---|
36 | }; |
---|
37 | |
---|
38 | |
---|
39 | double readFreqFromParmout(FILE* fp) |
---|
40 | { |
---|
41 | float freq=0.0; |
---|
42 | char chaine[80]; |
---|
43 | while ( fscanf(fp, " %s ", chaine) > 0 ) |
---|
44 | { |
---|
45 | if ( strcmp(chaine, "frequency") == 0 ) |
---|
46 | { |
---|
47 | printf(" chaone lue : %s \n", chaine); |
---|
48 | fscanf(fp, " %s ", chaine); |
---|
49 | fscanf(fp, " %e ", &freq); |
---|
50 | break; |
---|
51 | } |
---|
52 | } |
---|
53 | return freq; |
---|
54 | }; |
---|
55 | |
---|
56 | |
---|
57 | int main(int argc,char* argv[]) |
---|
58 | { |
---|
59 | |
---|
60 | double freqRef; |
---|
61 | FILE* filefais; |
---|
62 | FILE* fileParmout; |
---|
63 | FILE* fileFaisTrans; |
---|
64 | std::vector<struct particle> faisceau; |
---|
65 | int k; |
---|
66 | |
---|
67 | if (argc!=2 ) { |
---|
68 | printf(" usage : recupFaisceau <nom du probleme parmela> \n"); |
---|
69 | exit(0); |
---|
70 | } |
---|
71 | |
---|
72 | printf(" ATTENTION : PROGRAMME recupFaisceauParm NON VERIFIE PRECISEMENT \n"); |
---|
73 | |
---|
74 | |
---|
75 | // printf(" frequence de reference %e Mhz \n", freqRef); |
---|
76 | char* nomParm = argv[1]; |
---|
77 | char nomParmout[132]; |
---|
78 | strcpy(nomParmout, nomParm); |
---|
79 | strcat(nomParmout, ".out"); |
---|
80 | |
---|
81 | printf(" probleme parmela : %s \n", argv[1]); |
---|
82 | // |
---|
83 | fileParmout = fopen(nomParmout, "r"); |
---|
84 | if ( fileParmout == (FILE*)0 ) |
---|
85 | { |
---|
86 | printf(" erreur a l'ouverture du fichier %s \n", nomParmout); |
---|
87 | exit(0); |
---|
88 | } |
---|
89 | else printf(" ouverture du fichier %s \n", nomParmout ); |
---|
90 | freqRef= readFreqFromParmout(fileParmout); |
---|
91 | fclose(fileParmout); |
---|
92 | // |
---|
93 | printf(" frequence relue : %e \n", freqRef); |
---|
94 | char nomfilefais[132]; |
---|
95 | strcpy(nomfilefais, nomParm); |
---|
96 | strcat(nomfilefais, ".desz"); |
---|
97 | |
---|
98 | filefais = fopen(nomfilefais, "r"); |
---|
99 | |
---|
100 | if ( filefais == (FILE*)0 ) |
---|
101 | { |
---|
102 | printf(" erreur a l'ouverture du fichier %s \n",nomfilefais ); |
---|
103 | exit(0); |
---|
104 | } |
---|
105 | else printf(" ouverture du fichier %s \n", nomfilefais); |
---|
106 | |
---|
107 | fileFaisTrans = fopen("faisceauTransp.dat", "w"); |
---|
108 | if ( fileFaisTrans == (FILE*)0 ) |
---|
109 | { |
---|
110 | printf(" erreur a l'ouverture du fichier faisceauTransp.dat \n"); |
---|
111 | exit(0); |
---|
112 | } |
---|
113 | |
---|
114 | struct particle partic; |
---|
115 | struct particle* partRefPtr=NULL; |
---|
116 | int numElem=-1; |
---|
117 | while( partic.readFromFile(filefais) > 0 ) |
---|
118 | |
---|
119 | { |
---|
120 | // on ne va conserver que les resultats a la sortie du dernier element |
---|
121 | if ( partic.ne != numElem ) |
---|
122 | { |
---|
123 | faisceau.clear(); |
---|
124 | partRefPtr=NULL; |
---|
125 | numElem = partic.ne; |
---|
126 | } |
---|
127 | faisceau.push_back(partic); |
---|
128 | if ( partic.np == 1 ) |
---|
129 | { |
---|
130 | if ( fabs(partic.xx) > EPSILON || fabs(partic.yy) > EPSILON || fabs(partic.xxp) > EPSILON || fabs(partic.yyp) > EPSILON) |
---|
131 | { |
---|
132 | printf(" ATTENTION part. reference douteuse \n"); |
---|
133 | partic.imprim(); |
---|
134 | } |
---|
135 | partRefPtr=&faisceau.back(); |
---|
136 | // printf("part. reference \n"); |
---|
137 | // partRefPtr->imprim(); |
---|
138 | } |
---|
139 | } |
---|
140 | // for (int k=0; k < faisceau.size(); k++) faisceau[k].imprim(); |
---|
141 | printf("dernier element %d \n", numElem); |
---|
142 | // pour preserver l'avenir, on definit un centroid (cf. TRANSPORT) |
---|
143 | double xg=0.; |
---|
144 | double xpg=0.; |
---|
145 | double yg=0.; |
---|
146 | double ypg=0.; |
---|
147 | double lg=0.; |
---|
148 | double delg=0.; |
---|
149 | |
---|
150 | // sigmas TRANSPORT |
---|
151 | double sig11=0.0; |
---|
152 | double sig22=0.0; |
---|
153 | double sig33=0.0; |
---|
154 | double sig44=0.0; |
---|
155 | double sig55=0.0; |
---|
156 | double sig66=0.0; |
---|
157 | |
---|
158 | // correlations |
---|
159 | double r21=0.0; |
---|
160 | double r31=0.0; |
---|
161 | double r32=0.0; |
---|
162 | double r41=0.0; |
---|
163 | double r42=0.0; |
---|
164 | double r43=0.0; |
---|
165 | double r51=0.0; |
---|
166 | double r52=0.0; |
---|
167 | double r53=0.0; |
---|
168 | double r54=0.0; |
---|
169 | double r61=0.0; |
---|
170 | double r62=0.0; |
---|
171 | double r63=0.0; |
---|
172 | double r64=0.0; |
---|
173 | double r65=0.0; |
---|
174 | |
---|
175 | |
---|
176 | double x,xp,y,yp,ll,del; |
---|
177 | // unites : TRANSPORT : cm, mrad, % |
---|
178 | // PARMELA : phase en degrés |
---|
179 | |
---|
180 | |
---|
181 | double gref = partRefPtr->wz/EMeV; |
---|
182 | double PrefMeVsc = sqrt( gref*(gref+2) ); |
---|
183 | |
---|
184 | |
---|
185 | for (int k=0; k < faisceau.size(); k++) |
---|
186 | // for (int k=0; k < 10; k++) |
---|
187 | { |
---|
188 | x=faisceau[k].xx; |
---|
189 | xp=faisceau[k].xxp; |
---|
190 | y=faisceau[k].yy; |
---|
191 | yp=faisceau[k].yyp; |
---|
192 | // dephasage par rapport a la reference |
---|
193 | double dephas = faisceau[k].phi - partRefPtr->phi; // degrés |
---|
194 | double g = faisceau[k].wz/EMeV; |
---|
195 | double betaz = faisceau[k].begamz/(g+1.0); |
---|
196 | // printf( " betaz= %e \n", betaz); |
---|
197 | // si j'ai bien compris, une phase positive indique une particule en avance dans l'espace |
---|
198 | double deltaz = FACTEUR * betaz * dephas / freqRef; |
---|
199 | // printf(" correction dz %e \n", deltaz); |
---|
200 | ll = deltaz; |
---|
201 | // printf("ancien x %e correction %e \n", x, xp * deltaz); |
---|
202 | x += xp * deltaz; |
---|
203 | faisceau[k].xx = x; |
---|
204 | // printf("ancien y %e correction %e \n", y, yp * deltaz); |
---|
205 | y += yp * deltaz; |
---|
206 | faisceau[k].yy =y; |
---|
207 | double PMeVsc = sqrt( g*(g+2) ); |
---|
208 | del = 100.0 * ( PMeVsc - PrefMeVsc ) / PrefMeVsc ; // en % |
---|
209 | // printf("defaut energie %e \n", del); |
---|
210 | sig11 += ( x - xg ) * ( x - xg ); |
---|
211 | sig22 += ( xp - xpg ) * ( xp - xpg ); |
---|
212 | sig33 += ( y - yg ) * ( y - yg ); |
---|
213 | sig44 += ( yp - ypg ) * ( yp - ypg ); |
---|
214 | sig55 += uniteLong*uniteLong*( ll - lg ) * ( ll - lg ); |
---|
215 | sig66 += ( del - delg ) * ( del - delg ); |
---|
216 | // |
---|
217 | |
---|
218 | |
---|
219 | r21 += ( xp - xpg ) * ( x - xg ); |
---|
220 | r31 += ( y - yg ) * ( x - xg ); |
---|
221 | r32 += ( y - yg ) * ( xp - xpg ); |
---|
222 | r41 += ( yp - ypg ) * ( x - xg ); |
---|
223 | r42 += ( yp - ypg ) * ( xp - xpg ); |
---|
224 | r43 += ( yp - ypg ) * ( y - yg ); |
---|
225 | r51 += ( ll - lg ) * ( x - xg ); |
---|
226 | r52 += ( ll - lg ) * ( xp - xpg ); |
---|
227 | r53 += ( ll - lg ) * ( y - yg ); |
---|
228 | r54 += ( ll - lg ) * ( yp - ypg ); |
---|
229 | r61 += ( del - delg ) * ( x - xg ); |
---|
230 | r62 += ( del - delg ) * ( xp - xpg ); |
---|
231 | r63 += ( del - delg ) * ( y - yg ); |
---|
232 | r64 += ( del - delg ) * ( yp - ypg ); |
---|
233 | r65 += ( del - delg ) * ( ll - lg ); |
---|
234 | |
---|
235 | } |
---|
236 | double facmoy = 1.0/double( faisceau.size() ); |
---|
237 | sig11 = sqrt(sig11*facmoy); |
---|
238 | sig22 = sqrt(sig22*facmoy); |
---|
239 | sig33 = sqrt(sig33*facmoy); |
---|
240 | sig44 = sqrt(sig44*facmoy); |
---|
241 | sig55 = sqrt(sig55*facmoy); |
---|
242 | sig66 = sqrt(sig66*facmoy); |
---|
243 | |
---|
244 | |
---|
245 | |
---|
246 | |
---|
247 | // |
---|
248 | |
---|
249 | r21= r21*facmoy/(sig22*sig11); |
---|
250 | r31= r31*facmoy/(sig33*sig11); |
---|
251 | r32= r32*facmoy/(sig33*sig22); |
---|
252 | r41= r41*facmoy/(sig44*sig11); |
---|
253 | r42= r42*facmoy/(sig44*sig22); |
---|
254 | r43= r43*facmoy/(sig44*sig33); |
---|
255 | r51= r51*facmoy/(sig55*sig11); |
---|
256 | r52= r52*facmoy/(sig55*sig22); |
---|
257 | r53= r53*facmoy/(sig55*sig33); |
---|
258 | r54= r54*facmoy/(sig55*sig44); |
---|
259 | r61= r61*facmoy/(sig66*sig11); |
---|
260 | r62= r62*facmoy/(sig66*sig22); |
---|
261 | r63= r63*facmoy/(sig66*sig33); |
---|
262 | r64= r64*facmoy/(sig66*sig44); |
---|
263 | r65= r65*facmoy/(sig66*sig55); |
---|
264 | |
---|
265 | // passage aux unites de sortie |
---|
266 | |
---|
267 | sig11 *= uniteLong; |
---|
268 | sig22 *= uniteAngle; |
---|
269 | sig33 *= uniteLong; |
---|
270 | sig44 *= uniteAngle; |
---|
271 | sig55 *= uniteLong; |
---|
272 | |
---|
273 | int dummy=0; |
---|
274 | fprintf(fileFaisTrans, " titre du probleme : %s \n", nomParm); |
---|
275 | fprintf(fileFaisTrans, " %d \n", dummy); |
---|
276 | fprintf(fileFaisTrans, " UTRANS \n"); |
---|
277 | |
---|
278 | |
---|
279 | fprintf(fileFaisTrans, " BEAM, X=%e, XP=%e, Y=%e, YP=%e, & \n", sig11, sig22, sig33, sig44 ); |
---|
280 | fprintf(fileFaisTrans, " DL=%e, DEL=%e, P0=%e ; \n", sig55, sig66, 1.0e-3*EMeV*PrefMeVsc ); |
---|
281 | |
---|
282 | |
---|
283 | fprintf(fileFaisTrans, " CORR, C21=%e, C31=%e, C32=%e, & \n", r21, r31, r32 ); |
---|
284 | fprintf(fileFaisTrans, " C41=%e, C42=%e, C43=%e, C51=%e, & \n", r41, r42, r43, r51 ); |
---|
285 | fprintf(fileFaisTrans, " C52=%e, C53=%e, C54=%e, C61=%e, & \n", r52, r53, r54, r61 ); |
---|
286 | fprintf(fileFaisTrans, " C62=%e, C63=%e, C64=%e, C65=%e ; \n", r62, r63, r64, r65 ); |
---|
287 | |
---|
288 | |
---|
289 | printf(" le faisceau TRANSPORT est cree sur le fichier : faisceauTransp.dat \n"); |
---|
290 | |
---|
291 | |
---|
292 | |
---|
293 | fclose(filefais); |
---|
294 | fclose(fileFaisTrans); |
---|
295 | } |
---|