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

#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
#include <string>
#include <cmath>
#include "CrystalCollimator.h"
using namespace std;

int CrystalCollimator::crystaloutflag = 0;


CrystalCollimator::CrystalCollimator(const double& ALFX, const double& ALFY, const double& APER_1, const double& APER_2, const double& APER_3, const double& APER_4, const string& APERTYPE, const double& BETX, const double& BETY, const double& DPX, const double& DPY, const double& DX, const double& DY, const string& KEYWORD, const double& L, const double& MUX, const double& MUY, const string& NAME, const double& PTC, const double& PXC, const double& PYC, const double& S, const double& TC, const double& XC, const double& YC, const double& K0L, const double& K0SL, const double& K1L, const double& K1SL, const double& K2L, const double& K2SL, const string& PARENT, const string& meth, const double& hgap, const double& hgap2, const double& collang, const double& pdepth, const double& pdepth2, const double& tcang, const double& nsig)
    : Collimator(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, meth, hgap, hgap2, collang, pdepth, pdepth2, tcang, nsig)
{};

CrystalCollimator::CrystalCollimator(Element elt, const double& tcang, const double& nsig, const string& meth)
    : Collimator(elt, tcang, nsig, meth)
{};

CrystalCollimator::CrystalCollimator(const CrystalCollimator& obj)
    : Collimator(obj)
{};

void CrystalCollimator::affiche()
{

    cout << "Crystal collimator: " << endl;

    this->Collimator::affiche();
};


void CrystalCollimator::collipassCrystal(vector <Particle>& bunch, const double& betgam, const int& pas, string outputpath)
{

    //Calculate proton energy
    double Mproton(0.93827231); //rest mass [GeV]

    vector <double> betgamProton;
    vector <double> betaProton;
    vector <double> Eproton;

    for (int j(0); j < bunch.size(); ++j) {
        betgamProton.push_back(betgam * (1 + bunch[j].coordonnees[0][4])); //relativistic beta * relativistic gamma
        betaProton.push_back(sqrt((betgamProton[j]*betgamProton[j]) / ((betgamProton[j]*betgamProton[j]) + 1))); //relativistic beta
        Eproton.push_back(Mproton * betgamProton[j] / betaProton[j]); //total energy [GeV]
    }


    ofstream sortie;
    string filename;
    filename = outputpath + "/ascii_before.csv";

    sortie.open(filename.c_str());

    if (sortie.fail()) {
        cerr << "Warning: problem with the file " << filename << "!" << endl;
    } else {

        int col(15);

        for (int k(0); k < bunch.size(); ++k) {

            sortie  << bunch[k].coordonnees[0][0] << "," << bunch[k].coordonnees[0][2] << "," << bunch[k].coordonnees[0][1] << "," << bunch[k].coordonnees[0][3] << "," << Eproton[k] << endl;
        }

        sortie.close();

    }

    ofstream out;
    string file;
    file = outputpath + "/crystal_output.out";

    if (crystaloutflag == 0) {
        out.open(file.c_str());
        crystaloutflag = 1;
    } else {
        out.open(file.c_str(), ios::out | ios::app);
    }

    if (out.fail()) {
        cerr << "Warning: problem writing in the file: " << file << "!" << endl;
    } else {
        out << this->NAME << endl;
        out.close();
    }


    SimCrys sim(Crystal (C_orient,  IS, C_xmax, C_ymax, Cry_length, Rcurv), Partcrys (0)); //simulation of the passage through the crystal

    sim.general(sim, pas, outputpath, C_rotation, C_aperture, C_offset, C_tilt, Cry_tilt);

    ifstream enter;
    string file_in;
    file_in = outputpath + "/ascii_after.csv";

    enter.open(file_in.c_str());

    if (enter.fail()) {
        cerr << "Warning: problem reading the file " << file_in << "!" << endl;
    } else {

        for (int k(0); k < bunch.size(); ++k) {

            string rest;

            getline(enter, rest, ',');
            bunch[k].coordonnees[1][0] = atof(rest.c_str());

            getline(enter, rest, ',');
            bunch[k].coordonnees[1][2] = atof(rest.c_str());

            getline(enter, rest, ',');
            bunch[k].coordonnees[1][1] = atof(rest.c_str());

            getline(enter, rest, ',');
            bunch[k].coordonnees[1][3] = atof(rest.c_str());

            getline(enter, rest);
            Eproton[k] = atof(rest.c_str());
        }

        enter.close();
    }

    vector <double> gamProton;
    vector <double> betProton;

    //we return to the initial referential
    for (int j(0); j < Eproton.size(); ++j) {
        gamProton.push_back(Eproton[j] / Mproton);
        betProton.push_back(sqrt(1 - (1 / (gamProton[j]*gamProton[j]))));
        bunch[j].coordonnees[1][4] = gamProton[j] * betProton[j] / betgam - 1;
    }

};