#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <string>
#include <string.h>

#include "TTree.h"
#include "TFile.h"
#include "TString.h"

#include "TSystem.h"
#include "THaEvent.h"
#include "THaRun.h"
#include "TH1.h"
#include "TF1.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TChain.h"
#include "TEllipse.h"
#include "TLine.h"
#include "TMarker.h"
#include "TVector3.h"
#include "TRotation.h"
#include "TMath.h"
#include "TGraph.h"


#include "TMath.h"
#include "TComplex.h"


#include <cstring>
#include <cstdio>
#include <iostream>

//******************* E-plane corrections ****************************************

  double EGainCorrections[24][4] = {	
{1.01282, 2.12441, 2.51691, 1067.13},
{0.998535, 2.531, 3.01, 1569.59},
{1.00527, 1.735, 2.04935, 1727.12},
{0.987726, 2.08, 1.84, 1856.0},
{1.00142, 1.72415, 1.38549, 1919.02}, 
{0.985367, 1.98209, 1.47693, 1847.42},//{0.985367, 2.896, 2.72, 746.4},
{0.980072, 1.73345, 1.62718, 1917.7}, //{0.980072, 2.896, 2.72, 777.9}, //6
{0.993866, 1.40852, 1.7899, 2020.53},
{1.0119, 1.93019, 1.90826, 1964.34}, //{1.0119, 2.896, 2.72, 777.9},
{0.979065, 1.75472, 1.90445, 1824.34}, //9
{0.987534, 1.72314, 1.77495, 1967.84}, //{0.987534, 2.896, 2.72, 777.9},
{1.01997, 1.85261, 1.66637, 1906.34}, //11
{0.99292, 1.92804, 1.55635, 1906.34}, //12
{0.983154, 1.84102, 1.45781, 1988.34}, //13
{0.998004, 1.55193, 1.26637, 1988.34}, //14
{0.996338, 2.394, 2.495, 807.6}, //{0.996338, 2.15438, 1.59494, 1967.84}, //15  // now aligned to MIP
{0.998474, 2.095, 2.58, 804.7}, //{0.998474, 1.79016, 1.49752, 2012.59}, //16 // aligned to MIP
{1.0009,   1.22551, 1.38512, 1992.19}, //17
{0.996709, 0.2995, 0.8678, 785.1}, //{0.996709, 1.31375, 1.50701, 1864.87}, //18
{0.97673, 0.7594, 0.5874, 734.7}, //19
{1.01855, 1.881, 3.014, 670.5}, //{1.01855, 1.71273, 1.82481, 1706.95}, //20
{0.996763, 1.312, 2.45, 830.1}, //21
{1.00377,  0.641, 1.133, 909.2}, //22
{0.842193, 0.0, 0.0, 0.0} //23
};


// upper bounary must be 2200. Lower boundary must be 340
double EFinalCorrections[24] = {	
0.91892,//1.0,
0.952177,//1.01939,
1.1458333, //1.05,//1.147028, //2
1.17901,
1.138716, //1.1387,
1.181525,
1.138716,
1.060369, //7
1.138373, // 1.126472, //8
1.186377,//1.200398, //1.22768, //9
1.13872,
1.126547,//1.14703, //11
1.149282,//1.15973, //12
1.10664, //13
1.093714, //1.10664 //14
1.10664, // 15
1.0697201,  //16
1.100939, //1.1184544, //17
1.1483270, //1.1180530, //18
1.21346, //19
1.253052, //1.3196916, //20
1.065265, //1.022056, //1.045382, //21 !!!!
0.952887, //1.056741, //22
1.0 //23
};

//******************* dE-plane corrections ****************************************

double dEGainCorrections[24][4] = {	
{0.942847, 3.159, 2.701, 659.7}, //0
{1.01363, 2.381, 2.358, 656.7}, //1
{0.957214, 2.178, 2.171, 555.6}, //2
{0.978934, 4.54085,  3.58637, 416.94}, //3
{0.993587, 3.96687,  3.9343, 305.725}, //4
{0.943484, 4.91193,  4.96377, 240.524}, //5
{0.938617, 4.40056,  5.01474, 284.225}, //6
{0.984372, 4.47317,  5.2921, 299.741}, //7
{0.973126, 5.76316,  5.72438, 257.804}, //8
{0.972052, 5.20898,  5.73479, 250.0}, //9
{0.966815, 5.14983,  6.51716, 270.598}, //10
{0.984943, 4.36701,  4.6097, 278.071}, //11
{0.959833, 4.42663,  5.12853, 265.518}, //12
{1.20093, 5.4342,  5.26873, 321.272}, //13
{0.866577, 4.05522,  3.93413, 234.327}, //14
{1.05817, 4.44272,  3.98357, 284.014}, //15
{1.04382, 4.44544,  5.72768, 273.388}, //16
{1.08008, 3.97472,  4.38101, 243.698}, //17
{1.06031, 3.148, 4.224, 324.5}, //18
{1.08027, 4.043, 3.413, 252.7}, //19
{0.919415, 2.486,  5.764, 291.4}, //20
{0.909647, 4.028,  5.470, 287.9}, //21
{0.968645, 3.414, 4.118, 249.9}, //22
{1.0, 0.0, 0.0,0.0} //23
};

// zgornja meja je 1250.
// spodnja meja pa je 760
double dEFinalCorrections[24] = {	
2.226926, //2.1378485,
1.996275,//1.8956627,
1.761272, //1.85004,//1.9595661, //1.74323, //2
1.77088,//1.67171,//1.78655,
1.74313,  //1.794589,
2.01694,
1.687588, //1.658722, //1.63292, //6
1.567452, //1.58877,
1.881716, //1.89677, //8
1.986606,//1.965945,//2.03672, //9
1.804876,//1.852525, //10
1.669922, //1.71000, //11
1.73213, //1.78063, //12
1.434891, //13
1.983233, //2.022898, //14
1.602513, //15
1.499816, //1.545205, //16
1.789397, // 1.861444, //1.786987, //17
1.535044, //1.465563, // 1.517672, //18
1.789642, //19
1.689952, //20
1.626228, //1.532991, //1.570805, //1.560333, //21
1.859684, //1.965066, //22
1.0
};

double dEPlusFinalCorrections[24] = {	
2.226926, //2.1378485,
1.996275,//1.8956627,
1.85004,//1.9595661, //1.74323,
1.621685, //1.77088,//1.67171,//1.78655, //3
1.74313,  //1.794589,
2.01694,
1.59937, //1.53412, //1.63292, //6
1.567452, //1.58877,
1.881716, //1.89677, // 8
1.986606,//1.965945,//2.03672, //9
1.804876,//1.852525, //10
1.689493, //1.739502,//1.669922, //1.71000, //11
1.73213, //1.78063, //12
1.434891, //13
1.983233, //2.022898, //14
1.602513, //15
1.751946,//1.806452, //16
1.745816, //1.814730, //17
1.562990, //1.492445, //1.545510, //18
1.789642, //19 
1.689952, //20
1.543371, //1.625718, //1.570805, //1.560333, //21
1.8047797, //1.787133, //22
1.0
};


//#################### Vertex Time & Momentum Reconstruction #####################
//Magnetic field density in [T] ---0.693---0.84--------------------------------------------
double magnetic_field_density = 0.92;//0.92;

//Distance from target to collimator [mm] --------------------------------
double DistanceToCollimator = 1267.03;

//Distance from target to front face of BB magnet in [mm] --------------------------------
double dist_to_target = 1267.0 + 350.0; //1267.0 + 350.0;

//Position of exit face with respect to to front face of BB magnet in [mm] 1.0-------------
double exit_face_xc = 800.0;
double exit_face_yc = 0.0; //10.0


//Length and Height of BB magnet in [mm] 753-- 818----------------------------------------
double magnet_length = 1000.1;
double magnet_height = 1200.0;

//Exit face angle in [deg] -------------------------------------------------------
double exit_face_angle = 20.0;

//Multi Wire Drift Chamber angle in [deg] ----------------------------------------
double MWDC_angle = 25.0;

//Position of First MWDC plane with respect to the entrance face in [mm] 1100 in 320---------------------
double MWDC_x = 1100.0 + 150.0; // 
double MWDC_y = 300.0  + 150.0*tan(25.0*3.1415/180);

//Distance from first MWDC plane fo second MWDC plane in [mm] --------------------
double MWDC_dist_to_second = 757.0;

//Scattering angle  [rad] --------------------
double scattering_angle = -75.0 * 3.1415/180.0;

//Distance from first MWDC plane to Scintillation plane in [mm] ------------------
double distance_to_EdE  = 973.0;
//################################################################################

using namespace std;

Double_t fitf(Double_t *x, Double_t *par)
{
  Double_t G = par[0];
  Double_t L1 = par[1];
  Double_t L2 = par[2];

  return G*(exp(L1*x[0]) + exp(-1.0*L2*x[0]))/2.0;

}

TComplex mpow(TComplex z, double n)
{
 double phi = atan2(z.Im(),z.Re());
 double r = sqrt(z.Re()*z.Re() + z.Im()*z.Im());
 return TComplex( pow(r,n)*cos(phi*n),pow(r,n)*sin(phi*n));
}

TComplex msqrt(TComplex z)
{
 return  mpow(z,0.5);
}

