HbbAnalysis/plugins/MCTreeMaker.cc

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

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

#include "DataFormats/Candidate/interface/Candidate.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "DataFormats/JetReco/interface/GenJet.h"
#include "DataFormats/JetReco/interface/Jet.h"
#include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"


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

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

MCTreeMaker::MCTreeMaker(const edm::ParameterSet & pset):
  debug_(pset.getParameter<int>("DEBUG")),
  doGen_(pset.getParameter<bool>("DOGEN")),
  genParticleSrc_(pset.getParameter<edm::InputTag>("GenParticles")),
  genJetSrc_(pset.getParameter<edm::InputTag>("GenJets")),
  partonJetSrc_(pset.getParameter<edm::InputTag>("PartonJets")),
  pfJetSrc_(pset.getParameter<edm::InputTag>("PFJets")),
  status1ToKeep_(pset.getParameter<std::vector<std::string> >("Status1ToKeep")),
  status2ToKeep_(pset.getParameter<std::vector<std::string> >("Status2ToKeep")),
  status3ToKeep_(pset.getParameter<std::vector<std::string> >("Status3ToKeep")),
  event_(0),
  tree_(0)
  //processVec_(pset.getParameter<std::vector<std::string> >("ProcessVec"))
{//constructor
  for (unsigned i = 0; i < status1ToKeep_.size(); ++i){
    status1ToKeepRegex_.push_back(boost::regex(status1ToKeep_[i]));
  }
  for (unsigned i = 0; i < status2ToKeep_.size(); ++i){
    status2ToKeepRegex_.push_back(boost::regex(status2ToKeep_[i]));
  }
  for (unsigned i = 0; i < status3ToKeep_.size(); ++i){
    status3ToKeepRegex_.push_back(boost::regex(status3ToKeep_[i]));
  }
}//constructor

MCTreeMaker::~MCTreeMaker()
{//destructor

}//destructor


void MCTreeMaker::beginJob(){//beginJob


  edm::Service<TFileService> lFileService;

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

  tree_ = lDir.make<TTree>("Tree","Tree");
  tree_->Branch("MCEvent","HbbAnalysis::MCEvent",&event_);
  event_ = new HbbAnalysis::MCEvent();
  
}//beginJob


void MCTreeMaker::endJob(){//endJob

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

void MCTreeMaker::beginRun(const edm::Run& aRun, const edm::EventSetup& aEvtSetup){//beginRun

}

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

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


  static bool firstEvent = true;

  event_->Clear();
  event_->event(aEvt.id().event());

  event_->run(aEvt.run());
  //used in jet method....
  edm::Handle<reco::GenParticleCollection> lGenParticles;

  //    edm::Handle<edm::HepMCProduct> mcProd;
  //     edm::Handle<GenRunInfoProduct> lRunProd;
  //     edm::Handle<GenEventInfoProduct> genEventInfoHandle;
  
  //     try {
  //       aEvt.getByLabel("generator", mcProd );
  //     } catch(cms::Exception& e)  {
  //       std::cout << "AMM: Collection HepMCProduct not available! Exception : " << e.what() << ". " << std::endl;
  //     }
  //     try {
  //       aEvt.getRun().getByLabel("generator",lRunProd);
  //     } catch(cms::Exception& e)  {
  //       std::cout << "AMM: Collection GenInfoProduct not available! Exception : " << e.what() << ". " << std::endl;
  //     }
  //     try {
  //       aEvt.getByLabel("generator",genEventInfoHandle);
  //     } catch(cms::Exception& e)  {
  //       std::cout << "AMM: Collection GenEventInfoProduct not available! Exception : " << e.what() << ". " << std::endl;
  //     }
  
  //     HbbGenInfo(mcProd,lRunProd,genEventInfoHandle);
  
  edm::Handle<reco::GenJetCollection> lGenJets;
  try {
    aEvt.getByLabel(genJetSrc_,lGenJets);  
    if (debug_) std::cout << "** ngenJets = " << lGenJets->size() << std::endl;
  } catch(cms::Exception& e)  {
    std::cout << "AMM: Collection GenJets not available! Exception : " << e.what() << ". " << std::endl;
  }
  
  //save components
  if (lGenJets.isValid()) HbbGenJets(lGenJets,event_->genjets(),true);

  edm::Handle<reco::GenJetCollection> lPartonJets;
  try {
    aEvt.getByLabel(partonJetSrc_,lPartonJets);  
    if (debug_) std::cout << "** npartonJets = " << lPartonJets->size() << std::endl;
  } catch(cms::Exception& e)  {
    std::cout << "AMM: Collection PartonJets not available! Exception : " << e.what() << ". " << std::endl;
  }
  
  //save components
  if (lPartonJets.isValid()) HbbGenJets(lPartonJets,event_->partonjets(),true);
   
  edm::Handle<LHEEventProduct> lLHEEvt;
  try {
    aEvt.getByType(lLHEEvt);
  } catch(cms::Exception& e)  {
    std::cout << "AMM: Collection of LHEEventProduct not available! Exception : " << e.what() << ". " << std::endl;
  }
  
  if (lLHEEvt.isValid()){
    HbbLHEInfo(lLHEEvt,event_->lheParticles());
  }
    
  try {
    aEvt.getByLabel(genParticleSrc_,lGenParticles);  
    if (debug_) std::cout << "** ngenParticles = " << lGenParticles->size() << std::endl;
  } catch(cms::Exception& e)  {
    std::cout << "AMM: Collection genParticles not available! Exception : " << e.what() << ". " << std::endl;
  }
  
  //if (doGen_) HbbParticles(lGenParticles,event_->particles());
  HbbParticles(lGenParticles,event_->particles(), event_->p4total());

  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,lGenParticles,event_->pfjets());
  

  tree_->Fill();
  firstEvent = false;

}//analyze


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


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

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

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

      if (keep){

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

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


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


}//genparticles

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

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

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


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

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


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

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


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

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

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

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

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


