#include "TTree.h"
#include "TFile.h"
#include <math.h>
#include <iostream>
//all lengths are in mm, angles in rad
const Double_t light_c = 299792458.;
const Double_t PI = 3.14159265358979323846;
const Double_t target_shift = 0.;
//new, I try to use now also the shift, to get more out of it....
const Double_t xMW1_shift = 0.344833+0.7626;
const Double_t xMW2_shift = 0.451628+0.9026;
//--------------------------------------------------------------
const Double_t xMW3_shift = 79.84;
const Double_t yMW1_shift = 1.37840e+01;
const Double_t yMW2_shift = 3.22703e+01;
const Double_t yMW3_shift = -6.70074e+00;
//original value of alpha_g is 14 degrees
const Double_t alpha_G = (14./180.)*PI;
//this is original value of zGm-----
const Double_t zGm = 2577.;
//----------------------------------
const Double_t zMW0 = -2520.;
//probably START and MW0 Detectors are switched, have to check...
//const Double_t zMW0 = -1867.75;
//new, I try to use now also the shift, to get more out of it....
const Double_t xMW0_shift = 5.8959;
const Double_t middle_zM3 = 5937;
const Double_t middle_xM3 = -(middle_zM3-zGm)*tan(PI/10.);
const Double_t middle_zToFW = (6100+760-2777)*cos(PI/10.)+2777;
const Double_t middle_xToFW = -(6100+760-2777)*sin(PI/10.);
const Double_t cutoff_ToFW = middle_xToFW - middle_zToFW*tan(2*PI/5.);
const Double_t zT = -684.5;
const Double_t zM1 = 279.;  //to check again....
const Double_t zM2 = 854.;
const Double_t psi_out_0 = PI/10.;
const Double_t start_to_target = 1183.25;
//const Double_t time_start_target = (1183.25/(0.714549*light_c))*pow(10,6);
const Double_t time_start_target = (1162./(0.714549*light_c))*pow(10,6);
const Double_t current = 1444;
const Double_t gamma_given = 1.42942;
//const Double_t current = 1498;
const Double_t zGm_cutoff = -tan(PI/2.-alpha_G)*(zGm);
Int_t entries_mw0 = 0.;
Int_t entries_mw1 = 0.;
Int_t entries_mw2 = 0.;
Int_t entries_mw3 = 0.;
Int_t entries_start = 0.;
Int_t entries_tofw = 0.;
Int_t entries_califa = 0.;
Double_t tot_length = 0.;
Double_t charge_val = 0.;
Double_t psi_out;
//Declare more functions here...
Double_t zC;
Double_t xC;
Double_t z_labMW3;
Double_t x_labMW3;
Double_t slope;
Double_t offset_slope;
Double_t psi_out_rec[50];
Double_t z_pos_shift[50];
Double_t x_pos_shift[50];
Double_t zB;
Double_t xB;
Double_t z_D;
Double_t x_D;
//original value of Leff!!----
const Double_t Leff = 2067;
//----------------------------
Double_t l_diff[50];
const Double_t zGLAD1 = (Leff/2.)/(cos(alpha_G));
const Double_t bGLAD_cutoff = -tan(PI/2.-alpha_G)*(zGm-zGLAD1);
const Double_t D_cutoff = -tan(PI/2.-alpha_G)*(zGm+zGLAD1);

void latest_automatic(char const count_i[50]){
char fname[200];
char f_out_name[200];
char hist_name[200];

//these are the parameters needed for the time of flight calculation
Long64_t nHits =0;
Long64_t nHitsStart = 0;
UShort_t iCh;//Channel number of the TOFD per Detector
UShort_t iDet;//detector number of the TOFD
Int_t NbDets = 28; //detector number of the TOFD
Int_t NbChs = 2; //number of preams per detector plastic
Double_t iRawTimeNs[NbDets * 2];
Double_t iRawTimeStartNs[2];
UShort_t mult[NbDets * NbChs];

////these are the new offsets when doing exact measurement for xD, zD and zTOFW for RUN 36
double tof_offs[27] = {
114.4,
113.588,
116.132,
114.358,
115.261,
115.886,
116.734,
115.982,
114.222,
115.801,
112.402,
115.786,
113.347,
114.917,
118.516,
117.823,
119.217,
116.913,
118.063,
117.591,
121.295,
121.690,
120.896,
120.931,
119.649,
119.225,
119.274
};

//these are the new offsets when doing exact measurement for xD, zD and zTOFW for RUN 39-61
//double tof_offs[27] = {
//114.916,
//113.592,
//116.172,
//114.396,
//115.279,
//115.910,
//116.751,
//116.002,
//114.242,
//115.822,
//112.410,
//115.801,
//113.355,
//114.925,
//113.514,
//112.817,
//114.214,
//111.899,
//113.050,
//112.578,
//116.283,
//116.685,
//115.884,
//115.922,
//114.642,
//114.207,
//114.212
//};


//offsets, I should remove them, not used...
double offsets[29] = {0,-2.31775e+00,-8.63789e-01,-1.28369e+00,2.02607e+00,6.30893e-01,5.29162e-01,-8.16594e-01,1.69442e-01,-3.52510e+00,-2.16979e+00,1.66105e+00,1.26320e+00,-2.26356e+00,-4.42463e+00,-6.51834e+00,-1.42568e-01,2.20837e+00,6.15364e+00,4.49112e-01,-4.82568e-01,-1.27471e+00,-3.10422e+00,-1.06742e+00,-2.47400e+00,-5.09690e-01,2.38830e-01,-7.83351e-01,0};

//Definition of histos-----------------
TH1F* h1_z_all_iso;
TH1F* h1_a_q_z5;
TH1F* h1_a_q_z6;

TH2F* h2_theta_out_vs_mw3_10b;
TH2F* h2_theta_out_vs_mw3_11b;
TH2F* h2_theta_out_vs_mw3_11c;
TH2F* h2_theta_out_vs_mw3_12c;

TH2F* h2_z_vs_a_q;
TH1F* h1_gamma_energyE_10B;
TH1F* h1_gamma_energyE_11B;
TH1F* h1_gamma_energyE_11B_no_border;
TH1F* h1_gamma_energyE_max_val_11B;
TH2F* h2_gamma_fst_vs_snd_11B;
TH1F* h1_theta_1_plus_theta_2_CALIFA_11b;
TH1F* h1_theta_1_plus_theta_2_CALIFA_10b;
TH2F* h2_E1_vs_E2_CALIFA_10b;
TH2F* h2_E1_vs_E2_CALIFA_11b;

TH1F* h1_E1_plus_E2_CALIFA_10b;
TH1F* h1_E1_plus_E2_CALIFA_11b;

TH2F* h2_long_mom_p1p2_long_mom11b;
TH2F* h2_long_mom_p1p2_long_mom10b;
TH1F* h1_transv_mom_difference_p1_p2_10b;
TH1F* h1_transv_mom_difference_p1_p2_11b;
TH2F* h2_tof_vs_aq_fix_g_11b;
TH2F* h2_tof_vs_aq_fix_g_10b;
TH2F* h2_tof_vs_aq_fix_g_11c;
TH2F* h2_tof_vs_aq_fix_g_12c;
TH2F* h2_E1_plus_E2_CALIFA_vs_full_mom11b;

TH2F* h2_E1_vs_E2_CALIFA_11b_cut;
TH2F* h2_theta_1_vs_theta_2_CALIFA_11b_cut;
TH1F* h1_binding_energy_11b;
TH1F* h1_cos_gamma_11b_p_i;
TH1F* h1_time_diff_gamma_11b;

TH2F* h2_theta1_theta2_11b;
TH2F* h2_psi1_psi2_11b;
TH2F* h2_psi1_psi2_11b_2pi;
TH1F* h1_tof_detnr_strange_12c;
TH1F* h1_tof_detnr_strange_12c_large_aq;
TH1F* h1_mw3_fy_12c;
TH1F* h1_mw3_fy_strange_12c;
TH1F* h1_mw3_fy_strange_12c_large_aq;
TH1F* h1_tof_detnr_12c;
TH2F* h2_z_vs_a_q_nocut;
TH1F* h1_theta_1_CALIFA_11b;
TH1F* h1_theta_2_CALIFA_11b;

TH2F* h2_gamma_energyE_psi_11B;
TH2F* h2_long_mom_p1p2_long_mom11b_no_border;
TH2F* h2_long_mom_p1p2_long_mom11b_high;
TH2F* h2_long_mom_p1p2_long_mom11b_low;
TH2F* h2_long_mom_p1p2_long_mom11b_no_border_high;
TH2F* h2_long_mom_p1p2_long_mom11b_no_border_low;
TH2F* h2_long_mom_p1p2_long_mom11b_small_opening;
TH2F* h2_long_mom_p1p2_long_mom11b_theta_phi_constr;
TH2F* h2_long_mom_p1p2_long_mom11b_large_opening;
TH1F* h1_cos_gamma_11b_p_i_cms;
TH1F* h1_cos_gamma_11b_p_i_cms_optimized;
TH1F* h1_psi_in_11B_cut_triangle;
TH1F* h1_psi_in_11B;
TH1F* h1_gamma_energyE_max_val_11B_angle_cut;

TH2F* h2_tof_path_strange12c;
TH2F* h2_tof_path_strange_12c_large_aq;
TH2F* h2_tof_path_12c;

TH2F* h2_radius_path_strange12c;
TH2F* h2_radius_path_strange_12c_large_aq;
TH2F* h2_radius_path_12c;
TH2F* h2_mw2x_vs_tof_12c;
TH2F* h2_mw3x_vs_tof_12c;
TH2F* h2_charge_vs_tof_12c;
TH2F* h2_psi_in_vs_tof_12c;
TH2F* h2_detnr_vs_tof_12c;

//further analyis histograms for mw2,3 and psi_in..
TH1F* h1_mw3_fy_11b_cut_triangle;
TH1F* h1_mw2_fy_11b_cut_triangle;
TH1F* h1_mw3_fy_11b;
TH1F* h1_mw2_fy_11b;
TH1F* h1_psi_in_11b_cut_triangle;
TH1F* h1_psi_in_11b;

TH1F* h1_angl_diff_y;
TH1F* h1_angl_diff_y_cut_triangle;
//-------------------------------------
TH2F* h2_long_mom_p1p2_long_mom11b_low_optimized;
TH2F* h2_long_mom_p1p2_long_mom11b_high_optimized;
TH2F* h2_long_mom_p1p2_long_mom11b_low_optimized_08;
TH1F* h1_p_missing_11B;
TH1F* h1_cos_low_gamma;
TH2F* h2_theta_sum_vs_phi_diff_11b;
//HISTOS----------------------------------------------------------------------

sprintf(hist_name, "Z for all isotopes");
h1_z_all_iso = new TH1F(hist_name,hist_name,1000,1,10);
h1_z_all_iso->GetXaxis()->SetTitle("Z (charge)");
h1_z_all_iso->GetYaxis()->SetTitle("# Events ");
h1_z_all_iso->GetXaxis()->CenterTitle(true);
h1_z_all_iso->GetYaxis()->CenterTitle(true);
h1_z_all_iso->GetYaxis()->SetLabelSize(0.045);
h1_z_all_iso->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "A/q for Z = 5");
h1_a_q_z5 = new TH1F(hist_name,hist_name,400,1,3);
h1_a_q_z5->GetXaxis()->SetTitle("A/q for Z = 5");
h1_a_q_z5->GetYaxis()->SetTitle(" # Events");
h1_a_q_z5->GetXaxis()->CenterTitle(true);
h1_a_q_z5->GetYaxis()->CenterTitle(true);
h1_a_q_z5->GetYaxis()->SetLabelSize(0.045);
h1_a_q_z5->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "A/q for Z = 6");
h1_a_q_z6 = new TH1F(hist_name,hist_name,400,1,3);
h1_a_q_z6->GetXaxis()->SetTitle("A/q for Z = 6");
h1_a_q_z6->GetYaxis()->SetTitle(" # Events");
h1_a_q_z6->GetXaxis()->CenterTitle(true);
h1_a_q_z6->GetYaxis()->CenterTitle(true);
h1_a_q_z6->GetYaxis()->SetLabelSize(0.045);
h1_a_q_z6->GetYaxis()->SetTitleSize(0.045);

//- Plots for the third step to get psi_in parameters
sprintf(hist_name, "Theta_out versus MWPC3.fX for 12C");
h2_theta_out_vs_mw3_12c = new TH2F(hist_name,hist_name,800,-400,400,1200,0,0.6);
h2_theta_out_vs_mw3_12c->GetXaxis()->SetTitle("MWPC3 f.X [mm]");
h2_theta_out_vs_mw3_12c->GetYaxis()->SetTitle("Theta_out [degrees]");
h2_theta_out_vs_mw3_12c->GetXaxis()->CenterTitle(true);
h2_theta_out_vs_mw3_12c->GetYaxis()->CenterTitle(true);
h2_theta_out_vs_mw3_12c->GetYaxis()->SetLabelSize(0.045);
h2_theta_out_vs_mw3_12c->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Theta_out versus MWPC3.fX for 11C");
h2_theta_out_vs_mw3_11c = new TH2F(hist_name,hist_name,800,-400,400,1200,0,0.6);
h2_theta_out_vs_mw3_11c->GetXaxis()->SetTitle("MWPC3 f.X [mm]");
h2_theta_out_vs_mw3_11c->GetYaxis()->SetTitle("Theta_out [degrees]");
h2_theta_out_vs_mw3_11c->GetXaxis()->CenterTitle(true);
h2_theta_out_vs_mw3_11c->GetYaxis()->CenterTitle(true);
h2_theta_out_vs_mw3_11c->GetYaxis()->SetLabelSize(0.045);
h2_theta_out_vs_mw3_11c->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Theta_out versus MWPC3.fX for 11B");
h2_theta_out_vs_mw3_11b = new TH2F(hist_name,hist_name,800,-400,400,1200,0,0.6);
h2_theta_out_vs_mw3_11b->GetXaxis()->SetTitle("MWPC3 f.X [mm]");
h2_theta_out_vs_mw3_11b->GetYaxis()->SetTitle("Theta_out [degrees]");
h2_theta_out_vs_mw3_11b->GetXaxis()->CenterTitle(true);
h2_theta_out_vs_mw3_11b->GetYaxis()->CenterTitle(true);
h2_theta_out_vs_mw3_11b->GetYaxis()->SetLabelSize(0.045);
h2_theta_out_vs_mw3_11b->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Theta_out versus MWPC3.fX for 10B");
h2_theta_out_vs_mw3_10b = new TH2F(hist_name,hist_name,800,-400,400,1200,0,0.6);
h2_theta_out_vs_mw3_10b->GetXaxis()->SetTitle("MWPC3 f.X [mm]");
h2_theta_out_vs_mw3_10b->GetYaxis()->SetTitle("Theta_out [degrees]");
h2_theta_out_vs_mw3_10b->GetXaxis()->CenterTitle(true);
h2_theta_out_vs_mw3_10b->GetYaxis()->CenterTitle(true);
h2_theta_out_vs_mw3_10b->GetYaxis()->SetLabelSize(0.045);
h2_theta_out_vs_mw3_10b->GetYaxis()->SetTitleSize(0.045);
//
//first serious try:
sprintf(hist_name, "Z versus A/q");
h2_z_vs_a_q = new TH2F(hist_name,hist_name,200,1,3,100,0,10);
h2_z_vs_a_q->GetXaxis()->SetTitle("A/q");
h2_z_vs_a_q->GetYaxis()->SetTitle("Z (charge) ");
h2_z_vs_a_q->GetXaxis()->CenterTitle(true);
h2_z_vs_a_q->GetYaxis()->CenterTitle(true);
h2_z_vs_a_q->GetYaxis()->SetLabelSize(0.045);
h2_z_vs_a_q->GetYaxis()->SetTitleSize(0.045);

//gamma analysis from CALIFA

sprintf(hist_name, "Gamma Energies E for 10B");
h1_gamma_energyE_10B = new TH1F(hist_name,hist_name,200,0,10);
h1_gamma_energyE_10B->GetXaxis()->SetTitle("Energy E [MeV]");
h1_gamma_energyE_10B->GetYaxis()->SetTitle("Number of entries");
h1_gamma_energyE_10B->GetXaxis()->CenterTitle(true);
h1_gamma_energyE_10B->GetYaxis()->CenterTitle(true);
h1_gamma_energyE_10B->GetYaxis()->SetLabelSize(0.045);
h1_gamma_energyE_10B->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Gamma Energies E for 11B");
h1_gamma_energyE_11B = new TH1F(hist_name,hist_name,200,0,10);
h1_gamma_energyE_11B->GetXaxis()->SetTitle("Energy E [MeV]");
h1_gamma_energyE_11B->GetYaxis()->SetTitle("Number of entries");
h1_gamma_energyE_11B->GetXaxis()->CenterTitle(true);
h1_gamma_energyE_11B->GetYaxis()->CenterTitle(true);
h1_gamma_energyE_11B->GetYaxis()->SetLabelSize(0.045);
h1_gamma_energyE_11B->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Highest Gamma Energy E for 11B");
h1_gamma_energyE_max_val_11B = new TH1F(hist_name,hist_name,200,0,10);
h1_gamma_energyE_max_val_11B->GetXaxis()->SetTitle("Energy E [MeV]");
h1_gamma_energyE_max_val_11B->GetYaxis()->SetTitle("Number of entries");
h1_gamma_energyE_max_val_11B->GetXaxis()->CenterTitle(true);
h1_gamma_energyE_max_val_11B->GetYaxis()->CenterTitle(true);
h1_gamma_energyE_max_val_11B->GetYaxis()->SetLabelSize(0.045);
h1_gamma_energyE_max_val_11B->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Highest Gamma Energy E  vs second highest for 11B");
h2_gamma_fst_vs_snd_11B = new TH2F(hist_name,hist_name,200,0,10,200,0,10);
h2_gamma_fst_vs_snd_11B->GetXaxis()->SetTitle("Gamma 2nd highestEnergy E [MeV]");
h2_gamma_fst_vs_snd_11B->GetYaxis()->SetTitle("Gamma highest Energy E [MeV]");
h2_gamma_fst_vs_snd_11B->GetXaxis()->CenterTitle(true);
h2_gamma_fst_vs_snd_11B->GetYaxis()->CenterTitle(true);
h2_gamma_fst_vs_snd_11B->GetYaxis()->SetLabelSize(0.045);
h2_gamma_fst_vs_snd_11B->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Highest Gamma Energy E for 11B with theta and phi cut");
h1_gamma_energyE_max_val_11B_angle_cut = new TH1F(hist_name,hist_name,200,0,10);
h1_gamma_energyE_max_val_11B_angle_cut->GetXaxis()->SetTitle("Energy E [MeV]");
h1_gamma_energyE_max_val_11B_angle_cut->GetYaxis()->SetTitle("Number of entries");
h1_gamma_energyE_max_val_11B_angle_cut->GetXaxis()->CenterTitle(true);
h1_gamma_energyE_max_val_11B_angle_cut->GetYaxis()->CenterTitle(true);
h1_gamma_energyE_max_val_11B_angle_cut->GetYaxis()->SetLabelSize(0.045);
h1_gamma_energyE_max_val_11B_angle_cut->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "Gamma Energies E for 11B without border");
h1_gamma_energyE_11B_no_border = new TH1F(hist_name,hist_name,200,0,10);
h1_gamma_energyE_11B_no_border->GetXaxis()->SetTitle("Energy E [MeV]");
h1_gamma_energyE_11B_no_border->GetYaxis()->SetTitle("Number of entries");
h1_gamma_energyE_11B_no_border->GetXaxis()->CenterTitle(true);
h1_gamma_energyE_11B_no_border->GetYaxis()->CenterTitle(true);
h1_gamma_energyE_11B_no_border->GetYaxis()->SetLabelSize(0.045);
h1_gamma_energyE_11B_no_border->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Gamma Energies versus Psi for  11B");
h2_gamma_energyE_psi_11B = new TH2F(hist_name,hist_name,200,0,10,60,0,360);
h2_gamma_energyE_psi_11B->GetXaxis()->SetTitle("Energy E [MeV]");
h2_gamma_energyE_psi_11B->GetYaxis()->SetTitle("Psi [degrees]");
h2_gamma_energyE_psi_11B->GetXaxis()->CenterTitle(true);
h2_gamma_energyE_psi_11B->GetYaxis()->CenterTitle(true);
h2_gamma_energyE_psi_11B->GetYaxis()->SetLabelSize(0.045);
h2_gamma_energyE_psi_11B->GetYaxis()->SetTitleSize(0.045);

//check psi_in for all 11B and for the triangle...
sprintf(hist_name, "Theta_in for 11B");
h1_psi_in_11B = new TH1F(hist_name,hist_name,200,-0.1,0.1);
h1_psi_in_11B->GetXaxis()->SetTitle("Theta_in [rad]");
h1_psi_in_11B->GetYaxis()->SetTitle("Number of entries");
h1_psi_in_11B->GetXaxis()->CenterTitle(true);
h1_psi_in_11B->GetYaxis()->CenterTitle(true);
h1_psi_in_11B->GetYaxis()->SetLabelSize(0.045);
h1_psi_in_11B->GetYaxis()->SetTitleSize(0.045);

//cut triangle:
//check psi_in for all 11B and for the triangle...
sprintf(hist_name, "Theta_in for 11B in cut triangle");
h1_psi_in_11B_cut_triangle = new TH1F(hist_name,hist_name,200,-0.1,0.1);
h1_psi_in_11B_cut_triangle->GetXaxis()->SetTitle("Theta_in [rad]");
h1_psi_in_11B_cut_triangle->GetYaxis()->SetTitle("Number of entries");
h1_psi_in_11B_cut_triangle->GetXaxis()->CenterTitle(true);
h1_psi_in_11B_cut_triangle->GetYaxis()->CenterTitle(true);
h1_psi_in_11B_cut_triangle->GetYaxis()->SetLabelSize(0.045);
h1_psi_in_11B_cut_triangle->GetYaxis()->SetTitleSize(0.045);

//here I do some Mw3 and mw2 and psi_in analysis for the cut triangle and the rest and compare it, to see if the strange cut triangle comes from straggling in Y direction...

sprintf(hist_name, "MW3.fY for the cut triangle , 11B");
h1_mw3_fy_11b_cut_triangle = new TH1F(hist_name,hist_name,400,-400,400);
h1_mw3_fy_11b_cut_triangle->GetXaxis()->SetTitle("MW3.fY [mm]");
h1_mw3_fy_11b_cut_triangle->GetYaxis()->SetTitle("No. of Events");
h1_mw3_fy_11b_cut_triangle->GetXaxis()->CenterTitle(true);
h1_mw3_fy_11b_cut_triangle->GetYaxis()->CenterTitle(true);
h1_mw3_fy_11b_cut_triangle->GetYaxis()->SetLabelSize(0.045);
h1_mw3_fy_11b_cut_triangle->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "MW2.fY for the cut triangle , 11B");
h1_mw2_fy_11b_cut_triangle = new TH1F(hist_name,hist_name,400,-400,400);
h1_mw2_fy_11b_cut_triangle->GetXaxis()->SetTitle("MW2.fY [mm]");
h1_mw2_fy_11b_cut_triangle->GetYaxis()->SetTitle("No. of Events");
h1_mw2_fy_11b_cut_triangle->GetXaxis()->CenterTitle(true);
h1_mw2_fy_11b_cut_triangle->GetYaxis()->CenterTitle(true);
h1_mw2_fy_11b_cut_triangle->GetYaxis()->SetLabelSize(0.045);
h1_mw2_fy_11b_cut_triangle->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "MW3.fY for anti-correlated line , 11B");
h1_mw3_fy_11b = new TH1F(hist_name,hist_name,400,-400,400);
h1_mw3_fy_11b->GetXaxis()->SetTitle("MW3.fY [mm]");
h1_mw3_fy_11b->GetYaxis()->SetTitle("No. of Events");
h1_mw3_fy_11b->GetXaxis()->CenterTitle(true);
h1_mw3_fy_11b->GetYaxis()->CenterTitle(true);
h1_mw3_fy_11b->GetYaxis()->SetLabelSize(0.045);
h1_mw3_fy_11b->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "MW2.fY for anti-correlated line , 11B");
h1_mw2_fy_11b = new TH1F(hist_name,hist_name,400,-400,400);
h1_mw2_fy_11b->GetXaxis()->SetTitle("MW2.fY [mm]");
h1_mw2_fy_11b->GetYaxis()->SetTitle("No. of Events");
h1_mw2_fy_11b->GetXaxis()->CenterTitle(true);
h1_mw2_fy_11b->GetYaxis()->CenterTitle(true);
h1_mw2_fy_11b->GetYaxis()->SetLabelSize(0.045);
h1_mw2_fy_11b->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "Psi_in for the cut triangle , 11B");
h1_psi_in_11b_cut_triangle = new TH1F(hist_name,hist_name,2000,-0.1,0.1);
h1_psi_in_11b_cut_triangle->GetXaxis()->SetTitle("Psi_in  [degrees]");
h1_psi_in_11b_cut_triangle->GetYaxis()->SetTitle("No. of Events");
h1_psi_in_11b_cut_triangle->GetXaxis()->CenterTitle(true);
h1_psi_in_11b_cut_triangle->GetYaxis()->CenterTitle(true);
h1_psi_in_11b_cut_triangle->GetYaxis()->SetLabelSize(0.045);
h1_psi_in_11b_cut_triangle->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "Psi_in for anti-correlated line, 11B");
h1_psi_in_11b = new TH1F(hist_name,hist_name,2000,-0.1,0.1);
h1_psi_in_11b->GetXaxis()->SetTitle("Psi_in [degrees]");
h1_psi_in_11b->GetYaxis()->SetTitle("No. of Events");
h1_psi_in_11b->GetXaxis()->CenterTitle(true);
h1_psi_in_11b->GetYaxis()->CenterTitle(true);
h1_psi_in_11b->GetYaxis()->SetLabelSize(0.045);
h1_psi_in_11b->GetYaxis()->SetTitleSize(0.045);


//-------------------------------------------------

//more CALIFA
sprintf(hist_name, "Theta1 plus Theta2 for CALIFA 11B");
h1_theta_1_plus_theta_2_CALIFA_11b = new TH1F(hist_name,hist_name,52,22.15,152.15);
h1_theta_1_plus_theta_2_CALIFA_11b->GetXaxis()->SetTitle("Theta1 plus Theta2 [degrees]");
h1_theta_1_plus_theta_2_CALIFA_11b->GetYaxis()->SetTitle("No. of Events");
h1_theta_1_plus_theta_2_CALIFA_11b->GetXaxis()->CenterTitle(true);
h1_theta_1_plus_theta_2_CALIFA_11b->GetYaxis()->CenterTitle(true);
h1_theta_1_plus_theta_2_CALIFA_11b->GetYaxis()->SetLabelSize(0.045);
h1_theta_1_plus_theta_2_CALIFA_11b->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Theta1 plus Theta2 versus phi difference for 11B");
h2_theta_sum_vs_phi_diff_11b = new TH2F(hist_name,hist_name,52,22.15,152.15,60,0,360);
h2_theta_sum_vs_phi_diff_11b->GetXaxis()->SetTitle("Theta sum [degrees]");
h2_theta_sum_vs_phi_diff_11b->GetYaxis()->SetTitle(" Phi difference [degrees]");
h2_theta_sum_vs_phi_diff_11b->GetXaxis()->CenterTitle(true);
h2_theta_sum_vs_phi_diff_11b->GetYaxis()->CenterTitle(true);
h2_theta_sum_vs_phi_diff_11b->GetYaxis()->SetLabelSize(0.045);
h2_theta_sum_vs_phi_diff_11b->GetYaxis()->SetTitleSize(0.045);



//here to check binning of angles:
sprintf(hist_name, "Theta1 for  CALIFA 11B");
h1_theta_1_CALIFA_11b = new TH1F(hist_name,hist_name,1500,0,150);
h1_theta_1_CALIFA_11b->GetXaxis()->SetTitle("Theta1 [degrees]");
h1_theta_1_CALIFA_11b->GetYaxis()->SetTitle("No. of Events");
h1_theta_1_CALIFA_11b->GetXaxis()->CenterTitle(true);
h1_theta_1_CALIFA_11b->GetYaxis()->CenterTitle(true);
h1_theta_1_CALIFA_11b->GetYaxis()->SetLabelSize(0.045);
h1_theta_1_CALIFA_11b->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Theta2 for  CALIFA 11B");
h1_theta_2_CALIFA_11b = new TH1F(hist_name,hist_name,1500,0,150);
h1_theta_2_CALIFA_11b->GetXaxis()->SetTitle("Theta2 [degrees]");
h1_theta_2_CALIFA_11b->GetYaxis()->SetTitle("No. of Events");
h1_theta_2_CALIFA_11b->GetXaxis()->CenterTitle(true);
h1_theta_2_CALIFA_11b->GetYaxis()->CenterTitle(true);
h1_theta_2_CALIFA_11b->GetYaxis()->SetLabelSize(0.045);
h1_theta_2_CALIFA_11b->GetYaxis()->SetTitleSize(0.045);



//--------------------
sprintf(hist_name, "Theta1 plus Theta2 for CALIFA 10B");
h1_theta_1_plus_theta_2_CALIFA_10b = new TH1F(hist_name,hist_name,52,22.15,152.15);
h1_theta_1_plus_theta_2_CALIFA_10b->GetXaxis()->SetTitle("Theta1 plus Theta2 [degrees]");
h1_theta_1_plus_theta_2_CALIFA_10b->GetYaxis()->SetTitle("No. of Events");
h1_theta_1_plus_theta_2_CALIFA_10b->GetXaxis()->CenterTitle(true);
h1_theta_1_plus_theta_2_CALIFA_10b->GetYaxis()->CenterTitle(true);
h1_theta_1_plus_theta_2_CALIFA_10b->GetYaxis()->SetLabelSize(0.045);
h1_theta_1_plus_theta_2_CALIFA_10b->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "E1 versus E2 for CALIFA for 11B");
h2_E1_vs_E2_CALIFA_11b = new TH2F(hist_name,hist_name,50,0,500,50,0,500);
h2_E1_vs_E2_CALIFA_11b->GetXaxis()->SetTitle("E1 [MeV]");
h2_E1_vs_E2_CALIFA_11b->GetYaxis()->SetTitle("E2 [MeV]");
h2_E1_vs_E2_CALIFA_11b->GetXaxis()->CenterTitle(true);
h2_E1_vs_E2_CALIFA_11b->GetYaxis()->CenterTitle(true);
h2_E1_vs_E2_CALIFA_11b->GetYaxis()->SetLabelSize(0.045);
h2_E1_vs_E2_CALIFA_11b->GetYaxis()->SetTitleSize(0.045);



sprintf(hist_name, "E1 versus E2 for CALIFA for 10B");
h2_E1_vs_E2_CALIFA_10b = new TH2F(hist_name,hist_name,50,0,500,50,0,500);
h2_E1_vs_E2_CALIFA_10b->GetXaxis()->SetTitle("E1 [MeV]");
h2_E1_vs_E2_CALIFA_10b->GetYaxis()->SetTitle("E2 [MeV]");
h2_E1_vs_E2_CALIFA_10b->GetXaxis()->CenterTitle(true);
h2_E1_vs_E2_CALIFA_10b->GetYaxis()->CenterTitle(true);
h2_E1_vs_E2_CALIFA_10b->GetYaxis()->SetLabelSize(0.045);
h2_E1_vs_E2_CALIFA_10b->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "E1 plus E2 for CALIFA 10B");
h1_E1_plus_E2_CALIFA_10b = new TH1F(hist_name,hist_name,100,0,1000);
h1_E1_plus_E2_CALIFA_10b->GetXaxis()->SetTitle("E1 plus E2 [MeV]");
h1_E1_plus_E2_CALIFA_10b->GetYaxis()->SetTitle("No. of Events");
h1_E1_plus_E2_CALIFA_10b->GetXaxis()->CenterTitle(true);
h1_E1_plus_E2_CALIFA_10b->GetYaxis()->CenterTitle(true);
h1_E1_plus_E2_CALIFA_10b->GetYaxis()->SetLabelSize(0.045);
h1_E1_plus_E2_CALIFA_10b->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "E1 plus E2 for CALIFA 11B");
h1_E1_plus_E2_CALIFA_11b = new TH1F(hist_name,hist_name,100,0,1000);
h1_E1_plus_E2_CALIFA_11b->GetXaxis()->SetTitle("E1 plus E2 [MeV]");
h1_E1_plus_E2_CALIFA_11b->GetYaxis()->SetTitle("No. of Events");
h1_E1_plus_E2_CALIFA_11b->GetXaxis()->CenterTitle(true);
h1_E1_plus_E2_CALIFA_11b->GetYaxis()->CenterTitle(true);
h1_E1_plus_E2_CALIFA_11b->GetYaxis()->SetLabelSize(0.045);
h1_E1_plus_E2_CALIFA_11b->GetYaxis()->SetTitleSize(0.045);

if (gamma_given == 1.42942){

//these histos are for checking strange behavior of tof vs A/q at 12C:
sprintf(hist_name, "TOF Wall Detector Number for tof < 35.13 ns , 12C");
h1_tof_detnr_strange_12c = new TH1F(hist_name,hist_name,30,0,30);
h1_tof_detnr_strange_12c->GetXaxis()->SetTitle("DetectorNumber TOFW");
h1_tof_detnr_strange_12c->GetYaxis()->SetTitle("No. of Events");
h1_tof_detnr_strange_12c->GetXaxis()->CenterTitle(true);
h1_tof_detnr_strange_12c->GetYaxis()->CenterTitle(true);
h1_tof_detnr_strange_12c->GetYaxis()->SetLabelSize(0.045);
h1_tof_detnr_strange_12c->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "TOF Wall Detector Number for 35.35 < tof < 35.61 ns , 12C");
h1_tof_detnr_12c = new TH1F(hist_name,hist_name,30,0,30);
h1_tof_detnr_12c->GetXaxis()->SetTitle("DetectorNumber TOFW");
h1_tof_detnr_12c->GetYaxis()->SetTitle("No. of Events");
h1_tof_detnr_12c->GetXaxis()->CenterTitle(true);
h1_tof_detnr_12c->GetYaxis()->CenterTitle(true);
h1_tof_detnr_12c->GetYaxis()->SetLabelSize(0.045);
h1_tof_detnr_12c->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "TOF Wall Detector Number for A/q > 1.997, 12C");
h1_tof_detnr_strange_12c_large_aq = new TH1F(hist_name,hist_name,30,0,30);
h1_tof_detnr_strange_12c_large_aq->GetXaxis()->SetTitle("DetectorNumber TOFW");
h1_tof_detnr_strange_12c_large_aq->GetYaxis()->SetTitle("No. of Events");
h1_tof_detnr_strange_12c_large_aq->GetXaxis()->CenterTitle(true);
h1_tof_detnr_strange_12c_large_aq->GetYaxis()->CenterTitle(true);
h1_tof_detnr_strange_12c_large_aq->GetYaxis()->SetLabelSize(0.045);
h1_tof_detnr_strange_12c_large_aq->GetYaxis()->SetTitleSize(0.045);

//Y inspection in MW3
sprintf(hist_name, "MW3.fY for tof < 35.13 ns , 12C");
h1_mw3_fy_strange_12c = new TH1F(hist_name,hist_name,400,-400,400);
h1_mw3_fy_strange_12c->GetXaxis()->SetTitle("MW3.fY [mm]");
h1_mw3_fy_strange_12c->GetYaxis()->SetTitle("No. of Events");
h1_mw3_fy_strange_12c->GetXaxis()->CenterTitle(true);
h1_mw3_fy_strange_12c->GetYaxis()->CenterTitle(true);
h1_mw3_fy_strange_12c->GetYaxis()->SetLabelSize(0.045);
h1_mw3_fy_strange_12c->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "MW3.fY for 35.35 < tof < 35.61 ns , 12C");
h1_mw3_fy_12c = new TH1F(hist_name,hist_name,400,-400,400);
h1_mw3_fy_12c->GetXaxis()->SetTitle("MW3.fY [mm]");
h1_mw3_fy_12c->GetYaxis()->SetTitle("No. of Events");
h1_mw3_fy_12c->GetXaxis()->CenterTitle(true);
h1_mw3_fy_12c->GetYaxis()->CenterTitle(true);
h1_mw3_fy_12c->GetYaxis()->SetLabelSize(0.045);
h1_mw3_fy_12c->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "MW3.fY for  A/q > 1.997, 12C");
h1_mw3_fy_strange_12c_large_aq = new TH1F(hist_name,hist_name,400,-400,400);
h1_mw3_fy_strange_12c_large_aq->GetXaxis()->SetTitle("MW3.fY [mm]");
h1_mw3_fy_strange_12c_large_aq->GetYaxis()->SetTitle("No. of Events");
h1_mw3_fy_strange_12c_large_aq->GetXaxis()->CenterTitle(true);
h1_mw3_fy_strange_12c_large_aq->GetYaxis()->CenterTitle(true);
h1_mw3_fy_strange_12c_large_aq->GetYaxis()->SetLabelSize(0.045);
h1_mw3_fy_strange_12c_large_aq->GetYaxis()->SetTitleSize(0.045);

//check tof vs pathlength...
sprintf(hist_name, "Time of flight vs pathlength, for tof < 35.13 ns , 12C");
h2_tof_path_strange12c = new TH2F(hist_name,hist_name,400,7300,7700,1000,30,40);
h2_tof_path_strange12c->GetXaxis()->SetTitle("Pathlength [mm]");
h2_tof_path_strange12c->GetYaxis()->SetTitle("ToF [ns]");
h2_tof_path_strange12c->GetXaxis()->CenterTitle(true);
h2_tof_path_strange12c->GetYaxis()->CenterTitle(true);
h2_tof_path_strange12c->GetYaxis()->SetLabelSize(0.045);
h2_tof_path_strange12c->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "Time of flight vs pathlength, for A/q > 1.997, 12C");
h2_tof_path_strange_12c_large_aq = new TH2F(hist_name,hist_name,400,7300,7700,1000,30,40);
h2_tof_path_strange_12c_large_aq->GetXaxis()->SetTitle("Pathlength [mm]");
h2_tof_path_strange_12c_large_aq->GetYaxis()->SetTitle("ToF [ns]");
h2_tof_path_strange_12c_large_aq->GetXaxis()->CenterTitle(true);
h2_tof_path_strange_12c_large_aq->GetYaxis()->CenterTitle(true);
h2_tof_path_strange_12c_large_aq->GetYaxis()->SetLabelSize(0.045);
h2_tof_path_strange_12c_large_aq->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "Time of flight vs pathlength, for 35.35 < tof < 35.61 ns , 12C");
h2_tof_path_12c = new TH2F(hist_name,hist_name,400,7300,7700,1000,30,40);
h2_tof_path_12c->GetXaxis()->SetTitle("Pathlength [mm]");
h2_tof_path_12c->GetYaxis()->SetTitle("ToF [ns]");
h2_tof_path_12c->GetXaxis()->CenterTitle(true);
h2_tof_path_12c->GetYaxis()->CenterTitle(true);
h2_tof_path_12c->GetYaxis()->SetLabelSize(0.045);
h2_tof_path_12c->GetYaxis()->SetTitleSize(0.045);

//check radius vs time of flight....
sprintf(hist_name, "Radius vs time of flight, for tof < 35.13 ns , 12C");
h2_radius_path_strange12c = new TH2F(hist_name,hist_name,400,6300,6700,1000,30,40);
h2_radius_path_strange12c->GetXaxis()->SetTitle("Radius [mm]");
h2_radius_path_strange12c->GetYaxis()->SetTitle("ToF [ns]");
h2_radius_path_strange12c->GetXaxis()->CenterTitle(true);
h2_radius_path_strange12c->GetYaxis()->CenterTitle(true);
h2_radius_path_strange12c->GetYaxis()->SetLabelSize(0.045);
h2_radius_path_strange12c->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Radius vs time of flight, for A/q > 1.997, 12C");
h2_radius_path_strange_12c_large_aq = new TH2F(hist_name,hist_name,400,6300,6700,1000,30,40);
h2_radius_path_strange_12c_large_aq->GetXaxis()->SetTitle("Radius [mm]");
h2_radius_path_strange_12c_large_aq->GetYaxis()->SetTitle("ToF [ns]");
h2_radius_path_strange_12c_large_aq->GetXaxis()->CenterTitle(true);
h2_radius_path_strange_12c_large_aq->GetYaxis()->CenterTitle(true);
h2_radius_path_strange_12c_large_aq->GetYaxis()->SetLabelSize(0.045);
h2_radius_path_strange_12c_large_aq->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Radius vs time of flight, for 35.35 < tof < 35.61 ns , 12C");
h2_radius_path_12c = new TH2F(hist_name,hist_name,400,6300,6700,1000,30,40);
h2_radius_path_12c->GetXaxis()->SetTitle("Radius [mm]");
h2_radius_path_12c->GetYaxis()->SetTitle("ToF [ns]");
h2_radius_path_12c->GetXaxis()->CenterTitle(true);
h2_radius_path_12c->GetYaxis()->CenterTitle(true);
h2_radius_path_12c->GetYaxis()->SetLabelSize(0.045);
h2_radius_path_12c->GetYaxis()->SetTitleSize(0.045);

//check mw2.x
sprintf(hist_name, "MW2.X vs time of flight for 12C");
h2_mw2x_vs_tof_12c = new TH2F(hist_name,hist_name,400,-400,400,1000,30,40);
h2_mw2x_vs_tof_12c->GetXaxis()->SetTitle("MW2.fX [mm]");
h2_mw2x_vs_tof_12c->GetYaxis()->SetTitle("ToF [ns]");
h2_mw2x_vs_tof_12c->GetXaxis()->CenterTitle(true);
h2_mw2x_vs_tof_12c->GetYaxis()->CenterTitle(true);
h2_mw2x_vs_tof_12c->GetYaxis()->SetLabelSize(0.045);
h2_mw2x_vs_tof_12c->GetYaxis()->SetTitleSize(0.045);

//check mw3.x
sprintf(hist_name, "MW3.X vs time of flight for 12C");
h2_mw3x_vs_tof_12c = new TH2F(hist_name,hist_name,400,-400,400,1000,30,40);
h2_mw3x_vs_tof_12c->GetXaxis()->SetTitle("MW3.fX [mm]");
h2_mw3x_vs_tof_12c->GetYaxis()->SetTitle("ToF [ns]");
h2_mw3x_vs_tof_12c->GetXaxis()->CenterTitle(true);
h2_mw3x_vs_tof_12c->GetYaxis()->CenterTitle(true);
h2_mw3x_vs_tof_12c->GetYaxis()->SetLabelSize(0.045);
h2_mw3x_vs_tof_12c->GetYaxis()->SetTitleSize(0.045);

//check charge...
sprintf(hist_name, "Charge vs time of flight for 12C");
h2_charge_vs_tof_12c = new TH2F(hist_name,hist_name,100,0,10,1000,30,40);
h2_charge_vs_tof_12c->GetXaxis()->SetTitle("Charge");
h2_charge_vs_tof_12c->GetYaxis()->SetTitle("ToF [ns]");
h2_charge_vs_tof_12c->GetXaxis()->CenterTitle(true);
h2_charge_vs_tof_12c->GetYaxis()->CenterTitle(true);
h2_charge_vs_tof_12c->GetYaxis()->SetLabelSize(0.045);
h2_charge_vs_tof_12c->GetYaxis()->SetTitleSize(0.045);

