#define cnt_cxx
#include "cnt.h"
#include "TH2.h"
#include "TStyle.h"
#include "TCanvas.h"
#include <multimap>

using std::multimap;

void cnt::Loop()
{

  ofstream ndfile;
  ndfile.open("neutrino.txt", ios::app);

  //initialising fitLengths[500][500] , needed for constructing justin's variable
	
  Int_t fitLengths[500][500];
	
    for (Int_t row=0; row<500; ++row){
    for (Int_t column=0; column<500; ++column){
      fitLengths[row][column] = 0;
    }
  }

  //Iteration 1
  fitLengths[0][0] = 0;
  fitLengths[1][0] = 0;
  fitLengths[2][0] = 0;

  for(Int_t tl=3; tl<8; ++tl){
    Int_t i=0;
    for(Int_t fl=tl; fl>1; --fl, ++i){
      fitLengths[tl][i] = fl;
      if (2==fl){fitLengths[tl][i] = 0;}
    }
  }

  for(Int_t tl=8; tl<11; ++tl){
    Int_t i=0;
    for(Int_t fl = 7; fl>3; --fl,++i){
      fitLengths[tl][i] = fl;
    }
    for(Int_t fl = 8; fl<=tl; ++fl,++i){
      fitLengths[tl][i] = fl;
    }
    fitLengths[tl][i] = 3;
    ++i;
    fitLengths[tl][i] = 0;
  }

  for (Int_t tl=11; tl<500; ++tl){
    fitLengths[tl][0] = 7;
    fitLengths[tl][1] = 8;
    fitLengths[tl][2] = 6;
    fitLengths[tl][3] = 5;
    fitLengths[tl][4] = 4;
    fitLengths[tl][5] = 0;

    //cout<<tl<<"  "<<fitLengths[tl][0]<<endl;
  }
	
  for (Int_t row=0; row<500; ++row){
    for (Int_t column=0; column<500; ++column){

      //ndfile<<row<<"  "<<column<<"  "<<fitLengths[row][column]<<endl;

    }
  }

////////  finished initialising the arrays /////////////////////////////////////

  if (fChain == 0) return;

  Long64_t nentries = fChain->GetEntries();       

  TFile *f2 = new TFile("nocuts_neutrino.root","recreate");
  TTree *neus = new TTree("neus","neus tree");

  TFile *f4 = new TFile("nocuts_antineu.root","recreate");
  TTree *aneus = new TTree("aneus","aneus tree");

  ndfile<<"no of entries "<<nentries<<endl;

  Int_t nbytes = 0, nb = 0;

  Int_t trk_ind;
  Int_t trk_ind_next;
  Int_t trk_ind_big;  
  Int_t trk_ind_little;

  Int_t header_run;
  Int_t header_subrun;
  Int_t header_snarl;

  Int_t   trkl_number;
  Float_t trkl_vtxx;
  Float_t trkl_vtxy;
  Float_t trkl_vtxz;
  Float_t trkl_vtxu;
  Float_t trkl_vtxv;
  Float_t trkl_vtxrad;
  Float_t trkl_endx;
  Float_t trkl_endy;
  Float_t trkl_endz;
  Float_t trkl_endu;
  Float_t trkl_endv;
  Float_t trkl_endrad;

  Int_t trkl_planes;
  Int_t trkl_uplanes;
  Int_t trkl_vplanes;
  Int_t trkl_begplane;
  Int_t trkl_ubegplane;
  Int_t trkl_vbegplane;
  Int_t trkl_endplane; 
  Int_t trkl_uendplane; 
  Int_t trkl_vendplane; 

  Float_t trkl_redchi2;
  Int_t trkl_pass;
  Int_t trkl_nstps;
  Int_t trkl_coilhole;
  Int_t trkl_exitstop;
  Int_t trkl_new_pflag;
  Int_t trkl_oldpos_stps;
  Int_t trkl_oldneg_stps;
  Int_t trkl_pos_stps;
  Int_t trkl_neg_stps;
  Int_t trkl_trklike;

  Float_t trkl_fitmom;
  Float_t trkl_qoverp;
  Float_t trkl_eqp;
  Float_t trkl_rangemom;
  Float_t trkl_energy;
  Float_t trkl_length;
  Float_t trkl_prob;

  Float_t trkl_dcosz;
  Float_t trkl_dcosx;
  Float_t trkl_dcosy;

  Int_t shw_ind;
  Int_t shwl_number;
  Float_t shwl_energy;
  Float_t shwl_vtxx;
  Float_t shwl_vtxy;
  Float_t shwl_vtxz;
  Float_t shwl_vtxu;
  Float_t shwl_vtxv;
  Int_t shwl_planes;
  Int_t shwl_uplanes;
  Int_t shwl_vplanes;
  Int_t shwl_begplane;
  Int_t shwl_ubegplane;
  Int_t shwl_vbegplane;
  Int_t shwl_endplane; 
  Int_t shwl_uendplane; 
  Int_t shwl_vendplane; 

  Float_t shwl_dcosz;
  Float_t shwl_dcosx;
  Float_t shwl_dcosy;
  Int_t shwl_count;

  Int_t  evt_number;
  Int_t  evt_ntracks;
  Int_t  evt_nshowers;

  Float_t evt_vtxx;
  Float_t evt_vtxy;
  Float_t evt_vtxz;
  Float_t evt_vtxu;
  Float_t evt_vtxv;
  Float_t evt_endx;
  Float_t evt_endy;
  Float_t evt_endz;
  Float_t evt_endu;
  Float_t evt_endv;

  Float_t evt_neu_energy;
  Float_t evt_y;
  Float_t evt_x;
  Float_t evt_q2;
  Float_t evt_pid;

  //MC variables

  Int_t mcl_index;

  Float_t mcl_mu_energy;
  Float_t mcl_mu_px;
  Float_t mcl_mu_py;
  Float_t mcl_mu_pz;

  Float_t mcl_shw_energy;
  Float_t mcl_shw_px;
  Float_t mcl_shw_py;
  Float_t mcl_shw_pz;

  Float_t mcl_nu_energy;
  Float_t mcl_nu_px;
  Float_t mcl_nu_py;
  Float_t mcl_nu_pz;

  Float_t mcl_tar_energy;
  Float_t mcl_tar_px;
  Float_t mcl_tar_py;
   Float_t mcl_tar_pz;

  Float_t   mcl_a;
  Int_t     mcl_inu;
  Int_t     mcl_iaction;
  Int_t     mcl_ires;
  Float_t   mcl_x;
  Float_t   mcl_y;
  Float_t   mcl_q2;
  Float_t   mcl_w2;
  Int_t     mcl_initstate;
  Int_t     mcl_nucleus;
  Int_t     mcl_hadfs;
  Int_t     mc_mul;

  //DECLARING THE NEUTRINO TREE  

  neus->Branch("header_run",&header_run,"header_run/I");
  neus->Branch("header_subrun",&header_subrun,"header_subrun/I");
  neus->Branch("header_snarl",&header_snarl,"header_snarl/I");

  neus->Branch("trkl_vtxx",&trkl_vtxx,"trkl_vtxx/F");
  neus->Branch("trkl_vtxy",&trkl_vtxy,"trkl_vtxy/F");
  neus->Branch("trkl_vtxz",&trkl_vtxz,"trkl_vtxz/F");
  neus->Branch("trkl_vtxu",&trkl_vtxu,"trkl_vtxu/F");
  neus->Branch("trkl_vtxv",&trkl_vtxv,"trkl_vtxv/F");
  neus->Branch("trkl_vtxrad",&trkl_vtxrad,"trkl_vtxrad/F");
  neus->Branch("trkl_endx",&trkl_endx,"trkl_endx/F");
  neus->Branch("trkl_endy",&trkl_endy,"trkl_endy/F");
  neus->Branch("trkl_endz",&trkl_endz,"trkl_endz/F");
  neus->Branch("trkl_endu",&trkl_endu,"trkl_endu/F");
  neus->Branch("trkl_endv",&trkl_endv,"trkl_endv/F");
  neus->Branch("trkl_endrad",&trkl_endrad,"trkl_endrad/F");

  neus->Branch("trkl_number",&trkl_number,"trkl_number/I");
  neus->Branch("trkl_planes",&trkl_planes,"trkl_planes/I");
  neus->Branch("trkl_uplanes",&trkl_uplanes,"trkl_uplanes/I");
  neus->Branch("trkl_vplanes",&trkl_vplanes,"trkl_vplanes/I");
  neus->Branch("trkl_begplane",&trkl_begplane,"trkl_begplane/I");
  neus->Branch("trkl_ubegplane",&trkl_ubegplane,"trkl_ubegplane/I");
  neus->Branch("trkl_vbegplane",&trkl_vbegplane,"trkl_vbegplane/I");
  neus->Branch("trkl_endplane",&trkl_endplane,"trkl_endplane/I");
  neus->Branch("trkl_uendplane",&trkl_uendplane,"trkl_uendplane/I");
  neus->Branch("trkl_vendplane",&trkl_vendplane,"trkl_vendplane/I");

  neus->Branch("trkl_redchi2",&trkl_redchi2,"trkl_redchi2/F");
  neus->Branch("trkl_pass",&trkl_pass,"trkl_pass/I");
  neus->Branch("trkl_nstps",&trkl_nstps,"trkl_nstps/I");
  neus->Branch("trkl_coilhole",&trkl_coilhole,"trkl_coilhole/I");
  neus->Branch("trkl_exitstop",&trkl_exitstop,"trkl_exitstop/I");
  neus->Branch("trkl_new_pflag",&trkl_new_pflag,"trkl_new_pflag/I");
  neus->Branch("trkl_trklike",&trkl_trklike,"trkl_trklike/I");
  neus->Branch("trkl_oldpos_stps",&trkl_oldpos_stps,"trkl_oldpos_stps/I");
  neus->Branch("trkl_oldneg_stps",&trkl_oldneg_stps,"trkl_oldneg_stps/I");
  neus->Branch("trkl_pos_stps",&trkl_pos_stps,"trkl_pos_stps/I");
  neus->Branch("trkl_neg_stps",&trkl_neg_stps,"trkl_neg_stps/I");

  neus->Branch("trkl_fitmom",&trkl_fitmom,"trkl_fitmom/F");
  neus->Branch("trkl_qoverp",&trkl_qoverp,"trkl_qoverp/F");
  neus->Branch("trkl_eqp",&trkl_eqp,"trkl_eqp/F");
  neus->Branch("trkl_rangemom",&trkl_rangemom,"trkl_rangemom/F");
  neus->Branch("trkl_energy",&trkl_energy,"trkl_energy/F");
  neus->Branch("trkl_length",&trkl_length,"trkl_length/F");
  neus->Branch("trkl_prob",&trkl_prob,"trkl_prob/F");

  neus->Branch("trkl_dcosz",&trkl_dcosz,"trkl_dcosz/F");
  neus->Branch("trkl_dcosx",&trkl_dcosx,"trkl_dcosx/F");
  neus->Branch("trkl_dcosy",&trkl_dcosy,"trkl_dcosy/F");

  neus->Branch("shwl_number",&shwl_number,"shwl_number/I");
  neus->Branch("shwl_energy",&shwl_energy,"shwl_energy/F");
  neus->Branch("shwl_vtxx",&shwl_vtxx,"shwl_vtxx/F");
  neus->Branch("shwl_vtxy",&shwl_vtxy,"shwl_vtxy/F");
  neus->Branch("shwl_vtxz",&shwl_vtxz,"shwl_vtxz/F");
  neus->Branch("shwl_vtxu",&shwl_vtxu,"shwl_vtxu/F");
  neus->Branch("shwl_vtxv",&shwl_vtxv,"shwl_vtxv/F");

  neus->Branch("shwl_planes",&shwl_planes,"shwl_planes/I");
  neus->Branch("shwl_uplanes",&shwl_uplanes,"shwl_uplanes/I");
  neus->Branch("shwl_vplanes",&shwl_vplanes,"shwl_vplanes/I");
  neus->Branch("shwl_begplane",&shwl_begplane,"shwl_begplane/I");
  neus->Branch("shwl_ubegplane",&shwl_ubegplane,"shwl_ubegplane/I");
  neus->Branch("shwl_vbegplane",&shwl_vbegplane,"shwl_vbegplane/I");
  neus->Branch("shwl_endplane",&shwl_endplane,"shwl_endplane/I");
  neus->Branch("shwl_uendplane",&shwl_uendplane,"shwl_uendplane/I");
  neus->Branch("shwl_vendplane",&shwl_vendplane,"shwl_vendplane/I");

  neus->Branch("shwl_dcosz",&shwl_dcosz,"shwl_dcosz/F");
  neus->Branch("shwl_dcosx",&shwl_dcosx,"shwl_dcosx/F");
  neus->Branch("shwl_dcosy",&shwl_dcosy,"shwl_dcosy/F");

  neus->Branch("shwl_count",&shwl_count,"shwl_count/I");

  neus->Branch("evt_number",&evt_number,"evt_number/I");
  neus->Branch("evt_ntracks",&evt_ntracks,"evt_ntracks/I");
  neus->Branch("evt_nshowers",&evt_nshowers,"evt_nshowers/I");
  neus->Branch("evt_vtxx",&evt_vtxx,"evt_vtxx/F");
  neus->Branch("evt_vtxy",&evt_vtxy,"evt_vtxy/F");
  neus->Branch("evt_vtxz",&evt_vtxz,"evt_vtxz/F");
  neus->Branch("evt_vtxu",&evt_vtxu,"evt_vtxu/F");
  neus->Branch("evt_vtxv",&evt_vtxv,"evt_vtxv/F");
  neus->Branch("evt_endx",&evt_endx,"evt_endx/F");
  neus->Branch("evt_endy",&evt_endy,"evt_endy/F");
  neus->Branch("evt_endz",&evt_endz,"evt_endz/F");
  neus->Branch("evt_endu",&evt_endu,"evt_endu/F");
  neus->Branch("evt_endv",&evt_endv,"evt_endv/F");
  
  neus->Branch("evt_neu_energy",&evt_neu_energy,"evt_neu_energy/F");
  neus->Branch("evt_q2",&evt_q2,"evt_q2/F");
  neus->Branch("evt_y",&evt_y,"evt_y/F");
  neus->Branch("evt_x",&evt_x,"evt_x/F");
  neus->Branch("evt_pid",&evt_pid,"evt_pid/F");

  neus->Branch("mcl_index",&mcl_index,"mcl_index/I");

  neus->Branch("mcl_mu_energy",&mcl_mu_energy,"mcl_mu_energy/F");
  neus->Branch("mcl_mu_px",&mcl_mu_px,"mcl_mu_px/F");
  neus->Branch("mcl_mu_py",&mcl_mu_py,"mcl_mu_py/F");
  neus->Branch("mcl_mu_pz",&mcl_mu_pz,"mcl_mu_pz/F");

  neus->Branch("mcl_shw_energy",&mcl_shw_energy,"mcl_shw_energy/F");
  neus->Branch("mcl_shw_px",&mcl_shw_px,"mcl_shw_px/F");
  neus->Branch("mcl_shw_py",&mcl_shw_py,"mcl_shw_py/F");
  neus->Branch("mcl_shw_pz",&mcl_shw_pz,"mcl_shw_pz/F");

  neus->Branch("mcl_nu_energy",&mcl_nu_energy,"mcl_nu_energy/F");
  neus->Branch("mcl_nu_px",&mcl_nu_px,"mcl_nu_px/F");
  neus->Branch("mcl_nu_py",&mcl_nu_py,"mcl_nu_py/F");
  neus->Branch("mcl_nu_pz",&mcl_nu_pz,"mcl_nu_pz/F");

  neus->Branch("mcl_tar_energy",&mcl_tar_energy,"mcl_tar_energy/F");
  neus->Branch("mcl_tar_px",&mcl_tar_px,"mcl_tar_px/F");
  neus->Branch("mcl_tar_py",&mcl_tar_py,"mcl_tar_py/F");
  neus->Branch("mcl_tar_pz",&mcl_tar_pz,"mcl_tar_pz/F");

  neus->Branch("mcl_a",&mcl_a,"mcl_a/F");
  neus->Branch("mcl_inu",&mcl_inu,"mcl_inu/I");
  neus->Branch("mcl_iaction",&mcl_iaction,"mcl_iaction/I");
  neus->Branch("mcl_ires",&mcl_ires,"mcl_res/I");
  neus->Branch("mcl_x",&mcl_x,"mcl_x/F");
  neus->Branch("mcl_y",&mcl_y,"mcl_y/F");
  neus->Branch("mcl_q2",&mcl_q2,"mcl_q2/F");
  neus->Branch("mcl_w2",&mcl_w2,"mcl_w2/F");
  neus->Branch("mcl_initstate",&mcl_initstate,"mcl_initstate/I");
  neus->Branch("mcl_nucleus",&mcl_nucleus,"mcl_nucleus/I");
  neus->Branch("mcl_hadfs",&mcl_hadfs,"mcl_hadfs/I");



  //DECLARING THE ANTI-NEUTRINO TREE  

  aneus->Branch("header_run",&header_run,"header_run/I");
  aneus->Branch("header_subrun",&header_subrun,"header_subrun/I");
  aneus->Branch("header_snarl",&header_snarl,"header_snarl/I");

  aneus->Branch("trkl_vtxx",&trkl_vtxx,"trkl_vtxx/F");
  aneus->Branch("trkl_vtxy",&trkl_vtxy,"trkl_vtxy/F");
  aneus->Branch("trkl_vtxz",&trkl_vtxz,"trkl_vtxz/F");
  aneus->Branch("trkl_vtxu",&trkl_vtxu,"trkl_vtxu/F");
  aneus->Branch("trkl_vtxv",&trkl_vtxv,"trkl_vtxv/F");
  aneus->Branch("trkl_vtxrad",&trkl_vtxrad,"trkl_vtxrad/F");
  aneus->Branch("trkl_endx",&trkl_endx,"trkl_endx/F");
  aneus->Branch("trkl_endy",&trkl_endy,"trkl_endy/F");
  aneus->Branch("trkl_endz",&trkl_endz,"trkl_endz/F");
  aneus->Branch("trkl_endu",&trkl_endu,"trkl_endu/F");
  aneus->Branch("trkl_endv",&trkl_endv,"trkl_endv/F");
  aneus->Branch("trkl_endrad",&trkl_endrad,"trkl_endrad/F");

  aneus->Branch("trkl_number",&trkl_number,"trkl_number/I");
  aneus->Branch("trkl_planes",&trkl_planes,"trkl_planes/I");
  aneus->Branch("trkl_uplanes",&trkl_uplanes,"trkl_uplanes/I");
  aneus->Branch("trkl_vplanes",&trkl_vplanes,"trkl_vplanes/I");
  aneus->Branch("trkl_begplane",&trkl_begplane,"trkl_begplane/I");
  aneus->Branch("trkl_ubegplane",&trkl_ubegplane,"trkl_ubegplane/I");
  aneus->Branch("trkl_vbegplane",&trkl_vbegplane,"trkl_vbegplane/I");
  aneus->Branch("trkl_endplane",&trkl_endplane,"trkl_endplane/I");
  aneus->Branch("trkl_uendplane",&trkl_uendplane,"trkl_uendplane/I");
  aneus->Branch("trkl_vendplane",&trkl_vendplane,"trkl_vendplane/I");

  aneus->Branch("trkl_redchi2",&trkl_redchi2,"trkl_redchi2/F");
  aneus->Branch("trkl_pass",&trkl_pass,"trkl_pass/I");
  aneus->Branch("trkl_nstps",&trkl_nstps,"trkl_nstps/I");
  aneus->Branch("trkl_coilhole",&trkl_coilhole,"trkl_coilhole/I");
  aneus->Branch("trkl_exitstop",&trkl_exitstop,"trkl_exitstop/I");
  aneus->Branch("trkl_new_pflag",&trkl_new_pflag,"trkl_new_pflag/I");
  aneus->Branch("trkl_trklike",&trkl_trklike,"trkl_trklike/I");
  aneus->Branch("trkl_oldpos_stps",&trkl_oldpos_stps,"trkl_oldpos_stps/I");
  aneus->Branch("trkl_oldneg_stps",&trkl_oldneg_stps,"trkl_oldneg_stps/I");
  aneus->Branch("trkl_pos_stps",&trkl_pos_stps,"trkl_pos_stps/I");
  aneus->Branch("trkl_neg_stps",&trkl_neg_stps,"trkl_neg_stps/I");

  aneus->Branch("trkl_fitmom",&trkl_fitmom,"trkl_fitmom/F");
  aneus->Branch("trkl_qoverp",&trkl_qoverp,"trkl_qoverp/F");
  aneus->Branch("trkl_eqp",&trkl_eqp,"trkl_eqp/F");
  aneus->Branch("trkl_rangemom",&trkl_rangemom,"trkl_rangemom/F");
  aneus->Branch("trkl_energy",&trkl_energy,"trkl_energy/F");
  aneus->Branch("trkl_length",&trkl_length,"trkl_length/F");
  aneus->Branch("trkl_prob",&trkl_prob,"trkl_prob/F");

  aneus->Branch("trkl_dcosz",&trkl_dcosz,"trkl_dcosz/F");
  aneus->Branch("trkl_dcosx",&trkl_dcosx,"trkl_dcosx/F");
  aneus->Branch("trkl_dcosy",&trkl_dcosy,"trkl_dcosy/F");

  aneus->Branch("shwl_number",&shwl_number,"shwl_number/I");
  aneus->Branch("shwl_energy",&shwl_energy,"shwl_energy/F");
  aneus->Branch("shwl_vtxx",&shwl_vtxx,"shwl_vtxx/F");
  aneus->Branch("shwl_vtxy",&shwl_vtxy,"shwl_vtxy/F");
  aneus->Branch("shwl_vtxz",&shwl_vtxz,"shwl_vtxz/F");
  aneus->Branch("shwl_vtxu",&shwl_vtxu,"shwl_vtxu/F");
  aneus->Branch("shwl_vtxv",&shwl_vtxv,"shwl_vtxv/F");

  aneus->Branch("shwl_planes",&shwl_planes,"shwl_planes/I");
  aneus->Branch("shwl_uplanes",&shwl_uplanes,"shwl_uplanes/I");
  aneus->Branch("shwl_vplanes",&shwl_vplanes,"shwl_vplanes/I");
  aneus->Branch("shwl_begplane",&shwl_begplane,"shwl_begplane/I");
  aneus->Branch("shwl_ubegplane",&shwl_ubegplane,"shwl_ubegplane/I");
  aneus->Branch("shwl_vbegplane",&shwl_vbegplane,"shwl_vbegplane/I");
  aneus->Branch("shwl_endplane",&shwl_endplane,"shwl_endplane/I");
  aneus->Branch("shwl_uendplane",&shwl_uendplane,"shwl_uendplane/I");
  aneus->Branch("shwl_vendplane",&shwl_vendplane,"shwl_vendplane/I");

  aneus->Branch("shwl_dcosz",&shwl_dcosz,"shwl_dcosz/F");
  aneus->Branch("shwl_dcosx",&shwl_dcosx,"shwl_dcosx/F");
  aneus->Branch("shwl_dcosy",&shwl_dcosy,"shwl_dcosy/F");

  aneus->Branch("shwl_count",&shwl_count,"shwl_count/I");

  aneus->Branch("evt_number",&evt_number,"evt_number/I");
  aneus->Branch("evt_ntracks",&evt_ntracks,"evt_ntracks/I");
  aneus->Branch("evt_nshowers",&evt_nshowers,"evt_nshowers/I");
  aneus->Branch("evt_vtxx",&evt_vtxx,"evt_vtxx/F");
  aneus->Branch("evt_vtxy",&evt_vtxy,"evt_vtxy/F");
  aneus->Branch("evt_vtxz",&evt_vtxz,"evt_vtxz/F");
  aneus->Branch("evt_vtxu",&evt_vtxu,"evt_vtxu/F");
  aneus->Branch("evt_vtxv",&evt_vtxv,"evt_vtxv/F");
  aneus->Branch("evt_endx",&evt_endx,"evt_endx/F");
  aneus->Branch("evt_endy",&evt_endy,"evt_endy/F");
  aneus->Branch("evt_endz",&evt_endz,"evt_endz/F");
  aneus->Branch("evt_endu",&evt_endu,"evt_endu/F");
  aneus->Branch("evt_endv",&evt_endv,"evt_endv/F");

  aneus->Branch("evt_neu_energy",&evt_neu_energy,"evt_neu_energy/F");
  aneus->Branch("evt_q2",&evt_q2,"evt_q2/F");
  aneus->Branch("evt_y",&evt_y,"evt_y/F");
  aneus->Branch("evt_x",&evt_x,"evt_x/F");
  aneus->Branch("evt_pid",&evt_pid,"evt_pid/F");

  aneus->Branch("mcl_index",&mcl_index,"mcl_index/I");

  aneus->Branch("mcl_mu_energy",&mcl_mu_energy,"mcl_mu_energy/F");
  aneus->Branch("mcl_mu_px",&mcl_mu_px,"mcl_mu_px/F");
  aneus->Branch("mcl_mu_py",&mcl_mu_py,"mcl_mu_py/F");
  aneus->Branch("mcl_mu_pz",&mcl_mu_pz,"mcl_mu_pz/F");

  aneus->Branch("mcl_shw_energy",&mcl_shw_energy,"mcl_shw_energy/F");
  aneus->Branch("mcl_shw_px",&mcl_shw_px,"mcl_shw_px/F");
  aneus->Branch("mcl_shw_py",&mcl_shw_py,"mcl_shw_py/F");
  aneus->Branch("mcl_shw_pz",&mcl_shw_pz,"mcl_shw_pz/F");

  aneus->Branch("mcl_nu_energy",&mcl_nu_energy,"mcl_nu_energy/F");
  aneus->Branch("mcl_nu_px",&mcl_nu_px,"mcl_nu_px/F");
  aneus->Branch("mcl_nu_py",&mcl_nu_py,"mcl_nu_py/F");
  aneus->Branch("mcl_nu_pz",&mcl_nu_pz,"mcl_nu_pz/F");

  aneus->Branch("mcl_tar_energy",&mcl_tar_energy,"mcl_tar_energy/F");
  aneus->Branch("mcl_tar_px",&mcl_tar_px,"mcl_tar_px/F");
  aneus->Branch("mcl_tar_py",&mcl_tar_py,"mcl_tar_py/F");
  aneus->Branch("mcl_tar_pz",&mcl_tar_pz,"mcl_tar_pz/F");

  aneus->Branch("mcl_a",&mcl_a,"mcl_a/F");
  aneus->Branch("mcl_inu",&mcl_inu,"mcl_inu/I");
  aneus->Branch("mcl_iaction",&mcl_iaction,"mcl_iaction/I");
  aneus->Branch("mcl_ires",&mcl_ires,"mcl_res/I");
  aneus->Branch("mcl_x",&mcl_x,"mcl_x/F");
  aneus->Branch("mcl_y",&mcl_y,"mcl_y/F");
  aneus->Branch("mcl_q2",&mcl_q2,"mcl_q2/F");
  aneus->Branch("mcl_w2",&mcl_w2,"mcl_w2/F");
  aneus->Branch("mcl_initstate",&mcl_initstate,"mcl_initstate/I");
  aneus->Branch("mcl_nucleus",&mcl_nucleus,"mcl_nucleus/I");
  aneus->Branch("mcl_hadfs",&mcl_hadfs,"mcl_hadfs/I");


  //ndfile<<"\n\n beginning a new run\n"<<endl;                  

  for (Int_t jentry=0;jentry<nentries;jentry++)         
  {
    //ndfile<<"entry no "<<jentry<<endl;
        
    Long64_t ientry = LoadTree(jentry);                       
    if (ientry < 0) break;

    nb = fChain->GetEntry(jentry);                              //loads all the branches

    if(evsum_ntracks==0) continue;                              //skip this event (only shw)
   		
    else if(evsum_ntracks<=2){ 

      //INITIALIZE THE VARIABLES

       trk_ind = 0;
       trk_ind_next = 1;


        header_run=0;
        header_subrun=0;
        header_snarl=0;

	trkl_number=0;
	trkl_vtxx=0;
	trkl_vtxy=0;
	trkl_vtxz=0;
	trkl_vtxu=0;
	trkl_vtxv=0;
	trkl_vtxrad=0;
	trkl_endx=0;
	trkl_endy=0;
	trkl_endz=0;
	trkl_endu=0;
	trkl_endv=0;
	trkl_endrad=0;

	trkl_planes=0;
	trkl_uplanes=0;
	trkl_vplanes=0;
	trkl_begplane=0;
	trkl_ubegplane=0;
	trkl_vbegplane=0;
	trkl_endplane=0; 
	trkl_uendplane=0; 
	trkl_vendplane=0; 
  
	trkl_redchi2=0;
	trkl_pass=0;
	trkl_nstps=0;
	trkl_coilhole=0;
	trkl_exitstop=0;
	trkl_new_pflag=0;
	trkl_oldpos_stps=0;
	trkl_oldneg_stps=0;
	trkl_pos_stps=0;
	trkl_neg_stps=0;
	trkl_trklike=0;

	trkl_fitmom=0;
	trkl_qoverp=0;
	trkl_eqp=0;
	trkl_rangemom=0;
        trkl_energy=0;
	trkl_length=0;
	trkl_prob=0;

	trkl_dcosz=0;
	trkl_dcosx=0;
	trkl_dcosy=0;

	shw_ind=-1;
	shwl_number=-1;
	shwl_energy=0;
	shwl_vtxx=1000;
	shwl_vtxy=1000;
	shwl_vtxz=1000;
	shwl_vtxu=1000;
	shwl_vtxv=1000;

	shwl_planes=1000;
	shwl_uplanes=1000;
	shwl_vplanes=1000;
	shwl_begplane=1000;
	shwl_ubegplane=1000;
	shwl_vbegplane=1000;
	shwl_endplane=1000;
	shwl_uendplane=1000;
	shwl_vendplane=1000;

	shwl_dcosz=1000;
	shwl_dcosx=1000;
	shwl_dcosy=1000;
	shwl_count=0;
  
	evt_number=0;
	evt_ntracks=0;
	evt_nshowers=0;

	evt_vtxx=0;
	evt_vtxy=0;
	evt_vtxz=0;
	evt_vtxu=0;
	evt_vtxv=0;
	evt_endx=0;
	evt_endy=0;
	evt_endz=0;
	evt_endu=0;
	evt_endv=0;

	evt_neu_energy=0;
	evt_y=0;
	evt_x=0;
	evt_q2=0;
	evt_pid=0;

        mcl_index=0;

	mcl_mu_energy=0;
	mcl_mu_px=0;
	mcl_mu_py=0;
	mcl_mu_pz=0;

	mcl_shw_energy=0;
	mcl_shw_px=0;
	mcl_shw_py=0;
	mcl_shw_pz=0;
	
	mcl_nu_energy=0;
	mcl_nu_px=0;
	mcl_nu_py=0;
	mcl_nu_pz=0;
	
	mcl_tar_energy=0;
	mcl_tar_px=0;
	mcl_tar_py=0;
	mcl_tar_pz=0;
	
	mcl_a=0;
	mcl_inu=0;
	mcl_iaction=0;
	mcl_ires=0;
	mcl_x=0;
	mcl_y=0;
	mcl_q2=0;
	mcl_w2=0;
	mcl_initstate=0;
	mcl_nucleus=0;
	mcl_hadfs=0;



	if(evsum_ntracks==1){
          trk_ind_big = trk_ind;

	}
        
      //selecting the longer track

        if(evsum_ntracks==2)
          {
              if(trkl_ds[trk_ind]>trkl_ds[trk_ind_next])
                {
		  trk_ind_big=trk_ind;
		  trk_ind_little=trk_ind_next;
                }else{
		  trk_ind_big=trk_ind_next;
                  trk_ind_little=trk_ind;
                }

	  }
          

	//FILLING IN THE HEADER INFORMATION

	header_run    = hdr_run;
	header_subrun = hdr_subrun;
	header_snarl  = hdr_snarl;        

	ndfile<<"\nrun and snarl "<<header_run<<"  "<<header_snarl<<endl;

	//FILLING IN THE TRACK INFORMATION

       if(trkl_qp[trk_ind_big]!=0)
       trkl_fitmom   = (1/trkl_qp[trk_ind_big]);
       else
       trkl_fitmom   = 0;

       
       trkl_qoverp   = trkl_qp[trk_ind_big];
       trkl_eqp      = trkl_dqp[trk_ind_big];
       trkl_rangemom = trkl_prange[trk_ind_big];
       trkl_length   = trkl_ds[trk_ind_big];

       trkl_prob     = TMath::Prob(trkl_chi2[trk_ind_big],trkl_ndof[trk_ind_big]);


       trkl_vtxx     = trkl_beg_x[trk_ind_big];
       trkl_vtxy     = trkl_beg_y[trk_ind_big];
       trkl_vtxz     = trkl_beg_z[trk_ind_big];
       trkl_vtxu     = trkl_beg_u[trk_ind_big];
       trkl_vtxv     = trkl_beg_v[trk_ind_big];
       trkl_vtxrad    = TMath::Sqrt((trkl_vtxx*trkl_vtxx)+(trkl_vtxy*trkl_vtxy));

       trkl_endx     = trkl_end_x[trk_ind_big];
       trkl_endy     = trkl_end_y[trk_ind_big];
       trkl_endz     = trkl_end_z[trk_ind_big];
       trkl_endu     = trkl_end_u[trk_ind_big];
       trkl_endv     = trkl_end_v[trk_ind_big];
       trkl_endrad    = TMath::Sqrt((trkl_endx*trkl_endx)+(trkl_endy*trkl_endy));

       trkl_number      = trkl_index[trk_ind_big];

       trkl_planes   = trkl_nplanes[trk_ind_big];
       trkl_uplanes  = trkl_nuplanes[trk_ind_big];
       trkl_vplanes  = trkl_nvplanes[trk_ind_big];

       trkl_begplane = trkl_beg_plane[trk_ind_big];
       trkl_ubegplane = trkl_begpu[trk_ind_big];
       trkl_vbegplane = trkl_begpv[trk_ind_big];

       trkl_endplane = trkl_end_plane[trk_ind_big];
       trkl_uendplane = trkl_endpu[trk_ind_big];
       trkl_vendplane = trkl_endpv[trk_ind_big];

       trkl_redchi2   = trkl_rchi2[trk_ind_big];
       trkl_pass      = trkl_fitpass[trk_ind_big];
       trkl_nstps     = trkl_nstrips[trk_ind_big];      
       trkl_exitstop  = trkl_pflag[trk_ind_big];
       trkl_trklike   = trkl_ntrklike[trk_ind_big];

       trkl_dcosz     = trkl_beg_dcosz[trk_ind_big];
       trkl_dcosx     = trkl_beg_dcosx[trk_ind_big];
       trkl_dcosy     = trkl_beg_dcosy[trk_ind_big];


       //ndfile<<"run "<<header_run<<" ,snarl "<<header_snarl<<" ,#of trks "<<evsum_ntracks<<" ,trk index "<<trkl_number<<" ,#of shws "<<evsum_nshowers<<endl;

       //FILLING IN THE SHOWER INFORMATION


      //adding the shower energy

      if(evsum_nshowers==0){
	shwl_energy=0;
	shwl_number=-1;
	shw_ind=-1;
	shwl_count=0;
	//ndfile<<"no shower at all"<<endl;

      }else {
	

      for(int count =0;count<evsum_nshowers;count++)
        {

	 //ndfile<<"trk "<<trkl_vtxz<<" ,shw "<<shwl_vtx_z[count]<<" ,diff "<<fabs(trkl_vtxz-shwl_vtx_z[count])<<endl;
	  
	 if(fabs(trkl_vtxz-shwl_vtx_z[count])<0.5){

         shwl_energy =  shwl_energy + shwl_shwphlinCC[count];
	 shwl_count++;

	 if(shwl_count==1)
	   {
	     shwl_number=shwl_index[count];
	     shw_ind = count;
	   }


	 /*if(count==0)
	 //ndfile<<"found a matching shw with first"<<endl;
	 if(count==1)
	 //ndfile<<"found a matching shw with second"<<endl;
	 if(count==2)
	 //ndfile<<"found a matching shw with third"<<endl;
	 if(count==3)
	 //ndfile<<"found a matching shw with fourth"<<endl;
	 if(count==4)
	 //ndfile<<"found a matching shw with fifth"<<endl;
	 if(count==5)
	 //ndfile<<"found a matching shw with sixth"<<endl;
	 if(count==6)
	 //ndfile<<"found a matching shw with seventh"<<endl;*/

	 }
	  


	}

	 if(shwl_count==0){
	  
            shwl_energy=0;
	    shwl_number=-1;
	    shw_ind=0;
	    shwl_count=0;
	    //ndfile<<"did not find a match"<<endl;
	 
      }
      }

      //ndfile<<"FINAL RESULT :#of shws "<<shwl_count<<",index "<<shw_ind<<",original "<<shwl_number<<endl;


									       
      //if(evsum_ntracks==2)
      //shwl_energy = shwl_energy + trkl_prange[trk_ind_little];

      if(evsum_nshowers!=0){

       shwl_vtxx     = shwl_vtx_x[shw_ind];
       shwl_vtxy     = shwl_vtx_y[shw_ind];
       shwl_vtxz     = shwl_vtx_z[shw_ind];
       shwl_vtxu     = shwl_vtx_u[shw_ind];
       shwl_vtxv     = shwl_vtx_v[shw_ind];

       //ndfile<<"assigining values "<<shwl_vtx_z[shw_ind]<<endl;

       shwl_planes   = shwl_nplanes[shw_ind];
       shwl_uplanes  = shwl_nuplanes[shw_ind];
       shwl_vplanes  = shwl_nvplanes[shw_ind];

       shwl_begplane = shwl_beg[shw_ind];
       shwl_ubegplane = shwl_begu[shw_ind];
       shwl_vbegplane = shwl_begv[shw_ind];

       shwl_endplane = shwl_end[shw_ind];
       shwl_uendplane = shwl_endu[shw_ind];
       shwl_vendplane = shwl_endv[shw_ind];       

       shwl_dcosz = shwl_vtx_dcosz[shw_ind];
       shwl_dcosx = shwl_vtx_dcosx[shw_ind];
       shwl_dcosy = shwl_vtx_dcosy[shw_ind];

      }
       
       //FILLING IN THE EVENT INFORMATION

       evt_number    = evsum_event;
       evt_ntracks   = evsum_ntracks;
       evt_nshowers  = evsum_nshowers;
       evt_pid       = evsum_pid;

       evt_vtxx      = svtx_x;
       evt_vtxy      = svtx_y;
       evt_vtxz      = svtx_z;
       evt_vtxu      = svtx_u;
       evt_vtxv      = svtx_v;

       evt_endx      = evtx_x;
       evt_endy      = evtx_y;
       evt_endz      = evtx_z;
       evt_endu      = evtx_u;
       evt_endv      = evtx_v;


       //determining if the track passes through the coil hole
                  
       for(int mm=0;mm<trkl_nstps;mm++){

         Float_t trkl_stprad = 
TMath::Sqrt((trkl_stpx[trk_ind_big][mm]*trkl_stpx[trk_ind_big][mm])+(trkl_stpy[trk_ind_big][mm]*trkl_stpy[trk_ind_big][mm]));

         if(trkl_stprad<0.5){

           trkl_coilhole=trkl_coilhole+1;

         }


       }

   	//determining the number of positive and negative q/p's
	//instead of looking at all the strips, I will look at 3/4th of the strips

        for(int stp=0;stp<trkl_nstps;stp++)
	  {

	    if(trkl_stpfitqp[trk_ind_big][stp] > 0) trkl_oldpos_stps++;
	   
	    if(trkl_stpfitqp[trk_ind_big][stp] < 0) trkl_oldneg_stps++;

	  }



/////////////////////////////////////////////////////////////////////////////////////////////////////////

  /* Float_t        *trkl_stpu[kMaxtrkl];   //[trkl_nstrips]
   Float_t        *trkl_stpv[kMaxtrkl];   //[trkl_nstrips]
   Float_t        *trkl_stpx[kMaxtrkl];   //[trkl_nstrips]
   Float_t        *trkl_stpy[kMaxtrkl];   //[trkl_nstrips]
   Float_t        *trkl_stpz[kMaxtrkl];   //[trkl_nstrips]*/


//	  ndfile<<trkl_begplane<<"  "<<trkl_endplane<<endl;

	  Float_t xCurv = XCurvature(trkl_begplane,trkl_endplane,fitLengths,trkl_nstps,trkl_stpu[trk_ind_big],trkl_stpv[trk_ind_big],trkl_stpx[trk_ind_big],trkl_stpy[trk_ind_big],trkl_stpz[trk_ind_big]);
	  Float_t yCurv = YCurvature(trkl_begplane,trkl_endplane,fitLengths,trkl_nstps,trkl_stpu[trk_ind_big],trkl_stpv[trk_ind_big],trkl_stpx[trk_ind_big],trkl_stpy[trk_ind_big],trkl_stpz[trk_ind_big]);

          ndfile<<"final xCurv "<<xCurv<<" yCurv "<<yCurv<<endl;

/////////////////////////////////////////////////////////////////////////////////////////////////////////


        //taking the events thrown out by the us stopping cut and putting them in the us exiting sample.

      if(trkl_exitstop==0 && trkl_endz<7 && trkl_endrad>0.5 && trkl_endrad<0.8)
        {
          trkl_exitstop=2;
        }


    //revising the stopping and exiting cuts


      //UNCHANGED FOR THE UPSTREAM TRACKS

      if(trkl_exitstop==1 || trkl_exitstop==2){

	trkl_new_pflag=trkl_exitstop;
      }


      //CHANGED FOR THE DOWNSTREAM TRACKS

      if(trkl_exitstop==3 || trkl_exitstop==4){


	if(trkl_endu>0.3 && trkl_endu<1.8 && trkl_endv>-1.8 && trkl_endv<-0.3 && trkl_endx<2.4 && trkl_endrad>0.8 && trkl_endz<15){

          trkl_new_pflag=3;
	  //ndfile<<"old "<<trkl_exitstop<<" new "<<trkl_new_pflag<<endl;
	}
	

        if(trkl_endrad>0.5 && ((trkl_endu<0.3 || trkl_endu>1.8 || trkl_endv<-1.8 || trkl_endv>-0.3 || trkl_endx>2.4) || trkl_endz>15)){

          trkl_new_pflag=4;
	  //ndfile<<"old "<<trkl_exitstop<<" new "<<trkl_new_pflag<<endl;

	}


	if(trkl_endrad>0.5 && trkl_endrad<0.8)
	  trkl_new_pflag=4;
	  //ndfile<<"old "<<trkl_exitstop<<" new "<<trkl_new_pflag<<endl;
      }


       //determining the kinematic variables for stopping and exiting

       if(trkl_new_pflag==1 || trkl_new_pflag==3)
	 {
          
           trkl_energy   = TMath::Sqrt((trkl_rangemom*trkl_rangemom)+(0.106*0.106));
           evt_neu_energy = trkl_energy + shwl_energy;

           //ndfile<<header_run<<" "<<header_snarl<<" "<<trkl_number<<" "<<trkl_energy<<" "<<evt_neu_energy<<endl;

	   evt_y = shwl_energy/evt_neu_energy;
	   evt_q2 = 2*(evt_neu_energy)*(trkl_energy)*(1-trkl_dcosz);

           if(shwl_energy!=0)
	     evt_x =  evt_q2/(2*0.938*shwl_energy);

	   }


         else if(trkl_new_pflag==2 || trkl_new_pflag==4)
	 {
          
           trkl_energy   = TMath::Sqrt((trkl_fitmom*trkl_fitmom)+(0.106*0.106));
           evt_neu_energy = trkl_energy + shwl_energy;
           

           if(evt_neu_energy!=0)           
	   evt_y = shwl_energy/evt_neu_energy;

	   evt_q2 = 2*(evt_neu_energy)*(trkl_energy)*(1-trkl_dcosz);

           if(shwl_energy!=0)
	   evt_x =  evt_q2/(2*0.938*shwl_energy);

	   }

         else if(trkl_exitstop==0 && trkl_endz>15)
	   {
	     trkl_new_pflag=0;

	     trkl_energy   = TMath::Sqrt((trkl_fitmom*trkl_fitmom)+(0.106*0.106));
	     evt_neu_energy=  trkl_energy + shwl_energy;
	     evt_neu_energy = trkl_energy + shwl_energy;
           
	     //ndfile<<header_run<<" "<<header_snarl<<" "<<trkl_number<<" "<<trkl_energy<<" "<<evt_neu_energy<<endl;

	     if(evt_neu_energy!=0)           
	       evt_y = shwl_energy/evt_neu_energy;

	     evt_q2 = 2*(evt_neu_energy)*(trkl_energy)*(1-trkl_dcosz);

	     if(shwl_energy!=0)
	       evt_x =  evt_q2/(2*0.938*shwl_energy);
             
	   }

         else if(trkl_exitstop==0 && trkl_endz<15)
	   {
	     trkl_new_pflag=0;

	     trkl_energy   = TMath::Sqrt((trkl_rangemom*trkl_rangemom)+(0.106*0.106));
	     evt_neu_energy= trkl_energy + shwl_energy;
	     evt_neu_energy = trkl_energy + shwl_energy;

	     //ndfile<<header_run<<" "<<header_snarl<<" "<<trkl_number<<" "<<trkl_energy<<" "<<evt_neu_energy<<endl;
	     
	     evt_y = shwl_energy/evt_neu_energy;
	     evt_q2 = 2*(evt_neu_energy)*(trkl_energy)*(1-trkl_dcosz);

	     if(shwl_energy!=0)
	       evt_x =  evt_q2/(2*0.938*shwl_energy);

	   }







       //FILLING IN THE MC INFORMATION

        mcl_index=mc_index;

	mcl_mu_energy=mc_mu1_E;
	mcl_mu_px=mc_mu1_Px;
	mcl_mu_py=mc_mu1_Py;
	mcl_mu_pz=mc_mu1_Pz;

	mcl_shw_energy=mc_shw_E;
	mcl_shw_px=mc_shw_Px;
	mcl_shw_py=mc_shw_Py;
	mcl_shw_pz=mc_shw_Pz;
	
	mcl_nu_energy=mc_nu_E;
	mcl_nu_px=mc_nu_Px;
	mcl_nu_py=mc_nu_Py;
	mcl_nu_pz=mc_nu_Pz;
	
	mcl_tar_energy=mc_tgt_E;
	mcl_tar_px=mc_tgt_Px;
	mcl_tar_py=mc_tgt_Py;
	mcl_tar_pz=mc_tgt_Pz;
	
	mcl_a=mc_A;
	mcl_inu=mc_nu_pdg;
	mcl_iaction=mc_curr;
	mcl_ires=mc_reac;
	mcl_x=mc_x;
	mcl_y=mc_y;
	mcl_q2=mc_Q2;
	mcl_w2=mc_W*mc_W;

	mcl_nucleus=0;
	mcl_hadfs=0;
        mc_mul=0;

	//ndfile<<"mc info "<<mcl_index<<"  "<<mcl_ires<<endl;

                  //calculating initial state, hadronic final state and nucleus


                  Int_t protneut = -1;  // 0 = neutron, 1 = proton, 2 = N, 3 = A
                  Int_t nubarnu = 0;    // +1 = neutrino, -1 = antineutrino

                      if(abs(mc_nu_pdg)==12 ||abs(mc_nu_pdg)==14 ||abs(mc_nu_pdg)==16){
                          nubarnu = mc_nu_pdg/abs(mc_nu_pdg);  //get sign
		      }


                        if(mc_tgt_pdg==2212) protneut = 1; //proton
                        else if(mc_tgt_pdg==2112) protneut = 0; //neutron
                        else if(abs(mc_tgt_pdg)>1000000000) protneut = 2; //nucleus
                        else protneut = 3; //atom - probably never get here since IdHEP A>N?
                      

		  if(protneut==1 && nubarnu==1)  mcl_initstate=1;  //p + v
                  if(protneut==0 && nubarnu==1)  mcl_initstate=2;  //n + v
                  if(protneut==1 && nubarnu==-1) mcl_initstate=3;  //p + vbar
                  if(protneut==0 && nubarnu==-1) mcl_initstate=4;  //n + vbar
                  if(protneut==2 && nubarnu==1)  mcl_initstate=5;  //N + v
                  if(protneut==3 && nubarnu==1)  mcl_initstate=6;  //A + v
                  if(protneut==2 && nubarnu==-1) mcl_initstate=7;  //N + vbar
                  if(protneut==3 && nubarnu==-1) mcl_initstate=8;  //A + vbar


		  //end( calculating initial state)



                  //calculating the nucleus
                  Int_t z=int(mc_Z);
                  mcl_nucleus=1;

                  switch (z) {
                  case 1:
                    mcl_nucleus=0;   // free nucleon
                    break;
                  case 6:
                    mcl_nucleus=274; // carbon
                    break;
                  case 8:
                    mcl_nucleus=284; // oxygen
		    break;
                  case 26:
                    mcl_nucleus=372; // iron
                    break;
                  default:
                    mcl_nucleus=1;   // unknown
                    break;
		  }


                  //calculating the hadronic final state
		  // stdhep_ is the total no of stdhep entries for a given mc info
			//ndfile<<"stdhep entries "<<stdhep_<<endl;

		  	for(int count=0;count<stdhep_;count++){
		        if(stdhep_ist[count]==3){				
			mcl_hadfs = stdhep_pdg[count];
			}
			}

                  /*if(mcl_ires==1002){
                        mcl_hadfs = mc_res_pdg;
			mcl_hadfs = mcl_hadfs%1000;
		  }
		  
                 
                  else {
                        mc_mul =  mc_np+mc_nn+mc_npi0+mc_npip+mc_npim+mc_nK0+mc_nKp+mc_nKm;
			mcl_hadfs = (4*(mc_mul+20)) + 128;
		  }*/

			mcl_hadfs = mcl_hadfs%1000;
		                 
   // reconstruction cuts
	
   if(trkl_vtxz<0.6 || trkl_vtxz>5) continue;
   if(trkl_vtxu<0.3 || trkl_vtxu>1.8) continue;
   if(trkl_vtxv<-1.8 || trkl_vtxv>-0.3) continue;
   if(trkl_vtxx>2.4) continue;
   if(trkl_vtxrad<0.8) continue;

   if(evt_vtxz>5) continue;


   //ndfile<<header_snarl<<"  "<<mcl_index<<"   "<<mcl_ires<<"  "<<mc_mul<<"  "<<mcl_hadfs<<endl;
   //if(trkl_redchi2 > 20 )                                           continue;
   //if(trkl_pass == 0 )                                              continue;
   //if(TMath::Abs(trkl_ubegplane - trkl_vbegplane)> 3 )              continue;
    /*if(trkl_eqp!=0){
	if ((TMath::Abs(trkl_qoverp/trkl_eqp))<3)                      continue; 
	}*/

      //energy cut
      //if (trkl_rangemom<1) continue;


      if(trkl_qoverp<0)
        neus->Fill(); 


      if(trkl_qoverp>0)
       aneus->Fill(); 


    }

 

      nbytes += nb;

  }                                       //end of loop over the chain

      f2->Write();
      f2->Close();

      f4->Write();
      f4->Close();

      ndfile.close();    

      ndfile<<"\n\n end of loop\n"<<endl;           
                    

}


