// Z->ee and W->enu actual data, does not look at plug electrons
// tight and loose electrons (CC and TCC)
//for run number 141544 through 163527, from /cdf/data10a/Datasets/non_express/electron/UCntuples/4.8.4a_remake4.11.1_tight/
// Not a fully debugged code.  Please replace this with a fully debugged code.
// August 22, 2003.
// use with GoodRunNoSi.C and load.C


#include <iostream>
#include <math.h>
#include "TChain.h"
#include "TTree.h"
#include "TFile.h"
#include "TLeaf.h"
#include "TBranch.h"
#include "TObjString.h"
#include "TClonesArray.h"
#include "UCNtuple/myuber_t.hh"

#include "TH1F.h"


void declarestructs(TTree* tree, myuber_t* uber);

// Function definitions. Implementations in GoodRunFcn.C, GoodRunFcnNoSi.C.
bool good_run_nosi(int run);
//Baseline track quality cut
bool GoodTrackQ(mytrack_t& mytrack, int i);
//Baseline electron cut
bool StandardElectronQ(myelectron_t& myelectron, int i, mytrack_t& mytrack);
bool LooseElectronQ(myelectron_t& myelectron, int i, mytrack_t& mytrack);
//Inv mass calculation.
Float_t invmass(Float_t et1, Float_t theta1, Float_t phi1, Float_t et2, Float_t theta2, Float_t phi2);



TFile *Zfile = new TFile("Zeewintsum3.root", "RECREATE");

void zsearch(TTree* oldtree) {
  //begin Z search


myuber_t uber;
  TTree* markchange;
  if (! oldtree) return;
  cout << "Expect to process " << oldtree->GetEntries() << " entries" << endl;
  oldtree->SetBranchStatus("*", 1);
  oldtree->GetEntry(0);


// forward backward electron asymmetry
 TH1F *Afb = new TH1F("Afb", "Forward Backward Electron Asymmetry Versus Z Mass", 7, 60, 200);
 TH1F *Afb_w = new TH1F("Afb_w", "Forward Backward Electron Asymmetry Versus Z Mass (Winter)", 7, 60, 200);
 TH1F *Afb_s = new TH1F("Afb_s", "Forward Backward Electron Asymmetry Versus Z Mass (Summer)", 7, 60, 200);
 TH1F *zee_lepton = new TH1F("zee_lepton", "Z to ee Lepton Mass 4.11.1", 100, 0, 200);
 TH1F *zee = new TH1F("zee", "Z to ee Invariant Z Mass 4.11.1 between 66 and 116", 100, 0, 200);
 TH1F *zee_winter = new TH1F("zee_winter", "Z to ee Invariant Z Mass (Winter 02-03) 4.11.1", 100, 0, 200);
 TH1F *zee_summer = new TH1F("zee_summer", "Z to ee Invariant Z Mass (Summer 03) 4.11.1", 100, 0, 200);
 TH1F *ZEta = new TH1F("ZEta", "Eta of Z Candidate Electrons", 100, -1.5, 1.5);
 TH1F *runnumber = new TH1F("runnumber", "The Run Numbers of Passing Events", 100,140000,164000);
 TH1F *runnumber_w = new TH1F("runnumber_w", "The Run Numbers of Passing Events (winter)", 100,140000,164000);
 TH1F *runnumber_s = new TH1F("runnumber_s", "The Run Numbers of Passing Events (summer)", 100,140000,164000);
 TH1F *PEta = new TH1F("PEta", "Number of Positive Electons Versus Eta", 100, -1.5, 1.5);
 TH1F *NEta = new TH1F("NEta", "Number of Negative Electrons Versus Eta", 100, -1.5, 1.5);
 TH1F *Phi = new TH1F("Phi", "Phi of Z Candidate Electrons", 100, -3.14159, 3.14159*2);

 TH1F *zeeTCC = new TH1F("zeeTCC", "Z to ee Invariant Z Mass 4.11.1 TCC", 100, 0, 200);
 TH1F *ZEtaTCC = new TH1F("ZEtaTCC", "Eta of Z Candidate Electrons TCC", 100, -1.5, 1.5);
 TH1F *PhiTCC = new TH1F("PhiTCC", "Phi of Z Candidate Electrons TCC", 100, -3.14159, 3.14159*2);

 TH1F *zeeCC = new TH1F("zeeCC", "Z to ee Invariant Z Mass 4.11.1 CC", 100, 0, 200);
 TH1F *ZEtaCC = new TH1F("ZEtaCC", "Eta of Z Candidate Electrons CC", 100, -1.5, 1.5);
 TH1F *PhiCC = new TH1F("PhiCC", "Phi of Z Candidate Electrons CC", 100, -3.14159, 3.14159*2);
 TH1F *WMasst = new TH1F("WMasst", "Transverse Mass of W->e nu", 200, 0, 200);
 TH1F *WMasst_winter = new TH1F("WMasst_winter", "Transverse Mass of W->e nu (Winter 02-03)", 200, 0, 200);
 TH1F *WMasst_summer = new TH1F("WMasst_summer", "Transverse Mass of W->e nu (Summer 03)", 200, 0, 200);

 TH1F *zeeTCC_winter = new TH1F("zeeTCC_winter", "Z to ee Invariant Z Mass 4.11.1 TCC (winter)", 100, 0, 200);
 TH1F *zeeCC_winter = new TH1F("zeeCC_winter", "Z to ee Invariant Z Mass 4.11.1 CC (winter)", 100, 0, 200);


// finding the first and last run numbers
   int r = 200000;
   int r2 = 0; 

// Counting the number of runs and events
   int runno_new = 0;
   int numruns = 0;
   int eventno_new = 0;
   int numevents = 0;

// forward backward electron asymmetry
   Float_t A[100];
   Float_t F[100];
   Float_t B[100];
   Float_t Aw[100];
   Float_t Fw[100];
   Float_t Bw[100];
   Float_t As[100];
   Float_t Fs[100];
   Float_t Bs[100];

for (int o=0; o<100; o++)
   {
     A[o] = 0.0;     
     F[o] = 0.0;
     B[o] = 0.0;
     Aw[o] = 0.0;     
     Fw[o] = 0.0;
     Bw[o] = 0.0;
     As[o] = 0.0;     
     Fs[o] = 0.0;
     Bs[o] = 0.0;
   }


  myuber_t uber;
  Int_t treeno = -1;

  if (! oldtree) return;
//  cout << "Expect to process " << oldtree->GetEntries() << " entries" << endl;
  oldtree->SetBranchStatus("*", 1);
  oldtree->GetEntry(0);
  
#ifdef __CINT__
  gInterpreter->LoadMacro("declarestructs.C");
  //  gInterpreter->LoadMacro("/cdf/data1e/ucntuple/goodrun.C");
#endif

  declarestructs(oldtree, &uber);
  
  int npass = 0; int npass_Z = 0; int nwenu = 0;

  zee.Reset(); Afb.Reset(); ZEta.Reset();

  for (int i = 0; i < oldtree->GetEntries(); i++) {
    // only way I can get this to be reliable
    oldtree->GetEntry(i);
    if(oldtree->IsA() == TChain::Class()){
      if (treeno != ((TChain*) oldtree)->GetTreeNumber()) {
	declarestructs(oldtree->GetTree(), &uber);
	oldtree->GetEntry(i);
	treeno = ((TChain*) oldtree)->GetTreeNumber();
      }
    }

// finding the first and last run numbers
   int runno = uber.header.run;

   if (runno < r)
   {
      r = runno;
   }
   
   if (runno > r2)
   {
      r2 = runno;
   }


   if (!(good_run_nosi(runno))) continue;

// Counting the number of runs and events
   int eventno = uber.header.event;
   if (runno_new != runno)
   {
      runno_new = runno;
      numruns++;
   }   

   if (eventno_new != eventno)
   {
      eventno_new = eventno;
      numevents++;
   }



    if (i % 1000 == 0) {
      cout << "entry: " << i << endl;
    }

    int nel = 0;
    int elind[MAXELE], elID[MAXELE];
    int ntcem=0,ncem=0;
    int tcemind[MAXELE],cemind[MAXELE];

    if (uber.electron.nele > 0) {
      for (int j = 0; j < uber.electron.nele; j++) {
	if(StandardElectronQ(uber.electron, j, uber.track)) {  //tight CEM
	  tcemind[nel] = j;
	  ntcem++;
	  elind[nel] = j;
	elID[nel] = 1;  
	nel++;
	   
	} else if(LooseElectronQ(uber.electron, j, uber.track)) {  //loose CEM
	  cemind[nel] = j;
	  ncem++;
	  elind[nel] = j;
	  
	  elID[nel] = 2;
	 nel++;
	} 

       }
     }
 
  
   if (nel == 2 && ntcem>0) {  //two electrons but at least one tight cem
      npass++;
      // cout << "Event passes: " << nel << " di-electrons" << endl; 


      // We use the Etcorr for inv. mass calculation.

      float e1_theta = 2.*atan2(exp(-(uber.electron.eveta[elind[0]]) ),1);
      float e1_phi = uber.electron.phi[elind[0]];
      float e1_et = uber.electron.etcorr[elind[0]];

      float e2_theta = 2.*atan2(exp(-(uber.electron.eveta[elind[1]]) ),1);
      float e2_phi = uber.electron.phi[elind[1]];
      float e2_et = uber.electron.etcorr[elind[1]];

      float massll = invmass(e1_et,e1_theta,e1_phi,e2_et,e2_theta,e2_phi);

      if ((massll > 66) && (massll < 116)) {
	npass_Z++;
        zee.Fill(massll);
	Phi.Fill(uber.electron.phi[elind[0]]);
	Phi.Fill(uber.electron.phi[elind[1]]);
	ZEta.Fill(uber.electron.eveta[elind[0]]);
        ZEta.Fill(uber.electron.eveta[elind[1]]);
        runnumber->Fill(uber.header.run);
	  if (uber.header.run <= 156487){
	   zee_winter.Fill(massll);
	   runnumber_w->Fill(uber.header.run);
	  }
	  if (uber.header.run > 156487){
	   zee_summer.Fill(massll);
	   runnumber_s->Fill(uber.header.run);
	  }
       }  

      zee_lepton.Fill(massll);

  if ((massll > 66) && (massll < 116)) {
      if(elID[0]==1 && elID[1]==1) {
	zeeTCC->Fill(massll);
	ZEtaTCC->Fill(uber.electron.eveta[elind[0]]);
	ZEtaTCC->Fill(uber.electron.eveta[elind[1]]);
        PhiTCC->Fill(uber.electron.phi[elind[0]]);
	PhiTCC->Fill(uber.electron.phi[elind[1]]);
	 if (uber.header.run <= 156487)
	   zeeTCC_winter.Fill(massll);


      } // TCEM+ TCEM
      if((elID[0]==1 && elID[1]==2) || (elID[0]==2 && elID[1]==1))  {
	zeeCC->Fill(massll);
	ZEtaCC->Fill(uber.electron.eveta[elind[0]]);
	ZEtaCC->Fill(uber.electron.eveta[elind[1]]);
	PhiCC->Fill(uber.electron.phi[elind[0]]);
	PhiCC->Fill(uber.electron.phi[elind[1]]);
	 if (uber.header.run <= 156487)
	   zeeCC_winter.Fill(massll);      

      } // TCEM+ CEM
   }


// To plot the number of positive and negative electrons versus Eta

	 if(uber.electron.charge[elind[0]]==1)
	     PEta->Fill(uber.electron.eveta[elind[0]]);
	
	 if(uber.electron.charge[elind[1]]==1)
	     PEta->Fill(uber.electron.eveta[elind[1]]);
	 
	 if(uber.electron.charge[elind[0]]==-1)
	     NEta->Fill(uber.electron.eveta[elind[0]]);
	 
	 if(uber.electron.charge[elind[1]]==-1)
	     NEta->Fill(uber.electron.eveta[elind[1]]);



// To plot asymmetry versus Z mass
     for(int o=1; o<=7; o++)
	{
	if((massll>=Afb->GetBinLowEdge(o))&&(massll<Afb->GetBinLowEdge(o+1)))
	{
	   if((uber.electron.charge[elind[0]]==-1)&&(uber.electron.eveta[elind[0]]>uber.electron.eveta[elind[1]]))
	   F[o-1]++;
	   if((uber.electron.charge[elind[1]]==-1)&&(uber.electron.eveta[elind[0]]>uber.electron.eveta[elind[1]]))
	   B[o-1]++;
	   if((uber.electron.charge[elind[0]]==-1)&&(uber.electron.eveta[elind[0]]<uber.electron.eveta[elind[1]]))
	   B[o-1]++;
	   if((uber.electron.charge[elind[1]]==-1)&&(uber.electron.eveta[elind[0]]<uber.electron.eveta[elind[1]]))
	   F[o-1]++;
	}
	if ((F[o-1]+B[o-1]) != 0)
	  A[o-1]=((F[o-1]-B[o-1])/(F[o-1]+B[o-1]));
	Afb->SetBinContent(o,A[o-1]);
        }

    

   if (uber.header.run <= 156487){
     for(int o=1; o<=7; o++)
       {
	 if((massll>=Afb_w->GetBinLowEdge(o))&&(massll<Afb_w->GetBinLowEdge(o+1)))
	   {
	     if((uber.electron.charge[elind[0]]==-1)&&(uber.electron.eveta[elind[0]]>uber.electron.eveta[elind[1]]))
	       Fw[o-1]++;
	     if((uber.electron.charge[elind[1]]==-1)&&(uber.electron.eveta[elind[0]]>uber.electron.eveta[elind[1]]))
	       Bw[o-1]++;
	     if((uber.electron.charge[elind[0]]==-1)&&(uber.electron.eveta[elind[0]]<uber.electron.eveta[elind[1]]))
	       Bw[o-1]++;
	     if((uber.electron.charge[elind[1]]==-1)&&(uber.electron.eveta[elind[0]]<uber.electron.eveta[elind[1]]))
	       Fw[o-1]++;
	   }
	 if ((Fw[o-1]+Bw[o-1]) != 0)
	   Aw[o-1]=((Fw[o-1]-Bw[o-1])/(Fw[o-1]+Bw[o-1]));
	 Afb_w->SetBinContent(o,Aw[o-1]);
       }
     
   }
   

   if (uber.header.run > 156487){
     for(int o=1; o<=7; o++)
       {
	 if((massll>=Afb_s->GetBinLowEdge(o))&&(massll<Afb_s->GetBinLowEdge(o+1)))
	   {
	     if((uber.electron.charge[elind[0]]==-1)&&(uber.electron.eveta[elind[0]]>uber.electron.eveta[elind[1]]))
	       Fs[o-1]++;
	     if((uber.electron.charge[elind[1]]==-1)&&(uber.electron.eveta[elind[0]]>uber.electron.eveta[elind[1]]))
	       Bs[o-1]++;
	     if((uber.electron.charge[elind[0]]==-1)&&(uber.electron.eveta[elind[0]]<uber.electron.eveta[elind[1]]))
	       Bs[o-1]++;
	     if((uber.electron.charge[elind[1]]==-1)&&(uber.electron.eveta[elind[0]]<uber.electron.eveta[elind[1]]))
	       Fs[o-1]++;
	   }
	 if ((Fs[o-1]+Bs[o-1]) != 0)
	   As[o-1]=((Fs[o-1]-Bs[o-1])/(Fs[o-1]+Bs[o-1]));
	 Afb_s->SetBinContent(o,As[o-1]);
       }
     
   }
   }





// Selecting the W's
   if (ntcem == 1 && uber.met.metraw > 25) {
     nwenu++;
     double transmass = pow(uber.electron.etcorr[tcemind[0]] + uber.met.metraw,2)-
				pow(uber.electron.px[tcemind[0]] + uber.met.metraw*cos(uber.met.metrawphi),2)-
				pow(uber.electron.py[tcemind[0]] + uber.met.metraw*sin(uber.met.metrawphi),2);
     transmass = sqrt(transmass);
     WMasst->Fill(transmass);
        if (uber.header.run <= 156487)
	  WMasst_winter.Fill(transmass);
	if (uber.header.run > 156487)
	  WMasst_summer.Fill(transmass);
    }


  }
  cout << "Number of passing Z events: " << npass << endl;
  cout << "Number of passing Z events between 66 and 116 GeV/c^2: " << npass_Z << endl;
  cout << "Number of passing W events: " << nwenu << endl;


// Counting the number of runs and events
  cout << "   The total number of runs: " << numruns << endl;
  cout << "   The total number of events: " << numevents << endl;

// finding the first and last run numbers
  cout << "   The first run number:  "<<r<<endl;
  cout << "   The last run number:  "<<r2<<endl;

// Forward backward muon asymmetry
  cout << "   # of forward events for Z's with a mass >60 and <80: " << F[0] << " +/- " << sqrt(F[0]) << endl;
  cout << "   # of backward events for Z's with a mass >60 and <80: " << B[0] << " +/- " << sqrt(B[0]) << endl;
  cout << "   # of forward events for Z's with a mass >80 and <100: " << F[1] << " +/- " << sqrt(F[1]) << endl;
  cout << "   # of backward events for Z's with a mass >80 and <100: " << B[1] << " +/- " << sqrt(B[1]) << endl;
  cout << "   # of forward events for Z's with a mass >100 and <120: " << F[2] << " +/- " << sqrt(F[2]) << endl;
  cout << "   # of backward events for Z's with a mass >100 and <120: " << B[2] << " +/- " << sqrt(B[2]) << endl;

  cout << "   The asymmetry for Z's with a mass >60 and <80: " << A[0] << endl;
  cout << "   The asymmetry for Z's with a mass >80 and <100: " << A[1] << endl;
  cout << "   The asymmetry for Z's with a mass >100 and <120: " << A[2] << endl;

  cout << "";
  cout << "   # of forward events for Z's (winter) >60 and <80: " << Fw[0] << " +/- " << sqrt(F[0]) << endl;
  cout << "   # of backward events for Z's (winter) >60 and <80: " << Bw[0] << " +/- " << sqrt(B[0]) << endl;
  cout << "   # of forward events for Z's (winter) >80 and <100: " << Fw[1] << " +/- " << sqrt(F[1]) << endl;
  cout << "   # of backward events for Z's (winter) >80 and <100: " << Bw[1] << " +/- " << sqrt(B[1]) << endl;
  cout << "   # of forward events for Z's (winter) >100 and <120: " << Fw[2] << " +/- " << sqrt(F[2]) << endl;
  cout << "   # of backward events for Z's (winter) >100 and <120: " << Bw[2] << " +/- " << sqrt(B[2]) << endl;
  cout << "";
  cout << "   The asymmetry (winter) >60 and <80: " << Aw[0] << endl;
  cout << "   The asymmetry (winter) >80 and <100: " << Aw[1] << endl;
  cout << "   The asymmetry (winter) >100 and <120: " << Aw[2] << endl;
  cout << "";
  cout << "   # of forward events for Z's (summer) >60 and <80: " << Fs[0] << " +/- " << sqrt(F[0]) << endl;
  cout << "   # of backward events for Z's (summer) >60 and <80: " << Bs[0] << " +/- " << sqrt(B[0]) << endl;
  cout << "   # of forward events for Z's (summer) >80 and <100: " << Fs[1] << " +/- " << sqrt(F[1]) << endl;
  cout << "   # of backward events for Z's (summer) >80 and <100: " << Bs[1] << " +/- " << sqrt(B[1]) << endl;
  cout << "   # of forward events for Z's (summer) >100 and <120: " << Fs[2] << " +/- " << sqrt(F[2]) << endl;
  cout << "   # of backward events for Z's (summer) >100 and <120: " << Bs[2] << " +/- " << sqrt(B[2]) << endl;
  cout << "";
  cout << "   The asymmetry (summer) >60 and <80: " << As[0] << endl;
  cout << "   The asymmetry (summer) >80 and <100: " << As[1] << endl;
  cout << "   The asymmetry (summer) >100 and <120: " << As[2] << endl;


  Zfile ->Write();
}




//_____________________________________________________________________________

bool StandardElectronQ(myelectron_t& myelectron, int i, mytrack_t& mytrack)
{

  bool ans = false;

  //with fidele==1, you don't need to put electron.detector==0 (CEM)
  bool fiducial = (myelectron.fidele[i]==1);  // |Xces|<21cm, 9<|Zces|<230cm, tower 9 is excluded.

  bool goodTrack = GoodTrackQ(mytrack, myelectron.trkindex[i]);

  bool highpt = ( (myelectron.etcorr[i] > 25.0) &&
                  (myelectron.pttrk[i] > 10.0) );

  if(!highpt)   return(ans);
 
  bool centralEleId = 
    (
     (myelectron.charge[i]*myelectron.delx[i] >-3.0) &&
     (myelectron.charge[i]*myelectron.delx[i] < 1.5) &&
     (fabs(myelectron.delz[i]) < 3.0)
     );

  bool centralEmId =
      (
       ( myelectron.hadem[i]<(0.055+0.00045*myelectron.e[i]) ) &&
       ( myelectron.lshr[i] < 0.2) &&
       ( myelectron.iso4[i]/myelectron.etcorr[i] < 0.10 ) &&
       ( myelectron.chi2strip[i] < 10.0) &&
       ((myelectron.ep[i] < 2.0)||(myelectron.etcorr[i] > 50.0))
       );

        ans = fiducial && goodTrack && highpt && centralEleId && centralEmId ;
  //        ans = fiducial && goodTrack && highpt && centralEmId ;
  
  return(ans);
}

