00001
00002
00003
00004
00005
00006
00007
00008
00009 #include <memory>
00010 #include <map>
00011
00012 #include "FWCore/Framework/interface/Frameworkfwd.h"
00013 #include "FWCore/Framework/interface/Event.h"
00014 #include "FWCore/Framework/interface/EventSetup.h"
00015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00016 #include "FWCore/Framework/interface/EDAnalyzer.h"
00017
00018 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
00019 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00020 #include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"
00021 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
00022
00023 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00024 #include "DataFormats/VertexReco/interface/Vertex.h"
00025
00026 #include "FWCore/ServiceRegistry/interface/Service.h"
00027 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00028 #include "TTree.h"
00029 #include "TVector2.h"
00030
00031 #include "DataFormats/Math/interface/deltaR.h"
00032
00033 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
00034 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
00035
00036 #include "RecoEcal/EgammaCoreTools/interface/Mustache.h"
00037 namespace MK = reco::MustacheKernel;
00038
00039 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h"
00040
00041 #include <algorithm>
00042 #include <memory>
00043
00044 typedef edm::ParameterSet PSet;
00045
00046 namespace {
00047 template<typename T>
00048 struct array_deleter{
00049 void operator () (T* arr) { delete [] arr; }
00050 };
00051
00052 typedef std::unary_function<const edm::Ptr<reco::PFCluster>&,
00053 double> ClusUnaryFunction;
00054
00055 struct GetSharedRecHitFraction : public ClusUnaryFunction {
00056 const edm::Ptr<reco::PFCluster> the_seed;
00057 double x_rechits_tot, x_rechits_match;
00058 GetSharedRecHitFraction(const edm::Ptr<reco::PFCluster>& s) :
00059 the_seed(s) {}
00060 double operator()(const edm::Ptr<reco::PFCluster>& x) {
00061
00062 const auto& seedHitsAndFractions =
00063 the_seed->hitsAndFractions();
00064 const auto& xHitsAndFractions =
00065 x->hitsAndFractions();
00066 x_rechits_tot = xHitsAndFractions.size();
00067 x_rechits_match = 0.0;
00068 for( const std::pair<DetId, float>& seedHit : seedHitsAndFractions ) {
00069 for( const std::pair<DetId, float>& xHit : xHitsAndFractions ) {
00070 if( seedHit.first == xHit.first ) {
00071 x_rechits_match += 1.0;
00072 }
00073 }
00074 }
00075 return x_rechits_match/x_rechits_tot;
00076 }
00077 };
00078 }
00079
00080 class PFSuperClusterTreeMaker : public edm::EDAnalyzer {
00081 typedef TTree* treeptr;
00082 public:
00083 PFSuperClusterTreeMaker(const PSet&);
00084 ~PFSuperClusterTreeMaker() {}
00085
00086 void analyze(const edm::Event&, const edm::EventSetup&);
00087 private:
00088 edm::Service<TFileService> _fs;
00089 bool _dogen;
00090 edm::InputTag _geninput;
00091 edm::InputTag _vtxsrc;
00092 edm::InputTag _scInputEB,_scInputEE;
00093 std::shared_ptr<PFEnergyCalibration> _calib;
00094 void processSuperClusterFillTree(const edm::Event&,
00095 const reco::SuperCluster&);
00096
00097 void setTreeArraysForSize(const size_t N_ECAL,const size_t N_PS);
00098 treeptr _tree;
00099 Int_t nVtx;
00100 Float_t scRawEnergy, scCalibratedEnergy, scPreshowerEnergy,
00101 scEta, scPhi, scR, scPhiWidth, scEtaWidth, scSeedRawEnergy,
00102 scSeedCalibratedEnergy, scSeedEta, scSeedPhi;
00103 Float_t genEnergy, genEta, genPhi, genDRToCentroid, genDRToSeed;
00104 Int_t N_ECALClusters;
00105 std::shared_ptr<Float_t> clusterRawEnergy, clusterCalibEnergy,
00106 clusterEta, clusterPhi, clusterDPhiToSeed, clusterDEtaToSeed,
00107 clusterDPhiToCentroid, clusterDEtaToCentroid,
00108 clusterDPhiToGen, clusterDEtaToGen, clusterHitFractionSharedWithSeed;
00109 std::shared_ptr<Int_t> clusterInMustache, clusterInDynDPhi;
00110 Int_t N_PSClusters;
00111 std::shared_ptr<Float_t> psClusterRawEnergy, psClusterEta, psClusterPhi;
00112 };
00113
00114 void PFSuperClusterTreeMaker::analyze(const edm::Event& e,
00115 const edm::EventSetup& es) {
00116 edm::Handle<reco::VertexCollection> vtcs;
00117 e.getByLabel(_vtxsrc,vtcs);
00118 if( vtcs.isValid() ) nVtx = vtcs->size();
00119 else nVtx = -1;
00120
00121 edm::Handle<reco::SuperClusterCollection> ebSCs, eeSCs;
00122 e.getByLabel(_scInputEB, ebSCs);
00123 e.getByLabel(_scInputEE, eeSCs);
00124
00125 if( ebSCs.isValid() ) {
00126 for( const auto& sc : *ebSCs ) processSuperClusterFillTree(e,sc);
00127 } else {
00128 throw cms::Exception("PFSuperClusterTreeMaker")
00129 << "Product ID for the EB SuperCluster collection was invalid!"
00130 << std::endl;
00131 }
00132
00133 if( eeSCs.isValid() ) {
00134 for( const auto& sc : *eeSCs ) processSuperClusterFillTree(e,sc);
00135 } else {
00136 throw cms::Exception("PFSuperClusterTreeMaker")
00137 << "Product ID for the EE SuperCluster collection was invalid!"
00138 << std::endl;
00139 }
00140 }
00141
00142 void PFSuperClusterTreeMaker::
00143 processSuperClusterFillTree(const edm::Event& e,
00144 const reco::SuperCluster& sc) {
00145 const int N_ECAL = sc.clustersEnd() - sc.clustersBegin();
00146 const int N_PS = ( sc.preshowerClustersEnd() -
00147 sc.preshowerClustersBegin() );
00148 if( sc.rawEnergy()/std::cosh(sc.position().Eta()) < 4.0 ) return;
00149 N_ECALClusters = std::max(0,N_ECAL - 1);
00150 N_PSClusters = N_PS;
00151 reco::GenParticleRef genmatch;
00152
00153 if( _dogen ) {
00154 edm::Handle<reco::GenParticleCollection> genp;
00155 e.getByLabel(_geninput,genp);
00156 if( genp.isValid() ) {
00157 double minDr = 1e6, this_dr;
00158 reco::GenParticleRef bestmatch;
00159 for(size_t i = 0; i < genp->size(); ++i) {
00160 this_dr = reco::deltaR(genp->at(i),sc);
00161
00162 if( this_dr < minDr && genp->at(i).pdgId() == 22) {
00163 minDr = this_dr;
00164 bestmatch = reco::GenParticleRef(genp,i);
00165 }
00166 }
00167 if( bestmatch.isNonnull() ) {
00168 genmatch = bestmatch;
00169 genEnergy = bestmatch->energy();
00170 genEta = bestmatch->eta();
00171 genPhi = bestmatch->phi();
00172 genDRToCentroid = minDr;
00173 genDRToSeed = reco::deltaR(*bestmatch,*(sc.seed()));
00174 } else {
00175 genEnergy = -1.0;
00176 genEta = 999.0;
00177 genPhi = 999.0;
00178 genDRToCentroid = 999.0;
00179 genDRToSeed = 999.0;
00180 }
00181 } else {
00182 throw cms::Exception("PFSuperClusterTreeMaker")
00183 << "Requested generator level information was not available!"
00184 << std::endl;
00185 }
00186 }
00187
00188 setTreeArraysForSize(N_ECALClusters,N_PSClusters);
00189 scRawEnergy = sc.rawEnergy();
00190 scCalibratedEnergy = sc.energy();
00191 scPreshowerEnergy = sc.preshowerEnergy();
00192 scEta = sc.position().Eta();
00193 scPhi = sc.position().Phi();
00194 scR = sc.position().R();
00195 scPhiWidth = sc.phiWidth();
00196 scEtaWidth = sc.etaWidth();
00197
00198 edm::Ptr<reco::PFCluster> theseed = edm::Ptr<reco::PFCluster>(sc.seed());
00199 GetSharedRecHitFraction fractionOfSeed(theseed);
00200 scSeedRawEnergy = theseed->energy();
00201 scSeedCalibratedEnergy = _calib->energyEm(*theseed,0.0,0.0,false);
00202 scSeedEta = theseed->eta();
00203 scSeedPhi = theseed->phi();
00204
00205 auto clusend = sc.clustersEnd();
00206 size_t iclus = 0;
00207 edm::Ptr<reco::PFCluster> pclus;
00208 for( auto clus = sc.clustersBegin(); clus != clusend; ++clus ) {
00209 pclus = edm::Ptr<reco::PFCluster>(*clus);
00210 if( theseed == pclus ) continue;
00211 clusterRawEnergy.get()[iclus] = pclus->energy();
00212 clusterCalibEnergy.get()[iclus] = _calib->energyEm(*pclus,0.0,0.0,false);
00213 clusterEta.get()[iclus] = pclus->eta();
00214 clusterPhi.get()[iclus] = pclus->phi();
00215 clusterDPhiToSeed.get()[iclus] =
00216 TVector2::Phi_mpi_pi(pclus->phi() - theseed->phi());
00217 clusterDEtaToSeed.get()[iclus] = pclus->eta() - theseed->eta();
00218 clusterDPhiToCentroid.get()[iclus] =
00219 TVector2::Phi_mpi_pi(pclus->phi() - sc.phi());
00220 clusterDEtaToCentroid.get()[iclus] = pclus->eta() - sc.eta();
00221 clusterDPhiToCentroid.get()[iclus] =
00222 TVector2::Phi_mpi_pi(pclus->phi() - sc.phi());
00223 clusterDEtaToCentroid.get()[iclus] = pclus->eta() - sc.eta();
00224 clusterHitFractionSharedWithSeed.get()[iclus] = fractionOfSeed(pclus);
00225 if( _dogen && genmatch.isNonnull() ) {
00226 clusterDPhiToGen.get()[iclus] =
00227 TVector2::Phi_mpi_pi(pclus->phi() - genmatch->phi());
00228 clusterDEtaToGen.get()[iclus] = pclus->eta() - genmatch->eta();
00229 }
00230 clusterInMustache.get()[iclus] = (Int_t) MK::inMustache(theseed->eta(),
00231 theseed->phi(),
00232 pclus->energy(),
00233 pclus->eta(),
00234 pclus->phi());
00235 clusterInDynDPhi.get()[iclus] = (Int_t)
00236 MK::inDynamicDPhiWindow(PFLayer::ECAL_BARREL == pclus->layer(),
00237 theseed->phi(),
00238 pclus->energy(),
00239 pclus->eta(),
00240 pclus->phi());
00241 ++iclus;
00242 }
00243
00244 auto psclusend = sc.preshowerClustersEnd();
00245 size_t ipsclus = 0;
00246 edm::Ptr<reco::PFCluster> ppsclus;
00247 for( auto psclus = sc.preshowerClustersBegin(); psclus != psclusend;
00248 ++psclus ) {
00249 ppsclus = edm::Ptr<reco::PFCluster>(*psclus);
00250 psClusterRawEnergy.get()[ipsclus] = ppsclus->energy();
00251 psClusterEta.get()[ipsclus] = ppsclus->eta();
00252 psClusterPhi.get()[ipsclus] = ppsclus->phi();
00253 ++ipsclus;
00254 }
00255 _tree->Fill();
00256 }
00257
00258 PFSuperClusterTreeMaker::PFSuperClusterTreeMaker(const PSet& p) {
00259 _calib.reset(new PFEnergyCalibration());
00260 N_ECALClusters = 1;
00261 N_PSClusters = 1;
00262 _tree = _fs->make<TTree>("SuperClusterTree","Dump of all available SC info");
00263 _tree->Branch("N_ECALClusters",&N_ECALClusters,"N_ECALClusters/I");
00264 _tree->Branch("N_PSClusters",&N_PSClusters,"N_PSClusters/I");
00265 _tree->Branch("nVtx",&nVtx,"nVtx/I");
00266 _tree->Branch("scRawEnergy",&scRawEnergy,"scRawEnergy/F");
00267 _tree->Branch("scCalibratedEnergy",&scCalibratedEnergy,
00268 "scCalibratedEnergy/F");
00269 _tree->Branch("scPreshowerEnergy",&scPreshowerEnergy,"scPreshowerEnergy/F");
00270 _tree->Branch("scEta",&scEta,"scEta/F");
00271 _tree->Branch("scPhi",&scPhi,"scPhi/F");
00272 _tree->Branch("scR",&scR,"scR/F");
00273 _tree->Branch("scPhiWidth",&scPhiWidth,"scPhiWidth/F");
00274 _tree->Branch("scEtaWidth",&scEtaWidth,"scEtaWidth/F");
00275 _tree->Branch("scSeedRawEnergy",&scSeedRawEnergy,"scSeedRawEnergy/F");
00276 _tree->Branch("scSeedCalibratedEnergy",&scSeedCalibratedEnergy,
00277 "scSeedCalibratedEnergy/F");
00278 _tree->Branch("scSeedEta",&scSeedEta,"scSeedEta/F");
00279 _tree->Branch("scSeedPhi",&scSeedPhi,"scSeedPhi/F");
00280
00281 clusterRawEnergy.reset(new Float_t[1],array_deleter<Float_t>());
00282 _tree->Branch("clusterRawEnergy",clusterRawEnergy.get(),
00283 "clusterRawEnergy[N_ECALClusters]/F");
00284 clusterCalibEnergy.reset(new Float_t[1],array_deleter<Float_t>());
00285 _tree->Branch("clusterCalibEnergy",clusterCalibEnergy.get(),
00286 "clusterCalibEnergy[N_ECALClusters]/F");
00287 clusterEta.reset(new Float_t[1],array_deleter<Float_t>());
00288 _tree->Branch("clusterEta",clusterEta.get(),
00289 "clusterEta[N_ECALClusters]/F");
00290 clusterPhi.reset(new Float_t[1],array_deleter<Float_t>());
00291 _tree->Branch("clusterPhi",clusterPhi.get(),
00292 "clusterPhi[N_ECALClusters]/F");
00293 clusterDPhiToSeed.reset(new Float_t[1],array_deleter<Float_t>());
00294 _tree->Branch("clusterDPhiToSeed",clusterDPhiToSeed.get(),
00295 "clusterDPhiToSeed[N_ECALClusters]/F");
00296 clusterDEtaToSeed.reset(new Float_t[1],array_deleter<Float_t>());
00297 _tree->Branch("clusterDEtaToSeed",clusterDEtaToSeed.get(),
00298 "clusterDEtaToSeed[N_ECALClusters]/F");
00299 clusterDPhiToCentroid.reset(new Float_t[1],array_deleter<Float_t>());
00300 _tree->Branch("clusterDPhiToCentroid",clusterDPhiToCentroid.get(),
00301 "clusterDPhiToCentroid[N_ECALClusters]/F");
00302 clusterDEtaToCentroid.reset(new Float_t[1],array_deleter<Float_t>());
00303 _tree->Branch("clusterDEtaToCentroid",clusterDEtaToCentroid.get(),
00304 "clusterDEtaToCentroid[N_ECALClusters]/F");
00305 clusterHitFractionSharedWithSeed.reset(new Float_t[1],
00306 array_deleter<Float_t>());
00307 _tree->Branch("clusterHitFractionSharedWithSeed",
00308 clusterHitFractionSharedWithSeed.get(),
00309 "clusterHitFractionSharedWithSeed[N_ECALClusters]/F");
00310 clusterInMustache.reset(new Int_t[1],array_deleter<Int_t>());
00311 _tree->Branch("clusterInMustache",clusterInMustache.get(),
00312 "clusterInMustache[N_ECALClusters]/I");
00313 clusterInDynDPhi.reset(new Int_t[1],array_deleter<Int_t>());
00314 _tree->Branch("clusterInDynDPhi",clusterInDynDPhi.get(),
00315 "clusterInDynDPhi[N_ECALClusters]/I");
00316
00317 psClusterRawEnergy.reset(new Float_t[1],array_deleter<Float_t>());
00318 _tree->Branch("psClusterRawEnergy",psClusterRawEnergy.get(),
00319 "psClusterRawEnergy[N_PSClusters]/F");
00320 psClusterEta.reset(new Float_t[1],array_deleter<Float_t>());
00321 _tree->Branch("psClusterEta",psClusterEta.get(),
00322 "psClusterEta[N_PSClusters]/F");
00323 psClusterPhi.reset(new Float_t[1],array_deleter<Float_t>());
00324 _tree->Branch("psClusterPhi",psClusterPhi.get(),
00325 "psClusterPhi[N_PSClusters]/F");
00326
00327
00328 if( (_dogen = p.getUntrackedParameter<bool>("doGen",false)) ) {
00329 _geninput = p.getParameter<edm::InputTag>("genSrc");
00330 _tree->Branch("genEta",&genEta,"genEta/F");
00331 _tree->Branch("genPhi",&genPhi,"genPhi/F");
00332 _tree->Branch("genEnergy",&genEnergy,"genEnergy/F");
00333 _tree->Branch("genDRToCentroid",&genDRToCentroid,"genDRToCentroid/F");
00334 _tree->Branch("genDRToSeed",&genDRToSeed,"genDRToSeed/F");
00335
00336 clusterDPhiToGen.reset(new Float_t[1],array_deleter<Float_t>());
00337 _tree->Branch("clusterDPhiToGen",clusterDPhiToGen.get(),
00338 "clusterDPhiToGen[N_ECALClusters]/F");
00339 clusterDEtaToGen.reset(new Float_t[1],array_deleter<Float_t>());
00340 _tree->Branch("clusterDEtaToGen",clusterDEtaToGen.get(),
00341 "clusterDPhiToGen[N_ECALClusters]/F");
00342 }
00343 _vtxsrc = p.getParameter<edm::InputTag>("primaryVertices");
00344 _scInputEB = p.getParameter<edm::InputTag>("superClusterSrcEB");
00345 _scInputEE = p.getParameter<edm::InputTag>("superClusterSrcEE");
00346 }
00347
00348
00349 void PFSuperClusterTreeMaker::setTreeArraysForSize(const size_t N_ECAL,
00350 const size_t N_PS) {
00351 Float_t* cRE_new = new Float_t[N_ECAL];
00352 clusterRawEnergy.reset(cRE_new,array_deleter<Float_t>());
00353 _tree->GetBranch("clusterRawEnergy")->SetAddress(clusterRawEnergy.get());
00354 Float_t* cCE_new = new Float_t[N_ECAL];
00355 clusterCalibEnergy.reset(cCE_new,array_deleter<Float_t>());
00356 _tree->GetBranch("clusterCalibEnergy")->SetAddress(clusterCalibEnergy.get());
00357 Float_t* cEta_new = new Float_t[N_ECAL];
00358 clusterEta.reset(cEta_new,array_deleter<Float_t>());
00359 _tree->GetBranch("clusterEta")->SetAddress(clusterEta.get());
00360 Float_t* cPhi_new = new Float_t[N_ECAL];
00361 clusterPhi.reset(cPhi_new,array_deleter<Float_t>());
00362 _tree->GetBranch("clusterPhi")->SetAddress(clusterPhi.get());
00363 Float_t* cDPhiSeed_new = new Float_t[N_ECAL];
00364 clusterDPhiToSeed.reset(cDPhiSeed_new,array_deleter<Float_t>());
00365 _tree->GetBranch("clusterDPhiToSeed")->SetAddress(clusterDPhiToSeed.get());
00366 Float_t* cDEtaSeed_new = new Float_t[N_ECAL];
00367 clusterDEtaToSeed.reset(cDEtaSeed_new,array_deleter<Float_t>());
00368 _tree->GetBranch("clusterDEtaToSeed")->SetAddress(clusterDEtaToSeed.get());
00369 Float_t* cDPhiCntr_new = new Float_t[N_ECAL];
00370 clusterDPhiToCentroid.reset(cDPhiCntr_new,array_deleter<Float_t>());
00371 _tree->GetBranch("clusterDPhiToCentroid")->SetAddress(clusterDPhiToCentroid.get());
00372 Float_t* cDEtaCntr_new = new Float_t[N_ECAL];
00373 clusterDEtaToCentroid.reset(cDEtaCntr_new,array_deleter<Float_t>());
00374 _tree->GetBranch("clusterDEtaToCentroid")->SetAddress(clusterDEtaToCentroid.get());
00375 Float_t* cHitFracShared_new = new Float_t[N_ECAL];
00376 clusterHitFractionSharedWithSeed.reset(cHitFracShared_new,
00377 array_deleter<Float_t>());
00378 _tree->GetBranch("clusterHitFractionSharedWithSeed")->SetAddress(clusterHitFractionSharedWithSeed.get());
00379
00380 if( _dogen ) {
00381 Float_t* cDPhiGen_new = new Float_t[N_ECAL];
00382 clusterDPhiToGen.reset(cDPhiGen_new,array_deleter<Float_t>());
00383 _tree->GetBranch("clusterDPhiToGen")->SetAddress(clusterDPhiToGen.get());
00384 Float_t* cDEtaGen_new = new Float_t[N_ECAL];
00385 clusterDEtaToGen.reset(cDEtaGen_new,array_deleter<Float_t>());
00386 _tree->GetBranch("clusterDEtaToGen")->SetAddress(clusterDEtaToGen.get());
00387 }
00388 Int_t* cInMust_new = new Int_t[N_ECAL];
00389 clusterInMustache.reset(cInMust_new,array_deleter<Int_t>());
00390 _tree->GetBranch("clusterInMustache")->SetAddress(clusterInMustache.get());
00391 Int_t* cInDynDPhi_new = new Int_t[N_ECAL];
00392 clusterInDynDPhi.reset(cInDynDPhi_new,array_deleter<Int_t>());
00393 _tree->GetBranch("clusterInDynDPhi")->SetAddress(clusterInDynDPhi.get());
00394 Float_t* psRE_new = new Float_t[N_PS];
00395 psClusterRawEnergy.reset(psRE_new,array_deleter<Float_t>());
00396 _tree->GetBranch("psClusterRawEnergy")->SetAddress(psClusterRawEnergy.get());
00397 Float_t* psEta_new = new Float_t[N_PS];
00398 psClusterEta.reset(psEta_new,array_deleter<Float_t>());
00399 _tree->GetBranch("psClusterEta")->SetAddress(psClusterEta.get());
00400 Float_t* psPhi_new = new Float_t[N_PS];
00401 psClusterPhi.reset(psPhi_new,array_deleter<Float_t>());
00402 _tree->GetBranch("psClusterPhi")->SetAddress(psClusterPhi.get());
00403 }
00404
00405 #include "FWCore/Framework/interface/MakerMacros.h"
00406 DEFINE_FWK_MODULE(PFSuperClusterTreeMaker);