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.
3 
9 
10 using namespace edm;
11 using namespace reco;
12 using namespace std;
13 
14 double TrackAssociatorByChi2::compareTracksParam ( TrackCollection::const_iterator rt,
15  SimTrackContainer::const_iterator st,
16  const math::XYZTLorentzVectorD vertexPosition,
17  GlobalVector magField,
18  TrackBase::CovarianceMatrix invertedCovariance,
19  const reco::BeamSpot& bs) const{
20 
21  Basic3DVector<double> momAtVtx(st->momentum().x(),st->momentum().y(),st->momentum().z());
22  Basic3DVector<double> vert = (Basic3DVector<double>) vertexPosition;
23 
24  std::pair<bool,reco::TrackBase::ParameterVector> params = parametersAtClosestApproach(vert, momAtVtx, st->charge(), bs);
25  if (params.first){
26  TrackBase::ParameterVector sParameters = params.second;
27  TrackBase::ParameterVector rParameters = rt->parameters();
28 
29  TrackBase::ParameterVector diffParameters = rParameters - sParameters;
30  diffParameters[2] = reco::deltaPhi(diffParameters[2],0.f);
31  double chi2 = ROOT::Math::Dot(diffParameters * invertedCovariance, diffParameters);
32 
33  return chi2;
34  } else {
35  return 10000000000.;
36  }
37 }
38 
39 
42  const SimTrackContainer& stColl,
43  const SimVertexContainer& svColl,
44  const reco::BeamSpot& bs) const{
45 
46  RecoToSimPairAssociation outputVec;
47 
48  for (TrackCollection::const_iterator track=rtColl.begin(); track!=rtColl.end(); track++){
49  Chi2SimMap outMap;
50 
51  TrackBase::ParameterVector rParameters = track->parameters();
52 
53  TrackBase::CovarianceMatrix recoTrackCovMatrix = track->covariance();
54  if (onlyDiagonal){
55  for (unsigned int i=0;i<5;i++){
56  for (unsigned int j=0;j<5;j++){
57  if (i!=j) recoTrackCovMatrix(i,j)=0;
58  }
59  }
60  }
61  recoTrackCovMatrix.Invert();
62 
63  for (SimTrackContainer::const_iterator st=stColl.begin(); st!=stColl.end(); st++){
64 
65  Basic3DVector<double> momAtVtx(st->momentum().x(),st->momentum().y(),st->momentum().z());
66  Basic3DVector<double> vert = (Basic3DVector<double>) svColl[st->vertIndex()].position();
67 
68  std::pair<bool,reco::TrackBase::ParameterVector> params = parametersAtClosestApproach(vert, momAtVtx, st->charge(), bs);
69  if (params.first){
70  TrackBase::ParameterVector sParameters = params.second;
71 
72  TrackBase::ParameterVector diffParameters = rParameters - sParameters;
73  diffParameters[2] = reco::deltaPhi(diffParameters[2],0.f);
74  double chi2 = ROOT::Math::Dot(diffParameters * recoTrackCovMatrix, diffParameters);
75  chi2/=5;
76  if (chi2<chi2cut) outMap[chi2]=*st;
77  }
78  }
79  outputVec.push_back(RecoToSimPair(*track,outMap));
80  }
81  return outputVec;
82 }
83 
85  TrackBase::CovarianceMatrix& recoTrackCovMatrix,
86  Basic3DVector<double>& momAtVtx,
88  int& charge,
89  const reco::BeamSpot& bs) const{
90 
91  double chi2;
92 
93  std::pair<bool,reco::TrackBase::ParameterVector> params = parametersAtClosestApproach(vert, momAtVtx, charge, bs);
94  if (params.first){
95  TrackBase::ParameterVector sParameters=params.second;
96 
97  TrackBase::ParameterVector diffParameters = rParameters - sParameters;
98  diffParameters[2] = reco::deltaPhi(diffParameters[2],0.f);
99  chi2 = ROOT::Math::Dot(diffParameters * recoTrackCovMatrix, diffParameters);
100  chi2 /= 5;
101 
102  LogDebug("TrackAssociator") << "====NEW RECO TRACK WITH PT=" << sin(rParameters[1])*float(charge)/rParameters[0] << "====\n"
103  << "qoverp sim: " << sParameters[0] << "\n"
104  << "lambda sim: " << sParameters[1] << "\n"
105  << "phi sim: " << sParameters[2] << "\n"
106  << "dxy sim: " << sParameters[3] << "\n"
107  << "dsz sim: " << sParameters[4] << "\n"
108  << ": " /*<< */ << "\n"
109  << "qoverp rec: " << rParameters[0] << "\n"
110  << "lambda rec: " << rParameters[1] << "\n"
111  << "phi rec: " << rParameters[2] << "\n"
112  << "dxy rec: " << rParameters[3] << "\n"
113  << "dsz rec: " << rParameters[4] << "\n"
114  << ": " /*<< */ << "\n"
115  << "chi2: " << chi2 << "\n";
116 
117  return chi2;
118  } else {
119  return 10000000000.;
120  }
121 }
122 
123 
124 double TrackAssociatorByChi2::associateRecoToSim( TrackCollection::const_iterator rt,
125  TrackingParticleCollection::const_iterator tp,
126  const reco::BeamSpot& bs) const{
127  TrackBase::ParameterVector rParameters = rt->parameters();
128  TrackBase::CovarianceMatrix recoTrackCovMatrix = rt->covariance();
129  if (onlyDiagonal){
130  for (unsigned int i=0;i<5;i++){
131  for (unsigned int j=0;j<5;j++){
132  if (i!=j) recoTrackCovMatrix(i,j)=0;
133  }
134  }
135  }
136 
137  recoTrackCovMatrix.Invert();
138  Basic3DVector<double> momAtVtx(tp->momentum().x(),tp->momentum().y(),tp->momentum().z());
139  Basic3DVector<double> vert(tp->vertex().x(),tp->vertex().y(),tp->vertex().z());
140  int charge = tp->charge();
141  return getChi2(rParameters,recoTrackCovMatrix,momAtVtx,vert,charge,bs);
142 }
143 
144 pair<bool,TrackBase::ParameterVector>
146  Basic3DVector<double> momAtVtx,
147  float charge,
148  const BeamSpot& bs) const{
149 
150  TrackBase::ParameterVector sParameters;
151  try {
152  FreeTrajectoryState ftsAtProduction(GlobalPoint(vertex.x(),vertex.y(),vertex.z()),
153  GlobalVector(momAtVtx.x(),momAtVtx.y(),momAtVtx.z()),
154  TrackCharge(charge),
155  theMF.product());
156  TSCBLBuilderNoMaterial tscblBuilder;
157  TrajectoryStateClosestToBeamLine tsAtClosestApproach = tscblBuilder(ftsAtProduction,bs);//as in TrackProducerAlgorithm
158 
159  GlobalPoint v = tsAtClosestApproach.trackStateAtPCA().position();
160  GlobalVector p = tsAtClosestApproach.trackStateAtPCA().momentum();
161  sParameters[0] = tsAtClosestApproach.trackStateAtPCA().charge()/p.mag();
162  sParameters[1] = Geom::halfPi() - p.theta();
163  sParameters[2] = p.phi();
164  sParameters[3] = (-v.x()*sin(p.phi())+v.y()*cos(p.phi()));
165  sParameters[4] = v.z()*p.perp()/p.mag() - (v.x()*p.x()+v.y()*p.y())/p.perp() * p.z()/p.mag();
166 
167  return pair<bool,TrackBase::ParameterVector>(true,sParameters);
168  } catch ( ... ) {
169  return pair<bool,TrackBase::ParameterVector>(false,sParameters);
170  }
171 }
172 
175  const edm::Event * e,
176  const edm::EventSetup *setup ) const{
177  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
178  e->getByLabel(bsSrc,recoBeamSpotHandle);
179  reco::BeamSpot bs = *recoBeamSpotHandle;
180 
181  RecoToSimCollection outputCollection;
182 
184  if (tPCH.size()!=0) tPC = *tPCH.product();
185 
186  int tindex=0;
187  for (RefToBaseVector<reco::Track>::const_iterator rt=tC.begin(); rt!=tC.end(); rt++, tindex++){
188 
189  LogDebug("TrackAssociator") << "=========LOOKING FOR ASSOCIATION===========" << "\n"
190  << "rec::Track #"<<tindex<<" with pt=" << (*rt)->pt() << "\n"
191  << "===========================================" << "\n";
192 
193  TrackBase::ParameterVector rParameters = (*rt)->parameters();
194 
195  TrackBase::CovarianceMatrix recoTrackCovMatrix = (*rt)->covariance();
196  if (onlyDiagonal){
197  for (unsigned int i=0;i<5;i++){
198  for (unsigned int j=0;j<5;j++){
199  if (i!=j) recoTrackCovMatrix(i,j)=0;
200  }
201  }
202  }
203 
204  recoTrackCovMatrix.Invert();
205 
206  int tpindex =0;
207  for (TrackingParticleCollection::const_iterator tp=tPC.begin(); tp!=tPC.end(); tp++, ++tpindex){
208 
209  //skip tps with a very small pt
210  //if (sqrt(tp->momentum().perp2())<0.5) continue;
211  int charge = tp->charge();
212  if (charge==0) continue;
213  Basic3DVector<double> momAtVtx(tp->momentum().x(),tp->momentum().y(),tp->momentum().z());
214  Basic3DVector<double> vert=(Basic3DVector<double>) tp->vertex();
215 
216  double chi2 = getChi2(rParameters,recoTrackCovMatrix,momAtVtx,vert,charge,bs);
217 
218  if (chi2<chi2cut) {
219  outputCollection.insert(tC[tindex],
220  std::make_pair(edm::Ref<TrackingParticleCollection>(tPCH, tpindex),
221  -chi2));//-chi2 because the Association Map is ordered using std::greater
222  }
223  }
224  }
225  outputCollection.post_insert();
226  return outputCollection;
227 }
228 
229 
232  const edm::Event * e,
233  const edm::EventSetup *setup ) const {
234  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
235  e->getByLabel(bsSrc,recoBeamSpotHandle);
236  reco::BeamSpot bs = *recoBeamSpotHandle;
237 
238  SimToRecoCollection outputCollection;
239 
241  if (tPCH.size()!=0) tPC = *tPCH.product();
242 
243  int tpindex =0;
244  for (TrackingParticleCollection::const_iterator tp=tPC.begin(); tp!=tPC.end(); tp++, ++tpindex){
245 
246  //skip tps with a very small pt
247  //if (sqrt(tp->momentum().perp2())<0.5) continue;
248  int charge = tp->charge();
249  if (charge==0) continue;
250 
251  LogDebug("TrackAssociator") << "=========LOOKING FOR ASSOCIATION===========" << "\n"
252  << "TrackingParticle #"<<tpindex<<" with pt=" << sqrt(tp->momentum().perp2()) << "\n"
253  << "===========================================" << "\n";
254 
255  Basic3DVector<double> momAtVtx(tp->momentum().x(),tp->momentum().y(),tp->momentum().z());
256  Basic3DVector<double> vert(tp->vertex().x(),tp->vertex().y(),tp->vertex().z());
257 
258  int tindex=0;
259  for (RefToBaseVector<reco::Track>::const_iterator rt=tC.begin(); rt!=tC.end(); rt++, tindex++){
260 
261  TrackBase::ParameterVector rParameters = (*rt)->parameters();
262  TrackBase::CovarianceMatrix recoTrackCovMatrix = (*rt)->covariance();
263  if (onlyDiagonal) {
264  for (unsigned int i=0;i<5;i++){
265  for (unsigned int j=0;j<5;j++){
266  if (i!=j) recoTrackCovMatrix(i,j)=0;
267  }
268  }
269  }
270  recoTrackCovMatrix.Invert();
271 
272  double chi2 = getChi2(rParameters,recoTrackCovMatrix,momAtVtx,vert,charge,bs);
273 
274  if (chi2<chi2cut) {
275  outputCollection.insert(edm::Ref<TrackingParticleCollection>(tPCH, tpindex),
276  std::make_pair(tC[tindex],
277  -chi2));//-chi2 because the Association Map is ordered using std::greater
278  }
279  }
280  }
281  outputCollection.post_insert();
282  return outputCollection;
283 }
284 
285 
286 
287 
290  const edm::Event * e,
291  const edm::EventSetup *setup ) const{
292  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
293  e->getByLabel(bsSrc,recoBeamSpotHandle);
294  reco::BeamSpot bs = *recoBeamSpotHandle;
295 
296  RecoToGenCollection outputCollection;
297 
299  if (tPCH.size()!=0) tPC = *tPCH.product();
300 
301  int tindex=0;
302  for (RefToBaseVector<reco::Track>::const_iterator rt=tC.begin(); rt!=tC.end(); rt++, tindex++){
303 
304  LogDebug("TrackAssociator") << "=========LOOKING FOR ASSOCIATION===========" << "\n"
305  << "rec::Track #"<<tindex<<" with pt=" << (*rt)->pt() << "\n"
306  << "===========================================" << "\n";
307 
308  TrackBase::ParameterVector rParameters = (*rt)->parameters();
309 
310  TrackBase::CovarianceMatrix recoTrackCovMatrix = (*rt)->covariance();
311  if (onlyDiagonal){
312  for (unsigned int i=0;i<5;i++){
313  for (unsigned int j=0;j<5;j++){
314  if (i!=j) recoTrackCovMatrix(i,j)=0;
315  }
316  }
317  }
318 
319  recoTrackCovMatrix.Invert();
320 
321  int tpindex =0;
322  for (GenParticleCollection::const_iterator tp=tPC.begin(); tp!=tPC.end(); tp++, ++tpindex){
323 
324  //skip tps with a very small pt
325  //if (sqrt(tp->momentum().perp2())<0.5) continue;
326  int charge = tp->charge();
327  if (charge==0) continue;
328  Basic3DVector<double> momAtVtx(tp->momentum().x(),tp->momentum().y(),tp->momentum().z());
329  Basic3DVector<double> vert=(Basic3DVector<double>) tp->vertex();
330 
331  double chi2 = getChi2(rParameters,recoTrackCovMatrix,momAtVtx,vert,charge,bs);
332 
333  if (chi2<chi2cut) {
334  outputCollection.insert(tC[tindex],
335  std::make_pair(edm::Ref<GenParticleCollection>(tPCH, tpindex),
336  -chi2));//-chi2 because the Association Map is ordered using std::greater
337  }
338  }
339  }
340  outputCollection.post_insert();
341  return outputCollection;
342 }
343 
344 
347  const edm::Event * e,
348  const edm::EventSetup *setup ) const {
349 
350  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
351  e->getByLabel(bsSrc,recoBeamSpotHandle);
352  reco::BeamSpot bs = *recoBeamSpotHandle;
353 
354  GenToRecoCollection outputCollection;
355 
357  if (tPCH.size()!=0) tPC = *tPCH.product();
358 
359  int tpindex =0;
360  for (GenParticleCollection::const_iterator tp=tPC.begin(); tp!=tPC.end(); tp++, ++tpindex){
361 
362  //skip tps with a very small pt
363  //if (sqrt(tp->momentum().perp2())<0.5) continue;
364  int charge = tp->charge();
365  if (charge==0) continue;
366 
367  LogDebug("TrackAssociator") << "=========LOOKING FOR ASSOCIATION===========" << "\n"
368  << "TrackingParticle #"<<tpindex<<" with pt=" << sqrt(tp->momentum().perp2()) << "\n"
369  << "===========================================" << "\n";
370 
371  Basic3DVector<double> momAtVtx(tp->momentum().x(),tp->momentum().y(),tp->momentum().z());
372  Basic3DVector<double> vert(tp->vertex().x(),tp->vertex().y(),tp->vertex().z());
373 
374  int tindex=0;
375  for (RefToBaseVector<reco::Track>::const_iterator rt=tC.begin(); rt!=tC.end(); rt++, tindex++){
376 
377  TrackBase::ParameterVector rParameters = (*rt)->parameters();
378  TrackBase::CovarianceMatrix recoTrackCovMatrix = (*rt)->covariance();
379  if (onlyDiagonal) {
380  for (unsigned int i=0;i<5;i++){
381  for (unsigned int j=0;j<5;j++){
382  if (i!=j) recoTrackCovMatrix(i,j)=0;
383  }
384  }
385  }
386  recoTrackCovMatrix.Invert();
387 
388  double chi2 = getChi2(rParameters,recoTrackCovMatrix,momAtVtx,vert,charge,bs);
389 
390  if (chi2<chi2cut) {
391  outputCollection.insert(edm::Ref<GenParticleCollection>(tPCH, tpindex),
392  std::make_pair(tC[tindex],
393  -chi2));//-chi2 because the Association Map is ordered using std::greater
394  }
395  }
396  }
397  outputCollection.post_insert();
398  return outputCollection;
399 }
#define LogDebug(id)
std::vector< GenParticle > GenParticleCollection
collection of GenParticles
int i
Definition: DBlmapReader.cc:9
T y() const
Cartesian y coordinate.
T perp() const
Definition: PV3DBase.h:72
T x() const
Cartesian x coordinate.
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:15
std::vector< TrackingParticle > TrackingParticleCollection
C const * product() const
Accessor for product collection.
Definition: RefVector.h:272
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
const_iterator end() const
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:10
T y() const
Definition: PV3DBase.h:63
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...
std::pair< reco::Track, Chi2SimMap > RecoToSimPair
math::Vector< dimension >::type ParameterVector
parameter vector
Definition: TrackBase.h:70
TrackCharge charge() const
double charge(const std::vector< uint8_t > &Ampls)
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 ...
Geom::Theta< T > theta() const
Definition: PV3DBase.h:75
int TrackCharge
Definition: TrackCharge.h:4
T mag() const
Definition: PV3DBase.h:67
T z() const
Cartesian z coordinate.
void post_insert()
post insert action
T sqrt(T t)
Definition: SSEVec.h:48
T z() const
Definition: PV3DBase.h:64
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
int j
Definition: DBlmapReader.cc:9
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.
double f[11][100]
double getChi2(reco::TrackBase::ParameterVector &rParameters, reco::TrackBase::CovarianceMatrix &recoTrackCovMatrix, Basic3DVector< double > &momAtVtx, Basic3DVector< double > &vert, int &charge, const reco::BeamSpot &) const
basic method where chi2 is computed
double halfPi()
Definition: Pi.h:33
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:361
GlobalVector momentum() const
double associateRecoToSim(reco::TrackCollection::const_iterator, TrackingParticleCollection::const_iterator, const reco::BeamSpot &) const
compare reco::TrackCollection and TrackingParticleCollection iterators: returns the chi2 ...
double deltaPhi(double phi1, double phi2)
Definition: deltaPhi.h:12
GlobalPoint position() const
void insert(const key_type &k, const data_type &v)
insert an association
reco::RecoToGenCollection associateRecoToGen(const edm::RefToBaseVector< reco::Track > &, const edm::RefVector< reco::GenParticleCollection > &, const edm::Event *event=0, const edm::EventSetup *setup=0) const
Association Sim To Reco with Collections (Gen Particle version)
std::vector< SimVertex > SimVertexContainer
reco::GenToRecoCollection associateGenToReco(const edm::RefToBaseVector< reco::Track > &, const edm::RefVector< reco::GenParticleCollection > &, const edm::Event *event=0, const edm::EventSetup *setup=0) const
Association Sim To Reco with Collections (Gen Particle version)
const_iterator begin() const
std::map< double, SimTrack > Chi2SimMap
std::vector< RecoToSimPair > RecoToSimPairAssociation
size_type size() const
Size of the RefVector.
Definition: RefVector.h:89
T x() const
Definition: PV3DBase.h:62
std::vector< SimTrack > SimTrackContainer
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
math::Error< dimension >::type CovarianceMatrix
5 parameter covariance matrix
Definition: TrackBase.h:72
Global3DVector GlobalVector
Definition: GlobalVector.h:10