HbbAnalysis/plugins/HbbTreeMaker.cc

#include "TLorentzVector.h"
#include <limits>

#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/Common/interface/TriggerResults.h"
#include "DataFormats/Common/interface/HLTenums.h"
#include "DataFormats/Common/interface/ValueMap.h"

#include "DataFormats/Candidate/interface/Candidate.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "DataFormats/Scalers/interface/DcsStatus.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "DataFormats/BeamSpot/interface/BeamSpot.h"

#include "DataFormats/MuonReco/interface/Muon.h"
#include "DataFormats/MuonReco/interface/MuonFwd.h"
#include "DataFormats/MuonReco/interface/MuonSelectors.h"
#include "DataFormats/MuonReco/interface/MuonChamberMatch.h"
#include "DataFormats/MuonReco/interface/MuonIsolation.h"
#include "DataFormats/MuonReco/interface/MuonEnergy.h"
#include "DataFormats/MuonReco/interface/MuonTime.h"
#include "DataFormats/MuonReco/interface/MuonQuality.h"

#include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
#include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
#include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"


#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/HLTReco/interface/TriggerEvent.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/JetReco/interface/GenJet.h"

#include "DataFormats/PatCandidates/interface/Lepton.h"
#include "DataFormats/PatCandidates/interface/Muon.h"
#include "DataFormats/PatCandidates/interface/Electron.h"
#include "DataFormats/PatCandidates/interface/Tau.h"
#include "DataFormats/PatCandidates/interface/Jet.h"
#include "DataFormats/PatCandidates/interface/JetCorrFactors.h"
#include "DataFormats/PatCandidates/interface/MET.h"

#include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"


#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Framework/interface/Run.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"

#include "PhysicsTools/SelectorUtils/interface/SimpleCutBasedElectronIDSelectionFunctor.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "MagneticField/Engine/interface/MagneticField.h"

#include "UserCode/HbbAnalysis/interface/Electron.hh"
#include "UserCode/HbbAnalysis/interface/Muon.hh"
//#include "UserCode/HbbAnalysis/interface/Tau.hh" //CURRENTLY BROKEN
#include "UserCode/HbbAnalysis/interface/Jet.hh"
#include "UserCode/HbbAnalysis/interface/Met.hh"
#include "UserCode/HbbAnalysis/interface/Trigger.hh"
#include "UserCode/HbbAnalysis/interface/Vertex.hh"
//#include "UserCode/HbbAnalysis/interface/L1Object.hh" //CURRENTLY BROKEN
//#include "UserCode/HbbAnalysis/interface/HLTObject.hh" //CURRENTLY BROKEN

//AJG
#include "DataFormats/JetReco/interface/JPTJetCollection.h"
#include "DataFormats/JetReco/interface/JPTJet.h"
#include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
#include "SimDataFormats/Vertex/interface/SimVertex.h"

#include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
#include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
#include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
#include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h" 
#include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"

#include "UserCode/HbbAnalysis/plugins/HbbTreeMaker.hh"
#include "boost/regex.hpp"
#include "boost/lexical_cast.hpp"

HbbTreeMaker::HbbTreeMaker(const edm::ParameterSet & pset):
  jecUncSet_(0),
  debug_(pset.getParameter<int>("DEBUG")),
  processData_(false),
  flavour_(pset.getParameter<unsigned int>("JetFlavour")),
  doGen_(pset.getParameter<bool>("DOGEN")),
  genParticleSrc_(pset.getParameter<edm::InputTag>("GenParticles")),
  genJetSrc_(pset.getParameter<edm::InputTag>("GenJets")),
  partonJetSrc_(pset.getParameter<edm::InputTag>("PartonJets")),
  ak7genJetSrc_(pset.getParameter<edm::InputTag>("AK7GenJets")),
  ak7partonJetSrc_(pset.getParameter<edm::InputTag>("AK7PartonJets")),
  trackSrc_(pset.getParameter<edm::InputTag>("Tracks")),
  dcsSrc_(pset.getParameter<edm::InputTag>("DCSInfo")),
  electronSrc_(pset.getParameter<edm::InputTag>("Electrons")),
  muonSrc_(pset.getParameter<edm::InputTag>("Muons")),
  caloTauSrc_(pset.getParameter<edm::InputTag>("CaloTaus")),
  pfTauSrc_(pset.getParameter<edm::InputTag>("PFTaus")),
  caloJetSrc_(pset.getParameter<edm::InputTag>("CaloJets")),
  jptJetSrc_(pset.getParameter<edm::InputTag>("JPTJets")),
  ak7JetSrc_(pset.getParameter<edm::InputTag>("AK7Jets")),
  pfJetSrc_(pset.getParameter<edm::InputTag>("PFJets")),
  ak7pfJetSrc_(pset.getParameter<edm::InputTag>("AK7PFJets")),
  caloMetSrc_(pset.getParameter<edm::InputTag>("CaloMET")),
  tcMetSrc_(pset.getParameter<edm::InputTag>("TCMET")),
  pfMetSrc_(pset.getParameter<edm::InputTag>("PFMET")),
  //pairSrc_(pset.getParameter<edm::InputTag>("Pair")),
  //mmPairSrc_(pset.getParameter<edm::InputTag>("MuMuPair")),
  //etPairSrc_(pset.getParameter<edm::InputTag>("ETauPair")),
  //mtPairSrc_(pset.getParameter<edm::InputTag>("MuTauPair")),
  vertexSrc_(pset.getParameter<edm::InputTag>("Vertex")),
  triggerSrc_(pset.getParameter<edm::InputTag>("Trigger")),
  hltPaths_(pset.getParameter<std::vector<std::string> >("HLTPaths")),
  l1CenJetSrc_(pset.getParameter<edm::InputTag>("L1CenJet")),
  l1TauJetSrc_(pset.getParameter<edm::InputTag>("L1TauJet")),
  l1FwdJetSrc_(pset.getParameter<edm::InputTag>("L1FwdJet")),
  hltSummarySrc_(pset.getParameter<edm::InputTag>("HLTSummary")),
  status1ToKeep_(pset.getParameter<std::vector<std::string> >("Status1ToKeep")),
  status2ToKeep_(pset.getParameter<std::vector<std::string> >("Status2ToKeep")),
  status3ToKeep_(pset.getParameter<std::vector<std::string> >("Status3ToKeep")),
  event_(0),
  tree_(0)
  //processVec_(pset.getParameter<std::vector<std::string> >("ProcessVec"))
{//constructor
  for (unsigned i = 0; i < status1ToKeep_.size(); ++i){
    status1ToKeepRegex_.push_back(boost::regex(status1ToKeep_[i]));
  }
  for (unsigned i = 0; i < status2ToKeep_.size(); ++i){
    status2ToKeepRegex_.push_back(boost::regex(status2ToKeep_[i]));
  }
  for (unsigned i = 0; i < status3ToKeep_.size(); ++i){
    status3ToKeepRegex_.push_back(boost::regex(status3ToKeep_[i]));
  }
}//constructor

HbbTreeMaker::~HbbTreeMaker()
{//destructor
  if (jecUncSet_){
  delete jecUncCalo;
  delete jecUncPF;
  delete jecUncJPT;
  }

}//destructor


void HbbTreeMaker::beginJob(){//beginJob


  edm::Service<TFileService> lFileService;

  TFileDirectory lDir = lFileService->mkdir("HbbAnalysis");

  tree_ = lDir.make<TTree>("Tree","Tree");
  tree_->Branch("HbbEvent","HbbAnalysis::HbbEvent",&event_);
  event_ = new HbbAnalysis::HbbEvent();
  
  if (debug_) std::cout << "Initialising JetFlavour : " << std::endl;

  lDir = lFileService->mkdir("JetFlavours");
  jetFlav_.Initialise(lDir, debug_, flavour_);


}//beginJob


void HbbTreeMaker::endJob(){//endJob
  if (!processData_) jetFlav_.printSummary();

  //tree_->Write();
  //delete tree_;
  //delete event_;
  
}//endJob

void HbbTreeMaker::beginRun(const edm::Run& aRun, const edm::EventSetup& aEvtSetup){//beginRun
    bool changed(true);
      if (hltConfig_.init(aRun,aEvtSetup,"HLT",changed)) {
            // if init returns TRUE, initialisation has succeeded!
            if (changed) {
            // The HLT config has actually changed wrt the previous Run, hence rebook your
            // histograms or do anything else dependent on the revised HLT config
            }
            } else {
            // if init returns FALSE, initialisation has NOT succeeded, which indicates a problem
            // with the file and/or code and needs to be investigated!
              std::cerr << " HLT config extraction failure with process name " << "HLT";
            // In this case, all access methods will return empty values!
            }

}

