#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. , 100.);
          b1_pt ->GetXaxis()->SetTitle("pt_{b1 }[Gev/c] ");
          b1_pt->GetYaxis()->SetTitle("Events");
          
          TH1D *b1_eta = new TH1D("b1_eta", "eta1_Eta distribution ",1000 ,-5.,9. );
          b1_eta ->GetXaxis()->SetTitle("eta_{b1 }");
          b1_eta->GetYaxis()->SetTitle("Events");
                                   
          TH1D *b2_pt = new TH1D("b2_pt", " b2_pt distribution ",3000 ,0. , 100.);
          b2_pt ->GetXaxis()->SetTitle("pt_{b2 [Gev/c]}");
          b2_pt->GetYaxis()->SetTitle("Events");
          
           TH1D *b2_eta = new TH1D("b2_eta", "eta2_Eta distribution ",1000 ,-3.,9. );
          b2_eta ->GetXaxis()->SetTitle("eta_{b2 }");
          b2_eta->GetYaxis()->SetTitle("Events");
       
       
         TH1D *h_mass = new TH1D("h_mass", " invariant mass distribution ",3000 ,0. , 100.);
          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. , 100.);
          h_pt ->GetXaxis()->SetTitle("pt_{b1 & b2 }[Gev/c] ");
          h_pt->GetYaxis()->SetTitle("Events");
                                       
          TH1D *h_eta = new TH1D("h_eta", "eta distribution ",1000 ,-5.,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. , 8.);
         DR_b ->GetXaxis()->SetTitle("DR_b_{b1 & b2} ");
         DR_b->GetYaxis()->SetTitle("Events");
       
          TH1D *scalarHt_Ht = new TH1D("scalarHt_Ht","scalarHt_Ht distribution",3000 ,0. , 100.);
          scalarHt_Ht->GetXaxis()->SetTitle("scalarHt_Ht[Gev/c]");
          scalarHt_Ht->GetYaxis()->SetTitle("Events");
                 
         int cutflow[10];          
          
                     
                          
                          
         
     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;
         
         
         
         
      
			  
           	for (int i=0;i<Particle_size;i++){           	        	
 	           	        	        
 	       
 	        if(fabs(Particle_PID[i])!=5)continue;
 	        b_quarkIndex.push_back(i); 
 	               	     
							 }         
 	           	           
            	  
       
         TLorentzVector  b1, b2, h;
        
          
                      
			  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]]);
        
           
                  
			 
             float  dR_b;
              h= b1+b2;     
        
          
       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));
      
          cutflow[0]+=cutflow[0];    
         if(Particle_PT[b_quarkIndex[0]] >40.)continue;
         cutflow[1]+=cutflow[1];  
         if(Particle_PT[b_quarkIndex[1]] >30.)continue;
         cutflow[2]+=cutflow[2];
         
          if(h.M() >40.) continue;
          cutflow[3]+=cutflow[3];      
          if(h.Pt() >40.) continue;
          cutflow[4]+=cutflow[4];      
          if(h.Eta() <10) continue;
         cutflow[5]+=cutflow[5];      
          if(ScalarHT_HT[0] >60)continue;
          cutflow[6]+=cutflow[6];
   
   
              
    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.M(), wt);
    h_pt->Fill(h.Pt(), wt);
    h_eta->Fill(h.Eta(), wt);
    DR_b ->Fill( dR_b,wt);
    scalarHt_Ht->Fill(ScalarHT_HT[0],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();
 scalarHt_Ht->Write();
  

  
   outputfile->Close(); 
 
  cout<<" Initial number of events = "<<cutflow[0]<<endl;
  cout<<" After b1 pT cut = "<<cutflow[1]<<endl;
  cout<<" After b2 pT cut = "<<cutflow[2]<<endl;
  cout<<" After higgs pT cut = "<<cutflow[3]<<endl;
  cout<<" After invariant Mass cut = "<<cutflow[4]<<endl;
  cout<<" After Higgs Eta cut = "<<cutflow[5]<<endl;
  cout<<" After scalarHt_ht cut = "<<cutflow[6]<<endl;

}