#include <dirent.h>
#include <TLorentzVector.h>
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];

#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 *);

void pi0anal1(){

  //char fnam[126]="pi0test_sum.root";
  //char fnam[126]="localdir/run30993/pi0test_run30993.root";
  //char fnam[126]="localdir/run31000/pi0test_run31000.root";
  //char fnam[126]="./bigfile.root";
  //char fnam[126]="pi0test_sum.root";

  TChain * myChain(0);
  myChain = new TChain("t3");

  DIR *dir;
  struct dirent *ent;
  TString directory = "./";
  
  vector <TString> FileList;
  if ((dir = opendir (directory)) != 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.BeginsWith("01")){
	FileList.push_back(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);

  TH1D *beamphotontime = new TH1D("beamphotontime", "Beam photon time minus RF-time", 1000, -25., 25.);
  TH1D *goodbeamphotons = new TH1D("goodbeamphotons", "Number of good beam photons", 10, 0., 10.);
  TH1D *goodbeamphotonsEcut = new TH1D("goodbeamphotonsEcut", "Number of good beam photons with E>8GeV", 10, 0., 10.);
  TH2D *sce = new TH2D("sce", "Start Counter Paddle vs SC-Energy NoCuts", 100, 0., 0.015, 32, 0., 32.);
  TH2D *sct = new TH2D("sct", "Start Counter Paddle vs SC-Time NoCuts", 3600, -400., 400., 32, 0., 32.);
  TH2D *scewc = new TH2D("scewc", "Start Counter Paddle vs SC-Energy With T Cut", 100, 0., 0.015, 32, 0., 32.);
  TH2D *sctwc = new TH2D("sctwc", "Start Counter Paddle vs SC-Time With E Cut", 3600, -400., 400., 32, 0., 32.);
  TH2D *sctwcTOK = new TH2D("sctwcTOK", "Start Counter Paddle vs SC-Time With E Cut", 3600, -400., 400., 32, 0., 32.);
  TH2D *sctwcTOK3 = new TH2D("sctwcTOK3", "Start Counter Paddle vs SC-Time With E Cut 3photons", 3600, -400., 400., 32, 0., 32.);
  TH1D *gammatimediff = new TH1D("gammatimediff", "2Gamma Times difference at target", 200, -10., 10.);
  TH1D *missingMass = new TH1D("missingMass", "BEST Missing Mass", 200, 10., 10.);
  TH1D *bestMM = new TH1D("bestMM", "Missing Mass", 200, -5., 5.);
  TH1D *missingP = new TH1D("missingP", "Missing Momentum", 200, 0., 10.);
  TH1D *missingEnergy = new TH1D("missingEnergy", "Missing Energy", 200, -10., 10.);
  TH1D *goodSChitsTOK = new TH1D("goodSChitsTOK", "Good SC hits in time with 2photons", 10, 0., 10.);
  TH1D *goodSChitsTOK3 = new TH1D("goodSChitsTOK3", "Good SC hits in time with 3photons", 10, 0., 10.);

  TH2D *mmPvsmmM = new TH2D("mmPvsmmM", "Missing P vs Missing M", 200, -5., 5., 200, 0., 10.);
  TH2D *bestmmPvsmmM = new TH2D("bestmmPvsmmM", "Missing P vs Missing M", 200, -5., 5., 200, 0., 10.);


  TH1D *missingMass3 = new TH1D("missingMass3", "Missing Mass three photon system", 200, 10., 10.);
  TH1D *bestMM3 = new TH1D("bestMM3", "BEST Missing Mass three photon system", 200, -5., 5.);
  TH1D *missingP3 = new TH1D("missingP3", "Missing Momentum three photon system", 200, 0., 10.);
  TH1D *missingEnergy3 = new TH1D("missingEnergy3", "Missing Energy three photon system", 200, -10., 20.);
  TH2D *mmPvsmmM3 = new TH2D("mmPvsmmM3", "Missing P vs Missing M 3gamma", 200, -5., 5., 200, 0., 10.);
  TH2D *bestmmPvsmmM3 = new TH2D("bestmmPvsmmM3", "Missing P vs Missing M 3gamma", 200, -5., 5., 200, 0., 10.);

  TH1D *twophotonInvariantMass = new TH1D("twophotonInvariantMass","Invariant Mass of two photon system", 900, 0., 3.5);
  TH1D *twophotonInvariantMassTOK = new TH1D("twophotonInvariantMassTOK","Invariant Mass of two photon system", 900, 0., 3.5);
  TH1D *twophotonInvariantMassTOKphi = new TH1D("twophotonInvariantMassTOKphi","Invariant Mass of two photon system TOK and Phi", 900, 0., 3.5);
  TH1D *twophotonInvariantMassTOKphiMM = new TH1D("twophotonInvariantMassTOKphiMM","Invariant Mass of two photon system TOK and Phi, MM<2.", 900, 0., 3.5);

  TH2D *twophotonPHIvsIMTOKphiMM = new TH2D("twophotonPHIvsIMTOKphiMM","Phi vs Invariant Mass of two photon system TOK and Phi, MM<2.", 
					    900, 0., 3.5, 30, 0., 3.1415926*2.);

  TH1D *threephotonInvariantMass = new TH1D("threephotonInvariantMass","Invariant Mass of three photon system TOK", 900, 0., 3.5);
  TH1D *threephotonInvariantMassTOK = new TH1D("threephotonInvariantMassTOK","Invariant Mass of three photon system TOK", 900, 0., 3.5);
  TH1D *threephotonInvariantMassTOKphi = new TH1D("threephotonInvariantMassTOKphi","Invariant Mass of three photon system TOK and Phi", 900, 0., 3.5);
  TH1D *threephotonInvariantMassTOKphiMM = new TH1D("threephotonInvariantMassTOKphiMM","Invariant Mass of three photon system TOK and Phi, MM<2.", 900, 0., 3.5);

  TH2D *threephotonPHIvsIMTOKphiMM = new TH2D("threephotonPHIvsIMTOKphiMM","Phi vs Invariant Mass of three photon system TOK and Phi, MM<2.", 
					      900, 0., 3.5, 30, 0., 3.1415926*2.);

  TH2D *scvsim = new TH2D("scvsim", "Start Counter Segment vs twophoton invariant mass", 200, 0., 3., 32, 0., 32.);
  TH2D *sctvsim = new TH2D("sctvsim", "Start Counter time vs twophoton invariant mass", 200, 0., 3., 400, -50., 50.);
  TH2D *scevsim = new TH2D("scevsim", "Start Counter energy vs twophoton invariant mass", 200, 0., 3., 100, 0., 0.015);
  TH2D *scevsimTOK = new TH2D("scevsimTOK", "Start Counter energy vs twophoton invariant mass", 200, 0., 3., 100, 0., 0.015);
  TH1D *screltime = new TH1D("screltime", "Start Counter relative time to photon time", 2000., -100., 100.);
  TH2D *sctvsgammat = new TH2D("sctvsgammat", "Start Counter time vs good photon time", 500, -50., 50., 400, -50., 50.);
  TH2D *sctvsgammatTOK = new TH2D("sctvsgammatTOK", "Start Counter time vs good photon time", 500, -50., 50., 400, -50., 50.);

  TH2D *scvsim3 = new TH2D("scvsim3", "Start Counter Segment vs threephoton invariant mass", 200, 0., 3., 32, 0., 32.);
  TH2D *sctvsim3 = new TH2D("sctvsim3", "Start Counter time vs threephoton invariant mass", 200, 0., 3., 400, -50., 50.);
  TH2D *scevsim3 = new TH2D("scevsim3", "Start Counter energy vs threephoton invariant mass", 200, 0., 3., 100, 0., 0.015);
  TH2D *scevsimTOK3 = new TH2D("scevsimTOK3", "Start Counter energy vs threephoton invariant mass", 200, 0., 3., 100, 0., 0.015);
  TH1D *screltime3 = new TH1D("screltime3", "Start Counter relative time to photon time 3gamma events", 2000., -100., 100.);

  TH1D *beamphotons = new TH1D("beamphotons", "Beam Photon Candidates", 400., 0., 12.);
  TH1D *bestbeamphoton = new TH1D("bestbeamphoton", "Best Beam Photon Candidate from Invariant mass", 400., 0., 12.);
  TH1D *twophotonphi = new TH1D("twophotonphi","Phi of the two photon system", 60, 0., 3.1415926*2.);
  TH1D *twophotontheta = new TH1D("twophotontheta","Theta of the two photon system", 60, 0., 1.5);
  TH1D *missingpphi = new TH1D("missingpphi","Phi of missing proton", 60, 0., 3.1415926*2.);
  TH1D *deltaPphi = new TH1D("deltaPphi","Phi of missing proton minus phi SC", 60, -3.1415926*2., 3.1415926*2.);

  TH1D *bestbeamphoton3 = new TH1D("bestbeamphoton3", "Best Beam Photon Candidate from 3 photon Invariant mass", 400., 0., 12.);
  TH1D *threephotonphi = new TH1D("threephotonphi","Phi of the two photon system", 60, 0., 3.1415926*2.);
  TH1D *threephotontheta = new TH1D("threephotontheta","Theta of the two photon system", 60, 0., 1.5);
  TH1D *missingpphi3 = new TH1D("missingpphi3","Phi of missing proton", 60, 0., 3.1415926*2.);
  TH1D *deltaPphi3 = new TH1D("deltaPphi3","Phi of missing proton minus phi SC", 60, -3.1415926*2., 3.1415926*2.);

  TH2D *twogvstwog = new TH2D("twogvstwog","#gamma #gamma M2 vs #gamma #gamma M2", 2000, 0., 2.5, 2000, 0., 2.5);
  TH1D *screltime4gammas = new TH1D("screltime4gammas", "Start Counter relative time to 4gammas ", 2000., -100., 100.);
  TH1D *fourphotonInvariantMass = new TH1D("fourphotonInvariantMass","Invariant Mass of four photon system", 900, 0., 3.5);

  TH1D *missingMass4 = new TH1D("missingMass4", "Missing Mass four photon system", 200, 10., 10.);
  TH1D *bestMM4 = new TH1D("bestMM4", "BEST Missing Mass four photon system", 200, -5., 5.);
  TH1D *missingP4 = new TH1D("missingP4", "Missing Momentum four photon system", 200, 0., 10.);
  TH1D *missingEnergy4 = new TH1D("missingEnergy4", "Missing Energy four photon system", 200, -10., 10.);
  TH1D *deltaPphi4 = new TH1D("deltaPphi4","Phi of missing proton minus phi SC", 60, -3.1415926*2., 3.1415926*2.);

  TH2D *mmPvsmmM4 = new TH2D("mmPvsmmM4", "Missing P vs Missing M 4gamma", 200, -5., 5., 200, 0., 10.);
  TH2D *bestmmPvsmmM4 = new TH2D("bestmmPvsmmM4", "Missing P vs Missing M 4gamma", 200, -5., 5., 200, 0., 10.);


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

    // loop over beam photons find intime photons
    vector <TLorentzVector> GoodBeamPhotons;
    int NBEcut = 0;
    for  (int n=0; n<Nbeam; n++){
      float dt = BeamT[n] - RFT;
      beamphotontime->Fill(dt);
      if (TMath::Abs(dt)<1.){
	if (BeamE[n]>8.){
	  TLorentzVector lv(0.,0.,BeamE[n],BeamE[n]);
	  GoodBeamPhotons.push_back(lv);
	  beamphotons->Fill(BeamE[n]);
	  NBEcut++;
	  
	}
      }
    }


    if (GoodBeamPhotons.size()<1){
      continue;
    }

    goodbeamphotons->Fill((float)GoodBeamPhotons.size());
    goodbeamphotonsEcut->Fill((float)NBEcut);
    int NGBP = GoodBeamPhotons.size();


    // Find START COUNTER Hits with energy larget than 0.5MeV AND resaonalble timing
    int FoundSCE = 0;
    int FoundSCT = 0;
    vector <int> SCIDX;
    for (int n=0; n<Nstart; n++){
      sct->Fill(SC_T[n],SC_S[n]);
      sce->Fill(SC_E[n],SC_S[n]);
      if (TMath::Abs(SC_T[n])<50.) {
	FoundSCT = 1;
	scewc->Fill(SC_E[n],SC_S[n]);
      }
      if (SC_E[n] > 0.0005) {
	FoundSCE = 1;
	sctwc->Fill(SC_T[n],SC_S[n]);
      }
      if ((TMath::Abs(SC_T[n])<50.) && (SC_E[n] > 0.0005)){
	SCIDX.push_back(n);
      }
    }
    if ((!FoundSCE) || (!FoundSCT)){
      continue;
    }

    
    // find photons with energy above 0.15GeV
    vector <TLorentzVector> Gammas;
    vector <double> GammaTimes;
    vector <int> IdxList;
    int NP = findphotons(Gammas, GammaTimes, 0.15, IdxList);

    // require at least two good photons in the detector
    if (NP<2){
      continue;
    }
    
    // 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	
	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];
	    twophotonInvariantMass->Fill(Ptot.M());
	    
	    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++){
	      scvsim->Fill(Ptot.M(),SC_S[SCIDX[l]]);
	      sctvsim->Fill(Ptot.M(),SC_T[SCIDX[l]]);
	      scevsim->Fill(Ptot.M(),SC_E[SCIDX[l]]);
	      double reltime = (GammaTimes[n]+GammaTimes[j])/2.  - SC_T[SCIDX[l]];
	      screltime->Fill(reltime);
	      sctvsgammat->Fill((GammaTimes[n]+GammaTimes[j])/2. , SC_T[SCIDX[l]]);
	      if ( (reltime>-7.) && (reltime<-1)){
		nhits++;
		TOK = 1;
		scE = SC_E[SCIDX[l]];
		scevsimTOK->Fill(Ptot.M(),SC_E[SCIDX[l]]);
		sctwcTOK->Fill(SC_T[SCIDX[l]],SC_S[SCIDX[l]]);
		sctvsgammatTOK->Fill((GammaTimes[n]+GammaTimes[j])/2. , SC_T[SCIDX[l]]);
		SC_phi = SC_S[SCIDX[l]] * 3.1414926/16. ; 		
	      }
	    }
	    if (TOK){
	      goodSChitsTOK->Fill(nhits);
	    }
	    if ((TOK) && (nhits == 1) && (scE>0.004)) {
	      twophotonInvariantMassTOK->Fill(Ptot.M());
	      goodSChitsTOK->Fill(nhits);
	      
	      double phi = GetPhiAngle(Ptot);
	      twophotonphi->Fill(phi);
	      
	      double theta = GetThetaAngle(Ptot);
	      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];
	      }
	      int IDX = FindBestBeamPhotonMatch(GoodBeamPhotons, Ptot, P3, missingMass, missingP, missingEnergy, mmPvsmmM);
	      bestbeamphoton->Fill(GoodBeamPhotons[IDX].E());
	      double BeamE = GoodBeamPhotons[IDX].E();
	      TLorentzVector BestBmP = GoodBeamPhotons[IDX] + ProtonTarget - Ptot - P3;
	      bestMM->Fill(BestBmP.M());
	      bestmmPvsmmM->Fill(BestBmP.M(), BestBmP.P());

	      // the two photons making the pi
	      
	      



	      double Bestphi = GetPhiAngle(BestBmP);
	      missingpphi->Fill(Bestphi);
	      
	      double deltaphi = Bestphi - SC_phi;
	      deltaPphi->Fill(deltaphi);
	      if (TMath::Abs(deltaphi)<1.) {
		twophotonInvariantMassTOKphi->Fill(Ptot.M());
		if ((BestBmP.M()<2.) && (BestBmP.M()>0.) && (BeamE>8.)) {
		  twophotonInvariantMassTOKphiMM->Fill(Ptot.M());
		  twophotonPHIvsIMTOKphiMM->Fill(Ptot.M(),phi);
		}
	      }
	      
	    }
	    
	    // if there are more than two photons test if these two photons match the pi0 mass
	    // if yes then look at the third photon
	    if (NP>2) {
	      if (TMath::Abs(Ptot.M()-0.133) < 0.026){
		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];
		  threephotonInvariantMass->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++){
		    scvsim3->Fill(Ptot3.M(),SC_S[SCIDX[l]]);
		    sctvsim3->Fill(Ptot3.M(),SC_T[SCIDX[l]]);
		    scevsim3->Fill(Ptot3.M(),SC_E[SCIDX[l]]);
		    double reltime3 =  (GammaTimes[0] + GammaTimes[1] + GammaTimes[2])/3. - SC_T[SCIDX[l]];
		    screltime3->Fill(reltime3);
		    if ( (reltime3>-6.) && (reltime3<-1.5)){
		      nhits3++;
		      TOK3 = 1;
		      scE = SC_E[SCIDX[l]];
		      scevsimTOK3->Fill(Ptot3.M(),SC_E[SCIDX[l]]);
		      sctwcTOK3->Fill(SC_T[SCIDX[l]],SC_S[SCIDX[l]]);
		      SC_phi3 = SC_S[SCIDX[l]] * 3.1414926/16. ;
		    }
		  }
		  if (TOK3){
		    goodSChitsTOK3->Fill(nhits3);
		  }
		  if ((TOK3) && (nhits3==1) && (scE>0.004)) {
		    threephotonInvariantMassTOK->Fill(Ptot3.M());
		    double phi = GetPhiAngle(Ptot3);
		    threephotonphi->Fill(phi);
		    double theta3 = GetThetaAngle(Ptot3);
		    threephotontheta->Fill(theta3);
		    
		    TLorentzVector P3(0.,0.,0.,0.);
		    int IDX = FindBestBeamPhotonMatch(GoodBeamPhotons, Ptot3, P3, missingMass3, missingP3, missingEnergy3,
						      mmPvsmmM3);
		    bestbeamphoton3->Fill(GoodBeamPhotons[IDX].E());
		    double BeamE = GoodBeamPhotons[IDX].E();

		    TLorentzVector BestBmP3 = GoodBeamPhotons[IDX] + ProtonTarget - Ptot3;
		    bestMM3->Fill(BestBmP3.M());
		    bestmmPvsmmM3->Fill(BestBmP3.M(), BestBmP3.P());
	      
		    double Bestphi = GetPhiAngle(BestBmP3);
		    missingpphi3->Fill(Bestphi);
		    
		    double deltaphi3 = Bestphi - SC_phi3;
		    deltaPphi3->Fill(deltaphi3);
		    if (TMath::Abs(deltaphi3)<1.) {
		      threephotonInvariantMassTOKphi->Fill(Ptot3.M());

		      if ((BestBmP3.M()<2) && (BestBmP3.M()>0.) && (BeamE>8.) ) {
			threephotonInvariantMassTOKphiMM->Fill(Ptot3.M());
			threephotonPHIvsIMTOKphiMM->Fill(Ptot3.M(), phi);
		      }

		    }
		    
		  }
		  
		}
		
	      }// 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;
	for (unsigned int l=0; l<Nstart; l++){
	  if (SC_E[l]>0.003) {
	    double reltime = meandt  - SC_T[l];
	    screltime4gammas->Fill(reltime);
	    if ( (reltime>-6.5) && (reltime<-1.5)){
	      nhits++;
	      TOK4 = 1;
	      scE = SC_E[l];
	      SC_phi4 = SC_S[l] * 3.1414926/16. ; 		
	    }
	  }
	}


	if (nhits==1){
	  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];
		    twogvstwog->Fill(Ptot1.M(),Ptot2.M());

		    if ( (TMath::Abs(Ptot1.M2()-0.018)<0.008) && (TMath::Abs(Ptot2.M2()-0.018)<0.008) ){
		      TLorentzVector Ptot4 = Ptot1 + Ptot2;

		      TLorentzVector P3(0.,0.,0.,0.);
		      if (NP>4){
			int idx = (1<<n) + (1<<j) + (1<<s) + (1<<i);
			int idx1 = -1;
			for (int l=0; l<5; l++){
			  if (((idx>>l) & 0x1) == 0){
			    idx1 = l;
			    break;
			  }
			}
			if (idx1>0){
			  P3 = Gammas[idx1];
			}
		      }
		      int IDX = FindBestBeamPhotonMatch(GoodBeamPhotons, Ptot4, P3, missingMass4, missingP4, missingEnergy4,
							mmPvsmmM4);
		      double BeamE = GoodBeamPhotons[IDX].E();
		      TLorentzVector BestBmP4 = GoodBeamPhotons[IDX] + ProtonTarget - Ptot4 - P3;
		      bestMM4->Fill(BestBmP4.M());
		      bestmmPvsmmM4->Fill(BestBmP4.M(), BestBmP4.P());

		      double Bestphi4 = GetPhiAngle(BestBmP4);
		      
		      double deltaphi4 = Bestphi4 - SC_phi4;
		      deltaPphi4->Fill(deltaphi4);
		      if ((BestBmP4.M()<1.5) && (BestBmP4.M()>0.) && (BestBmP4.P()<1.) && (TMath::Abs(deltaphi4)<1.)) {
			fourphotonInvariantMass->Fill(Ptot4.M());
		      }
		    }

		  }
		}
	      }	      
	    }
	  }
	}
      } 
    }
  }
  
  TFile *outf = new TFile("newroot.root","RECREATE");

  beamphotons->Write();
  goodbeamphotons->Write();
  goodbeamphotonsEcut->Write();
  sce->Write();
  sct->Write();
  scewc->Write();
  sctwc->Write();
  sctwcTOK->Write();
  gammatimediff->Write();
  missingMass->Write();
  bestMM->Write();
  missingP->Write(); 
  missingEnergy->Write();
  mmPvsmmM->Write();
  bestmmPvsmmM->Write();

  missingMass3->Write();
  bestMM3->Write();
  missingP3->Write(); 
  missingEnergy3->Write();
  mmPvsmmM3->Write();
  bestmmPvsmmM3->Write();

  twophotonInvariantMass->Write();
  scvsim->Write();
  sctvsim->Write();
  scevsim->Write();
  screltime->Write();
  scevsimTOK->Write();
  sctvsgammat->Write();
  sctvsgammatTOK->Write();
  bestbeamphoton->Write();
  twophotonphi->Write();
  twophotontheta->Write();
  missingpphi->Write();
  deltaPphi->Write();
  goodSChitsTOK->Write();

  bestbeamphoton3->Write();
  threephotonphi->Write();
  threephotontheta->Write();
  missingpphi3->Write();
  deltaPphi3->Write();


  threephotonInvariantMass->Write();
  threephotonInvariantMassTOK->Write();
  twophotonInvariantMassTOKphi->Write();
  twophotonInvariantMassTOKphiMM->Write();
  twophotonPHIvsIMTOKphiMM->Write();

  threephotonInvariantMassTOKphi->Write();
  threephotonInvariantMassTOKphiMM->Write();
  threephotonPHIvsIMTOKphiMM->Write();
 
  sctwcTOK3->Write();
  scvsim3->Write();
  sctvsim3->Write();
  scevsim3->Write();
  screltime3->Write();
  scevsimTOK3->Write();
  goodSChitsTOK3->Write();


  twogvstwog->Write();
  screltime4gammas->Write();
  fourphotonInvariantMass->Write();

  missingMass4->Write();
  bestMM4->Write();
  missingP4->Write(); 
  missingEnergy4->Write();
  deltaPphi4->Write();
  mmPvsmmM4->Write();
  bestmmPvsmmM4->Write();

  outf->Close();

}



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

  int cnt = 0;
  for  (int n=0; n<Nneut; n++){
    if (NeutralE[n]>EnergyCut) {
      double R2 = NeutralX[n]*NeutralX[n] + NeutralY[n]*NeutralY[n] + (NeutralZ[n] - 65.)*(NeutralZ[n] - 65.);
      double R = TMath::Sqrt(R2);
      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 rphig = P.Y() / P.X();
  double 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){

  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;
}