#define TwoG_jan_8_2021_cxx

// Author: J. McIntyre
// March 4, 2021
// This file contains cuts on real data for a 2 photon sample.
////////////////////////////////////////////////////////////////////////
// This was written to look at in the 2 gamma sample on Mar. 4, 2021. //
// __________________________________________________________________ //
// __________________________________________________________________ //
//                                                                    //
// THIS FILE HAS UPDATES THAT:                                        //
// --------------------------                                         //
//    - USES THE LATEST ENERGY RESCALING FACTOR (ER) = 0.97          //
// __________________________________________________________________ //
//                                                                    //
// This file is a re-do of ER_real_pi02gp_oct_21_2021.C,              //
// /nfs/direct/annex/mcintyre/analysisIII/ER/ER_0_955/5MeV_mass_bins  //
// since the ER & ASF were updated and a valid cross section for pi0  //
// was plotted. Now this file is the start of a final analysis of the //
// two gamma sample for the pi0, eta, and etap cross sections.        //
//                                                                    //
// The new ER equals 0.97, as determined by the various attempts in  //
// folder: /nfs/direct/annex/mcintyre/analysisIX/TwoGamma/inv_mass    //
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
// To modify this file start by update the file name listed below.    //
// This file's name --> TwoG_jan_8_2021                               //
////////////////////////////////////////////////////////////////////////

/*
   The following methods are defined in this file:
      Begin():        called every time a loop on the tree starts,
                      a convenient place to create your histograms.
      SlaveBegin():   called after Begin(), when on PROOF called only on the
                      slave servers.
      Process():      called for each event, in this function you decide what
                      to read and fill your histograms.
      SlaveTerminate: called at the end of the loop on the tree, when on PROOF
                      called only on the slave servers.
      Terminate():    called at the end of the loop on the tree,
                      a convenient place to draw/fit your histograms.
*/

#include <iostream>
#include "TwoG_jan_8_2021.h"
#include <TH2.h>
#include <TStyle.h>
#include <math.h>
#include <cmath>
#include <TLorentzVector.h>
#include <TVector3.h>
#include <TMath.h>
#include <TProofServ.h>
#include <TF1.h>
#include <vector>


#define sqr(A) ((A)*(A))     // Define a function to square a number
#define PI 3.14159265359
#define fwdgamma 2           // Defines the number of photons in 
                             // the LGD for the events of interest
#define bgvgamma 2           // Defines how many photons are permitted in the BGV
                             // [bgvgamma - fwdgamma] = photons_in_BGV

#define ResonanceEner 0.300  // Cut on missing energy (e.g. resonance energy of the proton)(GeV)

#define ASF 0.935198         // Accidental Scaling Factor -- **pass-9-2020 dataset** Average ASF = 0.935198 +/- 0.00249013 
                             // (*old* pass-1-2020 -> 0.935144 +/- 0.00162298) 
#define ER 0.97              // Energy Rescaler -- Shifts invariant mass plot to 
                             // match well known mesons
                             // Previous tries: 0.955, 0.953, 0.9668            

#define TYPECUT 1            // Cluster_cleanup cut

using namespace std;         // Removes the need for std::

TFile *fHistFile;            // Set Pointer to the output file


void TwoG_jan_8_2021::Begin(TTree * /*tree*/)
{
   TString option = GetOption();
}

void TwoG_jan_8_2021::SlaveBegin(TTree * /*tree*/)
{
   TString option = GetOption();

   // Histogram definition
   // Cut 1
   hcut1 = new TH1D("hcut1", ";M_{2#gamma} [GeV/c^{2}];Events per 5 MeV/c^{2}", 280, 0.0, 1.4);
   hcut1->Sumw2();
   GetOutputList()->Add(hcut1);

   // Cut 2
   hcut2 = new TH1D("hcut2", ";M_{2#gamma} [GeV/c^{2}];Events per 5 MeV/c^{2}", 280, 0.0, 1.4);
   hcut2->Sumw2();
   GetOutputList()->Add(hcut2);

   // Cut 3
   hcut3 = new TH1D("hcut3", ";M_{2#gamma} [GeV/c^{2}];Events per 5 MeV/c^{2}", 280, 0.0, 1.4);
   hcut3->Sumw2();
   GetOutputList()->Add(hcut3);

   // Cut 4
   hcut4 = new TH1D("hcut4", ";M_{2#gamma} [GeV/c^{2}];Events per 5 MeV/c^{2}", 280, 0.0, 1.4);
   hcut4->Sumw2();
   GetOutputList()->Add(hcut4);

   // Cut 5
   hcut5 = new TH1D("hcut5", ";M_{2#gamma} [GeV/c^{2}];Events per 5 MeV/c^{2}", 280, 0.0, 1.4);
   hcut5->Sumw2();
   GetOutputList()->Add(hcut5);

   // Cut 6
   hcut6 = new TH1D("hcut6", ";M_{2#gamma} [GeV/c^{2}];Events per 5 MeV/c^{2}", 280, 0.0, 1.4);
   hcut6->Sumw2();
   GetOutputList()->Add(hcut6);

   // Cut 7
   hcut7 = new TH1D("hcut7", ";M_{2#gamma} [GeV/c^{2}];Events per 5 MeV/c^{2}", 280, 0.0, 1.4);
   hcut7->Sumw2();
   GetOutputList()->Add(hcut7);

   // Cut 8
   hcut8 = new TH1D("hcut8", ";M_{2#gamma} [GeV/c^{2}];Events per 5 MeV/c^{2}", 280, 0.0, 1.4);
   hcut8->Sumw2();
   GetOutputList()->Add(hcut8);


   // Coincidence hits after cuts 1 - 8
   h_ncoin = new TH1D("h_ncoin", ";Coincidence Hits;Number of 2#gamma Events", 124, -1.0, 30.0);
   h_ncoin->Sumw2();
   GetOutputList()->Add(h_ncoin);

   // Time difference between recoil proton and the i'th tagger hit
   // (Time of best recoil proton candidate minus tagger hit time of hit[i])
   h_dt = new TH1D("h_dt", ";t_{recoil particle} - t_{tagger hit} [ns];Events per 0.1 ns", 1900, -90.0, 100.0);
   h_dt->Sumw2();
   GetOutputList()->Add(h_dt);

   // Cuts 1-8, channel diff. > 2, Peak at -40 ns, tag_{j} btw -53 to -27 ns;cotime_{i} - cotime_{j} [ns]
   h2g_all = new TH1D("h2g_all", ";t_{accidental tag (i)} - t_{accidental tag (j)} [ns];Events per 0.1 ns", 320, -16.0, 16.0);
   h2g_all->Sumw2();
   GetOutputList()->Add(h2g_all);

   // t_{tag_{i}} - t_{tag_{j}}, Cuts 1-8, channel diff. > 2, Peak at -42 ns, tag_{j} btw -61 to -25 ns
   h_dcotime = new TH1D("h_dcotime", ";t_{accidental tag (i)} - t_{accidental tag (j)} [ns];Events per 0.1 ns", 420, -22.0, 20.0);
   h_dcotime->Sumw2();
   GetOutputList()->Add(h_dcotime);

   // t_{tag_{i}} - t_{tag_{j}}, Cuts 1-8, channel diff. > 2, Peak at -44 ns, tag_{j} btw -61 to -25 ns
   h_dcotime2 = new TH1D("h_dcotime2", ";t_{accidental tag (i)} - t_{accidental tag (j)} [ns];Events per 0.1 ns", 420, -20.0, 22.0);
   h_dcotime2->Sumw2();
   GetOutputList()->Add(h_dcotime2);

   // |t| for 2#gamma multiplicity
   t_2gamma = new TH1D("t_2gamma", ";|t| [GeV^{2}];Events per 5 MeV^{2}", 1000, 0.0, 5.0);
   t_2gamma->Sumw2();
   GetOutputList()->Add(t_2gamma);

   // M_{2#gamma} for 2#gamma final state (all |t|)
   M_2gamma = new TH1D("M_2gamma", ";M_{2#gamma} [GeV/c^{2}];Events per 5 MeV/c^{2}", 280, 0.0, 1.4);
   M_2gamma->Sumw2();
   GetOutputList()->Add(M_2gamma);

   // M_{2#gamma} for 2#gamma final state (for all #pi^{0} |t| bins)
   M_Pi = new TH1D("M_Pi", ";M_{2#gamma} [GeV/c^{2}];Events per 5 MeV/c^{2}", 280, 0.0, 1.4);
   M_Pi->Sumw2();
   GetOutputList()->Add(M_Pi);

   // M_{2#gamma} for 2#gamma final state (for all #eta |t| bins)
   M_Eta = new TH1D("M_Eta", ";M_{2#gamma} [GeV/c^{2}];Events per 5 MeV/c^{2}", 280, 0.0, 1.4);
   M_Eta->Sumw2();
   GetOutputList()->Add(M_Eta);

   // M_{2#gamma} for 2#gamma final state (for all #eta' |t| bins)
   M_Etap = new TH1D("M_Etap", ";M_{2#gamma} [GeV/c^{2}];Events per 5 MeV/c^{2}", 280, 0.0, 1.4);
   M_Etap->Sumw2();
   GetOutputList()->Add(M_Etap);


   double XPi_bins[28][2] = {{0.100, 0.125}, {0.125, 0.150}, {0.150, 0.175}, {0.175, 0.200}, {0.200, 0.240},
                             {0.240, 0.280}, {0.280, 0.320}, {0.320, 0.360}, {0.360, 0.410}, {0.410, 0.460},
                             {0.460, 0.510}, {0.510, 0.585}, {0.585, 0.660}, {0.660, 0.735}, {0.735, 0.810},
                             {0.810, 0.885}, {0.885, 0.970}, {0.970, 1.055}, {1.055, 1.155}, {1.155, 1.255},
                             {1.255, 1.375}, {1.375, 1.500}, {1.500, 1.650}, {1.650, 1.825}, {1.825, 2.050},
                             {2.050, 2.350}, {2.350, 2.750}, {2.750, 3.500}};

   double XEta_bins[18][2] = {{0.100, 0.150}, {0.150, 0.200}, {0.200, 0.250}, {0.250, 0.300}, {0.300, 0.350},
                              {0.350, 0.400}, {0.400, 0.450}, {0.450, 0.525}, {0.525, 0.600}, {0.600, 0.700},
                              {0.700, 0.800}, {0.800, 0.950}, {0.950, 1.150}, {1.150, 1.400}, {1.400, 1.650},
                              {1.650, 1.950}, {1.950, 2.350}, {2.350, 3.000}};

   double XEtap_bins[6][2] = {{0.100, 0.250}, {0.250, 0.450}, {0.450, 0.700}, {0.700, 1.000}, {1.000, 1.400},
                              {1.400, 2.000}};

   TString title;
   TString name;
   for (int i = 0; i < 28; i++) {
      name.Form("Pi_mass%d", i);
      title.Form("M_{2#gamma} for #pi^{0} fit, %.3f <= |t| < %.3f;M_{2#gamma} [MeV/c^{2}];Events per 5 MeV/c^{2}", XPi_bins[i][0], XPi_bins[i][1]);
      Pi_mass[i] = new TH1D(name, title, 280, 0.0, 1.4);
      Pi_mass[i]->Sumw2();
      GetOutputList()->Add(Pi_mass[i]);
   }

   for (int i = 0; i < 18; i++) {
      name.Form("Eta_mass%d", i);
      title.Form("M_{2#gamma} for #eta fit, %.3f <= |t| < %.3f;M_{2#gamma} [MeV/c^{2}];Events per 5 MeV/c^{2}", XEta_bins[i][0], XEta_bins[i][1]);
      Eta_mass[i] = new TH1D(name, title, 280, 0.0, 1.4);
      Eta_mass[i]->Sumw2();
      GetOutputList()->Add(Eta_mass[i]);
   }

   for (int i = 0; i < 6; i++) {
      name.Form("Etap_mass%d", i);
      title.Form("M_{2#gamma} for #eta' fit, %.3f <= |t| < %.3f;M_{2#gamma} [MeV/c^{2}];Events per 5 MeV/c^{2}", XEtap_bins[i][0], XEtap_bins[i][1]);
      Etap_mass[i] = new TH1D(name, title, 280, 0.0, 1.4);
      Etap_mass[i]->Sumw2();
      GetOutputList()->Add(Etap_mass[i]);
   }

}


