Context Navigation

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

Last change on this file since 15 was 15, checked in by zhangj, 10 years ago
(1) Modify the bug to connect to the Fortran files. (2) Fix the bug of the definition type of the random functions used in Fortran file crystal_dan_FINAL....f.
File size: 62.9 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
1004	/*
1005	crys.Rcurv = crysRcurv;
1006	crys.C_xmax = crysC_xmax;
1007	crys.C_ymax = crysC_ymax;
1008	crys.Alayer = crysAlayer;
1009	crys.C_orient = crysC_orient;
1010	crys.miscut = crysmiscut;
1011	*/
1012
1013	// cout<< "ICI AHHHHHHHHHHHHHHHHHHHHHHHHHHHH"<<endl;
1014	// //cout<<" xp rel and xp crit " << part.xp_rel<< " " << crys.xpcrit<<endl;
1015	// bases(sim);
1016	// //cout << part.xp<< " XP "<<endl;
1017	// bool crossing(cross(sim));
1018	// //cout<<part.xp-xp0<<endl;
1019	// //cout<<" xp rel and xp crit " << part.xp_rel<< " " << crys.xpcrit<<endl;
1020	// bool layering;
1021	// if (crossing != 0) { ////////////////////////// IF NUMERO 8//////////////////
1022	// layering = layer(sim);
1023	// //cout<<"Crossing !!!!"<<endl;
1024	//
1025	// if (layering == true) { ////////////////////////// IF NUMERO 9//////////////////
1026	// parameters(sim);
1027	// //cout<<" Layering =true " <<endl;
1028	//
1029	// if ((abs(part.xp_rel) < crys.xpcrit)) { ////////////////////////// IF NUMERO 10//////////////////
1030	// channel(sim);
1031	// //cout<< " CHANNELING " << endl;
1032	// } ////////////////////////// fin IF NUMERO 10//////////////////
1033	// else { ////////////////////////// IF NUMERO 11//////////////////
1034	// reflection(sim);
1035	// //cout<<" REFLECTION "<<endl;
1036	// } ////////////////////////// fin IF NUMERO 11//////////////////
1037	//
1038	//
1039	// } ////////////////////////// fin IF NUMERO 9//////////////////
1040	// } ////////////////////////// fin IF NUMERO 8//////////////////
1041
1042	//--------------------------------------------------------------------------------------------------
1043
1044
1045	// xp = part.xp;
1046
1047	//}
1048
1049	s = crys.Rcurv * sin(crys.cry_bend);
1050	zlm = crys.Rcurv * sin(crys.cry_bend);
1051	//cout<<"process:"<<PROC<<endl;
1052	//cout<<"s after"<< s<<endl;
1053	//cout<<"part.xp"<<part.xp<<endl;
1054
1055	//cout<<"----------------------1111111111111111111111111-----------------"<<endl;
1056
1057
1058	if (strncmp(Proc,"out",3)==0)
1059	crossing = false;
1060	else
1061	crossing = true;
1062
1063	if (crossing == true){
1064	nhit = nhit + 1;
1065	// lhit = 100000000*ie + ITURN;
1066	impact = x_in0;
1067	indiv = xp_in0;
1068	}
1069	}
1070	// Particle not hit the crystal
1071	////////////////////////// fin IF NUMERO 7//////////////////
1072	else { ////////////////////////// IF NUMERO 12//////////////////
1073
1074	//cout<<"OU LA ABBBBBBBBBBBB" <<endl;
1075	xp_tangent = sqrt((-2 * x * crys.Rcurv + x * x) / (crys.Rcurv * crys.Rcurv));
1076	//cout<<"RCURV "<<crys.Rcurv<<endl;
1077	//cout<<"tangent"<<xp_tangent<<"angle"<<xp<<endl;
1078
1079	//cout<<"s tan " << crys.Rcurv*sin(xp_tangent)<<endl;
1080	//cout<< " s tot "<< c_length<< crys.Rcurv*sin(crys.cry_bend)<<endl;
1081	if ( xp >= xp_tangent) { ////////////////////////// IF NUMERO 13//////////////////
1082
1083	//if it hits the crystal, calculate in which point and apply the
1084	//transformation and drift to that point
1085	a_eq = (1 + xp * xp);
1086	b_eq = 2 * xp * (x - crys.Rcurv);
1087	c_eq = -2 * x * crys.Rcurv + x * x;
1088	Delta = b_eq * b_eq - 4 * (a_eq * c_eq);
1089	s_int = (-b_eq - sqrt(Delta)) / (2 * a_eq);
1090	//cout<< " s int "<<s_int <<endl;
1091	if (s_int < crys.Rcurv * sin(crys.cry_bend)) { ////////////////////////// fin IF NUMERO 14//////////////////
1092	// transform to a new ref system:shift and rotate
1093
1094	x_int = xp * s_int + x;
1095	xp_int = xp;
1096	z = z + zp * s_int;
1097	x = 0;
1098	s = 0;
1099
1100
1101	// tilt_int=2*X_int/S_int
1102	tilt_int = s_int / crys.Rcurv;
1103	xp = xp - tilt_int;
1104
1105	//---------------------------------------------------------------------------------------
1106	// Fortran Cyrstal routine for the protons
1107	//
1108	//---------------------------------------------------------------------------------------
1109
1110	part.x = x;
1111	part.xp = xp;
1112	part.y = z;
1113	part.yp = zp;
1114	part.PC = p;
1115	crys.Cry_length = crys.Cry_length - (tilt_int * crys.Rcurv);
1116
1117
1118
1119	/*
1120
1121	//Fortran crystal routine in SixTrack
1122	// Modified by Jianfeng Zhang @ LAL, 30/04/2014
1123	CRYST_(mat-7,x,xp,z,zp,p,cry_length);
1124	*/
1125
1126
1127	//------------------------------------------------------------------------------
1128	// The proton crystal routine in Fortran version
1129	//
1130	//------------------------------------------------------------------------------
1131
1132	//Fortran crystal routine in SixTrack
1133	// Modified by Jianfeng Zhang @ LAL, 30/04/2014
1134	double crysCry_length=crys.Cry_length;
1135	double crysRcurv=crys.Rcurv;
1136	double crysC_xmax=crys.C_xmax;
1137	double crysC_ymax=crys.C_ymax;
1138	double crysAlayer=crys.Alayer;
1139	int crysC_orient=crys.C_orient;
1140	double crysmiscut = crys.miscut;
1141	int IS =0;
1142	double length =0.0;
1143	IS = mat -7 ;
1144	length = crys.Cry_length;
1145	double eUmIS = crys.eUm[crys.IS];
1146	double AIIS = crys.AI[crys.IS];
1147
1148
1149	double crysparaDes[4];
1150	crysparaDes[0]=crys.DES[0],crysparaDes[1]=crys.DES[1],crysparaDes[2]=crys.DES[2],crysparaDes[3]=crys.DES[3];
1151	double crysparaDlyi[4];
1152	crysparaDlyi[0]=crys.DLYI[0],crysparaDlyi[1]=crys.DLYI[1],crysparaDlyi[2]=crys.DLYI[2],crysparaDlyi[3]=crys.DLYI[3];
1153	double crysparaDlri[4];
1154	crysparaDlri[0]=crys.DLRI[0],crysparaDlri[1]=crys.DLRI[1],crysparaDlri[2]=crys.DLRI[2],crysparaDlri[3]=crys.DLRI[3];
1155
1156	double crysparaeUmIS[4];
1157	crysparaeUmIS[0]=crys.eUm[0], crysparaeUmIS[1]=crys.eUm[1], crysparaeUmIS[2]=crys.eUm[2], crysparaeUmIS[3]=crys.eUm[3];
1158	double crysparaAIIS[4];
1159	crysparaAIIS[0]=crys.AI[0],crysparaAIIS[1]=crys.AI[1],crysparaAIIS[2]=crys.AI[2],crysparaAIIS[3]=crys.AI[3];
1160
1161
1162
1163
1164
1165	// cout << "before enter crystal routine...."<<endl;
1166
1167	cryst_(&length,&layerflag,Proc,&crys.L_chan,&crys.tchan,&crys.Alayer,&crys.ymin,&crys.ymax,&crys.Rcurv,
1168	&crys.s_length,&crys.IS,&crys.NAM,
1169	&crys.xpcrit0,&crys.xpcrit,&crys.Rcrit,&crys.ratio,&crys.C_orient,&crys.miscut,&crys.Cry_length,
1170	&crys.ZN,&crys.C_xmax,&crys.C_ymax,
1171	&part.x,&part.xp,&part.y,&part.yp,&part.PC,&part.s,&part.Chann,&part.xp_rel,&part.Ang_rms,
1172	& part.Ang_avr,&part.Vcapt,&part.Dechan,crysparaDes,crysparaDlri,crysparaDlyi,crysparaeUmIS,crysparaAIIS);
1173
1174	//convert back the Fortran variables to C++ variables
1175	//mat = IS +7;
1176	crys.Cry_length = length;
1177	crys.IS = crys.IS - 1;
1178
1179
1180
1181	/*
1182	crys.Rcurv = crysRcurv;
1183	crys.C_xmax = crysC_xmax;
1184	crys.C_ymax = crysC_ymax;
1185	crys.Alayer = crysAlayer;
1186	crys.C_orient = crysC_orient;
1187	crys.miscut = crysmiscut;
1188	*/
1189
1190	/*
1191	bases(sim);
1192
1193	//cout<<"part.xp"<<part.xp<<endl;
1194
1195	bool crossing(cross(sim));
1196	//cout<<"Crossing !!!!"<<endl;
1197	bool layering(layer(sim));
1198
1199	if (layering == true) {
1200	//cout<<" Layering =true " <<endl;
1201	parameters(sim);
1202	if ((abs(part.xp_rel) < crys.xpcrit)) {
1203	channel(sim);
1204	//cout<< " CHANNELING " << endl;
1205	} else {
1206	reflection(sim);
1207	//cout<< " REFLECTION " << endl;
1208	}
1209
1210	}
1211	*/
1212	//---------------------------------------------------------------------------------------
1213
1214
1215
1216	//cout<<"part.xp"<<part.xp<<endl;
1217	//cout<<"----------------------2222222222222222222-----------------"<<endl;
1218
1219	s = crys.Rcurv * sin(crys.cry_bend - tilt_int);
1220	zlm = crys.Rcurv * sin(crys.cry_bend - tilt_int);
1221
1222
1223	if(layerflag == 0){
1224	part.namor = part.namor + 1; /////////////COUNT/////////////
1225	part.effet = 2;
1226	}
1227	if (strncmp(Proc,"out",3)==0)
1228	crossing = false;
1229	else
1230	crossing = true;
1231
1232	if (crossing == true) { ////////////////////////// IF NUMERO 15//////////////////
1233	x_rot = x_int;
1234	s_rot = s_int;
1235	xp_rot = xp_int;
1236	s_shift = s_rot * cos(-Cry_tilt) + x_rot * sin(-Cry_tilt);
1237	x_shift = -s_rot * sin(-Cry_tilt) + x_rot * cos(-Cry_tilt);
1238	xp_shift = xp_rot + Cry_tilt;
1239	if (Cry_tilt < 0) { ////////////////////////// IF NUMERO 16//////////////////
1240	s_impact = s_shift;
1241	x_in0 = x_shift + shift;
1242	xp_in0 = xp_shift;
1243	} ////////////////////////// fin IF NUMERO 16//////////////////
1244	else { ////////////////////////// IF NUMERO 17//////////////////
1245	x_in0 = x_shift;
1246	s_impact = s_shift;
1247	xp_in0 = xp_shift;
1248	} ////////////////////////// fin IF NUMERO 17//////////////////
1249
1250	nhit = nhit + 1;
1251	//lhit = 100000000*ie + ITURN;
1252	impact = x_in0;
1253	indiv = xp_in0;
1254	} ////////////////////////// fin IF NUMERO 15//////////////////
1255
1256
1257
1258	x_temp = x;
1259	s_temp = s;
1260	xp_temp = xp;
1261	s = s_temp * cos(-tilt_int) + x_temp * sin(-tilt_int);
1262	x = -s_temp * sin(-tilt_int) + x_temp * cos(-tilt_int);
1263	xp = xp_temp + tilt_int;
1264
1265	// 2nd: shift back the 2 axis
1266	x = x + x_int;
1267	s = s + s_int;
1268
1269	} ////////////////////////// fin IF NUMERO 14//////////////////
1270	// Finish if (s_int < crys.Rcurv * sin(crys.cry_bend))
1271	else { ///////////////////////////////////////////////PROBLEM : TOO MUCH OF PARTICULES WITHOUT ANY CHANGE OF ANGLE///////////THE PROB IS HERE WITH THOSES TWO "ELSE"/////////////////////////
1272
1273	s = crys.Rcurv * sin(crys.Cry_length / crys.Rcurv);
1274	x = x + s * xp;
1275	z = z + s * zp;
1276	++part.nout;
1277	//cout<< "the particule does not hit the crystal (Trop basse dans neg.)"<<endl;
1278	}// Finish else if (s_int < crys.Rcurv * sin(crys.cry_bend))
1279
1280	}else {
1281	s = crys.Rcurv * sin(crys.Cry_length / crys.Rcurv);
1282	x = x + s * xp;
1283	z = z + s * zp;
1284	++part.nout;
1285	//cout<< "the particule does not hit the crystal (negatif et plus petit angle qu xp_tang)"<<endl;
1286	}
1287	//} ////////////////////////// fin IF NUMERO 13//////////////////
1288	// Finish else if ( xp >= xp_tangent)
1289
1290	}////////////////////// 12
1291	// Finish particle not hit the crystal else (x < 0)
1292
1293
1294
1295
1296	//**************************************************************************************************************************
1297	//trasform back from the crystal to the collimator reference system
1298	// 1st: un-rotate the coordinates
1299
1300
1301	x_rot = x;
1302	s_rot = s;
1303	xp_rot = xp;
1304	//cout<< "debug - S Cry RF 2" , s_rot <<endl;
1305	//cout<<"debug - X Cry RF 2" , x_rot <<endl;
1306	//cout<<"debug - XP Cry RF 2" , xp_rot <<endl;
1307	s_shift = s_rot * cos(-Cry_tilt) + x_rot * sin(-Cry_tilt);
1308	x_shift = -s_rot * sin(-Cry_tilt) + x_rot * cos(-Cry_tilt);
1309	xp_shift = xp_rot + Cry_tilt;
1310	// 2nd: shift back the reference frame
1311	if (Cry_tilt < 0) { ////////////////////////// IF NUMERO 18//////////////////
1312	s = s_shift;
1313	x = x_shift + shift;
1314	xp = xp_shift;
1315	} ////////////////////////// fin IF NUMERO 18//////////////////
1316	else {
1317	x = x_shift;
1318	s = s_shift;
1319	xp = xp_shift;
1320	}
1321	// 3rd: shift to new S=Length position
1322	x = xp * (c_length - s) + x;
1323	z = zp * (c_length - s) + z;
1324	s = c_length;
1325
1326
1327	nabs = 0;
1328	part.effet = nabs;
1329	//cout<< "debug - S Coll RF 2" , s_rot <<endl;
1330	//cout<<"debug - X Coll RF 2" , x_rot <<endl;
1331	//cout<<"debug - XP Coll RF 2" , xp_rot <<endl;
1332
1333	//cout<<"X1_coll"<<x<<"Z1_coll"<<z<<"XP1_coll"<<xp<<"ZP1_coll"<<zp <<"s" <<s<<endl;
1334
1335
1336	if(strncmp(Proc,"out",3)==0){
1337	part.nout = part.nout + 1; /////////////COUNT/////////////
1338	part.effet = 1;
1339	}
1340	if(layerflag == 0){
1341	part.namor = part.namor + 1; /////////////COUNT/////////////
1342	part.effet = 2;
1343	}
1344	if(strncmp(Proc,"AM",2)==0){
1345	part.namor = part.namor + 1; /////////////COUNT/////////////
1346	part.effet = 2;
1347	}
1348	if(strncmp(Proc,"CH",2)==0){
1349	part.nchann = part.nchann + 1;
1350	part.effet = 3;
1351	}
1352	if(strncmp(Proc,"DC",2)==0){
1353	part.ndechann = part.ndechann + 1; /////////////COUNT/////////////
1354	part.effet = 4;;
1355	}
1356	if(strncmp(Proc,"VR",2)==0){
1357	part.nvrefl = part.nvrefl + 1; /////////////COUNT/////////////
1358	part.effet = 5;
1359	}
1360	if(strncmp(Proc,"VC",2)==0){
1361	part.nvcapt = part.nvcapt + 1; /////////////COUNT/////////////
1362	part.effet = 6;}
1363
1364	// ?????????????????????????????????????
1365	if(part.effet == 2) {
1366	part_abs = 0;
1367	}else {
1368	part_abs = 1;
1369	}
1370
1371
1372	//=========================================================================================================================
1373	// Transform back to particle coordinates with opening and offset
1374
1375	if (part_abs == 0) { ////////////////////////// IF NUMERO 19//////////////////
1376	//
1377	// Include collimator tilt
1378
1379	if (tiltangle > 0) {
1380	x = x + tiltangle * c_length;
1381	xp = xp + tiltangle;
1382	} else if (tiltangle < 0) {
1383	x = x + tiltangle * c_length;
1384	xp = xp + tiltangle;
1385
1386	x = x - sin(tiltangle) * c_length;
1387	}
1388
1389	// Transform back to particle coordinates with opening and offset
1390
1391	z00 = z;
1392	x00 = x + mirror * c_offset;
1393	x = x + c_aperture / 2 + mirror * c_offset;
1394
1395	//+ Now mirror at the horizontal axis for negative X offset
1396
1397	x = mirror * x;
1398	xp = mirror * xp;
1399
1400	// Last do rotation into collimator frame
1401
1402	part.x = x * cos((-1) * c_rotation) + z * sin((-1) * c_rotation);
1403	part.y = z * cos((-1) * c_rotation) - x * sin((-1) * c_rotation);
1404	part.xp = xp * cos((-1) * c_rotation) + zp * sin((-1) * c_rotation);
1405	part.yp = zp * cos((-1) * c_rotation) - xp * sin((-1) * c_rotation);
1406
1407
1408
1409	//**********************************************************************************
1410	//For output.... by Zhang @ CERN, 24/04/2014
1411	//p_in = (1 + dpop) * p0;
1412	//From Dianiele
1413	//p_in = p;
1414	//s_in = s_in + s;
1415	part.PC = p;
1416	part.s = part.s + s;
1417
1418	if (part.effet == 1) { ////////////////////////// IF NUMERO 21//////////////////
1419
1420	part.x = 99.99e-3;
1421	part.y = 99.99e-3;
1422	cout<< "Mean that the particules is lost"<<endl;
1423	} ////////////////////////// fin IF NUMERO 21//////////////////
1424	} ////////////////////////// fin IF NUMERO 19//////////////////
1425
1426
1427	// End of check for particles not being lost before (see @330)
1428
1429	//} ////////////////////////// fin IF NUMERO 12//////////////////
1430	// End of loop over all particles
1431
1432
1433
1434	} //collimator with length = 0 ////////////////////////// fin IF NUMERO 1//////////////////
1435
1436	// =================================================================================FIN DE RETOUR AU COORDONEE DE ICOSIM====================================================================
1437	/*
1438	part.xp=xp;
1439	part.yp=zp;
1440	part.y=z;
1441	part.x=x;
1442	part.PC=p;
1443	*/
1444	xp1 = part.xp;
1445
1446
1447	//if(xp1-xp0!=0){
1448	//sortie << xp1 << "," << xp0<<"," << xfin<<"," << yfin<<endl;
1449	//}
1450	//cout<<"ICIIIIIIIIIIIIIIIIIIIIIIIII"<<part.PC<<endl;
1451
1452
1453	AV_xp1 = AV_xp1 + part.xp; //average x angle after interaction
1454	AV_yp1 = AV_yp1 + part.yp; //average y angle after interaction
1455	RMSxp1 = RMSxp1 + part.xp * part.xp; //rms x angle after interaction
1456	RMSyp1 = RMSyp1 + part.xp * part.xp; //rms y angle after interaction
1457	av_pc1 = av_pc1 + part.PC;
1458
1459
1460	cout<<" X apres :"<<part.x<<" XP apres :"<<part.xp<<" Y apres :"<<part.y<<" YP apres :"<<part.yp<<" PC apres :"<<part.PC<<endl;
1461	//cout<<endl;
1462	//sortieIco.setf(ios::scientific);
1463	sortieIco << part.x << "," << part.y << "," << part.xp << "," << part.yp << "," << part.PC << endl;
1464	//}//end of for....
1465
1466	}
1467	} //End of while
1468
1469
1470	} else {
1471	cerr<< "We can not open the file !" << endl;
1472	}
1473	entree.close();
1474
1475
1476
1477	AV_xp0 = AV_xp0 / count;
1478	AV_yp0 = AV_yp0 / count;
1479
1480	RMSxp0 = sqrt(RMSxp0 / (count));
1481	RMSyp0 = sqrt(RMSyp0 / (count));
1482
1483	SIGxp0 = sqrt((RMSxp0 * RMSxp0) - (AV_xp0 * AV_xp0));
1484	SIGyp0 = sqrt((RMSyp0 * RMSyp0) - (AV_yp0 * AV_yp0));
1485
1486	AV_xp1 = AV_xp1 / count;
1487	AV_yp1 = AV_yp1 / count;
1488	RMSxp1 = sqrt(RMSxp1 / (count));
1489	RMSyp1 = sqrt(RMSyp1 / (count));
1490
1491	SIGxp1 = sqrt((RMSxp1 * RMSxp1) - (AV_xp1 * AV_xp1));
1492	SIGyp1 = sqrt((RMSyp1 * RMSyp1) - (AV_yp1 * AV_yp1));
1493	av_pc0 = av_pc0 / count;
1494	av_pc1 = av_pc1 / count;
1495	//--------------------------------------------------------------------------
1496	//-------------write a summary file-----------------------------------------
1497	//cout<<endl;
1498	//cout<<endl;
1499	//cout<<endl;
1500
1501	//cout<<"avg xp_in: "<< AV_xp0 <<"// avg xp_out: "<< AV_xp1<<endl;
1502	//cout<<"rms xp_in: "<< RMSxp0 <<"// rms xp_out: "<< RMSxp1<<endl;
1503	//cout<<"sigma xp_in: "<< SIGxp0 <<"// sigma xp_out: "<< SIGxp1<<endl;
1504	//cout<<"avg pc_in: "<< av_pc0 <<"// avg pc_out: "<< av_pc1<<endl;
1505
1506
1507
1508	sortieIco.close();
1509	// cout<<endl;
1510	// affiche();
1511	file_out(pas, outputpath);
1512
1513	//cout<<endl;
1514	//cout<<endl;
1515	//cout<<endl;
1516	}
1517
1518
1519
1520	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1521	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1522	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1523	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1524
1525	//WE WRITE SOME FUNCTIONS IMPORTANT FUNCTION NOW
1526
1527	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1528	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1529	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1530	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1531
1532
1533	void SimCrys::move_am(int IS, int NAM, double DZ, double DEI, double DLY, double DLR, double& xp, double& yp, double& PC)
1534	{
1535	//cout<<"The move_am function is use !!!!! "<<endl;
1536	double DYA(0);
1537	double Wp(0);
1538	if (NAM == 0) {
1539	return;
1540	}
1541	xp = xp * 1000;
1542	yp = yp * 1000;
1543	//cout << "xp avt: "<<xp<<" yp avt: "<<yp<<endl;
1544	//DEI ---> dE/dx (stoping energy)
1545	PC = PC - DEI * DZ; //energy lost beacause of ionization process[GeV]
1546
1547	// Coulomb Scattering
1548	//DYA ---> rms of coloumb scattering
1549	DYA = (14 / PC) * sqrt(DZ / DLR); //rms scattering (mrad)
1550	xp = xp + DYA * random_gauss_cut();
1551	yp = yp + DYA * random_gauss_cut();
1552
1553	//Elastic scattering
1554	//cout<<"Plusieurs possibilitees : "<<endl;
1555	if (random() <= (DZ / crys.DLAI[IS])) {
1556	//cout<<"Poss. 1 !!!! "<<endl;
1557	xp = xp + 196 * random_gauss_cut() / PC; //angle elast. scattering in R plane
1558	yp = yp + 196 * random_gauss_cut() / PC;
1559	}
1560
1561	//Diffraction interaction
1562	if (random() <= (DZ / (DLY * 6.143))) {
1563	//cout<<"Poss. 2 !!!! "<<endl;
1564	xp = xp + 257 * random_gauss_cut() / PC; // angle elast. scattering in R plane[mr]
1565	yp = yp + 257 * random_gauss_cut() / PC;
1566	Wp = random();
1567	PC = PC - 0.5 * pow(0.3 * PC, Wp); // m0*c = 1 GeV/c for particle
1568
1569	}
1570
1571	//Inelastic interaction
1572	if (random() <= DZ / DLY) {
1573	//cout<<" Absorbed !! "<< endl;
1574	}
1575	//cout<<"XP APRES : "<<xp<<" YP APRES : "<<yp<<endl;
1576	xp = xp / 1000;
1577	yp = yp / 1000;
1578	}
1579
1580
1581	double SimCrys::random()
1582	{
1583	double scale = RAND_MAX;
1584	double base = rand() / scale;
1585	double fine = rand() / scale;
1586	return base + fine / scale;
1587	}
1588
1589	double SimCrys::random_dist()
1590	{
1591	double scale = RAND_MAX;
1592	double base = rand() / scale;
1593	double fine = rand() / scale;
1594	return (-2.5 * 1e-4 + (3.5 * 1e-4 ) * (base + fine / scale));
1595	}
1596
1597	double SimCrys::random_gauss()
1598	{
1599	const double PI (4.0 * atan(1.0));
1600	double U(random());
1601	double V(random());
1602	return (sqrt(-2 * log(U)) * cos(2 * PI * V));
1603	}
1604
1605	double SimCrys::random_gauss_cut()
1606	{
1607	double resu;
1608	do {
1609	resu = random_gauss();
1610	} while (abs(resu) >= 1);
1611	return resu;
1612	}
1613
1614	double SimCrys::random_gauss_dist() //POUR faire varier la moyenne de la gaussienne il faut mettre en argu la sim& !!!!!
1615	{
1616	const double PI (4.0 * atan(1.0));
1617	double U(random());
1618	double V(random());
1619	//moy=random();
1620	return 1e-4 * (sqrt(-2 * log(U)) * cos(2 * PI * V)); //+(0.0002*moy); //POUR faire varier la moyenne de la gaussienne aussi !!!!!
1621	}
1622
1623	/////////////////////////////////*******************************************************+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1624	//POUR 409 :
1625
1626	double SimCrys::random_gauss_dist_409() //POUR faire varier la moyenne de la gaussienne il faut mettre en argu la sim& !!!!!
1627	{
1628	const double PI (4.0 * atan(1.0));
1629	double U(random());
1630	double V(random());
1631	//moy=random();
1632	return (8.939203e-06) * (sqrt(-2 * log(U)) * cos(2 * PI * V)) + (-7.183296e-08); //POUR faire varier la moyenne de la gaussienne aussi !!!!!
1633	}
1634
1635	///////////////////////////////////************************************************************++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1636
1637	void SimCrys::affiche()
1638	{
1639	cout << "Crystal Curved Length [m] :" << crys.Cry_length << endl;
1640	cout << "Crystal Curvature Radius [m] :" << crys.Rcurv << endl;
1641	cout << "Crystal Bending Angle [urad] :" << crys.cry_bend << " , " << crys.Cry_length / crys.Rcurv << endl;
1642	cout << "Critical Radius :" << crys.Rcrit << endl;
1643	cout << "Ratio :" << crys.ratio << endl;
1644	cout << "Critical Angle for straight Crystal :" << crys.xpcrit0 << endl;
1645	cout << "Critical Angle for curved Crystal :" << crys.xpcrit << endl;
1646	cout << "Angle average reflexion : " << part.Ang_avr << endl;
1647	cout << "Full Channeling : " << crys.Cry_length / crys.Rcurv << endl;
1648	cout << endl;
1649	cout << endl;
1650	cout << "N.AMORPH :" << part.namor << endl;
1651	cout << "N.DECHANNELING:" << part.ndechann << endl;
1652	cout << "N.CHANNELING :" << part.nchann << endl;
1653	cout << "N.OUT :" << part.nout << endl;
1654	cout << "N.VREFLECTION :" << part.nvrefl << endl;
1655	cout << "N.VCAPTURE :" << part.nvcapt << endl;
1656
1657	//cout << "Si la somme ne donne pas le nbre de part., c est normale !! En effet, il est possible que plusieurs evenement se produisent!!!";
1658	}
1659
1660	void SimCrys::file_out(const int& pas, string outputpath)
1661	{
1662
1663	ofstream out;
1664	string file;
1665	file = outputpath + "/crystal_output.out";
1666
1667	out.open(file.c_str(), ios::out \| ios::app);
1668
1669	if (out.fail()) {
1670	cerr << "Warning: problem writing in the file " << file << "!" << endl;
1671	} else {
1672
1673	out << endl;
1674	out << "Crystal Curved Length [m] :" << crys.Cry_length << endl;
1675	out << "Crystal Curvature Radius [m] :" << crys.Rcurv << endl;
1676	out << "Crystal Bending Angle [urad] :" << crys.cry_bend << " , " << crys.Cry_length / crys.Rcurv << endl;
1677	out << "Critical Radius :" << crys.Rcrit << endl;
1678	out << "Ratio :" << crys.ratio << endl;
1679	out << "Critical Angle for straight Crystal :" << crys.xpcrit0 << endl;
1680	out << "Critical Angle for curved Crystal :" << crys.xpcrit << endl;
1681	out << "Angle average reflexion : " << part.Ang_avr << endl;
1682	out << "Full Channeling : " << crys.Cry_length / crys.Rcurv << endl;
1683	out << endl;
1684	out << endl;
1685	out << "N.AMORPH :" << part.namor << endl;
1686	out << "N.DECHANNELING:" << part.ndechann << endl;
1687	out << "N.CHANNELING :" << part.nchann << endl;
1688	out << "N.OUT :" << part.nout << endl;
1689	out << "N.VREFLECTION :" << part.nvrefl << endl;
1690	out << "N.VCAPTURE :" << part.nvcapt << endl << endl << endl;
1691
1692	out.close();
1693	}
1694
1695	}
1696

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

Download in other formats: