Reference Manual

Go to the documentation of this file.
00001 #include <iostream>
00002 #include <vector>
00003 #include <memory>
00004 
00005 // Framework
00006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00007 #include "FWCore/Utilities/interface/Exception.h"
00008 
00009 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00010 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00011 #include "Geometry/CaloTopology/interface/CaloTopology.h"
00012 #include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
00013 
00014 
00015 #include "DataFormats/VertexReco/interface/Vertex.h"
00016 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00017 #include "DataFormats/EgammaReco/interface/ClusterShape.h"
00018 #include "DataFormats/EgammaCandidates/interface/PhotonCore.h"
00019 #include "DataFormats/EgammaCandidates/interface/Photon.h"
00020 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
00021 #include "DataFormats/EgammaCandidates/interface/Conversion.h"
00022 
00023 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
00024 #include "RecoCaloTools/Selectors/interface/CaloConeSelector.h"
00025 #include "RecoEgamma/EgammaPhotonProducers/interface/PhotonProducer.h"
00026 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaTowerIsolation.h"
00027 
00028 
00029 PhotonProducer::PhotonProducer(const edm::ParameterSet& config) : 
00030   conf_(config)
00031 {
00032 
00033   // use onfiguration file to setup input/output collection names
00034 
00035   photonCoreProducer_   = conf_.getParameter<edm::InputTag>("photonCoreProducer");
00036   barrelEcalHits_   = conf_.getParameter<edm::InputTag>("barrelEcalHits");
00037   endcapEcalHits_   = conf_.getParameter<edm::InputTag>("endcapEcalHits");
00038   vertexProducer_   = conf_.getParameter<std::string>("primaryVertexProducer");
00039   hcalTowers_ = conf_.getParameter<edm::InputTag>("hcalTowers");
00040   hOverEConeSize_   = conf_.getParameter<double>("hOverEConeSize");
00041   highEt_        = conf_.getParameter<double>("highEt");
00042   // R9 value to decide converted/unconverted
00043   minR9Barrel_        = conf_.getParameter<double>("minR9Barrel");
00044   minR9Endcap_        = conf_.getParameter<double>("minR9Endcap");
00045   usePrimaryVertex_ = conf_.getParameter<bool>("usePrimaryVertex");
00046  
00047   edm::ParameterSet posCalcParameters = 
00048     config.getParameter<edm::ParameterSet>("posCalcParameters");
00049   posCalculator_ = PositionCalc(posCalcParameters);
00050 
00051   // Parameters for the position calculation:
00052   //  std::map<std::string,double> providedParameters;
00053   // providedParameters.insert(std::make_pair("LogWeighted",conf_.getParameter<bool>("posCalc_logweight")));
00054   //providedParameters.insert(std::make_pair("T0_barl",conf_.getParameter<double>("posCalc_t0_barl")));
00055   //providedParameters.insert(std::make_pair("T0_endc",conf_.getParameter<double>("posCalc_t0_endc")));
00056   //providedParameters.insert(std::make_pair("T0_endcPresh",conf_.getParameter<double>("posCalc_t0_endcPresh")));
00057   //providedParameters.insert(std::make_pair("W0",conf_.getParameter<double>("posCalc_w0")));
00058   //providedParameters.insert(std::make_pair("X0",conf_.getParameter<double>("posCalc_x0")));
00059   //posCalculator_ = PositionCalc(providedParameters);
00060   // cut values for pre-selection
00061   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("minSCEtBarrel")); 
00062   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("maxHoverEBarrel")); 
00063   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("ecalRecHitSumEtOffsetBarrel")); 
00064   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("ecalRecHitSumEtSlopeBarrel")); 
00065   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("hcalTowerSumEtOffsetBarrel"));
00066   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("hcalTowerSumEtSlopeBarrel"));
00067   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("nTrackSolidConeBarrel"));
00068   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("nTrackHollowConeBarrel"));     
00069   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("trackPtSumSolidConeBarrel"));     
00070   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("trackPtSumHollowConeBarrel"));     
00071   preselCutValuesBarrel_.push_back(conf_.getParameter<double>("sigmaIetaIetaCutBarrel"));     
00072   //  
00073   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("minSCEtEndcap")); 
00074   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("maxHoverEEndcap")); 
00075   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("ecalRecHitSumEtOffsetEndcap")); 
00076   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("ecalRecHitSumEtSlopeEndcap")); 
00077   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("hcalTowerSumEtOffsetEndcap"));
00078   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("hcalTowerSumEtSlopeEndcap"));
00079   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("nTrackSolidConeEndcap"));
00080   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("nTrackHollowConeEndcap"));     
00081   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("trackPtSumSolidConeEndcap"));     
00082   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("trackPtSumHollowConeEndcap"));     
00083   preselCutValuesEndcap_.push_back(conf_.getParameter<double>("sigmaIetaIetaCutEndcap"));     
00084   //
00085   energyCorrectionF = EcalClusterFunctionFactory::get()->create("EcalClusterEnergyCorrection", conf_);
00086 
00087   // Register the product
00088   produces< reco::PhotonCollection >(PhotonCollection_);
00089 
00090 }
00091 
00092 PhotonProducer::~PhotonProducer() {
00093 
00094   delete energyCorrectionF;
00095 
00096 }
00097 
00098 
00099 
00100 void  PhotonProducer::beginRun (edm::Run& r, edm::EventSetup const & theEventSetup) {
00101 
00102     thePhotonIsolationCalculator_ = new PhotonIsolationCalculator();
00103     edm::ParameterSet isolationSumsCalculatorSet = conf_.getParameter<edm::ParameterSet>("isolationSumsCalculatorSet"); 
00104     thePhotonIsolationCalculator_->setup(isolationSumsCalculatorSet);
00105 
00106 }
00107 
00108 void  PhotonProducer::endRun (edm::Run& r, edm::EventSetup const & theEventSetup) {
00109 
00110   delete thePhotonIsolationCalculator_;
00111 
00112 }
00113 
00114 
00115 
00116 
00117 void PhotonProducer::produce(edm::Event& theEvent, const edm::EventSetup& theEventSetup) {
00118 
00119   using namespace edm;
00120   //  nEvt_++;
00121 
00122   reco::PhotonCollection outputPhotonCollection;
00123   std::auto_ptr< reco::PhotonCollection > outputPhotonCollection_p(new reco::PhotonCollection);
00124 
00125 
00126   // Get the PhotonCore collection
00127   bool validPhotonCoreHandle=true;
00128   Handle<reco::PhotonCoreCollection> photonCoreHandle;
00129   theEvent.getByLabel(photonCoreProducer_,photonCoreHandle);
00130   if (!photonCoreHandle.isValid()) {
00131     edm::LogError("PhotonProducer") << "Error! Can't get the product "<<photonCoreProducer_.label();
00132     validPhotonCoreHandle=false;
00133   }
00134 
00135  // Get EcalRecHits
00136   bool validEcalRecHits=true;
00137   Handle<EcalRecHitCollection> barrelHitHandle;
00138   EcalRecHitCollection barrelRecHits;
00139   theEvent.getByLabel(barrelEcalHits_, barrelHitHandle);
00140   if (!barrelHitHandle.isValid()) {
00141     edm::LogError("PhotonProducer") << "Error! Can't get the product "<<barrelEcalHits_.label();
00142     validEcalRecHits=false; 
00143   }
00144   if (  validEcalRecHits)  barrelRecHits = *(barrelHitHandle.product());
00145 
00146   
00147   Handle<EcalRecHitCollection> endcapHitHandle;
00148   theEvent.getByLabel(endcapEcalHits_, endcapHitHandle);
00149   EcalRecHitCollection endcapRecHits;
00150   if (!endcapHitHandle.isValid()) {
00151     edm::LogError("PhotonProducer") << "Error! Can't get the product "<<endcapEcalHits_.label();
00152     validEcalRecHits=false; 
00153   }
00154   if( validEcalRecHits) endcapRecHits = *(endcapHitHandle.product());
00155 
00156 
00157 // get Hcal towers collection 
00158   Handle<CaloTowerCollection> hcalTowersHandle;
00159   theEvent.getByLabel(hcalTowers_, hcalTowersHandle);
00160 
00161 
00162   // get the geometry from the event setup:
00163   theEventSetup.get<CaloGeometryRecord>().get(theCaloGeom_);
00164 
00165   //
00166   // update energy correction function
00167   energyCorrectionF->init(theEventSetup);  
00168 
00169   edm::ESHandle<CaloTopology> pTopology;
00170   theEventSetup.get<CaloTopologyRecord>().get(theCaloTopo_);
00171   const CaloTopology *topology = theCaloTopo_.product();
00172 
00173   // Get the primary event vertex
00174   Handle<reco::VertexCollection> vertexHandle;
00175   reco::VertexCollection vertexCollection;
00176   bool validVertex=true;
00177   if ( usePrimaryVertex_ ) {
00178     theEvent.getByLabel(vertexProducer_, vertexHandle);
00179     if (!vertexHandle.isValid()) {
00180       edm::LogWarning("PhotonProducer") << "Error! Can't get the product primary Vertex Collection "<< "\n";
00181       validVertex=false;
00182     }
00183     if (validVertex) vertexCollection = *(vertexHandle.product());
00184   }
00185   math::XYZPoint vtx(0.,0.,0.);
00186   if (vertexCollection.size()>0) vtx = vertexCollection.begin()->position();
00187 
00188 
00189   int iSC=0; // index in photon collection
00190   // Loop over barrel and endcap SC collections and fill the  photon collection
00191   if ( validPhotonCoreHandle) 
00192     fillPhotonCollection(theEvent,
00193                          theEventSetup,
00194                          photonCoreHandle,
00195                          topology,
00196                          &barrelRecHits,
00197                          &endcapRecHits,
00198                          hcalTowersHandle,
00199                          vtx,
00200                          outputPhotonCollection,
00201                          iSC);
00202  
00203 
00204   // put the product in the event
00205   edm::LogInfo("PhotonProducer") << " Put in the event " << iSC << " Photon Candidates \n";
00206   outputPhotonCollection_p->assign(outputPhotonCollection.begin(),outputPhotonCollection.end());
00207   theEvent.put( outputPhotonCollection_p, PhotonCollection_);
00208 
00209 }
00210 
00211 void PhotonProducer::fillPhotonCollection(edm::Event& evt,
00212                                           edm::EventSetup const & es,
00213                                           const edm::Handle<reco::PhotonCoreCollection> & photonCoreHandle,
00214                                           const CaloTopology* topology,
00215                                           const EcalRecHitCollection* ecalBarrelHits,
00216                                           const EcalRecHitCollection* ecalEndcapHits,
00217                                           const edm::Handle<CaloTowerCollection> & hcalTowersHandle, 
00218                                           math::XYZPoint & vtx,
00219                                           reco::PhotonCollection & outputPhotonCollection, int& iSC) {
00220   
00221   const CaloGeometry* geometry = theCaloGeom_.product();
00222   const CaloSubdetectorGeometry* subDetGeometry =0 ;
00223   const CaloSubdetectorGeometry* geometryES = theCaloGeom_->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
00224   const EcalRecHitCollection* hits = 0 ;
00225   std::vector<double> preselCutValues;
00226   float minR9=0;
00227 
00228 
00229 
00230   for(unsigned int lSC=0; lSC < photonCoreHandle->size(); lSC++) {
00231 
00232     reco::PhotonCoreRef coreRef(reco::PhotonCoreRef(photonCoreHandle, lSC));
00233     reco::SuperClusterRef scRef=coreRef->superCluster();
00234     const reco::SuperCluster* pClus=&(*scRef);
00235     iSC++;
00236 
00237     int subdet = scRef->seed()->hitsAndFractions()[0].first.subdetId();
00238     subDetGeometry =  theCaloGeom_->getSubdetectorGeometry(DetId::Ecal, subdet);
00239 
00240     if (subdet==EcalBarrel) 
00241       { 
00242         preselCutValues = preselCutValuesBarrel_;
00243         minR9=minR9Barrel_;
00244         hits=  ecalBarrelHits;
00245       }
00246     else if  (subdet==EcalEndcap) 
00247       { 
00248         preselCutValues = preselCutValuesEndcap_;
00249         minR9=minR9Endcap_;
00250         hits=  ecalEndcapHits;
00251       }
00252     else
00253       { edm::LogWarning("")<<"PhotonProducer: do not know if it is a barrel or endcap SuperCluster" ; }
00254     
00255         
00256     // SC energy preselection
00257     if (scRef->energy()/cosh(scRef->eta()) <= preselCutValues[0] ) continue;
00258     // calculate HoE
00259     const CaloTowerCollection* hcalTowersColl = hcalTowersHandle.product();
00260     EgammaTowerIsolation towerIso1(hOverEConeSize_,0.,0.,1,hcalTowersColl) ;  
00261     EgammaTowerIsolation towerIso2(hOverEConeSize_,0.,0.,2,hcalTowersColl) ;  
00262     double HoE1=towerIso1.getTowerESum(&(*scRef))/scRef->energy();
00263     double HoE2=towerIso2.getTowerESum(&(*scRef))/scRef->energy(); 
00264     //  std::cout << " PhotonProducer " << HoE1  << "  HoE2 " << HoE2 << std::endl;
00265     // std::cout << " PhotonProducer calcualtion of HoE1 " << HoE1  << "  HoE2 " << HoE2 << std::endl;
00266 
00267     
00268     // recalculate position of seed BasicCluster taking shower depth for unconverted photon
00269     math::XYZPoint unconvPos = posCalculator_.Calculate_Location(scRef->seed()->hitsAndFractions(),hits,subDetGeometry,geometryES);
00270 
00271 
00272     float maxXtal =   EcalClusterTools::eMax( *(scRef->seed()), &(*hits) );
00273     float e1x5    =   EcalClusterTools::e1x5(  *(scRef->seed()), &(*hits), &(*topology)); 
00274     float e2x5    =   EcalClusterTools::e2x5Max(  *(scRef->seed()), &(*hits), &(*topology)); 
00275     float e3x3    =   EcalClusterTools::e3x3(  *(scRef->seed()), &(*hits), &(*topology)); 
00276     float e5x5    =   EcalClusterTools::e5x5( *(scRef->seed()), &(*hits), &(*topology)); 
00277     std::vector<float> cov =  EcalClusterTools::covariances( *(scRef->seed()), &(*hits), &(*topology), geometry); 
00278     float sigmaEtaEta = sqrt(cov[0]);
00279     std::vector<float> locCov =  EcalClusterTools::localCovariances( *(scRef->seed()), &(*hits), &(*topology)); 
00280     float sigmaIetaIeta = sqrt(locCov[0]);
00281 
00282     float r9 =e3x3/(scRef->rawEnergy());
00283     // compute position of ECAL shower
00284     math::XYZPoint caloPosition;
00285     double photonEnergy=0;
00286     if (r9>minR9) {
00287       caloPosition = unconvPos;
00288       // f(eta) correction to e5x5
00289       double deltaE = energyCorrectionF->getValue(*pClus, 1);
00290       if (subdet==EcalBarrel) e5x5 = e5x5 * (1.0 +  deltaE/scRef->rawEnergy() );
00291       photonEnergy=  e5x5    + scRef->preshowerEnergy() ;
00292     } else {
00293       caloPosition = scRef->position();
00294       photonEnergy=scRef->energy();
00295     }
00296     
00297 
00298     //std::cout << " dete " << subdet << " r9 " << r9 <<  " uncorrected e5x5 " <<  e5x5uncor << " corrected e5x5 " << e5x5 << " photon energy " << photonEnergy << std::endl;
00299     
00300     // compute momentum vector of photon from primary vertex and cluster position
00301     math::XYZVector direction = caloPosition - vtx;
00302     math::XYZVector momentum = direction.unit() * photonEnergy ;
00303 
00304     // Create candidate
00305     const reco::Particle::LorentzVector  p4(momentum.x(), momentum.y(), momentum.z(), photonEnergy );
00306     reco::Photon newCandidate(p4, caloPosition, coreRef, vtx);
00307 
00308     // Calculate fiducial flags and isolation variable. Blocked are filled from the isolationCalculator
00309     reco::Photon::FiducialFlags fiducialFlags;
00310     reco::Photon::IsolationVariables isolVarR03, isolVarR04;
00311     thePhotonIsolationCalculator_-> calculate ( &newCandidate,evt,es,fiducialFlags,isolVarR04, isolVarR03);
00312     newCandidate.setFiducialVolumeFlags( fiducialFlags );
00313     newCandidate.setIsolationVariables(isolVarR04, isolVarR03 );
00314   
00316     reco::Photon::ShowerShape  showerShape;
00317     showerShape.e1x5= e1x5;
00318     showerShape.e2x5= e2x5;
00319     showerShape.e3x3= e3x3;
00320     showerShape.e5x5= e5x5;
00321     showerShape.maxEnergyXtal =  maxXtal;
00322     showerShape.sigmaEtaEta =    sigmaEtaEta;
00323     showerShape.sigmaIetaIeta =  sigmaIetaIeta;
00324     showerShape.hcalDepth1OverEcal = HoE1;
00325     showerShape.hcalDepth2OverEcal = HoE2;
00326     newCandidate.setShowerShapeVariables ( showerShape ); 
00327 
00329     bool isLooseEM=true;
00330     if ( newCandidate.pt() < highEt_) { 
00331       if ( newCandidate.hadronicOverEm()                   >= preselCutValues[1] )                                            isLooseEM=false;
00332       if ( newCandidate.ecalRecHitSumEtConeDR04()          > preselCutValues[2]+ preselCutValues[3]*newCandidate.pt() )       isLooseEM=false;
00333       if ( newCandidate.hcalTowerSumEtConeDR04()           > preselCutValues[4]+ preselCutValues[5]*newCandidate.pt() )       isLooseEM=false;
00334       if ( newCandidate.nTrkSolidConeDR04()                > int(preselCutValues[6]) )                                        isLooseEM=false;
00335       if ( newCandidate.nTrkHollowConeDR04()               > int(preselCutValues[7]) )                                        isLooseEM=false;
00336       if ( newCandidate.trkSumPtSolidConeDR04()            > preselCutValues[8] )                                             isLooseEM=false;
00337       if ( newCandidate.trkSumPtHollowConeDR04()           > preselCutValues[9] )                                             isLooseEM=false;
00338       if ( newCandidate.sigmaIetaIeta()                    > preselCutValues[10] )                                            isLooseEM=false;
00339     } 
00340     
00341     
00342     if ( isLooseEM)  
00343       outputPhotonCollection.push_back(newCandidate);
00344       
00345         
00346   }
00347 }
00348
/data/refman/pasoursint/CMSSW_4_2_9_HLT1_bphpatch4/src/RecoEgamma/EgammaPhotonProducers/src/PhotonProducer.cc