void HbbTreeMaker::analyze(const edm::Event& aEvt, const edm::EventSetup& aEvtSetup){//analyze

  //dumpContent(aEvt,"","MuTauPAT");
  if (debug_) std::cout << "Processing event #" << aEvt.id().event() << std::endl;


  static bool firstEvent = true;
  if (!jecUncSet_){
    edm::ESHandle<JetCorrectorParametersCollection> JetCorParCollCalo;
    edm::ESHandle<JetCorrectorParametersCollection> JetCorParCollPF;
    edm::ESHandle<JetCorrectorParametersCollection> JetCorParCollJPT;
    aEvtSetup.get<JetCorrectionsRecord>().get("AK5Calo",JetCorParCollCalo); 
    aEvtSetup.get<JetCorrectionsRecord>().get("AK5PF",JetCorParCollPF); 
    aEvtSetup.get<JetCorrectionsRecord>().get("AK5JPT",JetCorParCollJPT); 
    JetCorrectorParameters const & JetCorParCalo = (*JetCorParCollCalo)["Uncertainty"];
    JetCorrectorParameters const & JetCorParPF = (*JetCorParCollPF)["Uncertainty"];
    JetCorrectorParameters const & JetCorParJPT = (*JetCorParCollJPT)["Uncertainty"];

    jecUncCalo = new JetCorrectionUncertainty(JetCorParCalo);
    jecUncPF = new JetCorrectionUncertainty(JetCorParPF);
    jecUncJPT = new JetCorrectionUncertainty(JetCorParJPT);
    jecUncSet_ = true;
  }

  event_->Clear();
  event_->event(aEvt.id().event());
  processData_ = aEvt.isRealData();

  event_->run(aEvt.run());
  event_->lumiBlock(aEvt.luminosityBlock());
  event_->isRealData(aEvt.isRealData());
  event_->bx(aEvt.bunchCrossing());

  edm::Handle<reco::BeamSpot> lBeamSpot;
  try {
    aEvt.getByLabel(edm::InputTag("offlineBeamSpot"), lBeamSpot);
  } catch(cms::Exception& e) {
    std::cout<< "AMM: No beam spot found. Exception : " << e.what() << "." << std::endl;
  }
  HbbBeamSpot(lBeamSpot,event_->beamSpot());

  //AJG - This will only work on RECO sample
  /*
  std::cout << "Beamspot = (" << lBeamSpot->x0() << "," << lBeamSpot->y0() << "," << lBeamSpot->z0() << ")" << std::endl;
  
  edm::Handle<edm::SimVertexContainer> simVertexCollection;
  aEvt.getByLabel("g4SimHits", simVertexCollection);
  const edm::SimVertexContainer simVC = *(simVertexCollection.product());
  std::cout << "SimVertexContainer size = " << simVC.size() << std::endl;
  std::cout << "SimVertex[0] = (" << simVC.at(0).position().x() << ","
                                  << simVC.at(0).position().y() << ","
                                  << simVC.at(0).position().z() << ")" << std::endl;
  */

  //used in jet method....
  edm::Handle<reco::GenParticleCollection> lGenParticles;

  if (!processData_) {
    edm::Handle<edm::HepMCProduct> mcProd;
    edm::Handle<GenRunInfoProduct> lRunProd;
    edm::Handle<GenEventInfoProduct> genEventInfoHandle;

    try {
      aEvt.getByLabel("generator", mcProd );
    } catch(cms::Exception& e)  {
      std::cout << "AMM: Collection HepMCProduct not available! Exception : " << e.what() << ". " << std::endl;
    }
    try {
      aEvt.getRun().getByLabel("generator",lRunProd);
    } catch(cms::Exception& e)  {
      std::cout << "AMM: Collection GenInfoProduct not available! Exception : " << e.what() << ". " << std::endl;
    }
    try {
      aEvt.getByLabel("generator",genEventInfoHandle);
    } catch(cms::Exception& e)  {
      std::cout << "AMM: Collection GenEventInfoProduct not available! Exception : " << e.what() << ". " << std::endl;
    }

    HbbGenInfo(mcProd,lRunProd,genEventInfoHandle);

    edm::Handle<reco::GenJetCollection> lGenJets;
    try {
      aEvt.getByLabel(genJetSrc_,lGenJets);  
      if (debug_) std::cout << "** ngenJets = " << lGenJets->size() << std::endl;
    } catch(cms::Exception& e)  {
      std::cout << "AMM: Collection GenJets not available! Exception : " << e.what() << ". " << std::endl;
    }

    //do not save components
    if (lGenJets.isValid()) HbbGenJets(lGenJets,event_->genjets(),false);

    edm::Handle<reco::GenJetCollection> lPartonJets;
    try {
      aEvt.getByLabel(partonJetSrc_,lPartonJets);  
      if (debug_) std::cout << "** npartonJets = " << lPartonJets->size() << std::endl;
    } catch(cms::Exception& e)  {
      std::cout << "AMM: Collection PartonJets not available! Exception : " << e.what() << ". " << std::endl;
    }

    //save components
    if (lPartonJets.isValid()) HbbGenJets(lPartonJets,event_->partonjets(),true);
    
    //AK7 
    edm::Handle<reco::GenJetCollection> lAk7GenJets;
    try {
      aEvt.getByLabel(ak7genJetSrc_,lAk7GenJets);  
      if (debug_) std::cout << "** ngenJets = " << lAk7GenJets->size() << std::endl;
    } catch(cms::Exception& e)  {
      std::cout << "AMM: Collection GenJets not available! Exception : " << e.what() << ". " << std::endl;
    }

    //do not save components
    if (lAk7GenJets.isValid()) HbbGenJets(lAk7GenJets,event_->ak7genjets(),false);

    edm::Handle<reco::GenJetCollection> lAk7PartonJets;
    try {
      aEvt.getByLabel(ak7partonJetSrc_,lAk7PartonJets);  
      if (debug_) std::cout << "** npartonJets = " << lAk7PartonJets->size() << std::endl;
    } catch(cms::Exception& e)  {
      std::cout << "AMM: Collection PartonJets not available! Exception : " << e.what() << ". " << std::endl;
    }

    //save components
    if (lAk7PartonJets.isValid()) HbbGenJets(lAk7PartonJets,event_->ak7partonjets(),true);

    edm::Handle<LHEEventProduct> lLHEEvt;
    try {
      aEvt.getByType(lLHEEvt);
    } catch(cms::Exception& e)  {
      std::cout << "AMM: Collection of LHEEventProduct not available! Exception : " << e.what() << ". " << std::endl;
    }

    if (lLHEEvt.isValid()){
      HbbLHEInfo(lLHEEvt,event_->lheParticles());
    }
    
    edm::Handle<std::vector<PileupSummaryInfo> > lPUInfo;
    try {
      aEvt.getByLabel("addPileupInfo",lPUInfo);
    } catch(cms::Exception& e)  {
      std::cout << "AMM: Collection of PileupSummaryInfo  not available! Exception : " << e.what() << ". " << std::endl;
    }

    if (lPUInfo.isValid()){
      HbbPUInfo(lPUInfo,event_->puVars());
    }

    try {
      aEvt.getByLabel(genParticleSrc_,lGenParticles);  
      if (debug_) std::cout << "** ngenParticles = " << lGenParticles->size() << std::endl;
    } catch(cms::Exception& e)  {
      std::cout << "AMM: Collection genParticles not available! Exception : " << e.what() << ". " << std::endl;
    }

    //if (doGen_) HbbParticles(lGenParticles,event_->particles());
    HbbParticles(lGenParticles,event_->particles(), event_->p4total());

    unsigned int lNPartons = 0;
    lNPartons = jetFlav_.fillPartons(lGenParticles);
  
    if (debug_) std::cout << "--- Number of partons = " << lNPartons << "." << std::endl;
  }

  edm::Handle<std::vector<reco::Vertex> > lRecoVertices;
  try {
    aEvt.getByLabel(vertexSrc_, lRecoVertices);
    if (!lRecoVertices.isValid()){
      edm::LogInfo("ERROR")<< "Error! Can't get vertex by label. ";
    }
    if (debug_) std::cout << "** vertexcollection = " << lRecoVertices->size() << " elements." << std::endl;
  } catch(cms::Exception& e)  {
    std::cout << "AMM: Collection " << vertexSrc_  << " not available! Exception : " << e.what() << ". " << std::endl;
  }

  HbbVertices(lRecoVertices,event_->vertices());

  edm::Handle<std::vector<pat::Electron> > lElectronCollection;
  
  try {
    aEvt.getByLabel(electronSrc_,lElectronCollection);
    if (!lElectronCollection.isValid()){
      edm::LogInfo("ERROR")<< "Error! Can't get electron by label. ";
    } 
    if (debug_) std::cout << "** electroncollection = " << lElectronCollection->size() << " elements." << std::endl;
  } catch(cms::Exception& e)  {
    std::cout << "AMM: Collection " << electronSrc_  << " not available! Exception : " << e.what() << ". " << std::endl;
  }

  HbbElectrons(lElectronCollection,event_->electrons(),lBeamSpot);

  edm::Handle<std::vector<pat::Muon> > lMuonCollection;
  
  try {
    aEvt.getByLabel(muonSrc_,lMuonCollection);
    if (!lMuonCollection.isValid()){
      edm::LogInfo("ERROR")<< "Error! Can't get muon by label. ";
    } 
    if (debug_) std::cout << "** muoncollection = " << lMuonCollection->size() << " elements." << std::endl;
  } catch(cms::Exception& e)  {
    std::cout << "AMM: Collection " << muonSrc_  << " not available! Exception : " << e.what() << ". " << std::endl;
  }

  HbbMuons(lMuonCollection,lRecoVertices,lBeamSpot,event_->muons());


  /*
  if (!( caloTauSrc_.label()=="" && caloTauSrc_.instance()=="" )) {
    edm::Handle<std::vector<pat::Tau> > lTauCollection;
  
    try {
      aEvt.getByLabel(caloTauSrc_,lTauCollection);
      if (!lTauCollection.isValid()){
        edm::LogInfo("ERROR")<< "Error! Can't get caloTau by label. ";
      } 
      if (debug_) std::cout << "** caloTaucollection = " << lTauCollection->size() << " elements." << std::endl;
    } catch(cms::Exception& e)  {
      std::cout << "AMM: Collection " << caloTauSrc_  << " not available! Exception : " << e.what() << ". " << std::endl;
    }

    HbbTaus(lTauCollection,lRecoVertices,event_->caloTaus());
  }

  edm::Handle<std::vector<pat::Tau> > lPFTauCollection;
  
  try {
    aEvt.getByLabel(pfTauSrc_,lPFTauCollection);
    if (!lPFTauCollection.isValid()){
      edm::LogInfo("ERROR")<< "Error! Can't get pfTau by label. ";
    } 
    if (debug_) std::cout << "** pfTaucollection = " << lPFTauCollection->size() << " elements." << std::endl;
  } catch(cms::Exception& e)  {
    std::cout << "AMM: Collection " << pfTauSrc_  << " not available! Exception : " << e.what() << ". " << std::endl;
  }

  HbbTaus(lPFTauCollection,lRecoVertices,event_->pfTaus());
  */

  edm::Handle<std::vector<pat::Jet> > lCaloJetCollection;
  
  try {
    aEvt.getByLabel(caloJetSrc_,lCaloJetCollection);
    if (!lCaloJetCollection.isValid()){
      edm::LogInfo("ERROR")<< "Error! Can't get caloJets by label. ";
    } 
    if (debug_) std::cout << "** caloJetcollection = " << lCaloJetCollection->size() << " elements." << std::endl;
  } catch(cms::Exception& e)  {
    std::cout << "AMM: Collection " << caloJetSrc_  << " not available! Exception : " << e.what() << ". " << std::endl;
  }

  HbbJets(lCaloJetCollection,jetFlav_,lGenParticles,event_->caloJets());
 
  //AJG - Process reco JPT Jet Corrections:
  //This section should be uncommented for testing purposes

  /* 
     edm::Handle<reco::JPTJetCollection > lJptRL2L3RESJetCollection;
     edm::Handle<reco::JPTJetCollection > lJptRL1L2L3RESJetCollection;

     edm::Handle<reco::JPTJetCollection > lJptRRAWJetCollection;
     try {
     aEvt.getByLabel("JetPlusTrackZSPCorJetAntiKt5",lJptRRAWJetCollection);
     if (!lJptRRAWJetCollection.isValid()){
     edm::LogInfo("ERROR")<< "Error! Can't get JetPlusZSPCorJetAntiKt5 by label. ";
     } 
     if (debug_) std::cout << "** jetPlusTrackZSPCorJetAntiKt5 = " << lJptRRAWJetCollection->size() << " elements." << std::endl;
     } catch(cms::Exception& e)  {
     std::cout << "AJG: Collection not available! Exception : " << e.what() << ". " << std::endl;
     }

     try {
     aEvt.getByLabel("ak5JPTJetsL2L3Residual",lJptRL2L3RESJetCollection);
     aEvt.getByLabel("JetPlusTrackZSPCorJetAntiKt5",lJptRRAWJetCollection);
     aEvt.getByLabel("ak5JPTJetsL1L2L3Residual",lJptRL1L2L3RESJetCollection);
     if (!lJptRRAWJetCollection.isValid()){
     edm::LogInfo("ERROR")<< "Error! Can't get jptRecoJets by label. ";
     } 
     if (debug_) std::cout << "** jptRecoJetcollection = " << lJptRRAWJetCollection->size() << " elements." << std::endl;
     } catch(cms::Excepftion& e)  {
     std::cout << "AJG: Collection " << jptJetSrc_  << " not available! Exception : " << e.what() << ". " << std::endl;
     }


     for(reco::JPTJetCollection::const_iterator jptjet = lJptRRAWJetCollection->begin(); jptjet != lJptRRAWJetCollection->end(); ++jptjet ) { 

     std::cout << "pT = " << jptjet->pt() << std::endl;

     reco::TrackRefVector tpionsInVInC = jptjet->getPionsInVertexInCalo();
     std::cout << "Pion Track Count (reco) = " << tpionsInVInC.size() << std::endl;
     }

     reco::JPTJetCollection::const_iterator jptRAWjet = lJptRRAWJetCollection->begin();
     reco::JPTJetCollection::const_iterator jptL2L3RESjet = lJptRL2L3RESJetCollection->begin();
     reco::JPTJetCollection::const_iterator jptL1L2L3RESjet = lJptRL1L2L3RESJetCollection->begin();
     std::cout << "pT RAW = " << jptRAWjet->pt() << std::endl;
     std::cout << "pT L2L3RES = " << jptL2L3RESjet->pt() << std::endl;
     std::cout << "pT L1L2L3RES = " << jptL1L2L3RESjet->pt() << std::endl;
     */


  edm::Handle<std::vector<pat::Jet> > lJptJetCollection;
  
  try {
    aEvt.getByLabel(jptJetSrc_,lJptJetCollection);
    if (!lJptJetCollection.isValid()){
      edm::LogInfo("ERROR")<< "Error! Can't get jptJets by label. ";
    } 
    if (debug_) std::cout << "** jptJetcollection = " << lJptJetCollection->size() << " elements." << std::endl;
  } catch(cms::Exception& e)  {
    std::cout << "AMM: Collection " << jptJetSrc_  << " not available! Exception : " << e.what() << ". " << std::endl;
  }

  //std::cout << "Processing JPT jets:" << std::endl;
  HbbJets(lJptJetCollection,jetFlav_,lGenParticles,event_->jptJets());
  
  edm::Handle<std::vector<pat::Jet> > lAk7JetCollection;
  
  try {
    aEvt.getByLabel(ak7JetSrc_,lAk7JetCollection);
    if (!lAk7JetCollection.isValid()){
      edm::LogInfo("ERROR")<< "Error! Can't get ak7Jets by label. ";
    } 
    if (debug_) std::cout << "** ak7Jetcollection = " << lAk7JetCollection->size() << " elements." << std::endl;
  } catch(cms::Exception& e)  {
    std::cout << "AMM: Collection " << ak7JetSrc_  << " not available! Exception : " << e.what() << ". " << std::endl;
  }

  //std::cout << "Processing JPT jets:" << std::endl;
  HbbJets(lAk7JetCollection,jetFlav_,lGenParticles,event_->ak7Jets());
  
  edm::Handle<std::vector<pat::Jet> > lPfJetCollection;
  
  try {
    aEvt.getByLabel(pfJetSrc_,lPfJetCollection);
    if (!lPfJetCollection.isValid()){
      edm::LogInfo("ERROR")<< "Error! Can't get pfJets by label. ";
    } 
    if (debug_) std::cout << "** pfJetcollection = " << lPfJetCollection->size() << " elements." << std::endl;
  } catch(cms::Exception& e)  {
    std::cout << "AMM: Collection " << pfJetSrc_  << " not available! Exception : " << e.what() << ". " << std::endl;
  }

  //std::cout << "Processing PF jets:" << std::endl;
  HbbJets(lPfJetCollection,jetFlav_,lGenParticles,event_->pfJets());
  
  edm::Handle<std::vector<pat::Jet> > lAk7PfJetCollection;
  
  try {
    aEvt.getByLabel(ak7pfJetSrc_,lAk7PfJetCollection);
    if (!lAk7PfJetCollection.isValid()){
      edm::LogInfo("ERROR")<< "Error! Can't get ak7pfJets by label. ";
    } 
    if (debug_) std::cout << "** ak7pfJetcollection = " << lAk7PfJetCollection->size() << " elements." << std::endl;
  } catch(cms::Exception& e)  {
    std::cout << "AMM: Collection " << ak7pfJetSrc_  << " not available! Exception : " << e.what() << ". " << std::endl;
  }

  //std::cout << "Processing PF jets:" << std::endl;
  HbbJets(lAk7PfJetCollection,jetFlav_,lGenParticles,event_->ak7pfJets());


  edm::Handle<std::vector<pat::MET> > lCaloMetCol;

  try {
    aEvt.getByLabel(caloMetSrc_,lCaloMetCol);
    if (!lCaloMetCol.isValid()){
      edm::LogInfo("ERROR")<< "Error! Can't get caloMet by label. ";
    } 
    if (debug_) std::cout << "** caloMet  = " << lCaloMetCol->size() << " elements." << std::endl;
  } catch(cms::Exception& e)  {
    std::cout << "AMM: Collection " << caloMetSrc_  << " not available! Exception : " << e.what() << ". " << std::endl;
  }

  HbbMet(lCaloMetCol,event_->caloMet());
  
  edm::Handle<std::vector<pat::MET> > lTcMetCol;

  try {
    aEvt.getByLabel(tcMetSrc_,lTcMetCol);
    if (!lTcMetCol.isValid()){
      edm::LogInfo("ERROR")<< "Error! Can't get tcMet by label. ";
    } 
    if (debug_) std::cout << "** tcMet  = " << lTcMetCol->size() << " elements." << std::endl;
  } catch(cms::Exception& e)  {
    std::cout << "AMM: Collection " << tcMetSrc_  << " not available! Exception : " << e.what() << ". " << std::endl;
  }

  HbbMet(lTcMetCol,event_->tcMet());
  
  edm::Handle<std::vector<pat::MET> > lPfMetCol;

  try {
    aEvt.getByLabel(pfMetSrc_,lPfMetCol);
    if (!lPfMetCol.isValid()){
      edm::LogInfo("ERROR")<< "Error! Can't get pfMet by label. ";
    } 
    if (debug_) std::cout << "** pfMet  = " << lPfMetCol->size() << " elements." << std::endl;
  } catch(cms::Exception& e)  {
    std::cout << "AMM: Collection " << pfMetSrc_  << " not available! Exception : " << e.what() << ". " << std::endl;
  }

  HbbMet(lPfMetCol,event_->pfMet());


  //triggers
  edm::Handle<edm::TriggerResults> lTrigCol;
  try {
    aEvt.getByLabel(triggerSrc_,lTrigCol);
    if (!lTrigCol.isValid()){
      edm::LogInfo("ERROR")<< "Error! Can't get triggers by label. ";
    } 
    if (debug_) std::cout << "** triggers = " << lTrigCol->size() << std::endl;
  } catch(cms::Exception& e)  {
    std::cout << "AMM: Collection " << triggerSrc_  << " not available! Exception : " << e.what() << ". " << std::endl;
  }
  if (lTrigCol.isValid()) {
    const edm::TriggerNames & lNames = aEvt.triggerNames(*lTrigCol);
    int prescaleSet = hltConfig_.prescaleSet(aEvt, aEvtSetup);
    HbbTrigger(lTrigCol,lNames,event_->triggers(), prescaleSet);
  }

  //L1 jet objects

  /*
  edm::Handle<l1extra::L1JetParticleCollection> lCenJet;
  edm::Handle<l1extra::L1JetParticleCollection> lFwdJet;
  edm::Handle<l1extra::L1JetParticleCollection> lTauJet;
  try {
    aEvt.getByLabel(l1CenJetSrc_,lCenJet);
    aEvt.getByLabel(l1TauJetSrc_,lTauJet);
    aEvt.getByLabel(l1FwdJetSrc_,lFwdJet);
  } catch(cms::Exception& e)  {
    std::cout << "AMM: L1 Collections " << l1CenJetSrc_  << ", " << l1TauJetSrc_ << ", " << l1FwdJetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
  }
  if (lCenJet.isValid()) {
    HbbL1Objects(lCenJet,event_->l1Jets(),1);
  }
  if (lTauJet.isValid()) {
    HbbL1Objects(lTauJet,event_->l1Jets(),2);
  }
  if (lFwdJet.isValid()) {
    HbbL1Objects(lFwdJet,event_->l1Jets(),3);
  }


  edm::Handle<trigger::TriggerEvent> lHltSummary;
  try {
    aEvt.getByLabel(hltSummarySrc_,lHltSummary);
  } catch(cms::Exception& e)  {
    std::cout << "AMM: collection " << hltSummarySrc_  << " not available! Exception : " << e.what() << ". " << std::endl;
  }
  if (lHltSummary.isValid()) {
    HbbHLTObjects(lHltSummary,event_->hltObjects());
  }
  */


  //event_->order();
  tree_->Fill();
  firstEvent = false;

}//analyze


void HbbTreeMaker::HbbBeamSpot(const edm::Handle<reco::BeamSpot> & aBeamSpot,
                               HbbAnalysis::BeamSpot & aBS)
{

  if (aBeamSpot.isValid()) {
    HbbAnalysis::BeamSpotVars lbs;
    lbs.x0 = aBeamSpot->x0();
    lbs.y0 = aBeamSpot->y0();
    lbs.z0 = aBeamSpot->z0();
    lbs.sigmaZ = aBeamSpot->sigmaZ();
    lbs.BeamWidthX = aBeamSpot->BeamWidthX();
    lbs.BeamWidthY = aBeamSpot->BeamWidthY();
    lbs.x0Error = aBeamSpot->x0Error();
    lbs.y0Error = aBeamSpot->y0Error();
    lbs.z0Error = aBeamSpot->z0Error();
    lbs.sigmaZ0Error = aBeamSpot->sigmaZ0Error();
    lbs.BeamWidthXError = aBeamSpot->BeamWidthXError();
    lbs.BeamWidthYError = aBeamSpot->BeamWidthYError();

    aBS.bsVars(lbs);
  }

}


void HbbTreeMaker::HbbElectrons(const edm::Handle<std::vector<pat::Electron> > & aCol,
                                std::vector<HbbAnalysis::Electron> & aVec, const edm::Handle<reco::BeamSpot> & aBeamSpot)
{//HbbElectrons

  if (aCol.isValid()){
    if (aCol->size() > 0) {
        
      unsigned int iEle = 0;
      for (std::vector<pat::Electron>::const_iterator iter = aCol->begin(); 
           iter != aCol->end(); 
           ++iter,++iEle) 
        {
          //if (debug_ > 1) std::cout << "**** Ele #" << iEle << ", pT,eta,phi =  " << (*iter).pt() << " " << (*iter).eta() << " " << (*iter).phi() << std::endl;
// get the track collection


          // define your function for  WPXX using relative or absolute isolations
          SimpleCutBasedElectronIDSelectionFunctor patSeleR95(SimpleCutBasedElectronIDSelectionFunctor::relIso95);
          SimpleCutBasedElectronIDSelectionFunctor patSeleR90(SimpleCutBasedElectronIDSelectionFunctor::relIso90);
          SimpleCutBasedElectronIDSelectionFunctor patSeleR85(SimpleCutBasedElectronIDSelectionFunctor::relIso85);
          SimpleCutBasedElectronIDSelectionFunctor patSeleR80(SimpleCutBasedElectronIDSelectionFunctor::relIso80);
          SimpleCutBasedElectronIDSelectionFunctor patSeleR70(SimpleCutBasedElectronIDSelectionFunctor::relIso70);
          SimpleCutBasedElectronIDSelectionFunctor patSeleR60(SimpleCutBasedElectronIDSelectionFunctor::relIso60);
          SimpleCutBasedElectronIDSelectionFunctor patSeleC95(SimpleCutBasedElectronIDSelectionFunctor::cIso95);
          SimpleCutBasedElectronIDSelectionFunctor patSeleC90(SimpleCutBasedElectronIDSelectionFunctor::cIso90);
          SimpleCutBasedElectronIDSelectionFunctor patSeleC85(SimpleCutBasedElectronIDSelectionFunctor::cIso85);
          SimpleCutBasedElectronIDSelectionFunctor patSeleC80(SimpleCutBasedElectronIDSelectionFunctor::cIso80);
          SimpleCutBasedElectronIDSelectionFunctor patSeleC70(SimpleCutBasedElectronIDSelectionFunctor::cIso70);
          SimpleCutBasedElectronIDSelectionFunctor patSeleC60(SimpleCutBasedElectronIDSelectionFunctor::cIso60);
   

    //----------------------------------------------------------------------
    HbbAnalysis::FourMomentum lFourVec;
    //----------------------------------------------------------------------
          lFourVec.e = (*iter).energy();
          lFourVec.px = (*iter).px();
          lFourVec.py = (*iter).py();
          lFourVec.pz = (*iter).pz();
    //----------------------------------------------------------------------
    HbbAnalysis::LeptonVars lLep;
    //----------------------------------------------------------------------
          lLep.vertexPCA.vx = (*iter).vx();
          lLep.vertexPCA.vy = (*iter).vy();
          lLep.vertexPCA.vz = (*iter).vz();
          lLep.charge = (*iter).charge();
    lLep.dBeamspot= (*iter).gsfTrack()->dxy(*aBeamSpot);
    //lLep.dBeamspot = (*iter).dB(pat::Electron::BS2D);
    lLep.dVertex = (*iter).dB(pat::Electron::PV2D);
    if ((*iter).triggerObjectMatches().size() > 0){
      lLep.hltMatchPaths = (*iter).triggerObjectMatches().at(0).pathNames(1,0);
      lLep.hltMatchFilters = (*iter).triggerObjectMatches().at(0).filterLabels();
    }
    //Construct a GenParticle if processing MC and a match has been found by PAT
          if (!processData_ && (*iter).genLepton()){
      HbbAnalysis::FourMomentum lGenFourVec;
            HbbAnalysis::GenVars lGen;
            lGenFourVec.px = (*iter).genLepton()->px();
            lGenFourVec.py = (*iter).genLepton()->py();
            lGenFourVec.pz = (*iter).genLepton()->pz();
            lGenFourVec.e = (*iter).genLepton()->energy();
            lGen.charge = (*iter).genLepton()->charge();
            lGen.pdgId = (*iter).genLepton()->pdgId();
            lGen.status = (*iter).genLepton()->status();
            lGen.vertex.vx = (*iter).genLepton()->vx();
            lGen.vertex.vy = (*iter).genLepton()->vy();
            lGen.vertex.vz = (*iter).genLepton()->vz();
      lLep.genMatches.push_back(HbbAnalysis::GenParticle(lGenFourVec, lGen));
          }

    
    //TESTING CODE
    /* 
    std::cout << "-----------------------------------------------" << std::endl;
    std::cout << "Electron: " << iEle << std::endl;
    std::cout << "--Found " << (*iter).triggerObjectMatches().size() << " matched trigger objects" << std::endl;
    std::cout << "--Path Names with Option (1,1): " << std::endl;
    for (unsigned i = 0; i < (*iter).triggerObjectMatches().size(); ++i){
      std::cout << "---Trigger object has pT: " << (*iter).triggerObjectMatches().at(i).pt() << std::endl;
      std::vector<std::string> HLTpathNames = (*iter).triggerObjectMatches().at(i).pathNames(1,1);
      std::cout << "---HLT Path Names:" << std::endl;
      for (unsigned j = 0; j < HLTpathNames.size(); ++j){
        std::cout << "----" << HLTpathNames.at(j) << std::endl;
      }
    }
    std::cout << "--Path Names with Option (1,0): " << std::endl;
    for (unsigned i = 0; i < (*iter).triggerObjectMatches().size(); ++i){
      std::vector<std::string> HLTpathNames = (*iter).triggerObjectMatches().at(i).pathNames(1,0);
      std::cout << "---HLT Path Names:" << std::endl;
      for (unsigned j = 0; j < HLTpathNames.size(); ++j){
        std::cout << "----" << HLTpathNames.at(j) << std::endl;
      }
    }
    std::cout << "--Path Names with Option (0,1): " << std::endl;
    for (unsigned i = 0; i < (*iter).triggerObjectMatches().size(); ++i){
      std::vector<std::string> HLTpathNames = (*iter).triggerObjectMatches().at(i).pathNames(0,1);
      std::cout << "---HLT Path Names:" << std::endl;
      for (unsigned j = 0; j < HLTpathNames.size(); ++j){
        std::cout << "----" << HLTpathNames.at(j) << std::endl;
      }
    }
    std::cout << "--Path Names with Option (0,0): " << std::endl;
    for (unsigned i = 0; i < (*iter).triggerObjectMatches().size(); ++i){
      std::vector<std::string> HLTpathNames = (*iter).triggerObjectMatches().at(i).pathNames(0,0);
      std::cout << "---HLT Path Names:" << std::endl;
      for (unsigned j = 0; j < HLTpathNames.size(); ++j){
        std::cout << "----" << HLTpathNames.at(j) << std::endl;
      }
    }
    for (unsigned i = 0; i < (*iter).triggerObjectMatches().size(); ++i){
      std::vector<std::string> HLTfilterLabels = (*iter).triggerObjectMatches().at(i).filterLabels();
      std::cout << "---HLT FilterLabels:" << std::endl;
      for (unsigned j = 0; j < HLTfilterLabels.size(); ++j){
        std::cout << "----" << HLTfilterLabels.at(j) << std::endl;
      }
    }
    */
    
    
    HbbAnalysis::SCVars lSC;
          if ( (*iter).superCluster().isAvailable() && (*iter).superCluster().isNonnull() ) 
            lSC.eta = (*iter).superCluster()->eta();
          else lSC.eta = (*iter).eta();
          lSC.sigmaEtaEta = (*iter).scSigmaEtaEta();
          lSC.sigmaIEtaIEta = (*iter).scSigmaIEtaIEta();
          lSC.e1x5 = (*iter).scE1x5();
          lSC.e2x5Max = (*iter).scE2x5Max();
          lSC.e5x5 = (*iter).scE5x5();
          lSC.eOverP = (*iter).eSuperClusterOverP();

          HbbAnalysis::EleIsoVars lIso;
          lIso.track = (*iter).trackIso();
          lIso.ecal = (*iter).ecalIso();
          lIso.hcal = (*iter).hcalIso();

          HbbAnalysis::EleIDVars lID;
          lID.idAndIso = 
            (patSeleR95(*iter)<<11) | (patSeleR90(*iter)<<10) |
            (patSeleR85(*iter)<<9) | (patSeleR80(*iter)<<8) |
            (patSeleR70(*iter)<<7) | (patSeleR60(*iter)<<6) |
            (patSeleC95(*iter)<<5) | (patSeleC90(*iter)<<4) |
            (patSeleC85(*iter)<<3) | (patSeleC80(*iter)<<2) |
            (patSeleC70(*iter)<<1) | (patSeleC60(*iter));
    std::vector<std::pair<std::string,float> > electronIDsTemp = (*iter).electronIDs();
    for (unsigned iIDs = 0; iIDs < electronIDsTemp.size(); ++iIDs){
      lID.electronIDs[electronIDsTemp.at(iIDs).first] = electronIDsTemp.at(iIDs).second;
    }
          lID.hOverE = (*iter).hadronicOverEm();
          lID.deltaPhiIn = (*iter).deltaPhiSuperClusterTrackAtVtx();
          lID.deltaEtaIn = (*iter).deltaEtaSuperClusterTrackAtVtx();
          lID.ecalDrivenSeed = (*iter).ecalDrivenSeed();
          lID.trackerDrivenSeed = (*iter).trackerDrivenSeed();

          aVec.push_back(HbbAnalysis::Electron(lFourVec,lLep,lSC,lIso,lID));

    //Testing
    std::vector<HbbAnalysis::Electron>::iterator it = aVec.end() - 1;
    if ( (int(it->idVars().electronIDs.find("simpleEleId85relIso")->second+0.5) != 7 &&  it->WPrel85())
        || (int(it->idVars().electronIDs.find("simpleEleId85relIso")->second+0.5) == 7 &&  !it->WPrel85()) ) {
      std::cout << "SHIT!" << std::endl;
      std::cout << "Direct: " << it->WPrel85() << "\tFrom PAT: "  << int(it->idVars().electronIDs.find("simpleEleId85relIso")->second+0.5) << std::endl;
    }

    
        }
    }
  }

}//HbbElectrons


void HbbTreeMaker::HbbMuons(const edm::Handle<std::vector<pat::Muon> > & aCol,
                            const edm::Handle<std::vector<reco::Vertex> > & aRecoVertices,
                            const edm::Handle<reco::BeamSpot> & aBeamSpot,
                            std::vector<HbbAnalysis::Muon> & aVec)
{//HbbMuons

  if (aCol.isValid()){
    if (aCol->size() > 0) {
        
      unsigned int iEle = 0;
      for (std::vector<pat::Muon>::const_iterator iter = aCol->begin(); 
           iter != aCol->end(); 
           ++iter,++iEle) 
        {
          const reco::Muon* recoMuon = dynamic_cast<const reco::Muon*>((*iter).originalObject());

          //select only global and tracker muons
          if (!recoMuon->isGlobalMuon() || !recoMuon->isTrackerMuon()) continue;

    //----------------------------------------------------------------------
    HbbAnalysis::FourMomentum lFourVec;
    //----------------------------------------------------------------------
          lFourVec.e = (*iter).energy();
          lFourVec.px = (*iter).px();
          lFourVec.py = (*iter).py();
          lFourVec.pz = (*iter).pz();
    //----------------------------------------------------------------------
    HbbAnalysis::LeptonVars lLep;
    //----------------------------------------------------------------------
          lLep.vertexPCA.vx = (*iter).vx();
          lLep.vertexPCA.vy = (*iter).vy();
          lLep.vertexPCA.vz = (*iter).vz();
          lLep.charge = (*iter).charge();
    lLep.dBeamspot = (*iter).dB(pat::Muon::BS2D);
    lLep.dVertex = (*iter).dB(pat::Muon::PV2D);
    if ((*iter).triggerObjectMatches().size() > 0){
      lLep.hltMatchPaths = (*iter).triggerObjectMatches().at(0).pathNames(1,0);
      lLep.hltMatchFilters = (*iter).triggerObjectMatches().at(0).filterLabels();
    }
    //Construct a GenParticle if processing MC and a match has been found by PAT
          if (!processData_ && (*iter).genLepton()){
      HbbAnalysis::FourMomentum lGenFourVec;
            HbbAnalysis::GenVars lGen;
            lGenFourVec.px = (*iter).genLepton()->px();
            lGenFourVec.py = (*iter).genLepton()->py();
            lGenFourVec.pz = (*iter).genLepton()->pz();
            lGenFourVec.e = (*iter).genLepton()->energy();
            lGen.charge = (*iter).genLepton()->charge();
            lGen.pdgId = (*iter).genLepton()->pdgId();
            lGen.status = (*iter).genLepton()->status();
            lGen.vertex.vx = (*iter).genLepton()->vx();
            lGen.vertex.vy = (*iter).genLepton()->vy();
            lGen.vertex.vz = (*iter).genLepton()->vz();
      lLep.genMatches.push_back(HbbAnalysis::GenParticle(lGenFourVec, lGen));
          }


    //TESTING CODE
    /* 
    std::cout << "-----------------------------------------------" << std::endl;
    std::cout << "Muon: " << iEle << std::endl;
    std::cout << "--Found " << (*iter).triggerObjectMatches().size() << " matched trigger objects" << std::endl;
    
    std::cout << "--Path Names with Option (1,1): " << std::endl;
    for (unsigned i = 0; i < (*iter).triggerObjectMatches().size(); ++i){
      std::vector<std::string> hltpathNames = (*iter).triggerObjectMatches().at(i).pathNames(1,1);
      std::cout << "---hlt path names:" << std::endl;
      for (unsigned j = 0; j < hltpathNames.size(); ++j){
        std::cout << "----" << hltpathNames.at(j) << std::endl;
      }
    }
    std::cout << "--Path Names with Option (1,0): " << std::endl;
    for (unsigned i = 0; i < (*iter).triggerObjectMatches().size(); ++i){
      std::vector<std::string> hltpathNames = (*iter).triggerObjectMatches().at(i).pathNames(1,0);
      std::cout << "---hlt path names:" << std::endl;
      for (unsigned j = 0; j < hltpathNames.size(); ++j){
        std::cout << "----" << hltpathNames.at(j) << std::endl;
      }
    }
    std::cout << "--Path Names with Option (0,1): " << std::endl;
    for (unsigned i = 0; i < (*iter).triggerObjectMatches().size(); ++i){
      std::vector<std::string> hltpathNames = (*iter).triggerObjectMatches().at(i).pathNames(0,1);
      std::cout << "---hlt path names:" << std::endl;
      for (unsigned j = 0; j < hltpathNames.size(); ++j){
        std::cout << "----" << hltpathNames.at(j) << std::endl;
      }
    }
    std::cout << "--Path Names with Option (0,0): " << std::endl;
    for (unsigned i = 0; i < (*iter).triggerObjectMatches().size(); ++i){
      std::vector<std::string> hltpathNames = (*iter).triggerObjectMatches().at(i).pathNames(0,0);
      std::cout << "---hlt path names:" << std::endl;
      for (unsigned j = 0; j < hltpathNames.size(); ++j){
        std::cout << "----" << hltpathNames.at(j) << std::endl;
      }
    }

    for (unsigned i = 0; i < (*iter).triggerObjectMatches().size(); ++i){
      std::vector<std::string> HLTfilterLabels = (*iter).triggerObjectMatches().at(i).filterLabels();
      std::cout << "---HLT FilterLabels:" << std::endl;
      for (unsigned j = 0; j < HLTfilterLabels.size(); ++j){
        std::cout << "----" << HLTfilterLabels.at(j) << std::endl;
      }
    }
    */
    
  
    HbbAnalysis::MuTrkVars lTrk;
          reco::TrackRef glmuon = recoMuon->globalTrack();
          reco::TrackRef tkmuon = recoMuon->innerTrack();

          if ( glmuon.isAvailable() && glmuon.isNonnull() ) {
            if (aBeamSpot.isValid()) {
              lTrk.dxy = glmuon->dxy(aBeamSpot->position());
              lTrk.dz = glmuon->dz(aBeamSpot->position());
            }
            else {
              lTrk.dxy = 0;
              lTrk.dz = 0;
            }
            lTrk.normalizedChi2 = glmuon->normalizedChi2();
            lTrk.muonHits = glmuon->hitPattern().numberOfValidMuonHits();
            lTrk.charge = glmuon->charge();
          }
          else {
            lTrk.dxy = 0;
            lTrk.dz = 0;
            lTrk.normalizedChi2 = 0;
            lTrk.muonHits = 0;
            lTrk.charge = 0;
          }
          if ( tkmuon.isAvailable() && tkmuon.isNonnull() ) {
            lTrk.trackerHits = tkmuon->hitPattern().numberOfValidTrackerHits();
            lTrk.pixelHits = tkmuon->hitPattern().numberOfValidPixelHits();
          }
          else {
            lTrk.trackerHits = 0;
            lTrk.pixelHits = 0;
          }

          HbbAnalysis::MuIsoVars lIsoR03;
          lIsoR03.sumPt = recoMuon->isolationR03().sumPt;
          lIsoR03.emEt = recoMuon->isolationR03().emEt;
          lIsoR03.hadEt = recoMuon->isolationR03().hadEt;
          lIsoR03.nTracks = recoMuon->isolationR03().nTracks;

          HbbAnalysis::MuIsoVars lIsoR05;
          lIsoR05.sumPt = recoMuon->isolationR05().sumPt;
          lIsoR05.emEt = recoMuon->isolationR05().emEt;
          lIsoR05.hadEt = recoMuon->isolationR05().hadEt;
          lIsoR05.nTracks = recoMuon->isolationR05().nTracks;

          HbbAnalysis::MuIDVars lID;
          //if (recoMuon->isGlobalMuon()) lID.type = 1;
          //else if (recoMuon->isTrackerMuon()) lID.type = 2;
          //else if (recoMuon->isStandAloneMuon()) lID.type = 3;
          //else if (recoMuon->isCaloMuon()) lID.type = 4;
          //else lID.type = 0;
          lID.type =  (recoMuon->isCaloMuon()<<3) | (recoMuon->isStandAloneMuon()<<2) | (recoMuon->isTrackerMuon()<<1) | recoMuon->isGlobalMuon();

          if (muon::isGoodMuon(*recoMuon,muon::AllGlobalMuons)) lID.ids.push_back(1);
          if (muon::isGoodMuon(*recoMuon,muon::AllStandAloneMuons)) lID.ids.push_back(2);
          if (muon::isGoodMuon(*recoMuon,muon::AllTrackerMuons)) lID.ids.push_back(3);
          if (muon::isGoodMuon(*recoMuon,muon::TrackerMuonArbitrated)) lID.ids.push_back(4);
          if (muon::isGoodMuon(*recoMuon,muon::AllArbitrated)) lID.ids.push_back(5);
          if (muon::isGoodMuon(*recoMuon,muon::GlobalMuonPromptTight)) lID.ids.push_back(6);
          if (muon::isGoodMuon(*recoMuon,muon::TMLastStationLoose)) lID.ids.push_back(7);
          if (muon::isGoodMuon(*recoMuon,muon::TMLastStationTight)) lID.ids.push_back(8);
          if (muon::isGoodMuon(*recoMuon,muon::TM2DCompatibilityLoose)) lID.ids.push_back(9);
          if (muon::isGoodMuon(*recoMuon,muon::TM2DCompatibilityTight)) lID.ids.push_back(10);
          if (muon::isGoodMuon(*recoMuon,muon::TMOneStationLoose)) lID.ids.push_back(11);
          if (muon::isGoodMuon(*recoMuon,muon::TMOneStationTight)) lID.ids.push_back(12);
          if (muon::isGoodMuon(*recoMuon,muon::TMLastStationOptimizedLowPtLoose)) lID.ids.push_back(13);
          if (muon::isGoodMuon(*recoMuon,muon::TMLastStationOptimizedLowPtTight)) lID.ids.push_back(14);
          if (muon::isGoodMuon(*recoMuon,muon::GMTkChiCompatibility)) lID.ids.push_back(15);
          if (muon::isGoodMuon(*recoMuon,muon::GMStaChiCompatibility)) lID.ids.push_back(16);
          if (muon::isGoodMuon(*recoMuon,muon::GMTkKinkTight)) lID.ids.push_back(17);
          if (muon::isGoodMuon(*recoMuon,muon::TMLastStationAngLoose)) lID.ids.push_back(18);
          if (muon::isGoodMuon(*recoMuon,muon::TMLastStationAngTight)) lID.ids.push_back(19);
          if (muon::isGoodMuon(*recoMuon,muon::TMOneStationAngLoose)) lID.ids.push_back(20);
          if (muon::isGoodMuon(*recoMuon,muon::TMOneStationAngTight)) lID.ids.push_back(21);
          if (muon::isGoodMuon(*recoMuon,muon::TMLastStationOptimizedBarrelLowPtLoose)) lID.ids.push_back(22);
          if (muon::isGoodMuon(*recoMuon,muon::TMLastStationOptimizedBarrelLowPtTight)) lID.ids.push_back(23);

          lID.caloCompat = recoMuon->caloCompatibility();
          lID.segCompat = muon::segmentCompatibility(*recoMuon);
          lID.nChambers = recoMuon->numberOfChambers();
          lID.nMatchesLoose = recoMuon->numberOfMatches(reco::Muon::NoArbitration);
          lID.nMatchesMedium = recoMuon->numberOfMatches(reco::Muon::SegmentArbitration);
          lID.nMatchesTight = recoMuon->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);

          aVec.push_back(HbbAnalysis::Muon(lFourVec,lLep,lTrk,lIsoR03,lIsoR05,lID));
        }
    }
  }

}//HbbMuons

/*
void HbbTreeMaker::HbbTaus(const edm::Handle<std::vector<pat::Tau> > & aCol, 
                           const edm::Handle<std::vector<reco::Vertex> > & aRecoVertices,
                           std::vector<HbbAnalysis::Tau> & aVec)
{//HbbTaus
  
  if (aCol.isValid()){
    if (aCol->size() > 0) {
        
      unsigned int iEle = 0;
      for (std::vector<pat::Tau>::const_iterator iter = aCol->begin(); 
           iter != aCol->end(); 
           ++iter,++iEle) 
        {
          //if (debug_ > 1) std::cout << "**** Ele #" << iEle << ", pT,eta,phi =  " << (*iter).pt() << " " << (*iter).eta() << " " << (*iter).phi() << std::endl;

          HbbAnalysis::GenVars lGen;
          if (!processData_ && (*iter).genLepton()){
            lGen.valid = true;
            lGen.E = (*iter).genLepton()->energy();
            lGen.pT = (*iter).genLepton()->pt();
            lGen.eta = (*iter).genLepton()->eta();
            lGen.y = (*iter).genLepton()->y();
            lGen.phi = (*iter).genLepton()->phi();
            lGen.charge = (*iter).genLepton()->charge();
            lGen.pdgId = (*iter).genLepton()->pdgId();
            lGen.status = (*iter).genLepton()->status();
            lGen.mass = (*iter).genLepton()->mass();
            lGen.vx = (*iter).genLepton()->vx();
            lGen.vy = (*iter).genLepton()->vy();
            lGen.vz = (*iter).genLepton()->vz();
          }
          else {
            lGen.valid = false;
            lGen.E = 0;
            lGen.pT = 0;
            lGen.eta = 0;
            lGen.y = 0;
            lGen.phi = 0;
            lGen.charge = 0;
            lGen.pdgId = 0;
            lGen.status = 0;
            lGen.mass = 0;
            lGen.vx = 0;
            lGen.vy = 0;
            lGen.vz = 0;
          }


          HbbAnalysis::GenVars lGenJet;
          if (!processData_ && (*iter).genJet()){
            lGenJet.valid = true;
            lGenJet.E = (*iter).genJet()->energy();
            lGenJet.pT = (*iter).genJet()->pt();
            lGenJet.eta = (*iter).genJet()->eta();
            lGenJet.y = (*iter).genJet()->y();
            lGenJet.phi = (*iter).genJet()->phi();
            lGenJet.charge = (*iter).genJet()->charge();
            lGenJet.pdgId = (*iter).genJet()->pdgId();
            lGenJet.status = (*iter).genJet()->status();
            lGenJet.mass = (*iter).genJet()->mass();
            lGenJet.vx = (*iter).genJet()->vx();
            lGenJet.vy = (*iter).genJet()->vy();
            lGenJet.vz = (*iter).genJet()->vz();
          }
          else {
            lGenJet.valid = false;
            lGenJet.E = 0;
            lGenJet.pT = 0;
            lGenJet.eta = 0;
            lGenJet.y = 0;
            lGenJet.phi = 0;
            lGenJet.charge = 0;
            lGenJet.pdgId = 0;
            lGenJet.status = 0;
            lGenJet.mass = 0;
            lGenJet.vx = 0;
            lGenJet.vy = 0;
            lGenJet.vz = 0;
          }

          HbbAnalysis::BaseVars lReco;
          lReco.E = (*iter).energy();
          lReco.pT = (*iter).pt();
          lReco.eta = (*iter).eta();
          lReco.y = (*iter).y();
          lReco.phi = (*iter).phi();
          lReco.charge = (*iter).charge();
          lReco.vx = (*iter).vx();
          lReco.vy = (*iter).vy();
          lReco.vz = (*iter).vz();

          bool isCalo = (*iter).isCaloTau();
          bool isPF = (*iter).isPFTau();

          assert (isCalo == !isPF);

          HbbAnalysis::TrkVars lLead;
          if (isCalo) {
            lLead.signedSipt = (*iter).leadTracksignedSipt();

            reco::TrackRef lCaloTrk = (*iter).leadTrack();

            if ( lCaloTrk.isAvailable() && lCaloTrk.isNonnull() ) {
              lLead.pT = lCaloTrk->pt();
              lLead.eta = lCaloTrk->eta();
              lLead.phi = lCaloTrk->phi();
      
              lLead.matchDist = reco::deltaR(lCaloTrk->momentum(), (*iter).p4());

              if ( aRecoVertices->size() >= 1 ) {
                const reco::Vertex& thePrimaryEventVertex = (*aRecoVertices->begin());
                lLead.IPxy = lCaloTrk->dxy(thePrimaryEventVertex.position());
                lLead.IPz = lCaloTrk->dz(thePrimaryEventVertex.position());
              }
              else {
                lLead.IPxy = 0;
                lLead.IPz = 0;
              }
            }
            else {
              lLead.pT = 0;
              lLead.eta = 0;
              lLead.phi = 0;
              lLead.matchDist = 0;
              lLead.IPxy = 0;
              lLead.IPz = 0;
              lLead.signedSipt = 0;
            }
          }//calo
          else {
            lLead.signedSipt = (*iter).leadPFChargedHadrCandsignedSipt();

            reco::PFCandidateRef lHadrCand = (*iter).leadPFChargedHadrCand();

            if ( lHadrCand.isAvailable() && lHadrCand.isNonnull() ) {
              lLead.pT = lHadrCand->pt();
              lLead.eta = lHadrCand->eta();
              lLead.phi = lHadrCand->phi();

              lLead.matchDist = reco::deltaR(lHadrCand->p4(), (*iter).p4());

              reco::TrackRef lTrk = lHadrCand->trackRef();

              if ( aRecoVertices->size() >= 1) {
                const reco::Vertex& thePrimaryEventVertex = (*aRecoVertices->begin());
                lLead.IPxy = lTrk->dxy(thePrimaryEventVertex.position());
                lLead.IPz = lTrk->dz(thePrimaryEventVertex.position());
              }
              else {
                lLead.IPxy = 0;
                lLead.IPz = 0;
              }
            }
            else {
              lLead.pT = 0;
              lLead.eta = 0;
              lLead.phi = 0;
              lLead.matchDist = 0;
              lLead.IPxy = 0;
              lLead.IPz = 0;
              lLead.signedSipt = 0;
            }

          }//pf


          const std::vector<std::pair<std::string,float> > &  lIDs = (*iter).tauIDs();
          bool lPrint = false;
          //a few warning if additional discriminants are found:
          if ((isPF && lIDs.size() != 16) || 
              (isCalo && lIDs.size() != 3))
            lPrint = true;

          if (lPrint) {
            std::cout << "!!!!!!! Discriminants changed, please update tree !!!!!!!" << std::endl;
            std::cout << "--- isCaloTau = " << isCalo << ", isPFTau = " << isPF << std::endl;
            std::cout << "------ ID names = " << std::endl;
          }


          HbbAnalysis::CaloTauIDVars lcaloId;
          lcaloId.byIsolation = 0;
          lcaloId.byLeadingTrackFinding = 0;
          lcaloId.byLeadingTrackPtCut = 0;

          HbbAnalysis::PFTauIDVars lpfId;
          lpfId.byLeadingTrackFinding = 0;
          lpfId.byLeadingTrackPtCut = 0;
          lpfId.byTrackIsolation = 0;
          lpfId.byECALIsolation = 0;
          lpfId.byIsolation = 0;
          lpfId.againstElectron = 0;
          lpfId.againstMuon = 0;
          lpfId.byIsolationUsingLeadingPion = 0;
          lpfId.byTaNC = 0;
          lpfId.byTaNCfrHalfPercent = 0;
          lpfId.byTaNCfrOnePercent = 0;
          lpfId.byTaNCfrQuarterPercent = 0;
          lpfId.byTaNCfrTenthPercent = 0;
          lpfId.ecalIsolationUsingLeadingPion = 0;
          lpfId.leadingPionPtCut = 0;
          lpfId.trackIsolationUsingLeadingPion = 0;


          for (unsigned int id(0); id<lIDs.size(); ++id){

            std::string lName = lIDs.at(id).first; 
            float lDiscri = lIDs.at(id).second;

            if (lPrint) std::cout << "--------- " << lName << " = " << lDiscri << std::endl;
            //pf

            if (isPF) {
              if (lName.find("leadingTrackFinding") != lName.npos) lpfId.byLeadingTrackFinding = lDiscri;
              if (lName.find("leadingTrackPtCut") != lName.npos) lpfId.byLeadingTrackPtCut = lDiscri;
              if (lName.find("trackIsolationUsingLeadingPion") != lName.npos) lpfId.trackIsolationUsingLeadingPion = lDiscri;
              if (lName.find("trackIsolation") != lName.npos && 
                  lName.find("trackIsolationUsingLeadingPion") == lName.npos) lpfId.byTrackIsolation = lDiscri;
              if (lName.find("ecalIsolationUsingLeadingPion") != lName.npos) lpfId.ecalIsolationUsingLeadingPion = lDiscri;
              if (lName.find("ecalIsolation") != lName.npos &&
                  lName.find("ecalIsolationUsingLeadingPion") == lName.npos) lpfId.byECALIsolation = lDiscri;
              if (lName.find("byIsolationUsingLeadingPion") != lName.npos) lpfId.byIsolationUsingLeadingPion = lDiscri;
              if (lName.find("byIsolation") != lName.npos &&
                  lName.find("byIsolationUsingLeadingPion") == lName.npos) lpfId.byIsolation = lDiscri;
              if (lName.find("againstElectron") != lName.npos) lpfId.againstElectron = lDiscri;
              if (lName.find("againstMuon") != lName.npos) lpfId.againstMuon = lDiscri;
              if (lName.find("byTaNCfrHalfPercent") != lName.npos) lpfId.byTaNCfrHalfPercent = lDiscri;
              if (lName.find("byTaNCfrOnePercent") != lName.npos) lpfId.byTaNCfrOnePercent = lDiscri;
              if (lName.find("byTaNCfrQuarterPercent") != lName.npos) lpfId.byTaNCfrQuarterPercent = lDiscri;
              if (lName.find("byTaNCfrTenthPercent") != lName.npos) lpfId.byTaNCfrTenthPercent = lDiscri;
              if (lName.find("byTaNC") != lName.npos && 
                  lName.find("byTaNCfr") == lName.npos) lpfId.byTaNC = lDiscri;
              if (lName.find("leadingPionPtCut") != lName.npos) lpfId.leadingPionPtCut = lDiscri;
            }
            if (isCalo){
              if (lName.find("byIsolation") != lName.npos) lcaloId.byIsolation = lDiscri;
              if (lName.find("leadingTrackFinding") != lName.npos) lcaloId.byLeadingTrackFinding = lDiscri;
              if (lName.find("leadingTrackPtCut") != lName.npos) lcaloId.byLeadingTrackPtCut = lDiscri;
            }


          }


          if (isCalo){

            HbbAnalysis::CaloTauIsoVars lcaloIso;
            lcaloIso.nIsoTracks = (*iter).isolationTracks().size();
            lcaloIso.nSigTracks = (*iter).signalTracks().size();
            lcaloIso.leadTrackHCAL3x3hitsEtSum = (*iter).leadTrackHCAL3x3hitsEtSum();
            lcaloIso.leadTrackHCAL3x3hottesthitDEta = (*iter).leadTrackHCAL3x3hottesthitDEta();
            lcaloIso.signalTracksInvariantMass = (*iter).signalTracksInvariantMass();
            lcaloIso.tracksInvariantMass = (*iter).TracksInvariantMass(); 
            lcaloIso.isolationTracksPtSum = (*iter).isolationTracksPtSum();
            lcaloIso.isolationECALhitsEtSum = (*iter).isolationECALhitsEtSum();
            lcaloIso.maximumHCALhitEt = (*iter).maximumHCALhitEt();

            HbbAnalysis::Tau lObj(lGen,lGenJet,lReco,lLead,lcaloIso,lcaloId);
            aVec.push_back(lObj);
          }
          else {

            HbbAnalysis::PFTauIsoVars lpfIso;
            lpfIso.nSigCands = (*iter).signalPFCands().size();
            lpfIso.nIsoCands = (*iter).isolationPFCands().size();
            lpfIso.maximumHCALPFClusterEt = (*iter).maximumHCALPFClusterEt();
            //lpfIso.emFraction = (*iter).emFraction();
            lpfIso.hcalTotOverPLead = (*iter).hcalTotOverPLead();
            lpfIso.hcalMaxOverPLead = (*iter).hcalMaxOverPLead();
            lpfIso.hcal3x3OverPLead = (*iter).hcal3x3OverPLead();
            lpfIso.ecalStripSumEOverPLead = (*iter).ecalStripSumEOverPLead();
            lpfIso.bremsRecoveryEOverPLead = (*iter).bremsRecoveryEOverPLead();

            HbbAnalysis::PFTauCandVars lHadr;
            lHadr.nSigCands = (*iter).signalPFChargedHadrCands().size();
            lHadr.nIsoCands = (*iter).isolationPFChargedHadrCands().size();
            lHadr.isolationPtSum = (*iter).isolationPFChargedHadrCandsPtSum();

            HbbAnalysis::PFTauCandVars lNeutr;
            lNeutr.nSigCands = (*iter).signalPFNeutrHadrCands().size();
            lNeutr.nIsoCands = (*iter).isolationPFNeutrHadrCands().size();
            lNeutr.isolationPtSum = (*iter).neutralHadronIso();

            HbbAnalysis::PFTauCandVars lGamma;
            lGamma.nSigCands = (*iter).signalPFGammaCands().size();
            lGamma.nIsoCands = (*iter).isolationPFGammaCands().size();
            lGamma.isolationPtSum = (*iter).isolationPFGammaCandsEtSum();

            HbbAnalysis::PFTauEleIDVars lEleId;
            if ( (*iter).electronPreIDTrack().isAvailable() && (*iter).electronPreIDTrack().isNonnull() ) {
              lEleId.pT = (*iter).electronPreIDTrack()->pt();
              lEleId.eta = (*iter).electronPreIDTrack()->eta();
              lEleId.phi = (*iter).electronPreIDTrack()->phi();
            }
            else {
              lEleId.pT = 0;
              lEleId.eta = 0;
              lEleId.phi = 0;
            }

            lEleId.output = (*iter).electronPreIDOutput();
            lEleId.decision = (*iter).electronPreIDDecision();

            HbbAnalysis::PFTauMuIDVars lMuId;
            lMuId.caloCompat = (*iter).caloComp();
            lMuId.segCompat = (*iter).segComp();
            lMuId.decision = (*iter).muonDecision();

            HbbAnalysis::Tau lObj(lGen,lGenJet,lReco,lLead,lpfIso,lHadr,lNeutr,lGamma,lpfId,lEleId,lMuId);
            aVec.push_back(lObj);
          }
        }
    }
  }

}//HbbTaus
*/

