#include <cstdlib>
#include <vector>
#include <iostream>
#include <map>
#include <string>

#include "TFile.h"
#include "TTree.h"
#include "TString.h"
#include "TSystem.h"
#include "TROOT.h"
#include "TPluginManager.h"
#include "TStopwatch.h"
//
#include "TROOT.h"
#include "TRint.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TChain.h"
#include "TVector3.h"
#include "TMath.h"
#include "TFile.h"
#include "TSystem.h"
#include "TGStatusBar.h"
#include "TSystem.h"
#include "TXMLEngine.h"
#include "TTree.h"
#include "TLorentzVector.h"
#include "TNtuple.h"
#include "TStyle.h"
#include "TGraphErrors.h"
#include "TGraph.h"
#include "TLegend.h"
#include "TPaveStats.h"
#include "TList.h"
#include "TMatrixT.h"
#include "TMatrixD.h"
#include "TLatex.h"

#include <cmath>
#include <iostream>
#include <string>
#include <cstdlib>
#include<iomanip>
#include <fstream>
#include <map>
#include <sstream>
#include <stdexcept>
#include <vector>

//#include "TMVAGui.C"

#if not defined(__CINT__) || defined(__MAKECINT__)
#include "TMVA/Tools.h"
#include "TMVA/Reader.h"
#include "TMVA/MethodCuts.h"
#endif

using namespace std;
#define SQR(a) ((a)*(a))
#define PI 3.14159
ofstream qin;

