#define h99_cxx
#include "h99.h"
#include "TH2.h"
#include "TStyle.h"
#include "TCanvas.h"

void h99::Loop()
{
//   In a ROOT session, you can do:
//      Root > .L h99.C
//      Root > h99 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:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(i);  // read all branches
//by  b_branchname->GetEntry(i); //read only this branch

   if (fChain == 0) return;

 // call define histograms
   padx = 2;
   pady = 3;
   define_histos();
   TCanvas *canvas1 = blankcanvas(); //initial track
   

   Float_t etemp, VNU[4][2];
   Float_t Ptsq, Ppsq, Ptele, Ptmuo, Ptlep, Etae, Etam, Etahw, PpsqE, PpsqM, Etalep;
   
   //Hadronic W variables
   Float_t pthw, Mhw, Etaw, Pthwbos_true;

   //cout << "fChain->GetEntries()" <<fChain->GetEntries() << endl
   
   Int_t nentries = Int_t(fChain->GetEntriesFast());
     
 //cout <<"Events: "<<nentries<<endl;

   Int_t nbytes = 0, nb = 0;

   for (Int_t jentry=0; jentry<nentries;jentry++) {
      Int_t ientry = LoadTree(jentry); //in case of a TChain, ientry is the entry number in the current file
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // if (Cut(ientry) < 0) continue;

      //EVENT LOOP
      //find and store highest pt lepton in event
   
      Float_t pthigh(0.0);
      Int_t glw(0), ghw(0), nlh(0); 
      Int_t count1(0);
      Int_t count2(0);

      if ((Nele_unsmeared || Nmuo_unsmeared) ==0 ) { 

	continue;
      }


      //Ptele = sqrt(pow(Pxele_unsmeared[nlh-1],2) + pow(Pyele_unsmeared[nlh-1],2));
      // Ptmuo = sqrt(pow(Pxmuo_unsmeared[nlh-1],2) + pow(Pymuo_unsmeared[nlh-1],2));

      // Ptlep = Ptele + Ptmuo;

      for (Int_t i = 0; i<(Nele_unsmeared + Nmuo_unsmeared); i++){  

	if (abs(Kfele_unsmeared[i]) == 11 || abs(Kfmuo_unsmeared[i]) == 13){ 

	  if (sqrt(pow(Pxele_unsmeared[i],2) + pow(Pyele_unsmeared[i],2)) +  sqrt(pow(Pxmuo_unsmeared[i],2) +  pow(Pymuo_unsmeared[i],2)) > pthigh){

	    //if (sqrt(pow(Pxele_unsmeared[i],2) + pow(Pyele_unsmeared[i],2)) >  sqrt(pow(Pxmuo_unsmeared[i],2) +  pow(Pymuo_unsmeared[i],2))) {
   
	    //pthigh = sqrt(pow(Pxele_unsmeared[i],2) + pow(Pyele_unsmeared[i],2)); 
	    //} else { 
	    //pthigh = sqrt(pow(Pxmuo_unsmeared[i],2) + pow(Pymuo_unsmeared[i],2));
	    
	  
	    pthigh = sqrt(pow(Pxele_unsmeared[i],2) + pow(Pyele_unsmeared[i],2)) +  sqrt(pow(Pxmuo_unsmeared[i],2) +  pow(Pymuo_unsmeared[i],2));
	    
	    nlh = i+1;

	  }
	}
      }
   
          
	  
        if (nlh == 0) { 
	continue;
	}

      etemp = pow(Ptmiss_unsmeared,2);
      if (etemp <0){

	cout << "Event " << jentry << ": missing ET < 0" << endl;
	continue;
      }

      Ptmiss_unsmeared = sqrt(etemp);

      //      if (Nele_unsmeared == 0 ) {
      //cout << "Event " << jentry << ": Nele_unsmeared = 0" << endl;
      //continue;
       //}

      //if (Nmuo_unsmeared == 0 ) {
      //cout << "Event " << jentry << ": Nmuo_unsmeared = 0" << endl;
      // continue;
      //}

      Ptele = sqrt(pow(Pxele_unsmeared[nlh-1],2) + pow(Pyele_unsmeared[nlh-1],2));
      Ptmuo = sqrt(pow(Pxmuo_unsmeared[nlh-1],2) + pow(Pymuo_unsmeared[nlh-1],2));

      Ptlep = Ptele + Ptmuo;

	    ///////////////////calculate eta of lepton.
     
	    
      PpsqE = pow((sqrt(Ptele*Ptele + pow(Pzele_unsmeared[nlh-1],2)) + abs(Pzele_unsmeared[nlh-1])),2);

      if ((Ptele>0) || (Ptele<0)){
	Etae = 0.5*log(PpsqE/(Ptele*Ptele));}

      //make Etae have same sign as Pz
      if ((Pzele_unsmeared[nlh-1]<0) || (Pzele_unsmeared[nlh-1]>0)){
      Etae = Etae*(Pzele_unsmeared[nlh-1]/fabs(Pzele_unsmeared[nlh-1]));
      }
    
      PpsqM = pow((sqrt(Ptmuo*Ptmuo + pow(Pzmuo_unsmeared[nlh-1],2)) + abs(Pzmuo_unsmeared[nlh- 1])),2);

      if ((Ptmuo>0) || (Ptmuo<0)){
      Etam = 0.5*log(PpsqM/(Ptmuo*Ptmuo));
      }

      //make Etam have same sign as Pz
      if ((Pzmuo_unsmeared[nlh-1]<0) || (Pzmuo_unsmeared[nlh-1]>0)){
      Etam = Etam*(Pzmuo_unsmeared[nlh-1]/fabs(Pzmuo_unsmeared[nlh-1]));
      }

      Etalep = Etae + Etam;
      
      //-------------------------------------------------//

      //histo[0]->Fill(Ptele);
      histo[0]->Fill(Ptlep);
      //histo[1]->Fill(Ptmuo);
      //histo[2]->Fill(Etalep);
      histo[1]->Fill(Ptmiss_unsmeared);
      histo[2]->Fill(Etalep);


/*****************************hadronic W****************************/
	    //GENERATE HADRONIC W PLOTS;
   
         if (Nwbos_true == 0 ) {
        cout << "Event " << jentry << ": Nwbos_true = 0" << endl;
       continue;
       }


      Pthwbos_true = sqrt(pow(Pxwbos_true[Nwbos_true - 1],2) + pow(Pywbos_true[Nwbos_true - 1],2));

      Mhw =  sqrt(pow(Eewbos_true[Nwbos_true- 1],2) - pow(Pxwbos_true[Nwbos_true - 1],2) - pow(Pywbos_true[Nwbos_true - 1],2) - pow(Pzwbos_true[Nwbos_true - 1],2));

      Ptsq= pow(Pthwbos_true,2);
	    
      Ppsq = pow((sqrt(Ptsq + pow(Pzwbos_true[Nwbos_true - 1],2)) + abs(Pzwbos_true[Nwbos_true - 1])),2);

      Etahw = 0.5*log(Ppsq/Ptsq);

      //make Etaw have same sign as Pz

      Etahw = Etahw*(Pzwbos_true[Nwbos_true - 1]/fabs(Pzwbos_true[Nwbos_true - 1]));


      //Fill hadronic W candidate plots (highest Et jet in event)

      histo[3]->Fill(Pthwbos_true);
      histo[4]->Fill(Mhw);
      histo[5]->Fill(Etahw);
      


      //MAKING HADRONIC W SELECTION;
      
      if (Pthwbos_true > 320.0) {
	if ((Mhw>70.0) && (Mhw<100.0) && (abs(Etahw)<4)) {
	  ghw = ghw + 1; }

	if (ghw>1) {cout << "More than one hadronic W candidates\n" << endl;}
      }

      if (ghw=0) {
	
      continue;
      }

      //Now have at least one hadronic W candidate, and have passed trigger and leptonic W cuts);

      //CANNOT DO Ycut BECAUSE IT ISN'T IMPLEMENTED IN ATLFAST//
	 
      
      //top quark veto, assuming highest Pt subjet had W candidate
      
      Float_t tc[4], wlep[4], tcm;
      
      for (i =0; i<Nwbos_true; i++){
	
	Mhw =  sqrt(pow(Eewbos_true[Nwbos_true- 1],2) - pow(Pxwbos_true[Nwbos_true - 1],2) - pow(Pywbos_true[Nwbos_true - 1],2) - pow(Pzwbos_true[Nwbos_true - 1],2));

	if (Mhw>160.0 && Mhw<200.0){
	  cout << "Event " << jentry << ": top quark veto" << endl;
	  continue;
	/}

	//Now combine jet with leptonic W

	tc[0]= Pxwbos_true[i] + wlep[0];
	tc[1]= Pywbos_true[i];
	tc[2]= Pzwbos_true[i];
	tc[3]= Eewbos_true[i];

	//tcm = sqrt((tc[3]*tc[3]) -(tc[0]*tc[0]) -(tc[1]*tc[1]) -(tc[2]*[2]));

	//histo[8]->Fill(tcm);

	//if ((tcm>60.0) && (tcm<240.0)){
	//continue;
	//}

      }
/************************** END OF HADRONIC W *******************************/


/*********************** lEPTONIC W *****************************************/
   
      //      search for leptonic W decay
      // reconstruct W from highest pt lepton and ptmiss

      // FIX FOR PHI
    
       Float_t Philep; //...........................//

      
       Philep = atan2(Pyele_unsmeared[Nele_unsmeared] + Pymuo_unsmeared[Nmuo_unsmeared],Pxele_unsmeared[Nele_unsmeared] + Pxmuo_unsmeared[Nmuo_unsmeared]);

       //first do Pt and rapidity cut on highest lepton
       // wtest (function)


      Float_t Pt, Phi, mass(80.1);
      Float_t nx, ny, Ptnu, t, lx, ly, lz, le;
      Float_t EtaNu1, EtaNu2, a, b, c, testz1, testz2, teste1, teste2;
      Float_t Ptlw, Mlw, Ptsqlw, Ptsqlw, Ppsqlw, Etalw;

      if (Nnue_true == 0 ) {
        cout << "Event " << jentry << ": Nnue_true = 0" << endl;
       continue;
       }

      nx = -Pxnue_true[Nnue_true];//I've included no of neutrinos
      ny = -Pynue_true[Nnue_true];
      Ptnu = sqrt(nx*nx + ny*ny);
      lx = Ptlep*cos(Philep);
      ly = Ptlep*sin(Philep); 
     
      t = exp(Etalep);      
      
      if (t==0) {
	  continue;
	}

      lz = Ptlep*(t*t-1)/(2.0*t);
      le = Ptlep*(1+t*t)/(2.0*t);

      a = 4.0*(le*le - lz*lz);
      b = 4.0*mass*mass*lz + 8.0*nx*lx*lz + 8.0*ny*ly*lz;
      c = 4.0*(nx*nx + ny*ny)*le*le - 4.0*nx*lx*nx*lx - 4.0*ny*ly*ny*ly - 4.0*mass*mass*nx*lx - 4.0*mass*mass*ny*ly - 8.0*nx*lx*ny*ly - mass*mass*mass*mass;

      if (b*b < 4*a*c){
	wlep[3]=-999.0;
	cout << "Event " << jentry << ": b*b < 4ac " << endl;
	continue;
      }

      wlep[0] = lx + nx;
      wlep[1] = ly + ny;

      if (a == 0.0){	
	cout << "Event " << jentry << ": Error: a = 0" << endl;
	continue;
      }

      testz1 = (b + sqrt(b*b - 4.0*a*c))/(2.0*a);
      testz2 = (b - sqrt(b*b - 4.0*a*c))/(2.0*a);
      
      teste1 = sqrt(nx*nx + ny*ny + testz1*testz1);
      teste2 = sqrt(nx*nx + ny*ny + testz2*testz2);

      //Calculate rapidity of neutrino

      if ((teste1 - testz1) == 0.0) {
	cout << "Event " << jentry << ": EtaNu1 is infinite " << endl;
	continue;
      }


      EtaNu1 = 0.5*log((teste1 + testz1)/(teste1 - testz1));

      if ((teste2 - testz2) == 0.0) {
       
	cout << "Event " << jentry << ": EtaNu2 is infinite " << endl;
	continue;
      }

      EtaNu2 = 0.5*log((teste2 + testz2)/(teste2 - testz2));

      //EtaN = sqrt(pow(EtaNu[0][Nnue_true],2) + pow(EtaNu[1][Nnue_true],2));

       //     histo[6]->Fill(EtaNu);

      VNU[0][0] = nx;
      VNU[0][1] = nx;
      VNU[1][0] = ny;
      VNU[1][1] = ny;
      
      if (testz1 > testz2) {

	VNU[2][0] = testz1;
	VNU[2][1] = testz2;
	VNU[3][0] = teste1;
	VNU[3][1] = teste2;
      
      } else {

	VNU[2][0] = testz2;
	VNU[2][1] = testz1;
	VNU[3][0] = teste2;
	VNU[3][1] = teste1;

      }

      wlep[3] = VNU[3][0] + le;


      //Jon bug fix 06/11

      wlep[2] = VNU[2][0] + lz;

      // max energy constraint

      if (wlep[3] > 7000.0) {
	wlep[3] = -999.0;
	  }

      // check error

      if (wlep[3] < -998.0) {continue;}


      // check for a good leptonic W
      // If found, carry on

      
      Ptlw = sqrt(pow(Pxwbos_true[Nwbos_true],2) + pow(Pywbos_true[Nwbos_true],2));

      
      Mlw = sqrt(pow(Eewbos_true[Nwbos_true],2) - pow(Pxwbos_true[Nwbos_true],2) -  pow(Pywbos_true[Nwbos_true],2) - pow(Pzwbos_true[Nwbos_true],2));

      histo[6]->Fill(Ptlw);
      histo[7]->Fill(Mlw);


      if (Ptlw >320.0){glw=1;}

      if (glw=0){continue;}

      //Pseudorapidity of leptonic W for later

      Ptsqlw = sqrt(wlep[0]*wlep[0] + wlep[1]*wlep[1]);

      Ppsqlw = pow((sqrt(Ptsqlw + wlep[2]*wlep[2]) + abs(wlep[2])),2);

      Etalw = 0.5*log(Ppsqlw/Ptsqlw);

      Etalw = 0.5*(wlep[2]/fabs(wlep[2]));



      //END OF LEPTONIC W//


      //................TAG JETS.....................//

      //Identify WW rapidity region

      Float_t Wetah, Wetal, Ptjet, Etajet, PpsqJ;
      Int_t tagf(0);
      Int_t tagb(0);

      Float_t tagfpt(0.0);
      Float_t tagbpt(0.0);
      Int_t tageta(0);
      Float_t tagfeta(100.0);
      Float_t tagbeta(-100.0);
      Int_t high(0);
      Int_t low(0);
      Float_t mjveto(0.0);

      if (Etahw > Etalw){

	Wetah = Etahw;
	Wetal = Etalw;
      } else

	{
	  Wetah = Etalw;
	  Wetal = Etahw;
	}


      if (Njet_cone_unsmeared == 0 ) {
        cout << "Event " << jentry << ": Njet_cone_unsmeared = 0" << endl;
	continue;
       }

     
      Ptjet = sqrt(pow(Pxjet_cone_unsmeared[Njet_cone_unsmeared-1],2) + pow(Pyjet_cone_unsmeared[Njet_cone_unsmeared-1],2));
     

	    //calculate eta of jet//.
   
      PpsqJ = pow((sqrt(Ptjet*Ptjet + pow(Pzjet_cone_unsmeared[Njet_cone_unsmeared-1],2)) + abs(Pzjet_cone_unsmeared[Njet_cone_unsmeared-1])),2);

      if ((Ptjet>0) || (Ptjet<0)){
	Etajet = 0.5*log(PpsqJ/(Ptjet*Ptjet));
      }



      //make Etajet have same sign as Pz
      
      if ((Pzjet_cone_unsmeared[Njet_cone_unsmeared-1]<0) || (Pzjet_cone_unsmeared[Njet_cone_unsmeared-1]>0)){
    
	Etajet = Etajet*(Pzjet_cone_unsmeared[Njet_cone_unsmeared-1]/fabs(Pzjet_cone_unsmeared[Njet_cone_unsmeared-1]));
      }
    

      //Find highest ET jets forward and backward of the W rapidity region


      for (Int_t i=0; i<Njet_cone_unsmeared-1; i++){

	if ((Ptjet[i]>20) && abs(Etajet[i]<4.5)){

	  if (Etajet[i] > Wetah + 0.01){

	    if (Ptjet[i]>tagfpt){
	      
	      tagfpt = Ptjet[i];
	      high = i+1;
	      tagfeta = Etajet[i];
	    }
	  }

	  if (etajet[i]<Wetal - 0.01){

	    if (Ptjet[i]>tagbpt){

	      tagbpt=Ptjet[i];
	      low=i+1;
	      tagbeta=Etajet[i];
	    }
	  }
	}
     
	//minijet veto

	if ((Ptjet[i]>20.0) && (abs(Etajet[i]<2))){

	  mjveto = mjveto + 1.0;
	}
      }

      histo[8]->Fill(mjveto);
      histo[9]->Fill(tagbeta);
      histo[10]->Fill(tagfeta);

//Now check that highest ET tag jet has abs(eta) > 2.0

      if ((high == 0) || (low == 0)){

	continue;
      }

      if ((tagfeta < -99.0) || (tagfeta > -99.0)){

	cout <<"bad tag eta" << endl;}
   
      if (tagfpt > tagbpt) {

	if (abs(tagfeta) > 2.0) {

	  tageta = 1;}
       
	if (abs(tagbeta > 2.0)) {

	  tageta = 1;}
      }
   
      if (tageta = 0) {

	continue;
      }



      Float_t PxWW, PxWW, PyWW, MWW, m[4]

      //Now we have good tag jets, calculate Pt of WW + tagjet system

      PxWW = Pxwbos_true[Nwbos_true -1] + wlep[0] + Pyjet[high-1] + Pyjet[low-1];
      PyWW =  Pywbos_true[Nwbos_true -1] + wlep[1] + Pzjet[high-1] + Pzjet[low-1];
      PtWW = sqrt(pow(PxWW, 2) + pow(PyWW,2));

      ////histo[]->Fill(PtWW);
	

      //Now do hard Pt cut

      if (PtWW < 50.0){

	continue;
	count1++
      }

      if (mjvet0 > 1.0){

	continue;
	count2++;
      }


      //WW Mass

      m[0] = Pxwbos_true[Nwbos_true-1] + wlep[0];
      m[1] = Pywbos_true[Nwbos_true-1] + wlep[1];
      m[2] = Pzwbos_true[Nwbos_true-1] + wlep[2];
      m[3] = Eewbos_true[Nwbos_true-1] + wlep[3];

      MWW = sqrt(pow(m[3],2) - pow(m[0],2) - pow(m[1],2) - pow(m[2],2));

      //histo[]->Fill(MWW);

 




   } //event Loop

// plot the histograms
   int colour = 8;

   for (int j=0; j<1; j++) {
      gStyle->SetOptStat(0);
      gStyle->SetOptStat(10111110);

     for(int i=0; i<11; i++) {
       pad[j]->cd(i+1);
       if (i==0 || i==1 || i==3){
	 gPad->SetLogy();
       }
       // if((i==4 ) && j==0) {
       ////currpad = (TPad*)gPad;
       //currpad->SetLogx();
       // }
       currpad = (TPad*)gPad;
       if (colour) histo[i+11*j]->SetFillColor(colour);

       histo[i+11*j]->DrawCopy();

       //  histo[i+6*j]->Write();
}
     pad[j]->Update();
   }
}
   

