Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-03-26 02:43:54

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