//|||||||||||||||||||||||||||||||||||||||||||||||||||||//
/* Create template class for Selector to run on a tree */
//|||||||||||||||||||||||||||||||||||||||||||||||||||||//

// open root file containing the Tree in ROOT
   TFile f("r8493-3.root");
// create TTree object from it
   TTree *t; f.GetObject("h1", t);
// this generates the files MySelector.h and MySelector.C
   t->MakeSelector("MySelector");

/*------------------------------------------------------------------------------------------------*/
/*                                             OR                                                 */
/* You can be in the directory where a file is (shown below) or just point to it and open in ROOT */
/*------------------------------------------------------------------------------------------------*/
mcintyre@gluey /nfs/direct/annex/mcintyre/RadPhi/test
$ root -l r8493-3.root 
root [0] 
Attaching file r8493-3.root as _file0...
(TFile *) 0x46cdc90
root [1] .ls
TFile**		r8493-3.root	HBOOK file: ntuple.hbook converted to ROOT
 TFile*		r8493-3.root	HBOOK file: ntuple.hbook converted to ROOT
  KEY: TTree	h1;9	Radphi
  KEY: TTree	h1;8	Radphi
// Since h1 is my TTree
root [2] h1->MakeSelector("MY_NEW_FILE");
/*------------------------------------------------------------------------------------------------*/


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


//||||||||||||||||||||||||||||||||||||||||||||||//
/*   YOUR C++ FILE CAN INCLUDE THE FOLLOWING:   */
/*       Placed after all your #includes        */
//||||||||||||||||||||||||||||||||||||||||||||||//

// Example include files
/* (depends on what's in your analysis code, as to what's needed) */
#include <iostream>
#include "YOUR_C++_FILE_NAME.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>

/* These are placed at the beginning of my C++ file after #include. */
// I define a function for squaring numbers.
#define sqr(x) ((x)*(x))

// Accidental Scaling Factor (ASF), value calculated previously.
// Used in determining a tag weighting factor.
#define ASF 0.935198

// Energy Rescaler (ER), value calculated previously.
// Used to shift invariant mass plot to match well
// known mesons. I determined its value using the
// pi0, eta, and etaprime mesons in the two gamma sample.
#define ER 0.97              

