CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/RecoTauTag/TauTagTools/src/TauTagTools.cc

Go to the documentation of this file.
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   //--- check functional form 
00020   //    given as configuration parameter for matching and signal cone sizes;
00021   //
00022   //    The size of a cone may depend on the energy "E" and/or transverse energy "ET" of the tau-jet candidate.
00023   //    Alternatively one can additionally use "JetOpeningDR", which specifies the opening angle of the seed jet.
00024   //    Any functional form that is supported by ROOT's TFormula class can be used (e.g. "3.0/E", "0.25/sqrt(ET)")
00025   //
00026   //    replace "E"  by TFormula variable "x"
00027   //            "ET"                      "y"
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()); // the function definition is actually stored in the "Title" data-member of the TFormula object
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   //--- protect replacement algorithm
00044   //    from case that oldSubStr and newSubStr are equal
00045   //    (nothing to be done anyway)
00046   if ( oldSubStr == newSubStr ) return;
00047   
00048   //--- protect replacement algorithm
00049   //    from case that oldSubStr contains no characters
00050   //    (i.e. matches everything)
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   //--- consecutively replace all occurences of oldSubStr by newSubStr;
00060   //    keep iterating until no occurence of oldSubStr left
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         // *** Whether the charged hadron candidate will be selected or not depends on its rec. tk properties. 
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         // *** Whether the charged hadron candidate will be selected or not depends on its rec. tk properties. 
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         // *** Whether the charged hadron candidate will be selected or not depends on its rec. tk properties. 
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         // *** Whether the neutral hadron candidate will be selected or not depends on its rec. HCAL cluster properties. 
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         // *** Whether the gamma candidate will be selected or not depends on its rec. ECAL cluster properties. 
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     // get the initial Track FreeTrajectoryState - at outermost point position if possible, else at innermost point position:
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     // FIX THIS !!!
00204     // need to convert from perigee to global or helix (curvilinear) frame
00205     // for now just an arbitrary matrix.
00206     CLHEP::HepSymMatrix covM_HepSM(6,1); covM_HepSM*=1e-6; // initialize to sigma=1e-3
00207     CartesianTrajectoryError cov_CTE(covM_HepSM);
00208     FreeTrajectoryState Track_initialFTS(theTrack_initialGTPs,cov_CTE);
00209     // ***
00210   
00211     // propagate to ECAL surface: 
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 }