Context Navigation

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

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

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

Download in other formats: