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/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 "FWCore/ServiceRegistry/interface/Service.h"
#include "PhysicsTools/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/plugins/HbbTreeMaker.hh"

using namespace HbbAnalysis;

HbbTreeMaker::HbbTreeMaker(const edm::ParameterSet & pset):
  debug_(pset.getParameter<int>("DEBUG")),
  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")),
  event_(0),
  tree_(0)
  //processVec_(pset.getParameter<std::vector<std::string> >("ProcessVec"))
{//constructor


}//constructor

HbbTreeMaker::~HbbTreeMaker(){//destructor
}//destructor


void HbbTreeMaker::beginJob(const edm::EventSetup&){//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
  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());

  edm::Handle<reco::GenParticleCollection> lGenParticles;
  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());

  unsigned int 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;
  }

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


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

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

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

  HbbTrigger(lTrigCol,event_->triggers());

  //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 ((*iter).genLepton()){
            lGen.valid = true;
            lGen.E = (*iter).genLepton()->energy();
            lGen.pT = (*iter).genLepton()->pt();
            lGen.eta = (*iter).genLepton()->eta();
            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.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.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();

          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 ((*iter).genLepton()){
            lGen.valid = true;
            lGen.E = (*iter).genLepton()->energy();
            lGen.pT = (*iter).genLepton()->pt();
            lGen.eta = (*iter).genLepton()->eta();
            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.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.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 (recoMuon->isGood(reco::Muon::AllGlobalMuons)) lID.ids.push_back(1);
          if (recoMuon->isGood(reco::Muon::AllStandAloneMuons)) lID.ids.push_back(2);
          if (recoMuon->isGood(reco::Muon::AllTrackerMuons)) lID.ids.push_back(3);
          if (recoMuon->isGood(reco::Muon::TrackerMuonArbitrated)) lID.ids.push_back(4);
          if (recoMuon->isGood(reco::Muon::AllArbitrated)) lID.ids.push_back(5);
          if (recoMuon->isGood(reco::Muon::GlobalMuonPromptTight)) lID.ids.push_back(6);
          if (recoMuon->isGood(reco::Muon::TMLastStationLoose)) lID.ids.push_back(7);
          if (recoMuon->isGood(reco::Muon::TMLastStationTight)) lID.ids.push_back(8);
          if (recoMuon->isGood(reco::Muon::TM2DCompatibilityLoose)) lID.ids.push_back(9);
          if (recoMuon->isGood(reco::Muon::TM2DCompatibilityTight)) lID.ids.push_back(10);
          if (recoMuon->isGood(reco::Muon::TMOneStationLoose)) lID.ids.push_back(11);
          if (recoMuon->isGood(reco::Muon::TMOneStationTight)) lID.ids.push_back(12);
          if (recoMuon->isGood(reco::Muon::TMLastStationOptimizedLowPtLoose)) lID.ids.push_back(13);
          if (recoMuon->isGood(reco::Muon::TMLastStationOptimizedLowPtTight)) lID.ids.push_back(14);

          lID.caloCompat = recoMuon->caloCompatibility();
          lID.segCompat = recoMuon->segmentCompatibility();
          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 ((*iter).genLepton()){
            lGen.valid = true;
            lGen.E = (*iter).genLepton()->energy();
            lGen.pT = (*iter).genLepton()->pt();
            lGen.eta = (*iter).genLepton()->eta();
            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.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 ((*iter).genJet()){
            lGenJet.valid = true;
            lGenJet.pT = (*iter).genJet()->pt();
            lGenJet.eta = (*iter).genJet()->eta();
            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.pT = 0;
            lGenJet.eta = 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.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::TauLeadTrkVars 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() != 7) || 
              (isCalo && lIDs.size() != 3))
            lPrint = true;

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


          HbbAnalysis::PFTauIDVars lpfId;
          HbbAnalysis::CaloTauIDVars lcaloId;

          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("trackIsolation") != lName.npos) lpfId.byTrackIsolation = lDiscri;
              if (lName.find("ecalIsolation") != lName.npos) lpfId.byECALIsolation = lDiscri;
              if (lName.find("byIsolation") != lName.npos) lpfId.byIsolation = lDiscri;
              if (lName.find("againstElectron") != lName.npos) lpfId.againstElectron = lDiscri;
              if (lName.find("againstMuon") != lName.npos) lpfId.againstMuon = 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 = 0;

            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 ((*iter).genParton()){
            lGen.valid = true;
            lGen.E = (*iter).genParton()->energy();
            lGen.pT = (*iter).genParton()->pt();
            lGen.eta = (*iter).genParton()->eta();
            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.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 ((*iter).genJet()){
            lGenJet.valid = true;
            lGenJet.pT = (*iter).genJet()->pt();
            lGenJet.eta = (*iter).genJet()->eta();
            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.pT = 0;
            lGenJet.eta = 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.phi = (*iter).phi();
          lReco.charge = (*iter).jetCharge();
          lReco.vx = (*iter).vx();
          lReco.vy = (*iter).vy();
          lReco.vz = (*iter).vz();

          unsigned int lFlav = 0;
          if (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).hasJetCorrFactors())
            lCommon.rawpT = (*iter).pt()/((*iter).jetCorrFactors().scaleDefault());
          //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.sSV = (*iter).bDiscriminator("simpleSecondaryVertexBJetTags");
          lBtag.softElectron = (*iter).bDiscriminator("softElectronBJetTags");
          lBtag.softMuon = (*iter).bDiscriminator("softMuonBJetTags");
          lBtag.softMuonNoIP = (*iter).bDiscriminator("softMuonNoIPBJetTags");
          lBtag.tCHE = (*iter).bDiscriminator("trackCountingHighEffBJetTags");
          lBtag.tCHP = (*iter).bDiscriminator("trackCountingHighPurBJetTags");
          
          bool isCalo = (*iter).isCaloJet();
          bool isPF = (*iter).isPFJet();

          assert (isCalo == !isPF);
          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);
            aVec.push_back(lObj);
          
          }
     
          if (isPF){
            HbbAnalysis::PFJetVars lPF;
            lPF.chargedHadronEnergy = (*iter).chargedHadronEnergy();
            lPF.chargedHadronEnergyFraction = (*iter).chargedHadronEnergyFraction();
            lPF.neutralHadronEnergy = (*iter).neutralHadronEnergy();
            lPF.neutralHadronEnergyFraction = (*iter).neutralHadronEnergyFraction();
            lPF.chargedEmEnergy = (*iter).chargedEmEnergy();
            lPF.chargedEmEnergyFraction = (*iter).chargedEmEnergyFraction();
            lPF.chargedMuEnergy = (*iter).chargedMuEnergy();
            lPF.chargedMuEnergyFraction = (*iter).chargedMuEnergyFraction();
            lPF.neutralEmEnergy = (*iter).neutralEmEnergy();
            lPF.neutralEmEnergyFraction = (*iter).neutralEmEnergyFraction();
            lPF.chargedMultiplicity = (*iter).chargedMultiplicity();
            lPF.neutralMultiplicity = (*iter).neutralMultiplicity();
            lPF.muonMultiplicity = (*iter).muonMultiplicity();

            HbbAnalysis::Jet lObj(lGen,lGenJet,lReco,lCommon,lPF,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 (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,
                              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 = lNames.triggerNames().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.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
Revision:	1.3
Committed:	Tue Oct 13 12:14:05 2009 UTC (15 years, 7 months ago) by amagnan
Content type:	text/plain
Branch:	MAIN
Changes since 1.2:	+28 -5 lines
Log Message:	add innerTrack IPd0, IPdz and nHits info for muons.