// CPP includes
#include <cstdlib>
#include <cstdio>
#include <iostream>
#include <cmath>
#include <vector>
#include <string>
#include <iostream>
#include <sstream>
#include <fstream>

// C includes
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

// ROOT includes
#include <TSystem.h>
#include <TROOT.h>
#include <TApplication.h>
#include <TRint.h>
#include <TDirectory.h>
#include <TFile.h>
#include <TChain.h>
#include <TTree.h>
#include <TNtuple.h>
#include <TBranch.h>
#include <TH1.h>
#include <TAxis.h>
#include <TCanvas.h>
#include <TColor.h>
#include <TPad.h>
#include <TStyle.h>
#include <TKey.h>
#include <TLegend.h>
#include <TText.h>
#include <TPaveText.h>
#include <TLine.h>
#include <TMarker.h>
#include <TLatex.h>
#include <TMath.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TF2.h>
#include <TString.h>
#include "TObjString.h"
#include "TRandom3.h"
#include "TGraph.h"
#include "TGraphSmooth.h"
#include "TSpline.h"

#include "./include/RungeKutta.hh"


using namespace std;

// Review the bunch period and the correlation with the position
// This code is tuned on a tollerance in the Runge Kutta method of 1.e-3.
// If you want to have a better tolerance, you need to change the code (namely
// to add a condition about rtol in the Runge-Kutta algorithm)



#define VERSION "v3.02-08-2016";

void run(int, double , double , double , double ,double ); 



