/// JpsiDumper
#define JpsiDumper_cxx
#include "JpsiDumper.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TFile.h>
#include <iostream>
#include "families.h"
#include <fstream>
#include <algorithm>
#include <map>
#include "math.h"


//////////////////////////////////////
//            F L A G S             //
//////////////////////////////////////
Bool_t doTrackTrack = kFALSE;
Bool_t doTrackMuon  = kFALSE;
Bool_t doMuonMuon   = kTRUE;
Bool_t doMuid       = kTRUE;
Bool_t doStaco      = kTRUE;
Bool_t doSMuon      = kFALSE; // activate flat Tree for single muon
Bool_t doDimuon     = kFALSE; // activate flat Tree for dimuons
Bool_t doMultiMuon  = kTRUE;  // activate flat Tree for multi muons
Float_t barrelEta = 1.05;
Float_t minInvmass = 2000.;
Float_t maxInvmass = 4000.;
Float_t muonPDGmass = 105.658367;
Bool_t IS_DEBUG = 1; 
Float_t PiGreek = 3.14159;
//////////////////////////////////////

typedef  pair<int, float> my_pair;

bool sort_pred(const my_pair& left, const my_pair& right)
{
return left.second < right.second;
}

void JpsiDumper::LoadROIMap(){
  if(IS_DEBUG) file << " ** Load Roi Map "<<endl;
  ifstream fileRoI;
  fileRoI.open("RoIMap.txt");
  if(!fileRoI) file << "Error: file could not be opened" << endl;
  
  int SubSysID;
  int sectorID;
  int roINumber;
  float etamin,etamax,phimin,phimax;
  
  //look on file
  while (! fileRoI.eof())
    {
      //if (!fileRoI.good()) break;
      
      fileRoI >> SubSysID>>sectorID>>roINumber>>etamin
	      >>etamax>>phimin>>phimax;
      
      int mySecId=sectorID+32*SubSysID;
      file  << "  sec = "<< mySecId  
	    << "  roi = "<< roINumber
	    << "  etamin = "<<etamin
	    << " etamax = "<<etamax   
	    << "  phimin = "<<phimin
	    << " phimax = "<<phimax   
	    <<endl;
      
      muctpiEta_min[mySecId][roINumber]= etamin;
      muctpiEta_max[mySecId][roINumber]= etamax;
      muctpiPhi_min[mySecId][roINumber]= phimin;
      muctpiPhi_max[mySecId][roINumber]= phimax;
      
    }
}


void JpsiDumper::Begin()
{
  file << " Beginning analysis..." << std::endl;
}

