30 long int ic, jc, it, jt;
36 long int NVar = (NTRK + 1) * 3;
47 for (it = 1; it <= NTRK; ++it) {
54 for (jt = 1; jt <= NTRK; ++jt) {
55 for (
int k = 0; k < 3; ++k) {
56 for (
int l = 0; l < 3; ++l) {
58 for (
int j = 0; j < 2; ++j) {
59 dpipj[k][l] += vk->
tmpArr[jt-1]->drdp[j][k] * drdpy[j][l];
63 for (
int k = 1; k <= 3; ++k) {
64 for (
int l = 1; l <= 3; ++l) {
65 ader_ref(it * 3 + k, jt * 3 + l) += dpipj[l-1][k-1];
84 for (it = 1; it <= NTRK; ++it) {
86 ader_ref(1, it * 3 + 1) += cnt * th0t.
X * tf0t.
X;
87 ader_ref(2, it * 3 + 1) += cnt * th0t.
Y * tf0t.
X;
88 ader_ref(3, it * 3 + 1) += cnt * th0t.
Z * tf0t.
X;
89 ader_ref(1, it * 3 + 2) += cnt * th0t.
X * tf0t.
Y;
90 ader_ref(2, it * 3 + 2) += cnt * th0t.
Y * tf0t.
Y;
91 ader_ref(3, it * 3 + 2) += cnt * th0t.
Z * tf0t.
Y;
92 ader_ref(1, it * 3 + 3) += cnt * th0t.
X * tf0t.
Z;
93 ader_ref(2, it * 3 + 3) += cnt * th0t.
Y * tf0t.
Z;
94 ader_ref(3, it * 3 + 3) += cnt * th0t.
Z * tf0t.
Z;
102 for (it = 1; it <= NTRK; ++it) {
103 for (jt = it; jt <= NTRK; ++jt) {
106 ader_ref(it*3 + 1, jt*3 + 1) += cnt * tf0ti.
X * tf0tj.
X;
107 ader_ref(it*3 + 2, jt*3 + 1) += cnt * tf0ti.
Y * tf0tj.
X;
108 ader_ref(it*3 + 3, jt*3 + 1) += cnt * tf0ti.
Z * tf0tj.
X;
109 ader_ref(it*3 + 1, jt*3 + 2) += cnt * tf0ti.
X * tf0tj.
Y;
110 ader_ref(it*3 + 2, jt*3 + 2) += cnt * tf0ti.
Y * tf0tj.
Y;
111 ader_ref(it*3 + 3, jt*3 + 2) += cnt * tf0ti.
Z * tf0tj.
Y;
112 ader_ref(it*3 + 1, jt*3 + 3) += cnt * tf0ti.
X * tf0tj.
Z;
113 ader_ref(it*3 + 2, jt*3 + 3) += cnt * tf0ti.
Y * tf0tj.
Z;
114 ader_ref(it*3 + 3, jt*3 + 3) += cnt * tf0ti.
Z * tf0tj.
Z;
120 for (i=1; i<=NVar-1; ++i) {
121 for (j = i+1; j<=NVar; ++j) {
144 for(i=0; i<NVar+1; i++){ ta[i] = tab.data() + i*(NVar+1);}
145 for (i=1; i<=NVar; ++i)
for (j = i; j<=NVar; ++j) ta[i][j] = ta[j][i] =
ader_ref(i,j);
153 for(i=0; i<NVar+1; i++){ tv[i] = tvb.data() + i*(NVar+1); tr[i] = trb.data() + i*(NVar+1);}
155 vkSVDCmp( ta.data(), NVar, NVar, tw.data(), tv.data());
158 for(i=1; i<NVar+1; i++)
if(fabs(tw[i])>tmax)tmax=fabs(tw[i]);
159 for(i=1; i<NVar+1; i++)
if(fabs(tw[i])/tmax < 1.e-18) tw[i]=0.;
160 for(i=1; i<=NVar; i++){
for(j=1; j<=NVar; j++){
161 tr[i][j]=0.;
for(
int k=1; k<=NVar; k++)
if(tw[k]!=0.) tr[i][j] += ta[i][k]*tv[j][k]/tw[k];
164 for (i=1; i<=NVar; ++i)
for (j=1; j<=NVar; ++j)
ader_ref(i,j)=tr[i][j];
170 for( i=1; i<=NVar; i++){
for(j=i; j<=NVar; j++){
171 double mcheck=0.;
for(
int k=1; k<=NVar; k++) mcheck+=ta[i][k]*
ader_ref(k,j);
172 if(i!=j) maxDiff = (maxDiff > std::abs(mcheck)) ? maxDiff : std::abs(mcheck);
173 if(i==j) maxDiff = (maxDiff > std::abs(1.-mcheck)) ? maxDiff : std::abs(1.-mcheck);
176 if(maxDiff>0.1)
return -1;
181 for (i=1; i<=NVar; i++) {
182 for (j=1; j<=NVar; j++) {
197 for (it=1; it<=NTRK; it++) {
198 t_trk=vk->
tmpArr[it-1].get();
200 - vcov[1] * t_trk->
wbci[1]
201 - vcov[3] * t_trk->
wbci[2];
203 - vcov[2] * t_trk->
wbci[1]
204 - vcov[4] * t_trk->
wbci[2];
206 - vcov[4] * t_trk->
wbci[1]
207 - vcov[5] * t_trk->
wbci[2];
209 - vcov[1] * t_trk->
wbci[4]
210 - vcov[3] * t_trk->
wbci[5];
212 - vcov[2] * t_trk->
wbci[4]
213 - vcov[4] * t_trk->
wbci[5];
215 - vcov[4] * t_trk->
wbci[4]
216 - vcov[5] * t_trk->
wbci[5];
218 - vcov[1] * t_trk->
wbci[7]
219 - vcov[3] * t_trk->
wbci[8];
221 - vcov[2] * t_trk->
wbci[7]
222 - vcov[4] * t_trk->
wbci[8];
224 - vcov[4] * t_trk->
wbci[7]
225 - vcov[5] * t_trk->
wbci[8];
238 for (it = 1; it<=NTRK; ++it) {
239 t_trk=vk->
tmpArr[it-1].get();
240 for (jt=1; jt<=NTRK; ++jt) {
272 std::vector<std::vector< Vect3DF> > tf0t;
273 std::vector< Vect3DF > th0t;
274 std::vector< double > taa;
275 std::vector< Vect3DF > tmpVec;
285 tf0t.push_back( tmpVec );
289 double **R =
new double*[totNC];
for(ic=0; ic<totNC; ic++) R[ic]=
new double[NVar];
290 double **RC=
new double*[totNC];
for(ic=0; ic<totNC; ic++)RC[ic]=
new double[NVar];
291 double *RCRt=
new double[totNC*totNC];
292 for(ic=0; ic<totNC; ic++){
296 for(it=1; it<=NTRK; it++){
297 R[ic][it*3+0]=tf0t[ic][it-1].X;
298 R[ic][it*3+1]=tf0t[ic][it-1].Y;
299 R[ic][it*3+2]=tf0t[ic][it-1].Z;
303 for(ic=0; ic<totNC; ic++){
304 for(j=0; j<NVar; j++){ RC[ic][j]=0;
305 for(i=0; i<NVar; i++) RC[ic][j] += R[ic][i]*
ader_ref(i+1,j+1);
309 for(ic=0; ic<totNC; ic++){
310 for(jc=0; jc<totNC; jc++){ RCRt[ic*totNC + jc] =0.;
311 for(i=0; i<NVar; i++) RCRt[ic*totNC + jc] += RC[ic][i]*R[jc][i];
314 dsinv(totNC, RCRt, totNC, &IERR);
315 if ( IERR != 0)
return IERR;
317 for(i=0; i<NVar; i++){
318 for(j=0; j<NVar; j++){
double COR=0.;
319 for(ic=0; ic<totNC; ic++){
320 for(jc=0; jc<totNC; jc++){
321 COR += RC[ic][i]*RC[jc][j]*RCRt[ic*totNC +jc];
328 for(ic=0; ic<totNC; ic++)
delete[] R[ic];
330 for(ic=0; ic<totNC; ic++)
delete[] RC[ic];
346 for (i = 1; i <= 6; ++i) {
347 for (j = 1; j <= 6; ++j) {
349 for (ic=1; ic<=NVar; ++ic) {
350 if(
dcv_ref(i, ic)==0.)
continue;
351 for (jc=1; jc<=NVar; ++jc) {
352 if(
dcv_ref(j, jc)==0.)
continue;