File indexing completed on 2024-06-22 02:24:05
0001 #include "RecoParticleFlow/PFClusterProducer/test/PFClusterComparator.h"
0002 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0003
0004 #include "FWCore/Framework/interface/ESHandle.h"
0005
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/Utilities/interface/Exception.h"
0008 #include "FWCore/Framework/interface/EventSetup.h"
0009
0010 #include <iomanip>
0011
0012 using namespace std;
0013 using namespace edm;
0014 using namespace reco;
0015
0016 PFClusterComparator::PFClusterComparator(const edm::ParameterSet& iConfig)
0017 : inputTokenPFClusters_(consumes<reco::PFClusterCollection>(iConfig.getParameter<InputTag>("PFClusters"))),
0018 inputTokenPFClustersCompare_(
0019 consumes<reco::PFClusterCollection>(iConfig.getParameter<InputTag>("PFClustersCompare"))),
0020 verbose_(iConfig.getUntrackedParameter<bool>("verbose", false)),
0021 printBlocks_(iConfig.getUntrackedParameter<bool>("printBlocks", false)) {
0022 usesResource("TFileService");
0023
0024 edm::Service<TFileService> fs;
0025 log10E_old = fs->make<TH1F>("log10E_old", "log10(E cluster)", 500, -5, 5);
0026 log10E_new = fs->make<TH1F>("log10E_new", "log10(E cluster)", 500, -5, 5);
0027 deltaEnergy = fs->make<TH1F>("delta_energy", "E_{old} - E_{new}", 5000, -5, 5);
0028
0029 posX_old = fs->make<TH1F>("posX_old", "log10(E cluster)", 50000, 0, 500);
0030 posX_new = fs->make<TH1F>("posX_new", "log10(E cluster)", 50000, 0, 500);
0031 deltaX = fs->make<TH1F>("delta_X", "X_{old} - X_{new}", 5000, -5, 5);
0032
0033 posY_old = fs->make<TH1F>("posY_old", "log10(E cluster)", 50000, 0, 500);
0034 posY_new = fs->make<TH1F>("posY_new", "log10(E cluster)", 50000, 0, 500);
0035 deltaY = fs->make<TH1F>("delta_Y", "Y_{old} - Y_{new}", 5000, -5, 5);
0036
0037 posZ_old = fs->make<TH1F>("posZ_old", "log10(E cluster)", 50000, 0, 500);
0038 posZ_new = fs->make<TH1F>("posZ_new", "log10(E cluster)", 50000, 0, 500);
0039 deltaZ = fs->make<TH1F>("delta_Z", "Z_{old} - Z_{new}", 5000, -5, 5);
0040
0041 LogDebug("PFClusterComparator") << " input collection : " << iConfig.getParameter<InputTag>("PFClusters");
0042 }
0043
0044 void PFClusterComparator::analyze(const Event& iEvent, const EventSetup& iSetup) {
0045 std::map<unsigned, unsigned> detId_count;
0046
0047
0048
0049 Handle<PFClusterCollection> pfClusters;
0050 fetchCandidateCollection(pfClusters, inputTokenPFClusters_, iEvent);
0051
0052 Handle<PFClusterCollection> pfClustersCompare;
0053 fetchCandidateCollection(pfClustersCompare, inputTokenPFClustersCompare_, iEvent);
0054
0055
0056
0057 std::cout << "There are " << pfClusters->size() << " PFClusters in the original cluster collection." << std::endl;
0058 std::cout << "There are " << pfClustersCompare->size() << " PFClusters in the new cluster collection." << std::endl;
0059
0060 std::cout << std::flush << "---- COMPARING OLD TO NEW ----" << std::endl << std::flush;
0061
0062 for (auto const& cluster : *pfClusters) {
0063 detId_count[cluster.seed().rawId()] += 1;
0064 log10E_old->Fill(std::log10(cluster.energy()));
0065 posX_old->Fill(std::abs(cluster.position().x()));
0066 posY_old->Fill(std::abs(cluster.position().y()));
0067 posZ_old->Fill(std::abs(cluster.position().z()));
0068 bool foundmatch = false;
0069 for (auto const& clustercomp : *pfClustersCompare) {
0070 if (cluster.seed().rawId() == clustercomp.seed().rawId()) {
0071 foundmatch = true;
0072 const double denergy = std::abs(cluster.energy() - clustercomp.energy());
0073 const double dcenergy = std::abs(cluster.correctedEnergy() - clustercomp.correctedEnergy());
0074 const double dx = std::abs(cluster.position().x() - clustercomp.position().x());
0075 const double dy = std::abs(cluster.position().y() - clustercomp.position().y());
0076 const double dz = std::abs(cluster.position().z() - clustercomp.position().z());
0077 deltaEnergy->Fill((cluster.energy() - clustercomp.energy()) / cluster.energy());
0078 deltaX->Fill((cluster.position().x() - clustercomp.position().x()) / cluster.position().x());
0079 deltaY->Fill((cluster.position().y() - clustercomp.position().y()) / cluster.position().y());
0080 deltaZ->Fill((cluster.position().z() - clustercomp.position().z()) / cluster.position().z());
0081
0082 if (denergy / std::abs(cluster.energy()) > 1e-5) {
0083 std::cout << " " << cluster.seed() << " Energies different by larger than tolerance! "
0084 << "( " << denergy << " )"
0085 << "[ " << detId_count[cluster.seed().rawId()] << " ]"
0086 << " Old: " << std::setprecision(7) << cluster.energy() << " GeV , New: " << clustercomp.energy()
0087 << " GeV" << std::endl;
0088 }
0089 if (dcenergy / std::abs(cluster.correctedEnergy()) > 1e-5) {
0090 std::cout << " " << cluster.seed() << " Corrected energies different by larger than tolerance! "
0091 << "( " << dcenergy << " )"
0092 << "[ " << detId_count[cluster.seed().rawId()] << " ]"
0093 << " Old: " << std::setprecision(7) << cluster.correctedEnergy()
0094 << " GeV , New: " << clustercomp.correctedEnergy() << " GeV" << std::endl;
0095 }
0096 std::cout << std::flush;
0097 if (dx / std::abs(cluster.position().x()) > 1e-5) {
0098 std::cout << "***" << cluster.seed() << " X's different by larger than tolerance! "
0099 << "( " << dx << " )"
0100 << "[ " << detId_count[cluster.seed().rawId()] << " ]"
0101 << " Old: " << std::setprecision(7) << cluster.position().x()
0102 << " , New: " << clustercomp.position().x() << std::endl;
0103 }
0104 std::cout << std::flush;
0105 if (dy / std::abs(cluster.position().y()) > 1e-5) {
0106 std::cout << "---" << cluster.seed() << " Y's different by larger than tolerance! "
0107 << "( " << dy << " )"
0108 << "[ " << detId_count[cluster.seed().rawId()] << " ]"
0109 << " Old: " << std::setprecision(7) << cluster.position().y()
0110 << " , New: " << clustercomp.position().y() << std::endl;
0111 }
0112 std::cout << std::flush;
0113 if (dz / std::abs(cluster.position().z()) > 1e-5) {
0114 std::cout << "+++" << cluster.seed() << " Z's different by larger than tolerance! "
0115 << "( " << dz << " )"
0116 << "[ " << detId_count[cluster.seed().rawId()] << " ]"
0117 << " Old: " << std::setprecision(7) << cluster.position().z()
0118 << " , New: " << clustercomp.position().z() << std::endl;
0119 }
0120 std::cout << std::flush;
0121 }
0122 }
0123 if (!foundmatch) {
0124 std::cout << "Seed in old clusters and not new: " << cluster << std::endl;
0125 }
0126 }
0127
0128 std::cout << std::flush << "---- COMPARING NEW TO OLD ----" << std::endl << std::flush;
0129
0130 for (auto const& cluster : *pfClustersCompare) {
0131 log10E_new->Fill(std::log10(cluster.energy()));
0132 posX_new->Fill(std::abs(cluster.position().x()));
0133 posY_new->Fill(std::abs(cluster.position().y()));
0134 posZ_new->Fill(std::abs(cluster.position().z()));
0135 bool foundmatch = false;
0136 for (auto const& clustercomp : *pfClusters) {
0137 if (cluster.seed() == clustercomp.seed()) {
0138 foundmatch = true;
0139 }
0140 }
0141 if (!foundmatch) {
0142 std::cout << "Seed in new clusters and not old: " << cluster << std::endl;
0143 }
0144 }
0145 std::cout << std::flush;
0146 }
0147
0148 void PFClusterComparator::fetchCandidateCollection(Handle<reco::PFClusterCollection>& c,
0149 const edm::EDGetTokenT<reco::PFClusterCollection>& token,
0150 const Event& iEvent) const {
0151 c = iEvent.getHandle(token);
0152
0153 if (!c.isValid()) {
0154 ostringstream err;
0155 err << " cannot get PFClusters " << endl;
0156 LogError("PFClusters") << err.str();
0157 throw cms::Exception("MissingProduct", err.str());
0158 }
0159 }
0160
0161 DEFINE_FWK_MODULE(PFClusterComparator);
