Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-28 05:58:08

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