#define MetSign_cxx

#include "MetSign.h"
#include <TH2.h>
#include <TStyle.h>
#include <iostream>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <TCanvas.h>
#include <TLorentzVector.h>
#include <TDirectoryFile.h>
#include "TFile.h"
#include "TTree.h"

using namespace std;
/**
    CS-11 Asn 2
    checks.cpp
    Purpose: Calculates the Met Significance

    @author Dilia Portillo
    @version 2. 28-11-2016
*/


/* 
***Direct
root -l /data/atlas0/dportill/DougNtuples/mc15_13TeV.363365.Sherpa_NNPDF30NNLO_Zmumu_Pt0_70_CFilterBVeto.merge.DAOD_SUSY5.e4716_s2726_r7725_r7676_p2666.root
root -l /data/atlas0/dportill/DougNtuples/mc15_13TeV.361068.Sherpa_CT10_llvv.merge.DAOD_SUSY5.e3836_s2608_s2183_r7725_r7676_p2666_TRUTH.root
CollectionTree->Process("MetSign.C+")


***Proof:
root -l

.L Loader.C+
TProof::Open("");
  TChain* mychain=new TChain("CollectionTree")
  mychain->Add("/data/atlas0/dportill/DougNtuples/mc15_13TeV.363364.Sherpa_NNPDF30NNLO_Zmumu_Pt0_70_CVetoBVeto.merge.DAOD_SUSY5.e4716_s2726_r7725_r7676_p2666.root")
  mychain->Add("/data/atlas0/dportill/DougNtuples/mc15_13TeV.363366.Sherpa_NNPDF30NNLO_Zmumu_Pt0_70_BFilter.merge.DAOD_SUSY5.e4716_s2726_r7725_r7676_p2666.root")
  mychain->Add("/data/atlas0/dportill/DougNtuples/mc15_13TeV.363365.Sherpa_NNPDF30NNLO_Zmumu_Pt0_70_CFilterBVeto.merge.DAOD_SUSY5.e4716_s2726_r7725_r7676_p2666.root")
  mychain->SetProof();
  mychain->Process("MetSign.C")
*/
Int_t fNumberOfEvents; // Total number of events

void MetSign::Begin(TTree * /*tree*/)
{
  fNumberOfEvents = 0 ;
  cout<<"BEGIN"<<endl;

   TString option = GetOption();

   h_MetSig = 0 ;
   h_MSgNum = 0 ;
   h_MSgDen = 0 ;
   h_MSgRho = 0 ;
   h_VarPar = 0 ;
   h_VarPer = 0 ;
   h_Covari = 0 ;
   h_SqrMSg = 0 ;
   h_Met_sm = 0 ;
   h_Met_Ht = 0 ;
   h_Met_si = 0 ;
   h_LepPt_ = 0 ;
   h_LepEta = 0 ;
   h_LepPhi = 0 ;
   h_Lep_Nu = 0 ;
   h_Jet_Pt = 0 ;
   h_JetEta = 0 ;
   h_JetPhi = 0 ;
   h_Jet_Nu = 0 ;
   h_JetRes = 0 ;
   h_MetPhi = 0 ;
   h____Met = 0 ;
   h_softTm = 0 ;
   h_metjet = 0 ;

   h_MCchnu = 0 ;
   h_MCweig = 0 ;
   h_sqrtMS = 0 ;
   t_tree  = 0;
   vesseltree = 0;
   MetSig_t = 0;
   h_dilePt = 0 ;

   h_phi_reso_pt20  = 0 ; 
   h_phi_reso_pt50  = 0 ;
   h_phi_reso_pt100 = 0 ;


}