Bool_t TwoG_jan_8_2021::Process(Long64_t entry)
{
   if (fChain == 0) {   // REMEMBER:
      return false;     // false is valid for ROOT 6 and above.
   }                    // kFALSE works for all ROOT, but not C++.

   GetEntry(entry);

   // Rescaling the LGD
   double EFRWD = ER * Efrwd;
   double PVECT[7][4];
   for (int ph = 0; ph < 7; ph++) {
      for (int fv = 0; fv < 4; fv++) {
         PVECT[ph][fv] = ER * pvect[ph][fv];
      }
   }


   // Setting TLorentzVectors for all photons
   TLorentzVector photon_0;
   TLorentzVector photon_1;

   photon_0.SetPxPyPzE(PVECT[0][1],PVECT[0][2],PVECT[0][3],PVECT[0][0]);
   photon_1.SetPxPyPzE(PVECT[1][1],PVECT[1][2],PVECT[1][3],PVECT[1][0]);


   //////////////////////////////////////////////////////////
   //////////// Determining the Invariant Masses ////////////
   //////////////////////////////////////////////////////////
   // 2 Gamma Invariant Mass (imass2)
   double imass2 = 0.0;
   imass2 = (photon_0 + photon_1).M();


   // Cut to look at cluster_cleanup turned on
   #if TYPECUT
   if (evtype > 10) {
      return false;
   }
   #endif

   ///////////////////////////////////////////////////////
   //////////// Performing the first few cuts ////////////
   ///////////////////////////////////////////////////////
   // Cut #1 - Number of showers reconstructed in the LGD
   if (nfrwd != fwdgamma) {
      return false;   
   }
   hcut1->Fill(imass2);

   // Cut #2 - Number of reconstructed recoil particles
   // Exactly 1 pixel cluster in the BSD
   if (nrec != 1) {
      return false;
   }
   hcut2->Fill(imass2);

   // Cut #3 - No showers unrelated to the recoil reconstructed in the BGD
   if (nphot != bgvgamma) {
      return false;
   }
   hcut3->Fill(imass2);

   // Cut #4 - Individual shower energy in LGD
   // Individual shower energy in the LGD must be at least 0.05 GeV
   int cut4 = 0;
   for (int i = 0; i < nphot; i++) {
      if ((PVECT[i][0]) < 0.050) {
         cut4 += 1;
      }
   }
   if (cut4 != 0) {
      return false;
   }
   hcut4->Fill(imass2);

   // Cut #5 - Total shower energy in the LGD
   // Total shower energy in the LGD must be at least 3.0 GeV
   // Use this if allowing a photon in BGD (e.g. nphot != nfrwd)
   // Otherwise i < nphot can be used.
   double tot_shower_e = 0.0;
   for (int i = 0; i < nfrwd; i++) {
      tot_shower_e += PVECT[i][0];
   }
   if (tot_shower_e < 3.0) {
      return false;
   }
   hcut5->Fill(imass2);

   // Cut #6 - No coincident hits in the forward veto scintillators (CPV)
   // No coincidence between a hit in the CPV & a recoil particle
   int cpv_recoil_coin = 0;
   // Cycling through the hits in the CPV
   for (int i = 0; i < nhcpv; i++) {
      // Cycling through number of times there were a CPV hit [i]
      for (int j = 0; j < ntcpv[i]; j++) {
         // Leading edge times minus recoil proton time
         double dt = le[t1cpv[i] + j] - trec0;
         if (dt >= -3.0 && dt <= 3.0) {
            cpv_recoil_coin += 1;
         }
      }
   }
   if (cpv_recoil_coin != 0) {
      return false;
   }
   hcut6->Fill(imass2);

   // Cut #7 - No low-energy showers around the beam hole in the LGD
   // Individual showers in the LGD that are < 0.5 GeV must satisfy
   // the following condition. This reduces the observed EM background
   // that peaks at low angles & low energies.
   int beam_hole_cut = 0;
   for (int i = 0; i < nphot; i++) {
      double lgd_xy = (sqrt((sqr(PVECT[i][1])) + (sqr(PVECT[i][2]))));
      double lgd_z = PVECT[i][3];
      // Polar angle in deg.
      double theta = ((atan2(lgd_xy, lgd_z))*(180/PI));
      // individual shower energy 
      double ishower_E = PVECT[i][0];
      double test_level = (1.033 - (0.13*theta));
      if ((ishower_E <= test_level) && (ishower_E < 0.5)) {
         beam_hole_cut += 1;
      }
   }
   if (beam_hole_cut != 0) {
      return false;
   }
   hcut7->Fill(imass2);

   // Cut #8 - Suppresses proton resonances (cleans up the signal)
   int summed_eloss = 0;
   double E_missing = 0.0;
   for (int i = 0; i < ncoin; i++) {
      E_missing = (coenergy[i] - EFRWD);
      if (E_missing < ResonanceEner) {
         summed_eloss += 1;
      }
   }
   if (summed_eloss == 0) {
      return false;
   }
   hcut8->Fill(imass2);

   // Histos showing the number of coincidence hits we are dealing with
   double ncoinplus = ncoin + 0.1;
   double ncoinminus = ncoin - 0.1;
   h_ncoin->Fill(ncoinplus);
   h_ncoin->Fill(ncoinminus);

   ///////////////////// CALCULATION SECTION /////////////////////
   // Summed shower components
   double ener_sum = 0.0;
   double px_sum = 0.0; 
   double py_sum = 0.0; 
   double pz_sum = 0.0; 
   for (int i = 0; i < nphot; i++) {
      ener_sum += PVECT[i][0];     
      px_sum   += PVECT[i][1];     
      py_sum   += PVECT[i][2];     
      pz_sum   += PVECT[i][3];     
   }

   // nphot Gamma Invariant Mass
   TLorentzVector decay_lab;
   decay_lab.SetPxPyPzE(px_sum, py_sum, pz_sum, ener_sum);

   /*
   // delta phi =  azimuthal angular difference between the recoil plane
   // & the Pi meson. Delta phi is used to suppress proton resonances.
   double phi_rec = 0.0;
   double rec_phi = 0.0;
   double two_phi = 0.0;
   double del_phi = 0.0;

   phi_rec = (phirec[0])*(180.0 / PI); 
   rec_phi = constrain_angle(phi_rec);   // Gives the angle - 180 deg.
   two_phi = decay_lab.Phi()*(180.0 / PI);
   del_phi = angle_diff(two_phi, rec_phi);
   */

   // Weighting factor for the fill statements
   double dt_tag2 = 0.0;
   double dt = 0.0;
   double dt3 = 0.0;
   int delta_cochan2 = 0;

   for (int j = 0; j < ncoin; j++) {
      // Time difference between recoil proton and the j'th tagger hit
      // Time of best recoil proton candidate minus tagger hit time of hit[j].
      dt = trec0 - cotime[j];
      h_dt->Fill(dt);

      if ((dt >= -41.0) && (dt <= -39.0)) {
         for (int k = 0; k < ncoin; k++) {
            delta_cochan2 = abs(cochan[j] - cochan[k]);
            if ((j != k) && (delta_cochan2 >= 3)) {
               dt_tag2 = cotime[j] - cotime[k];
               dt3 = trec0 - cotime[k];
               if ((dt3 >= -53.0) && (dt3 <= -27.0)) {
                  h2g_all->Fill(dt_tag2);
               }
            }
         }
      }
   }

   // reset variables to zero for reuse without bias
   dt_tag2 = 0.0;
   dt = 0.0;
   dt3 = 0.0;
   delta_cochan2 = 0;

   for (int l = 0; l < ncoin; l++) {
		// reset variables to zero for reuse without bias
		dt_tag2 = 0.0;
		dt = 0.0;
		dt3 = 0.0;
		delta_cochan2 = 0;

      // Time difference between recoil proton and the l'th tagger hit
      // Time of best recoil proton candidate minus tagger hit time of hit[l].
      dt = trec0 - cotime[l];

      if ((dt >= -43.0) && (dt <= -41.0)) {
         for (int m = 0; m < ncoin; m++) {
            delta_cochan2 = abs(cochan[l] - cochan[m]);
            if ((l != m) && (delta_cochan2 >= 3)) {
               dt_tag2 = cotime[l] - cotime[m];
               dt3 = trec0 - cotime[m];
               if ((dt3 >= -61.0) && (dt3 <= -25.0)) {
                  h_dcotime->Fill(dt_tag2);
               }
            }
         }
      }

		// reset variables to zero for reuse without bias
		dt_tag2 = 0.0;
		dt3 = 0.0;
		delta_cochan2 = 0;

      if ((dt >= -45.0) && (dt <= -43.0)) {
         for (int n = 0; n < ncoin; n++) {
            delta_cochan2 = abs(cochan[l] - cochan[n]);
            if ((l != n) && (delta_cochan2 >= 3)) {
               dt_tag2 = cotime[l] - cotime[n];
               dt3 = trec0 - cotime[n];
               if ((dt3 >= -61.0) && (dt3 <= -25.0)) {
                  h_dcotime2->Fill(dt_tag2);
               }
            }
         }
      }
   }

   // Filling t-bins
   double Pi_bins[28][2] = {{0.100, 0.125}, {0.125, 0.150}, {0.150, 0.175}, {0.175, 0.200}, {0.200, 0.240},
                            {0.240, 0.280}, {0.280, 0.320}, {0.320, 0.360}, {0.360, 0.410}, {0.410, 0.460},
                            {0.460, 0.510}, {0.510, 0.585}, {0.585, 0.660}, {0.660, 0.735}, {0.735, 0.810},
                            {0.810, 0.885}, {0.885, 0.970}, {0.970, 1.055}, {1.055, 1.155}, {1.155, 1.255},
                            {1.255, 1.375}, {1.375, 1.500}, {1.500, 1.650}, {1.650, 1.825}, {1.825, 2.050},
                            {2.050, 2.350}, {2.350, 2.750}, {2.750, 3.500}};

   double Eta_bins[18][2] = {{0.100, 0.150}, {0.150, 0.200}, {0.200, 0.250}, {0.250, 0.300}, {0.300, 0.350},
                             {0.350, 0.400}, {0.400, 0.450}, {0.450, 0.525}, {0.525, 0.600}, {0.600, 0.700},
                             {0.700, 0.800}, {0.800, 0.950}, {0.950, 1.150}, {1.150, 1.400}, {1.400, 1.650},
                             {1.650, 1.950}, {1.950, 2.350}, {2.350, 3.000}};

   double Etap_bins[6][2] = {{0.100, 0.250}, {0.250, 0.450}, {0.450, 0.700}, {0.700, 1.000}, {1.000, 1.400},
                             {1.400, 2.000}};

   /*
   double Eta_bins[23][2] = {{0.000, 0.050}, {0.050, 0.100}, {0.100, 0.150}, {0.150, 0.200}, {0.200, 0.250},
                             {0.250, 0.300}, {0.300, 0.350}, {0.350, 0.400}, {0.400, 0.450}, {0.450, 0.525},
                             {0.525, 0.600}, {0.600, 0.675}, {0.675, 0.750}, {0.750, 0.850}, {0.850, 0.950},
                             {0.950, 1.050}, {1.050, 1.200}, {1.200, 1.350}, {1.350, 1.550}, {1.550, 1.800},
                             {1.800, 2.100}, {2.100, 2.500}, {2.500, 3.000}};

   double Etap_bins[8][2] = {{0.000, 0.050}, {0.050, 0.150}, {0.150, 0.300}, {0.300, 0.450}, {0.450, 0.600},
                             {0.600, 0.750}, {0.750, 1.000}, {1.000, 2.000}};
   */

   for (int i = 0; i < ncoin; i++) {
      double tabs = 0.0;
      double dt = 0.0;
      double tag_wgt = 0.0;
      double Piwgt = 0.0;
      double Pidphiwgt = 0.0;

      // Time of best recoil proton candidate minus tagger hit time of hit[i].
      dt = trec0 - cotime[i];
      if (dt >= -1.0 && dt <= 5.0) {
         tag_wgt = 1.0;
      }

      else if (dt > 5.0 && dt <= 11.0) {
         tag_wgt = -(1.0 / ASF);
      }

      else {
         tag_wgt = 0.0;
      }

      tabs = -((imass2 * imass2) - (2 * coenergy[i] * ener_sum) + (2 * coenergy[i] * pz_sum));

      if (tabs > 0.1 && tabs < 3.5) {
         M_2gamma->Fill(imass2, tag_wgt);
      }

      for (int xp = 0; xp < 28; xp++) {
         if (tabs >= Pi_bins[xp][0] && tabs < Pi_bins[xp][1]) {
            Pi_mass[xp]->Fill(imass2, tag_wgt);
            M_Pi->Fill(imass2, tag_wgt);
         }
      }

      for (int xpa = 0; xpa < 18; xpa++) {
         if (tabs >= Eta_bins[xpa][0] && tabs < Eta_bins[xpa][1]) {
            Eta_mass[xpa]->Fill(imass2, tag_wgt);
            M_Eta->Fill(imass2, tag_wgt);
         }
      }

      for (int xpb = 0; xpb < 6; xpb++) {
         if (tabs >= Etap_bins[xpb][0] && tabs < Etap_bins[xpb][1]) {
            Etap_mass[xpb]->Fill(imass2, tag_wgt);
            M_Etap->Fill(imass2, tag_wgt);
         }
      }

   }

   return kTRUE;
}