void HbbTreeMaker::HbbJets(const edm::Handle<std::vector<pat::Jet> > & aCol, 
                           const HbbAnalysis::JetFlavour & aJetFlav, 
                           const edm::Handle<reco::GenParticleCollection> & aGenParticles,
                           std::vector<HbbAnalysis::Jet> & aVec)
{//HbbJets
  
  if (aCol.isValid()){//valid
    if (aCol->size() > 0) {//non empty
        
      unsigned int iEle = 0;
      for (std::vector<pat::Jet>::const_iterator iter = aCol->begin(); 
           iter != aCol->end(); 
           ++iter,++iEle) 
        {//loop on element
          //if (debug_ > 1) std::cout << "**** Ele #" << iEle << ", pT,eta,phi =  " << (*iter).pt() << " " << (*iter).eta() << " " << (*iter).phi() << std::endl;

    //----------------------------------------------------------------------
    HbbAnalysis::FourMomentum lFourVec;
    //----------------------------------------------------------------------
          lFourVec.e = (*iter).energy();
          lFourVec.px = (*iter).px();
          lFourVec.py = (*iter).py();
          lFourVec.pz = (*iter).pz();
    //----------------------------------------------------------------------
    HbbAnalysis::JetVars lJet;
    //----------------------------------------------------------------------
          unsigned int lFlav = 0;
          if (!processData_ && aJetFlav.nPartons()) {
            lFlav = aJetFlav.partonMatchingGenJet((*iter),aGenParticles,0.4).second;
          }
          lJet.mcFlavour = lFlav;
          lJet.partonFlavour = (*iter).partonFlavour();
          lJet.nAssociatedTracks = ((*iter).associatedTracks()).size();
          lJet.etaMean = (*iter).etaPhiStatistics().etaMean;
          lJet.phiMean = (*iter).etaPhiStatistics().phiMean;
          lJet.etaEtaMoment = (*iter).etaPhiStatistics().etaEtaMoment;
          lJet.phiPhiMoment = (*iter).etaPhiStatistics().phiPhiMoment;
          lJet.etaPhiMoment = (*iter).etaPhiStatistics().etaPhiMoment;
    //Construct a GenParticle if processing MC and a match has been found by PAT
          if (!processData_ && (*iter).genParton()){
      HbbAnalysis::FourMomentum lGenFourVec;
            HbbAnalysis::GenVars lGen;
            lGenFourVec.px = (*iter).genParton()->px();
            lGenFourVec.py = (*iter).genParton()->py();
            lGenFourVec.pz = (*iter).genParton()->pz();
            lGenFourVec.e = (*iter).genParton()->energy();
            lGen.charge = (*iter).genParton()->charge();
            lGen.pdgId = (*iter).genParton()->pdgId();
            lGen.status = (*iter).genParton()->status();
            lGen.vertex.vx = (*iter).genParton()->vx();
            lGen.vertex.vy = (*iter).genParton()->vy();
            lGen.vertex.vz = (*iter).genParton()->vz();
      lJet.genMatches.push_back(HbbAnalysis::GenParticle(lGenFourVec, lGen));
          }
    //Construct a GenJet if processing MC and a match has been found by PAT
          if (!processData_ && (*iter).genJet()){
      HbbAnalysis::FourMomentum lGenJetFourVec;
      HbbAnalysis::GenJetVars lGenJet;
            lGenJetFourVec.px = (*iter).genJet()->px();
            lGenJetFourVec.py = (*iter).genJet()->py();
            lGenJetFourVec.pz = (*iter).genJet()->pz();
            lGenJetFourVec.e = (*iter).genJet()->energy();
            lGenJet.flavour = 0;  //Flavour unknown in this context (will be filled for parton jets)
            lGenJet.mass = (*iter).genJet()->mass();
      lGenJet.nConstituents = (*iter).genJet()->getGenConstituents().size();
      lJet.genJetMatches.push_back(HbbAnalysis::GenJet(lGenJetFourVec, lGenJet));
          }
               
          HbbAnalysis::JetECorVars lEnergyCorrections;
          if ((*iter).jecSetsAvailable()){
            std::vector<std::string> lJECnames = (*iter).availableJECSets();
      for (unsigned int iSet(0); iSet < lJECnames.size(); ++iSet){
        if (debug_) std::cout << " -- JEC SET " << iSet << " " << lJECnames[iSet] << std::endl;
        std::vector<std::string> lJEClevels = (*iter).availableJECLevels(iSet);
        for (unsigned int iLev(0); iLev<lJEClevels.size(); ++iLev){
          if (debug_) std::cout << " ---- JEC level " << iLev << " " << lJEClevels[iLev] 
            << " , factor = " << (*iter).jecFactor(iLev,pat::JetCorrFactors::NONE,iSet)
              << std::endl;
          lEnergyCorrections.Levels[(lJEClevels[iLev])] = (*iter).jecFactor(iLev,pat::JetCorrFactors::NONE,iSet);
        }
      }
    }
    //Care must be taken here with getUncertainty calls, this appears to "reset" the corrector object and if called again
    //without re-setting Eta and Pt will cause a crash.
    if(jecUncCalo && jecUncPF && jecUncJPT){
      if ((*iter).isCaloJet()) {
        jecUncCalo->setJetEta((*iter).eta());
        jecUncCalo->setJetPt((*iter).pt());
        lEnergyCorrections.Levels["Uncertainty"] = jecUncCalo->getUncertainty(true);
      } else if ((*iter).isPFJet()) {
        jecUncPF->setJetEta((*iter).eta());
        jecUncPF->setJetPt((*iter).pt());
        lEnergyCorrections.Levels["Uncertainty"] = jecUncPF->getUncertainty(true);
      } else if ((*iter).isJPTJet()) {
        jecUncJPT->setJetEta((*iter).eta());
        jecUncJPT->setJetPt((*iter).pt());
        lEnergyCorrections.Levels["Uncertainty"] = jecUncJPT->getUncertainty(true);
      }
    }
          
          HbbAnalysis::JetBtagVars lBtag;
          lBtag.cSV = (*iter).bDiscriminator("combinedSecondaryVertexBJetTags");
          lBtag.cSVMVA = (*iter).bDiscriminator("combinedSecondaryVertexMVABJetTags");
          lBtag.iPMVA = (*iter).bDiscriminator("impactParameterMVABJetTags");
          lBtag.bProba = (*iter).bDiscriminator("jetBProbabilityBJetTags");
          lBtag.probability = (*iter).bDiscriminator("jetProbabilityBJetTags");
          lBtag.sSVHP = (*iter).bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
          lBtag.sSVHE = (*iter).bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
          lBtag.softElectronByPt = ((*iter).bDiscriminator("softElectronByPtBJetTags") < -100) ? -100 : (*iter).bDiscriminator("softElectronByPtBJetTags");
          lBtag.softElectronByIP3d = ((*iter).bDiscriminator("softElectronByIP3dBJetTags") < -100) ? -100 : (*iter).bDiscriminator("softElectronByIP3dBJetTags");
          lBtag.softMuon = ((*iter).bDiscriminator("softMuonBJetTags") < -100) ? -100 : (*iter).bDiscriminator("softMuonBJetTags");
          lBtag.softMuonByPt = ((*iter).bDiscriminator("softMuonByPtBJetTags") < -100) ? -100 : (*iter).bDiscriminator("softMuonByPtBJetTags");
          lBtag.softMuonByIP3d = ((*iter).bDiscriminator("softMuonByIP3dBJetTags") < -100) ? -100 : (*iter).bDiscriminator("softMuonByIP3dBJetTags");
          lBtag.tCHE = (*iter).bDiscriminator("trackCountingHighEffBJetTags");
          lBtag.tCHP = (*iter).bDiscriminator("trackCountingHighPurBJetTags");

          //std::cout << " -- New values of b-discri for jet : " << iEle << std::endl
          //        << " ---- softElecs: " << lBtag.softElectronByPt << " " << lBtag.softElectronByIP3d << std::endl
          //        << " ---- softMus: " << lBtag.softMuon << " " << lBtag.softMuonByPt << " " << lBtag.softMuonByIP3d << std::endl;


    std::vector<HbbAnalysis::JetSecVertexVars> jetSecVertexVars;
    for (unsigned iVtx = 0; iVtx < (*iter).tagInfoSecondaryVertex()->nVertices(); ++iVtx){
      HbbAnalysis::JetSecVertexVars jetSecVertexVarsTmp;
      math::XYZTLorentzVectorD vec = (*iter).tagInfoSecondaryVertex()->secondaryVertex(iVtx).p4();
      jetSecVertexVarsTmp.fourVec.e  = vec.E();
      jetSecVertexVarsTmp.fourVec.px = vec.Px();
      jetSecVertexVarsTmp.fourVec.py = vec.Py();
      jetSecVertexVarsTmp.fourVec.pz = vec.Pz();
      jetSecVertexVarsTmp.vertex.vx = (*iter).tagInfoSecondaryVertex()->secondaryVertex(iVtx).x();
      jetSecVertexVarsTmp.vertex.vy = (*iter).tagInfoSecondaryVertex()->secondaryVertex(iVtx).y();
      jetSecVertexVarsTmp.vertex.vz = (*iter).tagInfoSecondaryVertex()->secondaryVertex(iVtx).z();
      jetSecVertexVars.push_back(jetSecVertexVarsTmp);
    }


          bool isCalo = (*iter).isCaloJet();
          bool isPF = (*iter).isPFJet();
          bool isJPT = (*iter).isJPTJet();

          //std::cout << " -- isCalo = " << isCalo << ", isPF = " << isPF << std::endl;

          //assert (isCalo == !isPF);
          HbbAnalysis::JetIDVars lId;
          if (isCalo || isJPT) {
            lId.fHPD = (*iter).jetID().fHPD; 
            lId.fRBX = (*iter).jetID().fRBX;
            lId.n90Hits = (*iter).jetID().n90Hits;
            lId.hitsInN90 = (*iter).jetID().hitsInN90;
            lId.restrictedEMF = (*iter).jetID().restrictedEMF;
            lId.nHCALTowers = (*iter).jetID().nHCALTowers;
            lId.nECALTowers = (*iter).jetID().nECALTowers;
          }
          
    if (isCalo){
            HbbAnalysis::CaloJetVars lCalo;
            lCalo.maxEInEmTowers = (*iter).maxEInEmTowers();
            lCalo.maxEInHadTowers = (*iter).maxEInHadTowers();
            lCalo.energyFractionHadronic = (*iter).energyFractionHadronic();
            lCalo.emEnergyFraction = (*iter).emEnergyFraction();
            lCalo.hadEnergyInHB = (*iter).hadEnergyInHB();
            lCalo.hadEnergyInHO = (*iter).hadEnergyInHO();
            lCalo.hadEnergyInHE = (*iter).hadEnergyInHE();
            lCalo.hadEnergyInHF = (*iter).hadEnergyInHF();
            lCalo.emEnergyInEB = (*iter).emEnergyInEB();
            lCalo.emEnergyInEE = (*iter).emEnergyInEE();
            lCalo.emEnergyInHF = (*iter).emEnergyInHF();
            lCalo.towersArea = (*iter).towersArea();
            lCalo.n90 = (*iter).n90();
            lCalo.n60 = (*iter).n60();
            HbbAnalysis::Jet lObj(lFourVec,lJet,lCalo,lBtag,lId,lEnergyCorrections);
      lObj.secVertexVars(jetSecVertexVars);
            aVec.push_back(lObj);
          }
          
    HbbAnalysis::JPTPFJetVars lJPTPF;
          if (isPF || isJPT){
            lJPTPF.chargedEmEnergy = (*iter).chargedEmEnergy();
            lJPTPF.neutralEmEnergy = (*iter).neutralEmEnergy();
            lJPTPF.muonMultiplicity = (*iter).muonMultiplicity();
            lJPTPF.chargedHadronEnergy = (*iter).chargedHadronEnergy();
            lJPTPF.neutralHadronEnergy = (*iter).neutralHadronEnergy();
            lJPTPF.chargedMultiplicity = (*iter).chargedMultiplicity();
            lJPTPF.chargedHadronEnergyFraction = (*iter).chargedHadronEnergyFraction();
            lJPTPF.neutralHadronEnergyFraction = (*iter).neutralHadronEnergyFraction();
            lJPTPF.chargedEmEnergyFraction = (*iter).chargedEmEnergyFraction();
            lJPTPF.neutralEmEnergyFraction = (*iter).neutralEmEnergyFraction();
          }
          if (isJPT) {
            HbbAnalysis::JPTJetVars lJPT;
      reco::TrackRefVector pionsInVInC = (*iter).pionsInVertexInCalo();
      reco::TrackRefVector pionsInVOuC = (*iter).pionsInVertexOutCalo();
      reco::TrackRefVector pionsOuVInC = (*iter).pionsOutVertexInCalo();
      
      reco::TrackRefVector muonsInVInC = (*iter).muonsInVertexInCalo();
      reco::TrackRefVector muonsInVOuC = (*iter).muonsInVertexOutCalo();
      reco::TrackRefVector muonsOuVInC = (*iter).muonsOutVertexInCalo();

      reco::TrackRefVector elecsInVInC = (*iter).elecsInVertexInCalo();
      reco::TrackRefVector elecsInVOuC = (*iter).elecsInVertexOutCalo();
      reco::TrackRefVector elecsOuVInC = (*iter).elecsOutVertexInCalo();

      /*
      unsigned test_track_count = pionsInVInC.size() + pionsInVOuC.size() + pionsOuVInC.size();
      std::cout << "Total Pion Track Count = " << test_track_count << std::endl;
      */

      for (reco::TrackRefVector::const_iterator trk_iter = pionsInVInC.begin(); 
          trk_iter != pionsInVInC.end(); 
          ++trk_iter) {
        lJPT.particleStatusPt.push_back(std::pair<unsigned int,double>(1,(*trk_iter)->pt()));
      }
      for (reco::TrackRefVector::const_iterator trk_iter = pionsInVOuC.begin(); 
          trk_iter != pionsInVOuC.end(); 
          ++trk_iter) {
        lJPT.particleStatusPt.push_back(std::pair<unsigned int,double>(2,(*trk_iter)->pt()));
      }
      for (reco::TrackRefVector::const_iterator trk_iter = pionsOuVInC.begin(); 
          trk_iter != pionsOuVInC.end(); 
          ++trk_iter) {
        lJPT.particleStatusPt.push_back(std::pair<unsigned int,double>(3,(*trk_iter)->pt()));
      }

      for (reco::TrackRefVector::const_iterator trk_iter = muonsInVInC.begin(); 
          trk_iter != muonsInVInC.end(); 
          ++trk_iter) {
        lJPT.particleStatusPt.push_back(std::pair<unsigned int,double>(4,(*trk_iter)->pt()));
      }
      for (reco::TrackRefVector::const_iterator trk_iter = muonsInVOuC.begin(); 
          trk_iter != muonsInVOuC.end(); 
          ++trk_iter) {
        lJPT.particleStatusPt.push_back(std::pair<unsigned int,double>(5,(*trk_iter)->pt()));
      }
      for (reco::TrackRefVector::const_iterator trk_iter = muonsOuVInC.begin(); 
          trk_iter != muonsOuVInC.end(); 
          ++trk_iter) {
        lJPT.particleStatusPt.push_back(std::pair<unsigned int,double>(6,(*trk_iter)->pt()));
      }
      
      for (reco::TrackRefVector::const_iterator trk_iter = elecsInVInC.begin(); 
          trk_iter != elecsInVInC.end(); 
          ++trk_iter) {
        lJPT.particleStatusPt.push_back(std::pair<unsigned int,double>(7,(*trk_iter)->pt()));
      }
      for (reco::TrackRefVector::const_iterator trk_iter = elecsInVOuC.begin(); 
          trk_iter != elecsInVOuC.end(); 
          ++trk_iter) {
        lJPT.particleStatusPt.push_back(std::pair<unsigned int,double>(8,(*trk_iter)->pt()));
      }
      for (reco::TrackRefVector::const_iterator trk_iter = elecsOuVInC.begin(); 
          trk_iter != elecsOuVInC.end(); 
          ++trk_iter) {
        lJPT.particleStatusPt.push_back(std::pair<unsigned int,double>(9,(*trk_iter)->pt()));
      }

      //            const reco::TrackRefVector& pionsInVertexInCalo () const{return jptSpecific().pionsInVertexInCalo; }
//          const reco::TrackRefVector& pionsInVertexOutCalo() const{return jptSpecific().pionsInVertexOutCalo;}
//          const reco::TrackRefVector& pionsOutVertexInCalo() const{return jptSpecific().pionsOutVertexInCalo;}
//          const reco::TrackRefVector& muonsInVertexInCalo () const{return jptSpecific().muonsInVertexInCalo; }
//          const reco::TrackRefVector& muonsInVertexOutCalo() const{return jptSpecific().muonsInVertexOutCalo;}
//          const reco::TrackRefVector& muonsOutVertexInCalo() const{return jptSpecific().muonsOutVertexInCalo;}
//          const reco::TrackRefVector& elecsInVertexInCalo () const{return jptSpecific().elecsInVertexInCalo; }
//          const reco::TrackRefVector& elecsInVertexOutCalo() const{return jptSpecific().elecsInVertexOutCalo;}
//          const reco::TrackRefVector& elecsOutVertexInCalo() const{return jptSpecific().elecsOutVertexInCalo;}
            lJPT.zspCorrection = (*iter).zspCorrection();
            lJPT.elecMultiplicity = (*iter).elecMultiplicity();

            HbbAnalysis::CaloJetVars lCalo;
            const reco::JPTJet * lJetRef = dynamic_cast<const reco::JPTJet *>((*iter).originalObject());
            edm::RefToBase<reco::Jet> jptjetRef = lJetRef->getCaloJetRef();

            // Dynamic_cast to CaloJet object
            reco::CaloJet const * rawcalojet = 
              dynamic_cast<reco::CaloJet const *>( &* jptjetRef);
      lJPT.calopT = rawcalojet->pt();
      lJPT.caloEta = rawcalojet->eta();
      if (debug_) {
        std::cout << " -- Raw CaloJet pT = " << rawcalojet->pt() << std::endl;     
        std::cout << " -- Raw CaloJet Eta = " << rawcalojet->eta() << std::endl;
      }
            lCalo.maxEInEmTowers = rawcalojet->maxEInEmTowers();
            lCalo.maxEInHadTowers = rawcalojet->maxEInHadTowers();
            lCalo.energyFractionHadronic = rawcalojet->energyFractionHadronic();
            lCalo.emEnergyFraction = rawcalojet->emEnergyFraction();
            lCalo.hadEnergyInHB = rawcalojet->hadEnergyInHB();
            lCalo.hadEnergyInHO = rawcalojet->hadEnergyInHO();
            lCalo.hadEnergyInHE = rawcalojet->hadEnergyInHE();
            lCalo.hadEnergyInHF = rawcalojet->hadEnergyInHF();
            lCalo.emEnergyInEB = rawcalojet->emEnergyInEB();
            lCalo.emEnergyInEE = rawcalojet->emEnergyInEE();
            lCalo.emEnergyInHF = rawcalojet->emEnergyInHF();
            lCalo.towersArea = rawcalojet->towersArea();
            lCalo.n90 = rawcalojet->n90();
            lCalo.n60 = rawcalojet->n60();

            lJet.etaMean = rawcalojet->etaPhiStatistics().etaMean;
            lJet.phiMean = rawcalojet->etaPhiStatistics().phiMean;
            lJet.etaEtaMoment = rawcalojet->etaPhiStatistics().etaEtaMoment;
            lJet.phiPhiMoment = rawcalojet->etaPhiStatistics().phiPhiMoment;
            lJet.etaPhiMoment = rawcalojet->etaPhiStatistics().etaPhiMoment;

            HbbAnalysis::Jet lObj(lFourVec,lJet,lCalo,lJPT,lJPTPF,lBtag,lId,lEnergyCorrections);
      lObj.secVertexVars(jetSecVertexVars);
            aVec.push_back(lObj);
          }
          if (isPF) {
            HbbAnalysis::PFJetVars lPF;
            lPF.chargedMuEnergy = (*iter).chargedMuEnergy();
            lPF.chargedMuEnergyFraction = (*iter).chargedMuEnergyFraction();
            lPF.neutralMultiplicity = (*iter).neutralMultiplicity();
      lPF.numberOfDaughters = (*iter).numberOfDaughters();
      lPF.HFHadronEnergy = (*iter).HFHadronEnergy();

            HbbAnalysis::Jet lObj(lFourVec,lJet,lPF,lJPTPF,lBtag,lEnergyCorrections);
      lObj.secVertexVars(jetSecVertexVars);
            aVec.push_back(lObj);
          }

        }//loop on element
    }//non empty
  }//valid

}//HbbJets