void run(int tc, double ef, double ke, double sx, double sy,double sz){
    
  TH1::AddDirectory(kFALSE);
  TH1::AddDirectory(kFALSE);
 //=======================================================================================//
 //
 //                              DEFINE INITIAL VARIABLES:
 //
 // ======================================================================================//   
  
 // distance         = distance betwee the center of the IPM and the collection plate (meters)
 // c                = speed of light (m/s)
 // q                = elementary charge in Coulomb
 // m                = electron mass in Kg  
 // E                = applied electric field in V/m.  The field is along the y direction   
 // N_sigma_b        = NOT USED. Implemented to keep the same structure of Cyrille's code  
 // sigma_x          = bunch spread in x in lab sys in meters
 // sigma_y          = bunch spread in y in lab sys in meters
 // sigma_z          = bunch spread in z in lab sys in meters
 // sigma_bx         = bunch spread in x in bunch ref sys in meters
 // sigma_by         = bunch spread in y in bunch ref sys in meters
 // sigma_bz         = bunch spread in z in bunch ref sys in meters
 // Nb               = number of protons in the bunch
 // Ep0              = proton rest mass in MeV
 // Epi              = proton kinetic energy in MeV
 // dL               = longitudinal sampling region, i.e. the ionization charges 
 //                    are created in in a z poistion defined by 
 //                    -dL/2 <= z <= dL/2 
 // T                = bunch period
 // Qb               = unit charge of proton bunch
 // Q0               = unit charge of electrons   
 // dn               = dn/3 = dimension of the problem. In our case the  
 //                    problem is treated in 3-D, therefore dn = 9.
 //                    This means that you start by creating ionization
 //                    charges with an initial x, vx, ax, y, vy, ay, z, vz,
 //                    az and in the end of the code you will get the final
 //                    values of these same quantities   
 // N                = number of test particles followed in every loop. Here
 //                    the particles are followed 1 by one from creation to
 //                    the end of their history (when they reach a screen).  
 // tot_nb_test_part = total number of test particles created
 // MV               = MegaVolt unit 
 // mm               = scaling factor from mm to m   
 // scm              = scaling factor from meters to mm
 // GHz              = unit definition wrt Hz
 // ns               = unit definition wrt seconds
 // n                = number of time stamps used to solve the ODE equation
 
 
 double distance             = 5e-02;
 double c                    = 3e8 ;    
 double q                    = 1.6e-19 ;                   
 double m;
 if(tc<0) m                  = 9.10938356e-31;  //<--when the test charge is negative, you are following electrons, therefore m = mass of electron
 if(tc>0) m                  = 3.34605901e-27;  //<--when the test charge is negative, you are following H2+, therefore m = mass of electron
 double E;                   // Electric field between electrodes, when externally given //<-- Negative electric field if electrons are tracked as test charges (tc). Positive if positive test charges are tracked
 if(tc>0) E                  = abs(ef);                
 if(tc<0) E                  = ef;                
 double N_sigma_b            = 1 ;
 double sigma_x              = sx;
 double sigma_y              = sy; 
 double sigma_z              = sz; 
 double Nb                   = 1.1e+9;  
 double Ep0                  = 938 ;  
 double Epi                  = ke;      
 double gamma                = 1 + Epi / Ep0 ;
 double beta                 = sqrt(1 - 1./(gamma*gamma)) ;
 double dL                   = 5e-3 ;
 double T                    = 2.84e-9;
 double Qb                   = +1;
 double Q0                   = tc;  
 int dn                      = 9;
 int N                       = 1;
 double tot_nb_test_part     = 10;//0000;
 double MV                   = 1e6;
 double mm                   = 1e-3;
 double scm                  = 1e3; 
 double GHz                  = 1e9;
 double ns                   = 1e-9;
 TRandom3 r;
 
 // In the vectors of vectors you will store the solution of the differential equation for every particle and time
 // Using vector of vectors because you do not need to say in advance how many elements every vector will contain,
 // because you do not need to initialize elements to 0 and because you do not need to define as "new", with the
 // risk that you forget to delete it. 

 vector<vector<double> > history_ps_x_all_part;
 vector<vector<double> > history_vl_x_all_part;
 vector<vector<double> > history_ac_x_all_part;
 vector<vector<double> > history_ps_y_all_part;
 vector<vector<double> > history_vl_y_all_part;
 vector<vector<double> > history_ac_y_all_part;
 vector<vector<double> > history_ps_z_all_part;
 vector<vector<double> > history_vl_z_all_part;
 vector<vector<double> > history_ac_z_all_part;
 vector<vector<double> > history_tm_all_part;
 vector<vector<double> > history_Field_x_all_part;
 vector<vector<double> > history_Field_y_all_part;
 vector<vector<double> > history_Field_z_all_part;
 vector<vector<double> > history_Field_y_ext_all_part; 
 vector<vector<double> > history_Field_ps_z_bunch_all_part;
 
 
vector<double> history_ps_x;
vector<double> history_vl_x;
vector<double> history_ac_x;
vector<double> history_ps_y;
vector<double> history_vl_y;
vector<double> history_ac_y;
vector<double> history_ps_z;
vector<double> history_vl_z;
vector<double> history_ac_z;
vector<double> history_tm; 
vector<double> history_Field_x; 
vector<double> history_Field_y; 
vector<double> history_Field_z; 
vector<double> history_Field_y_ext;                            // "External" electric field, i.e. the one imposed between electrodes 
vector<double> history_ps_zbunch;     

 
 // Create the TFile and TTree where you are going to store the results 
 char fileName [200];
 TFile* f_res;
 if( Q0==-1) f_res = new TFile(Form("/local/home/fbelloni/WORK/MATLAB2C++/Step13TTree/tmp/Res_Electrons_%.1lf_field_%.1lf_MeV_%.2lf_x_%.2lf_y%.2lf_z.root",E,Epi,sigma_x*1000,sigma_y*1000,sigma_z*1000.),"recreate");
 if( Q0>0)   f_res = new TFile(Form("/local/home/fbelloni/WORK/MATLAB2C++/Step13TTree/tmp/Res_Ions_%.1lf_field_%.1lf_MeV_%.2lf_x_%.2lf_y%.2lf_z.root",E,Epi,sigma_x*1000,sigma_y*1000,sigma_z*1000.),"recreate");
 
 TTree *ttree = new TTree("tvec","Tree with vectors");
 ttree->Branch("x",&history_ps_x);
 ttree->Branch("y",&history_ps_y);
 ttree->Branch("z",&history_ps_z);
 ttree->Branch("vx",&history_vl_x);
 ttree->Branch("vy",&history_vl_y);
 ttree->Branch("vz",&history_vl_z);
 ttree->Branch("ax",&history_ac_x);
 ttree->Branch("ay",&history_ac_y);
 ttree->Branch("az",&history_ac_z);
 ttree->Branch("time",&history_tm);
 ttree->Branch("Field_x",&history_Field_x);
 ttree->Branch("Field_y",&history_Field_y);
 ttree->Branch("Field_z",&history_Field_z);
 ttree->Branch("Field_y_ext",&history_Field_y_ext);
 ttree->Branch("z_bunch",&history_ps_zbunch);

 
 // time                = time stamps at which you want the positions, speed, acceleration of the test particle and the field
 //                       it feels, to be stored. <-- changed. Now simply used to give a time and dt = 0 or different from 0 to Runge-Kutta
 // n                   = number of time stamps 
 // timen               = time stamps at which you want the positions, speed, acceleration of the test particle and the field
 //                       it feels, to be stored.
 // nn                  = number of time stamps 
 // dt                  = time step ("h" in the Runge Kutta solver) 
 // dtnew               = vector which will contain the new dt ( = h) 
 int    n                  = 87;//105;
 double *time              = new double[n]   {0 ,                    4.86782268547076e-21, 9.73564537094153e-21,  
                                              1.46034680564123e-20,  1.94712907418831e-20, 4.38104041692369e-20,
                                              6.81495175965907e-20,  9.24886310239445e-20, 1.16827744451298e-19,
                                              2.38523311588067e-19,  3.60218878724837e-19, 4.81914445861606e-19, 
                                              6.03610012998375e-19,  1.21208784868222e-18, 1.82056568436607e-18, 
                                              2.42904352004991e-18,  3.03752135573376e-18, 6.07991053415298e-18, 
                                              9.12229971257221e-18,  1.21646888909914e-17, 1.52070780694107e-17,
                                              3.04190239615068e-17,  4.56309698536029e-17, 6.08429157456991e-17,
                                              7.60548616377952e-17,  1.52114591098276e-16, 2.28174320558756e-16,
                                              3.04234050019237e-16,  3.80293779479718e-16, 7.60592426782121e-16,
                                              1.14089107408452e-15,  1.52118972138693e-15, 1.90148836868933e-15,
                                              3.80298160520135e-15,  5.70447484171336e-15, 7.60596807822538e-15,
                                              9.50746131473739e-15,  1.90149274972975e-14, 2.85223936798576e-14,
                                              3.80298598624176e-14,  4.75373260449777e-14, 9.50746569577781e-14,
                                              1.42611987870579e-13,  1.90149318783379e-13, 2.37686649696179e-13,
                                              4.75373304260181e-13,  7.13059958824183e-13, 9.50746613388185e-13,
                                              1.18843326795219e-12,  2.37686654077220e-12, 3.56529981359221e-12,
                                              4.75373308641222e-12,  5.94216635923223e-12, 1.18843327233323e-11,
                                              1.78264990874323e-11,  2.37686654515324e-11, 2.97108318156324e-11,
                                              5.94216636361327e-11,  8.91324954566329e-11, 1.18843327277133e-10,
                                              1.48554159097633e-10,  1.85053435546994e-10, 2.21552711996355e-10,
                                              2.58051988445716e-10,  2.94551264895077e-10, 3.31050541344438e-10,
                                              3.67549817793799e-10,  4.04049094243160e-10, 4.40548370692521e-10,
                                              4.93693760207255e-10,  5.46839149721989e-10, 5.99984539236723e-10,
                                              6.53129928751457e-10,  7.22238913724113e-10, 7.91347898696770e-10,
                                              8.60456883669427e-10,  9.29565868642083e-10, 1.03588250468420e-09,
                                              1.14219914072631e-09,  1.24851577676842e-09, 1.35483241281053e-09,
                                              1.47718121948677e-09,  1.59953002616300e-09, 1.72187883283923e-09,
                                              1.84422763951546e-09,  2.07104290640650e-09, 2.84e-09};//
                                              //2.29785817329753e-09,
                                              //2.52467344018857e-09,  2.75148870707960e-09, 3.05202392781177e-09,
                                              //3.35255914854393e-09,  3.65309436927610e-09, 3.95362959000827e-09,
                                              //4.03366776835209e-09,  4.11370594669591e-09, 4.19374412503973e-09,
                                              //4.27378230338355e-09,  5.e-09,               5.5e-09,
                                              //6e-09,                 7e-09,                8.e-09,
                                              //9.5e-09,               1.5e-08,              2.5e-08};
 
 int    nn                  = 30;//105;
 double *timen              = new double[n]   {0 ,                     
                                              1.46034680564123e-20,  
                                              6.81495175965907e-20,  
                                              2.38523311588067e-19,  
                                              6.03610012998375e-19,  
                                              2.42904352004991e-18,  
                                              9.12229971257221e-18,  
                                              3.04190239615068e-17,  
                                              7.60548616377952e-17,  
                                              3.04234050019237e-16,  
                                              1.14089107408452e-15,  
                                              3.80298160520135e-15,  
                                              9.50746131473739e-15,  
                                              3.80298598624176e-14,  
                                              1.42611987870579e-13,  
                                              4.75373304260181e-13,  
                                              1.18843326795219e-12,  
                                              4.75373308641222e-12,  
                                              1.78264990874323e-11,  
                                              5.94216636361327e-11,  
                                              1.48554159097633e-10,  
                                              2.58051988445716e-10,  
                                              3.67549817793799e-10,  
                                              4.93693760207255e-10,  
                                              6.53129928751457e-10,  
                                              8.60456883669427e-10,  
                                              1.14219914072631e-09,  
                                              1.47718121948677e-09,  
                                              1.84422763951546e-09,   
                                              2.84e-09};//

 long double *dt= new long double[n];
 dt[0]=0;
 for(int i=1;i<n;i++) dt[i]=time[i]-time[i-1]; 
 double tf=time[n-1];
 long double tdir=1;
 vector<double> dtnew;
 

 
//FILE* fff;
//fff = fopen("./tmp/Data.dat","w");
 
// Create the binary file where you will store the info you want. 
char BinFileName[256];
sprintf(BinFileName,"/local/home/fbelloni/WORK/MATLAB2C++/Step13TTree/tmp/Res_Electrons_%.1lf_field_%.1lf_MeV_%.2lf_x_%.2lf_y%.2lf_z.dat",E,Epi,sigma_x*1000,sigma_y*1000,sigma_z*1000.);
FILE *binaryfile = fopen(BinFileName,"wb");
double info[16];     
     
 // =============== CREATE AND FOLLOW PARTICLES ====================
 //
 // Loop over the number of particles
  for (int test=1; test<=tot_nb_test_part; test++)
   {
       dtnew.clear();
       cout << " particle nb " << test << endl;
    // initial beam size for the ODE (sets the distribution of ions/electrons in horizontal and vertical directions (x,y)
    // Convert the bunch width from lab system to bunch ref system 
    double sigma_bx=sigma_x; 
    double sigma_by=sigma_y; 
    double sigma_bz=sigma_z*gamma;
    // Beam kinetic energy  
    double Ep = Epi; 
    // Prepare a vector of 9 elements, ready to accept in the following order: x, vx, ax, y, vy, ay, z, vz, az
    double *P_init = new double[dn*N]; 
    for(int i=0;i<dn*N;i++) P_init[i]=0;
    // initial x and y position of the test charge
   for(int i =0;i<dn*N;i=i+dn) P_init[i] = sigma_x*r.Gaus(0.,1.);//(rand()/(float)RAND_MAX-1/2.);
    for(int i =3;i<dn*N;i=i+dn) P_init[i] = sigma_y*r.Gaus(0.,1.);//(rand()/(float)RAND_MAX-1/2.);
    // initial z position of the test charge 
    if (dn==9){for(int i =6;i<dn*N;i=i+dn) 
                 {P_init[i] =  dL*(rand()/(float)RAND_MAX-0.5);}//-6.0695e-04;}// dL*(rand()/(float)RAND_MAX-0.5);}
              }           
           

            
    // Some more initialisation for the analysis:
    
    
    // E_in_bunch_sys      = array which will contain the x, y and z component of the electric field in the bunch system
    // TotField_in_lab_sys = array which will contain the x, y and z component of the elm field in the lab reference frame
    // time_start          = since the time goes to 0 to a maximium and then starts over again looping from 0 to the same maxiumum value 
    //                       (periodic trend), time start will contain the period value times the loop number;
    // loop                = loop nuber, i.e. how many periods already passed 
    // t                   = time. It takes a periodic value. It is used only to have a t == 0 every T ns. 
    // t_time              = time variable (its value increases with the time passing by)
    // y                   = array which will contain P_init as the time passes by
    // tmp                 = array which wil contain the new time and dt. tmp[0] = dt, tmp[1] = t
     
    double* E_in_bunch_sys= new double[3];
    double* TotField_in_lab_sys= new double[5];
    double time_start         = 0;
    int    loop               = 0;
    double t                  = 0;
    double *y                 = new double [9]{P_init[0],P_init[1],P_init[2],P_init[3],P_init[4],P_init[5],P_init[6],P_init[7],P_init[8]};
    int k=0;
   


     long double *tmp = new long double[2];
     tmp[1]=0;
     while(y[3]<=distance){
    
      t=time[k]+time_start;     
      rk4(y, dt[k], dt[1], t, sigma_bx, sigma_by, sigma_bz, gamma,beta, N, Qb, Nb, Q0, m, E, E_in_bunch_sys,TotField_in_lab_sys, T, q, c,tf,  tdir,dtnew, tmp);
      dtnew.push_back(tmp[0]);
      history_ps_x.push_back (y[0]);
      history_vl_x.push_back (y[1]);
      history_ac_x.push_back (y[2]);
      history_ps_y.push_back (y[3]);
      history_vl_y.push_back (y[4]);
      history_ac_y.push_back (y[5]);
      history_ps_z.push_back (y[6]);
      history_vl_z.push_back (y[7]);
      history_ac_z.push_back (y[8]);
      history_tm.push_back (tmp[1]);
      history_Field_x.push_back(TotField_in_lab_sys[0]); 
      history_Field_y.push_back(TotField_in_lab_sys[1]); 
      history_Field_z.push_back(TotField_in_lab_sys[2]);
      history_Field_y_ext.push_back(TotField_in_lab_sys[3]);       
      history_ps_zbunch.push_back(TotField_in_lab_sys[4]);      
      //fprintf(fff,"%Le    %lf \n",tmp[1],y[3]);

      k++;
      if(k==n) {  k=0;
                  loop++;
                  time_start=loop*time[n-1];
          }  
          
       }
       // Here you pass the values you stored in vectors into array to be able to make a Graph and interpolate at the times "time[k]"
       int length=history_tm.size();
       double* x_val         = new double[length+1];
       double* vx_val        = new double[length+1];
       double* ax_val        = new double[length+1];
       double* y_val         = new double[length+1];
       double* vy_val        = new double[length+1];
       double* ay_val        = new double[length+1];
       double* z_val         = new double[length+1];
       double* vz_val        = new double[length+1];
       double* az_val        = new double[length+1];
       double* XField_val    = new double[length+1];
       double* YField_val    = new double[length+1];
       double* ZField_val    = new double[length+1];
       double* YExtField_val = new double[length+1];       
       double* zbunch_val    = new double[length+1];       
       double* vremia        = new double[length+1];       
       
       for(int i =0;i<length; i++){x_val[i]=history_ps_x.at(i);
                                   vx_val[i]=history_vl_x.at(i);
                                   ax_val[i]=history_ac_x.at(i);
                                   y_val[i]=history_ps_y.at(i);
                                   vy_val[i]=history_vl_y.at(i);
                                   ay_val[i]=history_ac_y.at(i);
                                   z_val[i]=history_ps_z.at(i);
                                   vz_val[i]=history_vl_z.at(i);
                                   az_val[i]=history_ac_z.at(i);
                                   vremia[i]=history_tm.at(i);
                                   XField_val[i]=history_Field_x.at(i);
                                   YField_val[i]=history_Field_y.at(i);
                                   ZField_val[i]=history_Field_z.at(i);
                                   YExtField_val[i]=history_Field_y_ext.at(i);
                                   zbunch_val[i] =history_ps_zbunch.at(i); 
                                   //cout << vremia[i] << " " << x_val[i]<< " " << y_val[i]<<endl;
                                   }
        history_ps_x.clear();
        history_vl_x.clear();
        history_ac_x.clear();
        history_ps_y.clear();
        history_vl_y.clear();
        history_ac_y.clear();
        history_ps_z.clear();
        history_vl_z.clear();
        history_ac_z.clear();
        history_tm.clear();
        history_Field_x.clear();
        history_Field_y.clear();
        history_Field_z.clear();
        history_Field_y_ext.clear();
        history_ps_zbunch.clear();
        double time_point=0;

        
        k=1;
        time_start=0;
        while(time_point<vremia[length-1])
          { 
            history_tm.push_back(time_point);
            time_point= time_start + timen[k];
            k++;
            if(k==n) {  k=0;
                       loop++;
                       time_start=loop*time[n-1];}
                       
          }      
        // TSpline does not work if, in the x & y data you want to fit, there are two y with the same value. 
        // For this reason you need to use TSmooth
        TGraph*g = new TGraph(length, y_val,vremia);  
        TGraph *grout=new TGraph(length);
        TGraphSmooth *gs = new TGraphSmooth("g");
        grout = gs->Approx(g,"linear",length);
        TSpline3 *s = new TSpline3("s",grout);
        vremia[length]=s->Eval(distance);
        if(s->Eval(distance)>history_tm.at(history_tm.size()-1)) history_tm.push_back(s->Eval(distance));
        else (history_tm.at(history_tm.size()-1) = s->Eval(distance));

        
        TGraph*g1 = new TGraph(length, vremia, x_val);
        TGraph *grout1=new TGraph(length);
        TGraphSmooth *gs1 = new TGraphSmooth("g1");
        grout1 = gs1->Approx(g1,"linear",length);
        TSpline5*s1 = new TSpline5("s1",grout1);
        
 
       
        TGraph*g2 = new TGraph(length, vremia,vx_val);
        TGraph *grout2=new TGraph(length);
        TGraphSmooth *gs2 = new TGraphSmooth("g2");
        grout2 = gs2->Approx(g2,"linear",length);
        TSpline5*s2 = new TSpline5("s2",grout2);

        
        TGraph*g3 = new TGraph(length, vremia,ax_val);
        TGraph *grout3=new TGraph(length);
        TGraphSmooth *gs3 = new TGraphSmooth("g3");
        grout3 = gs3->Approx(g3,"linear",length);
        TSpline5*s3 = new TSpline5("s3",grout3);
        
                
        TGraph*g4 = new TGraph(length, vremia,y_val);
        TGraph *grout4=new TGraph(length);
        TGraphSmooth *gs4 = new TGraphSmooth("g4");
        grout4 = gs4->Approx(g4,"linear",length);
        TSpline5*s4 = new TSpline5("s4",grout4);
        

        TGraph*g5 = new TGraph(length, vremia,vy_val);
        TGraph *grout5=new TGraph(length);
        TGraphSmooth *gs5 = new TGraphSmooth("g5");
        grout5 = gs5->Approx(g5,"linear",length);
        TSpline5*s5 = new TSpline5("s5",grout5);


        TGraph*g6 = new TGraph(length, vremia,ay_val);
        TGraph *grout6=new TGraph(length);
        TGraphSmooth *gs6 = new TGraphSmooth("g6");
        grout6 = gs6->Approx(g6,"linear",length);
        TSpline5*s6 = new TSpline5("s6",grout6);

  
        TGraph*g7 = new TGraph(length, vremia,z_val);
        TGraph *grout7=new TGraph(length);
        TGraphSmooth *gs7 = new TGraphSmooth("g7");
        grout7 = gs7->Approx(g7,"linear",length);
        TSpline5*s7 = new TSpline5("s7",grout7);


        TGraph*g8 = new TGraph(length, vremia,vz_val);
        TGraph *grout8=new TGraph(length);
        TGraphSmooth *gs8 = new TGraphSmooth("g8");
        grout8 = gs8->Approx(g8,"linear",length);
        TSpline5*s8 = new TSpline5("s8",grout8);


        
        TGraph*g9 = new TGraph(length, vremia,az_val);
        TGraph *grout9=new TGraph(length);
        TGraphSmooth *gs9 = new TGraphSmooth("g9");
        grout9 = gs9->Approx(g9,"linear",length);
        TSpline5*s9 = new TSpline5("s9",grout9);

        
        TGraph*g10 = new TGraph(length, vremia,XField_val);
        TGraph *grout10=new TGraph(length);
        TGraphSmooth *gs10 = new TGraphSmooth("g10");
        grout10 = gs10->Approx(g10,"linear",length);
        TSpline5*s10 = new TSpline5("s10",grout10);

        TGraph*g11 = new TGraph(length, vremia,YField_val);
        TGraph *grout11=new TGraph(length);
        TGraphSmooth *gs11 = new TGraphSmooth("g11");
        grout11 = gs11->Approx(g11,"linear",length);
        TSpline5*s11 = new TSpline5("s11",grout11);
        
          
        TGraph*g12 = new TGraph(length, vremia,ZField_val);
        TGraph *grout12=new TGraph(length);
        TGraphSmooth *gs12 = new TGraphSmooth("g12");
        grout12 = gs12->Approx(g12,"linear",length);
        TSpline5*s12 = new TSpline5("s12",grout12);


        TGraph*g13 = new TGraph(length, vremia,YExtField_val);
        TGraph *grout13=new TGraph(length);
        TGraphSmooth *gs13 = new TGraphSmooth("g13");
        grout13 = gs13->Approx(g13,"linear",length);
        TSpline5*s13 = new TSpline5("s13",grout13);

        
        TGraph*g14 = new TGraph(length, vremia,zbunch_val);
        TGraph *grout14=new TGraph(length);
        TGraphSmooth *gs14 = new TGraphSmooth("g14");
        grout14 = gs14->Approx(g14,"linear",length);
        TSpline5*s14 = new TSpline5("s14",grout14);



        
        for (int i=0;i<history_tm.size();i++) 
                                {history_ps_x.push_back(s1->Eval(history_tm.at(i)));
                                 history_vl_x.push_back(s2->Eval(history_tm.at(i)));
                                 history_ac_x.push_back(s3->Eval(history_tm.at(i)));     
                                 history_ps_y.push_back(s4->Eval(history_tm.at(i)));   
                                 history_vl_y.push_back(s5->Eval(history_tm.at(i)));
                                 history_ac_y.push_back(s6->Eval(history_tm.at(i)));             
                                 history_ps_z.push_back(s7->Eval(history_tm.at(i)));   
                                 history_vl_z.push_back(s8->Eval(history_tm.at(i)));
                                 history_ac_z.push_back(s9->Eval(history_tm.at(i)));  
                                 history_Field_x.push_back(s10->Eval(history_tm.at(i)));   
                                 history_Field_y.push_back(s11->Eval(history_tm.at(i)));
                                 history_Field_z.push_back(s12->Eval(history_tm.at(i)));                                        
                                 history_Field_y_ext.push_back(s13->Eval(history_tm.at(i)));
                                 history_ps_zbunch.push_back(s14->Eval(history_tm.at(i)));
                                 info[0]= test;
                                 info[1] = history_tm.at(i);
                                 info[2]=s1->Eval(history_tm.at(i));
                                 info[3]=s2->Eval(history_tm.at(i));
                                 info[4]=s3->Eval(history_tm.at(i));
                                 info[5]=s4->Eval(history_tm.at(i));
                                 info[6]=s5->Eval(history_tm.at(i));
                                 info[7]=s6->Eval(history_tm.at(i));
                                 info[8]=s7->Eval(history_tm.at(i));
                                 info[9]=s8->Eval(history_tm.at(i));
                                 info[10]=s9->Eval(history_tm.at(i));
                                 info[11]=s10->Eval(history_tm.at(i));
                                 info[12]=s11->Eval(history_tm.at(i));
                                 info[13]=s12->Eval(history_tm.at(i));
                                 info[14]=s13->Eval(history_tm.at(i));
                                 info[15]=s14->Eval(history_tm.at(i));
                                 //printf("%e %e  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e\n",history_tm.at(i), info[0], info[1], info[2], info[3], info[4], info[5],info[6], info[7], info[8], info[9], info[10], info[11], info[12], info[13]);
                                 fwrite(info,sizeof(double),16,binaryfile); 
                                
        }
        
        ttree->Fill();   
        delete g;
        delete grout;
        delete s;
        //delete gs;
        
 

        
                 
        delete g1;
        delete grout1;
        //delete gs1;
        delete s1; 
        delete g2;
        delete grout2;
        //delete gs2;
        delete s2;         
        delete g3;
        delete grout3;
        //delete gs3;
        delete s3;                 
        delete g4;
        delete grout4;
        //delete gs4;
        delete s4;
        delete g5;
        delete grout5;
        //delete gs5;
        delete s5; 
    
        delete g6;
        delete grout6;
        //delete gs6;
        delete s6;
        delete g7;
        delete grout7;
        //delete gs7;
        delete s7; 
        delete g8;
        delete grout8;
        //delete gs8;
        delete s8;         
        delete g9;
        delete grout9;
        //delete gs9;
        delete s9;                 
        delete g10;
        delete grout10;
        //delete gs10;
        delete s10;
        delete g11;
        delete grout11;
        //delete gs11;
        delete s11;
        delete g12;
        delete grout12;
        //delete gs12;
        delete s12;
        delete g13;
        delete grout13;
        //delete gs13;
        delete s13; 
        delete g14;
        delete grout14;
        //delete gs14;
        delete s14;
        
                         

        delete x_val;       
        delete vx_val;        
        delete ax_val;        
        delete y_val;         
        delete vy_val;        
        delete ay_val;        
        delete z_val;         
        delete vz_val;        
        delete az_val;        
        delete XField_val;    
        delete YField_val;    
        delete ZField_val;    
        delete YExtField_val; 
        delete zbunch_val;    
        delete vremia;         
        
        
        history_ps_x_all_part.push_back(history_ps_x);
        history_vl_x_all_part.push_back(history_vl_x);
        history_ac_x_all_part.push_back(history_ac_x);
        history_ps_y_all_part.push_back(history_ps_y);
        history_vl_y_all_part.push_back(history_vl_y);
        history_ac_y_all_part.push_back(history_ac_y);
        history_ps_z_all_part.push_back(history_ps_z);
        history_vl_z_all_part.push_back(history_vl_z);
        history_ac_z_all_part.push_back(history_ac_z);
        history_tm_all_part.push_back(history_tm);
        history_Field_x_all_part.push_back(history_Field_x); 
        history_Field_y_all_part.push_back(history_Field_y); 
        history_Field_z_all_part.push_back(history_Field_z); 
        history_Field_y_ext_all_part.push_back(history_Field_y_ext); 
        history_Field_ps_z_bunch_all_part.push_back(history_ps_zbunch);       

        delete y;

        delete E_in_bunch_sys;
        delete TotField_in_lab_sys;
        delete P_init;
        history_ps_x.clear();
        history_vl_x.clear();
        history_ac_x.clear();
        history_ps_y.clear();
        history_vl_y.clear();
        history_ac_y.clear();
        history_ps_z.clear();
        history_vl_z.clear();
        history_ac_z.clear();
        history_tm.clear(); 
        history_Field_x.clear(); 
        history_Field_y.clear(); 
        history_Field_z.clear(); 
        history_Field_y_ext.clear(); // External, imposed in the chamber
        history_ps_zbunch.clear();  
        delete tmp;
    
   }
   fclose(binaryfile);
   delete time;
   delete timen;
   delete dt;
   f_res->cd();
   ttree->Write();
//   f_res->Write();
   delete f_res;
}





//==========




//=====

int main(int argc, char **argv)
{
  int    tc             = atoi(argv[1]);     // test charge
  double ef             = atof(argv[2]);   // electric field
  double ke             = atof(argv[3]);   // proton kinetic energy
  double sx             = atof(argv[4]);
  double sy             = atof(argv[5]);
  double sz             = atof(argv[6]);
  
      
  TRint *theApp = new TRint("App", &argc, argv, NULL, 0);
  run(tc, ef, ke, sx, sy,sz); 
  //process_one_raw_data_file(TString::Format("rfio:/castor/cern.ch/ntof/%s/stream%d/run%d_%d_s%d.raw.finished",rawdata_dir, stream, nrunr, segm, stream), detector_name_list, Run_Type, output_dir, XYConfig, nrunr, segm, stream);
  //run();
  theApp->Run();
  return 0;

    }