#define Higgs_cxx
#include "Higgs.h"
#include "TH2.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TLorentzVector.h"
#include <iostream>
using std::cout; using std::endl;

void Higgs::Loop()
{
//   In a ROOT session, you can do:
//      Root > .L Higgs.C
//      Root > Higgs 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 where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch

  TH1F* pt_el  = new TH1F("pt_el","p_T of electrons",100,0.,200.);
  TH1F* pt_mu  = new TH1F("pt_mu","p_T of muons",    100,0.,200.);
  TH1F* eta_el = new TH1F("eta_el","eta of electrons",50,-5.,5.);
  TH1F* eta_mu = new TH1F("eta_mu","eta of muons",    50,-5.,5.);
  
  TH1F* mass_H     = new TH1F("mass_H","MC Higgs mass",           100,0.,500.);
  TH1F* mass_Heemm = new TH1F("mass_Heemm","e+e-mu+mu- inv. mass",100,0.,500.);
  TH1F* mass_Z     = new TH1F("mass_Z","MC Z mass",               100,0.,200.);
  TH1F* mass_Zee   = new TH1F("mass_Zee","e+e- inv. mass",        100,0.,200.);
  TH1F* mass_Zmm   = new TH1F("mass_Zmm","mu+mu- inv. mass",      100,0.,200.);
  if (fChain == 0) return;

  //   Int_t nentries = Int_t(fChain->GetEntriesFast()); // ???
  Int_t nentries = Int_t(fChain->GetEntries());


  
  // ... Particle MC IDs; PDG=HEPEVT standard: electron/positron = -11/+11 ...
  const int ID_el = -11, ID_mu = -13, ID_Z = 23, ID_Higgs = 25;
  
  // ... and Line numbers in HEPEVT block for e+/e-, mu+/mu-, etc. ...
  int L_eln, L_elp, L_mun, L_mup;
  int L_H, L_Z1, L_Z2;             // L_Zee, L_Zmm;
  
  // ... Set some counters to 0 ...
  int nev[20] = {0};
  
  Int_t nbytes = 0, nb = 0;


  for (Int_t jentry=0; jentry<nentries;jentry++) {
    Int_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    // if (Cut(ientry) < 0) continue;

    if(jentry%1000  ==  0) cout<<" Ievent = "<<jentry<<endl; // << "!";

    ++nev[0];
    L_eln = -1; L_elp = -1; L_mun = -1; L_mup = -1;
    // ... Loop over particles
    for(int i=0; i<Nhep; i++) {
      // ... Find a Higgs 
      Int_t status  = Jsmhep[i]/16000000;            // status code of particle 
      if(Idhep[i] == ID_Higgs && status == 2){
	Int_t last    = Jsdhep[i]%4000 -1 ;            // index of last daughter particle 
	Int_t first   = (Jsdhep[i]-last)/4000-1;       // index of first daughter particle 
	L_H = i;
	mass_H->Fill(Phep[i][4]);
	L_Z1 = first;
	L_Z2 = last;
	// ... Make sure this is Higgs -> XX two-body decay
	// ... Make sure daughter ptcls. are Z's
	if(Idhep[L_Z1] == ID_Z && Idhep[L_Z2] == ID_Z )
	  {
	    mass_Z->Fill(Phep[L_Z1][4]);
	    mass_Z->Fill(Phep[L_Z2][4]);
	    ++nev[1];
	    // ... Make sure Z1 daughters exist
	    last    = Jsdhep[L_Z1]%4000 -1 ;            // index of last daughter particle 
	    first   = (Jsdhep[L_Z1]-last)/4000-1;       // index of first daughter particle 
	    if(first>0 && last>0) {
	      ++nev[2];
	      for(int j=first; j<=last; j++) {
		if(Idhep[j] ==  ID_el) L_eln = j;
		if(Idhep[j] == -ID_el) L_elp = j;
		if(Idhep[j] ==  ID_mu) L_mun = j;
		if(Idhep[j] == -ID_mu) L_mup = j;
	      }
	    }
	    // ... Make sure Z2 daughters exist
	    last    = Jsdhep[L_Z2]%4000 -1 ;            // index of last daughter particle 
	    first   = (Jsdhep[L_Z2]-last)/4000-1;       // index of first daughter particle 
	    if(first>0 && last>0) {
	      ++nev[3];
	      for(int j=first; j<=last; j++) {
		if(Idhep[j] ==  ID_el) L_eln = j;
		if(Idhep[j] == -ID_el) L_elp = j;
		if(Idhep[j] ==  ID_mu) L_mun = j;
		if(Idhep[j] == -ID_mu) L_mup = j;
	      }
	    }
	  }
      }
    } // End Loop over particles/i


    // ... Are all four leptons found ?
    if(L_eln<=0 || L_elp<=0 || L_mun<=0 || L_mup<=0) continue;

    // ... Fill in the vectors (electrons/positrons first), histos, ...
    TLorentzVector mcpL[4];
    TLorentzVector Zee, Zmm, Higgs;

    mcpL[0].SetPxPyPzE( double(Phep[L_eln][0]),double(Phep[L_eln][1]),
                        double(Phep[L_eln][2]),double(Phep[L_eln][3]) );
    mcpL[1].SetPxPyPzE( double(Phep[L_elp][0]),double(Phep[L_elp][1]),
                        double(Phep[L_elp][2]),double(Phep[L_elp][3]) );
    Zee = mcpL[0] + mcpL[1];

    mcpL[2].SetPxPyPzE( double(Phep[L_mun][0]),double(Phep[L_mun][1]),
                        double(Phep[L_mun][2]),double(Phep[L_mun][3]) );
    mcpL[3].SetPxPyPzE( double(Phep[L_mup][0]),double(Phep[L_mup][1]),
                        double(Phep[L_mup][2]),double(Phep[L_mup][3]) );
    Zmm = mcpL[2] + mcpL[3];

    Higgs = Zee + Zmm;

    Float_t pt, eta, mass;

    pt = mcpL[0].Pt(); pt_el->Fill(pt); pt = mcpL[1].Pt(); pt_el->Fill(pt);
    pt = mcpL[2].Pt(); pt_mu->Fill(pt); pt = mcpL[3].Pt(); pt_mu->Fill(pt);

    eta=mcpL[0].Eta(); eta_el->Fill(eta); eta=mcpL[1].Eta(); eta_el->Fill(eta);
    eta=mcpL[2].Eta(); eta_mu->Fill(eta); eta=mcpL[3].Eta(); eta_mu->Fill(eta);

    mass  = Zee.M();   mass_Zee->Fill(mass);
    mass  = Zmm.M();   mass_Zmm->Fill(mass);
    mass  = Higgs.M(); mass_Heemm->Fill(mass);

  } // End Loop over events/jentry


}