//check psi_in
sprintf(hist_name, "Psi_in vs time of flight for 12C");
h2_psi_in_vs_tof_12c = new TH2F(hist_name,hist_name,200,-0.1,0.1,1000,30,40);
h2_psi_in_vs_tof_12c->GetXaxis()->SetTitle("Psi_in [degrees]");
h2_psi_in_vs_tof_12c->GetYaxis()->SetTitle("ToF [ns]");
h2_psi_in_vs_tof_12c->GetXaxis()->CenterTitle(true);
h2_psi_in_vs_tof_12c->GetYaxis()->CenterTitle(true);
h2_psi_in_vs_tof_12c->GetYaxis()->SetLabelSize(0.045);
h2_psi_in_vs_tof_12c->GetYaxis()->SetTitleSize(0.045);


//check detnr...
sprintf(hist_name, "Detnr TOFW vs time of flight for 12C");
h2_detnr_vs_tof_12c = new TH2F(hist_name,hist_name,30,0,30,1000,30,40);
h2_detnr_vs_tof_12c->GetXaxis()->SetTitle("Detnr TOFW");
h2_detnr_vs_tof_12c->GetYaxis()->SetTitle("ToF [ns]");
h2_detnr_vs_tof_12c->GetXaxis()->CenterTitle(true);
h2_detnr_vs_tof_12c->GetYaxis()->CenterTitle(true);
h2_detnr_vs_tof_12c->GetYaxis()->SetLabelSize(0.045);
h2_detnr_vs_tof_12c->GetYaxis()->SetTitleSize(0.045);



//

sprintf(hist_name, "Z versus A/q all isotopes");
h2_z_vs_a_q_nocut = new TH2F(hist_name,hist_name,200,1,3,100,0,10);
h2_z_vs_a_q_nocut->GetXaxis()->SetTitle("A/q");
h2_z_vs_a_q_nocut->GetYaxis()->SetTitle("Z (charge) ");
h2_z_vs_a_q_nocut->GetXaxis()->CenterTitle(true);
h2_z_vs_a_q_nocut->GetYaxis()->CenterTitle(true);
h2_z_vs_a_q_nocut->GetYaxis()->SetLabelSize(0.045);
h2_z_vs_a_q_nocut->GetYaxis()->SetTitleSize(0.045);


//

sprintf(hist_name, "Theta1 vs Theta2 in CALIFA");
h2_theta1_theta2_11b = new TH2F(hist_name,hist_name,25,22.15,84.65,25,22.15,84.65);
h2_theta1_theta2_11b->GetXaxis()->SetTitle("Theta1 [degrees]");
h2_theta1_theta2_11b->GetYaxis()->SetTitle("Theta2 [degrees]");
h2_theta1_theta2_11b->GetXaxis()->CenterTitle(true);
h2_theta1_theta2_11b->GetYaxis()->CenterTitle(true);
h2_theta1_theta2_11b->GetYaxis()->SetLabelSize(0.045);
h2_theta1_theta2_11b->GetYaxis()->SetTitleSize(0.045);

//


sprintf(hist_name, "Phi1 vs Phi2 in CALIFA");
h2_psi1_psi2_11b = new TH2F(hist_name,hist_name,60,0,360,60,0,360);
h2_psi1_psi2_11b->GetXaxis()->SetTitle("Psi1 [degrees]");
h2_psi1_psi2_11b->GetYaxis()->SetTitle("Psi2 [degrees]");
h2_psi1_psi2_11b->GetXaxis()->CenterTitle(true);
h2_psi1_psi2_11b->GetYaxis()->CenterTitle(true);
h2_psi1_psi2_11b->GetYaxis()->SetLabelSize(0.045);
h2_psi1_psi2_11b->GetYaxis()->SetTitleSize(0.045);


//2pi test for phi...

sprintf(hist_name, "Phi1 vs Phi2 in CALIFA with 2Pi");
h2_psi1_psi2_11b_2pi = new TH2F(hist_name,hist_name,60,0,360,60,0,360);
h2_psi1_psi2_11b_2pi->GetXaxis()->SetTitle("Psi1 [degrees]");
h2_psi1_psi2_11b_2pi->GetYaxis()->SetTitle("Psi2 [degrees]");
h2_psi1_psi2_11b_2pi->GetXaxis()->CenterTitle(true);
h2_psi1_psi2_11b_2pi->GetYaxis()->CenterTitle(true);
h2_psi1_psi2_11b_2pi->GetYaxis()->SetLabelSize(0.045);
h2_psi1_psi2_11b_2pi->GetYaxis()->SetTitleSize(0.045);
//


sprintf(hist_name, "cos(gamma) in the z-x plane for 11B and p_i");
h1_cos_gamma_11b_p_i = new TH1F(hist_name,hist_name,40,-2,2);
h1_cos_gamma_11b_p_i->GetXaxis()->SetTitle("cos(gamma)");
h1_cos_gamma_11b_p_i->GetYaxis()->SetTitle("No. of Events");
h1_cos_gamma_11b_p_i->GetXaxis()->CenterTitle(true);
h1_cos_gamma_11b_p_i->GetYaxis()->CenterTitle(true);
h1_cos_gamma_11b_p_i->GetYaxis()->SetLabelSize(0.045);
h1_cos_gamma_11b_p_i->GetYaxis()->SetTitleSize(0.045);


//angle between p_in and p_11b in 12C frame...
sprintf(hist_name, "cos(gamma) in the z-x plane for 11B and p_i in 12C rest frame");
h1_cos_gamma_11b_p_i_cms = new TH1F(hist_name,hist_name,50,-1,1);
h1_cos_gamma_11b_p_i_cms->GetXaxis()->SetTitle("cos(gamma)");
h1_cos_gamma_11b_p_i_cms->GetYaxis()->SetTitle("No. of Events");
h1_cos_gamma_11b_p_i_cms->GetXaxis()->CenterTitle(true);
h1_cos_gamma_11b_p_i_cms->GetYaxis()->CenterTitle(true);
h1_cos_gamma_11b_p_i_cms->GetYaxis()->SetLabelSize(0.045);
h1_cos_gamma_11b_p_i_cms->GetYaxis()->SetTitleSize(0.045);

//here the optimized version of h1_cos_gamma_11b_p_i_cms sweeping both phi1 and ph2 +- 2.8 degrees...
sprintf(hist_name, "cos(gamma) in the z-x plane for 11B and p_i in 12C rest frame  optimized version");
h1_cos_gamma_11b_p_i_cms_optimized = new TH1F(hist_name,hist_name,50,-1,1);
h1_cos_gamma_11b_p_i_cms_optimized->GetXaxis()->SetTitle("cos(gamma)");
h1_cos_gamma_11b_p_i_cms_optimized->GetYaxis()->SetTitle("No. of Events");
h1_cos_gamma_11b_p_i_cms_optimized->GetXaxis()->CenterTitle(true);
h1_cos_gamma_11b_p_i_cms_optimized->GetYaxis()->CenterTitle(true);
h1_cos_gamma_11b_p_i_cms_optimized->GetYaxis()->SetLabelSize(0.045);
h1_cos_gamma_11b_p_i_cms_optimized->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Time difference between START and CALIFA for the gammas");
h1_time_diff_gamma_11b = new TH1F(hist_name,hist_name,200,-100,100);
h1_time_diff_gamma_11b->GetXaxis()->SetTitle("Time difference [ns]");
h1_time_diff_gamma_11b->GetYaxis()->SetTitle("No. of Events");
h1_time_diff_gamma_11b->GetXaxis()->CenterTitle(true);
h1_time_diff_gamma_11b->GetYaxis()->CenterTitle(true);
h1_time_diff_gamma_11b->GetYaxis()->SetLabelSize(0.045);
h1_time_diff_gamma_11b->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "Bindungsenergie B_p = T_0 - (T_1+T_2+T_A-1)");
h1_binding_energy_11b = new TH1F(hist_name,hist_name,1000,-500,500);
h1_binding_energy_11b->GetXaxis()->SetTitle("Calc. Binding energy [MeV]");
h1_binding_energy_11b->GetYaxis()->SetTitle("No. of Events");
h1_binding_energy_11b->GetXaxis()->CenterTitle(true);
h1_binding_energy_11b->GetYaxis()->CenterTitle(true);
h1_binding_energy_11b->GetYaxis()->SetLabelSize(0.045);
h1_binding_energy_11b->GetYaxis()->SetTitleSize(0.045);
//here I look at the triangle cut, see if these events have special angles, Energies,...


sprintf(hist_name, "E1 versus E2 for CALIFA for 11B in the cut triangle");
h2_E1_vs_E2_CALIFA_11b_cut = new TH2F(hist_name,hist_name,50,0,500,50,0,500);
h2_E1_vs_E2_CALIFA_11b_cut->GetXaxis()->SetTitle("E1 [MeV]");
h2_E1_vs_E2_CALIFA_11b_cut->GetYaxis()->SetTitle("E2 [MeV]");
h2_E1_vs_E2_CALIFA_11b_cut->GetXaxis()->CenterTitle(true);
h2_E1_vs_E2_CALIFA_11b_cut->GetYaxis()->CenterTitle(true);
h2_E1_vs_E2_CALIFA_11b_cut->GetYaxis()->SetLabelSize(0.045);
h2_E1_vs_E2_CALIFA_11b_cut->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Theta1 vs Theta2 for CALIFA 11B in the cut triangle");
h2_theta_1_vs_theta_2_CALIFA_11b_cut = new TH2F(hist_name,hist_name,75,0,150,75,0,150);
h2_theta_1_vs_theta_2_CALIFA_11b_cut->GetXaxis()->SetTitle("Theta1 [degrees]");
h2_theta_1_vs_theta_2_CALIFA_11b_cut->GetYaxis()->SetTitle("Theta2 [degrees]");
h2_theta_1_vs_theta_2_CALIFA_11b_cut->GetXaxis()->CenterTitle(true);
h2_theta_1_vs_theta_2_CALIFA_11b_cut->GetYaxis()->CenterTitle(true);
h2_theta_1_vs_theta_2_CALIFA_11b_cut->GetYaxis()->SetLabelSize(0.045);
h2_theta_1_vs_theta_2_CALIFA_11b_cut->GetYaxis()->SetTitleSize(0.045);

//-----------------------------------------------------------------------------------------

sprintf(hist_name, "E1+E2 versus full momentum 11B");
h2_E1_plus_E2_CALIFA_vs_full_mom11b = new TH2F(hist_name,hist_name,100,0,1000,200,9500,11500);
h2_E1_plus_E2_CALIFA_vs_full_mom11b->GetXaxis()->SetTitle("E1+E2 [MeV]");
h2_E1_plus_E2_CALIFA_vs_full_mom11b->GetYaxis()->SetTitle("Full Momentum 11B");
h2_E1_plus_E2_CALIFA_vs_full_mom11b->GetXaxis()->CenterTitle(true);
h2_E1_plus_E2_CALIFA_vs_full_mom11b->GetYaxis()->CenterTitle(true);
h2_E1_plus_E2_CALIFA_vs_full_mom11b->GetYaxis()->SetLabelSize(0.045);
h2_E1_plus_E2_CALIFA_vs_full_mom11b->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Long.  mom. p1p2 vs long. mom. 11B");
h2_long_mom_p1p2_long_mom11b = new TH2F(hist_name,hist_name,200,0,2000,200,9500,11500);
h2_long_mom_p1p2_long_mom11b->GetXaxis()->SetTitle("Longitudinal mom. p1+p2");
h2_long_mom_p1p2_long_mom11b->GetYaxis()->SetTitle("Longitudinal mom. 11B");
h2_long_mom_p1p2_long_mom11b->GetXaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b->GetYaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b->GetYaxis()->SetLabelSize(0.045);
h2_long_mom_p1p2_long_mom11b->GetYaxis()->SetTitleSize(0.045);

//here a cut on opening angle > < 84 degreees..
sprintf(hist_name, "Long.  mom. p1p2 vs long. mom. 11B for opening angles smaller than 84 degrees");
h2_long_mom_p1p2_long_mom11b_small_opening = new TH2F(hist_name,hist_name,200,0,2000,200,9500,11500);
h2_long_mom_p1p2_long_mom11b_small_opening->GetXaxis()->SetTitle("Longitudinal mom. p1+p2");
h2_long_mom_p1p2_long_mom11b_small_opening->GetYaxis()->SetTitle("Longitudinal mom. 11B");
h2_long_mom_p1p2_long_mom11b_small_opening->GetXaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_small_opening->GetYaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_small_opening->GetYaxis()->SetLabelSize(0.045);
h2_long_mom_p1p2_long_mom11b_small_opening->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Long.  mom. p1p2 vs long. mom. 11B for opening angles larger than 84 degrees");
h2_long_mom_p1p2_long_mom11b_large_opening = new TH2F(hist_name,hist_name,200,0,2000,200,9500,11500);
h2_long_mom_p1p2_long_mom11b_large_opening->GetXaxis()->SetTitle("Longitudinal mom. p1+p2");
h2_long_mom_p1p2_long_mom11b_large_opening->GetYaxis()->SetTitle("Longitudinal mom. 11B");
h2_long_mom_p1p2_long_mom11b_large_opening->GetXaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_large_opening->GetYaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_large_opening->GetYaxis()->SetLabelSize(0.045);
h2_long_mom_p1p2_long_mom11b_large_opening->GetYaxis()->SetTitleSize(0.045);

//constraints both on theta and phi...
sprintf(hist_name, "Long.  mom. p1p2 vs long. mom. 11B for opening angles smaller than 84 degrees and phi-diff +-30degrees");
h2_long_mom_p1p2_long_mom11b_theta_phi_constr = new TH2F(hist_name,hist_name,200,0,2000,200,9500,11500);
h2_long_mom_p1p2_long_mom11b_theta_phi_constr->GetXaxis()->SetTitle("Longitudinal mom. p1+p2");
h2_long_mom_p1p2_long_mom11b_theta_phi_constr->GetYaxis()->SetTitle("Longitudinal mom. 11B");
h2_long_mom_p1p2_long_mom11b_theta_phi_constr->GetXaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_theta_phi_constr->GetYaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_theta_phi_constr->GetYaxis()->SetLabelSize(0.045);
h2_long_mom_p1p2_long_mom11b_theta_phi_constr->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "Long.  mom. p1p2 vs long. mom. 11B cos_high");
h2_long_mom_p1p2_long_mom11b_high = new TH2F(hist_name,hist_name,200,0,2000,200,9500,11500);
h2_long_mom_p1p2_long_mom11b_high->GetXaxis()->SetTitle("Longitudinal mom. p1+p2");
h2_long_mom_p1p2_long_mom11b_high->GetYaxis()->SetTitle("Longitudinal mom. 11B");
h2_long_mom_p1p2_long_mom11b_high->GetXaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_high->GetYaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_high->GetYaxis()->SetLabelSize(0.045);
h2_long_mom_p1p2_long_mom11b_high->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Long.  mom. p1p2 vs long. mom. 11B cos_low");
h2_long_mom_p1p2_long_mom11b_low = new TH2F(hist_name,hist_name,200,0,2000,200,9500,11500);
h2_long_mom_p1p2_long_mom11b_low->GetXaxis()->SetTitle("Longitudinal mom. p1+p2");
h2_long_mom_p1p2_long_mom11b_low->GetYaxis()->SetTitle("Longitudinal mom. 11B");
h2_long_mom_p1p2_long_mom11b_low->GetXaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_low->GetYaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_low->GetYaxis()->SetLabelSize(0.045);
h2_long_mom_p1p2_long_mom11b_low->GetYaxis()->SetTitleSize(0.045);

//do p1+p2 long.mom. vs 11b long. mom. plots with optimizations.... and restrictions on phi and theta....
sprintf(hist_name, "Long.  mom. p1p2 vs long. mom. 11B cos_low with optimized");
h2_long_mom_p1p2_long_mom11b_low_optimized = new TH2F(hist_name,hist_name,200,0,2000,200,9500,11500);
h2_long_mom_p1p2_long_mom11b_low_optimized->GetXaxis()->SetTitle("Longitudinal mom. p1+p2");
h2_long_mom_p1p2_long_mom11b_low_optimized->GetYaxis()->SetTitle("Longitudinal mom. 11B");
h2_long_mom_p1p2_long_mom11b_low_optimized->GetXaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_low_optimized->GetYaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_low_optimized->GetYaxis()->SetLabelSize(0.045);
h2_long_mom_p1p2_long_mom11b_low_optimized->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "Long.  mom. p1p2 vs long. mom. 11B with cos < -0.8");
h2_long_mom_p1p2_long_mom11b_low_optimized_08 = new TH2F(hist_name,hist_name,200,0,2000,200,9500,11500);
h2_long_mom_p1p2_long_mom11b_low_optimized_08->GetXaxis()->SetTitle("Longitudinal mom. p1+p2");
h2_long_mom_p1p2_long_mom11b_low_optimized_08->GetYaxis()->SetTitle("Longitudinal mom. 11B");
h2_long_mom_p1p2_long_mom11b_low_optimized_08->GetXaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_low_optimized_08->GetYaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_low_optimized_08->GetYaxis()->SetLabelSize(0.045);
h2_long_mom_p1p2_long_mom11b_low_optimized_08->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "Long.  mom. p1p2 vs long. mom. 11B cos_high with optimized");
h2_long_mom_p1p2_long_mom11b_high_optimized = new TH2F(hist_name,hist_name,200,0,2000,200,9500,11500);
h2_long_mom_p1p2_long_mom11b_high_optimized->GetXaxis()->SetTitle("Longitudinal mom. p1+p2");
h2_long_mom_p1p2_long_mom11b_high_optimized->GetYaxis()->SetTitle("Longitudinal mom. 11B");
h2_long_mom_p1p2_long_mom11b_high_optimized->GetXaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_high_optimized->GetYaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_high_optimized->GetYaxis()->SetLabelSize(0.045);
h2_long_mom_p1p2_long_mom11b_high_optimized->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "missing four momentum p_missing 11B");
h1_p_missing_11B = new TH1F(hist_name,hist_name,750,0,1500);
h1_p_missing_11B->GetXaxis()->SetTitle("p_missing [MeV/c]");
h1_p_missing_11B->GetYaxis()->SetTitle("No. of Events");
h1_p_missing_11B->GetXaxis()->CenterTitle(true);
h1_p_missing_11B->GetYaxis()->CenterTitle(true);
h1_p_missing_11B->GetYaxis()->SetLabelSize(0.045);
h1_p_missing_11B->GetYaxis()->SetTitleSize(0.045);


//------------------------------------------------------------------------------------------------

