Context Navigation

simcrys.cc @ 17

Last change on this file since 17 was 17, checked in by zhangj, 10 years ago
Added comments...
File size: 64.6 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, long int& nhit, 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
663
664	//----------------------------------END OF INITATING REFERNENTIAL PARAMETERS-------------------------------------------------------------------------
665	c_length = crys.Cry_length;
666
667	if (crys.IS == 0) {
668	mat = 8;
669	} else if (crys.IS == 1) {
670	mat = 9;
671	} else if (crys.IS == 2) {
672	mat = 10;
673	} else if (crys.IS == 3) {
674	mat = 11;
675	} else {
676	cout << "Material of the Crystal not fund !!!" << endl;
677	}
678
679
680	if (c_length > 0 ) {
681	p0 = part.PC;
682	c_tilt0[0] = c_tilt[0];
683	c_tilt0[1] = c_tilt[1];
684	tiltangle = c_tilt0[0];
685	}
686
687
688	//new..........................
689	// call scatin(p0) ??????????????
690
691	// scatin_(p0);
692
693	/* EVENT LOOP, initial distribution is here a flat distribution with
694	* xmin=x-, xmax=x+, etc. from the input file
695	*/
696
697
698	//initialize particle variables
699	part.n_part = 0;
700	part.nabsorbed = 0;
701	part.namor = 0;
702	part.nchann = 0;
703	part.ndechann = 0;
704	part.nout =0;
705	part.nvcapt = 0;
706	part.nvrefl = 0;
707	part.nhit = 0;
708	nhit = 0;
709	double fracab = 0;
710	// int n_chan = 0; // valentina :initialize to zero the counters for crystal effects
711	// int n_VR = 0;
712	// int n_amorphous = 0;
713
714
715	double impact = 0;
716	double indiv = 0;
717	double lint = 0;
718
719
720	// c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
721	//c WRITE(,) ITURN,ICOLL
722	/* do j = 1, nev
723	//c
724	impact(j) = -1.0;
725	lint(j) = -1.0;
726	indiv(j) = -1.0;
727
728
729	idx_proc = name(j)-100*samplenumber;
730
731	if (ITURN == 1){ // !daniele
732	bool_proc_old(idx_proc)=-1;} // !daniele
733	// elseif (bool_proc(idx_proc) != -1){
734	else{
735	bool_proc_old(idx_proc)=bool_proc(idx_proc); // !daniele
736	} // !daniele
737
738	PROC= "out"; // !the default process is 'out' !daniele
739	bool_proc(idx_proc)=-1; // !daniele
740	*/
741
742	int nabs = 0;
743
744	//ofstream sortie("results.csv");
745	//sortie << "xp1" <<","<< "xp0"<<","<< "x"<<","<< "y"<<endl;
746	ofstream sortieIco;
747	string name1;
748
749	//file with the particle coordinate after the crystal
750	name1 = outputpath + "/ascii_after.csv";
751
752	sortieIco.open(name1.c_str());
753
754	string rest;
755
756	ifstream entree;
757	string name2;
758	name2 = outputpath + "/ascii_before.csv";
759
760	entree.open(name2.c_str());
761	if (entree) {
762	while (!entree.eof()) {
763	count = count + 1;
764
765	//particle hit the crystal or not
766	bool hitflag = false;
767
768
769	getline(entree, rest, ',');
770	part.x = atof(rest.c_str());
771	//cout<<" X " << part.x<<endl;
772
773	getline(entree, rest, ',');
774	part.y = atof(rest.c_str());
775	//cout<<" Y " << part.y<<endl;
776
777	getline(entree, rest, ',');
778	part.xp = atof(rest.c_str());
779	//cout<<" XP " << part.xp<<endl;
780
781	getline(entree, rest, ',');
782	part.yp = atof(rest.c_str());
783	//cout<<" YP " << part.yp<<endl;
784
785	getline(entree, rest);
786	part.PC = atof(rest.c_str());
787
788	if (part.PC != 0) { // POUR OTER LA DERNIERE PART QUI VAUT 0 POUR TOUT
789
790	av_pc0 = av_pc0 + part.PC;
791
792	//part.x=random_gauss()*0.00115;
793	//part.xp=random_gauss()sqrt(51e-6);
794	//cout<<"X et Y avt :"<<part.x<<" , "<<part.y<<endl;
795
796	// IL FAUT PASSER DES CM EN M!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
797	//bdelta=tan(part.xp)*10;
798	part.x = (part.x); //-0.0045;
799
800
801	//bdelta=tan(part.yp)*10;
802	part.y = (part.y); //-0.0045;
803
804	//cout<<"X et Y apres :"<<part.x<<" , "<<part.y<<endl;
805	xp0 = part.xp;
806
807	//part.y=random_gauss()*sqrt(0.00095);
808	//part.yp=5*1e-6;
809
810	//part.PC=120.976;
811	//cout<<" X :"<<part.x<<" XP :"<<part.xp<<" Y :"<<part.y<<" YP :"<<part.yp<<" PC :"<<part.PC<<endl;
812	AV_xp0 = AV_xp0 + part.xp; //average x angle before impact
813	AV_yp0 = AV_yp0 + part.yp; //average y angle before impact
814	RMSxp0 = RMSxp0 + part.xp * part.xp; //rms x angle before impact
815	RMSyp0 = RMSyp0 + part.yp * part.yp; //rms y angle before impact
816
817	// NOW WE NEED TO CHANGE THE REFERENTIAL!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
818	//change the beam frame to the crystal frame !!!!!!!!!
819	// since the interaction between the particle and the crystal is observed in the
820	// crystal reference system
821
822	//initial coordinates in the beam frame
823	s = 0;
824	x = part.x;
825	xp = part.xp;
826	z = part.y;
827	zp = part.yp;
828	p = part.PC;
829
830	x_temp = 0;
831	x_int = 0;
832	x_rot = 0;
833	x_shift = 0;
834	s_temp = 0;
835	s_int = 0;
836	s_rot = 0;
837	s_shift = 0;
838	xp_int = 0;
839	xp_temp = 0;
840	xp_rot = 0;
841	xp_shift = 0;
842	shift = 0;
843	tilt_int = 0;
844	dpop = (p - p0) / p0; //momentum spread
845
846
847	/*---------------DANIELE-------------------------
848	! corrected position for variable association for FirstImpact.dat
849	! also for a vertical crystal case.
850	! Only x_in0 (i.e. b) have to be assigned after the change of reference frame
851	!-----------------------------------------------*/
852
853	s_in0 = s_in; //daniele
854	// x_in0(j) = x //!daniele
855	xp_in0 = xp; //daniele
856	y_in0 = z; //daniele
857	yp_in0 = zp; //daniele
858	p_in0 = p; //daniele
859
860	/----------------DANIELE------------------/
861
862
863	//**********************************************************************************
864
865	// ------------------------------transform particle coordinates to get into collimator coordinate system -----------------------------------------
866
867
868	// first check whether particle was lost before
869	totcount = totcount + 1;
870	if ((x < 99.0 * 1e-3) && (z < 99.0 * 1e-3)) { ////////////////////////// IF NUMERO 1//////////////////
871
872	// first do rotation into collimator frame
873
874	x = part.x * cos(c_rotation) + sin(c_rotation) * part.y;
875	z = part.y * cos(c_rotation) - sin(c_rotation) * part.x;
876	xp = part.xp * cos(c_rotation) + sin(c_rotation) * part.yp;
877	zp = part.yp * cos(c_rotation) - sin(c_rotation) * part.xp;
878
879	//cout<<"PARTICULE NON PERDU>>><<"<<endl;
880
881	//To define the correct value of the mirror for the install locations
882	// of the crystal in SPS & LHC, inner side or outside of the vacuum chamber????
883	// by Zhang @ CERN, 23/04/2014.....
884	//mirror = 1;
885
886
887	x = mirror * x;
888	xp = mirror * xp;
889
890	//Shift with opening and offset
891
892	x = x - c_aperture / 2 - c_offset;
893
894	// Include collimator tilt
895
896	if (tiltangle > 0) { ////////////////////////// IF NUMERO 2//////////////////
897	xp = xp - tiltangle;
898	} //////////////////////////FIN IF NUMERO 2//////////////////
899	else if (tiltangle < 0) { ////////////////////////// IF NUMERO 3//////////////////
900	x = x + sin(tiltangle) * c_length;
901	xp = xp - tiltangle;
902	} //////////////////////////FIN IF NUMERO 3//////////////////
903
904	if (Cry_tilt < 0) { ////////////////////////// IF NUMERO 4//////////////////
905
906	s_shift = s;
907	shift = crys.Rcurv * (1 - cos(Cry_tilt));
908	if (Cry_tilt < (-crys.cry_bend) ) { ////////////////////////// IF NUMERO 5//////////////////
909	shift = ( crys.Rcurv * ( cos ( - Cry_tilt) - cos( crys.cry_bend - Cry_tilt) ) );
910	} //////////////////////////FIN IF NUMERO 5//////////////////
911	x_shift = x - shift;
912	} //////////////////////////fin IF NUMERO 4//////////////////
913	else { //////////////////////////IF NUMERO 6//////////////////
914	//cout<< "CRY-TILT POSSITIVE ++ or nul"<<endl;
915	s_shift = s;
916	x_shift = x;
917	} //////////////////////////fin IF NUMERO 6//////////////////
918	// cout<<"debug - S shift" << s_shift<<endl;
919	//cout<<"debug - X shift" << x_shift <<endl;
920	//
921	// 2nd transformation: rotation
922	s_rot = x_shift * sin(Cry_tilt) + s_shift * cos(Cry_tilt);
923	x_rot = x_shift * cos(Cry_tilt) - s_shift * sin(Cry_tilt);
924	xp_rot = xp - Cry_tilt;
925	//cout<< "debug - S rot" << s_rot <<endl;
926	//cout<<"debug - X rot" << x_rot <<endl;
927	//cout<<"debug - XP rot" << xp_rot<<endl;
928
929	// 3rd transformation: drift to the new coordinate s=0
930	xp = xp_rot;
931	x = x_rot - xp_rot * s_rot;
932	z = z - zp * s_rot;
933	s = 0;
934	//cout<<"debug - S cryRF" << s_rot <<endl;
935	//cout<<"debug - X cryRF" << x_rot <<endl;
936	//cout<<"debug - XP cryRF" << xp_rot <<endl;
937
938
939
940	//********************************************************************************************************************
941	//----------------------------------NOW CHECK IF THE PARTICLE HIT the crystal---------------------------------------------------------------
942	//cout<<" X " << x<<endl;
943	//Particle hit the crystal
944	if (x >= 0) { ////////////////////////// IF NUMERO 7//////////////////
945	s_impact = s_in0; //!(for the first impact)
946
947	if(0){
948	cout<<"hit the cry entrance face"<<endl;
949	cout<<"impact at s = "<< s_impact<< ", x = "<<x_in0<<endl;
950	cout<<"with angle xp = "<<xp<<endl;
951	cout<<"s before: "<< s<<endl;
952	}
953
954	part.x = x;
955	part.xp = xp;
956	part.y = z;
957	part.yp = zp;
958	part.PC = p;
959
960	//------------------------------------------------------------------------------
961	// The proton crystal routine in Fortran version
962	//
963	//------------------------------------------------------------------------------
964
965	//Fortran crystal routine in SixTrack
966	// Modified by Jianfeng Zhang @ LAL, 30/04/2014
967	double crysCry_length=crys.Cry_length;
968	double crysRcurv=crys.Rcurv;
969	double crysC_xmax=crys.C_xmax;
970	double crysC_ymax=crys.C_ymax;
971	double crysAlayer=crys.Alayer;
972	int crysC_orient=crys.C_orient;
973	double crysmiscut = crys.miscut;
974	int IS =0;
975	double length =0.0;
976	IS = mat -7 ;
977	length = crys.Cry_length;
978	double eUmIS = crys.eUm[crys.IS];
979	double AIIS = crys.AI[crys.IS];
980
981
982
983	double crysparaDes[4];
984	crysparaDes[0]=crys.DES[0],crysparaDes[1]=crys.DES[1],crysparaDes[2]=crys.DES[2],crysparaDes[3]=crys.DES[3];
985	double crysparaDlyi[4];
986	crysparaDlyi[0]=crys.DLYI[0],crysparaDlyi[1]=crys.DLYI[1],crysparaDlyi[2]=crys.DLYI[2],crysparaDlyi[3]=crys.DLYI[3];
987	double crysparaDlri[4];
988	crysparaDlri[0]=crys.DLRI[0],crysparaDlri[1]=crys.DLRI[1],crysparaDlri[2]=crys.DLRI[2],crysparaDlri[3]=crys.DLRI[3];
989
990	double crysparaeUmIS[4];
991	crysparaeUmIS[0]=crys.eUm[0], crysparaeUmIS[1]=crys.eUm[1], crysparaeUmIS[2]=crys.eUm[2], crysparaeUmIS[3]=crys.eUm[3];
992	double crysparaAIIS[4];
993	crysparaAIIS[0]=crys.AI[0],crysparaAIIS[1]=crys.AI[1],crysparaAIIS[2]=crys.AI[2],crysparaAIIS[3]=crys.AI[3];
994
995	// cout << "before enter crystal routine...."<<endl;
996
997	// cryst_(&IS,&x,&xp,&z,&zp,&p,&length,&crysCry_length,&crysRcurv,&crysC_xmax,&crysC_ymax,&crysAlayer,&crysC_orient,&crysmiscut, &Proc,&layerflag, &part.Chann);
998
999	// cryst_(&length,&layerflag, crypara,partpara,crysparaDes,crysparaDlri,crysparaDlyi,crysparaeUmIS,crysparaAIIS,Proc);
1000
1001	cryst_(&length,&layerflag,Proc,&crys.L_chan,&crys.tchan,&crys.Alayer,&crys.ymin,&crys.ymax,&crys.Rcurv,
1002	&crys.s_length,&crys.IS,&crys.NAM,
1003	&crys.xpcrit0,&crys.xpcrit,&crys.Rcrit,&crys.ratio,&crys.C_orient,&crys.miscut,&crys.Cry_length,
1004	&crys.ZN,&crys.C_xmax,&crys.C_ymax,
1005	&part.x,&part.xp,&part.y,&part.yp,&part.PC,&part.s,&part.Chann,&part.xp_rel,&part.Ang_rms,
1006	& part.Ang_avr,&part.Vcapt,&part.Dechan,crysparaDes,crysparaDlri,crysparaDlyi,crysparaeUmIS,crysparaAIIS);
1007
1008
1009
1010	//convert back the Fortran variables to C++ variables
1011	//mat = IS +7;
1012
1013	// cout << "IS in Fortran = ..." << crys.IS << endl;
1014
1015	crys.IS = crys.IS - 1;
1016
1017	// cout << "IS in C++ = ..." << crys.IS << endl;
1018
1019	crys.Cry_length = length;
1020
1021	x=part.x;
1022	xp=part.xp;
1023	z=part.y;
1024	zp=part.yp;
1025	p=part.PC;
1026	s=part.s;
1027
1028
1029	/*
1030	crys.Rcurv = crysRcurv;
1031	crys.C_xmax = crysC_xmax;
1032	crys.C_ymax = crysC_ymax;
1033	crys.Alayer = crysAlayer;
1034	crys.C_orient = crysC_orient;
1035	crys.miscut = crysmiscut;
1036	*/
1037
1038	// cout<< "ICI AHHHHHHHHHHHHHHHHHHHHHHHHHHHH"<<endl;
1039	// //cout<<" xp rel and xp crit " << part.xp_rel<< " " << crys.xpcrit<<endl;
1040	// bases(sim);
1041	// //cout << part.xp<< " XP "<<endl;
1042	// bool crossing(cross(sim));
1043	// //cout<<part.xp-xp0<<endl;
1044	// //cout<<" xp rel and xp crit " << part.xp_rel<< " " << crys.xpcrit<<endl;
1045	// bool layering;
1046	// if (crossing != 0) { ////////////////////////// IF NUMERO 8//////////////////
1047	// layering = layer(sim);
1048	// //cout<<"Crossing !!!!"<<endl;
1049	//
1050	// if (layering == true) { ////////////////////////// IF NUMERO 9//////////////////
1051	// parameters(sim);
1052	// //cout<<" Layering =true " <<endl;
1053	//
1054	// if ((abs(part.xp_rel) < crys.xpcrit)) { ////////////////////////// IF NUMERO 10//////////////////
1055	// channel(sim);
1056	// //cout<< " CHANNELING " << endl;
1057	// } ////////////////////////// fin IF NUMERO 10//////////////////
1058	// else { ////////////////////////// IF NUMERO 11//////////////////
1059	// reflection(sim);
1060	// //cout<<" REFLECTION "<<endl;
1061	// } ////////////////////////// fin IF NUMERO 11//////////////////
1062	//
1063	//
1064	// } ////////////////////////// fin IF NUMERO 9//////////////////
1065	// } ////////////////////////// fin IF NUMERO 8//////////////////
1066
1067	//--------------------------------------------------------------------------------------------------
1068
1069
1070	// xp = part.xp;
1071
1072	//}
1073
1074	s = crys.Rcurv * sin(crys.cry_bend);
1075	zlm = crys.Rcurv * sin(crys.cry_bend);
1076	//cout<<"process:"<<PROC<<endl;
1077	//cout<<"s after"<< s<<endl;
1078	//cout<<"part.xp"<<part.xp<<endl;
1079
1080	//cout<<"----------------------1111111111111111111111111-----------------"<<endl;
1081
1082	if(layerflag == 0){
1083	part.namor = part.namor + 1; /////////////COUNT/////////////
1084	part.effet = 2;
1085	}
1086
1087	if (strncmp(Proc,"out",3)==0)
1088	crossing = false;
1089	else
1090	crossing = true;
1091
1092	if (crossing == true){
1093	hitflag = true;
1094	nhit = nhit + 1;
1095	part.nhit = nhit;
1096	// lhit = 100000000*ie + ITURN;
1097	impact = x_in0;
1098	indiv = xp_in0;
1099	}
1100	}
1101	// Particle not hit the crystal
1102	////////////////////////// fin IF NUMERO 7//////////////////
1103	else { ////////////////////////// IF NUMERO 12//////////////////
1104
1105	//cout<<"OU LA ABBBBBBBBBBBB" <<endl;
1106	xp_tangent = sqrt((-2 * x * crys.Rcurv + x * x) / (crys.Rcurv * crys.Rcurv));
1107	//cout<<"RCURV "<<crys.Rcurv<<endl;
1108	//cout<<"tangent"<<xp_tangent<<"angle"<<xp<<endl;
1109
1110	//cout<<"s tan " << crys.Rcurv*sin(xp_tangent)<<endl;
1111	//cout<< " s tot "<< c_length<< crys.Rcurv*sin(crys.cry_bend)<<endl;
1112	if ( xp >= xp_tangent) { ////////////////////////// IF NUMERO 13//////////////////
1113
1114	//if it hits the crystal, calculate in which point and apply the
1115	//transformation and drift to that point
1116	a_eq = (1 + xp * xp);
1117	b_eq = 2 * xp * (x - crys.Rcurv);
1118	c_eq = -2 * x * crys.Rcurv + x * x;
1119	Delta = b_eq * b_eq - 4 * (a_eq * c_eq);
1120	s_int = (-b_eq - sqrt(Delta)) / (2 * a_eq);
1121	//cout<< " s int "<<s_int <<endl;
1122	if (s_int < crys.Rcurv * sin(crys.cry_bend)) { ////////////////////////// fin IF NUMERO 14//////////////////
1123	// transform to a new ref system:shift and rotate
1124
1125	x_int = xp * s_int + x;
1126	xp_int = xp;
1127	z = z + zp * s_int;
1128	x = 0;
1129	s = 0;
1130
1131
1132	// tilt_int=2*X_int/S_int
1133	tilt_int = s_int / crys.Rcurv;
1134	xp = xp - tilt_int;
1135
1136	//---------------------------------------------------------------------------------------
1137	// Fortran Cyrstal routine for the protons
1138	//
1139	//---------------------------------------------------------------------------------------
1140
1141	part.x = x;
1142	part.xp = xp;
1143	part.y = z;
1144	part.yp = zp;
1145	part.PC = p;
1146	crys.Cry_length = crys.Cry_length - (tilt_int * crys.Rcurv);
1147
1148
1149
1150	/*
1151
1152	//Fortran crystal routine in SixTrack
1153	// Modified by Jianfeng Zhang @ LAL, 30/04/2014
1154	CRYST_(mat-7,x,xp,z,zp,p,cry_length);
1155	*/
1156
1157
1158	//------------------------------------------------------------------------------
1159	// The proton crystal routine in Fortran version
1160	//
1161	//------------------------------------------------------------------------------
1162
1163	//Fortran crystal routine in SixTrack
1164	// Modified by Jianfeng Zhang @ LAL, 30/04/2014
1165	double crysCry_length=crys.Cry_length;
1166	double crysRcurv=crys.Rcurv;
1167	double crysC_xmax=crys.C_xmax;
1168	double crysC_ymax=crys.C_ymax;
1169	double crysAlayer=crys.Alayer;
1170	int crysC_orient=crys.C_orient;
1171	double crysmiscut = crys.miscut;
1172	int IS =0;
1173	double length =0.0;
1174	IS = mat -7 ;
1175	length = crys.Cry_length;
1176	double eUmIS = crys.eUm[crys.IS];
1177	double AIIS = crys.AI[crys.IS];
1178
1179
1180	double crysparaDes[4];
1181	crysparaDes[0]=crys.DES[0],crysparaDes[1]=crys.DES[1],crysparaDes[2]=crys.DES[2],crysparaDes[3]=crys.DES[3];
1182	double crysparaDlyi[4];
1183	crysparaDlyi[0]=crys.DLYI[0],crysparaDlyi[1]=crys.DLYI[1],crysparaDlyi[2]=crys.DLYI[2],crysparaDlyi[3]=crys.DLYI[3];
1184	double crysparaDlri[4];
1185	crysparaDlri[0]=crys.DLRI[0],crysparaDlri[1]=crys.DLRI[1],crysparaDlri[2]=crys.DLRI[2],crysparaDlri[3]=crys.DLRI[3];
1186
1187	double crysparaeUmIS[4];
1188	crysparaeUmIS[0]=crys.eUm[0], crysparaeUmIS[1]=crys.eUm[1], crysparaeUmIS[2]=crys.eUm[2], crysparaeUmIS[3]=crys.eUm[3];
1189	double crysparaAIIS[4];
1190	crysparaAIIS[0]=crys.AI[0],crysparaAIIS[1]=crys.AI[1],crysparaAIIS[2]=crys.AI[2],crysparaAIIS[3]=crys.AI[3];
1191
1192
1193
1194
1195
1196	// cout << "before enter crystal routine...."<<endl;
1197
1198	cryst_(&length,&layerflag,Proc,&crys.L_chan,&crys.tchan,&crys.Alayer,&crys.ymin,&crys.ymax,&crys.Rcurv,
1199	&crys.s_length,&crys.IS,&crys.NAM,
1200	&crys.xpcrit0,&crys.xpcrit,&crys.Rcrit,&crys.ratio,&crys.C_orient,&crys.miscut,&crys.Cry_length,
1201	&crys.ZN,&crys.C_xmax,&crys.C_ymax,
1202	&part.x,&part.xp,&part.y,&part.yp,&part.PC,&part.s,&part.Chann,&part.xp_rel,&part.Ang_rms,
1203	& part.Ang_avr,&part.Vcapt,&part.Dechan,crysparaDes,crysparaDlri,crysparaDlyi,crysparaeUmIS,crysparaAIIS);
1204
1205	//convert back the Fortran variables to C++ variables
1206	//mat = IS +7;
1207	crys.Cry_length = length;
1208	crys.IS = crys.IS - 1;
1209
1210
1211
1212	x=part.x;
1213	xp=part.xp;
1214	z=part.y;
1215	zp=part.yp;
1216	p=part.PC;
1217	s=part.s;
1218
1219	/*
1220	crys.Rcurv = crysRcurv;
1221	crys.C_xmax = crysC_xmax;
1222	crys.C_ymax = crysC_ymax;
1223	crys.Alayer = crysAlayer;
1224	crys.C_orient = crysC_orient;
1225	crys.miscut = crysmiscut;
1226	*/
1227
1228	/*
1229	bases(sim);
1230
1231	//cout<<"part.xp"<<part.xp<<endl;
1232
1233	bool crossing(cross(sim));
1234	//cout<<"Crossing !!!!"<<endl;
1235	bool layering(layer(sim));
1236
1237	if (layering == true) {
1238	//cout<<" Layering =true " <<endl;
1239	parameters(sim);
1240	if ((abs(part.xp_rel) < crys.xpcrit)) {
1241	channel(sim);
1242	//cout<< " CHANNELING " << endl;
1243	} else {
1244	reflection(sim);
1245	//cout<< " REFLECTION " << endl;
1246	}
1247
1248	}
1249	*/
1250	//---------------------------------------------------------------------------------------
1251
1252
1253
1254	//cout<<"part.xp"<<part.xp<<endl;
1255	//cout<<"----------------------2222222222222222222-----------------"<<endl;
1256
1257	s = crys.Rcurv * sin(crys.cry_bend - tilt_int);
1258	zlm = crys.Rcurv * sin(crys.cry_bend - tilt_int);
1259
1260
1261	if(layerflag == 0){
1262	part.namor = part.namor + 1; /////////////COUNT/////////////
1263	part.effet = 2;
1264	}
1265	if (strncmp(Proc,"out",3)==0)
1266	crossing = false;
1267	else
1268	crossing = true;
1269
1270	if (crossing == true) { ////////////////////////// IF NUMERO 15//////////////////
1271	x_rot = x_int;
1272	s_rot = s_int;
1273	xp_rot = xp_int;
1274	s_shift = s_rot * cos(-Cry_tilt) + x_rot * sin(-Cry_tilt);
1275	x_shift = -s_rot * sin(-Cry_tilt) + x_rot * cos(-Cry_tilt);
1276	xp_shift = xp_rot + Cry_tilt;
1277	if (Cry_tilt < 0) { ////////////////////////// IF NUMERO 16//////////////////
1278	s_impact = s_shift;
1279	x_in0 = x_shift + shift;
1280	xp_in0 = xp_shift;
1281	} ////////////////////////// fin IF NUMERO 16//////////////////
1282	else { ////////////////////////// IF NUMERO 17//////////////////
1283	x_in0 = x_shift;
1284	s_impact = s_shift;
1285	xp_in0 = xp_shift;
1286	} ////////////////////////// fin IF NUMERO 17//////////////////
1287	hitflag = true;
1288	nhit = nhit + 1;
1289	part.nhit = nhit;
1290	//lhit = 100000000*ie + ITURN;
1291	impact = x_in0;
1292	indiv = xp_in0;
1293	} ////////////////////////// fin IF NUMERO 15//////////////////
1294
1295
1296
1297	x_temp = x;
1298	s_temp = s;
1299	xp_temp = xp;
1300	s = s_temp * cos(-tilt_int) + x_temp * sin(-tilt_int);
1301	x = -s_temp * sin(-tilt_int) + x_temp * cos(-tilt_int);
1302	xp = xp_temp + tilt_int;
1303
1304	// 2nd: shift back the 2 axis
1305	x = x + x_int;
1306	s = s + s_int;
1307
1308	} ////////////////////////// fin IF NUMERO 14//////////////////
1309	// Finish if (s_int < crys.Rcurv * sin(crys.cry_bend))
1310	else { ///////////////////////////////////////////////PROBLEM : TOO MUCH OF PARTICULES WITHOUT ANY CHANGE OF ANGLE///////////THE PROB IS HERE WITH THOSES TWO "ELSE"/////////////////////////
1311
1312	s = crys.Rcurv * sin(crys.Cry_length / crys.Rcurv);
1313	x = x + s * xp;
1314	z = z + s * zp;
1315	++part.nout;
1316	//cout<< "the particule does not hit the crystal (Trop basse dans neg.)"<<endl;
1317	}// Finish else if (s_int < crys.Rcurv * sin(crys.cry_bend))
1318
1319	}else {
1320	s = crys.Rcurv * sin(crys.Cry_length / crys.Rcurv);
1321	x = x + s * xp;
1322	z = z + s * zp;
1323	++part.nout;
1324	//cout<< "the particule does not hit the crystal (negatif et plus petit angle qu xp_tang)"<<endl;
1325	}
1326	//} ////////////////////////// fin IF NUMERO 13//////////////////
1327	// Finish else if ( xp >= xp_tangent)
1328
1329	}////////////////////// 12
1330	// Finish particle not hit the crystal else (x < 0)
1331
1332
1333
1334
1335	//**************************************************************************************************************************
1336	//trasform back from the crystal to the collimator reference system
1337	// 1st: un-rotate the coordinates
1338
1339
1340	x_rot = x;
1341	s_rot = s;
1342	xp_rot = xp;
1343	//cout<< "debug - S Cry RF 2" , s_rot <<endl;
1344	//cout<<"debug - X Cry RF 2" , x_rot <<endl;
1345	//cout<<"debug - XP Cry RF 2" , xp_rot <<endl;
1346	s_shift = s_rot * cos(-Cry_tilt) + x_rot * sin(-Cry_tilt);
1347	x_shift = -s_rot * sin(-Cry_tilt) + x_rot * cos(-Cry_tilt);
1348	xp_shift = xp_rot + Cry_tilt;
1349	// 2nd: shift back the reference frame
1350	if (Cry_tilt < 0) { ////////////////////////// IF NUMERO 18//////////////////
1351	s = s_shift;
1352	x = x_shift + shift;
1353	xp = xp_shift;
1354	} ////////////////////////// fin IF NUMERO 18//////////////////
1355	else {
1356	x = x_shift;
1357	s = s_shift;
1358	xp = xp_shift;
1359	}
1360	// 3rd: shift to new S=Length position
1361	x = xp * (c_length - s) + x;
1362	z = zp * (c_length - s) + z;
1363	s = c_length;
1364
1365
1366	nabs = 0; //particle lost or not lost.....
1367	part.effet = nabs;
1368	//cout<< "debug - S Coll RF 2" , s_rot <<endl;
1369	//cout<<"debug - X Coll RF 2" , x_rot <<endl;
1370	//cout<<"debug - XP Coll RF 2" , xp_rot <<endl;
1371
1372	//cout<<"X1_coll"<<x<<"Z1_coll"<<z<<"XP1_coll"<<xp<<"ZP1_coll"<<zp <<"s" <<s<<endl;
1373
1374
1375	if(strncmp(Proc,"out",3)==0){
1376	part.nout = part.nout + 1; /////////////COUNT/////////////
1377	part.effet = 1;
1378	}
1379	if(layerflag == 0){
1380	part.namor = part.namor + 1; /////////////COUNT/////////////
1381	part.effet = 2;
1382	}
1383	if(strncmp(Proc,"AM",2)==0){
1384	part.namor = part.namor + 1; /////////////COUNT/////////////
1385	part.effet = 2;
1386	}
1387	if(strncmp(Proc,"CH",2)==0){
1388	part.nchann = part.nchann + 1;
1389	part.effet = 3;
1390	}
1391	if(strncmp(Proc,"DC",2)==0){
1392	part.ndechann = part.ndechann + 1; /////////////COUNT/////////////
1393	part.effet = 4;;
1394	}
1395	if(strncmp(Proc,"VR",2)==0){
1396	part.nvrefl = part.nvrefl + 1; /////////////COUNT/////////////
1397	part.effet = 5;
1398	}
1399	if(strncmp(Proc,"VC",2)==0){
1400	part.nvcapt = part.nvcapt + 1; /////////////COUNT/////////////
1401	part.effet = 6;}
1402	if(strncmp(Proc,"absorbed",8==0)){
1403	part.nabsorbed = part.nabsorbed +1;
1404	part.effet = 7;
1405	}
1406	/*
1407	if (PROC(1:3).eq.'pne')
1408	PROC_dan=7
1409	if (PROC(1:3).eq.'ppe')
1410	PROC_dan=8
1411	if (PROC(1:4).eq.'diff')
1412	PROC_dan=9
1413	if (PROC(1:4).eq.'ruth')
1414	PROC_dan=10
1415	*/
1416
1417
1418	// ?????????????????????????????????????
1419	if(part.effet == 2) {
1420	part_abs = 0;
1421	}else {
1422	part_abs = 1;
1423	}
1424
1425
1426	//=========================================================================================================================
1427	// Transform back to particle coordinates with opening and offset
1428
1429	if (part_abs == 0) { ////////////////////////// IF NUMERO 19//////////////////
1430	//
1431	// Include collimator tilt
1432
1433	if (tiltangle > 0) {
1434	x = x + tiltangle * c_length;
1435	xp = xp + tiltangle;
1436	} else if (tiltangle < 0) {
1437	x = x + tiltangle * c_length;
1438	xp = xp + tiltangle;
1439
1440	x = x - sin(tiltangle) * c_length;
1441	}
1442
1443	// Transform back to particle coordinates with opening and offset
1444
1445	z00 = z;
1446	x00 = x + mirror * c_offset;
1447	x = x + c_aperture / 2 + mirror * c_offset;
1448
1449	//+ Now mirror at the horizontal axis for negative X offset
1450
1451	x = mirror * x;
1452	xp = mirror * xp;
1453
1454	// Last do rotation into collimator frame
1455
1456	part.x = x * cos((-1) * c_rotation) + z * sin((-1) * c_rotation);
1457	part.y = z * cos((-1) * c_rotation) - x * sin((-1) * c_rotation);
1458	part.xp = xp * cos((-1) * c_rotation) + zp * sin((-1) * c_rotation);
1459	part.yp = zp * cos((-1) * c_rotation) - xp * sin((-1) * c_rotation);
1460
1461
1462
1463	//**********************************************************************************
1464	//For output.... by Zhang @ CERN, 24/04/2014
1465	//p_in = (1 + dpop) * p0;
1466	//From Dianiele
1467	//p_in = p;
1468	//s_in = s_in + s;
1469	part.PC = p;
1470	part.s = part.s + s;
1471
1472	if (part.effet == 1) { ////////////////////////// IF NUMERO 21//////////////////
1473
1474	part.x = 99.99e-3;
1475	part.y = 99.99e-3;
1476	cout<< "Mean that the particules is lost"<<endl;
1477	} ////////////////////////// fin IF NUMERO 21//////////////////
1478	} ////////////////////////// fin IF NUMERO 19//////////////////
1479
1480
1481	// End of check for particles not being lost before (see @330)
1482
1483	//} ////////////////////////// fin IF NUMERO 12//////////////////
1484	// End of loop over all particles
1485
1486
1487
1488	} //collimator with length = 0 ////////////////////////// fin IF NUMERO 1//////////////////
1489
1490	// =================================================================================FIN DE RETOUR AU COORDONEE DE ICOSIM====================================================================
1491
1492	part.xp=xp;
1493	part.yp=zp;
1494	part.y=z;
1495	part.x=x;
1496	part.PC=p;
1497
1498	xp1 = part.xp;
1499
1500
1501	//if(xp1-xp0!=0){
1502	//sortie << xp1 << "," << xp0<<"," << xfin<<"," << yfin<<endl;
1503	//}
1504	//cout<<"ICIIIIIIIIIIIIIIIIIIIIIIIII"<<part.PC<<endl;
1505
1506
1507	AV_xp1 = AV_xp1 + part.xp; //average x angle after interaction
1508	AV_yp1 = AV_yp1 + part.yp; //average y angle after interaction
1509	RMSxp1 = RMSxp1 + part.xp * part.xp; //rms x angle after interaction
1510	RMSyp1 = RMSyp1 + part.xp * part.xp; //rms y angle after interaction
1511	av_pc1 = av_pc1 + part.PC;
1512
1513	if(0){
1514	cout<<" X apres :"<<part.x<<" XP apres :"<<part.xp<<" Y apres :"<<part.y<<" YP apres :"<<part.yp<<" PC apres :"<<part.PC<<endl;
1515	//cout<<endl;
1516	}
1517	//only record the particles that hit the collimator.
1518	if(hitflag != true){
1519	//sortieIco.setf(ios::scientific);
1520	sortieIco << part.x << "," << part.y << "," << part.xp << "," << part.yp << "," << part.PC << endl;
1521	//}//end of for....
1522	}
1523
1524	}
1525	} //End of while
1526
1527
1528	} else {
1529	cerr<< "We can not open the file !" << endl;
1530	}
1531	entree.close();
1532
1533
1534
1535	AV_xp0 = AV_xp0 / count;
1536	AV_yp0 = AV_yp0 / count;
1537
1538	RMSxp0 = sqrt(RMSxp0 / (count));
1539	RMSyp0 = sqrt(RMSyp0 / (count));
1540
1541	SIGxp0 = sqrt((RMSxp0 * RMSxp0) - (AV_xp0 * AV_xp0));
1542	SIGyp0 = sqrt((RMSyp0 * RMSyp0) - (AV_yp0 * AV_yp0));
1543
1544	AV_xp1 = AV_xp1 / count;
1545	AV_yp1 = AV_yp1 / count;
1546	RMSxp1 = sqrt(RMSxp1 / (count));
1547	RMSyp1 = sqrt(RMSyp1 / (count));
1548
1549	SIGxp1 = sqrt((RMSxp1 * RMSxp1) - (AV_xp1 * AV_xp1));
1550	SIGyp1 = sqrt((RMSyp1 * RMSyp1) - (AV_yp1 * AV_yp1));
1551	av_pc0 = av_pc0 / count;
1552	av_pc1 = av_pc1 / count;
1553	//--------------------------------------------------------------------------
1554	//-------------write a summary file-----------------------------------------
1555	//cout<<endl;
1556	//cout<<endl;
1557	//cout<<endl;
1558
1559	//cout<<"avg xp_in: "<< AV_xp0 <<"// avg xp_out: "<< AV_xp1<<endl;
1560	//cout<<"rms xp_in: "<< RMSxp0 <<"// rms xp_out: "<< RMSxp1<<endl;
1561	//cout<<"sigma xp_in: "<< SIGxp0 <<"// sigma xp_out: "<< SIGxp1<<endl;
1562	//cout<<"avg pc_in: "<< av_pc0 <<"// avg pc_out: "<< av_pc1<<endl;
1563
1564
1565
1566	sortieIco.close();
1567	// cout<<endl;
1568	// affiche();
1569	file_out(pas, outputpath);
1570
1571	if(1){
1572	cout << "number of particles before the crystal: " << count<< endl;
1573	cout << "number of particle hit the crystal: " << part.nhit << endl;
1574	cout << "number of particles after the crystal: " << count - part.nhit << endl;
1575	}
1576	//cout<<endl;
1577	//cout<<endl;
1578	//cout<<endl;
1579	}
1580
1581
1582
1583	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1584	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1585	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1586	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1587
1588	//WE WRITE SOME FUNCTIONS IMPORTANT FUNCTION NOW
1589
1590	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1591	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1592	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1593	///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1594
1595
1596	void SimCrys::move_am(int IS, int NAM, double DZ, double DEI, double DLY, double DLR, double& xp, double& yp, double& PC)
1597	{
1598	//cout<<"The move_am function is use !!!!! "<<endl;
1599	double DYA(0);
1600	double Wp(0);
1601	if (NAM == 0) {
1602	return;
1603	}
1604	xp = xp * 1000;
1605	yp = yp * 1000;
1606	//cout << "xp avt: "<<xp<<" yp avt: "<<yp<<endl;
1607	//DEI ---> dE/dx (stoping energy)
1608	PC = PC - DEI * DZ; //energy lost beacause of ionization process[GeV]
1609
1610	// Coulomb Scattering
1611	//DYA ---> rms of coloumb scattering
1612	DYA = (14 / PC) * sqrt(DZ / DLR); //rms scattering (mrad)
1613	xp = xp + DYA * random_gauss_cut();
1614	yp = yp + DYA * random_gauss_cut();
1615
1616	//Elastic scattering
1617	//cout<<"Plusieurs possibilitees : "<<endl;
1618	if (random() <= (DZ / crys.DLAI[IS])) {
1619	//cout<<"Poss. 1 !!!! "<<endl;
1620	xp = xp + 196 * random_gauss_cut() / PC; //angle elast. scattering in R plane
1621	yp = yp + 196 * random_gauss_cut() / PC;
1622	}
1623
1624	//Diffraction interaction
1625	if (random() <= (DZ / (DLY * 6.143))) {
1626	//cout<<"Poss. 2 !!!! "<<endl;
1627	xp = xp + 257 * random_gauss_cut() / PC; // angle elast. scattering in R plane[mr]
1628	yp = yp + 257 * random_gauss_cut() / PC;
1629	Wp = random();
1630	PC = PC - 0.5 * pow(0.3 * PC, Wp); // m0*c = 1 GeV/c for particle
1631
1632	}
1633
1634	//Inelastic interaction
1635	if (random() <= DZ / DLY) {
1636	//cout<<" Absorbed !! "<< endl;
1637	}
1638	//cout<<"XP APRES : "<<xp<<" YP APRES : "<<yp<<endl;
1639	xp = xp / 1000;
1640	yp = yp / 1000;
1641	}
1642
1643
1644	double SimCrys::random()
1645	{
1646	double scale = RAND_MAX;
1647	double base = rand() / scale;
1648	double fine = rand() / scale;
1649	return base + fine / scale;
1650	}
1651
1652	double SimCrys::random_dist()
1653	{
1654	double scale = RAND_MAX;
1655	double base = rand() / scale;
1656	double fine = rand() / scale;
1657	return (-2.5 * 1e-4 + (3.5 * 1e-4 ) * (base + fine / scale));
1658	}
1659
1660	double SimCrys::random_gauss()
1661	{
1662	const double PI (4.0 * atan(1.0));
1663	double U(random());
1664	double V(random());
1665	return (sqrt(-2 * log(U)) * cos(2 * PI * V));
1666	}
1667
1668	double SimCrys::random_gauss_cut()
1669	{
1670	double resu;
1671	do {
1672	resu = random_gauss();
1673	} while (abs(resu) >= 1);
1674	return resu;
1675	}
1676
1677	double SimCrys::random_gauss_dist() //POUR faire varier la moyenne de la gaussienne il faut mettre en argu la sim& !!!!!
1678	{
1679	const double PI (4.0 * atan(1.0));
1680	double U(random());
1681	double V(random());
1682	//moy=random();
1683	return 1e-4 * (sqrt(-2 * log(U)) * cos(2 * PI * V)); //+(0.0002*moy); //POUR faire varier la moyenne de la gaussienne aussi !!!!!
1684	}
1685
1686	/////////////////////////////////*******************************************************+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1687	//POUR 409 :
1688
1689	double SimCrys::random_gauss_dist_409() //POUR faire varier la moyenne de la gaussienne il faut mettre en argu la sim& !!!!!
1690	{
1691	const double PI (4.0 * atan(1.0));
1692	double U(random());
1693	double V(random());
1694	//moy=random();
1695	return (8.939203e-06) * (sqrt(-2 * log(U)) * cos(2 * PI * V)) + (-7.183296e-08); //POUR faire varier la moyenne de la gaussienne aussi !!!!!
1696	}
1697
1698	///////////////////////////////////************************************************************++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1699
1700	void SimCrys::affiche()
1701	{
1702	cout << "Crystal Curved Length [m] :" << crys.Cry_length << endl;
1703	cout << "Crystal Curvature Radius [m] :" << crys.Rcurv << endl;
1704	cout << "Crystal Bending Angle [urad] :" << crys.cry_bend << " , " << crys.Cry_length / crys.Rcurv << endl;
1705	cout << "Critical Radius :" << crys.Rcrit << endl;
1706	cout << "Ratio :" << crys.ratio << endl;
1707	cout << "Critical Angle for straight Crystal :" << crys.xpcrit0 << endl;
1708	cout << "Critical Angle for curved Crystal :" << crys.xpcrit << endl;
1709	cout << "Angle average reflexion : " << part.Ang_avr << endl;
1710	cout << "Full Channeling : " << crys.Cry_length / crys.Rcurv << endl;
1711	cout << endl;
1712	cout << endl;
1713	cout << "N.AMORPH :" << part.namor << endl;
1714	cout << "N.DECHANNELING:" << part.ndechann << endl;
1715	cout << "N.CHANNELING :" << part.nchann << endl;
1716	cout << "N.OUT :" << part.nout << endl;
1717	cout << "N.VREFLECTION :" << part.nvrefl << endl;
1718	cout << "N.VCAPTURE:" << part.nvcapt << endl;
1719	cout << "N.ABSORBED: "<< part.nabsorbed << endl;
1720
1721	//cout << "Si la somme ne donne pas le nbre de part., c est normale !! En effet, il est possible que plusieurs evenement se produisent!!!";
1722	}
1723
1724	void SimCrys::file_out(const int& pas, string outputpath)
1725	{
1726
1727	ofstream out;
1728	string file;
1729	file = outputpath + "/crystal_output.out";
1730
1731	out.open(file.c_str(), ios::out \| ios::app);
1732
1733	if (out.fail()) {
1734	cerr << "Warning: problem writing in the file " << file << "!" << endl;
1735	} else {
1736
1737	out << endl;
1738	out << "Crystal Curved Length [m] :" << crys.Cry_length << endl;
1739	out << "Crystal Curvature Radius [m] :" << crys.Rcurv << endl;
1740	out << "Crystal Bending Angle [urad] :" << crys.cry_bend << " , " << 1e6*crys.Cry_length / crys.Rcurv << endl;
1741	out << "Critical Radius :" << crys.Rcrit << endl;
1742	out << "Ratio :" << crys.ratio << endl;
1743	out << "Critical Angle for straight Crystal :" << crys.xpcrit0 << endl;
1744	out << "Critical Angle for curved Crystal :" << crys.xpcrit << endl;
1745	out << "Angle average reflexion : " << part.Ang_avr << endl;
1746	out << "Full Channeling : " << crys.Cry_length / crys.Rcurv << endl;
1747	out << endl;
1748	out << endl;
1749	out << "N.AMORPH :" << part.namor << endl;
1750	out << "N.DECHANNELING:" << part.ndechann << endl;
1751	out << "N.CHANNELING :" << part.nchann << endl;
1752	out << "N.OUT :" << part.nout << endl;
1753	out << "N.VREFLECTION :" << part.nvrefl << endl;
1754	out << "N.VCAPTURE :" << part.nvcapt << endl ;
1755	out << "N.ABSORBED: "<< part.nabsorbed << endl;
1756	out << "N. of particles that hit the crystal collimator: " << part.nhit << endl;
1757	cout << endl;
1758	cout << endl;
1759
1760	out.close();
1761	}
1762
1763	}
1764

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

Context Navigation

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

Download in other formats: