#include <iostream>
#include <vector>
#include <string>
#include <cmath>
#include "lattice.h"
using namespace std;

int Lattice::count = 0;
int Lattice::turn = -2;
int Lattice::choicePart = 1;
int Lattice::eltOutNber = 1561;



//the following two constructors are just usefull for the initial tests

Lattice::Lattice(Element element1, Element element2, StandardCollimator stdcolli, FlukaCollimator flukacolli, MagneticCollimator magnetcolli, const vector <int>& ips,  int size, int npcle)
    : size_res(size), npart(npcle), ips(ips)
{
    addElement(&element1);
    addElement(&element2);
    addElement(&stdcolli);
    addElement(&flukacolli);
    addElement(&magnetcolli);
    addCollimator(new StandardCollimator(stdcolli));
    addCollimator(new FlukaCollimator(flukacolli));
    addCollimator(new MagneticCollimator(magnetcolli));

}


Lattice::Lattice(Element start, const vector <int>& ips)
    : size_res(1), npart(1), ips(ips)
{
    addElement(&start);
}



Lattice::~Lattice()
{

    for (int i(0); i < resColli.size(); ++i) {
        delete resColli[i];
    }
    ipcoll.push_back(0);
    ipcoll.clear();
}



void Lattice::setsize_res(int s)
{
    size_res = s;
}


void Lattice::addCollimator(Collimator* colli)
{

    resColli.push_back(colli);
    cocount.push_back(0);
}

void Lattice::addElement(Element* elt)
{

    reseau.push_back(elt);

}


double Lattice::time(Particle& p, const double& l, const double& betgam, const int& elt)
{

    double c(2.99792458e8); //speed of light [m/s]
    double beta(sqrt(betgam * betgam / (betgam * betgam + 1))); //relativistic beta
    double gamma(betgam / beta); //relativistic gamma
    double dist;
    double dist2;


    dist = sqrt(l * l + (p.coordonnees[1][0] - p.coordonnees[0][0]) * (p.coordonnees[1][0] - p.coordonnees[0][0]) + (p.coordonnees[1][2] - p.coordonnees[0][2]) * (p.coordonnees[1][2] - p.coordonnees[0][2]));

    if (elt != reseau.size() - 1) {

        dist2 = sqrt(l * l + (reseau[elt + 1]->XC - reseau[elt]->XC) * (reseau[elt + 1]->XC - reseau[elt]->XC) + (reseau[elt + 1]->YC - reseau[elt]->YC) * (reseau[elt + 1]->YC - reseau[elt]->YC));
    } else {

        dist2 = sqrt(l * l + (reseau[0]->XC - reseau[elt]->XC) * (reseau[0]->XC - reseau[elt]->XC) + (reseau[0]->YC - reseau[elt]->YC) * (reseau[0]->YC - reseau[elt]->YC));
    }

    p.dt = (dist - dist2) / (beta * c);

    return (dist / (beta * c));
}



void Lattice::outCoord(const Particle& p, const int& indic, const string& fileOut)
{

    ofstream output;

    int col(15);

    if ((turn == -2) && (indic == 1)) {
        output.open(fileOut.c_str());
    } else {
        output.open(fileOut.c_str(), ios::out | ios::app);
    }

    if (output.fail()) {

        cerr << "Warning: problem openning the file " << fileOut << "!" << endl;
    } else {

        output << setw(col) << p.coordonnees[1][0] << setw(col) << p.coordonnees[1][1] << setw(col) << p.coordonnees[1][2] << setw(col) << p.coordonnees[1][3] << setw(col) << p.coordonnees[1][4] << setw(col) << p.coordonnees[1][5] <<  endl;

    }

    output.close();

}


void Lattice::outPunct(const int& elt, const Particle& p, const Particle& p2, const double& var, string outputpath)
{

    ofstream outstream;
    string file;
    file = outputpath + "/coordinates_punctual.dat";

    int col(40);

    if (elt == eltOutNber) {

        if ((turn == -2) || (turn == -1)) {
            outstream.open(file.c_str());
        } else {
            outstream.open(file.c_str(), ios::out | ios::app);
        }

        if (outstream.fail()) {
            cerr << "Warning: problem openning the file " << file << "!" << endl;
        } else {

            outstream << setw(col) << p.coordonnees[1][0] << setw(col) << p.coordonnees[1][1] << setw(col) << p.coordonnees[1][2] << setw(col) << p.coordonnees[1][3] << setw(col) << p.coordonnees[1][4] << setw(col) << p.coordonnees[1][5] <<  endl;
        }

        outstream.close();
    }

}

void Lattice::outElt(const int& elt, const Particle& p, string outputpath, int& indication)
{

    ofstream tusors;
    string file;
    file = outputpath + "/coordinates_elt.dat";

    int col(35);

    if (elt == eltOutNber) {

        if (((turn == -2) || (turn == -1)) && (indication == 1)) {
            tusors.open(file.c_str());
            indication = 0;
        } else {
            tusors.open(file.c_str(), ios::out | ios::app);
        }

        if (tusors.fail()) {
            cerr << "Warning: problem openning the file " << file << "!" << endl;
        } else {

            tusors << p.coordonnees[1][0] << setw(col) << p.coordonnees[1][1] << setw(col) << p.coordonnees[1][2] << setw(col) << p.coordonnees[1][3] << setw(col) << p.coordonnees[1][4] << setw(col) << p.coordonnees[1][5] <<  endl;

        }

        tusors.close();
    }
}



void Lattice::outrf(const double& x1, const double& x2, string outputpath)
{

    ofstream print;
    string fich;
    fich = outputpath + "/valuestestrf.dat";

    int col(30);

    if ((turn == -1) || (turn == 0)) {
        print.open(fich.c_str());
    } else {
        print.open(fich.c_str(), ios::out | ios::app);
    }

    if (print.fail()) {
        cerr << "Warning: problem openning the file " << fich << "!" << endl;
    } else {

        print << setw(col) << setprecision(15) << x1 << setw(col) << setprecision(15) << x2 << endl;
    }

    print.close();
}

