Context Navigation

simulation.cc @ 17

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

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

Context Navigation

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

Download in other formats: