Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:26:55

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