#include <PFJetBenchmark.h>

Public Member Functions
void	gettrue (const reco::GenJet *truth, double &true_ChargedHadEnergy, double &true_NeutralHadEnergy, double &true_NeutralEmEnergy)
	PFJetBenchmark ()
void	printGenJet (const reco::GenJet *)
void	printPFJet (const reco::PFJet *)
void	process (const reco::PFJetCollection &, const reco::GenJetCollection &)
double	resChargedHadEnergyMax () const
double	resNeutralEmEnergyMax () const
double	resNeutralHadEnergyMax () const
double	resPtMax () const
void	setup (std::string Filename, bool debug, bool plotAgainstReco=0, bool onlyTwoJets=1, double deltaRMax=0.1, std::string benchmarkLabel_="ParticleFlow", double recPt=-1, double maxEta=-1, DQMStore *dbe_store=NULL)
void	write ()
virtual	~PFJetBenchmark ()
Protected Attributes
PFBenchmarkAlgo *	algo_
DQMStore *	dbe_
bool	debug_
double	deltaRMax_
unsigned int	entry_
double	maxEta_cut
bool	onlyTwoJets_
bool	plotAgainstReco_
double	recPt_cut
double	resChargedHadEnergyMax_
double	resNeutralEmEnergyMax_
double	resNeutralHadEnergyMax_
double	resPtMax_
Private Attributes
TFile *	file_
TH2F *	hBDEtavsPt
TH2F *	hBDPhivsPt
TH1F *	hBNCH
TH2F *	hBNCH0vsPt
TH2F *	hBNCH1vsPt
TH2F *	hBNCH2vsPt
TH2F *	hBNCH3vsPt
TH2F *	hBNCH4vsPt
TH2F *	hBNCH5vsPt
TH2F *	hBNCH6vsPt
TH2F *	hBNCHvsPt
TH1F *	hBRCHE
TH2F *	hBRCHEvsPt
TH2F *	hBRHCALvsP
TH2F *	hBRHONLvsP
TH1F *	hBRNEE
TH2F *	hBRNEEvsPt
TH1F *	hBRneut
TH2F *	hBRNEUTvsP
TH2F *	hBRneutvsPt
TH1F *	hBRNHE
TH2F *	hBRNHEvsPt
TH2F *	hBRNONLvsP
TH1F *	hBRPt
TH1F *	hBRPt100_150
TH1F *	hBRPt1250_2000
TH1F *	hBRPt150_200
TH1F *	hBRPt2000_5000
TH1F *	hBRPt200_250
TH1F *	hBRPt20_40
TH1F *	hBRPt250_300
TH1F *	hBRPt300_400
TH1F *	hBRPt400_500
TH1F *	hBRPt40_60
TH1F *	hBRPt500_750
TH1F *	hBRPt60_80
TH1F *	hBRPt750_1250
TH1F *	hBRPt80_100
TH2F *	hBRPtvsPt
TH2F *	hDEtavsEta
TH2F *	hDPhivsEta
TH2F *	hEDEtavsPt
TH2F *	hEDPhivsPt
TH1F *	hENCH
TH2F *	hENCH0vsPt
TH2F *	hENCH1vsPt
TH2F *	hENCH2vsPt
TH2F *	hENCH3vsPt
TH2F *	hENCH4vsPt
TH2F *	hENCH5vsPt
TH2F *	hENCH6vsPt
TH2F *	hENCHvsPt
TH1F *	hERCHE
TH2F *	hERCHEvsPt
TH2F *	hERHCALvsP
TH2F *	hERHONLvsP
TH1F *	hERNEE
TH2F *	hERNEEvsPt
TH1F *	hERneut
TH2F *	hERNEUTvsP
TH2F *	hERneutvsPt
TH1F *	hERNHE
TH2F *	hERNHEvsPt
TH2F *	hERNONLvsP
TH1F *	hERPt
TH1F *	hERPt100_150
TH1F *	hERPt1250_2000
TH1F *	hERPt150_200
TH1F *	hERPt2000_5000
TH1F *	hERPt200_250
TH1F *	hERPt20_40
TH1F *	hERPt250_300
TH1F *	hERPt300_400
TH1F *	hERPt400_500
TH1F *	hERPt40_60
TH1F *	hERPt500_750
TH1F *	hERPt60_80
TH1F *	hERPt750_1250
TH1F *	hERPt80_100
TH2F *	hERPtvsPt
TH2F *	hFDEtavsPt
TH2F *	hFDPhivsPt
TH1F *	hFNCH
TH2F *	hFNCH0vsPt
TH2F *	hFNCH1vsPt
TH2F *	hFNCH2vsPt
TH2F *	hFNCH3vsPt
TH2F *	hFNCH4vsPt
TH2F *	hFNCH5vsPt
TH2F *	hFNCH6vsPt
TH2F *	hFNCHvsPt
TH1F *	hFRCHE
TH2F *	hFRCHEvsPt
TH2F *	hFRHCALvsP
TH2F *	hFRHONLvsP
TH1F *	hFRNEE
TH2F *	hFRNEEvsPt
TH1F *	hFRneut
TH2F *	hFRNEUTvsP
TH2F *	hFRneutvsPt
TH1F *	hFRNHE
TH2F *	hFRNHEvsPt
TH2F *	hFRNONLvsP
TH1F *	hFRPt
TH1F *	hFRPt100_150
TH1F *	hFRPt1250_2000
TH1F *	hFRPt150_200
TH1F *	hFRPt2000_5000
TH1F *	hFRPt200_250
TH1F *	hFRPt20_40
TH1F *	hFRPt250_300
TH1F *	hFRPt300_400
TH1F *	hFRPt400_500
TH1F *	hFRPt40_60
TH1F *	hFRPt500_750
TH1F *	hFRPt60_80
TH1F *	hFRPt750_1250
TH1F *	hFRPt80_100
TH2F *	hFRPtvsPt
TH1F *	hjetsEta
TH1F *	hjetsPt
TH2F *	hNCH0vsEta
TH2F *	hNCH1vsEta
TH2F *	hNCH2vsEta
TH2F *	hNCH3vsEta
TH2F *	hNCH4vsEta
TH2F *	hNCH5vsEta
TH2F *	hNCH6vsEta
TH2F *	hNCHvsEta
TH1F *	hNjets
TH2F *	hRCHEvsEta
TH2F *	hRHCALvsEta
TH2F *	hRHONLvsEta
TH2F *	hRNEUTvsEta
TH2F *	hRNeutvsEta
TH2F *	hRNONLvsEta
TH2F *	hRPtvsEta
std::string	outputFile_

Detailed Description

Definition at line 38 of file PFJetBenchmark.h.

Constructor & Destructor Documentation

PFJetBenchmark::PFJetBenchmark ( )

Definition at line 40 of file PFJetBenchmark.cc.

: file_(0), entry_(0) {}

PFJetBenchmark::~PFJetBenchmark ( ) [virtual]

Definition at line 42 of file PFJetBenchmark.cc.

References file_.

                                {
  if(file_) file_->Close();
}

Member Function Documentation

void PFJetBenchmark::gettrue	(	const reco::GenJet *	truth,
		double &	true_ChargedHadEnergy,
		double &	true_NeutralHadEnergy,
		double &	true_NeutralEmEnergy
	)

Definition at line 663 of file PFJetBenchmark.cc.

References abs, reco::LeafCandidate::energy(), reco::GenJet::getGenConstituents(), i, and reco::LeafCandidate::pdgId().

Referenced by process().

                                                                                          {
  std::vector <const GenParticle*> mcparts = truth->getGenConstituents ();
  true_NeutralEmEnergy = 0.;
  true_ChargedHadEnergy = 0.;
  true_NeutralHadEnergy = 0.;
  // for each MC particle in turn  
  for (unsigned i = 0; i < mcparts.size (); i++) {
    const GenParticle* mcpart = mcparts[i];
    int PDG = std::abs( mcpart->pdgId());
    double e = mcpart->energy(); 
    switch(PDG){  // start PDG switch
    case 22: // photon
      true_NeutralEmEnergy += e;
      break;
    case 211: // pi
    case 321: // K
    case 2212: // p
    case 11: //electrons (until recognised)
      true_ChargedHadEnergy += e;
      break;
    case 310: // K_S0
    case 130: // K_L0
    case 3122: // Lambda0
    case 2112: // n0
      true_NeutralHadEnergy += e;
    default:
      break;
    }  // end PDG switch                
  }  // end loop on constituents.
}

void PFJetBenchmark::printGenJet ( const reco::GenJet * truth )

Definition at line 716 of file PFJetBenchmark.cc.

References gather_cfg::cout, reco::LeafCandidate::eta(), reco::GenJet::getGenConstituents(), i, reco::LeafCandidate::p(), reco::LeafCandidate::pdgId(), reco::LeafCandidate::phi(), reco::LeafCandidate::pt(), reco::LeafCandidate::px(), reco::LeafCandidate::py(), and reco::LeafCandidate::pz().

Referenced by process().

                                                        {
  std::vector <const GenParticle*> mcparts = truth->getGenConstituents ();
  cout << "GenJet p/px/py/pz/pt: " << truth->p() << '/' << truth->px () 
       << '/' << truth->py() << '/' << truth->pz() << '/' << truth->pt() << endl
       << "    eta/phi: " << truth->eta () << '/' << truth->phi () << endl
       << "    # of constituents: " << mcparts.size() << endl;
  cout << "    constituents:" << endl;
  for (unsigned i = 0; i < mcparts.size (); i++) {
    const GenParticle* mcpart = mcparts[i];
    cout << "      #" << i << "  PDG code:" << mcpart->pdgId() 
         << ", p/pt/eta/phi: " << mcpart->p() << '/' << mcpart->pt() 
         << '/' << mcpart->eta() << '/' << mcpart->phi() << endl;       
  }    
}

void PFJetBenchmark::printPFJet ( const reco::PFJet * pfj )

Definition at line 695 of file PFJetBenchmark.cc.

References reco::PFJet::chargedEmEnergy(), reco::PFJet::chargedHadronEnergy(), reco::PFJet::chargedMuEnergy(), reco::PFJet::chargedMultiplicity(), gather_cfg::cout, reco::LeafCandidate::eta(), reco::PFJet::neutralEmEnergy(), reco::PFJet::neutralHadronEnergy(), reco::PFJet::neutralMultiplicity(), reco::LeafCandidate::p(), reco::LeafCandidate::phi(), reco::PFJet::print(), reco::LeafCandidate::pt(), reco::LeafCandidate::px(), reco::LeafCandidate::py(), and reco::LeafCandidate::pz().