void Lattice::read(vector <Particle>& bunch)
{

    ifstream lecture;
    string nomfich;
    nomfich = "../sample/output_data_120GeV_500pt.txt";
    //nomfich = "test.txt";

    lecture.open(nomfich.c_str());

    if (lecture.fail()) {

        cerr << "Warning: problem with the file " << nomfich << " !" << endl;
    } else {

        bunch.clear();
        double var;
        int id;
        Particle p;
        string phrase;

        lecture >> ws;

        for (int k(0); k < 500; ++k) {

            //getline(lecture, phrase);

            lecture >> var;
            p.coordonnees[0][0] = var;
            p.coordonnees[1][0] = var;
            lecture >> var;
            p.coordonnees[0][1] = var;
            p.coordonnees[1][1] = var;
            lecture >> var;
            p.coordonnees[0][2] = var;
            p.coordonnees[1][2] = var;
            lecture >> var;
            p.coordonnees[0][3] = var;
            p.coordonnees[1][3] = var;
            lecture >> var;
            p.coordonnees[0][4] = var;
            p.coordonnees[1][4] = var;
            lecture >> var;
            p.coordonnees[0][5] = var;
            p.coordonnees[1][5] = var;
            //lecture >> id;
            p.Ap0 = 1;
            p.Zp0 = 1;

            bunch.push_back(p);
            bunch[bunch.size() - 1].inabs = 1;
        }

        lecture.close();
    }
}


/*
 * 
 */
void Lattice::trackensemblelinearnew(vector <Particle>& bunch, vector <Particle>& bunchhit, const int& nrev, const double& blowup2, const int& blowupperiod)
{

  /*nip: number of primary collimators in the accelerator; niph: number of primary collimators hit by the particles; nrevhitp: number of turns when the particle hit the primary collimators */
    int nip=0, niph=0, nrevhitp=0; 
    vector <int> iph;
    double blowup=0.0; /* sqrt(blowup2), beam blow up strength*/

    //we only take the primary collimators into account here, as well as the first and the last elements of the lattice
    nip = ip.size();//the number of primary collimators
    niph = nip + 1;

    //ip --> ip-1

    iph.push_back(0);

    for (int j(0); j < ip.size(); ++j) {
        iph.push_back(ip[j] - 1);
    }
    iph.push_back(reseau.size() - 1);
    blowup = sqrt(blowup2);

    for (int q(0); q < bunch.size(); ++q) {
        bunch[q].in = 1;
    }

  //  cout <<" "<<endl;
   //   cout <<"********************************"<<endl;
    cout << "Initialising LINEAR R-Matrix (solutions of Hill's equations using the Floquet theory, based on the Twiss parameters)." << endl;

    for (int i(0); i < niph; ++i) { //making matrices for the Twiss transform

        long double cx, sx, cy, sy;

        cx = cos(2 * M_PI * (reseau[iph[i + 1]]->MUX - reseau[iph[i]]->MUX));
        sx = sin(2 * M_PI * (reseau[iph[i + 1]]->MUX - reseau[iph[i]]->MUX));
        R11X.push_back(sqrt(reseau[iph[i + 1]]->BETX / reseau[iph[i]]->BETX) * (cx + reseau[iph[i]]->ALFX * sx));
        R12X.push_back(sqrt(reseau[iph[i + 1]]->BETX * reseau[iph[i]]->BETX)*sx);
        R21X.push_back(((reseau[iph[i]]->ALFX - reseau[iph[i + 1]]->ALFX)*cx - (1 + reseau[iph[i]]->ALFX * reseau[iph[i + 1]]->ALFX)*sx) / sqrt(reseau[iph[i + 1]]->BETX * reseau[iph[i]]->BETX));
        R22X.push_back(sqrt(reseau[iph[i]]->BETX / reseau[iph[i + 1]]->BETX) * (cx - reseau[iph[i + 1]]->ALFX * sx));


        cy = cos(2 * M_PI * (reseau[iph[i + 1]]->MUY - reseau[iph[i]]->MUY));
        sy = sin(2 * M_PI * (reseau[iph[i + 1]]->MUY - reseau[iph[i]]->MUY));
        R11Y.push_back(sqrt(reseau[iph[i + 1]]->BETY / reseau[iph[i]]->BETY) * (cy + reseau[iph[i]]->ALFY * sy));
        R12Y.push_back(sqrt(reseau[iph[i + 1]]->BETY * reseau[iph[i]]->BETY)*sy);
        R21Y.push_back(((reseau[iph[i]]->ALFY - reseau[iph[i + 1]]->ALFY)*cy - (1 + reseau[iph[i]]->ALFY * reseau[iph[i + 1]]->ALFY)*sy) / sqrt(reseau[iph[i + 1]]->BETY * reseau[iph[i]]->BETY));
        R22Y.push_back(sqrt(reseau[iph[i]]->BETY / reseau[iph[i + 1]]->BETY) * (cy - reseau[iph[i + 1]]->ALFY * sy));
    }

    int count(0);
    int total(bunch.size());
    for (int k(0); k < nrev; ++k) { //loop over the number of revolution
        vector <Particle> bunchtemp;

        ++turn;

        for (int w(0); w < bunch.size(); ++w) {
            bunchtemp.push_back(bunch[w]);
        }
        for (int i(0); i < niph; ++i) { //loop through the primary collimators
            for (int p(0); p < bunch.size(); ++p) { //loop through the particles
                if (bunch[p].in == 1) { //to assure the particle is still remainding in the accelerator

                    long double pdepth=0.0, pdepth2=0.0;

                    bunch[p].coordonnees[1][0] = R11X[i] * bunch[p].coordonnees[0][0] + R12X[i] * bunch[p].coordonnees[0][1] + (reseau[iph[i + 1]]->DX - R11X[i] * reseau[iph[i]]->DX - R12X[i] * reseau[iph[i]]->DPX) * bunch[p].coordonnees[0][4];
                    bunch[p].coordonnees[1][1] = R21X[i] * bunch[p].coordonnees[0][0] + R22X[i] * bunch[p].coordonnees[0][1] + (reseau[iph[i + 1]]->DPX - R21X[i] * reseau[iph[i]]->DX - R22X[i] * reseau[iph[i]]->DPX) * bunch[p].coordonnees[0][4];
                    bunch[p].coordonnees[1][2] = R11Y[i] * bunch[p].coordonnees[0][2] + R12Y[i] * bunch[p].coordonnees[0][3] + (reseau[iph[i + 1]]->DY - R11Y[i] * reseau[iph[i]]->DY - R12Y[i] * reseau[iph[i]]->DPY) * bunch[p].coordonnees[0][4];
                    bunch[p].coordonnees[1][3] = R21Y[i] * bunch[p].coordonnees[0][2] + R22Y[i] * bunch[p].coordonnees[0][3] + (reseau[iph[i + 1]]->DPY - R21Y[i] * reseau[iph[i]]->DY - R22Y[i] * reseau[iph[i]]->DPY) * bunch[p].coordonnees[0][4];

                    if (i < niph - 1) {
                        long double lcoll=0.0, sa=0.0, ca=0.0;

                        sa = sin(resColli[ipcoll[i]]->tcang);//sinus of the collimator's angle
                        lcoll = resColli[ipcoll[i]]->L;//length of the collimator
                        

                        if (sa == 0) {
                            pdepth = abs(bunch[p].coordonnees[1][0]) - resColli[ipcoll[i]]->hgap;//impact parameter at the beginning of the collimaator
                            pdepth2 = abs(bunch[p].coordonnees[1][0] + lcoll * bunch[p].coordonnees[1][1]) - resColli[ipcoll[i]]->hgap2; //impact parameter at the end of the collimator

                        } else {

                            long double xl, xsl;
                            ca = cos(resColli[ipcoll[i]]->tcang);
                            xl = bunch[p].coordonnees[1][0] * ca + bunch[p].coordonnees[1][2] * sa;
                            xsl = bunch[p].coordonnees[1][1] * ca + bunch[p].coordonnees[1][3] * sa;
                            pdepth = abs(xl) - resColli[ipcoll[i]]->hgap;
                            pdepth2 = abs(xl + lcoll * xsl) - resColli[ipcoll[i]]->hgap2;
                        }

                        if ((pdepth <= 0) && (pdepth2 <= 0)) { //we only continue with the particles that have not disappeared
                            bunch[p].coordonnees[0][0] = bunch[p].coordonnees[1][0];
                            bunch[p].coordonnees[0][1] = bunch[p].coordonnees[1][1];
                            bunch[p].coordonnees[0][2] = bunch[p].coordonnees[1][2];
                            bunch[p].coordonnees[0][3] = bunch[p].coordonnees[1][3];
                        } else {
                            bunch[p].in = false;
                            count = count + 1;
                            bunchhit.push_back(bunchtemp[p]);
                            bunch[p].nrevhitp = k + 1;
                            apdepth.push_back(pdepth);
                            apdepth2.push_back(pdepth2);
                            cocount[i] = cocount[i] + 1;
                        }
                    }//end (if i<niph)
                    else {//we only continue with the particles that have not disappeared
                        bunch[p].coordonnees[0][0] = bunch[p].coordonnees[1][0];
                        bunch[p].coordonnees[0][1] = bunch[p].coordonnees[1][1];
                        bunch[p].coordonnees[0][2] = bunch[p].coordonnees[1][2];
                        bunch[p].coordonnees[0][3] = bunch[p].coordonnees[1][3];
                    }
                }//end if in

                //the following part can be uncommented to have an output of the x- and  y-coordinates of a particle in the file coordinates.dat

                /*if(bunch[p].getidentification() == choicePart){
                  outCoord(bunch[p], i+1, "coordinates.dat");
                  }*/

            }//end loop over the particles
        }//end loop over i

       // cout << "Revolution: " << k + 1 << ", number of particles left: " << total - count << endl;

        if (total - count == 0) { //we return if all the particles are gone
            return;
        }

        vector <Particle> tempinabs;

        //we continue just with particles that are not lost
        for (int w(0); w < bunch.size(); ++w) {
            if (bunch[w].in != 0) {
                tempinabs.push_back(bunch[w]);
            }
        }

        bunch.clear();

        for (int w(0); w < tempinabs.size(); ++w) {
            bunch.push_back(tempinabs[w]);
        }

        tempinabs.clear();

        for (int p(0); p < bunch.size(); ++p) {
            if ((k + 1) % blowupperiod == 0) { //blowup/diffusion
                if (bunch[p].in == 1) {
                //    cout << "Blow-up!!" << k+1 << endl;
                    int im=0;
                    im = reseau.size() - 1;
                    bunch[p].coordonnees[0][0] = (bunch[p].coordonnees[0][0] - reseau[im]->DX * bunch[p].coordonnees[0][4]) * blowup + reseau[im]->DX * bunch[p].coordonnees[0][4];
                    bunch[p].coordonnees[0][1] = (bunch[p].coordonnees[0][1] - reseau[im]->DPX * bunch[p].coordonnees[0][4]) * blowup + reseau[im]->DPX * bunch[p].coordonnees[0][4];
                    bunch[p].coordonnees[0][2] = (bunch[p].coordonnees[0][2] - reseau[im]->DY * bunch[p].coordonnees[0][4]) * blowup + reseau[im]->DY * bunch[p].coordonnees[0][4];
                    bunch[p].coordonnees[0][3] = (bunch[p].coordonnees[0][3] - reseau[im]->DPY * bunch[p].coordonnees[0][4]) * blowup + reseau[im]->DPY * bunch[p].coordonnees[0][4];
                }
            }
        }

    }//end loop over k

}