void find_p_and_L(double xD1 ,double yD1, double zD1, double xD2, double yD2, double zD2, double &TgMomentum, double &TgTh, double &TgPh, double &TgY, TVector3 &TgP, double &pathlength, int DoPlot=0)
{
  

  //cout<<"Let's start calculating!!!"<<endl;

  TgMomentum = -12345.0;
  TgTh = -12345.0;
  TgPh = -12345.0;
  TgY = -12345.0;
  TVector3 pp0(-12345,-12345,-12345);
  TgP = pp0;
  pathlength = -12345.0;
  
  double deg2rad = 3.1415/180.0;

  double x0 = MWDC_x;
  double z0 = MWDC_y; 

  double ph = MWDC_angle*deg2rad; 

  // now we need to transform coordinates that correspond to trajectory from MWDC coordinate system to
  // spectrometer coordinate system, which has its center in the middle of front face plane of BB's magnet.

  // (zD2,xD2) --> (x2, z2) corresponds to second MWDC plane
  Double_t x2 = zD2*cos(ph) + xD2*sin(ph) + x0; 
  Double_t y2 = yD2;
  Double_t z2 = (zD2*sin(ph) - xD2*cos(ph) + z0); 

  // (zD1,xD1) --> (x1, z1) corresponds to first MWDC plane
  Double_t x1 = zD1*cos(ph) + xD1*sin(ph) + x0; 
  Double_t y1 = yD1;
  Double_t z1 = zD1*sin(ph) - xD1*cos(ph) + z0; 

  
  // now we calculate initial gradient of trajectory in spec. c-s.
  double q0 = (y2-y1)/(x2-x1);

  // exit face of  the dipole
  Double_t xC = exit_face_xc; 
  Double_t zC = exit_face_yc; 

  // length and height of dipole
  Double_t LC= magnet_length; 
  Double_t hC= magnet_height; 


  Double_t  d1 = dist_to_target;
  Double_t  alphaB = exit_face_angle*deg2rad; 


    // calc point X where the particle gets out of the magnetic field
   double  aI = (z1 - z2)/(x1 - x2);
   double  bI = (z2*x1 - z1*x2)/(x1 - x2);
    
   double  aB = -1.0/tan(alphaB);
   double  bB = zC - aB*xC;
    
   double  xX = (bB - bI) / (aI - aB);
   double  zX = aI*xX + bI;


  double  low  = zC + (xC - LC)/tan(alphaB);
  double  high = low + hC;
 

  double zT=0.0;


    // calc point Y where the particle gets in to the magnetic field
    Double_t xY = 0.0;
    
    double a0, a1, a2;
    a0 = d1*(aI*zX*zX - aI*xX*xX + 2.0*xX*zX);
    a1 = zX*zX - xX*xX - 2.0*aI*xX*zX - 
         2.0*d1*xX - 2.0*aI*d1*zX;
    a2 = aI*d1 - 2.0*zX + 2.0*aI*xX;
 
    double p, q;
    p = a1 - pow(a2, 2.0)/3.0;
    q = 2.0*pow(a2, 3.0)/27.0 - a1*a2/3.0 + a0;
	

    double dis_sqr;
    dis_sqr = pow(q/2.0, 2.0) + pow(p/3.0, 3.0);



    TComplex offset(-a2/3.0,0.0);
    TComplex dis = msqrt(TComplex(dis_sqr)); 
    TComplex u = mpow(TComplex(-q/2.0) + dis, 1.0/3.0); 
    TComplex v = TComplex(p/3.0,0.0)/u;

    TComplex s[3];
    s[0] = u-v; 
    s[1] = (v-u)/TComplex(2.0,0.0) + TComplex(0.0, sqrt(3.0)/2.0)*(u+v);  
    s[2] = (v-u)/TComplex(2.0,0.0) - TComplex(0.0, sqrt(3.0)/2.0)*(u+v); 

    s[0] = s[0] + offset;
    s[1] = s[1] + offset;
    s[2] = s[2] + offset;

   // three solutions, from that we have to choose the right one
    bool flag[3]; // flag the right route.
    int num_true = 3;
    int i;
    for (i=0; i<3; i++) flag[i]=true;
	
    
    double _r_sqr, _xJ, _zJ, _beta, theta;
    
    for (i=0; i<3; i++) {
        // check if s[i] is real
        // due too roundoff errors, we check if s[i] is close to real
        if ((s[i].Re() == 0.0) && (fabs(s[i].Im()) > 0.000001)) {
            flag[i] = false;
            num_true--;
            continue;
        }
        else if ((s[i].Im()/s[i].Re()) > 0.000001) {
            flag[i] = false;
            num_true--;
            continue;
        }
 
        // check if s[i] is in the trapezium
        if ((s[i].Re() > high) || (s[i].Re() < low)) {
            flag[i] = false;
            num_true--;
            continue;
        }
        
        // check if the short way is in the particle direction
       
        _xJ = (s[i].Re() - zX - xX/aI)/(d1/s[i].Re() - 1.0/aI);
        _zJ = -_xJ/aI + zX + xX/aI;
        _r_sqr = pow(_xJ, 2.0) + pow(s[i].Re() - _zJ, 2.0);
        
        _beta = acos(((s[i].Re() - _zJ)*(zX - _zJ) - _xJ*(xX - _xJ))/_r_sqr);  // dolocimo velikost kota
        if (((s[i].Re() - _zJ)*(xX - _xJ) + _xJ*(zX - _zJ)) < 0.0) _beta*=-1.0; // dolocimo smer vrtenja

	// check if anything goes wrong.
        if (!( ( (_beta < 0.0) && (s[i].Re() < _zJ) &&        // clockwise
               (_zJ > ((_xJ - xX)*tan(alphaB) + zX) ) ) ||
             ( (_beta > 0.0) && (s[i].Re() > _zJ) &&        // counter clockwise
               (_zJ < ((_xJ - xX)*tan(alphaB) + zX) ) ) )) {
            flag[i] = false;
            num_true--;
            continue;
        }
	
	
    }
    if (num_true != 1) {
         //printf("ERROR - num_true = %d\n", num_true);
    }
    else {
    
        i=0;
        while (!flag[i]) i++;
	
        double zY = s[i].Re(); 
        double xJ = (zY - zX - xX/aI)/(d1/zY - 1/aI);
        double zJ = -xJ/aI + zX + xX/aI;
        _r_sqr = pow(xJ, 2.0) + pow(zY - zJ, 2.0);
        double R = sqrt(_r_sqr);
        double beta = acos(((zY - zJ)*(zX - zJ) - xJ*(xX - xJ))/_r_sqr); // radians
        double xT = -d1;

        theta = atan2(-zY,xT);
        double Lxz = sqrt(pow(xY-xT,2.)+pow(zY-zT,2)) + 
            sqrt(pow(x1-xX,2.)+pow(z1-zX,2)) + 
            R*beta;     // length in plane xz // R*beta  = (2*pi*R)*(beta/2*pi)

       	// calc the momentum according to P = eBR which is correct in MKS
	double Pxz = magnetic_field_density * R *300./1000.;
        TgTh = (zY/xT);
        TgPh = 0.90*(q0);
        //TgPh = (q0);
        TgY = y1 - tan(TgPh)*Lxz;


        pathlength = Lxz*sqrt(1.0 + pow(tan(TgPh),2));
 
        double TgPz = Pxz/sqrt(1+pow(tan(TgTh),2));
        double TgPy = TgPz*tan(TgPh);
        double TgPx = TgPz*tan(TgTh);
        TgMomentum = sqrt(pow(TgPx,2)+pow(TgPy,2)+pow(TgPz,2));

        TVector3 pp(TgPx, TgPy, TgPz);
        TgP = pp;  

    
  

        if (DoPlot >-1){
          TCanvas *c1 = new TCanvas("c1","c1");
          c1->cd(1);
          //gPad->Range(-2000,-1000,3000,2000);
          gPad->SetGrid();

          double tockiX[2] = {-2000,3000};
          double tockiY[2] = {-1000,2000};
          TGraph *g1 = new TGraph(2, tockiX, tockiY);
          g1->Draw("AP");

          TLine *lineEntranceFace = new TLine(0,-1000,0,3000);
          lineEntranceFace->Draw();
	  TMarker *pointMWDC = new TMarker(x0,z0, 20);
  	  pointMWDC->Draw();

  	  TMarker *point1 = new TMarker(x1,z1, 22);
  	  point1->Draw();
  	  TMarker *point2 = new TMarker(x2,z2, 22);
  	  point2->Draw();

  	  TMarker *pointExitFace = new TMarker(xC,zC, 20);
  	  pointExitFace->Draw();

          TMarker *pointTarget = new TMarker(xT,0, 20);
          pointTarget->Draw();
  
          TMarker *point3 = new TMarker(xJ,zJ, 22);
          point3->SetMarkerColor(kRed);
          point3->Draw();
 
          TMarker *point4 = new TMarker(xY,zY, 22);
          point4->Draw();

          TEllipse *circRadius = new TEllipse(xJ,zJ,R,R);
          circRadius->SetLineColor(kRed);
          circRadius->SetFillStyle(0);
          circRadius->Draw();

          TMarker *point5 = new TMarker(xX,zX, 22);
   	  point5->Draw();

  	  TLine *lineExitFace = new TLine(exit_face_xc-3000*tan(alphaB),3000+exit_face_yc, exit_face_xc+1000*tan(alphaB),-1000+exit_face_yc);
  	  lineExitFace->Draw();

          TLine *lineMWDC1 = new TLine(x0-700*sin(ph),z0+700*cos(ph), x0+700*sin(ph),z0-700*cos(ph));
          lineMWDC1->SetLineColor(kGreen);
  	  lineMWDC1->Draw();

          TLine *lineMWDC2 = new TLine(x0-1000*sin(ph) + MWDC_dist_to_second*cos(ph),z0+1000*cos(ph)+ MWDC_dist_to_second*sin(ph), x0+1000*sin(ph)+ MWDC_dist_to_second*cos(ph),z0-1000*cos(ph) + MWDC_dist_to_second*sin(ph));
          lineMWDC2->SetLineColor(kGreen);
  	  lineMWDC2->Draw();

          TLine *lineTrack2 = new TLine(xX,zX,x2,z2);
          lineTrack2->SetLineColor(kRed);
  	  lineTrack2->Draw();

          TLine *lineTrack1 = new TLine(xT,zT,xY,zY);
          lineTrack1->SetLineColor(kRed);
  	  lineTrack1->Draw();

          TLine *lineMagnetLow = new TLine(0,low, LC,low);
  	  lineMagnetLow->Draw();
      
          TLine *lineMagnetHigh = new TLine(0,high, xC,high);
  	  lineMagnetHigh->Draw();

 
          cout<<"Target Y: "<<TgY<<endl;
          cout<<"Target Th: "<<TgTh<<endl;
          cout<<"MagnetiFieldDensity: "<<magnetic_field_density<<endl;
          cout<<"Radius: "<<R<<endl;
          cout<<"Momentum: "<<TgMomentum<<endl;
       }
 







    }

  return; 
}