void MetSign::SlaveBegin(TTree * /*tree*/)
{
  TString option = GetOption();

  TFile *m_file = new TFile("/data/atlas0/dportill/PerformanceROC/data/jet_unc.root");

  if(m_file)
    {
      h_phi_reso_pt20 = static_cast<TH2F *>(m_file->Get("phi_reso_pt20"));
      if(!h_phi_reso_pt20) std::cerr << "Jet Uncertainty Histogram not valid" << std::endl;
      h_phi_reso_pt50 = static_cast<TH2F *>(m_file->Get("phi_reso_pt50"));
      if(!h_phi_reso_pt50) std::cerr << "Jet Uncertainty Histogram not valid" << std::endl;
      h_phi_reso_pt100 = static_cast<TH2F *>(m_file->Get("phi_reso_pt100"));
      if(!h_phi_reso_pt100) std::cerr << "Jet Uncertainty Histogram not valid" << std::endl;
    }



 //***************** Histogram definitions ***************************************8
 
  h_MetSig = new TH1D("h_Zmm_MetSig","",  900, 0.0 , 500.0   );
  h_MSgNum = new TH1D("h_Zmm_MSgNum","",  900, 0.0 , 500.0  );//700, 0.0 , 5000.0  );
  h_MSgDen = new TH1D("h_Zmm_MSgDen","",  900, 0.0 , 300.0   );
  h_MSgRho = new TH1D("h_Zmm_MSgRho","",  900, -2.0 , 2.0    );
  h_VarPar = new TH1D("h_Zmm_VarPar","",  900, 0.0 , 300.0   );
  h_VarPer = new TH1D("h_Zmm_VarPer","",  900, 0.0 , 300.0   );
  h_Covari = new TH1D("h_Zmm_Covari","",  900, -100.0 , 100.0);
  h_SqrMSg = new TH1D("h_Zmm_SqrMSg","",  900, 0.0 , 200.0   );
  h_Met_sm = new TH1D("h_Zmm_Met_sm","",  900, 0.0 , 100.0   );
  h_Met_Ht = new TH1D("h_Zmm_Met_Ht","",  900, 0.0 , 100.0   );
  h_Met_si = new TH1D("h_Zmm_Met_si","",  900, 0.0 , 100.0   );
  h_LepPt_ = new TH1D("h_Zmm_LepPt_","",  900, 0.0 , 300.0   );
  h_LepEta = new TH1D("h_Zmm_LepEta","",  900, -3.5 , 3.5    );
  h_LepPhi = new TH1D("h_Zmm_LepPhi","",  900, -3.5 , 3.5    );
  h_Jet_Pt = new TH1D("h_Zmm_Jet_Pt","",  900, 0.0 , 400.0   );
  h_JetEta = new TH1D("h_Zmm_JetEta","",  900, -5. , 5.      );
  h_JetPhi = new TH1D("h_Zmm_JetPhi","",  900, -4. , 4.      );
  h_JetRes = new TH1D("h_Zmm_JetRes","",  900, 0.0 , 0.3     );
  h_MetPhi = new TH1D("h_Zmm_MetPhi","",  900, -4. , 4.      );
  h_softTm = new TH1D("h_Zmm_softTm","",  900, 0., 400.      );
  h_metjet = new TH1D("h_Zmm_metjet","",  900, 0., 400.      );
  h_sqrtMS = new TH1D("h_Zmm_sqrtMS","",  900, 0.,100.); 
  h_dilePt = new TH1D("h_Zmm_dilePt","",  600, 0.0 , 2000.0   );
  h_Lep_Nu = new TH1D("h_Zmm_Lep_Nu","",  20, 0.0 , 10.0     );
  h_Jet_Nu = new TH1D("h_Zmm_Jet_Nu","",  20, 0.0 , 20.0     );
  h_MCchnu = new TH1D("h_Zmm_MCchnu","",  24, 363364.,363387.); //MC channel #
  h_MCweig = new TH1D("h_Zmm_MCweig","",  24, 363364.,363387.); //MC weight
  h____Met = new TH1D("h_Zmm____Met","",  1000, 0., 2000.      );//(600,0.,400.)
  
  fOutput->Add(h_MetSig);
  fOutput->Add(h_MSgNum);
  fOutput->Add(h_MSgDen);
  fOutput->Add(h_MSgRho);
  fOutput->Add(h_VarPar);
  fOutput->Add(h_VarPer);
  fOutput->Add(h_Covari);
  fOutput->Add(h_SqrMSg);
  fOutput->Add(h_Met_sm);
  fOutput->Add(h_Met_Ht);
  fOutput->Add(h_Met_si);
  fOutput->Add(h_LepPt_);
  fOutput->Add(h_LepEta);
  fOutput->Add(h_LepPhi);
  fOutput->Add(h_Lep_Nu);
  fOutput->Add(h_Jet_Pt);
  fOutput->Add(h_JetEta);
  fOutput->Add(h_JetPhi);
  fOutput->Add(h_Jet_Nu);
  fOutput->Add(h_JetRes);
  fOutput->Add(h_MetPhi);
  fOutput->Add(h____Met);
  fOutput->Add(h_softTm);
  fOutput->Add(h_metjet);
  fOutput->Add(h_MCchnu);
  fOutput->Add(h_MCweig);
  fOutput->Add(h_sqrtMS);
  fOutput->Add(h_dilePt);

  fOutput->Add(h_phi_reso_pt20);
  fOutput->Add(h_phi_reso_pt50);
  fOutput->Add(h_phi_reso_pt100);

  //***************** Tree definitions ***************************************
  //f1
  TString treefile = "/data/atlas0/dportill/DougNtuples/mc15_13TeV.363379.Sherpa_NNPDF30NNLO_Zmumu_Pt700_1000_CVetoBVeto.merge.DAOD_SUSY5.e4716_s2726_r7725_r7676_p2666.root";
  TFile *ftree = new TFile(treefile);
  vesseltree = (TTree*)ftree->Get("CollectionTree");
  //vesseltree->Print();

  //create the new file
  TFile *g = new TFile("./NTuples/test.root", "RECREATE");

  //TTree *newtree = vesseltree->CloneTree(0);

  t_tree = new TTree("testtree","Ntuples for Met Significance study");


  t_tree = vesseltree ->CloneTree(0);
  
  t_tree->Branch("Significance",&MetSig_t);
  t_tree->Branch("sqrtSignificance" , &sqrtMS_t );
  t_tree->Branch("NumeradorSignific" , &MSgNum_t );
  t_tree->Branch("DenominadorSignif" , &MSgDen_t );
  t_tree->Branch("Corr_fact_rho" , &MSgRho_t );
  t_tree->Branch("VarianceParalel" , &VarPar_t );
  t_tree->Branch("VariancePerpend" , &VarPer_t );
  t_tree->Branch("Covariance" , &Covari_t );
  t_tree->Branch("Met_sumet" , &Met_sm_t );
  t_tree->Branch("Met_Ht" , &Met_Ht_t );
  t_tree->Branch("Met_sigma" , &Met_si_t );
  t_tree->Branch("Met" , &Met____t );
  t_tree->Branch("MetPhi" , &MetPhi_t );
  t_tree->Branch("LepEta" ,   &LepEta_t );
  t_tree->Branch("LepPhi" ,   &LepPhi_t );
  t_tree->Branch("LepPt_" ,   &LepPt__t );
  t_tree->Branch("Lep_Number",&Lep_Nu_t );
  t_tree->Branch("softTermm" , &softTm_t );
  t_tree->Branch("ZPt" , &dilePt_t );
  t_tree->Branch("Jet_Pt" , &Jet_Pt_t );
  t_tree->Branch("JetEta" , &JetEta_t );
  t_tree->Branch("JetPhi" , &JetPhi_t );
  t_tree->Branch("Jet_Number" , &Jet_Nu_t );
  t_tree->Branch("JetResolution" , &JetRes_t );
  t_tree->Branch("JetMetTerm" , &JetMet_t );
  t_tree->Branch("JetPileUp" , &Jetpup_t );
  t_tree->Branch("Jet_PassOverlapRemoval" ,&Jet_OR_t );
  t_tree->Branch("JetJvt_t", &JetJvt_t );
  t_tree->Branch("JetFracUsedInMet" ,&JetFMt_t );
  t_tree->Branch("Muo_Pt" ,&Muo_Pt_t );
  t_tree->Branch("MuoEta" ,&MuoEta_t );
  t_tree->Branch("MuoPhi" ,&MuoPhi_t );
  t_tree->Branch("Muo_Number" ,&Muo_Nu_t );
  t_tree->Branch("MuoResolution" ,&MuoRes_t );
  t_tree->Branch("MuoType" , &MuoTyp_t );
  t_tree->Branch("MuoCharge" ,&Muocha_t );
  t_tree->Branch("MuoisoLooseTrackOnly" , &MuoITO_t );
  t_tree->Branch("Muo_PassOverlapRemoval",&Muo_OR_t );
  t_tree->Branch("Ele_Pt" ,   &Ele_Pt_t );
  t_tree->Branch("EleEta" ,   &EleEta_t );
  t_tree->Branch("ElePhi" ,   &ElePhi_t );
  t_tree->Branch("Ele_Number",&Ele_Nu_t );
  t_tree->Branch("EleResolution", &EleRes_t );
  t_tree->Branch("Ele_PassOverlapRemoval" ,  &Ele_OR_t);
  t_tree->Branch("weight" ,             &weight_t );              
  t_tree->Branch("MCEventweight" ,      &MCweig_t );
  t_tree->Branch("PileUp_weight" ,      &pup_we_t );
  t_tree->Branch("weightSherpa22_njets",&wshenj_t );
  t_tree->Branch("sf_jvt" ,&sf_jvt_t );
  t_tree->Branch("sf_el" , &sf_el__t );
  t_tree->Branch("sf_mu" , &sf_mu__t );
  t_tree->Branch("sf_muo_trigger" , &sfmtri_t );
  t_tree->Branch("sf_ele_trigger" , &sfetri_t );
  t_tree->Branch("TotatWeight" , &TotWGT_t );
  t_tree->Branch("Z_mass" ,      &Z_mass_t );
  t_tree->Branch("EleVarParal" , &EleVPa_t );
  t_tree->Branch("EleVarPerpe" , &EleVPe_t );
  t_tree->Branch("MuoVarParal" , &MuoVPa_t );
  t_tree->Branch("MuoVarPerpe" , &MuoVPe_t );
  t_tree->Branch("JetVarParal" , &JetVPa_t );
  t_tree->Branch("JetVarPerpe" , &JetVPe_t );
  t_tree->Branch("JetPhiReso" ,  &JetPre_t );



  
  fOutput->Add(t_tree);
}