// I defined Pi, useful for converting to/from radians(degrees), etc.
#define PI 3.14159265359


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


   ////////////// One method for calculating invariant mass //////////////
   ///////////////////////////////////////////////////////////////////////
   /*
      Remember that nphot is LGD photons (0, ..., nphot - BGV_photons) plus
      the photons allowed in the BGV (n - BGD_photons, ..., nphot).
      Where nphot = number of photons in LGD + BGD for an event.
      I implement a Barrel Gamma Veto (BGV), so no BGD photons therefore 
      nphot = LDG photons.
   */

   //// Summing the components of the photons' 4-vector. ////
   double e_lab = 0.0;
   double px_lab = 0.0;
   double py_lab = 0.0;
   double pz_lab = 0.0;
   // Below, you can also specify number of photons instead of using nphot.
   for (int m = 0; m < nphot; m++) {
      e_lab   +=  pvect[m][0];  // Energy
      px_lab  +=  pvect[m][1];  // p_x
      py_lab  +=  pvect[m][2];  // p_y
      pz_lab  +=  pvect[m][3];  // p_z
   }
 
   // Calculating the invariant mass
   double imass = 0.0;
   imass = sqrt(sqr(e_lab) - sqr(px_lab) - sqr(py_lab) - sqr(pz_lab));


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


   //////////// Another method for calculating invariant mass ////////////
   /////////// I use this method along with an energy rescaler ///////////
   ///////////////////////////////////////////////////////////////////////

   // At the beginning of my code I introduce an energy rescaler, as follows.
   // Rescaling the LGD. Efrwd is from the .root files.
   double EFRWD = ER * Efrwd;

   // Rescaling the photons. pvect[][] is from the .root files.
   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 up to six photons
   TLorentzVector photon_0;
   TLorentzVector photon_1;
   TLorentzVector photon_2;
   TLorentzVector photon_3;
   TLorentzVector photon_4;
   TLorentzVector photon_5;

   // Remember: pvect[][] (lower case) comes from the .root data file,
   // but since I used an energy rescaler (ER) I made new
   // variables PVECT[][] (upper case). See rescaling above.
   // If you don't need to rescale then just do the same as below using pvect[][].
   /* 
      pvect[photon][4-vector]
      where: 4-vector --> [0] = energy, [1] = p_x, 
                          [2] = p_y,    [3] = p_z
   */  
   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]);
   photon_2.SetPxPyPzE(PVECT[2][1],PVECT[2][2],PVECT[2][3],PVECT[2][0]);
   photon_3.SetPxPyPzE(PVECT[3][1],PVECT[3][2],PVECT[3][3],PVECT[3][0]);
   photon_4.SetPxPyPzE(PVECT[4][1],PVECT[4][2],PVECT[4][3],PVECT[4][0]);
   photon_5.SetPxPyPzE(PVECT[5][1],PVECT[5][2],PVECT[5][3],PVECT[5][0]);

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


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


   ////////////// Calculating Mandelstam variable |t| //////////////
   //////////////    The momentum transfer squared    //////////////
   //////////////  t = (p_1 - p_3)^2 = (p_4 - p_2)^2  //////////////
   //////////////  This is often written a -t or |t|  //////////////
   /////////////////////////////////////////////////////////////////
   double tabsolute = 0.0; 
   for (int taghit = 0; taghit < ncoin; taghit++) {
      // There is a separate |t| for each coincidence hit (ncoin)
      tabsolute = -(sqr(imass) - 2 * (coenergy[taghit] * (e_lab - pz_lab)));

      /* This is the spot within the ncoin loop where you can do a 
         histogram Fill. We use tag weighting (accidental subtraction)
         in the Fill statement, as follows. */

      // Tag weight example
      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[taghit].
      dt = trec0 - cotime[taghit];
      // These dt ranges are based off of looking at a timing histogram.
      if (dt >= -1.0 && dt <= 5.0) {
         // Range used for "true" hits
         tag_wgt = 1.0;
      }

      else if (dt > 5.0 && dt <= 11.0) {
         // Range used for "accidental" hits
         tag_wgt = -(1.0 / ASF);
      }

      else {
         // If coincidence hit is outside timing ranges, ignore it
         tag_wgt = 0.0;
      }

		// Example Fill statements: 
      // The code below will do a fill for each coincidence hit and will have a weighting
      // factor based on where the hit (ncoin) falls in dt (timing difference).
      // Remember: taghit cycles through all coincidence hits. taghit = 0, 1, 2, ..., ncoin.
      // If (dt >= -1.0 && dt <= 5.0): We consider this a "true" hit and give it 
      //                               a fill weight of +1.
      // If (dt > 5.0 && dt <= 11.0):  We consider this to be a range representative
      //                               of what the accidental background looks like 
      //                               in the 6 ns "true" hit range and therefore apply
      //                               a negative weight so that we can subtract out the
      //                               accidental hits.
      // If dt is outside the aforementioned ranges then we apply a weight of zero (0).

		tabs->Fill(tabsolute, tag_wgt);
      Inv_Mass->Fill(imass2, tag_wgt);

      /* 
         NOTE: The weighting factor for Fill statements can consist of multiple
               components. For example, we can include a mass sideband subtraction
               weight factor and/or a delta phi SB subtraction factor.
               tabs->Fill(tabsolute, tag_wgt * mass_wgt * dphi_wgt);
      */
   }



      // lab frame four momenta of the incident beam,incident proton,
      // omega meson, recoiling proton.

      // missing mass and the decay plane azimuthal angle.

      // Boosting to center of mass frame. cmboost is the boost vector  

      // COM four vectors. 


      /* rotate about the y axis by the polar angle */


      /* Now boost to omega rest frame */ 


      // we are now in helicity frame, quantization axis is defined
      // as the direction opposite to the recoiling nucleon

      //  alpha angle, angle the beam makes with the z axis
      //  (theta in helicity frame),Since GJ frame and helicity frame differ
      //  by just a rotation.A rotation about the y-axis by the polar angle 
      //  gets you to GJ frame , in this frame the quantization axis is 
      //  defined as the direction of the incident beam in the omega meson rest frame 


      /* omega meson decay products angular distribution */


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


   /////////// Getting Acceptance & Resolution with Radphi MC  ///////////
   ///////////////////////////////////////////////////////////////////////

   // Declare variables
   double t_gen = 0;
   double t_acc = 0;
   double t_res = 0;
   double mass_res = 0;
   double gen_mass = 0;
   double p[4] = {0.0, 0.0, 0.0, 0.0};
   // Below, you can also specify number of photons instead of using nphot.
   // This is the same as I did above with pvect, except momf is for MC data.
   for (int i = 0; i < nphot; i++) {
      p[0] += momf[i][0];
      p[1] += momf[i][1];
      p[2] += momf[i][2];
      p[3] += momf[i][3];
   }

   // Calculate invariant mass for MC data.
   gen_mass = sqrt(sqr(p[0]) - sqr(p[1]) - sqr(p[2]) - sqr(p[3]));

   // Calculate |t| generated from MC data. momi[0][0] = photon beam energy
   t_gen = -(sqr(gen_mass) - 2 * momi[0][0] * (p[0] - p[3]));

   // Calculations and fills for getting acceptace (used in calculating cross section) 
   // and resolution used to help determine bin widths.
   t_res =  tabsolute - t_gen;
   t_res1d->Fill(t_res);
   // Acceptance
   t_acc = tabsolute / t_gen;
   t_acc1d->Fill(t_acc);
   // Resolution
   tabs_res->Fill(t_gen, tabsolute);


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


   //////////// Calculating delta phi and transverse momentum ////////////
   //////////////////////// using TLorentzVector /////////////////////////

   TLorentzVector prot_lab;     // Proton
   TLorentzVector recprot_lab;  // Recoil proton
   TLorentzVector meson_lab;    // Meson of interest
   TLorentzVector B_lab;        // Photon beam

   prot_lab.SetPxPyPzE(0.0, 0.0, 0.0, 0.938);            // Proton at rest
   meson_lab.SetPxPyPzE(px_lab, py_lab, pz_lab, e_lab);  // Variables defined above
                                                         // using photons' summed pvect[][]
   /*
      delta phi (del_phi) =  the azimuthal angle difference between the recoil 
      proton plane and the meson of interest. Delta phi is used to suppress and 
      proton resonances. Gives the angle -180 degree different. Since recoil proton
      and meson phi should be 180 degree different.
   */

   double pTy = 0.0;
   double phi_rec = 0.0;
   double rec_phi = 0.0;
   double meson_phi = 0.0;
   double del_phi = 0.0;
   double pt = 0.0;
   phi_rec = phirec[0] * (180.0 / PI);  // Converts radians to degrees
   // constrain_angle() is a function defined at the end of my code.
   // Gives the angle -180 degree different. 
   /* Since recoil proton and meson phi should be 180 degree different. */
   rec_phi = constrain_angle(phi_rec);
   pTy = -(sin(phirec[0]) * px_lab) + cos(phirec[0]) * py_lab;
   pt = meson_lab.Perp();
   meson_phi = meson_lab.Phi() * (180.0 / PI);  // Gets phi angle of the meson then converts radians to degrees

   // angle_diff() is a function defined at the end of my code
   // Takes the difference and +/- 360 to keep the value between +/- 180 deg.
   del_phi = angle_diff(meson_phi, rec_phi);

   // ... more analysis code. Then, (see below)

   /*
      The following code placed after Terminate(), at the end of your file. 
      This function must be declared in your .h file and can be called within
      your analysis code, like I did above. Useful if you need to call multiple
      times using different inputs at different places in your code. 
   */

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

   double YOUR_C++_FILE_NAME::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;
   }


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