sprintf(hist_name, "Long.  mom. p1p2 vs long. mom. 11B no borders");
h2_long_mom_p1p2_long_mom11b_no_border = new TH2F(hist_name,hist_name,200,0,2000,200,9500,11500);
h2_long_mom_p1p2_long_mom11b_no_border->GetXaxis()->SetTitle("Longitudinal mom. p1+p2");
h2_long_mom_p1p2_long_mom11b_no_border->GetYaxis()->SetTitle("Longitudinal mom. 11B");
h2_long_mom_p1p2_long_mom11b_no_border->GetXaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_no_border->GetYaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_no_border->GetYaxis()->SetLabelSize(0.045);
h2_long_mom_p1p2_long_mom11b_no_border->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Long.  mom. p1p2 vs long. mom. 11B cos_high no borders");
h2_long_mom_p1p2_long_mom11b_no_border_high = new TH2F(hist_name,hist_name,200,0,2000,200,9500,11500);
h2_long_mom_p1p2_long_mom11b_no_border_high->GetXaxis()->SetTitle("Longitudinal mom. p1+p2");
h2_long_mom_p1p2_long_mom11b_no_border_high->GetYaxis()->SetTitle("Longitudinal mom. 11B");
h2_long_mom_p1p2_long_mom11b_no_border_high->GetXaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_no_border_high->GetYaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_no_border_high->GetYaxis()->SetLabelSize(0.045);
h2_long_mom_p1p2_long_mom11b_no_border_high->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Long.  mom. p1p2 vs long. mom. 11B cos_low no borders");
h2_long_mom_p1p2_long_mom11b_no_border_low = new TH2F(hist_name,hist_name,200,0,2000,200,9500,11500);
h2_long_mom_p1p2_long_mom11b_no_border_low->GetXaxis()->SetTitle("Longitudinal mom. p1+p2");
h2_long_mom_p1p2_long_mom11b_no_border_low->GetYaxis()->SetTitle("Longitudinal mom. 11B");
h2_long_mom_p1p2_long_mom11b_no_border_low->GetXaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_no_border_low->GetYaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b_no_border_low->GetYaxis()->SetLabelSize(0.045);
h2_long_mom_p1p2_long_mom11b_no_border_low->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "Long.  mom. p1p2 vs long. mom. 10B");
h2_long_mom_p1p2_long_mom10b = new TH2F(hist_name,hist_name,200,0,2000,200,8500,10500);
h2_long_mom_p1p2_long_mom10b->GetXaxis()->SetTitle("Longitudinal mom. p1+p2");
h2_long_mom_p1p2_long_mom10b->GetYaxis()->SetTitle("Longitudinal mom. 10B");
h2_long_mom_p1p2_long_mom10b->GetXaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom10b->GetYaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom10b->GetYaxis()->SetLabelSize(0.045);
h2_long_mom_p1p2_long_mom10b->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "B*rho/(beta*fixed_gamma) versus Time of Flight 11B");
h2_tof_vs_aq_fix_g_11b = new TH2F(hist_name,hist_name,900,1,3,1000,30,40);
h2_tof_vs_aq_fix_g_11b->GetXaxis()->SetTitle("B*rho/(beta*fixed_gamma)");
h2_tof_vs_aq_fix_g_11b->GetYaxis()->SetTitle("Time of Flight [ns]");
h2_tof_vs_aq_fix_g_11b->GetXaxis()->CenterTitle(true);
h2_tof_vs_aq_fix_g_11b->GetYaxis()->CenterTitle(true);
h2_tof_vs_aq_fix_g_11b->GetYaxis()->SetLabelSize(0.045);
h2_tof_vs_aq_fix_g_11b->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "B*rho/(beta*fixed_gamma) versus Time of Flight 10B");
h2_tof_vs_aq_fix_g_10b = new TH2F(hist_name,hist_name,900,1,3,1000,30,40);
h2_tof_vs_aq_fix_g_10b->GetXaxis()->SetTitle("B*rho/(beta*fixed_gamma)");
h2_tof_vs_aq_fix_g_10b->GetYaxis()->SetTitle("Time of Flight [ns]");
h2_tof_vs_aq_fix_g_10b->GetXaxis()->CenterTitle(true);
h2_tof_vs_aq_fix_g_10b->GetYaxis()->CenterTitle(true);
h2_tof_vs_aq_fix_g_10b->GetYaxis()->SetLabelSize(0.045);
h2_tof_vs_aq_fix_g_10b->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "B*rho/(beta*fixed_gamma) versus Time of Flight 11C");
h2_tof_vs_aq_fix_g_11c = new TH2F(hist_name,hist_name,900,1,3,1000,30,40);
h2_tof_vs_aq_fix_g_11c->GetXaxis()->SetTitle("B*rho/(beta*fixed_gamma)");
h2_tof_vs_aq_fix_g_11c->GetYaxis()->SetTitle("Time of Flight [ns]");
h2_tof_vs_aq_fix_g_11c->GetXaxis()->CenterTitle(true);
h2_tof_vs_aq_fix_g_11c->GetYaxis()->CenterTitle(true);
h2_tof_vs_aq_fix_g_11c->GetYaxis()->SetLabelSize(0.045);
h2_tof_vs_aq_fix_g_11c->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "B*rho/(beta*fixed_gamma) versus Time of Flight 12C");
h2_tof_vs_aq_fix_g_12c = new TH2F(hist_name,hist_name,900,1,3,1000,30,40);
h2_tof_vs_aq_fix_g_12c->GetXaxis()->SetTitle("B*rho/(beta*fixed_gamma)");
h2_tof_vs_aq_fix_g_12c->GetYaxis()->SetTitle("Time of Flight [ns]");
h2_tof_vs_aq_fix_g_12c->GetXaxis()->CenterTitle(true);
h2_tof_vs_aq_fix_g_12c->GetYaxis()->CenterTitle(true);
h2_tof_vs_aq_fix_g_12c->GetYaxis()->SetLabelSize(0.045);
h2_tof_vs_aq_fix_g_12c->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "Abs. Angle Difference before and after GLAD in y-Plane without cut triangle");
h1_angl_diff_y = new TH1F(hist_name,hist_name,100,0,25);
h1_angl_diff_y->GetXaxis()->SetTitle("Angle differnce (abs.) [degrees]");
h1_angl_diff_y->GetYaxis()->SetTitle("No. of Events");
h1_angl_diff_y->GetXaxis()->CenterTitle(true);
h1_angl_diff_y->GetYaxis()->CenterTitle(true);
h1_angl_diff_y->GetYaxis()->SetLabelSize(0.045);
h1_angl_diff_y->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Abs. Angle Difference before and after GLAD in y-Plane for cut triangle");
h1_angl_diff_y_cut_triangle = new TH1F(hist_name,hist_name,100,0,25);
h1_angl_diff_y_cut_triangle->GetXaxis()->SetTitle("Angle differnce (abs.) [degrees]");
h1_angl_diff_y_cut_triangle->GetYaxis()->SetTitle("No. of Events");
h1_angl_diff_y_cut_triangle->GetXaxis()->CenterTitle(true);
h1_angl_diff_y_cut_triangle->GetYaxis()->CenterTitle(true);
h1_angl_diff_y_cut_triangle->GetYaxis()->SetLabelSize(0.045);
h1_angl_diff_y_cut_triangle->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Gamma spectrum for 12C(p,2p)11B reaction  for cos < -0.6");
h1_cos_low_gamma = new TH1F(hist_name,hist_name,100,0,25);
h1_cos_low_gamma->GetXaxis()->SetTitle("Energy [MeV]");
h1_cos_low_gamma->GetYaxis()->SetTitle("No. of Events");
h1_cos_low_gamma->GetXaxis()->CenterTitle(true);
h1_cos_low_gamma->GetYaxis()->CenterTitle(true);
h1_cos_low_gamma->GetYaxis()->SetLabelSize(0.045);
h1_cos_low_gamma->GetYaxis()->SetTitleSize(0.045);

}
else {

sprintf(hist_name, "E1+E2 versus full momentum 11B");
h2_E1_plus_E2_CALIFA_vs_full_mom11b = new TH2F(hist_name,hist_name,100,0,1000,300,15000,18000);
h2_E1_plus_E2_CALIFA_vs_full_mom11b->GetXaxis()->SetTitle("E1+E2 [MeV]");
h2_E1_plus_E2_CALIFA_vs_full_mom11b->GetYaxis()->SetTitle("Full Momentum 11B");
h2_E1_plus_E2_CALIFA_vs_full_mom11b->GetXaxis()->CenterTitle(true);
h2_E1_plus_E2_CALIFA_vs_full_mom11b->GetYaxis()->CenterTitle(true);
h2_E1_plus_E2_CALIFA_vs_full_mom11b->GetYaxis()->SetLabelSize(0.045);
h2_E1_plus_E2_CALIFA_vs_full_mom11b->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "Long.  mom. p1p2 vs long. mom. 11B");
h2_long_mom_p1p2_long_mom11b = new TH2F(hist_name,hist_name,200,0,2000,300,15000,18000);
h2_long_mom_p1p2_long_mom11b->GetXaxis()->SetTitle("Longitudinal mom. p1+p2");
h2_long_mom_p1p2_long_mom11b->GetYaxis()->SetTitle("Longitudinal mom. 11B");
h2_long_mom_p1p2_long_mom11b->GetXaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b->GetYaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom11b->GetYaxis()->SetLabelSize(0.045);
h2_long_mom_p1p2_long_mom11b->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "Long.  mom. p1p2 vs long. mom. 10B");
h2_long_mom_p1p2_long_mom10b = new TH2F(hist_name,hist_name,200,0,2000,300,13500,16500);
h2_long_mom_p1p2_long_mom10b->GetXaxis()->SetTitle("Longitudinal mom. p1+p2");
h2_long_mom_p1p2_long_mom10b->GetYaxis()->SetTitle("Longitudinal mom. 10B");
h2_long_mom_p1p2_long_mom10b->GetXaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom10b->GetYaxis()->CenterTitle(true);
h2_long_mom_p1p2_long_mom10b->GetYaxis()->SetLabelSize(0.045);
h2_long_mom_p1p2_long_mom10b->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "B*rho/(beta*fixed_gamma) versus Time of Flight 11B");
h2_tof_vs_aq_fix_g_11b = new TH2F(hist_name,hist_name,900,1,3,1000,23,33);
h2_tof_vs_aq_fix_g_11b->GetXaxis()->SetTitle("B*rho/(beta*fixed_gamma)");
h2_tof_vs_aq_fix_g_11b->GetYaxis()->SetTitle("Time of Flight [ns]");
h2_tof_vs_aq_fix_g_11b->GetXaxis()->CenterTitle(true);
h2_tof_vs_aq_fix_g_11b->GetYaxis()->CenterTitle(true);
h2_tof_vs_aq_fix_g_11b->GetYaxis()->SetLabelSize(0.045);
h2_tof_vs_aq_fix_g_11b->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "B*rho/(beta*fixed_gamma) versus Time of Flight 10B");
h2_tof_vs_aq_fix_g_10b = new TH2F(hist_name,hist_name,900,1,3,1000,23,33);
h2_tof_vs_aq_fix_g_10b->GetXaxis()->SetTitle("B*rho/(beta*fixed_gamma)");
h2_tof_vs_aq_fix_g_10b->GetYaxis()->SetTitle("Time of Flight [ns]");
h2_tof_vs_aq_fix_g_10b->GetXaxis()->CenterTitle(true);
h2_tof_vs_aq_fix_g_10b->GetYaxis()->CenterTitle(true);
h2_tof_vs_aq_fix_g_10b->GetYaxis()->SetLabelSize(0.045);
h2_tof_vs_aq_fix_g_10b->GetYaxis()->SetTitleSize(0.045);

sprintf(hist_name, "B*rho/(beta*fixed_gamma) versus Time of Flight 11C");
h2_tof_vs_aq_fix_g_11c = new TH2F(hist_name,hist_name,900,1,3,1000,23,33);
h2_tof_vs_aq_fix_g_11c->GetXaxis()->SetTitle("B*rho/(beta*fixed_gamma)");
h2_tof_vs_aq_fix_g_11c->GetYaxis()->SetTitle("Time of Flight [ns]");
h2_tof_vs_aq_fix_g_11c->GetXaxis()->CenterTitle(true);
h2_tof_vs_aq_fix_g_11c->GetYaxis()->CenterTitle(true);
h2_tof_vs_aq_fix_g_11c->GetYaxis()->SetLabelSize(0.045);
h2_tof_vs_aq_fix_g_11c->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "B*rho/(beta*fixed_gamma) versus Time of Flight 12C");
h2_tof_vs_aq_fix_g_12c = new TH2F(hist_name,hist_name,900,1,3,1000,23,33);
h2_tof_vs_aq_fix_g_12c->GetXaxis()->SetTitle("B*rho/(beta*fixed_gamma)");
h2_tof_vs_aq_fix_g_12c->GetYaxis()->SetTitle("Time of Flight [ns]");
h2_tof_vs_aq_fix_g_12c->GetXaxis()->CenterTitle(true);
h2_tof_vs_aq_fix_g_12c->GetYaxis()->CenterTitle(true);
h2_tof_vs_aq_fix_g_12c->GetYaxis()->SetLabelSize(0.045);
h2_tof_vs_aq_fix_g_12c->GetYaxis()->SetTitleSize(0.045);

}

sprintf(hist_name, "Transv. mom. difference p1 p2 for 10B");
h1_transv_mom_difference_p1_p2_10b = new TH1F(hist_name,hist_name,100,0,500);
h1_transv_mom_difference_p1_p2_10b->GetXaxis()->SetTitle("Transv. mom. difference p1p2 [MeV/c]");
h1_transv_mom_difference_p1_p2_10b->GetYaxis()->SetTitle("Number of entries");
h1_transv_mom_difference_p1_p2_10b->GetXaxis()->CenterTitle(true);
h1_transv_mom_difference_p1_p2_10b->GetYaxis()->CenterTitle(true);
h1_transv_mom_difference_p1_p2_10b->GetYaxis()->SetLabelSize(0.045);
h1_transv_mom_difference_p1_p2_10b->GetYaxis()->SetTitleSize(0.045);


sprintf(hist_name, "Transv. mom. difference p1 p2 for 11B");
h1_transv_mom_difference_p1_p2_11b = new TH1F(hist_name,hist_name,100,0,500);
h1_transv_mom_difference_p1_p2_11b->GetXaxis()->SetTitle("Transv. mom. difference p1p2 [MeV/c]");
h1_transv_mom_difference_p1_p2_11b->GetYaxis()->SetTitle("Number of entries");
h1_transv_mom_difference_p1_p2_11b->GetXaxis()->CenterTitle(true);
h1_transv_mom_difference_p1_p2_11b->GetYaxis()->CenterTitle(true);
h1_transv_mom_difference_p1_p2_11b->GetYaxis()->SetLabelSize(0.045);
h1_transv_mom_difference_p1_p2_11b->GetYaxis()->SetTitleSize(0.045);


//----------------------------------------------------------------------------


//Path to input file
sprintf(fname,"/scratch8/ge37liw/workingspace/data/root_files/all_ts_unpack/sweep_target/ts_cone_cluster_%s_zero_cm.root",count_i);

TChain* chain = new TChain("evt");
chain->Reset();
chain->Add(fname);

//TClonesArrays containing the TObjects R3BSofToFWTcalData,... which allow access to data over function calls
//
TClonesArray* SofMwpc3HitData = new TClonesArray("R3BSofMwpcHitData",5);
R3BSofMwpcHitData** sofmwpc3hitdata;
TBranch *branchSofMwpc3HitData = chain->GetBranch("Mwpc3HitData");
branchSofMwpc3HitData->SetAddress(&SofMwpc3HitData);
//
TClonesArray* SofSciTcalData = new TClonesArray("R3BSofSciTcalData",2);
R3BSofSciTcalData** sofscitcaldata;
TBranch *branchSofSciTcalData = chain->GetBranch("SofSciTcalData");
branchSofSciTcalData->SetAddress(&SofSciTcalData);
//
TClonesArray* SofToFWTcalData = new TClonesArray("R3BSofToFWTcalData",2);
R3BSofToFWTcalData** softofwtcaldata;
TBranch *branchSofToFWTcalData = chain->GetBranch("SofToFWTcalData");
branchSofToFWTcalData->SetAddress(&SofToFWTcalData);
//
TClonesArray* SofMwpc0HitData = new TClonesArray("R3BSofMwpcHitData",5);
R3BSofMwpcHitData** sofmwpc0hitdata;
TBranch *branchSofMwpc0HitData = chain->GetBranch("Mwpc0HitData");
branchSofMwpc0HitData->SetAddress(&SofMwpc0HitData);
//
TClonesArray* SofMwpc1HitData = new TClonesArray("R3BSofMwpcHitData",5);
R3BSofMwpcHitData** sofmwpc1hitdata;
TBranch *branchSofMwpc1HitData = chain->GetBranch("Mwpc1HitData");
branchSofMwpc1HitData->SetAddress(&SofMwpc1HitData);
//
TClonesArray* SofMwpc2HitData = new TClonesArray("R3BSofMwpcHitData",5);
R3BSofMwpcHitData** sofmwpc2hitdata;
TBranch *branchSofMwpc2HitData = chain->GetBranch("Mwpc2HitData");
branchSofMwpc2HitData->SetAddress(&SofMwpc2HitData);
//
TClonesArray* SofToFWMappedData = new TClonesArray("R3BSofToFWMappedData",2);
R3BSofToFWMappedData** softofwmappeddata;
TBranch *branchSofToFWMappedData = chain->GetBranch("SofToFWMappedData");
branchSofToFWMappedData->SetAddress(&SofToFWMappedData);
//
TClonesArray* SofTwimHitData = new TClonesArray("R3BSofTwimHitData",2);
R3BSofTwimHitData** softwimhitdata;
TBranch *branchSofTwimHitData = chain->GetBranch("TwimHitData");
branchSofTwimHitData->SetAddress(&SofTwimHitData);
//
TClonesArray* CalifaHitData = new TClonesArray("R3BCalifaHitData",3);
R3BCalifaHitData** califahitdata;
TBranch *branchCalifaHitData = chain->GetBranch("CalifaHitData");
branchCalifaHitData->SetAddress(&CalifaHitData);
//
//---------------------------------------------------------------------------
Long64_t nevents = chain->GetEntries(); //number of events in file with name "fname"
//here I check the charge:-----------------
for(Long64_t i=0;i< nevents;i++){
    Long64_t evtnr = i;
	if (i%100000==0)
        cout<<"Processing event for charge analysis "<<i<<endl;
	SofMwpc1HitData->Clear();
    SofMwpc2HitData->Clear();
    SofMwpc3HitData->Clear();
    SofMwpc0HitData->Clear();
    SofSciTcalData->Clear();
    SofToFWTcalData->Clear();
    SofToFWMappedData->Clear();
    SofTwimHitData->Clear();

    chain->GetEvent(i); //event_i_call
	entries_mw0 = SofMwpc0HitData->GetEntries();
    entries_mw1 = SofMwpc1HitData->GetEntries();
    entries_mw2 = SofMwpc2HitData->GetEntries();
    entries_mw3 = SofMwpc3HitData->GetEntries();
    entries_start = SofSciTcalData->GetEntriesFast();
    entries_tofw = SofToFWTcalData->GetEntriesFast();
    entries_califa = CalifaHitData->GetEntries();
    softwimhitdata = new R3BSofTwimHitData*[1];
    softwimhitdata[0] = (R3BSofTwimHitData*)SofTwimHitData->At(0);

	if (entries_califa >= 2 && entries_start ==2 && entries_tofw == 2 && entries_mw0 == 1 && entries_mw1 == 1 && entries_mw2 == 1 && entries_mw3 == 1 && softwimhitdata[0]){
        charge_val = softwimhitdata[0]->GetZcharge();
        //softwimhitdata = new R3BSofTwimHitData*[1];
        sofmwpc0hitdata = new R3BSofMwpcHitData*[1];
        sofmwpc1hitdata = new R3BSofMwpcHitData*[1];
        sofmwpc2hitdata = new R3BSofMwpcHitData*[1];
        sofmwpc3hitdata = new R3BSofMwpcHitData*[1];
        softofwtcaldata = new R3BSofToFWTcalData*[1];
        softofwmappeddata = new R3BSofToFWMappedData*[2];
        califahitdata = new R3BCalifaHitData*[entries_califa];
        sofmwpc0hitdata[0] = (R3BSofMwpcHitData*)SofMwpc0HitData->At(0);
        sofmwpc1hitdata[0] = (R3BSofMwpcHitData*)SofMwpc1HitData->At(0);
        sofmwpc2hitdata[0] = (R3BSofMwpcHitData*)SofMwpc2HitData->At(0);
        sofmwpc3hitdata[0] = (R3BSofMwpcHitData*)SofMwpc3HitData->At(0);
        softofwtcaldata[0] = (R3BSofToFWTcalData*)SofToFWTcalData->At(0);
        //softwimhitdata[0] = (R3BSofTwimHitData*)SofTwimHitData->At(0);
        califahitdata[0] = (R3BCalifaHitData*)CalifaHitData->At(0);
        Double_t xMW0 = sofmwpc0hitdata[0]->GetX();
        Double_t xMW1 = sofmwpc1hitdata[0]->GetX();
        Double_t xMW2 = sofmwpc2hitdata[0]->GetX();
        Double_t xMW3 = sofmwpc3hitdata[0]->GetX();
        Double_t yMW1 = sofmwpc1hitdata[0]->GetY();
        Double_t yMW2 = sofmwpc2hitdata[0]->GetY();
        Double_t yMW3 = sofmwpc3hitdata[0]->GetY();
		
		if (xMW1 != -1000 && xMW2 != -1000 && xMW3 != -1000  && xMW0 != -1000 && yMW3 != -1000){
		//fill here a histogram with charge only....
		h1_z_all_iso->Fill(charge_val);
		}
		 //delete [] softwimhitdata;
         delete [] sofmwpc0hitdata;  
         delete [] sofmwpc1hitdata;  
         delete [] sofmwpc2hitdata;  
         delete [] sofmwpc3hitdata;  
         delete [] softofwtcaldata;  
         delete [] softofwmappeddata;
		 delete [] califahitdata;
		}
		delete [] softwimhitdata;
		}
Int_t binmax_Z6 = h1_z_all_iso->GetMaximumBin();
Double_t middle_Z6 = h1_z_all_iso->GetXaxis()->GetBinCenter(binmax_Z6);
Double_t middle_Z5 = middle_Z6 -1;
cout << "max bin x position" << middle_Z6 << endl;

//TCanvas* canvas = new TCanvas("canvas", "canvas");
//h1_z_all_iso->Draw();
//fit for Z = 6
TF1 *gauss_6 = new TF1("gauss_6","gaus",middle_Z6-0.5,middle_Z6+0.5);
h1_z_all_iso->Fit("gauss_6","R+");

////fit for Z = 5
TF1 *gauss_5 = new TF1("gauss_5","gaus",middle_Z5-0.5,middle_Z5+0.5);
h1_z_all_iso->Fit("gauss_5","R+");

TF1 *fit_g6 = h1_z_all_iso->GetFunction("gauss_6");
TF1 *fit_g5 = h1_z_all_iso->GetFunction("gauss_5");
Double_t mean_g6 = fit_g6->GetParameter(1);
Double_t mean_g5 = fit_g5->GetParameter(1);
Double_t offs_g6 = fit_g6->GetParameter(0);
Double_t offs_g5 = fit_g5->GetParameter(0);
Double_t sigma_g6 = fit_g6->GetParameter(2);
Double_t sigma_g5 = fit_g5->GetParameter(2);

TF1 *f1 = new TF1("f1","abs([0]*exp(-0.5*((x-[1])/[2])^2)-([3]*exp(-0.5*((x-[4])/[5])^2)))",mean_g5,mean_g6 );

f1->SetParameters(offs_g5,mean_g5,sigma_g5,offs_g6,mean_g6,sigma_g6);
f1->Print("V");
cout << "new min" << f1->GetMinimumX() << endl;
const Double_t charge_cut = f1->GetMinimumX();
//delete canvas;
delete gauss_6;
delete gauss_5;
delete f1;
cout << "this is the charge_cut\t" << charge_cut << endl;
cout << "this is mean value of Z = 6\t " << mean_g6 << endl;
cout << "this is mean value of Z = 5\t " << mean_g5 << endl;

//-----------------------------------------