void JpsiDumper::Loop()
{

  file.open("log.txt");
 
  LoadROIMap();
 
  outputfile = new TFile("H200.root", "RECREATE");
  
  L1_items["L1_MU0"]=0;
  L1_items["L1_MU6"]=0;
  L1_items["L1_MU10"]=0;
  L1_items["L1_MU0_COMM"]=0;

  EF_items["EF_mu4"]=0;
  EF_items["EF_mu6"]=0;
  EF_items["EF_mu10"]=0;
  EF_items["EF_mu4_MSonly_MB2_noL2_EFFS"]=0;

   // Here the definition of variables in DIMUON DUMPER TREE
  file << " Allocating the JpsiDIMUON output tree..." << std::endl;
  dump = new TTree("JpsiTree", "Jpsi flat tree");
  dump->Branch("runNumber", &m_runNumber, "runNumber/I");
  dump->Branch("eventNumber ", &m_eventNumber, "eventNumber/I");
  dump->Branch("lumiBlock", &m_lumiBlock, "lumiBlock/I");
  dump->Branch("bcId", &m_bcId, "bcId/I");
  dump->Branch("invmass_novtx", &m_invmass_novtx, "invmass_novtx/F");
  dump->Branch("invmass", &m_invmass, "invmass/F");
  dump->Branch("q", &m_q, "q/I");
  dump->Branch("pt", &m_pt, "pt/F");
  dump->Branch("deltaEta", &m_deltaEta, "deltaEta/F");
  dump->Branch("deltaPhi", &m_deltaPhi, "deltaPhi/F");
  dump->Branch("deltaPt", &m_deltaPt, "deltaPt/F");
  dump->Branch("family", &m_family, "family/I");
  dump->Branch("etarange", &m_etarange, "etarange/I"); /// start from here
  
  /*  
      dump->Branch("nmuctpi", &m_nmuctpi, "nmuctpi/I"); //
   dump->Branch("mu1index", &m_mu1index, "mu1index/I");
   
   dump->Branch("mu1hasMSinfo", &m_mu1hasMSinfo, "mu1hasMSinfo/I"); //fra
   dump->Branch("mu1hasExtra", &m_mu1hasExtra, "mu1hasExtra/I"); //fra
   dump->Branch("mu1RoiHitsEta", &m_mu1RoiHitsEta, "mu1RoiHitsEta/F"); //fra
   dump->Branch("mu1RoiHitsPhi", &m_mu1RoiHitsPhi, "mu1RoiHitsPhi/F"); //fra
   dump->Branch("mu1roisysId", &m_mu1roisysId,   "mu1roisysId/I");       // fra
   dump->Branch("mu1dr1l1", &m_mu1dr1l1,   "mu1dr1l1/F");       // fra
   dump->Branch("mu1dr1l1e", &m_mu1dr1l1e,   "mu1dr1l1e/F");       // fra
   dump->Branch("mu1roi1l1", &m_mu1roi1l1, "mu1roi1l1/I");    // fra
   dump->Branch("mu1thre1l1",&m_mu1thre1l1,"mu1thre1l1/I"); // fra
   dump->Branch("mu1dr2l1", &m_mu1dr2l1, "mu1dr2l1/F");    // fra
   dump->Branch("mu1dr2l1e", &m_mu1dr2l1e, "mu1dr2l1e/F");    // fra
   dump->Branch("mu1roi2l1", &m_mu1roi2l1, "mu1roi2l1/I");    // fra
   dump->Branch("mu1thre2l1",&m_mu1thre2l1,"mu1thre2l1/I"); // fra
   dump->Branch("mu1newchi",&m_mu1newchi,"mu1newchi/F"); // fra
   dump->Branch("mu1makesStacoCombined", &m_mu1makesStacoCombined, "mu1makesStacoCombined/I");
   dump->Branch("mu1makesMuidCombined", &m_mu1makesMuidCombined, "mu1makesMuidCombined/I");
  */
  
   dump->Branch("mu1index", &m_mu1index, "mu1index/I");
   dump->Branch("mu1author", &m_mu1author, "mu1author/I");
   dump->Branch("mu1pt", &m_mu1pt, "mu1pt/F");
   dump->Branch("mu1eta", &m_mu1eta, "mu1eta/F");
   dump->Branch("mu1etaex", &m_mu1etaex, "mu1etaex/F");
   dump->Branch("mu1etae", &m_mu1etae, "mu1etae/F");
   dump->Branch("mu1phie", &m_mu1phie, "mu1phie/F");
   dump->Branch("mu1etas", &m_mu1etas, "mu1etas/F");
   dump->Branch("mu1phi", &m_mu1phi, "mu1phi/F");
   dump->Branch("mu1phiex", &m_mu1phiex, "mu1phiex/F");
   dump->Branch("mu1phis", &m_mu1phis, "mu1phis/F");
   dump->Branch("mu1pt_vtx", &m_mu1pt_vtx, "mu1pt_vtx/F");
   dump->Branch("mu1eta_vtx", &m_mu1eta_vtx, "mu1eta_vtx/F");
   dump->Branch("mu1phi_vtx", &m_mu1phi_vtx, "mu1phi_vtx/F");
   dump->Branch("mu1pt_id", &m_mu1pt_id, "mu1pt_id/F");
   dump->Branch("mu1eta_id", &m_mu1eta_id, "mu1eta_id/F");
   dump->Branch("mu1phi_id", &m_mu1phi_id, "mu1phi_id/F");
   dump->Branch("mu1IDTrackChi2Reduct", &m_mu1IDTrackChi2Reduct, "mu1IDTrackChi2Reduct/F");
   dump->Branch("mu1IDProbChi2Fit", &m_mu1IDProbChi2Fit, "mu1IDProbChi2Fit/F");
   dump->Branch("mu1IDTrackDeltaChi2", &m_mu1IDTrackDeltaChi2, "mu1IDTrackDeltaChi2/F");

   dump->Branch("mu1_etIsol10", &m_mu1Etisol10, "mu1_etIsol10/F");   
   dump->Branch("mu1_etIsol20", &m_mu1Etisol20, "mu1_etIsol20/F");
   dump->Branch("mu1_etIsol30", &m_mu1Etisol30, "mu1_etIsol30/F");
   dump->Branch("mu1_etIsol40", &m_mu1Etisol40, "mu1_etIsol40/F");
   dump->Branch("mu1_ptIsol10", &m_mu1Ptisol10, "mu1_ptIsol10/F");
   dump->Branch("mu1_ptIsol20", &m_mu1Ptisol20, "mu1_ptIsol20/F");
   dump->Branch("mu1_ptIsol30", &m_mu1Ptisol30, "mu1_ptIsol30/F");
   dump->Branch("mu1_ptIsol40", &m_mu1Ptisol40, "mu1_ptIsol40/F");
   dump->Branch("mu1_nuIsol10", &m_mu1Nuisol10, "mu1_nuIsol10/F");
   dump->Branch("mu1_nuIsol20", &m_mu1Nuisol20, "mu1_nuIsol20/F");
   dump->Branch("mu1_nuIsol30", &m_mu1Nuisol30, "mu1_nuIsol30/F");
   dump->Branch("mu1_nuIsol40", &m_mu1Nuisol40, "mu1_nuIsol40/F");
   dump->Branch("mu1ProbChi2Fit", &m_mu1ProbChi2Fit, "mu1ProbChi2Fit/F");
   dump->Branch("mu1ProbChi2Match", &m_mu1ProbChi2Match, "mu1ProbChi2Match/F");
   dump->Branch("mu1PtRatio_ex", &m_mu1PtRatio_ex, "mu1PtRatio_ex/F");
   dump->Branch("mu1DeltaPtSign_ex", &m_mu1DeltaPtSign_ex, "mu1DeltaPtSign_ex/F");
   dump->Branch("mu1EnergyBalance", &m_mu1EnergyBalance, "mu1EnergyBalance/F");
   dump->Branch("mu1d0", &m_mu1d0, "mu1d0/F");
   dump->Branch("mu1z0", &m_mu1z0, "mu1z0/F");
   dump->Branch("mu1d0PVCorr", &m_mu1d0PVCorr, "mu1d0PVCorr/F");
   dump->Branch("mu1z0PVCorr", &m_mu1z0PVCorr, "mu1z0PVCorr/F");
   dump->Branch("mu1q", &m_mu1q, "mu1q/I");
   // mu2 vars
   dump->Branch("mu2index", &m_mu2index, "mu2index/I"); 
   dump->Branch("mu2author", &m_mu2author, "mu2author/I");
   dump->Branch("mu2pt", &m_mu2pt, "mu2pt/F");
   dump->Branch("mu2eta", &m_mu2eta, "mu2eta/F");
   dump->Branch("mu2etaex", &m_mu2etaex, "mu2etaex/F");
   dump->Branch("mu2etae", &m_mu2etae, "mu2etae/F");
   dump->Branch("mu2phie", &m_mu2phie, "mu2phie/F");
   dump->Branch("mu2etas", &m_mu2etas, "mu2etas/F");
   dump->Branch("mu2phi", &m_mu2phi, "mu2phi/F");
   dump->Branch("mu2phiex", &m_mu2phiex, "mu2phiex/F");
   dump->Branch("mu2phis", &m_mu2phis, "mu2phis/F");
   dump->Branch("mu2pt_vtx", &m_mu2pt_vtx, "mu2pt_vtx/F");
   dump->Branch("mu2eta_vtx", &m_mu2eta_vtx, "mu2eta_vtx/F");
   dump->Branch("mu2phi_vtx", &m_mu2phi_vtx, "mu2phi_vtx/F");
   dump->Branch("mu2pt_id", &m_mu2pt_id, "mu2pt_id/F");
   dump->Branch("mu2eta_id", &m_mu2eta_id, "mu2eta_id/F");
   dump->Branch("mu2phi_id", &m_mu2phi_id, "mu2phi_id/F");
   dump->Branch("mu2IDTrackChi2Reduct", &m_mu2IDTrackChi2Reduct, "mu2IDTrackChi2Reduct/F");
   dump->Branch("mu2IDProbChi2Fit", &m_mu2IDProbChi2Fit, "mu2IDProbChi2Fit/F");
   dump->Branch("mu2IDTrackDeltaChi2", &m_mu2IDTrackDeltaChi2, "mu2IDTrackDeltaChi2/F");

   dump->Branch("mu2_etIsol10", &m_mu2Etisol10, "mu2_etIsol10/F");   
   dump->Branch("mu2_etIsol20", &m_mu2Etisol20, "mu2_etIsol20/F");
   dump->Branch("mu2_etIsol30", &m_mu2Etisol30, "mu2_etIsol30/F");
   dump->Branch("mu2_etIsol40", &m_mu2Etisol40, "mu2_etIsol40/F");
   dump->Branch("mu2_ptIsol10", &m_mu2Ptisol10, "mu2_ptIsol10/F");
   dump->Branch("mu2_ptIsol20", &m_mu2Ptisol20, "mu2_ptIsol20/F");
   dump->Branch("mu2_ptIsol30", &m_mu2Ptisol30, "mu2_ptIsol30/F");
   dump->Branch("mu2_ptIsol40", &m_mu2Ptisol40, "mu2_ptIsol40/F");
   dump->Branch("mu2_nuIsol10", &m_mu2Nuisol10, "mu2_nuIsol10/F");
   dump->Branch("mu2_nuIsol20", &m_mu2Nuisol20, "mu2_nuIsol20/F");
   dump->Branch("mu2_nuIsol30", &m_mu2Nuisol30, "mu2_nuIsol30/F");
   dump->Branch("mu2_nuIsol40", &m_mu2Nuisol40, "mu2_nuIsol40/F");
   dump->Branch("mu2ProbChi2Fit", &m_mu2ProbChi2Fit, "mu2ProbChi2Fit/F");
   dump->Branch("mu2ProbChi2Match", &m_mu2ProbChi2Match, "mu2ProbChi2Match/F");
   dump->Branch("mu2PtRatio_ex", &m_mu2PtRatio_ex, "mu2PtRatio_ex/F");
   dump->Branch("mu2DeltaPtSign_ex", &m_mu2DeltaPtSign_ex, "mu2DeltaPtSign_ex/F");
   dump->Branch("mu2EnergyBalance", &m_mu2EnergyBalance, "mu2EnergyBalance/F");
   dump->Branch("mu2d0", &m_mu2d0, "mu2d0/F");
   dump->Branch("mu2z0", &m_mu2z0, "mu2z0/F");
   dump->Branch("mu2d0PVCorr", &m_mu2d0PVCorr, "mu2d0PVCorr/F");
   dump->Branch("mu2z0PVCorr", &m_mu2z0PVCorr, "mu2z0PVCorr/F");
   dump->Branch("mu2q", &m_mu2q, "mu2q/I");
   // trigger info vars
   dump->Branch("muL1MU0pass", &m_muL1MU0pass, "muL1MU0pass/I");
   dump->Branch("muL1MU6pass", &m_muL1MU6pass, "muL1MU6pass/I");
   dump->Branch("muL1MU10pass", &m_muL1MU10pass, "muL1MU10pass/I");
   dump->Branch("muL1MU0COMMpass", &m_muL1MU0COMMpass, "muL1MU0COMMpass/I");
   dump->Branch("muEFMU4pass", &m_muEFMU4pass, "muEFMU4pass/I");
   dump->Branch("muEFMU6pass", &m_muEFMU6pass, "muEFMU6pass/I");
   dump->Branch("muEFMU10pass", &m_muEFMU10pass, "muEFMU10pass/I");
   dump->Branch("muEFMU4MSONLYpass", &m_muEFMU4MSONLYpass, "muEFM4MSONLYpass/I");
  
   dump->Branch("mu1L1MU0dr", &m_mu1L1MU0dr, "mu1L1MU0dr/F");
   dump->Branch("mu1L1MU0index", &m_mu1L1MU0index, "mu1L1MU0index/I");
   dump->Branch("mu1L1MU6dr", &m_mu1L1MU6dr, "mu1L1MU6dr/F");
   dump->Branch("mu1L1MU6index", &m_mu1L1MU6index, "mu1L1MU6index/I");
   dump->Branch("mu1L1MU10dr", &m_mu1L1MU10dr, "mu1L1MU10dr/F");
   dump->Branch("mu1L1MU10index", &m_mu1L1MU10index, "mu1L1MU10index/I");
   dump->Branch("mu1L1MU0COMMdr", &m_mu1L1MU0COMMdr, "mu1L1MU0COMMdr/F");
   dump->Branch("mu1L1MU0COMMindex", &m_mu1L1MU0COMMindex, "mu1L1MU0COMMindex/I");

   dump->Branch("mu1L1MU0dr2", &m_mu1L1MU0dr2, "mu1L1MU0dr2/F");
   dump->Branch("mu1L1MU0index2", &m_mu1L1MU0index2, "mu1L1MU0index2/I");
   dump->Branch("mu1L1MU6dr2", &m_mu1L1MU6dr2, "mu1L1MU6dr2/F");
   dump->Branch("mu1L1MU6index2", &m_mu1L1MU6index2, "mu1L1MU6index2/I");
   dump->Branch("mu1L1MU10dr2", &m_mu1L1MU10dr2, "mu1L1MU10dr2/F");
   dump->Branch("mu1L1MU10index2", &m_mu1L1MU10index2, "mu1L1MU10index2/I");
   dump->Branch("mu1L1MU0COMMdr2", &m_mu1L1MU0COMMdr2, "mu1L1MU0COMMdr2/F");
   dump->Branch("mu1L1MU0COMMindex2", &m_mu1L1MU0COMMindex2, "mu1L1MU0COMMindex2/I");


   dump->Branch("mu1EFMU4dr", &m_mu1EFMU4dr, "mu1EFMU4dr/F");
   dump->Branch("mu1EFMU4index", &m_mu1EFMU4index, "mu1EFMU4index/I");
   dump->Branch("mu1EFMU6dr", &m_mu1EFMU6dr, "mu1EFMU6dr/F");
   dump->Branch("mu1EFMU6index", &m_mu1EFMU6index, "mu1EFMU6index/I");
   dump->Branch("mu1EFMU10dr", &m_mu1EFMU10dr, "mu1EFMU10dr/F");
   dump->Branch("mu1EFMU10index", &m_mu1EFMU10index, "mu1EFMU10index/I");
   dump->Branch("mu1EFMU4MSONLYdr", &m_mu1EFMU4MSONLYdr, "mu1EFMU4MSONLYdr/F");
   dump->Branch("mu1EFMU4MSONLYindex", &m_mu1EFMU4MSONLYindex, "mu1EFMU4MSONLYindex/I");
   
   dump->Branch("mu2L1MU0dr", &m_mu2L1MU0dr, "mu2L1MU0dr/F");
   dump->Branch("mu2L1MU0index", &m_mu2L1MU0index, "mu2L1MU0index/I");
   dump->Branch("mu2L1MU6dr", &m_mu2L1MU6dr, "mu2L1MU6dr/F");
   dump->Branch("mu2L1MU6index", &m_mu2L1MU6index, "mu2L1MU6index/I");
   dump->Branch("mu2L1MU10dr", &m_mu2L1MU10dr, "mu2L1MU10dr/F");
   dump->Branch("mu2L1MU10index", &m_mu2L1MU10index, "mu2L1MU10index/I");
   dump->Branch("mu2L1MU0COMMdr", &m_mu2L1MU0COMMdr, "mu2L1MU0COMMdr/F");
   dump->Branch("mu2L1MU0COMMindex", &m_mu2L1MU0COMMindex, "mu2L1MU0COMMindex/I");

   dump->Branch("mu2L1MU0dr2", &m_mu2L1MU0dr2, "mu2L1MU0dr2/F");
   dump->Branch("mu2L1MU0index2", &m_mu2L1MU0index2, "mu2L1MU0index2/I");
   dump->Branch("mu2L1MU6dr2", &m_mu2L1MU6dr2, "mu2L1MU6dr2/F");
   dump->Branch("mu2L1MU6index2", &m_mu2L1MU6index2, "mu2L1MU6index2/I");
   dump->Branch("mu2L1MU10dr2", &m_mu2L1MU10dr2, "mu2L1MU10dr2/F");
   dump->Branch("mu2L1MU10index2", &m_mu2L1MU10index2, "mu2L1MU10index2/I");
   dump->Branch("mu2L1MU0COMMdr2", &m_mu2L1MU0COMMdr2, "mu2L1MU0COMMdr2/F");
   dump->Branch("mu2L1MU0COMMindex2", &m_mu2L1MU0COMMindex2, "mu2L1MU0COMMindex2/I");

   dump->Branch("mu2EFMU4dr", &m_mu2EFMU4dr, "mu2EFMU4dr/F");
   dump->Branch("mu2EFMU4index", &m_mu2EFMU4index, "mu2EFMU4index/I");
   dump->Branch("mu2EFMU6dr", &m_mu2EFMU6dr, "mu2EFMU6dr/F");
   dump->Branch("mu2EFMU6index", &m_mu2EFMU6index, "mu2EFMU6index/I");
   dump->Branch("mu2EFMU10dr", &m_mu2EFMU10dr, "mu2EFMU10dr/F");
   dump->Branch("mu2EFMU10index", &m_mu2EFMU10index, "mu2EFMU10index/I");
   dump->Branch("mu2EFMU4MSONLYdr", &m_mu2EFMU4MSONLYdr, "mu2EFMU4MSONLYdr/F");
   dump->Branch("mu2EFMU4MSONLYindex", &m_mu2EFMU4MSONLYindex, "mu2EFMU4MSONLYindex/I");


  file << "  Done." << std::endl;
  
  // Here the definition of variables in SINLGE-MUON DUMPER TREE
  file << " Allocating the single  muons output tree " << std::endl;
  smu = new TTree("SingleMuonTree", "SingleMuon flat tree");
  smu->Branch("runNumber", &m_runNumber, "runNumber/I");
  smu->Branch("eventNumber ", &m_eventNumber, "eventNumber/I");
  smu->Branch("lumiBlock", &m_lumiBlock, "lumiBlock/I");
  smu->Branch("bcId", &m_bcId, "bcId/I");
  smu->Branch("smuHasIDTrk", &m_smuHasIDTrk, "smuHasIDTrk/I");
  smu->Branch("smuID_NSCTHits", &m_smuID_NSCTHits, "smuID_NSCTHits/I");
  smu->Branch("smuID_NPixHits", &m_smuID_NPixHits, "smuID_NPixHits/I");
  smu->Branch("smuID_NTRTHits", &m_smuID_NTRTHits, "smuID_NTRTHits/I");
  smu->Branch("smuq", &m_smuq, "smuq/I");
  smu->Branch("smuauthor", &m_smuauthor, "smuauthor/I");
  smu->Branch("smud0", &m_smud0, "smud0/F");
  smu->Branch("smuz0", &m_smuz0, "smuz0/F");
  smu->Branch("smueta", &m_smueta, "smueta/F");
  smu->Branch("smuphi", &m_smuphi, "smuphi/F");
  smu->Branch("smupt", &m_smupt, "smpupt/F");
  smu->Branch("smuchi2match", &m_smuchi2match, "smuchi2match/F");
  smu->Branch("smuchi2ID", &m_smuchi2ID, "smuchi2ID/F");
  smu->Branch("smumatchNDoF", &m_smumatchNDoF, "smumatchNDoF/F");
  smu->Branch("smufitchi2", &m_smufitchi2, "smufitchi2/F");
  smu->Branch("smufitNDoF", &m_smufitNDoF, "smufitNDoF/F");
  smu->Branch("smuetas", &m_smuetas, "smuetas/F");
  smu->Branch("smuphis", &m_smuphis, "smuphis/F");
  smu->Branch("smupts", &m_smupts, "smpupts/F");
  smu->Branch("smuetae", &m_smuetae, "smuetae/F");
  smu->Branch("smuphie", &m_smuphie, "smuphie/F");
  smu->Branch("smupte", &m_smupte, "smpupte/F");
  smu->Branch("smuetaID", &m_smuetaID, "smuetaID/F");
  smu->Branch("smuphiID", &m_smuphiID, "smuphiID/F");
  smu->Branch("smuptID", &m_smuptID, "smuptID/F");
  
  smu->Branch("smuIDTrackChi2Reduct", &m_smuIDTrackChi2Reduct, "smuIDTrackChi2Reduct/F");
  smu->Branch("smuIDProbChi2Fit", &m_smuIDProbChi2Fit, "smuIDProbChi2Fit/F");
  smu->Branch("smuIDTrackDeltaChi2", &m_smuIDTrackDeltaChi2, "smuIDTrackDeltaChi2/F");
  
  smu->Branch("smuetaex", &m_smuetaex, "smuetaex/F");
  smu->Branch("smuphiex", &m_smuphiex, "smuphiex/F");
  smu->Branch("smuptex", &m_smuptex, "smuptex/F");
  smu->Branch("smunMDT", &m_smunMDT, "smunMDT/I");
  smu->Branch("smunRPCeta", &m_smunRPCeta, "smunRPCeta/I");
  smu->Branch("smunRPCphi", &m_smunRPCphi, "smunRPCphi/I");
  smu->Branch("smunTGCeta", &m_smunTGCeta, "smunTGCeta/I");
  smu->Branch("smunTGCphi", &m_smunTGCphi, "smunTGCphi/I");
  smu->Branch("smunCSCeta", &m_smunCSCeta, "smunCSCeta/I");
  smu->Branch("smunCSCphi", &m_smunCSCphi, "smunCSCphi/I");
  smu->Branch("smuEtisol10", &m_smuEtisol10, "smuEtisol10/F");
  smu->Branch("smuEtisol20", &m_smuEtisol20, "smuEtisol20/F");
  smu->Branch("smuEtisol30", &m_smuEtisol30, "smuEtisol30/F");
  smu->Branch("smuEtisol40", &m_smuEtisol40, "smuEtisol40/F");  
  smu->Branch("smuPtisol10", &m_smuPtisol10, "smuPtisol10/F");
  smu->Branch("smuPtisol20", &m_smuPtisol20, "smuPtisol20/F");
  smu->Branch("smuPtisol30", &m_smuPtisol30, "smuPtisol30/F");
  smu->Branch("smuPtisol40", &m_smuPtisol40, "smuPtisol40/F");  
  smu->Branch("smuNuisol10", &m_smuNuisol10, "smuNuisol10/F");
  smu->Branch("smuNuisol20", &m_smuNuisol20, "smuNuisol20/F");
  smu->Branch("smuNuisol30", &m_smuNuisol30, "smuNuisol30/F");
  smu->Branch("smuNuisol40", &m_smuNuisol40, "smuNuisol40/F");  
  smu->Branch("smud0PVCorr", &m_smud0PVCorr, "smud0PVCorr/F");

  smu->Branch("smuL1MU0pass", &m_smuL1MU0pass, "smuL1MU0pass/I");
  smu->Branch("smuL1MU6pass", &m_smuL1MU6pass, "smuL1MU6pass/I");
  smu->Branch("smuL1MU10pass", &m_smuL1MU10pass, "smuL1MU10pass/I");
  smu->Branch("smuL1MU0COMMpass", &m_smuL1MU0COMMpass, "smuL1MU0COMMpass/I");
  smu->Branch("smuEFMU4pass", &m_smuEFMU4pass, "smuEFMU4pass/I");
  smu->Branch("smuEFMU6pass", &m_smuEFMU6pass, "smuEFMU6pass/I");
  smu->Branch("smuEFMU10pass", &m_smuEFMU10pass, "smuEFMU10pass/I");
  smu->Branch("smuEFMU4MSONLYpass", &m_smuEFMU4MSONLYpass, "smuEFM4MSONLYpass/I");
  
  
  smu->Branch("smuL1MU0dr", &m_smuL1MU0dr, "smuL1MU0dr/F");
  smu->Branch("smuL1MU0index", &m_smuL1MU0index, "smuL1MU0index/I");
  smu->Branch("smuL1MU6dr", &m_smuL1MU6dr, "smuL1MU6dr/F");
  smu->Branch("smuL1MU6index", &m_smuL1MU6index, "smuL1MU6index/I");
  smu->Branch("smuL1MU10dr", &m_smuL1MU10dr, "smuL1MU10dr/F");
  smu->Branch("smuL1MU10index", &m_smuL1MU10index, "smuL1MU10index/I");
  smu->Branch("smuL1MU0COMMdr", &m_smuL1MU0COMMdr, "smuL1MU0COMMdr/F");
  smu->Branch("smuL1MU0COMMindex", &m_smuL1MU0COMMindex, "smuL1MU0COMMindex/I");
  smu->Branch("smuEFMU4dr", &m_smuEFMU4dr, "smuEFMU4dr/F");
  smu->Branch("smuEFMU4index", &m_smuEFMU4index, "smuEFMU4index/I");
  smu->Branch("smuEFMU6dr", &m_smuEFMU6dr, "smuEFMU6dr/F");
  smu->Branch("smuEFMU6index", &m_smuEFMU6index, "smuEFMU6index/I");
  smu->Branch("smuEFMU10dr", &m_smuEFMU10dr, "smuEFMU10dr/F");
  smu->Branch("smuEFMU10index", &m_smuEFMU10index, "smuEFMU10index/I");
  smu->Branch("smuEFMU4MSONLYdr", &m_smuEFMU4MSONLYdr, "smuEFMU4MSONLYdr/F");
  smu->Branch("smuEFMU4MSONLYindex", &m_smuEFMU4MSONLYindex, "smuEFMU4MSONLYindex/I");

  smu->Branch("smuProbChi2Fit", &m_smuProbChi2Fit, "smuProbChi2Fit/F");
  smu->Branch("smuProbChi2Match", &m_smuProbChi2Match, "smuProbChi2Match/F");
  smu->Branch("smuPtRatio_ex", &m_smuPtRatio_ex, "smuPtRatio_ex/F");
  smu->Branch("smuDeltaPtSign_ex", &m_smuDeltaPtSign_ex, "smuDeltaPtSign_ex/F");
  smu->Branch("smuEnergyBalance", &m_smuEnergyBalance, "smuEnergyBalance/F");
  file << "  Done." << std::endl;

  // Here the definition of variables in MULTI-MUON DUMPER TREE
  file << " Allocating the single  muons output tree " << std::endl;
  multi = new TTree("MultiMuonTree", "MultiMuon flat tree");
  multi->Branch("runNumber", &m_runNumber, "runNumber/I");
  multi->Branch("eventNumber ", &m_eventNumber, "eventNumber/I");
  multi->Branch("lumiBlock", &m_lumiBlock, "lumiBlock/I");
  multi->Branch("bcId", &m_bcId, "bcId/I");
  multi->Branch("mnmuons", &m_mnmuons, "mnmuons/I");
  multi->Branch("mchargetot", &m_mchargetot, "mchargetot/I");
  multi->Branch("mIsStaco", &m_mIsStaco, "mIsStaco/I");
  multi->Branch("mmupairtype", &m_mmupairtype, "mmupairtype/I");
  multi->Branch("mmu1eta", &m_mmu1eta, "mmu1eta/F");
  multi->Branch("mmu1phi", &m_mmu1phi, "mmu1phi/F");
  multi->Branch("mmu1pt",  &m_mmu1pt,  "mmu1pt/F");
  multi->Branch("mmu1author",&m_mmu1author, "mmu1author/I");
  multi->Branch("mmu1vertex",&m_mmu1vertex, "mmu1vertex/I");
  multi->Branch("mmu1vxtype",&m_mmu1vxtype, "mmu1vxtype/I");
  multi->Branch("mmu2eta", &m_mmu2eta, "mmu2eta/F");
  multi->Branch("mmu2phi", &m_mmu2phi, "mmu2phi/F");
  multi->Branch("mmu2pt",  &m_mmu2pt,  "mmu2pt/F");
  multi->Branch("mmu2author",&m_mmu2author, "mmu2author/I");
  multi->Branch("mmu2vertex",&m_mmu2vertex, "mmu2vertex/I");
  multi->Branch("mmu2vxtype",&m_mmu2vxtype, "mmu2vxtype/I");
  multi->Branch("mmu3eta", &m_mmu3eta, "mmu3eta/F");
  multi->Branch("mmu3phi", &m_mmu3phi, "mmu3phi/F");
  multi->Branch("mmu3pt",  &m_mmu3pt,  "mmu3pt/F");
  multi->Branch("mmu3author",&m_mmu3author, "mmu3author/I");
  multi->Branch("mmu3vertex",&m_mmu3vertex, "mmu3vertex/I");
  multi->Branch("mmu3vxtype",&m_mmu3vxtype, "mmu3vxtype/I");
  multi->Branch("mmu4eta", &m_mmu4eta, "mmu4eta/F");
  multi->Branch("mmu4phi", &m_mmu4phi, "mmu4phi/F");
  multi->Branch("mmu4pt",  &m_mmu4pt,  "mmu4pt/F");
  multi->Branch("mmu4author",&m_mmu4author, "mmu4author/I");
  multi->Branch("mmu4vertex",&m_mmu4vertex, "mmu4vertex/I");
  multi->Branch("mmu4vxtype",&m_mmu4vxtype, "mmu4vxtype/I");  
  multi->Branch("mmotherAmass",  &m_mmotherAmass,  "mmotherAmass/F");
  multi->Branch("mmotherAp",  &m_mmotherAp,  "mmotherAp/F");
  multi->Branch("mmotherAcharge",  &m_mmotherAcharge,  "mmotherAcharge/I");
  multi->Branch("mmotherApt",  &m_mmotherApt,  "mmotherApt/F");
  multi->Branch("mmotherAeta",  &m_mmotherAeta,  "mmotherAeta/F");
  multi->Branch("mmotherAphi",  &m_mmotherAphi,  "mmotherAphi/F");
  multi->Branch("mmotherBmass",  &m_mmotherBmass,  "mmotherBmass/F");
  multi->Branch("mmotherBp",  &m_mmotherBp,  "mmotherBp/F");
  multi->Branch("mmotherBcharge",  &m_mmotherBcharge,  "mmotherBcharge/I");
  multi->Branch("mmotherBpt",  &m_mmotherBpt,  "mmotherBpt/F");
  multi->Branch("mmotherBeta",  &m_mmotherBeta,  "mmotherBeta/F");
  multi->Branch("mmotherBphi",  &m_mmotherBphi,  "mmotherBphi/F");
  multi->Branch("mgmothermass",  &m_mgmothermass,  "mgmothermass/F");
  multi->Branch("mgmotherp",  &m_mgmotherp,  "mgmotherp/F");
  multi->Branch("mgmothercharge",  &m_mgmothercharge,  "mgmothercharge/I");
  multi->Branch("mgmotherpt",  &m_mgmotherpt,  "mgmotherpt/F");
  multi->Branch("mgmothereta",  &m_mgmothereta,  "mgmothereta/F");
  multi->Branch("mgmotherphi",  &m_mgmotherphi,  "mgmotherphi/F");
   if (fChain == 0) return;

   Long64_t nentries = fChain->GetEntries();

   cout<< " tot nev=" <<nentries<< endl;

   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=0; jentry<nentries;jentry++) 
     {

       //       if(jentry>3000000) break;
       Long64_t ientry = LoadTree(jentry);
       if (ientry < 0) break;
       //nb = fChain->GetEntry(jentry);   nbytes += nb;

       if( jentry%1000==0) cout<< " evt=" << jentry<<endl; 
       file<< " jentry=" << jentry<<endl;

       
     //   // Get needed branches from JpsiIt ntuple
        b_evinfo_runNumber->GetEntry(ientry);
        b_evinfo_evtNumber->GetEntry(ientry);
        b_evinfo_lumiBlock->GetEntry(ientry);
        b_evinfo_bcId->GetEntry(ientry);
        b_VtxType->GetEntry(ientry);
        // ID
        //b_ID_TrackChi2Reduct->GetEntry(ientry);
	//b_ID_ProbChi2Fit->GetEntry(ientry);  
	//b_ID_TrackDeltaChi2->GetEntry(ientry); 
	b_ID_nTracks->GetEntry(ientry);
        b_ID_VxIndex->GetEntry(ientry);
        b_ID_nPixelHits->GetEntry(ientry);
        b_ID_nSCTHits->GetEntry(ientry);
        b_ID_nTRTHits->GetEntry(ientry);
        b_ID_TrackPhi->GetEntry(ientry);
        b_ID_TrackTheta->GetEntry(ientry);
        b_ID_TrackPt->GetEntry(ientry);
        // muons
        b_mu_n->GetEntry(ientry);
        b_mu_IDindex->GetEntry(ientry);
        b_mu_author->GetEntry(ientry);
        b_mu_Q->GetEntry(ientry);
        b_mu_author->GetEntry(ientry);
        b_mu_Eta->GetEntry(ientry);
        b_mu_Phi->GetEntry(ientry);
        b_mu_Pt->GetEntry(ientry);
        b_mu_matchChi2->GetEntry(ientry);
        b_mu_matchNDoF->GetEntry(ientry);
        b_mu_fitChi2->GetEntry(ientry);	 
        b_mu_fitNDoF->GetEntry(ientry); 
        //b_mu_ProbChi2Fit->GetEntry(ientry);
	//b_mu_ProbChi2Match->GetEntry(ientry); 
	//b_mu_PtRatio_ex->GetEntry(ientry);
	//b_mu_DeltaPtSign_ex->GetEntry(ientry);
	//b_mu_EnergyBalance->GetEntry(ientry); 

	b_mu_Eta_ms->GetEntry(ientry);       	   
        b_mu_Phi_ms->GetEntry(ientry);  		   
        b_mu_Pt_ms->GetEntry(ientry); 		   
        b_mu_IDEtaExtrapolatedToMS->GetEntry(ientry); 
        b_mu_IDPhiExtrapolatedToMS->GetEntry(ientry); 
        b_mu_IDPtExtrapolatedToMS->GetEntry(ientry); 
        b_mu_Eta_ex->GetEntry(ientry);        
        b_mu_Phi_ex->GetEntry(ientry);        
        b_mu_Pt_ex->GetEntry(ientry); 	   
        b_mu_nMDTHits->GetEntry(ientry); 
        b_mu_nRPCEtaHits->GetEntry(ientry);  
        b_mu_nRPCPhiHits->GetEntry(ientry);  
        b_mu_nTGCEtaHits->GetEntry(ientry);  
        b_mu_nTGCPhiHits->GetEntry(ientry);  
        b_mu_nCSCEtaHits->GetEntry(ientry);  
        b_mu_nCSCPhiHits->GetEntry(ientry);  
        b_mu_etIsol10->GetEntry(ientry); 
        b_mu_etIsol20->GetEntry(ientry); 
        b_mu_etIsol30->GetEntry(ientry); 
        b_mu_etIsol40->GetEntry(ientry); 
        b_mu_ptIsol10->GetEntry(ientry); 
        b_mu_ptIsol20->GetEntry(ientry); 
        b_mu_ptIsol30->GetEntry(ientry); 
        b_mu_ptIsol40->GetEntry(ientry); 
        b_mu_nuIsol10->GetEntry(ientry); 
        b_mu_nuIsol20->GetEntry(ientry); 
        b_mu_nuIsol30->GetEntry(ientry); 
        b_mu_nuIsol40->GetEntry(ientry); 
        b_mu_d0->GetEntry(ientry); 
        b_mu_z0->GetEntry(ientry); 
        b_mu_d0corrPV->GetEntry(ientry); 
        b_mu_z0corrPV->GetEntry(ientry); 
        // Trigger 
        b_L1_ispassed->GetEntry(ientry);
        b_L2_ispassed->GetEntry(ientry);
        b_EF_ispassed->GetEntry(ientry);
        b_L1_trigObj_names->GetEntry(ientry);
        b_L1_trigObj_eta->GetEntry(ientry);
        b_L1_trigObj_phi->GetEntry(ientry);
        b_L1_trigObj_pt->GetEntry(ientry);
        b_L1_trigObj_type->GetEntry(ientry);
        b_L1_trigObj_id->GetEntry(ientry);
        b_L1_trigObj_parentId->GetEntry(ientry);  
        b_EF_trigObj_names->GetEntry(ientry);
        b_EF_trigObj_eta->GetEntry(ientry);  
        b_EF_trigObj_phi->GetEntry(ientry); 
        b_EF_trigObj_pt->GetEntry(ientry);   
        b_EF_trigObj_type->GetEntry(ientry); 
        b_EF_trigObj_id->GetEntry(ientry);   
        b_EF_trigObj_parentId->GetEntry(ientry);  
        // Dimuons
        b_dimuon_n->GetEntry(ientry);   //!		  
        b_dimuon_invMass->GetEntry(ientry);   //!	  
        b_dimuon_pt->GetEntry(ientry);   //!		  
        b_dimuon_mu1Author->GetEntry(ientry);   //!	  
        b_dimuon_mu1Index->GetEntry(ientry);   //!	  
        b_dimuon_mu2Author->GetEntry(ientry);   //!	  
        b_dimuon_mu2Index->GetEntry(ientry);   //!	  
        b_dimuon_corrPrimVtxIndex->GetEntry(ientry);   //
        b_dimuon_corrPrimVtxPosX->GetEntry(ientry);   //!
        b_dimuon_corrPrimVtxPosY->GetEntry(ientry);   //!
        b_dimuon_corrPrimVtxPosZ->GetEntry(ientry);   //!
        b_dimuon_mu1CorrPx->GetEntry(ientry);   //!	  
        b_dimuon_mu1CorrPy->GetEntry(ientry);   //!	  
        b_dimuon_mu1CorrPz->GetEntry(ientry);   //!	  
        b_dimuon_mu2CorrPx->GetEntry(ientry);   //!	  
        b_dimuon_mu2CorrPy->GetEntry(ientry);   //!	  
        b_dimuon_mu2CorrPz->GetEntry(ientry);   //!	  
        b_dimuon_IDEtaExtrapolatedToMS->GetEntry(ientry);
        b_dimuon_IDPhiExtrapolatedToMS->GetEntry(ientry);
        b_dimuon_IDPtExtrapolatedToMS->GetEntry(ientry); 
        b_dimuon_mu1CorrPx->GetEntry(ientry); 
        b_dimuon_mu1CorrPy->GetEntry(ientry);
        b_dimuon_mu1CorrPz->GetEntry(ientry);
        b_dimuon_mu2CorrPx->GetEntry(ientry);
        b_dimuon_mu2CorrPy->GetEntry(ientry);
        b_dimuon_mu2CorrPz->GetEntry(ientry);

        b_muctpi_n->GetEntry(ientry);   //!
        b_muctpi_sector->GetEntry(ientry);   //!
        b_muctpi_sysId->GetEntry(ientry);   //!
        b_muctpi_pad->GetEntry(ientry);   //!
        b_muctpi_roiNumber->GetEntry(ientry);   //!
        b_muctpi_eta->GetEntry(ientry);   //!
        b_muctpi_phi->GetEntry(ientry);   //!
     
       /////
       m_runNumber   = evinfo_runNumber;
       m_eventNumber = evinfo_evtNumber; 
       m_lumiBlock   = evinfo_lumiBlock;
       m_bcId        = evinfo_bcId;


       //
       /// Here start to fill SINGLE-MUON TREE 
       //

       if (doSMuon)
	 {
	   //file <<endl << "** NEW EVENT WITH Number of MUON = "<<mu_n<<endl;
	   bool baseline_sel = InclusiveBaselineSelection(); // Ask for baseline inclusive selection
	   if (! baseline_sel ) continue;
	   //if (baseline_sel ) file <<" baseline is ok"<<endl;
	   for (int ind_m=0; ind_m < mu_n ; ind_m ++)
	     {

	       //file << endl << " Iteration on muon number = " << ind_m<<endl;
	       //file << " Author is "<< mu_author->at(ind_m) << endl;

	       int mu_id= mu_IDindex->at(ind_m); // index of ID track associated to muon 
	       
	       m_smuID_NSCTHits=-1;
	       m_smuID_NPixHits=-1;
	       m_smuID_NTRTHits=-1; 
	       m_smuchi2IDnew=-1;
	       m_smuHasIDTrk=1;
	       if (mu_id <0) m_smuHasIDTrk=-1; // if muon hasn't any associated ID track
	       if (m_smuHasIDTrk==1)
		 {
		   m_smuID_NSCTHits = ID_nSCTHits->at(mu_id);
		   m_smuID_NPixHits = ID_nPixelHits->at(mu_id);
		   m_smuID_NTRTHits = ID_nTRTHits->at(mu_id); 
		 }
	       
	       // CUT: Muon Author selection
	       if (mu_author->at(ind_m) != 6 && mu_author->at(ind_m) != 12  
		   && mu_author->at(ind_m) != 7  && mu_author->at(ind_m) != 13  
		   && mu_author->at(ind_m) != 18) continue;

	       m_smuq         = mu_Q->at(ind_m);      
	       m_smuauthor    = mu_author->at(ind_m); 
	       m_smueta       = mu_Eta->at(ind_m);    
	       m_smuphi       = mu_Phi->at(ind_m);    
	       m_smupt        = mu_Pt->at(ind_m);     
	       m_smuchi2match = mu_matchChi2->at(ind_m); 
	       m_smumatchNDoF = mu_matchNDoF->at(ind_m); 
	       m_smufitchi2   = mu_fitChi2->at(ind_m); 	 
	       m_smufitNDoF   = mu_fitNDoF->at(ind_m);   
	       //if (m_smuHasIDTrk==1) m_smuchi2ID  = ID_TrackChi2->at(mu_id);
	       
	       m_smuetas     = mu_Eta_ms->at(ind_m);      	   
	       m_smuphis     = mu_Phi_ms->at(ind_m); 		   
	       m_smupts      = mu_Pt_ms->at(ind_m); 		   
	       m_smuetae     = mu_IDEtaExtrapolatedToMS->at(ind_m);  // extrapolator vars    
	       m_smuphie     = mu_IDPhiExtrapolatedToMS->at(ind_m);      
	       m_smupte      = mu_IDPtExtrapolatedToMS->at(ind_m); 

	       if (m_smuHasIDTrk==1) // ID tracks vars
		  {
		    float eta_id_trk = -log(fabs(tan(ID_TrackTheta->at(mu_id)))/2.);
		    m_smuetaID    = eta_id_trk;
		    m_smuphiID    = ID_TrackPhi->at(mu_id);
		    m_smuptID     = ID_TrackPt->at(mu_id);
		    //Pi and K rejection
		    m_smuIDTrackChi2Reduct=ID_TrackChi2Reduct->at(mu_id);
		    m_smuIDProbChi2Fit=ID_ProbChi2Fit->at(mu_id);
		    m_smuIDTrackDeltaChi2=ID_TrackDeltaChi2->at(mu_id);
   

		  }

	       m_smuetaex    = mu_Eta_ex->at(ind_m);       
	       m_smuphiex    = mu_Phi_ex->at(ind_m);       
	       m_smuptex     = mu_Pt_ex->at(ind_m); 	   
	       m_smunMDT     = mu_nMDTHits->at(ind_m);        
	       m_smunRPCeta  = mu_nRPCEtaHits->at(ind_m);  
	       m_smunRPCphi  = mu_nRPCPhiHits->at(ind_m);  
	       m_smunTGCeta  = mu_nTGCEtaHits->at(ind_m);  
	       m_smunTGCphi  = mu_nTGCPhiHits->at(ind_m);  
	       m_smunCSCeta  = mu_nCSCEtaHits->at(ind_m);  
	       m_smunCSCphi  = mu_nCSCPhiHits->at(ind_m);  
	       
	       // isolation :
	       m_smuEtisol10 = mu_etIsol10->at(ind_m); 
	       m_smuEtisol20 = mu_etIsol20->at(ind_m);
	       m_smuEtisol30 = mu_etIsol30->at(ind_m);
	       m_smuEtisol40 = mu_etIsol40->at(ind_m);
	       m_smuPtisol10 = mu_ptIsol10->at(ind_m); 
	       m_smuPtisol20 = mu_ptIsol20->at(ind_m);
	       m_smuPtisol30 = mu_ptIsol30->at(ind_m);
	       m_smuPtisol40 = mu_ptIsol40->at(ind_m);
	       m_smuNuisol10 = mu_nuIsol10->at(ind_m); 
	       m_smuNuisol20 = mu_nuIsol20->at(ind_m);
	       m_smuNuisol30 = mu_nuIsol30->at(ind_m);
	       m_smuNuisol40 = mu_nuIsol40->at(ind_m);

	       // vertexing :
	       m_smud0       = mu_d0->at(ind_m);       
	       m_smuz0       = mu_z0->at(ind_m);       
	       m_smud0PVCorr = mu_d0corrPV->at(ind_m); 
	       m_smuz0PVCorr = mu_z0corrPV->at(ind_m); 
	       
	       // Trigger: Start here

	       m_smuL1MU0pass=0;
	       m_smuL1MU6pass=0;
	       m_smuL1MU10pass=0;
	       m_smuL1MU0COMMpass=0;
	       m_smuEFMU4pass=0;	   
	       m_smuEFMU6pass=0;	   
	       m_smuEFMU10pass=0;	   
	       m_smuEFMU4MSONLYpass=0;
	       m_smuL1MU0dr=0;	       
	       m_smuL1MU0index=0;      
	       m_smuL1MU6dr=0;	       
	       m_smuL1MU6index=0;      
	       m_smuL1MU10dr=0;	       
	       m_smuL1MU10index=0;     
	       m_smuL1MU0COMMdr=0;     
	       m_smuL1MU0COMMindex=0;  
	       m_smuEFMU4dr=0;	       
	       m_smuEFMU4index=0;      
	       m_smuEFMU6dr=0;	       
	       m_smuEFMU6index=0;      
	       m_smuEFMU10dr=0;	       
	       m_smuEFMU10index=0;     
	       m_smuEFMU4MSONLYdr=0;   
	       m_smuEFMU4MSONLYindex=0;
	       
	       //Pi and K rejection
	       /*
	       m_smuProbChi2Fit=mu_ProbChi2Fit->at(ind_m);
	       m_smuProbChi2Match=mu_ProbChi2Match->at(ind_m);
	       m_smuPtRatio_ex=mu_PtRatio_ex->at(ind_m);
	       m_smuDeltaPtSign_ex=mu_DeltaPtSign_ex->at(ind_m);
	       m_smuEnergyBalance=mu_EnergyBalance->at(ind_m);
	       */

	        L1_items["L1_MU0"]=0;
		L1_items["L1_MU6"]=0;
		L1_items["L1_MU10"]=0;
		L1_items["L1_MU0_COMM"]=0;
		
		EF_items["EF_mu4"]=0;
		EF_items["EF_mu6"]=0;
		EF_items["EF_mu10"]=0;
		EF_items["EF_mu4_MSonly_MB2_noL2_EFFS"]=0;

		 

		if (L1_ispassed->size() >0 )
		   {
		     mapType::const_iterator end = L1_items.end(); 
		     for(mapType::const_iterator it = L1_items.begin(); it != end; ++it)
		       {
			 //file << " requested item is "<<it->first<<endl;
			 int passed=0;
			 bool L1_item_isfound=0;
			 while (passed < L1_ispassed->size() &&  !L1_item_isfound)
			   {
			     //file << " Item is : "<< L1_ispassed->at(passed)<<endl;
			     if (L1_ispassed->at(passed) == it->first )
			     { 
			       //file << "Chain " << it->first << "is OK" <<endl;
			       L1_item_isfound = 1 ;
			       string t = it->first;
			       L1_items[t]=1;
			     } 
			     passed++;
			   }   
		       }
		   } // end of LVl1


	

	       if (EF_ispassed->size() >0 )
		   {
		     //file <<"EF is passed"<<endl;
		     //file <<"EF item trovati hanno size = "<<EF_ispassed->size()<<endl;
		     mapType::const_iterator end = EF_items.end(); 
		     for(mapType::const_iterator it = EF_items.begin(); it != end; ++it)
		       {
			 //file << " requested EF item is "<<it->first<<endl;
			 int passed=0;
			 bool EF_item_isfound=0;
			 while (passed < EF_ispassed->size() &&  !EF_item_isfound)
			   {
			     //file << " Item is : "<< EF_ispassed->at(passed)<<endl;
			     if (EF_ispassed->at(passed) == it->first )
			     { 
			       //file << "found rquested item "<<endl;
			       EF_item_isfound = 1 ;
			       string t = it->first;
			       EF_items[t]=1;
			     } 
			     passed++;
			   }   
		       }
		   }
	       
	       //
	       // Trigger: matching
	       //
	       int index_trigger_chain=1000;
	       int index_trigger_chain2=1000;
	       bool  IsNotCombined = 0;
	       float L1_DR=1000;
	       float L1_DR2=1000;
	       float EF_DR=1000;
	       if (mu_author->at(ind_m) != 12 || mu_author->at(ind_m) != 6) IsNotCombined = 1; 

	       ////// L1MU0
	       mapType::const_iterator it = L1_items.find("L1_MU0");
	       string chain = it->first;

	       if (it->second == 1) 
		 {
		   m_smuL1MU0pass = 1;  //L1MU0 is fired in this event
		   MatchL1Chain(ind_m, L1_DR, L1_DR2, IsNotCombined, index_trigger_chain, index_trigger_chain2, chain);
		   m_smuL1MU0dr    =  L1_DR;
		   m_smuL1MU0index =  index_trigger_chain;
		 }

	       ////// L1MU6
	       it = L1_items.find("L1_MU6");
	       chain = it->first;
	       L1_DR=1000;
	       index_trigger_chain=1000;

	       if (it->second == 1) 
		 {
		   m_smuL1MU6pass = 1;  //L1MU6 is fired in this event
		   MatchL1Chain(ind_m, L1_DR, L1_DR2, IsNotCombined, index_trigger_chain, index_trigger_chain2, chain);
		   m_smuL1MU6dr    =  L1_DR;
		   m_smuL1MU6index =  index_trigger_chain;
		 }

	       ////// L1MU10
	       it = L1_items.find("L1_MU10");
	       chain = it->first;
	       L1_DR=1000;
	       index_trigger_chain=1000;
	       if (it->second == 1) 
		 {
		   m_smuL1MU10pass = 1;  //L1MU10 is fired in this event
		   MatchL1Chain(ind_m, L1_DR, L1_DR2, IsNotCombined, index_trigger_chain, index_trigger_chain2, chain);
		   m_smuL1MU10dr    =  L1_DR;
		   m_smuL1MU10index =  index_trigger_chain;
		 }
	       
	       ////// L1MU0COMM
	      it = L1_items.find("L1_MU0_COMM");
	      chain = it->first;
	      L1_DR=1000;
	      index_trigger_chain=1000;
	      
	      if (it->second == 1) 
		 {
		   m_smuL1MU0COMMpass = 1;  //L1MU0COMM is fired in this event
		   MatchL1Chain(ind_m, L1_DR, L1_DR2, IsNotCombined, index_trigger_chain, index_trigger_chain2, chain);
		   m_smuL1MU0COMMdr    =  L1_DR;
		   m_smuL1MU0COMMindex =  index_trigger_chain;
		 }

	       
	      ////// EFMU4
	      mapType::const_iterator cc = EF_items.find("EF_mu4");
	      chain = cc->first;
	      L1_DR=1000;
	      index_trigger_chain=1000;
	       
	       if (cc->second == 1) 
		 {
		   m_smuEFMU4pass = 1;  //EFMU4 is fired in this event
		   MatchEFChain(ind_m, EF_DR, IsNotCombined, index_trigger_chain, chain);
		   m_smuEFMU4dr    =  EF_DR;
		   m_smuEFMU4index =  index_trigger_chain;
		 }
	      
	        ////// EFMU6
	       cc = EF_items.find("EF_mu6");
	       chain = cc->first;
	       L1_DR=1000;
	       index_trigger_chain=1000;
	       if (cc->second == 1) 
		 {
		   m_smuEFMU6pass = 1;  //EFMU6 is fired in this event
		   MatchEFChain(ind_m, EF_DR, IsNotCombined, index_trigger_chain, chain);
		   m_smuEFMU6dr    =  EF_DR;
		   m_smuEFMU6index =  index_trigger_chain;
		 }
	       
	        ////// EFMU10
	       cc = EF_items.find("EF_mu10");
	       chain = cc->first;
	       L1_DR=1000;
	       index_trigger_chain=1000;

	       if (cc->second == 1) 
		 {
		   m_smuEFMU10pass = 1;  //EFMU10 is fired in this event
		   MatchEFChain(ind_m, EF_DR, IsNotCombined, index_trigger_chain, chain);
		   m_smuEFMU10dr    =  EF_DR;
		   m_smuEFMU10index =  index_trigger_chain;
		 }
	       
	        ////// EFMU4MSONLY
	      cc = EF_items.find("EF_mu4_MSonly_MB2_noL2_EFFS");
	      chain = cc->first;
	      L1_DR=1000;
	      index_trigger_chain=1000;
	      
	       if (cc->second == 1) 
		 {
		   m_smuEFMU4MSONLYpass = 1;  //EFMU4MSONLY is fired in this event
		   MatchEFChain(ind_m, EF_DR, IsNotCombined, index_trigger_chain, chain);
		   m_smuEFMU4MSONLYdr    =  EF_DR;
		   m_smuEFMU4MSONLYindex =  index_trigger_chain;
		 }
	       
	       // end of Trigger: matching
	       
	       smu->Fill();
	     }
	 }
       
       // end of filling for SINGLE-MUON TREE 
       
       //
       /// Here start to fill MULTIMUON TREE 
       //
       
       if (doMultiMuon && mu_n >= 4)
	 {
	   // here reset all vars in MultiMuons Tree
	   m_mnmuons = 0;
	   m_mchargetot= 0;
	   m_mIsStaco = 0;
	   //
	   file << " ** Enter in Multimuon **" << endl;
	  
	   bool baseline_sel = InclusiveBaselineSelection(); // Ask for baseline inclusive 
	                                 
	   if (! baseline_sel ) continue;  // selection on primary vertex
	   file << " Pass BaselineInclusiveSelection" << endl;

	   vector < my_pair > staco_list;
	   vector < my_pair > muid_list;
	   
	   for (int im= 0; im<mu_n; im++)
	     {   
	       Bool_t good_mu = false;
	       float  mu_p ;
	       int    muid = mu_IDindex->at(im);
	       if ( muid < 0 ) continue;
	       if ( ID_nSCTHits->at(muid) >= 6 && ID_nPixelHits->at(muid) >= 1) good_mu = true; 
	       
	       file << " Muon n." << im << "pass quality selection on Pixel and SCT" << endl;
	       
	       // get p for this muon
	       //mu_p = ID_TrackPt->at(muid)/sin(ID_TrackTheta->at(muid));
	       mu_p = mu_Pt->at(im)/sin(2*atan(exp(-1*mu_Eta->at(im))));
	       // record this muon as a c++ pair with mu_index and p 
	       my_pair mu_pair (im,mu_p);
	       
	       if (mu_author->at(im) == 6  || mu_author->at(im) == 7 ) 	
		 staco_list.push_back(mu_pair);
	       if (mu_author->at(im) == 12 || mu_author->at(im) == 18)	
		 muid_list.push_back(mu_pair);
	       
	     } // end of for (int im = 0...)

	   if (staco_list.size() >= 4 || muid_list.size() >= 4)
	     {   
	       bool ReadStacoList = false;
	       bool ReadMuidList = false;
	       if (staco_list.size() >= 4) ReadStacoList = true;
	       if ( muid_list.size() >= 4) ReadMuidList = true;
	       
	       if (ReadStacoList)
		 { 
		   m_mnmuons=staco_list.size();  
		   int n_of_muons_to_be_erased = staco_list.size() - 4; 
		   // at max 4 muons are analyzed
		   // Ordering both vector of muons with respect the muons p
		   sort(staco_list.begin(), staco_list.end(), sort_pred);
		   if (n_of_muons_to_be_erased > 0)
		     staco_list.erase (staco_list.begin(),
				       staco_list.begin()+n_of_muons_to_be_erased);
		   
		   int staco_flag = 1;
		   
		   FillMulti(staco_list, staco_flag);
		 }

	       if (ReadMuidList)
		 {
		   m_mnmuons=muid_list.size();   
		   int n_of_muons_to_be_erased = muid_list.size() - 4; 
		   // at max 4 muons are analyzed
		   // Ordering both vector of muons with respect the muons p
		   sort(muid_list.begin(), muid_list.end(), sort_pred);
		   if (n_of_muons_to_be_erased > 0)
		     muid_list.erase (muid_list.begin(),
				       muid_list.begin()+n_of_muons_to_be_erased);
		   
		   int muid_flag = 0;
		  
		   FillMulti(muid_list, muid_flag);
		 }

	      
	       file << " Multimuons : Fill! **"<<endl;
	     }
	 }
       // end of filling for MULTI-MUON TREE 



       if (doDimuon)
	 {   
       //
       /// Here start to fill DIMUON TREE 
       //
	   
       if (IS_DEBUG) file <<endl<<" Run number is "<<evinfo_runNumber<<endl;
       if (IS_DEBUG) file <<endl<<" Event number is "<<evinfo_evtNumber<<endl;
       if (IS_DEBUG) file <<endl << "** NEW EVENT WITH NDIMUONS = "<<dimuon_n<<endl;
       bool baseline_sel = InclusiveBaselineSelection(); // Ask for baseline inclusive selection
       if (! baseline_sel ) continue;
       if (IS_DEBUG) file <<" Pass Baseline inclusive selection"<<endl;

       for (int i = 0; i < dimuon_n; i++) 
	 {

	   if (IS_DEBUG) file <<endl<<" Dimuon n." <<i<<endl;

	   m_invmass_novtx=0;
	   m_invmass=0;
	   m_q=0;
	   m_pt=0;
	   m_deltaEta=0;
	   m_deltaPhi=0;
	   m_deltaPt=0;
	   m_family=0;
	   m_etarange=0;
	   m_mu1index=0;
	   m_mu1author=0;
	   m_mu1pt=0;
	   m_mu1pts=0;
	   m_mu1pte=0;
	   m_mu1ptex=0;
	   m_mu1etae=0;
	   m_mu1phie=0;
	   m_mu1eta=0;
	   m_mu1etas=0;
	   m_mu1etaex=0;
	   m_mu1phi=0;
	   m_mu1phis=0;
	   m_mu1phiex=0;
	   m_mu1pt_vtx=0;
	   m_mu1eta_vtx=0;
	   m_mu1phi_vtx=0;
	   m_mu1pt_id=0;
	   m_mu1eta_id=0;
	   m_mu1phi_id=0;
	   m_mu1IDTrackChi2Reduct=0;
	   m_mu1IDProbChi2Fit=0;
	   m_mu1IDTrackDeltaChi2=0;
	   m_mu1Etisol10=0;
	   m_mu1Etisol20=0;
	   m_mu1Etisol30=0;
	   m_mu1Etisol40=0;
	   m_mu1Ptisol10=0;
	   m_mu1Ptisol20=0;
	   m_mu1Ptisol30=0;
	   m_mu1Ptisol40=0;
	   m_mu1Nuisol10=0;
	   m_mu1Nuisol20=0;
	   m_mu1Nuisol30=0;
	   m_mu1Nuisol40=0;
	   m_mu1d0=0;
	   m_mu1z0=0;
	   m_mu1d0PVCorr=0;
	   m_mu1z0PVCorr=0;
	   m_mu1q=0;
	   m_mu1ProbChi2Fit=0;
	   m_mu1ProbChi2Match=0;
	   m_mu1PtRatio_ex=0;
	   m_mu1DeltaPtSign_ex=0;
	   m_mu1EnergyBalance=0;
	   m_mu2index=0;
	   m_mu2author=0;
	   m_mu2pt=0;
	   m_mu2eta=0;
	   m_mu2etas=0;
	   m_mu2pts=0;
	   m_mu2pte=0;
	   m_mu2ptex=0;
	   m_mu2etae=0;
	   m_mu2phie=0;
	   m_mu2etaex=0;
	   m_mu2phi=0;
	   m_mu2phis=0;
	   m_mu2phiex=0;
	   m_mu2pt_vtx=0;
	   m_mu2eta_vtx=0;
	   m_mu2phi_vtx=0;
	   m_mu2pt_id=0;
	   m_mu2eta_id=0;
	   m_mu2phi_id=0;
	   m_mu2IDTrackChi2Reduct=0;
	   m_mu2IDProbChi2Fit=0;
	   m_mu2IDTrackDeltaChi2=0;
	   m_mu2Etisol10=0;
	   m_mu2Etisol20=0;
	   m_mu2Etisol30=0;
	   m_mu2Etisol40=0;
	   m_mu2Ptisol10=0;
	   m_mu2Ptisol20=0;
	   m_mu2Ptisol30=0;
	   m_mu2Ptisol40=0;
	   m_mu2Nuisol10=0;
	   m_mu2Nuisol20=0;
	   m_mu2Nuisol30=0;
	   m_mu2Nuisol40=0;
	   m_mu2d0=0;
	   m_mu2z0=0;
	   m_mu2d0PVCorr=0;
	   m_mu2z0PVCorr=0;
	   m_mu2q=0;
	   m_mu2ProbChi2Fit=0;
	   m_mu2ProbChi2Match=0;
	   m_mu2PtRatio_ex=0;
	   m_mu2DeltaPtSign_ex=0;
	   m_mu2EnergyBalance=0;
	   

	   int mu1_index = (*dimuon_mu1Index)[i];
	   int mu2_index = (*dimuon_mu2Index)[i];
	   
	   if (mu1_index <0 || mu2_index <0) continue; // Fix is still needed for this problem: mu_index = -1
	   if (IS_DEBUG) file <<" Pass check for muon index <0"<<endl;

	   int mu_vtx_type[2];
	  
	   mu_vtx_type[0] = GetPrimaryVertex(mu1_index);
	   mu_vtx_type[1] = GetPrimaryVertex(mu2_index);
	   
	   // CUT:: ON PRIMARY VERTEX FOR EACH MUON !!
	   if (mu_vtx_type[0] == 999 || mu_vtx_type[1] == 999 
	       || (mu_vtx_type[0] == 1 && mu_vtx_type[1] == 3)
	       || (mu_vtx_type[0] == 3 && mu_vtx_type[1] == 1))
	     {
	       file << "Discarding this dimuon: "<< endl;
	       file << "mu1 vtx type is "<<mu_vtx_type[0]<<endl;
	       file << "mu2 vtx type is "<<mu_vtx_type[1]<<endl;
	       continue;
	     }
	   
	   if (IS_DEBUG) file <<" Pass check for different vertexes"<<endl;

	   // check the family
	   Int_t mu1_author = (*dimuon_mu1Author)[i];
	   Int_t mu2_author = (*dimuon_mu2Author)[i];
	   Int_t family = -1;

	   Bool_t isStaco = kFALSE;
	   Bool_t isMuid  = kFALSE;
	   Bool_t isTrackMuon = kFALSE;
	   Bool_t isMuonMuon = kFALSE;

	   if ((mu1_author == 6 || mu1_author == 7) && (mu2_author == 6 || mu2_author == 7)) 
	     { // STACO family muons
	       isStaco = kTRUE;
	       isMuonMuon = kTRUE;
	       if (mu1_author + mu2_author == 12) family = STACO_COMBCOMB;
	       else if (mu1_author + mu2_author == 13) family = STACO_COMBTAG;
	       else family = STACO_TAGTAG;
	     }
	   else if ((mu1_author == 12 || mu1_author == 13 || mu1_author == 18) && (mu2_author == 12 || mu2_author == 13 || mu2_author == 18)) 
	     { // Muid family muons
	       isMuid = kTRUE;
	       isMuonMuon = kTRUE;
	       if (mu1_author + mu2_author == 24) family = MUID_COMBCOMB;
	       else if (mu1_author + mu2_author == 25 || mu1_author + mu2_author == 30) family = MUID_COMBTAG;
	       else family = MUID_TAGTAG;
	     }
	   else if ((mu1_author == -1 && (mu2_author == 6 || mu2_author == 7)) || ((mu1_author == 6 || mu1_author == 7) && mu2_author == -1)) 
	     { // trk + STACO family muon
	       isStaco = kTRUE;
	       isTrackMuon = kTRUE;
	       if (mu1_author + mu2_author == 5) family = ID_STACO_COMBINED;
	       else family = ID_STACO_TAGGED;
	     }
	   else if ((mu1_author == -1 && (mu2_author == 12 || mu2_author == 13 || mu2_author == 18)) || ((mu1_author == 12 || mu1_author == 13 || mu1_author == 18) && mu2_author == -1)) 
	     { // trk + Muid family muon
	       isMuid = kTRUE;
	       isTrackMuon = kTRUE;
	       if (mu1_author + mu2_author == 11) family = ID_MUID_COMBINED;
	       else family = ID_MUID_TAGGED;
	     }
	   else if (mu1_author == -1 && mu2_author == -1) family = ID_ID;
	   
	   
	   // if (!doStaco && isStaco) continue;
// 	   if (!doMuid && isMuid) continue;
// 	   if (!doTrackTrack && family == ID_ID) continue;
// 	   if (!doTrackMuon  && isTrackMuon) continue;
// 	   if (!doMuonMuon  && isMuonMuon) continue;
	   
	   if (IS_DEBUG) file <<" Pass at least 1 family is ok "<<endl;

	   // CUT:: Baseline Muon quality
	   Bool_t good_mu1 = false; 
	   Bool_t good_mu2 = false; 
	   int mu1id= mu_IDindex->at(mu1_index);
	   int mu2id= mu_IDindex->at(mu2_index);
	   if ( ID_nSCTHits->at(mu1id) >= 6 && ID_nPixelHits->at(mu1id) >= 1) good_mu1 = true; 
	   if ( ID_nSCTHits->at(mu2id) >= 6 && ID_nPixelHits->at(mu2id) >= 1) good_mu2 = true; 
	   if (!(good_mu1 && good_mu2)) continue;
	   
	   if (IS_DEBUG) file <<" Pass check on good muons inside the dimuons (baseline)"<<endl;

	   /*
	   // Set the right pointer (trk or muon) to have a symmetric access to (muon, IDtrack) parameters
	   std::vector<float> *mu1_d0  = mu_d0corrPV;
	   std::vector<float> *mu1_z0  = mu_z0;	
	   std::vector<int>   *mu1_Q   = mu_Q;	
	   std::vector<float> *mu1_pt  = mu_Pt;	
	   std::vector<float> *mu1_eta = mu_Eta;
	   std::vector<float> *mu1_phi = mu_Phi;
	   std::vector<float> *mu2_d0  = mu_d0corrPV;
	   std::vector<float> *mu2_z0  = mu_z0;	
	   std::vector<int>   *mu2_Q   = mu_Q;	
	   std::vector<float> *mu2_pt  = mu_Pt;	
	   std::vector<float> *mu2_eta = mu_Eta;
	   std::vector<float> *mu2_phi = mu_Phi;
	   
	   // create the objects
	 
	   TLorentzVector mu1, mu2, mu1_vtx, mu2_vtx, mu1_id, mu2_id;
	   TLorentzVector sum, sum_novtx;
   
	   mu1.SetPtEtaPhiM((*mu1_pt)[mu1_index], (*mu1_eta)[mu1_index], (*mu1_phi)[mu1_index], muonPDGmass);
	   mu1_vtx.SetXYZM((*dimuon_mu1CorrPx)[i], (*dimuon_mu1CorrPy)[i], (*dimuon_mu1CorrPz)[i], muonPDGmass);

	   if (mu1_author != -1)
	     mu1_id.SetPtEtaPhiM((*ID_TrackPt)[mu1id], (*ID_TrackTheta)[mu1id], (*ID_TrackPhi)[mu1id], muonPDGmass);
	   else mu1_id = mu1;
	   
	   mu2.SetPtEtaPhiM((*mu2_pt)[mu2_index], (*mu2_eta)[mu2_index], (*mu2_phi)[mu2_index], muonPDGmass);
	   mu2_vtx.SetXYZM((*dimuon_mu2CorrPx)[i], (*dimuon_mu2CorrPy)[i], (*dimuon_mu2CorrPz)[i], muonPDGmass);

	   if (mu2_author != -1)
	     mu2_id.SetPtEtaPhiM((*ID_TrackPt)[mu2id], (*ID_TrackPt)[mu2id], (*ID_TrackPt)[mu2id], muonPDGmass);
	   else mu2_id = mu2;
	   	   
	   sum = mu1_vtx + mu2_vtx;
	   sum_novtx = mu1 + mu2;
	   */

	   // Trigger: Start here
	   
	   m_muL1MU0pass=0;
	   m_muL1MU6pass=0;
	   m_muL1MU10pass=0;
	   m_muL1MU0COMMpass=0;
	   m_muEFMU4pass=0;	   
	   m_muEFMU6pass=0;	   
	   m_muEFMU10pass=0;	   
	   m_muEFMU4MSONLYpass=0;
	   
	   m_mu1L1MU0dr=100;	       
	   m_mu1L1MU0index=100;      
	   m_mu1L1MU6dr=100;	       
	   m_mu1L1MU6index=100;      
	   m_mu1L1MU10dr=100;	       
	   m_mu1L1MU10index=100;     
	   m_mu1L1MU0COMMdr=100;     
	   m_mu1L1MU0COMMindex=100; 
 
	   m_mu1L1MU0dr2=100;	       
	   m_mu1L1MU0index2=100;      
	   m_mu1L1MU6dr2=100;	       
	   m_mu1L1MU6index2=100;      
	   m_mu1L1MU10dr2=100;	       
	   m_mu1L1MU10index2=100;     
	   m_mu1L1MU0COMMdr2=100;     
	   m_mu1L1MU0COMMindex2=100;  
	   
	   m_mu1EFMU4dr=100;	       
	   m_mu1EFMU4index=100;      
	   m_mu1EFMU6dr=100;	       
	   m_mu1EFMU6index=100;      
	   m_mu1EFMU10dr=100;	       
	   m_mu1EFMU10index=100;     
	   m_mu1EFMU4MSONLYdr=100;   
	   m_mu1EFMU4MSONLYindex=100;
	   
	   m_mu2L1MU0dr=100;	       
	   m_mu2L1MU0index=100;      
	   m_mu2L1MU6dr=100;	       
	   m_mu2L1MU6index=100;      
	   m_mu2L1MU10dr=100;	       
	   m_mu2L1MU10index=100;     
	   m_mu2L1MU0COMMdr=100;     
	   m_mu2L1MU0COMMindex=100; 

	   m_mu2L1MU0dr2=100;	       
	   m_mu2L1MU0index2=100;      
	   m_mu2L1MU6dr2=100;	       
	   m_mu2L1MU6index2=100;      
	   m_mu2L1MU10dr2=100;	       
	   m_mu2L1MU10index2=100;     
	   m_mu2L1MU0COMMdr2=100;     
	   m_mu2L1MU0COMMindex2=100;

	   m_mu2EFMU4dr=100;	       
	   m_mu2EFMU4index=100;      
	   m_mu2EFMU6dr=100;	       
	   m_mu2EFMU6index=100;      
	   m_mu2EFMU10dr=100;	       
	   m_mu2EFMU10index=100;     
	   m_mu2EFMU4MSONLYdr=100;   
	   m_mu2EFMU4MSONLYindex=100;
	   
	   L1_items["L1_MU0"]=0;
	   L1_items["L1_MU6"]=0;
	   L1_items["L1_MU10"]=0;
	   L1_items["L1_MU0_COMM"]=0;
	   
	   EF_items["EF_mu4"]=0;
	   EF_items["EF_mu6"]=0;
	   EF_items["EF_mu10"]=0;
	   EF_items["EF_mu4_MSonly_MB2_noL2_EFFS"]=0;
	   
	   mapType::const_iterator it = L1_items.find("L1_MU0");
	   if (IS_DEBUG) file <<"Trigger: Reset L1MU0 status is: "<<it->second<<endl;
	   
	   
	   if (L1_ispassed->size() >0 )
	     {
	       if (IS_DEBUG) file <<"L1 is passed"<<endl;
	       mapType::const_iterator end = L1_items.end(); 
	       for( it = L1_items.begin(); it != end; ++it)
		 {
		   //file << " requested item is "<<it->first<<endl;
		   int passed=0;
		   bool L1_item_isfound=0;
		   while (passed < L1_ispassed->size() &&  !L1_item_isfound)
		     {
		       //file << " Item is : "<< L1_ispassed->at(passed)<<endl;
		       if (L1_ispassed->at(passed) == it->first )
			 { 
			   if (IS_DEBUG) file <<"Requested Chain " << it->first << " is OK" <<endl;
			   L1_item_isfound = 1 ;
			   string t = it->first;
			   L1_items[t]=1;
			 } 
		       passed++;
			   }   
		 }
	     } // end of LVl1
	   
	   if (EF_ispassed->size() >0 )
	     {
	       if (IS_DEBUG) file <<"EF is passed"<<endl;
	       //file <<"EF item trovati hanno size = "<<EF_ispassed->size()<<endl;
	       mapType::const_iterator end = EF_items.end(); 
	       for( it = EF_items.begin(); it != end; ++it)
		 {
		   //file << " requested EF item is "<<it->first<<endl;
		   int passed=0;
		   bool EF_item_isfound=0;
		   while (passed < EF_ispassed->size() &&  !EF_item_isfound)
		     {
		       //file << " Item is : "<< EF_ispassed->at(passed)<<endl;
		       if (EF_ispassed->at(passed) == it->first )
			 { 
			   if (IS_DEBUG) file <<"Requested Chain " << it->first << " is OK" <<endl;
			   EF_item_isfound = 1 ;
			   string t = it->first;
			   EF_items[t]=1;
			 } 
		       passed++;
		     }   
		 }
	     }
	       
	       //
	       // Trigger: matching
	       //
	       int index_trigger_chain=1000;
	       int index_trigger_chain2=1000;
	       bool  IsNotCombined[2];
	       IsNotCombined[0]=0;
	       IsNotCombined[1]=0;
	       float L1_DR=1000;
	       float L1_DR2=1000;
	       float EF_DR=1000;

	       if (mu1_author != 12 && mu1_author != 6) IsNotCombined[0] = 1; 
	       if (mu2_author != 12 && mu2_author != 6) IsNotCombined[1] = 1; 

		   ////// L1MU0
	       it = L1_items.find("L1_MU0");
	       string chain = it->first;
	       
	       if (it->second == 1) 
		 {
		   m_muL1MU0pass = 1;  //L1MU0 is fired in this event
		   if (IS_DEBUG) file <<" L1mu0 is fired in the event"<<endl;
		   bool IsNotComb = IsNotCombined[0];

		   if (IS_DEBUG) file <<" muon1 author is "<<mu1_author<<endl;
		   if (IS_DEBUG) file <<" IsNot Combined = "<<IsNotCombined[0]<<endl;
		   if (IS_DEBUG) file <<" Make Match with mu1_index = "<<mu1_index<<endl;
		   MatchL1Chain(mu1_index, L1_DR, L1_DR2, IsNotComb, index_trigger_chain, index_trigger_chain2, chain);
		   if (IS_DEBUG) file <<"Dr is "<<L1_DR<<endl;
		   if (IS_DEBUG) file <<"Dr2 is "<<L1_DR2<<endl;
		   m_mu1L1MU0dr     = L1_DR;
		   m_mu1L1MU0dr2    = L1_DR2;
		   m_mu1L1MU0index  = index_trigger_chain;
		   m_mu1L1MU0index2 = index_trigger_chain2;
		   
		   if (IS_DEBUG) file <<"Assigned index1= "<<m_mu1L1MU0index <<endl;
		   if (IS_DEBUG) file <<"Put  mu1dr =  "<<m_mu1L1MU0dr<<endl;
		   
		   L1_DR=1000;
		   L1_DR2=1000;
		   index_trigger_chain=1000;
		   index_trigger_chain2=1000;
		   
		   IsNotComb = IsNotCombined[1];
		   if (IS_DEBUG) file <<" muon2 author is "<<mu2_author<<endl;
		   if (IS_DEBUG) file <<" IsNot Combined = "<<IsNotCombined[1]<<endl;
		   if (IS_DEBUG) file <<" Make Match with mu1_index = "<<mu1_index<<endl;
		   MatchL1Chain(mu2_index, L1_DR, L1_DR2, IsNotComb, index_trigger_chain, index_trigger_chain2, chain);
		   if (IS_DEBUG) file <<"Dr is " << L1_DR<<endl;
		   if (IS_DEBUG) file <<"Dr2 is "<< L1_DR2<<endl;
		   m_mu2L1MU0dr     = L1_DR;
		   m_mu2L1MU0dr2    = L1_DR2;
		   m_mu2L1MU0index  = index_trigger_chain;
		   m_mu2L1MU0index2 = index_trigger_chain2;
		   if (IS_DEBUG) file <<"Assigned index2= "<<m_mu2L1MU0index <<endl;
		 }
		   
		   ////// L1MU6
		  it = L1_items.find("L1_MU6");
		  chain = it->first;
		  L1_DR=1000;
		  L1_DR2=1000;
		  index_trigger_chain=1000;
		  index_trigger_chain2=1000;
		   
		  if (it->second == 1) 
		     {
		       m_muL1MU6pass = 1;  //L1MU6 is fired in this event
		       bool IsNotComb = IsNotCombined[0];
		       
		         if (IS_DEBUG) file <<" muon1 author is "<<mu1_author<<endl;
			 if (IS_DEBUG) file <<" IsNot Combined = "<<IsNotCombined[0]<<endl;
			 if (IS_DEBUG) file <<" Make Match with mu1_index = "<<mu1_index<<endl;
			 MatchL1Chain(mu1_index, L1_DR, L1_DR2, IsNotComb, index_trigger_chain, index_trigger_chain2, chain);
			 if (IS_DEBUG) file <<"Dr is "<<L1_DR<<endl;
			 if (IS_DEBUG) file <<"Dr2 is "<<L1_DR2<<endl;
			 m_mu1L1MU6dr     = L1_DR;
			 m_mu1L1MU6dr2    = L1_DR2;
			 m_mu1L1MU6index  = index_trigger_chain;
			 m_mu1L1MU6index2 = index_trigger_chain2;
			 
			 if (IS_DEBUG) file <<"Assigned index1= "<<m_mu1L1MU6index <<endl;
			 if (IS_DEBUG) file <<"Put  mu1dr =  "<<m_mu1L1MU6dr<<endl;
			
			 L1_DR=1000;
			 L1_DR2=1000;
			 index_trigger_chain=1000;
			 index_trigger_chain2=1000;
			
			 
			 IsNotComb = IsNotCombined[1];
			 if (IS_DEBUG) file <<" muon2 author is "<<mu2_author<<endl;
			 if (IS_DEBUG) file <<" IsNot Combined = "<<IsNotCombined[1]<<endl;
			 if (IS_DEBUG) file <<" Make Match with mu1_index = "<<mu1_index<<endl;
			 MatchL1Chain(mu2_index, L1_DR, L1_DR2, IsNotComb, index_trigger_chain, index_trigger_chain2, chain);
			 if (IS_DEBUG) file <<"Dr is " << L1_DR<<endl;
			 if (IS_DEBUG) file <<"Dr2 is "<< L1_DR2<<endl;
			 m_mu2L1MU6dr     = L1_DR;
			 m_mu2L1MU6dr2    = L1_DR2;
			 m_mu2L1MU6index  = index_trigger_chain;
			 m_mu2L1MU6index2 = index_trigger_chain2;
			 if (IS_DEBUG) file <<"Assigned index2= "<<m_mu2L1MU6index <<endl;
		   
		     }
		   

		   ////// L1MU10
		    it = L1_items.find("L1_MU10");
		    chain = it->first;
		    L1_DR=1000;
		    L1_DR2=1000;
		    index_trigger_chain=1000;
		    index_trigger_chain2=1000;
		    
		   if (it->second == 1) 
		     {
		       m_muL1MU10pass = 1;  //L1MU10 is fired in this event
		       bool IsNotComb = IsNotCombined[0];

		       if (IS_DEBUG) file <<" muon1 author is "<<mu1_author<<endl;
			 if (IS_DEBUG) file <<" IsNot Combined = "<<IsNotCombined[0]<<endl;
			 if (IS_DEBUG) file <<" Make Match with mu1_index = "<<mu1_index<<endl;
			 MatchL1Chain(mu1_index, L1_DR, L1_DR2, IsNotComb, index_trigger_chain, index_trigger_chain2, chain);
			 if (IS_DEBUG) file <<"Dr is "<<L1_DR<<endl;
			 if (IS_DEBUG) file <<"Dr2 is "<<L1_DR2<<endl;
			 m_mu1L1MU10dr     = L1_DR;
			 m_mu1L1MU10dr2    = L1_DR2;
			 m_mu1L1MU10index  = index_trigger_chain;
			 m_mu1L1MU10index2 = index_trigger_chain2;
			 
			 if (IS_DEBUG) file <<"Assigned index1= "<<m_mu1L1MU10index <<endl;
			 if (IS_DEBUG) file <<"Put  mu1dr =  "<<m_mu1L1MU10dr<<endl;
			 
			 L1_DR=1000;
			 L1_DR2=1000;
			 index_trigger_chain=1000;
			 index_trigger_chain2=1000;
			
			 IsNotComb = IsNotCombined[1];
			 if (IS_DEBUG) file <<" muon2 author is "<<mu2_author<<endl;
			 if (IS_DEBUG) file <<" IsNot Combined = "<<IsNotCombined[1]<<endl;
			 if (IS_DEBUG) file <<" Make Match with mu1_index = "<<mu1_index<<endl;
			 MatchL1Chain(mu2_index, L1_DR, L1_DR2, IsNotComb, index_trigger_chain, index_trigger_chain2, chain);
			 if (IS_DEBUG) file <<"Dr is " << L1_DR<<endl;
			 if (IS_DEBUG) file <<"Dr2 is "<< L1_DR2<<endl;
			 m_mu2L1MU10dr     = L1_DR;
			 m_mu2L1MU10dr2    = L1_DR2;
			 m_mu2L1MU10index  = index_trigger_chain;
			 m_mu2L1MU10index2 = index_trigger_chain2;
			 if (IS_DEBUG) file <<"Assigned index2= "<<m_mu2L1MU10index <<endl;
		     }

		  ////// L1MU0_COMM
		  it = L1_items.find("L1_MU0_COMM");
		  chain = it->first;
		  L1_DR=1000;
		  L1_DR2=1000;
		  index_trigger_chain=1000;
		  index_trigger_chain2=1000;
		  
		  if (it->second == 1) 
		     {
		       m_muL1MU0COMMpass = 1;  //L1MU0_COMM is fired in this event
		       bool IsNotComb = IsNotCombined[0];
		       
		        if (IS_DEBUG) file <<" muon1 author is "<<mu1_author<<endl;
			 if (IS_DEBUG) file <<" IsNot Combined = "<<IsNotCombined[0]<<endl;
			 if (IS_DEBUG) file <<" Make Match with mu1_index = "<<mu1_index<<endl;
			 MatchL1Chain(mu1_index, L1_DR, L1_DR2, IsNotComb, index_trigger_chain, index_trigger_chain2, chain);
			 if (IS_DEBUG) file <<"Dr is "<<L1_DR<<endl;
			 if (IS_DEBUG) file <<"Dr2 is "<<L1_DR2<<endl;
			 m_mu1L1MU0COMMdr     = L1_DR;
			 m_mu1L1MU0COMMdr2    = L1_DR2;
			 m_mu1L1MU0COMMindex  = index_trigger_chain;
			 m_mu1L1MU0COMMindex2 = index_trigger_chain2;
			 
			 if (IS_DEBUG) file <<"Assigned index1= "<<m_mu1L1MU0COMMindex <<endl;
			 if (IS_DEBUG) file <<"Put  mu1dr =  "<<m_mu1L1MU0COMMdr<<endl;
			 
			 L1_DR=1000;
			 L1_DR2=1000;
			 index_trigger_chain=1000;
			 index_trigger_chain2=1000;
			 
			 IsNotComb = IsNotCombined[1];
			 if (IS_DEBUG) file <<" muon2 author is "<<mu2_author<<endl;
			 if (IS_DEBUG) file <<" IsNot Combined = "<<IsNotCombined[1]<<endl;
			 if (IS_DEBUG) file <<" Make Match with mu1_index = "<<mu1_index<<endl;
			 MatchL1Chain(mu2_index, L1_DR, L1_DR2, IsNotComb, index_trigger_chain, index_trigger_chain2, chain);
			 if (IS_DEBUG) file <<"Dr is " << L1_DR<<endl;
			 if (IS_DEBUG) file <<"Dr2 is "<< L1_DR2<<endl;
			 m_mu2L1MU0COMMdr     = L1_DR;
			 m_mu2L1MU0COMMdr2    = L1_DR2;
			 m_mu2L1MU0COMMindex  = index_trigger_chain;
			 m_mu2L1MU0COMMindex2 = index_trigger_chain2;
			 if (IS_DEBUG) file <<"Assigned index2= "<<m_mu2L1MU0COMMindex <<endl;


		     }

		   ////// EF_MU4
		   it = EF_items.find("EF_mu4");
		   chain = it->first;
		   L1_DR=1000;
		   L1_DR2=1000;
		   index_trigger_chain=1000;
		   index_trigger_chain2=1000;
			 
		   if (it->second == 1) 
		     {
		       m_muEFMU4pass = 1;  //EFMU0_COMM is fired in this event
		       bool IsNotComb = IsNotCombined[0];
		       MatchEFChain(mu1_index, EF_DR, IsNotComb, index_trigger_chain, chain);
		       m_mu1EFMU4dr= EF_DR;
		       m_mu1EFMU4index= index_trigger_chain;
		       
		       L1_DR=1000;
		       L1_DR2=1000;
		       index_trigger_chain=1000;
		       index_trigger_chain2=1000;
		       
		       IsNotComb = IsNotCombined[1];
		       MatchEFChain(mu2_index, EF_DR, IsNotComb, index_trigger_chain, chain);
		       m_mu2EFMU4dr= EF_DR;
		       m_mu2EFMU4index= index_trigger_chain;
		     }
		   
		   ////// EF_MU6
		   it = EF_items.find("EF_mu6");
		   chain = it->first;
		   L1_DR=1000;
		   L1_DR2=1000;
		   index_trigger_chain=1000;
		   index_trigger_chain2=1000;
		   
		   if (it->second == 1) 
		     {
		       file<<" Found Event Passing " << chain << endl;
		       m_muEFMU6pass = 1;  
		       bool IsNotComb = IsNotCombined[0];
		       MatchEFChain(mu1_index, EF_DR, IsNotComb, index_trigger_chain, chain);
		       file<<" First muon1 DR= " << EF_DR << " index=" << index_trigger_chain<< endl;
		       m_mu1EFMU6dr= EF_DR;
		       m_mu1EFMU6index= index_trigger_chain;
		       L1_DR=1000;
		       L1_DR2=1000;
		       index_trigger_chain=1000;
		       index_trigger_chain2=1000;
		       IsNotComb = IsNotCombined[1];
		       MatchEFChain(mu2_index, EF_DR, IsNotComb, index_trigger_chain, chain);
		       file<<" Second muon2 DR= " << EF_DR << " index=" << index_trigger_chain<< endl;
		       m_mu2EFMU6dr= EF_DR;
		       m_mu2EFMU6index= index_trigger_chain;
		     }
		   
		    ////// EF_MU10
		   it = EF_items.find("EF_mu10");
		   chain = it->first;
		   L1_DR=1000;
		   L1_DR2=1000;
		   index_trigger_chain=1000;
		   index_trigger_chain2=1000;
		   
		   if (it->second == 1) 
		     {
		       m_muEFMU10pass = 1;  //EFMU0_COMM is fired in this event
		       bool IsNotComb = IsNotCombined[0];
		       MatchEFChain(mu1_index, EF_DR, IsNotComb, index_trigger_chain, chain);
		       m_mu1EFMU10dr= EF_DR;
		       m_mu1EFMU10index= index_trigger_chain;
		       L1_DR=1000;
		       L1_DR2=1000;
		       index_trigger_chain=1000;
		       index_trigger_chain2=1000;
		       IsNotComb = IsNotCombined[1];
		       MatchEFChain(mu2_index, EF_DR, IsNotComb, index_trigger_chain, chain);
		       m_mu2EFMU10dr= EF_DR;
		       m_mu2EFMU10index= index_trigger_chain;
		     }
		   
		    ////// EF_MU4MSONLY
		   it = EF_items.find("EF_mu4_MSonly_MB2_noL2_EFFS");
		   chain = it->first;
		   L1_DR=1000;
		   L1_DR2=1000;
		   index_trigger_chain=1000;
		   index_trigger_chain2=1000;
		   
		   if (it->second == 1) 
		     {
		       m_muEFMU4MSONLYpass = 1;  //EFMU0_COMM is fired in this event
		       bool IsNotComb = IsNotCombined[0];
		       MatchEFChain(mu1_index, EF_DR, IsNotComb, index_trigger_chain, chain);
		       m_mu1EFMU4MSONLYdr= EF_DR;
		       m_mu1EFMU4MSONLYindex= index_trigger_chain;
		       L1_DR=1000;
		       L1_DR2=1000;
		       index_trigger_chain=1000;
		       index_trigger_chain2=1000;
		       IsNotComb = IsNotCombined[1];
		       MatchEFChain(mu2_index, EF_DR, IsNotComb, index_trigger_chain, chain);
		       m_mu2EFMU4MSONLYdr= EF_DR;
		       m_mu2EFMU4MSONLYindex= index_trigger_chain;
		     }
		   



		   // mu1id mu2id;
	   // additional variables
	   int etarange = -1;
	   float eta_mu1id = -log(fabs(tan(ID_TrackTheta->at(mu1id)))/2.);
	   float eta_mu2id = -log(fabs(tan(ID_TrackTheta->at(mu2id)))/2.);

	   if      (TMath::Abs(eta_mu1id) <  1.05 && TMath::Abs(eta_mu2id) <  1.05) etarange = BB;
	   else if (TMath::Abs(eta_mu1id) >= 1.05 && TMath::Abs(eta_mu2id) >= 1.05) etarange = EE;
	   else etarange = BE;
	   
	   // fill the tree
	   //m_invmass_novtx = sum_novtx.M();
	   m_invmass = (*dimuon_invMass)[i];
	   m_q = (mu_Q->at(mu1_index)+mu_Q->at(mu2_index));
	   m_pt = (*dimuon_pt)[i];
	   m_deltaEta = mu_Eta->at(mu1_index)-mu_Eta->at(mu2_index);
	   m_deltaPhi = mu_Phi->at(mu1_index)-mu_Phi->at(mu2_index);
	   m_deltaPt  = mu_Pt->at(mu1_index)-mu_Pt->at(mu2_index);
	   m_family = family;
	   m_etarange = etarange; // start from here

	   // Mu1 vars
	   m_mu1index = mu1_index;
	   m_mu1author = mu1_author;
	   m_mu1pt = mu_Pt->at(mu1_index);
	   m_mu1eta = mu_Eta->at(mu1_index);
	   m_mu1phi = mu_Phi->at(mu1_index);
	   //m_mu1pt_vtx = mu1_vtx.Pt();
	   //m_mu1eta_vtx = mu1_vtx.Eta();
	   //m_mu1phi_vtx = mu1_vtx.Phi();
	   m_mu1pt_id = ID_TrackPt->at(mu1id);
	   m_mu1eta_id = log(fabs(tan(ID_TrackTheta->at(mu1id)))/2.);
	   m_mu1phi_id = ID_TrackPhi->at(mu1id);
	   m_mu1IDTrackChi2Reduct=ID_TrackChi2Reduct->at(mu1id);
	   m_mu1IDProbChi2Fit=ID_ProbChi2Fit->at(mu1id);
	   m_mu1IDTrackDeltaChi2=ID_TrackDeltaChi2->at(mu1id);
	   m_mu1etas     = mu_Eta_ms->at(mu1_index);      	   
	   m_mu1phis     = mu_Phi_ms->at(mu1_index); 		   
	   m_mu1pts      = mu_Pt_ms->at(mu1_index); 		   
	   m_mu1etae     = mu_IDEtaExtrapolatedToMS->at(mu1_index);  // extrapolator vars    
	   m_mu1phie     = mu_IDPhiExtrapolatedToMS->at(mu1_index);      
	   m_mu1pte      = mu_IDPtExtrapolatedToMS->at(mu1_index); 
	   m_mu1etaex    = mu_Eta_ex->at(mu1_index);       
	   m_mu1phiex    = mu_Phi_ex->at(mu1_index);       
	   m_mu1ptex     = mu_Pt_ex->at(mu1_index); 
	   
	   /// isolation
	   m_mu1Etisol10 = mu_etIsol10->at(mu1_index); 
	   m_mu1Etisol20 = mu_etIsol20->at(mu1_index);
	   m_mu1Etisol30 = mu_etIsol30->at(mu1_index);
	   m_mu1Etisol40 = mu_etIsol40->at(mu1_index);
	   m_mu1Ptisol10 = mu_ptIsol10->at(mu1_index); 
	   m_mu1Ptisol20 = mu_ptIsol20->at(mu1_index);
	   m_mu1Ptisol30 = mu_ptIsol30->at(mu1_index);
	   m_mu1Ptisol40 = mu_ptIsol40->at(mu1_index);
	   m_mu1Nuisol10 = mu_nuIsol10->at(mu1_index); 
	   m_mu1Nuisol20 = mu_nuIsol20->at(mu1_index);
	   m_mu1Nuisol30 = mu_nuIsol30->at(mu1_index);
	   m_mu1Nuisol40 = mu_nuIsol40->at(mu1_index);
	   // d0 and z0
	   m_mu1d0       = mu_d0->at(mu1_index);       
	   m_mu1z0       = mu_z0->at(mu1_index);       
	   m_mu1d0PVCorr = mu_d0corrPV->at(mu1_index); 
	   m_mu1z0PVCorr = mu_z0corrPV->at(mu1_index);
	   // charge
	   m_mu1q = mu_Q->at(mu1_index);
	   
	   // Pi k vars
	   /*
	   m_mu1ProbChi2Fit=mu_ProbChi2Fit->at(mu1_index);
	   m_mu1ProbChi2Match=mu_ProbChi2Match->at(mu1_index);
	   m_mu1PtRatio_ex=mu_PtRatio_ex->at(mu1_index);
	   m_mu1DeltaPtSign_ex=mu_DeltaPtSign_ex->at(mu1_index);
	   m_mu1EnergyBalance=mu_EnergyBalance->at(mu1_index);
	   */
	   
	   // MU2 vars
	   m_mu2index = mu2_index;
	   m_mu2author = mu2_author;
	   m_mu2pt = mu_Pt->at(mu2_index);
	   m_mu2eta = mu_Eta->at(mu2_index);
	   m_mu2phi = mu_Phi->at(mu2_index);
	   //m_mu2pt_vtx = mu2_vtx.Pt();
	   //m_mu2eta_vtx = mu2_vtx.Eta();
	   //m_mu2phi_vtx = mu2_vtx.Phi();
	   m_mu2pt_id = ID_TrackPt->at(mu2id);
	   m_mu2eta_id = log(fabs(tan(ID_TrackTheta->at(mu2id)))/2.);
	   m_mu2phi_id = ID_TrackPhi->at(mu2id);
	   m_mu2IDTrackChi2Reduct=ID_TrackChi2Reduct->at(mu2id);
	   m_mu2IDProbChi2Fit=ID_ProbChi2Fit->at(mu2id);
	   m_mu2IDTrackDeltaChi2=ID_TrackDeltaChi2->at(mu2id);
	   m_mu2etas     = mu_Eta_ms->at(mu2_index);      	   
	   m_mu2phis     = mu_Phi_ms->at(mu2_index); 		   
	   m_mu2pts      = mu_Pt_ms->at(mu2_index); 		   
	   m_mu2etae     = mu_IDEtaExtrapolatedToMS->at(mu2_index);  // extrapolator vars    
	   m_mu2phie     = mu_IDPhiExtrapolatedToMS->at(mu2_index);      
	   m_mu2pte      = mu_IDPtExtrapolatedToMS->at(mu2_index); 
	   m_mu2etaex    = mu_Eta_ex->at(mu2_index);       
	   m_mu2phiex    = mu_Phi_ex->at(mu2_index);       
	   m_mu2ptex     = mu_Pt_ex->at(mu2_index);

	   // Pi k vars
	   /*
	   m_mu2ProbChi2Fit=mu_ProbChi2Fit->at(mu2_index);
	   m_mu2ProbChi2Match=mu_ProbChi2Match->at(mu2_index);
	   m_mu2PtRatio_ex=mu_PtRatio_ex->at(mu2_index);
	   m_mu2DeltaPtSign_ex=mu_DeltaPtSign_ex->at(mu2_index);
	   m_mu2EnergyBalance=mu_EnergyBalance->at(mu2_index);
	   */

	   /// isolation
	   m_mu2Etisol10 = mu_etIsol10->at(mu2_index); 
	   m_mu2Etisol20 = mu_etIsol20->at(mu2_index);
	   m_mu2Etisol30 = mu_etIsol30->at(mu2_index);
	   m_mu2Etisol40 = mu_etIsol40->at(mu2_index);
	   m_mu2Ptisol10 = mu_ptIsol10->at(mu2_index); 
	   m_mu2Ptisol20 = mu_ptIsol20->at(mu2_index);
	   m_mu2Ptisol30 = mu_ptIsol30->at(mu2_index);
	   m_mu2Ptisol40 = mu_ptIsol40->at(mu2_index);
	   m_mu2Nuisol10 = mu_nuIsol10->at(mu2_index); 
	   m_mu2Nuisol20 = mu_nuIsol20->at(mu2_index);
	   m_mu2Nuisol30 = mu_nuIsol30->at(mu2_index);
	   m_mu2Nuisol40 = mu_nuIsol40->at(mu2_index);
	   // d0 and z0
	   m_mu2d0       = mu_d0->at(mu2_index);       
	   m_mu2z0       = mu_z0->at(mu2_index);       
	   m_mu2d0PVCorr = mu_d0corrPV->at(mu2_index); 
	   m_mu2z0PVCorr = mu_z0corrPV->at(mu2_index);
	   // charge
	   m_mu2q = mu_Q->at(mu2_index);
	   
	   dump->Fill();
	   // end dimuon fill
	 }
	 }
     } 
   file << "  Done." << std::endl;
   outputfile->Write();
   outputfile->Close();
   file << " Output data written, analysis OK." << std::endl;
}

