BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
BmnMwpcHitFinder.cxx
Go to the documentation of this file.
1// Author: Vasilisa Lenivenko <vasilisa@jinr.ru> 2018-07-18
2
3////////////////////////////////////////////////////////////////////////////////
4// //
5// BmnMwpcHitFinder //
6// //
7// Implementation of an algorithm developed by //
8// Vasilisa Lenivenko and Vladimir Palchik //
9// to the BmnRoot software //
10// //
11// The algorithm serves for searching for hits //
12// in the Multi Wire Prop. Chambers of the BM@N experiment //
13// //
14////////////////////////////////////////////////////////////////////////////////
15#include "BmnMwpcHitFinder.h"
16
17#include "BmnMwpcSegment.h"
18#include "BmnMwpcTrack.h"
19#include "BmnTrack.h"
20#include "FairRootManager.h"
21#include "TCanvas.h"
22#include "TFile.h"
23
24#include <BmnEventHeader.h>
25#include <algorithm>
26#include <climits>
27#include <vector>
28
29static Float_t workTime = 0.0;
30
31using namespace std;
32using namespace TMath;
33
34struct segments
35{
36 Int_t Nhits = 0;
37 Double_t Chi2 = 999.;
38 Double_t coord[6] = {-999., -999., -999., -999., -999., -999.};
39 Double_t clust[6] = {0., 0., 0., 0., 0., 0.};
40 Double_t param[4] = {999., 999., 999., 999.};
41 Double_t sigma2[4][4] = {{0., 0., 0., 0.}, {0., 0., 0., 0.}, {0., 0., 0., 0.}, {0., 0., 0., 0.}};
42};
43
44bool compareSegments(const segments& a, const segments& b)
45{
46 return a.Chi2 < b.Chi2;
47}
48
49BmnMwpcHitFinder::BmnMwpcHitFinder(Bool_t isExp, Int_t runPeriod, Int_t runNumber)
50 : expData(isExp)
51 , fEventNo(0)
52{
53 if (expData) {
54 fRunPeriod = runPeriod;
55 fRunNumber = runNumber;
56 fInputBranchName = "MWPC";
57 fBmnEventHeaderBranchName = "EventHeader";
58
59 if (fRunPeriod == 6 || (fRunPeriod == 7 && fRunNumber > 3588)) // bmn
60 {
61 kNumPairs = 1;
62 kCh_max = 2;
63 } else if ((fRunPeriod == 7 && fRunNumber <= 3588) || (fRunPeriod == 8)) // SRC
64 {
65 kNumPairs = 2;
66 kCh_max = 4;
67 }
68 } else {
69 fInputBranchName = "BmnMwpcHit";
70 kNumPairs = 2;
71 kCh_max = 4;
72 fRunPeriod = 7;
73 fRunNumber = 0001;
74 }
75
76 fOutputBranchName = "BmnMwpcSegment";
77 nInputDigits = 3;
78 nTimeSamples = 3;
79 kBig = 100;
80 kCh_min = 0;
81 // fBmnEvQualityBranchName = "BmnEventQuality";
82}
83
84BmnMwpcHitFinder::~BmnMwpcHitFinder()
85{
86
87 fBmnMwpcSegmentsArray->Delete();
88 delete fBmnMwpcSegmentsArray;
89
90 for (Int_t iCh = 0; iCh < kNChambers; ++iCh) {
91 for (Int_t iPlane = 0; iPlane < kNPlanes; ++iPlane) {
92 delete[] Wires_Ch[iCh][iPlane];
93 delete[] clust_Ch[iCh][iPlane];
94 delete[] XVU_Ch[iCh][iPlane];
95 delete[] DigitsArray[iCh][iPlane];
96 delete[] ClusterSize[iCh][iPlane];
97 delete[] Coord_wire[iCh][iPlane];
98 delete[] Coord_xuv[iCh][iPlane];
99 delete[] XVU_coord[iCh][iPlane];
100 delete[] Cluster_coord[iCh][iPlane];
101 delete[] Coor_seg[iCh][iPlane];
102 delete[] Cluster_seg[iCh][iPlane];
103 }
104
105 for (Int_t iBig = 0; iBig < kBig; ++iBig) {
106 for (Int_t i = 0; i < 4; ++i) {
107 delete[] sigma2_seg[iCh][iBig][i];
108 }
109 delete[] sigma2_seg[iCh][iBig];
110 }
111
112 for (Int_t iWire = 0; iWire < kNWires; ++iWire) {
113 delete[] wire_Ch[iCh][iWire];
114 delete[] xuv_Ch[iCh][iWire];
115 }
116
117 for (Int_t i = 0; i < 4; ++i) {
118 delete[] par_ab_Ch[iCh][i];
119 delete[] par_ab_seg[iCh][i];
120 }
121
122 delete[] sigma2_seg[iCh];
123 delete[] Wires_Ch[iCh];
124 delete[] clust_Ch[iCh];
125 delete[] XVU_Ch[iCh];
126 delete[] DigitsArray[iCh];
127 delete[] ClusterSize[iCh];
128 delete[] Coord_wire[iCh];
129 delete[] Coord_xuv[iCh];
130 delete[] XVU_coord[iCh];
131 delete[] Cluster_coord[iCh];
132 delete[] Coor_seg[iCh];
133 delete[] Cluster_seg[iCh];
134 delete[] par_ab_Ch[iCh];
135 delete[] par_ab_seg[iCh];
136 delete[] kPln[iCh];
137 delete[] iw_Ch[iCh];
138 delete[] shift[iCh];
139 delete[] Nhits_Ch[iCh];
140 delete[] Beam_wires_min[iCh];
141 delete[] Beam_wires_max[iCh];
142 delete[] ind_best_Ch[iCh];
143 delete[] best_Ch_gl[iCh];
144 delete[] Chi2_ndf_Ch[iCh];
145 delete[] Chi2_ndf_best_Ch[iCh];
146 delete[] XVU[iCh];
147 delete[] XVU_cl[iCh];
148 delete[] kZ_loc[iCh];
149 delete[] z_gl[iCh];
150 delete[] wire_Ch[iCh];
151 delete[] xuv_Ch[iCh];
152 delete[] Chi2_ndf_seg[iCh];
153 delete[] Nhits_seg[iCh];
154 delete[] Nclust[iCh];
155 }
156 for (Int_t i = 0; i < 4; ++i) {
157 delete[] matrA[i];
158 delete[] matrb[i];
159 delete[] Amatr[i];
160 delete[] bmatr[i];
161 }
162 delete[] ChCent;
163 delete[] par_ab_Ch;
164 delete[] par_ab_seg;
165 delete[] sigma2_seg;
166 delete[] Wires_Ch;
167 delete[] clust_Ch;
168 delete[] XVU_Ch;
169 delete[] DigitsArray;
170 delete[] ClusterSize;
171 delete[] Coord_wire;
172 delete[] Coord_xuv;
173 delete[] XVU_coord;
174 delete[] Cluster_coord;
175 delete[] Coor_seg;
176 delete[] Cluster_seg;
177 delete[] kPln;
178 delete[] iw_Ch;
179 delete[] shift;
180 delete[] wire_Ch;
181 delete[] xuv_Ch;
182 delete[] Nhits_Ch;
183 delete[] Beam_wires_min;
184 delete[] Beam_wires_max;
185 delete[] Nseg_Ch;
186 delete[] Nbest_Ch;
187 delete[] ind_best_Ch;
188 delete[] best_Ch_gl;
189 delete[] Chi2_ndf_Ch;
190 delete[] Chi2_ndf_best_Ch;
191 delete[] XVU;
192 delete[] XVU_cl;
193 delete[] kZ_loc;
194 delete[] z_gl;
195
196 delete[] counter_pl;
197 delete[] ipl;
198 delete[] Nlay_w_wires;
199 delete[] sigm2;
200 delete[] z2;
201
202 delete[] matrA;
203 delete[] matrb;
204 delete[] Amatr;
205 delete[] bmatr;
206 delete[] Chi2_ndf_seg;
207 delete[] Nhits_seg;
208 delete[] Nclust;
209 delete[] Nbest_seg;
210 delete[] ChZ;
211 delete[] Zmid;
212}
213
214InitStatus BmnMwpcHitFinder::Init()
215{
216
217 if (!expData) {
218 if (fDebug)
219 cout << "BmnMwpcHitFinder::Init(): simulation data is reconstructed " << endl;
220 SetActive(kTRUE);
221 } else {
222 if (fDebug)
223 cout << " BmnMwpcHitFinder::Init() " << endl;
224 }
225
226 FairRootManager* ioman = FairRootManager::Instance();
227
228 if (expData) {
229 if (fDebug)
230 cout << " expData " << endl;
231 fBmnMwpcDigitArray = (TClonesArray*)ioman->GetObject(fInputBranchName);
232 if (!fBmnMwpcDigitArray) {
233 cout << "BmnMwpcHitFinder::Init(): branch " << fInputBranchName << " not found! Task will be deactivated"
234 << endl;
235 SetActive(kFALSE);
236 return kERROR;
237 }
238 } else {
239 if (fDebug)
240 cout << " !expData " << endl;
241 fBmnHitsArray = (TClonesArray*)ioman->GetObject(fInputBranchName);
242 if (fDebug)
243 cout << "fBmnHitsArray = " << fInputBranchName << "\n";
244 if (!fBmnHitsArray) {
245 cout << "BmnMwpcHitFinder::Init(): branch " << fInputBranchName << " not found! Task will be deactivated"
246 << endl;
247 SetActive(kFALSE);
248 return kERROR;
249 }
250 }
251
252 fBmnMwpcSegmentsArray =
253 new TClonesArray(fOutputBranchName); // fBmnMwpcSegmentsArray = new TClonesArray(fOutputBranchName.Data());
254 ioman->Register(fOutputBranchName.Data(), "MWPC", fBmnMwpcSegmentsArray, kTRUE);
255
256 // fMwpcGeo = new BmnMwpcGeometrySRC(fRunPeriod, fRunNumber);
257 cout << " fRunPeriod " << fRunPeriod << " fRunNumber " << fRunNumber << endl;
258 fMwpcGeometrySRC = new BmnMwpcGeometrySRC(fRunPeriod, fRunNumber);
259
260 kNChambers = fMwpcGeometrySRC->GetNChambers();
261 kNPlanes = fMwpcGeometrySRC->GetNPlanes();
262 kNWires = fMwpcGeometrySRC->GetNWires();
263
264 if (fDebug)
265 printf("C-P-W: %d %d %d\n", kNChambers, kNPlanes, kNWires);
266 // cout<<" MWPC runPeriod "<<fRunPeriod<<" fRunNumber "<<fRunNumber<<endl;
267
268 ChCent = new TVector3[kNChambers];
269 ChZ = new Float_t[kNChambers];
270 Zmid = new Float_t[kNChambers];
271
272 kChi2_Max = 20.;
273
274 for (int i = 0; i < kNChambers; ++i) {
275
276 if (fDebug) {
277 TH1D* h;
278
279 h = new TH1D(Form("Chisq_ndf_Ch%d", i), Form("Chisq_ndf_Ch%d", i), 100, 0., kChi2_Max);
280 fList.Add(h);
281 hChisq_ndf_Ch.push_back(h);
282 h = new TH1D(Form("Np_best_Ch%d", i), Form("Np_best_Ch%d", i), 6, 1.0, 7.0);
283 fList.Add(h);
284 hNp_best_Ch.push_back(h);
285 h = new TH1D(Form("Nbest_Ch%d", i), Form("Nbest_Ch%d", i), 6, 0.0, 6.0);
286 fList.Add(h);
287 hNbest_Ch.push_back(h);
288 h = new TH1D(Form("hpar_Ax_Ch%d", i), Form("par_Ax_Ch%d", i), 100, -0.5, 0.5);
289 fList.Add(h);
290 hpar_Ax_Ch.push_back(h);
291 h = new TH1D(Form("hpar_Bx_Ch%d", i), Form("par_Bx_Ch%d", i), 100, -10, 10);
292 fList.Add(h);
293 hpar_Bx_Ch.push_back(h);
294 h = new TH1D(Form("hpar_Ay_Ch%d", i), Form("par_Ay_Ch%d", i), 100, -0.5, 0.5);
295 fList.Add(h);
296 hpar_Ay_Ch.push_back(h);
297 h = new TH1D(Form("hpar_By_Ch%d", i), Form("par_By_Ch%d", i), 100, -10, 10);
298 fList.Add(h);
299 hpar_By_Ch.push_back(h);
300
301 h = new TH1D(Form("NFired_layers_Ch%d", i), Form("NFired_layers_Ch%d; NFired_layers; Events", i), 7, 0.0,
302 7.0);
303 fList.Add(h);
304 hNFired_layers_Ch.push_back(h);
305 h = new TH1D(Form("Num_layers_out_beam_Ch%d", i),
306 Form("Num_layers_out_beam_Ch%d; NFired_layers; Events", i), 7, 0.0, 7.0);
307 fList.Add(h);
308 hNum_layers_out_beam_Ch.push_back(h);
309 h = new TH1D(Form("Residuals_pl0_Ch%d", i), Form("Residuals_pl0_Ch%d", i), 10, 0., 10.);
310 fList.Add(h);
311 hResiduals_pl0_Ch.push_back(h);
312 h = new TH1D(Form("Residuals_pl1_Ch%d", i), Form("Residuals_pl1_Ch%d", i), 10, 0., 10.);
313 fList.Add(h);
314 hResiduals_pl1_Ch.push_back(h);
315 h = new TH1D(Form("Residuals_pl2_Ch%d", i), Form("Residuals_pl2_Ch%d", i), 10, 0., 10.);
316 fList.Add(h);
317 hResiduals_pl2_Ch.push_back(h);
318 h = new TH1D(Form("Residuals_pl3_Ch%d", i), Form("Residuals_pl3_Ch%d", i), 10, 0., 10.);
319 fList.Add(h);
320 hResiduals_pl3_Ch.push_back(h);
321 h = new TH1D(Form("Residuals_pl4_Ch%d", i), Form("Residuals_pl4_Ch%d", i), 10, 0., 10.);
322 fList.Add(h);
323 hResiduals_pl4_Ch.push_back(h);
324 h = new TH1D(Form("Residuals_pl5_Ch%d", i), Form("Residuals_pl5_Ch%d", i), 10, 0., 10.);
325 fList.Add(h);
326 hResiduals_pl5_Ch.push_back(h);
327 h = new TH1D(Form("occupancyXp%d", i), Form("occupancyXp%d", i), 250, -10., 10.);
328 fList.Add(h);
329 hoccupancyXp.push_back(h);
330 h = new TH1D(Form("occupancyUp%d", i), Form("occupancyUp%d", i), 250, -10., 10.);
331 fList.Add(h);
332 hoccupancyUp.push_back(h);
333 h = new TH1D(Form("occupancyVp%d", i), Form("occupancyVp%d", i), 250, -10., 10.);
334 fList.Add(h);
335 hoccupancyVp.push_back(h);
336 h = new TH1D(Form("occupancyXm%d", i), Form("occupancyXm%d", i), 250, -10., 10.);
337 fList.Add(h);
338 hoccupancyXm.push_back(h);
339 h = new TH1D(Form("occupancyUm%d", i), Form("occupancyUm%d", i), 250, -10., 10.);
340 fList.Add(h);
341 hoccupancyUm.push_back(h);
342 h = new TH1D(Form("occupancyVm%d", i), Form("occupancyVm%d", i), 250, -10., 10.);
343 fList.Add(h);
344 hoccupancyVm.push_back(h);
345 h = new TH1D(Form("firedWire_Ch%d", i), Form("firedWire_Ch%d", i), 100, 0., 100.);
346 fList.Add(h);
347 hfiredWire_Ch.push_back(h);
348 h = new TH1D(Form("ClusterSize_Ch%d", i), Form("ClusterSize_Ch%d", i), 15, 0., 15.);
349 fList.Add(h);
350 hClusterSize.push_back(h);
351
352 h = new TH1D(Form("fired_wire_Ch%d", i), Form("fired_wire_Ch%d;Wires;Events", i), kNWires * 6, 0.,
353 kNWires * 6);
354 fList.Add(h);
355 hfired_wire_Ch.push_back(h);
356
357 h = new TH1D(Form("WiresXm_Ch%d", i), Form("WiresXm_Ch%d", i), kNWires, 0., kNWires);
358 fList.Add(h);
359 hWiresXm.push_back(h);
360 h = new TH1D(Form("WiresVm_Ch%d", i), Form("WiresVm_Ch%d", i), kNWires, 0., kNWires);
361 fList.Add(h);
362 hWiresVm.push_back(h);
363 h = new TH1D(Form("WiresUp_Ch%d", i), Form("WiresUp_Ch%d", i), kNWires, 0., kNWires);
364 fList.Add(h);
365 hWiresUp.push_back(h);
366 h = new TH1D(Form("WiresXp_Ch%d", i), Form("WiresXp_Ch%d", i), kNWires, 0., kNWires);
367 fList.Add(h);
368 hWiresXp.push_back(h);
369 h = new TH1D(Form("WiresVp_Ch%d", i), Form("WiresVp_Ch%d", i), kNWires, 0., kNWires);
370 fList.Add(h);
371 hWiresVp.push_back(h);
372 h = new TH1D(Form("WiresUm_Ch%d", i), Form("WiresUm_Ch%d", i), kNWires, 0., kNWires);
373 fList.Add(h);
374 hWiresUm.push_back(h);
375
376 TH2D* h2;
377 h2 = new TH2D(Form("hY_X_Ch%d", i), Form("hY_X_Ch%d", i), 100, -20, 20, 100, -20, 20);
378 fList.Add(h2);
379 hY_X_Ch.push_back(h2);
380 h2 = new TH2D(Form("Event_display_Ch%d", i), Form("Event_display_Ch%d; Events; Wires", i), 100, 0., 700.,
381 100, 0., 100.);
382 fList.Add(h2);
383 hEvent_display_Ch.push_back(h2);
384 h2 = new TH2D(Form("all_pl_time_wire_Ch%d", i), Form("time_wire_allplane_Ch%d; Wires; Time", i),
385 kNWires * 6, 0., kNWires * 6, 100, 0, 500.);
386 fList.Add(h2);
387 htime_wire_Ch.push_back(h2);
388 }
389
390 ChCent[i] = fMwpcGeometrySRC->GetChamberCenter(i);
391 ChZ[i] = ChCent[i].Z();
392 }
393
394 if (fDebug) {
395
396 for (int i = 0; i < kNChambers * kNPlanes; ++i) {
397 TH1D* h;
398 h = new TH1D(Form("Time%d", i), Form("Time%d", i), 500, 0., 500.);
399 fList.Add(h);
400 hTime.push_back(h);
401 h = new TH1D(Form("Occupancy%d", i), Form("Occupancy%d", i), 100, 0., 100.);
402 fList.Add(h);
403 hOccupancy.push_back(h);
404 }
405
406 hNtrMC = new TH1D("NtrMC", "NtrMC: ch2 || ch3", 10, 0, 10);
407 fList.Add(hNtrMC);
408
409 hNp_MCtrue_ch2 = new TH1D("Np_MCtrue_ch2", "Np_MCtrue_ch2", 7, 0, 7);
410 fList.Add(hNp_MCtrue_ch2);
411 hNp_MCtrue_ch3 = new TH1D("Np_MCtrue_ch3", "Np_MCtrue_ch3", 7, 0, 7);
412 fList.Add(hNp_MCtrue_ch3);
413
414 hx_target = new TH1D("x_target", "x_target", 100, -30, 30);
415 fList.Add(hx_target);
416 hy_target = new TH1D("y_target", "y_target", 100, -30, 30);
417 fList.Add(hy_target);
418
419 for (int i = 0; i < kNChambers; ++i) {
420 TH1D* h;
421
422 h = new TH1D(Form("hx_target_best%d", i), Form("hx_target_best%d", i), 100, -30, 30);
423 fList.Add(h);
424 hx_target_best.push_back(h);
425 h = new TH1D(Form("hy_target_best%d", i), Form("hy_target_best%d", i), 100, -30, 30);
426 fList.Add(h);
427 hy_target_best.push_back(h);
428
429 h = new TH1D(Form("hAx_mc_ch%d", i), Form("hAx_mc_ch%d", i), 100, -0.1, 0.1);
430 fList.Add(h);
431 hAx_mc_ch.push_back(h);
432 h = new TH1D(Form("hBx_mc_ch%d", i), Form("hBx_mc_ch%d", i), 100, -10, 10);
433 fList.Add(h);
434 hBx_mc_ch.push_back(h);
435 h = new TH1D(Form("hAy_mc_ch%d", i), Form("hAy_mc_ch%d", i), 100, -0.1, 0.1);
436 fList.Add(h);
437 hAy_mc_ch.push_back(h);
438 h = new TH1D(Form("hBy_mc_ch%d", i), Form("hBy_mc_ch%d", i), 100, -10, 10);
439 fList.Add(h);
440 hBy_mc_ch.push_back(h);
441
442 h = new TH1D(Form("hdAx_mc_SegCh%d", i), Form("dAx_mc_SegCh%d", i), 100, -0.5, 0.5);
443 fList.Add(h);
444 hdAx_mc_SegCh.push_back(h);
445 h = new TH1D(Form("hdX_mc_SegCh%d", i), Form("dX_mc_SegCh%d", i), 100, -0.5, 0.5);
446 fList.Add(h);
447 hdX_mc_SegCh.push_back(h);
448 h = new TH1D(Form("hdAy_mc_SegCh%d", i), Form("dAy_mc_SegCh%d", i), 100, -0.5, 0.5);
449 fList.Add(h);
450 hdAy_mc_SegCh.push_back(h);
451 h = new TH1D(Form("hdY_mc_SegCh%d", i), Form("dY_mc_SegCh%d", i), 100, -0.5, 0.5);
452 fList.Add(h);
453 hdY_mc_SegCh.push_back(h);
454
455 h = new TH1D(Form("hdAx_mc_Seg_deltaCh%d", i), Form("dAx_mc_Seg_deltaCh%d", i), 100, -0.5, 0.5);
456 fList.Add(h);
457 hdAx_mc_Seg_deltaCh.push_back(h);
458 h = new TH1D(Form("hdX_mc_Seg_deltaCh%d", i), Form("dX_mc_Seg_deltaCh%d", i), 100, -0.5, 0.5);
459 fList.Add(h);
460 hdX_mc_Seg_deltaCh.push_back(h);
461 h = new TH1D(Form("hdAy_mc_Seg_deltaCh%d", i), Form("dAy_mc_Seg_deltaCh%d", i), 100, -0.5, 0.5);
462 fList.Add(h);
463 hdAy_mc_Seg_deltaCh.push_back(h);
464 h = new TH1D(Form("hdY_mc_Seg_deltaCh%d", i), Form("dY_mc_Seg_deltaCh%d", i), 100, -0.5, 0.5);
465 fList.Add(h);
466 hdY_mc_Seg_deltaCh.push_back(h);
467 }
468
469 hYvsX_mc_ch2 = new TH2D("hYvsX_mc_ch2", "YvsX_mc_ch2", 100, -20, 20, 100, -20, 20);
470 fList.Add(hYvsX_mc_ch2);
471 hYvsX_mc_ch3 = new TH2D("hYvsX_mc_ch3", "YvsX_mc_ch3", 100, -20, 20, 100, -20, 20);
472 fList.Add(hYvsX_mc_ch3);
473
474 hHalfDeadCh = new TH1D("hHalfDeadCh", "HalfDeadCh - Ch0+Ch1", 5, 0, 5);
475 hHalfDeadCh->GetXaxis()->SetBinLabel(1, "Ch0");
476 hHalfDeadCh->GetXaxis()->SetBinLabel(2, "Ch1");
477 hHalfDeadCh->GetXaxis()->SetBinLabel(3, "Ch2");
478 hHalfDeadCh->GetXaxis()->SetBinLabel(4, "Ch3");
479 hHalfDeadCh->GetXaxis()->SetBinLabel(5, "Ch2&Ch3");
480 fList.Add(hHalfDeadCh);
481
482 hDen_mc = new TH1D("hDen_mc", "Den_mc", 6, 0, 6);
483 fList.Add(hDen_mc);
484 hNum_mc = new TH1D("hNum_mc", "Num_mc", 6, 0, 6);
485 fList.Add(hNum_mc);
486 hEff_mc = new TH1D("hEff_mc", "Eff_mc", 6, 0, 6);
487 hEff_mc->GetXaxis()->SetBinLabel(1, "Ch2");
488 hEff_mc->GetXaxis()->SetBinLabel(2, "Ch3");
489 hEff_mc->GetXaxis()->SetBinLabel(3, "Ch2(3coord)");
490 hEff_mc->GetXaxis()->SetBinLabel(4, "Ch3(3coord)");
491 hEff_mc->GetXaxis()->SetBinLabel(5, "Ch2(3c/3cmc)");
492 hEff_mc->GetXaxis()->SetBinLabel(6, "Ch3(3c/3cmc)");
493 hDen_mcreaction2 = new TH1D("hDen_mcreaction2", "hDen_mcreaction2", 1, 0, 1);
494 fList.Add(hDen_mcreaction2);
495 hNum_mcreaction2 = new TH1D("hNum_mcreaction2", "hNum_mcreaction2", 1, 0, 1);
496 fList.Add(hNum_mcreaction2);
497 hEff_mcreaction2 = new TH1D("hEff_mcreaction2", "hEff_mcreaction2", 1, 0, 1);
498 fList.Add(hEff_mcreaction2);
499 hDen_mcreaction3 = new TH1D("hDen_mcreaction3", "hDen_mcreaction3", 1, 0, 1);
500 fList.Add(hDen_mcreaction3);
501 hNum_mcreaction3 = new TH1D("hNum_mcreaction3", "hNum_mcreaction3", 1, 0, 1);
502 fList.Add(hNum_mcreaction3);
503 hEff_mcreaction3 = new TH1D("hEff_mcreaction3", "hEff_mcreaction3", 1, 0, 1);
504 fList.Add(hEff_mcreaction3);
505
506 fList.Add(hEff_mc);
507
508 } // if (fDebug)
509
510 // Segment finding cuts
511 kMinHits = 4; //// for alignment kMinHits = 6;
512 // if(!expData) kMinHits = 5;
513 kMinHits_before_target = 5;
514 kmaxSeg = 15;
515 kChMaxAllWires = 200;
516 // Constants
517 dw = fMwpcGeometrySRC->GetWireStep(); // 0.25; // [cm] // wires step
518 dw_half = 0.5 * dw;
519 sq3 = sqrt(3.);
520 sq12 = sqrt(12.);
521 sigma = dw / sq12;
522 kMiddlePl = 47.25; // Center of wires plane
523 // Matrices
524 matrA = new Double_t*[4];
525 matrb = new Double_t*[4];
526 Amatr = new Double_t*[4];
527 bmatr = new Double_t*[4];
528
529 for (Int_t ii = 0; ii < 4; ii++) {
530 Amatr[ii] = new Double_t[4];
531 bmatr[ii] = new Double_t[4];
532 }
533
534 // Arrays
535 Nlay_w_wires = new Int_t[kNChambers];
536 kPln = new Int_t*[kNChambers];
537 iw_Ch = new Int_t*[kNChambers];
538 shift = new Float_t*[kNChambers];
539 wire_Ch = new Int_t**[kNChambers];
540 xuv_Ch = new Float_t**[kNChambers];
541 Wires_Ch = new Int_t**[kNChambers];
542 clust_Ch = new Int_t**[kNChambers];
543 XVU_Ch = new Float_t**[kNChambers];
544 Nhits_Ch = new Int_t*[kNChambers];
545 Beam_wires_min = new Int_t*[kNChambers];
546 Beam_wires_max = new Int_t*[kNChambers];
547 Nseg_Ch = new Int_t[kNChambers];
548 Nbest_Ch = new Int_t[kNChambers];
549 ind_best_Ch = new Int_t*[kNChambers];
550 best_Ch_gl = new Int_t*[kNChambers];
551 Chi2_ndf_Ch = new Double_t*[kNChambers];
552 Chi2_ndf_best_Ch = new Double_t*[kNChambers];
553 par_ab_Ch = new Double_t**[kNChambers];
554 XVU = new Float_t*[kNChambers];
555 XVU_cl = new Float_t*[kNChambers];
556 kZ_loc = new Float_t*[kNChambers];
557 z_gl = new Float_t*[kNChambers];
558 sigm2 = new Float_t[kNPlanes];
559 ipl = new Int_t[kNPlanes];
560 counter_pl = new Int_t[kNPlanes];
561 z2 = new Float_t[kNPlanes];
562 DigitsArray = new Double_t**[kNChambers];
563 ClusterSize = new Int_t**[kNChambers];
564 Nclust = new Int_t*[kNChambers];
565 Coord_wire = new Double_t**[kNChambers];
566 Coord_xuv = new Double_t**[kNChambers];
567 XVU_coord = new Double_t**[kNChambers];
568 Cluster_coord = new Double_t**[kNChambers];
569 Nhits_seg = new Int_t*[kNChambers];
570 Chi2_ndf_seg = new Double_t*[kNChambers];
571 Coor_seg = new Double_t**[kNChambers];
572 Cluster_seg = new Double_t**[kNChambers];
573 par_ab_seg = new Double_t**[kNChambers];
574 sigma2_seg = new Double_t***[kNChambers];
575 Nbest_seg = new Int_t[kNChambers];
576
577 for (Int_t i = 0; i < kNChambers; ++i) {
578
579 if (i == 0 || i == 2) {
580 Zmid[i] = (ChZ[i] - ChZ[i + 1]) * 0.5;
581 }
582 if (i == 1 || i == 3) {
583 Zmid[i] = (ChZ[i - 1] - ChZ[i]) * -0.5;
584 }
585 if (fDebug)
586 printf("Chamber %d Z: %f Zmid: %f\n", i, ChZ[i], Zmid[i]);
587
588 kPln[i] = new Int_t[kNPlanes];
589 iw_Ch[i] = new Int_t[kNPlanes];
590 kZ_loc[i] = new Float_t[kNPlanes];
591 z_gl[i] = new Float_t[kNPlanes];
592 Nhits_Ch[i] = new Int_t[kBig];
593 shift[i] = new Float_t[4];
594 wire_Ch[i] = new Int_t*[kNWires];
595 xuv_Ch[i] = new Float_t*[kNWires];
596 Wires_Ch[i] = new Int_t*[kNPlanes];
597 clust_Ch[i] = new Int_t*[kNPlanes];
598 XVU_Ch[i] = new Float_t*[kNPlanes];
599 par_ab_Ch[i] = new Double_t*[4];
600 XVU[i] = new Float_t[kNPlanes];
601 XVU_cl[i] = new Float_t[kNPlanes];
602 ind_best_Ch[i] = new Int_t[kmaxSeg];
603 best_Ch_gl[i] = new Int_t[kmaxSeg];
604 Chi2_ndf_Ch[i] = new Double_t[kBig];
605 Chi2_ndf_best_Ch[i] = new Double_t[kmaxSeg];
606 DigitsArray[i] = new Double_t*[kNPlanes];
607 ClusterSize[i] = new Int_t*[kNPlanes];
608 Nclust[i] = new Int_t[kNPlanes];
609 Coord_wire[i] = new Double_t*[kNPlanes];
610 Coord_xuv[i] = new Double_t*[kNPlanes];
611 XVU_coord[i] = new Double_t*[kNPlanes];
612 Cluster_coord[i] = new Double_t*[kNPlanes];
613 Nhits_seg[i] = new Int_t[kBig];
614 Chi2_ndf_seg[i] = new Double_t[kBig];
615 Coor_seg[i] = new Double_t*[kNPlanes];
616 Cluster_seg[i] = new Double_t*[kNPlanes];
617 par_ab_seg[i] = new Double_t*[4];
618 sigma2_seg[i] = new Double_t**[kBig];
619 Beam_wires_min[i] = new Int_t[kNPlanes];
620 Beam_wires_max[i] = new Int_t[kNPlanes];
621
622 shift[i][0] = fMwpcGeometrySRC->GetTx(i);
623 shift[i][2] = fMwpcGeometrySRC->GetTy(i);
624 shift[i][1] = ChCent[i].X();
625 shift[i][3] = ChCent[i].Y();
626
627 for (int iWire = 0; iWire < kNWires; iWire++) {
628 wire_Ch[i][iWire] = new Int_t[kNPlanes];
629 xuv_Ch[i][iWire] = new Float_t[kNPlanes];
630 }
631 for (int iPla = 0; iPla < kNPlanes; ++iPla) {
632 Wires_Ch[i][iPla] = new Int_t[kBig];
633 clust_Ch[i][iPla] = new Int_t[kBig];
634 XVU_Ch[i][iPla] = new Float_t[kBig];
635 DigitsArray[i][iPla] = new Double_t[kNWires];
636 ClusterSize[i][iPla] = new Int_t[kBig];
637 Coord_wire[i][iPla] = new Double_t[kBig];
638 Coord_xuv[i][iPla] = new Double_t[kBig];
639 XVU_coord[i][iPla] = new Double_t[kBig];
640 Cluster_coord[i][iPla] = new Double_t[kBig];
641 Coor_seg[i][iPla] = new Double_t[kBig];
642 Cluster_seg[i][iPla] = new Double_t[kBig];
643
644 if (fRunPeriod == 6 || (fRunPeriod == 7 && fRunNumber > 3588)) // if ( i == 2 || i == 3)
645 {
646 kPln[i][0] = 4;
647 kPln[i][1] = 5;
648 kPln[i][2] = 0;
649 kPln[i][3] = 1;
650 kPln[i][4] = 2;
651 kPln[i][5] = 3; //{4,5,0,1,2,3, 7,11,6,10,9,8, 0,0,0,0,0,0, 0,0,0,0,0,0};
652 kZ_loc[i][iPla] = -0.5 + iPla;
653 if (iPla == 4)
654 kZ_loc[i][iPla] = -2.5;
655 if (iPla == 5)
656 kZ_loc[i][iPla] = -1.5;
657 }
658 if ((fRunPeriod == 7 && fRunNumber <= 3588) || (fRunPeriod == 8)) // SRC
659 {
660 if (i == 0) {
661 kPln[i][0] = 5;
662 kZ_loc[i][0] = -1.5;
663 kPln[i][1] = 0;
664 kZ_loc[i][1] = -0.5;
665 kPln[i][2] = 1;
666 kZ_loc[i][2] = 0.5;
667 kPln[i][3] = 2;
668 kZ_loc[i][3] = 1.5;
669 kPln[i][4] = 3;
670 kZ_loc[i][4] = 2.5;
671 kPln[i][5] = 4;
672 kZ_loc[i][5] = -2.5;
673 }
674 if (i == 1) {
675 kPln[i][0] = 1;
676 kZ_loc[1][0] = -1.5;
677 kPln[i][1] = 0;
678 kZ_loc[i][1] = -2.5;
679 kPln[i][2] = 5;
680 kZ_loc[i][2] = 2.5;
681 kPln[i][3] = 4;
682 kZ_loc[i][3] = 1.5;
683 kPln[i][4] = 3;
684 kZ_loc[i][4] = 0.5;
685 kPln[i][5] = 2;
686 kZ_loc[i][5] = -0.5;
687 }
688 if (fRunPeriod == 7 && (i == 2 || i == 3)) {
689 kPln[i][0] = 4;
690 kPln[i][1] = 5;
691 kPln[i][2] = 0;
692 kPln[i][3] = 1;
693 kPln[i][4] = 2;
694 kPln[i][5] = 3;
695
696 kZ_loc[i][iPla] = -0.5 + iPla;
697 if (iPla == 4) {
698 kZ_loc[i][iPla] = -2.5;
699 }
700 if (iPla == 5) {
701 kZ_loc[i][iPla] = -1.5;
702 }
703 } // fRunPeriod == 7 ch 2 3
704 if (fRunPeriod == 8 && (i == 2 || i == 3)) {
705 kPln[i][0] = 3;
706 kZ_loc[i][0] = 0.5;
707 kPln[i][1] = 2;
708 kZ_loc[i][1] = -0.5;
709 kPln[i][2] = 1;
710 kZ_loc[i][2] = -1.5;
711 kPln[i][3] = 0;
712 kZ_loc[i][3] = -2.5;
713 kPln[i][4] = 5;
714 kZ_loc[i][4] = 2.5;
715 kPln[i][5] = 4;
716 kZ_loc[i][5] = 1.5;
717 } // fRunPeriod == 8 ch 2 3
718
719 } // SRC
720
721 if (fRunPeriod == 6)
722 Z0_SRC = 0.;
723 if (fRunPeriod == 7)
724 Z0_SRC = -647.476;
725 if (fRunPeriod == 8)
726 Z0_SRC = -574.91;
727
728 if (!expData) { // MC
729
730 kPln[i][0] = 4;
731 kPln[i][1] = 5;
732 kPln[i][2] = 0;
733 kPln[i][3] = 1;
734 kPln[i][4] = 2;
735 kPln[i][5] = 3;
736
737 kZ_loc[i][iPla] = -0.5 + iPla;
738 if (iPla == 4) {
739 kZ_loc[i][iPla] = -2.5;
740 }
741 if (iPla == 5) {
742 kZ_loc[i][iPla] = -1.5;
743 }
744 } // MC
745
746 } // iPla
747
748 for (int ii = 0; ii < 4; ++ii) { // 4 parameters: tan(x), tan(y), x ,y
749 par_ab_Ch[i][ii] = new Double_t[kBig];
750 par_ab_seg[i][ii] = new Double_t[kBig];
751 } // 4
752 } // kChamber
753 for (int ii = 0; ii < 4; ++ii) {
754 matrA[ii] = new Double_t[4];
755 matrb[ii] = new Double_t[4];
756 }
757 for (Int_t i = 0; i < kNChambers; ++i) {
758 for (Int_t j = 0; j < kBig; ++j) {
759 sigma2_seg[i][j] = new Double_t*[4];
760 }
761 }
762 for (Int_t i = 0; i < kNChambers; ++i) {
763 for (Int_t j = 0; j < kBig; ++j) {
764 for (Int_t k = 0; k < 4; ++k) {
765 sigma2_seg[i][j][k] = new Double_t[4];
766 }
767 }
768 }
769
770 // if (fRunPeriod == 6) {
771 // kPln[1][0] = 1; kZ_loc[1][0] =-2.5;// kPln[3][0] = 1;//run6-II
772 // kPln[1][3] = 4; kZ_loc[1][3] = 0.5;// kPln[3][3] = 4;//
773 // }
774
775 // if (fDebug) printf("Chamber Plane kZ_loc z_gl\n");
776 for (Int_t i = 0; i < kNChambers; ++i) {
777 for (int ii = 0; ii < kNPlanes; ii++) {
778 z_gl[i][ii] = Zmid[i] + kZ_loc[i][ii];
779 // if (fDebug) printf("%5d %5d %8.4f %8.4f\n", i, ii, kZ_loc[i][ii], z_gl[i][ii]);
780 if (i < 2)
781 Beam_wires_min[i][ii] = 0;
782 }
783 }
784 if ((fRunPeriod == 7 && fRunNumber <= 3588) || (fRunPeriod == 8)) // SRC
785 {
786 // beam area //run7
787 Beam_wires_min[2][0] = 42;
788 Beam_wires_max[2][0] = 58; // x-
789 Beam_wires_min[2][1] = 21;
790 Beam_wires_max[2][1] = 34; // v-
791 Beam_wires_min[2][2] = 21;
792 Beam_wires_max[2][2] = 37; // u+
793 Beam_wires_min[2][3] = 41;
794 Beam_wires_max[2][3] = 57; // x+
795 Beam_wires_min[2][4] = 58;
796 Beam_wires_max[2][4] = 73; // v+
797 Beam_wires_min[2][5] = 58;
798 Beam_wires_max[2][5] = 72; // u-
799
800 Beam_wires_min[3][0] = 40;
801 Beam_wires_max[3][0] = 62;
802 Beam_wires_min[3][1] = 15;
803 Beam_wires_max[3][1] = 35;
804 Beam_wires_min[3][2] = 19;
805 Beam_wires_max[3][2] = 38;
806 Beam_wires_min[3][3] = 40;
807 Beam_wires_max[3][3] = 60;
808 Beam_wires_min[3][4] = 56;
809 Beam_wires_max[3][4] = 77;
810 Beam_wires_min[3][5] = 58;
811 Beam_wires_max[3][5] = 81;
812 }
813
814 // fBmnEvQuality = (TClonesArray*) ioman->GetObject(fBmnEvQualityBranchName);
815 cout << " end init " << endl;
816 return kSUCCESS;
817}
818
819void BmnMwpcHitFinder::Exec(Option_t* opt)
820{
821 if (!IsActive())
822 return;
823 clock_t tStart = clock();
824 PrepareArraysToProcessEvent();
825
826 if (fDebug)
827 cout << "\n======================== MWPC hit finder exec started =====================\n" << endl;
828 if (fDebug)
829 printf("Event number: %d\n", fEventNo++);
830
831 // ---------- Digi-file reading-------------------------------------
832 ReadWires(DigitsArray, iw_Ch, vec_points);
833
834 //--------- Big chambers cycle -> looking for track-segments ---------
835 for (Int_t iChamber = 0; iChamber < kNChambers; iChamber++) {
836 Int_t counter = 0;
837
838 for (Int_t iplane = 0; iplane < kNPlanes; iplane++) {
839 counter += iw_Ch[iChamber][iplane];
840 } // iplane
841
842 if (counter < kMinHits || counter > kChMaxAllWires)
843 continue; // too many or too few wires per Chamber
844
845 if (fDebug)
846 cout << "-- Clustering --" << endl;
847 Clustering(iChamber, ClusterSize, DigitsArray, Coord_wire, Coord_xuv, Nclust,
848 counter_pl); // Cluster production function
849
850 if (fDebug)
851 cout << "-- SegmentFinder3coord --" << endl;
852 for (Int_t iCase = 1; iCase < 9; iCase++) {
853 SegmentFinder(iChamber, Nclust, Coord_xuv, ClusterSize, Nseg_Ch, XVU_coord, Cluster_coord, Nhits_Ch,
854 kMinHits, iCase, kBig); // Combinatorial segment selection
855 }
856 if (fDebug)
857 cout << "-- after SegmentFinder3coord: Nseg_[" << iChamber << "]= " << Nseg_Ch[iChamber] << endl;
858
859 ProcessSegments(iChamber, Nseg_Ch, XVU_coord, Cluster_coord, Nhits_Ch, kZ_loc, kMinHits, sigma, kChi2_Max,
860 Nhits_seg, Chi2_ndf_seg, Coor_seg, Cluster_seg, par_ab_seg, Nbest_seg, Nlay_w_wires,
861 sigma2_seg); // Segment fitting & finding the best track-candidate
862
863 if (fDebug) {
864 cout << " Nbest_seg[" << iChamber << "]= " << Nbest_seg[iChamber] << endl;
865 for (Int_t iseg = 0; iseg < Nbest_seg[iChamber]; iseg++) {
866 cout << " iseg " << iseg << " Nhits_Ch " << Nhits_seg[iChamber][iseg] << " Chi2 "
867 << Chi2_ndf_seg[iChamber][iseg] << " Ax " << par_ab_seg[iChamber][0][iseg] << " bx "
868 << par_ab_seg[iChamber][1][iseg] << " by " << par_ab_seg[iChamber][3][iseg] << endl;
869 }
870 }
871
872 if (fDebug)
873 cout << "-- SegmentFinder2coord --" << endl;
874 if (iChamber > 1) {
875 for (Int_t iCase = 1; iCase < 4; iCase++) {
876 SegmentFinder2coord(iChamber, Nclust, Coord_xuv, ClusterSize, Nseg_Ch, XVU_coord, Cluster_coord,
877 Nhits_Ch, kMinHits, iCase, kBig, Nbest_seg,
878 Coor_seg); // Combinatorial segment selection
879 }
880 ProcessSegments2coord(iChamber, Nseg_Ch, XVU_coord, Cluster_coord, Nhits_Ch, kZ_loc, kMinHits, sigma,
881 kChi2_Max, Nhits_seg, Chi2_ndf_seg, Coor_seg, Cluster_seg, par_ab_seg, Nbest_seg,
882 Nlay_w_wires, sigma2_seg); //
883 }
884 if (fDebug)
885 cout << "--after ProcessSegments: Nbest[" << iChamber << "] " << Nbest_seg[iChamber] << endl;
886 if (fDebug && Nbest_seg[iChamber] > 0)
887 hNbest_Ch.at(iChamber)->Fill(Nbest_seg[iChamber]);
888
889 if (fDebug) {
890 cout << " ----------------" << endl;
891 cout << " print all segments Nbest_seg[" << iChamber << "]= " << Nbest_seg[iChamber] << endl;
892 for (Int_t iseg = 0; iseg < Nbest_seg[iChamber]; iseg++) {
893 cout << " iseg " << iseg << " Nhits_Ch " << Nhits_seg[iChamber][iseg] << " Chi2 "
894 << Chi2_ndf_seg[iChamber][iseg] << " Ax " << par_ab_seg[iChamber][0][iseg] << " bx "
895 << par_ab_seg[iChamber][1][iseg] << " by " << par_ab_seg[iChamber][3][iseg] << endl;
896 if (fabs(par_ab_seg[iChamber][1][iseg]) > 10. && fabs(par_ab_seg[iChamber][3][iseg]) < 4.) {
897 cout << " !!! out of hexagon Nhits_Ch: " << Nhits_seg[iChamber][iseg] << endl;
898 }
899 hChisq_ndf_Ch.at(iChamber)->Fill(Chi2_ndf_seg[iChamber][iseg]);
900 if (Nhits_seg[iChamber][iseg] > 0)
901 hNp_best_Ch.at(iChamber)->Fill(Nhits_seg[iChamber][iseg]);
902 }
903 }
904
905 //--efficiency calculation for MC--
906 if (!expData)
907 MCefficiencyCalculation(iChamber, vec_points, par_ab_seg, Nbest_seg);
908 //-- Alignment--
909 SegmentParamAlignment(iChamber, Nbest_seg, par_ab_seg, shift);
910 } // iChamber
911
912 SegmentsStoring(Nbest_seg, par_ab_seg, Chi2_ndf_seg, Nhits_seg, Coor_seg, Cluster_seg, sigma2_seg);
913
914 clock_t tFinish = clock();
915 workTime += ((Float_t)(tFinish - tStart)) / CLOCKS_PER_SEC;
916 if (fDebug)
917 cout << "\n======================== MWPC hit finder exec finished ====================" << endl;
918}
919
920void BmnMwpcHitFinder::SegmentsStoring(Int_t* Nbest,
921 Double_t*** par_ab,
922 Double_t** Chi2_ndf,
923 Int_t** Nhits,
924 Double_t*** Coor,
925 Double_t*** Cluster,
926 Double_t**** sigma2)
927{
928 if (fDebug)
929 cout << "--SegmentsStoring--" << endl;
930 vector<Double_t> vtmpCoor;
931 vector<Double_t> vtmpClust;
932 for (Int_t iChamber = 0; iChamber < kNChambers; iChamber++) {
933 for (Int_t ise = 0; ise < Nbest[iChamber]; ise++) {
934
935 if (fDebug) {
936 hpar_Ax_Ch.at(iChamber)->Fill(par_ab[iChamber][0][ise]);
937 hpar_Bx_Ch.at(iChamber)->Fill(par_ab[iChamber][1][ise]);
938 hpar_Ay_Ch.at(iChamber)->Fill(par_ab[iChamber][2][ise]);
939 hpar_By_Ch.at(iChamber)->Fill(par_ab[iChamber][3][ise]);
940 hY_X_Ch.at(iChamber)->Fill(par_ab[iChamber][1][ise], par_ab[iChamber][3][ise]);
941 }
942
943 // if (Nhits[iChamber][ise] > 3) {
944
945 BmnMwpcSegment* pSeg =
946 new ((*fBmnMwpcSegmentsArray)[fBmnMwpcSegmentsArray->GetEntriesFast()]) BmnMwpcSegment();
947 pSeg->SetChi2(Chi2_ndf[iChamber][ise]);
948 pSeg->SetNHits(Nhits[iChamber][ise]);
949 pSeg->SetFlag(ise);
950
951 vtmpCoor.clear();
952 vtmpClust.clear();
953
954 for (Int_t i1 = 0; i1 < 6; i1++) {
955 vtmpCoor.push_back(Coor[iChamber][i1][ise]);
956 vtmpClust.push_back(Cluster[iChamber][i1][ise]);
957 }
958 pSeg->SetClust(vtmpClust);
959 pSeg->SetCoord(vtmpCoor);
960
961 FairTrackParam pSegParams;
962 pSegParams.SetPosition(TVector3(par_ab[iChamber][1][ise], par_ab[iChamber][3][ise], ChZ[iChamber]));
963 pSegParams.SetTx(par_ab[iChamber][0][ise]);
964 pSegParams.SetTy(par_ab[iChamber][2][ise]);
965 for (Int_t i1 = 0; i1 < 4; i1++) {
966 for (Int_t j1 = 0; j1 < 4; j1++) {
967 // cout<<" sigma2_seg["<<iChamber<<"]["<<ise<<"]["<<i1<<"]["<<j1<<"]
968 // "<<sigma2_seg[iChamber][ise][i1][j1]<<endl;
969 pSegParams.SetCovariance(i1, j1, sigma2[iChamber][ise][i1][j1]);
970 }
971 }
972 pSeg->SetParamFirst(pSegParams);
973 //}//if
974 } // ise
975 } //[iChamber]
976 //--------------------------------------------------------------------
977}
978
979void BmnMwpcHitFinder::MCefficiencyCalculation(Int_t iCh,
980 vector<MC_points>& vec,
981 Double_t*** par_ab_seg_,
982 Int_t* Nbest_seg_)
983{
984 if (!expData) {
985
986 if (fDebug)
987 cout << " ---MC tracks association--" << endl;
988 // ax, bx, ay, by
989 Double_t delta2[4] = {-99., -999., -99., -999.};
990 Double_t delta3[4] = {-99., -999., -99., -999.};
991
992 // ax, bx, ay, by
993 Double_t sig[4] = {0.04, 0.08, 0.05, 0.08};
994
995 Double_t dmatch = 0.;
996 Double_t dmc_match[kNChambers][kMaxMC];
997 Int_t mc_tr_assoc[kNChambers][kMaxMC];
998
999 for (Int_t i = 0; i < kNChambers; i++) {
1000 for (Int_t j = 0; j < kMaxMC; j++) {
1001 dmc_match[i][j] = 1000.;
1002 mc_tr_assoc[i][j] = -1;
1003 }
1004 }
1005 Double_t dax = -999.;
1006 Double_t day = -999.;
1007 Double_t dx = -999.;
1008 Double_t dy = -999.;
1009
1010 if (fDebug)
1011 cout << " Hit: Nmctr = " << vec.size() << " Nreco[" << iCh << "] = " << Nbest_seg_[iCh] << endl;
1012
1013 Bool_t wasLi7_ch2 = 0;
1014 Bool_t wasHe4_ch2 = 0;
1015 Bool_t wasLi7_ch3 = 0;
1016 Bool_t wasHe4_ch3 = 0;
1017 Int_t Nassoc2 = 0;
1018 Int_t Nassoc3 = 0;
1019
1020 for (size_t itr = 0; itr < vec.size(); itr++) { // mc_tr
1021
1022 if (fDebug)
1023 cout << " Np2 " << vec.at(itr).Np2 << " Np3 " << vec.at(itr).Np3 << endl;
1024
1025 if (vec.at(itr).Np2 >= 4 && vec.at(itr).Pdg == PDG_Li7)
1026 wasLi7_ch2 = 1;
1027 if (vec.at(itr).Np3 >= 4 && vec.at(itr).Pdg == PDG_Li7)
1028 wasLi7_ch3 = 1;
1029 if (vec.at(itr).Np2 >= 4 && vec.at(itr).Pdg == PDG_He4)
1030 wasHe4_ch2 = 1;
1031 if (vec.at(itr).Np3 >= 4 && vec.at(itr).Pdg == PDG_He4)
1032 wasHe4_ch3 = 1;
1033
1034 if (iCh == 2 && fDebug)
1035 hNp_MCtrue_ch2->Fill(vec.at(itr).Np2);
1036 if (iCh == 3 && fDebug)
1037 hNp_MCtrue_ch3->Fill(vec.at(itr).Np3);
1038 //---MC Eff ---
1039 //---Den
1040 if (vec.at(itr).Np2 >= 4) {
1041 // if (vec.at(itr).Np2 >= 4 && vec.at(itr).xWas2 && vec.at(itr).uWas2 && vec.at(itr).vWas2)
1042 if (iCh == 2 && fDebug) {
1043 hDen_mc->Fill(0);
1044 hDen_mc->Fill(2);
1045 if (vec.at(itr).xWas2 && vec.at(itr).uWas2 && vec.at(itr).vWas2) {
1046 hDen_mc->Fill(4);
1047 }
1048 cout << " Den_mcPC2 " << endl;
1049 }
1050 }
1051 // if (fDebug) cout<<" xWas3 "<<vec.at(itr).xWas3<<" uWas3 "<<vec.at(itr).uWas3<<" vWas3
1052 // "<<vec.at(itr).vWas3<<endl;
1053 if (vec.at(itr).Np3 >= 4) {
1054 // if (vec.at(itr).xWas3 && vec.at(itr).uWas3 && vec.at(itr).vWas3
1055 if (iCh == 3 && fDebug) {
1056 hDen_mc->Fill(1);
1057 hDen_mc->Fill(3);
1058 if (vec.at(itr).xWas3 && vec.at(itr).uWas3 && vec.at(itr).vWas3) {
1059 hDen_mc->Fill(5);
1060 }
1061 cout << " Den_mcPC3 " << endl;
1062 }
1063 }
1064
1065 if (wasLi7_ch2 && wasHe4_ch2) {
1066 if (fDebug)
1067 hDen_mcreaction2->Fill(0);
1068 }
1069 if (wasLi7_ch3 && wasHe4_ch3) {
1070 if (fDebug)
1071 hDen_mcreaction3->Fill(0);
1072 }
1073 // for (Int_t iChamber = 2; iChamber < kNChambers; iChamber++) {
1074
1075 if (iCh == 2) {
1076 for (Int_t iseg = 0; iseg < Nbest_seg_[iCh]; iseg++) {
1077 delta2[0] = vec.at(itr).param2[0] - par_ab_seg_[iCh][0][iseg];
1078 delta2[1] = vec.at(itr).param2[1] - par_ab_seg_[iCh][1][iseg];
1079 delta2[2] = vec.at(itr).param2[2] - par_ab_seg_[iCh][2][iseg];
1080 delta2[3] = vec.at(itr).param2[3] - par_ab_seg_[iCh][3][iseg];
1081
1082 if (fDebug) {
1083 // combinatorics
1084 if (delta2[0] > -90.)
1085 hdAx_mc_SegCh.at(iCh)->Fill(delta2[0]);
1086 if (delta2[1] > -900.)
1087 hdX_mc_SegCh.at(iCh)->Fill(delta2[1]);
1088 if (delta2[2] > -90.)
1089 hdAy_mc_SegCh.at(iCh)->Fill(delta2[2]);
1090 if (delta2[3] > -900.)
1091 hdY_mc_SegCh.at(iCh)->Fill(delta2[3]);
1092 }
1093
1094 dmatch = 0.;
1095 dmatch = (delta2[0] * delta2[0]) / (sig[0] * sig[0]) + (delta2[1] * delta2[1]) / (sig[1] * sig[1])
1096 + (delta2[2] * delta2[2]) / (sig[2] * sig[2])
1097 + (delta2[3] * delta2[3]) / (sig[3] * sig[3]);
1098
1099 if (dmc_match[iCh][itr] > dmatch) {
1100 dmc_match[iCh][itr] = dmatch;
1101 mc_tr_assoc[iCh][itr] = iseg;
1102 dax = delta2[0];
1103 dx = delta2[1];
1104 day = delta2[2];
1105 dy = delta2[3];
1106 }
1107 }
1108 if (fDebug) {
1109 cout << " itr " << itr << " Np2 " << vec.at(itr).Np2 << " mc_Id " << vec.at(itr).Id << " bx_mc "
1110 << vec.at(itr).param2[1] << " iseg_ind " << mc_tr_assoc[iCh][itr] << " bx "
1111 << par_ab_seg_[iCh][1][mc_tr_assoc[iCh][itr]] << " dmc_match " << dmc_match[iCh][itr] << endl;
1112 if (mc_tr_assoc[iCh][itr] > -1) {
1113 if (dax > -900.)
1114 hdAx_mc_Seg_deltaCh.at(iCh)->Fill(dax);
1115 if (dx > -900.)
1116 hdX_mc_Seg_deltaCh.at(iCh)->Fill(dx);
1117 if (day > -900.)
1118 hdAy_mc_Seg_deltaCh.at(iCh)->Fill(day);
1119 if (dy > -900.)
1120 hdY_mc_Seg_deltaCh.at(iCh)->Fill(dy);
1121 }
1122 }
1123 } // iChamber == 2
1124
1125 if (iCh == 3) {
1126 for (Int_t iseg = 0; iseg < Nbest_seg_[iCh]; iseg++) {
1127 delta3[0] = vec.at(itr).param3[0] - par_ab_seg_[iCh][0][iseg];
1128 delta3[1] = vec.at(itr).param3[1] - par_ab_seg_[iCh][1][iseg];
1129 delta3[2] = vec.at(itr).param3[2] - par_ab_seg_[iCh][2][iseg];
1130 delta3[3] = vec.at(itr).param3[3] - par_ab_seg_[iCh][3][iseg];
1131
1132 if (fDebug) {
1133 // combinatorics
1134 if (delta3[0] > -90.)
1135 hdAx_mc_SegCh.at(iCh)->Fill(delta3[0]);
1136 if (delta3[1] > -900.)
1137 hdX_mc_SegCh.at(iCh)->Fill(delta3[1]);
1138 if (delta3[2] > -90.)
1139 hdAy_mc_SegCh.at(iCh)->Fill(delta3[2]);
1140 if (delta3[3] > -900.)
1141 hdY_mc_SegCh.at(iCh)->Fill(delta3[3]);
1142 }
1143
1144 dmatch = 0.;
1145 dmatch = (delta3[0] * delta3[0]) / (sig[0] * sig[0]) + (delta3[1] * delta3[1]) / (sig[1] * sig[1])
1146 + (delta3[2] * delta3[2]) / (sig[2] * sig[2])
1147 + (delta3[3] * delta3[3]) / (sig[3] * sig[3]);
1148
1149 if (dmc_match[iCh][itr] > dmatch) {
1150 dmc_match[iCh][itr] = dmatch;
1151 mc_tr_assoc[iCh][itr] = iseg;
1152 dax = delta3[0];
1153 dx = delta3[1];
1154 day = delta3[2];
1155 dy = delta3[3];
1156 }
1157 }
1158
1159 if (fDebug) {
1160 cout << " itr " << itr << " Np3 " << vec.at(itr).Np3 << " mc_Id " << vec.at(itr).Id << " bx_mc "
1161 << vec.at(itr).param3[1] << " iseg_ind " << mc_tr_assoc[iCh][itr] << " bx "
1162 << par_ab_seg_[iCh][1][mc_tr_assoc[iCh][itr]] << " dmc_match " << dmc_match[iCh][itr] << endl;
1163
1164 if (mc_tr_assoc[iCh][itr] > -1) {
1165
1166 if (dax > -900.)
1167 hdAx_mc_Seg_deltaCh.at(iCh)->Fill(dax);
1168 if (dx > -900.)
1169 hdX_mc_Seg_deltaCh.at(iCh)->Fill(dx);
1170 if (day > -900.)
1171 hdAy_mc_Seg_deltaCh.at(iCh)->Fill(day);
1172 if (dy > -900.)
1173 hdY_mc_Seg_deltaCh.at(iCh)->Fill(dy);
1174 }
1175 }
1176 } // iChamber == 3
1177 //}//iChamber
1178
1179 } // vec_points.size
1180
1181 if (fDebug)
1182 cout << "reject poorly chosen association segments " << endl;
1183
1184 // for (Int_t iChamber = 2; iChamber < kNChambers; iChamber++) {
1185 for (size_t itr = 0; itr < vec.size(); itr++) { // mc_tr
1186 if (mc_tr_assoc[iCh][itr] == -1)
1187 continue;
1188
1189 for (size_t itr2 = 0; itr2 < vec.size(); itr2++) { // mc_tr
1190 if (itr2 == itr)
1191 continue;
1192 if (mc_tr_assoc[iCh][itr2] == -1)
1193 continue;
1194
1195 if (mc_tr_assoc[iCh][itr] == mc_tr_assoc[iCh][itr2]) {
1196 if (dmc_match[iCh][itr2] > 20.)
1197 mc_tr_assoc[iCh][itr2] = -1;
1198 else {
1199 // mc_tr_assoc[iCh][itr] = -1;
1200 continue;
1201 }
1202 }
1203 } // itr2
1204 //---MC Eff ---
1205 //---Num
1206 if (fDebug)
1207 cout << " mc_Id " << vec.at(itr).Id << " assoc " << mc_tr_assoc[iCh][itr] << " iChamber " << iCh
1208 << endl;
1209 if (fDebug && mc_tr_assoc[iCh][itr] > -1) {
1210 if (fDebug && iCh == 2 && vec.at(itr).Np2 >= 4) {
1211 if (vec.at(itr).xWas2 && vec.at(itr).uWas2 && vec.at(itr).vWas2) {
1212 hNum_mc->Fill(2);
1213 hNum_mc->Fill(4);
1214 }
1215 hNum_mc->Fill(0);
1216 cout << " Num_mcPC2 " << endl;
1217
1218 Nassoc2++;
1219 }
1220 if (fDebug && iCh == 3 && vec.at(itr).Np3 >= 4) {
1221 if (vec.at(itr).xWas3 && vec.at(itr).uWas3 && vec.at(itr).vWas3) {
1222 hNum_mc->Fill(3);
1223 hNum_mc->Fill(5);
1224 }
1225 hNum_mc->Fill(1);
1226 cout << " Num_mcPC3 " << endl;
1227
1228 Nassoc3++;
1229 }
1230 }
1231 } // itr
1232 // }//iChamber
1233 if (Nassoc2 == 2 && wasLi7_ch2 && wasHe4_ch2) {
1234 if (fDebug)
1235 cout << " hNum_mcreaction " << endl;
1236 hNum_mcreaction2->Fill(0);
1237 }
1238 if (Nassoc3 == 2 && wasLi7_ch3 && wasHe4_ch3) {
1239 if (fDebug)
1240 cout << " hNum_mcreaction " << endl;
1241 hNum_mcreaction3->Fill(0);
1242 }
1243
1244 } // if(!expData)
1245}
1246
1247void BmnMwpcHitFinder::ReadWires(Double_t*** DigitsArray_, Int_t** iw_Ch_, vector<MC_points>& vec)
1248{
1249 if (fDebug)
1250 cout << "--ReadWires--" << endl;
1251
1252 Int_t Fired_layer_[kNChambers][kNPlanes];
1253 int Npl_MC2[kMaxMC];
1254 int Npl_MC3[kMaxMC];
1255 Short_t st, wire, pn, pl;
1256 UInt_t ts;
1257
1258 for (Int_t iChamber = 0; iChamber < kNChambers; iChamber++) {
1259 for (Int_t ipll = 0; ipll < kNPlanes; ipll++) {
1260 Fired_layer_[iChamber][ipll] = 0;
1261 }
1262 }
1263
1264 if (!expData) {
1265
1266 Int_t mcTracksArray[kMaxMC];
1267 Int_t mcPdgCode[kMaxMC];
1268 Double_t X2mc[kMaxMC][kNPlanes];
1269 Double_t Y2mc[kMaxMC][kNPlanes];
1270 Double_t Z2mc[kMaxMC][kNPlanes];
1271 Double_t X3mc[kMaxMC][kNPlanes];
1272 Double_t Y3mc[kMaxMC][kNPlanes];
1273 Double_t Z3mc[kMaxMC][kNPlanes];
1274 for (Int_t Id = 0; Id < kMaxMC; Id++) {
1275 Npl_MC2[Id] = 0;
1276 Npl_MC3[Id] = 0;
1277 mcTracksArray[Id] = -1;
1278 mcPdgCode[Id] = -1;
1279 for (Int_t i = 0; i < kNPlanes; i++) {
1280 X2mc[Id][i] = -999.;
1281 Y2mc[Id][i] = -999.;
1282 Z2mc[Id][i] = -999.;
1283 X3mc[Id][i] = -999.;
1284 Y3mc[Id][i] = -999.;
1285 Z3mc[Id][i] = -999.;
1286 }
1287 }
1288
1289 int tr_before = -1;
1290 int Nmc_tracks = -1;
1291
1292 for (Int_t iMC = 0; iMC < fBmnHitsArray->GetEntriesFast(); ++iMC) {
1293 BmnMwpcHit* hit = (BmnMwpcHit*)fBmnHitsArray->UncheckedAt(iMC);
1294
1295 Int_t st_MC = hit->GetMwpcId();
1296 Int_t trackId_MC = hit->GetHitId();
1297 Int_t pdg_MC = hit->GetPdgId();
1298 Int_t pl_MC = hit->GetPlaneId();
1299 Short_t wire_MC = hit->GetWireNumber();
1300 Double_t time_MC = hit->GetWireTime();
1301
1302 if (fDebug)
1303 cout << " st_MC " << st_MC << " trackId_MC " << trackId_MC << " pdg_MC " << pdg_MC << " pl_MC " << pl_MC
1304 << " X " << hit->GetX() << " wire_MC " << wire_MC << endl;
1305
1306 if (tr_before != trackId_MC) {
1307 Nmc_tracks++;
1308 mcTracksArray[Nmc_tracks] = hit->GetHitId();
1309 }
1310 tr_before = trackId_MC;
1311 pn = kPln[st_MC][pl_MC]; // made for the canonical sequence / x- v- u+ x+ v+ u-/
1312 // if (fDebug)cout<<" pn "<<pn<<" pl_MC "<<pl_MC<<endl;
1313
1314 if (st_MC == 2) {
1315 X2mc[Nmc_tracks][pn] = hit->GetX();
1316 Y2mc[Nmc_tracks][pn] = hit->GetY();
1317 Z2mc[Nmc_tracks][pn] = hit->GetZ();
1318 Npl_MC2[Nmc_tracks]++;
1319 }
1320 if (st_MC == 3) {
1321 X3mc[Nmc_tracks][pn] = hit->GetX();
1322 Y3mc[Nmc_tracks][pn] = hit->GetY();
1323 Z3mc[Nmc_tracks][pn] = hit->GetZ();
1324 mcPdgCode[Nmc_tracks] = pdg_MC;
1325 Npl_MC3[Nmc_tracks]++;
1326 }
1327 // if (fDebug)cout<<" X2["<<Nmc_tracks<<"]["<<pn<<"] "<<X2mc[Nmc_tracks][pn]<<" Npl_MC2 "<<
1328 // Npl_MC2[Nmc_tracks]<<" X3["<<Nmc_tracks<<"]["<<pn<<"] "<<X3mc[Nmc_tracks][pn]<<" Npl_MC3 "<<
1329 // Npl_MC3[Nmc_tracks]<<endl;
1330
1331 DigitsArray_[st_MC][pn][wire_MC] = time_MC;
1332 Fired_layer_[st_MC][pn] = 1;
1333 iw_Ch_[st_MC][pn]++;
1334 } // iMC
1335 Nmc_tracks++;
1336
1337 MC_points tmpTr;
1338 if (fDebug)
1339 cout << " Nmc_tracks " << Nmc_tracks << endl;
1340 for (Int_t id = 0; id < Nmc_tracks; id++) {
1341 if (fDebug)
1342 cout << " Id_mc " << mcTracksArray[id] << " Npl2 " << Npl_MC2[id] << " Npl3 " << Npl_MC3[id] << endl;
1343 for (Int_t i = 0; i < 6; i++) {
1344 if (fDebug && X2mc[id][i] > -900.)
1345 cout << " ipl " << i << " X2 " << X2mc[id][i] << " Y " << Y2mc[id][i] << endl;
1346 tmpTr.x2[i] = X2mc[id][i];
1347 tmpTr.y2[i] = Y2mc[id][i];
1348 tmpTr.z2[i] = Z2mc[id][i];
1349 }
1350 if (fDebug)
1351 cout << "---" << endl;
1352 for (Int_t i = 0; i < 6; i++) {
1353 if (fDebug && X3mc[id][i] > -900.)
1354 cout << " ipl " << i << " X3 " << X3mc[id][i] << " Y " << Y3mc[id][i] << endl;
1355 tmpTr.x3[i] = X3mc[id][i];
1356 tmpTr.y3[i] = Y3mc[id][i];
1357 tmpTr.z3[i] = Z3mc[id][i];
1358 }
1359 if (fDebug)
1360 cout << "---" << endl;
1361 if (fDebug)
1362 cout << endl;
1363 tmpTr.Id = mcTracksArray[id];
1364 tmpTr.Np2 = Npl_MC2[id];
1365 tmpTr.Np3 = Npl_MC3[id];
1366 tmpTr.Pdg = mcPdgCode[id];
1367
1368 if (Npl_MC2[id] >= 4 || Npl_MC3[id] >= 4)
1369 vec.push_back(tmpTr);
1370 } // Nmc_tracks
1371
1372 if (fDebug)
1373 cout << " MC vec_points.size() " << vec.size() << endl;
1374 if (fDebug)
1375 hNtrMC->Fill(vec.size());
1376
1377 Double_t x_target_ch2, y_target_ch2, x_target_ch3, y_target_ch3;
1378
1379 for (size_t itr = 0; itr < vec.size(); itr++) {
1380
1381 if (vec.at(itr).x2[0] > -900. || vec.at(itr).x2[3] > -900.)
1382 vec.at(itr).xWas2 = 1;
1383 if (vec.at(itr).x2[1] > -900. || vec.at(itr).x2[4] > -900.)
1384 vec.at(itr).vWas2 = 1;
1385 if (vec.at(itr).x2[2] > -900. || vec.at(itr).x2[5] > -900.)
1386 vec.at(itr).uWas2 = 1;
1387
1388 if (vec.at(itr).x3[0] > -900. || vec.at(itr).x3[3] > -900.)
1389 vec.at(itr).xWas3 = 1;
1390 if (vec.at(itr).x3[1] > -900. || vec.at(itr).x3[4] > -900.)
1391 vec.at(itr).vWas3 = 1;
1392 if (vec.at(itr).x3[2] > -900. || vec.at(itr).x3[5] > -900.)
1393 vec.at(itr).uWas3 = 1;
1394
1395 int i2 = -1;
1396 if (vec.at(itr).x2[5] > -900.)
1397 i2 = 5;
1398 if (i2 < 0 && vec.at(itr).x2[4] > -900.)
1399 i2 = 4;
1400 if (i2 < 0 && vec.at(itr).x2[3] > -900.)
1401 i2 = 3;
1402
1403 int i20 = 0;
1404 if (vec.at(itr).x2[i20] < -900.)
1405 i20 = 1;
1406 if (vec.at(itr).x2[i20] < -900.)
1407 i20 = 2;
1408
1409 if (i2 > -1 && vec.at(itr).Np2 > 3) {
1410
1411 vec.at(itr).param2[0] =
1412 (vec.at(itr).x2[i20] - vec.at(itr).x2[i2]) / (vec.at(itr).z2[i20] - vec.at(itr).z2[i2]);
1413 vec.at(itr).param2[1] = (vec.at(itr).x2[i20] + vec.at(itr).x2[i2]) * 0.5;
1414 vec.at(itr).param2[2] =
1415 (vec.at(itr).y2[i20] - vec.at(itr).y2[i2]) / (vec.at(itr).z2[i20] - vec.at(itr).z2[i2]);
1416 vec.at(itr).param2[3] = (vec.at(itr).y2[i20] + vec.at(itr).y2[i2]) * 0.5;
1417 x_target_ch2 = vec.at(itr).param2[0] * (Z0_SRC - ChZ[2]) + vec.at(itr).param2[1];
1418 y_target_ch2 = vec.at(itr).param2[2] * (Z0_SRC - ChZ[2]) + vec.at(itr).param2[3];
1419 }
1420 int i3 = -1;
1421 if (vec.at(itr).x3[5] > -900.)
1422 i3 = 5;
1423 if (i3 < 0 && vec.at(itr).x3[4] > -900.)
1424 i3 = 4;
1425 if (i3 < 0 && vec.at(itr).x3[3] > -900.)
1426 i3 = 3;
1427
1428 int i30 = 0;
1429 if (vec.at(itr).x3[i30] < -900.)
1430 i30 = 1;
1431 if (vec.at(itr).x3[i30] < -900.)
1432 i30 = 2;
1433
1434 if (i3 > -1 && vec.at(itr).Np3 > 3) {
1435
1436 vec.at(itr).param3[0] =
1437 (vec.at(itr).x3[i30] - vec.at(itr).x3[i3]) / (vec.at(itr).z3[i30] - vec.at(itr).z3[i3]);
1438 vec.at(itr).param3[1] = (vec.at(itr).x3[i30] + vec.at(itr).x3[i3]) * 0.5;
1439 vec.at(itr).param3[2] =
1440 (vec.at(itr).y3[i30] - vec.at(itr).y3[i3]) / (vec.at(itr).z3[i30] - vec.at(itr).z3[i3]);
1441 vec.at(itr).param3[3] = (vec.at(itr).y3[i30] + vec.at(itr).y3[i3]) * 0.5;
1442 x_target_ch3 = vec.at(itr).param3[0] * (Z0_SRC - ChZ[3]) + vec.at(itr).param3[1];
1443 y_target_ch3 = vec.at(itr).param3[2] * (Z0_SRC - ChZ[3]) + vec.at(itr).param3[3];
1444
1445 vec.at(itr).xt = (x_target_ch2 + x_target_ch3) / 2;
1446 vec.at(itr).yt = (y_target_ch2 + y_target_ch3) / 2;
1447 }
1448 if (fDebug)
1449 cout << " itr " << itr << " Id_mc " << vec.at(itr).Id << " 2:Ax " << vec.at(itr).param2[0] << " bx "
1450 << vec.at(itr).param2[1] << " Ay " << vec.at(itr).param2[2] << " by " << vec.at(itr).param2[3]
1451 << " x_target " << x_target_ch2 << " y_target " << y_target_ch2 << " Np2 " << vec.at(itr).Np2
1452 << endl;
1453 // tmp
1454 // Double_t v = (x_target_ch2 - sq3*y_target_ch2)*0.5;
1455 // Double_t u = (x_target_ch2 + sq3*y_target_ch2)*0.5;
1456
1457 // if (fDebug) cout<<" vt_xv: "<< (2*x_target_ch2 - v)/sq3 <<" ut_xu: "<< (2*x_target_ch2 - u)/sq3 <<endl;
1458 // if (fDebug) cout<<" vt_uv: "<< (2*u + v)/sq3 <<" ut_uv: "<< (u + 2*v)/sq3 <<endl;
1459
1460 if (fDebug)
1461 cout << " itr " << itr << " Id_mc " << vec.at(itr).Id << " 3:Ax " << vec.at(itr).param3[0] << " bx "
1462 << vec.at(itr).param3[1] << " Ay " << vec.at(itr).param3[2] << " by " << vec.at(itr).param3[3]
1463 << " x_target " << x_target_ch3 << " y_target " << y_target_ch3 << " Np3 " << vec.at(itr).Np3
1464 << endl;
1465 if (fDebug)
1466 cout << " xt " << vec.at(itr).xt << " yt " << vec.at(itr).yt << " Pdg " << vec.at(itr).Pdg << endl;
1467
1468 if (fDebug) {
1469 if (vec.at(itr).Np2 > 2) {
1470 hAx_mc_ch.at(2)->Fill(vec.at(itr).param2[0]);
1471 hBx_mc_ch.at(2)->Fill(vec.at(itr).param2[1]);
1472 hAy_mc_ch.at(2)->Fill(vec.at(itr).param2[2]);
1473 hBy_mc_ch.at(2)->Fill(vec.at(itr).param2[3]);
1474 hYvsX_mc_ch2->Fill(vec.at(itr).param2[1], vec.at(itr).param2[3]);
1475 }
1476 if (vec.at(itr).Np3 > 2) {
1477 hAx_mc_ch.at(3)->Fill(vec.at(itr).param3[0]);
1478 hBx_mc_ch.at(3)->Fill(vec.at(itr).param3[1]);
1479 hAy_mc_ch.at(3)->Fill(vec.at(itr).param3[2]);
1480 hBy_mc_ch.at(3)->Fill(vec.at(itr).param3[3]);
1481 hYvsX_mc_ch3->Fill(vec.at(itr).param3[1], vec.at(itr).param3[3]);
1482 }
1483 }
1484
1485 } // vec_points
1486 if (fDebug)
1487 cout << endl;
1488
1489 } else {
1490 for (Int_t iDigit = 0; iDigit < fBmnMwpcDigitArray->GetEntries(); iDigit++) {
1491 BmnMwpcDigit* digit = (BmnMwpcDigit*)fBmnMwpcDigitArray->At(iDigit);
1492
1493 st = digit->GetStation();
1494 wire = digit->GetWireNumber();
1495 pl = digit->GetPlane();
1496 ts = digit->GetTime(); // ns
1497
1498 if (fRunPeriod == 7 && fRunNumber > 3588) { // BMN
1499 if (st > 1)
1500 st = st - 2;
1501 }
1502
1503 Fired_layer_[st][pl] = 1;
1504 Int_t ind = st * kNPlanes + pl;
1505 if (fDebug)
1506 hTime.at(ind)->Fill(ts);
1507 if (fDebug)
1508 hOccupancy.at(ind)->Fill(wire);
1509 if (ts < 80 || ts > 280)
1510 continue;
1511
1512 pn = kPln[st][pl]; // made for the canonical sequence / x- v- u+ x+ v+ u-/
1513
1514 // Loop over repeated wires
1515 Bool_t repeat_wire = 0;
1516 if (iw_Ch_[st][pn] > 0) {
1517 for (Int_t ix = 0; ix < iw_Ch[st][pn]; ix++) {
1518 if (wire == wire_Ch[st][ix][pn]) {
1519 repeat_wire = 1;
1520 break;
1521 }
1522 } // ix
1523 }
1524 wire_Ch[st][iw_Ch_[st][pn]][pn] = wire;
1525
1526 if (repeat_wire)
1527 continue;
1528 if (iw_Ch_[st][pn] >= 80)
1529 continue;
1530
1531 DigitsArray_[st][pn][wire] = ts;
1532 iw_Ch_[st][pn]++;
1533 } // iDigit
1534 } // else
1535 for (Int_t iChamber = 0; iChamber < kNChambers; iChamber++) {
1536 for (Int_t iplane = 0; iplane < kNPlanes; iplane++) {
1537 Nlay_w_wires[iChamber] += Fired_layer_[iChamber][iplane];
1538 counter_pl[iplane] = 0;
1539 counter_pl[iplane] += iw_Ch_[iChamber][iplane];
1540 if (fDebug)
1541 hfiredWire_Ch.at(iChamber)->Fill(counter_pl[iplane]);
1542 } // iplane
1543 if (fDebug) {
1544 hNFired_layers_Ch[iChamber]->Fill(Nlay_w_wires[iChamber]);
1545 if (fDebug && iChamber < 2 && Nlay_w_wires[iChamber] < kMinHits)
1546 hHalfDeadCh->Fill(iChamber);
1547 if (fDebug && iChamber >= 2 && Nlay_w_wires[0] > kMinHits && Nlay_w_wires[1] > kMinHits
1548 && Nlay_w_wires[iChamber] < kMinHits)
1549 hHalfDeadCh->Fill(iChamber);
1550 }
1551 } // ch
1552
1553 if (fDebug && Nlay_w_wires[0] > kMinHits && Nlay_w_wires[1] > kMinHits && Nlay_w_wires[2] < kMinHits
1554 && Nlay_w_wires[3] < kMinHits)
1555 hHalfDeadCh->Fill(4);
1556}
1557//--------------------------------------------------------------------
1558
1559//------------------ Cluster production function------------------------
1560void BmnMwpcHitFinder::Clustering(Int_t chNum,
1561 Int_t*** ClusterSize_,
1562 Double_t*** DigitsArray_,
1563 Double_t*** Coord_wire_,
1564 Double_t*** Coord_xuv_,
1565 Int_t** Nclust_,
1566 Int_t* counter_pl_)
1567{
1568 Int_t Nfirst[kCh_max][kNPlanes] /*, Nlast[kCh_max][kNPlanes]*/;
1569 // Int_t Min_time_wires, fast_wire ;
1570 // Int_t N_wires[kCh_max][kNPlanes];
1571 // Int_t N_1_cluster[kCh_max][kNPlanes];
1572 // Int_t Fired_layer_[kNChambers][kNPlanes];
1573 // Double_t earlyWires[kNChambers][kNPlanes][kNWires];
1574 // Double_t Coord_fast[kNChambers][kNPlanes][kBig];
1575 // Double_t Cut_time_wire = 16.;//ns
1576 // Double_t Num_layers_out_beam = 0;
1577 Bool_t wire_was = 0;
1578
1579 /*
1580 for (Int_t ipll = 0; ipll < kNPlanes; ipll++) {
1581 Fired_layer_[chNum][ipll] =0;
1582 N_wires[chNum][ipll] = 0;
1583 N_1_cluster[chNum][ipll] = 0;
1584 for (Int_t iwire = 0; iwire < kNWires; iwire++) {
1585 if ( DigitsArray_[chNum][ipll][iwire] > 0 ) Fired_layer_[chNum][ipll] = 1;
1586 }
1587 for (Int_t ib = 0; ib < kBig; ib++) {
1588 Coord_fast[chNum][ipll][ib] = 0.;
1589 }
1590 }*/
1591
1592 /*
1593 if ( chNum < 2 ) {
1594 for (Int_t ipll = 0; ipll < kNPlanes; ipll++) {
1595 Nfirst[chNum][ipll] = -1;
1596 Nlast[chNum][ipll] = -1;
1597 Min_time_wires = 999.;
1598 fast_wire = -1;
1599
1600 for (Int_t iwire = 0; iwire < kNWires; iwire++) {
1601 earlyWires[chNum][ipll][iwire] = 0.;
1602 if ( DigitsArray_[chNum][ipll][iwire] > 0 && DigitsArray_[chNum][ipll][iwire] < Min_time_wires ) {
1603 Min_time_wires = DigitsArray_[chNum][ipll][iwire];
1604 }
1605 }
1606 if (fDebug) cout<<" --ipll "<<ipll<<" Min_time_wires "<<Min_time_wires<<endl;
1607
1608 for (Int_t iwire = 0; iwire < kNWires; iwire++) {
1609 if (fDebug && DigitsArray_[chNum][ipll][iwire] > 0) cout<<" DigitsArray_["<<chNum<<"]["<<ipll<<"]["<<iwire<<"]
1610 "<<DigitsArray_[chNum][ipll][iwire]<<endl;
1611
1612 if ( fDebug && DigitsArray_[chNum][ipll][iwire] > 0 ){
1613 htime_wire_Ch.at(chNum)->Fill(iwire + kNWires*ipll, DigitsArray_[chNum][ipll][iwire]);
1614 hfired_wire_Ch.at(chNum)->Fill(iwire + kNWires*ipll);
1615
1616 if(ipll == 0) hWiresXm.at(chNum)->Fill(iwire);
1617 if(ipll == 1) hWiresVm.at(chNum)->Fill(iwire);
1618 if(ipll == 2) hWiresUp.at(chNum)->Fill(iwire);
1619 if(ipll == 3) hWiresXp.at(chNum)->Fill(iwire);
1620 if(ipll == 4) hWiresVp.at(chNum)->Fill(iwire);
1621 if(ipll == 5) hWiresUm.at(chNum)->Fill(iwire);
1622
1623 if ( ipll == 0 && fEventNo < 700) hEvent_display_Ch.at(chNum)->Fill(fEventNo,iwire);
1624 }
1625
1626 if( Nfirst[chNum][ipll] < 0 && DigitsArray_[chNum][ipll][iwire] == 0.) continue;
1627 if( Nfirst[chNum][ipll] < 0 && DigitsArray_[chNum][ipll][iwire] > (Min_time_wires + Cut_time_wire) ) continue;
1628
1629 if( Nfirst[chNum][ipll] < 0 && DigitsArray_[chNum][ipll][iwire] > 0.
1630 && DigitsArray_[chNum][ipll][iwire] <= (Min_time_wires + Cut_time_wire)) Nfirst[chNum][ipll] = iwire;
1631
1632 if( Nfirst[chNum][ipll] >=0 && (DigitsArray_[chNum][ipll][iwire+1] == 0.
1633 || DigitsArray_[chNum][ipll][iwire+ 1] > (Min_time_wires + Cut_time_wire) ) )
1634 Nlast[chNum][ipll] = iwire;
1635
1636 if (Nfirst[chNum][ipll] >= 0 && Nlast[chNum][ipll] >= Nfirst[chNum][ipll]) {
1637 Int_t min_time_in_clust = 999, num_fast_wire = -1;
1638
1639 for (Int_t iww = Nfirst[chNum][ipll] ; iww < Nlast[chNum][ipll] + 1; iww++) {
1640
1641 if ( DigitsArray_[chNum][ipll][iww] < min_time_in_clust ) {
1642 min_time_in_clust = DigitsArray_[chNum][ipll][iww] ;
1643 num_fast_wire = iww;
1644 }
1645 }
1646
1647 ClusterSize_[chNum][ipll][Nclust_[chNum][ipll]] = Nlast[chNum][ipll] - Nfirst[chNum][ipll] + 1;
1648
1649
1650 if (fDebug) hClusterSize.at(chNum)->Fill(ClusterSize_[chNum][ipll][Nclust_[chNum][ipll]]);
1651
1652 if ( min_time_in_clust < 900) {
1653 Coord_wire_[chNum][ipll][Nclust_[chNum][ipll] ] = Nfirst[chNum][ipll] +
1654 0.5*ClusterSize_[chNum][ipll][Nclust_[chNum][ipll]]; Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]] =
1655 (Coord_wire_[chNum][ipll][Nclust_[chNum][ipll]] - kMiddlePl)* dw;//cm
1656
1657 // made for the canonical sequence / x- v- u+ x+ v+ u-/ 0 1 5 have back reading
1658 if (ipll == 0 || ipll == 1 || ipll == 5 ) Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]] = -
1659 Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]];
1660 }
1661 if ( fDebug) cout<<" Nfirst "<<Nfirst[chNum][ipll]<<" Nlast "<<Nlast[chNum][ipll]<<"
1662 ClusterSize_["<<chNum<<"]["<<ipll<<"] "<<ClusterSize_[chNum][ipll][Nclust_[chNum][ipll]]<<" Coord_xuv
1663 "<<Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]]<<endl;
1664
1665 if (fDebug) {
1666 hoccupancyXm.at(chNum)->Fill(Coord_xuv_[chNum][0][Nclust_[chNum][ipll]]);
1667 hoccupancyVm.at(chNum)->Fill(Coord_xuv_[chNum][1][Nclust_[chNum][ipll]]);
1668 hoccupancyUp.at(chNum)->Fill(Coord_xuv_[chNum][2][Nclust_[chNum][ipll]]);
1669 hoccupancyXp.at(chNum)->Fill(Coord_xuv_[chNum][3][Nclust_[chNum][ipll]]);
1670 hoccupancyVp.at(chNum)->Fill(Coord_xuv_[chNum][4][Nclust_[chNum][ipll]]);
1671 hoccupancyUm.at(chNum)->Fill(Coord_xuv_[chNum][5][Nclust_[chNum][ipll]]);
1672 }
1673
1674 Nclust_[chNum][ipll]++;
1675
1676 Nfirst[chNum][ipll] = -1;
1677 Nlast[chNum][ipll] = -1;
1678
1679 } //if (Nfirst[chNum][ipll] >= 0
1680 }//iwire
1681 }// ipll
1682 } else {
1683 */
1684
1685 // Num_layers_out_beam = 0;
1686
1687 for (Int_t ipll = 0; ipll < kNPlanes; ipll++) {
1688
1689 // if (fDebug) cout<<" counter_pl["<<ipll<<"] "<<counter_pl_[ipll]<<endl;
1690
1691 if (counter_pl_[ipll] <= 8) { // for MC-case
1692 if (fDebug)
1693 cout << "<= 8 " << endl;
1694 for (Int_t iwire = 0; iwire < kNWires; iwire++) {
1695
1696 if (DigitsArray_[chNum][ipll][iwire] > 0.) {
1697 Nfirst[chNum][ipll] = iwire;
1698 ClusterSize_[chNum][ipll][Nclust_[chNum][ipll]] = 1;
1699 Coord_wire_[chNum][ipll][Nclust_[chNum][ipll]] =
1700 Nfirst[chNum][ipll] + 0.5 * ClusterSize_[chNum][ipll][Nclust_[chNum][ipll]];
1701 Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]] =
1702 (Coord_wire_[chNum][ipll][Nclust_[chNum][ipll]] - kMiddlePl) * dw;
1703 if (ipll == 0 || ipll == 1 || ipll == 5)
1704 Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]] = -Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]];
1705 if (fDebug)
1706 cout << " Coord_xuv_[" << chNum << "][" << ipll << "][" << Nclust_[chNum][ipll] << "] "
1707 << Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]] << " Coord_wire_ "
1708 << Coord_wire_[chNum][ipll][Nclust_[chNum][ipll]] << endl;
1709 Nclust_[chNum][ipll]++;
1710 }
1711 } // iwire
1712
1713 } else { // loaded area
1714
1715 Nfirst[chNum][ipll] = -1;
1716 // Nlast[chNum][ipll] = -1;
1717 Bool_t Next_next_wire = 0;
1718 // Int_t Last_wire_Digit = 1000;
1719
1720 for (Int_t iwire = 0; iwire < kNWires; iwire++) {
1721 wire_was = 0;
1722
1723 if ((ipll == 0 || ipll == 3) && (iwire < 33 || iwire > 63)) {
1724 // if ( DigitsArray_[chNum][ipll][iwire] > 0 ) Fired_layer_[chNum][ipll] = 1;
1725 }
1726
1727 if ((ipll == 1 || ipll == 2) && (iwire < 12 || iwire > 42)) {
1728 // if ( DigitsArray_[chNum][ipll][iwire] > 0 ) Fired_layer_[chNum][ipll] = 1;
1729 }
1730 if ((ipll == 4 || ipll == 5) && (iwire < 50 || iwire > 80)) {
1731 // if ( DigitsArray_[chNum][ipll][iwire] > 0 ) Fired_layer_[chNum][ipll] = 1;
1732 }
1733
1734 if (fDebug && DigitsArray_[chNum][ipll][iwire] > 0) {
1735 htime_wire_Ch.at(chNum)->Fill(iwire + kNWires * ipll, DigitsArray_[chNum][ipll][iwire]);
1736 hfired_wire_Ch.at(chNum)->Fill(iwire + kNWires * ipll);
1737
1738 if (ipll == 0)
1739 hWiresXm.at(chNum)->Fill(iwire);
1740 if (ipll == 1)
1741 hWiresVm.at(chNum)->Fill(iwire);
1742 if (ipll == 2)
1743 hWiresUp.at(chNum)->Fill(iwire);
1744 if (ipll == 3)
1745 hWiresXp.at(chNum)->Fill(iwire);
1746 if (ipll == 4)
1747 hWiresVp.at(chNum)->Fill(iwire);
1748 if (ipll == 5)
1749 hWiresUm.at(chNum)->Fill(iwire);
1750
1751 if (ipll == 0 && fEventNo < 700)
1752 hEvent_display_Ch.at(chNum)->Fill(fEventNo, iwire);
1753 } //
1754
1755 // if (fDebug ) cout<<" 1 first "<<Nfirst[chNum][ipll]<<" time
1756 // "<<DigitsArray_[chNum][ipll][iwire]<<endl;
1757 if (fDebug && DigitsArray_[chNum][ipll][iwire] > 0)
1758 cout << " DigitsArray_[" << chNum << "][" << ipll << "][" << iwire << "] "
1759 << DigitsArray_[chNum][ipll][iwire] << " first " << Nfirst[chNum][ipll] << endl;
1760
1761 // area without beam for ch2 & ch3
1762 // if (chNum > 1 ){
1763 if (chNum > 1 && DigitsArray_[chNum][ipll][iwire] > 0
1764 && (iwire < Beam_wires_min[chNum][ipll] || iwire > Beam_wires_max[chNum][ipll]))
1765 {
1766 wire_was = 1;
1767 if (fDebug)
1768 cout << " area without beam wire: " << iwire << " wire_was " << wire_was << endl;
1769
1770 Nfirst[chNum][ipll] = iwire;
1771 ClusterSize_[chNum][ipll][Nclust_[chNum][ipll]] = 1;
1772 Coord_wire_[chNum][ipll][Nclust_[chNum][ipll]] =
1773 Nfirst[chNum][ipll] + 0.5 * ClusterSize_[chNum][ipll][Nclust_[chNum][ipll]];
1774 Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]] =
1775 (Coord_wire_[chNum][ipll][Nclust_[chNum][ipll]] - kMiddlePl) * dw;
1776 if (ipll == 0 || ipll == 1 || ipll == 5)
1777 Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]] = -Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]];
1778 if (fDebug)
1779 cout << " Coord_xuv_[" << chNum << "][" << ipll << "][" << Nclust_[chNum][ipll] << "] "
1780 << Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]] << " Coord_wire_ "
1781 << Coord_wire_[chNum][ipll][Nclust_[chNum][ipll]] << endl;
1782 Nclust_[chNum][ipll]++;
1783
1784 } else { // area with beam
1785
1786 if (wire_was)
1787 continue;
1788
1789 if (Nfirst[chNum][ipll] < 0 && DigitsArray_[chNum][ipll][iwire] == 0) {
1790 Next_next_wire = 0;
1791 // Last_wire_Digit = 1000;
1792 // if (fDebug) cout<<" -------------0? Next_next_wire "<<Next_next_wire<<" wire "<<iwire<<endl;
1793 continue;
1794 }
1795
1796 // if (fDebug) cout<<" Next_next_wire "<<Next_next_wire<<endl;
1797 if (Nfirst[chNum][ipll] < 0 && DigitsArray_[chNum][ipll][iwire] > 0. && Next_next_wire) {
1798
1799 if (DigitsArray_[chNum][ipll][iwire] > DigitsArray_[chNum][ipll][iwire - 1]) {
1800 Next_next_wire = 0;
1801 // Last_wire_Digit = 1000;
1802 continue;
1803 } else {
1804 Nfirst[chNum][ipll] = iwire;
1805 ClusterSize_[chNum][ipll][Nclust_[chNum][ipll]] = 1;
1806 Coord_wire_[chNum][ipll][Nclust_[chNum][ipll]] =
1807 Nfirst[chNum][ipll] + 0.5 * ClusterSize_[chNum][ipll][Nclust_[chNum][ipll]];
1808 Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]] =
1809 (Coord_wire_[chNum][ipll][Nclust_[chNum][ipll]] - kMiddlePl) * dw;
1810 Next_next_wire = 1;
1811 // Last_wire_Digit = DigitsArray_[chNum][ipll][Nfirst[chNum][ipll]];
1812
1813 // made for the canonical sequence / x- v- u+ x+ v+ u-/ 0 1 5 have back reading
1814 if (ipll == 0 || ipll == 1 || ipll == 5)
1815 Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]] =
1816 -Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]];
1817
1818 if (fDebug) {
1819 hClusterSize.at(chNum)->Fill(ClusterSize_[chNum][ipll][Nclust_[chNum][ipll]]);
1820 // cout<<" Coord_xuv_["<< chNum<<"]["<<ipll<<"]["<<Nclust_[chNum][ipll]<<"]
1821 // "<<Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]]<<endl; cout<<"2 Nfirst "<<
1822 // Nfirst[chNum][ipll]<<" ClusterSize_
1823 // "<<ClusterSize_[chNum][ipll][Nclust_[chNum][ipll]]<<" Coord_wire_
1824 // "<<Coord_wire_[chNum][ipll][Nclust_[chNum][ipll]]<<endl;
1825
1826 hoccupancyXm.at(chNum)->Fill(Coord_xuv_[chNum][0][Nclust_[chNum][ipll]]);
1827 hoccupancyVm.at(chNum)->Fill(Coord_xuv_[chNum][1][Nclust_[chNum][ipll]]);
1828 hoccupancyUp.at(chNum)->Fill(Coord_xuv_[chNum][2][Nclust_[chNum][ipll]]);
1829 hoccupancyXp.at(chNum)->Fill(Coord_xuv_[chNum][3][Nclust_[chNum][ipll]]);
1830 hoccupancyVp.at(chNum)->Fill(Coord_xuv_[chNum][4][Nclust_[chNum][ipll]]);
1831 hoccupancyUm.at(chNum)->Fill(Coord_xuv_[chNum][5][Nclust_[chNum][ipll]]);
1832 }
1833
1834 Nclust_[chNum][ipll]++;
1835 Nfirst[chNum][ipll] = -1;
1836 } // else
1837 }
1838
1839 if (Nfirst[chNum][ipll] >= 0
1840 && (DigitsArray_[chNum][ipll][iwire] > 0
1841 && DigitsArray_[chNum][ipll][iwire] < DigitsArray_[chNum][ipll][Nfirst[chNum][ipll]]))
1842 {
1843 Nfirst[chNum][ipll] = iwire;
1844 if (fDebug)
1845 cout << " first Nfirst " << Nfirst[chNum][ipll] << endl;
1846 }
1847
1848 if (Nfirst[chNum][ipll] >= 0
1849 && (DigitsArray_[chNum][ipll][iwire] == 0
1850 || (DigitsArray_[chNum][ipll][iwire] > 0
1851 && DigitsArray_[chNum][ipll][iwire] > DigitsArray_[chNum][ipll][Nfirst[chNum][ipll]])))
1852 {
1853
1854 ClusterSize_[chNum][ipll][Nclust_[chNum][ipll]] = 1;
1855 Coord_wire_[chNum][ipll][Nclust_[chNum][ipll]] =
1856 Nfirst[chNum][ipll] + 0.5 * ClusterSize_[chNum][ipll][Nclust_[chNum][ipll]];
1857 Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]] =
1858 (Coord_wire_[chNum][ipll][Nclust_[chNum][ipll]] - kMiddlePl) * dw;
1859 if (DigitsArray_[chNum][ipll][iwire] > 0)
1860 Next_next_wire = 1;
1861
1862 // Last_wire_Digit = DigitsArray_[chNum][ipll][Nfirst[chNum][ipll]];
1863
1864 // made for the canonical sequence / x- v- u+ x+ v+ u-/ 0 1 5 have back reading
1865 if (ipll == 0 || ipll == 1 || ipll == 5)
1866 Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]] =
1867 -Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]];
1868 if (fDebug) {
1869 hClusterSize.at(chNum)->Fill(ClusterSize_[chNum][ipll][Nclust_[chNum][ipll]]);
1870 cout << " Coord_xuv_[" << chNum << "][" << ipll << "][" << Nclust_[chNum][ipll] << "] "
1871 << Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]] << " Coord_wire_ "
1872 << Coord_wire_[chNum][ipll][Nclust_[chNum][ipll]] << endl;
1873 // cout<<"1 Nfirst "<< Nfirst[chNum][ipll]<<" ClusterSize_
1874 // "<<ClusterSize_[chNum][ipll][Nclust_[chNum][ipll]]<<" Coord_wire_
1875 // "<<Coord_wire_[chNum][ipll][Nclust_[chNum][ipll]]<<endl; cout<<" iwire "<<iwire<<endl;
1876
1877 hoccupancyXm.at(chNum)->Fill(Coord_xuv_[chNum][0][Nclust_[chNum][ipll]]);
1878 hoccupancyVm.at(chNum)->Fill(Coord_xuv_[chNum][1][Nclust_[chNum][ipll]]);
1879 hoccupancyUp.at(chNum)->Fill(Coord_xuv_[chNum][2][Nclust_[chNum][ipll]]);
1880 hoccupancyXp.at(chNum)->Fill(Coord_xuv_[chNum][3][Nclust_[chNum][ipll]]);
1881 hoccupancyVp.at(chNum)->Fill(Coord_xuv_[chNum][4][Nclust_[chNum][ipll]]);
1882 hoccupancyUm.at(chNum)->Fill(Coord_xuv_[chNum][5][Nclust_[chNum][ipll]]);
1883 }
1884
1885 Nclust_[chNum][ipll]++;
1886 Nfirst[chNum][ipll] = -1;
1887 }
1888
1889 if (Nfirst[chNum][ipll] < 0 && DigitsArray_[chNum][ipll][iwire] > 0. && !Next_next_wire)
1890 Nfirst[chNum][ipll] = iwire;
1891
1892 if (iwire == 95 && Nfirst[chNum][ipll] > 0) {
1893 ClusterSize_[chNum][ipll][Nclust_[chNum][ipll]] = 1;
1894 Coord_wire_[chNum][ipll][Nclust_[chNum][ipll]] =
1895 Nfirst[chNum][ipll] + 0.5 * ClusterSize_[chNum][ipll][Nclust_[chNum][ipll]];
1896 Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]] =
1897 (Coord_wire_[chNum][ipll][Nclust_[chNum][ipll]] - kMiddlePl) * dw;
1898 if (fDebug)
1899 cout << " Coord_xuv_[" << chNum << "][" << ipll << "][" << Nclust_[chNum][ipll] << "] "
1900 << Coord_xuv_[chNum][ipll][Nclust_[chNum][ipll]] << endl;
1901 if (fDebug)
1902 cout << "3 Nfirst " << Nfirst[chNum][ipll] << " ClusterSize_ "
1903 << ClusterSize_[chNum][ipll][Nclust_[chNum][ipll]] << " Coord_wire_ "
1904 << Coord_wire_[chNum][ipll][Nclust_[chNum][ipll]] << endl;
1905 if (fDebug)
1906 cout << " iwire " << iwire << endl;
1907 Nclust_[chNum][ipll]++;
1908 Nfirst[chNum][ipll] = -1;
1909 }
1910 } // area with beam
1911 } // iwire
1912 // Num_layers_out_beam += Fired_layer_[chNum][ipll];
1913 } // else loaded area
1914 } // ipll
1915
1916 // if (fDebug) hNum_layers_out_beam_Ch.at(chNum)->Fill( Num_layers_out_beam);
1917 // }//else chamber 23
1918
1919} // Clustering
1920//----------------------------------------------------------------------
1921
1922//----------------- Combinatorial segment selection --------------------
1923void BmnMwpcHitFinder::SegmentFinder(Int_t chNum,
1924 Int_t** Nclust_,
1925 Double_t*** Coord_xuv_,
1926 Int_t*** ClusterSize_,
1927 Int_t* Nsegm,
1928 Double_t*** XVU_coor,
1929 Double_t*** Cluster_coor,
1930 Int_t** Nhits_Ch_,
1931 Int_t minHits,
1932 Short_t code,
1933 Int_t kBig_)
1934{
1935 if (fDebug)
1936 cout << "SegmentFinder: hexagon on the base of triad chNum: " << chNum << endl;
1937 // Coord_xuv_ - coordinates of all clusters
1938 Int_t minHits4_5;
1939 if (chNum > 1)
1940 minHits4_5 = minHits; // for run7
1941 else
1942 minHits4_5 = kMinHits_before_target;
1943 // if (fDebug) cout<<" chNum "<<chNum<<" minHits4_5 "<<minHits4_5<<endl;
1944
1945 Double_t min_for_triples = 5.; // minimum delta wires for tree planes
1946 Double_t min_for_conjugate = 3.; // minimum delta wires for conjugate plane
1947
1948 Double_t a_side = 12. * sq3 + 3 * dw; // wire half length //incorrect
1949
1950 // code : first triples is
1951 // 1 {X-, V-, U+}
1952 // 2 {X-, V+, U+}
1953 // 3 {X-, V-, U-}
1954 // 7 {X+, V-, U+}
1955 // 5 {X+, V+, U+}
1956 // 6 {X+, V-, U-}
1957 // 4 {X-, V+, U-}
1958 // 8 {X+, V+, U-}
1959
1960 Int_t x = 0, v = 1, u = 2, x1 = 3, v1 = 4, u1 = 5; // MK
1961
1962 switch (code) {
1963 case 1:
1964 x = 0;
1965 v = 1;
1966 u = 2;
1967 x1 = 3;
1968 v1 = 4;
1969 u1 = 5;
1970 break;
1971 case 2:
1972 x = 0;
1973 v = 4;
1974 u = 2;
1975 x1 = 3;
1976 v1 = 1;
1977 u1 = 5;
1978 break;
1979 case 3:
1980 x = 0;
1981 v = 1;
1982 u = 5;
1983 x1 = 3;
1984 v1 = 4;
1985 u1 = 2;
1986 break;
1987 case 7:
1988 x = 3;
1989 v = 1;
1990 u = 2;
1991 x1 = 0;
1992 v1 = 4;
1993 u1 = 5;
1994 break;
1995 case 5:
1996 x = 3;
1997 v = 4;
1998 u = 2;
1999 x1 = 0;
2000 v1 = 1;
2001 u1 = 5;
2002 break;
2003 case 6:
2004 x = 3;
2005 v = 1;
2006 u = 5;
2007 x1 = 0;
2008 v1 = 4;
2009 u1 = 2;
2010 break;
2011 case 4:
2012 x = 0;
2013 v = 4;
2014 u = 5;
2015 x1 = 3;
2016 v1 = 1;
2017 u1 = 2;
2018 break;
2019 case 8:
2020 x = 3;
2021 v = 4;
2022 u = 5;
2023 x1 = 0;
2024 v1 = 1;
2025 u1 = 2;
2026 break;
2027 }
2028
2029 if (Nsegm[chNum] > kBig_ - 2)
2030 return; // MP
2031
2032 for (Int_t ix = 0; ix < Nclust_[chNum][x]; ix++) {
2033
2034 for (Int_t iv = 0; iv < Nclust_[chNum][v]; iv++) {
2035
2036 for (Int_t iu = 0; iu < Nclust_[chNum][u]; iu++) {
2037
2038 if (Nsegm[chNum] > kBig_ - 2)
2039 return;
2040
2041 // --for repeated triples--
2042 Bool_t it_was = 0;
2043 if (Nsegm[chNum] > 0) {
2044 for (Int_t iseg = 0; iseg < Nsegm[chNum]; iseg++) {
2045 Bool_t it_was_x = 0;
2046 Bool_t it_was_v = 0;
2047 Bool_t it_was_u = 0;
2048
2049 if (XVU_coor[chNum][x][iseg] == Coord_xuv_[chNum][x][ix])
2050 it_was_x = 1;
2051 if (XVU_coor[chNum][v][iseg] == Coord_xuv_[chNum][v][iv])
2052 it_was_v = 1;
2053 if (XVU_coor[chNum][u][iseg] == Coord_xuv_[chNum][u][iu])
2054 it_was_u = 1;
2055
2056 it_was = it_was_x && it_was_v && it_was_u;
2057
2058 if (it_was) {
2059 break;
2060 }
2061 } // iseg
2062 } // Nsegm[chNum] > 0
2063
2064 if (it_was)
2065 continue;
2066 // --for repeated triples.
2067
2068 //--- main equation---// u + v - x < delta
2069 if (fabs(Coord_xuv_[chNum][u][iu] + Coord_xuv_[chNum][v][iv] - Coord_xuv_[chNum][x][ix])
2070 < min_for_triples * dw)
2071 {
2072
2073 // if (fDebug)cout<<" x " <<Coord_xuv_[chNum][x][ix]<<
2074 // " y_ux "<<(2*Coord_xuv_[chNum][u][iu] - Coord_xuv_[chNum][x][ix])/sq3<<
2075 // " vt " <<(2 *Coord_xuv_[chNum][u][iu] + Coord_xuv_[chNum][v][iv])/sq3<<
2076 // " ut " <<(2 *Coord_xuv_[chNum][x][ix] - Coord_xuv_[chNum][u][iu])/sq3<<" a_side
2077 // "<<a_side<<endl;
2078 //" ut " <<(Coord_xuv_[chNum][u][iu] + 2*Coord_xuv_[chNum][v][iv])/sq3<<endl;
2079
2080 if (fabs((2 * Coord_xuv_[chNum][u][iu] - Coord_xuv_[chNum][x][ix]) / sq3) < a_side
2081 && // y = (2u-x)/sq3
2082 fabs((2 * Coord_xuv_[chNum][u][iu] + Coord_xuv_[chNum][v][iv]) / sq3) < a_side
2083 && // vt = (2u+v)/sq3
2084 // fabs((2 *Coord_xuv_[chNum][x][ix] - Coord_xuv_[chNum][u][iu])/sq3) < a_side ){//
2085 fabs((Coord_xuv_[chNum][u][iu] + 2 * Coord_xuv_[chNum][v][iv]) / sq3) < a_side)
2086 { // ut = (u+2v)/sq3
2087
2088 if (fDebug)
2089 cout << " 3p-candidate new Nsegm " << endl;
2090 XVU_coor[chNum][x][Nsegm[chNum]] = Coord_xuv_[chNum][x][ix];
2091 XVU_coor[chNum][v][Nsegm[chNum]] = Coord_xuv_[chNum][v][iv];
2092 XVU_coor[chNum][u][Nsegm[chNum]] = Coord_xuv_[chNum][u][iu];
2093 Cluster_coor[chNum][x][Nsegm[chNum]] = ClusterSize_[chNum][x][ix];
2094 Cluster_coor[chNum][v][Nsegm[chNum]] = ClusterSize_[chNum][v][iv];
2095 Cluster_coor[chNum][u][Nsegm[chNum]] = ClusterSize_[chNum][u][iu];
2096
2097 XVU_coor[chNum][x1][Nsegm[chNum]] = -999.;
2098 XVU_coor[chNum][v1][Nsegm[chNum]] = -999.;
2099 XVU_coor[chNum][u1][Nsegm[chNum]] = -999.;
2100 Cluster_coor[chNum][x1][Nsegm[chNum]] = -1;
2101 Cluster_coor[chNum][v1][Nsegm[chNum]] = -1;
2102 Cluster_coor[chNum][u1][Nsegm[chNum]] = -1;
2103
2104 Nhits_Ch_[chNum][Nsegm[chNum]] = 3;
2105 // if (fDebug)cout<<" Nhits 3 = "<<Nhits_Ch_[chNum][Nsegm[chNum]]<<endl;
2106 // if (fDebug)cout<<endl;
2107
2108 //-- if 3p-candidate was look for conjugate coord
2109 if (XVU_coor[chNum][x][Nsegm[chNum]] != -999.) {
2110 Bool_t point_was = 0;
2111 Float_t Min_D = min_for_conjugate * dw;
2112 for (Int_t ix2 = 0; ix2 < Nclust_[chNum][x1]; ix2++) {
2113
2114 if (abs(XVU_coor[chNum][x][Nsegm[chNum]] - Coord_xuv_[chNum][x1][ix2]) < Min_D) {
2115 XVU_coor[chNum][x1][Nsegm[chNum]] = Coord_xuv_[chNum][x1][ix2];
2116 Cluster_coor[chNum][x1][Nsegm[chNum]] = ClusterSize_[chNum][x1][ix2];
2117 if (!point_was)
2118 Nhits_Ch_[chNum][Nsegm[chNum]]++; // 4 points
2119 point_was = 1;
2120 Min_D = abs(XVU_coor[chNum][x][Nsegm[chNum]] - Coord_xuv_[chNum][x1][ix2]);
2121 } // abs(
2122 } // ix2
2123 } // if it was
2124
2125 // if (fDebug)cout<<" Nhits 4 = "<<Nhits_Ch_[chNum][Nsegm[chNum]]<<endl;
2126 // if (fDebug)cout<<endl;
2127
2128 if (XVU_coor[chNum][v][Nsegm[chNum]] != -999.) {
2129 Bool_t point_was = 0;
2130 Float_t Min_D = min_for_conjugate * dw;
2131 for (Int_t iv2 = 0; iv2 < Nclust_[chNum][v1]; iv2++) {
2132
2133 if (abs(XVU_coor[chNum][v][Nsegm[chNum]] - Coord_xuv_[chNum][v1][iv2]) < Min_D) {
2134 XVU_coor[chNum][v1][Nsegm[chNum]] = Coord_xuv_[chNum][v1][iv2];
2135 Cluster_coor[chNum][v1][Nsegm[chNum]] = ClusterSize_[chNum][v1][iv2];
2136 if (!point_was)
2137 Nhits_Ch_[chNum][Nsegm[chNum]]++; // 5 points
2138 point_was = 1;
2139 Min_D = abs(XVU_coor[chNum][v][Nsegm[chNum]] - Coord_xuv_[chNum][v1][iv2]);
2140 } // abs(
2141 } // iv2
2142 } // if it was
2143 // if (fDebug)cout<<" Nhits 5 = "<<Nhits_Ch_[chNum][Nsegm[chNum]]<<endl;
2144 // if (fDebug)cout<<endl;
2145
2146 if (XVU_coor[chNum][u][Nsegm[chNum]] != -999.) {
2147 Bool_t point_was = 0;
2148 Float_t Min_D = min_for_conjugate * dw;
2149 for (Int_t iu2 = 0; iu2 < Nclust_[chNum][u1]; iu2++) {
2150 if (abs(XVU_coor[chNum][u][Nsegm[chNum]] - Coord_xuv_[chNum][u1][iu2]) < Min_D) {
2151 XVU_coor[chNum][u1][Nsegm[chNum]] = Coord_xuv_[chNum][u1][iu2];
2152 Cluster_coor[chNum][u1][Nsegm[chNum]] = ClusterSize_[chNum][u1][iu2];
2153 if (!point_was)
2154 Nhits_Ch_[chNum][Nsegm[chNum]]++; // 6 points
2155 point_was = 1;
2156 Min_D = abs(XVU_coor[chNum][u][Nsegm[chNum]] - Coord_xuv_[chNum][u1][iu2]);
2157 } // abs(
2158 } // iu2
2159 } // if it was
2160
2161 if (chNum > 1 && Nhits_Ch_[chNum][Nsegm[chNum]] > 4)
2162 minHits4_5 = minHits;
2163 if (Nsegm[chNum] > 15)
2164 minHits4_5 = 5;
2165 if (Nsegm[chNum] > 30)
2166 minHits4_5 = 6;
2167
2168 // if (fDebug) cout<<" chNum "<<chNum<<" minHits4_5 "<<minHits4_5<<" Nsegm "<<
2169 // Nsegm[chNum]<<endl;
2170
2171 if (Nhits_Ch_[chNum][Nsegm[chNum]] >= minHits4_5) {
2172 // if (fDebug) cout<<endl;
2173 // if (fDebug) cout<<" Nsegm "<< Nsegm[chNum] <<"
2174 // Nhits_Ch_["<<chNum<<"]["<<Nsegm[chNum]<<"] "<<Nhits_Ch_[chNum][Nsegm[chNum]]<<"
2175 // minHits4_5 "<< minHits4_5 <<endl; if (fDebug) cout<<" coord
2176 // "<<XVU_coor[chNum][0][Nsegm[chNum]]<<" "<<XVU_coor[chNum][1][Nsegm[chNum]]<<"
2177 // "<<XVU_coor[chNum][2][Nsegm[chNum]]<<" "<<XVU_coor[chNum][3][Nsegm[chNum]]<<"
2178 // "<<XVU_coor[chNum][4][Nsegm[chNum]]<<" "<<XVU_coor[chNum][5][Nsegm[chNum]]<<" "<<endl;
2179
2180 Nsegm[chNum]++; // go to next segment
2181 }
2182 if (Nhits_Ch_[chNum][Nsegm[chNum]] < minHits4_5) {
2183 //--segment deleting
2184 Nhits_Ch_[chNum][Nsegm[chNum]] = 0;
2185
2186 XVU_coor[chNum][x][Nsegm[chNum]] = -999.;
2187 XVU_coor[chNum][v][Nsegm[chNum]] = -999.;
2188 XVU_coor[chNum][u][Nsegm[chNum]] = -999.;
2189 XVU_coor[chNum][x1][Nsegm[chNum]] = -999.;
2190 XVU_coor[chNum][v1][Nsegm[chNum]] = -999.;
2191 XVU_coor[chNum][u1][Nsegm[chNum]] = -999.;
2192 Cluster_coor[chNum][x][Nsegm[chNum]] = -1;
2193 Cluster_coor[chNum][v][Nsegm[chNum]] = -1;
2194 Cluster_coor[chNum][u][Nsegm[chNum]] = -1;
2195 Cluster_coor[chNum][x1][Nsegm[chNum]] = -1;
2196 Cluster_coor[chNum][v1][Nsegm[chNum]] = -1;
2197 Cluster_coor[chNum][u1][Nsegm[chNum]] = -1;
2198 } // else
2199 } // +-a cm
2200 } // u + v - x = delta
2201 if (Nsegm[chNum] > kBig_ - 2)
2202 break;
2203 } // iu
2204 if (Nsegm[chNum] > kBig_ - 2)
2205 break;
2206 } // iv
2207 if (Nsegm[chNum] > kBig_ - 2)
2208 break;
2209 } // ix
2210 if (Nsegm[chNum] > kBig_ - 2)
2211 return;
2212 // if (fDebug) cout<<" --end SegmentFinder: Nsegm["<<chNum<<"] "<<Nsegm[chNum]<<endl;
2213
2214} // SegmentFinder1
2215//----------------------------------------------------------------------
2216
2217//----------------- Combinatorial segment selection --------------------
2218void BmnMwpcHitFinder::SegmentFinder2coord(Int_t chNum,
2219 Int_t** Nclust_,
2220 Double_t*** Coord_xuv_,
2221 Int_t*** ClusterSize_,
2222 Int_t* Nsegm,
2223 Double_t*** XVU_coor,
2224 Double_t*** Cluster_coor,
2225 Int_t** Nhits_Ch_,
2226 Int_t minHits,
2227 Short_t code,
2228 Int_t kBig_,
2229 Int_t* nbest_seg,
2230 Double_t*** coor_seg)
2231{
2232
2233 if (fDebug)
2234 cout << " SegmentFinder2coord " << endl;
2235 if (fDebug)
2236 cout << " Nsegm[" << chNum << "] " << Nsegm[chNum] << endl;
2237 Int_t minHits4 = 4;
2238 Double_t min_for_conjugate = 3.; // minimum delta wires for conjugate plane
2239 Double_t a_side = 12. * sq3 + 3 * dw; // wire half length
2240
2241 Int_t f, l, f1, l1;
2242 /* for 3-4 points
2243 switch (code) {
2244 case 1: f = 0; l = 1; f1 = 3; l1 = 4; break;
2245 case 2: f = 3; l = 1; f1 = 0; l1 = 4; break;
2246 case 3: f = 0; l = 4; f1 = 3; l1 = 1; break;
2247 case 4: f = 3; l = 4; f1 = 0; l1 = 1; break;
2248 case 5: f = 0; l = 2; f1 = 3; l1 = 5; break;
2249 case 6: f = 3; l = 2; f1 = 0; l1 = 5; break;
2250 case 7: f = 0; l = 5; f1 = 3; l1 = 2; break;
2251 case 8: f = 3; l = 5; f1 = 0; l1 = 2; break;
2252 case 9: f = 1; l = 2; f1 = 4; l1 = 5; break;
2253 case 10: f = 4; l = 2; f1 = 1; l1 = 5; break;
2254 case 11: f = 1; l = 5; f1 = 4; l1 = 2; break;
2255 case 12: f = 4; l = 5; f1 = 1; l1 = 2; break;
2256 }
2257 */
2258 // only 4 points
2259 switch (code) {
2260 case 1:
2261 f = 0;
2262 l = 4;
2263 f1 = 3;
2264 l1 = 1;
2265 break; // XV
2266 case 2:
2267 f = 0;
2268 l = 5;
2269 f1 = 3;
2270 l1 = 2;
2271 break; // XU
2272 case 3:
2273 f = 4;
2274 l = 5;
2275 f1 = 1;
2276 l1 = 2;
2277 break; // VU
2278 }
2279
2280 if (Nsegm[chNum] > kBig_ - 2)
2281 return; // MP
2282 // if (fDebug) cout<<" Nclustf "<<Nclust_[chNum][f]<<" Nclustl "<<Nclust_[chNum][l]<<endl;
2283
2284 for (Int_t jf = 0; jf < Nclust_[chNum][f]; jf++) {
2285
2286 Bool_t it_was_f = 0;
2287 for (Int_t ib = 0; ib < nbest_seg[chNum]; ib++) {
2288 if (coor_seg[chNum][f][ib] == Coord_xuv_[chNum][f][jf]) {
2289 it_was_f = 1;
2290 break;
2291 }
2292 }
2293 if (it_was_f)
2294 continue;
2295
2296 for (Int_t il = 0; il < Nclust_[chNum][l]; il++) {
2297
2298 if (Nsegm[chNum] > kBig_ - 2)
2299 return;
2300 Bool_t it_was_l = 0;
2301 for (Int_t ib = 0; ib < nbest_seg[chNum]; ib++) {
2302 if (coor_seg[chNum][l][ib] == Coord_xuv_[chNum][l][il]) {
2303 it_was_l = 1;
2304 break;
2305 }
2306 }
2307 if (it_was_l)
2308 continue;
2309
2310 // --for repeated duplets--
2311 /* // 30/04/2021
2312 Bool_t it_was = 0;
2313
2314 if (Nsegm[chNum] > 0) {
2315 for (Int_t iseg = 0; iseg < Nsegm[chNum]; iseg++) {
2316 it_was_f = 0;
2317 it_was_l = 0;
2318
2319 //cout<<" XVU "<<XVU_coor[chNum][l][iseg]<<" cand "<<Coord_xuv_[chNum][l][il]<<endl;
2320
2321 //if (XVU_coor[chNum][f][iseg] == Coord_xuv_[chNum][f][jf]) it_was_f = 1;
2322 //if (XVU_coor[chNum][l][iseg] == Coord_xuv_[chNum][l][il]) it_was_l = 1;
2323 // it_was = it_was_f ^ it_was_l;
2324 //cout<<" it_was "<<it_was<<" it_was_f "<<it_was_f<<" it_was_l "<<it_was_l<<endl;
2325
2326 if (it_was) {
2327 break;
2328 }
2329 } // iseg
2330 } // Nsegm[chNum] > 0
2331 if (it_was) continue;
2332 */
2333 // --for repeated duplets.
2334
2335 //--- equation---//
2336 Double_t ft, lt;
2337 if (code == 1) { //--xv--
2338 ft = (Coord_xuv_[chNum][f][jf] - 2 * Coord_xuv_[chNum][l][il]) / sq3;
2339 lt = (2 * Coord_xuv_[chNum][f][jf] - Coord_xuv_[chNum][l][il]) / sq3;
2340 }
2341 if (code == 2) { //--xu--
2342 ft = (-Coord_xuv_[chNum][f][jf] + 2 * Coord_xuv_[chNum][l][il]) / sq3;
2343 lt = (2 * Coord_xuv_[chNum][f][jf] - Coord_xuv_[chNum][l][il]) / sq3;
2344 }
2345 if (code == 3) { //--uv--
2346 ft = (Coord_xuv_[chNum][f][jf] + 2 * Coord_xuv_[chNum][l][il]) / sq3;
2347 lt = (2 * Coord_xuv_[chNum][f][jf] + Coord_xuv_[chNum][l][il]) / sq3;
2348 }
2349
2350 if (fDebug)
2351 cout << " pl f:" << f << " coord " << Coord_xuv_[chNum][f][jf] << " pl l: " << l << " coord "
2352 << Coord_xuv_[chNum][l][il] << " ft " << ft << " lt " << lt << " a_side " << a_side << endl;
2353
2354 if (fabs(ft) < a_side && fabs(lt) < a_side) {
2355 // 2p-candidate new Nsegm
2356 XVU_coor[chNum][f][Nsegm[chNum]] = Coord_xuv_[chNum][f][jf];
2357 XVU_coor[chNum][l][Nsegm[chNum]] = Coord_xuv_[chNum][l][il];
2358 Cluster_coor[chNum][f][Nsegm[chNum]] = ClusterSize_[chNum][f][jf];
2359 Cluster_coor[chNum][l][Nsegm[chNum]] = ClusterSize_[chNum][l][il];
2360
2361 Nhits_Ch_[chNum][Nsegm[chNum]] = 2;
2362 // if (fDebug)cout<<" Nhits 2 = "<<Nhits_Ch_[chNum][Nsegm[chNum]]<<endl;
2363 // if (fDebug)cout<<endl;
2364
2365 //-- if 3p-candidate was look for conjugate coord
2366 if (XVU_coor[chNum][f][Nsegm[chNum]] != -999.) {
2367 Bool_t point_was = 0;
2368 Float_t Min_D = min_for_conjugate * dw;
2369 for (Int_t ix2 = 0; ix2 < Nclust_[chNum][f1]; ix2++) {
2370
2371 if (abs(XVU_coor[chNum][f][Nsegm[chNum]] - Coord_xuv_[chNum][f1][ix2]) < Min_D) {
2372 XVU_coor[chNum][f1][Nsegm[chNum]] = Coord_xuv_[chNum][f1][ix2];
2373 Cluster_coor[chNum][f1][Nsegm[chNum]] = ClusterSize_[chNum][f1][ix2];
2374 if (!point_was)
2375 Nhits_Ch_[chNum][Nsegm[chNum]]++; // 3 points
2376 point_was = 1;
2377 Min_D = abs(XVU_coor[chNum][f][Nsegm[chNum]] - Coord_xuv_[chNum][f1][ix2]);
2378 } // abs(
2379 } // ix2
2380 } // if it was
2381 // if (fDebug)cout<<" Nhits 3 = "<<Nhits_Ch_[chNum][Nsegm[chNum]]<<endl;
2382 // if (fDebug)cout<<endl;
2383
2384 if (XVU_coor[chNum][l][Nsegm[chNum]] != -999.) {
2385 Bool_t point_was = 0;
2386 Float_t Min_D = min_for_conjugate * dw;
2387 for (Int_t iv2 = 0; iv2 < Nclust_[chNum][l1]; iv2++) {
2388
2389 if (abs(XVU_coor[chNum][l][Nsegm[chNum]] - Coord_xuv_[chNum][l1][iv2]) < Min_D) {
2390 XVU_coor[chNum][l1][Nsegm[chNum]] = Coord_xuv_[chNum][l1][iv2];
2391 Cluster_coor[chNum][l1][Nsegm[chNum]] = ClusterSize_[chNum][l1][iv2];
2392 if (!point_was)
2393 Nhits_Ch_[chNum][Nsegm[chNum]]++; // 4 points
2394 point_was = 1;
2395 Min_D = abs(XVU_coor[chNum][l][Nsegm[chNum]] - Coord_xuv_[chNum][l1][iv2]);
2396 } // abs(
2397 } // iv2
2398 } // if it was
2399 // if (fDebug)cout<<" Nhits 4 = "<<Nhits_Ch_[chNum][Nsegm[chNum]]<<endl;
2400 // if (fDebug)cout<<endl;
2401
2402 if (Nhits_Ch_[chNum][Nsegm[chNum]] == minHits4) {
2403 if (fDebug)
2404 cout << endl;
2405 if (fDebug)
2406 cout << " Nsegm " << Nsegm[chNum] << " Nhits_Ch_[" << chNum << "][" << Nsegm[chNum] << "] "
2407 << Nhits_Ch_[chNum][Nsegm[chNum]] << endl;
2408 if (fDebug)
2409 cout << " coord " << XVU_coor[chNum][0][Nsegm[chNum]] << " " << XVU_coor[chNum][1][Nsegm[chNum]]
2410 << " " << XVU_coor[chNum][2][Nsegm[chNum]] << " " << XVU_coor[chNum][3][Nsegm[chNum]]
2411 << " " << XVU_coor[chNum][4][Nsegm[chNum]] << " " << XVU_coor[chNum][5][Nsegm[chNum]]
2412 << " " << endl;
2413
2414 Nsegm[chNum]++; // go to next segment
2415 }
2416
2417 if (Nhits_Ch_[chNum][Nsegm[chNum]] < minHits4) {
2418 //--segment deleting
2419 Nhits_Ch_[chNum][Nsegm[chNum]] = 0;
2420
2421 XVU_coor[chNum][f][Nsegm[chNum]] = -999.;
2422 XVU_coor[chNum][l][Nsegm[chNum]] = -999.;
2423 XVU_coor[chNum][f1][Nsegm[chNum]] = -999.;
2424 XVU_coor[chNum][l1][Nsegm[chNum]] = -999.;
2425 Cluster_coor[chNum][f][Nsegm[chNum]] = -1;
2426 Cluster_coor[chNum][l][Nsegm[chNum]] = -1;
2427 Cluster_coor[chNum][f1][Nsegm[chNum]] = -1;
2428 Cluster_coor[chNum][l1][Nsegm[chNum]] = -1;
2429 } //
2430
2431 } // if (fabs(ft) < 12.5 && fabs(lt) < 12.5)
2432
2433 if (Nsegm[chNum] > kBig_ - 2)
2434 break;
2435 } // if
2436 if (Nsegm[chNum] > kBig_ - 2)
2437 break;
2438 } // il
2439 if (Nsegm[chNum] > kBig_ - 2)
2440 return;
2441
2442 if (fDebug)
2443 cout << " -- SegmentFinder2coord: Nsegm[" << chNum << "] " << Nsegm[chNum] << endl;
2444
2445} // SegmentFinder2coord
2446
2447void BmnMwpcHitFinder::ProcessSegments2coord(Int_t chNum,
2448 Int_t* Nsegm,
2449 Double_t*** XVU_coor,
2450 Double_t*** Cluster_coor,
2451 Int_t** Nhits_Ch_,
2452 Float_t** z_loc,
2453 Int_t Min_hits,
2454 Double_t sigma_,
2455 Double_t kChi2_Max_,
2456 Int_t** Nhits_seg_,
2457 Double_t** Chi2_ndf_seg_,
2458 Double_t*** coor_seg,
2459 Double_t*** cluster_seg,
2460 Double_t*** Par_ab_seg,
2461 Int_t* Nbest_seg_,
2462 Int_t* Nlay_w_wires_,
2463 Double_t**** sigma2_s)
2464{
2465 if (fDebug)
2466 cout << " -- ProcessSegments2coord -- " << endl;
2467 // Double_t dx_;
2468 Double_t Chi2[kNChambers][kBig];
2469 Double_t Par_ab[kNChambers][4][kBig];
2470 Double_t x1, x2, y1, y2, zx1, zx2, zy1, zy2;
2471 if (fDebug)
2472 cout << " Nbest_seg_[chNum] " << Nbest_seg_[chNum] << endl;
2473
2474 for (Int_t iseg = 0; iseg < Nsegm[chNum]; iseg++) { //---cycle by segments
2475 if (Nbest_seg_[chNum] > kBig - 2)
2476 return;
2477
2478 //---Chi2 calculating---
2479 Chi2[chNum][iseg] = 0.;
2480 x1 = -999.;
2481 x2 = -999.;
2482 y1 = -999.;
2483 y2 = -999.;
2484 zx1 = 0.;
2485 zx2 = 0.;
2486 zy1 = 0.;
2487 zy2 = 0.;
2488 Bool_t x_was = 0, u_was = 0, v_was = 0;
2489 //---x---
2490 if (XVU_coor[chNum][0][iseg] > -900. && XVU_coor[chNum][3][iseg] > -900.) { // x
2491
2492 // dx_ = XVU_coor[chNum][0][iseg] - XVU_coor[chNum][3][iseg];
2493 // Chi2[chNum][iseg] += dx_*dx_;
2494 if (fDebug)
2495 cout << " iseg " << iseg << " x1 " << XVU_coor[chNum][0][iseg] << " x2 " << XVU_coor[chNum][3][iseg]
2496 << endl;
2497 x1 = XVU_coor[chNum][0][iseg];
2498 x2 = XVU_coor[chNum][3][iseg];
2499 zx1 = z_loc[chNum][0];
2500 zx2 = z_loc[chNum][3];
2501 zy1 = z_loc[chNum][0];
2502 zy2 = z_loc[chNum][3];
2503 x_was = 1;
2504 }
2505 //---v---
2506 if (XVU_coor[chNum][4][iseg] > -999. && XVU_coor[chNum][1][iseg] > -999.) { // v
2507
2508 // dx_ = XVU_coor[chNum][4][iseg] - XVU_coor[chNum][1][iseg];
2509 // Chi2[chNum][iseg] += dx_*dx_;
2510 if (fDebug)
2511 cout << " iseg " << iseg << " v1 " << XVU_coor[chNum][4][iseg] << " v2 " << XVU_coor[chNum][1][iseg]
2512 << endl;
2513 v_was = 1;
2514
2515 zy1 += z_loc[chNum][4];
2516 zy2 += z_loc[chNum][1];
2517 if (x_was) {
2518 y1 = (x1 - 2 * XVU_coor[chNum][4][iseg]) / sq3;
2519 y2 = (x2 - 2 * XVU_coor[chNum][1][iseg]) / sq3;
2520 } else {
2521 zx1 = z_loc[chNum][4];
2522 zx2 = z_loc[chNum][1];
2523 }
2524 }
2525 //---u---
2526 if (XVU_coor[chNum][5][iseg] > -999. && XVU_coor[chNum][2][iseg] > -999.) { // u
2527
2528 // dx_ = XVU_coor[chNum][5][iseg] - XVU_coor[chNum][2][iseg];
2529 // Chi2[chNum][iseg] += dx_*dx_;
2530 if (fDebug)
2531 cout << " iseg " << iseg << " u1 " << XVU_coor[chNum][5][iseg] << " u2 " << XVU_coor[chNum][2][iseg]
2532 << endl;
2533 u_was = 1;
2534 zy1 += z_loc[chNum][5];
2535 zy2 += z_loc[chNum][2];
2536
2537 if (x_was) {
2538 y1 = (2 * XVU_coor[chNum][5][iseg] - x1) / sq3;
2539 y2 = (2 * XVU_coor[chNum][2][iseg] - x2) / sq3;
2540 }
2541 if (v_was) { // not x_was
2542 x1 = XVU_coor[chNum][5][iseg] + XVU_coor[chNum][4][iseg];
2543 x2 = XVU_coor[chNum][2][iseg] + XVU_coor[chNum][1][iseg];
2544
2545 y1 = (XVU_coor[chNum][5][iseg] - XVU_coor[chNum][4][iseg]) / sq3;
2546 y2 = (XVU_coor[chNum][2][iseg] - XVU_coor[chNum][1][iseg]) / sq3;
2547
2548 zx1 += z_loc[chNum][5];
2549 zx2 += z_loc[chNum][2];
2550 zx1 = zx1 / 2.;
2551 zx2 = zx2 / 2.;
2552 }
2553 }
2554
2555 zy1 = zy1 / 2.;
2556 zy2 = zy2 / 2.;
2557
2558 Double_t Ax = (x1 - x2) / (zx1 - zx2);
2559 Double_t Ay = (y1 - y2) / (zy1 - zy2);
2560 //--bisectrix--
2561 if (!x_was) {
2562 Double_t y_t = Ay * (Z0_SRC - ChZ[chNum]) + y1;
2563 if (fDebug)
2564 cout << " y_t " << y_t << endl;
2565 Chi2[chNum][iseg] = fabs(y_t + 6.1); // run7: in local system yv = -6.1 cm
2566 }
2567 if (!u_was) {
2568 Double_t u_t1 = (XVU_coor[chNum][0][iseg] + XVU_coor[chNum][4][iseg]) / sq3;
2569 Double_t u_t2 = (XVU_coor[chNum][3][iseg] + XVU_coor[chNum][1][iseg]) / sq3;
2570 Double_t Au = (u_t1 - u_t2) / (zy1 - zy2);
2571 Double_t u_t = Au * (Z0_SRC - ChZ[chNum]) + u_t1;
2572 Chi2[chNum][iseg] = fabs(u_t - 3.0); // run7: in local system uv = 3.0 cm
2573 if (fDebug)
2574 cout << " u_t " << u_t << endl;
2575 }
2576 if (!v_was) {
2577 Double_t v_t1 = (XVU_coor[chNum][0][iseg] + XVU_coor[chNum][5][iseg]) / sq3;
2578 Double_t v_t2 = (XVU_coor[chNum][3][iseg] + XVU_coor[chNum][2][iseg]) / sq3;
2579 Double_t Av = (v_t1 - v_t2) / (zy1 - zy2);
2580 Double_t v_t = Av * (Z0_SRC - ChZ[chNum]) + v_t1;
2581 Chi2[chNum][iseg] = fabs(v_t + 3.3); // run7: in local system vv = -3.3 cm
2582 if (fDebug)
2583 cout << " v_t " << v_t << endl;
2584 }
2585
2586 // if (fDebug) cout<<" zy1= "<<zy1<<" zy2= "<<zy2<<" zx1= "<<zx1<<" zx2= "<<zx2<<" x1= "<<x1<<" x2= "<<x2<<" y1=
2587 // "<<y1<<" y2= "<<y2<<endl;
2588 if (fDebug)
2589 cout << " iseg " << iseg << " Chi2 " << Chi2[chNum][iseg] << endl;
2590
2591 Par_ab[chNum][0][iseg] = Ax;
2592 Par_ab[chNum][1][iseg] = (x1 + x2) / 2.;
2593 Par_ab[chNum][2][iseg] = Ay;
2594 Par_ab[chNum][3][iseg] = (y1 + y2) / 2.;
2595
2596 } //---cycle by segments
2597
2598 //--choice of min chi2
2599 Int_t vec_ind_best_seg[kmaxSeg];
2600 Int_t Nbest_vec_ind = 0;
2601
2602 for (Int_t ind_best = 0; ind_best < kmaxSeg; ind_best++) {
2603 Double_t Chi2_best = 999.;
2604 Int_t iseg_best = -1;
2605 vec_ind_best_seg[ind_best] = -1;
2606
2607 for (Int_t iseg = 0; iseg < Nsegm[chNum]; iseg++) {
2608
2609 if (Chi2[chNum][iseg] >= Chi2_best)
2610 continue;
2611
2612 Bool_t best_was = 0;
2613 if (Nbest_vec_ind > 0) {
2614 for (Int_t ii = 0; ii < Nbest_vec_ind; ii++) {
2615 if (iseg == vec_ind_best_seg[ii]) {
2616 best_was = 1;
2617 break;
2618 }
2619 }
2620 }
2621 if (best_was)
2622 continue;
2623 Chi2_best = Chi2[chNum][iseg];
2624 iseg_best = iseg;
2625 } // iseg
2626
2627 if (iseg_best == -1)
2628 continue;
2629 if (fDebug)
2630 cout << " ind_best " << ind_best << " Chi2_best " << Chi2_best << " iseg_best " << iseg_best << endl;
2631 vec_ind_best_seg[ind_best] = iseg_best;
2632 Nbest_vec_ind++;
2633
2634 if (fDebug)
2635 cout << "-- reject(common points)" << endl;
2636 for (Int_t iseg = 0; iseg < Nsegm[chNum]; iseg++) {
2637 // if (iseg == iseg_best)continue;
2638 Bool_t best_was = 0;
2639 if (Nbest_vec_ind > 0) {
2640 for (Int_t ii = 0; ii < Nbest_vec_ind; ii++) {
2641 if (iseg == vec_ind_best_seg[ii]) {
2642 best_was = 1;
2643 break;
2644 }
2645 }
2646 }
2647 if (best_was)
2648 continue;
2649
2650 for (Int_t i1 = 0; i1 < kNPlanes; i1++) {
2651 if (XVU_coor[chNum][i1][iseg] > -999.) {
2652 if (fabs(XVU_coor[chNum][i1][iseg] - XVU_coor[chNum][i1][iseg_best]) < dw_half) { // 3*dw_half
2653 if (fDebug)
2654 cout << " reject point pl:" << i1 << " XUV " << XVU_coor[chNum][i1][iseg]
2655 << " iseg = " << iseg << endl;
2656 // Nhits_Ch_[chNum][iseg] = 0;
2657 XVU_coor[chNum][i1][iseg] = -999.;
2658 Chi2[chNum][iseg] = 999.; // Nhits_Ch_[chNum][iseg] = 0;
2659 break;
2660 }
2661 }
2662 } // i1
2663 } // iseg
2664 if (fDebug)
2665 cout << " Nbest_vec_ind " << Nbest_vec_ind << endl;
2666 if (fDebug)
2667 cout << " Chi2_best " << Chi2_best << endl;
2668 } // ind_best
2669
2670 //-----bests writing---
2671
2672 if (Nbest_seg_[chNum] < kmaxSeg) {
2673 for (Int_t iseg = 0; iseg < Nsegm[chNum]; iseg++) {
2674
2675 if (Chi2[chNum][iseg] < 900.) {
2676 Par_ab_seg[chNum][0][Nbest_seg_[chNum]] = Par_ab[chNum][0][iseg];
2677 Par_ab_seg[chNum][1][Nbest_seg_[chNum]] = Par_ab[chNum][1][iseg];
2678 Par_ab_seg[chNum][2][Nbest_seg_[chNum]] = Par_ab[chNum][2][iseg];
2679 Par_ab_seg[chNum][3][Nbest_seg_[chNum]] = Par_ab[chNum][3][iseg];
2680 // cout<<" Par_ab_seg[chNum][0][iseg] "<<Par_ab_seg[chNum][0][iseg]<<endl;
2681 if (fDebug)
2682 cout << " iseg " << iseg << " Chi2 " << Chi2[chNum][iseg] << endl;
2683
2684 if (fDebug) {
2685 hpar_Ax_Ch.at(chNum)->Fill(Par_ab_seg[chNum][0][Nbest_seg_[chNum]]);
2686 hpar_Bx_Ch.at(chNum)->Fill(Par_ab_seg[chNum][1][Nbest_seg_[chNum]]);
2687 hpar_Ay_Ch.at(chNum)->Fill(Par_ab_seg[chNum][2][Nbest_seg_[chNum]]);
2688 hpar_By_Ch.at(chNum)->Fill(Par_ab_seg[chNum][3][Nbest_seg_[chNum]]);
2689 }
2690 for (int ipla = 0; ipla < kNPlanes; ipla++) {
2691 // cout<<" XVU "<<XVU_coor[chNum][ipla][iseg]<<endl;
2692 coor_seg[chNum][ipla][Nbest_seg_[chNum]] = XVU_coor[chNum][ipla][iseg];
2693 cluster_seg[chNum][ipla][Nbest_seg_[chNum]] = Cluster_coor[chNum][ipla][iseg];
2694 }
2695
2696 Nhits_seg_[chNum][Nbest_seg_[chNum]] = Nhits_Ch_[chNum][iseg];
2697 Chi2_ndf_seg_[chNum][Nbest_seg_[chNum]] = Chi2[chNum][iseg];
2698 for (int i1 = 0; i1 < 4; i1++) {
2699 for (int j1 = 0; j1 < 4; j1++) {
2700 // sigma2_s[chNum][Nbest_seg_[chNum]][i1][j1] = sigma2[i1][j1];//?
2701 }
2702 }
2703 Nbest_seg_[chNum]++;
2704 }
2705 }
2706 } // kmaxSeg
2707
2708 if (fDebug)
2709 cout << "---end cycle by segments " << endl;
2710 if (fDebug)
2711 cout << " Nbest_seg_[chNum] " << Nbest_seg_[chNum] << endl;
2712 if (fDebug)
2713 cout << endl;
2714 for (Int_t iseg = 0; iseg < Nbest_seg_[chNum]; iseg++) { // print
2715 if (fDebug)
2716 cout << " iseg " << iseg << " Ax " << Par_ab_seg[chNum][0][iseg] << " bx " << Par_ab_seg[chNum][1][iseg]
2717 << " Ay " << Par_ab_seg[chNum][2][iseg] << " by " << Par_ab_seg[chNum][3][iseg] << " Chi2 "
2718 << Chi2_ndf_seg_[chNum][iseg] << " Nhits " << Nhits_seg_[chNum][iseg] << endl;
2719 }
2720}
2721
2722//------------ Segment fitting & finding the best track-candidate ------
2723void BmnMwpcHitFinder::ProcessSegments(Int_t chNum,
2724 Int_t* Nsegm,
2725 Double_t*** XVU_coor,
2726 Double_t*** Cluster_coor,
2727 Int_t** Nhits_Ch_,
2728 Float_t** z_loc,
2729 Int_t Min_hits,
2730 Double_t sigma_,
2731 Double_t kChi2_Max_,
2732 Int_t** Nhits_seg_,
2733 Double_t** Chi2_ndf_seg_,
2734 Double_t*** coor_seg,
2735 Double_t*** cluster_seg,
2736 Double_t*** Par_ab_seg,
2737 Int_t* Nbest_seg_,
2738 Int_t* Nlay_w_wires_,
2739 Double_t**** sigma2_s)
2740{
2741
2742 // segments seg[kNChambers][kBig];
2743
2744 Double_t par_ab[kNChambers][kNPlanes][kBig];
2745 Double_t dx_[kNPlanes];
2746 Double_t Chi2[kNChambers][kBig];
2747 Double_t Chi2_ndf[kNChambers][kBig];
2748 Double_t sigm[kNPlanes], sigm2_[kNPlanes];
2749 Int_t h1[kNPlanes];
2750 Double_t bmatr_seg[kBig][4][4];
2751 for (Int_t i = 0; i < kBig; i++) {
2752 for (Int_t i1 = 0; i1 < 4; i1++) {
2753 for (Int_t j1 = 0; j1 < 4; j1++) {
2754 bmatr_seg[i][i1][j1] = 0.;
2755 }
2756 }
2757 }
2758
2759 Int_t Min_hits6;
2760 // Double_t delta = 3*dw;
2761 // Double_t Chi2_Super_min = 0.001;
2762 Double_t x_target, y_target;
2763 Double_t xy_t_max[4] = {30., 35., 40., 45.};
2764
2765 if (chNum > 1)
2766 Min_hits6 = Min_hits; // for run7
2767 else
2768 Min_hits6 = kMinHits_before_target;
2769
2770 // cout<<" Nlay_w_wires["<<chNum<<"] "<<Nlay_w_wires_[chNum]<<endl;
2771 // TEMPORARY OUT SEGMENTS ARRY
2772 int OutSegCount = 0;
2773 segments OutSegArray[kmaxSeg];
2774 // cout<<" Nlay_w_wires_["<<chNum<<"] "<<Nlay_w_wires_[chNum]<<endl;
2775 if (Nlay_w_wires_[chNum] < 4)
2776 return;
2777
2778 if (chNum > 1 && Min_hits < 6) {
2779 if (Nlay_w_wires_[chNum] == 4)
2780 Min_hits6 = 4;
2781 }
2782
2783 if (Nsegm[chNum] > kBig - 2)
2784 return;
2785 /*
2786 if (fDebug) cout<<" ProcessSegments "<<endl;
2787 for (Int_t iseg = 0; iseg < Nsegm[chNum]; iseg++) { //---print number of hits
2788 if (fDebug) cout<<" iseg "<<iseg<<" Nhits "<<Nhits_Ch_[chNum][iseg]<<endl;
2789 }
2790 */
2791
2792 for (Int_t Nhitm = kNPlanes; Nhitm > Min_hits6 - 1; Nhitm--) { // -- first view 6 points
2793 Bool_t ifNhitm = 0;
2794
2795 if (Nhitm < Min_hits6)
2796 break;
2797
2798 for (Int_t iseg = 0; iseg < Nsegm[chNum]; iseg++) { //---cycle by segments
2799 if (Nbest_seg_[chNum] > kBig - 2)
2800 return;
2801 ifNhitm = 1;
2802
2803 if (Nhits_Ch_[chNum][iseg] != Nhitm)
2804 continue;
2805 // if (fDebug) cout<<" iseg "<<iseg<<" Nhits "<<Nhits_Ch_[chNum][iseg]<<endl;
2806 if (Nhits_Ch_[chNum][iseg] != 0) {
2807 Nhits_Ch_[chNum][iseg] = 0;
2808 for (Int_t i1 = 0; i1 < kNPlanes; i1++) {
2809 if (XVU_coor[chNum][i1][iseg] > -900.)
2810 Nhits_Ch_[chNum][iseg]++;
2811 }
2812 if (Nhits_Ch_[chNum][iseg] < Min_hits6)
2813 Nhits_Ch_[chNum][iseg] = 0; // control shot
2814 }
2815
2816 // if (fDebug) cout<<"after recalc nhits iseg "<<iseg<<" Nhits "<<Nhits_Ch_[chNum][iseg]<<endl;
2817 //-----!!!----
2818 if (Nhits_Ch_[chNum][iseg] == 0)
2819 continue; // go to next segment!
2820
2821 for (Int_t i = 0; i < 6; i++) {
2822 sigm[i] = 0.;
2823 h1[i] = 0;
2824 if (XVU_coor[chNum][i][iseg] > -900.) {
2825 h1[i] = 1;
2826 sigm[i] = (Cluster_coor[chNum][i][iseg] * dw) / sq12;
2827 sigm2_[i] = sigm[i] * sigm[i];
2828 } // if coord was
2829 // if (fDebug) cout<<" chNum "<<chNum<<" i "<<i<<" iseg "<<iseg<<" coor
2830 // "<<XVU_coor[chNum][i][iseg]<<endl; if (fDebug) cout<<" chNum "<<chNum<<" h "<<h1[i]<<" coor
2831 // "<<XVU_coor[chNum][i][iseg]<<" z_loc "<<z_loc[chNum][i]<<endl; if (fDebug) cout<<" h "<<h1[i]<<"
2832 // Cluster_coor["<<chNum<<"]["<<i<<"]["<<iseg<<"]= "<<Cluster_coor[chNum][i][iseg]<<" sigm
2833 // "<<sigm[i]<<endl;
2834 }
2835 if (fDebug)
2836 cout << endl;
2837
2838 for (Int_t im = 0; im < 4; im++) {
2839 for (Int_t ii = 0; ii < 4; ii++) {
2840 Amatr[im][ii] = 0.;
2841 bmatr[im][ii] = 0.;
2842 }
2843 }
2844 Double_t matrF[4] = {0, 0, 0, 0}; // free coef
2845
2846 FillFitMatrix(chNum, Amatr, z_loc, sigm2_, h1);
2847 FillFreeCoefVector(chNum, matrF, XVU_coor, iseg, z_loc, sigm2_, h1);
2848
2849 // Gaussian algorithm for 4x4 matrix inversion
2850 Double_t A0matr[4][4];
2851 for (Int_t i1 = 0; i1 < 4; i1++) {
2852 for (Int_t j1 = 0; j1 < 4; j1++) {
2853 A0matr[i1][j1] = Amatr[i1][j1];
2854 // if (fDebug) cout<<" Amatr["<<i1<<"]["<<j1<<"] "<<Amatr[i1][j1]<<endl;
2855 }
2856 }
2857
2858 InverseMatrix(Amatr, bmatr);
2859 Double_t sum;
2860 // Double_t A1[4][4] = {{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}};
2861 // cout<<" A1 "<<endl;
2862
2863 for (Int_t i1 = 0; i1 < 4; ++i1)
2864 for (Int_t j1 = 0; j1 < 4; ++j1) {
2865 sum = 0;
2866 for (Int_t k1 = 0; k1 < 4; ++k1) {
2867 Double_t a0 = A0matr[i1][k1];
2868 Double_t b0 = bmatr[k1][j1];
2869 sum += a0 * b0;
2870 // A1[i1][j1] = sum;
2871 }
2872 }
2873
2874 for (Int_t i1 = 0; i1 < 4; i1++) {
2875 par_ab[chNum][i1][iseg] = 0;
2876 for (Int_t j1 = 0; j1 < 4; j1++) {
2877 par_ab[chNum][i1][iseg] += bmatr[i1][j1] * matrF[j1];
2878 bmatr_seg[iseg][i1][j1] = 2 * bmatr[i1][j1];
2879 // if (fDebug) cout<<" i1 "<<i1<<" bmatr "<<bmatr[i1][j1]<<" F "<<matrF[j1] <<endl;
2880 }
2881 } // i1
2882
2883 //---Chi2 calculating---
2884 Chi2[chNum][iseg] = 0.;
2885 Chi2_ndf[chNum][iseg] = 999.;
2886 Double_t Max_dx = 0.;
2887
2888 //--- target area checking---
2889 x_target = par_ab[chNum][0][iseg] * (Z0_SRC - ChZ[chNum]) + par_ab[chNum][1][iseg];
2890 y_target = par_ab[chNum][2][iseg] * (Z0_SRC - ChZ[chNum]) + par_ab[chNum][3][iseg];
2891
2892 if (fabs(x_target) > xy_t_max[chNum] || fabs(y_target) > xy_t_max[chNum]) {
2893 Nhits_Ch_[chNum][iseg] = 0;
2894 continue; // go to next segment!
2895 }
2896
2897 if (fDebug)
2898 cout << " Nhits_Ch_[" << chNum << "][" << iseg << "] " << Nhits_Ch_[chNum][iseg] << endl;
2899
2900 for (Int_t i1 = 0; i1 < 6; i1++) {
2901 dx_[i1] = 0.;
2902
2903 if (XVU_coor[chNum][i1][iseg] > -999.) {
2904 if (i1 == 0 || i1 == 3) {
2905 dx_[i1] = XVU_coor[chNum][i1][iseg] - par_ab[chNum][0][iseg] * z_loc[chNum][i1]
2906 - par_ab[chNum][1][iseg];
2907 if (fabs(dx_[i1]) > Max_dx)
2908 Max_dx = fabs(dx_[i1]);
2909 // if (fDebug) cout<<" fitX
2910 // "<<par_ab[chNum][0][iseg]*z_loc[chNum][i1]+par_ab[chNum][1][iseg]<<endl;
2911 }
2912 if (fDebug && i1 == 0)
2913 hResiduals_pl0_Ch.at(chNum)->Fill(dx_[i1]);
2914 if (fDebug && i1 == 3)
2915 hResiduals_pl3_Ch.at(chNum)->Fill(dx_[i1]);
2916 if (i1 == 2 || i1 == 5) {
2917 dx_[i1] = XVU_coor[chNum][i1][iseg]
2918 - 0.5 * (par_ab[chNum][0][iseg] + sq3 * par_ab[chNum][2][iseg]) * z_loc[chNum][i1]
2919 - 0.5 * (par_ab[chNum][1][iseg] + sq3 * par_ab[chNum][3][iseg]);
2920 if (fabs(dx_[i1]) > Max_dx)
2921 Max_dx = fabs(dx_[i1]);
2922 // if (fDebug) cout<<" fitU
2923 // "<<0.5*(par_ab[chNum][0][iseg]+sq3*par_ab[chNum][2][iseg])*z_loc[chNum][i1]+0.5*(par_ab[chNum][1][iseg]+sq3*par_ab[chNum][3][iseg])<<endl;
2924 }
2925 if (fDebug && i1 == 2)
2926 hResiduals_pl2_Ch.at(chNum)->Fill(dx_[i1]);
2927 if (fDebug && i1 == 5)
2928 hResiduals_pl5_Ch.at(chNum)->Fill(dx_[i1]);
2929 if (i1 == 1 || i1 == 4) {
2930 dx_[i1] = XVU_coor[chNum][i1][iseg]
2931 - 0.5 * (par_ab[chNum][0][iseg] - sq3 * par_ab[chNum][2][iseg]) * z_loc[chNum][i1]
2932 - 0.5 * (par_ab[chNum][1][iseg] - sq3 * par_ab[chNum][3][iseg]);
2933 if (fabs(dx_[i1]) > Max_dx)
2934 Max_dx = fabs(dx_[i1]);
2935 // if (fDebug) cout<<" fitV
2936 // "<<0.5*(par_ab[chNum][0][iseg]-sq3*par_ab[chNum][2][iseg])*z_loc[chNum][i1]+0.5*(par_ab[chNum][1][iseg]-sq3*par_ab[chNum][3][iseg])<<endl;
2937 }
2938 if (fDebug && i1 == 1)
2939 hResiduals_pl1_Ch.at(chNum)->Fill(dx_[i1]);
2940 if (fDebug && i1 == 4)
2941 hResiduals_pl4_Ch.at(chNum)->Fill(dx_[i1]);
2942 Chi2[chNum][iseg] += dx_[i1] * dx_[i1] / (sigm2_[i1]);
2943 if (fDebug)
2944 cout << "iseg " << iseg << " i " << i1 << " dx " << dx_[i1] << " coor "
2945 << XVU_coor[chNum][i1][iseg] << " Chi2 " << Chi2[chNum][iseg] << " z " << z_loc[chNum][i1]
2946 << endl;
2947 }
2948 } // i1
2949 //---Chi2 calculating.---
2950 if (fDebug)
2951 cout << "---Chi2 calculating.--- Nhits[" << chNum << "][" << iseg << "] " << Nhits_Ch_[chNum][iseg]
2952 << " Chi2 " << Chi2[chNum][iseg] << endl;
2953 if (fDebug)
2954 cout << " Ax= " << par_ab[chNum][0][iseg] << " bx= " << par_ab[chNum][1][iseg]
2955 << " Ay= " << par_ab[chNum][2][iseg] << " by= " << par_ab[chNum][3][iseg] << endl;
2956 if (Nhits_Ch_[chNum][iseg] > 4) {
2957 Chi2_ndf[chNum][iseg] = Chi2[chNum][iseg] / (Nhits_Ch_[chNum][iseg] - 4);
2958 if (fDebug)
2959 cout << " Chi2_ndf[" << chNum << "][" << iseg << "] " << Chi2_ndf[chNum][iseg] << endl;
2960 } else {
2961 Chi2_ndf[chNum][iseg] = 0.;
2962 }
2963
2964 if (Chi2_ndf[chNum][iseg] > kChi2_Max_) { // --if bad chi2--
2965
2966 if (Nhits_Ch_[chNum][iseg] <= Min_hits6) { // --if no enough points--
2967 Nhits_Ch_[chNum][iseg] = 0;
2968 continue;
2969 } else { //--reject most distant point--
2970 if (fDebug)
2971 cout << " iseg " << iseg << " --reject or not bad point" << endl;
2972 for (Int_t i1 = 0; i1 < kNPlanes; i1++) {
2973 if (XVU_coor[chNum][i1][iseg] > -999. && fabs(dx_[i1]) >= Max_dx) {
2974 XVU_coor[chNum][i1][iseg] = -999.;
2975 Nhits_Ch_[chNum][iseg]--; // --reject bad point
2976 if (Nhits_Ch_[chNum][iseg] < Min_hits6)
2977 Nhits_Ch_[chNum][iseg] = 0;
2978 break; // continue;
2979 } // if point was
2980 } // i1
2981 } // else
2982
2983 } // if bad chi2.
2984
2985 if (fDebug)
2986 cout << "---Nhits calc--- Nhits[" << chNum << "][" << iseg << "] " << Nhits_Ch_[chNum][iseg]
2987 << " Chi2_ndf[" << chNum << "][" << iseg << "] " << Chi2_ndf[chNum][iseg] << endl;
2988
2989 if (Nhits_Ch_[chNum][iseg] != 0) {
2990 // if ( Chi2_ndf[chNum][iseg] > Chi2_Super_min && Chi2_ndf[chNum][iseg] < kChi2_Max_ &&
2991 // Nhits_Ch_[chNum][iseg] != 0){
2992 Nhits_Ch_[chNum][iseg] = 0;
2993 for (Int_t i1 = 0; i1 < kNPlanes; i1++) {
2994 if (XVU_coor[chNum][i1][iseg] > -900.)
2995 Nhits_Ch_[chNum][iseg]++;
2996 }
2997 if (fDebug)
2998 cout << "if ( Nhits_Ch_[chNum][iseg] != 0) Nhits_Ch " << Nhits_Ch_[chNum][iseg] << " Min_hits6 "
2999 << Min_hits6 << endl;
3000 if (Nhits_Ch_[chNum][iseg] < Min_hits6)
3001 Nhits_Ch_[chNum][iseg] = 0; // control shot
3002 }
3003
3004 } //--iseg--------------
3005
3006 if (fDebug)
3007 cout << endl;
3008 // if (fDebug ) cout<<" Nhitm "<<Nhitm<<endl;
3009
3010 for (Int_t iseg = 0; iseg < Nsegm[chNum]; iseg++) {
3011 x_target = par_ab[chNum][0][iseg] * (Z0_SRC - ChZ[chNum]) + par_ab[chNum][1][iseg];
3012 y_target = par_ab[chNum][2][iseg] * (Z0_SRC - ChZ[chNum]) + par_ab[chNum][3][iseg];
3013
3014 if (fDebug) {
3015 hx_target->Fill(x_target);
3016 hy_target->Fill(y_target);
3017 // cout<<Nhits_Ch_[chNum][iseg]<<" ";
3018 if (Nhits_Ch_[chNum][iseg] == Nhitm)
3019 cout << " chNum " << chNum << " iseg " << iseg << " Nhits_Ch " << Nhits_Ch_[chNum][iseg] << " Chi2 "
3020 << Chi2_ndf[chNum][iseg] << " x_target " << x_target << " y_target " << y_target << endl;
3021 // if (Nhits_Ch_[chNum][iseg] == Nhitm ) cout<<" ax "<<par_ab[chNum][0][iseg]<<" bx
3022 // "<<par_ab[chNum][1][iseg]<<" ay "<<par_ab[chNum][2][iseg]<<" by "<<par_ab[chNum][3][iseg]<<endl;
3023 }
3024 }
3025 if (fDebug)
3026 cout << endl;
3027 if (!ifNhitm)
3028 continue;
3029
3030 //--choice of min chi2
3031
3032 Int_t vec_ind_best_seg[kmaxSeg];
3033 Int_t Nbest_vec_ind = 0;
3034
3035 for (Int_t ind_best = 0; ind_best < kmaxSeg; ind_best++) {
3036 Double_t Chi2_best = 999.;
3037 Int_t iseg_best = -1;
3038 vec_ind_best_seg[ind_best] = -1;
3039
3040 for (Int_t iseg = 0; iseg < Nsegm[chNum]; iseg++) {
3041
3042 if (Nhits_Ch_[chNum][iseg] != Nhitm)
3043 continue;
3044 if (Nhitm > 4 && Chi2_ndf[chNum][iseg] >= Chi2_best)
3045 continue;
3046
3047 Bool_t best_was = 0;
3048 if (Nbest_vec_ind > 0) {
3049 for (Int_t ii = 0; ii < Nbest_vec_ind; ii++) {
3050 if (iseg == vec_ind_best_seg[ii]) {
3051 best_was = 1;
3052 break;
3053 }
3054 }
3055 }
3056 if (best_was)
3057 continue;
3058 Chi2_best = Chi2_ndf[chNum][iseg];
3059 iseg_best = iseg;
3060 } // iseg
3061
3062 if (iseg_best == -1)
3063 continue;
3064 if (fDebug)
3065 cout << " ind_best " << ind_best << " Chi2_best " << Chi2_best << " iseg_best " << iseg_best << endl;
3066 vec_ind_best_seg[ind_best] = iseg_best;
3067 Nbest_vec_ind++;
3068
3069 if (fDebug)
3070 cout << "-- reject(common points)" << endl;
3071 for (Int_t iseg = 0; iseg < Nsegm[chNum]; iseg++) {
3072 // if (iseg == iseg_best)continue;
3073 Bool_t best_was = 0;
3074 if (Nbest_vec_ind > 0) {
3075 for (Int_t ii = 0; ii < Nbest_vec_ind; ii++) {
3076 if (iseg == vec_ind_best_seg[ii]) {
3077 best_was = 1;
3078 break;
3079 }
3080 }
3081 }
3082 if (best_was)
3083 continue;
3084
3085 for (Int_t i1 = 0; i1 < kNPlanes; i1++) {
3086 if (XVU_coor[chNum][i1][iseg] > -999.) {
3087 if (fabs(XVU_coor[chNum][i1][iseg] - XVU_coor[chNum][i1][iseg_best]) < dw_half) { // 3*dw_half
3088 if (fDebug)
3089 cout << " reject point pl:" << i1 << " XUV " << XVU_coor[chNum][i1][iseg]
3090 << " iseg = " << iseg << endl;
3091 // Nhits_Ch_[chNum][iseg] = 0;
3092 XVU_coor[chNum][i1][iseg] = -999.;
3093 Nhits_Ch_[chNum][iseg]--; // --reject bad point
3094 if (Nhits_Ch_[chNum][iseg] < Min_hits6) {
3095 Nhits_Ch_[chNum][iseg] = 0;
3096 break;
3097 }
3098 }
3099 }
3100 } // i1
3101 } // iseg
3102 if (fDebug)
3103 cout << " Nbest_vec_ind " << Nbest_vec_ind << endl;
3104 if (fDebug)
3105 cout << " Chi2_best " << Chi2_best << " Nhits " << Nhits_Ch_[chNum][iseg_best] << endl;
3106 } // ind_best
3107
3108 vector<segments> vtmpSeg;
3109 segments tmpSeg;
3110
3111 for (int itSeg = 0; itSeg < Nsegm[chNum]; ++itSeg) {
3112 if (Nhits_Ch_[chNum][itSeg] == Nhitm && Chi2_ndf[chNum][itSeg] < kChi2_Max_) {
3113 tmpSeg.Nhits = Nhits_Ch_[chNum][itSeg];
3114 tmpSeg.Chi2 = Chi2_ndf[chNum][itSeg];
3115
3116 for (int ipla = 0; ipla < kNPlanes; ipla++) {
3117 tmpSeg.coord[ipla] = XVU_coor[chNum][ipla][itSeg];
3118 tmpSeg.clust[ipla] = Cluster_coor[chNum][ipla][itSeg];
3119 }
3120
3121 Double_t x_t = par_ab[chNum][0][itSeg] * (Z0_SRC - ChZ[chNum]) + par_ab[chNum][1][itSeg];
3122 // Double_t y_t = par_ab[chNum][2][itSeg]*( Z0_SRC - ChZ[chNum]) + par_ab[chNum][3][itSeg];
3123 if (fDebug) {
3124 hx_target_best.at(chNum)->Fill(x_t);
3125 hy_target_best.at(chNum)->Fill(x_t);
3126 }
3127
3128 for (int ipar = 0; ipar < 4; ipar++) {
3129 tmpSeg.param[ipar] = par_ab[chNum][ipar][itSeg];
3130 }
3131 for (Int_t i1 = 0; i1 < 4; i1++) {
3132 for (Int_t j1 = 0; j1 < 4; j1++) {
3133 tmpSeg.sigma2[i1][j1] = bmatr_seg[itSeg][i1][j1];
3134 }
3135 }
3136 vtmpSeg.push_back(tmpSeg);
3137 }
3138 }
3139
3140 // vector sorting
3141 if (Nhitm > 4)
3142 sort(vtmpSeg.begin(), vtmpSeg.end(), compareSegments);
3143
3144 // storing
3145 for (size_t iterOut = 0; iterOut < vtmpSeg.size(); iterOut++) {
3146 if (OutSegCount < kmaxSeg) {
3147 OutSegArray[OutSegCount] = vtmpSeg.at(iterOut);
3148 OutSegCount++;
3149 }
3150 }
3151 // !!! vector clear for next Nhitm
3152 // if ( Nhitm == 5)
3153 vtmpSeg.clear(); // clear for 4p-segments
3154
3155 // if (fDebug) cout<<" ---------Nhitm "<<Nhitm<<endl;
3156 // for (Int_t iseg = 0; iseg < Nsegm[chNum]; iseg++) { //---print number of hits
3157 // if (fDebug) cout<<" iseg "<<iseg<<" Nhits "<<Nhits_Ch_[chNum][iseg]<<endl;
3158 // }
3159 // if (fDebug) cout<<" ---------"<<endl;
3160 } // Nhitm--
3161
3162 for (int iterOut = 0; iterOut < kmaxSeg; iterOut++) {
3163 if (OutSegArray[iterOut].Chi2 < 900.) {
3164 Nbest_seg_[chNum]++;
3165 Nhits_seg_[chNum][iterOut] = OutSegArray[iterOut].Nhits;
3166 Chi2_ndf_seg_[chNum][iterOut] = OutSegArray[iterOut].Chi2;
3167
3168 for (int iterCoord = 0; iterCoord < kNPlanes; iterCoord++) {
3169 coor_seg[chNum][iterCoord][iterOut] = OutSegArray[iterOut].coord[iterCoord];
3170 cluster_seg[chNum][iterCoord][iterOut] = OutSegArray[iterOut].clust[iterCoord];
3171 }
3172 for (int ipar = 0; ipar < 4; ipar++) {
3173 Par_ab_seg[chNum][ipar][iterOut] = OutSegArray[iterOut].param[ipar];
3174 }
3175 for (int i1 = 0; i1 < 4; i1++) {
3176 for (int j1 = 0; j1 < 4; j1++) {
3177 sigma2_s[chNum][iterOut][i1][j1] = OutSegArray[iterOut].sigma2[i1][j1];
3178 }
3179 }
3180 }
3181 }
3182
3183 Nsegm[chNum] = 0; // go to the next case for 2coord
3184 for (Int_t i1 = 0; i1 < kNPlanes; i1++) {
3185 for (Int_t iBig = 0; iBig < kBig; iBig++) {
3186 XVU_coor[chNum][i1][iBig] = -999.;
3187 Cluster_coor[chNum][i1][iBig] = -1;
3188 }
3189 }
3190
3191 // printf(">>>TESTING\n");
3192 // for (int iterOut = 0; iterOut < kmaxSeg; iterOut++){
3193 // if(OutSegArray[iterOut].Chi2 < 900.){
3194 // printf("> Ch: %d, Seg %d, Nhits: %d, Chi2: %8.4f\n", chNum, iterOut, OutSegArray[iterOut].Nhits,
3195 // OutSegArray[iterOut].Chi2);
3196 // for (int iterCoord = 0; iterCoord < kNPlanes; iterCoord++)
3197 // printf("> %8.4f - %8.4f\n", OutSegArray[iterOut].coord[iterCoord], OutSegArray[iterOut].clust[iterCoord]);
3198 // }
3199 // }
3200
3201} // ProcessSegments
3202//----------------------------------------------------------------------
3203
3204// --------------local parameters to Global parameters------------------
3205void BmnMwpcHitFinder::SegmentParamAlignment(Int_t chNum, Int_t* Nbest, Double_t*** par_ab, Float_t** shiftt)
3206{
3207 if (fDebug)
3208 cout << endl;
3209 if (fDebug)
3210 cout << " SegmentParamAlignment " << endl;
3211
3212 for (Int_t iBest = 0; iBest < Nbest[chNum]; iBest++) {
3213 // cout<<" chNum "<<chNum<<" before Alignment: iBest "<<iBest<<" Ax "<< par_ab[chNum][0][iBest]<<" bx "<<
3214 // par_ab[chNum][1][iBest]<<" Ay "<< par_ab[chNum][1][iBest]<<" by "<< par_ab[chNum][3][iBest]<<endl;
3215 // ax alpha ax^2
3216 par_ab[chNum][0][iBest] +=
3217 shiftt[chNum][0] + shiftt[chNum][0] * par_ab[chNum][0][iBest] * par_ab[chNum][0][iBest];
3218 par_ab[chNum][2][iBest] +=
3219 shiftt[chNum][2] + shiftt[chNum][2] * par_ab[chNum][2][iBest] * par_ab[chNum][2][iBest];
3220 par_ab[chNum][1][iBest] += shiftt[chNum][1];
3221 par_ab[chNum][3][iBest] += shiftt[chNum][3];
3222 if (fDebug)
3223 cout << " ax " << par_ab[chNum][0][iBest] << " bx " << par_ab[chNum][1][iBest] << " ay "
3224 << par_ab[chNum][2][iBest] << " by " << par_ab[chNum][3][iBest] << endl;
3225
3226 } // iBest
3227} // SegmentParamAlignment
3228//----------------------------------------------------------------------
3229
3230//------ Arrays Initialization -----------------------------------------
3231void BmnMwpcHitFinder::PrepareArraysToProcessEvent()
3232{
3233 fBmnMwpcSegmentsArray->Delete();
3234 vec_points.clear();
3235
3236 for (Int_t iCh = 0; iCh < kNChambers; iCh++) {
3237 Nseg_Ch[iCh] = 0;
3238 Nbest_Ch[iCh] = 0;
3239 Nbest_seg[iCh] = 0;
3240 Nlay_w_wires[iCh] = 0;
3241
3242 for (Int_t iPl = 0; iPl < kNPlanes; iPl++) {
3243 iw_Ch[iCh][iPl] = 0;
3244 XVU[iCh][iPl] = 0;
3245 XVU_cl[iCh][iPl] = 0;
3246 Nclust[iCh][iPl] = 0;
3247
3248 for (Int_t iWire = 0; iWire < kNWires; iWire++) {
3249 DigitsArray[iCh][iPl][iWire] = 0.;
3250 }
3251
3252 for (Int_t iBig = 0; iBig < kBig; iBig++) {
3253 Wires_Ch[iCh][iPl][iBig] = -1;
3254 clust_Ch[iCh][iPl][iBig] = -1;
3255 XVU_Ch[iCh][iPl][iBig] = -999.;
3256 Coord_wire[iCh][iPl][iBig] = -999.;
3257 Coord_xuv[iCh][iPl][iBig] = -999.;
3258 ClusterSize[iCh][iPl][iBig] = 0;
3259 XVU_coord[iCh][iPl][iBig] = -999.;
3260 Cluster_coord[iCh][iPl][iBig] = -1;
3261 Coor_seg[iCh][iPl][iBig] = -999.;
3262 Cluster_seg[iCh][iPl][iBig] = -1;
3263 }
3264 for (Int_t iWire = 0; iWire < kNWires; iWire++) {
3265 wire_Ch[iCh][iWire][iPl] = 0;
3266 xuv_Ch[iCh][iWire][iPl] = 0.;
3267 }
3268 }
3269
3270 for (Int_t ii = 0; ii < 4; ii++) {
3271 for (Int_t jj = 0; jj < kBig; jj++) {
3272 par_ab_Ch[iCh][ii][jj] = 999.;
3273 par_ab_seg[iCh][ii][jj] = 999.;
3274 for (Int_t j = 0; j < 4; j++) {
3275 sigma2_seg[iCh][jj][ii][j] = 0;
3276 }
3277 }
3278 }
3279
3280 for (Int_t iBig = 0; iBig < kBig; iBig++) {
3281 Nhits_Ch[iCh][iBig] = 0;
3282 Nhits_seg[iCh][iBig] = 0;
3283 Chi2_ndf_Ch[iCh][iBig] = 0;
3284 }
3285
3286 for (Int_t i = 0; i < kmaxSeg; i++) {
3287 ind_best_Ch[iCh][i] = 0;
3288 best_Ch_gl[iCh][i] = -1;
3289 Chi2_ndf_best_Ch[iCh][i] = -999.;
3290 Chi2_ndf_seg[iCh][i] = -999.;
3291 }
3292 } // iCh
3293
3294 for (Int_t iPl = 0; iPl < kNPlanes; iPl++) {
3295 sigm2[iPl] = sigma * sigma;
3296 ipl[iPl] = 6;
3297 z2[iPl] = 0;
3298 }
3299
3300 for (Int_t ii = 0; ii < 4; ii++) {
3301 for (Int_t jj = 0; jj < 4; jj++) {
3302 matrA[ii][jj] = 0.;
3303 matrb[ii][jj] = 0.;
3304 }
3305 }
3306} // PrepareArraysToProcessEvent
3307//----------------------------------------------------------------------
3308
3309//----------------------Fill FitMatrix----------------------------------
3310void BmnMwpcHitFinder::FillFitMatrix(Int_t chN, Double_t** AA, Float_t** z, Double_t* sigm2_, Int_t* h_)
3311{
3312 // AA - matrix to be filledlayers)
3313 // sigm2 - square of sigma
3314 // h_ - array to include/exclude planes (h_[i] = 0 or 1)
3315 // Float_t z2_[nPlanes];
3316 Float_t z2_[6] = {z[chN][0] * z[chN][0], z[chN][1] * z[chN][1], z[chN][2] * z[chN][2],
3317 z[chN][3] * z[chN][3], z[chN][4] * z[chN][4], z[chN][5] * z[chN][5]}; // cm
3318
3319 AA[0][0] += 2 * z2_[0] * h_[0] / sigm2_[0] + z2_[2] * h_[2] / (2 * sigm2_[2]) + z2_[1] * h_[1] / (2 * sigm2_[1])
3320 + 2 * z2_[3] * h_[3] / sigm2_[3] + z2_[5] * h_[5] / (2 * sigm2_[5])
3321 + z2_[4] * h_[4] / (2 * sigm2_[4]); // Ax
3322 AA[0][1] += 2 * z[chN][0] * h_[0] / sigm2_[0] + z[chN][2] * h_[2] / (2 * sigm2_[2])
3323 + z[chN][1] * h_[1] / (2 * sigm2_[1]) + 2 * z[chN][3] * h_[3] / sigm2_[3]
3324 + z[chN][5] * h_[5] / (2 * sigm2_[5]) + z[chN][4] * h_[4] / (2 * sigm2_[4]); // Bx
3325 AA[0][2] += sq3
3326 * (z2_[2] * h_[2] / (2 * sigm2_[2]) - z2_[1] * h_[1] / (2 * sigm2_[1])
3327 + z2_[5] * h_[5] / (2 * sigm2_[5]) - z2_[4] * h_[4] / (2 * sigm2_[4])); // Ay
3328 AA[0][3] += sq3
3329 * (z[chN][2] * h_[2] / (2 * sigm2_[2]) - z[chN][1] * h_[1] / (2 * sigm2_[1])
3330 + z[chN][5] * h_[5] / (2 * sigm2_[5]) - z[chN][4] * h_[4] / (2 * sigm2_[4])); // By
3331 AA[1][0] = AA[0][1];
3332 AA[1][1] += 2 * h_[0] / sigm2_[0] + 0.5 * h_[2] / sigm2_[2] + 0.5 * h_[1] / sigm2_[1] + 2 * h_[3] / sigm2_[3]
3333 + 0.5 * h_[5] / sigm2_[5] + 0.5 * h_[4] / sigm2_[4];
3334 AA[1][2] += sq3
3335 * (z[chN][2] * h_[2] / sigm2_[2] - z[chN][1] * h_[1] / sigm2_[1] + z[chN][5] * h_[5] / sigm2_[5]
3336 - z[chN][4] * h_[4] / sigm2_[4])
3337 * 0.5;
3338 AA[1][3] += sq3 * (h_[2] / sigm2_[2] - h_[1] / sigm2_[1] + h_[5] / sigm2_[5] - h_[4] / sigm2_[4]) * 0.5;
3339 AA[2][0] = AA[0][2];
3340 AA[2][1] = AA[1][2];
3341 AA[2][2] += 3.0
3342 * (z2_[2] * h_[2] / sigm2_[2] + z2_[1] * h_[1] / sigm2_[1] + z2_[5] * h_[5] / sigm2_[5]
3343 + z2_[4] * h_[4] / sigm2_[4])
3344 * 0.5;
3345 AA[2][3] += 3.0
3346 * (z[chN][2] * h_[2] / sigm2_[2] + z[chN][1] * h_[1] / sigm2_[1] + z[chN][5] * h_[5] / sigm2_[5]
3347 + z[chN][4] * h_[4] / sigm2_[4])
3348 * 0.5;
3349 AA[3][0] = AA[0][3];
3350 AA[3][1] = AA[1][3];
3351 AA[3][2] = AA[2][3];
3352 AA[3][3] +=
3353 3.0 * (0.5 * h_[2] / sigm2_[2] + 0.5 * h_[1] / sigm2_[1] + 0.5 * h_[5] / sigm2_[5] + 0.5 * h_[4] / sigm2_[4]);
3354}
3355//----------------------------------------------------------------------
3356
3357//--------------------Matrix Free Coefficients Calculation--------------
3358void BmnMwpcHitFinder::FillFreeCoefVector(Int_t ichNum,
3359 Double_t* F,
3360 Double_t*** XVU_,
3361 Int_t ise,
3362 Float_t** z,
3363 Double_t* sigmm2,
3364 Int_t* h_)
3365{
3366 // F - vector to be filled
3367 // XVU_ - coordinates of segment in chamber (Is it correct definition?)
3368 // segIdx - index of current segment
3369 // z - local z-positions of planes(layers)
3370 // sigmm2 - square of sigma
3371 // h_ - array to include/exclude planes (h_[i] = 0 or 1)
3372
3373 F[0] += 2 * XVU_[ichNum][0][ise] * z[ichNum][0] * h_[0] / sigmm2[0]
3374 + XVU_[ichNum][1][ise] * z[ichNum][1] * h_[1] / sigmm2[1]
3375 + XVU_[ichNum][2][ise] * z[ichNum][2] * h_[2] / sigmm2[2]
3376 + 2 * XVU_[ichNum][3][ise] * z[ichNum][3] * h_[3] / sigmm2[3]
3377 + XVU_[ichNum][4][ise] * z[ichNum][4] * h_[4] / sigmm2[4]
3378 + XVU_[ichNum][5][ise] * z[ichNum][5] * h_[5] / sigmm2[5];
3379 F[1] += 2 * XVU_[ichNum][0][ise] * h_[0] / sigmm2[0] + XVU_[ichNum][1][ise] * h_[1] / sigmm2[1]
3380 + XVU_[ichNum][2][ise] * h_[2] / sigmm2[2] + 2 * XVU_[ichNum][3][ise] * h_[3] / sigmm2[3]
3381 + XVU_[ichNum][4][ise] * h_[4] / sigmm2[4] + XVU_[ichNum][5][ise] * h_[5] / sigmm2[5];
3382 F[2] += sq3
3383 * (-XVU_[ichNum][1][ise] * z[ichNum][1] * h_[1] / sigmm2[1]
3384 + XVU_[ichNum][2][ise] * z[ichNum][2] * h_[2] / sigmm2[2]
3385 - XVU_[ichNum][4][ise] * z[ichNum][4] * h_[4] / sigmm2[4]
3386 + XVU_[ichNum][5][ise] * z[ichNum][5] * h_[5] / sigmm2[5]);
3387 F[3] += sq3
3388 * (-XVU_[ichNum][1][ise] * h_[1] / sigmm2[1] + XVU_[ichNum][2][ise] * h_[2] / sigmm2[2]
3389 - XVU_[ichNum][4][ise] * h_[4] / sigmm2[4] + XVU_[ichNum][5][ise] * h_[5] / sigmm2[5]);
3390}
3391//----------------------------------------------------------------------
3392
3393//--------------------Matrix Coefficients Calculation--------------
3394void BmnMwpcHitFinder::FillFreeCoefVectorXUV(Int_t ichNum,
3395 Double_t* F,
3396 Float_t** XVU_,
3397 Float_t** z,
3398 Float_t* sigm2_,
3399 Int_t* h_)
3400{
3401 // F - vector to be filled
3402 // XVU_ - coordinates of segment in chamber (Is it correct definition?)
3403 // segIdx - index of current segment
3404 // z - local z-positions of planes(layers)
3405 // sigm2_ - square of sigma
3406 // h_ - array to include/exclude planes (h_[i] = 0 or 1)
3407 F[0] += 2 * XVU_[ichNum][0] * z[ichNum][0] * h_[0] / sigm2_[0] + XVU_[ichNum][1] * z[ichNum][1] * h_[1] / sigm2_[1]
3408 + XVU_[ichNum][2] * z[ichNum][2] * h_[2] / sigm2_[2]
3409 + 2 * XVU_[ichNum][3] * z[ichNum][3] * h_[3] / sigm2_[3]
3410 + XVU_[ichNum][4] * z[ichNum][4] * h_[4] / sigm2_[4] + XVU_[ichNum][5] * z[ichNum][5] * h_[5] / sigm2_[5];
3411 F[1] += 2 * XVU_[ichNum][0] * h_[0] / sigm2_[0] + XVU_[ichNum][1] * h_[1] / sigm2_[1]
3412 + XVU_[ichNum][2] * h_[2] / sigm2_[2] + 2 * XVU_[ichNum][3] * h_[3] / sigm2_[3]
3413 + XVU_[ichNum][4] * h_[4] / sigm2_[4] + XVU_[ichNum][5] * h_[5] / sigm2_[5];
3414 F[2] += (-XVU_[ichNum][1] * z[ichNum][1] * h_[1] / sigm2_[1] + XVU_[ichNum][2] * z[ichNum][2] * h_[2] / sigm2_[2]
3415 - XVU_[ichNum][4] * z[ichNum][4] * h_[4] / sigm2_[4] + XVU_[ichNum][5] * z[ichNum][5] * h_[5] / sigm2_[5]);
3416 F[3] += (-XVU_[ichNum][1] * h_[1] / sigm2_[1] + XVU_[ichNum][2] * h_[2] / sigm2_[2]
3417 - XVU_[ichNum][4] * h_[4] / sigm2_[4] + XVU_[ichNum][5] * h_[5] / sigm2_[5]);
3418 F[2] = F[2] * sq3;
3419 F[3] = F[3] * sq3;
3420}
3421//----------------------------------------------------------------------
3422
3423//--------------------Matrix Coefficients Calculation--------------
3424void BmnMwpcHitFinder::InverseMatrix(Double_t** AA, Double_t** bb)
3425{
3426 // Gaussian algorithm for 4x4 matrix inversion
3427 Double_t factor;
3428 Double_t temp[4];
3429 // Set b to I
3430 for (Int_t i1 = 0; i1 < 4; i1++) {
3431 for (Int_t j1 = 0; j1 < 4; j1++) {
3432 if (i1 == j1)
3433 bb[i1][j1] = 1.0;
3434 else
3435 bb[i1][j1] = 0.0;
3436 }
3437 }
3438 for (Int_t i1 = 0; i1 < 4; i1++) {
3439 for (Int_t j1 = i1 + 1; j1 < 4; j1++) {
3440 if (fabs(AA[i1][i1]) < fabs(AA[j1][i1])) {
3441 for (Int_t l1 = 0; l1 < 4; l1++)
3442 temp[l1] = AA[i1][l1];
3443 for (Int_t l1 = 0; l1 < 4; l1++)
3444 AA[i1][l1] = AA[j1][l1];
3445 for (Int_t l1 = 0; l1 < 4; l1++)
3446 AA[j1][l1] = temp[l1];
3447 for (Int_t l1 = 0; l1 < 4; l1++)
3448 temp[l1] = bb[i1][l1];
3449 for (Int_t l1 = 0; l1 < 4; l1++)
3450 bb[i1][l1] = bb[j1][l1];
3451 for (Int_t l1 = 0; l1 < 4; l1++)
3452 bb[j1][l1] = temp[l1];
3453 }
3454 }
3455 factor = AA[i1][i1];
3456 for (Int_t j1 = 4 - 1; j1 > -1; j1--) {
3457 bb[i1][j1] /= factor;
3458 AA[i1][j1] /= factor;
3459 }
3460 for (Int_t j1 = i1 + 1; j1 < 4; j1++) {
3461 factor = -AA[j1][i1];
3462 for (Int_t k1 = 0; k1 < 4; k1++) {
3463 AA[j1][k1] += AA[i1][k1] * factor;
3464 bb[j1][k1] += bb[i1][k1] * factor;
3465 }
3466 }
3467 } // i1
3468 for (Int_t i1 = 3; i1 > 0; i1--) {
3469 for (Int_t j1 = i1 - 1; j1 > -1; j1--) {
3470 factor = -AA[j1][i1];
3471 for (Int_t k1 = 0; k1 < 4; k1++) {
3472 AA[j1][k1] += AA[i1][k1] * factor;
3473 bb[j1][k1] += bb[i1][k1] * factor;
3474 }
3475 }
3476 } // i1
3477} // end inverse
3478//----------------------------------------------------------------------
3479
3480//----------------------------------------------------------------------
3481void BmnMwpcHitFinder::Finish()
3482{
3483
3484 if (fDebug) {
3485
3486 hEff_mc->Divide(hNum_mc, hDen_mc, 1, 1);
3487 hEff_mcreaction2->Divide(hNum_mcreaction2, hDen_mcreaction2, 1, 1);
3488 hEff_mcreaction3->Divide(hNum_mcreaction3, hDen_mcreaction3, 1, 1);
3489
3490 printf("MWPC hit finder: write hists to file... ");
3491 fOutputFileName = Form("hMWPChits_p%d_run%d.root", fRunPeriod, fRunNumber);
3492 cout << fOutputFileName << endl;
3493 TFile file(fOutputFileName, "RECREATE");
3494 fList.Write();
3495 file.Close();
3496 }
3497
3498 if (fDebug)
3499 printf("done\n");
3500 delete fMwpcGeometrySRC;
3501
3502 cout << "Work time of the MWPC hit finder: " << workTime << " s" << endl;
3503}
3504//----------------------------------------------------------------------
compareSegments
bool compareSegments(const segments &a, const segments &b)
Definition BmnMwpcHitFinder.cxx:44
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
v
__m128 v
Definition P4_F32vec4.h:1
i
int i
Definition P4_F32vec4.h:22
Bool_t
bool Bool_t
Int_t
int Int_t
Float_t
float Float_t
Short_t
short Short_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,...)
BmnMwpcDigit::GetPlane
Short_t GetPlane() const
Definition BmnMwpcDigit.h:26
BmnMwpcDigit::GetStation
Short_t GetStation() const
Definition BmnMwpcDigit.h:25
BmnMwpcDigit::GetTime
UInt_t GetTime() const
Definition BmnMwpcDigit.h:28
BmnMwpcDigit::GetWireNumber
Int_t GetWireNumber() const
Definition BmnMwpcDigit.h:27
BmnMwpcHitFinder::hBy_mc_ch
vector< TH1D * > hBy_mc_ch
Definition BmnMwpcHitFinder.h:127
BmnMwpcHitFinder::hWiresUp
vector< TH1D * > hWiresUp
Definition BmnMwpcHitFinder.h:124
BmnMwpcHitFinder::Cluster_coord
Double_t *** Cluster_coord
Definition BmnMwpcHitFinder.h:186
BmnMwpcHitFinder::hAx_mc_ch
vector< TH1D * > hAx_mc_ch
Definition BmnMwpcHitFinder.h:127
BmnMwpcHitFinder::~BmnMwpcHitFinder
virtual ~BmnMwpcHitFinder()
Destructor.
Definition BmnMwpcHitFinder.cxx:84
BmnMwpcHitFinder::par_ab_seg
Double_t *** par_ab_seg
Definition BmnMwpcHitFinder.h:184
BmnMwpcHitFinder::Coord_wire
Double_t *** Coord_wire
Definition BmnMwpcHitFinder.h:174
BmnMwpcHitFinder::XVU_coord
Double_t *** XVU_coord
Definition BmnMwpcHitFinder.h:173
BmnMwpcHitFinder::hY_X_Ch
vector< TH2D * > hY_X_Ch
Definition BmnMwpcHitFinder.h:131
BmnMwpcHitFinder::fBmnMwpcSegmentsArray
TClonesArray * fBmnMwpcSegmentsArray
Output array of MWPC hits.
Definition BmnMwpcHitFinder.h:96
BmnMwpcHitFinder::hWiresVm
vector< TH1D * > hWiresVm
Definition BmnMwpcHitFinder.h:125
BmnMwpcHitFinder::hoccupancyUm
vector< TH1D * > hoccupancyUm
Definition BmnMwpcHitFinder.h:124
BmnMwpcHitFinder::Init
virtual InitStatus Init()
Virtual method Init.
Definition BmnMwpcHitFinder.cxx:214
BmnMwpcHitFinder::hdY_mc_SegCh
vector< TH1D * > hdY_mc_SegCh
Definition BmnMwpcHitFinder.h:128
BmnMwpcHitFinder::hAy_mc_ch
vector< TH1D * > hAy_mc_ch
Definition BmnMwpcHitFinder.h:127
BmnMwpcHitFinder::hWiresXp
vector< TH1D * > hWiresXp
Definition BmnMwpcHitFinder.h:124
BmnMwpcHitFinder::hResiduals_pl5_Ch
vector< TH1D * > hResiduals_pl5_Ch
Definition BmnMwpcHitFinder.h:126
BmnMwpcHitFinder::hWiresVp
vector< TH1D * > hWiresVp
Definition BmnMwpcHitFinder.h:125
BmnMwpcHitFinder::hBx_mc_ch
vector< TH1D * > hBx_mc_ch
Definition BmnMwpcHitFinder.h:127
BmnMwpcHitFinder::Clustering
void Clustering(Int_t, Int_t ***, Double_t ***, Double_t ***, Double_t ***, Int_t **, Int_t *)
Definition BmnMwpcHitFinder.cxx:1560
BmnMwpcHitFinder::hdAy_mc_Seg_deltaCh
vector< TH1D * > hdAy_mc_Seg_deltaCh
Definition BmnMwpcHitFinder.h:128
BmnMwpcHitFinder::hResiduals_pl3_Ch
vector< TH1D * > hResiduals_pl3_Ch
Definition BmnMwpcHitFinder.h:126
BmnMwpcHitFinder::fMwpcGeometrySRC
BmnMwpcGeometrySRC * fMwpcGeometrySRC
Definition BmnMwpcHitFinder.h:104
BmnMwpcHitFinder::hoccupancyVm
vector< TH1D * > hoccupancyVm
Definition BmnMwpcHitFinder.h:124
BmnMwpcHitFinder::SegmentFinder2coord
void SegmentFinder2coord(Int_t, Int_t **, Double_t ***, Int_t ***, Int_t *, Double_t ***, Double_t ***, Int_t **, Int_t, Short_t, Int_t, Int_t *, Double_t ***)
Definition BmnMwpcHitFinder.cxx:2218
BmnMwpcHitFinder::hResiduals_pl1_Ch
vector< TH1D * > hResiduals_pl1_Ch
Definition BmnMwpcHitFinder.h:125
BmnMwpcHitFinder::Coord_xuv
Double_t *** Coord_xuv
Definition BmnMwpcHitFinder.h:172
BmnMwpcHitFinder::par_ab_Ch
Double_t *** par_ab_Ch
Definition BmnMwpcHitFinder.h:171
BmnMwpcHitFinder::fBmnHitsArray
TClonesArray * fBmnHitsArray
Definition BmnMwpcHitFinder.h:94
BmnMwpcHitFinder::hResiduals_pl2_Ch
vector< TH1D * > hResiduals_pl2_Ch
Definition BmnMwpcHitFinder.h:125
BmnMwpcHitFinder::hChisq_ndf_Ch
vector< TH1D * > hChisq_ndf_Ch
Definition BmnMwpcHitFinder.h:124
BmnMwpcHitFinder::Coor_seg
Double_t *** Coor_seg
Definition BmnMwpcHitFinder.h:182
BmnMwpcHitFinder::hdAx_mc_SegCh
vector< TH1D * > hdAx_mc_SegCh
Definition BmnMwpcHitFinder.h:128
BmnMwpcHitFinder::ProcessSegments2coord
void ProcessSegments2coord(Int_t, Int_t *, Double_t ***, Double_t ***, Int_t **, Float_t **, Int_t, Double_t, Double_t, Int_t **, Double_t **, Double_t ***, Double_t ***, Double_t ***, Int_t *, Int_t *, Double_t ****)
Definition BmnMwpcHitFinder.cxx:2447
BmnMwpcHitFinder::Chi2_ndf_Ch
Double_t ** Chi2_ndf_Ch
Definition BmnMwpcHitFinder.h:164
BmnMwpcHitFinder::hdX_mc_SegCh
vector< TH1D * > hdX_mc_SegCh
Definition BmnMwpcHitFinder.h:128
BmnMwpcHitFinder::htime_wire_Ch
vector< TH2D * > htime_wire_Ch
Definition BmnMwpcHitFinder.h:131
BmnMwpcHitFinder::BmnMwpcHitFinder
BmnMwpcHitFinder()
Default constructor.
Definition BmnMwpcHitFinder.h:43
BmnMwpcHitFinder::hx_target_best
vector< TH1D * > hx_target_best
Definition BmnMwpcHitFinder.h:127
BmnMwpcHitFinder::DigitsArray
Double_t *** DigitsArray
Definition BmnMwpcHitFinder.h:180
BmnMwpcHitFinder::hy_target_best
vector< TH1D * > hy_target_best
Definition BmnMwpcHitFinder.h:127
BmnMwpcHitFinder::Exec
virtual void Exec(Option_t *opt)
Virtual method Exec.
Definition BmnMwpcHitFinder.cxx:819
BmnMwpcHitFinder::hdY_mc_Seg_deltaCh
vector< TH1D * > hdY_mc_Seg_deltaCh
Definition BmnMwpcHitFinder.h:129
BmnMwpcHitFinder::hdAx_mc_Seg_deltaCh
vector< TH1D * > hdAx_mc_Seg_deltaCh
Definition BmnMwpcHitFinder.h:128
BmnMwpcHitFinder::hNum_layers_out_beam_Ch
vector< TH1D * > hNum_layers_out_beam_Ch
Definition BmnMwpcHitFinder.h:126
BmnMwpcHitFinder::hResiduals_pl4_Ch
vector< TH1D * > hResiduals_pl4_Ch
Definition BmnMwpcHitFinder.h:126
BmnMwpcHitFinder::ReadWires
void ReadWires(Double_t ***, Int_t **, vector< MC_points > &)
Definition BmnMwpcHitFinder.cxx:1247
BmnMwpcHitFinder::ProcessSegments
void ProcessSegments(Int_t, Int_t *, Double_t ***, Double_t ***, Int_t **, Float_t **, Int_t, Double_t, Double_t, Int_t **, Double_t **, Double_t ***, Double_t ***, Double_t ***, Int_t *, Int_t *, Double_t ****)
Definition BmnMwpcHitFinder.cxx:2723
BmnMwpcHitFinder::hpar_Bx_Ch
vector< TH1D * > hpar_Bx_Ch
Definition BmnMwpcHitFinder.h:127
BmnMwpcHitFinder::hfiredWire_Ch
vector< TH1D * > hfiredWire_Ch
Definition BmnMwpcHitFinder.h:124
BmnMwpcHitFinder::hNp_best_Ch
vector< TH1D * > hNp_best_Ch
Definition BmnMwpcHitFinder.h:123
BmnMwpcHitFinder::Finish
virtual void Finish()
Virtual method Finish.
Definition BmnMwpcHitFinder.cxx:3481
BmnMwpcHitFinder::SegmentParamAlignment
void SegmentParamAlignment(Int_t, Int_t *, Double_t ***, Float_t **)
Definition BmnMwpcHitFinder.cxx:3205
BmnMwpcHitFinder::FillFitMatrix
void FillFitMatrix(Int_t, Double_t **, Float_t **, Double_t *, Int_t *)
Definition BmnMwpcHitFinder.cxx:3310
BmnMwpcHitFinder::FillFreeCoefVector
void FillFreeCoefVector(Int_t, Double_t *, Double_t ***, Int_t, Float_t **, Double_t *, Int_t *)
Definition BmnMwpcHitFinder.cxx:3358
BmnMwpcHitFinder::hResiduals_pl0_Ch
vector< TH1D * > hResiduals_pl0_Ch
Definition BmnMwpcHitFinder.h:125
BmnMwpcHitFinder::vec_points
vector< MC_points > vec_points
Definition BmnMwpcHitFinder.h:141
BmnMwpcHitFinder::SegmentsStoring
void SegmentsStoring(Int_t *, Double_t ***, Double_t **, Int_t **, Double_t ***, Double_t ***, Double_t ****)
Definition BmnMwpcHitFinder.cxx:920
BmnMwpcHitFinder::hEvent_display_Ch
vector< TH2D * > hEvent_display_Ch
Definition BmnMwpcHitFinder.h:131
BmnMwpcHitFinder::fBmnEventHeaderBranchName
TString fBmnEventHeaderBranchName
Definition BmnMwpcHitFinder.h:90
BmnMwpcHitFinder::hoccupancyVp
vector< TH1D * > hoccupancyVp
Definition BmnMwpcHitFinder.h:123
BmnMwpcHitFinder::hoccupancyUp
vector< TH1D * > hoccupancyUp
Definition BmnMwpcHitFinder.h:123
BmnMwpcHitFinder::Chi2_ndf_best_Ch
Double_t ** Chi2_ndf_best_Ch
Definition BmnMwpcHitFinder.h:165
BmnMwpcHitFinder::InverseMatrix
void InverseMatrix(Double_t **, Double_t **)
Definition BmnMwpcHitFinder.cxx:3424
BmnMwpcHitFinder::hNbest_Ch
vector< TH1D * > hNbest_Ch
Definition BmnMwpcHitFinder.h:123
BmnMwpcHitFinder::hClusterSize
vector< TH1D * > hClusterSize
Definition BmnMwpcHitFinder.h:124
BmnMwpcHitFinder::hdAy_mc_SegCh
vector< TH1D * > hdAy_mc_SegCh
Definition BmnMwpcHitFinder.h:128
BmnMwpcHitFinder::hOccupancy
vector< TH1D * > hOccupancy
Definition BmnMwpcHitFinder.h:123
BmnMwpcHitFinder::hfired_wire_Ch
vector< TH1D * > hfired_wire_Ch
Definition BmnMwpcHitFinder.h:126
BmnMwpcHitFinder::fBmnMwpcDigitArray
TClonesArray * fBmnMwpcDigitArray
Input array of MWPC digits.
Definition BmnMwpcHitFinder.h:93
BmnMwpcHitFinder::hpar_By_Ch
vector< TH1D * > hpar_By_Ch
Definition BmnMwpcHitFinder.h:127
BmnMwpcHitFinder::sigma2_seg
Double_t **** sigma2_seg
Definition BmnMwpcHitFinder.h:185
BmnMwpcHitFinder::hpar_Ax_Ch
vector< TH1D * > hpar_Ax_Ch
Definition BmnMwpcHitFinder.h:126
BmnMwpcHitFinder::hWiresUm
vector< TH1D * > hWiresUm
Definition BmnMwpcHitFinder.h:125
BmnMwpcHitFinder::hpar_Ay_Ch
vector< TH1D * > hpar_Ay_Ch
Definition BmnMwpcHitFinder.h:127
BmnMwpcHitFinder::hWiresXm
vector< TH1D * > hWiresXm
Definition BmnMwpcHitFinder.h:125
BmnMwpcHitFinder::MCefficiencyCalculation
void MCefficiencyCalculation(Int_t, vector< MC_points > &, Double_t ***, Int_t *)
Definition BmnMwpcHitFinder.cxx:979
BmnMwpcHitFinder::hoccupancyXm
vector< TH1D * > hoccupancyXm
Definition BmnMwpcHitFinder.h:123
BmnMwpcHitFinder::SegmentFinder
void SegmentFinder(Int_t, Int_t **, Double_t ***, Int_t ***, Int_t *, Double_t ***, Double_t ***, Int_t **, Int_t, Short_t, Int_t)
Definition BmnMwpcHitFinder.cxx:1923
BmnMwpcHitFinder::Cluster_seg
Double_t *** Cluster_seg
Definition BmnMwpcHitFinder.h:183
BmnMwpcHitFinder::hNFired_layers_Ch
vector< TH1D * > hNFired_layers_Ch
Definition BmnMwpcHitFinder.h:124
BmnMwpcHitFinder::hoccupancyXp
vector< TH1D * > hoccupancyXp
Definition BmnMwpcHitFinder.h:123
BmnMwpcHitFinder::hTime
vector< TH1D * > hTime
Definition BmnMwpcHitFinder.h:123
BmnMwpcHitFinder::hdX_mc_Seg_deltaCh
vector< TH1D * > hdX_mc_Seg_deltaCh
Definition BmnMwpcHitFinder.h:129
BmnMwpcHitFinder::FillFreeCoefVectorXUV
void FillFreeCoefVectorXUV(Int_t, Double_t *, Float_t **, Float_t **, Float_t *, Int_t *)
Definition BmnMwpcHitFinder.cxx:3394
BmnMwpcHitFinder::Chi2_ndf_seg
Double_t ** Chi2_ndf_seg
Definition BmnMwpcHitFinder.h:181
BmnMwpcHit::GetPdgId
Int_t GetPdgId()
Definition BmnMwpcHit.h:54
BmnMwpcHit::GetMwpcId
Short_t GetMwpcId() const
Definition BmnMwpcHit.h:32
BmnMwpcHit::GetWireNumber
Int_t GetWireNumber()
Definition BmnMwpcHit.h:50
BmnMwpcHit::GetPlaneId
Int_t GetPlaneId()
Definition BmnMwpcHit.h:53
BmnMwpcHit::GetHitId
Int_t GetHitId() const
Definition BmnMwpcHit.h:34
BmnMwpcHit::GetWireTime
Int_t GetWireTime()
Definition BmnMwpcHit.h:52
BmnMwpcTrack::SetClust
void SetClust(vector< Double_t > Clust)
Definition BmnMwpcTrack.h:23
BmnMwpcTrack::SetCoord
void SetCoord(vector< Double_t > Coord)
Definition BmnMwpcTrack.h:21
BmnTrack::SetChi2
void SetChi2(Double_t chi2)
Definition BmnTrack.h:76
BmnTrack::SetFlag
void SetFlag(Int_t flag)
Definition BmnTrack.h:74
BmnTrack::SetNHits
void SetNHits(Int_t n)
Definition BmnTrack.h:80
BmnTrack::SetParamFirst
void SetParamFirst(FairTrackParam &par)
Modifiers.
Definition BmnTrack.h:70
TFile::Close
void Close(Option_t *option="") override
abs
RVec< PromoteType< T > > abs(const RVec< T > &v)
sqrt
RVec< PromoteType< T > > sqrt(const RVec< T > &v)
std
Definition json.hpp:4446
BmnMwpcHitFinder::MC_points
Definition BmnMwpcHitFinder.h:57
BmnMwpcHitFinder::MC_points::z2
Double_t z2[6]
Definition BmnMwpcHitFinder.h:74
BmnMwpcHitFinder::MC_points::Pdg
Int_t Pdg
Definition BmnMwpcHitFinder.h:59
BmnMwpcHitFinder::MC_points::y3
Double_t y3[6]
Definition BmnMwpcHitFinder.h:76
BmnMwpcHitFinder::MC_points::x2
Double_t x2[6]
Definition BmnMwpcHitFinder.h:72
BmnMwpcHitFinder::MC_points::Np2
Int_t Np2
Definition BmnMwpcHitFinder.h:60
BmnMwpcHitFinder::MC_points::Id
Int_t Id
Definition BmnMwpcHitFinder.h:58
BmnMwpcHitFinder::MC_points::x3
Double_t x3[6]
Definition BmnMwpcHitFinder.h:75
BmnMwpcHitFinder::MC_points::Np3
Int_t Np3
Definition BmnMwpcHitFinder.h:61
BmnMwpcHitFinder::MC_points::y2
Double_t y2[6]
Definition BmnMwpcHitFinder.h:73
BmnMwpcHitFinder::MC_points::z3
Double_t z3[6]
Definition BmnMwpcHitFinder.h:77
segments::param
Double_t param[4]
Definition BmnMwpcHitFinder.cxx:40
segments::clust
Double_t clust[6]
Definition BmnMwpcHitFinder.cxx:39
segments::coord
Double_t coord[6]
Definition BmnMwpcHitFinder.cxx:38
segments::sigma2
Double_t sigma2[4][4]
Definition BmnMwpcHitFinder.cxx:41