EcalClusters/src/EnergyScaleAnalyzer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    EnergyScaleAnalyzer
0004 // Class:      EnergyScaleAnalyzer
0005 //
0006 /**\class EnergyScaleAnalyzer EnergyScaleAnalyzer.cc
0007  Validation/EcalClusters/src/EnergyScaleAnalyzer.cc
0008
0009  Description: <one line class summary>
0010
0011  Implementation:
0012      <Notes on implementation>
0013 */
0014 // Original Author:  Keti Kaadze
0015 //         Created:  Thu Jun 21 08:59:42 CDT 2007
0016 //
0017
0018 #include "Validation/EcalClusters/interface/EnergyScaleAnalyzer.h"
0019
0020 // Framework
0021 #include "DataFormats/Common/interface/Handle.h"
0022 #include "FWCore/Framework/interface/ESHandle.h"
0023 #include "FWCore/Framework/interface/Event.h"
0024 #include "FWCore/Framework/interface/EventSetup.h"
0025 #include "FWCore/Framework/interface/Frameworkfwd.h"
0026 #include "FWCore/Framework/interface/MakerMacros.h"
0027 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029 #include "FWCore/Utilities/interface/Exception.h"
0030
0031 // Geometry
0032 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0033 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0034 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0035 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0036 #include "Geometry/CaloTopology/interface/EcalBarrelHardcodedTopology.h"
0037 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0038
0039 #include "TFile.h"
0040 // Reconstruction classes
0041 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0042 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0043 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0044 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0045 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
0046 #include "DataFormats/EgammaReco/interface/PreshowerCluster.h"
0047 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0048 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0049
0050 #include "DataFormats/EgammaReco/interface/BasicClusterShapeAssociation.h"
0051 #include "DataFormats/Math/interface/LorentzVector.h"
0052 #include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
0053
0054 //// Class header file
0055 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
0056 #include "DataFormats/EgammaReco/interface/ClusterShape.h"
0057 #include "DataFormats/EgammaReco/interface/ClusterShapeFwd.h"
0058 #include "RecoEcal/EgammaCoreTools/interface/ClusterShapeAlgo.h"
0059 #include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
0060
0061 #include "TString.h"
0062 #include "Validation/EcalClusters/interface/AnglesUtil.h"
0063 #include <fstream>
0064 #include <iomanip>
0065 #include <ios>
0066 #include <iostream>
0067 #include <map>
0068
0069 //========================================================================
0070 EnergyScaleAnalyzer::EnergyScaleAnalyzer(const edm::ParameterSet &ps)
0071 //========================================================================
0072 {
0073   hepMCLabel_ = consumes<edm::HepMCProduct>(ps.getParameter<std::string>("hepMCLabel"));
0074   hybridSuperClusters_token = consumes<reco::SuperClusterCollection>(
0075       ps.getUntrackedParameter<edm::InputTag>("hybridSuperClusters", edm::InputTag("hybridSuperClusters")));
0076   dynamicHybridSuperClusters_token = consumes<reco::SuperClusterCollection>(ps.getUntrackedParameter<edm::InputTag>(
0077       "dynamicHybridSuperClusters", edm::InputTag("dynamicHybridSuperClusters")));
0078   correctedHybridSuperClusters_token = consumes<reco::SuperClusterCollection>(ps.getUntrackedParameter<edm::InputTag>(
0079       "correctedHybridSuperClusters", edm::InputTag("correctedHybridSuperClusters")));
0080   correctedDynamicHybridSuperClusters_token =
0081       consumes<reco::SuperClusterCollection>(ps.getUntrackedParameter<edm::InputTag>(
0082           "correctedDynamicHybridSuperClusters", edm::InputTag("correctedDynamicHybridSuperClusters")));
0083   correctedFixedMatrixSuperClustersWithPreshower_token = consumes<reco::SuperClusterCollection>(
0084       ps.getUntrackedParameter<edm::InputTag>("correctedFixedMatrixSuperClustersWithPreshower",
0085                                               edm::InputTag("correctedFixedMatrixSuperClustersWithPreshower")));
0086   fixedMatrixSuperClustersWithPreshower_token =
0087       consumes<reco::SuperClusterCollection>(ps.getUntrackedParameter<edm::InputTag>(
0088           "fixedMatrixSuperClustersWithPreshower", edm::InputTag("fixedMatrixSuperClustersWithPreshower")));
0089
0090   outputFile_ = ps.getParameter<std::string>("outputFile");
0091   rootFile_ = TFile::Open(outputFile_.c_str(),
0092                           "RECREATE");  // open output file to store histograms
0093
0094   evtN = 0;
0095 }
0096
0097 //========================================================================
0098 EnergyScaleAnalyzer::~EnergyScaleAnalyzer()
0099 //========================================================================
0100 {
0101   delete rootFile_;
0102 }
0103
0104 //========================================================================
0105 void EnergyScaleAnalyzer::beginJob() {
0106   //========================================================================
0107
0108   mytree_ = new TTree("energyScale", "");
0109   TString treeVariables =
0110       "mc_npar/I:parID:mc_sep/"
0111       "F:mc_e:mc_et:mc_phi:mc_eta:mc_theta:";  // MC information
0112   treeVariables += "em_dR/F:";                 // MC <-> EM matching information
0113   treeVariables +=
0114       "em_isInCrack/I:em_scType:em_e/F:em_et:em_phi:em_eta:em_theta:em_nCell/"
0115       "I:em_nBC:";                                       // EM SC info
0116   treeVariables += "em_pet/F:em_pe:em_peta:em_ptheta:";  // EM SC physics (eta corrected
0117                                                          // information)
0118
0119   treeVariables += "emCorr_e/F:emCorr_et:emCorr_eta:emCorr_phi:emCorr_theta:";  // CMSSW
0120                                                                                 // standard
0121                                                                                 // corrections
0122   treeVariables += "emCorr_pet/F:emCorr_peta:emCorr_ptheta:";                   // CMSSW standard physics
0123
0124   treeVariables += "em_pw/F:em_ew:em_br";  // EM widths pw -- phiWidth, ew --
0125                                            // etaWidth, ratios of pw/ew
0126
0127   mytree_->Branch("energyScale", &(tree_.mc_npar), treeVariables);
0128 }
0129
0130 //========================================================================
0131 void EnergyScaleAnalyzer::analyze(const edm::Event &evt, const edm::EventSetup &es) {
0132   using namespace edm;  // needed for all fwk related classes
0133   using namespace std;
0134
0135   //  std::cout << "Proceccing event # " << ++evtN << std::endl;
0136
0137   // Get containers for MC truth, SC etc.
0138   // ===================================================
0139   // =======================================================================================
0140   // =======================================================================================
0141   Handle<HepMCProduct> hepMC;
0142   evt.getByToken(hepMCLabel_, hepMC);
0143
0144   Labels l;
0145   labelsForToken(hepMCLabel_, l);
0146
0147   [[clang::suppress]]
0148   const HepMC::GenEvent *genEvent = hepMC->GetEvent();
0149   if (!(hepMC.isValid())) {
0150     LogInfo("EnergyScaleAnalyzer") << "Could not get MC Product!";
0151     return;
0152   }
0153
0154   //=======================For Vertex correction
0155   std::vector<Handle<HepMCProduct>> evtHandles;
0156
0157   //evt.getManyByType(evtHandles);
0158   throw cms::Exception("UnsupportedFunction") << "EnergyScaleAnalyzer::analyze: "
0159                                               << "getManyByType has not been supported by the Framework since 2015. "
0160                                               << "This module has been broken since then. Maybe it should be deleted. "
0161                                               << "Another possibility is to upgrade to use GetterOfProducts instead.";
0162
0163   for (unsigned int i = 0; i < evtHandles.size(); ++i) {
0164     if (evtHandles[i].isValid()) {
0165       const HepMC::GenEvent *evt = evtHandles[i]->GetEvent();
0166
0167       // take only 1st vertex for now - it's been tested only of PGuns...
0168       //
0169       HepMC::GenEvent::vertex_const_iterator vtx = evt->vertices_begin();
0170       if (evtHandles[i].provenance()->moduleLabel() == std::string(l.module)) {
0171         // Corrdinates of Vertex w.r.o. the point (0,0,0)
0172         xVtx_ = 0.1 * (*vtx)->position().x();
0173         yVtx_ = 0.1 * (*vtx)->position().y();
0174         zVtx_ = 0.1 * (*vtx)->position().z();
0175       }
0176     }
0177   }
0178   //==============================================================================
0179   // Get handle to SC collections
0180
0181   Handle<reco::SuperClusterCollection> hybridSuperClusters;
0182   try {
0183     evt.getByToken(hybridSuperClusters_token, hybridSuperClusters);
0184   } catch (cms::Exception &ex) {
0185     edm::LogError("EnergyScaleAnalyzer") << "Can't get collection with producer hybridSuperClusters.";
0186   }
0187
0188   Handle<reco::SuperClusterCollection> dynamicHybridSuperClusters;
0189   try {
0190     evt.getByToken(dynamicHybridSuperClusters_token, dynamicHybridSuperClusters);
0191   } catch (cms::Exception &ex) {
0192     edm::LogError("EnergyScaleAnalyzer") << "Can't get collection with producer dynamicHybridSuperClusters.";
0193   }
0194
0195   Handle<reco::SuperClusterCollection> fixedMatrixSuperClustersWithPS;
0196   try {
0197     evt.getByToken(fixedMatrixSuperClustersWithPreshower_token, fixedMatrixSuperClustersWithPS);
0198   } catch (cms::Exception &ex) {
0199     edm::LogError("EnergyScaleAnalyzer") << "Can't get collection with producer "
0200                                             "fixedMatrixSuperClustersWithPreshower.";
0201   }
0202
0203   // Corrected collections
0204   Handle<reco::SuperClusterCollection> correctedHybridSC;
0205   try {
0206     evt.getByToken(correctedHybridSuperClusters_token, correctedHybridSC);
0207   } catch (cms::Exception &ex) {
0208     edm::LogError("EnergyScaleAnalyzer") << "Can't get collection with producer correctedHybridSuperClusters.";
0209   }
0210
0211   Handle<reco::SuperClusterCollection> correctedDynamicHybridSC;
0212   try {
0213     evt.getByToken(correctedDynamicHybridSuperClusters_token, correctedDynamicHybridSC);
0214   } catch (cms::Exception &ex) {
0215     edm::LogError("EnergyScaleAnalyzer") << "Can't get collection with producer "
0216                                             "correctedDynamicHybridSuperClusters.";
0217   }
0218
0219   Handle<reco::SuperClusterCollection> correctedFixedMatrixSCWithPS;
0220   try {
0221     evt.getByToken(correctedFixedMatrixSuperClustersWithPreshower_token, correctedFixedMatrixSCWithPS);
0222   } catch (cms::Exception &ex) {
0223     edm::LogError("EnergyScaleAnalyzer") << "Can't get collection with producer "
0224                                             "correctedFixedMatrixSuperClustersWithPreshower.";
0225   }
0226
0227   const reco::SuperClusterCollection *dSC = dynamicHybridSuperClusters.product();
0228   const reco::SuperClusterCollection *hSC = hybridSuperClusters.product();
0229   const reco::SuperClusterCollection *fmSC = fixedMatrixSuperClustersWithPS.product();
0230   const reco::SuperClusterCollection *chSC = correctedHybridSC.product();
0231   const reco::SuperClusterCollection *cdSC = correctedDynamicHybridSC.product();
0232   const reco::SuperClusterCollection *cfmSC = correctedFixedMatrixSCWithPS.product();
0233
0234   // -----------------------  Print outs for debugging
0235   /*
0236   std::cout << "MC truth" << std::endl;
0237   int counterI = 0;
0238   for(  HepMC::GenEvent::particle_const_iterator p =
0239   genEvent->particles_begin(); p != genEvent->particles_end(); ++p ) { if (
0240   fabs((*p)->momentum().eta()) < 1.5 ) { std::cout << ++counterI << " " <<
0241   (*p)->momentum().e() << " "
0242                 << (*p)->momentum().phi() << " " << (*p)->momentum().eta() <<
0243   std::endl;
0244     }
0245   }
0246
0247   std::cout << "Standard clusters" << std::endl;
0248   counterI = 0;
0249   for(reco::SuperClusterCollection::const_iterator em = hSC->begin();
0250       em != hSC->end(); ++em )
0251     std::cout << ++counterI << " " << em->energy() << " " <<
0252   em->position().phi() << " " << em->position().eta() << std::endl;
0253
0254   std::cout << "Dynamic clusters" << std::endl;
0255   counterI = 0;
0256   for(reco::SuperClusterCollection::const_iterator em = dSC->begin();
0257       em != dSC->end(); ++em )
0258     std::cout << ++counterI << " " << em->energy() << " " <<
0259   em->position().phi() << " " << em->position().eta() << std::endl;
0260
0261   std::cout << "FixedMatrix clusters with PS" << std::endl;
0262   counterI = 0;
0263   for(reco::SuperClusterCollection::const_iterator em = fmSC->begin();
0264       em != fmSC->end(); ++em )
0265     std::cout << ++counterI << " " << em->energy() << " " <<
0266   em->position().phi() << " " << em->position().eta() << std::endl;
0267   */
0268   // -----------------------------
0269   //=====================================================================
0270   // All containers are loaded, perform the analysis
0271   //====================================================================
0272
0273   // --------------------------- Store MC particles
0274   HepMC::GenEvent::particle_const_iterator p = genEvent->particles_begin();
0275
0276   // Search for MC electrons or photons that satisfy the criteria
0277   float min_eT = 5;
0278   float max_eta = 2.5;
0279
0280   std::vector<HepMC::GenParticle *> mcParticles;
0281   // int counter = 0;
0282   for (HepMC::GenEvent::particle_const_iterator p = genEvent->particles_begin(); p != genEvent->particles_end(); ++p) {
0283     // LogInfo("EnergyScaleAnalyzer") << "Particle " << ++counter
0284     //<< " PDG ID = " << (*p)->pdg_id() << " pT = " << (*p)->momentum().perp();
0285     // require photon or electron
0286     if ((*p)->pdg_id() != 22 && abs((*p)->pdg_id()) != 11)
0287       continue;
0288
0289     // require selection criteria
0290     bool satisfySelectionCriteria = (*p)->momentum().perp() > min_eT && fabs((*p)->momentum().eta()) < max_eta;
0291
0292     if (!satisfySelectionCriteria)
0293       continue;
0294
0295     // EM MC particle is found, save it in the vector
0296     mcParticles.push_back(*p);
0297   }
0298   // separation in dR between 2 first MC particles
0299   // should not be used for MC samples with > 2 em objects generated!
0300   if (mcParticles.size() == 2) {
0301     HepMC::GenParticle *mc1 = mcParticles[0];
0302     HepMC::GenParticle *mc2 = mcParticles[1];
0303     tree_.mc_sep =
0304         kinem::delta_R(mc1->momentum().eta(), mc1->momentum().phi(), mc2->momentum().eta(), mc2->momentum().phi());
0305   } else
0306     tree_.mc_sep = -100;
0307
0308   // now loop over MC particles, find the match with SC and do everything we
0309   // need then save info in the tree for every MC particle
0310   for (std::vector<HepMC::GenParticle *>::const_iterator p = mcParticles.begin(); p != mcParticles.end(); ++p) {
0311     HepMC::GenParticle *mc = *p;
0312
0313     // Fill MC information
0314     tree_.mc_npar = mcParticles.size();
0315     tree_.parID = mc->pdg_id();
0316     tree_.mc_e = mc->momentum().e();
0317     tree_.mc_et = mc->momentum().e() * sin(mc->momentum().theta());
0318     tree_.mc_phi = mc->momentum().phi();
0319     tree_.mc_eta = mc->momentum().eta();
0320     tree_.mc_theta = mc->momentum().theta();
0321
0322     // Call function to fill tree
0323     // scType coprreponds:
0324     // HybridSuperCluster                     -- 1
0325     // DynamicHybridSuperCluster              -- 2
0326     // FixedMatrixSuperClustersWithPreshower  -- 3
0327
0328     fillTree(hSC, chSC, mc, tree_, xVtx_, yVtx_, zVtx_, 1);
0329     //    std::cout << " TYPE " << 1 << " : " << tree_.em_e << " : " <<
0330     //    tree_.em_phi << " : " << tree_.em_eta << std::endl;
0331
0332     fillTree(dSC, cdSC, mc, tree_, xVtx_, yVtx_, zVtx_, 2);
0333     //    std::cout << " TYPE " << 2 << " : " << tree_.em_e << " : " <<
0334     //    tree_.em_phi << " : " << tree_.em_eta << std::endl;
0335
0336     fillTree(fmSC, cfmSC, mc, tree_, xVtx_, yVtx_, zVtx_, 3);
0337     //    std::cout << " TYPE " << 3 << " : " << tree_.em_e << " : " <<
0338     //    tree_.em_phi << " : " << tree_.em_eta << std::endl;
0339
0340     //   mytree_->Fill();
0341   }  // loop over particles
0342 }
0343
0344 void EnergyScaleAnalyzer::fillTree(const reco::SuperClusterCollection *scColl,
0345                                    const reco::SuperClusterCollection *corrSCColl,
0346                                    HepMC::GenParticle *mc,
0347                                    tree_structure_ &tree_,
0348                                    float xV,
0349                                    float yV,
0350                                    float zV,
0351                                    int scType) {
0352   // -----------------------------  SuperClusters before energy correction
0353   reco::SuperClusterCollection::const_iterator em = scColl->end();
0354   float energyMax = -100.0;  // dummy energy of the matched SC
0355   for (reco::SuperClusterCollection::const_iterator aClus = scColl->begin(); aClus != scColl->end(); ++aClus) {
0356     // check the matching
0357     float dR =
0358         kinem::delta_R(mc->momentum().eta(), mc->momentum().phi(), aClus->position().eta(), aClus->position().phi());
0359     if (dR < 0.4) {  // a rather loose matching cut
0360       float energy = aClus->energy();
0361       if (energy < energyMax)
0362         continue;
0363       energyMax = energy;
0364       em = aClus;
0365     }
0366   }
0367
0368   if (em == scColl->end()) {
0369     //    std::cout << "No matching SC with type " << scType << " was found for
0370     //    MC particle! "  << std::endl; std::cout << "Going to next type of SC.
0371     //    " << std::endl;
0372     return;
0373   }
0374   //  ------------
0375
0376   tree_.em_scType = scType;
0377
0378   tree_.em_isInCrack = 0;
0379   double emAbsEta = fabs(em->position().eta());
0380   // copied from
0381   // RecoEgama/EgammaElectronAlgos/src/EgammaElectronClassification.cc
0382   if (emAbsEta < 0.018 || (emAbsEta > 0.423 && emAbsEta < 0.461) || (emAbsEta > 0.770 && emAbsEta < 0.806) ||
0383       (emAbsEta > 1.127 && emAbsEta < 1.163) || (emAbsEta > 1.460 && emAbsEta < 1.558))
0384     tree_.em_isInCrack = 1;
0385
0386   tree_.em_dR = kinem::delta_R(mc->momentum().eta(), mc->momentum().phi(), em->position().eta(), em->position().phi());
0387   tree_.em_e = em->energy();
0388   tree_.em_et = em->energy() * sin(em->position().theta());
0389   tree_.em_phi = em->position().phi();
0390   tree_.em_eta = em->position().eta();
0391   tree_.em_theta = em->position().theta();
0392   tree_.em_nCell = em->size();
0393   tree_.em_nBC = em->clustersSize();
0394
0395   // Get physics e, et etc:
0396   // Coordinates of EM object with respect of the point (0,0,0)
0397   xClust_zero_ = em->position().x();
0398   yClust_zero_ = em->position().y();
0399   zClust_zero_ = em->position().z();
0400   // Coordinates of EM object w.r.o. the Vertex position
0401   xClust_vtx_ = xClust_zero_ - xV;
0402   yClust_vtx_ = yClust_zero_ - yV;
0403   zClust_vtx_ = zClust_zero_ - zV;
0404
0405   energyMax_ = em->energy();
0406   thetaMax_ = em->position().theta();
0407   etaMax_ = em->position().eta();
0408   phiMax_ = em->position().phi();
0409   eTMax_ = energyMax_ * sin(thetaMax_);
0410   if (phiMax_ < 0)
0411     phiMax_ += 2 * 3.14159;
0412
0413   rClust_vtx_ = sqrt(xClust_vtx_ * xClust_vtx_ + yClust_vtx_ * yClust_vtx_ + zClust_vtx_ * zClust_vtx_);
0414   thetaMaxVtx_ = acos(zClust_vtx_ / rClust_vtx_);
0415   etaMaxVtx_ = -log(tan(thetaMaxVtx_ / 2));
0416   eTMaxVtx_ = energyMax_ * sin(thetaMaxVtx_);
0417   phiMaxVtx_ = atan2(yClust_vtx_, xClust_vtx_);
0418   if (phiMaxVtx_ < 0)
0419     phiMaxVtx_ += 2 * 3.14159;
0420   //=============================
0421   // parametres of EM object after vertex correction
0422   tree_.em_pet = eTMaxVtx_;
0423   tree_.em_pe = tree_.em_pet / sin(thetaMaxVtx_);
0424   tree_.em_peta = etaMaxVtx_;
0425   tree_.em_ptheta = thetaMaxVtx_;
0426
0427   //-------------------------------   Get SC after energy correction
0428   em = corrSCColl->end();
0429   energyMax = -100.0;  // dummy energy of the matched SC
0430   for (reco::SuperClusterCollection::const_iterator aClus = corrSCColl->begin(); aClus != corrSCColl->end(); ++aClus) {
0431     // check the matching
0432     float dR =
0433         kinem::delta_R(mc->momentum().eta(), mc->momentum().phi(), aClus->position().eta(), aClus->position().phi());
0434     if (dR < 0.4) {
0435       float energy = aClus->energy();
0436       if (energy < energyMax)
0437         continue;
0438       energyMax = energy;
0439       em = aClus;
0440     }
0441   }
0442
0443   if (em == corrSCColl->end()) {
0444     //    std::cout << "No matching corrected SC with type " << scType << " was
0445     //    found for MC particle! "  << std::endl; std::cout << "Going to next
0446     //    type of SC. " << std::endl;
0447     return;
0448   }
0449   //  ------------
0450
0451   /// fill tree with kinematic variables of corrected Super Cluster
0452   tree_.emCorr_e = em->energy();
0453   tree_.emCorr_et = em->energy() * sin(em->position().theta());
0454   tree_.emCorr_phi = em->position().phi();
0455   tree_.emCorr_eta = em->position().eta();
0456   tree_.emCorr_theta = em->position().theta();
0457
0458   // =========== Eta and Theta wrt Vertex does not change after energy
0459   // corrections are applied
0460   // =========== So, no need to calculate them again
0461
0462   tree_.emCorr_peta = tree_.em_peta;
0463   tree_.emCorr_ptheta = tree_.em_ptheta;
0464   tree_.emCorr_pet = tree_.emCorr_e / cosh(tree_.emCorr_peta);
0465
0466   tree_.em_pw = em->phiWidth();
0467   tree_.em_ew = em->etaWidth();
0468   tree_.em_br = tree_.em_pw / tree_.em_ew;
0469
0470   mytree_->Fill();
0471 }
0472
0473 //==========================================================================
0474 void EnergyScaleAnalyzer::endJob() {
0475   //========================================================================
0476   // Fill ROOT tree
0477   rootFile_->Write();
0478 }
0479
0480 DEFINE_FWK_MODULE(EnergyScaleAnalyzer);