CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/Calibration/IsolatedParticles/src/ChargeIsolation.cc

Go to the documentation of this file.
00001 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
00002 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
00003 #include "Calibration/IsolatedParticles/interface/FindDistCone.h"
00004 #include "Calibration/IsolatedParticles/interface/MatrixECALDetIds.h"
00005 #include "Calibration/IsolatedParticles/interface/MatrixHCALDetIds.h"
00006 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00007 
00008 #include<iostream>
00009 
00010 namespace spr{
00011 
00012 
00013   double chargeIsolationEcal(unsigned int trkIndex, std::vector<spr::propagatedTrackID>& vdetIds, const CaloGeometry* geo, const CaloTopology* caloTopology, int ieta, int iphi, bool debug) {
00014   
00015     const DetId coreDet = vdetIds[trkIndex].detIdECAL;
00016     if (debug) {
00017       if (coreDet.subdetId() == EcalBarrel) 
00018         std::cout << "DetId " << (EBDetId)(coreDet) << " Flag " << vdetIds[trkIndex].okECAL << std::endl;
00019       else
00020         std::cout << "DetId " << (EEDetId)(coreDet) << " Flag " << vdetIds[trkIndex].okECAL << std::endl;
00021     }
00022 
00023     double maxNearP = -1.0;
00024     if (vdetIds[trkIndex].okECAL) {
00025       std::vector<DetId> vdets = spr::matrixECALIds(coreDet, ieta, iphi, geo, caloTopology, debug);
00026       if (debug) std::cout << "chargeIsolationEcal:: eta/phi/dets " << ieta << " " << iphi << " " << vdets.size() << std::endl;
00027 
00028       for (unsigned int indx=0; indx<vdetIds.size(); ++indx) {
00029         if (indx != trkIndex && vdetIds[indx].ok && vdetIds[indx].okECAL) {
00030           const DetId anyCell = vdetIds[indx].detIdECAL;
00031           if (!spr::chargeIsolation(anyCell,vdets)) {
00032             const reco::Track* pTrack = &(*(vdetIds[indx].trkItr));
00033             if (maxNearP < pTrack->p()) maxNearP = pTrack->p();
00034           }
00035         }
00036       }
00037     }
00038     return maxNearP;
00039   }
00040 
00041   //===========================================================================================================
00042 
00043   double chargeIsolationEcal(const DetId& coreDet, reco::TrackCollection::const_iterator trkItr, edm::Handle<reco::TrackCollection> trkCollection, const CaloGeometry* geo, const CaloTopology* caloTopology, const MagneticField* bField, int ieta, int iphi, std::string& theTrackQuality, bool debug) {
00044   
00045     const EcalBarrelGeometry *barrelGeom = (dynamic_cast< const EcalBarrelGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel)));
00046     const EcalEndcapGeometry *endcapGeom = (dynamic_cast< const EcalEndcapGeometry *> (geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap)));
00047 
00048     std::vector<DetId> vdets = spr::matrixECALIds(coreDet, ieta, iphi, geo, caloTopology, debug);
00049     if (debug) std::cout << "chargeIsolation:: eta/phi/dets " << ieta << " " << iphi << " " << vdets.size() << std::endl;
00050 
00051     double maxNearP = -1.0;
00052     reco::TrackBase::TrackQuality trackQuality_=reco::TrackBase::qualityByName(theTrackQuality);
00053 
00054     // const DetId anyCell,
00055     reco::TrackCollection::const_iterator trkItr2;
00056     for (trkItr2 = trkCollection->begin(); trkItr2 != trkCollection->end(); ++trkItr2) {
00057 
00058       const reco::Track* pTrack2 = &(*trkItr2);
00059 
00060       bool   trkQuality  = pTrack2->quality(trackQuality_);
00061       if ( (trkItr2 != trkItr) && trkQuality )  {
00062       
00063         std::pair<math::XYZPoint,bool> info = spr::propagateECAL(pTrack2,bField);
00064         const GlobalPoint point2(info.first.x(),info.first.y(),info.first.z());
00065 
00066         if (info.second) {
00067           if (std::abs(point2.eta())<1.479) {
00068             const DetId anyCell = barrelGeom->getClosestCell(point2);
00069             if (!spr::chargeIsolation(anyCell,vdets)) {
00070               if (maxNearP<pTrack2->p()) maxNearP=pTrack2->p();
00071             }
00072           } else {
00073             const DetId anyCell = endcapGeom->getClosestCell(point2);
00074             if (!spr::chargeIsolation(anyCell,vdets)) {
00075               if (maxNearP<pTrack2->p()) maxNearP=pTrack2->p();
00076             }
00077           }
00078         } //info.second
00079       }
00080     }
00081     return maxNearP;
00082   }
00083 
00084   //===========================================================================================================
00085 
00086   double chargeIsolationHcal(unsigned int trkIndex, std::vector<spr::propagatedTrackID> & vdetIds, const HcalTopology* topology, int ieta, int iphi, bool debug) {
00087   
00088     std::vector<DetId> dets(1,vdetIds[trkIndex].detIdHCAL);
00089     if (debug) {
00090       std::cout << "DetId " << (HcalDetId)(dets[0]) << " Flag " << vdetIds[trkIndex].okHCAL << std::endl;
00091     }
00092 
00093     double maxNearP = -1.0;
00094     if (vdetIds[trkIndex].okHCAL) {
00095       std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, false, debug);
00096       if (debug) std::cout << "chargeIsolationHcal:: eta/phi/dets " << ieta << " " << iphi << " " << vdets.size() << std::endl;
00097 
00098       for (unsigned indx = 0; indx<vdetIds.size(); ++indx) {
00099         if (indx != trkIndex && vdetIds[indx].ok && vdetIds[indx].okHCAL) {
00100           const DetId anyCell = vdetIds[indx].detIdHCAL;
00101           if (!spr::chargeIsolation(anyCell,vdets)) {
00102             const reco::Track* pTrack = &(*(vdetIds[indx].trkItr));
00103             if (maxNearP < pTrack->p()) maxNearP = pTrack->p();
00104           }
00105         }
00106       }
00107     }
00108     return maxNearP;
00109   }
00110 
00111   //===========================================================================================================
00112 
00113   double chargeIsolationHcal(reco::TrackCollection::const_iterator trkItr, edm::Handle<reco::TrackCollection> trkCollection, const DetId ClosestCell, const HcalTopology* topology, const CaloSubdetectorGeometry* gHB, const MagneticField* bField, int ieta, int iphi, std::string& theTrackQuality, bool debug) {
00114 
00115     std::vector<DetId> dets(1,ClosestCell);
00116 
00117     if (debug) std::cout << (HcalDetId) ClosestCell << std::endl;
00118 
00119     std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, false, debug);
00120     std::vector<DetId>::iterator it;  
00121   
00122     if (debug) {
00123       for (unsigned int i=0; i<vdets.size(); i++) {
00124         std::cout << "HcalDetId in " <<2*ieta+1 << "x" << 2*iphi+1 << " " << (HcalDetId) vdets[i] << std::endl;
00125       }
00126     }
00127 
00128     double maxNearP = -1.0;
00129     reco::TrackBase::TrackQuality trackQuality_=reco::TrackBase::qualityByName(theTrackQuality);
00130   
00131     reco::TrackCollection::const_iterator trkItr2;
00132     for (trkItr2 = trkCollection->begin(); trkItr2 != trkCollection->end(); ++trkItr2) {
00133     
00134       const reco::Track* pTrack2 = &(*trkItr2);
00135     
00136       bool   trkQuality  = pTrack2->quality(trackQuality_);
00137       if ( (trkItr2 != trkItr) && trkQuality )  {
00138         std::pair<math::XYZPoint,bool> info = spr::propagateHCAL(pTrack2,bField);
00139         const GlobalPoint point2(info.first.x(),info.first.y(),info.first.z());
00140 
00141         if (debug) {
00142           std::cout << "Track2 (p,eta,phi) " << pTrack2->p() << " " << pTrack2->eta() << " " << pTrack2->phi() << std::endl;
00143         }
00144 
00145         if (info.second) {
00146           const DetId anyCell = gHB->getClosestCell(point2);
00147           if (!spr::chargeIsolation(anyCell,vdets)) {   
00148             if(maxNearP<pTrack2->p())  maxNearP=pTrack2->p();
00149           }
00150           if (debug){
00151             std::cout << "maxNearP " << maxNearP << " thisCell " 
00152                       << (HcalDetId)anyCell << " (" 
00153                       << info.first.x() << "," << info.first.y() <<","
00154                       << info.first.z() << ")" << std::endl;
00155           }
00156         }
00157       }
00158     }
00159     return maxNearP;
00160   }
00161 
00162   //===========================================================================================================
00163 
00164   bool chargeIsolation(const DetId anyCell, std::vector<DetId>& vdets) {
00165     bool isIsolated = true;
00166     for (unsigned int i=0; i<vdets.size(); i++){
00167       if (anyCell == vdets[i] ) {
00168         isIsolated = false;
00169         break;
00170       }
00171     }
00172     return isIsolated;
00173   }
00174 
00175   //===========================================================================================================
00176 
00177   double coneChargeIsolation(const edm::Event& iEvent, const edm::EventSetup& iSetup, reco::TrackCollection::const_iterator trkItr, edm::Handle<reco::TrackCollection> trkCollection, TrackDetectorAssociator& associator, TrackAssociatorParameters& parameters_, std::string theTrackQuality, int &nNearTRKs, int &nLayers_maxNearP, int &trkQual_maxNearP, double &maxNearP_goodTrk, const GlobalPoint& hpoint1, const GlobalVector& trackMom, double dR) {
00178 
00179     nNearTRKs=0;
00180     nLayers_maxNearP=0;
00181     trkQual_maxNearP=-1; 
00182     maxNearP_goodTrk = -999.0;
00183     double maxNearP = -999.0;
00184     reco::TrackBase::TrackQuality trackQuality_=  reco::TrackBase::qualityByName(theTrackQuality);
00185 
00186     // Iterate over tracks
00187     reco::TrackCollection::const_iterator trkItr2;
00188     for( trkItr2 = trkCollection->begin(); 
00189          trkItr2 != trkCollection->end(); ++trkItr2){
00190 
00191       // Get track
00192       const reco::Track* pTrack2 = &(*trkItr2);
00193     
00194       // Get track qual, nlayers, and hit pattern
00195       if (pTrack2->quality(trackQuality_)) trkQual_maxNearP  = 1;
00196       const reco::HitPattern& hitp = pTrack2->hitPattern();
00197       nLayers_maxNearP = hitp.trackerLayersWithMeasurement() ;        
00198     
00199       // Skip if the neighboring track candidate is the iso-track
00200       // candidate
00201       if (trkItr2 != trkItr) {
00202     
00203         // Get propagator
00204         const FreeTrajectoryState fts2 = associator.getFreeTrajectoryState(iSetup, *pTrack2);
00205         TrackDetMatchInfo info2 = associator.associate(iEvent, iSetup, fts2, parameters_);
00206     
00207         // Make sure it reaches Hcal
00208         if (info2.isGoodHcal ) {
00209     
00210           const GlobalPoint point2(info2.trkGlobPosAtHcal.x(),
00211                                    info2.trkGlobPosAtHcal.y(),
00212                                    info2.trkGlobPosAtHcal.z());
00213     
00214           int isConeChargedIso = spr::coneChargeIsolation(hpoint1, point2, trackMom, dR);
00215     
00216           if (isConeChargedIso==0) {
00217             nNearTRKs++;
00218             if(maxNearP<pTrack2->p()) {
00219               maxNearP=pTrack2->p();
00220               if (trkQual_maxNearP>0 && nLayers_maxNearP>7 && maxNearP_goodTrk<pTrack2->p()) {
00221                 maxNearP_goodTrk=pTrack2->p();
00222               }
00223             }
00224           }
00225         }
00226       }
00227     } // Iterate over track loop
00228     
00229     return maxNearP;
00230   }
00231 
00232   double chargeIsolationCone(unsigned int trkIndex, std::vector<spr::propagatedTrackDirection> & trkDirs, double dR, int & nNearTRKs, bool debug) {
00233 
00234     double maxNearP = -1.0;
00235     nNearTRKs = 0;
00236     if (trkDirs[trkIndex].okHCAL) {
00237       for (unsigned int indx=0; indx<trkDirs.size(); ++indx) {
00238         if (indx != trkIndex && trkDirs[indx].ok && trkDirs[indx].okHCAL) {
00239           int isConeChargedIso = spr::coneChargeIsolation(trkDirs[trkIndex].pointHCAL, trkDirs[indx].pointHCAL, trkDirs[trkIndex].directionHCAL, dR);
00240           if (isConeChargedIso==0) {
00241             nNearTRKs++;
00242             const reco::Track* pTrack = &(*(trkDirs[indx].trkItr));
00243             if (maxNearP < pTrack->p()) maxNearP = pTrack->p();
00244           }
00245         }
00246       }
00247     }
00248     return maxNearP;
00249   }
00250 
00251 
00252   int coneChargeIsolation(const GlobalPoint& hpoint1, const GlobalPoint& point2, const GlobalVector& trackMom, double dR) {                      
00253 
00254     int isIsolated = 1;
00255     if (spr::getDistInPlaneTrackDir(hpoint1, trackMom, point2) > dR) isIsolated = 1;
00256     else isIsolated = 0;
00257     return isIsolated;
00258   } 
00259 
00260 }