1 | #include <iostream> |
---|
2 | #include <fstream> |
---|
3 | #include <cmath> |
---|
4 | #include <iomanip> |
---|
5 | #include <string> |
---|
6 | #include "simcrys.h" |
---|
7 | |
---|
8 | |
---|
9 | using namespace std; |
---|
10 | |
---|
11 | |
---|
12 | |
---|
13 | SimCrys::SimCrys(Crystal crys, Partcrys part) |
---|
14 | : crys(crys), part(part) //,moy(moy) //POUR faire varier la moyenne de la gaussienne aussi !!!!! |
---|
15 | { |
---|
16 | } |
---|
17 | |
---|
18 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
19 | ////////////////////////// FUNCTION BASES THAT MAKE SOME BASIC CALCULATION |
---|
20 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
21 | |
---|
22 | void SimCrys::bases(SimCrys& sim) |
---|
23 | { |
---|
24 | if ( (crys.miscut < 0 ) && ( part.x > 0) && (part.x < -crys.Cry_length * tan(crys.miscut))) { |
---|
25 | crys.L_chan = -part.x * tan(crys.miscut); |
---|
26 | crys.tchan = crys.L_chan / crys.Rcurv; |
---|
27 | } |
---|
28 | part.xp_rel = part.xp - crys.miscut; |
---|
29 | } |
---|
30 | |
---|
31 | |
---|
32 | |
---|
33 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
34 | ////////////////////////FUNCTIN CROSS TO KNOW IF THE PARTICULE CROSS THE CRYSTAL |
---|
35 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
36 | |
---|
37 | |
---|
38 | |
---|
39 | |
---|
40 | bool SimCrys::cross(SimCrys& sim) |
---|
41 | { |
---|
42 | bool out(true); |
---|
43 | if ((part.y < crys.ymin) || (part.y > crys.ymax) || (part.x > crys.xmax)) { |
---|
44 | //cout << "OUT of the Crystal !!" << part.x << " , "<<part.xp << endl; |
---|
45 | out = false; |
---|
46 | part.x = part.x + part.xp * crys.s_length; |
---|
47 | |
---|
48 | part.y = part.y + part.yp * crys.s_length; |
---|
49 | |
---|
50 | part.nout = part.nout + 1; /////////////COUNT///////////// |
---|
51 | part.effet = 1; |
---|
52 | //cout<<" BLAAAAAAAAa"<<endl; |
---|
53 | |
---|
54 | } |
---|
55 | |
---|
56 | return out; |
---|
57 | } |
---|
58 | |
---|
59 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
60 | ////////////////////////////FUNCTION LAYER THAT DETERMIN IF THE PARTICULE IS IN THE AMORPHOUS LAYER///////////////////////////// |
---|
61 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
62 | |
---|
63 | bool SimCrys::layer(SimCrys& sim) |
---|
64 | { |
---|
65 | //cout<<"La fonction layer est utilisee>>>>>>>>>>>>>>"<<endl; |
---|
66 | bool out(true); |
---|
67 | if ((part.x < crys.Alayer) || ((part.y - crys.ymin) < crys.Alayer) || ((crys.ymax - part.y) < crys.Alayer)) { |
---|
68 | //cout<<"Amorphous layer! " << part.x <<" , " <<part.xp<<endl; |
---|
69 | double a_eq; |
---|
70 | double b_eq; |
---|
71 | double c_eq; |
---|
72 | double delta; //a, b, c, Delta of the second order eq. |
---|
73 | double x0; |
---|
74 | double y0; |
---|
75 | double length_xs; //!Amorphous length |
---|
76 | double length_ys; //!Amorphous length |
---|
77 | double am_length; //!Amorphous length |
---|
78 | out = false; |
---|
79 | x0 = part.x; |
---|
80 | y0 = part.y; |
---|
81 | a_eq = (1 + part.xp * part.xp); |
---|
82 | b_eq = (2 * (part.x) * (part.xp) - 2 * (part.xp) * crys.Rcurv); |
---|
83 | c_eq = (part.x) * part.x - 2 * abs(part.x) * crys.Rcurv; |
---|
84 | delta = b_eq * b_eq - 4 * a_eq * c_eq; |
---|
85 | //cout<<"a : " << a_eq<<endl; //POUR TESTER FONCTION |
---|
86 | //cout<<"b: " << b_eq <<endl; //POUR TESTER FONCTION |
---|
87 | //cout<<"c : " << c_eq <<endl; //POUR TESTER FONCTION |
---|
88 | //cout<<"Delta : " << delta <<endl; //POUR TESTER FONCTION |
---|
89 | part.s = (-b_eq + sqrt(delta)) / ((2 * a_eq)) ; //in [m] |
---|
90 | //cout<<"s : " << part.s <<endl; //POUR TESTER FONCTION |
---|
91 | if (part.s >= crys.s_length) { |
---|
92 | part.s = crys.s_length; |
---|
93 | } |
---|
94 | //cout<<"s : " << part.s <<endl; //POUR TESTER FONCTION |
---|
95 | part.x = part.xp * part.s + x0 ; //in [m] |
---|
96 | //cout<<"x0 : " << x0 <<endl; //POUR TESTER FONCTION |
---|
97 | //cout<<"x : " << part.x <<endl; //POUR TESTER FONCTION |
---|
98 | |
---|
99 | length_xs = sqrt((part.x - x0) * (part.x - x0) + part.s * part.s); |
---|
100 | //cout<<"Length_XS : " << length_xs <<endl; //POUR TESTER FONCTION |
---|
101 | |
---|
102 | if ((((part.yp >= 0) && ((part.y + part.yp * part.s) <= crys.ymax))) || (((part.yp < 0) && ((part.y + part.yp * part.s) >= crys.ymin)))) { |
---|
103 | length_ys = part.yp * length_xs; |
---|
104 | } else { |
---|
105 | part.s = (crys.ymax - part.y) / part.yp; |
---|
106 | length_ys = sqrt((crys.ymax - part.y) * (crys.ymax - part.y) + part.s * part.s); |
---|
107 | part.x = x0 + part.xp * part.s; |
---|
108 | length_xs = sqrt((part.x - x0) * (part.x - x0) + part.s * part.s); |
---|
109 | } |
---|
110 | |
---|
111 | am_length = sqrt(length_xs * length_xs + length_ys * length_ys); |
---|
112 | //cout << "AM_Length :" << am_length << endl; // TESTER FONCTION |
---|
113 | part.s = part.s / 2; |
---|
114 | part.x = x0 + part.xp * part.s; |
---|
115 | part.y = y0 + part.yp * part.s; |
---|
116 | //cout<<" am_length "<<am_length<<endl; |
---|
117 | move_am(crys.IS, crys.NAM, am_length, crys.DES[crys.IS], crys.DLYI[crys.IS], crys.DLRI[crys.IS], part.xp, part.yp, part.PC); //We call the move_am function |
---|
118 | //cout<<"AM LAYER !!!! "<<part.PC<<endl; |
---|
119 | //cout<< "Amorphous"<< endl; |
---|
120 | part.x = part.x + part.xp * (crys.Cry_length - part.s); |
---|
121 | part.y = part.y + part.yp * (crys.Cry_length - part.s); |
---|
122 | part.namor = part.namor + 1; /////////////COUNT///////////// |
---|
123 | part.effet = 2; |
---|
124 | } else if ((part.x > (crys.xmax - crys.Alayer)) && (part.x < crys.xmax)) { |
---|
125 | out = false; |
---|
126 | |
---|
127 | move_am(crys.IS, crys.NAM, crys.s_length, crys.DES[crys.IS], crys.DLYI[crys.IS], crys.DLRI[crys.IS], part.xp, part.yp, part.PC); //We call the move_am function |
---|
128 | |
---|
129 | part.namor = part.namor + 1; /////////////COUNT///////////// |
---|
130 | part.effet = 2; |
---|
131 | //cout<<"Fix HERE"<<endl; |
---|
132 | } |
---|
133 | return out; |
---|
134 | } |
---|
135 | |
---|
136 | |
---|
137 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
138 | //SI LES DEUX FONCTIONS D AVT RETURNE TRUE, ALORS ON CALCULE LES PARA CRITIQUE ET D AUTRE PARA (POINT C - E) |
---|
139 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
140 | |
---|
141 | |
---|
142 | |
---|
143 | void SimCrys::parameters(SimCrys& sim) |
---|
144 | { |
---|
145 | //THE FIRST PART IS FOR THE (110) ORIENTATION |
---|
146 | //cout << "The parameters function is use !!!!! "<< endl; |
---|
147 | double c_v1; |
---|
148 | double c_v2; |
---|
149 | |
---|
150 | |
---|
151 | crys.xpcrit0 = sqrt(2.10e-9 * crys.eUm[crys.IS] / part.PC); |
---|
152 | crys.Rcrit = part.PC / (2.e-6 * crys.eUm[crys.IS]) * crys.AI[crys.IS]; //if R>Rcritical=>no channeling is possible |
---|
153 | crys.ratio = crys.Rcurv / crys.Rcrit; |
---|
154 | crys.xpcrit = crys.xpcrit0 * (crys.Rcurv - crys.Rcrit) / crys.Rcurv; |
---|
155 | c_v1 = 1.7; |
---|
156 | c_v2 = -1.5; |
---|
157 | if (crys.ratio <= 1) { |
---|
158 | part.Ang_rms = c_v1 * 0.42 * crys.xpcrit0 * sin(1.4 * crys.ratio); //rms scattering |
---|
159 | part.Ang_avr = c_v2 * crys.xpcrit0 * 0.05 * crys.ratio; //average angle reflection |
---|
160 | part.Vcapt = 0.0; |
---|
161 | } else if (crys.ratio <= 3) { |
---|
162 | part.Ang_rms = c_v1 * 0.42 * crys.xpcrit0 * sin(1.571 * 0.3 * crys.ratio + 0.85); //rms scattering |
---|
163 | part.Ang_avr = c_v2 * crys.xpcrit0 * (0.1972 * crys.ratio - 0.1472); //average angle reflection |
---|
164 | part.Vcapt = 0.0007 * (crys.ratio - 0.7) / pow(part.PC, 0.2); //K=0.00070 is taken based on simulations using CATCH.f (V.Biryukov) |
---|
165 | } else { |
---|
166 | part.Ang_rms = c_v1 * crys.xpcrit0 * (1 / crys.ratio); |
---|
167 | part.Ang_avr = c_v2 * crys.xpcrit0 * (1 - 1.6667 / crys.ratio); |
---|
168 | part.Vcapt = 0.0007 * (crys.ratio - 0.7) / pow(part.PC, 0.2); |
---|
169 | } |
---|
170 | //cout <<crys.xpcrit<<endl; |
---|
171 | if (crys.C_orient == 2) { |
---|
172 | part.Ang_avr = part.Ang_avr * 0.93; |
---|
173 | part.Ang_rms = part.Ang_rms * 1.05; // FOR (111) ORIENTATION |
---|
174 | crys.xpcrit = crys.xpcrit * 0.98; |
---|
175 | } |
---|
176 | //cout <<crys.xpcrit<<endl; |
---|
177 | } |
---|
178 | |
---|
179 | |
---|
180 | |
---|
181 | |
---|
182 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
183 | //MAINTENANT ON FAIT UNE FONCTION POUR LA PARTIE I CF SHEMA!!! C EST LA PARTIE CHANNELING D OU LE NOM |
---|
184 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
185 | |
---|
186 | |
---|
187 | |
---|
188 | |
---|
189 | void SimCrys::channel(SimCrys& sim) |
---|
190 | { |
---|
191 | //cout << "La fonction channel est utilisee!!!! "<<endl; |
---|
192 | |
---|
193 | |
---|
194 | part.Chann = sqrt(crys.xpcrit * crys.xpcrit - part.xp_rel * part.xp_rel) / crys.xpcrit0; //probability of channeling/volume capture |
---|
195 | |
---|
196 | //cout <<"Channeling Prob :" << part.Chann <<endl; //POUR REGARDER LA VALEUR!!!!!!!!!!!!!! |
---|
197 | //cout <<"Channeling ANGLE :" << crys.tchan<<endl; //POUR REGARDER LA VALEUR!!!!!!!!!!!!!! |
---|
198 | |
---|
199 | |
---|
200 | double n_atom(0.1); //probability for entering channeling near atomic planes |
---|
201 | double Dxp(0); //change angle from channeling [mrad] |
---|
202 | double TLdech1(0); //tipical dechanneling length(1) [m] |
---|
203 | double tdech(0); |
---|
204 | double Ldech(0); //Dechanneling. length |
---|
205 | double Sdech(0); |
---|
206 | //cout<<"Chann : "<<part.Chann<<endl; |
---|
207 | if (random() <= part.Chann) { |
---|
208 | TLdech1 = 0.0005 * part.PC * pow(1 - 1 / crys.ratio, 2); // calculate tipical dechanneling length(m) |
---|
209 | if (random() <= n_atom) { |
---|
210 | TLdech1 = 0.000002 * part.PC * pow(1. - 1. / crys.ratio, 2); // dechanneling length (m) for short crystal for high amplitude particles |
---|
211 | } |
---|
212 | //cout<<"@ energy " << part.PC <<" dechanneling length "<< TLdech1 << endl; //FAUDRA OTER CELA!!!!!!!!!! |
---|
213 | part.Dechan = log(random()); |
---|
214 | //cout<<"dechannn " << part.Dechan <<endl; //FAUDRA OTER CELA!!!!!!!!!! |
---|
215 | Ldech = -TLdech1 * part.Dechan; //careful: the dechanneling lentgh is along the trajectory of the particle -not along the longitudinal coordinate... |
---|
216 | tdech = Ldech / crys.Rcurv; |
---|
217 | Sdech = Ldech * cos(crys.miscut + 0.5 * tdech); |
---|
218 | // 2 option if the particule channeling !!! |
---|
219 | // (1) dechanneling |
---|
220 | // (2) full channeling |
---|
221 | |
---|
222 | // (1) : |
---|
223 | |
---|
224 | if (Ldech < crys.L_chan) { |
---|
225 | //cout <<"Dechanneling ! "<< endl; |
---|
226 | Dxp = Ldech / crys.Rcurv; //change angle from channeling [mrad] |
---|
227 | |
---|
228 | part.ndechann = part.ndechann + 1; /////////////COUNT///////////// |
---|
229 | part.effet = 4; |
---|
230 | part.x = part.x + Ldech * (sin(0.5 * Dxp + crys.miscut)); // trajectory at channeling exit |
---|
231 | part.xp = part.xp + Dxp + (random_gauss() - 0.5) * crys.xpcrit * 0.5; |
---|
232 | part.y = part.y + part.yp * Sdech; |
---|
233 | |
---|
234 | part.x = part.x + 0.5 * (crys.s_length - Sdech) * part.xp; |
---|
235 | part.y = part.y + 0.5 * (crys.s_length - Sdech) * part.yp; |
---|
236 | |
---|
237 | |
---|
238 | move_am(crys.IS, crys.NAM, crys.s_length - Sdech, crys.DES[crys.IS], crys.DLYI[crys.IS], crys.DLRI[crys.IS], part.xp, part.yp, part.PC); //We call the move_am function |
---|
239 | |
---|
240 | part.PC = part.PC - 0.5 * crys.DES[crys.IS] * part.y; //energy loss to ionization [GeV] |
---|
241 | part.x = part.x + 0.5 * (crys.s_length - Sdech) * part.xp; |
---|
242 | part.y = part.y + 0.5 * (crys.s_length - Sdech) * part.yp; |
---|
243 | } |
---|
244 | |
---|
245 | // (2) : |
---|
246 | |
---|
247 | else { |
---|
248 | //cout <<"Channeling ! "<< endl; |
---|
249 | Dxp = crys.L_chan / crys.Rcurv; //change angle [mrad] |
---|
250 | part.x = part.x + crys.L_chan * (sin(0.5 * Dxp + crys.miscut)); //trajectory at channeling exit |
---|
251 | |
---|
252 | |
---|
253 | // Removed dependence on incoming angle according to discussion with I.Yazynin on 24/10/2012 |
---|
254 | //part.xp = part.xp + Dxp + (random_gauss()-0.5)*crys.xpcrit*0.5; // [mrad] |
---|
255 | part.xp = /*part.xp*/ + Dxp + (random_gauss() - 0.5) * crys.xpcrit * 0.5; // [mrad] |
---|
256 | |
---|
257 | part.y = part.y + crys.s_length * part.yp; |
---|
258 | part.PC = part.PC - 0.5 * crys.DES[crys.IS] * crys.Cry_length; // energy loss to ionization [GeV] |
---|
259 | part.nchann = part.nchann + 1; /////////////COUNT///////////// |
---|
260 | part.effet = 3; |
---|
261 | } |
---|
262 | } else { |
---|
263 | //cout <<"Volume Reflection ! "<< endl; |
---|
264 | part.nvrefl = part.nvrefl + 1; /////////////COUNT///////////// |
---|
265 | part.effet = 5; |
---|
266 | Dxp = 0.5 * (part.xp_rel / crys.xpcrit + 1) * part.Ang_avr; |
---|
267 | part.xp = part.xp + Dxp + part.Ang_rms * random_gauss(); |
---|
268 | /* |
---|
269 | next line new from sasha |
---|
270 | part.xp=part.xp+0.45*(part.xp/crys.xpcrit+1)*part.Ang_avr; |
---|
271 | valentina fix here |
---|
272 | */ |
---|
273 | part.x = part.x + 0.5 * crys.s_length * part.xp; |
---|
274 | part.y = part.y + 0.5 * crys.s_length * part.yp; |
---|
275 | |
---|
276 | move_am(crys.IS, crys.NAM, crys.s_length, crys.DES[crys.IS], crys.DLYI[crys.IS], crys.DLRI[crys.IS], part.xp , part.yp, part.PC); //We call the move_am function |
---|
277 | |
---|
278 | part.x = part.x + 0.5 * crys.s_length * part.xp; |
---|
279 | part.y = part.y + 0.5 * crys.s_length * part.yp; |
---|
280 | |
---|
281 | } |
---|
282 | } |
---|
283 | |
---|
284 | |
---|
285 | |
---|
286 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
287 | //MAINTENANT ON FAIT UNE FONCTION POUR LA PARTIE II CF SHEMA!!! C EST LA PARTIE VOLUME REFLEXION D OU LE NOM |
---|
288 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
289 | |
---|
290 | |
---|
291 | |
---|
292 | |
---|
293 | void SimCrys::reflection(SimCrys& sim) |
---|
294 | { |
---|
295 | //cout<<"The function reflection is used" <<endl; |
---|
296 | double Dxp(0); //change angle from channeling [mrad] |
---|
297 | double Lrefl(0); //distance of the reflection point [m] |
---|
298 | double Srefl(0); //distance of the reflection point [m] |
---|
299 | double TLdech2(0); //tipical dechanneling length(2) [m] |
---|
300 | double tdech(0); |
---|
301 | double Ldech(0); //Dechanneling. length |
---|
302 | double Sdech(0); |
---|
303 | double Red_S(0); //Reduced length (in case of dechanneling) |
---|
304 | double Rlength(0); //Reduced length |
---|
305 | |
---|
306 | Lrefl = (part.xp_rel) * crys.Rcurv; |
---|
307 | Srefl = sin(part.xp) * Lrefl; |
---|
308 | |
---|
309 | |
---|
310 | if ((Lrefl > 0) && (Lrefl < crys.Cry_length)) { // IE : If the VR point is inside the crystal!!!! |
---|
311 | |
---|
312 | //cout<<"VRCapt : "<<part.Vcapt<<endl; |
---|
313 | if ((random() > part.Vcapt) || (crys.ZN == 0)) { |
---|
314 | //cout<< "Volume Reflexion! "<<endl; |
---|
315 | part.nvrefl = part.nvrefl + 1; /////////////COUNT///////////// |
---|
316 | part.effet = 5; |
---|
317 | part.x = part.x + part.xp * Srefl; |
---|
318 | part.y = part.y + part.yp * Srefl; |
---|
319 | Dxp = part.Ang_avr; |
---|
320 | part.xp = part.xp + Dxp + part.Ang_rms * random_gauss(); |
---|
321 | part.x = part.x + 0.5 * part.xp * (crys.s_length - Srefl); |
---|
322 | part.y = part.y + 0.5 * part.yp * (crys.s_length - Srefl); |
---|
323 | |
---|
324 | if (crys.ZN > 0) { |
---|
325 | move_am(crys.IS, crys.NAM, crys.s_length - Srefl, crys.DES[crys.IS], crys.DLYI[crys.IS], crys.DLRI[crys.IS], part.xp , part.yp, part.PC); //We call the move_am function |
---|
326 | } |
---|
327 | |
---|
328 | part.x = part.x + 0.5 * part.xp * (crys.s_length - Srefl); |
---|
329 | part.y = part.y + 0.5 * part.yp * (crys.s_length - Srefl); |
---|
330 | } else { |
---|
331 | //cout<< "Volume Capture! "; |
---|
332 | part.nvcapt = part.nvcapt + 1; /////////////COUNT///////////// |
---|
333 | part.effet = 6; |
---|
334 | part.x = part.x + part.xp * Srefl; |
---|
335 | part.y = part.y + part.yp * Srefl; |
---|
336 | /* |
---|
337 | TLdech2= 0.00011*pow(part.PC,0.25)*pow((1.-1./crys.ratio),2) //dechanneling length [m] |
---|
338 | part.Dechan = log(1.-random()) |
---|
339 | Ldech = -TLdech2*part.Dechan |
---|
340 | next 2 lines new from sasha - different dechanneling probability |
---|
341 | */ |
---|
342 | TLdech2 = 0.01 * part.PC * pow((1. - 1. / crys.ratio), 2); // typical dechanneling length [m] |
---|
343 | Ldech = 0.005 * TLdech2 * pow((sqrt(0.01 - log(random())) - 0.1), 2); // (dechanneling length) [m] |
---|
344 | tdech = Ldech / crys.Rcurv; |
---|
345 | Sdech = Ldech * cos(part.xp + 0.5 * tdech); |
---|
346 | |
---|
347 | |
---|
348 | |
---|
349 | if (Ldech < (crys.Cry_length - Lrefl)) { // for dechanneling part |
---|
350 | //cout<< "Dechanneling! "; |
---|
351 | part.ndechann = part.ndechann + 1; /////////////COUNT///////////// |
---|
352 | part.effet = 4; |
---|
353 | Dxp = Ldech / crys.Rcurv; |
---|
354 | part.x = part.x + Ldech * (sin(0.5 * Dxp + part.xp)); //trajectory at channeling exit |
---|
355 | part.y = part.y + Sdech * part.yp; |
---|
356 | part.xp = part.xp + Dxp + (random_gauss() - 0.5) * crys.xpcrit * 0.5; |
---|
357 | Red_S = crys.s_length - Srefl - Sdech; |
---|
358 | // move in amorphous substance---------------------------- |
---|
359 | part.x = part.x + 0.5 * part.xp * Red_S; |
---|
360 | part.y = part.y + 0.5 * part.yp * Red_S; |
---|
361 | |
---|
362 | if (crys.ZN > 0) { |
---|
363 | //cout<<"Amorphous! "<<endl; /////////////COUNT///////////// |
---|
364 | move_am(crys.IS, crys.NAM, Red_S, crys.DES[crys.IS], crys.DLYI[crys.IS], crys.DLRI[crys.IS], part.xp , part.yp, part.PC); //We call the move_am function |
---|
365 | } |
---|
366 | |
---|
367 | part.x = part.x + 0.5 * part.xp * Red_S; |
---|
368 | part.y = part.y + 0.5 * part.yp * Red_S; |
---|
369 | //cout<<endl; |
---|
370 | } else { |
---|
371 | //cout<< "Volume Capture! "<<endl; |
---|
372 | //Pas besoin d additionner 1 dans le compte de vol capt car deja fait!!! |
---|
373 | Rlength = crys.Cry_length - Lrefl; |
---|
374 | crys.tchan = Rlength / crys.Rcurv; |
---|
375 | Red_S = Rlength * cos(part.xp + 0.5 * crys.tchan); |
---|
376 | Dxp = (crys.Cry_length - Lrefl) / crys.Rcurv; |
---|
377 | part.x = part.x + sin(0.5 * Dxp + part.xp) * Rlength; // trajectory at channeling exit |
---|
378 | part.y = part.y + Red_S * part.yp; |
---|
379 | part.xp = part.xp + Dxp + (random_gauss() - 0.5) * crys.xpcrit * 0.5; //[mrad] |
---|
380 | } |
---|
381 | } |
---|
382 | } |
---|
383 | // move in amorphous substance for big input angles--------------- |
---|
384 | else { |
---|
385 | //cout<<"Amorphous! "<<endl; |
---|
386 | part.namor = part.namor + 1; /////////////COUNT///////////// |
---|
387 | part.effet = 2; |
---|
388 | part.x = part.x + 0.5 * crys.s_length * part.xp; |
---|
389 | part.y = part.y + 0.5 * crys.s_length * part.yp; |
---|
390 | |
---|
391 | if (crys.ZN > 0) { |
---|
392 | move_am(crys.IS, crys.NAM, crys.s_length, crys.DES[crys.IS], crys.DLYI[crys.IS], crys.DLRI[crys.IS], part.xp , part.yp, part.PC); //We call the move_am function |
---|
393 | } |
---|
394 | |
---|
395 | part.x = part.x + 0.5 * crys.s_length * part.xp; |
---|
396 | part.y = part.y + 0.5 * crys.s_length * part.yp; |
---|
397 | } |
---|
398 | |
---|
399 | } |
---|
400 | |
---|
401 | |
---|
402 | |
---|
403 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
404 | //THE FUNCTION THAT DO ALL THE SCHEMA USING THE PREVIOUS FUNCTION !!!!!!!!!!111 |
---|
405 | //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
406 | |
---|
407 | |
---|
408 | |
---|
409 | |
---|
410 | |
---|
411 | void SimCrys::general(SimCrys& sim, const int& pas, string outputpath, double C_rotation, double C_aperture, double C_offset, double C_tilt, double Crystal_tilt) |
---|
412 | { |
---|
413 | srand ( time(NULL) ); |
---|
414 | double AV_xp0(0); |
---|
415 | double AV_yp0(0); |
---|
416 | double RMSxp0(0); |
---|
417 | double RMSyp0(0); |
---|
418 | double AV_xp1(0); |
---|
419 | double AV_yp1(0); |
---|
420 | double long RMSxp1(0); |
---|
421 | double long RMSyp1(0); |
---|
422 | double av_pc1(0); |
---|
423 | double av_pc0(0); |
---|
424 | double SIGxp0(0); |
---|
425 | double SIGyp0(0); |
---|
426 | double SIGxp1(0); |
---|
427 | double SIGyp1(0); |
---|
428 | double xp0(0); |
---|
429 | double xp1(0); |
---|
430 | double count(0); |
---|
431 | double zpos(-0.1); |
---|
432 | double xfin(0); |
---|
433 | double yfin(0); |
---|
434 | int const col(15); |
---|
435 | |
---|
436 | |
---|
437 | //-----------------------------------------------NOW THE PARAMETERS FOR THE CHANGE OF REFERENTIAL ------------------------------------------------ |
---|
438 | int mat; |
---|
439 | |
---|
440 | double p0; |
---|
441 | double zlm; |
---|
442 | double shift; |
---|
443 | double x_shift, xp_shift, s_shift; //coordinates after shift/rotation |
---|
444 | double x_rot, xp_rot, s_rot; |
---|
445 | double x_temp, xp_temp, s_temp; //all the _temp variables are used when you hit the cry from below |
---|
446 | double tilt_int, x_int, xp_int, s_int; //all the _int variables are used when you hit the cry from below (int=interaction point) |
---|
447 | double x00; |
---|
448 | double z; |
---|
449 | double z00; |
---|
450 | double zp; |
---|
451 | double dpop; |
---|
452 | double s; |
---|
453 | double a_eq, b_eq, c_eq, Delta; |
---|
454 | int part_abs; |
---|
455 | double x; |
---|
456 | double xp; |
---|
457 | double y; |
---|
458 | double yp; |
---|
459 | double p; //be careful: [Gev] |
---|
460 | double s_in; |
---|
461 | |
---|
462 | // adding variables for the pencil beam. Variables in the absolute reference frame. |
---|
463 | |
---|
464 | double x_in0; |
---|
465 | double xp_in0; |
---|
466 | double y_in0; |
---|
467 | double yp_in0; |
---|
468 | double p_in0; //be careful: [Gev] |
---|
469 | double s_in0; |
---|
470 | double s_impact; |
---|
471 | double totcount(0); |
---|
472 | |
---|
473 | double mirror; |
---|
474 | double tiltangle; |
---|
475 | |
---|
476 | |
---|
477 | double c_length; //Length in m |
---|
478 | double c_rotation(C_rotation) ; //Rotation angle vs vertical in radian initiated and fixed at 0. Can BE CHANGED IF NEED |
---|
479 | double c_aperture(C_aperture) ; //Aperture in m |
---|
480 | double c_offset(C_offset) ; //Offset in m |
---|
481 | double c_tilt[2] = {C_tilt, 0} ; //Tilt in radian |
---|
482 | double c_tilt0[2] ; //Tilt in radian |
---|
483 | |
---|
484 | |
---|
485 | double xp_tangent; |
---|
486 | |
---|
487 | double Cry_tilt(Crystal_tilt); //crystal tilt [rad] ///////////////////////////////////++++++++++++++++++++++++(IF NEEDED, put it into the input file !!! See later)+++++++++++++++//////////////// |
---|
488 | |
---|
489 | |
---|
490 | |
---|
491 | //----------------------------------END OF INITATING REFERNENTIAL PARAMETERS------------------------------------------------------------------------- |
---|
492 | c_length = crys.Cry_length; |
---|
493 | |
---|
494 | if (crys.IS == 0) { |
---|
495 | mat = 8; |
---|
496 | } else if (crys.IS == 1) { |
---|
497 | mat = 9; |
---|
498 | } else if (crys.IS == 2) { |
---|
499 | mat = 10; |
---|
500 | } else if (crys.IS == 3) { |
---|
501 | mat = 11; |
---|
502 | } else { |
---|
503 | cout << "Material of the Crystal not fund !!!" << endl; |
---|
504 | } |
---|
505 | |
---|
506 | |
---|
507 | if (c_length > 0 ) { |
---|
508 | p0 = part.PC; |
---|
509 | c_tilt0[0] = c_tilt[0]; |
---|
510 | c_tilt0[1] = c_tilt[1]; |
---|
511 | tiltangle = c_tilt0[0]; |
---|
512 | } |
---|
513 | |
---|
514 | |
---|
515 | |
---|
516 | //ofstream sortie("results.csv"); |
---|
517 | //sortie << "xp1" <<","<< "xp0"<<","<< "x"<<","<< "y"<<endl; |
---|
518 | ofstream sortieIco; |
---|
519 | string name1; |
---|
520 | name1 = outputpath + "/ascii_after.csv"; |
---|
521 | |
---|
522 | sortieIco.open(name1.c_str()); |
---|
523 | |
---|
524 | string rest; |
---|
525 | |
---|
526 | ifstream entree; |
---|
527 | string name2; |
---|
528 | name2 = outputpath + "/ascii_before.csv"; |
---|
529 | |
---|
530 | entree.open(name2.c_str()); |
---|
531 | if (entree) { |
---|
532 | while (!entree.eof()) { |
---|
533 | count = count + 1; |
---|
534 | |
---|
535 | getline(entree, rest, ','); |
---|
536 | part.x = atof(rest.c_str()); |
---|
537 | |
---|
538 | |
---|
539 | //cout<<" X " << part.x<<endl; |
---|
540 | |
---|
541 | getline(entree, rest, ','); |
---|
542 | part.y = atof(rest.c_str()); |
---|
543 | |
---|
544 | |
---|
545 | //cout<<" Y " << part.y<<endl; |
---|
546 | |
---|
547 | getline(entree, rest, ','); |
---|
548 | part.xp = atof(rest.c_str()); |
---|
549 | |
---|
550 | //cout<<" XP " << part.xp<<endl; |
---|
551 | |
---|
552 | |
---|
553 | |
---|
554 | getline(entree, rest, ','); |
---|
555 | part.yp = atof(rest.c_str()); |
---|
556 | |
---|
557 | |
---|
558 | //cout<<" YP " << part.yp<<endl; |
---|
559 | |
---|
560 | getline(entree, rest); |
---|
561 | part.PC = atof(rest.c_str()); |
---|
562 | |
---|
563 | if (part.PC != 0) { // POUR OTER LA DERNIERE PART QUI VAUT 0 POUR TOUT |
---|
564 | |
---|
565 | av_pc0 = av_pc0 + part.PC; |
---|
566 | |
---|
567 | //part.x=random_gauss()*0.00115; |
---|
568 | //part.xp=random_gauss()*sqrt(5*1e-6); |
---|
569 | //cout<<"X et Y avt :"<<part.x<<" , "<<part.y<<endl; |
---|
570 | |
---|
571 | // IL FAUT PASSER DES CM EN M!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11 |
---|
572 | //bdelta=tan(part.xp)*10; |
---|
573 | part.x = (part.x); //-0.0045; |
---|
574 | |
---|
575 | |
---|
576 | //bdelta=tan(part.yp)*10; |
---|
577 | part.y = (part.y); //-0.0045; |
---|
578 | |
---|
579 | |
---|
580 | |
---|
581 | //cout<<"X et Y apres :"<<part.x<<" , "<<part.y<<endl; |
---|
582 | xp0 = part.xp; |
---|
583 | |
---|
584 | //part.y=random_gauss()*sqrt(0.00095); |
---|
585 | //part.yp=5*1e-6; |
---|
586 | |
---|
587 | |
---|
588 | //part.PC=120.976; |
---|
589 | //cout<<" X :"<<part.x<<" XP :"<<part.xp<<" Y :"<<part.y<<" YP :"<<part.yp<<" PC :"<<part.PC<<endl; |
---|
590 | AV_xp0 = AV_xp0 + part.xp; //average x angle before impact |
---|
591 | AV_yp0 = AV_yp0 + part.yp; //average y angle before impact |
---|
592 | RMSxp0 = RMSxp0 + part.xp * part.xp; //rms x angle before impact |
---|
593 | RMSyp0 = RMSyp0 + part.yp * part.yp; //rms y angle before impact |
---|
594 | |
---|
595 | |
---|
596 | // NOW WE NEED TO CHANGE THE REFERENTIAL!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
597 | |
---|
598 | |
---|
599 | |
---|
600 | s = 0; |
---|
601 | x = part.x; |
---|
602 | xp = part.xp; |
---|
603 | z = part.y; |
---|
604 | zp = part.yp; |
---|
605 | p = part.PC; |
---|
606 | |
---|
607 | x_temp = 0; |
---|
608 | x_int = 0; |
---|
609 | x_rot = 0; |
---|
610 | x_shift = 0; |
---|
611 | s_temp = 0; |
---|
612 | s_int = 0; |
---|
613 | s_rot = 0; |
---|
614 | s_shift = 0; |
---|
615 | xp_int = 0; |
---|
616 | xp_temp = 0; |
---|
617 | xp_rot = 0; |
---|
618 | xp_shift = 0; |
---|
619 | shift = 0; |
---|
620 | tilt_int = 0; |
---|
621 | dpop = (p - p0) / p0; |
---|
622 | |
---|
623 | // ------------------------------transform particle coordinates to get into collimator coordinate system ----------------------------------------- |
---|
624 | |
---|
625 | |
---|
626 | // first check whether particle was lost before |
---|
627 | totcount = totcount + 1; |
---|
628 | if ((x < 99.0 * 1e-3) && (z < 99.0 * 1e-3)) { ////////////////////////// IF NUMERO 1////////////////// |
---|
629 | |
---|
630 | // first do rotation into collimator frame |
---|
631 | |
---|
632 | x = part.x * cos(c_rotation) + sin(c_rotation) * part.y; |
---|
633 | z = part.y * cos(c_rotation) - sin(c_rotation) * part.x; |
---|
634 | xp = part.xp * cos(c_rotation) + sin(c_rotation) * part.yp; |
---|
635 | zp = part.yp * cos(c_rotation) - sin(c_rotation) * part.xp; |
---|
636 | |
---|
637 | //cout<<"PARTICULE NON PERDU>>><<"<<endl; |
---|
638 | |
---|
639 | mirror = 1; |
---|
640 | x = mirror * x; |
---|
641 | xp = mirror * xp; |
---|
642 | |
---|
643 | //Shift with opening and offset |
---|
644 | |
---|
645 | x = x - c_aperture / 2 - c_offset; |
---|
646 | |
---|
647 | // Include collimator tilt |
---|
648 | |
---|
649 | if (tiltangle > 0) { ////////////////////////// IF NUMERO 2////////////////// |
---|
650 | xp = xp - tiltangle; |
---|
651 | } //////////////////////////FIN IF NUMERO 2////////////////// |
---|
652 | else if (tiltangle < 0) { ////////////////////////// IF NUMERO 3////////////////// |
---|
653 | x = x + sin(tiltangle) * c_length; |
---|
654 | xp = xp - tiltangle; |
---|
655 | } //////////////////////////FIN IF NUMERO 3////////////////// |
---|
656 | |
---|
657 | |
---|
658 | |
---|
659 | if (Cry_tilt < 0) { ////////////////////////// IF NUMERO 4////////////////// |
---|
660 | |
---|
661 | s_shift = s; |
---|
662 | shift = crys.Rcurv * (1 - cos(Cry_tilt)); |
---|
663 | if (Cry_tilt < (-crys.cry_bend) ) { ////////////////////////// IF NUMERO 5////////////////// |
---|
664 | shift = ( crys.Rcurv * ( cos ( - Cry_tilt) - cos( crys.cry_bend - Cry_tilt) ) ); |
---|
665 | } //////////////////////////FIN IF NUMERO 5////////////////// |
---|
666 | x_shift = x - shift; |
---|
667 | } //////////////////////////fin IF NUMERO 4////////////////// |
---|
668 | else { //////////////////////////IF NUMERO 6////////////////// |
---|
669 | //cout<< "CRY-TILT POSSITIVE ++ or nul"<<endl; |
---|
670 | s_shift = s; |
---|
671 | x_shift = x; |
---|
672 | } //////////////////////////fin IF NUMERO 6////////////////// |
---|
673 | // cout<<"debug - S shift" << s_shift<<endl; |
---|
674 | //cout<<"debug - X shift" << x_shift <<endl; |
---|
675 | // |
---|
676 | // 2nd transformation: rotation |
---|
677 | s_rot = x_shift * sin(Cry_tilt) + s_shift * cos(Cry_tilt); |
---|
678 | x_rot = x_shift * cos(Cry_tilt) - s_shift * sin(Cry_tilt); |
---|
679 | xp_rot = xp - Cry_tilt; |
---|
680 | //cout<< "debug - S rot" << s_rot <<endl; |
---|
681 | //cout<<"debug - X rot" << x_rot <<endl; |
---|
682 | //cout<<"debug - XP rot" << xp_rot<<endl; |
---|
683 | |
---|
684 | // 3rd transformation: drift to the new coordinate s=0 |
---|
685 | xp = xp_rot; |
---|
686 | x = x_rot - xp_rot * s_rot; |
---|
687 | z = z - zp * s_rot; |
---|
688 | s = 0; |
---|
689 | //cout<<"debug - S cryRF" << s_rot <<endl; |
---|
690 | //cout<<"debug - X cryRF" << x_rot <<endl; |
---|
691 | //cout<<"debug - XP cryRF" << xp_rot <<endl; |
---|
692 | |
---|
693 | |
---|
694 | //----------------------------------NOW CHECK IF THE PARTICLE HIT the crystal--------------------------------------------------------------- |
---|
695 | //cout<<" X " << x<<endl; |
---|
696 | if (x >= 0) { ////////////////////////// IF NUMERO 7////////////////// |
---|
697 | s_impact = s_in0; //!(for the first impact) |
---|
698 | //cout<<"hit the cry entrance face"<<endl; |
---|
699 | //cout<<"impact at s,x = "<< s_impact<<x_in0<<endl; |
---|
700 | //cout<<"with angle xp = "<<xp<<endl; |
---|
701 | //cout<<"s before"<< s<<endl; |
---|
702 | |
---|
703 | part.x = x; |
---|
704 | part.xp = xp; |
---|
705 | part.y = z; |
---|
706 | part.yp = zp; |
---|
707 | part.PC = p; |
---|
708 | //cout<< "ICI AHHHHHHHHHHHHHHHHHHHHHHHHHHHH"<<endl; |
---|
709 | //cout<<" xp rel and xp crit " << part.xp_rel<< " " << crys.xpcrit<<endl; |
---|
710 | bases(sim); |
---|
711 | //cout << part.xp<< " XP "<<endl; |
---|
712 | bool crossing(cross(sim)); |
---|
713 | //cout<<part.xp-xp0<<endl; |
---|
714 | //cout<<" xp rel and xp crit " << part.xp_rel<< " " << crys.xpcrit<<endl; |
---|
715 | bool layering; |
---|
716 | if (crossing != 0) { ////////////////////////// IF NUMERO 8////////////////// |
---|
717 | layering = layer(sim); |
---|
718 | //cout<<"Crossing !!!!"<<endl; |
---|
719 | |
---|
720 | if (layering == true) { ////////////////////////// IF NUMERO 9////////////////// |
---|
721 | parameters(sim); |
---|
722 | //cout<<" Layering =true " <<endl; |
---|
723 | |
---|
724 | if ((abs(part.xp_rel) < crys.xpcrit)) { ////////////////////////// IF NUMERO 10////////////////// |
---|
725 | channel(sim); |
---|
726 | //cout<< " CHANNELING " << endl; |
---|
727 | } ////////////////////////// fin IF NUMERO 10////////////////// |
---|
728 | else { ////////////////////////// IF NUMERO 11////////////////// |
---|
729 | reflection(sim); |
---|
730 | //cout<<" REFLECTION "<<endl; |
---|
731 | } ////////////////////////// fin IF NUMERO 11////////////////// |
---|
732 | |
---|
733 | |
---|
734 | } ////////////////////////// fin IF NUMERO 9////////////////// |
---|
735 | } ////////////////////////// fin IF NUMERO 8////////////////// |
---|
736 | |
---|
737 | xp = part.xp; |
---|
738 | |
---|
739 | //} |
---|
740 | |
---|
741 | s = crys.Rcurv * sin(crys.cry_bend); |
---|
742 | zlm = crys.Rcurv * sin(crys.cry_bend); |
---|
743 | //cout<<"process:"<<PROC<<endl; |
---|
744 | //cout<<"s after"<< s<<endl; |
---|
745 | //cout<<"part.xp"<<part.xp<<endl; |
---|
746 | |
---|
747 | //cout<<"----------------------1111111111111111111111111-----------------"<<endl; |
---|
748 | |
---|
749 | |
---|
750 | } ////////////////////////// fin IF NUMERO 7////////////////// |
---|
751 | else { ////////////////////////// IF NUMERO 12////////////////// |
---|
752 | |
---|
753 | //cout<<"OU LA ABBBBBBBBBBBB" <<endl; |
---|
754 | xp_tangent = sqrt((-2 * x * crys.Rcurv + x * x) / (crys.Rcurv * crys.Rcurv)); |
---|
755 | //cout<<"RCURV "<<crys.Rcurv<<endl; |
---|
756 | //cout<<"tangent"<<xp_tangent<<"angle"<<xp<<endl; |
---|
757 | |
---|
758 | //cout<<"s tan " << crys.Rcurv*sin(xp_tangent)<<endl; |
---|
759 | //cout<< " s tot "<< c_length<< crys.Rcurv*sin(crys.cry_bend)<<endl; |
---|
760 | if ( xp >= xp_tangent) { ////////////////////////// IF NUMERO 13////////////////// |
---|
761 | |
---|
762 | //if it hits the crystal, calculate in which point and apply the |
---|
763 | //transformation and drift to that point |
---|
764 | a_eq = (1 + xp * xp); |
---|
765 | b_eq = 2 * xp * (x - crys.Rcurv); |
---|
766 | c_eq = -2 * x * crys.Rcurv + x * x; |
---|
767 | Delta = b_eq * b_eq - 4 * (a_eq * c_eq); |
---|
768 | s_int = (-b_eq - sqrt(Delta)) / (2 * a_eq); |
---|
769 | //cout<< " s int "<<s_int <<endl; |
---|
770 | if (s_int < crys.Rcurv * sin(crys.cry_bend)) { ////////////////////////// fin IF NUMERO 14////////////////// |
---|
771 | // transform to a new ref system:shift and rotate |
---|
772 | |
---|
773 | x_int = xp * s_int + x; |
---|
774 | xp_int = xp; |
---|
775 | z = z + zp * s_int; |
---|
776 | x = 0; |
---|
777 | s = 0; |
---|
778 | |
---|
779 | |
---|
780 | // tilt_int=2*X_int/S_int |
---|
781 | tilt_int = s_int / crys.Rcurv; |
---|
782 | xp = xp - tilt_int; |
---|
783 | |
---|
784 | |
---|
785 | |
---|
786 | part.x = x; |
---|
787 | part.xp = xp; |
---|
788 | part.y = z; |
---|
789 | part.yp = zp; |
---|
790 | part.PC = p; |
---|
791 | crys.Cry_length = crys.Cry_length - (tilt_int * crys.Rcurv); |
---|
792 | bases(sim); |
---|
793 | |
---|
794 | //cout<<"part.xp"<<part.xp<<endl; |
---|
795 | |
---|
796 | bool crossing(cross(sim)); |
---|
797 | //cout<<"Crossing !!!!"<<endl; |
---|
798 | bool layering(layer(sim)); |
---|
799 | |
---|
800 | |
---|
801 | |
---|
802 | if (layering == true) { |
---|
803 | //cout<<" Layering =true " <<endl; |
---|
804 | parameters(sim); |
---|
805 | if ((abs(part.xp_rel) < crys.xpcrit)) { |
---|
806 | channel(sim); |
---|
807 | //cout<< " CHANNELING " << endl; |
---|
808 | } else { |
---|
809 | reflection(sim); |
---|
810 | //cout<< " REFLECTION " << endl; |
---|
811 | } |
---|
812 | |
---|
813 | |
---|
814 | } |
---|
815 | |
---|
816 | //cout<<"part.xp"<<part.xp<<endl; |
---|
817 | //cout<<"----------------------2222222222222222222-----------------"<<endl; |
---|
818 | |
---|
819 | s = crys.Rcurv * sin(crys.cry_bend - tilt_int); |
---|
820 | zlm = crys.Rcurv * sin(crys.cry_bend - tilt_int); |
---|
821 | |
---|
822 | |
---|
823 | if (crossing != 0) { ////////////////////////// IF NUMERO 15////////////////// |
---|
824 | x_rot = x_int; |
---|
825 | s_rot = s_int; |
---|
826 | xp_rot = xp_int; |
---|
827 | s_shift = s_rot * cos(-Cry_tilt) + x_rot * sin(-Cry_tilt); |
---|
828 | x_shift = -s_rot * sin(-Cry_tilt) + x_rot * cos(-Cry_tilt); |
---|
829 | xp_shift = xp_rot + Cry_tilt; |
---|
830 | if (Cry_tilt < 0) { ////////////////////////// IF NUMERO 16////////////////// |
---|
831 | s_impact = s_shift; |
---|
832 | x_in0 = x_shift + shift; |
---|
833 | xp_in0 = xp_shift; |
---|
834 | } ////////////////////////// fin IF NUMERO 16////////////////// |
---|
835 | else { ////////////////////////// IF NUMERO 17////////////////// |
---|
836 | x_in0 = x_shift; |
---|
837 | s_impact = s_shift; |
---|
838 | xp_in0 = xp_shift; |
---|
839 | } ////////////////////////// fin IF NUMERO 17////////////////// |
---|
840 | |
---|
841 | } ////////////////////////// fin IF NUMERO 15////////////////// |
---|
842 | |
---|
843 | |
---|
844 | x_temp = x; |
---|
845 | s_temp = s; |
---|
846 | xp_temp = xp; |
---|
847 | s = s_temp * cos(-tilt_int) + x_temp * sin(-tilt_int); |
---|
848 | x = -s_temp * sin(-tilt_int) + x_temp * cos(-tilt_int); |
---|
849 | xp = xp_temp + tilt_int; |
---|
850 | |
---|
851 | // 2nd: shift back the 2 axis |
---|
852 | x = x + x_int; |
---|
853 | s = s + s_int; |
---|
854 | |
---|
855 | } ////////////////////////// fin IF NUMERO 14////////////////// |
---|
856 | else { ///////////////////////////////////////////////PROBLEM : TOO MUCH OF PARTICULES WITHOUT ANY CHANGE OF ANGLE///////////THE PROB IS HERE WITH THOSES TWO "ELSE"///////////////////////// |
---|
857 | |
---|
858 | s = crys.Rcurv * sin(crys.Cry_length / crys.Rcurv); |
---|
859 | x = x + s * xp; |
---|
860 | z = z + s * zp; |
---|
861 | ++part.nout; |
---|
862 | //cout<< "the particule does not hit the crystal (Trop basse dans neg.)"<<endl; |
---|
863 | } |
---|
864 | } else { |
---|
865 | s = crys.Rcurv * sin(crys.Cry_length / crys.Rcurv); |
---|
866 | x = x + s * xp; |
---|
867 | z = z + s * zp; |
---|
868 | ++part.nout; |
---|
869 | //cout<< "the particule does not hit the crystal (negatif et plus petit angle qu xp_tang)"<<endl; |
---|
870 | } |
---|
871 | //} ////////////////////////// fin IF NUMERO 13////////////////// |
---|
872 | |
---|
873 | }////////////////////// 12 |
---|
874 | |
---|
875 | |
---|
876 | //trasform back from the crystal to the collimator reference system |
---|
877 | // 1st: un-rotate the coordinates |
---|
878 | |
---|
879 | |
---|
880 | x_rot = x; |
---|
881 | s_rot = s; |
---|
882 | xp_rot = xp; |
---|
883 | //cout<< "debug - S Cry RF 2" , s_rot <<endl; |
---|
884 | //cout<<"debug - X Cry RF 2" , x_rot <<endl; |
---|
885 | //cout<<"debug - XP Cry RF 2" , xp_rot <<endl; |
---|
886 | s_shift = s_rot * cos(-Cry_tilt) + x_rot * sin(-Cry_tilt); |
---|
887 | x_shift = -s_rot * sin(-Cry_tilt) + x_rot * cos(-Cry_tilt); |
---|
888 | xp_shift = xp_rot + Cry_tilt; |
---|
889 | // 2nd: shift back the reference frame |
---|
890 | if (Cry_tilt < 0) { ////////////////////////// IF NUMERO 18////////////////// |
---|
891 | s = s_shift; |
---|
892 | x = x_shift + shift; |
---|
893 | xp = xp_shift; |
---|
894 | } ////////////////////////// fin IF NUMERO 18////////////////// |
---|
895 | else { |
---|
896 | x = x_shift; |
---|
897 | s = s_shift; |
---|
898 | xp = xp_shift; |
---|
899 | } |
---|
900 | // 3rd: shift to new S=Length position |
---|
901 | x = xp * (c_length - s) + x; |
---|
902 | z = zp * (c_length - s) + z; |
---|
903 | s = c_length; |
---|
904 | |
---|
905 | |
---|
906 | //cout<< "debug - S Coll RF 2" , s_rot <<endl; |
---|
907 | //cout<<"debug - X Coll RF 2" , x_rot <<endl; |
---|
908 | //cout<<"debug - XP Coll RF 2" , xp_rot <<endl; |
---|
909 | |
---|
910 | //cout<<"X1_coll"<<x<<"Z1_coll"<<z<<"XP1_coll"<<xp<<"ZP1_coll"<<zp <<"s" <<s<<endl; |
---|
911 | |
---|
912 | if (part.effet == 2) { |
---|
913 | part_abs = 0; |
---|
914 | } else { |
---|
915 | part_abs = 1; |
---|
916 | } |
---|
917 | |
---|
918 | |
---|
919 | //========================================================================================================================= |
---|
920 | // Transform back to particle coordinates with opening and offset |
---|
921 | |
---|
922 | if (part_abs == 0) { ////////////////////////// IF NUMERO 19////////////////// |
---|
923 | // |
---|
924 | // Include collimator tilt |
---|
925 | |
---|
926 | if (tiltangle > 0) { |
---|
927 | x = x + tiltangle * c_length; |
---|
928 | xp = xp + tiltangle; |
---|
929 | } else if (tiltangle < 0) { |
---|
930 | x = x + tiltangle * c_length; |
---|
931 | xp = xp + tiltangle; |
---|
932 | |
---|
933 | x = x - sin(tiltangle) * c_length; |
---|
934 | } |
---|
935 | |
---|
936 | // Transform back to particle coordinates with opening and offset |
---|
937 | |
---|
938 | z00 = z; |
---|
939 | x00 = x + mirror * c_offset; |
---|
940 | x = x + c_aperture / 2 + mirror * c_offset; |
---|
941 | |
---|
942 | //+ Now mirror at the horizontal axis for negative X offset |
---|
943 | |
---|
944 | x = mirror * x; |
---|
945 | xp = mirror * xp; |
---|
946 | |
---|
947 | // Last do rotation into collimator frame |
---|
948 | |
---|
949 | part.x = x * cos((-1) * c_rotation) + z * sin((-1) * c_rotation); |
---|
950 | part.y = z * cos((-1) * c_rotation) - x * sin((-1) * c_rotation); |
---|
951 | part.xp = xp * cos((-1) * c_rotation) + zp * sin((-1) * c_rotation); |
---|
952 | part.yp = zp * cos((-1) * c_rotation) - xp * sin((-1) * c_rotation); |
---|
953 | |
---|
954 | |
---|
955 | //p_in = (1 + dpop) * p0; |
---|
956 | s_in = s_in + s; |
---|
957 | |
---|
958 | if (part.effet == 1) { ////////////////////////// IF NUMERO 21////////////////// |
---|
959 | |
---|
960 | part.x = 99.99e-3; |
---|
961 | part.y = 99.99e-3; |
---|
962 | //cout<< "Mean that the particules is lost"<<endl; |
---|
963 | } ////////////////////////// fin IF NUMERO 21////////////////// |
---|
964 | } ////////////////////////// fin IF NUMERO 19////////////////// |
---|
965 | |
---|
966 | |
---|
967 | // End of check for particles not being lost before (see @330) |
---|
968 | |
---|
969 | //} ////////////////////////// fin IF NUMERO 12////////////////// |
---|
970 | // End of loop over all particles |
---|
971 | |
---|
972 | |
---|
973 | |
---|
974 | } //collimator with length = 0 ////////////////////////// fin IF NUMERO 1////////////////// |
---|
975 | |
---|
976 | // =================================================================================FIN DE RETOUR AU COORDONEE DE ICOSIM==================================================================== |
---|
977 | /* |
---|
978 | part.xp=xp; |
---|
979 | part.yp=zp; |
---|
980 | part.y=z; |
---|
981 | part.x=x; |
---|
982 | part.PC=p; |
---|
983 | */ |
---|
984 | xp1 = part.xp; |
---|
985 | |
---|
986 | |
---|
987 | //if(xp1-xp0!=0){ |
---|
988 | //sortie << xp1 << "," << xp0<<"," << xfin<<"," << yfin<<endl; |
---|
989 | //} |
---|
990 | //cout<<"ICIIIIIIIIIIIIIIIIIIIIIIIII"<<part.PC<<endl; |
---|
991 | |
---|
992 | |
---|
993 | AV_xp1 = AV_xp1 + part.xp; //average x angle after interaction |
---|
994 | AV_yp1 = AV_yp1 + part.yp; //average y angle after interaction |
---|
995 | RMSxp1 = RMSxp1 + part.xp * part.xp; //rms x angle after interaction |
---|
996 | RMSyp1 = RMSyp1 + part.xp * part.xp; //rms y angle after interaction |
---|
997 | av_pc1 = av_pc1 + part.PC; |
---|
998 | |
---|
999 | |
---|
1000 | //cout<<" X apres :"<<part.x<<" XP apres :"<<part.xp<<" Y apres :"<<part.y<<" YP apres :"<<part.yp<<" PC apres :"<<part.PC<<endl; |
---|
1001 | //cout<<endl; |
---|
1002 | //sortieIco.setf(ios::scientific); |
---|
1003 | sortieIco << part.x << "," << part.y << "," << part.xp << "," << part.yp << "," << part.PC << endl; |
---|
1004 | //}//end of for.... |
---|
1005 | |
---|
1006 | } |
---|
1007 | } //End of while |
---|
1008 | |
---|
1009 | |
---|
1010 | } else { |
---|
1011 | //cout<< "We can not open the file !" << endl; |
---|
1012 | } |
---|
1013 | entree.close(); |
---|
1014 | |
---|
1015 | |
---|
1016 | |
---|
1017 | AV_xp0 = AV_xp0 / count; |
---|
1018 | AV_yp0 = AV_yp0 / count; |
---|
1019 | |
---|
1020 | RMSxp0 = sqrt(RMSxp0 / (count)); |
---|
1021 | RMSyp0 = sqrt(RMSyp0 / (count)); |
---|
1022 | |
---|
1023 | SIGxp0 = sqrt((RMSxp0 * RMSxp0) - (AV_xp0 * AV_xp0)); |
---|
1024 | SIGyp0 = sqrt((RMSyp0 * RMSyp0) - (AV_yp0 * AV_yp0)); |
---|
1025 | |
---|
1026 | AV_xp1 = AV_xp1 / count; |
---|
1027 | AV_yp1 = AV_yp1 / count; |
---|
1028 | RMSxp1 = sqrt(RMSxp1 / (count)); |
---|
1029 | RMSyp1 = sqrt(RMSyp1 / (count)); |
---|
1030 | |
---|
1031 | SIGxp1 = sqrt((RMSxp1 * RMSxp1) - (AV_xp1 * AV_xp1)); |
---|
1032 | SIGyp1 = sqrt((RMSyp1 * RMSyp1) - (AV_yp1 * AV_yp1)); |
---|
1033 | av_pc0 = av_pc0 / count; |
---|
1034 | av_pc1 = av_pc1 / count; |
---|
1035 | //-------------------------------------------------------------------------- |
---|
1036 | //-------------write a summary file----------------------------------------- |
---|
1037 | //cout<<endl; |
---|
1038 | //cout<<endl; |
---|
1039 | //cout<<endl; |
---|
1040 | |
---|
1041 | //cout<<"avg xp_in: "<< AV_xp0 <<"// avg xp_out: "<< AV_xp1<<endl; |
---|
1042 | //cout<<"rms xp_in: "<< RMSxp0 <<"// rms xp_out: "<< RMSxp1<<endl; |
---|
1043 | //cout<<"sigma xp_in: "<< SIGxp0 <<"// sigma xp_out: "<< SIGxp1<<endl; |
---|
1044 | //cout<<"avg pc_in: "<< av_pc0 <<"// avg pc_out: "<< av_pc1<<endl; |
---|
1045 | |
---|
1046 | |
---|
1047 | |
---|
1048 | sortieIco.close(); |
---|
1049 | //cout<<endl; |
---|
1050 | //affiche(); |
---|
1051 | file_out(pas, outputpath); |
---|
1052 | |
---|
1053 | //cout<<endl; |
---|
1054 | //cout<<endl; |
---|
1055 | //cout<<endl; |
---|
1056 | } |
---|
1057 | |
---|
1058 | |
---|
1059 | |
---|
1060 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
1061 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
1062 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
1063 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
1064 | |
---|
1065 | //WE WRITE SOME FUNCTIONS IMPORTANT FUNCTION NOW |
---|
1066 | |
---|
1067 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
1068 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
1069 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
1070 | /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// |
---|
1071 | |
---|
1072 | |
---|
1073 | void SimCrys::move_am(int IS, int NAM, double DZ, double DEI, double DLY, double DLR, double& xp, double& yp, double& PC) |
---|
1074 | { |
---|
1075 | //cout<<"The move_am function is use !!!!! "<<endl; |
---|
1076 | double DYA(0); |
---|
1077 | double Wp(0); |
---|
1078 | if (NAM == 0) { |
---|
1079 | return; |
---|
1080 | } |
---|
1081 | xp = xp * 1000; |
---|
1082 | yp = yp * 1000; |
---|
1083 | //cout << "xp avt: "<<xp<<" yp avt: "<<yp<<endl; |
---|
1084 | //DEI ---> dE/dx (stoping energy) |
---|
1085 | PC = PC - DEI * DZ; //energy lost beacause of ionization process[GeV] |
---|
1086 | |
---|
1087 | // Coulomb Scattering |
---|
1088 | //DYA ---> rms of coloumb scattering |
---|
1089 | DYA = (14 / PC) * sqrt(DZ / DLR); //rms scattering (mrad) |
---|
1090 | xp = xp + DYA * random_gauss_cut(); |
---|
1091 | yp = yp + DYA * random_gauss_cut(); |
---|
1092 | |
---|
1093 | //Elastic scattering |
---|
1094 | //cout<<"Plusieurs possibilitees : "<<endl; |
---|
1095 | if (random() <= (DZ / crys.DLAI[IS])) { |
---|
1096 | //cout<<"Poss. 1 !!!! "<<endl; |
---|
1097 | xp = xp + 196 * random_gauss_cut() / PC; //angle elast. scattering in R plane |
---|
1098 | yp = yp + 196 * random_gauss_cut() / PC; |
---|
1099 | } |
---|
1100 | |
---|
1101 | //Diffraction interaction |
---|
1102 | if (random() <= (DZ / (DLY * 6.143))) { |
---|
1103 | //cout<<"Poss. 2 !!!! "<<endl; |
---|
1104 | xp = xp + 257 * random_gauss_cut() / PC; // angle elast. scattering in R plane[mr] |
---|
1105 | yp = yp + 257 * random_gauss_cut() / PC; |
---|
1106 | Wp = random(); |
---|
1107 | PC = PC - 0.5 * pow(0.3 * PC, Wp); // m0*c = 1 GeV/c for particle |
---|
1108 | |
---|
1109 | } |
---|
1110 | |
---|
1111 | //Inelastic interaction |
---|
1112 | if (random() <= DZ / DLY) { |
---|
1113 | //cout<<" Absorbed !! "<< endl; |
---|
1114 | } |
---|
1115 | //cout<<"XP APRES : "<<xp<<" YP APRES : "<<yp<<endl; |
---|
1116 | xp = xp / 1000; |
---|
1117 | yp = yp / 1000; |
---|
1118 | } |
---|
1119 | |
---|
1120 | |
---|
1121 | double SimCrys::random() |
---|
1122 | { |
---|
1123 | double scale = RAND_MAX; |
---|
1124 | double base = rand() / scale; |
---|
1125 | double fine = rand() / scale; |
---|
1126 | return base + fine / scale; |
---|
1127 | } |
---|
1128 | |
---|
1129 | double SimCrys::random_dist() |
---|
1130 | { |
---|
1131 | double scale = RAND_MAX; |
---|
1132 | double base = rand() / scale; |
---|
1133 | double fine = rand() / scale; |
---|
1134 | return (-2.5 * 1e-4 + (3.5 * 1e-4 ) * (base + fine / scale)); |
---|
1135 | } |
---|
1136 | |
---|
1137 | double SimCrys::random_gauss() |
---|
1138 | { |
---|
1139 | const double PI (4.0 * atan(1.0)); |
---|
1140 | double U(random()); |
---|
1141 | double V(random()); |
---|
1142 | return (sqrt(-2 * log(U)) * cos(2 * PI * V)); |
---|
1143 | } |
---|
1144 | |
---|
1145 | double SimCrys::random_gauss_cut() |
---|
1146 | { |
---|
1147 | double resu; |
---|
1148 | do { |
---|
1149 | resu = random_gauss(); |
---|
1150 | } while (abs(resu) >= 1); |
---|
1151 | return resu; |
---|
1152 | } |
---|
1153 | |
---|
1154 | double SimCrys::random_gauss_dist() //POUR faire varier la moyenne de la gaussienne il faut mettre en argu la sim& !!!!! |
---|
1155 | { |
---|
1156 | const double PI (4.0 * atan(1.0)); |
---|
1157 | double U(random()); |
---|
1158 | double V(random()); |
---|
1159 | //moy=random(); |
---|
1160 | return 1e-4 * (sqrt(-2 * log(U)) * cos(2 * PI * V)); //+(0.0002*moy); //POUR faire varier la moyenne de la gaussienne aussi !!!!! |
---|
1161 | } |
---|
1162 | |
---|
1163 | /////////////////////////////////*******************************************************+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
1164 | //POUR 409 : |
---|
1165 | |
---|
1166 | double SimCrys::random_gauss_dist_409() //POUR faire varier la moyenne de la gaussienne il faut mettre en argu la sim& !!!!! |
---|
1167 | { |
---|
1168 | const double PI (4.0 * atan(1.0)); |
---|
1169 | double U(random()); |
---|
1170 | double V(random()); |
---|
1171 | //moy=random(); |
---|
1172 | return (8.939203e-06) * (sqrt(-2 * log(U)) * cos(2 * PI * V)) + (-7.183296e-08); //POUR faire varier la moyenne de la gaussienne aussi !!!!! |
---|
1173 | } |
---|
1174 | |
---|
1175 | ///////////////////////////////////************************************************************++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
1176 | |
---|
1177 | void SimCrys::affiche() |
---|
1178 | { |
---|
1179 | cout << "Crystal Curved Length [m] :" << crys.Cry_length << endl; |
---|
1180 | cout << "Crystal Curvature Radius [m] :" << crys.Rcurv << endl; |
---|
1181 | cout << "Crystal Bending Angle [urad] :" << crys.cry_bend << " , " << crys.Cry_length / crys.Rcurv << endl; |
---|
1182 | cout << "Critical Radius :" << crys.Rcrit << endl; |
---|
1183 | cout << "Ratio :" << crys.ratio << endl; |
---|
1184 | cout << "Critical Angle for straight Crystal :" << crys.xpcrit0 << endl; |
---|
1185 | cout << "Critical Angle for curved Crystal :" << crys.xpcrit << endl; |
---|
1186 | cout << "Angle average reflexion : " << part.Ang_avr << endl; |
---|
1187 | cout << "Full Channeling : " << crys.Cry_length / crys.Rcurv << endl; |
---|
1188 | cout << endl; |
---|
1189 | cout << endl; |
---|
1190 | cout << "N.AMORPH :" << part.namor << endl; |
---|
1191 | cout << "N.DECHANNELING:" << part.ndechann << endl; |
---|
1192 | cout << "N.CHANNELING :" << part.nchann << endl; |
---|
1193 | cout << "N.OUT :" << part.nout << endl; |
---|
1194 | cout << "N.VREFLECTION :" << part.nvrefl << endl; |
---|
1195 | cout << "N.VCAPTURE :" << part.nvcapt << endl; |
---|
1196 | |
---|
1197 | //cout << "Si la somme ne donne pas le nbre de part., c est normale !! En effet, il est possible que plusieurs evenement se produisent!!!"; |
---|
1198 | } |
---|
1199 | |
---|
1200 | void SimCrys::file_out(const int& pas, string outputpath) |
---|
1201 | { |
---|
1202 | |
---|
1203 | ofstream out; |
---|
1204 | string file; |
---|
1205 | file = outputpath + "/crystal_output.out"; |
---|
1206 | |
---|
1207 | out.open(file.c_str(), ios::out | ios::app); |
---|
1208 | |
---|
1209 | if (out.fail()) { |
---|
1210 | cerr << "Warning: problem writing in the file " << file << "!" << endl; |
---|
1211 | } else { |
---|
1212 | |
---|
1213 | out << endl; |
---|
1214 | out << "Crystal Curved Length [m] :" << crys.Cry_length << endl; |
---|
1215 | out << "Crystal Curvature Radius [m] :" << crys.Rcurv << endl; |
---|
1216 | out << "Crystal Bending Angle [urad] :" << crys.cry_bend << " , " << crys.Cry_length / crys.Rcurv << endl; |
---|
1217 | out << "Critical Radius :" << crys.Rcrit << endl; |
---|
1218 | out << "Ratio :" << crys.ratio << endl; |
---|
1219 | out << "Critical Angle for straight Crystal :" << crys.xpcrit0 << endl; |
---|
1220 | out << "Critical Angle for curved Crystal :" << crys.xpcrit << endl; |
---|
1221 | out << "Angle average reflexion : " << part.Ang_avr << endl; |
---|
1222 | out << "Full Channeling : " << crys.Cry_length / crys.Rcurv << endl; |
---|
1223 | out << endl; |
---|
1224 | out << endl; |
---|
1225 | out << "N.AMORPH :" << part.namor << endl; |
---|
1226 | out << "N.DECHANNELING:" << part.ndechann << endl; |
---|
1227 | out << "N.CHANNELING :" << part.nchann << endl; |
---|
1228 | out << "N.OUT :" << part.nout << endl; |
---|
1229 | out << "N.VREFLECTION :" << part.nvrefl << endl; |
---|
1230 | out << "N.VCAPTURE :" << part.nvcapt << endl << endl << endl; |
---|
1231 | |
---|
1232 | out.close(); |
---|
1233 | } |
---|
1234 | |
---|
1235 | } |
---|
1236 | |
---|