BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
read_mwpc1_688.C
Go to the documentation of this file.
1#include "TBranch.h"
2#include "TCanvas.h"
3#include "TChain.h"
4#include "TClonesArray.h"
5#include "TFile.h"
6#include "TGeoManager.h"
7#include "TH1F.h"
8#include "TH2F.h"
9#include "TMath.h"
10#include "TTree.h"
11
12#include <fstream>
13#include <iostream>
14#include <stdio.h>
15using namespace std;
16using namespace TMath;
17
18const Int_t nWires = 96;
19const Int_t nPlanes = 6;
20const Int_t big = 100;
21const Float_t dw = 0.25; // cm
22const Double_t sq3 = Sqrt(3.);
23const Double_t Chi2_Max = 20;
24// Geometry for RUN3
25const Float_t Zmid = 16.5; // average plane of chamber 2
26const Float_t Zmid3 = -16.5; // average plane of chamber 3
27const Float_t Z_Ch1 = -150; // average plane of chamber 1
28const Float_t Z_Ch2 = 453; // global
29const Float_t Z_Ch3 = 486; // global
30
31void SegmentFinder(Int_t chNum,
32 Int_t wires_Ch[nPlanes][big],
33 Float_t XVU_Ch[nPlanes][big],
34 Int_t Nhits_Ch[big],
35 Int_t iw_Ch[nPlanes],
36 Int_t& Nseg,
37 Int_t wires_glob[nWires][nPlanes],
38 Float_t xuv_glob[nWires][nPlanes],
39 Int_t minHits,
40 Short_t code)
41{
42 // xuv_glob - coordinates of all hits
43 // wires_glob - wires of all hits
44 Float_t delta = (chNum == 1) ? 0.5 : 1.15; // condition for the nearest "lighted" wires
45 Int_t cntr_sumWW_x = (chNum == 1) ? 95 : 96; // center wires Wx1+Wx2 due to x-slope
46 Int_t min_sumWW_x = (chNum == 1) ? 2 : 3; // min sum of wires to 95 ^ 96
47 Int_t cntr_sumWW = 95; // center wires W1+W2 due to u,v-slope
48 Int_t min_sumWW = 2; // min sum of wires to 95
49 Int_t minHits4_5 = minHits;
50
51 // code :
52 // 1 {X-, V-, U+}
53 // 2 {X-, V+, U+}
54 // 3 {X-, V-, U-}
55 // 4 {X+, V-, U+}
56 // 5 {X+, V+, U+}
57 // 6 {X+, V-, U-}
58 // 7 {X+, V-, U+}
59 // 8 {X-, V-, U-}
60 Int_t x, v, u, x1, v1, u1;
61
62 switch (code) {
63 case 1:
64 x = 0;
65 v = 1;
66 u = 2;
67 x1 = 3;
68 v1 = 4;
69 u1 = 5;
70 break;
71 case 2:
72 x = 0;
73 v = 4;
74 u = 2;
75 x1 = 3;
76 v1 = 1;
77 u1 = 5;
78 break;
79 case 3:
80 x = 0;
81 v = 1;
82 u = 5;
83 x1 = 3;
84 v1 = 4;
85 u1 = 2;
86 break;
87 case 4:
88 x = 3;
89 v = 1;
90 u = 2;
91 x1 = 0;
92 v1 = 4;
93 u1 = 5;
94 break;
95 case 5:
96 x = 3;
97 v = 4;
98 u = 2;
99 x1 = 0;
100 v1 = 1;
101 u1 = 5;
102 break;
103 case 6:
104 x = 3;
105 v = 1;
106 u = 5;
107 x1 = 0;
108 v1 = 4;
109 u1 = 2;
110 break;
111 case 7:
112 x = 0;
113 v = 4;
114 u = 5;
115 x1 = 3;
116 v1 = 1;
117 u1 = 2;
118 break;
119 case 8:
120 x = 3;
121 v = 4;
122 u = 5;
123 x1 = 0;
124 v1 = 1;
125 u1 = 2;
126 break;
127 }
128
129 if (iw_Ch[x] > 0) {
130 for (Int_t ix = 0; ix < iw_Ch[x]; ix++) {
131 if (iw_Ch[v] > 0) {
132 for (Int_t iv = 0; iv < iw_Ch[v]; iv++) {
133 if (iw_Ch[u] > 0) {
134 for (Int_t iu = 0; iu < iw_Ch[u]; iu++) {
135 if (Nseg > big - 1)
136 break;
137 Bool_t it_was = 0;
138 if (Nseg > 0) {
139 for (Int_t iseg = 0; iseg < Nseg; iseg++) {
140 if (wires_Ch[x][iseg] == wires_glob[ix][x] && wires_Ch[v][iseg] == wires_glob[iv][v]
141 && wires_Ch[u][iseg] == wires_glob[iu][u])
142 {
143 it_was = 1;
144 break;
145 }
146 }
147 }
148 if (it_was)
149 continue;
150 if (fabs(xuv_glob[iu][u] + xuv_glob[iv][v] - xuv_glob[ix][x]) < delta) {
151 XVU_Ch[x][Nseg] = xuv_glob[ix][x];
152 XVU_Ch[v][Nseg] = xuv_glob[iv][v];
153 XVU_Ch[u][Nseg] = xuv_glob[iu][u];
154 wires_Ch[x1][Nseg] = -1;
155 wires_Ch[v1][Nseg] = -1;
156 wires_Ch[u1][Nseg] = -1;
157 wires_Ch[x][Nseg] = wires_glob[ix][x];
158 wires_Ch[v][Nseg] = wires_glob[iv][v];
159 wires_Ch[u][Nseg] = wires_glob[iu][u];
160 Nhits_Ch[Nseg] = 3;
161 if (iw_Ch[u1] > 0) {
162 for (Int_t iu2 = 0; iu2 < iw_Ch[u1]; iu2++) {
163
164 if (abs(wires_glob[iu][u] + wires_glob[iu2][u1] - cntr_sumWW) < min_sumWW) {
165 XVU_Ch[u1][Nseg] = xuv_glob[iu2][u1];
166 wires_Ch[u1][Nseg] = wires_glob[iu2][u1];
167 Nhits_Ch[Nseg] = 4;
168 }
169 }
170 }
171 if (iw_Ch[v1] > 0) {
172 for (Int_t iv2 = 0; iv2 < iw_Ch[v1]; iv2++) {
173
174 if (abs(wires_glob[iv][v] + wires_glob[iv2][v1] - cntr_sumWW) < min_sumWW) {
175 XVU_Ch[v1][Nseg] = xuv_glob[iv2][v1];
176 wires_Ch[v1][Nseg] = wires_glob[iv2][v1];
177 }
178 }
179 if (wires_Ch[v1][Nseg] > -1)
180 Nhits_Ch[Nseg] = Nhits_Ch[Nseg] + 1; // 5 points
181 }
182 if (iw_Ch[x1] > 0) {
183 for (Int_t ix2 = 0; ix2 < iw_Ch[x1]; ix2++) {
184
185 if (abs(wires_glob[ix][x] + wires_glob[ix2][x1] - cntr_sumWW_x) < min_sumWW_x) {
186 XVU_Ch[x1][Nseg] = xuv_glob[ix2][x1];
187 wires_Ch[x1][Nseg] = wires_glob[ix2][x1];
188 }
189 }
190 if (wires_Ch[x1][Nseg] > -1)
191 Nhits_Ch[Nseg] = Nhits_Ch[Nseg] + 1; // 6 points
192 }
193 if (Nhits_Ch[Nseg] >= minHits4_5) {
194 Nseg++;
195 if (Nseg > 25)
196 minHits4_5 = 5;
197 } else {
198 Nhits_Ch[Nseg] = 0;
199 wires_Ch[x1][Nseg] = -1;
200 wires_Ch[v1][Nseg] = -1;
201 wires_Ch[u1][Nseg] = -1;
202 wires_Ch[x][Nseg] = -1;
203 wires_Ch[v][Nseg] = -1;
204 wires_Ch[u][Nseg] = -1;
205 }
206 }
207 } // iu
208 } // u
209 } // iv
210 } // v
211 } // ix
212 } // x
213}
214
215void FillFitMatrix(Double_t A[4][4], Float_t z[nPlanes], Float_t sigm2[nPlanes], Int_t h[nPlanes])
216{
217
218 // A - matrix to be filled
219 // z - local z-positions of planes(layers)
220 // sigm2 - square of sigma
221 // h - array to include/exclude planes (h[i] = 0 or 1)
222
223 Float_t z2[nPlanes] = {z[0] * z[0], z[1] * z[1], z[2] * z[2], z[3] * z[3], z[4] * z[4], z[5] * z[5]}; // cm
224
225 A[0][0] = 2 * z2[0] * h[0] / sigm2[0] + z2[2] * h[2] / (2 * sigm2[2]) + z2[1] * h[1] / (2 * sigm2[1])
226 + 2 * z2[3] * h[3] / sigm2[3] + z2[5] * h[5] / (2 * sigm2[5]) + z2[4] * h[4] / (2 * sigm2[4]); // Ax
227 A[0][1] = 2 * z[0] * h[0] / sigm2[0] + z[2] * h[2] / (2 * sigm2[2]) + z[1] * h[1] / (2 * sigm2[1])
228 + 2 * z[3] * h[3] / (sigm2[3]) + z[5] * h[5] / (2 * sigm2[5]) + z[4] * h[4] / (2 * sigm2[4]); // Bx
229 A[0][2] = sq3
230 * (z2[2] * h[2] / (2 * sigm2[2]) - z2[1] * h[1] / (2 * sigm2[1]) + z2[5] * h[5] / (2 * sigm2[5])
231 - z2[4] * h[4] / (2 * sigm2[4])); // Ay
232 A[0][3] = sq3
233 * (z[2] * h[2] / (2 * sigm2[2]) - z[1] * h[1] / (2 * sigm2[1]) + z[5] * h[5] / (2 * sigm2[5])
234 - z[4] * h[4] / (2 * sigm2[4])); // By
235 A[1][0] = A[0][1];
236 A[1][1] = 2 * h[0] / sigm2[0] + 0.5 * h[2] / sigm2[2] + 0.5 * h[1] / sigm2[1] + 2 * h[3] / sigm2[3]
237 + 0.5 * h[5] / sigm2[5] + 0.5 * h[4] / sigm2[4];
238 A[1][2] =
239 sq3 * (z[2] * h[2] / sigm2[2] - z[1] * h[1] / sigm2[1] + z[5] * h[5] / sigm2[5] - z[4] * h[4] / sigm2[4]) * 0.5;
240 A[1][3] = sq3 * (h[2] / sigm2[2] - h[1] / sigm2[1] + h[5] / sigm2[5] - h[4] / sigm2[4]) * 0.5;
241 A[2][0] = A[0][2];
242 A[2][1] = A[1][2];
243 A[2][2] = 3.0
244 * (z2[2] * h[2] / sigm2[2] + z2[1] * h[1] / sigm2[1] + z2[5] * h[5] / sigm2[5] + z2[4] * h[4] / sigm2[4])
245 * 0.5;
246 A[2][3] =
247 3.0 * (z[2] * h[2] / sigm2[2] + z[1] * h[1] / sigm2[1] + z[5] * h[5] / sigm2[5] + z[4] * h[4] / sigm2[4]) * 0.5;
248 A[3][0] = A[0][3];
249 A[3][1] = A[1][3];
250 A[3][2] = A[2][3];
251 A[3][3] = 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]);
252}
253
254void FillFreeCoefVector(Double_t F[4],
255 Float_t XVU_Ch[nPlanes][big],
256 Int_t segIdx,
257 Float_t z[nPlanes],
258 Float_t sigm2[nPlanes],
259 Int_t h[nPlanes])
260{
261 // F - vector to be filled
262 // XVU_Ch - coordinates of segment in chamber (Is it correct definition?)
263 // segIdx - index of current segment
264 // z - local z-positions of planes(layers)
265 // sigm2 - square of sigma
266 // h - array to include/exclude planes (h[i] = 0 or 1)
267
268 F[0] = 2 * XVU_Ch[0][segIdx] * z[0] * h[0] / sigm2[0] + XVU_Ch[1][segIdx] * z[1] * h[1] / sigm2[1]
269 + XVU_Ch[2][segIdx] * z[2] * h[2] / sigm2[2] + 2 * XVU_Ch[3][segIdx] * z[3] * h[3] / sigm2[3]
270 + XVU_Ch[4][segIdx] * z[4] * h[4] / sigm2[4] + XVU_Ch[5][segIdx] * z[5] * h[5] / sigm2[5];
271 F[1] = 2 * XVU_Ch[0][segIdx] * h[0] / sigm2[0] + XVU_Ch[1][segIdx] * h[1] / sigm2[1]
272 + XVU_Ch[2][segIdx] * h[2] / sigm2[2] + 2 * XVU_Ch[3][segIdx] * h[3] / sigm2[3]
273 + XVU_Ch[4][segIdx] * h[4] / sigm2[4] + XVU_Ch[5][segIdx] * h[5] / sigm2[5];
274 F[2] = sq3
275 * (-XVU_Ch[1][segIdx] * z[1] * h[1] / sigm2[1] + XVU_Ch[2][segIdx] * z[2] * h[2] / sigm2[2]
276 - XVU_Ch[4][segIdx] * z[4] * h[4] / sigm2[4] + XVU_Ch[5][segIdx] * z[5] * h[5] / sigm2[5]);
277 F[3] = sq3
278 * (-XVU_Ch[1][segIdx] * h[1] / sigm2[1] + XVU_Ch[2][segIdx] * h[2] / sigm2[2]
279 - XVU_Ch[4][segIdx] * h[4] / sigm2[4] + XVU_Ch[5][segIdx] * h[5] / sigm2[5]);
280}
281
282void InverseMatrix(Double_t A[4][4], Double_t b[4][4])
283{
284 /**** Gaussian algorithm for 4x4 matrix inversion ****/
285
286 Double_t factor;
287 Double_t temp[4];
288
289 // Set b to I
290 for (Int_t i1 = 0; i1 < 4; i1++)
291 for (Int_t j1 = 0; j1 < 4; j1++)
292 if (i1 == j1)
293 b[i1][j1] = 1.0;
294 else
295 b[i1][j1] = 0.0;
296
297 for (Int_t i1 = 0; i1 < 4; i1++) {
298 for (Int_t j1 = i1 + 1; j1 < 4; j1++) {
299 if (fabs(A[i1][i1]) < fabs(A[j1][i1])) {
300 for (Int_t l1 = 0; l1 < 4; l1++)
301 temp[l1] = A[i1][l1];
302 for (Int_t l1 = 0; l1 < 4; l1++)
303 A[i1][l1] = A[j1][l1];
304 for (Int_t l1 = 0; l1 < 4; l1++)
305 A[j1][l1] = temp[l1];
306 for (Int_t l1 = 0; l1 < 4; l1++)
307 temp[l1] = b[i1][l1];
308 for (Int_t l1 = 0; l1 < 4; l1++)
309 b[i1][l1] = b[j1][l1];
310 for (Int_t l1 = 0; l1 < 4; l1++)
311 b[j1][l1] = temp[l1];
312 }
313 }
314 factor = A[i1][i1];
315 for (Int_t j1 = 4 - 1; j1 > -1; j1--) {
316 b[i1][j1] /= factor;
317 A[i1][j1] /= factor;
318 }
319 for (Int_t j1 = i1 + 1; j1 < 4; j1++) {
320 factor = -A[j1][i1];
321 for (Int_t k1 = 0; k1 < 4; k1++) {
322 A[j1][k1] += A[i1][k1] * factor;
323 b[j1][k1] += b[i1][k1] * factor;
324 }
325 }
326 } // i1
327 for (Int_t i1 = 3; i1 > 0; i1--) {
328 for (Int_t j1 = i1 - 1; j1 > -1; j1--) {
329 factor = -A[j1][i1];
330 for (Int_t k1 = 0; k1 < 4; k1++) {
331 A[j1][k1] += A[i1][k1] * factor;
332 b[j1][k1] += b[i1][k1] * factor;
333 }
334 }
335 } // i1
336 // end inverse
337}
338
339void ProcessSegments(Int_t chNum,
340 Float_t sigma,
341 Float_t z_loc[nPlanes],
342 Int_t Min_hits,
343 Int_t& Nseg,
344 Int_t Nhits_Ch[big],
345 Int_t Wires_Ch[nPlanes][big],
346 Float_t XVU_Ch[nPlanes][big],
347 Int_t& Nbest_Ch,
348 Int_t ind_best_Ch[5],
349 Double_t Chi2_ndf_Ch[big],
350 Double_t Chi2_ndf_best_Ch[5],
351 Double_t par_ab_Ch[4][big])
352{
353 Float_t sigma2 = sigma * sigma;
354 Float_t sigm2[nPlanes] = {sigma2, sigma2, sigma2, sigma2, sigma2, sigma2};
355
356 Int_t h6[nPlanes] = {1, 1, 1, 1, 1, 1};
357
358 for (Int_t Nhitm = nPlanes; Nhitm > Min_hits - 1; Nhitm--) {
359 Int_t ifNhitm = 0;
360
361 for (Int_t iseg = 0; iseg < Nseg; iseg++) {
362
363 if (Nhits_Ch[iseg] != Nhitm)
364 continue;
365 ifNhitm = 1;
366
367 Int_t h[nPlanes] = {0, 0, 0, 0, 0, 0};
368 for (Int_t i = 0; i < nPlanes; i++)
369 if (Wires_Ch[i][iseg] > -1)
370 h[i] = 1;
371
372 // linear fit
373 Double_t A[4][4]; // coef matrix
374 Double_t F[4]; // free coef
375
376 if (Nhits_Ch[iseg] == nPlanes) // case 6-point segment
377 FillFitMatrix(A, z_loc, sigm2, h6);
378 else
379 FillFitMatrix(A, z_loc, sigm2, h);
380
381 FillFreeCoefVector(F, XVU_Ch, iseg, z_loc, sigm2, h);
382
383 Double_t A0[4][4];
384 for (Int_t i1 = 0; i1 < 4; i1++)
385 for (Int_t j1 = 0; j1 < 4; j1++)
386 A0[i1][j1] = A[i1][j1];
387
388 Double_t b[4][4];
389 InverseMatrix(A, b);
390
391 Double_t sum;
392 Double_t A1[4][4];
393
394 for (Int_t i1 = 0; i1 < 4; ++i1)
395 for (j1 = 0; j1 < 4; ++j1) {
396 sum = 0;
397 for (Int_t k1 = 0; k1 < 4; ++k1)
398 sum += A0[i1][k1] * b[k1][j1];
399 A1[i1][j1] = sum;
400 }
401
402 for (Int_t i1 = 0; i1 < 4; i1++) {
403 par_ab_Ch[i1][iseg] = 0;
404 for (j1 = 0; j1 < 4; j1++) {
405 par_ab_Ch[i1][iseg] += b[i1][j1] * F[j1];
406 }
407 } // i1
408
409 Float_t dX[nPlanes];
410 Chi2_ndf_Ch[iseg] = 0;
411 for (Int_t i1 = 0; i1 < nPlanes; i1++) {
412 if (Wires_Ch[i1][iseg] > -1) {
413 if (i1 == 0 || i1 == 3)
414 dX[i1] = XVU_Ch[i1][iseg] - par_ab_Ch[0][iseg] * z_loc[i1] - par_ab_Ch[1][iseg];
415 if (i1 == 2 || i1 == 5)
416 dX[i1] = XVU_Ch[i1][iseg] - 0.5 * (par_ab_Ch[0][iseg] + sq3 * par_ab_Ch[2][iseg]) * z_loc[i1]
417 - 0.5 * (par_ab_Ch[1][iseg] + sq3 * par_ab_Ch[3][iseg]);
418 if (i1 == 1 || i1 == 4)
419 dX[i1] = XVU_Ch[i1][iseg] - 0.5 * (par_ab_Ch[0][iseg] - sq3 * par_ab_Ch[2][iseg]) * z_loc[i1]
420 - 0.5 * (par_ab_Ch[1][iseg] - sq3 * par_ab_Ch[3][iseg]);
421 Chi2_ndf_Ch[iseg] = Chi2_ndf_Ch[iseg] + dX[i1] * dX[i1] / sigm2[i1];
422 }
423 }
424
425 if (Nhits_Ch[iseg] > 4)
426 Chi2_ndf_Ch[iseg] = Chi2_ndf_Ch[iseg] / (Nhits_Ch[iseg] - 4);
427
428 if (Chi2_ndf_Ch[iseg] > Chi2_Max) {
429 if (Nhits_Ch[iseg] <= Min_hits)
430 Nhits_Ch[iseg] = 0;
431 else { // reject most distant point
432 Float_t Max_dev = 0;
433 Int_t irej = -1;
434 for (Int_t i1 = 0; i1 < nPlanes; i1++)
435 if (Wires_Ch[i1][iseg] > -1)
436 if (fabs(dX[i1]) > Max_dev) {
437 irej = i1;
438 Max_dev = fabs(dX[i1]);
439 }
440 Wires_Ch[irej][iseg] = -1;
441 Nhits_Ch[iseg]--; // reject most distant point
442 Chi2_ndf_Ch[iseg] = 0; // recalculate Chi2 with old(?) parameters
443 for (Int_t ii = 0; ii < nPlanes; ii++) {
444 if (Wires_Ch[ii][iseg] > -1) {
445 if (ii == 0 || ii == 3)
446 dX[ii] = XVU_Ch[ii][iseg] - par_ab_Ch[0][iseg] * z_loc[ii] - par_ab_Ch[1][iseg];
447 if (ii == 2 || ii == 5)
448 dX[ii] = XVU_Ch[ii][iseg]
449 - 0.5 * (par_ab_Ch[0][iseg] + sq3 * par_ab_Ch[2][iseg]) * z_loc[ii]
450 - 0.5 * (par_ab_Ch[1][iseg] + sq3 * par_ab_Ch[3][iseg]);
451 if (ii == 1 || ii == 4)
452 dX[ii] = XVU_Ch[ii][iseg]
453 - 0.5 * (par_ab_Ch[0][iseg] - sq3 * par_ab_Ch[2][iseg]) * z_loc[ii]
454 - 0.5 * (par_ab_Ch[1][iseg] - sq3 * par_ab_Ch[3][iseg]);
455 Chi2_ndf_Ch[iseg] = Chi2_ndf_Ch[iseg] + dX[ii] * dX[ii] / sigm2[ii];
456 }
457 }
458 } // else
459 }
460 } // iseg
461
462 if (!ifNhitm)
463 continue;
464
465 if (Nbest_Ch == 0) {
466
467 Double_t Chi2_best = 9999;
468 Int_t iseg_best = -1;
469 for (Int_t iseg = 0; iseg < Nseg; iseg++) {
470 if (Nhits_Ch[iseg] != Nhitm)
471 continue;
472 if (Chi2_ndf_Ch[iseg] >= Chi2_best)
473 continue;
474 Chi2_best = Chi2_ndf_Ch[iseg];
475 iseg_best = iseg;
476 } // iseg
477
478 if (iseg_best == -1)
479 continue;
480
481 ind_best_Ch[0] = iseg_best;
482 Chi2_ndf_best_Ch[0] = Chi2_best;
483 Nbest_Ch = 1;
484
485 // reject(common points)
486 for (Int_t iseg = 0; iseg < Nseg; iseg++) {
487 if (iseg == iseg_best)
488 continue;
489 Int_t u_cross = 0;
490 Int_t v_cross = 0;
491 Int_t x_cross = 0;
492 for (Int_t i1 = 0; i1 < nPlanes; i1++) {
493 if (Wires_Ch[i1][iseg] > -1) {
494 if (Wires_Ch[i1][iseg] == Wires_Ch[i1][iseg_best])
495 Nhits_Ch[iseg] = 0;
496 }
497 }
498 }
499 }
500
501 if (Nbest_Ch == 1) {
502 Double_t Chi2_best = 9999;
503 Int_t iseg_best = -1;
504 for (Int_t iseg = 0; iseg < Nseg; iseg++) {
505 if (iseg == ind_best_Ch[0])
506 continue;
507 if (Nhits_Ch[iseg] != Nhitm)
508 continue;
509 if (Chi2_ndf_Ch[iseg] > Chi2_best)
510 continue;
511 Chi2_best = Chi2_ndf_Ch[iseg];
512 iseg_best = iseg;
513 } // iseg
514
515 if (iseg_best > -1) {
516 ind_best_Ch[Nbest_Ch] = iseg_best;
517 Chi2_ndf_best_Ch[Nbest_Ch] = Chi2_best;
518 Nbest_Ch++;
519
520 if (chNum == 2 || chNum == 3) {
521 // reject(common points)
522 for (Int_t iseg = 0; iseg < Nseg; iseg++) {
523 if (iseg == iseg_best)
524 continue;
525 Int_t u_cross = 0;
526 Int_t v_cross = 0;
527 Int_t x_cross = 0;
528 for (Int_t i1 = 0; i1 < nPlanes; i1++) {
529 if (Wires_Ch[i1][iseg] > -1) {
530 if (Wires_Ch[i1][iseg] == Wires_Ch[i1][iseg_best])
531 Nhits_Ch[iseg] = 0;
532 }
533 }
534 } // iseg
535 }
536 }
537 }
538
539 if (chNum == 2 || chNum == 3) {
540 if (Nbest_Ch == 2) {
541
542 Double_t Chi2_best = 20;
543 Int_t iseg_best = -1;
544 for (Int_t iseg = 0; iseg < Nseg; iseg++) {
545 if (iseg == ind_best_Ch[0] || iseg == ind_best_Ch[1])
546 continue;
547 if (Nhits_Ch[iseg] != Nhitm)
548 continue;
549 if (Chi2_ndf_Ch[iseg] > Chi2_best)
550 continue;
551 Chi2_best = Chi2_ndf_Ch[iseg];
552 iseg_best = iseg;
553 } // iseg
554
555 if (iseg_best > -1) {
556
557 ind_best_Ch[Nbest_Ch] = iseg_best;
558 Chi2_ndf_best_Ch[Nbest_Ch] = Chi2_best;
559 Nbest_Ch++;
560
561 // reject(common points)
562 for (Int_t iseg = 0; iseg < Nseg; iseg++) {
563 if (iseg == iseg_best)
564 continue;
565 Int_t u_cross = 0;
566 Int_t v_cross = 0;
567 Int_t x_cross = 0;
568 for (Int_t i1 = 0; i1 < nPlanes; i1++) {
569 if (Wires_Ch[i1][iseg] > -1) {
570 if (Wires_Ch[i1][iseg] == Wires_Ch[i1][iseg_best])
571 Nhits_Ch[iseg] = 0;
572 }
573 }
574 }
575 }
576 }
577 }
578 } // Nhitm
579}
580void FillEfficiency(Int_t Wires_Ch[nPlanes][big],
581 Int_t Nhits_Ch[big],
582 Int_t segIdx,
583 TH1F* Denom_Ch,
584 TH1F* Nomin_Ch,
585 Int_t min_distX,
586 Int_t min_distY)
587{
588 // segIdx - index of current segment ch2 or ch3 // Int_t ind_best_Ch[5]
589 // 4p&4p -> all matched / Efficiency per layer
590
591 if (fabs(min_distX) < 4 && fabs(min_distY) < 3.5) {
592 for (int i1 = 0; i1 < 6; i1++) {
593
594 if (Wires_Ch[i1][segIdx] > -1 && Nhits_Ch[segIdx] == 4)
595 continue;
596
597 Denom_Ch->Fill(i1);
598
599 if (Wires_Ch[i1][segIdx] > -1) {
600 Nomin_Ch->Fill(i1);
601 }
602 }
603 } // min_distX
604} // FillEfficiency
605
606void FillResid(Int_t Wires_Ch[nPlanes][big],
607 Int_t segIdx,
608 Float_t XVU_Ch[nPlanes][big],
609 Double_t par_ab_Ch[4][big],
610 Float_t z_loc[nPlanes],
611 TH1F* hResid_Ch[nPlanes])
612{
613 Float_t dX[nPlanes];
614 for (Int_t i1 = 0; i1 < nPlanes; i1++) {
615 if (Wires_Ch[i1][segIdx] > -1) {
616 Float_t xvuBest = XVU_Ch[i1][segIdx];
617 if (i1 == 0 || i1 == 3)
618 dX[i1] = xvuBest - par_ab_Ch[0][segIdx] * z_loc[i1] - par_ab_Ch[1][segIdx];
619 if (i1 == 2 || i1 == 5)
620 dX[i1] = xvuBest - 0.5 * (par_ab_Ch[0][segIdx] + sq3 * par_ab_Ch[2][segIdx]) * z_loc[i1]
621 - 0.5 * (par_ab_Ch[1][segIdx] + sq3 * par_ab_Ch[3][segIdx]);
622 if (i1 == 1 || i1 == 4)
623 dX[i1] = xvuBest - 0.5 * (par_ab_Ch[0][segIdx] - sq3 * par_ab_Ch[2][segIdx]) * z_loc[i1]
624 - 0.5 * (par_ab_Ch[1][segIdx] - sq3 * par_ab_Ch[3][segIdx]);
625 hResid_Ch[i1]->Fill(dX[i1]);
626 }
627 }
628} // FillResid
629
630Int_t read_mwpc1_688()
631{
632
633 gROOT->LoadMacro("$VMCWORKDIR/macro/run/bmnloadlibs.C");
634 bmnloadlibs(); // load libraries
635
636 // Use :: modificator before global name
637
638 TFile in("/dataBMN/bmndata1/run3/root/bmn_run0688_hrb.root", "READ");
639
640 TTree* T = (TTree*)in.Get("BMN_RAW");
641
642 BmnMwpcDigit* bmndigit = new BmnMwpcDigit();
643 TClonesArray* clar = new TClonesArray("BmnMwpcDigit", 100);
644 TBranch* BMN = T->GetBranch("bmn_mwpc");
645 BMN->SetAddress(&clar);
646
647 TFile* fReco = new TFile("bmndst_run688.root", "RECREATE");
648 TTree* tReco = new TTree("cbmsim", "bmndst_run688");
649
650 // TClonesArray* dchHits = new TClonesArray("BmnDchHit");
651 TClonesArray* mwpcHits = new TClonesArray("BmnMwpcHit");
652 // TClonesArray* dchTracks = new TClonesArray("CbmTrack");
653 TClonesArray* mwpcTracks1 = new TClonesArray("CbmTrack");
654 TClonesArray* mwpcTracks2 = new TClonesArray("CbmTrack");
655 TClonesArray* mwpcTracks3 = new TClonesArray("CbmTrack");
656 TClonesArray* mwpcMatchedTracks = new TClonesArray("CbmTrack");
657
658 // tReco->Branch("BmnDchHit", &dchHits);
659 tReco->Branch("BmnMwpcHit", &mwpcHits);
660 // tReco->Branch("DchTracks", &dchTracks);
661 tReco->Branch("Mwpc1Segmants", &mwpcTracks1);
662 tReco->Branch("Mwpc2Segmants", &mwpcTracks2);
663 tReco->Branch("Mwpc3Segmants", &mwpcTracks3);
664 tReco->Branch("MwpcMatchedTracks", &mwpcMatchedTracks);
665
666 TH1F* Nomin_Ch2 = new TH1F("Nomin23_Ch2", "Nomin23_Ch2", 6, 0, 6);
667 TH1F* Nomin_Ch3 = new TH1F("Nomin23_Ch3", "Nomin23_Ch3", 6, 0, 6);
668
669 TH1F* Denom_Ch2 = new TH1F("Denom23_Ch2", "Denom23_Ch2", 6, 0, 6);
670 TH1F* Denom_Ch3 = new TH1F("Denom23_Ch3", "Denom23_Ch3", 6, 0, 6);
671
672 TH1F* Effic23_Ch2 = new TH1F("Effic23_Ch2", "Layer Efficiency23 (MWPC2)", 6, 0, 6);
673 TH1F* Effic23_Ch3 = new TH1F("Effic23_Ch3", "Layer Efficiency23 (MWPC3)", 6, 0, 6);
674
675 TH1F* Nomin_Ch2_9p = new TH1F("Nomin_Ch2_9p", "Nomin_Ch2_9p", 6, 0, 6);
676 TH1F* Nomin_Ch3_9p = new TH1F("Nomin_Ch3_9p", "Nomin_Ch3_9p", 6, 0, 6);
677
678 TH1F* Denom_Ch2_9p = new TH1F("Denom_Ch2_9p", "Denom_Ch2_9p", 6, 0, 6);
679 TH1F* Denom_Ch3_9p = new TH1F("Denom_Ch3_9p", "Denom_Ch3_9p", 6, 0, 6);
680
681 TH1F* Effic23_Ch2_9p = new TH1F("Effic23_Ch2_9p", "Layer Efficiency ch2 ch3 min9p (MWPC2)", 6, 0, 6);
682 TH1F* Effic23_Ch3_9p = new TH1F("Effic23_Ch3_9p", "Layer Efficiency ch2 ch3 min9p (MWPC3)", 6, 0, 6);
683
684 TH1F* Nomin23_overlap = new TH1F("Nomin23_overlap", "Nomin23_overlap", 2, 1, 3);
685 TH1F* Denom23_overlap = new TH1F("Denom23_overlap", "Denom23_overlap", 2, 1, 3);
686 TH1F* Effic23_overlap = new TH1F("Effic23_overlap", "Effic23_overlap", 2, 1, 3);
687
688 TH1F* hChi2_ndf_Ch23 = new TH1F("hChi2_ndf_Ch23", "Chi2_ndf_Ch23", 100, 0, 100);
689
690 TH1F* hResid_Ch1[6];
691 TH1F* hResid_Ch2[6];
692 TH1F* hResid_Ch3[6];
693 TH1F* hResid23_Ch2[6];
694 TH1F* hResid23_Ch3[6];
695
696 TString namHist;
697
698 for (int i1 = 0; i1 < 6; i1++) {
699 namHist = "hResid_Ch1_";
700 namHist += i1 + 1;
701 hResid_Ch1[i1] = new TH1F(namHist, namHist, 100, -0.50, 0.50);
702 namHist = "hResid_Ch2_";
703 namHist += i1 + 1;
704 hResid_Ch2[i1] = new TH1F(namHist, namHist, 100, -0.50, 0.50);
705 namHist = "hResid_Ch3_";
706 namHist += i1 + 1;
707 hResid_Ch3[i1] = new TH1F(namHist, namHist, 100, -0.50, 0.50);
708 namHist = "hResid23_Ch2_";
709 namHist += i1 + 1;
710 hResid23_Ch2[i1] = new TH1F(namHist, namHist, 100, -0.50, 0.50);
711 namHist = "hResid23_Ch3_";
712 namHist += i1 + 1;
713 hResid23_Ch3[i1] = new TH1F(namHist, namHist, 100, -0.50, 0.50);
714 }
715
716 Int_t Nevents = T->GetEntries();
717
718 //
719 // Nevents = 1000;
720
721 Int_t nentries;
722
723 Short_t wn, pn, ts, pl;
724 Int_t wire_Ch1[nWires][nPlanes] = {0};
725 Float_t xuv_Ch1[nWires][nPlanes] = {0};
726 Int_t wire_Ch2[nWires][nPlanes] = {0};
727 Float_t xuv_Ch2[nWires][nPlanes] = {0};
728 Int_t wire_Ch3[nWires][nPlanes] = {0};
729 Float_t xuv_Ch3[nWires][nPlanes] = {0};
730 Float_t dw = 0.25; // cm
731 Double_t sq3 = sqrt(3.);
732
733 for (Int_t iev = 1; iev < Nevents; ++iev) {
734
735 T->GetEntry(iev, 1);
736 if (iev % 100 == 0)
737 cout << " --- Ev " << iev << endl;
738
739 // dchHits->Clear();
740 mwpcHits->Clear();
741 // dchTracks->Clear();
742 mwpcTracks1->Clear();
743 mwpcTracks2->Clear();
744 mwpcTracks3->Clear();
745 mwpcMatchedTracks->Clear();
746
747 // Ch : 3 2 2 2 2 2 2 3 3 3 3 3 1 1 1 1 1 1
748 Int_t pln[18] = {15, 8, 7, 6, 11, 10, 9, 14, 13, 12, 17, 16, 0, 1, 2, 3, 4, 5}; // 9.IX.15
749
750 // x- v- u+ x+ v+ u-
751 Float_t z1_loc[nPlanes] = {-2.5, -1.5, -0.5, 0.5, 1.5, 2.5}; // cm
752 Float_t z2_loc[nPlanes] = {1.5, 0.5, -0.5, -1.5, -2.5, 2.5}; // cm
753 Float_t z3_loc[nPlanes] = {1.5, 2.5, -2.5, -1.5, -0.5, 0.5}; // cm
754
755 Int_t iw[18] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
756 Int_t iw_Ch1[nPlanes] = {0, 0, 0, 0, 0, 0};
757 Int_t iw_Ch2[nPlanes] = {0, 0, 0, 0, 0, 0};
758 Int_t iw_Ch3[nPlanes] = {0, 0, 0, 0, 0, 0};
759
760 nentries = clar->GetEntries();
761
762 for (Int_t j = 0; j < nentries; ++j) {
763 bmndigit = (BmnMwpcDigit*)clar->At(j);
764 wn = bmndigit->GetWireNumber();
765 pl = bmndigit->GetPlane();
766
767 pn = pln[pl];
768 Int_t pn0 = 6;
769 if (pn < 6)
770 pn0 = 0;
771 if (pn > 11)
772 pn0 = 12;
773 ts = bmndigit->GetTime();
774 bool repeat = 0;
775 if (iw[pn] > 0) {
776 for (Int_t ix = 0; ix < iw[pn]; ix++) {
777 if (wn == wire_Ch1[ix][pn - pn0] && pn0 == 0) {
778 repeat = 1;
779 break;
780 }
781 if (wn == wire_Ch2[ix][pn - pn0] && pn0 == 6) {
782 repeat = 1;
783 break;
784 }
785 if (wn == wire_Ch3[ix][pn - pn0] && pn0 == 12) {
786 repeat = 1;
787 break;
788 }
789 }
790 }
791 if (repeat)
792 continue;
793 if (pn < 6) {
794 // isCh1Work = kTRUE;
795 wire_Ch1[iw[pn]][pn] = wn;
796 xuv_Ch1[iw[pn]][pn] = (wn - 46.75) * dw;
797 } else if (pn > 11) {
798 // isCh2Work = kTRUE;
799 wire_Ch3[iw[pn]][pn - pn0] = wn;
800 xuv_Ch3[iw[pn]][pn - pn0] = (wn - 46.75) * dw;
801 } else {
802 // isCh3Work = kTRUE;
803 wire_Ch2[iw[pn]][pn - pn0] = wn;
804 xuv_Ch2[iw[pn]][pn - pn0] = (wn - 46.75) * dw;
805 }
806
807 if (pn == 0 || pn == 1 || pn == 5)
808 xuv_Ch1[iw[pn]][pn] = -xuv_Ch1[iw[pn]][pn];
809 if (pn == 6 || pn == 7 || pn == 11)
810 xuv_Ch2[iw[pn]][pn - pn0] = -xuv_Ch2[iw[pn]][pn - pn0];
811 if (pn == 12 || pn == 13 || pn == 17)
812 xuv_Ch3[iw[pn]][pn - pn0] = -xuv_Ch3[iw[pn]][pn - pn0];
813
814 iw[pn]++;
815 if (pn0 < 6)
816 iw_Ch1[pn - pn0] = iw[pn];
817 else if (pn0 > 11)
818 iw_Ch3[pn - pn0] = iw[pn];
819 else
820 iw_Ch2[pn - pn0] = iw[pn];
821 }
822
823 Int_t Min_hits = 5;
824 Double_t Chi2_Max = 20;
825
826 //--------- Ch1 ---------//
827 // if (isCh1Work) {
828 Float_t Z_to_target = 150; // cm
829
830 // slope and shift Ch1 for parametrs
831 Float_t x1_sh = 0;
832 Float_t y1_sh = 8.5; // cm
833 Float_t x1_slope_sh = -0.03; // + -> -
834 Float_t y1_slope_sh = 0.08;
835
836 Int_t Nseg_Ch1 = 0;
837 Int_t Wires_Ch1[nPlanes][big] = {-1};
838 Float_t XVU_Ch1[nPlanes][big];
839 Int_t Nhits_Ch1[big] = {0};
840
841 SegmentFinder(1, Wires_Ch1, XVU_Ch1, Nhits_Ch1, iw_Ch1, Nseg_Ch1, wire_Ch1, xuv_Ch1, Min_hits, 1);
842 SegmentFinder(1, Wires_Ch1, XVU_Ch1, Nhits_Ch1, iw_Ch1, Nseg_Ch1, wire_Ch1, xuv_Ch1, Min_hits, 2);
843 SegmentFinder(1, Wires_Ch1, XVU_Ch1, Nhits_Ch1, iw_Ch1, Nseg_Ch1, wire_Ch1, xuv_Ch1, Min_hits, 3);
844 SegmentFinder(1, Wires_Ch1, XVU_Ch1, Nhits_Ch1, iw_Ch1, Nseg_Ch1, wire_Ch1, xuv_Ch1, Min_hits, 4);
845 SegmentFinder(1, Wires_Ch1, XVU_Ch1, Nhits_Ch1, iw_Ch1, Nseg_Ch1, wire_Ch1, xuv_Ch1, Min_hits, 5);
846 SegmentFinder(1, Wires_Ch1, XVU_Ch1, Nhits_Ch1, iw_Ch1, Nseg_Ch1, wire_Ch1, xuv_Ch1, Min_hits, 6);
847 SegmentFinder(1, Wires_Ch1, XVU_Ch1, Nhits_Ch1, iw_Ch1, Nseg_Ch1, wire_Ch1, xuv_Ch1, Min_hits, 7);
848 SegmentFinder(1, Wires_Ch1, XVU_Ch1, Nhits_Ch1, iw_Ch1, Nseg_Ch1, wire_Ch1, xuv_Ch1, Min_hits, 8);
849
850 Int_t Nbest_Ch1 = 0;
851 Int_t ind_best_Ch1[5] = {0, 0, 0, 0, 0};
852 Double_t Chi2_ndf_best_Ch1[5] = {-999, -999, -999, -999, -999};
853 Double_t Chi2_ndf_Ch1[big];
854 Double_t par_ab_Ch1[4][big] = {999, 999, 999, 999};
855
856 if (Nseg_Ch1 > 0)
857 ProcessSegments(1, 0.0722, z1_loc, Min_hits, Nseg_Ch1, Nhits_Ch1, Wires_Ch1, XVU_Ch1, Nbest_Ch1,
858 ind_best_Ch1, Chi2_ndf_Ch1, Chi2_ndf_best_Ch1, par_ab_Ch1);
859
860 // create track and put it into TClonesArray
861 if (Nbest_Ch1 > 0) {
862 // for (Int_t iBest = 0; iBest < Nbest_Ch1; iBest++) {
863 Int_t iBest = 0;
864 if (Nhits_Ch1[ind_best_Ch1[iBest]] < 5)
865 continue;
866 Float_t dX[nPlanes];
867 // Int_t iBest = ind_best_Ch1[iBest];
868
869 FillResid(Wires_Ch1, ind_best_Ch1[iBest], XVU_Ch1, par_ab_Ch1, z1_loc, hResid_Ch1);
870
871 // local parameters to Global parameters
872 // ax alpha ax^2
873 par_ab_Ch1[0][ind_best_Ch1[iBest]] +=
874 x1_slope_sh + x1_slope_sh * par_ab_Ch1[0][ind_best_Ch1[iBest]] * par_ab_Ch1[0][ind_best_Ch1[iBest]];
875 par_ab_Ch1[2][ind_best_Ch1[iBest]] +=
876 y1_slope_sh + y1_slope_sh * par_ab_Ch1[2][ind_best_Ch1[iBest]] * par_ab_Ch1[2][ind_best_Ch1[iBest]];
877 par_ab_Ch1[1][ind_best_Ch1[iBest]] += x1_sh;
878 par_ab_Ch1[3][ind_best_Ch1[iBest]] += y1_sh;
879
880 Float_t X1_v = par_ab_Ch1[0][ind_best_Ch1[iBest]] * Z_to_target + par_ab_Ch1[1][ind_best_Ch1[iBest]];
881 Float_t Y1_v = par_ab_Ch1[2][ind_best_Ch1[iBest]] * Z_to_target + par_ab_Ch1[3][ind_best_Ch1[iBest]];
882
883 CbmTrack TracksCh1;
884 FairTrackParam ParamsTrackCh1;
885 Float_t z0 = Z_Ch1;
886 Float_t x0 = par_ab_Ch1[1][ind_best_Ch1[iBest]]; // FIXME!!! Check formula x0 - x coordinate of point in
887 // which line crosses mid-plane
888 Float_t y0 = par_ab_Ch1[3][ind_best_Ch1[iBest]]; // FIXME!!! Check formula y0 - y coordinate of point in
889 // which line crosses mid-plane
890
891 ParamsTrackCh1.SetPosition(TVector3(x0, y0, z0));
892 ParamsTrackCh1.SetTx(par_ab_Ch1[0][ind_best_Ch1[iBest]]);
893 ParamsTrackCh1.SetTy(par_ab_Ch1[2][ind_best_Ch1[iBest]]);
894 TracksCh1.SetParamFirst(&ParamsTrackCh1);
895 TracksCh1.SetChiSq(Chi2_ndf_Ch1[ind_best_Ch1[iBest]]);
896 TracksCh1.SetNDF(Nhits_Ch1[ind_best_Ch1[iBest]]); // WARNING!!! now it is number of hits in matched track!
897 new ((*mwpcTracks1)[mwpcTracks1->GetEntriesFast()]) CbmTrack(TracksCh1);
898
899 //}//iBest
900 } // Nbest_Ch1 > 0
901
902 // }
903 //--------- Ch2 ---------//
904 // if (isCh2Work) {
905
906 Float_t x2_sh = 1.0 - 73.0; //-0.6 //cm
907 Float_t y2_sh = 0;
908 Float_t x2_slope_sh = 0.030; // 0; + -> -//-0.0025//0.005; 0.015;//0.015;//0.02//.025// .032
909 Float_t y2_slope_sh = -0.01; // 0.006;//0.01// - 0.02;
910
911 Int_t Nbest_Ch2 = 0;
912 Int_t ind_best_Ch2[5] = {0, 0, 0, 0, 0};
913 Double_t par_ab_Ch2[4][100] = {999, 999, 999, 999};
914
915 Int_t Nseg_Ch2 = 0;
916 Int_t Wires_Ch2[nPlanes][big] = {-1};
917 Float_t XVU_Ch2[nPlanes][big];
918 Int_t Nhits_Ch2[big] = {0};
919 Min_hits = 4;
920
921 SegmentFinder(2, Wires_Ch2, XVU_Ch2, Nhits_Ch2, iw_Ch2, Nseg_Ch2, wire_Ch2, xuv_Ch2, Min_hits, 1);
922 SegmentFinder(2, Wires_Ch2, XVU_Ch2, Nhits_Ch2, iw_Ch2, Nseg_Ch2, wire_Ch2, xuv_Ch2, Min_hits, 2);
923 SegmentFinder(2, Wires_Ch2, XVU_Ch2, Nhits_Ch2, iw_Ch2, Nseg_Ch2, wire_Ch2, xuv_Ch2, Min_hits, 3);
924 SegmentFinder(2, Wires_Ch2, XVU_Ch2, Nhits_Ch2, iw_Ch2, Nseg_Ch2, wire_Ch2, xuv_Ch2, Min_hits, 4);
925 SegmentFinder(2, Wires_Ch2, XVU_Ch2, Nhits_Ch2, iw_Ch2, Nseg_Ch2, wire_Ch2, xuv_Ch2, Min_hits, 5);
926 SegmentFinder(2, Wires_Ch2, XVU_Ch2, Nhits_Ch2, iw_Ch2, Nseg_Ch2, wire_Ch2, xuv_Ch2, Min_hits, 6);
927 SegmentFinder(2, Wires_Ch2, XVU_Ch2, Nhits_Ch2, iw_Ch2, Nseg_Ch2, wire_Ch2, xuv_Ch2, Min_hits, 7);
928 SegmentFinder(2, Wires_Ch2, XVU_Ch2, Nhits_Ch2, iw_Ch2, Nseg_Ch2, wire_Ch2, xuv_Ch2, Min_hits, 8);
929
930 Double_t Chi2_ndf_best_Ch2[5] = {-999, -999, -999, -999, -999};
931 Double_t Chi2_ndf_Ch2[100];
932
933 if (Nseg_Ch2 > 0)
934 ProcessSegments(2, 0.0722, z2_loc, Min_hits, Nseg_Ch2, Nhits_Ch2, Wires_Ch2, XVU_Ch2, Nbest_Ch2,
935 ind_best_Ch2, Chi2_ndf_Ch2, Chi2_ndf_best_Ch2, par_ab_Ch2);
936
937 // create track and put it into TClonesArray
938 if (Nbest_Ch2 > 0) {
939 for (Int_t iBest = 0; iBest < Nbest_Ch2; iBest++) {
940
941 if (Nhits_Ch2[ind_best_Ch2[iBest]] > 4) {
942
943 FillResid(Wires_Ch2, ind_best_Ch2[iBest], XVU_Ch2, par_ab_Ch2, z2_loc, hResid_Ch2);
944 } // > 4
945 // local parameters to Global parameters
946 par_ab_Ch2[0][ind_best_Ch2[iBest]] +=
947 x2_slope_sh + x2_slope_sh * par_ab_Ch2[0][ind_best_Ch2[iBest]] * par_ab_Ch2[0][ind_best_Ch2[iBest]];
948 par_ab_Ch2[2][ind_best_Ch2[iBest]] +=
949 y2_slope_sh + y2_slope_sh * par_ab_Ch2[2][ind_best_Ch2[iBest]] * par_ab_Ch2[2][ind_best_Ch2[iBest]];
950 par_ab_Ch2[1][ind_best_Ch2[iBest]] += x2_sh;
951 par_ab_Ch2[3][ind_best_Ch2[iBest]] += y2_sh;
952
953 CbmTrack TracksCh2;
954 FairTrackParam ParamsTrackCh2;
955 Float_t z0 = Z_Ch2;
956 Float_t x0 = par_ab_Ch2[1][ind_best_Ch2[iBest]]; // Check formula x0 - x coordinate of point in which
957 // line crosses mid-plane
958 Float_t y0 = par_ab_Ch2[3][ind_best_Ch2[iBest]]; // Check formula y0 - y coordinate of point in which
959 // line crosses mid-plane
960 ParamsTrackCh2.SetPosition(TVector3(x0, y0, z0));
961 ParamsTrackCh2.SetTx(par_ab_Ch2[0][ind_best_Ch2[iBest]]);
962 ParamsTrackCh2.SetTy(par_ab_Ch2[2][ind_best_Ch2[iBest]]);
963 TracksCh2.SetParamFirst(&ParamsTrackCh2);
964 TracksCh2.SetChiSq(Chi2_ndf_Ch2[ind_best_Ch2[iBest]]);
965 TracksCh2.SetNDF(
966 Nhits_Ch2[ind_best_Ch2[iBest]]); // WARNING!!! now it is number of hits in matched track!
967 new ((*mwpcTracks2)[mwpcTracks2->GetEntriesFast()]) CbmTrack(TracksCh2);
968
969 } // iBest
970 } // Nbest_Ch2 > 0
971
972 //--------- Ch3 ---------//
973 // if (isCh3Work) {
974
975 Float_t x3_sh = -2.0 - 73.0; // 1.6 //cm
976 Float_t y3_sh = 0.1; // cm
977 Float_t x3_slope_sh = 0.038; // 0;// 0.005// 0.015// 0.028//.033 //.040
978 Float_t y3_slope_sh = -0.01; // 0.006// 0.01//-0.02
979
980 Int_t Nbest_Ch3 = 0;
981 Int_t ind_best_Ch3[5] = {0, 0, 0, 0, 0};
982 Double_t par_ab_Ch3[4][100] = {999, 999, 999, 999};
983
984 Int_t Nseg_Ch3 = 0;
985 Int_t Wires_Ch3[nPlanes][big] = {-1};
986 Float_t XVU_Ch3[nPlanes][big];
987 Int_t Nhits_Ch3[big] = {0};
988 Min_hits = 4;
989
990 SegmentFinder(3, Wires_Ch3, XVU_Ch3, Nhits_Ch3, iw_Ch3, Nseg_Ch3, wire_Ch3, xuv_Ch3, Min_hits, 1);
991 SegmentFinder(3, Wires_Ch3, XVU_Ch3, Nhits_Ch3, iw_Ch3, Nseg_Ch3, wire_Ch3, xuv_Ch3, Min_hits, 2);
992 SegmentFinder(3, Wires_Ch3, XVU_Ch3, Nhits_Ch3, iw_Ch3, Nseg_Ch3, wire_Ch3, xuv_Ch3, Min_hits, 3);
993 SegmentFinder(3, Wires_Ch3, XVU_Ch3, Nhits_Ch3, iw_Ch3, Nseg_Ch3, wire_Ch3, xuv_Ch3, Min_hits, 4);
994 SegmentFinder(3, Wires_Ch3, XVU_Ch3, Nhits_Ch3, iw_Ch3, Nseg_Ch3, wire_Ch3, xuv_Ch3, Min_hits, 5);
995 SegmentFinder(3, Wires_Ch3, XVU_Ch3, Nhits_Ch3, iw_Ch3, Nseg_Ch3, wire_Ch3, xuv_Ch3, Min_hits, 6);
996 SegmentFinder(3, Wires_Ch3, XVU_Ch3, Nhits_Ch3, iw_Ch3, Nseg_Ch3, wire_Ch3, xuv_Ch3, Min_hits, 7);
997 SegmentFinder(3, Wires_Ch3, XVU_Ch3, Nhits_Ch3, iw_Ch3, Nseg_Ch3, wire_Ch3, xuv_Ch3, Min_hits, 8);
998
999 Double_t Chi2_ndf_best_Ch3[5] = {-999, -999, -999, -999, -999};
1000 Double_t Chi2_ndf_Ch3[100];
1001
1002 if (Nseg_Ch3 > 0)
1003 ProcessSegments(3, 0.0722, z3_loc, Min_hits, Nseg_Ch3, Nhits_Ch3, Wires_Ch3, XVU_Ch3, Nbest_Ch3,
1004 ind_best_Ch3, Chi2_ndf_Ch3, Chi2_ndf_best_Ch3, par_ab_Ch3);
1005
1006 if (Nbest_Ch3 > 0) {
1007 for (Int_t iBest = 0; iBest < Nbest_Ch3; iBest++) {
1008 if (Nhits_Ch3[ind_best_Ch3[iBest]] > 4) {
1009
1010 FillResid(Wires_Ch3, ind_best_Ch3[iBest], XVU_Ch3, par_ab_Ch3, z3_loc, hResid_Ch3);
1011 } // > 4
1012
1013 // local parameters to Global parameters
1014 par_ab_Ch3[0][ind_best_Ch3[iBest]] +=
1015 x3_slope_sh + x3_slope_sh * par_ab_Ch3[0][ind_best_Ch3[iBest]] * par_ab_Ch3[0][ind_best_Ch3[iBest]];
1016 par_ab_Ch3[2][ind_best_Ch3[iBest]] +=
1017 y3_slope_sh + y3_slope_sh * par_ab_Ch3[2][ind_best_Ch3[iBest]] * par_ab_Ch3[2][ind_best_Ch3[iBest]];
1018 par_ab_Ch3[1][ind_best_Ch3[iBest]] += x3_sh;
1019 par_ab_Ch3[3][ind_best_Ch3[iBest]] += y3_sh;
1020
1021 CbmTrack TracksCh3;
1022 FairTrackParam ParamsTrackCh3;
1023 Float_t z0 = Z_Ch3;
1024 Float_t x0 = par_ab_Ch3[1][ind_best_Ch3[iBest]]; // Check formula x0 - x coordinate of point in which
1025 // line crosses mid-plane
1026 Float_t y0 = par_ab_Ch3[3][ind_best_Ch3[iBest]]; // Check formula y0 - y coordinate of point in which
1027 // line crosses mid-plane
1028 ParamsTrackCh3.SetPosition(TVector3(x0, y0, z0));
1029 ParamsTrackCh3.SetTx(par_ab_Ch3[0][ind_best_Ch3[iBest]]);
1030 ParamsTrackCh3.SetTy(par_ab_Ch3[2][ind_best_Ch3[iBest]]);
1031 TracksCh3.SetParamFirst(&ParamsTrackCh3);
1032 TracksCh3.SetChiSq(Chi2_ndf_Ch3[ind_best_Ch3[iBest]]);
1033 TracksCh3.SetNDF(
1034 Nhits_Ch3[ind_best_Ch3[iBest]]); // WARNING!!! now it is number of hits in matched track!
1035 new ((*mwpcTracks3)[mwpcTracks3->GetEntriesFast()]) CbmTrack(TracksCh3);
1036
1037 } // iBest++
1038 } // Nbest_Ch3 > 0
1039
1040 // === Ch2-Ch3 Matching
1041
1042 Int_t best = -1;
1043 Int_t best3 = -1;
1044
1045 Int_t best_Ch2_gl[2] = {-1, -1};
1046 Int_t best_Ch3_gl[2] = {-1, -1};
1047 Int_t Nbest_Ch23_gl = 0;
1048 Float_t Chi2_match[2] = {999, 999};
1049 Float_t sig_dx = 0.85;
1050 Float_t sig_dy = 0.76;
1051 Float_t sig_dax = 0.063;
1052 Float_t sig_day = 0.045;
1053
1054 if (Nbest_Ch2 > 0 && Nbest_Ch3 > 0) {
1055
1056 // for (int pair1 = 0; pair1 < 2; pair1++) {
1057 Int_t pair1 = 0;
1058
1059 if (pair1 == 1 && Nbest_Ch2 < 2)
1060 break;
1061 if (pair1 == 1 && Nbest_Ch3 < 2)
1062 break;
1063
1064 Float_t min_Chi2m = 100; // 400
1065 Float_t min_distX = 99;
1066 Float_t min_distY = 99;
1067 Float_t dAx23;
1068 Float_t dAy23;
1069
1070 for (Int_t bst = 0; bst < Nbest_Ch2; bst++) {
1071 if (pair1 == 1 && best_Ch2_gl[0] == bst)
1072 continue;
1073
1074 Float_t x2mid = par_ab_Ch2[0][ind_best_Ch2[bst]] * Zmid + par_ab_Ch2[1][ind_best_Ch2[bst]];
1075 Float_t y2mid = par_ab_Ch2[2][ind_best_Ch2[bst]] * Zmid + par_ab_Ch2[3][ind_best_Ch2[bst]];
1076
1077 for (Int_t bst3 = 0; bst3 < Nbest_Ch3; bst3++) {
1078 if (pair1 == 1 && best_Ch3_gl[0] == bst3)
1079 continue;
1080
1081 Float_t x3mid = par_ab_Ch3[0][ind_best_Ch3[bst3]] * Zmid3 + par_ab_Ch3[1][ind_best_Ch3[bst3]];
1082 Float_t y3mid = par_ab_Ch3[2][ind_best_Ch3[bst3]] * Zmid3 + par_ab_Ch3[3][ind_best_Ch3[bst3]];
1083
1084 dAx23 = par_ab_Ch2[0][ind_best_Ch2[bst]] - par_ab_Ch3[0][ind_best_Ch3[bst3]];
1085 dAy23 = par_ab_Ch2[2][ind_best_Ch2[bst]] - par_ab_Ch3[2][ind_best_Ch3[bst3]];
1086 min_distX = x2mid - x3mid; // min
1087 min_distY = y2mid - y3mid; // min
1088
1089 Float_t Chi2_m =
1090 (min_distX * min_distX / (sig_dx * sig_dx) + min_distY * min_distY / (sig_dy * sig_dy)
1091 + dAx23 * dAx23 / (sig_dax * sig_dax) + dAy23 * dAy23 / (sig_day * sig_day));
1092 Float_t DAx23;
1093 Float_t DAy23;
1094 Float_t Min_distX;
1095 Float_t Min_distY;
1096
1097 if (Chi2_m < min_Chi2m) {
1098 min_Chi2m = Chi2_m; // min
1099 best = bst;
1100 best3 = bst3;
1101
1102 Nbest_Ch23_gl = pair1 + 1;
1103
1104 best_Ch2_gl[pair1] = best;
1105 best_Ch3_gl[pair1] = best3;
1106 Chi2_match[pair1] = Chi2_m;
1107
1108 DAx23 = dAx23;
1109 DAy23 = dAy23;
1110 Min_distX = min_distX;
1111 Min_distY = min_distY;
1112
1113 } // Chi2_m
1114 } // bst3++
1115 } // bst++
1116
1117 if (Nbest_Ch23_gl == pair1 + 1) {
1118
1119 FillEfficiency(Wires_Ch2, Nhits_Ch2, ind_best_Ch2[best], Denom_Ch2, Nomin_Ch2, min_distX,
1120 min_distY); // ch2 8point
1121 FillEfficiency(Wires_Ch3, Nhits_Ch3, ind_best_Ch3[best3], Denom_Ch3, Nomin_Ch3, min_distX,
1122 min_distY); // ch3 8point
1123
1124 // min 9p (4+5)
1125 if (Nhits_Ch2[ind_best_Ch2[best]] + Nhits_Ch3[ind_best_Ch3[best3]] > 8) { // proviso 9p
1126
1127 FillEfficiency(Wires_Ch2, Nhits_Ch2, ind_best_Ch2[best], Denom_Ch2_9p, Nomin_Ch2_9p, min_distX,
1128 min_distY); // ch2 9point
1129 FillEfficiency(Wires_Ch3, Nhits_Ch3, ind_best_Ch3[best3], Denom_Ch3_9p, Nomin_Ch3_9p, min_distX,
1130 min_distY); // ch3 9point
1131
1132 } //>8
1133 } // Nbest_Ch23_gl == pair1+1
1134 // }//pair1++
1135
1136 } // Nbest_Ch2>0 Nbest_Ch3>0
1137
1138 if (Nbest_Ch23_gl > 0) {
1139
1140 Float_t Z_mid23 = 469.5; // cm
1141 Float_t Z_to_target = -469.5; // cm
1142 Float_t X23_v[2];
1143 Float_t Y23_v[2];
1144
1145 for (Int_t bst = 0; bst < Nbest_Ch23_gl; bst++) {
1146
1147 Int_t bst2 = best_Ch2_gl[bst];
1148 Int_t bst3 = best_Ch3_gl[bst];
1149
1150 Int_t best2 = ind_best_Ch2[bst2]; // it's index!
1151 Int_t best3 = ind_best_Ch3[bst3]; // it's index!
1152
1153 Float_t X2_v = par_ab_Ch2[0][best2] * (-Z_Ch2) + par_ab_Ch2[1][best2];
1154 Float_t Y2_v = par_ab_Ch2[2][best2] * (-Z_Ch2) + par_ab_Ch2[3][best2];
1155
1156 Float_t X3_v = par_ab_Ch3[0][best3] * (-Z_Ch3) + par_ab_Ch3[1][best3];
1157 Float_t Y3_v = par_ab_Ch3[2][best3] * (-Z_Ch3) + par_ab_Ch3[3][best3];
1158
1159 X23_v[bst] = 0.5 * (X2_v + X3_v);
1160 Y23_v[bst] = 0.5 * (Y2_v + Y3_v);
1161 /*
1162 //create track and put it into TClonesArray
1163 CbmTrack matchedTrack;
1164 FairTrackParam matchedParams;
1165 Float_t z0 = Z_mid23;
1166 Float_t x0 = par_ab_Ch23[1][bst]; //Check formula x0 - x coordinate of point in which line crosses
1167 mid-plane Float_t y0 = par_ab_Ch23[3][bst]; //Check formula y0 - y coordinate of point in which line
1168 crosses mid-plane matchedParams.SetPosition(TVector3(x0, y0, z0));
1169 matchedParams.SetTx(par_ab_Ch23[0][bst]);
1170 matchedParams.SetTy(par_ab_Ch23[2][bst]);
1171 matchedTrack.SetParamFirst(&matchedParams);
1172 matchedTrack.SetChiSq(Chi2_ndf_Ch23[bst]);
1173 matchedTrack.SetNDF(Nhits_Ch2[best2] + Nhits_Ch3[best3]); //WARNING!!! now it is number of hits
1174 in matched track! new((*mwpcMatchedTracks)[mwpcMatchedTracks->GetEntriesFast()]) CbmTrack(matchedTrack);
1175 */
1176 } // bst++
1177
1178 } // Nbest_Ch23_gl>0
1179
1180 tReco->Fill();
1181
1182 } // iev
1183
1184 Nomin_Ch2->Write();
1185 Nomin_Ch3->Write();
1186
1187 Denom_Ch2->Write();
1188 Denom_Ch3->Write();
1189
1190 Effic23_Ch2->Divide(Nomin_Ch2, Denom_Ch2, 1, 1);
1191 Effic23_Ch3->Divide(Nomin_Ch3, Denom_Ch3, 1, 1);
1192
1193 Effic23_Ch2->Write();
1194 Effic23_Ch3->Write();
1195
1196 Nomin_Ch2_9p->Write();
1197 Nomin_Ch3_9p->Write();
1198
1199 Denom_Ch2_9p->Write();
1200 Denom_Ch3_9p->Write();
1201
1202 Effic23_Ch2_9p->Divide(Nomin_Ch2_9p, Denom_Ch2_9p, 1, 1);
1203 Effic23_Ch3_9p->Divide(Nomin_Ch3_9p, Denom_Ch3_9p, 1, 1);
1204
1205 Effic23_Ch2_9p->Write();
1206 Effic23_Ch3_9p->Write();
1207
1208 for (int i1 = 0; i1 < 6; i1++) {
1209 hResid_Ch1[i1]->Write();
1210 hResid_Ch2[i1]->Write();
1211 hResid_Ch3[i1]->Write();
1212 hResid23_Ch2[i1]->Write();
1213 hResid23_Ch3[i1]->Write();
1214 }
1215
1216 tReco->Write();
1217 fReco->Close();
1218 return 0;
1219}
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
Double_t
double Double_t
gROOT
#define gROOT
BmnMwpcDigit::GetPlane
Short_t GetPlane() const
Definition BmnMwpcDigit.h:26
BmnMwpcDigit::GetTime
UInt_t GetTime() const
Definition BmnMwpcDigit.h:28
BmnMwpcDigit::GetWireNumber
Int_t GetWireNumber() const
Definition BmnMwpcDigit.h:27
TBranch::SetAddress
virtual void SetAddress(void *add)
TClonesArray::Clear
void Clear(Option_t *option="") override
TDirectoryFile::Get
T * Get(const char *namecycle)
TFile::Close
void Close(Option_t *option="") override
TH1::Fill
virtual Int_t Fill(const char *name, Double_t w)
TH1::Divide
virtual Bool_t Divide(const TH1 *h1)
TObjArray::GetEntriesFast
Int_t GetEntriesFast() const
TObjArray::GetEntries
Int_t GetEntries() const override
TObjArray::At
TObject * At(Int_t idx) const override
TObject::Write
virtual Int_t Write(const char *name=nullptr, Int_t option=0, Int_t bufsize=0)
TTree::Fill
virtual Int_t Fill()
TTree::Write
Int_t Write(const char *name=nullptr, Int_t option=0, Int_t bufsize=0) const override
TTree::Branch
virtual Int_t Branch(const char *folder, Int_t bufsize=32000, Int_t splitlevel=99)
abs
RVec< PromoteType< T > > abs(const RVec< T > &v)
sqrt
RVec< PromoteType< T > > sqrt(const RVec< T > &v)
dw
const Float_t dw
Definition read_mwpc1_688.C:21
FillFreeCoefVector
void FillFreeCoefVector(Double_t F[4], Float_t XVU_Ch[nPlanes][big], Int_t segIdx, Float_t z[nPlanes], Float_t sigm2[nPlanes], Int_t h[nPlanes])
Definition read_mwpc1_688.C:254
Z_Ch3
const Float_t Z_Ch3
Definition read_mwpc1_688.C:29
Z_Ch1
const Float_t Z_Ch1
Definition read_mwpc1_688.C:27
InverseMatrix
void InverseMatrix(Double_t A[4][4], Double_t b[4][4])
Definition read_mwpc1_688.C:282
sq3
const Double_t sq3
Definition read_mwpc1_688.C:22
Z_Ch2
const Float_t Z_Ch2
Definition read_mwpc1_688.C:28
read_mwpc1_688
Int_t read_mwpc1_688()
Definition read_mwpc1_688.C:630
FillFitMatrix
void FillFitMatrix(Double_t A[4][4], Float_t z[nPlanes], Float_t sigm2[nPlanes], Int_t h[nPlanes])
Definition read_mwpc1_688.C:215
nWires
const Int_t nWires
Definition read_mwpc1_688.C:18
SegmentFinder
void SegmentFinder(Int_t chNum, Int_t wires_Ch[nPlanes][big], Float_t XVU_Ch[nPlanes][big], Int_t Nhits_Ch[big], Int_t iw_Ch[nPlanes], Int_t &Nseg, Int_t wires_glob[nWires][nPlanes], Float_t xuv_glob[nWires][nPlanes], Int_t minHits, Short_t code)
Definition read_mwpc1_688.C:31
Zmid3
const Float_t Zmid3
Definition read_mwpc1_688.C:26
Chi2_Max
const Double_t Chi2_Max
Definition read_mwpc1_688.C:23
ProcessSegments
void ProcessSegments(Int_t chNum, Float_t sigma, Float_t z_loc[nPlanes], Int_t Min_hits, Int_t &Nseg, Int_t Nhits_Ch[big], Int_t Wires_Ch[nPlanes][big], Float_t XVU_Ch[nPlanes][big], Int_t &Nbest_Ch, Int_t ind_best_Ch[5], Double_t Chi2_ndf_Ch[big], Double_t Chi2_ndf_best_Ch[5], Double_t par_ab_Ch[4][big])
Definition read_mwpc1_688.C:339
FillEfficiency
void FillEfficiency(Int_t Wires_Ch[nPlanes][big], Int_t Nhits_Ch[big], Int_t segIdx, TH1F *Denom_Ch, TH1F *Nomin_Ch, Int_t min_distX, Int_t min_distY)
Definition read_mwpc1_688.C:580
big
const Int_t big
Definition read_mwpc1_688.C:20
FillResid
void FillResid(Int_t Wires_Ch[nPlanes][big], Int_t segIdx, Float_t XVU_Ch[nPlanes][big], Double_t par_ab_Ch[4][big], Float_t z_loc[nPlanes], TH1F *hResid_Ch[nPlanes])
Definition read_mwpc1_688.C:606
Zmid
const Float_t Zmid
Definition read_mwpc1_688.C:25
nPlanes
const Int_t nPlanes
Definition read_mwpc1_688.C:19
TMath::Sqrt
Double_t Sqrt(Double_t x)
std
Definition json.hpp:4446