//||||||||||||||||||||||||||||||||||||||||||||||//
/* YOUR HEADER FILE WILL INCLUDE THE FOLLOWING: */
//||||||||||||||||||||||||||||||||||||||||||||||//

#ifndef YOUR_C++_FILE_NAME_h
#define YOUR_C++_FILE_NAME_h

#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TSelector.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TString.h>
#include <TProfile.h>
#include <TTree.h>

class YOUR_C++_FILE_NAME : public TSelector {
public :
   TTree          *fChain;   //!pointer to the analyzed TTree or TChain
   Int_t           fCurrent; //!current Tree number in a TChain

   double constrain_angle(double x);
   double angle_diff(double angle1, double angle2);

   // Declare histograms
   TH1D *tPi;
   TH1D *dphi_Pi;
   TH1D *Pi_dphi[28];  // This one is 28 TH1D histos named:
                       // Pi_dphi0, Pi_dphi1, ..., Pi_dphi27

   // Fixed size dimensions of array or collections stored in the TTree, if any.

   // Declaration of leaf types
      /*
         The rest comes from a class that you can automatically generated from a TTree
         in a .root file you specify. When doing analysis you'll generate it once,
         modify to work well with PROOF, then reuse with subsequent analysis
         (obviously modified to your analysis code,
          i.e. histo names, function names, etc.)
      */

//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^//
/* YOUR HEADER FILE WILL INCLUDE THE ABOVE CODE */
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^//


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


