CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoEcal/EgammaClusterProducers/src/PFSuperClusterTreeMaker.cc

Go to the documentation of this file.
00001 //
00002 // Class: PFSuperClusterTreeMaker.cc
00003 //
00004 // Info: Processes a track into histograms of delta-phis and such
00005 //
00006 // Author: L. Gray (FNAL)
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       // now see if the clusters overlap in rechits
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   // the tree  
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); // minus 1 because of seed
00150   N_PSClusters = N_PS;
00151   reco::GenParticleRef genmatch;
00152   // gen information (if needed)
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           // oh hai, this is hard coded to photons for the time being
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   // supercluster information
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   // sc seed information
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   // loop over all clusters that aren't the seed
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   // loop over all preshower clusters 
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   // ecal cluster information
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   // preshower information
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);