void HbbTreeMaker::HbbMet(const edm::Handle<std::vector<pat::MET> > & aCol,
                          HbbAnalysis::Met & aMet)
{//HbbMet
  if (aCol.isValid()) {
    HbbAnalysis::MetVars lGen;
    HbbAnalysis::MetVars lReco;

    assert(aCol->size() == 1);
    const pat::MET & lMet = *(aCol->begin());

    lReco.mET = lMet.pt();
    lReco.mEx = lMet.px();
    lReco.mEy = lMet.py();
    lReco.sumET = lMet.sumEt();
    lReco.phi = lMet.phi();
    lReco.mEtSig = lMet.mEtSig();
    
    const reco::GenMET *lGenMET = lMet.genMET();
    
    if (!processData_ && lGenMET){
      lGen.mET = lGenMET->pt();
      lGen.mEx = lGenMET->px();
      lGen.mEy = lGenMET->py();
      lGen.sumET = lGenMET->sumEt();
      lGen.phi = lGenMET->phi();
      lGen.mEtSig = lGenMET->mEtSig();
    }
    else {
      lGen.mET = 0;
      lGen.mEx = 0;
      lGen.mEy = 0;
      lGen.sumET = 0;
      lGen.phi = 0;
      lGen.mEtSig = 0;
    }
    
    aMet.genVars(lGen);
    aMet.recoVars(lReco);
  }
}//HbbMet

void HbbTreeMaker::HbbTrigger(const edm::Handle<edm::TriggerResults> & aCol,
                              const edm::TriggerNames & aNames,
                              std::vector<HbbAnalysis::Trigger> & aVec, int prescaleSet)
{//HbbTrigger

  if (aCol.isValid()){
    //edm::TriggerNames lNames;
    //lNames.init(*aCol);

    //for (unsigned int j=0; j<hltPaths_.size(); j++) {
    //bool valid=false;
    for (unsigned int i=0; i<aCol->size(); ++i){
      std::string trueName = aNames.triggerName(i);//.at(i);
      if (aCol->accept(i)){
        HbbAnalysis::TriggerVars lTrigVars;
        lTrigVars.name = trueName;
        lTrigVars.index = i;
        lTrigVars.accept = aCol->accept(i);
        lTrigVars.prescale = hltConfig_.prescaleValue(prescaleSet, trueName);
        HbbAnalysis::Trigger lObj(lTrigVars);
        aVec.push_back(lObj);
      }//If trigger path is passed
    }//Loop through trigger paths
  }

}

void HbbTreeMaker::HbbParticles(const edm::Handle<reco::GenParticleCollection> & aCol,
                                std::vector<HbbAnalysis::GenParticle> & aVec, HbbAnalysis::P4Total & p4total)
{//genparticles


  const reco::GenParticle *ptrToFirst  = 0;
  if (aCol->size() > 0) ptrToFirst = &((*aCol)[0]);
  //add genParticle inside vector
  p4total.comE = 0.0;
  p4total.sumPx = 0.0;
  p4total.sumPy = 0.0;
  p4total.sumPz = 0.0;
  p4total.sumE = 0.0;
  for (unsigned int mccount = 0;mccount<aCol->size();++mccount)
    {//loop on particles
      const reco::GenParticle & p = (*aCol)[mccount];

      if (p.status() == 3 && p.pdgId() == 2212){
        p4total.comE += p.energy();
      }
      if (p.status() == 1){
        p4total.sumPx += p.px();
        p4total.sumPy += p.py();
        p4total.sumPz += p.pz();
        p4total.sumE += p.energy();
      }

      bool keep = false;
      if (p.status() == 1){
        for (unsigned i = 0; i < status1ToKeepRegex_.size(); ++i){
          if (boost::regex_match(boost::lexical_cast<std::string>(p.pdgId()),status1ToKeepRegex_[i]) && p.pt() > 5.0) keep = true;
          if (abs(p.pdgId()) <= 16) keep = true; //Keep all status 1 leptons, regardless of matches passed through config
        }
      }
      if (p.status() == 2){
        for (unsigned i = 0; i < status2ToKeepRegex_.size(); ++i){
          if (boost::regex_match(boost::lexical_cast<std::string>(p.pdgId()),status2ToKeepRegex_[i])) keep = true;
        }
      }
      if (p.status() == 3){
        for (unsigned i = 0; i < status3ToKeepRegex_.size(); ++i){
          if (boost::regex_match(boost::lexical_cast<std::string>(p.pdgId()),status3ToKeepRegex_[i])) keep = true;
        }
      }

      if (keep){

        HbbAnalysis::FourMomentum lGenFourVec;
        HbbAnalysis::GenVars lGen;
        lGenFourVec.px = p.px();
        lGenFourVec.py = p.py();
        lGenFourVec.pz = p.pz();
        lGenFourVec.e = p.energy();
        lGen.charge = p.charge();
        lGen.pdgId = p.pdgId();
        lGen.status = p.status();
        lGen.index = mccount;
        lGen.vertex.vx = p.vx();
        lGen.vertex.vy = p.vy();
        lGen.vertex.vz = p.vz();
              
        HbbAnalysis::GenParticle lObj(lGenFourVec,lGen);

        for (unsigned i = 0; i < p.numberOfMothers(); ++i){
          lObj.setParent((reco::GenParticle*)p.mother(i) - ptrToFirst);
          //std::cout << ((reco::GenParticle*)p.mother(i) - ptrToFirst) << " " << p.mother(i) << "\t";
        }
        //std::cout << std::endl;
        //std::cout << "Daughters: ";
        for (unsigned i = 0; i < p.numberOfDaughters(); ++i){
          lObj.setChild((reco::GenParticle*)p.daughter(i) - ptrToFirst);
          //std::cout << ((reco::GenParticle*)p.daughter(i) - ptrToFirst) << " " << p.daughter(i) << "\t";
        }
                aVec.push_back(lObj);
          //if (mccount < 50) std::cout << "------Particle " << mccount << std::endl;
          //if (mccount < 50) lObj.Print();
        //std::cout << std::endl;
      }
/*
      unsigned int nMum = p.numberOfMothers();
      //std::cout << mccount << " index = " << const_cast<reco::GenParticle *>(&p)-lBegin << " has " << nMum << " mother(s)." << std::endl;
 
      if ( (!doGen_ && 
            (p.status()==3 || 
             (p.status()==2 && (abs(p.pdgId())==4 || abs(p.pdgId())==5)) || 
             (p.status()==1 && ( (p.pt()>5 && abs(p.pdgId()) > 16) || abs(p.pdgId()) <= 16) ) 
             )) || doGen_){
        for(unsigned int mum=0;mum<nMum;++mum){//loop on mothers
          const reco::Candidate * pMother = p.mother(mum);
          //std::cout << "-- mother index = " << pMother-lBegin << std::endl;
          for (unsigned int mccount2 = 0;mccount2<mccount;++mccount2)
            {//loop on particles
              const reco::Candidate & p2 = (*aCol)[mccount2];
              if (pMother == &p2) {
                lObj.setParent(mccount2);
                break;
              }
            }
        }//loop on mothers
  */
    }//loop on particles


  if (debug_>1){
    for (unsigned int imc(0); imc < aVec.size(); ++imc){
      HbbAnalysis::GenParticle & p = aVec.at(imc);
      p.Print();
    }
  }


}//genparticles

