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