#define emstudy_cxx
#include "emstudy.h"
#include "TH2.h"  // added by me

void emstudy::Loop()
{
//   In a Root session, you can do:
//      Root > .L emstudy.C
//      Root > emstudy t
//      Root > t.GetEntry(12); // Fill t data members with entry number 12
//      Root > t.Show();       // Show values of entry 12
//      Root > t.Show(16);     // Read and show values of entry 16
//      Root > t.Loop();       // Loop on all entries
//

//     This is the loop skeleton
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fTree->SetBranchStatus("*",0);  // disable all branches
//    fTree->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fTree->GetEntry(i);  // read all branches
//by  b_branchname->GetEntry(i); //read only this branch


  // define some variables

  Int_t em1 ; // added by me
  Int_t em2 ; // added by me
  Int_t k;    // added by me
  Real_t invmass; // added by my
   

   // define some histograms
   TH1F *h_nempart = new TH1F("h_nempart","n em particles",10,0.,10.);
   TH1F *h_emE = new TH1F("h_emE","emE",200,0.,200.);
   TH1F *h_empx = new TH1F("h_empx","empx",200,0.,200.);
   TH1F *h_empy = new TH1F("h_empy","empy",200,0.,200.);
   TH1F *h_empz = new TH1F("h_empz","empz",200,0.,200.);
   TH1F *h_eemass = new TH1F("h_eemass","electron mass",200,0.,200.);
   TH1F *h_zemfrac = new TH1F("h_zemfrac","emfrac for peak",100.,0.,1.);
   TH1F *h_semfrac = new TH1F("h_semfrac","emfrac for side",100.,0.,1.);


   if (fTree == 0) return;

   Int_t nentries = Int_t(fTree->GetEntries());


   Int_t nbytes = 0, nb = 0;


   // event loop;

   for (Int_t i=0; i<nentries;i++) {
      if (LoadTree(i) < 0) break;
      nb = fTree->GetEntry(i);   nbytes += nb;
      h_nempart->Fill(EMPART_nempart);


      // emobject loop rest is added by me

     for(em1 =0; em1 <EMPART_nempart; em1++){
       if(emobject(em1)){
	h_emE->Fill(EMPART_emE[em1]);
	h_empx->Fill(EMPART_empz[em1]);
	h_empy->Fill(EMPART_empy[em1]);
	h_empz->Fill(EMPART_empz[em1]);
	for(em2 = 0; em2 < em1; em2++){
	  if(emobject(em2)){
	    invmass = (EMPART_emE[em1]+EMPART_emE[em2])*(EMPART_emE[em1]+EMPART_emE[em2]);
	    invmass -= (EMPART_empx[em1]+EMPART_empx[em2])*(EMPART_empx[em1]+EMPART_empx[em2]);
	    invmass -= (EMPART_empy[em1]+EMPART_empy[em2])*(EMPART_empy[em1]+EMPART_empy[em2]);
	    invmass -= (EMPART_empz[em1]+EMPART_empz[em2])*(EMPART_empz[em1]+EMPART_empz[em2]);

	    invmass = sqrt(invmass);
	    h_eemass->Fill(invmass);
	    if(zpeak(invmass)){ 
	      h_zemfrac->Fill(EMPART_emfrac[em1]);
	      h_zemfrac->Fill(EMPART_emfrac[em2]);
	    }
	    elseif(sideband(invmass)){
	      h_semfrac->Fill(EMPART_emfrac[em1]);
	      h_semfrac->Fill(EMPART_emfrac[em2]);
	    }    
	  }
	}
       }
     }  


   }
   
}

/// added by me

bool emstudy::emobject(const Int_t index){
  return 1;
}

bool zpeak(const Real_t invmass){
  if(invmass>96.)return 0;
  if(invmass<86.) return 0;
  return 1;
}
bool lowband(const Real_t invmass){
  if(invmass>85.)return 0;
  if(invmass<80.)return 0;
  return 1;
}
bool highband(const Real_t invmass){
  if(invmass>97.)return 0;
  if(invmass<102.)return 0;
  return 1;
}
bool sideband(const Real_t invmass){
   return lowband(invmass)|| highband(invmass);
}