BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
BmnMwpcTrackFinder.cxx
Go to the documentation of this file.
1// Author: Vasilisa Lenivenko <vasilisa@jinr.ru> 2018-07-18
2
3////////////////////////////////////////////////////////////////////////////////
4// //
5// BmnMwpcTrackFinder //
6// //
7// //
8// Implementation of an algorithm developed by //
9// Vasilisa Lenivenko and Vladimir Palchik //
10// to the BmnRoot software //
11// //
12// The algorithm serves for searching for track segments //
13// in the MWPC of the BM@N experiment //
14// //
15////////////////////////////////////////////////////////////////////////////////
16#include "BmnMwpcTrackFinder.h"
17
18#include "BmnMwpcSegment.h"
19#include "BmnMwpcTrack.h"
20#include "FairRootManager.h"
21#include "TFile.h"
22
23using namespace std;
24static Float_t workTime = 0.0;
25
26struct match
27{
28 Double_t Chi2m = 0.;
29 Int_t Ind1 = -1;
30 Int_t Ind2 = -1;
31 Int_t Nhits1 = 0;
32 Int_t Nhits2 = 0;
33 Double_t param1[4] = {999., 999., 999., 999.};
34 Double_t param2[4] = {999., 999., 999., 999.};
35 Float_t distX = -1.;
36 Float_t distY = -1.;
37 Float_t distAX = -1.;
38 Float_t distAY = -1.;
39};
40
41bool compareSegments(const match& a, const match& b)
42{
43 return a.Chi2m < b.Chi2m;
44}
45
46BmnMwpcTrackFinder::BmnMwpcTrackFinder(Bool_t isExp, Int_t runP, Int_t runNumber)
47 : expData(isExp)
48 , fEventNo(0)
49{
50 if (expData) {
51 fRunPeriod = runP;
52 fRunNumber = runNumber;
53
54 if (fRunPeriod == 6 || (fRunPeriod == 7 && fRunNumber > 3588)) // BMN
55 {
56 kNumPairs = 1;
57 kCh_max = 2;
58 } else if ((fRunPeriod == 7 && fRunNumber <= 3588) || (fRunPeriod == 8)) // SRC
59 {
60 kNumPairs = 2;
61 kCh_max = 4;
62 }
63 } else {
64 fInputBranchNameHit = "BmnMwpcHit";
65 kNumPairs = 2;
66 kCh_max = 4;
67 fRunPeriod = 7;
68 fRunNumber = 0001;
69 }
70
71 fMwpcGeo = new BmnMwpcGeometrySRC(fRunPeriod, fRunNumber);
72 fInputBranchName = "BmnMwpcSegment";
73 fOutputBranchName = "BmnMwpcTrack";
74 kBig = 100;
75 kCh_min = 0;
76}
77
78BmnMwpcTrackFinder::~BmnMwpcTrackFinder()
79{
80
81 for (Int_t iCh = 0; iCh < kNChambers; ++iCh) {
82 for (Int_t iPlane = 0; iPlane < kNPlanes; ++iPlane) {
83 delete[] XVU_Ch[iCh][iPlane];
84 delete[] Clust_Ch[iCh][iPlane];
85 }
86 for (Int_t i = 0; i < 4; ++i) {
87 delete[] par_ab_Ch[iCh][i];
88 }
89 delete[] XVU_Ch[iCh];
90 delete[] Clust_Ch[iCh];
91 delete[] par_ab_Ch[iCh];
92 delete[] shift[iCh];
93 delete[] kPln[iCh];
94 delete[] kZ_loc[iCh];
95 delete[] z_gl[iCh];
96 delete[] Chi2_ndf_Ch[iCh];
97 delete[] Nhits_Ch[iCh];
98 delete[] Nhits_match[iCh];
99 delete[] ind_best_Ch[iCh];
100 }
101
102 for (Int_t iPair = 0; iPair < kNumPairs; ++iPair) {
103 for (Int_t iBig = 0; iBig < kBig; ++iBig) {
104 for (Int_t i = 0; i < 4; ++i) {
105 delete[] sigma2_pair[iPair][iBig][i];
106 }
107 delete[] sigma2_pair[iPair][iBig];
108 }
109
110 for (Int_t i = 0; i < 4; ++i) {
111 delete[] par_ab_pair[iPair][i];
112 }
113
114 delete[] sigma2_pair[iPair];
115 delete[] par_ab_pair[iPair];
116 delete[] shift_pair[iPair];
117 delete[] Chi2_match_pair[iPair];
118 delete[] Chi2_ndf_pair[iPair];
119 delete[] ind_best_pair[iPair];
120 delete[] Nhits_pair[iPair];
121 delete[] sigma_delta[iPair];
122 }
123
124 for (Int_t ii = 0; ii < 4; ii++) {
125 delete[] Amatr[ii];
126 delete[] bmatr[ii];
127 }
128
129 delete[] Amatr;
130 delete[] bmatr;
131
132 delete[] kPln;
133 delete[] kZ_loc;
134 delete[] z_gl;
135 delete[] XVU_Ch;
136 delete[] Clust_Ch;
137 delete[] par_ab_Ch;
138 delete[] par_ab_pair;
139 delete[] sigma2_pair;
140 delete[] shift;
141 delete[] shift_pair;
142 delete[] kZ_midle_pair;
143
144 delete[] ChCent;
145 delete[] ZCh;
146 delete[] kZmid;
147 delete[] Chi2_match_pair;
148 delete[] Chi2_ndf_pair;
149 delete[] ind_best_pair;
150 delete[] Chi2_ndf_Ch;
151 delete[] Nhits_Ch;
152 delete[] Nhits_match;
153 delete[] Nhits_pair;
154 delete[] sigm2;
155 delete[] ipl;
156 delete[] Nbest_pair;
157 delete[] Nbest_Ch;
158 delete[] ind_best_Ch;
159 delete[] sigma_delta;
160}
161
162void BmnMwpcTrackFinder::Exec(Option_t* opt)
163{
164 if (!IsActive())
165 return;
166 clock_t tStart = clock();
167 PrepareArraysToProcessEvent();
168 if (fDebug)
169 cout << "\n======================== MWPC track finder exec started ===================\n" << endl;
170 if (fDebug)
171 cout << "Event number: " << fEventNo++ << endl;
172
173 ReadSegments(par_ab_Ch, Nhits_Ch, Chi2_ndf_Ch, Nbest_Ch, XVU_Ch, Clust_Ch, vec_points);
174
175 //--------Track-Segment matching between chambers--------------------
176 if (Nbest_Ch[0] > 0 && Nbest_Ch[1] > 0) {
177 if (fDebug)
178 cout << "Reco: N0= " << Nbest_Ch[0] << " N1= " << Nbest_Ch[1] << endl;
179 SegmentMatching(0, Nbest_Ch, par_ab_Ch, kZmid, ind_best_Ch, Nbest_pair, Chi2_match_pair, XVU_Ch, Nhits_Ch,
180 Nhits_match);
181 }
182
183 if (kNChambers > 2 && Nbest_Ch[2] > 0 && Nbest_Ch[3] > 0) {
184 if (fDebug)
185 cout << "Reco: N2= " << Nbest_Ch[2] << " N3= " << Nbest_Ch[3] << endl;
186 SegmentMatchingAfterTarget(2, Nbest_Ch, par_ab_Ch, kZmid, ind_best_Ch, Nbest_pair, Chi2_match_pair, XVU_Ch,
187 Nhits_Ch, Nhits_match);
188 }
189
190 if (fDebug) {
191 // printf("kNumPairs: %d\n", kNumPairs);
192 for (Int_t p = 0; p < kNumPairs; p++) {
193 for (Int_t se = 0; se < kmaxPairs; se++) {
194 if (Chi2_match_pair[p][se] != 999.) {
195 hChi2_match_pair.at(p)->Fill(Chi2_match_pair[p][se]);
196 // cout<<" pair "<<p<<" se "<<se<<" Chi2_match "<<Chi2_match_pair[p][se]<<endl;
197 }
198 }
199 // cout<<"Nbest_pair["<<p<<"]= "<<Nbest_pair[p]<<endl;
200 if (Nbest_pair[p]) {
201 hNbest_pair.at(p)->Fill(Nbest_pair[p]);
202 }
203 } // p
204 }
205 //--------------------------------------------------------------------
206
207 // ----------------Segment from 2-stations Fit-------------------------
208 Int_t First_Chamber = 0;
209 if (fDebug)
210 cout << " kNumPairs " << kNumPairs << " Nbest_pair[0] " << Nbest_pair[0] << " Nbest_pair[1] " << Nbest_pair[1]
211 << endl;
212 if (kNumPairs > 0 && Nbest_pair[0] > 0) {
213 if (fDebug)
214 cout << " Nbest_pair[0] " << Nbest_pair[0] << endl;
215 SegmentFit(First_Chamber, z_gl, sigm2, Nbest_pair, ind_best_Ch, par_ab_pair, Chi2_ndf_pair, XVU_Ch, Clust_Ch,
216 ind_best_pair, Nhits_Ch, Nhits_pair, sigma2_pair);
217 }
218
219 if (kNumPairs == 2 && Nbest_pair[1] > 0) {
220 // if (fDebug) cout<<" Nbest_pair[1] "<<Nbest_pair[1]<<endl;
221 SegmentFit(First_Chamber + 2, z_gl, sigm2, Nbest_pair, ind_best_Ch, par_ab_pair, Chi2_ndf_pair, XVU_Ch,
222 Clust_Ch, ind_best_pair, Nhits_Ch, Nhits_pair, sigma2_pair);
223 }
224 //--------------------------------------------------------------------
225 if (!expData)
226 MCefficiencyCalculation(vec_points, Nbest_pair, par_ab_pair);
227 //------------------Segment Parameters Alignment----------------------
228 for (Int_t iPair = 0; iPair < kNumPairs; iPair++) {
229 if (Nbest_pair[iPair] > 0)
230 SegmentParamAlignment(iPair, Nbest_pair, par_ab_pair, shift_pair);
231 }
232 //--------------------------------------------------------------------
233
234 //--------------------MWPC pairs matching ----------------------------
235 // if ( Nbest_pair[0] > 0 && Nbest_pair[1] > 0 ) PairMatching(Nbest_pair, par_ab_pair, kZ_midle_pair);
236 //--------------------------------------------------------------------
237
238 TracksStoring(Nbest_pair, par_ab_pair, Chi2_ndf_pair, Nhits_pair, sigma2_pair);
239
240 if (fDebug)
241 cout << "\n======================== MWPC track finder exec finished ==================" << endl;
242 clock_t tFinish = clock();
243 workTime += ((Float_t)(tFinish - tStart)) / CLOCKS_PER_SEC;
244} // Exec
245//----------------------------------------------------------------------------------
246
247void BmnMwpcTrackFinder::TracksStoring(Int_t* Nbest,
248 Double_t*** par_ab,
249 Double_t** Chi2_ndf,
250 Int_t** Nhits,
251 Double_t**** sigma2)
252{
253 Double_t theta, phi;
254 Float_t X_par_to_target, Y_par_to_target;
255 for (Int_t iPair = 0; iPair < kNumPairs; iPair++) {
256 if (Nbest[iPair] > 0) {
257 if (fDebug)
258 cout << " Nbest[" << iPair << "] " << Nbest[iPair] << endl;
259 for (Int_t itr = 0; itr < Nbest[iPair]; itr++) {
260 if (fDebug)
261 cout << " pair " << itr << " Chi2_ndf_pair " << Chi2_ndf[iPair][itr] << endl;
262
263 if (Chi2_ndf[iPair][itr] > 1000.)
264 continue;
265
266 if (fDebug) {
267 cout << " Chi2_ndf_pair " << Chi2_ndf[iPair][itr] << " Nhits_pair[" << iPair << "][" << itr << "] "
268 << Nhits[iPair][itr] << endl;
269 X_par_to_target = par_ab[iPair][0][itr] * (Z0_SRC - kZ_midle_pair[iPair]) + par_ab[iPair][1][itr];
270 Y_par_to_target = par_ab[iPair][2][itr] * (Z0_SRC - kZ_midle_pair[iPair]) + par_ab[iPair][3][itr];
271 phi = TMath::ATan2(par_ab[iPair][2][itr], par_ab[iPair][0][itr]); // phi = arctan(tgy/tgx)
272 theta = TMath::ATan2(par_ab[iPair][0][itr], TMath::Cos(phi)); // theta = arctan(tgx/cos(phi))
273 hpar_Ax_pair.at(iPair)->Fill(TMath::RadToDeg() * par_ab[iPair][0][itr]);
274 hpar_Bx_pair.at(iPair)->Fill(par_ab[iPair][1][itr]);
275 hpar_Ay_pair.at(iPair)->Fill(TMath::RadToDeg() * par_ab[iPair][2][itr]);
276 hpar_By_pair.at(iPair)->Fill(par_ab[iPair][3][itr]);
277 hYvsX_pair.at(iPair)->Fill(par_ab[iPair][1][itr], par_ab[iPair][3][itr]);
278
279 hpar_theta_pair.at(iPair)->Fill(TMath::RadToDeg() * theta);
280 hpar_phi_pair.at(iPair)->Fill(TMath::RadToDeg() * phi);
281 hX_in_target_pair.at(iPair)->Fill(X_par_to_target);
282 hY_in_target_pair.at(iPair)->Fill(Y_par_to_target);
283 hY_X_in_target_pair.at(iPair)->Fill(X_par_to_target, Y_par_to_target);
284 hAx_bx_in_target->Fill(X_par_to_target, TMath::RadToDeg() * par_ab[iPair][0][itr]);
285 hAy_by_in_target->Fill(Y_par_to_target, TMath::RadToDeg() * par_ab[iPair][2][itr]);
286 }
287
288 BmnTrack* Tr = new ((*fBmnMwpcTracksArray)[fBmnMwpcTracksArray->GetEntriesFast()]) BmnTrack();
289 Tr->SetChi2(Chi2_ndf[iPair][itr]);
290 Tr->SetNHits(Nhits[iPair][itr]);
291 Tr->SetFlag(itr);
292 FairTrackParam TrParams;
293 TrParams.SetPosition(TVector3(par_ab[iPair][1][itr], par_ab[iPair][3][itr], kZ_midle_pair[iPair]));
294 TrParams.SetTx(par_ab[iPair][0][itr]);
295 TrParams.SetTy(par_ab[iPair][2][itr]);
296
297 for (Int_t i1 = 0; i1 < 4; i1++) {
298 for (Int_t j1 = 0; j1 < 4; j1++) {
299 // if (fDebug) cout<<" sigma2["<<iPair<<"]["<<itr<<"]["<<i1<<"]["<<j1<<"]
300 // "<<sigma2[iPair][itr][i1][j1]<<endl;
301 TrParams.SetCovariance(i1, j1, sigma2[iPair][itr][i1][j1]);
302 }
303 }
304 Tr->SetParamFirst(TrParams);
305
306 } // Nbest_pair[iPair]
307
308 } //> 0
309 } // iPair
310}
311
312void BmnMwpcTrackFinder::MCefficiencyCalculation(vector<MC_points>& vec, Int_t* Nbest, Double_t*** par_ab)
313{
314
315 if (fDebug)
316 cout << " ---MC tracks association--" << endl;
317 // ax, bx, ay, by
318 Double_t delta2[4] = {-999., -999., -999., -999.};
319 Double_t sig[4] = {0.04, 0.08, 0.05, 0.08};
320
321 Double_t dmatch = 0.;
322 Double_t dmc_match[kMaxMC];
323 Int_t mc_tr_assoc[kMaxMC];
324
325 for (Int_t j = 0; j < kMaxMC; j++) {
326 dmc_match[j] = 1000.;
327 mc_tr_assoc[j] = -1;
328 }
329
330 Double_t dax = -999.;
331 Double_t day = -999.;
332 Double_t dx = -999.;
333 Double_t dy = -999.;
334 Int_t Ngood_mc_tracks = 0;
335 Int_t Ngood_reco_tracks = 0;
336
337 for (size_t itr = 0; itr < vec.size(); itr++) { // mc_tr
338
339 //---MC Eff ---
340 //---Den
341 if (fDebug && vec.at(itr).Np >= 8) { // min 4 + 4
342 // if ( vec.at(itr).xWas3 && vec.at(itr).uWas3 && vec.at(itr).vWas3 &&
343 // vec.at(itr).xWas2 && vec.at(itr).uWas2 && vec.at(itr).vWas2)
344 hDen_mctr->Fill(1);
345 hDen_mctr->Fill(0);
346 cout << " PCDen " << endl;
347 Ngood_mc_tracks++;
348 }
349
350 for (Int_t iseg = 0; iseg < Nbest[1]; iseg++) { // reco-tr
351 delta2[0] = vec.at(itr).param[0] - par_ab[1][0][iseg];
352 delta2[1] = vec.at(itr).param[1] - par_ab[1][1][iseg];
353 delta2[2] = vec.at(itr).param[2] - par_ab[1][2][iseg];
354 delta2[3] = vec.at(itr).param[3] - par_ab[1][3][iseg];
355
356 dmatch = 0.;
357 dmatch = (delta2[0] * delta2[0]) / (sig[0] * sig[0]) + (delta2[1] * delta2[1]) / (sig[1] * sig[1])
358 + (delta2[2] * delta2[2]) / (sig[2] * sig[2]) + (delta2[3] * delta2[3]) / (sig[3] * sig[3]);
359
360 if (dmc_match[itr] > dmatch) {
361 dmc_match[itr] = dmatch;
362 mc_tr_assoc[itr] = iseg;
363 dax = delta2[0];
364 dx = delta2[1];
365 day = delta2[2];
366 dy = delta2[3];
367 }
368 } // Nbest_pair[1]
369 if (fDebug) {
370 if (mc_tr_assoc[itr] > -1) {
371 cout << " index_mc " << itr << " Np " << vec.at(itr).Np << " mc_Id " << vec.at(itr).Id
372 << " index_track " << mc_tr_assoc[itr] <<
373 // " ax_mc "<<vec.at(itr).param[0]<<" reco_ind "<<mc_tr_assoc[itr]<<" ax "<<
374 // par_ab[1][0][mc_tr_assoc[itr]] <<
375 " dmc_match " << dmc_match[itr] << endl;
376 if (dax > -900.)
377 hdAx_tr_mc->Fill(dax);
378 if (dx > -900.)
379 hdX_tr_mc->Fill(dx);
380 if (day > -900.)
381 hdAy_tr_mc->Fill(day);
382 if (dy > -900.)
383 hdY_tr_mc->Fill(dy);
384 }
385 }
386
387 } // vec_points.size
388
389 if (fDebug)
390 cout << "reject poorly chosen association PC tracks " << endl;
391 for (size_t itr = 0; itr < vec.size(); itr++) { // mc_tr
392 if (mc_tr_assoc[itr] == -1)
393 continue;
394
395 for (size_t itr2 = 0; itr2 < vec.size(); itr2++) { // mc_tr
396 if (itr2 == itr)
397 continue;
398 if (mc_tr_assoc[itr2] == -1)
399 continue;
400
401 if (mc_tr_assoc[itr] == mc_tr_assoc[itr2]) {
402 if (dmc_match[itr2] > 20.)
403 mc_tr_assoc[itr2] = -1;
404 else {
405 // mc_tr_assoc[itr] = -1;
406 continue; // break;
407 }
408 }
409 } // itr2
410 //---MC Eff ---
411 //---Num
412 if (fDebug)
413 cout << " mc_Id " << vec.at(itr).Id << " assoc " << mc_tr_assoc[itr] << endl;
414 if (fDebug && mc_tr_assoc[itr] > -1 && vec.at(itr).Np >= 8) {
415 if (vec.at(itr).xWas3 && vec.at(itr).uWas3 && vec.at(itr).vWas3 && vec.at(itr).xWas2 && vec.at(itr).uWas2
416 && vec.at(itr).vWas2)
417 hNum_mctr->Fill(1);
418 hNum_mctr->Fill(0);
419 cout << " PCNum " << endl;
420 Ngood_reco_tracks++;
421 }
422 } // itr
423 if (fDebug) {
424 hNtrpc_mc->Fill(Ngood_mc_tracks);
425 hNtrpc_reco->Fill(Ngood_reco_tracks);
426 hNtrpc_mc_vs_reco->Fill(Ngood_reco_tracks, Ngood_mc_tracks);
427 cout << "PCp: Ngood_mc_tracks " << Ngood_mc_tracks << " Ngood_reco_tracks " << Ngood_reco_tracks << endl;
428 }
429}
430
431void BmnMwpcTrackFinder::ReadSegments(Double_t*** par_ab,
432 Int_t** Nhits,
433 Double_t** Chi2_ndf,
434 Int_t* Nbest,
435 Double_t*** XVU,
436 Double_t*** Clust,
437 vector<MC_points>& vec)
438{
439 if (fDebug)
440 cout << "--ReadSegments--" << endl;
441 int Npl_MC2[kMaxMC];
442 int Npl_MC3[kMaxMC];
443 if (!expData) {
444
445 Int_t mcTracksArray[kMaxMC];
446 Double_t X2mc[kMaxMC][kNPlanes];
447 Double_t Y2mc[kMaxMC][kNPlanes];
448 Double_t Z2mc[kMaxMC][kNPlanes];
449 Double_t X3mc[kMaxMC][kNPlanes];
450 Double_t Y3mc[kMaxMC][kNPlanes];
451 Double_t Z3mc[kMaxMC][kNPlanes];
452 for (Int_t Id = 0; Id < kMaxMC; Id++) {
453 Npl_MC2[Id] = 0;
454 Npl_MC3[Id] = 0;
455 mcTracksArray[Id] = -1;
456 for (Int_t i = 0; i < kNPlanes; i++) {
457 X2mc[Id][i] = -999.;
458 Y2mc[Id][i] = -999.;
459 Z2mc[Id][i] = -999.;
460 X3mc[Id][i] = -999.;
461 Y3mc[Id][i] = -999.;
462 Z3mc[Id][i] = -999.;
463 }
464 }
465
466 int tr_before = -1;
467 int Nmc_tracks = -1;
468
469 for (Int_t iMC = 0; iMC < fBmnHitsArray->GetEntriesFast(); ++iMC) {
470 BmnMwpcHit* hit = (BmnMwpcHit*)fBmnHitsArray->UncheckedAt(iMC);
471
472 Int_t st_MC = hit->GetMwpcId();
473 Int_t trackId_MC = hit->GetHitId();
474 Int_t pl_MC = hit->GetPlaneId();
475 // Short_t wire_MC = hit->GetWireNumber();
476 // Double_t time_MC = hit->GetWireTime();
477
478 // if (fDebug)cout<<" st_MC "<<st_MC<<" trackId_MC "<<trackId_MC<<" pl_MC "<<pl_MC<<" X "<<hit->GetX()<<"
479 // wire_MC "<<wire_MC<<endl;
480
481 if (tr_before != trackId_MC) {
482 Nmc_tracks++;
483 mcTracksArray[Nmc_tracks] = hit->GetHitId();
484 }
485 tr_before = trackId_MC;
486
487 if (st_MC == 2) {
488 X2mc[Nmc_tracks][pl_MC] = hit->GetX();
489 Y2mc[Nmc_tracks][pl_MC] = hit->GetY();
490 Z2mc[Nmc_tracks][pl_MC] = hit->GetZ();
491 Npl_MC2[Nmc_tracks]++;
492 }
493 if (st_MC == 3) {
494 X3mc[Nmc_tracks][pl_MC] = hit->GetX();
495 Y3mc[Nmc_tracks][pl_MC] = hit->GetY();
496 Z3mc[Nmc_tracks][pl_MC] = hit->GetZ();
497 Npl_MC3[Nmc_tracks]++;
498 }
499 // if (fDebug)cout<<" X2["<<Nmc_tracks<<"]["<<pl_MC<<"] "<<X2mc[Nmc_tracks][pl_MC]<<" Npl_MC2 "<<
500 // Npl_MC2[Nmc_tracks]<<" X3["<<Nmc_tracks<<"]["<<pl_MC<<"] "<<X3mc[Nmc_tracks][pl_MC]<<" Npl_MC3 "<<
501 // Npl_MC3[Nmc_tracks]<<endl;
502
503 } // iMC
504 Nmc_tracks++;
505
506 MC_points tmpTr;
507 if (fDebug)
508 cout << " Nmc_tracks " << Nmc_tracks << endl;
509 for (Int_t id = 0; id < Nmc_tracks; id++) {
510 if (fDebug)
511 cout << " id " << id << " Id_mc " << mcTracksArray[id] << " Npl2 " << Npl_MC2[id] << " Npl3 "
512 << Npl_MC3[id] << endl;
513 for (Int_t i = 0; i < 6; i++) {
514 if (fDebug && X2mc[id][i] > -900.)
515 cout << " ipl " << i << " X2 " << X2mc[id][i] << endl;
516 tmpTr.x2[i] = X2mc[id][i];
517 tmpTr.y2[i] = Y2mc[id][i];
518 tmpTr.z2[i] = Z2mc[id][i];
519 }
520 for (Int_t i = 0; i < 6; i++) {
521 if (fDebug && X3mc[id][i] > -900.)
522 cout << " ipl " << i << " X3 " << X3mc[id][i] << endl;
523 tmpTr.x3[i] = X3mc[id][i];
524 tmpTr.y3[i] = Y3mc[id][i];
525 tmpTr.z3[i] = Z3mc[id][i];
526 }
527 if (fDebug)
528 cout << endl;
529 tmpTr.Id = mcTracksArray[id];
530 tmpTr.Np2 = Npl_MC2[id];
531 tmpTr.Np3 = Npl_MC3[id];
532 tmpTr.Np = Npl_MC2[id] + Npl_MC3[id];
533
534 if (Npl_MC2[id] >= 4 || Npl_MC3[id] >= 4)
535 vec.push_back(tmpTr);
536 }
537
538 if (fDebug)
539 cout << " MC vec_points.size() " << vec.size() << endl;
540
541 // Double_t x_target_ch2, y_target_ch2, x_target_ch3, y_target_ch3;
542
543 for (size_t itr = 0; itr < vec.size(); itr++) {
544 // ch2
545 int i2 = -1;
546 if (vec.at(itr).x2[5] > -900.)
547 i2 = 5;
548 if (i2 < 0 && vec.at(itr).x2[4] > -900.)
549 i2 = 4;
550 if (i2 < 0 && vec.at(itr).x2[3] > -900.)
551 i2 = 3;
552 int i20 = 0;
553 if (vec.at(itr).x2[i20] < -900.)
554 i20 = 1;
555 if (vec.at(itr).x2[i20] < -900.)
556 i20 = 2;
557 if (i2 > -1 && vec.at(itr).Np2 > 3) {
558 vec.at(itr).param2[0] =
559 (vec.at(itr).x2[i20] - vec.at(itr).x2[i2]) / (vec.at(itr).z2[i20] - vec.at(itr).z2[i2]);
560 vec.at(itr).param2[1] = (vec.at(itr).x2[i20] + vec.at(itr).x2[i2]) * 0.5;
561 vec.at(itr).param2[2] =
562 (vec.at(itr).y2[i20] - vec.at(itr).y2[i2]) / (vec.at(itr).z2[i20] - vec.at(itr).z2[i2]);
563 vec.at(itr).param2[3] = (vec.at(itr).y2[i20] + vec.at(itr).y2[i2]) * 0.5;
564 }
565 // ch3
566 int i3 = -1;
567 if (vec.at(itr).x3[5] > -900.)
568 i3 = 5;
569 if (i3 < 0 && vec.at(itr).x3[4] > -900.)
570 i3 = 4;
571 if (i3 < 0 && vec.at(itr).x3[3] > -900.)
572 i3 = 3;
573 int i30 = 0;
574 if (vec.at(itr).x3[i30] < -900.)
575 i30 = 1;
576 if (vec.at(itr).x3[i30] < -900.)
577 i30 = 2;
578
579 if (i3 > -1 && vec.at(itr).Np3 > 3) {
580 vec.at(itr).param3[0] =
581 (vec.at(itr).x3[i30] - vec.at(itr).x3[i3]) / (vec.at(itr).z3[i30] - vec.at(itr).z3[i3]);
582 vec.at(itr).param3[1] = (vec.at(itr).x3[i30] + vec.at(itr).x3[i3]) * 0.5;
583 vec.at(itr).param3[2] =
584 (vec.at(itr).y3[i30] - vec.at(itr).y3[i3]) / (vec.at(itr).z3[i30] - vec.at(itr).z3[i3]);
585 vec.at(itr).param3[3] = (vec.at(itr).y3[i30] + vec.at(itr).y3[i3]) * 0.5;
586 }
587 // pair
588 if (i2 > -1 && vec.at(itr).Np2 > 3 && i3 > -1 && vec.at(itr).Np3 > 3) {
589 vec.at(itr).param[0] =
590 (vec.at(itr).x2[i20] - vec.at(itr).x3[i3]) / (vec.at(itr).z2[i20] - vec.at(itr).z3[i3]);
591 vec.at(itr).param[1] = (vec.at(itr).x2[i20] + vec.at(itr).x3[i3]) * 0.5;
592 vec.at(itr).param[2] =
593 (vec.at(itr).y2[i20] - vec.at(itr).y3[i3]) / (vec.at(itr).z2[i20] - vec.at(itr).z3[i3]);
594 vec.at(itr).param[3] = (vec.at(itr).y2[i20] + vec.at(itr).y3[i3]) * 0.5;
595 }
596 // triplet check
597 if (vec.at(itr).x2[0] > -900. || vec.at(itr).x2[3] > -900.)
598 vec.at(itr).xWas2 = 1;
599 if (vec.at(itr).x2[1] > -900. || vec.at(itr).x2[4] > -900.)
600 vec.at(itr).vWas2 = 1;
601 if (vec.at(itr).x2[2] > -900. || vec.at(itr).x2[5] > -900.)
602 vec.at(itr).uWas2 = 1;
603
604 if (vec.at(itr).x3[0] > -900. || vec.at(itr).x3[3] > -900.)
605 vec.at(itr).xWas3 = 1;
606 if (vec.at(itr).x3[1] > -900. || vec.at(itr).x3[4] > -900.)
607 vec.at(itr).vWas3 = 1;
608 if (vec.at(itr).x3[2] > -900. || vec.at(itr).x3[5] > -900.)
609 vec.at(itr).uWas3 = 1;
610
611 if (fDebug && vec.at(itr).Np > 7)
612 hYvsX_mc_pair->Fill(vec.at(itr).param[1], vec.at(itr).param[3]);
613 }
614 } // !expData
615
616 //----------------------- Read MWPC-Segmets---------------------------
617 for (Int_t iSegment = 0; iSegment < fBmnMwpcSegmentsArray->GetEntries(); iSegment++) {
618 BmnMwpcSegment* segment = (BmnMwpcSegment*)fBmnMwpcSegmentsArray->At(iSegment);
619 Int_t iCh;
620 Double_t Z = segment->GetParamFirst()->GetZ();
621 Int_t ise = segment->GetFlag(); // iSegmentID
622
623 if (ZCh[0] == Z)
624 iCh = 0;
625 if (ZCh[1] == Z)
626 iCh = 1;
627 if (ZCh[2] == Z)
628 iCh = 2;
629 if (ZCh[3] == Z)
630 iCh = 3;
631
632 Nhits[iCh][ise] = segment->GetNHits();
633 Chi2_ndf[iCh][ise] = segment->GetChi2();
634 par_ab[iCh][0][ise] = segment->GetParamFirst()->GetTx();
635 par_ab[iCh][1][ise] = segment->GetParamFirst()->GetX();
636 par_ab[iCh][2][ise] = segment->GetParamFirst()->GetTy();
637 par_ab[iCh][3][ise] = segment->GetParamFirst()->GetY();
638
639 for (Int_t i1 = 0; i1 < 6; i1++) {
640 XVU[iCh][i1][ise] = segment->GetCoord().at(i1);
641 Clust[iCh][i1][ise] = segment->GetClust().at(i1);
642 }
643 Nbest[iCh]++;
644 } // iSegment
645 //--------------------------------------------------------------------
646
647 for (Int_t iChamber = 0; iChamber < kNChambers; iChamber++) {
648 if (fDebug && Nbest[iChamber] > 0) {
649 for (Int_t ise = 0; ise < Nbest[iChamber]; ise++) {
650 hChi2_ndf_Ch.at(iChamber)->Fill(Chi2_ndf[iChamber][ise]);
651 hpar_Ax_Ch.at(iChamber)->Fill(par_ab[iChamber][0][ise]);
652 hpar_Bx_Ch.at(iChamber)->Fill(par_ab[iChamber][1][ise]);
653 hpar_Ay_Ch.at(iChamber)->Fill(par_ab[iChamber][2][ise]);
654 hpar_By_Ch.at(iChamber)->Fill(par_ab[iChamber][3][ise]);
655 // if (fDebug) cout<<" iChamber "<<iChamber<<" ax= "<<par_ab[iChamber][0][ise]<<" bx=
656 // "<<par_ab[iChamber][1][ise]<<" ay= "<<par_ab[iChamber][2][ise]<<" by= "<<par_ab[iChamber][3][ise]<<"
657 // Chi2 "<<Chi2_ndf[iChamber][ise]<<endl;
658 for (Int_t i1 = 0; i1 < 6; i1++) {
659 // if (fDebug) cout<<" Coord "<<i1 <<" " <<XVU[iChamber][i1][ise] <<endl;
660 }
661 }
662 } // if (fDebug)
663 } // iChamber
664}
665
666//------------------Segment Parameters Alignment----------------------
667void BmnMwpcTrackFinder::SegmentParamAlignment(Int_t chNum, Int_t* Nbest, Double_t*** par_ab, Float_t** shiftt)
668{
669 if (fDebug)
670 cout << endl;
671 for (Int_t iBest = 0; iBest < Nbest[chNum]; iBest++) {
672 // ax alpha ax^2
673 par_ab[chNum][0][iBest] +=
674 shiftt[chNum][0] + shiftt[chNum][0] * par_ab[chNum][0][iBest] * par_ab[chNum][0][iBest];
675 par_ab[chNum][2][iBest] +=
676 shiftt[chNum][2] + shiftt[chNum][2] * par_ab[chNum][2][iBest] * par_ab[chNum][2][iBest];
677 par_ab[chNum][1][iBest] += shiftt[chNum][1];
678 par_ab[chNum][3][iBest] += shiftt[chNum][3];
679 // if (fDebug) cout<<" chNum "<<chNum<<" after Alignment: iBest "<<iBest<<" Ax "<< par_ab[chNum][0][iBest]<<" bx
680 // "<< par_ab[chNum][1][iBest]<<" Ay "<< par_ab[chNum][1][iBest]<<" by "<< par_ab[chNum][3][iBest]<<endl;
681
682 } // iBest
683} // SegmentParamAlignment
684//----------------------------------------------------------------------
685
686//--------------Track Matching------------------------------------------
687void BmnMwpcTrackFinder::SegmentMatching(Int_t first_Ch,
688 Int_t* Nbest,
689 Double_t*** par_ab,
690 Float_t* Zmid,
691 Int_t** best_Ch,
692 Int_t* Nbest_pair_,
693 Double_t** Chi2_match_,
694 Double_t*** XVU_Ch_,
695 Int_t** Nhits_,
696 Int_t** Nhits_m)
697{
698 // cout<<" SegmentMatching "<<endl;
699 Float_t min_Chi2m = 100;
700 // Float_t min_distX = 99;
701 // Float_t min_distY = 99;
702 Double_t dx_loc = 99;
703 Double_t dy_loc = 99;
704 Float_t dAx12 = 0;
705 Float_t dAy12 = 0;
706 Float_t Min_distX[kmaxPairs];
707 Float_t Min_distY[kmaxPairs];
708
709 for (Int_t i = 0; i < kmaxPairs; i++) {
710 Min_distX[i] = -1;
711 Min_distY[i] = -1;
712 }
713
714 Int_t Nvariat = 0;
715 Int_t Nvariations = 100;
716
717 vector<match> vtmpSeg;
718 vector<match> OutVector;
719 match tmpSeg;
720 // match OutSegArray[kmaxPairs];
721
722 Int_t Pairr = 0;
723 if ((fRunPeriod == 7 && first_Ch == 2) || (fRunPeriod == 8))
724 Pairr = 1;
725 Int_t Secon_Ch = first_Ch + 1;
726
727 if (Nbest_Ch[first_Ch] > 0 && Nbest_Ch[Secon_Ch] > 0) {
728
729 for (Int_t bst1 = 0; bst1 < Nbest[first_Ch]; bst1++) {
730 // ch1 zloc0 -z_i
731 // Float_t x1mid = par_ab[first_Ch][0][bst1] *( 0 - kZmid[first_Ch]) + par_ab[first_Ch][1][bst1] ;
732 // Float_t y1mid = par_ab[first_Ch][2][bst1] *( 0 - kZmid[first_Ch]) + par_ab[first_Ch][3][bst1] ;
733
734 for (Int_t bst2 = 0; bst2 < Nbest[Secon_Ch]; bst2++) {
735 // ch2
736 // Float_t x2mid = par_ab[Secon_Ch][0][bst2] *( 0 - kZmid[Secon_Ch]) + par_ab[Secon_Ch][1][bst2] ;
737 // Float_t y2mid = par_ab[Secon_Ch][2][bst2] *( 0 - kZmid[Secon_Ch]) + par_ab[Secon_Ch][3][bst2] ;
738 dAx12 = par_ab[first_Ch][0][bst1] - par_ab[Secon_Ch][0][bst2];
739 dAy12 = par_ab[first_Ch][2][bst1] - par_ab[Secon_Ch][2][bst2];
740 // min_distX = x1mid - x2mid; //min
741 // min_distY = y1mid - y2mid; //min
742 dx_loc = par_ab[first_Ch][1][bst1] - par_ab[Secon_Ch][1][bst2];
743 dy_loc = par_ab[first_Ch][3][bst1] - par_ab[Secon_Ch][3][bst2];
744 // if (fDebug) cout<<" Pairr "<<Pairr<<" sigma ax "<<sigma_delta[Pairr][0]<<" ay
745 // "<<sigma_delta[Pairr][2]<<" sigma x "<<sigma_delta[Pairr][1]<<" y "<<sigma_delta[Pairr][3]<<endl;
746
747 Float_t Chi2_m = (dx_loc * dx_loc / (sigma_delta[Pairr][1] * sigma_delta[Pairr][1])
748 + dy_loc * dy_loc / (sigma_delta[Pairr][3] * sigma_delta[Pairr][3]));
749 // + dAx12*dAx12 /(sigma_delta[Pairr][0]*sigma_delta[Pairr][0]) + dAy12*dAy12
750 // /(sigma_delta[Pairr][2]*sigma_delta[Pairr][2]) ); //13.09.2019
751 if (fDebug) {
752 hdX_pair.at(Pairr)->Fill(dx_loc);
753 hdY_pair.at(Pairr)->Fill(dy_loc);
754 hdAx_pair.at(Pairr)->Fill(dAx12);
755 hdAy_pair.at(Pairr)->Fill(dAy12);
756 }
757 // if (fDebug)cout<<" Pairr "<<Pairr<<" bst1 "<<bst1<<" Nhits "<<Nhits_[first_Ch][bst1]<<" bst2
758 // "<<bst2<<" Nhits "<<Nhits_[Secon_Ch][bst2]<<" Chi2_m "<<Chi2_m<<endl;
759
760 if (Chi2_m < min_Chi2m && Nvariat < Nvariations && fabs(dx_loc) < 1.5 && fabs(dy_loc) < 1.5)
761 { // 09.2019
762
763 tmpSeg.Chi2m = Chi2_m;
764 tmpSeg.Ind1 = bst1;
765 tmpSeg.Ind2 = bst2;
766 tmpSeg.Nhits1 = Nhits_[first_Ch][bst1];
767 tmpSeg.Nhits2 = Nhits_[Secon_Ch][bst2];
768 tmpSeg.distX = dx_loc; // min_distX;
769 tmpSeg.distY = dy_loc; // min_distY;
770
771 for (int ipar = 0; ipar < 4; ipar++) {
772 tmpSeg.param1[ipar] = par_ab[first_Ch][ipar][bst1];
773 tmpSeg.param2[ipar] = par_ab[Secon_Ch][ipar][bst2];
774 }
775
776 vtmpSeg.push_back(tmpSeg);
777 Nvariat++;
778 } // if (Chi2_m
779
780 } // bst2++
781 } // bst1++
782
783 if (vtmpSeg.size() < 1)
784 return;
785
786 // vector sorting
787 sort(vtmpSeg.begin(), vtmpSeg.end(), compareSegments);
788 OutVector.clear();
789
790 // first best
791 OutVector.push_back(vtmpSeg.at(0));
792
793 // reject repeat index
794 Bool_t isMatch;
795 for (size_t iter = 0; iter < vtmpSeg.size(); ++iter) {
796 // printf("vtmpSeg.at(%d): %8.4f | %d - %d\n", iter, vtmpSeg.at(iter).Chi2m, vtmpSeg.at(iter).Ind1,
797 // vtmpSeg.at(iter).Ind2 );
798 isMatch = 0;
799 for (size_t InIter = 0; InIter < OutVector.size(); ++InIter) {
800 if (vtmpSeg.at(iter).Ind1 == OutVector.at(InIter).Ind1
801 || vtmpSeg.at(iter).Ind2 == OutVector.at(InIter).Ind2)
802 {
803 isMatch = 1;
804 continue;
805 }
806 }
807 // writing unique index
808 if (isMatch == 0)
809 OutVector.push_back(vtmpSeg.at(iter));
810 } // iter
811
812 if (fDebug && vtmpSeg.size() > 1)
813 hChi2best_Chi2fake_before_target->Fill(OutVector.at(0).Chi2m, vtmpSeg.at(1).Chi2m);
814
815 for (size_t iter = 0; iter < OutVector.size(); ++iter) {
816 // printf("OutVector.at(%d): %8.4f | %d - %d\n", iter, OutVector.at(iter).Chi2m, OutVector.at(iter).Ind1,
817 // OutVector.at(iter).Ind2);
818 if (Nbest_pair_[Pairr] < kmaxPairs) {
819 Chi2_match_[Pairr][Nbest_pair_[Pairr]] = OutVector.at(iter).Chi2m;
820 best_Ch[first_Ch][Nbest_pair_[Pairr]] = OutVector.at(iter).Ind1;
821 best_Ch[Secon_Ch][Nbest_pair_[Pairr]] = OutVector.at(iter).Ind2;
822 Nhits_m[first_Ch][Nbest_pair_[Pairr]] = OutVector.at(iter).Nhits1;
823 Nhits_m[Secon_Ch][Nbest_pair_[Pairr]] = OutVector.at(iter).Nhits2;
824 Min_distX[Nbest_pair_[Pairr]] = OutVector.at(iter).distX;
825 Min_distY[Nbest_pair_[Pairr]] = OutVector.at(iter).distY;
826
827 Nbest_pair_[Pairr]++;
828 } // < kmaxPairs){
829 } // iter
830
831 if (fDebug) {
832 for (Int_t ii = 0; ii < Nbest_pair_[Pairr]; ii++) {
833 FillEfficiency(first_Ch, XVU_Ch, Nhits_m, kMinHits, best_Ch[first_Ch][ii], Min_distX[ii],
834 Min_distY[ii]);
835 FillEfficiency(Secon_Ch, XVU_Ch, Nhits_m, kMinHits, best_Ch[Secon_Ch][ii], Min_distX[ii],
836 Min_distY[ii]);
837 }
838 }
839
840 } // if (Nbest_Ch[first_Ch] > 0 && Nbest_Ch[Secon_Ch] > 0)
841} // SegmentMatching
842//----------------------------------------------------------------------
843
844//-----------------Segment Matching after target------------------------
845void BmnMwpcTrackFinder::SegmentMatchingAfterTarget(Int_t first_Ch,
846 Int_t* Nbest,
847 Double_t*** par_ab,
848 Float_t* Zmid,
849 Int_t** best_Ch,
850 Int_t* Nbest_pair_,
851 Double_t** Chi2_match_,
852 Double_t*** XVU_Ch_,
853 Int_t** Nhits_,
854 Int_t** Nhits_m)
855{
856 if (fDebug)
857 cout << " SegmentMatching AfterTarget" << endl;
858 Float_t max_Chi2m = 100;
859 Float_t dAx = 0;
860 Float_t dAy = 0;
861 Double_t dx_loc = 10;
862 Double_t dy_loc = 10;
863 Float_t Min_distX[kmaxPairs];
864 Float_t Min_distY[kmaxPairs];
865
866 for (Int_t i = 0; i < kmaxPairs; i++) {
867 Min_distX[i] = -1;
868 Min_distY[i] = -1;
869 }
870
871 Int_t Nvariat = 0;
872 Int_t Nvariations = 100;
873 vector<match> vtmpSeg;
874 vector<match> OutVector;
875 match tmpSeg;
876 // match OutSegArray[kmaxPairs];
877
878 Int_t Pairr = 0; // doesn't work
879 if ((fRunPeriod == 7 && first_Ch == 2) || (fRunPeriod == 8))
880 Pairr = 1; // main stream
881 Int_t Secon_Ch = first_Ch + 1;
882
883 if (Nbest_Ch[first_Ch] > 0 && Nbest_Ch[Secon_Ch] > 0) {
884 if (fDebug)
885 hDen->Fill(0);
886 // if (fDebug)cout<<" Nbest[ "<<first_Ch<<"] = "<<Nbest_Ch[first_Ch]<<endl;
887 // if (fDebug)cout<<" Nbest[ "<<Secon_Ch<<"] = "<<Nbest_Ch[Secon_Ch]<<endl;
888
889 for (Int_t bst1 = 0; bst1 < Nbest[first_Ch]; bst1++) {
890 // ch1 zloc0 -z_i
891 // Float_t x1mid = par_ab[first_Ch][0][bst1] *( 0 - kZmid[first_Ch]) + par_ab[first_Ch][1][bst1] ;
892 // Float_t y1mid = par_ab[first_Ch][2][bst1] *( 0 - kZmid[first_Ch]) + par_ab[first_Ch][3][bst1] ;
893 // cout<<" bst1 " <<bst1<<" x1mid "<<x1mid<<" y1mid "<<y1mid<<endl;
894
895 for (Int_t bst2 = 0; bst2 < Nbest[Secon_Ch]; bst2++) {
896 // ch2
897 // Float_t x2mid = par_ab[Secon_Ch][0][bst2] *( 0 - kZmid[Secon_Ch]) + par_ab[Secon_Ch][1][bst2] ;
898 // Float_t y2mid = par_ab[Secon_Ch][2][bst2] *( 0 - kZmid[Secon_Ch]) + par_ab[Secon_Ch][3][bst2] ;
899 // if (fDebug)cout<<" bst2 " <<bst2<<" x2mid "<<x2mid<<" y2mid "<<y2mid<<endl;
900
901 dx_loc = par_ab[first_Ch][1][bst1] - par_ab[Secon_Ch][1][bst2];
902 dy_loc = par_ab[first_Ch][3][bst1] - par_ab[Secon_Ch][3][bst2];
903 dAx = par_ab[first_Ch][0][bst1] - par_ab[Secon_Ch][0][bst2];
904 dAy = par_ab[first_Ch][2][bst1] - par_ab[Secon_Ch][2][bst2];
905
906 Double_t Ax_23 =
907 (par_ab[Secon_Ch][1][bst2] - par_ab[first_Ch][1][bst1]) / (kZmid[Secon_Ch] - kZmid[first_Ch]);
908 Double_t Ay_23 =
909 (par_ab[Secon_Ch][3][bst2] - par_ab[first_Ch][3][bst1]) / (kZmid[Secon_Ch] - kZmid[first_Ch]);
910 Double_t x_target = Ax_23 * (Z0_SRC - ZCh[2]) + par_ab[first_Ch][1][bst1];
911 Double_t y_target = Ay_23 * (Z0_SRC - ZCh[2]) + par_ab[first_Ch][3][bst1];
912
913 // if (fDebug) cout<<" x_target "<< x_target<<" y_target "<< y_target<<endl;
914 // if (fDebug) cout<<" Pairr "<<Pairr<<" sigma ax "<<sigma_delta[Pairr][0]<<" ay
915 // "<<sigma_delta[Pairr][2]<<" sigma x "<<sigma_delta[Pairr][1]<<" y "<<sigma_delta[Pairr][3]<<endl;
916
917 Float_t Chi2_m = (dx_loc * dx_loc / (sigma_delta[Pairr][1] * sigma_delta[Pairr][1])
918 + dy_loc * dy_loc / (sigma_delta[Pairr][3] * sigma_delta[Pairr][3]));
919 //+ dAx*dAx /(sigma_delta[Pairr][0]*sigma_delta[Pairr][0]) + dAy*dAy
920 ///(sigma_delta[Pairr][2]*sigma_delta[Pairr][2]) ); //13.09.19
921 if (fDebug) {
922 hdX_pair.at(Pairr)->Fill(dx_loc); // Fill(min_distX);
923 hdY_pair.at(Pairr)->Fill(dy_loc); // Fill(min_distY);
924 hdAx_pair.at(Pairr)->Fill(dAx);
925 hdAy_pair.at(Pairr)->Fill(dAy);
926 }
927 // if (fDebug)cout<<" Pairr "<<Pairr<<" bst1 "<<bst1<<" Nhits "<<Nhits_[first_Ch][bst1]<<" bst2
928 // "<<bst2<<" Nhits "<<Nhits_[Secon_Ch][bst2]<<" Chi2_m "<<Chi2_m<<endl;
929
930 if (fDebug)
931 cout << " x_target " << x_target << " fabs(x_target - kX_target) " << fabs(x_target - kX_target)
932 << " fabs(y_target - kX_target) " << fabs(y_target - kX_target) << endl;
933 if (fDebug)
934 cout << " bst1 " << bst1 << " bst2 " << bst2 << " Chi2_m " << Chi2_m << endl;
935 if (Chi2_m < max_Chi2m && Nvariat < Nvariations && fabs(x_target - kX_target) < ktarget_region
936 && fabs(y_target - kY_target) < ktarget_regionY)
937 {
938
939 tmpSeg.Chi2m = Chi2_m;
940 tmpSeg.Ind1 = bst1;
941 tmpSeg.Ind2 = bst2;
942 tmpSeg.Nhits1 = Nhits_[first_Ch][bst1];
943 tmpSeg.Nhits2 = Nhits_[Secon_Ch][bst2];
944 tmpSeg.distX = dx_loc; // min_distX;
945 tmpSeg.distY = dy_loc; // min_distY;
946 tmpSeg.distAX = dAx;
947 tmpSeg.distAY = dAy;
948
949 for (int ipar = 0; ipar < 4; ipar++) {
950 tmpSeg.param1[ipar] = par_ab[first_Ch][ipar][bst1];
951 tmpSeg.param2[ipar] = par_ab[Secon_Ch][ipar][bst2];
952 }
953
954 vtmpSeg.push_back(tmpSeg);
955 Nvariat++;
956 } // if (Chi2_m
957 } // bst2++
958 } // bst1++
959
960 if (vtmpSeg.size() < 1)
961 return;
962
963 // vector sorting
964 sort(vtmpSeg.begin(), vtmpSeg.end(), compareSegments);
965 OutVector.clear();
966
967 Bool_t exist_pair[Nvariations];
968 for (size_t im = 0; im < vtmpSeg.size(); ++im)
969 exist_pair[im] = 1;
970
971 // cout<<" kmaxPairs "<<kmaxPairs<<endl;
972 for (size_t im = 0; im < vtmpSeg.size(); ++im) {
973 if (!exist_pair[im])
974 continue;
975 OutVector.push_back(vtmpSeg.at(im));
976 int InIter = OutVector.size() - 1;
977 if (fDebug)
978 cout << " im " << im << " InIter " << InIter << endl;
979
980 // reject repeat index
981 if (im + 1 < vtmpSeg.size()) {
982 for (size_t iter = im + 1; iter < vtmpSeg.size(); ++iter) {
983 if (fDebug)
984 cout << " iter " << iter << endl;
985 // printf("vtmpSeg.at(%d): %8.4f | %d - %d\n", iter, vtmpSeg.at(iter).Chi2m, vtmpSeg.at(iter).Ind1,
986 // vtmpSeg.at(iter).Ind2 );
987 if (!exist_pair[iter])
988 continue;
989 if (fDebug)
990 cout << " vtmpSeg.at(iter).Ind1 " << vtmpSeg.at(iter).Ind1 << " OutVector.at(InIter).Ind1 "
991 << OutVector.at(InIter).Ind1 << endl;
992 if (vtmpSeg.at(iter).Ind1 == OutVector.at(InIter).Ind1
993 || vtmpSeg.at(iter).Ind2 == OutVector.at(InIter).Ind2)
994 exist_pair[iter] = 0; // VVP 17.07.2019 -> for matching with silicon!!!
995 } // iter
996 } // if ( im + 1 < vtmpSeg.size())
997 } // im
998
999 if (fDebug)
1000 cout << " OutVector.at(0).Chi2m " << OutVector.at(0).Chi2m << " vtmpSeg.at(0).Chi2m " << vtmpSeg.at(0).Chi2m
1001 << endl;
1002 if (fDebug && vtmpSeg.size() > 1)
1003 hChi2best_Chi2fake_after_target->Fill(OutVector.at(0).Chi2m, vtmpSeg.at(1).Chi2m);
1004
1005 for (size_t iter = 0; iter < OutVector.size(); ++iter) {
1006 if (fDebug)
1007 printf("OutVector.at(%zu): %8.4f | %d - %d\n", iter, OutVector.at(iter).Chi2m, OutVector.at(iter).Ind1,
1008 OutVector.at(iter).Ind2);
1009 if (Nbest_pair_[Pairr] < kmaxPairs) {
1010
1011 Chi2_match_[Pairr][Nbest_pair_[Pairr]] = OutVector.at(iter).Chi2m;
1012 best_Ch[first_Ch][Nbest_pair_[Pairr]] = OutVector.at(iter).Ind1;
1013 best_Ch[Secon_Ch][Nbest_pair_[Pairr]] = OutVector.at(iter).Ind2;
1014 Nhits_m[first_Ch][Nbest_pair_[Pairr]] = OutVector.at(iter).Nhits1;
1015 Nhits_m[Secon_Ch][Nbest_pair_[Pairr]] = OutVector.at(iter).Nhits2;
1016 Min_distX[Nbest_pair_[Pairr]] = OutVector.at(iter).distX;
1017 Min_distY[Nbest_pair_[Pairr]] = OutVector.at(iter).distY;
1018
1019 Double_t Ax = (OutVector.at(iter).param2[1] - OutVector.at(iter).param1[1]) / (ZCh[3] - ZCh[2]);
1020 Double_t Xv = Ax * (Z0_SRC - ZCh[2]) + OutVector.at(iter).param1[1];
1021 Double_t Ay = (OutVector.at(iter).param2[3] - OutVector.at(iter).param1[3]) / (ZCh[3] - ZCh[2]);
1022 Double_t Yv = Ay * (Z0_SRC - ZCh[2]) + OutVector.at(iter).param1[3];
1023 if (fDebug) {
1024 hdX_pair_1->Fill(OutVector.at(iter).distX);
1025 hdY_pair_1->Fill(OutVector.at(iter).distY);
1026 hdAx_pair_1->Fill(OutVector.at(iter).distAX);
1027 hdAy_pair_1->Fill(OutVector.at(iter).distAY);
1028 hChi2m_pair_1->Fill(OutVector.at(iter).Chi2m);
1029 hXv_pair_1->Fill(Xv);
1030 hYv_pair_1->Fill(Yv);
1031 }
1032
1033 Nbest_pair_[Pairr]++;
1034 } // < kmaxPairs){
1035 } // iter
1036
1037 if (fDebug) {
1038 for (Int_t ii = 0; ii < Nbest_pair_[Pairr]; ii++) {
1039 FillEfficiency(first_Ch, XVU_Ch, Nhits_m, kMinHits, best_Ch[first_Ch][ii], Min_distX[ii],
1040 Min_distY[ii]);
1041 FillEfficiency(Secon_Ch, XVU_Ch, Nhits_m, kMinHits, best_Ch[Secon_Ch][ii], Min_distX[ii],
1042 Min_distY[ii]);
1043 }
1044 }
1045 if (fDebug)
1046 cout << "--Matching efficiency--" << endl;
1047 if (fDebug)
1048 cout << " Nbest_pair_[" << Pairr << "] " << Nbest_pair_[Pairr] << endl;
1049 if (fDebug && Nbest_pair_[Pairr] > 0)
1050 hNum->Fill(0);
1051
1052 } // if (Nbest_Ch[first_Ch] > 0 && Nbest_Ch[Secon_Ch] > 0)
1053} // SegmentMatching
1054//----------------------------------------------------------------------
1055
1056//------------------Pair Matching---------------------------------------
1057void BmnMwpcTrackFinder::PairMatching(Int_t* Nbest_p, Double_t*** par_ab_p, Float_t* kZ_midle_pair_)
1058{
1059 Int_t Npair_dist = 0;
1060 Double_t sig_dx = 1., sig_dy = 1., sig_dax = .006, sig_day = .006;
1061 Double_t min_Chi2m = 50.;
1062 Double_t dXm, dYm, dAx12m, dAy12m;
1063 Float_t X_pair0_to_pair1, Y_pair0_to_pair1;
1064 Double_t phi0, theta0, phi1, theta1;
1065
1066 for (Int_t pair0 = 0; pair0 < Nbest_p[0]; pair0++) {
1067
1068 // pair 0 zloc0 -z_i
1069 Double_t x1 = par_ab_p[0][0][pair0] * (Z0_SRC - kZ_midle_pair_[0]) + par_ab_p[0][1][pair0];
1070 Double_t y1 = par_ab_p[0][2][pair0] * (Z0_SRC - kZ_midle_pair_[0]) + par_ab_p[0][3][pair0];
1071 // if (fDebug) cout<<" pair "<<0<<" pair0 "<<pair0<<" ax "<<par_ab_p[0][0][pair0]<<" bx "<<
1072 // par_ab_p[0][1][pair0]<<" ay "<<par_ab_p[0][2][pair0]<<" by "<<par_ab_p[0][3][pair0]<<" Z0_SRC "<<Z0_SRC<<"
1073 // kZ_midle_pair "<<kZ_midle_pair_[0]<<endl;
1074
1075 for (Int_t pair1 = 0; pair1 < Nbest_p[1]; pair1++) {
1076
1077 // pair 1
1078 Double_t x2 = par_ab_p[1][0][pair1] * (Z0_SRC - kZ_midle_pair_[1]) + par_ab_p[1][1][pair1];
1079 Double_t y2 = par_ab_p[1][2][pair1] * (Z0_SRC - kZ_midle_pair_[1]) + par_ab_p[1][3][pair1];
1080 // if (fDebug) cout<<" pair "<<1<<" pair1 "<<pair1<<" ax "<<par_ab_p[1][0][pair1]<<" bx "<<
1081 // par_ab_p[1][1][pair1]<<" ay "<<par_ab_p[1][2][pair1]<<" by "<<par_ab_p[1][3][pair1]<<" Z0_SRC
1082 // "<<Z0_SRC<<" kZ_midle_pair "<<kZ_midle_pair_[1]<<endl;
1083 Float_t dAx12 = par_ab_p[0][0][pair0] - par_ab_p[1][0][pair1];
1084 Float_t dAy12 = par_ab_p[0][2][pair0] - par_ab_p[1][2][pair1];
1085 Float_t dX = x1 - x2;
1086 Float_t dY = y1 - y2;
1087
1088 Float_t Chi2_m = dX / (sig_dx * sig_dx) + dY / (sig_dy * sig_dy) + dAx12 * dAx12 / (sig_dax * sig_dax)
1089 + dAy12 * dAy12 / (sig_day * sig_day);
1090 // Pair
1091 X_pair0_to_pair1 = par_ab_p[0][0][pair0] * (kZ_midle_pair[0] - kZ_midle_pair[1]) + par_ab_p[0][1][pair0];
1092 Y_pair0_to_pair1 = par_ab_p[0][2][pair0] * (kZ_midle_pair[0] - kZ_midle_pair[1]) + par_ab_p[0][3][pair0];
1093
1094 if (Chi2_m < min_Chi2m) {
1095 min_Chi2m = Chi2_m; // min
1096 dXm = dX;
1097 dYm = dY;
1098 dAx12m = dAx12;
1099 dAy12m = dAy12;
1100 Npair_dist++;
1101
1102 if (fDebug) {
1103 hdX_pair01_vs_x1->Fill(x2, dXm);
1104 hdY_pair01_vs_y1->Fill(y2, dYm);
1105 hdX_pair01_inZpair1->Fill(X_pair0_to_pair1 - par_ab_p[1][1][pair1]);
1106 hdY_pair01_inZpair1->Fill(Y_pair0_to_pair1 - par_ab_p[1][3][pair1]);
1107 // hY_X_in_target_pair.at(0) -> Fill(x1, y1);
1108 // hY_X_in_target_pair.at(1) -> Fill(x2, y2);
1109 hAx_bx_in_target_pair.at(0)->Fill(x1, TMath::RadToDeg() * par_ab_p[0][0][pair0]);
1110 hAx_bx_in_target_pair.at(1)->Fill(x2, TMath::RadToDeg() * par_ab_p[1][0][pair1]);
1111 hAy_by_in_target_pair.at(0)->Fill(y1, TMath::RadToDeg() * par_ab_p[0][2][pair0]);
1112 hAy_by_in_target_pair.at(1)->Fill(y2, TMath::RadToDeg() * par_ab_p[1][2][pair1]);
1113 }
1114
1115 phi0 = TMath::ATan2(par_ab_pair[0][2][pair0], par_ab_pair[0][0][pair0]); // phi = arctan(tgy/tgx)
1116 theta0 = TMath::ATan2(par_ab_pair[0][0][pair0], TMath::Cos(phi0)); // theta = arctan(tgx/cos(phi))
1117 phi1 = TMath::ATan2(par_ab_pair[1][2][pair1], par_ab_pair[1][0][pair1]); // phi = arctan(tgy/tgx)
1118 theta1 = TMath::ATan2(par_ab_pair[1][0][pair1], TMath::Cos(phi1)); // theta = arctan(tgx/cos(phi))
1119
1120 if (fDebug)
1121 htheta_p1vsp0->Fill(TMath::RadToDeg() * theta0, TMath::RadToDeg() * theta1);
1122 }
1123
1124 } // pair1
1125 } // pair0
1126
1127 if (fDebug)
1128 hChi2_m_target->Fill(min_Chi2m);
1129 if (fDebug && min_Chi2m < 50) {
1130 hdX_target->Fill(dXm);
1131 hdY_target->Fill(dYm);
1132 hdAx_target->Fill(dAx12m);
1133 hdAy_target->Fill(dAy12m);
1134 }
1135
1136} // PairMatching
1137//----------------------------------------------------------------------
1138
1139//-------------Segment Fit----------------------------------------------
1140void BmnMwpcTrackFinder::SegmentFit(Int_t First_Ch,
1141 Float_t** z_gl_,
1142 Float_t* sig,
1143 Int_t* Nbest_pair_,
1144 Int_t** ind_best_Ch_,
1145 Double_t*** par_ab_pair_,
1146 Double_t** Chi2_ndf_pair_,
1147 Double_t*** XVU_Ch_,
1148 Double_t*** Clust_Ch_,
1149 Int_t** ind_best_pair_,
1150 Int_t** Nhits_Ch_,
1151 Int_t** Nhits_pair_,
1152 Double_t**** sigma2_p)
1153{
1154 int chiEl = 0;
1155 if (First_Ch == 2)
1156 chiEl = 1;
1157 if (Nbest_pair_[chiEl] >= kmaxPairs) {
1158 printf("!!! ERROR: Nbest_pair_[%d] > 10\n", chiEl);
1159 return;
1160 }
1161
1162 Double_t sigm_1[kNPlanes], sigm2_1[kNPlanes];
1163 Double_t sigm_2[kNPlanes], sigm2_2[kNPlanes];
1164 Int_t Pair1 = 0;
1165 if ((fRunPeriod == 7 && First_Ch == 2) || (fRunPeriod == 8))
1166 Pair1 = 1;
1167
1168 for (Int_t bst = 0; bst < Nbest_pair_[Pair1]; bst++) {
1169 if (fDebug)
1170 cout << "SegmentFit bst " << bst << endl;
1171
1172 Int_t fir = First_Ch;
1173 Int_t sec = First_Ch + 1;
1174 Int_t best1 = ind_best_Ch_[fir][bst];
1175 Int_t best2 = ind_best_Ch_[sec][bst];
1176
1177 int h1[6] = {0, 0, 0, 0, 0, 0};
1178 int h2[6] = {0, 0, 0, 0, 0, 0};
1179
1180 for (Int_t i = 0; i < 6; i++) {
1181 sigm_1[i] = 1.;
1182 sigm_2[i] = 1.;
1183 sigm2_1[i] = 1.;
1184 sigm2_2[i] = 1.;
1185 h1[i] = 0;
1186 h2[i] = 0;
1187
1188 if (XVU_Ch_[fir][i][best1] > -900.) {
1189 h1[i] = 1;
1190 sigm_1[i] = (Clust_Ch_[fir][i][best1] * dw) / sq12;
1191 sigm2_1[i] = sigm_1[i] * sigm_1[i];
1192 } // if coord was
1193
1194 if (XVU_Ch_[sec][i][best2] > -900.) {
1195 h2[i] = 1;
1196 sigm_2[i] = (Clust_Ch_[sec][i][best2] * dw) / sq12;
1197 sigm2_2[i] = sigm_2[i] * sigm_2[i];
1198 } // if coord was
1199 } // i6
1200
1201 for (Int_t im = 0; im < 4; im++) {
1202 for (Int_t ii = 0; ii < 4; ii++) {
1203 Amatr[im][ii] = 0.;
1204 bmatr[im][ii] = 0.;
1205 }
1206 }
1207
1208 FillFitMatrix(fir, Amatr, z_gl_, sigm2_1, h1);
1209 FillFitMatrix(sec, Amatr, z_gl_, sigm2_2, h2);
1210
1211 Double_t matrF[4] = {0, 0, 0, 0}; // free coef
1212
1213 FillFreeCoefVector(fir, matrF, XVU_Ch, best1, z_gl_, sigm2_1, h1);
1214 FillFreeCoefVector(sec, matrF, XVU_Ch, best2, z_gl_, sigm2_2, h2);
1215
1216 if (fDebug) {
1217 for (Int_t i = 0; i < kNPlanes; i++) {
1218 // cout<<" h1 "<<h1[i]<<" XVU_Ch1 "<<XVU_Ch_[fir][i][best1]<<" sigm2 "<< sigm2_1[i]<<endl;
1219 //<<" z "<<z_gl_[fir][i]<<endl;
1220 }
1221 cout << endl;
1222 for (Int_t i = 0; i < 6; i++) {
1223 // cout<<" h2 "<<h2[i]<<" XVU_Ch2 "<<XVU_Ch_[sec][i][best2]<<" sigm2 "<< sigm2_2[i]<<endl;
1224 //<<" z "<<z_gl_[sec][i]<<endl;
1225 }
1226 }
1227 // Gaussian algorithm for 4x4 matrix inversion
1228 Double_t A0matr[4][4];
1229 for (Int_t i1 = 0; i1 < 4; i1++) {
1230 for (Int_t j1 = 0; j1 < 4; j1++) {
1231 A0matr[i1][j1] = Amatr[i1][j1];
1232 // if (fDebug) cout<<" Amatr["<<i1<<"]["<<j1<<"] "<<Amatr[i1][j1]<<endl;
1233 }
1234 }
1235
1236 InverseMatrix(Amatr, bmatr);
1237 Double_t sum;
1238 // Double_t A1[4][4] = {{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}};
1239
1240 for (Int_t i1 = 0; i1 < 4; ++i1)
1241 for (Int_t j1 = 0; j1 < 4; ++j1) {
1242 sum = 0;
1243 for (Int_t k1 = 0; k1 < 4; ++k1) {
1244 Double_t a0 = A0matr[i1][k1];
1245 Double_t b0 = bmatr[k1][j1];
1246 sum += a0 * b0;
1247 // A1[i1][j1] = sum;
1248 }
1249 }
1250
1251 for (Int_t i1 = 0; i1 < 4; i1++) {
1252 par_ab_pair_[Pair1][i1][bst] = 0;
1253 for (Int_t j1 = 0; j1 < 4; j1++) {
1254 par_ab_pair_[Pair1][i1][bst] += bmatr[i1][j1] * matrF[j1];
1255 // if (fDebug)cout<<"bmatr["<<i1<<"]["<<j1<<"] "<<bmatr[i1][j1]<<endl;
1256 }
1257 } // i1
1258
1259 Double_t Xtarget = par_ab_pair_[Pair1][0][bst] * (Z0_SRC - kZ_midle_pair[Pair1])
1260 + par_ab_pair_[Pair1][1][bst]; //? + shift_pair[Pair1][1];
1261 Double_t Ytarget = par_ab_pair_[Pair1][2][bst] * (Z0_SRC - kZ_midle_pair[Pair1])
1262 + par_ab_pair_[Pair1][3][bst]; //? + shift_pair[Pair1][3];
1263
1264 if (fDebug) {
1265 cout << endl;
1266 cout << " Xtarget " << Xtarget << " Ytarget " << Ytarget << endl;
1267 cout << " Pair " << Pair1 << " par [0] " << par_ab_pair_[Pair1][0][bst] << " par [1] "
1268 << par_ab_pair_[Pair1][1][bst] << " par [2] " << par_ab_pair_[Pair1][2][bst] << " par [3] "
1269 << par_ab_pair_[Pair1][3][bst] << endl;
1270 cout << endl;
1271 }
1272
1273 Float_t dx_[kNPlanes];
1274 Chi2_ndf_pair_[Pair1][bst] = 0;
1275
1276 for (Int_t i1 = 0; i1 < kNPlanes; i1++) {
1277 dx_[i1] = 0.;
1278
1279 if (XVU_Ch_[fir][i1][best1] > -900.) {
1280 if (i1 == 0 || i1 == 3)
1281 dx_[i1] = XVU_Ch_[fir][i1][best1] - par_ab_pair_[Pair1][0][bst] * z_gl_[fir][i1]
1282 - par_ab_pair_[Pair1][1][bst];
1283 if (i1 == 2 || i1 == 5)
1284 dx_[i1] = XVU_Ch_[fir][i1][best1]
1285 - 0.5 * (par_ab_pair_[Pair1][0][bst] + sq3 * par_ab_pair_[Pair1][2][bst]) * z_gl_[fir][i1]
1286 - 0.5 * (par_ab_pair_[Pair1][1][bst] + sq3 * par_ab_pair_[Pair1][3][bst]);
1287 if (i1 == 1 || i1 == 4)
1288 dx_[i1] = XVU_Ch_[fir][i1][best1]
1289 - 0.5 * (par_ab_pair_[Pair1][0][bst] - sq3 * par_ab_pair_[Pair1][2][bst]) * z_gl_[fir][i1]
1290 - 0.5 * (par_ab_pair_[Pair1][1][bst] - sq3 * par_ab_pair_[Pair1][3][bst]);
1291 Chi2_ndf_pair_[Pair1][bst] = Chi2_ndf_pair_[Pair1][bst] + dx_[i1] * dx_[i1] / (sigm2_1[i1]);
1292 // if ( fDebug) cout<<"best1 "<<best1 <<" i1 "<<i1<<" dx_ "<<dx_[i1]<<" XVU_Ch1
1293 // "<<XVU_Ch_[fir][i1][best1]<<" Chi2_ndf_Ch1_2 "<<Chi2_ndf_pair_[Pair1][bst]<<" z_gl1
1294 // "<<z_gl_[fir][i1]<<endl;
1295
1296 } // if( Wires_Ch1[i1][best2]>-1){
1297 }
1298
1299 for (Int_t i2 = 0; i2 < kNPlanes; i2++) {
1300
1301 if (XVU_Ch_[sec][i2][best2] > -900.) { // if(Wires_Ch2_[i2][best2]>-1){
1302 if (i2 == 0 || i2 == 3)
1303 dx_[i2] = XVU_Ch_[sec][i2][best2] - par_ab_pair_[Pair1][0][bst] * z_gl_[sec][i2]
1304 - par_ab_pair_[Pair1][1][bst];
1305 if (i2 == 2 || i2 == 5)
1306 dx_[i2] = XVU_Ch_[sec][i2][best2]
1307 - 0.5 * (par_ab_pair_[Pair1][0][bst] + sq3 * par_ab_pair_[Pair1][2][bst]) * z_gl_[sec][i2]
1308 - 0.5 * (par_ab_pair_[Pair1][1][bst] + sq3 * par_ab_pair_[Pair1][3][bst]);
1309 if (i2 == 1 || i2 == 4)
1310 dx_[i2] = XVU_Ch_[sec][i2][best2]
1311 - 0.5 * (par_ab_pair_[Pair1][0][bst] - sq3 * par_ab_pair_[Pair1][2][bst]) * z_gl_[sec][i2]
1312 - 0.5 * (par_ab_pair_[Pair1][1][bst] - sq3 * par_ab_pair_[Pair1][3][bst]);
1313 Chi2_ndf_pair_[Pair1][bst] = Chi2_ndf_pair_[Pair1][bst] + dx_[i2] * dx_[i2] / (sigm2_2[i2]);
1314 // if ( fDebug) cout<<"best2 "<<best2 <<" i2 "<<i2<<" dx_ "<<dx_[i2]<<" XVU_Ch2
1315 // "<<XVU_Ch_[sec][i2][best2]<<" Chi2_ndf_Ch1_2 "<<Chi2_ndf_pair_[Pair1][bst]<<" z_gl2
1316 // "<<z_gl_[sec][i2]<<endl;
1317 } // if( Wires_Ch2[i2][best2]>-1){
1318 }
1319
1320 if (Nhits_Ch_[fir][best1] + Nhits_Ch_[sec][best2] > 4)
1321 Chi2_ndf_pair_[Pair1][bst] =
1322 Chi2_ndf_pair_[Pair1][bst] / (Nhits_Ch_[fir][best1] + Nhits_Ch_[sec][best2] - 4);
1323 if (fDebug)
1324 hChi2_ndf_pair.at(Pair1)->Fill(Chi2_ndf_pair_[Pair1][bst]);
1325 if (fDebug)
1326 cout << " in fun Chi2_ndf_pair[" << Pair1 << "][" << bst << "] " << Chi2_ndf_pair_[Pair1][bst] << endl;
1327
1328 ind_best_pair_[Pair1][bst] = bst;
1329 Nhits_pair_[Pair1][bst] = Nhits_Ch_[fir][best1] + Nhits_Ch_[sec][best2];
1330 for (Int_t i1 = 0; i1 < 4; i1++) {
1331 for (Int_t j1 = 0; j1 < 4; j1++) {
1332 sigma2_p[Pair1][bst][i1][j1] = 2 * bmatr[i1][j1];
1333 }
1334 }
1335
1336 } //< Nbest_Ch12_g_l
1337} // SegmentFit
1338//----------------------------------------------------------------------
1339
1340//--------------------Matrix Coefficients Calculation--------------
1341void BmnMwpcTrackFinder::FillFitMatrix(Int_t chN, Double_t** AA, Float_t** z, Double_t* sigmm2, Int_t* h_)
1342{
1343
1344 // AA - matrix to be filledlayers)
1345 // sigm2 - square of sigma
1346 // h_ - array to include/exclude planes (h_[i] = 0 or 1)
1347
1348 Float_t z2_[6] = {z[chN][0] * z[chN][0], z[chN][1] * z[chN][1], z[chN][2] * z[chN][2],
1349 z[chN][3] * z[chN][3], z[chN][4] * z[chN][4], z[chN][5] * z[chN][5]}; // cm
1350
1351 AA[0][0] += 2 * z2_[0] * h_[0] / sigmm2[0] + z2_[2] * h_[2] / (2 * sigmm2[2]) + z2_[1] * h_[1] / (2 * sigmm2[1])
1352 + 2 * z2_[3] * h_[3] / sigmm2[3] + z2_[5] * h_[5] / (2 * sigmm2[5])
1353 + z2_[4] * h_[4] / (2 * sigmm2[4]); // Ax
1354 AA[0][1] += 2 * z[chN][0] * h_[0] / sigmm2[0] + z[chN][2] * h_[2] / (2 * sigmm2[2])
1355 + z[chN][1] * h_[1] / (2 * sigmm2[1]) + 2 * z[chN][3] * h_[3] / sigmm2[3]
1356 + z[chN][5] * h_[5] / (2 * sigmm2[5]) + z[chN][4] * h_[4] / (2 * sigmm2[4]); // Bx
1357 AA[0][2] += sq3
1358 * (z2_[2] * h_[2] / (2 * sigmm2[2]) - z2_[1] * h_[1] / (2 * sigmm2[1])
1359 + z2_[5] * h_[5] / (2 * sigmm2[5]) - z2_[4] * h_[4] / (2 * sigmm2[4])); // Ay
1360 AA[0][3] += sq3
1361 * (z[chN][2] * h_[2] / (2 * sigmm2[2]) - z[chN][1] * h_[1] / (2 * sigmm2[1])
1362 + z[chN][5] * h_[5] / (2 * sigmm2[5]) - z[chN][4] * h_[4] / (2 * sigmm2[4])); // By
1363 AA[1][0] = AA[0][1];
1364 AA[1][1] += 2 * h_[0] / sigmm2[0] + 0.5 * h_[2] / sigmm2[2] + 0.5 * h_[1] / sigmm2[1] + 2 * h_[3] / sigmm2[3]
1365 + 0.5 * h_[5] / sigmm2[5] + 0.5 * h_[4] / sigmm2[4];
1366 AA[1][2] += sq3
1367 * (z[chN][2] * h_[2] / sigmm2[2] - z[chN][1] * h_[1] / sigmm2[1] + z[chN][5] * h_[5] / sigmm2[5]
1368 - z[chN][4] * h_[4] / sigmm2[4])
1369 * 0.5;
1370 AA[1][3] += sq3 * (h_[2] / sigmm2[2] - h_[1] / sigmm2[1] + h_[5] / sigmm2[5] - h_[4] / sigmm2[4]) * 0.5;
1371 AA[2][0] = AA[0][2];
1372 AA[2][1] = AA[1][2];
1373 AA[2][2] += 3.0
1374 * (z2_[2] * h_[2] / sigmm2[2] + z2_[1] * h_[1] / sigmm2[1] + z2_[5] * h_[5] / sigmm2[5]
1375 + z2_[4] * h_[4] / sigmm2[4])
1376 * 0.5;
1377 AA[2][3] += 3.0
1378 * (z[chN][2] * h_[2] / sigmm2[2] + z[chN][1] * h_[1] / sigmm2[1] + z[chN][5] * h_[5] / sigmm2[5]
1379 + z[chN][4] * h_[4] / sigmm2[4])
1380 * 0.5;
1381 AA[3][0] = AA[0][3];
1382 AA[3][1] = AA[1][3];
1383 AA[3][2] = AA[2][3];
1384 AA[3][3] +=
1385 3.0 * (0.5 * h_[2] / sigmm2[2] + 0.5 * h_[1] / sigmm2[1] + 0.5 * h_[5] / sigmm2[5] + 0.5 * h_[4] / sigmm2[4]);
1386}
1387//----------------------------------------------------------------------
1388
1389//--------------------Matrix Coefficients Calculation--------------
1390void BmnMwpcTrackFinder::FillFreeCoefVector(Int_t ichNum,
1391 Double_t* F,
1392 Double_t*** XVU_,
1393 Int_t ise,
1394 Float_t** z,
1395 Double_t* sigmm2,
1396 Int_t* h_)
1397{
1398 // F - vector to be filled
1399 // XVU_ - coordinates of segment in chamber (Is it correct definition?)
1400 // segIdx - index of current segment
1401 // z - local z-positions of planes(layers)
1402 // sigmm2 - square of sigma
1403 // h_ - array to include/exclude planes (h_[i] = 0 or 1)
1404
1405 F[0] += 2 * XVU_[ichNum][0][ise] * z[ichNum][0] * h_[0] / sigmm2[0]
1406 + XVU_[ichNum][1][ise] * z[ichNum][1] * h_[1] / sigmm2[1]
1407 + XVU_[ichNum][2][ise] * z[ichNum][2] * h_[2] / sigmm2[2]
1408 + 2 * XVU_[ichNum][3][ise] * z[ichNum][3] * h_[3] / sigmm2[3]
1409 + XVU_[ichNum][4][ise] * z[ichNum][4] * h_[4] / sigmm2[4]
1410 + XVU_[ichNum][5][ise] * z[ichNum][5] * h_[5] / sigmm2[5];
1411 F[1] += 2 * XVU_[ichNum][0][ise] * h_[0] / sigmm2[0] + XVU_[ichNum][1][ise] * h_[1] / sigmm2[1]
1412 + XVU_[ichNum][2][ise] * h_[2] / sigmm2[2] + 2 * XVU_[ichNum][3][ise] * h_[3] / sigmm2[3]
1413 + XVU_[ichNum][4][ise] * h_[4] / sigmm2[4] + XVU_[ichNum][5][ise] * h_[5] / sigmm2[5];
1414 F[2] += sq3
1415 * (-XVU_[ichNum][1][ise] * z[ichNum][1] * h_[1] / sigmm2[1]
1416 + XVU_[ichNum][2][ise] * z[ichNum][2] * h_[2] / sigmm2[2]
1417 - XVU_[ichNum][4][ise] * z[ichNum][4] * h_[4] / sigmm2[4]
1418 + XVU_[ichNum][5][ise] * z[ichNum][5] * h_[5] / sigmm2[5]);
1419 F[3] += sq3
1420 * (-XVU_[ichNum][1][ise] * h_[1] / sigmm2[1] + XVU_[ichNum][2][ise] * h_[2] / sigmm2[2]
1421 - XVU_[ichNum][4][ise] * h_[4] / sigmm2[4] + XVU_[ichNum][5][ise] * h_[5] / sigmm2[5]);
1422}
1423//----------------------------------------------------------------------
1424
1425//--------------------Matrix Coefficients Calculation--------------
1426void BmnMwpcTrackFinder::InverseMatrix(Double_t** AA, Double_t** bb)
1427{
1428 // Gaussian algorithm for 4x4 matrix inversion
1429 Double_t factor;
1430 Double_t temp[4];
1431
1432 // Set b to I
1433 for (Int_t i1 = 0; i1 < 4; i1++)
1434 for (Int_t j1 = 0; j1 < 4; j1++)
1435 if (i1 == j1)
1436 bb[i1][j1] = 1.0;
1437 else
1438 bb[i1][j1] = 0.0;
1439
1440 for (Int_t i1 = 0; i1 < 4; i1++) {
1441 for (Int_t j1 = i1 + 1; j1 < 4; j1++) {
1442 if (fabs(AA[i1][i1]) < fabs(AA[j1][i1])) {
1443 for (Int_t l1 = 0; l1 < 4; l1++)
1444 temp[l1] = AA[i1][l1];
1445 for (Int_t l1 = 0; l1 < 4; l1++)
1446 AA[i1][l1] = AA[j1][l1];
1447 for (Int_t l1 = 0; l1 < 4; l1++)
1448 AA[j1][l1] = temp[l1];
1449 for (Int_t l1 = 0; l1 < 4; l1++)
1450 temp[l1] = bb[i1][l1];
1451 for (Int_t l1 = 0; l1 < 4; l1++)
1452 bb[i1][l1] = bb[j1][l1];
1453 for (Int_t l1 = 0; l1 < 4; l1++)
1454 bb[j1][l1] = temp[l1];
1455 }
1456 }
1457 factor = AA[i1][i1];
1458 for (Int_t j1 = 4 - 1; j1 > -1; j1--) {
1459 bb[i1][j1] /= factor;
1460 AA[i1][j1] /= factor;
1461 }
1462 for (Int_t j1 = i1 + 1; j1 < 4; j1++) {
1463 factor = -AA[j1][i1];
1464 for (Int_t k1 = 0; k1 < 4; k1++) {
1465 AA[j1][k1] += AA[i1][k1] * factor;
1466 bb[j1][k1] += bb[i1][k1] * factor;
1467 }
1468 }
1469 } // i1
1470 for (Int_t i1 = 3; i1 > 0; i1--) {
1471 for (Int_t j1 = i1 - 1; j1 > -1; j1--) {
1472 factor = -AA[j1][i1];
1473 for (Int_t k1 = 0; k1 < 4; k1++) {
1474 AA[j1][k1] += AA[i1][k1] * factor;
1475 bb[j1][k1] += bb[i1][k1] * factor;
1476 // if ( fDebug) cout<<"AA["<<j1<<"]["<<k1<<"] "<<AA[j1][k1]<<endl;
1477 // if ( fDebug) cout<<"bb["<<j1<<"]["<<k1<<"] "<<bb[j1][k1]<<endl;
1478 }
1479 }
1480 } // i1
1481 // end inverse
1482}
1483//----------------------------------------------------------------------
1484
1485//--------------------Calculation Efficiency----------------------------
1486void BmnMwpcTrackFinder::FillEfficiency(Int_t ChN,
1487 Double_t*** XVU_Ch_,
1488 Int_t** Nhits_C,
1489 Int_t MinHits,
1490 Int_t ind_best,
1491 Float_t min_distX,
1492 Float_t min_distY)
1493{
1494 // cout<<" FillEfficiency : ChN "<<ChN<<" min_distX "<<min_distX<<" min_distY "<<min_distY<<"
1495 // Nhits_Ch["<<ChN<<"]["<<ind_best<<"] "<<Nhits_Ch[ChN][ ind_best ]<<endl;
1496 // segIdx - index of current segment ch2 or ch3 // Int_t ind_best_Ch[5]
1497 // 4p&4p -> all matched / Efficiency per layer
1498
1499 if (fabs(min_distX) < 5. && fabs(min_distY) < 5.5) {
1500 for (int i1 = 0; i1 < 6; i1++) {
1501 // if (fDebug) cout<<" XVU_Ch_["<<ChN<<"]["<<i1<<"]["<< ind_best<<"] = "<<XVU_Ch_[ChN][i1][ ind_best]<<endl;
1502 if (XVU_Ch_[ChN][i1][ind_best] > -999. && Nhits_Ch[ChN][ind_best] == MinHits)
1503 continue; // segIdx[ChN][j]
1504 if (fDebug)
1505 Denom_Ch.at(ChN)->Fill(i1);
1506 if (XVU_Ch_[ChN][i1][ind_best] > -999.) {
1507 if (fDebug)
1508 Nomin_Ch.at(ChN)->Fill(i1);
1509 }
1510 } // i1
1511 } // min_distX
1512} // FillEfficiency
1513//----------------------------------------------------------------------
1514
1515//----------------------------------------------------------------------
1516InitStatus BmnMwpcTrackFinder::Init()
1517{
1518
1519 FairRootManager* ioman = FairRootManager::Instance();
1520
1521 if (!expData) {
1522 if (fDebug)
1523 cout << " !expData " << endl;
1524 fBmnHitsArray = (TClonesArray*)ioman->GetObject(fInputBranchNameHit);
1525 if (fVerbose)
1526 cout << "fBmnHitsArray = " << fInputBranchNameHit << "\n";
1527 if (!fBmnHitsArray) {
1528 cout << "BmnMwpcTrackFinder::Init(): branchHits " << fInputBranchNameHit
1529 << " not found! Task will be deactivated" << endl;
1530 SetActive(kFALSE);
1531 return kERROR;
1532 }
1533 } else {
1534 if (fDebug)
1535 cout << " BmnMwpcTrackFinder::Init() " << endl;
1536 }
1537
1538 fBmnMwpcSegmentsArray = (TClonesArray*)ioman->GetObject(fInputBranchName);
1539 if (!fBmnMwpcSegmentsArray) {
1540 cout << "BmnMwpcTrackFinder::Init(): branch " << fInputBranchName << " not found! Task will be deactivated"
1541 << endl;
1542 SetActive(kFALSE);
1543 return kERROR;
1544 }
1545
1546 fBmnMwpcTracksArray = new TClonesArray(fOutputBranchName.Data());
1547 ioman->Register(fOutputBranchName.Data(), "MWPC", fBmnMwpcTracksArray, kTRUE);
1548
1549 fMwpcGeo = new BmnMwpcGeometrySRC(fRunPeriod, fRunNumber);
1550 kNChambers = fMwpcGeo->GetNChambers();
1551 kNPlanes = fMwpcGeo->GetNPlanes(); // 6
1552 if (fDebug)
1553 printf("fRunPeriod: %d %d %d\n", fRunPeriod, kNChambers, kNPlanes);
1554
1555 kMinHits = 4;
1556 kChi2_Max = 20.;
1557 kmaxPairs = 15; // 10;//5;
1558
1559 dw = fMwpcGeo->GetWireStep(); // 0.25; // [cm] // wires step
1560 dw_half = 0.5 * dw;
1561 sq3 = sqrt(3.);
1562 sq2 = sqrt(2.);
1563 sq12 = sqrt(12.);
1564 sigma = dw / sq12;
1565 kMiddlePl = 47.25;
1566 ktarget_region = 10.; // 04.05.21 VVP // 4.;
1567 ktarget_regionY = 10.; // 4.*sq2;
1568 kX_target = 0;
1569 kY_target = 0;
1570 Z0_SRC = 0.;
1571 if (fRunPeriod == 7) {
1572 kX_target = 0.5;
1573 kY_target = -4.5;
1574 Z0_SRC = -647.476;
1575 }
1576 if (fRunPeriod == 8) {
1577 kX_target = 0.;
1578 kY_target = 0.;
1579 Z0_SRC = -574.91;
1580 }
1581
1582 ChCent = new TVector3[kNChambers];
1583 ZCh = new Float_t[kNChambers];
1584 kZmid = new Float_t[kNChambers];
1585 shift = new Float_t*[kNChambers];
1586 shift_pair = new Float_t*[kNumPairs];
1587 kZ_midle_pair = new Float_t[kNumPairs];
1588 XVU_Ch = new Double_t**[kNChambers];
1589 Clust_Ch = new Double_t**[kNChambers];
1590 par_ab_Ch = new Double_t**[kNChambers];
1591 par_ab_pair = new Double_t**[kNumPairs];
1592 sigma2_pair = new Double_t***[kNumPairs];
1593 kPln = new Int_t*[kNChambers];
1594 kZ_loc = new Float_t*[kNChambers];
1595 z_gl = new Float_t*[kNChambers];
1596 Chi2_match_pair = new Double_t*[kNumPairs];
1597 Chi2_ndf_pair = new Double_t*[kNumPairs];
1598 ind_best_pair = new Int_t*[kNumPairs];
1599 Chi2_ndf_Ch = new Double_t*[kNChambers];
1600 Nhits_Ch = new Int_t*[kNChambers];
1601 Nhits_match = new Int_t*[kNChambers];
1602 Nhits_pair = new Int_t*[kNumPairs];
1603 sigm2 = new Float_t[kNPlanes];
1604 ipl = new Int_t[kNPlanes];
1605 Nbest_pair = new Int_t[kNumPairs];
1606 Nbest_Ch = new Int_t[kNChambers];
1607 ind_best_Ch = new Int_t*[kNChambers];
1608 sigma_delta = new Float_t*[kNumPairs];
1609
1610 Amatr = new Double_t*[4];
1611 bmatr = new Double_t*[4];
1612
1613 for (Int_t ii = 0; ii < 4; ii++) {
1614 Amatr[ii] = new Double_t[4];
1615 bmatr[ii] = new Double_t[4];
1616 }
1617
1618 for (Int_t i = 0; i < kNChambers; i++) {
1619 if (fDebug) {
1620 TH1D* h;
1621 h = new TH1D(Form("par_Ax_Ch%d", i), Form("par_Ax_Ch%d", i), 100, -.4, .4);
1622 fList.Add(h);
1623 hpar_Ax_Ch.push_back(h);
1624 h = new TH1D(Form("par_Bx_Ch%d", i), Form("par_Bx_Ch%d", i), 100, -10., 10.0);
1625 fList.Add(h);
1626 hpar_Bx_Ch.push_back(h);
1627 h = new TH1D(Form("par_Ay_Ch%d", i), Form("par_Ay_Ch%d", i), 100, -.4, .4);
1628 fList.Add(h);
1629 hpar_Ay_Ch.push_back(h);
1630 h = new TH1D(Form("par_By_Ch%d", i), Form("par_By_Ch%d", i), 100, -10., 10.0);
1631 fList.Add(h);
1632 hpar_By_Ch.push_back(h);
1633 h = new TH1D(Form("Chi2_ndf_Ch%d", i), Form("Chi2_ndf_Ch%d", i), 100, 0., 20.0);
1634 fList.Add(h);
1635 hChi2_ndf_Ch.push_back(h);
1636 h = new TH1D(Form("Nomin_Ch%d", i), Form("Nomin_Ch%d", i), 6, 0., 6.);
1637 fList.Add(h);
1638 Nomin_Ch.push_back(h);
1639 h = new TH1D(Form("Denom_Ch%d", i), Form("Denom_Ch%d", i), 6, 0., 6.);
1640 fList.Add(h);
1641 Denom_Ch.push_back(h);
1642 h = new TH1D(Form("Efficiency_Ch%d", i), Form("Efficiency_Ch%d", i), 6, 0., 6.);
1643 fList.Add(h);
1644 Eff_Ch.push_back(h);
1645 }
1646
1647 shift[i] = new Float_t[4];
1648 ChCent[i] = fMwpcGeo->GetChamberCenter(i);
1649 ZCh[i] = ChCent[i].Z();
1650 shift[i][0] = fMwpcGeo->GetTx(i);
1651 shift[i][2] = fMwpcGeo->GetTy(i);
1652 shift[i][1] = ChCent[i].X();
1653 shift[i][3] = ChCent[i].Y();
1654 kPln[i] = new Int_t[kNPlanes];
1655 kZ_loc[i] = new Float_t[kNPlanes];
1656 z_gl[i] = new Float_t[kNPlanes];
1657 XVU_Ch[i] = new Double_t*[kNPlanes];
1658 Clust_Ch[i] = new Double_t*[kNPlanes];
1659 par_ab_Ch[i] = new Double_t*[4];
1660 Nhits_Ch[i] = new Int_t[kBig];
1661 Chi2_ndf_Ch[i] = new Double_t[kBig];
1662 ind_best_Ch[i] = new Int_t[kmaxPairs];
1663 Nhits_match[i] = new Int_t[kmaxPairs];
1664 }
1665
1666 //----- hists booking -----
1667 for (Int_t i = 0; i < kNChambers; i++) {
1668 if (i == 0 || i == 2) {
1669 kZmid[i] = (ZCh[i] - ZCh[i + 1]) * 0.5;
1670 }
1671 if (i == 1 || i == 3) {
1672 kZmid[i] = (ZCh[i - 1] - ZCh[i]) * -0.5;
1673 }
1674 if (fDebug)
1675 printf("Chamber %d Z: %f Zmid: %f\n", i, ZCh[i], kZmid[i]);
1676 }
1677
1678 for (int i = 0; i < kNumPairs; ++i) {
1679 if (fDebug) {
1680 TH1D* h;
1681 h = new TH1D(Form("Chi2_match_pair%d", i), Form("Chi2_match_pair%d", i), 100, 0., 100.0);
1682 fList.Add(h);
1683 hChi2_match_pair.push_back(h);
1684 h = new TH1D(Form("par_Ax_pair%d", i), Form("slopeX pair%d; [degrees]; Events", i), 100, -2.3, 2.3);
1685 fList.Add(h);
1686 hpar_Ax_pair.push_back(h);
1687 h = new TH1D(Form("par_Bx_pair%d", i), Form("posX pair%d; [cm]; Events", i), 100, -10., 10.0);
1688 fList.Add(h);
1689 hpar_Bx_pair.push_back(h);
1690 h = new TH1D(Form("par_Ay_pair%d", i), Form("slopeY pair%d; [degrees]; Events", i), 100, -2.3, 2.3);
1691 fList.Add(h);
1692 hpar_Ay_pair.push_back(h);
1693 h = new TH1D(Form("par_By_pair%d", i), Form("posY pair%d; [cm]; Events", i), 100, -10., 10.0);
1694 fList.Add(h);
1695 hpar_By_pair.push_back(h);
1696
1697 h = new TH1D(Form("Chi2_ndf_pair%d", i), Form("Chi2_ndf_pair%d", i), 30, 0., 30.0);
1698 fList.Add(h);
1699 hChi2_ndf_pair.push_back(h);
1700 h = new TH1D(Form("Nbest_pair%d", i), Form("Nbest_pair%d; Ntracks; Events", i), 5, 0., 5.);
1701 fList.Add(h);
1702 hNbest_pair.push_back(h);
1703 h = new TH1D(Form("theta_pair%d", i), Form("theta_pair%d; [degrees]; Events", i), 160, 0., 8.);
1704 fList.Add(h);
1705 hpar_theta_pair.push_back(h);
1706 h = new TH1D(Form("phi_pair%d", i), Form("phi_pair%d; [degrees]; Events", i), 380, -190., 190.);
1707 fList.Add(h);
1708 hpar_phi_pair.push_back(h);
1709
1710 h = new TH1D(Form("dX_pair%d", i), Form("dX_pair%d; cm; Events", i), 100, -5., 5.);
1711 fList.Add(h);
1712 hdX_pair.push_back(h);
1713 h = new TH1D(Form("dY_pair%d", i), Form("dY_pair%d; cm; Events", i), 100, -5., 5.);
1714 fList.Add(h);
1715 hdY_pair.push_back(h);
1716 h = new TH1D(Form("dAx_pair%d", i), Form("dAx_pair%d; rad;Events", i), 100, -.5, .5);
1717 fList.Add(h);
1718 hdAx_pair.push_back(h);
1719 h = new TH1D(Form("dAy_pair%d", i), Form("dAy_pair%d; rad;Events", i), 100, -.5, .5);
1720 fList.Add(h);
1721 hdAy_pair.push_back(h);
1722
1723 h = new TH1D(Form("X_in_target_pair%d", i), Form(" posX_pair%d in target;[cm]; Events ", i), 100, -10.,
1724 10.);
1725 fList.Add(h);
1726 hX_in_target_pair.push_back(h);
1727 h = new TH1D(Form("Y_in_target_pair%d", i), Form(" posY_pair%d in target;[cm]; Events ", i), 100, -10.,
1728 10.);
1729 fList.Add(h);
1730 hY_in_target_pair.push_back(h);
1731
1732 TH2D* h1;
1733 h1 = new TH2D(Form("YvsX_pair%d", i), Form("YvsX_pair%d; X[cm]; Y[cm]", i), 100, -25, 25, 100, -21, 21);
1734 fList.Add(h1);
1735 hYvsX_pair.push_back(h1);
1736
1737 h1 = new TH2D(Form("Y_X_in_target_pair%d", i), Form("posY vs posX pair%d in target; X[cm]; Y[cm]", i), 150,
1738 -7.5, 7.5, 150, -13., 2.);
1739 fList.Add(h1);
1740 hY_X_in_target_pair.push_back(h1);
1741 h1 = new TH2D(Form("Ax_bx_in_target_pair%d", i),
1742 Form("slopeX vs posX in target pair%d; posX[cm]; slopeX", i), 100, -10., 10., 100, -2.3, 2.3);
1743 fList.Add(h1);
1744 hAx_bx_in_target_pair.push_back(h1);
1745 h1 = new TH2D(Form("Ay_by_in_target_pair%d", i),
1746 Form("slopeY vs posY in target pair%d; posY[cm]; slopeY", i), 100, -10., 10., 100, -2.3, 2.3);
1747 fList.Add(h1);
1748 hAy_by_in_target_pair.push_back(h1);
1749 } // if (fDebug)
1750 Chi2_match_pair[i] = new Double_t[kmaxPairs];
1751 Chi2_ndf_pair[i] = new Double_t[kmaxPairs];
1752 ind_best_pair[i] = new Int_t[kmaxPairs];
1753 Nhits_pair[i] = new Int_t[kmaxPairs];
1754 par_ab_pair[i] = new Double_t*[4];
1755 sigma2_pair[i] = new Double_t**[kBig];
1756 shift_pair[i] = new Float_t[4];
1757 sigma_delta[i] = new Float_t[4];
1758 }
1759
1760 for (Int_t i = 0; i < kNumPairs; i++) {
1761 for (Int_t j = 0; j < kBig; j++) {
1762 sigma2_pair[i][j] = new Double_t*[4];
1763 }
1764 }
1765 for (Int_t i = 0; i < kNumPairs; i++) {
1766 for (Int_t j = 0; j < kBig; j++) {
1767 for (Int_t k = 0; k < 4; k++) {
1768 sigma2_pair[i][j][k] = new Double_t[4];
1769 }
1770 }
1771 }
1772
1773 if (fDebug) {
1774 hAx_bx_in_target =
1775 new TH2D("Ax_bx_in_target", "slopeX vs posX in target; posX[cm]; slopeX", 100, -10., 10., 100, -2.3, 2.3);
1776 fList.Add(hAx_bx_in_target);
1777 hAy_by_in_target =
1778 new TH2D("Ay_by_in_target", "slopeY vs posY in target; posY[cm]; slopeY", 100, -10., 10., 100, -2.3, 2.3);
1779 fList.Add(hAy_by_in_target);
1780 hY_X_in_target =
1781 new TH2D("Y_X_in_target", "posY vs posX (pair0) in target; X[cm]; Y[cm]", 100, -10., 10., 100, -10., 10.);
1782 fList.Add(hY_X_in_target);
1783
1784 htheta_p1vsp0 = new TH2D("theta_p1vsp0", "theta pair1 vs pair0", 160, 0., 3., 160, 0., 3.);
1785 fList.Add(htheta_p1vsp0);
1786 hdX_pair01_vs_x1 = new TH2D("dX_pair01_vs_x1", "dX(pair0- pair1)_vs_Xpair1;Xpair1[cm];dX(pair0- pair1)[cm]",
1787 100, -10., 10., 100, -10., 10.);
1788 fList.Add(hdX_pair01_vs_x1);
1789 hdY_pair01_vs_y1 = new TH2D("dY_pair01_vs_y1", "dY(pair0- pair1)_vs_Ypair1;Ypair1[cm];dY(pair0- pair1)[cm]",
1790 100, -10., 10., 100, -10., 10.);
1791 fList.Add(hdY_pair01_vs_y1);
1792 hdX_pair01_inZpair1 = new TH1D("dX_pair01_inZpair1", "dX(pair0- pair1) inZpair1;X[cm]; ", 100, -10., 10.);
1793 fList.Add(hdX_pair01_inZpair1);
1794 hdY_pair01_inZpair1 = new TH1D("dY_pair01_inZpair1", "dY(pair0- pair1) inZpair1;Y[cm]; ", 100, -10., 10.);
1795 fList.Add(hdY_pair01_inZpair1);
1796
1797 hdX_target = new TH1D("dX_target", " PosX(pair0-pair1) in target;[cm]; Events ", 200, -10., 10.);
1798 fList.Add(hdX_target);
1799 hdY_target = new TH1D("dY_target", " PosY(pair0-pair1) in target;[cm]; Events ", 200, -10., 10.);
1800 fList.Add(hdY_target);
1801 hdAx_target = new TH1D("dAx_target", "SlopeX(pair0-pair1); rad; Events ", 100, -.05, .05);
1802 fList.Add(hdAx_target);
1803 hdAy_target = new TH1D("dAy_target", "SlopeY(pair0-pair1); rad; Events ", 100, -.05, .05);
1804 fList.Add(hdAy_target);
1805 hChi2_m_target = new TH1D("Chi2_m_target", "Chi2_m in target;; Events ", 100, 0., 100.);
1806 fList.Add(hChi2_m_target);
1807
1808 hChi2best_Chi2fake_before_target =
1809 new TH2D("Chi2best_Chi2fake_before_target", "Chi2best_Chi2fake_before_target; Chi2_best; Chi2_second", 20,
1810 0., 20., 20, 0., 20.);
1811 fList.Add(hChi2best_Chi2fake_before_target);
1812 hChi2best_Chi2fake_after_target =
1813 new TH2D("Chi2best_Chi2fake_after_target", "Chi2best_Chi2fake_after_target; Chi2_best; Chi2_second", 100,
1814 0., 100., 100, 0., 100.);
1815 fList.Add(hChi2best_Chi2fake_after_target);
1816
1817 hdX_pair_1 = new TH1D("dX_pair_1", "dX_pair_1", 200, -20., 20.);
1818 fList.Add(hdX_pair_1);
1819 hdY_pair_1 = new TH1D("dY_pair_1", "dY_pair_1", 200, -20., 20.);
1820 fList.Add(hdY_pair_1);
1821 hdAx_pair_1 = new TH1D("dAx_pair_1", "dAx_pair_1", 100, -.3, .3);
1822 fList.Add(hdAx_pair_1);
1823 hdAy_pair_1 = new TH1D("dAy_pair_1", "dAy_pair_1", 100, -.3, .3);
1824 fList.Add(hdAy_pair_1);
1825 hChi2m_pair_1 = new TH1D("Chi2m_pair_1", "Chi2_xy_pair_1", 100, 0, 50);
1826 fList.Add(hChi2m_pair_1);
1827 hXv_pair_1 = new TH1D("Xv_pair_1", "Xv_pair_1", 200, -10., 10.);
1828 fList.Add(hXv_pair_1);
1829 hYv_pair_1 = new TH1D("Yv_pair_1", "Yv_pair_1", 200, -10., 10.);
1830 fList.Add(hYv_pair_1);
1831
1832 hDen = new TH1D("hDen", "hDen", 1, 0, 1);
1833 fList.Add(hDen);
1834 hNum = new TH1D("hNum", "hNum", 1, 0, 1);
1835 fList.Add(hNum);
1836 hEff = new TH1D("hEff", "hEff", 1, 0, 1);
1837 fList.Add(hEff);
1838
1839 hYvsX_mc_pair = new TH2D("hYvsX_mc_pair", "YvsX_mc_pair", 100, -25, 25, 100, -21, 21);
1840 fList.Add(hYvsX_mc_pair);
1841
1842 hdAx_tr_mc = new TH1D("hdAx_tr_mc", "dAx_tr_mc", 200, -0.005, 0.005);
1843 fList.Add(hdAx_tr_mc);
1844 hdX_tr_mc = new TH1D("hdX_tr_mc", "dX_tr_mc", 200, -0.2, 0.2);
1845 fList.Add(hdX_tr_mc);
1846 hdAy_tr_mc = new TH1D("hdAy_tr_mc", "dAy_tr_mc", 200, -0.005, 0.005);
1847 fList.Add(hdAy_tr_mc);
1848 hdY_tr_mc = new TH1D("hdY_tr_mc", "dY_tr_mc", 200, -0.2, 0.2);
1849 fList.Add(hdY_tr_mc);
1850
1851 hDen_mctr = new TH1D("hDen_mctr", "Den_mctr", 2, 0, 2);
1852 fList.Add(hDen_mctr);
1853 hNum_mctr = new TH1D("hNum_mctr", "Num_mctr", 2, 0, 2);
1854 fList.Add(hNum_mctr);
1855 hEff_mctr = new TH1D("hEff_mctr", "Eff_mctr", 2, 0, 2);
1856 fList.Add(hEff_mctr);
1857 hEff_mctr->GetXaxis()->SetBinLabel(1, "2+3 coord");
1858 hEff_mctr->GetXaxis()->SetBinLabel(2, "only 3coord");
1859
1860 hNtrpc_reco = new TH1D("Ntrpc_reco", ";N of pair1 tr reco", 10, 0, 10);
1861 fList.Add(hNtrpc_reco);
1862 hNtrpc_mc = new TH1D("Ntrpc_mc", ";N of tr pair1 mc", 10, 0, 10);
1863 fList.Add(hNtrpc_mc);
1864 hNtrpc_mc_vs_reco = new TH2D("Ntrpc_mc_vs_reco", "Ntr_mc_vs_reco; Nreco;Nmctrue", 10, 0, 10, 10, 0, 10);
1865 fList.Add(hNtrpc_mc_vs_reco);
1866
1867 } // if (fDebug)
1868 Int_t i1 = 0;
1869 for (Int_t i = 0; i < kNumPairs; i++) {
1870
1871 i1 = i;
1872 sigma_delta[0][0] = 3 * .14; //.0624;// sigm_dax
1873 sigma_delta[0][2] = 3 * .14; //.066; // sigm_day
1874 sigma_delta[0][1] = 2 * .35; // 4.08;// sigm_dx
1875 sigma_delta[0][3] = 2 * .35; // 4.30;// sigm_dy
1876
1877 if (i == 1) {
1878 i1 = 2;
1879 // pair 0
1880 sigma_delta[0][0] = 3 * .14;
1881 sigma_delta[0][1] = 1.2; // 2*.35;
1882 sigma_delta[0][2] = 3 * .14;
1883 sigma_delta[0][3] = 1.3; // 2*.35;
1884 // pair 1
1885 sigma_delta[1][0] = 6 * .1; // sigm_dax //VP
1886 sigma_delta[1][1] = 1.; // 1.;//6.;//4.08;//2*.35;//4.08;// sigm_dx
1887 sigma_delta[1][2] = 6 * .15; // sigm_day
1888 sigma_delta[1][3] = 1.5; // 1.//2*7.;//shift +10cm? 4.30;// 2*.35;//4.30;// sigm_dy
1889 }
1890 kZ_midle_pair[i] = ZCh[i1] + kZmid[i1 + 1];
1891
1892 shift_pair[i][0] = (shift[i1 + 1][1] - shift[i1][1]) / (ZCh[i1 + 1] - ZCh[i1]);
1893 shift_pair[i][2] = (shift[i1 + 1][3] - shift[i1][3]) / (ZCh[i1 + 1] - ZCh[i1]);
1894 shift_pair[i][1] = 0.5 * (shift[i1 + 1][1] + shift[i1][1]);
1895 shift_pair[i][3] = 0.5 * (shift[i1 + 1][3] + shift[i1][3]);
1896 // Additional alignment
1897 /*
1898 if (i == 0){// alignment with GEM detectors
1899 shift_pair[i][0] += 0.;//Ax
1900 shift_pair[i][1] += 0.05;//Bx
1901 shift_pair[i][2] += -0.0008;//Ay
1902 shift_pair[i][3] += -0.4;//By
1903 }
1904 */
1905 if (i == 0) {
1906 shift_pair[i][0] += 0.; // Ax
1907 shift_pair[i][1] += 0.; // Bx
1908 shift_pair[i][2] += 0.; // Ay
1909 shift_pair[i][3] += 0.; // By
1910 }
1911 if (i == 1 && expData) {
1912 shift_pair[i][0] += 0.; // Ax
1913 shift_pair[i][1] += -0.05; // Bx
1914 shift_pair[i][2] += 0.0009; // Ay
1915 shift_pair[i][3] += 0.7; // By
1916 }
1917
1918 // if (fDebug) cout<<" i "<<i<<" kZ_midle_pair[i] "<<kZ_midle_pair[i]<<" i1 "<<i1<<" i1+1 "<<i1+1<<" -( ZCh[i1]-
1919 // ZCh[i1+1] )= "<<-( ZCh[i1]- ZCh[i1+1] )<<endl;
1920 }
1921
1922 for (Int_t ichh = 0; ichh < kNChambers; ichh++) {
1923 for (int ii = 0; ii < 6; ii++) {
1924
1925 if (fRunPeriod == 6) {
1926
1927 if (ichh == 0 || ichh == 1) {
1928 kZ_loc[ichh][ii] = -0.5 + ii;
1929 if (ii == 4) {
1930 kZ_loc[ichh][ii] = -2.5;
1931 }
1932 if (ii == 5) {
1933 kZ_loc[ichh][ii] = -1.5;
1934 }
1935 }
1936 }
1937
1938 if (fRunPeriod == 6 || (fRunPeriod == 7 && fRunNumber > 3588)) // BMN
1939 {
1940 if (ichh == 0 || ichh == 1) {
1941 kZ_loc[ichh][ii] = -0.5 + ii;
1942 if (ii == 4) {
1943 kZ_loc[ichh][ii] = -2.5;
1944 }
1945 if (ii == 5) {
1946 kZ_loc[ichh][ii] = -1.5;
1947 }
1948 }
1949 }
1950 if ((fRunPeriod == 7 && fRunNumber <= 3588) || (fRunPeriod == 8)) // SRC
1951 {
1952 if (ichh == 0) {
1953 kZ_loc[ichh][0] = -1.5;
1954 kZ_loc[ichh][1] = -0.5;
1955 kZ_loc[ichh][2] = 0.5;
1956 kZ_loc[ichh][3] = 1.5;
1957 kZ_loc[ichh][4] = 2.5;
1958 kZ_loc[ichh][5] = -2.5;
1959 }
1960 if (ichh == 1) {
1961 kZ_loc[ichh][0] = -1.5;
1962 kZ_loc[ichh][1] = -2.5;
1963 kZ_loc[ichh][2] = 2.5;
1964 kZ_loc[ichh][3] = 1.5;
1965 kZ_loc[ichh][4] = 0.5;
1966 kZ_loc[ichh][5] = -0.5;
1967 }
1968 if (fRunPeriod == 7 && (ichh == 2 || ichh == 3)) {
1969
1970 kZ_loc[ichh][ii] = -0.5 + ii;
1971 if (ii == 4) {
1972 kZ_loc[ichh][ii] = -2.5;
1973 }
1974 if (ii == 5) {
1975 kZ_loc[ichh][ii] = -1.5;
1976 }
1977 } // fRunPeriod == 7 ch 2 3
1978 if (fRunPeriod == 8 && (ichh == 2 || ichh == 3)) {
1979 kZ_loc[ichh][0] = 0.5;
1980 kZ_loc[ichh][1] = -0.5;
1981 kZ_loc[ichh][2] = -1.5;
1982 kZ_loc[ichh][3] = -2.5;
1983 kZ_loc[ichh][4] = 2.5;
1984 kZ_loc[ichh][5] = 1.5;
1985 } // fRunPeriod == 8 ch 2 3
1986
1987 } // SRC
1988
1989 z_gl[ichh][ii] = kZmid[ichh] + kZ_loc[ichh][ii];
1990
1991 if (fDebug)
1992 cout << " ich " << ichh << " ii " << ii << " kZ_loc " << kZ_loc[ichh][ii] << " z_gl " << z_gl[ichh][ii]
1993 << endl;
1994 }
1995 } // ich
1996
1997 for (Int_t ii = 0; ii < kNChambers; ii++) {
1998 for (Int_t iPla = 0; iPla < kNPlanes; iPla++) {
1999 XVU_Ch[ii][iPla] = new Double_t[kBig];
2000 Clust_Ch[ii][iPla] = new Double_t[kBig];
2001 }
2002 for (Int_t i = 0; i < 4; i++) {
2003 par_ab_Ch[ii][i] = new Double_t[kBig];
2004 }
2005 }
2006
2007 for (Int_t ip = 0; ip < kNumPairs; ip++) {
2008 for (Int_t i4 = 0; i4 < 4; i4++) {
2009 par_ab_pair[ip][i4] = new Double_t[kmaxPairs];
2010 }
2011 }
2012
2013 return kSUCCESS;
2014} // Init
2015//----------------------------------------------------------------------
2016
2017//------ Arrays Initialization -----------------------------------------
2018void BmnMwpcTrackFinder::PrepareArraysToProcessEvent()
2019{
2020 fBmnMwpcTracksArray->Clear();
2021 vec_points.clear();
2022
2023 // Clean and initialize arrays:
2024
2025 for (Int_t iCh = 0; iCh < kNChambers; iCh++) {
2026 Nbest_Ch[iCh] = 0;
2027
2028 for (Int_t iPlane = 0; iPlane < kNPlanes; iPlane++) {
2029 for (Int_t iBig = 0; iBig < kBig; iBig++) {
2030 XVU_Ch[iCh][iPlane][iBig] = -999.;
2031 Clust_Ch[iCh][iPlane][iBig] = -1.;
2032 } // iBig
2033 } // iPlane
2034
2035 for (Int_t ii = 0; ii < 4; ii++) {
2036 for (Int_t jj = 0; jj < kBig; jj++) {
2037 par_ab_Ch[iCh][ii][jj] = 999.;
2038 }
2039 }
2040
2041 for (Int_t iBig = 0; iBig < kBig; iBig++) {
2042 Nhits_Ch[iCh][iBig] = 0;
2043 Chi2_ndf_Ch[iCh][iBig] = 0;
2044 }
2045
2046 for (Int_t i = 0; i < kmaxPairs; i++) {
2047 ind_best_Ch[iCh][i] = -1;
2048 Nhits_match[iCh][i] = 0;
2049 }
2050 } // iCh
2051
2052 for (Int_t iPl = 0; iPl < kNPlanes; iPl++) {
2053 sigm2[iPl] = sigma * sigma;
2054 ipl[iPl] = 6;
2055 }
2056
2057 for (Int_t ip = 0; ip < kNumPairs; ip++) {
2058 Nbest_pair[ip] = 0; // Nbest_Ch12_gl = 0;
2059
2060 for (Int_t i4 = 0; i4 < 4; i4++) {
2061 for (Int_t i5 = 0; i5 < kmaxPairs; i5++) {
2062 par_ab_pair[ip][i4][i5] = 999.; // par_ab_Ch1_2[ii][jj] = 999.;
2063 }
2064 }
2065 for (Int_t i = 0; i < kmaxPairs; i++) {
2066 Chi2_match_pair[ip][i] = 999.;
2067 Chi2_ndf_pair[ip][i] = 999.;
2068 ind_best_pair[ip][i] = -1;
2069 }
2070 for (Int_t j = 0; j < kBig; j++) {
2071 for (Int_t k = 0; k < 4; k++) {
2072 for (Int_t i4 = 0; i4 < 4; i4++) {
2073 sigma2_pair[ip][j][k][i4] = 0.;
2074 }
2075 }
2076 }
2077 }
2078
2079} // PrepareArraysToProcessEvent
2080//----------------------------------------------------------------------
2081
2082//----------------------------------------------------------------------
2083void BmnMwpcTrackFinder::Finish()
2084{
2085 delete fMwpcGeo;
2086
2087 if (fDebug) {
2088 printf("MWPC track finder: write hists to file... ");
2089 fOutputFileName = Form("hMWPCtracks_p%d_run%d.root", fRunPeriod, fRunNumber);
2090 cout << fOutputFileName << endl;
2091 TFile file(fOutputFileName, "RECREATE");
2092
2093 for (Int_t iCh = 0; iCh < kNChambers; iCh++) {
2094 Eff_Ch.at(iCh)->Sumw2();
2095 Eff_Ch.at(iCh)->Divide(Nomin_Ch.at(iCh), Denom_Ch.at(iCh), 1, 1);
2096 }
2097
2098 hEff->Divide(hNum, hDen, 1, 1);
2099 hEff_mctr->Divide(hNum_mctr, hDen_mctr, 1, 1);
2100 fList.Write();
2101 file.Close();
2102 }
2103 if (fDebug)
2104 printf("done\n");
2105
2106 cout << "Work time of the MWPC track finder: " << workTime << " s" << endl;
2107} // Finish
2108//----------------------------------------------------------------------
compareSegments
bool compareSegments(const segments &a, const segments &b)
Definition BmnMwpcHitFinder.cxx:44
compareSegments
bool compareSegments(const match &a, const match &b)
Definition BmnMwpcTrackFinder.cxx:41
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
i
int i
Definition P4_F32vec4.h:22
Bool_t
bool Bool_t
Int_t
int Int_t
Float_t
float Float_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,...)
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
BmnMwpcTrackFinder::FillFreeCoefVector
void FillFreeCoefVector(Int_t, Double_t *, Double_t ***, Int_t, Float_t **, Double_t *, Int_t *)
Definition BmnMwpcTrackFinder.cxx:1390
BmnMwpcTrackFinder::Exec
virtual void Exec(Option_t *opt)
Definition BmnMwpcTrackFinder.cxx:162
BmnMwpcTrackFinder::hChi2_ndf_Ch
vector< TH1D * > hChi2_ndf_Ch
Definition BmnMwpcTrackFinder.h:92
BmnMwpcTrackFinder::hpar_Bx_pair
vector< TH1D * > hpar_Bx_pair
Definition BmnMwpcTrackFinder.h:93
BmnMwpcTrackFinder::fBmnHitsArray
TClonesArray * fBmnHitsArray
Definition BmnMwpcTrackFinder.h:85
BmnMwpcTrackFinder::Eff_Ch
vector< TH1D * > Eff_Ch
Definition BmnMwpcTrackFinder.h:93
BmnMwpcTrackFinder::Nomin_Ch
vector< TH1D * > Nomin_Ch
Definition BmnMwpcTrackFinder.h:93
BmnMwpcTrackFinder::hpar_Ax_pair
vector< TH1D * > hpar_Ax_pair
Definition BmnMwpcTrackFinder.h:92
BmnMwpcTrackFinder::hpar_Ax_Ch
vector< TH1D * > hpar_Ax_Ch
Definition BmnMwpcTrackFinder.h:92
BmnMwpcTrackFinder::SegmentMatching
void SegmentMatching(Int_t, Int_t *, Double_t ***, Float_t *, Int_t **, Int_t *, Double_t **, Double_t ***, Int_t **, Int_t **)
Definition BmnMwpcTrackFinder.cxx:687
BmnMwpcTrackFinder::fBmnMwpcTracksArray
TClonesArray * fBmnMwpcTracksArray
Output array of MWPC tracks.
Definition BmnMwpcTrackFinder.h:87
BmnMwpcTrackFinder::hY_in_target_pair
vector< TH1D * > hY_in_target_pair
Definition BmnMwpcTrackFinder.h:94
BmnMwpcTrackFinder::hpar_Ay_pair
vector< TH1D * > hpar_Ay_pair
Definition BmnMwpcTrackFinder.h:93
BmnMwpcTrackFinder::hdX_pair
vector< TH1D * > hdX_pair
Definition BmnMwpcTrackFinder.h:94
BmnMwpcTrackFinder::hAx_bx_in_target_pair
vector< TH2D * > hAx_bx_in_target_pair
Definition BmnMwpcTrackFinder.h:97
BmnMwpcTrackFinder::hChi2_ndf_pair
vector< TH1D * > hChi2_ndf_pair
Definition BmnMwpcTrackFinder.h:95
BmnMwpcTrackFinder::hY_X_in_target_pair
vector< TH2D * > hY_X_in_target_pair
Definition BmnMwpcTrackFinder.h:97
BmnMwpcTrackFinder::hdAy_pair
vector< TH1D * > hdAy_pair
Definition BmnMwpcTrackFinder.h:94
BmnMwpcTrackFinder::Denom_Ch
vector< TH1D * > Denom_Ch
Definition BmnMwpcTrackFinder.h:93
BmnMwpcTrackFinder::hdY_pair
vector< TH1D * > hdY_pair
Definition BmnMwpcTrackFinder.h:94
BmnMwpcTrackFinder::hYvsX_pair
vector< TH2D * > hYvsX_pair
Definition BmnMwpcTrackFinder.h:97
BmnMwpcTrackFinder::hAy_by_in_target_pair
vector< TH2D * > hAy_by_in_target_pair
Definition BmnMwpcTrackFinder.h:97
BmnMwpcTrackFinder::fBmnMwpcSegmentsArray
TClonesArray * fBmnMwpcSegmentsArray
Input array of MWPC hits.
Definition BmnMwpcTrackFinder.h:84
BmnMwpcTrackFinder::hpar_phi_pair
vector< TH1D * > hpar_phi_pair
Definition BmnMwpcTrackFinder.h:93
BmnMwpcTrackFinder::PairMatching
void PairMatching(Int_t *, Double_t ***, Float_t *)
Definition BmnMwpcTrackFinder.cxx:1057
BmnMwpcTrackFinder::TracksStoring
void TracksStoring(Int_t *, Double_t ***, Double_t **, Int_t **, Double_t ****)
Definition BmnMwpcTrackFinder.cxx:247
BmnMwpcTrackFinder::SegmentMatchingAfterTarget
void SegmentMatchingAfterTarget(Int_t, Int_t *, Double_t ***, Float_t *, Int_t **, Int_t *, Double_t **, Double_t ***, Int_t **, Int_t **)
Definition BmnMwpcTrackFinder.cxx:845
BmnMwpcTrackFinder::vec_points
vector< MC_points > vec_points
Definition BmnMwpcTrackFinder.h:113
BmnMwpcTrackFinder::hpar_By_Ch
vector< TH1D * > hpar_By_Ch
Definition BmnMwpcTrackFinder.h:92
BmnMwpcTrackFinder::Init
virtual InitStatus Init()
Definition BmnMwpcTrackFinder.cxx:1516
BmnMwpcTrackFinder::fMwpcGeo
BmnMwpcGeometrySRC * fMwpcGeo
Definition BmnMwpcTrackFinder.h:89
BmnMwpcTrackFinder::ReadSegments
void ReadSegments(Double_t ***, Int_t **, Double_t **, Int_t *, Double_t ***, Double_t ***, vector< MC_points > &)
Definition BmnMwpcTrackFinder.cxx:431
BmnMwpcTrackFinder::hX_in_target_pair
vector< TH1D * > hX_in_target_pair
Definition BmnMwpcTrackFinder.h:94
BmnMwpcTrackFinder::hpar_By_pair
vector< TH1D * > hpar_By_pair
Definition BmnMwpcTrackFinder.h:93
BmnMwpcTrackFinder::hNbest_pair
vector< TH1D * > hNbest_pair
Definition BmnMwpcTrackFinder.h:95
BmnMwpcTrackFinder::hpar_Ay_Ch
vector< TH1D * > hpar_Ay_Ch
Definition BmnMwpcTrackFinder.h:92
BmnMwpcTrackFinder::InverseMatrix
void InverseMatrix(Double_t **, Double_t **)
Definition BmnMwpcTrackFinder.cxx:1426
BmnMwpcTrackFinder::SegmentParamAlignment
void SegmentParamAlignment(Int_t, Int_t *, Double_t ***, Float_t **)
Definition BmnMwpcTrackFinder.cxx:667
BmnMwpcTrackFinder::hpar_theta_pair
vector< TH1D * > hpar_theta_pair
Definition BmnMwpcTrackFinder.h:93
BmnMwpcTrackFinder::hdAx_pair
vector< TH1D * > hdAx_pair
Definition BmnMwpcTrackFinder.h:94
BmnMwpcTrackFinder::SegmentFit
void SegmentFit(Int_t, Float_t **, Float_t *, Int_t *, Int_t **, Double_t ***, Double_t **, Double_t ***, Double_t ***, Int_t **, Int_t **, Int_t **, Double_t ****)
Definition BmnMwpcTrackFinder.cxx:1140
BmnMwpcTrackFinder::hChi2_match_pair
vector< TH1D * > hChi2_match_pair
Definition BmnMwpcTrackFinder.h:92
BmnMwpcTrackFinder::FillFitMatrix
void FillFitMatrix(Int_t, Double_t **, Float_t **, Double_t *, Int_t *)
Definition BmnMwpcTrackFinder.cxx:1341
BmnMwpcTrackFinder::MCefficiencyCalculation
void MCefficiencyCalculation(vector< MC_points > &, Int_t *, Double_t ***)
Definition BmnMwpcTrackFinder.cxx:312
BmnMwpcTrackFinder::hpar_Bx_Ch
vector< TH1D * > hpar_Bx_Ch
Definition BmnMwpcTrackFinder.h:92
BmnMwpcTrackFinder::FillEfficiency
void FillEfficiency(Int_t, Double_t ***, Int_t **, Int_t, Int_t, Float_t, Float_t)
Definition BmnMwpcTrackFinder.cxx:1486
BmnMwpcTrack::GetCoord
vector< Double_t > GetCoord()
Definition BmnMwpcTrack.h:17
BmnMwpcTrack::GetClust
vector< Double_t > GetClust()
Definition BmnMwpcTrack.h:19
BmnTrack::SetChi2
void SetChi2(Double_t chi2)
Definition BmnTrack.h:76
BmnTrack::SetFlag
void SetFlag(Int_t flag)
Definition BmnTrack.h:74
BmnTrack::GetParamFirst
FairTrackParam * GetParamFirst()
Definition BmnTrack.h:63
BmnTrack::GetChi2
Float_t GetChi2() const
Definition BmnTrack.h:55
BmnTrack::GetNHits
Int_t GetNHits() const
Accessors.
Definition BmnTrack.h:49
BmnTrack::SetNHits
void SetNHits(Int_t n)
Definition BmnTrack.h:80
BmnTrack::SetParamFirst
void SetParamFirst(FairTrackParam &par)
Modifiers.
Definition BmnTrack.h:70
BmnTrack::GetFlag
Int_t GetFlag() const
Definition BmnTrack.h:53
TFile::Close
void Close(Option_t *option="") override
sqrt
RVec< PromoteType< T > > sqrt(const RVec< T > &v)
TMath::ATan2
Double_t ATan2(Double_t y, Double_t x)
TMath::Cos
Double_t Cos(Double_t)
TMath::RadToDeg
constexpr Double_t RadToDeg()
std
Definition json.hpp:4446
BmnMwpcTrackFinder::MC_points
Definition BmnMwpcTrackFinder.h:47
BmnMwpcTrackFinder::MC_points::Id
Int_t Id
Definition BmnMwpcTrackFinder.h:68
BmnMwpcTrackFinder::MC_points::z2
Double_t z2[6]
Definition BmnMwpcTrackFinder.h:50
BmnMwpcTrackFinder::MC_points::Np2
Int_t Np2
Definition BmnMwpcTrackFinder.h:65
BmnMwpcTrackFinder::MC_points::y3
Double_t y3[6]
Definition BmnMwpcTrackFinder.h:52
BmnMwpcTrackFinder::MC_points::y2
Double_t y2[6]
Definition BmnMwpcTrackFinder.h:49
BmnMwpcTrackFinder::MC_points::z3
Double_t z3[6]
Definition BmnMwpcTrackFinder.h:53
BmnMwpcTrackFinder::MC_points::x2
Double_t x2[6]
Definition BmnMwpcTrackFinder.h:48
BmnMwpcTrackFinder::MC_points::x3
Double_t x3[6]
Definition BmnMwpcTrackFinder.h:51
BmnMwpcTrackFinder::MC_points::Np3
Int_t Np3
Definition BmnMwpcTrackFinder.h:66
BmnMwpcTrackFinder::MC_points::Np
Int_t Np
Definition BmnMwpcTrackFinder.h:67
match::Chi2m
Double_t Chi2m
Definition BmnMwpcTrackFinder.cxx:28
match::param1
Double_t param1[4]
Definition BmnMwpcTrackFinder.cxx:33
match::param2
Double_t param2[4]
Definition BmnMwpcTrackFinder.cxx:34