bool JpsiDumper::InclusiveBaselineSelection()
{
  int HowmanyVertexes = VtxType->size();
  int ind_v = 0;
  bool HasGoodPrimary = 0;
  while (! HasGoodPrimary && ind_v < HowmanyVertexes)
    {
      if (VtxType->at(ind_v) == 1) // Is Primary Vertex !!
	 {
	   int n_id_good_tracks = 0;
	   for (Int_t i = 0; i < ID_nTracks; i++)
	   {
	     if ( ID_VxIndex->at(i) == ind_v && ID_nSCTHits->at(i) >= 6
		  && ID_nPixelHits->at(i) >= 1) 
	       {
		 n_id_good_tracks++;
	       }
	   }
	   if (n_id_good_tracks >= 3)  HasGoodPrimary = 1;
	 }
      ind_v++;
    }
  return HasGoodPrimary ;
}

int JpsiDumper::GetPrimaryVertex(int mu_index)
{
  int type_of_vertex;
  int ind_of_vertex;
  int mu_id= mu_IDindex->at(mu_index); // index of ID track associated to muon
  if (mu_id <0) type_of_vertex = 999;
  if (mu_id >= 0) 
    {
      ind_of_vertex = ID_VxIndex->at(mu_id);
      type_of_vertex = VtxType->at(ind_of_vertex);
    }  
  return  type_of_vertex;
}

void JpsiDumper::MatchL1Chain(int imuon, float& deltaR, float& deltaR2, bool isNotCombined, int& itri, int& itri2, string ThisChain)
{
  
  float etaMS=999;
  float phiMS=999;
  float eta_vertex = mu_Eta->at(imuon);
  float phi_vertex = mu_Phi->at(imuon);

  if(! isNotCombined){
    etaMS=mu_Eta_ms->at(imuon);
    phiMS=mu_Phi_ms->at(imuon);
  }  else {
    if (etaMS=mu_IDEtaExtrapolatedToMS->at(imuon) != -999 ){	
      etaMS=mu_IDEtaExtrapolatedToMS->at(imuon);
      phiMS=mu_IDPhiExtrapolatedToMS->at(imuon);
    }    else      {	
      etaMS=mu_Eta->at(imuon);
      phiMS=mu_Phi->at(imuon);
    }
  }
  
  if (IS_DEBUG) file <<"MatchL1:: muon eta= "<< etaMS << " phi=" << phiMS << endl;
  if (IS_DEBUG) file <<"MatchL1:: Nr Objects L1 are "<<L1_trigObj_names->size() <<endl;

  float maxDR=1000;
  int ThisChainLevel = L1ThMap(ThisChain);
  if (IS_DEBUG) file <<"MatchL1: This Chain has level "<<ThisChainLevel<<endl;

  for (int ind_TriggerObj_names=0; ind_TriggerObj_names< L1_trigObj_names->size(); ind_TriggerObj_names++)
    {
      vector<string> TriggerChain_string =  L1_trigObj_names->at(ind_TriggerObj_names);
      bool FoundItem = 0 ;
      int ind= 0;
      if (IS_DEBUG) file <<"MatchL1 index::"<<ind_TriggerObj_names<< "  etaL1="<< L1_trigObj_eta->at(ind_TriggerObj_names) 
			 << " phiL1="<< L1_trigObj_phi->at(ind_TriggerObj_names)
			 << endl;
      
      while (!FoundItem && ind < TriggerChain_string.size())
	{
	  if (IS_DEBUG) file <<"Object n. "<<ind_TriggerObj_names<<" chain in = "<<TriggerChain_string[ind]<<endl;

	  if (L1ThMap(TriggerChain_string[ind]) >= ThisChainLevel)
	 	{   
		  FoundItem=1;
		  if (IS_DEBUG) file <<"Found Good Item in object n."<<ind_TriggerObj_names
				     <<" L1 item is "<<TriggerChain_string[ind]<<endl;
		}
	  ind++;
	}
      if (FoundItem) {
	// find the associated muctpi
	int imuctpi=-1;
	for(int jj=0; jj<muctpi_n; ++jj){
	  if(muctpi_sysId->at(jj)!=0 ) continue;
	  if(
	     (L1_trigObj_eta->at(ind_TriggerObj_names)==muctpi_eta->at(jj))
	     &&
	     (L1_trigObj_phi->at(ind_TriggerObj_names)==muctpi_phi->at(jj))
	     ){
	    imuctpi=jj;
	    break;
	  }
	}
	//determine the new eta/phi

	 float dr =1000;
	 
	 if(imuctpi>-1)
	   {
	     int sl=muctpi_sector->at(imuctpi);
	     int roiN=muctpi_roiNumber->at(imuctpi);
	     float newTrigEta= ( muctpiEta_max[sl][roiN]+ muctpiEta_min[sl][roiN] )/2;
	     float newTrigPhi = ( muctpiPhi_max[sl][roiN]+ muctpiPhi_min[sl][roiN] )/2;
	     
	     dr = DR(etaMS, phiMS,newTrigEta, newTrigPhi);
	     // for non combine dimuons dr is calculated also 
	     // with info @ IP at the end we take the min between the two DR
	     float dr_vertex = 10000;
	     if (! isNotCombined) dr_vertex = dr = DR(eta_vertex, phi_vertex,newTrigEta, newTrigPhi);
	     if (dr_vertex < dr) dr = dr_vertex;

	     if(IS_DEBUG) file<< "After the remapping: eta= " << newTrigEta<< " phi=" << newTrigPhi 	     << "** DR= "<< dr 
			      << endl;
	   } 
	 else 
	   {
	     dr = DR(etaMS, phiMS,L1_trigObj_eta->at(ind_TriggerObj_names) ,L1_trigObj_phi->at(ind_TriggerObj_names) );
	   }

	 if(IS_DEBUG) file<<  "** DR= "<< dr << "index=" << ind_TriggerObj_names<< endl;
	 if(dr<maxDR) 
	  {
	    deltaR2 = maxDR;
	    itri2   = itri; 
	    maxDR   = dr;
	    itri    = ind_TriggerObj_names;
	  }
      }
    }
  if(IS_DEBUG) file<<  "** Final DR= "<< maxDR << " index= " << itri<< endl;
  deltaR=maxDR;
} 

