CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/RecoEgamma/EgammaPhotonProducers/src/ConversionTrackCandidateProducer.cc

Go to the documentation of this file.
00001 #include <iostream>
00002 #include <vector>
00003 #include <memory>
00004 
00005 // Framework
00006 #include "FWCore/Framework/interface/Event.h"
00007 #include "FWCore/Framework/interface/EventSetup.h"
00008 #include "DataFormats/Common/interface/Handle.h"
00009 #include "FWCore/Framework/interface/ESHandle.h"
00010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00011 #include "FWCore/Utilities/interface/Exception.h"
00012 //
00013 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
00014 #include "DataFormats/EgammaTrackReco/interface/TrackCandidateSuperClusterAssociation.h"
00015 #include "DataFormats/EgammaTrackReco/interface/TrackCandidateCaloClusterAssociation.h"
00016 //
00017 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00018 //
00019 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00020 //  Abstract classes for the conversion tracking components
00021 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionSeedFinder.h"
00022 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionTrackFinder.h"
00023 // Class header file
00024 #include "RecoEgamma/EgammaPhotonProducers/interface/ConversionTrackCandidateProducer.h"
00025 //
00026 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
00027 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
00028 #include "RecoTracker/Record/interface/NavigationSchoolRecord.h"
00029 #include "RecoEgamma/EgammaPhotonAlgos/interface/OutInConversionSeedFinder.h"
00030 #include "RecoEgamma/EgammaPhotonAlgos/interface/InOutConversionSeedFinder.h"
00031 #include "RecoEgamma/EgammaPhotonAlgos/interface/OutInConversionTrackFinder.h"
00032 #include "RecoEgamma/EgammaPhotonAlgos/interface/InOutConversionTrackFinder.h"
00033 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaTowerIsolation.h"
00034 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaRecHitIsolation.h"
00035 
00036 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
00037 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00038 #include "CommonTools/Utils/interface/StringToEnumValue.h"
00039 #include "DataFormats/EgammaCandidates/interface/Photon.h"
00040 #include "DataFormats/DetId/interface/DetId.h"
00041 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00042 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00043 
00044 
00045 ConversionTrackCandidateProducer::ConversionTrackCandidateProducer(const edm::ParameterSet& config) : 
00046   conf_(config), 
00047   theNavigationSchool_(0), 
00048   theOutInSeedFinder_(0), 
00049   theOutInTrackFinder_(0), 
00050   theInOutSeedFinder_(0),
00051   theInOutTrackFinder_(0)
00052 {  
00053   //std::cout << "ConversionTrackCandidateProducer CTOR " << "\n";
00054   nEvt_=0;  
00055    
00056   // use onfiguration file to setup input/output collection names
00057  
00058 
00059   bcBarrelCollection_     = conf_.getParameter<edm::InputTag>("bcBarrelCollection");
00060   bcEndcapCollection_     = conf_.getParameter<edm::InputTag>("bcEndcapCollection");
00061   
00062   scHybridBarrelProducer_       = conf_.getParameter<edm::InputTag>("scHybridBarrelProducer");
00063   scIslandEndcapProducer_       = conf_.getParameter<edm::InputTag>("scIslandEndcapProducer");
00064   
00065   OutInTrackCandidateCollection_ = conf_.getParameter<std::string>("outInTrackCandidateCollection");
00066   InOutTrackCandidateCollection_ = conf_.getParameter<std::string>("inOutTrackCandidateCollection");
00067 
00068 
00069   OutInTrackSuperClusterAssociationCollection_ = conf_.getParameter<std::string>("outInTrackCandidateSCAssociationCollection");
00070   InOutTrackSuperClusterAssociationCollection_ = conf_.getParameter<std::string>("inOutTrackCandidateSCAssociationCollection");
00071 
00072   barrelecalCollection_ = conf_.getParameter<edm::InputTag>("barrelEcalRecHitCollection");
00073   endcapecalCollection_ = conf_.getParameter<edm::InputTag>("endcapEcalRecHitCollection");
00074   
00075   hcalTowers_        = conf_.getParameter<edm::InputTag>("hcalTowers");
00076   hOverEConeSize_    = conf_.getParameter<double>("hOverEConeSize");
00077   maxHOverE_         = conf_.getParameter<double>("maxHOverE");
00078   minSCEt_           = conf_.getParameter<double>("minSCEt");
00079   isoConeR_          = conf_.getParameter<double>("isoConeR");
00080   isoInnerConeR_     = conf_.getParameter<double>("isoInnerConeR");
00081   isoEtaSlice_       = conf_.getParameter<double>("isoEtaSlice");
00082   isoEtMin_          = conf_.getParameter<double>("isoEtMin");
00083   isoEMin_           = conf_.getParameter<double>("isoEMin");
00084   vetoClusteredHits_ = conf_.getParameter<bool>("vetoClusteredHits");
00085   useNumXtals_       = conf_.getParameter<bool>("useNumXstals");
00086    ecalIsoCut_offset_ = conf_.getParameter<double>("ecalIsoCut_offset");
00087   ecalIsoCut_slope_  = conf_.getParameter<double>("ecalIsoCut_slope");
00088 
00089   //Flags and Severities to be excluded from photon calculations
00090   const std::vector<std::string> flagnamesEB = 
00091     config.getParameter<std::vector<std::string> >("RecHitFlagToBeExcludedEB");
00092 
00093   const std::vector<std::string> flagnamesEE =
00094     config.getParameter<std::vector<std::string> >("RecHitFlagToBeExcludedEE");
00095 
00096   flagsexclEB_= 
00097     StringToEnumValue<EcalRecHit::Flags>(flagnamesEB);
00098 
00099   flagsexclEE_=
00100     StringToEnumValue<EcalRecHit::Flags>(flagnamesEE);
00101 
00102   const std::vector<std::string> severitynamesEB = 
00103     config.getParameter<std::vector<std::string> >("RecHitSeverityToBeExcludedEB");
00104 
00105   severitiesexclEB_= 
00106     StringToEnumValue<EcalSeverityLevel::SeverityLevel>(severitynamesEB);
00107 
00108   const std::vector<std::string> severitynamesEE = 
00109     config.getParameter<std::vector<std::string> >("RecHitSeverityToBeExcludedEE");
00110 
00111   severitiesexclEE_= 
00112     StringToEnumValue<EcalSeverityLevel::SeverityLevel>(severitynamesEE);
00113 
00114 
00115   // Register the product
00116   produces< TrackCandidateCollection > (OutInTrackCandidateCollection_);
00117   produces< TrackCandidateCollection > (InOutTrackCandidateCollection_);
00118 
00119   produces< reco::TrackCandidateCaloClusterPtrAssociation > ( OutInTrackSuperClusterAssociationCollection_);
00120   produces< reco::TrackCandidateCaloClusterPtrAssociation > ( InOutTrackSuperClusterAssociationCollection_);
00121   
00122 
00123 }
00124 
00125 ConversionTrackCandidateProducer::~ConversionTrackCandidateProducer() {}
00126 
00127 void  ConversionTrackCandidateProducer::setEventSetup (const edm::EventSetup & theEventSetup) {
00128 
00129 
00130   theOutInSeedFinder_->setEventSetup(theEventSetup);
00131   theInOutSeedFinder_->setEventSetup(theEventSetup);
00132   theOutInTrackFinder_->setEventSetup(theEventSetup);
00133   theInOutTrackFinder_->setEventSetup(theEventSetup);
00134 
00135 
00136 }
00137 
00138 
00139 void  ConversionTrackCandidateProducer::beginRun (edm::Run& r , edm::EventSetup const & theEventSetup) {
00140 
00141   edm::ESHandle<NavigationSchool> nav;
00142   theEventSetup.get<NavigationSchoolRecord>().get("SimpleNavigationSchool", nav);
00143   theNavigationSchool_ = nav.product();
00144 
00145   // get the Out In Seed Finder  
00146   theOutInSeedFinder_ = new OutInConversionSeedFinder (  conf_ );
00147   
00148   // get the Out In Track Finder
00149   theOutInTrackFinder_ = new OutInConversionTrackFinder ( theEventSetup, conf_  );
00150 
00151   
00152   // get the In Out Seed Finder  
00153   theInOutSeedFinder_ = new InOutConversionSeedFinder ( conf_ );
00154   
00155   
00156   // get the In Out Track Finder
00157   theInOutTrackFinder_ = new InOutConversionTrackFinder ( theEventSetup, conf_  );
00158 }
00159 
00160 
00161 void  ConversionTrackCandidateProducer::endRun (edm::Run& r , edm::EventSetup const & theEventSetup) {
00162   delete theOutInSeedFinder_; 
00163   delete theOutInTrackFinder_;
00164   delete theInOutSeedFinder_;  
00165   delete theInOutTrackFinder_;
00166 }
00167 
00168 
00169 
00170 
00171 void ConversionTrackCandidateProducer::produce(edm::Event& theEvent, const edm::EventSetup& theEventSetup) {
00172   
00173   using namespace edm;
00174   nEvt_++;
00175   //  std::cout << "ConversionTrackCandidateProducer Analyzing event number " <<   theEvent.id() <<  " Global Counter " << nEvt_ << "\n";
00176   
00177 
00178   
00179   setEventSetup( theEventSetup );
00180   theOutInSeedFinder_->setEvent(theEvent);
00181   theInOutSeedFinder_->setEvent(theEvent);
00182   theOutInTrackFinder_->setEvent(theEvent);
00183   theInOutTrackFinder_->setEvent(theEvent);
00184 
00185 // Set the navigation school  
00186   NavigationSetter setter(*theNavigationSchool_);  
00187 
00188   //
00189   // create empty output collections
00190   //
00191   //  Out In Track Candidates
00192   std::auto_ptr<TrackCandidateCollection> outInTrackCandidate_p(new TrackCandidateCollection); 
00193   //  In Out  Track Candidates
00194   std::auto_ptr<TrackCandidateCollection> inOutTrackCandidate_p(new TrackCandidateCollection); 
00195   //   Track Candidate  calo  Cluster Association
00196   std::auto_ptr<reco::TrackCandidateCaloClusterPtrAssociation> outInAssoc_p(new reco::TrackCandidateCaloClusterPtrAssociation);
00197   std::auto_ptr<reco::TrackCandidateCaloClusterPtrAssociation> inOutAssoc_p(new reco::TrackCandidateCaloClusterPtrAssociation);
00198     
00199   // Get the basic cluster collection in the Barrel 
00200   bool validBarrelBCHandle=true;
00201   edm::Handle<edm::View<reco::CaloCluster> > bcBarrelHandle;
00202   theEvent.getByLabel(bcBarrelCollection_, bcBarrelHandle);
00203   if (!bcBarrelHandle.isValid()) {
00204     edm::LogError("ConversionTrackCandidateProducer") << "Error! Can't get the product "<<bcBarrelCollection_.label();
00205     validBarrelBCHandle=false;
00206   }
00207   
00208   
00209   // Get the basic cluster collection in the Endcap 
00210   bool validEndcapBCHandle=true;
00211   edm::Handle<edm::View<reco::CaloCluster> > bcEndcapHandle;
00212   theEvent.getByLabel(bcEndcapCollection_, bcEndcapHandle);
00213   if (!bcEndcapHandle.isValid()) {
00214     edm::LogError("CoonversionTrackCandidateProducer") << "Error! Can't get the product "<<bcEndcapCollection_.label();
00215     validEndcapBCHandle=false; 
00216   }
00217   
00218   
00219 
00220   // Get the Super Cluster collection in the Barrel
00221   bool validBarrelSCHandle=true;
00222   edm::Handle<edm::View<reco::CaloCluster> > scBarrelHandle;
00223   theEvent.getByLabel(scHybridBarrelProducer_,scBarrelHandle);
00224   if (!scBarrelHandle.isValid()) {
00225     edm::LogError("CoonversionTrackCandidateProducer") << "Error! Can't get the product "<<scHybridBarrelProducer_.label();
00226     validBarrelSCHandle=false;
00227   }
00228 
00229 
00230   // Get the Super Cluster collection in the Endcap
00231   bool validEndcapSCHandle=true;
00232   edm::Handle<edm::View<reco::CaloCluster> > scEndcapHandle;
00233   theEvent.getByLabel(scIslandEndcapProducer_,scEndcapHandle);
00234   if (!scEndcapHandle.isValid()) {
00235     edm::LogError("CoonversionTrackCandidateProducer") << "Error! Can't get the product "<<scIslandEndcapProducer_.label();
00236     validEndcapSCHandle=false;
00237   }
00238 
00239 
00240   // get the geometry from the event setup:
00241   theEventSetup.get<CaloGeometryRecord>().get(theCaloGeom_);
00242 
00243   // get Hcal towers collection 
00244   Handle<CaloTowerCollection> hcalTowersHandle;
00245   theEvent.getByLabel(hcalTowers_, hcalTowersHandle);
00246 
00247   edm::Handle<EcalRecHitCollection> ecalhitsCollEB;
00248   edm::Handle<EcalRecHitCollection> ecalhitsCollEE;
00249 
00250   theEvent.getByLabel(endcapecalCollection_, ecalhitsCollEE);
00251   theEvent.getByLabel(barrelecalCollection_, ecalhitsCollEB);
00252 
00253   edm::ESHandle<EcalSeverityLevelAlgo> sevlv;
00254   theEventSetup.get<EcalSeverityLevelAlgoRcd>().get(sevlv);
00255   const EcalSeverityLevelAlgo* sevLevel = sevlv.product();
00256 
00257   std::auto_ptr<CaloRecHitMetaCollectionV> RecHitsEE(0); 
00258   RecHitsEE = std::auto_ptr<CaloRecHitMetaCollectionV>(new EcalRecHitMetaCollection(ecalhitsCollEE.product()));
00259  
00260   std::auto_ptr<CaloRecHitMetaCollectionV> RecHitsEB(0); 
00261   RecHitsEB = std::auto_ptr<CaloRecHitMetaCollectionV>(new EcalRecHitMetaCollection(ecalhitsCollEB.product()));
00262 
00263 
00264   caloPtrVecOutIn_.clear();
00265   caloPtrVecInOut_.clear();
00266 
00267   bool isBarrel=true;
00268   if ( validBarrelBCHandle && validBarrelSCHandle ) 
00269     buildCollections(isBarrel, scBarrelHandle, bcBarrelHandle, ecalhitsCollEB, &(*RecHitsEB), sevLevel, hcalTowersHandle, *outInTrackCandidate_p, *inOutTrackCandidate_p, caloPtrVecOutIn_, caloPtrVecInOut_);
00270 
00271   if ( validEndcapBCHandle && validEndcapSCHandle ) {
00272     isBarrel=false; 
00273     buildCollections(isBarrel, scEndcapHandle, bcEndcapHandle, ecalhitsCollEE, &(*RecHitsEE), sevLevel, hcalTowersHandle, *outInTrackCandidate_p, *inOutTrackCandidate_p, caloPtrVecOutIn_, caloPtrVecInOut_);
00274   }
00275 
00276 
00277   //  std::cout  << "  ConversionTrackCandidateProducer  caloPtrVecOutIn_ size " <<  caloPtrVecOutIn_.size() << " caloPtrVecInOut_ size " << caloPtrVecInOut_.size()  << "\n"; 
00278   
00279 
00280 
00281   // put all products in the event
00282  // Barrel
00283  //std::cout  << "ConversionTrackCandidateProducer Putting in the event " << (*outInTrackCandidate_p).size() << " Out In track Candidates " << "\n";
00284  const edm::OrphanHandle<TrackCandidateCollection> refprodOutInTrackC = theEvent.put( outInTrackCandidate_p, OutInTrackCandidateCollection_ );
00285  //std::cout  << "ConversionTrackCandidateProducer  refprodOutInTrackC size  " <<  (*(refprodOutInTrackC.product())).size()  <<  "\n";
00286  //
00287  //std::cout  << "ConversionTrackCandidateProducer Putting in the event  " << (*inOutTrackCandidate_p).size() << " In Out track Candidates " <<  "\n";
00288  const edm::OrphanHandle<TrackCandidateCollection> refprodInOutTrackC = theEvent.put( inOutTrackCandidate_p, InOutTrackCandidateCollection_ );
00289  //std::cout  << "ConversionTrackCandidateProducer  refprodInOutTrackC size  " <<  (*(refprodInOutTrackC.product())).size()  <<  "\n";
00290 
00291 
00292  edm::ValueMap<reco::CaloClusterPtr>::Filler fillerOI(*outInAssoc_p);
00293  fillerOI.insert(refprodOutInTrackC, caloPtrVecOutIn_.begin(), caloPtrVecOutIn_.end());
00294  fillerOI.fill();
00295  edm::ValueMap<reco::CaloClusterPtr>::Filler fillerIO(*inOutAssoc_p);
00296  fillerIO.insert(refprodInOutTrackC, caloPtrVecInOut_.begin(), caloPtrVecInOut_.end());
00297  fillerIO.fill();
00298 
00299 
00300   
00301  // std::cout  << "ConversionTrackCandidateProducer Putting in the event   OutIn track - SC association: size  " <<  (*outInAssoc_p).size() << "\n";  
00302  theEvent.put( outInAssoc_p, OutInTrackSuperClusterAssociationCollection_);
00303  
00304  // std::cout << "ConversionTrackCandidateProducer Putting in the event   InOut track - SC association: size  " <<  (*inOutAssoc_p).size() << "\n";  
00305  theEvent.put( inOutAssoc_p, InOutTrackSuperClusterAssociationCollection_);
00306 
00307  theOutInSeedFinder_->clear();
00308  theInOutSeedFinder_->clear();
00309  
00310 
00311   
00312 }
00313 
00314 
00315 void ConversionTrackCandidateProducer::buildCollections(bool isBarrel, 
00316                                                         const edm::Handle<edm::View<reco::CaloCluster> > & scHandle,
00317                                                         const edm::Handle<edm::View<reco::CaloCluster> > & bcHandle,
00318                                                         edm::Handle<EcalRecHitCollection> ecalRecHitHandle, 
00319                                                         CaloRecHitMetaCollectionV* ecalRecHits,
00320                                                         const EcalSeverityLevelAlgo* sevLevel,
00321                                                         //edm::ESHandle<EcalChannelStatus>  chStatus,
00322                                                         //const EcalChannelStatus* chStatus,
00323                                                         const edm::Handle<CaloTowerCollection> & hcalTowersHandle,
00324                                                         TrackCandidateCollection& outInTrackCandidates,
00325                                                         TrackCandidateCollection& inOutTrackCandidates,
00326                                                         std::vector<edm::Ptr<reco::CaloCluster> >& vecRecOI,
00327                                                         std::vector<edm::Ptr<reco::CaloCluster> >& vecRecIO )
00328 
00329 {
00330 
00331   //  std::cout << "ConversionTrackCandidateProducer builcollections bc size " << bcHandle->size() <<  "\n";
00332   //const CaloGeometry* geometry = theCaloGeom_.product();
00333 
00334   //  Loop over SC in the barrel and reconstruct converted photons
00335   for (unsigned i = 0; i < scHandle->size(); ++i ) {
00336 
00337     reco::CaloClusterPtr aClus= scHandle->ptrAt(i);
00338   
00339     // preselection based in Et and H/E cut. 
00340     if (aClus->energy()/cosh(aClus->eta()) <= minSCEt_) continue;
00341     const reco::CaloCluster*   pClus=&(*aClus);
00342     const reco::SuperCluster*  sc=dynamic_cast<const reco::SuperCluster*>(pClus);
00343     double scEt = sc->energy()/cosh(sc->eta());  
00344     const CaloTowerCollection* hcalTowersColl = hcalTowersHandle.product();
00345     EgammaTowerIsolation towerIso(hOverEConeSize_,0.,0.,-1,hcalTowersColl) ;
00346     double HoE = towerIso.getTowerESum(sc)/sc->energy();
00347     if (HoE >= maxHOverE_)  continue;
00348 
00350     EgammaRecHitIsolation ecalIso(isoConeR_,     
00351                                   isoInnerConeR_, 
00352                                   isoEtaSlice_,  
00353                                   isoEtMin_,    
00354                                   isoEMin_,    
00355                                   theCaloGeom_,
00356                                   &(*ecalRecHits),
00357                                   sevLevel,
00358                                   DetId::Ecal);
00359 
00360     ecalIso.setVetoClustered(vetoClusteredHits_);
00361     ecalIso.setUseNumCrystals(useNumXtals_);
00362     if (isBarrel) {
00363       ecalIso.doFlagChecks(flagsexclEB_);
00364       ecalIso.doSeverityChecks(ecalRecHitHandle.product(), severitiesexclEB_);
00365     } else {
00366       ecalIso.doFlagChecks(flagsexclEE_);
00367       ecalIso.doSeverityChecks(ecalRecHitHandle.product(), severitiesexclEE_);
00368     }
00369 
00370     double ecalIsolation = ecalIso.getEtSum(sc);
00371     if ( ecalIsolation >   ecalIsoCut_offset_ + ecalIsoCut_slope_*scEt ) continue;
00372 
00373     // Now launch the seed finding
00374     theOutInSeedFinder_->setCandidate(pClus->energy(), GlobalPoint(pClus->position().x(),pClus->position().y(),pClus->position().z() ) );
00375     theOutInSeedFinder_->makeSeeds( bcHandle );
00376 
00377     std::vector<Trajectory> theOutInTracks= theOutInTrackFinder_->tracks(theOutInSeedFinder_->seeds(),  outInTrackCandidates);    
00378 
00379     theInOutSeedFinder_->setCandidate(pClus->energy(), GlobalPoint(pClus->position().x(),pClus->position().y(),pClus->position().z() ) );  
00380     theInOutSeedFinder_->setTracks(  theOutInTracks );   
00381     theInOutSeedFinder_->makeSeeds(  bcHandle);
00382     
00383     std::vector<Trajectory> theInOutTracks= theInOutTrackFinder_->tracks(theInOutSeedFinder_->seeds(),  inOutTrackCandidates); 
00384 
00385     // Debug
00386     //   std::cout  << "ConversionTrackCandidateProducer  theOutInTracks.size() " << theOutInTracks.size() << " theInOutTracks.size() " << theInOutTracks.size() <<  " Event pointer to out in track size barrel " << outInTrackCandidates.size() << " in out track size " << inOutTrackCandidates.size() <<   "\n";
00387 
00388 
00390     for (std::vector<Trajectory>::const_iterator it = theOutInTracks.begin(); it !=  theOutInTracks.end(); ++it) {
00391       caloPtrVecOutIn_.push_back(aClus);
00392       //     std::cout  << "ConversionTrackCandidateProducer Barrel OutIn Tracks Number of hits " << (*it).foundHits() << "\n"; 
00393     }
00394 
00395     for (std::vector<Trajectory>::const_iterator it = theInOutTracks.begin(); it !=  theInOutTracks.end(); ++it) {
00396       caloPtrVecInOut_.push_back(aClus);
00397       //     std::cout  << "ConversionTrackCandidateProducer Barrel InOut Tracks Number of hits " << (*it).foundHits() << "\n"; 
00398     }
00399   }
00400 }
00401