void HbbTreeMaker::HbbGenJets(const edm::Handle<reco::GenJetCollection> & aCol,
                              std::vector<HbbAnalysis::GenJet> & aVec,
                              const bool aSaveConstituants)
{

  if (aCol.isValid()){//valid
    if (aCol->size() > 0) {//non empty
        
      unsigned int iEle = 0;
      for (std::vector<reco::GenJet>::const_iterator iter = aCol->begin(); 
           iter != aCol->end(); 
           ++iter,++iEle) 
        {//loop on element
          //if (debug_ > 1) std::cout << "**** Ele #" << iEle << ", pT,eta,phi =  " << (*iter).pt() << " " << (*iter).eta() << " " << (*iter).phi() << std::endl;

    HbbAnalysis::FourMomentum lGenJetFourVec;
    HbbAnalysis::GenJetVars lGenJet;
    lGenJetFourVec.px = (*iter).px();
    lGenJetFourVec.py = (*iter).py();
    lGenJetFourVec.pz = (*iter).pz();
    lGenJetFourVec.e = (*iter).energy();
    lGenJet.mass = (*iter).mass();
    lGenJet.nConstituents = (*iter).getGenConstituents().size();
    lGenJet.flavour = 0;  //Flavour will be set if constituents are added 


            //only status 1 saved !!!!
            //get flavour from list of constituants.
            //std::vector <const reco::GenParticle*> lConst = (*iter).getGenConstituents();
            //unsigned int nb=0, nc=0;
            // for (unsigned int iMC(0); iMC<lConst.size(); ++iMC){
//            const reco::GenParticle* lPart = lConst[iMC];
//            //std::cout << lPart->pdgId() << " " << fabs(lPart->pdgId()) << " " << abs(lPart->pdgId()) << std::endl;
//            if (fabs(lPart->pdgId()) == 5){
//              ++nb;
//              break;
//            }
//            if (fabs(lPart->pdgId()) == 4) ++nc;
//          }
//          if (nb > 0) lGenJet.pdgId = 5;
//          else if (nc > 0) lGenJet.pdgId = 4;
//          else lGenJet.pdgId = 3;

          HbbAnalysis::GenJet lObj(lGenJetFourVec, lGenJet);
          if (aSaveConstituants){
            for (unsigned int iC(0); iC<(*iter).getGenConstituents().size(); ++iC){
              const reco::GenParticle* p = (*iter).getGenConstituent(iC);
        HbbAnalysis::FourMomentum lGenFourVec;
        HbbAnalysis::GenVars lGen;
        lGenFourVec.px = p->px();
        lGenFourVec.py = p->py();
        lGenFourVec.pz = p->pz();
        lGenFourVec.e = p->energy();
        lGen.charge = p->charge();
        lGen.pdgId = p->pdgId();
        lGen.status = p->status();
        lGen.index = -1; //We don't know this now
        lGen.vertex.vx = p->vx();
        lGen.vertex.vy = p->vy();
        lGen.vertex.vz = p->vz();
              lObj.addConstituent(HbbAnalysis::GenParticle(lGenFourVec,lGen));
            }
          }
    //std::cout << "--------------------------------------------" << std::endl;
    //std::cout << "GenJet: " << std::endl;
    //std::cout << "--------------------------------------------" << std::endl;
    //lObj.Print();
          aVec.push_back(lObj);
        }
    }
  }
}

void HbbTreeMaker::HbbVertices(const edm::Handle<std::vector<reco::Vertex> > & aCol,
                               std::vector<HbbAnalysis::Vertex> & aVec)
{

  if (aCol.isValid()) {
    for (std::vector<reco::Vertex>::const_iterator iter = aCol->begin(); 
         iter != aCol->end(); 
         ++iter) 
      {

        if (!((*iter).isValid()) || (*iter).isFake()) continue;

        HbbAnalysis::VertexVars lVtx;
  double sumPtSqTmp = 0;

        if ((*iter).tracksSize() > 0) {
          for (std::vector<reco::TrackBaseRef>::const_iterator lTrk = (*iter).tracks_begin();
               lTrk != (*iter).tracks_end();
               ++lTrk) {
      sumPtSqTmp += ( ((*(*lTrk)).pt()) * ((*(*lTrk)).pt()) );
          }
        }

  lVtx.sumPtSq = sumPtSqTmp;
        lVtx.chi2 = (*iter).chi2();
        lVtx.ndof = (*iter).ndof();
        lVtx.x = (*iter).x();
        lVtx.y = (*iter).y();
        lVtx.z = (*iter).z();
        lVtx.rho = (*iter).position().Rho();
        lVtx.xError = (*iter).xError();
        lVtx.yError = (*iter).yError();
        lVtx.zError = (*iter).zError();
        lVtx.cov01 = (*iter).covariance(0,1);
        lVtx.cov02 = (*iter).covariance(0,2);
        lVtx.cov12 = (*iter).covariance(1,2);

  //AJG
  /*
  std::cout << "Vertex (" << (*iter).x() << ","
                          << (*iter).y() << ","
                          << (*iter).z() << ")" << std::endl;
  */

        aVec.push_back(HbbAnalysis::Vertex(lVtx));
      }
  }

}

/*
void HbbTreeMaker::HbbL1Objects(const edm::Handle<l1extra::L1JetParticleCollection> & aCol,
                                std::vector<HbbAnalysis::L1Object> & aVec,
                                const unsigned int aType)
{

  for (l1extra::L1JetParticleCollection::const_iterator iter = aCol->begin(); 
         iter != aCol->end(); 
         ++iter) 
    {
        
      HbbAnalysis::L1Vars lVars;
      lVars.pT = (*iter).pt();
      lVars.ET = (*iter).et();
      lVars.eta = (*iter).eta();
      lVars.phi = (*iter).phi();
      lVars.bx = (*iter).bx();
      lVars.type = aType;

      HbbAnalysis::L1Object lObj(lVars);
      aVec.push_back(lObj);
        
    }

}


void HbbTreeMaker::HbbHLTObjects(const edm::Handle<trigger::TriggerEvent> & aCol,
                                 std::vector<HbbAnalysis::HLTObject> & aVec)
{


  //cout
  //const trigger::TriggerObjectCollection & toc(aCol->getObjects());
  //for (unsigned i=0; i<aCol->sizeFilters(); i++){
  //const edm::InputTag & lTag = aCol->filterTag(i);
    //std::cout << " -- ele " << i << ": " << lTag << ", key ids = " ;
    ////const int index = aCol->filterIndex(lTag);
    ////if (index >= aCol->sizeFilters()) {
      ////std::cout << " -- Filter " << aTag << " not found !" << std::endl;
      ////continue;
    ////}
    //const trigger::Keys & k = aCol->filterKeys(i);
    //for (trigger::Keys::const_iterator ki = k.begin(); ki !=k.end(); ++ki ) {
    //std::cout <<  toc[*ki].id() << ", " ;
    //}
    //std::cout << std::endl;
  //}

}

void HbbTreeMaker::fillHLTVector(const edm::InputTag & aTag,
                                 const edm::Handle<trigger::TriggerEvent> & aCol,
                                 std::vector<HbbAnalysis::HLTObject> & aVec,
                                 const unsigned int aType){

  const trigger::TriggerObjectCollection & toc(aCol->getObjects());
  const int index = aCol->filterIndex(aTag);
  if (index >= aCol->sizeFilters()) {
    //std::cout << " -- Filter " << aTag << " not found !" << std::endl;
    return;
  }

  const trigger::Keys & k = aCol->filterKeys(index);
         
  for (trigger::Keys::const_iterator ki = k.begin(); ki !=k.end(); ++ki ) {
    HbbAnalysis::HLTVars lVars;
    lVars.pT = toc[*ki].pt();
    lVars.eta = toc[*ki].eta();
    lVars.phi = toc[*ki].phi();
    lVars.id = toc[*ki].id();
    lVars.mass = toc[*ki].mass();
    lVars.type = aType;
    HbbAnalysis::HLTObject lObj(lVars);
    aVec.push_back(lObj);           
  }
}
*/
void HbbTreeMaker::HbbGenInfo(const edm::Handle<edm::HepMCProduct> & amcProd,
                              const edm::Handle<GenRunInfoProduct> & aRunProd,
            const edm::Handle<GenEventInfoProduct> & aInfoProd)
{

  HbbAnalysis::GenInfo lInfo;

  if (amcProd.isValid()) {
    const HepMC::GenEvent * genEvt = amcProd->GetEvent();
    lInfo.processID(genEvt->signal_process_id());
    lInfo.ptHat(genEvt->event_scale());
  }

  if (aRunProd.isValid()) {

    lInfo.internalXS(aRunProd->internalXSec().value(),
        aRunProd->internalXSec().error());

    lInfo.externalXSLO(aRunProd->externalXSecLO().value(),
        aRunProd->externalXSecLO().error());

    lInfo.externalXSNLO(aRunProd->externalXSecNLO().value(),
        aRunProd->externalXSecNLO().error());

    lInfo.filterEff(aRunProd->filterEfficiency());
  }

  if (aInfoProd.isValid()) {
    lInfo.weight(aInfoProd->weight());
  }
    
  event_->generatorInfo(lInfo);

}

void HbbTreeMaker::HbbPUInfo(const edm::Handle<std::vector<PileupSummaryInfo> > & aPUInfo,
    std::vector<HbbAnalysis::PUVars> & aVec) {
  std::vector<PileupSummaryInfo>::const_iterator PVI;
  for (PVI = aPUInfo->begin(); PVI != aPUInfo->end(); ++PVI) {
    HbbAnalysis::PUVars aVars;
    aVars.bunchCrossing = PVI->getBunchCrossing();
    aVars.numInteractions = PVI->getPU_NumInteractions();
    aVec.push_back(aVars);
  }
}

void HbbTreeMaker::HbbLHEInfo(const edm::Handle<LHEEventProduct> & aLHEEvt,
                              std::vector<HbbAnalysis::GenParticle> & aVec) {

  const lhef::HEPEUP lhepeup = aLHEEvt->hepeup();
 
  //4-momentum (px,py,pz,E)+mass
  const std::vector<lhef::HEPEUP::FiveVector> lpup = lhepeup.PUP;
  //pdgid int
  const std::vector<int> lidup = lhepeup.IDUP;
  //status int
  const std::vector<int> listup = lhepeup.ISTUP;
  //mother: pair<int,int>
  const std::vector<std::pair<int,int> > lmothup = lhepeup.MOTHUP;


  //add genParticle inside vector
  for (unsigned int imc = 0;imc<lpup.size();++imc)
    {//loop on particles
      
      HbbAnalysis::FourMomentum lFourVec;
      HbbAnalysis::GenVars lMC;
      lMC.index = imc;

      TLorentzVector lVec((lpup[imc])[0],(lpup[imc])[1],(lpup[imc])[2],(lpup[imc])[3]);

      lFourVec.e = lVec.E();
      lFourVec.px = lVec.Px();
      lFourVec.py = lVec.Py();
      lFourVec.pz = lVec.Pz();
      /*
      // Special case if pT = 0 (otherwise ROOT complains)
      if (fabs(0.0 - lVec.Pt()) < std::numeric_limits<double>::epsilon()){
        if (lVec.Pz() >= 0) {
          lMC.eta = 1E10;
          lMC.y = 1E10;
        } else {
          lMC.eta = -1E10;
          lMC.y = -1E10;
        }
      } else {
      lMC.eta = lVec.Eta();
      lMC.y = lVec.Rapidity();
      }
      lMC.phi = lVec.Phi();
      */
      lMC.pdgId = lidup[imc];
      lMC.status = listup[imc];

      HbbAnalysis::GenParticle lObj(lFourVec,lMC);
      lObj.setParent(lmothup[imc].first);
      lObj.setParent(lmothup[imc].second);
      aVec.push_back(lObj);

    }//loop on particles
  
  if (debug_>1){
    for (unsigned int imc(0); imc < aVec.size(); ++imc){
      //HbbAnalysis::GenParticle & p = aVec.at(imc);
      //p.print();
    }
  }
}

#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(HbbTreeMaker);


Revision:	1.48
Committed:	Wed Jan 25 18:45:49 2012 UTC (13 years, 3 months ago) by agilbert
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.47:	+53 -18 lines
Log Message:	Lots of updates