CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/RecoTauTag/RecoTau/src/CaloRecoTauTagInfoAlgorithm.cc

Go to the documentation of this file.
00001 #include "RecoTauTag/RecoTau/interface/CaloRecoTauTagInfoAlgorithm.h"
00002 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00003 
00004 using namespace reco;
00005 
00006 CaloRecoTauTagInfoAlgorithm::CaloRecoTauTagInfoAlgorithm(const edm::ParameterSet& parameters){
00007   // parameters of the considered rec. Tracks (catched through a JetTracksAssociation object) :
00008   tkminPt_                            = parameters.getParameter<double>("tkminPt");
00009   tkminPixelHitsn_                    = parameters.getParameter<int>("tkminPixelHitsn");
00010   tkminTrackerHitsn_                  = parameters.getParameter<int>("tkminTrackerHitsn");
00011   tkmaxipt_                           = parameters.getParameter<double>("tkmaxipt");
00012   tkmaxChi2_                          = parameters.getParameter<double>("tkmaxChi2");
00013   // 
00014   UsePVconstraint_                    = parameters.getParameter<bool>("UsePVconstraint");
00015   tkPVmaxDZ_                          = parameters.getParameter<double>("tkPVmaxDZ");
00016   //
00017   UseTrackQuality_                    = parameters.getParameter<bool>("UseTrackQuality");
00018   if (UseTrackQuality_) {
00019     tkQuality_ = reco::TrackBase::qualityByName(parameters.getParameter<std::string>("tkQuality"));
00020   }
00021   // parameters of the considered EcalRecHits 
00022   BarrelBasicClusters_                = parameters.getParameter<edm::InputTag>("BarrelBasicClustersSource"); 
00023   EndcapBasicClusters_                = parameters.getParameter<edm::InputTag>("EndcapBasicClustersSource"); 
00024   // parameters of the considered neutral ECAL BasicClusters
00025   ECALBasicClustersAroundCaloJet_DRConeSize_      = parameters.getParameter<double>("ECALBasicClustersAroundCaloJet_DRConeSize");
00026   ECALBasicClusterminE_                           = parameters.getParameter<double>("ECALBasicClusterminE");
00027   ECALBasicClusterpropagTrack_matchingDRConeSize_ = parameters.getParameter<double>("ECALBasicClusterpropagTrack_matchingDRConeSize");
00028 }
00029   
00030 CaloTauTagInfo CaloRecoTauTagInfoAlgorithm::buildCaloTauTagInfo(edm::Event& theEvent,const edm::EventSetup& theEventSetup,const CaloJetRef& theCaloJet,const TrackRefVector& theTracks,const Vertex& thePV){
00031   CaloTauTagInfo resultExtended;
00032   resultExtended.setcalojetRef(theCaloJet);
00033 
00034   TrackRefVector theFilteredTracks;
00035   if (UsePVconstraint_) theFilteredTracks=TauTagTools::filteredTracks(theTracks,tkminPt_,tkminPixelHitsn_,tkminTrackerHitsn_,tkmaxipt_,tkmaxChi2_,tkPVmaxDZ_,thePV, thePV.z());
00036   else theFilteredTracks=TauTagTools::filteredTracks(theTracks,tkminPt_,tkminPixelHitsn_,tkminTrackerHitsn_,tkmaxipt_,tkmaxChi2_,thePV);
00037   if (UseTrackQuality_) theFilteredTracks = filterTracksByQualityBit(theFilteredTracks,tkQuality_);
00038   resultExtended.setTracks(theFilteredTracks);
00039   
00040   //resultExtended.setpositionAndEnergyECALRecHits(getPositionAndEnergyEcalRecHits(theEvent,theEventSetup,theCaloJet));
00041 
00042   std::vector<BasicClusterRef> theNeutralEcalBasicClusters=getNeutralEcalBasicClusters(theEvent,theEventSetup,theCaloJet,theFilteredTracks,ECALBasicClustersAroundCaloJet_DRConeSize_,ECALBasicClusterminE_,ECALBasicClusterpropagTrack_matchingDRConeSize_);
00043   resultExtended.setneutralECALBasicClusters(theNeutralEcalBasicClusters);
00044   
00045   return resultExtended; 
00046 }
00047 CaloTauTagInfo CaloRecoTauTagInfoAlgorithm::buildCaloTauTagInfo(edm::Event& theEvent,const edm::EventSetup& theEventSetup,const JetBaseRef& theJet,const TrackRefVector& theTracks,const Vertex& thePV){
00048   CaloTauTagInfo resultExtended;
00049   resultExtended.setJetRef(theJet);
00050 
00051   TrackRefVector theFilteredTracks;
00052   if (UsePVconstraint_) theFilteredTracks=TauTagTools::filteredTracks(theTracks,tkminPt_,tkminPixelHitsn_,tkminTrackerHitsn_,tkmaxipt_,tkmaxChi2_,tkPVmaxDZ_,thePV, thePV.z());
00053   else theFilteredTracks=TauTagTools::filteredTracks(theTracks,tkminPt_,tkminPixelHitsn_,tkminTrackerHitsn_,tkmaxipt_,tkmaxChi2_,thePV);
00054   if (UseTrackQuality_) theFilteredTracks = filterTracksByQualityBit(theFilteredTracks,tkQuality_);
00055   resultExtended.setTracks(theFilteredTracks);
00056 
00057   //resultExtended.setpositionAndEnergyECALRecHits(getPositionAndEnergyEcalRecHits(theEvent,theEventSetup,theCaloJet));
00058 
00059   //reco::JPTJetRef const theJPTJetRef = theJet.castTo<reco::JPTJetRef>();
00060   //reco::CaloJetRef const theCaloJet = (theJPTJetRef->getCaloJetRef()).castTo<reco::CaloJetRef>();
00061   reco::CaloJetRef const theCaloJet = resultExtended.calojetRef();
00062   std::vector<BasicClusterRef> theNeutralEcalBasicClusters=getNeutralEcalBasicClusters(theEvent,theEventSetup,theCaloJet,theFilteredTracks,ECALBasicClustersAroundCaloJet_DRConeSize_,ECALBasicClusterminE_,ECALBasicClusterpropagTrack_matchingDRConeSize_);
00063   resultExtended.setneutralECALBasicClusters(theNeutralEcalBasicClusters);
00064 
00065   return resultExtended;
00066 }
00067 /*
00068 std::vector<std::pair<math::XYZPoint,float> > CaloRecoTauTagInfoAlgorithm::getPositionAndEnergyEcalRecHits(edm::Event& theEvent,const edm::EventSetup& theEventSetup,const CaloJetRef& theCaloJet){
00069   std::vector<std::pair<math::XYZPoint,float> > thePositionAndEnergyEcalRecHits;
00070   std::vector<CaloTowerPtr> theCaloTowers=theCaloJet->getCaloConstituents();
00071   ESHandle<CaloGeometry> theCaloGeometry;
00072   theEventSetup.get<CaloGeometryRecord>().get(theCaloGeometry);
00073   const CaloSubdetectorGeometry* theCaloSubdetectorGeometry;  
00074   edm::Handle<EBRecHitCollection> EBRecHits;
00075   edm::Handle<EERecHitCollection> EERecHits; 
00076   edm::Handle<ESRecHitCollection> ESRecHits; 
00077   theEvent.getByLabel(EBRecHitsLabel_,EBRecHits);
00078   theEvent.getByLabel(EERecHitsLabel_,EERecHits);
00079   theEvent.getByLabel(ESRecHitsLabel_,ESRecHits);
00080   for(std::vector<CaloTowerPtr>::const_iterator i_Tower=theCaloTowers.begin();i_Tower!=theCaloTowers.end();i_Tower++){
00081     size_t numRecHits = (**i_Tower).constituentsSize();
00082     for(size_t j=0;j<numRecHits;j++) {
00083       DetId RecHitDetID=(**i_Tower).constituent(j);      
00084 
00085 
00086       DetId::Detector DetNum=RecHitDetID.det();     
00087       if(DetNum==DetId::Ecal){
00088         if((EcalSubdetector)RecHitDetID.subdetId()==EcalBarrel){
00089           theCaloSubdetectorGeometry = theCaloGeometry->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
00090           EBDetId EcalID=RecHitDetID;
00091           EBRecHitCollection::const_iterator theRecHit=EBRecHits->find(EcalID);
00092           const CaloCellGeometry* theRecHitCell=theCaloSubdetectorGeometry->getGeometry(RecHitDetID);
00093           math::XYZPoint theRecHitCell_XYZPoint(theRecHitCell->getPosition().x(),theRecHitCell->getPosition().y(),theRecHitCell->getPosition().z());
00094           pair<math::XYZPoint,float> thePositionAndEnergyEcalRecHit(theRecHitCell_XYZPoint,theRecHit->energy());
00095           thePositionAndEnergyEcalRecHits.push_back(thePositionAndEnergyEcalRecHit);
00096         }else if((EcalSubdetector)RecHitDetID.subdetId()==EcalEndcap){
00097           theCaloSubdetectorGeometry = theCaloGeometry->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);
00098           EEDetId EcalID = RecHitDetID;
00099           EERecHitCollection::const_iterator theRecHit=EERecHits->find(EcalID);     
00100           const CaloCellGeometry* theRecHitCell=theCaloSubdetectorGeometry->getGeometry(RecHitDetID);
00101           math::XYZPoint theRecHitCell_XYZPoint(theRecHitCell->getPosition().x(),theRecHitCell->getPosition().y(),theRecHitCell->getPosition().z());
00102           pair<math::XYZPoint,float> thePositionAndEnergyEcalRecHit(theRecHitCell_XYZPoint,theRecHit->energy());
00103           thePositionAndEnergyEcalRecHits.push_back(thePositionAndEnergyEcalRecHit);
00104         }else if((EcalSubdetector)RecHitDetID.subdetId()==EcalPreshower){
00105           theCaloSubdetectorGeometry = theCaloGeometry->getSubdetectorGeometry(DetId::Ecal,EcalPreshower);
00106           ESDetId EcalID = RecHitDetID;
00107           ESRecHitCollection::const_iterator theRecHit=ESRecHits->find(EcalID);     
00108           const CaloCellGeometry* theRecHitCell=theCaloSubdetectorGeometry->getGeometry(RecHitDetID);
00109           math::XYZPoint theRecHitCell_XYZPoint(theRecHitCell->getPosition().x(),theRecHitCell->getPosition().y(),theRecHitCell->getPosition().z());
00110           pair<math::XYZPoint,float> thePositionAndEnergyEcalRecHit(theRecHitCell_XYZPoint,theRecHit->energy());
00111           thePositionAndEnergyEcalRecHits.push_back(thePositionAndEnergyEcalRecHit);
00112         }        
00113       } 
00114     }
00115   }
00116   return thePositionAndEnergyEcalRecHits;
00117 }
00118 */
00119 
00120 std::vector<DetId> CaloRecoTauTagInfoAlgorithm::getVectorDetId(const CaloJetRef& theCaloJet){
00121   std::vector<CaloTowerPtr> theCaloTowers=theCaloJet->getCaloConstituents();
00122   std::vector<DetId> myDetIds;
00123   myDetIds.clear();
00124 
00125   for(std::vector<CaloTowerPtr>::const_iterator i_Tower=theCaloTowers.begin();i_Tower!=theCaloTowers.end();i_Tower++){
00126     size_t numRecHits = (**i_Tower).constituentsSize();
00127     for(size_t j=0;j<numRecHits;j++) {
00128       DetId RecHitDetID=(**i_Tower).constituent(j);      
00129 
00130       myDetIds.push_back(RecHitDetID);
00131     }
00132   }
00133   return myDetIds;
00134 }
00135 
00136 
00137 std::vector<BasicClusterRef> CaloRecoTauTagInfoAlgorithm::getNeutralEcalBasicClusters(edm::Event& theEvent,const edm::EventSetup& theEventSetup,const CaloJetRef& theCaloJet,const TrackRefVector& theTracks,float theECALBasicClustersAroundCaloJet_DRConeSize,float theECALBasicClusterminE,float theECALBasicClusterpropagTrack_matchingDRConeSize){
00138   std::vector<math::XYZPoint> thepropagTracksECALSurfContactPoints;
00139   edm::ESHandle<MagneticField> theMF;
00140   theEventSetup.get<IdealMagneticFieldRecord>().get(theMF);
00141   const MagneticField* theMagField=theMF.product();
00142   for(TrackRefVector::const_iterator i_Track=theTracks.begin();i_Track!=theTracks.end();i_Track++){
00143     math::XYZPoint thepropagTrackECALSurfContactPoint=TauTagTools::propagTrackECALSurfContactPoint(theMagField,*i_Track);
00144     if(thepropagTrackECALSurfContactPoint.R()!=0.) thepropagTracksECALSurfContactPoints.push_back(thepropagTrackECALSurfContactPoint);
00145   }
00146   
00147   math::XYZPoint aCaloJetFakePosition((*theCaloJet).px(),(*theCaloJet).py(),(*theCaloJet).pz());
00148     
00149   std::vector<BasicClusterRef> theBasicClusters; 
00150   
00151   edm::Handle<BasicClusterCollection> theBarrelBCCollection;
00152   //  theEvent.getByLabel("islandBasicClusters","islandBarrelBasicClusters",theBarrelBCCollection);
00153   theEvent.getByLabel(BarrelBasicClusters_,theBarrelBCCollection);
00154   for(unsigned int i_BC=0;i_BC!=theBarrelBCCollection->size();i_BC++) { 
00155     BasicClusterRef theBasicClusterRef(theBarrelBCCollection,i_BC);    
00156     if (theBasicClusterRef.isNull()) continue;  
00157     if (ROOT::Math::VectorUtil::DeltaR(aCaloJetFakePosition,(*theBasicClusterRef).position())<=theECALBasicClustersAroundCaloJet_DRConeSize && (*theBasicClusterRef).energy()>=theECALBasicClusterminE) theBasicClusters.push_back(theBasicClusterRef);
00158   }
00159   edm::Handle<BasicClusterCollection> theEndcapBCCollection;
00160   //  theEvent.getByLabel("islandBasicClusters","islandEndcapBasicClusters",theEndcapBCCollection);
00161   theEvent.getByLabel(EndcapBasicClusters_,theEndcapBCCollection);
00162   for(unsigned int j_BC=0;j_BC!=theEndcapBCCollection->size();j_BC++) { 
00163     BasicClusterRef theBasicClusterRef(theEndcapBCCollection,j_BC); 
00164     if (theBasicClusterRef.isNull()) continue;  
00165     if (ROOT::Math::VectorUtil::DeltaR(aCaloJetFakePosition,(*theBasicClusterRef).position())<=theECALBasicClustersAroundCaloJet_DRConeSize && (*theBasicClusterRef).energy()>=theECALBasicClusterminE) theBasicClusters.push_back(theBasicClusterRef);
00166   }  
00167   
00168   std::vector<BasicClusterRef> theNeutralBasicClusters=theBasicClusters;  
00169   std::vector<BasicClusterRef>::iterator kmatchedBasicCluster;
00170   for (std::vector<math::XYZPoint>::iterator ipropagTrackECALSurfContactPoint=thepropagTracksECALSurfContactPoints.begin();ipropagTrackECALSurfContactPoint!=thepropagTracksECALSurfContactPoints.end();ipropagTrackECALSurfContactPoint++) {
00171     double theMatchedEcalBasicClusterpropagTrack_minDR=theECALBasicClusterpropagTrack_matchingDRConeSize;
00172     bool Track_matchedwithEcalBasicCluster=false;
00173     for (std::vector<BasicClusterRef>::iterator jBasicCluster=theNeutralBasicClusters.begin();jBasicCluster!=theNeutralBasicClusters.end();jBasicCluster++) {
00174       if(ROOT::Math::VectorUtil::DeltaR((*ipropagTrackECALSurfContactPoint),(**jBasicCluster).position())<theMatchedEcalBasicClusterpropagTrack_minDR){
00175         Track_matchedwithEcalBasicCluster=true;
00176         theMatchedEcalBasicClusterpropagTrack_minDR=ROOT::Math::VectorUtil::DeltaR((*ipropagTrackECALSurfContactPoint),(**jBasicCluster).position());
00177         kmatchedBasicCluster=jBasicCluster;
00178       }
00179     }
00180     if(Track_matchedwithEcalBasicCluster) kmatchedBasicCluster=theNeutralBasicClusters.erase(kmatchedBasicCluster);
00181   }
00182   return theNeutralBasicClusters;
00183 }
00184 
00185 TrackRefVector CaloRecoTauTagInfoAlgorithm::filterTracksByQualityBit(const TrackRefVector& tracks, reco::TrackBase::TrackQuality quality) const
00186 {
00187   TrackRefVector filteredTracks;
00188   for (TrackRefVector::const_iterator iTrack = tracks.begin(); iTrack != tracks.end(); iTrack++) {
00189     if ((*iTrack)->quality(quality)) filteredTracks.push_back(*iTrack);
00190   }
00191   return filteredTracks;
00192 }