#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_h1_e3e3bar/Events/run_01/h1_e3e3bar.root","READ");
     // TFile *f = new TFile("/home/ahmed/Downloads/MG5_aMC_v3_4_2/pp_h_tata/Events/run_01/h_tata.root","READ");
    //TFile *f = new TFile("/home/ahmed/Downloads/MG5_aMC_v3_4_2/pp_z_tata/Events/run_01/z_tata.root","READ");
    // TFile *f = new TFile("/home/ahmed/Downloads/MG5_aMC_v3_4_2/pp_tata/Events/run_01/tata.root","READ");
                                                      //***************
    
                                                                                                            
           
       
       // **************************************                           
                                   
   
   // weight calculation 
     
     
      
   //float cross_section =  11.87; //(in fb)- pp_h1_e3e3bar  (100 TeV)
  //float cross_section =  0.09916; //(in fb)- pp_h_tata(100 TeV)
    //  float cross_section =  5482 ; //(in fb)- pp_z_tata  (100 TeV)
     float cross_section =  8046 ; //(in fb)- pp_tata (100 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 *tau1_pt = new TH1D("tau1_pt", " tau1_pt distribution ",1000 ,0. , 1000.);
          tau1_pt ->GetXaxis()->SetTitle("pt_{tau1 }[Gev/c] ");
          tau1_pt->GetYaxis()->SetTitle("Events");
          
          TH1D *tau1_eta = new TH1D("tau1_eta", "eta1_Eta distribution ",100 ,-5.,5. );
          tau1_eta ->GetXaxis()->SetTitle("eta_{tau1 }");
          tau1_eta->GetYaxis()->SetTitle("Events");
                                   
          TH1D *tau2_pt = new TH1D("tau2_pt", " tau2_pt distribution ",1000 ,0. , 1000.);
          tau2_pt ->GetXaxis()->SetTitle("pt_{tau2 [Gev/c]}");
          tau2_pt->GetYaxis()->SetTitle("Events");
          
           TH1D *tau2_eta = new TH1D("tau2_eta", "eta2_Eta distribution ",100 ,-5.,5. );
          tau2_eta ->GetXaxis()->SetTitle("eta_{tau2 }");
          tau2_eta->GetYaxis()->SetTitle("Events");
       
       
         TH1D *h_mass = new TH1D("h_mass", " invariant mass distribution ",1000 ,0. , 1000.);
          h_mass ->GetXaxis()->SetTitle("M_{tau1 & tau2} [Gev]");
          h_mass->GetYaxis()->SetTitle("Events");
         
         
          TH1D *h_pt = new TH1D("h_pt", " h_pt distribution ",1000 ,0. , 1000.);
          h_pt ->GetXaxis()->SetTitle("pt_{tau1 & tau2 }[Gev/c] ");
          h_pt->GetYaxis()->SetTitle("Events");
                                       
          TH1D *h_eta = new TH1D("h_eta", "eta distribution ",100 ,-5.,5.  );
          h_eta ->GetXaxis()->SetTitle("eta_{tau1 & tau2}");
          h_eta->GetYaxis()->SetTitle("Events");
       
        TH1D *DR_tau = new TH1D("DR_tau", " DR_tau distribution ",1000 ,0 , 10.);
         DR_tau ->GetXaxis()->SetTitle("DR_tau_{tau1 & tau2} ");
         DR_tau->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> tauIndex;
                 
           	for (int i=0;i<Particle_size;i++){
 	      
 	        if(fabs(Particle_PID[i])!=15)continue;
 	        tauIndex.push_back(i);
         		  
		  }  
     
		 
	           
        
         TLorentzVector  tau1, tau2, h;        
          
          
            tau1.SetPtEtaPhiE(Particle_PT[ tauIndex[0]],Particle_Eta[tauIndex[0]],Particle_Phi[tauIndex[0]],Particle_E[tauIndex[0]]);
			 tau2.SetPtEtaPhiE(Particle_PT[tauIndex[1]],Particle_Eta[tauIndex[1]],Particle_Phi[tauIndex[1]],Particle_E[tauIndex[1]]);  
			      
          float  dR_tau;                  
            h=tau1+tau2;
            dR_tau= TMath::Sqrt(TMath::Power((Particle_Eta[tauIndex[0]] - Particle_Eta[tauIndex[1]]), 2) + TMath::Power((Particle_Phi[tauIndex[0]] - Particle_Phi[tauIndex[1]]), 2));
      
           
			
              
    tau1_pt->Fill(Particle_PT[tauIndex[0]], wt);
    tau1_eta->Fill(Particle_Eta[tauIndex[0]], wt);
    tau2_pt->Fill(Particle_PT[tauIndex[1]], wt);  
    tau2_eta->Fill(Particle_Eta[tauIndex[1]], wt); 
    h_mass->Fill(h.Mt(), wt);
    h_pt->Fill(h.Pt(), wt);
    h_eta->Fill(h.Eta(), wt);
    DR_tau ->Fill( dR_tau,wt);


  
   
  
}//end of event loop.



//**************outputfile
    
   
    
    
                         //*********************************
   //  TFile *outputfile = new TFile("h1_e3e3bar.root", "RECREATE");
   // TFile *outputfile = new TFile("h_tata.root", "RECREATE");
   //TFile *outputfile = new TFile("z_tata.root", "RECREATE");
     TFile *outputfile = new TFile("tata.root", "RECREATE");
     
                          //**********************8   
   

 
 
  tau1_pt-> Write();
  tau1_eta->Write();
  tau2_pt-> Write();
  tau2_eta->Write();
  h_mass->Write();
  h_pt ->Write();
  h_eta->Write();
  DR_tau->Write();

  
  
   outputfile->Close();  
   
}