27 double& chirz,
double& chirf)
const {
28 double pi = 4.*atan(1.);
29 double twopi=8.*atan(1.);
32 std::cout<<
" HICMuonUpdator::update:: can not start::initial tsos is not valid " <<std::endl;
39 std::cout<<
" HICMuonUpdator::update::MTM size "<<MTM.size()<<
" vertex "<<zvert<<std::endl;
40 std::cout<<
" HICMuonUpdator::update::charge from trajectory"<<(MTM.back()).updatedState().freeTrajectoryState()->parameters().charge()<<std::endl;
42 std::cout<<
" HICMuonUpdator::update::momentum "<<(MTM.back()).updatedState().freeTrajectoryState()->parameters().momentum()<<std::endl;
44 vector<double> phihit,rhit,zhit,dphihit,drhit,dzhit,dzhitl,ehitphi,dehitphi,ehitstrip;
54 GlobalVector pold=(MTM.back()).updatedState().freeTrajectoryState()->parameters().momentum();
57 for(vector<TrajectoryMeasurement>::const_iterator ihit=MTM.begin();ihit!=MTM.end();ihit++){
60 phihit.push_back((*ihit).recHit()->globalPosition().phi());
61 rhit.push_back((*ihit).recHit()->globalPosition().perp());
62 zhit.push_back((*ihit).recHit()->globalPosition().z());
63 GlobalPoint realhit = (*ihit).recHit()->globalPosition();
65 double phierror=
sqrt((*ihit).recHit()->globalPositionError().phierr(realhit));
67 if(fabs(phierror)<0.0000001) {
70 ehitphi.push_back(phierror);
72 #ifdef UPDATOR_BARREL_DEBUG
73 cout<<
" Errors "<<phierror<<
" "<<(*ihit).recHit()->globalPositionError().rerr(realhit)<<
" "<<
tan(theta)<<endl;
78 ehitstrip.push_back(
sqrt((*ihit).recHit()->globalPositionError().czz()));
80 ehitstrip.push_back(
sqrt((*ihit).recHit()->globalPositionError().rerr(realhit)/fabs(
tan(theta))));
85 phihit.push_back(nRecHit->globalPosition().phi());
86 rhit.push_back(nRecHit->globalPosition().perp());
87 zhit.push_back(nRecHit->globalPosition().z());
88 ehitphi.push_back(
sqrt(nRecHit->globalPositionError().phierr(nRecHit->globalPosition())));
95 ehitstrip.push_back(
sqrt(nRecHit->globalPositionError().czz()));
97 ehitstrip.push_back(
sqrt(nRecHit->globalPositionError().rerr(nRecHit->globalPosition()))/fabs(
tan(theta)));
102 rhit.push_back(rvert);
103 zhit.push_back(zvert);
104 ehitstrip.push_back(ezvert);
106 for(vector<double>::const_iterator iphi=phihit.begin();iphi!=phihit.end()-1;iphi++){
107 double dpnew=fabs(*iphi-*(iphi+1));
108 if(dpnew>pi) dpnew=twopi-dpnew;
110 #ifdef UPDATOR_BARREL_DEBUG
111 cout<<
" dphi=dpnew="<<dpnew<<
" "<<*iphi<<
" "<<*(iphi+1)<<endl;
114 dphihit.push_back(dpnew);
117 for(vector<double>::const_iterator ir=rhit.begin();ir!=rhit.end()-2;ir++){
118 double dpnew=fabs(*ir-*(ir+1));
120 #ifdef UPDATOR_BARREL_DEBUG
121 cout<<
" dr=dpnew="<<dpnew<<
" "<<*ir<<
" "<<*(ir+1)<<endl;
124 drhit.push_back(dpnew);
126 for(vector<double>::const_iterator iz=zhit.begin();iz!=zhit.end()-2;iz++){
127 double dpnew=*iz-*(iz+1);
129 #ifdef UPDATOR_BARREL_DEBUG
130 cout<<
" dZ=dpnew="<<dpnew<<
" "<<*iz<<
" "<<*(iz+1)<<endl;
133 dzhit.push_back(fabs(dpnew));
134 dzhitl.push_back(dpnew);
137 dzhitl.push_back(*(zhit.end()-1)-zvert);
139 for(vector<double>::const_iterator ie=ehitphi.begin();ie!=ehitphi.end()-1;ie++){
140 double dpnew=(*ie)*1.14;
141 dehitphi.push_back(dpnew);
150 nTsos, chirz, chirf, tType);
154 std::cout<<
" Trajectory is not valid "<<std::endl;
164 nTsos, chirz, chirf, tType);
171 const vector <double>& err,
double& co1,
double& co2,
172 double & chi2)
const{
185 cout<<
" linefit2::sizes="<<x.size()<<
" "<<y.size()<<
" "<<err.size()<<endl;
186 cout<<
" linefit2::a00 "<<a00<<
" "<<a01<<
" "<<a11<<
" "<<b0<<
" "<<b1<<endl;
189 if(x.size() != y.size())
return fit;
190 if(x.size() != err.size())
return fit;
193 for (
unsigned int i=0;
i<x.size();
i++){
195 cout<<
" line2fit "<<a00<<
" "<<x[
i]/err[
i]<<
" "<<x[
i]<<
" "<<err[
i]<<
" second try "<<err[
i]<<endl;
197 a00=a00+(x[
i]/err[
i])*(x[
i]/err[
i]);
198 a01=a01+(x[
i]/err[
i])/err[i];
200 a11=a11+1./(err[
i]*err[
i]);
201 b0=b0+(x[
i]/err[
i])*(y[i]/err[i]);
202 b1=b1+(y[
i]/err[
i])/err[i];
205 cout<<
"linefit2="<<x[
i]<<
" "<<y[
i]<<
" "<<err[
i]<<
" "<<a00<<
" "<<a01<<
" "<<a11<<
" "<<b0<<
" "<<b1<<endl;
209 double det=a00*a11-a01*a01;
212 cout<<
" linefit2::det="<<det<<endl;
215 if(fabs(det)<0.00000001)
return fit;
216 co1=(b0*a11-b1*a01)/det;
217 co2=(b1*a00-b0*a10)/det;
221 cout<<
" linefit2::Previous element= "<<y[x.size()-3]<<
" "<<x[x.size()-3]<<endl;
224 if(y[x.size()-2]<14.) {
227 cout<<
" check 90 degree track "<<endl;
230 if(fabs(x[x.size()-2]-x[x.size()-1])<0.1){
233 cout<<
" Redetermine line - 90 degree "<<endl;
236 if(fabs(x[x.size()-2]-x[x.size()-1])>0.0001){
237 co1=(y[x.size()-2]-y[x.size()-1])/(x[x.size()-2]-x[x.size()-1]);
238 co2=y[x.size()-2]-co1*x[x.size()-2];
249 cout<<
"linefit2::co1,co2="<<co1<<
" "<<co2<<
" size "<<x.size()<<endl;
253 for (
unsigned int i=0;
i<x.size();
i++){
254 double zdet=(y[
i]-co2)/co1;
255 chi2=chi2+(x[
i]-zdet)*(x[
i]-zdet)/(err[
i]*err[
i]);
258 cout<<
"linefit2::chi2="<<chi2<<
" err="<<err[
i]<<
" x[i]="<<x[
i]<<
" teor="<<zdet<<endl;
265 cout<<
" linefit2::chi2= "<<chi2<<endl;
271 double& co1,
double& chi2)
280 cout<<
" linefit1::sizes="<<x.size()<<
" "<<y.size()<<
" "<<err.size()<<endl;
283 if(x.size() != y.size())
return fit;
284 if(x.size() != err.size())
return fit;
286 for (
unsigned int i=0;
i<x.size();
i++){
287 s1=s1+(x[
i]/err[
i])*(y[
i]/err[
i]);
288 s2=s2+(x[
i]/err[
i])*(x[i]/err[i]);
291 cout<<
"linefit1="<<x[
i]<<
" "<<y[
i]<<
" "<<err[
i]<<
" "<<s1<<
" "<<s2<<endl;
299 cout<<
"linefit1::co1,co2="<<co1<<endl;
303 for (
unsigned int i=0;
i<x.size();
i++){
304 chi2=chi2+(y[
i]-co1*x[
i])*(y[
i]-co1*x[
i])/(err[
i]*err[
i]);
308 cout<<
"linefit1::chi2="<<chi2<<endl;
317 double pi = 4.*atan(1.);
325 double xdouble=xrclus/(2.*rc);
326 psi= phiclus+acharge*asin(xdouble);
330 psi=phiclus+acharge*0.006*fabs(zclus)/fabs(pl);
332 double phic = psi-acharge*pi/2.;
338 cout<<
"radius of track="<<rc<<endl;
339 cout<<
"acharge="<<acharge<<endl;
340 cout<<
"psi="<<psi<<endl;
341 cout<<
"phic="<<phic<<
" pi="<<pi<<
" pi2="<<
pi2<<endl;
348 vector<double>& dphihit, vector<double>& drhit,
349 vector<double>& ehitstrip, vector<double>& dehitphi,
353 double& chirz,
double& chirf,
int& tType)
const{
357 double pi = 4.*atan(1.);
358 double twopi=8.*atan(1.);
362 double ch1,dphi,dr,ptnew;
368 fitrz=this->linefit2(rhit,zhit,ehitstrip,co1,co2,chirz);
370 #ifdef UPDATOR_BARREL_DEBUG
371 cout<<
"UPDATE::BARREL::line fit rz= "<<fitrz<<
" chirz="<<chirz<<endl;
372 cout<<
" co1="<<co1<<
" co2="<<co2<<endl;
377 cout<<
"UPDATE::BARREL::line fit failed rz= "<<fitrz<<
" chirz="<<chirz<<endl;
382 if(dphihit.size()>1){
383 fitrf=this->linefit1(dphihit,drhit,dehitphi,ch1,chirf);
385 #ifdef UPDATOR_BARREL_DEBUG
386 cout<<
"UPDATE::BARREL::line fit dphi= "<<fitrf<<
" chirf="<<chirf<<endl;
387 cout<<
" ch1="<<ch1<<endl;
392 cout<<
"UPDATE::BARREL::line fit failed dphi= "<<fitrf<<
" chirz="<<chirf<<endl;
399 dphi=fabs(dphihit.back());
400 dr=fabs(drhit.back());
401 if(dphi > pi) dphi = twopi-dphi;
404 #ifdef UPDATOR_BARREL_DEBUG
405 cout<<
"UPDATE::BARREL::line calc dphi= "<<dphi<<
" dr="<<dr<<
" chirf="<<chirf<<endl;
406 cout<<
" ch1="<<ch1<<endl;
415 double phiclus=xrhit.
phi();
416 double xrclus=xrhit.
perp();
418 double rc=100.*ptnew/(0.3*4.);
419 double xdouble=xrclus/(2.*rc);
420 double psi=phiclus-aCharge*asin(xdouble);
421 double pznew=ptnew/co1;
422 double znew=(xrclus-co2)/co1;
423 double delphinew=fabs(0.006*drhit.back()/ptnew);
424 double phinew=pRecHit->globalPosition().phi()+aCharge*delphinew;
428 m(0,0)=ptnew;
m(1,1)=thePhiWin,
433 #ifdef UPDATOR_BARREL_DEBUG
434 cout<<
"MuUpdator::xnew=" << xnew<<endl;
446 vector<double>& dphihit, vector<double>& drhit,
447 vector<double>& ehitstrip, vector<double>& dehitphi,
451 double& chirz,
double& chirf,
int& tType)
const{
454 double pi = 4.*atan(1.);
455 double twopi=8.*atan(1.);
466 #ifdef UPDATOR_ENDCAP_DEBUG
467 cout<<
"updateEndcap switched on"<<endl;
472 fitrz=this->linefit2(rhit,zhit,ehitstrip,co1,co2,chirz);
474 #ifdef UPDATOR_ENDCAP_DEBUG
475 cout<<
"line fit rz= "<<fitrz<<
" chirz="<<chirz<<endl;
476 cout<<
" co1="<<co1<<
" co2="<<co2<<endl;
481 cout<<
"UPDATE::ENDCAP::line fit failed rz= "<<fitrz<<
" chirz="<<chirz<<endl;
485 if(dphihit.size()>1){
486 fitrf=this->linefit1(dphihit,drhit,dehitphi,ch1,chirf);
487 if(zhit.front()<0.) ch1=-1.*ch1;
488 #ifdef UPDATOR_ENDCAP_DEBUG
489 cout<<
"MuUpdate::barrel::line fit dphi= "<<fitrf<<
" chirf="<<chirf<<endl;
490 cout<<
" ch1="<<ch1<<endl;
494 cout<<
"UPDATE::ENDCAP::line fit failed dphi= "<<fitrf<<
" chirz="<<chirf<<endl;
500 dphi=fabs(dphihit.back());
501 dr=fabs(drhit.back());
502 if(dphi > pi) dphi = twopi-dphi;
504 if(zhit.front()<0.) ch1=-1.*ch1;
506 #ifdef UPDATOR_ENDCAP_DEBUG
507 cout<<
"MuUpdate::barrel::linedphi=" << dphi <<
" dz= "<< dr <<
" ch1= " <<ch1<<
" pznew= "<<0.006/ch1<<endl;
512 double pznew=0.006/ch1;
518 double phiclus=xrhit.
phi();
520 double xzclus=xrhit.
z();
522 double psi=phiclus-aCharge*0.006*fabs(xzclus-zvert)/fabs(pznew);
524 double ptnew=pznew*co1;
526 double xrclus=co1*xzclus+co2;
528 double delphinew=fabs(0.006*drhit.back()/pznew);
530 double phinew=pRecHit->globalPosition().phi()+aCharge*delphinew;
539 #ifdef UPDATOR_ENDCAP_DEBUG
540 cout<<
"MuUpdator::xnew=" << xnew<<endl;
545 m(0,0)=pznew;
m(1,1)=thePhiWin,
bool linefit1(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &err, double &a, double &chi) const
virtual Location location() const =0
Which part of the detector (barrel, endcap)
const GlobalTrajectoryParameters & parameters() const
FreeTrajectoryState * freeTrajectoryState(bool withErrors=true) const
bool linefit2(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &err, double &a, double &b, double &chi) const
Sin< T >::type sin(const T &t)
Geom::Phi< T > phi() const
Geom::Theta< T > theta() const
ConstRecHitPointer recHit() const
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
TrackCharge charge() const
std::map< std::string, int, std::less< std::string > > psi
Geom::Theta< T > theta() const
DataContainer const & measurements() const
double findPhiInVertex(const FreeTrajectoryState &fts, const double &rc, const GeomDetEnumerators::Location) const
GlobalVector momentum() const
Cos< T >::type cos(const T &t)
Tan< T >::type tan(const T &t)
const DetLayer * layer() const
TrajectoryStateOnSurface updateBarrel(std::vector< double > &rhit, std::vector< double > &zhit, std::vector< double > &dphihit, std::vector< double > &drhit, std::vector< double > &ehitstrip, std::vector< double > &dehitphi, const TransientTrackingRecHit::ConstRecHitPointer &pRecHit, const TransientTrackingRecHit::ConstRecHitPointer &nRecHit, const TrajectoryStateOnSurface &nTsos, double &, double &, int &) const
GlobalPoint position() const
const Surface & surface() const
TrajectoryStateOnSurface update(const Trajectory &mt, const TrajectoryStateOnSurface &, const TrajectoryMeasurement &, const DetLayer *, double &, double &) const
TrackCharge charge() const
TrajectoryStateOnSurface updateEndcap(std::vector< double > &rhit, std::vector< double > &zhit, std::vector< double > &dphihit, std::vector< double > &drhit, std::vector< double > &ehitstrip, std::vector< double > &dehitphi, const TransientTrackingRecHit::ConstRecHitPointer &pRecHit, const TransientTrackingRecHit::ConstRecHitPointer &nRecHit, const TrajectoryStateOnSurface &nTsos, double &, double &, int &) const