for(Long64_t i=0;i< nevents;i++){
    Long64_t evtnr = i;
    //cout << "This is eventnr.:\t" << evtnr <<endl;
    if (i%100000==0)
        cout<<"Processing event "<<i<<endl;
    SofMwpc1HitData->Clear();
    SofMwpc2HitData->Clear();
    SofMwpc3HitData->Clear();
    SofMwpc0HitData->Clear();
    SofSciTcalData->Clear();
    SofToFWTcalData->Clear();
    SofToFWMappedData->Clear();
    SofTwimHitData->Clear();

    chain->GetEvent(i); //event_i_call
    //checking the number of entries/multiplicity for various detectors...
    entries_mw0 = SofMwpc0HitData->GetEntries();
    entries_mw1 = SofMwpc1HitData->GetEntries();
    entries_mw2 = SofMwpc2HitData->GetEntries();
    entries_mw3 = SofMwpc3HitData->GetEntries();
    entries_start = SofSciTcalData->GetEntriesFast();
    entries_tofw = SofToFWTcalData->GetEntriesFast();
    entries_califa = CalifaHitData->GetEntries();
    softwimhitdata = new R3BSofTwimHitData*[1];
    softwimhitdata[0] = (R3BSofTwimHitData*)SofTwimHitData->At(0);

    //
    if (entries_califa >= 2 && entries_start ==2 && entries_tofw == 2 && entries_mw0 == 1 && entries_mw1 == 1 && entries_mw2 == 1 && entries_mw3 == 1 && softwimhitdata[0]){
        charge_val = softwimhitdata[0]->GetZcharge();
        //softwimhitdata = new R3BSofTwimHitData*[1];
        sofmwpc0hitdata = new R3BSofMwpcHitData*[1];
        sofmwpc1hitdata = new R3BSofMwpcHitData*[1];
        sofmwpc2hitdata = new R3BSofMwpcHitData*[1];
        sofmwpc3hitdata = new R3BSofMwpcHitData*[1];
        softofwtcaldata = new R3BSofToFWTcalData*[1];
        softofwmappeddata = new R3BSofToFWMappedData*[2];
        califahitdata = new R3BCalifaHitData*[entries_califa];
        sofmwpc0hitdata[0] = (R3BSofMwpcHitData*)SofMwpc0HitData->At(0);
        sofmwpc1hitdata[0] = (R3BSofMwpcHitData*)SofMwpc1HitData->At(0);
        sofmwpc2hitdata[0] = (R3BSofMwpcHitData*)SofMwpc2HitData->At(0);
        sofmwpc3hitdata[0] = (R3BSofMwpcHitData*)SofMwpc3HitData->At(0);
        softofwtcaldata[0] = (R3BSofToFWTcalData*)SofToFWTcalData->At(0);
        //softwimhitdata[0] = (R3BSofTwimHitData*)SofTwimHitData->At(0);
        califahitdata[0] = (R3BCalifaHitData*)CalifaHitData->At(0);
        Double_t xMW0 = sofmwpc0hitdata[0]->GetX();
        Double_t xMW1 = sofmwpc1hitdata[0]->GetX();
        Double_t xMW2 = sofmwpc2hitdata[0]->GetX();
        Double_t xMW3 = sofmwpc3hitdata[0]->GetX();
        Double_t yMW1 = sofmwpc1hitdata[0]->GetY();
        Double_t yMW2 = sofmwpc2hitdata[0]->GetY();
        Double_t yMW3 = sofmwpc3hitdata[0]->GetY();
        //here without cuts on y direction:
if (xMW1 != -1000 && xMW2 != -1000 && xMW3 != -1000  && xMW0 != -1000){
        //original formula-----------------------------------------------------
        //Double_t psi_in = atan(((xMW2+xMW2_shift)-(xMW1+xMW1_shift))/(zM2-zT));
        //I think I had a typo in my original formula, it should be zM2-zM1?!:
        Double_t psi_in = atan(((xMW2+xMW2_shift)-(xMW1+xMW1_shift))/(zM2-zM1));
        //------------------------------------------------------------------------------
        Double_t T_to_M2 = (zM2-zT)/(cos(psi_in));
        //new: also here I try to use the shift for better calculation:
        Double_t b1_cutoff = (xMW1+xMW1_shift)-tan(psi_in)*zM1;
        Double_t z_labMW3 = middle_zM3+cos(2*PI/5.)*(xMW3-xMW3_shift);
        Double_t x_labMW3 = middle_xM3+sin(2*PI/5.)*(xMW3-xMW3_shift);
        Double_t zB = (bGLAD_cutoff-b1_cutoff)/(tan(psi_in)-tan(PI/2.-alpha_G));
        Double_t xB = tan(psi_in)*zB+b1_cutoff;

        for (Int_t rec = 0; rec < 50; rec++){
            zC = (zGm_cutoff-b1_cutoff)/(tan(psi_in)-tan(PI/2.-alpha_G)) + 2*rec;
            xC = tan(psi_in)*zC+b1_cutoff;
            slope = (x_labMW3-xC)/(z_labMW3-zC);
            offset_slope = x_labMW3-slope*z_labMW3;
            psi_out_rec[rec] = -atan(slope);
            z_D = (D_cutoff-offset_slope)/(slope-tan(PI/2.-alpha_G));
            x_D = tan(PI/2.-alpha_G)*z_D +D_cutoff;
            l_diff[rec] = (sqrt((xC-xB)*(xC-xB)+(zC-zB)*(zC-zB))-sqrt((x_D-xC)*(x_D-xC)+(z_D-zC)*(z_D-zC)));
            z_pos_shift[rec] = zC;
            x_pos_shift[rec] = xC;
            //cout <<"this is the pathdiff" << l_diff[rec] << endl;
        }
        Int_t n = 50;
        TGraph *gr1 = new TGraph(n,psi_out_rec,l_diff);
        gr1->Fit("pol1","Q");
        TF1 *f3 =gr1->GetFunction("pol1");
        Double_t f3_slope = f3->GetParameter(1);
        Double_t f3_offset = f3->GetParameter(0);
        Double_t psi_out = -f3_offset/f3_slope;
		TGraph *gr2 = new TGraph(n,z_pos_shift,l_diff);
        gr2->Fit("pol1", "Q");
        TF1 *f4 = gr2->GetFunction("pol1");
        Double_t f4_slope = f4->GetParameter(1);
        Double_t f4_offset = f4->GetParameter(0);
        Double_t new_z_C = -f4_offset/f4_slope;

        TGraph *gr3 = new TGraph(n,x_pos_shift,l_diff);
        gr3->Fit("pol1","Q");
        TF1 *f5 = gr3->GetFunction("pol1");
        Double_t f5_slope = f5->GetParameter(1);
        Double_t f5_offset = f5->GetParameter(0);
        Double_t new_x_C = -f5_offset/f5_slope;
//      cout << "this is x" << new_x_C << endl; 
        Double_t z_diff = new_z_C-z_pos_shift[0];
        //cout <<"this is z-diff, shoul be negative" << z_diff << endl;

        //Here I write the exact value for zD and xD:----------------
        Double_t offset_new = new_x_C+tan(psi_out)*new_z_C;
        z_D = (D_cutoff-offset_new)/(-tan(psi_out)-tan(PI/2.-alpha_G));
        x_D = -tan(psi_out)*z_D + offset_new;

        //----calculation of TOFW position---
        Double_t zToFW = (cutoff_ToFW-offset_new)/(-tan(psi_out)-tan(2*PI/5.));
        Double_t xToFW = -tan(psi_out)*zToFW+offset_new;
        Double_t path_D_to_TOFW = (zToFW-z_D)/(cos(psi_out));

        //calculation of distance between B and D:
        Double_t BD = sqrt((x_D-xB)*(x_D-xB)+(z_D-zB)*(z_D-zB));

		Double_t m = (cos(psi_out)-cos(psi_in))/(sin(psi_out)+sin(psi_in));
        Double_t rho =  (Leff/(2*sin((psi_in+psi_out)/2)))*sqrt(pow(((m+tan(alpha_G))/(1-m*tan(alpha_G))),2)+1);
		Double_t w = 2*abs(asin(BD/(2*rho)));


		tot_length = start_to_target+ abs((zB-zT)/(cos(psi_in)))+rho*w+abs((zToFW-z_D)/cos(psi_out));


		if (SofToFWTcalData && SofToFWTcalData->GetEntriesFast())
            {
            //this for-loop initializes the matrix elementys mult[i][j] and set them to 0 beforehand
            for (UShort_t i = 0; i < NbDets; i++)
                {
                for (UShort_t j = 0; j < NbChs; j++)
                    {
                    mult[i * NbChs + j] = 0;
                    }
                }
            //this for-loop fills the two matrix elements which stand for the Preamps with signal
            for (Int_t ihit = 0; ihit < entries_tofw; ihit++)
                {
                R3BSofToFWMappedData* hitmapped = (R3BSofToFWMappedData*)SofToFWMappedData->At(ihit);
                if (!hitmapped)
                    continue;
                iDet = hitmapped->GetDetector() - 1;
                iCh = hitmapped->GetPmt() - 1;
                mult[iDet * NbChs + iCh]++;
                }
            //now we have to get the raw times for the signals at the TOFW using the cal-level
            for (Int_t ihit = 0; ihit < entries_tofw; ihit++)
                {
                R3BSofToFWTcalData* hittcal = (R3BSofToFWTcalData*)SofToFWTcalData->At(ihit);
                if (!hittcal)
                    continue;
                if (hittcal->GetPmt() == 3)
                    continue;
                iDet = hittcal->GetDetector() - 1;
                iCh = hittcal->GetPmt() - 1;
                iRawTimeNs[iDet * 2 + iCh] = hittcal->GetRawTimeNs();
                }
            //now we have to get the raw times for the signals at the START detector using the cal-level
            for (Int_t s_hit = 0 ; s_hit < entries_start; s_hit++)
                {
                R3BSofSciTcalData* start_hits = (R3BSofSciTcalData*)SofSciTcalData->At(s_hit);
                if(!start_hits)
                    continue;
                iRawTimeStartNs[s_hit] = start_hits->GetRawTimeNs();

                }
                sofmwpc3hitdata[0]=(R3BSofMwpcHitData*)SofMwpc3HitData->At(0);
            if ((sofmwpc3hitdata[0]->GetY() > -500.) && (sofmwpc3hitdata[0]->GetX() > -500.))  //cut on the mwpc3
            {
            Double_t raw_mwpc3 = sofmwpc3hitdata[0]->GetY();
            Double_t tofpos_ns = -1000;
            Double_t raw_tofpos_ns = -1000;
            Double_t raw_t_start = -1000000.;
            Double_t raw_time_of_flight = 0.;
			
			for (UShort_t i = 0; i < NbDets; i++)
                {
                if ((mult[i * NbChs] == 1) && (mult[i * NbChs + 1] == 1))
                    {
                    tofpos_ns = 0.5*(iRawTimeNs[i * NbChs + 1] - iRawTimeNs[i * NbChs]-offsets[i+1]);
                    raw_tofpos_ns = 0.5*(iRawTimeNs[i * NbChs + 1] - iRawTimeNs[i * NbChs]);
                    raw_t_start = 0.5*(iRawTimeStartNs[0]-iRawTimeStartNs[1]);
                    if (raw_t_start !=-1000000. && raw_tofpos_ns != -1000 && tofpos_ns != -1000)
                        {
						raw_time_of_flight=(0.5*(iRawTimeNs[i*NbChs+1]+iRawTimeNs[i*NbChs]))-0.5*(iRawTimeStartNs[0]+iRawTimeStartNs[1]);
						Double_t time_target_tof = raw_time_of_flight+tof_offs[i] -time_start_target;
						Double_t path_from_target = abs((zB-zT)/(cos(psi_in)))+rho*w+abs((zToFW-z_D)/cos(psi_out));
						Double_t beta = ((path_from_target/time_target_tof)*pow(10,6))/light_c;
						Double_t mag_field = current*0.0006527728074785267;
						Double_t a_q = (pow(10,-3)*((((mag_field*rho)/(beta))*1.602176634*pow(10,-19))/(1.66053906660*pow(10,-27)*light_c)))/gamma_given;
				//here starts separation of isotopes.....
						if (charge_val< charge_cut && charge_val> charge_cut-1.) //Z = 5
						{
						h1_a_q_z5->Fill(a_q);
						}
						if (charge_val< charge_cut+0.8 && charge_val> charge_cut) //Z = 6
						{
						h1_a_q_z6->Fill(a_q);
						}
						h2_z_vs_a_q_nocut->Fill(a_q,charge_val);
}
}
}
}
}
delete gr1;
delete gr2;
delete gr3;
}
 //delete [] softwimhitdata;
 delete [] sofmwpc0hitdata;  
 delete [] sofmwpc1hitdata;  
 delete [] sofmwpc2hitdata;  
 delete [] sofmwpc3hitdata;  
 delete [] softofwtcaldata;  
 delete [] softofwmappeddata;
 delete [] califahitdata;
}
delete [] softwimhitdata;
}

TF1* gauss_10B = new TF1("gauss_10B","gaus",1.85,2.03);
TF1* gauss_11B = new TF1("gauss_11B","gaus",2.1,2.3);
TF1* gauss_11C = new TF1("gauss_11C","gaus",1.7,1.88);
TF1* gauss_12C = new TF1("gauss_12C","gaus",1.9,2.1);

//
h1_a_q_z6->Fit("gauss_11C","R+");
h1_a_q_z6->Fit("gauss_12C","R+");

TF1* fit_11C = h1_a_q_z6->GetFunction("gauss_11C");
TF1* fit_12C = h1_a_q_z6->GetFunction("gauss_12C");
Double_t offs_11C = fit_11C->GetParameter(0);
Double_t offs_12C = fit_12C->GetParameter(0);
Double_t mean_11C = fit_11C->GetParameter(1);
Double_t mean_12C = fit_12C->GetParameter(1);
Double_t sigma_11C = fit_11C->GetParameter(2);
Double_t sigma_12C = fit_12C->GetParameter(2);

TF1 *f2 = new TF1("f2","abs([0]*exp(-0.5*((x-[1])/[2])^2)-([3]*exp(-0.5*((x-[4])/[5])^2)))",mean_11C,mean_12C);
f2->SetParameters(offs_11C,mean_11C,sigma_11C,offs_12C,mean_12C,sigma_12C);
const Double_t aq_cut_z_6 = f2->GetMinimumX();
const Double_t width_11c = abs(mean_11C-aq_cut_z_6);
const Double_t width_12c = abs(mean_12C-aq_cut_z_6);
cout << "cut on z6" << aq_cut_z_6 << endl;
//

h1_a_q_z5->Fit("gauss_10B","R+");
h1_a_q_z5->Fit("gauss_11B","R+");
TF1* fit_10B = h1_a_q_z5->GetFunction("gauss_10B");
TF1* fit_11B = h1_a_q_z5->GetFunction("gauss_11B");
Double_t offs_10B = fit_10B->GetParameter(0);
Double_t offs_11B = fit_11B->GetParameter(0);
Double_t mean_10B = fit_10B->GetParameter(1);
Double_t mean_11B = fit_11B->GetParameter(1);
Double_t sigma_10B = fit_10B->GetParameter(2);
Double_t sigma_11B = fit_11B->GetParameter(2);
TF1 *f3 = new TF1("f3","abs([0]*exp(-0.5*((x-[1])/[2])^2)-([3]*exp(-0.5*((x-[4])/[5])^2)))",mean_10B,mean_11B);
f3->SetParameters(offs_10B,mean_10B,sigma_10B,offs_11B,mean_11B,sigma_11B);
f3->Print("V");
const Double_t aq_cut_z_5 = f3->GetMinimumX();
const Double_t width_10b = abs(mean_10B-aq_cut_z_5);
const Double_t width_11b = abs(mean_11B-aq_cut_z_5);
cout << "cut on z5" << aq_cut_z_5 << endl;


delete gauss_10B;
delete gauss_11B;
delete gauss_11C;
delete gauss_12C;
delete f3;
delete f2;

cout << "----------------------------------------------------------" << endl;
cout << "----------------------------------------------------------" << endl;
cout << "cut on z5\t" << aq_cut_z_5 << endl;
cout << "cut on z6\t" << aq_cut_z_6 << endl;
cout << "----------------------------------------------------------" << endl;
cout << "witth 10b \t" << width_10b << endl;
cout << "width 11b \t" << width_11b << endl;
cout << "width 11c \t" << width_11c << endl;
cout << "width 12c \t" << width_12c << endl;
cout << "----------------------------------------------------------" << endl;
cout << "mean_10b \t" << mean_10B << endl;
cout << "mean_10b -cut on z5\t " << mean_10B-aq_cut_z_5 << endl;
cout << "mean_11b \t" << mean_11B << endl; 
cout << "mean_11b -cut on z5 \t" << mean_11B-aq_cut_z_5 << endl;
cout << "mean_11c \t" << mean_11C << endl;
cout << "mean_12c \t" << mean_12C << endl;


cout << "----------------------------------------------------------" << endl;
cout << "----------------------------------------------------------" << endl;

//third step to get variables of psi_in:------------------------------------------------------------------------

Int_t entries_10b = 0;
Int_t entries_11b = 0;
Int_t entries_11c = 0;
Int_t entries_12c = 0;

