#include <dirent.h>
#include <TLorentzVector.h>
#include <TLorentzRotation.h>
#include <TFile.h>
#include <TTree.h>
#include <TChain.h>
#include <TH1.h>
#include <TH2.h>
#include <TMinuit.h>
#include <TMath.h>

#include <vector>
#include<iostream>

using namespace std;

int EventNUM;
float RFT;
int Nbeam;
float BeamE[500];
float BeamT[500];
int Nneut;
float NeutralX[40];
float NeutralY[40];
float NeutralZ[40];
float NeutralT[40];
float NeutralE[40];
int   DetSys[40];
int Nstart;
float SC_E[60];
float SC_T[60];
int SC_S[60];


vector <int> Good_SC_idx;
vector <int> Good_gamma_idx;
vector <int> Good_beam_idx;

double ThePhiCut = 0.5;
double BeamEnergyCutCenter = 9.2;
double BeamEnergyCutWidth = 2.2;
double StartCounterTimePeak = 10.;
double MinimumPhotonEnergy = 0.15;


#define ProtonMass  .938272
TLorentzVector ProtonTarget(0.,0.,0., ProtonMass); // Proton Target 4-vector

int findphotons(vector <TLorentzVector> &, vector <double> &, double, vector <int> & );
double GetPhiAngle(TLorentzVector );
double GetThetaAngle(TLorentzVector );
int FindBestBeamPhotonMatch(vector <TLorentzVector> &, TLorentzVector &, TLorentzVector &, TH1D *, TH1D *, TH1D *, TH2D *);
TVector3 NormalizeVector(TVector3);

double FCAL_R_CUT = 65.;

