HbbAnalysis/plugins/HbbTreeMaker.cc

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

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

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

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

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


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

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

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

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

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

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

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

#include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
#include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h" 

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


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")),
  genJetSrc_(pset.getParameter<edm::InputTag>("GenJets")),
  partonJetSrc_(pset.getParameter<edm::InputTag>("PartonJets")),
  trackSrc_(pset.getParameter<edm::InputTag>("Tracks")),
  dcsSrc_(pset.getParameter<edm::InputTag>("DCSInfo")),
  electronSrc_(pset.getParameter<edm::InputTag>("Electrons")),
  muonSrc_(pset.getParameter<edm::InputTag>("Muons")),
  caloTauSrc_(pset.getParameter<edm::InputTag>("CaloTaus")),
  pfTauSrc_(pset.getParameter<edm::InputTag>("PFTaus")),
  caloJetSrc_(pset.getParameter<edm::InputTag>("CaloJets")),
  jptJetSrc_(pset.getParameter<edm::InputTag>("JPTJets")),
  ak7JetSrc_(pset.getParameter<edm::InputTag>("AK7Jets")),
  pfJetSrc_(pset.getParameter<edm::InputTag>("PFJets")),
  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)
  //processVec_(pset.getParameter<std::vector<std::string> >("ProcessVec"))
{//constructor
 
  std::string JEC_PATH("CondFormats/JetMETObjects/data/");
  try {
    edm::FileInPath fipCalo(JEC_PATH+"Jec10V3_Uncertainty_AK5Calo.txt");
    edm::FileInPath fipPF(JEC_PATH+"Jec10V3_Uncertainty_AK5PF.txt");
    edm::FileInPath fipJPT(JEC_PATH+"Jec10V3_Uncertainty_AK5JPT.txt");
    jecUncCalo = new JetCorrectionUncertainty(fipCalo.fullPath());
    jecUncPF = new JetCorrectionUncertainty(fipPF.fullPath());
    jecUncJPT = new JetCorrectionUncertainty(fipJPT.fullPath());
  } catch(edm::Exception& e)  {
    std::cout << "JEC Uncertainty files unavailable! Exception : " << e.what() << ". " << std::endl;
    jecUncCalo = 0;
    jecUncPF = 0;
    jecUncJPT = 0;
  }

}//constructor

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

}//destructor


void HbbTreeMaker::beginJob(){//beginJob


  edm::Service<TFileService> lFileService;

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

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

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


}//beginJob

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

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

void HbbTreeMaker::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());
  event_->bx(aEvt.bunchCrossing());

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

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

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

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

    if (mcProd.isValid() && lRunProd.isValid()){
    HbbGenInfo(mcProd,lRunProd);
    }

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

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

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

    //save components
    if (lPartonJets.isValid()) HbbGenJets(lPartonJets,event_->partonjets(),true);

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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


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

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

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

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


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

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

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

  //std::cout << "Processing PF jets:" << std::endl;
  HbbJets(lPfJetCollection,jetFlav_,lGenParticles,event_->pfJets());
  
  edm::Handle<std::vector<pat::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_->hltObjects());
  }


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

}//analyze


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

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

    aBS.bsVars(lbs);
  }

}


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

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


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

      
          HbbAnalysis::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();
    lReco.hltMatch = (*iter).triggerObjectMatches().size();
          
          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.idAndIso = 
            (patSeleR95(*iter)<<11) | (patSeleR90(*iter)<<10) |
            (patSeleR85(*iter)<<9) | (patSeleR80(*iter)<<8) |
            (patSeleR70(*iter)<<7) | (patSeleR60(*iter)<<6) |
            (patSeleC95(*iter)<<5) | (patSeleC90(*iter)<<4) |
            (patSeleC85(*iter)<<3) | (patSeleC80(*iter)<<2) |
            (patSeleC70(*iter)<<1) | (patSeleC60(*iter));
          lID.electronIDs = (*iter).electronIDs();
          lID.hOverE = (*iter).hadronicOverEm();
          lID.deltaPhiIn = (*iter).deltaPhiSuperClusterTrackAtVtx();
          lID.deltaEtaIn = (*iter).deltaEtaSuperClusterTrackAtVtx();
          lID.ecalDrivenSeed = (*iter).ecalDrivenSeed();
          lID.trackerDrivenSeed = (*iter).trackerDrivenSeed();
    lID.dB = (*iter).dB();

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

}//HbbElectrons


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

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

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

          HbbAnalysis::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();
    lReco.hltMatch = (*iter).triggerObjectMatches().size();

          HbbAnalysis::MuTrkVars lTrk;
          reco::TrackRef glmuon = recoMuon->globalTrack();
          reco::TrackRef tkmuon = recoMuon->innerTrack();

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

          HbbAnalysis::MuIsoVars lIsoR03;
          lIsoR03.sumPt = recoMuon->isolationR03().sumPt;
          lIsoR03.emEt = recoMuon->isolationR03().emEt;
          lIsoR03.hadEt = recoMuon->isolationR03().hadEt;
          lIsoR03.nTracks = recoMuon->isolationR03().nTracks;
          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;
          lID.type =  (recoMuon->isCaloMuon()<<3) | (recoMuon->isStandAloneMuon()<<2) | (recoMuon->isTrackerMuon()<<1) | recoMuon->isGlobalMuon();

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

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

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

}//HbbMuons

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

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


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

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

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

          assert (isCalo == !isPF);

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

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

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

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

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

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

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

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

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

          }//pf


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

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


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

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


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

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

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

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


          }


          if (isCalo){

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

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

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

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

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

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

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

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

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

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

}//HbbTaus

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

          HbbAnalysis::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.hltMatch = (*iter).triggerObjectMatches().size();
          //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"));