Referenced by process().

                                                   {
  cout<<setiosflags(ios::right);
  cout<<setiosflags(ios::fixed);
  cout<<setprecision(3);

  cout << "PFJet  p/px/py/pz/pt: " << pfj->p() << "/" << pfj->px () 
       << "/" << pfj->py() << "/" << pfj->pz() << "/" << pfj->pt() << endl
       << "    eta/phi: " << pfj->eta () << "/" << pfj->phi () << endl                  
       << "    PFJet specific:" << std::endl
       << "      charged/neutral hadrons energy: " << pfj->chargedHadronEnergy () << '/' << pfj->neutralHadronEnergy () << endl
       << "      charged/neutral em energy: " << pfj->chargedEmEnergy () << '/' << pfj->neutralEmEnergy () << endl
       << "      charged muon energy: " << pfj->chargedMuEnergy () << '/' << endl
       << "      charged/neutral multiplicity: " << pfj->chargedMultiplicity () << '/' << pfj->neutralMultiplicity () << endl;
  
  // And print the constituents
  std::cout << pfj->print() << std::endl;

  cout<<resetiosflags(ios::right|ios::fixed);
}

void PFJetBenchmark::process	(	const reco::PFJetCollection &	pfJets,
		const reco::GenJetCollection &	genJets
	)

Definition at line 228 of file PFJetBenchmark.cc.

References abs, algo_, sipixelsummary::Barrel, reco::PFJet::chargedHadronEnergy(), reco::PFJet::chargedMultiplicity(), gather_cfg::cout, debug_, deltaR(), PFBenchmarkAlgo::deltaR(), deltaRMax_, entry_, reco::LeafCandidate::eta(), reco::PFJet::getPFConstituents(), gettrue(), hBDEtavsPt, hBDPhivsPt, hBNCH, hBNCH0vsPt, hBNCH1vsPt, hBNCH2vsPt, hBNCH3vsPt, hBNCH4vsPt, hBNCH5vsPt, hBNCH6vsPt, hBNCHvsPt, hBRCHE, hBRCHEvsPt, hBRHCALvsP, hBRHONLvsP, hBRNEE, hBRNEEvsPt, hBRneut, hBRNEUTvsP, hBRneutvsPt, hBRNHE, hBRNHEvsPt, hBRNONLvsP, hBRPt, hBRPt100_150, hBRPt1250_2000, hBRPt150_200, hBRPt2000_5000, hBRPt200_250, hBRPt20_40, hBRPt250_300, hBRPt300_400, hBRPt400_500, hBRPt40_60, hBRPt500_750, hBRPt60_80, hBRPt750_1250, hBRPt80_100, hBRPtvsPt, hDEtavsEta, hDPhivsEta, hEDEtavsPt, hEDPhivsPt, hENCH, hENCH0vsPt, hENCH1vsPt, hENCH2vsPt, hENCH3vsPt, hENCH4vsPt, hENCH5vsPt, hENCH6vsPt, hENCHvsPt, hERCHE, hERCHEvsPt, hERHCALvsP, hERHONLvsP, hERNEE, hERNEEvsPt, hERneut, hERNEUTvsP, hERneutvsPt, hERNHE, hERNHEvsPt, hERNONLvsP, hERPt, hERPt100_150, hERPt1250_2000, hERPt150_200, hERPt2000_5000, hERPt200_250, hERPt20_40, hERPt250_300, hERPt300_400, hERPt400_500, hERPt40_60, hERPt500_750, hERPt60_80, hERPt750_1250, hERPt80_100, hERPtvsPt, hFDEtavsPt, hFDPhivsPt, hFRCHE, hFRCHEvsPt, hFRHCALvsP, hFRHONLvsP, hFRNEE, hFRNEEvsPt, hFRneut, hFRNEUTvsP, hFRneutvsPt, hFRNHE, hFRNHEvsPt, hFRNONLvsP, hFRPt, hFRPt100_150, hFRPt1250_2000, hFRPt150_200, hFRPt2000_5000, hFRPt200_250, hFRPt20_40, hFRPt250_300, hFRPt300_400, hFRPt400_500, hFRPt40_60, hFRPt500_750, hFRPt60_80, hFRPt750_1250, hFRPt80_100, hFRPtvsPt, hjetsEta, hjetsPt, hNCH0vsEta, hNCH1vsEta, hNCH2vsEta, hNCH3vsEta, hNCH4vsEta, hNCH5vsEta, hNCH6vsEta, hNCHvsEta, hNjets, hRCHEvsEta, hRHCALvsEta, hRHONLvsEta, hRNEUTvsEta, hRNeutvsEta, hRNONLvsEta, hRPtvsEta, i, edm::Ref< C, T, F >::isNull(), j, PFBenchmarkAlgo::matchByDeltaR(), maxEta_cut, reco::PFJet::neutralEmEnergy(), reco::PFJet::neutralHadronEnergy(), onlyTwoJets_, reco::LeafCandidate::p(), reco::LeafCandidate::phi(), plotAgainstReco_, printGenJet(), printPFJet(), reco::LeafCandidate::pt(), recPt_cut, resChargedHadEnergyMax_, resNeutralEmEnergyMax_, resNeutralHadEnergyMax_, and resPtMax_.