for(Long64_t i=0;i< nevents;i++){
    Long64_t evtnr = i;
    //cout << "This is eventnr.:\t" << evtnr <<endl;
    if (i%100000==0)
        cout<<"Processing event "<<i<<endl;
    SofMwpc1HitData->Clear();
    SofMwpc2HitData->Clear();
    SofMwpc3HitData->Clear();
    SofMwpc0HitData->Clear();
    SofSciTcalData->Clear();
    SofToFWTcalData->Clear();
    SofToFWMappedData->Clear();
    SofTwimHitData->Clear();

    chain->GetEvent(i); //event_i_call
    //checking the number of entries/multiplicity for various detectors...
    entries_mw0 = SofMwpc0HitData->GetEntries();
    entries_mw1 = SofMwpc1HitData->GetEntries();
    entries_mw2 = SofMwpc2HitData->GetEntries();
    entries_mw3 = SofMwpc3HitData->GetEntries();
    entries_start = SofSciTcalData->GetEntriesFast();
    entries_tofw = SofToFWTcalData->GetEntriesFast();
    entries_califa = CalifaHitData->GetEntries();
    softwimhitdata = new R3BSofTwimHitData*[1];
    softwimhitdata[0] = (R3BSofTwimHitData*)SofTwimHitData->At(0);

    //
    if (entries_califa >= 2 && entries_start ==2 && entries_tofw == 2 && entries_mw0 == 1 && entries_mw1 == 1 && entries_mw2 == 1 && entries_mw3 == 1 && softwimhitdata[0]){
        charge_val = softwimhitdata[0]->GetZcharge();
        //softwimhitdata = new R3BSofTwimHitData*[1];
        sofmwpc0hitdata = new R3BSofMwpcHitData*[1];
        sofmwpc1hitdata = new R3BSofMwpcHitData*[1];
        sofmwpc2hitdata = new R3BSofMwpcHitData*[1];
        sofmwpc3hitdata = new R3BSofMwpcHitData*[1];
        softofwtcaldata = new R3BSofToFWTcalData*[1];
        softofwmappeddata = new R3BSofToFWMappedData*[2];
        califahitdata = new R3BCalifaHitData*[entries_califa];
        sofmwpc0hitdata[0] = (R3BSofMwpcHitData*)SofMwpc0HitData->At(0);
        sofmwpc1hitdata[0] = (R3BSofMwpcHitData*)SofMwpc1HitData->At(0);
        sofmwpc2hitdata[0] = (R3BSofMwpcHitData*)SofMwpc2HitData->At(0);
        sofmwpc3hitdata[0] = (R3BSofMwpcHitData*)SofMwpc3HitData->At(0);
        softofwtcaldata[0] = (R3BSofToFWTcalData*)SofToFWTcalData->At(0);
        //softwimhitdata[0] = (R3BSofTwimHitData*)SofTwimHitData->At(0);
        califahitdata[0] = (R3BCalifaHitData*)CalifaHitData->At(0);
        Double_t xMW0 = sofmwpc0hitdata[0]->GetX();
        Double_t xMW1 = sofmwpc1hitdata[0]->GetX();
        Double_t xMW2 = sofmwpc2hitdata[0]->GetX();
        Double_t xMW3 = sofmwpc3hitdata[0]->GetX();
        Double_t yMW1 = sofmwpc1hitdata[0]->GetY();
        Double_t yMW2 = sofmwpc2hitdata[0]->GetY();
        Double_t yMW3 = sofmwpc3hitdata[0]->GetY();
        //here without cuts on y direction:
if (xMW1 != -1000 && xMW2 != -1000 && xMW3 != -1000  && xMW0 != -1000){
        //original formula-----------------------------------------------------
        //Double_t psi_in = atan(((xMW2+xMW2_shift)-(xMW1+xMW1_shift))/(zM2-zT));
        //I think I had a typo in my original formula, it should be zM2-zM1?!:
        Double_t psi_in = atan(((xMW2+xMW2_shift)-(xMW1+xMW1_shift))/(zM2-zM1));
        //------------------------------------------------------------------------------
        Double_t T_to_M2 = (zM2-zT)/(cos(psi_in));
        //new: also here I try to use the shift for better calculation:
        Double_t b1_cutoff = (xMW1+xMW1_shift)-tan(psi_in)*zM1;
        Double_t z_labMW3 = middle_zM3+cos(2*PI/5.)*(xMW3-xMW3_shift);
        Double_t x_labMW3 = middle_xM3+sin(2*PI/5.)*(xMW3-xMW3_shift);
        Double_t zB = (bGLAD_cutoff-b1_cutoff)/(tan(psi_in)-tan(PI/2.-alpha_G));
        Double_t xB = tan(psi_in)*zB+b1_cutoff;

        for (Int_t rec = 0; rec < 50; rec++){
            zC = (zGm_cutoff-b1_cutoff)/(tan(psi_in)-tan(PI/2.-alpha_G)) + 2*rec;
            xC = tan(psi_in)*zC+b1_cutoff;
            slope = (x_labMW3-xC)/(z_labMW3-zC);
            offset_slope = x_labMW3-slope*z_labMW3;
            psi_out_rec[rec] = -atan(slope);
            z_D = (D_cutoff-offset_slope)/(slope-tan(PI/2.-alpha_G));
            x_D = tan(PI/2.-alpha_G)*z_D +D_cutoff;
            l_diff[rec] = (sqrt((xC-xB)*(xC-xB)+(zC-zB)*(zC-zB))-sqrt((x_D-xC)*(x_D-xC)+(z_D-zC)*(z_D-zC)));
            z_pos_shift[rec] = zC;
            x_pos_shift[rec] = xC;
            //cout <<"this is the pathdiff" << l_diff[rec] << endl;
        }
        Int_t n = 50;
        TGraph *gr1 = new TGraph(n,psi_out_rec,l_diff);
        gr1->Fit("pol1","Q");
        TF1 *f3 =gr1->GetFunction("pol1");
        Double_t f3_slope = f3->GetParameter(1);
        Double_t f3_offset = f3->GetParameter(0);
        Double_t psi_out = -f3_offset/f3_slope;
		TGraph *gr2 = new TGraph(n,z_pos_shift,l_diff);
        gr2->Fit("pol1", "Q");
        TF1 *f4 = gr2->GetFunction("pol1");
        Double_t f4_slope = f4->GetParameter(1);
        Double_t f4_offset = f4->GetParameter(0);
        Double_t new_z_C = -f4_offset/f4_slope;

        TGraph *gr3 = new TGraph(n,x_pos_shift,l_diff);
        gr3->Fit("pol1","Q");
        TF1 *f5 = gr3->GetFunction("pol1");
        Double_t f5_slope = f5->GetParameter(1);
        Double_t f5_offset = f5->GetParameter(0);
        Double_t new_x_C = -f5_offset/f5_slope;
//      cout << "this is x" << new_x_C << endl; 
        Double_t z_diff = new_z_C-z_pos_shift[0];
        //cout <<"this is z-diff, shoul be negative" << z_diff << endl;

        //Here I write the exact value for zD and xD:----------------
        Double_t offset_new = new_x_C+tan(psi_out)*new_z_C;
        z_D = (D_cutoff-offset_new)/(-tan(psi_out)-tan(PI/2.-alpha_G));
        x_D = -tan(psi_out)*z_D + offset_new;

        //----calculation of TOFW position---
        Double_t zToFW = (cutoff_ToFW-offset_new)/(-tan(psi_out)-tan(2*PI/5.));
        Double_t xToFW = -tan(psi_out)*zToFW+offset_new;
        Double_t path_D_to_TOFW = (zToFW-z_D)/(cos(psi_out));

        //calculation of distance between B and D:
        Double_t BD = sqrt((x_D-xB)*(x_D-xB)+(z_D-zB)*(z_D-zB));

		Double_t m = (cos(psi_out)-cos(psi_in))/(sin(psi_out)+sin(psi_in));
        Double_t rho =  (Leff/(2*sin((psi_in+psi_out)/2)))*sqrt(pow(((m+tan(alpha_G))/(1-m*tan(alpha_G))),2)+1);
		Double_t w = 2*abs(asin(BD/(2*rho)));


		tot_length = start_to_target+ abs((zB-zT)/(cos(psi_in)))+rho*w+abs((zToFW-z_D)/cos(psi_out));


		if (SofToFWTcalData && SofToFWTcalData->GetEntriesFast())
            {
            //this for-loop initializes the matrix elementys mult[i][j] and set them to 0 beforehand
            for (UShort_t i = 0; i < NbDets; i++)
                {
                for (UShort_t j = 0; j < NbChs; j++)
                    {
                    mult[i * NbChs + j] = 0;
                    }
                }
            //this for-loop fills the two matrix elements which stand for the Preamps with signal
            for (Int_t ihit = 0; ihit < entries_tofw; ihit++)
                {
                R3BSofToFWMappedData* hitmapped = (R3BSofToFWMappedData*)SofToFWMappedData->At(ihit);
                if (!hitmapped)
                    continue;
                iDet = hitmapped->GetDetector() - 1;
                iCh = hitmapped->GetPmt() - 1;
                mult[iDet * NbChs + iCh]++;
                }
            //now we have to get the raw times for the signals at the TOFW using the cal-level
            for (Int_t ihit = 0; ihit < entries_tofw; ihit++)
                {
                R3BSofToFWTcalData* hittcal = (R3BSofToFWTcalData*)SofToFWTcalData->At(ihit);
                if (!hittcal)
                    continue;
                if (hittcal->GetPmt() == 3)
                    continue;
                iDet = hittcal->GetDetector() - 1;
                iCh = hittcal->GetPmt() - 1;
                iRawTimeNs[iDet * 2 + iCh] = hittcal->GetRawTimeNs();
                }
            //now we have to get the raw times for the signals at the START detector using the cal-level
            for (Int_t s_hit = 0 ; s_hit < entries_start; s_hit++)
                {
                R3BSofSciTcalData* start_hits = (R3BSofSciTcalData*)SofSciTcalData->At(s_hit);
                if(!start_hits)
                    continue;
                iRawTimeStartNs[s_hit] = start_hits->GetRawTimeNs();

                }
                sofmwpc3hitdata[0]=(R3BSofMwpcHitData*)SofMwpc3HitData->At(0);
            if ((sofmwpc3hitdata[0]->GetY() > -500.) && (sofmwpc3hitdata[0]->GetX() > -500.))  //cut on the mwpc3
            {
            Double_t raw_mwpc3 = sofmwpc3hitdata[0]->GetY();
            Double_t tofpos_ns = -1000;
            Double_t raw_tofpos_ns = -1000;
            Double_t raw_t_start = -1000000.;
            Double_t raw_time_of_flight = 0.;
			
			for (UShort_t i = 0; i < NbDets; i++)
                {
                if ((mult[i * NbChs] == 1) && (mult[i * NbChs + 1] == 1))
                    {
                    tofpos_ns = 0.5*(iRawTimeNs[i * NbChs + 1] - iRawTimeNs[i * NbChs]-offsets[i+1]);
                    raw_tofpos_ns = 0.5*(iRawTimeNs[i * NbChs + 1] - iRawTimeNs[i * NbChs]);
                    raw_t_start = 0.5*(iRawTimeStartNs[0]-iRawTimeStartNs[1]);
                    if (raw_t_start !=-1000000. && raw_tofpos_ns != -1000 && tofpos_ns != -1000)
                        {
						raw_time_of_flight=(0.5*(iRawTimeNs[i*NbChs+1]+iRawTimeNs[i*NbChs]))-0.5*(iRawTimeStartNs[0]+iRawTimeStartNs[1]);
						Double_t time_target_tof = raw_time_of_flight+tof_offs[i] -time_start_target;
						Double_t path_from_target = abs((zB-zT)/(cos(psi_in)))+rho*w+abs((zToFW-z_D)/cos(psi_out));
						Double_t beta = ((path_from_target/time_target_tof)*pow(10,6))/light_c;
						Double_t mag_field = current*0.0006527728074785267;
						Double_t a_q = (pow(10,-3)*((((mag_field*rho)/(beta))*1.602176634*pow(10,-19))/(1.66053906660*pow(10,-27)*light_c)))/gamma_given;
				//here starts separation of isotopes.....
						if (charge_val< charge_cut && charge_val> charge_cut-1.) //Z = 5 before it was charge_val< charge_cut && charge_val> charge_cut-0.8, too straight....
						{
						if (a_q > aq_cut_z_5 && a_q < (aq_cut_z_5+2*width_11b)){  //11B
							h2_theta_out_vs_mw3_11b->Fill(xMW3,psi_out);
							entries_11b++;
						}						
						if (a_q < aq_cut_z_5 && a_q > (aq_cut_z_5-0.15)){  //10B
							h2_theta_out_vs_mw3_10b->Fill(xMW3,psi_out);
							entries_10b++;
                        }
						}
						if (charge_val< charge_cut+0.8 && charge_val> charge_cut) //Z = 6
						{
						if (a_q > aq_cut_z_6 && a_q < (aq_cut_z_6+2*width_12c)){  //12C
							h2_theta_out_vs_mw3_12c->Fill(xMW3,psi_out);
							entries_12c++;
                        }
						if (a_q < aq_cut_z_6-0.06 && a_q > (aq_cut_z_6-1.6*width_11c)){  //11C before it wasa_q < aq_cut_z_6 && a_q > (aq_cut_z_6-2*width_11c), bad for 11C ...
							h2_theta_out_vs_mw3_11c->Fill(xMW3,psi_out);
							entries_11c++;

                        }

						}
}
}
}
}
}
delete gr1;
delete gr2;
delete gr3;
}
 delete [] sofmwpc0hitdata;  
 delete [] sofmwpc1hitdata;  
 delete [] sofmwpc2hitdata;  
 delete [] sofmwpc3hitdata;  
 delete [] softofwtcaldata;  
 delete [] softofwmappeddata;
 delete [] califahitdata;
}
delete [] softwimhitdata;
}

//now get ranges for the fits:
//10B
Double_t st_dev_mw3_10b = h2_theta_out_vs_mw3_10b->GetStdDev(1);
Double_t st_dev_psi_out_10b = h2_theta_out_vs_mw3_10b->GetStdDev(2);
Double_t mean_mw3_10b = h2_theta_out_vs_mw3_10b->GetMean(1);
const Double_t mean_psi_out_10b = h2_theta_out_vs_mw3_10b->GetMean(2);
//
//11B
Double_t st_dev_mw3_11b = h2_theta_out_vs_mw3_11b->GetStdDev(1);
Double_t st_dev_psi_out_11b = h2_theta_out_vs_mw3_11b->GetStdDev(2);
Double_t mean_mw3_11b = h2_theta_out_vs_mw3_11b->GetMean(1);
const Double_t mean_psi_out_11b = h2_theta_out_vs_mw3_11b->GetMean(2);
//

//12C
Double_t st_dev_mw3_12c = h2_theta_out_vs_mw3_12c->GetStdDev(1);
Double_t st_dev_psi_out_12c = h2_theta_out_vs_mw3_12c->GetStdDev(2);
Double_t mean_mw3_12c = h2_theta_out_vs_mw3_12c->GetMean(1);
const Double_t mean_psi_out_12c = h2_theta_out_vs_mw3_12c->GetMean(2);
//

//11C

Double_t st_dev_mw3_11c = h2_theta_out_vs_mw3_11c->GetStdDev(1);
Double_t st_dev_psi_out_11c = h2_theta_out_vs_mw3_11c->GetStdDev(2);
Double_t mean_mw3_11c = h2_theta_out_vs_mw3_11c->GetMean(1);
const Double_t mean_psi_out_11c = h2_theta_out_vs_mw3_11c->GetMean(2);
//

TF1 *angle_fit_10b = new TF1("angle_fit_10b","pol1",mean_mw3_10b-2*st_dev_mw3_10b,mean_mw3_10b+2*st_dev_mw3_10b);
TF1 *angle_fit_11b = new TF1("angle_fit_11b","pol1",mean_mw3_11b-2*st_dev_mw3_11b,mean_mw3_11b+2*st_dev_mw3_11b);
TF1 *angle_fit_11c = new TF1("angle_fit_11c","pol1",mean_mw3_11c-2*st_dev_mw3_11c,mean_mw3_11c+2*st_dev_mw3_11c);
TF1 *angle_fit_12c = new TF1("angle_fit_12c","pol1",mean_mw3_12c-2*st_dev_mw3_12c,mean_mw3_12c+2*st_dev_mw3_12c);


//psi_parameters for 10B:
h2_theta_out_vs_mw3_10b->Fit("angle_fit_10b", "R+");
h2_theta_out_vs_mw3_10b->Draw();
TF1* fit_angle_10b = h2_theta_out_vs_mw3_10b->GetFunction("angle_fit_10b");
const Double_t angle_offs_10b = fit_angle_10b->GetParameter(0);
const Double_t slope_10b = fit_angle_10b->GetParameter(1);
cout << "this is the offs_angle of 10B" << angle_offs_10b << endl;
cout << "this is the mean theta_out of 10B" << mean_psi_out_10b << endl;
cout << "this is the slope of 10B" << slope_10b << endl;
//

//psi_parameters for 11B:
h2_theta_out_vs_mw3_11b->Fit("angle_fit_11b", "R+");
TF1* fit_angle_11b = h2_theta_out_vs_mw3_11b->GetFunction("angle_fit_11b");
const Double_t angle_offs_11b = fit_angle_11b->GetParameter(0); 
const Double_t slope_11b = fit_angle_11b->GetParameter(1);
//

//psi_parameters for 11C:
h2_theta_out_vs_mw3_11c->Fit("angle_fit_11c", "R+");
TF1* fit_angle_11c = h2_theta_out_vs_mw3_11c->GetFunction("angle_fit_11c");
const Double_t angle_offs_11c = fit_angle_11c->GetParameter(0);
const Double_t slope_11c = fit_angle_11c->GetParameter(1);
//

//psi_parameters for 12C:
h2_theta_out_vs_mw3_12c->Fit("angle_fit_12c", "R+");
TF1* fit_angle_12c = h2_theta_out_vs_mw3_12c->GetFunction("angle_fit_12c");
const Double_t angle_offs_12c = fit_angle_12c->GetParameter(0);
const Double_t slope_12c = fit_angle_12c->GetParameter(1);
//

//

//-----------------------------------------------------------------------------------------------------------------------------

//Finally the precise calculation with the parameters from Step1 -> Step3