Bool_t MetSign::Process(Long64_t entry)
{
  //Initializacion cov. matrix components
  Double_t Var_L = 0 ;
  Double_t Var_T = 0 ;
  Double_t Cv_LT = 0 ;

  //Selection number of jets:
  int Nj = 0; 

  // Counter of jets
  int NJets = 0; 
  // Counter of muons
  int NMuon = 0; 
  // Counter of Electrons
  int NElec = 0; 
  
  // 1 if the event has a jet that pass the selection and 0 to veto the event
  int Veto = 0 ;

  //****************** Branches *******************************

  vesseltree ->GetEntry(fNumberOfEvents);
  
  b_met       ->GetEntry(entry);
  b_met_jet   ->GetEntry(entry); 
  b_met_soft  ->GetEntry(entry);

  b_met_phi   ->GetEntry(entry);
  b_n_mu      ->GetEntry(entry);
  b_n_dilep   ->GetEntry(entry);  
  b_met_sumet ->GetEntry(entry);
  b_ht        ->GetEntry(entry);
  b_ht_sig    ->GetEntry(entry);
  b_mu_trigger->GetEntry(entry);

  b_dilep_m   ->GetEntry(entry);   //Vector [n_dilep]
  b_dilep_pt  ->GetEntry(entry);

  b_n_jet     ->GetEntry(entry);
  b_jet_pTreso->GetEntry(entry);   //Vector [n_jet] 
  b_jet_phi   ->GetEntry(entry);   //Vector [n_jet]
  b_jet_pt    ->GetEntry(entry);
  b_jet_eta   ->GetEntry(entry);

  b_n_lep    ->GetEntry(entry);
  b_lep_pt   ->GetEntry(entry);
  b_lep_eta  ->GetEntry(entry);
  b_lep_phi  ->GetEntry(entry);

  b_n_basejet     ->GetEntry(entry);
  b_basejet_pTreso->GetEntry(entry);//[n_basejet]
  b_basejet_phi   ->GetEntry(entry);//[n_basejet]
  b_basejet_pt    ->GetEntry(entry);//[n_basejet]
  b_basejet_eta   ->GetEntry(entry);
  b_basejet_pileup->GetEntry(entry);
  b_basejet_passOR->GetEntry(entry);
  b_basejet_Jvt   ->GetEntry(entry);  //0.64--> 0.59
  b_basejet_fracUsedInMET->GetEntry(entry);

  b_n_basemu       ->GetEntry(entry);
  b_basemu_pt	   ->GetEntry(entry);
  b_basemu_passOR ->GetEntry(entry);
  //b_basemu_muType ->GetEntry(entry);  //Not needed for selection
  b_basemu_pTreso ->GetEntry(entry);
  b_basemu_phi   ->GetEntry(entry);//[n_basejet]
  b_basemu_pt    ->GetEntry(entry);//[n_basejet]
  b_basemu_eta   ->GetEntry(entry);
  b_basemu_charge->GetEntry(entry);

  b_n_mu       ->GetEntry(entry);
  b_mu_muType  ->GetEntry(entry);
  // b_mu_passOR ->GetEntry(entry); //not defined for muo
  b_mu_pTreso ->GetEntry(entry);
  b_mu_phi   ->GetEntry(entry);
  b_mu_pt    ->GetEntry(entry);
  b_mu_eta   ->GetEntry(entry);
  b_mu_charge->GetEntry(entry);
  b_mu_isoLooseTrackOnly->GetEntry(entry);

  b_n_baseel       ->GetEntry(entry);
  b_baseel_pt	   ->GetEntry(entry);
  b_baseel_passOR ->GetEntry(entry);
  b_baseel_Enreso ->GetEntry(entry);
  b_baseel_phi   ->GetEntry(entry);//[n_basejet]
  b_baseel_pt    ->GetEntry(entry);//[n_basejet]
  b_baseel_eta   ->GetEntry(entry);

  //b_truth_Z_pt->GetEntry(entry);  //No filled

  b_weight               ->GetEntry(entry);
  b_mc_event_weight      ->GetEntry(entry);
  b_pileup_weight        ->GetEntry(entry);
  b_weight_sherpa22_njets->GetEntry(entry);
  b_sf_jvt                 ->GetEntry(entry);
  b_sf_el      	           ->GetEntry(entry);
  b_sf_mu    	           ->GetEntry(entry);
  b_sf_mu_trigger            ->GetEntry(entry);
  b_sf_el_trigger            ->GetEntry(entry);
  b_mc_channel_number ->GetEntry(entry);

  //Existing Tree
  b_n_trackjet ->GetEntry(entry);
  b_trackjet_pt->GetEntry(entry);//?
  b_n_bph          ->GetEntry(entry);
  b_n_bjet	       ->GetEntry(entry);
  b_n_lep	       ->GetEntry(entry);
  b_n_baseel       ->GetEntry(entry);
  b_n_el	       ->GetEntry(entry);
  b_n_basemu       ->GetEntry(entry);
  b_n_mu	       ->GetEntry(entry);
  b_n_ph	       ->GetEntry(entry);
  b_n_truth_photon ->GetEntry(entry);
  b_n_truth_top    ->GetEntry(entry);
  b_n_truth_el     ->GetEntry(entry);
  b_n_truth_jet    ->GetEntry(entry);
  b_n_truth_mu     ->GetEntry(entry);
  b_n_truth_ph     ->GetEntry(entry);
  b_n_truth_tau    ->GetEntry(entry);
  b_n_truth_W      ->GetEntry(entry);
  b_n_truth_H      ->GetEntry(entry);
  b_n_truth_Z      ->GetEntry(entry);
  b_n_truth_B      ->GetEntry(entry);
  b_n_truth_C      ->GetEntry(entry);
  b_n_dilep        ->GetEntry(entry);
  b_met	       ->GetEntry(entry);
  b_met_x	       ->GetEntry(entry);
  b_met_y	       ->GetEntry(entry);
  b_met_phi        ->GetEntry(entry);
  b_met_sumet      ->GetEntry(entry);
  b_met_noEL       ->GetEntry(entry);
  b_met_noEL_x     ->GetEntry(entry);
  b_met_noEL_y     ->GetEntry(entry);
  b_met_noEL_phi   ->GetEntry(entry);
  b_met_noEL_sumet ->GetEntry(entry);
  b_met_noMU       ->GetEntry(entry);
  b_met_noMU_x     ->GetEntry(entry);
  b_met_noMU_y     ->GetEntry(entry);
  b_met_noMU_phi   ->GetEntry(entry);
  b_met_noMU_sumet ->GetEntry(entry);
  b_met_noJVT      ->GetEntry(entry);
  b_met_noJVT_x    ->GetEntry(entry);
  b_met_noJVT_y    ->GetEntry(entry);
  b_met_noJVT_phi  ->GetEntry(entry);
  b_met_noJVT_sumet->GetEntry(entry);
  b_met_invis      ->GetEntry(entry);
  b_met_invis_x    ->GetEntry(entry);
  b_met_invis_y    ->GetEntry(entry);
  b_met_invis_phi  ->GetEntry(entry);
  b_met_invis_sumet->GetEntry(entry);
  b_met_el         ->GetEntry(entry);
  b_met_el_x       ->GetEntry(entry);    
  b_met_el_y       ->GetEntry(entry);    
  b_met_el_phi     ->GetEntry(entry);
  b_met_el_sumet   ->GetEntry(entry);
  b_met_mu	       ->GetEntry(entry);
  b_met_mu_x       ->GetEntry(entry);    
  b_met_mu_y       ->GetEntry(entry);    
  b_met_mu_phi     ->GetEntry(entry);
  b_met_mu_sumet   ->GetEntry(entry);
  b_met_jet	       ->GetEntry(entry);
  b_met_jet_x      ->GetEntry(entry);    
  b_met_jet_y      ->GetEntry(entry);    
  b_met_jet_phi    ->GetEntry(entry);
  b_met_jet_sumet  ->GetEntry(entry);
  b_met_soft       ->GetEntry(entry);    
  b_met_soft_x     ->GetEntry(entry);
  b_met_soft_y     ->GetEntry(entry);
  b_met_soft_phi   ->GetEntry(entry);
  b_met_soft_sumet ->GetEntry(entry);
  b_met_cst        ->GetEntry(entry);
  b_met_cst_x      ->GetEntry(entry);    
  b_met_cst_y      ->GetEntry(entry);    
  b_met_cst_phi    ->GetEntry(entry);
  b_met_cst_sumet  ->GetEntry(entry);
  b_met_cst_soft   ->GetEntry(entry);
  b_met_photons    ->GetEntry(entry);
  b_met_photons_x  ->GetEntry(entry);
  b_met_photons_y  ->GetEntry(entry);
  b_met_photons_phi->GetEntry(entry);
  b_met_photons_sumet->GetEntry(entry);
  










  
  // ********* Weights *************************
 
  Double_t Lumi_0 = 1. ; //?
  Double_t NumEvnets = 1. ; //? MetSign::GetXSxFiltEff(mc_channel_number) includes the number of events for the Zmumu

  Double_t XSecWeight = (MetSign::GetXSxFiltEff(mc_channel_number) *  Lumi_0)   ;

  Double_t WEIGHT = XSecWeight* weight *weight_sherpa22_njets * sf_jvt * sf_el * sf_mu * sf_mu_trigger * sf_el_trigger  ;

  Int_t Bin = h_MCweig->GetXaxis()->FindBin(mc_channel_number);
  h_MCweig->SetBinContent(Bin, mc_event_weight);

 
  if(mu_trigger!=1) return kFALSE;
  if(n_mu<2) return kFALSE;//if(n_basemu<2) return kFALSE;

  TLorentzVector mu1;
  TLorentzVector mu2;

  Int_t Countm1=0;
  Int_t Countm2=-1;
  
  Double_t Zcand_mass=0;
  Int_t chargemuons;
  
  // ***Z mass**

  for (int i=0; i<n_mu; i++)
    {
      if(mu_pt[i] > 25.0e3 && mu_muType[i] < 2.5 && mu_isoLooseTrackOnly[i]>0.5 && Countm2!= i)//&& mu_muType[i] >2.5
	{
	  mu1.SetPtEtaPhiM(mu_pt[i]*0.001,mu_eta[i], mu_phi[i],0.1);
	  Countm1 = i;
	  for (int j= i+1; j<n_mu; j++)
	    {
	      if(mu_pt[j] > 25.0e3 && mu_muType[j] < 2.5  && mu_isoLooseTrackOnly[j]>0.5)//&& mu_muType[j] >2.5
		{
		  mu2.SetPtEtaPhiM(mu_pt[j]*0.001,mu_eta[j], mu_phi[j],0.1);
		  Countm2 = j;
		  chargemuons = mu_charge[i]*mu_charge[j];
		  if(chargemuons<0)
		    {
		      Zcand_mass = ( mu1 + mu2 ).M();
		      if(Zcand_mass > 66 && Zcand_mass < 116)
			{
			  h_dilePt->Fill(( mu1 + mu2 ).Pt(), WEIGHT );
			  Z_mass_t = ( mu1 + mu2 ).Pt();
			  break;
			}
		    }//Charge Selection
		}//Muon 2 Selection
	    } //Loop Muon 2
	}//Muon 1 Selection (baseline muons)
    } //Loop on Muons 

      //***************** Selection ******************************
      //***************** Selection ******************************

  cout << "Zcand_mass: "<<Zcand_mass<<endl;
  
  if(Zcand_mass > 66. && Zcand_mass < 116.)//&& (met*0.001)>300.)//fabs(dilep_m[0]*0.001-91.)<=25.) Baseline Event selection
    {
      cout<< "***************Event #"<<fNumberOfEvents<< " Pass Selection*****************"<<endl;
      Double_t Phi = met_phi; //Angle for rotation of the cov matrix
      Double_t Met = met*0.001 ; // Numerator
      // cout<<"Event with Jets! ---> n_basejet != 0"<<endl;
      cout<< "- WEIGHT = mc_event_weight*XSecWeight: "<<WEIGHT<<endl;
      cout<<"- Met [GeV]: "<<met*0.001<<endl;
      cout<<"- Met Phi: "<<met_phi<<endl;
      
      //************************************* //
      //*************** JETS  *************** //
      //Initializacion cov. matrix components for the container = JET
      Double_t Var_L_j = 0 ;
      Double_t Var_T_j = 0 ;
      Double_t Cv_LT_j = 0 ;
      tuple<double,double,double> CovMatrix_j; //Cov matrix entries jets
      //Loop over number of Jets
      cout<<"___________________________________"<<endl;
      cout<<"          JETS         "<<endl;
      cout<<"___________________________________"<<endl;
      cout<<"- Number of Jets (n_basejet): "<<n_basejet<<endl;
      cout<<"*Loop over Jets:*"<<endl;
      cout<<"Selecting Jets with basejet_fracUsedInMET[i]>0"<<endl;
      // cout<<"___________________________________"<<endl;
      for (int i=0; i<n_basejet; i++)
	{
	  if( basejet_fracUsedInMET[i]>0  )
	    {
	      cout<< "Jet #"<<i<<" in MET"<<endl;
	      NJets = NJets + 1;
	    }//Jet Selection (baseline jets)
	}
     cout<<"- Number of Jets that pass selection: "<<NJets<<endl;
      for ( int i=0; i<n_basejet; i++)
	{
	  if( NJets>=Nj  && basejet_fracUsedInMET[i]>0.0 )//!basejet_pileup[i] && basejet_passOR[i]>0.5)
	    { 
	      Double_t VarPhi = MetSign::GetPhiUnc(basejet_eta[i], basejet_phi[i], basejet_pt[i]*basejet_fracUsedInMET[i]*0.001);//PHI RESOLUTION Calculate the delta phi(reco-true)->angular variance 
	      Double_t Pt_unc = basejet_pTreso[i]*basejet_fracUsedInMET[i]*basejet_pt[i]*0.001; //sqrt of Entry 1,1 of Cov matrix in Jet direction
	      Double_t Phi_unc = VarPhi * basejet_fracUsedInMET[i]*basejet_pt[i]*0.001; //sqrt of Entry 2,2 of Cov matrix in Jet direction
	      Double_t var_L, var_T,cv_LT; //Cov for each jet in parallel-transverse
	      ////xy Resolutions
	      Double_t PxReso = Pt_unc *cos(basejet_phi[i]);
	      Double_t PyReso = Pt_unc *sin(basejet_phi[i]);
	      cout<<"Jet #"<<i<<" pass the selection"<<endl;
	      cout<<" - basejet_phi : "<< basejet_phi[i] <<endl;
	      cout<<" - basejet_pt [GeV] : "<<basejet_pt[i]*0.001 <<" - basejet_pt*fracUsedInMET [GeV] : " <<basejet_pt[i]*basejet_fracUsedInMET[i]*0.001<<endl;
	      cout<<" - basejet_pTreso : "<<basejet_pTreso[i]*basejet_pt[i]*0.001<<" - basejet_pTreso*fracUsedInMET : " <<basejet_pTreso[i]*basejet_fracUsedInMET[i] <<endl;
	      cout<<" - basejet_VarPhi : "<<VarPhi<<endl;
	      cout<<" - PxResoJet [GeV]: "<< PxReso <<endl;
	      cout<<" - PyResoJet [GeV]: "<< PyReso <<endl;
	      ////Covariance matrix in x-y coordinate system
	      Double_t Var2_x, Var2_y, Cv2_xy;
	      tie(Var2_x, Var2_y, Cv2_xy) = MetSign::CovMatrixRotation( pow(Pt_unc,2),pow( Phi_unc,2), 0,basejet_phi[i] );
	      cout<<"Pt_unc: "<<Pt_unc<<", Phi_unc: "<<Phi_unc<<endl;
	      Double_t Rhoxy = Cv2_xy / sqrt( Var2_x * Var2_y ); //For one jet Rho can be negative or positive in xy
	      ////Covariance matrix in x-y coordinate system without Phi Reso 
	      Double_t Var_x = pow(PxReso,2);
	      Double_t Var_y = pow(PyReso,2);
	      Double_t Cv_xy = PxReso * PyReso;  // 100 % correlated -> rho=1
	      cout<<" - Variance in x [GeV]^2: "<< Var_x <<" with Phi Reso: " << Var2_x <<endl;
	      cout<<" - Variance in y [GeV]^2: "<< Var_y <<" with Phi Reso: " << Var2_y <<endl;
	      cout<<" - Covariance xy [GeV]^2: "<< Cv_xy <<" with Phi Reso: " << Cv2_xy <<endl;
	      ////Rotation (components)
	      //tie(var_L, var_T,cv_LT) = MetSign::CovMatrixRotation( Var_x, Var_y, Cv_xy, Phi); //Whitout Phi Reso
	      tie(var_L, var_T,cv_LT) = MetSign::CovMatrixRotation( Var2_x, Var2_y, Cv2_xy, Phi);
	      ////Covariance matrix rotaded for all jets
	      Var_L_j = var_L + Var_L_j ;
	      Var_T_j = var_T + Var_T_j ;
	      Cv_LT_j = cv_LT + Cv_LT_j ;
	       cout<<" - Jet[i] Variance parallel  [GeV]^2: "<< var_L <<endl;
	       cout<<" - Jet[i] Variance perpendicular  [GeV]^2: "<< var_T <<endl;
	       cout<<" - Jet[i] Covariance   [GeV]^2: "<< cv_LT <<endl;
	      // *** HISTOGRAM TREE FILLING --> Jets***
	      h_Jet_Pt-> Fill( basejet_pt[i]*0.001  ,weight);	      
	      h_JetEta-> Fill( basejet_eta[i] ,weight);
	      h_JetPhi-> Fill( basejet_phi[i] ,weight);
	      h_JetRes-> Fill( basejet_pTreso[i] ,weight);
	      h_Jet_Nu -> Fill( NJets  ,WEIGHT);
	      Jet_Pt_t = basejet_pt[i]*0.001 ;
	      JetEta_t = basejet_eta[i] ;
	      JetPhi_t = basejet_phi[i] ;
	      JetRes_t = basejet_pTreso[i] ;
	      Jetpup_t = basejet_pileup[i]  ;
	      Jet_OR_t = basejet_passOR[i]  ;
	      JetJvt_t = basejet_Jvt[i]  ;
	      JetFMt_t = basejet_fracUsedInMET[i]  ;
	      JetPre_t = VarPhi;
	      Veto = 1; //Good Event
	      // cout<<"___________________________________"<<endl;
	    }// Jet Selection and Cov Selection for baseline jets
	  else
	    {
	      cout<<"Jet #"<<i<<" not pass the selection"<<endl;
	      //cout<<"___________________________________"<<endl;
	    }
	}//loop over Jets
      cout<<" - - - - - - - - - - - - - - - - - - - -  - - -  "<<endl;
      cout<<" - Jet Variance parallel  [GeV]^2: "<< Var_L_j <<endl;
      cout<<" - Jet Variance perpendicular  [GeV]^2: "<< Var_T_j <<endl;
      cout<<" - Jet Covariance   [GeV]^2: "<< Cv_LT_j <<endl;
      //// Sum the jet contribution to the total Cov Matrix:
      Var_L = Var_L_j + Var_L ;
      Var_T = Var_T_j + Var_T ;
      Cv_LT = Cv_LT_j + Cv_LT ;

      //*******************************************//
      //*************** SOFT TERMS  *************** //  
      if(NJets>= Nj )
	{
	  Double_t SoftTermReso = pow(13.4,2.);
	  Double_t var_L_st=0, var_T_st=0,cv_LT_st=0; //Cov for soft Term
	  tie(var_L_st, var_T_st,cv_LT_st) = MetSign::CovMatrixRotation(SoftTermReso ,SoftTermReso , 0., Phi);
	  cout<<"___________________________________"<<endl;
	  cout<<"          Soft Terms         "<<endl;
	  cout<<"___________________________________"<<endl;
	  cout<<" - ST Variance parallel  [GeV]^2: "<< var_L_st <<endl;
	  cout<<" - ST Variance perpendicular  [GeV]^2: "<< var_T_st <<endl;
	  cout<<" - ST Covariance   [GeV]^2: "<< cv_LT_st <<endl;
	  // Sum the ST contribution to the total Cov Matrix:
	  Var_L = var_L_st + Var_L ;
	  Var_T = var_T_st + Var_T ;
	  Cv_LT = cv_LT_st + Cv_LT ;
        }


      //*******************************************//
      //***************   MUONS     *************** // 
      ////Initializacion cov. matrix components for the container = Muon
      Double_t Var_L_m = 0 ;
      Double_t Var_T_m = 0 ;
      Double_t Cv_LT_m = 0 ;
      cout<<"___________________________________"<<endl;
      cout<<"            MUONS           "<<endl;
      cout<<"___________________________________"<<endl;
      cout<<"- Number of Muons (n_basemu): "<<n_basemu<<endl;
      cout<<"*Loop over Muons:*"<<endl;
      cout<<"Selecting Muons: mu_pt[i] > 7.0e3 && mu_passOR[i] > 0.5"<<endl;
      for (int i=0; i<n_mu; i++)
	{
	  if(mu_pt[i] > 7.0e3 && mu_muType[i] <2.5  )
	    {
	      cout<< "Muo #"<<i<<" is good"<<endl;
	      NMuon = NMuon +1;
	    }//Muon Selection (baseline muons)
	}
      Muo_Nu_t = NMuon;//n_mu  ;
      cout<<"- Number of muons that pass selection: "<<NMuon<<endl;  
      for ( int i=0; i<n_mu; i++)
	{
	  if( NJets>= Nj && mu_pt[i] > 7.0e3 && mu_muType[i] <2.5 )
	    {
	      Double_t VarPhi = 0.01;//PHI RESOLUTION 
	      Double_t Pt_unc = mu_pTreso[i]*mu_pt[i]*0.001; //sqrt of Entry 1,1 of Cov matrix in Ele direction
	      Double_t Phi_unc = VarPhi*mu_pt[i]*0.001;//sqrt of Entry 2,2 of Cov matrix in Ele direction
	      Double_t var_L, var_T,cv_LT; //Cov for each muon
	      ////xy Resolutions
	      Double_t PxReso = mu_pTreso[i]*cos(mu_phi[i])*mu_pt[i]*0.001;
	      Double_t PyReso = mu_pTreso[i]*sin(mu_phi[i])*mu_pt[i]*0.001;
	      cout<<"Muo #"<<i<<" pass the selection"<<endl;	     
	      cout<<" - mu_pTreso : "<<mu_pTreso[i]  <<endl;
	      cout<<" - mu_phi : "<< mu_phi[i] <<endl;
	      cout<<" - mu_pt [GeV] : "<<mu_pt[i]*0.001  <<endl;
	      cout<<" - PxResoMuo [GeV]: "<< PxReso <<endl;
	      cout<<" - PyResoMuo [GeV]: "<< PyReso <<endl;
	      ////Covariance matrix in x-y coordinate system 
	      Double_t Var2_x, Var2_y, Cv2_xy; //With PhiReso
	      tie(Var2_x, Var2_y, Cv2_xy) = MetSign::CovMatrixRotation( pow(Pt_unc,2),pow( Phi_unc,2), 0,mu_phi[i] );
	      ////Covariance matrix in x-y coordinate system without Phi Reso 
	      Double_t Var_x = pow(PxReso,2);
	      Double_t Var_y = pow(PyReso,2);
	      Double_t Cv_xy = PxReso * PyReso;  // 100 % correlated-> rho=1
	      cout<<" - Variance in x [GeV]^2: "<< Var_x <<" with Phi Reso: " << Var2_x <<endl;
	      cout<<" - Variance in y [GeV]^2: "<< Var_y <<" with Phi Reso: " << Var2_y <<endl;
	      cout<<" - Covariance xy [GeV]^2: "<< Cv_xy <<" with Phi Reso: " << Cv2_xy <<endl;
	      ////Rotation (components)
      	      tie(var_L, var_T,cv_LT) = MetSign::CovMatrixRotation( Var2_x, Var2_y, Cv2_xy, Phi);
	      //tie(var_L, var_T,cv_LT) = MetSign::CovMatrixRotation( Var_x, Var_y, Cv_xy, Phi);
	      Var_L_m = var_L + Var_L_m ;
	      Var_T_m = var_T + Var_T_m ;
	      Cv_LT_m = cv_LT + Cv_LT_m ;
	      Veto = 1; //Good Event
	      // *** HISTOGRAM FILLING *** //Histograms to define!
	      // h_MuoPt_-> Fill( mu_pt[i]*0.001  ,weight);	      
	      // h_MuoEta-> Fill( mu_eta[i] ,weight);
	      // h_MuoPhi-> Fill( mu_phi[i] ,weight);
	      // h_MuoRes-> Fill( mu_pTreso[i] ,weight);
	      // h_Muo_Nu -> Fill( NMuon  ,WEIGHT);
	      Muo_Pt_t = mu_pt[i]  ;
	      MuoEta_t = mu_eta[i]  ;
	      MuoPhi_t = mu_phi[i]  ;
	      MuoRes_t = mu_pTreso[i]  ;
	      MuoTyp_t = mu_muType[i]  ;
	      Muocha_t = mu_charge[i]  ;
	      MuoITO_t = mu_isoLooseTrackOnly[i] ;
	    }//Selection muons
	  else
	    {
	      cout<<"Muon #"<<i<<" not pass the selection"<<endl;
	    }
	}//loop over muons
      // Sum the Muon contribution to the total Cov Matrix:
      Var_L = Var_L_m + Var_L ;
      Var_T = Var_T_m + Var_T ;
      Cv_LT = Cv_LT_m + Cv_LT ;
      // *** Printing ***
      cout<<" - - - - - - - - - - - - - - - - - - - -  - - -  "<<endl;
      cout<<" - *Muon Variance parallel  [GeV]^2: "<< Var_L_m <<endl;
      cout<<" - *Muon Variance perpendicular  [GeV]^2: "<< Var_T_m <<endl;
      cout<<" - *Muon Covariance   [GeV]^2: "<< Cv_LT_m <<endl;

      //*******************************************//
      //*************** ELECTRONS  *************** // 
      ////Initializacion cov. matrix components for the container = Muon
      Double_t Var_L_e = 0 ;
      Double_t Var_T_e = 0 ;
      Double_t Cv_LT_e = 0 ;
      ////Loop over number of Electrons
      cout<<"___________________________________"<<endl;
      cout<<"           ELECTRONS           "<<endl;
      cout<<"___________________________________"<<endl;
      cout<<"- Number of Electrons (n_baseel): "<<n_baseel<<endl;
      cout<<"*Loop over Electrons:*"<<endl;
      cout<<"Selecting Electrons: baseel_pt[i] > 7.0e3 && baseel_passOR[i] > 0.5"<<endl;
      for (int i=0; i<n_baseel; i++)
      	{
      	  if(baseel_pt[i] > 7.0e3 && baseel_passOR[i] > 0.5 )
      	    {
      	      cout<< "Ele #"<<i<<" is good"<<endl;
      	      NElec = NElec +1;
      	    }//Ele Selection (baseline ele)
      	}
      Ele_Nu_t = NElec ;
      cout<<"- Number of electrons that pass selection: "<< NElec <<endl;
      for ( int i=0; i<n_baseel; i++)
      	{
      	  if( NJets>= Nj && baseel_pt[i] > 7.0e3 && baseel_passOR[i] > 0.5 )
      	    {
      	      Double_t theta = 2 * atan ( exp(baseel_eta[i] ) );
	      Double_t VarPhi = 0.02;//PHI RESOLUTION 
	      Double_t Pt_unc = baseel_Enreso[i]*sin(theta)*baseel_pt[i]*0.001; //sqrt of Entry 1,1 of Cov matrix in Ele direction
	      Double_t Phi_unc = VarPhi*baseel_pt[i]*0.001;//sqrt of Entry 2,2 of Cov matrix in Ele direction
      	      Double_t var_L, var_T,cv_LT; //Cov for each muon
      	      ////xy Resolutions
      	      Double_t PxReso = baseel_Enreso[i]*cos(baseel_phi[i])*baseel_pt[i]*sin(theta)*0.001;
      	      Double_t PyReso = baseel_Enreso[i]*sin(baseel_phi[i])*baseel_pt[i]*sin(theta)*0.001;
      	      cout<<"Ele #"<<i<<" pass the selection"<<endl;	     
      	      cout<<" - baseel_pTreso : "<<baseel_Enreso[i]  <<endl;
      	      cout<<" - baseel_phi : "<< baseel_phi[i] <<endl;
      	      cout<<" - baseel_pt [GeV] : "<<baseel_pt[i]*0.001  <<endl;
      	      cout<<" - PxResoEle [GeV]: "<< PxReso <<endl;
      	      cout<<" - PyResoEle [GeV]: "<< PyReso <<endl;
      	      ////Covariance matrix in x-y coordinate system 
	      Double_t Var2_x, Var2_y, Cv2_xy; //With PhiReso
	      tie(Var2_x, Var2_y, Cv2_xy) = MetSign::CovMatrixRotation( pow(Pt_unc,2),pow( Phi_unc,2), 0,baseel_phi[i] );
	      ////Covariance matrix in x-y coordinate system without Phi Reso 
      	      Double_t Var_x = pow(PxReso,2);  //No PhiReso
      	      Double_t Var_y = pow(PyReso,2);  //No PhiReso
      	      Double_t Cv_xy = PxReso * PyReso;//No PhiReso  // 100 % correlated -> rho=1
	      cout<<" - Variance in x [GeV]^2: "<< Var_x <<" with Phi Reso: " << Var2_x <<endl;
	      cout<<" - Variance in y [GeV]^2: "<< Var_y <<" with Phi Reso: " << Var2_y <<endl;
	      cout<<" - Covariance xy [GeV]^2: "<< Cv_xy <<" with Phi Reso: " << Cv2_xy <<endl;
      	      ////Rotation (components)
      	      tie(var_L, var_T,cv_LT) = MetSign::CovMatrixRotation( Var2_x, Var2_y, Cv2_xy, Phi);
      	      Var_L_e = var_L + Var_L_e ;
      	      Var_T_e = var_T + Var_T_e ;
      	      Cv_LT_e = cv_LT + Cv_LT_e ;
      	      Veto = 1; //Good Event
      	      // // *** HISTOGRAM-TREE FILLING ***
      	      // h_ElePt_-> Fill( baseel_pt[i]*0.001  ,weight);	      
      	      // h_EleEta-> Fill( baseel_eta[i]     ,weight);
      	      // h_ElePhi-> Fill( baseel_phi[i]     ,weight);
      	      // h_EleRes-> Fill( baseel_Enreso[i]  ,weight);
      	      // h_Ele_Nu -> Fill( NElec  ,WEIGHT);
	      Ele_Pt_t =  baseel_pt[i]*0.001 ;
	      EleEta_t =  baseel_eta[i]      ;
	      ElePhi_t =  baseel_phi[i]      ;
	      EleRes_t =  baseel_Enreso[i] ;
	      Ele_OR_t =  baseel_passOR[i] ;
      	    }//Selection electrons
      	}//loop electrons
      // Sum the Electron contribution to the total Cov Matrix:
      Var_L = Var_L_e + Var_L ;
      Var_T = Var_T_e + Var_T ;
      Cv_LT = Cv_LT_e + Cv_LT ;
      // *** Printing ***
      cout<<" - - - - - - - - - - - - - - - - - - - -  - - -  "<<endl;
      cout<<" - *Elec Variance parallel  [GeV]^2: "<< Var_L_e <<endl;
      cout<<" - *Elec Variance perpendicular  [GeV]^2: "<< Var_T_e <<endl;
      cout<<" - *Elec Covariance   [GeV]^2: "<< Cv_LT_e <<endl;


      //**************************************************//
      //***************   SIGNIFICANCE    *************** //
      if( Var_L != 0 )//&& Veto ==1) --> Filling only for events with good jets (&& is not needed) 
	{
	  Double_t Significance = MetSign::Significance_LT(Met,Var_L,Var_T,Cv_LT );
    
	  
	  Double_t Rho =  Cv_LT / sqrt( Var_L * Var_T ) ;
	  
	  // *** Printing ***
	  cout<<"  "<<endl;
	  cout<<"_________________________________________________"<<endl;
	  cout<<" - - - - - - - - - - - - - - - - - - - -  - - -  "<<endl;
	  cout<<" -> * Total Variance parallel  [GeV]^2: "<< Var_L <<endl;
	  cout<<" -> * Total Variance perpendicular  [GeV]^2: "<< Var_T <<endl;
	  cout<<" -> * Total Covariance   [GeV]^2: "<< Cv_LT <<endl;
	  cout<<" -> * Total correlationfactor rho: "<< Rho <<endl;
	  cout<<" -> * SIGNIFICANCE: "<< Significance <<endl;
	  cout<<" - - - - - - - - - - - - - - - - - - - -  - - -  "<<endl;
	  cout<<"_________________________________________________"<<endl;
	  
	  // *** HISTOGRAM FILLING ***
	  h_MetSig -> Fill( Significance  ,WEIGHT);
	  h_sqrtMS -> Fill( sqrt(Significance) ,WEIGHT);
	  h_MSgNum -> Fill( pow(Met,2)  ,WEIGHT);           
	  h_VarPar -> Fill( Var_L  ,WEIGHT);
	  h_VarPer -> Fill( Var_T  ,WEIGHT);
	  h_Covari -> Fill( Cv_LT  ,WEIGHT);
	  h_SqrMSg -> Fill( sqrt(Significance)  ,WEIGHT);
	  h_Met_sm -> Fill( Met/sqrt(met_sumet*0.001)  ,WEIGHT);
	  h_Met_Ht -> Fill( Met/sqrt(ht*0.001)  ,WEIGHT);
	  h_Met_si -> Fill( Met/sqrt(Var_L)  ,WEIGHT); 
	  h_MetPhi -> Fill( met_phi  ,WEIGHT);
	  h____Met -> Fill( Met ,WEIGHT);
	  h_softTm -> Fill( met_soft*0.001 ,WEIGHT);
	  h_metjet -> Fill( met_jet *0.001 ,WEIGHT);
	  h_MCchnu -> Fill( mc_channel_number,WEIGHT);
	  //h_dilePt -> Fill( dilep_pt[0]*0.001 ,WEIGHT);
	  // cout<<"..............................................................."<<endl;
	  // cout<<"  Met"<<Met<<"  met_soft"<<met_soft<<"  met_jet"<< met_jet<<endl;
	  // cout<<"..............................................................."<<endl;

	  // *** Ntuples Filling
	  MetSig_t = Significance ;
	  sqrtMS_t = sqrt(Significance) ;
	  MSgNum_t = pow(Met,2) ;
	  VarPar_t = Var_L ;
	  VarPer_t = Var_T ;
	  Covari_t = Cv_LT ;
	  Met_sm_t = Met/sqrt(met_sumet*0.001) ;
	  Met_Ht_t = Met/sqrt(ht*0.001) ;
	  Met_si_t = Met/sqrt(Var_L) ;
	  Jet_Nu_t = NJets ;
	  JetMet_t = met_jet *0.001 ;
	  Met____t = Met ;
	  MetPhi_t = met_phi ;
	  softTm_t = met_soft*0.001 ;
	  dilePt_t = dilep_pt[0]*0.001 ;
	  weight_t = weight               ;
	  MCweig_t = mc_event_weight      ;
	  pup_we_t = pileup_weight        ;
	  wshenj_t = weight_sherpa22_njets;
	  sf_jvt_t = sf_jvt               ;
	  sf_el__t = sf_el                ; 
	  sf_mu__t = sf_mu    	        ; 
	  sfmtri_t = sf_mu_trigger        ;
	  sfetri_t = sf_el_trigger        ;
	  TotWGT_t =  WEIGHT ;
	  EleVPa_t =  Var_L_e;
	  EleVPe_t = Var_T_e;
	  MuoVPa_t = Var_L_m ;
	  MuoVPe_t = Var_T_m;
	  JetVPa_t = Var_L_j;
	  JetVPe_t = Var_T_j;

	  if(fabs(Rho) >= 0.9)
	    {
	      //cout<<"NOT INVERTIBLE!!! -->rho="<<Rho<<endl;
	      Rho = 0; //In order to fill the Denominator distribution
	    }
	  // *** HISTOGRAM-TREE FILLING ***
	  h_MSgDen-> Fill( Var_L* ( 1 - pow(Rho,2) )  ,WEIGHT);
	  h_MSgRho -> Fill( Rho  , WEIGHT);
	  MSgDen_t = Var_L* ( 1 - pow(Rho,2) )  ;
	  MSgRho_t = Rho ;
	  //Loop over leptons
	  for (int i=0; i<n_lep; i++)
	    {
	      // *** HISTOGRAM-TREE FILLING ***
	      h_LepEta  -> Fill( lep_eta[i]     ,WEIGHT);
	      h_LepPhi  -> Fill( lep_phi[i]     ,WEIGHT);
	      h_LepPt_  -> Fill( lep_pt[i]*0.001  ,WEIGHT);
	      LepEta_t =  lep_eta[i]        ;
	      LepPhi_t =  lep_phi[i]        ;
	      LepPt__t =  lep_pt[i]*0.001   ;
	    }//loop over leptons
	  h_Lep_Nu  -> Fill( n_lep  ,WEIGHT);
	  Lep_Nu_t =  n_lep ;
	}//Var_L != 0 -  veto ==1
      else
      	{
      	  cout<<"No object passed selection in the event=> No Significance calculated"<<endl;
      	}//

    
      	  
      
    }//End baseline Event selection
  else
    {
      cout<< "***************Event #"<<fNumberOfEvents<< " did NOT Pass Selection*****************"<<endl;
      Jet_Pt_t = -100;
      JetEta_t = -100;
      JetPhi_t = -100;
      JetRes_t = -100;
      Jetpup_t = -100;
      Jet_OR_t = -100;
      JetJvt_t = -100;
      JetFMt_t = -100; 
      JetPre_t = -100;
      Muo_Pt_t = -100;
      MuoEta_t = -100;
      MuoPhi_t = -100;
      MuoRes_t = -100;
      MuoTyp_t = -100;
      Muocha_t = -100;
      MuoITO_t = -100;
      Ele_Pt_t = -100;
      EleEta_t = -100;
      ElePhi_t = -100;
      EleRes_t = -100;
      Ele_OR_t = -100;
      Jet_Nu_t = -100;
      JetMet_t = -100;
      Met____t = -100;
      MetPhi_t = -100;
      softTm_t = -100;
      dilePt_t = -100;
      weight_t = weight               ;
      MCweig_t = mc_event_weight      ;
      pup_we_t = pileup_weight        ;
      wshenj_t = weight_sherpa22_njets;
      sf_jvt_t = sf_jvt               ;
      sf_el__t = sf_el                ; 
      sf_mu__t = sf_mu    	        ; 
      sfmtri_t = sf_mu_trigger        ;
      sfetri_t = sf_el_trigger        ;
      TotWGT_t =  WEIGHT ;
     
      MetSig_t = -100;
      sqrtMS_t = -100;
      MSgNum_t = -100;
      VarPar_t = -100;
      VarPer_t = -100;
      Covari_t = -100;
      Met_sm_t = -100;
      Met_Ht_t = -100;
      Met_si_t = -100;
      EleVPa_t = -100;
      EleVPe_t = -100;
      MuoVPa_t = -100;
      MuoVPe_t = -100;
      JetVPa_t = -100;
      JetVPe_t = -100;
    }
  

  cout<<" Check->  Veto: "<<Veto<<endl;
  cout<<" ..........................................................................."<<endl;

   ++fNumberOfEvents;
   // Int_t EventMax = 100;
   // if(fNumberOfEvents>EventMax)
   //   {
   //     Abort("aborting...", kAbortProcess);
   //   }


  t_tree->Fill();
  
  return kTRUE;

}