void Lattice::trackensemblechrom(vector <Particle>& bunch, const int& irev, const int& i0, const int& im, const double& Apr, const double& Zpr, const double& wecolli, const double& betgam, const int& nonlinflag, const int& scaleorbit, double& attr1, const int& idpart, const int& idelt, const int& outcoord, const string& plotflag, vector <vector <double> >& xco, vector <vector <double> >& yco, string outputpath, int RFflag, int& indication)
{

    if (idpart >= 0) {
        if (idpart > bunch.size()) {
            cerr << "Warning, the particle that you want to spy using IDPART does not exist!" << "  idpart = " << idpart << ",  bunch.size = " << bunch.size() <<endl;
        } else {
            choicePart = idpart;
        }
    }

    if (idelt >= 0) {
        if (idelt >= reseau.size()) {
            cerr << "Warning, the element after which you want to spy using IDELT does not exist!" << endl;
        } else {
            eltOutNber = idelt;
        }
    }

    ++turn;

    if (nhitcolli.size() < resColli.size()) {
        for (int k(0); k < resColli.size(); ++k) {
            nhitcolli.push_back(0);
        }
    }

    if (plotflag == "Yes") {
        xco.clear();
        yco.clear();
        vector <double> temp1, temp2;
        for (int k(0); k < bunch.size(); ++k) {
            temp1.push_back(bunch[k].coordonnees[0][0]);
            temp2.push_back(bunch[k].coordonnees[0][2]);
        }

        xco.push_back(temp1);
        yco.push_back(temp2);

        temp1.clear();
        temp2.clear();
    }

    //we set the revolution number (trackensemblechrom is called two times during the first turn)
    int rev;

    if (irev == 0) {
        rev = 1;
    } else {
        rev = irev;
    }

    if (irev == 0) {
        cout << "Initialising NONLINEAR R-matrix." << endl;

        double K=0.0;
        double cx=0.0, sx=0.0, cy=0.0, sy=0.0;

        R11X.clear();
        R12X.clear();
        R21X.clear();
        R22X.clear();
        R11Y.clear();
        R12Y.clear();
        R21Y.clear();
        R22Y.clear();
        turn = -2;

        Lelem.push_back(0);
        R11X.push_back(0);
        R12X.push_back(0);
        R21X.push_back(0);
        R22X.push_back(0);
        R11Y.push_back(0);
        R12Y.push_back(0);
        R21Y.push_back(0);
        R22Y.push_back(0);
        R11XC.push_back(0);
        R12XC.push_back(0);
        R21XC.push_back(0);
        R22XC.push_back(0);
        R11YC.push_back(0);
        R12YC.push_back(0);
        R21YC.push_back(0);
        R22YC.push_back(0);

        for (int i(1); i < size_res; ++i) {

            Lelem.push_back(reseau[i]->S - reseau[i - 1]->S);
	    
            if (reseau[i]->K1L > 0) {
                if (Lelem[i] == 0) {
                    Lelem[i] = 0.001;
                }

                K = reseau[i]->K1L / Lelem[i];
                cx = cos(sqrt(K) * Lelem[i]);
                sx = sin(sqrt(K) * Lelem[i]);
                R11XC.push_back(0.5 * sqrt(K)*Lelem[i]*sx);
                R12XC.push_back(-0.5 * Lelem[i]*cx + 0.5 / sqrt(K)*sx);
                R21XC.push_back(0.5 * K * Lelem[i]*cx + 0.5 * sqrt(K)*sx);
                R22XC.push_back(0.5 * sqrt(K)*Lelem[i]*sx);
                cy = cosh(sqrt(K) * Lelem[i]);
                sy = sinh(sqrt(K) * Lelem[i]);
                R11YC.push_back(-0.5 * sqrt(K)*Lelem[i]*sy);
                R12YC.push_back(-0.5 * Lelem[i]*cy + 0.5 / sqrt(K)*sy);
                R21YC.push_back(-0.5 * K * Lelem[i]*cy - 0.5 * sqrt(K)*sy);
                R22YC.push_back(-0.5 * sqrt(K)*Lelem[i]*sy);
            } else if (reseau[i]->K1L < 0) {
                if (Lelem[i] == 0) {
                    Lelem[i] = 0.001;
                }

                K = -reseau[i]->K1L / Lelem[i];
                cx = cosh(sqrt(K) * Lelem[i]);
                sx = sinh(sqrt(K) * Lelem[i]);
                R11XC.push_back(-0.5 * sqrt(K)*Lelem[i]*sx);
                R12XC.push_back(-0.5 * Lelem[i]*cx + 0.5 / sqrt(K)*sx);
                R21XC.push_back(-0.5 * K * Lelem[i]*cx - 0.5 * sqrt(K)*sx);
                R22XC.push_back(-0.5 * sqrt(K)*Lelem[i]*sx);
                cy = cos(sqrt(K) * Lelem[i]);
                sy = sin(sqrt(K) * Lelem[i]);
                R11YC.push_back(0.5 * sqrt(K)*Lelem[i]*sy);
                R12YC.push_back(-0.5 * Lelem[i]*cy + 0.5 / sqrt(K)*sy);
                R21YC.push_back(0.5 * K * Lelem[i]*cy + 0.5 * sqrt(K)*sy);
                R22YC.push_back(0.5 * sqrt(K)*Lelem[i]*sy);
            } else {
                R11XC.push_back(0);
                R12XC.push_back(0);
                R21XC.push_back(0);
                R22XC.push_back(0);
                R11YC.push_back(0);
                R12YC.push_back(0);
                R21YC.push_back(0);
                R22YC.push_back(0);
            }

            cx = cos(2 * M_PI * (reseau[i]->MUX - reseau[i - 1]->MUX));
            sx = sin(2 * M_PI * (reseau[i]->MUX - reseau[i - 1]->MUX));
            R11X.push_back(sqrt(reseau[i]->BETX / reseau[i - 1]->BETX) * (cx + reseau[i - 1]->ALFX * sx));
            R12X.push_back(sqrt(reseau[i]->BETX * reseau[i - 1]->BETX)*sx);
            R21X.push_back(((reseau[i - 1]->ALFX - reseau[i]->ALFX)*cx - (1 + reseau[i - 1]->ALFX * reseau[i]->ALFX)*sx) / sqrt(reseau[i]->BETX * reseau[i - 1]->BETX));
            R22X.push_back(sqrt(reseau[i - 1]->BETX / reseau[i]->BETX) * (cx - reseau[i]->ALFX * sx));


            cy = cos(2 * M_PI * (reseau[i]->MUY - reseau[i - 1]->MUY));
            sy = sin(2 * M_PI * (reseau[i]->MUY - reseau[i - 1]->MUY));
            R11Y.push_back(sqrt(reseau[i]->BETY / reseau[i - 1]->BETY) * (cy + reseau[i - 1]->ALFY * sy));
            R12Y.push_back(sqrt(reseau[i]->BETY * reseau[i - 1]->BETY)*sy);
            R21Y.push_back(((reseau[i - 1]->ALFY - reseau[i]->ALFY)*cy - (1 + reseau[i - 1]->ALFY * reseau[i]->ALFY)*sy) / sqrt(reseau[i]->BETY * reseau[i - 1]->BETY));
            R22Y.push_back(sqrt(reseau[i - 1]->BETY / reseau[i]->BETY) * (cy - reseau[i]->ALFY * sy));

        }

    }


    // loop througt the elements in the machine

    for (int i(i0 + 1); i <= im; ++i) { //i is the element number along the ring (similar numbering as s)

        /*if(i == 2822){
          read(bunch);
          }*/

        for (int p(0); p < bunch.size(); ++p) { //loop througt the particles

            double dpopeff;

            if (bunch[p].inabs == 1) {

                //cout <<"Element " << i << endl;

                int icop(-1);

                for (int k(0); k < ips.size(); ++k) {
                    if (i == ips[k]) { //element is a collimator
                        icop = k;//The element is the 'icop'th collimator
                    }
                }

                if (icop != -1) { //the element we consider is a collimator

//cout << "icop = "<<icop << endl;

                    if (resColli[icop]->method == "standard") {

                        resColli[icop]->collipass(bunch[p], dpopeff, scaleorbit, R11X[i], R12X[i], R21X[i], R22X[i], R11Y[i], R12Y[i], R21Y[i], R22Y[i], reseau[i - 1]->DX, reseau[i - 1]->DPX, reseau[i - 1]->DY, reseau[i - 1]->DPY, reseau[i]->S - reseau[i - 1]->S, Apr, Zpr, betgam);

                        if (bunch[p].Ap0 != 0) {
                            bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5] + time(bunch[p], reseau[i]->L, betgam, i);
                        }

                    }
#if defined(FLUKA)
                    else if (resColli[icop]->method == "fluka") {

                        vector <Particle> temp;

                        if (p == 0) { //we send all the bunch at the same time to Fluka only one time per element

                            vector <Particle> bunchend;//bunch after the passage through Fluka

                            for (int y(0); y < bunch.size(); ++y) {
                                temp.push_back(bunch[y]);
                            }

                            resColli[icop]->collipassfluka(bunch, bunchend, conn, turn, momentum);


                            nhitcolli[icop] = bunch.size() - bunchend.size(); //nhitcolli(icop) gives the number particles getting lost in collimator number icop,

                            int creation(0);

                            bunch.clear();

                            for (int g(0); g < bunchend.size(); ++g) {
                                if (bunchend[g].getidentification() == 0) {
                                    bunch.push_back(bunchend[g]);
                                    break;
                                }
                            }

                            for (int g(0); g < bunchend.size(); ++g) {
                                if (bunchend[g].getidentification() == bunch[bunch.size() - 1].getidentification() + 1) {
                                    bunch.push_back(bunchend[g]);
                                } else if (bunchend[g].getidentification() == bunch[bunch.size() - 1].getidentification() + 10000) {
                                    bunch.push_back(bunchend[g]);
                                    ++creation;
                                } else if ((bunchend[g].getidentification() == bunch[bunch.size() - 1].getidentification()) && (bunchend[g].getidentification() != 0)) {
                                    bunch.push_back(bunchend[g]);
                                    ++creation;
                                } else if (bunchend[g].getidentification() == bunch[bunch.size() - 1].getidentification() - 10000 + 1) {
                                    bunch.push_back(bunchend[g]);
                                }
                            }

                            nhitcolli[icop] = nhitcolli[icop] + creation;

                            cout << creation << " particles created." << endl;

                            int uu(1);
                            int rr0, rr1;
                            long double tps;
                            bunch[0].coordonnees[0][5] = temp[0].coordonnees[0][5];
                            rr0 = bunch[0].getidentification();
                            bunch[0].setidentification(temp[0].getidentification());
                            for (int g(1); g < bunch.size(); ++g) {
                                rr1 = bunch[g].getidentification();
                                if (bunch[g].getidentification() == rr0 + 1) {
                                    bunch[g].setidentification(temp[uu].getidentification());
                                    bunch[g].coordonnees[0][5] = temp[uu].coordonnees[0][5];
                                    rr0 = rr1;
                                } else {
                                    tps = temp[uu - 1].coordonnees[0][5];
                                    while (bunch[g].getidentification() == rr1) {
                                        bunch[g].coordonnees[0][5] = tps;
                                        ++g;
                                    }
                                    --g;
                                    --uu;
                                }
                                ++uu;
                            }

                            for (int g(0); g < bunch.size(); ++g) {
                                bunch[g].coordonnees[1][0] = bunch[g].coordonnees[0][0];
                                bunch[g].coordonnees[1][1] = bunch[g].coordonnees[0][1];
                                bunch[g].coordonnees[1][2] = bunch[g].coordonnees[0][2];
                                bunch[g].coordonnees[1][3] = bunch[g].coordonnees[0][3];
                                bunch[g].coordonnees[1][4] = bunch[g].coordonnees[0][4];
                            }

                        }

                        bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5] + time(bunch[p], reseau[i]->L, betgam, i);

                        for (int k(0); k < bunch.size(); ++k) {
                            bunch[k].inabs = 1;
                        }

                        temp.clear();

                    }
#endif
                    else if (resColli[icop]->method == "magnetic") {

                        if (p == 0) {
                            resColli[icop]->hgap = resColli[icop]->hgap + resColli[icop]->deltaGap;
                            resColli[icop]->hgap2 = resColli[icop]->hgap2 + resColli[icop]->deltaGap;
                        }


                        resColli[icop]->collipass(bunch[p], dpopeff, scaleorbit, R11X[i], R12X[i], R21X[i], R22X[i], R11Y[i], R12Y[i], R21Y[i], R22Y[i], reseau[i - 1]->DX, reseau[i - 1]->DPX, reseau[i - 1]->DY, reseau[i - 1]->DPY, reseau[i]->S - reseau[i - 1]->S, Apr, Zpr, betgam);


                        if (bunch[p].Ap0 != 0) {
                            bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5] + time(bunch[p], reseau[i]->L, betgam, i);
                        }

                    } else if (resColli[icop]->method == "crystal") {

                        if (p == 0) { //we send all the bunch at the same time through collipassCrystal

                            int lost(bunch.size());

                            resColli[icop]->collipassCrystal(bunch, betgam, turn, outputpath);

                            nhitcolli[icop] = lost - bunch.size(); //nhitcolli(icop) gives the number particles getting lost in collimator number icop,

                        }

                        bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5] + time(bunch[p], reseau[i]->L, betgam, i);

                    } else {

                        cerr << "Error: unknown type of collimator, the method is not good defined!!" << endl;
                    }

                    //we test if the particle is lost in the preceeding collimator
                    if (bunch[p].Ap0 == 0) {
                        nhitcolli[icop] = nhitcolli[icop] + 1;
                        bunch[p].inabs = 0;
                    }

                } else {//element not collimator

                    if ((reseau[i]->KEYWORD == "RFCAVITY") && (RFflag == 1)) { //attention: voir jusqu ou aller avec le else (est-ce qu on controle l aperture hit?? pour le moment oui...)

                        //Note thate the phase is taken here to be equal to pi.

                        double period(rfharmonic / freqrf);
                        double omega;
                        double phase(0);
                        double phi;
                        double beta(sqrt(betgam * betgam / (betgam * betgam + 1))); //relativistic beta
                        long double c(2.99792458e8);//speed of light [m/s]
                        long double e(1.60218e-19);//elementary charge [C]


                        omega = 2 * M_PI * (freqrf / rfharmonic);
                        phi = phase + omega * bunch[p].coordonnees[0][5];

                        bunch[p].coordonnees[1][0] = bunch[p].coordonnees[0][0];
                        bunch[p].coordonnees[1][1] = bunch[p].coordonnees[0][1];
                        bunch[p].coordonnees[1][2] = bunch[p].coordonnees[0][2];
                        bunch[p].coordonnees[1][3] = bunch[p].coordonnees[0][3];

                        double attr;
                        attr = bunch[p].Zp0 * sin(phi) * (rfvoltage) / (beta * momentum);
                        attr1 =  attr;
                        bunch[p].coordonnees[1][4] = bunch[p].dpoporiginal + attr;//attention vraiment pas sur des parametre (surtout p dans la formule)

                        bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5];

                        //uncomment the following line to have output related to the rf-cavity (cf lattice.h)
                        //outrf(bunch[p].coordonnees[1][5], phi);


                    }


                    dpopeff = (bunch[p].Ap0 * Zpr) / (bunch[p].Zp0 * Apr) * (1 + bunch[p].coordonnees[0][4]) - 1;

                    if (Lelem[i] == 0) {

                        bunch[p].coordonnees[1][0] = bunch[p].coordonnees[0][0];
                        bunch[p].coordonnees[1][1] = bunch[p].coordonnees[0][1];
                        bunch[p].coordonnees[1][2] = bunch[p].coordonnees[0][2];
                        bunch[p].coordonnees[1][3] = bunch[p].coordonnees[0][3];
                        bunch[p].coordonnees[1][4] = bunch[p].coordonnees[0][4];
                        bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5];

                    } else if ((reseau[i]->K1L != 0) && (nonlinflag == 1)) {

                        double R11Xh, R12Xh, R21Xh, R22Xh, R11Yh, R12Yh, R21Yh, R22Yh;


                        R11Xh = R11X[i] + R11XC[i] * dpopeff;
                        R12Xh = R12X[i] + R12XC[i] * dpopeff;
                        R21Xh = R21X[i] + R21XC[i] * dpopeff;
                        R22Xh = R22X[i] + R22XC[i] * dpopeff;
                        R11Yh = R11Y[i] + R11YC[i] * dpopeff;
                        R12Yh = R12Y[i] + R12YC[i] * dpopeff;
                        R21Yh = R21Y[i] + R21YC[i] * dpopeff;
                        R22Yh = R22Y[i] + R22YC[i] * dpopeff;

                        bunch[p].coordonnees[1][0] = R11Xh * bunch[p].coordonnees[0][0] + R12Xh * bunch[p].coordonnees[0][1] + (reseau[i]->DX - R11Xh * reseau[i - 1]->DX - R12Xh * reseau[i - 1]->DPX) * dpopeff;
                        bunch[p].coordonnees[1][1] = R21Xh * bunch[p].coordonnees[0][0] + R22Xh * bunch[p].coordonnees[0][1] + (reseau[i]->DPX - R21Xh * reseau[i - 1]->DX - R22Xh * reseau[i - 1]->DPX) * dpopeff;
                        bunch[p].coordonnees[1][2] = R11Yh * bunch[p].coordonnees[0][2] + R12Yh * bunch[p].coordonnees[0][3] + (reseau[i]->DY - R11Yh * reseau[i - 1]->DY - R12Yh * reseau[i - 1]->DPY) * dpopeff;
                        bunch[p].coordonnees[1][3] = R21Yh * bunch[p].coordonnees[0][2] + R22Yh * bunch[p].coordonnees[0][3] + (reseau[i]->DPY - R21Yh * reseau[i - 1]->DY - R22Yh * reseau[i - 1]->DPY) * dpopeff;
                        bunch[p].coordonnees[1][4] = bunch[p].coordonnees[0][4];

                        bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5] + time(bunch[p], reseau[i]->L, betgam, i);

                    } else if ((reseau[i]->K2L != 0) && (nonlinflag == 1)) {
                        double Lelemha(Lelem[i] / 2);
                        double dxha(reseau[i - 1]->DX + reseau[i - 1]->DPX * Lelemha);
                        double dyha(reseau[i - 1]->DY + reseau[i - 1]->DPY * Lelemha);

                        bunch[p].coordonnees[1][0] = bunch[p].coordonnees[0][0] + Lelemha * bunch[p].coordonnees[0][1] + (dxha - reseau[i - 1]->DX - Lelemha * reseau[i - 1]->DPX) * dpopeff;
                        bunch[p].coordonnees[1][2] = bunch[p].coordonnees[0][2] + Lelemha * bunch[p].coordonnees[0][3] + (dyha - reseau[i - 1]->DY - Lelemha * reseau[i - 1]->DPY) * dpopeff;

                        bunch[p].coordonnees[1][1] = bunch[p].coordonnees[0][1] - 0.5 * reseau[i]->K2L * (bunch[p].coordonnees[1][0] * bunch[p].coordonnees[1][0] - bunch[p].coordonnees[1][2] * bunch[p].coordonnees[1][2]);
                        bunch[p].coordonnees[1][3] = bunch[p].coordonnees[0][3] + reseau[i]->K2L * (bunch[p].coordonnees[1][0] * bunch[p].coordonnees[1][2]);
                        bunch[p].coordonnees[1][0] = bunch[p].coordonnees[1][0] + Lelemha * bunch[p].coordonnees[1][1] + (reseau[i]->DX - dxha - Lelemha * reseau[i - 1]->DPX) * dpopeff;
                        bunch[p].coordonnees[1][2] = bunch[p].coordonnees[1][2] + Lelemha * bunch[p].coordonnees[1][3] + (reseau[i]->DY - dyha - Lelemha * reseau[i - 1]->DPY) * dpopeff;
                        bunch[p].coordonnees[1][4] = bunch[p].coordonnees[0][4];

                        bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5] + time(bunch[p], reseau[i]->L, betgam, i);
                    } else {

                        bunch[p].coordonnees[1][0] = R11X[i] * bunch[p].coordonnees[0][0] + R12X[i] * bunch[p].coordonnees[0][1] + (reseau[i]->DX - R11X[i] * reseau[i - 1]->DX - R12X[i] * reseau[i - 1]->DPX) * dpopeff;
                        bunch[p].coordonnees[1][1] = R21X[i] * bunch[p].coordonnees[0][0] + R22X[i] * bunch[p].coordonnees[0][1] + (reseau[i]->DPX - R21X[i] * reseau[i - 1]->DX - R22X[i] * reseau[i - 1]->DPX) * dpopeff;
                        bunch[p].coordonnees[1][2] = R11Y[i] * bunch[p].coordonnees[0][2] + R12Y[i] * bunch[p].coordonnees[0][3] + (reseau[i]->DY - R11Y[i] * reseau[i - 1]->DY - R12Y[i] * reseau[i - 1]->DPY) * dpopeff;
                        bunch[p].coordonnees[1][3] = R21Y[i] * bunch[p].coordonnees[0][2] + R22Y[i] * bunch[p].coordonnees[0][3] + (reseau[i]->DPY - R21Y[i] * reseau[i - 1]->DY - R22Y[i] * reseau[i - 1]->DPY) * dpopeff;

                        if ((reseau[i]->KEYWORD != "RFCAVITY") && (RFflag == 1)) {
                            bunch[p].coordonnees[1][4] = bunch[p].coordonnees[0][4];
                        }

                        bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5] + time(bunch[p], reseau[i]->L, betgam, i);

                    }

                }

                if ((nonlinflag == 1) && ((icop == -1) || (resColli[icop]->method == "standard") || (resColli[icop]->method == "magnetic")) && (bunch[p].inabs != 0)) {

                    int ixcor(-1), iycor (-1);

                    for (int k(0); k < ixcormag.size(); ++k) {
                        if (ixcormag[k] == i) {
                            ixcor = k;
                        }
                    }

                    if (ixcor != -1) {
                        bunch[p].coordonnees[1][0] = bunch[p].coordonnees[1][0] + scaleorbit * (reseau[i]->S - reseau[i - 1]->S) / 2 * xcormag[ixcor] * (1 + dpopeff);
                        bunch[p].coordonnees[1][1] = bunch[p].coordonnees[1][1] + scaleorbit * xcormag[ixcor] * (1 + dpopeff);
                    }

                    for (int k(0); k < iycormag.size(); ++k) {
                        if (iycormag[k] == i) {
                            iycor = k;
                        }
                    }

                    if (iycor != -1) {
                        bunch[p].coordonnees[1][2] = bunch[p].coordonnees[1][2] + scaleorbit * (reseau[i]->S - reseau[i - 1]->S) / 2 * ycormag[iycor] * (1 + dpopeff);
                        bunch[p].coordonnees[1][3] = bunch[p].coordonnees[1][3] + scaleorbit * ycormag[iycor] * (1 + dpopeff);
                    }

                }



                //checking for aperture hits
                if (bunch[p].inabs == 1) {

                    if ((icop == -1) || (resColli[icop]->method == "standard") || (resColli[icop]->method == "magnetic") || (resColli[icop]->method == "crystal")) {

                        bool inside(false);

                        if (reseau[i]->APERTYPE == "RECTANGLE") {
                            if ((abs(bunch[p].coordonnees[1][0]) < reseau[i]->aperx) && (abs(bunch[p].coordonnees[1][2]) < reseau[i]->apery) && (abs(bunch[p].coordonnees[0][0]) < reseau[i]->aperx) && (abs(bunch[p].coordonnees[0][2]) < reseau[i]->apery)) {
                                inside = true;
                            }
                        } else if ((reseau[i]->APERTYPE == "ELLIPSE") || (reseau[i]->APERTYPE == "CIRCLE") || (reseau[i]->APERTYPE == "RECTELLIPSE")) {
                            if ((((bunch[p].coordonnees[1][0] / reseau[i]->aperx) * (bunch[p].coordonnees[1][0] / reseau[i]->aperx) + (bunch[p].coordonnees[1][2] / reseau[i]->apery) * (bunch[p].coordonnees[1][2] / reseau[i]->apery)) < 1) && (((bunch[p].coordonnees[0][0] / reseau[i]->aperx) * (bunch[p].coordonnees[0][0] / reseau[i]->aperx) + (bunch[p].coordonnees[0][2] / reseau[i]->apery) * (bunch[p].coordonnees[0][2] / reseau[i]->apery)) < 1)) {
                                inside = true;

                            }
                        } else {
                            cerr << "Error for the " << i << "th element: unknown aperture type!" << endl;
                        }

                        double L;

                        L = reseau[i]->S - reseau[i - 1]->S; //distance to next element in the accelerator

                        double Ap0h, Zp0h, xh0, xh, xhs, yh0, yh, yhs;
                        int nh(0);//number of lost particles

                        if (inside == false) {
                            cout <<"The particle is lost at element " << reseau[i]->NAME << endl;
                            xh0 = bunch[p].coordonnees[0][0];
                            bunch[p].inabs = 0;

                            if (L == 0) {
                                nh = nh + 1;
                                hits.push_back(reseau[i]->S);//saving the value of s where the particles get lost
                                Aphit.push_back(bunch[p].Ap0);//saving the mass of the lost particles
                                Zphit.push_back(bunch[p].Zp0);//saving the charge of the lost particles
                            } else {

                                xh = bunch[p].coordonnees[1][0];
                                xhs = (xh - xh0) / L;

                                yh0 = bunch[p].coordonnees[0][2];
                                yh = bunch[p].coordonnees[1][2];
                                yhs = (yh - yh0) / L;

                                Aphit.push_back(bunch[p].Ap0);
                                Zphit.push_back(bunch[p].Zp0);

                                double slostx, slosty, splus, sminus;

                                if (reseau[i]->APERTYPE == "RECTANGLE") {

                                    //find hit position in x

                                    if ((abs(xh0) < reseau[i]->aperx) && (abs(xh) < reseau[i]->aperx)) { //both points inside aperture - particle lost in y instead
                                        slostx = reseau[i]->S;

                                    } else if (abs(xh0) > reseau[i]->aperx) { //first point outside - particle lost already at the entrance of the element
                                        slostx = reseau[i - 1]->S;

                                    } else {//first point inside, second point outside - do a linear interpolation

                                        splus = (reseau[i]->aperx - xh0) / xhs;
                                        sminus = (-reseau[i]->aperx - xh0) / xhs;

                                        if (splus > sminus) { //choose the correct point where the straight line hits the aperture, at +-a. The highest s-value is correct
                                            slostx = splus + reseau[i - 1]->S;

                                        } else {
                                            slostx = sminus + reseau[i - 1]->S;

                                        }
                                    }

                                    //find hit position in y

                                    if ((abs(yh0) < reseau[i]->apery) && (abs(yh) < reseau[i]->apery)) { //both points inside aperture - particle lost in x instead
                                        slosty = reseau[i]->S;
                                    } else if (abs(yh0) > reseau[i]->apery) { //first point outside - particle lost already at the entrance of the element
                                        slosty = reseau[i - 1]->S;
                                    } else {//first point inside, second point outside - do a linear interpolation

                                        splus = (reseau[i]->apery - yh0) / yhs;
                                        sminus = (-reseau[i]->apery - yh0) / yhs;

                                        if (splus > sminus) { //choose the right point where the straight line hits the aperture, at +-b. The highest s-value is correct
                                            slosty = splus + reseau[i - 1]->S;
                                        } else {
                                            slosty = sminus + reseau[i - 1]->S;
                                        }
                                    }

                                    //choose the point where the particle hits first, x or y

                                    if (slostx < slosty) {
                                        hits.push_back(slostx);
                                    } else {
                                        hits.push_back(slosty);
                                    }
                                } else if (reseau[i]->APERTYPE == "ELLIPSE") {

                                    double sqrarg;
                                    sqrarg = reseau[i]->apery * reseau[i]->apery * xhs * xhs - xhs * xhs * yh0 * yh0 + 2 * xh0 * xhs * yh0 * yhs + reseau[i]->aperx * reseau[i]->aperx * yhs * yhs - xh0 * xh0 * yhs * yhs;
                                    if (sqrarg < 0) {
                                        cout << "Warning, sqrarg < 0" << endl;
                                        hits.push_back(reseau[i - 1]->S);
                                    } else {
                                        double stry;

                                        stry = (-reseau[i]->apery * reseau[i]->apery * xh0 * xhs - reseau[i]->aperx * reseau[i]->aperx * yh0 * yhs - reseau[i]->aperx * reseau[i]->apery * sqrt(sqrarg)) / (reseau[i]->apery * reseau[i]->apery * xhs * xhs + reseau[i]->aperx * reseau[i]->aperx * yhs * yhs);

                                        if ((stry > 0) && (stry < L)) {
                                            hits.push_back(reseau[i - 1]->S + stry);
                                        } else {

                                            double stry2;

                                            stry2 = (-reseau[i]->apery * reseau[i]->apery * xh0 * xhs - reseau[i]->aperx * reseau[i]->aperx * yh0 * yhs + reseau[i]->aperx * reseau[i]->aperx * sqrt(sqrarg)) / (reseau[i]->apery * reseau[i]->apery * xhs * xhs + reseau[i]->aperx * reseau[i]->aperx * yhs * yhs);

                                            if ((stry2 > 0) && (stry2 < L)) {
                                                hits.push_back(reseau[i - 1]->S + stry2);
                                            } else {
                                                hits.push_back(reseau[i - 1]->S);
                                            }
                                        }
                                    }
                                }//end if ELLIPSE
                                else if (reseau[i]->APERTYPE == "RECTELLIPSE") {

                                    double sqrarg;
                                    sqrarg = reseau[i]->apery * reseau[i]->apery * xhs * xhs - xhs * xhs * yh0 * yh0 + 2 * xh0 * xhs * yh0 * yhs + reseau[i]->aperx * reseau[i]->aperx * yhs * yhs - xh0 * xh0 * yhs * yhs;
                                    if (sqrarg < 0) {
                                        cout << "Warning, sqrarg < 0" << endl;
                                        hits.push_back(reseau[i - 1]->S);
                                    } else {
                                        double stry;
                                        stry = (-reseau[i]->apery * reseau[i]->apery * xh0 * xhs - reseau[i]->aperx * reseau[i]->aperx * yh0 * yhs - reseau[i]->aperx * reseau[i]->apery * sqrt(sqrarg)) / (reseau[i]->apery * reseau[i]->apery * xhs * xhs + reseau[i]->aperx * reseau[i]->aperx * yhs * yhs);

                                        if ((stry > 0) && (stry < L)) {
                                            hits.push_back(reseau[i - 1]->S + stry);
                                        } else {
                                            double stry2;

                                            stry2 = (-reseau[i]->apery * reseau[i]->apery * xh0 * xhs - reseau[i]->aperx * reseau[i]->aperx * yh0 * yhs + reseau[i]->aperx * reseau[i]->apery * sqrt(sqrarg)) / (reseau[i]->apery * reseau[i]->apery * xhs * xhs + reseau[i]->aperx * reseau[i]->aperx * yhs * yhs);

                                            if ((stry2 > 0) && (stry2 < L)) {
                                                hits.push_back(reseau[i - 1]->S + stry2);
                                            } else {
                                                hits.push_back(reseau[i - 1]->S);
                                            }
                                        }
                                    }
                                }//end RECTELLIPSE
                            }//end else after if L == 0

                        }//end si particle lost
                        else {
                            //cout << "No losses, all the particles stay in the experience." << endl;
                        }
                    }//end if icop==vide ||...

                }
            }

            //Different ways to print the coordinates of the particle through the parameters IDPART and OUTCOORD (see manual)

            if ((outcoord != 0) && (bunch[p].inabs != 0)) {

                if (outcoord == 1) {

                    cout << "Particle " <<  p + 1 << " after the passage through element: " << i << " " << reseau[i]->NAME << endl;
                    bunch[p].afficheCoordonnees();
                    cout << "Turn number " << turn << endl;
                    cout << "Massnumber: " << bunch[p].Ap0 << " Chargestate: " << bunch[p].Zp0 << endl;
                    cout << "Particle ID: " << bunch[p].getidentification() << endl;
                } else if (outcoord == 3) {

                    if (bunch[p].getidentification() == choicePart) {
                        outCoord(bunch[p], i, outputpath + "/coordinates.dat");
                    }
                } else if (outcoord == 4) {

                    if (bunch[p].getidentification() == choicePart) {
                        outPunct(i, bunch[p], bunch[0], attr1, outputpath);
                    }
                }

                else if (outcoord == 2) {
                    outElt(i, bunch[p], outputpath, indication);
                }
            }


        }//end loop over particles


        //we prepare the particles for the next element
        for (int b(0); b < bunch.size(); ++b) {
            if (bunch[b].inabs != 0) {
                for (int r(0); r < 6; ++r) {
                    bunch[b].coordonnees[0][r] = bunch[b].coordonnees[1][r];
                }
            }
        }

        if (plotflag == "Yes") {

            vector <double> temp1, temp2;

            for (int h(0); h < bunch.size(); ++h) {
                if (bunch[h].inabs != 0) {
                    temp1.push_back(bunch[h].coordonnees[0][0]);
                    temp2.push_back(bunch[h].coordonnees[0][2]);
                } else {
                    temp1.push_back(100);
                    temp2.push_back(100);
                }
            }

            xco.push_back(temp1);
            yco.push_back(temp2);

            temp1.clear();
            temp2.clear();
        }


        //if there are no more particles, we return
        if (bunch.size() == 0) {
            return;
        }

    }//end loop elements

    vector <Particle> tempinabs;

    //we contine just with particles that are not lost
    for (int w(0); w < bunch.size(); ++w) {
        if (bunch[w].inabs != 0) {
            tempinabs.push_back(bunch[w]);
        }
    }

    bunch.clear();

    for (int w(0); w < tempinabs.size(); ++w) {
        bunch.push_back(tempinabs[w]);
    }

    tempinabs.clear();

}