for(Long64_t i=0;i< nevents;i++){
    Long64_t evtnr = i;
    //cout << "This is eventnr.:\t" << evtnr <<endl;
    if (i%100000==0)
        cout<<"Processing event "<<i<<endl;
    SofMwpc1HitData->Clear();
    SofMwpc2HitData->Clear();
    SofMwpc3HitData->Clear();
    SofMwpc0HitData->Clear();
    SofSciTcalData->Clear();
    SofToFWTcalData->Clear();
    SofToFWMappedData->Clear();
    SofTwimHitData->Clear();

    chain->GetEvent(i); //event_i_call
    //checking the number of entries/multiplicity for various detectors...
    entries_mw0 = SofMwpc0HitData->GetEntries();
    entries_mw1 = SofMwpc1HitData->GetEntries();
    entries_mw2 = SofMwpc2HitData->GetEntries();
    entries_mw3 = SofMwpc3HitData->GetEntries();
    entries_start = SofSciTcalData->GetEntriesFast();
    entries_tofw = SofToFWTcalData->GetEntriesFast();
    entries_califa = CalifaHitData->GetEntries();
    softwimhitdata = new R3BSofTwimHitData*[1];
    softwimhitdata[0] = (R3BSofTwimHitData*)SofTwimHitData->At(0);

    //
    if (entries_start ==2 && entries_tofw == 2 && entries_mw0 == 1 && entries_mw1 == 1 && entries_mw2 == 1 && entries_mw3 == 1 && softwimhitdata[0]){
        charge_val = softwimhitdata[0]->GetZcharge();
        //softwimhitdata = new R3BSofTwimHitData*[1];
        sofmwpc0hitdata = new R3BSofMwpcHitData*[1];
        sofmwpc1hitdata = new R3BSofMwpcHitData*[1];
        sofmwpc2hitdata = new R3BSofMwpcHitData*[1];
        sofmwpc3hitdata = new R3BSofMwpcHitData*[1];
        softofwtcaldata = new R3BSofToFWTcalData*[1];
        softofwmappeddata = new R3BSofToFWMappedData*[2];
        //califahitdata = new R3BCalifaHitData*[entries_califa];
        sofmwpc0hitdata[0] = (R3BSofMwpcHitData*)SofMwpc0HitData->At(0);
        sofmwpc1hitdata[0] = (R3BSofMwpcHitData*)SofMwpc1HitData->At(0);
        sofmwpc2hitdata[0] = (R3BSofMwpcHitData*)SofMwpc2HitData->At(0);
        sofmwpc3hitdata[0] = (R3BSofMwpcHitData*)SofMwpc3HitData->At(0);
        softofwtcaldata[0] = (R3BSofToFWTcalData*)SofToFWTcalData->At(0);
        //softwimhitdata[0] = (R3BSofTwimHitData*)SofTwimHitData->At(0);
        //califahitdata[0] = (R3BCalifaHitData*)CalifaHitData->At(0);
        Double_t xMW0 = sofmwpc0hitdata[0]->GetX();
        Double_t xMW1 = sofmwpc1hitdata[0]->GetX();
        Double_t xMW2 = sofmwpc2hitdata[0]->GetX();
        Double_t xMW3 = sofmwpc3hitdata[0]->GetX();
        Double_t yMW1 = sofmwpc1hitdata[0]->GetY();
        Double_t yMW2 = sofmwpc2hitdata[0]->GetY();
        Double_t yMW3 = sofmwpc3hitdata[0]->GetY();

        //here without cuts on y direction:
if (xMW1 != -1000 && xMW2 != -1000 && xMW3 != -1000  && xMW0 != -1000){
        //original formula-----------------------------------------------------
        //Double_t psi_in = atan(((xMW2+xMW2_shift)-(xMW1+xMW1_shift))/(zM2-zT));
        //I think I had a typo in my original formula, it should be zM2-zM1?!:
        Double_t psi_in = atan(((xMW2+xMW2_shift)-(xMW1+xMW1_shift))/(zM2-zM1));
        //------------------------------------------------------------------------------
        Double_t T_to_M2 = (zM2-zT)/(cos(psi_in));
        //new: also here I try to use the shift for better calculation:
        Double_t b1_cutoff = (xMW1+xMW1_shift)-tan(psi_in)*zM1;
        Double_t z_labMW3 = middle_zM3+cos(2*PI/5.)*(xMW3-xMW3_shift);
        Double_t x_labMW3 = middle_xM3+sin(2*PI/5.)*(xMW3-xMW3_shift);
        Double_t zB = (bGLAD_cutoff-b1_cutoff)/(tan(psi_in)-tan(PI/2.-alpha_G));
        Double_t xB = tan(psi_in)*zB+b1_cutoff;

        for (Int_t rec = 0; rec < 50; rec++){
            zC = (zGm_cutoff-b1_cutoff)/(tan(psi_in)-tan(PI/2.-alpha_G)) + 2*rec;
            xC = tan(psi_in)*zC+b1_cutoff;
            slope = (x_labMW3-xC)/(z_labMW3-zC);
            offset_slope = x_labMW3-slope*z_labMW3;
            psi_out_rec[rec] = -atan(slope);
            z_D = (D_cutoff-offset_slope)/(slope-tan(PI/2.-alpha_G));
            x_D = tan(PI/2.-alpha_G)*z_D +D_cutoff;
            l_diff[rec] = (sqrt((xC-xB)*(xC-xB)+(zC-zB)*(zC-zB))-sqrt((x_D-xC)*(x_D-xC)+(z_D-zC)*(z_D-zC)));
            z_pos_shift[rec] = zC;
            x_pos_shift[rec] = xC;
            //cout <<"this is the pathdiff" << l_diff[rec] << endl;
        }
        Int_t n = 50;
        TGraph *gr1 = new TGraph(n,psi_out_rec,l_diff);
        gr1->Fit("pol1","Q");
        TF1 *f3 =gr1->GetFunction("pol1");
        Double_t f3_slope = f3->GetParameter(1);
        Double_t f3_offset = f3->GetParameter(0);
        Double_t psi_out = -f3_offset/f3_slope;
		TGraph *gr2 = new TGraph(n,z_pos_shift,l_diff);
        gr2->Fit("pol1", "Q");
        TF1 *f4 = gr2->GetFunction("pol1");
        Double_t f4_slope = f4->GetParameter(1);
        Double_t f4_offset = f4->GetParameter(0);
        Double_t new_z_C = -f4_offset/f4_slope;

        TGraph *gr3 = new TGraph(n,x_pos_shift,l_diff);
        gr3->Fit("pol1","Q");
        TF1 *f5 = gr3->GetFunction("pol1");
        Double_t f5_slope = f5->GetParameter(1);
        Double_t f5_offset = f5->GetParameter(0);
        Double_t new_x_C = -f5_offset/f5_slope;
//      cout << "this is x" << new_x_C << endl; 
        Double_t z_diff = new_z_C-z_pos_shift[0];
        //cout <<"this is z-diff, shoul be negative" << z_diff << endl;

        //Here I write the exact value for zD and xD:----------------
        Double_t offset_new = new_x_C+tan(psi_out)*new_z_C;
        z_D = (D_cutoff-offset_new)/(-tan(psi_out)-tan(PI/2.-alpha_G));
        x_D = -tan(psi_out)*z_D + offset_new;

        //----calculation of TOFW position---
        Double_t zToFW = (cutoff_ToFW-offset_new)/(-tan(psi_out)-tan(2*PI/5.));
        Double_t xToFW = -tan(psi_out)*zToFW+offset_new;
        Double_t path_D_to_TOFW = (zToFW-z_D)/(cos(psi_out));

        //calculation of distance between B and D:
        Double_t BD = sqrt((x_D-xB)*(x_D-xB)+(z_D-zB)*(z_D-zB));

		Double_t m = (cos(psi_out)-cos(psi_in))/(sin(psi_out)+sin(psi_in));
        Double_t rho =  (Leff/(2*sin((psi_in+psi_out)/2)))*sqrt(pow(((m+tan(alpha_G))/(1-m*tan(alpha_G))),2)+1);
		Double_t w = 2*abs(asin(BD/(2*rho)));


		tot_length = start_to_target+ abs((zB-zT)/(cos(psi_in)))+rho*w+abs((zToFW-z_D)/cos(psi_out));


		if (SofToFWTcalData && SofToFWTcalData->GetEntriesFast())
            {
            //this for-loop initializes the matrix elementys mult[i][j] and set them to 0 beforehand
            for (UShort_t i = 0; i < NbDets; i++)
                {
                for (UShort_t j = 0; j < NbChs; j++)
                    {
                    mult[i * NbChs + j] = 0;
                    }
                }
            //this for-loop fills the two matrix elements which stand for the Preamps with signal
            for (Int_t ihit = 0; ihit < entries_tofw; ihit++)
                {
                R3BSofToFWMappedData* hitmapped = (R3BSofToFWMappedData*)SofToFWMappedData->At(ihit);
                if (!hitmapped)
                    continue;
                iDet = hitmapped->GetDetector() - 1;
                iCh = hitmapped->GetPmt() - 1;
                mult[iDet * NbChs + iCh]++;
                }
            //now we have to get the raw times for the signals at the TOFW using the cal-level
            for (Int_t ihit = 0; ihit < entries_tofw; ihit++)
                {
                R3BSofToFWTcalData* hittcal = (R3BSofToFWTcalData*)SofToFWTcalData->At(ihit);
                if (!hittcal)
                    continue;
                if (hittcal->GetPmt() == 3)
                    continue;
                iDet = hittcal->GetDetector() - 1;
                iCh = hittcal->GetPmt() - 1;
                iRawTimeNs[iDet * 2 + iCh] = hittcal->GetRawTimeNs();
                }
            //now we have to get the raw times for the signals at the START detector using the cal-level
            for (Int_t s_hit = 0 ; s_hit < entries_start; s_hit++)
                {
                R3BSofSciTcalData* start_hits = (R3BSofSciTcalData*)SofSciTcalData->At(s_hit);
                if(!start_hits)
                    continue;
                iRawTimeStartNs[s_hit] = start_hits->GetRawTimeNs();

                }
                sofmwpc3hitdata[0]=(R3BSofMwpcHitData*)SofMwpc3HitData->At(0);
            if ((sofmwpc3hitdata[0]->GetY() > -500.) && (sofmwpc3hitdata[0]->GetX() > -500.))  //cut on the mwpc3
            {
            Double_t raw_mwpc3 = sofmwpc3hitdata[0]->GetY();
            Double_t tofpos_ns = -1000;
            Double_t raw_tofpos_ns = -1000;
            Double_t raw_t_start = -1000000.;
            Double_t raw_time_of_flight = 0.;
			
			for (UShort_t i = 0; i < NbDets; i++)
                {
                if ((mult[i * NbChs] == 1) && (mult[i * NbChs + 1] == 1))
                    {
                    tofpos_ns = 0.5*(iRawTimeNs[i * NbChs + 1] - iRawTimeNs[i * NbChs]-offsets[i+1]);
                    raw_tofpos_ns = 0.5*(iRawTimeNs[i * NbChs + 1] - iRawTimeNs[i * NbChs]);
                    raw_t_start = 0.5*(iRawTimeStartNs[0]-iRawTimeStartNs[1]);
                    if (raw_t_start !=-1000000. && raw_tofpos_ns != -1000 && tofpos_ns != -1000)
                        {
						raw_time_of_flight=(0.5*(iRawTimeNs[i*NbChs+1]+iRawTimeNs[i*NbChs]))-0.5*(iRawTimeStartNs[0]+iRawTimeStartNs[1]);
						Double_t time_target_tof = raw_time_of_flight+tof_offs[i] -time_start_target;
						Double_t path_from_target = abs((zB-zT)/(cos(psi_in)))+rho*w+abs((zToFW-z_D)/cos(psi_out));
						Double_t beta = ((path_from_target/time_target_tof)*pow(10,6))/light_c;
						Double_t mag_field = current*0.0006527728074785267;
						Double_t a_q = (pow(10,-3)*((((mag_field*rho)/(beta))*1.602176634*pow(10,-19))/(1.66053906660*pow(10,-27)*light_c)))/gamma_given;
				//here starts separation of isotopes.....
						if (charge_val< charge_cut && charge_val> charge_cut-1.) //Z = 5 before it was charge_val< charge_cut && charge_val> charge_cut-0.8, too straight ....
						{

						if (a_q > aq_cut_z_5 && a_q < (aq_cut_z_5+2*width_11b) && entries_califa >= 2){  //11B
						psi_in = mean_psi_out_11b -slope_11b*xMW3 - angle_offs_11b;
						m = (cos(psi_out)-cos(psi_in))/(sin(psi_out)+sin(psi_in));
						rho =  (Leff/(2*sin((psi_in+psi_out)/2)))*sqrt(pow(((m+tan(alpha_G))/(1-m*tan(alpha_G))),2)+1);
                        w = 2*abs(asin(BD/(2*rho)));
						Double_t pathlength_precise = abs((zB-zT)/(cos(psi_in)))+rho*w+abs((zToFW-z_D)/cos(psi_out));
						Double_t beta_precise = ((pathlength_precise/time_target_tof)*pow(10,6))/light_c;
						Double_t gamma_precise = 1/(sqrt(1-beta_precise*beta_precise));
						Double_t a_q_precise = (pow(10,-3)*((((mag_field*rho)/(beta_precise))*1.602176634*pow(10,-19))/(1.66053906660*pow(10,-27)*light_c)))/gamma_given;
						//cout << "beta is:\t " << beta << endl;
						//Double_t mom_11b = a_q_precise*5*(beta/(sqrt(1-(beta*beta))))*938.272;
						//TJ some test, should I do it in that way?
						Double_t mom_11b = 11*931.5*(beta_precise/(sqrt(1-(beta_precise*beta_precise))));
						//cout << "Long mom 11B:\t " << mom_11b*cos(psi_in) << endl;
						h2_z_vs_a_q->Fill(a_q_precise,charge_val);
                        h2_tof_vs_aq_fix_g_11b->Fill(a_q_precise,time_target_tof);
						
						//CALIFA handling
						Double_t psi_CALIFA;
                        Double_t theta_CALIFA;
						Double_t time_CALIFA;
                        Double_t E_Califa;
                        vector<double> vE1;
                        vector<double> vtheta1CALIFA;
                        vector<double> vpsi1CALIFA;
                        vector<double> vE2;
                        vector<double> vtheta2CALIFA;
                        vector<double> vpsi2CALIFA;
						//TJ time measurement CALIFA
						vector<double> vCALIFA_time_1;
						vector<double> vCALIFA_time_2;
						vector<double> v_E_gamma;
						//
						//TJ do the vector handling in a different way...
						vector<double> vpsi;
						vector<double> vtheta;
						vector<double> vE;
						vector<double> v_gamma_sorted;
						vector<double> vtime;
						vector<double> v_origin_E;
						Double_t ind_high_energy;
						Double_t ind_second_high_energy;
						Double_t ind_first_gamma;
						Double_t ind_second_gamma;
						Double_t theta_gamma;
						Double_t psi1;
						Double_t psi2;
						Double_t psi1_2pi;
						Double_t psi2_2pi;
						Double_t psi1_2pi_degr;
						Double_t psi2_2pi_degr;
						Double_t theta1_degr;
						Double_t theta2_degr;
						//Double_t first_gamma;
						//Double_t second_gamma;
						//Int_t no_third_gamma = 0;
						//Double_t energy_gamma_first;
						//Double_t energy_gamma_second;
						Double_t complete_gamma_deposited = 0.;
						for (Int_t j = 0.;j<entries_califa;j++){
							califahitdata[j] = (R3BCalifaHitData*)CalifaHitData->At(j);
							psi_CALIFA = califahitdata[j]->GetPhi();
							theta_CALIFA = califahitdata[j]->GetTheta();
							E_Califa = califahitdata[j]->GetEnergy();
							time_CALIFA = califahitdata[j]->GetTime();
							vpsi.push_back(psi_CALIFA);
							vtheta.push_back(theta_CALIFA);
							vE.push_back(E_Califa);
							v_origin_E.push_back(E_Califa);
							vtime.push_back(time_CALIFA);
							}
						sort(vE.begin(),vE.end());
						if (vE[vE.size()-1] >30000 && vE[vE.size()-2] > 30000)
						{
						for (Int_t k =0;k< vE.size();k++)
						{
							if (v_origin_E[k] == vE[vE.size()-1]) {
							ind_high_energy = k;
							}
							else if (v_origin_E[k] == vE[vE.size()-2]){
							ind_second_high_energy = k;
							}
							
							else if(v_origin_E[k] < 10000)
							{
							theta_gamma = vtheta[k];
							Double_t energy_gamma = (v_origin_E[k]/1000)*(1/(sqrt(1-beta_precise*beta_precise)))*(1-beta_precise*cos(theta_gamma));
							//here I just want to take the gamma_energy with highest energy, not all gammas. Therefore I store all doppler-corrected energies in a vector...
							v_E_gamma.push_back(energy_gamma);
							h1_gamma_energyE_11B->Fill(energy_gamma);
							complete_gamma_deposited +=v_origin_E[k]/1000;
							if (vpsi[k]< 0.){
							h2_gamma_energyE_psi_11B->Fill(energy_gamma,((2*PI+vpsi[k])/PI)*180);
							}
							if (vpsi[k] > 0.){
							h2_gamma_energyE_psi_11B->Fill(energy_gamma,(vpsi[k]/PI)*180);
	
							}
							}
						

						
						}
						Double_t E1 = v_origin_E[ind_high_energy]/1000.;
						Double_t E2 = v_origin_E[ind_second_high_energy]/1000.;
						Double_t theta1 = vtheta[ind_high_energy];
						Double_t theta2 = vtheta[ind_second_high_energy];
						cout << "------------------------------------------------------" << endl;
						cout << "phi1 = \t" << (vpsi[ind_high_energy]/PI)*180. << endl;
						cout << "phi2 = \t" << (vpsi[ind_second_high_energy]/PI)*180. << endl;
						cout << "------------------------------------------------------" << endl;
						if (vpsi[ind_high_energy] < 0.){
						psi1 = PI - vpsi[ind_high_energy];	
						psi1_2pi = 2*PI+vpsi[ind_high_energy];
						}
						if (vpsi[ind_high_energy] > 0.){
						psi1_2pi = vpsi[ind_high_energy];
						}
						if (vpsi[ind_second_high_energy] < 0.){
						psi2 = PI - vpsi[ind_second_high_energy];
						psi2_2pi = 2*PI+vpsi[ind_second_high_energy];
						}
						if (vpsi[ind_second_high_energy] > 0.){
						psi2_2pi = vpsi[ind_second_high_energy];
						}
						
						psi1_2pi_degr = (psi1_2pi/PI)*180.;
						psi2_2pi_degr = (psi2_2pi/PI)*180.; 
						theta1_degr = (theta1/PI)*180.;
						theta2_degr = (theta2/PI)*180.;

						
						Double_t mom_p1 = sqrt(pow((E1+938.272),2)-pow(938.272,2));
                        Double_t mom_p2 = sqrt(pow((E2+938.272),2)-pow(938.272,2));
						//here I do my calculation for the excitation energy:
						Double_t beta_1 = mom_p1/(sqrt(mom_p1*mom_p1+938.272*938.272));
						Double_t gamma_1 = 1/(sqrt(1-beta_1*beta_1));
						Double_t beta_2 = mom_p2/(sqrt(mom_p2*mom_p2+938.272*938.272));						
						Double_t gamma_2 = 1/(sqrt(1-beta_2*beta_2));
						Double_t angle_p1_11b = abs(atan(tan(theta1)*cos(psi1_2pi) )-psi_in);
						Double_t angle_p2_11b = abs(atan(tan(theta2)*cos(psi2_2pi) )-psi_in);
						Double_t angle_p2_p1 = abs(atan(tan(theta1)*cos(psi1_2pi) )-atan(tan(theta2)*cos(psi2_2pi)));
						Double_t exc_energy = sqrt(2*938.272*938.272+(11*931.5)*(11*931.5)+2*gamma_1*gamma_2*938.272*938.272*(1-beta_1*beta_2*cos(angle_p2_p1))+gamma_1*gamma_precise*(11*931.5)*938.272*(1-beta_precise*beta_1*cos(angle_p1_11b))+gamma_2*gamma_precise*(11*931.5)*938.272*(1-beta_precise*beta_2*cos(angle_p2_11b)) + 938.272*0.714549*1.42942*938.272*0.714549*1.42942)+complete_gamma_deposited-(12*931.5)-938.272-1.42942*938;
						cout << "exec energy\t "<< exc_energy << endl;
						h1_psi_in_11B->Fill(psi_in);
						//

						TVector2 k_0 =   TVector2(11414.5,0);
                        TVector2 k_1 =   TVector2(mom_p1*cos(theta1),mom_p1*sin(theta1));
                        TVector2 k_2 =   TVector2(mom_p2*cos(theta2),mom_p2*sin(theta2));
                        TVector2 k_a_1 = TVector2(mom_11b*cos(psi_in),mom_11b*sin(psi_in));
                        Double_t cos_distr = ((k_0 - k_1 - k_2)*k_a_1)/(sqrt((k_0 - k_1 - k_2)*(k_0 - k_1 - k_2))*sqrt(k_a_1*k_a_1));

                        //CALIFA histos..
                        if (v_E_gamma.size()){
                        h1_gamma_energyE_max_val_11B->Fill(*max_element(v_E_gamma.begin(),v_E_gamma.end()));
						if(v_E_gamma.size() > 1){
						sort(v_E_gamma.begin(),v_E_gamma.end());
						h2_gamma_fst_vs_snd_11B->Fill(v_E_gamma[v_E_gamma.size()-2],v_E_gamma[v_E_gamma.size()-1]);
						}
						}
                        h1_cos_gamma_11b_p_i->Fill(cos_distr);
                        //cout << "cos_distr:\t" << cos_distr << endl;
                        h1_theta_1_plus_theta_2_CALIFA_11b->Fill(((theta1+theta2)/PI)*180.);
                        h2_E1_vs_E2_CALIFA_11b->Fill(E1,E2);
                        h1_E1_plus_E2_CALIFA_11b->Fill(E1+E2);
                        h2_long_mom_p1p2_long_mom11b->Fill(mom_p1*cos(theta1)+mom_p2*cos(theta2),mom_11b*cos(psi_in));

                        h1_transv_mom_difference_p1_p2_11b->Fill(abs(mom_p1*sin(theta1)-mom_p2*sin(theta2)));
                        h2_E1_plus_E2_CALIFA_vs_full_mom11b->Fill(E1+E2,mom_11b);
                        h1_binding_energy_11b->Fill(11*400-(E1+E2+sqrt(pow(mom_11b,2)+pow(11*938.272,2))-(11*938.272)));
                        h2_theta1_theta2_11b->Fill((theta1/PI)*180.,(theta2/PI)*180.);
                        h2_psi1_psi2_11b->Fill((psi1/PI)*180.,(psi2/PI)*180.);
						h2_psi1_psi2_11b_2pi->Fill((psi1_2pi/PI)*180.,(psi2_2pi/PI)*180.);
						h1_theta_1_CALIFA_11b->Fill((theta1/PI)*180.);
						h1_theta_2_CALIFA_11b->Fill((theta2/PI)*180.);
						h2_theta_sum_vs_phi_diff_11b->Fill(theta1_degr+theta2_degr,abs(psi1_2pi_degr-psi2_2pi_degr));

						//for my analysis in y direction..
						//Double_t y1 = yMW1+yMW1_shift;
        				//Double_t y2 = -(yMW2+yMW2_shift);
        				//Double_t y3 = yMW3+yMW3_shift;
						//Double_t angle_before_glad_y = atan((y2-y1)/(zM2-zM1));
						//Double_t mw1_to_mw3_z = pathlength_precise -(760./cos(psi_out-(PI/10.))) - (1408-445);
						//Double_t angle_after_glad_y = atan((y3-y1)/mw1_to_mw3_z);
						//cout << "distance between MW1 and MW2:\t" << zM2-zM1 << endl;
						//cout << "y1:\t" << y1 << endl;
						//cout << "y2:\t" << y2 << endl;
						//cout << "this is mw2.fy:\t" << yMW2 << "------------------------<<<" << endl;
						//cout << "y difference y2 and y1:\t" << y2-y1 << endl;
						//cout << "angle before glad:\t" << (angle_before_glad_y/PI)*180. << endl;
						//cout << "angle after glad:\t" << (angle_after_glad_y/PI)*180. << endl;
						//cout << "angle diff before GLAD - after GLAD:\t" << ((angle_before_glad_y-angle_after_glad_y)/PI)*180. << endl;
                        //looking at the cut triangle
                        if (-0.9973*(mom_p1*cos(theta1)+mom_p2*cos(theta2))+10913 > (mom_11b*cos(psi_in)) && (mom_11b*cos(psi_in)) > 9820 && (mom_p1*cos(theta1)+mom_p2*cos(theta2)) > 750) {
                        h2_theta_1_vs_theta_2_CALIFA_11b_cut->Fill((theta1/PI)*180.,(theta2/PI)*180.);
                        h2_E1_vs_E2_CALIFA_11b_cut->Fill(E1,E2);
						h1_psi_in_11B_cut_triangle->Fill(psi_in);
						//h1_angl_diff_y_cut_triangle->Fill((abs(angle_after_glad_y-angle_before_glad_y)/PI)*180.);

						}
						else {
							//h1_angl_diff_y->Fill((abs(angle_after_glad_y-angle_before_glad_y)/PI)*180.);
						}
						//here I try to plot cos(p_i vs p_11b) in the 12C frame...
						TVector3 p_11b (11*931.5*gamma_precise*beta_precise*sin(psi_in),0.,-gamma_given*0.714549*11*931.5*gamma_precise+11*931.5*gamma_precise*beta_precise*cos(psi_in)*gamma_given);
                        TVector3 p_p1(938.272*gamma_1*beta_1*sin(theta1)*cos(psi1_2pi),0.,-gamma_given*0.714549*938.272*gamma_1+gamma_given*938.272*gamma_1*beta_1*cos(theta1));
                        TVector3 p_p2(938.272*gamma_2*beta_2*sin(theta2)*cos(psi2_2pi),0.,-gamma_given*0.714549*938.272*gamma_2+gamma_given*938.272*gamma_2*beta_2*cos(theta2));
                        TVector3 p_tar(0.,0.,-gamma_given*0.714549*938.272);
                        TVector3 p_initial(p_p1+p_p2-p_tar);
						cout << "this is momentum of initial projectile\t" << p_initial.Mag() << endl;
						h1_p_missing_11B->Fill(p_initial.Mag());

						//here I will do my analysis with phi angle sweeping...+-2.8degrees = 0.048869219056 rad------------
						Double_t my_sweep = 0.048869219056;
						vector<double> cos_sweep;
						for (Int_t sweep_1 = -1; sweep_1 < 2; sweep_1++){
							for (Int_t sweep_2 = -1; sweep_2 < 2; sweep_2++){
								
						TVector3 p_11b_sweep (11*931.5*gamma_precise*beta_precise*sin(psi_in),0.,-gamma_given*0.714549*11*931.5*gamma_precise+11*931.5*gamma_precise*beta_precise*cos(psi_in)*gamma_given);
                        TVector3 p_p1_sweep(938.272*gamma_1*beta_1*sin(theta1)*cos(psi1_2pi+sweep_1*my_sweep),0.,-gamma_given*0.714549*938.272*gamma_1+gamma_given*938.272*gamma_1*beta_1*cos(theta1));
                        TVector3 p_p2_sweep(938.272*gamma_2*beta_2*sin(theta2)*cos(psi2_2pi+sweep_2*my_sweep),0.,-gamma_given*0.714549*938.272*gamma_2+gamma_given*938.272*gamma_2*beta_2*cos(theta2));
                        TVector3 p_tar_sweep(0.,0.,-gamma_given*0.714549*938.272);
                        TVector3 p_initial_sweep(p_p1_sweep+p_p2_sweep-p_tar_sweep);
						Double_t angle_cms_sweep = p_initial_sweep.Angle(p_11b_sweep);
						cos_sweep.push_back(cos(angle_cms_sweep));
							}


						}	
						Double_t min_cos_angle_cms = *min_element(cos_sweep.begin(),cos_sweep.end());
						cos_sweep.clear();
						//----------------------------------------------------------------
                        Double_t angle_cms = p_initial.Angle(p_11b);	
						cout << "cos(gamma) questioning....." << endl;	
						if (cos(angle_cms) < -0.6){
						h2_long_mom_p1p2_long_mom11b_low->Fill(mom_p1*cos(theta1)+mom_p2*cos(theta2),mom_11b*cos(psi_in));	
						}
						else if (cos(angle_cms)>-0.6){
						h2_long_mom_p1p2_long_mom11b_high->Fill(mom_p1*cos(theta1)+mom_p2*cos(theta2),mom_11b*cos(psi_in));
						}
						if ((theta1_degr+theta2_degr) < 84 ){
						h2_long_mom_p1p2_long_mom11b_small_opening->Fill(mom_p1*cos(theta1)+mom_p2*cos(theta2),mom_11b*cos(psi_in));
						}
						else if ((theta1_degr+theta2_degr) > 84 ){
						h2_long_mom_p1p2_long_mom11b_large_opening->Fill(mom_p1*cos(theta1)+mom_p2*cos(theta2),mom_11b*cos(psi_in));
                        }
						cout << "asking for angular restrictions...." << endl;
						if ((theta1_degr+theta2_degr) < 84 && abs(180-abs(psi1_2pi_degr-psi2_2pi_degr)) < 30){
					
							h2_long_mom_p1p2_long_mom11b_theta_phi_constr->Fill(mom_p1*cos(theta1)+mom_p2*cos(theta2),mom_11b*cos(psi_in));	
							if (min_cos_angle_cms < -0.6){
								h2_long_mom_p1p2_long_mom11b_low_optimized->Fill(mom_p1*cos(theta1)+mom_p2*cos(theta2),mom_11b*cos(psi_in));	
								if (v_E_gamma.size()){
									h1_cos_low_gamma->Fill(*max_element(v_E_gamma.begin(),v_E_gamma.end()));	
								}
							}
							if (min_cos_angle_cms > 0.6){
								h2_long_mom_p1p2_long_mom11b_high_optimized->Fill(mom_p1*cos(theta1)+mom_p2*cos(theta2),mom_11b*cos(psi_in));
							}

							if (min_cos_angle_cms < -0.8){
								h2_long_mom_p1p2_long_mom11b_low_optimized_08->Fill(mom_p1*cos(theta1)+mom_p2*cos(theta2),mom_11b*cos(psi_in));
							}
							//triangle cut vs straight anti-correlated line....
							if (-0.9973*(mom_p1*cos(theta1)+mom_p2*cos(theta2))+10913 > (mom_11b*cos(psi_in)) && (mom_11b*cos(psi_in)) > 9820 && (mom_p1*cos(theta1)+mom_p2*cos(theta2)) > 750) {
							h1_mw3_fy_11b_cut_triangle->Fill(raw_mwpc3);
							h1_mw2_fy_11b_cut_triangle->Fill(yMW2);
							h1_psi_in_11b_cut_triangle->Fill(psi_in);	

							if (v_E_gamma.size()){
							h1_gamma_energyE_max_val_11B_angle_cut->Fill(v_E_gamma[v_E_gamma.size()-1]);
							}
							v_E_gamma.clear();

							}
							else {
							h1_mw3_fy_11b->Fill(raw_mwpc3);
							h1_mw2_fy_11b->Fill(yMW2);
							h1_psi_in_11b->Fill(psi_in);
							if (v_E_gamma.size()){
							h1_gamma_energyE_max_val_11B_angle_cut->Fill(v_E_gamma[v_E_gamma.size()-1]);
							}
							v_E_gamma.clear();
							}
						}
	
						//here I make strict constraints on the angles, so that I don't have borders
						cout << "before border constraints...." << endl;
						cout << "psi1 in degrees:\t" << psi1_2pi_degr<< endl;
						cout << "psi2 in degrees:\t" << psi2_2pi_degr << endl;
						cout << "theta1 in degrees:\t" << theta1_degr << endl;
						cout << "theta2 in degrees:\t" << theta2_degr << endl;
						 if( (psi1_2pi_degr > 84 && psi1_2pi_degr < 96) || (psi1_2pi_degr > 264 && psi1_2pi_degr < 276) || (psi2_2pi_degr > 84 && psi2_2pi_degr < 96) || (psi2_2pi_degr > 264 && psi2_2pi_degr)){
							cout << "hello in first condition" << endl;
							continue;
							}
						if((psi1_2pi_degr> 42 &&  psi1_2pi_degr < 160.5 && theta1_degr < 46) || (psi1_2pi_degr> 199.5 &&  psi1_2pi_degr < 318 && theta1_degr < 46)){
							cout << "hello in second contidtion" << endl;
							continue;
							}
						if((psi2_2pi_degr> 42 &&  psi2_2pi_degr < 160.5 && theta2_degr < 46) || (psi2_2pi_degr> 199.5 &&  psi2_2pi_degr < 318 && theta2_degr < 46)){
							cout << "hello in third condition " << endl;
							continue;
							}
						if (theta1_degr < 25 || theta2_degr < 25){
							cout << "hello in forth condition" << endl;
							continue;
							}
						cout << "After border constraints ...." << endl;
						h2_long_mom_p1p2_long_mom11b_no_border->Fill(mom_p1*cos(theta1)+mom_p2*cos(theta2),mom_11b*cos(psi_in));
						h1_cos_gamma_11b_p_i_cms->Fill(cos(angle_cms));
						h1_cos_gamma_11b_p_i_cms_optimized->Fill(min_cos_angle_cms);
						cout << "after filling some histograms..." << endl;	
						if (cos(angle_cms) < -0.6){
                        h2_long_mom_p1p2_long_mom11b_no_border_low->Fill(mom_p1*cos(theta1)+mom_p2*cos(theta2),mom_11b*cos(psi_in));
                        }
                        else if (cos(angle_cms)>0.6){
                        h2_long_mom_p1p2_long_mom11b_no_border_high->Fill(mom_p1*cos(theta1)+mom_p2*cos(theta2),mom_11b*cos(psi_in));
                        }
						cout << "After asking cos(gamma) <> -0.6...."<< endl;

					
						//some calculations without border....
						cout << "last calculation....." << endl;
						for (Int_t k =0;k< vE.size();k++)
                        {
                            if(v_origin_E[k] < 10000)
                            {
                            Double_t theta_gamma_noborder = vtheta[k];
                            Double_t energy_gamma = (v_origin_E[k]/1000)*(1/(sqrt(1-beta*beta)))*(1-beta*cos(theta_gamma_noborder));
                            h1_gamma_energyE_11B_no_border->Fill(energy_gamma);
                            
                            }



                        }
						
						}
						//End of CALIFA handling---------------------------------------------------	
						//                                                
						
						}						
						if (a_q < aq_cut_z_5 && a_q > (aq_cut_z_5-0.15) && entries_califa >= 2){  //10B
						psi_in = mean_psi_out_10b-slope_10b*xMW3 - angle_offs_10b;
						m = (cos(psi_out)-cos(psi_in))/(sin(psi_out)+sin(psi_in));
						rho =  (Leff/(2*sin((psi_in+psi_out)/2)))*sqrt(pow(((m+tan(alpha_G))/(1-m*tan(alpha_G))),2)+1);
                        w = 2*abs(asin(BD/(2*rho)));
						Double_t pathlength_precise = abs((zB-zT)/(cos(psi_in)))+rho*w+abs((zToFW-z_D)/cos(psi_out));
                        Double_t beta_precise = ((pathlength_precise/time_target_tof)*pow(10,6))/light_c;
                        Double_t a_q_precise = (pow(10,-3)*((((mag_field*rho)/(beta_precise))*1.602176634*pow(10,-19))/(1.66053906660*pow(10,-27)*light_c)))/gamma_given;
						Double_t mom_10b = a_q_precise*5*(beta/(sqrt(1-(beta*beta))))*938.272;
						//cout << "Long mom 10b:\t" << mom_10b*cos(psi_in) << endl;

						//CALIFA handling
						Double_t psi_CALIFA;
                        Double_t theta_CALIFA;
						Double_t time_CALIFA;
                        Double_t E_Califa;
                        vector<double> vE1;
                        vector<double> vtheta1CALIFA;
                        vector<double> vpsi1CALIFA;
                        vector<double> vE2;
                        vector<double> vtheta2CALIFA;
                        vector<double> vpsi2CALIFA;

						//TJ do the vector handling in a different way...
                        vector<double> vpsi;
                        vector<double> vtheta;
                        vector<double> vE;
                        vector<double> vtime;
                        vector<double> v_origin_E;
                        Double_t ind_high_energy;
                        Double_t ind_second_high_energy;
                        Double_t theta_gamma;
						Double_t psi1;
						Double_t psi2;
                        for (Int_t j = 0.;j<entries_califa;j++){
                            califahitdata[j] = (R3BCalifaHitData*)CalifaHitData->At(j);
                            psi_CALIFA = califahitdata[j]->GetPhi();
                            theta_CALIFA = califahitdata[j]->GetTheta();
                            E_Califa = califahitdata[j]->GetEnergy();
                            time_CALIFA = califahitdata[j]->GetTime();
                            vpsi.push_back(psi_CALIFA);
                            vtheta.push_back(theta_CALIFA);
                            vE.push_back(E_Califa);
                            v_origin_E.push_back(E_Califa);
                            vtime.push_back(time_CALIFA);
                            }
                        sort(vE.begin(),vE.end());
                        if (vE[vE.size()-1] >30000 && vE[vE.size()-2] > 30000)
                        {
                        for (Int_t k =0;k<vE.size();k++)
                        {
                            if (v_origin_E[k] == vE[vE.size()-1]) {
                            ind_high_energy = k;
                            }
                            else if (v_origin_E[k] == vE[vE.size()-2]){
                            ind_second_high_energy = k;
                            }
                            else if(v_origin_E[k] < 10000)
                            {
                            theta_gamma = vtheta[k];
                            Double_t energy_gamma = (v_origin_E[k]/1000)*(1/(sqrt(1-beta*beta)))*(1-beta*cos(theta_gamma));
                            h1_gamma_energyE_10B->Fill(energy_gamma);

                            }



                        }	
						Double_t E1 = v_origin_E[ind_high_energy]/1000.;
                        Double_t E2 = v_origin_E[ind_second_high_energy]/1000.;
                        Double_t theta1 = vtheta[ind_high_energy];
                        Double_t theta2 = vtheta[ind_second_high_energy];
		
						if (vpsi[ind_high_energy] < 0.){
                        psi1 = 2*PI + vpsi[ind_high_energy];
                        }   
						if (vpsi[ind_high_energy] > 0.){
						psi1 = vpsi[ind_high_energy];
						}
                        if (vpsi[ind_second_high_energy] < 0.){
                        psi2 = 2*PI + vpsi[ind_second_high_energy];
                        }
						if (vpsi[ind_second_high_energy] > 0.){
						psi2 = vpsi[ind_second_high_energy];
						}
						
                        Double_t mom_p1 = sqrt(pow((E1+938.272),2)-pow(938.272,2));
                        Double_t mom_p2 = sqrt(pow((E2+938.272),2)-pow(938.272,2));
						h1_theta_1_plus_theta_2_CALIFA_10b->Fill(((theta1+theta2)/PI)*180.);
                        h2_E1_vs_E2_CALIFA_10b->Fill(E1,E2);
                        h1_E1_plus_E2_CALIFA_10b->Fill(E1+E2);
                        h2_long_mom_p1p2_long_mom10b->Fill(mom_p1*cos(theta1)+mom_p2*cos(theta2),mom_10b*cos(psi_in));
                        h1_transv_mom_difference_p1_p2_10b->Fill(abs(mom_p1*sin(theta1)-mom_p2*sin(theta2)));
						}
					//	//End of CALIFA handling--------------------------------------------

						//histos
                        h2_z_vs_a_q->Fill(a_q_precise,charge_val);
						h2_tof_vs_aq_fix_g_10b->Fill(a_q_precise,time_target_tof);
                        }
						}
						if (charge_val< charge_cut+0.8 && charge_val> charge_cut) //Z = 6
						{
						if (a_q > aq_cut_z_6 && a_q < (aq_cut_z_6+2*width_12c)){  //12C
						psi_in = mean_psi_out_12c-slope_12c*xMW3 - angle_offs_12c;
						m = (cos(psi_out)-cos(psi_in))/(sin(psi_out)+sin(psi_in));
						rho =  (Leff/(2*sin((psi_in+psi_out)/2)))*sqrt(pow(((m+tan(alpha_G))/(1-m*tan(alpha_G))),2)+1);
                        w = 2*abs(asin(BD/(2*rho)));
						Double_t pathlength_precise = abs((zB-zT)/(cos(psi_in)))+rho*w+abs((zToFW-z_D)/cos(psi_out));
                        Double_t beta_precise = ((pathlength_precise/time_target_tof)*pow(10,6))/light_c;
                        Double_t a_q_precise = (pow(10,-3)*((((mag_field*rho)/(beta_precise))*1.602176634*pow(10,-19))/(1.66053906660*pow(10,-27)*light_c)))/gamma_given;
						if (time_target_tof < 35.13){
						h1_tof_detnr_strange_12c->Fill(i+1);							
						h1_mw3_fy_strange_12c->Fill(raw_mwpc3);
						h2_tof_path_strange12c->Fill(pathlength_precise,time_target_tof);
						h2_radius_path_strange12c->Fill(rho,time_target_tof);
						
						}
						if (time_target_tof> 35.35 && time_target_tof<35.61){
						h1_tof_detnr_12c->Fill(i+1);	
						h1_mw3_fy_12c->Fill(raw_mwpc3);
						h2_tof_path_12c->Fill(pathlength_precise,time_target_tof);
						h2_radius_path_12c->Fill(rho,time_target_tof);
						}
						if (a_q_precise > 1.997){
						h1_tof_detnr_strange_12c_large_aq->Fill(i+1);	
						h1_mw3_fy_strange_12c_large_aq->Fill(raw_mwpc3);
						h2_tof_path_strange_12c_large_aq->Fill(pathlength_precise,time_target_tof);
						h2_radius_path_strange_12c_large_aq->Fill(rho,time_target_tof);
						}
						//histos
                        h2_z_vs_a_q->Fill(a_q_precise,charge_val);
						h2_tof_vs_aq_fix_g_12c->Fill(a_q_precise,time_target_tof);
						h2_mw2x_vs_tof_12c->Fill(xMW2,time_target_tof);
						h2_mw3x_vs_tof_12c->Fill(xMW3,time_target_tof);
						h2_charge_vs_tof_12c->Fill(charge_val,time_target_tof);
						h2_psi_in_vs_tof_12c->Fill(psi_in,time_target_tof);
						h2_detnr_vs_tof_12c->Fill(i+1,time_target_tof);
                        }
						if (a_q < aq_cut_z_6-0.06 && a_q > (aq_cut_z_6-1.65*width_11c) && entries_califa >= 2){  //11C before it was a_q < aq_cut_z_6 && a_q > (aq_cut_z_6-2*width_11c), too bad for 11c...
						psi_in = mean_psi_out_11c-slope_11c*xMW3 - angle_offs_11c;
						m = (cos(psi_out)-cos(psi_in))/(sin(psi_out)+sin(psi_in));
						rho =  (Leff/(2*sin((psi_in+psi_out)/2)))*sqrt(pow(((m+tan(alpha_G))/(1-m*tan(alpha_G))),2)+1);
                        w = 2*abs(asin(BD/(2*rho)));
						Double_t pathlength_precise = abs((zB-zT)/(cos(psi_in)))+rho*w+abs((zToFW-z_D)/cos(psi_out));
                        Double_t beta_precise = ((pathlength_precise/time_target_tof)*pow(10,6))/light_c;
                        Double_t a_q_precise = (pow(10,-3)*((((mag_field*rho)/(beta_precise))*1.602176634*pow(10,-19))/(1.66053906660*pow(10,-27)*light_c)))/gamma_given;
						//histos
                        h2_z_vs_a_q->Fill(a_q_precise,charge_val);
						h2_tof_vs_aq_fix_g_11c->Fill(a_q_precise,time_target_tof);
                        }

						}
}
}
}
}
}
delete gr1;
delete gr2;
delete gr3;
}
 delete [] sofmwpc0hitdata;   
 delete [] sofmwpc1hitdata;   
 delete [] sofmwpc2hitdata;   
 delete [] sofmwpc3hitdata;   
 delete [] softofwtcaldata;   
 delete [] softofwmappeddata; 

}
delete [] softwimhitdata;
}

