CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_2_7_hltpatch1/src/Validation/RecoJets/plugins/PFJetTester.cc

Go to the documentation of this file.
00001 // Producer for validation histograms for CaloJet objects
00002 // F. Ratnikov, Sept. 7, 2006
00003 // Modified by J F Novak July 10, 2008
00004 // $Id: PFJetTester.cc,v 1.11.2.2 2010/08/07 20:42:15 wmtan Exp $
00005 
00006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00007 #include "FWCore/Framework/interface/Event.h"
00008 #include "DataFormats/Common/interface/Handle.h"
00009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00010 #include "FWCore/ServiceRegistry/interface/Service.h"
00011 
00012 #include "DataFormats/Math/interface/deltaR.h"
00013 
00014 #include "DQMServices/Core/interface/DQMStore.h"
00015 #include "DQMServices/Core/interface/MonitorElement.h"
00016 
00017 #include "DataFormats/JetReco/interface/PFJet.h"
00018 #include "DataFormats/JetReco/interface/GenJet.h"
00019 #include "DataFormats/JetReco/interface/GenJetCollection.h"
00020 
00021 //I don't know which of these I actually need yet
00022 #include "DataFormats/METReco/interface/CaloMET.h"
00023 #include "DataFormats/METReco/interface/CaloMETCollection.h"
00024 #include "DataFormats/METReco/interface/GenMET.h"
00025 #include "DataFormats/METReco/interface/GenMETCollection.h"
00026 #include "DataFormats/METReco/interface/MET.h"
00027 #include "DataFormats/METReco/interface/METCollection.h"
00028 
00029 #include "RecoJets/JetProducers/interface/JetMatchingTools.h"
00030 
00031 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
00032 
00033 #include "PFJetTester.h"
00034 
00035 #include <cmath>
00036 
00037 using namespace edm;
00038 using namespace reco;
00039 using namespace std;
00040 
00041 namespace {
00042   bool is_B (const reco::Jet& fJet) {return fabs (fJet.eta()) < 1.3;}
00043   bool is_E (const reco::Jet& fJet) {return fabs (fJet.eta()) >= 1.3 && fabs (fJet.eta()) < 3.;}
00044   bool is_F (const reco::Jet& fJet) {return fabs (fJet.eta()) >= 3.;}
00045 }
00046 
00047 PFJetTester::PFJetTester(const edm::ParameterSet& iConfig)
00048   : mInputCollection (iConfig.getParameter<edm::InputTag>( "src" )),
00049     mInputGenCollection (iConfig.getParameter<edm::InputTag>( "srcGen" )),
00050     mOutputFile (iConfig.getUntrackedParameter<std::string>("outputFile", "")),
00051     mMatchGenPtThreshold (iConfig.getParameter<double>("genPtThreshold")),
00052     mGenEnergyFractionThreshold (iConfig.getParameter<double>("genEnergyFractionThreshold")),
00053     mReverseEnergyFractionThreshold (iConfig.getParameter<double>("reverseEnergyFractionThreshold")),
00054     mRThreshold (iConfig.getParameter<double>("RThreshold")),
00055     mTurnOnEverything (iConfig.getUntrackedParameter<std::string>("TurnOnEverything",""))
00056 {
00057     numberofevents
00058     = mEta = mEtaFineBin = mPhi = mPhiFineBin = mE = mE_80 = mE_3000
00059     = mP = mP_80 = mP_3000 = mPt = mPt_80 = mPt_3000
00060     = mMass = mMass_80 = mMass_3000 = mConstituents = mConstituents_80
00061     = mEtaFirst = mPhiFirst = mEFirst = mEFirst_80 = mEFirst_3000 = mPtFirst = mPtFirst_80 = mPtFirst_3000
00062     = mMjj = mMjj_3000 = mDelEta = mDelPhi = mDelPt 
00063       //    = mMaxEInEmTowers = mMaxEInHadTowers 
00064       //    = mHadEnergyInHO = mHadEnergyInHB = mHadEnergyInHF = mHadEnergyInHE 
00065       //    = mHadEnergyInHO_80 = mHadEnergyInHB_80 = mHadEnergyInHE_80 
00066       //    = mHadEnergyInHO_3000 = mHadEnergyInHB_3000 = mHadEnergyInHE_3000 
00067       //    = mEmEnergyInEB = mEmEnergyInEE = mEmEnergyInHF 
00068       //    = mEmEnergyInEB_80 = mEmEnergyInEE_80
00069       //    = mEmEnergyInEB_3000 = mEmEnergyInEE_3000
00070       //    = mEnergyFractionHadronic = mEnergyFractionEm 
00071       //    = mHFLong = mHFTotal = mHFLong_80 = mHFLong_3000 = mHFShort = mHFShort_80 = mHFShort_3000
00072       //    = mN90
00073     = mChargedEmEnergy_80 = mChargedHadronEnergy_80 = mNeutralEmEnergy_80 = mNeutralHadronEnergy_80
00074     = mChargedEmEnergy_3000 = mChargedHadronEnergy_3000 = mNeutralEmEnergy_3000 = mNeutralHadronEnergy_3000
00075     = mChargedEmEnergyFraction = mChargedHadronEnergyFraction = mNeutralEmEnergyFraction = mNeutralHadronEnergyFraction
00076     = mCaloMEx = mCaloMEx_3000 = mCaloMEy = mCaloMEy_3000 = mCaloMETSig = mCaloMETSig_3000
00077     = mCaloMET = mCaloMET_3000 =  mCaloMETPhi = mCaloSumET  = mCaloSumET_3000   
00078     = mHadTiming = mEmTiming 
00079     = mNJetsEtaC = mNJetsEtaF = mNJets1 = mNJets2
00080     = mAllGenJetsPt = mMatchedGenJetsPt = mAllGenJetsEta = mMatchedGenJetsEta 
00081     = mGenJetMatchEnergyFraction = mReverseMatchEnergyFraction = mRMatch
00082     = mDeltaEta = mDeltaPhi = mEScale = mlinEScale = mDeltaE
00083     = mHadEnergyProfile = mEmEnergyProfile = mJetEnergyProfile 
00084       //    = mHadJetEnergyProfile = mEMJetEnergyProfile
00085     = mEScale_pt10 = mEScaleFineBin
00086     = mpTScaleB_s = mpTScaleE_s = mpTScaleF_s 
00087     = mpTScaleB_d = mpTScaleE_d = mpTScaleF_d
00088     = mpTScale_60_120_s = mpTScale_200_300_s = mpTScale_600_900_s = mpTScale_2700_3500_s
00089     = mpTScale_60_120_d = mpTScale_200_300_d = mpTScale_600_900_d = mpTScale_2700_3500_d
00090     = mpTScale1DB_60_120    = mpTScale1DE_60_120    = mpTScale1DF_60_120 
00091     = mpTScale1DB_200_300   = mpTScale1DE_200_300   = mpTScale1DF_200_300 
00092     = mpTScale1DB_600_900   = mpTScale1DE_600_900   = mpTScale1DF_600_900 
00093     = mpTScale1DB_2700_3500 = mpTScale1DE_2700_3500 = mpTScale1DF_2700_3500
00094     = mpTScale1D_60_120 = mpTScale1D_200_300 = mpTScale1D_600_900 = mpTScale1D_2700_3500
00095     = mHBEne = mHBTime = mHEEne = mHETime = mHFEne = mHFTime = mHOEne = mHOTime
00096     = mEBEne = mEBTime = mEEEne = mEETime 
00097     = mPthat_80 = mPthat_3000
00098     = 0;
00099   
00100   DQMStore* dbe = &*edm::Service<DQMStore>();
00101   if (dbe) {
00102     dbe->setCurrentFolder("JetMET/RecoJetsV/PFJetTask_" + mInputCollection.label());
00103     //
00104     numberofevents    = dbe->book1D("numberofevents","numberofevents", 3, 0 , 2);
00105     //
00106     mEta              = dbe->book1D("Eta", "Eta", 100, -5, 5); 
00107     mEtaFineBin       = dbe->book1D("EtaFineBin_Pt10", "EtaFineBin_Pt10", 500, -5, 5); 
00108     mEtaFineBin1p     = dbe->book1D("EtaFineBin1p_Pt10", "EtaFineBin1p_Pt10", 100, 0, 1.3); 
00109     mEtaFineBin2p     = dbe->book1D("EtaFineBin2p_Pt10", "EtaFineBin2p_Pt10", 100, 1.3, 3); 
00110     mEtaFineBin3p     = dbe->book1D("EtaFineBin3p_Pt10", "EtaFineBin3p_Pt10", 100, 3, 5); 
00111     mEtaFineBin1m     = dbe->book1D("EtaFineBin1m_Pt10", "EtaFineBin1m_Pt10", 100, -1.3, 0); 
00112     mEtaFineBin2m     = dbe->book1D("EtaFineBin2m_Pt10", "EtaFineBin2m_Pt10", 100, -3, -1.3); 
00113     mEtaFineBin3m     = dbe->book1D("EtaFineBin3m_Pt10", "EtaFineBin3m_Pt10", 100, -5, -3); 
00114     //
00115     mPhi              = dbe->book1D("Phi", "Phi", 70, -3.5, 3.5); 
00116     mPhiFineBin       = dbe->book1D("PhiFineBin_Pt10", "PhiFineBin_Pt10", 350, -3.5, 3.5); 
00117     //
00118     mE                = dbe->book1D("E", "E", 100, 0, 500); 
00119     mE_80             = dbe->book1D("E_80", "E_80", 100, 0, 4500); 
00120     mE_3000           = dbe->book1D("E_3000", "E_3000", 100, 0, 6000); 
00121     //
00122     mP                = dbe->book1D("P", "P", 100, 0, 500); 
00123     mP_80             = dbe->book1D("P_80", "P_80", 100, 0, 4500); 
00124     mP_3000           = dbe->book1D("P_3000", "P_3000", 100, 0, 6000); 
00125     //
00126     mPt               = dbe->book1D("Pt", "Pt", 100, 0, 50); 
00127     mPt_80            = dbe->book1D("Pt_80", "Pt_80", 100, 0, 140); 
00128     mPt_3000          = dbe->book1D("Pt_3000", "Pt_3000", 100, 0, 4000); 
00129     //
00130     mMass             = dbe->book1D("Mass", "Mass", 100, 0, 25); 
00131     mMass_80          = dbe->book1D("Mass_80", "Mass_80", 100, 0, 120); 
00132     mMass_3000        = dbe->book1D("Mass_3000", "Mass_3000", 100, 0, 500); 
00133     //
00134     mConstituents     = dbe->book1D("Constituents", "# of Constituents", 100, 0, 100); 
00135     mConstituents_80  = dbe->book1D("Constituents_80", "# of Constituents_80", 40, 0, 40); 
00136     //
00137     mEtaFirst         = dbe->book1D("EtaFirst", "EtaFirst", 100, -5, 5); 
00138     mPhiFirst         = dbe->book1D("PhiFirst", "PhiFirst", 70, -3.5, 3.5);     
00139     mEFirst           = dbe->book1D("EFirst", "EFirst", 100, 0, 1000); 
00140     mEFirst_80        = dbe->book1D("EFirst_80", "EFirst_80", 100, 0, 180); 
00141     mEFirst_3000      = dbe->book1D("EFirst_3000", "EFirst_3000", 100, 0, 4000); 
00142     mPtFirst          = dbe->book1D("PtFirst", "PtFirst", 100, 0, 50); 
00143     mPtFirst_80       = dbe->book1D("PtFirst_80", "PtFirst_80", 100, 0, 140);
00144     mPtFirst_3000     = dbe->book1D("PtFirst_3000", "PtFirst_3000", 100, 0, 4000);
00145     //
00146     mMjj              = dbe->book1D("Mjj", "Mjj", 100, 0, 2000); 
00147     mMjj_3000         = dbe->book1D("Mjj_3000", "Mjj_3000", 100, 0, 10000); 
00148     mDelEta           = dbe->book1D("DelEta", "DelEta", 100, -.5, .5); 
00149     mDelPhi           = dbe->book1D("DelPhi", "DelPhi", 100, -.5, .5); 
00150     mDelPt            = dbe->book1D("DelPt", "DelPt", 100, -1, 1); 
00151     //
00152     //    mMaxEInEmTowers   = dbe->book1D("MaxEInEmTowers", "MaxEInEmTowers", 100, 0, 100); 
00153     //    mMaxEInHadTowers  = dbe->book1D("MaxEInHadTowers", "MaxEInHadTowers", 100, 0, 100); 
00154     //    mHadEnergyInHO    = dbe->book1D("HadEnergyInHO", "HadEnergyInHO", 100, 0, 10); 
00155     //    mHadEnergyInHB    = dbe->book1D("HadEnergyInHB", "HadEnergyInHB", 100, 0, 50); 
00156     //    mHadEnergyInHF    = dbe->book1D("HadEnergyInHF", "HadEnergyInHF", 100, 0, 50); 
00157     //    mHadEnergyInHE    = dbe->book1D("HadEnergyInHE", "HadEnergyInHE", 100, 0, 100); 
00158     //
00159     //    mHadEnergyInHO_80    = dbe->book1D("HadEnergyInHO_80", "HadEnergyInHO_80", 100, 0, 50); 
00160     //    mHadEnergyInHB_80    = dbe->book1D("HadEnergyInHB_80", "HadEnergyInHB_80", 100, 0, 200); 
00161     //    mHadEnergyInHE_80    = dbe->book1D("HadEnergyInHE_80", "HadEnergyInHE_80", 100, 0, 1000); 
00162     //    mHadEnergyInHO_3000  = dbe->book1D("HadEnergyInHO_3000", "HadEnergyInHO_3000", 100, 0, 500); 
00163     //    mHadEnergyInHB_3000  = dbe->book1D("HadEnergyInHB_3000", "HadEnergyInHB_3000", 100, 0, 3000); 
00164     //    mHadEnergyInHE_3000  = dbe->book1D("HadEnergyInHE_3000", "HadEnergyInHE_3000", 100, 0, 2000); 
00165     //
00166     //    mEmEnergyInEB     = dbe->book1D("EmEnergyInEB", "EmEnergyInEB", 100, 0, 50); 
00167     //    mEmEnergyInEE     = dbe->book1D("EmEnergyInEE", "EmEnergyInEE", 100, 0, 50); 
00168     //    mEmEnergyInHF     = dbe->book1D("EmEnergyInHF", "EmEnergyInHF", 120, -20, 100); 
00169     //    mEmEnergyInEB_80  = dbe->book1D("EmEnergyInEB_80", "EmEnergyInEB_80", 100, 0, 200); 
00170     //    mEmEnergyInEE_80  = dbe->book1D("EmEnergyInEE_80", "EmEnergyInEE_80", 100, 0, 1000); 
00171     //    mEmEnergyInEB_3000= dbe->book1D("EmEnergyInEB_3000", "EmEnergyInEB_3000", 100, 0, 3000); 
00172     //    mEmEnergyInEE_3000= dbe->book1D("EmEnergyInEE_3000", "EmEnergyInEE_3000", 100, 0, 2000); 
00173     //    mEnergyFractionHadronic = dbe->book1D("EnergyFractionHadronic", "EnergyFractionHadronic", 120, -0.1, 1.1); 
00174     //    mEnergyFractionEm = dbe->book1D("EnergyFractionEm", "EnergyFractionEm", 120, -0.1, 1.1); 
00175     //
00176     //    mHFTotal          = dbe->book1D("HFTotal", "HFTotal", 100, 0, 500);
00177     //    mHFTotal_80       = dbe->book1D("HFTotal_80", "HFTotal_80", 100, 0, 3000);
00178     //    mHFTotal_3000     = dbe->book1D("HFTotal_3000", "HFTotal_3000", 100, 0, 6000);
00179     //    mHFLong           = dbe->book1D("HFLong", "HFLong", 100, 0, 500);
00180     //    mHFLong_80        = dbe->book1D("HFLong_80", "HFLong_80", 100, 0, 200);
00181     //    mHFLong_3000      = dbe->book1D("HFLong_3000", "HFLong_3000", 100, 0, 1500);
00182     //    mHFShort          = dbe->book1D("HFShort", "HFShort", 100, 0, 500);
00183     //    mHFShort_80       = dbe->book1D("HFShort_80", "HFShort_80", 100, 0, 200);
00184     //    mHFShort_3000     = dbe->book1D("HFShort_3000", "HFShort_3000", 100, 0, 1500);
00185     //
00186     //    mN90              = dbe->book1D("N90", "N90", 50, 0, 50); 
00187     //
00188 
00189     mChargedEmEnergy_80 = dbe->book1D("ChargedEmEnergy_80","ChargedEmEnergy_80",100,0.,500.);
00190     mChargedHadronEnergy_80 = dbe->book1D("ChargedHadronEnergy_80","ChargedHadronEnergy_80",100,0.,500.);
00191     mNeutralEmEnergy_80 = dbe->book1D("NeutralEmEnergy_80","NeutralEmEnergy_80",100,0.,500.);
00192     mNeutralHadronEnergy_80 = dbe->book1D("NeutralHadronEnergy_80","NeutralHadronEnergy_80",100,0.,500.);
00193 
00194     mChargedEmEnergy_3000 = dbe->book1D("ChargedEmEnergy_3000","ChargedEmEnergy_3000",100,0.,2000.);
00195     mChargedHadronEnergy_3000 = dbe->book1D("ChargedHadronEnergy_3000","ChargedHadronEnergy_3000",100,0.,2000.);
00196     mNeutralEmEnergy_3000 = dbe->book1D("NeutralEmEnergy_3000","NeutralEmEnergy_3000",100,0.,2000.);
00197     mNeutralHadronEnergy_3000 = dbe->book1D("NeutralHadronEnergy_3000","NeutralHadronEnergy_3000",100,0.,2000.);
00198 
00199     mChargedEmEnergyFraction = dbe->book1D("ChargedEmEnergyFraction","ChargedEmEnergyFraction",120,-0.1,1.1);
00200     mChargedHadronEnergyFraction = dbe->book1D("ChargedHadronEnergyFraction","ChargedHadronEnergyFraction",120,-0.1,1.1);
00201     mNeutralEmEnergyFraction = dbe->book1D("NeutralEmEnergyFraction","NeutralEmEnergyFraction",120,-0.1,1.1);
00202     mNeutralHadronEnergyFraction = dbe->book1D("NeutralHadronEnergyFraction","NeutralHadronEnergyFraction",120,-0.1,1.1);
00203 
00204 
00205     mGenEta           = dbe->book1D("GenEta", "GenEta", 100, -5, 5);
00206     mGenPhi           = dbe->book1D("GenPhi", "GenPhi", 70, -3.5, 3.5);
00207     mGenPt            = dbe->book1D("GenPt", "GenPt", 100, 0, 50);
00208     mGenPt_80         = dbe->book1D("GenPt_80", "GenPt_80", 100, 0, 140);
00209     mGenPt_3000       = dbe->book1D("GenPt_3000", "GenPt_3000", 100, 0, 1500);
00210     //
00211     mGenEtaFirst      = dbe->book1D("GenEtaFirst", "GenEtaFirst", 100, -5, 5);
00212     mGenPhiFirst      = dbe->book1D("GenPhiFirst", "GenPhiFirst", 70, -3.5, 3.5);
00213     //
00214     mCaloMEx          = dbe->book1D("CaloMEx","CaloMEx",200,-150,150);
00215     mCaloMEx_3000     = dbe->book1D("CaloMEx_3000","CaloMEx_3000",100,-500,500);
00216     mCaloMEy          = dbe->book1D("CaloMEy","CaloMEy",200,-150,150);
00217     mCaloMEy_3000     = dbe->book1D("CaloMEy_3000","CaloMEy_3000",100,-500,500);
00218     mCaloMETSig       = dbe->book1D("CaloMETSig","CaloMETSig",100,0,15);
00219     mCaloMETSig_3000  = dbe->book1D("CaloMETSig_3000","CaloMETSig_3000",100,0,50);
00220     mCaloMET          = dbe->book1D("CaloMET","CaloMET",100,0,150);
00221     mCaloMET_3000     = dbe->book1D("CaloMET_3000","CaloMET_3000",100,0,1000);
00222     mCaloMETPhi       = dbe->book1D("CaloMETPhi","CaloMETPhi",70, -3.5, 3.5);
00223     mCaloSumET        = dbe->book1D("CaloSumET","CaloSumET",100,0,500);
00224     mCaloSumET_3000   = dbe->book1D("CaloSumET_3000","CaloSumET_3000",100,3000,8000);
00225     //
00226     mHadTiming        = dbe->book1D("HadTiming", "HadTiming", 75, -50, 100);
00227     mEmTiming         = dbe->book1D("EMTiming", "EMTiming", 75, -50, 100);
00228     //
00229     mNJetsEtaC        = dbe->book1D("NJetsEtaC_Pt10", "NJetsEtaC_Pt10", 15, 0, 15);
00230     mNJetsEtaF        = dbe->book1D("NJetsEtaF_Pt10", "NJetsEtaF_Pt10", 15, 0, 15);
00231     //
00232     mNJets1           = dbe->bookProfile("NJets1", "NJets1", 100, 0, 200,  100, 0, 50, "s");
00233     mNJets2           = dbe->bookProfile("NJets2", "NJets2", 100, 0, 4000, 100, 0, 50, "s");
00234     //
00235     mHBEne     = dbe->book1D( "HBEne",  "HBEne", 1000, -20, 100 );
00236     mHBTime    = dbe->book1D( "HBTime", "HBTime", 200, -200, 200 );
00237     mHEEne     = dbe->book1D( "HEEne",  "HEEne", 1000, -20, 100 );
00238     mHETime    = dbe->book1D( "HETime", "HETime", 200, -200, 200 );
00239     mHOEne     = dbe->book1D( "HOEne",  "HOEne", 1000, -20, 100 );
00240     mHOTime    = dbe->book1D( "HOTime", "HOTime", 200, -200, 200 );
00241     mHFEne     = dbe->book1D( "HFEne",  "HFEne", 1000, -20, 100 );
00242     mHFTime    = dbe->book1D( "HFTime", "HFTime", 200, -200, 200 );
00243     mEBEne     = dbe->book1D( "EBEne",  "EBEne", 1000, -20, 100 );
00244     mEBTime    = dbe->book1D( "EBTime", "EBTime", 200, -200, 200 );
00245     mEEEne     = dbe->book1D( "EEEne",  "EEEne", 1000, -20, 100 );
00246     mEETime    = dbe->book1D( "EETime", "EETime", 200, -200, 200 );
00247     //
00248     mPthat_80            = dbe->book1D("Pthat_80", "Pthat_80", 100, 40.0, 160.0); 
00249     mPthat_3000          = dbe->book1D("Pthat_3000", "Pthat_3000", 100, 2500.0, 4000.0); 
00250     //
00251     double log10PtMin = 0.5; //=3.1622766
00252     double log10PtMax = 3.75; //=5623.41325
00253     int log10PtBins = 26; 
00254     double etaMin = -5.;
00255     double etaMax = 5.;
00256     int etaBins = 50;
00257 
00258     double linPtMin = 5;
00259     double linPtMax = 155;
00260     int linPtBins = 15;
00261 
00262     int log10PtFineBins = 50;
00263 
00264     mAllGenJetsPt = dbe->book1D("GenJetLOGpT", "GenJet LOG(pT_gen)", 
00265                                 log10PtBins, log10PtMin, log10PtMax);
00266     mMatchedGenJetsPt = dbe->book1D("MatchedGenJetLOGpT", "MatchedGenJet LOG(pT_gen)", 
00267                                     log10PtBins, log10PtMin, log10PtMax);
00268     mAllGenJetsEta = dbe->book2D("GenJetEta", "GenJet Eta vs LOG(pT_gen)", 
00269                                  log10PtBins, log10PtMin, log10PtMax, etaBins, etaMin, etaMax);
00270     mMatchedGenJetsEta = dbe->book2D("MatchedGenJetEta", "MatchedGenJet Eta vs LOG(pT_gen)", 
00271                                      log10PtBins, log10PtMin, log10PtMax, etaBins, etaMin, etaMax);
00272     //
00273     if (mTurnOnEverything.compare("yes")==0) {
00274       mHadEnergyProfile = dbe->bookProfile2D("HadEnergyProfile", "HadEnergyProfile", 82, -41, 41, 73, 0, 73, 100, 0, 10000, "s");
00275       mEmEnergyProfile  = dbe->bookProfile2D("EmEnergyProfile", "EmEnergyProfile", 82, -41, 41, 73, 0, 73, 100, 0, 10000, "s");
00276     }
00277     mJetEnergyProfile = dbe->bookProfile2D("JetEnergyProfile", "JetEnergyProfile", 50, -5, 5, 36, -3.1415987, 3.1415987, 100, 0, 10000, "s");
00278     //    mHadJetEnergyProfile = dbe->bookProfile2D("HadJetEnergyProfile", "HadJetEnergyProfile", 50, -5, 5, 36, -3.1415987, 3.1415987, 100, 0, 10000, "s");
00279     //    mEMJetEnergyProfile = dbe->bookProfile2D("EMJetEnergyProfile", "EMJetEnergyProfile", 50, -5, 5, 36, -3.1415987, 3.1415987, 100, 0, 10000, "s");
00280     //
00281     if (mTurnOnEverything.compare("yes")==0) {
00282     mGenJetMatchEnergyFraction  = dbe->book3D("GenJetMatchEnergyFraction", "GenJetMatchEnergyFraction vs LOG(pT_gen) vs eta", 
00283                                               log10PtBins, log10PtMin, log10PtMax, etaBins, etaMin, etaMax, 101, 0, 1.01);
00284     mReverseMatchEnergyFraction  = dbe->book3D("ReverseMatchEnergyFraction", "ReverseMatchEnergyFraction vs LOG(pT_gen) vs eta", 
00285                                                log10PtBins, log10PtMin, log10PtMax, etaBins, etaMin, etaMax, 101, 0, 1.01);
00286     mRMatch  = dbe->book3D("RMatch", "delta(R)(Gen-Calo) vs LOG(pT_gen) vs eta", 
00287                            log10PtBins, log10PtMin, log10PtMax, etaBins, etaMin, etaMax, 60, 0, 3);
00288     mDeltaEta = dbe->book3D("DeltaEta", "DeltaEta vs LOG(pT_gen) vs eta", 
00289                               log10PtBins, log10PtMin, log10PtMax, etaBins, etaMin, etaMax, 100, -1, 1);
00290     mDeltaPhi = dbe->book3D("DeltaPhi", "DeltaPhi vs LOG(pT_gen) vs eta", 
00291                               log10PtBins, log10PtMin, log10PtMax, etaBins, etaMin, etaMax, 100, -1, 1);
00292     mEScale = dbe->book3D("EScale", "EnergyScale vs LOG(pT_gen) vs eta", 
00293                             log10PtBins, log10PtMin, log10PtMax, etaBins, etaMin, etaMax, 100, 0, 2);
00294     mlinEScale = dbe->book3D("linEScale", "EnergyScale vs LOG(pT_gen) vs eta", 
00295                             linPtBins, linPtMin, linPtMax, etaBins, etaMin, etaMax, 100, 0, 2);
00296     mDeltaE = dbe->book3D("DeltaE", "DeltaE vs LOG(pT_gen) vs eta", 
00297                             log10PtBins, log10PtMin, log10PtMax, etaBins, etaMin, etaMax, 2000, -200, 200);
00298     //
00299     mEScale_pt10 = dbe->book3D("EScale_pt10", "EnergyScale vs LOG(pT_gen) vs eta", 
00300                             log10PtBins, log10PtMin, log10PtMax, etaBins, etaMin, etaMax, 100, 0, 2);
00301     mEScaleFineBin = dbe->book3D("EScaleFineBins", "EnergyScale vs LOG(pT_gen) vs eta", 
00302                             log10PtFineBins, log10PtMin, log10PtMax, etaBins, etaMin, etaMax, 100, 0, 2);
00303     }
00304 
00305     mpTScaleB_s = dbe->bookProfile("pTScaleB_s", "pTScale_s_0<|eta|<1.3",
00306                                     log10PtBins, log10PtMin, log10PtMax, 0, 2, "s");
00307     mpTScaleE_s = dbe->bookProfile("pTScaleE_s", "pTScale_s_1.3<|eta|<3.0",
00308                                     log10PtBins, log10PtMin, log10PtMax, 0, 2, "s");
00309     mpTScaleF_s = dbe->bookProfile("pTScaleF_s", "pTScale_s_3.0<|eta|<5.0",
00310                                     log10PtBins, log10PtMin, log10PtMax, 0, 2, "s");
00311     mpTScaleB_d = dbe->bookProfile("pTScaleB_d", "pTScale_d_0<|eta|<1.3",
00312                                    log10PtBins, log10PtMin, log10PtMax, 0, 2, " ");
00313     mpTScaleE_d = dbe->bookProfile("pTScaleE_d", "pTScale_d_1.3<|eta|<3.0",
00314                                    log10PtBins, log10PtMin, log10PtMax, 0, 2, " ");
00315     mpTScaleF_d = dbe->bookProfile("pTScaleF_d", "pTScale_d_3.0<|eta|<5.0",
00316                                    log10PtBins, log10PtMin, log10PtMax, 0, 2, " ");
00317 
00318     mpTScale_60_120_s    = dbe->bookProfile("pTScale_60_120_s", "pTScale_s_60<pT<120",
00319                                           etaBins, etaMin, etaMax, 0., 2., "s");
00320     mpTScale_200_300_s   = dbe->bookProfile("pTScale_200_300_s", "pTScale_s_200<pT<300",
00321                                           etaBins, etaMin, etaMax, 0., 2., "s");
00322     mpTScale_600_900_s   = dbe->bookProfile("pTScale_600_900_s", "pTScale_s_600<pT<900",
00323                                           etaBins, etaMin, etaMax, 0., 2., "s");
00324     mpTScale_2700_3500_s = dbe->bookProfile("pTScale_2700_3500_s", "pTScale_s_2700<pt<3500",
00325                                           etaBins, etaMin, etaMax, 0., 2., "s");
00326     mpTScale_60_120_d    = dbe->bookProfile("pTScale_60_120_d", "pTScale_d_60<pT<120",
00327                                           etaBins, etaMin, etaMax, 0., 2., " ");
00328     mpTScale_200_300_d   = dbe->bookProfile("pTScale_200_300_d", "pTScale_d_200<pT<300",
00329                                           etaBins, etaMin, etaMax, 0., 2., " ");
00330     mpTScale_600_900_d   = dbe->bookProfile("pTScale_600_900_d", "pTScale_d_600<pT<900",
00331                                           etaBins, etaMin, etaMax, 0., 2., " ");
00332     mpTScale_2700_3500_d = dbe->bookProfile("pTScale_2700_3500_d", "pTScale_d_2700<pt<3500",
00333                                           etaBins, etaMin, etaMax, 0., 2., " ");
00334 
00335     mpTScale1DB_60_120 = dbe->book1D("pTScale1DB_60_120", "pTScale_distribution_for_0<|eta|<1.3_60_120",
00336                                    100, 0, 2);
00337     mpTScale1DE_60_120 = dbe->book1D("pTScale1DE_60_120", "pTScale_distribution_for_1.3<|eta|<3.0_60_120",
00338                                    100, 0, 2);
00339     mpTScale1DF_60_120 = dbe->book1D("pTScale1DF_60_120", "pTScale_distribution_for_3.0<|eta|<5.0_60_120",
00340                                    100, 0, 2);
00341 
00342     mpTScale1DB_200_300 = dbe->book1D("pTScale1DB_200_300", "pTScale_distribution_for_0<|eta|<1.3_200_300",
00343                                    100, 0, 2);
00344     mpTScale1DE_200_300 = dbe->book1D("pTScale1DE_200_300", "pTScale_distribution_for_1.3<|eta|<3.0_200_300",
00345                                    100, 0, 2);
00346     mpTScale1DF_200_300 = dbe->book1D("pTScale1DF_200_300", "pTScale_distribution_for_3.0<|eta|<5.0_200_300",
00347                                    100, 0, 2);
00348 
00349     mpTScale1DB_600_900 = dbe->book1D("pTScale1DB_600_900", "pTScale_distribution_for_0<|eta|<1.3_600_900",
00350                                    100, 0, 2);
00351     mpTScale1DE_600_900 = dbe->book1D("pTScale1DE_600_900", "pTScale_distribution_for_1.3<|eta|<3.0_600_900",
00352                                    100, 0, 2);
00353     mpTScale1DF_600_900 = dbe->book1D("pTScale1DF_600_900", "pTScale_distribution_for_3.0<|eta|<5.0_600_900",
00354                                    100, 0, 2);
00355 
00356     mpTScale1DB_2700_3500 = dbe->book1D("pTScale1DB_2700_3500", "pTScale_distribution_for_0<|eta|<1.3_2700_3500",
00357                                    100, 0, 2);
00358     mpTScale1DE_2700_3500 = dbe->book1D("pTScale1DE_2700_3500", "pTScale_distribution_for_1.3<|eta|<3.0_2700_3500",
00359                                    100, 0, 2);
00360     mpTScale1DF_2700_3500 = dbe->book1D("pTScale1DF_2700_3500", "pTScale_distribution_for_3.0<|eta|<5.0_2700_3500",
00361                                    100, 0, 2);
00362 
00363     mpTScale1D_60_120    = dbe->book1D("pTScale1D_60_120", "pTScale_distribution_for_60<pT<120",
00364                                             100, 0, 2);
00365     mpTScale1D_200_300    = dbe->book1D("pTScale1D_200_300", "pTScale_distribution_for_200<pT<300",
00366                                             100, 0, 2);
00367     mpTScale1D_600_900    = dbe->book1D("pTScale1D_600_900", "pTScale_distribution_for_600<pT<900",
00368                                             100, 0, 2);
00369     mpTScale1D_2700_3500 = dbe->book1D("pTScale1D_2700_3500", "pTScale_distribution_for_2700<pt<3500",
00370                                             100, 0, 2);
00371 
00372   }
00373 
00374   if (mOutputFile.empty ()) {
00375     LogInfo("OutputInfo") << " PFJet histograms will NOT be saved";
00376   } 
00377   else {
00378     LogInfo("OutputInfo") << " PFJethistograms will be saved to file:" << mOutputFile;
00379   }
00380 }
00381    
00382 PFJetTester::~PFJetTester()
00383 {
00384 }
00385 
00386 void PFJetTester::beginJob(){
00387 }
00388 
00389 void PFJetTester::endJob() {
00390  if (!mOutputFile.empty() && &*edm::Service<DQMStore>()) edm::Service<DQMStore>()->save (mOutputFile);
00391 }
00392 
00393 
00394 void PFJetTester::analyze(const edm::Event& mEvent, const edm::EventSetup& mSetup)
00395 {
00396   double countsfornumberofevents = 1;
00397   numberofevents->Fill(countsfornumberofevents);
00398   // *********************************
00399   // *** Get pThat
00400   // *********************************
00401 
00402   edm::Handle<HepMCProduct> evt;
00403   mEvent.getByLabel("generator", evt);
00404   if (evt.isValid()) {
00405   HepMC::GenEvent * myGenEvent = new HepMC::GenEvent(*(evt->GetEvent()));
00406   
00407   double pthat = myGenEvent->event_scale();
00408 
00409   mPthat_80->Fill(pthat);
00410   mPthat_3000->Fill(pthat);
00411 
00412   delete myGenEvent; 
00413   }
00414 
00415   // ***********************************
00416   // *** Get CaloMET
00417   // ***********************************
00418 
00419   const CaloMET *calomet;
00420   edm::Handle<CaloMETCollection> calo;
00421   mEvent.getByLabel("met", calo);
00422   if (!calo.isValid()) {
00423     edm::LogInfo("OutputInfo") << " failed to retrieve data required by MET Task";
00424     edm::LogInfo("OutputInfo") << " MET Task cannot continue...!";
00425   } else {
00426     const CaloMETCollection *calometcol = calo.product();
00427     calomet = &(calometcol->front());
00428     
00429     double caloSumET = calomet->sumEt();
00430     double caloMETSig = calomet->mEtSig();
00431     double caloMET = calomet->pt();
00432     double caloMEx = calomet->px();
00433     double caloMEy = calomet->py();
00434     double caloMETPhi = calomet->phi();
00435 
00436     mCaloMEx->Fill(caloMEx);
00437     mCaloMEx_3000->Fill(caloMEx);
00438     mCaloMEy->Fill(caloMEy);
00439     mCaloMEy_3000->Fill(caloMEy);
00440     mCaloMET->Fill(caloMET);
00441     mCaloMET_3000->Fill(caloMET);
00442     mCaloMETPhi->Fill(caloMETPhi);
00443     mCaloSumET->Fill(caloSumET);
00444     mCaloSumET_3000->Fill(caloSumET);
00445     mCaloMETSig->Fill(caloMETSig);
00446     mCaloMETSig_3000->Fill(caloMETSig);
00447   }
00448 
00449   // ***********************************
00450   // *** Get the CaloTower collection
00451   // ***********************************
00452   Handle<CaloTowerCollection> caloTowers;
00453   mEvent.getByLabel( "towerMaker", caloTowers );
00454   if (caloTowers.isValid()) {
00455   for( CaloTowerCollection::const_iterator cal = caloTowers->begin(); cal != caloTowers->end(); ++ cal ){
00456 
00457     //To compensate for the index
00458     if (mTurnOnEverything.compare("yes")==0) {
00459       if (cal->ieta() >> 0 ){mHadEnergyProfile->Fill (cal->ieta()-1, cal->iphi(), cal->hadEnergy());
00460       mEmEnergyProfile->Fill (cal->ieta()-1, cal->iphi(), cal->emEnergy());}
00461       mHadEnergyProfile->Fill (cal->ieta(), cal->iphi(), cal->hadEnergy());
00462       mEmEnergyProfile->Fill (cal->ieta(), cal->iphi(), cal->emEnergy());
00463     }
00464 
00465     mHadTiming->Fill (cal->hcalTime());
00466     mEmTiming->Fill (cal->ecalTime());    
00467   }
00468   }  
00469 
00470   // ***********************************
00471   // *** Get the RecHits collection
00472   // ***********************************
00473   try {
00474     std::vector<edm::Handle<HBHERecHitCollection> > colls;
00475     mEvent.getManyByType(colls);
00476     std::vector<edm::Handle<HBHERecHitCollection> >::iterator i;
00477     for (i=colls.begin(); i!=colls.end(); i++) {
00478       for (HBHERecHitCollection::const_iterator j=(*i)->begin(); j!=(*i)->end(); j++) {
00479         //      std::cout << *j << std::endl;
00480         if (j->id().subdet() == HcalBarrel) {
00481           mHBEne->Fill(j->energy()); 
00482           mHBTime->Fill(j->time()); 
00483         }
00484         if (j->id().subdet() == HcalEndcap) {
00485           mHEEne->Fill(j->energy()); 
00486           mHETime->Fill(j->time()); 
00487         }
00488 
00489       }
00490     }
00491   } catch (...) {
00492     edm::LogInfo("OutputInfo") << " No HB/HE RecHits.";
00493   }
00494 
00495   try {
00496     std::vector<edm::Handle<HFRecHitCollection> > colls;
00497     mEvent.getManyByType(colls);
00498     std::vector<edm::Handle<HFRecHitCollection> >::iterator i;
00499     for (i=colls.begin(); i!=colls.end(); i++) {
00500       for (HFRecHitCollection::const_iterator j=(*i)->begin(); j!=(*i)->end(); j++) {
00501         //      std::cout << *j << std::endl;
00502         if (j->id().subdet() == HcalForward) {
00503           mHFEne->Fill(j->energy()); 
00504           mHFTime->Fill(j->time()); 
00505         }
00506       }
00507     }
00508   } catch (...) {
00509     edm::LogInfo("OutputInfo") << " No HF RecHits.";
00510   }
00511 
00512   try {
00513     std::vector<edm::Handle<HORecHitCollection> > colls;
00514     mEvent.getManyByType(colls);
00515     std::vector<edm::Handle<HORecHitCollection> >::iterator i;
00516     for (i=colls.begin(); i!=colls.end(); i++) {
00517       for (HORecHitCollection::const_iterator j=(*i)->begin(); j!=(*i)->end(); j++) {
00518         if (j->id().subdet() == HcalOuter) {
00519           mHOEne->Fill(j->energy()); 
00520           mHOTime->Fill(j->time()); 
00521         }
00522       }
00523     }
00524   } catch (...) {
00525     edm::LogInfo("OutputInfo") << " No HO RecHits.";
00526   }
00527   try {
00528     std::vector<edm::Handle<EBRecHitCollection> > colls;
00529     mEvent.getManyByType(colls);
00530     std::vector<edm::Handle<EBRecHitCollection> >::iterator i;
00531     for (i=colls.begin(); i!=colls.end(); i++) {
00532       for (EBRecHitCollection::const_iterator j=(*i)->begin(); j!=(*i)->end(); j++) {
00533         //      if (j->id() == EcalBarrel) {
00534         mEBEne->Fill(j->energy()); 
00535         mEBTime->Fill(j->time()); 
00536         //    }
00537         //      std::cout << *j << std::endl;
00538         //      std::cout << j->id() << std::endl;
00539       }
00540     }
00541   } catch (...) {
00542     edm::LogInfo("OutputInfo") << " No EB RecHits.";
00543   }
00544 
00545   try {
00546     std::vector<edm::Handle<EERecHitCollection> > colls;
00547     mEvent.getManyByType(colls);
00548     std::vector<edm::Handle<EERecHitCollection> >::iterator i;
00549     for (i=colls.begin(); i!=colls.end(); i++) {
00550       for (EERecHitCollection::const_iterator j=(*i)->begin(); j!=(*i)->end(); j++) {
00551         //      if (j->id().subdet() == EcalEndcap) {
00552         mEEEne->Fill(j->energy()); 
00553         mEETime->Fill(j->time()); 
00554         //    }
00555         //      std::cout << *j << std::endl;
00556       }
00557     }
00558   } catch (...) {
00559     edm::LogInfo("OutputInfo") << " No EE RecHits.";
00560   }
00561 
00562 
00563   //***********************************
00564   //*** Get the Jet collection
00565   //***********************************
00566   math::XYZTLorentzVector p4tmp[2];
00567   Handle<PFJetCollection> pfJets;
00568   mEvent.getByLabel(mInputCollection, pfJets);
00569   if (!pfJets.isValid()) return;
00570   PFJetCollection::const_iterator jet = pfJets->begin ();
00571   int jetIndex = 0;
00572   int nJet = 0;
00573   int nJetF = 0;
00574   int nJetC = 0;
00575   for (; jet != pfJets->end (); jet++, jetIndex++) {
00576     if (mEta) mEta->Fill (jet->eta());
00577 
00578     if (jet->pt() > 10.) {
00579       if (fabs(jet->eta()) > 1.3) 
00580         nJetF++;
00581       else 
00582         nJetC++;          
00583     }
00584     if (jet->pt() > 10.) {
00585       if (mEtaFineBin) mEtaFineBin->Fill (jet->eta());
00586       if (mEtaFineBin1p) mEtaFineBin1p->Fill (jet->eta());
00587       if (mEtaFineBin2p) mEtaFineBin2p->Fill (jet->eta());
00588       if (mEtaFineBin3p) mEtaFineBin3p->Fill (jet->eta());
00589       if (mEtaFineBin1m) mEtaFineBin1m->Fill (jet->eta());
00590       if (mEtaFineBin2m) mEtaFineBin2m->Fill (jet->eta());
00591       if (mEtaFineBin3m) mEtaFineBin3m->Fill (jet->eta());
00592       if (mPhiFineBin) mPhiFineBin->Fill (jet->phi());
00593     }
00594     if (mPhi) mPhi->Fill (jet->phi());
00595     if (mE) mE->Fill (jet->energy());
00596     if (mE_80) mE_80->Fill (jet->energy());
00597     if (mE_3000) mE_3000->Fill (jet->energy());
00598     if (mP) mP->Fill (jet->p());
00599     if (mP_80) mP_80->Fill (jet->p());
00600     if (mP_3000) mP_3000->Fill (jet->p());
00601     if (mPt) mPt->Fill (jet->pt());
00602     if (mPt_80) mPt_80->Fill (jet->pt());
00603     if (mPt_3000) mPt_3000->Fill (jet->pt());
00604     if (mMass) mMass->Fill (jet->mass());
00605     if (mMass_80) mMass_80->Fill (jet->mass());
00606     if (mMass_3000) mMass_3000->Fill (jet->mass());
00607     if (mConstituents) mConstituents->Fill (jet->nConstituents());
00608     if (mConstituents_80) mConstituents_80->Fill (jet->nConstituents());
00609     if (jet == pfJets->begin ()) { // first jet
00610       if (mEtaFirst) mEtaFirst->Fill (jet->eta());
00611       if (mPhiFirst) mPhiFirst->Fill (jet->phi());
00612       if (mEFirst) mEFirst->Fill (jet->energy());
00613       if (mEFirst_80) mEFirst_80->Fill (jet->energy());
00614       if (mEFirst_3000) mEFirst_3000->Fill (jet->energy());
00615       if (mPtFirst) mPtFirst->Fill (jet->pt());
00616       if (mPtFirst_80) mPtFirst_80->Fill (jet->pt());
00617       if (mPtFirst_3000) mPtFirst_3000->Fill (jet->pt());
00618     }
00619     if (jetIndex == 0) {
00620       nJet++;
00621       p4tmp[0] = jet->p4();     
00622     }
00623     if (jetIndex == 1) {
00624       nJet++;
00625       p4tmp[1] = jet->p4();     
00626     }
00627 
00628     //    if (mMaxEInEmTowers) mMaxEInEmTowers->Fill (jet->maxEInEmTowers());
00629     //    if (mMaxEInHadTowers) mMaxEInHadTowers->Fill (jet->maxEInHadTowers());
00630     //    if (mHadEnergyInHO) mHadEnergyInHO->Fill (jet->hadEnergyInHO());
00631     //    if (mHadEnergyInHO_80)   mHadEnergyInHO_80->Fill (jet->hadEnergyInHO());
00632     //    if (mHadEnergyInHO_3000) mHadEnergyInHO_3000->Fill (jet->hadEnergyInHO());
00633     //    if (mHadEnergyInHB) mHadEnergyInHB->Fill (jet->hadEnergyInHB());
00634     //    if (mHadEnergyInHB_80)   mHadEnergyInHB_80->Fill (jet->hadEnergyInHB());
00635     //    if (mHadEnergyInHB_3000) mHadEnergyInHB_3000->Fill (jet->hadEnergyInHB());
00636     //    if (mHadEnergyInHF) mHadEnergyInHF->Fill (jet->hadEnergyInHF());
00637     //    if (mHadEnergyInHE) mHadEnergyInHE->Fill (jet->hadEnergyInHE());
00638     //    if (mHadEnergyInHE_80)   mHadEnergyInHE_80->Fill (jet->hadEnergyInHE());
00639     //    if (mHadEnergyInHE_3000) mHadEnergyInHE_3000->Fill (jet->hadEnergyInHE());
00640     //    if (mEmEnergyInEB) mEmEnergyInEB->Fill (jet->emEnergyInEB());
00641     //    if (mEmEnergyInEB_80)   mEmEnergyInEB_80->Fill (jet->emEnergyInEB());
00642     //    if (mEmEnergyInEB_3000) mEmEnergyInEB_3000->Fill (jet->emEnergyInEB());
00643     //    if (mEmEnergyInEE) mEmEnergyInEE->Fill (jet->emEnergyInEE());
00644     //    if (mEmEnergyInEE_80)   mEmEnergyInEE_80->Fill (jet->emEnergyInEE());
00645     //    if (mEmEnergyInEE_3000) mEmEnergyInEE_3000->Fill (jet->emEnergyInEE());
00646     //    if (mEmEnergyInHF) mEmEnergyInHF->Fill (jet->emEnergyInHF());
00647     //    if (mEnergyFractionHadronic) mEnergyFractionHadronic->Fill (jet->energyFractionHadronic());
00648     //    if (mEnergyFractionEm) mEnergyFractionEm->Fill (jet->emEnergyFraction());
00649     //
00650     //    if (mHFTotal)      mHFTotal->Fill (jet->hadEnergyInHF()+jet->emEnergyInHF());
00651     //    if (mHFTotal_80)   mHFTotal_80->Fill (jet->hadEnergyInHF()+jet->emEnergyInHF());
00652     //    if (mHFTotal_3000) mHFTotal_3000->Fill (jet->hadEnergyInHF()+jet->emEnergyInHF());
00653     //    if (mHFLong)       mHFLong->Fill (jet->hadEnergyInHF()*0.5+jet->emEnergyInHF());
00654     //    if (mHFLong_80)    mHFLong_80->Fill (jet->hadEnergyInHF()*0.5+jet->emEnergyInHF());
00655     //    if (mHFLong_3000)  mHFLong_3000->Fill (jet->hadEnergyInHF()*0.5+jet->emEnergyInHF());
00656     //    if (mHFShort)      mHFShort->Fill (jet->hadEnergyInHF()*0.5);
00657     //    if (mHFShort_80)   mHFShort_80->Fill (jet->hadEnergyInHF()*0.5);
00658     //    if (mHFShort_3000) mHFShort_3000->Fill (jet->hadEnergyInHF()*0.5);
00659 
00660     if (mChargedEmEnergy_80) mChargedEmEnergy_80->Fill (jet->chargedEmEnergy());
00661     if (mChargedHadronEnergy_80) mChargedHadronEnergy_80->Fill (jet->chargedHadronEnergy());
00662     if (mNeutralEmEnergy_80) mNeutralEmEnergy_80->Fill (jet->neutralEmEnergy());
00663     if (mNeutralHadronEnergy_80) mNeutralHadronEnergy_80->Fill (jet->neutralHadronEnergy());
00664 
00665     if (mChargedEmEnergy_3000) mChargedEmEnergy_3000->Fill (jet->chargedEmEnergy());
00666     if (mChargedHadronEnergy_3000) mChargedHadronEnergy_3000->Fill (jet->chargedHadronEnergy());
00667     if (mNeutralEmEnergy_3000) mNeutralEmEnergy_3000->Fill (jet->neutralEmEnergy());
00668     if (mNeutralHadronEnergy_3000) mNeutralHadronEnergy_3000->Fill (jet->neutralHadronEnergy());
00669 
00670     if (mChargedEmEnergyFraction) mChargedEmEnergyFraction->Fill (jet->chargedEmEnergyFraction());
00671     if (mChargedHadronEnergyFraction) mChargedHadronEnergyFraction->Fill (jet->chargedHadronEnergyFraction());
00672     if (mNeutralEmEnergyFraction) mNeutralEmEnergyFraction->Fill (jet->neutralEmEnergyFraction());
00673     if (mNeutralHadronEnergyFraction) mNeutralHadronEnergyFraction->Fill (jet->neutralHadronEnergyFraction());
00674 
00675     //    if (mN90) mN90->Fill (jet->n90());
00676     mJetEnergyProfile->Fill (jet->eta(), jet->phi(), jet->energy());
00677     //    mHadJetEnergyProfile->Fill (jet->eta(), jet->phi(), jet->hadEnergyInHO()+jet->hadEnergyInHB()+jet->hadEnergyInHF()+jet->hadEnergyInHE());
00678     //    mEMJetEnergyProfile->Fill (jet->eta(), jet->phi(), jet->emEnergyInEB()+jet->emEnergyInEE()+jet->emEnergyInHF());
00679   }
00680 
00681 
00682 
00683 
00684   if (mNJetsEtaC) mNJetsEtaC->Fill( nJetC );
00685   if (mNJetsEtaF) mNJetsEtaF->Fill( nJetF );
00686 
00687   if (nJet == 2) {
00688     if (mMjj) mMjj->Fill( (p4tmp[0]+p4tmp[1]).mass() );
00689     if (mMjj_3000) mMjj_3000->Fill( (p4tmp[0]+p4tmp[1]).mass() );
00690   }
00691 
00692   // Count Jets above Pt cut
00693   for (int istep = 0; istep < 100; ++istep) {
00694     int     njet = 0;
00695     float ptStep = (istep * (200./100.));
00696 
00697     for ( PFJetCollection::const_iterator cal = pfJets->begin(); cal != pfJets->end(); ++ cal ) {
00698       if ( cal->pt() > ptStep ) njet++;
00699     }
00700     mNJets1->Fill( ptStep, njet );
00701   }
00702 
00703   for (int istep = 0; istep < 100; ++istep) {
00704     int     njet = 0;
00705     float ptStep = (istep * (4000./100.));
00706     for ( PFJetCollection::const_iterator cal = pfJets->begin(); cal != pfJets->end(); ++ cal ) {
00707       if ( cal->pt() > ptStep ) njet++;
00708     }
00709     mNJets2->Fill( ptStep, njet );
00710   }
00711 
00712   // Gen jet analysis
00713   Handle<GenJetCollection> genJets;
00714   mEvent.getByLabel(mInputGenCollection, genJets);
00715   if (!genJets.isValid()) return;
00716   GenJetCollection::const_iterator gjet = genJets->begin ();
00717   int gjetIndex = 0;
00718   for (; gjet != genJets->end (); gjet++, gjetIndex++) {
00719     if (mGenEta) mGenEta->Fill (gjet->eta());
00720     if (mGenPhi) mGenPhi->Fill (gjet->phi());
00721     if (mGenPt) mGenPt->Fill (gjet->pt());
00722     if (mGenPt_80) mGenPt_80->Fill (gjet->pt());
00723     if (mGenPt_3000) mGenPt_3000->Fill (gjet->pt());
00724     if (gjet == genJets->begin ()) { // first jet
00725       if (mGenEtaFirst) mGenEtaFirst->Fill (gjet->eta());
00726       if (mGenPhiFirst) mGenPhiFirst->Fill (gjet->phi());
00727     }
00728   }
00729 
00730 
00731   // now match PFJets to GenJets
00732   JetMatchingTools jetMatching (mEvent);
00733   if (!(mInputGenCollection.label().empty())) {
00734     //    Handle<GenJetCollection> genJets;
00735     //    mEvent.getByLabel(mInputGenCollection, genJets);
00736 
00737     std::vector <std::vector <const reco::GenParticle*> > genJetConstituents (genJets->size());
00738     std::vector <std::vector <const reco::GenParticle*> > pfJetConstituents (pfJets->size());
00739     /*
00740     if (mRThreshold > 0) { 
00741     }
00742     else {
00743       for (unsigned iGenJet = 0; iGenJet < genJets->size(); ++iGenJet) {
00744         genJetConstituents [iGenJet] = jetMatching.getGenParticles ((*genJets) [iGenJet]);
00745       }
00746       
00747       for (unsigned iPFJet = 0; iPFJet < pfJets->size(); ++iPFJet) {
00748         pfJetConstituents [iPFJet] = jetMatching.getGenParticles ((*pfJets) [iPFJet], false);
00749       }
00750     }
00751     */
00752     
00753     for (unsigned iGenJet = 0; iGenJet < genJets->size(); ++iGenJet) {               //****************************************************************
00754     //for (unsigned iGenJet = 0; iGenJet < 1; ++iGenJet) {                           // only FIRST Jet !!!!
00755       const GenJet& genJet = (*genJets) [iGenJet];
00756       double genJetPt = genJet.pt();
00757 
00758       //std::cout << iGenJet <<". Genjet: pT = " << genJetPt << "GeV" << std::endl;  //  *****************************************************
00759 
00760       if (fabs(genJet.eta()) > 5.) continue; // out of detector 
00761       if (genJetPt < mMatchGenPtThreshold) continue; // no low momentum 
00762       double logPtGen = log10 (genJetPt);
00763       mAllGenJetsPt->Fill (logPtGen);
00764       mAllGenJetsEta->Fill (logPtGen, genJet.eta());
00765       if (pfJets->size() <= 0) continue; // no PFJets - nothing to match
00766       if (mRThreshold > 0) {
00767         unsigned iPFJetBest = 0;
00768         double deltaRBest = 999.;
00769         for (unsigned iPFJet = 0; iPFJet < pfJets->size(); ++iPFJet) {
00770           double dR = deltaR (genJet.eta(), genJet.phi(), (*pfJets) [iPFJet].eta(), (*pfJets) [iPFJet].phi());
00771           if (deltaRBest < mRThreshold && dR < mRThreshold && genJet.pt() > 5.) {
00772             /*
00773             std::cout << "Yet another matched jet for GenJet pt=" << genJet.pt()
00774                       << " previous PFJet pt/dr: " << (*pfJets) [iPFJetBest].pt() << '/' << deltaRBest
00775                       << " new PFJet pt/dr: " << (*pfJets) [iPFJet].pt() << '/' << dR
00776                       << std::endl;
00777             */
00778           }
00779           if (dR < deltaRBest) {
00780             iPFJetBest = iPFJet;
00781             deltaRBest = dR;
00782           }
00783         }
00784         if (mTurnOnEverything.compare("yes")==0) {
00785           mRMatch->Fill (logPtGen, genJet.eta(), deltaRBest);
00786         }
00787         if (deltaRBest < mRThreshold) { // Matched
00788           fillMatchHists (genJet, (*pfJets) [iPFJetBest]);
00789         }
00790       }
00791       /*
00792       else {
00793         unsigned iPFJetBest = 0;
00794         double energyFractionBest = 0.;
00795         for (unsigned iPFJet = 0; iPFJet < pfJets->size(); ++iPFJet) {
00796           double energyFraction = jetMatching.overlapEnergyFraction (genJetConstituents [iGenJet], 
00797                                                                      pfJetConstituents [iPFJet]);
00798           if (energyFraction > energyFractionBest) {
00799             iPFJetBest = iPFJet;
00800             energyFractionBest = energyFraction;
00801           }
00802         }
00803         if (mTurnOnEverything.compare("yes")==0) {
00804           mGenJetMatchEnergyFraction->Fill (logPtGen, genJet.eta(), energyFractionBest);
00805         }
00806         if (energyFractionBest > mGenEnergyFractionThreshold) { // good enough
00807           double reverseEnergyFraction = jetMatching.overlapEnergyFraction (pfJetConstituents [iPFJetBest], 
00808                                                                             genJetConstituents [iGenJet]);
00809           if (mTurnOnEverything.compare("yes")==0) {
00810             mReverseMatchEnergyFraction->Fill (logPtGen, genJet.eta(), reverseEnergyFraction);
00811           }
00812           if (reverseEnergyFraction > mReverseEnergyFractionThreshold) { // Matched
00813             fillMatchHists (genJet, (*pfJets) [iPFJetBest]);
00814           }
00815         }
00816       }
00817       */
00818     }
00819   }
00820 }
00821 
00822 
00823 void PFJetTester::fillMatchHists (const reco::GenJet& fGenJet, const reco::PFJet& fPFJet) {
00824   double logPtGen = log10 (fGenJet.pt());
00825   double PtGen = fGenJet.pt();
00826   double PtPF = fPFJet.pt();
00827   mMatchedGenJetsPt->Fill (logPtGen);
00828   mMatchedGenJetsEta->Fill (logPtGen, fGenJet.eta());
00829 
00830   double PtThreshold = 10.;
00831 
00832   if (mTurnOnEverything.compare("yes")==0) {
00833     mDeltaEta->Fill (logPtGen, fGenJet.eta(), fPFJet.eta()-fGenJet.eta());
00834     mDeltaPhi->Fill (logPtGen, fGenJet.eta(), fPFJet.phi()-fGenJet.phi());
00835     mEScale->Fill (logPtGen, fGenJet.eta(), fPFJet.energy()/fGenJet.energy());
00836     mlinEScale->Fill (fGenJet.pt(), fGenJet.eta(), fPFJet.energy()/fGenJet.energy());
00837     mDeltaE->Fill (logPtGen, fGenJet.eta(), fPFJet.energy()-fGenJet.energy());
00838 
00839     mEScaleFineBin->Fill (logPtGen, fGenJet.eta(), fPFJet.energy()/fGenJet.energy());
00840   
00841     if (fGenJet.pt()>PtThreshold) {
00842       mEScale_pt10->Fill (logPtGen, fGenJet.eta(), fPFJet.energy()/fGenJet.energy());
00843 
00844     }
00845 
00846   }
00847   if (fPFJet.pt() > PtThreshold) {
00848     mDelEta->Fill (fGenJet.eta()-fPFJet.eta());
00849     mDelPhi->Fill (fGenJet.phi()-fPFJet.phi());
00850     mDelPt->Fill  ((fGenJet.pt()-fPFJet.pt())/fGenJet.pt());
00851   }
00852 
00853   if (fabs(fGenJet.eta())<1.3) {
00854 
00855     mpTScaleB_s->Fill (log10(PtGen), PtPF/PtGen);
00856     mpTScaleB_d->Fill (log10(PtGen), PtPF/PtGen);
00857     if (PtGen>60.0 && PtGen<120.0) {
00858       mpTScale1DB_60_120->Fill (fPFJet.pt()/fGenJet.pt());
00859     }
00860     if (PtGen>200.0 && PtGen<300.0) {
00861       mpTScale1DB_200_300->Fill (fPFJet.pt()/fGenJet.pt());
00862     }
00863     if (PtGen>600.0 && PtGen<900.0) {
00864       mpTScale1DB_600_900->Fill (fPFJet.pt()/fGenJet.pt());
00865     }
00866     if (PtGen>2700.0 && PtGen<3500.0) {
00867       mpTScale1DB_2700_3500->Fill (fPFJet.pt()/fGenJet.pt());
00868     }
00869   }
00870 
00871   if (fabs(fGenJet.eta())>1.3 && fabs(fGenJet.eta())<3.0) {
00872 
00873     mpTScaleE_s->Fill (log10(PtGen), PtPF/PtGen);
00874     mpTScaleE_d->Fill (log10(PtGen), PtPF/PtGen);
00875     if (PtGen>60.0 && PtGen<120.0) {
00876       mpTScale1DE_60_120->Fill (fPFJet.pt()/fGenJet.pt());
00877     }
00878     if (PtGen>200.0 && PtGen<300.0) {
00879       mpTScale1DE_200_300->Fill (fPFJet.pt()/fGenJet.pt());
00880     }
00881     if (PtGen>600.0 && PtGen<900.0) {
00882       mpTScale1DE_600_900->Fill (fPFJet.pt()/fGenJet.pt());
00883     }
00884     if (PtGen>2700.0 && PtGen<3500.0) {
00885       mpTScale1DE_2700_3500->Fill (fPFJet.pt()/fGenJet.pt());
00886     }
00887   }
00888 
00889   if (fabs(fGenJet.eta())>3.0 && fabs(fGenJet.eta())<5.0) {
00890 
00891     mpTScaleF_s->Fill (log10(PtGen), PtPF/PtGen);
00892     mpTScaleF_d->Fill (log10(PtGen), PtPF/PtGen);
00893     if (PtGen>60.0 && PtGen<120.0) {
00894       mpTScale1DF_60_120->Fill (fPFJet.pt()/fGenJet.pt());
00895     }
00896     if (PtGen>200.0 && PtGen<300.0) {
00897       mpTScale1DF_200_300->Fill (fPFJet.pt()/fGenJet.pt());
00898     }
00899     if (PtGen>600.0 && PtGen<900.0) {
00900       mpTScale1DF_600_900->Fill (fPFJet.pt()/fGenJet.pt());
00901     }
00902     if (PtGen>2700.0 && PtGen<3500.0) {
00903       mpTScale1DF_2700_3500->Fill (fPFJet.pt()/fGenJet.pt());
00904     }
00905   }
00906 
00907   if (fGenJet.pt()>60.0 && fGenJet.pt()<120.0) {
00908     mpTScale_60_120_s->Fill (fGenJet.eta(),fPFJet.pt()/fGenJet.pt());
00909     mpTScale_60_120_d->Fill (fGenJet.eta(),fPFJet.pt()/fGenJet.pt());
00910     mpTScale1D_60_120->Fill (fPFJet.pt()/fGenJet.pt());
00911   }
00912 
00913   if (fGenJet.pt()>200.0 && fGenJet.pt()<300.0) {
00914     mpTScale_200_300_s->Fill (fGenJet.eta(),fPFJet.pt()/fGenJet.pt());
00915     mpTScale_200_300_d->Fill (fGenJet.eta(),fPFJet.pt()/fGenJet.pt());
00916     mpTScale1D_200_300->Fill (fPFJet.pt()/fGenJet.pt());
00917   }
00918 
00919   if (fGenJet.pt()>600.0 && fGenJet.pt()<900.0) {
00920     mpTScale_600_900_s->Fill (fGenJet.eta(),fPFJet.pt()/fGenJet.pt());
00921     mpTScale_600_900_d->Fill (fGenJet.eta(),fPFJet.pt()/fGenJet.pt());
00922     mpTScale1D_600_900->Fill (fPFJet.pt()/fGenJet.pt());
00923   }
00924 
00925   if (fGenJet.pt()>2700.0 && fGenJet.pt()<3500.0) {
00926     mpTScale_2700_3500_s->Fill (fGenJet.eta(),fPFJet.pt()/fGenJet.pt());
00927     mpTScale_2700_3500_d->Fill (fGenJet.eta(),fPFJet.pt()/fGenJet.pt());
00928     mpTScale1D_2700_3500->Fill (fPFJet.pt()/fGenJet.pt());
00929   }
00930 
00931 
00932 
00933 }