void pi0anal3(int DataPeriod, int OR){

  if (OR<0 || OR>4){
    cout<<"OR out of range: 0==135DEG, 1==45DEG, 2==PARA, 3==PERP, 4==AMO"<<endl;
    return;
  }


  TChain * myChain(0);
  myChain = new TChain("TreeDir/t3");
  /*
  TString datatype(whichone);
  TString directory = "./spring2017/";
  directory += datatype;
  directory += "/";
  */
  
  TString directories[5] = {"./spring2017/135DEG/", "./spring2017/45DEG/","./spring2017/PARA/","./spring2017/PERP/",
			    "./spring2017/AMO/"};

  

  char outputrootfile[128];
  int Ndirs = 4;
  if (DataPeriod==1){
    sprintf(outputrootfile,"newroot_2017_lowint_%d.root",OR); // change >30795 to <30796 below
  } else if (DataPeriod==2){
    sprintf(outputrootfile,"newroot_2017_higint_%d.root",OR); // change <30796 to >30795 below
  } else if (DataPeriod==3){
    sprintf(outputrootfile,"newroot_2016_golden_%d.root",OR); 
    directories[0] = "./spring2016/PARA/";
    directories[1] = "./spring2016/PERP/";
    directories[2] = "dummy";
    directories[3] = "dummy";
    Ndirs = 2;
    StartCounterTimePeak = -40;
  }

  DIR *dir;
  struct dirent *ent;

  vector <TString> FileList;
  for (int id=OR; id<OR+1; id++) {
    if ((dir = opendir (directories[id])) != NULL) {
      /* print all the files and directories within directory */
      while ((ent = readdir (dir)) != NULL) {
	//printf ("%s\n", ent->d_name);
	TString s(ent->d_name);
	if (s.EndsWith(".root")){
	  TString r = s(0,6);
	  int rnum = atoi(r.Data());
	  cout<<rnum<<endl;
	  if (DataPeriod==1) {
	    if (rnum<30796){ // low intensity spring 2017 data  CHANGE IF NEEDED!!!!!
	      FileList.push_back(directories[id]+s);
	    }
	  } else if (DataPeriod==2){
	    if (rnum>30795){ // low intensity spring 2017 data  CHANGE IF NEEDED!!!!!
	      FileList.push_back(directories[id]+s);
	    }
	  } else if (DataPeriod==3){
	    if ((rnum>11365) && (rnum<11556)){ // low intensity spring 2017 data  CHANGE IF NEEDED!!!!!
	      FileList.push_back(directories[id]+s);
	    }
	  }
	}
      }
      closedir (dir);
    } else {
      /* could not open directory */
      cout<<"error open directory"<<endl;
      return;      
    }
  }
  std::sort(FileList.begin(), FileList.end());

  int nfiles = FileList.size();
  for (int k=0; k<nfiles; k++){
    TString newf = FileList[k];
    cout<<newf<<endl;
    myChain->Add(newf);
  }

  myChain->SetBranchAddress("EventNUM",&EventNUM);
  myChain->SetBranchAddress("RFT",&RFT);
  myChain->SetBranchAddress("Nbeam",&Nbeam);
  myChain->SetBranchAddress("BeamE",BeamE);
  myChain->SetBranchAddress("BeamT",BeamT);
  myChain->SetBranchAddress("Nneut",&Nneut);
  myChain->SetBranchAddress("NeutralX",NeutralX);
  myChain->SetBranchAddress("NeutralY",NeutralY);
  myChain->SetBranchAddress("NeutralZ",NeutralZ);
  myChain->SetBranchAddress("NeutralT",NeutralT);
  myChain->SetBranchAddress("NeutralE",NeutralE);
  myChain->SetBranchAddress("DetSys", DetSys);
  myChain->SetBranchAddress("Nstart",&Nstart);
  myChain->SetBranchAddress("SC_E",SC_E);
  myChain->SetBranchAddress("SC_T",SC_T);
  myChain->SetBranchAddress("SC_S",SC_S);

  cout<<"create histograms.......";

  TH1D *beamphotontime = new TH1D("beamphotontime", "Beam photon time minus RF-time", 1000, -25., 25.);
  TH1D *beamphotons[2];
  beamphotons[0] = new TH1D("beamphotons0", "Beam photon Energy out of time", 100, 0., 12.);
  beamphotons[1] = new TH1D("beamphotons1", "Beam photon Energy in of time", 100, 0., 12.);
  TH1D *NumberOfBeamPhotons[2];
  NumberOfBeamPhotons[1] = new TH1D("NumberOfBeamPhotons1", "Number of Beam photon in of time", 15, 0., 15.);
  NumberOfBeamPhotons[0] = new TH1D("NumberOfBeamPhotons0", "Number of Beam photon out of time", 15, 0., 15.);
  TH2D *StartCounterTime = new TH2D("StartCounterTime", "StartCounterNumber vs Start Counter Time", 3000, -150., 150., 32, 0., 32.);
  TH2D *StartCounterEnergy = new TH2D("StartCounterEnergy", "StartCounteRNumber vs Start Counter Energy", 150, 0., 0.015, 32, 0., 32.);
  TH2D *StarCounterEnergyWithTimingCut = new TH2D("StarCounterEnergyWithTimingCut", "StartCounterNumber vs Start Counter Energy",  150, 0., 0.015, 32, 0., 32.);
  TH2D *StartCountertimeWithEnergyCut = new TH2D("StartCountertimeWithEnergyCut", "StartCounterNumber vs Start Counter Time", 3000, -150., 150., 32, 0., 32.);
  
  TH1D *gammatimediff = new TH1D("gammatimediff", "2Gamma Times difference at target", 200, -10., 10.);
  TH1D *TwoPhotonOnlyInvMass = new TH1D("TwoPhotonOnlyInvMass","Invariant Mass of two photon system", 900, 0., 3.5);
  TH1D *TwoPhotonInvMass = new TH1D("TwoPhotonInvMass","Invariant Mass of two photon system", 900, 0., 3.5);
  TH2D *StartCounterNumVsInvMass = new TH2D("StartCounterNumVsInvMass", "SC# vs Inv Mass of 2-photon system", 900, 0., 3.5, 32, 0., 32);
  TH2D *StartCounterTimeVsInvMass = new TH2D("StartCounterTimeVsInvMass", "SC Time vs Inv Mass of 2-photon system", 900, 0., 3.5, 900, -150., 150.);
  TH2D *StartCounterEnergyVsInvMass = new TH2D("StartCounterEnergyVsInvMass", "SC Energy vs Inv Mass of 2-photon system", 900, 0., 3.5, 150, 0., 0.015);
  TH1D *StarCounterRelativeTime = new TH1D("StarCounterRelativeTime", "Mean Photon time minus SC time", 500, -50., 50.);
  TH2D *StartCounterTimeVsMeanGammaTime = new TH2D("StartCounterTimeVsMeanGammaTime","SC time vs mean gamma time", 2000, -100., 100., 3000, -150., 150.);
  TH2D *StartCounterEnergyVsInvMassAndTimeCut = new TH2D("StartCounterEnergyVsInvMassAndTimeCut", "SC Energy Vs Inv Mass with Time cut ", 800, 0., 4., 150, 0., 0.015);
  TH2D *StartCounterNumberVsStartCounterTimeAndTimeCut = new TH2D("StartCounterNumberVsStartCounterTimeAndTimeCut", "SC# vs SC time with Time Cut", 3000, -150., 150., 32, 0., 32);
  TH2D *StartCounterTimeVsInvMassAndTimeCut = new TH2D("StartCounterTimeVsInvMassAndTimeCut", "SC Time vs Inv. Mass with Time cut",  800, 0., 4., 3000, -150., 150.);
  TH1D *GoodStartCounterHitsAndTimeCut = new TH1D("GoodStartCounterHitsAndTimeCut", "Number of SC Hits with good time", 10, 0., 10.) ;

  TH1D *TwoPhotonOnlyInvMassAndTimeCut = new TH1D("TwoPhotonOnlyInvMassAndTimeCut", "Invariant 2-photon Mass (only2) with time cut",  800, 0., 4.);
  TH1D *TwoPhotonInvMassAndTimeCut = new TH1D("TwoPhotonInvMassAndTimeCut", "Invariant 2-photon Mass with time cut",  800, 0., 4.);
  TH1D *TwoPhotonOnlyPhi = new TH1D("TwoPhotonOnlyPhi", "Phi of 2-photon system (2only)",36,  0., 2.*TMath::Pi() );
  TH1D *TwoPhotonPhi = new TH1D("TwoPhotonPhi", "Phi of 2-photon system",36, 0., 2.*TMath::Pi() );
  TH1D *TwoPhotonOnlyTheta = new TH1D("TwoPhotonOnlyTheta", "Theta of 2-photon system (2-only)", 400, 0., 0.5);    
  TH1D *TwoPhotonTheta = new TH1D("TwoPhotonTheta", "Theta of 2-photon system", 400, -TMath::Pi(), TMath::Pi());    
  TH1D *MissingMass[3];
  MissingMass[0] = new TH1D("MissingMass0",  "Missing Mass with 2photons, Background", 100, -4., 8.);
  MissingMass[1] = new TH1D("MissingMass1",  "Missing Mass with 2photons, Signal Only", 100, -4., 8.);
  MissingMass[2] = new TH1D("MissingMass2",  "Missing Mass with 2photons, Final Bgr. subtracted", 100, -4., 8.);
  TH1D *MissingMomentum[3];
  MissingMomentum[0] = new TH1D("MissingMomentum0", "Missing Momentum with 2photons, Background", 100, -4., 8.);
  MissingMomentum[1] = new TH1D("MissingMomentum1", "Missing Momentum with 2photons, Signal", 100, -4., 8.);
  MissingMomentum[2] = new TH1D("MissingMomentum2", "Missing Momentum with 2photons, Final Bgr. subtracted", 100, -4., 8.);
  TH1D *MissingEnergy[3];
  MissingEnergy[0] = new TH1D("MissingEnergy0", "Missing Energy with 2photons, Background", 100, -4., 8. );
  MissingEnergy[1] = new TH1D("MissingEnergy1", "Missing Energy with 2photons, Signal", 100, -4., 8. );
  MissingEnergy[2] = new TH1D("MissingEnergy2", "Missing Energy with 2photons, Final Bgr. subtracted", 100, -4., 8. );
  TH2D *MomentumVsMissingMass[3];
  MomentumVsMissingMass[0] = new TH2D("MomentumVsMissingMass0", "Momentum Vs Missing Mass with 2 photongs", 100, -4., 8. , 100, -4., 8. );
  MomentumVsMissingMass[1] = new TH2D("MomentumVsMissingMass1", "Momentum Vs Missing Mass with 2 photongs", 100, -4., 8. , 100, -4., 8. );
  MomentumVsMissingMass[2] = new TH2D("MomentumVsMissingMass2", "Momentum Vs Missing Mass with 2 photongs", 100, -4., 8. , 100, -4., 8. );
  TH1D *MissingPhi[3];
  MissingPhi[0] = new TH1D("MissingPhi0","Missing Phi", 36, 0., 2.*TMath::Pi());
  MissingPhi[1] = new TH1D("MissingPhi1","Missing Phi", 36, 0., 2.*TMath::Pi());
  MissingPhi[2] = new TH1D("MissingPhi2","Missing Phi", 36, 0., 2.*TMath::Pi());
  TH1D *DeltaPhi[3];
  DeltaPhi[0] = new TH1D("DeltaPhi0", "phi 2photon system minus phi SC", 72, -TMath::Pi(), TMath::Pi());
  DeltaPhi[1] = new TH1D("DeltaPhi1", "phi 2photon system minus phi SC", 72, -TMath::Pi(), TMath::Pi());
  DeltaPhi[2] = new TH1D("DeltaPhi2", "phi 2photon system minus phi SC", 72, -TMath::Pi(), TMath::Pi());

  TH1D *TwoPhotonOnlyInvMassAndTimeCutAndPhiCut[3];
  TwoPhotonOnlyInvMassAndTimeCutAndPhiCut[0] = new TH1D("TwoPhotonOnlyInvMassAndTimeCutAndPhiCut0", "Inv. Mass 2gonly with Cuts in Time and phi", 800, 0., 4.);
  TwoPhotonOnlyInvMassAndTimeCutAndPhiCut[1] = new TH1D("TwoPhotonOnlyInvMassAndTimeCutAndPhiCut1", "Inv. Mass 2gonly with Cuts in Time and phi", 800, 0., 4.);
  TwoPhotonOnlyInvMassAndTimeCutAndPhiCut[2] = new TH1D("TwoPhotonOnlyInvMassAndTimeCutAndPhiCut2", "Inv. Mass 2gonly with Cuts in Time and phi", 800, 0., 4.);
  TH1D *TwoPhotonOnlyMissingMassAndTimeCutAndPhiCut[3];
  TwoPhotonOnlyMissingMassAndTimeCutAndPhiCut[0] = new TH1D("TwoPhotonOnlyMissingMassAndTimeCutAndPhiCut0", "Missing Mass 2gonly with Cuts in Time and phi", 100, -4., 8. );
  TwoPhotonOnlyMissingMassAndTimeCutAndPhiCut[1] = new TH1D("TwoPhotonOnlyMissingMassAndTimeCutAndPhiCut1", "Missing Mass 2gonly with Cuts in Time and phi", 100, -4., 8. );
  TwoPhotonOnlyMissingMassAndTimeCutAndPhiCut[2] = new TH1D("TwoPhotonOnlyMissingMassAndTimeCutAndPhiCut2", "Missing Mass 2gonly with Cuts in Time and phi", 100, -4., 8. );
  TH1D *TwoPhotonOnlyInvMassAndTimeCutAndPhiCutAndMMCut[3];
  TwoPhotonOnlyInvMassAndTimeCutAndPhiCutAndMMCut[0] = new TH1D("TwoPhotonOnlyInvMassAndTimeCutAndPhiCutAndMMCut0", "Inv. Mass 2gonly with Cuts in Time, phi and MM", 800, 0., 4.);
  TwoPhotonOnlyInvMassAndTimeCutAndPhiCutAndMMCut[1] = new TH1D("TwoPhotonOnlyInvMassAndTimeCutAndPhiCutAndMMCut1", "Inv. Mass 2gonly with Cuts in Time, phi and MM", 800, 0., 4.);
  TwoPhotonOnlyInvMassAndTimeCutAndPhiCutAndMMCut[2] = new TH1D("TwoPhotonOnlyInvMassAndTimeCutAndPhiCutAndMMCut2", "Inv. Mass 2gonly with Cuts in Time, phi and MM", 800, 0., 4.);
  TH2D *TwoPhotonOnlyPHIvsTwoPhotonInvMassAndTimeCutAndPhiCutAndMMCut[3];
  TwoPhotonOnlyPHIvsTwoPhotonInvMassAndTimeCutAndPhiCutAndMMCut[0] = new TH2D("TwoPhotonOnlyPHIvsTwoPhotonInvMassAndTimeCutAndPhiCutAndMMCut0", "2g only Phi vs Inv. Mass w.all cuts",
									      100, 0., 4., 36, 0., 2.*TMath::Pi() );
  TwoPhotonOnlyPHIvsTwoPhotonInvMassAndTimeCutAndPhiCutAndMMCut[1] = new TH2D("TwoPhotonOnlyPHIvsTwoPhotonInvMassAndTimeCutAndPhiCutAndMMCut1", "2g only Phi vs Inv. Mass w.all cuts",
									      100, 0., 4., 36, 0., 2.*TMath::Pi() );
  TwoPhotonOnlyPHIvsTwoPhotonInvMassAndTimeCutAndPhiCutAndMMCut[2] = new TH2D("TwoPhotonOnlyPHIvsTwoPhotonInvMassAndTimeCutAndPhiCutAndMMCut2", "2g only Phi vs Inv. Mass w.all cuts",
									      100, 0., 4., 36, 0., 2.*TMath::Pi());

  TH1D *TwoPhotonOnlyPi0PhiAndAllCuts[3];
  TwoPhotonOnlyPi0PhiAndAllCuts[0] = new TH1D("TwoPhotonOnlyPi0PhiAndAllCuts0", "Pi0 phi distribution Background", 16, 0., 2.*TMath::Pi());
  TwoPhotonOnlyPi0PhiAndAllCuts[1] = new TH1D("TwoPhotonOnlyPi0PhiAndAllCuts1", "Pi0 phi distribution Signal", 16, 0., 2.*TMath::Pi());
  TwoPhotonOnlyPi0PhiAndAllCuts[2] = new TH1D("TwoPhotonOnlyPi0PhiAndAllCuts2", "Pi0 phi distribution Signal Background subtracted", 16, 0., 2.*TMath::Pi());

  TH1D *TwoPhotonOnlyPi0ThetaAndAllCuts[3];
  TwoPhotonOnlyPi0ThetaAndAllCuts[0] = new TH1D("TwoPhotonOnlyPi0ThetaAndAllCuts0", "Pi0 Theta after all cuts Background", 100, 0., 0.1); 
  TwoPhotonOnlyPi0ThetaAndAllCuts[1] = new TH1D("TwoPhotonOnlyPi0ThetaAndAllCuts1", "Pi0 Theta after all cuts Signal", 100, 0., 0.1); 
  TwoPhotonOnlyPi0ThetaAndAllCuts[2] = new TH1D("TwoPhotonOnlyPi0ThetaAndAllCuts2", "Pi0 Theta after all cuts Signal Background subtracted", 100, 0., 0.1); 
  
  TH1D *TwoPhotonOnlyPi0_t_AndAllCuts[3];
  TwoPhotonOnlyPi0_t_AndAllCuts[0] = new TH1D("TwoPhotonOnlyPi0_t_AndAllCuts0", "Exclusive Pi0 -t  after all cuts Background", 50, 0., 0.2);
  TwoPhotonOnlyPi0_t_AndAllCuts[1] = new TH1D("TwoPhotonOnlyPi0_t_AndAllCuts1", "Exclusive Pi0 -t  after all cuts Signal", 50, 0., 0.2);
  TwoPhotonOnlyPi0_t_AndAllCuts[2] = new TH1D("TwoPhotonOnlyPi0_t_AndAllCuts2", "Exclusive Pi0 -t  after all cuts Signal Background subtracted", 50, 0., 0.2);

  TH1D *TwoPhotonPi0PhiAndAllCuts[3];
  TwoPhotonPi0PhiAndAllCuts[0] = new TH1D("TwoPhotonPi0PhiAndAllCuts0", "Pi0 phi distribution Background", 16, 0., 2.*TMath::Pi());
  TwoPhotonPi0PhiAndAllCuts[1] = new TH1D("TwoPhotonPi0PhiAndAllCuts1", "Pi0 phi distribution Signal", 16, 0., 2.*TMath::Pi());
  TwoPhotonPi0PhiAndAllCuts[2] = new TH1D("TwoPhotonPi0PhiAndAllCuts2", "Pi0 phi distribution Signal Background subtracted", 16, 0., 2.*TMath::Pi());


  TH1D *TwoPhotonInvMassAndTimeCutAndPhiCut[3];
  TwoPhotonInvMassAndTimeCutAndPhiCut[0] = new TH1D("TwoPhotonInvMassAndTimeCutAndPhiCut0", "Inv. Mass 2g with Cuts in Time and phi", 800, 0., 4.);
  TwoPhotonInvMassAndTimeCutAndPhiCut[1] = new TH1D("TwoPhotonInvMassAndTimeCutAndPhiCut1", "Inv. Mass 2g with Cuts in Time and phi", 800, 0., 4.);
  TwoPhotonInvMassAndTimeCutAndPhiCut[2] = new TH1D("TwoPhotonInvMassAndTimeCutAndPhiCut2", "Inv. Mass 2g with Cuts in Time and phi", 800, 0., 4.);
  TH1D *TwoPhotonMissingMassAndTimeCutAndPhiCut[3];
  TwoPhotonMissingMassAndTimeCutAndPhiCut[0] = new TH1D("TwoPhotonMissingMassAndTimeCutAndPhiCut0", "Missing Mass 2g with Cuts in Time and phi", 100, -4., 8. );
  TwoPhotonMissingMassAndTimeCutAndPhiCut[1] = new TH1D("TwoPhotonMissingMassAndTimeCutAndPhiCut1", "Missing Mass 2g with Cuts in Time and phi", 100, -4., 8. );
  TwoPhotonMissingMassAndTimeCutAndPhiCut[2] = new TH1D("TwoPhotonMissingMassAndTimeCutAndPhiCut2", "Missing Mass 2g with Cuts in Time and phi", 100, -4., 8. );
  TH1D *TwoPhotonInvMassAndTimeCutAndPhiCutAndMMCut[3];
  TwoPhotonInvMassAndTimeCutAndPhiCutAndMMCut[0] = new TH1D("TwoPhotonInvMassAndTimeCutAndPhiCutAndMMCut0", "Inv. Mass 2g with Cuts in Time, phi and MM", 800, 0., 4.);
  TwoPhotonInvMassAndTimeCutAndPhiCutAndMMCut[1] = new TH1D("TwoPhotonInvMassAndTimeCutAndPhiCutAndMMCut1", "Inv. Mass 2g with Cuts in Time, phi and MM", 800, 0., 4.);
  TwoPhotonInvMassAndTimeCutAndPhiCutAndMMCut[2] = new TH1D("TwoPhotonInvMassAndTimeCutAndPhiCutAndMMCut2", "Inv. Mass 2g with Cuts in Time, phi and MM", 800, 0., 4.);
  TH2D *TwoGPHIvsTwoGInvMassTimeCutPhiCutMMCut[3];
  TwoGPHIvsTwoGInvMassTimeCutPhiCutMMCut[0] = new TH2D("TwoGPHIvsTwoGInvMassTimeCutPhiCutMMCut0", "2g Phi vs Inv. Mass w.all cuts",
	   								      100, 0., 4., 36, 0., 2.*TMath::Pi());
  TwoGPHIvsTwoGInvMassTimeCutPhiCutMMCut[1] = new TH2D("TwoGPHIvsTwoGInvMassTimeCutPhiCutMMCut1", "2g Phi vs Inv. Mass w.all cuts",
	   								      100, 0., 4., 36, 0., 2.*TMath::Pi());
  TwoGPHIvsTwoGInvMassTimeCutPhiCutMMCut[2] = new TH2D("TwoGPHIvsTwoGInvMassTimeCutPhiCutMMCut2", "2g Phi vs Inv. Mass w.all cuts",
									      100, 0., 4., 36, 0., 2.*TMath::Pi());
    

  TH1D *Pi0GammaInvMass = new TH1D("Pi0GammaInvMass", "pi0 gamma invariant mass", 800, 0., 4.);
  TH2D *StartCounterIDvsInvMassPi0Gamma = new TH2D("StartCounterIDvsInvMassPi0Gamma", "SC ID vs Inv. Mass of pi0#gamma system", 800, 0., 4., 32, 0., 32.);
  TH2D *StartCounterTimeVsInvMassPi0Gamma = new TH2D("StartCounterTimeVsInvMassPi0Gamma", "SC time vs Inv. Mass of pi0#gamma system", 800, 0., 4., 500, -100., 100. );
  TH2D *StartCounterEnergyVsInvMassPi0Gamma = new TH2D("StartCounterEnergyVsInvMassPi0Gamma", "SC energy vs Inv. Mass of pi0#gamma system", 800, 0., 4., 100, 0., 0.015);
  TH1D *StartCounterRelTimeToPi0Gamma = new TH1D("StartCounterRelTimeToPi0Gamma", "SC time relative to pi0#gamma system", 500, -100., 100.);
  TH2D *StartCounterEnergyVsInvMassPi0GammaAndTimeCut = new TH2D("StartCounterEnergyVsInvMassPi0GammaAndTimeCut", 
								 "SC Energy vs Inv.Mass of pi0#gamma system with time cut", 800, 0., 4., 100, 0., 0.015);
  TH2D *StartCounterTimeVsInvMassPi0GammaAndTimeCut = new TH2D("StartCounterTimeVsInvMassPi0GammaAndTimeCut", 
							       "SC time vs Inv. Mass of pi0#gamma system with time cut ", 800, 0., 4., 500, -100., 100. );
  TH1D *NumberOfStartCounterHitsWithPi0GammaTimeCut = new TH1D("NumberOfStartCounterHitsWithPi0GammaTimeCut","SC Hit multiplicity with good time for pi0#gamma system", 10, 0., 10.);
  TH1D *Pi0GammaInvMassAndTimeCutAndOneSCHitAndSCEnergyCut = new TH1D("Pi0GammaInvMassAndTimeCutAndOneSCHitAndSCEnergyCut", "Inv. Mass of pi0#gamma system with 1 SC Hit and ECut", 800, 0., 4.);
  TH1D *Pi0GammaPhiAngleAndTimeCutAndOneSCHitAndSCEnergyCut = new TH1D("Pi0GammaPhiAngleAndTimeCutAndOneSCHitAndSCEnergyCut",
								       "Phi of pi0#gamma system with 1 SC Hit and ECut", 800, 0., 2.*TMath::Pi());
  TH1D *Pi0GammaThetaAngleAndTimeCutAndOneSCHitAndSCEnergyCut = new TH1D("Pi0GammaThetaAngleAndTimeCutAndOneSCHitAndSCEnergyCut", 
									 "Theta of pi0#gamma system with 1 SC Hit and ECut", 400, 0, 0.5);

  TH1D *MissingMassPi0Gamma[3];
  MissingMassPi0Gamma[0] = new TH1D("MissingMassPi0Gamma0", "Missing Mass of pi0#gamma system Background", 800, -4., 8.);
  MissingMassPi0Gamma[1] = new TH1D("MissingMassPi0Gamma1", "Missing Mass of pi0#gamma system Signal", 800, -4., 8.);
  MissingMassPi0Gamma[2] = new TH1D("MissingMassPi0Gamma2", "Missing Mass of pi0#gamma system Signal w. Bgr. subtracted", 800, -4., 8.);
  TH1D *MissingMomentumPi0Gamma[3];
  MissingMomentumPi0Gamma[0] = new TH1D("MissingMomentumPi0Gamma0", "Missing Momentum of pi0#gamma system Background", 800, -4., 8.);
  MissingMomentumPi0Gamma[1] = new TH1D("MissingMomentumPi0Gamma1", "Missing Momentum of pi0#gamma system Signal", 800, -4., 8.);
  MissingMomentumPi0Gamma[2] = new TH1D("MissingMomentumPi0Gamma2", "Missing Momentum of pi0#gamma system Signal w. Bgr. subtracted", 800, -4., 8.);
  TH1D *MissingEnergyPi0Gamma[3];
  MissingEnergyPi0Gamma[0] = new TH1D("MissingEnergyPi0Gamma0", "Missing Energy of pi0#gamma system Background", 800, -4., 8.);
  MissingEnergyPi0Gamma[1] = new TH1D("MissingEnergyPi0Gamma1", "Missing Energy of pi0#gamma system Signal", 800, -4., 8.);
  MissingEnergyPi0Gamma[2] = new TH1D("MissingEnergyPi0Gamma2", "Missing Energy of pi0#gamma system Signal w. Bgr. subtracted", 800, -4., 8.);

  TH2D *MomentumVsMissingMassPi0Gamma[3];
  MomentumVsMissingMassPi0Gamma[0] = new TH2D("MomentumVsMissingMassPi0Gamma0","Momentum vs MM of pi0#gamma system Background", 800, -4., 8., 800, -4., 8. );
  MomentumVsMissingMassPi0Gamma[1] = new TH2D("MomentumVsMissingMassPi0Gamma1","Momentum vs MM of pi0#gamma system Signal", 800, -4., 8., 800, -4., 8. );
  MomentumVsMissingMassPi0Gamma[2] = new TH2D("MomentumVsMissingMassPi0Gamma2","Momentum vs MM of pi0#gamma system Signal w. Bgr. subtracted", 800, -4., 8., 800, -4., 8. );
  TH1D *MissingpPhiPi0Gamma[3];
  MissingpPhiPi0Gamma[0] = new TH1D("MissingpPhiPi0Gamma0", "Phi of missing momentum pi0#gamma system Background", 36, 0., 2.*TMath::Pi());
  MissingpPhiPi0Gamma[1] = new TH1D("MissingpPhiPi0Gamma1", "Phi of missing momentum pi0#gamma system Signal", 36, 0., 2.*TMath::Pi());
  MissingpPhiPi0Gamma[2] = new TH1D("MissingpPhiPi0Gamma2", "Phi of missing momentum pi0#gamma system Signal w. Bgr. subtracted", 36, 0., 2.*TMath::Pi());
  TH1D *DeltaPhiPi0Gamma[3];
  DeltaPhiPi0Gamma[0] = new TH1D("DeltaPhiPi0Gamma0", "Missing Pphi minus SC phi for pi0#gamma system Background", 72, -TMath::Pi(), TMath::Pi());
  DeltaPhiPi0Gamma[1] = new TH1D("DeltaPhiPi0Gamma1", "Missing Pphi minus SC phi for pi0#gamma system Signal", 72, -TMath::Pi(), TMath::Pi());
  DeltaPhiPi0Gamma[2] = new TH1D("DeltaPhiPi0Gamma2", "Missing Pphi minus SC phi for pi0#gamma system Signal w. Bgr. subtracted", 72, -TMath::Pi(), TMath::Pi());

  TH1D *Pi0GammaInvMassAndTimeCutAndPhiCut[3];
  Pi0GammaInvMassAndTimeCutAndPhiCut[0] = new TH1D("Pi0GammaInvMassAndTimeCutAndPhiCut0", "pi0#gamma Inv. Mass with cuts on Time,phi Background", 800, 0., 4.);
  Pi0GammaInvMassAndTimeCutAndPhiCut[1] = new TH1D("Pi0GammaInvMassAndTimeCutAndPhiCut1", "pi0#gamma Inv. Mass with cuts on Time,phi Signal", 800, 0., 4.);
  Pi0GammaInvMassAndTimeCutAndPhiCut[2] = new TH1D("Pi0GammaInvMassAndTimeCutAndPhiCut2", "pi0#gamma Inv. Mass with cuts on Time,phi Signal w. Bgr. subtracted", 800, 0., 4.);

  TH1D *Pi0GammaInvMassAndTimeCutAndPhiCutAndMMCut[3];
  Pi0GammaInvMassAndTimeCutAndPhiCutAndMMCut[0] = new TH1D("Pi0GammaInvMassAndTimeCutAndPhiCutAndMMCut0", "pi0#gamma Inv. Mass with cuts on Time,phi,MM Background", 800, 0., 4.);
  Pi0GammaInvMassAndTimeCutAndPhiCutAndMMCut[1] = new TH1D("Pi0GammaInvMassAndTimeCutAndPhiCutAndMMCut1", "pi0#gamma Inv. Mass with cuts on Time,phi,MM Signal", 800, 0., 4.);
  Pi0GammaInvMassAndTimeCutAndPhiCutAndMMCut[2] = new TH1D("Pi0GammaInvMassAndTimeCutAndPhiCutAndMMCut2", "pi0#gamma Inv. Mass with cuts on Time,phi,MM Signal w. Bgr. subtracted", 800, 0., 4.);
  TH2D *Pi0GammaInvMassVsPi0GammaPhiWithTimeCutAndPhiCutAndMMCut[3];
  Pi0GammaInvMassVsPi0GammaPhiWithTimeCutAndPhiCutAndMMCut[0] = new TH2D("Pi0GammaInvMassVsPi0GammaPhiWithTimeCutAndPhiCutAndMMCut0", 
									 "phi vs Inv. Mass of pi0#gamma system Background", 400, 0., 4., 36, 0., 2.*TMath::Pi());
  Pi0GammaInvMassVsPi0GammaPhiWithTimeCutAndPhiCutAndMMCut[1] = new TH2D("Pi0GammaInvMassVsPi0GammaPhiWithTimeCutAndPhiCutAndMMCut1", 
									 "phi vs Inv. Mass of pi0#gamma system Background", 400, 0., 4., 36, 0., 2.*TMath::Pi());
  Pi0GammaInvMassVsPi0GammaPhiWithTimeCutAndPhiCutAndMMCut[2] = new TH2D("Pi0GammaInvMassVsPi0GammaPhiWithTimeCutAndPhiCutAndMMCut2", 
									 "phi vs Inv. Mass of pi0#gamma system Background", 400, 0., 4., 36, 0., 2.*TMath::Pi());


  TH1D *omegaDecayPlaneAllCuts[3];
  omegaDecayPlaneAllCuts[0] = new TH1D("omegaDecayPlaneAllCuts0", "Decay plane of omega #phi Background", 36., 0., 2.*TMath::Pi());
  omegaDecayPlaneAllCuts[1] = new TH1D("omegaDecayPlaneAllCuts1", "Decay plane of omega #phi Signal", 36., 0., 2.*TMath::Pi());
  omegaDecayPlaneAllCuts[2] = new TH1D("omegaDecayPlaneAllCuts2", "Decay plane of omega #phi Signal Bgr. subtr.", 36., 0., 2.*TMath::Pi());

  TH1D *omegapi0DecayAngleAllCuts[3]; 
  omegapi0DecayAngleAllCuts[0] = new TH1D("omegapi0DecayAngleAllCuts0", "#pi^{0} decay angle #theta in omega CM Frame (helicity F)  Background", 180, 0., TMath::Pi());
  omegapi0DecayAngleAllCuts[1] = new TH1D("omegapi0DecayAngleAllCuts1", "#pi^{0} decay angle #theta in omega CM Frame (helicity F)  Signal", 180, 0., TMath::Pi());
  omegapi0DecayAngleAllCuts[2] = new TH1D("omegapi0DecayAngleAllCuts2", "#pi^{0} decay angle #theta in omega CM Frame (helicity F)  Signal Bgr. subtr.", 180, 0., TMath::Pi());

  TH1D *omegapi0PSIAllCuts[3]; 
  omegapi0PSIAllCuts[0] = new TH1D("omegapi0PSIAllCuts0", "#psi=#Phi-#phi Background", 16, 0., 2.*TMath::Pi());
  omegapi0PSIAllCuts[1] = new TH1D("omegapi0PSIAllCuts1", "#psi=#Phi-#phi Signal", 16, 0., 2.*TMath::Pi());
  omegapi0PSIAllCuts[2] = new TH1D("omegapi0PSIAllCuts2", "#psi=#Phi-#phi Signal Bgr. subtr.", 16, 0., 2.*TMath::Pi());

  TH1D *omegaTHETAAllCuts[3]; 
  omegaTHETAAllCuts[0] = new TH1D("omegaTHETAAllCuts0", "#theta omega Background", 100, 0., 0.1);
  omegaTHETAAllCuts[1] = new TH1D("omegaTHETAAllCuts1", "#theta omega Signal", 100, 0., 0.1);
  omegaTHETAAllCuts[2] = new TH1D("omegaTHETAAllCuts2", "#theta omega Signal Bgr. subtr.", 100, 0., 0.1);

  TH1D *omegapi0XSIAllCuts[3]; 
  omegapi0XSIAllCuts[0] = new TH1D("omegapi0XSIAllCuts0", "#psi Background", 36, 0., 2.*TMath::Pi());
  omegapi0XSIAllCuts[1] = new TH1D("omegapi0XSIAllCuts1", "#psi Signal", 36, 0., 2.*TMath::Pi());
  omegapi0XSIAllCuts[2] = new TH1D("omegapi0XSIAllCuts2", "#psi Signal Bgr. subtr.", 36, 0., 2.*TMath::Pi());

  TH1D *omegaPHIAllCuts[3]; 
  omegaPHIAllCuts[0] = new TH1D("omegaPHIAllCuts0", "#Phi omega production plane Background", 36, 0., 2.*TMath::Pi());
  omegaPHIAllCuts[1] = new TH1D("omegaPHIAllCuts1", "#Phi omega production plane  Signal", 36, 0., 2.*TMath::Pi());
  omegaPHIAllCuts[2] = new TH1D("omegaPHIAllCuts2", "#Phi omega production plane  Signal Bgr. subtr.", 36, 0., 2.*TMath::Pi());

  TH1D *omega_t_AndAllCuts[3];
  omega_t_AndAllCuts[0] = new TH1D("omega_t_AndAllCuts0", "Exclusive #omega -t  after all cuts Background", 50, 0., 0.2);
  omega_t_AndAllCuts[1] = new TH1D("omega_t_AndAllCuts1", "Exclusive #omega -t  after all cuts Signal", 50, 0., 0.2);
  omega_t_AndAllCuts[2] = new TH1D("omega_t_AndAllCuts2", "Exclusive #omega -t  after all cuts Signal Background subtracted", 50, 0., 0.2);


  TH1D *StartCounterRelTimeWith4Gammas = new TH1D("StartCounterRelTimeWith4Gammas", "SC relative time with 4#gamma system", 500, -50., 50.);
  TH2D *TwoGammaVsTwoGammaInvMass = new TH2D("TwoGammaVsTwoGammaInvMass", "2#gamma vs 2#gamma inv. mass", 400, 0., 2., 400, 0., 2.);
  TH1D *MissingMass4Gamma[3];
  MissingMass4Gamma[0] = new TH1D("MissingMass4Gamma0", "Missing Mass of 4#gamma system Background", 100, -4., 8.);
  MissingMass4Gamma[1] = new TH1D("MissingMass4Gamma1", "Missing Mass of 4#gamma system Signal", 100, -4., 8.);
  MissingMass4Gamma[2] = new TH1D("MissingMass4Gamma2", "Missing Mass of 4#gamma system Signal w. Bgr. subtracted", 100, -4., 8.);
  TH1D *MissingMomentum4Gamma[3];
  MissingMomentum4Gamma[0] = new TH1D("MissingMomentum4Gamma0", "Missing Momentum of 4#gamma system Background", 100, -4., 8. );
  MissingMomentum4Gamma[1] = new TH1D("MissingMomentum4Gamma1", "Missing Momentum of 4#gamma system Signal", 100, -4., 8. );
  MissingMomentum4Gamma[2] = new TH1D("MissingMomentum4Gamma2", "Missing Momentum of 4#gamma system Signal w. Bgr. subtracted", 100, -4., 8. );
  TH1D *MissingEnergy4Gamma[3];
  MissingEnergy4Gamma[0] = new TH1D("MissingEnergy4Gamma0", "Missing Energy of 4#gamma system Background", 100, -4., 8. );
  MissingEnergy4Gamma[1] = new TH1D("MissingEnergy4Gamma1", "Missing Energy of 4#gamma system Signal", 100, -4., 8. );
  MissingEnergy4Gamma[2] = new TH1D("MissingEnergy4Gamma2", "Missing Energy of 4#gamma system Signal w. Bgr. subtracted", 100, -4., 8. );
  TH2D *MomentumVsMissingMass4Gamma[3];
  MomentumVsMissingMass4Gamma[0] = new TH2D("MomentumVsMissingMass4Gamma0", "MP vs MM of 4#gamma system Background", 100, -4., 8., 100, -4., 8.);
  MomentumVsMissingMass4Gamma[1] = new TH2D("MomentumVsMissingMass4Gamma1", "MP vs MM of 4#gamma system Signal", 100, -4., 8., 100, -4., 8.);
  MomentumVsMissingMass4Gamma[2] = new TH2D("MomentumVsMissingMass4Gamma2", "MP vs MM of 4#gamma system Signal w. Bgr. subtracted", 100, -4., 8., 100, -4., 8.);
  TH1D *MissingMomentumPhi4Gamma[3];
  MissingMomentumPhi4Gamma[0] = new TH1D("MissingMomentumPhi4Gamma0", "MPPhi of 4#gamma system Background", 36, 0., 2.*TMath::Pi());
  MissingMomentumPhi4Gamma[1] = new TH1D("MissingMomentumPhi4Gamma1", "MPPhi of 4#gamma system Signal", 36, 0., 2.*TMath::Pi());
  MissingMomentumPhi4Gamma[2] = new TH1D("MissingMomentumPhi4Gamma2", "MPPhi of 4#gamma system Signal w. Bgr. subtracted", 36, 0., 2.*TMath::Pi());
  TH1D *DeltaPhi4Gammas[3];
  DeltaPhi4Gammas[0] = new TH1D("DeltaPhi4Gammas0","Missing momentum phi minus SC phi Background", 72,  -TMath::Pi(), TMath::Pi());
  DeltaPhi4Gammas[1] = new TH1D("DeltaPhi4Gammas1","Missing momentum phi minus SC phi Signal", 72,  -TMath::Pi(), TMath::Pi());
  DeltaPhi4Gammas[2] = new TH1D("DeltaPhi4Gammas2","Missing momentum phi minus SC phi Signal w. Bgr. subtracted", 72,  -TMath::Pi(), TMath::Pi());
  TH1D *FourGammaInvMassAndMMCutAndPhiCut[3];
  FourGammaInvMassAndMMCutAndPhiCut[0] = new TH1D("FourGammaInvMassAndMMCutAndPhiCut0", "Inv. Mass of 4#gamma system with Cut on MM, and phi Background", 800, 0., 4.);
  FourGammaInvMassAndMMCutAndPhiCut[1] = new TH1D("FourGammaInvMassAndMMCutAndPhiCut1", "Inv. Mass of 4#gamma system with Cut on MM, and phi Signal", 800, 0., 4.);
  FourGammaInvMassAndMMCutAndPhiCut[2] = new TH1D("FourGammaInvMassAndMMCutAndPhiCut2", "Inv. Mass of 4#gamma system with Cut on MM, and phi Signal w. Bgr. subtracted", 800, 0., 4.);


   TH1D *f1270PHI[3];
   f1270PHI[0] = new TH1D("f1270PHI0", "azimuthal angular distribution of f1270", 16, 0., 2.*TMath::Pi());
   f1270PHI[1] = new TH1D("f1270PHI1", "azimuthal angular distribution of f1270", 16, 0., 2.*TMath::Pi());
   f1270PHI[2] = new TH1D("f1270PHI2", "azimuthal angular distribution of f1270", 16, 0., 2.*TMath::Pi());

   TH1D *f1270Theta[3];
   f1270Theta[0] = new TH1D("f1270Theta0", "#theta distribution of f1270", 100, 0., 0.1);
   f1270Theta[1] = new TH1D("f1270Theta1", "#theta angular distribution of f1270", 100, 0., 0.1);
   f1270Theta[2] = new TH1D("f1270Theta2", "#theta angular distribution of f1270", 100, 0., 0.1);

   TH1D *f1270t[3];
   f1270t[0] = new TH1D("f1270t0", "-t of f1270 region", 50, 0., 0.2);
   f1270t[1] = new TH1D("f1270t1", "-t of f1270 region", 50, 0., 0.2);
   f1270t[2] = new TH1D("f1270t2", "-t of f1270 region", 50, 0., 0.2);


   TH1D *f1270DecayPlaneAllCuts[3];
   f1270DecayPlaneAllCuts[0] = new TH1D("f1270DecayPlaneAllCuts0", "#phi f1270 decay plane Bg.",  36., 0., 2.*TMath::Pi());
   f1270DecayPlaneAllCuts[1] = new TH1D("f1270DecayPlaneAllCuts1", "#phi f1270 decay plane Signal",  36., 0., 2.*TMath::Pi());
   f1270DecayPlaneAllCuts[2] = new TH1D("f1270DecayPlaneAllCuts2", "#phi f1270 decay plane Signal-Bg.",  36., 0., 2.*TMath::Pi());

   TH1D *f1270PSIAllCuts[3];
   f1270PSIAllCuts[0] = new TH1D("f1270PSIAllCuts0","#psi=#Phi-#phi f1270 Bg.",  16., 0., 2.*TMath::Pi());
   f1270PSIAllCuts[1] = new TH1D("f1270PSIAllCuts1","#psi=#Phi-#phi f1270 Signal",  16., 0., 2.*TMath::Pi());
   f1270PSIAllCuts[2] = new TH1D("f1270PSIAllCuts2","#psi=#Phi-#phi f1270 Singla-Bg.",  16., 0., 2.*TMath::Pi());

   
   TH1D *f1270PHIAllCuts[3];
   f1270PHIAllCuts[0] = new TH1D("f1270PHIAllCuts0","#Phi f1270 Production plane Bg.",  16., 0., 2.*TMath::Pi());
   f1270PHIAllCuts[1] = new TH1D("f1270PHIAllCuts1","#Phi f1270 Production plane Signal",  16., 0., 2.*TMath::Pi());
   f1270PHIAllCuts[2] = new TH1D("f1270PHIAllCuts2","#Phi f1270 Production plane Signal-Bg.",  16., 0., 2.*TMath::Pi());
   
   TH1D *f1270pi0DecayAngleAllCuts[3];
   f1270pi0DecayAngleAllCuts[0] = new TH1D("f1270pi0DecayAngleAllCuts0", "#theta f1270 #pi0 decay angle Bg.", 32, 0., TMath::Pi() );
   f1270pi0DecayAngleAllCuts[1] = new TH1D("f1270pi0DecayAngleAllCuts1", "#theta f1270 #pi0 decay angle Signal", 32, 0., TMath::Pi() );
   f1270pi0DecayAngleAllCuts[2] = new TH1D("f1270pi0DecayAngleAllCuts2", "#theta f1270 #pi0 decay angle Signal-Bg.", 32, 0., TMath::Pi() );


  cout<<"...done"<<endl;

  // loop over events in the tree
  unsigned int nentries = (unsigned int) myChain->GetEntries();
  cout<<"Start reading root tree chain with "<<nentries<<" events"<<endl;
  //nentries = 100000;
  for (unsigned int k=0; k<nentries; k++){
    myChain->GetEntry(k);
    
    //cout<<"Entry: "<<k<<endl;

    // loop over beam photons find intime photons
    Good_beam_idx.clear();
    vector <TLorentzVector> BeamPhotons;
    vector <int> InTimePhotons;
    int InTimeCounter = 0;
    int OutOfTimeCounter = 0;

    
    for  (int n=0; n<Nbeam; n++){
      float dt = BeamT[n] - RFT;
      beamphotontime->Fill(dt);

      if (TMath::Abs(BeamE[n]-BeamEnergyCutCenter)<BeamEnergyCutWidth){

	if (TMath::Abs(dt)<1.){

	  TLorentzVector lv(0.,0.,BeamE[n],BeamE[n]);
	  BeamPhotons.push_back(lv);
	  InTimePhotons.push_back(1);
	  beamphotons[1]->Fill(BeamE[n]);
	  InTimeCounter++;
	} else {
	  // only select tagger hits +/- 1ns round timing peaks
	  int ok =0;
	  double abs_dt = TMath::Abs(dt);
	  for (int j=1;j<5;j++) {
	    double tloc = 4.*(double)j;
	    if (TMath::Abs(abs_dt-tloc)<1.)      { // accidentals +/- four peaks
	      ok = 1;
	    }
	  }
	  if (ok) {
	    TLorentzVector lv(0.,0.,BeamE[n],BeamE[n]);
	    BeamPhotons.push_back(lv);
	    InTimePhotons.push_back(0);
	    beamphotons[0]->Fill(BeamE[n]);
	    OutOfTimeCounter++;
	  }
	}
      }
    }
    
    if (BeamPhotons.size()<1){
      //cout<<"no good beam photons"<<endl;
      continue;
    }

    NumberOfBeamPhotons[1]->Fill((float)InTimeCounter);      // beam photon E>BeamEnergyCut
    NumberOfBeamPhotons[0]->Fill((float)OutOfTimeCounter);   // beam photon E>BeamEnergyCut 

    int GoodBeamPhotons = 0;
    if (InTimeCounter){
      GoodBeamPhotons =1;    // Beam photon candidate with good timing and Energy over BeamEnergyCut
    }

    // Find START COUNTER Hits with energy larger than 0.5MeV AND resaonalble timing
    int FoundSCE = 0;
    int FoundSCT = 0;
    vector <int> SCIDX;
    for (int n=0; n<Nstart; n++){
      StartCounterTime->Fill(SC_T[n],SC_S[n]);
      StartCounterEnergy->Fill(SC_E[n],SC_S[n]);
      if (TMath::Abs(SC_T[n]-StartCounterTimePeak)<50.) {
	FoundSCT = 1;
	StarCounterEnergyWithTimingCut->Fill(SC_E[n],SC_S[n]);
      }
      if (SC_E[n] > 0.0005) {
	FoundSCE = 1;
	StartCountertimeWithEnergyCut->Fill(SC_T[n],SC_S[n]);
      }
      if ((TMath::Abs(SC_T[n]-StartCounterTimePeak)<50.) && (SC_E[n] > 0.0005)){
	SCIDX.push_back(n);
      }
    }
    if ((!FoundSCE) || (!FoundSCT)){
      //cout<<"not good SC hit"<<endl;
      continue;
    }

    
    // find photons with energy above 0.05GeV
    vector <TLorentzVector> Gammas;    // final state Photon 4-Vectors
    vector <double> GammaTimes;        // timing for these photons
    vector <int> IdxList;              // photon id number in the event-tree
    int NP = findphotons(Gammas, GammaTimes, MinimumPhotonEnergy, IdxList);
    
    // require at least two good photons in the detector
    if (NP<2){
      //cout<<"less than two good photons"<<endl;
      continue;
    }
    
    //cout<<"found event with more than 2 photons"<<endl;
    // look only at events with no more than 3 photons (and at least two)
    if (NP <4) {
      
      for  (int n=0; n<NP-1; n++){ // 1st loop over photons	
	
	Good_SC_idx.clear();
	Good_gamma_idx.clear();
	
	for  (int j=n+1; j<NP; j++){ // 2nd loop over photons
	  
	  double dt = GammaTimes[n]-GammaTimes[j];
	  gammatimediff->Fill(dt);
	  
	  // require that the time difference of the two photons at the target is smaller than 1.5 ns
	  if (TMath::Abs(dt)<1.5) {
	    TLorentzVector Ptot = Gammas[n] + Gammas[j];
	    if (NP ==2 ) {
	      TwoPhotonOnlyInvMass->Fill(Ptot.M());
	    } 

	    TwoPhotonInvMass->Fill(Ptot.M());
	    
	    
	    Good_gamma_idx.push_back(n);
	    Good_gamma_idx.push_back(j);
	    
	    int TOK = 0;
	    int nhits=0;
	    double SC_phi = 99999.;
	    // loop over START COUNTER hits and find the ones that are in time with
	    // the average target time of the two photons
	    // this time is called reltime
	    double scE = 0.;
	    for (unsigned int l=0; l<SCIDX.size(); l++){
	      StartCounterNumVsInvMass->Fill(Ptot.M(),SC_S[SCIDX[l]]);
	      StartCounterTimeVsInvMass->Fill(Ptot.M(),SC_T[SCIDX[l]]);
	      StartCounterEnergyVsInvMass->Fill(Ptot.M(),SC_E[SCIDX[l]]);
	      double reltime = (GammaTimes[n]+GammaTimes[j])/2.  - SC_T[SCIDX[l]];
	      StarCounterRelativeTime->Fill(reltime);
	      StartCounterTimeVsMeanGammaTime->Fill((GammaTimes[n]+GammaTimes[j])/2. , SC_T[SCIDX[l]]);
	      if ( (reltime>-6.) && (reltime<-1.) ) {
		nhits++;
		TOK = 1;
		scE = SC_E[SCIDX[l]];
		StartCounterEnergyVsInvMassAndTimeCut->Fill(Ptot.M(),SC_E[SCIDX[l]]);
		StartCounterNumberVsStartCounterTimeAndTimeCut->Fill(SC_T[SCIDX[l]],SC_S[SCIDX[l]]);
		StartCounterTimeVsInvMassAndTimeCut->Fill(Ptot.M() , SC_T[SCIDX[l]]);
		SC_phi = SC_S[SCIDX[l]] * 3.1414926/16. ; 		
		Good_SC_idx.push_back(SCIDX[l]);
	      }
	    }
	    
	    if (TOK){
	      GoodStartCounterHitsAndTimeCut->Fill(nhits);
	    }
	    
	    if ((TOK) && (nhits == 1) && (scE>0.003)) {
	      if (NP == 2 ) {
		TwoPhotonOnlyInvMassAndTimeCut->Fill(Ptot.M());
	      } 
	      TwoPhotonInvMassAndTimeCut->Fill(Ptot.M());
	      
	      
	      double phi = GetPhiAngle(Ptot);
	      if (NP == 2 ) {
		TwoPhotonOnlyPhi->Fill(phi);
	      }
	      TwoPhotonPhi->Fill(phi);
	      
	      double theta = GetThetaAngle(Ptot);
	      if (NP == 2 ) {
		TwoPhotonOnlyTheta->Fill(theta);
	      }
	      TwoPhotonTheta->Fill(theta);
	      
	      TLorentzVector P3(0.,0.,0.,0.);
	      if (NP>2){
		int idx = 2;
		if ((n==0) && (j==2)){
		  idx = 1;
		}
		if (n==1){
		  idx = 0;
		}
		P3 = Gammas[idx];
	      }

	      //std::cout.flush();
	      //cout<<"2photon system, loop over all beam photons ...."<<endl;
	      for (unsigned int IDX=0;IDX<BeamPhotons.size();IDX++) {  // loop over all possible BeamPhotons
		
		double Weight = -0.125; // Eight beam bunches +/- 4 off zero peak
		if (InTimePhotons[IDX]){
		  Weight = 1.;
		}
		
		TLorentzVector BmP = BeamPhotons[IDX] + ProtonTarget - Ptot - P3;
		MissingMass[InTimePhotons[IDX]]->Fill(BmP.M());
		MissingMomentum[InTimePhotons[IDX]]->Fill(BmP.P());
		MissingEnergy[InTimePhotons[IDX]]->Fill(BmP.E());    
		MomentumVsMissingMass[InTimePhotons[IDX]]->Fill(BmP.M(),BmP.P());

		MissingMass[2]->Fill(BmP.M(), Weight);
		MissingMomentum[2]->Fill(BmP.P(), Weight);
		MissingEnergy[2]->Fill(BmP.E(), Weight);    
		MomentumVsMissingMass[2]->Fill(BmP.M(),BmP.P(),Weight);
		
		double Bestphi = GetPhiAngle(BmP);
		MissingPhi[InTimePhotons[IDX]]->Fill(Bestphi);
		MissingPhi[2]->Fill(Bestphi, Weight);
		
		double deltaphi = Bestphi - SC_phi;
		
		DeltaPhi[InTimePhotons[IDX]]->Fill(deltaphi);
		DeltaPhi[2]->Fill(deltaphi, Weight);
		if (TMath::Abs(deltaphi)<ThePhiCut) {
		  if (NP == 2 ) {
		    TwoPhotonOnlyInvMassAndTimeCutAndPhiCut[InTimePhotons[IDX]]->Fill(Ptot.M());
		    TwoPhotonOnlyMissingMassAndTimeCutAndPhiCut[InTimePhotons[IDX]]->Fill(BmP.M());

		    TwoPhotonOnlyInvMassAndTimeCutAndPhiCut[2]->Fill(Ptot.M(), Weight);
		    TwoPhotonOnlyMissingMassAndTimeCutAndPhiCut[2]->Fill(BmP.M(), Weight);
		    if ((BmP.M()<1.5) && (BmP.M()>0.) ) {
		      TwoPhotonOnlyInvMassAndTimeCutAndPhiCutAndMMCut[InTimePhotons[IDX]]->Fill(Ptot.M());
		      TwoPhotonOnlyPHIvsTwoPhotonInvMassAndTimeCutAndPhiCutAndMMCut[InTimePhotons[IDX]]->Fill(Ptot.M(),phi);

		      TwoPhotonOnlyInvMassAndTimeCutAndPhiCutAndMMCut[2]->Fill(Ptot.M(), Weight);
		      TwoPhotonOnlyPHIvsTwoPhotonInvMassAndTimeCutAndPhiCutAndMMCut[2]->Fill(Ptot.M(),phi, Weight);
		      
		      if (TMath::Abs(Ptot.M()-0.132)<0.02){
			TwoPhotonOnlyPi0PhiAndAllCuts[InTimePhotons[IDX]]->Fill(phi);
			TwoPhotonOnlyPi0PhiAndAllCuts[2]->Fill(phi, Weight);

			TwoPhotonOnlyPi0ThetaAndAllCuts[InTimePhotons[IDX]]->Fill(theta);
			TwoPhotonOnlyPi0ThetaAndAllCuts[2]->Fill(theta, Weight);

			double t = -(BeamPhotons[IDX]  - Ptot).M2();
			TwoPhotonOnlyPi0_t_AndAllCuts[InTimePhotons[IDX]]->Fill(t);
			TwoPhotonOnlyPi0_t_AndAllCuts[2]->Fill(t, Weight);

		      }
		    }
		  }
		  TwoPhotonInvMassAndTimeCutAndPhiCut[InTimePhotons[IDX]]->Fill(Ptot.M());
		  TwoPhotonMissingMassAndTimeCutAndPhiCut[InTimePhotons[IDX]]->Fill(BmP.M());
		  
		  TwoPhotonInvMassAndTimeCutAndPhiCut[2]->Fill(Ptot.M(), Weight);
		  TwoPhotonMissingMassAndTimeCutAndPhiCut[2]->Fill(BmP.M(), Weight);
		  if ((BmP.M()<1.5) && (BmP.M()>0.) ) {
		    TwoPhotonInvMassAndTimeCutAndPhiCutAndMMCut[InTimePhotons[IDX]]->Fill(Ptot.M());
		    TwoGPHIvsTwoGInvMassTimeCutPhiCutMMCut[InTimePhotons[IDX]]->Fill(Ptot.M(),phi);
		    
		    TwoPhotonInvMassAndTimeCutAndPhiCutAndMMCut[2]->Fill(Ptot.M(), Weight);
		    TwoGPHIvsTwoGInvMassTimeCutPhiCutMMCut[2]->Fill(Ptot.M(),phi, Weight);

		    if (TMath::Abs(Ptot.M()-0.135)<0.02){
		      TwoPhotonPi0PhiAndAllCuts[InTimePhotons[IDX]]->Fill(phi);
		      TwoPhotonPi0PhiAndAllCuts[2]->Fill(phi, Weight);
		    }
		  }
		}		
	      }
	      //cout<<"...done"<<endl;
	    }
	    
	    // if there are more than two photons test if these two photons match the pi0 mass
	    // if yes then look at the third photon
	    //cout<<"case of more than two photons"<<endl;
	    if (NP > 2) {
	      if (TMath::Abs(Ptot.M()-0.135) < 0.02){
		int i3 = 2;
		if ((n == 0) && (j==n+2))  {
		  i3 = 1;
		}
		if ((n == 1) && (j==n+1))  {
		  i3 = 0;
		}
		
		// check if the third photon is in time with the average time of the other two
		double dt3 =  (GammaTimes[n]+GammaTimes[j])/2. - GammaTimes[i3];
		if (TMath::Abs(dt3)<1.5) {
		  TLorentzVector Ptot3 = Gammas[n] + Gammas[j] + Gammas[i3];
		  Pi0GammaInvMass->Fill(Ptot3.M());
		  
		  int TOK3 = 0;
		  int nhits3=0;
		  double SC_phi3 = 9999.;
		  // find START COUNTER hits that are in time with the average of all three photons
		  double scE = 0;
		  for (unsigned int l=0; l<SCIDX.size(); l++){
		    StartCounterIDvsInvMassPi0Gamma->Fill(Ptot3.M(),SC_S[SCIDX[l]]);
		    StartCounterTimeVsInvMassPi0Gamma->Fill(Ptot3.M(),SC_T[SCIDX[l]]);
		    StartCounterEnergyVsInvMassPi0Gamma->Fill(Ptot3.M(),SC_E[SCIDX[l]]);
		    double reltime3 =  (GammaTimes[0] + GammaTimes[1] + GammaTimes[2])/3. - SC_T[SCIDX[l]];
		    StartCounterRelTimeToPi0Gamma->Fill(reltime3);
		    if ( (reltime3>-6.) && (reltime3<-0.5)){
		      nhits3++;
		      TOK3 = 1;
		      scE = SC_E[SCIDX[l]];
		      StartCounterEnergyVsInvMassPi0GammaAndTimeCut->Fill(Ptot3.M(),SC_E[SCIDX[l]]);
		      StartCounterTimeVsInvMassPi0GammaAndTimeCut->Fill(Ptot3.M(),SC_T[SCIDX[l]]);
		      SC_phi3 = SC_S[SCIDX[l]] * 3.1414926/16. ;
		    }
		  }
		  if (TOK3){
		    NumberOfStartCounterHitsWithPi0GammaTimeCut->Fill(nhits3);
		  }
		  if ((TOK3) && (nhits3==1) && (scE>0.004)) {
		    Pi0GammaInvMassAndTimeCutAndOneSCHitAndSCEnergyCut->Fill(Ptot3.M());
		    double phi = GetPhiAngle(Ptot3);
		    Pi0GammaPhiAngleAndTimeCutAndOneSCHitAndSCEnergyCut->Fill(phi);
		    double theta3 = GetThetaAngle(Ptot3);
		    Pi0GammaThetaAngleAndTimeCutAndOneSCHitAndSCEnergyCut->Fill(theta3);
		    
		    for (unsigned int IDX=0;IDX<BeamPhotons.size();IDX++) {  // loop over all possible BeamPhotons
		      

		      double Weight = -0.125;
		      if (InTimePhotons[IDX]){
			Weight = 1.;
		      }

		      TLorentzVector BmP = BeamPhotons[IDX] + ProtonTarget - Ptot3;
		      MissingMassPi0Gamma[InTimePhotons[IDX]]->Fill(BmP.M());
		      MissingMomentumPi0Gamma[InTimePhotons[IDX]]->Fill(BmP.P());
		      MissingEnergyPi0Gamma[InTimePhotons[IDX]]->Fill(BmP.E());    
		      MomentumVsMissingMassPi0Gamma[InTimePhotons[IDX]]->Fill(BmP.M(),BmP.P());
		      

		      MissingMassPi0Gamma[2]->Fill(BmP.M(), Weight);
		      MissingMomentumPi0Gamma[2]->Fill(BmP.P(), Weight);
		      MissingEnergyPi0Gamma[2]->Fill(BmP.E(), Weight);    
		      MomentumVsMissingMassPi0Gamma[2]->Fill(BmP.M(),BmP.P(), Weight);

		      double Bestphi = GetPhiAngle(BmP);
		      MissingpPhiPi0Gamma[InTimePhotons[IDX]]->Fill(Bestphi);
		      MissingpPhiPi0Gamma[2]->Fill(Bestphi, Weight);
		      
		      double deltaphi3 = Bestphi - SC_phi3;
		      DeltaPhiPi0Gamma[InTimePhotons[IDX]]->Fill(deltaphi3);
		      DeltaPhiPi0Gamma[2]->Fill(deltaphi3, Weight);
		      if (TMath::Abs(deltaphi3)<ThePhiCut) {
			Pi0GammaInvMassAndTimeCutAndPhiCut[InTimePhotons[IDX]]->Fill(Ptot3.M());	
			Pi0GammaInvMassAndTimeCutAndPhiCut[2]->Fill(Ptot3.M(), Weight);	
			if ( (BmP.M()<1.5) && (BmP.M()>0.) ) {
			  Pi0GammaInvMassAndTimeCutAndPhiCutAndMMCut[InTimePhotons[IDX]]->Fill(Ptot3.M());
			  Pi0GammaInvMassVsPi0GammaPhiWithTimeCutAndPhiCutAndMMCut[InTimePhotons[IDX]]->Fill(Ptot3.M(), phi);
			  Pi0GammaInvMassAndTimeCutAndPhiCutAndMMCut[2]->Fill(Ptot3.M(), Weight);
			  Pi0GammaInvMassVsPi0GammaPhiWithTimeCutAndPhiCutAndMMCut[2]->Fill(Ptot3.M(), phi, Weight);

			  if (TMath::Abs(Ptot3.M()-0.768)<0.036*2.5){ // cut on omega mass

			    int BENI = 1;
			    if (BENI) {
			      
			      //cout<<"do some boosting..."<<endl;
			      // First Calculate CM energy and Boost Vector
			      double Ecm2 = ProtonMass*ProtonMass + 2.*ProtonMass*BeamPhotons[IDX].E();
			      double Ecm = TMath::Sqrt(Ecm2);
			      double protonCMp2 = Ecm2 - ProtonMass*ProtonMass;
			      double PprotonCM = TMath::Sqrt(protonCMp2);
			      TLorentzVector ProtonCM = TLorentzVector(0., 0., -PprotonCM, Ecm);
			      TVector3 Boost = ProtonCM.BoostVector();
			      TLorentzRotation BOOST2CM(Boost);
			      
			      // Boost particles from LAB to CM
			      TLorentzVector pi0CM = BOOST2CM * Ptot;
			      TLorentzVector g1CM = BOOST2CM * Gammas[i3];
			      TLorentzVector omegaCM = BOOST2CM * Ptot3;
			      TLorentzVector BeamCM = BOOST2CM * BeamPhotons[IDX];

			      TVector3 Boost2Rest = -omegaCM.BoostVector();
			      TLorentzRotation BOOST2REST(Boost2Rest);
			      // Boost particles to omega rest frame
			      TLorentzVector pi0CM_Rest = BOOST2REST * pi0CM;
			      TLorentzVector g1CM_Rest = BOOST2REST * g1CM;
			      TLorentzVector BeamCM_Rest = BOOST2REST * BeamCM;

			      TVector3 y = (BeamPhotons[IDX].Vect().Unit().Cross(omegaCM.Vect().Unit())).Unit();
			      TVector3 z = omegaCM.Vect().Unit();
			      TVector3 x = y.Cross(z).Unit();

			      // now express pi0CM_Rest vector in this Helicity frame
			      TVector3 Pi0Vector( pi0CM_Rest.Vect().Dot(x),
						  pi0CM_Rest.Vect().Dot(y),
						  pi0CM_Rest.Vect().Dot(z));
			      
			      double phipol = 0;                           // *** phipol=0 means horizontal polarization plane.
			      TVector3 eps(cos(phipol), sin(phipol), 0.0); // beam polarization vector in lab

			      double Theta = Pi0Vector.Theta();
			      double phi = Pi0Vector.Phi(); // angle between produciton plane and decay plane
			      double PHI = TMath::ATan2( y.Dot(eps), BeamCM.Vect().Unit().Dot(eps.Cross(y)));
			      
			      double psi = PHI - phi;
			      if (psi < -TMath::Pi()){
				psi += 2.*TMath::Pi();
			      }
			      if (psi > TMath::Pi()){
				psi -= 2.*TMath::Pi();
			      }
			      
			      omegaTHETAAllCuts[InTimePhotons[IDX]]->Fill(theta3);
			      omegaTHETAAllCuts[2]->Fill(theta3, Weight);
			      
			      omegaDecayPlaneAllCuts[InTimePhotons[IDX]]->Fill(phi+TMath::Pi());
			      omegaDecayPlaneAllCuts[2]->Fill(phi+TMath::Pi(), Weight);
			      
			      omegapi0PSIAllCuts[InTimePhotons[IDX]]->Fill(psi+TMath::Pi());
			      omegapi0PSIAllCuts[2]->Fill(psi+TMath::Pi(), Weight);
			      
			      omegaPHIAllCuts[InTimePhotons[IDX]]->Fill(PHI+TMath::Pi());
			      omegaPHIAllCuts[2]->Fill(PHI+TMath::Pi(), Weight);
			      
			      omegapi0DecayAngleAllCuts[InTimePhotons[IDX]]->Fill(Theta);
			      omegapi0DecayAngleAllCuts[2]->Fill(Theta, Weight);


			      double t = -(BeamPhotons[IDX]  - Ptot3).M2();
			      omega_t_AndAllCuts[InTimePhotons[IDX]]->Fill(t);
			      omega_t_AndAllCuts[2]->Fill(t, Weight);

		
			    } else { // this is from Elton

			      // boost all into omega rest frame
			      TLorentzRotation ResonanceBoost( -Ptot3.BoostVector()); // Boost into omega rest frame
			      TLorentzVector beam_rest   = ResonanceBoost * BeamPhotons[IDX];
			      TLorentzVector recoil_rest = ResonanceBoost * BmP;
			      TLorentzVector pi0_rest    = ResonanceBoost * Ptot;
			      TLorentzVector g3_rest     = ResonanceBoost * Gammas[i3];

			      double phipol = 0;                           // *** phipol=0 means horizontal polarization plane.
			      TVector3 eps(cos(phipol), sin(phipol), 0.0); // beam polarization vector in lab

			      // Helicity Frame: y || to lab production plane
			      // z-axis is antiprallele to recoil_proton in the resonance rest frame
			      TVector3 y = (BeamPhotons[IDX].Vect().Unit().Cross(BmP.Vect().Unit())).Unit();
			      TVector3 z = -1.*recoil_rest.Vect().Unit();
			      TVector3 x = y.Cross(z).Unit();

			      // now express pi0_rest vector in this Helicity frame
			      TVector3 Pi0Vector( pi0_rest.Vect().Dot(x),
						  pi0_rest.Vect().Dot(y),
						  pi0_rest.Vect().Dot(z));

			      double Theta = Pi0Vector.Theta();
			      double phi = Pi0Vector.Phi(); // angle between produciton plane and decay plane
			      double PHI = TMath::ATan2( y.Dot(eps), BeamPhotons[IDX].Vect().Unit().Dot(eps.Cross(y)));
			      
			      double psi = PHI - phi;
			      if (psi < -TMath::Pi()){
				psi += 2.*TMath::Pi();
			      }
			      if (psi > TMath::Pi()){
				psi -= 2.*TMath::Pi();
			      }

			      omegaTHETAAllCuts[InTimePhotons[IDX]]->Fill(theta3);
			      omegaTHETAAllCuts[2]->Fill(theta3, Weight);

			      omegaDecayPlaneAllCuts[InTimePhotons[IDX]]->Fill(phi+TMath::Pi());
			      omegaDecayPlaneAllCuts[2]->Fill(phi+TMath::Pi(), Weight);
			      
			      omegapi0PSIAllCuts[InTimePhotons[IDX]]->Fill(psi+TMath::Pi());
			      omegapi0PSIAllCuts[2]->Fill(psi+TMath::Pi(), Weight);
			      
			      omegaPHIAllCuts[InTimePhotons[IDX]]->Fill(PHI+TMath::Pi());
			      omegaPHIAllCuts[2]->Fill(PHI+TMath::Pi(), Weight);
			      
			      omegapi0DecayAngleAllCuts[InTimePhotons[IDX]]->Fill(Theta);
			      omegapi0DecayAngleAllCuts[2]->Fill(Theta, Weight);


			      double t = -(BeamPhotons[IDX]  - Ptot3).M2();
			      omega_t_AndAllCuts[InTimePhotons[IDX]]->Fill(t);
			      omega_t_AndAllCuts[2]->Fill(t, Weight);


			    }

			  }			    
			}
		      }
		    }
		    
		  }
		  
		}
		
	      } // end of found good pi0
	    } // end of have more than two photons 
	  }
	  
	}
      } // end loop over all good gammas
    } else if (NP ==4) {
            
      int NgoodT = 0;
      double meandt=0;
      for  (int n=0; n<NP-1; n++){ // loop over all photons	
	for  (int j=n+1; j<NP; j++){ // loop over all photons	
	  double dt = TMath::Abs(GammaTimes[n]-GammaTimes[j]);
	  if (dt<1.5){
	    NgoodT++;
	    meandt += GammaTimes[n]+GammaTimes[j];	 
	  }	
	}
      }
      
      if (NgoodT>3) { //  4 photons are in time with each other.

	meandt /= (double)(NgoodT) * 2.0;
	
	int TOK4=0;
	double SC_phi4=99999.;
	double scE=0.;
	int nhits=0;
	double RELTIME=99999.;
	for (int l=0; l<Nstart; l++){
	  if (SC_E[l]>0.003) {
	    double reltime = meandt  - SC_T[l];
	    StartCounterRelTimeWith4Gammas->Fill(reltime);
	    if ( (reltime>-6.5) && (reltime<-1.5)){
	      nhits++;
	      TOK4 = 1;
	      scE = SC_E[l];
	      SC_phi4 = SC_S[l] * 3.1414926/16. ; 		
	      RELTIME = reltime;
	    }
	  }
	}


	if (nhits==1){

	  int foundone = 0; // make sure I count only one!
	  for  (int n=0; n<NP-1; n++){
	    for  (int j=n+1; j<NP; j++){
	      TLorentzVector Ptot1 = Gammas[n] + Gammas[j];
	      
	      for (int s=0; s<NP-1; s++){
		for  (int i=s+1; i<NP; i++){
		  if ((s!=n) && (j!=i)){
		    TLorentzVector Ptot2 = Gammas[s] + Gammas[i];
		    TwoGammaVsTwoGammaInvMass->Fill(Ptot1.M(),Ptot2.M());

		    if ( (TMath::Abs(Ptot1.M2()-0.018)<0.008) && (TMath::Abs(Ptot2.M2()-0.018)<0.008) && !foundone){
		      TLorentzVector Ptot4 = Ptot1 + Ptot2;
		      foundone = 1;
		      for (unsigned int IDX=0;IDX<BeamPhotons.size();IDX++) {  // loop over all possible BeamPhotons
			
			TLorentzVector BmP = BeamPhotons[IDX] + ProtonTarget - Ptot4;
			MissingMass4Gamma[InTimePhotons[IDX]]->Fill(BmP.M());
			MissingMomentum4Gamma[InTimePhotons[IDX]]->Fill(BmP.P());
			MissingEnergy4Gamma[InTimePhotons[IDX]]->Fill(BmP.E());    
			MomentumVsMissingMass4Gamma[InTimePhotons[IDX]]->Fill(BmP.M(),BmP.P());

			double Weight = -0.125;
			if (InTimePhotons[IDX]){
			  Weight = 1;
			}
			MissingMass4Gamma[2]->Fill(BmP.M(), Weight);
			MissingMomentum4Gamma[2]->Fill(BmP.P(), Weight);
			MissingEnergy4Gamma[2]->Fill(BmP.E(), Weight);    
			MomentumVsMissingMass4Gamma[2]->Fill(BmP.M(),BmP.P(), Weight);
			
			double Bestphi = GetPhiAngle(BmP);
			MissingMomentumPhi4Gamma[InTimePhotons[IDX]]->Fill(Bestphi);
			MissingMomentumPhi4Gamma[2]->Fill(Bestphi, Weight);
			
			double deltaphi4 = Bestphi - SC_phi4;
			DeltaPhi4Gammas[InTimePhotons[IDX]]->Fill(deltaphi4);
			DeltaPhi4Gammas[2]->Fill(deltaphi4, Weight);
			if ((BmP.M()<1.5) && (BmP.M()>0.) && (BmP.P()<1.) && (TMath::Abs(deltaphi4)<ThePhiCut)) {
			  FourGammaInvMassAndMMCutAndPhiCut[InTimePhotons[IDX]]->Fill(Ptot4.M());
			  FourGammaInvMassAndMMCutAndPhiCut[2]->Fill(Ptot4.M(), Weight);

			  if (TMath::Abs(Ptot4.M()-1.25)<0.15){

			    f1270PHI[InTimePhotons[IDX]]->Fill(GetPhiAngle(Ptot4));
			    f1270PHI[2]->Fill(GetPhiAngle(Ptot4), Weight);

			    f1270Theta[InTimePhotons[IDX]]->Fill(GetThetaAngle(Ptot4));
			    f1270Theta[2]->Fill(GetThetaAngle(Ptot4), Weight);



			    if (1){

			      // boost all into 2pi0 rest frame
			      TLorentzRotation ResonanceBoost( -Ptot4.BoostVector()); // Boost into omega rest frame
			      TLorentzVector beam_rest   = ResonanceBoost * BeamPhotons[IDX];
			      TLorentzVector recoil_rest = ResonanceBoost * BmP;
			      TLorentzVector pi0A_rest    = ResonanceBoost * Ptot1;
			      TLorentzVector pi0B_rest     = ResonanceBoost * Ptot2;

			      double phipol = 0;                           // *** phipol=0 means horizontal polarization plane.
			      TVector3 eps(cos(phipol), sin(phipol), 0.0); // beam polarization vector in lab

			      // Helicity Frame: y || to lab production plane
			      // z-axis is antiprallele to recoil_proton in the resonance rest frame
			      TVector3 y = (BeamPhotons[IDX].Vect().Unit().Cross(BmP.Vect().Unit())).Unit();
			      TVector3 z = -1.*recoil_rest.Vect().Unit();
			      TVector3 x = y.Cross(z).Unit();

			      // now express pi0_rest vector in this Helicity frame
			      TVector3 Pi0AVector( pi0A_rest.Vect().Dot(x),
						   pi0A_rest.Vect().Dot(y),
						   pi0A_rest.Vect().Dot(z));
			      TVector3 Pi0BVector( pi0B_rest.Vect().Dot(x),
						   pi0B_rest.Vect().Dot(y),
						   pi0B_rest.Vect().Dot(z));

			      double ThetaA = Pi0AVector.Theta();
			      double ThetaB = Pi0BVector.Theta();
			      double phi = Pi0AVector.Phi(); // angle between produciton plane and decay plane
			      double PHI = TMath::ATan2( y.Dot(eps), BeamPhotons[IDX].Vect().Unit().Dot(eps.Cross(y)));

			      double psi = PHI - phi;
			      if (psi < -TMath::Pi()){
				psi += 2.*TMath::Pi();
			      }
			      if (psi > TMath::Pi()){
				psi -= 2.*TMath::Pi();
			      }

			      f1270DecayPlaneAllCuts[InTimePhotons[IDX]]->Fill(phi+TMath::Pi());
			      f1270DecayPlaneAllCuts[2]->Fill(phi+TMath::Pi(), Weight);
			      
			      f1270PSIAllCuts[InTimePhotons[IDX]]->Fill(psi+TMath::Pi());
			      f1270PSIAllCuts[2]->Fill(psi+TMath::Pi(), Weight);
			      
			      f1270PHIAllCuts[InTimePhotons[IDX]]->Fill(PHI+TMath::Pi());
			      f1270PHIAllCuts[2]->Fill(PHI+TMath::Pi(), Weight);
			      
			      f1270pi0DecayAngleAllCuts[InTimePhotons[IDX]]->Fill(ThetaA);
			      f1270pi0DecayAngleAllCuts[2]->Fill(ThetaA, Weight);

			      //f1270pi0DecayAngleAllCuts[InTimePhotons[IDX]]->Fill(ThetaB);
			      //f1270pi0DecayAngleAllCuts[2]->Fill(ThetaB, Weight);

			      double t = -(BeamPhotons[IDX]  - Ptot4).M2();
			      f1270t[InTimePhotons[IDX]]->Fill(t);
			      f1270t[2]->Fill(t, Weight);

			    }

			  }

			}
		      }
		    }
		  }
		}
	      }	      
	    }
	  }
	}
      } 
    }
  }
  
  char of1[128];
  sprintf(of1, "%s", outputrootfile);
  TFile *outf = new TFile(of1,"RECREATE");

  beamphotontime->Write();
  beamphotons[0]->Write();
  beamphotons[1]->Write();
  
  NumberOfBeamPhotons[1]->Write();
  NumberOfBeamPhotons[0]->Write();
  StartCounterTime->Write();
  StartCounterEnergy->Write();
  StarCounterEnergyWithTimingCut->Write();
  StartCountertimeWithEnergyCut->Write();
  
  gammatimediff->Write();
  TwoPhotonOnlyInvMass->Write();
  StartCounterNumVsInvMass->Write();
  StartCounterTimeVsInvMass->Write();
  StartCounterEnergyVsInvMass->Write();
  StarCounterRelativeTime->Write();
  StartCounterTimeVsMeanGammaTime->Write();
  StartCounterEnergyVsInvMassAndTimeCut->Write();
  StartCounterNumberVsStartCounterTimeAndTimeCut->Write();
  StartCounterTimeVsInvMassAndTimeCut->Write();
  GoodStartCounterHitsAndTimeCut->Write();

  TwoPhotonOnlyInvMassAndTimeCut->Write();
  TwoPhotonOnlyPhi->Write();
  TwoPhotonOnlyTheta->Write();    
  MissingMass[0]->Write();
  MissingMass[1]->Write();
  MissingMass[2]->Write();
  MissingMomentum[0]->Write();
  MissingMomentum[1]->Write();
  MissingMomentum[2]->Write();
  MissingEnergy[0]->Write();
  MissingEnergy[1]->Write();
  MissingEnergy[2]->Write();
  MomentumVsMissingMass[0]->Write();
  MomentumVsMissingMass[1]->Write();
  MomentumVsMissingMass[2]->Write();
  MissingPhi[0]->Write();
  MissingPhi[1]->Write();
  MissingPhi[2]->Write();
  DeltaPhi[0]->Write();
  DeltaPhi[1]->Write();
  DeltaPhi[2]->Write();
  TwoPhotonOnlyInvMassAndTimeCutAndPhiCut[0]->Write();
  TwoPhotonOnlyInvMassAndTimeCutAndPhiCut[1]->Write();
  TwoPhotonOnlyInvMassAndTimeCutAndPhiCut[2]->Write();

  TwoPhotonOnlyInvMassAndTimeCutAndPhiCutAndMMCut[0]->Write();
  TwoPhotonOnlyInvMassAndTimeCutAndPhiCutAndMMCut[1]->Write();
  TwoPhotonOnlyInvMassAndTimeCutAndPhiCutAndMMCut[2]->Write();
  
  TwoPhotonOnlyMissingMassAndTimeCutAndPhiCut[0]->Write();
  TwoPhotonOnlyMissingMassAndTimeCutAndPhiCut[1]->Write();
  TwoPhotonOnlyMissingMassAndTimeCutAndPhiCut[2]->Write();

  TwoPhotonInvMassAndTimeCutAndPhiCut[0]->Write();
  TwoPhotonInvMassAndTimeCutAndPhiCut[1]->Write();
  TwoPhotonInvMassAndTimeCutAndPhiCut[2]->Write();

  TwoPhotonInvMassAndTimeCutAndPhiCutAndMMCut[0]->Write();
  TwoPhotonInvMassAndTimeCutAndPhiCutAndMMCut[1]->Write();
  TwoPhotonInvMassAndTimeCutAndPhiCutAndMMCut[2]->Write();

  TwoGPHIvsTwoGInvMassTimeCutPhiCutMMCut[0]->Write();
  TwoGPHIvsTwoGInvMassTimeCutPhiCutMMCut[1]->Write();
  TwoGPHIvsTwoGInvMassTimeCutPhiCutMMCut[2]->Write();

  TwoPhotonPi0PhiAndAllCuts[0]->Write();
  TwoPhotonPi0PhiAndAllCuts[1]->Write();
  TwoPhotonPi0PhiAndAllCuts[2]->Write();

  TwoPhotonOnlyPi0PhiAndAllCuts[0]->Write();
  TwoPhotonOnlyPi0PhiAndAllCuts[1]->Write();
  TwoPhotonOnlyPi0PhiAndAllCuts[2]->Write();

  TwoPhotonOnlyPi0ThetaAndAllCuts[0]->Write();   
  TwoPhotonOnlyPi0ThetaAndAllCuts[1]->Write();   
  TwoPhotonOnlyPi0ThetaAndAllCuts[2]->Write();   

  TwoPhotonOnlyPi0_t_AndAllCuts[0]->Write();
  TwoPhotonOnlyPi0_t_AndAllCuts[1]->Write();
  TwoPhotonOnlyPi0_t_AndAllCuts[2]->Write();

  Pi0GammaInvMass->Write();
  StartCounterIDvsInvMassPi0Gamma->Write();
  StartCounterTimeVsInvMassPi0Gamma->Write();
  StartCounterEnergyVsInvMassPi0Gamma->Write();
  StartCounterRelTimeToPi0Gamma->Write();
  StartCounterEnergyVsInvMassPi0GammaAndTimeCut->Write();
  StartCounterTimeVsInvMassPi0GammaAndTimeCut->Write();
  NumberOfStartCounterHitsWithPi0GammaTimeCut->Write();
  Pi0GammaInvMassAndTimeCutAndOneSCHitAndSCEnergyCut->Write();
  Pi0GammaInvMassAndTimeCutAndOneSCHitAndSCEnergyCut->Write();
  Pi0GammaPhiAngleAndTimeCutAndOneSCHitAndSCEnergyCut->Write();
  Pi0GammaPhiAngleAndTimeCutAndOneSCHitAndSCEnergyCut->Write();
  Pi0GammaThetaAngleAndTimeCutAndOneSCHitAndSCEnergyCut->Write();
  Pi0GammaThetaAngleAndTimeCutAndOneSCHitAndSCEnergyCut->Write();
  MissingMassPi0Gamma[0]->Write();
  MissingMassPi0Gamma[1]->Write();
  MissingMassPi0Gamma[2]->Write();
  MissingMomentumPi0Gamma[0]->Write();
  MissingMomentumPi0Gamma[1]->Write();
  MissingMomentumPi0Gamma[2]->Write();
  MissingEnergyPi0Gamma[0]->Write();
  MissingEnergyPi0Gamma[1]->Write();
  MissingEnergyPi0Gamma[2]->Write();
  MomentumVsMissingMassPi0Gamma[0]->Write();
  MomentumVsMissingMassPi0Gamma[1]->Write();
  MomentumVsMissingMassPi0Gamma[2]->Write();
  MissingpPhiPi0Gamma[0]->Write();
  MissingpPhiPi0Gamma[1]->Write();
  MissingpPhiPi0Gamma[2]->Write();
  DeltaPhiPi0Gamma[0]->Write();
  DeltaPhiPi0Gamma[1]->Write();
  DeltaPhiPi0Gamma[2]->Write();
  Pi0GammaInvMassAndTimeCutAndPhiCut[0]->Write();
  Pi0GammaInvMassAndTimeCutAndPhiCut[1]->Write();
  Pi0GammaInvMassAndTimeCutAndPhiCut[2]->Write();
  Pi0GammaInvMassAndTimeCutAndPhiCutAndMMCut[0]->Write();
  Pi0GammaInvMassAndTimeCutAndPhiCutAndMMCut[1]->Write();
  Pi0GammaInvMassAndTimeCutAndPhiCutAndMMCut[2]->Write();
  Pi0GammaInvMassVsPi0GammaPhiWithTimeCutAndPhiCutAndMMCut[0]->Write();
  Pi0GammaInvMassVsPi0GammaPhiWithTimeCutAndPhiCutAndMMCut[1]->Write();
  Pi0GammaInvMassVsPi0GammaPhiWithTimeCutAndPhiCutAndMMCut[2]->Write();

  
  omegaDecayPlaneAllCuts[0]->Write();
  omegaDecayPlaneAllCuts[1]->Write();
  omegaDecayPlaneAllCuts[2]->Write();
  omegapi0DecayAngleAllCuts[0]->Write();
  omegapi0DecayAngleAllCuts[1]->Write();
  omegapi0DecayAngleAllCuts[2]->Write();
  omegapi0PSIAllCuts[0]->Write();
  omegapi0PSIAllCuts[1]->Write();
  omegapi0PSIAllCuts[2]->Write();
  omegapi0XSIAllCuts[0]->Write();
  omegapi0XSIAllCuts[1]->Write();
  omegapi0XSIAllCuts[2]->Write();
  omegaPHIAllCuts[0]->Write();
  omegaPHIAllCuts[1]->Write();
  omegaPHIAllCuts[2]->Write();

  omegaTHETAAllCuts[0]->Write();
  omegaTHETAAllCuts[1]->Write();
  omegaTHETAAllCuts[2]->Write();

  omega_t_AndAllCuts[0]->Write();
  omega_t_AndAllCuts[1]->Write();
  omega_t_AndAllCuts[2]->Write();

  StartCounterRelTimeWith4Gammas->Write();
  TwoGammaVsTwoGammaInvMass->Write();
  MissingMass4Gamma[0]->Write();
  MissingMass4Gamma[1]->Write();
  MissingMass4Gamma[2]->Write();
  MissingMomentum4Gamma[0]->Write();
  MissingMomentum4Gamma[1]->Write();
  MissingMomentum4Gamma[2]->Write();
  MissingEnergy4Gamma[0]->Write();
  MissingEnergy4Gamma[1]->Write();
  MissingEnergy4Gamma[2]->Write();
  MomentumVsMissingMass4Gamma[0]->Write();
  MomentumVsMissingMass4Gamma[1]->Write();
  MomentumVsMissingMass4Gamma[2]->Write();
  MissingMomentumPhi4Gamma[0]->Write();
  MissingMomentumPhi4Gamma[1]->Write();
  MissingMomentumPhi4Gamma[2]->Write();
  DeltaPhi4Gammas[0]->Write();
  DeltaPhi4Gammas[1]->Write();
  DeltaPhi4Gammas[2]->Write();
  FourGammaInvMassAndMMCutAndPhiCut[0]->Write();
  FourGammaInvMassAndMMCutAndPhiCut[1]->Write();
  FourGammaInvMassAndMMCutAndPhiCut[2]->Write();

  f1270PHI[0]->Write();
  f1270PHI[1]->Write();
  f1270PHI[2]->Write();
 
  f1270Theta[0]->Write();
  f1270Theta[1]->Write();
  f1270Theta[2]->Write();

  f1270t[0]->Write();
  f1270t[1]->Write();
  f1270t[2]->Write();
  
  f1270DecayPlaneAllCuts[0]->Write();
  f1270DecayPlaneAllCuts[1]->Write();
  f1270DecayPlaneAllCuts[2]->Write();

  f1270PSIAllCuts[0]->Write();
  f1270PSIAllCuts[1]->Write();
  f1270PSIAllCuts[2]->Write();
  
  f1270PHIAllCuts[0]->Write();
  f1270PHIAllCuts[1]->Write();
  f1270PHIAllCuts[2]->Write();

  f1270pi0DecayAngleAllCuts[0]->Write();
  f1270pi0DecayAngleAllCuts[1]->Write();
  f1270pi0DecayAngleAllCuts[2]->Write();




  outf->Close();

  cout<<"ALL DONE."<<endl;
}



