CMSSW Reference Manual

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