EgammaCoreTools/test/testEcalClusterSeverityAlgo.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    testEcalClusterSeverityAlgo
0004 // Class:      testEcalClusterSeverityAlgo
0005 //
0006 /**\class testEcalClusterSeverityAlgo testEcalClusterSeverityAlgo.cc
0007
0008 Description: <one line class summary>
0009
0010 Implementation:
0011 <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  "Paolo Meridiani CERN CMG"
0015
0016 // system include files
0017 #include <memory>
0018
0019 // user include files
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 // to access recHits and BasicClusters
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     //MC Variables
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   // Initialize Tree
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   //const reco::SuperClusterCollection *eeClusters = pEEClusters.product();
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   //const EcalRecHitCollection *eeRecHits = pEERecHits.product();
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     //apply an et Cut
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       // select electrons
0207       if (std::abs((*mcIter)->pdg_id()) == 11) {
0208         // single primary electrons or electrons from Zs or Ws
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 //define this as a plug-in
0267 DEFINE_FWK_MODULE(testEcalClusterSeverityAlgo);