HbbAnalysis/plugins/HbbTreeMaker.cc

#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/MuonReco/interface/Muon.h"
#include "DataFormats/MuonReco/interface/MuonFwd.h"
#include "DataFormats/MuonReco/interface/MuonSelectors.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/HLTReco/interface/TriggerEvent.h"
#include "DataFormats/TrackReco/interface/TrackFwd.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/MET.h"

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

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

#include "UserCode/HbbAnalysis/interface/Electron.hh"
#include "UserCode/HbbAnalysis/interface/Muon.hh"
#include "UserCode/HbbAnalysis/interface/Tau.hh"
#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"
#include "UserCode/HbbAnalysis/interface/HLTObject.hh"

#include "UserCode/HbbAnalysis/plugins/HbbTreeMaker.hh"

using namespace HbbAnalysis;

HbbTreeMaker::HbbTreeMaker(const edm::ParameterSet & pset):
  debug_(pset.getParameter<int>("DEBUG")),
  processData_(false),
  flavour_(pset.getParameter<unsigned int>("JetFlavour")),
  doGen_(pset.getParameter<bool>("DOGEN")),
  genParticleSrc_(pset.getParameter<edm::InputTag>("GenParticles")),
  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")),
  pfJetSrc_(pset.getParameter<edm::InputTag>("PFJets")),
  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")),
  event_(0),
  tree_(0),
  hltTags_(pset.getParameter<std::vector<edm::InputTag> >("HLTTags"))
  //processVec_(pset.getParameter<std::vector<std::string> >("ProcessVec"))
{//constructor


}//constructor

HbbTreeMaker::~HbbTreeMaker(){//destructor
}//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 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::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;

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

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


  if (!processData_) {
    edm::Handle<edm::HepMCProduct> mcProd;
    edm::Handle<GenRunInfoProduct> lRunProd;
    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;
    }

    HbbGenInfo(mcProd,lRunProd);
  }

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

  if (doGen_ && !processData_) HbbParticles(lGenParticles,event_->particles());

  unsigned int lNPartons = 0;
  if (!processData_) lNPartons = jetFlav_.fillPartons(lGenParticles);
  
  if (debug_ && !processData_) 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());


  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,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;
  }

  //std::cout << "Processing calo jets:" << std::endl;
  HbbJets(lCaloJetCollection,jetFlav_,lGenParticles,event_->caloJets());
  
  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> > 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::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);
    HbbTrigger(lTrigCol,lNames,event_->triggers());
  }

  //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_->hltJets());
  }


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

}//analyze


void HbbTreeMaker::HbbElectrons(const edm::Handle<std::vector<pat::Electron> > & aCol,
                                std::vector<HbbAnalysis::Electron> & aVec)
{//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++) 
        {
          //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::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();
          
          HbbAnalysis::SCVars lSC;
          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.calo = (*iter).caloIso();
          lIso.track = (*iter).trackIso();
          lIso.ecal = (*iter).ecalIso();
          lIso.hcal = (*iter).hcalIso();

          HbbAnalysis::EleIDVars lID;
          lID.electronIDs = (*iter).electronIDs();
          lID.hOverE = (*iter).hadronicOverEm();
          lID.deltaPhiIn = (*iter).deltaPhiSuperClusterTrackAtVtx();
          lID.deltaEtaIn = (*iter).deltaEtaSuperClusterTrackAtVtx();
          lID.ecalDrivenSeed = (*iter).ecalDrivenSeed();
          lID.trackerDrivenSeed = (*iter).trackerDrivenSeed();

          HbbAnalysis::Electron lObj(lGen,lReco,lSC,lIso,lID);
          aVec.push_back(lObj);
          iEle++;
        }
    }
  }

}//HbbElectrons


