Reference Manual

Go to the documentation of this file.
00001 #include "RecoTauTag/HLTProducers/interface/L2TauModularIsolationProducer.h"
00002 #include "RecoTauTag/HLTProducers/interface/L2TauIsolationAlgs.h"
00003 #include "RecoTauTag/HLTProducers/interface/L2TauSimpleClustering.h"
00004 
00005 
00006 using namespace reco;
00007 using namespace edm;
00008 
00009 
00010 L2TauModularIsolationProducer::L2TauModularIsolationProducer(const edm::ParameterSet& iConfig):
00011   l2CaloJets_(iConfig.getParameter<edm::InputTag>("L2TauJetCollection")),
00012   EBRecHits_(iConfig.getParameter<edm::InputTag>("EBRecHits")),
00013   EERecHits_(iConfig.getParameter<edm::InputTag>("EERecHits")),
00014   caloTowers_(iConfig.getParameter<edm::InputTag>("CaloTowers")),
00015   pfClustersECAL_(iConfig.getParameter<edm::InputTag>("pfClustersECAL")),
00016   pfClustersHCAL_(iConfig.getParameter<edm::InputTag>("pfClustersHCAL")),
00017   ecalIsolationAlg_(iConfig.getParameter<std::string>("ecalIsolationAlgorithm")),
00018   hcalIsolationAlg_(iConfig.getParameter<std::string>("hcalIsolationAlgorithm")),
00019   ecalClusteringAlg_(iConfig.getParameter<std::string>("ecalClusteringAlgorithm")),
00020   hcalClusteringAlg_(iConfig.getParameter<std::string>("hcalClusteringAlgorithm")),
00021   associationRadius_(iConfig.getParameter<double>("associationRadius")),
00022   simpleClusterRadiusECAL_(iConfig.getParameter<double>("simpleClusterRadiusEcal")),
00023   simpleClusterRadiusHCAL_(iConfig.getParameter<double>("simpleClusterRadiusHcal")),
00024   innerConeECAL_(iConfig.getParameter<double>("innerConeECAL")),
00025   outerConeECAL_(iConfig.getParameter<double>("outerConeECAL")),
00026   innerConeHCAL_(iConfig.getParameter<double>("innerConeHCAL")),
00027   outerConeHCAL_(iConfig.getParameter<double>("outerConeHCAL")),
00028   crystalThresholdE_(iConfig.getParameter<double>("crystalThresholdEE")),
00029   crystalThresholdB_(iConfig.getParameter<double>("crystalThresholdEB")),
00030   towerThreshold_(iConfig.getParameter<double>("towerThreshold"))
00031 
00032  {
00033        
00034   //Add the products
00035   produces<L2TauInfoAssociation>();
00036 
00037 }
00038 
00039 
00040 L2TauModularIsolationProducer::~L2TauModularIsolationProducer()
00041 {
00042   //Destruction
00043 
00044 }
00045 
00046 
00047 
00048 void
00049 L2TauModularIsolationProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00050 {
00051    edm::Handle<CaloJetCollection> l2CaloJets; //Handle to the input (L2TauCaloJets);
00052    iEvent.getByLabel(l2CaloJets_ ,l2CaloJets);//get the handle
00053 
00054    //Create the Association
00055    std::auto_ptr<L2TauInfoAssociation> l2InfoAssoc( new L2TauInfoAssociation);
00056 
00057    //If the JetCrystalsAssociation exists -> RUN The Producer
00058    if(l2CaloJets->size()>0)
00059      {
00060       CaloJetCollection::const_iterator jcStart = l2CaloJets->begin();
00061        //Loop on Jets
00062        for(CaloJetCollection::const_iterator jc = jcStart ;jc!=l2CaloJets->end();++jc)
00063          {
00064            //Create Algorithm Object
00065            L2TauIsolationAlgs alg;
00066            //Create Info Object
00067            L2TauIsolationInfo info; 
00068 
00069            //Define objects to be loaded from file:
00070            math::PtEtaPhiELorentzVectorCollection hitsECAL;
00071            math::PtEtaPhiELorentzVectorCollection hitsHCAL;
00072            math::PtEtaPhiELorentzVectorCollection pfClustersECAL;
00073            math::PtEtaPhiELorentzVectorCollection pfClustersHCAL;
00074 
00075            if(ecalIsolationAlg_=="recHits" || ecalClusteringAlg_=="recHits" || ecalIsolationAlg_ =="simpleClusters" ||ecalClusteringAlg_ =="simpleClusters")
00076              hitsECAL =getECALHits(*jc,iEvent,iSetup); 
00077 
00078            if(hcalIsolationAlg_=="recHits" || hcalClusteringAlg_=="recHits" || hcalIsolationAlg_ =="simpleClusters" ||hcalClusteringAlg_ =="simpleClusters")
00079              hitsHCAL =getHCALHits(*jc,iEvent); 
00080 
00081            if(ecalIsolationAlg_=="particleFlow" || ecalClusteringAlg_=="particleFlow")
00082              pfClustersECAL =getPFClusters(*jc,iEvent,pfClustersECAL_);
00083 
00084            if(hcalIsolationAlg_=="particleFlow" || hcalClusteringAlg_=="particleFlow")
00085              pfClustersHCAL =getPFClusters(*jc,iEvent,pfClustersHCAL_);
00086 
00087 
00088 
00089            //Do ECAL Isolation
00090 
00091            if(ecalIsolationAlg_ == "recHits")
00092              {
00093                //Use Rechits
00094                info.setEcalIsolEt( alg.isolatedEt(hitsECAL , jc->p4().Vect(), innerConeECAL_,outerConeECAL_) );
00095              }
00096            else if(ecalIsolationAlg_ == "simpleClusters")
00097              {
00098                //create the simple clusters
00099                L2TauSimpleClustering clustering(simpleClusterRadiusECAL_);
00100                math::PtEtaPhiELorentzVectorCollection clusters = clustering.clusterize(hitsECAL);
00101                info.setEcalIsolEt( alg.isolatedEt(clusters , jc->p4().Vect(), innerConeECAL_,outerConeECAL_) );
00102 
00103              }
00104            else if(ecalIsolationAlg_ == "particleFlow")
00105              {
00106                //Use ParticleFlow
00107                info.setEcalIsolEt( alg.isolatedEt(pfClustersECAL , jc->p4().Vect(), innerConeECAL_,outerConeECAL_) );
00108              }
00109 
00110            //Do ECAL Clustering
00111 
00112            if(ecalClusteringAlg_ == "recHits")
00113              {
00114                //Use Rechits
00115                info.setEcalClusterShape(alg.clusterShape(hitsECAL,jc->p4().Vect(),0.,outerConeECAL_) );
00116                info.setNEcalHits( alg.nClustersAnnulus(hitsECAL , jc->p4().Vect(), innerConeECAL_,outerConeECAL_));
00117                if(hitsECAL.size()>0)
00118                  info.setSeedEcalHitEt(hitsECAL[0].pt());
00119 
00120              }
00121            else if(ecalClusteringAlg_ == "simpleClusters")
00122              {
00123                //create the simple clusters
00124                L2TauSimpleClustering clustering(simpleClusterRadiusECAL_);
00125                math::PtEtaPhiELorentzVectorCollection clusters = clustering.clusterize(hitsECAL);
00126                info.setEcalClusterShape(alg.clusterShape(clusters,jc->p4().Vect(),0.,outerConeECAL_) );
00127                info.setNEcalHits( alg.nClustersAnnulus(clusters, jc->p4().Vect(), innerConeECAL_,outerConeECAL_));
00128 
00129                if(clusters.size()>0)
00130                  info.setSeedEcalHitEt(clusters[0].pt());
00131              }
00132            else if(ecalClusteringAlg_ == "particleFlow")
00133              {
00134                //Use ParticleFlow
00135                info.setEcalClusterShape(alg.clusterShape(pfClustersECAL,jc->p4().Vect(),0.,outerConeECAL_) );
00136                info.setNEcalHits( alg.nClustersAnnulus(pfClustersECAL, jc->p4().Vect(), innerConeECAL_,outerConeECAL_));
00137                if(pfClustersECAL.size()>0)
00138                  info.setSeedEcalHitEt(pfClustersECAL[0].pt());
00139              }
00140 
00141 
00142            //Do HCAL Isolation
00143 
00144            if(hcalIsolationAlg_ == "recHits")
00145              {
00146                //Use Rechits
00147                info.setHcalIsolEt( alg.isolatedEt(hitsHCAL , jc->p4().Vect(), innerConeHCAL_,outerConeHCAL_) );
00148              }
00149            else if(hcalIsolationAlg_ == "simpleClusters")
00150              {
00151                //create the simple clusters
00152                L2TauSimpleClustering clustering(simpleClusterRadiusHCAL_);
00153                math::PtEtaPhiELorentzVectorCollection clusters = clustering.clusterize(hitsHCAL);
00154                info.setHcalIsolEt( alg.isolatedEt(clusters , jc->p4().Vect(), innerConeHCAL_,outerConeHCAL_) );
00155 
00156              }
00157            else if(hcalIsolationAlg_ == "particleFlow")
00158              {
00159                //Use ParticleFlow
00160                info.setHcalIsolEt( alg.isolatedEt(pfClustersHCAL , jc->p4().Vect(), innerConeHCAL_,outerConeHCAL_) );
00161              }
00162 
00163 
00164 
00165            //Do HCAL Clustering
00166 
00167            if(hcalClusteringAlg_ == "recHits")
00168              {
00169                //Use Rechits
00170                info.setHcalClusterShape(alg.clusterShape(hitsHCAL,jc->p4().Vect(),0.,outerConeHCAL_) );
00171                info.setNHcalHits( alg.nClustersAnnulus(hitsHCAL, jc->p4().Vect(), innerConeHCAL_,outerConeHCAL_));
00172                if(hitsHCAL.size()>0)
00173                  info.setSeedHcalHitEt(hitsHCAL[0].pt());
00174              }
00175            else if(hcalClusteringAlg_ == "simpleClusters")
00176              {
00177                //create the simple clusters
00178                L2TauSimpleClustering clustering(simpleClusterRadiusHCAL_);
00179                math::PtEtaPhiELorentzVectorCollection clusters = clustering.clusterize(hitsHCAL);
00180                info.setHcalClusterShape(alg.clusterShape(clusters,jc->p4().Vect(),0.,outerConeHCAL_) );
00181                info.setNHcalHits( alg.nClustersAnnulus(clusters, jc->p4().Vect(), innerConeHCAL_,outerConeHCAL_));
00182                if(clusters.size()>0)
00183                  info.setSeedHcalHitEt(clusters[0].pt());
00184              }
00185            else if(hcalClusteringAlg_ == "particleFlow")
00186              {
00187                //Use ParticleFlow
00188                info.setHcalClusterShape(alg.clusterShape(pfClustersHCAL,jc->p4().Vect(),0.,outerConeHCAL_) );
00189                info.setNHcalHits( alg.nClustersAnnulus(pfClustersHCAL, jc->p4().Vect(), innerConeHCAL_,outerConeHCAL_));
00190                if(pfClustersHCAL.size()>0)
00191                  info.setSeedHcalHitEt(pfClustersHCAL[0].pt());
00192              }
00193 
00194            //Store the info Class
00195            edm::Ref<CaloJetCollection> jcRef(l2CaloJets,jc-jcStart);
00196            l2InfoAssoc->insert(jcRef, info);
00197          }
00198 
00199      } //end of if(*jetCrystalsObj)
00200 
00201     iEvent.put(l2InfoAssoc);
00202 }
00203 
00204 
00205 void 
00206 L2TauModularIsolationProducer::beginJob()
00207 {
00208 }
00209 
00210 
00211 void 
00212 L2TauModularIsolationProducer::endJob() {
00213 }
00214 
00215 
00216 
00217 
00218 math::PtEtaPhiELorentzVectorCollection 
00219 L2TauModularIsolationProducer::getHCALHits(const CaloJet& jet,const edm::Event& iEvent)
00220 {
00221   edm::Handle<CaloTowerCollection> towers;
00222 
00223   math::PtEtaPhiELorentzVectorCollection towers2;
00224 
00225   if(iEvent.getByLabel(caloTowers_,towers))
00226     if(towers->size()>0)
00227     for(size_t i=0;i<towers->size();++i)
00228       {
00229         math::PtEtaPhiELorentzVector tower((*towers)[i].et(),(*towers)[i].eta(),(*towers)[i].phi(),(*towers)[i].energy());
00230         if(ROOT::Math::VectorUtil::DeltaR(tower,jet.p4()) <associationRadius_)
00231           {
00232             if(tower.pt()>towerThreshold_)
00233               towers2.push_back(tower);
00234           }       
00235       }
00236 
00237   if(towers2.size()>0)
00238     std::sort(towers2.begin(),towers2.end(),comparePt);
00239   return towers2;
00240 }
00241 
00242 
00243 math::PtEtaPhiELorentzVectorCollection 
00244 L2TauModularIsolationProducer::getECALHits(const CaloJet& jet,const edm::Event& iEvent,const edm::EventSetup& iSetup)
00245 {
00246   //Init Geometry
00247   ESHandle<CaloGeometry> geometry;
00248   iSetup.get<CaloGeometryRecord>().get(geometry);
00249 
00250   //Create ECAL Geometry
00251   const CaloSubdetectorGeometry* EB = geometry->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
00252   const CaloSubdetectorGeometry* EE = geometry->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);
00253 
00254   //Handle To the ECAL
00255   edm::Handle<EBRecHitCollection> EBRecHits;
00256   edm::Handle<EERecHitCollection> EERecHits;
00257 
00258   //Create a container for the hits
00259   math::PtEtaPhiELorentzVectorCollection jetRecHits;
00260 
00261   //Loop on the barrel hits
00262   if(iEvent.getByLabel( EBRecHits_, EBRecHits))
00263      for(EBRecHitCollection::const_iterator hit = EBRecHits->begin();hit!=EBRecHits->end();++hit)
00264        {
00265          //get Detector Geometry
00266          const CaloCellGeometry* this_cell = EB->getGeometry(hit->detid());
00267          GlobalPoint posi = this_cell->getPosition();
00268          double energy = hit->energy();
00269          double eta = posi.eta();
00270          double phi = posi.phi();
00271          double theta = posi.theta();
00272          if(theta > M_PI) theta = 2 * M_PI- theta;
00273          double et = energy * sin(theta);
00274          math::PtEtaPhiELorentzVector p(et, eta, phi, energy);
00275          if(ROOT::Math::VectorUtil::DeltaR(p,jet.p4()) <associationRadius_)
00276            if(p.pt()>crystalThresholdB_)
00277              jetRecHits.push_back(p);
00278        }
00279 
00280  if(iEvent.getByLabel( EERecHits_, EERecHits))
00281      for(EERecHitCollection::const_iterator hit = EERecHits->begin();hit!=EERecHits->end();++hit)
00282        {
00283          //get Detector Geometry
00284          const CaloCellGeometry* this_cell = EE->getGeometry(hit->detid());
00285          GlobalPoint posi = this_cell->getPosition();
00286          double energy = hit->energy();
00287          double eta = posi.eta();
00288          double phi = posi.phi();
00289          double theta = posi.theta();
00290          if(theta > M_PI) theta = 2 * M_PI- theta;
00291          double et = energy * sin(theta);
00292          math::PtEtaPhiELorentzVector p(et, eta, phi, energy);
00293          if(ROOT::Math::VectorUtil::DeltaR(p,jet.p4()) < associationRadius_)
00294            if(p.pt()>crystalThresholdE_)
00295              jetRecHits.push_back(p);
00296        }
00297 
00298  if(jetRecHits.size()>0)
00299    std::sort(jetRecHits.begin(),jetRecHits.end(),comparePt);
00300 
00301  return jetRecHits;
00302 }
00303 
00304 
00305 
00306 math::PtEtaPhiELorentzVectorCollection 
00307 L2TauModularIsolationProducer::getPFClusters(const CaloJet& jet,const edm::Event& iEvent,const edm::InputTag& input)
00308 {
00309   edm::Handle<PFClusterCollection> clusters;
00310   math::PtEtaPhiELorentzVectorCollection clusters2;
00311 
00312   //get Clusters near the jet
00313   if(iEvent.getByLabel(input,clusters))
00314     if(clusters->size()>0)
00315     for(PFClusterCollection::const_iterator c = clusters->begin();c!=clusters->end();++c)
00316     {
00317       double energy = c->energy();
00318       double eta = c->eta();
00319       double phi = c->phi();
00320       double theta = c->position().theta();
00321       if(theta > M_PI) theta = 2 * M_PI- theta;
00322       double et = energy * sin(theta);
00323       math::PtEtaPhiELorentzVector p(et, eta, phi, energy);
00324 
00325       if(ROOT::Math::VectorUtil::DeltaR(p,jet.p4()) < associationRadius_)
00326         clusters2.push_back(p);
00327     }
00328 
00329   if(clusters2.size()>0)
00330     std::sort(clusters2.begin(),clusters2.end(),comparePt);
00331   return clusters2;
00332 }
00333
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/RecoTauTag/HLTProducers/src/L2TauModularIsolationProducer.cc
