#define dataAnalysis_cxx
#include "dataAnalysis.h"
#include <TH2.h>
#include <TH1.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <vector>
#include<TMath.h>
#include<iostream>
#include<cmath>
using namespace std;

 void dataAnalysis::Loop()
{
//   In a ROOT session, you can do:
//      root> .L dataAnalysis.C 
//      root> dataAnalysis t
//      root> t.GetEntry(12); // Fill t dataAnalysis 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
   if (fChain == 0) return;
   
   
                                                   
       //TFile *f = new TFile("/home/ahmed/Downloads/MG5_aMC_v3_4_2/pp_h90BP1_bb/Events/run_01/h90BP1_bb.root","READ");
      // TFile *f = new TFile("/home/ahmed/Downloads/MG5_aMC_v3_4_2/pp_h90BP2_bb/Events/run_01/h90BP2_bb.root","READ");
      
      
      //TFile *f = new TFile("/home/ahmed/Downloads/MG5_aMC_v3_4_2/pp_h95BP1_bb/Events/run_01/h95BP1_bb.root","READ");
     // TFile *f = new TFile("/home/ahmed/Downloads/MG5_aMC_v3_4_2/pp_h95BP3_bb/Events/run_01/h95BP3.root","READ");
                                                     
      //  TFile *f = new TFile("/home/ahmed/Downloads/MG5_aMC_v3_4_2/pp_h96BP1_bb/Events/run_01/h96BP1_bb.root","READ");
     // TFile *f = new TFile("/home/ahmed/Downloads/MG5_aMC_v3_4_2/pp_h96BP2_bb/Events/run_01/h96BP2_bb.root","READ");
     //TFile *f = new TFile("/home/ahmed/Downloads/MG5_aMC_v3_4_2/pp_h96BP3_bb/Events/run_01/h96BP3_bb.root","READ");
    // TFile *f = new TFile("/home/ahmed/Downloads/MG5_aMC_v3_4_2/pp_h96BP6_bb/Events/run_01/h96BP6_bb.root","READ");                                                
                                                    
                                              
                                // **************************************                   
                                                                                                                 
                                                  
                                    
     
     // TFile *f = new TFile("/home/ahmed/Downloads/MG5_aMC_v3_4_2/pp_h_bb/Events/run_01/h_bb.root","READ");
    //TFile *f = new TFile("/home/ahmed/Downloads/MG5_aMC_v3_4_2/pp_z_bb/Events/run_01/z_bb.root","READ");
   // TFile *f = new TFile("/home/ahmed/Downloads/MG5_aMC_v3_4_2/pp_bb/Events/run_01/bb.root","READ");
                                                      //***************
    
    
                                                                                                            
           
       
       // **************************************                           
                                   
   
   // weight calculation 
      
        float cross_section = 2.881e+03  ;   //(in fb)- pp_h90BP1_bb  (14 TeV)
      // float cross_section =   2.718e+03  ;    //(in fb)- pp_h90BP2_bb  (14 TeV)
 
 
        //float cross_section = 3.959e+03   ;   //(in fb)- pp_95BP1_bb (14TeV)
       // float cross_section =  2.54e+03 ; //(in fb)- pp_95BP3_bb (14 TeV)
      
      
       //float cross_section = 4.122e+03   ; //(in fb)- pp_h96BP1_bb  (14 TeV)
       //float cross_section = 4.181e+03    ; //(in fb)- pp_h96BP2_bb  (14 TeV)
       //float cross_section =  2.64e+03    ; //(in fb)- pp_h96BP3_bb  (14 TeV)
      // float cross_section =   4.031e+03  ; //(in fb)- pp_h96BP6_bb  (14 TeV)  
      
      
                             //*****************************
         
     
       //float cross_section =  255.2 ; //(in fb)- pp_h_bb  (14 TeV)
      // float cross_section =  7.297e+6; //(in fb)- pp_z_bb (14TeV)
      //float cross_section = 6.039e+09; //(in fb)- pp_bb (14 TeV)
             
              //*****************************
        
            
               
   
  
       
    
      float N_generation= 1.e+06;
      float lum_mc = N_generation/cross_section;
      float lum_Data =300.0;
      float wt= lum_Data /lum_mc; 

     
     
       //**************************************************************                                 
    
    
       // define histogram
     
                               
        TH1D *b1_pt = new TH1D("b1_pt", " b1_pt distribution ",3000 ,0. , 3000.);
          b1_pt ->GetXaxis()->SetTitle("pt_{b1 }[Gev/c] ");
          b1_pt->GetYaxis()->SetTitle("Events");
          
          TH1D *b1_eta = new TH1D("b1_eta", "eta1_Eta distribution ",1000 ,-20.,20. );
          b1_eta ->GetXaxis()->SetTitle("eta_{b1 }");
          b1_eta->GetYaxis()->SetTitle("Events");
                                   
          TH1D *b2_pt = new TH1D("b2_pt", " b2_pt distribution ",3000 ,0. , 3000.);
          b2_pt ->GetXaxis()->SetTitle("pt_{b2 [Gev/c]}");
          b2_pt->GetYaxis()->SetTitle("Events");
          
           TH1D *b2_eta = new TH1D("b2_eta", "eta2_Eta distribution ",1000 ,-20.,20. );
          b2_eta ->GetXaxis()->SetTitle("eta_{b2 }");
          b2_eta->GetYaxis()->SetTitle("Events");
       
       
         TH1D *h_mass = new TH1D("h_mass", " invariant mass distribution ",3000 ,0. , 3000.);
          h_mass ->GetXaxis()->SetTitle("M_{b1 & b2} [Gev]");
          h_mass->GetYaxis()->SetTitle("Events");
         
         
          TH1D *h_pt = new TH1D("h_pt", " h_pt distribution ",3000 ,0. , 3000.);
          h_pt ->GetXaxis()->SetTitle("pt_{b1 & b2 }[Gev/c] ");
          h_pt->GetYaxis()->SetTitle("Events");
                                       
          TH1D *h_eta = new TH1D("h_eta", "eta distribution ",1000 ,-20.,20.  );
          h_eta ->GetXaxis()->SetTitle("eta_{b1 & b2}");
          h_eta->GetYaxis()->SetTitle("Events");
       
        TH1D *DR_b = new TH1D("DR_b", " DR_b distribution ",1000 ,0 , 40.);
         DR_b ->GetXaxis()->SetTitle("DR_b_{b1 & b2} ");
         DR_b->GetYaxis()->SetTitle("Events");
       
                   //************************************************************
                   
          TH1D *hfJet_mass = new TH1D("hfJet_mass", " hfJet_mass distribution ",3000 ,0. , 3000.);
         hfJet_mass->GetXaxis()->SetTitle("M[Gev]");
         hfJet_mass->GetYaxis()->SetTitle("Events");
                   
                
            TH1D *  hfJet_eta= new TH1D("  hfJet_eta", " hfJet_eta distribution ",1000 ,-20.,20. );
             hfJet_eta ->GetXaxis()->SetTitle("eta");
            hfJet_eta->GetYaxis()->SetTitle("Events");
                                     
                   
          TH1D *  hfJet_pt= new TH1D("  hfJet_pt", "  hfJet_pt distribution ",3000 ,0. , 3000.);
          hfJet_pt->GetXaxis()->SetTitle("pt[Gev/c] ");
          hfJet_pt->GetYaxis()->SetTitle("Events");
                   
                   
         
           TH1D *fjet_deltaEta = new TH1D("fjet_deltaEta", "fjet_deltaEta distribution ",1000 ,-20.,20. );
          fjet_deltaEta ->GetXaxis()->SetTitle("deltaEta");
          fjet_deltaEta->GetYaxis()->SetTitle("Events");
              
               
           TH1D *fjet_deltaPhi = new TH1D("fjet_deltaPhi", "fjet_deltaPhi distribution ",1000 ,-20.,20. );
          fjet_deltaPhi ->GetXaxis()->SetTitle("deltaPhi");
          fjet_deltaPhi->GetYaxis()->SetTitle("Events");
                         
                         
           TH1D *fjet_EhadOverEem = new TH1D(" fjet_EhadOverEem", " fjet_EhadOverEem distribution ",1000 ,-20.,20. );
           fjet_EhadOverEem ->GetXaxis()->SetTitle("fjet_EhadOverEem");
          fjet_EhadOverEem->GetYaxis()->SetTitle("Events"); 
                  
                          
          TH1D *missingET_MET= new TH1D("missingET_MET", " missingET_MET distribution ",3000 ,0. , 3000.);
         missingET_MET ->GetXaxis()->SetTitle("missingET_MET[Gev]");
          missingET_MET->GetYaxis()->SetTitle("Events");
                      
         
         
           TH1D *missingET_phi = new TH1D("missingET_phi", "missingET_phi distribution ",1000 ,-20.,20. );
          missingET_phi->GetXaxis()->SetTitle("missingET_phi");
          missingET_phi->GetYaxis()->SetTitle("Events");
                         
          
          
          TH1D *scalarHt = new TH1D("scalarHt", " scalarHt distribution ",3000 ,0. , 3000.);
          scalarHt ->GetXaxis()->SetTitle("scalarHt[Gev]");
         scalarHt->GetYaxis()->SetTitle("Events");
                     
                          
                          
         
     Long64_t nentries = fChain->GetEntriesFast();
     Long64_t nbytes = 0, nb = 0;
    for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // if (Cut(ientry) < 0) continue;
     
     vector <int> b_quarkIndex,b_fatIndex;
         
          
			  
           	for (int i=0;i<Particle_size;i++){           	        	
 	        
 	        if(fabs(Particle_PID[i])!=5)continue;
 	         b_quarkIndex.push_back(i);  
 	               	                          	     	         	            
 	               	         	            }  
  
  
  	       	for (int i=0;i<Particle_size;i++){           	        	
 	        
 	          if(fabs(Particle_PID[i])!=5)continue;
 	          b_fatIndex.push_back(i);
 	               	         	                         	        
 	               	         	            }  
  
  
        
       
         TLorentzVector  b1, b2, h;
         TLorentzVector  bf1,bf2,hf;
         
                      
			 b1.SetPtEtaPhiE(Particle_PT[b_quarkIndex[0]],Particle_Eta[b_quarkIndex[0]],Particle_Phi[b_quarkIndex[0]],Particle_E[b_quarkIndex[0]]);
			 b2.SetPtEtaPhiE(Particle_PT[b_quarkIndex[1]],Particle_Eta[b_quarkIndex[1]],Particle_Phi[b_quarkIndex[1]],Particle_E[b_quarkIndex[1]]);
        
        
             bf1.SetPtEtaPhiE(Particle_PT[b_fatIndex[0]],Particle_Eta[b_fatIndex[0]],Particle_Phi[b_fatIndex[0]],Particle_E[b_fatIndex[0]]);
			 bf2.SetPtEtaPhiE(Particle_PT[b_fatIndex[1]],Particle_Eta[b_fatIndex[1]],Particle_Phi[b_fatIndex[1]],Particle_E[b_fatIndex[1]]);
               
			 
          float  dR_b, dEta_fJet,dPhi_fJet;     
         h= b1+b2;
         hf=bf1+bf2;
         
       dR_b= TMath::Sqrt(TMath::Power((Particle_Eta[b_quarkIndex[0]] - Particle_Eta[b_quarkIndex[1]]), 2) + TMath::Power((Particle_Phi[b_quarkIndex[0]] - Particle_Phi[b_quarkIndex[1]]), 2));
      
       dEta_fJet = (Particle_Eta[b_fatIndex[0]] - Particle_Eta[b_fatIndex[1]]);
       dPhi_fJet =(Particle_Phi[b_fatIndex[0]] - Particle_Phi[b_fatIndex[1]]);
               
               //********************************************
               
                // calculate Energy components -> caluclate EhadoverEem 
         float energy_hadronic = hf.E();
        float  energy_bf1 = bf1.E();
         float  energy_bf2 = bf2.E();
			 
             // Calculate total hadronic and electromagnetic energies
     
       float hadronicEnergy = energy_hadronic;
       float electromagneticEnergy = energy_bf1 + energy_bf2;
    
    // Calculate the ratio of hadronic energy to electromagnetic energy
   
    float ratio = electromagneticEnergy != 0 ? hadronicEnergy / electromagneticEnergy : 0;
              
                
              //**********************************
              
        // Calculate the total transverse momentum of visible particles (bf1, bf2) -> calculte MissingET 
       float px = bf1.Px() + bf2.Py() ;
      float py = bf1.Py() + bf2.Py() ;

     // Calculate the magnitude of the missing transverse energy
       float met = TMath::Sqrt (px *px +py *py);

     // Calculate the phi of the missing transverse energy
     float  met_phi = TMath::ATan2(py, px);      
              
           //***************************************   
          //calculate  ScalarHT, denoted as _HT, typically refers to the scalar sum of transverse momenta (pT) of all visible objects in an event
              
              // Calculate transverse momentum components
              
      float pt_bf1 = bf1.Pt();
     float pt_bf2 = bf2.Pt();
     float pt_hf = hf.Pt();

    // Calculate ScalarHT
    float scalarht_ht = pt_bf1 + pt_bf2 + pt_hf;
        
          

              
    b1_pt->Fill(Particle_PT[b_quarkIndex[0]], wt);
    b1_eta->Fill(Particle_Eta[b_quarkIndex[0]], wt);
    b2_pt->Fill(Particle_PT[b_quarkIndex[1]], wt);  
    b2_eta->Fill(Particle_Eta[b_quarkIndex[1]], wt); 
    h_mass->Fill(h.Mt(), wt);
    h_pt->Fill(h.Pt(), wt);
    h_eta->Fill(h.Eta(), wt);
    DR_b ->Fill( dR_b,wt);
  

    
    
     
   hfJet_mass->Fill(hf.Mt(), wt);
    hfJet_pt->Fill(hf.Pt(), wt);
   hfJet_eta->Fill(hf.Eta(), wt);
    
   
    fjet_deltaEta->Fill(dEta_fJet, wt);
    fjet_deltaPhi->Fill(dPhi_fJet, wt);
    fjet_EhadOverEem->Fill(ratio, wt);
    
    
    missingET_MET->Fill(met, wt); 
    missingET_phi->Fill(met_phi, wt);

     scalarHt->Fill(scalarht_ht, wt);
     
   
  
}//end of event loop.

