#include <PhotonProducer.h>

Inheritance diagram for PhotonProducer:

Public Member Functions
virtual void	beginRun (edm::Run &r, edm::EventSetup const &es)

virtual void	endRun (edm::Run &, edm::EventSetup const &)

	PhotonProducer (const edm::ParameterSet &ps)

virtual void	produce (edm::Event &evt, const edm::EventSetup &es)

	~PhotonProducer ()

Public Member Functions inherited from edm::EDProducer
	EDProducer ()

virtual	~EDProducer ()

Public Member Functions inherited from edm::ProducerBase
	ProducerBase ()

void	registerProducts (ProducerBase , ProductRegistry , ModuleDescription const &)

boost::function< void(const BranchDescription &)>	registrationCallback () const
	used by the fwk to register list of products More...

virtual	~ProducerBase ()

Private Member Functions
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::PhotonCollection &outputCollection, int &iSC)

Private Attributes
edm::InputTag	barrelEcalHits_

edm::ParameterSet	conf_

std::string	conversionCollection_

std::string	conversionProducer_

edm::InputTag	endcapEcalHits_

EcalClusterFunctionBaseClass *	energyCorrectionF

edm::InputTag	hcalTowers_

double	highEt_

double	hOverEConeSize_

double	maxHOverE_

double	minR9Barrel_

double	minR9Endcap_

double	minSCEt_

std::string	PhotonCollection_

edm::InputTag	photonCoreProducer_

std::string	pixelSeedProducer_

PositionCalc	posCalculator_

std::vector< double >	preselCutValuesBarrel_

std::vector< double >	preselCutValuesEndcap_

edm::ESHandle< CaloGeometry >	theCaloGeom_

edm::ESHandle< CaloTopology >	theCaloTopo_

PhotonIsolationCalculator *	thePhotonIsolationCalculator_

bool	usePrimaryVertex_

bool	validConversions_

bool	validPixelSeeds_

std::string	vertexProducer_

Additional Inherited Members
Public Types inherited from edm::EDProducer
typedef EDProducer	ModuleType

typedef WorkerT< EDProducer >	WorkerType

Public Types inherited from edm::ProducerBase
typedef ProductRegistryHelper::TypeLabelList	TypeLabelList

Static Public Member Functions inherited from edm::EDProducer
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

Protected Member Functions inherited from edm::EDProducer
CurrentProcessingContext const *	currentContext () const

Protected Member Functions inherited from edm::ProducerBase
template<class TProducer , class TMethod >
void	callWhenNewProductsRegistered (TProducer *iProd, TMethod iMethod)

Detailed Description

Id:: PhotonProducer.h,v 1.36 2009/05/29 18:11:15 nancy Exp

Date:: 2009/05/29 18:11:15

Revision:: 1.36

Author: Nancy Marinelli, U. of Notre Dame, US

Definition at line 38 of file PhotonProducer.h.

Constructor & Destructor Documentation

PhotonProducer::PhotonProducer ( const edm::ParameterSet & ps )

Definition at line 29 of file PhotonProducer.cc.

References barrelEcalHits_, conf_, endcapEcalHits_, energyCorrectionF, reco::get(), edm::ParameterSet::getParameter(), hcalTowers_, highEt_, hOverEConeSize_, minR9Barrel_, minR9Endcap_, PhotonCollection_, photonCoreProducer_, posCalculator_, preselCutValuesBarrel_, preselCutValuesEndcap_, usePrimaryVertex_, and vertexProducer_.

                                                             : 
   conf_(config)
 {
 
   // use onfiguration file to setup input/output collection names
 
   photonCoreProducer_   = conf_.getParameter<edm::InputTag>("photonCoreProducer");
   barrelEcalHits_   = conf_.getParameter<edm::InputTag>("barrelEcalHits");
   endcapEcalHits_   = conf_.getParameter<edm::InputTag>("endcapEcalHits");
   vertexProducer_   = conf_.getParameter<std::string>("primaryVertexProducer");
   hcalTowers_ = conf_.getParameter<edm::InputTag>("hcalTowers");
   hOverEConeSize_   = conf_.getParameter<double>("hOverEConeSize");
   highEt_        = conf_.getParameter<double>("highEt");
   // R9 value to decide converted/unconverted
   minR9Barrel_        = conf_.getParameter<double>("minR9Barrel");
   minR9Endcap_        = conf_.getParameter<double>("minR9Endcap");
   usePrimaryVertex_ = conf_.getParameter<bool>("usePrimaryVertex");
  
   edm::ParameterSet posCalcParameters = 
     config.getParameter<edm::ParameterSet>("posCalcParameters");
   posCalculator_ = PositionCalc(posCalcParameters);
 
   // Parameters for the position calculation:
   //  std::map<std::string,double> providedParameters;
   // providedParameters.insert(std::make_pair("LogWeighted",conf_.getParameter<bool>("posCalc_logweight")));
   //providedParameters.insert(std::make_pair("T0_barl",conf_.getParameter<double>("posCalc_t0_barl")));
   //providedParameters.insert(std::make_pair("T0_endc",conf_.getParameter<double>("posCalc_t0_endc")));
   //providedParameters.insert(std::make_pair("T0_endcPresh",conf_.getParameter<double>("posCalc_t0_endcPresh")));
   //providedParameters.insert(std::make_pair("W0",conf_.getParameter<double>("posCalc_w0")));
   //providedParameters.insert(std::make_pair("X0",conf_.getParameter<double>("posCalc_x0")));
   //posCalculator_ = PositionCalc(providedParameters);
   // cut values for pre-selection
   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("minSCEtBarrel")); 
   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("maxHoverEBarrel")); 
   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("ecalRecHitSumEtOffsetBarrel")); 
   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("ecalRecHitSumEtSlopeBarrel")); 
   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("hcalTowerSumEtOffsetBarrel"));
   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("hcalTowerSumEtSlopeBarrel"));
   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("nTrackSolidConeBarrel"));
   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("nTrackHollowConeBarrel"));     
   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("trackPtSumSolidConeBarrel"));     
   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("trackPtSumHollowConeBarrel"));     
   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("sigmaIetaIetaCutBarrel"));     
   //  
   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("minSCEtEndcap")); 
   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("maxHoverEEndcap")); 
   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("ecalRecHitSumEtOffsetEndcap")); 
   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("ecalRecHitSumEtSlopeEndcap")); 
   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("hcalTowerSumEtOffsetEndcap"));
   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("hcalTowerSumEtSlopeEndcap"));
   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("nTrackSolidConeEndcap"));
   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("nTrackHollowConeEndcap"));     
   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("trackPtSumSolidConeEndcap"));     
   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("trackPtSumHollowConeEndcap"));     
   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("sigmaIetaIetaCutEndcap"));     
   //
   energyCorrectionF = EcalClusterFunctionFactory::get()->create("EcalClusterEnergyCorrection", conf_);
 
   // Register the product
   produces< reco::PhotonCollection >(PhotonCollection_);
 
 }

PhotonProducer::~PhotonProducer ( )

Definition at line 92 of file PhotonProducer.cc.

References energyCorrectionF.

                                 {
 
   delete energyCorrectionF;
 
 }

Member Function Documentation

void PhotonProducer::beginRun	(	edm::Run &	r,
		edm::EventSetup const &	es
	)

virtual

Reimplemented from edm::EDProducer.

Definition at line 100 of file PhotonProducer.cc.

References conf_, edm::ParameterSet::getParameter(), PhotonIsolationCalculator::setup(), and thePhotonIsolationCalculator_.

                                                                               {
 
     thePhotonIsolationCalculator_ = new PhotonIsolationCalculator();
     edm::ParameterSet isolationSumsCalculatorSet = conf_.getParameter<edm::ParameterSet>("isolationSumsCalculatorSet"); 
     thePhotonIsolationCalculator_->setup(isolationSumsCalculatorSet);
 
 }

void PhotonProducer::endRun	(	edm::Run &	r,
		edm::EventSetup const &	theEventSetup
	)

virtual

Reimplemented from edm::EDProducer.

Definition at line 108 of file PhotonProducer.cc.