   double recprot_0 = 0.0;
   double recprot_1 = 0.0;
   double recprot_2 = 0.0;
   double recprot_3 = 0.0;
   double t_abs = 0.0;
   TVector3 cmboost;
   TLorentzVector B_lab;        // Photon beam
   TLorentzVector prot_lab;     // Incident proton
   TLorentzVector meson_lab;    // Meson of interest
   TLorentzVector recprot_lab;  // Recoil proton

   // In the analysis you can have kinematic cuts. For example with the 3 gamma final state 
   // sample: omega -> pi0 * gamma, pi0 -> 2 gammas. Most of the time the pi0 is the
   // lowest mass 2 gamma combo. So if we look at the three different 2 gamma combos
   // that can be uniquely made, then select the lowest mass combo that falls within
   // the pi0 mass range, the other photon (third) will be your "gamma_lab" photon.
   // To get the mass range to cut on: Plot a TH2D of 3 gamma invariant mass vs all 
   // 2 gamma combos. This TH2D will have three entries per event. Then look at the 
   // histogram to find a suitable mass range to cut on for the pi0 invariant mass.      
   TLorentzVector Pi_lab;       // 2 photon combo in pi0 mass range
   TLorentzVector gamma_lab;    // Photon in final state not associated to a meson
   //-------------------------------------------------------------------------------//
   //-- Analysis code to fill the two TLorentzVectors above is not included here. --//
   //-------------------------------------------------------------------------------//

   /* Summing the components of the photons' 4-vectors. */
   double e_lab = 0.0;
   double px_lab = 0.0;
   double py_lab = 0.0;
   double pz_lab = 0.0;
   // Below, you can also specify number of photons instead of using nphot.
   for (int m = 0; m < nphot; m++) {
      e_lab   +=  pvect[m][0];  // Energy
      px_lab  +=  pvect[m][1];  // p_x
      py_lab  +=  pvect[m][2];  // p_y
      pz_lab  +=  pvect[m][3];  // p_z
   }