//**************outputfile
       
                                //*********************************
            TFile *outputfile = new TFile("h90BP1_bb.root", "RECREATE");
           //TFile *outputfile = new TFile("h90BP2_bb.root", "RECREATE");
         
         //  TFile *outputfile = new TFile("h95BP1_bb.root", "RECREATE");
         // TFile *outputfile = new TFile("h95BP3_bb.root", "RECREATE");
    
    
       //TFile *outputfile = new TFile("h96BP1_bb.root", "RECREATE");
      //TFile *outputfile = new TFile("h96BP2_bb.root", "RECREATE");                 
      //TFile *outputfile = new TFile("h96BP3_bb.root", "RECREATE");
     // TFile *outputfile = new TFile("h96BP6_bb.root", "RECREATE");
 
                             //**********************
    
    //TFile *outputfile = new TFile("h_bb.root", "RECREATE");
    //  TFile *outputfile = new TFile("z_bb.root", "RECREATE");
    //TFile *outputfile = new TFile("bb.root", "RECREATE");
 
 
  b1_pt-> Write();
  b1_eta->Write();
  b2_pt-> Write();
  b2_eta->Write();
  h_mass->Write();
  h_pt->Write();
  h_eta->Write();
  DR_b->Write();

 
 
 
  hfJet_mass-> Write();
  hfJet_pt-> Write();
  hfJet_eta-> Write();
 
  fjet_deltaEta-> Write();
  fjet_deltaPhi -> Write();
  fjet_EhadOverEem-> Write();
  
   missingET_MET->Write();
   missingET_phi->Write();
   
   scalarHt->Write();
   
   
   outputfile->Close();  
   
}