void HbbTreeMaker::HbbMuons(const edm::Handle<std::vector<pat::Muon> > & aCol,
                            const edm::Handle<std::vector<reco::Vertex> > & aRecoVertices,
                            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++) 
        {
          const reco::Muon* recoMuon = dynamic_cast<const reco::Muon*>((*iter).originalObject());

          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::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();

          HbbAnalysis::MuTrkVars lTrk;

          reco::TrackRef lTrackerTrk = (*iter).innerTrack();
          if ( lTrackerTrk.isAvailable() && lTrackerTrk.isNonnull() ) {
            if ( aRecoVertices->size() >= 1 ) {
              const reco::Vertex& thePrimaryEventVertex = (*aRecoVertices->begin());
              lTrk.IPd0 = -lTrackerTrk->dxy(thePrimaryEventVertex.position());
              lTrk.IPdz = lTrackerTrk->dz(thePrimaryEventVertex.position());
            }
            else {
              lTrk.IPd0 = 0;
              lTrk.IPdz = 0;
            }

            lTrk.nHits =  lTrackerTrk->numberOfValidHits();
          }
          else {
            lTrk.IPd0 = 0;
            lTrk.IPdz = 0;
            lTrk.nHits = 0;
          }

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

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

          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;

          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);

          HbbAnalysis::Muon lObj(lGen,lReco,lTrk,lIsoR03,lIsoR05,lID);
          aVec.push_back(lObj);
          iEle++;
        }
    }
  }

}//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++) 
        {
          //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);
          }
          iEle++;
        }
    }
  }

}//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++) 
        {//loop on element
          //if (debug_ > 1) std::cout << "**** Ele #" << iEle << ", pT,eta,phi =  " << (*iter).pt() << " " << (*iter).eta() << " " << (*iter).phi() << std::endl;

          HbbAnalysis::GenVars lGen;
          if (!processData_ && (*iter).genParton()){
            lGen.valid = true;
            lGen.E = (*iter).genParton()->energy();
            lGen.pT = (*iter).genParton()->pt();
            lGen.eta = (*iter).genParton()->eta();
            lGen.y = (*iter).genParton()->y();
            lGen.phi = (*iter).genParton()->phi();
            lGen.charge = (*iter).genParton()->charge();
            lGen.pdgId = (*iter).genParton()->pdgId();
            lGen.status = (*iter).genParton()->status();
            lGen.mass = (*iter).genParton()->mass();
            lGen.vx = (*iter).genParton()->vx();
            lGen.vy = (*iter).genParton()->vy();
            lGen.vz = (*iter).genParton()->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.etaDet = (*iter).detectorEta((*iter).vz(),(*iter).eta());
          lReco.y = (*iter).y();
          lReco.phi = (*iter).phi();
          lReco.charge = (*iter).jetCharge();
          lReco.vx = (*iter).vx();
          lReco.vy = (*iter).vy();
          lReco.vz = (*iter).vz();
          //lReco.myEtaDet = HbbAnalysis::EtaDetector(lReco);

          unsigned int lFlav = 0;
          if (!processData_ && aJetFlav.nPartons()) {
            lFlav = aJetFlav.partonMatchingGenJet((*iter),aGenParticles,0.4).second;
          }

          HbbAnalysis::JetVars lCommon;
          lCommon.flavour = lFlav;
          lCommon.partonFlavour = (*iter).partonFlavour();
          lCommon.nAssociatedTracks = ((*iter).associatedTracks()).size();
          lCommon.rawpT = (*iter).pt();
          if ((*iter).hasCorrFactors() && (*iter).corrFactor("raw")>0)
            lCommon.rawpT = (*iter).pt()/((*iter).corrFactor("raw"));

          lCommon.etaMean = (*iter).etaPhiStatistics().etaMean;
          lCommon.phiMean = (*iter).etaPhiStatistics().phiMean;
          lCommon.etaEtaMoment = (*iter).etaPhiStatistics().etaEtaMoment;
          lCommon.phiPhiMoment = (*iter).etaPhiStatistics().phiPhiMoment;
          lCommon.etaPhiMoment = (*iter).etaPhiStatistics().etaPhiMoment;
          
          lCommon.l1Match = false;//TODO
          lCommon.hltMatch = false;//TODO

          //lCommon.rawEta = (*iter).;
          //lCommon.rawPhi = (*iter).;
          //    p_hasJetCorrFactors->Fill(aJet.hasJetCorrFactors());

          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") < -1000) ? -1000 : (*iter).bDiscriminator("softElectronByPtBJetTags");
          lBtag.softElectronByIP3d = ((*iter).bDiscriminator("softElectronByIP3dBJetTags") < -1000) ? -1000 : (*iter).bDiscriminator("softElectronByIP3dBJetTags");
          lBtag.softMuon = ((*iter).bDiscriminator("softMuonBJetTags") < -1000) ? -1000 : (*iter).bDiscriminator("softMuonBJetTags");
          lBtag.softMuonByPt = ((*iter).bDiscriminator("softMuonByPtBJetTags") < -1000) ? -1000 : (*iter).bDiscriminator("softMuonByPtBJetTags");
          lBtag.softMuonByIP3d = ((*iter).bDiscriminator("softMuonByIP3dBJetTags") < -1000) ? -1000 : (*iter).bDiscriminator("softMuonByIP3dBJetTags");

          //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;


          lBtag.tCHE = (*iter).bDiscriminator("trackCountingHighEffBJetTags");
          lBtag.tCHP = (*iter).bDiscriminator("trackCountingHighPurBJetTags");


          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;
          lId.fHPD = 0;
          lId.fRBX = 0;
          lId.n90Hits = 0;
          lId.hitsInN90 = 0;
          lId.restrictedEMF = 0;
          lId.nHCALTowers = 0;
          lId.nECALTowers = 0;
          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(lGen,lGenJet,lReco,lCommon,lCalo,lBtag,lId);
            aVec.push_back(lObj);
          }
          HbbAnalysis::JPTPFJetVars lJPTPF;
          lJPTPF.chargedHadronEnergy = 0;
          lJPTPF.chargedHadronEnergyFraction = 0;
          lJPTPF.neutralHadronEnergy = 0;
          lJPTPF.neutralHadronEnergyFraction = 0;
          lJPTPF.chargedEmEnergy = 0;
          lJPTPF.chargedEmEnergyFraction = 0;
          lJPTPF.neutralEmEnergy = 0;
          lJPTPF.neutralEmEnergyFraction = 0;
          lJPTPF.chargedMultiplicity = 0;
          lJPTPF.muonMultiplicity = 0;

          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;
//          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);
                                   
            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();

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

            HbbAnalysis::Jet lObj(lGen,lGenJet,lReco,lCommon,lCalo,lJPT,lJPTPF,lBtag,lId);
            aVec.push_back(lObj);
          }
          if (isPF) {
            HbbAnalysis::PFJetVars lPF;
            lPF.chargedMuEnergy = (*iter).chargedMuEnergy();
            lPF.chargedMuEnergyFraction = (*iter).chargedMuEnergyFraction();
            lPF.neutralMultiplicity = (*iter).neutralMultiplicity();

            HbbAnalysis::Jet lObj(lGen,lGenJet,lReco,lCommon,lPF,lJPTPF,lBtag);
            aVec.push_back(lObj);
          }

          iEle++;
        }//loop on element
    }//non empty
  }//valid

}//HbbJets