const Float_t cnt::XCurvature(const Int_t begplane,const Int_t endplane,
                         const Int_t fitLengths[500][500],
                         const Int_t trkl_nstps,
                         const Float_t trkl_stpu[],
                         const Float_t trkl_stpv[],
                         const Float_t trkl_stpx[],
                         const Float_t trkl_stpy[],
	    	         const Float_t trkl_stpz[]) const{

  ofstream ndfile;
  ndfile.open("neutrino.txt", ios::app);

	const Int_t trkLthPln =  endplane -  begplane;

	if (-1 == endplane || -1 == begplane){return -1;}

        ndfile<<"inside xcurvature "<<trkLthPln<<endl;	

        Int_t fitLength = fitLengths[trkLthPln][0];
        if (!fitLength) {return 0.0;}

        ndfile<<"fitLength is "<<fitLength<<" , sntrips  "<<trkl_nstps<<endl;

	for(int k=0;k<trkl_nstps;k++){
       
	  ndfile<<k<<"  "<<trkl_stpx[k]<<"  "<<trkl_stpy[k]<<"  "<<trkl_stpz[k]<<endl;

	}

        Float_t xCurv = CalcXCurv(fitLength,trkl_nstps,trkl_stpu,trkl_stpv,trkl_stpx,trkl_stpy,trkl_stpz);

	return 0;

}

