Reference Manual

Go to the documentation of this file.
00001 #include "RecoTauTag/TauTagTools/interface/TauElementsOperators.h"
00002 
00003 using namespace reco;
00004 using std::string;
00005 
00006 TauElementsOperators::TauElementsOperators(BaseTau& theBaseTau) : BaseTau_(theBaseTau),AreaMetric_recoElements_maxabsEta_(2.5){
00007   IsolTracks_=theBaseTau.isolationTracks();
00008 }
00009   
00010 double TauElementsOperators::computeConeSize(const TFormula& ConeSizeTFormula,double ConeSizeMin,double ConeSizeMax){
00011   double x=BaseTau_.energy();
00012   double y=BaseTau_.et();
00013   double ConeSize=ConeSizeTFormula.Eval(x,y);
00014   if (ConeSize<ConeSizeMin)ConeSize=ConeSizeMin;
00015   if (ConeSize>ConeSizeMax)ConeSize=ConeSizeMax;
00016   return ConeSize;
00017 }
00018 
00019 double TauElementsOperators::computeConeSize(const TFormula& ConeSizeTFormula,double ConeSizeMin,double ConeSizeMax, double transverseEnergy, double energy, double jetOpeningAngle){
00020   double ConeSize=ConeSizeTFormula.Eval(energy, transverseEnergy, jetOpeningAngle);
00021   if (ConeSize<ConeSizeMin)ConeSize=ConeSizeMin;
00022   if (ConeSize>ConeSizeMax)ConeSize=ConeSizeMax;
00023   return ConeSize;
00024 }
00025 
00026 /*
00027  * DEPRECATED????
00028 TFormula  TauElementsOperators::computeConeSizeTFormula(const string& ConeSizeFormula,const char* errorMessage){
00029   //--- check functional form 
00030   //    given as configuration parameter for matching and signal cone sizes;
00031   //
00032   //    The size of a cone may depend on the energy "E" and/or transverse energy "ET" of the tau-jet candidate.
00033   //    Any functional form that is supported by ROOT's TFormula class can be used (e.g. "3.0/E", "0.25/sqrt(ET)")
00034   //
00035   //    replace "E"  by TFormula variable "x"
00036   //            "ET"                      "y"
00037   string ConeSizeFormulaStr = ConeSizeFormula;
00038   replaceSubStr(ConeSizeFormulaStr,"ET","y");
00039   replaceSubStr(ConeSizeFormulaStr,"E","x");
00040   ConeSizeTFormula.SetName("ConeSize");
00041   ConeSizeTFormula.SetTitle(ConeSizeFormulaStr.data()); // the function definition is actually stored in the "Title" data-member of the TFormula object
00042   int errorFlag = ConeSizeTFormula.Compile();
00043   if (errorFlag!= 0) {
00044     throw cms::Exception("") << "\n unsupported functional Form for " << errorMessage << " " << ConeSizeFormula << endl
00045                              << "Please check that the Definition in \"" << ConeSizeTFormula.GetName() << "\" only contains the variables \"E\" or \"ET\""
00046                              << " and Functions that are supported by ROOT's TFormular Class." << endl;
00047   }else return ConeSizeTFormula;
00048 }
00049 */
00050 
00051 
00052 
00053 void TauElementsOperators::replaceSubStr(string& s,const string& oldSubStr,const string& newSubStr){
00054   //--- protect replacement algorithm
00055   //    from case that oldSubStr and newSubStr are equal
00056   //    (nothing to be done anyway)
00057   if ( oldSubStr == newSubStr ) return;
00058   
00059   //--- protect replacement algorithm
00060   //    from case that oldSubStr contains no characters
00061   //    (i.e. matches everything)
00062   if ( oldSubStr.empty() ) return;
00063   
00064   const string::size_type lengthOldSubStr = oldSubStr.size();
00065   const string::size_type lengthNewSubStr = newSubStr.size();
00066   
00067   string::size_type positionPreviousMatch = 0;
00068   string::size_type positionNextMatch = 0;
00069   
00070   //--- consecutively replace all occurences of oldSubStr by newSubStr;
00071   //    keep iterating until no occurence of oldSubStr left
00072   while ( (positionNextMatch = s.find(oldSubStr, positionPreviousMatch)) != string::npos ) {
00073     s.replace(positionNextMatch, lengthOldSubStr, newSubStr);
00074     positionPreviousMatch = positionNextMatch + lengthNewSubStr;
00075   } 
00076 }
00077 
00078 const TrackRefVector TauElementsOperators::tracksInCone(const math::XYZVector& coneAxis,const string coneMetric,const double coneSize,const double ptTrackMin) const{
00079   TrackRefVector theFilteredTracks;
00080   for (TrackRefVector::const_iterator iTrack=Tracks_.begin();iTrack!=Tracks_.end();++iTrack) {
00081     if ((**iTrack).pt()>ptTrackMin)theFilteredTracks.push_back(*iTrack);
00082   }  
00083   TrackRefVector theFilteredTracksInCone;
00084   if (coneMetric=="DR"){
00085     theFilteredTracksInCone=TracksinCone_DRmetric_(coneAxis,metricDR_,coneSize,theFilteredTracks);
00086   }else if(coneMetric=="angle"){
00087     theFilteredTracksInCone=TracksinCone_Anglemetric_(coneAxis,metricAngle_,coneSize,theFilteredTracks);
00088   }else if(coneMetric=="area"){
00089     int errorFlag = 0;
00090     FixedAreaIsolationCone fixedAreaCone;
00091     fixedAreaCone.setAcceptanceLimit(AreaMetric_recoElements_maxabsEta_);
00092     double coneAngle=fixedAreaCone(coneAxis.theta(),coneAxis.phi(),0,coneSize,errorFlag);
00093     if (errorFlag!=0) return TrackRefVector();
00094     theFilteredTracksInCone=TracksinCone_Anglemetric_(coneAxis,metricAngle_,coneAngle,theFilteredTracks);
00095   }else return TrackRefVector(); 
00096   return theFilteredTracksInCone;
00097 }
00098 const TrackRefVector TauElementsOperators::tracksInCone(const math::XYZVector& coneAxis,const string coneMetric,const double coneSize,const double ptTrackMin,const double tracktorefpoint_maxDZ,const double refpoint_Z, const Vertex &myPV) const{
00099   TrackRefVector theFilteredTracks;
00100   for (TrackRefVector::const_iterator iTrack=Tracks_.begin();iTrack!=Tracks_.end();++iTrack) {
00101     if ((**iTrack).pt()>ptTrackMin && fabs((**iTrack).dz(myPV.position())-refpoint_Z)<=tracktorefpoint_maxDZ)theFilteredTracks.push_back(*iTrack);
00102   }  
00103   TrackRefVector theFilteredTracksInCone;
00104   if (coneMetric=="DR"){
00105     theFilteredTracksInCone=TracksinCone_DRmetric_(coneAxis,metricDR_,coneSize,theFilteredTracks);
00106   }else if(coneMetric=="angle"){
00107     theFilteredTracksInCone=TracksinCone_Anglemetric_(coneAxis,metricAngle_,coneSize,theFilteredTracks);
00108   }else if(coneMetric=="area"){
00109     int errorFlag = 0;
00110     FixedAreaIsolationCone fixedAreaCone;
00111     fixedAreaCone.setAcceptanceLimit(AreaMetric_recoElements_maxabsEta_);
00112     double coneAngle=fixedAreaCone(coneAxis.theta(),coneAxis.phi(),0,coneSize,errorFlag);
00113     if (errorFlag!=0) return TrackRefVector();
00114     theFilteredTracksInCone=TracksinCone_Anglemetric_(coneAxis,metricAngle_,coneAngle,theFilteredTracks);
00115   }else return TrackRefVector(); 
00116   return theFilteredTracksInCone;
00117 }
00118 const TrackRefVector TauElementsOperators::tracksInAnnulus(const math::XYZVector& myVector,const string innercone_metric,const double innercone_size,const string outercone_metric,const double outercone_size,const double minPt)const{     
00119   TrackRefVector theFilteredTracks;
00120   for (TrackRefVector::const_iterator iTrack=Tracks_.begin();iTrack!=Tracks_.end();++iTrack) {
00121     if ((**iTrack).pt()>minPt)theFilteredTracks.push_back(*iTrack);
00122   }  
00123   TrackRefVector theFilteredTracksInAnnulus;
00124   if (outercone_metric=="DR"){
00125     if (innercone_metric=="DR"){
00126       theFilteredTracksInAnnulus=TracksinAnnulus_innerDRouterDRmetrics_(myVector,metricDR_,innercone_size,metricDR_,outercone_size,theFilteredTracks);
00127     }else if(innercone_metric=="angle"){
00128       theFilteredTracksInAnnulus=TracksinAnnulus_innerAngleouterDRmetrics_(myVector,metricAngle_,innercone_size,metricDR_,outercone_size,theFilteredTracks);
00129     }else if(innercone_metric=="area"){
00130       int errorFlag=0;
00131       FixedAreaIsolationCone theFixedAreaSignalCone;
00132       theFixedAreaSignalCone.setAcceptanceLimit(AreaMetric_recoElements_maxabsEta_);
00133       double innercone_angle=theFixedAreaSignalCone(myVector.theta(),myVector.phi(),0,innercone_size,errorFlag);
00134       if (errorFlag!=0)return TrackRefVector();
00135       theFilteredTracksInAnnulus=TracksinAnnulus_innerAngleouterDRmetrics_(myVector,metricAngle_,innercone_angle,metricDR_,outercone_size,theFilteredTracks);
00136     }else return TrackRefVector();
00137   }else if(outercone_metric=="angle"){
00138     if (innercone_metric=="DR"){
00139       theFilteredTracksInAnnulus=TracksinAnnulus_innerDRouterAnglemetrics_(myVector,metricDR_,innercone_size,metricAngle_,outercone_size,theFilteredTracks);
00140     }else if(innercone_metric=="angle"){
00141       theFilteredTracksInAnnulus=TracksinAnnulus_innerAngleouterAnglemetrics_(myVector,metricAngle_,innercone_size,metricAngle_,outercone_size,theFilteredTracks);
00142     }else if(innercone_metric=="area"){
00143       int errorFlag=0;
00144       FixedAreaIsolationCone theFixedAreaSignalCone;
00145       theFixedAreaSignalCone.setAcceptanceLimit(AreaMetric_recoElements_maxabsEta_);
00146       double innercone_angle=theFixedAreaSignalCone(myVector.theta(),myVector.phi(),0,innercone_size,errorFlag);
00147       if (errorFlag!=0)return theFilteredTracksInAnnulus;
00148       theFilteredTracksInAnnulus=TracksinAnnulus_innerAngleouterAnglemetrics_(myVector,metricAngle_,innercone_angle,metricAngle_,outercone_size,theFilteredTracks);
00149     }else return TrackRefVector();
00150   }else if(outercone_metric=="area"){
00151     int errorFlag=0;
00152     FixedAreaIsolationCone theFixedAreaSignalCone;
00153     theFixedAreaSignalCone.setAcceptanceLimit(AreaMetric_recoElements_maxabsEta_);
00154     if (innercone_metric=="DR"){
00155       // not implemented yet
00156     }else if(innercone_metric=="angle"){
00157       double outercone_angle=theFixedAreaSignalCone(myVector.theta(),myVector.phi(),innercone_size,outercone_size,errorFlag);    
00158       if (errorFlag!=0)return theFilteredTracksInAnnulus;
00159       theFilteredTracksInAnnulus=TracksinAnnulus_innerAngleouterAnglemetrics_(myVector,metricAngle_,innercone_size,metricAngle_,outercone_angle,theFilteredTracks);
00160     }else if(innercone_metric=="area"){
00161       double innercone_angle=theFixedAreaSignalCone(myVector.theta(),myVector.phi(),0,innercone_size,errorFlag);    
00162       if (errorFlag!=0)return theFilteredTracksInAnnulus;
00163       double outercone_angle=theFixedAreaSignalCone(myVector.theta(),myVector.phi(),innercone_angle,outercone_size,errorFlag);
00164       if (errorFlag!=0)return theFilteredTracksInAnnulus;
00165       theFilteredTracksInAnnulus=TracksinAnnulus_innerAngleouterAnglemetrics_(myVector,metricAngle_,innercone_angle,metricAngle_,outercone_angle,theFilteredTracks);
00166     }else return TrackRefVector();
00167   }
00168   return theFilteredTracksInAnnulus;
00169 }
00170 const TrackRefVector TauElementsOperators::tracksInAnnulus(const math::XYZVector& myVector,const string innercone_metric,const double innercone_size,const string outercone_metric,const double outercone_size,const double minPt,const double tracktorefpoint_maxDZ,const double refpoint_Z, const Vertex &myPV)const{     
00171   TrackRefVector theFilteredTracks;
00172   for (TrackRefVector::const_iterator iTrack=Tracks_.begin();iTrack!=Tracks_.end();++iTrack) {
00173     if ((**iTrack).pt()>minPt && fabs((**iTrack).dz(myPV.position())-refpoint_Z)<=tracktorefpoint_maxDZ)theFilteredTracks.push_back(*iTrack);
00174   }  
00175   TrackRefVector theFilteredTracksInAnnulus;
00176   if (outercone_metric=="DR"){
00177     if (innercone_metric=="DR"){
00178       theFilteredTracksInAnnulus=TracksinAnnulus_innerDRouterDRmetrics_(myVector,metricDR_,innercone_size,metricDR_,outercone_size,theFilteredTracks);
00179     }else if(innercone_metric=="angle"){
00180       theFilteredTracksInAnnulus=TracksinAnnulus_innerAngleouterDRmetrics_(myVector,metricAngle_,innercone_size,metricDR_,outercone_size,theFilteredTracks);
00181     }else if(innercone_metric=="area"){
00182       int errorFlag=0;
00183       FixedAreaIsolationCone theFixedAreaSignalCone;
00184       theFixedAreaSignalCone.setAcceptanceLimit(AreaMetric_recoElements_maxabsEta_);
00185       double innercone_angle=theFixedAreaSignalCone(myVector.theta(),myVector.phi(),0,innercone_size,errorFlag);
00186       if (errorFlag!=0)return TrackRefVector();
00187       theFilteredTracksInAnnulus=TracksinAnnulus_innerAngleouterDRmetrics_(myVector,metricAngle_,innercone_angle,metricDR_,outercone_size,theFilteredTracks);
00188     }else return TrackRefVector();
00189   }else if(outercone_metric=="angle"){
00190     if (innercone_metric=="DR"){
00191       theFilteredTracksInAnnulus=TracksinAnnulus_innerDRouterAnglemetrics_(myVector,metricDR_,innercone_size,metricAngle_,outercone_size,theFilteredTracks);
00192     }else if(innercone_metric=="angle"){
00193       theFilteredTracksInAnnulus=TracksinAnnulus_innerAngleouterAnglemetrics_(myVector,metricAngle_,innercone_size,metricAngle_,outercone_size,theFilteredTracks);
00194     }else if(innercone_metric=="area"){
00195       int errorFlag=0;
00196       FixedAreaIsolationCone theFixedAreaSignalCone;
00197       theFixedAreaSignalCone.setAcceptanceLimit(AreaMetric_recoElements_maxabsEta_);
00198       double innercone_angle=theFixedAreaSignalCone(myVector.theta(),myVector.phi(),0,innercone_size,errorFlag);
00199       if (errorFlag!=0)return theFilteredTracksInAnnulus;
00200       theFilteredTracksInAnnulus=TracksinAnnulus_innerAngleouterAnglemetrics_(myVector,metricAngle_,innercone_angle,metricAngle_,outercone_size,theFilteredTracks);
00201     }else return TrackRefVector();
00202   }else if(outercone_metric=="area"){
00203     int errorFlag=0;
00204     FixedAreaIsolationCone theFixedAreaSignalCone;
00205     theFixedAreaSignalCone.setAcceptanceLimit(AreaMetric_recoElements_maxabsEta_);
00206     if (innercone_metric=="DR"){
00207       // not implemented yet
00208     }else if(innercone_metric=="angle"){
00209       double outercone_angle=theFixedAreaSignalCone(myVector.theta(),myVector.phi(),innercone_size,outercone_size,errorFlag);    
00210       if (errorFlag!=0)return theFilteredTracksInAnnulus;
00211       theFilteredTracksInAnnulus=TracksinAnnulus_innerAngleouterAnglemetrics_(myVector,metricAngle_,innercone_size,metricAngle_,outercone_angle,theFilteredTracks);
00212     }else if(innercone_metric=="area"){
00213       double innercone_angle=theFixedAreaSignalCone(myVector.theta(),myVector.phi(),0,innercone_size,errorFlag);    
00214       if (errorFlag!=0)return theFilteredTracksInAnnulus;
00215       double outercone_angle=theFixedAreaSignalCone(myVector.theta(),myVector.phi(),innercone_angle,outercone_size,errorFlag);
00216       if (errorFlag!=0)return theFilteredTracksInAnnulus;
00217       theFilteredTracksInAnnulus=TracksinAnnulus_innerAngleouterAnglemetrics_(myVector,metricAngle_,innercone_angle,metricAngle_,outercone_angle,theFilteredTracks);
00218     }else return TrackRefVector();
00219   }
00220   return theFilteredTracksInAnnulus;
00221 }
00222 
00223 const TrackRef TauElementsOperators::leadTk(string matchingConeMetric,double matchingConeSize,double ptTrackMin)const{
00224   return leadTk(BaseTau_.momentum(),matchingConeMetric,matchingConeSize,ptTrackMin);
00225 }
00226 
00227 const TrackRef TauElementsOperators::leadTk(const math::XYZVector& jetAxis,string matchingConeMetric,double matchingConeSize,double ptTrackMin)const{
00228   const TrackRefVector matchingConeTracks=tracksInCone(jetAxis,matchingConeMetric,matchingConeSize,ptTrackMin);
00229   if ((int)matchingConeTracks.size()==0) return TrackRef();
00230   TrackRef leadingTrack;
00231   double leadingTrackPt=0.;
00232   for (TrackRefVector::const_iterator track=matchingConeTracks.begin();track!=matchingConeTracks.end();++track) {
00233     if ((*track)->pt()>ptTrackMin && (*track)->pt()>leadingTrackPt){
00234       leadingTrack=(*track);
00235       leadingTrackPt=leadingTrack->pt();
00236     }
00237   }  
00238   return leadingTrack;
00239 }
00240 // ***
00241 double TauElementsOperators::discriminatorByIsolTracksN(unsigned int isolationAnnulus_Tracksmaxn)const{
00242   if ((unsigned int)IsolTracks_.size()>isolationAnnulus_Tracksmaxn)return 0.;
00243   else return 1.;
00244 }
00245 double TauElementsOperators::discriminatorByIsolTracksN(const math::XYZVector& jetAxis, 
00246                                                         string matchingConeMetric,double matchingConeSize,double ptLeadingTrackMin,double ptOtherTracksMin,
00247                                                         string signalConeMetric,double signalConeSize,string isolationConeMetric,double isolationConeSize, 
00248                                                         unsigned int isolationAnnulus_Tracksmaxn)const{
00249   const TrackRef leadingTrack=leadTk(jetAxis,matchingConeMetric,matchingConeSize,ptLeadingTrackMin);
00250   if(!leadingTrack)return 0.; 
00251   math::XYZVector coneAxis=leadingTrack->momentum();
00252   TrackRefVector isolationAnnulusTracks=tracksInAnnulus(coneAxis,signalConeMetric,signalConeSize,isolationConeMetric,isolationConeSize,ptOtherTracksMin);
00253   if ((unsigned int)isolationAnnulusTracks.size()>isolationAnnulus_Tracksmaxn)return 0.;
00254   else return 1.;
00255 }
00256 double TauElementsOperators::discriminatorByIsolTracksN(string matchingConeMetric,double matchingConeSize,double ptLeadingTrackMin,double ptOtherTracksMin, 
00257                                                         string signalConeMetric,double signalConeSize,string isolationConeMetric,double isolationConeSize, 
00258                                                         unsigned int isolationAnnulus_Tracksmaxn)const{
00259   return discriminatorByIsolTracksN(BaseTau_.momentum(),matchingConeMetric,matchingConeSize,ptLeadingTrackMin,ptOtherTracksMin,signalConeMetric,signalConeSize,isolationConeMetric,isolationConeSize,isolationAnnulus_Tracksmaxn);
00260 }
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/RecoTauTag/TauTagTools/src/TauElementsOperators.cc
