ATLAS Offline Software
Derclc2.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 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 
14 namespace Trk {
15 //
16 // Pointing constraint calculation in new data model.
17 //
18 // cnstV and cnstP values are used!!!
19 //-----------------------------------------------
20 void calcPointConstraint( VKPointConstraint * cnst )
21 {
22  VKConstraintBase * base_cnst = (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 
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
CommonPars.h
Trk::vkalMagFld::getMagFld
static void getMagFld(const double, const double, const double, double &, double &, double &, VKalVrtControlBase *)
Definition: VKalVrtBMag.cxx:24
Trk::VKPointConstraint::onlyZ
bool onlyZ() const
Definition: Derivt.h:92
skel.it
it
Definition: skel.GENtoEVGEN.py:396
M_PI
#define M_PI
Definition: ActiveFraction.h:11
mc.diff
diff
Definition: mc.SFGenPy8_MuMu_DD.py:14
Trk::z0
@ z0
Definition: ParamDefs.h:64
Trk::VKVertex::refIterV
double refIterV[3]
Definition: TrkVKalVrtCoreBase.h:146
Trk::VKTrack
Definition: TrkVKalVrtCoreBase.h:64
Trk::VKTrack::Charge
int Charge
Definition: TrkVKalVrtCoreBase.h:73
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
Trk::getCnstParticleMom
std::array< double, 4 > getCnstParticleMom(const VKTrack *trk, const VKVertex *vk)
Definition: cfMomentum.cxx:92
VKalVrtBMag.h
Trk::VKConstraintBase
Definition: Derivt.h:24
TrkVKalVrtCoreBase.h
Trk::VKVertex::TrackList
std::vector< std::unique_ptr< VKTrack > > TrackList
Definition: TrkVKalVrtCoreBase.h:167
ParticleGun_EoverP_Config.mom
mom
Definition: ParticleGun_EoverP_Config.py:63
Trk::VKPointConstraint
Definition: Derivt.h:87
vkalMagCnvCst
#define vkalMagCnvCst
Definition: CommonPars.h:23
Trk::theta
@ theta
Definition: ParamDefs.h:66
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
Trk::calcPointConstraint
void calcPointConstraint(VKPointConstraint *cnst)
Definition: Derclc2.cxx:20
a0
double a0
Definition: globals.cxx:27
Trk::VKConstraintBase::f0t
std::vector< std::vector< Vect3DF > > f0t
Definition: Derivt.h:37
Derivt.h
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
Trk::VKVertex::vk_fitterControl
std::unique_ptr< VKalVrtControl > vk_fitterControl
Definition: TrkVKalVrtCoreBase.h:194
Trk::VKConstraintBase::aa
std::vector< double > aa
Definition: Derivt.h:36
Trk::VKVertex::cnstV
double cnstV[3]
Definition: TrkVKalVrtCoreBase.h:143
Prompt::Def::Pt
@ Pt
Definition: VarHolder.h:76
Trk::VKConstraintBase::getOriginVertex
const VKVertex * getOriginVertex() const
Definition: Derivt.h:30
Trk::VKConstraintBase::h0t
std::vector< Vect3DF > h0t
Definition: Derivt.h:38
Trk::VKPointConstraint::getTargetVertex
const double * getTargetVertex() const
Definition: Derivt.h:98
Trk::VKVertex
Definition: TrkVKalVrtCoreBase.h:128
Trk::phi
@ phi
Definition: ParamDefs.h:75
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
cfMomentum.h
Generated on Thu Nov 7 2024 21:13:36 for ATLAS Offline Software by   doxygen 1.8.18