BmnRoot
Loading...
Searching...
No Matches
BmnBeamTracking.cxx
Go to the documentation of this file.
1#include "BmnBeamTracking.h"
2
3#include "BmnKalmanFilter.h"
4#include "BmnMath.h"
5#include "BmnSiBTHit.h"
6#include "FairRunAna.h"
7#include "FairTrackParam.h"
8#include "TFile.h"
9#include "TH1F.h"
10#include "TStyle.h"
11
12#include <TMath.h>
13#include <map>
14#include <vector>
15
16using namespace std;
17using namespace TMath;
18
19BmnBeamTracking::BmnBeamTracking(Int_t run)
20{
21
22 fEventNo = 0;
23 fSiBTHitsArray = nullptr;
24 fBeamTracksArray = nullptr;
25 fField = nullptr;
26 fSiBTHitsBranchName = "BmnSiBTHit";
27 fBeamTracksBranchName = "BmnBeamTrack";
28 fSiBTDetector = nullptr;
29}
30
31BmnBeamTracking::~BmnBeamTracking()
32{
33 delete fKalman;
34 delete fSiBTDetector;
35}
36
37InitStatus BmnBeamTracking::Init()
38{
39 if (fVerbose > 1)
40 cout << "======================== BEAM tracking init started ====================" << endl;
41
42 fKalman = new BmnKalmanFilter();
43
44 // Get ROOT Manager
45 FairRootManager* ioman = FairRootManager::Instance();
46 if (!ioman)
47 Fatal("Init", "FairRootManager is not instantiated");
48
49 // SILICON
50 fSiBTHitsArray = (TClonesArray*)ioman->GetObject(fSiBTHitsBranchName); // in
51
52 if (!fSiBTHitsArray) {
53 cout << "BmnBeamTracking::Init(): branch " << fSiBTHitsBranchName << " not found! Task will be deactivated"
54 << endl;
55 SetActive(kFALSE);
56 return InitStatus::kERROR;
57 }
58
59 fBeamTracksArray = new TClonesArray("BmnTrack", 100); // out
60 ioman->Register(fBeamTracksBranchName, "BEAM", fBeamTracksArray, kTRUE);
61
62 fField = FairRunAna::Instance()->GetField();
63 if (!fField)
64 Fatal("Init", "No Magnetic Field found");
65
66 TString gPathConfig = gSystem->Getenv("VMCWORKDIR");
67
68 TString gPathGemConfig = gPathConfig + "/parameters/SiBT/XMLConfigs/";
69 TString confSiBT = "SiBTRun8.xml";
70 fSiBTDetector = new BmnSiBTStationSet(gPathGemConfig + confSiBT);
71
72 if (fVerbose > 1)
73 cout << "======================== BEAM tracking init finished ===================" << endl;
74
75 return kSUCCESS;
76}
77
78void BmnBeamTracking::Exec(Option_t* opt)
79{
80 if (fVerbose > 1)
81 cout << "\n======================== BEAM tracking exec started ====================" << endl;
82 if (fVerbose > 1)
83 cout << "\n Event number: " << fEventNo << endl;
84
85 if (!IsActive())
86 return;
87
88 fBeamTracksArray->Delete();
89
90 fEventNo++;
91 FindTracks();
92 if (fVerbose > 0)
93 cout << "BmnBeamTracking: " << fBeamTracksArray->GetEntriesFast() << " tracks" << endl;
94 if (fVerbose > 1)
95 cout << "\n======================== BEAM tracking exec finished ===================" << endl;
96}
97
98BmnStatus BmnBeamTracking::FindTracks()
99{
100
101 vector<BmnSiBTHit*> sortedHits[fSiBTDetector->GetNStations()];
102 SortHits(sortedHits);
103 Double_t Xv = 0.0;
104 Double_t Yv = 0.0;
105 Double_t Zv = 0.0;
106
107 const Int_t nTrMax = 10;
108 Int_t nTr = 0;
109
110 BmnTrack sortedCandidates[nTrMax];
111
112 for (BmnHit* hit2 : sortedHits[2]) {
113 for (BmnHit* hit1 : sortedHits[1]) {
114 for (BmnHit* hit0 : sortedHits[0]) {
115 BmnTrack cand;
116 cand.AddHit(hit0->GetIndex(), hit0);
117 cand.AddHit(hit1->GetIndex(), hit1);
118 cand.AddHit(hit2->GetIndex(), hit2);
119 cand.SortHits();
120
121 if (CalculateTrackParams(&cand) != kBMNERROR) {
122 if (nTr == nTrMax) {
123 for (Int_t i = 0; i < nTr; ++i) {
124 if (cand.GetChi2() < sortedCandidates[i].GetChi2()) {
125 sortedCandidates[i] = cand;
126 }
127 }
128 } else {
129 sortedCandidates[nTr] = cand;
130 nTr++;
131 }
132 }
133 }
134 }
135 }
136
137 for (Int_t i = 0; i < nTr; ++i) {
138
139 BmnTrack tr = sortedCandidates[i];
140
141 TrackUpdateByKalman(&tr);
142 FairTrackParam par = *(tr.GetParamFirst());
143
144 for (Int_t iHit = 0; iHit < tr.GetNHits(); ++iHit) {
145 BmnHit* hit = (BmnHit*)fSiBTHitsArray->At(tr.GetHitIndex(iHit));
146 fKalman->TGeoTrackPropagate(&par, hit->GetZ(), 2212, nullptr, nullptr, 1);
147 hit->SetResXY(par.GetX() - hit->GetX(), par.GetY() - hit->GetY());
148 }
149 fKalman->TGeoTrackPropagate(&par, Zv, 2212, nullptr, nullptr, 1);
150 tr.SetB(Sqrt(Sq(par.GetX() - Xv) + Sq(par.GetY() - Yv)));
151 tr.SetParamLast(par); // ParamLast of beamtrack is now on target
152 new ((*fBeamTracksArray)[fBeamTracksArray->GetEntriesFast()]) BmnTrack(tr);
153 }
154
155 return kBMNSUCCESS;
156}
157
158void BmnBeamTracking::Finish() {}
159
160BmnStatus BmnBeamTracking::SortHits(vector<BmnSiBTHit*>* sortedHits)
161{
162 for (Int_t iHit = 0; iHit < fSiBTHitsArray->GetEntriesFast(); ++iHit) {
163 BmnSiBTHit* hit = (BmnSiBTHit*)fSiBTHitsArray->At(iHit);
164 if (!hit)
165 continue;
166 if (hit->GetStripTotalSignalInLowerLayer() < 200 || hit->GetStripTotalSignalInUpperLayer() < 200)
167 continue;
168 Int_t station = hit->GetStation();
169 sortedHits[station].push_back(hit);
170 }
171
172 return kBMNSUCCESS;
173}
174
175BmnStatus BmnBeamTracking::TrackUpdateByKalman(BmnTrack* cand)
176{
177
178 FairTrackParam par = *(cand->GetParamFirst());
179 Double_t chiTot = 0.0;
180 for (Int_t iHit = 0; iHit < cand->GetNHits(); ++iHit) {
181 BmnHit* hit = (BmnHit*)fSiBTHitsArray->At(cand->GetHitIndex(iHit));
182 Double_t chi = 0.0;
183 fKalman->TGeoTrackPropagate(&par, hit->GetZ(), 2212, nullptr, nullptr, 1);
184 fKalman->Update(&par, hit, chi);
185 chiTot += chi;
186 }
187 cand->SetParamLast(par);
188
189 for (Int_t iHit = cand->GetNHits() - 1; iHit >= 0; iHit--) {
190 BmnHit* hit = (BmnHit*)fSiBTHitsArray->At(cand->GetHitIndex(iHit));
191 Double_t chi = 0.0;
192 fKalman->TGeoTrackPropagate(&par, hit->GetZ(), 2212, nullptr, nullptr, 1);
193 fKalman->Update(&par, hit, chi);
194 chiTot += chi;
195 }
196 cand->SetParamFirst(par);
197 cand->SetChi2(chiTot);
198 //---to vertex for online monitoring------
199 // fKalman->TGeoTrackPropagate(&par, 0, 2212, nullptr, nullptr);
200 // cand->SetParamLast(par);
201
202 return kBMNSUCCESS;
203}
204
205BmnStatus BmnBeamTracking::CalcCovMatrix(BmnTrack* tr)
206{
207 const UInt_t nHits = 3;
208 Double_t chi2circ = 0.0;
209 TVector3 CircParZX = CircleFit(tr, fSiBTHitsArray, chi2circ);
210
211 Double_t Xc = CircParZX.Y(); // x-coordinate of fit-circle center
212 Double_t Zc = CircParZX.X(); // z-coordinate of fit-circle center
213 fField = FairRunAna::Instance()->GetField();
214 const Double_t B = (LineFit(tr, fSiBTHitsArray, "ZY")).X(); // angle coefficient for helicoid
215
216 // Double_t Q = (Xc > 0) ? +1 : -1;
217
218 Double_t sumX(0.0), sumY(0.0), sumTx(0.0), sumTy(0.0), sumQP(0.0);
219 Double_t sumXX(0.0), sumXY(0.0), sumXTx(0.0), sumXTy(0.0), sumXQP(0.0);
220 Double_t sumYY(0.0), sumYTx(0.0), sumYTy(0.0), sumYQP(0.0);
221 Double_t sumTxTx(0.0), sumTxTy(0.0), sumTxQP(0.0);
222 Double_t sumTyTy(0.0), sumTyQP(0.0);
223 Double_t sumQPQP(0.0);
224
225 for (UInt_t i = 0; i < nHits; ++i) {
226 BmnHit* hit = (BmnHit*)fSiBTHitsArray->At(tr->GetHitIndex(i));
227 if (!hit)
228 continue;
229 Double_t Xi = hit->GetX();
230 Double_t Yi = hit->GetY();
231 Double_t Zi = hit->GetZ();
232
233 Double_t Txi = -1.0 * (Zi - Zc) / (Xi - Xc);
234 Double_t Tyi = B;
235 // Double_t Ri = Sqrt(Sq(Xi - Xc) + Sq(Zi - Zc));
236 Double_t QPi = 1.0 / 10.678; // Q / (0.0003 * Abs(fField->GetBy(Xi, Yi, Zi)) * Ri);
237
238 sumX += Xi;
239 sumY += Yi;
240 sumTx += Txi;
241 sumTy += Tyi;
242 sumQP += QPi;
243 sumXX += Xi * Xi;
244 sumXY += Xi * Yi;
245 sumXTx += Xi * Txi;
246 sumXTy += Xi * Tyi;
247 sumXQP += Xi * QPi;
248 sumYY += Yi * Yi;
249 sumYTx += Yi * Txi;
250 sumYTy += Yi * Tyi;
251 sumYQP += Yi * QPi;
252 sumTxTx += Txi * Txi;
253 sumTxTy += Txi * Tyi;
254 sumTxQP += Txi * QPi;
255 sumTyTy += Tyi * Tyi;
256 sumTyQP += Tyi * QPi;
257 sumQPQP += QPi * QPi;
258 }
259
260 Double_t meanX = sumX / nHits;
261 Double_t meanY = sumY / nHits;
262 Double_t meanTx = sumTx / nHits;
263 Double_t meanTy = sumTy / nHits;
264 Double_t meanQP = sumQP / nHits;
265
266 Double_t Cov_X_X = sumXX / nHits - Sq(meanX);
267 Double_t Cov_X_Y = sumXY / nHits - meanX * meanY;
268 Double_t Cov_X_Tx = sumXTx / nHits - meanX * meanTx;
269 Double_t Cov_X_Ty = sumXTy / nHits - meanX * meanTy;
270 Double_t Cov_X_Qp = sumXQP / nHits - meanX * meanQP;
271 Double_t Cov_Y_Y = sumYY / nHits - Sq(meanY);
272 Double_t Cov_Y_Tx = sumYTx / nHits - meanY * meanTx;
273 Double_t Cov_Y_Ty = sumYTy / nHits - meanY * meanTy;
274 Double_t Cov_Y_Qp = sumYQP / nHits - meanY * meanQP;
275 Double_t Cov_Tx_Tx = sumTxTx / nHits - Sq(meanTx);
276 Double_t Cov_Tx_Ty = sumTxTy / nHits - meanTx * meanTy;
277 Double_t Cov_Tx_Qp = sumTxQP / nHits - meanTx * meanQP;
278 Double_t Cov_Ty_Ty = sumTyTy / nHits - Sq(meanTy);
279 Double_t Cov_Ty_Qp = sumTyQP / nHits - meanTy * meanQP;
280 Double_t Cov_Qp_Qp = sumQPQP / nHits - Sq(meanQP);
281
282 FairTrackParam par;
283 par.SetCovariance(0, 0, Cov_X_X);
284 par.SetCovariance(0, 1, Cov_X_Y);
285 par.SetCovariance(0, 2, Cov_X_Tx);
286 par.SetCovariance(0, 3, Cov_X_Ty);
287 par.SetCovariance(0, 4, Cov_X_Qp);
288 par.SetCovariance(1, 1, Cov_Y_Y);
289 par.SetCovariance(1, 2, Cov_Y_Tx);
290 par.SetCovariance(1, 3, Cov_Y_Ty);
291 par.SetCovariance(1, 4, Cov_Y_Qp);
292 par.SetCovariance(2, 2, Cov_Tx_Tx);
293 par.SetCovariance(2, 3, Cov_Tx_Ty);
294 par.SetCovariance(2, 4, Cov_Tx_Qp);
295 par.SetCovariance(3, 3, Cov_Ty_Ty);
296 par.SetCovariance(3, 4, Cov_Ty_Qp);
297 par.SetCovariance(4, 4, Cov_Qp_Qp);
298
299 // если ставить ненулевые ковариации в недиагональных членах, то трек начинает сильно уводить в сторону. Проверить!
300 // par.SetCovariance(0, 0, Cov_X_X);
301 // par.SetCovariance(0, 1, 0.0);
302 // par.SetCovariance(0, 2, 0.0);
303 // par.SetCovariance(0, 3, 0.0);
304 // par.SetCovariance(0, 4, 0.0);
305 // par.SetCovariance(1, 1, Cov_Y_Y);
306 // par.SetCovariance(1, 2, 0.0);
307 // par.SetCovariance(1, 3, 0.0);
308 // par.SetCovariance(1, 4, 0.0);
309 // par.SetCovariance(2, 2, Cov_Tx_Tx);
310 // par.SetCovariance(2, 3, 0.0);
311 // par.SetCovariance(2, 4, 0.0);
312 // par.SetCovariance(3, 3, Cov_Ty_Ty);
313 // par.SetCovariance(3, 4, 0.0);
314 // par.SetCovariance(4, 4, Cov_Qp_Qp);
315
316 tr->SetParamFirst(par);
317 tr->SetParamLast(par);
318
319 return kBMNSUCCESS;
320}
321
322BmnStatus BmnBeamTracking::CalculateTrackParams(BmnTrack* tr)
323{
324 // Estimation of track parameters for events with magnetic field
325
326 TVector3 lineParZY = LineFit(tr, fSiBTHitsArray, "ZY");
327 const Double_t B = lineParZY.X(); // angle coefficient for helicoid
328
329 Double_t Tx_first =
330 CalcTx((BmnHit*)fSiBTHitsArray->At(tr->GetHitIndex(0)), (BmnHit*)fSiBTHitsArray->At(tr->GetHitIndex(1)),
331 (BmnHit*)fSiBTHitsArray->At(tr->GetHitIndex(2)));
332 Double_t Tx_last =
333 CalcTx((BmnHit*)fSiBTHitsArray->At(tr->GetHitIndex(2)), (BmnHit*)fSiBTHitsArray->At(tr->GetHitIndex(1)),
334 (BmnHit*)fSiBTHitsArray->At(tr->GetHitIndex(0)));
335
336 CalcCovMatrix(tr);
337 TVector3 firstPos;
338 TVector3 lastPos;
339 ((BmnHit*)fSiBTHitsArray->At(tr->GetHitIndex(0)))->Position(firstPos);
340 ((BmnHit*)fSiBTHitsArray->At(tr->GetHitIndex(2)))->Position(lastPos);
341 tr->GetParamFirst()->SetPosition(firstPos);
342 tr->GetParamFirst()->SetTx(Tx_first);
343 tr->GetParamFirst()->SetTy(B);
344 tr->GetParamLast()->SetPosition(lastPos);
345 tr->GetParamLast()->SetTx(Tx_last);
346 tr->GetParamLast()->SetTy(B);
347 Double_t QP = 1.0 / 10.678; // CalcQp(tr); 10.678=4.65*124/54 Xe124
348 tr->GetParamFirst()->SetQp(QP);
349 tr->GetParamLast()->SetQp(QP);
350 tr->SetChi2(lineParZY.Z());
351 return kBMNSUCCESS;
352}
353
354Double_t BmnBeamTracking::CalcQp(BmnTrack* track)
355{
356
357 vector<BmnHit*> hits;
358 for (Int_t iHit = 0; iHit < track->GetNHits(); iHit++)
359 hits.push_back((BmnHit*)fSiBTHitsArray->At(track->GetHitIndex(iHit)));
360
361 // Get q/p info from all track segments
362 BmnHit* first = hits[0];
363 BmnHit* second = hits[1];
364 BmnHit* third = hits[2];
365
366 TVector3 CircParZX = CircleBy3Hit(first, second, third);
367 Double_t R = CircParZX.Z();
368 Double_t Xc = CircParZX.Y();
369 Double_t Zc = CircParZX.X();
370
371 Double_t Q = (Xc > 0) ? +1. : -1.;
372 Double_t S = 0.0003
373 * (Abs(fField->GetBy(third->GetX(), third->GetY(), third->GetZ()))
374 + Abs(fField->GetBy(second->GetX(), second->GetY(), second->GetZ()))
375 + Abs(fField->GetBy(first->GetX(), first->GetY(), first->GetZ())))
376 / 3.;
377
378 Double_t Pt = S * R; // actually Pt/Q, but it doesn't matter
379 Double_t fX = first->GetX();
380 Double_t fZ = first->GetZ();
381
382 Double_t h = -1.0;
383
384 Double_t Tx_first = h * (fZ - Zc) / (fX - Xc);
385 TVector3 lineParZY = LineFit(track, fSiBTHitsArray, "ZY");
386 const Double_t B = lineParZY.X(); // angle coefficient for helicoid
387 Double_t Ty_first = B; // / (fX - Xc);
388
389 Double_t Pz = Pt / Sqrt(1 + Sq(Tx_first));
390 Double_t Px = Tx_first * Pz;
391 Double_t Py = Ty_first * Pz;
392 Double_t P = Sqrt(Sq(Px) + Sq(Py) + Sq(Pz));
393 Double_t QP = Q / P;
394
395 return QP;
396}
LineFit
TVector3 LineFit(BmnTrack *track, const TClonesArray *arr, TString type)
Definition BmnMath.cxx:113
CircleBy3Hit
TVector3 CircleBy3Hit(BmnTrack *track, const TClonesArray *arr)
Definition BmnMath.cxx:470
CircleFit
TVector3 CircleFit(BmnTrack *track, const TClonesArray *arr, Double_t &chi2)
Definition BmnMath.cxx:312
CalcTx
Double_t CalcTx(const BmnHit *h0, const BmnHit *h1, const BmnHit *h2)
Definition BmnMath.cxx:422
i
int i
Definition P4_F32vec4.h:22
BmnStatus
BmnStatus
Definition BmnEnums.h:24
kBMNERROR
@ kBMNERROR
Definition BmnEnums.h:26
kBMNSUCCESS
@ kBMNSUCCESS
Definition BmnEnums.h:25
BmnBeamTracking::Finish
virtual void Finish()
Definition BmnBeamTracking.cxx:158
BmnBeamTracking::~BmnBeamTracking
virtual ~BmnBeamTracking()
Definition BmnBeamTracking.cxx:31
BmnBeamTracking::Init
virtual InitStatus Init()
Definition BmnBeamTracking.cxx:37
BmnBeamTracking::Exec
virtual void Exec(Option_t *opt)
Definition BmnBeamTracking.cxx:78
BmnHit::SetResXY
void SetResXY(Double_t resX, Double_t resY)
Definition BmnHit.h:93
BmnHit::GetStation
Short_t GetStation() const
Definition BmnHit.h:45
BmnKalmanFilter::TGeoTrackPropagate
BmnStatus TGeoTrackPropagate(FairTrackParam *par, Double_t zOut, Int_t pdg, vector< Double_t > *F, Double_t *length, Bool_t isField)
Definition BmnKalmanFilter.cxx:480
BmnKalmanFilter::Update
BmnStatus Update(FairTrackParam *par, const BmnHit *hit, Double_t &chiSq)
Definition BmnKalmanFilter.cxx:270
BmnSiBTHit::GetStripTotalSignalInUpperLayer
Double_t GetStripTotalSignalInUpperLayer()
Definition BmnSiBTHit.h:52
BmnSiBTHit::GetStripTotalSignalInLowerLayer
Double_t GetStripTotalSignalInLowerLayer()
Definition BmnSiBTHit.h:48
BmnTrack::SetChi2
void SetChi2(Double_t chi2)
Definition BmnTrack.h:97
BmnTrack::SetParamLast
void SetParamLast(FairTrackParam &par)
Definition BmnTrack.h:89
BmnTrack::GetParamFirst
FairTrackParam * GetParamFirst()
Definition BmnTrack.h:72
BmnTrack::GetChi2
Float_t GetChi2() const
Definition BmnTrack.h:56
BmnTrack::GetNHits
Int_t GetNHits() const
Definition BmnTrack.h:44
BmnTrack::GetHitIndex
Int_t GetHitIndex(Int_t iHit) const
Definition BmnTrack.h:48
BmnTrack::AddHit
void AddHit(Int_t hitIndex, FairHit *Hit)
Definition BmnTrack.cxx:26
BmnTrack::SetB
void SetB(Double_t b)
Definition BmnTrack.h:109
BmnTrack::SetParamFirst
void SetParamFirst(FairTrackParam &par)
Definition BmnTrack.h:85
BmnTrack::GetParamLast
FairTrackParam * GetParamLast()
Definition BmnTrack.h:76
BmnTrack::SortHits
void SortHits()
Definition BmnTrack.cxx:41
run
-clang-format
std
STL namespace.