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