CMS 3D CMS Logo

MuonErrorMatrixAdjuster.cc
Go to the documentation of this file.
2 
3 #include "TString.h"
4 #include "TMath.h"
7 
10 
12 
15 
18 
20 {
21  theCategory="MuonErrorMatrixAdjuster";
22  theInstanceName = iConfig.getParameter<std::string>("instanceName");
23  //register your products
24  produces<reco::TrackCollection>( theInstanceName);
25  produces<TrackingRecHitCollection>();
26  produces<reco::TrackExtraCollection>();
27 
28  theTrackLabel = iConfig.getParameter<edm::InputTag>("trackLabel");
29  theRescale = iConfig.getParameter<bool>("rescale");
30 
31  theMatrixProvider_pset = iConfig.getParameter<edm::ParameterSet>("errorMatrix_pset");
32 }
33 
34 
36 {
37 
38  // do anything here that needs to be done at desctruction time
39  // (e.g. close files, deallocate resources etc.)
40 
41 }
42 
43 
44 //
45 // member functions
46 //
47 
48 //take the error matrix and rescale it or just replace it
50 {
51  //CovarianceMatrix is template for SMatrix
52  reco::TrackBase::CovarianceMatrix revised_matrix( theMatrixProvider->get(momentum));
53 
54  if( theRescale){
55  //rescale old error matrix up by a factor
56  multiply(revised_matrix,error_matrix);
57  }
58  return revised_matrix;
59 }
60 
62  //scale term by term the matrix
63  // the true type of the matrix is such that [i][j] is the same memory object as [j][i]: looping i:0-5, j:0-5 double multiply the terms
64  // need to loop only on i:0-5, j:i-5
65  for(int i = 0;i!=5;i++){for(int j = i;j!=5;j++){
66  revised_matrix(i,j)*=scale_matrix(i,j);
67  }}
68 }
70  //divide term by term the matrix
71  // the true type of the matrix is such that [i][j] is the same memory object as [j][i]: looping i:0-5, j:0-5 double multiply the terms
72  // need to loop only on i:0-5, j:i-5
73  for(int i = 0;i!=5;i++){for(int j = i;j!=5;j++){
74  if (denom_matrix(i,j)==0) return false;
75  num_matrix(i,j)/=denom_matrix(i,j);
76  }}
77  return true;
78 }
79 
81  const FreeTrajectoryState & PCAstate){
82  //get the parameters of the track so I can reconstruct it
83  double chi2 = recotrack_orig.chi2();
84  double ndof = recotrack_orig.ndof();
85  const math::XYZPoint& refpos = recotrack_orig.referencePoint();
86  const math::XYZVector& mom = recotrack_orig.momentum();
87  int charge = recotrack_orig.charge();
88 
89  reco::TrackBase::CovarianceMatrix covariance_matrix = fix_cov_matrix(recotrack_orig.covariance(),PCAstate.momentum());
90 
91  LogDebug(theCategory)<<"chi2: "<<chi2
92  <<"\n ndof: "<<ndof
93  <<"\n refpos: "<<refpos
94  <<"\n mom: "<<mom
95  <<"\n charge: "<<charge
96  <<"\n covariance:\n"<<recotrack_orig.covariance()
97  <<"\n replaced by:\n"<<covariance_matrix;
98 
99  return reco::Track(chi2,ndof,refpos,mom,charge,covariance_matrix);
100 }
101 
103  reco::Track & recotrack,
105  //get the 5x5 matrix of recotrack/recotrack_orig
106  reco::TrackBase::CovarianceMatrix scale_matrix(recotrack.covariance());
107  if (!divide(scale_matrix,recotrack_orig.covariance())){
108  edm::LogError( theCategory)<<"original track error matrix has term ==0... skipping.";
109  return nullptr; }
110 
111  const reco::TrackExtraRef & trackExtra_orig = recotrack_orig.extra();
112  if (trackExtra_orig.isNull()) {
113  edm::LogError( theCategory)<<"original track has no track extra... skipping.";
114  return nullptr;}
115 
116  //copy the outer state. rescaling the error matrix
117  reco::TrackBase::CovarianceMatrix outerCov(trackExtra_orig->outerStateCovariance());
118  multiply(outerCov,scale_matrix);
119 
120  //copy the inner state, rescaling the error matrix
121  reco::TrackBase::CovarianceMatrix innerCov(trackExtra_orig->innerStateCovariance());
122  multiply(innerCov,scale_matrix);
123 
124  //put the trackExtra
125  TEcol.push_back(reco::TrackExtra (trackExtra_orig->outerPosition(), trackExtra_orig->outerMomentum(), true,
126  trackExtra_orig->innerPosition(), trackExtra_orig->innerMomentum(), true,
127  outerCov, trackExtra_orig->outerDetId(),
128  innerCov, trackExtra_orig->innerDetId(),
129  trackExtra_orig->seedDirection()));
130 
131  //add a reference to the trackextra on the track
133 
134  //return the reference to the last inserted then
135  return &(TEcol.back());
136 }
137 
138 
140  reco::Track & recotrack,
141  reco::TrackExtra & trackextra,
143  const TrackerTopology& ttopo){
144  //loop over the hits of the original track
145  trackingRecHit_iterator recHit = recotrack_orig.recHitsBegin();
146  auto const firstHitIndex = theRHi;
147  for (; recHit!=recotrack_orig.recHitsEnd();++recHit){
148  //clone it. this is meandatory
149  TrackingRecHit * hit = (*recHit)->clone();
150 
151  //put it on the new track
152  recotrack.appendHitPattern(*hit, ttopo);
153  //copy them in the new collection
154  RHcol.push_back(hit);
155  ++theRHi;
156 
157  }//loop over original rechits
158  //do something with the trackextra
159  trackextra.setHits(theRefprodRH, firstHitIndex, theRHi-firstHitIndex);
160 
161  return true; //if nothing fails
162 }
163 
164 bool MuonErrorMatrixAdjuster::selectTrack(const reco::Track & recotrack_orig){return true;}
165 
166 // ------------ method called to produce the data ------------
167 void
169 {
170  using namespace edm;
171 
172  //open a collection of track
174  iEvent.getByLabel( theTrackLabel,tracks);
175  LogDebug( theCategory)<<"considering: "<<tracks->size()<<" uncorrected reco::Track from the event.("<< theTrackLabel<<")";
176 
177  //get the mag field
178  iSetup.get<IdealMagneticFieldRecord>().get( theField);
179 
181  iSetup.get<TrackerTopologyRcd>().get(httopo);
182  const TrackerTopology& ttopo = *httopo;
183 
184  //prepare the output collection
185  auto Toutput = std::make_unique<reco::TrackCollection>();
186  auto TRHoutput = std::make_unique<TrackingRecHitCollection>();
187  auto TEoutput = std::make_unique<reco::TrackExtraCollection>();
189  theTEi=0;
191  theRHi=0;
192 
193 
194 
195  for(unsigned int it=0;it!=tracks->size();it++)
196  {
197  const reco::Track & recotrack_orig = (*tracks)[it];
199  if (PCAstate.position().mag()==0)
200  {edm::LogError( theCategory)<<"invalid state from track initial state in "<< theTrackLabel <<". skipping.";continue;}
201 
202  //create a reco::Track
203  reco::Track recotrack = makeTrack(recotrack_orig,PCAstate);
204 
205  //make a selection on the create reco::Track
206  if (!selectTrack(recotrack)) continue;
207 
208  Toutput->push_back(recotrack);
209  reco::Track & recotrackref = Toutput->back();
210 
211  //build the track extra
212  reco::TrackExtra * extra = makeTrackExtra(recotrack_orig,recotrackref,*TEoutput);
213  if (!extra){
214  edm::LogError( theCategory)<<"cannot create the track extra for this track.";
215  //pop the inserted track
216  Toutput->pop_back();
217  continue;}
218 
219  //attach the collection of rechits
220  if (!attachRecHits(recotrack_orig,recotrackref,*extra,*TRHoutput, ttopo)){
221  edm::LogError( theCategory)<<"cannot attach any rechits on this track";
222  //pop the inserted track
223  Toutput->pop_back();
224  //pop the track extra
225  TEoutput->pop_back();
226  theTEi--;
227  continue;}
228 
229  }//loop over the original reco tracks
230 
231  LogDebug( theCategory)<<"writing: "<<Toutput->size()<<" corrected reco::Track to the event.";
232 
233  iEvent.put(std::move(Toutput), theInstanceName);
234  iEvent.put(std::move(TEoutput));
235  iEvent.put(std::move(TRHoutput));
236 
237 }
238 
239 // ------------ method called once each job just before starting event loop ------------
240 void
242 {
244 }
245 
246 // ------------ method called once each job just after ending the event loop ------------
247 void
250 }
#define LogDebug(id)
void produce(edm::Event &, const edm::EventSetup &) override
T getParameter(std::string const &) const
const Point & referencePoint() const
Reference point on the track.
Definition: TrackBase.h:681
void beginJob() override
framework method
edm::RefProd< reco::TrackExtraCollection > theRefprodTE
get reference before put track extra to the event, in order to create edm::Ref
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:137
const TrackExtraRef & extra() const
reference to "extra" object
Definition: Track.h:189
reco::TrackBase::CovarianceMatrix fix_cov_matrix(const reco::TrackBase::CovarianceMatrix &error_matrix, const GlobalVector &momentum)
return a corrected error matrix
void setHits(TrackingRecHitRefProd const &prod, unsigned firstH, unsigned int nH)
CurvilinearTrajectoryError get(GlobalVector momentum, bool convolute=true)
main method to be used. Retrieve a 5x5 symetrical matrix according to parametrization of error or sca...
std::string theInstanceName
instrance name of the created track collecion. rechit and trackextra have no instance name ...
const Vector & momentum() const
track momentum vector
Definition: TrackBase.h:675
edm::InputTag theTrackLabel
input tag of the reco::Track collection to be corrected
bool attachRecHits(const reco::Track &recotrack_orig, reco::Track &recotrack, reco::TrackExtra &trackextra, TrackingRecHitCollection &RHcol, const TrackerTopology &ttopo)
attached rechits to the newly created reco::Track and reco::TrackExtra
std::string theCategory
log category: MuonErrorMatrixAdjuster
void multiply(reco::TrackBase::CovarianceMatrix &revised_matrix, const reco::TrackBase::CovarianceMatrix &scale_matrix)
mutliply revised_matrix (first argument) by second matrix TERM by TERM
edm::ESHandle< MagneticField > theField
hold on to the magnetic field
MuonErrorMatrixAdjuster(const edm::ParameterSet &)
constructor
void push_back(D *&d)
Definition: OwnVector.h:290
MuonErrorMatrix * theMatrixProvider
int iEvent
Definition: GenABIO.cc:230
T mag() const
Definition: PV3DBase.h:67
double chi2() const
chi-squared of the fit
Definition: TrackBase.h:549
bool divide(reco::TrackBase::CovarianceMatrix &num_matrix, const reco::TrackBase::CovarianceMatrix &denom_matrix)
divide the num_matrix (first argument) by second matrix, TERM by TERM
double ndof() const
number of degrees of freedom of the fit
Definition: TrackBase.h:555
CovarianceMatrix covariance() const
return track covariance matrix
Definition: TrackBase.h:731
~MuonErrorMatrixAdjuster() override
destructor
edm::Ref< reco::TrackExtraCollection >::key_type theTEi
edm::RefProd< TrackingRecHitCollection > theRefprodRH
get reference before put rechit to the event, in order to create edm::Ref
trackingRecHit_iterator recHitsBegin() const
Iterator to first hit on the track.
Definition: Track.h:104
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:535
bool isNull() const
Checks for null.
Definition: Ref.h:250
GlobalVector momentum() const
RefProd< PROD > getRefBeforePut()
Definition: Event.h:167
edm::ParameterSet theMatrixProvider_pset
holds the error matrix parametrization
GlobalPoint position() const
std::vector< TrackExtra > TrackExtraCollection
collection of TrackExtra objects
Definition: TrackExtraFwd.h:11
void setExtra(const TrackExtraRef &ref)
set reference to "extra" object
Definition: Track.h:184
reco::TrackExtra * makeTrackExtra(const reco::Track &recotrack_orig, reco::Track &recotrack, reco::TrackExtraCollection &TEcol)
create a track extra for the newly created recotrack, scaling the outer/inner measurment error matrix...
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
bool selectTrack(const reco::Track &recotrack_orig)
make a selection on the reco:Track. (dummy for the moment)
bool appendHitPattern(const TrackingRecHit &hit, const TrackerTopology &ttopo)
append a single hit to the HitPattern
Definition: TrackBase.h:456
reco::Track makeTrack(const reco::Track &recotrack_orig, const FreeTrajectoryState &PCAstate)
create a corrected reco::Track from itself and trajectory state (redundant information) ...
bool theRescale
select the rescaling or replacing method to correct the error matrix
HLT enums.
T get() const
Definition: EventSetup.h:63
int charge() const
track electric charge
Definition: TrackBase.h:567
FreeTrajectoryState initialFreeState(const reco::Track &tk, const MagneticField *field, bool withErr=true)
T const * product() const
Definition: ESHandle.h:86
def move(src, dest)
Definition: eostools.py:510
math::Error< dimension >::type CovarianceMatrix
5 parameter covariance matrix
Definition: TrackBase.h:77
edm::Ref< TrackingRecHitCollection >::key_type theRHi
trackingRecHit_iterator recHitsEnd() const
Iterator to last hit on the track.
Definition: Track.h:109