void MCTreeMaker::HbbJets(const edm::Handle<std::vector<pat::Jet> > & aCol, 
                           const edm::Handle<reco::GenParticleCollection> & aGenParticles,
                           std::vector<HbbAnalysis::Jet> & aVec)
{//HbbJets
  
  if (aCol.isValid()){//valid
    if (aCol->size() > 0) {//non empty
      
      unsigned int iEle = 0;
      for (std::vector<pat::Jet>::const_iterator iter = aCol->begin(); 
           iter != aCol->end(); 
           ++iter,++iEle) 
        {//loop on element
          if (debug_ > 1) std::cout << "**** Ele #" << iEle << ", pT,eta,phi =  " << (*iter).pt() << " " << (*iter).eta() << " " << (*iter).phi() << std::endl;
          
          //----------------------------------------------------------------------
          HbbAnalysis::FourMomentum lFourVec;
          //----------------------------------------------------------------------
          lFourVec.e = (*iter).energy();
          lFourVec.px = (*iter).px();
          lFourVec.py = (*iter).py();
          lFourVec.pz = (*iter).pz();
          //----------------------------------------------------------------------
          HbbAnalysis::JetVars lJet;
          //----------------------------------------------------------------------
          lJet.partonFlavour = (*iter).partonFlavour();
          lJet.nAssociatedTracks = ((*iter).associatedTracks()).size();
          lJet.etaMean = (*iter).etaPhiStatistics().etaMean;
          lJet.phiMean = (*iter).etaPhiStatistics().phiMean;
          lJet.etaEtaMoment = (*iter).etaPhiStatistics().etaEtaMoment;
          lJet.phiPhiMoment = (*iter).etaPhiStatistics().phiPhiMoment;
          lJet.etaPhiMoment = (*iter).etaPhiStatistics().etaPhiMoment;
          //Construct a GenParticle if processing MC and a match has been found by PAT
          if ((*iter).genParton()){
            HbbAnalysis::FourMomentum lGenFourVec;
            HbbAnalysis::GenVars lGen;
            lGenFourVec.px = (*iter).genParton()->px();
            lGenFourVec.py = (*iter).genParton()->py();
            lGenFourVec.pz = (*iter).genParton()->pz();
            lGenFourVec.e = (*iter).genParton()->energy();
            lGen.charge = (*iter).genParton()->charge();
            lGen.pdgId = (*iter).genParton()->pdgId();
            lGen.status = (*iter).genParton()->status();
            lGen.vertex.vx = (*iter).genParton()->vx();
            lGen.vertex.vy = (*iter).genParton()->vy();
            lGen.vertex.vz = (*iter).genParton()->vz();
            lJet.genMatches.push_back(HbbAnalysis::GenParticle(lGenFourVec, lGen));
          }
          //Construct a GenJet if processing MC and a match has been found by PAT
          if ((*iter).genJet()){
            HbbAnalysis::FourMomentum lGenJetFourVec;
            HbbAnalysis::GenJetVars lGenJet;
            lGenJetFourVec.px = (*iter).genJet()->px();
            lGenJetFourVec.py = (*iter).genJet()->py();
            lGenJetFourVec.pz = (*iter).genJet()->pz();
            lGenJetFourVec.e = (*iter).genJet()->energy();
            lGenJet.flavour = 0;  //Flavour unknown in this context (will be filled for parton jets)
            lGenJet.mass = (*iter).genJet()->mass();
            lGenJet.nConstituents = (*iter).genJet()->getGenConstituents().size();
            lJet.genJetMatches.push_back(HbbAnalysis::GenJet(lGenJetFourVec, lGenJet));
          }
          
          HbbAnalysis::JetECorVars lEnergyCorrections;
          if ((*iter).jecSetsAvailable()){
            std::vector<std::string> lJECnames = (*iter).availableJECSets();
            for (unsigned int iSet(0); iSet < lJECnames.size(); ++iSet){
              if (debug_>1) std::cout << " -- JEC SET " << iSet << " " << lJECnames[iSet] << std::endl;
              std::vector<std::string> lJEClevels = (*iter).availableJECLevels(iSet);
              for (unsigned int iLev(0); iLev<lJEClevels.size(); ++iLev){
                if (debug_>1) std::cout << " ---- JEC level " << iLev << " " << lJEClevels[iLev] 
                                      << " , factor = " << (*iter).jecFactor(iLev,pat::JetCorrFactors::NONE,iSet)
                                      << std::endl;
                lEnergyCorrections.Levels[(lJEClevels[iLev])] = (*iter).jecFactor(iLev,pat::JetCorrFactors::NONE,iSet);
              }
            }
          }
          
          HbbAnalysis::JetBtagVars lBtag;
          if (debug_>1) {
            const std::vector<std::pair<std::string, float> > & lDiscris = (*iter).getPairDiscri();
            std::cout << " -- btaggers: " << lDiscris.size() << " elements"
                      << std::endl;
            for (unsigned int iD(0); iD<lDiscris.size(); ++iD){
              std::cout << "---- " << iD << " " << lDiscris[iD].first << std::endl;
            }
          }

          lBtag.cSV = (*iter).bDiscriminator("combinedSecondaryVertexBJetTags");
          lBtag.cSVMVA = (*iter).bDiscriminator("combinedSecondaryVertexMVABJetTags");
          //lBtag.iPMVA = (*iter).bDiscriminator("impactParameterMVABJetTags");
          lBtag.bProba = (*iter).bDiscriminator("jetBProbabilityBJetTags");
          lBtag.probability = (*iter).bDiscriminator("jetProbabilityBJetTags");
          lBtag.sSVHP = (*iter).bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
          lBtag.sSVHE = (*iter).bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
          //lBtag.softElectronByPt = ((*iter).bDiscriminator("softElectronByPtBJetTags") < -100) ? -100 : (*iter).bDiscriminator("softElectronByPtBJetTags");
          //lBtag.softElectronByIP3d = ((*iter).bDiscriminator("softElectronByIP3dBJetTags") < -100) ? -100 : (*iter).bDiscriminator("softElectronByIP3dBJetTags");
          lBtag.softMuon = ((*iter).bDiscriminator("softMuonBJetTags") < -100) ? -100 : (*iter).bDiscriminator("softMuonBJetTags");
          lBtag.softMuonByPt = ((*iter).bDiscriminator("softMuonByPtBJetTags") < -100) ? -100 : (*iter).bDiscriminator("softMuonByPtBJetTags");
          lBtag.softMuonByIP3d = ((*iter).bDiscriminator("softMuonByIP3dBJetTags") < -100) ? -100 : (*iter).bDiscriminator("softMuonByIP3dBJetTags");
          lBtag.tCHE = (*iter).bDiscriminator("trackCountingHighEffBJetTags");
          lBtag.tCHP = (*iter).bDiscriminator("trackCountingHighPurBJetTags");

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

          if (debug_>1) {
            std::cout << " -- Filling tag info variables : " 
                      << (*iter).hasTagInfo("SecondaryVertex")
                      << std::endl;
          }

          std::vector<HbbAnalysis::JetSecVertexVars> jetSecVertexVars;

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


          bool isPF = (*iter).isPFJet();
          if (debug_>1) {
            std::cout << " -- Filling JPTPF jet variables " << std::endl;
          }

          HbbAnalysis::JPTPFJetVars lJPTPF;
          if (isPF){
            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 (debug_>1) {
            std::cout << " -- Filling PF jet variables " << std::endl;
          }
          if (isPF) {
            HbbAnalysis::PFJetVars lPF;
            lPF.chargedMuEnergy = (*iter).chargedMuEnergy();
            lPF.chargedMuEnergyFraction = (*iter).chargedMuEnergyFraction();
            lPF.neutralMultiplicity = (*iter).neutralMultiplicity();
            lPF.numberOfDaughters = (*iter).numberOfDaughters();
            lPF.HFHadronEnergy = (*iter).HFHadronEnergy();

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

          if (debug_>1) {
            std::cout << " -- Element added. " << std::endl;
          }

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

}//HbbJets


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


Revision:	1.2
Committed:	Wed Jun 6 13:42:42 2012 UTC (12 years, 11 months ago) by amagnan
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.1:	+193 -2 lines
Log Message:	add reco jets