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