void MetSign::SlaveTerminate()
{
}

void MetSign::Terminate()
{

  cout<<"TOTAL NUMBER OF EVENTS: "<<fNumberOfEvents<<endl;

  /// ** Output file ***********************************************************************
  //TFile f1("./HISTOS/Z_inclusive_J.root","RECREATE");// Nj
  TFile f1("./HISTOS/Zmumu_Pt700_1000_CVetoBVeto.root","RECREATE");//
  h_MetSig -> Write();
  h_sqrtMS -> Write();
  h_MSgNum -> Write();
  h_MSgDen -> Write();
  h_MSgRho -> Write();
  h_VarPar -> Write();
  h_VarPer -> Write();
  h_Covari -> Write();
  h_SqrMSg -> Write();
  h_Met_sm -> Write();
  h_Met_Ht -> Write();
  h_Met_si -> Write();
  h_Jet_Pt -> Write();
  h_JetEta -> Write();
  h_JetPhi -> Write();
  h_Jet_Nu -> Write();
  h_JetRes -> Write();
  h_MetPhi -> Write();
  h_LepEta -> Write();
  h_LepPhi -> Write();
  h_LepPt_ -> Write();
  h_Lep_Nu -> Write();
  h____Met -> Write();
  h_softTm -> Write();
  h_metjet -> Write();
  h_MCchnu -> Write();
  h_MCweig -> Write();
  h_dilePt -> Write();
  t_tree->Write();
  f1.Close(); 

  cout<<"File Written"<<endl;



  // TString treefile = "/data/atlas0/dportill/DougNtuples/mc15_13TeV.363387.Sherpa_NNPDF30NNLO_Zmumu_Pt2000_E_CMS_BFilter.merge.DAOD_SUSY5.e4716_s2726_r7725_r7676_p2666.root";
  // TFile *ftree = new TFile(treefile);
  // TTree *vesseltree = (TTree*)ftree->Get("CollectionTree");
  // //vesseltree->Print();

  // //create the new file
  // TFile *g = new TFile("./NTuples/test.root", "RECREATE");
  // TTree *newtree = vesseltree->CloneTree(0);

  // //int nentries = vesseltree->GetEntries();
  
}