//        std::cout << " -- Current JEC SET " << (*iter).currentJECSet() 
//                  << " Level " << (*iter).currentJECLevel() 
//                  << std::endl;
               
          HbbAnalysis::JetECorVars lEnergyCorrections;
          if ((*iter).jecSetsAvailable()){
            float lJecFactor = (*iter).jecFactor("Uncorrected","NONE",(*iter).currentJECSet());
            if (lJecFactor > 0) lCommon.rawpT = (*iter).pt()*lJecFactor;
            if (debug_){ 
              std::cout << " ---------------------------------------------- " << std::endl
                        << " - Raw to Cor Jet # " << iEle << " = " << std::endl
                        << " -- Raw pT = " << lCommon.rawpT << " eta " << lReco.eta << " phi " << lReco.phi << std::endl;
            }

            std::vector<std::string> lJECnames = (*iter).availableJECSets();
      for (unsigned int iSet(0); iSet < lJECnames.size(); ++iSet){
        if (debug_) std::cout << " -- JEC SET " << iSet << " " << lJECnames[iSet] << std::endl;
        std::vector<std::string> lJEClevels = (*iter).availableJECLevels(iSet);
        for (unsigned int iLev(0); iLev<lJEClevels.size(); ++iLev){
          if (debug_) std::cout << " ---- JEC level " << iLev << " " << lJEClevels[iLev] 
            << " , factor = " << (*iter).jecFactor(iLev,pat::JetCorrFactors::NONE,iSet)
              << std::endl;
          lEnergyCorrections.Levels.push_back(std::pair<std::string,double>(lJEClevels[iLev],(*iter).jecFactor(iLev,pat::JetCorrFactors::NONE,iSet)));
        }
      }
      if (debug_) std::cout << " -- Reco pT = " << lReco.pT
        << std::endl;
    }

    if(jecUncCalo && jecUncPF && jecUncJPT){
      if ((*iter).isCaloJet()) {
        jecUncCalo->setJetEta((*iter).eta());
        jecUncCalo->setJetPt((*iter).pt());
        lEnergyCorrections.Levels.push_back(std::pair<std::string,double>("Uncertainty",jecUncCalo->getUncertainty(true)));
      } else if ((*iter).isPFJet()) {
        jecUncPF->setJetEta((*iter).eta());
        jecUncPF->setJetPt((*iter).pt());
        lEnergyCorrections.Levels.push_back(std::pair<std::string,double>("Uncertainty",jecUncPF->getUncertainty(true)));
      } else if ((*iter).isJPTJet()) {
        jecUncJPT->setJetEta((*iter).eta());
        jecUncJPT->setJetPt((*iter).pt());
        lEnergyCorrections.Levels.push_back(std::pair<std::string,double>("Uncertainty",jecUncJPT->getUncertainty(true)));
      }
    }

          //lEnergyCorrections.L1 = (*iter).jecFactor("L1Offset","NONE",(*iter).currentJECSet());
          //lEnergyCorrections.L1JPT = (*iter).jecFactor("L1JPTOffset","NONE",(*iter).currentJECSet());
          //lEnergyCorrections.L2 = (*iter).jecFactor("L2Relative","NONE",(*iter).currentJECSet());
          //lEnergyCorrections.L3 = (*iter).jecFactor("L3Absolute","NONE",(*iter).currentJECSet());
          //if (processData_) lEnergyCorrections.L2L3Residual = (*iter).jecFactor("L2L3Residual","NONE",(*iter).currentJECSet());
          
          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.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") < -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");

          //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,lEnergyCorrections);
            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;
      reco::TrackRefVector pionsInVInC = (*iter).pionsInVertexInCalo();
      reco::TrackRefVector pionsInVOuC = (*iter).pionsInVertexOutCalo();
      reco::TrackRefVector pionsOuVInC = (*iter).pionsOutVertexInCalo();
      
      reco::TrackRefVector muonsInVInC = (*iter).muonsInVertexInCalo();
      reco::TrackRefVector muonsInVOuC = (*iter).muonsInVertexOutCalo();
      reco::TrackRefVector muonsOuVInC = (*iter).muonsOutVertexInCalo();

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

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

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

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

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

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

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

            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,lEnergyCorrections);
            aVec.push_back(lObj);
          }
          if (isPF) {
            HbbAnalysis::PFJetVars lPF;
            lPF.chargedMuEnergy = (*iter).chargedMuEnergy();
            lPF.chargedMuEnergyFraction = (*iter).chargedMuEnergyFraction();
            lPF.neutralMultiplicity = (*iter).neutralMultiplicity();
      lPF.numberOfDaughters = (*iter).numberOfDaughters();
      lPF.HFHadronEnergy = (*iter).HFHadronEnergy();

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

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

}//HbbJets

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

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

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