int findphotons(vector <TLorentzVector> &Gammas, vector <double> &GammaTimes, double EnergyCut, vector <int> &IdxList){

  // find all Neutrals with energies above EnergyCut and construct 4-vectors
  // assuming the vertex at the center of the target = 65cm

  int cnt = 0;
  for  (int n=0; n<Nneut; n++){
    if (NeutralE[n]>EnergyCut) {

      if (DetSys[n]){
	double trans = TMath::Sqrt(NeutralX[n]*NeutralX[n] + NeutralY[n]*NeutralY[n]);
	if (trans>FCAL_R_CUT){
	  continue;
	}
      }

      double R2 = NeutralX[n]*NeutralX[n] + NeutralY[n]*NeutralY[n] + (NeutralZ[n] - 65.)*(NeutralZ[n] - 65.);
      if (R2 < 0.)
	continue;
      double R = TMath::Sqrt(R2);
      if (R == 0.)
	continue;
      double f = NeutralE[n]/R;
      TLorentzVector gamma(NeutralX[n]*f, NeutralY[n]*f, (NeutralZ[n]-65.)*f, NeutralE[n]);
      double TargetTime = NeutralT[n] - R/29.97; // hit time at target center	
      Gammas.push_back(gamma);
      GammaTimes.push_back(TargetTime);
      IdxList.push_back(n);
    }
  }

  return (int)Gammas.size();

}