TList *l = new TList();
l->Add(h2_z_vs_a_q);
l->Add(h1_a_q_z6);
l->Add(h1_a_q_z5);
l->Add(h1_gamma_energyE_10B);
l->Add(h1_gamma_energyE_11B);
l->Add(h1_theta_1_plus_theta_2_CALIFA_10b);
l->Add(h1_theta_1_plus_theta_2_CALIFA_11b);
l->Add(h2_E1_vs_E2_CALIFA_10b);
l->Add(h2_E1_vs_E2_CALIFA_11b);
l->Add(h1_E1_plus_E2_CALIFA_10b);
l->Add(h1_E1_plus_E2_CALIFA_11b);
l->Add(h2_long_mom_p1p2_long_mom11b);
l->Add(h2_long_mom_p1p2_long_mom10b);
l->Add(h1_transv_mom_difference_p1_p2_10b);
l->Add(h1_transv_mom_difference_p1_p2_11b);
l->Add(h2_tof_vs_aq_fix_g_11b);
l->Add(h2_tof_vs_aq_fix_g_10b);
l->Add(h2_tof_vs_aq_fix_g_11c);
l->Add(h2_tof_vs_aq_fix_g_12c);
l->Add(h2_theta_out_vs_mw3_10b);
l->Add(h2_theta_out_vs_mw3_11b);
l->Add(h2_theta_out_vs_mw3_11c);
l->Add(h2_theta_out_vs_mw3_12c);
l->Add(h2_E1_plus_E2_CALIFA_vs_full_mom11b);
l->Add(h2_theta_1_vs_theta_2_CALIFA_11b_cut);
l->Add(h2_E1_vs_E2_CALIFA_11b_cut);
l->Add(h1_binding_energy_11b);
l->Add(h1_cos_gamma_11b_p_i);
l->Add(h1_time_diff_gamma_11b);
l->Add(h2_theta1_theta2_11b);
l->Add(h2_psi1_psi2_11b);
l->Add(h2_psi1_psi2_11b_2pi);
l->Add(h1_tof_detnr_strange_12c);
l->Add(h1_tof_detnr_strange_12c_large_aq);
l->Add(h1_tof_detnr_12c);
l->Add(h1_mw3_fy_strange_12c);
l->Add(h1_mw3_fy_12c);
l->Add(h1_mw3_fy_strange_12c_large_aq);
l->Add(h2_z_vs_a_q_nocut);
l->Add(h1_theta_1_CALIFA_11b);
l->Add(h1_theta_2_CALIFA_11b);
l->Add(h2_gamma_energyE_psi_11B);
l->Add(h2_long_mom_p1p2_long_mom11b_no_border);
l->Add(h1_cos_gamma_11b_p_i_cms);
l->Add(h1_cos_gamma_11b_p_i_cms_optimized);
l->Add(h1_gamma_energyE_11B_no_border);
l->Add(h1_psi_in_11B_cut_triangle);
l->Add(h1_psi_in_11B);
l->Add(h2_long_mom_p1p2_long_mom11b_no_border_high);
l->Add(h2_long_mom_p1p2_long_mom11b_no_border_low);
l->Add(h2_long_mom_p1p2_long_mom11b_high);
l->Add(h2_long_mom_p1p2_long_mom11b_low);
l->Add(h1_gamma_energyE_max_val_11B);
l->Add(h2_gamma_fst_vs_snd_11B);
l->Add(h2_long_mom_p1p2_long_mom11b_small_opening);
l->Add(h2_long_mom_p1p2_long_mom11b_theta_phi_constr);
l->Add(h2_long_mom_p1p2_long_mom11b_large_opening);
l->Add(h1_gamma_energyE_max_val_11B_angle_cut);
//---------------------------------------------------
l->Add(h1_mw3_fy_11b_cut_triangle);
l->Add(h1_mw2_fy_11b_cut_triangle);
l->Add(h1_psi_in_11b_cut_triangle);
l->Add(h1_mw3_fy_11b);
l->Add(h1_mw2_fy_11b);
l->Add(h1_psi_in_11b);

l->Add(h2_tof_path_strange12c);
l->Add(h2_tof_path_strange_12c_large_aq);
l->Add(h2_tof_path_12c);
l->Add(h2_radius_path_strange12c);
l->Add(h2_radius_path_strange_12c_large_aq);
l->Add(h2_radius_path_12c);
l->Add(h2_mw2x_vs_tof_12c);
l->Add(h2_mw3x_vs_tof_12c);
l->Add(h2_charge_vs_tof_12c);
l->Add(h2_psi_in_vs_tof_12c);
l->Add(h2_detnr_vs_tof_12c);
//l->Add(h1_angl_diff_y_cut_triangle);
//l->Add(h1_angl_diff_y);
l->Add(h2_long_mom_p1p2_long_mom11b_low_optimized);
l->Add(h2_long_mom_p1p2_long_mom11b_high_optimized);
l->Add(h2_long_mom_p1p2_long_mom11b_low_optimized_08);
l->Add(h1_p_missing_11B);
l->Add(h1_cos_low_gamma);
l->Add(h2_theta_sum_vs_phi_diff_11b);

sprintf(f_out_name,"/home/ge37liw/plots/histos/zero_with_miss/phi_test_automatic_src_califa_%s.root",count_i);


TFile *f = new TFile(f_out_name,"RECREATE");
l->Write("histlist", TObject::kSingleKey);
l->Delete();
//write all cut parameters to file
char par_file_name[200];
sprintf(par_file_name,"/home/ge37liw/plots/par_files/parameters_run%s.txt",count_i);
ofstream par_file;
par_file.open(par_file_name);
par_file << " ----------------- CUT PARAMETERS for RUN" << count_i << "----------"  <<endl;
par_file << "---> 10B:\t " << endl;
par_file << "charge:\t" << charge_cut-1 << "< Z < " << charge_cut << endl;
par_file << "A/q:\t" << aq_cut_z_5-0.15 << "< A/q < " <<  aq_cut_z_5 << endl;
par_file << "psi_in = " << mean_psi_out_10b << "+" << -slope_10b << "*xMW3" << "-" << angle_offs_10b << endl;
par_file << "-----------------------------------------------------" << endl;
par_file << "---> 11B:\t " << endl;
par_file << "charge:\t" << charge_cut-1 << "< Z < " << charge_cut << endl;
par_file << "A/q:\t" << aq_cut_z_5 << "< A/q < " << aq_cut_z_5+2*width_11b << endl;
par_file << "psi_in = " << mean_psi_out_11b << "+" << -slope_11b << "*xMW3" << "-" << angle_offs_11b << endl;
par_file << "-----------------------------------------------------" << endl;
par_file << "---> 11C:\t " << endl;
par_file << "charge:\t" << charge_cut << "< Z < " << charge_cut + 0.8 << endl;
par_file << "A/q:\t" << aq_cut_z_6-1.65*width_11c << "< A/q < " << aq_cut_z_6-0.06 << endl;
par_file << "psi_in = " << mean_psi_out_11c << "+" << -slope_11c << "*xMW3" << "-" << angle_offs_11c << endl;
par_file << "-----------------------------------------------------" << endl; 
par_file << "---> 12C:\t " << endl;
par_file << "charge:\t" << charge_cut << "< Z < " << charge_cut + 0.8 << endl;
par_file << "A/q:\t" << aq_cut_z_6 << "< A/q < " << aq_cut_z_6+2*width_12c << endl;
par_file << "psi_in = " << mean_psi_out_12c << "+" << -slope_12c << "*xMW3" << "-" << angle_offs_12c << endl;
par_file << "-----------------------------------------------------" << endl;
par_file.close();

//
f->Close();
cout << "end of process, successful" <<endl;
}