double MetSign::Significance_LT(double Numerator, double var_parall, double var_perpen, double cov)
{
  Double_t rho = cov / sqrt( var_parall * var_perpen ) ;
  Double_t Significance = 0;
  if (abs( rho ) >= 0.9 )  //Cov Max not invertible -> Significance diverges
    {
      Significance = pow( Numerator , 2 ) / (  var_parall  ) ;
    }
  else
    {
      Significance = pow( Numerator , 2 ) / (  var_parall * ( 1 - pow(rho,2) ) ) ;
    }

  if( abs(Significance) >= 10e+15)
    {
      cout<<"warning -->"<< Significance <<endl;
    }
  
  return Significance;
}
tuple<double,double,double> MetSign::CovMatrixRotation(double var_x, double var_y, double cv_xy, double Phi)
{
  //Covariance matrix parallel and transverse to the Phi direction 
  Double_t V11 = pow(cos(Phi),2)*var_x + 2*sin(Phi)*cos(Phi)*cv_xy + pow(sin(Phi),2)*var_y;
  Double_t V22 = pow(sin(Phi),2)*var_x - 2*sin(Phi)*cos(Phi)*cv_xy + pow(cos(Phi),2)*var_y;
  Double_t V12 = pow(cos(Phi),2)*cv_xy -sin(Phi)*cos(Phi)*var_x + sin(Phi)*cos(Phi)*var_y - pow(sin(Phi),2)*cv_xy;   // rho is equal to one for just one jet
  return  make_tuple( V11, V22, V12);
}
TMatrixD MetSign::MatrixRotation(TMatrixD CovM_xy , double Phi)
{
  TMatrixD r(2,2);
  TMatrixD rInv(2,2);
  r(0,0) = cos(Phi); r(0,1) = sin(Phi); r(1,0) = -sin(Phi); r(1,1) = cos(Phi);
  rInv = r;
  rInv.Invert();
  //return rInv * CovM_xy * r;
  return r * CovM_xy * rInv;
}


