void pi0anal(){

  //char fnam[126]="pi0test_sum.root";
  char fnam[126]="localdir/run30992/pi0test_run30992_full.root";

  TFile *RF = new TFile(fnam,"READ");

  TTree *t3 = (TTree*)RF->Get("t3");
  
  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];

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

  TH1D *goodbeamphotons = new TH1D("goodbeamphotons", "Number of good beam photons", 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", "Missing Mass", 200, 10., 10.);
  TH1D *missingP = new TH1D("missingP", "Missing Momentum", 200, 0., 10.);
  TH1D *missingEnergy = new TH1D("missingEnergy", "Missing Mass", 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.);


  TH1D *missingMass3 = new TH1D("missingMass3", "Missing Mass three photon system", 200, 10., 10.);
  TH1D *missingP3 = new TH1D("missingP3", "Missing Momentum three photon system", 200, 0., 10.);
  TH1D *missingEnergy3 = new TH1D("missingEnergy3", "Missing Mass three photon system", 200, -10., 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 *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);

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

  TLorentzVector ProtonTarget(0.,0.,0., 0.938272); // Proton Target 4-vector

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

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

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

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

    int NGBP = GoodBeamPhotons.size();

    if (NGBP<1) {
      continue;
    }
    
    vector <TLorentzVector> Gammas;
    vector <double> GammaTimes;
    for  (int n=0; n<Nneut; n++){
      if (NeutralE[n]>0.15) {
	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);
      }
    }

    int NP = Gammas.size();
    double SC_phi = 99999.;
    TLorentzVector Ptot;
    TLorentzVector Ptot3;
    int TOK = 0;
    int TOK3 = 0;
    if (NP <4) {
      for  (int n=0; n<NP-1; n++){
	
	for  (int j=n+1; j<NP; j++){
	  
	  double dt = GammaTimes[n]-GammaTimes[j];
	  gammatimediff->Fill(dt);
	  if (TMath::Abs(dt)<1.5) {
	    Ptot = Gammas[n] + Gammas[j];
	    twophotonInvariantMass->Fill(Ptot.M());

	    TOK = 0;
	    int nhits=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] - SC_T[SCIDX[l]];
	      screltime->Fill(reltime);
	      if ( (reltime>-7.) && (reltime<-1)){
		nhits++;
		TOK = 1;
		scevsimTOK->Fill(Ptot.M(),SC_E[SCIDX[l]]);
		sctwcTOK->Fill(SC_T[SCIDX[l]],SC_S[SCIDX[l]]);
		SC_phi = SC_S[SCIDX[l]] * 3.1414926/16. ; 		
	      }
	    }
	    if (TOK){
	      twophotonInvariantMassTOK->Fill(Ptot.M());
	      goodSChitsTOK->Fill(nhits);
	    }

	    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;
		}
		
		double dt3 =  GammaTimes[n] - GammaTimes[i3];
		if (TMath::Abs(dt3)<1.5) {
		  Ptot3 = Gammas[n] + Gammas[j] + Gammas[i3];
		  threephotonInvariantMass->Fill(Ptot3.M());

		  TOK3 = 0;
		  int nhits3=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;
		      scevsimTOK3->Fill(Ptot3.M(),SC_E[SCIDX[l]]);
		      sctwcTOK3->Fill(SC_T[SCIDX[l]],SC_S[SCIDX[l]]);
		    }
		  }
		  if (TOK3){
		    threephotonInvariantMassTOK->Fill(Ptot3.M());
		    goodSChitsTOK3->Fill(nhits);
		  }
		}
		
	      }
	    }

	    //if ((TOK) && (TMath::Abs(Ptot.M()-0.133) < 0.026)){
	    if ((TOK)){
	      
	      double rphig = Ptot.Y() / Ptot.X();
	      double phig = TMath::ATan(rphig);
	      if (Ptot.X()>0){
		if (Ptot.Y()<0){
		  phig += 3.1415926*2.;
		} 
	      } else {
		if (Ptot.Y()>0){
		  phig += 3.1415926;
		} else if (Ptot.Y()<0){
		  phig += 3.1415926;
		}
	      }
	      twophotonphi->Fill(phig);


	      double rtheta =  TMath::Sqrt(Ptot.X()*Ptot.X() + Ptot.Y()*Ptot.Y()) / Ptot.Z();
	      double theta = TMath::ATan(rtheta);
	      twophotontheta->Fill(theta);

	      // ok got 2photon particle
	      // compare with beam photon and find closed to missing proton particle
	      // find photon with closest Missing Mass to 0.938GeV
	      double dmass = 999999.;
	      int bestIDX=0;
	      for (int s=0; s<NGBP; s++) {	      
		TLorentzVector BmP = GoodBeamPhotons[s] + ProtonTarget - Ptot;
		missingMass->Fill(BmP.M());
		missingP->Fill(BmP.P());
		missingEnergy->Fill(BmP.E());
		
		double dM = TMath::Abs(BmP.M()-0.938);
		if (dM<dmass){
		  dmass = dM;
		  bestIDX = s;
		}
	      } // end loop over all beam photons
	      bestbeamphoton->Fill(GoodBeamPhotons[bestIDX].E());
	      TLorentzVector BestBmP = GoodBeamPhotons[bestIDX] + ProtonTarget - Ptot;
	      double rphi =  BestBmP.Y() / BestBmP.X();
	      double phi = TMath::ATan(rphi);
	      if (BestBmP.X()>0){
		if (BestBmP.Y()<0){
		  phi += 3.1415926*2.;
		} 
	      } else {
		if (BestBmP.Y()>0){
		  phi += 3.1415926;
		} else if (BestBmP.Y()<0){
		  phi += 3.1415926;
		}
	      }
	      missingpphi->Fill(phi);

	      double deltaphi = phi - SC_phi;
	      deltaPphi->Fill(deltaphi);
	      if (TMath::Abs(deltaphi)<1.) {
		twophotonInvariantMassTOKphi->Fill(Ptot.M());
	      }



	      // 3 photons
	      if (TOK3) {
		rphig = Ptot3.Y() / Ptot3.X();
		phig = TMath::ATan(rphig);
		if (Ptot3.X()>0){
		  if (Ptot3.Y()<0){
		    phig += 3.1415926*2.;
		  } 
		} else {
		  phig += 3.1415926;
		}
		threephotonphi->Fill(phig);

		double rtheta3 =  TMath::Sqrt( Ptot3.X()*Ptot3.X() +  Ptot3.Y()*Ptot3.Y()) / Ptot3.Z();
		double theta3 = TMath::ATan(rtheta3);
		threephotontheta->Fill(theta3);
		dmass = 999999.;
		bestIDX=0;
		for (int s=0; s<NGBP; s++) {	      
		  TLorentzVector BmP3 = GoodBeamPhotons[s] + ProtonTarget - Ptot3;
		  missingMass3->Fill(BmP3.M());
		  missingP3->Fill(BmP3.P());
		  missingEnergy3->Fill(BmP3.E());
		  
		  double dM = TMath::Abs(BmP3.M()-0.938);
		  if (dM<dmass){
		    dmass = dM;
		    bestIDX = s;
		  }
		} // end loop over all beam photons
		bestbeamphoton3->Fill(GoodBeamPhotons[bestIDX].E());
		TLorentzVector BestBmP3 = GoodBeamPhotons[bestIDX] + ProtonTarget - Ptot3;
		rphi =  BestBmP3.Y() / BestBmP3.X();
		phi = TMath::ATan(rphi);
		if (BestBmP3.X()>0){
		  if (BestBmP3.Y()<0){
		    phi += 3.1415926*2.;
		  } 
		} else {
		  phi += 3.1415926;
		}
		missingpphi3->Fill(phi);
		
		double deltaphi3 = phi - SC_phi;
		deltaPphi3->Fill(deltaphi3);
		if (TMath::Abs(deltaphi3)<1.) {
		  threephotonInvariantMassTOKphi->Fill(Ptot3.M());
		}
	      }
	    }
	  }
	}
      }
    } // end loop over all good gammas
  }


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

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

  missingMass3->Write();
  missingP3->Write(); 
  missingEnergy3->Write();

  twophotonInvariantMass->Write();
  scvsim->Write();
  sctvsim->Write();
  scevsim->Write();
  screltime->Write();
  scevsimTOK->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();
  threephotonInvariantMassTOKphi->Write();

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

  outf->Close();

}