BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
BmnDchTrackFinder.cxx
Go to the documentation of this file.
1// @(#)bmnroot/dch:$Id$
2// Author: Pavel Batyuk (VBLHEP) <pavel.batyuk@jinr.ru> 2017-01-21
3
4////////////////////////////////////////////////////////////////////////////////
5// //
6// BmnDchTrackFinder //
7// //
8// Implementation of an algorithm developed by //
9// N.Voytishin and V.Paltchik (LIT) //
10// to the BmnRoot software //
11// //
12// The algorithm serves for searching for track segments //
13// in the Drift Chambers of the BM@N experiment //
14// //
15////////////////////////////////////////////////////////////////////////////////
16
17#include "BmnDchTrackFinder.h"
18
19#include "FairLogger.h"
20#include "FairMCPoint.h"
21#include "FairRootManager.h"
22#include "FitWLSQ.h"
23#include "TFile.h"
24#include "TH1F.h"
25#include "TH2F.h"
26#include "TProfile.h"
27
28#include <TMath.h>
29#include <cmath>
30
31static Double_t workTime = 0.0;
32
33BmnDchTrackFinder::BmnDchTrackFinder(Int_t period, Int_t number, Bool_t isExp)
34 : fPeriod(period)
35 , fRunId(number)
36 , expData(isExp)
37 , fhTestFlnm("test.BmnDCHTracking.root")
38{
39 fEventNo = 0;
40 N = (fPeriod == 7 && fRunId > 3589) ? 100 : 15; // needs to be adjusted according to hit multiplicity
41 tracksDch = "BmnDchTrack";
42 if (expData) {
43 InputDigitsBranchName = "DCH";
44 fTransferFunctionName = "transfer_func.txt";
45
46 prev_wire = -1;
47 prev_time = -1;
48 } else {
49 InputDigitsBranchName = "BmnDchHit";
50 }
51
52 layers_with = 0;
53 layers_with2 = 0;
54
55 prev_wire = -1;
56 prev_time = -1;
57
58 nChambers = 2;
59 nWires = 4;
60 nLayers = 2;
61 nSegmentsMax = 100;
62
63 Z_dch1 = 510.2;
64 Z_dch2 = 709.5;
65 Z_dch_mid = (Z_dch1 + Z_dch2) / 2.;
66 dZ_dch_mid = Z_dch2 - Z_dch_mid; // > 0
67 dZ_dch = Z_dch2 - Z_dch1;
68
69 // Some alignment corrections
70 if (fPeriod == 7) {
71
72 x1_sh = 0; // 10.3;//0; //11.2;//-8.56;//run6 10.9;// 5.34;//5.14;//10.31; //
73 y1_sh = 0; //-3;//0; //1.76;//5.76;//run6 -1.19;//0;//-1.03;
74 x2_sh = -5.99;
75 ; //-5.4;//4.7;//8.60; //9.-3.5;//0.0;//run6 5.56;//0;//7.95;// //
76 y2_sh = -.24; //-.47;//-3.3;//-12.40; //1.8;//-7.1;//run6 -2.47;//-1.38;//-3.7;
77
78 x1_slope_sh = 0; //-.070;//0;
79 y1_slope_sh = 0; //.060;//0;
80 x2_slope_sh = -.005; //-.034;//-.071;//0.;
81 y2_slope_sh = .005; //.014;//.061;//0.007;
82 }
83
84 if (!expData) { // no alignment shifts for MC data
85
86 x1_sh = 0.;
87 y1_sh = 0.;
88 x2_sh = 0.;
89 y2_sh = 0.;
90
91 x1_slope_sh = 0.;
92 y1_slope_sh = 0.;
93 x2_slope_sh = 0.;
94 y2_slope_sh = 0.;
95 }
96 scale = 0.5;
97 // cout<<" x2 sh "<<x2_sh<<endl;
98}
99
100void BmnDchTrackFinder::initHists()
101{
102 // added temporary hists
103 hNLayWithFiredWiresDc1 = new TH1F("hNLayWithFiredWiresDc1", "layers with fired wires in dc1", 9, 0, 9);
104 hNLayWithFiredWiresDc2 = new TH1F("hNLayWithFiredWiresDc2", "layers with fired wires in dc2", 9, 0, 9);
105 hNoFiredWiresOnLayerDc1 = new TH1F("hNoFiredWiresOnLayerDc1", "no fired wires on a layer in dc1", 8, 0, 8);
106 hNoFiredWiresOnLayerDc2 = new TH1F("hNoFiredWiresOnLayerDc2", "no fired wires on a layer in dc2", 8, 0, 8);
107 hNomin1 = new TH1F("hNomin1", " eff nominator dc1", 8, 0, 8);
108 hDenom1 = new TH1F("hDenom1", " eff denominator dc1", 8, 0, 8);
109 hNomin2 = new TH1F("hNomin2", " eff nominator dc2", 8, 0, 8);
110 hDenom2 = new TH1F("hDenom2", " eff denominator dc2", 8, 0, 8);
111
112 hEff1 = new TH1F("hEff1", " eff dc1", 8, 0, 8);
113 hEff2 = new TH1F("hEff2", " eff dc2", 8, 0, 8);
114 // hNomin1->Sumw2();
115 // hDenom1->Sumw2();
116 // hNomin2->Sumw2();
117 // hDenom2->Sumw2();
118 hEff1->Sumw2();
119 hEff2->Sumw2();
120
121 hXDC1_atZ0 = new TH1F("hXDC1_atZ0", " XDC1 extrapol Z = 0", 200, -100, 100);
122 hYDC1_atZ0 = new TH1F("hYDC1_atZ0", " YDC1 extrapol Z = 0", 200, -100, 100);
123 hXDC2_atZ0 = new TH1F("hXDC2_atZ0", " XDC2 extrapol Z = 0", 200, -100, 100);
124 hYDC2_atZ0 = new TH1F("hYDC2_atZ0", " YDC2 extrapol Z = 0", 200, -100, 100);
125
126 hSlot_1xa_time8p = new TH1F("hSlot_1xa_time8p", "time8ps_for_plane_DC1_xa", 2000, 0, 2000);
127 hSlot_1xb_time8p = new TH1F("hSlot_1xb_time8p", "time8ps_for_plane_DC1_xb", 2000, 0, 2000);
128 hSlot_1ya_time8p = new TH1F("hSlot_1ya_time8p", "time8ps_for_plane_DC1_ya", 2000, 0, 2000);
129 hSlot_1yb_time8p = new TH1F("hSlot_1yb_time8p", "time8ps_for_plane_DC1_yb", 2000, 0, 2000);
130 hSlot_1ua_time8p = new TH1F("hSlot_1ua_time8p", "time8ps_for_plane_DC1_ua", 2000, 0, 2000);
131 hSlot_1ub_time8p = new TH1F("hSlot_1ub_time8p", "time8ps_for_plane_DC1_ub", 2000, 0, 2000);
132 hSlot_1va_time8p = new TH1F("hSlot_1va_time8p", "time8ps_for_plane_DC1_va", 2000, 0, 2000);
133 hSlot_1vb_time8p = new TH1F("hSlot_1vb_time8p", "time8ps_for_plane_DC1_vb", 2000, 0, 2000);
134 hSlot_2xa_time8p = new TH1F("hSlot_2xa_time8p", "time8ps_for_plane_DC2_xa", 2000, 0, 2000);
135 hSlot_2xb_time8p = new TH1F("hSlot_2xb_time8p", "time8ps_for_plane_DC2_xb", 2000, 0, 2000);
136 hSlot_2ya_time8p = new TH1F("hSlot_2ya_time8p", "time8ps_for_plane_DC2_ya", 2000, 0, 2000);
137 hSlot_2yb_time8p = new TH1F("hSlot_2yb_time8p", "time8ps_for_plane_DC2_yb", 2000, 0, 2000);
138 hSlot_2ua_time8p = new TH1F("hSlot_2ua_time8p", "time8ps_for_plane_DC2_ua", 2000, 0, 2000);
139 hSlot_2ub_time8p = new TH1F("hSlot_2ub_time8p", "time8ps_for_plane_DC2_ub", 2000, 0, 2000);
140 hSlot_2va_time8p = new TH1F("hSlot_2va_time8p", "time8ps_for_plane_DC2_va", 2000, 0, 2000);
141 hSlot_2vb_time8p = new TH1F("hSlot_2vb_time8p", "time8ps_for_plane_DC2_vb", 2000, 0, 2000);
142
143 hResidx1a = new TH1F("hResidx1a", "dc1 xa_mes - xa_from_fit", 500, -0.1, 0.1);
144 hResidu1a = new TH1F("hResidu1a", "dc1 ua_mes - ua_from_fit", 500, -0.1, 0.1);
145 hResidy1a = new TH1F("hResidy1a", "dc1 ya_mes - ya_from_fit", 500, -0.1, 0.1);
146 hResidv1a = new TH1F("hResidv1a", "dc1 va_mes - va_from_fit", 500, -0.1, 0.1);
147 hResidx2a = new TH1F("hResidx2a", "dc2 xa_mes - xa_from_fit", 500, -0.1, 0.1);
148 hResidu2a = new TH1F("hResidu2a", "dc2 ua_mes - ua_from_fit", 500, -0.1, 0.1);
149 hResidy2a = new TH1F("hResidy2a", "dc2 ya_mes - ya_from_fit", 500, -0.1, 0.1);
150 hResidv2a = new TH1F("hResidv2a", "dc2 va_mes - va_from_fit", 500, -0.1, 0.1);
151
152 hResidx1b = new TH1F("hResidx1b", "dc1 xb_mes - xb_from_fit", 500, -0.1, 0.1);
153 hResidu1b = new TH1F("hResidu1b", "dc1 ub_mes - ub_from_fit", 500, -0.1, 0.1);
154 hResidy1b = new TH1F("hResidy1b", "dc1 yb_mes - yb_from_fit", 500, -0.1, 0.1);
155 hResidv1b = new TH1F("hResidv1b", "dc1 vb_mes - vb_from_fit", 500, -0.1, 0.1);
156 hResidx2b = new TH1F("hResidx2b", "dc2 xb_mes - xb_from_fit", 500, -0.1, 0.1);
157 hResidu2b = new TH1F("hResidu2b", "dc2 ub_mes - ub_from_fit", 500, -0.1, 0.1);
158 hResidy2b = new TH1F("hResidy2b", "dc2 yb_mes - yb_from_fit", 500, -0.1, 0.1);
159 hResidv2b = new TH1F("hResidv2b", "dc2 vb_mes - vb_from_fit", 500, -0.1, 0.1);
160
161 hResidx1a_0p1 = new TH1F(
162 "hResidx1a_0p1", "dc1 xa_mes - x_from_fit(w/o x in fit) for hit at (0; 0.1) away from wire", 500, -0.1, 0.1);
163 hResidu1a_0p1 = new TH1F(
164 "hResidu1a_0p1", "dc1 ua_mes - u_from_fit(w/o u in fit) for hit at (0; 0.1) away from wire", 500, -0.1, 0.1);
165 hResidy1a_0p1 = new TH1F(
166 "hResidy1a_0p1", "dc1 ya_mes - y_from_fit(w/o y in fit) for hit at (0; 0.1) away from wire", 500, -0.1, 0.1);
167 hResidv1a_0p1 = new TH1F(
168 "hResidv1a_0p1", "dc1 va_mes - v_from_fit(w/o v in fit) for hit at (0; 0.1) away from wire", 500, -0.1, 0.1);
169 hResidx2a_0p1 = new TH1F(
170 "hResidx2a_0p1", "dc2 xa_mes - x_from_fit(w/o x in fit) for hit at (0; 0.1) away from wire", 500, -0.1, 0.1);
171 hResidu2a_0p1 = new TH1F(
172 "hResidu2a_0p1", "dc2 ua_mes - u_from_fit(w/o u in fit) for hit at (0; 0.1) away from wire", 500, -0.1, 0.1);
173 hResidy2a_0p1 = new TH1F(
174 "hResidy2a_0p1", "dc2 ya_mes - y_from_fit(w/o y in fit) for hit at (0; 0.1) away from wire", 500, -0.1, 0.1);
175 hResidv2a_0p1 = new TH1F(
176 "hResidv2a_0p1", "dc2 va_mes - v_from_fit(w/o v in fit) for hit at (0; 0.1) away from wire", 500, -0.1, 0.1);
177
178 hResidx1a_0p1_0p4 =
179 new TH1F("hResidx1a_0p1_0p4", "dc1 xa_mes - x_from_fit(w/o x in fit) for hit at (0.1; 0.4) away from wire", 500,
180 -0.1, 0.1);
181 hResidu1a_0p1_0p4 =
182 new TH1F("hResidu1a_0p1_0p4", "dc1 ua_mes - u_from_fit(w/o u in fit) for hit at (0.1; 0.4) away from wire", 500,
183 -0.1, 0.1);
184 hResidy1a_0p1_0p4 =
185 new TH1F("hResidy1a_0p1_0p4", "dc1 ya_mes - y_from_fit(w/o y in fit) for hit at (0.1; 0.4) away from wire", 500,
186 -0.1, 0.1);
187 hResidv1a_0p1_0p4 =
188 new TH1F("hResidv1a_0p1_0p4", "dc1 va_mes - v_from_fit(w/o v in fit) for hit at (0.1; 0.4) away from wire", 500,
189 -0.1, 0.1);
190 hResidx2a_0p1_0p4 =
191 new TH1F("hResidx2a_0p1_0p4", "dc2 xa_mes - x_from_fit(w/o x in fit) for hit at (0.1; 0.4) away from wire", 500,
192 -0.1, 0.1);
193 hResidu2a_0p1_0p4 =
194 new TH1F("hResidu2a_0p1_0p4", "dc2 ua_mes - u_from_fit(w/o u in fit) for hit at (0.1; 0.4) away from wire", 500,
195 -0.1, 0.1);
196 hResidy2a_0p1_0p4 =
197 new TH1F("hResidy2a_0p1_0p4", "dc2 ya_mes - y_from_fit(w/o y in fit) for hit at (0.1; 0.4) away from wire", 500,
198 -0.1, 0.1);
199 hResidv2a_0p1_0p4 =
200 new TH1F("hResidv2a_0p1_0p4", "dc2 va_mes - v_from_fit(w/o v in fit) for hit at (0.1; 0.4) away from wire", 500,
201 -0.1, 0.1);
202
203 hResidx1a_0p4_0p5 =
204 new TH1F("hResidx1a_0p4_0p5", "dc1 xa_mes - x_from_fit(w/o x in fit) for hit at (0.4; 0.5) away from wire", 500,
205 -0.1, 0.1);
206 hResidu1a_0p4_0p5 =
207 new TH1F("hResidu1a_0p4_0p5", "dc1 ua_mes - u_from_fit(w/o u in fit) for hit at (0.4; 0.5) away from wire", 500,
208 -0.1, 0.1);
209 hResidy1a_0p4_0p5 =
210 new TH1F("hResidy1a_0p4_0p5", "dc1 ya_mes - y_from_fit(w/o y in fit) for hit at (0.4; 0.5) away from wire", 500,
211 -0.1, 0.1);
212 hResidv1a_0p4_0p5 =
213 new TH1F("hResidv1a_0p4_0p5", "dc1 va_mes - v_from_fit(w/o v in fit) for hit at (0.4; 0.5) away from wire", 500,
214 -0.1, 0.1);
215 hResidx2a_0p4_0p5 =
216 new TH1F("hResidx2a_0p4_0p5", "dc2 xa_mes - x_from_fit(w/o x in fit) for hit at (0.4; 0.5) away from wire", 500,
217 -0.1, 0.1);
218 hResidu2a_0p4_0p5 =
219 new TH1F("hResidu2a_0p4_0p5", "dc2 ua_mes - u_from_fit(w/o u in fit) for hit at (0.4; 0.5) away from wire", 500,
220 -0.1, 0.1);
221 hResidy2a_0p4_0p5 =
222 new TH1F("hResidy2a_0p4_0p5", "dc2 ya_mes - y_from_fit(w/o y in fit) for hit at (0.4; 0.5) away from wire", 500,
223 -0.1, 0.1);
224 hResidv2a_0p4_0p5 =
225 new TH1F("hResidv2a_0p4_0p5", "dc2 va_mes - v_from_fit(w/o v in fit) for hit at (0.4; 0.5) away from wire", 500,
226 -0.1, 0.1);
227
228 hResidx1b_0p1 = new TH1F(
229 "hResidx1b_0p1", "dc1 xb_mes - x_from_fit(w/o x in fit) for hit at (0; 0.1) away from wire", 500, -0.1, 0.1);
230 hResidu1b_0p1 = new TH1F(
231 "hResidu1b_0p1", "dc1 ub_mes - u_from_fit(w/o u in fit) for hit at (0; 0.1) away from wire", 500, -0.1, 0.1);
232 hResidy1b_0p1 = new TH1F(
233 "hResidy1b_0p1", "dc1 yb_mes - y_from_fit(w/o y in fit) for hit at (0; 0.1) away from wire", 500, -0.1, 0.1);
234 hResidv1b_0p1 = new TH1F(
235 "hResidv1b_0p1", "dc1 vb_mes - v_from_fit(w/o v in fit) for hit at (0; 0.1) away from wire", 500, -0.1, 0.1);
236 hResidx2b_0p1 = new TH1F(
237 "hResidx2b_0p1", "dc2 xb_mes - x_from_fit(w/o x in fit) for hit at (0; 0.1) away from wire", 500, -0.1, 0.1);
238 hResidu2b_0p1 = new TH1F(
239 "hResidu2b_0p1", "dc2 ub_mes - u_from_fit(w/o u in fit) for hit at (0; 0.1) away from wire", 500, -0.1, 0.1);
240 hResidy2b_0p1 = new TH1F(
241 "hResidy2b_0p1", "dc2 yb_mes - y_from_fit(w/o y in fit) for hit at (0; 0.1) away from wire", 500, -0.1, 0.1);
242 hResidv2b_0p1 = new TH1F(
243 "hResidv2b_0p1", "dc2 vb_mes - v_from_fit(w/o v in fit) for hit at (0; 0.1) away from wire", 500, -0.1, 0.1);
244
245 hResidx1b_0p1_0p4 =
246 new TH1F("hResidx1b_0p1_0p4", "dc1 xb_mes - x_from_fit(w/o x in fit) for hit at (0.1; 0.4) away from wire", 500,
247 -0.1, 0.1);
248 hResidu1b_0p1_0p4 =
249 new TH1F("hResidu1b_0p1_0p4", "dc1 ub_mes - u_from_fit(w/o u in fit) for hit at (0.1; 0.4) away from wire", 500,
250 -0.1, 0.1);
251 hResidy1b_0p1_0p4 =
252 new TH1F("hResidy1b_0p1_0p4", "dc1 yb_mes - y_from_fit(w/o y in fit) for hit at (0.1; 0.4) away from wire", 500,
253 -0.1, 0.1);
254 hResidv1b_0p1_0p4 =
255 new TH1F("hResidv1b_0p1_0p4", "dc1 vb_mes - v_from_fit(w/o v in fit) for hit at (0.1; 0.4) away from wire", 500,
256 -0.1, 0.1);
257 hResidx2b_0p1_0p4 =
258 new TH1F("hResidx2b_0p1_0p4", "dc2 xb_mes - x_from_fit(w/o x in fit) for hit at (0.1; 0.4) away from wire", 500,
259 -0.1, 0.1);
260 hResidu2b_0p1_0p4 =
261 new TH1F("hResidu2b_0p1_0p4", "dc2 ub_mes - u_from_fit(w/o u in fit) for hit at (0.1; 0.4) away from wire", 500,
262 -0.1, 0.1);
263 hResidy2b_0p1_0p4 =
264 new TH1F("hResidy2b_0p1_0p4", "dc2 yb_mes - y_from_fit(w/o y in fit) for hit at (0.1; 0.4) away from wire", 500,
265 -0.1, 0.1);
266 hResidv2b_0p1_0p4 =
267 new TH1F("hResidv2b_0p1_0p4", "dc2 vb_mes - v_from_fit(w/o v in fit) for hit at (0.1; 0.4) away from wire", 500,
268 -0.1, 0.1);
269
270 hResidx1b_0p4_0p5 =
271 new TH1F("hResidx1b_0p4_0p5", "dc1 xb_mes - x_from_fit(w/o x in fit) for hit at (0.4; 0.5) away from wire", 500,
272 -0.1, 0.1);
273 hResidu1b_0p4_0p5 =
274 new TH1F("hResidu1b_0p4_0p5", "dc1 ub_mes - u_from_fit(w/o u in fit) for hit at (0.4; 0.5) away from wire", 500,
275 -0.1, 0.1);
276 hResidy1b_0p4_0p5 =
277 new TH1F("hResidy1b_0p4_0p5", "dc1 yb_mes - y_from_fit(w/o y in fit) for hit at (0.4; 0.5) away from wire", 500,
278 -0.1, 0.1);
279 hResidv1b_0p4_0p5 =
280 new TH1F("hResidv1b_0p4_0p5", "dc1 vb_mes - v_from_fit(w/o v in fit) for hit at (0.4; 0.5) away from wire", 500,
281 -0.1, 0.1);
282 hResidx2b_0p4_0p5 =
283 new TH1F("hResidx2b_0p4_0p5", "dc2 xb_mes - x_from_fit(w/o x in fit) for hit at (0.4; 0.5) away from wire", 500,
284 -0.1, 0.1);
285 hResidu2b_0p4_0p5 =
286 new TH1F("hResidu2b_0p4_0p5", "dc2 ub_mes - u_from_fit(w/o u in fit) for hit at (0.4; 0.5) away from wire", 500,
287 -0.1, 0.1);
288 hResidy2b_0p4_0p5 =
289 new TH1F("hResidy2b_0p4_0p5", "dc2 yb_mes - y_from_fit(w/o y in fit) for hit at (0.4; 0.5) away from wire", 500,
290 -0.1, 0.1);
291 hResidv2b_0p4_0p5 =
292 new TH1F("hResidv2b_0p4_0p5", "dc2 vb_mes - v_from_fit(w/o v in fit) for hit at (0.4; 0.5) away from wire", 500,
293 -0.1, 0.1);
294
295 // pulls
296
297 hPullsx1a_0p1 = new TH1F("hPullsx1a_0p1", "dc1 xa pull (0 - 0.1) away from wire", 100, -3, 3);
298 hPullsu1a_0p1 = new TH1F("hPullsu1a_0p1", "dc1 ua pull (0 - 0.1) away from wire", 100, -3, 3);
299 hPullsy1a_0p1 = new TH1F("hPullsy1a_0p1", "dc1 ya pull (0 - 0.1) away from wire", 100, -3, 3);
300 hPullsv1a_0p1 = new TH1F("hPullsv1a_0p1", "dc1 va pull (0 - 0.1) away from wire", 100, -3, 3);
301 hPullsx2a_0p1 = new TH1F("hPullsx2a_0p1", "dc2 xa pull (0 - 0.1) away from wire", 100, -3, 3);
302 hPullsu2a_0p1 = new TH1F("hPullsu2a_0p1", "dc2 ua pull (0 - 0.1) away from wire", 100, -3, 3);
303 hPullsy2a_0p1 = new TH1F("hPullsy2a_0p1", "dc2 ya pull (0 - 0.1) away from wire", 100, -3, 3);
304 hPullsv2a_0p1 = new TH1F("hPullsv2a_0p1", "dc2 va pull (0 - 0.1) away from wire", 100, -3, 3);
305
306 hPullsx1a_0p1_0p4 = new TH1F("hPullsx1a_0p1_0p4", "dc1 xa pull (0.1 - 0.4) away from wire", 100, -3, 3);
307 hPullsu1a_0p1_0p4 = new TH1F("hPullsu1a_0p1_0p4", "dc1 ua pull (0.1 - 0.4) away from wire", 100, -3, 3);
308 hPullsy1a_0p1_0p4 = new TH1F("hPullsy1a_0p1_0p4", "dc1 ya pull (0.1 - 0.4) away from wire", 100, -3, 3);
309 hPullsv1a_0p1_0p4 = new TH1F("hPullsv1a_0p1_0p4", "dc1 va pull (0.1 - 0.4) away from wire", 100, -3, 3);
310 hPullsx2a_0p1_0p4 = new TH1F("hPullsx2a_0p1_0p4", "dc2 xa pull (0.1 - 0.4) away from wire", 100, -3, 3);
311 hPullsu2a_0p1_0p4 = new TH1F("hPullsu2a_0p1_0p4", "dc2 ua pull (0.1 - 0.4) away from wire", 100, -3, 3);
312 hPullsy2a_0p1_0p4 = new TH1F("hPullsy2a_0p1_0p4", "dc2 ya pull (0.1 - 0.4) away from wire", 100, -3, 3);
313 hPullsv2a_0p1_0p4 = new TH1F("hPullsv2a_0p1_0p4", "dc2 va pull (0.1 - 0.4) away from wire", 100, -3, 3);
314
315 hPullsx1a_0p4_0p5 = new TH1F("hPullsx1a_0p4_0p5", "dc1 xa pull (0.4 - 0.5) away from wire", 100, -3, 3);
316 hPullsu1a_0p4_0p5 = new TH1F("hPullsu1a_0p4_0p5", "dc1 ua pull (0.4 - 0.5) away from wire", 100, -3, 3);
317 hPullsy1a_0p4_0p5 = new TH1F("hPullsy1a_0p4_0p5", "dc1 ya pull (0.4 - 0.5) away from wire", 100, -3, 3);
318 hPullsv1a_0p4_0p5 = new TH1F("hPullsv1a_0p4_0p5", "dc1 va pull (0.4 - 0.5) away from wire", 100, -3, 3);
319 hPullsx2a_0p4_0p5 = new TH1F("hPullsx2a_0p4_0p5", "dc2 xa pull (0.4 - 0.5) away from wire", 100, -3, 3);
320 hPullsu2a_0p4_0p5 = new TH1F("hPullsu2a_0p4_0p5", "dc2 ua pull (0.4 - 0.5) away from wire", 100, -3, 3);
321 hPullsy2a_0p4_0p5 = new TH1F("hPullsy2a_0p4_0p5", "dc2 ya pull (0.4 - 0.5) away from wire", 100, -3, 3);
322 hPullsv2a_0p4_0p5 = new TH1F("hPullsv2a_0p4_0p5", "dc2 va pull (0.4 - 0.5) away from wire", 100, -3, 3);
323
324 hPullsx1b_0p1 = new TH1F("hPullsx1b_0p1", "dc1 xb pull (0 - 0.1) away from wire", 100, -3, 3);
325 hPullsu1b_0p1 = new TH1F("hPullsu1b_0p1", "dc1 ub pull (0 - 0.1) away from wire", 100, -3, 3);
326 hPullsy1b_0p1 = new TH1F("hPullsy1b_0p1", "dc1 yb pull (0 - 0.1) away from wire", 100, -3, 3);
327 hPullsv1b_0p1 = new TH1F("hPullsv1b_0p1", "dc1 vb pull (0 - 0.1) away from wire", 100, -3, 3);
328 hPullsx2b_0p1 = new TH1F("hPullsx2b_0p1", "dc2 xb pull (0 - 0.1) away from wire", 100, -3, 3);
329 hPullsu2b_0p1 = new TH1F("hPullsu2b_0p1", "dc2 ub pull (0 - 0.1) away from wire", 100, -3, 3);
330 hPullsy2b_0p1 = new TH1F("hPullsy2b_0p1", "dc2 yb pull (0 - 0.1) away from wire", 100, -3, 3);
331 hPullsv2b_0p1 = new TH1F("hPullsv2b_0p1", "dc2 vb pull (0 - 0.1) away from wire", 100, -3, 3);
332
333 hPullsx1b_0p1_0p4 = new TH1F("hPullsx1b_0p1_0p4", "dc1 xb pull (0.1 - 0.4) away from wire", 100, -3, 3);
334 hPullsu1b_0p1_0p4 = new TH1F("hPullsu1b_0p1_0p4", "dc1 ub pull (0.1 - 0.4) away from wire", 100, -3, 3);
335 hPullsy1b_0p1_0p4 = new TH1F("hPullsy1b_0p1_0p4", "dc1 yb pull (0.1 - 0.4) away from wire", 100, -3, 3);
336 hPullsv1b_0p1_0p4 = new TH1F("hPullsv1b_0p1_0p4", "dc1 vb pull (0.1 - 0.4) away from wire", 100, -3, 3);
337 hPullsx2b_0p1_0p4 = new TH1F("hPullsx2b_0p1_0p4", "dc2 xb pull (0.1 - 0.4) away from wire", 100, -3, 3);
338 hPullsu2b_0p1_0p4 = new TH1F("hPullsu2b_0p1_0p4", "dc2 ub pull (0.1 - 0.4) away from wire", 100, -3, 3);
339 hPullsy2b_0p1_0p4 = new TH1F("hPullsy2b_0p1_0p4", "dc2 yb pull (0.1 - 0.4) away from wire", 100, -3, 3);
340 hPullsv2b_0p1_0p4 = new TH1F("hPullsv2b_0p1_0p4", "dc2 vb pull (0.1 - 0.4) away from wire", 100, -3, 3);
341
342 hPullsx1b_0p4_0p5 = new TH1F("hPullsx1b_0p4_0p5", "dc1 xb pull (0.4 - 0.5) away from wire", 100, -3, 3);
343 hPullsu1b_0p4_0p5 = new TH1F("hPullsu1b_0p4_0p5", "dc1 ub pull (0.4 - 0.5) away from wire", 100, -3, 3);
344 hPullsy1b_0p4_0p5 = new TH1F("hPullsy1b_0p4_0p5", "dc1 yb pull (0.4 - 0.5) away from wire", 100, -3, 3);
345 hPullsv1b_0p4_0p5 = new TH1F("hPullsv1b_0p4_0p5", "dc1 vb pull (0.4 - 0.5) away from wire", 100, -3, 3);
346 hPullsx2b_0p4_0p5 = new TH1F("hPullsx2b_0p4_0p5", "dc2 xb pull (0.4 - 0.5) away from wire", 100, -3, 3);
347 hPullsu2b_0p4_0p5 = new TH1F("hPullsu2b_0p4_0p5", "dc2 ub pull (0.4 - 0.5) away from wire", 100, -3, 3);
348 hPullsy2b_0p4_0p5 = new TH1F("hPullsy2b_0p4_0p5", "dc2 yb pull (0.4 - 0.5) away from wire", 100, -3, 3);
349 hPullsv2b_0p4_0p5 = new TH1F("hPullsv2b_0p4_0p5", "dc2 vb pull (0.4 - 0.5) away from wire", 100, -3, 3);
350
351 hSlot_1xa_time = new TH1F("hSlot_1xa_time", "times_for_plane_DC1_xa", 2000, 0, 2000);
352 hSlot_1xb_time = new TH1F("hSlot_1xb_time", "times_for_plane_DC1_xb", 2000, 0, 2000);
353 Ntrack1 = new TH1D("Ntrack1", "N segs in dch1", 20, 0, 20);
354 Ntrack2 = new TH1D("Ntrack2", "N segs in dch2", 20, 0, 20);
355 NGtracks = new TH1D("NGtracks", "N global dch tracks per event", 10, 0, 10);
356 hXa_wireOccupancy = new TH1F("hXa_wireOccupancy", " xa wire nr ", 250, 0, 250);
357 hXb_wireOccupancy = new TH1F("hXb_wireOccupancy", " xb wire nr ", 250, 0, 250);
358 hYa_wireOccupancy = new TH1F("hYa_wireOccupancy", " ya wire nr ", 250, 0, 250);
359 hYb_wireOccupancy = new TH1F("hYb_wireOccupancy", " yb wire nr ", 250, 0, 250);
360 hUa_wireOccupancy = new TH1F("hUa_wireOccupancy", " ua wire nr ", 250, 0, 250);
361 hUb_wireOccupancy = new TH1F("hUb_wireOccupancy", " ub wire nr ", 250, 0, 250);
362 hVa_wireOccupancy = new TH1F("hVa_wireOccupancy", " va wire nr ", 250, 0, 250);
363 hVb_wireOccupancy = new TH1F("hVb_wireOccupancy", " vb wire nr ", 250, 0, 250);
364 hXa2_wireOccupancy = new TH1F("hXa2_wireOccupancy", " xa2 wire nr ", 250, 0, 250);
365 hXb2_wireOccupancy = new TH1F("hXb2_wireOccupancy", " xb2 wire nr ", 250, 0, 250);
366 hYa2_wireOccupancy = new TH1F("hYa2_wireOccupancy", " ya2 wire nr ", 250, 0, 250);
367 hYb2_wireOccupancy = new TH1F("hYb2_wireOccupancy", " yb2 wire nr ", 250, 0, 250);
368 hUa2_wireOccupancy = new TH1F("hUa2_wireOccupancy", " ua2 wire nr ", 250, 0, 250);
369 hUb2_wireOccupancy = new TH1F("hUb2_wireOccupancy", " ub2 wire nr ", 250, 0, 250);
370 hVa2_wireOccupancy = new TH1F("hVa2_wireOccupancy", " va2 wire nr ", 250, 0, 250);
371 hVb2_wireOccupancy = new TH1F("hVb2_wireOccupancy", " vb2 wire nr ", 250, 0, 250);
372 haX2_aX1 = new TH1F("haX2_aX1", "aX2 - aX1", 200, -0.25, 0.25);
373 haY2_aY1 = new TH1F("haY2_aY1", "aY2 - aY1", 200, -0.25, 0.25);
374 hX1_X2_glob_all = new TH1F("hX1_X2_glob_all", " Xdc2 - Xdc1 extrapol Zmid, all closest", 300, -30, 30);
375 hY1_Y2_glob_all = new TH1F("hY1_Y2_glob_all", " Ydc2 - Ydc1 extrapol Zmid, all closest", 300, -30, 30);
376 haX2_aX1m = new TH1F("haX2_aX1m", "aX2 - aX1m", 200, -0.1, 0.1);
377 haY2_aY1m = new TH1F("haY2_aY1m", "aY2 - aY1m", 200, -0.1, 0.1);
378 hX1_X2_glob_m = new TH1F("hX1_X2_glob_m", " Xdc2 - Xdc1 extrapol Zmid, m closest", 300, -15, 15);
379 hY1_Y2_glob_m = new TH1F("hY1_Y2_glob_m", " Ydc2 - Ydc1 extrapol Zmid, m closest", 300, -15, 15);
380 hx1Da_Db = new TH1F("hx1Da_Db", "d_a + d_b - 5mm dc1", 200, -0.5, 0.5);
381 hx2Da_Db = new TH1F("hx2Da_Db", "d_a + d_b - 5mm dc2", 200, -0.5, 0.5);
382 hy1Da_Db = new TH1F("hy1Da_Db", "d_a + d_b - 5mm dc1", 200, -0.5, 0.5);
383 hy2Da_Db = new TH1F("hy2Da_Db", "d_a + d_b - 5mm dc2", 200, -0.5, 0.5);
384 hu1Da_Db = new TH1F("hu1Da_Db", "d_a + d_b - 5mm dc1", 200, -0.5, 0.5);
385 hu2Da_Db = new TH1F("hu2Da_Db", "d_a + d_b - 5mm dc2", 200, -0.5, 0.5);
386 hv1Da_Db = new TH1F("hv1Da_Db", "d_a + d_b - 5mm dc1", 200, -0.5, 0.5);
387 hv2Da_Db = new TH1F("hv2Da_Db", "d_a + d_b - 5mm dc2", 200, -0.5, 0.5);
388 hx1Da_Db_best = new TH1F("hx1Da_Db_best", "d_a + d_b - 5mm for best dc1", 200, -0.5, 0.5);
389 hx2Da_Db_best = new TH1F("hx2Da_Db_best", "d_a + d_b - 5mm for best dc2", 200, -0.5, 0.5);
390 hy1Da_Db_best = new TH1F("hy1Da_Db_best", "d_a + d_b - 5mm for best dc1", 200, -0.5, 0.5);
391 hy2Da_Db_best = new TH1F("hy2Da_Db_best", "d_a + d_b - 5mm for best dc2", 200, -0.5, 0.5);
392 hu1Da_Db_best = new TH1F("hu1Da_Db_best", "d_a + d_b - 5mm for best dc1", 200, -0.5, 0.5);
393 hu2Da_Db_best = new TH1F("hu2Da_Db_best", "d_a + d_b - 5mm for best dc2", 200, -0.5, 0.5);
394 hv1Da_Db_best = new TH1F("hv1Da_Db_best", "d_a + d_b - 5mm for best dc1", 200, -0.5, 0.5);
395 hv2Da_Db_best = new TH1F("hv2Da_Db_best", "d_a + d_b - 5mm for best dc2", 200, -0.5, 0.5);
396 haX1 = new TH1F("haX1", "aX1", 300, -0.05, 0.15);
397 haY1 = new TH1F("haY1", "aY1", 200, -0.1, 0.1);
398 haX2 = new TH1F("haX2", "aX2", 300, -0.05, 0.15);
399 haY2 = new TH1F("haY2", "aY2", 200, -0.1, 0.1);
400 hx1 = new TH1F("hx1", "X1", 130, -80, 50);
401 hy1 = new TH1F("hy1", "Y1", 200, -30, 50);
402 hx2 = new TH1F("hx2", "X2", 130, -80, 50);
403 hy2 = new TH1F("hy2", "Y2", 200, -30, 50);
404 hX = new TH1F("hX", "X", 80, -80, 0);
405 hY = new TH1F("hY", "Y", 90, -10, 20);
406 haX = new TH1F("haX", "aX", 300, -0.05, 0.15);
407 haY = new TH1F("haY", "aY", 200, -0.1, 0.1);
408 hlocXY1 = new TH2F("hlocXY1", " (X,Y) local coord of a seg in dc1", 240, -120, 120, 240, -120, 120);
409 hlocXY2 = new TH2F("hlocXY2", " (X,Y) local coord of a seg in dc2", 240, -120, 120, 240, -120, 120);
410 hXvsAx = new TH2F("hXvsAx", " X coord of a seg in dc1 vs x slope", 240, -120, 120, 100, -10, 10);
411
412 hAx1_loc = new TH1F("hAx1_loc", "slope x dc1 in loc coord system", 400, -0.2, 0.2);
413 hAy1_loc = new TH1F("hAy1_loc", "slope y dc1 in loc coord system", 400, -0.2, 0.2);
414 hAu1_loc = new TH1F("hAu1_loc", "slope u dc1 in loc coord system", 400, -0.2, 0.2);
415 hAv1_loc = new TH1F("hAv1_loc", "slope v dc1 in loc coord system", 400, -0.2, 0.2);
416
417 hAx12_loc = new TH1F("hAx12_loc", "slope x dc1 in loc coord system", 200, -0.2, 0.2);
418 hAy12_loc = new TH1F("hAy12_loc", "slope y dc1 in loc coord system", 200, -0.2, 0.2);
419 hAu12_loc = new TH1F("hAu12_loc", "slope u dc1 in loc coord system", 200, -0.2, 0.2);
420 hAv12_loc = new TH1F("hAv12_loc", "slope v dc1 in loc coord system", 200, -0.2, 0.2);
421
422 hAx2_loc = new TH1F("hAx2_loc", "slope x dc2 in loc coord system", 400, -0.2, 0.2);
423 hAy2_loc = new TH1F("hAy2_loc", "slope y dc2 in loc coord system", 400, -0.2, 0.2);
424 hAu2_loc = new TH1F("hAu2_loc", "slope u dc2 in loc coord system", 400, -0.2, 0.2);
425 hAv2_loc = new TH1F("hAv2_loc", "slope v dc2 in loc coord system", 400, -0.2, 0.2);
426
427 hAx122loc = new TH1F("hAx122loc", "slope x dc2 in loc coord system", 200, -0.2, 0.2);
428 hAy122loc = new TH1F("hAy122loc", "slope y dc2 in loc coord system", 200, -0.2, 0.2);
429 hAu122loc = new TH1F("hAu122loc", "slope u dc2 in loc coord system", 200, -0.2, 0.2);
430 hAv122loc = new TH1F("hAv122loc", "slope v dc2 in loc coord system", 200, -0.2, 0.2);
431
432 hXm_Xe_loc = new TH1F("hXm_Xe_loc", " x_mes - x_est(u,v)", 100, -1.5, 1.5);
433 hYm_Ye_loc = new TH1F("hYm_Ye_loc", " y_mes - y_est(u,v)", 100, -1.5, 1.5);
434 hUm_Ue_loc = new TH1F("hUm_Ue_loc", " u_mes - u_est(x,y)", 100, -1.5, 1.5);
435 hVm_Ve_loc = new TH1F("hVm_Ve_loc", " v_mes - v_est(x,y)", 100, -1.5, 1.5);
436
437 hDc1XaResVsDa = new TProfile("hDc1XaResVsDa", " res x1a vs da", 100, 0, .5, 0, 5);
438 hDc1XbResVsDb = new TProfile("hDc1XbResVsDb", " res x1b vs db", 100, 0, .5, 0, 5);
439 hDc1XaResVsDa_p = new TProfile("hDc1XaResVsDa_p", " res x1a vs da +", 100, 0, .5, 0, 5);
440 hDc1XbResVsDb_p = new TProfile("hDc1XbResVsDb_p", " res x1b vs db +", 100, 0, .5, 0, 5);
441 hDc1XaResVsDa_m = new TProfile("hDc1XaResVsDa_m", " res x1a vs da -", 100, 0, .5, 0, 5);
442 hDc1XbResVsDb_m = new TProfile("hDc1XbResVsDb_m", " res x1b vs db -", 100, 0, .5, 0, 5);
443
444 hUvsXV1 = new TH1F("hUvsXV1", "u1 diff", 100, -2, 2);
445 hYvsXV1 = new TH1F("hYvsXV1", "y1 diff", 100, -2, 2);
446 hUvsXV2 = new TH1F("hUvsXV2", "u2 diff", 100, -2, 2);
447 hYvsXV2 = new TH1F("hYvsXV2", "y2 diff", 100, -2, 2);
448 hChi2ndf_dc1 = new TH1F("hChi2ndf_dc1", "chi2/nDoF for dc1 segments", 300, 0, 300);
449 hChi2ndf_dc2 = new TH1F("hChi2ndf_dc2", "chi2/nDoF for dc2 segments", 300, 0, 300);
450 hSegSize_dc1 = new TH1F("hSegSize_dc1", " segment size for dc1", 9, 0, 9);
451 hSegSize_dc2 = new TH1F("hSegSize_dc2", " segment size for dc2", 9, 0, 9);
452 hChi2Matched = new TH1F("hChi2Matched", " chi2 for matched pair", 2000, 0, 2000);
453 hP = new TH1F("hP", " momentum = .3*Int(BL)/[sin(alphaX_out-alphaX_in)+C]", 200, 5, 15);
454
455 fhList.Add(hNLayWithFiredWiresDc1);
456 fhList.Add(hNLayWithFiredWiresDc2);
457 fhList.Add(hNoFiredWiresOnLayerDc1);
458 fhList.Add(hNoFiredWiresOnLayerDc2);
459
460 fhList.Add(hNomin1);
461 fhList.Add(hDenom1);
462 fhList.Add(hNomin2);
463 fhList.Add(hDenom2);
464
465 fhList.Add(hEff1);
466 fhList.Add(hEff2);
467 fhList.Add(hSlot_1xa_time8p);
468 fhList.Add(hSlot_1xb_time8p);
469 fhList.Add(hSlot_1ya_time8p);
470 fhList.Add(hSlot_1yb_time8p);
471 fhList.Add(hSlot_1ua_time8p);
472 fhList.Add(hSlot_1ub_time8p);
473 fhList.Add(hSlot_1va_time8p);
474 fhList.Add(hSlot_1vb_time8p);
475 fhList.Add(hSlot_2xa_time8p);
476 fhList.Add(hSlot_2xb_time8p);
477 fhList.Add(hSlot_2ya_time8p);
478 fhList.Add(hSlot_2yb_time8p);
479 fhList.Add(hSlot_2ua_time8p);
480 fhList.Add(hSlot_2ub_time8p);
481 fhList.Add(hSlot_2va_time8p);
482 fhList.Add(hSlot_2vb_time8p);
483
484 fhList.Add(hUvsXV1);
485 fhList.Add(hYvsXV1);
486 fhList.Add(hUvsXV2);
487 fhList.Add(hYvsXV2);
488 fhList.Add(hDc1XaResVsDa);
489 fhList.Add(hDc1XbResVsDb);
490
491 fhList.Add(hDc1XaResVsDa_m);
492 fhList.Add(hDc1XbResVsDb_m);
493 fhList.Add(hDc1XaResVsDa_p);
494 fhList.Add(hDc1XbResVsDb_p);
495
496 fhList.Add(hXm_Xe_loc);
497 fhList.Add(hYm_Ye_loc);
498 fhList.Add(hUm_Ue_loc);
499 fhList.Add(hVm_Ve_loc);
500
501 fhList.Add(hAx1_loc);
502 fhList.Add(hAy1_loc);
503 fhList.Add(hAu1_loc);
504 fhList.Add(hAv1_loc);
505 fhList.Add(hAx12_loc);
506 fhList.Add(hAy12_loc);
507 fhList.Add(hAu12_loc);
508 fhList.Add(hAv12_loc);
509 fhList.Add(hAx2_loc);
510 fhList.Add(hAy2_loc);
511 fhList.Add(hAu2_loc);
512 fhList.Add(hAv2_loc);
513 fhList.Add(hAx122loc);
514 fhList.Add(hAy122loc);
515 fhList.Add(hAu122loc);
516 fhList.Add(hAv122loc);
517 fhList.Add(hResidx1a);
518 fhList.Add(hResidy1a);
519 fhList.Add(hResidu1a);
520 fhList.Add(hResidv1a);
521 fhList.Add(hResidx2a);
522 fhList.Add(hResidy2a);
523 fhList.Add(hResidu2a);
524 fhList.Add(hResidv2a);
525
526 fhList.Add(hResidx1b);
527 fhList.Add(hResidy1b);
528 fhList.Add(hResidu1b);
529 fhList.Add(hResidv1b);
530 fhList.Add(hResidx2b);
531 fhList.Add(hResidy2b);
532 fhList.Add(hResidu2b);
533 fhList.Add(hResidv2b);
534
535 fhList.Add(hResidx1a_0p1);
536 fhList.Add(hResidy1a_0p1);
537 fhList.Add(hResidu1a_0p1);
538 fhList.Add(hResidv1a_0p1);
539 fhList.Add(hResidx2a_0p1);
540 fhList.Add(hResidy2a_0p1);
541 fhList.Add(hResidu2a_0p1);
542 fhList.Add(hResidv2a_0p1);
543
544 fhList.Add(hResidx1a_0p1_0p4);
545 fhList.Add(hResidy1a_0p1_0p4);
546 fhList.Add(hResidu1a_0p1_0p4);
547 fhList.Add(hResidv1a_0p1_0p4);
548 fhList.Add(hResidx2a_0p1_0p4);
549 fhList.Add(hResidy2a_0p1_0p4);
550 fhList.Add(hResidu2a_0p1_0p4);
551 fhList.Add(hResidv2a_0p1_0p4);
552
553 fhList.Add(hResidx1a_0p4_0p5);
554 fhList.Add(hResidy1a_0p4_0p5);
555 fhList.Add(hResidu1a_0p4_0p5);
556 fhList.Add(hResidv1a_0p4_0p5);
557 fhList.Add(hResidx2a_0p4_0p5);
558 fhList.Add(hResidy2a_0p4_0p5);
559 fhList.Add(hResidu2a_0p4_0p5);
560 fhList.Add(hResidv2a_0p4_0p5);
561
562 fhList.Add(hResidx1b_0p1);
563 fhList.Add(hResidy1b_0p1);
564 fhList.Add(hResidu1b_0p1);
565 fhList.Add(hResidv1b_0p1);
566 fhList.Add(hResidx2b_0p1);
567 fhList.Add(hResidy2b_0p1);
568 fhList.Add(hResidu2b_0p1);
569 fhList.Add(hResidv2b_0p1);
570
571 fhList.Add(hResidx1b_0p1_0p4);
572 fhList.Add(hResidy1b_0p1_0p4);
573 fhList.Add(hResidu1b_0p1_0p4);
574 fhList.Add(hResidv1b_0p1_0p4);
575 fhList.Add(hResidx2b_0p1_0p4);
576 fhList.Add(hResidy2b_0p1_0p4);
577 fhList.Add(hResidu2b_0p1_0p4);
578 fhList.Add(hResidv2b_0p1_0p4);
579
580 fhList.Add(hResidx1b_0p4_0p5);
581 fhList.Add(hResidy1b_0p4_0p5);
582 fhList.Add(hResidu1b_0p4_0p5);
583 fhList.Add(hResidv1b_0p4_0p5);
584 fhList.Add(hResidx2b_0p4_0p5);
585 fhList.Add(hResidy2b_0p4_0p5);
586 fhList.Add(hResidu2b_0p4_0p5);
587 fhList.Add(hResidv2b_0p4_0p5);
588
589 // puls
590 fhList.Add(hP);
591 fhList.Add(hPullsx1a_0p1);
592 fhList.Add(hPullsy1a_0p1);
593 fhList.Add(hPullsu1a_0p1);
594 fhList.Add(hPullsv1a_0p1);
595 fhList.Add(hPullsx2a_0p1);
596 fhList.Add(hPullsy2a_0p1);
597 fhList.Add(hPullsu2a_0p1);
598 fhList.Add(hPullsv2a_0p1);
599
600 fhList.Add(hPullsx1a_0p1_0p4);
601 fhList.Add(hPullsy1a_0p1_0p4);
602 fhList.Add(hPullsu1a_0p1_0p4);
603 fhList.Add(hPullsv1a_0p1_0p4);
604 fhList.Add(hPullsx2a_0p1_0p4);
605 fhList.Add(hPullsy2a_0p1_0p4);
606 fhList.Add(hPullsu2a_0p1_0p4);
607 fhList.Add(hPullsv2a_0p1_0p4);
608
609 fhList.Add(hPullsx1a_0p4_0p5);
610 fhList.Add(hPullsy1a_0p4_0p5);
611 fhList.Add(hPullsu1a_0p4_0p5);
612 fhList.Add(hPullsv1a_0p4_0p5);
613 fhList.Add(hPullsx2a_0p4_0p5);
614 fhList.Add(hPullsy2a_0p4_0p5);
615 fhList.Add(hPullsu2a_0p4_0p5);
616 fhList.Add(hPullsv2a_0p4_0p5);
617
618 fhList.Add(hPullsx1b_0p1);
619 fEventNo = 0;
620
621 // cout<<" x2 sh "<<x2_sh<<endl;
622 fhList.Add(hPullsy1b_0p1);
623 fhList.Add(hPullsu1b_0p1);
624 fhList.Add(hPullsv1b_0p1);
625 fhList.Add(hPullsx2b_0p1);
626 fhList.Add(hPullsy2b_0p1);
627 fhList.Add(hPullsu2b_0p1);
628 fhList.Add(hPullsv2b_0p1);
629
630 fhList.Add(hPullsx1b_0p1_0p4);
631 fhList.Add(hPullsy1b_0p1_0p4);
632 fhList.Add(hPullsu1b_0p1_0p4);
633 fhList.Add(hPullsv1b_0p1_0p4);
634 fhList.Add(hPullsx2b_0p1_0p4);
635 fhList.Add(hPullsy2b_0p1_0p4);
636 fhList.Add(hPullsu2b_0p1_0p4);
637 fhList.Add(hPullsv2b_0p1_0p4);
638
639 fhList.Add(hPullsx1b_0p4_0p5);
640 fhList.Add(hPullsy1b_0p4_0p5);
641 fhList.Add(hPullsu1b_0p4_0p5);
642 fhList.Add(hPullsv1b_0p4_0p5);
643 fhList.Add(hPullsx2b_0p4_0p5);
644 fhList.Add(hPullsy2b_0p4_0p5);
645 fhList.Add(hPullsu2b_0p4_0p5);
646 fhList.Add(hPullsv2b_0p4_0p5);
647
648 fhList.Add(hXDC1_atZ0);
649 fhList.Add(hYDC1_atZ0);
650 fhList.Add(hXDC2_atZ0);
651 fhList.Add(hYDC2_atZ0);
652 fhList.Add(hSlot_1xa_time);
653 fhList.Add(hSlot_1xb_time);
654 fhList.Add(haX2_aX1);
655 fhList.Add(haY2_aY1);
656 fhList.Add(hX1_X2_glob_all);
657 fhList.Add(hY1_Y2_glob_all);
658 fhList.Add(haX2_aX1m);
659 fhList.Add(haY2_aY1m);
660 fhList.Add(hX1_X2_glob_m);
661 fhList.Add(hY1_Y2_glob_m);
662
663 fhList.Add(hXa_wireOccupancy);
664 fhList.Add(hXb_wireOccupancy);
665 fhList.Add(hYa_wireOccupancy);
666 fhList.Add(hYb_wireOccupancy);
667 fhList.Add(hUa_wireOccupancy);
668 fhList.Add(hUb_wireOccupancy);
669 fhList.Add(hVa_wireOccupancy);
670 fhList.Add(hVb_wireOccupancy);
671 fhList.Add(hXa2_wireOccupancy);
672 fhList.Add(hXb2_wireOccupancy);
673 fhList.Add(hYa2_wireOccupancy);
674 fhList.Add(hYb2_wireOccupancy);
675 fhList.Add(hUa2_wireOccupancy);
676 fhList.Add(hUb2_wireOccupancy);
677 fhList.Add(hVa2_wireOccupancy);
678 fhList.Add(hVb2_wireOccupancy);
679
680 fhList.Add(hChi2ndf_dc1);
681 fhList.Add(hChi2ndf_dc2);
682 fhList.Add(hChi2Matched);
683 fhList.Add(hSegSize_dc1);
684 fhList.Add(hSegSize_dc2);
685 fhList.Add(Ntrack1);
686 fhList.Add(Ntrack2);
687 fhList.Add(NGtracks);
688 fhList.Add(hx1Da_Db);
689 fhList.Add(hy1Da_Db);
690 fhList.Add(hu1Da_Db);
691 fhList.Add(hv1Da_Db);
692
693 fhList.Add(hx2Da_Db);
694 fhList.Add(hy2Da_Db);
695 fhList.Add(hu2Da_Db);
696 fhList.Add(hv2Da_Db);
697 fhList.Add(hx1Da_Db_best);
698 fhList.Add(hy1Da_Db_best);
699 fhList.Add(hu1Da_Db_best);
700 fhList.Add(hv1Da_Db_best);
701
702 fhList.Add(hx2Da_Db_best);
703 fhList.Add(hy2Da_Db_best);
704 fhList.Add(hu2Da_Db_best);
705 fhList.Add(hv2Da_Db_best);
706 fhList.Add(haX1);
707 fhList.Add(haY1);
708 fhList.Add(haX2);
709 fhList.Add(haY2);
710 fhList.Add(hx1);
711 fhList.Add(hy1);
712 fhList.Add(hx2);
713 fhList.Add(hy2);
714 fhList.Add(hX);
715 fhList.Add(hY);
716 fhList.Add(haX);
717 fhList.Add(haY);
718 fhList.Add(hlocXY1);
719 fhList.Add(hlocXY2);
720 fhList.Add(hXvsAx);
721 // end
722}
723
724BmnDchTrackFinder::~BmnDchTrackFinder()
725{
726 if (isInitialized) {
727 // Delete 1d-arrays
728 delete[] nSegments;
729 delete[] has7DC;
730 delete[] x_mid;
731 delete[] y_mid;
732 delete[] nhits;
733 delete[] pairID;
734 delete[] aX;
735 delete[] aY;
736 delete[] leng;
737 delete[] imp;
738 // Delete 2d-arrays and 3d-arrays
739 for (Int_t iChamber = 0; iChamber < nChambers; iChamber++) {
740 delete[] x_global[iChamber];
741 delete[] y_global[iChamber];
742 delete[] Chi2[iChamber];
743 delete[] SegMCId[iChamber];
744 delete[] SegMCIdCount[iChamber];
745 delete[] pairs[iChamber];
746 delete[] segment_size[iChamber];
747 for (Int_t iDim = 0; iDim < N; iDim++) {
748 delete[] v[iChamber][iDim];
749 delete[] u[iChamber][iDim];
750 delete[] y[iChamber][iDim];
751 delete[] x[iChamber][iDim];
752 delete[] vId[iChamber][iDim];
753 delete[] uId[iChamber][iDim];
754 delete[] yId[iChamber][iDim];
755 delete[] xId[iChamber][iDim];
756 delete[] v_time[iChamber][iDim];
757 delete[] u_time[iChamber][iDim];
758 delete[] y_time[iChamber][iDim];
759 delete[] x_time[iChamber][iDim];
760 delete[] sigm_v[iChamber][iDim];
761 delete[] sigm_u[iChamber][iDim];
762 delete[] sigm_y[iChamber][iDim];
763 delete[] sigm_x[iChamber][iDim];
764 delete[] v_Single[iChamber][iDim];
765 delete[] u_Single[iChamber][iDim];
766 delete[] y_Single[iChamber][iDim];
767 delete[] x_Single[iChamber][iDim];
768 delete[] v_SingleId[iChamber][iDim];
769 delete[] u_SingleId[iChamber][iDim];
770 delete[] y_SingleId[iChamber][iDim];
771 delete[] x_SingleId[iChamber][iDim];
772 delete[] Sigm_v_single[iChamber][iDim];
773 delete[] Sigm_u_single[iChamber][iDim];
774 delete[] Sigm_y_single[iChamber][iDim];
775 delete[] Sigm_x_single[iChamber][iDim];
776 }
777 for (Int_t iDim = 0; iDim < 4; iDim++) {
778 delete[] params[iChamber][iDim];
779 delete[] params_sigmas[iChamber][iDim];
780 }
781 for (Int_t iDim = 0; iDim < 8; iDim++) {
782 delete[] rh_segment[iChamber][iDim];
783 delete[] rhId_segment[iChamber][iDim];
784 delete[] rh_segment_time[iChamber][iDim];
785 delete[] rh_sigm_segment[iChamber][iDim];
786 }
787 delete[] v[iChamber];
788 delete[] u[iChamber];
789 delete[] y[iChamber];
790 delete[] x[iChamber];
791 delete[] vId[iChamber];
792 delete[] uId[iChamber];
793 delete[] yId[iChamber];
794 delete[] xId[iChamber];
795 delete[] v_time[iChamber];
796 delete[] u_time[iChamber];
797 delete[] y_time[iChamber];
798 delete[] x_time[iChamber];
799 delete[] sigm_v[iChamber];
800 delete[] sigm_u[iChamber];
801 delete[] sigm_y[iChamber];
802 delete[] sigm_x[iChamber];
803 delete[] v_Single[iChamber];
804 delete[] u_Single[iChamber];
805 delete[] y_Single[iChamber];
806 delete[] x_Single[iChamber];
807 delete[] v_SingleId[iChamber];
808 delete[] u_SingleId[iChamber];
809 delete[] y_SingleId[iChamber];
810 delete[] x_SingleId[iChamber];
811 delete[] Sigm_v_single[iChamber];
812 delete[] Sigm_u_single[iChamber];
813 delete[] Sigm_y_single[iChamber];
814 delete[] Sigm_x_single[iChamber];
815 delete[] params[iChamber];
816 delete[] params_sigmas[iChamber];
817 delete[] rh_segment[iChamber];
818 delete[] rhId_segment[iChamber];
819 delete[] rh_segment_time[iChamber];
820 delete[] rh_sigm_segment[iChamber];
821 for (Int_t iWire = 0; iWire < nWires; iWire++)
822 delete[] singles[iChamber][iWire];
823 delete[] singles[iChamber];
824 }
825 }
826
827 delete[] x_global;
828 delete[] y_global;
829 delete[] Chi2;
830 delete[] pairs;
831 delete[] SegMCId;
832 delete[] SegMCIdCount;
833 delete[] segment_size;
834 delete[] v;
835 delete[] u;
836 delete[] y;
837 delete[] x;
838 delete[] v_time;
839 delete[] u_time;
840 delete[] y_time;
841 delete[] x_time;
842 delete[] sigm_v;
843 delete[] sigm_u;
844 delete[] sigm_y;
845 delete[] sigm_x;
846 delete[] v_Single;
847 delete[] u_Single;
848 delete[] y_Single;
849 delete[] x_Single;
850 delete[] v_SingleId;
851 delete[] u_SingleId;
852 delete[] y_SingleId;
853 delete[] x_SingleId;
854 delete[] Sigm_v_single;
855 delete[] Sigm_u_single;
856 delete[] Sigm_y_single;
857 delete[] Sigm_x_single;
858 delete[] params;
859 delete[] params_sigmas;
860 delete[] rh_segment;
861 delete[] rhId_segment;
862 delete[] rh_segment_time;
863 delete[] rh_sigm_segment;
864 delete[] singles;
865}
866
867void BmnDchTrackFinder::Exec(Option_t* opt)
868{
869 if (!IsActive())
870 return;
871 fEventNo++;
872 clock_t tStart = clock();
873 PrepareArraysToProcessEvent();
874 if (fVerbose)
875 cout << "======================== DCH track finder exec started ====================" << endl;
876 if (fVerbose)
877 cout << "Event number: " << fEventNo << endl;
878 // temporary containers
879 // Order used: va1, vb1, ua1, ub1, ya1, yb1, xa1, xb1 (dch1, 0 - 7) - va2, vb2, ua2, ub2, ya2, yb2, xa2, xb2 (dch2,
880 // 8 - 15)
881 const Int_t nDim = 80; // max number of fired wires per plane
882 const Int_t nPlanes = 16; // Total number of planes in both DCHs (0-7, 8-15)
883 // cout << "------------------------------------------------------------Event number: " << fEventNo << endl;
884 Double_t driftLength[nPlanes][nDim];
885 Double_t wires[nPlanes][nDim];
886 Int_t hitId[nPlanes][nDim];
887 // Double_t xMC[nPlanes][nDim];
888 // Double_t yMC[nPlanes][nDim];
889 // Double_t zMC[nPlanes][nDim];
890
891 Int_t it[nPlanes];
892 Bool_t used[nPlanes][nDim];
893 for (Int_t iPlanes = 0; iPlanes < nPlanes; iPlanes++) {
894 it[iPlanes] = 0;
895 for (Int_t iDim = 0; iDim < nDim; iDim++) {
896 used[iPlanes][iDim] = kFALSE;
897 driftLength[iPlanes][iDim] = 0.;
898 hitId[iPlanes][iDim] = -999;
899 wires[iPlanes][iDim] = 0.;
900 // xMC[iPlanes][iDim] = 999.;
901 // yMC[iPlanes][iDim] = 999.;
902 // zMC[iPlanes][iDim] = 999.;
903 }
904 }
905
906 if (expData) {
907 Bool_t goodEv = kTRUE;
908 Bool_t written = kFALSE;
909
910 for (Int_t iDig = 0; iDig < fBmnDchDigitsArray->GetEntriesFast(); ++iDig) {
911 BmnDchDigit* digit = (BmnDchDigit*)fBmnDchDigitsArray->UncheckedAt(iDig);
912 // skip identical events
913 if (!written) {
914 written = kTRUE;
915 if (digit->GetTime() == prev_time && digit->GetWireNumber() == prev_wire) {
916 goodEv = kFALSE;
917 } else {
918 prev_time = Int_t(digit->GetTime());
919 prev_wire = Int_t(digit->GetWireNumber());
920 }
921 } //! written
922
923 if (!goodEv)
924 return;
925
926 // Order used: vb1(0), va1(1), ub1(2), ua1(3), yb1(4), ya1(5), xb1(6), xa1(7) ->
927 // vb2(8), va2(9), ub2(10), ua2(11), yb2(12), ya2(13), xb2(14), xa2(15)
928 Short_t plane = digit->GetPlane();
929 Short_t wire = digit->GetWireNumber();
930 Double_t time = digit->GetTime();
931 Bool_t secondaries = kFALSE;
932 if (wire > 239)
933 wire = wire - 128; // recalculate last amplifier
934 if (plane == 4) {
935 if (wire > 143 && wire < 176) {
936 if (wire > 159 && wire < 176) {
937 wire = wire - 16;
938 } else {
939 wire = wire + 16;
940 }
941 }
942 }
943
944 for (Int_t sec = 0; sec < it[plane] - 1; sec++)
945 if (wire == wires[plane][sec]) {
946 secondaries = kTRUE;
947 break;
948 }
949
950 if (it[plane] == (nDim - 1) || secondaries)
951 continue;
952
953 wires[plane][it[plane]] = wire;
954 driftLength[plane][it[plane]] = time;
955
956 it[plane]++;
957 if (plane == 6)
958 hSlot_1xa_time->Fill(time);
959 if (plane == 7)
960 hSlot_1xb_time->Fill(time);
961
962 if (plane == 1)
963 hVa_wireOccupancy->Fill(wire);
964 if (plane == 0)
965 hVb_wireOccupancy->Fill(wire);
966 if (plane == 3)
967 hUa_wireOccupancy->Fill(wire);
968 if (plane == 2)
969 hUb_wireOccupancy->Fill(wire);
970 if (plane == 5)
971 hYa_wireOccupancy->Fill(wire);
972 if (plane == 4)
973 hYb_wireOccupancy->Fill(wire);
974 if (plane == 7)
975 hXa_wireOccupancy->Fill(wire);
976 if (plane == 6)
977 hXb_wireOccupancy->Fill(wire);
978 if (plane == 9)
979 hVa2_wireOccupancy->Fill(wire);
980 if (plane == 8)
981 hVb2_wireOccupancy->Fill(wire);
982 if (plane == 11)
983 hUa2_wireOccupancy->Fill(wire);
984 if (plane == 10)
985 hUb2_wireOccupancy->Fill(wire);
986 if (plane == 13)
987 hYa2_wireOccupancy->Fill(wire);
988 if (plane == 12)
989 hYb2_wireOccupancy->Fill(wire);
990 if (plane == 15)
991 hXa2_wireOccupancy->Fill(wire);
992 if (plane == 14)
993 hXb2_wireOccupancy->Fill(wire);
994
995 } // iDig
996 } // isExp
997 else
998 {
999
1000 // read MC data
1001 // Double_t x1_MC1 = 999;
1002 // Double_t y1_MC1 = 999;
1003 // Double_t z1_MC1 = 999;
1004 // Double_t x1_MC2 = 999;
1005 // Double_t y1_MC2 = 999;
1006 // Double_t z1_MC2 = 999;
1007
1008 // bool _MC1 = false;
1009 // bool _MC2 = false;
1010 for (Int_t iMC = 0; iMC < fBmnHitsArray->GetEntriesFast(); ++iMC) {
1011 BmnDchHit* hit = (BmnDchHit*)fBmnHitsArray->UncheckedAt(iMC);
1012
1013 // Order used: vb1(0), va1(1), ub1(2), ua1(3), yb1(4), ya1(5), xb1(6), xa1(7) ->
1014 // vb2(8), va2(9), ub2(10), ua2(11), yb2(12), ya2(13), xb2(14), xa2(15)
1015 Short_t plane = hit->GetLayerNumber();
1016 Short_t wire = hit->GetWirePosition();
1017 Double_t time = hit->GetDrift();
1018 // Double_t x_MC = hit->GetX();
1019 // Double_t y_MC = hit->GetY();
1020 // Double_t z_MC = hit->GetZ();
1021
1022 hitId[plane][it[plane]] = hit->GetHitId();
1023 wires[plane][it[plane]] = wire;
1024 driftLength[plane][it[plane]] = time;
1025 // xMC[plane][it[plane]] = x_MC;
1026 // yMC[plane][it[plane]] = y_MC;
1027 // zMC[plane][it[plane]] = z_MC;
1028 it[plane]++;
1029 if (plane == 7)
1030 hSlot_1xa_time->Fill(time);
1031 if (plane == 6)
1032 hSlot_1xb_time->Fill(time);
1033
1034 } // end MC read
1035 }
1036
1037 const Int_t nWires_ = 4;
1038 const Int_t nLayers_ = 2;
1039 TString wireNames[nWires_] = {"v", "u", "y", "x"};
1040 TString layNames[nLayers_] = {"a", "b"};
1041 Int_t cntr = 0;
1042
1043 for (Int_t iDch = 0; iDch < nChambers; iDch++) {
1044 for (Int_t iWire = 0; iWire < nWires; iWire++) {
1045 Int_t start = 2 * iWire + (nPlanes / 2) * iDch;
1046 Int_t finish = start + 1;
1047 Double_t*** coord = (iWire == 0) ? v : (iWire == 1) ? u : (iWire == 2) ? y : x;
1048 Int_t*** coordId = (iWire == 0) ? vId : (iWire == 1) ? uId : (iWire == 2) ? yId : xId;
1049 Double_t*** sigma = (iWire == 0) ? sigm_v : (iWire == 1) ? sigm_u : (iWire == 2) ? sigm_y : sigm_x;
1050 Double_t*** coord_time = (iWire == 0) ? v_time : (iWire == 1) ? u_time : (iWire == 2) ? y_time : x_time;
1051 if (expData) {
1052 pairs[iDch][iWire] =
1053 Reconstruction(iDch + 1, wireNames[iWire], pairs[iDch][iWire], it[start], it[finish], wires[start],
1054 wires[finish], driftLength[start], driftLength[finish], used[start], used[finish],
1055 coord[iDch], sigma[iDch], coord_time[iDch]);
1056
1057 for (Int_t iLayer = 0; iLayer < nLayers; iLayer++) {
1058 Double_t*** single_coord = (iWire == 0) ? v_Single
1059 : (iWire == 1) ? u_Single
1060 : (iWire == 2) ? y_Single
1061 : x_Single;
1062 Double_t*** single_sigma = (iWire == 0) ? Sigm_v_single
1063 : (iWire == 1) ? Sigm_u_single
1064 : (iWire == 2) ? Sigm_y_single
1065 : Sigm_x_single;
1066 singles[iDch][iWire][iLayer] = ReconstructionSingle(
1067 iDch + 1, wireNames[iWire], layNames[iLayer], singles[iDch][iWire][iLayer], it[cntr],
1068 wires[cntr], driftLength[cntr], used[cntr], single_coord[iDch], single_sigma[iDch]);
1069 cntr++;
1070 }
1071 } else {
1072 pairs[iDch][iWire] = ReconstructionMC(
1073 iDch + 1, wireNames[iWire], pairs[iDch][iWire], it[start], it[finish], wires[start], wires[finish],
1074 driftLength[start], driftLength[finish], hitId[start], hitId[finish], used[start], used[finish],
1075 coord[iDch], coordId[iDch], sigma[iDch], coord_time[iDch]);
1076
1077 for (Int_t iLayer = 0; iLayer < nLayers; iLayer++) {
1078 Double_t*** single_coord = (iWire == 0) ? v_Single
1079 : (iWire == 1) ? u_Single
1080 : (iWire == 2) ? y_Single
1081 : x_Single;
1082 Int_t*** single_coordId = (iWire == 0) ? v_SingleId
1083 : (iWire == 1) ? u_SingleId
1084 : (iWire == 2) ? y_SingleId
1085 : x_SingleId;
1086 Double_t*** single_sigma = (iWire == 0) ? Sigm_v_single
1087 : (iWire == 1) ? Sigm_u_single
1088 : (iWire == 2) ? Sigm_y_single
1089 : Sigm_x_single;
1090 singles[iDch][iWire][iLayer] = ReconstructionSingleMC(
1091 iDch + 1, wireNames[iWire], layNames[iLayer], singles[iDch][iWire][iLayer], it[cntr],
1092 wires[cntr], driftLength[cntr], hitId[cntr], used[cntr], single_coord[iDch],
1093 single_coordId[iDch], single_sigma[iDch]);
1094 cntr++;
1095 }
1096 }
1097 }
1098
1099 nSegments[iDch] = BuildUVSegments(
1100 iDch + 1, pairs[iDch][1], pairs[iDch][0], pairs[iDch][3], pairs[iDch][2], singles[iDch][1][0],
1101 singles[iDch][1][1], singles[iDch][0][0], singles[iDch][0][1], x[iDch], y[iDch], u[iDch], v[iDch],
1102 xId[iDch], yId[iDch], uId[iDch], vId[iDch], x_time[iDch], y_time[iDch], u_time[iDch], v_time[iDch],
1103 sigm_x[iDch], sigm_y[iDch], sigm_u[iDch], sigm_v[iDch], rh_segment[iDch], rhId_segment[iDch],
1104 rh_segment_time[iDch], rh_sigm_segment[iDch], u_Single[iDch], v_Single[iDch], u_SingleId[iDch],
1105 v_SingleId[iDch], Sigm_u_single[iDch], Sigm_v_single[iDch]);
1106 nSegments[iDch] = BuildXYSegments(
1107 iDch + 1, pairs[iDch][1], pairs[iDch][0], pairs[iDch][3], pairs[iDch][2], singles[iDch][3][0],
1108 singles[iDch][3][1], singles[iDch][2][0], singles[iDch][2][1], x[iDch], y[iDch], u[iDch], v[iDch],
1109 xId[iDch], yId[iDch], uId[iDch], vId[iDch], x_time[iDch], y_time[iDch], u_time[iDch], v_time[iDch],
1110 sigm_x[iDch], sigm_y[iDch], sigm_u[iDch], sigm_v[iDch], rh_segment[iDch], rhId_segment[iDch],
1111 rh_segment_time[iDch], rh_sigm_segment[iDch], x_Single[iDch], y_Single[iDch], x_SingleId[iDch],
1112 y_SingleId[iDch], Sigm_x_single[iDch], Sigm_y_single[iDch]);
1113
1114 has7DC[iDch] = FitDchSegments(iDch + 1, segment_size[iDch], rh_segment[iDch], rh_sigm_segment[iDch],
1115 params[iDch], Chi2[iDch], x_global[iDch], y_global[iDch]); // Fit found segments
1116 SelectLongestAndBestSegments(iDch + 1, segment_size[iDch], rh_segment[iDch],
1117 Chi2[iDch]); // Save only longest and best chi2 segments
1118 FillSegmentParametersSigmas(iDch + 1, rh_segment[iDch], rh_sigm_segment[iDch], Chi2[iDch], params_sigmas[iDch]);
1119 FindSegmentTrackMCId(iDch + 1, rhId_segment[iDch], Chi2[iDch], SegMCId[iDch], SegMCIdCount[iDch]);
1120
1121 // count wire multiplicity per layer
1122 if (iDch == 0) {
1123 layers_with = 0;
1124 for (Int_t i = 0; i < 8; i++) {
1125 if (it[i] > 0)
1126 layers_with++;
1127 if (it[i] == 0)
1128 hNoFiredWiresOnLayerDc1->Fill(i);
1129 }
1130 hNLayWithFiredWiresDc1->Fill(layers_with);
1131 } else {
1132 layers_with2 = 0;
1133 for (Int_t i = 8; i < 16; i++) {
1134 if (it[i] > 0)
1135 layers_with2++;
1136 if (it[i] == 0)
1137 hNoFiredWiresOnLayerDc2->Fill(i - 8);
1138 }
1139 hNLayWithFiredWiresDc2->Fill(layers_with2);
1140 }
1141 // calculate segment parameters errors
1142 CreateDchTrack(iDch + 1, Chi2[iDch], params[iDch], segment_size[iDch], rh_segment[iDch], rh_segment_time[iDch],
1143 params_sigmas[iDch], SegMCId[iDch], SegMCIdCount[iDch]);
1144 // ! Fill segment parameters
1145 // - Next function works if the local shifts are set to 0.
1146 // - Be careful - the refit brackes segment parameters.
1147 // - This function must be used only for residual calculation. Everything else becomes wrong!
1148
1149 // FillPlaneResiduals(iDch + 1, segment_size[iDch], rh_segment[iDch], rh_sigm_segment[iDch], params[iDch],
1150 // Chi2[iDch]);
1151 }
1152
1153 // Try to match the reconstructed segments from the two chambers
1154 SegmentsToBeMatched();
1155 CreateDchTrack();
1156
1157 if (fVerbose)
1158 cout << "======================== DCH track finder exec finished ===================" << endl;
1159 clock_t tFinish = clock();
1160 workTime += ((Double_t)(tFinish - tStart)) / CLOCKS_PER_SEC;
1161}
1162
1163void BmnDchTrackFinder::SegmentsToBeMatched()
1164{
1165 while (true) {
1166 // Int_t match_dc1_seg = -1;
1167 // Double_t ax(0.), ay(0.), xMean(0.), yMean(0.);
1168
1169 int best1 = -1;
1170 int best2 = -1;
1171 double chi2_Match_min = 1200.;
1172
1173 Double_t dx = -999;
1174 Double_t dy = -999;
1175 Double_t daX = -999;
1176 Double_t daY = -999;
1177
1178 for (Int_t segdc2Nr = 0; segdc2Nr < nSegments[1]; segdc2Nr++) {
1179 if (Chi2[1][segdc2Nr] > 998.)
1180 continue;
1181
1182 // Double_t chi2_match = 0;
1183 for (int segdc1Nr = 0; segdc1Nr < nSegments[0]; segdc1Nr++) {
1184 if (Chi2[0][segdc1Nr] > 998.)
1185 continue;
1186
1187 dx = x_global[1][segdc2Nr] - x_global[0][segdc1Nr];
1188 dy = y_global[1][segdc2Nr] - y_global[0][segdc1Nr];
1189 // xMean = 0.5 * (x_global[0][segdc1Nr] + x_global[1][segdc2Nr]);
1190 // yMean = 0.5 * (y_global[0][segdc1Nr] + y_global[1][segdc2Nr]);
1191 // ax = (params[1][1][segdc2Nr] - params[0][1][segdc1Nr]) / dZ_dch;
1192 // ay = (params[1][3][segdc2Nr] - params[0][3][segdc1Nr]) / dZ_dch;
1193 // match_dc1_seg = segdc1Nr;
1194 daX = params[1][0][segdc2Nr] - params[0][0][segdc1Nr];
1195 daY = params[1][2][segdc2Nr] - params[0][2][segdc1Nr];
1196 hX1_X2_glob_all->Fill(dx);
1197 hY1_Y2_glob_all->Fill(dy);
1198 haX2_aX1->Fill(daX);
1199 haY2_aY1->Fill(daY);
1200
1201 double sigma2_dx = ((params_sigmas[1][0][segdc2Nr] * params_sigmas[1][0][segdc2Nr])
1202 + (params_sigmas[0][0][segdc1Nr] * params_sigmas[0][0][segdc1Nr]))
1203 * 99.75 * 99.75
1204 + ((params_sigmas[0][1][segdc1Nr]) * (params_sigmas[0][1][segdc1Nr]))
1205 + ((params_sigmas[1][1][segdc2Nr]) * (params_sigmas[1][1][segdc2Nr]));
1206
1207 double sigma2_dy = ((params_sigmas[1][2][segdc2Nr] * params_sigmas[1][2][segdc2Nr])
1208 + (params_sigmas[0][2][segdc1Nr] * params_sigmas[0][2][segdc1Nr]))
1209 * 99.75 * 99.75
1210 + ((params_sigmas[0][3][segdc1Nr]) * (params_sigmas[0][3][segdc1Nr]))
1211 + ((params_sigmas[1][3][segdc2Nr]) * (params_sigmas[1][3][segdc2Nr]));
1212
1213 double sigma2_dax = 9.
1214 * ((params_sigmas[1][0][segdc2Nr] * params_sigmas[1][0][segdc2Nr])
1215 + (params_sigmas[0][0][segdc1Nr] * params_sigmas[0][0][segdc1Nr]));
1216
1217 double sigma2_day = 9.
1218 * ((params_sigmas[1][2][segdc2Nr] * params_sigmas[1][2][segdc2Nr])
1219 + (params_sigmas[0][2][segdc1Nr] * params_sigmas[0][2][segdc1Nr]));
1220
1221 double chi2_Match = (((dx * dx) / sigma2_dx) + ((dy * dy) / sigma2_dy) + ((daX * daX) / sigma2_dax)
1222 + ((daY * daY) / sigma2_day))
1223 / (segment_size[0][segdc1Nr] + segment_size[1][segdc2Nr] - 8);
1224
1225 if ((segment_size[0][segdc1Nr] + segment_size[1][segdc2Nr]) < 12)
1226 continue;
1227 if (fabs(dx) > 20. || fabs(dy) > 20. || fabs(daX) > 0.25 || fabs(daY) > 0.25)
1228 continue;
1229 if (chi2_Match < chi2_Match_min) {
1230 chi2_Match_min = chi2_Match;
1231 best1 = segdc1Nr;
1232 best2 = segdc2Nr;
1233 hX1_X2_glob_m->Fill(dx);
1234 hY1_Y2_glob_m->Fill(dy);
1235 haX2_aX1m->Fill(daX);
1236 haY2_aY1m->Fill(daY);
1237 }
1238 } // segdc1Nr
1239
1240 } // segdc2Nr
1241
1242 if (chi2_Match_min > 1199.)
1243 break;
1244 nSegmentsToBeMatched++;
1245 pairID[nSegmentsToBeMatched] = 1000 * (best1 + 1) + best2 + 1;
1246 nhits[nSegmentsToBeMatched] = segment_size[0][best1] + segment_size[1][best2];
1247 x_mid[nSegmentsToBeMatched] =
1248 0.5 * (x_global[0][best1] + x_global[1][best2]); // par_ab1[1][best1];//xDC1_glob[best1];//xMean;
1249 y_mid[nSegmentsToBeMatched] = 0.5 * (y_global[0][best1] + y_global[1][best2]); // yDC1_glob[best1];//yMean;
1250 aX[nSegmentsToBeMatched] = -(params[1][1][best2] - params[0][1][best1]) / dZ_dch;
1251 aY[nSegmentsToBeMatched] = (params[1][3][best2] - params[0][3][best1]) / dZ_dch;
1252 Chi2_match[nSegmentsToBeMatched] = chi2_Match_min;
1253 haX->Fill(aX[nSegmentsToBeMatched]);
1254 haY->Fill(aY[nSegmentsToBeMatched]);
1255 hX->Fill(x_mid[nSegmentsToBeMatched]);
1256 hY->Fill(y_mid[nSegmentsToBeMatched]);
1257 // hP->Fill(0.3*2.84/(atan(aX[nSegmentsToBeMatched])));
1258 // set of field integrals to choose from
1259 // float intBL = 0.984;//600A 818
1260 // float intBL = 1.455;//900A 816 817
1261 // float intBL = 2.725;//1750A
1262 // float intBL = 1.93;//1200A 815
1263 // float intBL = 2.38;//1500A
1264 // float intBL = 3.38;//2100A
1265 // float intBL = 2.55;//1600A 814 834
1266 // float intBL = 2.84;//1800A 813
1267 hChi2Matched->Fill(chi2_Match_min);
1268 Chi2[0][best1] = 999;
1269 Chi2[1][best2] = 999;
1270
1271 } // while
1272}
1273
1274Int_t BmnDchTrackFinder::BuildXYSegments(Int_t dchID,
1275 Int_t pairU,
1276 Int_t pairV,
1277 Int_t pairX,
1278 Int_t pairY,
1279 Int_t single_xa,
1280 Int_t single_xb,
1281 Int_t single_ya,
1282 Int_t single_yb,
1283 Double_t** x_ab,
1284 Double_t** y_ab,
1285 Double_t** u_ab,
1286 Double_t** v_ab,
1287 Int_t** x_abId,
1288 Int_t** y_abId,
1289 Int_t** u_abId,
1290 Int_t** v_abId,
1291 Double_t** x_ab_time,
1292 Double_t** y_ab_time,
1293 Double_t** u_ab_time,
1294 Double_t** v_ab_time,
1295 Double_t** sigm_x_ab,
1296 Double_t** sigm_y_ab,
1297 Double_t** sigm_u_ab,
1298 Double_t** sigm_v_ab,
1299 Double_t** rh_seg,
1300 Int_t** rhId_seg,
1301 Double_t** rh_seg_time,
1302 Double_t** rh_sigm_seg,
1303 Double_t** x_single,
1304 Double_t** y_single,
1305 Int_t** x_singleId,
1306 Int_t** y_singleId,
1307 Double_t** sigm_x_single,
1308 Double_t** sigm_y_single)
1309{
1310 Double_t sqrt_2 = sqrt(2.);
1311 Double_t isqrt_2 = 1. / sqrt_2;
1312
1313 Int_t nDC_segments = (dchID == 1) ? nSegments[0] : nSegments[1];
1314
1315 for (Int_t i = 0; i < pairU; i++) {
1316 if (nDC_segments > 25 * N - 1)
1317 break;
1318 Double_t u_coord = (u_ab[0][i] + u_ab[1][i]) / 2;
1319 Double_t u_slope = (u_ab[0][i] - u_ab[1][i]) / (z_loc[4] - z_loc[5]);
1320 Double_t UX = u_coord + u_slope * 0.5 * (z_loc[1] + z_loc[0] - z_loc[5] - z_loc[4]);
1321 Double_t UY = u_coord + u_slope * 0.5 * (z_loc[3] + z_loc[2] - z_loc[5] - z_loc[4]);
1322
1323 for (Int_t j = 0; j < pairV; j++) {
1324 if (nDC_segments > 25 * N - 1)
1325 break;
1326
1327 Double_t v_coord = (v_ab[0][j] + v_ab[1][j]) / 2;
1328 if (fabs(u_coord) < 12. && fabs(v_coord) < 12.)
1329 continue;
1330 if (fPeriod == 7 && fRunId > 3589) {
1331 if (dchID == 1) {
1332 if (fabs(u_coord - 30) < 12 && fabs(v_coord + 22) < 12)
1333 continue; // cut on beam for Ar
1334 }
1335 if (dchID == 2) {
1336 if (fabs(u_coord - 37) < 14 && fabs(v_coord + 25) < 16)
1337 continue; // cut on beam for Ar
1338 }
1339 }
1340 Double_t v_slope = (v_ab[0][j] - v_ab[1][j]) / (z_loc[6] - z_loc[7]);
1341 Double_t VX = v_coord + v_slope * 0.5 * (z_loc[1] + z_loc[0] - z_loc[6] - z_loc[7]);
1342 Double_t VY = v_coord + v_slope * 0.5 * (z_loc[3] + z_loc[2] - z_loc[6] - z_loc[7]);
1343 Bool_t foundX = kFALSE;
1344 Double_t x_est = isqrt_2 * (VX - UX);
1345 Double_t y_est = isqrt_2 * (UY + VY);
1346 if (pairX > 0) {
1347 Double_t dX_thresh = 1.2;
1348 for (Int_t k = 0; k < pairX; k++) {
1349 Double_t x_coord = (x_ab[0][k] + x_ab[1][k]) / 2;
1350 if (nDC_segments > 25 * N - 1)
1351 break;
1352 if (dchID == 1)
1353 hXm_Xe_loc->Fill(x_coord - x_est);
1354 if (Abs(x_coord - x_est) > dX_thresh)
1355 continue;
1356 dX_thresh = Abs(x_coord - x_est);
1357
1358 rh_seg[0][nDC_segments] = x_ab[0][k];
1359 rh_seg[1][nDC_segments] = x_ab[1][k];
1360 rh_seg[4][nDC_segments] = u_ab[0][i];
1361 rh_seg[5][nDC_segments] = u_ab[1][i];
1362 rh_seg[6][nDC_segments] = v_ab[0][j];
1363 rh_seg[7][nDC_segments] = v_ab[1][j];
1364 rhId_seg[0][nDC_segments] = x_abId[0][k];
1365 rhId_seg[1][nDC_segments] = x_abId[1][k];
1366 rhId_seg[4][nDC_segments] = u_abId[0][i];
1367 rhId_seg[5][nDC_segments] = u_abId[1][i];
1368 rhId_seg[6][nDC_segments] = v_abId[0][j];
1369 rhId_seg[7][nDC_segments] = v_abId[1][j];
1370 rh_seg_time[0][nDC_segments] = x_ab_time[0][k];
1371 rh_seg_time[1][nDC_segments] = x_ab_time[1][k];
1372 rh_seg_time[4][nDC_segments] = u_ab_time[0][i];
1373 rh_seg_time[5][nDC_segments] = u_ab_time[1][i];
1374 rh_seg_time[6][nDC_segments] = v_ab_time[0][j];
1375 rh_seg_time[7][nDC_segments] = v_ab_time[1][j];
1376 rh_sigm_seg[0][nDC_segments] = sigm_x_ab[0][k];
1377 rh_sigm_seg[1][nDC_segments] = sigm_x_ab[1][k];
1378 rh_sigm_seg[4][nDC_segments] = sigm_u_ab[0][i];
1379 rh_sigm_seg[5][nDC_segments] = sigm_u_ab[1][i];
1380 rh_sigm_seg[6][nDC_segments] = sigm_v_ab[0][j];
1381 rh_sigm_seg[7][nDC_segments] = sigm_v_ab[1][j];
1382
1383 foundX = kTRUE;
1384 if (nDC_segments > 25 * N - 1)
1385 break;
1386 } // k
1387 } //(pair_x2>0)
1388
1389 Bool_t foundY = kFALSE;
1390 if (pairY > 0) {
1391 Double_t dY_thresh = 1.2;
1392 for (Int_t m = 0; m < pairY; m++) {
1393 if (nDC_segments > 25 * N - 1)
1394 break;
1395 Double_t y_coord = (y_ab[0][m] + y_ab[1][m]) / 2;
1396 if (dchID == 1)
1397 hYm_Ye_loc->Fill(y_coord - y_est);
1398 if (Abs(y_coord - y_est) > dY_thresh)
1399 continue;
1400 dY_thresh = Abs(y_coord - y_est);
1401 foundY = kTRUE;
1402 rh_seg[2][nDC_segments] = y_ab[0][m];
1403 rh_seg[3][nDC_segments] = y_ab[1][m];
1404 rh_seg[4][nDC_segments] = u_ab[0][i];
1405 rh_seg[5][nDC_segments] = u_ab[1][i];
1406 rh_seg[6][nDC_segments] = v_ab[0][j];
1407 rh_seg[7][nDC_segments] = v_ab[1][j];
1408 rhId_seg[2][nDC_segments] = y_abId[0][m];
1409 rhId_seg[3][nDC_segments] = y_abId[1][m];
1410 rhId_seg[4][nDC_segments] = u_abId[0][i];
1411 rhId_seg[5][nDC_segments] = u_abId[1][i];
1412 rhId_seg[6][nDC_segments] = v_abId[0][j];
1413 rhId_seg[7][nDC_segments] = v_abId[1][j];
1414 rh_seg_time[2][nDC_segments] = y_ab_time[0][m];
1415 rh_seg_time[3][nDC_segments] = y_ab_time[1][m];
1416 rh_seg_time[4][nDC_segments] = u_ab_time[0][i];
1417 rh_seg_time[5][nDC_segments] = u_ab_time[1][i];
1418 rh_seg_time[6][nDC_segments] = v_ab_time[0][j];
1419 rh_seg_time[7][nDC_segments] = v_ab_time[1][j];
1420 rh_sigm_seg[2][nDC_segments] = sigm_y_ab[0][m];
1421 rh_sigm_seg[3][nDC_segments] = sigm_y_ab[1][m];
1422 rh_sigm_seg[4][nDC_segments] = sigm_u_ab[0][i];
1423 rh_sigm_seg[5][nDC_segments] = sigm_u_ab[1][i];
1424 rh_sigm_seg[6][nDC_segments] = sigm_v_ab[0][j];
1425 rh_sigm_seg[7][nDC_segments] = sigm_v_ab[1][j];
1426 if (!foundX) {
1427 Double_t min_a = 999;
1428 Double_t min_b = 999;
1429 for (Int_t kk = 0; kk < single_xa; kk++) {
1430 if (Abs(x_single[0][kk] - x_est) > 1.2)
1431 continue; //????? 0.5 needs to be reviewed
1432
1433 if (Abs(x_single[0][kk] - x_est) < min_a) {
1434 min_a = Abs(x_single[0][kk] - x_est);
1435 rh_seg[0][nDC_segments] = x_single[0][kk];
1436 rhId_seg[0][nDC_segments] = x_singleId[0][kk];
1437 rh_sigm_seg[0][nDC_segments] = sigm_x_single[0][kk];
1438 foundX = kTRUE;
1439 }
1440 } // for kk
1441 for (Int_t kk = 0; kk < single_xb; kk++) {
1442 if (Abs(x_single[1][kk] - x_est) > 1.2)
1443 continue; //????? 0.5 needs to be reviewed
1444 if (Abs(x_single[1][kk] - x_est) < min_b) {
1445 min_b = Abs(x_single[1][kk] - x_est);
1446 rh_seg[1][nDC_segments] = x_single[1][kk];
1447 rhId_seg[1][nDC_segments] = x_singleId[1][kk];
1448 rh_sigm_seg[1][nDC_segments] = sigm_x_single[1][kk];
1449 foundX = kTRUE;
1450 }
1451 } // for kk
1452 if (nDC_segments > 25 * N - 1)
1453 break;
1454 } //! foundX
1455 } // m
1456 if (foundX && !foundY) {
1457 Double_t min_a = 999;
1458 Double_t min_b = 999;
1459 for (Int_t kk = 0; kk < single_ya; kk++) {
1460 if (Abs(y_single[0][kk] - y_est) > 1.2)
1461 continue; //????? 0.5 needs to be reviewed
1462 if (Abs(y_single[0][kk] - y_est) < min_a) {
1463 min_a = Abs(y_single[0][kk] - y_est);
1464 rh_seg[2][nDC_segments] = y_single[0][kk];
1465 rhId_seg[2][nDC_segments] = y_singleId[0][kk];
1466 rh_sigm_seg[2][nDC_segments] = sigm_y_single[0][kk];
1467 foundY = kTRUE;
1468 }
1469 } // for kk
1470 for (Int_t kk = 0; kk < single_yb; kk++) {
1471 if (Abs(y_single[1][kk] - y_est) > 1.2)
1472 continue; //????? 0.5 needs to be reviewed
1473 if (Abs(y_single[1][kk] - y_est) < min_b) {
1474 min_b = Abs(y_single[1][kk] - y_est);
1475 rh_seg[3][nDC_segments] = y_single[1][kk];
1476 rhId_seg[3][nDC_segments] = y_singleId[1][kk];
1477 rh_sigm_seg[3][nDC_segments] = sigm_y_single[1][kk];
1478 foundY = kTRUE;
1479 }
1480 } // for kk
1481 }
1482 } //(pair_y2>0)
1483 if (foundX || foundY)
1484 nDC_segments++;
1485 }
1486 }
1487 return nDC_segments;
1488}
1489
1490Int_t BmnDchTrackFinder::BuildUVSegments(Int_t dchID,
1491 Int_t pairU,
1492 Int_t pairV,
1493 Int_t pairX,
1494 Int_t pairY,
1495 Int_t single_ua,
1496 Int_t single_ub,
1497 Int_t single_va,
1498 Int_t single_vb,
1499 Double_t** x_ab,
1500 Double_t** y_ab,
1501 Double_t** u_ab,
1502 Double_t** v_ab,
1503 Int_t** x_abId,
1504 Int_t** y_abId,
1505 Int_t** u_abId,
1506 Int_t** v_abId,
1507 Double_t** x_ab_time,
1508 Double_t** y_ab_time,
1509 Double_t** u_ab_time,
1510 Double_t** v_ab_time,
1511 Double_t** sigm_x_ab,
1512 Double_t** sigm_y_ab,
1513 Double_t** sigm_u_ab,
1514 Double_t** sigm_v_ab,
1515 Double_t** rh_seg,
1516 Int_t** rhId_seg,
1517 Double_t** rh_seg_time,
1518 Double_t** rh_sigm_seg,
1519 Double_t** u_single,
1520 Double_t** v_single,
1521 Int_t** u_singleId,
1522 Int_t** v_singleId,
1523 Double_t** sigm_u_single,
1524 Double_t** sigm_v_single)
1525{
1526 Double_t sqrt_2 = sqrt(2.);
1527 Double_t isqrt_2 = 1. / sqrt_2;
1528
1529 Int_t nDC_segments = (dchID == 1) ? nSegments[0] : nSegments[1];
1530
1531 for (Int_t i = 0; i < pairX; i++) {
1532 if (nDC_segments > 25 * N - 1)
1533 break;
1534 Double_t x_coord = (x_ab[0][i] + x_ab[1][i]) / 2;
1535 Double_t x_slope = (x_ab[0][i] - x_ab[1][i]) / (z_loc[0] - z_loc[1]);
1536 Double_t XU = x_coord + x_slope * 0.5 * (z_loc[5] + z_loc[4] - z_loc[1] - z_loc[0]);
1537 Double_t XV = x_coord + x_slope * 0.5 * (z_loc[7] + z_loc[6] - z_loc[1] - z_loc[0]);
1538
1539 for (Int_t j = 0; j < pairY; j++) {
1540 Double_t y_coord = (y_ab[0][j] + y_ab[1][j]) / 2;
1541 if (fPeriod == 7 && fRunId > 3589) {
1542 if (dchID == 1) {
1543 if (fabs(x_coord + 37) < 12 && fabs(y_coord - 5) < 12)
1544 continue; // cut on beam for Ar
1545 }
1546 if (dchID == 2) {
1547 if (fabs(x_coord + 45) < 12 && fabs(y_coord - 8) < 12)
1548 continue; // cut on beam for Ar
1549 }
1550 }
1551 if (fabs(y_coord) < 12. && fabs(x_coord) < 12.)
1552 continue;
1553 Double_t y_slope = (y_ab[0][j] - y_ab[1][j]) / (z_loc[2] - z_loc[3]);
1554 Double_t YU = y_coord + y_slope * 0.5 * (z_loc[5] + z_loc[4] - z_loc[3] - z_loc[2]);
1555 Double_t YV = y_coord + y_slope * 0.5 * (z_loc[7] + z_loc[6] - z_loc[3] - z_loc[2]);
1556 Bool_t foundU = kFALSE;
1557 Double_t u_est = isqrt_2 * (YU - XU);
1558 Double_t v_est = isqrt_2 * (YV + XV);
1559
1560 Double_t dU_thresh = 1.2;
1561 for (Int_t k = 0; k < pairU; k++) {
1562 Double_t u_coord = (u_ab[0][k] + u_ab[1][k]) / 2;
1563 if (dchID == 1)
1564 hUm_Ue_loc->Fill(u_coord - u_est);
1565 if (Abs(u_coord - u_est) > dU_thresh)
1566 continue;
1567 dU_thresh = Abs(u_coord - u_est);
1568
1569 rh_seg[0][nDC_segments] = x_ab[0][i];
1570 rh_seg[1][nDC_segments] = x_ab[1][i];
1571 rh_seg[2][nDC_segments] = y_ab[0][j];
1572 rh_seg[3][nDC_segments] = y_ab[1][j];
1573 rh_seg[4][nDC_segments] = u_ab[0][k];
1574 rh_seg[5][nDC_segments] = u_ab[1][k];
1575 rhId_seg[0][nDC_segments] = x_abId[0][i];
1576 rhId_seg[1][nDC_segments] = x_abId[1][i];
1577 rhId_seg[2][nDC_segments] = y_abId[0][j];
1578 rhId_seg[3][nDC_segments] = y_abId[1][j];
1579 rhId_seg[4][nDC_segments] = u_abId[0][k];
1580 rhId_seg[5][nDC_segments] = u_abId[1][k];
1581 rh_seg_time[0][nDC_segments] = x_ab_time[0][i];
1582 rh_seg_time[1][nDC_segments] = x_ab_time[1][i];
1583 rh_seg_time[2][nDC_segments] = y_ab_time[0][j];
1584 rh_seg_time[3][nDC_segments] = y_ab_time[1][j];
1585 rh_seg_time[4][nDC_segments] = u_ab_time[0][k];
1586 rh_seg_time[5][nDC_segments] = u_ab_time[1][k];
1587 rh_sigm_seg[0][nDC_segments] = sigm_x_ab[0][i];
1588 rh_sigm_seg[1][nDC_segments] = sigm_x_ab[1][i];
1589 rh_sigm_seg[2][nDC_segments] = sigm_y_ab[0][j];
1590 rh_sigm_seg[3][nDC_segments] = sigm_y_ab[1][j];
1591 rh_sigm_seg[4][nDC_segments] = sigm_u_ab[0][k];
1592 rh_sigm_seg[5][nDC_segments] = sigm_u_ab[1][k];
1593
1594 foundU = kTRUE;
1595 if (nDC_segments > 25 * N - 1)
1596 break;
1597 }
1598
1599 Bool_t foundV = kFALSE;
1600
1601 Double_t dV_thresh = 1.2;
1602 for (Int_t m = 0; m < pairV; m++) {
1603 if (nDC_segments > 25 * N - 1)
1604 break;
1605 Double_t v_coord = (v_ab[0][m] + v_ab[1][m]) / 2;
1606 if (dchID == 1)
1607 hVm_Ve_loc->Fill(v_coord - v_est);
1608 if (Abs(v_coord - v_est) > dV_thresh)
1609 continue;
1610 dV_thresh = Abs(v_coord - v_est);
1611
1612 foundV = kTRUE;
1613 rh_seg[0][nDC_segments] = x_ab[0][i];
1614 rh_seg[1][nDC_segments] = x_ab[1][i];
1615 rh_seg[2][nDC_segments] = y_ab[0][j];
1616 rh_seg[3][nDC_segments] = y_ab[1][j];
1617 rh_seg[6][nDC_segments] = v_ab[0][m];
1618 rh_seg[7][nDC_segments] = v_ab[1][m];
1619 rhId_seg[0][nDC_segments] = x_abId[0][i];
1620 rhId_seg[1][nDC_segments] = x_abId[1][i];
1621 rhId_seg[2][nDC_segments] = y_abId[0][j];
1622 rhId_seg[3][nDC_segments] = y_abId[1][j];
1623 rhId_seg[6][nDC_segments] = v_abId[0][m];
1624 rhId_seg[7][nDC_segments] = v_abId[1][m];
1625 rh_seg_time[0][nDC_segments] = x_ab_time[0][i];
1626 rh_seg_time[1][nDC_segments] = x_ab_time[1][i];
1627 rh_seg_time[2][nDC_segments] = y_ab_time[0][j];
1628 rh_seg_time[3][nDC_segments] = y_ab_time[1][j];
1629 rh_seg_time[6][nDC_segments] = v_ab_time[0][m];
1630 rh_seg_time[7][nDC_segments] = v_ab_time[1][m];
1631 rh_sigm_seg[0][nDC_segments] = sigm_x_ab[0][i];
1632 rh_sigm_seg[1][nDC_segments] = sigm_x_ab[1][i];
1633 rh_sigm_seg[2][nDC_segments] = sigm_y_ab[0][j];
1634 rh_sigm_seg[3][nDC_segments] = sigm_y_ab[1][j];
1635 rh_sigm_seg[6][nDC_segments] = sigm_v_ab[0][m];
1636 rh_sigm_seg[7][nDC_segments] = sigm_v_ab[1][m];
1637
1638 if (!foundU) {
1639 Double_t min_a = 999;
1640 Double_t min_b = 999;
1641 for (Int_t kk = 0; kk < single_ua; kk++) {
1642 if (Abs(u_single[0][kk] - u_est) > 1.2)
1643 continue; //????? 0.5 needs to be reviewed
1644 if (Abs(u_single[0][kk] - u_est) < min_a) {
1645 min_a = Abs(u_single[0][kk] - u_est);
1646 rh_seg[4][nDC_segments] = u_single[0][kk];
1647 rhId_seg[4][nDC_segments] = u_singleId[0][kk];
1648 rh_sigm_seg[4][nDC_segments] = sigm_u_single[0][kk];
1649 foundU = kTRUE;
1650 }
1651 } // for kk
1652 for (Int_t kk = 0; kk < single_ub; kk++) {
1653 if (Abs(u_single[1][kk] - u_est) > 1.2)
1654 continue; //????? 0.5 needs to be reviewed
1655 if (Abs(u_single[1][kk] - u_est) < min_b) {
1656 min_b = Abs(u_single[1][kk] - u_est);
1657 rh_seg[5][nDC_segments] = u_single[1][kk];
1658 rhId_seg[5][nDC_segments] = u_singleId[1][kk];
1659 rh_sigm_seg[5][nDC_segments] = sigm_u_single[1][kk];
1660 foundU = kTRUE;
1661 }
1662 } // for kk
1663 if (nDC_segments > 25 * N - 1)
1664 break;
1665 } //! foundU
1666
1667 if (nDC_segments > 25 * N - 1)
1668 break;
1669 } // m
1670 if (!foundV && foundU) {
1671 Double_t min_a = 999;
1672 Double_t min_b = 999;
1673 for (Int_t kk = 0; kk < single_va; kk++) {
1674 if (Abs(v_single[0][kk] - v_est) > 1.2)
1675 continue; //????? 0.5 needs to be reviewed
1676 if (Abs(v_single[0][kk] - v_est) < min_a) {
1677 min_a = Abs(v_single[0][kk] - v_est);
1678 rh_seg[6][nDC_segments] = v_single[0][kk];
1679 rhId_seg[6][nDC_segments] = v_singleId[0][kk];
1680 rh_sigm_seg[6][nDC_segments] = sigm_v_single[0][kk];
1681 foundV = kTRUE;
1682 }
1683 } // for kk
1684 for (Int_t kk = 0; kk < single_vb; kk++) {
1685 if (Abs(v_single[1][kk] - v_est) > 1.2)
1686 continue; //????? 0.5 needs to be reviewed
1687 if (Abs(v_single[1][kk] - v_est) < min_b) {
1688 min_b = Abs(v_single[1][kk] - v_est);
1689 rh_seg[7][nDC_segments] = v_single[1][kk];
1690 rhId_seg[7][nDC_segments] = v_singleId[1][kk];
1691 rh_sigm_seg[7][nDC_segments] = sigm_v_single[1][kk];
1692 foundV = kTRUE;
1693 }
1694 }
1695 }
1696 if (foundV || foundU)
1697 nDC_segments++;
1698 } // y
1699 }
1700 return nDC_segments;
1701}
1702
1703Bool_t BmnDchTrackFinder::FitDchSegments(Int_t dchID,
1704 Int_t* size_seg,
1705 Double_t** rh_seg,
1706 Double_t** rh_sigm_seg,
1707 Double_t** par_ab,
1708 Double_t* chi2,
1709 Double_t* x_glob,
1710 Double_t* y_glob)
1711{
1712 Int_t nDC_segments = (dchID == 1) ? nSegments[0] : nSegments[1];
1713 Double_t Z_dch = (dchID == 1) ? Z_dch1 : Z_dch2;
1714 Bool_t hasSuffNumberOfSegments = kFALSE;
1715 // Double_t sqrt_2 = sqrt(2.);
1716 // Double_t isqrt_2 = 1. / sqrt_2;
1717 for (Int_t j = 0; j < nDC_segments; j++) {
1718 Int_t worst_hit = -1;
1719 Double_t max_resid = 0;
1720
1721 Double_t _rh_seg[8];
1722 Double_t _rh_sigm_seg[8];
1723 Double_t _par_ab[4];
1724
1725 for (Int_t i = 0; i < 8; i++)
1726 if (Abs(rh_seg[i][j] + 999.) > FLT_EPSILON)
1727 size_seg[j]++;
1728
1729 for (Int_t rej = 0; rej < 2; rej++) { // allow 2 passes max 8->7 & 7->6
1730 for (Int_t i = 0; i < 8; i++) {
1731 _rh_seg[i] = rh_seg[i][j];
1732 _rh_sigm_seg[i] = rh_sigm_seg[i][j];
1733 }
1734
1735 fit_seg(z_loc, _rh_seg, _rh_sigm_seg, _par_ab, -1, -1); // usual fit without skipping any plane
1736 for (Int_t i = 0; i < 4; i++)
1737 par_ab[i][j] = _par_ab[i];
1738
1739 chi2[j] = 0;
1740
1741 Double_t resid(LDBL_MAX);
1742 for (Int_t i = 0; i < 8; i++) {
1743 if (Abs(rh_seg[i][j] + 999.) < FLT_EPSILON)
1744 continue;
1745
1746 resid = CalculateResidual(i, j, rh_seg, par_ab);
1747 chi2[j] += (resid * resid) / rh_sigm_seg[i][j];
1748 if (Abs(resid) > max_resid) {
1749 worst_hit = i;
1750 max_resid = Abs(resid);
1751 }
1752 }
1753
1754 chi2[j] /= (size_seg[j] - 4);
1755 // if chi2 is big and seg_size = min erase this seg
1756 if (chi2[j] > 180.) {
1757 if (size_seg[j] == 6) {
1758 chi2[j] = 999.;
1759 break;
1760 } else {
1761 rh_seg[worst_hit][j] = -999.; // erase worst hit and refit
1762 size_seg[j]--;
1763 max_resid = 0;
1764 continue;
1765 }
1766 }
1767 }
1768
1769 // Add shifts to slopes and coords
1770 Double_t x_slope_sh = (dchID == 1) ? x1_slope_sh : x2_slope_sh;
1771 Double_t y_slope_sh = (dchID == 1) ? y1_slope_sh : y2_slope_sh;
1772 Double_t x_sh = (dchID == 1) ? x1_sh : x2_sh;
1773 Double_t y_sh = (dchID == 1) ? y1_sh : y2_sh;
1774 // the following 4 lines to be commented for residual calculation
1775 par_ab[0][j] += x_slope_sh + x_slope_sh * par_ab[0][j] * par_ab[0][j];
1776 par_ab[2][j] += y_slope_sh + y_slope_sh * par_ab[2][j] * par_ab[2][j];
1777 par_ab[1][j] += x_sh;
1778 par_ab[3][j] += y_sh;
1779
1780 // cout<<" par_ab[1][j] "<< par_ab[1][j]<<endl;
1781 if (dchID == 1) {
1782 hXDC1_atZ0->Fill(par_ab[0][j] * (-Z_dch + 249) + par_ab[1][j]);
1783 hYDC1_atZ0->Fill(par_ab[2][j] * (-Z_dch + 249) + par_ab[3][j]);
1784 } else {
1785 hXDC2_atZ0->Fill(par_ab[0][j] * (-Z_dch + 249) + par_ab[1][j]);
1786 hYDC2_atZ0->Fill(par_ab[2][j] * (-Z_dch + 249) + par_ab[3][j]);
1787 }
1788 // cut on x0 and y0 at back of the magnet
1789 if (expData) {
1790 if ((par_ab[2][j] * (-Z_dch + 249) + par_ab[3][j]) < -50
1791 || (par_ab[2][j] * (-Z_dch + 249) + par_ab[3][j]) > 50
1792 || (par_ab[0][j] * (-Z_dch + 249) + par_ab[1][j]) < -85
1793 || (par_ab[0][j] * (-Z_dch + 249) + par_ab[1][j]) > 85)
1794 {
1795 chi2[j] = 999;
1796 }
1797 }
1798 Int_t coeff = (dchID == 1) ? 1 : -1;
1799 x_glob[j] = coeff * par_ab[0][j] * dZ_dch_mid + par_ab[1][j];
1800 y_glob[j] = coeff * par_ab[2][j] * dZ_dch_mid + par_ab[3][j];
1801 if (size_seg[j] > 5)
1802 hasSuffNumberOfSegments = kTRUE;
1803 }
1804 return hasSuffNumberOfSegments;
1805}
1806
1807void BmnDchTrackFinder::CompareDaDb(Double_t d, Double_t& ele)
1808{
1809 ele = (d < 0.02) ? (0.04 * 0.04)
1810 : (d >= 0.02 && d < 0.1) ? (0.03 * 0.03)
1811 : (d >= 0.1 && d < 0.4) ? (0.015 * 0.015)
1812 : (d >= 0.4 && d < 0.41) ? (0.04 * 0.04)
1813 : (0.05 * 0.05); // 2019.11.5
1814}
1815
1816void BmnDchTrackFinder::CompareDaDb(Double_t d, Double_t& ele1, Double_t& ele2)
1817{
1818 ele1 = (d < 0.02) ? (0.04 * 0.04)
1819 : (d >= 0.02 && d < 0.1) ? (0.03 * 0.03)
1820 : (d >= 0.1 && d < 0.4) ? (0.015 * 0.015)
1821 : (d >= 0.4 && d < 0.41) ? (0.04 * 0.04)
1822 : (0.05 * 0.05); // 2019.11.5
1823 ele2 = ele1;
1824}
1825
1826void BmnDchTrackFinder::SelectLongestAndBestSegments(Int_t dchID, Int_t* size_seg, Double_t** rh_seg, Double_t* chi2)
1827{
1828 Int_t nDC_segments = (dchID == 1) ? nSegments[0] : nSegments[1];
1829 for (Int_t max_size = 8; max_size > 5; max_size--) {
1830 for (Int_t it1 = 0; it1 < nDC_segments; it1++) {
1831 if (size_seg[it1] != max_size || chi2[it1] > 990.)
1832 continue;
1833 for (Int_t it2 = 0; it2 < nDC_segments; it2++) {
1834 if (it2 == it1 || chi2[it2] > 990.)
1835 continue;
1836 for (Int_t hit = 0; hit < 4; hit++) {
1837 if ((rh_seg[2 * hit][it1] == rh_seg[2 * hit][it2]
1838 || rh_seg[2 * hit + 1][it1] == rh_seg[2 * hit + 1][it2])
1839 && (rh_seg[2 * hit][it1] > -998 && rh_seg[2 * hit + 1][it1] > -998))
1840 {
1841 if (size_seg[it1] == size_seg[it2]) {
1842 if (chi2[it1] <= chi2[it2]) {
1843 chi2[it2] = 999.;
1844 break;
1845 } else {
1846 chi2[it1] = 999.;
1847 break;
1848 }
1849 } // SAME SIZE
1850 if (size_seg[it1] > size_seg[it2]) {
1851 chi2[it2] = 999.;
1852 break;
1853 }
1854 } // have common hits
1855 } // rh cycle
1856 } // SECOND CYCLE
1857 } // FIRST CYCLE
1858 } // SCAN SEG SIZE
1859}
1860
1861void BmnDchTrackFinder::FillPlaneResiduals(Int_t dchID,
1862 Int_t* size_seg,
1863 Double_t** rh_seg,
1864 Double_t** rh_sigm_seg,
1865 Double_t** par_ab,
1866 Double_t* chi2)
1867{
1868 Int_t nDC_segments = (dchID == 1) ? nSegments[0] : nSegments[1];
1869 for (Int_t it1 = 0; it1 < nDC_segments; it1++) {
1870 if (chi2[it1] > 990.)
1871 continue;
1872 Double_t _rh_seg[8];
1873 Double_t _rh_sigm_seg[8];
1874 Double_t _par_ab[4];
1875
1876 for (Int_t i = 0; i < 8; i++) {
1877 _rh_seg[i] = rh_seg[i][it1];
1878 _rh_sigm_seg[i] = rh_sigm_seg[i][it1];
1879 }
1880
1881 fit_seg(z_loc, _rh_seg, _rh_sigm_seg, _par_ab, -1, -1); // usual fit without skipping any plane
1882 for (Int_t i = 0; i < 4; i++)
1883 par_ab[i][it1] = _par_ab[i];
1884 // fill residuals
1885 if (size_seg[it1] == 8) {
1886 if (dchID == 1) {
1887
1888 double res1 = CalculateResidual(0, it1, rh_seg, par_ab);
1889 double res2 = CalculateResidual(1, it1, rh_seg, par_ab);
1890
1891 hResidx1a->Fill(res1);
1892 hResidx1b->Fill(res2);
1893
1894 double int1;
1895 // double int2;
1896 double dist_a = 999;
1897 double dist_b = 999;
1898 if (fabs(modf(rh_seg[0][it1], &int1)) < .5) {
1899 dist_a = fabs(modf(rh_seg[0][it1], &int1));
1900 hDc1XaResVsDa->Fill(dist_a, fabs(res1), 1);
1901 if (res1 < 0.)
1902 hDc1XaResVsDa_m->Fill(dist_a, -res1, 1);
1903 if (res1 > 0.)
1904 hDc1XaResVsDa_p->Fill(dist_a, res1, 1);
1905 } else {
1906 dist_a = 1 - fabs(modf(rh_seg[0][it1], &int1));
1907 hDc1XaResVsDa->Fill(dist_a, fabs(res1), 1);
1908 if (res1 < 0.)
1909 hDc1XaResVsDa_m->Fill(dist_a, -res1, 1);
1910 if (res1 > 0.)
1911 hDc1XaResVsDa_p->Fill(dist_a, res1, 1);
1912 }
1913 if (fabs(modf(rh_seg[1][it1], &int1)) < .5) {
1914 dist_b = 0.5 - fabs(modf(rh_seg[1][it1], &int1));
1915 hDc1XbResVsDb->Fill(dist_b, fabs(res2), 1);
1916 if (res2 < 0.)
1917 hDc1XbResVsDb_m->Fill(dist_b, -res2, 1);
1918 if (res2 > 0.)
1919 hDc1XbResVsDb_p->Fill(dist_b, res2, 1);
1920 } else {
1921 dist_b = fabs(modf(rh_seg[1][it1], &int1)) - .5;
1922 hDc1XbResVsDb->Fill(dist_b, fabs(res2), 1);
1923 if (res2 < 0.)
1924 hDc1XbResVsDb_m->Fill(dist_b, -res2, 1);
1925 if (res2 > 0.)
1926 hDc1XbResVsDb_p->Fill(dist_b, res2, 1);
1927 }
1928 if (dist_a < .1) {
1929 hResidx1a_0p1->Fill(res1);
1930 hPullsx1a_0p1->Fill(res1 / sqrt(rh_sigm_seg[0][it1]));
1931 }
1932 if (dist_a > .1 && dist_a < .4) {
1933 hResidx1a_0p1_0p4->Fill(res1);
1934 hPullsx1a_0p1_0p4->Fill(res1 / sqrt(rh_sigm_seg[0][it1]));
1935 }
1936 if (dist_a > .4) {
1937 hResidx1a_0p4_0p5->Fill(res1);
1938 hPullsx1a_0p4_0p5->Fill(res1 / sqrt(rh_sigm_seg[0][it1]));
1939 }
1940 if (dist_b < .1) {
1941 hResidx1b_0p1->Fill(res2);
1942 hPullsx1b_0p1->Fill(res2 / sqrt(rh_sigm_seg[1][it1]));
1943 }
1944 if (dist_b > .1 && dist_b < .4) {
1945 hResidx1b_0p1_0p4->Fill(res2);
1946 hPullsx1b_0p1_0p4->Fill(res2 / sqrt(rh_sigm_seg[1][it1]));
1947 }
1948 if (dist_b > .4) {
1949 hResidx1b_0p4_0p5->Fill(res2);
1950 hPullsx1b_0p4_0p5->Fill(res2 / sqrt(rh_sigm_seg[1][it1]));
1951 }
1952 hx1Da_Db_best->Fill(dist_a + dist_b - 0.5);
1953
1954 // y
1955 dist_a = 999;
1956 dist_b = 999;
1957 res1 = CalculateResidual(2, it1, rh_seg, par_ab);
1958 res2 = CalculateResidual(3, it1, rh_seg, par_ab);
1959 hResidy1a->Fill(res1);
1960 hResidy1b->Fill(res2);
1961 if (fabs(modf(rh_seg[2][it1], &int1)) < .5) {
1962 dist_a = fabs(modf(rh_seg[2][it1], &int1));
1963 } else {
1964 dist_a = 1 - fabs(modf(rh_seg[2][it1], &int1));
1965 }
1966 if (fabs(modf(rh_seg[3][it1], &int1)) < .5) {
1967 dist_b = 0.5 - fabs(modf(rh_seg[3][it1], &int1));
1968 } else {
1969 dist_b = fabs(modf(rh_seg[3][it1], &int1)) - .5;
1970 }
1971 if (dist_a < .1)
1972 hResidy1a_0p1->Fill(res1);
1973 if (dist_a > .1 && dist_a < .4)
1974 hResidy1a_0p1_0p4->Fill(res1);
1975 if (dist_a > .4)
1976 hResidy1a_0p4_0p5->Fill(res1);
1977 if (dist_b < .1)
1978 hResidy1b_0p1->Fill(res2);
1979 if (dist_b > .1 && dist_b < .4)
1980 hResidy1b_0p1_0p4->Fill(res2);
1981 if (dist_b > .4)
1982 hResidy1b_0p4_0p5->Fill(res2);
1983 hy1Da_Db_best->Fill(dist_a + dist_b - 0.5);
1984 if (dist_a < .1) {
1985 hPullsy1a_0p1->Fill(res1 / sqrt(rh_sigm_seg[2][it1]));
1986 }
1987 if (dist_a > .1 && dist_a < .4) {
1988 hPullsy1a_0p1_0p4->Fill(res1 / sqrt(rh_sigm_seg[2][it1]));
1989 }
1990 if (dist_a > .4) {
1991 hPullsy1a_0p4_0p5->Fill(res1 / sqrt(rh_sigm_seg[2][it1]));
1992 }
1993 if (dist_b < .1) {
1994 hPullsy1b_0p1->Fill(res2 / sqrt(rh_sigm_seg[3][it1]));
1995 }
1996 if (dist_b > .1 && dist_b < .4) {
1997 hPullsy1b_0p1_0p4->Fill(res2 / sqrt(rh_sigm_seg[3][it1]));
1998 }
1999 if (dist_b > .4) {
2000 hPullsy1b_0p4_0p5->Fill(res2 / sqrt(rh_sigm_seg[3][it1]));
2001 }
2002
2003 // u
2004 dist_a = 999;
2005 dist_b = 999;
2006 res1 = CalculateResidual(4, it1, rh_seg, par_ab);
2007 res2 = CalculateResidual(5, it1, rh_seg, par_ab);
2008 hResidu1a->Fill(res1);
2009 hResidu1b->Fill(res2);
2010
2011 if (fabs(modf(rh_seg[4][it1], &int1)) < .5) {
2012 dist_a = fabs(modf(rh_seg[4][it1], &int1));
2013 } else {
2014 dist_a = 1 - fabs(modf(rh_seg[4][it1], &int1));
2015 }
2016 if (fabs(modf(rh_seg[5][it1], &int1)) < .5) {
2017 dist_b = 0.5 - fabs(modf(rh_seg[5][it1], &int1));
2018 } else {
2019 dist_b = fabs(modf(rh_seg[5][it1], &int1)) - .5;
2020 }
2021 if (dist_a < .1)
2022 hResidu1a_0p1->Fill(res1);
2023 if (dist_a > .1 && dist_a < .4)
2024 hResidu1a_0p1_0p4->Fill(res1);
2025 if (dist_a > .4)
2026 hResidu1a_0p4_0p5->Fill(res1);
2027 if (dist_b < .1)
2028 hResidu1b_0p1->Fill(res2);
2029 if (dist_b > .1 && dist_b < .4)
2030 hResidu1b_0p1_0p4->Fill(res2);
2031 if (dist_b > .4)
2032 hResidu1b_0p4_0p5->Fill(res2);
2033 hu1Da_Db_best->Fill(dist_a + dist_b - 0.5);
2034 if (dist_a < .1) {
2035 hPullsu1a_0p1->Fill(res1 / sqrt(rh_sigm_seg[4][it1]));
2036 }
2037 if (dist_a > .1 && dist_a < .4) {
2038 hPullsu1a_0p1_0p4->Fill(res1 / sqrt(rh_sigm_seg[4][it1]));
2039 }
2040 if (dist_a > .4) {
2041 hPullsu1a_0p4_0p5->Fill(res1 / sqrt(rh_sigm_seg[4][it1]));
2042 }
2043 if (dist_b < .1) {
2044 hPullsu1b_0p1->Fill(res2 / sqrt(rh_sigm_seg[5][it1]));
2045 }
2046 if (dist_b > .1 && dist_b < .4) {
2047 hPullsu1b_0p1_0p4->Fill(res2 / sqrt(rh_sigm_seg[5][it1]));
2048 }
2049 if (dist_b > .4) {
2050 hPullsu1b_0p4_0p5->Fill(res2 / sqrt(rh_sigm_seg[5][it1]));
2051 }
2052 // v
2053 dist_a = 999;
2054 dist_b = 999;
2055 res1 = CalculateResidual(6, it1, rh_seg, par_ab);
2056 res2 = CalculateResidual(7, it1, rh_seg, par_ab);
2057 hResidv1a->Fill(res1);
2058 hResidv1b->Fill(res2);
2059 if (fabs(modf(rh_seg[6][it1], &int1)) < .5) {
2060
2061 dist_a = fabs(modf(rh_seg[6][it1], &int1));
2062 } else {
2063 dist_a = 1 - fabs(modf(rh_seg[6][it1], &int1));
2064 }
2065 if (fabs(modf(rh_seg[7][it1], &int1)) < .5) {
2066 dist_b = 0.5 - fabs(modf(rh_seg[7][it1], &int1));
2067 } else {
2068 dist_b = fabs(modf(rh_seg[7][it1], &int1)) - .5;
2069 }
2070 if (dist_a < .1)
2071 hResidv1a_0p1->Fill(res1);
2072 if (dist_a > .1 && dist_a < .4)
2073 hResidv1a_0p1_0p4->Fill(res1);
2074 if (dist_a > .4)
2075 hResidv1a_0p4_0p5->Fill(res1);
2076 if (dist_b < .1)
2077 hResidv1b_0p1->Fill(res2);
2078 if (dist_b > .1 && dist_b < .4)
2079 hResidv1b_0p1_0p4->Fill(res2);
2080 if (dist_b > .4)
2081 hResidv1b_0p4_0p5->Fill(res2);
2082 hv1Da_Db_best->Fill(dist_a + dist_b - 0.5);
2083
2084 if (dist_a < .1) {
2085 hPullsv1a_0p1->Fill(res1 / sqrt(rh_sigm_seg[6][it1]));
2086 }
2087 if (dist_a > .1 && dist_a < .4) {
2088 hPullsv1a_0p1_0p4->Fill(res1 / sqrt(rh_sigm_seg[6][it1]));
2089 }
2090 if (dist_a > .4) {
2091 hPullsv1a_0p4_0p5->Fill(res1 / sqrt(rh_sigm_seg[6][it1]));
2092 }
2093 if (dist_b < .1) {
2094 hPullsv1b_0p1->Fill(res2 / sqrt(rh_sigm_seg[7][it1]));
2095 }
2096 if (dist_b > .1 && dist_b < .4) {
2097 hPullsv1b_0p1_0p4->Fill(res2 / sqrt(rh_sigm_seg[7][it1]));
2098 }
2099 if (dist_b > .4) {
2100 hPullsv1b_0p4_0p5->Fill(res2 / sqrt(rh_sigm_seg[7][it1]));
2101 }
2102
2103 } else {
2104 double res1 = CalculateResidual(0, it1, rh_seg, par_ab);
2105 double res2 = CalculateResidual(1, it1, rh_seg, par_ab);
2106 hResidx2a->Fill(res1);
2107 hResidx2b->Fill(res2);
2108
2109 double int1;
2110 // double int2;
2111 double dist_a = 999;
2112 double dist_b = 999;
2113 if (fabs(modf(rh_seg[0][it1], &int1)) < .5) {
2114 dist_a = fabs(modf(rh_seg[0][it1], &int1));
2115 } else {
2116 dist_a = 1 - fabs(modf(rh_seg[0][it1], &int1));
2117 }
2118 if (fabs(modf(rh_seg[1][it1], &int1)) < .5) {
2119 dist_b = 0.5 - fabs(modf(rh_seg[1][it1], &int1));
2120 } else {
2121 dist_b = fabs(modf(rh_seg[1][it1], &int1)) - .5;
2122 }
2123 if (dist_a < .1)
2124 hResidx2a_0p1->Fill(res1);
2125 if (dist_a > .1 && dist_a < .4)
2126 hResidx2a_0p1_0p4->Fill(res1);
2127 if (dist_a > .4)
2128 hResidx2a_0p4_0p5->Fill(res1);
2129 if (dist_b < .1)
2130 hResidx2b_0p1->Fill(res2);
2131 if (dist_b > .1 && dist_b < .4)
2132 hResidx2b_0p1_0p4->Fill(res2);
2133 if (dist_b > .4)
2134 hResidx2b_0p4_0p5->Fill(res2);
2135 if (dist_a < .1) {
2136 hPullsx2a_0p1->Fill(res1 / sqrt(rh_sigm_seg[0][it1]));
2137 }
2138 if (dist_a > .1 && dist_a < .4) {
2139 hPullsx2a_0p1_0p4->Fill(res1 / sqrt(rh_sigm_seg[0][it1]));
2140 }
2141 if (dist_a > .4) {
2142 hPullsx2a_0p4_0p5->Fill(res1 / sqrt(rh_sigm_seg[0][it1]));
2143 }
2144 if (dist_b < .1) {
2145 hPullsx2b_0p1->Fill(res2 / sqrt(rh_sigm_seg[1][it1]));
2146 }
2147 if (dist_b > .1 && dist_b < .4) {
2148 hPullsx2b_0p1_0p4->Fill(res2 / sqrt(rh_sigm_seg[1][it1]));
2149 }
2150 if (dist_b > .4) {
2151 hPullsx2b_0p4_0p5->Fill(res2 / sqrt(rh_sigm_seg[1][it1]));
2152 }
2153 hx2Da_Db_best->Fill(dist_a + dist_b - 0.5);
2154
2155 // y
2156 dist_a = 999;
2157 dist_b = 999;
2158 res1 = CalculateResidual(2, it1, rh_seg, par_ab);
2159 res2 = CalculateResidual(3, it1, rh_seg, par_ab);
2160 hResidy2a->Fill(res1);
2161 hResidy2b->Fill(res2);
2162
2163 if (fabs(modf(rh_seg[2][it1], &int1)) < .5) {
2164 dist_a = fabs(modf(rh_seg[2][it1], &int1));
2165 } else {
2166 dist_a = 1 - fabs(modf(rh_seg[2][it1], &int1));
2167 }
2168 if (fabs(modf(rh_seg[3][it1], &int1)) < .5) {
2169 dist_b = 0.5 - fabs(modf(rh_seg[3][it1], &int1));
2170 } else {
2171 dist_b = fabs(modf(rh_seg[3][it1], &int1)) - .5;
2172 }
2173 if (dist_a < .1)
2174 hResidy2a_0p1->Fill(res1);
2175 if (dist_a > .1 && dist_a < .4)
2176 hResidy2a_0p1_0p4->Fill(res1);
2177 if (dist_a > .4)
2178 hResidy2a_0p4_0p5->Fill(res1);
2179 if (dist_b < .1)
2180 hResidy2b_0p1->Fill(res2);
2181 if (dist_b > .1 && dist_b < .4)
2182 hResidy2b_0p1_0p4->Fill(res2);
2183 if (dist_b > .4)
2184 hResidy2b_0p4_0p5->Fill(res2);
2185 if (dist_a < .1) {
2186 hPullsy2a_0p1->Fill(res1 / sqrt(rh_sigm_seg[2][it1]));
2187 }
2188 if (dist_a > .1 && dist_a < .4) {
2189 hPullsy2a_0p1_0p4->Fill(res1 / sqrt(rh_sigm_seg[2][it1]));
2190 }
2191 if (dist_a > .4) {
2192 hPullsy2a_0p4_0p5->Fill(res1 / sqrt(rh_sigm_seg[2][it1]));
2193 }
2194 if (dist_b < .1) {
2195 hPullsy2b_0p1->Fill(res2 / sqrt(rh_sigm_seg[3][it1]));
2196 }
2197 if (dist_b > .1 && dist_b < .4) {
2198 hPullsy2b_0p1_0p4->Fill(res2 / sqrt(rh_sigm_seg[3][it1]));
2199 }
2200 if (dist_b > .4) {
2201 hPullsy2b_0p4_0p5->Fill(res2 / sqrt(rh_sigm_seg[3][it1]));
2202 }
2203 hy2Da_Db_best->Fill(dist_a + dist_b - 0.5);
2204
2205 // u
2206 dist_a = 999;
2207 dist_b = 999;
2208 res1 = CalculateResidual(4, it1, rh_seg, par_ab);
2209 res2 = CalculateResidual(5, it1, rh_seg, par_ab);
2210 hResidu2a->Fill(res1);
2211 hResidu2b->Fill(res2);
2212
2213 if (fabs(modf(rh_seg[4][it1], &int1)) < .5) {
2214 dist_a = fabs(modf(rh_seg[4][it1], &int1));
2215 } else {
2216 dist_a = 1 - fabs(modf(rh_seg[4][it1], &int1));
2217 }
2218 if (fabs(modf(rh_seg[5][it1], &int1)) < .5) {
2219 dist_b = 0.5 - fabs(modf(rh_seg[5][it1], &int1));
2220 } else {
2221 dist_b = fabs(modf(rh_seg[5][it1], &int1)) - .5;
2222 }
2223 if (dist_a < .1)
2224 hResidu2a_0p1->Fill(res1);
2225 if (dist_a > .1 && dist_a < .4)
2226 hResidu2a_0p1_0p4->Fill(res1);
2227 if (dist_a > .4)
2228 hResidu2a_0p4_0p5->Fill(res1);
2229 if (dist_b < .1)
2230 hResidu2b_0p1->Fill(res2);
2231 if (dist_b > .1 && dist_b < .4)
2232 hResidu2b_0p1_0p4->Fill(res2);
2233 if (dist_b > .4)
2234 hResidu2b_0p4_0p5->Fill(res2);
2235 if (dist_a < .1) {
2236 hPullsu2a_0p1->Fill(res1 / sqrt(rh_sigm_seg[4][it1]));
2237 }
2238 if (dist_a > .1 && dist_a < .4) {
2239 hPullsu2a_0p1_0p4->Fill(res1 / sqrt(rh_sigm_seg[4][it1]));
2240 }
2241 if (dist_a > .4) {
2242 hPullsu2a_0p4_0p5->Fill(res1 / sqrt(rh_sigm_seg[4][it1]));
2243 }
2244 if (dist_b < .1) {
2245 hPullsu2b_0p1->Fill(res2 / sqrt(rh_sigm_seg[5][it1]));
2246 }
2247 if (dist_b > .1 && dist_b < .4) {
2248 hPullsu2b_0p1_0p4->Fill(res2 / sqrt(rh_sigm_seg[5][it1]));
2249 }
2250 if (dist_b > .4) {
2251 hPullsu2b_0p4_0p5->Fill(res2 / sqrt(rh_sigm_seg[5][it1]));
2252 }
2253 hu2Da_Db_best->Fill(dist_a + dist_b - 0.5);
2254
2255 // v
2256 dist_a = 999;
2257 dist_b = 999;
2258 res1 = CalculateResidual(6, it1, rh_seg, par_ab);
2259 res2 = CalculateResidual(7, it1, rh_seg, par_ab);
2260 hResidv2a->Fill(res1);
2261 hResidv2b->Fill(res2);
2262
2263 if (fabs(modf(rh_seg[6][it1], &int1)) < .5) {
2264 dist_a = fabs(modf(rh_seg[6][it1], &int1));
2265 } else {
2266 dist_a = 1 - fabs(modf(rh_seg[6][it1], &int1));
2267 }
2268 if (fabs(modf(rh_seg[7][it1], &int1)) < .5) {
2269 dist_b = 0.5 - fabs(modf(rh_seg[7][it1], &int1));
2270 } else {
2271 dist_b = fabs(modf(rh_seg[7][it1], &int1)) - .5;
2272 }
2273 if (dist_a < .1)
2274 hResidv2a_0p1->Fill(res1);
2275 if (dist_a > .1 && dist_a < .4)
2276 hResidv2a_0p1_0p4->Fill(res1);
2277 if (dist_a > .4)
2278 hResidv2a_0p4_0p5->Fill(res1);
2279 if (dist_b < .1)
2280 hResidv2b_0p1->Fill(res2);
2281 if (dist_b > .1 && dist_b < .4)
2282 hResidv2b_0p1_0p4->Fill(res2);
2283 if (dist_b > .4)
2284 hResidv2b_0p4_0p5->Fill(res2);
2285 if (dist_a < .1) {
2286 hPullsv2a_0p1->Fill(res1 / sqrt(rh_sigm_seg[6][it1]));
2287 }
2288 if (dist_a > .1 && dist_a < .4) {
2289 hPullsv2a_0p1_0p4->Fill(res1 / sqrt(rh_sigm_seg[6][it1]));
2290 }
2291 if (dist_a > .4) {
2292 hPullsv2a_0p4_0p5->Fill(res1 / sqrt(rh_sigm_seg[6][it1]));
2293 }
2294 if (dist_b < .1) {
2295 hPullsv2b_0p1->Fill(res2 / sqrt(rh_sigm_seg[7][it1]));
2296 }
2297 if (dist_b > .1 && dist_b < .4) {
2298 hPullsv2b_0p1_0p4->Fill(res2 / sqrt(rh_sigm_seg[7][it1]));
2299 }
2300 if (dist_b > .4) {
2301 hPullsv2b_0p4_0p5->Fill(res2 / sqrt(rh_sigm_seg[7][it1]));
2302 }
2303 hv2Da_Db_best->Fill(dist_a + dist_b - 0.5);
2304
2305 } // else
2306
2307 // end fill residuals
2308 }
2309 }
2310}
2311
2312void BmnDchTrackFinder::FindSegmentTrackMCId(Int_t dchID,
2313 Int_t** rhId_seg,
2314 Double_t* chi2,
2315 Int_t* SegMC_Id,
2316 Int_t* SegMC_IdCount)
2317{
2318 Int_t nDC_segments = (dchID == 1) ? nSegments[0] : nSegments[1];
2319 for (Int_t it1 = 0; it1 < nDC_segments; it1++) {
2320 if (chi2[it1] > 990.)
2321 continue;
2322 // calculate segment param errors
2323
2324 // Int_t mcId = -999;
2325 // int prev = -9;
2326 int count = 0;
2327 Int_t maxHits = 0;
2328 Int_t bestId = -99;
2329
2330 for (int i = 0; i < 8; i++) {
2331 if (rhId_seg[i][it1] == -999 || rhId_seg[i][it1] == bestId)
2332 continue;
2333 count = 1;
2334 for (int j = i + 1; j < 8; j++) {
2335 if (rhId_seg[j][it1] == -999)
2336 continue;
2337 if (rhId_seg[i][it1] == rhId_seg[j][it1])
2338 count++;
2339 }
2340 if (count > maxHits) {
2341 maxHits = count;
2342 bestId = rhId_seg[i][it1];
2343 }
2344 }
2345
2346 if (maxHits > 3) {
2347 SegMC_Id[it1] = bestId;
2348 SegMC_IdCount[it1] = maxHits;
2349 }
2350 }
2351}
2352
2353void BmnDchTrackFinder::FillSegmentParametersSigmas(Int_t dchID,
2354 Double_t** rh_seg,
2355 Double_t** rh_sigm_seg,
2356 Double_t* chi2,
2357 Double_t** sigm_seg_par)
2358{
2359 Int_t nDC_segments = (dchID == 1) ? nSegments[0] : nSegments[1];
2360 for (Int_t it1 = 0; it1 < nDC_segments; it1++) {
2361 if (chi2[it1] > 990.)
2362 continue;
2363 // calculate segment param errors
2364 Double_t XSum = 0.;
2365 Double_t XSumZ = 0.;
2366 Double_t XSumZZ = 0.;
2367 Double_t YSum = 0.;
2368 Double_t YSumZ = 0.;
2369 Double_t YSumZZ = 0.;
2370 for (int i = 0; i < 8; i++) {
2371 if (rh_seg[i][it1] == -999 || i == 2 || i == 3)
2372 continue;
2373 if (i == 0 || i == 1) {
2374 XSum += 1. / rh_sigm_seg[i][it1];
2375 XSumZ += z_loc[i] / rh_sigm_seg[i][it1];
2376 XSumZZ += z_loc[i] * z_loc[i] / rh_sigm_seg[i][it1];
2377 } else {
2378 XSum += 1. / (2. * rh_sigm_seg[i][it1]);
2379 XSumZ += z_loc[i] / (2. * rh_sigm_seg[i][it1]);
2380 XSumZZ += z_loc[i] * z_loc[i] / (2. * rh_sigm_seg[i][it1]);
2381 }
2382 }
2383 for (int i = 0; i < 8; i++) {
2384 if (rh_seg[i][it1] == -999 || i == 0 || i == 1)
2385 continue;
2386 if (i == 2 || i == 3) {
2387 YSum += 1. / rh_sigm_seg[i][it1];
2388 YSumZ += z_loc[i] / rh_sigm_seg[i][it1];
2389 YSumZZ += z_loc[i] * z_loc[i] / rh_sigm_seg[i][it1];
2390 } else {
2391 YSum += 1. / (2. * rh_sigm_seg[i][it1]);
2392 YSumZ += z_loc[i] / (2. * rh_sigm_seg[i][it1]);
2393 YSumZZ += z_loc[i] * z_loc[i] / (2. * rh_sigm_seg[i][it1]);
2394 }
2395 }
2396 Double_t Xdelta = XSum * XSumZZ - XSumZ * XSumZ;
2397 Double_t Ydelta = YSum * YSumZZ - YSumZ * YSumZ;
2398 sigm_seg_par[0][it1] = sqrt(XSum / Xdelta); // ax
2399 sigm_seg_par[1][it1] = sqrt(XSumZZ / Xdelta); // bx
2400 sigm_seg_par[2][it1] = sqrt(YSum / Ydelta); // ay
2401 sigm_seg_par[3][it1] = sqrt(YSumZZ / Ydelta); // by
2402 }
2403}
2404
2405void BmnDchTrackFinder::CreateDchTrack(Int_t dchID,
2406 Double_t* chi2Arr,
2407 Double_t** parArr,
2408 Int_t* sizeArr,
2409 Double_t** rh_seg,
2410 Double_t** rh_seg_time,
2411 Double_t** sigm_seg_par,
2412 Int_t* SegMC_Id,
2413 Int_t* SegMC_IdCount)
2414{
2415 Int_t nDC_segments = (dchID == 1) ? nSegments[0] : nSegments[1];
2416 Int_t good_dc_segs = 0;
2417 Int_t best_seg1 = -1;
2418 Int_t best_chi1 = 1000.;
2419 Double_t sqrt_2 = sqrt(2.);
2420 Double_t isqrt_2 = 1. / sqrt_2;
2421 for (Int_t iSegment = 0; iSegment < nDC_segments; iSegment++) {
2422 if (chi2Arr[iSegment] > 990.)
2423 continue;
2424 good_dc_segs++;
2425
2426 if (sizeArr[iSegment] == 8) {
2427 if (chi2Arr[iSegment] < best_chi1) {
2428 best_chi1 = chi2Arr[iSegment];
2429 best_seg1 = iSegment;
2430 }
2431 }
2432 FairTrackParam trackParam;
2433 Double_t z0 = (dchID == 1) ? Z_dch1 : Z_dch2;
2434 Double_t x0 = parArr[1][iSegment];
2435 Double_t y0 = parArr[3][iSegment];
2436 Double_t x_align = (dchID == 1) ? 0.0 : 0.0; //-10.89 : +4.93;
2437 Double_t y_align = (dchID == 1) ? 0.0 : 0.0; //-2.95 : +9.12;
2438 Double_t tx_align = (dchID == 1) ? 0.0 : 0.0; //+0.073 + 0.049 : +0.070 + 0.052;
2439 Double_t ty_align = (dchID == 1) ? 0.0 : 0.0; //+0.054 : +0.062 - 0.047;
2440 trackParam.SetPosition(TVector3(-x0 + x_align, y0 + y_align, z0));
2441 trackParam.SetTx(-parArr[0][iSegment] + tx_align);
2442 trackParam.SetTy(parArr[2][iSegment] + ty_align);
2443 trackParam.SetCovariance(0, 0, sigm_seg_par[1][iSegment] * sigm_seg_par[1][iSegment]); // bx^2
2444 trackParam.SetCovariance(1, 1, sigm_seg_par[3][iSegment] * sigm_seg_par[3][iSegment]); // by^2
2445 trackParam.SetCovariance(2, 2, sigm_seg_par[0][iSegment] * sigm_seg_par[0][iSegment]); // tx^2
2446 trackParam.SetCovariance(3, 3, sigm_seg_par[2][iSegment] * sigm_seg_par[2][iSegment]); // ty^2
2447 BmnDchTrack* track = new ((*fDchTracks)[fDchTracks->GetEntriesFast()]) BmnDchTrack();
2448 track->SetChi2(chi2Arr[iSegment]);
2449 if (dchID == 1) {
2450 hChi2ndf_dc1->Fill(chi2Arr[iSegment]);
2451 } else {
2452 hChi2ndf_dc2->Fill(chi2Arr[iSegment]);
2453 }
2454 track->SetTrackId(SegMC_Id[iSegment]);
2455 track->SetLength(SegMC_IdCount[iSegment]);
2456 track->SetFlag(iSegment);
2457 track->SetNHits(sizeArr[iSegment]);
2458 track->SetParamFirst(trackParam);
2459 if (dchID == 1) {
2460 haX1->Fill(-parArr[0][iSegment]);
2461 haY1->Fill(parArr[2][iSegment]);
2462 hx1->Fill(parArr[1][iSegment]);
2463 hy1->Fill(parArr[3][iSegment]);
2464 hlocXY1->Fill(rh_seg[0][iSegment], rh_seg[2][iSegment]);
2465 hXvsAx->Fill(parArr[1][iSegment], -parArr[0][iSegment]);
2466 hSegSize_dc1->Fill(sizeArr[iSegment]);
2467 for (int k = 0; k < 8; k++) {
2468 if (sizeArr[iSegment] == 6 && rh_seg[k][iSegment] == -999)
2469 hDenom1->Fill(k, 1);
2470 if (sizeArr[iSegment] != 6) {
2471 if (rh_seg[k][iSegment] == -999) {
2472 hDenom1->Fill(k, 1);
2473 } else {
2474 hNomin1->Fill(k, 1);
2475 hDenom1->Fill(k, 1);
2476 }
2477 }
2478 } // for
2479
2480 } else {
2481 haX2->Fill(-parArr[0][iSegment]);
2482 haY2->Fill(parArr[2][iSegment]);
2483 hx2->Fill(parArr[1][iSegment]);
2484 hy2->Fill(parArr[3][iSegment]);
2485 hlocXY2->Fill(rh_seg[0][iSegment], rh_seg[2][iSegment]);
2486 hSegSize_dc2->Fill(sizeArr[iSegment]);
2487 for (int k = 0; k < 8; k++) {
2488 if (sizeArr[iSegment] == 6 && rh_seg[k][iSegment] == -999)
2489 hDenom2->Fill(k, 1);
2490 if (sizeArr[iSegment] != 6) {
2491 if (rh_seg[k][iSegment] == -999) {
2492 hDenom2->Fill(k, 1);
2493 } else {
2494 hNomin2->Fill(k, 1);
2495 hDenom2->Fill(k, 1);
2496 }
2497 }
2498 } // for
2499 }
2500 } // for
2501 if (dchID == 1 && best_seg1 != -1) {
2502 hSlot_1xa_time8p->Fill(rh_seg_time[0][best_seg1]);
2503 hSlot_1xb_time8p->Fill(rh_seg_time[1][best_seg1]);
2504 hSlot_1ya_time8p->Fill(rh_seg_time[2][best_seg1]);
2505 hSlot_1yb_time8p->Fill(rh_seg_time[3][best_seg1]);
2506 hSlot_1ua_time8p->Fill(rh_seg_time[4][best_seg1]);
2507 hSlot_1ub_time8p->Fill(rh_seg_time[5][best_seg1]);
2508 hSlot_1va_time8p->Fill(rh_seg_time[6][best_seg1]);
2509 hSlot_1vb_time8p->Fill(rh_seg_time[7][best_seg1]);
2510
2511 double vAtu = 4.8 * (-.074) + (0.5 * (rh_seg[6][best_seg1] + rh_seg[7][best_seg1]));
2512 double xAtu = (-9.6) * (-.107) + (0.5 * (rh_seg[0][best_seg1] + rh_seg[1][best_seg1]));
2513 double uEst = vAtu - sqrt(2.) * xAtu;
2514 double diffU = 0.5 * (rh_seg[4][best_seg1] + rh_seg[5][best_seg1]) - uEst;
2515 hUvsXV1->Fill(diffU);
2516
2517 double vAty = 9.6 * (-.074) + (0.5 * (rh_seg[6][best_seg1] + rh_seg[7][best_seg1])); //-.0802
2518 double xAty = (-4.8) * (-.107) + (0.5 * (rh_seg[0][best_seg1] + rh_seg[1][best_seg1])); //-.1135
2519 double yEst = sqrt(2.) * vAty - xAty;
2520 double diffY = 0.5 * (rh_seg[2][best_seg1] + rh_seg[3][best_seg1]) - yEst;
2521 hYvsXV1->Fill(diffY);
2522
2523 hAx12_loc->Fill((rh_seg[0][best_seg1] - rh_seg[1][best_seg1]) / (z_loc[0] - z_loc[1]));
2524 hAy12_loc->Fill((rh_seg[2][best_seg1] - rh_seg[3][best_seg1]) / (z_loc[2] - z_loc[3]));
2525 hAu12_loc->Fill((rh_seg[4][best_seg1] - rh_seg[5][best_seg1]) / (z_loc[4] - z_loc[5]));
2526 hAv12_loc->Fill((rh_seg[6][best_seg1] - rh_seg[7][best_seg1]) / (z_loc[6] - z_loc[7]));
2527
2528 hAx1_loc->Fill(parArr[0][best_seg1]);
2529 hAy1_loc->Fill(parArr[2][best_seg1]); //
2530 hAu1_loc->Fill(isqrt_2 * (parArr[2][best_seg1] - parArr[0][best_seg1]));
2531 hAv1_loc->Fill(isqrt_2 * (parArr[2][best_seg1] + parArr[0][best_seg1]));
2532 }
2533
2534 if (dchID == 2 && best_seg1 != -1) {
2535 hSlot_2xa_time8p->Fill(rh_seg_time[0][best_seg1]);
2536 hSlot_2xb_time8p->Fill(rh_seg_time[1][best_seg1]);
2537 hSlot_2ya_time8p->Fill(rh_seg_time[2][best_seg1]);
2538 hSlot_2yb_time8p->Fill(rh_seg_time[3][best_seg1]);
2539 hSlot_2ua_time8p->Fill(rh_seg_time[4][best_seg1]);
2540 hSlot_2ub_time8p->Fill(rh_seg_time[5][best_seg1]);
2541 hSlot_2va_time8p->Fill(rh_seg_time[6][best_seg1]);
2542 hSlot_2vb_time8p->Fill(rh_seg_time[7][best_seg1]);
2543 double vAtu = 4.8 * (-.075) + (0.5 * (rh_seg[6][best_seg1] + rh_seg[7][best_seg1])); //-.0632
2544 double xAtu = (-9.6) * (-.103) + (0.5 * (rh_seg[0][best_seg1] + rh_seg[1][best_seg1])); //-.0781
2545 double uEst = vAtu - sqrt(2.) * xAtu;
2546 double diffU = 0.5 * (rh_seg[4][best_seg1] + rh_seg[5][best_seg1]) - uEst;
2547 hUvsXV2->Fill(diffU);
2548
2549 double vAty = 9.6 * (-.075) + (0.5 * (rh_seg[6][best_seg1] + rh_seg[7][best_seg1]));
2550 double xAty = (-4.8) * (-.103) + (0.5 * (rh_seg[0][best_seg1] + rh_seg[1][best_seg1]));
2551 double yEst = sqrt(2.) * vAty - xAty;
2552 double diffY = 0.5 * (rh_seg[2][best_seg1] + rh_seg[3][best_seg1]) - yEst;
2553 hYvsXV2->Fill(diffY);
2554
2555 hAx122loc->Fill((rh_seg[0][best_seg1] - rh_seg[1][best_seg1]) / (z_loc[0] - z_loc[1]));
2556 hAy122loc->Fill((rh_seg[2][best_seg1] - rh_seg[3][best_seg1]) / (z_loc[2] - z_loc[3]));
2557 hAu122loc->Fill((rh_seg[4][best_seg1] - rh_seg[5][best_seg1]) / (z_loc[4] - z_loc[5]));
2558 hAv122loc->Fill((rh_seg[6][best_seg1] - rh_seg[7][best_seg1]) / (z_loc[6] - z_loc[7]));
2559
2560 hAx2_loc->Fill(parArr[0][best_seg1]);
2561 hAy2_loc->Fill(parArr[2][best_seg1]); //
2562 hAu2_loc->Fill(isqrt_2 * (parArr[2][best_seg1] - parArr[0][best_seg1]));
2563 hAv2_loc->Fill(isqrt_2 * (parArr[2][best_seg1] + parArr[0][best_seg1]));
2564 }
2565
2566 if (dchID == 1 && layers_with > 5) {
2567 Ntrack1->Fill(good_dc_segs);
2568 }
2569 if (dchID == 2 && layers_with2 > 5) {
2570 Ntrack2->Fill(good_dc_segs);
2571 }
2572}
2573
2574void BmnDchTrackFinder::CreateDchTrack()
2575{
2576
2577 if (layers_with > 5 && layers_with2 > 5)
2578 NGtracks->Fill(nSegmentsToBeMatched + 1);
2579 for (Int_t iSeg = 0; iSeg < nSegmentsToBeMatched + 1; iSeg++) {
2580 FairTrackParam trackParam;
2581 Double_t z0 = Z_dch_mid;
2582 Double_t x0 = -x_mid[iSeg];
2583 Double_t y0 = y_mid[iSeg];
2584 trackParam.SetPosition(TVector3(x0, y0, z0)); // Go to right reference frame
2585 trackParam.SetTx(aX[iSeg]); // Go to right reference frame
2586 trackParam.SetTy(aY[iSeg]);
2587 hP->Fill(0.3 * 2.84 / (atan(aX[iSeg])));
2588 Int_t iseg1 = int((pairID[iSeg] - 1000) / 1000);
2589 Int_t iseg2 = (pairID[iSeg] % 1000) - 1;
2590 Double_t sigma2_dx = ((params_sigmas[1][1][iseg2] * params_sigmas[1][1][iseg2])
2591 + (params_sigmas[0][1][iseg1] * params_sigmas[0][1][iseg1]))
2592 * 99.75 * 99.75
2593 + ((params_sigmas[0][0][iseg1]) * (params_sigmas[0][0][iseg1]))
2594 + ((params_sigmas[1][0][iseg2]) * (params_sigmas[1][0][iseg2]));
2595
2596 Double_t sigma2_dy = ((params_sigmas[1][3][iseg2] * params_sigmas[1][3][iseg2])
2597 + (params_sigmas[0][3][iseg1] * params_sigmas[0][3][iseg1]))
2598 * 99.75 * 99.75
2599 + ((params_sigmas[0][2][iseg1]) * (params_sigmas[0][2][iseg1]))
2600 + ((params_sigmas[1][2][iseg2]) * (params_sigmas[1][2][iseg2]));
2601
2602 Double_t sigma2_dax = 9.
2603 * ((params_sigmas[1][0][iseg2] * params_sigmas[1][0][iseg2])
2604 + (params_sigmas[0][0][iseg1] * params_sigmas[0][0][iseg1]));
2605
2606 Double_t sigma2_day = 9.
2607 * ((params_sigmas[1][2][iseg2] * params_sigmas[1][2][iseg2])
2608 + (params_sigmas[0][2][iseg1] * params_sigmas[0][2][iseg1]));
2609 trackParam.SetCovariance(0, 0, sigma2_dx);
2610 trackParam.SetCovariance(1, 1, sigma2_dy);
2611 trackParam.SetCovariance(2, 2, sigma2_dax);
2612 trackParam.SetCovariance(3, 3, sigma2_day);
2613 BmnDchTrack* track = new ((*fDchTracks)[fDchTracks->GetEntriesFast()]) BmnDchTrack();
2614 track->SetChi2(Chi2_match[iSeg]);
2615 track->SetNHits(nhits[iSeg]);
2616 track->SetFlag(pairID[iSeg]);
2617 track->SetParamFirst(trackParam);
2618
2619 if (SegMCId[0][iseg1] != SegMCId[1][iseg2]) {
2620 if (SegMCIdCount[0][iseg1] != SegMCIdCount[1][iseg2]) {
2621 if (SegMCIdCount[0][iseg1] > SegMCIdCount[1][iseg2]) {
2622 track->SetTrackId(SegMCId[0][iseg1]);
2623 track->SetLength(SegMCIdCount[0][iseg1]);
2624 if (layers_with > 5 && layers_with2 > 5) {
2625 // if(SegMCId[0][iseg1]==0)NGtracks_beam->Fill(1);
2626 // if(SegMCId[0][iseg1]>0)NGtracks_bkg->Fill(1);
2627 }
2628 }
2629 if (SegMCIdCount[0][iseg1] < SegMCIdCount[1][iseg2]) {
2630 track->SetTrackId(SegMCId[1][iseg2]);
2631 track->SetLength(SegMCIdCount[1][iseg2]);
2632 if (layers_with > 5 && layers_with2 > 5) {
2633 // if(SegMCId[1][iseg2]==0)NGtracks_beam->Fill(1);
2634 // if(SegMCId[1][iseg2]>0)NGtracks_bkg->Fill(1);
2635 }
2636 }
2637 } else {
2638 if (Chi2[0][iseg1] < Chi2[1][iseg2]) {
2639 track->SetTrackId(SegMCId[0][iseg1]);
2640 track->SetLength(SegMCIdCount[0][iseg1]);
2641 if (layers_with > 5 && layers_with2 > 5) {
2642 // if(SegMCId[0][iseg1]==0)NGtracks_beam->Fill(1);
2643 // if(SegMCId[0][iseg1]>0)NGtracks_bkg->Fill(1);
2644 }
2645 } else {
2646
2647 track->SetTrackId(SegMCId[1][iseg2]);
2648 track->SetLength(SegMCIdCount[1][iseg2]);
2649 if (layers_with > 5 && layers_with2 > 5) {
2650 // if(SegMCId[1][iseg2]==0)NGtracks_beam->Fill(1);
2651 // if(SegMCId[1][iseg2]>0)NGtracks_bkg->Fill(1);
2652 }
2653 }
2654 }
2655 } else {
2656 if (layers_with > 5 && layers_with2 > 5) {
2657 // if(SegMCId[0][iseg1]==0)NGtracks_beam->Fill(1);
2658 // if(SegMCId[0][iseg1]>0)NGtracks_bkg->Fill(1);
2659 }
2660
2661 track->SetTrackId(SegMCId[0][iseg1]);
2662 track->SetLength(SegMCIdCount[0][iseg1] + SegMCIdCount[1][iseg2]);
2663 }
2664 }
2665}
2666Double_t BmnDchTrackFinder::CalculateResidualMatch(Int_t dchID,
2667 Int_t i,
2668 Int_t j,
2669 Double_t*** rh_seg,
2670 Double_t*** par_ab)
2671{
2672 Double_t sqrt_2 = sqrt(2.);
2673 Double_t isqrt_2 = 1 / sqrt_2;
2674
2675 return (i < 2) ? rh_seg[dchID][i][j] - z_loc[i] * par_ab[dchID][0][j] - par_ab[dchID][1][j]
2676 : (i >= 2 && i < 4) ? rh_seg[dchID][i][j] - z_loc[i] * par_ab[dchID][2][j] - par_ab[dchID][3][j]
2677 : (i >= 4 && i < 6) ? rh_seg[dchID][i][j] - isqrt_2 * z_loc[i] * (par_ab[dchID][2][j] - par_ab[dchID][0][j])
2678 - isqrt_2 * (par_ab[dchID][3][j] - par_ab[dchID][1][j])
2679 : rh_seg[dchID][i][j] - isqrt_2 * z_loc[i] * (par_ab[dchID][2][j] + par_ab[dchID][0][j])
2680 - isqrt_2 * (par_ab[dchID][3][j] + par_ab[dchID][1][j]);
2681}
2682
2683Double_t BmnDchTrackFinder::CalculateResidual(Int_t i, Int_t j, Double_t** rh_seg, Double_t** par_ab)
2684{
2685 Double_t sqrt_2 = sqrt(2.);
2686 Double_t isqrt_2 = 1 / sqrt_2;
2687
2688 return (i < 2) ? rh_seg[i][j] - z_loc[i] * par_ab[0][j] - par_ab[1][j]
2689 : (i >= 2 && i < 4) ? rh_seg[i][j] - z_loc[i] * par_ab[2][j] - par_ab[3][j]
2690 : (i >= 4 && i < 6) ? rh_seg[i][j] - isqrt_2 * z_loc[i] * (par_ab[2][j] - par_ab[0][j])
2691 - isqrt_2 * (par_ab[3][j] - par_ab[1][j])
2692 : rh_seg[i][j] - isqrt_2 * z_loc[i] * (par_ab[2][j] + par_ab[0][j])
2693 - isqrt_2 * (par_ab[3][j] + par_ab[1][j]);
2694}
2695
2696InitStatus BmnDchTrackFinder::Init()
2697{
2698 if (!expData) {
2699 // cout << "BmnDchTrackFinder::Init(): simulation data is reconstructed " << endl;
2700 SetActive(kTRUE);
2701 }
2702
2703 if (fVerbose)
2704 cout << "BmnDchTrackFinder::Init()" << endl;
2705 FairRootManager* ioman = FairRootManager::Instance();
2706
2707 if (expData) {
2708 fBmnDchDigitsArray = (TClonesArray*)ioman->GetObject(InputDigitsBranchName);
2709 if (!fBmnDchDigitsArray) {
2710 cout << "BmnDchTrackFinder::Init(): branch " << InputDigitsBranchName
2711 << " not found! Task will be deactivated" << endl;
2712 SetActive(kFALSE);
2713 return kERROR;
2714 }
2715 } else {
2716 fBmnHitsArray = (TClonesArray*)ioman->GetObject(InputDigitsBranchName);
2717 if (!fBmnHitsArray) {
2718 cout << "BmnDchTrackFinder::Init(): branch " << InputDigitsBranchName
2719 << " not found! Task will be deactivated" << endl;
2720 SetActive(kFALSE);
2721 return kERROR;
2722 }
2723 }
2724
2725 initHists();
2726
2727 // Create and register track arrays
2728 fDchTracks = new TClonesArray(tracksDch.Data());
2729 ioman->Register(tracksDch.Data(), "DCH", fDchTracks, kTRUE);
2730
2731 if (expData) { // transfer function read for experimental data
2732 ifstream fin;
2733 TString dir = getenv("VMCWORKDIR");
2734 dir += "/input/";
2735 fin.open((TString(dir + fTransferFunctionName)).Data(), ios::in);
2736 // cout << "Transfer function for DCH: " << fTransferFunctionName << endl;
2737 for (Int_t fi = 0; fi < 16; fi++) {
2738 fin >> t_dc[0][fi] >> t_dc[1][fi] >> t_dc[2][fi] >> t_dc[3][fi] >> pol_par_dc[0][0][fi]
2739 >> pol_par_dc[0][1][fi] >> pol_par_dc[0][2][fi] >> pol_par_dc[0][3][fi] >> pol_par_dc[0][4][fi]
2740 >> pol_par_dc[1][0][fi] >> pol_par_dc[1][1][fi] >> pol_par_dc[1][2][fi] >> pol_par_dc[1][3][fi]
2741 >> pol_par_dc[1][4][fi] >> pol_par_dc[2][0][fi] >> pol_par_dc[2][1][fi] >> pol_par_dc[2][2][fi]
2742 >> pol_par_dc[2][3][fi] >> pol_par_dc[2][4][fi] >> scaling[fi];
2743 if (fVerbose) {
2744 cout << t_dc[0][fi] << t_dc[1][fi] << t_dc[2][fi] << t_dc[3][fi] << pol_par_dc[0][0][fi]
2745 << pol_par_dc[0][1][fi] << pol_par_dc[0][2][fi] << pol_par_dc[0][3][fi] << pol_par_dc[0][4][fi]
2746 << pol_par_dc[1][0][fi] << pol_par_dc[1][1][fi] << pol_par_dc[1][2][fi] << pol_par_dc[1][3][fi]
2747 << pol_par_dc[1][4][fi] << pol_par_dc[2][0][fi] << pol_par_dc[2][1][fi] << pol_par_dc[2][2][fi]
2748 << pol_par_dc[2][3][fi] << pol_par_dc[2][4][fi] << endl;
2749 }
2750 }
2751 fin.close();
2752 }
2753
2754 // z local xa->vb (cm)
2755 Double_t arr1[8] = {7.8, 6.6, 3.0, 1.8, -1.8, -3.0, -6.6, -7.8};
2756 for (Int_t iSize = 0; iSize < 8; iSize++)
2757 z_loc[iSize] = arr1[iSize];
2758
2759 // z global dc 1 & dc 2 (cm)
2760 Double_t arr2[16] = {-45.7, -46.9, -51.5, -52.7, -57.3, -58.5, -63.1, -64.3,
2761 64.3, 63.1, 58.5, 57.3, 52.7, 51.5, 46.9, 45.7};
2762 for (Int_t iSize = 0; iSize < 16; iSize++)
2763 z_glob[iSize] = arr2[iSize];
2764
2765 has7DC = new Bool_t[nChambers];
2766 x_mid = new Double_t[25 * N];
2767 y_mid = new Double_t[25 * N];
2768 pairID = new Int_t[25 * N];
2769 nhits = new Int_t[25 * N];
2770 aX = new Double_t[25 * N];
2771 aY = new Double_t[25 * N];
2772 imp = new Double_t[25 * N];
2773 leng = new Double_t[25 * N];
2774 Chi2_match = new Double_t[25 * N];
2775 v = new Double_t**[nChambers];
2776 u = new Double_t**[nChambers];
2777 y = new Double_t**[nChambers];
2778 x = new Double_t**[nChambers];
2779 vId = new Int_t**[nChambers];
2780 uId = new Int_t**[nChambers];
2781 yId = new Int_t**[nChambers];
2782 xId = new Int_t**[nChambers];
2783 v_time = new Double_t**[nChambers];
2784 u_time = new Double_t**[nChambers];
2785 y_time = new Double_t**[nChambers];
2786 x_time = new Double_t**[nChambers];
2787 v_Single = new Double_t**[nChambers];
2788 u_Single = new Double_t**[nChambers];
2789 y_Single = new Double_t**[nChambers];
2790 x_Single = new Double_t**[nChambers];
2791 v_SingleId = new Int_t**[nChambers];
2792 u_SingleId = new Int_t**[nChambers];
2793 y_SingleId = new Int_t**[nChambers];
2794 x_SingleId = new Int_t**[nChambers];
2795 sigm_v = new Double_t**[nChambers];
2796 sigm_u = new Double_t**[nChambers];
2797 sigm_y = new Double_t**[nChambers];
2798 sigm_x = new Double_t**[nChambers];
2799 Sigm_v_single = new Double_t**[nChambers];
2800 Sigm_u_single = new Double_t**[nChambers];
2801 Sigm_y_single = new Double_t**[nChambers];
2802 Sigm_x_single = new Double_t**[nChambers];
2803 segment_size = new Int_t*[nChambers];
2804 Chi2 = new Double_t*[nChambers];
2805 SegMCId = new Int_t*[nChambers];
2806 SegMCIdCount = new Int_t*[nChambers];
2807 x_global = new Double_t*[nChambers];
2808 y_global = new Double_t*[nChambers];
2809 params = new Double_t**[nChambers];
2810 params_sigmas = new Double_t**[nChambers];
2811 rh_segment = new Double_t**[nChambers];
2812 rhId_segment = new Int_t**[nChambers];
2813 rh_segment_time = new Double_t**[nChambers];
2814 rh_sigm_segment = new Double_t**[nChambers];
2815 for (Int_t iChamber = 0; iChamber < nChambers; iChamber++) {
2816 v[iChamber] = new Double_t*[N];
2817 u[iChamber] = new Double_t*[N];
2818 y[iChamber] = new Double_t*[N];
2819 x[iChamber] = new Double_t*[N];
2820 vId[iChamber] = new Int_t*[N];
2821 uId[iChamber] = new Int_t*[N];
2822 yId[iChamber] = new Int_t*[N];
2823 xId[iChamber] = new Int_t*[N];
2824 v_time[iChamber] = new Double_t*[N];
2825 u_time[iChamber] = new Double_t*[N];
2826 y_time[iChamber] = new Double_t*[N];
2827 x_time[iChamber] = new Double_t*[N];
2828 v_Single[iChamber] = new Double_t*[N];
2829 u_Single[iChamber] = new Double_t*[N];
2830 y_Single[iChamber] = new Double_t*[N];
2831 x_Single[iChamber] = new Double_t*[N];
2832 v_SingleId[iChamber] = new Int_t*[N];
2833 u_SingleId[iChamber] = new Int_t*[N];
2834 y_SingleId[iChamber] = new Int_t*[N];
2835 x_SingleId[iChamber] = new Int_t*[N];
2836 sigm_v[iChamber] = new Double_t*[N];
2837 sigm_u[iChamber] = new Double_t*[N];
2838 sigm_y[iChamber] = new Double_t*[N];
2839 sigm_x[iChamber] = new Double_t*[N];
2840 Sigm_v_single[iChamber] = new Double_t*[N];
2841 Sigm_u_single[iChamber] = new Double_t*[N];
2842 Sigm_y_single[iChamber] = new Double_t*[N];
2843 Sigm_x_single[iChamber] = new Double_t*[N];
2844 segment_size[iChamber] = new Int_t[75 * N];
2845 Chi2[iChamber] = new Double_t[75 * N];
2846 SegMCId[iChamber] = new Int_t[75 * N];
2847 SegMCIdCount[iChamber] = new Int_t[75 * N];
2848 x_global[iChamber] = new Double_t[75 * N];
2849 y_global[iChamber] = new Double_t[75 * N];
2850 params[iChamber] = new Double_t*[4];
2851 params_sigmas[iChamber] = new Double_t*[4];
2852 rh_segment[iChamber] = new Double_t*[8];
2853 rhId_segment[iChamber] = new Int_t*[8];
2854 rh_segment_time[iChamber] = new Double_t*[8];
2855 rh_sigm_segment[iChamber] = new Double_t*[8];
2856 for (Int_t iDim = 0; iDim < N; iDim++) {
2857 v[iChamber][iDim] = new Double_t[75 * N];
2858 u[iChamber][iDim] = new Double_t[75 * N];
2859 y[iChamber][iDim] = new Double_t[75 * N];
2860 x[iChamber][iDim] = new Double_t[75 * N];
2861 vId[iChamber][iDim] = new Int_t[75 * N];
2862 uId[iChamber][iDim] = new Int_t[75 * N];
2863 yId[iChamber][iDim] = new Int_t[75 * N];
2864 xId[iChamber][iDim] = new Int_t[75 * N];
2865 v_time[iChamber][iDim] = new Double_t[75 * N];
2866 u_time[iChamber][iDim] = new Double_t[75 * N];
2867 y_time[iChamber][iDim] = new Double_t[75 * N];
2868 x_time[iChamber][iDim] = new Double_t[75 * N];
2869 v_Single[iChamber][iDim] = new Double_t[20 * N];
2870 u_Single[iChamber][iDim] = new Double_t[20 * N];
2871 y_Single[iChamber][iDim] = new Double_t[20 * N];
2872 x_Single[iChamber][iDim] = new Double_t[20 * N];
2873 v_SingleId[iChamber][iDim] = new Int_t[20 * N];
2874 u_SingleId[iChamber][iDim] = new Int_t[20 * N];
2875 y_SingleId[iChamber][iDim] = new Int_t[20 * N];
2876 x_SingleId[iChamber][iDim] = new Int_t[20 * N];
2877 sigm_v[iChamber][iDim] = new Double_t[75 * N];
2878 sigm_u[iChamber][iDim] = new Double_t[75 * N];
2879 sigm_y[iChamber][iDim] = new Double_t[75 * N];
2880 sigm_x[iChamber][iDim] = new Double_t[75 * N];
2881 Sigm_v_single[iChamber][iDim] = new Double_t[20 * N];
2882 Sigm_u_single[iChamber][iDim] = new Double_t[20 * N];
2883 Sigm_y_single[iChamber][iDim] = new Double_t[20 * N];
2884 Sigm_x_single[iChamber][iDim] = new Double_t[20 * N];
2885 }
2886 for (Int_t iDim = 0; iDim < 4; iDim++) {
2887 params[iChamber][iDim] = new Double_t[75 * N];
2888 params_sigmas[iChamber][iDim] = new Double_t[75 * N];
2889 }
2890 for (Int_t iDim = 0; iDim < 8; iDim++) {
2891 rh_segment[iChamber][iDim] = new Double_t[75 * N];
2892 rhId_segment[iChamber][iDim] = new Int_t[75 * N];
2893 rh_segment_time[iChamber][iDim] = new Double_t[75 * N];
2894 rh_sigm_segment[iChamber][iDim] = new Double_t[75 * N];
2895 }
2896 }
2897 pairs = new Int_t*[nChambers];
2898 for (Int_t iChamber = 0; iChamber < nChambers; iChamber++)
2899 pairs[iChamber] = new Int_t[nWires];
2900 singles = new Int_t**[nChambers];
2901 for (Int_t iChamber = 0; iChamber < nChambers; iChamber++) {
2902 singles[iChamber] = new Int_t*[nWires];
2903 for (Int_t iWire = 0; iWire < nWires; iWire++)
2904 singles[iChamber][iWire] = new Int_t[nLayers];
2905 }
2906 nSegments = new Int_t[nSegmentsMax];
2907
2908 isInitialized = true;
2909
2910 return kSUCCESS;
2911}
2912
2913void BmnDchTrackFinder::PrepareArraysToProcessEvent()
2914{
2915 nSegmentsToBeMatched = -1;
2916 fDchTracks->Clear();
2917
2918 // Array cleaning and initializing
2919 for (Int_t iChamber = 0; iChamber < nChambers; iChamber++) {
2920 for (Int_t iWire = 0; iWire < nWires; iWire++) {
2921 pairs[iChamber][iWire] = 0;
2922 for (Int_t iLayer = 0; iLayer < nLayers; iLayer++)
2923 singles[iChamber][iWire][iLayer] = 0;
2924 }
2925
2926 has7DC[iChamber] = kFALSE;
2927
2928 for (Int_t iDim1 = 0; iDim1 < 4; iDim1++)
2929 for (Int_t iDim2 = 0; iDim2 < 75 * N; iDim2++) {
2930 params[iChamber][iDim1][iDim2] = -999.;
2931 params_sigmas[iChamber][iDim1][iDim2] = -999.;
2932 }
2933 for (Int_t iDim1 = 0; iDim1 < 75 * N; iDim1++) {
2934 segment_size[iChamber][iDim1] = 0;
2935 Chi2[iChamber][iDim1] = 50.;
2936 SegMCId[iChamber][iDim1] = -999;
2937 SegMCIdCount[iChamber][iDim1] = 0;
2938 x_global[iChamber][iDim1] = -999.;
2939 y_global[iChamber][iDim1] = -999.;
2940 }
2941 for (Int_t iDim1 = 0; iDim1 < N; iDim1++) {
2942 for (Int_t iDim2 = 0; iDim2 < 75 * N; iDim2++) {
2943 v[iChamber][iDim1][iDim2] = -999.;
2944 u[iChamber][iDim1][iDim2] = -999.;
2945 y[iChamber][iDim1][iDim2] = -999.;
2946 x[iChamber][iDim1][iDim2] = -999.;
2947 vId[iChamber][iDim1][iDim2] = -999;
2948 uId[iChamber][iDim1][iDim2] = -999;
2949 yId[iChamber][iDim1][iDim2] = -999;
2950 xId[iChamber][iDim1][iDim2] = -999;
2951 v_time[iChamber][iDim1][iDim2] = -999.;
2952 u_time[iChamber][iDim1][iDim2] = -999.;
2953 y_time[iChamber][iDim1][iDim2] = -999.;
2954 x_time[iChamber][iDim1][iDim2] = -999.;
2955 sigm_v[iChamber][iDim1][iDim2] = 1.;
2956 sigm_u[iChamber][iDim1][iDim2] = 1.;
2957 sigm_y[iChamber][iDim1][iDim2] = 1.;
2958 sigm_x[iChamber][iDim1][iDim2] = 1.;
2959 }
2960 for (Int_t iDim3 = 0; iDim3 < 20 * N; iDim3++) {
2961 v_Single[iChamber][iDim1][iDim3] = -999.;
2962 u_Single[iChamber][iDim1][iDim3] = -999.;
2963 y_Single[iChamber][iDim1][iDim3] = -999.;
2964 x_Single[iChamber][iDim1][iDim3] = -999.;
2965 v_SingleId[iChamber][iDim1][iDim3] = -999;
2966 u_SingleId[iChamber][iDim1][iDim3] = -999;
2967 y_SingleId[iChamber][iDim1][iDim3] = -999;
2968 x_SingleId[iChamber][iDim1][iDim3] = -999;
2969 Sigm_v_single[iChamber][iDim1][iDim3] = 1.;
2970 Sigm_u_single[iChamber][iDim1][iDim3] = 1.;
2971 Sigm_y_single[iChamber][iDim1][iDim3] = 1.;
2972 Sigm_x_single[iChamber][iDim1][iDim3] = 1.;
2973 }
2974 }
2975 for (Int_t iDim1 = 0; iDim1 < 8; iDim1++)
2976 for (Int_t iDim2 = 0; iDim2 < 75 * N; iDim2++) {
2977 rh_segment[iChamber][iDim1][iDim2] = -999.;
2978 rhId_segment[iChamber][iDim1][iDim2] = -999;
2979 rh_segment_time[iChamber][iDim1][iDim2] = -999.;
2980 rh_sigm_segment[iChamber][iDim1][iDim2] = 1.;
2981 }
2982 }
2983 for (Int_t iSegment = 0; iSegment < nSegmentsMax; iSegment++)
2984 nSegments[iSegment] = 0;
2985 for (Int_t iDim = 0; iDim < 25 * N; iDim++) {
2986 x_mid[iDim] = -999.;
2987 y_mid[iDim] = -999.;
2988 pairID[iDim] = -999.;
2989 nhits[iDim] = 0.;
2990 aX[iDim] = -999.;
2991 aY[iDim] = -999.;
2992 leng[iDim] = -999.;
2993 imp[iDim] = -999.;
2994 }
2995}
2996
2997void BmnDchTrackFinder::Finish()
2998{
2999 // added
3000 //===============================================================================================================
3001
3002 if (fVerbose) {
3003 TFile file(fhTestFlnm.Data(), "RECREATE");
3004
3005 hEff1->Divide(hNomin1, hDenom1, 1., 1.);
3006 hEff2->Divide(hNomin2, hDenom2, 1., 1.);
3007
3008 fhList.Write();
3009 file.Close();
3010 }
3011 //===============================================================================================================
3012
3013 /// end
3014
3015 cout << "Work time of the DCH track finder: " << workTime << " s" << endl;
3016}
3017
3018Int_t BmnDchTrackFinder::Reconstruction(Int_t dchID,
3019 TString wire,
3020 Int_t pair,
3021 Int_t it_a,
3022 Int_t it_b,
3023 Double_t* wirenr_a,
3024 Double_t* wirenr_b,
3025 Double_t* time_a,
3026 Double_t* time_b,
3027 Bool_t* used_a,
3028 Bool_t* used_b,
3029 Double_t** _ab,
3030 Double_t** sigm_ab,
3031 Double_t** ab_time)
3032{
3033 // cout<<"61"<<endl;
3034 const Int_t arrIdxShift = (dchID == 2) ? 8 : 0;
3035 const Int_t arrIdxStart = (wire == "x") ? 0 : (wire == "y") ? 2 : (wire == "u") ? 4 : 6;
3036
3037 Double_t a_pm[2], b_pm[2];
3038
3039 for (Int_t i = 0; i < it_a; ++i)
3040 for (Int_t j = 0; j < it_b; ++j) {
3041 if (pair > 25 * N - 1)
3042 break;
3043 if ((wirenr_a[i] != wirenr_b[j] && wirenr_a[i] != wirenr_b[j] + 1))
3044 continue;
3045 Int_t func_nr_a = -1;
3046 Int_t func_nr_b = -1;
3047 for (Int_t t_it = 0; t_it < 3; t_it++) {
3048 if (time_a[i] >= t_dc[t_it][0 + arrIdxStart + arrIdxShift]
3049 && time_a[i] < t_dc[t_it + 1][0 + arrIdxStart + arrIdxShift])
3050 {
3051 func_nr_a = t_it;
3052 break;
3053 }
3054 }
3055 if (func_nr_a == -1)
3056 break;
3057 Double_t time = time_a[i] - t_dc[0][0 + arrIdxStart + arrIdxShift];
3058 Double_t d_a = 0;
3059 Double_t d_b = 0;
3060 d_a = scale
3061 * (pol_par_dc[func_nr_a][0][0 + arrIdxStart + arrIdxShift]
3062 + pol_par_dc[func_nr_a][1][0 + arrIdxStart + arrIdxShift] * time
3063 + pol_par_dc[func_nr_a][2][0 + arrIdxStart + arrIdxShift] * Power(time, 2)
3064 + pol_par_dc[func_nr_a][3][0 + arrIdxStart + arrIdxShift] * Power(time, 3)
3065 + pol_par_dc[func_nr_a][4][0 + arrIdxStart + arrIdxShift] * Power(time, 4));
3066 for (Int_t t_it = 0; t_it < 3; t_it++)
3067 if (time_b[j] >= t_dc[t_it][1 + arrIdxStart + arrIdxShift]
3068 && time_b[j] < t_dc[t_it + 1][1 + arrIdxStart + arrIdxShift])
3069 {
3070 func_nr_b = t_it;
3071 break;
3072 }
3073 if (func_nr_b == -1)
3074 continue;
3075 time = time_b[j] - t_dc[0][1 + arrIdxStart + arrIdxShift];
3076
3077 d_b = scale
3078 * (pol_par_dc[func_nr_b][0][1 + arrIdxStart + arrIdxShift]
3079 + pol_par_dc[func_nr_b][1][1 + arrIdxStart + arrIdxShift] * time
3080 + pol_par_dc[func_nr_b][2][1 + arrIdxStart + arrIdxShift] * Power(time, 2)
3081 + pol_par_dc[func_nr_b][3][1 + arrIdxStart + arrIdxShift] * Power(time, 3)
3082 + pol_par_dc[func_nr_b][4][1 + arrIdxStart + arrIdxShift] * Power(time, 4));
3083
3084 a_pm[0] = wirenr_a[i] - 119 + d_a;
3085 a_pm[1] = wirenr_a[i] - 119 - d_a;
3086 b_pm[0] = wirenr_b[j] - 118.5 + d_b;
3087 b_pm[1] = wirenr_b[j] - 118.5 - d_b;
3088
3089 if (fPeriod == 7 && fRunId > 3589 && fabs(d_a + d_b - .5) > .35)
3090 continue; // get rid of bad timing scale
3091
3092 Double_t dmin2 = LDBL_MAX;
3093 Double_t dmin1 = Abs(a_pm[0] - b_pm[0]); // set first pair as one of the minima
3094 Double_t a_coord1 = a_pm[0];
3095 Double_t b_coord1 = b_pm[0];
3096 Double_t a_coord2 = LDBL_MAX;
3097 Double_t b_coord2 = LDBL_MAX;
3098 for (Int_t k = 0; k < 2; k++) {
3099 for (Int_t m = 0; m < 2; m++) {
3100 if (k == 0 && m == 0)
3101 continue;
3102 if (Abs(a_pm[k] - b_pm[m]) < dmin1) {
3103 dmin2 = dmin1;
3104 dmin1 = Abs(a_pm[k] - b_pm[m]);
3105 a_coord2 = a_coord1;
3106 b_coord2 = b_coord1;
3107 a_coord1 = a_pm[k];
3108 b_coord1 = b_pm[m];
3109
3110 } else if (Abs(a_pm[k] - b_pm[m]) < dmin2 && Abs(a_pm[k] - b_pm[m]) != dmin1) {
3111 dmin2 = Abs(a_pm[k] - b_pm[m]);
3112 a_coord2 = a_pm[k];
3113 b_coord2 = b_pm[m];
3114 }
3115 }
3116 }
3117
3118 _ab[0][pair] = a_coord1;
3119 _ab[1][pair] = b_coord1;
3120 _ab[0][pair + 1] = a_coord2;
3121 _ab[1][pair + 1] = b_coord2;
3122
3123 ab_time[0][pair] = time_a[i];
3124 ab_time[1][pair] = time_b[j];
3125
3126 if (arrIdxStart == 0) {
3127 if (arrIdxShift == 0)
3128 hx1Da_Db->Fill(d_a + d_b - .5);
3129 if (arrIdxShift == 8)
3130 hx2Da_Db->Fill(d_a + d_b - .5);
3131 }
3132 if (arrIdxStart == 4) {
3133 if (arrIdxShift == 0)
3134 hu1Da_Db->Fill(d_a + d_b - .5);
3135 if (arrIdxShift == 8)
3136 hu2Da_Db->Fill(d_a + d_b - .5);
3137 }
3138 if (arrIdxStart == 6) {
3139 if (arrIdxShift == 0)
3140 hv1Da_Db->Fill(d_a + d_b - .5);
3141 if (arrIdxShift == 8)
3142 hv2Da_Db->Fill(d_a + d_b - .5);
3143 }
3144
3145 if (arrIdxStart == 2) {
3146 if (arrIdxShift == 0)
3147 hy1Da_Db->Fill(d_a + d_b - .5);
3148 if (arrIdxShift == 8)
3149 hy2Da_Db->Fill(d_a + d_b - .5);
3150 }
3151
3152 // add shifts for y and u
3153 // dc1
3154 if (arrIdxStart == 2 && arrIdxShift == 0) {
3155 _ab[0][pair] -= .461;
3156 _ab[1][pair] -= .461;
3157 }
3158 if (arrIdxStart == 4 && arrIdxShift == 0) {
3159 _ab[0][pair] -= 1.056;
3160 _ab[1][pair] -= 1.056;
3161 }
3162 // dc2
3163 if (arrIdxStart == 2 && arrIdxShift == 8) {
3164 _ab[0][pair] -= .346;
3165 _ab[1][pair] -= .346;
3166 }
3167 if (arrIdxStart == 4 && arrIdxShift == 8) {
3168 _ab[0][pair] -= .966;
3169 _ab[1][pair] -= .966;
3170 }
3171 // dc1
3172 if (arrIdxStart == 2 && arrIdxShift == 0) {
3173 _ab[0][pair + 1] -= .461;
3174 _ab[1][pair + 1] -= .461;
3175 }
3176 if (arrIdxStart == 4 && arrIdxShift == 0) {
3177 _ab[0][pair + 1] -= 1.056;
3178 _ab[1][pair + 1] -= 1.056;
3179 }
3180 // dc2
3181 if (arrIdxStart == 2 && arrIdxShift == 8) {
3182 _ab[0][pair + 1] -= .346;
3183 _ab[1][pair + 1] -= .346;
3184 }
3185 if (arrIdxStart == 4 && arrIdxShift == 8) {
3186 _ab[0][pair + 1] -= .966;
3187 _ab[1][pair + 1] -= .966;
3188 }
3189
3190 // end
3191 CompareDaDb(d_a, sigm_ab[0][pair]);
3192 CompareDaDb(d_b, sigm_ab[1][pair]);
3193 CompareDaDb(d_a, sigm_ab[0][pair + 1]);
3194 CompareDaDb(d_b, sigm_ab[1][pair + 1]);
3195
3196 pair += 2;
3197
3198 used_a[i] = kTRUE;
3199 used_b[j] = kTRUE;
3200 }
3201 return pair;
3202}
3203
3204Int_t BmnDchTrackFinder::ReconstructionSingle(Int_t dchID,
3205 TString wire,
3206 TString lay,
3207 Int_t single,
3208 Int_t it,
3209 Double_t* wirenr,
3210 Double_t* time_,
3211 Bool_t* used,
3212 Double_t** _single,
3213 Double_t** sigm_single)
3214{
3215 const Int_t arrIdxStart = (wire == "x") ? 0 : (wire == "y") ? 2 : (wire == "u") ? 4 : 6;
3216
3217 const Int_t arrIdx1 = (lay == "a") ? 0 : 1;
3218 const Int_t arrIdx2 = (dchID == 2) ? 8 : 0;
3219 const Double_t coeff = (lay == "a") ? 119 : 118.5;
3220
3221 for (Int_t i = 0; i < it; ++i) {
3222 if (used[i])
3223 continue;
3224
3225 Int_t func_nr = -1;
3226 for (Int_t t_it = 0; t_it < 3; t_it++) {
3227 if (time_[i] >= t_dc[t_it][0 + arrIdxStart + arrIdx1 + arrIdx2]
3228 && time_[i] < t_dc[t_it + 1][0 + arrIdxStart + arrIdx1 + arrIdx2])
3229 {
3230 func_nr = t_it;
3231 break;
3232 }
3233 }
3234 if (func_nr == -1)
3235 continue;
3236 Double_t time = time_[i] - t_dc[0][0 + arrIdxStart + arrIdx1 + arrIdx2];
3237 Double_t d = 0;
3238
3239 d = scale
3240 * (pol_par_dc[func_nr][0][0 + arrIdx1 + arrIdx2] + pol_par_dc[func_nr][1][0 + arrIdx1 + arrIdx2] * time
3241 + pol_par_dc[func_nr][2][0 + arrIdxStart + arrIdx1 + arrIdx2] * Power(time, 2)
3242 + pol_par_dc[func_nr][3][0 + arrIdxStart + arrIdx1 + arrIdx2] * Power(time, 3)
3243 + pol_par_dc[func_nr][4][0 + arrIdxStart + arrIdx1 + arrIdx2] * Power(time, 4));
3244
3245 _single[0 + arrIdx1][single] = wirenr[i] - coeff + d;
3246 _single[0 + arrIdx1][single + 1] = wirenr[i] - coeff - d;
3247
3248 // add shifts for y and u
3249 // dc1
3250 if (arrIdxStart == 2 && arrIdx2 == 0) {
3251 _single[0 + arrIdx1][single] -= .461;
3252 _single[0 + arrIdx1][single + 1] -= .461;
3253 }
3254 if (arrIdxStart == 4 && arrIdx2 == 0) {
3255 _single[0 + arrIdx1][single] -= 1.056;
3256 _single[0 + arrIdx1][single + 1] -= 1.056;
3257 }
3258 // dc2
3259 if (arrIdxStart == 2 && arrIdx2 == 8) {
3260 _single[0 + arrIdx1][single] -= .346;
3261 _single[0 + arrIdx1][single + 1] -= .346;
3262 }
3263 if (arrIdxStart == 4 && arrIdx2 == 8) {
3264 _single[0 + arrIdx1][single] -= .966;
3265 _single[0 + arrIdx1][single + 1] -= .966;
3266 }
3267
3268 CompareDaDb(d, sigm_single[0 + arrIdx1][single], sigm_single[0 + arrIdx1][single + 1]);
3269
3270 single += 2;
3271 }
3272 return single;
3273}
3274
3275Int_t BmnDchTrackFinder::ReconstructionMC(Int_t dchID,
3276 TString wire,
3277 Int_t pair,
3278 Int_t it_b,
3279 Int_t it_a,
3280 Double_t* wirenr_b,
3281 Double_t* wirenr_a,
3282 Double_t* time_b,
3283 Double_t* time_a,
3284 Int_t* hitId_b,
3285 Int_t* hitId_a,
3286 Bool_t* used_b,
3287 Bool_t* used_a,
3288 Double_t** _ab,
3289 Int_t** _abId,
3290 Double_t** sigm_ab,
3291 Double_t** ab_time)
3292{
3293 // cout<<"61"<<endl;
3294 const Int_t arrIdxShift = (dchID == 2) ? 8 : 0;
3295 const Int_t arrIdxStart = (wire == "x") ? 0 : (wire == "y") ? 2 : (wire == "u") ? 4 : 6;
3296
3297 Double_t a_pm[2], b_pm[2];
3298
3299 for (Int_t i = 0; i < it_a; ++i)
3300 for (Int_t j = 0; j < it_b; ++j) {
3301 if (pair > 25 * N - 1)
3302 break;
3303 if ((wirenr_a[i] != wirenr_b[j] && wirenr_a[i] != wirenr_b[j] + 1))
3304 continue;
3305
3306 Double_t d_a = time_a[i];
3307 Double_t d_b = time_b[j];
3308 a_pm[0] = wirenr_a[i] - 119 + d_a;
3309 a_pm[1] = wirenr_a[i] - 119 - d_a;
3310 b_pm[0] = wirenr_b[j] - 118.5 + d_b;
3311 b_pm[1] = wirenr_b[j] - 118.5 - d_b;
3312
3313 Double_t dmin2 = LDBL_MAX;
3314 Double_t dmin1 = Abs(a_pm[0] - b_pm[0]); // set first pair as one of the minima
3315 Double_t a_coord1 = a_pm[0];
3316 Double_t b_coord1 = b_pm[0];
3317 Double_t a_coord2 = LDBL_MAX;
3318 Double_t b_coord2 = LDBL_MAX;
3319 for (Int_t k = 0; k < 2; k++) {
3320 for (Int_t m = 0; m < 2; m++) {
3321 if (k == 0 && m == 0)
3322 continue;
3323 if (Abs(a_pm[k] - b_pm[m]) < dmin1) {
3324 dmin2 = dmin1;
3325 dmin1 = Abs(a_pm[k] - b_pm[m]);
3326 a_coord2 = a_coord1;
3327 b_coord2 = b_coord1;
3328 a_coord1 = a_pm[k];
3329 b_coord1 = b_pm[m];
3330
3331 } else if (Abs(a_pm[k] - b_pm[m]) < dmin2 && Abs(a_pm[k] - b_pm[m]) != dmin1) {
3332 dmin2 = Abs(a_pm[k] - b_pm[m]);
3333 a_coord2 = a_pm[k];
3334 b_coord2 = b_pm[m];
3335 }
3336 }
3337 }
3338
3339 _ab[0][pair] = a_coord1;
3340 _ab[1][pair] = b_coord1;
3341 _ab[0][pair + 1] = a_coord2;
3342 _ab[1][pair + 1] = b_coord2;
3343
3344 _abId[0][pair] = hitId_a[i];
3345 _abId[1][pair] = hitId_b[j];
3346 _abId[0][pair + 1] = hitId_a[i];
3347 _abId[1][pair + 1] = hitId_b[j];
3348
3349 ab_time[0][pair] = time_a[i];
3350 ab_time[1][pair] = time_b[j];
3351
3352 if (arrIdxStart == 0) {
3353 if (arrIdxShift == 0)
3354 hx1Da_Db->Fill(d_a + d_b - .5);
3355 if (arrIdxShift == 8)
3356 hx2Da_Db->Fill(d_a + d_b - .5);
3357 }
3358 if (arrIdxStart == 4) {
3359 if (arrIdxShift == 0)
3360 hu1Da_Db->Fill(d_a + d_b - .5);
3361 if (arrIdxShift == 8)
3362 hu2Da_Db->Fill(d_a + d_b - .5);
3363 }
3364 if (arrIdxStart == 6) {
3365 if (arrIdxShift == 0)
3366 hv1Da_Db->Fill(d_a + d_b - .5);
3367 if (arrIdxShift == 8)
3368 hv2Da_Db->Fill(d_a + d_b - .5);
3369 }
3370
3371 if (arrIdxStart == 2) {
3372 if (arrIdxShift == 0)
3373 hy1Da_Db->Fill(d_a + d_b - .5);
3374 if (arrIdxShift == 8)
3375 hy2Da_Db->Fill(d_a + d_b - .5);
3376 }
3377
3378 CompareDaDb(d_a, sigm_ab[0][pair]);
3379 CompareDaDb(d_b, sigm_ab[1][pair]);
3380 CompareDaDb(d_a, sigm_ab[0][pair + 1]);
3381 CompareDaDb(d_b, sigm_ab[1][pair + 1]);
3382
3383 pair += 2;
3384
3385 used_a[i] = kTRUE;
3386 used_b[j] = kTRUE;
3387 }
3388 return pair;
3389}
3390
3391Int_t BmnDchTrackFinder::ReconstructionSingleMC(Int_t dchID,
3392 TString wire,
3393 TString lay,
3394 Int_t single,
3395 Int_t it,
3396 Double_t* wirenr,
3397 Double_t* time_,
3398 Int_t* hitId,
3399 Bool_t* used,
3400 Double_t** _single,
3401 Int_t** _singleId,
3402 Double_t** sigm_single)
3403{
3404 // const Int_t arrIdxStart = (wire == "x") ? 0 : (wire == "y") ? 2 : (wire == "u") ? 4 : 6;
3405
3406 const Int_t arrIdx1 = (lay == "a") ? 0 : 1;
3407 // const Int_t arrIdx2 = (dchID == 2) ? 8 : 0;
3408 const Double_t coeff = (lay == "a") ? 119 : 118.5;
3409
3410 for (Int_t i = 0; i < it; ++i) {
3411 if (used[i])
3412 continue;
3413
3414 Double_t d = time_[i];
3415
3416 _single[0 + arrIdx1][single] = wirenr[i] - coeff + d;
3417 _single[0 + arrIdx1][single + 1] = wirenr[i] - coeff - d;
3418 _singleId[0 + arrIdx1][single] = hitId[i];
3419 _singleId[0 + arrIdx1][single + 1] = hitId[i];
3420
3421 CompareDaDb(d, sigm_single[0 + arrIdx1][single], sigm_single[0 + arrIdx1][single + 1]);
3422
3423 single += 2;
3424 }
3425 return single;
3426}
fit_seg
void fit_seg(Double_t *z_loc, Double_t *rh_seg, Double_t *rh_sigm_seg, Double_t *par_ab, Int_t skip_first, Int_t skip_second)
Definition BmnMath.cxx:589
workTime
static Double_t workTime
Definition BmnDchTrackFinder.cxx:31
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
i
int i
Definition P4_F32vec4.h:22
m
__m128 m
Definition P4_F32vec4.h:27
Bool_t
bool Bool_t
Int_t
int Int_t
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
BmnDchDigit::GetTime
Double_t GetTime() const
Definition BmnDchDigit.h:20
BmnDchDigit::GetWireNumber
Short_t GetWireNumber() const
Definition BmnDchDigit.h:19
BmnDchDigit::GetPlane
UChar_t GetPlane() const
Definition BmnDchDigit.h:18
BmnDchHit::GetDrift
Double_t GetDrift(void) const
Definition BmnDchHit.h:48
BmnDchHit::GetLayerNumber
UShort_t GetLayerNumber() const
Definition BmnDchHit.h:34
BmnDchHit::GetWirePosition
Double_t GetWirePosition(void) const
Definition BmnDchHit.h:50
BmnDchHit::GetHitId
Int_t GetHitId() const
Definition BmnDchHit.h:45
BmnDchTrackFinder::FindSegmentTrackMCId
void FindSegmentTrackMCId(Int_t, Int_t **, Double_t *, Int_t *, Int_t *)
Definition BmnDchTrackFinder.cxx:2312
BmnDchTrackFinder::fhList
TList fhList
chi2 of the matched seg
Definition BmnDchTrackFinder.h:175
BmnDchTrackFinder::FillPlaneResiduals
void FillPlaneResiduals(Int_t, Int_t *, Double_t **, Double_t **, Double_t **, Double_t *)
Definition BmnDchTrackFinder.cxx:1861
BmnDchTrackFinder::fBmnDchDigitsArray
TClonesArray * fBmnDchDigitsArray
Input array.
Definition BmnDchTrackFinder.h:72
BmnDchTrackFinder::y_mid
Double_t * y_mid
x glob of matched segment in the z situated between the two DCH chambers
Definition BmnDchTrackFinder.h:165
BmnDchTrackFinder::ReconstructionSingleMC
Int_t ReconstructionSingleMC(Int_t, TString, TString, Int_t, Int_t, Double_t *, Double_t *, Int_t *, Bool_t *, Double_t **, Int_t **, Double_t **)
Definition BmnDchTrackFinder.cxx:3391
BmnDchTrackFinder::pol_par_dc
Double_t pol_par_dc[3][5][16]
Definition BmnDchTrackFinder.h:89
BmnDchTrackFinder::Sigm_y_single
Double_t *** Sigm_y_single
Definition BmnDchTrackFinder.h:150
BmnDchTrackFinder::nhits
Int_t * nhits
segments numbers that were matche into global track in form of 1000 * (DC1_segIt + 1) + DC_segIt + 1
Definition BmnDchTrackFinder.h:168
BmnDchTrackFinder::FitDchSegments
Bool_t FitDchSegments(Int_t, Int_t *, Double_t **, Double_t **, Double_t **, Double_t *, Double_t *, Double_t *)
Definition BmnDchTrackFinder.cxx:1703
BmnDchTrackFinder::FillSegmentParametersSigmas
void FillSegmentParametersSigmas(Int_t, Double_t **, Double_t **, Double_t *, Double_t **)
Definition BmnDchTrackFinder.cxx:2353
BmnDchTrackFinder::rh_sigm_segment
Double_t *** rh_sigm_segment
Definition BmnDchTrackFinder.h:157
BmnDchTrackFinder::Exec
virtual void Exec(Option_t *opt)
Definition BmnDchTrackFinder.cxx:867
BmnDchTrackFinder::aX
Double_t * aX
total number of hits from the two segments that were matched
Definition BmnDchTrackFinder.h:169
BmnDchTrackFinder::ReconstructionMC
Int_t ReconstructionMC(Int_t, TString, Int_t, Int_t, Int_t, Double_t *, Double_t *, Double_t *, Double_t *, Int_t *, Int_t *, Bool_t *, Bool_t *, Double_t **, Int_t **, Double_t **, Double_t **)
Definition BmnDchTrackFinder.cxx:3275
BmnDchTrackFinder::SelectLongestAndBestSegments
void SelectLongestAndBestSegments(Int_t, Int_t *, Double_t **, Double_t *)
Definition BmnDchTrackFinder.cxx:1826
BmnDchTrackFinder::CompareDaDb
void CompareDaDb(Double_t, Double_t &)
Definition BmnDchTrackFinder.cxx:1807
BmnDchTrackFinder::aY
Double_t * aY
x slope of the matched segment
Definition BmnDchTrackFinder.h:170
BmnDchTrackFinder::CalculateResidual
Double_t CalculateResidual(Int_t, Int_t, Double_t **, Double_t **)
Definition BmnDchTrackFinder.cxx:2683
BmnDchTrackFinder::leng
Double_t * leng
reconstructed particle trajectory momentum
Definition BmnDchTrackFinder.h:172
BmnDchTrackFinder::ReconstructionSingle
Int_t ReconstructionSingle(Int_t, TString, TString, Int_t, Int_t, Double_t *, Double_t *, Bool_t *, Double_t **, Double_t **)
Definition BmnDchTrackFinder.cxx:3204
BmnDchTrackFinder::BuildUVSegments
Int_t BuildUVSegments(Int_t, Int_t, Int_t, Int_t, Int_t, Int_t, Int_t, Int_t, Int_t, Double_t **, Double_t **, Double_t **, Double_t **, Int_t **, Int_t **, Int_t **, Int_t **, Double_t **, Double_t **, Double_t **, Double_t **, Double_t **, Double_t **, Double_t **, Double_t **, Double_t **, Int_t **, Double_t **, Double_t **, Double_t **, Double_t **, Int_t **, Int_t **, Double_t **, Double_t **)
Definition BmnDchTrackFinder.cxx:1490
BmnDchTrackFinder::CalculateResidualMatch
Double_t CalculateResidualMatch(Int_t, Int_t, Int_t, Double_t ***, Double_t ***)
Definition BmnDchTrackFinder.cxx:2666
BmnDchTrackFinder::Sigm_x_single
Double_t *** Sigm_x_single
Definition BmnDchTrackFinder.h:151
BmnDchTrackFinder::Init
virtual InitStatus Init()
Definition BmnDchTrackFinder.cxx:2696
BmnDchTrackFinder::Reconstruction
Int_t Reconstruction(Int_t, TString, Int_t, Int_t, Int_t, Double_t *, Double_t *, Double_t *, Double_t *, Bool_t *, Bool_t *, Double_t **, Double_t **, Double_t **)
Definition BmnDchTrackFinder.cxx:3018
BmnDchTrackFinder::t_dc
Double_t t_dc[4][16]
Definition BmnDchTrackFinder.h:88
BmnDchTrackFinder::Sigm_u_single
Double_t *** Sigm_u_single
Definition BmnDchTrackFinder.h:149
BmnDchTrackFinder::fDchTracks
TClonesArray * fDchTracks
Definition BmnDchTrackFinder.h:77
BmnDchTrackFinder::params_sigmas
Double_t *** params_sigmas
Definition BmnDchTrackFinder.h:153
BmnDchTrackFinder::rh_segment_time
Double_t *** rh_segment_time
Definition BmnDchTrackFinder.h:156
BmnDchTrackFinder::Sigm_v_single
Double_t *** Sigm_v_single
Definition BmnDchTrackFinder.h:148
BmnDchTrackFinder::fBmnHitsArray
TClonesArray * fBmnHitsArray
Definition BmnDchTrackFinder.h:73
BmnDchTrackFinder::BuildXYSegments
Int_t BuildXYSegments(Int_t, Int_t, Int_t, Int_t, Int_t, Int_t, Int_t, Int_t, Int_t, Double_t **, Double_t **, Double_t **, Double_t **, Int_t **, Int_t **, Int_t **, Int_t **, Double_t **, Double_t **, Double_t **, Double_t **, Double_t **, Double_t **, Double_t **, Double_t **, Double_t **, Int_t **, Double_t **, Double_t **, Double_t **, Double_t **, Int_t **, Int_t **, Double_t **, Double_t **)
Definition BmnDchTrackFinder.cxx:1274
BmnDchTrackFinder::imp
Double_t * imp
y slope of the matched segment
Definition BmnDchTrackFinder.h:171
BmnDchTrackFinder::pairID
Int_t * pairID
y glob of matched segment in the z situated between the two DCH chambers
Definition BmnDchTrackFinder.h:166
BmnDchTrackFinder::Chi2_match
Double_t * Chi2_match
the distance from z = 0 to the global point of the matched segment
Definition BmnDchTrackFinder.h:173
BmnDchTrack::SetTrackId
void SetTrackId(Int_t itr)
Definition BmnDchTrack.h:19
BmnTrack::SetChi2
void SetChi2(Double_t chi2)
Definition BmnTrack.h:76
BmnTrack::SetFlag
void SetFlag(Int_t flag)
Definition BmnTrack.h:74
BmnTrack::SetNHits
void SetNHits(Int_t n)
Definition BmnTrack.h:80
BmnTrack::SetParamFirst
void SetParamFirst(FairTrackParam &par)
Modifiers.
Definition BmnTrack.h:70
BmnTrack::SetLength
void SetLength(Double_t length)
Definition BmnTrack.h:84
TFile::Close
void Close(Option_t *option="") override
sqrt
RVec< PromoteType< T > > sqrt(const RVec< T > &v)
atan
RVec< PromoteType< T > > atan(const RVec< T > &v)
TMath::Power
Double_t Power(Double_t x, Double_t y)
TMath::Abs
Double_t Abs(Double_t d)