References thePhotonIsolationCalculator_.

                                                                             {
 
   delete thePhotonIsolationCalculator_;
 
 }

void PhotonProducer::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::PhotonCollection &	outputCollection,
		int &	iSC
	)

private

fill shower shape block

Pre-selection loose isolation cuts

Definition at line 211 of file PhotonProducer.cc.

Referenced by produce().

                                                                                {
   
   const CaloGeometry* geometry = theCaloGeom_.product();
   const CaloSubdetectorGeometry* subDetGeometry =0 ;
   const CaloSubdetectorGeometry* geometryES = theCaloGeom_->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
   const EcalRecHitCollection* hits = 0 ;
   std::vector<double> preselCutValues;
   float minR9=0;
 
 
 
   for(unsigned int lSC=0; lSC < photonCoreHandle->size(); lSC++) {
 
     reco::PhotonCoreRef coreRef(reco::PhotonCoreRef(photonCoreHandle, lSC));
     reco::SuperClusterRef scRef=coreRef->superCluster();
     const reco::SuperCluster* pClus=&(*scRef);
     iSC++;
 
     int subdet = scRef->seed()->hitsAndFractions()[0].first.subdetId();
     subDetGeometry =  theCaloGeom_->getSubdetectorGeometry(DetId::Ecal, subdet);
 
     if (subdet==EcalBarrel) 
       { 
     preselCutValues = preselCutValuesBarrel_;
         minR9=minR9Barrel_;
         hits=  ecalBarrelHits;
       }
     else if  (subdet==EcalEndcap) 
       { 
     preselCutValues = preselCutValuesEndcap_;
         minR9=minR9Endcap_;
     hits=  ecalEndcapHits;
       }
     else
       { edm::LogWarning("")<<"PhotonProducer: do not know if it is a barrel or endcap SuperCluster" ; }
     
         
     // SC energy preselection
     if (scRef->energy()/cosh(scRef->eta()) <= preselCutValues[0] ) continue;
     // calculate HoE
     const CaloTowerCollection* hcalTowersColl = hcalTowersHandle.product();
     EgammaTowerIsolation towerIso1(hOverEConeSize_,0.,0.,1,hcalTowersColl) ;  
     EgammaTowerIsolation towerIso2(hOverEConeSize_,0.,0.,2,hcalTowersColl) ;  
     double HoE1=towerIso1.getTowerESum(&(*scRef))/scRef->energy();
     double HoE2=towerIso2.getTowerESum(&(*scRef))/scRef->energy(); 
     //  std::cout << " PhotonProducer " << HoE1  << "  HoE2 " << HoE2 << std::endl;
     // std::cout << " PhotonProducer calcualtion of HoE1 " << HoE1  << "  HoE2 " << HoE2 << std::endl;
 
     
     // recalculate position of seed BasicCluster taking shower depth for unconverted photon
     math::XYZPoint unconvPos = posCalculator_.Calculate_Location(scRef->seed()->hitsAndFractions(),hits,subDetGeometry,geometryES);
 
 
     float maxXtal =   EcalClusterTools::eMax( *(scRef->seed()), &(*hits) );
     float e1x5    =   EcalClusterTools::e1x5(  *(scRef->seed()), &(*hits), &(*topology)); 
     float e2x5    =   EcalClusterTools::e2x5Max(  *(scRef->seed()), &(*hits), &(*topology)); 
     float e3x3    =   EcalClusterTools::e3x3(  *(scRef->seed()), &(*hits), &(*topology)); 
     float e5x5    =   EcalClusterTools::e5x5( *(scRef->seed()), &(*hits), &(*topology)); 
     std::vector<float> cov =  EcalClusterTools::covariances( *(scRef->seed()), &(*hits), &(*topology), geometry); 
     float sigmaEtaEta = sqrt(cov[0]);
     std::vector<float> locCov =  EcalClusterTools::localCovariances( *(scRef->seed()), &(*hits), &(*topology)); 
     float sigmaIetaIeta = sqrt(locCov[0]);
 
     float r9 =e3x3/(scRef->rawEnergy());
     // compute position of ECAL shower
     math::XYZPoint caloPosition;
     double photonEnergy=0;
     if (r9>minR9) {
       caloPosition = unconvPos;
       // f(eta) correction to e5x5
       double deltaE = energyCorrectionF->getValue(*pClus, 1);
       if (subdet==EcalBarrel) e5x5 = e5x5 * (1.0 +  deltaE/scRef->rawEnergy() );
       photonEnergy=  e5x5    + scRef->preshowerEnergy() ;
     } else {
       caloPosition = scRef->position();
       photonEnergy=scRef->energy();
     }
     
 
     //std::cout << " dete " << subdet << " r9 " << r9 <<  " uncorrected e5x5 " <<  e5x5uncor << " corrected e5x5 " << e5x5 << " photon energy " << photonEnergy << std::endl;
     
     // compute momentum vector of photon from primary vertex and cluster position
     math::XYZVector direction = caloPosition - vtx;
     math::XYZVector momentum = direction.unit() * photonEnergy ;
 
     // Create candidate
     const reco::Particle::LorentzVector  p4(momentum.x(), momentum.y(), momentum.z(), photonEnergy );
     reco::Photon newCandidate(p4, caloPosition, coreRef, vtx);
 
     // Calculate fiducial flags and isolation variable. Blocked are filled from the isolationCalculator
     reco::Photon::FiducialFlags fiducialFlags;
     reco::Photon::IsolationVariables isolVarR03, isolVarR04;
     thePhotonIsolationCalculator_-> calculate ( &newCandidate,evt,es,fiducialFlags,isolVarR04, isolVarR03);
     newCandidate.setFiducialVolumeFlags( fiducialFlags );
     newCandidate.setIsolationVariables(isolVarR04, isolVarR03 );
   
     reco::Photon::ShowerShape  showerShape;
     showerShape.e1x5= e1x5;
     showerShape.e2x5= e2x5;
     showerShape.e3x3= e3x3;
     showerShape.e5x5= e5x5;
     showerShape.maxEnergyXtal =  maxXtal;
     showerShape.sigmaEtaEta =    sigmaEtaEta;
     showerShape.sigmaIetaIeta =  sigmaIetaIeta;
     showerShape.hcalDepth1OverEcal = HoE1;
     showerShape.hcalDepth2OverEcal = HoE2;
     newCandidate.setShowerShapeVariables ( showerShape ); 
 
     bool isLooseEM=true;
     if ( newCandidate.pt() < highEt_) { 
       if ( newCandidate.hadronicOverEm()                   >= preselCutValues[1] )                                            isLooseEM=false;
       if ( newCandidate.ecalRecHitSumEtConeDR04()          > preselCutValues[2]+ preselCutValues[3]*newCandidate.pt() )       isLooseEM=false;
       if ( newCandidate.hcalTowerSumEtConeDR04()           > preselCutValues[4]+ preselCutValues[5]*newCandidate.pt() )       isLooseEM=false;
       if ( newCandidate.nTrkSolidConeDR04()                > int(preselCutValues[6]) )                                        isLooseEM=false;
       if ( newCandidate.nTrkHollowConeDR04()               > int(preselCutValues[7]) )                                        isLooseEM=false;
       if ( newCandidate.trkSumPtSolidConeDR04()            > preselCutValues[8] )                                             isLooseEM=false;
       if ( newCandidate.trkSumPtHollowConeDR04()           > preselCutValues[9] )                                             isLooseEM=false;
       if ( newCandidate.sigmaIetaIeta()                    > preselCutValues[10] )                                            isLooseEM=false;
     } 
     
     
     if ( isLooseEM)  
       outputPhotonCollection.push_back(newCandidate);
       
         
   }
 }

void PhotonProducer::produce	(	edm::Event &	evt,
		const edm::EventSetup &	es
	)

virtual

Implements edm::EDProducer.

Definition at line 117 of file PhotonProducer.cc.

References barrelEcalHits_, endcapEcalHits_, energyCorrectionF, fillPhotonCollection(), edm::EventSetup::get(), edm::Event::getByLabel(), hcalTowers_, EcalClusterFunctionBaseClass::init(), edm::InputTag::label(), PhotonCollection_, photonCoreProducer_, edm::ESHandle< class >::product(), edm::Event::put(), theCaloGeom_, theCaloTopo_, usePrimaryVertex_, GoodVertex_cfg::vertexCollection, and vertexProducer_.

                                                                                    {
 
   using namespace edm;
   //  nEvt_++;
 
   reco::PhotonCollection outputPhotonCollection;
   std::auto_ptr< reco::PhotonCollection > outputPhotonCollection_p(new reco::PhotonCollection);
 
 
   // Get the PhotonCore collection
   bool validPhotonCoreHandle=true;
   Handle<reco::PhotonCoreCollection> photonCoreHandle;
   theEvent.getByLabel(photonCoreProducer_,photonCoreHandle);
   if (!photonCoreHandle.isValid()) {
     edm::LogError("PhotonProducer") << "Error! Can't get the product "<<photonCoreProducer_.label();
     validPhotonCoreHandle=false;
   }
 
  // Get EcalRecHits
   bool validEcalRecHits=true;
   Handle<EcalRecHitCollection> barrelHitHandle;
   EcalRecHitCollection barrelRecHits;
   theEvent.getByLabel(barrelEcalHits_, barrelHitHandle);
   if (!barrelHitHandle.isValid()) {
     edm::LogError("PhotonProducer") << "Error! Can't get the product "<<barrelEcalHits_.label();
     validEcalRecHits=false; 
   }
   if (  validEcalRecHits)  barrelRecHits = *(barrelHitHandle.product());
 
   
   Handle<EcalRecHitCollection> endcapHitHandle;
   theEvent.getByLabel(endcapEcalHits_, endcapHitHandle);
   EcalRecHitCollection endcapRecHits;
   if (!endcapHitHandle.isValid()) {
     edm::LogError("PhotonProducer") << "Error! Can't get the product "<<endcapEcalHits_.label();
     validEcalRecHits=false; 
   }
   if( validEcalRecHits) endcapRecHits = *(endcapHitHandle.product());
 
 
 // get Hcal towers collection 
   Handle<CaloTowerCollection> hcalTowersHandle;
   theEvent.getByLabel(hcalTowers_, hcalTowersHandle);
 
 
   // get the geometry from the event setup:
   theEventSetup.get<CaloGeometryRecord>().get(theCaloGeom_);
 
   //
   // update energy correction function
   energyCorrectionF->init(theEventSetup);  
 
   edm::ESHandle<CaloTopology> pTopology;
   theEventSetup.get<CaloTopologyRecord>().get(theCaloTopo_);
   const CaloTopology *topology = theCaloTopo_.product();
 
   // Get the primary event vertex
   Handle<reco::VertexCollection> vertexHandle;
   reco::VertexCollection vertexCollection;
   bool validVertex=true;
   if ( usePrimaryVertex_ ) {
     theEvent.getByLabel(vertexProducer_, vertexHandle);
     if (!vertexHandle.isValid()) {
       edm::LogWarning("PhotonProducer") << "Error! Can't get the product primary Vertex Collection "<< "\n";
       validVertex=false;
     }
     if (validVertex) vertexCollection = *(vertexHandle.product());
   }
   math::XYZPoint vtx(0.,0.,0.);
   if (vertexCollection.size()>0) vtx = vertexCollection.begin()->position();
 
 
   int iSC=0; // index in photon collection
   // Loop over barrel and endcap SC collections and fill the  photon collection
   if ( validPhotonCoreHandle) 
     fillPhotonCollection(theEvent,
              theEventSetup,
              photonCoreHandle,
              topology,
              &barrelRecHits,
              &endcapRecHits,
              hcalTowersHandle,
              vtx,
              outputPhotonCollection,
              iSC);
  
 
   // put the product in the event
   edm::LogInfo("PhotonProducer") << " Put in the event " << iSC << " Photon Candidates \n";
   outputPhotonCollection_p->assign(outputPhotonCollection.begin(),outputPhotonCollection.end());
   theEvent.put( outputPhotonCollection_p, PhotonCollection_);
 
 }