#define TreeReader_cxx
#include "TreeReader.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include "include/Muon.h"
#include "include/Track.h"
#include "include/PrimaryVertex.h"
#include "include/PFJet.h"
#include <iostream>
#include <cmath>
#include <TLorentzVector.h>
#include <stdio.h>
#include <TBranch.h>

void TreeReader::Loop()
{
//   In a ROOT session, you can do:
//      Root > .L TreeReader.C
//      Root > TreeReader 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 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

createHistograms();
  
  //Set branches of Signal and Bkg trees

   treeS = new TTree("Signal",     "signal tracks");
   treeB = new TTree("Background", "Background tracks");

//    treeS->Branch("mva_track_firstHit_x",        &mva_track_firstHit_x);//exclude variables that higly depend on the flight distance
//    treeS->Branch("mva_track_firstHit_y",        &mva_track_firstHit_y);//exclude variables that higly depend on the flight distance
//    treeS->Branch("mva_track_firstHit_z",        &mva_track_firstHit_z);//exclude variables that higly depend on the flight distance
//    treeS->Branch("mva_track_dxy",               &mva_track_dxy);
//    treeS->Branch("mva_track_dxyError",          &mva_track_dxyError);
//    treeS->Branch("mva_track_dz",                &mva_track_dz);
//    treeS->Branch("mva_track_dzError",           &mva_track_dzError);
//    treeS->Branch("mva_track_pt",                &mva_track_pt);
//    treeS->Branch("mva_track_eta",               &mva_track_eta);
//    treeS->Branch("mva_track_phi",               &mva_track_phi);
//    treeS->Branch("mva_track_nchi2",             &mva_track_nchi2);
//    treeS->Branch("mva_track_nhits",             &mva_track_nhits);
//    treeS->Branch("mva_track_numberOfValidPixelHits", &mva_track_numberOfValidPixelHits);//exclude variables that higly depend on the flight distance
//    treeS->Branch("mva_track_numberOfValidStripHits", &mva_track_numberOfValidStripHits);//exclude variables that higly depend on the flight distance
//    treeS->Branch("mva_track_isinjet",                &mva_track_isinjet);
//    treeS->Branch("mva_track_algo",              &mva_track_algo);
//    //added by Paul
//    treeS->Branch("mva_ntrk10",              &mva_ntrk10);
//    treeS->Branch("mva_ntrk20",              &mva_ntrk20);
//    treeS->Branch("mva_ntrk30",              &mva_ntrk30);
//    treeS->Branch("mva_ntrk40",              &mva_ntrk40);
//    treeS->Branch("mva_ntrk1020",              &mva_ntrk1020);
//    treeS->Branch("mva_ntrk2030",              &mva_ntrk2030);
//    treeS->Branch("mva_drSig",              &mva_drSig);
//    treeS->Branch("mva_dzSig",              &mva_dzSig);
//    treeS->Branch("mva_ddSig",              &mva_ddSig);
//    treeS->Branch("mva_ValTIBHit",              &mva_ValTIBHit);//exclude variables that higly depend on the flight distance
//    treeS->Branch("mva_ValTOBHit",              &mva_ValTOBHit);//exclude variables that higly depend on the flight distance
//    treeS->Branch("mva_ValPixBarHit",              &mva_ValPixBarHit);//exclude variables that higly depend on the flight distance
//    treeS->Branch("mva_nValTECHHit",              &mva_nValTECHHit);//exclude variables that higly depend on the flight distance
//    treeS->Branch("mva_ntrk40XX",              &mva_ntrk40XX);
//    treeS->Branch("mva_track_dR",&mva_track_dR);
//    treeS->Branch("mva_track_dRmax",&mva_track_dRmax);
// //    treeS->Branch("mva_track_isTight",&mva_track_isTight);
//    treeS->Branch("mva_track_isHighQuality",&mva_track_isHighQuality);
// // treeS->Branch("mva_track_region",&mva_track_region);
// // treeS->Branch("mva_track_hitpattern",&mva_track_hitpattern);
//    //Added 10/02/2023 on MiniAOD datatier (1st year thesis)
//    //HT
//    treeS->Branch("mva_HT",&mva_HT);
   //PFMet
   treeS->Branch("mva_MET_et",&mva_MET_et);
   //jet-------------------------
//    treeS->Branch("mva_jet_nbr",&mva_jet_nbr);
//    treeS->Branch("mva_jet_btag",&mva_jet_btag);
   //tracks---------------
   treeS->Branch("mva_track_lost",&mva_track_lost);
   //muon---------
   treeS->Branch("mva_muon_pt",&mva_muon_pt);
   treeS->Branch("mva_muon_eta",&mva_muon_eta);
   treeS->Branch("mva_muon_phi",&mva_muon_phi);
   treeS->Branch("mva_muon_energy",&mva_muon_energy);
   treeS->Branch("mva_muon_isGlobal",&mva_muon_isGlobal);
   treeS->Branch("mva_muon_isoR3",&mva_muon_isoR3);
   treeS->Branch("mva_muon_Mmumu",&mva_muon_Mmumu);

//    ///BKG
//    treeB->Branch("mva_track_firstHit_x",        &mva_track_firstHit_x);//exclude variables that higly depend on the flight distance
//    treeB->Branch("mva_track_firstHit_y",        &mva_track_firstHit_y);//exclude variables that higly depend on the flight distance
//    treeB->Branch("mva_track_firstHit_z",        &mva_track_firstHit_z);//exclude variables that higly depend on the flight distance
//    treeB->Branch("mva_track_dxy",               &mva_track_dxy);
//    treeB->Branch("mva_track_dxyError",          &mva_track_dxyError);
//    treeB->Branch("mva_track_dz",                &mva_track_dz);
//    treeB->Branch("mva_track_dzError",           &mva_track_dzError);
//    treeB->Branch("mva_track_pt",                &mva_track_pt);
//    treeB->Branch("mva_track_eta",               &mva_track_eta);
//    treeB->Branch("mva_track_phi",               &mva_track_phi);
//    treeB->Branch("mva_track_nchi2",             &mva_track_nchi2);
//    treeB->Branch("mva_track_nhits",             &mva_track_nhits);//exclude variables that higly depend on the flight distance
//    treeB->Branch("mva_track_numberOfValidPixelHits", &mva_track_numberOfValidPixelHits);//exclude variables that higly depend on the flight distance
//    treeB->Branch("mva_track_numberOfValidStripHits", &mva_track_numberOfValidStripHits);//exclude variables that higly depend on the flight distance
//    treeB->Branch("mva_track_isinjet",                &mva_track_isinjet);
//    treeB->Branch("mva_track_algo",              &mva_track_algo);
// //    treeB->Branch("mva_track_isTight",&mva_track_isTight);
//    treeB->Branch("mva_track_isHighQuality",&mva_track_isHighQuality);


//    treeB->Branch("mva_ntrk10",              &mva_ntrk10);
//    treeB->Branch("mva_ntrk20",              &mva_ntrk20);
//    treeB->Branch("mva_ntrk30",              &mva_ntrk30);
//    treeB->Branch("mva_ntrk40",              &mva_ntrk40);
//    treeB->Branch("mva_ntrk1020",              &mva_ntrk1020);
//    treeB->Branch("mva_ntrk2030",              &mva_ntrk2030);
//    treeB->Branch("mva_drSig",              &mva_drSig);
//    treeB->Branch("mva_dzSig",              &mva_dzSig);
//    treeB->Branch("mva_ddSig",              &mva_ddSig);
//    treeB->Branch("mva_ValTIBHit",              &mva_ValTIBHit);//exclude variables that higly depend on the flight distance
//    treeB->Branch("mva_ValTOBHit",              &mva_ValTOBHit);//exclude variables that higly depend on the flight distance
//    treeB->Branch("mva_ValPixBarHit",              &mva_ValPixBarHit);//exclude variables that higly depend on the flight distance
//    treeB->Branch("mva_nValTECHHit",              &mva_nValTECHHit);//exclude variables that higly depend on the flight distance
//    treeB->Branch("mva_ntrk40XX",              &mva_ntrk40XX);
//    treeB->Branch("mva_track_dR",&mva_track_dR);//new
//    treeB->Branch("mva_track_dRmax",&mva_track_dRmax);

//    //Added 10/02/2023 on MiniAOD datatier (1st year thesis)
//    //HT
//    treeB->Branch("mva_HT",&mva_HT);
   //PFMet
   treeB->Branch("mva_MET_et",&mva_MET_et);
   //jet-------------------------
//    treeB->Branch("mva_jet_nbr",&mva_jet_nbr);
   //tracks---------------
//    treeB->Branch("mva_track_lost",&mva_track_lost);
   //muon---------
   treeB->Branch("mva_muon_pt",&mva_muon_pt);
   treeB->Branch("mva_muon_eta",&mva_muon_eta);
   treeB->Branch("mva_muon_phi",&mva_muon_phi);
   treeB->Branch("mva_muon_energy",&mva_muon_energy);
   treeB->Branch("mva_muon_isGlobal",&mva_muon_isGlobal);
   treeB->Branch("mva_muon_isoR3",&mva_muon_isoR3);
   treeB->Branch("mva_muon_Mmumu",&mva_muon_Mmumu);
   treeB->Branch("mva_jet_btag",&mva_jet_btag);
   //---------------------------------//

   if (fChain == 0) return;

   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;

      LoadTreeSizes(theFormat);//this line is bugged
    //   track_size = theFormat.tree_track_pt->size();
    //   muon_size = theFormat.tree_muon_pt->size();

      if (jentry%500 == 0) cout << "# of processed events " << jentry << endl;
	std::cout<< "jentry : "<< jentry <<std::endl;

         //MUONS/////////////////////////////////////////////////////////
    mva_MET_et = -10 ;
    mva_MET_et = theFormat.tree_PFMet_et;
    nSBr1++;

    treeS->GetBranch("mva_MET_et")->Fill();

    treeB->GetBranch("mva_MET_et")->Fill();

    int muon_idx=-1;
    // bool ZMuRec = false;
    // int PFiso, MVAiso, CutBasedId;
    float dVr, dVz;
    float Mmumu = 0.;
    int nmurec = 0, imu1 = -1, imu2 = -1;
    float mupt1, mueta1, muphi1, mupt2, mueta2, muphi2;
    TLorentzVector vmuon[2], v1,v2,v;
    float mu_mass = 0.1057;
    PrimaryVertex thepv; thepv.read(0, &theFormat);
    float PVx = thepv.mini_x();
    float PVy = thepv.mini_y();
    float PVz = thepv.mini_z();
	Muon themuon;
	Muon themuon2;
    bool Minvar = true;
   nSBr0 += muon_size;
//    treeS->GetBranch("mva_muon_pt")->SetEntries(nSBr0);
//    treeS->GetBranch("mva_muon_eta")->SetEntries(nSBr0);
//    treeS->GetBranch("mva_muon_phi")->SetEntries(nSBr0);
//    treeS->GetBranch("mva_muon_energy")->SetEntries(nSBr0);
//    treeS->GetBranch("mva_muon_isGlobal")->SetEntries(nSBr0);
//    treeS->GetBranch("mva_muon_isoR3")->SetEntries(nSBr0);

//    treeB->GetBranch("mva_muon_pt")->SetEntries(nSBr0);
//    treeB->GetBranch("mva_muon_eta")->SetEntries(nSBr0);
//    treeB->GetBranch("mva_muon_phi")->SetEntries(nSBr0);
//    treeB->GetBranch("mva_muon_energy")->SetEntries(nSBr0);
//    treeB->GetBranch("mva_muon_isGlobal")->SetEntries(nSBr0);
//    treeB->GetBranch("mva_muon_isoR3")->SetEntries(nSBr0);

    for (int midx=0; midx<muon_size ; midx++)//loop on muons
    {
        
    	themuon.read(midx, &theFormat);
        mva_muon_pt=themuon.mini_pt();
        mva_muon_eta=themuon.mini_eta();
        mva_muon_phi=themuon.mini_phi();
        mva_muon_energy=themuon.mini_energy();
        mva_muon_isGlobal=themuon.mini_isGlobal();
        mva_muon_isoR3=themuon.mini_isoR3();

        treeS->GetBranch("mva_muon_pt")->Fill();
        treeS->GetBranch("mva_muon_eta")->Fill();
        treeS->GetBranch("mva_muon_phi")->Fill();
        treeS->GetBranch("mva_muon_energy")->Fill();
        treeS->GetBranch("mva_muon_isGlobal")->Fill();
        treeS->GetBranch("mva_muon_isoR3")->Fill();
  
        treeB->GetBranch("mva_muon_pt")->Fill();
        treeB->GetBranch("mva_muon_eta")->Fill();
        treeB->GetBranch("mva_muon_phi")->Fill();
        treeB->GetBranch("mva_muon_energy")->Fill();
        treeB->GetBranch("mva_muon_isGlobal")->Fill();
        treeB->GetBranch("mva_muon_isoR3")->Fill();

        if ( themuon.mini_pt() > 3. ) 
            {
                 muon_idx++;
                int nmu = muon_size;
                
                //if ( !themuon.mini_isGlobalMuon() ) continue; //Need global muons
                if ( themuon.mini_isGlobal() )
                {
                    mupt1  = themuon.mini_pt();
                    //if ( mupt1 < 10. ) continue; //Zmu filter
                    if ( mupt1 > 10. )
                    {
                        dVr = TMath::Sqrt( (themuon.mini_vx()-PVx)*(themuon.mini_vx()-PVx) + (themuon.mini_vy()-PVy)*(themuon.mini_vy()-PVy) );
                        dVz = themuon.mini_vz()-PVz;
                        // if ( dVr > 0.1 || abs(dVz) > 0.2 ) continue;// on veut un bon fit pour nos PV d'ou un seuil maximum sur les distances
                        if ( dVr < 0.1 && abs(dVz) < 0.2 )
                            {
                              mueta1 = themuon.mini_eta();
                              muphi1 = themuon.mini_phi();
        
                                v1.SetPtEtaPhiM(mupt1,mueta1,muphi1,mu_mass);
                                vmuon[0].SetPtEtaPhiM(mupt1,mueta1,muphi1,0);
                                nmurec++;
	                            //if ( muon_idx == nmu-1 ) continue;
	                            if ( muon_idx < nmu-1 ) 
	                                {
                                        for ( int muon2_idx=muon_idx+1;muon2_idx<nmu;muon2_idx++) 
                                            {       // Loop on other reco muons

	                                            themuon2.read(muon2_idx, &theFormat);
                                                
                                                if ( themuon2.mini_isGlobal() )
                                                {
                                                    
                                                    if ( themuon.mini_charge() != themuon2.mini_charge() )
                                                        {
                                                            dVr = TMath::Sqrt( (themuon2.mini_vx()-PVx)*(themuon2.mini_vx()-PVx) + (themuon2.mini_vy()-PVy)*(themuon2.mini_vy()-PVy) );
                                                            dVz = themuon2.mini_vz()-PVz;
                                                            
                                                            if ( dVr < 0.1 && abs(dVz) < 0.2 ) 
                                                                {
                                                                    mupt2  = themuon2.mini_pt();
                                                                    
                                                                    if ( mupt2 > 10. )
                                                                        {
                                                                           
                                                                             if ( mupt1 > 28. || mupt2 > 28. )
                                                                             {
                                                                                mueta2 = themuon2.mini_eta();
                                                                                muphi2 = themuon2.mini_phi();
                                                                                //muon2_idx++;
                                                                                v2.SetPtEtaPhiM(mupt2,mueta2,muphi2,mu_mass);
                                                                                vmuon[1].SetPtEtaPhiM(mupt2,mueta2,muphi2,0);
                                                                                v = v1 + v2;
                                                                                if ( v.Mag() > Mmumu )
                                                                                    { //Mag pour masse invariante (magnitude)
                                                                                    Mmumu = v.Mag();
                                                                                    if ( Mmumu > 60. ) 
                                                                                        {
                                                                                        mva_muon_Mmumu=Mmumu;
                                                                                        Minvar =true;
                                                                                        nSBr2=nSBr2+1;
                                                                                        
                                                                                        treeS->GetBranch("mva_muon_Mmumu")->Fill();
                                                                                        treeB->GetBranch("mva_muon_Mmumu")->Fill();
                                                                                        imu1 = muon_idx;
                                                                                        imu2 = muon2_idx;
                                                                                        }
                                                                                    }
                                                                             }

                                                                        }

                                                                }

                                                        }

                                                }

                                           }//second loop on muons
	                                }

                            }

                    }

                }
            }
    }//end loop on muons

     if (imu1 != -1 && imu2 != -1 && imu1 <muon_size && imu2 <muon_size) 
     	{
        themuon.read(imu1, &theFormat);
        themuon2.read(imu2, &theFormat);

       if ( themuon2.mini_pt() > themuon.mini_pt()  && Minvar==true) 
       		{
            int imu0 = imu2;
            imu2 = imu1; //muons reco with imu1 having the highest pt
            imu1 = imu0;
            v = vmuon[0];
            vmuon[0] = vmuon[1];
            vmuon[1] = v;
        	}
	}

   }
   writeHistograms();
   treeS->SetEntries(nSBr0);
   treeB->SetEntries(nSBr0);
   treeS->Write();
   treeB->Write();


}