void HbbTreeMaker::HbbMet(const edm::Handle<std::vector<pat::MET> > & aCol,
                          HbbAnalysis::Met & aMet)
{//HbbMet
  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)
{//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);

      HbbAnalysis::TriggerVars lTrigVars;
      lTrigVars.name = trueName;
      lTrigVars.index = i;
      lTrigVars.accept = aCol->accept(i);

      HbbAnalysis::Trigger lObj(lTrigVars);

      aVec.push_back(lObj);
    }
    
  }

}

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

  //add genParticle inside vector
  for (unsigned int mccount = 0;mccount<aCol->size();mccount++)
    {//loop on particles
      const reco::Candidate & p = (*aCol)[mccount];

      HbbAnalysis::MCVars lMC;
      lMC.index = mccount;
      lMC.E = p.energy();
      lMC.pT = p.pt();
      lMC.eta = p.eta();
      lMC.y = p.y();
      lMC.phi = p.phi();
      lMC.pdgId = p.pdgId();
      lMC.status = p.status();
      HbbAnalysis::GenParticle lObj(lMC);
      aVec.push_back(lObj);
      
    }//loop on particles

  //now need to get the list of parents....

  for (unsigned int mccount = 0;mccount<aCol->size();mccount++)
    {//loop on particles
      const reco::Candidate & p = (*aCol)[mccount];
 
      unsigned int nD = p.numberOfDaughters();
      //std::cout << mccount << " has " << nD << " daughter(s) : " ;
      for(unsigned int dau=0;dau<nD;dau++){//loop on daughters
        const reco::Candidate * pDau = p.daughter(dau);
        //std::cout << pDau << " ("<<std::flush;
        for (unsigned int mccount2 = 0;mccount2<aCol->size();mccount2++)
          {//loop on particles
            const reco::Candidate & p2 = (*aCol)[mccount2];
            //std::cout<<"DBG: mccount2 = "<<mccount2<<"  gen  size = "<<data_->gen().size()<<std::endl;
            
            HbbAnalysis::GenParticle & part2 =  aVec.at(mccount2);
            if (pDau == &p2) {
              part2.setParent(mccount);
              //std::cout << &p2 << ", index = " << mccount2 << "), "<<std::flush;
              break;
            }
            //else std::cout << std::endl << "****no match " << mccount2 << " " << &p2 << std::endl;
            //if (p2.pdgId() == pDau->pdgId() && fabs(p2.energy()-pDau->energy())<0.000001 && fabs(p2.eta()-pDau->eta())<0.000001 && fabs(p2.phi()-pDau->phi())<0.000001){
            //part2.parent(mccount);
            //}
          }//loop on particles
      }//loop on daughters
      //std::cout << std::endl;

    }//loop on particles

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


}//genparticles


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

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

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

      HbbAnalysis::VertexVars lVtx;

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

      if (lVtx.trackWeights.size() != (*iter).tracksSize()) {
        std::cout<< " -- Problem with tracks, size is not as expected ! "
                 << std::endl
                 << " --- Size of recoVertex tracks = " << (*iter).tracksSize()
                 << std::endl
                 << " --- Size of trackWeights = " << lVtx.trackWeights.size()
                 << std::endl;
      }

      lVtx.chi2 = (*iter).chi2();
      lVtx.ndof = (*iter).ndof();
      lVtx.x = (*iter).x();
      lVtx.y = (*iter).y();
      lVtx.z = (*iter).z();
      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);

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

}

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)
{

  const trigger::TriggerObjectCollection & toc(aCol->getObjects());

  //cout
  //for (unsigned i=0; i<aCol->sizeFilters(); i++){
  //const edm::InputTag & lTag = aCol->filterTag(i);
  //std::cout << "ele " << i << ": " << lTag << std::endl;
  //}

  for (unsigned i(0); i<hltTags_.size(); i++){

    const int index = aCol->filterIndex(hltTags_[i]);
    if (index >= aCol->sizeFilters()) {
      //std::cout << " -- Filter " << hltTags_[i] << " not found !" << std::endl;
      continue;
    }

    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();
      HbbAnalysis::HLTObject lObj(lVars);
      aVec.push_back(lObj);           
    }
  }

}

void HbbTreeMaker::HbbGenInfo(const edm::Handle<edm::HepMCProduct> & amcProd,
                              const edm::Handle<GenRunInfoProduct> & aRunProd) 
{
  const HepMC::GenEvent * genEvt = amcProd->GetEvent();
  HbbAnalysis::GenInfo lInfo;

  lInfo.processID(genEvt->signal_process_id());
  lInfo.ptHat(genEvt->event_scale());

  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());
    
  event_->generatorInfo(lInfo);

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


Revision:	1.19
Committed:	Thu Jun 24 13:41:56 2010 UTC (14 years, 10 months ago) by amagnan
Content type:	text/plain
Branch:	MAIN
Changes since 1.18:	+47 -1 lines
Log Message:	add generator information to tree