bool LooseElectronQ(myelectron_t& myelectron, int i, mytrack_t& mytrack)
{

  bool ans = false;

  //with fidele==1, you don't need to put electron.detector==0 (CEM)
  bool fiducial = (myelectron.fidele[i]==1);  // |Xces|<21cm, 9<|Zces|<230cm, tower 9 is excluded.

  bool goodTrack = GoodTrackQ(mytrack, myelectron.trkindex[i]);

  bool highpt = ( (myelectron.etcorr[i] > 25.0) &&
                  (myelectron.pttrk[i] > 10.0) );

  if(!highpt)   return(ans);
 
  bool centralEleId = true;
   

  bool centralEmId =
      (
       ( myelectron.hadem[i]<(0.055+0.00045*myelectron.e[i]) ) &&
       ( myelectron.iso4[i]/myelectron.etcorr[i] < 0.10 )  );

        ans = fiducial && goodTrack && highpt && centralEleId && centralEmId ;
  //        ans = fiducial && goodTrack && highpt && centralEmId ;
  
  return(ans);
}


//_____________________________________________________________________________
bool GoodTrackQ(mytrack_t& mytrack, int i)
{

  if(i<0) return(false);

  bool ans = false;

  Int_t num_ax=0, num_axhits=0;
  Int_t num_st=0, num_sthits=0;

  //old method

  /*
    ans =

    (mytrack.naxialhits[i] > 24 && 
     mytrack.nstereohits[i] > 24 && 
     mytrack.naxialsl[i] >= 3 && 
     mytrack.nstereosl[i] >= 3 );

  */
  for (int j = 0; j < 8; j++) {
    Int_t num_hit_sl = (mytrack.slhitcount[i]>>j*4)&0xf;

//    cout << "check " << num_hit_sl << endl;

    if (fmod(j*1.0,2)>0) {
      num_axhits +=num_hit_sl; 
    } else {
      num_sthits +=num_hit_sl; 
    }

    if(num_hit_sl<7) continue;

    if (fmod(j*1.0,2)>0) {
      ++num_ax;}
    else {
      ++num_st;
    }
  }  

  /*
  if( (mytrack.naxialhits[i]!=num_axhits) ||
      (mytrack.nstereohits[i]!=num_sthits) ) {

    cout << "Problem in couting COT hit  " 
	 << setw(5) << mytrack.naxialhits[i]
	 << setw(5) << num_axhits
	 << setw(5) << mytrack.nstereohits[i]
	 << setw(5) << num_sthits << endl;

  }
  */

  ans = (num_ax>=3&&num_st>=3);
  ans &= (fabs(mytrack.z0[i])<60.);

  return(ans);

}
//_____________________________________________________________________________

Float_t invmass(Float_t et1, Float_t theta1, Float_t phi1, Float_t et2, Float_t theta2, Float_t phi2)
{

  //  if(theta1==0.0||theta2==0.0) {
  //    cout << "Warning in theta value " << endl;
  //  }

  Float_t E_1 = et1/sin(theta1);
  Float_t Ex_1 = E_1*sin(theta1)*cos(phi1);
  Float_t Ey_1 = E_1*sin(theta1)*sin(phi1);
  Float_t Ez_1 = E_1*cos(theta1);

  Float_t E_2 = et2/sin(theta2);
  Float_t Ex_2 = E_2*sin(theta2)*cos(phi2);
  Float_t Ey_2 = E_2*sin(theta2)*sin(phi2);
  Float_t Ez_2 = E_2*cos(theta2);
  
  Float_t mass  = sqrt(pow((E_1+E_2),2) - pow((Ex_1+Ex_2),2) 
		     - pow((Ey_1+Ey_2),2) - pow((Ez_1+Ez_2),2));
  
  return mass;

}