double GetPhiAngle(TLorentzVector P){

  double phig = 99999.;

  if (P.X() == 0.){
    //cout<<"Error Px is zero"<<endl;
    phig = 3.1415926/2.;
    if (P.Y()<0.){
      phig += 3.1415926;
    }
  }
  
  double rphig = P.Y() / P.X();
  phig = TMath::ATan(rphig);
  if (P.X()>0){
    if (P.Y()<0){
      phig += 3.1415926*2.;
    } 
  } else {
    if (P.Y()>0){
      phig += 3.1415926;
    } else if (P.Y()<0){
      phig += 3.1415926;
    }
  }
  
  return phig;

}

double GetThetaAngle(TLorentzVector P){

  if (P.Z() == 0.){
    //cout<<"Error Pz is zero"<<endl;
    return TMath::Pi()/2.;
  }

  double rtheta =  TMath::Sqrt(P.X()*P.X() + P.Y()*P.Y()) / P.Z();
  double theta = TMath::ATan(rtheta);
  
  return theta;
}


int FindBestBeamPhotonMatch(vector <TLorentzVector> &Beam, TLorentzVector &P, TLorentzVector &P3, TH1D *missingM, TH1D *missingP, TH1D *missingE, TH2D* PvM){

  
  double dmass = 999999.;
  int bestIDX=0;
  int Size = Beam.size();

  for (int s=0; s<Size; s++) {	      
    TLorentzVector BmP = Beam[s] + ProtonTarget - P - P3;
    missingM->Fill(BmP.M());
    missingP->Fill(BmP.P());
    missingE->Fill(BmP.E());    
    PvM->Fill(BmP.M(),BmP.P());

    double dM = TMath::Abs(BmP.M() - ProtonMass);
    if (dM<dmass){
      dmass = dM;
      bestIDX = s;
    }
  } // end loop over all beam photons
  
  return bestIDX;
}







TVector3 NormalizeVector(TVector3 v){

  double mag = v.Mag();
  TVector3 vn = TVector3(v.X()/mag, v.Y()/mag, v.Z()/mag);
  return vn;

}