Referenced by PFJetBenchmarkAnalyzer::analyze(), and PFRootEventManager::processEntry().

                                                                                                   {
  // loop over reco  pf  jets
  resPtMax_ = 0.;
  resChargedHadEnergyMax_ = 0.;
  resNeutralHadEnergyMax_ = 0.;
  resNeutralEmEnergyMax_ = 0.; 
  int NPFJets = 0;
        
  for(unsigned i=0; i<pfJets.size(); i++) {   

    // Count the number of jets with a larger energy
    unsigned highJets = 0;
    for(unsigned j=0; j<pfJets.size(); j++) { 
      if ( j != i && pfJets[j].pt() > pfJets[i].pt() ) highJets++;
    }
    if ( onlyTwoJets_ && highJets > 1 ) continue;
                
                
    const reco::PFJet& pfj = pfJets[i];
    double rec_pt = pfj.pt();
    double rec_eta = pfj.eta();
    double rec_phi = pfj.phi();

    // skip PFjets with pt < recPt_cut GeV
    if (rec_pt<recPt_cut and recPt_cut != -1.) continue;
    // skip PFjets with eta > maxEta_cut
    if (fabs(rec_eta)>maxEta_cut and maxEta_cut != -1.) continue;

    NPFJets++;
                
    // fill inclusive PFjet distribution pt > 20 GeV
    hNjets->Fill(NPFJets);
    hjetsPt->Fill(rec_pt);
    hjetsEta->Fill(rec_eta);

    // separate Barrel PFJets from Endcap PFJets
    bool Barrel = false;
    bool Endcap = false;
    bool Forward = false;
    if (std::abs(rec_eta) < 1.4 ) Barrel = true;
    if (std::abs (rec_eta) > 1.6 && std::abs (rec_eta) < 2.4 ) Endcap = true;
    if (std::abs (rec_eta) > 2.5 && std::abs (rec_eta) < 2.9 ) Forward = true;
    if (std::abs (rec_eta) > 3.1 && std::abs (rec_eta) < 4.7 ) Forward = true;

    // do only barrel for now
    //  if(!Barrel) continue;

    // look for the closets gen Jet : truth
    const GenJet *truth = algo_->matchByDeltaR(&pfj,&genJets);
    if(!truth) continue;   
    double deltaR = algo_->deltaR(&pfj, truth);
    // check deltaR is small enough
    if(deltaR < deltaRMax_ || (abs(rec_eta)>2.5 && deltaR < 0.2) || deltaRMax_ == -1.0 ) {//start case deltaR < deltaRMax

      // generate histograms comparing the reco and truth candidate (truth = closest in delta-R) 
      // get the quantities to place on the denominator and/or divide by
      double pt_denom;
      double true_E = truth->p();
      double true_pt = truth->pt();
      double true_eta = truth->eta();
      double true_phi = truth->phi();

      if (plotAgainstReco_) {pt_denom = rec_pt;}
      else {pt_denom = true_pt;}
      // get true specific quantities
      double true_ChargedHadEnergy;
      double true_NeutralHadEnergy;
      double true_NeutralEmEnergy;
      gettrue (truth, true_ChargedHadEnergy, true_NeutralHadEnergy, true_NeutralEmEnergy);
      double true_NeutralEnergy = true_NeutralHadEnergy + true_NeutralEmEnergy;
      double rec_ChargedHadEnergy = pfj.chargedHadronEnergy();
      double rec_NeutralHadEnergy = pfj.neutralHadronEnergy();
      double rec_NeutralEmEnergy = pfj.neutralEmEnergy();
      double rec_NeutralEnergy = rec_NeutralHadEnergy + rec_NeutralEmEnergy;
      double rec_ChargedMultiplicity = pfj.chargedMultiplicity();
      std::vector <PFCandidatePtr> constituents = pfj.getPFConstituents ();
      std::vector <unsigned int> chMult(7, static_cast<unsigned int>(0)); 
      for (unsigned ic = 0; ic < constituents.size (); ++ic) {
        if ( constituents[ic]->particleId() > 3 ) continue;
        reco::TrackRef trackRef = constituents[ic]->trackRef();
        if ( trackRef.isNull() ) {
          //std::cout << "Warning in entry " << entry_ 
          //        << " : a track with Id " << constituents[ic]->particleId() 
          //        << " has no track ref.." << std::endl;
          continue;
        }
        unsigned int iter = 0; 
        switch (trackRef->algo()) {
        case TrackBase::ctf:
        case TrackBase::iter0:
          iter = 0;
          break;
        case TrackBase::iter1:
          iter = 1;
          break;
        case TrackBase::iter2:
          iter = 2;
          break;
        case TrackBase::iter3:
          iter = 3;
          break;
        case TrackBase::iter4:
          iter = 4;
          break;
        case TrackBase::iter5:
          iter = 5;
          break;
        default:
          iter = 6;
          std::cout << "Warning in entry " << entry_ << " : iter = 6... " << std::endl;
          break;
        }
        ++(chMult[iter]);
      }

      bool plot1 = false;
      bool plot2 = false;
      bool plot3 = false;
      bool plot4 = false;
      bool plot5 = false;
      bool plot6 = false;
      bool plot7 = false;
      double cut1 = 0.0001;
      double cut2 = 0.0001;
      double cut3 = 0.0001;
      double cut4 = 0.0001;
      double cut5 = 0.0001;
      double cut6 = 0.0001;
      double cut7 = 0.0001;
      double resPt =0.;
      double resChargedHadEnergy= 0.;
      double resNeutralHadEnergy= 0.;
      double resNeutralEmEnergy= 0.;
      double resNeutralEnergy= 0.;

      double resHCALEnergy = 0.;
      double resNEUTEnergy = 0.;
      if ( rec_NeutralHadEnergy > cut6 && rec_ChargedHadEnergy < cut1 ) { 
        double true_NEUTEnergy = true_NeutralHadEnergy + true_NeutralEmEnergy;
        double true_HCALEnergy = true_NEUTEnergy - rec_NeutralEmEnergy;
        double rec_NEUTEnergy = rec_NeutralHadEnergy+rec_NeutralEmEnergy; 
        double rec_HCALEnergy = rec_NeutralHadEnergy; 
        resHCALEnergy = (rec_HCALEnergy-true_HCALEnergy)/rec_HCALEnergy;
        resNEUTEnergy = (rec_NEUTEnergy-true_NEUTEnergy)/rec_NEUTEnergy;
        if ( rec_NeutralEmEnergy > cut7 ) {
          plot6 = true;
        } else {
          plot7 = true;
        }
      }

      // get relative delta quantities (protect against division by zero!)
      if (true_pt > 0.0001){
        resPt = (rec_pt -true_pt)/true_pt ; 
        plot1 = true;}
      if (true_ChargedHadEnergy > cut1){
        resChargedHadEnergy = (rec_ChargedHadEnergy- true_ChargedHadEnergy)/true_ChargedHadEnergy;
        plot2 = true;}
      if (true_NeutralHadEnergy > cut2){
        resNeutralHadEnergy = (rec_NeutralHadEnergy- true_NeutralHadEnergy)/true_NeutralHadEnergy;
        plot3 = true;}
      else 
        if (rec_NeutralHadEnergy > cut3){
          resNeutralHadEnergy = (rec_NeutralHadEnergy- true_NeutralHadEnergy)/rec_NeutralHadEnergy;
          plot3 = true;}
      if (true_NeutralEmEnergy > cut4){
        resNeutralEmEnergy = (rec_NeutralEmEnergy- true_NeutralEmEnergy)/true_NeutralEmEnergy;
        plot4 = true;}
      if (true_NeutralEnergy > cut5){
        resNeutralEnergy = (rec_NeutralEnergy- true_NeutralEnergy)/true_NeutralEnergy;
        plot5 = true;}
      
      //double deltaEta = algo_->deltaEta(&pfj, truth);
      //double deltaPhi = algo_->deltaPhi(&pfj, truth);

      // Print outliers for further debugging
      if ( ( resPt > 0.2 && true_pt > 100. ) || 
           ( resPt < -0.5 && true_pt > 100. ) ) {
        //if ( ( true_pt > 50. && 
        //     ( ( truth->eta()>3.0 && rec_eta-truth->eta() < -0.1 ) || 
        //       ( truth->eta()<-3.0 && rec_eta-truth->eta() > 0.1 ) ))) {
        std::cout << "Entry " << entry_ 
                  << " resPt = " << resPt
                  <<" resCharged  " << resChargedHadEnergy
                  <<" resNeutralHad  " << resNeutralHadEnergy
                  << " resNeutralEm  " << resNeutralEmEnergy
                  << " pT (T/R) " << true_pt << "/" << rec_pt 
                  << " Eta (T/R) " << truth->eta() << "/" << rec_eta 
                  << " Phi (T/R) " << truth->phi() << "/" << rec_phi 
                  << std::endl;

        // check overlapping PF jets
        const reco::PFJet* pfoj = 0; 
        double dRo = 1E9;
        for(unsigned j=0; j<pfJets.size(); j++) { 
          const reco::PFJet& pfo = pfJets[j];
          if ( j != i &&  algo_->deltaR(&pfj,&pfo) < dRo && pfo.pt() > 0.25*pfj.pt()) { 
            dRo = algo_->deltaR(&pfj,&pfo);     
            pfoj = &pfo;
          }
        }
        
        // Check overlapping Gen Jet 
        math::XYZTLorentzVector overlappinGenJet(0.,0.,0.,0.);
        const reco::GenJet* genoj = 0;
        double dRgo = 1E9;
        for(unsigned j=0; j<genJets.size(); j++) { 
          const reco::GenJet* gjo = &(genJets[j]);
          if ( gjo != truth && algo_->deltaR(truth,gjo) < dRgo && gjo->pt() > 0.25*truth->pt() ) { 
            dRgo = algo_->deltaR(truth,gjo);
            genoj = gjo;
          }
        }
        
        if ( dRo < 0.8 && dRgo < 0.8 && algo_->deltaR(genoj,pfoj) < 2.*deltaRMax_ ) 
          std::cout << "Excess probably due to overlapping jets (DR = " <<   algo_->deltaR(genoj,pfoj) << "),"
                    << " at DeltaR(T/R) = " << dRgo << "/" << dRo  
                    << " with pT(T/R) " << genoj->pt() << "/" << pfoj->pt()
                    << " and Eta (T/R) " << genoj->eta() << "/" << pfoj->eta()
                    << " and Phi (T/R) " << genoj->phi() << "/" << pfoj->phi()
                    << std::endl;
      }

      if(std::abs(resPt) > std::abs(resPtMax_)) resPtMax_ = resPt;
      if(std::abs(resChargedHadEnergy) > std::abs(resChargedHadEnergyMax_) ) resChargedHadEnergyMax_ = resChargedHadEnergy;
      if(std::abs(resNeutralHadEnergy) > std::abs(resNeutralHadEnergyMax_) ) resNeutralHadEnergyMax_ = resNeutralHadEnergy;
      if(std::abs(resNeutralEmEnergy) > std::abs(resNeutralEmEnergyMax_) ) resNeutralEmEnergyMax_ = resNeutralEmEnergy;
      if (debug_) {
        cout << i <<"  =========PFJet Pt "<< rec_pt
             << " eta " << rec_eta
             << " phi " << rec_phi
             << " Charged Had Energy " << rec_ChargedHadEnergy
             << " Neutral Had Energy " << rec_NeutralHadEnergy
             << " Neutral elm Energy " << rec_NeutralEmEnergy << endl;
        cout << " matching Gen Jet Pt " << true_pt
             << " eta " << truth->eta()
             << " phi " << truth->phi()
             << " Charged Had Energy " << true_ChargedHadEnergy
             << " Neutral Had Energy " << true_NeutralHadEnergy
             << " Neutral elm Energy " << true_NeutralEmEnergy << endl;
        printPFJet(&pfj);
        //      cout<<pfj.print()<<endl;
        printGenJet(truth);
        //cout <<truth->print()<<endl;
                                
        cout << "==============deltaR " << deltaR << "  resPt " << resPt
             << " resChargedHadEnergy " << resChargedHadEnergy
             << " resNeutralHadEnergy " << resNeutralHadEnergy
             << " resNeutralEmEnergy " << resNeutralEmEnergy
             << endl;
      }
                        

      if(plot1) {
        if ( rec_eta > 0. ) 
          hDEtavsEta->Fill(true_eta,rec_eta-true_eta);
        else
          hDEtavsEta->Fill(true_eta,-rec_eta+true_eta);
        hDPhivsEta->Fill(true_eta,rec_phi-true_phi);

        hRPtvsEta->Fill(true_eta, resPt);
        hNCHvsEta->Fill(true_eta, rec_ChargedMultiplicity);
        hNCH0vsEta->Fill(true_eta,chMult[0]);
        hNCH1vsEta->Fill(true_eta,chMult[1]);
        hNCH2vsEta->Fill(true_eta,chMult[2]);
        hNCH3vsEta->Fill(true_eta,chMult[3]);
        hNCH4vsEta->Fill(true_eta,chMult[4]);
        hNCH5vsEta->Fill(true_eta,chMult[5]);
        hNCH6vsEta->Fill(true_eta,chMult[6]);
      }
      if(plot2)hRCHEvsEta->Fill(true_eta, resChargedHadEnergy);
      if(plot5)hRNeutvsEta->Fill(true_eta, resNeutralEnergy);
      if(plot6) { 
        hRHCALvsEta->Fill(true_eta, resHCALEnergy);
        hRNEUTvsEta->Fill(true_eta, resNEUTEnergy);
      }
      if(plot7) {  
        hRHONLvsEta->Fill(true_eta, resHCALEnergy);
        hRNONLvsEta->Fill(true_eta, resNEUTEnergy);
      }

      // fill histograms for relative delta quantitites of matched jets
      // delta Pt or E quantities for Barrel
      if (Barrel){
        if(plot1) { 
          hBRPt->Fill (resPt);
          if ( pt_denom >  20. && pt_denom <  40. ) hBRPt20_40->Fill (resPt);
          if ( pt_denom >  40. && pt_denom <  60. ) hBRPt40_60->Fill (resPt);
          if ( pt_denom >  60. && pt_denom <  80. ) hBRPt60_80->Fill (resPt);
          if ( pt_denom >  80. && pt_denom < 100. ) hBRPt80_100->Fill (resPt);
          if ( pt_denom > 100. && pt_denom < 150. ) hBRPt100_150->Fill (resPt);
          if ( pt_denom > 150. && pt_denom < 200. ) hBRPt150_200->Fill (resPt);
          if ( pt_denom > 200. && pt_denom < 250. ) hBRPt200_250->Fill (resPt);
          if ( pt_denom > 250. && pt_denom < 300. ) hBRPt250_300->Fill (resPt);
          if ( pt_denom > 300. && pt_denom < 400. ) hBRPt300_400->Fill (resPt);
          if ( pt_denom > 400. && pt_denom < 500. ) hBRPt400_500->Fill (resPt);
          if ( pt_denom > 500. && pt_denom < 750. ) hBRPt500_750->Fill (resPt);
          if ( pt_denom > 750. && pt_denom < 1250. ) hBRPt750_1250->Fill (resPt);
          if ( pt_denom > 1250. && pt_denom < 2000. ) hBRPt1250_2000->Fill (resPt);
          if ( pt_denom > 2000. && pt_denom < 5000. ) hBRPt2000_5000->Fill (resPt);
          hBNCH->Fill(rec_ChargedMultiplicity);
          hBNCH0vsPt->Fill(pt_denom,chMult[0]);
          hBNCH1vsPt->Fill(pt_denom,chMult[1]);
          hBNCH2vsPt->Fill(pt_denom,chMult[2]);
          hBNCH3vsPt->Fill(pt_denom,chMult[3]);
          hBNCH4vsPt->Fill(pt_denom,chMult[4]);
          hBNCH5vsPt->Fill(pt_denom,chMult[5]);
          hBNCH6vsPt->Fill(pt_denom,chMult[6]);
          hBNCHvsPt->Fill(pt_denom,rec_ChargedMultiplicity);
          if ( rec_eta > 0. ) 
            hBDEtavsPt->Fill(pt_denom,rec_eta-true_eta);
          else
            hBDEtavsPt->Fill(pt_denom,-rec_eta+true_eta);
          hBDPhivsPt->Fill(pt_denom,rec_phi-true_phi);
        }
        if(plot2)hBRCHE->Fill(resChargedHadEnergy);
        if(plot3)hBRNHE->Fill(resNeutralHadEnergy);
        if(plot4)hBRNEE->Fill(resNeutralEmEnergy);
        if(plot5)hBRneut->Fill(resNeutralEnergy);
        if(plot1)hBRPtvsPt->Fill(pt_denom, resPt);
        if(plot2)hBRCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
        if(plot3)hBRNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
        if(plot4)hBRNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
        if(plot5)hBRneutvsPt->Fill(pt_denom, resNeutralEnergy);
        if(plot6) { 
          hBRHCALvsP->Fill(true_E, resHCALEnergy);
          hBRNEUTvsP->Fill(true_E, resNEUTEnergy);
        }
        if(plot7) { 
          hBRHONLvsP->Fill(true_E, resHCALEnergy);
          hBRNONLvsP->Fill(true_E, resNEUTEnergy);
        }

      }
      // delta Pt or E quantities for Endcap
      if (Endcap){
        if(plot1) {
          hERPt->Fill (resPt);
          if ( pt_denom >  20. && pt_denom <  40. ) hERPt20_40->Fill (resPt);
          if ( pt_denom >  40. && pt_denom <  60. ) hERPt40_60->Fill (resPt);
          if ( pt_denom >  60. && pt_denom <  80. ) hERPt60_80->Fill (resPt);
          if ( pt_denom >  80. && pt_denom < 100. ) hERPt80_100->Fill (resPt);
          if ( pt_denom > 100. && pt_denom < 150. ) hERPt100_150->Fill (resPt);
          if ( pt_denom > 150. && pt_denom < 200. ) hERPt150_200->Fill (resPt);
          if ( pt_denom > 200. && pt_denom < 250. ) hERPt200_250->Fill (resPt);
          if ( pt_denom > 250. && pt_denom < 300. ) hERPt250_300->Fill (resPt);
          if ( pt_denom > 300. && pt_denom < 400. ) hERPt300_400->Fill (resPt);
          if ( pt_denom > 400. && pt_denom < 500. ) hERPt400_500->Fill (resPt);
          if ( pt_denom > 500. && pt_denom < 750. ) hERPt500_750->Fill (resPt);
          if ( pt_denom > 750. && pt_denom < 1250. ) hERPt750_1250->Fill (resPt);
          if ( pt_denom > 1250. && pt_denom < 2000. ) hERPt1250_2000->Fill (resPt);
          if ( pt_denom > 2000. && pt_denom < 5000. ) hERPt2000_5000->Fill (resPt);
          hENCH->Fill(rec_ChargedMultiplicity);
          hENCHvsPt->Fill(pt_denom,rec_ChargedMultiplicity);
          hENCH0vsPt->Fill(pt_denom,chMult[0]);
          hENCH1vsPt->Fill(pt_denom,chMult[1]);
          hENCH2vsPt->Fill(pt_denom,chMult[2]);
          hENCH3vsPt->Fill(pt_denom,chMult[3]);
          hENCH4vsPt->Fill(pt_denom,chMult[4]);
          hENCH5vsPt->Fill(pt_denom,chMult[5]);
          hENCH6vsPt->Fill(pt_denom,chMult[6]);
          if ( rec_eta > 0. ) 
            hEDEtavsPt->Fill(pt_denom,rec_eta-true_eta);
          else
            hEDEtavsPt->Fill(pt_denom,-rec_eta+true_eta);
          hEDPhivsPt->Fill(pt_denom,rec_phi-true_phi);
        }
        if(plot2)hERCHE->Fill(resChargedHadEnergy);
        if(plot3)hERNHE->Fill(resNeutralHadEnergy);
        if(plot4)hERNEE->Fill(resNeutralEmEnergy);
        if(plot5)hERneut->Fill(resNeutralEnergy);
        if(plot1)hERPtvsPt->Fill(pt_denom, resPt);
        if(plot2)hERCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
        if(plot3)hERNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
        if(plot4)hERNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
        if(plot5)hERneutvsPt->Fill(pt_denom, resNeutralEnergy);
        if(plot6) {
          hERHCALvsP->Fill(true_E, resHCALEnergy);
          hERNEUTvsP->Fill(true_E, resNEUTEnergy);
        }
        if(plot7) {
          hERHONLvsP->Fill(true_E, resHCALEnergy);
          hERNONLvsP->Fill(true_E, resNEUTEnergy);
        }
      }                                         
      // delta Pt or E quantities for Forward
      if (Forward){
        if(plot1) {
          hFRPt->Fill (resPt);
          if ( pt_denom >  20. && pt_denom <  40. ) hFRPt20_40->Fill (resPt);
          if ( pt_denom >  40. && pt_denom <  60. ) hFRPt40_60->Fill (resPt);
          if ( pt_denom >  60. && pt_denom <  80. ) hFRPt60_80->Fill (resPt);
          if ( pt_denom >  80. && pt_denom < 100. ) hFRPt80_100->Fill (resPt);
          if ( pt_denom > 100. && pt_denom < 150. ) hFRPt100_150->Fill (resPt);
          if ( pt_denom > 150. && pt_denom < 200. ) hFRPt150_200->Fill (resPt);
          if ( pt_denom > 200. && pt_denom < 250. ) hFRPt200_250->Fill (resPt);
          if ( pt_denom > 250. && pt_denom < 300. ) hFRPt250_300->Fill (resPt);
          if ( pt_denom > 300. && pt_denom < 400. ) hFRPt300_400->Fill (resPt);
          if ( pt_denom > 400. && pt_denom < 500. ) hFRPt400_500->Fill (resPt);
          if ( pt_denom > 500. && pt_denom < 750. ) hFRPt500_750->Fill (resPt);
          if ( pt_denom > 750. && pt_denom < 1250. ) hFRPt750_1250->Fill (resPt);
          if ( pt_denom > 1250. && pt_denom < 2000. ) hFRPt1250_2000->Fill (resPt);
          if ( pt_denom > 2000. && pt_denom < 5000. ) hFRPt2000_5000->Fill (resPt);
          if ( rec_eta > 0. ) 
            hFDEtavsPt->Fill(pt_denom,rec_eta-true_eta);
          else
            hFDEtavsPt->Fill(pt_denom,-rec_eta+true_eta);
          hFDPhivsPt->Fill(pt_denom,rec_phi-true_phi);
        }
        if(plot2)hFRCHE->Fill(resChargedHadEnergy);
        if(plot3)hFRNHE->Fill(resNeutralHadEnergy);
        if(plot4)hFRNEE->Fill(resNeutralEmEnergy);
        if(plot5)hFRneut->Fill(resNeutralEnergy);
        if(plot1)hFRPtvsPt->Fill(pt_denom, resPt);
        if(plot2)hFRCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
        if(plot3)hFRNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
        if(plot4)hFRNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
        if(plot5)hFRneutvsPt->Fill(pt_denom, resNeutralEnergy);
        if(plot6) {
          hFRHCALvsP->Fill(true_E, resHCALEnergy);
          hFRNEUTvsP->Fill(true_E, resNEUTEnergy);
        }
        if(plot7) {
          hFRHONLvsP->Fill(true_E, resHCALEnergy);
          hFRNONLvsP->Fill(true_E, resNEUTEnergy);
        }
      }                                         
    } // end case deltaR < deltaRMax
                
  } // i loop on pf Jets        

  // Increment counter
  entry_++;
}

