File indexing completed on 2024-10-08 05:12:07
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018 #include "Validation/EcalClusters/interface/EnergyScaleAnalyzer.h"
0019
0020
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
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
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
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");
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:";
0112 treeVariables += "em_dR/F:";
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:";
0116 treeVariables += "em_pet/F:em_pe:em_peta:em_ptheta:";
0117
0118
0119 treeVariables += "emCorr_e/F:emCorr_et:emCorr_eta:emCorr_phi:emCorr_theta:";
0120
0121
0122 treeVariables += "emCorr_pet/F:emCorr_peta:emCorr_ptheta:";
0123
0124 treeVariables += "em_pw/F:em_ew:em_br";
0125
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;
0133 using namespace std;
0134
0135
0136
0137
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
0155 std::vector<Handle<HepMCProduct>> evtHandles;
0156
0157
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
0168
0169 HepMC::GenEvent::vertex_const_iterator vtx = evt->vertices_begin();
0170 if (evtHandles[i].provenance()->moduleLabel() == std::string(l.module)) {
0171
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
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
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
0235
0236
0237
0238
0239
0240
0241
0242
0243
0244
0245
0246
0247
0248
0249
0250
0251
0252
0253
0254
0255
0256
0257
0258
0259
0260
0261
0262
0263
0264
0265
0266
0267
0268
0269
0270
0271
0272
0273
0274 HepMC::GenEvent::particle_const_iterator p = genEvent->particles_begin();
0275
0276
0277 float min_eT = 5;
0278 float max_eta = 2.5;
0279
0280 std::vector<HepMC::GenParticle *> mcParticles;
0281
0282 for (HepMC::GenEvent::particle_const_iterator p = genEvent->particles_begin(); p != genEvent->particles_end(); ++p) {
0283
0284
0285
0286 if ((*p)->pdg_id() != 22 && abs((*p)->pdg_id()) != 11)
0287 continue;
0288
0289
0290 bool satisfySelectionCriteria = (*p)->momentum().perp() > min_eT && fabs((*p)->momentum().eta()) < max_eta;
0291
0292 if (!satisfySelectionCriteria)
0293 continue;
0294
0295
0296 mcParticles.push_back(*p);
0297 }
0298
0299
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
0309
0310 for (std::vector<HepMC::GenParticle *>::const_iterator p = mcParticles.begin(); p != mcParticles.end(); ++p) {
0311 HepMC::GenParticle *mc = *p;
0312
0313
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
0323
0324
0325
0326
0327
0328 fillTree(hSC, chSC, mc, tree_, xVtx_, yVtx_, zVtx_, 1);
0329
0330
0331
0332 fillTree(dSC, cdSC, mc, tree_, xVtx_, yVtx_, zVtx_, 2);
0333
0334
0335
0336 fillTree(fmSC, cfmSC, mc, tree_, xVtx_, yVtx_, zVtx_, 3);
0337
0338
0339
0340
0341 }
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
0353 reco::SuperClusterCollection::const_iterator em = scColl->end();
0354 float energyMax = -100.0;
0355 for (reco::SuperClusterCollection::const_iterator aClus = scColl->begin(); aClus != scColl->end(); ++aClus) {
0356
0357 float dR =
0358 kinem::delta_R(mc->momentum().eta(), mc->momentum().phi(), aClus->position().eta(), aClus->position().phi());
0359 if (dR < 0.4) {
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
0370
0371
0372 return;
0373 }
0374
0375
0376 tree_.em_scType = scType;
0377
0378 tree_.em_isInCrack = 0;
0379 double emAbsEta = fabs(em->position().eta());
0380
0381
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
0396
0397 xClust_zero_ = em->position().x();
0398 yClust_zero_ = em->position().y();
0399 zClust_zero_ = em->position().z();
0400
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
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
0428 em = corrSCColl->end();
0429 energyMax = -100.0;
0430 for (reco::SuperClusterCollection::const_iterator aClus = corrSCColl->begin(); aClus != corrSCColl->end(); ++aClus) {
0431
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
0445
0446
0447 return;
0448 }
0449
0450
0451
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
0459
0460
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
0477 rootFile_->Write();
0478 }
0479
0480 DEFINE_FWK_MODULE(EnergyScaleAnalyzer);