#include <iostream>
#include <fstream>
#include <cmath>
#include <iomanip>
#include <string.h>
#include <stdio.h>
#include "simcrys.h"

//crystal routine of the proton beam 
#include "crystalproton.h"
#include "crystalprotonfuns.h"

extern "C"
{
   struct {
  double Cry_length, Rcurv, C_xmax, C_ymax, Alayer, C_orient;   
  }Par_Cry1_;
};


using namespace std;

//Fortran crystal routine for the proton beam
// By Jianfeng Zhang @ LAL, 05/2014
// extern "C" {
// // void cryst_(int,double,double,double,double,double,double);
// void cryst_(int IS, double x, double xp, double y,double yp,double PC,double Length);  
// }



SimCrys::SimCrys(Crystal crys, Partcrys part)
    : crys(crys), part(part) //,moy(moy)                   //POUR faire varier la moyenne de la gaussienne aussi !!!!!
{
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////// FUNCTION BASES THAT MAKE SOME BASIC CALCULATION
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

void SimCrys::bases(SimCrys& sim)
{
    if ( (crys.miscut < 0 ) &&  ( part.x > 0) && (part.x < -crys.Cry_length * tan(crys.miscut))) {
        crys.L_chan = -part.x * tan(crys.miscut);
    }
    crys.tchan = crys.L_chan / crys.Rcurv;
    part.xp_rel = part.xp - crys.miscut;
}



////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////FUNCTIN CROSS TO KNOW IF THE PARTICULE CROSS THE CRYSTAL
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////




bool SimCrys::cross(SimCrys& sim)
{
    bool crossing(true);
    if ((part.y < crys.ymin) || (part.y > crys.ymax) || (part.x > crys.xmax)) {
        //cout << "OUT of the Crystal !!" << part.x << " , "<<part.xp << endl;
        crossing = false;
        part.x = part.x + part.xp * crys.s_length;

        part.y = part.y + part.yp * crys.s_length;

        part.nout = part.nout + 1;     /////////////COUNT/////////////
        part.effet = 1;
        //cout<<" BLAAAAAAAAa"<<endl;

    }

    return crossing;
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////FUNCTION LAYER THAT DETERMIN IF THE PARTICULE IS IN THE AMORPHOUS LAYER/////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

bool SimCrys::layer(SimCrys& sim)
{
    //cout<<"La fonction layer est utilisee>>>>>>>>>>>>>>"<<endl;
    bool out(true);
    if ((part.x < crys.Alayer) || ((part.y - crys.ymin) < crys.Alayer) || ((crys.ymax - part.y) < crys.Alayer)) {
        cout<<"Amorphous layer! " << part.x <<" , " <<part.xp<<endl;
        double a_eq;
        double b_eq;
        double c_eq;
        double delta;         //a, b, c, Delta of the second order eq.
        double x0;
        double y0;
        double length_xs;     //!Amorphous length
        double length_ys;     //!Amorphous length
        double am_length;     //!Amorphous length
        out = false;
        x0 = part.x;
        y0 = part.y;
        a_eq = (1 + part.xp * part.xp);
        b_eq = (2 * (part.x) * (part.xp) - 2 * (part.xp) * crys.Rcurv);
        c_eq = (part.x) * part.x - 2 * abs(part.x) * crys.Rcurv;
        delta = b_eq * b_eq - 4 * a_eq * c_eq;
        //cout<<"a : " <<  a_eq<<endl;      //POUR TESTER FONCTION
        //cout<<"b: " << b_eq <<endl;      //POUR TESTER FONCTION
        //cout<<"c : " << c_eq <<endl;      //POUR TESTER FONCTION
        //cout<<"Delta : " << delta <<endl;      //POUR TESTER FONCTION
        part.s = (-b_eq + sqrt(delta)) / ((2 * a_eq)) ; //in [m]
        //cout<<"s : " << part.s <<endl;      //POUR TESTER FONCTION
        if (part.s >= crys.s_length) {
            part.s = crys.s_length;
        }
        //cout<<"s : " << part.s <<endl;      //POUR TESTER FONCTION
        part.x = part.xp * part.s + x0   ; //in [m]
        //cout<<"x0 : " << x0 <<endl;      //POUR TESTER FONCTION
        //cout<<"x : " << part.x <<endl;      //POUR TESTER FONCTION

        length_xs = sqrt((part.x - x0) * (part.x - x0) + part.s * part.s);
        //cout<<"Length_XS : " << length_xs <<endl;      //POUR TESTER FONCTION

        if ((((part.yp >= 0) && ((part.y + part.yp * part.s) <= crys.ymax))) || (((part.yp < 0) && ((part.y + part.yp * part.s) >= crys.ymin)))) {
            length_ys = part.yp * length_xs;
        } else {
            part.s = (crys.ymax - part.y) / part.yp;
            length_ys = sqrt((crys.ymax - part.y) * (crys.ymax - part.y) + part.s * part.s);
            part.x = x0 + part.xp * part.s;
            length_xs = sqrt((part.x - x0) * (part.x - x0) + part.s * part.s);
        }

        am_length = sqrt(length_xs * length_xs + length_ys * length_ys);
        //cout << "AM_Length :" << am_length << endl;  // TESTER FONCTION
        part.s = part.s / 2;
        part.x = x0 + part.xp * part.s;
        part.y = y0 + part.yp * part.s;
        //cout<<"  am_length      "<<am_length<<endl;
        move_am(crys.IS, crys.NAM, am_length, crys.DES[crys.IS], crys.DLYI[crys.IS], crys.DLRI[crys.IS], part.xp, part.yp, part.PC); //We call the move_am function
        //cout<<"AM LAYER !!!!       "<<part.PC<<endl;
        //cout<< "Amorphous"<< endl;
        part.x = part.x + part.xp * (crys.Cry_length - part.s);
        part.y = part.y + part.yp * (crys.Cry_length - part.s);
        part.namor = part.namor + 1;     /////////////COUNT/////////////
        part.effet = 2;
    } else if ((part.x > (crys.xmax - crys.Alayer)) && (part.x < crys.xmax)) {
        out = false;

        move_am(crys.IS, crys.NAM, crys.s_length, crys.DES[crys.IS], crys.DLYI[crys.IS], crys.DLRI[crys.IS], part.xp, part.yp, part.PC); //We call the move_am function

        part.namor = part.namor + 1;     /////////////COUNT/////////////
        part.effet = 2;
        //cout<<"Fix HERE"<<endl;
    }
    return out;
}


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//SI LES DEUX FONCTIONS D AVT RETURNE TRUE, ALORS ON CALCULE LES PARA CRITIQUE ET D AUTRE PARA (POINT C - E)
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////



void SimCrys::parameters(SimCrys& sim)
{
    //THE FIRST PART IS FOR THE (110) ORIENTATION
    //cout << "The parameters function is use !!!!! "<< endl;
    double c_v1;
    double c_v2;


    crys.xpcrit0 = sqrt(2.10e-9 * crys.eUm[crys.IS] / part.PC);
    crys.Rcrit  = part.PC / (2.e-6 * crys.eUm[crys.IS]) * crys.AI[crys.IS]; //if R>Rcritical=>no channeling is possible
    crys.ratio = crys.Rcurv / crys.Rcrit;
    crys.xpcrit = crys.xpcrit0 * (crys.Rcurv - crys.Rcrit) / crys.Rcurv;
    c_v1 = 1.7;
    c_v2 = -1.5;
    if (crys.ratio <= 1) {
        part.Ang_rms = c_v1 * 0.42 * crys.xpcrit0 * sin(1.4 * crys.ratio); //rms scattering
        part.Ang_avr = c_v2 * crys.xpcrit0 * 0.05 * crys.ratio;    //average angle reflection
        part.Vcapt = 0.0;
    } else if (crys.ratio <= 3) {
        part.Ang_rms = c_v1 * 0.42 * crys.xpcrit0 * sin(1.571 * 0.3 * crys.ratio + 0.85); //rms scattering
        part.Ang_avr = c_v2 * crys.xpcrit0 * (0.1972 * crys.ratio - 0.1472);       //average angle reflection
        part.Vcapt   = 0.0007 * (crys.ratio - 0.7) / pow(part.PC, 0.2);           //K=0.00070 is taken based on simulations using CATCH.f (V.Biryukov)
    } else {
        part.Ang_rms = c_v1 * crys.xpcrit0 * (1 / crys.ratio);
        part.Ang_avr = c_v2 * crys.xpcrit0 * (1 - 1.6667 / crys.ratio);
        part.Vcapt = 0.0007 * (crys.ratio - 0.7) / pow(part.PC, 0.2);
    }
    //cout <<crys.xpcrit<<endl;
    if (crys.C_orient == 2) {
        part.Ang_avr = part.Ang_avr * 0.93;
        part.Ang_rms = part.Ang_rms * 1.05;                 // FOR (111) ORIENTATION
        crys.xpcrit  = crys.xpcrit * 0.98;
    }
    //cout <<crys.xpcrit<<endl;
}




////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//MAINTENANT ON FAIT UNE FONCTION POUR LA PARTIE I CF SHEMA!!! C EST LA PARTIE CHANNELING D OU LE NOM
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////




void SimCrys::channel(SimCrys& sim)
{
    //cout << "La fonction channel est utilisee!!!! "<<endl;


    part.Chann  = sqrt(crys.xpcrit * crys.xpcrit - part.xp_rel * part.xp_rel) / crys.xpcrit0; //probability of channeling/volume capture

    //cout <<"Channeling Prob :" << part.Chann <<endl;   //POUR REGARDER LA VALEUR!!!!!!!!!!!!!!
    //cout <<"Channeling ANGLE :" << crys.tchan<<endl;   //POUR REGARDER LA VALEUR!!!!!!!!!!!!!!


    double n_atom(0.1);                                //probability for entering channeling near atomic planes
    double Dxp(0);					//change angle from channeling [mrad]
    double TLdech1(0);       //tipical dechanneling length(1) [m]
    double tdech(0);
    double Ldech(0);	      //Dechanneling. length
    double Sdech(0);
    //cout<<"Chann     : "<<part.Chann<<endl;
    if (random() <= part.Chann) {
        TLdech1 = 0.0005 * part.PC * pow(1 - 1 / crys.ratio, 2); // calculate tipical dechanneling length(m)
        if (random() <= n_atom) {
            TLdech1 = 0.000002 * part.PC * pow(1. - 1. / crys.ratio, 2); // dechanneling length (m) for short crystal for high amplitude particles
        }
        //cout<<"@ energy " << part.PC <<" dechanneling length "<< TLdech1 << endl; //FAUDRA OTER CELA!!!!!!!!!!
        part.Dechan = log(random());
        //cout<<"dechannn " << part.Dechan <<endl; //FAUDRA OTER CELA!!!!!!!!!!
        Ldech = -TLdech1 * part.Dechan; //careful: the dechanneling lentgh is along the trajectory of the particle -not along the longitudinal coordinate...
        tdech = Ldech / crys.Rcurv;
        Sdech = Ldech * cos(crys.miscut + 0.5 * tdech);
        // 2 option if the particule channeling !!!
        //      (1) dechanneling
        //      (2) full channeling

        //    (1) :

        if (Ldech < crys.L_chan) {
            //cout <<"Dechanneling ! "<< endl;
            Dxp = Ldech / crys.Rcurv;    //change angle from channeling [mrad]

            part.ndechann = part.ndechann + 1;     /////////////COUNT/////////////
            part.effet = 4;
            part.x  = part.x + Ldech * (sin(0.5 * Dxp + crys.miscut)); // trajectory at channeling exit
            part.xp = part.xp + Dxp + (random_gauss() - 0.5) * crys.xpcrit * 0.5;
            part.y = part.y + part.yp * Sdech;

            part.x = part.x + 0.5 * (crys.s_length - Sdech) * part.xp;
            part.y = part.y + 0.5 * (crys.s_length - Sdech) * part.yp;


            move_am(crys.IS, crys.NAM, crys.s_length - Sdech, crys.DES[crys.IS], crys.DLYI[crys.IS], crys.DLRI[crys.IS], part.xp, part.yp, part.PC); //We call the move_am function

            part.PC = part.PC - 0.5 * crys.DES[crys.IS] * part.y; //energy loss to ionization [GeV]
            part.x = part.x + 0.5 * (crys.s_length - Sdech) * part.xp;
            part.y = part.y + 0.5 * (crys.s_length - Sdech) * part.yp;
        }

        //      (2) :

        else {
            //cout <<"Channeling ! "<< endl;
            Dxp = crys.L_chan / crys.Rcurv;                             //change angle [mrad]
            part.x  = part.x + crys.L_chan * (sin(0.5 * Dxp + crys.miscut)); //trajectory at channeling exit


            // Removed dependence on incoming angle according to discussion with I.Yazynin on 24/10/2012
            //part.xp = part.xp + Dxp + (random_gauss()-0.5)*crys.xpcrit*0.5;        // [mrad]
            part.xp = /*part.xp*/ + Dxp + (random_gauss() - 0.5) * crys.xpcrit * 0.5;  // [mrad]

            part.y = part.y + crys.s_length * part.yp;
            part.PC = part.PC - 0.5 * crys.DES[crys.IS] * crys.Cry_length;    // energy loss to ionization [GeV]
            part.nchann = part.nchann + 1;     /////////////COUNT/////////////
            part.effet = 3;
        }
    } else {
        //cout <<"Volume Reflection ! "<< endl;
        part.nvrefl = part.nvrefl + 1;     /////////////COUNT/////////////
        part.effet = 5;
        Dxp = 0.5 * (part.xp_rel / crys.xpcrit + 1) * part.Ang_avr;
        part.xp = part.xp + Dxp + part.Ang_rms * random_gauss();
        /*
        next line new from sasha
        part.xp=part.xp+0.45*(part.xp/crys.xpcrit+1)*part.Ang_avr;
        valentina fix here
        */
        part.x = part.x + 0.5 * crys.s_length * part.xp;
        part.y = part.y + 0.5 * crys.s_length * part.yp;

        move_am(crys.IS, crys.NAM, crys.s_length, crys.DES[crys.IS], crys.DLYI[crys.IS], crys.DLRI[crys.IS], part.xp , part.yp, part.PC); //We call the move_am function

        part.x = part.x + 0.5 * crys.s_length * part.xp;
        part.y = part.y + 0.5 * crys.s_length * part.yp;

    }
}



////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//MAINTENANT ON FAIT UNE FONCTION POUR LA PARTIE II CF SHEMA!!! C EST LA PARTIE VOLUME REFLEXION D OU LE NOM
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////




void SimCrys::reflection(SimCrys& sim)
{
    //cout<<"The function reflection is used" <<endl;
    double Dxp(0);		//change angle from channeling [mrad]
    double Lrefl(0);         //distance of the reflection point [m]
    double Srefl(0);         //distance of the reflection point [m]
    double TLdech2(0);       //tipical dechanneling length(2) [m]
    double tdech(0);
    double Ldech(0);	      //Dechanneling. length
    double Sdech(0);
    double Red_S(0);         //Reduced length (in case of dechanneling)
    double Rlength(0);       //Reduced length

    Lrefl =  (part.xp_rel) * crys.Rcurv;
    Srefl = sin(part.xp) * Lrefl;


    if ((Lrefl > 0) && (Lrefl < crys.Cry_length)) {   // IE : If the VR point is inside the crystal!!!!

        //cout<<"VRCapt     : "<<part.Vcapt<<endl;
        if ((random() > part.Vcapt) || (crys.ZN == 0)) {
            //cout<< "Volume Reflexion! "<<endl;
            part.nvrefl = part.nvrefl + 1;     /////////////COUNT/////////////
            part.effet = 5;
            part.x = part.x + part.xp * Srefl;
            part.y = part.y + part.yp * Srefl;
            Dxp = part.Ang_avr;
            part.xp = part.xp + Dxp + part.Ang_rms * random_gauss();
            part.x = part.x + 0.5 * part.xp * (crys.s_length - Srefl);
            part.y = part.y + 0.5 * part.yp * (crys.s_length - Srefl);

            if (crys.ZN > 0) {
                move_am(crys.IS, crys.NAM, crys.s_length - Srefl, crys.DES[crys.IS], crys.DLYI[crys.IS], crys.DLRI[crys.IS], part.xp , part.yp, part.PC); //We call the move_am function
            }

            part.x = part.x + 0.5 * part.xp * (crys.s_length - Srefl);
            part.y = part.y + 0.5 * part.yp * (crys.s_length - Srefl);
        } else {
            //cout<< "Volume Capture! ";
            part.nvcapt = part.nvcapt + 1;     /////////////COUNT/////////////
            part.effet = 6;
            part.x = part.x + part.xp * Srefl;
            part.y = part.y + part.yp * Srefl;
            /*
            TLdech2= 0.00011*pow(part.PC,0.25)*pow((1.-1./crys.ratio),2)    //dechanneling length [m]
            part.Dechan = log(1.-random())
            Ldech  = -TLdech2*part.Dechan
            next 2 lines new from sasha - different dechanneling probability
            */
            TLdech2 = 0.01 * part.PC * pow((1. - 1. / crys.ratio), 2); // typical dechanneling length [m]
            Ldech  = 0.005 * TLdech2 * pow((sqrt(0.01 - log(random())) - 0.1), 2); // (dechanneling length)  [m]
            tdech = Ldech / crys.Rcurv;
            Sdech = Ldech * cos(part.xp + 0.5 * tdech);



            if (Ldech < (crys.Cry_length - Lrefl)) {    // for dechanneling part
                //cout<< "Dechanneling! ";
                part.ndechann = part.ndechann + 1;     /////////////COUNT/////////////
                part.effet = 4;
                Dxp = Ldech / crys.Rcurv;
                part.x  = part.x + Ldech * (sin(0.5 * Dxp + part.xp)); //trajectory at channeling exit
                part.y = part.y + Sdech * part.yp;
                part.xp = part.xp + Dxp + (random_gauss() - 0.5) * crys.xpcrit * 0.5;
                Red_S = crys.s_length - Srefl - Sdech;
                //  move in amorphous substance----------------------------
                part.x = part.x + 0.5 * part.xp * Red_S;
                part.y = part.y + 0.5 * part.yp * Red_S;

                if (crys.ZN > 0) {
                    //cout<<"Amorphous! "<<endl;        /////////////COUNT/////////////
                    move_am(crys.IS, crys.NAM, Red_S, crys.DES[crys.IS], crys.DLYI[crys.IS], crys.DLRI[crys.IS], part.xp , part.yp, part.PC); //We call the move_am function
                }

                part.x = part.x + 0.5 * part.xp * Red_S;
                part.y = part.y + 0.5 * part.yp * Red_S;
                //cout<<endl;
            } else {
                //cout<< "Volume Capture! "<<endl;
                //Pas besoin d additionner 1 dans le compte de vol capt car deja fait!!!
                Rlength = crys.Cry_length - Lrefl;
                crys.tchan = Rlength / crys.Rcurv;
                Red_S = Rlength * cos(part.xp + 0.5 * crys.tchan);
                Dxp = (crys.Cry_length - Lrefl) / crys.Rcurv;
                part.x  = part.x + sin(0.5 * Dxp + part.xp) * Rlength; // trajectory at channeling exit
                part.y = part.y + Red_S * part.yp;
                part.xp = part.xp + Dxp +  (random_gauss() - 0.5) * crys.xpcrit * 0.5; //[mrad]
            }
        }
    }
    // move in amorphous substance for big input angles---------------
    else {
        //cout<<"Amorphous! "<<endl;
        part.namor = part.namor + 1;     /////////////COUNT/////////////
        part.effet = 2;
        part.x = part.x + 0.5 * crys.s_length * part.xp;
        part.y = part.y + 0.5 * crys.s_length * part.yp;

        if (crys.ZN > 0) {
            move_am(crys.IS, crys.NAM, crys.s_length, crys.DES[crys.IS], crys.DLYI[crys.IS], crys.DLRI[crys.IS], part.xp , part.yp, part.PC); //We call the move_am function
        }

        part.x = part.x + 0.5 * crys.s_length * part.xp;
        part.y = part.y + 0.5 * crys.s_length * part.yp;
    }

}



////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//THE FUNCTION THAT DO ALL THE SCHEMA USING THE PREVIOUS FUNCTION !!!!!!!!!!111
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////




//called the FORTRAN cyrstal-particle routine from SixTrack
/*
void SimCrys::general(SimCrys& sim, const int& pas, char name, double emitx0, double emity0, double enom, string outputpath, double Mirror, double C_rotation, double C_aperture, double C_offset, double C_tilt, double Crystal_tilt)
{
 
  long int max_npart = 1500000;
  
   char name_coll[50] = name;
   
  int NP = ;
  double ENOM = enom*1e-3;  // [GeV]
  double EMITX0 = emitx0;
  double EMITY0 = emity0;
  double AX = , AY = , BX = , BY = ;
  
  double X_IN[max_npart] ;
  double XP_IN[max_npart];
  double Y_IN[max_npart];
  double YP_IN[max_npart];
  double P_IN[max_npart];
  double S_IN[max_npart] = ;
  
  int flagsec[max_npart] = ;
  bool dowrite_impact = true;
  int name[max_npart] = ;
  
  char C_MATERIAL[6];
  double C_LENGTH = crys.s_length;
  double C_ROTATION = C_rotation;
  double C_APERTURE = C_aperture;
  double C_OFFSET = C_offset;
  double C_TILT[2] = {C_tilt, 0} ;
  
  int LHIT[max_npart] = ;  //position of the particle been hit
  int PART_ABS[max_npart] = ;
  
  double IMPACT[max_npart] = ;
  double INDIV[max_npart] = ;
  double LINT[max_npart] = ;
  
  */
 /* -----------------
!++  Copy particle data to 1-dim array and go back to meters
*/
  /*  
  for(j = 1, napx)
                        rcx(j)  = (xv(1,j)-torbx(ie))/1d3
                        rcxp(j) = (yv(1,j)-torbxp(ie))/1d3
                        rcy(j)  = (xv(2,j)-torby(ie))/1d3
                        rcyp(j) = (yv(2,j)-torbyp(ie))/1d3
                        rcp(j)  = ejv(j)/1d3
                        rcs(j)  = 0d0
                        part_hit_before(j) = part_hit(j)
                        rcx0(j)  = rcx(j)
                        rcxp0(j) = rcxp(j)
                        rcy0(j)  = rcy(j)
                        rcyp0(j) = rcyp(j)
                        rcp0(j)  = rcp(j)
                        ejf0v(j) = ejfv(j)
!
!++  For zero length element track back half collimator length
!
                        if (stracki.eq.0.) then
                                rcx(j) = rcx(j) - 0.5d0*c_length*rcxp(j)
                                rcy(j)= rcy(j) - 0.5d0*c_length*rcyp(j)
                        else
                                Write(*,*) 
     &                          "ERROR: Non-zero length collimator!"
                                STOP
                        endif
                        flukaname(j) = ipart(j)+100*samplenumber
!
                end do
!
!++  Do the collimation tracking
!
                enom_gev = myenom*1d-3
!
!++  Allow



  
   if (crys.IS == 0) {
        C_MATERIAL = "CRY-Si";
    } else if (crys.IS == 1) {
	C_MATERIAL = "CRY-W";
    } else if (crys.IS == 2) {
	C_MATERIAL = "CRY-C";
    } else if (crys.IS == 3) {
	C_MATERIAL = "CRY-Ge";
    }else{
        cerr << "Material of the Crystal not fund !!!" << endl;
    }

      
	
  int count = 0;
  string rest;
    ifstream entree;
    string name2;
    name2 = outputpath + "/ascii_before.csv";
    entree.open(name2.c_str());
    if (entree) {
        while (!entree.eof()) {
            count = count + 1;

            getline(entree, rest, ',');
            X_IN[count] = atof(rest.c_str());
            //cout<<" X    " << part.x<<endl;

            getline(entree, rest, ',');
            Y_IN[count] = atof(rest.c_str());
            //cout<<" Y    " << part.y<<endl;

            getline(entree, rest, ',');
            XP_IN[count] = atof(rest.c_str());
            //cout<<" XP    " << part.xp<<endl;

            getline(entree, rest, ',');
            YP_IN[count] = atof(rest.c_str());
            //cout<<" YP    " << part.yp<<endl;

            getline(entree, rest);
            P_IN[count] = atof(rest.c_str());

	}
    }
    entree.close();
    
	
  //from the file "crystaltrack_dan.f", 08/2014, by Jianfeng Zhang @ LAL
  COLLIMATE_CRY_(name_coll,C_MATERIAL, C_LENGTH,
                        C_ROTATION, 
                        C_APERTURE, C_OFFSET, C_TILT, 
                        X_IN, XP_IN, Y_IN,  
                        YP_IN, P_IN, S_IN, NP, ENOM, LHIT, 
                        PART_ABS, IMPACT, INDIV, LINT,
                        BX,BY,AX,
                        AY,EMITX0,EMITY0,
                        name,flagsec,dowrite_impact);
}

*/

/********************************************************************************
 * 
* original crystal routine
* *******************************************************************************/
void SimCrys::general(SimCrys& sim, const int& pas, long int& nhit, string outputpath, double Mirror, double C_rotation, double C_aperture, double C_offset, double C_tilt, double Crystal_tilt)
{
    srand ( time(NULL) );
    double AV_xp0(0);
    double AV_yp0(0);
    double RMSxp0(0);
    double RMSyp0(0);
    double AV_xp1(0);
    double AV_yp1(0);
    double long RMSxp1(0);
    double long RMSyp1(0);
    double av_pc1(0);
    double av_pc0(0);
    double SIGxp0(0);
    double SIGyp0(0);
    double SIGxp1(0);
    double SIGyp1(0);
    double xp0(0);
    double xp1(0);
    double count(0);
    double zpos(-0.1);
    double xfin(0);
    double yfin(0);
    int const col(15);

    
    //parameter to check whether the particle hit the crystal
    char Proc[50] = " ";  // process of the particle-crystal interaction
    bool crossing = false;  //the particle hit the collimator or not
    int layerflag = 1;  //check the particle exit the layer in "AM"
    //-----------------------------------------------NOW THE PARAMETERS FOR THE CHANGE OF REFERENTIAL ------------------------------------------------
    int mat=0;

    double p0=0.0;
    double zlm=0.0;
    double shift=0.0;
    double x_shift=0.0, xp_shift=0.0, s_shift=0.0; //coordinates after shift/rotation
    double x_rot=0.0, xp_rot=0.0, s_rot=0.0;
    double x_temp=0.0, xp_temp=0.0, s_temp=0.0;  //all the _temp variables are used when you hit the cry from below
    double tilt_int=0.0, x_int=0.0, xp_int=0.0, s_int=0.0;//all the _int variables are used when you hit the cry from below (int=interaction point)
    double x00=0.0;
    double z=0.0;
    double z00=0.0;
    double zp=0.0;
    double dpop=0.0;
    double s=0.0;
    double a_eq=0.0, b_eq=0.0, c_eq=0.0, Delta=0.0;
    int part_abs=-1;
    double x=0.0;
    double xp=0.0;
    double y=0.0;
    double yp=0.0;
    double p=0.0;    //be careful: [Gev]
    double s_in=0.0;

    //        adding variables for the pencil beam. Variables in the absolute reference frame.

    double x_in0=0.0;
    double xp_in0=0.0;
    double y_in0=0.0;
    double yp_in0=0.0;
    double p_in0=0.0;    //be careful: [Gev]
    double s_in0=0.0;
    double s_impact=0.0;
    double totcount(0);

    double mirror=Mirror;
    double tiltangle=0.0;


    double c_length=0.0;      //Length in m
    double c_rotation(C_rotation) ;    //Rotation angle vs vertical in radian initiated and fixed at 0. Can BE CHANGED IF NEED
    double c_aperture(C_aperture) ;   //Aperture in m
    double c_offset(C_offset) ;      //Offset in m
    double c_tilt[2] = {C_tilt, 0} ;  //Tilt in radian
    double c_tilt0[2]  ;    //Tilt in radian
    

    double xp_tangent=0.0;

    double Cry_tilt(Crystal_tilt);     //crystal tilt [rad] ///////////////////////////////////++++++++++++++++++++++++(IF NEEDED, put it into the input file !!! See later)+++++++++++++++////////////////


    
     
     
    //----------------------------------END OF INITATING REFERNENTIAL PARAMETERS-------------------------------------------------------------------------
    c_length = crys.Cry_length;

    if (crys.IS == 0) {
      mat = 8;
    } else if (crys.IS == 1) {
      mat = 9;
    } else if (crys.IS == 2) {
      mat = 10;
    } else if (crys.IS == 3) {
      mat = 11;
    } else {
      cout << "Material of the Crystal not fund !!!" << endl;
    }

    
    if (c_length > 0 ) {
      p0 = part.PC;
      c_tilt0[0] = c_tilt[0];
      c_tilt0[1] = c_tilt[1];
      tiltangle = c_tilt0[0];
    }
    
    
    //new..........................
   // call scatin(p0)   ??????????????

   //  scatin_(p0);
     
      /* EVENT LOOP,  initial distribution is here a flat distribution with
       * xmin=x-, xmax=x+, etc. from the input file
       */
   

      //initialize particle variables
      part.n_part = 0;
      part.nabsorbed = 0;
      part.namor = 0;
      part.nchann = 0;
      part.ndechann = 0;
      part.nout =0;
      part.nvcapt = 0;
      part.nvrefl = 0; 
      part.nhit = 0;
      nhit    = 0;
    double fracab  = 0;
  //  int   n_chan  = 0;          // valentina :initialize to zero the counters for crystal effects 
  //  int  n_VR    = 0;          
   // int    n_amorphous = 0;   
  
  
    double impact = 0;
    double indiv = 0;
    double lint = 0;
  
  
    //  c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    //c        WRITE(*,*) ITURN,ICOLL  
 /*   do j = 1, nev 
	 //c 
	 impact(j) = -1.0;
    lint(j)   = -1.0;
    indiv(j)  = -1.0;
    

    idx_proc = name(j)-100*samplenumber;

    if (ITURN == 1){                       //                   !daniele
      bool_proc_old(idx_proc)=-1;} //                        !daniele
    //  elseif (bool_proc(idx_proc) != -1){
    else{
      bool_proc_old(idx_proc)=bool_proc(idx_proc); //                  !daniele
    } //                                          !daniele
    
    PROC= "out"; // !the default process is 'out'                 !daniele
    bool_proc(idx_proc)=-1;         //     !daniele
*/

  int  nabs = 0;
        
    //ofstream sortie("results.csv");
    //sortie << "xp1" <<","<<  "xp0"<<","<<  "x"<<","<<  "y"<<endl;
    ofstream sortieIco;
    string name1;
    
    //file with the particle coordinate after the crystal
    name1 = outputpath + "/ascii_after.csv";

    sortieIco.open(name1.c_str());

    string rest;

    ifstream entree;
    string name2;
    name2 = outputpath + "/ascii_before.csv";

    entree.open(name2.c_str());
    if (entree) {
      while (!entree.eof()) {
	count = count + 1;
	
	//particle hit the crystal or not
	bool hitflag = false;
	  
	
	getline(entree, rest, ',');
	part.x = atof(rest.c_str());
	//cout<<" X    " << part.x<<endl;

	getline(entree, rest, ',');
	part.y = atof(rest.c_str());
	//cout<<" Y    " << part.y<<endl;

	getline(entree, rest, ',');
	part.xp = atof(rest.c_str());    
	//cout<<" XP    " << part.xp<<endl;

	getline(entree, rest, ',');
	part.yp = atof(rest.c_str());
	//cout<<" YP    " << part.yp<<endl;

	getline(entree, rest);
	part.PC = atof(rest.c_str());
	
	if (part.PC != 0) { // POUR OTER LA DERNIERE PART QUI VAUT 0 POUR TOUT

	  av_pc0 = av_pc0 + part.PC;
	  
	  //part.x=random_gauss()*0.00115;
	  //part.xp=random_gauss()*sqrt(5*1e-6);
	  //cout<<"X et Y avt :"<<part.x<<" , "<<part.y<<endl;

	  // IL FAUT PASSER DES CM EN M!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
	  //bdelta=tan(part.xp)*10;
	  part.x = (part.x); //-0.0045;
	  
	  
	  //bdelta=tan(part.yp)*10;
	  part.y = (part.y); //-0.0045;

	  //cout<<"X et Y apres :"<<part.x<<" , "<<part.y<<endl;
	  xp0 = part.xp;

	  //part.y=random_gauss()*sqrt(0.00095);
	  //part.yp=5*1e-6;

	  //part.PC=120.976;
	  //cout<<"  X  :"<<part.x<<"  XP  :"<<part.xp<<"  Y  :"<<part.y<<"  YP  :"<<part.yp<<"  PC  :"<<part.PC<<endl;
	  AV_xp0 = AV_xp0 + part.xp;               //average x angle before impact
	  AV_yp0 = AV_yp0 + part.yp;               //average y angle before impact
	  RMSxp0 = RMSxp0 + part.xp * part.xp;        //rms x angle before impact
	  RMSyp0 = RMSyp0 + part.yp * part.yp;   //rms y angle before impact

	  // NOW WE NEED TO CHANGE THE REFERENTIAL!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	  //change the beam frame to the crystal frame !!!!!!!!!
	  // since the interaction between the particle and the crystal is observed in the 
	  // crystal reference system

	  //initial coordinates in the beam frame
	  s   = 0;
	  x   = part.x;
	  xp  = part.xp;
	  z   = part.y;
	  zp  = part.yp;
	  p   = part.PC;

	  x_temp = 0;
	  x_int = 0;
	  x_rot = 0;
	  x_shift = 0;
	  s_temp = 0;
	  s_int = 0;
	  s_rot = 0;
	  s_shift = 0;
	  xp_int = 0;
	  xp_temp = 0;
	  xp_rot = 0;
	  xp_shift = 0;
	  shift = 0;
	  tilt_int = 0;
	  dpop = (p - p0) / p0;  //momentum spread
                
        
	  /*---------------DANIELE-------------------------
	    ! corrected position for variable association for FirstImpact.dat
	    ! also for a vertical crystal case.
	    ! Only x_in0 (i.e. b) have to be assigned after the change of reference frame
	    !-----------------------------------------------*/

	  s_in0   = s_in;                    //daniele
	  //        x_in0(j)   = x                         //!daniele
	  xp_in0  = xp;                         //daniele
	  y_in0   = z;                         //daniele
	  yp_in0  = zp;                         //daniele
	  p_in0   = p;                         //daniele
        
	  /*----------------DANIELE------------------*/

                
	  //**********************************************************************************                
	  
	  // ------------------------------transform particle coordinates to get into collimator coordinate system -----------------------------------------
	  

	  // first check whether particle was lost before
	  totcount = totcount + 1;
	  if ((x < 99.0 * 1e-3) && (z < 99.0 * 1e-3)) {             ////////////////////////// IF NUMERO 1//////////////////
	    
	    // first do rotation into collimator frame
	    
	    x  = part.x * cos(c_rotation) + sin(c_rotation) * part.y;
	    z  = part.y * cos(c_rotation) - sin(c_rotation) * part.x;
	    xp = part.xp * cos(c_rotation) + sin(c_rotation) * part.yp;
	    zp = part.yp * cos(c_rotation) - sin(c_rotation) * part.xp;
	    
	    //cout<<"PARTICULE NON PERDU>>><<"<<endl;
	    
	    //To define the correct value of the mirror for the install locations
	    // of the crystal in SPS & LHC, inner side or outside of the vacuum chamber????
	    //  by Zhang @ CERN, 23/04/2014.....
	    //mirror = 1;
		    
	    
	    x  = mirror * x;
	    xp = mirror * xp;

	    //Shift with opening and offset

	    x  = x - c_aperture / 2 - c_offset;

	    // Include collimator tilt

	    if (tiltangle > 0) {                ////////////////////////// IF NUMERO 2//////////////////
	      xp = xp - tiltangle;
	    }                               //////////////////////////FIN IF NUMERO 2//////////////////
	    else if (tiltangle < 0) {                 ////////////////////////// IF NUMERO 3//////////////////
	      x  = x + sin(tiltangle) * c_length;
	      xp = xp - tiltangle;
	    }                 //////////////////////////FIN IF NUMERO 3//////////////////

	    if (Cry_tilt < 0) {                ////////////////////////// IF NUMERO 4//////////////////
	      
	      s_shift = s;
	      shift = crys.Rcurv * (1 - cos(Cry_tilt));
	      if (Cry_tilt < (-crys.cry_bend) ) {                ////////////////////////// IF NUMERO 5//////////////////
		shift = ( crys.Rcurv * ( cos ( - Cry_tilt) - cos( crys.cry_bend - Cry_tilt) ) );
	      }                 //////////////////////////FIN IF NUMERO 5//////////////////
	      x_shift = x - shift;
	    }                                  //////////////////////////fin  IF NUMERO 4//////////////////
	    else {                //////////////////////////IF NUMERO 6//////////////////
	      //cout<< "CRY-TILT POSSITIVE ++ or nul"<<endl;
	      s_shift = s;
	      x_shift = x;
	    }                 //////////////////////////fin IF NUMERO 6//////////////////
	    // cout<<"debug - S shift" <<  s_shift<<endl;
	    //cout<<"debug - X shift" <<  x_shift <<endl;
	    //
	    //    2nd transformation: rotation
	    s_rot = x_shift * sin(Cry_tilt) + s_shift * cos(Cry_tilt);
	    x_rot = x_shift * cos(Cry_tilt) - s_shift * sin(Cry_tilt);
	    xp_rot = xp - Cry_tilt;
	    //cout<< "debug - S rot" <<  s_rot <<endl;
	    //cout<<"debug - X rot" <<  x_rot <<endl;
	    //cout<<"debug - XP rot" <<  xp_rot<<endl;
	    
	    //    3rd transformation: drift to the new coordinate s=0
	    xp = xp_rot;
	    x = x_rot - xp_rot * s_rot;
	    z = z - zp * s_rot;
	    s = 0;
	    //cout<<"debug - S cryRF" <<  s_rot <<endl;
	    //cout<<"debug - X cryRF" <<  x_rot <<endl;
	    //cout<<"debug - XP cryRF" <<  xp_rot  <<endl;


	    
	    //********************************************************************************************************************
	    //----------------------------------NOW CHECK IF THE PARTICLE HIT the crystal---------------------------------------------------------------
	    //cout<<" X    " << x<<endl;
	    //Particle hit the crystal
	    if (x >= 0) {                ////////////////////////// IF NUMERO 7//////////////////
	      s_impact = s_in0;       //!(for the first impact)
              
	      if(0){
		cout<<"hit the cry entrance face"<<endl;
		cout<<"impact at s = "<< s_impact<< ",  x = "<<x_in0<<endl;
		cout<<"with angle xp = "<<xp<<endl;
		cout<<"s before: "<< s<<endl;
	      }
	      
	      part.x = x;
	      part.xp = xp;
	      part.y = z;
	      part.yp = zp;
	      part.PC = p;

	      //------------------------------------------------------------------------------
	      // The proton crystal routine in Fortran version 
	      //
	      //------------------------------------------------------------------------------
	      
	      //Fortran crystal routine in SixTrack
	      // Modified by Jianfeng Zhang @ LAL, 30/04/2014
	      double crysCry_length=crys.Cry_length;
	      double crysRcurv=crys.Rcurv;
	      double crysC_xmax=crys.C_xmax;
	      double crysC_ymax=crys.C_ymax;
	      double crysAlayer=crys.Alayer;
	      int crysC_orient=crys.C_orient;   
	      double crysmiscut = crys.miscut;
	      int IS =0;
	      double length =0.0;
	      IS = mat -7 ;
	      length = crys.Cry_length;
	      double eUmIS = crys.eUm[crys.IS];
	      double AIIS = crys.AI[crys.IS];
	
	   
	      
	      double crysparaDes[4];
	      crysparaDes[0]=crys.DES[0],crysparaDes[1]=crys.DES[1],crysparaDes[2]=crys.DES[2],crysparaDes[3]=crys.DES[3];
	      double crysparaDlyi[4];
	      crysparaDlyi[0]=crys.DLYI[0],crysparaDlyi[1]=crys.DLYI[1],crysparaDlyi[2]=crys.DLYI[2],crysparaDlyi[3]=crys.DLYI[3];
	      double crysparaDlri[4];
	      crysparaDlri[0]=crys.DLRI[0],crysparaDlri[1]=crys.DLRI[1],crysparaDlri[2]=crys.DLRI[2],crysparaDlri[3]=crys.DLRI[3];
	      
	      double crysparaeUmIS[4];
	      crysparaeUmIS[0]=crys.eUm[0], crysparaeUmIS[1]=crys.eUm[1], crysparaeUmIS[2]=crys.eUm[2], crysparaeUmIS[3]=crys.eUm[3];
	      double crysparaAIIS[4];
	      crysparaAIIS[0]=crys.AI[0],crysparaAIIS[1]=crys.AI[1],crysparaAIIS[2]=crys.AI[2],crysparaAIIS[3]=crys.AI[3];
	            
	      //	cout << "before enter crystal routine...."<<endl;
			
//	      cryst_(&IS,&x,&xp,&z,&zp,&p,&length,&crysCry_length,&crysRcurv,&crysC_xmax,&crysC_ymax,&crysAlayer,&crysC_orient,&crysmiscut, &Proc,&layerflag, &part.Chann);
	      
//	     cryst_(&length,&layerflag, crypara,partpara,crysparaDes,crysparaDlri,crysparaDlyi,crysparaeUmIS,crysparaAIIS,Proc);
	      
	      cryst_(&length,&layerflag,Proc,&crys.L_chan,&crys.tchan,&crys.Alayer,&crys.ymin,&crys.ymax,&crys.Rcurv,
		     &crys.s_length,&crys.IS,&crys.NAM, 
                    &crys.xpcrit0,&crys.xpcrit,&crys.Rcrit,&crys.ratio,&crys.C_orient,&crys.miscut,&crys.Cry_length, 
                    &crys.ZN,&crys.C_xmax,&crys.C_ymax,
                    &part.x,&part.xp,&part.y,&part.yp,&part.PC,&part.s,&part.Chann,&part.xp_rel,&part.Ang_rms,
    	            & part.Ang_avr,&part.Vcapt,&part.Dechan,crysparaDes,crysparaDlri,crysparaDlyi,crysparaeUmIS,crysparaAIIS);
	      
			 
	     
	      //convert back the Fortran variables to C++ variables
	      //mat = IS +7; 
	      
	 //     cout << "IS  in Fortran = ..." << crys.IS << endl;

	      crys.IS = crys.IS - 1; 
	      
	   //   	      cout << "IS in C++ = ..." << crys.IS << endl;
		      
	      crys.Cry_length = length;
	      
	          x=part.x;
		  xp=part.xp;
		  z=part.y;
		  zp=part.yp;
		  p=part.PC;
	          s=part.s;
	            
		
	      /*
	      crys.Rcurv = crysRcurv;
	      crys.C_xmax = crysC_xmax;
	      crys.C_ymax = crysC_ymax;
	      crys.Alayer = crysAlayer;
	      crys.C_orient = crysC_orient;
	      crys.miscut = crysmiscut;
	      */
	      
	      //              cout<< "ICI AHHHHHHHHHHHHHHHHHHHHHHHHHHHH"<<endl;
	      //                         //cout<<" xp rel and xp crit " << part.xp_rel<< " " << crys.xpcrit<<endl;
	      //                         bases(sim);
	      //                         //cout << part.xp<< " XP "<<endl;
	      //                         bool crossing(cross(sim));
	      //                         //cout<<part.xp-xp0<<endl;
	      //                         //cout<<" xp rel and xp crit " << part.xp_rel<< " " << crys.xpcrit<<endl;
	      //                         bool layering;
	      //                         if (crossing != 0) {             ////////////////////////// IF NUMERO 8//////////////////
	      //                             layering = layer(sim);
	      //                             //cout<<"Crossing !!!!"<<endl;
	      // 
	      //                             if (layering == true) {             ////////////////////////// IF NUMERO 9//////////////////
	      //                                 parameters(sim);
	      //                                 //cout<<" Layering =true " <<endl;
	      // 
	      //                                 if ((abs(part.xp_rel) < crys.xpcrit)) {               ////////////////////////// IF NUMERO 10//////////////////
	      //                                     channel(sim);
	      //                                     //cout<< " CHANNELING " << endl;
	      //                                 }                 ////////////////////////// fin IF NUMERO 10//////////////////
	      //                                 else {                ////////////////////////// IF NUMERO 11//////////////////
	      //                                     reflection(sim);
	      //                                     //cout<<" REFLECTION "<<endl;
	      //                                 }                 ////////////////////////// fin IF NUMERO 11//////////////////
	      // 
	      // 
	      //                             }                 ////////////////////////// fin IF NUMERO 9//////////////////
	      //                         }                 ////////////////////////// fin IF NUMERO 8//////////////////
	      
	      //--------------------------------------------------------------------------------------------------
	      
                        
	      // xp = part.xp;
	      
	      //}

	      s = crys.Rcurv * sin(crys.cry_bend);
	      zlm = crys.Rcurv * sin(crys.cry_bend);
	      //cout<<"process:"<<PROC<<endl;
	      //cout<<"s after"<< s<<endl;
	      //cout<<"part.xp"<<part.xp<<endl;
	      
	      //cout<<"----------------------1111111111111111111111111-----------------"<<endl;
	      
	   if(layerflag == 0){
	      part.namor = part.namor + 1;     /////////////COUNT/////////////
              part.effet = 2;
	      }
	      
	      if (strncmp(Proc,"out",3)==0) 
		crossing  = false;
	      else
		crossing = true;
	      
	      if (crossing == true){
 		hitflag = true;
 		nhit = nhit + 1;
 		part.nhit = nhit; 
	//	lhit = 100000000*ie + ITURN;
		impact = x_in0;            
		indiv = xp_in0;    
	      }	
	    }
	    // Particle not hit the crystal
	    ////////////////////////// fin IF NUMERO 7//////////////////
	    else {                //////////////////////////  IF NUMERO 12//////////////////
	      
	      //cout<<"OU LA ABBBBBBBBBBBB" <<endl;
	      xp_tangent = sqrt((-2 * x * crys.Rcurv + x * x) / (crys.Rcurv * crys.Rcurv));
	      //cout<<"RCURV      "<<crys.Rcurv<<endl;
	      //cout<<"tangent"<<xp_tangent<<"angle"<<xp<<endl;
	      
	      //cout<<"s tan " << crys.Rcurv*sin(xp_tangent)<<endl;
	      //cout<< " s tot "<< c_length<< crys.Rcurv*sin(crys.cry_bend)<<endl;
	      if ( xp >= xp_tangent) {                //////////////////////////  IF NUMERO 13//////////////////
		
		//if it hits the crystal, calculate in which point and apply the
		//transformation and drift to that point
		a_eq = (1 + xp * xp);
		b_eq = 2 * xp * (x - crys.Rcurv);
		c_eq = -2 * x * crys.Rcurv + x * x;
		Delta = b_eq * b_eq - 4 * (a_eq * c_eq);
		s_int = (-b_eq - sqrt(Delta)) / (2 * a_eq);
		//cout<< " s int "<<s_int <<endl;
		if (s_int < crys.Rcurv * sin(crys.cry_bend)) {               ////////////////////////// fin IF NUMERO 14//////////////////
		  //  transform to a new ref system:shift and rotate
		  
		  x_int = xp * s_int + x;
		  xp_int = xp;
		  z = z + zp * s_int;
		  x = 0;
		  s = 0;
		  

		  //               tilt_int=2*X_int/S_int
		  tilt_int = s_int / crys.Rcurv;
		  xp = xp - tilt_int;
		  
		  //---------------------------------------------------------------------------------------
		  //  Fortran Cyrstal routine for the protons 
		  //
		  //---------------------------------------------------------------------------------------
		  
		  part.x = x;
		  part.xp = xp;
		  part.y = z;
		  part.yp = zp;
		  part.PC = p;
		  crys.Cry_length = crys.Cry_length - (tilt_int * crys.Rcurv);
		  
				
				
		  /*				
				
		  //Fortran crystal routine in SixTrack
		  // Modified by Jianfeng Zhang @ LAL, 30/04/2014
		  CRYST_(mat-7,x,xp,z,zp,p,cry_length);
		  */				
		  
				 
		  //------------------------------------------------------------------------------
		  // The proton crystal routine in Fortran version 
		  //
		  //------------------------------------------------------------------------------
			
		  //Fortran crystal routine in SixTrack
		  // Modified by Jianfeng Zhang @ LAL, 30/04/2014
		  double crysCry_length=crys.Cry_length;
		  double crysRcurv=crys.Rcurv;
		  double crysC_xmax=crys.C_xmax;
		  double crysC_ymax=crys.C_ymax;
		  double crysAlayer=crys.Alayer;
		  int crysC_orient=crys.C_orient;   
		  double crysmiscut = crys.miscut;
		  int IS =0;
		  double length =0.0;
		  IS = mat -7 ;
		  length = crys.Cry_length;
		  double eUmIS = crys.eUm[crys.IS];
		  double AIIS = crys.AI[crys.IS];
		  
	      
	      double crysparaDes[4];
	      crysparaDes[0]=crys.DES[0],crysparaDes[1]=crys.DES[1],crysparaDes[2]=crys.DES[2],crysparaDes[3]=crys.DES[3];
	      double crysparaDlyi[4];
	      crysparaDlyi[0]=crys.DLYI[0],crysparaDlyi[1]=crys.DLYI[1],crysparaDlyi[2]=crys.DLYI[2],crysparaDlyi[3]=crys.DLYI[3];
	      double crysparaDlri[4];
	      crysparaDlri[0]=crys.DLRI[0],crysparaDlri[1]=crys.DLRI[1],crysparaDlri[2]=crys.DLRI[2],crysparaDlri[3]=crys.DLRI[3];
	      
	      double crysparaeUmIS[4];
	      crysparaeUmIS[0]=crys.eUm[0], crysparaeUmIS[1]=crys.eUm[1], crysparaeUmIS[2]=crys.eUm[2], crysparaeUmIS[3]=crys.eUm[3];
	      double crysparaAIIS[4];
	      crysparaAIIS[0]=crys.AI[0],crysparaAIIS[1]=crys.AI[1],crysparaAIIS[2]=crys.AI[2],crysparaAIIS[3]=crys.AI[3];
	      
	
	      
	      
	 
	      //	cout << "before enter crystal routine...."<<endl;
				
	      cryst_(&length,&layerflag,Proc,&crys.L_chan,&crys.tchan,&crys.Alayer,&crys.ymin,&crys.ymax,&crys.Rcurv,
		     &crys.s_length,&crys.IS,&crys.NAM, 
                    &crys.xpcrit0,&crys.xpcrit,&crys.Rcrit,&crys.ratio,&crys.C_orient,&crys.miscut,&crys.Cry_length, 
                    &crys.ZN,&crys.C_xmax,&crys.C_ymax,
                    &part.x,&part.xp,&part.y,&part.yp,&part.PC,&part.s,&part.Chann,&part.xp_rel,&part.Ang_rms,
    	            & part.Ang_avr,&part.Vcapt,&part.Dechan,crysparaDes,crysparaDlri,crysparaDlyi,crysparaeUmIS,crysparaAIIS);
	       
			  //convert back the Fortran variables to C++ variables
	      //mat = IS +7; 
	      crys.Cry_length = length;
	      crys.IS = crys.IS - 1;
	      

	          
	          x=part.x;
		  xp=part.xp;
		  z=part.y;
		  zp=part.yp;
		  p=part.PC;
	           s=part.s;
	      
			/*
		  crys.Rcurv = crysRcurv;
		  crys.C_xmax = crysC_xmax;
		  crys.C_ymax = crysC_ymax;
		  crys.Alayer = crysAlayer;
		  crys.C_orient = crysC_orient;
		  crys.miscut = crysmiscut;
		  */

		  /*
		    bases(sim);
		    
		    //cout<<"part.xp"<<part.xp<<endl;

		    bool crossing(cross(sim));
		    //cout<<"Crossing !!!!"<<endl;
		    bool layering(layer(sim));

		    if (layering == true) {
		    //cout<<" Layering =true " <<endl;
		    parameters(sim);
		    if ((abs(part.xp_rel) < crys.xpcrit)) {
		    channel(sim);
		    //cout<< " CHANNELING " << endl;
		    } else {
		    reflection(sim);
		    //cout<< " REFLECTION " << endl;
		    }
		    
		    }
		  */
		  //---------------------------------------------------------------------------------------


		  
                                //cout<<"part.xp"<<part.xp<<endl;
                                //cout<<"----------------------2222222222222222222-----------------"<<endl;
		  
		  s = crys.Rcurv * sin(crys.cry_bend - tilt_int);
		  zlm = crys.Rcurv * sin(crys.cry_bend - tilt_int);
				

	      if(layerflag == 0){
	      part.namor = part.namor + 1;     /////////////COUNT/////////////
              part.effet = 2;
	      }
		  if (strncmp(Proc,"out",3)==0)  
		    crossing  = false;
		  else
		    crossing = true;
		  
		  if (crossing == true) {              //////////////////////////  IF NUMERO 15//////////////////
		    x_rot = x_int;
		    s_rot = s_int;
		    xp_rot = xp_int;
		    s_shift = s_rot * cos(-Cry_tilt) + x_rot * sin(-Cry_tilt);
		    x_shift = -s_rot * sin(-Cry_tilt) + x_rot * cos(-Cry_tilt);
		    xp_shift = xp_rot + Cry_tilt;
		    if (Cry_tilt < 0) {                //////////////////////////  IF NUMERO 16//////////////////
		      s_impact = s_shift;
		      x_in0 = x_shift + shift;
		      xp_in0 = xp_shift;
		    }                 ////////////////////////// fin IF NUMERO 16//////////////////
		    else {                //////////////////////////  IF NUMERO 17//////////////////
		      x_in0 = x_shift;
		      s_impact = s_shift;
		      xp_in0 = xp_shift;
		    }                 ////////////////////////// fin IF NUMERO 17//////////////////
		    hitflag = true;
		    nhit = nhit + 1;
		    part.nhit = nhit; 
		    //lhit = 100000000*ie + ITURN;
		    impact = x_in0;            
		    indiv = xp_in0;  
		  }                 ////////////////////////// fin IF NUMERO 15//////////////////
		  
  

		  x_temp = x;
		  s_temp = s;
		  xp_temp = xp;
		  s = s_temp * cos(-tilt_int) + x_temp * sin(-tilt_int);
		  x = -s_temp * sin(-tilt_int) + x_temp * cos(-tilt_int);
		  xp = xp_temp + tilt_int;
		  
		  //     2nd: shift back the 2 axis
		  x = x + x_int;
		  s = s + s_int;
		  
		}                 ////////////////////////// fin IF NUMERO 14//////////////////
		// Finish if (s_int < crys.Rcurv * sin(crys.cry_bend))
		else { ///////////////////////////////////////////////PROBLEM : TOO MUCH OF PARTICULES WITHOUT ANY CHANGE OF ANGLE///////////THE PROB IS HERE WITH THOSES TWO "ELSE"/////////////////////////
		  
		  s = crys.Rcurv * sin(crys.Cry_length / crys.Rcurv);
		  x = x + s * xp;
		  z = z + s * zp;
		  ++part.nout;
		  //cout<< "the particule does not hit the crystal (Trop basse dans neg.)"<<endl;
		}// Finish else if (s_int < crys.Rcurv * sin(crys.cry_bend))
                
	      }else {
		s = crys.Rcurv * sin(crys.Cry_length / crys.Rcurv);
		x = x + s * xp;
		z = z + s * zp;
		++part.nout;
		//cout<< "the particule does not hit the crystal (negatif et plus petit angle qu xp_tang)"<<endl;
	      }
	      //}                 ////////////////////////// fin IF NUMERO 13//////////////////
	      // Finish else  if ( xp >= xp_tangent)
              
	    }////////////////////// 12
	    // Finish particle not hit the crystal  else (x < 0)
	    
       
       
	    
	    //**************************************************************************************************************************
	    //trasform back from the crystal to the collimator reference system
	    //   1st: un-rotate the coordinates
	    

	    x_rot = x;
	    s_rot = s;
	    xp_rot = xp;
	    //cout<< "debug - S Cry RF 2" ,  s_rot <<endl;
	    //cout<<"debug - X Cry RF 2" ,  x_rot <<endl;
	    //cout<<"debug - XP Cry RF 2" ,  xp_rot <<endl;
	    s_shift = s_rot * cos(-Cry_tilt) + x_rot * sin(-Cry_tilt);
	    x_shift = -s_rot * sin(-Cry_tilt) + x_rot * cos(-Cry_tilt);
	    xp_shift = xp_rot + Cry_tilt;
	    //     2nd: shift back the reference frame
	    if (Cry_tilt < 0) {                //////////////////////////  IF NUMERO 18//////////////////
	      s = s_shift;
	      x = x_shift + shift;
	      xp = xp_shift;
	    }                 ////////////////////////// fin IF NUMERO 18//////////////////
	    else {
	      x = x_shift;
	      s = s_shift;
	      xp = xp_shift;
	    }
	    //     3rd: shift to new S=Length position
	    x = xp * (c_length - s) + x;
	    z = zp * (c_length - s) + z;
	    s = c_length;


	    nabs = 0;   //particle lost or not lost.....
	    part.effet = nabs;
	    //cout<< "debug - S Coll RF 2" ,  s_rot <<endl;
	    //cout<<"debug - X Coll RF 2" ,  x_rot <<endl;
	    //cout<<"debug - XP Coll RF 2" ,  xp_rot <<endl;

	    //cout<<"X1_coll"<<x<<"Z1_coll"<<z<<"XP1_coll"<<xp<<"ZP1_coll"<<zp <<"s" <<s<<endl;

	
	    if(strncmp(Proc,"out",3)==0){
	    part.nout = part.nout + 1;     /////////////COUNT/////////////
		part.effet = 1; 
	    }
	      if(layerflag == 0){
		part.namor = part.namor + 1;     /////////////COUNT/////////////
		part.effet = 2;
	      }
	      if(strncmp(Proc,"AM",2)==0){
		part.namor = part.namor + 1;     /////////////COUNT/////////////
        part.effet = 2;
	      }
	      if(strncmp(Proc,"CH",2)==0){
		part.nchann = part.nchann + 1;
		part.effet = 3; 
	      }
	   if(strncmp(Proc,"DC",2)==0){
		part.ndechann = part.ndechann + 1;     /////////////COUNT/////////////
		part.effet = 4;;
	      }
	      if(strncmp(Proc,"VR",2)==0){
	       part.nvrefl = part.nvrefl + 1;     /////////////COUNT/////////////
        part.effet = 5;
	      }
	      if(strncmp(Proc,"VC",2)==0){
	     part.nvcapt = part.nvcapt + 1;     /////////////COUNT/////////////
            part.effet = 6;}
             if(strncmp(Proc,"absorbed",8==0)){
 	     part.nabsorbed = part.nabsorbed +1;
 	      part.effet = 7;
 	    }
	     /*           
 if (PROC(1:3).eq.'pne')
                 PROC_dan=7
             if (PROC(1:3).eq.'ppe')
                 PROC_dan=8
               if (PROC(1:4).eq.'diff')
                 PROC_dan=9
               if (PROC(1:4).eq.'ruth')
                 PROC_dan=10
	     */
    
            
              //  ?????????????????????????????????????
	      if(part.effet == 2) {
		part_abs = 0;
	      }else {
		part_abs = 1;
	      }


	      //=========================================================================================================================
	      //  Transform back to particle coordinates with opening and offset

	      if (part_abs == 0) {               //////////////////////////  IF NUMERO 19//////////////////
		//
		//  Include collimator tilt

		if (tiltangle > 0) {
		  x = x  + tiltangle * c_length;
		  xp = xp + tiltangle;
		} else if (tiltangle < 0) {
		  x = x + tiltangle * c_length;
		  xp = xp + tiltangle;
		  
		  x = x - sin(tiltangle) * c_length;
		}

		//  Transform back to particle coordinates with opening and offset
		
		z00 = z;
		x00 = x + mirror * c_offset;
		x = x + c_aperture / 2 + mirror * c_offset;

		//+  Now mirror at the horizontal axis for negative X offset

		x = mirror * x;
		xp = mirror * xp;

		//  Last do rotation into collimator frame
		
		part.x  = x  * cos((-1) * c_rotation) + z * sin((-1) * c_rotation);
		part.y  = z  * cos((-1) * c_rotation) - x  * sin((-1) * c_rotation);
		part.xp = xp * cos((-1) * c_rotation) + zp * sin((-1) * c_rotation);
		part.yp = zp * cos((-1) * c_rotation) - xp * sin((-1) * c_rotation);

		
		
		//**********************************************************************************
		//For output.... by Zhang @ CERN, 24/04/2014
		//p_in = (1 + dpop) * p0;
		//From Dianiele
		//p_in = p;
		//s_in = s_in + s;
               part.PC = p;
	       part.s = part.s + s;
		
		if (part.effet == 1) {              //////////////////////////  IF NUMERO 21//////////////////

		  part.x = 99.99e-3;
		  part.y = 99.99e-3;
		  cout<< "Mean that the particules is lost"<<endl;
		}              ////////////////////////// fin IF NUMERO 21//////////////////
	      }                  ////////////////////////// fin IF NUMERO 19//////////////////


	      //  End of check for particles not being lost before   (see @330)

	      //}                  ////////////////////////// fin IF NUMERO 12//////////////////
	      //  End of loop over all particles

	      
	      
	  }  //collimator with length = 0                  ////////////////////////// fin IF NUMERO 1//////////////////

	  // =================================================================================FIN DE RETOUR AU COORDONEE DE ICOSIM====================================================================
                
		  part.xp=xp;
		  part.yp=zp;
		  part.y=z;
		  part.x=x;
		  part.PC=p;
                
	  xp1 = part.xp;

	  
	  //if(xp1-xp0!=0){
	  //sortie << xp1 << "," << xp0<<"," << xfin<<"," << yfin<<endl;
	  //}
	  //cout<<"ICIIIIIIIIIIIIIIIIIIIIIIIII"<<part.PC<<endl;

	  
	  AV_xp1 = AV_xp1 + part.xp;               //average x angle after interaction
	  AV_yp1 = AV_yp1 + part.yp;               //average y angle after interaction
	  RMSxp1 = RMSxp1 + part.xp * part.xp;      //rms x angle after interaction
	  RMSyp1 = RMSyp1 + part.xp * part.xp;      //rms y angle after interaction
	  av_pc1 = av_pc1 + part.PC;

      if(0){
	  cout<<"  X apres :"<<part.x<<"  XP apres :"<<part.xp<<"  Y apres :"<<part.y<<"  YP apres :"<<part.yp<<"  PC apres :"<<part.PC<<endl;
	  //cout<<endl;
      }
      //only record the particles that hit the collimator.
      if(hitflag != true){
	  //sortieIco.setf(ios::scientific);
	  sortieIco << part.x << "," << part.y << "," << part.xp << "," << part.yp << "," << part.PC  << endl;
	  //}//end of for....
      }
      
	}
      } //End of while


    } else {
      cerr<< "We can not open the file !" << endl;
    }
    entree.close();



    AV_xp0 = AV_xp0 / count;
    AV_yp0 = AV_yp0 / count;
    
    RMSxp0 = sqrt(RMSxp0 / (count));
    RMSyp0 = sqrt(RMSyp0 / (count));
    
    SIGxp0 = sqrt((RMSxp0 * RMSxp0) - (AV_xp0 * AV_xp0));
    SIGyp0 = sqrt((RMSyp0 * RMSyp0) - (AV_yp0 * AV_yp0));

    AV_xp1 = AV_xp1 / count;
    AV_yp1 = AV_yp1 / count;
    RMSxp1 = sqrt(RMSxp1 / (count));
    RMSyp1 = sqrt(RMSyp1 / (count));
    
    SIGxp1 = sqrt((RMSxp1 * RMSxp1) - (AV_xp1 * AV_xp1));
    SIGyp1 = sqrt((RMSyp1 * RMSyp1) - (AV_yp1 * AV_yp1));
    av_pc0 = av_pc0 / count;
    av_pc1 = av_pc1 / count;
    //--------------------------------------------------------------------------
    //-------------write a summary file-----------------------------------------
    //cout<<endl;
    //cout<<endl;
    //cout<<endl;
    
    //cout<<"avg xp_in: "<< AV_xp0 <<"//  avg xp_out: "<< AV_xp1<<endl;
    //cout<<"rms xp_in: "<< RMSxp0 <<"//  rms xp_out: "<< RMSxp1<<endl;
    //cout<<"sigma xp_in: "<< SIGxp0 <<"//  sigma xp_out: "<< SIGxp1<<endl;
    //cout<<"avg pc_in: "<< av_pc0 <<"//  avg pc_out: "<< av_pc1<<endl;
    


    sortieIco.close();
    //     cout<<endl;
    //     affiche();
    file_out(pas, outputpath);
    
    if(1){
      cout << "number of particles before the crystal: " << count<< endl;
      cout << "number of particle hit the crystal: " << part.nhit << endl;
      cout << "number of particles after the crystal: " << count - part.nhit << endl;
    }
    //cout<<endl;
    //cout<<endl;
    //cout<<endl;
}



///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

//WE WRITE SOME FUNCTIONS IMPORTANT FUNCTION NOW

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


void SimCrys::move_am(int IS, int NAM, double DZ, double DEI, double DLY, double DLR, double& xp, double& yp, double& PC)
{
    //cout<<"The move_am function is use !!!!! "<<endl;
    double DYA(0);
    double Wp(0);
    if (NAM == 0) {
        return;
    }
    xp = xp * 1000;
    yp = yp * 1000;
    //cout << "xp avt: "<<xp<<" yp avt: "<<yp<<endl;
    //DEI ---> dE/dx (stoping energy)
    PC = PC - DEI * DZ;    //energy lost beacause of ionization process[GeV]

    // Coulomb Scattering
    //DYA ---> rms of coloumb scattering
    DYA = (14 / PC) * sqrt(DZ / DLR);    //rms scattering (mrad)
    xp = xp + DYA * random_gauss_cut();
    yp = yp + DYA * random_gauss_cut();

    //Elastic scattering
    //cout<<"Plusieurs possibilitees : "<<endl;
    if (random() <= (DZ / crys.DLAI[IS])) {
        //cout<<"Poss. 1 !!!! "<<endl;
        xp    = xp + 196 * random_gauss_cut() / PC;        //angle elast. scattering in R plane
        yp    = yp + 196 * random_gauss_cut() / PC;
    }

    //Diffraction interaction
    if (random() <= (DZ / (DLY * 6.143))) {
        //cout<<"Poss. 2 !!!! "<<endl;
        xp = xp + 257 * random_gauss_cut() / PC;          // angle elast. scattering in R plane[mr]
        yp = yp + 257 * random_gauss_cut() / PC;
        Wp = random();
        PC = PC - 0.5 * pow(0.3 * PC, Wp);         // m0*c = 1 GeV/c for particle

    }

    //Inelastic interaction
    if (random() <= DZ / DLY) {
        //cout<<" Absorbed   !! "<< endl;
    }
    //cout<<"XP APRES : "<<xp<<" YP APRES : "<<yp<<endl;
    xp = xp / 1000;
    yp = yp / 1000;
}


double SimCrys::random()
{
    double scale = RAND_MAX;
    double base = rand() / scale;
    double fine = rand() / scale;
    return base + fine / scale;
}

double SimCrys::random_dist()
{
    double scale = RAND_MAX;
    double base = rand() / scale;
    double fine = rand() / scale;
    return (-2.5 * 1e-4  + (3.5 * 1e-4 ) * (base + fine / scale));
}

double SimCrys::random_gauss()
{
    const double PI (4.0 * atan(1.0));
    double U(random());
    double V(random());
    return (sqrt(-2 * log(U)) * cos(2 * PI * V));
}

double SimCrys::random_gauss_cut()
{
    double resu;
    do {
        resu = random_gauss();
    } while (abs(resu) >= 1);
    return resu;
}

double SimCrys::random_gauss_dist()                     //POUR faire varier la moyenne de la gaussienne il faut mettre en argu la sim& !!!!!
{
    const double PI (4.0 * atan(1.0));
    double U(random());
    double V(random());
    //moy=random();
    return 1e-4 * (sqrt(-2 * log(U)) * cos(2 * PI * V)); //+(0.0002*moy);                   //POUR faire varier la moyenne de la gaussienne aussi !!!!!
}

/////////////////////////////////*******************************************************+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//POUR 409 :

double SimCrys::random_gauss_dist_409()                     //POUR faire varier la moyenne de la gaussienne il faut mettre en argu la sim& !!!!!
{
    const double PI (4.0 * atan(1.0));
    double U(random());
    double V(random());
    //moy=random();
    return (8.939203e-06) * (sqrt(-2 * log(U)) * cos(2 * PI * V)) + (-7.183296e-08);       //POUR faire varier la moyenne de la gaussienne aussi !!!!!
}

///////////////////////////////////************************************************************++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

void SimCrys::affiche()
{
    cout << "Crystal Curved Length [m] :" << crys.Cry_length << endl;
    cout << "Crystal Curvature Radius [m] :" << crys.Rcurv << endl;
    cout << "Crystal Bending Angle [urad] :" << crys.cry_bend << " , " << crys.Cry_length / crys.Rcurv << endl;
    cout << "Critical Radius :" << crys.Rcrit << endl;
    cout << "Ratio :" << crys.ratio << endl;
    cout << "Critical Angle for straight Crystal :" << crys.xpcrit0 << endl;
    cout << "Critical Angle for curved Crystal :" << crys.xpcrit << endl;
    cout << "Angle average reflexion : " << part.Ang_avr << endl;
    cout << "Full Channeling : " << crys.Cry_length / crys.Rcurv << endl;
    cout << endl;
    cout << endl;
    cout << "N.AMORPH :" << part.namor << endl;
    cout << "N.DECHANNELING:" << part.ndechann << endl;
    cout << "N.CHANNELING :" << part.nchann << endl;
    cout << "N.OUT :" << part.nout << endl;
    cout << "N.VREFLECTION :" << part.nvrefl << endl;
    cout << "N.VCAPTURE:" << part.nvcapt << endl;
    cout << "N.ABSORBED: "<< part.nabsorbed << endl;

    //cout << "Si la somme ne donne pas le nbre de part., c est normale !! En effet, il est possible que plusieurs evenement se produisent!!!";
}

void SimCrys::file_out(const int& pas, string outputpath)
{

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

    out.open(file.c_str(), ios::out | ios::app);

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

        out << endl;
        out << "Crystal Curved Length [m] :" << crys.Cry_length << endl;
        out << "Crystal Curvature Radius [m] :" << crys.Rcurv << endl;
        out << "Crystal Bending Angle [urad] :" << crys.cry_bend << " , " << 1e6*crys.Cry_length / crys.Rcurv << endl;
        out << "Critical Radius :" << crys.Rcrit << endl;
        out << "Ratio :" << crys.ratio << endl;
        out << "Critical Angle for straight Crystal :" << crys.xpcrit0 << endl;
        out << "Critical Angle for curved Crystal :" << crys.xpcrit << endl;
        out << "Angle average reflexion : " << part.Ang_avr << endl;
        out << "Full Channeling : " << crys.Cry_length / crys.Rcurv << endl;
        out << endl;
        out << endl;
        out << "N.AMORPH :" << part.namor << endl;
        out << "N.DECHANNELING:" << part.ndechann << endl;
        out << "N.CHANNELING :" << part.nchann << endl;
        out << "N.OUT :" << part.nout << endl;
        out << "N.VREFLECTION :" << part.nvrefl << endl;
        out << "N.VCAPTURE :" << part.nvcapt << endl ;
	out << "N.ABSORBED: "<< part.nabsorbed << endl;
	out << "N. of particles that hit the crystal collimator: " << part.nhit << endl;
	cout << endl;
	cout << endl;
	
        out.close();
    }

}