TDirectoryFile *FitDir[500];
void dbeta_study2()
{   

TLorentzVector MyProton, MyPip, MyPim, MyProton_Rest, My_Initial, MyPi0, MyPhoton, MM, MyKaon, MySigma;
         TVector3 omega_pvec;

         Float_t m_p = 0.938;
Float_t mp = 0.93827201;
  float omega_p;
  float kaon_mass;
  float sigma_mass;

Float_t  p_tof, pim_tof, pip_tof, beamE;
Float_t  p_tofp, pim_tofp, pip_tofp;
Float_t p_sec, pip_sec, pim_sec;
Float_t bitcon; 
//Int_t trigbit;
Float_t trigbit;


  int prefix;
  //Int_t NRun = 0;
  Float_t event;
  
  Float_t kin_beamE;
  Float_t mm_m;

  Float_t pip_pim_mm;

  //things not in gamp2root

  Float_t pip_pim, run;
  Float_t p_pip_pim;

  Float_t pip_mm;
  Float_t pim_mm;
  Float_t p_pip;
  Float_t p_pim;
  Float_t p_mm;
  Float_t p_pip_mm;
  Float_t p_pim_mm;

  Float_t z;
  Float_t x;
  Float_t y;

  Float_t mmsq;

  Float_t betapim;
  Float_t betapim_tof;
  Float_t betapip,betaproton;
  Float_t betapip_tof,betaproton_tof;
  Float_t dbetapim;
  Float_t dbetapip, dbetap;

  Float_t vtd;

  Float_t p_ngrf;
  Float_t pip_ngrf;
  Float_t pim_ngrf;
  Float_t p_tagr_id;
  Float_t pip_tagr_id;
  Float_t pim_tagr_id;

  Float_t mm_E;


  Float_t pip_px;
  Float_t pip_py;
  Float_t pip_pz;
  Float_t pim_px;
  Float_t pim_py;
  Float_t pim_pz;
  Float_t kin_pip_px;
  Float_t kin_pip_py;
  Float_t kin_pip_pz;
  Float_t kin_pim_px;
  Float_t kin_pim_py;
  Float_t kin_pim_pz;
  Float_t mm_px;
  Float_t mm_py;
  Float_t mm_pz;
  Float_t p_px;
  Float_t p_py;
  Float_t p_pz;
  Float_t kin_p_px;
  Float_t kin_p_py;
  Float_t kin_p_pz;
  Float_t pip_E;
  Float_t pim_E;
  Float_t p_E;
  Float_t kin_pip_E;
  Float_t kin_pim_E;
  Float_t kin_p_E;
  
  /*Float_t pip_p;
  Float_t pim_p;
  Float_t p_p;*/

  Float_t theta_mesons;
  Float_t t;
  Float_t p_m;

  Float_t beta_cm;
  Float_t gamma_cm,gamma_cm2;
  Float_t x_pz,x_px,x_py;
  Float_t x_pz_cm,p0_E_cm,p0_pz_cm,beamE_cm,p_pz_cm,p_E_cm;
  Float_t x_E;
  Float_t x_pt,theta_cm,cos_theta_cm;
  Int_t pi0_fail,mm_E_fail,dbeta_fail,vtd_fail;

  Float_t bec_pi0_prob;
  Float_t bec_nop_prob;

Float_t p_stvtime;
Float_t pip_stvtime;
Float_t pim_stvtime;

Float_t p_scvtime;
Float_t pip_scvtime;
Float_t pim_scvtime;

Float_t p_sc_id_part;
Float_t pip_sc_id_part;
Float_t pim_sc_id_part;

Float_t d_pim1;
Float_t d_pim2;


char name1[500];

float phi_p =0;
//float phi_rad_p=0;


float phi_pip=0;
//float phi_rad_pip=0;


float phi_pim=0;
//float phi_rad_pim=0;


TLorentzVector MyProton, MyPip, MyPim, MyProton_Rest, My_Initial, MyPi0, MyPhoton, MM, MyDeltazero, MyDeltaplus;
         TVector3 omega_pvec;
         Float_t m_p = 0.938;
TFile f ("beta_study_MC.root", "RECREATE", "Histograms from ntuples" );   
TH1D *three_sec = new TH1D("three_sec","energy",86,1.1,5.4);
TH1D *two_sec = new TH1D("two_sec","energy",86,1.1,5.4);

TH2D *proton = new TH2D("beta_proton","beta_proton",100,0,1.0,100,0,0.2);
TH2D *pip = new TH2D("beta_pip","beta_pip",100,0,1.0,100,0,0.2);
TH2D *pim = new TH2D("beta_pim","beta_pim",100,0,1.0,100,0,0.2);

TH2D *time = new TH2D("time","time",40,-10,10,40,-10,10);
TH2D *dbeta= new TH2D("dbeta","dbeta",40,-0.2,0.2,40,-0.2,0.2);

TH1D *mass_good_1= new TH1D("mass_good_1","mass",120,0.4,1.0);
TH1D *mass_bad_1 = new TH1D("mass_bad_1","mass",120,0.4,1.0);


TH1D *eta_good_1 = new TH1D("eta_good_1","mass",40,0.45,0.65);
TH1D *eta_bad_1 = new TH1D("eta_bad_1","mass",40,0.45,0.65);

TH1D *kaon_good_1 = new TH1D("kaon_good_1","mass",40,0.4,0.6);
TH1D *kaon_bad_1 = new TH1D("kaon_bad_1","mass",40,0.4,0.6);

TH1D *kaon_good_low = new TH1D("kaon_good_low","mass",40,0.4,0.6);
TH1D *kaon_good_high = new TH1D("kaon_good_high","mass",40,0.4,0.6);

TH1D *kaon_bad_low = new TH1D("kaon_bad_low","mass",40,0.4,0.6);
TH1D *kaon_bad_high = new TH1D("kaon_bad_high","mass",40,0.4,0.6);

TH1D *omega_good_low = new TH1D("omega_good_low","mass",40,0.4,0.6);
TH1D *omega_good_high = new TH1D("omega_good_high","mass",40,0.4,0.6);

TH1D *omega_bad_low = new TH1D("omega_bad_low","mass",40,0.4,0.6);
TH1D *omega_bad_high = new TH1D("omega_bad_high","mass",40,0.4,0.6);

TH1D *omega_good_1 = new TH1D("omega_good_1","mass",40,0.65,0.85);
TH1D *omega_bad_1 = new TH1D("omega_bad_1","mass",40,0.65,0.85);

TH1D *deltazero_good_1 = new TH1D("deltazero_good_1","mass",80,1.0,1.4);
TH1D *deltazero_bad_1 = new TH1D("deltazero_bad_1","mass",80,1.0,1.4);

TH1D *deltaplus_good_1 = new TH1D("deltaplus_good_1","mass",80,1.0,1.4);
TH1D *deltaplus_bad_1 = new TH1D("deltaplus_bad_1","mass",80,1.0,1.4);
//
TH1D *mass_good_2= new TH1D("mass_good_2","mass",60,0.4,1.0);
TH1D *mass_bad_2 = new TH1D("mass_bad_2","mass",60,0.4,1.0);


TH1D *eta_good_2 = new TH1D("eta_good_2","mass",40,0.45,0.65);
TH1D *eta_bad_2 = new TH1D("eta_bad_2","mass",40,0.45,0.65);

TH1D *omega_good_2 = new TH1D("omega_good_2","mass",40,0.65,0.85);
TH1D *omega_bad_2 = new TH1D("omega_bad_2","mass",40,0.65,0.85);
//
TH1D *mass_good_3= new TH1D("mass_good_3","mass",60,0.4,1.0);
TH1D *mass_bad_3 = new TH1D("mass_bad_3","mass",60,0.4,1.0);


TH1D *eta_good_3 = new TH1D("eta_good_3","mass",40,0.45,0.65);
TH1D *eta_bad_3 = new TH1D("eta_bad_3","mass",40,0.45,0.65);

TH1D *omega_good_3 = new TH1D("omega_good_3","mass",40,0.65,0.85);
TH1D *omega_bad_3 = new TH1D("omega_bad_3","mass",40,0.65,0.85);

TH1D *kaon_good_3 = new TH1D("kaon_good_3","mass",40,0.4,0.6);
TH1D *kaon_bad_3 = new TH1D("kaon_bad_3","mass",40,0.4,0.6);

//

TH1D *mass_good_4= new TH1D("mass_good_4","mass",60,0.4,1.0);
TH1D *mass_bad_4 = new TH1D("mass_bad_4","mass",60,0.4,1.0);

TH1D *eta_good_4 = new TH1D("eta_good_4","mass",40,0.45,0.65);
TH1D *eta_bad_4 = new TH1D("eta_bad_4","mass",40,0.45,0.65);

TH1D *kaon_good_4 = new TH1D("kaon_good_4","mass",40,0.4,0.6);
TH1D *kaon_bad_4 = new TH1D("kaon_bad_4","mass",40,0.4,0.6);

TH1D *sigma_good_4 = new TH1D("sigma_good_4","mass",40,1.1,1.3);
TH1D *sigma_bad_4 = new TH1D("sigma_bad_4","mass",40,1.1,1.3);

//
TH1D *kaon_good_5 = new TH1D("kaon_good_5","mass",40,0.4,0.6);
TH1D *kaon_bad_5 = new TH1D("kaon_bad_5","mass",40,0.4,0.6);

TH1D *sigma_good_5 = new TH1D("sigma_good_5","mass",40,1.1,1.3);
TH1D *sigma_bad_5 = new TH1D("sigma_bad_5","mass",40,1.1,1.3);TH1D *eta_good_3 = new TH1D("eta_good_3","mass",40,0.45,0.65);

TH1D *eta_good_5 = new TH1D("eta_good_5","mass",40,0.45,0.65);
TH1D *eta_bad_5 = new TH1D("eta_bad_5","mass",40,0.45,0.65);

TH1D *omega_good_5 = new TH1D("omega_good_5","mass",40,0.65,0.85);
TH1D *omega_bad_5 = new TH1D("omega_bad_5","mass",40,0.65,0.85);



TChain *ttree;
ttree = new TChain("ntp");
ttree->Add("/d/grid13/romines/3pion/sim/3pionsim*.root");
ttree->Add("/d/grid6/akbar/clas/eta/sim_1.3_1.5.root"); //eta mc
ttree->Add("/d/grid6/akbar/clas/eta/sim_1.5_1.7.root"); 
ttree->Add("/d/grid6/akbar/clas/eta/sim_1.7_1.9.root"); 
ttree->Add("/d/grid6/akbar/clas/eta/sim_1.9_2.2.root"); 
ttree->Add("/d/grid6/akbar/clas/eta/sim_2.2_2.8.root"); 
ttree->Add("/d/grid6/akbar/clas/eta/sim_2.8_3.8.root");
//ttree->Add("/d/grid6/akbar/g12/ROOT-new/56525.*.ppm.bos.root");
//ttree->Add("/d/grid11/akbar/g12/ROOT-new/56525.*.ppm.bos.root"); //this is 2pion
//ttree->Add("/d/grid13/akbar/new-files-may/ROOT-new/56525.*.ppm.bos.root");
//ttree->Add("/d/grid13/akbar/new-files-may/ROOT-new/5652*.*.ppm.bos.root");
//ttree->Add("/d/grid13/akbar/new-files-may/ROOT-new/5653*.*.ppm.bos.root");

//ttree->Add("/d/grid6/akbar/g12/ROOT-new/56525.*.ppm.bos.root");
//ttree->Add("/d/grid6/akbar/g12/ROOT-new/56527.*.ppm.bos.root");
//ttree->Add("/d/grid6/akbar/g12/ROOT-new/5652*.*.ppm.bos.root");
//ttree->Add("/d/grid6/akbar/g12/ROOT-new/5653*.*.ppm.bos.root");
//ttree->Add("/d/grid6/akbar/g12/ROOT-new/5654*.*.ppm.bos.root");
/*ttree->Add("/d/grid6/akbar/g12/ROOT-new/5655*.*.ppm.bos.root");
ttree->Add("/d/grid6/akbar/g12/ROOT-new/5656*.*.ppm.bos.root");

ttree->Add("/d/grid6/akbar/g12/ROOT-new/5657*.*.ppm.bos.root");
ttree->Add("/d/grid6/akbar/g12/ROOT-new/5658*.*.ppm.bos.root");
ttree->Add("/d/grid6/akbar/g12/ROOT-new/5659*.*.ppm.bos.root");*/



    


TTree *tree = (TTree*)gROOT->FindObject("ntp");

tree->SetBranchAddress( "p_sec", &p_sec);
tree->SetBranchAddress( "pip_sec", &pip_sec);
tree->SetBranchAddress( "pim_sec", &pim_sec);

tree->SetBranchAddress( "p_tof", &p_tof);
tree->SetBranchAddress( "pim_tof", &pim_tof);
tree->SetBranchAddress( "pip_tof", &pip_tof);

tree->SetBranchAddress( "p_tofp", &p_tofp);
tree->SetBranchAddress( "pim_tofp", &pim_tofp);
tree->SetBranchAddress( "pip_tofp", &pip_tofp);

tree->SetBranchAddress( "beamE", &beamE);
//tree->SetBranchAddress( "trigbit", &trigbit);
tree->SetBranchAddress("mm_m",&mm_m);
  tree->SetBranchAddress("mm_E",&mm_E);
  tree->SetBranchAddress("beamE",&beamE);


  tree->SetBranchAddress("pip_pim_mm",&pip_pim_mm);

  tree->SetBranchAddress("runNumber",&run);
  tree->SetBranchAddress("eventNumber",&event);
  tree->SetBranchAddress("pip_pim",&pip_pim);
  tree->SetBranchAddress("p_pip_pim",&p_pip_pim);

  tree->SetBranchAddress("pip_mm",&pip_mm);
  tree->SetBranchAddress("pim_mm",&pim_mm);
  tree->SetBranchAddress("p_pip_mm",&p_pip_mm);
  tree->SetBranchAddress("p_pim_mm",&p_pim_mm);
  tree->SetBranchAddress("p_pip",&p_pip);
  tree->SetBranchAddress("p_pim",&p_pim);
  tree->SetBranchAddress("p_mm",&p_mm);

  tree->SetBranchAddress("betapim",&betapim);
  tree->SetBranchAddress("betapim_tof",&betapim_tof);
  tree->SetBranchAddress("betapip",&betapip);
  tree->SetBranchAddress("betapip_tof",&betapip_tof);
  tree->SetBranchAddress("betaproton",&betaproton);
  tree->SetBranchAddress("betaproton_tof",&betaproton_tof);

  tree->SetBranchAddress("mmsqorig",&mmsq);


  tree->SetBranchAddress("z",&z);
  tree->SetBranchAddress("x",&x);
  tree->SetBranchAddress("y",&y);

  tree->SetBranchAddress("stVtimeDiff",&vtd);
  tree->SetBranchAddress("kin_pip_pz",&kin_pip_pz);
  tree->SetBranchAddress("kin_pip_px",&kin_pip_px);
  tree->SetBranchAddress("kin_pip_py",&kin_pip_py);
  tree->SetBranchAddress("kin_pim_pz",&kin_pim_pz);
  tree->SetBranchAddress("kin_pim_px",&kin_pim_px);
  tree->SetBranchAddress("kin_pim_py",&kin_pim_py);
  tree->SetBranchAddress("kin_p_pz",&kin_p_pz);
  tree->SetBranchAddress("kin_p_px",&kin_p_px);
  tree->SetBranchAddress("kin_p_py",&kin_p_py);
  tree->SetBranchAddress("kin_pip_E",&kin_pip_E);
  tree->SetBranchAddress("kin_pim_E",&kin_pim_E);
  tree->SetBranchAddress("kin_p_E",&kin_p_E);
  tree->SetBranchAddress("kin_beamE",&kin_beamE);
  tree->SetBranchAddress("p_ngrf",&p_ngrf);
  tree->SetBranchAddress("pip_ngrf",&pip_ngrf);
  tree->SetBranchAddress("pim_ngrf",&pim_ngrf);
  tree->SetBranchAddress("p_tagr_id",&p_tagr_id);
  tree->SetBranchAddress("pip_tagr_id",&pip_tagr_id);
  tree->SetBranchAddress("pim_tagr_id",&pim_tagr_id);

   tree->SetBranchAddress("p_m",&p_m);


  tree->SetBranchAddress("pip_pz",&pip_pz);
  tree->SetBranchAddress("pip_px",&pip_px);
  tree->SetBranchAddress("pip_py",&pip_py);
  tree->SetBranchAddress("pim_pz",&pim_pz);
  tree->SetBranchAddress("pim_px",&pim_px);
  tree->SetBranchAddress("pim_py",&pim_py);


  tree->SetBranchAddress("mm_pz",&mm_pz);
  tree->SetBranchAddress("mm_px",&mm_px);
  tree->SetBranchAddress("mm_py",&mm_py);
  tree->SetBranchAddress("p_pz",&p_pz);
  tree->SetBranchAddress("p_px",&p_px);
  tree->SetBranchAddress("p_py",&p_py);

  tree->SetBranchAddress("pip_E",&pip_E);
  tree->SetBranchAddress("pim_E",&pim_E);
  tree->SetBranchAddress("p_E",&p_E);

  
    tree->SetBranchAddress("bec_pi0_prob",&bec_pi0_prob);
    tree->SetBranchAddress("bec_nop_prob",&bec_nop_prob);


tree->SetBranchAddress("pim_stvtime",&pim_stvtime);
tree->SetBranchAddress("pip_stvtime",&pip_stvtime);
tree->SetBranchAddress("p_stvtime",&p_stvtime);

tree->SetBranchAddress("pim_scvtime",&pim_scvtime);
tree->SetBranchAddress("pip_scvtime",&pip_scvtime);
tree->SetBranchAddress("p_scvtime",&p_scvtime);

tree->SetBranchAddress("stVtimeDiff",&vtd);
tree->SetBranchAddress("p_sc_id_part",&p_sc_id_part);
tree->SetBranchAddress("pip_sc_id_part",&pip_sc_id_part);
tree->SetBranchAddress("pim_sc_id_part",&pim_sc_id_part);



//
Int_t nentries = (Int_t)tree->GetEntries();

Int_t counting;
counting = 0;



for (Int_t zz=0; zz<nentries ;zz++) {
tree->GetEntry(zz);
if (bec_pi0_prob > 0.01 && p_ngrf ==1 && pip_ngrf ==1 && pim_ngrf ==1)
{

counting = counting + 1;

cerr<<counting<<" "<<nentries<<endl;

dbetapim = betapim - betapim_tof;
      dbetapip = betapip - betapip_tof;
      dbetap = betaproton - betaproton_tof;

/*MyProton.SetPxPyPzE( p_px,p_py, p_pz, p_E );
                          MyPip.SetPxPyPzE( pip_px,pip_py, pip_pz, pip_E );
                          MyPim.SetPxPyPzE( pim_px,pim_py, pim_pz, pim_E);
                          MyProton_Rest.SetPxPyPzE( 0.0, 0.0, 0.0, m_p );
                          My_Initial.SetPxPyPzE(0.0,0.0,beamE, beamE+0.938);
                          MyPhoton.SetPxPyPzE(0.0,0.0,beamE, beamE);
                          MyPi0 = My_Initial - MyProton - MyPip - MyPim;
                          MM = MyPip + MyPim + MyPi0;
			  //omega_pvec = MM.Vect();
			  omega_p = MM.Mag();*/

MyProton.SetPxPyPzE( kin_p_px,kin_p_py, kin_p_pz, kin_p_E );
                          MyPip.SetPxPyPzE( kin_pip_px,kin_pip_py, kin_pip_pz, kin_pip_E );
                          MyPim.SetPxPyPzE( kin_pim_px,kin_pim_py, kin_pim_pz, kin_pim_E);
                          MyProton_Rest.SetPxPyPzE( 0.0, 0.0, 0.0, m_p );
                          My_Initial.SetPxPyPzE(0.0,0.0,kin_beamE, kin_beamE+0.938);
                          MyPhoton.SetPxPyPzE(0.0,0.0,kin_beamE, kin_beamE);
                          MyPi0 = My_Initial - MyProton - MyPip - MyPim;
                          MM = MyPip + MyPim + MyPi0;
			  //omega_pvec = MM.Vect();
			  omega_p = MM.Mag();

MyKaon = MyPip + MyPim;
kaon_mass = MyKaon.Mag();

MySigma = MyProton + MyPi0;
sigma_mass = MySigma.Mag();

double deltazero_mass = (MyProton + MyPim).Mag();
double deltaplus_mass = (MyProton + MyPip).Mag();

/*if (kin_beamE < 2.5)
{
deltazero_good_1->Fill(deltazero_mass);
deltaplus_good_1->Fill(deltaplus_mass);
}*/

if (fabs(dbetap) < 0.04 && fabs(dbetapip) < 0.04 && fabs(dbetapim) < 0.04)
//if (fabs(dbetap) < 0.04 || fabs(dbetapip) < 0.04 )
{
  mass_good_1->Fill(omega_p);
  omega_good_1->Fill(omega_p); 
  eta_good_1->Fill(omega_p);

  deltazero_good_1->Fill(deltazero_mass);
  deltaplus_good_1->Fill(deltaplus_mass);

    if (kin_beamE <= 1.55 )
    {
     omega_good_low->Fill(omega_p);
    }
    if (kin_beamE >= 1.55 )
    {
     kaon_good_high->Fill(omega_p);
    }

  if (sigma_mass >= 1.169 && sigma_mass <= 1.209)
  {
    kaon_good_1->Fill(kaon_mass);
    if (kin_beamE >= 1.6 && kin_beamE <= 1.8)
    {
     kaon_good_low->Fill(kaon_mass);
    }
    if (kin_beamE >= 2.6 && kin_beamE <= 2.8)
    {
     kaon_good_high->Fill(kaon_mass);
    }
  }
}
else
{
  mass_bad_1->Fill(omega_p);
  omega_bad_1->Fill(omega_p); 
  eta_bad_1->Fill(omega_p);

  if (kin_beamE <= 1.55 )
    {
     omega_bad_low->Fill(omega_p);
    }
    if (kin_beamE >= 1.55 )
    {
     kaon_bad_high->Fill(omega_p);
    }

deltazero_bad_1->Fill(deltazero_mass);
  deltaplus_bad_1->Fill(deltaplus_mass);  

  if (sigma_mass >= 1.169 && sigma_mass <= 1.209)
  {
    kaon_bad_1->Fill(kaon_mass);
    if (kin_beamE >= 1.6 && kin_beamE <= 1.8)
    {
     kaon_bad_low->Fill(kaon_mass);
    }
    if (kin_beamE >= 2.6 && kin_beamE <= 2.8)
    {
     kaon_bad_high->Fill(kaon_mass);
    }

  }
}

//if (fabs(dbetap) < 0.04 && fabs(dbetapip) < 0.04 && fabs(dbetapim) < 0.04)
if (fabs(dbetap) < 0.04 || fabs(dbetapip) < 0.04 )
{
  mass_good_2->Fill(omega_p);
  omega_good_2->Fill(omega_p); 
  eta_good_2->Fill(omega_p);
  }
else
{
  mass_bad_2->Fill(omega_p);
  omega_bad_2->Fill(omega_p); 
  eta_bad_2->Fill(omega_p);
}

if (fabs(p_tof-p_tofp) < 1.0 || fabs(pip_tof-pip_tofp) < 1.0 )
{



  mass_good_3->Fill(omega_p);
  omega_good_3->Fill(omega_p); 
  eta_good_3->Fill(omega_p);

  if (sigma_mass >= 1.169 && sigma_mass <= 1.209)
  {
    kaon_good_3->Fill(kaon_mass);
  }
}
else
{
  mass_bad_3->Fill(omega_p);
  omega_bad_3->Fill(omega_p); 
  eta_bad_3->Fill(omega_p);

  if (sigma_mass >= 1.169 && sigma_mass <= 1.209)
  {
    kaon_bad_3->Fill(kaon_mass);
  }
}

if (fabs(vtd) < 1.0 )
{


  mass_good_4->Fill(omega_p);
  eta_good_4->Fill(omega_p);

  if (sigma_mass >= 1.169 && sigma_mass <= 1.209)
  {
    kaon_good_4->Fill(kaon_mass);
  }

  if (kaon_mass >= 0.44 && kaon_mass <= 0.55)
  {
    sigma_good_4->Fill(sigma_mass);
  }
}
else
{
  mass_bad_4->Fill(omega_p);
  eta_bad_4->Fill(omega_p);

  if (sigma_mass >= 1.169 && sigma_mass <= 1.209)
  {
    kaon_bad_4->Fill(kaon_mass);
  }

  if (kaon_mass >= 0.44 && kaon_mass <= 0.55)
  {
    sigma_bad_4->Fill(sigma_mass);
  }

}

if ((x*x + y*y) < 4 )
{

omega_good_5->Fill(omega_p); 
eta_good_5->Fill(omega_p);
  

  if (sigma_mass >= 1.169 && sigma_mass <= 1.209)
  {
    kaon_good_5->Fill(kaon_mass);
  }

  if (kaon_mass >= 0.44 && kaon_mass <= 0.55)
  {
    sigma_good_5->Fill(sigma_mass);
  }
}
else
{
  
omega_bad_5->Fill(omega_p); 
eta_bad_5->Fill(omega_p);

  if (sigma_mass >= 1.169 && sigma_mass <= 1.209)
  {
    kaon_bad_5->Fill(kaon_mass);
  }

  if (kaon_mass >= 0.44 && kaon_mass <= 0.55)
  {
    sigma_bad_5->Fill(sigma_mass);
  }

}



 

//cerr<<pim_scvtime<<" "<<pim_stvtime<<" "<<pim_scvtime-pim_stvtime<<endl;

time->Fill(p_tof-p_tofp,pip_tof-pip_tofp);
dbeta->Fill(dbetap, dbetapip);
/*
if ( fabs(p_scvtime-p_stvtime) < 1 || fabs(pip_scvtime-pip_stvtime) < 1 )
{
 mass_p_good->Fill(omega_p);
}
else
{
 mass_p_bad->Fill(omega_p);
}*/

proton->Fill(omega_p, fabs(dbetap));
pip->Fill(omega_p, fabs(dbetapip));
pim->Fill(omega_p, fabs(dbetapim));



//}
//
}
}

proton->Write();
pip->Write();
pim->Write();


time->Write();
//dbeta->Write();

mass_good_1->Write();
mass_bad_1->Write();

omega_good_1->Write();
omega_bad_1->Write();

eta_good_1->Write();
eta_bad_1->Write();
//
deltazero_good_1->Write();
deltazero_bad_1->Write();

deltaplus_good_1->Write();
deltaplus_bad_1->Write();
//
mass_good_2->Write();
mass_bad_2->Write();

omega_good_2->Write();
omega_bad_2->Write();

eta_good_2->Write();
eta_bad_2->Write();

mass_good_3->Write();
mass_bad_3->Write();

omega_good_3->Write();
omega_bad_3->Write();

eta_good_3->Write();
eta_bad_3->Write();

kaon_good_3->Write();
kaon_bad_3->Write();

kaon_good_1->Write();
kaon_bad_1->Write();

kaon_good_high->Write();
kaon_bad_high->Write();

kaon_good_low->Write();
kaon_bad_low->Write();

omega_good_high->Write();
omega_bad_high->Write();

omega_good_low->Write();
omega_bad_low->Write();

mass_good_4->Write();
mass_bad_4->Write();


eta_good_4->Write();
eta_bad_4->Write();

kaon_good_4->Write();
kaon_bad_4->Write();

sigma_good_4->Write();
sigma_bad_4->Write();

kaon_good_5->Write();
kaon_bad_5->Write();

sigma_good_5->Write();
sigma_bad_5->Write();

omega_good_5->Write();
omega_bad_5->Write();

eta_good_5->Write();
eta_bad_5->Write();

cerr<<"FINISH FIRST LOOP";


}