| 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 |
|
|---|
| 14 | using namespace std;
|
|---|
| 15 |
|
|---|
| 16 | SolGridCov::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 |
|
|---|
| 36 | SolGridCov::~SolGridCov()
|
|---|
| 37 | {
|
|---|
| 38 | delete[] fCov;
|
|---|
| 39 | delete fAcc;
|
|---|
| 40 | }
|
|---|
| 41 |
|
|---|
| 42 | void 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
|
|---|
| 64 | fAcc = new AcceptanceClx(G);
|
|---|
| 65 | }
|
|---|
| 66 |
|
|---|
| 67 |
|
|---|
| 68 | //
|
|---|
| 69 | Bool_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 | //
|
|---|
| 80 | Bool_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 | //
|
|---|
| 93 | Bool_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
|
|---|
| 108 | Int_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
|
|---|
| 117 | TMatrixDSym 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
|
|---|
| 161 | TMatrixDSym 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 | }
|
|---|