const Float_t cnt::CalcXCurv(const Int_t fitLength,
			const Int_t trkl_nstps,
			const Float_t trkl_stpu[],
			const Float_t trkl_stpv[],
			const Float_t trkl_stpx[],
			const Float_t trkl_stpy[],
			const Float_t trkl_stpz[]) const{

     const std::multimap<Float_t , Float_t> md ;

     //ArrayXAlongTrack(trkl_nstps,trkl_stpu,trkl_stpv,trkl_stpx,trkl_stpy,trkl_stpz);

    // return CalcCurv(md, fitLength);

  return 0;
}


const multimap<Float_t, Float_t>  cnt::ArrayXAlongTrack(const Int_t trkl_nstps,
						   const Float_t trkl_stpu[],
						   const Float_t trkl_stpv[],
						   const Float_t trkl_stpx[],
						   const Float_t trkl_stpy[],
						   const Float_t trkl_stpz[]) const{
  multimap<Float_t, Float_t> md;
  md.clear();

  for (Int_t k=0; k<trkl_nstps; ++k){

    md.insert(pair<Float_t,Float_t>(trkl_stpz[k],
                                    trkl_stpx[k]));
  }
  return md;
}



const Float_t cnt::YCurvature(const Int_t begplane,const Int_t endplane,
                         const Int_t fitLengths[500][500],
			 const Int_t trkl_nstps,
		         const Float_t trkl_stpu[],
                         const Float_t trkl_stpv[],
                         const Float_t trkl_stpx[],
                         const Float_t trkl_stpy[],
                         const Float_t trkl_stpz[]) const{

  ofstream ndfile;
  ndfile.open("neutrino.txt", ios::app);

        const Int_t trkLthPln =  endplane -  begplane;

        if (-1 == endplane || -1 == begplane){return -1;}

        ndfile<<"inside ycurvature "<<trkLthPln<<endl;

        Int_t fitLength = fitLengths[trkLthPln][0];
        if (!fitLength) {return 0.0;}

	ndfile<<"fitLength is "<<fitLength<<endl;

        return 0;

}