Double_t DeuteronLine(Double_t *x, Double_t *par)
{
  Double_t a = par[0];
  Double_t b = par[1];
  Double_t c = par[2];
  Double_t k = par[2];

  return a + b*exp(-1.0*(x[0]-c)/k);

}


void BigBiteAnalyticalModelTrack( double trackX = -0.0, double trackY = 0.0, double trackPhi = 0.0, double trackTheta = 0.0){

/*
  double trackX = -0.5; //-0.032;
  double trackY = 0.0;
  double trackPhi = 0.0;
  double trackTheta = -0.11;//-0.011;
*/
  // now that we have all variables that we need, we can proceed and calculate
  // the momentum, pathlength and vertex time

  double x1 = trackX*1000.0;
  double y1 = trackY*1000.0;
  double z1 = 0.0;
  double th = trackTheta;
  double ph = trackPhi;
  double z2 = MWDC_dist_to_second; // 0.36
  double x2 = x1 + th*(z2-z1); /// th - means tan(th)
  double y2 = y1 + ph*(z2-z1); /// ph - means tan(ph)   


  double TgMomentum, TgY, TgTh, TgPh;
  double pathlength;
  TVector3 p;

  find_p_and_L(x1 , y1, z1, x2, y2, z2,  TgMomentum, TgTh, TgPh, TgY, p, pathlength, 0);



}



