// 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 <TH3.h>
#include <TF1.h>
#include <TF2.h>
#include <TF3.h>
#include <TString.h>
#include "TObjString.h"
#include "TRandom3.h"
#include "TObjArray.h"
#include "TFrame.h"


using namespace std;

#ifdef __MAKECINT__
#pragma link C++ class vector<float>+;
#endif





void run(){
    
    
   TFile *f = TFile::Open("Results.root","READ");
   TTree *t= (TTree*)f->Get("tvec");
   int nEntries=t->GetEntries();

  
  vector<vector<double> > history_ps_x_all_part;
  vector<vector<double> > history_ps_y_all_part;
  vector<vector<double> > history_ps_z_all_part;
  vector<vector<double> > history_vl_x_all_part;
  vector<vector<double> > history_vl_y_all_part;
  vector<vector<double> > history_vl_z_all_part;
  vector<vector<double> > history_ac_x_all_part;
  vector<vector<double> > history_ac_y_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_ps_z_bunch_all_part;  
   
  vector<double> *x           =0;
  vector<double> *y           =0;
  vector<double> *z           =0;
  vector<double> *vx          =0;
  vector<double> *vy          =0;
  vector<double> *vz          =0;
  vector<double> *ax          =0;
  vector<double> *ay          =0;
  vector<double> *az          =0;
  vector<double> *time        =0;   
  vector<double> *Field_x     =0;   
  vector<double> *Field_y     =0;   
  vector<double> *Field_z     =0;   
  vector<double> *Field_y_ext =0;   
  vector<double> *z_bunch     =0; 
  
  
  TBranch *bx = 0;
  TBranch *by = 0;
  TBranch *bz = 0;
  TBranch *bvx = 0;
  TBranch *bvy = 0;
  TBranch *bvz = 0;
  TBranch *bax = 0;
  TBranch *bay = 0;
  TBranch *baz = 0;
  TBranch *btime        = 0;
  TBranch *bField_x     = 0;
  TBranch *bField_y     = 0;
  TBranch *bField_z     = 0;
  TBranch *bField_y_ext = 0;
  TBranch *bz_bunch     = 0;  
  
  t->SetBranchAddress("x",&x,&bx);
  t->SetBranchAddress("y",&y,&by);
  t->SetBranchAddress("z",&z,&bz);
  t->SetBranchAddress("vx",&vx,&bvx);
  t->SetBranchAddress("vy",&vy,&bvy);
  t->SetBranchAddress("vz",&vz,&bvz);
  t->SetBranchAddress("ax",&ax,&bax);
  t->SetBranchAddress("ay",&ay,&bay);
  t->SetBranchAddress("az",&az,&baz);
  t->SetBranchAddress("time",&time,&btime);  
  t->SetBranchAddress("Field_x",&Field_x,&bField_x);
  t->SetBranchAddress("Field_y",&Field_y,&bField_y); 
  t->SetBranchAddress("Field_z",&Field_z,&bField_z);
  t->SetBranchAddress("Field_y_ext",&Field_y_ext,&bField_y_ext); 
  t->SetBranchAddress("z_bunch",&z_bunch,&bz_bunch);   


   for (Int_t i = 0; i < nEntries; i++) {
      
      Long64_t tentry = t->LoadTree(i);
      bx->GetEntry(tentry);
      history_ps_x_all_part.push_back(*x);
      by->GetEntry(tentry);
      history_ps_y_all_part.push_back(*y);
      bz->GetEntry(tentry);
      history_ps_z_all_part.push_back(*z);
      bvx->GetEntry(tentry);
      history_vl_x_all_part.push_back(*vx);
      bvy->GetEntry(tentry);
      history_vl_y_all_part.push_back(*vy);
      bvz->GetEntry(tentry);
      history_vl_z_all_part.push_back(*vz);

       
      bax->GetEntry(tentry);
      history_ac_x_all_part.push_back(*ax);
      bay->GetEntry(tentry);
      history_ac_y_all_part.push_back(*ay);
      baz->GetEntry(tentry);
      history_ac_z_all_part.push_back(*az);
      btime->GetEntry(tentry);
      history_tm_all_part.push_back(*time);    
      bField_x->GetEntry(tentry);
      history_Field_x_all_part.push_back(*Field_x);    
      bField_y->GetEntry(tentry);
      history_Field_y_all_part.push_back(*Field_y);    
      bField_z->GetEntry(tentry);
      history_Field_z_all_part.push_back(*Field_z);    
      bField_y_ext->GetEntry(tentry);
      history_Field_y_ext_all_part.push_back(*Field_y_ext);    
      bz_bunch->GetEntry(tentry);
      history_ps_z_bunch_all_part.push_back(*z_bunch);   

}


   
   // Since we passed the address of a local variable we need
   // to remove it.
   t->ResetBranchAddresses();  // It deletes bx, by, ...
   delete f; // automatically deletes "t", too
   delete x;
   delete y;
   delete z;
   delete vx;
   delete vy;
   delete vz;   
   delete ax;
   delete ay;
   delete az;
   delete time;   
   delete Field_x;   
   delete Field_y;   
   delete Field_z;   
   delete Field_y_ext;   
   delete z_bunch;

   double sigma_x  = 3.8e-03;      
   double distance = 5e-02;
   
/*   TFile*fout = new TFile("3DPlots.root","recreate");
   TObjArray HList(0);
   for(int i=0;i<history_ps_y_all_part.size();i++){
    for(int j=0;j<history_ps_y_all_part[i].size();j++){
      // x vs y vs y
      TString histoname =  TString::Format("x_vs_y__vs_z_time_indeks_%d", j);
      TString histotitle = TString::Format("x vs y vs z,  time = %e ns",history_tm_all_part[i][j]*1.e+9); 
      TH3D *myhist = ((TH3D *)(HList.FindObject(histoname)));
      if (!myhist)
        {
          myhist = new TH3D(histoname,
                            "x vs y vs z",
                            300, -8*sigma_x*1000, 8*sigma_x*1000, 300,-distance*1000.-50,distance*1000.+50, 300,-10,+10 );
          myhist->GetXaxis()->SetTitle("x (mm)");
          myhist->GetYaxis()->SetTitle("y (mm)");
          myhist->GetZaxis()->SetTitle("z (mm)");
          myhist->SetTitle(histotitle);
          myhist->SetMarkerStyle(20);
          myhist->SetMarkerColor(2);
          myhist->SetMarkerSize(0.4);
          
          HList.Add(myhist);
        } 
       myhist->Fill(history_ps_x_all_part[i][j]*1000.,history_ps_y_all_part[i][j]*1000.,history_ps_z_all_part[i][j]*1000.); 
        
       } // for (loop over events, elements, tracks, hits) ...
   }
 HList.Write();
 fout->Write();    
 
 */
        
  TFile*f_test= new TFile("3DPlots.root","recreate");
  TObjArray HList(0);
  for(int i=0;i<history_ps_y_all_part.size();i++){
    for(int j=0;j<history_ps_y_all_part[i].size();j++){
      // x vs y 
      TString histoname =  TString::Format("x_vs_yz_time_indeks_%d", j);
      TString histotitle = TString::Format("x vs yz,  time = %e ns",history_tm_all_part[i][j]*1.e+9); 
      TH3D *myhist = ((TH3D *)(HList.FindObject(histoname)));
      if (!myhist)
        {
          myhist = new TH3D(histoname,
                            "x vs y",
                            300, -8*sigma_x*1000, 8*sigma_x*1000, 300,-distance*1000.-50,distance*1000.+50,300, -8*sigma_x*1000, 8*sigma_x*1000);
          myhist->GetXaxis()->SetTitle("x (mm)");
          myhist->GetYaxis()->SetTitle("y (mm)");
          myhist->GetZaxis()->SetTitle("z (mm)");
          myhist->SetTitle(histotitle);
          myhist->SetMarkerStyle(20);
          myhist->SetMarkerColor(2);
          myhist->SetMarkerSize(0.4);
          
          HList.Add(myhist);
        } 
     
      myhist->Fill(history_ps_x_all_part[i][j]*1000.,history_ps_y_all_part[i][j]*1000.,history_ps_z_all_part[i][j]*1000.);
      // vx vs x
      TString histoname2 =  TString::Format("speedx_vs_x_time_indeks_%d", j);
      TString histotitle2 = TString::Format("vx vs x,  time = %e ns",history_tm_all_part[i][j]*1.e+9); 
      TH2D *myhist2 = ((TH2D *)(HList.FindObject(histoname2)));
       if (!myhist2)
         {
           myhist2 = new TH2D(histoname2,
                             "vx vs v",
                             300, -8*sigma_x*1000, 8*sigma_x*1000, 300,-1e8,1e8);
           myhist2->GetXaxis()->SetTitle("x (mm)");
           myhist2->GetYaxis()->SetTitle("v_x (m/s)");
           myhist2->SetTitle(histotitle2);
           myhist2->SetMarkerStyle(20);
           myhist2->SetMarkerColor(2);
           myhist2->SetMarkerSize(0.4);
           
           HList.Add(myhist2);
         } 
      
       myhist2->Fill(history_ps_x_all_part[i][j]*1000.,history_vl_x_all_part[i][j]);        
      // vy vs y
      TString histoname3 =  TString::Format("speedy_vs_y_time_indeks_%d", j);
      TString histotitle3 = TString::Format("vy vs y,  time = %e ns",history_tm_all_part[i][j]*1.e+9); 
      TH2D *myhist3 = ((TH2D *)(HList.FindObject(histoname3)));
       if (!myhist3)
         {
           myhist3 = new TH2D(histoname3,
                             "vy vs y",
                             300, -10., distance*1000.+50, 300,-1e5,3e8);
           myhist3->GetXaxis()->SetTitle("y (mm)");
           myhist3->GetYaxis()->SetTitle("v_y (m/s)");
           myhist3->SetTitle(histotitle3);
           myhist3->SetMarkerStyle(20);
           myhist3->SetMarkerColor(2);
           myhist3->SetMarkerSize(0.4);
           
           HList.Add(myhist3);
         } 
      
       myhist3->Fill(history_ps_y_all_part[i][j]*1000.,history_vl_y_all_part[i][j]);        
         // vz vs z
      TString histoname4 =  TString::Format("speedz_vs_z_time_indeks_%d", j);
      TString histotitle4 = TString::Format("vz vs z,  time = %e ns",history_tm_all_part[i][j]*1.e+9); 
      TH2D *myhist4 = ((TH2D *)(HList.FindObject(histoname4)));
       if (!myhist4)
         {
           myhist4 = new TH2D(histoname4,
                             "vz vs z",
                             300, -10.,10., 300,-1e5,1e5);
           myhist4->GetXaxis()->SetTitle("z (mm)");
           myhist4->GetYaxis()->SetTitle("v_z (m/s)");
           myhist4->SetTitle(histotitle4);
           myhist4->SetMarkerStyle(20);
           myhist4->SetMarkerColor(2);
           myhist4->SetMarkerSize(0.4);
           
           HList.Add(myhist4);
         } 
      
       myhist4->Fill(history_ps_z_all_part[i][j]*1000.,history_vl_z_all_part[i][j]);        
        // x counts
      TString histoname5 =  TString::Format("x_counts_time_indeks_%d", j);
      TString histotitle5 = TString::Format("x distribution,  time = %e ns",history_tm_all_part[i][j]*1.e+9); 
      TH1D *myhist5 = ((TH1D *)(HList.FindObject(histoname5)));
       if (!myhist5)
         {
           myhist5 = new TH1D(histoname5,
                             "x distribution",
                             300, -8*sigma_x*1000, 8*sigma_x*1000);
           myhist5->GetXaxis()->SetTitle("x (mm)");
           myhist5->GetYaxis()->SetTitle("Counts");
           myhist5->SetTitle(histotitle5);
           myhist5->SetMarkerStyle(20);
           myhist5->SetMarkerColor(2);
           myhist5->SetMarkerSize(0.4);
           
           HList.Add(myhist5);
         } 
      
       myhist5->Fill(history_ps_x_all_part[i][j]*1000.);        
        // y counts
      TString histoname6 =  TString::Format("y_counts_time_indeks_%d", j);
      TString histotitle6 = TString::Format("y distribution,  time = %e ns",history_tm_all_part[i][j]*1.e+9); 
      TH1D *myhist6 = ((TH1D *)(HList.FindObject(histoname6)));
       if (!myhist6)
         {
           myhist6 = new TH1D(histoname6,
                             "y distribution",
                             300.,-10., distance*1000+50);
           myhist6->GetXaxis()->SetTitle("y (mm)");
           myhist6->GetYaxis()->SetTitle("Counts");
           myhist6->SetTitle(histotitle6);
           myhist6->SetMarkerStyle(20);
           myhist6->SetMarkerColor(2);
           myhist6->SetMarkerSize(0.4);
           
           HList.Add(myhist6);
         } 
      
       myhist6->Fill(history_ps_y_all_part[i][j]*1000.);          
     // Ratio of electric field in the lab system x to y components vs x position.
      TString histoname7 =  TString::Format("Fieldx_over_Fieldy_time_indeks_%d", j);
      TString histotitle7 = TString::Format("Field_x/Field_y in lab ref frame vs x,  time = %e ns",history_tm_all_part[i][j]*1.e+9); 
      TH2D *myhist7 = ((TH2D *)(HList.FindObject(histoname7)));
       if (!myhist7)
         {
           myhist7 = new TH2D(histoname7,
                             "Field ratios",
                              300, -8*sigma_x*1000, 8*sigma_x*1000, 300, 1e-8, 1.);
           myhist7->GetXaxis()->SetTitle("x (mm)");
           myhist7->GetYaxis()->SetTitle("Ex/Ey");
           myhist7->SetTitle(histotitle7);
           myhist7->SetMarkerStyle(20);
           myhist7->SetMarkerColor(2);
           myhist7->SetMarkerSize(0.4);
           
           HList.Add(myhist7);
         } 
      
       myhist7->Fill(history_ps_x_all_part[i][j]*1000. ,abs(history_Field_x_all_part[i][j]/(history_Field_y_all_part[i][j]+history_Field_y_ext_all_part[i][j]) )); 
       
     //  cout <<i << " " << j << " " <<history_ps_x_all_part[i][j]*1000. <<" " << history_Field_x_all_part[i][j]/(history_Field_y_all_part[i][j]+history_Field_y_ext_all_part[i][j]) <<endl;
    } // for (loop over events, elements, tracks, hits) ...
   }
 HList.Write();
 f_test->Write();    
    
       
        
        
        
//         for (Int_t i = 0; i < 2; i++) {
//           for (UInt_t j = 0; j <history_ps_x_all_part[i].size() ; ++j) {
//                  cout << i << " " <<history_ps_x_all_part[i][j]<< " " <<   history_ps_y_all_part[i][j]<<" " << history_ps_z_all_part[i][j]<<endl;
//                  //cout << history_vl_x_all_part[i][j]<< " " <<   history_vl_y_all_part[i][j]<<" " << history_vl_z_all_part[i][j]<<endl;
//                  //cout << history_ac_x_all_part[i][j]<< " " <<   history_ac_y_all_part[i][j]<<" " << history_ac_z_all_part[i][j]<<endl;
//                  
//           }}

   
   
}