void JpsiDumper::MatchEFChain(int imuon, float& deltaR, bool isNotCombined, int& itri, string ThisChain)
{
  
  float  etaMS=mu_Eta->at(imuon);
  float  phiMS=mu_Phi->at(imuon);
  
 if(IS_DEBUG) file <<"****** muon eta= "<< etaMS << " phi=" << phiMS << endl;

 if (IS_DEBUG) file <<"*****MatchEF:: Nr Objects EF are "<<EF_trigObj_names->size() <<endl;

  float maxDR=1000;

  for (int ind_TriggerObj_names=0; ind_TriggerObj_names< EF_trigObj_names->size(); ind_TriggerObj_names++)
    {
     if(IS_DEBUG) file << "******** EF eta="<< EF_trigObj_eta->at(ind_TriggerObj_names) 
	   << " phi="<< EF_trigObj_phi->at(ind_TriggerObj_names)
	   << endl;
      vector<string> TriggerChain_string =  EF_trigObj_names->at(ind_TriggerObj_names);
      bool FoundItem = 0 ;
      int ind= 0;
      while (!FoundItem && ind < TriggerChain_string.size())
	{
	  if(IS_DEBUG) file <<"*** testing chain " << TriggerChain_string[ind]<< endl;
	  if (TriggerChain_string[ind] == ThisChain)
	 	{   
		  FoundItem=1;
		}
	  ind++;
	}
      if (FoundItem) 
	{
	  float dr = DR(etaMS, phiMS, EF_trigObj_eta->at(ind_TriggerObj_names), 
			EF_trigObj_phi->at(ind_TriggerObj_names));
	  
	  file 	   << " DR= "<< dr << endl;
	  
	  if(dr<maxDR) 
	    {
	      maxDR=dr;
	      itri=ind_TriggerObj_names;
	    }
	}
    }
  
  deltaR=maxDR;
} 

