#include <iostream>
#include <fstream>
#include <iomanip>
#include <signal.h>
#include <TROOT.h>
#include <TFile.h>
#include <TNtuple.h>
#include "TCanvas.h"
#include "TGaxis.h"
#include "TF1.h"
#include "TH1.h"
#include <clasEvent.h>
#include <eloss.h>
#include "mfuncs.h"
#include "TNtupleUtil.h"
#include "SIMULATION/KinLine.C"
#include "SIMULATION/Kstream.C"
#include "alghoulstuff/g12_pcor_may11.h"
#include "g12ecor/g12_ecor.hpp"



int StartRun(int runNo)
{  
  int static CurrentRun = -1;
  if (CurrentRun != runNo) {
    vertex_brun(runNo);
    CurrentRun = runNo;
    //Comment out line below for removing eloss
    InitEloss(CurrentRun);
  }
  return 0;
}

void PrintUsage(char *processName)
{
  cerr << processName << " <options> <filename>\n";
  cerr << "\toptions are:\n";
  cerr << "\t-o<name>.root\tRoot output file.\n";
  cerr << "\t-h\tprint the above" << endl;
}


static int realexit = 0;
void ctrlCHandle(int x)
{
  signal(SIGINT, ctrlCHandle);
  signal(SIGHUP, ctrlCHandle);
  if(realexit > 0) {
    cerr << endl << endl << "\t\t\t***  INTERRUPTED!!!  Exiting... ***"
	 << endl << endl;
    exit(1);
  }
  realexit++;
  cerr << endl << endl << "\t\t\t***  INTERRUPTED!!! Ending Event processing ***"
	 << endl << endl;
}


TROOT ClasNtuple("ClasNtuple","A CLAS Ntuple");