double MetSign::GetPhiUnc(double jet_eta, double jet_phi,double jet_pt)
{
  unsigned xbin=0, ybin=0;
  if(-4.5<jet_eta && -3.8>=jet_eta)      xbin=1;
  else if(-3.8<jet_eta && -3.5>=jet_eta) xbin=2;
  else if(-3.5<jet_eta && -3.0>=jet_eta) xbin=3;
  else if(-3.0<jet_eta && -2.7>=jet_eta) xbin=4;
  else if(-2.7<jet_eta && -2.4>=jet_eta) xbin=5;
  else if(-2.4<jet_eta && -1.5>=jet_eta) xbin=6;
  else if(-1.5<jet_eta && -0.5>=jet_eta) xbin=7;
  else if(-0.5<jet_eta &&  0.0>=jet_eta) xbin=8;
  else if(0.0<jet_eta  &&  0.5>=jet_eta) xbin=9;
  else if(0.5<jet_eta  &&  1.5>=jet_eta) xbin=10;
  else if(1.5<jet_eta  &&  2.4>=jet_eta) xbin=11;
  else if(2.4<jet_eta  &&  2.7>=jet_eta) xbin=12;
  else if(2.7<jet_eta  &&  3.0>=jet_eta) xbin=13;
  else if(3.0<jet_eta  &&  3.5>=jet_eta) xbin=14;
  else if(3.5<jet_eta  &&  3.8>=jet_eta) xbin=15;
  else if(3.8<jet_eta                  ) xbin=16;

  ybin = jet_phi>0.0? int(jet_phi/0.4)+9:int(jet_phi/0.4)+8;

  // Stored as bin content = Mean, error = RMS, add them in quadrature
  if(jet_pt<50.0)
    return sqrt(h_phi_reso_pt20->GetBinContent(xbin, ybin)*h_phi_reso_pt20->GetBinContent(xbin, ybin)+h_phi_reso_pt20->GetBinError(xbin, ybin)*h_phi_reso_pt20->GetBinContent(xbin, ybin));
  else if(jet_pt<100.0)
    return sqrt(h_phi_reso_pt50->GetBinContent(xbin, ybin)*h_phi_reso_pt50->GetBinContent(xbin, ybin)+h_phi_reso_pt50->GetBinError(xbin, ybin)*h_phi_reso_pt50->GetBinContent(xbin, ybin));
  return sqrt(h_phi_reso_pt100->GetBinContent(xbin, ybin)*h_phi_reso_pt100->GetBinContent(xbin, ybin)+h_phi_reso_pt100->GetBinError(xbin, ybin)*h_phi_reso_pt100->GetBinContent(xbin, ybin));
}