float JpsiDumper::DR(float eta1, float phi1, float eta2, float phi2){
  return sqrt(
          (eta1-eta2)*(eta1-eta2)
          +
          (phi1-phi2)*(phi1-phi2)
          );
} 

int JpsiDumper::L1ThMap(string chain)
{
  if(chain=="L1_MU0") return 1;
  if(chain=="L1_MU6") return 2;
  if(chain=="L1_MU10") return 3;
  if(chain=="L1_MU0_COMM") return 4;
  if(chain=="L1_MU15") return 5;
  if(chain=="L1_MU20") return 6; 
  return -1;
}

void JpsiDumper::FillMulti(vector< my_pair> muons_list, int algo_flag)
{
  
  m_mchargetot = 0;   
  m_mmu1eta = 0;
  m_mmu1phi = 0;
  m_mmu1pt  = 0;
  m_mmu2eta = 0;
  m_mmu2phi = 0;
  m_mmu2pt = 0;
  m_mmu3eta = 0;
  m_mmu3phi = 0;
  m_mmu3pt = 0;
  m_mmu4eta = 0;
  m_mmu4phi = 0;
  m_mmu4pt = 0;
  //  m_mmu1vertex  = 0;
  //m_mmu2vertex  = 0;
  // m_mmu3vertex  = 0;
  //m_mmu4vertex  = 0;
  m_mmu1author = 0;
  m_mmu2author = 0;
  m_mmu3author = 0;
  m_mmu4author =0;
  
     my_pair pair1 = muons_list[0];
     my_pair pair2 = muons_list[1];
     my_pair pair3 = muons_list[2];
     my_pair pair4 = muons_list[3];

     for (vector<my_pair>::iterator it = muons_list.begin(); 
       it!=muons_list.end(); ++it) 
    {
      int ind_muon = (*it).first;
      m_mchargetot += mu_Q->at(ind_muon);
    }
     
     float dmu1eta = mu_Eta->at(muons_list[0].first);
     float dmu1phi = mu_Phi->at(muons_list[0].first);
     float dmu1pt  =  mu_Pt->at(muons_list[0].first);
    
     float dmu2eta = mu_Eta->at(muons_list[1].first);
     float dmu2phi = mu_Phi->at(muons_list[1].first);
     float dmu2pt  =  mu_Pt->at(muons_list[1].first);
   
     float dmu3eta = mu_Eta->at(muons_list[2].first);
     float dmu3phi = mu_Phi->at(muons_list[2].first);
     float dmu3pt  =  mu_Pt->at(muons_list[2].first);

     float dmu4eta = mu_Eta->at(muons_list[3].first);
     float dmu4phi = mu_Phi->at(muons_list[3].first);
     float dmu4pt  =  mu_Pt->at(muons_list[3].first);

     float E1 = sqrt(muonPDGmass*muonPDGmass+muons_list[0].second*muons_list[0].second);
     float E2 = sqrt(muonPDGmass*muonPDGmass+muons_list[1].second*muons_list[1].second);
     float E3 = sqrt(muonPDGmass*muonPDGmass+muons_list[2].second*muons_list[2].second);
     float E4 = sqrt(muonPDGmass*muonPDGmass+muons_list[3].second*muons_list[3].second);
  
     float InvMassMu1=sqrt(E1*E1-((dmu1pt*cos(dmu1phi)*dmu1pt*cos(dmu1phi))+(dmu1pt*sin(dmu1phi)*dmu1pt*sin(dmu1phi))+(dmu1pt*sinh(dmu1eta)*dmu1pt*sinh(dmu1eta))));
     //cout<<"InvMassMu1 = " << InvMassMu1<<endl;
     
     
     float InvMassMu2=sqrt(E2*E2-((dmu2pt*cos(dmu2phi)*dmu2pt*cos(dmu2phi))+(dmu2pt*sin(dmu2phi)*dmu2pt*sin(dmu2phi))+(dmu2pt*sinh(dmu2eta)*dmu2pt*sinh(dmu2eta))));
     //cout<<"InvMassMu2 = " << InvMassMu2<<endl;

     float Etot=(E1+E2);
     float pxtot=(dmu1pt*cos(dmu1phi))+(dmu2pt*cos(dmu2phi));
     float pytot=(dmu1pt*sin(dmu1phi))+(dmu2pt*sin(dmu2phi));
     float pztot=(dmu1pt*sinh(dmu1eta))+(dmu2pt*sinh(dmu2eta));

     float InvMassG12=sqrt(Etot*Etot-(pxtot*pxtot+pytot*pytot+pztot*pztot));
     //cout<<"InvMassG12 = "<<InvMassG12<<endl;
     
     TLorentzVector p1(dmu1pt*cos(dmu1phi),dmu1pt*sin(dmu1phi),
		    dmu1pt*sinh(dmu1eta),E1);
     
     TLorentzVector p2(dmu2pt*cos(dmu2phi),dmu2pt*sin(dmu2phi),
		    dmu2pt*sinh(dmu2eta),E2);
     
     TLorentzVector p3(dmu3pt*cos(dmu3phi),dmu3pt*sin(dmu3phi),
		       dmu3pt*sinh(dmu3eta),E3);
     
     TLorentzVector p4(dmu4pt*cos(dmu4phi),dmu4pt*sin(dmu4phi),
		       dmu4pt*sinh(dmu4eta),E4);
         
     //m_mnmuons  = mu_n;
     m_mIsStaco = 0   ;
     if (algo_flag == 1) m_mIsStaco = 1;
     bool fill1 = true;
     bool fill2 = true;
     bool fill3 = true;
     m_mmupairtype=1;  
     GetMother(p1,p2,p3,p4,pair1,pair2,pair3,pair4, fill1);//,mu1q,mu2q,mu3q,mu4q);
     multi->Fill();
     m_mmupairtype=2;      
     GetMother(p1,p3,p2,p4,pair1,pair3,pair2,pair4, fill2);//,mu1q,mu3q,mu2q,mu4q);
     multi->Fill();
     m_mmupairtype=3;       
     GetMother(p1,p4,p2,p3,pair1,pair4,pair2,pair3, fill3);//,mu1q,mu4q,mu2q,mu3q);
     multi->Fill();
}