/*const Float_t cnt::CalcCurv(const multimap<Float_t,Float_t> md,
                       const Int_t fitLength) const
{

  //Set up some arrays.
  const Int_t size = md.size();
  Float_t* xerr = new Float_t[fitLength];
  Float_t* zerr = new Float_t[fitLength];
  if (size<fitLength){return 0.0;}

  multimap<Float_t, Float_t> mfits; //<abs(parabolic fit),(parabolic fit)>
  multimap<Float_t,Float_t>::const_iterator it = md.begin();
  //k loops over the fit segments.
  for (Int_t k=0; k<(size-fitLength); ++k){
    Float_t* z = new Float_t[fitLength];
    Float_t* x = new Float_t[fitLength];
    Float_t* absx = new Float_t[fitLength];
    Bool_t changesign = false;
    Bool_t positive = false;
    if (it->second > 0){positive = true;}
    //l loops over the strips in each fit segment
    //The variable 'positive' says whether the first strip in a
    //segment has a positive co-ordinate.
    for (Int_t l=0; l<fitLength; ++l)
      {
        Bool_t thispos = false;
        z[l] = it->first;
        zerr[l] = 0.0594;
        x[l] = it->second;
        xerr[l] = 0.041;
        absx[l] = fabs(it->second);
        if (x[l] > 0){thispos = true;}
        if (thispos != positive){changesign = true;}
        ++it;
      }
      if (!changesign){
      TGraphErrors g(fitLength,z,absx,zerr,xerr);
      //TCanvas cd("cd","",0,0,500,500);
      //g.Draw("A*");
      g.Fit("pol2","Q");
      Float_t curviness = 2*(g.GetFunction("pol2")->GetParameter(2));
      //      hCurv.Fill(curviness);
      mfits.insert(pair<Float_t,Float_t>(fabs(curviness),
                                         curviness));
    }
    for (Int_t l=0; l<fitLength; ++l)
      {
        --it;
      }
    ++it;
    if(z){delete z; z=0;}
    if(x){delete x; x=0;}
    if(absx){delete absx; absx=0;}
  }

  multimap<Float_t,Float_t>::const_iterator itfits = mfits.end();
  Int_t numRemove = 0;
  if (((Int_t) mfits.size()) >= numRemove){
    for (Int_t rem = 0; rem < numRemove; ++rem){
      itfits--;
      mfits.erase(itfits->first);
      itfits = mfits.end();
    }
  }
  TH1F hCurv("hCurv","",100,-5,5);
  for (itfits = mfits.begin();
       itfits != mfits.end();
       ++itfits){
    hCurv.Fill(itfits->second);
  }

  Float_t neg = hCurv.Integral(0,50);
  Float_t pos = hCurv.Integral(51,101);
  hvariances->Fill(hCurv.GetRMS());
  //  flastVariance = hCurv.GetRMS();
  if (pos > neg) {if (neg) {return pos/neg;} else {return 10.0;}}
  if (neg > pos) {if (pos) {return -1.0*neg/pos;} else {return -10.0;}}
  if (neg == pos) {return 0.0;}
  else {return 0.0;}

}
*/