double MetSign::GetXSxFiltEff(int mc_id ) // Xsec[pb]  * kFact * filtereff
{
  //Zmumu
  if(mc_id==363364) return  2077.0     *  0.9751  *   0.811 / 1.99687e+06;  // -> DerivationStat_Weights->GetBinContent(1)//1.167897625e+06;     -> Nov 23//2962000. ;//AOD AMI 
  if(mc_id==363365) return  2075.9     *  0.9751  *   0.1186/ 1.16676e+06;  // -> DerivationStat_Weights->GetBinContent(1)//8.034832500e+05;     -> Nov 23//1971000. ;//AOD AMI 
  if(mc_id==363366) return  2077.6     *  0.9751  *   0.0704/ 1.05937e+06;  // -> DerivationStat_Weights->GetBinContent(1)//6.958489375e+05;     -> Nov 23//1969000. ;//AOD AMI 
  if(mc_id==363367) return  71.72      *  0.9751  *   0.6669/ 572819.    ;  // -> DerivationStat_Weights->GetBinContent(1)//5.163017500e+05;     -> Nov 23//1237000. ;//AOD AMI 
  if(mc_id==363368) return  71.743     *  0.9751  *   0.2003/ 401147.    ;  // -> DerivationStat_Weights->GetBinContent(1)//3.690490625e+05;     -> Nov 23// 788000. ;//AOD AMI 
  if(mc_id==363369) return  71.574     *  0.9751  *   0.1276/ 544569.    ;  // -> DerivationStat_Weights->GetBinContent(1)//5.071501875e+05;     -> Nov 23//3236000. ;//AOD AMI 
  if(mc_id==363370) return  11.105     *  0.9751  *   0.6258/ 436614.    ;  // -> DerivationStat_Weights->GetBinContent(1)//4.1741959375e+05;    -> Nov 23// 789000. ;//AOD AMI 
  if(mc_id==363371) return  11.099     *  0.9751  *   0.2225/ 509500.    ;  // -> DerivationStat_Weights->GetBinContent(1)//4.9134059375e+05;    -> Nov 23// 791000. ;//AOD AMI 
  if(mc_id==363372) return  11.078     *  0.9751  *   0.1442/ 675932.    ;  // -> DerivationStat_Weights->GetBinContent(1)//6.5679400000e+05;    -> Nov 23//2956000. ;//AOD AMI 
  if(mc_id==363373) return  0.83396    *  0.9751  *   0.5923/ 317060.    ;  // -> DerivationStat_Weights->GetBinContent(1)//3.1410346875e+05;    -> Nov 23// 497000. ;//AOD AMI 
  if(mc_id==363374) return  0.83155    *  0.9751  *   0.2408/ 223157.    ;  // -> DerivationStat_Weights->GetBinContent(1)//2.2119806250e+05;    -> Nov 23// 297000. ;//AOD AMI 
  if(mc_id==363375) return  0.83216    *  0.9751  *   0.1523/ 190941.    ;  // -> DerivationStat_Weights->GetBinContent(1)//1.89595640625e+05;   -> Nov 23// 230000. ;//AOD AMI 
  if(mc_id==363376) return  0.053138   *  0.9751  *   0.5836/ 203933.    ;  // -> DerivationStat_Weights->GetBinContent(1)//2.03707656250e+05;   -> Nov 23// 293000. ;//AOD AMI 
  if(mc_id==363377) return  0.052965   *  0.9751  *   0.2511/ 156399.    ;  // -> DerivationStat_Weights->GetBinContent(1)//1.56253734375e+05;   -> Nov 23// 195000. ;//AOD AMI 
  if(mc_id==363378) return  0.053414   *  0.9751  *   0.147 / 4534.20    ;  // -> DerivationStat_Weights->GetBinContent(1)//4.53520312500e+03;   -> Nov 23//   5000. ;//AOD AMI 
  if(mc_id==363379) return  0.0095435  *  0.9751  *   0.5764/ 36438.9    ;  // -> DerivationStat_Weights->GetBinContent(1)//3.640801953125e+04;  -> Nov 23 //  49000.;//AOD AMI 
  if(mc_id==363380) return  0.0095439  *  0.9751  *   0.2548/ 24678.7    ;  // -> DerivationStat_Weights->GetBinContent(1)//2.4653939453125e+04; -> Nov 23 //  29000.;//AOD AMI 
  if(mc_id==363381) return  0.0095915  *  0.9751  *   0.1599/ 26621.7    ;  // -> DerivationStat_Weights->GetBinContent(1)//2.6591181640625e+04; -> Nov 23 //  29000.;//AOD AMI 
  if(mc_id==363382) return  0.0012698  *  0.9751  *   0.5579/ 22557.0    ;  // -> DerivationStat_Weights->GetBinContent(1)//2.2509687500000e+04; -> Nov 23 //  30000.;//AOD AMI 
  if(mc_id==363383) return  0.0012528  *  0.9751  *   0.2611/ 16776.2    ;  // -> DerivationStat_Weights->GetBinContent(1)//1.6758078125000e+04; -> Nov 23 //  20000.;//AOD AMI 
  if(mc_id==363384) return  0.0012306  *  0.9751  *   0.1539/ 17662.0    ;  // -> DerivationStat_Weights->GetBinContent(1)//1.7643875000000e+04; -> Nov 23 //  20000.;//AOD AMI 
  if(mc_id==363385) return  4.4846e-06 *  0.9751  *   0.5606/ 11097.3    ;  // -> DerivationStat_Weights->GetBinContent(1)//1.1052402343750e+04; -> Nov 23 //  20000.;//AOD AMI 
  if(mc_id==363386) return  5.2027e-06 *  0.9751  *   0.272 / 8526.28    ;  // -> DerivationStat_Weights->GetBinContent(1)//8.4850507812500e+03; -> Nov 23 //  10000.;//AOD AMI 
  if(mc_id==363387) return  4.7171e-06 *  0.9751  *   0.1455/ 7544.75    ;  // -> DerivationStat_Weights->GetBinContent(1)//7.52532373046875e+03;-> Nov 23  // 10000 ;//AOD AMI 
  //Zee
  if(mc_id==363388) return 2076.4      *  0.9751  *   0.8107  ;
  if(mc_id==363389) return 2078.5      *  0.9751  *   0.1191  ;
  if(mc_id==363390) return 2074.0*       0.9751*      0.0;       
  if(mc_id==363391) return 71.681      *  0.9751  *   0.6694  ;
  if(mc_id==363392) return 71.341      *  0.9751  *   0.202   ;
  if(mc_id==363393) return 71.744      *  0.9751  *   0.1267  ;
  if(mc_id==363394) return 11.095      *  0.9751  *   0.6276  ;
  if(mc_id==363395) return 11.051      *  0.9751  *   0.2252  ;
  if(mc_id==363396) return 10.952      *  0.9751  *   0.146   ;
  if(mc_id==363397) return 0.83477     *  0.9751  *   0.5985  ;
  if(mc_id==363398) return 0.8291      *  0.9751  *   0.2453  ;
  if(mc_id==363399) return 0.83251     *  0.9751  *   0.1515  ;
  if(mc_id==363400) return 0.053181    *  0.9751  *   0.5545  ;
  if(mc_id==363401) return 0.052912    *  0.9751  *   0.2507  ;
  if(mc_id==363402) return 0.052964    *  0.9751  *   0.1543  ;
  if(mc_id==363403) return 0.0096374   *  0.9751  *   0.5762  ;
  if(mc_id==363404) return 0.009431    *  0.9751  *   0.2611  ;
  if(mc_id==363405) return 0.0096257   *  0.9751  *   0.1593  ;
  if(mc_id==363406) return 0.0012529   *  0.9751  *   0.5712  ;
  if(mc_id==363407) return 0.0012711   *  0.9751  *   0.2685  ;
  if(mc_id==363408) return 0.0012502   *  0.9751  *   0.1554  ;
  if(mc_id==363409) return 4.9183e-06  *  0.9751  *   0.5264  ;
  if(mc_id==363410) return 4.7405e-06  *  0.9751  *   0.2635  ;
  if(mc_id==363411) return 4.6041e-06  *  0.9751  *   0.1556  ;
  
  // ZZ->llVV
  if(mc_id==361068) return 14.022*0.91/4.04033e+06;
  
  // ttbar + 1lep
  if(mc_id==410000) return 1.;//696.11*0.54383; 
}