00001 #include "RecoTauTag/TauTagTools/interface/TauTagTools.h"
00002
00003 using namespace reco;
00004 using std::string;
00005
00006 namespace TauTagTools{
00007
00008 double computeDeltaR(const math::XYZVector& vec1, const math::XYZVector& vec2)
00009 {
00010 DeltaR<math::XYZVector> myMetricDR_;
00011 return myMetricDR_(vec1, vec2);
00012 }
00013 double computeAngle(const math::XYZVector& vec1, const math::XYZVector& vec2)
00014 {
00015 Angle<math::XYZVector> myMetricAngle_;
00016 return myMetricAngle_(vec1, vec2);
00017 }
00018 TFormula computeConeSizeTFormula(const string& ConeSizeFormula,const char* errorMessage){
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 TFormula ConeSizeTFormula;
00029 string ConeSizeFormulaStr = ConeSizeFormula;
00030 replaceSubStr(ConeSizeFormulaStr,"JetOpeningDR", "z");
00031 replaceSubStr(ConeSizeFormulaStr,"ET","y");
00032 replaceSubStr(ConeSizeFormulaStr,"E","x");
00033 ConeSizeTFormula.SetName("ConeSize");
00034 ConeSizeTFormula.SetTitle(ConeSizeFormulaStr.data());
00035 int errorFlag = ConeSizeTFormula.Compile();
00036 if (errorFlag!= 0) {
00037 throw cms::Exception("") << "\n unsupported functional Form for " << errorMessage << " " << ConeSizeFormula << std::endl
00038 << "Please check that the Definition in \"" << ConeSizeTFormula.GetName() << "\" only contains the variables \"E\" or \"ET\""
00039 << " and Functions that are supported by ROOT's TFormular Class." << std::endl;
00040 }else return ConeSizeTFormula;
00041 }
00042 void replaceSubStr(string& s,const string& oldSubStr,const string& newSubStr){
00043
00044
00045
00046 if ( oldSubStr == newSubStr ) return;
00047
00048
00049
00050
00051 if ( oldSubStr.empty() ) return;
00052
00053 const string::size_type lengthOldSubStr = oldSubStr.size();
00054 const string::size_type lengthNewSubStr = newSubStr.size();
00055
00056 string::size_type positionPreviousMatch = 0;
00057 string::size_type positionNextMatch = 0;
00058
00059
00060
00061 while ( (positionNextMatch = s.find(oldSubStr, positionPreviousMatch)) != string::npos ) {
00062 s.replace(positionNextMatch, lengthOldSubStr, newSubStr);
00063 positionPreviousMatch = positionNextMatch + lengthNewSubStr;
00064 }
00065 }
00066
00067
00068 TrackRefVector filteredTracksByNumTrkHits(TrackRefVector theInitialTracks, int tkminTrackerHitsn){
00069 TrackRefVector filteredTracks;
00070 for(TrackRefVector::const_iterator iTk=theInitialTracks.begin();iTk!=theInitialTracks.end();iTk++){
00071 if ( (**iTk).numberOfValidHits() >= tkminTrackerHitsn )
00072 filteredTracks.push_back(*iTk);
00073 }
00074 return filteredTracks;
00075 }
00076
00077 TrackRefVector filteredTracks(TrackRefVector theInitialTracks,double tkminPt,int tkminPixelHitsn,int tkminTrackerHitsn,double tkmaxipt,double tkmaxChi2, Vertex pv){
00078 TrackRefVector filteredTracks;
00079 for(TrackRefVector::const_iterator iTk=theInitialTracks.begin();iTk!=theInitialTracks.end();iTk++){
00080 if ((**iTk).pt()>=tkminPt &&
00081 (**iTk).normalizedChi2()<=tkmaxChi2 &&
00082 fabs((**iTk).dxy(pv.position()))<=tkmaxipt &&
00083 (**iTk).numberOfValidHits()>=tkminTrackerHitsn &&
00084 (**iTk).hitPattern().numberOfValidPixelHits()>=tkminPixelHitsn)
00085 filteredTracks.push_back(*iTk);
00086 }
00087 return filteredTracks;
00088 }
00089 TrackRefVector filteredTracks(TrackRefVector theInitialTracks,double tkminPt,int tkminPixelHitsn,int tkminTrackerHitsn,double tkmaxipt,double tkmaxChi2,double tktorefpointmaxDZ,Vertex pv, double refpoint_Z){
00090 TrackRefVector filteredTracks;
00091 for(TrackRefVector::const_iterator iTk=theInitialTracks.begin();iTk!=theInitialTracks.end();iTk++){
00092 if(pv.isFake()) tktorefpointmaxDZ=30.;
00093 if ((**iTk).pt()>=tkminPt &&
00094 (**iTk).normalizedChi2()<=tkmaxChi2 &&
00095 fabs((**iTk).dxy(pv.position()))<=tkmaxipt &&
00096 (**iTk).numberOfValidHits()>=tkminTrackerHitsn &&
00097 (**iTk).hitPattern().numberOfValidPixelHits()>=tkminPixelHitsn &&
00098 fabs((**iTk).dz(pv.position()))<=tktorefpointmaxDZ)
00099 filteredTracks.push_back(*iTk);
00100 }
00101 return filteredTracks;
00102 }
00103
00104 PFCandidateRefVector filteredPFChargedHadrCandsByNumTrkHits(PFCandidateRefVector theInitialPFCands, int ChargedHadrCand_tkminTrackerHitsn){
00105 PFCandidateRefVector filteredPFChargedHadrCands;
00106 for(PFCandidateRefVector::const_iterator iPFCand=theInitialPFCands.begin();iPFCand!=theInitialPFCands.end();iPFCand++){
00107 if (PFCandidate::ParticleType((**iPFCand).particleId())==PFCandidate::h || PFCandidate::ParticleType((**iPFCand).particleId())==PFCandidate::mu || PFCandidate::ParticleType((**iPFCand).particleId())==PFCandidate::e){
00108
00109 TrackRef PFChargedHadrCand_rectk = (**iPFCand).trackRef();
00110
00111 if (!PFChargedHadrCand_rectk)continue;
00112 if ( (*PFChargedHadrCand_rectk).numberOfValidHits()>=ChargedHadrCand_tkminTrackerHitsn )
00113 filteredPFChargedHadrCands.push_back(*iPFCand);
00114 }
00115 }
00116 return filteredPFChargedHadrCands;
00117 }
00118
00119 PFCandidateRefVector filteredPFChargedHadrCands(PFCandidateRefVector theInitialPFCands,double ChargedHadrCand_tkminPt,int ChargedHadrCand_tkminPixelHitsn,int ChargedHadrCand_tkminTrackerHitsn,double ChargedHadrCand_tkmaxipt,double ChargedHadrCand_tkmaxChi2, Vertex pv){
00120 PFCandidateRefVector filteredPFChargedHadrCands;
00121 for(PFCandidateRefVector::const_iterator iPFCand=theInitialPFCands.begin();iPFCand!=theInitialPFCands.end();iPFCand++){
00122 if (PFCandidate::ParticleType((**iPFCand).particleId())==PFCandidate::h || PFCandidate::ParticleType((**iPFCand).particleId())==PFCandidate::mu || PFCandidate::ParticleType((**iPFCand).particleId())==PFCandidate::e){
00123
00124 TrackRef PFChargedHadrCand_rectk = (**iPFCand).trackRef();
00125
00126
00127 if (!PFChargedHadrCand_rectk)continue;
00128 if ((*PFChargedHadrCand_rectk).pt()>=ChargedHadrCand_tkminPt &&
00129 (*PFChargedHadrCand_rectk).normalizedChi2()<=ChargedHadrCand_tkmaxChi2 &&
00130 fabs((*PFChargedHadrCand_rectk).dxy(pv.position()))<=ChargedHadrCand_tkmaxipt &&
00131 (*PFChargedHadrCand_rectk).numberOfValidHits()>=ChargedHadrCand_tkminTrackerHitsn &&
00132 (*PFChargedHadrCand_rectk).hitPattern().numberOfValidPixelHits()>=ChargedHadrCand_tkminPixelHitsn)
00133 filteredPFChargedHadrCands.push_back(*iPFCand);
00134 }
00135 }
00136 return filteredPFChargedHadrCands;
00137 }
00138 PFCandidateRefVector filteredPFChargedHadrCands(PFCandidateRefVector theInitialPFCands,double ChargedHadrCand_tkminPt,int ChargedHadrCand_tkminPixelHitsn,int ChargedHadrCand_tkminTrackerHitsn,double ChargedHadrCand_tkmaxipt,double ChargedHadrCand_tkmaxChi2,double ChargedHadrCand_tktorefpointmaxDZ,Vertex pv, double refpoint_Z){
00139 if(pv.isFake()) ChargedHadrCand_tktorefpointmaxDZ = 30.;
00140 PFCandidateRefVector filteredPFChargedHadrCands;
00141 for(PFCandidateRefVector::const_iterator iPFCand=theInitialPFCands.begin();iPFCand!=theInitialPFCands.end();iPFCand++){
00142 if (PFCandidate::ParticleType((**iPFCand).particleId())==PFCandidate::h || PFCandidate::ParticleType((**iPFCand).particleId())==PFCandidate::mu || PFCandidate::ParticleType((**iPFCand).particleId())==PFCandidate::e){
00143
00144 TrackRef PFChargedHadrCand_rectk = (**iPFCand).trackRef();
00145 if (!PFChargedHadrCand_rectk)continue;
00146 if ((*PFChargedHadrCand_rectk).pt()>=ChargedHadrCand_tkminPt &&
00147 (*PFChargedHadrCand_rectk).normalizedChi2()<=ChargedHadrCand_tkmaxChi2 &&
00148 fabs((*PFChargedHadrCand_rectk).dxy(pv.position()))<=ChargedHadrCand_tkmaxipt &&
00149 (*PFChargedHadrCand_rectk).numberOfValidHits()>=ChargedHadrCand_tkminTrackerHitsn &&
00150 (*PFChargedHadrCand_rectk).hitPattern().numberOfValidPixelHits()>=ChargedHadrCand_tkminPixelHitsn &&
00151 fabs((*PFChargedHadrCand_rectk).dz(pv.position()))<=ChargedHadrCand_tktorefpointmaxDZ)
00152 filteredPFChargedHadrCands.push_back(*iPFCand);
00153 }
00154 }
00155 return filteredPFChargedHadrCands;
00156 }
00157
00158 PFCandidateRefVector filteredPFNeutrHadrCands(PFCandidateRefVector theInitialPFCands,double NeutrHadrCand_HcalclusMinEt){
00159 PFCandidateRefVector filteredPFNeutrHadrCands;
00160 for(PFCandidateRefVector::const_iterator iPFCand=theInitialPFCands.begin();iPFCand!=theInitialPFCands.end();iPFCand++){
00161 if (PFCandidate::ParticleType((**iPFCand).particleId())==PFCandidate::h0){
00162
00163 if ((**iPFCand).et()>=NeutrHadrCand_HcalclusMinEt){
00164 filteredPFNeutrHadrCands.push_back(*iPFCand);
00165 }
00166 }
00167 }
00168 return filteredPFNeutrHadrCands;
00169 }
00170
00171 PFCandidateRefVector filteredPFGammaCands(PFCandidateRefVector theInitialPFCands,double GammaCand_EcalclusMinEt){
00172 PFCandidateRefVector filteredPFGammaCands;
00173 for(PFCandidateRefVector::const_iterator iPFCand=theInitialPFCands.begin();iPFCand!=theInitialPFCands.end();iPFCand++){
00174 if (PFCandidate::ParticleType((**iPFCand).particleId())==PFCandidate::gamma){
00175
00176 if ((**iPFCand).et()>=GammaCand_EcalclusMinEt){
00177 filteredPFGammaCands.push_back(*iPFCand);
00178 }
00179 }
00180 }
00181 return filteredPFGammaCands;
00182 }
00183
00184 math::XYZPoint propagTrackECALSurfContactPoint(const MagneticField* theMagField,TrackRef theTrack){
00185 AnalyticalPropagator thefwdPropagator(theMagField,alongMomentum);
00186 math::XYZPoint propTrack_XYZPoint(0.,0.,0.);
00187
00188
00189 GlobalVector theTrack_initialGV(0.,0.,0.);
00190 GlobalPoint theTrack_initialGP(0.,0.,0.);
00191 if(theTrack->outerOk()){
00192 GlobalVector theTrack_initialoutermostGV(theTrack->outerMomentum().x(),theTrack->outerMomentum().y(),theTrack->outerMomentum().z());
00193 GlobalPoint theTrack_initialoutermostGP(theTrack->outerPosition().x(),theTrack->outerPosition().y(),theTrack->outerPosition().z());
00194 theTrack_initialGV=theTrack_initialoutermostGV;
00195 theTrack_initialGP=theTrack_initialoutermostGP;
00196 } else if(theTrack->innerOk()){
00197 GlobalVector theTrack_initialinnermostGV(theTrack->innerMomentum().x(),theTrack->innerMomentum().y(),theTrack->innerMomentum().z());
00198 GlobalPoint theTrack_initialinnermostGP(theTrack->innerPosition().x(),theTrack->innerPosition().y(),theTrack->innerPosition().z());
00199 theTrack_initialGV=theTrack_initialinnermostGV;
00200 theTrack_initialGP=theTrack_initialinnermostGP;
00201 } else return (propTrack_XYZPoint);
00202 GlobalTrajectoryParameters theTrack_initialGTPs(theTrack_initialGP,theTrack_initialGV,theTrack->charge(),&*theMagField);
00203
00204
00205
00206 CLHEP::HepSymMatrix covM_HepSM(6,1); covM_HepSM*=1e-6;
00207 CartesianTrajectoryError cov_CTE(covM_HepSM);
00208 FreeTrajectoryState Track_initialFTS(theTrack_initialGTPs,cov_CTE);
00209
00210
00211
00212 double ECALcorner_tantheta=ECALBounds::barrel_innerradius()/ECALBounds::barrel_halfLength();
00213 TrajectoryStateOnSurface Track_propagatedonECAL_TSOS=thefwdPropagator.propagate((Track_initialFTS),ECALBounds::barrelBound());
00214 if(!Track_propagatedonECAL_TSOS.isValid() || fabs(Track_propagatedonECAL_TSOS.globalParameters().position().perp()/Track_propagatedonECAL_TSOS.globalParameters().position().z())<ECALcorner_tantheta) {
00215 if(Track_propagatedonECAL_TSOS.isValid() && fabs(Track_propagatedonECAL_TSOS.globalParameters().position().perp()/Track_propagatedonECAL_TSOS.globalParameters().position().z())<ECALcorner_tantheta){
00216 if(Track_propagatedonECAL_TSOS.globalParameters().position().eta()>0.){
00217 Track_propagatedonECAL_TSOS=thefwdPropagator.propagate((Track_initialFTS),ECALBounds::positiveEndcapDisk());
00218 }else{
00219 Track_propagatedonECAL_TSOS=thefwdPropagator.propagate((Track_initialFTS),ECALBounds::negativeEndcapDisk());
00220 }
00221 }
00222 if(!Track_propagatedonECAL_TSOS.isValid()){
00223 if((Track_initialFTS).position().eta()>0.){
00224 Track_propagatedonECAL_TSOS=thefwdPropagator.propagate((Track_initialFTS),ECALBounds::positiveEndcapDisk());
00225 }else{
00226 Track_propagatedonECAL_TSOS=thefwdPropagator.propagate((Track_initialFTS),ECALBounds::negativeEndcapDisk());
00227 }
00228 }
00229 }
00230 if(Track_propagatedonECAL_TSOS.isValid()){
00231 math::XYZPoint validpropTrack_XYZPoint(Track_propagatedonECAL_TSOS.globalPosition().x(),
00232 Track_propagatedonECAL_TSOS.globalPosition().y(),
00233 Track_propagatedonECAL_TSOS.globalPosition().z());
00234 propTrack_XYZPoint=validpropTrack_XYZPoint;
00235 }
00236 return (propTrack_XYZPoint);
00237 }
00238
00239
00240
00241 void
00242 sortRefVectorByPt(PFCandidateRefVector& vec)
00243 {
00244
00245 std::vector<size_t> indices;
00246 indices.reserve(vec.size());
00247 for(unsigned int i=0;i<vec.size();++i)
00248 indices.push_back(i);
00249
00250 refVectorPtSorter sorter(vec);
00251 std::sort(indices.begin(),indices.end(),sorter);
00252
00253
00254 reco::PFCandidateRefVector sorted;
00255 sorted.reserve(vec.size());
00256
00257 for(unsigned int i=0;i<indices.size();++i)
00258 sorted.push_back(vec.at(indices.at(i)));
00259
00260 vec = sorted;
00261 }
00262
00263
00264 }