int main(int argc,char **argv){
  
  signal(SIGINT, ctrlCHandle);
  signal(SIGHUP, ctrlCHandle);
  unsigned long EventCount=0, bookCount=0, diffbeam=0 ;
  int partbank0 = 1;
  
  int ret = 1;
  long Nevents = 0;
  char *argptr;
  
  // bos stuff
  int MaxBanks = 1000; 
  //tambahan

 
  Char_t outputfile[50], ntpTitle[50]; 
  sprintf(outputfile,"ClasNtuple.root");
  sprintf(ntpTitle,"Hall  Ntuple");
 
  for (int i = 1; i < argc; i++) {
    argptr = argv[i];
    if (*argptr == '-') {
      argptr++;
      switch (*argptr) {
      case 'o':
	sprintf(outputfile,"%s.root",++argptr);
	cerr<<"Saving output in file: "<<outputfile<<endl;
	break;
      case 'h':
	PrintUsage(argv[0]);
	exit(0);
	break;
      default:
	cerr << "Unrecognized argument: " << argptr << endl;;
	break;
      }
    }
  }///////////////////////////////////////////////////////////


  //
  //  Open the output file.
  //
  Int_t comp   = 1; // compress levels from 0 to 9 (none to most)
  TFile *ntpfile = new TFile(outputfile,"RECREATE","A CLAS ntuple");
  ntpfile->SetCompressionLevel(comp);

  Bool_t firstEvent=kTRUE;

  //
  // create histograms
  //
  TH1F *pull_2pi[10], *pull_3pi[10];
  TH1F *cl_2pi = new TH1F("cl_2pi","cl_2pi",220,0.0,1.1);
  TH1F *cl_3pi = new TH1F("cl_3pi","cl_3pi",220,0.0,1.1);
  pull_2pi[0] = new TH1F("pull_photon_2pi","pull_photon_2pi",100,-5.0,5.0);
  pull_2pi[1] = new TH1F("pull_pip_mom_2pi","pull_pip_mom_2pi",100,-5.0,5.0);
  pull_2pi[2] = new TH1F("pull_pip_lambda_2pi","pull_pip_lambda_2pi",100,-5.0,5.0);
  pull_2pi[3] = new TH1F("pull_pip_phi_2pi","pull_pip_phi_2pi",100,-5.0,5.0);
  pull_2pi[4] = new TH1F("pull_pim_mom_2pi","pull_pim_mom_2pi",100,-5.0,5.0);
  pull_2pi[5] = new TH1F("pull_pim_lambda_2pi","pull_pim_lambda_2pi",100,-5.0,5.0);
  pull_2pi[6] = new TH1F("pull_pim_phi_2pi","pull_pim_phi_2pi",100,-5.0,5.0);
  pull_2pi[7] = new TH1F("pull_proton_mom_2pi","pull_proton_mom_2pi",100,-5.0,5.0);
  pull_2pi[8] = new TH1F("pull_proton_lambda_2pi","pull_proton_lambda_2pi",100,-5.0,5.0);
  pull_2pi[9] = new TH1F("pull_proton_phi_2pi","pull_proton_phi_2pi",100,-5.0,5.0);
  pull_3pi[0] = new TH1F("pull_photon_3pi","pull_photon_3pi",100,-5.0,5.0);
  pull_3pi[1] = new TH1F("pull_pip_mom_3pi","pull_pip_mom_3pi",100,-5.0,5.0);
  pull_3pi[2] = new TH1F("pull_pip_lambda_3pi","pull_pip_lambda_3pi",100,-5.0,5.0);
  pull_3pi[3] = new TH1F("pull_pip_phi_3pi","pull_pip_phi_3pi",100,-5.0,5.0);
  pull_3pi[4] = new TH1F("pull_pim_mom_3pi","pull_pim_mom_3pi",100,-5.0,5.0);
  pull_3pi[5] = new TH1F("pull_pim_lambda_3pi","pull_pim_lambda_3pi",100,-5.0,5.0);
  pull_3pi[6] = new TH1F("pull_pim_phi_3pi","pull_pim_phi_3pi",100,-5.0,5.0);
  pull_3pi[7] = new TH1F("pull_proton_mom_3pi","pull_proton_mom_3pi",100,-5.0,5.0);
  pull_3pi[8] = new TH1F("pull_proton_lambda_3pi","pull_proton_lambda_3pi",100,-5.0,5.0);
  pull_3pi[9] = new TH1F("pull_proton_phi_3pi","pull_proton_phi_3pi",100,-5.0,5.0);

  //
  // setup a TNtuple
  //
  TNtuple *ntp=0;
  TntpLables *vnames= new TntpLables();
  Int_t n_vectors;
  Float_t values[500];
  Char_t label[50];

  //
  // Initialize BOS
  bnames_(&MaxBanks);
  initbos();
  configure_banks(stderr,0);

  read_pcor( pcor_array );

  for (int i = 1; i < argc; ++i) {
    argptr = argv[i];
    // process all arguments on command line.
    if (*argptr != '-' && realexit==0 ) {
      // we have a file to process
      EventCount = 0;
      diffbeam = 0;
      cerr<<"\nWe have a file to process ...\n\n\n\n\n\n";
      bookCount =0;
      clasEvent event(argptr,&bcs_,partbank0,0);
      
      if (event.status() ) 
	ret = 1;
      
      while ( ret && realexit==0 ) {
	// process every event in the file
	ret = event.read(partbank0);
	//cerr<<"Event status: "<<event.status()<<endl;
	if(DISIO_DATA == ret && event.status() ){
	  Nevents++;
	  int runno = event.run();
          int eventno = event.event();	
	  ConfigEvent(runno,0); 
	  StartRun(runno);

	  //
	  // ProcessOneEvent
	  n_vectors=0;
	  int nReturnValue; nReturnValue = 0;
	  if(event.load(partbank0)==0) cerr<<"Error: event.load(partbank0)\n";
	  
	  //  clasTAGR_t *tagR = 0;
	  //if(tagR = event.TAGR()) {
	  // }
	  
	  
	  //
	  // Book selected event channel
	  // 
          if(EventCount<10)
            {
	      //cerr<<event.N(PiMinus)<<event.N(-1)<<event.N(PiPlus)<<event.N(Proton)<<event.N(1)<<"\n\n\n\n\n\n\n\n\n";
            }
	  if(event.N(PiMinus)==1 && event.N(-1)==1   
	      && event.N(PiPlus)==1 && event.N(Proton)==1
	      && event.N(1)==2){

            //cerr<<"It found an event \n\n\n\n\n";

	    fourVec pim, pip, proton, beam,beam2, mm,mm2,gammas;
	    fourVec target( PROTON_MASS, threeVec(0.0,0.0,0.0) );
	    threeVec vert = event.V();
            fourVec gamma0(0.0, threeVec(0.0,0.0,0.0) );
            fourVec temp_proton, temp_pip, temp_pim, temp_miss4v, temp_beam,temp_beam2,temp_miss4v2;

            fourVec wo_proton, wo_pip, wo_pim, wo_mm, wo_beam,wo_beam2, el_proton,el_pip,el_pim,el_mm,el_mm2,el_beam,el_beam2;
	    
            //get clas part
            clasParticle partpim,partpip,partproton,partbeam;
            partpim = event.cp(PiMinus,1);
            partpip = event.cp(PiPlus,1);
            partproton = event.cp(Proton,1);

	    
	    //  Fill 4-vector of each particle.
	    //
	    wo_pim = partpim.p();
	    wo_pip = partpip.p();
	    wo_proton = partproton.p();
	    wo_beam = event.beam().get4P();


            //beam2 = event.getBeam(4.4,5.7);
            int ediff = 0;
            if(fabs(wo_beam.t() - wo_beam2.t()) > .001)
            {
              diffbeam++;
              ediff = 1;
            }
            ////Calling Eloss - Temporary Storage////
            vector3_t vtx = makevector3_t(vert);
            temp_beam = wo_beam;

            temp_pim = makefourVec(c_momcor(makevector4_t(wo_pim), PI_CHARGED_MASS, vtx, 1, 7));
            temp_pip = makefourVec(c_momcor(makevector4_t(wo_pip), PI_CHARGED_MASS, vtx, 1, 7));
            temp_proton = makefourVec(c_momcor(makevector4_t(wo_proton), PROTON_MASS, vtx, 1, 7));
            temp_miss4v = temp_beam + target - temp_proton - temp_pim - temp_pip;

            el_pim = temp_pim;
            el_pip = temp_pip;
            el_proton = temp_proton;
            el_beam = temp_beam;
            el_mm = temp_miss4v;




            ///////////////// Momentum corrections

            TLorentzVector pc_P4vector[6],el_lproton, el_lpip, el_lpim;

            el_lproton.SetXYZT(el_proton.x(),el_proton.y(),el_proton.z(),el_proton.t());
            el_lpip.SetXYZT(el_pip.x(),el_pip.y(),el_pip.z(),el_pip.t());
            el_lpim.SetXYZT(el_pim.x(),el_pim.y(),el_pim.z(),el_pim.t());

            //pc_P4vector[0] = Calc_new_fourVec(el_lproton, 0);
            //pc_P4vector[1] = Calc_new_fourVec(el_lpip, 3);
            //pc_P4vector[2] = Calc_new_fourVec(el_lpim, 4);

            pc_P4vector[0].SetXYZT(el_proton.x(),el_proton.y(),el_proton.z(),el_proton.t());
	    pc_P4vector[1].SetXYZT(el_pip.x(),el_pip.y(),el_pip.z(),el_pip.t());
            pc_P4vector[2].SetXYZT(el_pim.x(),el_pim.y(),el_pim.z(),el_pim.t());

            //cout<<pc_P4vector[0].T()<<pc_P4vector[0].X()<<pc_P4vector[0].Y()<<pc_P4vector[0].Z()<<endl;

            proton.set(pc_P4vector[0].T(),pc_P4vector[0].X(),pc_P4vector[0].Y(),pc_P4vector[0].Z());
            pip.set(pc_P4vector[1].T(),pc_P4vector[1].X(),pc_P4vector[1].Y(),pc_P4vector[1].Z());
            pim.set(pc_P4vector[2].T(),pc_P4vector[2].X(),pc_P4vector[2].Y(),pc_P4vector[2].Z());
            beam = el_beam;

            mm = beam + target - proton - pim - pip;

           // cout<<proton.x()<<proton.y()<<proton.z()<<proton.t()<<endl;
           // cout<<el_proton.x()<<el_proton.y()<<el_proton.z()<<el_proton.t()<<endl;


            ////////////////// BEC

            double bec_beame = clas6::g12::corrected_beam_energy(runno, el_beam.t());
            double bec_beame2 = clas6::g12::corrected_beam_energy(runno, el_beam2.t());

            fourVec bec_beam(bec_beame,threeVec(0.0,0.0,bec_beame));


            //if(abs(bec_beam.z() - bec_beam.t())>.0001) cout<<"WRONG BEC BEAM ASSIGNMENT!";

            mm = bec_beam + target -proton -pim - pip;

            beam = bec_beam;

            gammas = gamma0;
            clasParticle CPg1, CPg2;
            if (event.N(Gamma)==1)
            {
              CPg1 = event.cp(Gamma, 1);
              gammas = CPg1.p();
            }
            if(event.N(Gamma)>1)
            {
              CPg1 = event.cp(Gamma, 1);
              CPg2 = event.cp(Gamma, 2);
              gammas = CPg1.p() + CPg2.p();
            }

            //////////////////////////////////////////////////////////////////
            ///////////////////////////Start Kinematic Fitting////////////////
            //////////////////////////////////////////////////////////////////
            const int DK_nparts = 3;  //number of particles in fit
            vector<TLorentzVector> fourVin(DK_nparts);
            vector<TLorentzVector> wo_eloss_fourVin(DK_nparts);
            vector<TLorentzVector> w_eloss_fourVin(DK_nparts);
            vector<TLorentzVector> bec_fourVin(DK_nparts);
            vector<TLorentzVector> bec_fourVin2(DK_nparts);
            vector<TLorentzVector> nop_fourVcor(DK_nparts);
            vector<TLorentzVector> pi0_fourVcor(DK_nparts);
            vector<TLorentzVector> bec_pi0_fourVcor(DK_nparts);

            vector<TLorentzVector> eta_fourVcor(DK_nparts);
            std::vector<TVector3> Vert(DK_nparts);
            std::vector<string> p_names(DK_nparts);
            bool Multi = false;
            bool MC = false;
            Double_t m_targ = PROTON_MASS;
            Double_t e_gamma = wo_beam.t();
            const int dim = 3 * DK_nparts + 1;
            TMatrixD covMatrix(dim, dim), woecovMatrix(dim,dim),wecovMatrix(dim,dim), beccovMatrix(dim,dim), beccovMatrix2(dim,dim);
            TMatrixD covTrack(dim, dim);
            clasTBER_t* TBER(NULL);
            clasPART_t* PART(NULL);
            clasTBID_t* TBID(NULL);
            clasSCRC_t* SCRC(NULL);
clasTGBI_t* TGBI(NULL);
clasHEAD_t* HEAD(NULL);
            if (!(PART = static_cast<clasPART_t*>(getGroup(&bcs_, "PART", 1)))) continue;
            if (!(TBID = static_cast<clasTBID_t*>(getBank(&bcs_, "TBID")))) continue;
            if (!(TBER = static_cast<clasTBER_t*>(getBank(&bcs_, "TBER")))) continue;
            if (!(SCRC = static_cast<clasSCRC_t*>(getBank(&bcs_, "SCRC")))) continue;
if (!(TGBI = static_cast<clasTGBI_t*>(getBank(&bcs_, "TGBI")))) continue;
if (!(HEAD = static_cast<clasHEAD_t*>(getBank(&bcs_, "HEAD")))) continue;
            for (int i = 0; i < TBID->bank.nrow; i++)
            {
              int track = TBID->tbid[i].track - 1;
              if (!(PART->part[i].q)) continue;
	      if (PART->part[i].pid == 8)
	      {
		double p_mag = pip.V().len();
		float mq = PART->part[i].q;
		covTrack(0, 0) = pow((0.001*5.715),2)/3;
		covTrack(1, 1) = (TBER->tber[track].c11) * pow(p_mag, 4);
		covTrack(1, 2) = -(TBER->tber[track].c12) * mq * p_mag * p_mag;
		covTrack(1, 3) = -(TBER->tber[track].c13) * mq * p_mag * p_mag;
		covTrack(2, 1) = -(TBER->tber[track].c12) * mq * p_mag * p_mag;
		covTrack(2, 2) = (TBER->tber[track].c22);
		covTrack(2, 3) = (TBER->tber[track].c23);
		covTrack(3, 1) = -(TBER->tber[track].c13) * mq * p_mag * p_mag;
		covTrack(3, 2) = (TBER->tber[track].c23);
		covTrack(3, 3) = (TBER->tber[track].c33);
	      }
	      if (PART->part[i].pid == 9)
	      {
		double p_mag = pim.V().len();
		float mq = PART->part[i].q;
		covTrack(4, 4) = (TBER->tber[track].c11) * pow(p_mag, 4);
		covTrack(4, 5) = -(TBER->tber[track].c12) * mq * p_mag * p_mag;
		covTrack(4, 6) = -(TBER->tber[track].c13) * mq * p_mag * p_mag;
		covTrack(5, 4) = -(TBER->tber[track].c12) * mq * p_mag * p_mag;
		covTrack(5, 5) = (TBER->tber[track].c22);
		covTrack(5, 6) = (TBER->tber[track].c23);
		covTrack(6, 4) = -(TBER->tber[track].c13) * mq * p_mag * p_mag;
		covTrack(6, 5) = (TBER->tber[track].c23);
		covTrack(6, 6) = (TBER->tber[track].c33);
	      }
	      if (PART->part[i].pid == 14)
	      {
		double p_mag = proton.V().len();
		float mq = PART->part[i].q;
		covTrack(7, 7) = (TBER->tber[track].c11) * pow(p_mag, 4);
		covTrack(7, 8) = -(TBER->tber[track].c12) * mq * p_mag * p_mag;
		covTrack(7, 9) = -(TBER->tber[track].c13) * mq * p_mag * p_mag;
		covTrack(8, 7) = -(TBER->tber[track].c12) * mq * p_mag * p_mag;
		covTrack(8, 8) = (TBER->tber[track].c22);
		covTrack(8, 9) = (TBER->tber[track].c23);
		covTrack(9, 7) = -(TBER->tber[track].c13) * mq * p_mag * p_mag;
		covTrack(9, 8) = (TBER->tber[track].c23);
		covTrack(9, 9) = (TBER->tber[track].c33);
	      }
	    }
            /*
	    clasPART_t* PARTy(NULL);
	    int seclist[3], sc_id[3], sc_id2[3];
	    char trypart[4][10] = { "Proton", "PiPlus", "PiMinus", "Neutral" };
	    PARTy = static_cast<clasPART_t*>(getGroup(&bcs_, "PART", 1));
	    for (int qq = 0; qq < TBID->bank.nrow; qq++)
	    {
	      if (PARTy->part[qq].pid == 14)
	      {
		seclist[0] = TBID->tbid[qq].sec;
                sc_id[0] = TBID->tbid[qq].sc_id;
                sc_id2[0] = SCRC->scrc[qq].id;

                if(sc_id[0] != sc_id2[0]) cout<<"CHECK SC_ID";
	      }
	      if (PARTy->part[qq].pid == 8)
	      {
		seclist[1] = TBID->tbid[qq].sec;
                sc_id[1] = TBID->tbid[qq].sc_id;
                sc_id2[1] = SCRC->scrc[qq].id;

                if(sc_id[1] != sc_id2[1]) cout<<"CHECK SC_ID";
	      }
	      if (PARTy->part[qq].pid == 9)
	      {
		seclist[2] = TBID->tbid[qq].sec;
                sc_id[2] = TBID->tbid[qq].sc_id;
                sc_id2[3] = SCRC->scrc[qq].id;

                if(sc_id[2] != sc_id2[2]) cout<<"CHECK SC_ID";
	      }
	    }
            */

            int sec[3],sc_id_part[3];
            sc_id_part[0] = partproton.scPaddleId();
            sc_id_part[1] = partpip.scPaddleId();
            sc_id_part[2] = partpim.scPaddleId();
            sec[0] = partproton.sec();
            sec[1] = partpip.sec();
            sec[2] = partpim.sec();
int pim_test = partpim.scPaddleId();
            

	    fourVin[0].SetXYZT(pip.x(), pip.y(), pip.z(), pip.t());
	    fourVin[1].SetXYZT(pim.x(), pim.y(), pim.z(), pim.t());
	    fourVin[2].SetXYZT(proton.x(), proton.y(), proton.z(), proton.t());
            wo_eloss_fourVin[0].SetXYZT(wo_pip.x(), wo_pip.y(), wo_pip.z(), wo_pip.t());
	    wo_eloss_fourVin[1].SetXYZT(wo_pim.x(), wo_pim.y(), wo_pim.z(), wo_pim.t());
	    wo_eloss_fourVin[2].SetXYZT(wo_proton.x(), wo_proton.y(), wo_proton.z(), wo_proton.t());
            w_eloss_fourVin[0].SetXYZT(el_pip.x(), el_pip.y(), el_pip.z(), el_pip.t());
	    w_eloss_fourVin[1].SetXYZT(el_pim.x(), el_pim.y(), el_pim.z(), el_pim.t());
	    w_eloss_fourVin[2].SetXYZT(el_proton.x(), el_proton.y(), el_proton.z(), el_proton.t());

	    bec_fourVin[0].SetXYZT(pip.x(), pip.y(), pip.z(), pip.t());
	    bec_fourVin[1].SetXYZT(pim.x(), pim.y(), pim.z(), pim.t());
            bec_fourVin[2].SetXYZT(proton.x(), proton.y(), proton.z(), proton.t());


	    std::vector<bool> set(DK_nparts);
	    set[0] = true;
	    set[1] = true;
	    set[2] = false;
	    Vert[0].SetXYZ(partpip.x(), partpip.y(), partpip.z());
	    Vert[1].SetXYZ(partpim.x(), partpim.y(), partpim.z());
	    Vert[2].SetXYZ(partproton.x(), partproton.y(), partproton.z());
	    p_names[0] = "pi+";
	    p_names[1] = "pi-";
	    p_names[2] = "p";
	    string RP = "g12";
	    covMatrix = CorrectCLAS_V(covTrack, p_names, fourVin, Vert, Multi, MC, RP);
            woecovMatrix = CorrectCLAS_V(covTrack,p_names, wo_eloss_fourVin, Vert, Multi, MC,RP);
            wecovMatrix = CorrectCLAS_V(covTrack,p_names, w_eloss_fourVin, Vert, Multi, MC,RP);
	    beccovMatrix = CorrectCLAS_V(covTrack, p_names, fourVin, Vert, Multi, MC, RP);

	    Kstream PhiFit;

            Kstream woelossFit,welossFit,becFit;

            //cout<< "covTrack(0,0) : covMatrix(0,0) ::"<<covTrack(0,0)<<" : "<<covMatrix(0,0)<<endl;

            //covTrack(0,0) = covMatrix(0,0);

	    PhiFit.StringNames(p_names);
	    welossFit.StringNames(p_names);
	    woelossFit.StringNames(p_names);
	    becFit.StringNames(p_names);


	    PhiFit.FitInput(e_gamma, fourVin, covMatrix, m_targ);
	    welossFit.FitInput(e_gamma, w_eloss_fourVin, wecovMatrix, m_targ);
	    woelossFit.FitInput(e_gamma, wo_eloss_fourVin, woecovMatrix, m_targ);
            becFit.FitInput(bec_beame,bec_fourVin,beccovMatrix,m_targ);



	    double pi0_pulls[10] = {};
	    double wo_eloss_nop_pulls[10] = {};
	    double w_eloss_nop_pulls[10] = {};
	    double bec_nop_pulls[10] = {};
	    double wo_eloss_pi_pulls[10] = {};
	    double w_eloss_pi_pulls[10] = {};
	    double bec_pi0_pulls[10] = {};

	    double nop_pulls[10] = {};
	    double eta_pulls[10] = {};
	    bool include = false;

	    int mm_pid=3;
	    float wo_eloss_nop_prob = 0.0;
	    float w_eloss_nop_prob = 0.0;
	    float bec_nop_prob = 0.0;
            float wo_eloss_pi_prob = 0.0;
            float w_eloss_pi_prob = 0.0;
	    float nop_prob = 0.0;
	    float pi0_prob = 0.0;
	    float eta_prob = 0.0;
	    float bec_pi0_prob = 0.0;
 
            float bec_pi0_fitEnergy = 0.0;

            woelossFit.Fit();
            wo_eloss_nop_prob = woelossFit.Prob();
            for (int i = 0; i < 10; i++) wo_eloss_nop_pulls[i] = woelossFit.GetPull(i);

            becFit.Fit();
            bec_nop_prob = becFit.Prob();
            cl_2pi->Fill(bec_nop_prob);
            for (int i = 0; i < 10; i++) {
	      bec_nop_pulls[i] = becFit.GetPull(i);

              if (partproton.nGammaRF() == 1 && partpip.nGammaRF() == 1 && partpim.nGammaRF() == 1 && bec_nop_prob > 0.005) {
              //if (bec_nop_prob > 0.005) {
		pull_2pi[i]->Fill(becFit.GetPull(i));
	      }
	    }

            becFit.Fit(0.134976);
            bec_pi0_prob = becFit.Prob();
            cl_3pi->Fill(bec_pi0_prob);
            for (int i = 0; i < 10; i++) {
	      bec_pi0_pulls[i] = becFit.GetPull(i);
              pull_3pi[i]->Fill(becFit.GetPull(i));
	    }
            for (int i = 0; i < 3; i++)
            {
              bec_pi0_fourVcor[i] = becFit.FitP4(i);
            }
            bec_pi0_fitEnergy = becFit.FitPhotonEnergy();

            woelossFit.Fit(0.134976);
            wo_eloss_pi_prob = woelossFit.Prob();
            for (int i = 0; i < 10; i++) wo_eloss_pi_pulls[i] = woelossFit.GetPull(i);

            welossFit.Fit();
            w_eloss_nop_prob = welossFit.Prob();
            for (int i = 0; i < 10; i++) w_eloss_nop_pulls[i] = welossFit.GetPull(i);

            welossFit.Fit(0.134976);
            w_eloss_pi_prob = welossFit.Prob();
            for (int i = 0; i < 10; i++) w_eloss_pi_pulls[i] = welossFit.GetPull(i);

            PhiFit.Fit();
	    nop_prob = PhiFit.Prob();
	    for (int i = 0; i < 10; i++) nop_pulls[i] = PhiFit.GetPull(i);
	    for (int i = 0; i < 3; i++)
	    {
	      nop_fourVcor[i] = PhiFit.FitP4(i);
	    }
	    
	    PhiFit.Fit(0.134976);
	    pi0_prob = PhiFit.Prob();
	    for (int i = 0; i < 10; i++) pi0_pulls[i] = PhiFit.GetPull(i);
	    for (int i = 0; i < 3; i++)
	    {
	      pi0_fourVcor[i] = PhiFit.FitP4(i);
	    }

            PhiFit.Fit(0.547853);
            eta_prob = PhiFit.Prob();
	    for (int i = 0; i < 10; i++) eta_pulls[i] = PhiFit.GetPull(i);
	    for (int i = 0; i < 3; i++)
	    {
	      eta_fourVcor[i] = PhiFit.FitP4(i);
	    }

	    //cout << "temprob1 :: "<< tempprob1 << endl;
	    //cout << "temprob2 :: "<< tempprob2 << endl;

            //cout << "MM : prob none : prob pi0 : prob eta :: " << ~(mm) << " : " << nop_prob << " : "<<
            //pi0_prob << " : "<< eta_prob << endl;

           // cout<<" pull0 wo : w : nop : pi0 : eta :: "<< wo_eloss_nop_pulls[0] <<" : "/*<< w_eloss_nop_pulls[0]*/ <<" : "
            //<< nop_pulls[0] <<" : "<< pi0_pulls[0] <<" : "<< eta_pulls[0] <<endl;
  
            //TEST
/*            if(fabs(mm.t() +proton.t() +pim.t() +pip.t() -beam.t() -target.t())>0.001 ){
		cerr<<"Something is wrong with the energies!"<<endl;
                break;
            }
            if(firstEvent){
                cout<<"Proton Energy: "<<proton.t()<<"Proton rest mass: "<<target.t()<<endl;
            }
*/












            /////////////////////////////////// SC Paddle
            /*
	    values[n_vectors++] = seclist[0];
	    if(firstEvent){
	      sprintf(label,"p_sec_alt");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }


	    values[n_vectors++] = sc_id[0];
	    if(firstEvent){
	      sprintf(label,"p_sc_id");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }


	    values[n_vectors++] = seclist[1];
	    if(firstEvent){
	      sprintf(label,"pip_sec_alt");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }


	    values[n_vectors++] = sc_id[1];
	    if(firstEvent){
	      sprintf(label,"pip_sc_id");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }

	    values[n_vectors++] = seclist[2];
	    if(firstEvent){
	      sprintf(label,"pim_sec_alt");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }


	    values[n_vectors++] = sc_id[2];
	    if(firstEvent){
	      sprintf(label,"pim_sc_id");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }

	    values[n_vectors++] = sc_id2[0];
	    if(firstEvent){
	      sprintf(label,"p_sc_id2");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }


	    values[n_vectors++] = sc_id2[1];
	    if(firstEvent){
	      sprintf(label,"pip_sc_id2");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }

	    values[n_vectors++] = sc_id2[2];
	    if(firstEvent){
	      sprintf(label,"pim_sc_id2");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }
            */

values[n_vectors++] = pim_test;
	    if(firstEvent){
	      sprintf(label,"pim_Testt");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }

	    values[n_vectors++] = sc_id_part[0];
	    if(firstEvent){
	      sprintf(label,"p_sc_id_part");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }


	    values[n_vectors++] = sc_id_part[1];
	    if(firstEvent){
	      sprintf(label,"pip_sc_id_part");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }

	    values[n_vectors++] = sc_id_part[2];
	    if(firstEvent){
	      sprintf(label,"pim_sc_id_part");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }

	    values[n_vectors++] = sec[0];
	    if(firstEvent){
	      sprintf(label,"p_sec");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }


	    values[n_vectors++] = sec[1];
	    if(firstEvent){
	      sprintf(label,"pip_sec");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }

	    values[n_vectors++] = sec[2];
	    if(firstEvent){
	      sprintf(label,"pim_sec");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }

            /////////////////////////////// Photon info
          
            int ngrf, p_ngrf, pip_ngrf, pim_ngrf, trigbit;
//Int_t trigbits = TGBI->tgbi[0].latch1;

            /*clasGPID_t* GPID(NULL);
            clasTAGR_t* TAGR(NULL);
            GPID = static_cast<clasGPID_t*>(getBank(&bcs_, "GPID"));
            TAGR = static_cast<clasTAGR_t*>(getBank(&bcs_, "TAGR"));
            if(GPID){
              cout<<"HAVE GPID!"<<endl;
            }
            int p_pho_ind;
            if(TAGR){
              if(partproton.tagr_id()>0){
                p_pho_ind = partproton.tagr_id() - 1;
              }
              else{
                cout<<"TAGR ID: "<<partproton.tagr_id()<<endl;
                p_pho_ind = 0;
              }
              cout<<p_pho_ind<<endl;
              int tagr_rows = TAGR->bank.nrow;
              cout<<tagr_rows<<endl;
              double tpho_p = TAGR->tagr[p_pho_ind].tpho;
              cout<<tpho_p<<endl;
            }
            else{
              cout<<"NO TAGR!"<<endl;
            }*/

            p_ngrf = partproton.nGammaRF();
            pip_ngrf = partpip.nGammaRF();
            pim_ngrf = partpim.nGammaRF();
            //trigbit = HEAD[O].TRIGBIT;

            if( p_ngrf == pip_ngrf && p_ngrf == pim_ngrf) ngrf = p_ngrf;
            else ngrf= -1;

            /*else{
            p_ngrf = -1;
            pim_ngrf = -1;
            pip_ngrf = -1;
            ngrf = -1;
            }*/

	    values[n_vectors++] = p_ngrf;
	    if(firstEvent){
	      sprintf(label,"p_ngrf");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }

	    values[n_vectors++] = pip_ngrf;
	    if(firstEvent){
	      sprintf(label,"pip_ngrf");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }

	    values[n_vectors++] = pim_ngrf;
	    if(firstEvent){
	      sprintf(label,"pim_ngrf");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }

	    values[n_vectors++] = ngrf;
	    if(firstEvent){
	      sprintf(label,"ngrf");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }
            
            //values[n_vectors++] = TGBI->tgbi[0].latch1;
values[n_vectors++] = HEAD->head[0].trigbits;
	    if(firstEvent){
	      sprintf(label,"trigbit");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }
 //
//values[n_vectors++] = TGBI[0].latch1;
values[n_vectors++] = TGBI->tgbi[0].latch1;
	    if(firstEvent){
	      sprintf(label,"beamPol");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }
//           
	    values[n_vectors++] = partpip.tagr_id();
	    if(firstEvent){
	      sprintf(label,"pip_tagr_id");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }


	    values[n_vectors++] = partpim.tagr_id();
	    if(firstEvent){
	      sprintf(label,"pim_tagr_id");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }


	    values[n_vectors++] = partproton.tagr_id();
	    if(firstEvent){
	      sprintf(label,"p_tagr_id");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
              }
            
            values[n_vectors++] = ediff;
            if(firstEvent){
              sprintf(label,"ediff");
              vnames->Add(label);
              }

            values[n_vectors++] = runno;;
            if(firstEvent){
              sprintf(label,"runNumber");
              vnames->Add(label);
	    }

            values[n_vectors++] = eventno;;
            if(firstEvent){
              sprintf(label,"eventNumber");
              vnames->Add(label);
            }

	    //
	    //  Book event info into ntuple
	    //
	    ////////////////////////////////////////////

            //Kinematic Fitter results
            for(int i = 0; i<10;i++){


              
	      values[n_vectors++] = bec_nop_pulls[i]; 
	      if(firstEvent){
	        sprintf(label,"bec_nop_pull%d",i);
	        vnames->Add(label);
	        //cout<<"beamE:\t\t  Photon beam energy\n";
              }
              
	      values[n_vectors++] = bec_pi0_pulls[i]; 
	      if(firstEvent){
	        sprintf(label,"bec_pi0_pull%d",i);
	        vnames->Add(label);
	        //cout<<"beamE:\t\t  Photon beam energy\n";
              }
              /*values[n_vectors++] = bec_pi0_pulls2[i]; 
              if(firstEvent){
                sprintf(label,"bec_pi0_pull2%d",i);
                vnames->Add(label);
                //cout<<"beamE:\t\t  Photon beam energy\n";
              } */         
              /*
	      values[n_vectors++] = wo_eloss_nop_pulls[i]; 
	      if(firstEvent){
	        sprintf(label,"wo_eloss_nop_pull%d",i);
	        vnames->Add(label);
	        //cout<<"beamE:\t\t  Photon beam energy\n";
              }
          
	      values[n_vectors++] = w_eloss_nop_pulls[i]; 
	      if(firstEvent){
	        sprintf(label,"w_eloss_nop_pull%d",i);
	        vnames->Add(label);
	        //cout<<"beamE:\t\t  Photon beam energy\n";
              }
          
	      values[n_vectors++] = wo_eloss_pi_pulls[i]; 
	      if(firstEvent){
	        sprintf(label,"wo_eloss_pi_pull%d",i);
	        vnames->Add(label);
	        //cout<<"beamE:\t\t  Photon beam energy\n";
              }
          
	      values[n_vectors++] = w_eloss_pi_pulls[i]; 
	      if(firstEvent){
	        sprintf(label,"w_eloss_pi_pull%d",i);
	        vnames->Add(label);
	        //cout<<"beamE:\t\t  Photon beam energy\n";
              }
          
	      values[n_vectors++] = nop_pulls[i]; 
	      if(firstEvent){
	        sprintf(label,"nop_pull%d",i);
	        vnames->Add(label);
	        //cout<<"beamE:\t\t  Photon beam energy\n";
              }

	      values[n_vectors++] = pi0_pulls[i]; 
	      if(firstEvent){
	        sprintf(label,"pi0_pull%d",i);
	        vnames->Add(label);
	        //cout<<"beamE:\t\t  Photon beam energy\n";
              }

	      values[n_vectors++] = eta_pulls[i]; 
	      if(firstEvent){
	        sprintf(label,"eta_pull%d",i);
	        vnames->Add(label);
	        //cout<<"beamE:\t\t  Photon beam energy\n";
              }
              */
            }
            
	    values[n_vectors++] = bec_nop_prob; 
	    if(firstEvent){
	      sprintf(label,"bec_nop_prob");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }
            

	    values[n_vectors++] = bec_pi0_prob; 
	    if(firstEvent){
	      sprintf(label,"bec_pi0_prob");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }        

            /*values[n_vectors++] = bec_pi0_prob2;
            if(firstEvent){
              sprintf(label,"bec_pi0_prob2");
              vnames->Add(label);
            }*/
            /*
	    values[n_vectors++] = wo_eloss_nop_prob; 
	    if(firstEvent){
	      sprintf(label,"wo_eloss_nop_prob");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }
        
	    values[n_vectors++] = w_eloss_nop_prob; 
	    if(firstEvent){
	      sprintf(label,"w_eloss_nop_prob");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }

	    values[n_vectors++] = wo_eloss_pi_prob; 
	    if(firstEvent){
	      sprintf(label,"wo_eloss_pi_prob");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }
        
	    values[n_vectors++] = w_eloss_pi_prob; 
	    if(firstEvent){
	      sprintf(label,"w_eloss_pi_prob");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }
        
	    values[n_vectors++] = nop_prob; 
	    if(firstEvent){
	      sprintf(label,"nop_prob");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }

	    values[n_vectors++] = pi0_prob; 
	    if(firstEvent){
	      sprintf(label,"pi0_prob");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }

	    values[n_vectors++] = eta_prob;
	    if(firstEvent){
	      sprintf(label,"eta_prob");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }
            */
            //test

	    values[n_vectors++] = partpim.scTOF(); 
	    if(firstEvent){
	      sprintf(label,"pim_tof");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }

	    values[n_vectors++] = partpim.scTOFexpected(); 
	    if(firstEvent){
	      sprintf(label,"pim_tofp");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }

	    values[n_vectors++] = partpip.scTOF(); 
	    if(firstEvent){
	      sprintf(label,"pip_tof");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }

	    values[n_vectors++] = partpip.scTOFexpected(); 
	    if(firstEvent){
	      sprintf(label,"pip_tofp");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }

	    values[n_vectors++] = partproton.scTOF(); 
	    if(firstEvent){
	      sprintf(label,"p_tof");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }

	    values[n_vectors++] = partproton.scTOFexpected(); 
	    if(firstEvent){
	      sprintf(label,"p_tofp");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }

	    values[n_vectors++] = partproton.scVtime(); 
	    if(firstEvent){
	      sprintf(label,"p_scvtime");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }

	    values[n_vectors++] = partproton.stVtime(); 
	    if(firstEvent){
	      sprintf(label,"p_stvtime");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }


	    values[n_vectors++] = partpip.scVtime(); 
	    if(firstEvent){
	      sprintf(label,"pip_scvtime");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }

	    values[n_vectors++] = partpip.stVtime(); 
	    if(firstEvent){
	      sprintf(label,"pip_stvtime");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }

	    values[n_vectors++] = partpim.scVtime(); 
	    if(firstEvent){
	      sprintf(label,"pim_scvtime");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }

	    values[n_vectors++] = partpim.stVtime(); 
	    if(firstEvent){
	      sprintf(label,"pim_stvtime");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }


	    values[n_vectors++] = beam.t(); 
	    if(firstEvent){
	      sprintf(label,"beamE");
	      vnames->Add(label);
	      //cout<<"beamE:\t\t  Photon beam energy\n";
	    }

            values[n_vectors++] = bec_pi0_fitEnergy;
            if(firstEvent){
              sprintf(label,"kin_beamE");
              vnames->Add(label);
              //cout<<"beamE:\t\t  Photon beam energy\n";
            }

	    values[n_vectors++] = event.vtime()- event.stVtime() ;
	    if(firstEvent){
	      sprintf(label,"stVtimeDiff");
	      vnames->Add(label);
	      //cout<<"stVtimeDiff:\t\t  vtime - stVtime x\n";
	    }

	    values[n_vectors++] = event.vtime()- event.scVtime() ;
	    if(firstEvent){
	      sprintf(label,"scVtimeDiff");
	      vnames->Add(label);
	      //cout<<"scVtimeDiff:\t\t  vtime - scVtime x\n";
	    }
	    values[n_vectors++] = event.stVtime();
	    if(firstEvent){
	      sprintf(label,"stVtime");
	      vnames->Add(label);
	      //cout<<"stVtime:\t\t Vertex time from Start Counter\n";
	    }
	    values[n_vectors++] = event.scVtime();
	    if(firstEvent){
	      sprintf(label,"scVtime");
	      vnames->Add(label);
	      //cout<<"scVtime:\t\t Vertex time from TOF\n";
	    }
	    values[n_vectors++] = event.vtime();
	    if(firstEvent){
	      sprintf(label,"vtime");
	      vnames->Add(label);
	      //cout<<"vtime:\t\t Vertex time from photon\n";
	    }

	    values[n_vectors++] = vert.x();
	    if(firstEvent){
	      sprintf(label,"x");
	      vnames->Add(label);
	      //cout<<"x:\t\t  Vertex x\n";
	    }
	    values[n_vectors++] = vert.y();
	    if(firstEvent){
	      sprintf(label,"y");
	      vnames->Add(label);
	      //cout<<"y:\t\t  Vertex y\n";
	    }
	    values[n_vectors++] = vert.z();
	    if(firstEvent){
	      sprintf(label,"z");
	      vnames->Add(label);
	      //cout<<"z:\t\t  Vertex z\n";
	    }
	    
	    //
	    // Book particle info
	    //
	    values[n_vectors++] = pim.x();
	    if(firstEvent){
	      sprintf(label,"pim_px");
	      vnames->Add(label);
	      //cout<<"pim_px:\t\t PiMinus x momentum \n";
	    }
	    values[n_vectors++] = pim.y();
	    if(firstEvent){
	      sprintf(label,"pim_py");
	      vnames->Add(label);
	      //cout<<"pim_py:\t\t PiMinus y momentum \n";
	    }
	    values[n_vectors++] = pim.r();
	    if(firstEvent){
	      sprintf(label,"pim_p");
	      vnames->Add(label);
	      //cout<<"pim_p:\t\t PiMinus momentum \n";
	    }
	    values[n_vectors++] = pim.z();
	    if(firstEvent){
	      sprintf(label,"pim_pz");
	      vnames->Add(label);
	      //cout<<"pim_pz:\t\t PiMinus z momentum \n";
	    }
	    values[n_vectors++] = pim.t();
	    if(firstEvent){
	      sprintf(label,"pim_E");
	      vnames->Add(label);
	      //cout<<"pim_E:\t\t PiMinus Energy \n";
	    }
            values[n_vectors++] = bec_pi0_fourVcor[1].X();
            if(firstEvent){
              sprintf(label,"kin_pim_px");
              vnames->Add(label);
              //cout<<"pim_px:\t\t PiMinus x momentum \n";
            }
            values[n_vectors++] = bec_pi0_fourVcor[1].Y();
            if(firstEvent){
              sprintf(label,"kin_pim_py");
              vnames->Add(label);
              //cout<<"pim_py:\t\t PiMinus y momentum \n";
            }
            values[n_vectors++] = bec_pi0_fourVcor[1].Z();
            if(firstEvent){
              sprintf(label,"kin_pim_pz");
              vnames->Add(label);
              //cout<<"pim_pz:\t\t PiMinus z momentum \n";
            }
            values[n_vectors++] = bec_pi0_fourVcor[1].T();
            if(firstEvent){
              sprintf(label,"kin_pim_E");
              vnames->Add(label);
              //cout<<"pim_E:\t\t PiMinus Energy \n";
            }

	    values[n_vectors++] = ~(pim);
	    if(firstEvent){
	      sprintf(label,"pim_m");
	      vnames->Add(label);
	      //cout<<"pim_m:\t\t PiMinus Mass \n";
	    }
	    values[n_vectors++] = sqrt( pim.x()*pim.x() + pim.y()*pim.y() );
	    if(firstEvent){
	      sprintf(label,"pim_pt");
	      vnames->Add(label);
	      //cout<<"pim_pt:\t\t PiMinus transverse momentum \n";
	    }
	    values[n_vectors++] = pim.V().theta();
	    if(firstEvent){
	      sprintf(label,"pim_theta");
	      vnames->Add(label);
	      //cout<<"pim_theta:\t\t PiMinus lab theta \n";
	    }
	    values[n_vectors++] = pim.V().phi();
	    if(firstEvent){
	      sprintf(label,"pim_phi");
	      vnames->Add(label);
	      //cout<<"pim_phi:\t\t PiMinus lab phi \n";
	    }
	    values[n_vectors++] = mm.x();
	    if(firstEvent){
	      sprintf(label,"mm_px");
	      vnames->Add(label);
	      //cout<<"mm_px:\t\t Missing x momentum \n";
	    }
	    values[n_vectors++] = mm.y();
	    if(firstEvent){
	      sprintf(label,"mm_py");
	      vnames->Add(label);
	      //cout<<"mm_py:\t\t Missing y momentum \n";
	    }
	    values[n_vectors++] = mm.r();
	    if(firstEvent){
	      sprintf(label,"mm_p");
	      vnames->Add(label);
	      //cout<<"mm_p:\t\t Missing momentum \n";
	    }
	    values[n_vectors++] = mm.z();
	    if(firstEvent){
	      sprintf(label,"mm_pz");
	      vnames->Add(label);
	      //cout<<"mm_pz:\t\t mm z momentum \n";
	    }
	    values[n_vectors++] = mm.t();
	    if(firstEvent){
	      sprintf(label,"mm_E");
	      vnames->Add(label);
	      //cout<<"mm_E:\t\t Missing Energy \n";
	    }
	    values[n_vectors++] = ~(mm);
	    if(firstEvent){
	      sprintf(label,"mm_m");
	      vnames->Add(label);
	      //cout<<"mm_pz:\t\t Missing Mass \n";
	    }
	    values[n_vectors++] = sqrt( mm.x()*mm.x() + mm.y()*mm.y() );
	    if(firstEvent){
	      sprintf(label,"mm_pt");
	      vnames->Add(label);
	      //cout<<"mm_pt:\t\t mm transverse momentum \n";
	    }
	    values[n_vectors++] = mm.V().theta();
	    if(firstEvent){
	      sprintf(label,"mm_theta");
	      vnames->Add(label);
	      //cout<<"mm_theta:\t\t Mm lab theta \n";
	    }
	    values[n_vectors++] = mm.V().phi();
	    if(firstEvent){
	      sprintf(label,"mm_phi");
	      vnames->Add(label);
	      //cout<<"mm_phi:\t\t mm lab phi \n";
	    }
	    values[n_vectors++] = mm2.t();
	    if(firstEvent){
	      sprintf(label,"mm2_E");
	      vnames->Add(label);
	      //cout<<"mm2_E:\t\t Missing Energy \n";
	    }
	    values[n_vectors++] = ~(mm2);
	    if(firstEvent){
	      sprintf(label,"mm2_m");
	      vnames->Add(label);
	      //cout<<"mm2_pz:\t\t Missing Mass \n";
	    }
	    values[n_vectors++] = sqrt( mm2.x()*mm2.x() + mm2.y()*mm2.y() );
	    if(firstEvent){
	      sprintf(label,"mm2_pt");
	      vnames->Add(label);
	      //cout<<"mm2_pt:\t\t mm2 transverse momentum \n";
	    }
	    values[n_vectors++] = mm2.V().theta();
	    if(firstEvent){
	      sprintf(label,"mm2_theta");
	      vnames->Add(label);
	      //cout<<"mm2_theta:\t\t Mm lab theta \n";
	    }
	    values[n_vectors++] = mm2.V().phi();
	    if(firstEvent){
	      sprintf(label,"mm2_phi");
	      vnames->Add(label);
	      //cout<<"mm2_phi:\t\t mm2 lab phi \n";
	    }

	    values[n_vectors++] = pip.x();
	    if(firstEvent){
	      sprintf(label,"pip_px");
	      vnames->Add(label);
	      //cout<<"pip_px:\t\t PiPlus x momentum \n";
	    }
	    values[n_vectors++] = pip.y();
	    if(firstEvent){
	      sprintf(label,"pip_py");
	      vnames->Add(label);
	      //cout<<"pip_py:\t\t PiPlus y momentum \n";
	    }
	    values[n_vectors++] = pip.r();
	    if(firstEvent){
	      sprintf(label,"pip_p");
	      vnames->Add(label);
	      //cout<<"pip_p:\t\t PiPlus momentum \n";
	    }

	    values[n_vectors++] = pip.z();
	    if(firstEvent){
	      sprintf(label,"pip_pz");
	      vnames->Add(label);
	      //cout<<"pip_pz:\t\t PiPlus z momentum \n";
	    }
	    values[n_vectors++] = pip.t();
	    if(firstEvent){
	      sprintf(label,"pip_E");
	      vnames->Add(label);
	      //cout<<"pip_E:\t\t PiPlus Energy \n";
	    }
            values[n_vectors++] = bec_pi0_fourVcor[0].X();
            if(firstEvent){
              sprintf(label,"kin_pip_px");
              vnames->Add(label);
              //cout<<"pip_px:\t\t PiPlus x momentum \n";
            }
            values[n_vectors++] = bec_pi0_fourVcor[0].Y();
            if(firstEvent){
              sprintf(label,"kin_pip_py");
              vnames->Add(label);
              //cout<<"pip_py:\t\t PiPlus y momentum \n";
            }
            values[n_vectors++] = bec_pi0_fourVcor[0].Z();
            if(firstEvent){
              sprintf(label,"kin_pip_pz");
              vnames->Add(label);
              //cout<<"pip_pz:\t\t PiPlus z momentum \n";
            }
            values[n_vectors++] = bec_pi0_fourVcor[0].T();
            if(firstEvent){
              sprintf(label,"kin_pip_E");
              vnames->Add(label);
              //cout<<"pip_E:\t\t PiPlus Energy \n";
            }

	    values[n_vectors++] = ~(pip);
	    if(firstEvent){
	      sprintf(label,"pip_m");
	      vnames->Add(label);
	      //cout<<"pip_m:\t\t PiPlus Mass \n";
	    }
	    values[n_vectors++] = sqrt( pip.x()*pip.x() + pip.y()*pip.y() );
	    if(firstEvent){
	      sprintf(label,"pip_pt");
	      vnames->Add(label);
	      //cout<<"pip_pt:\t\t PiPlus transverse momentum \n";
	    }
	    values[n_vectors++] = pip.V().theta();
	    if(firstEvent){
	      sprintf(label,"pip_theta");
	      vnames->Add(label);
	      //cout<<"pip_theta:\t\t PiPlus lab theta \n";
	    }
	    values[n_vectors++] = pip.V().phi();
	    if(firstEvent){
	      sprintf(label,"pip_phi");
	      vnames->Add(label);
	      //cout<<"pip_phi:\t\t PiPlus1 lab phi \n";
	    }

	    // proton
	    values[n_vectors++] = proton.x();
	    if(firstEvent){
	      sprintf(label,"p_px");
	      vnames->Add(label);
	      //cout<<"p_px:\t\t Proton x momentum \n";
	    }
	    values[n_vectors++] = proton.y();
	    if(firstEvent){
	      sprintf(label,"p_py");
	      vnames->Add(label);
	      //cout<<"p_py:\t\t Proton y momentum \n";
	    }
	    values[n_vectors++] = proton.r();
	    if(firstEvent){
	      sprintf(label,"p_p");
	      vnames->Add(label);
	      //cout<<"p_p:\t\t Proton momentum \n";
	    }

	    values[n_vectors++] = proton.z();
	    if(firstEvent){
	      sprintf(label,"p_pz");
	      vnames->Add(label);
	      //cout<<"p_pz:\t\t Proton z momentum \n";
	    }
	    values[n_vectors++] = proton.t();
	    if(firstEvent){
	      sprintf(label,"p_E");
	      vnames->Add(label);
	      //cout<<"p_E:\t\t Proton Energy \n";
	    }
            values[n_vectors++] = bec_pi0_fourVcor[2].X();
            if(firstEvent){
              sprintf(label,"kin_p_px");
              vnames->Add(label);
              //cout<<"p_px:\t\t Proton x momentum \n";
            }
            values[n_vectors++] = bec_pi0_fourVcor[2].Y();
            if(firstEvent){
              sprintf(label,"kin_p_py");
              vnames->Add(label);
              //cout<<"p_py:\t\t Proton y momentum \n";
            }
            values[n_vectors++] = bec_pi0_fourVcor[2].Z();
            if(firstEvent){
              sprintf(label,"kin_p_pz");
              vnames->Add(label);
              //cout<<"p_pz:\t\t Proton z momentum \n";
            }
            values[n_vectors++] = bec_pi0_fourVcor[2].T();
            if(firstEvent){
              sprintf(label,"kin_p_E");
              vnames->Add(label);
              //cout<<"p_E:\t\t Proton Energy \n";
            }

	    values[n_vectors++] = ~(proton);
	    if(firstEvent){
	      sprintf(label,"p_m");
	      vnames->Add(label);
	      //cout<<"p_m:\t\t Proton Mass \n";
	    }
	    values[n_vectors++] = sqrt( proton.x()*proton.x() + proton.y()*proton.y() );
	    if(firstEvent){
	      sprintf(label,"proton_pt");
	      vnames->Add(label);
	      //cout<<"proton_pt:\t\t Proton transverse momentum \n";
	    }
	    values[n_vectors++] = proton.V().theta();
	    if(firstEvent){
	      sprintf(label,"proton_theta");
	      vnames->Add(label);
	      //cout<<"proton_theta:\t\t Proton lab theta \n";
	    }
	    values[n_vectors++] = proton.V().phi();
	    if(firstEvent){
	      sprintf(label,"proton_phi");
	      vnames->Add(label);
	      //cout<<"proton_phi:\t\t Proton lab phi \n";
	    }
	    

	    //
	    // particle betas
	    //
	    values[n_vectors++] = event.cp(PiMinus, 1).beta();
	    if(firstEvent){
	      sprintf(label,"betapim_tof");
	      vnames->Add(label);
	      //cout<<"betapim:\t\t  Beta via TOF\n";
	    }
	    values[n_vectors++] = event.cp(PiPlus, 1).beta();
	    if(firstEvent){
	      sprintf(label,"betapip_tof");
	      vnames->Add(label);
	      //cout<<"betapip:\t\t  Beta via TOF\n";
	    }
	    values[n_vectors++] = event.cp(Proton, 1).beta();
	    if(firstEvent){
	      sprintf(label,"betaproton_tof");
	      vnames->Add(label);
	      //cout<<"betaproton:\t\t  Beta via TOF\n";
	    }
	    values[n_vectors++] = event.cp(PiMinus, 1).Beta();
	    if(firstEvent){
	      sprintf(label,"betapim");
	      vnames->Add(label);
	      //cout<<"betapim:\t\t  Beta via p/E\n";
	    }
	    values[n_vectors++] = event.cp(PiPlus, 1).Beta();
	    if(firstEvent){
	      sprintf(label,"betapip");
	      vnames->Add(label);
	      //cout<<"betapip:\t\t  Beta via p/E\n";
	    }
	    values[n_vectors++] = event.cp(Proton, 1).Beta();
	    if(firstEvent){
	      sprintf(label,"betaproton");
	      vnames->Add(label);
	      //cout<<"betaproton:\t\t  Beta via p/E\n";
	    }

	    //
	    //  Book Invariant Masses
	    //

	    values[n_vectors++] = ~(mm + pip); 
	    if(firstEvent){
	      sprintf(label,"pip_mm");
	      vnames->Add(label);
	      //cout<<"mass2pi:\t\t  Mass(piopip)\n";
	    }
	    values[n_vectors++] = ~(mm + pim); 
	    if(firstEvent){
	      sprintf(label,"pim_mm");
	      vnames->Add(label);
	      //cout<<"mass2pi:\t\t  Mass(piopim)\n";
	    }
	    values[n_vectors++] = ~(pip + pim); 
	    if(firstEvent){
	      sprintf(label,"pip_pim");
	      vnames->Add(label);
	      //cout<<"mass2pi2:\t\t  Mass(pippim)\n";
	    }
	    values[n_vectors++] = ~(pip + pim + mm); 
	    if(firstEvent){
	      sprintf(label,"pip_pim_mm");
	      vnames->Add(label);
	      //cout<<"mass3pi:\t\t  Mass(3pi)\n";
	    }
	    


	    //Mass( p pi )
	    values[n_vectors++] = ~(proton + pim+mm); 
	    if(firstEvent){
	      sprintf(label,"p_pim_mm");
	      vnames->Add(label);
	    }
	    values[n_vectors++] = ~(proton + pim); 
	    if(firstEvent){
	      sprintf(label,"p_pim");
	      vnames->Add(label);
	    }
	    values[n_vectors++] = ~(proton + pip); 
	    if(firstEvent){
	      sprintf(label,"p_pip");
	      vnames->Add(label);
	    }
	    values[n_vectors++] = ~(proton + pip+mm); 
	    if(firstEvent){
	      sprintf(label,"p_pip_mm");
	      vnames->Add(label);
	    }
	    values[n_vectors++] = ~(proton + pip+pim); 
	    if(firstEvent){
	      sprintf(label,"p_pip_pim");
	      vnames->Add(label);
	    }
	    values[n_vectors++] = ~(proton + mm); 
	    if(firstEvent){
	      sprintf(label,"p_mm");
	      vnames->Add(label);
	    }
	    values[n_vectors++] = ~(proton + mm + pim + pip); 
	    if(firstEvent){
	      sprintf(label,"p_pip_pim_mm");
	      vnames->Add(label);
	    }

	    values[n_vectors++] = mm.lenSq();
	    if(firstEvent){
	      sprintf(label,"mmsqorig");
	      vnames->Add(label);
	      //cout<<"mm:\t\t  Missing mass squared\n";
	    }
	    values[n_vectors++] = sqrt(mm.lenSq());
	    if(firstEvent){
	      sprintf(label,"mmorig");
	      vnames->Add(label);
	      //cout<<"mm:\t\t  Missing mass 2\n";
	    }

	    values[n_vectors++] = (~(mm))*(~(mm));
	    if(firstEvent){
	      sprintf(label,"mmsq2");
	      vnames->Add(label);
	      //cout<<"mm:\t\t  Missing mass squared\n";
	    }
	    
	    //////////////
	    if(firstEvent){
	      if(vnames->GetNlabels() != n_vectors){
		cerr<<"Error! vnames->GetNlabels()  !=  n_vectors\nExiting\n"
		    <<vnames->GetNlabels()<<"="<<n_vectors<<endl;
		
		exit(-1);
	      }
	      // Now create the ntuple
	      Char_t *vectorNames = vnames->GetLables();
	      cerr<<"Creating Ntuple\n\tTitle: "
		  <<ntpTitle<<"\n\tVectorNames: "
		  <<vectorNames<<endl<<endl;
	      ntpfile->cd();// create the ntuple in this file
	      ntp = new TNtuple("ntp",ntpTitle,vectorNames);
	      firstEvent=kFALSE;
	      
	    }
	    ntp->Fill(values);
	    
	    bookCount++; 
	  }
	  
	  
	  if(EventCount++%100==0)
	    cerr<<"event no: "<<EventCount<<"\tBooked: " << bookCount <<"\r";
	  
	  // clear buffers
	  dropAllBanks(&bcs_,"E");
	  cleanBanks(&bcs_);
	  
	  
	} 
	
      }
      
      if(EventCount)
        cerr << endl << setw(30) << "Events processed from file: "
             << EventCount << endl << setw(30) << "Total: " << Nevents << endl << setw(30) << "Booked: "<< bookCount << endl << setw(30) << "Different beam: "<< diffbeam << endl << endl;
    }// end of processing bos files
  }
  cerr << "\nTotal number of events read:\t" << Nevents << endl;

  ntpfile->Write();
  ntpfile->Close();

  return (0);
  }