   double min_delta_kin_energy = 0.05;
   for (int taghit = 0; taghit < ncoin; taghit++) {
      // Get the recoil momentum using kinetic energy conservation
      // and the polar and azimuthal angles
      double rec_prot_ke = 0.0;
      double rec_prot_pmag = 0.0;
      double rec_prot_pt = 0.0;
      double prot_mass = 0.938;
      double proton_binding_energy = 0.03;

      // Recoil proton KE = Beam Energy - Energy of decay photons - Proton binding energy
      rec_prot_ke = coenergy[taghit] - e_lab - proton_binding_energy;

      // Gives the recoiling proton a minimum KE of 50 MeV
      if (rec_prot_ke > min_delta_kin_energy) {
         recprot_0 =  rec_prot_ke + prot_mass;
      }
      else {
         recprot_0 = min_delta_kin_energy + prot_mass;
      }

      // Below comes from -> E^2 = (m * c^2)^2 + (p * c)^2 
      //              where: c = 1
      //                     p = rec_prot_pmag,
      //           remember: p is a 3-vector (lab frame)
      rec_prot_pmag = sqrt(sqr(recprot_0) - sqr(prot_mass));

      recprot_3 = rec_prot_pmag * cos(therec[0]);  // z-direction

      // pt is transverse momentum (x-y plane)
      // p^2 = px^2 + py^2 + pz^2
      // (rec_prot_pt)^2 = px^2 + py^2 = p^2 - pz^2
      rec_prot_pt = sqrt(sqr(rec_prot_pmag) - sqr(recprot_3));

      recprot_1 = rec_prot_pt * cos(phirec[0]);    // x-direction
      recprot_2 = rec_prot_pt * sin(phirec[0]);    // y-direction

      // The difference in energy between the incident beam and the meson
      // produced = energy transfered to the proton.
      // t_abs is the mandelstam variable |t|, the four-momentum transfer squared.
      t_abs = -(sqr(imass) - 2 * coenergy[taghit] * (e_lab- pz_lab));

      // Two example of an entry for a new .root inside your code.
      // Useful for later analysis.
      pwa_ntuple.tabs = -(sqr(imass) - 2 * (coenergy[taghit] * (e_lab- pz_lab)));
      pwa_ntuple.Ebeam = coenergy[taghit];

      // Lab frame four momenta of the incident beam, incident proton, meson, & recoiling proton.
      /* Incident photon beam is in + z-direction */
      B_lab.SetPxPyPzE(0.0, 0.0, coenergy[taghit], coenergy[taghit]);
      prot_lab.SetPxPyPzE(0.0, 0.0, 0.0, 0.938);
      meson_lab.SetPxPyPzE(px_lab, py_lab, pz_lab, e_lab);
      recprot_lab.SetPxPyPzE(recprot_1, recprot_2, recprot_3, recprot_0);

      // Decay plane polar & azimuthal angles
      double thetalab = meson_lab.Theta();
      double philab = meson_lab.Phi();

      // Boosting to center of mass frame. cmboost is the boost vector 
      // since the initial state proton is not measured(nuclear target),define boost vector 
      // using meson and recoiling proton 
      TLorentzVector cms = recprot_lab + meson_lab;
      cmboost = -cms.BoostVector(); 
         //double bcm = B_lab.E()/(B_lab.E() + prot_lab.M()); 
         //cmboost.SetXYZ(0.0, 0.0, bcm);

      // Center of mass 4-vectors 
      TLorentzVector B_cm(B_lab);
      TLorentzVector prot_cm(prot_lab);
      TLorentzVector meson_cm(meson_lab);
      TLorentzVector recprot_cm(recprot_lab);
      TLorentzVector Pi_cm(Pi_lab);
      TLorentzVector gamma_cm(gamma_lab);

      // Boosting to COM
      B_cm.Boost(cmboost);
      prot_cm.Boost(cmboost);
      meson_cm.Boost(cmboost);
      recprot_cm.Boost(cmboost);
      Pi_cm.Boost(cmboost);
      gamma_cm.Boost(cmboost);

      // COM plane polar & azimuthal angles
      double thetacm = meson_cm.Theta();
      double phicm = meson_cm.Phi();

      /* Rotate about the z-axis by phicm */
      B_cm.RotateZ(-phicm);
      prot_cm.RotateZ(-phicm);
      meson_cm.RotateZ(-phicm);
      recprot_cm.RotateZ(-phicm);
      Pi_cm.RotateZ(-phicm);
      gamma_cm.RotateZ(-phicm);

      /* Rotate about the y-axis by the polar angle */
      B_cm.RotateY(-thetacm);
      prot_cm.RotateY(-thetacm);
      meson_cm.RotateY(-thetacm);
      recprot_cm.RotateY(-thetacm);
      Pi_cm.RotateY(-thetacm);
      gamma_cm.RotateY(-thetacm);
      
      /* Now boost to meson rest frame */
      // Boot to helicity frame 
      double meson_rest = meson_cm.Pz() / meson_cm.E();

      TVector3 rfboost;
      rfboost.SetXYZ(0.0, 0.0, meson_rest);

      TLorentzVector B_hf(B_cm);
      TLorentzVector prot_hf(prot_cm);
      TLorentzVector meson_hf(meson_cm);
      TLorentzVector recprot_hf(recprot_cm);
      TLorentzVector Pi_hf(Pi_cm);
      TLorentzVector gamma_hf(gamma_cm);

      B_hf.Boost(-rfboost);
      prot_hf.Boost(-rfboost);
      meson_hf.Boost(-rfboost);
      recprot_hf.Boost(-rfboost);
      Pi_hf.Boost(-rfboost);
      gamma_hf.Boost(-rfboost);

      // We are now in helicity frame, quantization axis is defined
      // as the direction opposite to the recoiling nucleon

      double hfcosthe = Pi_hf.Dot(-recprot_hf);
      double helcosthe = cos(Pi_hf.Theta());
      double helphi    = (Pi_hf.Phi());
      double helthe    = (Pi_hf.Theta());
      // double helcosthe = cos(gamma_hf.Theta());
      // double helphi    = (gamma_hf.Phi());
      // double helthe    = (gamma_hf.Theta());

      // Alpha angle is the angle the beam makes with the z axis
      // (theta in helicity frame),Since GJ frame and helicity frame differ
      // by just a rotation.A rotation about the y-axis by the polar angle 
      // gets you to GJ frame , in this frame the quantization axis is 
      // defined as the direction of the incident beam in the omega meson rest frame 

      double alpha = B_hf.Theta();

      TLorentzVector B_gj(B_hf);
      TLorentzVector prot_gj(prot_hf);
      TLorentzVector meson_gj(meson_hf);
      TLorentzVector recprot_gj(recprot_hf);
      TLorentzVector Pi_gj(Pi_hf);
      TLorentzVector gamma_gj(gamma_hf);

      B_gj.RotateY(alpha);
      prot_gj.RotateY(alpha);
      meson_gj.RotateY(alpha);
      recprot_gj.RotateY(alpha);
      Pi_gj.RotateY(alpha);
      gamma_gj.RotateY(alpha);

      /* Meson decay products angular distribution */
      double gjcosthe = cos(Pi_gj.Theta());
      double gjphi    = Pi_gj.Phi();
      double gjthe    = Pi_gj.Theta();
      // double gjcosthe = cos(gamma_gj.Theta());
      // double gjphi    = gamma_gj.Phi();
      // double gjthe    = gamma_gj.Theta();
   }


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