void HbbTreeMaker::HbbTrigger(const edm::Handle<edm::TriggerResults> & aCol,
                              const edm::TriggerNames & aNames,
                              std::vector<HbbAnalysis::Trigger> & aVec)
{//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();

      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_){
        HbbAnalysis::GenParticle lObj(lMC);
        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
        aVec.push_back(lObj);
      }
    }//loop on particles


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


}//genparticles

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

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

          HbbAnalysis::GenVars lGenJet;
          if (!processData_ ){
            lGenJet.valid = true;
            lGenJet.E = (*iter).energy();
            lGenJet.pT = (*iter).pt();
            lGenJet.eta = (*iter).eta();
            lGenJet.y = (*iter).y();
            lGenJet.phi = (*iter).phi();
            lGenJet.charge = (*iter).charge();

            //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;
            lGenJet.pdgId = (*iter).pdgId();
            lGenJet.status = (*iter).status();
            lGenJet.mass = (*iter).mass();
            lGenJet.vx = (*iter).vx();
            lGenJet.vy = (*iter).vy();
            lGenJet.vz = (*iter).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::GenJet lObj(lGenJet);
          if (aSaveConstituants){
            for (unsigned int iC(0); iC<(*iter).getGenConstituents().size(); ++iC){
              const reco::GenParticle* lPart = (*iter).getGenConstituent(iC);

              HbbAnalysis::MCVars lMC;
              lMC.pT = lPart->pt();
              lMC.eta = lPart->eta();
              lMC.phi = lPart->phi();
              lMC.pdgId = lPart->pdgId();
              lObj.addConstituent(lMC);
            }
          }
          else lObj.nConstituents((*iter).getGenConstituents().size());
          aVec.push_back(lObj);
        }
    }
  }
}

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

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

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

        HbbAnalysis::VertexVars lVtx;

        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));
      lVtx.trackpT.push_back((*(*lTrk)).pt());
          }
        }

        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;
        }
        if (lVtx.trackpT.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 trackpT = " << lVtx.trackpT.size()
                   << std::endl;
        }

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

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

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


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

}

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

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

  const trigger::Keys & k = aCol->filterKeys(index);
         
  for (trigger::Keys::const_iterator ki = k.begin(); ki !=k.end(); ++ki ) {
    HbbAnalysis::HLTVars lVars;
    lVars.pT = toc[*ki].pt();
    lVars.eta = toc[*ki].eta();
    lVars.phi = toc[*ki].phi();
    lVars.id = toc[*ki].id();
    lVars.mass = toc[*ki].mass();
    lVars.type = aType;
    HbbAnalysis::HLTObject lObj(lVars);
    aVec.push_back(lObj);           
  }
}

void HbbTreeMaker::HbbGenInfo(const edm::Handle<edm::HepMCProduct> & amcProd,
                              const edm::Handle<GenRunInfoProduct> & aRunProd) 
{
  const 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);

}

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

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

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


  //add genParticle inside vector
  for (unsigned int imc = 0;imc<lpup.size();++imc)
    {//loop on particles
      
      HbbAnalysis::MCVars lMC;
      lMC.index = imc;
      lMC.E = (lpup[imc])[3];

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

      lMC.pT = lVec.Pt();
      // 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(lMC);
      lObj.setParent(lmothup[imc].first);
      lObj.setParent(lmothup[imc].second);
      aVec.push_back(lObj);

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

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


Revision:	1.42
Committed:	Thu Sep 1 14:07:31 2011 UTC (13 years, 8 months ago) by amagnan
Content type:	text/plain
Branch:	MAIN
Changes since 1.41:	+23 -40 lines
Log Message:	minor changes related to genParticles and gen info