double PFJetBenchmark::resChargedHadEnergyMax ( ) const [inline]

Definition at line 62 of file PFJetBenchmark.h.

References resChargedHadEnergyMax_.

Referenced by PFRootEventManager::processEntry().

{return resChargedHadEnergyMax_;};

double PFJetBenchmark::resNeutralEmEnergyMax ( ) const [inline]

Definition at line 64 of file PFJetBenchmark.h.

References resNeutralEmEnergyMax_.

Referenced by PFRootEventManager::processEntry().

{return resNeutralEmEnergyMax_;};

double PFJetBenchmark::resNeutralHadEnergyMax ( ) const [inline]

Definition at line 63 of file PFJetBenchmark.h.

References resNeutralHadEnergyMax_.

Referenced by PFRootEventManager::processEntry().

{return resNeutralHadEnergyMax_;};

double PFJetBenchmark::resPtMax ( ) const [inline]

Definition at line 61 of file PFJetBenchmark.h.

References resPtMax_.

Referenced by PFRootEventManager::processEntry().

{return resPtMax_;};

void PFJetBenchmark::setup	(	std::string	Filename,
		bool	debug,
		bool	plotAgainstReco = `0`,
		bool	onlyTwoJets = `1`,
		double	deltaRMax = `0.1`,
		std::string	benchmarkLabel_ = `"ParticleFlow"`,
		double	recPt = `-1`,
		double	maxEta = `-1`,
		DQMStore *	dbe_store = `NULL`
	)

Referenced by PFJetBenchmarkAnalyzer::PFJetBenchmarkAnalyzer(), and PFRootEventManager::readOptions().

void PFJetBenchmark::write ( void )

Definition at line 46 of file PFJetBenchmark.cc.

References gather_cfg::cout, dbe_, file_, outputFile_, and DQMStore::save().

Referenced by PFJetBenchmarkAnalyzer::endJob(), and PFRootEventManager::write().

                           {
   // Store the DAQ Histograms 
  if (outputFile_.size() != 0) {
    if (dbe_)
          dbe_->save(outputFile_.c_str());
    // use bare Root if no DQM (FWLite applications)
    else if (file_) {
       file_->Write(outputFile_.c_str());
       cout << "Benchmark output written to file " << outputFile_.c_str() << endl;
       file_->Close();
       }
  }
  else 
    cout << "No output file specified ("<<outputFile_<<"). Results will not be saved!" << endl;
    
}