Context Navigation

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

Last change on this file since 16 was 16, checked in by zhangj, 10 years ago
Fix the bug to generate random number move_am_() in the Fortran file.
File size: 63.1 KB

Line
1	#include <iostream>
2	#include <fstream>
3	#include <cmath>
4	#include <iomanip>
5	#include <string.h>
6	#include <stdio.h>
7	#include "simcrys.h"
8
9	//crystal routine of the proton beam
10	#include "crystalproton.h"
11	#include "crystalprotonfuns.h"
12
13	extern "C"
14	{
15	struct {
16	double Cry_length, Rcurv, C_xmax, C_ymax, Alayer, C_orient;
17	}Par_Cry1_;
18	};
19
20
21	using namespace std;
22
23	//Fortran crystal routine for the proton beam
24	// By Jianfeng Zhang @ LAL, 05/2014
25	// extern "C" {
26	// // void cryst_(int,double,double,double,double,double,double);
27	// void cryst_(int IS, double x, double xp, double y,double yp,double PC,double Length);
28	// }
29
30
31
32	SimCrys::SimCrys(Crystal crys, Partcrys part)
33	: crys(crys), part(part) //,moy(moy) //POUR faire varier la moyenne de la gaussienne aussi !!!!!
34	{
35	}
36
37	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
38	////////////////////////// FUNCTION BASES THAT MAKE SOME BASIC CALCULATION
39	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
40
41	void SimCrys::bases(SimCrys& sim)
42	{
43	if ( (crys.miscut < 0 ) && ( part.x > 0) && (part.x < -crys.Cry_length * tan(crys.miscut))) {
44	crys.L_chan = -part.x * tan(crys.miscut);
45	}
46	crys.tchan = crys.L_chan / crys.Rcurv;
47	part.xp_rel = part.xp - crys.miscut;
48	}
49
50
51
52	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
53	////////////////////////FUNCTIN CROSS TO KNOW IF THE PARTICULE CROSS THE CRYSTAL
54	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
55
56
57
58
59	bool SimCrys::cross(SimCrys& sim)
60	{
61	bool crossing(true);
62	if ((part.y < crys.ymin) \|\| (part.y > crys.ymax) \|\| (part.x > crys.xmax)) {
63	//cout << "OUT of the Crystal !!" << part.x << " , "<<part.xp << endl;
64	crossing = false;
65	part.x = part.x + part.xp * crys.s_length;
66
67	part.y = part.y + part.yp * crys.s_length;
68
69	part.nout = part.nout + 1; /////////////COUNT/////////////
70	part.effet = 1;
71	//cout<<" BLAAAAAAAAa"<<endl;
72
73	}
74
75	return crossing;
76	}
77
78	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
79	////////////////////////////FUNCTION LAYER THAT DETERMIN IF THE PARTICULE IS IN THE AMORPHOUS LAYER/////////////////////////////
80	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
81
82	bool SimCrys::layer(SimCrys& sim)
83	{
84	//cout<<"La fonction layer est utilisee>>>>>>>>>>>>>>"<<endl;
85	bool out(true);
86	if ((part.x < crys.Alayer) \|\| ((part.y - crys.ymin) < crys.Alayer) \|\| ((crys.ymax - part.y) < crys.Alayer)) {
87	cout<<"Amorphous layer! " << part.x <<" , " <<part.xp<<endl;
88	double a_eq;
89	double b_eq;
90	double c_eq;
91	double delta; //a, b, c, Delta of the second order eq.
92	double x0;
93	double y0;
94	double length_xs; //!Amorphous length
95	double length_ys; //!Amorphous length
96	double am_length; //!Amorphous length
97	out = false;
98	x0 = part.x;
99	y0 = part.y;
100	a_eq = (1 + part.xp * part.xp);
101	b_eq = (2 * (part.x) * (part.xp) - 2 * (part.xp) * crys.Rcurv);
102	c_eq = (part.x) * part.x - 2 * abs(part.x) * crys.Rcurv;
103	delta = b_eq * b_eq - 4 * a_eq * c_eq;
104	//cout<<"a : " << a_eq<<endl; //POUR TESTER FONCTION
105	//cout<<"b: " << b_eq <<endl; //POUR TESTER FONCTION
106	//cout<<"c : " << c_eq <<endl; //POUR TESTER FONCTION
107	//cout<<"Delta : " << delta <<endl; //POUR TESTER FONCTION
108	part.s = (-b_eq + sqrt(delta)) / ((2 * a_eq)) ; //in [m]
109	//cout<<"s : " << part.s <<endl; //POUR TESTER FONCTION
110	if (part.s >= crys.s_length) {
111	part.s = crys.s_length;
112	}
113	//cout<<"s : " << part.s <<endl; //POUR TESTER FONCTION
114	part.x = part.xp * part.s + x0 ; //in [m]
115	//cout<<"x0 : " << x0 <<endl; //POUR TESTER FONCTION
116	//cout<<"x : " << part.x <<endl; //POUR TESTER FONCTION
117
118	length_xs = sqrt((part.x - x0) * (part.x - x0) + part.s * part.s);
119	//cout<<"Length_XS : " << length_xs <<endl; //POUR TESTER FONCTION
120
121	if ((((part.yp >= 0) && ((part.y + part.yp * part.s) <= crys.ymax))) \|\| (((part.yp < 0) && ((part.y + part.yp * part.s) >= crys.ymin)))) {
122	length_ys = part.yp * length_xs;
123	} else {
124	part.s = (crys.ymax - part.y) / part.yp;
125	length_ys = sqrt((crys.ymax - part.y) * (crys.ymax - part.y) + part.s * part.s);
126	part.x = x0 + part.xp * part.s;
127	length_xs = sqrt((part.x - x0) * (part.x - x0) + part.s * part.s);
128	}
129
130	am_length = sqrt(length_xs * length_xs + length_ys * length_ys);
131	//cout << "AM_Length :" << am_length << endl; // TESTER FONCTION
132	part.s = part.s / 2;
133	part.x = x0 + part.xp * part.s;
134	part.y = y0 + part.yp * part.s;
135	//cout<<" am_length "<<am_length<<endl;
136	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
137	//cout<<"AM LAYER !!!! "<<part.PC<<endl;
138	//cout<< "Amorphous"<< endl;
139	part.x = part.x + part.xp * (crys.Cry_length - part.s);
140	part.y = part.y + part.yp * (crys.Cry_length - part.s);
141	part.namor = part.namor + 1; /////////////COUNT/////////////
142	part.effet = 2;
143	} else if ((part.x > (crys.xmax - crys.Alayer)) && (part.x < crys.xmax)) {
144	out = false;
145
146	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
147
148	part.namor = part.namor + 1; /////////////COUNT/////////////
149	part.effet = 2;
150	//cout<<"Fix HERE"<<endl;
151	}
152	return out;
153	}
154
155
156	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
157	//SI LES DEUX FONCTIONS D AVT RETURNE TRUE, ALORS ON CALCULE LES PARA CRITIQUE ET D AUTRE PARA (POINT C - E)
158	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
159
160
161
162	void SimCrys::parameters(SimCrys& sim)
163	{
164	//THE FIRST PART IS FOR THE (110) ORIENTATION
165	//cout << "The parameters function is use !!!!! "<< endl;
166	double c_v1;
167	double c_v2;
168
169
170	crys.xpcrit0 = sqrt(2.10e-9 * crys.eUm[crys.IS] / part.PC);
171	crys.Rcrit = part.PC / (2.e-6 * crys.eUm[crys.IS]) * crys.AI[crys.IS]; //if R>Rcritical=>no channeling is possible
172	crys.ratio = crys.Rcurv / crys.Rcrit;
173	crys.xpcrit = crys.xpcrit0 * (crys.Rcurv - crys.Rcrit) / crys.Rcurv;
174	c_v1 = 1.7;
175	c_v2 = -1.5;
176	if (crys.ratio <= 1) {
177	part.Ang_rms = c_v1 * 0.42 * crys.xpcrit0 * sin(1.4 * crys.ratio); //rms scattering
178	part.Ang_avr = c_v2 * crys.xpcrit0 * 0.05 * crys.ratio; //average angle reflection
179	part.Vcapt = 0.0;
180	} else if (crys.ratio <= 3) {
181	part.Ang_rms = c_v1 * 0.42 * crys.xpcrit0 * sin(1.571 * 0.3 * crys.ratio + 0.85); //rms scattering
182	part.Ang_avr = c_v2 * crys.xpcrit0 * (0.1972 * crys.ratio - 0.1472); //average angle reflection
183	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)
184	} else {
185	part.Ang_rms = c_v1 * crys.xpcrit0 * (1 / crys.ratio);
186	part.Ang_avr = c_v2 * crys.xpcrit0 * (1 - 1.6667 / crys.ratio);
187	part.Vcapt = 0.0007 * (crys.ratio - 0.7) / pow(part.PC, 0.2);
188	}
189	//cout <<crys.xpcrit<<endl;
190	if (crys.C_orient == 2) {
191	part.Ang_avr = part.Ang_avr * 0.93;
192	part.Ang_rms = part.Ang_rms * 1.05; // FOR (111) ORIENTATION
193	crys.xpcrit = crys.xpcrit * 0.98;
194	}
195	//cout <<crys.xpcrit<<endl;
196	}
197
198
199
200
201	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
202	//MAINTENANT ON FAIT UNE FONCTION POUR LA PARTIE I CF SHEMA!!! C EST LA PARTIE CHANNELING D OU LE NOM
203	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
204
205
206
207
208	void SimCrys::channel(SimCrys& sim)
209	{
210	//cout << "La fonction channel est utilisee!!!! "<<endl;
211
212
213	part.Chann = sqrt(crys.xpcrit * crys.xpcrit - part.xp_rel * part.xp_rel) / crys.xpcrit0; //probability of channeling/volume capture
214
215	//cout <<"Channeling Prob :" << part.Chann <<endl; //POUR REGARDER LA VALEUR!!!!!!!!!!!!!!
216	//cout <<"Channeling ANGLE :" << crys.tchan<<endl; //POUR REGARDER LA VALEUR!!!!!!!!!!!!!!
217
218
219	double n_atom(0.1); //probability for entering channeling near atomic planes
220	double Dxp(0); //change angle from channeling [mrad]
221	double TLdech1(0); //tipical dechanneling length(1) [m]
222	double tdech(0);
223	double Ldech(0); //Dechanneling. length
224	double Sdech(0);
225	//cout<<"Chann : "<<part.Chann<<endl;
226	if (random() <= part.Chann) {
227	TLdech1 = 0.0005 * part.PC * pow(1 - 1 / crys.ratio, 2); // calculate tipical dechanneling length(m)
228	if (random() <= n_atom) {
229	TLdech1 = 0.000002 * part.PC * pow(1. - 1. / crys.ratio, 2); // dechanneling length (m) for short crystal for high amplitude particles
230	}
231	//cout<<"@ energy " << part.PC <<" dechanneling length "<< TLdech1 << endl; //FAUDRA OTER CELA!!!!!!!!!!
232	part.Dechan = log(random());
233	//cout<<"dechannn " << part.Dechan <<endl; //FAUDRA OTER CELA!!!!!!!!!!
234	Ldech = -TLdech1 * part.Dechan; //careful: the dechanneling lentgh is along the trajectory of the particle -not along the longitudinal coordinate...
235	tdech = Ldech / crys.Rcurv;
236	Sdech = Ldech * cos(crys.miscut + 0.5 * tdech);
237	// 2 option if the particule channeling !!!
238	// (1) dechanneling
239	// (2) full channeling
240
241	// (1) :
242
243	if (Ldech < crys.L_chan) {
244	//cout <<"Dechanneling ! "<< endl;
245	Dxp = Ldech / crys.Rcurv; //change angle from channeling [mrad]
246
247	part.ndechann = part.ndechann + 1; /////////////COUNT/////////////
248	part.effet = 4;
249	part.x = part.x + Ldech * (sin(0.5 * Dxp + crys.miscut)); // trajectory at channeling exit
250	part.xp = part.xp + Dxp + (random_gauss() - 0.5) * crys.xpcrit * 0.5;
251	part.y = part.y + part.yp * Sdech;
252
253	part.x = part.x + 0.5 * (crys.s_length - Sdech) * part.xp;
254	part.y = part.y + 0.5 * (crys.s_length - Sdech) * part.yp;
255
256
257	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
258
259	part.PC = part.PC - 0.5 * crys.DES[crys.IS] * part.y; //energy loss to ionization [GeV]
260	part.x = part.x + 0.5 * (crys.s_length - Sdech) * part.xp;
261	part.y = part.y + 0.5 * (crys.s_length - Sdech) * part.yp;
262	}
263
264	// (2) :
265
266	else {
267	//cout <<"Channeling ! "<< endl;
268	Dxp = crys.L_chan / crys.Rcurv; //change angle [mrad]
269	part.x = part.x + crys.L_chan * (sin(0.5 * Dxp + crys.miscut)); //trajectory at channeling exit
270
271
272	// Removed dependence on incoming angle according to discussion with I.Yazynin on 24/10/2012
273	//part.xp = part.xp + Dxp + (random_gauss()-0.5)crys.xpcrit0.5; // [mrad]
274	part.xp = /part.xp/ + Dxp + (random_gauss() - 0.5) * crys.xpcrit * 0.5; // [mrad]
275
276	part.y = part.y + crys.s_length * part.yp;
277	part.PC = part.PC - 0.5 * crys.DES[crys.IS] * crys.Cry_length; // energy loss to ionization [GeV]
278	part.nchann = part.nchann + 1; /////////////COUNT/////////////
279	part.effet = 3;
280	}
281	} else {
282	//cout <<"Volume Reflection ! "<< endl;
283	part.nvrefl = part.nvrefl + 1; /////////////COUNT/////////////
284	part.effet = 5;
285	Dxp = 0.5 * (part.xp_rel / crys.xpcrit + 1) * part.Ang_avr;
286	part.xp = part.xp + Dxp + part.Ang_rms * random_gauss();
287	/*
288	next line new from sasha
289	part.xp=part.xp+0.45(part.xp/crys.xpcrit+1)part.Ang_avr;
290	valentina fix here
291	*/
292	part.x = part.x + 0.5 * crys.s_length * part.xp;
293	part.y = part.y + 0.5 * crys.s_length * part.yp;
294
295	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
296
297	part.x = part.x + 0.5 * crys.s_length * part.xp;
298	part.y = part.y + 0.5 * crys.s_length * part.yp;
299
300	}
301	}
302
303
304
305	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
306	//MAINTENANT ON FAIT UNE FONCTION POUR LA PARTIE II CF SHEMA!!! C EST LA PARTIE VOLUME REFLEXION D OU LE NOM
307	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
308
309
310
311
312	void SimCrys::reflection(SimCrys& sim)
313	{
314	//cout<<"The function reflection is used" <<endl;
315	double Dxp(0); //change angle from channeling [mrad]
316	double Lrefl(0); //distance of the reflection point [m]
317	double Srefl(0); //distance of the reflection point [m]
318	double TLdech2(0); //tipical dechanneling length(2) [m]
319	double tdech(0);
320	double Ldech(0); //Dechanneling. length
321	double Sdech(0);
322	double Red_S(0); //Reduced length (in case of dechanneling)
323	double Rlength(0); //Reduced length
324
325	Lrefl = (part.xp_rel) * crys.Rcurv;
326	Srefl = sin(part.xp) * Lrefl;
327
328
329	if ((Lrefl > 0) && (Lrefl < crys.Cry_length)) { // IE : If the VR point is inside the crystal!!!!
330
331	//cout<<"VRCapt : "<<part.Vcapt<<endl;
332	if ((random() > part.Vcapt) \|\| (crys.ZN == 0)) {
333	//cout<< "Volume Reflexion! "<<endl;
334	part.nvrefl = part.nvrefl + 1; /////////////COUNT/////////////
335	part.effet = 5;
336	part.x = part.x + part.xp * Srefl;
337	part.y = part.y + part.yp * Srefl;
338	Dxp = part.Ang_avr;
339	part.xp = part.xp + Dxp + part.Ang_rms * random_gauss();
340	part.x = part.x + 0.5 * part.xp * (crys.s_length - Srefl);
341	part.y = part.y + 0.5 * part.yp * (crys.s_length - Srefl);
342
343	if (crys.ZN > 0) {
344	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
345	}
346
347	part.x = part.x + 0.5 * part.xp * (crys.s_length - Srefl);
348	part.y = part.y + 0.5 * part.yp * (crys.s_length - Srefl);
349	} else {
350	//cout<< "Volume Capture! ";
351	part.nvcapt = part.nvcapt + 1; /////////////COUNT/////////////
352	part.effet = 6;
353	part.x = part.x + part.xp * Srefl;
354	part.y = part.y + part.yp * Srefl;
355	/*
356	TLdech2= 0.00011pow(part.PC,0.25)pow((1.-1./crys.ratio),2) //dechanneling length [m]
357	part.Dechan = log(1.-random())
358	Ldech = -TLdech2*part.Dechan
359	next 2 lines new from sasha - different dechanneling probability
360	*/
361	TLdech2 = 0.01 * part.PC * pow((1. - 1. / crys.ratio), 2); // typical dechanneling length [m]
362	Ldech = 0.005 * TLdech2 * pow((sqrt(0.01 - log(random())) - 0.1), 2); // (dechanneling length) [m]
363	tdech = Ldech / crys.Rcurv;
364	Sdech = Ldech * cos(part.xp + 0.5 * tdech);
365
366
367
368	if (Ldech < (crys.Cry_length - Lrefl)) { // for dechanneling part
369	//cout<< "Dechanneling! ";
370	part.ndechann = part.ndechann + 1; /////////////COUNT/////////////
371	part.effet = 4;
372	Dxp = Ldech / crys.Rcurv;
373	part.x = part.x + Ldech * (sin(0.5 * Dxp + part.xp)); //trajectory at channeling exit
374	part.y = part.y + Sdech * part.yp;
375	part.xp = part.xp + Dxp + (random_gauss() - 0.5) * crys.xpcrit * 0.5;
376	Red_S = crys.s_length - Srefl - Sdech;
377	// move in amorphous substance----------------------------
378	part.x = part.x + 0.5 * part.xp * Red_S;
379	part.y = part.y + 0.5 * part.yp * Red_S;
380
381	if (crys.ZN > 0) {
382	//cout<<"Amorphous! "<<endl; /////////////COUNT/////////////
383	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
384	}
385
386	part.x = part.x + 0.5 * part.xp * Red_S;
387	part.y = part.y + 0.5 * part.yp * Red_S;
388	//cout<<endl;
389	} else {
390	//cout<< "Volume Capture! "<<endl;
391	//Pas besoin d additionner 1 dans le compte de vol capt car deja fait!!!
392	Rlength = crys.Cry_length - Lrefl;
393	crys.tchan = Rlength / crys.Rcurv;
394	Red_S = Rlength * cos(part.xp + 0.5 * crys.tchan);
395	Dxp = (crys.Cry_length - Lrefl) / crys.Rcurv;
396	part.x = part.x + sin(0.5 * Dxp + part.xp) * Rlength; // trajectory at channeling exit
397	part.y = part.y + Red_S * part.yp;
398	part.xp = part.xp + Dxp + (random_gauss() - 0.5) * crys.xpcrit * 0.5; //[mrad]
399	}
400	}
401	}
402	// move in amorphous substance for big input angles---------------
403	else {
404	//cout<<"Amorphous! "<<endl;
405	part.namor = part.namor + 1; /////////////COUNT/////////////
406	part.effet = 2;
407	part.x = part.x + 0.5 * crys.s_length * part.xp;
408	part.y = part.y + 0.5 * crys.s_length * part.yp;
409
410	if (crys.ZN > 0) {
411	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
412	}
413
414	part.x = part.x + 0.5 * crys.s_length * part.xp;
415	part.y = part.y + 0.5 * crys.s_length * part.yp;
416	}
417
418	}
419
420
421
422	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
423	//THE FUNCTION THAT DO ALL THE SCHEMA USING THE PREVIOUS FUNCTION !!!!!!!!!!111
424	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
425
426
427
428
429	//called the FORTRAN cyrstal-particle routine from SixTrack
430	/*
431	void SimCrys::general(SimCrys& sim, const int& pas, char name, double emitx0, double emity0, double enom, string outputpath, double Mirror, double C_rotation, double C_aperture, double C_offset, double C_tilt, double Crystal_tilt)
432	{
433
434	long int max_npart = 1500000;
435
436	char name_coll[50] = name;
437
438	int NP = ;
439	double ENOM = enom*1e-3; // [GeV]
440	double EMITX0 = emitx0;
441	double EMITY0 = emity0;
442	double AX = , AY = , BX = , BY = ;
443
444	double X_IN[max_npart] ;
445	double XP_IN[max_npart];
446	double Y_IN[max_npart];
447	double YP_IN[max_npart];
448	double P_IN[max_npart];
449	double S_IN[max_npart] = ;
450
451	int flagsec[max_npart] = ;
452	bool dowrite_impact = true;
453	int name[max_npart] = ;
454
455	char C_MATERIAL[6];
456	double C_LENGTH = crys.s_length;
457	double C_ROTATION = C_rotation;
458	double C_APERTURE = C_aperture;
459	double C_OFFSET = C_offset;
460	double C_TILT[2] = {C_tilt, 0} ;
461
462	int LHIT[max_npart] = ; //position of the particle been hit
463	int PART_ABS[max_npart] = ;
464
465	double IMPACT[max_npart] = ;
466	double INDIV[max_npart] = ;
467	double LINT[max_npart] = ;
468
469	*/
470	/* -----------------
471	!++ Copy particle data to 1-dim array and go back to meters
472	*/
473	/*
474	for(j = 1, napx)
475	rcx(j) = (xv(1,j)-torbx(ie))/1d3
476	rcxp(j) = (yv(1,j)-torbxp(ie))/1d3
477	rcy(j) = (xv(2,j)-torby(ie))/1d3
478	rcyp(j) = (yv(2,j)-torbyp(ie))/1d3
479	rcp(j) = ejv(j)/1d3
480	rcs(j) = 0d0
481	part_hit_before(j) = part_hit(j)
482	rcx0(j) = rcx(j)
483	rcxp0(j) = rcxp(j)
484	rcy0(j) = rcy(j)
485	rcyp0(j) = rcyp(j)
486	rcp0(j) = rcp(j)
487	ejf0v(j) = ejfv(j)
488	!
489	!++ For zero length element track back half collimator length
490	!
491	if (stracki.eq.0.) then
492	rcx(j) = rcx(j) - 0.5d0c_lengthrcxp(j)
493	rcy(j)= rcy(j) - 0.5d0c_lengthrcyp(j)
494	else
495	Write(,)
496	& "ERROR: Non-zero length collimator!"
497	STOP
498	endif
499	flukaname(j) = ipart(j)+100*samplenumber
500	!
501	end do
502	!
503	!++ Do the collimation tracking
504	!
505	enom_gev = myenom*1d-3
506	!
507	!++ Allow
508
509
510
511
512	if (crys.IS == 0) {
513	C_MATERIAL = "CRY-Si";
514	} else if (crys.IS == 1) {
515	C_MATERIAL = "CRY-W";
516	} else if (crys.IS == 2) {
517	C_MATERIAL = "CRY-C";
518	} else if (crys.IS == 3) {
519	C_MATERIAL = "CRY-Ge";
520	}else{
521	cerr << "Material of the Crystal not fund !!!" << endl;
522	}
523
524
525
526	int count = 0;
527	string rest;
528	ifstream entree;
529	string name2;
530	name2 = outputpath + "/ascii_before.csv";
531	entree.open(name2.c_str());
532	if (entree) {
533	while (!entree.eof()) {
534	count = count + 1;
535
536	getline(entree, rest, ',');
537	X_IN[count] = atof(rest.c_str());
538	//cout<<" X " << part.x<<endl;
539
540	getline(entree, rest, ',');
541	Y_IN[count] = atof(rest.c_str());
542	//cout<<" Y " << part.y<<endl;
543
544	getline(entree, rest, ',');
545	XP_IN[count] = atof(rest.c_str());
546	//cout<<" XP " << part.xp<<endl;
547
548	getline(entree, rest, ',');
549	YP_IN[count] = atof(rest.c_str());
550	//cout<<" YP " << part.yp<<endl;
551
552	getline(entree, rest);
553	P_IN[count] = atof(rest.c_str());
554
555	}
556	}
557	entree.close();
558
559
560	//from the file "crystaltrack_dan.f", 08/2014, by Jianfeng Zhang @ LAL
561	COLLIMATE_CRY_(name_coll,C_MATERIAL, C_LENGTH,
562	C_ROTATION,
563	C_APERTURE, C_OFFSET, C_TILT,
564	X_IN, XP_IN, Y_IN,
565	YP_IN, P_IN, S_IN, NP, ENOM, LHIT,
566	PART_ABS, IMPACT, INDIV, LINT,
567	BX,BY,AX,
568	AY,EMITX0,EMITY0,
569	name,flagsec,dowrite_impact);
570	}
571
572	*/
573
574	/********************************************************************************
575	*
576	* original crystal routine
577	* *******************************************************************************/
578	void SimCrys::general(SimCrys& sim, const int& pas, string outputpath, double Mirror, double C_rotation, double C_aperture, double C_offset, double C_tilt, double Crystal_tilt)
579	{
580	srand ( time(NULL) );
581	double AV_xp0(0);
582	double AV_yp0(0);
583	double RMSxp0(0);
584	double RMSyp0(0);
585	double AV_xp1(0);
586	double AV_yp1(0);
587	double long RMSxp1(0);
588	double long RMSyp1(0);
589	double av_pc1(0);
590	double av_pc0(0);
591	double SIGxp0(0);
592	double SIGyp0(0);
593	double SIGxp1(0);
594	double SIGyp1(0);
595	double xp0(0);
596	double xp1(0);
597	double count(0);
598	double zpos(-0.1);
599	double xfin(0);
600	double yfin(0);
601	int const col(15);
602
603
604	//parameter to check whether the particle hit the crystal
605	char Proc[50] = " "; // process of the particle-crystal interaction
606	bool crossing = false; //the particle hit the collimator or not
607	int layerflag = 1; //check the particle exit the layer in "AM"
608	//-----------------------------------------------NOW THE PARAMETERS FOR THE CHANGE OF REFERENTIAL ------------------------------------------------
609	int mat=0;
610
611	double p0=0.0;
612	double zlm=0.0;
613	double shift=0.0;
614	double x_shift=0.0, xp_shift=0.0, s_shift=0.0; //coordinates after shift/rotation
615	double x_rot=0.0, xp_rot=0.0, s_rot=0.0;
616	double x_temp=0.0, xp_temp=0.0, s_temp=0.0; //all the _temp variables are used when you hit the cry from below
617	double tilt_int=0.0, x_int=0.0, xp_int=0.0, s_int=0.0;//all the _int variables are used when you hit the cry from below (int=interaction point)
618	double x00=0.0;
619	double z=0.0;
620	double z00=0.0;
621	double zp=0.0;
622	double dpop=0.0;
623	double s=0.0;
624	double a_eq=0.0, b_eq=0.0, c_eq=0.0, Delta=0.0;
625	int part_abs=-1;
626	double x=0.0;
627	double xp=0.0;
628	double y=0.0;
629	double yp=0.0;
630	double p=0.0; //be careful: [Gev]
631	double s_in=0.0;
632
633	// adding variables for the pencil beam. Variables in the absolute reference frame.
634
635	double x_in0=0.0;
636	double xp_in0=0.0;
637	double y_in0=0.0;
638	double yp_in0=0.0;
639	double p_in0=0.0; //be careful: [Gev]
640	double s_in0=0.0;
641	double s_impact=0.0;
642	double totcount(0);
643
644	double mirror=Mirror;
645	double tiltangle=0.0;
646
647
648	double c_length=0.0; //Length in m
649	double c_rotation(C_rotation) ; //Rotation angle vs vertical in radian initiated and fixed at 0. Can BE CHANGED IF NEED
650	double c_aperture(C_aperture) ; //Aperture in m
651	double c_offset(C_offset) ; //Offset in m
652	double c_tilt[2] = {C_tilt, 0} ; //Tilt in radian
653	double c_tilt0[2] ; //Tilt in radian
654
655
656	double xp_tangent=0.0;
657
658	double Cry_tilt(Crystal_tilt); //crystal tilt [rad] ///////////////////////////////////++++++++++++++++++++++++(IF NEEDED, put it into the input file !!! See later)+++++++++++++++////////////////
659
660
661
662	//----------------------------------END OF INITATING REFERNENTIAL PARAMETERS-------------------------------------------------------------------------
663	c_length = crys.Cry_length;
664
665	if (crys.IS == 0) {
666	mat = 8;
667	} else if (crys.IS == 1) {
668	mat = 9;
669	} else if (crys.IS == 2) {
670	mat = 10;
671	} else if (crys.IS == 3) {
672	mat = 11;
673	} else {
674	cout << "Material of the Crystal not fund !!!" << endl;
675	}
676
677
678	if (c_length > 0 ) {
679	p0 = part.PC;
680	c_tilt0[0] = c_tilt[0];
681	c_tilt0[1] = c_tilt[1];
682	tiltangle = c_tilt0[0];
683	}
684
685
686	//new..........................
687	// call scatin(p0) ??????????????
688
689	// scatin_(p0);
690
691	/* EVENT LOOP, initial distribution is here a flat distribution with
692	* xmin=x-, xmax=x+, etc. from the input file
693	*/
694	long int nhit = 0;
695	double fracab = 0;
696	// int n_chan = 0; // valentina :initialize to zero the counters for crystal effects
697	// int n_VR = 0;
698	// int n_amorphous = 0;
699
700
701	double impact = 0;
702	double indiv = 0;
703	double lint = 0;
704
705
706	// c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
707	//c WRITE(,) ITURN,ICOLL
708	/* do j = 1, nev
709	//c
710	impact(j) = -1.0;
711	lint(j) = -1.0;
712	indiv(j) = -1.0;
713
714
715	idx_proc = name(j)-100*samplenumber;
716
717	if (ITURN == 1){ // !daniele
718	bool_proc_old(idx_proc)=-1;} // !daniele
719	// elseif (bool_proc(idx_proc) != -1){
720	else{
721	bool_proc_old(idx_proc)=bool_proc(idx_proc); // !daniele
722	} // !daniele
723
724	PROC= "out"; // !the default process is 'out' !daniele
725	bool_proc(idx_proc)=-1; // !daniele
726	*/
727
728	int nabs = 0;
729
730	//ofstream sortie("results.csv");
731	//sortie << "xp1" <<","<< "xp0"<<","<< "x"<<","<< "y"<<endl;
732	ofstream sortieIco;
733	string name1;
734
735	//file with the particle coordinate after the crystal
736	name1 = outputpath + "/ascii_after.csv";
737
738	sortieIco.open(name1.c_str());
739
740	string rest;
741
742	ifstream entree;
743	string name2;
744	name2 = outputpath + "/ascii_before.csv";
745
746	entree.open(name2.c_str());
747	if (entree) {
748	while (!entree.eof()) {
749	count = count + 1;
750
751	getline(entree, rest, ',');
752	part.x = atof(rest.c_str());
753	//cout<<" X " << part.x<<endl;
754
755	getline(entree, rest, ',');
756	part.y = atof(rest.c_str());
757	//cout<<" Y " << part.y<<endl;
758
759	getline(entree, rest, ',');
760	part.xp = atof(rest.c_str());
761	//cout<<" XP " << part.xp<<endl;
762
763	getline(entree, rest, ',');
764	part.yp = atof(rest.c_str());
765	//cout<<" YP " << part.yp<<endl;
766
767	getline(entree, rest);
768	part.PC = atof(rest.c_str());
769
770	if (part.PC != 0) { // POUR OTER LA DERNIERE PART QUI VAUT 0 POUR TOUT
771
772	av_pc0 = av_pc0 + part.PC;
773
774	//part.x=random_gauss()*0.00115;
775	//part.xp=random_gauss()sqrt(51e-6);
776	//cout<<"X et Y avt :"<<part.x<<" , "<<part.y<<endl;
777
778	// IL FAUT PASSER DES CM EN M!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
779	//bdelta=tan(part.xp)*10;
780	part.x = (part.x); //-0.0045;
781
782
783	//bdelta=tan(part.yp)*10;
784	part.y = (part.y); //-0.0045;
785
786	//cout<<"X et Y apres :"<<part.x<<" , "<<part.y<<endl;
787	xp0 = part.xp;
788
789	//part.y=random_gauss()*sqrt(0.00095);
790	//part.yp=5*1e-6;
791
792	//part.PC=120.976;
793	//cout<<" X :"<<part.x<<" XP :"<<part.xp<<" Y :"<<part.y<<" YP :"<<part.yp<<" PC :"<<part.PC<<endl;
794	AV_xp0 = AV_xp0 + part.xp; //average x angle before impact
795	AV_yp0 = AV_yp0 + part.yp; //average y angle before impact
796	RMSxp0 = RMSxp0 + part.xp * part.xp; //rms x angle before impact
797	RMSyp0 = RMSyp0 + part.yp * part.yp; //rms y angle before impact
798
799	// NOW WE NEED TO CHANGE THE REFERENTIAL!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
800	//change the beam frame to the crystal frame !!!!!!!!!
801	// since the interaction between the particle and the crystal is observed in the
802	// crystal reference system
803
804	//initial coordinates in the beam frame
805	s = 0;
806	x = part.x;
807	xp = part.xp;
808	z = part.y;
809	zp = part.yp;
810	p = part.PC;
811
812	x_temp = 0;
813	x_int = 0;
814	x_rot = 0;
815	x_shift = 0;
816	s_temp = 0;
817	s_int = 0;
818	s_rot = 0;
819	s_shift = 0;
820	xp_int = 0;
821	xp_temp = 0;
822	xp_rot = 0;
823	xp_shift = 0;
824	shift = 0;
825	tilt_int = 0;
826	dpop = (p - p0) / p0; //momentum spread
827
828
829	/*---------------DANIELE-------------------------
830	! corrected position for variable association for FirstImpact.dat
831	! also for a vertical crystal case.
832	! Only x_in0 (i.e. b) have to be assigned after the change of reference frame
833	!-----------------------------------------------*/
834
835	s_in0 = s_in; //daniele
836	// x_in0(j) = x //!daniele
837	xp_in0 = xp; //daniele
838	y_in0 = z; //daniele
839	yp_in0 = zp; //daniele
840	p_in0 = p; //daniele
841
842	/----------------DANIELE------------------/
843
844
845	//**********************************************************************************
846
847	// ------------------------------transform particle coordinates to get into collimator coordinate system -----------------------------------------
848
849
850	// first check whether particle was lost before
851	totcount = totcount + 1;
852	if ((x < 99.0 * 1e-3) && (z < 99.0 * 1e-3)) { ////////////////////////// IF NUMERO 1//////////////////
853
854	// first do rotation into collimator frame
855
856	x = part.x * cos(c_rotation) + sin(c_rotation) * part.y;
857	z = part.y * cos(c_rotation) - sin(c_rotation) * part.x;
858	xp = part.xp * cos(c_rotation) + sin(c_rotation) * part.yp;
859	zp = part.yp * cos(c_rotation) - sin(c_rotation) * part.xp;
860
861	//cout<<"PARTICULE NON PERDU>>><<"<<endl;
862
863	//To define the correct value of the mirror for the install locations
864	// of the crystal in SPS & LHC, inner side or outside of the vacuum chamber????
865	// by Zhang @ CERN, 23/04/2014.....
866	//mirror = 1;
867
868
869	x = mirror * x;
870	xp = mirror * xp;
871
872	//Shift with opening and offset
873
874	x = x - c_aperture / 2 - c_offset;
875
876	// Include collimator tilt
877
878	if (tiltangle > 0) { ////////////////////////// IF NUMERO 2//////////////////
879	xp = xp - tiltangle;
880	} //////////////////////////FIN IF NUMERO 2//////////////////
881	else if (tiltangle < 0) { ////////////////////////// IF NUMERO 3//////////////////
882	x = x + sin(tiltangle) * c_length;
883	xp = xp - tiltangle;
884	} //////////////////////////FIN IF NUMERO 3//////////////////
885
886	if (Cry_tilt < 0) { ////////////////////////// IF NUMERO 4//////////////////
887
888	s_shift = s;
889	shift = crys.Rcurv * (1 - cos(Cry_tilt));
890	if (Cry_tilt < (-crys.cry_bend) ) { ////////////////////////// IF NUMERO 5//////////////////
891	shift = ( crys.Rcurv * ( cos ( - Cry_tilt) - cos( crys.cry_bend - Cry_tilt) ) );
892	} //////////////////////////FIN IF NUMERO 5//////////////////
893	x_shift = x - shift;
894	} //////////////////////////fin IF NUMERO 4//////////////////
895	else { //////////////////////////IF NUMERO 6//////////////////
896	//cout<< "CRY-TILT POSSITIVE ++ or nul"<<endl;
897	s_shift = s;
898	x_shift = x;
899	} //////////////////////////fin IF NUMERO 6//////////////////
900	// cout<<"debug - S shift" << s_shift<<endl;
901	//cout<<"debug - X shift" << x_shift <<endl;
902	//
903	// 2nd transformation: rotation
904	s_rot = x_shift * sin(Cry_tilt) + s_shift * cos(Cry_tilt);
905	x_rot = x_shift * cos(Cry_tilt) - s_shift * sin(Cry_tilt);
906	xp_rot = xp - Cry_tilt;
907	//cout<< "debug - S rot" << s_rot <<endl;
908	//cout<<"debug - X rot" << x_rot <<endl;
909	//cout<<"debug - XP rot" << xp_rot<<endl;
910
911	// 3rd transformation: drift to the new coordinate s=0
912	xp = xp_rot;
913	x = x_rot - xp_rot * s_rot;
914	z = z - zp * s_rot;
915	s = 0;
916	//cout<<"debug - S cryRF" << s_rot <<endl;
917	//cout<<"debug - X cryRF" << x_rot <<endl;
918	//cout<<"debug - XP cryRF" << xp_rot <<endl;
919
920
921
922	//********************************************************************************************************************
923	//----------------------------------NOW CHECK IF THE PARTICLE HIT the crystal---------------------------------------------------------------
924	//cout<<" X " << x<<endl;
925	//Particle hit the crystal
926	if (x >= 0) { ////////////////////////// IF NUMERO 7//////////////////
927	s_impact = s_in0; //!(for the first impact)
928
929	if(1){
930	cout<<"hit the cry entrance face"<<endl;
931	cout<<"impact at s = "<< s_impact<< ", x = "<<x_in0<<endl;
932	cout<<"with angle xp = "<<xp<<endl;
933	cout<<"s before: "<< s<<endl;
934	}
935
936	part.x = x;
937	part.xp = xp;
938	part.y = z;
939	part.yp = zp;
940	part.PC = p;
941
942	//------------------------------------------------------------------------------
943	// The proton crystal routine in Fortran version
944	//
945	//------------------------------------------------------------------------------
946
947	//Fortran crystal routine in SixTrack
948	// Modified by Jianfeng Zhang @ LAL, 30/04/2014
949	double crysCry_length=crys.Cry_length;
950	double crysRcurv=crys.Rcurv;
951	double crysC_xmax=crys.C_xmax;
952	double crysC_ymax=crys.C_ymax;
953	double crysAlayer=crys.Alayer;
954	int crysC_orient=crys.C_orient;
955	double crysmiscut = crys.miscut;
956	int IS =0;
957	double length =0.0;
958	IS = mat -7 ;
959	length = crys.Cry_length;
960	double eUmIS = crys.eUm[crys.IS];
961	double AIIS = crys.AI[crys.IS];
962
963
964
965	double crysparaDes[4];
966	crysparaDes[0]=crys.DES[0],crysparaDes[1]=crys.DES[1],crysparaDes[2]=crys.DES[2],crysparaDes[3]=crys.DES[3];
967	double crysparaDlyi[4];
968	crysparaDlyi[0]=crys.DLYI[0],crysparaDlyi[1]=crys.DLYI[1],crysparaDlyi[2]=crys.DLYI[2],crysparaDlyi[3]=crys.DLYI[3];
969	double crysparaDlri[4];
970	crysparaDlri[0]=crys.DLRI[0],crysparaDlri[1]=crys.DLRI[1],crysparaDlri[2]=crys.DLRI[2],crysparaDlri[3]=crys.DLRI[3];
971
972	double crysparaeUmIS[4];
973	crysparaeUmIS[0]=crys.eUm[0], crysparaeUmIS[1]=crys.eUm[1], crysparaeUmIS[2]=crys.eUm[2], crysparaeUmIS[3]=crys.eUm[3];
974	double crysparaAIIS[4];
975	crysparaAIIS[0]=crys.AI[0],crysparaAIIS[1]=crys.AI[1],crysparaAIIS[2]=crys.AI[2],crysparaAIIS[3]=crys.AI[3];
976
977	// cout << "before enter crystal routine...."<<endl;
978
979	// cryst_(&IS,&x,&xp,&z,&zp,&p,&length,&crysCry_length,&crysRcurv,&crysC_xmax,&crysC_ymax,&crysAlayer,&crysC_orient,&crysmiscut, &Proc,&layerflag, &part.Chann);
980
981	// cryst_(&length,&layerflag, crypara,partpara,crysparaDes,crysparaDlri,crysparaDlyi,crysparaeUmIS,crysparaAIIS,Proc);
982
983	cryst_(&length,&layerflag,Proc,&crys.L_chan,&crys.tchan,&crys.Alayer,&crys.ymin,&crys.ymax,&crys.Rcurv,
984	&crys.s_length,&crys.IS,&crys.NAM,
985	&crys.xpcrit0,&crys.xpcrit,&crys.Rcrit,&crys.ratio,&crys.C_orient,&crys.miscut,&crys.Cry_length,
986	&crys.ZN,&crys.C_xmax,&crys.C_ymax,
987	&part.x,&part.xp,&part.y,&part.yp,&part.PC,&part.s,&part.Chann,&part.xp_rel,&part.Ang_rms,
988	& part.Ang_avr,&part.Vcapt,&part.Dechan,crysparaDes,crysparaDlri,crysparaDlyi,crysparaeUmIS,crysparaAIIS);
989
990
991
992	//convert back the Fortran variables to C++ variables
993	//mat = IS +7;
994
995	// cout << "IS in Fortran = ..." << crys.IS << endl;
996
997	crys.IS = crys.IS - 1;
998
999	// cout << "IS in C++ = ..." << crys.IS << endl;
1000
1001	crys.Cry_length = length;
1002
1003	x=part.x;
1004	xp=part.xp;
1005	z=part.y;
1006	zp=part.yp;
1007	p=part.PC;
1008	s=part.s;
1009
1010
1011	/*
1012	crys.Rcurv = crysRcurv;
1013	crys.C_xmax = crysC_xmax;
1014	crys.C_ymax = crysC_ymax;
1015	crys.Alayer = crysAlayer;
1016	crys.C_orient = crysC_orient;
1017	crys.miscut = crysmiscut;
1018	*/
1019
1020	// cout<< "ICI AHHHHHHHHHHHHHHHHHHHHHHHHHHHH"<<endl;
1021	// //cout<<" xp rel and xp crit " << part.xp_rel<< " " << crys.xpcrit<<endl;
1022	// bases(sim);
1023	// //cout << part.xp<< " XP "<<endl;
1024	// bool crossing(cross(sim));
1025	// //cout<<part.xp-xp0<<endl;
1026	// //cout<<" xp rel and xp crit " << part.xp_rel<< " " << crys.xpcrit<<endl;
1027	// bool layering;
1028	// if (crossing != 0) { ////////////////////////// IF NUMERO 8//////////////////
1029	// layering = layer(sim);
1030	// //cout<<"Crossing !!!!"<<endl;
1031	//
1032	// if (layering == true) { ////////////////////////// IF NUMERO 9//////////////////
1033	// parameters(sim);
1034	// //cout<<" Layering =true " <<endl;
1035	//
1036	// if ((abs(part.xp_rel) < crys.xpcrit)) { ////////////////////////// IF NUMERO 10//////////////////
1037	// channel(sim);
1038	// //cout<< " CHANNELING " << endl;
1039	// } ////////////////////////// fin IF NUMERO 10//////////////////
1040	// else { ////////////////////////// IF NUMERO 11//////////////////
1041	// reflection(sim);
1042	// //cout<<" REFLECTION "<<endl;
1043	// } ////////////////////////// fin IF NUMERO 11//////////////////
1044	//
1045	//
1046	// } ////////////////////////// fin IF NUMERO 9//////////////////
1047	// } ////////////////////////// fin IF NUMERO 8//////////////////
1048
1049	//--------------------------------------------------------------------------------------------------
1050
1051
1052	// xp = part.xp;
1053
1054	//}
1055
1056	s = crys.Rcurv * sin(crys.cry_bend);
1057	zlm = crys.Rcurv * sin(crys.cry_bend);
1058	//cout<<"process:"<<PROC<<endl;
1059	//cout<<"s after"<< s<<endl;
1060	//cout<<"part.xp"<<part.xp<<endl;
1061
1062	//cout<<"----------------------1111111111111111111111111-----------------"<<endl;
1063
1064
1065	if (strncmp(Proc,"out",3)==0)
1066	crossing = false;
1067	else
1068	crossing = true;
1069
1070	if (crossing == true){
1071	nhit = nhit + 1;
1072	// lhit = 100000000*ie + ITURN;
1073	impact = x_in0;
1074	indiv = xp_in0;
1075	}
1076	}
1077	// Particle not hit the crystal
1078	////////////////////////// fin IF NUMERO 7//////////////////
1079	else { ////////////////////////// IF NUMERO 12//////////////////
1080
1081	//cout<<"OU LA ABBBBBBBBBBBB" <<endl;
1082	xp_tangent = sqrt((-2 * x * crys.Rcurv + x * x) / (crys.Rcurv * crys.Rcurv));
1083	//cout<<"RCURV "<<crys.Rcurv<<endl;
1084	//cout<<"tangent"<<xp_tangent<<"angle"<<xp<<endl;
1085
1086	//cout<<"s tan " << crys.Rcurv*sin(xp_tangent)<<endl;
1087	//cout<< " s tot "<< c_length<< crys.Rcurv*sin(crys.cry_bend)<<endl;
1088	if ( xp >= xp_tangent) { ////////////////////////// IF NUMERO 13//////////////////
1089
1090	//if it hits the crystal, calculate in which point and apply the
1091	//transformation and drift to that point
1092	a_eq = (1 + xp * xp);
1093	b_eq = 2 * xp * (x - crys.Rcurv);
1094	c_eq = -2 * x * crys.Rcurv + x * x;
1095	Delta = b_eq * b_eq - 4 * (a_eq * c_eq);
1096	s_int = (-b_eq - sqrt(Delta)) / (2 * a_eq);
1097	//cout<< " s int "<<s_int <<endl;
1098	if (s_int < crys.Rcurv * sin(crys.cry_bend)) { ////////////////////////// fin IF NUMERO 14//////////////////
1099	// transform to a new ref system:shift and rotate
1100
1101	x_int = xp * s_int + x;
1102	xp_int = xp;
1103	z = z + zp * s_int;
1104	x = 0;
1105	s = 0;
1106
1107
1108	// tilt_int=2*X_int/S_int
1109	tilt_int = s_int / crys.Rcurv;
1110	xp = xp - tilt_int;
1111
1112	//---------------------------------------------------------------------------------------
1113	// Fortran Cyrstal routine for the protons
1114	//
1115	//---------------------------------------------------------------------------------------
1116
1117	part.x = x;
1118	part.xp = xp;
1119	part.y = z;
1120	part.yp = zp;
1121	part.PC = p;
1122	crys.Cry_length = crys.Cry_length - (tilt_int * crys.Rcurv);
1123
1124
1125
1126	/*
1127
1128	//Fortran crystal routine in SixTrack
1129	// Modified by Jianfeng Zhang @ LAL, 30/04/2014
1130	CRYST_(mat-7,x,xp,z,zp,p,cry_length);
1131	*/
1132
1133
1134	//------------------------------------------------------------------------------
1135	// The proton crystal routine in Fortran version
1136	//
1137	//------------------------------------------------------------------------------
1138
1139	//Fortran crystal routine in SixTrack
1140	// Modified by Jianfeng Zhang @ LAL, 30/04/2014
1141	double crysCry_length=crys.Cry_length;
1142	double crysRcurv=crys.Rcurv;
1143	double crysC_xmax=crys.C_xmax;
1144	double crysC_ymax=crys.C_ymax;
1145	double crysAlayer=crys.Alayer;
1146	int crysC_orient=crys.C_orient;
1147	double crysmiscut = crys.miscut;
1148	int IS =0;
1149	double length =0.0;
1150	IS = mat -7 ;
1151	length = crys.Cry_length;
1152	double eUmIS = crys.eUm[crys.IS];
1153	double AIIS = crys.AI[crys.IS];
1154
1155
1156	double crysparaDes[4];
1157	crysparaDes[0]=crys.DES[0],crysparaDes[1]=crys.DES[1],crysparaDes[2]=crys.DES[2],crysparaDes[3]=crys.DES[3];
1158	double crysparaDlyi[4];
1159	crysparaDlyi[0]=crys.DLYI[0],crysparaDlyi[1]=crys.DLYI[1],crysparaDlyi[2]=crys.DLYI[2],crysparaDlyi[3]=crys.DLYI[3];
1160	double crysparaDlri[4];
1161	crysparaDlri[0]=crys.DLRI[0],crysparaDlri[1]=crys.DLRI[1],crysparaDlri[2]=crys.DLRI[2],crysparaDlri[3]=crys.DLRI[3];
1162
1163	double crysparaeUmIS[4];
1164	crysparaeUmIS[0]=crys.eUm[0], crysparaeUmIS[1]=crys.eUm[1], crysparaeUmIS[2]=crys.eUm[2], crysparaeUmIS[3]=crys.eUm[3];
1165	double crysparaAIIS[4];
1166	crysparaAIIS[0]=crys.AI[0],crysparaAIIS[1]=crys.AI[1],crysparaAIIS[2]=crys.AI[2],crysparaAIIS[3]=crys.AI[3];
1167
1168
1169
1170
1171
1172	// cout << "before enter crystal routine...."<<endl;
1173
1174	cryst_(&length,&layerflag,Proc,&crys.L_chan,&crys.tchan,&crys.Alayer,&crys.ymin,&crys.ymax,&crys.Rcurv,
1175	&crys.s_length,&crys.IS,&crys.NAM,
1176	&crys.xpcrit0,&crys.xpcrit,&crys.Rcrit,&crys.ratio,&crys.C_orient,&crys.miscut,&crys.Cry_length,
1177	&crys.ZN,&crys.C_xmax,&crys.C_ymax,
1178	&part.x,&part.xp,&part.y,&part.yp,&part.PC,&part.s,&part.Chann,&part.xp_rel,&part.Ang_rms,
1179	& part.Ang_avr,&part.Vcapt,&part.Dechan,crysparaDes,crysparaDlri,crysparaDlyi,crysparaeUmIS,crysparaAIIS);
1180
1181	//convert back the Fortran variables to C++ variables
1182	//mat = IS +7;
1183	crys.Cry_length = length;
1184	crys.IS = crys.IS - 1;
1185
1186
1187
1188	x=part.x;
1189	xp=part.xp;
1190	z=part.y;
1191	zp=part.yp;
1192	p=part.PC;
1193	s=part.s;
1194
1195	/*
1196	crys.Rcurv = crysRcurv;
1197	crys.C_xmax = crysC_xmax;
1198	crys.C_ymax = crysC_ymax;
1199	crys.Alayer = crysAlayer;
1200	crys.C_orient = crysC_orient;
1201	crys.miscut = crysmiscut;
1202	*/
1203
1204	/*
1205	bases(sim);
1206
1207	//cout<<"part.xp"<<part.xp<<endl;
1208
1209	bool crossing(cross(sim));
1210	//cout<<"Crossing !!!!"<<endl;
1211	bool layering(layer(sim));
1212
1213	if (layering == true) {
1214	//cout<<" Layering =true " <<endl;
1215	parameters(sim);
1216	if ((abs(part.xp_rel) < crys.xpcrit)) {
1217	channel(sim);
1218	//cout<< " CHANNELING " << endl;
1219	} else {
1220	reflection(sim);
1221	//cout<< " REFLECTION " << endl;
1222	}
1223
1224	}
1225	*/
1226	//---------------------------------------------------------------------------------------
1227
1228
1229
1230	//cout<<"part.xp"<<part.xp<<endl;
1231	//cout<<"----------------------2222222222222222222-----------------"<<endl;
1232
1233	s = crys.Rcurv * sin(crys.cry_bend - tilt_int);
1234	zlm = crys.Rcurv * sin(crys.cry_bend - tilt_int);
1235
1236
1237	if(layerflag == 0){
1238	part.namor = part.namor + 1; /////////////COUNT/////////////
1239	part.effet = 2;
1240	}
1241	if (strncmp(Proc,"out",3)==0)
1242	crossing = false;
1243	else
1244	crossing = true;
1245
1246	if (crossing == true) { ////////////////////////// IF NUMERO 15//////////////////
1247	x_rot = x_int;
1248	s_rot = s_int;
1249	xp_rot = xp_int;
1250	s_shift = s_rot * cos(-Cry_tilt) + x_rot * sin(-Cry_tilt);
1251	x_shift = -s_rot * sin(-Cry_tilt) + x_rot * cos(-Cry_tilt);
1252	xp_shift = xp_rot + Cry_tilt;
1253	if (Cry_tilt < 0) { ////////////////////////// IF NUMERO 16//////////////////
1254	s_impact = s_shift;
1255	x_in0 = x_shift + shift;
1256	xp_in0 = xp_shift;
1257	} ////////////////////////// fin IF NUMERO 16//////////////////
1258	else { ////////////////////////// IF NUMERO 17//////////////////
1259	x_in0 = x_shift;
1260	s_impact = s_shift;
1261	xp_in0 = xp_shift;
1262	} ////////////////////////// fin IF NUMERO 17//////////////////
1263
1264	nhit = nhit + 1;
1265	//lhit = 100000000*ie + ITURN;
1266	impact = x_in0;
1267	indiv = xp_in0;
1268	} ////////////////////////// fin IF NUMERO 15//////////////////
1269
1270
1271
1272	x_temp = x;
1273	s_temp = s;
1274	xp_temp = xp;
1275	s = s_temp * cos(-tilt_int) + x_temp * sin(-tilt_int);
1276	x = -s_temp * sin(-tilt_int) + x_temp * cos(-tilt_int);
1277	xp = xp_temp + tilt_int;
1278
1279	// 2nd: shift back the 2 axis
1280	x = x + x_int;
1281	s = s + s_int;
1282
1283	} ////////////////////////// fin IF NUMERO 14//////////////////
1284	// Finish if (s_int < crys.Rcurv * sin(crys.cry_bend))
1285	else { ///////////////////////////////////////////////PROBLEM : TOO MUCH OF PARTICULES WITHOUT ANY CHANGE OF ANGLE///////////THE PROB IS HERE WITH THOSES TWO "ELSE"/////////////////////////
1286
1287	s = crys.Rcurv * sin(crys.Cry_length / crys.Rcurv);
1288	x = x + s * xp;
1289	z = z + s * zp;
1290	++part.nout;
1291	//cout<< "the particule does not hit the crystal (Trop basse dans neg.)"<<endl;
1292	}// Finish else if (s_int < crys.Rcurv * sin(crys.cry_bend))
1293
1294	}else {
1295	s = crys.Rcurv * sin(crys.Cry_length / crys.Rcurv);
1296	x = x + s * xp;
1297	z = z + s * zp;
1298	++part.nout;
1299	//cout<< "the particule does not hit the crystal (negatif et plus petit angle qu xp_tang)"<<endl;
1300	}
1301	//} ////////////////////////// fin IF NUMERO 13//////////////////
1302	// Finish else if ( xp >= xp_tangent)
1303
1304	}////////////////////// 12
1305	// Finish particle not hit the crystal else (x < 0)
1306
1307
1308
1309
1310	//**************************************************************************************************************************
1311	//trasform back from the crystal to the collimator reference system
1312	// 1st: un-rotate the coordinates
1313
1314
1315	x_rot = x;
1316	s_rot = s;
1317	xp_rot = xp;
1318	//cout<< "debug - S Cry RF 2" , s_rot <<endl;
1319	//cout<<"debug - X Cry RF 2" , x_rot <<endl;
1320	//cout<<"debug - XP Cry RF 2" , xp_rot <<endl;
1321	s_shift = s_rot * cos(-Cry_tilt) + x_rot * sin(-Cry_tilt);
1322	x_shift = -s_rot * sin(-Cry_tilt) + x_rot * cos(-Cry_tilt);
1323	xp_shift = xp_rot + Cry_tilt;
1324	// 2nd: shift back the reference frame
1325	if (Cry_tilt < 0) { ////////////////////////// IF NUMERO 18//////////////////
1326	s = s_shift;
1327	x = x_shift + shift;
1328	xp = xp_shift;
1329	} ////////////////////////// fin IF NUMERO 18//////////////////
1330	else {
1331	x = x_shift;
1332	s = s_shift;
1333	xp = xp_shift;
1334	}
1335	// 3rd: shift to new S=Length position
1336	x = xp * (c_length - s) + x;
1337	z = zp * (c_length - s) + z;
1338	s = c_length;
1339
1340
1341	nabs = 0;
1342	part.effet = nabs;
1343	//cout<< "debug - S Coll RF 2" , s_rot <<endl;
1344	//cout<<"debug - X Coll RF 2" , x_rot <<endl;
1345	//cout<<"debug - XP Coll RF 2" , xp_rot <<endl;
1346
1347	//cout<<"X1_coll"<<x<<"Z1_coll"<<z<<"XP1_coll"<<xp<<"ZP1_coll"<<zp <<"s" <<s<<endl;
1348
1349
1350	if(strncmp(Proc,"out",3)==0){
1351	part.nout = part.nout + 1; /////////////COUNT/////////////
1352	part.effet = 1;
1353	}
1354	if(layerflag == 0){
1355	part.namor = part.namor + 1; /////////////COUNT/////////////
1356	part.effet = 2;
1357	}
1358	if(strncmp(Proc,"AM",2)==0){
1359	part.namor = part.namor + 1; /////////////COUNT/////////////
1360	part.effet = 2;
1361	}
1362	if(strncmp(Proc,"CH",2)==0){
1363	part.nchann = part.nchann + 1;
1364	part.effet = 3;
1365	}
1366	if(strncmp(Proc,"DC",2)==0){
1367	part.ndechann = part.ndechann + 1; /////////////COUNT/////////////
1368	part.effet = 4;;
1369	}
1370	if(strncmp(Proc,"VR",2)==0){
1371	part.nvrefl = part.nvrefl + 1; /////////////COUNT/////////////
1372	part.effet = 5;
1373	}
1374	if(strncmp(Proc,"VC",2)==0){
1375	part.nvcapt = part.nvcapt + 1; /////////////COUNT/////////////
1376	part.effet = 6;}
1377
1378	// ?????????????????????????????????????
1379	if(part.effet == 2) {
1380	part_abs = 0;
1381	}else {
1382	part_abs = 1;
1383	}
1384
1385
1386	//=========================================================================================================================
1387	// Transform back to particle coordinates with opening and offset
1388
1389	if (part_abs == 0) { ////////////////////////// IF NUMERO 19//////////////////
1390	//
1391	// Include collimator tilt
1392
1393	if (tiltangle > 0) {
1394	x = x + tiltangle * c_length;
1395	xp = xp + tiltangle;
1396	} else if (tiltangle < 0) {
1397	x = x + tiltangle * c_length;
1398	xp = xp + tiltangle;
1399
1400	x = x - sin(tiltangle) * c_length;
1401	}
1402
1403	// Transform back to particle coordinates with opening and offset
1404
1405	z00 = z;
1406	x00 = x + mirror * c_offset;
1407	x = x + c_aperture / 2 + mirror * c_offset;
1408
1409	//+ Now mirror at the horizontal axis for negative X offset
1410
1411	x = mirror * x;
1412	xp = mirror * xp;
1413
1414	// Last do rotation into collimator frame
1415
1416	part.x = x * cos((-1) * c_rotation) + z * sin((-1) * c_rotation);
1417	part.y = z * cos((-1) * c_rotation) - x * sin((-1) * c_rotation);
1418	part.xp = xp * cos((-1) * c_rotation) + zp * sin((-1) * c_rotation);
1419	part.yp = zp * cos((-1) * c_rotation) - xp * sin((-1) * c_rotation);
1420
1421
1422
1423	//**********************************************************************************
1424	//For output.... by Zhang @ CERN, 24/04/2014
1425	//p_in = (1 + dpop) * p0;
1426	//From Dianiele
1427	//p_in = p;
1428	//s_in = s_in + s;
1429	part.PC = p;
1430	part.s = part.s + s;
1431
1432	if (part.effet == 1) { ////////////////////////// IF NUMERO 21//////////////////
1433
1434	part.x = 99.99e-3;
1435	part.y = 99.99e-3;
1436	cout<< "Mean that the particules is lost"<<endl;
1437	} ////////////////////////// fin IF NUMERO 21//////////////////
1438	} ////////////////////////// fin IF NUMERO 19//////////////////
1439
1440
1441	// End of check for particles not being lost before (see @330)
1442
1443	//} ////////////////////////// fin IF NUMERO 12//////////////////
1444	// End of loop over all particles
1445
1446
1447
1448	} //collimator with length = 0 ////////////////////////// fin IF NUMERO 1//////////////////
1449
1450	// =================================================================================FIN DE RETOUR AU COORDONEE DE ICOSIM====================================================================
1451
1452	part.xp=xp;
1453	part.yp=zp;
1454	part.y=z;
1455	part.x=x;
1456	part.PC=p;
1457
1458	xp1 = part.xp;
1459
1460
1461	//if(xp1-xp0!=0){
1462	//sortie << xp1 << "," << xp0<<"," << xfin<<"," << yfin<<endl;
1463	//}
1464	//cout<<"ICIIIIIIIIIIIIIIIIIIIIIIIII"<<part.PC<<endl;
1465
1466
1467	AV_xp1 = AV_xp1 + part.xp; //average x angle after interaction
1468	AV_yp1 = AV_yp1 + part.yp; //average y angle after interaction
1469	RMSxp1 = RMSxp1 + part.xp * part.xp; //rms x angle after interaction
1470	RMSyp1 = RMSyp1 + part.xp * part.xp; //rms y angle after interaction
1471	av_pc1 = av_pc1 + part.PC;
1472
1473
1474	cout<<" X apres :"<<part.x<<" XP apres :"<<part.xp<<" Y apres :"<<part.y<<" YP apres :"<<part.yp<<" PC apres :"<<part.PC<<endl;
1475	//cout<<endl;
1476	//sortieIco.setf(ios::scientific);
1477	sortieIco << part.x << "," << part.y << "," << part.xp << "," << part.yp << "," << part.PC << endl;
1478	//}//end of for....
1479
1480	}
1481	} //End of while
1482
1483
1484	} else {
1485	cerr<< "We can not open the file !" << endl;
1486	}
1487	entree.close();
1488
1489
1490
1491	AV_xp0 = AV_xp0 / count;
1492	AV_yp0 = AV_yp0 / count;
1493
1494	RMSxp0 = sqrt(RMSxp0 / (count));
1495	RMSyp0 = sqrt(RMSyp0 / (count));
1496
1497	SIGxp0 = sqrt((RMSxp0 * RMSxp0) - (AV_xp0 * AV_xp0));
1498	SIGyp0 = sqrt((RMSyp0 * RMSyp0) - (AV_yp0 * AV_yp0));
1499
1500	AV_xp1 = AV_xp1 / count;
1501	AV_yp1 = AV_yp1 / count;
1502	RMSxp1 = sqrt(RMSxp1 / (count));
1503	RMSyp1 = sqrt(RMSyp1 / (count));
1504
1505	SIGxp1 = sqrt((RMSxp1 * RMSxp1) - (AV_xp1 * AV_xp1));
1506	SIGyp1 = sqrt((RMSyp1 * RMSyp1) - (AV_yp1 * AV_yp1));
1507	av_pc0 = av_pc0 / count;
1508	av_pc1 = av_pc1 / count;
1509	//--------------------------------------------------------------------------
1510	//-------------write a summary file-----------------------------------------
1511	//cout<<endl;
1512	//cout<<endl;
1513	//cout<<endl;
1514
1515	//cout<<"avg xp_in: "<< AV_xp0 <<"// avg xp_out: "<< AV_xp1<<endl;
1516	//cout<<"rms xp_in: "<< RMSxp0 <<"// rms xp_out: "<< RMSxp1<<endl;
1517	//cout<<"sigma xp_in: "<< SIGxp0 <<"// sigma xp_out: "<< SIGxp1<<endl;
1518	//cout<<"avg pc_in: "<< av_pc0 <<"// avg pc_out: "<< av_pc1<<endl;
1519
1520
1521
1522	sortieIco.close();
1523	// cout<<endl;
1524	// affiche();
1525	file_out(pas, outputpath);
1526
1527	//cout<<endl;
1528	//cout<<endl;
1529	//cout<<endl;
1530	}
1531
1532
1533
1534	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1535	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1536	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1537	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1538
1539	//WE WRITE SOME FUNCTIONS IMPORTANT FUNCTION NOW
1540
1541	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1542	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1543	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1544	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1545
1546
1547	void SimCrys::move_am(int IS, int NAM, double DZ, double DEI, double DLY, double DLR, double& xp, double& yp, double& PC)
1548	{
1549	//cout<<"The move_am function is use !!!!! "<<endl;
1550	double DYA(0);
1551	double Wp(0);
1552	if (NAM == 0) {
1553	return;
1554	}
1555	xp = xp * 1000;
1556	yp = yp * 1000;
1557	//cout << "xp avt: "<<xp<<" yp avt: "<<yp<<endl;
1558	//DEI ---> dE/dx (stoping energy)
1559	PC = PC - DEI * DZ; //energy lost beacause of ionization process[GeV]
1560
1561	// Coulomb Scattering
1562	//DYA ---> rms of coloumb scattering
1563	DYA = (14 / PC) * sqrt(DZ / DLR); //rms scattering (mrad)
1564	xp = xp + DYA * random_gauss_cut();
1565	yp = yp + DYA * random_gauss_cut();
1566
1567	//Elastic scattering
1568	//cout<<"Plusieurs possibilitees : "<<endl;
1569	if (random() <= (DZ / crys.DLAI[IS])) {
1570	//cout<<"Poss. 1 !!!! "<<endl;
1571	xp = xp + 196 * random_gauss_cut() / PC; //angle elast. scattering in R plane
1572	yp = yp + 196 * random_gauss_cut() / PC;
1573	}
1574
1575	//Diffraction interaction
1576	if (random() <= (DZ / (DLY * 6.143))) {
1577	//cout<<"Poss. 2 !!!! "<<endl;
1578	xp = xp + 257 * random_gauss_cut() / PC; // angle elast. scattering in R plane[mr]
1579	yp = yp + 257 * random_gauss_cut() / PC;
1580	Wp = random();
1581	PC = PC - 0.5 * pow(0.3 * PC, Wp); // m0*c = 1 GeV/c for particle
1582
1583	}
1584
1585	//Inelastic interaction
1586	if (random() <= DZ / DLY) {
1587	//cout<<" Absorbed !! "<< endl;
1588	}
1589	//cout<<"XP APRES : "<<xp<<" YP APRES : "<<yp<<endl;
1590	xp = xp / 1000;
1591	yp = yp / 1000;
1592	}
1593
1594
1595	double SimCrys::random()
1596	{
1597	double scale = RAND_MAX;
1598	double base = rand() / scale;
1599	double fine = rand() / scale;
1600	return base + fine / scale;
1601	}
1602
1603	double SimCrys::random_dist()
1604	{
1605	double scale = RAND_MAX;
1606	double base = rand() / scale;
1607	double fine = rand() / scale;
1608	return (-2.5 * 1e-4 + (3.5 * 1e-4 ) * (base + fine / scale));
1609	}
1610
1611	double SimCrys::random_gauss()
1612	{
1613	const double PI (4.0 * atan(1.0));
1614	double U(random());
1615	double V(random());
1616	return (sqrt(-2 * log(U)) * cos(2 * PI * V));
1617	}
1618
1619	double SimCrys::random_gauss_cut()
1620	{
1621	double resu;
1622	do {
1623	resu = random_gauss();
1624	} while (abs(resu) >= 1);
1625	return resu;
1626	}
1627
1628	double SimCrys::random_gauss_dist() //POUR faire varier la moyenne de la gaussienne il faut mettre en argu la sim& !!!!!
1629	{
1630	const double PI (4.0 * atan(1.0));
1631	double U(random());
1632	double V(random());
1633	//moy=random();
1634	return 1e-4 * (sqrt(-2 * log(U)) * cos(2 * PI * V)); //+(0.0002*moy); //POUR faire varier la moyenne de la gaussienne aussi !!!!!
1635	}
1636
1637	/////////////////////////////////*******************************************************+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1638	//POUR 409 :
1639
1640	double SimCrys::random_gauss_dist_409() //POUR faire varier la moyenne de la gaussienne il faut mettre en argu la sim& !!!!!
1641	{
1642	const double PI (4.0 * atan(1.0));
1643	double U(random());
1644	double V(random());
1645	//moy=random();
1646	return (8.939203e-06) * (sqrt(-2 * log(U)) * cos(2 * PI * V)) + (-7.183296e-08); //POUR faire varier la moyenne de la gaussienne aussi !!!!!
1647	}
1648
1649	///////////////////////////////////************************************************************++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1650
1651	void SimCrys::affiche()
1652	{
1653	cout << "Crystal Curved Length [m] :" << crys.Cry_length << endl;
1654	cout << "Crystal Curvature Radius [m] :" << crys.Rcurv << endl;
1655	cout << "Crystal Bending Angle [urad] :" << crys.cry_bend << " , " << crys.Cry_length / crys.Rcurv << endl;
1656	cout << "Critical Radius :" << crys.Rcrit << endl;
1657	cout << "Ratio :" << crys.ratio << endl;
1658	cout << "Critical Angle for straight Crystal :" << crys.xpcrit0 << endl;
1659	cout << "Critical Angle for curved Crystal :" << crys.xpcrit << endl;
1660	cout << "Angle average reflexion : " << part.Ang_avr << endl;
1661	cout << "Full Channeling : " << crys.Cry_length / crys.Rcurv << endl;
1662	cout << endl;
1663	cout << endl;
1664	cout << "N.AMORPH :" << part.namor << endl;
1665	cout << "N.DECHANNELING:" << part.ndechann << endl;
1666	cout << "N.CHANNELING :" << part.nchann << endl;
1667	cout << "N.OUT :" << part.nout << endl;
1668	cout << "N.VREFLECTION :" << part.nvrefl << endl;
1669	cout << "N.VCAPTURE :" << part.nvcapt << endl;
1670
1671	//cout << "Si la somme ne donne pas le nbre de part., c est normale !! En effet, il est possible que plusieurs evenement se produisent!!!";
1672	}
1673
1674	void SimCrys::file_out(const int& pas, string outputpath)
1675	{
1676
1677	ofstream out;
1678	string file;
1679	file = outputpath + "/crystal_output.out";
1680
1681	out.open(file.c_str(), ios::out \| ios::app);
1682
1683	if (out.fail()) {
1684	cerr << "Warning: problem writing in the file " << file << "!" << endl;
1685	} else {
1686
1687	out << endl;
1688	out << "Crystal Curved Length [m] :" << crys.Cry_length << endl;
1689	out << "Crystal Curvature Radius [m] :" << crys.Rcurv << endl;
1690	out << "Crystal Bending Angle [urad] :" << crys.cry_bend << " , " << crys.Cry_length / crys.Rcurv << endl;
1691	out << "Critical Radius :" << crys.Rcrit << endl;
1692	out << "Ratio :" << crys.ratio << endl;
1693	out << "Critical Angle for straight Crystal :" << crys.xpcrit0 << endl;
1694	out << "Critical Angle for curved Crystal :" << crys.xpcrit << endl;
1695	out << "Angle average reflexion : " << part.Ang_avr << endl;
1696	out << "Full Channeling : " << crys.Cry_length / crys.Rcurv << endl;
1697	out << endl;
1698	out << endl;
1699	out << "N.AMORPH :" << part.namor << endl;
1700	out << "N.DECHANNELING:" << part.ndechann << endl;
1701	out << "N.CHANNELING :" << part.nchann << endl;
1702	out << "N.OUT :" << part.nout << endl;
1703	out << "N.VREFLECTION :" << part.nvrefl << endl;
1704	out << "N.VCAPTURE :" << part.nvcapt << endl << endl << endl;
1705
1706	out.close();
1707	}
1708
1709	}
1710

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

Download in other formats: