BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
BmnUpstreamTracking.cxx
Go to the documentation of this file.
1// Author: Vasilisa Lenivenko (VBLHEP) <vasilisa@jinr.ru> 2019-07-18
2
3////////////////////////////////////////////////////////////////////////////////
4// //
5// BmnUpstreamTracking //
6// //
7// Implementation of an algorithm developed by //
8// V.Lenivenko and V.Paltchik //
9// to the BmnRoot software //
10// //
11// The algorithm serves for searching for tracks //
12// in the Silicon detectors & MWPCs //
13// of the configuration SRC of BM@N experiment //
14// //
15////////////////////////////////////////////////////////////////////////////////
16
17#include "BmnUpstreamTracking.h"
18
19#include "FairLogger.h"
20#include "FairRootManager.h"
21#include "TFile.h"
22
23#include <TMath.h>
24#include <vector>
25
26static Double_t workTime = 0.0;
27TList fhList;
28
29//----------------------------------------------------------------------
30void BmnUpstreamTracking::Exec(Option_t* opt)
31{
32 if (!IsActive())
33 return;
34
35 clock_t tStart = clock();
36 PrepareArraysToProcessEvent();
37 if (fDebug)
38 cout << "\n======================== Upstream track finder exec started ===================\n" << endl;
39 if (fDebug)
40 printf("Event number: %d\n", fEventNo);
41 fEventNo++;
42
43 if (fDebug)
44 cout << "--------------Read-------------" << endl;
45 ReadSiliconTracks(par_ab_SiTr, par_SiTr_z, NSiTracks, vec_points);
46 ReadSiliconHits(SiXYhits, NSiXYhits);
47
48 if (fDebug)
49 cout << "--------------Reco: MWPCSegments-----" << endl;
50 ReadMWPCSegments(par_ab_Ch, par_Seg_z, Nseg_Ch, vec_points);
51 if (fDebug)
52 cout << "--------------Reco: MWPCTracks" << endl;
53 ReadMWPCTracks(par_ab_tr, par_ab_trz, NPCTracks);
54
55 //-----------------Matching-------------------------------------------
56 if (fDebug)
57 cout << "--------------Matching Si-tracks & PC-tracks-----" << endl;
58 PCTracksSiTracksMatching(par_ab_SiTr, par_SiTr_z, NSiTracks, par_ab_tr, par_ab_trz, NPCTracks, vtmpSeg,
59 OutVector); // Step1
60 if (fDebug)
61 cout << "--------------Matching Si-tracks & PC-segments---" << endl;
62 PCSegmentsSiTracksMatching(par_ab_SiTr, par_SiTr_z, NSiTracks, par_ab_Ch, par_Seg_z, Nseg_Ch, vtmpSeg,
63 OutVector); // Step2
64 if (fDebug)
65 cout << "--------------Adding hits-------------------------" << endl;
66 GetAddSiHits(OutVector, SiXYhits, NSiXYhits);
67 GetHitsToPoints(OutVector, par_ab_tr, NPCTracks, par_ab_Ch, Nseg_Ch, Points, NumPoints);
68
69 if (fDebug)
70 cout << "--------------Fit--------------------------------" << endl;
71 CalculationOfChi2(Points, NumPoints, vecUpTracks);
72
73 if (fDebug)
74 PrintAllTracks(vecUpTracks);
75 if (!expData)
76 MCefficiencyCalculation(vec_points, vecUpTracks);
77
78 if (fDebug)
79 cout << "----TrackRecording--------------------------------" << endl;
80 TrackRecording(vecUpTracks);
81
82 if (fDebug)
83 cout << endl;
84 if (fDebug)
85 cout << "======================== Upstream track finder exec finished ===================" << endl;
86 clock_t tFinish = clock();
87 workTime += ((Double_t)(tFinish - tStart)) / CLOCKS_PER_SEC;
88}
89//----------------------------------------------------------------------
90
91void BmnUpstreamTracking::MCefficiencyCalculation(vector<MC_points>& points, vector<UpTracks>& vecUp)
92{
93
94 if (fDebug)
95 cout << " ---UpMC tracks association--" << endl;
96 // ax, bx, ay, by
97 Double_t delta2[4] = {-999., -999., -999., -999.};
98 Double_t sig[4] = {0.03, 4., 0.03, 4.};
99
100 Double_t dmatch = 0.;
101 Double_t dmc_match[kMaxMC];
102 Int_t mc_tr_assoc[kMaxMC];
103
104 for (Int_t j = 0; j < kMaxMC; j++) {
105 dmc_match[j] = 1000.;
106 mc_tr_assoc[j] = -1;
107 }
108
109 Double_t dax = -999.;
110 Double_t day = -999.;
111 Double_t dx = -999.;
112 Double_t dy = -999.;
113
114 if (fDebug)
115 cout << " MC " << points.size() << " UpTracks " << vecUp.size() << endl;
116
117 Int_t Ngood_mc_tracks = 0;
118 Int_t Ngood_reco_tracks = 0;
119 for (size_t itr = 0; itr < points.size(); itr++) { // mc_tr
120
121 //---MC Eff ---
122 //---Den
123 if (fDebug)
124 cout << " Np " << points.at(itr).Np << " pdg " << points.at(itr).Pdg << endl;
125 if (fDebug && points.at(itr).Np >= 7 //>=10
126 && points.at(itr).np_3si == 0 && (points.at(itr).Np_ch2 > 3 || points.at(itr).Np_ch3 > 3))
127 {
128 // if ( points.at(itr).xWas3 && points.at(itr).uWas3 && points.at(itr).vWas3 &&
129 // points.at(itr).xWas2 && points.at(itr).uWas2 && points.at(itr).vWas2 )
130 hDen_mcuptr->Fill(1);
131 hDen_mcuptr->Fill(0);
132 hY_vs_Xmctrue->Fill(points.at(itr).param[1], points.at(itr).param[3]);
133
134 Ngood_mc_tracks++;
135 cout << " itr " << itr << " Ngood_mc_tracks " << Ngood_mc_tracks << endl;
136 cout << "UpDen" << endl;
137 }
138
139 for (size_t iup = 0; iup < vecUp.size(); iup++) { // reco
140 delta2[0] = points.at(itr).param[0] - vecUp.at(iup).param[0];
141 delta2[1] = points.at(itr).param[1] - vecUp.at(iup).param[1];
142 delta2[2] = points.at(itr).param[2] - vecUp.at(iup).param[2];
143 delta2[3] = points.at(itr).param[3] - vecUp.at(iup).param[3];
144
145 dmatch = 0.;
146 dmatch = (delta2[0] * delta2[0]) / (sig[0] * sig[0]) + (delta2[1] * delta2[1]) / (sig[1] * sig[1])
147 + (delta2[2] * delta2[2]) / (sig[2] * sig[2]) + (delta2[3] * delta2[3]) / (sig[3] * sig[3]);
148
149 // combinatorics
150 if (fDebug && delta2[0] > -900.)
151 hdAx_tr_mc_comb->Fill(delta2[0]);
152 if (fDebug && delta2[1] > -900.)
153 hdX_tr_mc_comb->Fill(delta2[1]);
154 if (fDebug && delta2[2] > -900.)
155 hdAy_tr_mc_comb->Fill(delta2[2]);
156 if (fDebug && delta2[3] > -900.)
157 hdY_tr_mc_comb->Fill(delta2[3]);
158
159 if (fDebug)
160 cout << " imc " << itr << " iup " << iup << " Nhits " << vecUp.at(iup).Nhits << " ax_mc "
161 << points.at(itr).param[0] << " ay_mc " << points.at(itr).param[2] << " ax_up "
162 << vecUp.at(iup).param[0] << " ay_up " << vecUp.at(iup).param[2] << " dmc_match " << dmc_match[itr]
163 << " " << dmatch << " delta2[0] " << delta2[0] << " delta2[1] " << delta2[1] << " delta2[2] "
164 << delta2[2] << " delta2[3] " << delta2[3] << endl;
165
166 if (dmc_match[itr] > dmatch) {
167 dmc_match[itr] = dmatch;
168 mc_tr_assoc[itr] = iup;
169 // vecUp.at(iup).Pdg = points.at(itr).Pdg;
170 dax = delta2[0];
171 dx = delta2[1];
172 day = delta2[2];
173 dy = delta2[3];
174 }
175
176 } // vecUp.size()
177 if (fDebug) {
178 // cout<<" itr "<<itr<<" mc_Id "<<points.at(itr).Id<<
179 // " ax_mc "<<points.at(itr).paramSi[0]<<" reco_ind "<<mc_tr_assoc[itr]<<" ax "<<
180 // vecUp.at(mc_tr_assoc[itr]).param[0] << " dmc_match "<<dmc_match[itr]<<endl;
181 if (mc_tr_assoc[itr] > -1) {
182 cout << " mc itr " << itr << " rec " << mc_tr_assoc[itr] << endl;
183 if (dax > -900.)
184 hdAx_uptr_mc->Fill(dax);
185 if (dx > -900.)
186 hdX_uptr_mc->Fill(dx);
187 if (day > -900.)
188 hdAy_uptr_mc->Fill(day);
189 if (dy > -900.)
190 hdY_uptr_mc->Fill(dy);
191
192 } // if mc_tr_assoc[itr] > -1
193 } // if
194 } // points.size
195
196 if (fDebug)
197 cout << "reject poorly chosen association segments " << endl;
198 for (size_t itr = 0; itr < points.size(); itr++) { // mc_tr
199 if (mc_tr_assoc[itr] == -1)
200 continue;
201
202 for (size_t itr2 = 0; itr2 < points.size(); itr2++) { // mc_tr
203 if (itr2 == itr)
204 continue;
205 if (mc_tr_assoc[itr2] == -1)
206 continue;
207
208 if (mc_tr_assoc[itr] == mc_tr_assoc[itr2]) {
209 if (dmc_match[itr2] > dmc_match[itr])
210 mc_tr_assoc[itr2] = -1;
211 else {
212 mc_tr_assoc[itr] = -1;
213 break;
214 }
215 }
216 } // itr2
217 //---MC Eff ---
218 //---Num
219 if (fDebug)
220 cout << " mc_Id " << points.at(itr).Id << " pdg " << points.at(itr).Pdg << " assoc " << mc_tr_assoc[itr]
221 << endl;
222 if (fDebug && mc_tr_assoc[itr] > -1 && fDebug && points.at(itr).Np >= 7 //>= 10
223 && points.at(itr).np_3si == 0 && (points.at(itr).Np_ch2 > 3 || points.at(itr).Np_ch3 > 3))
224 {
225 if (points.at(itr).xWas3 && points.at(itr).uWas3 && points.at(itr).vWas3 && points.at(itr).xWas2
226 && points.at(itr).uWas2 && points.at(itr).vWas2)
227 hNum_mcuptr->Fill(1);
228
229 hNum_mcuptr->Fill(0);
230 Ngood_reco_tracks++;
231 // cout<<" itr "<<itr<<" Ngood_reco_tracks "<<Ngood_reco_tracks<<endl;
232 cout << "UpNum" << endl;
233 }
234 } // itr
235 if (fDebug)
236 cout << " Ngood_mc_tracks " << Ngood_mc_tracks << " Ngood_reco_tracks " << Ngood_reco_tracks << endl;
237 if (fDebug) {
238 hNtr_reco->Fill(Ngood_reco_tracks);
239 hNtr_mc->Fill(Ngood_mc_tracks);
240 hNtr_mc_vs_reco->Fill(Ngood_reco_tracks, Ngood_mc_tracks);
241 if (Ngood_mc_tracks == 2)
242 hNrecoTrif2mc->Fill(Ngood_reco_tracks);
243 }
244 // if (fDebug && points.size() == 2) hNrecoTrif2mc->Fill(vecUp.size());
245 for (size_t itr = 0; itr < points.size(); itr++) { // mc_tr
246 if (mc_tr_assoc[itr] > -1)
247 vecUp.at(mc_tr_assoc[itr]).Pdg = points.at(itr).Pdg;
248 }
249
250 twotracks tt;
251 Bool_t track1_was = 0;
252 Bool_t ifNmc2 = 0;
253 Bool_t ifSelectedCases = 0;
254 if (fDebug)
255 cout << " reco tracks: " << endl;
256
257 if (Ngood_mc_tracks == 2)
258 ifNmc2 = 1;
259
260 if (vecUp.size() == 2 && Ngood_mc_tracks == 2) {
261 for (size_t iup = 0; iup < vecUp.size(); iup++) { // reco
262 if (fDebug)
263 cout << " iup " << iup << " found code pdg " << vecUp.at(iup).Pdg << endl;
264 if (track1_was == 0) {
265 track1_was = 1;
266 tt.track1_pdg = vecUp.at(0).Pdg;
267 tt.track1[0] = vecUp.at(iup).param[0];
268 tt.track1[1] = vecUp.at(iup).param[1];
269 tt.track1[2] = vecUp.at(iup).param[2];
270 tt.track1[3] = vecUp.at(iup).param[3];
271 }
272 if (track1_was == 1) {
273 tt.track2_pdg = vecUp.at(1).Pdg;
274 tt.track2[0] = vecUp.at(iup).param[0];
275 tt.track2[1] = vecUp.at(iup).param[1];
276 tt.track2[2] = vecUp.at(iup).param[2];
277 tt.track2[3] = vecUp.at(iup).param[3];
278 }
279 } // reco
280 if (track1_was == 1)
281 vec_twotracks.push_back(tt);
282 if (fDebug)
283 cout << " vec_twotracks " << vec_twotracks.size() << endl;
284
285 // for(int itr = 0; itr< vec_twotracks.size(); ++itr){
286 if (vec_twotracks.size() == 1) {
287 Int_t itr = 0;
288 if ((vec_twotracks.at(itr).track1_pdg == PDG_Be7 && vec_twotracks.at(itr).track2_pdg == PDG_p)
289 || (vec_twotracks.at(itr).track2_pdg == PDG_Be7 && vec_twotracks.at(itr).track1_pdg == PDG_p)
290 || (vec_twotracks.at(itr).track1_pdg == PDG_Be7 && vec_twotracks.at(itr).track2_pdg == PDG_H2)
291 || (vec_twotracks.at(itr).track2_pdg == PDG_Be7 && vec_twotracks.at(itr).track1_pdg == PDG_H2)
292 || (vec_twotracks.at(itr).track1_pdg == PDG_Be7 && vec_twotracks.at(itr).track2_pdg == PDG_He4)
293 || (vec_twotracks.at(itr).track2_pdg == PDG_Be7 && vec_twotracks.at(itr).track1_pdg == PDG_He4)
294 || (vec_twotracks.at(itr).track1_pdg == PDG_Be7 && vec_twotracks.at(itr).track2_pdg == PDG_He3)
295 || (vec_twotracks.at(itr).track2_pdg == PDG_Be7 && vec_twotracks.at(itr).track1_pdg == PDG_He3)
296 || (vec_twotracks.at(itr).track1_pdg == PDG_Li7 && vec_twotracks.at(itr).track2_pdg == PDG_He4)
297 || (vec_twotracks.at(itr).track2_pdg == PDG_Li7 && vec_twotracks.at(itr).track1_pdg == PDG_He4)
298 || (vec_twotracks.at(itr).track1_pdg == PDG_Li7 && vec_twotracks.at(itr).track2_pdg == PDG_He3)
299 || (vec_twotracks.at(itr).track2_pdg == PDG_Li7 && vec_twotracks.at(itr).track1_pdg == PDG_He3)
300 || (vec_twotracks.at(itr).track1_pdg == PDG_Li6 && vec_twotracks.at(itr).track2_pdg == PDG_He4)
301 || (vec_twotracks.at(itr).track2_pdg == PDG_Li6 && vec_twotracks.at(itr).track1_pdg == PDG_He4)
302 || (vec_twotracks.at(itr).track1_pdg == PDG_Li6 && vec_twotracks.at(itr).track2_pdg == PDG_He3)
303 || (vec_twotracks.at(itr).track2_pdg == PDG_Li6 && vec_twotracks.at(itr).track1_pdg == PDG_He3)
304 || (vec_twotracks.at(itr).track1_pdg == PDG_He4 && vec_twotracks.at(itr).track2_pdg == PDG_He4)
305
306 )
307 ifSelectedCases = 1;
308
309 if (fDebug)
310 cout << " pdg1 " << vec_twotracks.at(itr).track1_pdg << " pdg2 " << vec_twotracks.at(itr).track2_pdg
311 << " ifSelectedCases " << ifSelectedCases << endl;
312 Double_t CosA = (vec_twotracks.at(itr).track1[0] * vec_twotracks.at(itr).track2[0]
313 + vec_twotracks.at(itr).track1[2] * vec_twotracks.at(itr).track2[2] + 1)
314 / (sqrt((vec_twotracks.at(itr).track1[0] * vec_twotracks.at(itr).track1[0]
315 + vec_twotracks.at(itr).track1[2] * vec_twotracks.at(itr).track1[2] + 1)
316 * (vec_twotracks.at(itr).track2[0] * vec_twotracks.at(itr).track2[0]
317 + vec_twotracks.at(itr).track2[2] * vec_twotracks.at(itr).track2[2] + 1)));
318 Double_t alpha = TMath::ACos(CosA);
319 Double_t alphaDeg = TMath::RadToDeg() * alpha;
320 if (fDebug)
321 cout << " alpha " << alpha << " CosA " << CosA << " ToDeg " << alphaDeg << endl;
322 hAngle_reco->Fill(alphaDeg);
323 if (ifNmc2)
324 hAngle_recoifNmc2->Fill(alphaDeg);
325 if (ifNmc2 && ifSelectedCases)
326 hAngle_recoifNmc2Cases->Fill(alphaDeg);
327 }
328
329 } // Ngood_mc_tracks == 2
330
331} // MCefficiencyCalculation
332
333//-----------------------------ReadSiliconTracks------------------------
334void BmnUpstreamTracking::ReadSiliconTracks(Double_t** par_ab, Double_t* par_z, Int_t& NTracks, vector<MC_points>& vec)
335{
336 MC_points tmpTr;
337 if (!expData) {
338 if (fDebug)
339 cout << " MC_SiliconTracks:" << endl;
340 int tr_before = -1;
341 int Nmc_tracks = -1;
342
343 int Nst_mc[kMaxMC];
344 int Np_in_3st[kMaxMC];
345 int mcTracksArray[kMaxMC];
346 double Xmc[kMaxMC][fNstations];
347 double Ymc[kMaxMC][fNstations];
348 double Zmc[kMaxMC][fNstations];
349 for (Int_t Id = 0; Id < kMaxMC; Id++) {
350 Nst_mc[Id] = 0;
351 Np_in_3st[Id] = 0;
352 mcTracksArray[Id] = -1;
353 for (Int_t i = 0; i < fNstations; i++) {
354 Xmc[Id][i] = -999.;
355 Ymc[Id][i] = -999.;
356 Zmc[Id][i] = -999.;
357 }
358 }
359
360 for (Int_t iMC = 0; iMC < fSiTracksSim->GetEntriesFast(); ++iMC) {
361 BmnSiliconHit* hit = (BmnSiliconHit*)fSiTracksSim->UncheckedAt(iMC);
362
363 Double_t x_MC = hit->GetX();
364 Double_t y_MC = hit->GetY();
365 Double_t z_MC = hit->GetZ();
366 Int_t trackId_MC = hit->GetIndex();
367 Int_t st_mc = -1;
368
369 if (tr_before != trackId_MC) {
370 Nmc_tracks++;
371 mcTracksArray[Nmc_tracks] = trackId_MC;
372 }
373 tr_before = trackId_MC;
374
375 if (z_MC > -320.) {
376 st_mc = 3;
377 Np_in_3st[Nmc_tracks]++;
378 }
379 if (z_MC < -434. && z_MC > -436.)
380 st_mc = 2;
381 if (z_MC < -438.)
382 st_mc = 1;
383 if (fDebug)
384 cout << " Id " << trackId_MC << " x " << x_MC << " y " << y_MC << " z " << z_MC << " st_mc " << st_mc
385 << endl;
386
387 Xmc[Nmc_tracks][st_mc] = x_MC;
388 Ymc[Nmc_tracks][st_mc] = y_MC;
389 Zmc[Nmc_tracks][st_mc] = z_MC;
390 Nst_mc[Nmc_tracks]++;
391
392 } // iMC
393 Nmc_tracks++;
394
395 for (Int_t Id = 0; Id < kMaxMC; Id++) {
396
397 for (Int_t i = 0; i < fNstations; i++) {
398 tmpTr.x_si[i] = Xmc[Id][i];
399 tmpTr.y_si[i] = Ymc[Id][i];
400 tmpTr.z_si[i] = Zmc[Id][i];
401 }
402
403 tmpTr.Id = mcTracksArray[Id];
404 tmpTr.Np_Si = Nst_mc[Id];
405 if (Np_in_3st[Id] == tmpTr.Np_Si)
406 tmpTr.np_3si = 1;
407 if (fDebug && Np_in_3st[Id] > 0)
408 cout << " Np_in_3st " << Np_in_3st[Id] << " tmpTr.Np_Si " << tmpTr.Np_Si << " tmpTr.np_3si "
409 << tmpTr.np_3si << endl;
410 if (Nst_mc[Id] >= 2 && Np_in_3st[Id] > 0)
411 vec.push_back(tmpTr);
412 }
413
414 if (fDebug)
415 cout << " vec_points.size() " << vec.size() << endl;
416
417 } // !expData
418 if (fDebug)
419 cout << endl;
420
421 if (fDebug)
422 cout << "----Reco_SiliconTracks:" << endl;
423 for (Int_t iTrack = 0; iTrack < fSiTracks->GetEntriesFast(); iTrack++) {
424 BmnTrack* track = (BmnTrack*)fSiTracks->At(iTrack);
425
426 FairTrackParam* par = track->GetParamLast();
427 // FairTrackParam* par = track->GetParamFirst();
428
429 if (fDebug)
430 cout << " SiTr " << iTrack << " SiTx " << par->GetTx() << " SiX " << par->GetX() << " SiTy " << par->GetTy()
431 << " SiY " << par->GetY() << " SiZ " << par->GetZ() << " Nhits " << track->GetNHits() << endl;
432
433 par_ab[0][NTracks] = par->GetTx();
434 par_ab[1][NTracks] = par->GetX();
435 par_ab[2][NTracks] = par->GetTy();
436 par_ab[3][NTracks] = par->GetY();
437 par_ab[4][NTracks] = track->GetNHits();
438 par_z[NTracks] = par->GetZ();
439 NTracks++;
440 }
441
442 if (fDebug)
443 cout << endl;
444
445} // ReadSiliconTracks
446//----------------------------------------------------------------------
447
448//----------------------------------------------------------------------
449void BmnUpstreamTracking::ReadSiliconHits(Double_t*** hits, Int_t* NSihits)
450{
451 //-------------Read Silicon Hits------------------------------------
452 if (fDebug)
453 cout << " Hits from Tracks " << endl;
454 for (Int_t iSihit = 0; iSihit < fSilHits->GetEntriesFast(); iSihit++) {
455 BmnSiliconHit* hit = (BmnSiliconHit*)fSilHits->At(iSihit);
456
457 // if ( fDebug ) cout<<" iSihit "<<iSihit<<endl;
458 Int_t st = hit->GetStation();
459
460 if (fDebug)
461 cout << " st " << st << " hit: X " << hit->GetX() << " Y " << hit->GetY() << " Sigx " << hit->GetDx()
462 << " Sigy " << hit->GetDy() << endl;
463 hits[st][NSihits[st]][0] = hit->GetX();
464 hits[st][NSihits[st]][1] = hit->GetY();
465 hits[st][NSihits[st]][2] = hit->GetZ();
466 hits[st][NSihits[st]][3] = hit->GetDx();
467 hits[st][NSihits[st]][4] = hit->GetDy();
468 hits[st][NSihits[st]][5] = hit->GetModule();
469
470 NSihits[st]++;
471 }
472} // ReadSiliconHits
473//----------------------------------------------------------------------
474
475//----------------------------------------------------------------------
476void BmnUpstreamTracking::ReadMWPCSegments(Double_t*** par_ab, Double_t** par_z, Int_t* Nseg, vector<MC_points>& vec)
477{
478 for (Int_t iSeg = 0; iSeg < fMWPCSegments->GetEntries(); ++iSeg) {
479 BmnTrack* segment = (BmnTrack*)fMWPCSegments->At(iSeg);
480
481 Int_t ich = -1;
482 Double_t Z_mwpc = segment->GetParamFirst()->GetZ();
483 // Double_t Chi2 = segment->GetChi2();
484
485 if (Z_mwpc >= -859. && Z_mwpc <= -857.)
486 ich = 0;
487 if (Z_mwpc >= -760. && Z_mwpc <= -758.)
488 ich = 1;
489 if (Z_mwpc >= -365. && Z_mwpc <= -363.)
490 ich = 2;
491 if (Z_mwpc >= -215. && Z_mwpc <= -213.)
492 ich = 3;
493
494 if (ich < 0)
495 continue;
496
497 if (fDebug)
498 cout << " MWPCSeg: ich " << ich << " Tx " << segment->GetParamFirst()->GetTx() << " X "
499 << segment->GetParamFirst()->GetX() << " Ty " << segment->GetParamFirst()->GetTy() << " Y "
500 << segment->GetParamFirst()->GetY() << " Z " << Z_mwpc << " Nhits " << segment->GetNHits() << endl;
501 par_ab[ich][0][Nseg[ich]] = segment->GetParamFirst()->GetTx();
502 par_ab[ich][1][Nseg[ich]] = segment->GetParamFirst()->GetX();
503 par_ab[ich][2][Nseg[ich]] = segment->GetParamFirst()->GetTy();
504 par_ab[ich][3][Nseg[ich]] = segment->GetParamFirst()->GetY();
505 par_ab[ich][4][Nseg[ich]] = segment->GetNHits();
506 par_ab[ich][5][Nseg[ich]] = segment->GetParamFirst()->GetCovariance(0, 0);
507 par_ab[ich][6][Nseg[ich]] = segment->GetParamFirst()->GetCovariance(1, 1);
508 par_ab[ich][7][Nseg[ich]] = segment->GetParamFirst()->GetCovariance(2, 2);
509 par_ab[ich][8][Nseg[ich]] = segment->GetParamFirst()->GetCovariance(3, 3);
510 par_z[ich][Nseg[ich]] = segment->GetParamFirst()->GetZ();
511
512 if (ich > 1) {
513 par_ab[ich][1][Nseg[ich]] += -X_shift - X_shift_seg[ich];
514 par_ab[ich][3][Nseg[ich]] += -Y_shift;
515 }
516 Nseg[ich]++;
517 }
518
519 if (!expData) {
520 if (fDebug)
521 cout << "----MC: MWPCSegments" << endl;
522
523 Int_t mcTracksArray[kMaxMC];
524 Int_t mcTracksPdg[kMaxMC];
525 Double_t X2mc[kMaxMC][kNPlanes];
526 // Double_t Y2mc[kMaxMC][kNPlanes];
527 // Double_t Z2mc[kMaxMC][kNPlanes];
528 Double_t X3mc[kMaxMC][kNPlanes];
529 // Double_t Y3mc[kMaxMC][kNPlanes];
530 // Double_t Z3mc[kMaxMC][kNPlanes];
531 int Npl_MC2[kMaxMC];
532 int Npl_MC3[kMaxMC];
533
534 for (Int_t Id = 0; Id < kMaxMC; Id++) {
535 Npl_MC2[Id] = 0;
536 Npl_MC3[Id] = 0;
537 mcTracksArray[Id] = -1;
538 mcTracksPdg[Id] = -1;
539 for (Int_t i = 0; i < kNPlanes; i++) {
540 X2mc[Id][i] = -999.;
541 // Y2mc[Id][i] = -999.;
542 // Z2mc[Id][i] = -999.;
543 X3mc[Id][i] = -999.;
544 // Y3mc[Id][i] = -999.;
545 // Z3mc[Id][i] = -999.;
546 }
547 }
548
549 int tr_before = -1;
550 int Nmc_tracks = -1;
551
552 for (Int_t iMC = 0; iMC < fBmnHitsArray->GetEntriesFast(); ++iMC) {
553 BmnMwpcHit* hit = (BmnMwpcHit*)fBmnHitsArray->UncheckedAt(iMC);
554
555 Int_t st_MC = hit->GetMwpcId();
556 Int_t trackId_MC = hit->GetHitId();
557 // Int_t pdg_MC = hit->GetPdgId();
558 Int_t pl_MC = hit->GetPlaneId();
559 // Short_t wire_MC = hit->GetWireNumber();
560 // Double_t time_MC = hit->GetWireTime();
561
562 // if (fDebug)cout<<" st_MC "<<st_MC<<" trackId_MC "<<trackId_MC<<" pl_MC "<<pl_MC<<" X "<<hit->GetX()<<"
563 // wire_MC "<<wire_MC<<endl;
564
565 if (tr_before != trackId_MC) {
566 Nmc_tracks++;
567 mcTracksArray[Nmc_tracks] = hit->GetHitId();
568 mcTracksPdg[Nmc_tracks] = hit->GetPdgId();
569 }
570 tr_before = trackId_MC;
571
572 if (st_MC == 2) {
573 X2mc[Nmc_tracks][pl_MC] = hit->GetX();
574 // Y2mc[Nmc_tracks][pl_MC] = hit->GetY();
575 // Z2mc[Nmc_tracks][pl_MC] = hit->GetZ();
576 Npl_MC2[Nmc_tracks]++;
577 }
578 if (st_MC == 3) {
579 X3mc[Nmc_tracks][pl_MC] = hit->GetX();
580 // Y3mc[Nmc_tracks][pl_MC] = hit->GetY();
581 // Z3mc[Nmc_tracks][pl_MC] = hit->GetZ();
582 Npl_MC3[Nmc_tracks]++;
583 }
584 // if (fDebug)cout<<" X2["<<Nmc_tracks<<"]["<<pl_MC<<"] "<<X2mc[Nmc_tracks][pl_MC]<<" Npl_MC2 "<<
585 // Npl_MC2[Nmc_tracks]<<" X3["<<Nmc_tracks<<"]["<<pl_MC<<"] "<<X3mc[Nmc_tracks][pl_MC]<<" Npl_MC3 "<<
586 // Npl_MC3[Nmc_tracks]<<" z3 "<<Z3mc[Nmc_tracks][pl_MC]<<endl;
587
588 } // iMC
589 Nmc_tracks++;
590
591 if (fDebug)
592 cout << " Nmc_tracks " << Nmc_tracks << " MC vec_points.size() " << vec.size() << endl;
593 MC_points tmpTr;
594 // bool mc_new = 0;
595 for (Int_t id = 0; id < Nmc_tracks; id++) { // mc_tr_curr
596 if (fDebug)
597 cout << " id " << id << " Id_mc " << mcTracksArray[id] << " pdg " << mcTracksPdg[id] << " Npl2 "
598 << Npl_MC2[id] << " Npl3 " << Npl_MC3[id] << endl;
599
600 for (size_t itr = 0; itr < vec.size(); itr++) { // mc_tr_before
601
602 if (Npl_MC2[id] >= 4 || Npl_MC3[id] >= 4) {
603 if (vec.at(itr).Id == mcTracksArray[id]) {
604 vec.at(itr).Pdg = mcTracksPdg[id];
605 vec.at(itr).Np_ch2 = Npl_MC2[id];
606 vec.at(itr).Np_ch3 = Npl_MC3[id];
607 vec.at(itr).Np_P = Npl_MC2[id] + Npl_MC3[id];
608 vec.at(itr).Np = vec.at(itr).Np_P + vec.at(itr).Np_Si;
609
610 int i3 = -1;
611 if (X3mc[id][5] > -900.)
612 i3 = 5;
613 if (i3 < 0 && X3mc[id][4] > -900.)
614 i3 = 4;
615 if (i3 < 0 && X3mc[id][3] > -900.)
616 i3 = 3;
617
618 int isi = -1;
619 if (isi < 0 && vec.at(itr).x_si[1] > -900.)
620 isi = 1;
621 if (isi < 0 && vec.at(itr).x_si[2] > -900.)
622 isi = 2;
623
624 vec.at(itr).param[0] = (vec.at(itr).x_si[isi] - vec.at(itr).x_si[3])
625 / (vec.at(itr).z_si[isi] - vec.at(itr).z_si[3]);
626 vec.at(itr).param[1] =
627 vec.at(itr).param[0] * (-350 - vec.at(itr).z_si[3]) + vec.at(itr).x_si[3];
628 vec.at(itr).param[2] = (vec.at(itr).y_si[isi] - vec.at(itr).y_si[3])
629 / (vec.at(itr).z_si[isi] - vec.at(itr).z_si[3]);
630 vec.at(itr).param[3] =
631 vec.at(itr).param[2] * (-350 - vec.at(itr).z_si[3]) + vec.at(itr).y_si[3];
632
633 if (fDebug)
634 cout << " si mc: x3 " << vec.at(itr).x_si[3] << " z3 " << vec.at(itr).z_si[3] << " UpX "
635 << vec.at(itr).param[1] << endl;
636
637 if (X2mc[id][0] > -900. || X2mc[id][3] > -900.)
638 vec.at(itr).xWas2 = 1;
639 if (X2mc[id][1] > -900. || X2mc[id][4] > -900.)
640 vec.at(itr).vWas2 = 1;
641 if (X2mc[id][2] > -900. || X2mc[id][5] > -900.)
642 vec.at(itr).uWas2 = 1;
643
644 if (X3mc[id][0] > -900. || X3mc[id][3] > -900.)
645 vec.at(itr).xWas3 = 1;
646 if (X3mc[id][1] > -900. || X3mc[id][4] > -900.)
647 vec.at(itr).vWas3 = 1;
648 if (X3mc[id][2] > -900. || X3mc[id][5] > -900.)
649 vec.at(itr).uWas3 = 1;
650
651 } // vec_points.at(itr).Id == mcTracksArray[id]
652 } // Npl_MC2[id] >= 4 || Npl_MC3[id] >= 4
653
654 } // vec_points.size
655 } // Nmc_tracksId = -1;
656
657 if (fDebug) {
658 cout << endl;
659 cout << "----MCtrue all systems, size: " << vec.size() << endl;
660 for (size_t itr = 0; itr < vec.size(); itr++) {
661 cout << " itr " << itr << " Id " << vec.at(itr).Id << " Np_Si " << vec.at(itr).Np_Si << " Np_P "
662 << vec.at(itr).Np_P << " Np_ch2 " << vec.at(itr).Np_ch2 << " Np_ch3 " << vec.at(itr).Np_ch3
663 << endl;
664 cout << " itr " << itr << " Id " << vec.at(itr).Id << " Up: ax " << vec.at(itr).param[0] << " bx "
665 << vec.at(itr).param[1] << " ay " << vec.at(itr).param[2] << " by " << vec.at(itr).param[3]
666 << endl;
667
668 cout << " Xt " << vec.at(itr).param[0] * (kZ_target - Zcentr) + vec.at(itr).param[1] << " Yt "
669 << vec.at(itr).param[2] * (kZ_target - Zcentr) + vec.at(itr).param[3] << endl;
670
671 hAx_upmc->Fill(vec.at(itr).param[0]);
672 hX_upmc->Fill(vec.at(itr).param[1]);
673 hAy_upmc->Fill(vec.at(itr).param[2]);
674 hY_upmc->Fill(vec.at(itr).param[3]);
675 }
676 } // if
677
678 } // !expData
679 if (fDebug)
680 cout << endl;
681} // ReadMWPCSegments
682//----------------------------------------------------------------------
683
684void BmnUpstreamTracking::ReadMWPCTracks(Double_t*** par_ab, Double_t** par_z, Int_t* NTracks)
685{
686 for (Int_t iTr = 0; iTr < fMWPCTracks->GetEntriesFast(); ++iTr) {
687 BmnTrack* track = (BmnTrack*)fMWPCTracks->At(iTr);
688
689 Int_t ip = -1;
690 Double_t Z = track->GetParamFirst()->GetZ();
691
692 if (Z >= -809. && Z <= -807.)
693 ip = 0;
694 if (Z >= -290. && Z <= -288.)
695 ip = 1;
696 if (ip < 0)
697 continue;
698
699 if (fDebug)
700 cout << " MWPCTr: ip " << iTr << " Tx " << track->GetParamFirst()->GetTx() << " X "
701 << track->GetParamFirst()->GetX() << " Ty " << track->GetParamFirst()->GetTy() << " Y "
702 << track->GetParamFirst()->GetY() << " Z " << Z << endl;
703
704 par_ab[ip][0][NTracks[ip]] = track->GetParamFirst()->GetTx();
705 par_ab[ip][1][NTracks[ip]] = track->GetParamFirst()->GetX();
706 par_ab[ip][2][NTracks[ip]] = track->GetParamFirst()->GetTy();
707 par_ab[ip][3][NTracks[ip]] = track->GetParamFirst()->GetY();
708 par_ab[ip][4][NTracks[ip]] = track->GetNHits();
709 par_ab[ip][5][NTracks[ip]] = track->GetParamFirst()->GetCovariance(0, 0);
710 par_ab[ip][6][NTracks[ip]] = track->GetParamFirst()->GetCovariance(1, 1);
711 par_ab[ip][7][NTracks[ip]] = track->GetParamFirst()->GetCovariance(2, 2);
712 par_ab[ip][8][NTracks[ip]] = track->GetParamFirst()->GetCovariance(3, 3);
713 par_z[ip][NTracks[ip]] = track->GetParamFirst()->GetZ();
714
715 par_ab[1][1][NTracks[ip]] += -X_shift;
716 par_ab[1][3][NTracks[ip]] += -Y_shift;
717 /*
718 hAxPCtr->Fill(par_ab[1][0][NTracks[ip]]);
719 hxPCtr->Fill(par_ab[1][1][NTracks[ip]]);
720 hAyPCtr->Fill(par_ab[1][2][NTracks[ip]]);
721 hyPCtr->Fill(par_ab[1][3][NTracks[ip]]);
722 hbeam_profile_pair_[ip]->Fill(track->GetParamFirst()->GetX(),track->GetParamFirst()->GetY());
723 */
724 NTracks[ip]++;
725
726 } // itr
727
728} // ReadMWPCTracks
729//----------------------------------------------------------------------
730
731//----------------------------------------------------------------------
732void BmnUpstreamTracking::PCTracksSiTracksMatching(Double_t** par_Si,
733 Double_t* par_Siz,
734 Int_t& NTracksSi,
735 Double_t*** par_PC,
736 Double_t** par_PCz,
737 Int_t* NTracksPC,
738 vector<Smatch>& vec_m,
739 vector<Smatch>& Out)
740{
741
742 Smatch tmpSeg;
743 vec_m.clear();
744
745 if (fDebug)
746 cout << " NTracksSi " << NTracksSi << " NPCTracks[1] " << NPCTracks[1] << endl;
747
748 for (Int_t iSi = 0; iSi < NTracksSi; ++iSi) {
749 for (Int_t iPCtr = 0; iPCtr < NPCTracks[1]; ++iPCtr) {
750
751 Double_t dz = par_PCz[1][iPCtr] - par_Siz[iSi];
752
753 Double_t dAx = par_Si[0][iSi] - par_PC[1][0][iPCtr];
754 Double_t dAy = par_Si[2][iSi] - par_PC[1][2][iPCtr];
755 Double_t dX = (par_Si[0][iSi] * dz + par_Si[1][iSi]) - par_PC[1][1][iPCtr];
756 Double_t dY = (par_Si[2][iSi] * dz + par_Si[3][iSi]) - par_PC[1][3][iPCtr];
757 /*
758 hdX_MWPCpair1_Sitr->Fill(dX);
759 hdAx_MWPCpair1_Sitr->Fill(dAx);
760 hdY_MWPCpair1_Sitr->Fill(dY);
761 hdAy_MWPCpair1_Sitr->Fill(dAy);
762 */
763 if (fDebug)
764 cout << " Si(" << iSi << ")-P1(" << iPCtr << "): dAx(" << cutAx << ") " << dAx << " dAy(" << cutAy
765 << ") " << dAy << " dX(" << cutX << ") " << dX << " dY(" << cutY << ") " << dY << endl;
766 if (fabs(dAx) < cutAx && fabs(dAy) < cutAy && fabs(dX) < cutX && fabs(dY) < cutY) {
767 if (fDebug)
768 cout << "after cuts: Si " << iSi << "- P1 " << iPCtr << endl;
769 Double_t Chi2_match = dX * dX + 0.25 * dY * dY;
770
771 tmpSeg.Ind1 = iSi;
772 tmpSeg.Ind2 = iPCtr;
773 tmpSeg.Inds2 = -1;
774 tmpSeg.Inds3 = -1;
775 tmpSeg.Nhits1 = par_Si[4][iSi];
776 tmpSeg.Nhits2 = par_PC[1][4][iPCtr];
777 tmpSeg.Z1 = par_Siz[iSi];
778 tmpSeg.Z2 = par_PCz[1][iPCtr];
779 tmpSeg.distX = dX;
780 tmpSeg.distY = dY;
781 tmpSeg.Chi2m = Chi2_match;
782 // if (fDebug) cout<<" iSi "<<iSi<<" iPCtr "<<iPCtr<<" dX "<<dx<<" dY "<<dy<<" Chi2_match
783 // "<<Chi2_match<<endl;
784
785 for (int ipar = 0; ipar < 4; ipar++) {
786 tmpSeg.param1[ipar] = par_Si[ipar][iSi];
787 tmpSeg.param2[ipar] = par_PC[1][ipar][iPCtr];
788 }
789 vec_m.push_back(tmpSeg);
790
791 /*
792 hdX_MWPCpair1_Sitr_cut->Fill(dX);
793 hdAx_MWPCpair1_Sitr_cut->Fill(dAx);
794 hdY_MWPCpair1_Sitr_cut->Fill(dY);
795 hdAy_MWPCpair1_Sitr_cut->Fill(dAy);
796 */
797 } // if
798 } // iPCtr
799 } // iSi
800 if (fDebug)
801 cout << " --End PCtracks: size " << vtmpSeg.size() << endl;
802
803 if (vtmpSeg.size() > 0) {
804 Int_t PCtrack = 1;
805 RecordingBest(PCtrack, vec_m, Out);
806 }
807
808 if (fDebug)
809 cout << "PCTracksSiTracksMatching: OutVector.size() " << Out.size() << endl;
810 for (size_t iter = 0; iter < Out.size(); ++iter) {
811 if (fDebug)
812 printf(" OutVector.at(%zu): %8.4f | Si %d - PC %d Seg2 %d - Seg3 %d\n", iter, Out.at(iter).Chi2m,
813 Out.at(iter).Ind1, Out.at(iter).Ind2, Out.at(iter).Inds2, Out.at(iter).Inds3);
814 }
815} // PCTracksSiTracksMatching
816//----------------------------------------------------------------------
817
818//----------------------------------------------------------------------
819void BmnUpstreamTracking::PCSegmentsSiTracksMatching(Double_t** par_Si,
820 Double_t* par_Siz,
821 Int_t& NTracksSi,
822 Double_t*** par_PC,
823 Double_t** par_PCz,
824 Int_t* NSegPC,
825 vector<Smatch>& vec_m,
826 vector<Smatch>& Out)
827{
828
829 vec_m.clear();
830 if (fDebug)
831 cout << endl;
832 if (fDebug)
833 cout << " ----Matching: PC segments Ch2 ---" << endl;
834 if (fDebug)
835 cout << " NTracksSi " << NTracksSi << " NSegPC[2] " << NSegPC[2] << endl;
836 for (int iSi = 0; iSi < NTracksSi; ++iSi) {
837 for (int iseg2 = 0; iseg2 < NSegPC[2]; ++iseg2) {
838
839 Double_t dz2 = par_PCz[2][iseg2] - par_Siz[iSi];
840 // if ( fDebug ) cout<<" Zpc "<<par_PCz[2][iseg2]<<" Zsi "<<par_Siz[iSi]<<endl;
841 Double_t X_extr2 = par_Si[0][iSi] * dz2 + par_Si[1][iSi];
842 Double_t Y_extr2 = par_Si[2][iSi] * dz2 + par_Si[3][iSi];
843
844 Double_t dX = X_extr2 - par_PC[2][1][iseg2];
845 Double_t dY = Y_extr2 - par_PC[2][3][iseg2];
846
847 // hdX_MWPCSeg2_Sitr->Fill(dX);
848 // hdY_MWPCSeg2_Sitr->Fill(dY);
849
850 if (fDebug)
851 cout << " Si[" << iSi << "]-Seg2[" << iseg2 << "] : dX(" << cutX << ") " << dX << " dY(" << cutY << ") "
852 << dY << endl;
853 if (fabs(dX) < cutX && fabs(dY) < cutY) {
854
855 Int_t iseg3 = -1;
856 if (fDebug)
857 cout << " cuts ok:iseg2 " << iseg2 << endl;
858 RecordCandidateSeg(iSi, iseg2, iseg3, par_Si, par_PC, par_Siz, par_PCz, dX, dY, vec_m);
859
860 } // if
861 } // iseg2
862 } // iSi
863
864 if (fDebug)
865 cout << " --End PCiseg2 : now vec_m.size " << vec_m.size() << endl;
866 if (fDebug)
867 cout << endl;
868 //-----------------------------------------------------------------
869 if (fDebug)
870 cout << endl;
871 if (fDebug)
872 cout << " ----Matching: PC segments Ch3 ---" << endl;
873 if (fDebug)
874 cout << " NTracksSi " << NTracksSi << " NSegPC[3] " << NSegPC[3] << endl;
875 for (int iSi = 0; iSi < NSiTracks; ++iSi) {
876 for (int iseg3 = 0; iseg3 < NSegPC[3]; ++iseg3) {
877
878 Double_t dz3 = par_PCz[3][iseg3] - par_Siz[iSi];
879 Double_t X_extr3 = par_Si[0][iSi] * dz3 + par_Si[1][iSi];
880 Double_t Y_extr3 = par_Si[2][iSi] * dz3 + par_Si[3][iSi];
881
882 Double_t dX = X_extr3 - par_PC[3][1][iseg3];
883 Double_t dY = Y_extr3 - par_PC[3][3][iseg3];
884 // hdX_MWPCSeg3_Sitr->Fill(dX);
885 // hdY_MWPCSeg3_Sitr->Fill(dY);
886
887 if (fDebug)
888 cout << " Si[" << iSi << "]-Seg3[" << iseg3 << "] : dX(" << cutX << ") " << dX << " dY(" << cutY << ") "
889 << dY << endl;
890 if (fabs(dX) < cutX && fabs(dY) < cutY) {
891
892 Int_t iseg2 = -1;
893 RecordCandidateSeg(iSi, iseg2, iseg3, par_Si, par_PC, par_Siz, par_PCz, dX, dY, vec_m);
894
895 } // if
896 } // iseg3
897 } // iSi
898
899 if (fDebug)
900 cout << " --End PCiseg3 : now vec_m.size " << vec_m.size() << endl;
901
902 if (vec_m.size() > 0) {
903 Int_t PCtrack = -1;
904 RecordingBest(PCtrack, vec_m, Out);
905 }
906 if (fDebug)
907 cout << " vec_m.size() " << vec_m.size() << " Out.size() " << Out.size() << endl;
908 if (fDebug)
909 cout << endl;
910} // PCSegmentsSiTracksMatching
911//----------------------------------------------------------------------
912
913//----------------------------------------------------------------------
914void BmnUpstreamTracking::RecordCandidateSeg(Int_t& iSi,
915 Int_t& is2,
916 Int_t& is3,
917 Double_t** par_Si,
918 Double_t*** par_PC,
919 Double_t* par_Siz,
920 Double_t** ZPC,
921 Double_t& dx,
922 Double_t& dy,
923 vector<Smatch>& vtmp)
924{
925
926 if (fDebug)
927 cout << " RecordCandidateSeg: Inds: seg2 " << is2 << " seg3 " << is3 << endl;
928 Smatch tmpSeg;
929 Int_t iseg = is2;
930 Int_t ich = 2;
931 if (is2 < 0) {
932 iseg = is3;
933 ich = 3;
934 }
935
936 Double_t Chi2_match = dx * dx + 0.25 * dy * dy;
937
938 tmpSeg.Ind1 = iSi;
939 tmpSeg.Ind2 = -1;
940 tmpSeg.Inds2 = is2;
941 tmpSeg.Inds3 = is3;
942 tmpSeg.Z1 = par_Siz[iSi];
943 tmpSeg.Zs2 = ZPC[ich][iseg];
944 tmpSeg.distX = dx;
945 tmpSeg.distY = dy;
946 tmpSeg.Chi2m = Chi2_match;
947 if (fDebug)
948 cout << " iSi " << iSi << " iseg2 " << is2 << " iseg3 " << is3 << " dX " << dx << " dY " << dy << " Chi2_match "
949 << Chi2_match << endl;
950
951 for (int ipar = 0; ipar < 4; ipar++) {
952 tmpSeg.param1[ipar] = par_Si[ipar][iSi];
953 tmpSeg.params2[ipar] = par_PC[ich][ipar][iseg];
954 }
955 vtmp.push_back(tmpSeg);
956}
957//----------------------------------------------------------------------
958
959//----------------------------------------------------------------------
960void BmnUpstreamTracking::RecordingBest(Int_t& iPCtrack, vector<Smatch>& vec, vector<Smatch>& OutVec)
961{
962 if (iPCtrack > -1) {
963 if (fDebug)
964 cout << endl;
965 if (fDebug)
966 cout << "--RecordingBest: tracks --" << endl;
967 // vector sorting
968 sort(vec.begin(), vec.end(), compareSegments);
969
970 // first best
971 OutVec.push_back(vec.at(0));
972 if (fDebug)
973 cout << "--isMatch: Chi2 " << vec.at(0).Chi2m << endl;
974
975 // reject repeat index
976 Bool_t isMatch;
977 for (size_t iter = 0; iter < vec.size(); ++iter) {
978 if (fDebug)
979 printf(" vtmpSeg.at(%zu): %8.4f | %d - %d\n", iter, vec.at(iter).Chi2m, vec.at(iter).Ind1,
980 vec.at(iter).Ind2);
981 isMatch = 0;
982 for (size_t InIter = 0; InIter < OutVec.size(); ++InIter) {
983 if (vec.at(iter).Ind1 == OutVec.at(InIter).Ind1 || vec.at(iter).Ind2 == OutVec.at(InIter).Ind2) {
984 isMatch = 1;
985 continue;
986 }
987 }
988 // writing unique index
989 if (isMatch == 0) {
990 if (fDebug)
991 cout << "--isMatch: Chi2 " << vec.at(iter).Chi2m << endl;
992 OutVec.push_back(vec.at(iter));
993 }
994 } // iter
995
996 if (fDebug)
997 cout << " Now OutVec.size " << OutVec.size() << endl;
998 vec.clear();
999 } else { //---segs--
1000 if (fDebug)
1001 cout << "--RecordingBest: segs --" << endl;
1002 // vector sorting
1003 sort(vec.begin(), vec.end(), compareSegments);
1004
1005 // reject repeat index
1006 Bool_t isMatch;
1007 for (size_t iter = 0; iter < vec.size(); ++iter) {
1008 if (fDebug)
1009 printf(" vec.at(%zu): %8.4f | Si %d - Seg2 %d - Seg3 %d\n", iter, vec.at(iter).Chi2m, vec.at(iter).Ind1,
1010 vec.at(iter).Inds2, vec.at(iter).Inds3);
1011 isMatch = 0;
1012 for (size_t InIter = 0; InIter < OutVec.size(); ++InIter) {
1013 if (vec.at(iter).Ind1 == OutVec.at(InIter).Ind1 && OutVec.at(InIter).Ind2 > -1) {
1014 isMatch = 1;
1015 continue;
1016 }
1017 if ((vec.at(iter).Ind1 == OutVec.at(InIter).Ind1 && OutVec.at(InIter).Inds3 > -1)
1018 || (vec.at(iter).Ind1 == OutVec.at(InIter).Ind1 && vec.at(iter).Inds2 == OutVec.at(InIter).Inds2))
1019 {
1020 isMatch = 1;
1021 continue;
1022 }
1023 if ((vec.at(iter).Ind1 == OutVec.at(InIter).Ind1 && OutVec.at(InIter).Inds2 > -1)
1024 || (vec.at(iter).Ind1 == OutVec.at(InIter).Ind1 && vec.at(iter).Inds3 == OutVec.at(InIter).Inds3))
1025 {
1026 isMatch = 1;
1027 continue;
1028 }
1029 if ((OutVec.at(InIter).Inds2 > -1 && vec.at(iter).Inds2 == OutVec.at(InIter).Inds2)
1030 || (OutVec.at(InIter).Inds3 > -1 && vec.at(iter).Inds3 == OutVec.at(InIter).Inds3))
1031 {
1032 isMatch = 1;
1033 continue;
1034 }
1035 }
1036 if (fDebug)
1037 cout << " writing unique index " << endl;
1038 if (isMatch == 0) {
1039 if (fDebug)
1040 cout << "--isMatch: Chi2 " << vec.at(iter).Chi2m << " Inds2 " << vec.at(iter).Inds2 << " Inds3 "
1041 << vec.at(iter).Inds3 << endl;
1042 OutVec.push_back(vec.at(iter));
1043 }
1044 } // iter
1045
1046 if (fDebug)
1047 cout << " RecordingBest: Now OutVec.size " << OutVec.size() << endl;
1048 for (size_t iter = 0; iter < OutVec.size(); ++iter) {
1049 if (fDebug)
1050 cout << " iter " << iter << " index: Si:" << OutVec.at(iter).Ind1 << " PCtr: " << OutVec.at(iter).Ind2
1051 << " PCseg2: " << OutVec.at(iter).Inds2 << " PCseg3: " << OutVec.at(iter).Inds3 << endl;
1052 }
1053
1054 vec.clear();
1055 }
1056
1057} // RecordingBest
1058//----------------------------------------------------------------------
1059
1060//----------------------------------------------------------------------
1061void BmnUpstreamTracking::GetAddSiHits(vector<Smatch>& Out, Double_t*** hits, Int_t* Nh)
1062{
1063 Bool_t was_ind;
1064 Double_t tmp_array[kNumStSi][6];
1065 if (fDebug)
1066 cout << "--GetAddSiHits--" << endl;
1067
1068 for (size_t iter = 0; iter < Out.size(); ++iter) {
1069 if (fDebug)
1070 printf(" Out.at(%zu): %8.4f | Si %d - PC %d Seg2 %d - Seg3 %d\n", iter, Out.at(iter).Chi2m,
1071 Out.at(iter).Ind1, Out.at(iter).Ind2, Out.at(iter).Inds2, Out.at(iter).Inds3);
1072
1073 was_ind = 0;
1074 for (int st = 1; st < kNumStSi; ++st) {
1075 for (int i = 0; i < 5; ++i) {
1076 tmp_array[st][i] = -999.;
1077 }
1078 for (int ihit = 0; ihit < Nh[st]; ++ihit) {
1079 if (Out.at(iter).Ind1 == ihit) {
1080 was_ind = 1;
1081 tmp_array[st][0] = hits[st][ihit][0];
1082 tmp_array[st][1] = hits[st][ihit][1];
1083 tmp_array[st][2] = hits[st][ihit][2];
1084 tmp_array[st][3] = hits[st][ihit][3];
1085 tmp_array[st][4] = hits[st][ihit][4];
1086 tmp_array[st][5] = hits[st][ihit][5];
1087 }
1088 }
1089 } // st
1090 if (was_ind) {
1091 for (int st = 1; st < kNumStSi; ++st) {
1092 Out.at(iter).SiCoordX[st] = tmp_array[st][0];
1093 Out.at(iter).SiCoordY[st] = tmp_array[st][1];
1094 Out.at(iter).SiCoordZ[st] = tmp_array[st][2];
1095 Out.at(iter).SiSigmaX[st] = tmp_array[st][3];
1096 Out.at(iter).SiSigmaY[st] = tmp_array[st][4];
1097 Out.at(iter).SiMod[st] = tmp_array[st][5];
1098 if (fDebug)
1099 cout << " st " << st << " X " << tmp_array[st][0] << " Y " << tmp_array[st][1] << " Z "
1100 << tmp_array[st][2] << " mod " << tmp_array[st][5] << endl;
1101 }
1102 }
1103 } // iter
1104
1105} // GetAddSiHits
1106//----------------------------------------------------------------------
1107
1108//----------------------------------------------------------------------
1109void BmnUpstreamTracking::GetHitsToPoints(vector<Smatch>& Out,
1110 Double_t*** ab_tr,
1111 Int_t* NPCTracks_,
1112 Double_t*** ab_Ch,
1113 Int_t* Nseg_ch_,
1114 Double_t*** Points_,
1115 Int_t& NumPoints_)
1116{
1117 if (fDebug)
1118 cout << " ---GetHitsToPoints---" << endl;
1119 for (size_t iter = 0; iter < Out.size(); iter++) {
1120 if (fDebug)
1121 cout << " silicon points " << endl;
1122 for (int ist = 1; ist < kNumStSi; ++ist) {
1123 if (ist < 3) {
1124 Points_[iter][ist - 1][0] = Out.at(iter).SiCoordX[ist]; // x
1125 Points_[iter][ist - 1][1] = Out.at(iter).SiCoordY[ist]; // y
1126 Points_[iter][ist - 1][2] = Out.at(iter).SiCoordZ[ist]; // z
1127 Points_[iter][ist - 1][3] = Out.at(iter).SiSigmaX[ist]; // sigx
1128 Points_[iter][ist - 1][4] = Out.at(iter).SiSigmaY[ist]; // sigy
1129 if (fDebug)
1130 cout << " X[" << ist << "] " << Points_[iter][ist - 1][0] << " Y[" << ist << "] "
1131 << Points_[ist - 1][ist][1] << endl;
1132 } else {
1133 Points_[iter][ist][0] = Out.at(iter).SiCoordX[ist]; // x
1134 Points_[iter][ist][1] = Out.at(iter).SiCoordY[ist]; // y
1135 Points_[iter][ist][2] = Out.at(iter).SiCoordZ[ist]; // z
1136 Points_[iter][ist][3] = Out.at(iter).SiSigmaX[ist]; // sigx
1137 Points_[iter][ist][4] = Out.at(iter).SiSigmaY[ist]; // sigy
1138 if (fDebug)
1139 cout << " X[" << ist << "] " << Points_[iter][ist][0] << " Y[" << ist << "] "
1140 << Points_[iter][ist][1] << endl;
1141 }
1142 }
1143 for (int ipctr = 0; ipctr < NPCTracks_[1]; ++ipctr) {
1144 // if was MWPC-track
1145 if (Out.at(iter).Ind2 == ipctr) {
1146 // x = kx*(dz) + bx //Z_pair[kNumPairs] = {-808.816, -289.401}; //Z_Chamber[kNumChambers] = {-858.,
1147 // -761., -364., -214.};
1148 Double_t deltaAx = ab_tr[1][7][ipctr];
1149 Double_t deltaAy = ab_tr[1][8][ipctr];
1150
1151 Points_[iter][2][0] =
1152 Out.at(iter).param2[0] * (Z_Chamber[2] - Z_pair[1]) + Out.at(iter).param2[1]; // x3
1153 Points_[iter][4][0] = Out.at(iter).param2[0] * (Z_Chamber[3] - Z_pair[1])
1154 + Out.at(iter).param2[1]; // x4 //(Z_Chamber[3]-Z_pair[1])
1155
1156 Points_[iter][2][1] =
1157 Out.at(iter).param2[2] * (Z_Chamber[2] - Z_pair[1]) + Out.at(iter).param2[3]; // y3
1158 Points_[iter][4][1] = Out.at(iter).param2[2] * (Z_Chamber[3] - Z_pair[1])
1159 + Out.at(iter).param2[3]; // y4 //Z_pair[1] - Z_Chamber[3]
1160 // Ax*deltaX + X*deltaAx + deltaX
1161 Points_[iter][2][3] =
1162 Out.at(iter).param2[0] * SigmXmwpc + Out.at(iter).param2[1] * deltaAx + SigmXmwpc; // sigx
1163 Points_[iter][2][4] =
1164 Out.at(iter).param2[2] * SigmYmwpc + Out.at(iter).param2[3] * deltaAy + SigmYmwpc; // sigy
1165 Points_[iter][4][3] = SigmXmwpc; // sigx// Points_[iter][3][3];
1166 Points_[iter][4][4] = SigmYmwpc; // sigy// Points_[iter][3][4];
1167 if (fDebug) {
1168 cout << " MWPCs from track" << endl;
1169 cout << " X3 " << Points_[iter][2][0] << " X4 " << Points_[iter][4][0]
1170 << endl; //" sigX "<<Points_[iter][3][3]<<endl;
1171 cout << " Y3 " << Points_[iter][2][1] << " X4 " << Points_[iter][4][1]
1172 << endl; //" sigX "<<Points_[iter][3][3]<<endl;
1173 //" sigXcov "<<ab_tr[1][5][ipctr]<<endl;
1174 }
1175 }
1176 } // NPCTracks_[1]
1177 //--if was segment--
1178 for (int iseg2 = 0; iseg2 < Nseg_ch_[2]; ++iseg2) {
1179 if (Out.at(iter).Inds2 == iseg2) {
1180 Points_[iter][2][0] = ab_Ch[2][1][iseg2]; // x
1181 Points_[iter][2][1] = ab_Ch[2][3][iseg2]; // y
1182 Points_[iter][2][3] = SigmXmwpc; // sigx
1183 Points_[iter][2][4] = SigmYmwpc; // sigy
1184 if (fDebug) {
1185 cout << " MWPCs from segment" << endl;
1186 cout << " X3 " << Points_[iter][2][0] << " Y3 " << Points_[iter][2][1]
1187 << endl; //" sigX "<<Points_[iter][3][3]<<endl;
1188 //" sigXcov "<<ab_tr[1][5][ipctr]<<endl;
1189 }
1190 }
1191 } // Nseg_ch_[2]
1192 for (int iseg3 = 0; iseg3 < Nseg_ch_[3]; ++iseg3) {
1193 if (Out.at(iter).Inds3 == iseg3) {
1194 Points_[iter][4][0] = ab_Ch[3][1][iseg3]; // x
1195 Points_[iter][4][1] = ab_Ch[3][3][iseg3]; // y
1196 Points_[iter][4][3] = SigmXmwpc; // sigx
1197 Points_[iter][4][4] = SigmYmwpc; // sigy
1198 if (fDebug) {
1199 cout << " MWPCs from segment" << endl;
1200 cout << " X4 " << Points_[iter][4][0] << " Y4 " << Points_[iter][4][1] << endl;
1201 }
1202 }
1203 } // Nseg_ch_[2]
1204
1205 Points_[iter][2][2] = Z_Chamber[2];
1206 Points_[iter][4][2] = Z_Chamber[3];
1207 NumPoints_ = iter + 1;
1208 } // Out.size()
1209 //????????
1210}
1211//----------------------------------------------------------------------
1212
1213//----------------------------------------------------------------------
1214void BmnUpstreamTracking::CalculationOfChi2(Double_t*** Points_, Int_t& NumPoints_, vector<UpTracks>& vec)
1215{
1216 if (fDebug)
1217 cout << " Check points before fit " << " Zcentr " << Zcentr << endl;
1218 Int_t Nhits_before;
1219 for (int ip = 0; ip < NumPoints_; ++ip) {
1220 Nhits_before = 0;
1221 for (int ist = 0; ist < kPoints; ++ist) {
1222 Points_[ip][ist][0] += X_shiftUp[ist];
1223 Points_[ip][ist][1] += Y_shiftUp[ist];
1224 Points_[ip][ist][2] -= Zcentr;
1225 // if (fDebug) cout<<"ip "<<ip<<" X["<<ist<<"] "<< Points_[ip][ist][0]<<" Y "<<Points_[ip][ist][1]<<" Z
1226 // "<<Points_[ip][ist][2]<<" sigx "<<Points_[ip][ist][3]<<" sigy "<<Points_[ip][ist][4]<<endl;
1227 if (Points_[ip][ist][0] > -900.)
1228 Nhits_before++;
1229 if (Points_[ip][ist][1] > -900.)
1230 Nhits_before++;
1231 }
1232 }
1233 if (fDebug)
1234 cout << " Nhits_before " << Nhits_before << endl;
1235 UpTracks tmptr;
1236
1237 for (Int_t i = 0; i < kPoints; ++i) {
1238 for (Int_t j = 0; j < 5; ++j) {
1239 Xforglfit[i][j] = -999.;
1240 }
1241 }
1242 Double_t ChiSquare_;
1243 for (Int_t itr = 0; itr < NumPoints_; ++itr) {
1244
1245 if (Nhits_before < 6)
1246 continue;
1247
1248 ChiSquare_ = 0.;
1249 if (fDebug)
1250 cout << endl;
1251 if (fDebug)
1252 cout << "--GlobalFit--" << endl;
1253 // setting global fit matrices to zero
1254
1255 for (Int_t i = 0; i < 4; ++i) {
1256 for (Int_t j = 0; j < 4; ++j) {
1257 Amatr[i][j] = 0.;
1258 }
1259 }
1260 for (Int_t i = 0; i < 4; ++i) {
1261 rhs[i] = 0;
1262 }
1263 for (Int_t i = 0; i < 16; ++i) {
1264 AmatrArray[i] = 0;
1265 }
1266 for (Int_t i = 0; i < 16; ++i) {
1267 AmatrInverted[i] = 0;
1268 }
1269
1270 for (Int_t i = 0; i < kPoints; ++i) {
1271 for (Int_t j = 0; j < kPoin_Par; ++j) {
1272 Xforglfit[i][j] = Points_[itr][i][j];
1273 // if ( fDebug && j == 0) cout<<" x ["<<i<<"]["<<j<<"] "<<Xforglfit[i][j]<<endl;
1274 }
1275 for (Int_t j = 0; j < kPoin_Par; ++j) {
1276 // if ( fDebug && j == 1) cout<<" y ["<<i<<"]["<<j<<"] "<<Xforglfit[i][j]<<endl;
1277 }
1278 }
1279
1280 GlobalFit(Xforglfit, Amatr, rhs);
1281
1282 // here we do smth with Amatr & rhs (solving)
1283 // Amatr is row-major
1284 for (Int_t row = 0; row < 4; ++row) {
1285 for (Int_t col = 0; col < 4; ++col) {
1286 AmatrArray[4 * col + row] = Amatr[row][col];
1287 }
1288 } // AmatrArray is col-major
1289
1290 bool ifSuccess = InvertMatrix(AmatrArray, AmatrInverted);
1291
1292 if (!ifSuccess) {
1293 cout << "???? InvertMatrix is unsuccessfull ????" << endl;
1294 }
1295 if (!ifSuccess) {
1296 cout << "???? InvertMatrix is unsuccessfull ????" << endl;
1297 return;
1298 }
1299 // line_ = { ax, bx, ay, by}
1300 Double_t line_[4] = {0., 0., 0., 0.};
1301 for (Int_t coeff = 0; coeff < 4; ++coeff) {
1302 for (Int_t iCol = 0; iCol < 4; ++iCol) {
1303 // if ( fDebug ) cout<<" AmatrInverted["<<4*iCol + coeff<<"] "<<AmatrInverted[4*iCol + coeff]<<"
1304 // rhs["<<iCol<<"] "<<rhs[iCol]<<endl;
1305 line_[coeff] += AmatrInverted[4 * iCol + coeff] * rhs[iCol];
1306 }
1307 }
1308
1309 if (fDebug)
1310 cout << " ax " << line_[0] << " bx " << line_[1] << " ay " << line_[2] << " by " << line_[3] << endl;
1311 if (fDebug)
1312 cout << endl;
1313
1314 //---------------------------Chi2---------------------------------
1315 Int_t nhits = 0;
1316 for (Int_t ist = 0; ist < kPoints; ++ist) {
1317 if (Xforglfit[ist][0] < -900.)
1318 continue; // x
1319
1320 Double_t hit = Xforglfit[ist][0];
1321 Double_t fit = line_[0] * Xforglfit[ist][2] + line_[1];
1322 nhits++;
1323
1324 Double_t sigmaX_2 = Xforglfit[ist][3] * Xforglfit[ist][3];
1325 ChiSquare_ += ((hit - fit) * (hit - fit)) / (sigmaX_2);
1326 // if ( fDebug )
1327 // if ((ist == 2 || ist == 4) && fabs(hit-fit) > 0.2 ) cout<<" fEventNo "<<fEventNo<<" PC ist "<<ist<<" hitx
1328 // "<<hit<<" fit "<<fit<<" dx "<<hit-fit<<" Z "<<Xforglfit[ist][2]<<" sigma "<<Xforglfit[ist][3]<<"
1329 // ChiSquare_ "<<ChiSquare_<<endl; if ((ist == 0 || ist == 1 || ist == 3) && fabs(hit-fit) > 0.02 ) cout<<"
1330 // fEventNo "<<fEventNo<<" Si ist "<<ist<<" hitx "<<hit<<" fit "<<fit<<" dx "<<hit-fit<<" Z
1331 // "<<Xforglfit[ist][2]<<" sigma "<<Xforglfit[ist][3]<<" ChiSquare_ "<<ChiSquare_<<endl;
1332 if (fDebug)
1333 cout << " ist " << ist << " hitx " << hit << " fit " << fit << " dx " << hit - fit << " Z "
1334 << Xforglfit[ist][2] << " sigma " << Xforglfit[ist][3] << " ChiSquare_ " << ChiSquare_ << endl;
1335 }
1336 for (Int_t ist = 0; ist < kPoints; ++ist) {
1337 if (Xforglfit[ist][1] < -900.)
1338 continue; // y
1339
1340 Double_t hit = Xforglfit[ist][1];
1341 Double_t fit = line_[2] * Xforglfit[ist][2] + line_[3];
1342 nhits++;
1343
1344 Double_t sigmaY_2 = Xforglfit[ist][4] * Xforglfit[ist][4];
1345 ChiSquare_ += ((hit - fit) * (hit - fit)) / (sigmaY_2);
1346 // if ( fDebug)
1347 // if (fabs(hit-fit) > 0.2 )
1348 if (fDebug)
1349 cout << " ist " << ist << " hity " << hit << " fit " << fit << " dy " << hit - fit << " Z "
1350 << Xforglfit[ist][2] << " sigma " << Xforglfit[ist][4] << " ChiSquare_ " << ChiSquare_ << endl;
1351 } // ist
1352
1353 Double_t ChiSquareNdf = ChiSquare_ / (nhits - 4);
1354
1355 if (ChiSquareNdf < kChi2_Max && nhits > 5) {
1356
1357 tmptr.Chi2 = ChiSquareNdf;
1358
1359 if (fDebug)
1360 cout << " Chi2 " << tmptr.Chi2 << " nhits " << nhits << endl;
1361 if (fDebug)
1362 cout << " ------" << endl;
1363
1364 for (Int_t iCol = 0; iCol < 4; ++iCol) {
1365 tmptr.param[iCol] = line_[iCol];
1366 }
1367 // if ( fDebug) cout<<" Ax "<<tmptr.param[0]<<" bx "<<tmptr.param[1]<<" Ay "<<tmptr.param[2]<<" by
1368 // "<<tmptr.param[3]<<endl; if ( fDebug) cout<<" X coor "<<endl;
1369 for (Int_t ist = 0; ist < kPoints; ++ist) {
1370 tmptr.CoordZ[ist] = Xforglfit[ist][2];
1371 tmptr.CoordX[ist] = Xforglfit[ist][0];
1372 if (tmptr.CoordX[ist] < -900.)
1373 continue;
1374 // if ( fDebug) cout<<" ist "<<ist<<" hitx "<<tmptr.CoordX[ist]<<endl;
1375 }
1376 // if ( fDebug) cout<<" Y coor "<<endl;
1377 for (Int_t ist = 0; ist < kPoints; ++ist) {
1378 tmptr.CoordY[ist] = Xforglfit[ist][1];
1379 if (tmptr.CoordY[ist] < -900.)
1380 continue;
1381 // if ( fDebug) cout<<" ist "<<ist<<" hity "<<tmptr.CoordY[ist]<<endl;
1382 }
1383 tmptr.Nhits = nhits;
1384 vec.push_back(tmptr);
1385 }
1386 } // itr
1387
1388} // calculationOfChi2
1389//----------------------------------------------------------------------
1390
1391//----------------------------------------------------------------------
1392void BmnUpstreamTracking::PrintAllTracks(vector<UpTracks>& vecUp)
1393{
1394
1395 cout << " ------PrintAllTracks: " << vecUp.size() << endl;
1396 for (size_t itr = 0; itr < vecUp.size(); ++itr) {
1397 cout << endl;
1398 cout << " itr " << itr << " Chi2 " << vecUp.at(itr).Chi2 << " Ax " << vecUp.at(itr).param[0] << " bx "
1399 << vecUp.at(itr).param[1] << " Ay " << vecUp.at(itr).param[2] << " by " << vecUp.at(itr).param[3]
1400 << " Nhits " << vecUp.at(itr).Nhits << endl;
1401
1402 cout << "--X--" << endl;
1403 for (Int_t st = 0; st < kPoints; ++st) {
1404 if (vecUp.at(itr).CoordX[st] < -900.)
1405 continue;
1406 Double_t Hit = vecUp.at(itr).CoordX[st];
1407 Double_t Fit = vecUp.at(itr).param[0] * vecUp.at(itr).CoordZ[st] + vecUp.at(itr).param[1];
1408 hResXst[st]->Fill(Hit - Fit);
1409 cout << " st " << st << " Xhit " << Hit << " Xfit " << Fit << " Z " << vecUp.at(itr).CoordZ[st] << endl;
1410 }
1411 cout << "--Y--" << endl;
1412 for (Int_t st = 0; st < kPoints; ++st) {
1413 if (vecUp.at(itr).CoordY[st] < -900.)
1414 continue;
1415 Double_t Hit = vecUp.at(itr).CoordY[st];
1416 Double_t Fit = vecUp.at(itr).param[2] * vecUp.at(itr).CoordZ[st] + vecUp.at(itr).param[3];
1417 hResYst[st]->Fill(Hit - Fit);
1418 cout << " st " << st << " Yhit " << Hit << " Yfit " << Fit << " Z " << vecUp.at(itr).CoordZ[st] << endl;
1419 }
1420 Double_t x_target = vecUp.at(itr).param[0] * (kZ_target - Zcentr) + vecUp.at(itr).param[1];
1421 Double_t y_target = vecUp.at(itr).param[2] * (kZ_target - Zcentr) + vecUp.at(itr).param[3];
1422
1423 hAx_fitUp->Fill(vecUp.at(itr).param[0]);
1424 hAy_fitUp->Fill(vecUp.at(itr).param[2]);
1425 hx_fitUp->Fill(vecUp.at(itr).param[1]);
1426 hy_fitUp->Fill(vecUp.at(itr).param[3]);
1427 hvertexXYUp->Fill(x_target, y_target);
1428 hy_vs_x_Up->Fill(vecUp.at(itr).param[1], vecUp.at(itr).param[3]);
1429 hchi2_fitUp->Fill(vecUp.at(itr).Chi2);
1430 hNhitsUp->Fill(vecUp.at(itr).Nhits);
1431 hTyTx_Up->Fill(vecUp.at(itr).param[0], vecUp.at(itr).param[2]);
1432 }
1433 cout << "---" << endl;
1434}
1435//----------------------------------------------------------------------
1436
1437//----------------------------------------------------------------------
1438void BmnUpstreamTracking::TrackRecording(vector<UpTracks>& vecUp)
1439{
1440 for (size_t InIter = 0; InIter < vecUp.size(); ++InIter) {
1441 /*
1442 Double_t z0 = vecUp.at(InIter).CoordZ[1] + Zcentr;
1443 Double_t x0 = vecUp.at(InIter).param[0] * vecUp.at(InIter).CoordZ[0] + vecUp.at(InIter).param[1]; //ax(z)+bx
1444 Double_t y0 = vecUp.at(InIter).param[2] * vecUp.at(InIter).CoordZ[0] + vecUp.at(InIter).param[3]; //y
1445 Double_t z4 = vecUp.at(InIter).CoordZ[4] + Zcentr;
1446 Double_t x4 = vecUp.at(InIter).param[0] * vecUp.at(InIter).CoordZ[4] + vecUp.at(InIter).param[1]; //ax(z)+bx
1447 Double_t y4 = vecUp.at(InIter).param[2] * vecUp.at(InIter).CoordZ[4] + vecUp.at(InIter).param[3]; //y
1448
1449 FairTrackParam trackParamf;
1450 trackParamf.SetPosition(TVector3(x0, y0, z0));
1451 trackParamf.SetTx(vecUp.at(InIter).param[0]);
1452 trackParamf.SetTy(vecUp.at(InIter).param[2]);
1453
1454 FairTrackParam trackParaml;
1455 trackParaml.SetPosition(TVector3(x4, y4, z4));
1456 trackParaml.SetTx(vecUp.at(InIter).param[0]);
1457 trackParaml.SetTy(vecUp.at(InIter).param[2]);
1458 */
1459
1460 FairTrackParam trackParaml;
1461 trackParaml.SetPosition(TVector3(vecUp.at(InIter).param[1] + Shift_toCenterOfMagnetX,
1462 vecUp.at(InIter).param[3] + Shift_toCenterOfMagnetY, Zcentr));
1463 trackParaml.SetTx(vecUp.at(InIter).param[0] + Shift_toCenterOfMagnetAX);
1464 trackParaml.SetTy(vecUp.at(InIter).param[2] + Shift_toCenterOfMagnetAY);
1465
1466 BmnTrack* track = new ((*fBmnUpstreamTracksArray)[fBmnUpstreamTracksArray->GetEntriesFast()]) BmnTrack();
1467 track->SetChi2(vecUp.at(InIter).Chi2);
1468 track->SetNHits(vecUp.at(InIter).Nhits);
1469 track->SetFlag(InIter);
1470 // track->SetParamFirst(trackParamf);
1471 track->SetParamLast(trackParaml);
1472
1473 } // InIter
1474}
1475//----------------------------------------------------------------------
1476
1477//----------------------------------------------------------------------
1478bool BmnUpstreamTracking::InvertMatrix(Double_t* m, Double_t* invOut)
1479{
1480 Double_t inv[16], det;
1481 int i;
1482
1483 inv[0] = m[5] * m[10] * m[15] - m[5] * m[11] * m[14] - m[9] * m[6] * m[15] + m[9] * m[7] * m[14]
1484 + m[13] * m[6] * m[11] - m[13] * m[7] * m[10];
1485 inv[4] = -m[4] * m[10] * m[15] + m[4] * m[11] * m[14] + m[8] * m[6] * m[15] - m[8] * m[7] * m[14]
1486 - m[12] * m[6] * m[11] + m[12] * m[7] * m[10];
1487 inv[8] = m[4] * m[9] * m[15] - m[4] * m[11] * m[13] - m[8] * m[5] * m[15] + m[8] * m[7] * m[13]
1488 + m[12] * m[5] * m[11] - m[12] * m[7] * m[9];
1489 inv[12] = -m[4] * m[9] * m[14] + m[4] * m[10] * m[13] + m[8] * m[5] * m[14] - m[8] * m[6] * m[13]
1490 - m[12] * m[5] * m[10] + m[12] * m[6] * m[9];
1491 inv[1] = -m[1] * m[10] * m[15] + m[1] * m[11] * m[14] + m[9] * m[2] * m[15] - m[9] * m[3] * m[14]
1492 - m[13] * m[2] * m[11] + m[13] * m[3] * m[10];
1493 inv[5] = m[0] * m[10] * m[15] - m[0] * m[11] * m[14] - m[8] * m[2] * m[15] + m[8] * m[3] * m[14]
1494 + m[12] * m[2] * m[11] - m[12] * m[3] * m[10];
1495 inv[9] = -m[0] * m[9] * m[15] + m[0] * m[11] * m[13] + m[8] * m[1] * m[15] - m[8] * m[3] * m[13]
1496 - m[12] * m[1] * m[11] + m[12] * m[3] * m[9];
1497 inv[13] = m[0] * m[9] * m[14] - m[0] * m[10] * m[13] - m[8] * m[1] * m[14] + m[8] * m[2] * m[13]
1498 + m[12] * m[1] * m[10] - m[12] * m[2] * m[9];
1499 inv[2] = m[1] * m[6] * m[15] - m[1] * m[7] * m[14] - m[5] * m[2] * m[15] + m[5] * m[3] * m[14] + m[13] * m[2] * m[7]
1500 - m[13] * m[3] * m[6];
1501 inv[6] = -m[0] * m[6] * m[15] + m[0] * m[7] * m[14] + m[4] * m[2] * m[15] - m[4] * m[3] * m[14]
1502 - m[12] * m[2] * m[7] + m[12] * m[3] * m[6];
1503 inv[10] = m[0] * m[5] * m[15] - m[0] * m[7] * m[13] - m[4] * m[1] * m[15] + m[4] * m[3] * m[13]
1504 + m[12] * m[1] * m[7] - m[12] * m[3] * m[5];
1505 inv[14] = -m[0] * m[5] * m[14] + m[0] * m[6] * m[13] + m[4] * m[1] * m[14] - m[4] * m[2] * m[13]
1506 - m[12] * m[1] * m[6] + m[12] * m[2] * m[5];
1507 inv[3] = -m[1] * m[6] * m[11] + m[1] * m[7] * m[10] + m[5] * m[2] * m[11] - m[5] * m[3] * m[10] - m[9] * m[2] * m[7]
1508 + m[9] * m[3] * m[6];
1509 inv[7] = m[0] * m[6] * m[11] - m[0] * m[7] * m[10] - m[4] * m[2] * m[11] + m[4] * m[3] * m[10] + m[8] * m[2] * m[7]
1510 - m[8] * m[3] * m[6];
1511 inv[11] = -m[0] * m[5] * m[11] + m[0] * m[7] * m[9] + m[4] * m[1] * m[11] - m[4] * m[3] * m[9] - m[8] * m[1] * m[7]
1512 + m[8] * m[3] * m[5];
1513 inv[15] = m[0] * m[5] * m[10] - m[0] * m[6] * m[9] - m[4] * m[1] * m[10] + m[4] * m[2] * m[9] + m[8] * m[1] * m[6]
1514 - m[8] * m[2] * m[5];
1515 det = m[0] * inv[0] + m[1] * inv[4] + m[2] * inv[8] + m[3] * inv[12];
1516
1517 if (det == 0)
1518 return false;
1519 det = 1.0 / det;
1520 for (i = 0; i < 16; i++)
1521 invOut[i] = inv[i] * det;
1522 return true;
1523}
1524//----------------------------------------------------------------------
1525
1526//----------------------------------------------------------------------
1527void BmnUpstreamTracking::GlobalFit(Double_t** XYZsigXY, Double_t** Amatr_, Double_t* rhs_)
1528{
1529
1530 Double_t sigmaX_2 = 999.;
1531 Double_t sigmaY_2 = 999.;
1532 Double_t X = 999., Y = 999., Z = 999.;
1533
1534 for (Int_t iSt = 0; iSt < kPoints; ++iSt) {
1535
1536 Z = XYZsigXY[iSt][2];
1537
1538 if (XYZsigXY[iSt][0] > -900.) {
1539 X = XYZsigXY[iSt][0];
1540 sigmaX_2 = XYZsigXY[iSt][3] * XYZsigXY[iSt][3];
1541
1542 Amatr_[0][0] += 2 * Z * Z / sigmaX_2;
1543 Amatr_[0][1] += 2 * Z / sigmaX_2;
1544 rhs_[0] += 2 * Z * X / sigmaX_2;
1545
1546 Amatr_[1][0] += 2 * Z / sigmaX_2;
1547 Amatr_[1][1] += 2 / sigmaX_2;
1548 rhs_[1] += 2 * X / sigmaX_2;
1549 }
1550
1551 if (XYZsigXY[iSt][1] > -900.) {
1552 Y = XYZsigXY[iSt][1];
1553 sigmaY_2 = XYZsigXY[iSt][4] * XYZsigXY[iSt][4];
1554
1555 Amatr_[2][2] += 2 * Z * Z / sigmaY_2;
1556 Amatr_[2][3] += 2 * Z / sigmaY_2;
1557 rhs_[2] += 2 * Z * Y / sigmaY_2;
1558
1559 Amatr_[3][2] += 2 * Z / sigmaY_2;
1560 Amatr_[3][3] += 2 / sigmaY_2;
1561 rhs_[3] += 2 * Y / sigmaY_2;
1562 }
1563
1564 for (Int_t coeff = 0; coeff < 4; ++coeff) {
1565 for (Int_t iCol = 0; iCol < 4; ++iCol) {
1566 // if ( fDebug ) cout << "GLOBALFIT: station = " << iSt << " " << "Amatr_[" << coeff << "][" << iCol <<
1567 // "] = " << Amatr_[coeff][iCol] << "\n";
1568 }
1569 // if ( fDebug ) cout<<" rhs_["<<coeff<<"] "<< rhs_[coeff]<<endl;
1570 }
1571 }
1572
1573} // GlobalFit
1574//----------------------------------------------------------------------
1575
1576//----------------------------------------------------------------------
1577InitStatus BmnUpstreamTracking::Init()
1578{
1579
1580 if (fVerbose)
1581 cout << "BmnUpstreamTracking == Init ==" << endl;
1582 FairRootManager* ioman = FairRootManager::Instance();
1583
1584 //----Data & MC---
1585 fSiTracks = (TClonesArray*)ioman->GetObject(fInputBranchName1);
1586 if (!fSiTracks) {
1587 cout << "BmnUpstreamTracking::Init(): branch " << fInputBranchName1 << " not found! Task will be deactivated"
1588 << endl;
1589 SetActive(kFALSE);
1590 return kERROR;
1591 }
1592
1593 fSilHits = (TClonesArray*)ioman->GetObject(fInputBranchHits);
1594 if (!fSilHits) {
1595 cout << "BmnUpstreamTracking::Init(): branch " << fInputBranchHits << " not found! Task will be deactivated"
1596 << endl;
1597 SetActive(kFALSE);
1598 return kERROR;
1599 }
1600
1601 fMWPCTracks = (TClonesArray*)ioman->GetObject(fInputBranchName2);
1602 if (!fMWPCTracks) {
1603 cout << "BmnUpstreamTracking::Init(): branch " << fInputBranchName2 << " not found! Task will be deactivated"
1604 << endl;
1605 SetActive(kFALSE);
1606 return kERROR;
1607 }
1608 fMWPCSegments = (TClonesArray*)ioman->GetObject(fInputBranchName3);
1609 if (!fMWPCSegments) {
1610 cout << "BmnUpstreamTracking::Init(): branch " << fInputBranchName3 << " not found! Task will be deactivated"
1611 << endl;
1612 SetActive(kFALSE);
1613 return kERROR;
1614 }
1615
1616 //----MC true---
1617 if (!expData) {
1618 cout << " !expData " << endl;
1619 fSiTracksSim = (TClonesArray*)ioman->GetObject(fInputBranchNameSimTrue);
1620 if (!fSiTracksSim) {
1621 cout << "BmnUpstreamTracking::Init(): branch " << fInputBranchNameSimTrue
1622 << " not found! Task will be deactivated" << endl;
1623 SetActive(kFALSE);
1624 return kERROR;
1625 }
1626 fBmnHitsArray = (TClonesArray*)ioman->GetObject(fInputBranchName);
1627 cout << "fBmnHitsArray = " << fInputBranchName << "\n";
1628 if (!fBmnHitsArray) {
1629 cout << "BmnMwpcHitFinder::Init(): branch " << fInputBranchName << " not found! Task will be deactivated"
1630 << endl;
1631 SetActive(kFALSE);
1632 return kERROR;
1633 }
1634 } // !expData
1635
1636 // Create and register tracks arrays
1637 fBmnUpstreamTracksArray = new TClonesArray(fOutputTracksBranchName);
1638 ioman->Register("BmnUpstreamTrack", "UpstreamTrack", fBmnUpstreamTracksArray, kTRUE);
1639
1640 par_ab_Ch = new Double_t**[kNumChambers];
1641 par_ab_tr = new Double_t**[kNumPairs];
1642 SiXYhits = new Double_t**[kNumStSi];
1643 Points = new Double_t**[kBig];
1644 par_ab_SiTr = new Double_t*[kNumPars];
1645 par_ab_trz = new Double_t*[kNumPairs];
1646 par_Seg_z = new Double_t*[kNumChambers];
1647 Xforglfit = new Double_t*[kPoints];
1648 Amatr = new Double_t*[4];
1649 par_SiTr_z = new Double_t[kBig];
1650 X_shift_seg = new Double_t[kNumChambers];
1651 Z_pair = new Double_t[kNumPairs];
1652 Z_Chamber = new Double_t[kNumChambers];
1653 AmatrInverted = new Double_t[16];
1654 rhs = new Double_t[4];
1655 AmatrArray = new Double_t[16];
1656 line = new Double_t[4];
1657 Nseg_Ch = new Int_t[kNumChambers];
1658 NPCTracks = new Int_t[kNumPairs];
1659 NSiXYhits = new Int_t[kNumStSi];
1660 X_shiftUp = new Double_t[kPoints];
1661 Y_shiftUp = new Double_t[kPoints];
1662
1663 for (Int_t i = 0; i < 4; i++) {
1664 Amatr[i] = new Double_t[4];
1665 }
1666
1667 for (Int_t iPars = 0; iPars < kNumPars; iPars++) {
1668 par_ab_SiTr[iPars] = new Double_t[kBig];
1669 }
1670
1671 for (Int_t iCh = 0; iCh < kNumChambers; iCh++) {
1672 par_Seg_z[iCh] = new Double_t[kBig];
1673 par_ab_Ch[iCh] = new Double_t*[kNumPars];
1674 for (Int_t iPars = 0; iPars < kNumPars; iPars++) {
1675 par_ab_Ch[iCh][iPars] = new Double_t[kBig];
1676 }
1677 }
1678 for (Int_t ip = 0; ip < kNumPairs; ip++) {
1679 par_ab_tr[ip] = new Double_t*[kNumPars];
1680 par_ab_trz[ip] = new Double_t[kBig];
1681 for (Int_t iPars = 0; iPars < kNumPars; iPars++) {
1682 par_ab_tr[ip][iPars] = new Double_t[kBig];
1683 }
1684 }
1685 for (Int_t ist = 0; ist < kNumStSi; ist++) {
1686 SiXYhits[ist] = new Double_t*[kBig];
1687 for (Int_t ib = 0; ib < kBig; ib++) {
1688 SiXYhits[ist][ib] = new Double_t[kNumPars];
1689 }
1690 }
1691 for (Int_t ib = 0; ib < kBig; ib++) {
1692 Points[ib] = new Double_t*[kPoints];
1693 for (Int_t ip = 0; ip < kPoints; ip++) {
1694 Points[ib][ip] = new Double_t[kPoin_Par];
1695 }
1696 }
1697 for (Int_t ip = 0; ip < kPoints; ip++) {
1698 Xforglfit[ip] = new Double_t[kPoin_Par];
1699 }
1700
1701 //--some hists--
1702 if (fDebug) {
1703
1704 hAx_fitUp = new TH1D("Ax_fitUp", "Tx_fitUp;[rad]", 200, -0.03, 0.03);
1705 fList.Add(hAx_fitUp);
1706 hAy_fitUp = new TH1D("Ay_fitUp", "Ty_fitUp;[rad]", 200, -0.05, 0.05);
1707 fList.Add(hAy_fitUp);
1708 hx_fitUp = new TH1D("x_fitUp", "x_fitUp;[cm]", 200, -10, 10);
1709 fList.Add(hx_fitUp);
1710 hy_fitUp = new TH1D("y_fitUp", "y_fitUp;[cm]", 200, -15, 5);
1711 fList.Add(hy_fitUp);
1712 hy_vs_x_Up = new TH2D("y_vs_x_Up", ";X[cm];Y[cm]", 200, -10, 10, 200, -15, 5);
1713 fList.Add(hy_vs_x_Up);
1714 hvertexXYUp = new TH2D("vertexXYUp", "vertexUp;X[cm];Y[cm]", 200, -10, 10, 200, -10, 10);
1715 fList.Add(hvertexXYUp);
1716 hTyTx_Up = new TH2D("hTyTx_Up", "hTyTx_Up", 100, -0.04, 0.04, 100, -0.04, 0.04);
1717 fList.Add(hTyTx_Up);
1718
1719 hchi2_fitUp = new TH1D("chi2_fitUp", "Chi2", 100, 0, 100);
1720 fList.Add(hchi2_fitUp);
1721 hNhitsUp = new TH1D("NhitsUp", "Number of hits", 10, 0.5, 10.5);
1722 fList.Add(hNhitsUp);
1723
1724 hAx_upmc = new TH1D("Ax_upmc", "Ax_upmc;[rad]", 100, -0.05, 0.05);
1725 fList.Add(hAx_upmc);
1726 hAy_upmc = new TH1D("Ay_upmc", "Ay_upmc;[rad]", 100, -0.05, 0.05);
1727 fList.Add(hAy_upmc);
1728 hX_upmc = new TH1D("X_upmc", "X_upmc;[rad]", 100, -10, 10);
1729 fList.Add(hX_upmc);
1730 hY_upmc = new TH1D("Y_upmc", "Y_upmc;[rad]", 100, -10, 10);
1731 fList.Add(hY_upmc);
1732
1733 hdAx_tr_mc_comb = new TH1D("dAx_tr_mc_comb", "dAx_tr_mc_comb", 100, -0.005, 0.005);
1734 fList.Add(hdAx_tr_mc_comb);
1735 hdAy_tr_mc_comb = new TH1D("dAy_tr_mc_comb", "dAy_tr_mc_comb", 100, -0.005, 0.005);
1736 fList.Add(hdAy_tr_mc_comb);
1737 hdX_tr_mc_comb = new TH1D("dX_tr_mc_comb", "dX_tr_mc_comb", 100, -10, 10);
1738 fList.Add(hdX_tr_mc_comb);
1739 hdY_tr_mc_comb = new TH1D("dY_tr_mc_comb", "dY_tr_mc_comb", 100, -10, 10);
1740 fList.Add(hdY_tr_mc_comb);
1741
1742 hdAx_uptr_mc = new TH1D("dAx_uptr_mc", "dAx_uptr_mc", 200, -0.002, 0.002);
1743 fList.Add(hdAx_uptr_mc);
1744 hdAy_uptr_mc = new TH1D("dAy_uptr_mc", "dAy_uptr_mc", 200, -0.005, 0.005);
1745 fList.Add(hdAy_uptr_mc);
1746 hdX_uptr_mc = new TH1D("dX_uptr_mc", "dX_uptr_mc", 200, -2, 2);
1747 fList.Add(hdX_uptr_mc);
1748 hdY_uptr_mc = new TH1D("dY_uptr_mc", "dY_uptr_mc", 200, -2, 2);
1749 fList.Add(hdY_uptr_mc);
1750 hY_vs_Xmctrue = new TH2D("Y_vs_Xmctrue", "mctrue;X[cm];Y[cm]", 200, -10, 10, 200, -15, 5);
1751 fList.Add(hY_vs_Xmctrue);
1752
1753 hNtr_reco = new TH1D("Ntr_reco", "N up tr reco", 10, 0, 10);
1754 fList.Add(hNtr_reco);
1755 hNtr_mc = new TH1D("Ntr_mc", "N up tr mc", 10, 0, 10);
1756 fList.Add(hNtr_mc);
1757 hNtr_mc_vs_reco = new TH2D("Ntr_mc_vs_reco", "; Nupreco; Nupmctrue", 10, 0, 10, 10, 0, 10);
1758 fList.Add(hNtr_mc_vs_reco);
1759 hNrecoTrif2mc = new TH1D("NrecoTrif2mc", "N Up tracks if Nmc=2", 3, 0, 3);
1760 fList.Add(hNrecoTrif2mc);
1761 hAngle_reco = new TH1D("Angle_reco", "Angle_reco if Nreco = 2 ;[deg];N/0.02^{o}", 200, 0., 4.);
1762 fList.Add(hAngle_reco);
1763 hAngle_recoifNmc2 =
1764 new TH1D("Angle_recoifNmc2", "Angle_reco if Nreco = 2 & Nmc = 2;[deg];N/0.02^{o}", 200, 0., 4.);
1765 fList.Add(hAngle_recoifNmc2);
1766 hAngle_recoifNmc2Cases =
1767 new TH1D("Angle_recoifNmc2Cases", "Angle_reco if Nreco = 2 & Nmc = 2 & SelectedCases;[deg];N/0.02^{o}",
1768 200, 0., 4.);
1769 fList.Add(hAngle_recoifNmc2Cases);
1770
1771 hNum_mcuptr = new TH1D("hNum_mcuptr", "Num_mcuptr", 2, 0, 2);
1772 fList.Add(hNum_mcuptr);
1773 hDen_mcuptr = new TH1D("hDen_mcuptr", "Den_mcuptr", 2, 0, 2);
1774 fList.Add(hDen_mcuptr);
1775 hEff_mcuptr = new TH1D("hEff_mcuptr", "Eff_mc_uptr", 2, 0, 2);
1776 fList.Add(hEff_mcuptr);
1777 hEff_mcuptr->GetXaxis()->SetBinLabel(1, "2+3 coord");
1778 hEff_mcuptr->GetXaxis()->SetBinLabel(2, "only 3coord");
1779
1780 TH1D* h;
1781 for (Int_t i = 0; i < kPoints; ++i) {
1782 if (i == 2 || i == 4)
1783 h = new TH1D(Form("ResXst%d", i), Form("ResXst%d", i), 100, -0.5, 0.5);
1784 else
1785 h = new TH1D(Form("ResXst%d", i), Form("ResXst%d", i), 150, -0.015, 0.015);
1786 fList.Add(h);
1787 hResXst.push_back(h);
1788 }
1789 for (Int_t i = 0; i < kPoints; ++i) {
1790 if (i == 2 || i == 4)
1791 h = new TH1D(Form("ResYst%d", i), Form("ResYst%d", i), 100, -0.25, 0.25);
1792 // else if(i == 4) h = new TH1D(Form("ResYst%d", i), Form("ResYst%d", i), 100, -0.2, 0.2);
1793 else
1794 h = new TH1D(Form("ResYst%d", i), Form("ResYst%d", i), 100, -0.4, 0.4);
1795 fList.Add(h);
1796 hResYst.push_back(h);
1797 }
1798 }
1799
1800 Z_pair[0] = -808.816;
1801 Z_pair[1] = -289.401;
1802 Z_Chamber[0] = -858.155579;
1803 Z_Chamber[1] = -761.476562;
1804 Z_Chamber[2] = -364.255554;
1805 Z_Chamber[3] = -214.545547;
1806
1807 //--Shift--
1808 if (expData) {
1809 X_shift_seg[0] = 0.;
1810 X_shift_seg[1] = 0.;
1811 X_shift_seg[2] = -0.04;
1812 X_shift_seg[3] = -0.05;
1813 Shift_toCenterOfMagnetX = 0.039;
1814 Shift_toCenterOfMagnetY = 0.9; // 0.13;// drawing
1815 Shift_toCenterOfMagnetAX = 0.;
1816 Shift_toCenterOfMagnetAY = 0.0019; //
1817 X_shift = 0.02;
1818 Y_shift = 0.3;
1819
1820 } else {
1821 X_shift_seg[0] = 0.;
1822 X_shift_seg[1] = 0.;
1823 X_shift_seg[2] = 0.;
1824 X_shift_seg[3] = 0.;
1825 Shift_toCenterOfMagnetX = 0.;
1826 Shift_toCenterOfMagnetY = 0.;
1827 Shift_toCenterOfMagnetAX = 0.;
1828 Shift_toCenterOfMagnetAY = 0.;
1829 X_shift = 0.;
1830 Y_shift = 0.;
1831 }
1832
1833 //-- 5 detectors Shift---
1834 X_shiftUp[0] = 0.; //
1835 X_shiftUp[1] = 0.;
1836 X_shiftUp[2] = 0.;
1837 X_shiftUp[3] = 0.;
1838 X_shiftUp[4] = 0.; //-0.015;
1839
1840 Y_shiftUp[0] = 0.; //-0.005 + 0.02-0.005-0.008+0.03;
1841 Y_shiftUp[1] = 0.; //+0.03;
1842 Y_shiftUp[2] = 0.; //+0.006;
1843 Y_shiftUp[3] = 0.; //+.1;
1844 Y_shiftUp[4] = 0.; //+0.002;
1845
1846 return kSUCCESS;
1847}
1848//----------------------------------------------------------------------
1849
1850//---------------------Initialisation-----------------------------------
1851void BmnUpstreamTracking::PrepareArraysToProcessEvent()
1852{
1853 fBmnUpstreamTracksArray->Clear();
1854 OutVector.clear();
1855 vecUpTracks.clear();
1856 vec_points.clear();
1857 vec_twotracks.clear();
1858
1859 NSiTracks = 0;
1860 for (Int_t iPars = 0; iPars < kNumPars; iPars++) {
1861 for (Int_t j = 0; j < kBig; j++) {
1862 par_ab_SiTr[iPars][j] = -999.;
1863 par_SiTr_z[j] = -999.;
1864 }
1865 }
1866
1867 for (Int_t iCh = 0; iCh < kNumChambers; iCh++) {
1868 Nseg_Ch[iCh] = 0;
1869
1870 for (Int_t iPars = 0; iPars < kNumPars; iPars++) {
1871 for (Int_t j = 0; j < kBig; j++) {
1872 par_ab_Ch[iCh][iPars][j] = -999.;
1873 par_Seg_z[iCh][j] = -999.;
1874 }
1875 }
1876 }
1877 for (Int_t ip = 0; ip < kNumPairs; ip++) {
1878 NPCTracks[ip] = 0;
1879 for (Int_t iPars = 0; iPars < kNumPars; iPars++) {
1880 for (Int_t j = 0; j < kBig; j++) {
1881 par_ab_tr[ip][iPars][j] = -999.;
1882 }
1883 }
1884 }
1885 for (Int_t ist = 0; ist < kNumStSi; ist++) {
1886 NSiXYhits[ist] = 0;
1887 for (Int_t ib = 0; ib < kBig; ib++) {
1888 for (Int_t iPars = 0; iPars < kNumPars; iPars++) {
1889 SiXYhits[ist][ib][iPars] = -999.;
1890 }
1891 }
1892 }
1893 for (Int_t ib = 0; ib < kBig; ib++) {
1894 for (Int_t ip = 0; ip < kPoints; ip++) {
1895 for (Int_t i = 0; i < kPoin_Par; i++) {
1896 Points[ib][ip][i] = -999.;
1897 }
1898 }
1899 }
1900}
1901//----------------------------------------------------------------------
1902
1903BmnUpstreamTracking::BmnUpstreamTracking()
1904 : hAx_fitUp(nullptr)
1905 , hAy_fitUp(nullptr)
1906 , hx_fitUp(nullptr)
1907 , hy_fitUp(nullptr)
1908 , hchi2_fitUp(nullptr)
1909 , hNhitsUp(nullptr)
1910 , hdAx_tr_mc_comb(nullptr)
1911 , hdX_tr_mc_comb(nullptr)
1912 , hdAy_tr_mc_comb(nullptr)
1913 , hdY_tr_mc_comb(nullptr)
1914 , hdAx_uptr_mc(nullptr)
1915 , hdX_uptr_mc(nullptr)
1916 , hdAy_uptr_mc(nullptr)
1917 , hdY_uptr_mc(nullptr)
1918 , hDen_mcuptr(nullptr)
1919 , hNum_mcuptr(nullptr)
1920 , hEff_mcuptr(nullptr)
1921 , hAx_upmc(nullptr)
1922 , hAy_upmc(nullptr)
1923 , hX_upmc(nullptr)
1924 , hY_upmc(nullptr)
1925 , hNtr_reco(nullptr)
1926 , hNtr_mc(nullptr)
1927 , hNrecoTrif2mc(nullptr)
1928 , hAngle_reco(nullptr)
1929 , hAngle_recoifNmc2(nullptr)
1930 , hAngle_recoifNmc2Cases(nullptr)
1931 , hNtr_mc_vs_reco(nullptr)
1932 , hy_vs_x_Up(nullptr)
1933 , hY_vs_Xmctrue(nullptr)
1934 , hvertexXYUp(nullptr)
1935 , hTyTx_Up(nullptr)
1936{}
1937
1938//-------------------------constructor----------------------------------
1939BmnUpstreamTracking::BmnUpstreamTracking(Bool_t isExp, Int_t runNumber)
1940{
1941 fRunNumber = runNumber;
1942 expData = isExp;
1943 fInputBranchName1 = "BmnSiliconTrack";
1944 fInputBranchName2 = "BmnMwpcTrack";
1945 fInputBranchName3 = "BmnMwpcSegment";
1946 fInputBranchHits = "BmnSiliconHits";
1947 fOutputTracksBranchName = "BmnTrack";
1948 if (!expData) {
1949 fRunNumber = 0001;
1950 fInputBranchNameSimTrue = "BmnSiliconHitClean";
1951 fInputBranchName = "BmnMwpcHit";
1952 }
1953}
1954//----------------------------------------------------------------------
1955
1956//-------------------Destructor-----------------------------------------
1957BmnUpstreamTracking::~BmnUpstreamTracking()
1958{
1959
1960 for (Int_t iPars = 0; iPars < kNumPars; iPars++) {
1961 delete[] par_ab_SiTr[iPars];
1962 }
1963 for (Int_t iCh = 0; iCh < kNumChambers; iCh++) {
1964 for (Int_t iPars = 0; iPars < kNumPars; iPars++) {
1965 delete[] par_ab_Ch[iCh][iPars];
1966 }
1967 delete[] par_ab_Ch[iCh];
1968 delete[] par_Seg_z[iCh];
1969 }
1970 for (Int_t ip = 0; ip < kNumPairs; ip++) {
1971 for (Int_t iPars = 0; iPars < kNumPars; iPars++) {
1972 delete[] par_ab_tr[ip][iPars];
1973 }
1974 delete[] par_ab_tr[ip];
1975 delete[] par_ab_trz[ip];
1976 }
1977 for (Int_t ist = 0; ist < kNumStSi; ist++) {
1978 for (Int_t ib = 0; ib < kBig; ib++) {
1979 delete[] SiXYhits[ist][ib];
1980 }
1981 delete[] SiXYhits[ist];
1982 }
1983 for (Int_t ib = 0; ib < kBig; ib++) {
1984 for (Int_t ip = 0; ip < kPoints; ip++) {
1985 delete[] Points[ib][ip];
1986 }
1987 delete[] Points[ib];
1988 }
1989 for (Int_t ip = 0; ip < kPoints; ip++) {
1990 delete[] Xforglfit[ip];
1991 }
1992 for (Int_t i = 0; i < 4; i++) {
1993 delete[] Amatr[i];
1994 }
1995 delete[] X_shiftUp;
1996 delete[] AmatrInverted;
1997 delete[] rhs;
1998 delete[] AmatrArray;
1999 delete[] Amatr;
2000 delete[] line;
2001 delete[] Xforglfit;
2002 delete[] Points;
2003 delete[] NSiXYhits;
2004 delete[] SiXYhits;
2005 delete[] par_ab_trz;
2006 delete[] par_ab_tr;
2007 delete[] par_ab_SiTr;
2008 delete[] par_SiTr_z;
2009 delete[] par_ab_Ch;
2010 delete[] par_Seg_z;
2011 delete[] Nseg_Ch;
2012 delete[] X_shift_seg;
2013 delete[] Z_pair;
2014 delete[] Z_Chamber;
2015 delete[] NPCTracks;
2016 delete[] Y_shiftUp;
2017}
2018//----------------------------------------------------------------------
2019
2020//----------------------------------------------------------------------
2021void BmnUpstreamTracking::Finish()
2022{
2023
2024 if (fDebug) {
2025
2026 hEff_mcuptr->Divide(hNum_mcuptr, hDen_mcuptr, 1, 1);
2027 printf("BmnUpstreamTracking: write hists to file... ");
2028 fOutputFileName = Form("hUpstreamTracks_run%d.root", fRunNumber);
2029 cout << fOutputFileName << endl;
2030 TFile file(fOutputFileName, "RECREATE");
2031 fList.Write();
2032 file.Close();
2033 }
2034 cout << "Work time of the UpstreamTrack finder: " << workTime << " s" << endl;
2035}
2036//----------------------------------------------------------------------
fhList
TList fhList
Definition BmnUpstreamTracking.cxx:27
workTime
static Double_t workTime
Definition BmnUpstreamTracking.cxx:26
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
i
int i
Definition P4_F32vec4.h:22
m
__m128 m
Definition P4_F32vec4.h:27
Bool_t
bool Bool_t
Int_t
int Int_t
kFALSE
constexpr Bool_t kFALSE
Double_t
double Double_t
kTRUE
constexpr Bool_t kTRUE
Option_t
const char Option_t
Form
char * Form(const char *fmt,...)
BmnHit::GetIndex
Int_t GetIndex() const
Definition BmnHit.h:32
BmnHit::GetModule
Int_t GetModule()
Definition BmnHit.h:36
BmnHit::GetStation
Short_t GetStation() const
Definition BmnHit.h:35
BmnMwpcHit::GetPdgId
Int_t GetPdgId()
Definition BmnMwpcHit.h:54
BmnMwpcHit::GetMwpcId
Short_t GetMwpcId() const
Definition BmnMwpcHit.h:32
BmnMwpcHit::GetPlaneId
Int_t GetPlaneId()
Definition BmnMwpcHit.h:53
BmnMwpcHit::GetHitId
Int_t GetHitId() const
Definition BmnMwpcHit.h:34
BmnTrack::SetChi2
void SetChi2(Double_t chi2)
Definition BmnTrack.h:76
BmnTrack::SetParamLast
void SetParamLast(FairTrackParam &par)
Definition BmnTrack.h:72
BmnTrack::SetFlag
void SetFlag(Int_t flag)
Definition BmnTrack.h:74
BmnTrack::GetParamFirst
FairTrackParam * GetParamFirst()
Definition BmnTrack.h:63
BmnTrack::GetNHits
Int_t GetNHits() const
Accessors.
Definition BmnTrack.h:49
BmnTrack::SetNHits
void SetNHits(Int_t n)
Definition BmnTrack.h:80
BmnTrack::GetParamLast
FairTrackParam * GetParamLast()
Definition BmnTrack.h:65
BmnUpstreamTracking::vecUpTracks
vector< UpTracks > vecUpTracks
Definition BmnUpstreamTracking.h:183
BmnUpstreamTracking::PCTracksSiTracksMatching
void PCTracksSiTracksMatching(Double_t **, Double_t *, Int_t &, Double_t ***, Double_t **, Int_t *, vector< Smatch > &, vector< Smatch > &)
Definition BmnUpstreamTracking.cxx:732
BmnUpstreamTracking::OutVector
vector< Smatch > OutVector
Definition BmnUpstreamTracking.h:182
BmnUpstreamTracking::RecordCandidateSeg
void RecordCandidateSeg(Int_t &, Int_t &, Int_t &, Double_t **, Double_t ***, Double_t *, Double_t **, Double_t &, Double_t &, vector< Smatch > &)
Definition BmnUpstreamTracking.cxx:914
BmnUpstreamTracking::InvertMatrix
bool InvertMatrix(Double_t *, Double_t *)
Definition BmnUpstreamTracking.cxx:1478
BmnUpstreamTracking::RecordingBest
void RecordingBest(Int_t &, vector< Smatch > &, vector< Smatch > &)
Definition BmnUpstreamTracking.cxx:960
BmnUpstreamTracking::PrintAllTracks
void PrintAllTracks(vector< UpTracks > &)
Definition BmnUpstreamTracking.cxx:1392
BmnUpstreamTracking::vec_points
vector< MC_points > vec_points
Definition BmnUpstreamTracking.h:184
BmnUpstreamTracking::PCSegmentsSiTracksMatching
void PCSegmentsSiTracksMatching(Double_t **, Double_t *, Int_t &, Double_t ***, Double_t **, Int_t *, vector< Smatch > &, vector< Smatch > &)
Definition BmnUpstreamTracking.cxx:819
BmnUpstreamTracking::GetHitsToPoints
void GetHitsToPoints(vector< Smatch > &, Double_t ***, Int_t *, Double_t ***, Int_t *, Double_t ***, Int_t &)
Definition BmnUpstreamTracking.cxx:1109
BmnUpstreamTracking::vec_twotracks
vector< twotracks > vec_twotracks
Definition BmnUpstreamTracking.h:185
BmnUpstreamTracking::ReadSiliconTracks
void ReadSiliconTracks(Double_t **, Double_t *, Int_t &, vector< MC_points > &)
Definition BmnUpstreamTracking.cxx:334
BmnUpstreamTracking::ReadMWPCSegments
void ReadMWPCSegments(Double_t ***, Double_t **, Int_t *, vector< MC_points > &)
Definition BmnUpstreamTracking.cxx:476
BmnUpstreamTracking::ReadSiliconHits
void ReadSiliconHits(Double_t ***, Int_t *)
Definition BmnUpstreamTracking.cxx:449
BmnUpstreamTracking::compareSegments
static bool compareSegments(const Smatch &a, const Smatch &b)
Definition BmnUpstreamTracking.h:92
BmnUpstreamTracking::Exec
virtual void Exec(Option_t *opt)
Definition BmnUpstreamTracking.cxx:30
BmnUpstreamTracking::Init
virtual InitStatus Init()
Definition BmnUpstreamTracking.cxx:1577
BmnUpstreamTracking::ReadMWPCTracks
void ReadMWPCTracks(Double_t ***, Double_t **, Int_t *)
Definition BmnUpstreamTracking.cxx:684
BmnUpstreamTracking::MCefficiencyCalculation
void MCefficiencyCalculation(vector< MC_points > &, vector< UpTracks > &)
Definition BmnUpstreamTracking.cxx:91
BmnUpstreamTracking::GlobalFit
void GlobalFit(Double_t **, Double_t **, Double_t *)
Definition BmnUpstreamTracking.cxx:1527
BmnUpstreamTracking::fBmnUpstreamTracksArray
TClonesArray * fBmnUpstreamTracksArray
Output array of Silicon Hits.
Definition BmnUpstreamTracking.h:146
BmnUpstreamTracking::GetAddSiHits
void GetAddSiHits(vector< Smatch > &, Double_t ***, Int_t *)
Definition BmnUpstreamTracking.cxx:1061
BmnUpstreamTracking::TrackRecording
void TrackRecording(vector< UpTracks > &)
Definition BmnUpstreamTracking.cxx:1438
BmnUpstreamTracking::CalculationOfChi2
void CalculationOfChi2(Double_t ***, Int_t &, vector< UpTracks > &)
Definition BmnUpstreamTracking.cxx:1214
TFile::Close
void Close(Option_t *option="") override
sqrt
RVec< PromoteType< T > > sqrt(const RVec< T > &v)
TMath::ACos
Double_t ACos(Double_t)
TMath::RadToDeg
constexpr Double_t RadToDeg()
BmnUpstreamTracking::MC_points
Definition BmnUpstreamTracking.h:95
BmnUpstreamTracking::MC_points::x_si
Double_t x_si[4]
Definition BmnUpstreamTracking.h:113
BmnUpstreamTracking::MC_points::Id
Int_t Id
Definition BmnUpstreamTracking.h:96
BmnUpstreamTracking::MC_points::Np_Si
Int_t Np_Si
Definition BmnUpstreamTracking.h:99
BmnUpstreamTracking::MC_points::y_si
Double_t y_si[4]
Definition BmnUpstreamTracking.h:114
BmnUpstreamTracking::MC_points::np_3si
Bool_t np_3si
Definition BmnUpstreamTracking.h:110
BmnUpstreamTracking::MC_points::z_si
Double_t z_si[4]
Definition BmnUpstreamTracking.h:115
Smatch::params2
Double_t params2[4]
Definition BmnUpstreamTracking.h:59
Smatch::param1
Double_t param1[4]
Definition BmnUpstreamTracking.h:57
Smatch::param2
Double_t param2[4]
Definition BmnUpstreamTracking.h:58
Smatch::Chi2m
Double_t Chi2m
Definition BmnUpstreamTracking.h:56
Smatch::distX
Double_t distX
Definition BmnUpstreamTracking.h:54
Smatch::distY
Double_t distY
Definition BmnUpstreamTracking.h:55
UpTracks::CoordX
Double_t CoordX[5]
Definition BmnUpstreamTracking.h:73
UpTracks::CoordY
Double_t CoordY[5]
Definition BmnUpstreamTracking.h:74
UpTracks::CoordZ
Double_t CoordZ[5]
Definition BmnUpstreamTracking.h:75
UpTracks::param
Double_t param[4]
Definition BmnUpstreamTracking.h:72