CMS 3D CMS Logo

TrackHelixVertexFitter.cc
Go to the documentation of this file.
1 /* From SimpleFits Package
2  * Designed an written by
3  * author: Ian M. Nugent
4  * Humboldt Foundations
5  */
8 #include "TDecompBK.h"
9 #include <iostream>
10 
11 using namespace tauImpactParameter;
12 
13 TrackHelixVertexFitter::TrackHelixVertexFitter(const std::vector<TrackParticle>& particles, const TVector3& vguess)
14  : isFit_(false),
15  isConfigured_(false),
16  nParticles_(particles.size()),
17  nPar_((NFreeTrackPar-NFreeVertexPar)*particles.size()+NFreeVertexPar),
18  nVal_(TrackParticle::NHelixPar*particles.size())
19 {
21  par_.ResizeTo(nPar_);
22  parcov_.ResizeTo(nPar_,nPar_);
23  val_.ResizeTo(nVal_);
24  cov_.ResizeTo(nVal_,nVal_);
25  for(unsigned int p=0; p<particles.size();p++){
26  for(unsigned int j=0; j<TrackParticle::NHelixPar;j++){
27  val_(measuredValueIndex(j,p))=particles[p].parameter(j);
28  for(unsigned int k=0; k<TrackParticle::NHelixPar;k++){
29  cov_(measuredValueIndex(j,p),measuredValueIndex(k,p))=particles[p].covariance(j,k);
30  }
31  }
32  }
33  TDecompBK Inverter(cov_);
34  double det = cov_.Determinant();
35  if(!Inverter.Decompose()){
36  edm::LogWarning("TrackHelixVertexFitter::TrackHelixVertexFitter") << "Fit failed: unable to invert SYM gain matrix " << det << " \n" << std::endl;
37  return;
38  }
39 
40  cov_inv_.ResizeTo(nVal_,nVal_);
41  cov_inv_=Inverter.Invert();
43  // Set Initial conditions within reason
44  par_(x0) = vguess.X(); parcov_(x0,x0)=1.0;
45  par_(y0) = vguess.Y(); parcov_(y0,y0)=1.0;
46  par_(z0) = vguess.Z(); parcov_(z0,z0)=1.0;
47  for(unsigned int p=0; p<particles_.size();p++){
51 
55  }
56  isConfigured_=true;
57 }
58 
60 
61 double TrackHelixVertexFitter::updateChisquare(const TVectorT<double>& inpar){
62  TVectorT<double> vprime=computePar(inpar);
63  TVectorT<double> dalpha=vprime-val_;
64  double c2=dalpha*(cov_inv_*dalpha);
65  return c2;
66 }
67 
68 std::vector<TrackParticle> TrackHelixVertexFitter::getRefitTracks(){
69  std::vector<TrackParticle> refitParticles;
70  for(unsigned int p=0;p<particles_.size();p++){
71  TVectorT<double> FreePar(NFreeTrackPar);
72  TMatrixTSym<double> FreeParCov(NFreeTrackPar);
73  for(int i=0;i<FreeParCov.GetNrows();i++){
74  FreePar(i)=par_(freeParIndex(i,p));
75  for(int j=0;j<FreeParCov.GetNrows();j++){
76  FreeParCov(i,j)=parcov_(freeParIndex(i,p),freeParIndex(j,p));
77  }
78  }
79  TVectorT<double> TrackPar=computePar(FreePar);
80  TMatrixTSym<double> TrackCov=ErrorMatrixPropagator::propagateError(&TrackHelixVertexFitter::computePar,FreePar,FreeParCov);
81  refitParticles.push_back(TrackParticle(TrackPar,TrackCov,particles_[p].pdgId(),particles_[p].mass(),particles_[p].charge(),particles_[p].bField()));
82  }
83  return refitParticles;
84 }
85 
86 std::vector<LorentzVectorParticle> TrackHelixVertexFitter::getRefitLorentzVectorParticles(){
87  std::vector<LorentzVectorParticle> refitParticles;
88  for(unsigned int p=0;p<particles_.size();p++){
89  TVectorT<double> FreePar(NFreeTrackPar+NExtraPar+MassOffSet);
90  TMatrixTSym<double> FreeParCov(NFreeTrackPar+NExtraPar+MassOffSet);
91  for(int i=0;i<NFreeTrackPar;i++){
92  FreePar(i)=par_(freeParIndex(i,p));
93  for(int j=0;j<NFreeTrackPar;j++){
94  FreeParCov(i,j)=parcov_(freeParIndex(i,p),freeParIndex(j,p));
95  }
96  }
97  FreePar(NFreeTrackPar+MassOffSet)=particles_[p].mass();
98  FreePar(NFreeTrackPar+BField0)=particles_[p].bField();
99  TVectorT<double> LVPar=computeLorentzVectorPar(FreePar);
100  TMatrixTSym<double> LVCov=ErrorMatrixPropagator::propagateError(&TrackHelixVertexFitter::computeLorentzVectorPar,FreePar,FreeParCov);
101  refitParticles.push_back(LorentzVectorParticle(LVPar,LVCov,particles_[p].pdgId(),particles_[p].charge(),particles_[p].bField()));
102  }
103  return refitParticles;
104 }
105 
107  double c(0),b(0);
108  TVectorT<double> FreePar(par_.GetNrows()+NExtraPar+particles_.size());
109  TMatrixTSym<double> FreeParCov(par_.GetNrows()+NExtraPar+particles_.size());
110  for(int i=0;i<par_.GetNrows();i++){
111  FreePar(i)=par_(i);
112  for(int j=0;j<par_.GetNrows();j++){FreeParCov(i,j)=parcov_(i,j);}
113  }
114  for(unsigned int p=0; p<particles_.size();p++){
115  b=particles_[p].bField();
116  c+=particles_[p].charge();
117  FreePar(par_.GetNrows()+MassOffSet+p)=particles_[p].mass();
118  }
119  FreePar(par_.GetNrows()+BField0)=b;
120  TVectorT<double> mpar=computeMotherLorentzVectorPar(FreePar);
122  return LorentzVectorParticle(mpar,mcov,pdgid,c,b);
123 }
124 
126  return TVector3(par_(freeParIndex(x0,0)),par_(freeParIndex(y0,0)),par_(freeParIndex(z0,0)));
127 }
128 
130  TMatrixTSym<double> c(NFreeVertexPar);
131  for(unsigned int i=0;i<NFreeVertexPar;i++){
132  for(unsigned int j=0;j<NFreeVertexPar;j++){c(freeParIndex(i,0),freeParIndex(j,0))=parcov_(freeParIndex(i,0),freeParIndex(j,0));}
133  }
134  return c;
135 }
136 
137 void TrackHelixVertexFitter::computedxydz(const TVectorT<double>& inpar,int p,double& kappa,double& lam,double& phi,double& x,double& y,double& z,double& s,double& dxy,double& dz){
138  kappa=inpar(freeParIndex(kappa0,p));
139  lam=inpar(freeParIndex(lambda0,p));
140  phi=inpar(freeParIndex(phi0,p));
141  x=inpar(freeParIndex(x0,p));
142  y=inpar(freeParIndex(y0,p));
143  z=inpar(freeParIndex(z0,p));
144  double v=(2.0*kappa*(x*cos(phi)+y*sin(phi)));
145  double arcsinv=0;
146  if(v>=1.0){arcsinv=TMath::Pi()/2;}
147  else if(v<=-1.0){arcsinv=-TMath::Pi()/2;}
148  else{arcsinv=asin(v);}
149  s=1.0/(2.0*kappa)*arcsinv;
150  dxy=y*cos(phi)-x*sin(phi)-(1/kappa)*sin(kappa*s)*sin(kappa*s);
151  dz=z-s*tan(lam);
152 }
153 
154 TVectorT<double> TrackHelixVertexFitter::computePar(const TVectorT<double>& inpar){
155  int nparticles=(inpar.GetNrows()-NFreeVertexPar)/(NFreeTrackPar-NFreeVertexPar);
156  TVectorT<double> helices(nparticles*TrackParticle::NHelixPar);
157  for(int p=0;p<nparticles;p++){
158  TVectorT<double> TrackPar=computeTrackPar(inpar,p);
159  for(int i=0;i<TrackParticle::NHelixPar;i++){helices(measuredValueIndex(i,p))=TrackPar(i);}
160  }
161  return helices;
162 }
163 
164 TVectorT<double> TrackHelixVertexFitter::computeTrackPar(const TVectorT<double>& inpar, int p){
165  TVectorT<double> helix(TrackParticle::NHelixPar);
166  // copy parameters that are 1 to 1
167  double kappa,lam,phi,x,y,z,s,dxy,dz;
168  TrackHelixVertexFitter::computedxydz(inpar,p,kappa,lam,phi,x,y,z,s,dxy,dz);
169  helix(TrackParticle::kappa) = kappa;
170  helix(TrackParticle::lambda) = lam;
171  helix(TrackParticle::phi) = phi;
172  helix(TrackParticle::dxy) = dxy;
173  helix(TrackParticle::dz) = dz;
174  return helix;
175 }
176 
177 TVectorT<double> TrackHelixVertexFitter::computeLorentzVectorPar(const TVectorT<double>& inpar){
178  int np(0), parsize(0); parSizeInfo(inpar,np,parsize,true);
179  double B=inpar(parsize+BField0);
180  double massHypothesis=inpar(parsize+MassOffSet);
182  double kappa,lam,phi,x,y,z,s,dxy,dz;
183  int p=0;
184  TrackHelixVertexFitter::computedxydz(inpar,p,kappa,lam,phi,x,y,z,s,dxy,dz);
185  double phi1 = 2*s*kappa+phi;
186  double bOverK = B*(1.0/fabs(kappa));
187  LV(LorentzVectorParticle::px) = bOverK*cos(phi1);
188  LV(LorentzVectorParticle::py) = bOverK*sin(phi1);
189  LV(LorentzVectorParticle::pz) = bOverK*tan(lam) ;
194  return LV;
195 }
196 
197 TVectorT<double> TrackHelixVertexFitter::computeMotherLorentzVectorPar(const TVectorT<double>& inpar){
198  TVectorT<double> mother(LorentzVectorParticle::NLorentzandVertexPar);
199  double E(0);
200  int np(0), parsize(0); parSizeInfo(inpar,np,parsize,true);
201  for(int p=0;p<np;p++){
202  TVectorT<double> particlepar(NFreeTrackPar+NExtraPar+MassOffSet);
203  for(int i=0;i<NFreeTrackPar;i++){particlepar(i)=inpar(freeParIndex(i,p));}
204  particlepar(NFreeTrackPar+BField0)=inpar(parsize+BField0);
205  particlepar(NFreeTrackPar+MassOffSet)=inpar(parsize+MassOffSet+p);
206  TVectorT<double> daughter=TrackHelixVertexFitter::computeLorentzVectorPar(particlepar);
217  }
218  double P2=(mother(LorentzVectorParticle::px)*mother(LorentzVectorParticle::px)+
221  mother(LorentzVectorParticle::m)=(E*E-P2)/sqrt(fabs(E*E-P2));
222  return mother;
223 }
224 
226  int p(0);
227  if(Par==x0) return "x0";
228  if(Par==y0) return "y0";
229  if(Par==z0) return "z0";
230  for(p=0;p<nParticles_;p++){
231  if((Par-NFreeVertexPar)<(p+1)*(NFreeTrackPar-NFreeVertexPar))break;
232  }
233  TString n;
234  int index=Par-p*(NFreeTrackPar-NFreeVertexPar);
235  if(index==kappa0) n="kappa0";
236  if(index==lambda0) n="lambda0";
237  if(index==phi0) n="phi0";
238  n+="_particle";n+=p;
239  return n;
240 }
241 
242 void TrackHelixVertexFitter::parSizeInfo(const TVectorT<double>& inpar, int& np, int& parsize, bool hasextras){
243  if(hasextras)np=(inpar.GetNrows()-NFreeVertexPar-NExtraPar)/(NFreeTrackPar+MassOffSet-NFreeVertexPar);
244  else np=(inpar.GetNrows()-NFreeVertexPar)/(NFreeTrackPar-NFreeVertexPar);
246 }
size
Write out results.
const double Pi
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
static TVectorT< double > computeMotherLorentzVectorPar(const TVectorT< double > &inpar)
virtual double updateChisquare(const TVectorT< double > &inpar)
virtual std::vector< TrackParticle > getRefitTracks()
T x() const
Cartesian x coordinate.
math::XYZTLorentzVectorD LV
int np
Definition: AMPTWrapper.h:33
static TVectorT< double > computePar(const TVectorT< double > &inpar)
T sqrt(T t)
Definition: SSEVec.h:18
static TVectorT< double > computeLorentzVectorPar(const TVectorT< double > &inpar)
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
static int measuredValueIndex(int TrackPar, int Particle)
int k[5][pyjets_maxn]
virtual TMatrixTSym< double > getVertexError()
std::vector< LorentzVectorParticle > getRefitLorentzVectorParticles()
TrackHelixVertexFitter(const std::vector< TrackParticle > &particles, const TVector3 &vguess)
static TVectorT< double > computeTrackPar(const TVectorT< double > &inpar, int p=0)
static int freeParIndex(int Par, int Particle)
double b
Definition: hdecay.h:120
static void computedxydz(const TVectorT< double > &inpar, int particle, double &kappa, double &lam, double &phi, double &x, double &y, double &z, double &s, double &dxy, double &dz)
static void parSizeInfo(const TVectorT< double > &inpar, int &np, int &parsize, bool hasextras=0)
static const G4double kappa
static TMatrixTSym< double > propagateError(std::function< TVectorT< double >(const TVectorT< double > &)> f, const TVectorT< double > &inPar, TMatrixTSym< double > &inCov, double epsilon=0.001, double errorEpsilonRatio=1000)