ATLAS Offline Software
Loading...
Searching...
No Matches
Derclc2.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3*/
4
5#include "TrkVKalVrtCore/CommonPars.h"
6#include "TrkVKalVrtCore/cfMomentum.h"
7#include "TrkVKalVrtCore/Derivt.h"
8#include "TrkVKalVrtCore/TrkVKalVrtCoreBase.h"
9#include "TrkVKalVrtCore/VKalVrtBMag.h"
10#include <array>
11#include <algorithm>
12#include <cmath>
13
14namespace Trk {
15//
16// Pointing constraint calculation in new data model.
17//
18// cnstV and cnstP values are used!!!
19//-----------------------------------------------
20void calcPointConstraint( VKPointConstraint * cnst )
21{
22 VKConstraintBase * base_cnst = static_cast<VKConstraintBase*>( cnst);
23 const VKVertex * vk=cnst->getOriginVertex();
24 // Magnetic field in fitted vertex position
25 double cnstPos[3];
26 cnstPos[0]=vk->refIterV[0]+vk->cnstV[0];
27 cnstPos[1]=vk->refIterV[1]+vk->cnstV[1];
28 cnstPos[2]=vk->refIterV[2]+vk->cnstV[2];
29 double localField = Trk::vkalMagFld::getMagFld(cnstPos,(vk->vk_fitterControl).get());
30 double magConst = localField * vkalMagCnvCst;
31 //
32 int it,Charge=0;
33 double ptot[4]={0.,0.,0.,0.};
34 int NTRK = vk->TrackList.size();
35
36 std::vector< std::array<double, 4> > pp(NTRK);
37 for( it=0; it<NTRK; it++){
38 VKTrack * trk = vk->TrackList[it].get();
39 pp[it]=getCnstParticleMom( trk, localField );
40 ptot[0] += pp[it][0];
41 ptot[1] += pp[it][1];
42 ptot[2] += pp[it][2];
43 ptot[3] += pp[it][3];
44 Charge += trk->Charge;
45 }
46 //Summary particle at fitted vertex
47 double Pt= sqrt(ptot[0]*ptot[0] + ptot[1]*ptot[1]) ;
48 double mom= sqrt(ptot[2]*ptot[2] + Pt*Pt) ;
49 double phi=atan2(ptot[1],ptot[0]);
50 double theta=acos( ptot[2]/mom );
51 if(theta<1.e-5) theta=1.e-5;
52 if(theta>M_PI-1.e-5) theta=M_PI-1.e-5;
53 double invR = magConst / Pt; if(Charge) invR *= Charge;
54 //
55 // Fitted vertex position in global reference frame
56 double curV[3] = {vk->refIterV[0]+vk->cnstV[0], vk->refIterV[1]+vk->cnstV[1],vk->refIterV[2]+vk->cnstV[2]};
57 double Px=ptot[0];
58 double Py=ptot[1];
59 double Pz=ptot[2];
60 double sinp = sin(phi);
61 double cosp = cos(phi);
62 double cott = 1./tan(theta);
63 //
64 //Derivatives
65 double dPhiDpx = -Py/Pt/Pt; //dPhi/dPx
66 double dPhiDpy = Px/Pt/Pt; //dPhi/dPy
67 //double dPhiDpz = 0; //dPhi/dPz
68 double dThetaDpx = Px*Pz/(Pt*mom*mom); //dTheta/dPx
69 double dThetaDpy = Py*Pz/(Pt*mom*mom); //dTheta/dPy
70 double dThetaDpz = -Pt/(mom*mom); //dTheta/dPz
71 double dRDpx = Px/(Pt*magConst); //dR/dPx
72 double dRDpy = Py/(Pt*magConst); //dR/dPy
73 //double dRDpz = 0; //dR/dPz
74 if(Charge){dRDpx /=Charge; dRDpy /=Charge;}
75
76// Main part
77
78 double a0,t,z0;
79 double dA0dPhi,dZdPhi,dZdTheta,dA0dR;
80 if( Charge==0 ) {
81 a0 =-sinp * (cnst->getTargetVertex()[0] - curV[0]) + //Perigee with respect to targetVertex
82 cosp * (cnst->getTargetVertex()[1] - curV[1]);
83 t = cosp * (cnst->getTargetVertex()[0] - curV[0]) +
84 sinp * (cnst->getTargetVertex()[1] - curV[1]);
85 z0 = curV[2] + t*cott - cnst->getTargetVertex()[2];
86//----------------------------------------------------------------------
87// First constraint a0=0
88//
89 dA0dPhi = -t;
90 for( it=0; it<NTRK; it++){
91 double invRi = vk->TrackList[it]->cnstP[2];
92 base_cnst->f0t[it][0].X = 0.; // dPx/dTheta==0, dPy/dTheta==0
93 base_cnst->f0t[it][0].Y = dA0dPhi * (dPhiDpx*(-pp[it][1]) + dPhiDpy*( pp[it][0])) ;
94 base_cnst->f0t[it][0].Z = dA0dPhi * (dPhiDpx*(-pp[it][0]) + dPhiDpy*(-pp[it][1]))/invRi ;
95 }
96 base_cnst->h0t[0].X = sinp; // dA0/dFitVertex==dA0/dCurV
97 base_cnst->h0t[0].Y = -cosp;
98 base_cnst->h0t[0].Z = 0.;
99//----------------------------------------------------------------------
100// Second constraint z=0
101//
102 dZdPhi = a0*cott;
103 dZdTheta = -t/(sin(theta)*sin(theta));
104 for( it=0; it<NTRK; it++){
105 double invRi = vk->TrackList[it]->cnstP[2];
106 double p2i = pp[it][0]*pp[it][0]+pp[it][1]*pp[it][1]+pp[it][2]*pp[it][2];
107 double pti = sqrt(pp[it][0]*pp[it][0]+pp[it][1]*pp[it][1]);
108 base_cnst->f0t[it][1].X = 0.; // dPx/dTheta==0, dPy/dTheta==0 // ../dTheta_i
109 base_cnst->f0t[it][1].Y = dZdPhi * (dPhiDpx*(-pp[it][1]) + dPhiDpy*( pp[it][0])) ; // ../dPhi_i
110 base_cnst->f0t[it][1].Z = dZdPhi * (dPhiDpx*(-pp[it][0]) + dPhiDpy*(-pp[it][1]))/invRi; // ../dInvR_i
111 base_cnst->f0t[it][1].X += dZdTheta * ( dThetaDpz*(-p2i/pti)); // ../dTheta_i
112 base_cnst->f0t[it][1].Y += dZdTheta * (dThetaDpx*(-pp[it][1]) + dThetaDpy*( pp[it][0]) + 0.); // ../dPhi_; dPz/dPhi_i=0
113 base_cnst->f0t[it][1].Z += dZdTheta * (dThetaDpx*(-pp[it][0]) + dThetaDpy*(-pp[it][1]) + dThetaDpz*(-pp[it][2]))/invRi;// ../dInvR_i
114 }
115 base_cnst->h0t[1].X = -cosp*cott; // dZ/dCurV
116 base_cnst->h0t[1].Y = -sinp*cott;
117 base_cnst->h0t[1].Z = 1.;
118//std::cout<<" test0="<<Charge<<", "<<a0<<", "<<z0<<", "<<t<<", "<<Px<<", "<<Py<<", "<<Pz<<'\n';
119//
120// The same for charged particle
121//
122
123 }else{
124 int Sign = 1; if(Charge<0)Sign=-1;
125 double R = 1./invR;
126 double xc = -R*sinp + curV[0]; // centre of circle
127 double yc = R*cosp + curV[1];
128 double diffx = cnst->getTargetVertex()[0] - xc;
129 double diffy = cnst->getTargetVertex()[1] - yc;
130 double diff = sqrt( diffx*diffx + diffy*diffy );
131 double sindphi = (curV[0]-xc)*diffy - (curV[1]-yc)*diffx;
132 sindphi = sindphi/std::abs(R)/diff;
133 double dphi = asin( sindphi );
134 double cosdphi=sqrt(std::abs(1.-sindphi*sindphi)); if(cosdphi<1.e-10)cosdphi=1.e-10;
135 a0 = R*(1.-diff/std::abs(R));
136 t = dphi*R;
137 z0 = curV[2] + t*cott - cnst->getTargetVertex()[2];
138 double dTdXcur = std::abs(R)/cosdphi * (cosp/diff-diffx*sindphi/(diff*diff)) * (-1.);
139 double dTdYcur = std::abs(R)/cosdphi * (sinp/diff-diffy*sindphi/(diff*diff)) * (-1.);
140//----------------------------------------------------------------------
141// First constraint a0=0
142//
143 dA0dPhi = - (diffx*cosp + diffy*sinp)/diff*std::abs(R);
144 dA0dR = 1. - (diffx*sinp - diffy*cosp)/diff*R/std::abs(R);
145 for( it=0; it<NTRK; it++){
146 double invRi = vk->TrackList[it]->cnstP[2];
147 base_cnst->f0t[it][0].X = 0.; // dPx/dTheta==0, dPy/dTheta==0
148 base_cnst->f0t[it][0].Y = dA0dPhi * (dPhiDpx*(-pp[it][1]) + dPhiDpy*( pp[it][0])) ;
149 base_cnst->f0t[it][0].Z = dA0dPhi * (dPhiDpx*(-pp[it][0]) + dPhiDpy*(-pp[it][1]))/invRi ;
150 base_cnst->f0t[it][0].X += 0.;
151 base_cnst->f0t[it][0].Y += dA0dR * (dRDpx*(-pp[it][1]) + dRDpy*( pp[it][0]));
152 base_cnst->f0t[it][0].Z += dA0dR * (dRDpx*(-pp[it][0]) + dRDpy*(-pp[it][1]))/invRi ;
153 }
154 base_cnst->h0t[0].X = Sign*diffx/diff; // dA0/dFitVertex==dA0/dCurV
155 base_cnst->h0t[0].Y = Sign*diffy/diff;
156 base_cnst->h0t[0].Z = 0.;
157//----------------------------------------------------------------------
158// Second constraint z=0
159//
160 dZdPhi = -a0*cott;
161 dZdTheta = -t/sin(theta)/sin(theta);
162 for( it=0; it<NTRK; it++){
163 double invRi = vk->TrackList[it]->cnstP[2];
164 double p2i = pp[it][0]*pp[it][0]+pp[it][1]*pp[it][1]+pp[it][2]*pp[it][2];
165 double pti = sqrt(pp[it][0]*pp[it][0]+pp[it][1]*pp[it][1]);
166 base_cnst->f0t[it][1].X = 0.; // dPx/dTheta==0, dPy/dTheta==0 // ../dTheta_i
167 base_cnst->f0t[it][1].Y = dZdPhi * (dPhiDpx*(-pp[it][1]) + dPhiDpy*( pp[it][0])) ; // ../dPhi_i
168 base_cnst->f0t[it][1].Z = dZdPhi * (dPhiDpx*(-pp[it][0]) + dPhiDpy*(-pp[it][1]))/invRi; // ../dInvR_i
169 base_cnst->f0t[it][1].X += dZdTheta * ( dThetaDpz*(-p2i/pti)); // ../dTheta_i
170 base_cnst->f0t[it][1].Y += dZdTheta * (dThetaDpx*(-pp[it][1]) + dThetaDpy*( pp[it][0]) + 0.); // ../dPhi_; dPz/dPhi_i=0
171 base_cnst->f0t[it][1].Z += dZdTheta * (dThetaDpx*(-pp[it][0]) + dThetaDpy*(-pp[it][1]) + dThetaDpz*(-pp[it][2]))/invRi;// ../dInvR_i
172 }
173 base_cnst->h0t[1].X = dTdXcur*cott; // dZ/dCurV
174 base_cnst->h0t[1].Y = dTdYcur*cott;
175 base_cnst->h0t[1].Z = 1.;
176//std::cout<<" test1="<<Charge<<", "<<a0<<", "<<z0<<" dst="<<t<<" NTRK="<<NTRK<<'\n';
177 }
178 base_cnst->aa[0]=a0;
179 base_cnst->aa[1]=z0;
180
181
182// -----------------------Check of propagation
183// double paro[5],parn[5],s,ref[3],pere[3]; long int ichg=0;
184// paro[0]=0.; paro[1]=0.; paro[2]=theta; paro[3]=phi; paro[4]=invR;
185// ref[0]=cnst->getTargetVertex()[0]-curV[0]; ref[1]=cnst->getTargetVertex()[1]-curV[1]; ref[2]=cnst->getTargetVertex()[2]-curV[2];
186// cfnewp(&ichg, paro, ref, &s, parn, pere);
187// std::cout<<" testp="<<Charge<<", "<<parn[0]<<", "<<parn[1]<<", "<<s<<'\n';
188// -----------------------Derivatives
189//std::cout<<base_cnst->f0t[0][0].X<<", "<<base_cnst->f0t[0][0].Y<<", "<<base_cnst->f0t[0][0].Z<<'\n';
190//std::cout<<base_cnst->f0t[0][1].X<<", "<<base_cnst->f0t[0][1].Y<<", "<<base_cnst->f0t[0][1].Z<<'\n';
191//std::cout<<base_cnst->f0t[1][0].X<<", "<<base_cnst->f0t[1][0].Y<<", "<<base_cnst->f0t[1][0].Z<<'\n';
192//std::cout<<base_cnst->f0t[1][1].X<<", "<<base_cnst->f0t[1][1].Y<<", "<<base_cnst->f0t[1][1].Z<<'\n';
193//std::cout<<base_cnst->h0t[0].X<<", "<<base_cnst->h0t[0].Y<<", "<<base_cnst->h0t[0].Z<<'\n';
194//std::cout<<base_cnst->h0t[1].X<<", "<<base_cnst->h0t[1].Y<<", "<<base_cnst->h0t[1].Z<<'\n';
195//std::cout<<"-------"<<'\n';
196//--------------------------------------------------------------------------------------------------
197 if( cnst->onlyZ() ) { //Z only constraint <= all Rphi ralated stuff is 0
198 for( it=0; it<NTRK; it++){
199 base_cnst->f0t[it][0].X = 0.001;
200 base_cnst->f0t[it][0].Y = 0.001;
201 base_cnst->f0t[it][0].Z = 0.001;
202 }
203 base_cnst->h0t[0].X=0.; base_cnst->h0t[0].Y=0.; base_cnst->h0t[0].Z=0.;
204 base_cnst->aa[0]=0.;
205 }
206}
207
208} /* End of namespace */
209
M_PI
#define M_PI
Definition ActiveFraction.h:11
vkalMagCnvCst
#define vkalMagCnvCst
Definition CommonPars.h:23
Sign
int Sign(int in)
Definition JSSTaggerUtils.cxx:23
diff
void diff(const Jet &rJet1, const Jet &rJet2, std::map< std::string, double > varDiff)
Difference between jets - Non-Class function required by trigger.
Definition Jet.cxx:631
Trk::VKConstraintBase::f0t
std::vector< std::vector< Vect3DF > > f0t
Definition Derivt.h:39
Trk::VKConstraintBase::getOriginVertex
const VKVertex * getOriginVertex() const
Definition Derivt.h:32
Trk::VKConstraintBase::h0t
std::vector< Vect3DF > h0t
Definition Derivt.h:40
Trk::VKConstraintBase::aa
std::vector< double > aa
Definition Derivt.h:38
Trk::VKPointConstraint::onlyZ
bool onlyZ() const
Definition Derivt.h:94
Trk::VKPointConstraint::getTargetVertex
const double * getTargetVertex() const
Definition Derivt.h:100
Trk::VKVertex::TrackList
std::vector< std::unique_ptr< VKTrack > > TrackList
Definition TrkVKalVrtCoreBase.h:167
Trk::VKVertex::vk_fitterControl
std::unique_ptr< VKalVrtControl > vk_fitterControl
Definition TrkVKalVrtCoreBase.h:194
Trk::vkalMagFld::getMagFld
static void getMagFld(const double, const double, const double, double &, double &, double &, const VKalVrtControlBase *)
Definition VKalVrtBMag.cxx:24
a0
double a0
Definition globals.cxx:27
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition FakeTrackBuilder.h:9
Trk::calcPointConstraint
void calcPointConstraint(VKPointConstraint *cnst)
Definition Derclc2.cxx:20
Trk::theta
@ theta
Definition ParamDefs.h:66
Trk::phi
@ phi
Definition ParamDefs.h:75
Trk::z0
@ z0
Definition ParamDefs.h:64
Trk::getCnstParticleMom
std::array< double, 4 > getCnstParticleMom(const VKTrack *trk, const VKVertex *vk)
Definition cfMomentum.cxx:96
Generated on for ATLAS Offline Software by doxygen 1.14.0