Fork me on GitHub

source: git/external/TrackCovariance/SolGridCov.cc@ e1e72fd

Last change on this file since e1e72fd was 170a11d, checked in by michele <michele.selvaggi@…>, 4 years ago

included acceptance and hits in TrackCovariance

  • Property mode set to 100644
File size: 5.9 KB
Line 
1#include <iostream>
2
3#include <TMath.h>
4#include <TVectorD.h>
5#include <TVector3.h>
6#include <TMatrixDSym.h>
7#include <TDecompChol.h>
8#include <TMatrixDSymEigen.h>
9
10#include "SolGridCov.h"
11#include "SolGeom.h"
12#include "SolTrack.h"
13
14using namespace std;
15
16SolGridCov::SolGridCov()
17{
18 // Define pt-polar angle grid
19 fNpt = 22;
20 fPta.ResizeTo(fNpt);
21 Double_t p[] = { 0.1, 0.2, 0.5, 0.7, 1., 2., 3., 4., 6., 8., 10., 15.,
22 20., 25., 30., 40., 50., 60., 80., 100., 150., 200. };
23 for (Int_t ip = 0; ip < fNpt; ip++) fPta(ip) = p[ip];
24
25 fNang = 13;
26 fAnga.ResizeTo(fNang);
27 Double_t a[] = { 10., 15., 20., 25., 30., 35., 40., 45., 50., 60., 70., 80., 90. };
28 for (Int_t ia = 0; ia < fNang; ia++) fAnga(ia) = a[ia];
29 fCov = new TMatrixDSym[fNpt * fNang];
30 for (Int_t ip = 0; ip < fNpt; ip++)
31 {
32 for (Int_t ia = 0; ia < fNang; ia++) fCov[ip * fNang + ia].ResizeTo(5, 5);
33 }
34}
35
36SolGridCov::~SolGridCov()
37{
38 delete[] fCov;
39 delete fAcc;
40}
41
42void SolGridCov::Calc(SolGeom *G)
43{
44 TVectorD pta = fPta;
45 TVectorD anga = fAnga;
46 Bool_t Res = kTRUE; Bool_t MS = kTRUE; // Resolution and multiple scattering flags
47 for (Int_t ip = 0; ip < fNpt; ip++) // Loop on pt grid
48 {
49 Int_t ipt = TMath::Nint(10 * pta(ip));
50 for (Int_t ia = 0; ia < fNang; ia++) // Loop on angle grid
51 {
52 Double_t th = TMath::Pi() * (anga(ia)) / 180.;
53 Double_t x[3], p[3];
54 x[0] = 0; x[1] = 0; x[2] = 0; // Set origin
55 p[0] = pta(ip); p[1] = 0; p[2] = pta(ip) / TMath::Tan(th);
56 //
57 SolTrack *tr = new SolTrack(x, p, G); // Initialize track
58 tr->CovCalc(Res, MS); // Calculate covariance
59 fCov[ip * fNang + ia] = tr->Cov(); // Get covariance
60 }
61 }
62
63// Now make acceptance
64fAcc = new AcceptanceClx(G);
65}
66
67
68//
69Bool_t SolGridCov::IsAccepted(Double_t pt, Double_t Theta)
70{
71 //
72 // pt in GeV, Theta in degrees
73 //
74 Bool_t Accept = kFALSE;
75 if (fAcc->HitNumber(pt, Theta) >= fNminHits)Accept = kTRUE;
76 //
77 return Accept;
78}
79//
80Bool_t SolGridCov::IsAccepted(Double_t *p)
81{
82 //
83 // pt in GeV, Theta in degrees
84 //
85 Bool_t Accept = kFALSE;
86 Double_t pt = TMath::Sqrt(p[0] * p[0] + p[1] * p[1]);
87 Double_t th = 180. * TMath::ATan2(pt, p[2])/TMath::Pi();
88 if (fAcc->HitNumber(pt,th) >= fNminHits)Accept = kTRUE;
89 //
90 return Accept;
91}
92//
93Bool_t SolGridCov::IsAccepted(TVector3 p)
94{
95 //
96 // pt in GeV, Theta in degrees
97 //
98 Bool_t Accept = kFALSE;
99 Double_t pt = p.Pt();
100 Double_t th = 180.*TMath::ACos(p.CosTheta())/TMath::Pi();
101 if (fAcc->HitNumber(pt,th) >= fNminHits)Accept = kTRUE;
102 //
103 return Accept;
104}
105
106
107// Find bin in grid
108Int_t SolGridCov::GetMinIndex(Double_t xval, Int_t N, TVectorD x)
109{
110 Int_t min = -1; // default for xval below the lower limit
111 if (xval < x(0))return min;
112 if (xval>x(N - 1)){ min = N; return min; }
113 for (Int_t i = 0; i < N; i++) if (xval>x(i))min = i;
114 return min;
115}
116// Force positive definitness in normalized matrix
117TMatrixDSym SolGridCov::MakePosDef(TMatrixDSym NormMat)
118{
119 // Input: symmetric matrix with 1's on diagonal
120 // Output: positive definite matrix with 1's on diagonal
121
122 // Default return value
123 TMatrixDSym rMatN = NormMat;
124 // Check the diagonal
125 Bool_t Check = kFALSE;
126 Int_t Size = NormMat.GetNcols();
127 for (Int_t i = 0; i < Size; i++)if (TMath::Abs(NormMat(i, i) - 1.0)>1.0E-15)Check = kTRUE;
128 if (Check)
129 {
130 cout << "SolGridCov::MakePosDef: input matrix doesn ot have 1 on diagonal. Abort." << endl;
131 return rMatN;
132 }
133 // Diagonalize matrix
134 TMatrixDSymEigen Eign(NormMat);
135 TMatrixD U = Eign.GetEigenVectors();
136 TVectorD lambda = Eign.GetEigenValues();
137 // Reset negative eigenvalues to small positive value
138 TMatrixDSym D(Size); D.Zero(); Double_t eps = 1.0e-13;
139 for (Int_t i = 0; i < Size; i++)
140 {
141 D(i, i) = lambda(i);
142 if (lambda(i) <= 0) D(i, i) = eps;
143 }
144 // Rebuild matrix
145 TMatrixD Ut(TMatrixD::kTransposed, U);
146 TMatrixD rMat = (U*D)*Ut; // Now it is positive defite
147 // Restore all ones on diagonal
148 for (Int_t i1 = 0; i1 < Size; i1++)
149 {
150 Double_t rn1 = TMath::Sqrt(rMat(i1, i1));
151 for (Int_t i2 = 0; i2 <= i1; i2++)
152 {
153 Double_t rn2 = TMath::Sqrt(rMat(i2, i2));
154 rMatN(i1, i2) = 0.5*(rMat(i1, i2) + rMat(i2, i1)) / (rn1*rn2);
155 rMatN(i2, i1) = rMatN(i1, i2);
156 }
157 }
158 return rMatN;
159}
160// Interpolate covariance matrix: Bi-linear interpolation
161TMatrixDSym SolGridCov::GetCov(Double_t pt, Double_t ang)
162{
163 // pt in GeV and ang in degrees
164 Int_t minPt = GetMinIndex(pt, fNpt, fPta);
165 if (minPt == -1)minPt = 0;
166 if (minPt == fNpt - 1)minPt = fNpt - 2;
167 Double_t dpt = fPta(minPt + 1) - fPta(minPt);
168 // Put ang in 0-90 range
169 ang = TMath::Abs(ang);
170 while (ang > 90.)ang -= 90.; // Needs to be fixed
171 Int_t minAng = GetMinIndex(ang, fNang, fAnga);
172 if (minAng == -1)minAng = 0;
173 if (minAng == fNang - 1)minAng = fNang - 2;
174 Double_t dang = fAnga(minAng + 1) - fAnga(minAng);
175 //
176 Double_t tpt = (pt - fPta(minPt)) / dpt;
177 Double_t tang = (ang - fAnga(minAng)) / dang;
178 //
179 TMatrixDSym C11 = fCov[minPt * fNang + minAng];
180 TMatrixDSym C12 = fCov[minPt * fNang + minAng + 1];
181 TMatrixDSym C21 = fCov[(minPt + 1) * fNang + minAng];
182 TMatrixDSym C22 = fCov[(minPt + 1) * fNang + minAng + 1];
183 TMatrixDSym Cv = ((1-tpt) * (1-tang)) * C11 +
184 ((1-tpt) * tang ) * C12 +
185 ( tpt * (1-tang)) * C21 +
186 ( tpt * tang ) * C22;
187 // Check for positive definiteness
188 TMatrixDSym CvN = Cv;
189 TMatrixDSym DCvInv(5); DCvInv.Zero();
190 for (Int_t id = 0; id < 5; id++) DCvInv(id, id) = 1.0 / TMath::Sqrt(Cv(id, id));
191 CvN.Similarity(DCvInv); // Normalize diagonal to 1
192 TDecompChol Chl(CvN);
193 if (!Chl.Decompose())
194 {
195 cout << "SolGridCov::GetCov: Interpolated matrix not positive definite. Recovering ...." << endl;
196 TMatrixDSym rCv = MakePosDef(CvN); CvN = rCv;
197 TMatrixDSym DCv(5); DCv.Zero();
198 for (Int_t id = 0; id < 5; id++) DCv(id, id) = TMath::Sqrt(Cv(id, id));
199 Cv = CvN.Similarity(DCv); // Normalize diagonal to 1
200 }
201
202 return Cv;
203}
Note: See TracBrowser for help on using the repository browser.