void BigBiteAnalyticalModelAnalysis(int MAXnum = -1){

  Double_t BB_tr_x[10];
  Double_t BB_tr_y[10];
  Double_t BB_tr_th[10];
  Double_t BB_tr_ph[10];
  Double_t L_gold_p;
  Double_t scattering_angle;
  Double_t rb_x;
  Double_t rb_y;
  Double_t ReactPt_L_z;
  Double_t DL_t3[100];
  Double_t DL_t1[100];

  Double_t PriKineLHe3_qx;
  Double_t PriKineLHe3_qy;
  Double_t PriKineLHe3_qz;

  Double_t BB_tp_e_LE[24];
  Double_t BB_tp_e_RE[24];  
  Int_t Ndata_BB_tp_e_LE;
  Int_t Ndata_BB_tp_e_RE;  
  Double_t BB_tp_de_LE[24];
  Double_t BB_tp_de_RE[24];  
  Int_t Ndata_BB_tp_de_LE;
  Int_t Ndata_BB_tp_de_RE;  
  Double_t BB_tp_trx[100];
  Double_t BB_tp_try[100];
  Double_t BB_tp_trHitIndex[100];
  Double_t BB_tp_e_hit_bar[100];
  Double_t BB_tp_de_hit_bar[100];
  Int_t Ndata_DL_t3;  
  Int_t Ndata_BB_tp_e_hit_bar;  
  Int_t Ndata_BB_tp_de_hit_bar;  



  TChain *t1 =new TChain("tout");
  //t1->Add("./ROOTResults/eed_Test_GoodEvents_3265.root");
  //t1->Add("./ROOTResults/eed_Test_GoodEvents_3608.root");
  //t1->Add("./ROOTResults/eed_Test_GoodEvents_3151.root");
  //t1->Add("./ROOTResults/eed_Test_GoodEvents_2413.root");
  //t1->Add("./ROOTResults/eed_Test_GoodEvents_1518.root");
  //t1->Add("./ROOTResults/eed_Test_GoodEvents_2488.root");
  t1->Add("./ROOTResults/eed_Test2_GoodEvents_3488.root");
  //t1->Add("./ROOTResults/eed_Test_GoodEvents_3491.root");
  //t1->Add("./ROOTResults/eed_Test_GoodEvents_3490.root");
  //t1->Add("./ROOTResults/eed_Test_GoodEvents_3489.root");
  //t1->Add("./ROOTResults/eed_Test_GoodEvents_2164.root");

  //t1->Add("./ROOTResults/eed_Test2_GoodEvents_2240.root");
  //t1->Add("./ROOTResults/eed_Test2_GoodEvents_2241.root");
  //t1->Add("./ROOTResults/eed_Test2_GoodEvents_2242.root");
  //t1->Add("./ROOTResults/eed_Test2_GoodEvents_2243.root");

  //t1->Add("./ROOTResults/eed_Test2_GoodEvents_2258.root");
  //t1->Add("./ROOTResults/eed_Test2_GoodEvents_2259.root");
  //t1->Add("./ROOTResults/eed_Test2_GoodEvents_2260.root");


  assert(t1);

  t1->SetBranchAddress( "L.gold.p", &L_gold_p);
  t1->SetBranchAddress( "BB.tr.x", &BB_tr_x);
  t1->SetBranchAddress( "BB.tr.y", &BB_tr_y);
  t1->SetBranchAddress( "BB.tr.th", &BB_tr_th);
  t1->SetBranchAddress( "BB.tr.ph", &BB_tr_ph);
  t1->SetBranchAddress( "scattering_angle", &scattering_angle);
  t1->SetBranchAddress( "ReactPt_L.z", &ReactPt_L_z);
  t1->SetBranchAddress( "rb.x", &rb_x);
  t1->SetBranchAddress( "rb.y", &rb_y);
  t1->SetBranchAddress( "DL.t3", &DL_t3); 
  t1->SetBranchAddress( "DL.t1", &DL_t1);
  t1->SetBranchAddress( "Ndata.DL.t3", &Ndata_DL_t3);  

  t1->SetBranchAddress( "BB.tp.e.LE", &BB_tp_e_LE);
  t1->SetBranchAddress( "BB.tp.e.RE", &BB_tp_e_RE);
  t1->SetBranchAddress( "Ndata.BB.tp.e.LE", &Ndata_BB_tp_e_LE);
  t1->SetBranchAddress( "Ndata.BB.tp.e.RE", &Ndata_BB_tp_e_RE);
  t1->SetBranchAddress( "BB.tp.de.LE", &BB_tp_de_LE);
  t1->SetBranchAddress( "BB.tp.de.RE", &BB_tp_de_RE);
  t1->SetBranchAddress( "Ndata.BB.tp.de.LE", &Ndata_BB_tp_de_LE);
  t1->SetBranchAddress( "Ndata.BB.tp.de.RE", &Ndata_BB_tp_de_RE);
  t1->SetBranchAddress( "BB.tp.trx", &BB_tp_trx);
  t1->SetBranchAddress( "BB.tp.try", &BB_tp_try); 
  t1->SetBranchAddress( "BB.tp.e.hit_bar", &BB_tp_e_hit_bar);
  t1->SetBranchAddress( "BB.tp.de.hit_bar", &BB_tp_de_hit_bar);
  t1->SetBranchAddress( "BB.tp.trHitIndex", &BB_tp_trHitIndex);
  t1->SetBranchAddress( "Ndata.BB.tp.de.hit_bar", &Ndata_BB_tp_de_hit_bar);
  t1->SetBranchAddress( "Ndata.BB.tp.e.hit_bar", &Ndata_BB_tp_e_hit_bar);



  t1->SetBranchAddress( "PriKineLHe3.q_x", &PriKineLHe3_qx);
  t1->SetBranchAddress( "PriKineLHe3.q_y", &PriKineLHe3_qy);
  t1->SetBranchAddress( "PriKineLHe3.q_z", &PriKineLHe3_qz);

  int num = t1->GetEntries();
  cout<<"Number of entries : "<<num<<endl;
  if (MAXnum > 0 &&  MAXnum < num ) num = MAXnum;



  TH1F *hMomentum = new TH1F("hMomentum", "hMomentum", 400,0,1000);
  TH1F *hpx = new TH1F("hpx", "hpy", 400,-200,200);
  TH1F *hpy = new TH1F("hpy", "hpy", 400,-200,200);
  TH1F *hpz = new TH1F("hpz", "hpz", 400, 200,600);

  TH1F *hqVec = new TH1F("hqVec", "hqVec", 400,0,1000);
  TH1F *hqx = new TH1F("hqx", "hqx", 400,-200,200);
  TH1F *hqy = new TH1F("hqy", "hqy", 400,-200,200);
  TH1F *hqz = new TH1F("hqz", "hqz", 400,200.0,600);

  TH1F *hTgY = new TH1F("hTgY", "hTgY", 400, -400, 400);
  TH1F *hTgTh = new TH1F("hTgTh", "hTgTh", 400, -0.4, 0.4);
  TH1F *hTgPh = new TH1F("hTgPh", "hTgPh", 400, -0.2, 0.2);

  TH1F *hQY = new TH1F("hQY", "hQY", 400, -400, 400);
  TH1F *hQTh = new TH1F("hQTh", "hQTh", 400, -0.4, 0.4);
  TH1F *hQPh = new TH1F("hQPh", "hQPh", 400, -0.2, 0.2);

  TH2F *hTgTh_vs_QTh = new TH2F("hTgTh_vs_QTh", "hTgTh_vs_QTh", 400, -0.4, 0.4, 400, -0.4, 0.4);
  TH2F *hTgPh_vs_QPh = new TH2F("hTgPh_vs_QPh", "hTgPh_vs_QPh", 400, -0.2, 0.2, 400, -0.2, 0.2);
  TH2F *hTgY_vs_QY = new TH2F("hTgY_vs_QY", "hTgY_vs_QY", 400, -400, 400, 400, -400, 400);
  TH2F *hTgMomentum_vs_DLt3 = new TH2F("hTgMomentum_vs_DLt3", "hTgMomentum_vs_DLt3", 400, 0.0, 1000.0, 140, 760, 900);

  TH2F *hColliXY = new TH2F("hColliXY","hColliXY", 500,-150.0,150.0,500,-400.0, 400.0);

  TH2F *hADCEdE = new TH2F("hADCEdE", "hADCEdE", 500,0,3500,500,0,3500);
  TH2F *hADCE_vs_p = new TH2F("hADCE_vs_p", "hADCE_vs_p", 500,0,3200,500,0,1000);
  TH2F *hADCE_vs_p_Deuterons = new TH2F("hADCE_vs_p_Deuterons", "hADCE_vs_p_Deuterons", 500,0,3200,500,0,1000);

  TH2F *hADCE_vs_DLt3 = new TH2F("hADCE_vs_DLt3", "hADCE_vs_DLt3", 500,0,3200,200,700,900);
  TH2F *hADCdE_vs_p = new TH2F("hADCdE_vs_p", "hADCdE_vs_p", 500,0,3200,500,0,1000);
  TH2F *hADCdE_vs_DLt3 = new TH2F("hADCdE_vs_DLt3", "hADCdE_vs_DLt3", 500,0,3200,200,700,900);

  TH1F *hMassBigBite = new TH1F("hMassBigBite", "hMassBigBite", 100,600,2200);
  TH1F *hMissingMass = new TH1F("hMissingMass", "hMissingMass", 300,600,2200);
  TH2F *hMissingMass_vs_momentum = new TH2F("hMissingMass_vs_momentum", "hMissingMass_vs_momentum", 300,600,2200, 300,200,800);

  TH2F *hFpPh_vs_QPh = new TH2F("hFpPh_vs_QPh", "hFpPh_vs_QPh", 400, -0.2, 0.2, 400, -0.2, 0.2);

  TH1F *hTgPathLength = new TH1F("hTgPathLength", "hTgPathLength", 400, 0, 4000);
  TH2F *hTgPathLength_vs_FpX = new TH2F("hTgPathLength_vs_FpX", "hTgPathLength_vs_FpX", 200, 2500, 3500, 200,-800, 800);
  TH1F *hTotalPathLength = new TH1F("hTotalPathLength", "hTotalPathLength", 400, 0, 4500);
  TH2F *hTotalPathLength_vs_EPlane = new TH2F("hTotalPathLength_vs_EPlane", "hTotalPathLength_vs_EPlane", 200, 3500, 4500, 25,0, 24);
  
  TH1F *hTimeOfFlight = new TH1F("hTimeOfFlight", "hTimeOfFlight", 400, 0, 150);
  TH2F *hTgMomentum_vs_TimeOfFlight = new TH2F("hTgMomentum_vs_TimeOfFlight", "hTgMomentum_vs_TimeOfFlight", 400, 0.0, 1000.0, 400, 0, 150);
  TH2F *hT3_vs_TimeOfFlight = new TH2F("hT3_vs_TimeOfFlight", "hT3_vs_TimeOfFlight", 180, 350.0, 440.0, 400, 0, 150);

  for (int i = 0; i<num; i++){
    t1->GetEntry(i);
    if ((i%1000==0) && i>0) cout<<i<<endl; 
 
    if ( DL_t3[0]< 780.0 || DL_t3[0] > 880.0 ) continue;
    if (Ndata_BB_tp_e_hit_bar !=1 || Ndata_BB_tp_de_hit_bar !=1) continue;

    double FpY = BB_tp_try[0];
    double FpX = BB_tp_try[0];
    int trHitIndex = (int) BB_tp_trHitIndex[0];
    int EplaneHit = (int)BB_tp_e_hit_bar[trHitIndex];
    int dEplaneHit = (int)BB_tp_de_hit_bar[trHitIndex];
    if ((EplaneHit!=dEplaneHit-1) && (EplaneHit!=dEplaneHit)) continue; 
    if (EplaneHit>22 || dEplaneHit>22 ) continue; 
    double EmeanADCValue= 0.5*(EGainCorrections[EplaneHit][0]*BB_tp_e_RE[EplaneHit]+BB_tp_e_LE[EplaneHit]);
    double ResParams[3];
    ResParams[0] = EGainCorrections[EplaneHit][3];
    ResParams[1] = EGainCorrections[EplaneHit][1];
    ResParams[2] = EGainCorrections[EplaneHit][2];
    double EmeanADCValueCorrected = EFinalCorrections[EplaneHit]*(EmeanADCValue - (fitf(&FpY, ResParams) - ResParams[0]));

    double dEmeanADCValue = 0.5*(dEGainCorrections[dEplaneHit][0]*BB_tp_de_RE[dEplaneHit]+BB_tp_de_LE[dEplaneHit]);

    ResParams[0] = dEGainCorrections[dEplaneHit][3];
    ResParams[1] = dEGainCorrections[dEplaneHit][1];
    ResParams[2] = dEGainCorrections[dEplaneHit][2];

    double dEmeanADCValueCorrected = 0.0;
    if ((EplaneHit==dEplaneHit)) dEmeanADCValueCorrected = dEFinalCorrections[dEplaneHit]*(dEmeanADCValue - (fitf(&FpY, ResParams) - ResParams[0]));
    if ((EplaneHit==dEplaneHit-1)) dEmeanADCValueCorrected = dEPlusFinalCorrections[dEplaneHit]*(dEmeanADCValue - (fitf(&FpY, ResParams) - ResParams[0]));



    TVector3 TCSX(0,-1,0);
    TVector3 TCSZ(TMath::Sin(scattering_angle),0,TMath::Cos(scattering_angle));
    TVector3 TCSY = TCSZ.Cross(TCSX);
    TRotation  fTCSInHCS;
    fTCSInHCS.RotateAxes(TCSX,TCSY,TCSZ);

    TVector3 BeamSpotHCS( rb_x, rb_y, ReactPt_L_z );
    TVector3 BeamSpotTCS=fTCSInHCS.Inverse()*(BeamSpotHCS);
    

    TVector3 qVectorHCS(PriKineLHe3_qx, PriKineLHe3_qy, PriKineLHe3_qz);
    TVector3 qVectorTCS = fTCSInHCS.Inverse()*(qVectorHCS);

    double trackX = BB_tr_x[0];
    double trackY = BB_tr_y[0] + 0.024;
    double trackPhi = BB_tr_ph[0];
    double trackTheta = BB_tr_th[0];

    // now that we have all variables that we need, we can proceed and calculate
    // the momentum, pathlength and vertex time

    double x1 = trackX*1000.0;
    double y1 = trackY*1000.0;
    double z1 = 0.0;
    double th = trackTheta;
    double ph = trackPhi;
    double z2 = MWDC_dist_to_second; // 0.36
    double x2 = x1 + th*(z2-z1); /// th - means tan(th)
    double y2 = y1 + ph*(z2-z1); /// ph - means tan(ph)   
	
    double TgMomentum, TgY, TgTh, TgPh, PathLength;
    TVector3 p;

    find_p_and_L(x1 , y1, z1, x2, y2, z2,  TgMomentum, TgTh, TgPh, TgY, p, PathLength, -1);


    if (fabs(1000.0*BeamSpotTCS.Y() - TgY)>32.0) continue;

    double MBB = 0.0;

    // 3He mass = 2809.41 MeV/c2
    double DeuteronPar[4];
    DeuteronPar[0]  = 400.0;
    DeuteronPar[1]  = 434.874;
    DeuteronPar[2]  = 428.0;
    DeuteronPar[3]  = 408.636;
    
    if (TgMomentum> DeuteronLine(&dEmeanADCValueCorrected, DeuteronPar) ) MBB = 1875.614; // deuterons
    else  MBB = 938.27; //protons
    

    double EBB = sqrt(pow(MBB,2) + pow(p.Mag(),2));

    TVector3 pmm(1000.0*qVectorTCS.X() - p.X(), 1000.0*qVectorTCS.Y() - p.Y(), 1000.0*qVectorTCS.Z() - p.Z());
    double omega = 2425.0 - 1000.0*L_gold_p;
    double Mmm = sqrt(pow(omega + 2809.41- EBB,2) - pow(pmm.Mag(),2));


    double path_from_MWDC_to_EdE = sqrt(pow(trackX - FpX,2)+pow(trackY - FpY,2)+pow(distance_to_EdE/1000.0,2));
    double total_path_length = 1000.0*path_from_MWDC_to_EdE + PathLength; 
    
    double beta = p.Mag()/EBB;
    double velocity = beta*3e8;
    double time_of_flight = 1e9*total_path_length/1000.0/velocity;
    

    hMomentum->Fill(TgMomentum);
    hpx->Fill(p.X()); 
    hpy->Fill(p.Y()); 
    hpz->Fill(p.Z()); 


    hTgY->Fill(TgY);
    hTgTh->Fill(TgTh);
    hTgPh->Fill(TgPh);

    double QTh = qVectorTCS.X()/qVectorTCS.Z()+0.022;
    double QPh = qVectorTCS.Y()/qVectorTCS.Z()-0.022;

    hQY->Fill(BeamSpotTCS.Y()*1000.0);
    hQTh->Fill(QTh);
    hQPh->Fill(QPh);
    

    hqVec->Fill(1000*qVectorTCS.Mag());
    hqx->Fill(1000*qVectorTCS.Mag()*QTh/(sqrt(1 + pow(QTh,2) + pow(QPh,2)))); 
    hqy->Fill(1000*qVectorTCS.Mag()*QPh/(sqrt(1 + pow(QTh,2) + pow(QPh,2)))); 
    hqz->Fill(1000*qVectorTCS.Mag()/(sqrt(1 + pow(QTh,2) + pow(QPh,2)))); 


    hTgTh_vs_QTh->Fill(qVectorTCS.X()/qVectorTCS.Z(), TgTh) ;
    hTgPh_vs_QPh->Fill(qVectorTCS.Y()/qVectorTCS.Z(), TgPh) ;
    hTgY_vs_QY->Fill( 1000.0*BeamSpotTCS.Y(), TgY) ;


    if (fabs(qVectorTCS.X()/qVectorTCS.Z())<0.01) hFpPh_vs_QPh->Fill(qVectorTCS.Y()/qVectorTCS.Z(), TgPh);

    hTgMomentum_vs_DLt3->Fill(TgMomentum, DL_t3[0]);

    double ColliX = DistanceToCollimator*tan(TgTh);
    double ColliY = TgY + DistanceToCollimator*tan(TgPh);

    hColliXY->Fill(-ColliY, -ColliX);

    hADCE_vs_p->Fill(EmeanADCValueCorrected, TgMomentum);
    if (MBB >1000) hADCE_vs_p_Deuterons->Fill(EmeanADCValueCorrected, TgMomentum);
    hADCE_vs_DLt3->Fill(EmeanADCValueCorrected, DL_t3[0]);

    hADCdE_vs_p->Fill(dEmeanADCValueCorrected, TgMomentum);
    hADCdE_vs_DLt3->Fill(dEmeanADCValueCorrected, DL_t3[0]);

    hADCEdE->Fill(EmeanADCValueCorrected, dEmeanADCValueCorrected);

    hMassBigBite->Fill(MBB);
    hMissingMass->Fill(Mmm);
    hMissingMass_vs_momentum->Fill(Mmm, TgMomentum);

    hTgPathLength->Fill(PathLength);
    hTgPathLength_vs_FpX->Fill(PathLength, x1);

    hTotalPathLength->Fill(total_path_length);
    hTotalPathLength_vs_EPlane->Fill(total_path_length,  EplaneHit );

    hTimeOfFlight->Fill(time_of_flight);

    hTgMomentum_vs_TimeOfFlight->Fill(TgMomentum, time_of_flight);
    hT3_vs_TimeOfFlight->Fill(0.5*DL_t3[0], time_of_flight);

  

  }

  TCanvas *c2 = new TCanvas("c2","c2");
  c2->Divide(2,2);
  c2->cd(1);
  hqVec->SetLineColor(kRed);
  hqVec->Draw();
  hMomentum->Draw("same");
  c2->cd(2); 
  hqx->SetLineColor(kRed);
  hqx->Draw();
  hpx->Draw("same");
  c2->cd(3); 
  hqy->SetLineColor(kRed);
  hqy->Draw();
  hpy->Draw("same");
  c2->cd(4); 
  hqz->SetLineColor(kRed);
  hqz->Draw();
  hpz->Draw("same");

  c2->Print("figure_BigBiteAnalyticalModelAnalysis_2.png");

  TCanvas *c3 = new TCanvas("c3","c3");
  c3->Divide(2,3);
  c3->cd(1);
   hTgTh->Draw();
  hQTh->SetLineColor(kRed);
  hQTh->Draw("same");
 
  c3->cd(2);
  gPad->SetGrid();
  gStyle->SetPalette(1);
  hTgTh_vs_QTh->Draw("colz");
  TLine *l1 = new TLine(-0.4,-0.4, 0.4, 0.4);
  l1->SetLineColor(kRed);
  l1->SetLineWidth(2);
  l1->Draw("same");

  c3->cd(3); 
  hQPh->SetLineColor(kRed);
  hQPh->Draw();
  hTgPh->Draw("same");
  c3->cd(4);
  gPad->SetGrid();
  gStyle->SetPalette(1);
  hTgPh_vs_QPh->Draw("colz");
  TLine *l2 = new TLine(-0.2,-0.2, 0.2, 0.2);
  l2->SetLineColor(kRed);
  l2->SetLineWidth(2);
  l2->Draw("same");


  c3->cd(5);
  hQY->SetLineColor(kRed);
  hQY->Draw();
  hTgY->Draw("same");
  c3->cd(6);
  gPad->SetGrid();
  gStyle->SetPalette(1);
  hTgY_vs_QY->Draw("colz");
  TLine *l3 = new TLine(-400.0,-400.0, 400., 400.);
  l3->SetLineColor(kRed);
  l3->SetLineWidth(2);
  l3->Draw("same");


  c3->Print("figure_BigBiteAnalyticalModelAnalysis_3.png");

  TCanvas *c4 = new TCanvas("c4","c4");
  c4->Divide(1,1);
  c4->cd(1);
  gPad->SetGrid();
  gStyle->SetPalette(1);
  hTgMomentum_vs_DLt3->Draw("colz");

  c4->Print("figure_BigBiteAnalyticalModelAnalysis_4.png");

  TCanvas *c5 = new TCanvas("c5","c5");
  c5->Divide(1,1);
  c5->cd(1);
  gPad->SetGrid();
  gStyle->SetPalette(1);
  hColliXY->Draw("colz");
  
  double offsetY = 0.0;
  double offsetX = 0.0;
  double in2mm = 0.0254*1000.0;
  

  TLine *HorizontalLines[13];
 
  for (int i = 0; i<13; i++)
  { 
    HorizontalLines[i] = new TLine(-7.375*in2mm + offsetX, (i-6)*1.938*in2mm + offsetY, 7.375*in2mm + offsetX, (i-6)*1.938*in2mm + offsetY);
    HorizontalLines[i]->SetLineColor(kRed);
    HorizontalLines[i]->SetLineWidth(1.5);
    HorizontalLines[i]->Draw("same");
  }

  TLine *VerticalLines[7];

  for (int i = 0; i<7; i++)
  { 
    VerticalLines[i] = new TLine((i-3)*1.5*in2mm + offsetX, -13.75*in2mm + offsetY, (i-3)*1.5*in2mm + offsetX, 13.75*in2mm + offsetY);
    VerticalLines[i]->SetLineColor(kRed);
    VerticalLines[i]->SetLineWidth(1.5);
    VerticalLines[i]->Draw("same");
  }

  double Tocke[91][2];
  TEllipse *Holes[91];

  for (int i = 0; i<13; i++)
  { 
    for (int j = 0; j<7; j++)
    { 
      Tocke[j+ i*7][0] =  (j-3)*1.5*in2mm + offsetX;
      Tocke[j+ i*7][1] =  (i-6)*1.938*in2mm + offsetY;
      Holes[j+ i*7] = new TEllipse(Tocke[j+ i*7][0], Tocke[j+ i*7][1], 0.009525*1000.0, 0.009525*1000.0);
      Holes[j+ i*7]->SetLineColor(kRed);
      Holes[j+ i*7]->SetLineWidth(1);
      Holes[j+ i*7]->SetFillStyle(0);
      Holes[j+ i*7]->Draw("same");
    }
  }

  c5->Print("figure_BigBiteAnalyticalModelAnalysis_5.png");

  TCanvas *c6 = new TCanvas("c6","c6");
  c6->Divide(2,1);
  c6->cd(1);
  gPad->SetGrid();
  gStyle->SetPalette(1);
  hADCE_vs_p->Draw("colz");
  c6->cd(2);
  gPad->SetGrid();
  gStyle->SetPalette(1);
  hADCdE_vs_p->Draw("colz");

  TF1 *fDeuterons = new TF1("fDeuterons", DeuteronLine, 400, 2000, 4);
  fDeuterons->SetParNames("a","b","c","k");
  fDeuterons->FixParameter(0, 400.0);
  fDeuterons->FixParameter(1, 434.874);
  fDeuterons->FixParameter(2, 428.0);
  fDeuterons->FixParameter(3, 408.636);
  fDeuterons->SetLineColor(kRed);
  fDeuterons->SetLineWidth(1.3);
  fDeuterons->Draw("same");

  c6->Print("figure_BigBiteAnalyticalModelAnalysis_6.png");
  
  TCanvas *c6a = new TCanvas("c6a","c6a");
  c6a->Divide(1,1);
  c6a->cd(1);
  gPad->SetGrid();
  gStyle->SetPalette(1);
  hADCE_vs_p->Draw("");
  hADCE_vs_p_Deuterons->SetMarkerColor(kRed);
  hADCE_vs_p_Deuterons->Draw("same");

  c6a->Print("figure_BigBiteAnalyticalModelAnalysis_6a.png");

  TCanvas *c7 = new TCanvas("c7","c7");
  c7->Divide(2,1);
  c7->cd(1);
  gPad->SetGrid();
  gStyle->SetPalette(1);
  hADCE_vs_DLt3->Draw("colz");
  c7->cd(2);
  gPad->SetGrid();
  gStyle->SetPalette(1);
  hADCdE_vs_DLt3->Draw("colz");

  c7->Print("figure_BigBiteAnalyticalModelAnalysis_7.png");

  TCanvas *c8 = new TCanvas("c8","c8");
  c8->Divide(1,1);
  c8->cd(1);
  gPad->SetGrid();
  gStyle->SetPalette(1);
  hADCEdE->Draw("colz");

  c8->Print("figure_BigBiteAnalyticalModelAnalysis_8.png");

  TCanvas *c9 = new TCanvas("c9","c9");
  c9->Divide(1,2);
  c9->cd(1);
  gPad->SetGrid();
  hMassBigBite->Draw();
  c9->cd(2);
  gPad->SetGrid();
  hMissingMass->Draw();

  c9->Print("figure_BigBiteAnalyticalModelAnalysis_9.png");

  TCanvas *c10 = new TCanvas("c10","c10");
  c10->Divide(1,2);
  c10->cd(1);
  gPad->SetGrid();
  hMissingMass->Draw();
  c10->cd(2);
  gPad->SetGrid();
  gStyle->SetPalette(1);
  hMissingMass_vs_momentum->Draw("colz");

  c10->Print("figure_BigBiteAnalyticalModelAnalysis_10.png");

  TCanvas *c11 = new TCanvas("c11","c11");
  c11->Divide(1,1);
  c11->cd(1);
  gPad->SetGrid();
  gStyle->SetPalette(1);
  hFpPh_vs_QPh->Draw("colz");

  c11->Print("figure_BigBiteAnalyticalModelAnalysis_11.png");

  TCanvas *c12 = new TCanvas("c12","c12");
  c12->Divide(2,2);
  c12->cd(1);
  gPad->SetGrid();
  hTgPathLength->Draw("");
  c12->cd(2);
  gPad->SetGrid();
  gStyle->SetPalette(1);
  hTgPathLength_vs_FpX->Draw("colz");
  c12->cd(3);
  gPad->SetGrid();
  hTotalPathLength->Draw("");
  c12->cd(4);
  gPad->SetGrid();
  gStyle->SetPalette(1);
  hTotalPathLength_vs_EPlane->Draw("colz");

  c12->Print("figure_BigBiteAnalyticalModelAnalysis_12.png");

  TCanvas *c13 = new TCanvas("c13","c13");
  c13->Divide(1,1);
  c13->cd(1);
  gPad->SetGrid();
  hTimeOfFlight->Draw("");

  c13->Print("figure_BigBiteAnalyticalModelAnalysis_13.png");

  TCanvas *c14 = new TCanvas("c14","c14");
  c14->Divide(1,1);
  c14->cd(1);
  gPad->SetGrid();
  gStyle->SetPalette(1);
  hTgMomentum_vs_TimeOfFlight->Draw("colz");  

  c14->Print("figure_BigBiteAnalyticalModelAnalysis_14.png");

  TCanvas *c15 = new TCanvas("c15","c15");
  c15->Divide(1,1);
  c15->cd(1);
  gPad->SetGrid();
  gStyle->SetPalette(1);
  hT3_vs_TimeOfFlight->Draw("colz");  
  
  c15->Print("figure_BigBiteAnalyticalModelAnalysis_15.png");

}