CMS 3D CMS Logo

TauTagTools.cc
Go to the documentation of this file.
2 
3 using namespace reco;
4 using std::string;
5 
6 namespace TauTagTools{
7 
9  {
10  DeltaR<math::XYZVector> myMetricDR_;
11  return myMetricDR_(vec1, vec2);
12  }
14  {
15  Angle<math::XYZVector> myMetricAngle_;
16  return myMetricAngle_(vec1, vec2);
17  }
18 TFormula computeConeSizeTFormula(const string& ConeSizeFormula,const char* errorMessage){
19  //--- check functional form
20  // given as configuration parameter for matching and signal cone sizes;
21  //
22  // The size of a cone may depend on the energy "E" and/or transverse energy "ET" of the tau-jet candidate.
23  // Alternatively one can additionally use "JetOpeningDR", which specifies the opening angle of the seed jet.
24  // Any functional form that is supported by ROOT's TFormula class can be used (e.g. "3.0/E", "0.25/sqrt(ET)")
25  //
26  // replace "E" by TFormula variable "x"
27  // "ET" "y"
28  TFormula ConeSizeTFormula;
29  string ConeSizeFormulaStr = ConeSizeFormula;
30  replaceSubStr(ConeSizeFormulaStr,"JetOpeningDR", "z");
31  replaceSubStr(ConeSizeFormulaStr,"ET","y");
32  replaceSubStr(ConeSizeFormulaStr,"E","x");
33  ConeSizeTFormula.SetName("ConeSize");
34  ConeSizeTFormula.SetTitle(ConeSizeFormulaStr.data()); // the function definition is actually stored in the "Title" data-member of the TFormula object
35  int errorFlag = ConeSizeTFormula.Compile();
36  if (errorFlag!= 0) {
37  throw cms::Exception("") << "\n unsupported functional Form for " << errorMessage << " " << ConeSizeFormula << std::endl
38  << "Please check that the Definition in \"" << ConeSizeTFormula.GetName() << "\" only contains the variables \"E\" or \"ET\""
39  << " and Functions that are supported by ROOT's TFormular Class." << std::endl;
40  }else return ConeSizeTFormula;
41 }
42 void replaceSubStr(string& s,const string& oldSubStr,const string& newSubStr){
43  //--- protect replacement algorithm
44  // from case that oldSubStr and newSubStr are equal
45  // (nothing to be done anyway)
46  if ( oldSubStr == newSubStr ) return;
47 
48  //--- protect replacement algorithm
49  // from case that oldSubStr contains no characters
50  // (i.e. matches everything)
51  if ( oldSubStr.empty() ) return;
52 
53  const string::size_type lengthOldSubStr = oldSubStr.size();
54  const string::size_type lengthNewSubStr = newSubStr.size();
55 
56  string::size_type positionPreviousMatch = 0;
57  string::size_type positionNextMatch = 0;
58 
59  //--- consecutively replace all occurences of oldSubStr by newSubStr;
60  // keep iterating until no occurence of oldSubStr left
61  while ( (positionNextMatch = s.find(oldSubStr, positionPreviousMatch)) != string::npos ) {
62  s.replace(positionNextMatch, lengthOldSubStr, newSubStr);
63  positionPreviousMatch = positionNextMatch + lengthNewSubStr;
64  }
65 }
66 
67 
70  for(TrackRefVector::const_iterator iTk=theInitialTracks.begin();iTk!=theInitialTracks.end();iTk++){
71  if ( (**iTk).numberOfValidHits() >= tkminTrackerHitsn )
72  filteredTracks.push_back(*iTk);
73  }
74  return filteredTracks;
75  }
76 
79  for(TrackRefVector::const_iterator iTk=theInitialTracks.begin();iTk!=theInitialTracks.end();iTk++){
80  if ((**iTk).pt()>=tkminPt &&
81  (**iTk).normalizedChi2()<=tkmaxChi2 &&
82  fabs((**iTk).dxy(pv.position()))<=tkmaxipt &&
83  (**iTk).numberOfValidHits()>=tkminTrackerHitsn &&
84  (**iTk).hitPattern().numberOfValidPixelHits()>=tkminPixelHitsn)
85  filteredTracks.push_back(*iTk);
86  }
87  return filteredTracks;
88  }
89  TrackRefVector filteredTracks(const TrackRefVector& theInitialTracks,double tkminPt,int tkminPixelHitsn,int tkminTrackerHitsn,double tkmaxipt,double tkmaxChi2,double tktorefpointmaxDZ,Vertex pv, double refpoint_Z){
91  for(TrackRefVector::const_iterator iTk=theInitialTracks.begin();iTk!=theInitialTracks.end();iTk++){
92  if(pv.isFake()) tktorefpointmaxDZ=30.;
93  if ((**iTk).pt()>=tkminPt &&
94  (**iTk).normalizedChi2()<=tkmaxChi2 &&
95  fabs((**iTk).dxy(pv.position()))<=tkmaxipt &&
96  (**iTk).numberOfValidHits()>=tkminTrackerHitsn &&
97  (**iTk).hitPattern().numberOfValidPixelHits()>=tkminPixelHitsn &&
98  fabs((**iTk).dz(pv.position()))<=tktorefpointmaxDZ)
99  filteredTracks.push_back(*iTk);
100  }
101  return filteredTracks;
102  }
103 
104  std::vector<reco::CandidatePtr> filteredPFChargedHadrCandsByNumTrkHits(const std::vector<reco::CandidatePtr>& theInitialPFCands, int ChargedHadrCand_tkminTrackerHitsn){
105  std::vector<reco::CandidatePtr> filteredPFChargedHadrCands;
106  for(std::vector<reco::CandidatePtr>::const_iterator iPFCand=theInitialPFCands.begin();iPFCand!=theInitialPFCands.end();iPFCand++){
107  if (std::abs((*iPFCand)->pdgId()) == 211 || std::abs((*iPFCand)->pdgId()) == 13 || std::abs((*iPFCand)->pdgId()) == 11){
108  // *** Whether the charged hadron candidate will be selected or not depends on its rec. tk properties.
109  const reco::Track* PFChargedHadrCand_rectk = (*iPFCand)->bestTrack();
110  if (PFChargedHadrCand_rectk != nullptr) {
111  if ( (*PFChargedHadrCand_rectk).numberOfValidHits()>=ChargedHadrCand_tkminTrackerHitsn )
112  filteredPFChargedHadrCands.push_back(*iPFCand);
113  }
114  }
115  }
117  }
118 
119  std::vector<reco::CandidatePtr> filteredPFChargedHadrCands(const std::vector<reco::CandidatePtr>& theInitialPFCands,double ChargedHadrCand_tkminPt,int ChargedHadrCand_tkminPixelHitsn,int ChargedHadrCand_tkminTrackerHitsn,double ChargedHadrCand_tkmaxipt,double ChargedHadrCand_tkmaxChi2, Vertex pv){
120  std::vector<reco::CandidatePtr> filteredPFChargedHadrCands;
121  for(std::vector<reco::CandidatePtr>::const_iterator iPFCand=theInitialPFCands.begin();iPFCand!=theInitialPFCands.end();iPFCand++){
122  if (std::abs((*iPFCand)->pdgId()) == 211 || std::abs((*iPFCand)->pdgId()) == 13 || std::abs((*iPFCand)->pdgId()) == 11){
123  // *** Whether the charged hadron candidate will be selected or not depends on its rec. tk properties.
124  const reco::Track* PFChargedHadrCand_rectk = (*iPFCand)->bestTrack();
125  if (PFChargedHadrCand_rectk != nullptr) {
126  if ((*PFChargedHadrCand_rectk).pt()>=ChargedHadrCand_tkminPt &&
127  (*PFChargedHadrCand_rectk).normalizedChi2()<=ChargedHadrCand_tkmaxChi2 &&
128  fabs((*PFChargedHadrCand_rectk).dxy(pv.position()))<=ChargedHadrCand_tkmaxipt &&
129  (*PFChargedHadrCand_rectk).numberOfValidHits()>=ChargedHadrCand_tkminTrackerHitsn &&
130  (*PFChargedHadrCand_rectk).hitPattern().numberOfValidPixelHits()>=ChargedHadrCand_tkminPixelHitsn)
131  filteredPFChargedHadrCands.push_back(*iPFCand);
132  }
133  }
134  }
136  }
137  std::vector<reco::CandidatePtr> filteredPFChargedHadrCands(const std::vector<reco::CandidatePtr>& theInitialPFCands,double ChargedHadrCand_tkminPt,int ChargedHadrCand_tkminPixelHitsn,int ChargedHadrCand_tkminTrackerHitsn,double ChargedHadrCand_tkmaxipt,double ChargedHadrCand_tkmaxChi2,double ChargedHadrCand_tktorefpointmaxDZ,Vertex pv, double refpoint_Z){
138  if(pv.isFake()) ChargedHadrCand_tktorefpointmaxDZ = 30.;
139  std::vector<reco::CandidatePtr> filteredPFChargedHadrCands;
140  for(std::vector<reco::CandidatePtr>::const_iterator iPFCand=theInitialPFCands.begin();iPFCand!=theInitialPFCands.end();iPFCand++){
141  if (std::abs((*iPFCand)->pdgId()) == 211 || std::abs((*iPFCand)->pdgId()) == 13 || std::abs((*iPFCand)->pdgId()) == 11){
142  // *** Whether the charged hadron candidate will be selected or not depends on its rec. tk properties.
143  const reco::Track* PFChargedHadrCand_rectk = (*iPFCand)->bestTrack();
144  if (PFChargedHadrCand_rectk != nullptr) {
145  if ((*PFChargedHadrCand_rectk).pt()>=ChargedHadrCand_tkminPt &&
146  (*PFChargedHadrCand_rectk).normalizedChi2()<=ChargedHadrCand_tkmaxChi2 &&
147  fabs((*PFChargedHadrCand_rectk).dxy(pv.position()))<=ChargedHadrCand_tkmaxipt &&
148  (*PFChargedHadrCand_rectk).numberOfValidHits()>=ChargedHadrCand_tkminTrackerHitsn &&
149  (*PFChargedHadrCand_rectk).hitPattern().numberOfValidPixelHits()>=ChargedHadrCand_tkminPixelHitsn &&
150  fabs((*PFChargedHadrCand_rectk).dz(pv.position()))<=ChargedHadrCand_tktorefpointmaxDZ)
151  filteredPFChargedHadrCands.push_back(*iPFCand);
152  }
153  }
154  }
156  }
157 
158  std::vector<reco::CandidatePtr> filteredPFNeutrHadrCands(const std::vector<reco::CandidatePtr>& theInitialPFCands,double NeutrHadrCand_HcalclusMinEt){
159  std::vector<reco::CandidatePtr> filteredPFNeutrHadrCands;
160  for(std::vector<reco::CandidatePtr>::const_iterator iPFCand=theInitialPFCands.begin();iPFCand!=theInitialPFCands.end();iPFCand++){
161  if (std::abs((*iPFCand)->pdgId()) == 130){
162  // *** Whether the neutral hadron candidate will be selected or not depends on its rec. HCAL cluster properties.
163  if ((**iPFCand).et()>=NeutrHadrCand_HcalclusMinEt){
164  filteredPFNeutrHadrCands.push_back(*iPFCand);
165  }
166  }
167  }
169  }
170 
171  std::vector<reco::CandidatePtr> filteredPFGammaCands(const std::vector<reco::CandidatePtr>& theInitialPFCands,double GammaCand_EcalclusMinEt){
172  std::vector<reco::CandidatePtr> filteredPFGammaCands;
173  for(std::vector<reco::CandidatePtr>::const_iterator iPFCand=theInitialPFCands.begin();iPFCand!=theInitialPFCands.end();iPFCand++){
174  if (std::abs((*iPFCand)->pdgId()) == 22){
175  // *** Whether the gamma candidate will be selected or not depends on its rec. ECAL cluster properties.
176  if ((**iPFCand).et()>=GammaCand_EcalclusMinEt){
177  filteredPFGammaCands.push_back(*iPFCand);
178  }
179  }
180  }
181  return filteredPFGammaCands;
182  }
183 
185  AnalyticalPropagator thefwdPropagator(theMagField,alongMomentum);
186  math::XYZPoint propTrack_XYZPoint(0.,0.,0.);
187 
188  // get the initial Track FreeTrajectoryState - at outermost point position if possible, else at innermost point position:
189  GlobalVector theTrack_initialGV(0.,0.,0.);
190  GlobalPoint theTrack_initialGP(0.,0.,0.);
191  if(theTrack->outerOk()){
192  GlobalVector theTrack_initialoutermostGV(theTrack->outerMomentum().x(),theTrack->outerMomentum().y(),theTrack->outerMomentum().z());
193  GlobalPoint theTrack_initialoutermostGP(theTrack->outerPosition().x(),theTrack->outerPosition().y(),theTrack->outerPosition().z());
194  theTrack_initialGV=theTrack_initialoutermostGV;
195  theTrack_initialGP=theTrack_initialoutermostGP;
196  } else if(theTrack->innerOk()){
197  GlobalVector theTrack_initialinnermostGV(theTrack->innerMomentum().x(),theTrack->innerMomentum().y(),theTrack->innerMomentum().z());
198  GlobalPoint theTrack_initialinnermostGP(theTrack->innerPosition().x(),theTrack->innerPosition().y(),theTrack->innerPosition().z());
199  theTrack_initialGV=theTrack_initialinnermostGV;
200  theTrack_initialGP=theTrack_initialinnermostGP;
201  } else return (propTrack_XYZPoint);
202  GlobalTrajectoryParameters theTrack_initialGTPs(theTrack_initialGP,theTrack_initialGV,theTrack->charge(),&*theMagField);
203  // FIX THIS !!!
204  // need to convert from perigee to global or helix (curvilinear) frame
205  // for now just an arbitrary matrix.
206  AlgebraicSymMatrix66 covM=AlgebraicMatrixID(); covM *= 1e-6; // initialize to sigma=1e-3
207  CartesianTrajectoryError cov_CTE(covM);
208  FreeTrajectoryState Track_initialFTS(theTrack_initialGTPs,cov_CTE);
209  // ***
210 
211  // propagate to ECAL surface:
212  double ECALcorner_tantheta=ECALBounds::barrel_innerradius()/ECALBounds::barrel_halfLength();
213  TrajectoryStateOnSurface Track_propagatedonECAL_TSOS=thefwdPropagator.propagate((Track_initialFTS),ECALBounds::barrelBound());
214  if(!Track_propagatedonECAL_TSOS.isValid() || fabs(Track_propagatedonECAL_TSOS.globalParameters().position().perp()/Track_propagatedonECAL_TSOS.globalParameters().position().z())<ECALcorner_tantheta) {
215  if(Track_propagatedonECAL_TSOS.isValid() && fabs(Track_propagatedonECAL_TSOS.globalParameters().position().perp()/Track_propagatedonECAL_TSOS.globalParameters().position().z())<ECALcorner_tantheta){
216  if(Track_propagatedonECAL_TSOS.globalParameters().position().eta()>0.){
217  Track_propagatedonECAL_TSOS=thefwdPropagator.propagate((Track_initialFTS),ECALBounds::positiveEndcapDisk());
218  }else{
219  Track_propagatedonECAL_TSOS=thefwdPropagator.propagate((Track_initialFTS),ECALBounds::negativeEndcapDisk());
220  }
221  }
222  if(!Track_propagatedonECAL_TSOS.isValid()){
223  if((Track_initialFTS).position().eta()>0.){
224  Track_propagatedonECAL_TSOS=thefwdPropagator.propagate((Track_initialFTS),ECALBounds::positiveEndcapDisk());
225  }else{
226  Track_propagatedonECAL_TSOS=thefwdPropagator.propagate((Track_initialFTS),ECALBounds::negativeEndcapDisk());
227  }
228  }
229  }
230  if(Track_propagatedonECAL_TSOS.isValid()){
231  math::XYZPoint validpropTrack_XYZPoint(Track_propagatedonECAL_TSOS.globalPosition().x(),
232  Track_propagatedonECAL_TSOS.globalPosition().y(),
233  Track_propagatedonECAL_TSOS.globalPosition().z());
234  propTrack_XYZPoint=validpropTrack_XYZPoint;
235  }
236  return (propTrack_XYZPoint);
237  }
238 
239 
240 
241  void
242  sortRefVectorByPt(std::vector<reco::CandidatePtr>& vec)
243  {
244 
245  std::vector<size_t> indices;
246  indices.reserve(vec.size());
247  for(unsigned int i=0;i<vec.size();++i)
248  indices.push_back(i);
249 
251  std::sort(indices.begin(),indices.end(),sorter);
252 
253 
254  std::vector<reco::CandidatePtr> sorted;
255  sorted.reserve(vec.size());
256 
257  for(unsigned int i=0;i<indices.size();++i)
258  sorted.push_back(vec.at(indices.at(i)));
259 
260  vec = sorted;
261  }
262 
263 
264 }
std::vector< reco::CandidatePtr > filteredPFChargedHadrCands(const std::vector< reco::CandidatePtr > &theInitialPFCands, double ChargedHadrCand_tkminPt, int ChargedHadrCand_tkminPixelHitsn, int ChargedHadrCand_tkminTrackerHitsn, double ChargedHadrCand_tkmaxipt, double ChargedHadrCand_tkmaxChi2, double ChargedHadrCand_tktorefpointmaxDZ, Vertex pv, double refpoint_Z)
Definition: TauTagTools.cc:137
T perp() const
Definition: PV3DBase.h:72
double computeAngle(const math::XYZVector &vec1, const math::XYZVector &vec2)
Definition: TauTagTools.cc:13
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
void replaceSubStr(string &s, const string &oldSubStr, const string &newSubStr)
Definition: TauTagTools.cc:42
ROOT::Math::SMatrix< double, 6, 6, ROOT::Math::MatRepSym< double, 6 > > AlgebraicSymMatrix66
math::XYZPoint propagTrackECALSurfContactPoint(const MagneticField *, reco::TrackRef)
Definition: TauTagTools.cc:184
T y() const
Definition: PV3DBase.h:63
GlobalPoint globalPosition() const
const_iterator end() const
Termination of iteration.
Definition: RefVector.h:253
uint16_t size_type
const Point & position() const
position
Definition: Vertex.h:109
static float barrel_innerradius()
Definition: ECALBounds.h:25
const_iterator begin() const
Initialize an iterator over the RefVector.
Definition: RefVector.h:248
static const BoundCylinder & barrelBound()
Definition: ECALBounds.h:19
static const BoundDisk & positiveEndcapDisk()
Definition: ECALBounds.h:21
std::vector< reco::CandidatePtr > filteredPFNeutrHadrCands(const std::vector< reco::CandidatePtr > &theInitialPFCands, double NeutrHadrCand_HcalclusMinEt)
Definition: TauTagTools.cc:158
static const BoundDisk & negativeEndcapDisk()
Definition: ECALBounds.h:20
std::vector< double > vec1
Definition: HCALResponse.h:15
T z() const
Definition: PV3DBase.h:64
def pv(vc)
Definition: MetAnalyzer.py:7
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
void sortRefVectorByPt(std::vector< reco::CandidatePtr > &vec)
Definition: TauTagTools.cc:242
TFormula computeConeSizeTFormula(const string &ConeSizeFormula, const char *errorMessage)
Definition: TauTagTools.cc:18
double computeDeltaR(const math::XYZVector &vec1, const math::XYZVector &vec2)
Definition: TauTagTools.cc:8
std::vector< reco::CandidatePtr > filteredPFGammaCands(const std::vector< reco::CandidatePtr > &theInitialPFCands, double GammaCand_EcalclusMinEt)
Definition: TauTagTools.cc:171
bool isFake() const
Definition: Vertex.h:72
const GlobalTrajectoryParameters & globalParameters() const
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:30
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
std::vector< reco::CandidatePtr > filteredPFChargedHadrCandsByNumTrkHits(const std::vector< reco::CandidatePtr > &theInitialPFCands, int ChargedHadrCand_tkminTrackerHitsn)
Definition: TauTagTools.cc:104
T eta() const
Definition: PV3DBase.h:76
reco::TrackRefVector filteredTracks(const reco::TrackRefVector &theInitialTracks, double tkminPt, int tkminPixelHitsn, int tkminTrackerHitsn, double tkmaxipt, double tkmaxChi2, reco::Vertex pV)
Definition: TauTagTools.cc:77
fixed size matrix
static float barrel_halfLength()
Definition: ECALBounds.h:27
reco::TrackRefVector filteredTracksByNumTrkHits(const reco::TrackRefVector &theInitialTracks, int tkminTrackerHitsn)
Definition: TauTagTools.cc:68
void push_back(value_type const &ref)
Add a Ref<C, T> to the RefVector.
Definition: RefVector.h:69
std::vector< vec1 > vec2
Definition: HCALResponse.h:16
T x() const
Definition: PV3DBase.h:62