First commit

RecoEgamma/EgammaPhotonProducers/interface/GEDPhotonCoreProducer.h

@@ -0,0 +1,49 @@
+#ifndef RecoEgamma_EgammaPhotonProducers_GEDPhotonCoreProducer_h
+#define RecoEgamma_EgammaPhotonProducers_GEDPhotonCoreProducer_h
+/** \class GEDPhotonCoreProducer
+ **  
+ **
+ **  $Id: $ 
+ **  $Date: $ 
+ **  $Revision: $
+ **  \author Nancy Marinelli, U. of Notre Dame, US
+ **
+ ***/
+
+#include "FWCore/Framework/interface/EDProducer.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "DataFormats/Common/interface/Handle.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "DataFormats/EgammaCandidates/interface/PhotonCore.h"
+//#include "RecoEgamma/EgammaTools/interface/HoECalculator.h"
+//#include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
+//#include "RecoEgamma/PhotonIdentification/interface/PhotonIsolationCalculator.h"
+
+// GEDPhotonCoreProducer inherits from EDProducer, so it can be a module:
+class GEDPhotonCoreProducer : public edm::EDProducer {
+
+ public:
+
+  GEDPhotonCoreProducer (const edm::ParameterSet& ps);
+  ~GEDPhotonCoreProducer();
+
+  virtual void produce(edm::Event& evt, const edm::EventSetup& es);
+
+ private:
+
+  void createSingleLegConversions( reco::CaloClusterPtr, std::vector<reco::TrackRef>, std::vector<float>,  reco::ConversionCollection &oneLegConversions  );
+
+  std::string GEDPhotonCoreCollection_;
+  std::string PFConversionCollection_;
+  edm::InputTag pfEgammaCandidates_;
+
+  double minSCEt_;
+  bool validConversions_;
+  std::string pixelSeedProducer_;
+  edm::ParameterSet conf_;
+  bool validPixelSeeds_;
+
+};
+#endif

RecoEgamma/EgammaPhotonProducers/interface/GEDPhotonProducer.h

@@ -0,0 +1,120 @@
+#ifndef RecoEgamma_EgammaPhotonProducers_GEDPhotonProducer_h
+#define RecoEgamma_EgammaPhotonProducers_GEDPhotonProducer_h
+/** \class GEDPhotonProducer
+ **  
+ **
+ **  $Id: GEDPhotonProducer.h $ 
+ **  $Date:  $ 
+ **  $Revision: $
+ **  \author Nancy Marinelli, U. of Notre Dame, US
+ **
+ ***/
+
+#include "FWCore/Framework/interface/EDProducer.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "DataFormats/Common/interface/Handle.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
+#include "Geometry/CaloTopology/interface/CaloTopology.h"
+#include "DataFormats/EgammaCandidates/interface/Conversion.h"
+#include "DataFormats/EgammaReco/interface/BasicCluster.h"
+#include "DataFormats/EgammaReco/interface/SuperCluster.h"
+#include "DataFormats/EgammaCandidates/interface/Photon.h"
+#include "DataFormats/EgammaCandidates/interface/PhotonCore.h"
+#include "DataFormats/EgammaReco/interface/BasicClusterShapeAssociation.h"
+#include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
+#include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
+#include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h"
+#include "RecoCaloTools/MetaCollections/interface/CaloRecHitMetaCollections.h"
+#include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
+#include "RecoEgamma/PhotonIdentification/interface/PhotonIsolationCalculator.h"
+#include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
+#include "RecoEgamma/PhotonIdentification/interface/PhotonMIPHaloTagger.h"
+#include "RecoEcal/EgammaCoreTools/interface/EcalClusterFunctionFactory.h"
+#include "RecoEcal/EgammaCoreTools/interface/EcalClusterFunctionBaseClass.h" 
+#include "CondFormats/EcalObjects/interface/EcalFunctionParameters.h" 
+#include "RecoEgamma/EgammaPhotonAlgos/interface/PhotonEnergyCorrector.h"
+
+// GEDPhotonProducer inherits from EDProducer, so it can be a module:
+class GEDPhotonProducer : public edm::EDProducer {
+
+ public:
+
+  GEDPhotonProducer (const edm::ParameterSet& ps);
+  ~GEDPhotonProducer();
+
+  virtual void beginRun (edm::Run const& r, edm::EventSetup const & es) override final;
+  virtual void endRun(edm::Run const&,  edm::EventSetup const&) override final;
+  virtual void produce(edm::Event& evt, const edm::EventSetup& es);
+
+ private:
+
+  void fillPhotonCollection(edm::Event& evt,
+			    edm::EventSetup const & es,
+                            const edm::Handle<reco::PhotonCoreCollection> & photonCoreHandle,
+                            const CaloTopology *topology,
+			    const EcalRecHitCollection* ecalBarrelHits,
+			    const EcalRecHitCollection* ecalEndcapHits,
+			    const edm::Handle<CaloTowerCollection> & hcalTowersHandle,
+			    //math::XYZPoint & vtx,
+			    reco::VertexCollection& pvVertices,
+			    reco::PhotonCollection & outputCollection,
+			    int& iSC,
+			    const EcalSeverityLevelAlgo * sevLv);
+
+  // std::string PhotonCoreCollection_;
+  std::string PhotonCollection_;
+  edm::InputTag photonCoreProducer_;
+  edm::InputTag barrelEcalHits_;
+  edm::InputTag endcapEcalHits_;
+
+  edm::InputTag hcalTowers_;
+
+  std::string conversionProducer_;
+  std::string conversionCollection_;
+
+  //AA
+  //Flags and severities to be excluded from calculations
+
+  std::vector<int> flagsexclEB_;
+  std::vector<int> flagsexclEE_;
+  std::vector<int> severitiesexclEB_;
+  std::vector<int> severitiesexclEE_;
+
+
+  double hOverEConeSize_;
+  double maxHOverE_;
+  double minSCEt_;
+  double highEt_;
+  double  minR9Barrel_;
+  double  minR9Endcap_;
+  bool   runMIPTagger_;
+
+  bool validConversions_;
+  std::string pixelSeedProducer_;
+  std::string vertexProducer_;
+  bool usePrimaryVertex_;
+  edm::ParameterSet conf_;
+
+  PositionCalc posCalculator_;
+
+  edm::ESHandle<CaloGeometry> theCaloGeom_;
+  edm::ESHandle<CaloTopology> theCaloTopo_;
+
+  bool validPixelSeeds_;
+  PhotonIsolationCalculator* thePhotonIsolationCalculator_;
+
+  //MIP
+  PhotonMIPHaloTagger* thePhotonMIPHaloTagger_;
+
+  std::vector<double>  preselCutValuesBarrel_; 
+  std::vector<double>  preselCutValuesEndcap_; 
+
+  EcalClusterFunctionBaseClass* energyCorrectionF;
+  PhotonEnergyCorrector* thePhotonEnergyCorrector_;
+  std::string  candidateP4type_;
+
+};
+#endif

RecoEgamma/EgammaPhotonProducers/src/GEDPhotonCoreProducer.cc