   // Find two lightest photons that reconstruct a pion, the third photon 
   // becomes the bachelor photon in the omega -> pi0 * gamma decay 

   TLorentzVector photon_1;
   TLorentzVector photon_2;
   TLorentzVector photon_3;
   TLorentzVector gamma_lab;
   TLorentzVector Pi;
   TLorentzVector Pi1;
   TLorentzVector Pi2;
   TLorentzVector Pi_lab;
   int i, j, k;
   int i_lgt, j_lgt;
   for (i = 0; i < nphot; i++) {
      photon_1.SetPxPyPzE(pvect[i][1], pvect[i][2], pvect[i][3], pvect[i][0]);
      for (j =i + 1; j < nphot; j++) {
         photon_2.SetPxPyPzE(pvect[j][1], pvect[j][2], pvect[j][3], pvect[j][0]);
         for (k = j + 1; k < nphot; k++) {
            photon_3.SetPxPyPzE(pvect[k][1], pvect[k][2], pvect[k][3], pvect[k][0]);
            Pi  = photon_1 + photon_2;
            Pi1 = photon_1 + photon_3;
            Pi2 = photon_2 + photon_3;
         }
      }
   }        

   if (Pi.M() < Pi1.M() && Pi.M() < Pi2.M()) {
	   Pi_lab = Pi;
      gamma_lab = photon_3;
   }

   else if (Pi1.M() < Pi.M() && Pi1.M() < Pi2.M()) {
	   Pi_lab = Pi1;            
      gamma_lab = photon_2;
   }

   else if (Pi2.M() < Pi.M() && Pi2.M() < Pi1.M()) {
	   Pi_lab = Pi2;            
      gamma_lab = photon_1;
   }