File indexing completed on 2024-04-06 12:24:46
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017 #include <memory>
0018
0019
0020 #include "FWCore/Framework/interface/Frameworkfwd.h"
0021 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0022
0023 #include "FWCore/Framework/interface/Event.h"
0024 #include "FWCore/Framework/interface/MakerMacros.h"
0025
0026 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0027
0028
0029 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0030 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
0031 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0032 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0033
0034 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterSeverityLevelAlgo.h"
0035 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
0036
0037 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0038 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0039
0040 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0041 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0042
0043 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0044 #include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"
0045
0046 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0047 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0048
0049 #include "CLHEP/Units/PhysicalConstants.h"
0050
0051 #include "TTree.h"
0052 #include "TFile.h"
0053
0054 class testEcalClusterSeverityAlgo : public edm::one::EDAnalyzer<> {
0055 public:
0056 explicit testEcalClusterSeverityAlgo(const edm::ParameterSet&);
0057 ~testEcalClusterSeverityAlgo() override;
0058
0059 void analyze(const edm::Event&, const edm::EventSetup&) override;
0060 void endJob() override;
0061
0062 const edm::InputTag barrelClusterCollection_;
0063 const edm::InputTag endcapClusterCollection_;
0064 const edm::InputTag reducedBarrelRecHitCollection_;
0065 const edm::InputTag reducedEndcapRecHitCollection_;
0066 const edm::InputTag mcTruthCollection_;
0067
0068 struct ClusterSeverityTreeContent {
0069 #define NMAXOBJ 50
0070
0071 int nSc;
0072 float scE[NMAXOBJ];
0073 float scEta[NMAXOBJ];
0074 float scPhi[NMAXOBJ];
0075 float scGoodFraction[NMAXOBJ];
0076 float scFracAroundClosProb[NMAXOBJ];
0077 float scClosProbEta[NMAXOBJ];
0078 float scClosProbPhi[NMAXOBJ];
0079 float mcE[NMAXOBJ];
0080 float mcEta[NMAXOBJ];
0081 float mcPhi[NMAXOBJ];
0082 };
0083
0084 static void setBranchAddresses(TTree* chain, ClusterSeverityTreeContent& treeVars) {
0085 chain->SetBranchAddress("nSc", &treeVars.nSc);
0086 chain->SetBranchAddress("scE", treeVars.scE);
0087 chain->SetBranchAddress("scEta", treeVars.scEta);
0088 chain->SetBranchAddress("scPhi", treeVars.scPhi);
0089 chain->SetBranchAddress("scGoodFraction", treeVars.scGoodFraction);
0090 chain->SetBranchAddress("scFracAroundClosProb", treeVars.scFracAroundClosProb);
0091 chain->SetBranchAddress("scClosProbEta", treeVars.scClosProbEta);
0092 chain->SetBranchAddress("scClosProbPhi", treeVars.scClosProbPhi);
0093 chain->SetBranchAddress("mcE", treeVars.mcE);
0094 chain->SetBranchAddress("mcEta", treeVars.mcEta);
0095 chain->SetBranchAddress("mcPhi", treeVars.mcPhi);
0096 }
0097
0098 static void setBranches(TTree* chain, ClusterSeverityTreeContent& treeVars) {
0099 chain->Branch("nSc", &treeVars.nSc, "nSc/I");
0100 chain->Branch("scE", treeVars.scE, "scE[nSc]/F");
0101 chain->Branch("scEta", treeVars.scEta, "scEta[nSc]/F");
0102 chain->Branch("scPhi", treeVars.scPhi, "scPhi[nSc]/F");
0103 chain->Branch("scGoodFraction", treeVars.scGoodFraction, "scGoodFraction[nSc]/F");
0104 chain->Branch("scFracAroundClosProb", treeVars.scFracAroundClosProb, "scFracAroundClosProb[nSc]/F");
0105 chain->Branch("scClosProbEta", treeVars.scClosProbEta, "scClosProbEta[nSc]/F");
0106 chain->Branch("scClosProbPhi", treeVars.scClosProbPhi, "scClosProbPhi[nSc]/F");
0107 chain->Branch("mcE", treeVars.mcE, "mcE[nSc]/F");
0108 chain->Branch("mcEta", treeVars.mcEta, "mcEta[nSc]/F");
0109 chain->Branch("mcPhi", treeVars.mcPhi, "mcPhi[nSc]/F");
0110 }
0111
0112 void initializeBranches(TTree* chain, ClusterSeverityTreeContent& treeVars) {
0113 treeVars.nSc = 0;
0114 for (int i = 0; i < NMAXOBJ; ++i) {
0115 treeVars.scE[i] = 0;
0116 treeVars.scEta[i] = 0;
0117 treeVars.scPhi[i] = 0;
0118 treeVars.scGoodFraction[i] = 0;
0119 treeVars.scFracAroundClosProb[i] = 0;
0120 treeVars.scClosProbEta[i] = 0;
0121 treeVars.scClosProbPhi[i] = 0;
0122 treeVars.mcE[i] = 0;
0123 treeVars.mcEta[i] = 0;
0124 treeVars.mcPhi[i] = 0;
0125 }
0126 }
0127
0128 private:
0129 const edm::EDGetTokenT<reco::SuperClusterCollection> ebSCToken_;
0130 const edm::EDGetTokenT<reco::SuperClusterCollection> eeSCToken_;
0131 const edm::EDGetTokenT<EcalRecHitCollection> ebRecHitsToken_;
0132 const edm::EDGetTokenT<EcalRecHitCollection> eeRecHitsToken_;
0133 const edm::EDGetTokenT<edm::HepMCProduct> mcTruthToken_;
0134 const edm::ESGetToken<CaloTopology, CaloTopologyRecord> topologyToken_;
0135 const edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> severityLevelAlgoToken_;
0136
0137 std::string outputFile_;
0138 TFile* treeFile_;
0139 TTree* tree_;
0140 ClusterSeverityTreeContent myTreeVariables_;
0141 };
0142
0143 testEcalClusterSeverityAlgo::testEcalClusterSeverityAlgo(const edm::ParameterSet& ps)
0144 : barrelClusterCollection_(ps.getParameter<edm::InputTag>("barrelClusterCollection")),
0145 endcapClusterCollection_(ps.getParameter<edm::InputTag>("endcapClusterCollection")),
0146 reducedBarrelRecHitCollection_(ps.getParameter<edm::InputTag>("reducedBarrelRecHitCollection")),
0147 reducedEndcapRecHitCollection_(ps.getParameter<edm::InputTag>("reducedEndcapRecHitCollection")),
0148 mcTruthCollection_(ps.getParameter<edm::InputTag>("mcTruthCollection")),
0149 ebSCToken_(consumes<reco::SuperClusterCollection>(barrelClusterCollection_)),
0150 eeSCToken_(consumes<reco::SuperClusterCollection>(endcapClusterCollection_)),
0151 ebRecHitsToken_(consumes<EcalRecHitCollection>(reducedBarrelRecHitCollection_)),
0152 eeRecHitsToken_(consumes<EcalRecHitCollection>(reducedEndcapRecHitCollection_)),
0153 mcTruthToken_(consumes<edm::HepMCProduct>(mcTruthCollection_)),
0154 topologyToken_(esConsumes()),
0155 severityLevelAlgoToken_(esConsumes()),
0156 outputFile_(ps.getParameter<std::string>("outputFile")) {
0157 treeFile_ = new TFile(outputFile_.c_str(), "RECREATE");
0158 treeFile_->cd();
0159
0160 tree_ = new TTree("ClusterSeverityAnalysis", "ClusterSeverityAnalysis");
0161 setBranches(tree_, myTreeVariables_);
0162 }
0163
0164 testEcalClusterSeverityAlgo::~testEcalClusterSeverityAlgo() {
0165 delete treeFile_;
0166 delete tree_;
0167 }
0168
0169 void testEcalClusterSeverityAlgo::analyze(const edm::Event& ev, const edm::EventSetup& es) {
0170 initializeBranches(tree_, myTreeVariables_);
0171
0172 edm::Handle<edm::HepMCProduct> hepMC;
0173 ev.getByToken(mcTruthToken_, hepMC);
0174 const HepMC::GenEvent* myGenEvent = hepMC->GetEvent();
0175
0176 const auto& sevLv = es.getData(severityLevelAlgoToken_);
0177
0178 edm::Handle<reco::SuperClusterCollection> pEBClusters;
0179 ev.getByToken(ebSCToken_, pEBClusters);
0180 const reco::SuperClusterCollection* ebClusters = pEBClusters.product();
0181
0182 edm::Handle<reco::SuperClusterCollection> pEEClusters;
0183 ev.getByToken(eeSCToken_, pEEClusters);
0184
0185
0186 edm::Handle<EcalRecHitCollection> pEBRecHits;
0187 ev.getByToken(ebRecHitsToken_, pEBRecHits);
0188 const EcalRecHitCollection* ebRecHits = pEBRecHits.product();
0189
0190 edm::Handle<EcalRecHitCollection> pEERecHits;
0191 ev.getByToken(eeRecHitsToken_, pEERecHits);
0192
0193
0194 const auto& topology = es.getData(topologyToken_);
0195
0196 int problematicSC = 0;
0197 for (reco::SuperClusterCollection::const_iterator it = ebClusters->begin(); it != ebClusters->end(); ++it) {
0198
0199 if ((*it).energy() / cosh((*it).eta()) < 15.)
0200 continue;
0201
0202 HepMC::GenParticle* bestMcMatch = 0;
0203 for (HepMC::GenEvent::particle_const_iterator mcIter = myGenEvent->particles_begin();
0204 mcIter != myGenEvent->particles_end();
0205 mcIter++) {
0206
0207 if (std::abs((*mcIter)->pdg_id()) == 11) {
0208
0209 HepMC::GenParticle* mother = 0;
0210 if ((*mcIter)->production_vertex()) {
0211 if ((*mcIter)->production_vertex()->particles_begin(HepMC::parents) !=
0212 (*mcIter)->production_vertex()->particles_end(HepMC::parents))
0213 mother = *((*mcIter)->production_vertex()->particles_begin(HepMC::parents));
0214 }
0215 if (((mother == 0) || ((mother != 0) && (std::abs(mother->pdg_id()) == 23)) ||
0216 ((mother != 0) && (std::abs(mother->pdg_id()) == 32)) ||
0217 ((mother != 0) && (std::abs(mother->pdg_id()) == 24)))) {
0218 HepMC::GenParticle* genPc = (*mcIter);
0219 HepMC::FourVector pAssSim = genPc->momentum();
0220
0221 double ScOkRatio = 999999.;
0222
0223 double dphi = (*it).phi() - pAssSim.phi();
0224 if (fabs(dphi) > CLHEP::pi)
0225 dphi = dphi < 0 ? (CLHEP::twopi) + dphi : dphi - CLHEP::twopi;
0226 double deltaR = sqrt(pow(((*it).eta() - pAssSim.eta()), 2) + pow(dphi, 2));
0227 if (deltaR < 0.15) {
0228 double tmpScRatio = (*it).energy() / pAssSim.t();
0229 if (fabs(tmpScRatio - 1) < fabs(ScOkRatio - 1)) {
0230 ScOkRatio = tmpScRatio;
0231 bestMcMatch = genPc;
0232 }
0233 }
0234 }
0235 }
0236 }
0237
0238 if (!bestMcMatch)
0239 continue;
0240
0241 myTreeVariables_.scE[problematicSC] = (*it).energy();
0242 myTreeVariables_.scEta[problematicSC] = (*it).eta();
0243 myTreeVariables_.scPhi[problematicSC] = (*it).phi();
0244 myTreeVariables_.scGoodFraction[problematicSC] = EcalClusterSeverityLevelAlgo::goodFraction(*it, *ebRecHits, sevLv);
0245 myTreeVariables_.scFracAroundClosProb[problematicSC] =
0246 EcalClusterSeverityLevelAlgo::fractionAroundClosestProblematic(*it, *ebRecHits, &topology, sevLv);
0247 std::pair<int, int> distanceClosestProblematic =
0248 EcalClusterSeverityLevelAlgo::etaphiDistanceClosestProblematic(*it, *ebRecHits, &topology, sevLv);
0249 myTreeVariables_.scClosProbEta[problematicSC] = distanceClosestProblematic.first;
0250 myTreeVariables_.scClosProbPhi[problematicSC] = distanceClosestProblematic.second;
0251 myTreeVariables_.mcE[problematicSC] = bestMcMatch->momentum().t();
0252 myTreeVariables_.mcEta[problematicSC] = bestMcMatch->momentum().eta();
0253 myTreeVariables_.mcPhi[problematicSC] = bestMcMatch->momentum().phi();
0254 ++problematicSC;
0255 }
0256 myTreeVariables_.nSc = problematicSC;
0257 tree_->Fill();
0258 }
0259
0260 void testEcalClusterSeverityAlgo::endJob() {
0261 treeFile_->cd();
0262 tree_->Write();
0263 treeFile_->Close();
0264 }
0265
0266
0267 DEFINE_FWK_MODULE(testEcalClusterSeverityAlgo);