void JpsiDumper::GetMother(TLorentzVector mu1, TLorentzVector mu2, TLorentzVector mu3, TLorentzVector mu4,  my_pair pair1, my_pair pair2, my_pair pair3, my_pair pair4, bool& ToFill)
{
  m_mmu1eta = mu_Eta->at(pair1.first);
  m_mmu1phi = mu_Phi->at(pair1.first);
  m_mmu1pt  =  mu_Pt->at(pair1.first);

  m_mmu2eta = mu_Eta->at(pair2.first);
  m_mmu2phi = mu_Phi->at(pair2.first);
  m_mmu2pt  =  mu_Pt->at(pair2.first);
  
  m_mmu3eta = mu_Eta->at(pair3.first);
  m_mmu3phi = mu_Phi->at(pair3.first);
  m_mmu3pt  =  mu_Pt->at(pair3.first);
  
  m_mmu4eta = mu_Eta->at(pair4.first);
  m_mmu4phi = mu_Phi->at(pair4.first);
  m_mmu4pt  =  mu_Pt->at(pair4.first);

  int mu1q = mu_Q->at(pair1.first);
  int mu2q = mu_Q->at(pair2.first);
  int mu3q = mu_Q->at(pair3.first);
  int mu4q = mu_Q->at(pair4.first);
  
  int mu_id1= mu_IDindex->at(pair1.first);
  int mu_id2= mu_IDindex->at(pair2.first);
  int mu_id3= mu_IDindex->at(pair3.first);
  int mu_id4= mu_IDindex->at(pair4.first);

  m_mmu1vertex  = ID_VxIndex->at(mu_id1);
  m_mmu2vertex  = ID_VxIndex->at(mu_id2);
  m_mmu3vertex  = ID_VxIndex->at(mu_id3);
  m_mmu4vertex  = ID_VxIndex->at(mu_id4);
  
  m_mmu1vxtype  = VtxType->at(m_mmu1vertex);
  m_mmu2vxtype  = VtxType->at(m_mmu2vertex);
  m_mmu3vxtype  = VtxType->at(m_mmu3vertex);
  m_mmu4vxtype  = VtxType->at(m_mmu4vertex);
  
  m_mmu1author = mu_author->at(pair1.first);
  m_mmu2author = mu_author->at(pair2.first);
  m_mmu3author = mu_author->at(pair3.first);
  m_mmu4author = mu_author->at(pair4.first);

  TLorentzVector MotherA = mu1+mu2;
  m_mmotherAmass = sqrt(MotherA*MotherA);
  m_mmotherAcharge = mu1q+ mu2q;
  m_mmotherAp = sqrt(MotherA.Px()*MotherA.Px()+MotherA.Py()*MotherA.Py()+MotherA.Pz()*MotherA.Pz());
  float thetaA = acos(MotherA.Pz()/m_mmotherAp);
  m_mmotherApt = m_mmotherAp*sin(thetaA);
  m_mmotherAeta = -1*log(tan(thetaA/2));
  m_mmotherAphi = atan(MotherA.Py()/MotherA.Px());
  if (MotherA.Py() >0 && MotherA.Px() < 0) m_mmotherAphi = atan(MotherA.Py()/MotherA.Px()) + PiGreek;
  if (MotherA.Py() <0 && MotherA.Px() < 0) m_mmotherAphi = atan(MotherA.Py()/MotherA.Px()) - PiGreek;

  // if(m_mmupairtype==1)cout<<"InvMAss MotherA = "<<m_mmotherAmass<<endl;

  TLorentzVector MotherB = mu3+mu4;
  m_mmotherBmass = sqrt(MotherB*MotherB);
  m_mmotherBcharge = mu3q+ mu4q;
  m_mmotherBp = sqrt(MotherB.Px()*MotherB.Px()+MotherB.Py()*MotherB.Py()+MotherB.Pz()*MotherB.Pz());
  float thetaB = acos(MotherB.Pz()/m_mmotherBp);
  m_mmotherBpt = m_mmotherBp*sin(thetaB);
  m_mmotherBeta = -1*log(tan(thetaB/2));
  m_mmotherBphi = atan(MotherB.Py()/MotherB.Px());
  if (MotherB.Py() >0 && MotherB.Px() < 0) m_mmotherBphi = atan(MotherB.Py()/MotherB.Px()) + PiGreek;
  if (MotherB.Py() <0 && MotherB.Px() < 0) m_mmotherBphi = atan(MotherB.Py()/MotherB.Px()) - PiGreek;

  TLorentzVector GranMother = MotherA+ MotherB;
  m_mgmothermass = sqrt(GranMother*GranMother);
  m_mgmothercharge = m_mmotherAcharge+ m_mmotherBcharge;
  m_mgmotherp = sqrt(GranMother.Px()*GranMother.Px()+GranMother.Py()*GranMother.Py()+GranMother.Pz()*GranMother.Pz());
  float thetaG = acos(GranMother.Pz()/m_mgmotherp);
  m_mgmotherpt = m_mgmotherp*sin(thetaG);
  m_mgmothereta = -1*log(tan(thetaG/2));
  m_mgmotherphi = atan(GranMother.Py()/GranMother.Px());

  if (GranMother.Py() >0 && GranMother.Px() < 0) m_mgmotherphi = 
						   atan(GranMother.Py()/GranMother.Px()) + PiGreek;
  if (GranMother.Py() <0 && GranMother.Px() < 0) m_mgmotherphi = 
						   atan(GranMother.Py()/GranMother.Px()) - PiGreek;

  //selection of the Z mass particles
  //if (m_mmotherAmass < 80000 || m_mmotherAmass > 100000 || m_mmotherBmass < 80000 || m_mmotherBmass > 100000 ) ToFill = false;
  
} 

bool JpsiDumper::MCselection()
{
  bool isHiggs=false;
  int nMu_fromZ=0;       
  for(int mumc=0;mumc<mc_NMuons;mumc++)
    {
      int ID_Z = mc_MuonMotherPdgID->at(mumc);
      int ID_Higgs =  mc_MuonGrandMotherPdgID->at(mumc);
      if(ID_Z==23&&ID_Higgs==25)
	{
	nMu_fromZ++;
	}
    }

  if(nMu_fromZ>=4)isHiggs=true;
  return isHiggs; 

}