ATLAS Offline Software
VtCFitC.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/VtCFitC.h"
6 #include "TrkVKalVrtCore/CommonPars.h"
7 #include "TrkVKalVrtCore/Matrix.h"
8 #include "TrkVKalVrtCore/Derivt.h"
9 #include "TrkVKalVrtCore/ForCFT.h"
10 #include "TrkVKalVrtCore/TrkVKalVrtCoreBase.h"
11 #include <cmath>
12 #include <iostream>
13 
14 namespace Trk {
15 
16 int vtcfitc( VKVertex * vk )
17 {
18 
19  int ii,ic,it,jc,index;
20 
21  double tmp[3];
22 
23 /* -------------------------------------------------------- */
24 /* Start of constraint treatment */
25 /* -------------------------------------------------------- */
26  long int IERR = 0;
27  long int totNC=0; //total number of constraints
28  long int NTRK = vk->TrackList.size();
29  if (vk->ConstraintList.empty()) {return 0;}
30 //
31  double dxyz[3]={vk->dxyz0[0],vk->dxyz0[1],vk->dxyz0[2]}; //nonconstraint vertex shift
32 //
33 // Extruction of derivatives
34 //
35  std::vector<std::vector< const Vect3DF*> > tf0t; // derivative collectors
36  std::vector< const Vect3DF* > th0t; // derivative collectors
37  std::vector< double > taa; // derivative collectors
38  th0t.reserve(vk->ConstraintList.size() * vk->ConstraintList[0]->NCDim);
39  tf0t.reserve(vk->ConstraintList.size());
40 
41  for(ii=0; ii<(int)vk->ConstraintList.size();ii++){
42  totNC += vk->ConstraintList[ii]->NCDim;
43  for(ic=0; ic<(int)vk->ConstraintList[ii]->NCDim; ic++){
44  taa.push_back( vk->ConstraintList[ii]->aa[ic] );
45  th0t.push_back( &(vk->ConstraintList[ii]->h0t[ic]) );
46  std::vector< const Vect3DF* > tmpVec;
47  tmpVec.reserve(vk->ConstraintList[ii]->f0t.size());
48  for(it=0; it<(int)vk->ConstraintList[ii]->f0t.size(); it++){
49  tmpVec.push_back( &(vk->ConstraintList[ii]->f0t[it][ic]) );
50  }
51  tf0t.push_back( std::move(tmpVec) );
52  }
53  }
54  if(totNC==0)return 0;
55 //
56  std::vector< std::vector<double> > denom;
57  denom.reserve(totNC);
58  for(int ic=0; ic<totNC; ic++) denom.emplace_back( totNC, 0 );
59 //
60 // Start of calc
61 //
62  //This is deliberately written without make_unqiue to bypass auto intialization!!
63  std::unique_ptr< Vect3DF[] > al0( new Vect3DF[ totNC ]);
64  for(ic=0;ic<totNC; ic++){ al0[ic].X=0.; al0[ic].Y=0.; al0[ic].Z=0.;}
65  std::vector<double> anum(totNC,0.);
66  for (ic=0; ic<totNC; ++ic) {
67  for (it = 0; it<NTRK; ++it) {
68 
69 /* summation of WBCI * F0T into vector AL0 */
70  al0[ic].X += vk->tmpArr[it]->wbci[0] * tf0t[ic][it]->X
71  + vk->tmpArr[it]->wbci[3] * tf0t[ic][it]->Y
72  + vk->tmpArr[it]->wbci[6] * tf0t[ic][it]->Z;
73  al0[ic].Y += vk->tmpArr[it]->wbci[1] * tf0t[ic][it]->X
74  + vk->tmpArr[it]->wbci[4] * tf0t[ic][it]->Y
75  + vk->tmpArr[it]->wbci[7] * tf0t[ic][it]->Z;
76  al0[ic].Z += vk->tmpArr[it]->wbci[2] * tf0t[ic][it]->X
77  + vk->tmpArr[it]->wbci[5] * tf0t[ic][it]->Y
78  + vk->tmpArr[it]->wbci[8] * tf0t[ic][it]->Z;
79 
80  anum[ic] += vk->tmpArr[it]->parf0[0] * tf0t[ic][it]->X
81  + vk->tmpArr[it]->parf0[1] * tf0t[ic][it]->Y
82  + vk->tmpArr[it]->parf0[2] * tf0t[ic][it]->Z;
83 
84  }
85  al0[ic].X -= th0t[ic]->X;
86  al0[ic].Y -= th0t[ic]->Y;
87  al0[ic].Z -= th0t[ic]->Z;
88  anum[ic] = 2.*anum[ic] + dxyz[0] * 2. * th0t[ic]->X
89  + dxyz[1] * 2. * th0t[ic]->Y
90  + dxyz[2] * 2. * th0t[ic]->Z
91  + taa[ic];
92  }
93 
94 //std::cout<<" vertex="<<dxyz[0]<<", "<<dxyz[1]<<", "<<dxyz[2]<<'\n';
95 //std::cout<<" h0t="<<th0t[0].X<<", "<<th0t[0].Y<<", "<<th0t[0].Z<<'\n';
96 //std::cout<<" f0t="<<tf0t[0][0].X<<", "<<tf0t[0][0].Y<<", "<<tf0t[0][0].Z<<'\n';
97 //std::cout<<" f0t="<<tf0t[0][1].X<<", "<<tf0t[0][1].Y<<", "<<tf0t[0][1].Z<<'\n';
98 
99 
100 
101 /* calculation (AL0)t * VCOV * AL0 and (F0T)t * WCI * F0T */
102  for (ic=0; ic<totNC; ++ic) {
103  for (jc=0; jc<totNC; ++jc) {
104  denom[ic][jc] = al0[ic].X * vk->fitVcov[0] * al0[jc].X
105  + al0[ic].Y * vk->fitVcov[1] * al0[jc].X
106  + al0[ic].Z * vk->fitVcov[3] * al0[jc].X
107  + al0[ic].X * vk->fitVcov[1] * al0[jc].Y
108  + al0[ic].Y * vk->fitVcov[2] * al0[jc].Y
109  + al0[ic].Z * vk->fitVcov[4] * al0[jc].Y
110  + al0[ic].X * vk->fitVcov[3] * al0[jc].Z
111  + al0[ic].Y * vk->fitVcov[4] * al0[jc].Z
112  + al0[ic].Z * vk->fitVcov[5] * al0[jc].Z
113  + al0[jc].X * vk->fitVcov[0] * al0[ic].X
114  + al0[jc].Y * vk->fitVcov[1] * al0[ic].X
115  + al0[jc].Z * vk->fitVcov[3] * al0[ic].X
116  + al0[jc].X * vk->fitVcov[1] * al0[ic].Y
117  + al0[jc].Y * vk->fitVcov[2] * al0[ic].Y
118  + al0[jc].Z * vk->fitVcov[4] * al0[ic].Y
119  + al0[jc].X * vk->fitVcov[3] * al0[ic].Z
120  + al0[jc].Y * vk->fitVcov[4] * al0[ic].Z
121  + al0[jc].Z * vk->fitVcov[5] * al0[ic].Z;
122 
123  for (it=0; it<NTRK; ++it) {
124  TWRK * t_trk=vk->tmpArr[it].get();
125  denom[ic][jc] += tf0t[ic][it]->X * t_trk->wci[0] * tf0t[jc][it]->X
126  + tf0t[ic][it]->Y * t_trk->wci[1] * tf0t[jc][it]->X
127  + tf0t[ic][it]->Z * t_trk->wci[3] * tf0t[jc][it]->X
128  + tf0t[ic][it]->X * t_trk->wci[1] * tf0t[jc][it]->Y
129  + tf0t[ic][it]->Y * t_trk->wci[2] * tf0t[jc][it]->Y
130  + tf0t[ic][it]->Z * t_trk->wci[4] * tf0t[jc][it]->Y
131  + tf0t[ic][it]->X * t_trk->wci[3] * tf0t[jc][it]->Z
132  + tf0t[ic][it]->Y * t_trk->wci[4] * tf0t[jc][it]->Z
133  + tf0t[ic][it]->Z * t_trk->wci[5] * tf0t[jc][it]->Z
134  + tf0t[jc][it]->X * t_trk->wci[0] * tf0t[ic][it]->X
135  + tf0t[jc][it]->Y * t_trk->wci[1] * tf0t[ic][it]->X
136  + tf0t[jc][it]->Z * t_trk->wci[3] * tf0t[ic][it]->X
137  + tf0t[jc][it]->X * t_trk->wci[1] * tf0t[ic][it]->Y
138  + tf0t[jc][it]->Y * t_trk->wci[2] * tf0t[ic][it]->Y
139  + tf0t[jc][it]->Z * t_trk->wci[4] * tf0t[ic][it]->Y
140  + tf0t[jc][it]->X * t_trk->wci[3] * tf0t[ic][it]->Z
141  + tf0t[jc][it]->Y * t_trk->wci[4] * tf0t[ic][it]->Z
142  + tf0t[jc][it]->Z * t_trk->wci[5] * tf0t[ic][it]->Z;
143  }
144  }
145  }
146 
147 /* Solving of system DENOM(i,j)*COEF(j)=ANUM(i) */
148 /* for lagrange couplings */
149 /*-------------------------------------------------*/
150  //This is deliberately written without make_unqiue to bypass auto intialization!!
151  auto coef = std::unique_ptr<double[]>(new double[totNC]);
152  if (totNC == 1) {
153  if (denom[0][0] != 0.) {
154  coef[0] = anum[0] / denom[0][0];
155  } else {
156  coef[0] = 1.e3;
157  }
158  } else {
159  //This is deliberately written without make_unqiue to bypass auto intialization!!
160  std::unique_ptr<double[]> adenom( new double[totNC*totNC] );
161 
162  for (ic=0; ic<totNC; ic++) {
163  for (jc=0; jc<totNC; jc++) {
164  adenom[ic*totNC + jc] = denom[ic][jc];
165  }
166  }
167 /* -- INVERT */
168  dsinv(totNC, adenom.get(), totNC, &IERR);
169  if (IERR) {
170  return IERR;
171  }
172  for (ic=0; ic<totNC; ++ic) {
173  coef[ic] = 0.;
174  for (jc=0; jc<totNC; ++jc) {
175  index = ic*totNC + jc;
176  coef[ic] += adenom[index] * anum[jc];
177  }
178  }
179  }
180 
181 
182 /* new vertex */
183  dxyz[0] = 0.;
184  dxyz[1] = 0.;
185  dxyz[2] = 0.;
186 
187  //This is deliberately written without make_unqiue to bypass auto intialization!!
188  auto tmpVec = std::unique_ptr<Vect3DF[]>(new Vect3DF[totNC]);
189  for (ic=0; ic<totNC; ++ic) {
190  tmpVec[ic].X = vk->fitVcov[0] * al0[ic].X
191  + vk->fitVcov[1] * al0[ic].Y
192  + vk->fitVcov[3] * al0[ic].Z;
193  tmpVec[ic].Y = vk->fitVcov[1] * al0[ic].X
194  + vk->fitVcov[2] * al0[ic].Y
195  + vk->fitVcov[4] * al0[ic].Z;
196  tmpVec[ic].Z = vk->fitVcov[3] * al0[ic].X
197  + vk->fitVcov[4] * al0[ic].Y
198  + vk->fitVcov[5] * al0[ic].Z;
199  dxyz[0] += coef[ic] * tmpVec[ic].X;
200  dxyz[1] += coef[ic] * tmpVec[ic].Y;
201  dxyz[2] += coef[ic] * tmpVec[ic].Z;
202  }
203  vk->fitV[0] = vk->iniV[0] + vk->dxyz0[0] + dxyz[0];
204  vk->fitV[1] = vk->iniV[1] + vk->dxyz0[1] + dxyz[1];
205  vk->fitV[2] = vk->iniV[2] + vk->dxyz0[2] + dxyz[2];
206 
207 //std::cout <<"New vrt="<<vk->dxyz0[0]<<", "<<vk->dxyz0[1]<<", "<<vk->dxyz0[2]<<'\n';
208 //std::cout<<"Ncntvrt shift="<<dxyz[0]<<", "<<dxyz[1]<<", "<<dxyz[2]<<'\n';
209 
210 /* new momenta */
211  for (it=0; it<NTRK; ++it) {
212  TWRK * t_trk=vk->tmpArr[it].get();
213  tmp[0] = 0.;
214  tmp[1] = 0.;
215  tmp[2] = 0.;
216  for (ic=0; ic<totNC; ++ic) {
217  tmp[0] += coef[ic] * ( tf0t[ic][it]->X + t_trk->wb[0]*tmpVec[ic].X
218  + t_trk->wb[1]*tmpVec[ic].Y
219  + t_trk->wb[2]*tmpVec[ic].Z);
220  tmp[1] += coef[ic] * ( tf0t[ic][it]->Y + t_trk->wb[3]*tmpVec[ic].X
221  + t_trk->wb[4]*tmpVec[ic].Y
222  + t_trk->wb[5]*tmpVec[ic].Z);
223  tmp[2] += coef[ic] * ( tf0t[ic][it]->Z + t_trk->wb[6]*tmpVec[ic].X
224  + t_trk->wb[7]*tmpVec[ic].Y
225  + t_trk->wb[8]*tmpVec[ic].Z);
226  }
227 
228 /* calculation WCI * TMP */
229  VKTrack * trk = vk->TrackList[it].get();
230  trk->fitP[0] -= t_trk->wci[0]*tmp[0] + t_trk->wci[1]*tmp[1] + t_trk->wci[3]*tmp[2];
231  trk->fitP[1] -= t_trk->wci[1]*tmp[0] + t_trk->wci[2]*tmp[1] + t_trk->wci[4]*tmp[2];
232  double concor = t_trk->wci[3]*tmp[0] + t_trk->wci[4]*tmp[1] + t_trk->wci[5]*tmp[2];
233  if( (trk->fitP[2] - concor)*trk->fitP[2] < 0) { concor=trk->fitP[2]/4.; } // VK 15.12.2009 - Protection against change of sign
234  trk->fitP[2] -= concor;
235 //std::cout<<" Npar="<<trk->fitP[0] <<", "<<trk->fitP[1]<<", "<<trk->fitP[2]<<'\n';
236  }
237 /* =================================================================== */
238  return 0;
239 }
240 
241 //
242 // Get sum of squared aa[ic] for all constraints. It's needed for postfit.
243 //
244 double getCnstValues2( VKVertex * vk ) noexcept
245 {
246  if (vk->ConstraintList.empty()) return 0.;
247  double sumSQ=0.;
248  for(int ii=0; ii<(int)vk->ConstraintList.size();ii++){
249  for(int ic=0; ic<(int)vk->ConstraintList[ii]->NCDim; ic++){
250  sumSQ += (vk->ConstraintList[ii]->aa[ic])*(vk->ConstraintList[ii]->aa[ic]);
251  }
252  }
253  return sumSQ;
254 }
255 
256 
257 } /* End of Trk namespace*/
258 
CommonPars.h
Trk::VKVertex::tmpArr
std::vector< std::unique_ptr< TWRK > > tmpArr
Definition: TrkVKalVrtCoreBase.h:168
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
VtCFitC.h
Trk::TWRK
Definition: TrkVKalVrtCoreBase.h:44
index
Definition: index.py:1
skel.it
it
Definition: skel.GENtoEVGEN.py:396
Trk::VKTrack
Definition: TrkVKalVrtCoreBase.h:64
Trk::vtcfitc
int vtcfitc(VKVertex *vk)
Definition: VtCFitC.cxx:16
Trk::VKVertex::ConstraintList
std::vector< std::unique_ptr< VKConstraintBase > > ConstraintList
Definition: TrkVKalVrtCoreBase.h:169
Trk::VKTrack::fitP
double fitP[3]
Definition: TrkVKalVrtCoreBase.h:76
TrkVKalVrtCoreBase.h
Trk::VKVertex::TrackList
std::vector< std::unique_ptr< VKTrack > > TrackList
Definition: TrkVKalVrtCoreBase.h:167
Trk::VKVertex::fitV
double fitV[3]
Definition: TrkVKalVrtCoreBase.h:140
Trk::dsinv
void dsinv(long int n, double *a, long int DIM, long int *ifail) noexcept
Definition: Tracking/TrkVertexFitter/TrkVKalVrtCore/src/Matrix.cxx:517
Trk::TWRK::wb
double wb[9]
Definition: TrkVKalVrtCoreBase.h:51
Matrix.h
Trk::VKVertex::fitVcov
double fitVcov[6]
Definition: TrkVKalVrtCoreBase.h:141
DeMoUpdate.tmp
string tmp
Definition: DeMoUpdate.py:1167
compute_lumi.denom
denom
Definition: compute_lumi.py:76
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::Vect3DF
Definition: TrkVKalUtils.h:14
Trk::getCnstValues2
double getCnstValues2(VKVertex *vk) noexcept
Definition: VtCFitC.cxx:244
DeMoScan.index
string index
Definition: DeMoScan.py:364
Trk::TWRK::wci
double wci[6]
Definition: TrkVKalVrtCoreBase.h:53
Trk::VKVertex::dxyz0
double dxyz0[3]
Definition: TrkVKalVrtCoreBase.h:165
Trk::VKVertex
Definition: TrkVKalVrtCoreBase.h:128
Trk::VKVertex::iniV
double iniV[3]
Definition: TrkVKalVrtCoreBase.h:142
ForCFT.h
WriteCalibToCool.coef
coef
Definition: WriteCalibToCool.py:582
Generated on Thu Nov 7 2024 21:31:49 for ATLAS Offline Software by   doxygen 1.8.18