void TwoG_jan_8_2021::SlaveTerminate()
{
  
}

void TwoG_jan_8_2021::Terminate()
{
   // Write histograms to the output file
   fHistFile = new TFile("TwoG_jan_8_2021Analysis.root", "RECREATE");   // Declare output file
   std::cout << "pushing result histograms to a file..." << std::endl;
   std::cout << "my outputlist has " << GetOutputList()->GetEntries() << " entries." << std::endl;
   for (TIter iter = GetOutputList()->begin(); iter != GetOutputList()->end(); ++iter) {
      (*iter)->Write();
   }
   fHistFile->Close();
   cout << "Output file has been written" << endl;
}


//////////////////////////////////////////////////////


   double TwoG_jan_8_2021::constrain_angle(double x) {
      double anglex;
      x = remainder(x, 360.0);
      if (x < 0) {
         x += 360.0;
      }
      anglex = (x - 180.0);
      return anglex;
   }

   double TwoG_jan_8_2021::angle_diff(double angle1, double angle2) {
      double diff = angle1 - angle2;
      while (diff < -180.0) {
         diff += 360.0;
      }
      while (diff > 180.0) {
         diff -= 360.0;
      }
      return diff;
   }

   double TwoG_jan_8_2021::delta_energy(int n_coin, float tagger_energy[30], double particle_ener) {
      double del_ener = 0.0;
      for (int i = 0; i < n_coin; i++) {
         del_ener = tagger_energy[i] - particle_ener;
      }
   return del_ener;
   }