CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TrackAssociatorByChi2.cc
Go to the documentation of this file.
1 #include "SimTracker/TrackAssociation/interface/TrackAssociatorByChi2.h"
2 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
3 
4 #include "DataFormats/Math/interface/deltaPhi.h"
5 #include "DataFormats/GeometrySurface/interface/Line.h"
6 #include "DataFormats/GeometryVector/interface/Pi.h"
7 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
8 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
9 using namespace edm;
10 using namespace reco;
11 using namespace std;
12 
13 double TrackAssociatorByChi2::compareTracksParam ( TrackCollection::const_iterator rt,
14  SimTrackContainer::const_iterator st,
15  const math::XYZTLorentzVectorD vertexPosition,
16  GlobalVector magField,
17  TrackBase::CovarianceMatrix invertedCovariance,
18  const reco::BeamSpot& bs) const{
19 
20  Basic3DVector<double> momAtVtx(st->momentum().x(),st->momentum().y(),st->momentum().z());
21  Basic3DVector<double> vert = (Basic3DVector<double>) vertexPosition;
22 
23  std::pair<bool,reco::TrackBase::ParameterVector> params = parametersAtClosestApproach(vert, momAtVtx, st->charge(), bs);
24  if (params.first){
25  TrackBase::ParameterVector sParameters = params.second;
26  TrackBase::ParameterVector rParameters = rt->parameters();
27 
28  TrackBase::ParameterVector diffParameters = rParameters - sParameters;
29  diffParameters[2] = reco::deltaPhi(diffParameters[2],0.f);
30  double chi2 = ROOT::Math::Dot(diffParameters * invertedCovariance, diffParameters);
31 
32  return chi2;
33  } else {
34  return 10000000000.;
35  }
36 }
37 
38 
39 TrackAssociatorByChi2::RecoToSimPairAssociation
40 TrackAssociatorByChi2::compareTracksParam(const TrackCollection& rtColl,
41  const SimTrackContainer& stColl,
42  const SimVertexContainer& svColl,
43  const reco::BeamSpot& bs) const{
44 
45  RecoToSimPairAssociation outputVec;
46 
47  for (TrackCollection::const_iterator track=rtColl.begin(); track!=rtColl.end(); track++){
48  Chi2SimMap outMap;
49 
50  TrackBase::ParameterVector rParameters = track->parameters();
51 
52  TrackBase::CovarianceMatrix recoTrackCovMatrix = track->covariance();
53  if (onlyDiagonal){
54  for (unsigned int i=0;i<5;i++){
55  for (unsigned int j=0;j<5;j++){
56  if (i!=j) recoTrackCovMatrix(i,j)=0;
57  }
58  }
59  }
60  recoTrackCovMatrix.Invert();
61 
62  for (SimTrackContainer::const_iterator st=stColl.begin(); st!=stColl.end(); st++){
63 
64  Basic3DVector<double> momAtVtx(st->momentum().x(),st->momentum().y(),st->momentum().z());
65  Basic3DVector<double> vert = (Basic3DVector<double>) svColl[st->vertIndex()].position();
66 
67  std::pair<bool,reco::TrackBase::ParameterVector> params = parametersAtClosestApproach(vert, momAtVtx, st->charge(), bs);
68  if (params.first){
69  TrackBase::ParameterVector sParameters = params.second;
70 
71  TrackBase::ParameterVector diffParameters = rParameters - sParameters;
72  diffParameters[2] = reco::deltaPhi(diffParameters[2],0.f);
73  double chi2 = ROOT::Math::Dot(diffParameters * recoTrackCovMatrix, diffParameters);
74  chi2/=5;
75  if (chi2<chi2cut) outMap[chi2]=*st;
76  }
77  }
78  outputVec.push_back(RecoToSimPair(*track,outMap));
79  }
80  return outputVec;
81 }
82 
83 
84 double TrackAssociatorByChi2::associateRecoToSim( TrackCollection::const_iterator rt,
85  TrackingParticleCollection::const_iterator tp,
86  const reco::BeamSpot& bs) const{
87 
88  double chi2;
89 
90  TrackBase::ParameterVector rParameters = rt->parameters();
91 
92  TrackBase::CovarianceMatrix recoTrackCovMatrix = rt->covariance();
93  if (onlyDiagonal) {
94  for (unsigned int i=0;i<5;i++){
95  for (unsigned int j=0;j<5;j++){
96  if (i!=j) recoTrackCovMatrix(i,j)=0;
97  }
98  }
99  }
100  recoTrackCovMatrix.Invert();
101 
102  Basic3DVector<double> momAtVtx(tp->momentum().x(),tp->momentum().y(),tp->momentum().z());
103  Basic3DVector<double> vert(tp->vertex().x(),tp->vertex().y(),tp->vertex().z());
104 
105  std::pair<bool,reco::TrackBase::ParameterVector> params = parametersAtClosestApproach(vert, momAtVtx, tp->charge(), bs);
106  if (params.first){
107  TrackBase::ParameterVector sParameters=params.second;
108 
109  TrackBase::ParameterVector diffParameters = rParameters - sParameters;
110  diffParameters[2] = reco::deltaPhi(diffParameters[2],0.f);
111  chi2 = ROOT::Math::Dot(diffParameters * recoTrackCovMatrix, diffParameters);
112  chi2 /= 5;
113 
114  LogDebug("TrackAssociator") << "====NEW RECO TRACK WITH PT=" << rt->pt() << "====\n"
115  << "qoverp sim: " << sParameters[0] << "\n"
116  << "lambda sim: " << sParameters[1] << "\n"
117  << "phi sim: " << sParameters[2] << "\n"
118  << "dxy sim: " << sParameters[3] << "\n"
119  << "dsz sim: " << sParameters[4] << "\n"
120  << ": " /*<< */ << "\n"
121  << "qoverp rec: " << rt->qoverp()/*rParameters[0]*/ << "\n"
122  << "lambda rec: " << rt->lambda()/*rParameters[1]*/ << "\n"
123  << "phi rec: " << rt->phi()/*rParameters[2]*/ << "\n"
124  << "dxy rec: " << rt->dxy(bs.position())/*rParameters[3]*/ << "\n"
125  << "dsz rec: " << rt->dsz(bs.position())/*rParameters[4]*/ << "\n"
126  << ": " /*<< */ << "\n"
127  << "chi2: " << chi2 << "\n";
128 
129  return chi2;
130  } else {
131  return 10000000000.;
132  }
133 }
134 
135 pair<bool,TrackBase::ParameterVector>
136 TrackAssociatorByChi2::parametersAtClosestApproach(Basic3DVector<double> vertex,
137  Basic3DVector<double> momAtVtx,
138  float charge,
139  const BeamSpot& bs) const{
140 
141  TrackBase::ParameterVector sParameters;
142  try {
143  FreeTrajectoryState ftsAtProduction(GlobalPoint(vertex.x(),vertex.y(),vertex.z()),
144  GlobalVector(momAtVtx.x(),momAtVtx.y(),momAtVtx.z()),
145  TrackCharge(charge),
146  theMF.product());
147  TSCBLBuilderNoMaterial tscblBuilder;
148  TrajectoryStateClosestToBeamLine tsAtClosestApproach = tscblBuilder(ftsAtProduction,bs);//as in TrackProducerAlgorithm
149 
150  GlobalPoint v = tsAtClosestApproach.trackStateAtPCA().position();
151  GlobalVector p = tsAtClosestApproach.trackStateAtPCA().momentum();
152  sParameters[0] = tsAtClosestApproach.trackStateAtPCA().charge()/p.mag();
153  sParameters[1] = Geom::halfPi() - p.theta();
154  sParameters[2] = p.phi();
155  sParameters[3] = (-v.x()*sin(p.phi())+v.y()*cos(p.phi()));
156  sParameters[4] = v.z()*p.perp()/p.mag() - (v.x()*p.x()+v.y()*p.y())/p.perp() * p.z()/p.mag();
157 
158  return pair<bool,TrackBase::ParameterVector>(true,sParameters);
159  } catch ( ... ) {
160  return pair<bool,TrackBase::ParameterVector>(false,sParameters);
161  }
162 }
163 
164 RecoToSimCollection TrackAssociatorByChi2::associateRecoToSim(const edm::RefToBaseVector<reco::Track>& tC,
165  const edm::RefVector<TrackingParticleCollection>& tPCH,
166  const edm::Event * e,
167  const edm::EventSetup *setup ) const{
168  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
169  e->getByLabel(bsSrc,recoBeamSpotHandle);
170  reco::BeamSpot bs = *recoBeamSpotHandle;
171 
172  RecoToSimCollection outputCollection;
173  double chi2;
174 
175  TrackingParticleCollection tPC;
176  if (tPCH.size()!=0) tPC = *tPCH.product();
177 
178  int tindex=0;
179  for (RefToBaseVector<reco::Track>::const_iterator rt=tC.begin(); rt!=tC.end(); rt++, tindex++){
180 
181  LogDebug("TrackAssociator") << "=========LOOKING FOR ASSOCIATION===========" << "\n"
182  << "rec::Track #"<<tindex<<" with pt=" << (*rt)->pt() << "\n"
183  << "===========================================" << "\n";
184 
185  TrackBase::ParameterVector rParameters = (*rt)->parameters();
186 
187  TrackBase::CovarianceMatrix recoTrackCovMatrix = (*rt)->covariance();
188  if (onlyDiagonal){
189  for (unsigned int i=0;i<5;i++){
190  for (unsigned int j=0;j<5;j++){
191  if (i!=j) recoTrackCovMatrix(i,j)=0;
192  }
193  }
194  }
195 
196  recoTrackCovMatrix.Invert();
197 
198  int tpindex =0;
199  for (TrackingParticleCollection::const_iterator tp=tPC.begin(); tp!=tPC.end(); tp++, ++tpindex){
200 
201  //skip tps with a very small pt
202  //if (sqrt(tp->momentum().perp2())<0.5) continue;
203  Basic3DVector<double> momAtVtx(tp->momentum().x(),tp->momentum().y(),tp->momentum().z());
204  Basic3DVector<double> vert=(Basic3DVector<double>) tp->vertex();
205 
206  std::pair<bool,reco::TrackBase::ParameterVector> params = parametersAtClosestApproach(vert, momAtVtx, tp->charge(), bs);
207  if (params.first){
208  TrackBase::ParameterVector sParameters=params.second;
209 
210  TrackBase::ParameterVector diffParameters = rParameters - sParameters;
211  diffParameters[2] = reco::deltaPhi(diffParameters[2],0.f);
212 
213  chi2 = ROOT::Math::Similarity(diffParameters, recoTrackCovMatrix);
214  chi2 /= 5;
215 
216  LogTrace("TrackAssociator") << "====TP index=" << tpindex << " RECO TRACK index="<<tindex<<" WITH PT=" << (*rt)->pt() << "====\n"
217  << "qoverp sim: " << sParameters[0] << "\n"
218  << "lambda sim: " << sParameters[1] << "\n"
219  << "phi sim: " << sParameters[2] << "\n"
220  << "dxy sim: " << sParameters[3] << "\n"
221  << "dsz sim: " << sParameters[4] << "\n"
222  << ": " /*<< */ << "\n"
223  << "qoverp rec: " << (*rt)->qoverp() << " err: " << (*rt)->error(0) << "\n"
224  << "lambda rec: " << (*rt)->lambda() << " err: " << (*rt)->error(1) << "\n"
225  << "phi rec: " << (*rt)->phi() << " err: " << (*rt)->error(2) << "\n"
226  << "dxy rec: " << (*rt)->dxy(bs.position()) << " err: " << (*rt)->error(3) << "\n"
227  << "dsz rec: " << (*rt)->dsz(bs.position()) << " err: " << (*rt)->error(4) << "\n"
228  << ": " /*<< */ << "\n"
229  << "chi2: " << chi2 << "\n";
230 
231  if (chi2<chi2cut) {
232  outputCollection.insert(tC[tindex],
233  std::make_pair(edm::Ref<TrackingParticleCollection>(tPCH, tpindex),
234  -chi2));//-chi2 because the Association Map is ordered using std::greater
235  }
236  }
237  }
238  }
239  outputCollection.post_insert();
240  return outputCollection;
241 }
242 
243 
244 
245 SimToRecoCollection TrackAssociatorByChi2::associateSimToReco(const edm::RefToBaseVector<reco::Track>& tC,
246  const edm::RefVector<TrackingParticleCollection>& tPCH,
247  const edm::Event * e,
248  const edm::EventSetup *setup ) const {
249  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
250  e->getByLabel(bsSrc,recoBeamSpotHandle);
251  reco::BeamSpot bs = *recoBeamSpotHandle;
252 
253  SimToRecoCollection outputCollection;
254  double chi2;
255 
256  TrackingParticleCollection tPC;
257  if (tPCH.size()!=0) tPC = *tPCH.product();
258 
259  int tpindex =0;
260  for (TrackingParticleCollection::const_iterator tp=tPC.begin(); tp!=tPC.end(); tp++, ++tpindex){
261 
262  //skip tps with a very small pt
263  //if (sqrt(tp->momentum().perp2())<0.5) continue;
264 
265  LogDebug("TrackAssociator") << "=========LOOKING FOR ASSOCIATION===========" << "\n"
266  << "TrackingParticle #"<<tpindex<<" with pt=" << sqrt(tp->momentum().perp2()) << "\n"
267  << "===========================================" << "\n";
268 
269  Basic3DVector<double> momAtVtx(tp->momentum().x(),tp->momentum().y(),tp->momentum().z());
270  Basic3DVector<double> vert(tp->vertex().x(),tp->vertex().y(),tp->vertex().z());
271 
272  std::pair<bool,reco::TrackBase::ParameterVector> params = parametersAtClosestApproach(vert, momAtVtx, tp->charge(), bs);
273  if (params.first){
274  TrackBase::ParameterVector sParameters=params.second;
275 
276  int tindex=0;
277  for (RefToBaseVector<reco::Track>::const_iterator rt=tC.begin(); rt!=tC.end(); rt++, tindex++){
278 
279  TrackBase::ParameterVector rParameters = (*rt)->parameters();
280 
281  TrackBase::CovarianceMatrix recoTrackCovMatrix = (*rt)->covariance();
282  if (onlyDiagonal) {
283  for (unsigned int i=0;i<5;i++){
284  for (unsigned int j=0;j<5;j++){
285  if (i!=j) recoTrackCovMatrix(i,j)=0;
286  }
287  }
288  }
289 
290  recoTrackCovMatrix.Invert();
291 
292  TrackBase::ParameterVector diffParameters = rParameters - sParameters;
293  diffParameters[2] = reco::deltaPhi(diffParameters[2],0.f);
294 
295  chi2 = ROOT::Math::Similarity(recoTrackCovMatrix, diffParameters);
296  chi2 /= 5;
297 
298  LogTrace("TrackAssociator") << "====TP index=" << tpindex << " RECO TRACK index="<<tindex<<" WITH PT=" << (*rt)->pt() << "====\n"
299  << "qoverp sim: " << sParameters[0] << "\n"
300  << "lambda sim: " << sParameters[1] << "\n"
301  << "phi sim: " << sParameters[2] << "\n"
302  << "dxy sim: " << sParameters[3] << "\n"
303  << "dsz sim: " << sParameters[4] << "\n"
304  << ": " /*<< */ << "\n"
305  << "qoverp rec: " << (*rt)->qoverp() << " err: " << (*rt)->error(0) << "\n"
306  << "lambda rec: " << (*rt)->lambda() << " err: " << (*rt)->error(1) << "\n"
307  << "phi rec: " << (*rt)->phi() << " err: " << (*rt)->error(2) << "\n"
308  << "dxy rec: " << (*rt)->dxy(bs.position()) << " err: " << (*rt)->error(3) << "\n"
309  << "dsz rec: " << (*rt)->dsz(bs.position()) << " err: " << (*rt)->error(4) << "\n"
310  << ": " /*<< */ << "\n"
311  << "chi2: " << chi2 << "\n";
312 
313  if (chi2<chi2cut) {
314  outputCollection.insert(edm::Ref<TrackingParticleCollection>(tPCH, tpindex),
315  std::make_pair(tC[tindex],
316  -chi2));//-chi2 because the Association Map is ordered using std::greater
317  }
318  }
319  }
320  }
321  outputCollection.post_insert();
322  return outputCollection;
323 }
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:501
i
int i
Definition: DBlmapReader.cc:9
Basic3DVector::y
T y() const
Cartesian y coordinate.
Definition: newBasic3DVector.h:95
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:71
Basic3DVector::x
T x() const
Cartesian x coordinate.
Definition: newBasic3DVector.h:92
math::XYZTLorentzVectorD
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:15
TrackingParticleCollection
std::vector< TrackingParticle > TrackingParticleCollection
Definition: HiEvtPlaneFlatProducer.cc:116
edm::RefVector::product
C const * product() const
Accessor for product collection.
Definition: RefVector.h:272
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:68
edm::RefToBaseVector::end
const_iterator end() const
Definition: RefToBaseVector.h:230
GlobalPoint
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
reco::TrackCollection
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:10
PV3DBase::y
T y() const
Definition: PV3DBase.h:62
TrackAssociatorByChi2::parametersAtClosestApproach
std::pair< bool, reco::TrackBase::ParameterVector > parametersAtClosestApproach(Basic3DVector< double >, Basic3DVector< double >, float, const reco::BeamSpot &) const
propagate the track parameters of TrackinParticle from production vertex to the point of closest appr...
Definition: TrackAssociatorByChi2.cc:136
TrackAssociatorByChi2::RecoToSimPair
std::pair< reco::Track, Chi2SimMap > RecoToSimPair
Definition: TrackAssociatorByChi2.h:34
reco::TrackBase::ParameterVector
math::Vector< dimension >::type ParameterVector
parameter vector
Definition: TrackBase.h:70
FreeTrajectoryState::charge
TrackCharge charge() const
Definition: FreeTrajectoryState.h:81
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:105
DeDxDiscriminatorTools::charge
double charge(const std::vector< uint8_t > &Ampls)
Definition: DeDxDiscriminatorTools.cc:43
TrackAssociatorByChi2::compareTracksParam
double compareTracksParam(reco::TrackCollection::const_iterator, edm::SimTrackContainer::const_iterator, const math::XYZTLorentzVectorD, GlobalVector, reco::TrackBase::CovarianceMatrix, const reco::BeamSpot &) const
compare reco::TrackCollection and edm::SimTrackContainer iterators: returns the chi2 ...
PV3DBase::theta
Geom::Theta< T > theta() const
Definition: PV3DBase.h:74
TrackCharge
int TrackCharge
Definition: TrackCharge.h:4
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:66
Basic3DVector::z
T z() const
Cartesian z coordinate.
Definition: newBasic3DVector.h:98
edm::AssociationMap::post_insert
void post_insert()
post insert action
Definition: AssociationMap.h:133
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
PV3DBase::z
T z() const
Definition: PV3DBase.h:63
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
j
int j
Definition: DBlmapReader.cc:9
TrackAssociatorByChi2::associateSimToReco
reco::SimToRecoCollection associateSimToReco(const edm::RefToBaseVector< reco::Track > &, const edm::RefVector< TrackingParticleCollection > &, const edm::Event *event=0, const edm::EventSetup *setup=0) const
Association Sim To Reco with Collections.
Definition: TrackAssociatorByChi2.cc:245
f
double f[11][100]
Definition: MuScleFitUtils.cc:79
Geom::halfPi
double halfPi()
Definition: Pi.h:33
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:356
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:502
FreeTrajectoryState::momentum
GlobalVector momentum() const
Definition: FreeTrajectoryState.h:78
TrackAssociatorByChi2::associateRecoToSim
double associateRecoToSim(reco::TrackCollection::const_iterator, TrackingParticleCollection::const_iterator, const reco::BeamSpot &) const
compare reco::TrackCollection and TrackingParticleCollection iterators: returns the chi2 ...
reco::deltaPhi
double deltaPhi(double phi1, double phi2)
Definition: deltaPhi.h:12
FreeTrajectoryState::position
GlobalPoint position() const
Definition: FreeTrajectoryState.h:75
edm::AssociationMap::insert
void insert(const key_type &k, const data_type &v)
insert an association
Definition: AssociationMap.h:81
edm::SimVertexContainer
std::vector< SimVertex > SimVertexContainer
Definition: SimVertexContainer.h:13
edm::RefToBaseVector::begin
const_iterator begin() const
Definition: RefToBaseVector.h:222
TrackAssociatorByChi2::Chi2SimMap
std::map< double, SimTrack > Chi2SimMap
Definition: TrackAssociatorByChi2.h:33
TrackAssociatorByChi2::RecoToSimPairAssociation
std::vector< RecoToSimPair > RecoToSimPairAssociation
Definition: TrackAssociatorByChi2.h:35
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
edm::RefVector::size
size_type size() const
Size of the RefVector.
Definition: RefVector.h:89
reco::BeamSpot::position
const Point & position() const
position
Definition: BeamSpot.h:63
PV3DBase::x
T x() const
Definition: PV3DBase.h:61
edm::SimTrackContainer
std::vector< SimTrack > SimTrackContainer
Definition: SimTrackContainer.h:12
HcalObjRepresent::setup
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
Definition: HcalObjRepresent.h:628
v
mathSSE::Vec4< T > v
Definition: newBasic3DVector.h:354
reco::TrackBase::CovarianceMatrix
math::Error< dimension >::type CovarianceMatrix
5 parameter covariance matrix
Definition: TrackBase.h:72
GlobalVector
Global3DVector GlobalVector
Definition: GlobalVector.h:10
Generated for CMSSW Reference Manual by doxygen 1.8.5