Fork me on GitHub

source: git/external/TrackCovariance/ObsTrk.cc@ 7bca620

Last change on this file since 7bca620 was 7bca620, checked in by Franco BEDESCHI <bed@…>, 3 years ago

Add feature to force starting radius in vertex fit
  • Property mode set to 100644
File size: 4.5 KB
Line 
1#include <TMath.h>
2#include <TVector3.h>
3#include <TMatrixD.h>
4#include <TMatrixDSym.h>
5#include <TDecompChol.h>
6#include <TRandom.h>
7#include <iostream>
8#include "SolGeom.h"
9#include "SolGridCov.h"
10#include "ObsTrk.h"
11//
12// Constructors
13//
14// x(3) track origin, p(3) track momentum at origin, Q charge, B ma gnetic field in Tesla
15ObsTrk::ObsTrk(TVector3 x, TVector3 p, Double_t Q, SolGridCov *GC, SolGeom *G)
16{
17 fB = G->B();
18 SetB(fB);
19 fG = G;
20 fGC = GC;
21 fGenX = x;
22 fGenP = p;
23 fGenQ = Q;
24 fGenPar.ResizeTo(5);
25 fGenParMm.ResizeTo(5);
26 fGenParACTS.ResizeTo(6);
27 fGenParILC.ResizeTo(5);
28 fObsPar.ResizeTo(5);
29 fObsParMm.ResizeTo(5);
30 fObsParACTS.ResizeTo(6);
31 fObsParILC.ResizeTo(5);
32 fCov.ResizeTo(5, 5);
33 fCovMm.ResizeTo(5, 5);
34 fCovACTS.ResizeTo(6, 6);
35 fCovILC.ResizeTo(5, 5);
36 fGenPar = XPtoPar(x,p,Q);
37 fGenParMm = ParToMm(fGenPar);
38 fGenParACTS = ParToACTS(fGenPar);
39 fGenParILC = ParToILC(fGenPar);
40 //
41 fObsPar = GenToObsPar(fGenPar);
42 fObsParMm = ParToMm(fObsPar);
43 fObsParACTS = ParToACTS(fObsPar);
44 fObsParILC = ParToILC(fObsPar);
45 fObsX = ParToX(fObsPar);
46 fObsP = ParToP(fObsPar);
47 fObsQ = ParToQ(fObsPar);
48 fCovMm = CovToMm(fCov);
49 fCovACTS = CovToACTS(fObsPar, fCov);
50 fCovILC = CovToILC(fCov);
51}
52//
53// x[3] track origin, p[3] track momentum at origin, Q charge, B magnetic field in Tesla
54ObsTrk::ObsTrk(Double_t *x, Double_t *p, Double_t Q, SolGridCov* GC, SolGeom *G)
55{
56 fB = G->B();
57 SetB(fB);
58 fG = G;
59 fGC = GC;
60 fGenX.SetXYZ(x[0],x[1],x[2]);
61 fGenP.SetXYZ(p[0],p[1],p[2]);
62 fGenQ = Q;
63 fGenPar.ResizeTo(5);
64 fGenParMm.ResizeTo(5);
65 fGenParACTS.ResizeTo(6);
66 fGenParILC.ResizeTo(5);
67 fObsPar.ResizeTo(5);
68 fObsParMm.ResizeTo(5);
69 fObsParACTS.ResizeTo(6);
70 fObsParILC.ResizeTo(5);
71 fCov.ResizeTo(5, 5);
72 fCovMm.ResizeTo(5, 5);
73 fCovACTS.ResizeTo(6, 6);
74 fCovILC.ResizeTo(5, 5);
75 fGenPar = XPtoPar(fGenX, fGenP, Q);
76 fGenParMm = ParToMm(fGenPar);
77 fGenParACTS = ParToACTS(fGenPar);
78 fGenParILC = ParToILC(fGenPar);
79 //
80 fObsPar = GenToObsPar(fGenPar);
81 fObsParMm = ParToMm(fObsPar);
82 fObsParACTS = ParToACTS(fObsPar);
83 fObsParILC = ParToILC(fObsPar);
84 fObsX = ParToX(fObsPar);
85 fObsP = ParToP(fObsPar);
86 fObsQ = ParToQ(fObsPar);
87 fCovMm = CovToMm(fCov);
88 fCovACTS = CovToACTS(fObsPar, fCov);
89 fCovILC = CovToILC(fCov);
90}
91//
92// Destructor
93ObsTrk::~ObsTrk()
94{
95 fGenX.Clear();
96 fGenP.Clear();
97 fGenPar.Clear();
98 fGenParMm.Clear();
99 fGenParACTS.Clear();
100 fGenParILC.Clear();
101 fObsX.Clear();
102 fObsP.Clear();
103 fObsPar.Clear();
104 fObsParMm.Clear();
105 fObsParACTS.Clear();
106 fObsParILC.Clear();
107 fCov.Clear();
108 fCovMm.Clear();
109 fCovACTS.Clear();
110 fCovILC.Clear();
111}
112//
113TVectorD ObsTrk::GenToObsPar(TVectorD gPar)
114{
115 //
116 // Check ranges
117 Double_t minPt = fGC->GetMinPt();
118 //if (pt < minPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is below grid range of " << minPt << std::endl;
119 Double_t maxPt = fGC->GetMaxPt();
120 //if (pt > maxPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is above grid range of " << maxPt << std::endl;
121 Double_t minAn = fGC->GetMinAng();
122 //if (angd < minAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd
123 // << " is below grid range of " << minAn << std::endl;
124 Double_t maxAn = fGC->GetMaxAng();
125 //if (angd > maxAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd
126 // << " is above grid range of " << maxAn << std::endl;
127 //
128 TMatrixDSym Cov(5);
129 //
130 // Check if track origin is inside beampipe and betwen the first disks
131 //
132 Double_t Rin = fG->GetRmin();
133 Double_t ZinPos = fG->GetZminPos();
134 Double_t ZinNeg = fG->GetZminNeg();
135 Bool_t inside = TrkUtil::IsInside(fGenX, Rin, ZinNeg, ZinPos); // Check if in inner box
136 if (inside)
137 {
138 //std::cout<<"ObsTrk:: inside: x= "<<fGenX(0)<<", y= "<<fGenX(1)
139 // <<", z= "<<fGenX(2)<<std::endl;
140 // Observed track parameters
141 Double_t pt = fGenP.Pt();
142 Double_t angd = fGenP.Theta() * 180. / TMath::Pi();
143 Cov = fGC->GetCov(pt, angd); // Track covariance
144 }
145 else
146 {
147 //std::cout<<"ObsTrk:: outside: x= "<<fGenX(0)<<", y= "<<fGenX(1)
148 // <<", z= "<<fGenX(2)<<std::endl;
149 SolTrack* trk = new SolTrack(fGenX, fGenP, fG);
150 Bool_t Res = kTRUE; Bool_t MS = kTRUE;
151 trk->CovCalc(Res, MS); // Calculate covariance matrix
152 Cov = trk->Cov(); // Track covariance
153 delete trk;
154 }
155 //
156 fCov = Cov;
157 //
158 // Now do Choleski decomposition and random number extraction, with appropriate stabilization
159 //
160 TVectorD oPar = TrkUtil::CovSmear(gPar, Cov);
161 //
162 return oPar;
163}
Note: See TracBrowser for help on using the repository browser.