Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 09:46:14

0001 #include "FastSimulation/Validation/test/JetComparison.h"
0002 
0003 //#define EDM_ML_DEBUG
0004 
0005 JetComparison::JetComparison(edm::ParameterSet const& conf)
0006     : jetToken_(consumes<reco::CaloJetCollection>(edm::InputTag("iterativeCone5CaloJets"))),
0007       genjetToken_(consumes<reco::GenJetCollection>(edm::InputTag("ak4GenJets"))),
0008       fMinEnergy(conf.getParameter<double>("MinEnergy")) {
0009   usesResource(TFileService::kSharedResource);
0010 
0011   nEvent = 0;
0012 
0013   edm::Service<TFileService> tf;
0014   meNumberJet = tf->make<TH1F>("number_of_jets", "number_of_jets", 500, 0., 500.);
0015   meEtJet = tf->make<TH1F>("jet_et", "jet_et", 1000, 0., 1000.);
0016   meEtGen = tf->make<TH1F>("gen_et", "gen_et", 1000, 0., 1000.);
0017   meEtJetMatched = tf->make<TH1F>("jet_et_two_matched", "jet_et_two_matched", 1000, 0., 1000.);
0018   meEtGenMatched = tf->make<TH1F>("gen_et_two_matched", "gen_et_two_matched", 1000, 0., 1000.);
0019   meEtaJet = tf->make<TH1F>("jet_eta", "jet_eta", 500, 0., 5.);
0020   meEtaGen = tf->make<TH1F>("gen_eta", "gen_eta", 500, 0., 5.);
0021   meRatio = tf->make<TH2F>("ratio_jet_gen_vs_eta", "ratio_jet_gen_vs_eta", 50, 0., 5., 20000, 0., 5.);
0022   meHadronicFrac_vs_eta =
0023       tf->make<TH2F>("hadronic_fraction_vs_eta", "hadronic_fraction_vs_eta", 50, 0., 5., 1000, 0., 1.);
0024   meNTowers60_vs_eta = tf->make<TH2F>("number_of_towers_with_60_percents_energy_vs_eta",
0025                                       "number_of_towers_with_60_percents_energy_vs_eta",
0026                                       50,
0027                                       0.,
0028                                       5.,
0029                                       100,
0030                                       -0.5,
0031                                       99.5);
0032   meNTowers90_vs_eta = tf->make<TH2F>("number_of_towers_with_90_percents_energy_vs_eta",
0033                                       "number_of_towers_with_90_percents_energy_vs_eta",
0034                                       50,
0035                                       0.,
0036                                       5.,
0037                                       100,
0038                                       -0.5,
0039                                       99.5);
0040   meDistR = tf->make<TH1F>("distance_R_between_matched_jets", "distance_R_between_matched_jets", 300, 0., 0.3);
0041   meDistR_vs_eta = tf->make<TH2F>(
0042       "distance_R_between_matched_jets_vs_eta", "stance_R_between_matched_jets_vs_eta", 50, 0., 5., 300, 0., 0.3);
0043   meNTowers_vs_eta =
0044       tf->make<TH2F>("number_of_calotowers_vs_eta", "number_of_calotowers_vs_eta", 50, 0., 5., 51, -0.5, 50.5);
0045 }
0046 
0047 void JetComparison::endJob() {
0048   meNumberJet->Scale(1. / nEvent);
0049   Float_t x[50];
0050   Float_t y[50];
0051   Float_t y_er[50];
0052   Float_t x_er[50];
0053   Float_t mean;
0054   Float_t sigm;
0055 
0056   for (Int_t i = 1; i < 51; i++) {
0057     mean = (meRatio->ProjectionY("", i, i, "")->GetMean());
0058     sigm = (meRatio->ProjectionY("", i, i, "")->GetRMS());
0059     edm::LogVerbatim("JetComparison") << " Mean " << mean;
0060     x[i - 1] = (i - 0.5) * 0.1;
0061     y[i - 1] = mean;
0062     y_er[i - 1] = sigm;
0063     x_er[i - 1] = 0.;
0064   }
0065 
0066   edm::Service<TFileService> tf;
0067   gr = tf->make<TGraphErrors>(50, x, y, x_er, y_er);
0068   gr->SetMarkerStyle(21);
0069   gr->SetMarkerColor(44);
0070 }
0071 
0072 void JetComparison::beginRun(edm::Run const& run, edm::EventSetup const& iSetup) {}
0073 
0074 void JetComparison::analyze(edm::Event const& event, edm::EventSetup const& c) {
0075   nEvent++;
0076 
0077   const edm::Handle<reco::CaloJetCollection>& jets = event.getHandle(jetToken_);
0078   meNumberJet->Fill(jets->size());
0079 
0080   for (unsigned int ijet = 0; ijet < jets->size(); ijet++) {
0081 #ifdef EDM_ML_DEBUG
0082     edm::LogVerbatim("JetComparison") << " JETS energy = " << (*jets)[ijet].et();
0083 #endif
0084     double etJet = (*jets)[ijet].et();
0085     double etaJet = (*jets)[ijet].eta();
0086 
0087     meEtJet->Fill(etJet);
0088     meEtaJet->Fill(etaJet);
0089   }
0090 
0091   const edm::Handle<reco::GenJetCollection>& jetsgen = event.getHandle(genjetToken_);
0092   for (unsigned int igen = 0; igen < jetsgen->size(); igen++) {
0093 #ifdef EDM_ML_DEBUG
0094     edm::LogVerbatim("JetComparison") << " GENS energy = " << (*jetsgen)[igen].et();
0095 #endif
0096     double etGen = (*jetsgen)[igen].et();
0097     double etaGen = (*jetsgen)[igen].eta();
0098     meEtGen->Fill(etGen);
0099     meEtaGen->Fill(etaGen);
0100   }
0101 
0102   for (unsigned int ijet = 0; ijet < jets->size(); ijet++) {
0103     double etaJet = (*jets)[ijet].eta();
0104     double phiJet = (*jets)[ijet].phi();
0105     double etJet = (*jets)[ijet].et();
0106     double fracHadr = (*jets)[ijet].energyFractionHadronic();
0107     int n60 = (*jets)[ijet].n60();
0108     int n90 = (*jets)[ijet].n90();
0109     int nCaloTower = (*jets)[ijet].nConstituents();
0110     double etaGen;
0111     double phiGen;
0112     double etGen;
0113 
0114     unsigned int isize = jetsgen->size();
0115     if (isize > 2)
0116       isize = 2;
0117     for (unsigned int igen = 0; igen < isize; igen++) {
0118       etaGen = (*jetsgen)[igen].eta();
0119       phiGen = (*jetsgen)[igen].phi();
0120       etGen = (*jetsgen)[igen].et();
0121 
0122       double rDist = std::sqrt(deltaR2(etaJet, phiJet, etaGen, phiGen));
0123 
0124       if (rDist < 0.3 && etGen > fMinEnergy) {
0125         meNTowers_vs_eta->Fill(fabs(etaGen), nCaloTower);
0126         meDistR->Fill(rDist);
0127         meDistR_vs_eta->Fill(fabs(etaGen), rDist);
0128         edm::LogVerbatim("JetComparison")
0129             << "  JETS energy = " << (*jets)[ijet].et() << " ETA = " << etaJet << " PHI = " << phiJet;
0130         edm::LogVerbatim("JetComparison")
0131             << "  GENS energy = " << (*jetsgen)[igen].et() << " ETA = " << etaGen << " PHI = " << phiGen;
0132         edm::LogVerbatim("JetComparison") << "Jets matched ";
0133         meRatio->Fill(fabs(etaGen), etJet / etGen);
0134         meEtJetMatched->Fill(etJet);
0135         meEtGenMatched->Fill(etGen);
0136         meNTowers60_vs_eta->Fill(fabs(etaJet), n60);
0137         meNTowers90_vs_eta->Fill(fabs(etaJet), n90);
0138         meHadronicFrac_vs_eta->Fill(fabs(etaJet), fracHadr);
0139       }
0140     }
0141   }
0142 }
0143 
0144 double JetComparison::deltaR2(double eta0, double phi0, double eta, double phi) {
0145   double dphi = phi - phi0;
0146   if (dphi > M_PI)
0147     dphi -= 2 * M_PI;
0148   else if (dphi <= -M_PI)
0149     dphi += 2 * M_PI;
0150   return dphi * dphi + (eta - eta0) * (eta - eta0);
0151 }
0152 
0153 #include "FWCore/PluginManager/interface/ModuleDef.h"
0154 #include "FWCore/Framework/interface/MakerMacros.h"
0155 
0156 DEFINE_FWK_MODULE(JetComparison);