void h99::define_histos()
{

histo[0]  = new TH1F("Ptlep","Pt of leptons",100, 0, 100);
  //histo[0] = new TH1F("Ptele","Pt of electron",100, 0, 50);
  //histo[1] = new TH1F("Ptmuo","Pt of muon", 100, 0, 50);
// histo[2] = new TH1F("Etalep", "Eta of lepton", 50, -5, 5);
histo[1] = new TH1F("Ptmiss_unsmeared","Missing energy", 100, 0, 700);
histo[2] = new TH1F("Etalep","Eta of lepton", 50, -4, 4);
histo[3] = new TH1F("Pthwbos_true","Pt of hadronic W", 1600, 0, 800);
histo[4] = new TH1F("Mhw","Mass of hadronic W", 40, 0, 200);
histo[5] = new TH1F("Etaw","Eta of hadronic W", 50, -4, 4);
histo[6] = new TH1F("Ptlw", "Pt of leptonic W", 200, 0, 1000);
histo[7] = new TH1F("Mlw", "Mass of leptonic W", 40, 0, 500);
//histo[8] = new TH1F("tcm", "Top mass", 50, 0, 200);
histo[8] = new TH1F("mjveto", "Minijet veto", 10, 0, 10);
histo[9] = new TH1F("tagbeta", "Eta of b tag jet", 10, 0, 10);
histo[10] = new TH1F("tagfeta", "Eta of f tag jet", 10, 0, 10);

}