@@ -0,0 +1,209 @@
+#include <iostream>
+#include <vector>
+#include <memory>
+/** \class GEDPhotonCoreProducer
+ **  
+ **
+ **  $Id: $ 
+ **  $Date: $ 
+ **  $Revision: $
+ **  \author Nancy Marinelli, U. of Notre Dame, US
+ **
+ ***/
+
+// Framework
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "FWCore/Utilities/interface/Exception.h"
+#include "DataFormats/EgammaReco/interface/ClusterShape.h"
+#include "DataFormats/EgammaCandidates/interface/Photon.h"
+#include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
+#include "DataFormats/EgammaCandidates/interface/Conversion.h"
+#include "RecoEgamma/EgammaPhotonProducers/interface/GEDPhotonCoreProducer.h"
+#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
+#include "DataFormats/ParticleFlowCandidate/interface/PFCandidateEGammaExtra.h"
+#include "DataFormats/ParticleFlowCandidate/interface/PFCandidateEGammaExtraFwd.h"
+#include "DataFormats/VertexReco/interface/Vertex.h"
+#include "DataFormats/VertexReco/interface/VertexFwd.h"
+#include "DataFormats/TrackReco/interface/Track.h"
+#include "DataFormats/TrackReco/interface/TrackFwd.h"
+#include "DataFormats/EgammaReco/interface/ElectronSeed.h"
+#include "DataFormats/Math/interface/Vector3D.h"
+#include "DataFormats/Math/interface/LorentzVector.h"
+#include "DataFormats/GeometryVector/interface/GlobalVector.h"
+#include "DataFormats/Common/interface/RefToPtr.h"
+#include <Math/VectorUtil.h>
+#include <vector>
+#include "TLorentzVector.h"
+#include "TMath.h"
+
+GEDPhotonCoreProducer::GEDPhotonCoreProducer(const edm::ParameterSet& config) : 
+  conf_(config)
+
+{
+
+  // use onfiguration file to setup input/output collection names
+  pfEgammaCandidates_      = conf_.getParameter<edm::InputTag>("pfEgammaCandidates");
+  GEDPhotonCoreCollection_ = conf_.getParameter<std::string>("gedPhotonCoreCollection");
+
+  // Register the product
+  produces< reco::PhotonCoreCollection >(GEDPhotonCoreCollection_);
+  produces<reco::ConversionCollection>(PFConversionCollection_);
+
+}
+
+GEDPhotonCoreProducer::~GEDPhotonCoreProducer() {}
+
+
+
+
+void GEDPhotonCoreProducer::produce(edm::Event &theEvent, const edm::EventSetup& theEventSetup) {
+
+
+  using namespace edm;
+  //  nEvt_++;
+
+  reco::PhotonCoreCollection outputPhotonCoreCollection;
+  std::auto_ptr< reco::PhotonCoreCollection > outputPhotonCoreCollection_p(new reco::PhotonCoreCollection);
+
+  reco::ConversionCollection outputOneLegConversionCollection;
+  std::auto_ptr<reco::ConversionCollection> SingleLeg_p(new reco::ConversionCollection(outputOneLegConversionCollection));  
+
+  // Get the  PF refined cluster  collection
+  Handle<reco::PFCandidateCollection> pfCandidateHandle;
+  theEvent.getByLabel(pfEgammaCandidates_,pfCandidateHandle);
+  if (!pfCandidateHandle.isValid()) {
+    edm::LogError("GEDPhotonCoreProducer") << "Error! Can't get the product "<<pfEgammaCandidates_.label();
+  }
+
+
+  //  std::cout <<  "  GEDPhotonCoreProducer::produce input PFcandidate size " <<   pfCandidateHandle->size() << std::endl;
+
+
+  // Loop over PF candidates and get only photons
+
+  for(unsigned int lCand=0; lCand < pfCandidateHandle->size(); lCand++) {
+    reco::PFCandidateRef candRef (reco::PFCandidateRef(pfCandidateHandle,lCand));
+    if(candRef->particleId()!=reco::PFCandidate::gamma) continue;
+
+    // Retrieve stuff from the pfPhoton
+    reco::PFCandidateEGammaExtraRef pfPhoRef =  candRef->egammaExtraRef();
+    reco::SuperClusterRef  refinedSC= pfPhoRef->superClusterRef();
+    reco::SuperClusterRef  boxSC= pfPhoRef->superClusterBoxRef();
+    const reco::ConversionRefVector & doubleLegConv = pfPhoRef->conversionRef();
+    reco::CaloClusterPtr refinedSCPtr= edm::refToPtr(refinedSC);
+
+    //    std::cout << "newCandidate  doubleLegConv="<<doubleLegConv.size()<< std::endl;
+    //std::cout << "newCandidate  singleLegConv="<<  pfPhoRef->singleLegConvTrackRef().size()<< std::endl;
+
+    //////////
+    reco::PhotonCore newCandidate;
+    newCandidate.setPFlowPhoton(true);
+    newCandidate.setStandardPhoton(false);
+    newCandidate.setSuperCluster(refinedSC);
+    newCandidate.setPflowSuperCluster(boxSC);
+    // fill conversion infos
+
+
+    for(unsigned int lConv=0; lConv < doubleLegConv.size(); lConv++) {
+      newCandidate.addConversion(doubleLegConv[lConv]);
+    } 
+
+    //    std::cout << "newCandidate pf refined SC energy="<< newCandidate.superCluster()->energy()<<std::endl;
+    //std::cout << "newCandidate pf SC energy="<< newCandidate.pfSuperCluster()->energy()<<std::endl;
+    //std::cout << "newCandidate  nconv2leg="<<newCandidate.conversions().size()<< std::endl;
+
+
+
+    //
+    createSingleLegConversions(refinedSCPtr,  pfPhoRef->singleLegConvTrackRef(), pfPhoRef->singleLegConvMva(), outputOneLegConversionCollection); 
+    outputPhotonCoreCollection.push_back(newCandidate);
+  }
+
+
+  const edm::OrphanHandle<reco::ConversionCollection> singleLegConvOrhpHandle = theEvent.put(SingleLeg_p,PFConversionCollection_);
+
+  //  std::cout <<  "  GEDPhotonCoreProducer::produce orphanHandle to single legs " <<  singleLegConvOrhpHandle->size() << std::endl;
+  //std::cout <<  "  GEDPhotonCoreProducer::produce photon size " <<  outputPhotonCoreCollection.size() << std::endl;
+
+
+
+  int ipho=0;
+  for (reco::PhotonCoreCollection::iterator gamIter=outputPhotonCoreCollection.begin(); gamIter != outputPhotonCoreCollection.end(); ++gamIter){
+
+    for( unsigned int icp = 0;  icp <  singleLegConvOrhpHandle->size(); icp++) {
+      const reco::ConversionRef cpRef(reco::ConversionRef(singleLegConvOrhpHandle,icp));
+      if ( !cpRef->caloCluster().size()) continue; 
+      if (!( gamIter->superCluster().id() == cpRef->caloCluster()[0].id() && gamIter->superCluster().key() == cpRef->caloCluster()[0].key() )) continue; 
+      gamIter->addOneLegConversion(cpRef);     
+    }
+    // debug
+    //    std::cout << "PhotonCoreCollection i="<<ipho<<" pf refined SC energy="<<gamIter->superCluster()->energy()<<std::endl;
+    //std::cout << "PhotonCoreCollection i="<<ipho<<" pf SC energy="<<gamIter->pfSuperCluster()->energy()<<std::endl;
+    //std::cout << "PhotonCoreCollection i="<<ipho<<" nconv2leg="<<gamIter->conversions().size()<<" nconv1leg="<<gamIter->conversionsOneLeg().size()<<std::endl;
+    ipho++;
+  }
+
+
+
+  // put the product in the event
+  //  edm::LogInfo("GEDPhotonCoreProducer") << " Put in the event " << iSC << " Photon Candidates \n";
+  outputPhotonCoreCollection_p->assign(outputPhotonCoreCollection.begin(),outputPhotonCoreCollection.end());
+  theEvent.put( outputPhotonCoreCollection_p, GEDPhotonCoreCollection_);
+
+
+
+
+} 
+
+
+void GEDPhotonCoreProducer::createSingleLegConversions(reco::CaloClusterPtr sc,  std::vector<reco::TrackRef>  conv, std::vector<float> mva,  reco::ConversionCollection &oneLegConversions) {
+  // this method translates the single track into the Conversion Data Format
+
+  math::Error<3>::type error;
+  for (unsigned int itk=0; itk<conv.size(); itk++){
+    const reco::Vertex  * convVtx = new reco::Vertex(conv[itk]->innerPosition(), error);
+    std::vector<reco::TrackRef> OneLegConvVector;
+    OneLegConvVector.push_back(conv[itk]);
+    std::vector< float > OneLegMvaVector;
+    OneLegMvaVector.push_back(mva[itk]);
+    std::vector<reco::CaloClusterPtr> dummymatchingBC;
+    reco::CaloClusterPtrVector scPtrVec;
+    scPtrVec.push_back(sc);
+
+
+    std::vector<math::XYZPointF>trackPositionAtEcalVec;
+    std::vector<math::XYZPointF>innPointVec;
+    std::vector<math::XYZVectorF>trackPinVec;
+    std::vector<math::XYZVectorF>trackPoutVec;
+    math::XYZPointF trackPositionAtEcal(conv[itk]->outerPosition().X(), conv[itk]->outerPosition().Y(), conv[itk]->outerPosition().Z());
+    trackPositionAtEcalVec.push_back(trackPositionAtEcal);
+
+    math::XYZPointF innPoint(conv[itk]->innerPosition().X(), conv[itk]->innerPosition().Y(), conv[itk]->innerPosition().Z());
+    innPointVec.push_back(innPoint);
+
+    math::XYZVectorF trackPin(conv[itk]->innerMomentum().X(), conv[itk]->innerMomentum().Y(), conv[itk]->innerMomentum().Z());
+    trackPinVec.push_back(trackPin);
+
+    math::XYZVectorF trackPout(conv[itk]->outerMomentum().X(), conv[itk]->outerMomentum().Y(), conv[itk]->outerMomentum().Z());
+    trackPoutVec.push_back( trackPout );
+
+    float DCA = conv[itk]->d0() ;
+    reco::Conversion singleLegConvCandidate(scPtrVec, 
+					OneLegConvVector,
+					trackPositionAtEcalVec,
+					*convVtx,
+					dummymatchingBC,
+					DCA,
+					innPointVec,
+					trackPinVec,
+					trackPoutVec,
+					mva[itk],			  
+					reco::Conversion::pflow);
+    singleLegConvCandidate.setOneLegMVA(OneLegMvaVector); 
+    oneLegConversions.push_back(singleLegConvCandidate);
+
+  } 
+
+
+
+}