Context Navigation

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

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

Line
1	#include <iostream>
2	#include <iomanip>
3	#include <fstream>
4	#include <sstream>
5	#include <algorithm>
6	#include <vector>
7	#include <string>
8	#include <cmath>
9	#include "simulation.h"
10
11	//Port for the link with Fluka with the script run.sh
12	#if defined(FLUKA)
13	#define PORT 12345
14	#endif
15	using namespace std;
16
17
18	//===============================Constructors, destructor======================
19
20
21	Simulation::Simulation(int size_bunch, int nbre_elts, string genere)
22	: sizeBunch(size_bunch), nbre_elts(nbre_elts)
23	{
24	double ax(-2.316794521);
25	double bx(103.4788975);
26	double dx(1.231007397);
27	double dpx(0.027602708);
28	double by(20.8994661);
29	double ay(0.537534425);
30	double dy(0);
31	double dpy(0);
32	srand(time(0));
33
34	Particle p1;
35
36	for (int l(0); l < sizeBunch; ++l) {
37
38	p1.genpartdist(1, 1, genere, 1.00095, 1.2e-8, 1.2e-8, 0.00011, bx, ax, dx, dpx, by, ay, dy, dpy, 6);
39	p1.afficheFirstCoordonnees();
40
41	bunch.push_back(p1);
42	}
43	}
44
45	Simulation::~Simulation()
46	{
47
48	for (int j(0); j < nbre_elts; ++j) {
49	delete lat.reseau[j];
50	}
51	}
52
53
54	void Simulation::die(const char* msg)
55	{
56	cerr << msg << endl;
57
58	#if defined(FLUKA)
59	if (lat.conn) {
60	flukaio_disconnect(lat.conn);
61	lat.conn = NULL;
62	}
63	#endif
64
65	exit(1);
66	}
67
68
69	void Simulation::readParticle(const double& Apr, const double& Zpr, string inputfile)
70	{
71
72	string particlefile;
73
74	particlefile = inputfile + "initial.dat";
75
76	//if the name of the file is not 'initial.dat', uncomment the following two lines
77
78	//cout << "Enter the name of the file containing the particles informations: " << endl;
79
80	//cin >> particlefile;
81
82	ifstream enterPart;
83
84	enterPart.open(particlefile.c_str());
85
86	if (enterPart.fail()) {
87	cerr << "Warning: problem reading the file " << particlefile << "!" << endl;
88	} else {
89
90	double temp2;
91
92	enterPart >> temp2;
93
94	while (!enterPart.eof()) {
95
96	Particle p;
97
98	p.coordonnees[0][0] = temp2;
99	enterPart >> p.coordonnees[0][1];
100	enterPart >> p.coordonnees[0][2];
101	enterPart >> p.coordonnees[0][3];
102	enterPart >> p.coordonnees[0][4];
103	p.Ap0 = Apr;
104	p.Zp0 = Zpr;
105	p.coordonnees[0][5] = 0;
106	p.dt = 0;
107	p.dpoporiginal = p.coordonnees[0][4];
108
109
110	bunch.push_back(p);
111
112	enterPart >> temp2;
113	}
114	}
115	}
116
117	//read information from collimator file
118
119	void Simulation::readInputFile(string inputfile, double& Apr, double& Zpr, double& mass, double& energyPerIon, double& emix, double& emiy, int& npart, int& kbunch, double& sigdpp, int& taubeam, int& nparti, string& partdistr, double& r1r2skin, int& nrev1, int& nrev2, double& blowup2, int& blowupperiod, string& stoplineartracking, int& scaleorbit, double& wecolli, double& nsigi, string& opticsFile, double& thicknessmagneticfield, double& Bmax, double& deltaGap, int& beamflag, int& RunWithFluka, int& RunWithCrystal, double& freqrf, double& rfvoltage, double& rfharmonic, int& RunningFlag, int& idpart, int& idelt, int& outcoord, double& momentum, string& plotflag, string& inputpath, int& RFflag)
120	{
121
122	ifstream entree;
123	string crosssectionspath;
124	entree.open(inputfile.c_str());
125
126	if (entree.fail()) {
127	cerr << "Warning: problem to read the collimator file!!" << endl;
128	} else {
129
130	//reading the parameter values in the top part of the collimator file
131
132	string name, comment, date, strtemp;
133	entree >> ws;
134	getline(entree, comment);
135	getline(entree, date);
136
137	while (!entree.eof()) {
138
139	entree >> ws;
140
141	getline(entree, name, ',');
142
143	if (name == "MASSNUMBER") {
144	getline(entree, strtemp);
145	Apr = atof(strtemp.c_str());
146	} else if (name == "CHARGESTATE") {
147	getline(entree, strtemp);
148	Zpr = atof(strtemp.c_str());
149	} else if (name == "MASS") {
150	getline(entree, strtemp);
151	mass = atof(strtemp.c_str());
152	} else if (name == "ENERGY") {
153	getline(entree, strtemp);
154	energyPerIon = atof(strtemp.c_str());
155	} else if (name == "EX") {
156	getline(entree, strtemp);
157	emix = atof(strtemp.c_str());
158	} else if (name == "EY") {
159	getline(entree, strtemp);
160	emiy = atof(strtemp.c_str());
161	} else if (name == "NPART") {
162	getline(entree, strtemp);
163	npart = atoi(strtemp.c_str());
164	} else if (name == "KBUNCH") {
165	getline(entree, strtemp);
166	kbunch = atoi(strtemp.c_str());
167	} else if (name == "SIGDPP") {
168	getline(entree, strtemp);
169	sigdpp = atof(strtemp.c_str());
170	} else if (name == "TAUBEAM") {
171	getline(entree, strtemp);
172	taubeam = atoi(strtemp.c_str());
173	} else if (name == "NPARTI") {
174	getline(entree, strtemp);
175	nparti = atoi(strtemp.c_str());
176	} else if (name == "PARTDISTR") {
177	getline(entree, partdistr);
178	} else if (name == "R1R2SKIN") {
179	if (partdistr == "r1r2") {
180	getline(entree, strtemp);
181	r1r2skin = atof(strtemp.c_str());
182	} else {
183	getline(entree, strtemp);
184	}
185	} else if (name == "NREV1") {
186	getline(entree, strtemp);
187	nrev1 = atoi(strtemp.c_str());
188	} else if (name == "NREV2") {
189	getline(entree, strtemp);
190	nrev2 = atoi(strtemp.c_str());
191	} else if (name == "BLOWUP2") {
192	getline(entree, strtemp);
193	blowup2 = atof(strtemp.c_str());
194	} else if (name == "BLOWUPPERIOD") {
195	getline(entree, strtemp);
196	blowupperiod = atoi(strtemp.c_str());
197	} else if (name == "STOPLINEARTRACKING") {
198	getline(entree, stoplineartracking);
199	} else if (name == "SCALEORBIT") {
200	getline(entree, strtemp);
201	scaleorbit = atoi(strtemp.c_str());
202	} else if (name == "WECOLLI") {
203	getline(entree, strtemp);
204	wecolli = atof(strtemp.c_str());
205	} else if (name == "NSIGI") {
206	getline(entree, strtemp);
207	nsigi = atof(strtemp.c_str());
208	} else if (name == "CROSSSECTIONSPATH") {
209	getline(entree, crosssectionspath);
210	} else if (name == "OPTICSPATH") {
211	getline(entree, opticsFile);
212	} else if (name == "INPUTPATH") {
213	getline(entree, inputpath);
214	} else if (name == "THICKNESSMAGNETICFIELD") {
215	getline(entree, strtemp);
216	thicknessmagneticfield = atof(strtemp.c_str());
217	} else if (name == "BMAX") {
218	getline(entree, strtemp);
219	Bmax = atof(strtemp.c_str());
220	} else if (name == "DELTAGAP") {
221	getline(entree, strtemp);
222	deltaGap = atof(strtemp.c_str());
223	} else if (name == "BEAMFLAG") {
224	getline(entree, strtemp);
225	beamflag = atoi(strtemp.c_str());
226	} else if (name == "RUNNINGFLAG") {
227	getline(entree, strtemp);
228	RunningFlag = atoi(strtemp.c_str());
229	} else if (name == "IDPART") {
230	getline(entree, strtemp);
231	idpart = atoi(strtemp.c_str());
232	} else if (name == "OUTCOORD") {
233	getline(entree, strtemp);
234	outcoord = atoi(strtemp.c_str());
235	} else if (name == "IDELT") {
236	getline(entree, strtemp);
237	idelt = atoi(strtemp.c_str());
238	} else if (name == "FREQRF") {
239	getline(entree, strtemp);
240	freqrf = atof(strtemp.c_str());
241	RFflag = 1;
242	} else if (name == "RFVOLTAGE") {
243	getline(entree, strtemp);
244	rfvoltage = atof(strtemp.c_str());
245	RFflag = 1;
246	} else if (name == "RFHARMONBER") {
247	getline(entree, strtemp);
248	rfharmonic = atof(strtemp.c_str());
249	RFflag = 1;
250	} else if (name == "MOMENTUM") {
251	getline(entree, strtemp);
252	momentum = atof(strtemp.c_str());
253	} else if (name == "PLOTFLAG") {
254	getline(entree, plotflag);
255	} else if (name == "NAME") {
256	break;
257	} else {
258	getline(entree, strtemp);
259	}
260
261	}
262	}
263	entree.close();
264
265	//reads optics file
266	extractOpticsParameters(opticsFile, beamflag);
267
268
269	//order of the headers in collimatorfile: NAME,ANGLE,LENGTH,PHASE,NSIG,MATERIAL,METHOD
270
271	ifstream entree1;
272	vector <Collimator*> temp;
273
274	entree1.open(inputfile.c_str());
275
276	if (entree1.fail()) {
277	cerr << "Warning: problem reading the second part of the collimator file!" << endl;
278	} else {
279
280	string NAME, MATERIAL, METHOD, strtemp;
281	double ANGLE, LENGTH;
282	int PHASE;
283	double NSIG2;
284
285
286	getline(entree1, NAME, ',');
287	while (NAME != "NAME") {
288	getline(entree1, strtemp);
289	getline(entree1, NAME, ',');
290	}
291	getline(entree1, strtemp);
292
293	while (!entree1.eof()) {
294	getline(entree1, NAME, ',');
295	getline(entree1, strtemp, ',');
296	ANGLE = atof(strtemp.c_str());
297	getline(entree1, strtemp, ',');
298	LENGTH = atof(strtemp.c_str());
299	getline(entree1, strtemp, ',');
300	PHASE = atoi(strtemp.c_str());
301	getline(entree1, strtemp, ',');
302	NSIG2 = atof(strtemp.c_str());
303	getline(entree1, MATERIAL, ',');
304	getline(entree1, METHOD);
305
306	//we construct some special vectors related to collimators
307
308
309	for (int i(0); i < nbre_elts; ++i) {
310	if (lat.reseau[i]->NAME == NAME) {
311	//we skip collimators that are not found in the optics file and collimators with names ending in the wrong letter
312	if ((atoi(&NAME[NAME.size() - 1]) == beamflag) \|\| (NAME == "TCTV.4R2.B2") \|\| (NAME == "TCTV.4R8.B2") \|\| (NAME == "TCTV.4L2.B1") \|\| (NAME == "TCTV.4L8.B1")) { //attention, il faut peut etre egalement supprimer ces elements dans le reseau
313
314	lat.ips.push_back(i);
315
316	if (NAME.substr(0, 3) == "TCI") {
317	lat.ii.push_back(i);
318	} else if ((NAME.substr(0, 3) == "TCP") \|\| (NAME.substr(0, 7) == "CRYSTAL") \|\| (NAME.substr(0, 5) == "TCRYO")) {
319	lat.ip.push_back(i);
320	} else if (NAME.substr(0, 3) == "TCS") {
321	lat.is.push_back(i);
322	} else if (NAME.substr(0, 3) == "TCT") {
323	lat.it.push_back(i);
324	} else if (NAME.substr(0, 4) == "TCLA") {
325	lat.itcla.push_back(i);
326	}
327
328
329	if (METHOD == "standard") {
330
331	StandardCollimator stdcolli(*lat.reseau[i], ANGLE, NSIG2, METHOD, MATERIAL);
332	temp.push_back(new StandardCollimator(stdcolli));
333
334	} else if (METHOD == "magnetic") {
335
336	MagneticCollimator magnetcolli(*lat.reseau[i], ANGLE, NSIG2, METHOD, MATERIAL);
337	temp.push_back(new MagneticCollimator(magnetcolli));
338
339	} else if (METHOD == "fluka") {
340	#if defined(FLUKA)
341	RunWithFluka = 1;
342	FlukaCollimator flukacolli(*lat.reseau[i], ANGLE, NSIG2, METHOD);
343	temp.push_back(new FlukaCollimator(flukacolli));
344	#else
345	die("Fluka coupling not supported in this version, please compile with the appropriate flags");
346	#endif
347	} else if (METHOD == "crystal") {
348	RunWithCrystal = 1;
349	CrystalCollimator crystalcolli(*lat.reseau[i], ANGLE, NSIG2, METHOD);
350	temp.push_back(new CrystalCollimator(crystalcolli));
351
352	} else {
353
354	cerr << "Warning: unknown method (" << METHOD << ") for the collimator " << NAME << endl;
355	}
356	}
357	}
358	}
359	}
360	}
361
362
363	sort(lat.ip.begin(), lat.ip.end());
364	sort(lat.ii.begin(), lat.ii.end());
365	sort(lat.is.begin(), lat.is.end());
366	sort(lat.it.begin(), lat.it.end());
367	sort(lat.itcla.begin(), lat.itcla.end());
368
369	//sorting all the information according to position along the ring
370
371	sort(lat.ips.begin(), lat.ips.end());
372
373	vector <int> indice;
374	vector <double> pos, pos1;
375
376	for (int k(0); k < temp.size(); ++k) {
377	pos.push_back(temp[k]->S);
378	pos1.push_back(temp[k]->S);
379	}
380
381	sort(pos.begin(), pos.end());
382
383	for (int j(0); j < pos.size(); ++j) {
384	int m(0);
385	for (int k(0); k < pos1.size(); ++k) {
386	if (pos[j] == pos1[k]) {
387	indice.push_back(m);
388	pos1[k] = -1;
389	break;
390	} else {
391	++m;
392	}
393	}
394	}
395
396	for (int j(0); j < temp.size(); ++j) {
397	if (temp[indice[j]]->method == "standard") {
398
399	lat.addCollimator(temp[indice[j]]);
400	lat.resColli[j]->crossSectionPath = crosssectionspath + '/' + lat.resColli[j]->material + '/';
401
402	} else if (temp[indice[j]]->method == "magnetic") {
403
404	lat.addCollimator(temp[indice[j]]);
405	lat.resColli[j]->crossSectionPath = crosssectionspath + '/' + lat.resColli[j]->material + '/';
406	lat.resColli[j]->Bmax = Bmax;
407	lat.resColli[j]->thicknessMagneticField = thicknessmagneticfield;
408	lat.resColli[j]->energyPerIon = energyPerIon;
409	lat.resColli[j]->mass = mass;
410	lat.resColli[j]->deltaGap = deltaGap;
411
412	} else if (temp[indice[j]]->method == "fluka") {
413
414	lat.addCollimator(temp[indice[j]]);
415
416	} else if (temp[indice[j]]->method == "crystal") {
417
418	lat.addCollimator(temp[indice[j]]);
419
420	ifstream readinfocrystal;
421	string inputtemp;
422	inputtemp = inputpath + "crystalinfo.csv";
423	readinfocrystal.open(inputtemp.c_str());
424
425	if (readinfocrystal.fail()) {
426	cerr << "Warning: problem reading the file crystalinfo.csv!" << endl;
427	} else {
428	string temp1;
429	getline(readinfocrystal, temp1, ',');
430	while (!readinfocrystal.eof()) {
431	if (temp[indice[j]]->NAME == temp1) {
432	getline(readinfocrystal, temp1, ',');
433	lat.resColli[j]->C_orient = atoi(temp1.c_str());
434	getline(readinfocrystal, temp1, ',');
435	lat.resColli[j]->IS = atoi(temp1.c_str());
436	getline(readinfocrystal, temp1, ',');
437	lat.resColli[j]->C_xmax = atof(temp1.c_str());
438	getline(readinfocrystal, temp1, ',');
439	lat.resColli[j]->C_ymax = atof(temp1.c_str());
440	getline(readinfocrystal, temp1, ',');
441	lat.resColli[j]->Cry_length = atof(temp1.c_str());
442	getline(readinfocrystal, temp1, ',');
443	lat.resColli[j]->Rcurv = atof(temp1.c_str());
444	getline(readinfocrystal, temp1, ',');
445	lat.resColli[j]->C_rotation = atof(temp1.c_str());
446	getline(readinfocrystal, temp1, ',');
447	lat.resColli[j]->C_aperture = atof(temp1.c_str());
448	getline(readinfocrystal, temp1, ',');
449	lat.resColli[j]->C_offset = atof(temp1.c_str());
450	getline(readinfocrystal, temp1, ',');
451	lat.resColli[j]->C_tilt = atof(temp1.c_str());
452	getline(readinfocrystal, temp1);
453	lat.resColli[j]->Cry_tilt = atof(temp1.c_str());
454	} else {
455	getline(readinfocrystal, temp1);
456	}
457	getline(readinfocrystal, temp1, ',');
458	}
459	readinfocrystal.close();
460	}
461
462	} else {
463
464	cerr << "Warning: problem to construct resColli " << endl;
465	}
466	}
467
468	for (int j(0); j < lat.resColli.size(); ++j) {
469	lat.resColli[j]->tcang = lat.resColli[j]->tcang * M_PI / 180;
470	}
471
472	vector <double> sx, sy;
473
474	//In this calculation we do not take into account that there is variation in the momentum! In that case the formulas would be:
475	//sx=nsigsqrt(bxemix+dx^2sigdpp^2) and sy=nsigsqrt(byemiy+dy^2sigdpp^2)
476
477	for (int k(0); k < lat.resColli.size(); ++k) {
478
479	sx.push_back(lat.resColli[k]->nsig * sqrt(lat.resColli[k]->BETX * emix));
480	sy.push_back(lat.resColli[k]->nsig * sqrt(lat.resColli[k]->BETY * emiy));
481	}
482
483	double b[lat.resColli.size()];
484
485	for (int j(0); j < lat.resColli.size(); ++j) {
486	if (sin(lat.resColli[j]->tcang) == 0) {
487	b[j] = sx[j];
488	} else if (sin(lat.resColli[j]->tcang) > 0) { //collimator is straight
489	b[j] = sin(lat.resColli[j]->tcang) * sqrt(sy[j] * sy[j] + sx[j] * sx[j] / (tan(lat.resColli[j]->tcang) * tan(lat.resColli[j]->tcang)));
490	} else {
491	b[j] = 0;
492	}
493	}
494
495	if (wecolli >= 0) { //rotate collimator around face
496	for (int k(0); k < lat.resColli.size(); ++k) {
497	lat.resColli[k]->hgap = b[k];
498	lat.resColli[k]->hgap2 = b[k] + lat.resColli[k]->L * wecolli;
499	}
500	} else {//rotate collimator around end
501	for (int k(0); k < lat.resColli.size(); ++k) {
502	lat.resColli[k]->hgap = b[k] + lat.resColli[k]->L * (-wecolli);
503	lat.resColli[k]->hgap2 = b[k];
504	}
505	}
506
507	for (int k(0); k < lat.resColli.size(); ++k) {
508	lat.resColli[k]->phi = (lat.resColli[k]->hgap2 - lat.resColli[k]->hgap) / lat.resColli[k]->L; //Calculate absolute angle of collimator
509	}
510
511
512	for (int n(0); n < lat.ips.size(); ++n) {
513	for (int p(0); p < lat.ip.size(); ++p) {
514	if (lat.ips[n] == lat.ip[p]) {
515	lat.ipcoll.push_back(n);
516	}
517	}
518	}
519
520
521	lat.freqrf = freqrf;
522	lat.rfvoltage = rfvoltage;
523	lat.rfharmonic = rfharmonic;
524	lat.momentum = momentum;
525
526	}
527
528
529
530	void Simulation::run(string inputfile, string outputpath)
531	{
532
533	cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl;
534	cout << "++++++++++++++++++++++++++++++++++++++ NEW SIMULATION +++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl;
535	cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl;
536
537
538	double Apr, Zpr, mass, energyPerIon, emix, emiy, sigdpp, r1r2skin, blowup2, thicknessmagneticfield, Bmax, deltaGap, gamma, betgam, W, PlossPb, freqrf, rfvoltage, rfharmonic, wecolli, nsigi, momentum;
539	int npart, kbunch, taubeam, nparti, nrev1, nrev2, blowupperiod, scaleorbit, beamflag, stopLinearTrackingNo, RunWithFluka, RunWithCrystal, RunningFlag, outcoord, idpart, idelt, RFflag;
540	string partdistr, stoplineartracking, crosssectionspath, opticsFile, plotflag, inputpath;
541
542	int indication(1);
543	outcoord = 0;
544
545	RunWithFluka = 0;
546	RunWithCrystal = 0;
547	RFflag = 0;
548
549	cout << "READING OF THE INPUT FILES." << endl;
550
551	//getting information about the optics, apertures and collimators from the input file
552	try {
553	readInputFile(inputfile, Apr, Zpr, mass, energyPerIon, emix, emiy, npart, kbunch, sigdpp, taubeam, nparti, partdistr, r1r2skin, nrev1, nrev2, blowup2, blowupperiod, stoplineartracking, scaleorbit, wecolli, nsigi, opticsFile, thicknessmagneticfield, Bmax, deltaGap, beamflag, RunWithFluka, RunWithCrystal, freqrf, rfvoltage, rfharmonic, RunningFlag, idpart, idelt, outcoord, momentum, plotflag, inputpath, RFflag);
554	}
555
556
557	catch (int erreur) {
558	cerr << "Error 1: no aperture infos!!" << endl;
559	cerr << "The program has ended without stopping the fluka server. You have to kill it by yourself (pkill flukaserver)." << endl;
560	}
561
562	momentum = sqrt(energyPerIon * energyPerIon - (mass * mass));
563	momentum = momentum * 1e09;
564
565	if ((RunningFlag != 1) && (RunningFlag != 2) && (RunningFlag != 3) && (RunningFlag != 4)) {
566	cerr << "Warning: Running Flag must be equal to 1, 2, 3 or 4!!" << endl;
567	cerr << "Check it in the collimator file, if it as no RUNNINGFLAG parameter, you must add it" << endl;
568	return;
569	}
570
571	gamma = energyPerIon / mass;
572	betgam = sqrt(gamma * gamma - 1);
573
574	//The beginning of the tracking with 'trackensemblechrom' is linear. We decide when we start doing the nonlinear tracking.
575
576	for (int i(0); i < lat.reseau.size(); ++i) {
577	if (stoplineartracking == lat.reseau[i]->NAME) {
578	stopLinearTrackingNo = i;
579	break;
580	}
581	}
582
583	if (stoplineartracking.size() == 0) {
584
585	cerr << "Warning: the element where we want to stop the linear tracking does not exist!" << endl;
586	}
587
588	W = npart;
589	W = W * kbunch;
590	W = W * energyPerIon;
591	W = W * 1000000000;
592	W = W * 1.602e-19; //total energy of beam
593	PlossPb = W / taubeam; // average power loss of beam
594
595	lat.setsize_res(lat.reseau.size());
596
597	if ((RunWithCrystal == 1) && ((Apr != 1) \|\| (Zpr != 1))) {
598	cerr << "Warning: Crystal collimation can only be used with protons!!" << endl;
599	return;
600	}
601
602	int i0(0);
603
604	cout << "GENERATION OF THE BUNCH OF PARTICLES." << endl;
605
606	//===============================================================================================================================================================================================
607	//THE FOLLOWING PART CAN BE COMMENTED/UNCOMMENTED TO COMPARE EXACT VALUES WITHOUT RANDOMNESS WITH ICOSIM (SEE THE FILE RANDOMNESS.TXT FOR MORE EXPANATION ON WHAT TO DO TO ELIMINATE RANDOMNESS)
608
609	//WE ALSO HAVE THE POSSIBILITY HERE TO READ THE PARTICLES FROM A FILE WITH THE FUNCTION readParticle
610	//===============================================================================================================================================================================================
611
612	if (RunningFlag == 2) {
613
614	srand(time(0));
615
616	readParticle(Apr, Zpr, inputfile);
617
618	//we set the id of the particles
619	for (int y(0); y < bunch.size(); ++y) {
620	bunch[y].setidentification(y);
621	}
622
623	} else if (RunningFlag == 1) {
624
625	genbunch(i0, nparti, partdistr, r1r2skin, nsigi * nsigi * emix, nsigi * nsigi * emiy, sigdpp, lat.reseau[i0]->BETX, lat.reseau[i0]->ALFX, lat.reseau[i0]->DX, lat.reseau[i0]->DPX, lat.reseau[i0]->BETY, lat.reseau[i0]->ALFY, lat.reseau[i0]->DY, lat.reseau[i0]->DPY, nsigi);
626
627	//we set the id of the particles
628	for (int y(0); y < bunch.size(); ++y) {
629	bunch[y].setidentification(y);
630	}
631
632	}
633
634	/*Particle p1(-9.62058e-07, -0.000153021, -4.97476e-05, -4.80766e-05, -8.21051e-05, 0, 208, 82, -0.00012159388136310807);
635	Particle p2(-2.29177e-05, -0.000135618, -6.02201e-05, 3.8913e-05, -6.28236e-05, 0, 208, 82, -6.569201361245118e-05);
636	Particle p3(1e-06, -0.00015, -5e-05, 2e-05, -7e-05, 0, 208, 82, -7e-05);
637	Particle p4(2e-05, -0.00012, 2e-07, 4e-06, 5e-05, 0, 208, 82, 5e-05);
638	Particle p5(-3e-05, -0.00016, 3e-05, -2e-05, -3e-06, 0, 208, 82, -3e-06);
639	Particle p6(-1e-06, -0.00013, -2e-05, -2e-05, 2e-05, 0, 208, 82, 2e-05);
640	Particle p7(1.5e-05, -0.00015, -3e-06, 1.5e-05, 6e-05, 0, 208, 82, 6e-05);
641	Particle p8(8e-06, -0.00013, 4e-06, 2.3e-05, -6e-05, 0, 208, 82, -6e-05);
642	Particle p9(-4e-05, -0.00015, 5e-05, 4e-05, 5e-05, 0, 208, 82, 5e-05);
643	Particle p10(5e-05, -0.00016, 6e-07, -5e-05, -3e-05, 0, 208, 82, -3e-05);
644
645
646
647	bunch.push_back(p1);
648	bunch.push_back(p2);
649	bunch.push_back(p3);
650	bunch.push_back(p4);
651	bunch.push_back(p5);
652	bunch.push_back(p6);
653	bunch.push_back(p7);
654	bunch.push_back(p8);
655	bunch.push_back(p9);
656	bunch.push_back(p10);*/
657
658
659	//we set the id of the particles
660	/*for(int y(0); y<bunch.size(); ++y){
661	bunch[y].setidentification(y);
662	}*/
663
664	//================================================================= END OF TEMP. MODIFICATIONS ===========================================================================================
665	int particlehit;
666
667	if ((RunningFlag == 1) \|\| (RunningFlag == 2)) {
668	cout << "FIRST TRACKING." << endl;
669
670	//We track the particles linearly between the primary collimators
671	lat.trackensemblelinearnew(bunch, bunchhit, nrev1, blowup2, blowupperiod);
672	particlehit = bunchhit.size();
673
674	for (int i(0); i < bunchhit.size(); ++i) {
675	bunchhit[i].Ap0 = Apr;
676	bunchhit[i].Zp0 = Zpr;
677	}
678	}
679
680
681	if (RunningFlag == 4) {
682	srand(time(0));
683	readParticle(Apr, Zpr, inputfile);
684	} else if (RunningFlag == 3) {
685	genbunch(i0, nparti, partdistr, r1r2skin, nsigi * nsigi * emix, nsigi * nsigi * emiy, sigdpp, lat.reseau[i0]->BETX, lat.reseau[i0]->ALFX, lat.reseau[i0]->DX, lat.reseau[i0]->DPX, lat.reseau[i0]->BETY, lat.reseau[i0]->ALFY, lat.reseau[i0]->DY, lat.reseau[i0]->DPY, nsigi);
686	}
687
688	if ((RunningFlag == 3) \|\| (RunningFlag == 4)) {
689	for (int y(0); y < bunch.size(); ++y) {
690	bunch[y].setidentification(y);
691	}
692
693	for (int k(0); k < bunch.size(); ++k) {
694	bunchhit.push_back(bunch[k]);
695	}
696	particlehit = bunchhit.size();
697
698	for (int i(0); i < bunchhit.size(); ++i) {
699	bunchhit[i].Ap0 = Apr;
700	bunchhit[i].Zp0 = Zpr;
701	}
702	}
703	//#################################################################################################################################
704
705	#if defined(FLUKA)
706	//if we have a Fluka element
707	if (RunWithFluka == 1) { //connection to Fluka
708
709	int port = PORT;
710
711	lat.conn = flukaio_connect("pcen33148", port);//Warning: the hostname (pcen33148) must be changed to yours. Your hostname can be find in the file network.nfo under the directory ICOSIM++/<name_of_run.nber>/1
712
713	if (lat.conn == NULL) {
714	die("Error connecting to server");
715	}
716
717	cout << "Connected to server" << endl;
718	}
719	#endif
720
721
722	if (bunchhit.size() == 0) { //no hit during the first tracking
723	string outputfile;
724
725	cout << "All the particles go througt the accelerator. No hit!!" << endl;
726
727	#if defined(FLUKA)
728	if (RunWithFluka == 1) { //if we have a Fluka Collimator, we disconnect
729
730	cout << "End of the connection." << endl;
731
732	int n;
733	flukaio_message_t msg;
734	n = flukaio_send_eoc(lat.conn);
735	if (n < 0) {
736	die("Error occurred sending");
737	}
738
739	n = flukaio_wait_message(lat.conn, &msg);
740	if (n == -1) {
741	die("Server timeout when waiting end of computation");
742	}
743	if (msg.type != N_END) {
744	die("Unexpected message received");
745	}
746
747	flukaio_disconnect(lat.conn);
748	lat.conn = NULL;
749	}
750	#endif
751
752	return;
753
754
755	}
756
757
758	//We take the particles hitting collimators during the previous tracking (we take their coordinates at the beginning of the last turn before the hit), and do a more precise tracking
759	cout << "SECOND TRACKING (LINEAR IN THE BEGINNING OF THE FIRST REVOLUTION)" << endl;
760
761	vector <double> Aphito, Zphito, hitso;
762	vector <int> nhitcollio, nhitspoilero;
763	vector <vector <double> > xco, yco;
764	int irev(0);
765	int nonlinflag;
766	double attr1(0);
767	nonlinflag = 0;//the second tracking is linear
768
769	lat.trackensemblechrom(bunchhit, irev, i0, stopLinearTrackingNo, Apr, Zpr, wecolli, betgam, nonlinflag, scaleorbit, attr1, idpart, idelt, outcoord, plotflag, xco, yco, outputpath, RFflag, indication);
770
771	if (bunchhit.size() == 0) {
772	cout << "No more particles!!" << endl;
773	#if defined(FLUKA)
774	if (RunWithFluka == 1) {
775
776	cout << "End of the connection." << endl;
777
778	int n;
779	flukaio_message_t msg;
780	n = flukaio_send_eoc(lat.conn);
781	if (n < 0) {
782	die("Error occurred sending");
783	}
784
785	n = flukaio_wait_message(lat.conn, &msg);
786	if (n == -1) {
787	die("Server timeout when waiting end of computation");
788	}
789	if (msg.type != N_END) {
790	die("Unexpected message received");
791	}
792
793	flukaio_disconnect(lat.conn);
794	lat.conn = NULL;
795	}
796	#endif
797
798	return;
799	}
800
801
802
803	//tracking the particles to the end of the first revolution
804	cout << "THIRD TRACKING, END OF THE FIRST REVOLUTION." << endl;
805
806	int im;
807
808	im = lat.reseau.size() - 1;
809	irev = 1;
810	nonlinflag = 1;//this tracking is non linear
811
812	lat.trackensemblechrom(bunchhit, irev, stopLinearTrackingNo, im, Apr, Zpr, wecolli, betgam, nonlinflag, scaleorbit, attr1, idpart, idelt, outcoord, plotflag, xco, yco, outputpath, RFflag, indication);
813
814	for (int i(0); i < lat.hits.size(); ++i) {
815	hitso.push_back(lat.hits[i]);//vector of s values where a particle hits the aperture
816	}
817
818	for (int i(0); i < lat.Aphit.size(); ++i) {
819	Aphito.push_back(lat.Aphit[i]);//the mass number of the particles when they hit the aperture
820	}
821
822	for (int i(0); i < lat.Zphit.size(); ++i) {
823	Zphito.push_back(lat.Zphit[i]);//the charge of the particles when they hit the aperture
824	}
825
826
827	for (int i(0); i < lat.nhitcolli.size(); ++i) {
828	nhitcollio.push_back(lat.nhitcolli[i]);//nhitcolli(icop) gives the number of particles that are absorbed at collimator number 'icop'
829	}
830
831	for (int i(0); i < lat.nhitspoiler.size(); ++i) {
832	nhitspoilero.push_back(lat.nhitspoiler[i]);//same as nhitcollio, but for ion spoilers
833	}
834
835
836	vector <int> asumrem, asumhitcolli, asumhitspoiler, asumhits;
837
838
839	asumhitcolli.push_back(0);//how many particles that have been absorbed in a collimator, as a function of the number of revolutions
840	asumhits.push_back(0);//number of particles that have hitted the aperture, as a function of the number of revolutions
841	asumrem.push_back(particlehit);//how many particles there are left in the beam after each revolution
842	asumhits.push_back(lat.hits.size());
843	lat.hits.clear();
844
845	int count1(0);
846	int count(0);
847
848	for (int i(0); i < nhitcollio.size(); ++i) {
849	count = count + nhitcollio[i];
850	}
851
852	count1 = count;
853
854	asumhitcolli.push_back(count);
855
856	lat.nhitcolli.clear();
857
858	count = 0;
859
860	for (int i(0); i < lat.nhitspoiler.size(); ++i) {
861	count = count + lat.nhitspoiler[i];
862	}
863
864	count1 = count1 + count;
865
866	asumhitspoiler.push_back(count);//how many particles that have been absorbed by an ion spoiner, as a function of the number of revolutions
867	lat.nhitspoiler.clear();
868	asumrem.push_back(bunchhit.size());
869
870	if (bunchhit.size() == 0) { //if we have no more particle, we prepare a final output
871	cout << "No more particles!!" << endl;
872	string outputfile;
873	string inputtemp;
874	inputtemp = inputfile.substr(10, inputfile.size() - 1);
875	outputfile = outputpath + "/output.txt";
876
877	writeoutput(outputfile, emix, emiy, Apr, Zpr, nparti, beamflag, PlossPb, i0, im, npart, kbunch, asumrem, asumhitcolli, asumhits, taubeam, hitso, nhitcollio, Aphito, Zphito, i0, stopLinearTrackingNo, xco, yco, plotflag, outputpath, inputpath);
878
879	#if defined(FLUKA)
880	if (RunWithFluka == 1) {
881
882	cout << "End of the connection." << endl;
883
884	int n;
885	flukaio_message_t msg;
886	n = flukaio_send_eoc(lat.conn);
887	if (n < 0) {
888	die("Error occurred sending");
889	}
890
891	n = flukaio_wait_message(lat.conn, &msg);
892	if (n == -1) {
893	die("Server timeout when waiting end of computation");
894	}
895	if (msg.type != N_END) {
896	die("Unexpected message received");
897	}
898
899	flukaio_disconnect(lat.conn);
900	lat.conn = NULL;
901	}
902	#endif
903
904	return;
905	}
906
907	if (asumrem[asumrem.size() - 1] == 0) {
908
909	cout << "No more particle!!" << endl;
910
911	string outputfile;
912	string inputtemp;
913	inputtemp = inputfile.substr(10, inputfile.size() - 1);
914	outputfile = outputpath + "/output.txt";
915
916	writeoutput(outputfile, emix, emiy, Apr, Zpr, nparti, beamflag, PlossPb, i0, im, npart, kbunch, asumrem, asumhitcolli, asumhits, taubeam, hitso, nhitcollio, Aphito, Zphito, i0, stopLinearTrackingNo, xco, yco, plotflag, outputpath, inputpath);
917
918	#if defined(FLUKA)
919	if (RunWithFluka == 1) {
920
921	cout << "End of the connection." << endl;
922
923	int n;
924	flukaio_message_t msg;
925	n = flukaio_send_eoc(lat.conn);
926	if (n < 0) {
927	die("Error occurred sending");
928	}
929
930	n = flukaio_wait_message(lat.conn, &msg);
931	if (n == -1) {
932	die("Server timeout when waiting end of computation");
933	}
934	if (msg.type != N_END) {
935	die("Unexpected message received");
936	}
937
938	flukaio_disconnect(lat.conn);
939	lat.conn = NULL;
940	}
941	#endif
942
943	return;
944	}
945
946
947
948	lat.nhitspoiler.clear();
949
950
951	//tracking the particles through the machine
952	cout << "FOURTH TRACKING, TRACKING THROUGHT THE MACHINE." << endl;
953
954
955	for (int i(1); i < nrev2; ++i) {
956
957	cout << i + 1 << "th revolution." << endl;
958
959
960	bool empty;
961	empty = true;
962
963	for (int j(0); j < bunchhit.size(); ++j) {
964	if (bunchhit[j].inabs == 1) {
965	empty = false;
966	}
967	}
968
969	if (empty == true) {
970
971	cout << "No more particle! The tracking through the machine ends." << endl;
972
973	string outputfile;
974	string inputtemp;
975	inputtemp = inputfile.substr(10, inputfile.size() - 1);
976	outputfile = outputpath + "/output.txt";
977
978
979	writeoutput(outputfile, emix, emiy, Apr, Zpr, nparti, beamflag, PlossPb, i0, im, npart, kbunch, asumrem, asumhitcolli, asumhits, taubeam, hitso, nhitcollio, Aphito, Zphito, i0, stopLinearTrackingNo, xco, yco, plotflag, outputpath, inputpath);
980
981	#if defined(FLUKA)
982	if (RunWithFluka == 1) {
983
984	cout << "End of the connection." << endl;
985
986	int n;
987	flukaio_message_t msg;
988	n = flukaio_send_eoc(lat.conn);
989	if (n < 0) {
990	die("Error occurred sending");
991	}
992
993	n = flukaio_wait_message(lat.conn, &msg);
994	if (n == -1) {
995	die("Server timeout when waiting end of computation");
996	}
997	if (msg.type != N_END) {
998	die("Unexpected message received");
999	}
1000
1001	flukaio_disconnect(lat.conn);
1002	lat.conn = NULL;
1003	}
1004	#endif
1005
1006	return;
1007	}
1008
1009	lat.trackensemblechrom(bunchhit, i, i0, im, Apr, Zpr, wecolli, betgam, nonlinflag, scaleorbit, attr1, idpart, idelt, outcoord, plotflag, xco, yco, outputpath, RFflag, indication);
1010
1011
1012	for (int m(0); m < lat.hits.size(); ++m) {
1013	hitso.push_back(lat.hits[m]);
1014	}
1015
1016	for (int m(0); m < lat.Aphit.size(); ++m) {
1017	Aphito.push_back(lat.Aphit[m]);
1018	}
1019
1020	for (int m(0); m < lat.Zphit.size(); ++m) {
1021	Zphito.push_back(lat.Zphit[m]);
1022	}
1023
1024	for (int m(0); m < lat.nhitcolli.size(); ++m) {
1025	nhitcollio[m] = nhitcollio[m] + lat.nhitcolli[m];
1026	}
1027
1028	for (int m(0); m < lat.nhitspoiler.size(); ++m) {
1029	nhitspoilero[m] = nhitspoilero[m] + lat.nhitspoiler[m];
1030	}
1031
1032	asumhits.push_back(asumhits[asumhits.size() - 1] + lat.hits.size());
1033	lat.hits.clear();
1034	lat.Aphit.clear();
1035	lat.Zphit.clear();
1036
1037	int count1(0);
1038	int count(0);
1039
1040	for (int l(0); l < nhitcollio.size(); ++l) {
1041	count = count + nhitcollio[l];
1042	}
1043
1044	count1 = count;
1045	asumhitcolli.push_back(count);
1046
1047	lat.nhitcolli.clear();
1048
1049	count = 0;
1050
1051	for (int l(0); l < lat.nhitspoiler.size(); ++l) {
1052	count = count + lat.nhitspoiler[l];
1053	}
1054
1055	count1 = count1 + count;
1056
1057	asumhitspoiler.push_back(count);
1058	asumrem.push_back(bunchhit.size());
1059	lat.nhitspoiler.clear();
1060
1061	if (bunchhit.size() == 0) {
1062	cout << "No more particles!!" << endl;
1063
1064	string outputfile;
1065	string inputtemp;
1066	inputtemp = inputfile.substr(10, inputfile.size() - 1);
1067	outputfile = outputpath + "/output.txt";
1068
1069	writeoutput(outputfile, emix, emiy, Apr, Zpr, nparti, beamflag, PlossPb, i0, im, npart, kbunch, asumrem, asumhitcolli, asumhits, taubeam, hitso, nhitcollio, Aphito, Zphito, i0, im, xco, yco, plotflag, outputpath, inputpath);
1070
1071	#if defined(FLUKA)
1072	if (RunWithFluka == 1) {
1073
1074	cout << "End of the connction." << endl;
1075
1076	int n;
1077	flukaio_message_t msg;
1078	n = flukaio_send_eoc(lat.conn);
1079	if (n < 0) {
1080	die("Error occurred sending");
1081	}
1082
1083	n = flukaio_wait_message(lat.conn, &msg);
1084	if (n == -1) {
1085	die("Server timeout when waiting end of computation");
1086	}
1087	if (msg.type != N_END) {
1088	die("Unexpected message received");
1089	}
1090
1091	flukaio_disconnect(lat.conn);
1092	lat.conn = NULL;
1093	}
1094	#endif
1095
1096	return;
1097	}
1098
1099	if (asumrem[i] == 0) {
1100	string outputfile;
1101	string inputtemp;
1102	inputtemp = inputfile.substr(10, inputfile.size() - 1);
1103	outputfile = outputpath + "/output.txt";
1104
1105	writeoutput(outputfile, emix, emiy, Apr, Zpr, nparti, beamflag, PlossPb, i0, im, npart, kbunch, asumrem, asumhitcolli, asumhits, taubeam, hitso, nhitcollio, Aphito, Zphito, i0, im, xco, yco, plotflag, outputpath, inputpath);
1106
1107	#if defined(FLUKA)
1108	if (RunWithFluka == 1) {
1109
1110	cout << "End of the connection." << endl;
1111
1112	int n;
1113	flukaio_message_t msg;
1114	n = flukaio_send_eoc(lat.conn);
1115	if (n < 0) {
1116	die("Error occurred sending");
1117	}
1118
1119	n = flukaio_wait_message(lat.conn, &msg);
1120	if (n == -1) {
1121	die("Server timeout when waiting end of computation");
1122	}
1123	if (msg.type != N_END) {
1124	die("Unexpected message received");
1125	}
1126
1127	flukaio_disconnect(lat.conn);
1128	lat.conn = NULL;
1129	}
1130	#endif
1131
1132	return;
1133	}
1134
1135	}
1136
1137	string outputfile;
1138	string inputtemp;
1139	inputtemp = inputfile.substr(10, inputfile.size() - 1);
1140	outputfile = outputpath + "/output.txt";
1141
1142	writeoutput(outputfile, emix, emiy, Apr, Zpr, nparti, beamflag, PlossPb, i0, im, npart, kbunch, asumrem, asumhitcolli, asumhits, taubeam, hitso, nhitcollio, Aphito, Zphito, i0, im, xco, yco, plotflag, outputpath, inputpath);
1143
1144	#if defined(FLUKA)
1145	if (RunWithFluka == 1) {
1146
1147	cout << "End of the connection." << endl;
1148
1149	int n;
1150	flukaio_message_t msg;
1151	n = flukaio_send_eoc(lat.conn);
1152	if (n < 0) {
1153	die("Error occurred sending");
1154	}
1155
1156	n = flukaio_wait_message(lat.conn, &msg);
1157	if (n == -1) {
1158	die("Server timeout when waiting end of computation");
1159	}
1160	if (msg.type != N_END) {
1161	die("Unexpected message received");
1162	}
1163
1164	flukaio_disconnect(lat.conn);
1165	lat.conn = NULL;
1166	}
1167	#endif
1168
1169	}
1170
1171
1172	void Simulation::extractOpticsParameters(string opticsFile, const int& beamflag)
1173	{
1174
1175
1176	int JohnJaperFlag;
1177
1178	defineParameters(opticsFile, beamflag, JohnJaperFlag);
1179
1180	bool apertureIsZero(false);
1181	int temp1(0);
1182	int temp2(0);
1183	int temp3;
1184
1185	if (JohnJaperFlag == 1) { //JohnJaperflag=1 if we have information about the aperture
1186
1187	for (int w(0); w < nbre_elts; ++w) {
1188	if (lat.reseau[w]->APERTYPE == "RECTELLIPSE") {
1189	++temp1;
1190	if ((lat.reseau[w]->APER_1 == 0) && (lat.reseau[w]->APER_2 == 0) && (lat.reseau[w]->APER_3 == 0) && (lat.reseau[w]->APER_4 == 0)) {
1191	++temp2;
1192	} else {
1193	break;
1194	}
1195	} else {
1196	break;
1197	}
1198	}
1199
1200	//If all apertures are RECTELLIPSES:
1201	//If an apeture variable is 0, it takes the value of the preceeding aperture variable.
1202	//If all the apertures are 0, we define them to be 1 instead.
1203
1204	if (temp1 == nbre_elts) {
1205	if (temp2 == nbre_elts) {
1206	for (int w(0); w < nbre_elts; ++w) {
1207	lat.reseau[w]->APER_1 = 1;
1208	lat.reseau[w]->APER_2 = 1;
1209	lat.reseau[w]->APER_3 = 1;
1210	lat.reseau[w]->APER_4 = 1;
1211	}
1212	}
1213	//if some apertures of the first elements are 0, we define them equal to those of the element nearest of the end of the lattice with apertures not equal to 0
1214	if (lat.reseau[0]->APER_1 == 0) {
1215	apertureIsZero = true;
1216	temp3 = nbre_elts;
1217	while (lat.reseau[temp3]->APER_1 == 0) {
1218	temp3 = temp3 - 1;
1219	}
1220	lat.reseau[0]->APER_1 = lat.reseau[temp3]->APER_1;
1221	}
1222	if (lat.reseau[0]->APER_2 == 0) {
1223	apertureIsZero = true;
1224	temp3 = nbre_elts;
1225	while (lat.reseau[temp3]->APER_2 == 0) {
1226	temp3 = temp3 - 1;
1227	}
1228	lat.reseau[0]->APER_2 = lat.reseau[temp3]->APER_2;
1229	}
1230	if (lat.reseau[0]->APER_3 == 0) {
1231	apertureIsZero = true;
1232	temp3 = nbre_elts;
1233	while (lat.reseau[temp3]->APER_3 == 0) {
1234	temp3 = temp3 - 1;
1235	}
1236	lat.reseau[0]->APER_3 = lat.reseau[temp3]->APER_3;
1237	}
1238	if (lat.reseau[0]->APER_4 == 0) {
1239	apertureIsZero = true;
1240	temp3 = nbre_elts;
1241	while (lat.reseau[temp3]->APER_4 == 0) {
1242	temp3 = temp3 - 1;
1243	}
1244	lat.reseau[0]->APER_4 = lat.reseau[temp3]->APER_4;
1245	}
1246
1247	//if some other apertures are equal to 0, we define them equal to the preceding ones
1248	for (int p(1); p < nbre_elts; ++p) {
1249	if (lat.reseau[p]->APER_1 == 0) {
1250	apertureIsZero = true;
1251	lat.reseau[p]->APER_1 = lat.reseau[p - 1]->APER_1;
1252	}
1253	if (lat.reseau[p]->APER_2 == 0) {
1254	apertureIsZero = true;
1255	lat.reseau[p]->APER_2 = lat.reseau[p - 1]->APER_2;
1256	}
1257	if (lat.reseau[p]->APER_3 == 0) {
1258	apertureIsZero = true;
1259	lat.reseau[p]->APER_3 = lat.reseau[p - 1]->APER_3;
1260	}
1261	if (lat.reseau[p]->APER_4 == 0) {
1262	apertureIsZero = true;
1263	lat.reseau[p]->APER_4 = lat.reseau[p - 1]->APER_4;
1264	}
1265	}
1266	}
1267
1268	//if the apertures are not only RECTANGLE:
1269	//if the aperture type is known, we just take the preceding apertures if we don't know the current ones (not the ideal solution if the aperture type changes);
1270	//if the aperture type is NONE, we also take the aperture type of the preceding element
1271	for (int r(0); r < nbre_elts; ++r) {
1272	if (lat.reseau[r]->APERTYPE != "NONE") {
1273	if ((lat.reseau[r]->APER_1 == 0) && (lat.reseau[r]->APER_2 == 0) && (lat.reseau[r]->APER_3 == 0) && (lat.reseau[r]->APER_4 == 0)) {
1274	apertureIsZero = true;
1275	lat.reseau[r]->APER_1 = 1;;
1276	lat.reseau[r]->APER_2 = 1;
1277	lat.reseau[r]->APER_3 = 1;
1278	lat.reseau[r]->APER_4 = 1;
1279	}
1280	} else if (lat.reseau[r]->APERTYPE == "NONE") {
1281	lat.reseau[r]->APER_1 = lat.reseau[r - 1]->APER_1;
1282	lat.reseau[r]->APER_2 = lat.reseau[r - 1]->APER_2;
1283	lat.reseau[r]->APER_3 = lat.reseau[r - 1]->APER_3;
1284	lat.reseau[r]->APER_4 = lat.reseau[r - 1]->APER_4;
1285	lat.reseau[r]->APERTYPE = lat.reseau[r - 1]->APERTYPE;
1286	}
1287	}
1288
1289	if (apertureIsZero == true) {
1290	cout << "At least one of the apertures in the file is 0." << endl;
1291	}
1292
1293	//we compute the final apertures
1294	for (int i(0); i < lat.reseau.size(); ++i) {
1295	if (lat.reseau[i]->APERTYPE == "CIRCLE") {
1296	lat.reseau[i]->aperx = lat.reseau[i]->APER_1;
1297	lat.reseau[i]->apery = lat.reseau[i]->APER_1;
1298	} else if (lat.reseau[i]->APERTYPE == "ELLIPSE") {
1299	lat.reseau[i]->aperx = lat.reseau[i]->APER_1;
1300	lat.reseau[i]->apery = lat.reseau[i]->APER_2;
1301	} else if (lat.reseau[i]->APERTYPE == "RECTANGLE") {
1302	lat.reseau[i]->aperx = lat.reseau[i]->APER_1;
1303	lat.reseau[i]->apery = lat.reseau[i]->APER_2;
1304	} else if (lat.reseau[i]->APERTYPE == "LHCSCREEN") {
1305	if (lat.reseau[i]->APER_1 == 0) {
1306	lat.reseau[i]->aperx = lat.reseau[i]->APER_3;
1307	lat.reseau[i]->apery = lat.reseau[i]->APER_2;
1308	} else if (lat.reseau[i]->APER_2 == 0) {
1309	lat.reseau[i]->aperx = lat.reseau[i]->APER_1;
1310	lat.reseau[i]->apery = lat.reseau[i]->APER_3;
1311	} else {
1312	cerr << "ERROR READMAD: invalid LHCSCREEN aperture definition!" << endl;
1313	}
1314	} else if (lat.reseau[i]->APERTYPE == "RECTELLIPSE") {
1315	if (lat.reseau[i]->APER_1 == 0) {
1316	lat.reseau[i]->APER_1 = 1;
1317	} else if (lat.reseau[i]->APER_2 == 0) {
1318	lat.reseau[i]->APER_2 = 1;
1319	} else if (lat.reseau[i]->APER_3 == 0) {
1320	lat.reseau[i]->APER_3 = 1;
1321	} else if (lat.reseau[i]->APER_4 == 0) {
1322	lat.reseau[i]->APER_4 = 1;
1323	}
1324	if (lat.reseau[i]->APER_1 > lat.reseau[i]->APER_3) {
1325	lat.reseau[i]->aperx = lat.reseau[i]->APER_3;
1326	} else {
1327	lat.reseau[i]->aperx = lat.reseau[i]->APER_1;
1328	}
1329	if (lat.reseau[i]->APER_2 > lat.reseau[i]->APER_4) {
1330	lat.reseau[i]->apery = lat.reseau[i]->APER_4;
1331	} else {
1332	lat.reseau[i]->apery = lat.reseau[i]->APER_2;
1333	}
1334	}
1335	if ((lat.reseau[i]->NAME[0] == 'T') && (lat.reseau[i]->NAME[1] == 'C') && (lat.reseau[i]->NAME[2] == 'P')) {
1336	lat.reseau[i]->aperx = 1;
1337	lat.reseau[i]->apery = 1;
1338	}
1339	if ((lat.reseau[i]->NAME[0] == 'T') && (lat.reseau[i]->NAME[1] == 'C') && (lat.reseau[i]->NAME[2] == 'S')) {
1340	lat.reseau[i]->aperx = 1;
1341	lat.reseau[i]->apery = 1;
1342	}
1343	if ((lat.reseau[i]->NAME[0] == 'T') && (lat.reseau[i]->NAME[1] == 'C') && (lat.reseau[i]->NAME[2] == 'T')) {
1344	lat.reseau[i]->aperx = 1;
1345	lat.reseau[i]->apery = 1;
1346	}
1347	if ((lat.reseau[i]->NAME[0] == 'T') && (lat.reseau[i]->NAME[1] == 'C') && (lat.reseau[i]->NAME[2] == 'L') && (lat.reseau[i]->NAME[3] == 'A')) {
1348	lat.reseau[i]->aperx = 1;
1349	lat.reseau[i]->apery = 1;
1350	}
1351	if ((lat.reseau[i]->NAME[0] == 'T') && (lat.reseau[i]->NAME[1] == 'C') && (lat.reseau[i]->NAME[2] == 'R') && (lat.reseau[i]->NAME[3] == 'Y') && (lat.reseau[i]->NAME[4] == 'O')) {
1352	lat.reseau[i]->aperx = 1;
1353	lat.reseau[i]->apery = 1;
1354	}
1355	}
1356
1357	//we take into account that the elements have a length not equal to 0 to determine the real apertures, using interpolation
1358
1359	vector <double> Lelem, sh1, aperx2, apery2;
1360	vector <int> indice;
1361
1362	Lelem.push_back(0);
1363
1364	for (int k(0); k < lat.reseau.size() - 1; ++k) {
1365	Lelem.push_back(lat.reseau[k + 1]->S - lat.reseau[k]->S);
1366	}
1367
1368	for (int i(0); i < Lelem.size(); ++i) {
1369	if (Lelem[i] > 0) {
1370
1371	sh1.push_back(lat.reseau[i - 1]->S + 1e10 * pow(2., -52));
1372	lat.aperx1.push_back(lat.reseau[i]->aperx);
1373	lat.apery1.push_back(lat.reseau[i]->apery);
1374	}
1375	}
1376	for (int j(0); j < lat.reseau.size(); ++j) {
1377	lat.sh.push_back(lat.reseau[j]->S);
1378	}
1379	for (int j(0); j < sh1.size(); ++j) {
1380	lat.sh.push_back(sh1[j]);
1381	}
1382
1383	for (int j(0); j < lat.reseau.size(); ++j) {
1384	aperx2.push_back(lat.reseau[j]->aperx);
1385	}
1386	for (int j(0); j < sh1.size(); ++j) {
1387	aperx2.push_back(lat.aperx1[j]);
1388	}
1389
1390	for (int j(0); j < lat.reseau.size(); ++j) {
1391	apery2.push_back(lat.reseau[j]->apery);
1392	}
1393	for (int j(0); j < sh1.size(); ++j) {
1394	apery2.push_back(lat.apery1[j]);
1395	}
1396
1397	sh1.clear();
1398
1399	for (int j(0); j < lat.sh.size(); ++j) {
1400	sh1.push_back(lat.sh[j]);
1401	}
1402	sort(lat.sh.begin(), lat.sh.end());
1403
1404	for (int j(0); j < lat.sh.size(); ++j) {
1405	int m(0);
1406	for (int k(0); k < sh1.size(); ++k) {
1407	if (lat.sh[j] == sh1[k]) {
1408	indice.push_back(m);
1409	sh1[k] = -10;
1410	break;
1411	} else {
1412	++m;
1413	}
1414	}
1415	}
1416
1417	lat.aperx1.clear();
1418	lat.apery1.clear();
1419
1420	for (int k(0); k < indice.size(); ++k) {
1421	for (int j(0); j < aperx2.size(); ++j) {
1422	if (j == indice[k]) {
1423	lat.aperx1.push_back(aperx2[j]);
1424	lat.apery1.push_back(apery2[j]);
1425	break;
1426	}
1427	}
1428	}
1429
1430	aperx2.clear();
1431	apery2.clear();
1432
1433
1434	for (int j(0); j < lat.aperx1.size(); ++j) {
1435	if (lat.aperx1[j] < 1) {
1436	aperx2.push_back(lat.aperx1[j]);
1437	}
1438	}
1439
1440	double mahx;
1441
1442	mahx = *max_element(aperx2.begin(), aperx2.end());
1443
1444	for (int j(0); j < lat.aperx1.size(); ++j) {
1445	if (lat.aperx1[j] == 1) {
1446	lat.aperx1[j] = mahx;
1447	}
1448	}
1449
1450
1451	for (int j(0); j < lat.apery1.size(); ++j) {
1452	if (lat.apery1[j] < 1) {
1453	apery2.push_back(lat.apery1[j]);
1454	}
1455	}
1456
1457	double mahy;
1458
1459	mahy = *max_element(apery2.begin(), apery2.end());
1460
1461	for (int j(0); j < lat.apery1.size(); ++j) {
1462	if (lat.apery1[j] == 1) {
1463	lat.apery1[j] = mahy;
1464	}
1465	}
1466
1467	for (int j(0); j < lat.sh.size() - 1; ++j) {
1468	if (lat.sh[j + 1] - lat.sh[j] == 0) {
1469	lat.sh[j] = lat.sh[j] - 0.0001;
1470	}
1471	}
1472
1473	for (int j(0); j < lat.reseau.size(); ++j) {
1474	int m(0);
1475	for (int k(0); k < lat.sh.size(); ++k) {
1476	if ((lat.reseau[j]->S < lat.sh[k]) && (lat.reseau[j]->S >= lat.sh[k - 1])) {
1477	lat.reseau[j]->aperx = interp1(lat.sh[k - 1], lat.aperx1[k - 1], lat.sh[k], lat.aperx1[k], lat.reseau[j]->S);
1478	lat.reseau[j]->apery = interp1(lat.sh[k - 1], lat.apery1[k - 1], lat.sh[k], lat.apery1[k], lat.reseau[j]->S);
1479	}
1480	}
1481	}
1482
1483	indice.clear();
1484	} else {//if JohnJaperflag=0, we don't have apertures information, so we throw an error
1485
1486	cout << "We don't have apertures informations!" << endl;
1487
1488	throw 1;
1489
1490	}
1491
1492	vector <int> indice;
1493
1494	for (int j(0); j < lat.reseau.size(); ++j) {
1495	if (lat.reseau[j]->KEYWORD == "HKICKER") {
1496	indice.push_back(j);
1497	}
1498	}
1499
1500	for (int j(0); j < lat.reseau.size(); ++j) {
1501	if (lat.reseau[j]->KEYWORD == "KICKER") {
1502	indice.push_back(j);
1503	}
1504	}
1505
1506	double hk;
1507
1508	for (int k(0); k < indice.size(); ++k) {
1509
1510	hk = lat.reseau[indice[k]]->PXC - lat.reseau[indice[k] - 1]->PXC;
1511	if (hk != 0) {
1512	lat.xcormag.push_back(hk);
1513	lat.ixcormag.push_back(indice[k]);
1514	}
1515	}
1516
1517	indice.clear();
1518
1519	for (int j(0); j < lat.reseau.size(); ++j) {
1520	if (lat.reseau[j]->KEYWORD == "VKICKER") {
1521	indice.push_back(j);
1522	}
1523	}
1524
1525	for (int j(0); j < lat.reseau.size(); ++j) {
1526	if (lat.reseau[j]->KEYWORD == "KICKER") {
1527	indice.push_back(j);
1528	}
1529	}
1530
1531	double vk;
1532
1533	for (int k(0); k < indice.size(); ++k) {
1534
1535	vk = lat.reseau[indice[k]]->PYC - lat.reseau[indice[k] - 1]->PYC;
1536	if (vk != 0) {
1537	lat.ycormag.push_back(vk);
1538	lat.iycormag.push_back(indice[k]);
1539	}
1540	}
1541
1542	}
1543
1544	void Simulation::defineParameters(string opticsFile, const int& beamflag, int& JohnJaperFlag)
1545	{
1546
1547	ifstream entree;
1548	string APERTYPE, name, KEYWORD, PARENT;
1549	entree.open(opticsFile.c_str());
1550
1551	if (entree.fail()) {
1552	cerr << "Warning: problem with the optics file " << opticsFile << "!!" << endl;
1553	} else {
1554
1555	string name, temp;
1556	string strtemp, strtemp1;
1557
1558	getline(entree, strtemp, ',');
1559	while ((strtemp != "NAME") && (strtemp != "ALFX") && (strtemp != "S") && (strtemp != "L") && (strtemp != "K0L") && (strtemp != "K0SL") && (strtemp != "K1L") && (strtemp != "K1SL") && (strtemp != "K2L") && (strtemp != "K2SL") && (strtemp != "XC") && (strtemp != "PXC") && (strtemp != "KEYWORD") && (strtemp != "PARENT") && (strtemp != "YC") && (strtemp != "PYC") && (strtemp != "TC") && (strtemp != "PTC") && (strtemp != "BETX") && (strtemp != "MUX") && (strtemp != "DX") && (strtemp != "DPX") && (strtemp != "BETY") && (strtemp != "ALFY") && (strtemp != "MUY") && (strtemp != "DY") && (strtemp != "DPY") && (strtemp != "APERTYPE") && (strtemp != "APER_1") && (strtemp != "APER_2") && (strtemp != "APER_3") && (strtemp != "APER_4")) {
1560	getline(entree, strtemp1);
1561	getline(entree, strtemp, ',');
1562	}
1563
1564	double S, L, K0L, K0SL, K1L, K1SL, K2L, K2SL, XC, PXC, YC, PYC, TC, PTC, BETX, ALFX, MUX, DX, DPX, BETY, ALFY, MUY, DY, DPY, APER_1, APER_2, APER_3, APER_4;
1565
1566	int taille(0);
1567
1568	vector <string> header;
1569	vector <string> header_base;
1570	vector <string> vectemp;
1571	vector <int> order;
1572
1573	//header_base is: NAME,KEYWORD,PARENT,S,L,K0L,K0SL,K1L,K1SL,K2L,K2SL,XC,PXC,YC,PYC,TC,PTC,BETX,ALFX,MUX,DX,DPX,BETY,ALFY,MUY,DY,DPY,APERTYPE,APER_1,APER_2,APER_3,APER_4
1574	//we can have any order of the headers in the opticsfile
1575	header_base.push_back("NAME");
1576	header_base.push_back("KEYWORD");
1577	header_base.push_back("PARENT");
1578	header_base.push_back("S");
1579	header_base.push_back("L");
1580	header_base.push_back("K0L");
1581	header_base.push_back("K0SL");
1582	header_base.push_back("K1L");
1583	header_base.push_back("K1SL");
1584	header_base.push_back("K2L");
1585	header_base.push_back("K2SL");
1586	header_base.push_back("XC");
1587	header_base.push_back("PXC");
1588	header_base.push_back("YC");
1589	header_base.push_back("PYC");
1590	header_base.push_back("TC");
1591	header_base.push_back("PTC");
1592	header_base.push_back("BETX");
1593	header_base.push_back("ALFX");
1594	header_base.push_back("MUX");
1595	header_base.push_back("DX");
1596	header_base.push_back("DPX");
1597	header_base.push_back("BETY");
1598	header_base.push_back("ALFY");
1599	header_base.push_back("MUY");
1600	header_base.push_back("DY");
1601	header_base.push_back("DPY");
1602	header_base.push_back("APERTYPE");
1603	header_base.push_back("APER_1");
1604	header_base.push_back("APER_2");
1605	header_base.push_back("APER_3");
1606	header_base.push_back("APER_4");
1607
1608	header.push_back(strtemp);
1609	for (int n(0); n < 30; ++n) {
1610	getline(entree, strtemp, ',');
1611	header.push_back(strtemp);
1612	}
1613	getline(entree, strtemp);
1614	header.push_back(strtemp);
1615
1616	for (int t(0); t < header.size(); ++t) {
1617	for (int q(0); q < header_base.size(); ++q) {
1618	if (header[t] == header_base[q]) {
1619	order.push_back(q);
1620	break;
1621	}
1622	}
1623	}
1624
1625	while (!entree.eof()) {
1626	++taille;
1627	for (int n(0); n < 31; ++n) {
1628	getline(entree, strtemp, ',');
1629	vectemp.push_back(strtemp);
1630	}
1631	getline(entree, strtemp);
1632	vectemp.push_back(strtemp);
1633
1634	for (int h(0); h < order.size(); ++h) {
1635
1636	if (order[h] == 0) {
1637	name = vectemp[h];
1638	} else if (order[h] == 1) {
1639	KEYWORD = vectemp[h];
1640	} else if (order[h] == 2) {
1641	PARENT = vectemp[h];
1642	} else if (order[h] == 3) {
1643	S = atof(vectemp[h].c_str());
1644	} else if (order[h] == 4) {
1645	L = atof(vectemp[h].c_str());
1646	} else if (order[h] == 5) {
1647	K0L = atof(vectemp[h].c_str());
1648	} else if (order[h] == 6) {
1649	K0SL = atof(vectemp[h].c_str());
1650	} else if (order[h] == 7) {
1651	K1L = atof(vectemp[h].c_str());
1652	} else if (order[h] == 8) {
1653	K1SL = atof(vectemp[h].c_str());
1654	} else if (order[h] == 9) {
1655	K2L = atof(vectemp[h].c_str());
1656	} else if (order[h] == 10) {
1657	K2SL = atof(vectemp[h].c_str());
1658	} else if (order[h] == 11) {
1659	XC = atof(vectemp[h].c_str());
1660	} else if (order[h] == 12) {
1661	PXC = atof(vectemp[h].c_str());
1662	} else if (order[h] == 13) {
1663	YC = atof(vectemp[h].c_str());
1664	} else if (order[h] == 14) {
1665	PYC = atof(vectemp[h].c_str());
1666	} else if (order[h] == 15) {
1667	TC = atof(vectemp[h].c_str());
1668	} else if (order[h] == 16) {
1669	PTC = atof(vectemp[h].c_str());
1670	} else if (order[h] == 17) {
1671	BETX = atof(vectemp[h].c_str());
1672	} else if (order[h] == 18) {
1673	ALFX = atof(vectemp[h].c_str());
1674	} else if (order[h] == 19) {
1675	MUX = atof(vectemp[h].c_str());
1676	} else if (order[h] == 20) {
1677	DX = atof(vectemp[h].c_str());
1678	} else if (order[h] == 21) {
1679	DPX = atof(vectemp[h].c_str());
1680	} else if (order[h] == 22) {
1681	BETY = atof(vectemp[h].c_str());
1682	} else if (order[h] == 23) {
1683	ALFY = atof(vectemp[h].c_str());
1684	} else if (order[h] == 24) {
1685	MUY = atof(vectemp[h].c_str());
1686	} else if (order[h] == 25) {
1687	DY = atof(vectemp[h].c_str());
1688	} else if (order[h] == 26) {
1689	DPY = atof(vectemp[h].c_str());
1690	} else if (order[h] == 27) {
1691	APERTYPE = vectemp[h];
1692	} else if (order[h] == 28) {
1693	APER_1 = atof(vectemp[h].c_str());
1694	} else if (order[h] == 29) {
1695	APER_2 = atof(vectemp[h].c_str());
1696	} else if (order[h] == 30) {
1697	APER_3 = atof(vectemp[h].c_str());
1698	} else if (order[h] == 31) {
1699	APER_4 = atof(vectemp[h].c_str());
1700	}
1701	}
1702
1703	Element elt(ALFX, ALFY, APER_1, APER_2, APER_3, APER_4, APERTYPE, BETX, BETY, DPX, DPY, DX, DY, KEYWORD, L, MUX, MUY, name, PTC, PXC, PYC, S, TC, XC, YC, K0L, K0SL, K1L, K1SL, K2L, K2SL, PARENT);
1704
1705	lat.addElement(new Element(elt));
1706
1707	vectemp.clear();
1708
1709	if ((name.substr(0, 3) == "END") \|\| (name.substr(name.size() - 3, 3) == "END")) {
1710	break;
1711	}
1712
1713	}
1714	this ->nbre_elts = taille;
1715	}
1716	entree.close();
1717
1718	if (beamflag == 2) {
1719
1720	Lattice lattemp(lat);
1721
1722	for (int i(0); i < nbre_elts; ++i) {
1723
1724	lat.reseau[nbre_elts - 1 - i] = lattemp.reseau[i];
1725	}
1726	lattemp = lat;
1727	lat.reseau[0] = lattemp.reseau[nbre_elts - 1];
1728
1729	for (int j(1); j < nbre_elts; ++j) {
1730	lat.reseau[j] = lattemp.reseau[j - 1];
1731	}
1732
1733	}
1734
1735	JohnJaperFlag = 1;//for the moment, ICOSIM++ cannot run without aperture parameters, so JohnJaperflag is fixed to 1; the variable is still here to facilitate further developpment
1736
1737	}
1738
1739
1740
1741	double Simulation::interp1(double x1, double y1, double x2, double y2, double x)
1742	{
1743
1744	double d1(x1 - x2);
1745	double d2(y1 - y2);
1746
1747	double y;
1748
1749	y = (x - x1) * (d2 / d1) + y1;
1750
1751	return y;
1752	}
1753
1754
1755	void Simulation::genbunch(const int& i0, const int& npart, const string& partdistr, const double& r1r2skin, const double& emx, const double& emy, const double& sigdpp, const double& bx, const double& ax, const double& dx, const double& dpx, const double& by, const double& ay, const double& dy, const double& dpy, const double& nsigi)
1756	{
1757
1758	srand(time(0));
1759
1760	Particle p1;
1761
1762	for (int l(0); l < npart; ++l) {
1763
1764	p1.genpartdist(i0, 1, partdistr, r1r2skin, emx, emy, sigdpp, bx, ax, dx, dpx, by, ay, dy, dpy, nsigi);
1765
1766	bunch.push_back(p1);
1767	}
1768
1769	}
1770
1771
1772	void Simulation::affichebunch()
1773	{
1774
1775	for (int i(0); i < bunch.size(); ++i) {
1776
1777	cout << "Particle " << i + 1 << endl;
1778	bunch[i].afficheFirstCoordonnees();
1779	}
1780	}
1781
1782
1783	void Simulation::writeoutput(string outputfile, const double& emix, const double& emiy, const double& Apr, const double& Zpr, const double& nparti, const double& beamflag, const double& PlossPb, const int& i0, const int& im, const int& npart, const int& kbunch, const vector <int>& asumrem, const vector <int>& asumhitcolli, const vector <int>& asumhits, const int& taubeam, vector <double>& hitso, const vector <int>& nhitcollio, vector <double>& Aphito, vector <double>& Zphito, const int& ibeg, const int& iend, vector <vector <double> >& xco, vector <vector <double> >& yco, const string& plotflag, string outputpath, string inputpath)
1784	{
1785
1786	ofstream sortie(outputfile.c_str());
1787
1788	if (sortie.fail()) {
1789	cerr << "Warning: problem with the output file!!" << endl;
1790	} else {
1791
1792	sortie << "Position of TCP, CRYSTAL, TCRYO collimators (ip):";
1793
1794	for (int i(0); i < lat.ip.size(); ++i) {
1795	sortie << lat.ip[i] << ", ";
1796	}
1797
1798	sortie << endl;
1799
1800	sortie << "Position of TCS collimators (is):";
1801
1802	for (int i(0); i < lat.is.size(); ++i) {
1803	sortie << lat.is[i] << ", ";
1804	}
1805
1806	sortie << endl;
1807
1808	sortie << "Position of all the collimators (ips):";
1809
1810	for (int i(0); i < lat.ips.size(); ++i) {
1811	sortie << lat.ips[i] << ", ";
1812	}
1813
1814	sortie << endl;
1815
1816	sortie << "Initial size of the beam (nsig):";
1817
1818	for (int i(0); i < lat.resColli.size(); ++i) {
1819	sortie << lat.resColli[i]->nsig << ", ";
1820	}
1821
1822	sortie << endl;
1823
1824
1825	sortie << "Mass number of the reference particle (Apr): " << Apr << endl;
1826	sortie << "Charge state of the reference particle (Zpr): " << Zpr << endl;
1827	sortie << "Number of particles in the simulation (nparti): " << nparti << endl;
1828	sortie << "Average power loss of the beam (PlossPb): " << PlossPb << endl;
1829
1830	sortie << "Mass for the particles hitting the aperture (Aphit):";
1831
1832	for (int i(0); i < Aphito.size(); ++i) {
1833	sortie << Aphito[i] << ", ";
1834	}
1835
1836	sortie << endl;
1837
1838	sortie << "Charge of the particles hitting the aperture (Zphit):";
1839
1840	for (int i(0); i < Zphito.size(); ++i) {
1841	sortie << Zphito[i] << ", ";
1842	}
1843
1844	sortie << endl;
1845
1846	sortie << "Number of particle staying in the accelerator for each turn (asumrem):";
1847
1848	for (int i(0); i < asumrem.size(); ++i) {
1849	sortie << asumrem[i] << ", ";
1850	}
1851
1852	sortie << endl;
1853
1854	sortie << "Number of particles that have hit a collimator during a preceding turn (asumhitcolli):";
1855
1856	for (int i(0); i < asumhitcolli.size(); ++i) {
1857	sortie << asumhitcolli[i] << ", ";
1858	}
1859
1860	sortie << endl;
1861
1862	sortie << "Number of particles that have hit another element during a preceding turn (asumhits):";
1863
1864	for (int i(0); i < asumhits.size(); ++i) {
1865	sortie << asumhits[i] << ", ";
1866	}
1867
1868	sortie << endl;
1869
1870	sortie << "Number of particles that have hit the collimators, for each collimator (nhitcolli):";
1871
1872	for (int i(0); i < nhitcollio.size(); ++i) {
1873	sortie << nhitcollio[i] << ", ";
1874	}
1875
1876	sortie << endl;
1877
1878
1879	sortie << "Hit position along the accelerator for the particles that have hit another element than a collimator (hits):";
1880
1881	for (int i(0); i < hitso.size(); ++i) {
1882	sortie << hitso[i] << ", ";
1883	}
1884
1885	sortie << endl;
1886
1887
1888	sortie << "Expected lifetime of the beam (taubeam): " << taubeam;
1889
1890	sortie.close();
1891
1892	}
1893
1894	PlotRunSummary(asumrem, asumhitcolli, asumhits, nhitcollio, hitso, PlossPb, outputpath);
1895	PlotLossSpectra(beamflag, nparti, asumrem, Apr, Zpr, hitso, Aphito, Zphito, PlossPb, npart, kbunch, taubeam, nhitcollio, outputpath, inputpath);
1896
1897	if (plotflag == "Yes") {
1898	PlotTrajectory(emix, emiy, ibeg, iend, xco, yco, outputpath);
1899	}
1900
1901	}
1902
1903
1904
1905	void Simulation::PlotRunSummary(const vector <int>& asumrem, const vector <int>& asumhitcolli, const vector <int>& asumhits, const vector <int>& nhitcollio, vector <double>& hitso, const double& PlossPb, string outputpath)
1906	{
1907
1908	//Preparing files for the plotting of the loss map
1909
1910	int NLostTotal(0);
1911	int irevloc, i0loc, imloc;
1912
1913	NLostTotal = asumhitcolli[asumhitcolli.size() - 1] + asumhits[asumhits.size() - 1];
1914
1915	irevloc = asumrem.size() - 1;
1916	i0loc = 0;
1917	imloc = lat.reseau.size() - 1;
1918
1919	sort(hitso.begin(), hitso.end());
1920
1921	vector <double> sep;
1922	vector <double> plosslocal, zlosslocal;
1923
1924	for (double k(lat.reseau[i0loc]->S); k < lat.reseau[imloc]->S; ++k) {
1925	sep.push_back(k);
1926	}
1927
1928	for (int i(0); i < sep.size() - 1; ++i) {
1929	zlosslocal.push_back(sep[i] + (sep[i + 1] - sep[i]) / 2);
1930	}
1931
1932	int count(0);
1933
1934	for (int i(0); i < sep.size() - 1; ++i) {
1935	for (int j(0); j < hitso.size(); ++j) {
1936	if ((hitso[j] >= sep[i]) && (hitso[j] < sep[i + 1])) {
1937	++count;
1938	}
1939	}
1940	plosslocal.push_back(count);
1941	count = 0;
1942	}
1943
1944	ofstream outplot1;
1945	string plot1;
1946	plot1 = outputpath + "/PlotRunSummary1.dat";
1947
1948	outplot1.open(plot1.c_str());
1949
1950	if (outplot1.fail()) {
1951	cerr << "Warning: problem with the file " << plot1 << "!" << endl;
1952	} else {
1953
1954	for (int i(0); i < plosslocal.size(); ++i) {
1955	outplot1 << setw(10) << plosslocal[i]*PlossPb / NLostTotal << setw(10) << zlosslocal[i] / 1000 << endl;
1956	}
1957
1958	outplot1.close();
1959	}
1960
1961	vector <double> sip;
1962	vector <int> pos;
1963
1964	for (int i(0); i < lat.reseau.size(); ++i) {
1965	if (lat.reseau[i]->NAME.substr(0, 2) == "IP") {
1966	sip.push_back(lat.reseau[i]->S);
1967	pos.push_back(i);
1968	}
1969	}
1970
1971	if (sip.size() != 0) {
1972
1973	ofstream outplot2;
1974	string plot2;
1975	plot2 = outputpath + "/PlotRunSummary2.dat";
1976
1977	outplot2.open(plot2.c_str());
1978
1979	if (outplot2.fail()) {
1980	cerr << "Warning: problem with the file " << plot2 << "!" << endl;
1981	} else {
1982
1983	for (int i(0); i < sip.size(); ++i) {
1984	outplot2 << setw(10) << 0.5 << setw(10) << sip[i] / 1000 << setw(10) << lat.reseau[pos[i]]->NAME << endl;
1985	}
1986	outplot2.close();
1987	}
1988
1989	ofstream gnuout;
1990	string gnuscript;
1991	gnuscript = outputpath + "/labelsIP.p";
1992
1993	gnuout.open(gnuscript.c_str());
1994
1995	if (gnuout.fail()) {
1996	cerr << "Warning: problem with the file " << gnuscript << "!" << endl;
1997	} else {
1998
1999	for (int i(0); i < sip.size(); ++i) {
2000	gnuout << "set label " << i + 1 << " " << '"' << lat.reseau[pos[i]]->NAME << '"' << """ at " << sip[i] / 1000 << ",0.6 front" << endl;
2001	}
2002	gnuout.close();
2003	}
2004
2005
2006	ofstream gnuoutbis;
2007	string gnuscriptbis;
2008	gnuscriptbis = outputpath + "/labelsIPend.p";
2009
2010	gnuoutbis.open(gnuscriptbis.c_str());
2011
2012	if (gnuoutbis.fail()) {
2013	cerr << "Warning: problem with the file " << gnuscriptbis << "!" << endl;
2014	} else {
2015
2016	for (int i(0); i < sip.size(); ++i) {
2017	gnuoutbis << "set label " << i + 1 << " " << '"' << '"' << endl;
2018	}
2019	gnuoutbis.close();
2020	}
2021	}
2022
2023	//Preparing the files for the plotting of the load on the collimators
2024
2025	vector <int> removedCollimators;
2026
2027	for (int i(0); i < lat.resColli.size(); ++i) {
2028	if (lat.resColli[i]->nsig > 100) {
2029	removedCollimators.push_back(1);
2030	} else {
2031	removedCollimators.push_back(0);
2032	}
2033	}
2034
2035	vector <int> newnhitcolli, ipsnew;
2036
2037	for (int i(0); i < nhitcollio.size(); ++i) {
2038	if (removedCollimators[i] == 0) {
2039	newnhitcolli.push_back(nhitcollio[i]);
2040	ipsnew.push_back(lat.ips[i]);
2041	}
2042	}
2043
2044	ofstream outplot3;
2045	string plot3;
2046	plot3 = outputpath + "/PlotRunSummary3.dat";
2047
2048	outplot3.open(plot3.c_str());
2049
2050	if (outplot3.fail()) {
2051	cerr << "Warning: problem with the file " << plot3 << "!" << endl;
2052	} else {
2053
2054	for (int i(0); i < newnhitcolli.size(); ++i) {
2055	outplot3 << setw(10) << i + 1 << setw(10) << newnhitcolli[i]*PlossPb / NLostTotal << setw(10) << ipsnew[i] << setw(20) << lat.reseau[ipsnew[i]]->NAME << endl;
2056	}
2057	outplot3.close();
2058	}
2059
2060	ofstream gnuout2;
2061	string gnuscript2;
2062	gnuscript2 = outputpath + "/labelscolli.p";
2063
2064	gnuout2.open(gnuscript2.c_str());
2065
2066	if (gnuout2.fail()) {
2067	cerr << "Warning: problem with the file " << gnuscript2 << "!" << endl;
2068	} else {
2069
2070	for (int i(0); i < newnhitcolli.size(); ++i) {
2071	gnuout2 << "set label " << i + 1 << " " << '"' << lat.reseau[ipsnew[i]]->NAME << '"' << """ at " << i << ",13 rotate front" << endl;
2072	}
2073	gnuout2.close();
2074	}
2075
2076	ofstream gnuout3;
2077	string gnuscript3;
2078	gnuscript3 = outputpath + "/labelscolliend.p";
2079
2080	gnuout3.open(gnuscript3.c_str());
2081
2082	if (gnuout3.fail()) {
2083	cerr << "Warning: problem with the file " << gnuscript3 << "!" << endl;
2084	} else {
2085
2086	for (int i(0); i < newnhitcolli.size(); ++i) {
2087	gnuout3 << "set label " << i + 1 << " " << '"' << '"' << endl;
2088	}
2089	gnuout3.close();
2090	}
2091
2092	//Preparing files for the plotting of the number of particles left in the beam, the number of particles that have hit a collimator and the number of particles that have hit the aperture
2093
2094	ofstream outplot4;
2095	string plot4;
2096	plot4 = outputpath + "/PlotRunSummary4.dat";
2097
2098	outplot4.open(plot4.c_str());
2099
2100	if (outplot4.fail()) {
2101	cerr << "Warning: problem with the file " << plot4 << "!" << endl;
2102	} else {
2103
2104	for (int i(0); i < asumrem.size(); ++i) {
2105	outplot4 << setw(10) << i << setw(10) << asumrem[i] << setw(10) << asumhitcolli[i] << setw(10) << asumhits[i] << endl;
2106	}
2107	outplot4.close();
2108	}
2109
2110	//Writing PlotRunSummary1.p
2111
2112	ofstream PlotRun1;
2113	string plotrunname1;
2114	plotrunname1 = outputpath + "/PlotRunSummary1.p";
2115
2116	PlotRun1.open(plotrunname1.c_str());
2117
2118	if (PlotRun1.fail()) {
2119	cerr << "Warning: problem with the file " << PlotRun1 << "!" << endl;
2120	} else {
2121
2122	PlotRun1 << "#Gnuplot script file for the PlotRunSummary (Plot of the loss map)" << endl << endl;
2123	PlotRun1 << "set autoscale xmax" << endl;
2124	PlotRun1 << "set autoscale y" << endl;
2125	PlotRun1 << "set key default" << endl;
2126	PlotRun1 << "set xrange [-1:]" << endl;
2127	PlotRun1 << "set ylabel \"P\'(W/m)\"" << endl;
2128	PlotRun1 << "set xlabel \"Distance from starting point (km)\"" << endl;
2129	PlotRun1 << "set title \"Loss map\"" << endl << endl;
2130	if (sip.size() != 0) {
2131	PlotRun1 << "load \'labelsIP.p\'" << endl << endl;
2132	}
2133	PlotRun1 << "plot \"PlotRunSummary1.dat\" using 2:1 title \'Losses\' with boxes";
2134	if (sip.size() != 0) {
2135	PlotRun1 << ",\"PlotRunSummary2.dat\" using 2:1 title \'IP\' with impulses" << endl << endl;
2136	PlotRun1 << "load \'labelsIPend.p\'" << endl;
2137	}
2138	PlotRun1.close();
2139	}
2140
2141
2142	//Writing PlotRunSummary2.p
2143
2144	ofstream PlotRun2;
2145	string plotrunname2;
2146	plotrunname2 = outputpath + "/PlotRunSummary2.p";
2147
2148	PlotRun2.open(plotrunname2.c_str());
2149
2150	if (PlotRun2.fail()) {
2151	cerr << "Warning: problem with the file " << PlotRun2 << "!" << endl;
2152	} else {
2153
2154	PlotRun2 << "#Gnuplot script file for the PlotRunSummary (Plot of the collimator load)" << endl << endl;
2155	PlotRun2 << "set autoscale" << endl;
2156	PlotRun2 << "set key default" << endl;
2157	PlotRun2 << "set style fill solid 0.5" << endl;
2158	PlotRun2 << "set ylabel \"P_coll (W)\"" << endl;
2159	PlotRun2 << "set xlabel \"Collimators\"" << endl;
2160	PlotRun2 << "set title \"Collimator load distribution\"" << endl << endl;
2161	PlotRun2 << "load \'labelscolli.p\'" << endl << endl;
2162	PlotRun2 << " plot \"PlotRunSummary3.dat\" using 1:2 title \'Collimator load\' with boxes" << endl << endl;
2163	PlotRun2 << "load \'labelscolliend.p\'" << endl;
2164	PlotRun2.close();
2165	}
2166
2167	//Writing PlotRunSummary3.p
2168
2169	ofstream PlotRun3;
2170	string plotrunname3;
2171	plotrunname3 = outputpath + "/PlotRunSummary3.p";
2172
2173	PlotRun3.open(plotrunname3.c_str());
2174
2175	if (PlotRun3.fail()) {
2176	cerr << "Warning: problem with the file " << PlotRun3 << "!" << endl;
2177	} else {
2178
2179	PlotRun3 << "#Gnuplot script file for the PlotRunSummary (Plot of the development after the first impact on collimator)" << endl << endl;
2180	PlotRun3 << "set autoscale" << endl;
2181	PlotRun3 << "set key default" << endl;
2182	PlotRun3 << "set ylabel \"Number of particles\"" << endl;
2183	PlotRun3 << "set xlabel \"N_turn\"" << endl;
2184	PlotRun3 << "set title \"Time Development after first impact on collimator\"" << endl << endl;
2185	PlotRun3 << "set logscale y" << endl << endl;
2186	PlotRun3 << " plot \"PlotRunSummary4.dat\" using 1:2 title \'particles in beam\' with line, \"PlotRunSummary4.dat\" using 1:3 title \'particles lost on collimators\' with line, \"PlotRunSummary4.dat\" using 1:4 title \'particles lost elsewhere\' with line" << endl << endl;
2187	PlotRun3 << "unset logscale y" << endl;
2188	PlotRun3.close();
2189	}
2190
2191	}
2192
2193	void Simulation::PlotTrajectory(const double& emix, const double& emiy, const int& ibeg, const int& iend, vector <vector <double> >& xco, vector <vector <double> >& yco, string outputpath)
2194	{
2195
2196	vector <double> skmn, varplotx, varploty;
2197
2198	for (int i(ibeg); i <= iend; ++i) {
2199	skmn.push_back(lat.reseau[i]->S / 1000);
2200	varplotx.push_back(1000 * 6 * sqrt(lat.reseau[i]->BETX * emix));
2201	varploty.push_back(1000 * 6 * sqrt(lat.reseau[i]->BETY * emiy));
2202	}
2203
2204	ofstream outplot5;
2205	string plot5;
2206	plot5 = outputpath + "/PlotTrajectory1.dat";
2207
2208	outplot5.open(plot5.c_str());
2209
2210	if (outplot5.fail()) {
2211	cerr << "Warning: problem with the file " << plot5 << "!" << endl;
2212	} else {
2213
2214	for (int i(0); i < skmn.size(); ++i) {
2215	outplot5 << setw(10) << skmn[i] << setw(10) << varplotx[i] << setw(10) << -varplotx[i] << setw(10) << varploty[i] << setw(10) << -varploty[i] << endl;
2216	}
2217	outplot5.close();
2218	}
2219
2220
2221	ofstream outplotxco;
2222	string plotxco;
2223	plotxco = outputpath + "/PlotTrajectory5.dat";
2224
2225	outplotxco.open(plotxco.c_str());
2226
2227	if (outplotxco.fail()) {
2228	cerr << "Warning: problem with the file " << plotxco << "!" << endl;
2229	} else {
2230
2231	for (int i(ibeg); i < iend; ++i) {
2232	outplotxco << setw(16) << lat.reseau[i]->S / 1000;
2233	for (int j(0); j < xco[i].size(); ++j) {
2234	if (xco[i][j] != 100) {
2235	outplotxco << setw(16) << 1000 * xco[i][j];
2236	} else {
2237	outplotxco << setw(16) << 0;
2238	}
2239	}
2240	outplotxco << endl;
2241	}
2242	outplotxco.close();
2243	}
2244
2245	ofstream gnuout4;
2246	string gnuscript4;
2247	gnuscript4 = outputpath + "/plotxco.p";
2248
2249	gnuout4.open(gnuscript4.c_str());
2250
2251	if (gnuout4.fail()) {
2252	cerr << "Warning: problem with the file " << gnuscript4 << "!" << endl;
2253	} else {
2254
2255	gnuout4 << "plot " << '"' << "PlotTrajectory1.dat" << '"' << " using 1:2 title " << '"' << "Beam envelope" << '"' << " with lines ls 1, ";
2256	gnuout4 << '"' << "PlotTrajectory1.dat" << '"' << " using 1:3 title " << '"' << '"' << " with lines ls 1,";
2257
2258
2259
2260	for (int i(0); i < xco[0].size() - 1; ++i) {
2261	gnuout4 << '"' << "PlotTrajectory5.dat" << '"' << " using 1:" << i + 2 << " title " << '"' << '"' << " with lines,";
2262	}
2263	gnuout4 << '"' << "PlotTrajectory5.dat" << '"' << " using 1:" << xco[0].size() + 1 << " title " << '"' << '"' << " with lines, ";
2264	gnuout4.close();
2265	}
2266
2267	ofstream outplotyco;
2268	string plotyco;
2269	plotyco = outputpath + "/PlotTrajectory6.dat";
2270
2271	outplotyco.open(plotyco.c_str());
2272
2273	if (outplotyco.fail()) {
2274	cerr << "Warning: problem with the file " << plotyco << "!" << endl;
2275	} else {
2276
2277	for (int i(ibeg); i < iend; ++i) {
2278	outplotyco << setw(16) << lat.reseau[i]->S / 1000;
2279	for (int j(0); j < yco[i].size(); ++j) {
2280	if (yco[i][j] != 100) {
2281	outplotyco << setw(16) << 1000 * yco[i][j];
2282	} else {
2283	outplotyco << setw(16) << 0;
2284	}
2285	}
2286	outplotyco << endl;
2287	}
2288	outplotyco.close();
2289	}
2290
2291	ofstream gnuout5;
2292	string gnuscript5;
2293	gnuscript5 = outputpath + "/plotyco.p";
2294
2295	gnuout5.open(gnuscript5.c_str());
2296
2297	if (gnuout5.fail()) {
2298	cerr << "Warning: problem with the file " << gnuscript5 << "!" << endl;
2299	} else {
2300
2301	gnuout5 << "plot " << '"' << "PlotTrajectory1.dat" << '"' << " using 1:4 title " << '"' << "Beam envelope" << '"' << " with lines ls 1, ";
2302	gnuout5 << '"' << "PlotTrajectory1.dat" << '"' << " using 1:5 title " << '"' << '"' << " with lines ls 1,";
2303
2304
2305
2306	for (int i(0); i < yco[0].size() - 1; ++i) {
2307	gnuout5 << '"' << "PlotTrajectory6.dat" << '"' << " using 1:" << i + 2 << " title " << '"' << '"' << " with lines,";
2308	}
2309	gnuout5 << '"' << "PlotTrajectory6.dat" << '"' << " using 1:" << yco[0].size() + 1 << " title " << '"' << '"' << " with lines, ";
2310	gnuout5.close();
2311	}
2312
2313	ofstream outplot6;
2314	string plot6;
2315	plot6 = outputpath + "/PlotTrajectory2.dat";
2316
2317	outplot6.open(plot6.c_str());
2318
2319	if (outplot6.fail()) {
2320	cerr << "Warning: problem with the file " << plot6 << "!" << endl;
2321	} else {
2322
2323	for (int i(0); i < lat.ips.size(); ++i) {
2324	outplot6 << setw(15) << 15 << setw(15) << -15 << setw(15) << lat.reseau[lat.ips[i]]->S / 1000 << endl;
2325	}
2326	outplot6.close();
2327	}
2328
2329	ofstream outplot7;
2330	string plot7;
2331	plot7 = outputpath + "/PlotTrajectory3.dat";
2332
2333	outplot7.open(plot7.c_str());
2334
2335	if (outplot7.fail()) {
2336	cerr << "Warning: problem with the file " << plot7 << "!" << endl;
2337	} else {
2338
2339	for (int i(0); i < lat.ip.size(); ++i) {
2340	outplot7 << setw(15) << 15 << setw(15) << -15 << setw(15) << lat.reseau[lat.ip[i]]->S / 1000 << endl;
2341	}
2342	outplot7.close();
2343	}
2344
2345	ofstream outplotis;
2346	string plotis;
2347	plotis = outputpath + "/PlotTrajectory8.dat";
2348
2349	outplotis.open(plotis.c_str());
2350
2351	if (outplotis.fail()) {
2352	cerr << "Warning: problem with the file " << plotis << "!" << endl;
2353	} else {
2354
2355	for (int i(0); i < lat.is.size(); ++i) {
2356	outplotis << setw(15) << 15 << setw(15) << -15 << setw(15) << lat.reseau[lat.is[i]]->S / 1000 << endl;
2357	}
2358	outplotis.close();
2359	}
2360
2361	vector <int> indice;
2362	vector <double> aperx1bis, apery1bis, shbis;
2363
2364	for (int k(0); k < lat.sh.size(); ++k) {
2365	if ((lat.sh[k] > lat.reseau[ibeg]->S) && (lat.sh[k] < lat.reseau[iend]->S)) {
2366	indice.push_back(1);
2367	} else {
2368	indice.push_back(0);
2369	}
2370	}
2371
2372	for (int k(0); k < lat.sh.size(); ++k) {
2373	if (indice[k] == 1) {
2374	aperx1bis.push_back(lat.aperx1[k]);
2375	apery1bis.push_back(lat.apery1[k]);
2376	shbis.push_back(lat.sh[k]);
2377	}
2378	}
2379
2380
2381	ofstream outplot8;
2382	string plot8;
2383	plot8 = outputpath + "/PlotTrajectory4.dat";
2384
2385	outplot8.open(plot8.c_str());
2386
2387	if (outplot8.fail()) {
2388	cerr << "Warning: problem with the file " << plot8 << "!" << endl;
2389	} else {
2390
2391	for (int i(0); i < shbis.size(); ++i) {
2392	outplot8 << setw(15) << aperx1bis[i] * 1000 << setw(15) << -aperx1bis[i] * 1000 << setw(15) << apery1bis[i] * 1000 << setw(15) << -apery1bis[i] * 1000 << setw(15) << shbis[i] / 1000 << endl;
2393	}
2394	outplot8.close();
2395	}
2396
2397	ofstream outplot9;
2398	string plot9;
2399	plot9 = outputpath + "/PlotTrajectory7.dat";
2400
2401	outplot9.open(plot9.c_str());
2402
2403	if (outplot9.fail()) {
2404	cerr << "Warning: problem with the file " << plot9 << "!" << endl;
2405	} else {
2406
2407	for (int i(ibeg); i < iend; ++i) {
2408	outplot9 << setw(15) << skmn[i] << setw(15) << lat.reseau[i]->XC * 1000 << setw(15) << lat.reseau[i]->YC * 1000 << endl;
2409	}
2410	outplot9.close();
2411	}
2412
2413
2414	gnuout4.open(gnuscript4.c_str(), ios::out \| ios::app);
2415
2416	if (gnuout4.fail()) {
2417	cerr << "Warning: problem with the file " << gnuscript4 << "!" << endl;
2418	} else {
2419	gnuout4 << '"' << "PlotTrajectory2.dat" << '"' << " using 3:1 title " << '"' << "ips" << '"' << " with impulses, ";
2420	gnuout4 << '"' << "PlotTrajectory3.dat" << '"' << " using 3:1 title " << '"' << "ip" << '"' << " with impulses lt 3 lw 6, ";
2421	gnuout4 << '"' << "PlotTrajectory8.dat" << '"' << " using 3:1 title " << '"' << "is" << '"' << " with impulses lt 4 lw 2, ";
2422	gnuout4 << '"' << "PlotTrajectory4.dat" << '"' << " using 5:1 title " << '"' << "Aperture" << '"' << " with lines ls 2, ";
2423	gnuout4 << '"' << "PlotTrajectory4.dat" << '"' << " using 5:2 title " << '"' << '"' << " with lines ls 2, ";
2424	gnuout4 << '"' << "PlotTrajectory7.dat" << '"' << " using 1:2 title " << '"' << "Reference particle" << '"' << " with lines lt -1 lw 2" << endl;
2425
2426	gnuout4.close();
2427	}
2428
2429	gnuout5.open(gnuscript5.c_str(), ios::out \| ios::app);
2430
2431	if (gnuout5.fail()) {
2432	cerr << "Warning: problem with the file " << gnuscript5 << "!" << endl;
2433	} else {
2434	gnuout5 << '"' << "PlotTrajectory2.dat" << '"' << " using 3:1 title " << '"' << "ips" << '"' << " with impulses, ";
2435	gnuout5 << '"' << "PlotTrajectory3.dat" << '"' << " using 3:1 title " << '"' << "ip" << '"' << " with impulses lt 3 lw 6, ";
2436	gnuout5 << '"' << "PlotTrajectory8.dat" << '"' << " using 3:1 title " << '"' << "is" << '"' << " with impulses lt 4 lw 2, ";
2437	gnuout5 << '"' << "PlotTrajectory4.dat" << '"' << " using 5:3 title " << '"' << "Aperture" << '"' << " with lines ls 2, ";
2438	gnuout5 << '"' << "PlotTrajectory4.dat" << '"' << " using 5:4 title " << '"' << '"' << " with lines ls 2, ";
2439	gnuout5 << '"' << "PlotTrajectory7.dat" << '"' << " using 1:3 title " << '"' << "Reference particle" << '"' << " with lines lt -1 lw 2" << endl;
2440
2441	gnuout5.close();
2442	}
2443
2444	//Writing PlotTrajectory1.p
2445
2446	ofstream PlotTraj1;
2447	string plottrajname1;
2448	plottrajname1 = outputpath + "/PlotTrajectory1.p";
2449
2450	PlotTraj1.open(plottrajname1.c_str());
2451
2452	if (PlotTraj1.fail()) {
2453	cerr << "Warning: problem with the file " << plottrajname1 << "!" << endl;
2454	} else {
2455
2456	PlotTraj1 << "#Gnuplot script file for the PlotTrajectory in x" << endl << endl;
2457	PlotTraj1 << "set autoscale x" << endl;
2458	PlotTraj1 << "set yrange [-25:25]" << endl;
2459	PlotTraj1 << " set style line 1 lw 1" << endl;
2460	PlotTraj1 << "set ylabel \"x [mm]\"" << endl;
2461	PlotTraj1 << "set xlabel \"S-S_IP1 [km]\"" << endl;
2462	PlotTraj1 << "set title \"Trajectory in x\"" << endl << endl;
2463	PlotTraj1 << "set style line 2 lw 1 lt 6" << endl << endl;
2464	PlotTraj1 << "set key outside below" << endl << endl;
2465	PlotTraj1 << "load \'plotxco.p\'" << endl;
2466	PlotTraj1 << "unset key" << endl;
2467	PlotTraj1.close();
2468	}
2469
2470	//Writing PlotTrajectory2.p
2471
2472	ofstream PlotTraj2;
2473	string plottrajname2;
2474	plottrajname2 = outputpath + "/PlotTrajectory2.p";
2475
2476	PlotTraj2.open(plottrajname2.c_str());
2477
2478	if (PlotTraj2.fail()) {
2479	cerr << "Warning: problem with the file " << plottrajname2 << "!" << endl;
2480	} else {
2481
2482	PlotTraj2 << "#Gnuplot script file for the PlotTrajectory in y" << endl << endl;
2483	PlotTraj2 << "set autoscale x" << endl;
2484	PlotTraj2 << "set yrange [-25:25]" << endl;
2485	PlotTraj2 << " set style line 1 lw 1" << endl;
2486	PlotTraj2 << "set ylabel \"y [mm]\"" << endl;
2487	PlotTraj2 << "set xlabel \"S-S_IP1 [km]\"" << endl;
2488	PlotTraj2 << "set title \"Trajectory in x\"" << endl << endl;
2489	PlotTraj2 << "set style line 2 lw 1 lt 6" << endl << endl;
2490	PlotTraj2 << "set key outside below" << endl << endl;
2491	PlotTraj2 << "load \'plotyco.p\'" << endl;
2492	PlotTraj2 << "unset key" << endl;
2493	PlotTraj2.close();
2494	}
2495
2496	}
2497
2498
2499	void Simulation::PlotLossSpectra(const double& beamflag, const double& nparti, const vector <int>& asumrem, const double& Apr, const double& Zpr, vector <double>& hitso, vector <double>& Aphito, vector <double>& Zphito, const double& PlossPb, const int& npart, const int& kbunch, const int& taubeam, const vector <int>& nhitcollio, string outputpath, string inputpath)
2500	{
2501
2502	vector <double> sip;//vector of s coordinates corresponding to element names starting with 'ip'
2503	vector <int> inip;
2504	vector <string> NAMES_IP;
2505
2506	for (int i(0); i < lat.reseau.size(); ++i) {
2507	if (lat.reseau[i]->NAME.substr(0, 2) == "IP") {
2508	NAMES_IP.push_back(lat.reseau[i]->NAME);
2509	sip.push_back(lat.reseau[i]->S);
2510	inip.push_back(i);
2511	}
2512	}
2513
2514	string interactionPoint, selflag, PowerOrPartFlag, decdPoP;
2515	double fragcutoff, d1, d2, s1, s2;
2516
2517	ifstream entree;
2518
2519	string plotname;
2520	plotname = inputpath + "infoPlotLossSpectra.csv";
2521
2522	entree.open(plotname.c_str());
2523
2524	if (entree.fail()) {
2525	cerr << "Warning: problem with the file 'infoPlotLossSpectra.csv'." << endl;
2526	} else {
2527
2528	string word;
2529	double nbre;
2530
2531	while (!entree.eof()) {
2532
2533	getline(entree, word, ',');
2534
2535	if (word == "interactionPoint") {
2536	getline(entree, interactionPoint);
2537	} else if (word == "d1") {
2538	getline(entree, word);
2539	d1 = atof(word.c_str());
2540	} else if (word == "d2") {
2541	getline(entree, word);
2542	d2 = atof(word.c_str());
2543	} else if (word == "selflag") {
2544	getline(entree, selflag);
2545	} else if (word == "fragcutoff") {
2546	getline(entree, word);
2547	fragcutoff = atof(word.c_str());
2548	} else if (word == "PowerOrPartFlag") {
2549	getline(entree, PowerOrPartFlag);
2550	} else if (word == "decdPoP") {
2551	getline(entree, decdPoP);
2552	} else if (word == "s1") {
2553	getline(entree, word);
2554	s1 = atof(word.c_str());
2555	} else if (word == "s2") {
2556	getline(entree, word);
2557	s2 = atof(word.c_str());
2558	} else {
2559	getline(entree, word);
2560	}
2561	}//end while(!entree.eof())
2562	entree.close();
2563	}
2564
2565	double offx;
2566	int pos(0);
2567
2568	for (int j(0); j < NAMES_IP.size(); ++j) {
2569	if (NAMES_IP[j] != interactionPoint) {
2570	++pos;
2571	} else {
2572	break;
2573	}
2574	}
2575
2576
2577	if (inip.size() != 0) {
2578	offx = sip[pos];//s coordinate of the interaction point
2579	s1 = offx + d1;//starting point along the accelerator [m]
2580	s2 = offx + d2;//ending point along the accelerator [m]
2581	} else {
2582	offx = 0;
2583	}
2584
2585
2586	vector <double> dPoPeff;
2587	vector <double> hitsnew, Aphitnew, Zphitnew;
2588
2589	if (decdPoP == "Bigger than reference particle") {
2590
2591	for (int i(0); i < Aphito.size(); ++i) {
2592	dPoPeff.push_back(Zpr / Apr * Aphito[i] / Zphito[i] - 1);
2593	if (dPoPeff[i] < 0) {
2594	hitsnew.push_back(hitso[i]);
2595	Aphitnew.push_back(Aphito[i]);
2596	Zphitnew.push_back(Zphito[i]);
2597	}
2598	}
2599	} else if (decdPoP == "Smaller than reference particle") {
2600	for (int i(0); i < Aphito.size(); ++i) {
2601	dPoPeff.push_back(Zpr / Apr * Aphito[i] / Zphito[i] - 1);
2602	if (dPoPeff[i] > 0) {
2603	hitsnew.push_back(hitso[i]);
2604	Aphitnew.push_back(Aphito[i]);
2605	Zphitnew.push_back(Zphito[i]);
2606	}
2607	}
2608	}
2609
2610	if (decdPoP != "Anything") {
2611	hitso.clear();
2612	Aphito.clear();
2613	Zphito.clear();
2614
2615	for (int i(0); i < hitsnew.size(); ++i) {
2616	hitso.push_back(hitsnew[i]);
2617	Aphito.push_back(Aphitnew[i]);
2618	Zphito.push_back(Zphitnew[i]);
2619	}
2620	}
2621
2622
2623	int NLostTotal(0);
2624
2625	for (int i(0); i < nhitcollio.size(); ++i) {
2626	NLostTotal += nhitcollio[i];
2627	}
2628
2629	NLostTotal += hitso.size();
2630
2631	double NLostPb;
2632
2633	NLostPb = npart * kbunch / taubeam;
2634
2635	vector <int> asel, zsel;
2636	vector <string> legarr;
2637
2638	if ((Apr == 29) && (Zpr == 63)) {
2639	asel.push_back(63);
2640	asel.push_back(62);
2641	asel.push_back(61);
2642	asel.push_back(60);
2643	asel.push_back(59);
2644	asel.push_back(58);
2645	asel.push_back(62);
2646	asel.push_back(61);
2647	asel.push_back(60);
2648	asel.push_back(59);
2649	asel.push_back(58);
2650	asel.push_back(57);
2651	asel.push_back(61);
2652	asel.push_back(60);
2653	asel.push_back(59);
2654	asel.push_back(58);
2655	asel.push_back(57);
2656	asel.push_back(56);
2657	asel.push_back(60);
2658	asel.push_back(59);
2659	asel.push_back(58);
2660	asel.push_back(57);
2661	asel.push_back(56);
2662	asel.push_back(55);
2663
2664	zsel.push_back(29);
2665	zsel.push_back(29);
2666	zsel.push_back(29);
2667	zsel.push_back(29);
2668	zsel.push_back(29);
2669	zsel.push_back(29);
2670	zsel.push_back(28);
2671	zsel.push_back(28);
2672	zsel.push_back(28);
2673	zsel.push_back(28);
2674	zsel.push_back(28);
2675	zsel.push_back(28);
2676	zsel.push_back(27);
2677	zsel.push_back(27);
2678	zsel.push_back(27);
2679	zsel.push_back(27);
2680	zsel.push_back(27);
2681	zsel.push_back(27);
2682	zsel.push_back(26);
2683	zsel.push_back(26);
2684	zsel.push_back(26);
2685	zsel.push_back(26);
2686	zsel.push_back(26);
2687	zsel.push_back(26);
2688
2689	legarr.push_back("Cu^{63}");
2690	legarr.push_back("Cu^{62}");
2691	legarr.push_back("Cu^{61}");
2692	legarr.push_back("Cu^{60}");
2693	legarr.push_back("Cu^{59}");
2694	legarr.push_back("Cu^{58}");
2695	legarr.push_back("Ni^{62}");
2696	legarr.push_back("Ni^{61}");
2697	legarr.push_back("Ni^{60}");
2698	legarr.push_back("Ni^{59}");
2699	legarr.push_back("Ni^{58}");
2700	legarr.push_back("Ni^{57}");
2701	legarr.push_back("Co^{61}");
2702	legarr.push_back("Co^{60}");
2703	legarr.push_back("Co^{59}");
2704	legarr.push_back("Co^{58}");
2705	legarr.push_back("Co^{57}");
2706	legarr.push_back("Co^{56}");
2707	legarr.push_back("Fe^{60}");
2708	legarr.push_back("Fe^{59}");
2709	legarr.push_back("Fe^{58}");
2710	legarr.push_back("Fe^{57}");
2711	legarr.push_back("Fe^{56}");
2712	legarr.push_back("Fe^{55}");
2713	} else if ((Apr == 197) && (Zpr == 79)) {
2714	asel.push_back(197);
2715	asel.push_back(196);
2716	asel.push_back(195);
2717	asel.push_back(194);
2718	asel.push_back(193);
2719	asel.push_back(192);
2720	asel.push_back(195);
2721	asel.push_back(194);
2722	asel.push_back(193);
2723	asel.push_back(192);
2724	asel.push_back(191);
2725	asel.push_back(190);
2726	asel.push_back(193);
2727	asel.push_back(192);
2728	asel.push_back(191);
2729	asel.push_back(190);
2730	asel.push_back(189);
2731	asel.push_back(188);
2732	asel.push_back(191);
2733	asel.push_back(190);
2734	asel.push_back(189);
2735	asel.push_back(188);
2736	asel.push_back(187);
2737	asel.push_back(186);
2738
2739	zsel.push_back(79);
2740	zsel.push_back(79);
2741	zsel.push_back(79);
2742	zsel.push_back(79);
2743	zsel.push_back(79);
2744	zsel.push_back(79);
2745	zsel.push_back(78);
2746	zsel.push_back(78);
2747	zsel.push_back(78);
2748	zsel.push_back(78);
2749	zsel.push_back(78);
2750	zsel.push_back(78);
2751	zsel.push_back(77);
2752	zsel.push_back(77);
2753	zsel.push_back(77);
2754	zsel.push_back(77);
2755	zsel.push_back(77);
2756	zsel.push_back(77);
2757	zsel.push_back(76);
2758	zsel.push_back(76);
2759	zsel.push_back(76);
2760	zsel.push_back(76);
2761	zsel.push_back(76);
2762	zsel.push_back(76);
2763
2764	legarr.push_back("Au^{197}");
2765	legarr.push_back("Au^{196}");
2766	legarr.push_back("Au^{195}");
2767	legarr.push_back("Au^{194}");
2768	legarr.push_back("Au^{193}");
2769	legarr.push_back("Au^{192}");
2770	legarr.push_back("Pt^{195}");
2771	legarr.push_back("Pt^{194}");
2772	legarr.push_back("Pt^{193}");
2773	legarr.push_back("Pt^{192}");
2774	legarr.push_back("Pt^{191}");
2775	legarr.push_back("Pt^{190}");
2776	legarr.push_back("Ir^{193}");
2777	legarr.push_back("Ir^{192}");
2778	legarr.push_back("Ir^{191}");
2779	legarr.push_back("Ir^{190}");
2780	legarr.push_back("Ir^{189}");
2781	legarr.push_back("Ir^{188}");
2782	legarr.push_back("Os^{191}");
2783	legarr.push_back("Os^{190}");
2784	legarr.push_back("Os^{189}");
2785	legarr.push_back("Os^{188}");
2786	legarr.push_back("Os^{187}");
2787	legarr.push_back("Os^{186}");
2788	} else if ((Apr == 208) && (Zpr == 82)) {
2789
2790	asel.push_back(208);
2791	asel.push_back(207);
2792	asel.push_back(206);
2793	asel.push_back(205);
2794	asel.push_back(204);
2795	asel.push_back(203);
2796	asel.push_back(202);
2797	asel.push_back(201);
2798	asel.push_back(207);
2799	asel.push_back(206);
2800	asel.push_back(205);
2801	asel.push_back(204);
2802	asel.push_back(203);
2803	asel.push_back(202);
2804	asel.push_back(201);
2805	asel.push_back(200);
2806	asel.push_back(199);
2807	asel.push_back(205);
2808	asel.push_back(204);
2809	asel.push_back(203);
2810	asel.push_back(202);
2811	asel.push_back(201);
2812	asel.push_back(200);
2813	asel.push_back(199);
2814	asel.push_back(198);
2815	asel.push_back(197);
2816	asel.push_back(202);
2817	asel.push_back(201);
2818	asel.push_back(200);
2819	asel.push_back(199);
2820	asel.push_back(198);
2821	asel.push_back(197);
2822	asel.push_back(196);
2823	asel.push_back(195);
2824
2825	zsel.push_back(82);
2826	zsel.push_back(82);
2827	zsel.push_back(82);
2828	zsel.push_back(82);
2829	zsel.push_back(82);
2830	zsel.push_back(82);
2831	zsel.push_back(82);
2832	zsel.push_back(82);
2833	zsel.push_back(81);
2834	zsel.push_back(81);
2835	zsel.push_back(81);
2836	zsel.push_back(81);
2837	zsel.push_back(81);
2838	zsel.push_back(81);
2839	zsel.push_back(81);
2840	zsel.push_back(81);
2841	zsel.push_back(81);
2842	zsel.push_back(80);
2843	zsel.push_back(80);
2844	zsel.push_back(80);
2845	zsel.push_back(80);
2846	zsel.push_back(80);
2847	zsel.push_back(80);
2848	zsel.push_back(80);
2849	zsel.push_back(80);
2850	zsel.push_back(80);
2851	zsel.push_back(79);
2852	zsel.push_back(79);
2853	zsel.push_back(79);
2854	zsel.push_back(79);
2855	zsel.push_back(79);
2856	zsel.push_back(79);
2857	zsel.push_back(79);
2858	zsel.push_back(79);
2859
2860	legarr.push_back("Pb^{208}");
2861	legarr.push_back("Pb^{207}");
2862	legarr.push_back("Pb^{206}");
2863	legarr.push_back("Pb^{205}");
2864	legarr.push_back("Pb^{204}");
2865	legarr.push_back("Pb^{203}");
2866	legarr.push_back("Pb^{202}");
2867	legarr.push_back("Pb^{201}");
2868	legarr.push_back("Tl^{207}");
2869	legarr.push_back("Tl^{206}");
2870	legarr.push_back("Tl^{205}");
2871	legarr.push_back("Tl^{204}");
2872	legarr.push_back("Tl^{203}");
2873	legarr.push_back("Tl^{202}");
2874	legarr.push_back("Tl^{201}");
2875	legarr.push_back("Tl^{200}");
2876	legarr.push_back("Tl^{199}");
2877	legarr.push_back("Hg^{205}");
2878	legarr.push_back("Hg^{204}");
2879	legarr.push_back("Hg^{203}");
2880	legarr.push_back("Hg^{202}");
2881	legarr.push_back("Hg^{201}");
2882	legarr.push_back("Hg^{200}");
2883	legarr.push_back("Hg^{199}");
2884	legarr.push_back("Hg^{198}");
2885	legarr.push_back("Hg^{197}");
2886	legarr.push_back("Au^{202}");
2887	legarr.push_back("Au^{201}");
2888	legarr.push_back("Au^{200}");
2889	legarr.push_back("Au^{199}");
2890	legarr.push_back("Au^{198}");
2891	legarr.push_back("Au^{197}");
2892	legarr.push_back("Au^{196}");
2893	legarr.push_back("Au^{195}");
2894	} else {
2895
2896	cerr << "Warning! It is not possible to mark which isotopes are lost when the reference particle has mass number " << Apr << " and charge " << Zpr << "!" << endl;
2897	selflag = "No";
2898	}
2899
2900	double aymax;
2901	vector <vector <double> > plosslocalt;
2902	vector <double> plosslocal, zlosslocal;
2903	vector <int> nasel, naselr;
2904
2905	if (selflag == "Yes") {
2906
2907	for (int i(0); i < asel.size(); ++i) {
2908	nasel.push_back(i);
2909	}
2910
2911	hitsnew.clear();
2912	Aphitnew.clear();
2913	Zphitnew.clear();
2914
2915	for (int i(0); i < hitso.size(); ++i) {
2916	if ((hitso[i] > s1) && (hitso[i] < s2)) {
2917	hitsnew.push_back(hitso[i]);
2918	Aphitnew.push_back(Aphito[i]);
2919	Zphitnew.push_back(Zphito[i]);
2920	}
2921	}
2922
2923	vector <int> alhitsp;
2924	vector <double> hitsp, hitsprest, Aphitrest, Zphitrest;
2925	vector <int> inh;
2926
2927	for (int k(0); k < nasel.size(); ++k) {
2928	int ia;
2929	ia = nasel[k];
2930	for (int i(0); i < hitso.size(); ++i) {
2931	if ((hitso[i] > s1) && (hitso[i] < s2) && (Aphito[i] == asel[ia]) && (Zphito[i] == zsel[ia])) {
2932	hitsp.push_back(hitso[i]);
2933	}
2934	}
2935	alhitsp.push_back(hitsp.size());
2936	hitsp.clear();
2937	for (int i(0); i < Zphitnew.size(); ++i) {
2938	if ((Zphitnew[i] == zsel[ia]) && (Aphitnew[i] == asel[ia])) {
2939	inh.push_back(1);
2940	} else {
2941	inh.push_back(0);
2942	}
2943	}
2944	for (int j(0); j < hitsnew.size(); ++j) {
2945	if (inh[j] == 0) {
2946	hitsprest.push_back(hitsnew[j]);
2947	Aphitrest.push_back(Aphitnew[j]);
2948	Zphitrest.push_back(Zphitnew[j]);
2949	}
2950	}
2951	inh.clear();
2952	hitsnew.clear();
2953	Aphitnew.clear();
2954	Zphitnew.clear();
2955
2956	for (int j(0); j < hitsprest.size(); ++j) {
2957	hitsnew.push_back(hitsprest[j]);
2958	Aphitnew.push_back(Aphitrest[j]);
2959	Zphitrest.push_back(Zphitrest[j]);
2960	}
2961	hitsprest.clear();
2962	Aphitrest.clear();
2963	Zphitrest.clear();
2964	}
2965	double maxalhitsp(alhitsp[0]);
2966
2967	for (int i(1); i < alhitsp.size(); ++i) {
2968	if (maxalhitsp < alhitsp[i]) {
2969	maxalhitsp = alhitsp[i];
2970	}
2971	}
2972
2973	for (int i(0); i < alhitsp.size(); ++i) {
2974	if (alhitsp[i] > maxalhitsp * fragcutoff) {
2975	naselr.push_back(nasel[i]);
2976	}
2977	}
2978
2979	hitsnew.clear();
2980	Aphitnew.clear();
2981	Zphitnew.clear();
2982
2983	for (int i(0); i < hitso.size(); ++i) {
2984	if ((hitso[i] > s1) && (hitso[i] < s2)) {
2985	hitsnew.push_back(hitso[i]);
2986	Aphitnew.push_back(Aphito[i]);
2987	Zphitnew.push_back(Zphito[i]);
2988	}
2989	}
2990
2991
2992	alhitsp.clear();
2993	hitsprest.clear();
2994	Aphitrest.clear();
2995	Zphitrest.clear();
2996	inh.clear();
2997
2998	for (int k(0); k < naselr.size(); ++k) {
2999	int ia;
3000	ia = naselr[k];
3001	hitsp.clear();
3002	for (int i(0); i < hitso.size(); ++i) {
3003	if ((hitso[i] > s1) && (hitso[i] < s2) && (Aphito[i] == asel[ia]) && (Zphito[i] == zsel[ia])) {
3004	hitsp.push_back(hitso[i]);
3005	}
3006	}
3007	alhitsp.push_back(hitsp.size());
3008
3009	vector <double> sep;
3010	zlosslocal.clear();
3011	for (double k(s1); k < s2; ++k) {
3012	sep.push_back(k);
3013	}
3014
3015	for (int i(0); i < sep.size() - 1; ++i) {
3016	zlosslocal.push_back(sep[i] + (sep[i + 1] - sep[i]) / 2);
3017	}
3018
3019	int count(0);
3020
3021	for (int i(0); i < sep.size() - 1; ++i) {
3022	for (int j(0); j < hitsp.size(); ++j) {
3023	if ((hitsp[j] >= sep[i]) && (hitsp[j] < sep[i + 1])) {
3024	++count;
3025	}
3026	}
3027	plosslocal.push_back(count);
3028	count = 0;
3029	}
3030
3031	plosslocalt.push_back(plosslocal);
3032	plosslocal.clear();
3033
3034	for (int i(0); i < Zphitnew.size(); ++i) {
3035	if ((Zphitnew[i] == zsel[ia]) && (Aphitnew[i] == asel[ia])) {
3036	inh.push_back(1);
3037	} else {
3038	inh.push_back(0);
3039	}
3040	}
3041	for (int j(0); j < hitsnew.size(); ++j) {
3042	if (inh[j] == 0) {
3043	hitsprest.push_back(hitsnew[j]);
3044	Aphitrest.push_back(Aphitnew[j]);
3045	Zphitrest.push_back(Zphitnew[j]);
3046	}
3047	}
3048
3049	inh.clear();
3050	hitsnew.clear();
3051	Aphitnew.clear();
3052	Zphitnew.clear();
3053
3054	for (int j(0); j < hitsprest.size(); ++j) {
3055	hitsnew.push_back(hitsprest[j]);
3056	Aphitnew.push_back(Aphitrest[j]);
3057	Zphitrest.push_back(Zphitrest[j]);
3058	}
3059	hitsprest.clear();
3060	Aphitrest.clear();
3061	Zphitrest.clear();
3062
3063	}//end for(int k(0); k<naslr.size(); ++k)
3064
3065	vector <double> sep1;
3066	zlosslocal.clear();
3067	for (double k(s1); k < s2; ++k) {
3068	sep1.push_back(k);
3069	}
3070
3071	for (int i(0); i < sep1.size() - 1; ++i) {
3072	zlosslocal.push_back(sep1[i] + (sep1[i + 1] - sep1[i]) / 2);
3073	}
3074
3075	int count(0);
3076
3077	for (int i(0); i < sep1.size() - 1; ++i) {
3078	for (int j(0); j < hitsprest.size(); ++j) {
3079	if ((hitsprest[j] >= sep1[i]) && (hitsprest[j] < sep1[i + 1])) {
3080	++count;
3081	}
3082	}
3083	plosslocal.push_back(count);
3084	count = 0;
3085	}
3086
3087	plosslocalt.push_back(plosslocal);
3088	plosslocal.clear();
3089
3090	double tempsum;
3091	vector <double> sumploss;
3092
3093	for (int l(0); l < plosslocalt.size(); ++l) {
3094	tempsum = plosslocalt[l][0];
3095	for (int h(1); h < plosslocalt[l].size(); ++h) {
3096	tempsum += plosslocalt[l][h];
3097	}
3098	sumploss.push_back(tempsum);
3099	tempsum = 0;
3100	}
3101
3102	if (PowerOrPartFlag == "Power") {
3103	double aymaxtemp;
3104	aymax = 1.5 * sumploss[0] * PlossPb / NLostTotal;
3105	for (int i(1); i < sumploss.size(); ++i) {
3106	aymaxtemp = 1.5 * sumploss[i] * PlossPb / NLostTotal;
3107	if (aymaxtemp > aymax) {
3108	aymax = aymaxtemp;
3109	}
3110	}
3111	} else if (PowerOrPartFlag == "Particles") {
3112	double aymaxtemp;
3113	aymax = (1.5 * sumploss[0] * PlossPb / NLostTotal) / 1000;
3114	for (int i(1); i < sumploss.size(); ++i) {
3115	aymaxtemp = (1.5 * sumploss[i] * PlossPb / NLostTotal) / 1000;
3116	if (aymaxtemp > aymax) {
3117	aymax = aymaxtemp;
3118	}
3119	}
3120	}
3121
3122	}//end if selflag == 'Yes'
3123	else if (selflag == "No") {
3124
3125	hitsnew.clear();
3126
3127	for (int i(0); i < hitso.size(); ++i) {
3128	if ((hitso[i] > s1) && (hitso[i] < s2)) {
3129	hitsnew.push_back(hitso[i]);
3130	}
3131	}
3132
3133	vector <double> sep2;
3134	zlosslocal.clear();
3135	for (double k(s1); k < s2; ++k) {
3136	sep2.push_back(k);
3137	}
3138
3139	for (int i(0); i < sep2.size() - 1; ++i) {
3140	zlosslocal.push_back(sep2[i] + (sep2[i + 1] - sep2[i]) / 2);
3141	}
3142
3143	int count(0);
3144
3145	for (int i(0); i < sep2.size() - 1; ++i) {
3146	for (int j(0); j < hitsnew.size(); ++j) {
3147	if ((hitsnew[j] >= sep2[i]) && (hitsnew[j] < sep2[i + 1])) {
3148	++count;
3149	}
3150	}
3151	plosslocal.push_back(count);
3152	count = 0;
3153	}
3154
3155	plosslocalt.push_back(plosslocal);
3156
3157	if (PowerOrPartFlag == "Power") {
3158	double aymaxtemp;
3159	aymax = 1.5 * plosslocal[0] * PlossPb / NLostTotal;
3160	for (int i(1); i < plosslocal.size(); ++i) {
3161	aymaxtemp = 1.5 * plosslocal[i] * PlossPb / NLostTotal;
3162	if (aymaxtemp > aymax) {
3163	aymax = aymaxtemp;
3164	}
3165	}
3166	} else if (PowerOrPartFlag == "Particles") {
3167	double aymaxtemp;
3168	aymax = (1.5 * plosslocal[0] * PlossPb / NLostTotal) / 1000;
3169	for (int i(1); i < plosslocal.size(); ++i) {
3170	aymaxtemp = (1.5 * plosslocal[i] * PlossPb / NLostTotal) / 1000;
3171	if (aymaxtemp > aymax) {
3172	aymax = aymaxtemp;
3173	}
3174	}
3175	}
3176	}//end if selflag == 'No'
3177
3178	if (aymax <= 0) {
3179	aymax = 50;
3180	}
3181
3182	double aymin(-0.4 * aymax);
3183
3184	vector <vector <double> > plosslocaltransp;
3185
3186	vector <double> temp;
3187	for (int j(0); j < plosslocalt.size(); ++j) {
3188	temp.push_back(0);
3189	}
3190
3191	for (int i(0); i < plosslocalt[0].size(); ++i) {
3192	plosslocaltransp.push_back(temp);
3193	}
3194
3195	temp.clear();
3196
3197	for (int i(0); i < plosslocaltransp.size(); ++i) {
3198	for (int j(0); j < plosslocaltransp[0].size(); ++j) {
3199	plosslocaltransp[i][j] = plosslocalt [j][i];
3200	}
3201	}
3202
3203
3204	ofstream gnuout5;
3205	string gnuscript5;
3206	gnuscript5 = outputpath + "/PlotLossSpectra.p";
3207
3208	gnuout5.open(gnuscript5.c_str());
3209
3210	if (gnuout5.fail()) {
3211	cerr << "Warning: problem with the file " << gnuscript5 << "!" << endl;
3212	} else {
3213
3214	gnuout5 << "#Gnuplot script file for the PlotLossSpectra(Plot of the loss map)" << endl << endl << endl;
3215	gnuout5 << "set autoscale y" << endl;
3216	gnuout5 << "set key default" << endl;
3217	gnuout5 << "set xrange [" << s1 << ':' << s2 << ']' << endl;
3218	gnuout5 << "set title " << '"' << "Loss map" << '"' << endl;
3219	if (PowerOrPartFlag == "Power") {
3220	gnuout5 << "set ylabel " << '"' << "Power load (W/m)(axis*20) / Optical function" << '"' << endl;
3221	} else {
3222	gnuout5 << "set ylabel " << '"' << "Particle loss rates (1/(ms)) / Optical function (axis /20)" << '"' << endl;
3223	}
3224	gnuout5 << "set xlabel " << '"' << "Distance along the accelerator" << '"' << endl;
3225	gnuout5 << "load " << "'labelsIPloss.p'" << endl;
3226	gnuout5 << "load " << "'labelsbendloss.p'" << endl;
3227	gnuout5 << "set style fill solid 0.5" << endl << endl << endl;
3228
3229	gnuout5.close();
3230	}
3231
3232	ofstream outplot10;
3233	string plot10;
3234	plot10 = outputpath + "/PlotLossSpectra1.dat";
3235
3236	outplot10.open(plot10.c_str());
3237
3238	if (outplot10.fail()) {
3239	cerr << "Warning: problem with the file " << plot10 << "!" << endl;
3240	} else {
3241
3242	for (int i(0); i < zlosslocal.size(); ++i) {
3243	outplot10 << setw(15) << zlosslocal[i] - offx;
3244	for (int j(0); j < plosslocaltransp[0].size(); ++j) {
3245	if (PowerOrPartFlag == "Power") {
3246	outplot10 << setw(15) << plosslocaltransp[i][j]*PlossPb / NLostTotal;
3247	} else {
3248	outplot10 << setw(15) << plosslocaltransp[i][j]*PlossPb / NLostTotal;
3249	}
3250	}
3251	outplot10 << endl;
3252	}
3253	outplot10.close();
3254	}
3255
3256
3257	ofstream outplot11;
3258	string plot11;
3259	plot11 = outputpath + "/PlotLossSpectra2.dat";
3260
3261	outplot11.open(plot11.c_str());
3262
3263	if (outplot11.fail()) {
3264	cerr << "Warning: problem with the file " << plot11 << "!" << endl;
3265	} else {
3266
3267	for (int i(0); i < lat.reseau.size(); ++i) {
3268	if ((lat.reseau[i]->S > s1) && (lat.reseau[i]->S < s2)) {
3269	outplot11 << setw(15) << lat.reseau[i]->S - offx << setw(15) << lat.reseau[i]->aperx * 50 << setw(15) << abs(lat.reseau[i]->DX) << setw(15) << lat.reseau[i]->BETX / 120 << endl;
3270	}
3271	}
3272	outplot11.close();
3273	}
3274
3275	ofstream outplot12;
3276	string plot12;
3277	plot12 = outputpath + "/PlotLossSpectra3.dat";
3278
3279	outplot12.open(plot12.c_str());
3280
3281	if (outplot12.fail()) {
3282	cerr << "Warning: problem with the file " << plot12 << "!" << endl;
3283	} else {
3284
3285	for (int i(0); i < sip.size(); ++i) {
3286	outplot12 << setw(10) << 0.5 << setw(10) << sip[i] - offx << endl;
3287	}
3288	outplot12.close();
3289	}
3290
3291	ofstream gnuout7;
3292	string gnuscript7;
3293	gnuscript7 = outputpath + "/labelsIPloss.p";
3294
3295	gnuout7.open(gnuscript7.c_str());
3296
3297	if (gnuout7.fail()) {
3298	cerr << "Warning: problem with the file " << gnuscript7 << "!" << endl;
3299	} else {
3300
3301	for (int i(0); i < sip.size(); ++i) {
3302	gnuout7 << "set label " << i + 1 << " " << '"' << NAMES_IP[i] << '"' << """ at " << sip[i] << ",0.6 front" << endl;
3303	}
3304	gnuout7.close();
3305	}
3306
3307	ofstream gnuout6;
3308	string gnuscript6;
3309	gnuscript6 = outputpath + "/PlotLossSpectra.p";
3310
3311	gnuout6.open(gnuscript6.c_str(), ios::out \| ios::app);
3312
3313	if (gnuout6.fail()) {
3314	cerr << "Warning: problem with the file " << gnuscript6 << "!" << endl;
3315	} else {
3316	if (selflag == "Yes") {
3317
3318	if (naselr. size() != 0) {
3319	gnuout6 << "plot " << '"' << "PlotLossSpectra1.dat" << '"' << " using 1:2 title " << '"' << legarr[naselr[0]] << '"' << " with boxes, ";
3320
3321	for (int i(1); i < naselr.size(); ++i) {
3322	gnuout6 << '"' << "PlotLossSpectra1.dat" << '"' << " using 1:" << i + 2 << " title " << '"' << legarr[naselr[i]] << '"' << " with boxes, ";
3323	}
3324	gnuout6 << '"' << "PlotLossSpectra1.dat" << '"' << " using 1:" << plosslocalt.size() + 1 << " title " << '"' << "others" << '"' << " with boxes, ";
3325	} else {
3326	gnuout6 << "plot" << '"' << "PlotLossSpectra1.dat" << '"' << " using 1:" << plosslocalt.size() + 1 << " title " << '"' << "others" << '"' << " with boxes, ";
3327	}
3328	} else {
3329	gnuout6 << "plot " << '"' << "PlotLossSpectra1.dat" << '"' << " using 1:2" << " title " << '"' << "others" << '"' << " with boxes, ";
3330	}
3331
3332	gnuout6 << '"' << "PlotLossSpectra2.dat" << '"' << " using 1:2" << " title " << '"' << "Aperture in x direction [m/50]" << '"' << " with lines, ";
3333	gnuout6 << '"' << "PlotLossSpectra2.dat" << '"' << " using 1:3" << " title " << '"' << "Dispersion function in x direction [m]" << '"' << " with lines, ";
3334	gnuout6 << '"' << "PlotLossSpectra2.dat" << '"' << " using 1:4" << " title " << '"' << "Beta function in x direction [120m]" << '"' << " with lines, ";
3335	if (sip.size() != 0) {
3336	gnuout6 << '"' << "PlotLossSpectra3.dat" << '"' << " using 2:1" << " title " << '"' << "IP" << '"' << " with impulses, ";
3337	}
3338
3339	gnuout6.close();
3340	}
3341
3342	vector <double> posx, precposx;
3343	vector <string> names;
3344
3345	if (s2 - s1 < 2000) {
3346
3347	vector <int> which;
3348
3349	for (int i(0); i < lat.reseau.size(); ++i) {
3350	if ((lat.reseau[i]->S > s1) && (lat.reseau[i]->S < s2)) {
3351	which.push_back(1);
3352	} else {
3353	which.push_back(0);
3354	}
3355	}
3356
3357
3358
3359	for (int i(0); i < which.size(); ++i) {
3360	if (which[i] == 1) {
3361	if ((lat.reseau[i]->KEYWORD == "RBEND") \|\| (lat.reseau[i]->KEYWORD == "SBEND") \|\| (lat.reseau[i]->KEYWORD == "QUADRUPOLE") \|\| (lat.reseau[i]->KEYWORD == "RCOLLIMATOR")) {
3362
3363	posx.push_back(lat.reseau[i]->S);
3364	precposx.push_back(lat.reseau[i - 1]->S);
3365	names.push_back(lat.reseau[i]->NAME);
3366
3367	}
3368	}
3369	}
3370
3371
3372
3373	ofstream outplot14;
3374	string plot14;
3375	plot14 = outputpath + "/PlotLossSpectra4.dat";
3376
3377	outplot14.open(plot14.c_str());
3378
3379	if (outplot14.fail()) {
3380	cerr << "Warning: problem with the file " << plot14 << "!" << endl;
3381	} else {
3382
3383	for (int i(0); i < posx.size(); ++i) {
3384	outplot14 << setw(15) << precposx[i] + (posx[i] - precposx[i]) / 2 << setw(15) << (posx[i] - precposx[i]) / 2 << setw(10) << -0.2 << setw(10) << 0.1 << endl;
3385	}
3386	outplot14.close();
3387	}
3388
3389	}
3390
3391	ofstream gnuout9;
3392	string gnuscript9;
3393	gnuscript9 = outputpath + "/labelsbendloss.p";
3394
3395	gnuout9.open(gnuscript9.c_str());
3396
3397	if (gnuout9.fail()) {
3398	cerr << "Warning: problem with the file " << gnuscript9 << "!" << endl;
3399	} else {
3400
3401	for (int i(0); i < names.size(); ++i) {
3402	gnuout9 << "set label " << i + 1 + sip.size() + 2 << " " << '"' << names[i] << '"' << """ at " << precposx[i] + (posx[i] - precposx[i]) / 2 << ",-0.3 rotate front " << endl;
3403	}
3404	gnuout9.close();
3405	}
3406
3407	ofstream gnuout8;
3408	string gnuscript8;
3409	gnuscript8 = outputpath + "/PlotLossSpectra.p";
3410
3411	gnuout8.open(gnuscript8.c_str(), ios::out \| ios::app);
3412
3413	if (gnuout8.fail()) {
3414	cerr << "Warning: problem with the file " << gnuscript8 << "!" << endl;
3415	} else {
3416	gnuout8 << '"' << "PlotLossSpectra4.dat" << '"' << " using 1:3:2:4" << " title " << '"' << '"' << " with boxxyerrorbars " << endl;
3417
3418
3419	gnuout8.close();
3420	}
3421
3422	}

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

Download in other formats: