// Simplified calculation of magnetic field effect on propagation of particle with 
// wrong momentum assumption (background) in PANDA mag.field
// Based on Runge-Kutta algorithm with 4 order 
// [refence to paper by I.Kisel & S.Gorbunov about this topic:
// http://www-linux.gsi.de/~ikisel/reco/CBM/DOC-2005-Mar-212-1.pdf]

// v.1.0: Bx=Bz=0, By!=0
// date: 23/09/2013
// v.2.0: Bx=Bz=0, By!=0, bug fixed
//date: 25/11/2013

// code author: A. Karavdina

#include "TMath.h"
#include "TF1.h"
#include "TVector3.h"
#include <iostream>
#include <fstream>
#include "TH2D.h"
#include "TTree.h"
#include "TFile.h"
#include "TString.h"
#include "TRandom3.h"
#include "TGraph.h"
#include "TNtuple.h"

using namespace std;
const double kappa = 2.99792458e-4; // [(GeV/c)/(kG^{−1}*cm^{−1})] 
//Description for params in functions ay&by:
//p[0] = momentum [GeV]
//p[1] = charge [e]
//x[0] = t_x, direction in x
//x[1] = t_y, direction in y

//TTree *Byintegttbl;
double gl_z1[2641500],gl_z2[2641500],gl_By[2641500],gl_By1[2641500],gl_By2[2641500],gl_By3[2641500];
int stField;
//TTree *T;
int FieldTable(char inputFilename[500] ="15_GeV_FieldIntegTable.txt"){
ifstream inFile(inputFilename);

  // std::ifstream inFile;
  // inFile.open(inputFilename,std::ios::in | std::ios::binary);
  //  char inputFilename[] = fname;
  //  inFile.open(inputFilename, ios::in);
  // std::ifstream inFile(inputFilename, std::ifstream::in);
  int j=0;
  while (!inFile.eof()) {
    inFile >> gl_z1[j] >>  gl_z2[j]  >> gl_By[j]  >>  gl_By1[j]  >>  gl_By2[j]  >>  gl_By3[j] ;
//cout<<"gl_z1["<<j<<"] = "<<gl_z1[j] <<" gl_By3[j]"<<endl;
    j++;
  }
stField = j-1;

  // TTree *Tk = new TTree("Tk","data from ascii file");
  // Long64_t nlines = Tk->ReadFile(fname,"z1:z2:By:By1:By2:By3");
  // stField = nlines;
  // //  T->Write();
  // // T->Print();
  // Tk->Print();
  // // T->Draw("z1");
  // printf(" found %lld points\n",nlines);
  // TBranch *br_z1 = Tk->GetBranch("z1");
  // TBranch *br_z2 = Tk->GetBranch("z2");
  // TBranch *br_By = Tk->GetBranch("By"); //field values
  // TBranch *br_By1 = Tk->GetBranch("By1");//integral values
  // TBranch *br_By2 = Tk->GetBranch("By2");
  // TBranch *br_By3 = Tk->GetBranch("By3");
  // double z1in,z2in,Byin,By1in,By2in,By3in;
  // br_z1->SetAddress(&z1in);
  // br_z2->SetAddress(&z2in);
  // br_By->SetAddress(&Byin);
  // br_By1->SetAddress(&By1in);
  // br_By2->SetAddress(&By2in);
  // br_By3->SetAddress(&By3in);

  // // T->SetBranchAddress("z1",&z1in);
  // // T->SetBranchAddress("z2",&z2in);
  // // T->SetBranchAddress("By",&Byin);
  // // T->SetBranchAddress("By1",&By1in);
  // // T->SetBranchAddress("By2",&By2in);
  // // T->SetBranchAddress("By3",&By3in);
  

  // //  int nevent = T->GetEntries();
  // for(int itree=0;itree<stField;itree++){
  //   Tk->GetEvent(itree);
  //   //    br_z1->GetEvent(itree);
  //   gl_z1[itree] = z1in;
  //   gl_z2[itree] = z2in;
  //   cout<<" z1in = "<<z1in<<endl;
  //   cout<<"gl_z1[itree] = "<<gl_z1[itree]<<endl;
  //   gl_By[itree] = Byin;
  //   gl_By1[itree] = By1in;
  //   gl_By2[itree] = By3in;
  //   gl_By3[itree] = By3in;
  // }
  return 0;
}

// tx' = ay*By
double ay(double *x, double *p){ 
  double tx_2 = TMath::Power(x[0],2);
  double ty_2 = TMath::Power(x[1],2);
  return (-1)*TMath::Power(kappa*p[1]/p[0],1)*TMath::Power((1+tx_2+ty_2),0.5)*(1+tx_2);
}

//ty' = by*By
double by(double *x, double *p){ 
  double tx_2 = TMath::Power(x[0],2);
  double ty_2 = TMath::Power(x[1],2);
  return (-1)*TMath::Power(kappa*p[1]/p[0],1)*TMath::Power((1+tx_2+ty_2),0.5)*(x[0]*x[1]);
}


int By_read(double *z, double &By, double &By1, double &By2, double &By3){
  By = 0;
  By1 = 0;
  By2 = 0;
  By3 = 0;
  if( (z[0]>324.9 &&  z[0]<621.8) || (z[1]>324.9 &&  z[1]<621.8)){
    if(z[0]>324.9 &&  z[0]<621.8){
      if(z[1]<324.9) z[1]=324.9;
      if(z[1]>621.8) z[1]=621.8;
    }
    if(z[1]>324.9 &&  z[1]<621.8){
      if(z[0]<324.9) z[0]=324.9;
      if(z[0]>621.8) z[0]=621.8;
    }  

    double glz1 = gl_z1[0];
    double glz2 = gl_z2[0];
    int i=0;
    // cout<<"z[0] = "<<z[0]<<" z[1] = "<<z[1]<<" glz1 = "<<glz1<<" glz2 = "<<glz2<<endl;
    while(z[0]>glz1){
      i++;
      glz1 = gl_z1[i];
      glz2 = gl_z2[i];
    }

    double fBy = 0;
    double fBy1 = 0;
    double fBy2 = 0;
    double fBy3 = 0;
    while(z[0]<=glz1 && z[1]>glz2){
      //   cout<<"z[0] = "<<z[0]<<" z[1] = "<<z[1]<<" glz1 = "<<glz1<<" glz2 = "<<glz2<<endl;
      fBy += gl_By[i]; //not correct, but also doesn't used here
      fBy1 += gl_By1[i]; 
      fBy2 += gl_By2[i]; 
      fBy3 += gl_By3[i]; 
      i++;
      glz1 = gl_z1[i];
      glz2 = gl_z2[i];
    }

    // for(int i=0;i<stField;i++){
    //   cout<<"z[0] = "<<z[0]<<" z[1] = "<<z[1]<<" gl_z1[i] = "<<gl_z1[i]<<endl;
    //   if(z[0]<=gl_z1[i]){
    // 	cout<<"gl_z1[i] = "<<gl_z1[i]<<" gl_z2[i] = "<<gl_z2[i]<<endl;
    // 	int j=i;
    // 	while(z[1]>gl_z2[j]){
    // 	  cout<<"gl_z1[j] = "<<gl_z1[j]<<" gl_z2[j] = "<<gl_z2[j]<<endl;
    // 	  fBy += gl_By[j]; 
    // 	  fBy1 += gl_By1[j]; 
    // 	  fBy2 += gl_By2[j]; 
    // 	  fBy3 += gl_By3[j]; 
    // 	  j++;
    // 	  i++;
    // 	}
    // 	continue;
    //   }
    // }
    By = fBy;
    By1 = fBy1;
    By2 = fBy2;
    By3 = fBy3;
    // cout<<"between z[0]="<<z[0]<<" and z[1]="<<z[1]<<" By1 = "<<By1<<" By2 = "<<By2<<" By3 = "<<By3<<endl;
  }
  return 0;
}

//! very approximate approximation TODO: change it for better one
double By(double *z){
  if( (z[0]>324.9 &&  z[0]<621.8) || (z[1]>324.9 &&  z[1]<621.8)){
    if(z[0]>324.9 &&  z[0]<621.8){
      if(z[1]<324.9) z[1]=324.9;
      if(z[1]>621.8) z[1]=621.8;
    }
    if(z[1]>324.9 &&  z[1]<621.8){
      if(z[0]<324.9) z[0]=324.9;
      if(z[0]>621.8) z[0]=621.8;
    }



  
    //    double par[7] = {1022.18, -13.5362, 0.0734285, -0.000209178, 3.30658e-07, -2.75369e-10, 9.44655e-14};
    double par[7] = {10*1022.18, 10*-13.5362, 10*0.0734285, 10*-0.000209178, 10*3.30658e-07, 10*-2.75369e-10, 10*9.44655e-14};
    double fBy = 0;
    for(int i=0;i<7;i++)
      fBy +=par[i]*TMath::Power(z[0],i);
    //   cout<<"somebody called By("<<z[0]<<")="<<fBy<<endl;

    return fBy;
  }
  else return 0;
}

Double_t Bydraw(double *x, double *par) {
  double z1 = x[0]-par[0];
  double z2 = x[0]+par[0];
  double z[2]={z1,z2};
  return By(z);
}

void Draw_Field(){
  TF1 *fa3 = new TF1("fa3",Bydraw,0,1124,1);
  fa3->SetParameters(0,0.01);
  fa3->SetParNames("step");
  fa3->Draw();
}

double By_integ_1(double *z){
  if( (z[0]>324.9 &&  z[0]<621.8) || (z[1]>324.9 &&  z[1]<621.8)){
    if(z[0]>324.9 &&  z[0]<621.8){
      if(z[1]<324.9) z[1]=324.9;
      if(z[1]>621.8) z[1]=621.8;
    }
    if(z[1]>324.9 &&  z[1]<621.8){
      if(z[0]<324.9) z[0]=324.9;
      if(z[0]>621.8) z[0]=621.8;
    }
    // double par[7]={101.276, -1.324.931, 0.00652559,
    // 		   -1.54572e-05,  1.77642e-08, -7.99135e-12}; 
    //    double par[4]={-22.44, 0.133067,   -0.00024342, 1.4136e-07};
    double par[7] = {10*1022.18, 10*-13.5362, 10*0.0734285, 10*-0.000209178, 10*3.30658e-07, 10*-2.75369e-10, 10*9.44655e-14};
    double fBy = 0;
    for(int i=0;i<7;i++)
      fBy +=(par[i]/(i+1))*(TMath::Power(z[1],i+1)-TMath::Power(z[0],i+1));
    // if(z[0]==324.9 || z[0]==621.8 || z[1]==324.9 || z[1]==621.8)
    //   cout<<"somebody called By_integ_1("<<z[0]<<", "<<z[1]<<")="<<fBy<<endl;
    return fBy;
  }
 else return 0;
}


double By_integ_2(double *z){
  if( (z[0]>324.9 &&  z[0]<621.8) || (z[1]>324.9 &&  z[1]<621.8)){
    if(z[0]>324.9 &&  z[0]<621.8){
      if(z[1]<324.9) z[1]=324.9;
      if(z[1]>621.8) z[1]=621.8;
    }
    if(z[1]>324.9 &&  z[1]<621.8){
      if(z[0]<324.9) z[0]=324.9;
      if(z[0]>621.8) z[0]=621.8;
    }
    // double par[7]={101.276, -1.324.931, 0.00652559,
    // 		   -1.54572e-05,  1.77642e-08, -7.99135e-12}; 
    //    double par[4]={-22.44, 0.133067,   -0.00024342, 1.4136e-07};
    double par[7] = {10*1022.18, 10*-13.5362, 10*0.0734285, 10*-0.000209178, 10*3.30658e-07, 10*-2.75369e-10, 10*9.44655e-14};
    double fBy = 0;
    for(int i=0;i<7;i++)
      for(int j=0;j<7;j++)
	fBy +=(par[i]*par[j]/(j+1))*((TMath::Power(z[1],i+j+2)-TMath::Power(z[0],i+j+2))/(i+j+2) 
				     -(TMath::Power(z[0],j+1)*(TMath::Power(z[1],i+1)-TMath::Power(z[0],i+1))/(i+1)));
    // if(z[0]==324.9 || z[0]==621.8 || z[1]==324.9 || z[1]==621.8)
    //   cout<<"somebody called By_integ_2("<<z[0]<<", "<<z[1]<<")="<<fBy<<endl;
    //    cout<<"@z="<<z[0]<<"integ By="<<fBy<<endl;
    return fBy;
  }
  else return 0;
}

double By_integ_3(double *z){
  if( (z[0]>324.9 &&  z[0]<621.8) || (z[1]>324.9 &&  z[1]<621.8)){
    if(z[0]>324.9 &&  z[0]<621.8){
      if(z[1]<324.9) z[1]=324.9;
      if(z[1]>621.8) z[1]=621.8;
    }
    if(z[1]>324.9 &&  z[1]<621.8){
      if(z[0]<324.9) z[0]=324.9;
      if(z[0]>621.8) z[0]=621.8;
    }
    // double par[7]={101.276, -1.324.931, 0.00652559,
    // 		   -1.54572e-05,  1.77642e-08, -7.99135e-12}; 
    //    double par[4]={-22.44, 0.133067,   -0.00024342, 1.4136e-07};
    double par[7] = {10*1022.18, 10*-13.5362, 10*0.0734285, 10*-0.000209178, 10*3.30658e-07, 10*-2.75369e-10, 10*9.44655e-14};
    double fBy = 0;
    for(int i=0;i<7;i++)
      for(int j=0;j<7;j++)
	for(int k=0;k<7;k++)
	  fBy +=(par[i]*par[j]*par[k]/(j+1))*((TMath::Power(z[1],k+i+j+3)-TMath::Power(z[0],k+i+j+3))/((k+j+i+3)*(i+j+2))
					      -TMath::Power(z[0],i+j+2)*(TMath::Power(z[1],k+1)-TMath::Power(z[0],k+1))/((i+j+2)*(k+1))
					      -TMath::Power(z[0],j+1)*(TMath::Power(z[1],k+i+2)-TMath::Power(z[0],k+i+2))/((i+1)*(k+i+2))
					      +TMath::Power(z[0],i+j+2)*(TMath::Power(z[1],k+1)-TMath::Power(z[0],k+1))/((i+1)*(k+1)));
    // if(z[0]==324.9 || z[0]==621.8 || z[1]==324.9 || z[1]==621.8)
    //   cout<<"somebody called By_integ_3("<<z[0]<<", "<<z[1]<<")="<<fBy<<endl;
    return fBy;
  }
  else return 0;
}

double By_integ_4(double *z){
  return 0;

  // if( (z[0]>324.9 &&  z[0]<621.8) || (z[1]>324.9 &&  z[1]<621.8)){
  //   if(z[0]>324.9 &&  z[0]<621.8){
  //     if(z[1]<324.9) z[1]=324.9;
  //     if(z[1]>621.8) z[1]=621.8;
  //   }
  //   if(z[1]>324.9 &&  z[1]<621.8){
  //     if(z[0]<324.9) z[0]=324.9;
  //     if(z[0]>621.8) z[0]=621.8;
  //   }
  //   //    double par[4]={-22.44, 0.133067,   -0.00024342, 1.4136e-07};
  //   double par[7] = {10*1022.18, 10*-13.5362, 10*0.0734285, 10*-0.000209178, 10*3.30658e-07, 10*-2.75369e-10, 10*9.44655e-14};
  //   double fBy = 0;
  //   for(int i=0;i<7;i++)
  //     for(int j=0;j<7;j++)
  // 	for(int k=0;k<7;k++)
  // 	  for(int m=0;m<7;m++){
  // 	    int kmji4 = k+m+j+i+4;
  // 	    int kji3 = k+j+i+3;
  // 	    int ji2 = j+i+2;
  // 	    int km2 = k+m+2;
  // 	    int kmi3 = k+m+i+3;
  // 	    double A = (par[i]*par[j]*par[k]*par[m]/(j+1));
  // 	    //    cout<<"A = "<<A<<endl;
  // 	    //	    double B = (TMath::Power(z[1],kmji4)-TMath::Power(z[0],kmji4))/((kmji4)*(kji3)*(ji2));
  // 	    double C =  -TMath::Power(z[0],kji3)*(((TMath::Power(z[1],m+1)-TMath::Power(z[0],m+1))/((kji3)*(ji2)*(m+1))-(TMath::Power(z[1],m+1)-TMath::Power(z[0],m+1))/((ji2)*(k+1)*(m+1)))+(TMath::Power(z[1],m+1)-TMath::Power(z[0],m+1))/((i+1)*(k+1)*(m+1)))+(TMath::Power(z[1],kmji4)-TMath::Power(z[0],kmji4))/((kmji4)*(kji3)*(ji2));
  // 	    double D = -TMath::Power(z[0],ji2)*((TMath::Power(z[1],km2)-TMath::Power(z[0],km2))/((ji2)*(k+1)*(km2)) - (TMath::Power(z[1],km2)-TMath::Power(z[0],km2))/((i+1)*(k+1)*(km2))) -TMath::Power(z[0],j+1)*(TMath::Power(z[1],kmi3)-TMath::Power(z[0],kmi3))/((i+1)*(k+i+2)*(kmi3));
  // 	    //	    double F = TMath::Power(z[0],kji3)*(TMath::Power(z[1],m+1)-TMath::Power(z[0],m+1))/((ji2)*(k+1)*(m+1));
  // 	    // double G =-TMath::Power(z[0],j+1)*(TMath::Power(z[1],kmi3)-TMath::Power(z[0],kmi3))/((i+1)*(k+i+2)*(kmi3));
  // 	    //	    double E = TMath::Power(z[0],ji2)*(TMath::Power(z[1],km2)-TMath::Power(z[0],km2))/((i+1)*(k+1)*(km2));
  // 	    //	    double H = -TMath::Power(z[0],kji3)*(TMath::Power(z[1],m+1)-TMath::Power(z[0],m+1))/((i+1)*(k+1)*(m+1));

  // 	    //   cout<<" C = "<<C<<" D = "<<D<<endl;
  // 	    fBy +=A*C+A*D;
  // 	  }
  //   // if(z[0]==324.9 || z[0]==621.8 || z[1]==324.9 || z[1]==621.8)
  //   //   cout<<"somebody called By_integ_4("<<z[0]<<", "<<z[1]<<")="<<fBy<<endl;
  //   return fBy;
  // }
  // else return 0;
}



// //T functions of k-order; k=1
// double T1_x(double *x, double *p){
//   return ay(x,p);
// }
// double T1_y(double *x, double *p){
//   return by(x,p);
// }

// //T functions of k-order; k=2
// double T2_x(double *x, double *p){
//   double tx_2 = TMath::Power(x[0],2);
//   double ty_2 = TMath::Power(x[1],2);
//   return TMath::Power(kappa*p[1]/p[0],2)*(1+tx_2+ty_2)*(1+tx_2)*x[0];
// }
// double T2_y(double *x, double *p){
//   double tx_2 = TMath::Power(x[0],2);
//   double ty_2 = TMath::Power(x[1],2);
//   return TMath::Power(kappa*p[1]/p[0],2)*(1+tx_2+ty_2)*(1+tx_2)*x[1];
// }

// //T functions of k-order; k=3
// double T3_x(double *x, double *p){
//   double tx_2 = TMath::Power(x[0],2);
//   double ty_2 = TMath::Power(x[1],2);
//   return (-1)*TMath::Power(kappa*p[1]/p[0],3)*TMath::Power((1+tx_2+ty_2),3./2)*(1+tx_2)*(1+5*tx_2);
// }
// double T3_y(double *x, double *p){
//   double tx_2 = TMath::Power(x[0],2);
//   double ty_2 = TMath::Power(x[1],2);
//   return (-5)*TMath::Power(kappa*p[1]/p[0],3)*TMath::Power((1+tx_2+ty_2),3./2)*(1+tx_2)*x[0]*x[1];
// }

// //T functions of k-order; k=4
// double T4_x(double *x, double *p){
//   double tx_2 = TMath::Power(x[0],2);
//   double ty_2 = TMath::Power(x[1],2);
//   return 5*TMath::Power(kappa*p[1]/p[0],4)*TMath::Power((1+tx_2+ty_2),2)*(1+tx_2)*x[0]*(3+7*tx_2);
// }
// double T4_y(double *x, double *p){
//   double tx_2 = TMath::Power(x[0],2);
//   double ty_2 = TMath::Power(x[1],2);
//   return 5*TMath::Power(kappa*p[1]/p[0],4)*TMath::Power((1+tx_2+ty_2),2)*(1+7*tx_2)*x[1];
// }



//From I.Kisel paper --------------------------
//T functions of k-order; k=1
double T1_x(double *x, double *p){
  return ay(x,p);
}
double T1_y(double *x, double *p){
  return by(x,p);
}

//T functions of k-order; k=2
double T2_x(double *x, double *p){
  //cout<<"p[1] = "<<p[1]<<" p[0] = "<<p[0]<<endl;
  double tx_2 = TMath::Power(x[0],2);
  double ty_2 = TMath::Power(x[1],2);
  double h = (kappa*p[1]/p[0])*sqrt(1+tx_2+ty_2);
  return 3*TMath::Power(h,2)*(1+tx_2)*x[0];
}
double T2_y(double *x, double *p){
  double tx_2 = TMath::Power(x[0],2);
  double ty_2 = TMath::Power(x[1],2);
  double h = (kappa*p[1]/p[0])*sqrt(1+tx_2+ty_2);
  return TMath::Power(h,2)*(1+3*tx_2)*x[1];
}

//T functions of k-order; k=3
double T3_x(double *x, double *p){
  double tx_2 = TMath::Power(x[0],2);
  double ty_2 = TMath::Power(x[1],2);
  double h = (kappa*p[1]/p[0])*sqrt(1+tx_2+ty_2);
  return (-3)*TMath::Power(h,3)*(1+5*tx_2*tx_2+6*tx_2);
}
double T3_y(double *x, double *p){
  double tx_2 = TMath::Power(x[0],2);
  double ty_2 = TMath::Power(x[1],2);
 double h = (kappa*p[1]/p[0])*sqrt(1+tx_2+ty_2);
  return (-1)*TMath::Power(h,3)*(9+15*tx_2)*x[0]*x[1];
}

//extrapol(momentum,icharge,x1,y1,tx1,ty1,z1,z2,x2,y2,tx2,ty2);
int extrapol(double p, double q, double x0, double y0, double tx0, double ty0,double z1,double z2, double &x2, double &y2, double &tx2, double &ty2){
  double x[2]={tx0,ty0};
  double par[2]={p,q};
  int chsg = 1;
  if(z1>z2){
    chsg = -1;
    double z1_tmp = z1;     double z2_tmp = z2;
    z1 = z2_tmp;  z2 = z1_tmp;
  }
  double z[2]={z1,z2};
  //  double fdx_sum = tx0*(z2-z1)+T1_x(x,par)*By_integ_2(z)+ T2_x(x,par)*By_integ_3(z)+ T3_x(x,par)*By_integ_4(z);
  double fBy, fBy1, fBy2, fBy3;
  int okfl = By_read(z,fBy,fBy1,fBy2,fBy3);
  double dx0 = tx0*(z2-z1)+T1_x(x,par)*fBy2+ T2_x(x,par)*fBy3;
  double dy0 = ty0*(z2-z1) +T1_y(x,par)*fBy2+ T2_y(x,par)*fBy3;
  double dtx0 = T1_x(x,par)*fBy1+ T2_x(x,par)*fBy2+ T3_x(x,par)*fBy3;
  double dty0 = T1_y(x,par)*fBy1+ T2_y(x,par)*fBy2+ T3_y(x,par)*fBy3;
  x2 = x0+chsg*dx0;
  y2 = y0+chsg*dy0;
  tx2 = tx0 + chsg*dtx0;
  ty2 = ty0 + chsg*dty0;
}

double dx_sum(double p, double q, double tx0, double ty0,double z1,double z2){
  double x[2]={tx0,ty0};
  double par[2]={p,q};
  double z[2]={z1,z2};
  double fdx_sum = tx0*(z2-z1)+T1_x(x,par)*By_integ_2(z)+ T2_x(x,par)*By_integ_3(z)+ T3_x(x,par)*By_integ_4(z);
  //  cout<<" By_integ_2(z) = "<<By_integ_2(z)<<" By_integ_3(z) = "<<By_integ_3(z)<<endl;
  // double fBy, fBy1, fBy2, fBy3;
  // int okfl = By_read(z,fBy,fBy1,fBy2,fBy3);
  // double fdx_sum = tx0*(z2-z1)+T1_x(x,par)*fBy2+ T2_x(x,par)*fBy3;
  return fdx_sum;
}

double dy_sum(double p, double q, double tx0, double ty0,double z1,double z2){
  double x[2]={tx0,ty0};
  double par[2]={p,q};
  double z[2]={z1,z2};
  double fdy_sum = ty0*(z2-z1) +T1_y(x,par)*By_integ_2(z)+ T2_y(x,par)*By_integ_3(z)+ T3_y(x,par)*By_integ_4(z);
  // double fBy, fBy1, fBy2, fBy3;
  // int okfl = By_read(z,fBy,fBy1,fBy2,fBy3);
  // double fdy_sum = ty0*(z2-z1) +T1_y(x,par)*fBy2+ T2_y(x,par)*fBy3;
  return fdy_sum;
}
//[END] From I.Kisel paper --------------------------

double tx_sum(double p, double q, double tx0, double ty0,double z1,double z2){
  double x[2]={tx0,ty0};
  double par[2]={p,q};
  double z[2]={z1,z2};
  double ftx_sum = tx0 + T1_x(x,par)*By_integ_1(z)+ T2_x(x,par)*By_integ_2(z)+ T3_x(x,par)*By_integ_3(z);
  // double fBy, fBy1, fBy2, fBy3;
  // int okfl = By_read(z,fBy,fBy1,fBy2,fBy3);
  // double ftx_sum = tx0 + T1_x(x,par)*fBy1+ T2_x(x,par)*fBy2+ T3_x(x,par)*fBy3;
  return ftx_sum;
}

double ty_sum(double p, double q, double tx0, double ty0,double z1,double z2){
  double x[2]={tx0,ty0};
  double par[2]={p,q};
  double z[2]={z1,z2};
  double fty_sum = ty0 + T1_y(x,par)*By_integ_1(z) + T2_y(x,par)*By_integ_2(z) + T3_y(x,par)*By_integ_3(z);

  // double fBy, fBy1, fBy2, fBy3;
  // int okfl = By_read(z,fBy,fBy1,fBy2,fBy3);
  // double fty_sum = ty0 + T1_y(x,par)*fBy1+ T2_y(x,par)*fBy2+ T3_y(x,par)*fBy3;
  return fty_sum;
}


//Generate theta,phi,x,y initial parameters (@ IP)
void Generate_pars_IP(TRandom3 *rand,double& th_in, double& ph_in, double& x0, double& y0, double& z0, double thmin=3e-3, double thmax=2*TMath::Pi(), double phmin=0, double phmax=2*TMath::Pi()){
  //Sample uniform theta and phi  in given range  ---------------------------

  double posth = rand->Rndm(); // Produces uniformly-distributed floating points in ]0,1]
  double posph = rand->Rndm();
  double costhmin = TMath::Cos(thmax);
  double costhmax = TMath::Cos(thmin);
  double costhcur = costhmin+posth*(costhmax-costhmin);
  if(fabs(costhcur)>1.)  cout<<"costhcur = "<<costhcur<<endl;
  double thcur = TMath::ACos(costhcur);
  double phcur = phmin+posph*(phmax-phmin);
  // cout<<"rand th = "<< posth <<" rand ph = "<<posph<<endl;
  // cout<<" thcur = "<<thcur<<" phcur = "<<phcur<<endl;
  TVector3 dir(0,0,1);
  dir.SetTheta(thcur);
  dir.SetPhi(phcur);
  // if(phcur>0) ph_in =phcur-TMath::Pi();
  // if(phcur<0) ph_in =phcur+TMath::Pi();
  // cout<<"phcur = "<<phcur<<" ph_in = "<<ph_in<<endl;
  ph_in = phcur;
  th_in = thcur;
  x0 = 0.; y0 = 0.; z0 = 0.;
 
}

//Read results of previuos propagation from TTree
//check are params in expected range or not
int Read_propag_pars(bool xycut, bool thphcut, TString treeName, TFile *f_in, int itree, double &x_sim, double &y_sim, double &z_sim, double &th_sim, double &ph_sim, double &th_in, double &ph_in, double &x0, double &y0, double &z0, bool &gutProp, double &rx0, double &rd0, double xacc=25, double yacc=0, double racc1=3, double racc2=9, double thacc=0.038, double phacc=0., double thlim = 0.008, double phlim=0.210,double zacc=1110.){
  /// Read results from forward propagation from tree ------------------------------
  TTree *tByProp_in = (TTree*)f_in->Get(treeName);
  TBranch *br_x_in = tByProp_in->GetBranch("x_out");
  TBranch *br_y_in = tByProp_in->GetBranch("y_out");
  TBranch *br_z_in = tByProp_in->GetBranch("z_out");
  TBranch *br_th_in = tByProp_in->GetBranch("th_out");
  TBranch *br_ph_in = tByProp_in->GetBranch("ph_out");
  TBranch *br_x_sim = tByProp_in->GetBranch("x_in");
  TBranch *br_y_sim = tByProp_in->GetBranch("y_in");
  TBranch *br_z_sim = tByProp_in->GetBranch("z_in");
  TBranch *br_th_sim = tByProp_in->GetBranch("th_in");
  TBranch *br_ph_sim = tByProp_in->GetBranch("ph_in");
  double x_in_b,y_in_b,z_in_b,th_in_b,ph_in_b; //input for backward propagation
  br_x_in->SetAddress(&x_in_b);
  br_y_in->SetAddress(&y_in_b);
  br_z_in->SetAddress(&z_in_b);
  br_th_in->SetAddress(&th_in_b);
  br_ph_in->SetAddress(&ph_in_b);
  //input for simulation
  br_x_sim->SetAddress(&x_sim);
  br_y_sim->SetAddress(&y_sim);
  br_z_sim->SetAddress(&z_sim);
  br_th_sim->SetAddress(&th_sim);
  br_ph_sim->SetAddress(&ph_sim);

  int nevent = tByProp_in->GetEntries();
  if(itree>nevent){
    cout<<"Tree has "<<nevent<<" entries! Entry #"<<itree<<" doesn't exist!"<<endl;
    return 1;
  }
  //  for(int itree=0;itree<nevent;itree++){
    tByProp_in->GetEvent(itree);
    x0=x_in_b; 
    y0=y_in_b;
    z0=z_in_b;
    th_in = th_in_b;
    ph_in =ph_in_b;
    // if(ph_in_b>0) ph_in =ph_in_b-TMath::Pi();
    // if(ph_in_b<0) ph_in =ph_in_b+TMath::Pi();
    // cout<<"ph_in_b = "<<ph_in_b<<" ph_in = "<<ph_in<<endl;
    ph_in =ph_in_b;
    //check if params within acceptance
    bool ifProp=false;
    if(z_in_b>zacc){
      ifProp=true;
      // cout<<"Trk was forward extraploted!"<<endl;
    }
    
    bool ifPropxy = true;    
    bool ifPropthph = true;
    double r = sqrt((x0-xacc)*(x0-xacc)+(y0-yacc)*(y0-yacc));
    rx0 = r;
    if(xycut){
      if(r<racc1 || r>racc2) ifPropxy=false;
    }
    //   if(!ifPropxy && ifProp) cout<<" x0 = "<<x0<<" y0 = "<<y0<<endl;
    double ph_cur = ph_in_b;
    // if(ph_in_b>=0) ph_cur =ph_in_b-TMath::Pi();
    // if(ph_in_b<0) ph_cur =ph_in_b+TMath::Pi();

    // double a = ph_cur-phacc;
    // double b = th_in-thacc;
    // TVector2 rphth(a,b);//radius-vector in ellipse
    // double phi_rphth = rphth.Phi();
    // double rd = a*b/sqrt(TMath::Power(b*TMath::Cos(phi_rphth),2)+TMath::Power(a*TMath::Sin(phi_rphth),2));
    // rd0 = rd;
    // if(thphcut){
    // 	if(fabs(rd)<rdacc1 || fabs(rd)>rdacc2) ifPropthph=false;
    // }

    if(thphcut){
      //   cout<<"ph_cur-phacc = "<<ph_cur-phacc<<" th_in-thacc = "<<th_in-thacc<<endl;
      if(fabs(th_in-thacc)>thlim || fabs(ph_cur-phacc)>phlim) ifPropthph=false;
    }

    //  if(!ifPropthph) cout<<"rd0"<<rd0<<endl;
    if(ifPropthph && ifPropxy && ifProp)
      gutProp=true;
    else gutProp=false;
  /// [END] Read results from forward propagation from tree ---------------------
 return 0;
}
     

//void RK_By(double momentum=3.0, int charge=-1)
void RK_By(int nEv=1000, double Pbeam=15, int imom=0, int icharge=-1, int th1=3, int th2=8, bool fl_in_prop=true, bool forwAss=true, bool xycut=true, bool thphcut=true, int seed=0, TString path="/home/karavdin/pandaRoot12/karavdina/RungeKuttaSimplified/")
//int main(double momentum=3.0, int charge=-1)
{
#if defined(__CINT__) && !defined(__MAKECINT__) 
  cout << "WARNING: This macro can run only using ACliC, you must run it by doing: " << endl;
  cout << "\t .x RK_By.C(double momentum, int dmomentum%, int charge, isForward)+" << endl; 
  //  return -1;
  return;
#endif

    // Draw_Field();
    // return;
   // bool fl_in_prop = true;// if TRUE start from the beggining, usually forward, FALSE = read results of previous propagation
  //  bool fl_in_prop = false;// if TRUE start from the beggining, usually forward, FALSE = read results of previous propagation
  TRandom3 rand(seed);

  //read mag.field and integrals
  int okfl = FieldTable("15_GeV_FieldIntegTable.txt");
  if(okfl<1) cout<<"Read lookup table for mag.field!"<<endl;

  //  return;

  // read data from forward propagation -----------------------------
  
  TString ftoysimname = path+"/forward/Prop_By_th_";
  ftoysimname+=th1;
  ftoysimname+="_";
  ftoysimname+=th2;
  ftoysimname+="_mrad_";
  ftoysimname +=Pbeam;
  ftoysimname +="_PLUS_";
  ftoysimname+=imom;
  ftoysimname+="procent_q_";
  ftoysimname+=icharge;
 // ftoysimname+="_isFfB_";
 // ftoysimname+=forwAss;
 //  ftoysimname+="_isXYcut_";
 //  ftoysimname+=xycut;
 //  ftoysimname+="_isTHPHcut_";
 //  ftoysimname+=thphcut;
  ftoysimname+=".root";

  TFile *f_in;
      if(fl_in_prop) f_in = new TFile(ftoysimname,"RECREATE");
      else  f_in = new TFile(ftoysimname,"READ");
      double   momentum = Pbeam+Pbeam*imom/100.;
      if(fl_in_prop) cout<<"We create data for: momentum = "<<momentum<<" charge = "<<icharge<<endl;
      else cout<<"We read data for: momentum = "<<momentum<<" charge = "<<icharge<<endl;


      TString momName = path+"/backward/Prop_By_th_";
      momName+=th1;
      momName+="_";
      momName+=th2;
      momName+="_mrad_";
      momName +=Pbeam;
      momName +="_PLUS_";
      momName+=imom;
      momName+="procent_q_";
      momName+=icharge;
      momName+="_isFfB_";
      momName+=forwAss;
      momName+="_isXYcut_";
      momName+=xycut;
      momName+="_isTHPHcut_";
      momName+=thphcut;
      momName+=".root";


      TFile *f_out;
      if(!fl_in_prop) f_out= new TFile(momName,"RECREATE");

    
      // const double mom=momentum;//GeV //TEST
      // const double q = icharge; //barp=-1!!! //TEST

      // TH2D *hthphi_in = new TH2D("hthphi_in_b","initial params;#phi, rad;#theta, mrad",1e3,-TMath::Pi(),TMath::Pi(),1e3,0,100);
      // TH2D *hthphi_out = new TH2D("hthphi_out_b","final params ;#phi, rad;#theta, mrad",1e3,-TMath::Pi(),TMath::Pi(),1e3,0,100);
      // TH2D *hxy_out = new TH2D("hxy_out_b","final params ;x, cm;y, cm",1e3,-10,10,1e3,-10,10);

      TString treeName("tByProp");
      TString treeNameout = treeName+"out";
     
      //      const int nsteps=1.1e5; //number of steps in propagation
    
      //  double z0=0;
     
      TTree *tByProp;
      if(fl_in_prop) tByProp = new TTree(treeName, "trk params before and after propagation in By [forward]");
      else
	tByProp = new TTree(treeNameout, "trk params before and after propagation in By [backward]");
      Double_t x_in,y_in,z_in,th_in,ph_in;
      Double_t x_out,y_out,z_out,th_out,ph_out;
      tByProp->Branch("x_in", &x_in);
      tByProp->Branch("y_in", &y_in);
      tByProp->Branch("z_in", &z_in);
      tByProp->Branch("th_in", &th_in);
      tByProp->Branch("ph_in", &ph_in);
      tByProp->Branch("x_out", &x_out);
      tByProp->Branch("y_out", &y_out);
      tByProp->Branch("z_out", &z_out);
      tByProp->Branch("th_out", &th_out);
      tByProp->Branch("ph_out", &ph_out);
      Double_t x_sim,y_sim,z_sim,th_sim,ph_sim;
      Double_t rd_in, r_in;
      if(!fl_in_prop){
	tByProp->Branch("x_sim", &x_sim);
	tByProp->Branch("y_sim", &y_sim);
	tByProp->Branch("z_sim", &z_sim);
	tByProp->Branch("th_sim", &th_sim);
	tByProp->Branch("ph_sim", &ph_sim);
	tByProp->Branch("rd_in", &rd_in);
	tByProp->Branch("r_in", &r_in);
      }


      const int nsampst = nEv;//number of th&ph samples 
      //      const int nsampst = 1e3;//number of th&ph samples 
      if(fl_in_prop){ //forward extrapolation
	const double factor = TMath::Power((fabs(icharge)/momentum),0);
	//	const int nsteps=factor*1e8;
	//	const int nsteps=factor*1e4;
       
	const int nsteps=factor*5e3;
	double zst = 1124./nsteps;
	//	cout<<"factor for number of internal steps = "<<factor<<" zst = "<<zst<<endl;
	/// Fill tree with results from forward propagation from ---------------------------
	//	double thmin=0; double thmax=0.25*TMath::Pi(); 
	//	double thmin=2e-3; double thmax=9e-3; 
	//	double thmin=3e-3; double thmax=8e-3; 
	double thmin=th1*1e-3; double thmax=th2*1e-3; 
	double phmin=0; double phmax=2*TMath::Pi();
	double th0; double ph0; double x0; double y0; double z0;
	for(int ithph=0;ithph<nsampst;ithph++){// loop over diff theta&phi, which are sampled uniformly 
	  if((ithph%200)==0) cout<<"th&ph st#"<<ithph<<endl;
	  Generate_pars_IP(&rand,th0,ph0,x0,y0,z0,thmin,thmax,phmin,phmax);
	  th_in = th0;  ph_in = ph0;
	  TVector3 dir(0,0,1);
	  dir.SetTheta(th0);
	  dir.SetPhi(ph0);
	  double dirXYZ[2];
	  dir.GetXYZ(dirXYZ);
	  double tx0 = dirXYZ[0];
	  double ty0 = dirXYZ[1];
	  double tz0 = dirXYZ[2];
	  x_in = x0;  y_in = y0;
	  x_out = x0;  y_out = y0;

	  //propagate ----------------------------
	    // double gBy[nsteps]; double zBy[nsteps]; //TEST By field
	    // gBy[0] = 0; zBy[0] = 0;
	  //	  int acst = 0;
	  for(int i=1;i<nsteps+1;i++){//forward

	    //	for(int i=nsteps;i>-1;i--){//backfward
	    double z1 = z0+(i-1)*zst;
	    double z2 = z0+i*zst;
	    double zpair[2] = {z1,z2};
	    // gBy[i] = By(zpair);
	    // zBy[i] = z1;
	    //   acst++;
	    // if(By(zpair)!=0){
            // cout.precision(15);
	    //   cout<<z1<<" "<<z2<<" "<<By(zpair)<<" "<<By_integ_1(zpair)<<" "<<By_integ_2(zpair)<<" "<<By_integ_3(zpair)<<endl;
	    // }


	    //   cout<<" mom = "<<mom<<" q = "<<q<<endl;
	    // double tx0_pr = tx0;
	    // double ty0_pr = ty0;

	    double x1=x_out; double y1=y_out; double tx1=tx0; double ty1=ty0;
	    double x2; double y2; double tx2; double ty2;
	    int okextr = extrapol(momentum,icharge,x1,y1,tx1,ty1,z1,z2,x2,y2,tx2,ty2);
            x_out=x2; y_out = y2; tx0 = tx2; ty0 = ty2;
	    z_out = z2;
	    // tx0 =dtx0; 	    ty0 =dty0;
	    // x_out += dx0;  	    y_out += dy0; 
	    // z_out = z2;

	    // // tx0 = tx_sum(momentum,icharge,tx0,ty0,z1,z2);
	    // // ty0 = ty_sum(momentum,icharge,tx0,ty0,z1,z2);
	    // // //   cout<<"tx0 = "<<tx0<<" ty0 = "<<ty0<<endl;
	    // double x_out_exact = x1+dx_sum(momentum,icharge,tx1,ty1,z1,z2);
	    // if(fabs(x_out-x_out_exact)>1e-1){
	    // cout<<"x0 = "<<x_out<<" z_out = "<<z_out<<endl;
	    // cout<<"x_out_exact = "<<x_out_exact<<endl;
	    // break;
	    // }

	    // y_out += dy_sum(momentum,icharge,tx0_pr,ty0_pr,z1,z2);
	    // z_out = z1;
	    if(fabs(x_out)>5e1 || fabs(y_out)>5e1) break;
	    // ///TEST!!!!!! ---------------
	    // if((tx0*tx0+ty0*ty0)<1) tz0 = TMath::Sqrt(1-tx0*tx0-ty0*ty0);
	    // else tz0 = -999;
	    // //	    if(tz0!=-999){
	    //   TVector3 dirFin(tx0,ty0,tz0);
	    //   th_out = dirFin.Theta(); 
	    //   //ph_out = dirFin.Phi();
	    //   if(dirFin.Phi()>0) ph_out = dirFin.Phi()-TMath::Pi();
	    //   if(dirFin.Phi()<0) ph_out = dirFin.Phi()+TMath::Pi();
	    //   //	    }
	    //   tByProp->Fill();
	    //   ///[END] TEST!!!!!! ---------------
	  }
	  // TGraph *gr = new TGraph(acst,zBy,gBy);
	  // gr->Draw("A*");
	  //  [END]propagate ----------------------------
	  if((tx0*tx0+ty0*ty0)<1) tz0 = TMath::Sqrt(1-tx0*tx0-ty0*ty0);
	  else tz0 = -999;
	  if(tz0!=-999){
	    TVector3 dirFin(tx0,ty0,tz0);
	    th_out = dirFin.Theta(); 
  	    ph_out = dirFin.Phi();
	    // if(dirFin.Phi()>0) ph_out = dirFin.Phi()-TMath::Pi();
	    // if(dirFin.Phi()<0) ph_out = dirFin.Phi()+TMath::Pi();
	    tByProp->Fill();
	  }
	}
	/// [END] Fill tree with results from forward propagation from ------------------
	f_in->Write();
      }
      else{
	double mom; double q;
	if(forwAss){
	  //use true values
	  mom=momentum;//GeV ,
	  q = icharge; //[e]
	}
	else{	  //barp assumption for propagation
	  mom=Pbeam;//GeV
	  //	  mom=1.5;//GeV
	  q = -1; //barp=-1!!!
	}
	
	const double factor = TMath::Power((fabs(q)/mom),0);
	const int nsteps=factor*5e3;
	//	const int nsteps=factor*5e3;
	double zst = 1124./nsteps;
	cout<<"assumption for propagation: mom="<<mom<<" q="<<q<<", factor for number of internal steps = "<<factor<<endl;
	/// Read results from forward propagation from tree ------------------------------
	double th0; double ph0; double x0; double y0; double z0; bool gutProp;
	//	double xacc1=15.5; double xacc2=33.8; double yacc=9; double thacc1=0.0295; double thacc2=0.046; double phacc=0.22; double zacc=1110.;
	double xacc=25; double yacc=0; double racc1=3.1; double racc2=9.3; double thacc=0.038; double phacc=0.; double thlim=0.011; double phlim=0.25; double zacc=1120.;
	for(int isamp=0;isamp<nsampst;isamp++){//go over diff. started th&ph
	//	for(int isamp=nsampst-1;isamp>-1;isamp--){//invert go over diff. started th&ph 
	//	for(int isamp=nsteps-1;isamp>nsteps-2;isamp--){//TEST
	  //	  cout<<"we read event from tree #"<<isamp<<endl;
	  int resread = Read_propag_pars(xycut,thphcut,treeName,f_in,isamp,x_sim,y_sim,z_sim,th_sim,ph_sim,th0,ph0,x0,y0,z0,gutProp,r_in,rd_in,xacc,yacc,racc1,racc2,thacc,phacc,thlim,phlim,zacc);


	  if(!gutProp) continue;
	  if(resread>0) break;
	  // cout<<"resread = "<<resread<<endl;
	  //!!!TEST: write everything!!! done :)

	  //	  double x_out_tmp, y_out_tmp, z_out_tmp;
	  th_in = th0;  ph_in = ph0;
	  TVector3 dir(0,0,1);
	  dir.SetTheta(th_in);
	  dir.SetPhi(ph_in);
	  double dirXYZ[2];
	  dir.GetXYZ(dirXYZ);
	  double tx0 = dirXYZ[0];
	  double ty0 = dirXYZ[1];
	  double tz0 = dirXYZ[2];
	  x_in = x0;  y_in = y0;
	  z_in = z0;
	  //	  x_out_tmp = x0;  y_out_tmp = y0;
	
	  //propagate ----------------------------
	  // for(int i=1;i<nsteps;i++){
	  //   double z1 = z_in-(i+1)*zst;//backfward
	  //   double z2 = z_in-i*zst;
	  for(int i=0;i<nsteps;i++){ //TEST !!!
	    double z1 = z_in-(i+1)*zst;//backfward
	    double z2 = z_in-i*zst;
	
	    // // cout<<"z2-z1 = "<<z2-z1<<endl;
	    // // double tx0_pr = tx0;
	    // // double ty0_pr = ty0;

	    // tx0 = tx_sum(mom,q,tx0,ty0,z2,z1);
	    // ty0 = ty_sum(mom,q,tx0,ty0,z2,z1);
	   

	    // x_out_tmp += dx_sum(mom,q,tx0_pr,ty0_pr,z2,z1);
	    // y_out_tmp += dy_sum(mom,q,tx0_pr,ty0_pr,z2,z1);

	    // z_out_tmp = z2;

	    double x2 = x0; double y2 = y0; double tx2 = tx0; double ty2 = ty0;
	    double x1; double y1; double tx1; double ty1;
	    int okextr = extrapol(momentum,icharge,x2,y2,tx2,ty2,z2,z1,x1,y1,tx1,ty1);
	    x0 = x1; y0 = y1; tx0 = tx1; ty0 = ty1;
	    //	    tx0 =dtx0; 	    ty0 =dty0;
	    x_out = x0;  	    y_out = y0; 
	    z_out = z1;
	    //   cout<<"x_out = "<<x_out<<" y_out = "<<y_out<<" z_out = "<<z_out<<endl;

	    //	    if((i%1000)==0)  cout<<"OUT x,y,z = "<<x_out_tmp<<", "<<y_out_tmp<<", "<<z_out_tmp<<endl;
	    //  cout<<"OUT tx0,ty0 = "<<tx0<<", "<<ty0<<" @z=["<<z1<<" -- "<<z2<<"]"<<endl;
	    //    if(fabs(x_out)>1e3 || fabs(y_out)>1e3) break;
	    // ///TEST!!!!!! ---------------
	    // if((tx0*tx0+ty0*ty0)<1) tz0 = TMath::Sqrt(1-tx0*tx0-ty0*ty0);
	    // else tz0 = -999;
	    // if(tz0!=-999){
	    //   TVector3 dirFin(tx0,ty0,tz0);
	    //   th_out = dirFin.Theta(); 
	    //   //ph_out = dirFin.Phi();
	    //   if(dirFin.Phi()>0) ph_out = dirFin.Phi()-TMath::Pi();
	    //   if(dirFin.Phi()<0) ph_out = dirFin.Phi()+TMath::Pi();
	    //   // x_out = x_out_tmp;
	    //   // y_out = y_out_tmp;
	    //   // z_out = z_out_tmp;
	    //   tByProp->Fill();
	    // }
	    //   ///[END] TEST!!!!!! ---------------
	  }
	  //[END]propagate ----------------------------
	  if((tx0*tx0+ty0*ty0)<1) tz0 = TMath::Sqrt(1-tx0*tx0-ty0*ty0);
	  else tz0 = -999;
	  if(tz0!=-999){
	    TVector3 dirFin(tx0,ty0,tz0);
	    th_out = dirFin.Theta(); 
  	    //ph_out = dirFin.Phi();
	    // if(dirFin.Phi()>0) ph_out = dirFin.Phi()-TMath::Pi();
	    // if(dirFin.Phi()<0) ph_out = dirFin.Phi()+TMath::Pi();
	    ph_out = dirFin.Phi();
	    // x_out = x_out_tmp;
	    // y_out = y_out_tmp;
	    // z_out = z_out_tmp;
	    tByProp->Fill();
	    //	    if((isamp%1000)==0)  
	    cout<<"For sample#"<<isamp<<" OUT x,y,z = "<<x_out<<", "<<y_out<<", "<<z_out<<endl;
	  }
	}
	/// [END] Read results from forward propagation from tree ---------------------
	f_out->Write();
      }
      //    }
    //}
}