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   const HepMC::GenEvent *genEvent = hepMC->GetEvent();
0148   if (!(hepMC.isValid())) {
0149     LogInfo("EnergyScaleAnalyzer") << "Could not get MC Product!";
0150     return;
0151   }
0152
0153   //=======================For Vertex correction
0154   std::vector<Handle<HepMCProduct>> evtHandles;
0155
0156   //evt.getManyByType(evtHandles);
0157   throw cms::Exception("UnsupportedFunction") << "EnergyScaleAnalyzer::analyze: "
0158                                               << "getManyByType has not been supported by the Framework since 2015. "
0159                                               << "This module has been broken since then. Maybe it should be deleted. "
0160                                               << "Another possibility is to upgrade to use GetterOfProducts instead.";
0161
0162   for (unsigned int i = 0; i < evtHandles.size(); ++i) {
0163     if (evtHandles[i].isValid()) {
0164       const HepMC::GenEvent *evt = evtHandles[i]->GetEvent();
0165
0166       // take only 1st vertex for now - it's been tested only of PGuns...
0167       //
0168       HepMC::GenEvent::vertex_const_iterator vtx = evt->vertices_begin();
0169       if (evtHandles[i].provenance()->moduleLabel() == std::string(l.module)) {
0170         // Corrdinates of Vertex w.r.o. the point (0,0,0)
0171         xVtx_ = 0.1 * (*vtx)->position().x();
0172         yVtx_ = 0.1 * (*vtx)->position().y();
0173         zVtx_ = 0.1 * (*vtx)->position().z();
0174       }
0175     }
0176   }
0177   //==============================================================================
0178   // Get handle to SC collections
0179
0180   Handle<reco::SuperClusterCollection> hybridSuperClusters;
0181   try {
0182     evt.getByToken(hybridSuperClusters_token, hybridSuperClusters);
0183   } catch (cms::Exception &ex) {
0184     edm::LogError("EnergyScaleAnalyzer") << "Can't get collection with producer hybridSuperClusters.";
0185   }
0186
0187   Handle<reco::SuperClusterCollection> dynamicHybridSuperClusters;
0188   try {
0189     evt.getByToken(dynamicHybridSuperClusters_token, dynamicHybridSuperClusters);
0190   } catch (cms::Exception &ex) {
0191     edm::LogError("EnergyScaleAnalyzer") << "Can't get collection with producer dynamicHybridSuperClusters.";
0192   }
0193
0194   Handle<reco::SuperClusterCollection> fixedMatrixSuperClustersWithPS;
0195   try {
0196     evt.getByToken(fixedMatrixSuperClustersWithPreshower_token, fixedMatrixSuperClustersWithPS);
0197   } catch (cms::Exception &ex) {
0198     edm::LogError("EnergyScaleAnalyzer") << "Can't get collection with producer "
0199                                             "fixedMatrixSuperClustersWithPreshower.";
0200   }
0201
0202   // Corrected collections
0203   Handle<reco::SuperClusterCollection> correctedHybridSC;
0204   try {
0205     evt.getByToken(correctedHybridSuperClusters_token, correctedHybridSC);
0206   } catch (cms::Exception &ex) {
0207     edm::LogError("EnergyScaleAnalyzer") << "Can't get collection with producer correctedHybridSuperClusters.";
0208   }
0209
0210   Handle<reco::SuperClusterCollection> correctedDynamicHybridSC;
0211   try {
0212     evt.getByToken(correctedDynamicHybridSuperClusters_token, correctedDynamicHybridSC);
0213   } catch (cms::Exception &ex) {
0214     edm::LogError("EnergyScaleAnalyzer") << "Can't get collection with producer "
0215                                             "correctedDynamicHybridSuperClusters.";
0216   }
0217
0218   Handle<reco::SuperClusterCollection> correctedFixedMatrixSCWithPS;
0219   try {
0220     evt.getByToken(correctedFixedMatrixSuperClustersWithPreshower_token, correctedFixedMatrixSCWithPS);
0221   } catch (cms::Exception &ex) {
0222     edm::LogError("EnergyScaleAnalyzer") << "Can't get collection with producer "
0223                                             "correctedFixedMatrixSuperClustersWithPreshower.";
0224   }
0225
0226   const reco::SuperClusterCollection *dSC = dynamicHybridSuperClusters.product();
0227   const reco::SuperClusterCollection *hSC = hybridSuperClusters.product();
0228   const reco::SuperClusterCollection *fmSC = fixedMatrixSuperClustersWithPS.product();
0229   const reco::SuperClusterCollection *chSC = correctedHybridSC.product();
0230   const reco::SuperClusterCollection *cdSC = correctedDynamicHybridSC.product();
0231   const reco::SuperClusterCollection *cfmSC = correctedFixedMatrixSCWithPS.product();
0232
0233   // -----------------------  Print outs for debugging
0234   /*
0235   std::cout << "MC truth" << std::endl;
0236   int counterI = 0;
0237   for(  HepMC::GenEvent::particle_const_iterator p =
0238   genEvent->particles_begin(); p != genEvent->particles_end(); ++p ) { if (
0239   fabs((*p)->momentum().eta()) < 1.5 ) { std::cout << ++counterI << " " <<
0240   (*p)->momentum().e() << " "
0241                 << (*p)->momentum().phi() << " " << (*p)->momentum().eta() <<
0242   std::endl;
0243     }
0244   }
0245
0246   std::cout << "Standard clusters" << std::endl;
0247   counterI = 0;
0248   for(reco::SuperClusterCollection::const_iterator em = hSC->begin();
0249       em != hSC->end(); ++em )
0250     std::cout << ++counterI << " " << em->energy() << " " <<
0251   em->position().phi() << " " << em->position().eta() << std::endl;
0252
0253   std::cout << "Dynamic clusters" << std::endl;
0254   counterI = 0;
0255   for(reco::SuperClusterCollection::const_iterator em = dSC->begin();
0256       em != dSC->end(); ++em )
0257     std::cout << ++counterI << " " << em->energy() << " " <<
0258   em->position().phi() << " " << em->position().eta() << std::endl;
0259
0260   std::cout << "FixedMatrix clusters with PS" << std::endl;
0261   counterI = 0;
0262   for(reco::SuperClusterCollection::const_iterator em = fmSC->begin();
0263       em != fmSC->end(); ++em )
0264     std::cout << ++counterI << " " << em->energy() << " " <<
0265   em->position().phi() << " " << em->position().eta() << std::endl;
0266   */
0267   // -----------------------------
0268   //=====================================================================
0269   // All containers are loaded, perform the analysis
0270   //====================================================================
0271
0272   // --------------------------- Store MC particles
0273   HepMC::GenEvent::particle_const_iterator p = genEvent->particles_begin();
0274
0275   // Search for MC electrons or photons that satisfy the criteria
0276   float min_eT = 5;
0277   float max_eta = 2.5;
0278
0279   std::vector<HepMC::GenParticle *> mcParticles;
0280   // int counter = 0;
0281   for (HepMC::GenEvent::particle_const_iterator p = genEvent->particles_begin(); p != genEvent->particles_end(); ++p) {
0282     // LogInfo("EnergyScaleAnalyzer") << "Particle " << ++counter
0283     //<< " PDG ID = " << (*p)->pdg_id() << " pT = " << (*p)->momentum().perp();
0284     // require photon or electron
0285     if ((*p)->pdg_id() != 22 && abs((*p)->pdg_id()) != 11)
0286       continue;
0287
0288     // require selection criteria
0289     bool satisfySelectionCriteria = (*p)->momentum().perp() > min_eT && fabs((*p)->momentum().eta()) < max_eta;
0290
0291     if (!satisfySelectionCriteria)
0292       continue;
0293
0294     // EM MC particle is found, save it in the vector
0295     mcParticles.push_back(*p);
0296   }
0297   // separation in dR between 2 first MC particles
0298   // should not be used for MC samples with > 2 em objects generated!
0299   if (mcParticles.size() == 2) {
0300     HepMC::GenParticle *mc1 = mcParticles[0];
0301     HepMC::GenParticle *mc2 = mcParticles[1];
0302     tree_.mc_sep =
0303         kinem::delta_R(mc1->momentum().eta(), mc1->momentum().phi(), mc2->momentum().eta(), mc2->momentum().phi());
0304   } else
0305     tree_.mc_sep = -100;
0306
0307   // now loop over MC particles, find the match with SC and do everything we
0308   // need then save info in the tree for every MC particle
0309   for (std::vector<HepMC::GenParticle *>::const_iterator p = mcParticles.begin(); p != mcParticles.end(); ++p) {
0310     HepMC::GenParticle *mc = *p;
0311
0312     // Fill MC information
0313     tree_.mc_npar = mcParticles.size();
0314     tree_.parID = mc->pdg_id();
0315     tree_.mc_e = mc->momentum().e();
0316     tree_.mc_et = mc->momentum().e() * sin(mc->momentum().theta());
0317     tree_.mc_phi = mc->momentum().phi();
0318     tree_.mc_eta = mc->momentum().eta();
0319     tree_.mc_theta = mc->momentum().theta();
0320
0321     // Call function to fill tree
0322     // scType coprreponds:
0323     // HybridSuperCluster                     -- 1
0324     // DynamicHybridSuperCluster              -- 2
0325     // FixedMatrixSuperClustersWithPreshower  -- 3
0326
0327     fillTree(hSC, chSC, mc, tree_, xVtx_, yVtx_, zVtx_, 1);
0328     //    std::cout << " TYPE " << 1 << " : " << tree_.em_e << " : " <<
0329     //    tree_.em_phi << " : " << tree_.em_eta << std::endl;
0330
0331     fillTree(dSC, cdSC, mc, tree_, xVtx_, yVtx_, zVtx_, 2);
0332     //    std::cout << " TYPE " << 2 << " : " << tree_.em_e << " : " <<
0333     //    tree_.em_phi << " : " << tree_.em_eta << std::endl;
0334
0335     fillTree(fmSC, cfmSC, mc, tree_, xVtx_, yVtx_, zVtx_, 3);
0336     //    std::cout << " TYPE " << 3 << " : " << tree_.em_e << " : " <<
0337     //    tree_.em_phi << " : " << tree_.em_eta << std::endl;
0338
0339     //   mytree_->Fill();
0340   }  // loop over particles
0341 }
0342
0343 void EnergyScaleAnalyzer::fillTree(const reco::SuperClusterCollection *scColl,
0344                                    const reco::SuperClusterCollection *corrSCColl,
0345                                    HepMC::GenParticle *mc,
0346                                    tree_structure_ &tree_,
0347                                    float xV,
0348                                    float yV,
0349                                    float zV,
0350                                    int scType) {
0351   // -----------------------------  SuperClusters before energy correction
0352   reco::SuperClusterCollection::const_iterator em = scColl->end();
0353   float energyMax = -100.0;  // dummy energy of the matched SC
0354   for (reco::SuperClusterCollection::const_iterator aClus = scColl->begin(); aClus != scColl->end(); ++aClus) {
0355     // check the matching
0356     float dR =
0357         kinem::delta_R(mc->momentum().eta(), mc->momentum().phi(), aClus->position().eta(), aClus->position().phi());
0358     if (dR < 0.4) {  // a rather loose matching cut
0359       float energy = aClus->energy();
0360       if (energy < energyMax)
0361         continue;
0362       energyMax = energy;
0363       em = aClus;
0364     }
0365   }
0366
0367   if (em == scColl->end()) {
0368     //    std::cout << "No matching SC with type " << scType << " was found for
0369     //    MC particle! "  << std::endl; std::cout << "Going to next type of SC.
0370     //    " << std::endl;
0371     return;
0372   }
0373   //  ------------
0374
0375   tree_.em_scType = scType;
0376
0377   tree_.em_isInCrack = 0;
0378   double emAbsEta = fabs(em->position().eta());
0379   // copied from
0380   // RecoEgama/EgammaElectronAlgos/src/EgammaElectronClassification.cc
0381   if (emAbsEta < 0.018 || (emAbsEta > 0.423 && emAbsEta < 0.461) || (emAbsEta > 0.770 && emAbsEta < 0.806) ||
0382       (emAbsEta > 1.127 && emAbsEta < 1.163) || (emAbsEta > 1.460 && emAbsEta < 1.558))
0383     tree_.em_isInCrack = 1;
0384
0385   tree_.em_dR = kinem::delta_R(mc->momentum().eta(), mc->momentum().phi(), em->position().eta(), em->position().phi());
0386   tree_.em_e = em->energy();
0387   tree_.em_et = em->energy() * sin(em->position().theta());
0388   tree_.em_phi = em->position().phi();
0389   tree_.em_eta = em->position().eta();
0390   tree_.em_theta = em->position().theta();
0391   tree_.em_nCell = em->size();
0392   tree_.em_nBC = em->clustersSize();
0393
0394   // Get physics e, et etc:
0395   // Coordinates of EM object with respect of the point (0,0,0)
0396   xClust_zero_ = em->position().x();
0397   yClust_zero_ = em->position().y();
0398   zClust_zero_ = em->position().z();
0399   // Coordinates of EM object w.r.o. the Vertex position
0400   xClust_vtx_ = xClust_zero_ - xV;
0401   yClust_vtx_ = yClust_zero_ - yV;
0402   zClust_vtx_ = zClust_zero_ - zV;
0403
0404   energyMax_ = em->energy();
0405   thetaMax_ = em->position().theta();
0406   etaMax_ = em->position().eta();
0407   phiMax_ = em->position().phi();
0408   eTMax_ = energyMax_ * sin(thetaMax_);
0409   if (phiMax_ < 0)
0410     phiMax_ += 2 * 3.14159;
0411
0412   rClust_vtx_ = sqrt(xClust_vtx_ * xClust_vtx_ + yClust_vtx_ * yClust_vtx_ + zClust_vtx_ * zClust_vtx_);
0413   thetaMaxVtx_ = acos(zClust_vtx_ / rClust_vtx_);
0414   etaMaxVtx_ = -log(tan(thetaMaxVtx_ / 2));
0415   eTMaxVtx_ = energyMax_ * sin(thetaMaxVtx_);
0416   phiMaxVtx_ = atan2(yClust_vtx_, xClust_vtx_);
0417   if (phiMaxVtx_ < 0)
0418     phiMaxVtx_ += 2 * 3.14159;
0419   //=============================
0420   // parametres of EM object after vertex correction
0421   tree_.em_pet = eTMaxVtx_;
0422   tree_.em_pe = tree_.em_pet / sin(thetaMaxVtx_);
0423   tree_.em_peta = etaMaxVtx_;
0424   tree_.em_ptheta = thetaMaxVtx_;
0425
0426   //-------------------------------   Get SC after energy correction
0427   em = corrSCColl->end();
0428   energyMax = -100.0;  // dummy energy of the matched SC
0429   for (reco::SuperClusterCollection::const_iterator aClus = corrSCColl->begin(); aClus != corrSCColl->end(); ++aClus) {
0430     // check the matching
0431     float dR =
0432         kinem::delta_R(mc->momentum().eta(), mc->momentum().phi(), aClus->position().eta(), aClus->position().phi());
0433     if (dR < 0.4) {
0434       float energy = aClus->energy();
0435       if (energy < energyMax)
0436         continue;
0437       energyMax = energy;
0438       em = aClus;
0439     }
0440   }
0441
0442   if (em == corrSCColl->end()) {
0443     //    std::cout << "No matching corrected SC with type " << scType << " was
0444     //    found for MC particle! "  << std::endl; std::cout << "Going to next
0445     //    type of SC. " << std::endl;
0446     return;
0447   }
0448   //  ------------
0449
0450   /// fill tree with kinematic variables of corrected Super Cluster
0451   tree_.emCorr_e = em->energy();
0452   tree_.emCorr_et = em->energy() * sin(em->position().theta());
0453   tree_.emCorr_phi = em->position().phi();
0454   tree_.emCorr_eta = em->position().eta();
0455   tree_.emCorr_theta = em->position().theta();
0456
0457   // =========== Eta and Theta wrt Vertex does not change after energy
0458   // corrections are applied
0459   // =========== So, no need to calculate them again
0460
0461   tree_.emCorr_peta = tree_.em_peta;
0462   tree_.emCorr_ptheta = tree_.em_ptheta;
0463   tree_.emCorr_pet = tree_.emCorr_e / cosh(tree_.emCorr_peta);
0464
0465   tree_.em_pw = em->phiWidth();
0466   tree_.em_ew = em->etaWidth();
0467   tree_.em_br = tree_.em_pw / tree_.em_ew;
0468
0469   mytree_->Fill();
0470 }
0471
0472 //==========================================================================
0473 void EnergyScaleAnalyzer::endJob() {
0474   //========================================================================
0475   // Fill ROOT tree
0476   rootFile_->Write();
0477 }
0478
0479 DEFINE_FWK_MODULE(EnergyScaleAnalyzer);