ATLAS Offline Software
FullMtx.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/TrkVKalVrtCoreBase.h"
7 #include "TrkVKalVrtCore/Derivt.h"
8 
9 namespace Trk {
10 
11 /* --------------------------------------------------- */
12 /* ADER is a full matrix of chi2 problem */
13 /* It inversion gives the same solution as usual */
14 /* Author: V.Kostyukhin */
15 /* --------------------------------------------------- */
16 
17 #define ader_ref(a_1,a_2) ader[(a_2)*(vkalNTrkM*3+3) + (a_1)]
18 
19 void FullMTXfill(VKVertex * vk, double * ader)
20 {
21  int i,j,it;
22 
23  int NTRK = vk->TrackList.size();
24  int NPar=3*NTRK+3;
25 
26  for (i=0; i<NPar; i++) { for (j=0; j<NPar; j++) { ader_ref(i, j)=0.;} }
27 
28  ader_ref(0, 0) += vk->wa[0];
29  ader_ref(0, 1) += vk->wa[1];
30  ader_ref(1, 1) += vk->wa[2];
31  ader_ref(0, 2) += vk->wa[3];
32  ader_ref(1, 2) += vk->wa[4];
33  ader_ref(2, 2) += vk->wa[5];
34  ader_ref(1, 0) = ader_ref(0, 1);
35  ader_ref(2, 0) = ader_ref(0, 2);
36  ader_ref(2, 1) = ader_ref(1, 2);
37 
38 
39  for (it=0; it<NTRK; ++it) {
40  ader_ref(0, it*3 + 3) += vk->tmpArr[it]->wb[0];
41  ader_ref(1, it*3 + 3) += vk->tmpArr[it]->wb[1];
42  ader_ref(2, it*3 + 3) += vk->tmpArr[it]->wb[2];
43  ader_ref(0, it*3 + 4) += vk->tmpArr[it]->wb[3];
44  ader_ref(1, it*3 + 4) += vk->tmpArr[it]->wb[4];
45  ader_ref(2, it*3 + 4) += vk->tmpArr[it]->wb[5];
46  ader_ref(0, it*3 + 5) += vk->tmpArr[it]->wb[6];
47  ader_ref(1, it*3 + 5) += vk->tmpArr[it]->wb[7];
48  ader_ref(2, it*3 + 5) += vk->tmpArr[it]->wb[8];
49  }
50 
51  for (it=0; it<NTRK; ++it) {
52  ader_ref(it*3 + 3, it*3 + 3) += vk->tmpArr[it]->wc[0];
53  ader_ref(it*3 + 3, it*3 + 4) += vk->tmpArr[it]->wc[1];
54  ader_ref(it*3 + 4, it*3 + 4) += vk->tmpArr[it]->wc[2];
55  ader_ref(it*3 + 3, it*3 + 5) += vk->tmpArr[it]->wc[3];
56  ader_ref(it*3 + 4, it*3 + 5) += vk->tmpArr[it]->wc[4];
57  ader_ref(it*3 + 5, it*3 + 5) += vk->tmpArr[it]->wc[5];
58  ader_ref(it*3 + 4, it*3 + 3) += vk->tmpArr[it]->wc[1];
59  ader_ref(it*3 + 5, it*3 + 3) += vk->tmpArr[it]->wc[3];
60  ader_ref(it*3 + 5, it*3 + 4) += vk->tmpArr[it]->wc[4];
61  }
62 /* symmetrisation */
63  for (i=0; i<NPar-1; i++) {
64  for (j=i+1; j<NPar; ++j) {
65  ader_ref(j, i) = ader_ref(i, j);
66  }
67  }
68 //std::cout<<" FULLMTX NEW="<<ader_ref(0, 0)<<", "<<ader_ref(1, 1)<<", "<<ader_ref(2, 2)<<", "
69 // <<ader_ref(3, 3)<<", "<<ader_ref(4, 4)<<", "<<ader_ref(5, 5)<<", "
70 // <<ader_ref(0, 3)<<", "<<ader_ref(1, 4)<<", "<<ader_ref(3, 5)<<'\n';
71 }
72 
73 #undef ader_ref
74 
75 //---------------------------------------------------------------------------------------------
76 // Fill full matrix of vertex with constraints problem
77 // It's used for cascade fit
78 //
79 // Total number of parameters Vertex(3)+NTRK*3+N_constr
80 
81 int FullMCNSTfill(VKVertex * vk, double * ader, double * LSide)
82 {
83  int i,j,k,l,ii,ic,it,jt;
84 
85  int totNC=0; //total number of constraints
86  for(i=0; i<(int)vk->ConstraintList.size();i++)totNC += vk->ConstraintList[i]->NCDim;
87 
88  int NTRK = vk->TrackList.size();
89  int NPar = 3*NTRK+3+totNC;
90 // int NPar=3*NTrkM+3;
91 
92  for (i=0; i<NPar; i++) { for (j=0; j<NPar; j++) ader[i+j*NPar]=0.; }
93  std::vector<std::vector< const Vect3DF*> > tf0t; // derivative collectors
94  std::vector< const Vect3DF* > th0t; // derivative collectors
95  std::vector< double > taa; // derivative collectors
96  std::vector< const Vect3DF* > tmpVec;
97  for(ii=0; ii<(int)vk->ConstraintList.size();ii++){
98  for(ic=0; ic<(int)vk->ConstraintList[ii]->NCDim; ic++){
99  taa.push_back( vk->ConstraintList[ii]->aa[ic] );
100  th0t.push_back( &(vk->ConstraintList[ii]->h0t[ic]) );
101  tmpVec.clear();
102  tmpVec.reserve(vk->ConstraintList[ii]->f0t.size());
103  for(it=0; it<(int)vk->ConstraintList[ii]->f0t.size(); it++){
104  tmpVec.push_back( &(vk->ConstraintList[ii]->f0t[it][ic]) );
105  }
106  tf0t.push_back( std::move(tmpVec) );
107  }
108  }
109 //-----------------------------------------------------------------------------
110  ader[0*NPar+ 0] += vk->wa[0];
111  ader[0*NPar+ 1] += vk->wa[1];
112  ader[1*NPar+ 1] += vk->wa[2];
113  ader[0*NPar+ 2] += vk->wa[3];
114  ader[1*NPar+ 2] += vk->wa[4];
115  ader[2*NPar+ 2] += vk->wa[5];
116  ader[1*NPar+ 0] = ader[0*NPar+ 1];
117  ader[2*NPar+ 0] = ader[0*NPar+ 2];
118  ader[2*NPar+ 1] = ader[1*NPar+ 2];
119 
120 
121  for (it=0; it<NTRK; ++it) {
122  ader[0 + (it*3 + 3)*NPar] += vk->tmpArr[it]->wb[0];
123  ader[1 + (it*3 + 3)*NPar] += vk->tmpArr[it]->wb[1];
124  ader[2 + (it*3 + 3)*NPar] += vk->tmpArr[it]->wb[2];
125  ader[0 + (it*3 + 4)*NPar] += vk->tmpArr[it]->wb[3];
126  ader[1 + (it*3 + 4)*NPar] += vk->tmpArr[it]->wb[4];
127  ader[2 + (it*3 + 4)*NPar] += vk->tmpArr[it]->wb[5];
128  ader[0 + (it*3 + 5)*NPar] += vk->tmpArr[it]->wb[6];
129  ader[1 + (it*3 + 5)*NPar] += vk->tmpArr[it]->wb[7];
130  ader[2 + (it*3 + 5)*NPar] += vk->tmpArr[it]->wb[8];
131  }
132 
133  for (it=0; it<NTRK; ++it) {
134  ader[(it*3 + 3) + (it*3 + 3)*NPar] += vk->tmpArr[it]->wc[0];
135  ader[(it*3 + 3) + (it*3 + 4)*NPar] += vk->tmpArr[it]->wc[1];
136  ader[(it*3 + 4) + (it*3 + 4)*NPar] += vk->tmpArr[it]->wc[2];
137  ader[(it*3 + 3) + (it*3 + 5)*NPar] += vk->tmpArr[it]->wc[3];
138  ader[(it*3 + 4) + (it*3 + 5)*NPar] += vk->tmpArr[it]->wc[4];
139  ader[(it*3 + 5) + (it*3 + 5)*NPar] += vk->tmpArr[it]->wc[5];
140  ader[(it*3 + 4) + (it*3 + 3)*NPar] += vk->tmpArr[it]->wc[1];
141  ader[(it*3 + 5) + (it*3 + 3)*NPar] += vk->tmpArr[it]->wc[3];
142  ader[(it*3 + 5) + (it*3 + 4)*NPar] += vk->tmpArr[it]->wc[4];
143  }
144 /* symmetrisation */
145  for (i=0; i<NPar-1; i++) {
146  for (j=i+1; j<NPar; ++j) {
147  ader[j + i*NPar] = ader[i + j*NPar];
148  }
149  }
150 //
151 // For "passing near" constraint
152  if (vk->passNearVertex){
153  double drdpy[2][3],dpipj[3][3];
154  for (it = 0; it < NTRK; ++it) {
155  drdpy[0][0] = vk->tmpArr[it]->drdp[0][0] * vk->FVC.ywgt[0] + vk->tmpArr[it]->drdp[1][0] * vk->FVC.ywgt[1];
156  drdpy[1][0] = vk->tmpArr[it]->drdp[0][0] * vk->FVC.ywgt[1] + vk->tmpArr[it]->drdp[1][0] * vk->FVC.ywgt[2];
157  drdpy[0][1] = vk->tmpArr[it]->drdp[0][1] * vk->FVC.ywgt[0] + vk->tmpArr[it]->drdp[1][1] * vk->FVC.ywgt[1];
158  drdpy[1][1] = vk->tmpArr[it]->drdp[0][1] * vk->FVC.ywgt[1] + vk->tmpArr[it]->drdp[1][1] * vk->FVC.ywgt[2];
159  drdpy[0][2] = vk->tmpArr[it]->drdp[0][2] * vk->FVC.ywgt[0] + vk->tmpArr[it]->drdp[1][2] * vk->FVC.ywgt[1];
160  drdpy[1][2] = vk->tmpArr[it]->drdp[0][2] * vk->FVC.ywgt[1] + vk->tmpArr[it]->drdp[1][2] * vk->FVC.ywgt[2];
161  for (jt = 0; jt < NTRK; ++jt) { /* Matrix */
162  for (k=0; k<3; ++k) {
163  for (l=0; l<3; ++l) {
164  dpipj[k][l] = 0.;
165  for (j=0; j<2; ++j) {
166  dpipj[k][l] += vk->tmpArr[jt]->drdp[j][k] * drdpy[j][l];
167  }
168  }
169  }
170  for (k=0; k<3; ++k) {
171  for (l=0; l<3; ++l) {
172  ader[(it*3+3 +k) + (jt*3+3 +l)*NPar] += dpipj[l][k];
173  }
174  }
175  } // jt cycle
176  } // it cycle
177  }
178 //
179 // Part for lagrange multipliers
180 //
181  int NTrP=NTRK*3 + 3; // track part of matrix
182  for(ic=0; ic<totNC; ic++){
183  ader[(0) + (NTrP+ic)*NPar] = -2.*th0t[ic]->X;
184  ader[(1) + (NTrP+ic)*NPar] = -2.*th0t[ic]->Y;
185  ader[(2) + (NTrP+ic)*NPar] = -2.*th0t[ic]->Z;
186  ader[(0)*NPar + (NTrP+ic) ] = -2.*th0t[ic]->X;
187  ader[(1)*NPar + (NTrP+ic) ] = -2.*th0t[ic]->Y;
188  ader[(2)*NPar + (NTrP+ic) ] = -2.*th0t[ic]->Z;
189  for (it=0; it<NTRK; ++it) {
190  ader[(it*3+3+0) + (NTrP+ic)*NPar] = - 2.*tf0t[ic][it]->X;
191  ader[(it*3+3+1) + (NTrP+ic)*NPar] = - 2.*tf0t[ic][it]->Y;
192  ader[(it*3+3+2) + (NTrP+ic)*NPar] = - 2.*tf0t[ic][it]->Z;
193  ader[(it*3+3+0)*NPar + (NTrP+ic)] = - 2.*tf0t[ic][it]->X;
194  ader[(it*3+3+1)*NPar + (NTrP+ic)] = - 2.*tf0t[ic][it]->Y;
195  ader[(it*3+3+2)*NPar + (NTrP+ic)] = - 2.*tf0t[ic][it]->Z;
196  }
197  }
198 //
199 // Left part
200 //
201 // double *LSide = new double[NPar];
202  LSide[0]=vk->T[0];
203  LSide[1]=vk->T[1];
204  LSide[2]=vk->T[2];
205  for (it=0; it<NTRK; ++it) {
206  LSide[it*3+3+0]=vk->tmpArr[it]->tt[0];
207  LSide[it*3+3+1]=vk->tmpArr[it]->tt[1];
208  LSide[it*3+3+2]=vk->tmpArr[it]->tt[2];
209  }
210  for(ic=0; ic<totNC; ic++){
211  LSide[NTRK*3+3+ic]=taa[ic];
212  }
213 //--------------------------------------------------
214  //int IERR = vkMSolve(ader, LSide, NPar);
215  //std::cout<<"global="<<NPar<<" err="<<IERR<<'\n';
216  //std::cout<<" Vrt="<<LSide[0]<<", "<<LSide[1]<<", "<<LSide[2]<<'\n';
217  //for(it=0; it<NTRK; it++)std::cout<<" trk="<<LSide[it*3+3+0]<<", "
218  // <<LSide[it*3+3+1]<<", "<<LSide[it*3+3+2]<<'\n';
219  return NPar;
220 }
221 
222 
223 } /* End of namespace */
224 
CommonPars.h
Trk::FullMTXfill
void FullMTXfill(VKVertex *vk, double *ader)
Definition: FullMtx.cxx:19
Trk::VKVertex::tmpArr
std::vector< std::unique_ptr< TWRK > > tmpArr
Definition: TrkVKalVrtCoreBase.h:168
ader_ref
#define ader_ref(a_1, a_2)
Definition: FullMtx.cxx:17
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
Trk::VKVertex::T
double T[3]
Definition: TrkVKalVrtCoreBase.h:163
skel.it
it
Definition: skel.GENtoEVGEN.py:423
Trk::VKVertex::FVC
ForVrtClose FVC
Definition: TrkVKalVrtCoreBase.h:160
UploadAMITag.l
list l
Definition: UploadAMITag.larcaf.py:158
Trk::VKVertex::ConstraintList
std::vector< std::unique_ptr< VKConstraintBase > > ConstraintList
Definition: TrkVKalVrtCoreBase.h:169
TrkVKalVrtCoreBase.h
Trk::VKVertex::TrackList
std::vector< std::unique_ptr< VKTrack > > TrackList
Definition: TrkVKalVrtCoreBase.h:167
lumiFormat.i
int i
Definition: lumiFormat.py:92
Trk::VKVertex::passNearVertex
bool passNearVertex
Definition: TrkVKalVrtCoreBase.h:156
grepfile.ic
int ic
Definition: grepfile.py:33
Derivt.h
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
Trk::FullMCNSTfill
int FullMCNSTfill(VKVertex *vk, double *ader, double *LSide)
Definition: FullMtx.cxx:81
Trk::ForVrtClose::ywgt
double ywgt[3]
Definition: ForVrtClose.h:29
Trk::VKVertex::wa
double wa[6]
Definition: TrkVKalVrtCoreBase.h:164
Trk::VKVertex
Definition: TrkVKalVrtCoreBase.h:128
fitman.k
k
Definition: fitman.py:528
Generated on Fri Jul 5 2024 21:16:10 for ATLAS Offline Software by   doxygen 1.8.18