BmnRoot
Loading...
Searching...
No Matches
BmnHypNuclPairFinder.cxx
Go to the documentation of this file.
1
2#include "BmnHypNuclPairFinder.h"
3
4BmnHypNuclPairFinder::BmnHypNuclPairFinder() {}
5
6BmnHypNuclPairFinder::BmnHypNuclPairFinder(TString hyp, Bool_t exp)
7{
8 if (hyp.Contains("H3L") || hyp.Contains("h3l")) {
9 massPos = fMassHe3;
10 massNeg = fMassPi;
11 } else if (hyp.Contains("H4L") || hyp.Contains("h4l")) {
12 massPos = fMassHe4;
13 massNeg = fMassPi;
14 } else {
15 massPos = fMassHe3;
16 massNeg = fMassPi;
17 }
18 fEventCounter = 0;
19 fKalman = nullptr;
20 fParticlePair = nullptr;
21 fVertex = nullptr;
22 fTrigger = nullptr;
23 fStsHits = nullptr;
24 fUpperGemClusters = nullptr;
25 fLowerGemClusters = nullptr;
26 fUpperFsdClusters = nullptr;
27 fLowerFsdClusters = nullptr;
28 fStsTracks = nullptr;
29 fGlobTracks = nullptr;
30 fTof700Hits = nullptr;
31 isExp = exp;
32}
33
34BmnHypNuclPairFinder::~BmnHypNuclPairFinder()
35{
36 // delete fDetector;
37}
38
39void BmnHypNuclPairFinder::Analysis()
40{
41 TLorentzVector lPos, lNeg;
42 Float_t VPX = fVertex->GetX();
43 Float_t VPY = fVertex->GetY();
44 Float_t VPZ = fVertex->GetZ();
45 Float_t dVPXsq = fVertex->GetCovariance(0, 0);
46 Float_t dVPYsq = fVertex->GetCovariance(1, 1);
47
48 vector<FairTrackParam> parInVPPos;
49 vector<FairTrackParam> parInVPNeg;
50 vector<Float_t> m2Pos700;
51 vector<Float_t> chi2NDFPos;
52 vector<Float_t> chi2NDFNeg;
53 vector<Int_t> nHitsPos;
54 vector<Float_t> dedxPos;
55 vector<Int_t> nGemHitsNeg;
56 vector<Int_t> nGemHitsPos;
57 vector<Int_t> nFsdHitsNeg;
58 vector<Int_t> nFsdHitsPos;
59
60 TF1* hypCut = new TF1("hypCut", "20000*TMath::Exp(-2.0*TMath::Sqrt(x)) + 600.0", 0, 15);
61
62 // m^2 cut to select He-4 (low)
63 TF1* m2CutLow = nullptr;
64 // m^2 cut to select He-4 (up)
65 TF1* m2CutUp = nullptr;
66
67 if (massPos > 3) { // He4
68 m2CutLow = new TF1("m2CutLow", "-0.0236*x*x+0.0240*x+3.1896", 0, 30);
69 m2CutUp = new TF1("m2CutUp", "0.0169*x*x+0.1306*x+4.0525", 0, 30);
70 } else { // He3
71 m2CutLow = new TF1("m2CutLow", "-0.0034*x*x-0.0360*x+1.7426", 0, 30);
72 m2CutUp = new TF1("m2CutUp", "0.0128*x*x-0.1060*x+2.6261", 0, 30);
73 }
74
75 for (Int_t i = 0; i < fGlobTracks->GetEntriesFast(); i++) {
76 BmnGlobalTrack* gl = (BmnGlobalTrack*)fGlobTracks->UncheckedAt(i);
77 CbmStsTrack* sts = (CbmStsTrack*)fStsTracks->At(gl->GetGemTrackIndex());
78
79 if (!sts->FromVertex()) // хороший ли это кат для пионов? Гелий скорее всего из вершины летят, а вот пионы -
80 // не ясно...
81 continue;
82
83 if (sts->GetNStsHits() < 4)
84 continue;
85
86 FairTrackParam par = *(sts->GetParamFirst());
87 const Float_t tZ = -0.01942; // z-coordinate of target for RUN-8
88 if (fKalman->TGeoTrackPropagate(&par, tZ, fPDG1, nullptr, nullptr) == kBMNERROR)
89 continue;
90
91 Int_t nGEM = 0;
92 Int_t nFSD = 0;
93 GetNHits(sts, nGEM, nFSD);
94
95 Float_t mom = gl->GetP();
96 if (mom > 0) { // He3/He4
97
98 // ДОБАВИТЬ КОРРЕКЦИЮ ИМПУЛЬСА ГЕЛИЯ С УЧЕТОМ ПОПРАВОК БЕТЕ-БЛОХА, КОТОРЫЕ НЕ УЧИТЫВАЕТ КАЛМАН!
99
100 Float_t dedx = gl->GetdQdNLower();
101 if (dedx < hypCut->Eval(mom))
102 continue;
103 Float_t m2 = gl->GetMass2(2);
104 if (m2 > m2CutUp->Eval(mom) || m2 < m2CutLow->Eval(mom))
105 continue;
106
107 dedxPos.push_back(dedx);
108 m2Pos700.push_back(m2);
109
110 parInVPPos.push_back(par);
111 chi2NDFPos.push_back(sts->GetChi2() / sts->GetNDF());
112 nGemHitsPos.push_back(nGEM);
113 nFsdHitsPos.push_back(nFSD);
114 } else { // pi-ROR)
115 parInVPNeg.push_back(par);
116 chi2NDFNeg.push_back(sts->GetChi2() / sts->GetNDF());
117 nGemHitsNeg.push_back(nGEM);
118 nFsdHitsNeg.push_back(nFSD);
119 }
120 }
121
122 for (size_t i = 0; i < parInVPPos.size(); ++i) {
123 for (size_t j = 0; j < parInVPNeg.size(); ++j) {
124
125 FairTrackParam VP1 = parInVPPos[i];
126 FairTrackParam VP2 = parInVPNeg[j];
127 FairTrackParam V01 = parInVPPos[i];
128 FairTrackParam V02 = parInVPNeg[j];
129
130 Float_t V0Z = FindV0ByVirtualPlanes(&V01, &V02, VPZ, 50.0);
131 if (V0Z > 30 || V0Z < 0)
132 continue;
133
134 V01 = parInVPPos[i];
135 V02 = parInVPNeg[j];
136
137 fKalman->TGeoTrackPropagate(&V01, V0Z, fPDG1, nullptr, nullptr, kTRUE);
138 fKalman->TGeoTrackPropagate(&V02, V0Z, fPDG2, nullptr, nullptr, kTRUE);
139
140 Float_t V0X1 = V01.GetX();
141 Float_t V0X2 = V02.GetX();
142 Float_t V0Y1 = V01.GetY();
143 Float_t V0Y2 = V02.GetY();
144 Float_t VPX1 = VP1.GetX();
145 Float_t VPX2 = VP2.GetX();
146 Float_t VPY1 = VP1.GetY();
147 Float_t VPY2 = VP2.GetY();
148
149 Float_t dV0X1sq = V01.GetCovariance(0, 0);
150 Float_t dV0Y1sq = V01.GetCovariance(1, 1);
151 Float_t dV0X2sq = V02.GetCovariance(0, 0);
152 Float_t dV0Y2sq = V02.GetCovariance(1, 1);
153 Float_t dVPX1sq = VP1.GetCovariance(0, 0);
154 Float_t dVPY1sq = VP1.GetCovariance(1, 1);
155 Float_t dVPX2sq = VP2.GetCovariance(0, 0);
156 Float_t dVPY2sq = VP2.GetCovariance(1, 1);
157
158 Float_t DCA12 = Sqrt(Sq(V0X2 - V0X1) + Sq(V0Y2 - V0Y1));
159 Float_t dDCA12 =
160 1.0 / DCA12 * Sqrt(Sq(V0X2 - V0X1) * (dV0X1sq + dV0X2sq) + Sq(V0Y2 - V0Y1) * (dV0Y1sq + dV0Y2sq));
161
162 Float_t DCA1 = Sqrt(Sq(VPX1 - VPX) + Sq(VPY1 - VPY));
163 Float_t dDCA1 =
164 1.0 / DCA1 * Sqrt(Sq(VPX1 - VPX) * (dVPX1sq + dVPXsq) + Sq(VPY1 - VPY) * (dVPY1sq + dVPYsq));
165 Float_t DCA2 = Sqrt(Sq(VPX2 - VPX) + Sq(VPY2 - VPY));
166 Float_t dDCA2 =
167 1.0 / DCA2 * Sqrt(Sq(VPX2 - VPX) * (dVPX2sq + dVPXsq) + Sq(VPY2 - VPY) * (dVPY2sq + dVPYsq));
168
169 Float_t V0X = .5 * (V0X1 + V0X2);
170 Float_t V0Y = .5 * (V0Y1 + V0Y2);
171 Float_t path = Sqrt(Sq(V0X - VPX) + Sq(V0Y - VPY) + Sq(V0Z - VPZ));
172
173 BmnHypNuclPair partPair;
174
175 partPair.SetM2(m2Pos700[i]);
176 partPair.SetChi2NDF(chi2NDFPos[i], chi2NDFNeg[j]);
177
178 partPair.SetDCA1(DCA1);
179 partPair.SetDCA2(DCA2);
180 partPair.SetDCA12(DCA12);
181 partPair.SetdDCA1(dDCA1);
182 partPair.SetdDCA2(dDCA2);
183 partPair.SetdDCA12(dDCA12);
184 partPair.SetPath(path);
185
186 vector<Int_t> nHits1 = {nFsdHitsPos[i], nGemHitsPos[i]};
187 vector<Int_t> nHits2 = {nFsdHitsNeg[j], nGemHitsNeg[j]};
188 partPair.SetNHitsPair(nHits1, nHits2);
189
190 Float_t p1 = Abs(1. / V01.GetQp());
191 p1 *= 2.0; // p1 is actually p1/Q, Q is 2 for helium
192 Float_t p2 = Abs(1. / V02.GetQp());
193 Float_t Tx1 = V01.GetTx();
194 Float_t Ty1 = V01.GetTy();
195 Float_t Tx2 = V02.GetTx();
196 Float_t Ty2 = V02.GetTy();
197 Float_t A1 = 1. / Sqrt(Tx1 * Tx1 + Ty1 * Ty1 + 1);
198 Float_t A2 = 1. / Sqrt(Tx2 * Tx2 + Ty2 * Ty2 + 1);
199 Float_t pz1 = p1 * A1;
200 Float_t pz2 = p2 * A2;
201 Float_t px1 = Tx1 * pz1;
202 Float_t px2 = Tx2 * pz2;
203 Float_t py1 = Ty1 * pz1;
204 Float_t py2 = Ty2 * pz2;
205
206 partPair.SetMomXYZ1(px1, py1, pz1);
207 partPair.SetMomXYZ2(px2, py2, pz2);
208
209 Float_t PzPart0 = pz1 + pz2;
210 Float_t PxPart0 = px1 + px2;
211 Float_t PyPart0 = py1 + py2;
212 Float_t TxPart0 = PxPart0 / PzPart0;
213 Float_t TyPart0 = PyPart0 / PzPart0;
214 Float_t xPart0 = TxPart0 * (VPZ - V0Z) + V0X;
215 Float_t yPart0 = TyPart0 * (VPZ - V0Z) + V0Y;
216 Float_t dca0 = Sqrt(Sq(xPart0 - VPX) + Sq(yPart0 - VPY));
217
218 partPair.SetDCA0(dca0);
219
220 lPos.SetXYZM(px1, py1, pz1, massPos);
221 lNeg.SetXYZM(px2, py2, pz2, massNeg);
222
223 partPair.SetDedx(dedxPos[i]);
224 partPair.SetInvMass(TLorentzVector((lPos + lNeg)).Mag());
225
226 partPair.SetNTrInEvent(fVertex->GetNTracks());
227
228 new ((*fParticlePair)[fParticlePair->GetEntriesFast()]) BmnHypNuclPair(partPair);
229 }
230 }
231 delete hypCut;
232 delete m2CutLow;
233 delete m2CutUp;
234}
235// -------------------------------------------------------------------
236
237// // анализ через sts треки (новы мэтчинг, пересчет длины и массы)
238// void BmnHypNuclPairFinder::Analysis()
239// {
240// TLorentzVector lPos, lNeg;
241// Float_t VPX = fVertex->GetX();
242// Float_t VPY = fVertex->GetY();
243// Float_t VPZ = fVertex->GetZ();
244// Float_t dVPXsq = fVertex->GetCovariance(0, 0);
245// Float_t dVPYsq = fVertex->GetCovariance(1, 1);
246
247// vector<FairTrackParam> parInVPPos;
248// vector<FairTrackParam> parInVPNeg;
249// vector<Float_t> m2Pos400;
250// vector<Float_t> m2Pos700;
251// vector<Float_t> chi2NDFPos;
252// vector<Float_t> chi2NDFNeg;
253// vector<Int_t> nHitsPos;
254// vector<Float_t> dedxPos;
255// vector<Int_t> nGemHitsNeg;
256// vector<Int_t> nGemHitsPos;
257// vector<Int_t> nFsdHitsNeg;
258// vector<Int_t> nFsdHitsPos;
259
260// TF1* hypCut = new TF1("hypCut", "20000*TMath::Exp(-2.0*TMath::Sqrt(x)) + 600.0", 0, 15);
261
262// // m^2 cut to select He-4 (low)
263// TF1* m2CutLow = nullptr;
264// // m^2 cut to select He-4 (up)
265// TF1* m2CutUp = nullptr;
266
267// if (massPos > 3) { // He4
268// m2CutLow = new TF1("m2CutLow", "-0.0236*x*x+0.0240*x+3.1896", 0, 30);
269// m2CutUp = new TF1("m2CutUp", "0.0169*x*x+0.1306*x+4.0525", 0, 30);
270// } else { // He3
271// m2CutLow = new TF1("m2CutLow", "-0.0034*x*x-0.0360*x+1.7426", 0, 30);
272// m2CutUp = new TF1("m2CutUp", "0.0128*x*x-0.1060*x+2.6261", 0, 30);
273// }
274
275// // создаем истинные пары между самыми близкими треками и хитами.
276// // Но только для положительных, так как pid нужна для гелия
277// map<Float_t, pair<Int_t, Int_t>> distPairs = FindPairs(fStsTracks, fTof700Hits);
278
279// for (auto const& distPair : distPairs) {
280// // Float_t dist = distPair.first;
281// pair<Int_t, Int_t> idxes = distPair.second;
282// CbmStsTrack* track = (CbmStsTrack*)fStsTracks->At(idxes.first);
283// BmnTofHit* hit = (BmnTofHit*)fTof700Hits->At(idxes.second);
284// if (track->GetFlag() == 777)
285// continue;
286// if (hit->IsUsed())
287// continue;
288
289// track->SetFlag(777);
290// hit->SetUsing(kTRUE);
291
292// if (track->GetNStsHits() < 4)
293// continue;
294
295// Float_t mom = 1.0 / track->GetParamFirst()->GetQp();
296// Float_t dedx = GetdEdx(track);
297// if (dedx < hypCut->Eval(mom))
298// continue;
299
300// Float_t m2Tof700 = -1000.0;
301
302// Double_t length1 = 0.0; // from paramFirst to vertexZ
303// FairTrackParam par = *(track->GetParamFirst());
304// if (fKalman->TGeoTrackPropagate(&par, VPZ, fPDG1, nullptr, &length1) == kBMNERROR)
305// continue;
306
307// Double_t length2 = 0.0; // from paramFirst to paramLast
308// // ADD LOOP OVER HITS TO CALCULATE CBMSTSTRACK LENGTH ====>
309// CbmStsHit* cbmhit = (CbmStsHit*)fStsHits->At(track->GetStsHitIndex(0));
310// Double_t xPrev = cbmhit->GetX();
311// Double_t yPrev = cbmhit->GetY();
312// Double_t zPrev = cbmhit->GetZ();
313// for (Int_t iHit = 1; iHit < track->GetNStsHits(); ++iHit) {
314// cbmhit = (CbmStsHit*)fStsHits->At(track->GetStsHitIndex(iHit));
315// length2 += Sqrt(Sq(cbmhit->GetX() - xPrev) + Sq(cbmhit->GetY() - yPrev) + Sq(cbmhit->GetZ() - zPrev));
316// xPrev = cbmhit->GetX();
317// yPrev = cbmhit->GetY();
318// zPrev = cbmhit->GetZ();
319// }
320// //<====
321
322// Double_t length3_701 = 0.0; // from paramLast to TOF 701 hit
323// FairTrackParam parL = *(track->GetParamLast());
324// if (fKalman->TGeoTrackPropagate(&parL, hit->GetZ(), fPDG1, nullptr, &length3_701) == kBMNERROR)
325// continue;
326// Float_t beta = (length1 + length2 + length3_701) / hit->GetTimeStamp() / (TMath::C() * 1e-7);
327// m2Tof700 = mom * mom * (1.0 / beta / beta - 1.0);
328
329// if (m2Tof700 > m2CutUp->Eval(mom) || m2Tof700 < m2CutLow->Eval(mom))
330// continue;
331
332// Int_t nGEM = 0;
333// Int_t nFSD = 0;
334// GetNHits(track, nGEM, nFSD);
335
336// parInVPPos.push_back(par);
337// m2Pos400.push_back(-1000);
338// m2Pos700.push_back(m2Tof700);
339// chi2NDFPos.push_back(track->GetChi2() / track->GetNDF());
340// dedxPos.push_back(dedx);
341// nGemHitsPos.push_back(nGEM);
342// nFsdHitsPos.push_back(nFSD);
343// }
344
345// // далее, обрабатываем отрицательные треки
346// for (Int_t i = 0; i < fStsTracks->GetEntriesFast(); i++) {
347// CbmStsTrack* track = (CbmStsTrack*)fStsTracks->UncheckedAt(i);
348
349// if (track->GetParamFirst()->GetQp() > 0.0)
350// continue;
351
352// if (track->GetNStsHits() < 4)
353// continue;
354// Int_t nGEM = 0;
355// Int_t nFSD = 0;
356// GetNHits(track, nGEM, nFSD);
357
358// FairTrackParam par = *(track->GetParamFirst());
359// if (fKalman->TGeoTrackPropagate(&par, VPZ, fPDG2, nullptr, nullptr) == kBMNERROR)
360// continue;
361// parInVPNeg.push_back(par);
362// chi2NDFNeg.push_back(track->GetChi2() / track->GetNDF());
363// nGemHitsNeg.push_back(nGEM);
364// nFsdHitsNeg.push_back(nFSD);
365// }
366
367// // cout << "size1 = " << parInVPPos.size() << " " << "size2 = " << parInVPNeg.size() << endl;
368
369// // for (FairTrackParam par1 : parInVPPos) {
370// // for (FairTrackParam par2 : parInVPNeg) {
371// for (size_t i = 0; i < parInVPPos.size(); ++i) {
372// for (size_t j = 0; j < parInVPNeg.size(); ++j) {
373
374// FairTrackParam VP1 = parInVPPos[i];
375// FairTrackParam VP2 = parInVPNeg[j];
376// FairTrackParam V01 = parInVPPos[i];
377// FairTrackParam V02 = parInVPNeg[j];
378
379// Float_t V0Z = FindV0ByVirtualPlanes(&V01, &V02, VPZ, 50.0);
380// if (V0Z > 30 || V0Z < -1)
381// continue;
382
383// V01 = parInVPPos[i];
384// V02 = parInVPNeg[j];
385
386// fKalman->TGeoTrackPropagate(&V01, V0Z, fPDG1, nullptr, nullptr, kTRUE);
387// fKalman->TGeoTrackPropagate(&V02, V0Z, fPDG2, nullptr, nullptr, kTRUE);
388
389// Float_t V0X1 = V01.GetX();
390// Float_t V0X2 = V02.GetX();
391// Float_t V0Y1 = V01.GetY();
392// Float_t V0Y2 = V02.GetY();
393// Float_t VPX1 = VP1.GetX();
394// Float_t VPX2 = VP2.GetX();
395// Float_t VPY1 = VP1.GetY();
396// Float_t VPY2 = VP2.GetY();
397
398// Float_t dV0X1sq = V01.GetCovariance(0, 0);
399// Float_t dV0Y1sq = V01.GetCovariance(1, 1);
400// Float_t dV0X2sq = V02.GetCovariance(0, 0);
401// Float_t dV0Y2sq = V02.GetCovariance(1, 1);
402// Float_t dVPX1sq = VP1.GetCovariance(0, 0);
403// Float_t dVPY1sq = VP1.GetCovariance(1, 1);
404// Float_t dVPX2sq = VP2.GetCovariance(0, 0);
405// Float_t dVPY2sq = VP2.GetCovariance(1, 1);
406
407// Float_t DCA12 = Sqrt(Sq(V0X2 - V0X1) + Sq(V0Y2 - V0Y1));
408// Float_t dDCA12 =
409// 1.0 / DCA12 * Sqrt(Sq(V0X2 - V0X1) * (dV0X1sq + dV0X2sq) + Sq(V0Y2 - V0Y1) * (dV0Y1sq + dV0Y2sq));
410
411// Float_t DCA1 = Sqrt(Sq(VPX1 - VPX) + Sq(VPY1 - VPY));
412// Float_t dDCA1 =
413// 1.0 / DCA1 * Sqrt(Sq(VPX1 - VPX) * (dVPX1sq + dVPXsq) + Sq(VPY1 - VPY) * (dVPY1sq + dVPYsq));
414// Float_t DCA2 = Sqrt(Sq(VPX2 - VPX) + Sq(VPY2 - VPY));
415// Float_t dDCA2 =
416// 1.0 / DCA2 * Sqrt(Sq(VPX2 - VPX) * (dVPX2sq + dVPXsq) + Sq(VPY2 - VPY) * (dVPY2sq + dVPYsq));
417
418// Float_t V0X = .5 * (V0X1 + V0X2);
419// Float_t V0Y = .5 * (V0Y1 + V0Y2);
420// Float_t path = Sqrt(Sq(V0X - VPX) + Sq(V0Y - VPY) + Sq(V0Z - VPZ));
421
422// BmnHypNuclPair partPair;
423
424// // put mass2 for positive particle in field beta (let it be)
425// // for nevative one put -100, we don't need PID for pion
426// partPair.SetM2(m2Pos700[i]);
427// partPair.SetChi2NDF(chi2NDFPos[i], chi2NDFNeg[j]);
428
429// partPair.SetDCA1(DCA1);
430// partPair.SetDCA2(DCA2);
431// partPair.SetDCA12(DCA12);
432// partPair.SetdDCA1(dDCA1);
433// partPair.SetdDCA2(dDCA2);
434// partPair.SetdDCA12(dDCA12);
435// partPair.SetPath(path);
436
437// vector<Int_t> nHits1 = {nFsdHitsPos[i], nGemHitsPos[i]};
438// vector<Int_t> nHits2 = {nFsdHitsNeg[j], nGemHitsNeg[j]};
439// partPair.SetNHitsPair(nHits1, nHits2);
440
441// Float_t p1 = Abs(1. / V01.GetQp());
442// p1 *= 2.0; // p1 is actually p1/Q, Q is 2 for helium
443// Float_t p2 = Abs(1. / V02.GetQp());
444// Float_t Tx1 = V01.GetTx();
445// Float_t Ty1 = V01.GetTy();
446// Float_t Tx2 = V02.GetTx();
447// Float_t Ty2 = V02.GetTy();
448// Float_t A1 = 1. / Sqrt(Tx1 * Tx1 + Ty1 * Ty1 + 1);
449// Float_t A2 = 1. / Sqrt(Tx2 * Tx2 + Ty2 * Ty2 + 1);
450// Float_t pz1 = p1 * A1;
451// Float_t pz2 = p2 * A2;
452// Float_t px1 = Tx1 * pz1;
453// Float_t px2 = Tx2 * pz2;
454// Float_t py1 = Ty1 * pz1;
455// Float_t py2 = Ty2 * pz2;
456
457// partPair.SetMomXYZ1(px1, py1, pz1);
458// partPair.SetMomXYZ2(px2, py2, pz2);
459
460// Float_t PzPart0 = pz1 + pz2;
461// Float_t PxPart0 = px1 + px2;
462// Float_t PyPart0 = py1 + py2;
463// Float_t TxPart0 = PxPart0 / PzPart0;
464// Float_t TyPart0 = PyPart0 / PzPart0;
465// Float_t xPart0 = TxPart0 * (VPZ - V0Z) + V0X;
466// Float_t yPart0 = TyPart0 * (VPZ - V0Z) + V0Y;
467// Float_t dca0 = Sqrt(Sq(xPart0 - VPX) + Sq(yPart0 - VPY));
468
469// partPair.SetDCA0(dca0);
470
471// // partPair.SetEtaPair(-100, -100/*_eta1, _eta2*/);
472
473// lPos.SetXYZM(px1, py1, pz1, massPos);
474// lNeg.SetXYZM(px2, py2, pz2, massNeg);
475
476// partPair.SetDedx1(dedxPos[i]);
477// partPair.SetInvMass(TLorentzVector((lPos + lNeg)).Mag());
478
479// // partPair.SetTxPair(Tx1, Tx2);
480// // partPair.SetTyPair(Ty1, Ty2);
481// partPair.SetMass1(massPos);
482// partPair.SetMass2(massNeg);
483
484// partPair.SetNTrInEvent(fVertex->GetNTracks());
485
486// new ((*fParticlePair)[fParticlePair->GetEntriesFast()]) BmnHypNuclPair(partPair);
487// }
488// }
489// delete hypCut;
490// delete m2CutLow;
491// delete m2CutUp;
492// }
493// -------------------------------------------------------------------
494
495InitStatus BmnHypNuclPairFinder::Init()
496{
497 cout << "\nBmnHypNuclPairFinder::Init()" << endl;
498
499 FairRootManager* ioman = FairRootManager::Instance();
500
501 fStsTracks = (TClonesArray*)ioman->GetObject("StsVector");
502 fGlobTracks = (TClonesArray*)ioman->GetObject("BmnGlobalTrack");
503 fStsHits = (TClonesArray*)ioman->GetObject("StsHit");
504 // fVertex = (CbmVertex*)ioman->GetObject("PrimaryVertex.");
505 fVertex = (CbmVertex*)ioman->GetObject("MpdVertex.");
506 fTrigger = (BmnTrigInfoDst*)ioman->GetObject("BmnTrigInfo.");
507 fUpperGemClusters = (TClonesArray*)ioman->GetObject("BmnGemUpperCluster");
508 fLowerGemClusters = (TClonesArray*)ioman->GetObject("BmnGemLowerCluster");
509 fUpperFsdClusters = (TClonesArray*)ioman->GetObject("BmnSiliconUpperCluster");
510 fLowerFsdClusters = (TClonesArray*)ioman->GetObject("BmnSiliconLowerCluster");
511 fTof700Hits = (TClonesArray*)ioman->GetObject("BmnTof700Hit");
512 // fTof700Hits = (TClonesArray*)ioman->GetObject("BmnTof701Hit");
513
514 fParticlePair = new TClonesArray("BmnHypNuclPair");
515 ioman->Register("pairs", "HYP", fParticlePair, kTRUE);
516
517 fKalman = new BmnKalmanFilter();
518
519 return kSUCCESS;
520}
521
522// -------------------------------------------------------------------
523
524void BmnHypNuclPairFinder::Exec(Option_t* option)
525{
526 fParticlePair->Delete();
527
528 fEventCounter++;
529
530 // if (fEventCounter % 100 == 0)
531 cout << "Event# " << fEventCounter << endl;
532
533 if (isExp) {
534 if (!(fTrigger->IsTriggerBitTrueAR(7))) // bit #7 is for CCT2
535 return;
536 if (fVertex->GetNTracks() < 2)
537 return;
538 if (fVertex->GetZ() < -0.5 || fVertex->GetZ() > +0.5)
539 return;
540
541 Float_t centerOfTargetX = 0.40;
542 Float_t centerOfTargetY = 0.15;
543 Float_t radiusOfTarget = 1.2;
544
545 if (TMath::Sqrt(TMath::Sq(fVertex->GetX() - centerOfTargetX) + TMath::Sq(fVertex->GetY() - centerOfTargetY))
546 > radiusOfTarget)
547 return;
548 }
549 Analysis();
550}
551// -------------------------------------------------------------------
552
553void BmnHypNuclPairFinder::Finish()
554{
555 delete fKalman;
556 // delete fMagField;
557 cout << "\n-I- [BmnHypNuclPairFinder::Finish] " << endl;
558}
559
560Float_t BmnHypNuclPairFinder::FindV0ByVirtualPlanes(FairTrackParam* par1,
561 FairTrackParam* par2,
562 Float_t z_0,
563 Float_t range)
564{
565 const Int_t nPlanes = 5; // FIXME
566
567 while (range >= 0.1) {
568 Float_t zMax = z_0 + range;
569 Float_t zMin = z_0 - range;
570 Float_t zStep = (zMax - zMin) / nPlanes;
571
572 Float_t minDist = 1000.0;
573 Float_t minZ = -1000.0;
574
575 for (Int_t iPlane = 0; iPlane < nPlanes; ++iPlane) {
576 Float_t z = zMax - iPlane * zStep;
577 BmnStatus stat1 = fKalman->TGeoTrackPropagate(par1, z, fPDG1, nullptr, nullptr, kTRUE);
578 BmnStatus stat2 = fKalman->TGeoTrackPropagate(par2, z, fPDG2, nullptr, nullptr, kTRUE);
579
580 if (stat1 == kBMNERROR || stat2 == kBMNERROR)
581 return DBL_MAX;
582
583 Float_t d = Sqrt(Sq(par1->GetX() - par2->GetX()) + Sq(par1->GetY() - par2->GetY()));
584 if (d < minDist) {
585 minDist = d;
586 minZ = z;
587 }
588 }
589
590 z_0 = minZ;
591 range *= 0.5;
592 }
593 return z_0;
594}
595
596Float_t BmnHypNuclPairFinder::GetdEdx(CbmStsTrack* tr)
597{
598 if (isExp) {
599 // set<Float_t> signalsUpp;
600 set<Float_t> signalsLow;
601 Int_t nHits = tr->GetNStsHits();
602 Int_t nGemHits = 0;
603 for (Int_t hitIdx = 0; hitIdx < nHits; ++hitIdx) {
604 CbmStsHit* hit = (CbmStsHit*)fStsHits->At(tr->GetStsHitIndex(hitIdx));
605
606 if (hit->GetSystemId() != kGEM)
607 continue;
608
609 nGemHits++;
610
611 // StripCluster* upp = (StripCluster*)fUpperGemClusters->At(hit->fDigiB % 1000000);
612 StripCluster* low = (StripCluster*)fLowerGemClusters->At(hit->fDigiF % 1000000);
613
614 // signalsUpp.insert(upp->TotalSignal * aUpp[hit->GetStationNr() - 1] + bUpp[hit->GetStationNr() - 1]);
615 signalsLow.insert(low->TotalSignal * aLow[hit->GetStationNr() - 1] + bLow[hit->GetStationNr() - 1]);
616 }
617
618 if (nGemHits < 3)
619 return 0.0;
620
621 Int_t usedHits = (nGemHits == 3) ? 2
622 : (nGemHits == 4) ? 2
623 : (nGemHits == 5) ? 3
624 : (nGemHits == 6) ? 4
625 : (nGemHits == 7) ? 4
626 : 0;
627
628 // auto itUpp = signalsUpp.begin();
629 auto itLow = signalsLow.begin();
630
631 Float_t totSigLow = 0.0;
632 // Float_t totSigUpp = 0.0;
633 for (Int_t i = 0; i < usedHits; ++i) {
634 // totSigUpp += (*itUpp);
635 totSigLow += (*itLow);
636 // itUpp++;
637 itLow++;
638 }
639
640 totSigLow /= usedHits;
641 // totSigUpp /= usedHits;
642 // return 0.5 * (totSigLow + totSigUpp);
643 return totSigLow;
644 } else {
645 return 10000;
646 }
647}
648
649void BmnHypNuclPairFinder::GetNHits(CbmStsTrack* tr, Int_t& nGEM, Int_t& nFSD)
650{
651 nGEM = 0;
652 nFSD = 0;
653 for (Int_t hitIdx = 0; hitIdx < tr->GetNStsHits(); ++hitIdx) {
654 CbmStsHit* hit = (CbmStsHit*)fStsHits->At(tr->GetStsHitIndex(hitIdx));
655
656 if (hit->GetSystemId() == kGEM)
657 nGEM++;
658 else if (hit->GetSystemId() == kSILICON)
659 nFSD++;
660 }
661}
662
663Float_t BmnHypNuclPairFinder::GetDyTof700(FairTrackParam* par)
664{
665 if (isExp) {
666 Float_t mom = 1.0 / par->GetQp();
667 return (mom > 0) ? 0.17 * mom - 0.40 : 0.0;
668 } else
669 return 0.0;
670}
671
672Float_t BmnHypNuclPairFinder::GetDxTof700(FairTrackParam* par)
673{
674 if (isExp) {
675 return 0.0;
676 } else
677 return 0.0;
678}
679
680Float_t BmnHypNuclPairFinder::GetSigXTof700He3(FairTrackParam* par)
681{
682 if (isExp) {
683 Float_t mom = 1.0 / par->GetQp();
684 return (mom > 0) ? 6.58 * TMath::Exp(-1.73 * mom) + 0.69 : 1.0;
685 } else
686 return 1.0;
687}
688
689Float_t BmnHypNuclPairFinder::GetSigYTof700He3(FairTrackParam* par)
690{
691 if (isExp) {
692 Float_t mom = 1.0 / par->GetQp();
693 return (mom > 0) ? 7.52 * TMath::Exp(-2.08 * mom) + 0.53 : 1.0;
694 } else
695 return 1.0;
696}
697
698map<Float_t, pair<Int_t, Int_t>> BmnHypNuclPairFinder::FindPairs(TClonesArray* tracks, TClonesArray* hits)
699{
700
701 map<Float_t, pair<Int_t, Int_t>> distancesPairs;
702
703 // Double_t minZ = 10000.0;
704 // for (Int_t hitIdx = 0; hitIdx < hits->GetEntriesFast(); ++hitIdx) {
705 // BmnHit* hit = (BmnHit*)hits->At(hitIdx);
706 // if (hit->GetZ() < minZ)
707 // minZ = hit->GetZ();
708 // }
709
710 Float_t zTof701HitsMin = 616; // median z position for tof701
711
712 for (Int_t trIdx = 0; trIdx < tracks->GetEntriesFast(); ++trIdx) {
713 CbmStsTrack* tr = (CbmStsTrack*)tracks->At(trIdx);
714 if (tr->GetParamFirst()->GetQp() < 0.0)
715 continue;
716 FairTrackParam parMinZ(*(tr->GetParamLast()));
717 if (fKalman->TGeoTrackPropagate(&parMinZ, zTof701HitsMin, fPDG1, nullptr, nullptr) == kBMNERROR)
718 continue;
719 for (Int_t hitIdx = 0; hitIdx < hits->GetEntriesFast(); ++hitIdx) {
720 BmnHit* hit = (BmnHit*)hits->At(hitIdx);
721 FairTrackParam param = parMinZ;
722 if (fKalman->TGeoTrackPropagate(&param, hit->GetZ(), fPDG1, nullptr, nullptr) == kBMNERROR)
723 continue;
724 Float_t dX = Abs(param.GetX() - GetDxTof700(&param) - hit->GetX());
725 Float_t dY = Abs(param.GetY() - GetDyTof700(&param) - hit->GetY());
726 Float_t xCut = 3.0 * GetSigXTof700He3(&param);
727 Float_t yCut = 3.0 * GetSigYTof700He3(&param);
728 if (dX < xCut && dY < yCut) {
729 Float_t dist = Sqrt(dX * dX + dY * dY);
730 hit->SetUsing(kFALSE);
731 distancesPairs.insert(pair<Float_t, pair<Int_t, Int_t>>(dist, pair<Int_t, Int_t>(trIdx, hitIdx)));
732 }
733 }
734 }
735
736 return distancesPairs;
737}
738
739ClassImp(BmnHypNuclPairFinder);
ClassImp
ClassImp(BmnHypNuclPairFinder)
dist
vector< Double_t > dist(vector< Double_t > qp, Double_t mu)
Definition BmnMath.cxx:869
d
const Float_t d
Z-ccordinate of the first GEM-station.
Definition BmnMwpcHit.cxx:7
i
int i
Definition P4_F32vec4.h:22
exp
friend F32vec4 exp(const F32vec4 &a)
Definition P4_F32vec4.h:122
kSILICON
@ kSILICON
Definition BmnDetectorList.h:25
kGEM
@ kGEM
Definition BmnDetectorList.h:15
BmnStatus
BmnStatus
Definition BmnEnums.h:24
kBMNERROR
@ kBMNERROR
Definition BmnEnums.h:26
BmnGlobalTrack::GetMass2
Double_t GetMass2(Int_t tofID)
Definition BmnGlobalTrack.cxx:60
BmnGlobalTrack::GetGemTrackIndex
Int_t GetGemTrackIndex() const
Definition BmnGlobalTrack.h:36
BmnGlobalTrack::GetdQdNLower
Double_t GetdQdNLower() const
Definition BmnGlobalTrack.h:53
BmnHit::SetUsing
void SetUsing(Bool_t use)
Definition BmnHit.h:53
BmnHypNuclPairFinder::Init
virtual InitStatus Init()
Definition BmnHypNuclPairFinder.cxx:495
BmnHypNuclPairFinder::Exec
virtual void Exec(Option_t *option)
Definition BmnHypNuclPairFinder.cxx:524
BmnHypNuclPair::SetdDCA12
void SetdDCA12(Float_t val)
Definition BmnHypNuclPair.h:32
BmnHypNuclPair::SetDCA0
void SetDCA0(Float_t val)
Definition BmnHypNuclPair.h:20
BmnHypNuclPair::SetM2
void SetM2(Float_t val)
Definition BmnHypNuclPair.h:58
BmnHypNuclPair::SetDCA1
void SetDCA1(Float_t val)
Definition BmnHypNuclPair.h:22
BmnHypNuclPair::SetInvMass
void SetInvMass(Float_t val)
Definition BmnHypNuclPair.h:36
BmnHypNuclPair::SetDCA2
void SetDCA2(Float_t val)
Definition BmnHypNuclPair.h:24
BmnHypNuclPair::SetChi2NDF
void SetChi2NDF(Float_t val1, Float_t val2)
Definition BmnHypNuclPair.h:52
BmnHypNuclPair::SetNHitsPair
void SetNHitsPair(vector< Int_t > part1, vector< Int_t > part2)
Definition BmnHypNuclPair.h:60
BmnHypNuclPair::SetDCA12
void SetDCA12(Float_t val)
Definition BmnHypNuclPair.h:26
BmnHypNuclPair::SetdDCA1
void SetdDCA1(Float_t val)
Definition BmnHypNuclPair.h:28
BmnHypNuclPair::SetNTrInEvent
void SetNTrInEvent(Int_t n)
Definition BmnHypNuclPair.h:78
BmnHypNuclPair::SetPath
void SetPath(Float_t val)
Definition BmnHypNuclPair.h:34
BmnHypNuclPair::SetDedx
void SetDedx(Float_t val)
Definition BmnHypNuclPair.h:160
BmnHypNuclPair::SetdDCA2
void SetdDCA2(Float_t val)
Definition BmnHypNuclPair.h:30
BmnHypNuclPair::SetMomXYZ2
void SetMomXYZ2(Float_t px, Float_t py, Float_t pz)
Definition BmnHypNuclPair.h:45
BmnHypNuclPair::SetMomXYZ1
void SetMomXYZ1(Float_t px, Float_t py, Float_t pz)
Definition BmnHypNuclPair.h:38
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
BmnTrack::GetP
Double_t GetP()
Definition BmnTrack.h:80
BmnTrigInfoDst::IsTriggerBitTrueAR
Bool_t IsTriggerBitTrueAR(Short_t bit)
Definition BmnTrigInfoDst.h:71
CbmStsHit::GetStationNr
virtual Int_t GetStationNr() const
Definition CbmStsHit.h:66
CbmStsHit::GetSystemId
Int_t GetSystemId() const
Definition CbmStsHit.h:64
CbmStsHit::fDigiF
Int_t fDigiF
Definition CbmStsHit.h:92
CbmStsTrack::FromVertex
int FromVertex() const
Definition CbmStsTrack.h:76
CbmStsTrack::GetParamLast
FairTrackParam * GetParamLast()
Definition CbmStsTrack.h:70
CbmStsTrack::GetNStsHits
Int_t GetNStsHits() const
Definition CbmStsTrack.h:60
CbmStsTrack::GetChi2
Double_t GetChi2() const
Definition CbmStsTrack.h:66
CbmStsTrack::GetStsHitIndex
Int_t GetStsHitIndex(Int_t iHit) const
Definition CbmStsTrack.h:62
CbmStsTrack::GetParamFirst
FairTrackParam * GetParamFirst()
Definition CbmStsTrack.h:69
CbmStsTrack::GetNDF
Int_t GetNDF() const
Definition CbmStsTrack.h:67
CbmVertex::GetZ
Double_t GetZ() const
Definition CbmVertex.h:60
CbmVertex::GetX
Double_t GetX() const
Definition CbmVertex.h:58
CbmVertex::GetNTracks
Int_t GetNTracks() const
Definition CbmVertex.h:63
CbmVertex::GetCovariance
Double_t GetCovariance(Int_t i, Int_t j) const
Definition CbmVertex.cxx:115
CbmVertex::GetY
Double_t GetY() const
Definition CbmVertex.h:59
StripCluster::TotalSignal
Double_t TotalSignal
Definition BmnStripData.h:181