ATLAS Offline Software
CascadeUtils.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/CascadeUtils.h"
6 #include "TrkVKalVrtCore/Utilities.h"
7 #include "TrkVKalVrtCore/Derivt.h"
8 #include "TrkVKalVrtCore/cfMomentum.h"
9 #include "TrkVKalVrtCore/TrkVKalVrtCoreBase.h"
10 #include "TrkVKalVrtCore/VKalVrtBMag.h"
11 #include <array>
12 #include <cmath>
13 #include <iostream>
14 
15 namespace Trk {
16 
17 // Add to system matrix the derivatives due to pseudotrack constraints
18 int fixPseudoTrackPt(long int NPar, double * fullMtx, double * LSide, CascadeEvent & cascadeEvent_)
19 {
20 
21  int iv,it,ivnext;
22  //Deliberately not make_unique to bypass inititalization
23  std::unique_ptr<double[]> DerivC( new double[NPar] );
24  std::unique_ptr<double[]> DerivP( new double[NPar] );
25  std::unique_ptr<double[]> DerivT( new double[NPar] );
26 //
27  std::vector<double> vMagFld; double vBx,vBy,vBz;
28  for( iv=0; iv<cascadeEvent_.cascadeNV; iv++){
29  VKVertex * vk = cascadeEvent_.cascadeVertexList[iv].get();
30  Trk::vkalMagFld::getMagFld(vk->refIterV[0]+vk->iniV[0], vk->refIterV[1]+vk->iniV[1], vk->refIterV[2]+vk->iniV[2],vBx,vBy,vBz,(vk->vk_fitterControl).get());
31  vMagFld.push_back(vBz); // fill mag.fields for all vertices
32  }
33 //
34  for( iv=0; iv<cascadeEvent_.cascadeNV; iv++){
35  int indCombTrk=-1;
36  int iniPosTrk=0; /* Start of track part of vertex in global matrix */
37  int posCombTrk=0; /* Conbined track position in global matrix */
38  VKVertex* vk = cascadeEvent_.cascadeVertexList[iv].get();
39  if(vk->nextCascadeVrt){ //next vertex exists
40  ivnext=-1; //index of next vertex in common structure
41  for(int ivt=0;ivt<cascadeEvent_.cascadeNV;ivt++)if(vk->nextCascadeVrt==cascadeEvent_.cascadeVertexList[ivt].get())ivnext=ivt;
42  if(ivnext<0){return -1;}; //error in cascade
43 //
44  int NV=vk->nextCascadeVrt->includedVrt.size();
45  if(NV>0){
46  for(it=0; it<NV; it++)
47  if(vk->nextCascadeVrt->includedVrt[it] == vk)
48  indCombTrk=vk->nextCascadeVrt->TrackList.size() - NV + it; // index of combined track in next vertex track list
49  }
50  if(indCombTrk>=0){
51  iniPosTrk =cascadeEvent_.matrixPnt[iv]+3; /*Start of track part of vertex in global matrix */
52  posCombTrk =cascadeEvent_.matrixPnt[ivnext]+3+3*indCombTrk; /*Get position in global matrix */
53  }
54  if(posCombTrk==0 || iniPosTrk==0) {return -1;} //ERROR in cascade structure somewhere....
55 //
56 // Momentum of pseudo track in next vertex
57  for(int ivt=0; ivt<NPar; ivt++)DerivC[ivt]=DerivP[ivt]=DerivT[ivt]=0.;
58  DerivC[posCombTrk+2]=-1.;
59  DerivT[posCombTrk+0]=-1.;
60  DerivP[posCombTrk+1]=-1.;
61  std::array<double, 4> ppsum = getIniParticleMom( vk->nextCascadeVrt->TrackList[indCombTrk].get(), vMagFld[ivnext] ); // INI for pseudo
62  double csum=vk->nextCascadeVrt->TrackList[indCombTrk]->iniP[2]; // INI for pseudo
63  double ptsum=sqrt(ppsum[0]*ppsum[0] + ppsum[1]*ppsum[1]);
64  double sinth2sum=(ppsum[0]*ppsum[0] + ppsum[1]*ppsum[1])/(ppsum[0]*ppsum[0] + ppsum[1]*ppsum[1] + ppsum[2]*ppsum[2]);
65 
66 //
67 // Momenta+Derivatives of tracks in vertex itself
68  for(it=0; it<(int)vk->TrackList.size(); it++){
69  std::array<double, 4> pp = getIniParticleMom( vk->TrackList[it].get(), vMagFld[iv]);
70  double curv=vk->TrackList[it]->iniP[2];
71  double pt=sqrt(pp[0]*pp[0] + pp[1]*pp[1]);
72  double cth=pp[2]/pt;
73  double sinth2=(pp[0]*pp[0] + pp[1]*pp[1])/(pp[0]*pp[0] + pp[1]*pp[1] + pp[2]*pp[2]);
74  DerivC[iniPosTrk+it*3+1] = csum/ptsum/ptsum*(ppsum[0]*pp[1]-ppsum[1]*pp[0]); // dC/dPhi_i
75  DerivC[iniPosTrk+it*3+2] = csum/ptsum/ptsum*(ppsum[0]*pp[0]+ppsum[1]*pp[1])/curv; // dC/dC_i
76  DerivP[iniPosTrk+it*3+1] = (ppsum[0]*pp[0]+ppsum[1]*pp[1])/ptsum/ptsum; // dPhi/dPhi_i
77  DerivP[iniPosTrk+it*3+2] = (ppsum[1]*pp[0]-ppsum[0]*pp[1])/ptsum/ptsum/curv; // dPhi/dC_i
78  DerivT[iniPosTrk+it*3+0] = (sinth2sum*pt)/(sinth2*ptsum); // dTheta/dTheta_i
79  DerivT[iniPosTrk+it*3+2] = (sinth2sum*pt*cth)/(curv*ptsum); // dTheta/dC_i
80  }
81 // double iniV0Curv=myMagFld.getCnvCst()*vMagFld[iv]/sqrt(tpx*tpx+tpy*tpy); //initial PseudoTrack Curvature
82 // if(csum<0)iniV0Curv *= -1.;
83 // iniV0Curv *= vMagFld[ivnext]/vMagFld[iv]; //magnetic field correction
84 //
85 //fill Full Matrix and left side vector
86 //
87 //---- Momentum only
88 // for(it=0; it<NPar; it++) fullMtx[ (NPar-1*(cascadeEvent_.cascadeNV-1)+1*iv )*NPar + it] = DerivC[it];
89 // for(it=0; it<NPar; it++) fullMtx[ (NPar-1*(cascadeEvent_.cascadeNV-1)+1*iv ) + it*NPar] = DerivC[it];
90 // LSide[ NPar-1*(cascadeEvent_.cascadeNV-1)+1*iv] = -iniV0Curv+csum;
91 //---- Momentum+phi+theta //VK seems overshooting because direction is fixed by vertex-vertex pointing. Returns wrong error matrix
92  for(it=0; it<NPar; it++) fullMtx[ (NPar-3*(cascadeEvent_.cascadeNV-1)+3*iv )*NPar + it] = DerivT[it];
93  for(it=0; it<NPar; it++) fullMtx[ (NPar-3*(cascadeEvent_.cascadeNV-1)+3*iv+1)*NPar + it] = DerivP[it];
94  for(it=0; it<NPar; it++) fullMtx[ (NPar-3*(cascadeEvent_.cascadeNV-1)+3*iv+2)*NPar + it] = DerivC[it];
95  for(it=0; it<NPar; it++) fullMtx[ (NPar-3*(cascadeEvent_.cascadeNV-1)+3*iv ) + it*NPar] = DerivT[it];
96  for(it=0; it<NPar; it++) fullMtx[ (NPar-3*(cascadeEvent_.cascadeNV-1)+3*iv+1) + it*NPar] = DerivP[it];
97  for(it=0; it<NPar; it++) fullMtx[ (NPar-3*(cascadeEvent_.cascadeNV-1)+3*iv+2) + it*NPar] = DerivC[it];
98  VKTrack* cmbt=vk->nextCascadeVrt->TrackList[indCombTrk].get();
99  LSide[ NPar-3*(cascadeEvent_.cascadeNV-1)+3*iv ] = cmbt->iniP[0]-cmbt->Perig[2];
100  LSide[ NPar-3*(cascadeEvent_.cascadeNV-1)+3*iv+1] = cmbt->iniP[1]-cmbt->Perig[3];
101  LSide[ NPar-3*(cascadeEvent_.cascadeNV-1)+3*iv+2] = cmbt->iniP[2]-cmbt->Perig[4];
102  }
103 
104  } //end of vertex cycle
105  return 0; //All ok
106 }
107 //---------------------------------------------------------------------------
108 // Returns address of VTrack of combined track for given vertex
109 //
110 VKTrack * getCombinedVTrack(VKVertex * vk)
111 {
112  if(!vk->nextCascadeVrt) return nullptr; //nonpointing vertex
113  int NV=vk->nextCascadeVrt->includedVrt.size();
114  if(NV==0) return nullptr; //Error in structure
115 
116  int itv=-1;
117  for(int it=0; it<NV; it++) if(vk->nextCascadeVrt->includedVrt[it] == vk) {itv=it; break;};
118  if(itv<0) return nullptr; // Not found in list
119 
120  int totNT = vk->nextCascadeVrt->TrackList.size();
121  return vk->nextCascadeVrt->TrackList[totNT - NV + itv].get(); // pointer to combined track in next vertex
122 }
123 
124 
125 //---------------------------------------------------------------------------
126 // Returns dimension of full matrix for cascade fit
127 // At the end the place for pseudotrack(all!) momenta(3) constraints
128 // By default (Type=0) return full cascade NPar.
129 // For Type=1 returns amount of physics parameters only (without all constraints)
130 //
131 // MUST BE CONSISTENT WITH fixPseudoTrackPt(...)!!!
132 //
133 int getCascadeNPar(CascadeEvent & cascadeEvent_, int Type/*=0*/)
134 {
135  int NV=cascadeEvent_.cascadeNV;
136  int NTrk=0;
137  int NCnst=0;
138  for( int iv=0; iv<cascadeEvent_.cascadeNV; iv++){
139  VKVertex *vk = cascadeEvent_.cascadeVertexList[iv].get();
140  NTrk += vk->TrackList.size();
141  for(int ic=0; ic<(int)vk->ConstraintList.size();ic++) NCnst += vk->ConstraintList[ic]->NCDim;
142  }
143  if(Type==1) return 3*(NV+NTrk); // Return amount of physics parameters
144  return 3*(NV+NTrk)+NCnst + 3*(cascadeEvent_.cascadeNV-1); //additional 3 momentum constraints
145  //return 3*(NV+NTrk)+NCnst + 1*(cascadeEvent_.cascadeNV-1); //additional 1 momentum constraints
146 }
147 
148 //
149 // Track parameters are translated at each iteration so iniV==(0,0,0)
150 //
151 void setFittedParameters(const double * result, std::vector<int> & matrixPnt, CascadeEvent & cascadeEvent_)
152 {
153  int iv,it,Pnt;
154  for( iv=0; iv<cascadeEvent_.cascadeNV; iv++){
155  VKVertex *vk = cascadeEvent_.cascadeVertexList[iv].get();
156  Pnt=matrixPnt[iv]; // start of vertex parameters
157  vk->fitV[0]=result[Pnt]; vk->fitV[1]=result[Pnt+1]; vk->fitV[2]=result[Pnt+2];
158  for( it=0; it<(int)vk->TrackList.size(); it++){
159  VKTrack * trk=vk->TrackList[it].get();
160  trk->fitP[0]=trk->iniP[0]+result[Pnt+3+it*3 + 0];
161  trk->fitP[1]=trk->iniP[1]+result[Pnt+3+it*3 + 1];
162  trk->fitP[2]=trk->iniP[2]+result[Pnt+3+it*3 + 2];
163  trk->Chi2 = cfchi2(vk->fitV, trk->fitP, trk );
164  }
165  }
166 }
167 
168 void setFittedMatrices(const double * COVFIT, long int MATRIXSIZE,
169  std::vector<int> & matrixPnt,
170  std::vector< std::vector<double> > & covarCascade,
171  CascadeEvent & cascadeEvent_)
172 {
173  int iv, Pnt, ic, ir, vrtMtxSize, count;
174  for( iv=0; iv<cascadeEvent_.cascadeNV; iv++){
175  VKVertex *vk = cascadeEvent_.cascadeVertexList[iv].get();
176  Pnt=matrixPnt[iv]; // start of vertex parameters
177  vrtMtxSize=3+vk->TrackList.size()*3; //size of matrix for given vertex
178  std::vector<double> Res(vrtMtxSize*(vrtMtxSize+1)/2);
179  count=0;
180  for (int col=0; col<vrtMtxSize; col++){
181  for (int row=0; row<=col; row++){
182  ic=Pnt+col; ir=Pnt+row;
183  Res[count]=COVFIT[ic*MATRIXSIZE + ir]; count++;
184  }
185  }
186  covarCascade.emplace_back(std::move(Res));
187  }
188 }
189 //
190 // Symmetrical indexing (I*(I+1)/2+J) is valid ONLY if I>=J
191 std::vector<double> transformCovar(int NPar, double **Deriv, const std::vector<double> &covarI )
192 {
193  std::vector<double> covarO(NPar*(NPar+1)/2, 0.);
194  int ii,ij,oi,oj, indexO, indexI;
195  for(oi=0; oi<NPar; oi++){
196  for(oj=0; oj<=oi; oj++){
197  indexO = oi*(oi+1)/2+oj;
198  for(ii=0; ii<NPar; ii++){
199  for(ij=0; ij<NPar; ij++){
200  indexI = ii>=ij ? ii*(ii+1)/2+ij : ij*(ij+1)/2+ii;
201  covarO[indexO] += Deriv[oi][ii]*covarI[indexI]*Deriv[oj][ij];
202  } } } }
203  return covarO;
204 }
205 
206 
207 void addCrossVertexDeriv(CascadeEvent & cascadeEvent_, double * ader, long int MATRIXSIZE, const std::vector<int> & matrixPnt)
208 {
209  int iv,ivn;
210  //for( iv=0; iv<cascadeEvent_.cascadeNV; iv++)std::cout<<matrixPnt[iv]<<", ";std::cout<<'\n';
211 
212  for( iv=0; iv<cascadeEvent_.cascadeNV; iv++){
213  VKVertex *vk = cascadeEvent_.cascadeVertexList[iv].get();
214  if(!vk->nextCascadeVrt)continue; //no pointing
215  for( ivn=iv; ivn<cascadeEvent_.cascadeNV; ivn++){
216  if(vk->nextCascadeVrt == cascadeEvent_.cascadeVertexList[ivn].get()) break; //vertex found
217  }
218  if( ivn == cascadeEvent_.cascadeNV ) continue; // no found vertex
219 //
220 // Now we have vertex pair "from"(iv) "to"(ivn)
221  int From=matrixPnt[iv]; // start of "from" information
222  int Next=matrixPnt[iv+1]; // start of "next vertex" information. Pnt.constraints are 2 prevous param.!!!
223  int Into=matrixPnt[ivn]; // start of "to" information
224 //
225 // The same constraints, but derivatives are with respect to other ("to") vertex
226  ader[(0+Into) + (Next-1)*MATRIXSIZE] = - ader[(0+From) + (Next-1)*MATRIXSIZE];
227  ader[(1+Into) + (Next-1)*MATRIXSIZE] = - ader[(1+From) + (Next-1)*MATRIXSIZE];
228  ader[(2+Into) + (Next-1)*MATRIXSIZE] = - ader[(2+From) + (Next-1)*MATRIXSIZE];
229  ader[(0+Into) + (Next-2)*MATRIXSIZE] = - ader[(0+From) + (Next-2)*MATRIXSIZE];
230  ader[(1+Into) + (Next-2)*MATRIXSIZE] = - ader[(1+From) + (Next-2)*MATRIXSIZE];
231  ader[(2+Into) + (Next-2)*MATRIXSIZE] = - ader[(2+From) + (Next-2)*MATRIXSIZE];
232  ader[(0+Into)*MATRIXSIZE + (Next-1)] = - ader[(0+From) + (Next-1)*MATRIXSIZE];
233  ader[(1+Into)*MATRIXSIZE + (Next-1)] = - ader[(1+From) + (Next-1)*MATRIXSIZE];
234  ader[(2+Into)*MATRIXSIZE + (Next-1)] = - ader[(2+From) + (Next-1)*MATRIXSIZE];
235  ader[(0+Into)*MATRIXSIZE + (Next-2)] = - ader[(0+From) + (Next-2)*MATRIXSIZE];
236  ader[(1+Into)*MATRIXSIZE + (Next-2)] = - ader[(1+From) + (Next-2)*MATRIXSIZE];
237  ader[(2+Into)*MATRIXSIZE + (Next-2)] = - ader[(2+From) + (Next-2)*MATRIXSIZE];
238  }
239 }
240 
241 
242 //--------------------------------------------------------------------
243 // Copy matrix Input with dimension IDIM to predefined place(TStart)
244 // into matrix Target with dimension TDIM
245 //
246 void copyFullMtx(const double *Input, long int IPar, long int IDIM,
247  double *Target, long int TStart, long int TDIM)
248 {
249  int i,j,it,jt;
250  for( i=0; i<IPar; i++){
251  for( j=0; j<IPar; j++){
252  it=i+TStart; jt=j+TStart;
253  Target[it*TDIM+jt] = Input[i*IDIM+j];
254  }
255  }
256 }
257 
258 //--------------------------------------------------------------------
259 // Make the convolution Cov=D*OldCov*Dt
260 //
261 void getNewCov(const double *OldCov, const double* Der, double* Cov, long int DIM) noexcept
262 {
263  int i,j,it,jt;
264  for( i=0; i<DIM; i++){
265  for( j=0; j<DIM; j++){
266  Cov[i*DIM+j]=0.;
267  for( it=0; it<DIM; it++){
268  for( jt=0; jt<DIM; jt++){
269  Cov[i*DIM+j] += Der[i*DIM+it]*OldCov[it*DIM+jt]*Der[j*DIM+jt];
270  }
271  }
272  }
273  }
274 }
275 
276 
277 
278 } /* End of namespace Trk*/
query_example.row
row
Definition: query_example.py:24
get_generator_info.result
result
Definition: get_generator_info.py:21
Trk::copyFullMtx
void copyFullMtx(const double *Input, long int IPar, long int IDIM, double *Target, long int TStart, long int TDIM)
Definition: CascadeUtils.cxx:246
Trk::vkalMagFld::getMagFld
static void getMagFld(const double, const double, const double, double &, double &, double &, VKalVrtControlBase *)
Definition: VKalVrtBMag.cxx:24
Trk::transformCovar
std::vector< double > transformCovar(int NPar, double **Deriv, const std::vector< double > &covarI)
Definition: CascadeUtils.cxx:191
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
Trk::VKTrack::Chi2
double Chi2
Definition: TrkVKalVrtCoreBase.h:77
Trk::Target
@ Target
Definition: TargetSurfaces.h:40
skel.it
it
Definition: skel.GENtoEVGEN.py:423
test_pyathena.pt
pt
Definition: test_pyathena.py:11
Trk::VKVertex::refIterV
double refIterV[3]
Definition: TrkVKalVrtCoreBase.h:146
Trk::addCrossVertexDeriv
void addCrossVertexDeriv(CascadeEvent &cascadeEvent_, double *ader, long int MATRIXSIZE, const std::vector< int > &matrixPnt)
Definition: CascadeUtils.cxx:207
Trk::VKTrack
Definition: TrkVKalVrtCoreBase.h:64
Trk::VKTrack::Perig
double Perig[5]
Definition: TrkVKalVrtCoreBase.h:86
VKalVrtBMag.h
Trk::getCascadeNPar
int getCascadeNPar(CascadeEvent &cascadeEvent_, int Type)
Definition: CascadeUtils.cxx:133
Trk::VKVertex::ConstraintList
std::vector< std::unique_ptr< VKConstraintBase > > ConstraintList
Definition: TrkVKalVrtCoreBase.h:169
XMLtoHeader.count
count
Definition: XMLtoHeader.py:85
Trk::VKTrack::fitP
double fitP[3]
Definition: TrkVKalVrtCoreBase.h:76
Trk::CascadeEvent::matrixPnt
std::vector< int > matrixPnt
Definition: TrkVKalVrtCoreBase.h:31
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::VKVertex::includedVrt
std::vector< VKVertex * > includedVrt
Definition: TrkVKalVrtCoreBase.h:180
lumiFormat.i
int i
Definition: lumiFormat.py:92
Trk::getIniParticleMom
std::array< double, 4 > getIniParticleMom(const VKTrack *trk, const VKVertex *vk)
Definition: cfMomentum.cxx:60
Trk::getCombinedVTrack
VKTrack * getCombinedVTrack(VKVertex *vk)
Definition: CascadeUtils.cxx:110
xAODType
Definition: ObjectType.h:13
grepfile.ic
int ic
Definition: grepfile.py:33
Trk::setFittedParameters
void setFittedParameters(const double *result, std::vector< int > &matrixPnt, CascadeEvent &cascadeEvent_)
Definition: CascadeUtils.cxx:151
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
query_example.col
col
Definition: query_example.py:7
Trk::getNewCov
void getNewCov(const double *OldCov, const double *Der, double *Cov, long int DIM) noexcept
Definition: CascadeUtils.cxx:261
ir
int ir
counter of the current depth
Definition: fastadd.cxx:49
CascadeUtils.h
Trk::CascadeEvent
Definition: TrkVKalVrtCoreBase.h:21
Trk::CascadeEvent::cascadeNV
int cascadeNV
Definition: TrkVKalVrtCoreBase.h:23
Trk::cfchi2
double cfchi2(double *xyzt, const long int ich, double *part, const double *par0, double *wgt, double *rmnd)
Definition: Tracking/TrkVertexFitter/TrkVKalVrtCore/src/Utilities.cxx:27
Trk::VKVertex::nextCascadeVrt
VKVertex * nextCascadeVrt
Definition: TrkVKalVrtCoreBase.h:179
Trk::VKVertex
Definition: TrkVKalVrtCoreBase.h:128
Trk::VKVertex::iniV
double iniV[3]
Definition: TrkVKalVrtCoreBase.h:142
Trk::VKTrack::iniP
double iniP[3]
Definition: TrkVKalVrtCoreBase.h:83
Trk::fixPseudoTrackPt
int fixPseudoTrackPt(long int NPar, double *fullMtx, double *LSide, CascadeEvent &cascadeEvent_)
Definition: CascadeUtils.cxx:18
NswErrorCalibData::Input
Helper struct to be parsed to the object to derive the specific error of the cluster.
Definition: NswErrorCalibData.h:25
Utilities.h
Trk::CascadeEvent::cascadeVertexList
std::vector< std::unique_ptr< VKVertex > > cascadeVertexList
Definition: TrkVKalVrtCoreBase.h:30
Trk::setFittedMatrices
void setFittedMatrices(const double *COVFIT, long int MATRIXSIZE, std::vector< int > &matrixPnt, std::vector< std::vector< double > > &covarCascade, CascadeEvent &cascadeEvent_)
Definition: CascadeUtils.cxx:168
cfMomentum.h
Generated on Fri Jul 5 2024 21:12:17 for ATLAS Offline Software by   doxygen 1.8.18