Context Navigation

source: ZHANGProjects/ICOSIM/CPP/trunk/source/simcrys.cc @ 2

Last change on this file since 2 was 2, checked in by zhangj, 10 years ago
Initial import
File size: 51.8 KB

Line
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.xpcrit0.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.00011pow(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(51e-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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: