BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
recoRunDCH_VP.C
Go to the documentation of this file.
1/*-------------------------------*/
2// macro that adds DCH segments to the root file
3// indicated at line number 310;
4// the files required at the moment:
5// bmn_run00%d_digi.root - containing the DCH digis
6// transfer_func.txt - containing the coeficients of the transfer functions for each layer
7//
8// the pol_coord0%d.txt that correspond to each file will be added soon after the raw data is available
9// the script doesn't compile for now and runs only in the interpreter mode (i.e. root -b recoRunDCH.C)
10// more optimization and feature adding to be done further
11/*-------------------------------*/
12
13#include "TCanvas.h"
14#include "TChain.h"
15#include "TClonesArray.h"
16#include "TFile.h"
17#include "TGeoManager.h"
18#include "TH1F.h"
19#include "TH2F.h"
20#include "TMath.h"
21#include "TVector3.h"
22
23#include <vector>
24// #include "TROOT.h"
25// MK
26
27// #include "FairTrackParam.h"
28// #include "BmnMwpcDigit.h"
29// #include "CbmTrack.h"
30
31#include "TString.h"
32// #include "TTree.h"
33// #include "TProfile2D.h"
34
35#include "TArrayS.h"
36#include "TBenchmark.h"
37#include "TClonesArray.h"
38#include "TLine.h"
39#include "TStyle.h"
40
41#include <Riostream.h>
42#include <algorithm>
43#include <ctime>
44#include <fstream>
45#include <iostream>
46#include <stdlib.h>
47// #include "eig3.cpp"
48// #include "eig3.h"
49
50#include <TPostScript.h>
51#include <TProfile.h>
52// #include <Math/Vector3D.h>
53// #include "eig3.h"
54// #include "eig3.cpp"
55// #include "mylib.cpp"
56// #include "6x6inv.cpp"
57
58// #include "BmnDchDigit.h"
59
60#include <vector>
61// #include <set>
62
63using namespace std;
64using namespace TMath;
65
66//
67// runType - "run1", "run2", "run3"
68//
69
70// some global constants
71
72const Float_t z_loc[8] = {9.3, 8.1, 3.5, 2.3, -2.3, -3.5, -8.1, -9.3}; // z local xa->vb (cm)
73const Float_t z_glob[16] = {-45.7, -46.9, -51.5, -52.7, -57.3, -58.5, -63.1, -64.3,
74 64.3, 63.1, 58.5, 57.3, 52.7, 51.5, 46.9, 45.7}; // z global dc 1 & dc 2 (cm)
75const Float_t Z_dch1 = 529.65; // z coord in the midle of dch1
76const Float_t Z_dch2 = 727.65; // z coord in the midle of dch2
77const Float_t Z_dch_mid =
78 628.65; // z coord between the two chambers, this z is considered the global z for the matched dch segment
79
80double sqrt_2 = sqrt(2.);
81double isqrt_2 = 1 / sqrt_2;
82// shift variables
83float intBL = 2; // switch for different fields
84
85/*
86// run2
87const Float_t x1_sh = -5.2 ;
88const Float_t x2_sh = -7.8 ;
89const Float_t y1_sh = -0.4;
90const Float_t y2_sh = -0.7;
91
92// run3
93const Float_t x1_sh = -4.6 - 21;//-5.2
94const Float_t x2_sh = -8.4 - 25;//-7.8
95const Float_t y1_sh = -5.7;
96const Float_t y2_sh = -5.8;
97
98// run4
99
100const Float_t x1_sh = 15.4;//
101const Float_t x2_sh = 11.8;//
102const Float_t y1_sh = 0.2;
103const Float_t y2_sh = 7.03;
104
105*/
106const Float_t x1_sh = 15.0; //
107const Float_t x2_sh = 12.2; //
108const Float_t y1_sh = 0.2;
109const Float_t y2_sh = -0.3;
110
111const Float_t u1_sh = isqrt_2 * (y1_sh - x1_sh);
112const Float_t u2_sh = isqrt_2 * (y2_sh - x2_sh);
113const Float_t v1_sh = isqrt_2 * (x1_sh + y1_sh);
114const Float_t v2_sh = isqrt_2 * (x2_sh + y2_sh);
115
116const Float_t x1_slope_sh = -0.08;
117const Float_t y1_slope_sh = 0.051; // was .05
118const Float_t x2_slope_sh = -0.08;
119const Float_t y2_slope_sh = 0.055; //
120
121const Float_t u1_slope_sh = isqrt_2 * (y1_slope_sh - x1_slope_sh);
122const Float_t u2_slope_sh = isqrt_2 * (y2_slope_sh - x2_slope_sh);
123const Float_t v1_slope_sh = isqrt_2 * (y1_slope_sh + x1_slope_sh);
124const Float_t v2_slope_sh = isqrt_2 * (y2_slope_sh + x2_slope_sh);
125
126void fit_seg(float* rh_seg, float* rh_sigm_seg, float* par_ab, int skip_first, int skip_second)
127{
128
129 // linear fit
130 float z_loc[8] = {9.3, 8.1, 3.5, 2.3, -2.3, -3.5, -8.1, -9.3}; // z local xa->vb (cm)
131 float A[4][4] = {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}; // coef matrix
132 float f[4] = {0}; // free coef
133 // float sigm_sq[8] = {1,1,1,1,1,1,1,1};
134 int h[8] = {0, 0, 0, 0, 0, 0, 0, 0};
135
136 for (int i = 0; i < 4; i++) {
137 par_ab[i] = 999;
138 }
139 float xDC_glob, yDC_glob = -999;
140
141 for (int i = 0; i < 8; i++) {
142 h[i] = 1;
143 // out1<<"setting h[i]"<<endl;
144 if (i == skip_first || i == skip_second || rh_seg[i] == -999) {
145 h[i] = 0;
146 }
147 } // i
148
149 A[0][0] = 2 * z_loc[0] * z_loc[0] * h[0] / rh_sigm_seg[0] + z_loc[4] * z_loc[4] * h[4] / rh_sigm_seg[4]
150 + z_loc[6] * z_loc[6] * h[6] / rh_sigm_seg[6] + 2 * z_loc[1] * z_loc[1] * h[1] / rh_sigm_seg[1]
151 + z_loc[5] * z_loc[5] * h[5] / rh_sigm_seg[5] + z_loc[7] * z_loc[7] * h[7] / rh_sigm_seg[7]; // aX_a
152
153 A[0][1] = 2 * z_loc[0] * h[0] / rh_sigm_seg[0] + z_loc[4] * h[4] / rh_sigm_seg[4] + z_loc[6] * h[6] / rh_sigm_seg[6]
154 + 2 * z_loc[1] * h[1] / rh_sigm_seg[1] + z_loc[5] * h[5] / rh_sigm_seg[5]
155 + z_loc[7] * h[7] / rh_sigm_seg[7]; // bX_a
156
157 A[0][2] = z_loc[6] * z_loc[6] * h[6] / rh_sigm_seg[6] - z_loc[4] * z_loc[4] * h[4] / rh_sigm_seg[4]
158 + z_loc[7] * z_loc[7] * h[7] / rh_sigm_seg[7] - z_loc[5] * z_loc[5] * h[5] / rh_sigm_seg[5]; // aY
159
160 A[0][3] = z_loc[6] * h[6] / rh_sigm_seg[6] - z_loc[4] * h[4] / rh_sigm_seg[4] + z_loc[7] * h[7] / rh_sigm_seg[7]
161 - z_loc[5] * h[5] / rh_sigm_seg[5]; // bY_a
162
163 // dChi2/d_b_x
164
165 A[1][0] = 2 * z_loc[0] * h[0] / rh_sigm_seg[0] + z_loc[4] * h[4] / rh_sigm_seg[4] + z_loc[6] * h[6] / rh_sigm_seg[6]
166 + 2 * z_loc[1] * h[1] / rh_sigm_seg[1] + z_loc[5] * h[5] / rh_sigm_seg[5]
167 + z_loc[7] * h[7] / rh_sigm_seg[7];
168
169 A[1][1] = 2 * h[0] / rh_sigm_seg[0] + 1 * h[4] / rh_sigm_seg[4] + 1 * h[6] / rh_sigm_seg[6]
170 + 2 * h[1] / rh_sigm_seg[1] + 1 * h[5] / rh_sigm_seg[5] + 1 * h[7] / rh_sigm_seg[7]; // bX_a
171
172 A[1][2] = z_loc[6] * h[6] / rh_sigm_seg[6] - z_loc[4] * h[4] / rh_sigm_seg[4] + z_loc[7] * h[7] / rh_sigm_seg[7]
173 - z_loc[5] * h[5] / rh_sigm_seg[5]; // aY_a
174
175 A[1][3] = 1 * h[7] / rh_sigm_seg[7] - 1 * h[5] / rh_sigm_seg[5] + 1 * h[6] / rh_sigm_seg[6]
176 - 1 * h[4] / rh_sigm_seg[4]; // bY_a
177
178 // dChi2/da_y
179
180 A[2][0] = z_loc[6] * z_loc[6] * h[6] / rh_sigm_seg[6] - z_loc[4] * z_loc[4] * h[4] / rh_sigm_seg[4]
181 + z_loc[7] * z_loc[7] * h[7] / rh_sigm_seg[7] - z_loc[5] * z_loc[5] * h[5] / rh_sigm_seg[5]; // aX_a
182
183 A[2][1] = z_loc[6] * h[6] / rh_sigm_seg[6] - z_loc[4] * h[4] / rh_sigm_seg[4] + z_loc[7] * h[7] / rh_sigm_seg[7]
184 - z_loc[5] * h[5] / rh_sigm_seg[5]; // bX_a
185
186 A[2][2] = 2 * z_loc[2] * z_loc[2] * h[2] / rh_sigm_seg[2] + z_loc[4] * z_loc[4] * h[4] / rh_sigm_seg[4]
187 + z_loc[6] * z_loc[6] * h[6] / rh_sigm_seg[6] + 2 * z_loc[3] * z_loc[3] * h[3] / rh_sigm_seg[3]
188 + z_loc[5] * z_loc[5] * h[5] / rh_sigm_seg[5] + z_loc[7] * z_loc[7] * h[7] / rh_sigm_seg[7]; // aY_a
189
190 A[2][3] = 2 * z_loc[2] * h[2] / rh_sigm_seg[2] + z_loc[4] * h[4] / rh_sigm_seg[4] + z_loc[6] * h[6] / rh_sigm_seg[6]
191 + 2 * z_loc[3] * h[3] / rh_sigm_seg[3] + z_loc[5] * h[5] / rh_sigm_seg[5]
192 + z_loc[7] * h[7] / rh_sigm_seg[7];
193
194 ////dChi2/db_y
195
196 A[3][0] = z_loc[6] * h[6] / rh_sigm_seg[6] - z_loc[4] * h[4] / rh_sigm_seg[4] + z_loc[7] * h[7] / rh_sigm_seg[7]
197 - z_loc[5] * h[5] / rh_sigm_seg[5]; // aX_a
198
199 A[3][1] = 1 * h[6] / rh_sigm_seg[6] - 1 * h[4] / rh_sigm_seg[4] + 1 * h[7] / rh_sigm_seg[7]
200 - 1 * h[5] / rh_sigm_seg[5]; // bX_a
201
202 A[3][2] = 2 * z_loc[2] * h[2] / rh_sigm_seg[2] + z_loc[4] * h[4] / rh_sigm_seg[4] + z_loc[6] * h[6] / rh_sigm_seg[6]
203 + 2 * z_loc[3] * h[3] / rh_sigm_seg[3] + z_loc[5] * h[5] / rh_sigm_seg[5]
204 + z_loc[7] * h[7] / rh_sigm_seg[7]; // aY_a
205
206 A[3][3] = 2 * h[2] / rh_sigm_seg[2] + 1 * h[4] / rh_sigm_seg[4] + 1 * h[6] / rh_sigm_seg[6]
207 + 2 * h[3] / rh_sigm_seg[3] + 1 * h[5] / rh_sigm_seg[5] + 1 * h[7] / rh_sigm_seg[7]; // bY_a
208
209 // free coef
210
211 // dChi2/da_x
212
213 f[0] = 2 * z_loc[0] * rh_seg[0] * h[0] / rh_sigm_seg[0] + sqrt_2 * z_loc[6] * rh_seg[6] * h[6] / rh_sigm_seg[6]
214 - sqrt_2 * z_loc[4] * rh_seg[4] * h[4] / rh_sigm_seg[4] + 2 * z_loc[1] * rh_seg[1] * h[1] / rh_sigm_seg[1]
215 + sqrt_2 * z_loc[7] * rh_seg[7] * h[7] / rh_sigm_seg[7]
216 - sqrt_2 * z_loc[5] * rh_seg[5] * h[5] / rh_sigm_seg[5]; // j = nr of seg
217 // dChi2/db_x
218 f[1] = 2 * rh_seg[0] * h[0] / rh_sigm_seg[0] + sqrt_2 * rh_seg[6] * h[6] / rh_sigm_seg[6]
219 - sqrt_2 * rh_seg[4] * h[4] / rh_sigm_seg[4] + 2 * rh_seg[1] * h[1] / rh_sigm_seg[1]
220 + sqrt_2 * rh_seg[7] * h[7] / rh_sigm_seg[7] - sqrt_2 * rh_seg[5] * h[5] / rh_sigm_seg[5]; // j = nr of seg
221 // dChi2/da_y
222 f[2] = 2 * z_loc[2] * rh_seg[2] * h[2] / rh_sigm_seg[2] + sqrt_2 * z_loc[6] * rh_seg[6] * h[6] / rh_sigm_seg[6]
223 + sqrt_2 * z_loc[4] * rh_seg[4] * h[4] / rh_sigm_seg[4] + 2 * z_loc[3] * rh_seg[3] * h[3] / rh_sigm_seg[3]
224 + sqrt_2 * z_loc[7] * rh_seg[7] * h[7] / rh_sigm_seg[7]
225 + sqrt_2 * z_loc[5] * rh_seg[5] * h[5] / rh_sigm_seg[5]; // j = nr of seg
226 ////dChi2/db_y
227 f[3] = 2 * rh_seg[2] * h[2] / rh_sigm_seg[2] + sqrt_2 * rh_seg[6] * h[6] / rh_sigm_seg[6]
228 + sqrt_2 * rh_seg[4] * h[4] / rh_sigm_seg[4] + 2 * rh_seg[3] * h[3] / rh_sigm_seg[3]
229 + sqrt_2 * rh_seg[7] * h[7] / rh_sigm_seg[7] + sqrt_2 * rh_seg[5] * h[5] / rh_sigm_seg[5]; // j = nr of seg
230
231 // inverse the matrix
232
233 /**** Gaussian algorithm for 4x4 matrix inversion ****/
234 Int_t i1, j1, k1, l1;
235 Double_t factor;
236 Double_t temp[4];
237 Double_t b[4][4];
238 Double_t A0[4][4];
239
240 for (i1 = 0; i1 < 4; i1++)
241 for (j1 = 0; j1 < 4; j1++)
242 A0[i1][j1] = A[i1][j1];
243
244 // Set b to I
245 for (i1 = 0; i1 < 4; i1++)
246 for (j1 = 0; j1 < 4; j1++)
247 if (i1 == j1)
248 b[i1][j1] = 1.0;
249 else
250 b[i1][j1] = 0.0;
251
252 for (i1 = 0; i1 < 4; i1++) {
253 for (j1 = i1 + 1; j1 < 4; j1++)
254 if (fabs(A[i1][i1]) < fabs(A[j1][i1])) {
255 for (l1 = 0; l1 < 4; l1++)
256 temp[l1] = A[i1][l1];
257 for (l1 = 0; l1 < 4; l1++)
258 A[i1][l1] = A[j1][l1];
259 for (l1 = 0; l1 < 4; l1++)
260 A[j1][l1] = temp[l1];
261 for (l1 = 0; l1 < 4; l1++)
262 temp[l1] = b[i1][l1];
263 for (l1 = 0; l1 < 4; l1++)
264 b[i1][l1] = b[j1][l1];
265 for (l1 = 0; l1 < 4; l1++)
266 b[j1][l1] = temp[l1];
267 }
268 factor = A[i1][i1];
269 for (j1 = 4 - 1; j1 > -1; j1--) {
270 b[i1][j1] /= factor;
271 A[i1][j1] /= factor;
272 }
273 for (j1 = i1 + 1; j1 < 4; j1++) {
274 factor = -A[j1][i1];
275 for (k1 = 0; k1 < 4; k1++) {
276 A[j1][k1] += A[i1][k1] * factor;
277 b[j1][k1] += b[i1][k1] * factor;
278 }
279 }
280 }
281 for (i1 = 3; i1 > 0; i1--) {
282 for (j1 = i1 - 1; j1 > -1; j1--) {
283 factor = -A[j1][i1];
284 for (k1 = 0; k1 < 4; k1++) {
285 A[j1][k1] += A[i1][k1] * factor;
286 b[j1][k1] += b[i1][k1] * factor;
287 }
288 }
289 }
290
291 float sum;
292
293 float A1[4][4] = {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}};
294
295 for (i1 = 0; i1 < 4; ++i1)
296 for (j1 = 0; j1 < 4; ++j1) {
297 sum = 0;
298
299 for (k1 = 0; k1 < 4; ++k1)
300 sum += A0[i1][k1] * b[k1][j1];
301 A1[i1][j1] = sum;
302 }
303
304 for (i1 = 0; i1 < 4; i1++) {
305 par_ab[i1] = 0;
306 for (j1 = 0; j1 < 4; j1++) {
307 par_ab[i1] += b[i1][j1] * f[j1];
308 }
309 }
310} // fit_seg
311
312void recoRunDCH_VP(Int_t runId = 810)
313{
314
315 gROOT->LoadMacro("$VMCWORKDIR/macro/run/bmnloadlibs.C");
316 bmnloadlibs(); // load bmn libraries
317
318 TChain* bmnTree = new TChain("cbmsim");
319 bmnTree->Add(TString::Format("bmn_run0%d_digi.root", runId));
320
321 TClonesArray* dchDigits;
322 bmnTree->SetBranchAddress("DCH", &dchDigits);
323 Int_t startEvent = 1;
324 Int_t nEvents = bmnTree->GetEntries();
325 nEvents = 1000;
326
327 TString outName = Form("bmn_reco_DCH_run%04d.root", runId);
328 TFile* fReco = new TFile(outName, "RECREATE"); // the root file where the dch hits and tracks will be filled
329 TTree* tReco = new TTree("cbmsim", "test_bmn");
330
331 TClonesArray* dchHits = new TClonesArray("BmnDchHit");
332 TClonesArray* dchTracks = new TClonesArray("BmnTrack");
333 TClonesArray* dchMatchedTracks = new TClonesArray("BmnTrack");
334
335 tReco->Branch("BmnDchHit", &dchHits);
336 tReco->Branch("DchTracks", &dchTracks);
337 tReco->Branch("DchMatchedTracks", &dchMatchedTracks);
338
339 // read and memorize the transfer fuctions' parameters
340
341 float p1_xa, p0_xa, p1_xb, p0_xb, p1_ya, p0_ya, p1_yb, p0_yb, p1_ua, p0_ua, p1_ub, p0_ub, p1_va, p0_va, p1_vb,
342 p0_vb;
343
344 double t_dc[5][16]; //[time interval][plane number]
345 double pol_par_dc[3][5][16]; //[polinom number][param number][plane number]
346 double scaling[16]; //[plane number]
347
348 ifstream fin;
349 // fin.open(TString::Format("pol_coord00%d.txt", runId),ios::in);
350 fin.open("transfer_func.txt", ios::in);
351 for (int fi = 0; fi < 16; fi++) {
352 fin >> t_dc[0][fi] >> t_dc[1][fi] >> t_dc[2][fi] >> t_dc[3][fi] >> t_dc[4][fi] >> pol_par_dc[0][0][fi]
353 >> pol_par_dc[0][1][fi] >> pol_par_dc[0][2][fi] >> pol_par_dc[0][3][fi] >> pol_par_dc[0][4][fi]
354 >> pol_par_dc[1][0][fi] >> pol_par_dc[1][1][fi] >> pol_par_dc[1][2][fi] >> pol_par_dc[1][3][fi]
355 >> pol_par_dc[1][4][fi] >> pol_par_dc[2][0][fi] >> pol_par_dc[2][1][fi] >> pol_par_dc[2][2][fi]
356 >> pol_par_dc[2][3][fi] >> pol_par_dc[2][4][fi] >> scaling[fi];
357 }
358
359 fin.close();
360
361 float scale[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
362 for (Int_t sc = 0; sc < 16; sc++) {
363 scale[sc] = 0.5; // / scaling[sc];
364 } // scaling for transfer function coeffs
365
366 int prev_wire = -1;
367 int prev_time = -1;
368 int skipped_ev = 0;
369 int poss = 0;
370 int imposs = 0;
371 for (Int_t iEv = startEvent; iEv < startEvent + nEvents; iEv++) {
372
373 bmnTree->GetEntry(iEv);
374
375 dchHits->Clear();
376 dchTracks->Clear();
377 dchMatchedTracks->Clear();
378
379 if (iEv % 1000 == 0)
380 cout << "Event/last : " << iEv + 1 << "/" << startEvent + nEvents << endl;
381
382 BmnDchDigit* digit = NULL;
383
384 // temporary containers
385 double time_xa1[20] = {-99};
386 double time_xa2[20] = {-99};
387 double time_xb1[20] = {-99};
388 double time_xb2[20] = {-99};
389 double time_ya1[20] = {-99};
390 double time_ya2[20] = {-99};
391 double time_yb1[20] = {-99};
392 double time_yb2[20] = {-99};
393 double time_ua1[20] = {-99};
394 double time_ua2[20] = {-99};
395 double time_ub1[20] = {-99};
396 double time_ub2[20] = {-99};
397 double time_va1[20] = {-99};
398 double time_va2[20] = {-99};
399 double time_vb1[20] = {-99};
400 double time_vb2[20] = {-99};
401
402 float wirenr_xa1[20] = {-99};
403 float wirenr_xa2[20] = {-99};
404 float wirenr_xb1[20] = {-99};
405 float wirenr_xb2[20] = {-99};
406 float wirenr_ya1[20] = {-99};
407 float wirenr_ya2[20] = {-99};
408 float wirenr_yb1[20] = {-99};
409 float wirenr_yb2[20] = {-99};
410 float wirenr_ua1[20] = {-99};
411 float wirenr_ua2[20] = {-99};
412 float wirenr_ub1[20] = {-99};
413 float wirenr_ub2[20] = {-99};
414 float wirenr_va1[20] = {-99};
415 float wirenr_va2[20] = {-99};
416 float wirenr_vb1[20] = {-99};
417 float wirenr_vb2[20] = {-99};
418
419 bool used_xa1[20] = {0};
420 bool used_xa2[20] = {0};
421 bool used_xb1[20] = {0};
422 bool used_xb2[20] = {0};
423 bool used_ya1[20] = {0};
424 bool used_ya2[20] = {0};
425 bool used_yb1[20] = {0};
426 bool used_yb2[20] = {0};
427 bool used_ua1[20] = {0};
428 bool used_ua2[20] = {0};
429 bool used_ub1[20] = {0};
430 bool used_ub2[20] = {0};
431 bool used_va1[20] = {0};
432 bool used_va2[20] = {0};
433 bool used_vb1[20] = {0};
434 bool used_vb2[20] = {0};
435
436 int it_xa1 = 0;
437 int it_xa2 = 0;
438 int it_xb1 = 0;
439 int it_xb2 = 0;
440 int it_ya1 = 0;
441 int it_ya2 = 0;
442 int it_yb1 = 0;
443 int it_yb2 = 0;
444 int it_ua1 = 0;
445 int it_ua2 = 0;
446 int it_ub1 = 0;
447 int it_ub2 = 0;
448 int it_va1 = 0;
449 int it_va2 = 0;
450 int it_vb1 = 0;
451 int it_vb2 = 0;
452
453 bool goodEv = true;
454 bool written = false;
455
456 for (Int_t iDig = 0; iDig < dchDigits->GetEntriesFast(); ++iDig) {
457
458 digit = (BmnDchDigit*)dchDigits->At(iDig);
459 Short_t plane = digit->GetPlane();
460
461 // skip identical events
462 if (!written) {
463 written = true;
464 if (digit->GetTime() == prev_time && digit->GetWireNumber() == prev_wire) {
465 goodEv = false;
466 skipped_ev++;
467 } else {
468 prev_time = int(digit->GetTime());
469 prev_wire = int(digit->GetWireNumber());
470 }
471 } //! written
472
473 double times = digit->GetTime();
474 bool secondaries = false;
475 switch (plane) {
476 case 0:
477 for (int sec = 0; sec < it_va1 - 1; sec++) {
478 if (digit->GetWireNumber() == wirenr_va1[sec]) {
479 secondaries = true;
480 break;
481 }
482 } // skip secondary hits
483 if (it_va1 == 19 || secondaries)
484 break;
485 wirenr_va1[it_va1] = digit->GetWireNumber();
486 time_va1[it_va1] = times;
487 it_va1++;
488 break;
489 case 1:
490
491 for (int sec = 0; sec < it_vb1 - 1; sec++) {
492 if (digit->GetWireNumber() == wirenr_vb1[sec]) {
493 secondaries = true;
494 break;
495 }
496 } // skip secondary hits
497 if (it_vb1 == 19 || secondaries)
498 break;
499 wirenr_vb1[it_vb1] = digit->GetWireNumber();
500 time_vb1[it_vb1] = times;
501 it_vb1++;
502 break;
503 case 2:
504 for (int sec = 0; sec < it_ua1 - 1; sec++) {
505 if (digit->GetWireNumber() == wirenr_ua1[sec]) {
506 secondaries = true;
507 break;
508 }
509 } // skip secondary hits
510
511 if (it_ua1 == 19 || secondaries)
512 break;
513 wirenr_ua1[it_ua1] = digit->GetWireNumber();
514 time_ua1[it_ua1] = times;
515 it_ua1++;
516 break;
517 case 3:
518 for (int sec = 0; sec < it_ub1 - 1; sec++) {
519 if (digit->GetWireNumber() == wirenr_ub1[sec]) {
520 secondaries = true;
521 break;
522 }
523 } // skip secondary hits
524
525 if (it_ub1 == 19 || secondaries)
526 break;
527 wirenr_ub1[it_ub1] = digit->GetWireNumber();
528 time_ub1[it_ub1] = times;
529 it_ub1++;
530 break;
531 case 4:
532
533 for (int sec = 0; sec < it_ya1 - 1; sec++) {
534 if (digit->GetWireNumber() == wirenr_ya1[sec]) {
535 secondaries = true;
536 break;
537 }
538 } // skip secondary hits
539 if (it_ya1 == 19 || secondaries)
540 break;
541
542 wirenr_ya1[it_ya1] = digit->GetWireNumber();
543 time_ya1[it_ya1] = times;
544 it_ya1++;
545 break;
546 case 5:
547 for (int sec = 0; sec < it_yb1 - 1; sec++) {
548 if (digit->GetWireNumber() == wirenr_yb1[sec]) {
549 secondaries = true;
550 break;
551 }
552 } // skip secondary hits
553 if (it_yb1 == 19 || secondaries)
554 break;
555 wirenr_yb1[it_yb1] = digit->GetWireNumber();
556 time_yb1[it_yb1] = times;
557 it_yb1++;
558 break;
559 case 6:
560 for (int sec = 0; sec < it_xa1 - 1; sec++) {
561 if (digit->GetWireNumber() == wirenr_xa1[sec]) {
562 secondaries = true;
563 break;
564 }
565 } // skip secondarx hits
566
567 if (it_xa1 == 19 || secondaries)
568 break;
569 wirenr_xa1[it_xa1] = digit->GetWireNumber();
570 time_xa1[it_xa1] = times;
571 it_xa1++;
572 break;
573 case 7:
574 for (int sec = 0; sec < it_xb1 - 1; sec++) {
575 if (digit->GetWireNumber() == wirenr_xb1[sec]) {
576 secondaries = true;
577 break;
578 }
579 } // skip secondary hits
580
581 if (it_xb1 == 19 || secondaries)
582 break;
583 wirenr_xb1[it_xb1] = digit->GetWireNumber();
584 time_xb1[it_xb1] = times;
585 it_xb1++;
586 break;
587 case 8:
588 for (int sec = 0; sec < it_va2 - 1; sec++) {
589 if (digit->GetWireNumber() == wirenr_va2[sec]) {
590 secondaries = true;
591 break;
592 }
593 } // skip secondary hits
594
595 if (it_va2 == 19 || secondaries)
596 break;
597 wirenr_va2[it_va2] = digit->GetWireNumber();
598 time_va2[it_va2] = times;
599 it_va2++;
600 break;
601 case 9:
602 for (int sec = 0; sec < it_vb2 - 1; sec++) {
603 if (digit->GetWireNumber() == wirenr_vb2[sec]) {
604 secondaries = true;
605 break;
606 }
607 } // skip secondary hits
608
609 if (it_vb2 == 19 || secondaries)
610 break;
611 wirenr_vb2[it_vb2] = digit->GetWireNumber();
612 time_vb2[it_vb2] = times;
613 it_vb2++;
614 break;
615 case 10:
616 for (int sec = 0; sec < it_ua2 - 1; sec++) {
617 if (digit->GetWireNumber() == wirenr_ua2[sec]) {
618 secondaries = true;
619 break;
620 }
621 } // skip secondary hits
622 if (it_ua2 == 19 || secondaries)
623 break;
624 wirenr_ua2[it_ua2] = digit->GetWireNumber();
625 time_ua2[it_ua2] = times;
626 it_ua2++;
627 break;
628 case 11:
629 for (int sec = 0; sec < it_ub2 - 1; sec++) {
630 if (digit->GetWireNumber() == wirenr_ub2[sec]) {
631 secondaries = true;
632 break;
633 }
634 } // skip secondary hits
635
636 if (it_ub2 == 19 || secondaries)
637 break;
638 wirenr_ub2[it_ub2] = digit->GetWireNumber();
639 time_ub2[it_ub2] = times;
640 it_ub2++;
641 break;
642 case 12:
643 for (int sec = 0; sec < it_ya2 - 1; sec++) {
644 if (digit->GetWireNumber() == wirenr_ya2[sec]) {
645 secondaries = true;
646 break;
647 }
648 } // skip secondary hits
649 if (digit->GetWireNumber() == 111)
650 break;
651 if (it_ya2 == 19 || secondaries)
652 break;
653 wirenr_ya2[it_ya2] = digit->GetWireNumber();
654 time_ya2[it_ya2] = times;
655 it_ya2++;
656 break;
657 case 13:
658 for (int sec = 0; sec < it_yb2 - 1; sec++) {
659 if (digit->GetWireNumber() == wirenr_yb2[sec]) {
660 secondaries = true;
661 break;
662 }
663 } // skip secondary hits
664
665 if (it_yb2 == 19 || secondaries)
666 break;
667 wirenr_yb2[it_yb2] = digit->GetWireNumber();
668 time_yb2[it_yb2] = times;
669 it_yb2++;
670 break;
671 case 14:
672
673 for (int sec = 0; sec < it_xa2 - 1; sec++) {
674 if (digit->GetWireNumber() == wirenr_xa2[sec]) {
675 secondaries = true;
676 break;
677 }
678 } // skip secondarx hits
679 if (it_xa2 == 19 || secondaries)
680 break;
681 wirenr_xa2[it_xa2] = digit->GetWireNumber();
682 time_xa2[it_xa2] = times;
683 it_xa2++;
684 break;
685 case 15:
686 for (int sec = 0; sec < it_xb2 - 1; sec++) {
687 if (digit->GetWireNumber() == wirenr_xb2[sec]) {
688 secondaries = true;
689 break;
690 }
691 } // skip secondarx hits
692 if (it_xb2 == 19 || secondaries)
693 break;
694 wirenr_xb2[it_xb2] = digit->GetWireNumber();
695 time_xb2[it_xb2] = times;
696 it_xb2++;
697 break;
698 }
699 } // for digis in event iDig
700
701 if (!goodEv)
702 continue;
703 float rh_segDC1[8][150];
704 float rh_segDC2[8][150];
705 float rh_sigm_segDC1[8][150];
706 float rh_sigm_segDC2[8][150];
707 float xDC1_glob[150];
708 float yDC1_glob[150];
709 float xDC2_glob[150];
710 float yDC2_glob[150];
711 float chi2_DC1[150];
712 float chi2_DC2[150];
713 int size_segDC1[150];
714 int size_segDC2[150];
715 int nDC1_segments = 0;
716 int nDC2_segments = 0;
717 bool has7DC1 = false;
718 bool has7DC2 = false;
719
720 float x1_ab[2][150];
721 float y1_ab[2][150];
722 float u1_ab[2][150];
723 float v1_ab[2][150];
724 float sigm_x1_ab[2][150];
725 float sigm_y1_ab[2][150];
726 float sigm_u1_ab[2][150];
727 float sigm_v1_ab[2][150];
728 float x2_ab[2][150];
729 float y2_ab[2][150];
730 float u2_ab[2][150];
731 float v2_ab[2][150];
732 float sigm_x2_ab[2][150];
733 float sigm_y2_ab[2][150];
734 float sigm_u2_ab[2][150];
735 float sigm_v2_ab[2][150];
736
737 // single hits on ab-plane
738 float x1_single[2][40];
739 float y1_single[2][40];
740 float u1_single[2][40];
741 float v1_single[2][40];
742 float sigm_x1_single[2][40];
743 float sigm_y1_single[2][40];
744 float sigm_u1_single[2][40];
745 float sigm_v1_single[2][40];
746 float x2_single[2][40];
747 float y2_single[2][40];
748 float u2_single[2][40];
749 float v2_single[2][40];
750 float sigm_x2_single[2][40];
751 float sigm_y2_single[2][40];
752 float sigm_u2_single[2][40];
753 float sigm_v2_single[2][40];
754
755 for (int i = 0; i < 8; i++) {
756 for (int j = 0; j < 150; j++) {
757 rh_segDC1[i][j] = -999;
758 rh_segDC2[i][j] = -999;
759 rh_sigm_segDC1[i][j] = 1;
760 rh_sigm_segDC2[i][j] = 1;
761 if (i > 1)
762 continue;
763 x1_ab[i][j] = -999;
764 y1_ab[i][j] = -999;
765 u1_ab[i][j] = -999;
766 v1_ab[i][j] = -999;
767 x2_ab[i][j] = -999;
768 y2_ab[i][j] = -999;
769 u2_ab[i][j] = -999;
770 v2_ab[i][j] = -999;
771 sigm_x1_ab[i][j] = 1;
772 sigm_y1_ab[i][j] = 1;
773 sigm_u1_ab[i][j] = 1;
774 sigm_v1_ab[i][j] = 1;
775 sigm_x2_ab[i][j] = 1;
776 sigm_y2_ab[i][j] = 1;
777 sigm_u2_ab[i][j] = 1;
778 sigm_v2_ab[i][j] = 1;
779
780 if (j > 39)
781 continue;
782 x1_single[i][j] = -999;
783 y1_single[i][j] = -999;
784 u1_single[i][j] = -999;
785 v1_single[i][j] = -999;
786 x2_single[i][j] = -999;
787 y2_single[i][j] = -999;
788 u2_single[i][j] = -999;
789 v2_single[i][j] = -999;
790 sigm_x1_single[i][j] = 1;
791 sigm_y1_single[i][j] = 1;
792 sigm_u1_single[i][j] = 1;
793 sigm_v1_single[i][j] = 1;
794 sigm_x2_single[i][j] = 1;
795 sigm_y2_single[i][j] = 1;
796 sigm_u2_single[i][j] = 1;
797 sigm_v2_single[i][j] = 1;
798
799 if (i > 0)
800 continue;
801
802 size_segDC1[j] = 0;
803 size_segDC2[j] = 0;
804 chi2_DC1[j] = 50;
805 chi2_DC2[j] = 50;
806 if (j > 39)
807 continue;
808
809 xDC1_glob[j] = -999;
810 yDC1_glob[j] = -999;
811 xDC2_glob[j] = -999;
812 yDC2_glob[j] = -999;
813 }
814 }
815 int pair_x1 = 0;
816 int pair_y1 = 0;
817 int pair_u1 = 0;
818 int pair_v1 = 0;
819 int single_xa1 = 0;
820 int single_ya1 = 0;
821 int single_ua1 = 0;
822 int single_va1 = 0;
823 int single_xb1 = 0;
824 int single_yb1 = 0;
825 int single_ub1 = 0;
826 int single_vb1 = 0;
827
828 float xa1_pm[2] = {-999};
829 float xb1_pm[2] = {-999};
830 float ya1_pm[2] = {-999};
831 float yb1_pm[2] = {-999};
832 float ua1_pm[2] = {-999};
833 float ub1_pm[2] = {-999};
834 float va1_pm[2] = {-999};
835 float vb1_pm[2] = {-999};
836
837 // --- X ---
838
839 for (Int_t i = 0; i < it_xa1; ++i) {
840 for (Int_t j = 0; j < it_xb1; ++j) {
841 if (pair_x1 > 48)
842 break;
843 if ((wirenr_xa1[i] != wirenr_xb1[j] && wirenr_xa1[i] != wirenr_xb1[j] + 1))
844 continue;
845 int func_nr_a = -1;
846 int func_nr_b = -1;
847 for (int t_it = 0; t_it < 4; t_it++) {
848 if (time_xa1[i] >= t_dc[t_it][0] && time_xa1[i] < t_dc[t_it + 1][0]) {
849 func_nr_a = t_it;
850 break;
851 }
852 }
853 double time_xa = time_xa1[i];
854 double d_a = 0;
855 double d_b = 0;
856 if (func_nr_a == 1 || func_nr_a == 2)
857 d_a =
858 scale[0]
859 * (pol_par_dc[1][0][0] + pol_par_dc[1][1][0] * time_xa + pol_par_dc[1][2][0] * time_xa * time_xa
860 + pol_par_dc[1][3][0] * time_xa * time_xa * time_xa
861 + pol_par_dc[1][4][0] * time_xa * time_xa * time_xa * time_xa);
862 if (func_nr_a == 0)
863 d_a = 0;
864 if (func_nr_a == 3)
865 d_a =
866 scale[0]
867 * (pol_par_dc[2][0][0] + pol_par_dc[2][1][0] * time_xa + pol_par_dc[2][2][0] * time_xa * time_xa
868 + pol_par_dc[2][3][0] * time_xa * time_xa * time_xa
869 + pol_par_dc[2][4][0] * time_xa * time_xa * time_xa * time_xa);
870
871 for (int t_it = 0; t_it < 4; t_it++) {
872 if (time_xb1[j] >= t_dc[t_it][1] && time_xb1[j] < t_dc[t_it + 1][1]) {
873 func_nr_b = t_it;
874 break;
875 }
876 }
877 time_xa = time_xb1[j];
878
879 if (func_nr_b == 1 || func_nr_b == 2)
880 d_b =
881 scale[1]
882 * (pol_par_dc[1][0][1] + pol_par_dc[1][1][1] * time_xa + pol_par_dc[1][2][1] * time_xa * time_xa
883 + pol_par_dc[1][3][1] * time_xa * time_xa * time_xa
884 + pol_par_dc[1][4][1] * time_xa * time_xa * time_xa * time_xa);
885 if (func_nr_b == 0)
886 d_b = 0;
887 if (func_nr_b == 3)
888 d_b =
889 scale[1]
890 * (pol_par_dc[2][0][1] + pol_par_dc[2][1][1] * time_xa + pol_par_dc[2][2][1] * time_xa * time_xa
891 + pol_par_dc[2][3][1] * time_xa * time_xa * time_xa
892 + pol_par_dc[2][4][1] * time_xa * time_xa * time_xa * time_xa);
893
894 xa1_pm[0] = wirenr_xa1[i] - 119 + d_a;
895 xa1_pm[1] = wirenr_xa1[i] - 119 - d_a;
896 xb1_pm[0] = wirenr_xb1[j] - 118.5 + d_b;
897 xb1_pm[1] = wirenr_xb1[j] - 118.5 - d_b;
898
899 double dmin1 = 999;
900 for (int k = 0; k < 2; k++) {
901 for (int m = 0; m < 2; m++) {
902 if (fabs(xa1_pm[k] - xb1_pm[m]) < dmin1) {
903 dmin1 = fabs(xa1_pm[k] - xb1_pm[m]);
904 x1_ab[0][pair_x1] = xa1_pm[k];
905 x1_ab[1][pair_x1] = xb1_pm[m];
906 }
907 }
908 }
909 if (d_a < 0.02) {
910 sigm_x1_ab[0][pair_x1] = 0.08 * 0.08;
911 }
912 if (d_a >= 0.02 && d_a < 0.1) {
913 sigm_x1_ab[0][pair_x1] = 0.06 * 0.06;
914 }
915 if (d_a >= 0.1 && d_a < 0.4) {
916 sigm_x1_ab[0][pair_x1] = 0.025 * 0.025;
917 }
918 if (d_a >= 0.4 && d_a < 0.41) {
919 sigm_x1_ab[0][pair_x1] = 0.08 * 0.08;
920 }
921 if (d_a >= 0.41) {
922 sigm_x1_ab[0][pair_x1] = 0.10 * 0.10;
923 }
924
925 if (d_b < 0.02) {
926 sigm_x1_ab[1][pair_x1] = 0.08 * 0.08;
927 }
928 if (d_b >= 0.02 && d_b < 0.1) {
929 sigm_x1_ab[1][pair_x1] = 0.06 * 0.06;
930 }
931 if (d_b >= 0.1 && d_b < 0.4) {
932 sigm_x1_ab[1][pair_x1] = 0.025 * 0.025;
933 }
934 if (d_b >= 0.4 && d_b < 0.41) {
935 sigm_x1_ab[1][pair_x1] = 0.08 * 0.08;
936 }
937 if (d_b >= 0.41) {
938 sigm_x1_ab[1][pair_x1] = 0.10 * 0.10;
939 }
940 pair_x1++;
941
942 used_xa1[i] = true;
943 used_xb1[j] = true;
944
945 } // j for2
946 } // i for1 X.
947
948 // reconstruct single X-plane hits
949
950 for (Int_t i = 0; i < it_xa1; ++i) {
951 if (used_xa1[i])
952 continue;
953
954 int func_nr_a = -1;
955 for (int t_it = 0; t_it < 4; t_it++) {
956 if (time_xa1[i] >= t_dc[t_it][0] && time_xa1[i] < t_dc[t_it + 1][0]) {
957 func_nr_a = t_it;
958 break;
959 }
960 }
961 double time_xa = time_xa1[i];
962 double d_a = 0;
963 double d_b = 0;
964
965 if (func_nr_a == 1 || func_nr_a == 2)
966 d_a = scale[0]
967 * (pol_par_dc[1][0][0] + pol_par_dc[1][1][0] * time_xa + pol_par_dc[1][2][0] * time_xa * time_xa
968 + pol_par_dc[1][3][0] * time_xa * time_xa * time_xa
969 + pol_par_dc[1][4][0] * time_xa * time_xa * time_xa * time_xa);
970 if (func_nr_a == 0)
971 d_a = 0;
972 if (func_nr_a == 3)
973 d_a = scale[0]
974 * (pol_par_dc[2][0][0] + pol_par_dc[2][1][0] * time_xa + pol_par_dc[2][2][0] * time_xa * time_xa
975 + pol_par_dc[2][3][0] * time_xa * time_xa * time_xa
976 + pol_par_dc[2][4][0] * time_xa * time_xa * time_xa * time_xa);
977
978 x1_single[0][single_xa1] = wirenr_xa1[i] - 119 + d_a;
979 x1_single[0][single_xa1 + 1] = wirenr_xa1[i] - 119 - d_a;
980
981 if (d_a < 0.02) {
982 sigm_x1_single[0][single_xa1] = 0.08 * 0.08;
983 sigm_x1_single[0][single_xa1 + 1] = 0.08 * 0.08;
984 }
985 if (d_a >= 0.02 && d_a < 0.1) {
986 sigm_x1_single[0][single_xa1] = 0.06 * 0.06;
987 sigm_x1_single[0][single_xa1 + 1] = 0.06 * 0.06;
988 }
989 if (d_a >= 0.1 && d_a < 0.4) {
990 sigm_x1_single[0][single_xa1] = 0.025 * 0.025;
991 sigm_x1_single[0][single_xa1 + 1] = 0.025 * 0.025;
992 }
993 if (d_a >= 0.4 && d_a < 0.41) {
994 sigm_x1_single[0][single_xa1] = 0.08 * 0.08;
995 sigm_x1_single[0][single_xa1 + 1] = 0.08 * 0.08;
996 }
997 if (d_a >= 0.41) {
998 sigm_x1_single[0][single_xa1] = 0.10 * 0.10;
999 sigm_x1_single[0][single_xa1 + 1] = 0.10 * 0.10;
1000 }
1001
1002 single_xa1 += 2;
1003 } // for single xa
1004
1005 for (Int_t j = 0; j < it_xb1; ++j) {
1006 if (used_xb1[j])
1007 continue;
1008 int func_nr_b = -1;
1009
1010 double time_xa = time_xb1[j];
1011 double d_a = 0;
1012 double d_b = 0;
1013 for (int t_it = 0; t_it < 4; t_it++) {
1014 if (time_xb1[j] >= t_dc[t_it][1] && time_xb1[j] < t_dc[t_it + 1][1]) {
1015 func_nr_b = t_it;
1016 break;
1017 }
1018 }
1019 time_xa = time_xb1[j];
1020 if (func_nr_b == 1 || func_nr_b == 2)
1021 d_b = scale[1]
1022 * (pol_par_dc[1][0][1] + pol_par_dc[1][1][1] * time_xa + pol_par_dc[1][2][1] * time_xa * time_xa
1023 + pol_par_dc[1][3][1] * time_xa * time_xa * time_xa
1024 + pol_par_dc[1][4][1] * time_xa * time_xa * time_xa * time_xa);
1025 if (func_nr_b == 0)
1026 d_b = 0;
1027 if (func_nr_b == 3)
1028 d_b = scale[1]
1029 * (pol_par_dc[2][0][1] + pol_par_dc[2][1][1] * time_xa + pol_par_dc[2][2][1] * time_xa * time_xa
1030 + pol_par_dc[2][3][1] * time_xa * time_xa * time_xa
1031 + pol_par_dc[2][4][1] * time_xa * time_xa * time_xa * time_xa);
1032 x1_single[1][single_xb1] = wirenr_xb1[j] - 118.5 + d_b;
1033 x1_single[1][single_xb1 + 1] = wirenr_xb1[j] - 118.5 - d_b;
1034
1035 if (d_b < 0.02) {
1036 sigm_x1_single[1][single_xb1] = 0.08 * 0.08;
1037 sigm_x1_single[1][single_xb1 + 1] = 0.08 * 0.08;
1038 }
1039 if (d_b >= 0.02 && d_b < 0.1) {
1040 sigm_x1_single[1][single_xb1] = 0.06 * 0.06;
1041 sigm_x1_single[1][single_xb1 + 1] = 0.06 * 0.06;
1042 }
1043 if (d_b >= 0.1 && d_b < 0.4) {
1044 sigm_x1_single[1][single_xb1] = 0.0250 * 0.0250;
1045 sigm_x1_single[1][single_xb1 + 1] = 0.0250 * 0.0250;
1046 }
1047 if (d_b >= 0.4 && d_b < 0.41) {
1048 sigm_x1_single[1][single_xb1] = 0.08 * 0.08;
1049 sigm_x1_single[1][single_xb1 + 1] = 0.08 * 0.08;
1050 }
1051 if (d_b >= 0.41) {
1052 sigm_x1_single[1][single_xb1] = 0.10 * 0.10;
1053 sigm_x1_single[1][single_xb1 + 1] = 0.10 * 0.10;
1054 }
1055
1056 single_xb1 += 2;
1057
1058 } // i for1 X.
1059
1060 // ---- Y ----
1061
1062 for (Int_t i = 0; i < it_ya1; ++i) {
1063 for (Int_t j = 0; j < it_yb1; ++j) {
1064 if (pair_y1 > 48)
1065 break;
1066 if ((wirenr_ya1[i] != wirenr_yb1[j] && wirenr_ya1[i] != wirenr_yb1[j] + 1))
1067 continue;
1068 int func_nr_a = -1;
1069 int func_nr_b = -1;
1070 for (int t_it = 0; t_it < 4; t_it++) {
1071 if (time_ya1[i] >= t_dc[t_it][2] && time_ya1[i] < t_dc[t_it + 1][2]) {
1072 func_nr_a = t_it;
1073 break;
1074 }
1075 }
1076 double d_a = 0;
1077 double d_b = 0;
1078
1079 if (func_nr_a == 1 || func_nr_a == 2)
1080 d_a = scale[2]
1081 * (pol_par_dc[1][0][2] + pol_par_dc[1][1][2] * time_ya1[i]
1082 + pol_par_dc[1][2][2] * time_ya1[i] * time_ya1[i]
1083 + pol_par_dc[1][3][2] * time_ya1[i] * time_ya1[i] * time_ya1[i]
1084 + pol_par_dc[1][4][2] * time_ya1[i] * time_ya1[i] * time_ya1[i] * time_ya1[i]);
1085 if (func_nr_a == 0)
1086 d_a = 0;
1087 if (func_nr_a == 3)
1088 d_a = scale[2]
1089 * (pol_par_dc[2][0][2] + pol_par_dc[2][1][2] * time_ya1[i]
1090 + pol_par_dc[2][2][2] * time_ya1[i] * time_ya1[i]
1091 + pol_par_dc[2][3][2] * time_ya1[i] * time_ya1[i] * time_ya1[i]
1092 + pol_par_dc[2][4][2] * time_ya1[i] * time_ya1[i] * time_ya1[i] * time_ya1[i]);
1093 for (int t_it = 0; t_it < 4; t_it++) {
1094 if (time_yb1[j] >= t_dc[t_it][3] && time_yb1[j] < t_dc[t_it + 1][3]) {
1095 func_nr_b = t_it;
1096 break;
1097 }
1098 }
1099 if (func_nr_b == 1 || func_nr_b == 2)
1100 d_b = scale[3]
1101 * (pol_par_dc[1][0][3] + pol_par_dc[1][1][3] * time_yb1[j]
1102 + pol_par_dc[1][2][3] * time_yb1[j] * time_yb1[j]
1103 + pol_par_dc[1][3][3] * time_yb1[j] * time_yb1[j] * time_yb1[j]
1104 + pol_par_dc[1][4][3] * time_yb1[j] * time_yb1[j] * time_yb1[j] * time_yb1[j]);
1105 if (func_nr_b == 0)
1106 d_b = 0;
1107 if (func_nr_b == 3)
1108 d_b = scale[3]
1109 * (pol_par_dc[2][0][3] + pol_par_dc[2][1][3] * time_yb1[j]
1110 + pol_par_dc[2][2][3] * time_yb1[j] * time_yb1[j]
1111 + pol_par_dc[2][3][3] * time_yb1[j] * time_yb1[j] * time_yb1[j]
1112 + pol_par_dc[2][4][3] * time_yb1[j] * time_yb1[j] * time_yb1[j] * time_yb1[j]);
1113
1114 ya1_pm[0] = wirenr_ya1[i] - 119 + d_a;
1115 ya1_pm[1] = wirenr_ya1[i] - 119 - d_a;
1116 yb1_pm[0] = wirenr_yb1[j] - 118.5 + d_b;
1117 yb1_pm[1] = wirenr_yb1[j] - 118.5 - d_b;
1118
1119 double dmin1 = 999;
1120
1121 for (int k = 0; k < 2; k++) {
1122 for (int m = 0; m < 2; m++) {
1123 if (fabs(ya1_pm[k] - yb1_pm[m]) < dmin1) {
1124 dmin1 = fabs(ya1_pm[k] - yb1_pm[m]);
1125 y1_ab[0][pair_y1] = ya1_pm[k];
1126 y1_ab[1][pair_y1] = yb1_pm[m];
1127 }
1128 }
1129 }
1130 if (d_a < 0.02) {
1131 sigm_y1_ab[0][pair_y1] = 0.08 * 0.08;
1132 }
1133 if (d_a >= 0.02 && d_a < 0.1) {
1134 sigm_y1_ab[0][pair_y1] = 0.06 * 0.06;
1135 }
1136 if (d_a >= 0.1 && d_a < 0.4) {
1137 sigm_y1_ab[0][pair_y1] = 0.025 * 0.025;
1138 }
1139 if (d_a >= 0.4 && d_a < 0.41) {
1140 sigm_y1_ab[0][pair_y1] = 0.08 * 0.08;
1141 }
1142 if (d_a >= 0.41) {
1143 sigm_y1_ab[0][pair_y1] = 0.10 * 0.10;
1144 }
1145
1146 if (d_b < 0.02) {
1147 sigm_y1_ab[1][pair_y1] = 0.08 * 0.08;
1148 }
1149 if (d_b >= 0.02 && d_b < 0.1) {
1150 sigm_y1_ab[1][pair_y1] = 0.06 * 0.06;
1151 }
1152 if (d_b >= 0.1 && d_b < 0.4) {
1153 sigm_y1_ab[1][pair_y1] = 0.025 * 0.025;
1154 }
1155 if (d_b >= 0.4 && d_b < 0.41) {
1156 sigm_y1_ab[1][pair_y1] = 0.08 * 0.08;
1157 }
1158 if (d_b >= 0.41) {
1159 sigm_y1_ab[1][pair_y1] = 0.10 * 0.10;
1160 }
1161 pair_y1++;
1162
1163 used_ya1[i] = true;
1164 used_yb1[j] = true;
1165
1166 } // j for2
1167 } // i for1 Y.
1168
1169 // reconstruct single Y-plane hits
1170
1171 for (Int_t i = 0; i < it_ya1; ++i) {
1172 if (used_ya1[i])
1173 continue;
1174 double d_a = 0;
1175 double d_b = 0;
1176 int func_nr_a = -1;
1177 for (int t_it = 0; t_it < 4; t_it++) {
1178 if (time_ya1[i] >= t_dc[t_it][2] && time_ya1[i] < t_dc[t_it + 1][2]) {
1179 func_nr_a = t_it;
1180 break;
1181 }
1182 }
1183 if (func_nr_a == 1 || func_nr_a == 2)
1184 d_a = scale[2]
1185 * (pol_par_dc[1][0][2] + pol_par_dc[1][1][2] * time_ya1[i]
1186 + pol_par_dc[1][2][2] * time_ya1[i] * time_ya1[i]
1187 + pol_par_dc[1][3][2] * time_ya1[i] * time_ya1[i] * time_ya1[i]
1188 + pol_par_dc[1][4][2] * time_ya1[i] * time_ya1[i] * time_ya1[i] * time_ya1[i]);
1189 if (func_nr_a == 0)
1190 d_a = 0;
1191 if (func_nr_a == 3)
1192 d_a = scale[2]
1193 * (pol_par_dc[2][0][2] + pol_par_dc[2][1][2] * time_ya1[i]
1194 + pol_par_dc[2][2][2] * time_ya1[i] * time_ya1[i]
1195 + pol_par_dc[2][3][2] * time_ya1[i] * time_ya1[i] * time_ya1[i]
1196 + pol_par_dc[2][4][2] * time_ya1[i] * time_ya1[i] * time_ya1[i] * time_ya1[i]);
1197
1198 y1_single[0][single_ya1] = wirenr_ya1[i] - 119 + d_a;
1199 y1_single[0][single_ya1 + 1] = wirenr_ya1[i] - 119 - d_a;
1200
1201 if (d_a < 0.02) {
1202 sigm_y1_single[0][single_ya1] = 0.08 * 0.08;
1203 sigm_y1_single[0][single_ya1 + 1] = 0.08 * 0.08;
1204 }
1205 if (d_a >= 0.02 && d_a < 0.1) {
1206 sigm_y1_single[0][single_ya1] = 0.06 * 0.06;
1207 sigm_y1_single[0][single_ya1 + 1] = 0.06 * 0.06;
1208 }
1209 if (d_a >= 0.1 && d_a < 0.4) {
1210 sigm_y1_single[0][single_ya1] = 0.025 * 0.025;
1211 sigm_y1_single[0][single_ya1 + 1] = 0.025 * 0.025;
1212 }
1213 if (d_a >= 0.4 && d_a < 0.41) {
1214 sigm_y1_single[0][single_ya1] = 0.08 * 0.08;
1215 sigm_y1_single[0][single_ya1 + 1] = 0.08 * 0.08;
1216 }
1217 if (d_a >= 0.41) {
1218 sigm_y1_single[0][single_ya1] = 0.10 * 0.10;
1219 sigm_y1_single[0][single_ya1 + 1] = 0.10 * 0.10;
1220 }
1221 single_ya1 += 2;
1222 } // for single ya
1223
1224 for (Int_t j = 0; j < it_yb1; ++j) {
1225 if (used_yb1[j])
1226 continue;
1227 double time_ya = time_yb1[j];
1228 int func_nr_b = -1;
1229 for (int t_it = 0; t_it < 4; t_it++) {
1230 if (time_yb1[j] >= t_dc[t_it][3] && time_yb1[j] < t_dc[t_it + 1][3]) {
1231 func_nr_b = t_it;
1232 break;
1233 }
1234 }
1235 double d_a = 0;
1236 double d_b = 0;
1237 if (func_nr_b == 1 || func_nr_b == 2)
1238 d_b = scale[3]
1239 * (pol_par_dc[1][0][3] + pol_par_dc[1][1][3] * time_yb1[j]
1240 + pol_par_dc[1][2][3] * time_yb1[j] * time_yb1[j]
1241 + pol_par_dc[1][3][3] * time_yb1[j] * time_yb1[j] * time_yb1[j]
1242 + pol_par_dc[1][4][3] * time_yb1[j] * time_yb1[j] * time_yb1[j] * time_yb1[j]);
1243 if (func_nr_b == 0)
1244 d_b = 0;
1245 if (func_nr_b == 3)
1246 d_b = scale[3]
1247 * (pol_par_dc[2][0][3] + pol_par_dc[2][1][3] * time_yb1[j]
1248 + pol_par_dc[2][2][3] * time_yb1[j] * time_yb1[j]
1249 + pol_par_dc[2][3][3] * time_yb1[j] * time_yb1[j] * time_yb1[j]
1250 + pol_par_dc[2][4][3] * time_yb1[j] * time_yb1[j] * time_yb1[j] * time_yb1[j]);
1251 y1_single[1][single_yb1] = wirenr_yb1[j] - 118.5 + d_b;
1252 y1_single[1][single_yb1 + 1] = wirenr_yb1[j] - 118.5 - d_b;
1253
1254 if (d_b < 0.02) {
1255 sigm_y1_single[1][single_yb1] = 0.08 * 0.08;
1256 sigm_y1_single[1][single_yb1 + 1] = 0.08 * 0.08;
1257 }
1258 if (d_b >= 0.02 && d_b < 0.1) {
1259 sigm_y1_single[1][single_yb1] = 0.06 * 0.06;
1260 sigm_y1_single[1][single_yb1 + 1] = 0.06 * 0.06;
1261 }
1262 if (d_b >= 0.1 && d_b < 0.4) {
1263 sigm_y1_single[1][single_yb1] = 0.025 * 0.025;
1264 sigm_y1_single[1][single_yb1 + 1] = 0.025 * 0.025;
1265 }
1266 if (d_b >= 0.4 && d_b < 0.41) {
1267 sigm_y1_single[1][single_yb1] = 0.10 * 0.10;
1268 sigm_y1_single[1][single_yb1 + 1] = 0.10 * 0.10;
1269 }
1270 if (d_b >= 0.41) {
1271 sigm_y1_single[1][single_yb1] = 0.10 * 0.10;
1272 sigm_y1_single[1][single_yb1 + 1] = 0.10 * 0.10;
1273 }
1274 single_yb1 += 2;
1275
1276 } // i for1 Y.
1277
1278 // ---- U ---
1279
1280 for (Int_t i = 0; i < it_ua1; ++i) {
1281 for (Int_t j = 0; j < it_ub1; ++j) {
1282 if (pair_u1 > 48)
1283 break;
1284 if ((wirenr_ua1[i] != wirenr_ub1[j] && wirenr_ua1[i] != wirenr_ub1[j] + 1))
1285 continue;
1286 double d_a = 0;
1287 double d_b = 0;
1288
1289 int func_nr_a = -1;
1290 int func_nr_b = -1;
1291 for (int t_it = 0; t_it < 4; t_it++) {
1292 if (time_ua1[i] >= t_dc[t_it][4] && time_ua1[i] < t_dc[t_it + 1][4]) {
1293 func_nr_a = t_it;
1294 break;
1295 }
1296 }
1297 if (func_nr_a == 0)
1298 d_a = 0;
1299 if (func_nr_a == 3)
1300 d_a = scale[4]
1301 * (pol_par_dc[2][0][4] + pol_par_dc[2][1][4] * time_ua1[i]
1302 + pol_par_dc[2][2][4] * time_ua1[i] * time_ua1[i]
1303 + pol_par_dc[2][3][4] * time_ua1[i] * time_ua1[i] * time_ua1[i]
1304 + pol_par_dc[2][4][4] * time_ua1[i] * time_ua1[i] * time_ua1[i] * time_ua1[i]);
1305 if (func_nr_a == 1 || func_nr_a == 2)
1306 d_a = scale[4]
1307 * (pol_par_dc[1][0][4] + pol_par_dc[1][1][4] * time_ua1[i]
1308 + pol_par_dc[1][2][4] * time_ua1[i] * time_ua1[i]
1309 + pol_par_dc[1][3][4] * time_ua1[i] * time_ua1[i] * time_ua1[i]
1310 + pol_par_dc[1][4][4] * time_ua1[i] * time_ua1[i] * time_ua1[i] * time_ua1[i]);
1311 for (int t_it = 0; t_it < 4; t_it++) {
1312 if (time_ub1[j] >= t_dc[t_it][5] && time_ub1[j] < t_dc[t_it + 1][5]) {
1313 func_nr_b = t_it;
1314 break;
1315 }
1316 }
1317
1318 if (func_nr_b == 0)
1319 d_b = 0;
1320 if (func_nr_b == 3)
1321 d_b = scale[5]
1322 * (pol_par_dc[2][0][5] + pol_par_dc[2][1][5] * time_ub1[j]
1323 + pol_par_dc[2][2][5] * time_ub1[j] * time_ub1[j]
1324 + pol_par_dc[2][3][5] * time_ub1[j] * time_ub1[j] * time_ub1[j]
1325 + pol_par_dc[2][4][5] * time_ub1[j] * time_ub1[j] * time_ub1[j] * time_ub1[j]);
1326 if (func_nr_b == 1 || func_nr_b == 2)
1327 d_b = scale[5]
1328 * (pol_par_dc[1][0][5] + pol_par_dc[1][1][5] * time_ub1[j]
1329 + pol_par_dc[1][2][5] * time_ub1[j] * time_ub1[j]
1330 + pol_par_dc[1][3][5] * time_ub1[j] * time_ub1[j] * time_ub1[j]
1331 + pol_par_dc[1][4][5] * time_ub1[j] * time_ub1[j] * time_ub1[j] * time_ub1[j]);
1332 ua1_pm[0] = wirenr_ua1[i] - 119 + d_a;
1333 ua1_pm[1] = wirenr_ua1[i] - 119 - d_a;
1334 ub1_pm[0] = wirenr_ub1[j] - 118.5 + d_b;
1335 ub1_pm[1] = wirenr_ub1[j] - 118.5 - d_b;
1336 double dmin1 = 999;
1337
1338 for (int k = 0; k < 2; k++) {
1339 for (int m = 0; m < 2; m++) {
1340 if (fabs(ua1_pm[k] - ub1_pm[m]) < dmin1) {
1341 dmin1 = fabs(ua1_pm[k] - ub1_pm[m]);
1342 u1_ab[0][pair_u1] = ua1_pm[k];
1343 u1_ab[1][pair_u1] = ub1_pm[m];
1344 }
1345 }
1346 }
1347 if (d_a < 0.02) {
1348 sigm_u1_ab[0][pair_u1] = 0.08 * 0.08;
1349 }
1350 if (d_a >= 0.02 && d_a < 0.1) {
1351 sigm_u1_ab[0][pair_u1] = 0.06 * 0.06;
1352 }
1353 if (d_a >= 0.1 && d_a < 0.4) {
1354 sigm_u1_ab[0][pair_u1] = 0.025 * 0.025;
1355 }
1356 if (d_a >= 0.4 && d_a < 0.41) {
1357 sigm_u1_ab[0][pair_u1] = 0.08 * 0.08;
1358 }
1359 if (d_a >= 0.41) {
1360 sigm_u1_ab[0][pair_u1] = 0.10 * 0.10;
1361 }
1362
1363 if (d_b < 0.02) {
1364 sigm_u1_ab[1][pair_u1] = 0.08 * 0.08;
1365 }
1366 if (d_b >= 0.02 && d_b < 0.1) {
1367 sigm_u1_ab[1][pair_u1] = 0.06 * 0.06;
1368 }
1369 if (d_b >= 0.1 && d_b < 0.4) {
1370 sigm_u1_ab[1][pair_u1] = 0.025 * 0.025;
1371 }
1372 if (d_b >= 0.4 && d_b < 0.41) {
1373 sigm_u1_ab[1][pair_u1] = 0.08 * 0.08;
1374 }
1375 if (d_b >= 0.41) {
1376 sigm_u1_ab[1][pair_u1] = 0.10 * 0.10;
1377 }
1378 pair_u1++;
1379 used_ua1[i] = true;
1380 used_ub1[j] = true;
1381 } // j for2
1382 } // i for1 U.
1383
1384 // reconstruct single U-plane hits
1385
1386 for (Int_t i = 0; i < it_ua1; ++i) {
1387 if (used_ua1[i])
1388 continue;
1389
1390 double d_a = 0;
1391 double d_b = 0;
1392 int func_nr_a = -1;
1393 for (int t_it = 0; t_it < 4; t_it++) {
1394 if (time_ua1[i] >= t_dc[t_it][4] && time_ua1[i] < t_dc[t_it + 1][4]) {
1395 func_nr_a = t_it;
1396 break;
1397 }
1398 }
1399 if (func_nr_a == 0)
1400 d_a = 0;
1401 if (func_nr_a == 3)
1402 d_a = scale[4]
1403 * (pol_par_dc[2][0][4] + pol_par_dc[2][1][4] * time_ua1[i]
1404 + pol_par_dc[2][2][4] * time_ua1[i] * time_ua1[i]
1405 + pol_par_dc[2][3][4] * time_ua1[i] * time_ua1[i] * time_ua1[i]
1406 + pol_par_dc[2][4][4] * time_ua1[i] * time_ua1[i] * time_ua1[i] * time_ua1[i]);
1407 if (func_nr_a == 1 || func_nr_a == 2)
1408 d_a = scale[4]
1409 * (pol_par_dc[1][0][4] + pol_par_dc[1][1][4] * time_ua1[i]
1410 + pol_par_dc[1][2][4] * time_ua1[i] * time_ua1[i]
1411 + pol_par_dc[1][3][4] * time_ua1[i] * time_ua1[i] * time_ua1[i]
1412 + pol_par_dc[1][4][4] * time_ua1[i] * time_ua1[i] * time_ua1[i] * time_ua1[i]);
1413
1414 u1_single[0][single_ua1] = (wirenr_ua1[i] - 119 + d_a);
1415 u1_single[0][single_ua1 + 1] = (wirenr_ua1[i] - 119 - d_a);
1416 if (d_a < 0.02) {
1417 sigm_u1_single[0][single_ua1] = 0.08 * 0.08;
1418 sigm_u1_single[0][single_ua1 + 1] = 0.08 * 0.08;
1419 }
1420 if (d_a >= 0.02 && d_a < 0.1) {
1421 sigm_u1_single[0][single_ua1] = 0.06 * 0.06;
1422 sigm_u1_single[0][single_ua1 + 1] = 0.06 * 0.06;
1423 }
1424 if (d_a >= 0.1 && d_a < 0.4) {
1425 sigm_u1_single[0][single_ua1] = 0.025 * 0.025;
1426 sigm_u1_single[0][single_ua1 + 1] = 0.025 * 0.025;
1427 }
1428 if (d_a >= 0.4 && d_a < 0.41) {
1429 sigm_u1_single[0][single_ua1] = 0.08 * 0.08;
1430 sigm_u1_single[0][single_ua1 + 1] = 0.08 * 0.08;
1431 }
1432 if (d_a >= 0.41) {
1433 sigm_u1_single[0][single_ua1] = 0.10 * 0.10;
1434 sigm_u1_single[0][single_ua1 + 1] = 0.10 * 0.10;
1435 }
1436 single_ua1 += 2;
1437 } // for single ua
1438
1439 for (Int_t j = 0; j < it_ub1; ++j) {
1440 if (used_ub1[j])
1441 continue;
1442 int func_nr_b = -1;
1443 double d_a = 0;
1444 double d_b = 0;
1445 for (int t_it = 0; t_it < 4; t_it++) {
1446 if (time_ub1[j] >= t_dc[t_it][5] && time_ub1[j] < t_dc[t_it + 1][5]) {
1447 func_nr_b = t_it;
1448 break;
1449 }
1450 }
1451 if (func_nr_b == 0)
1452 d_b = 0;
1453 if (func_nr_b == 3)
1454 d_b = scale[5]
1455 * (pol_par_dc[2][0][5] + pol_par_dc[2][1][5] * time_ub1[j]
1456 + pol_par_dc[2][2][5] * time_ub1[j] * time_ub1[j]
1457 + pol_par_dc[2][3][5] * time_ub1[j] * time_ub1[j] * time_ub1[j]
1458 + pol_par_dc[2][4][5] * time_ub1[j] * time_ub1[j] * time_ub1[j] * time_ub1[j]);
1459 if (func_nr_b == 1 || func_nr_b == 2)
1460 d_b = scale[5]
1461 * (pol_par_dc[1][0][5] + pol_par_dc[1][1][5] * time_ub1[j]
1462 + pol_par_dc[1][2][5] * time_ub1[j] * time_ub1[j]
1463 + pol_par_dc[1][3][5] * time_ub1[j] * time_ub1[j] * time_ub1[j]
1464 + pol_par_dc[1][4][5] * time_ub1[j] * time_ub1[j] * time_ub1[j] * time_ub1[j]);
1465 u1_single[1][single_ub1] = (wirenr_ub1[j] - 118.5 + d_b);
1466 u1_single[1][single_ub1 + 1] = (wirenr_ub1[j] - 118.5 - d_b);
1467
1468 // hdev_occup[1]->Fill(d_b);
1469 if (d_b < 0.02) {
1470 sigm_u1_single[1][single_ub1] = 0.08 * 0.08;
1471 sigm_u1_single[1][single_ub1 + 1] = 0.08 * 0.08;
1472 }
1473 if (d_b >= 0.02 && d_b < 0.1) {
1474 sigm_u1_single[1][single_ub1] = 0.06 * 0.06;
1475 sigm_u1_single[1][single_ub1 + 1] = 0.06 * 0.06;
1476 }
1477 if (d_b >= 0.1 && d_b < 0.4) {
1478 sigm_u1_single[1][single_ub1] = 0.025 * 0.025;
1479 sigm_u1_single[1][single_ub1 + 1] = 0.025 * 0.025;
1480 }
1481 if (d_b >= 0.4 && d_b < 0.41) {
1482 sigm_u1_single[1][single_ub1] = 0.08 * 0.08;
1483 sigm_u1_single[1][single_ub1 + 1] = 0.08 * 0.08;
1484 }
1485 if (d_b >= 0.41) {
1486 sigm_u1_single[1][single_ub1] = 0.10 * 0.10;
1487 sigm_u1_single[1][single_ub1 + 1] = 0.10 * 0.10;
1488 }
1489 single_ub1 += 2;
1490 } // i for1 U.
1491
1492 // ---- V ---
1493
1494 for (Int_t i = 0; i < it_va1; ++i) {
1495 for (Int_t j = 0; j < it_vb1; ++j) {
1496
1497 if (pair_v1 > 48)
1498 break;
1499 if ((wirenr_va1[i] != wirenr_vb1[j] && wirenr_va1[i] != wirenr_vb1[j] + 1))
1500 continue;
1501
1502 double d_a = 0;
1503 double d_b = 0;
1504
1505 int func_nr_a = -1;
1506 int func_nr_b = -1;
1507 for (int t_it = 0; t_it < 4; t_it++) {
1508 if (time_va1[i] >= t_dc[t_it][6] && time_va1[i] < t_dc[t_it + 1][6]) {
1509 func_nr_a = t_it;
1510 break;
1511 }
1512 }
1513
1514 if (func_nr_a == 0)
1515 d_a = 0;
1516 if (func_nr_a == 3)
1517 d_a = scale[6]
1518 * (pol_par_dc[2][0][6] + pol_par_dc[2][1][6] * time_va1[i]
1519 + pol_par_dc[2][2][6] * time_va1[i] * time_va1[i]
1520 + pol_par_dc[2][3][6] * time_va1[i] * time_va1[i] * time_va1[i]
1521 + pol_par_dc[2][4][6] * time_va1[i] * time_va1[i] * time_va1[i] * time_va1[i]);
1522 if (func_nr_a == 1 || func_nr_a == 2)
1523 d_a = scale[6]
1524 * (pol_par_dc[1][0][6] + pol_par_dc[1][1][6] * time_va1[i]
1525 + pol_par_dc[1][2][6] * time_va1[i] * time_va1[i]
1526 + pol_par_dc[1][3][6] * time_va1[i] * time_va1[i] * time_va1[i]
1527 + pol_par_dc[1][4][6] * time_va1[i] * time_va1[i] * time_va1[i] * time_va1[i]);
1528 for (int t_it = 0; t_it < 4; t_it++) {
1529 if (time_vb1[j] >= t_dc[t_it][7] && time_vb1[j] < t_dc[t_it + 1][7]) {
1530 func_nr_b = t_it;
1531 break;
1532 }
1533 }
1534
1535 if (func_nr_b == 1 || func_nr_b == 2)
1536 d_b = scale[7]
1537 * (pol_par_dc[1][0][7] + pol_par_dc[1][1][7] * time_vb1[j]
1538 + pol_par_dc[1][2][7] * time_vb1[j] * time_vb1[j]
1539 + pol_par_dc[1][3][7] * time_vb1[j] * time_vb1[j] * time_vb1[j]
1540 + pol_par_dc[1][4][7] * time_vb1[j] * time_vb1[j] * time_vb1[j] * time_vb1[j]);
1541
1542 if (func_nr_b == 0)
1543 d_b = 0;
1544 if (func_nr_b == 3)
1545 d_b = scale[7]
1546 * (pol_par_dc[2][0][7] + pol_par_dc[2][1][7] * time_vb1[j]
1547 + pol_par_dc[2][2][7] * time_vb1[j] * time_vb1[j]
1548 + pol_par_dc[2][3][7] * time_vb1[j] * time_vb1[j] * time_vb1[j]
1549 + pol_par_dc[2][4][7] * time_vb1[j] * time_vb1[j] * time_vb1[j] * time_vb1[j]);
1550
1551 va1_pm[0] = wirenr_va1[i] - 119 + d_a;
1552 va1_pm[1] = wirenr_va1[i] - 119 - d_a;
1553 vb1_pm[0] = wirenr_vb1[j] - 118.5 + d_b;
1554 vb1_pm[1] = wirenr_vb1[j] - 118.5 - d_b;
1555 double dmin1 = 999;
1556
1557 for (int k = 0; k < 2; k++) {
1558 for (int m = 0; m < 2; m++) {
1559 if (fabs(va1_pm[k] - vb1_pm[m]) < dmin1) {
1560 dmin1 = fabs(va1_pm[k] - vb1_pm[m]);
1561 v1_ab[0][pair_v1] = va1_pm[k];
1562 v1_ab[1][pair_v1] = vb1_pm[m];
1563 }
1564 }
1565 }
1566 if (d_a < 0.02) {
1567 sigm_v1_ab[0][pair_v1] = 0.08 * 0.08;
1568 }
1569 if (d_a >= 0.02 && d_a < 0.1) {
1570 sigm_v1_ab[0][pair_v1] = 0.06 * 0.06;
1571 }
1572 if (d_a >= 0.1 && d_a < 0.4) {
1573 sigm_v1_ab[0][pair_v1] = 0.025 * 0.025;
1574 }
1575 if (d_a >= 0.4 && d_a < 0.41) {
1576 sigm_v1_ab[0][pair_v1] = 0.08 * 0.08;
1577 }
1578 if (d_a >= 0.41) {
1579 sigm_v1_ab[0][pair_v1] = 0.10 * 0.10;
1580 }
1581
1582 if (d_b < 0.02) {
1583 sigm_v1_ab[1][pair_v1] = 0.08 * 0.08;
1584 }
1585 if (d_b >= 0.02 && d_b < 0.1) {
1586 sigm_v1_ab[1][pair_v1] = 0.06 * 0.06;
1587 }
1588 if (d_b >= 0.1 && d_b < 0.4) {
1589 sigm_v1_ab[1][pair_v1] = 0.025 * 0.025;
1590 }
1591 if (d_b >= 0.4 && d_b < 0.41) {
1592 sigm_v1_ab[1][pair_v1] = 0.08 * 0.08;
1593 }
1594 if (d_b >= 0.41) {
1595 sigm_v1_ab[1][pair_v1] = 0.10 * 0.10;
1596 }
1597 pair_v1++;
1598 used_va1[i] = true;
1599 used_vb1[j] = true;
1600
1601 } // j for2
1602 } // i for1 V.
1603
1604 // reconstruct single V-plane hits
1605
1606 for (Int_t i = 0; i < it_va1; ++i) {
1607 if (used_va1[i])
1608 continue;
1609
1610 double d_a = 0;
1611 double d_b = 0;
1612 int func_nr_a = -1;
1613 for (int t_it = 0; t_it < 4; t_it++) {
1614 if (time_va1[i] >= t_dc[t_it][6] && time_va1[i] < t_dc[t_it + 1][6]) {
1615 func_nr_a = t_it;
1616 break;
1617 }
1618 }
1619 if (func_nr_a == 0)
1620 d_a = 0;
1621 if (func_nr_a == 3)
1622 d_a = scale[6]
1623 * (pol_par_dc[2][0][6] + pol_par_dc[2][1][6] * time_va1[i]
1624 + pol_par_dc[2][2][6] * time_va1[i] * time_va1[i]
1625 + pol_par_dc[2][3][6] * time_va1[i] * time_va1[i] * time_va1[i]
1626 + pol_par_dc[2][4][6] * time_va1[i] * time_va1[i] * time_va1[i] * time_va1[i]);
1627 if (func_nr_a == 1 || func_nr_a == 2)
1628 d_a = scale[6]
1629 * (pol_par_dc[1][0][6] + pol_par_dc[1][1][6] * time_va1[i]
1630 + pol_par_dc[1][2][6] * time_va1[i] * time_va1[i]
1631 + pol_par_dc[1][3][6] * time_va1[i] * time_va1[i] * time_va1[i]
1632 + pol_par_dc[1][4][6] * time_va1[i] * time_va1[i] * time_va1[i] * time_va1[i]);
1633
1634 v1_single[0][single_va1] = wirenr_va1[i] - 119 + d_a;
1635 v1_single[0][single_va1 + 1] = wirenr_va1[i] - 119 - d_a;
1636
1637 if (d_a < 0.02) {
1638 sigm_v1_single[0][single_va1] = 0.08 * 0.08;
1639 sigm_v1_single[0][single_va1 + 1] = 0.08 * 0.08;
1640 }
1641 if (d_a >= 0.02 && d_a < 0.1) {
1642 sigm_v1_single[0][single_va1] = 0.06 * 0.06;
1643 sigm_v1_single[0][single_va1 + 1] = 0.06 * 0.06;
1644 }
1645 if (d_a >= 0.1 && d_a < 0.4) {
1646 sigm_v1_single[0][single_va1] = 0.025 * 0.025;
1647 sigm_v1_single[0][single_va1 + 1] = 0.025 * 0.025;
1648 }
1649 if (d_a >= 0.4 && d_a < 0.41) {
1650 sigm_v1_single[0][single_va1] = 0.08 * 0.08;
1651 sigm_v1_single[0][single_va1 + 1] = 0.08 * 0.08;
1652 }
1653 if (d_a >= 0.41) {
1654 sigm_v1_single[0][single_va1] = 0.10 * 0.10;
1655 sigm_v1_single[0][single_va1 + 1] = 0.10 * 0.10;
1656 }
1657 single_va1 += 2;
1658 } // for single va
1659
1660 for (Int_t j = 0; j < it_vb1; ++j) {
1661 if (used_vb1[j])
1662 continue;
1663 double d_a = 0;
1664 double d_b = 0;
1665 int func_nr_b = -1;
1666 for (int t_it = 0; t_it < 4; t_it++) {
1667 if (time_vb1[j] >= t_dc[t_it][7] && time_vb1[j] < t_dc[t_it + 1][7]) {
1668 func_nr_b = t_it;
1669 break;
1670 }
1671 }
1672 if (func_nr_b == 1 || func_nr_b == 2)
1673 d_b = scale[7]
1674 * (pol_par_dc[1][0][7] + pol_par_dc[1][1][7] * time_vb1[j]
1675 + pol_par_dc[1][2][7] * time_vb1[j] * time_vb1[j]
1676 + pol_par_dc[1][3][7] * time_vb1[j] * time_vb1[j] * time_vb1[j]
1677 + pol_par_dc[1][4][7] * time_vb1[j] * time_vb1[j] * time_vb1[j] * time_vb1[j]);
1678
1679 if (func_nr_b == 0)
1680 d_b = 0;
1681 if (func_nr_b == 3)
1682 d_b = scale[7]
1683 * (pol_par_dc[2][0][7] + pol_par_dc[2][1][7] * time_vb1[j]
1684 + pol_par_dc[2][2][7] * time_vb1[j] * time_vb1[j]
1685 + pol_par_dc[2][3][7] * time_vb1[j] * time_vb1[j] * time_vb1[j]
1686 + pol_par_dc[2][4][7] * time_vb1[j] * time_vb1[j] * time_vb1[j] * time_vb1[j]);
1687
1688 v1_single[1][single_vb1] = wirenr_vb1[j] - 118.5 + d_b;
1689 v1_single[1][single_vb1 + 1] = wirenr_vb1[j] - 118.5 - d_b;
1690 if (d_b < 0.02) {
1691 sigm_v1_single[1][single_vb1] = 0.08 * 0.08;
1692 sigm_v1_single[1][single_vb1 + 1] = 0.08 * 0.08;
1693 }
1694 if (d_b >= 0.02 && d_b < 0.1) {
1695 sigm_v1_single[1][single_vb1] = 0.06 * 0.06;
1696 sigm_v1_single[1][single_vb1 + 1] = 0.06 * 0.06;
1697 }
1698 if (d_b >= 0.1 && d_b < 0.4) {
1699 sigm_v1_single[1][single_vb1] = 0.0250 * 0.0250;
1700 sigm_v1_single[1][single_vb1 + 1] = 0.0250 * 0.0250;
1701 }
1702 if (d_b >= 0.4 && d_b < 0.41) {
1703 sigm_v1_single[1][single_vb1] = 0.08 * 0.08;
1704 sigm_v1_single[1][single_vb1 + 1] = 0.08 * 0.08;
1705 }
1706 if (d_b >= 0.41) {
1707 sigm_v1_single[1][single_vb1] = 0.10 * 0.10;
1708 sigm_v1_single[1][single_vb1 + 1] = 0.10 * 0.10;
1709 }
1710 single_vb1 += 2;
1711
1712 } // i for1 V.
1713
1714 // --- DC1 Segment Building ------
1715
1716 bool found_8p_seg = false;
1717 float x_slope = 0.0;
1718 float y_slope = 0.0;
1719
1720 if (pair_x1 * pair_y1 > 0) { //(x,y)first
1721 for (int i = 0; i < pair_x1; i++) {
1722 if (nDC1_segments > 48)
1723 break;
1724 float x_coord = (x1_ab[0][i] + x1_ab[1][i]) / 2;
1725 float XU = x_coord + x_slope * 0.5 * (z_loc[5] + z_loc[4] - z_loc[1] - z_loc[0]);
1726 float XV = x_coord + x_slope * 0.5 * (z_loc[7] + z_loc[6] - z_loc[1] - z_loc[0]);
1727 for (int j = 0; j < pair_y1; j++) {
1728 if (nDC1_segments > 48)
1729 break;
1730 float y_coord = (y1_ab[0][j] + y1_ab[1][j]) / 2;
1731 float YU = y_coord + y_slope * 0.5 * (z_loc[5] + z_loc[4] - z_loc[3] - z_loc[2]);
1732 float YV = y_coord + y_slope * 0.5 * (z_loc[7] + z_loc[6] - z_loc[3] - z_loc[2]);
1733 float u_est = isqrt_2 * (YU - XU);
1734 float v_est = isqrt_2 * (YV + XV);
1735 bool foundU = false;
1736
1737 if (pair_u1 > 0) {
1738 double dU_thresh = 1.3;
1739 for (int k = 0; k < pair_u1; k++) {
1740 float u_coord = (u1_ab[0][k] + u1_ab[1][k]) / 2;
1741 if (fabs(u_coord - u_est) > dU_thresh)
1742 continue;
1743 dU_thresh = fabs(u_coord - u_est);
1744
1745 rh_segDC1[0][nDC1_segments] = x1_ab[0][i];
1746 rh_segDC1[1][nDC1_segments] = x1_ab[1][i];
1747 rh_segDC1[2][nDC1_segments] = y1_ab[0][j];
1748 rh_segDC1[3][nDC1_segments] = y1_ab[1][j];
1749 rh_segDC1[4][nDC1_segments] = u1_ab[0][k];
1750 rh_segDC1[5][nDC1_segments] = u1_ab[1][k];
1751 rh_sigm_segDC1[0][nDC1_segments] = sigm_x1_ab[0][i];
1752 rh_sigm_segDC1[1][nDC1_segments] = sigm_x1_ab[1][i];
1753 rh_sigm_segDC1[2][nDC1_segments] = sigm_y1_ab[0][j];
1754 rh_sigm_segDC1[3][nDC1_segments] = sigm_y1_ab[1][j];
1755 rh_sigm_segDC1[4][nDC1_segments] = sigm_u1_ab[0][k];
1756 rh_sigm_segDC1[5][nDC1_segments] = sigm_u1_ab[1][k];
1757 foundU = true;
1758 if (nDC1_segments > 48)
1759 break;
1760 } // k
1761 } //(pair_u1>0)
1762 if (found_8p_seg && !foundU)
1763 continue;
1764 bool foundV = false;
1765 if (pair_v1 > 0) {
1766 double dV_thresh = 1.3;
1767 for (int m = 0; m < pair_v1; m++) {
1768 if (nDC1_segments > 48)
1769 break;
1770 float v_coord = (v1_ab[0][m] + v1_ab[1][m]) / 2;
1771 if (fabs(v_coord - v_est) > dV_thresh)
1772 continue;
1773 dV_thresh = fabs(v_coord - v_est);
1774
1775 foundV = true;
1776 rh_segDC1[0][nDC1_segments] = x1_ab[0][i];
1777 rh_segDC1[1][nDC1_segments] = x1_ab[1][i];
1778 rh_segDC1[2][nDC1_segments] = y1_ab[0][j];
1779 rh_segDC1[3][nDC1_segments] = y1_ab[1][j];
1780 rh_segDC1[6][nDC1_segments] = v1_ab[0][m];
1781 rh_segDC1[7][nDC1_segments] = v1_ab[1][m];
1782 rh_sigm_segDC1[0][nDC1_segments] = sigm_x1_ab[0][i];
1783 rh_sigm_segDC1[1][nDC1_segments] = sigm_x1_ab[1][i];
1784 rh_sigm_segDC1[2][nDC1_segments] = sigm_y1_ab[0][j];
1785 rh_sigm_segDC1[3][nDC1_segments] = sigm_y1_ab[1][j];
1786 rh_sigm_segDC1[6][nDC1_segments] = sigm_v1_ab[0][m];
1787 rh_sigm_segDC1[7][nDC1_segments] = sigm_v1_ab[1][m];
1788 if (!foundU) {
1789
1790 float min_a = 999;
1791 float min_b = 999;
1792
1793 for (int kk = 0; kk < single_ua1; kk++) {
1794 if (fabs(u1_single[0][kk] - u_est) > 1.5)
1795 continue; //????? 0.5 needs to be reviewed
1796 if (fabs(u1_single[0][kk] - u_est) < min_a) {
1797 min_a = fabs(u1_single[0][kk] - u_est);
1798 rh_segDC1[4][nDC1_segments] = u1_single[0][kk];
1799 rh_sigm_segDC1[4][nDC1_segments] = sigm_u1_single[0][kk];
1800
1801 foundU = true;
1802 }
1803 } // for kk
1804 for (int kk = 0; kk < single_ub1; kk++) {
1805 if (fabs(u1_single[1][kk] - u_est) > 1.5)
1806 continue; //????? 0.5 needs to be reviewed
1807 if (fabs(u1_single[1][kk] - u_est) < min_b) {
1808 min_a = fabs(u1_single[1][kk] - u_est);
1809 rh_segDC1[5][nDC1_segments] = u1_single[1][kk];
1810 rh_sigm_segDC1[5][nDC1_segments] = sigm_u1_single[1][kk];
1811
1812 foundU = true;
1813 }
1814 } // for kk
1815
1816 if (nDC1_segments > 48)
1817 break;
1818 } //! foundU
1819
1820 if (nDC1_segments > 48)
1821 break;
1822
1823 } // m
1824
1825 } //(pair_v1>0)
1826
1827 if (!foundV && foundU) {
1828 float min_a = 999;
1829 float min_b = 999;
1830 for (int kk = 0; kk < single_va1; kk++) {
1831 if (fabs(v1_single[0][kk] - v_est) > 1.5)
1832 continue; //????? 0.5 needs to be reviewed
1833
1834 if (fabs(v1_single[0][kk] - v_est) < min_a) {
1835 min_a = fabs(v1_single[0][kk] - v_est);
1836 rh_segDC1[6][nDC1_segments] = v1_single[0][kk];
1837 rh_sigm_segDC1[6][nDC1_segments] = sigm_v1_single[0][kk];
1838 foundV = true;
1839 }
1840 } // for kk
1841 for (int kk = 0; kk < single_vb1; kk++) {
1842 if (fabs(v1_single[1][kk] - v_est) > 1.5)
1843 continue; //????? 0.5 needs to be reviewed
1844 if (fabs(v1_single[1][kk] - v_est) < min_b) {
1845 min_b = fabs(v1_single[1][kk] - v_est);
1846 rh_segDC1[7][nDC1_segments] = v1_single[1][kk];
1847 rh_sigm_segDC1[7][nDC1_segments] = sigm_v1_single[1][kk];
1848 foundV = true;
1849 }
1850 } // for kk
1851
1852 } //! foundV
1853 if (foundV || foundU)
1854 nDC1_segments++;
1855 } // j
1856 } // i
1857
1858 } //(x,y)first
1859
1860 if (pair_u1 * pair_v1 > 0 && !found_8p_seg) { // (u,v) first
1861
1862 for (int i = 0; i < pair_u1; i++) {
1863 if (nDC1_segments > 48)
1864 break;
1865 float u_coord = (u1_ab[0][i] + u1_ab[1][i]) / 2;
1866 for (int j = 0; j < pair_v1; j++) {
1867 if (nDC1_segments > 48)
1868 break;
1869 float v_coord = (v1_ab[0][j] + v1_ab[1][j]) / 2;
1870 bool foundX = false;
1871 float y_est = isqrt_2 * (u_coord + v_coord);
1872 float x_est = isqrt_2 * (v_coord - u_coord);
1873 if (pair_x1 > 0) {
1874
1875 double dX_thresh = 1.5;
1876 for (int k = 0; k < pair_x1; k++) {
1877 float x_coord = (x1_ab[0][k] + x1_ab[1][k]) / 2;
1878 if (nDC1_segments > 48)
1879 break;
1880 if (fabs(x_coord - x_est) > dX_thresh)
1881 continue;
1882 dX_thresh = fabs(x_coord - x_est);
1883
1884 rh_segDC1[0][nDC1_segments] = x1_ab[0][k];
1885 rh_segDC1[1][nDC1_segments] = x1_ab[1][k];
1886 rh_segDC1[4][nDC1_segments] = u1_ab[0][i];
1887 rh_segDC1[5][nDC1_segments] = u1_ab[1][i];
1888 rh_segDC1[6][nDC1_segments] = v1_ab[0][j];
1889 rh_segDC1[7][nDC1_segments] = v1_ab[1][j];
1890 rh_sigm_segDC1[0][nDC1_segments] = sigm_x1_ab[0][k];
1891 rh_sigm_segDC1[1][nDC1_segments] = sigm_x1_ab[1][k];
1892 rh_sigm_segDC1[4][nDC1_segments] = sigm_u1_ab[0][i];
1893 rh_sigm_segDC1[5][nDC1_segments] = sigm_u1_ab[1][i];
1894 rh_sigm_segDC1[6][nDC1_segments] = sigm_v1_ab[0][j];
1895 rh_sigm_segDC1[7][nDC1_segments] = sigm_v1_ab[1][j];
1896
1897 foundX = true;
1898 if (nDC1_segments > 48)
1899 break;
1900 } // k
1901
1902 } //(pair_x1>0)
1903 if (found_8p_seg && !foundX)
1904 continue;
1905 bool foundY = false;
1906 if (pair_y1 > 0) {
1907
1908 double dY_thresh = 1.0;
1909 for (int m = 0; m < pair_y1; m++) {
1910 if (nDC1_segments > 48)
1911 break;
1912 float y_coord = (y1_ab[0][m] + y1_ab[1][m]) / 2;
1913 if (fabs(y_coord - y_est) > dY_thresh)
1914 continue;
1915 dY_thresh = fabs(y_coord - y_est);
1916
1917 foundY = true;
1918 if (nDC1_segments > 48)
1919 break;
1920 rh_segDC1[2][nDC1_segments] = y1_ab[0][m];
1921 rh_segDC1[3][nDC1_segments] = y1_ab[1][m];
1922 rh_segDC1[4][nDC1_segments] = u1_ab[0][i];
1923 rh_segDC1[5][nDC1_segments] = u1_ab[1][i];
1924 rh_segDC1[6][nDC1_segments] = v1_ab[0][j];
1925 rh_segDC1[7][nDC1_segments] = v1_ab[1][j];
1926 rh_sigm_segDC1[2][nDC1_segments] = sigm_y1_ab[0][m];
1927 rh_sigm_segDC1[3][nDC1_segments] = sigm_y1_ab[1][m];
1928 rh_sigm_segDC1[4][nDC1_segments] = sigm_u1_ab[0][i];
1929 rh_sigm_segDC1[5][nDC1_segments] = sigm_u1_ab[1][i];
1930 rh_sigm_segDC1[6][nDC1_segments] = sigm_v1_ab[0][j];
1931 rh_sigm_segDC1[7][nDC1_segments] = sigm_v1_ab[1][j];
1932
1933 if (!foundX) {
1934
1935 float min_a = 999;
1936 float min_b = 999;
1937 for (int kk = 0; kk < single_xa1; kk++) {
1938 if (fabs(x1_single[0][kk] - x_est) > 1.5)
1939 continue; //????? 0.5 needs to be reviewed
1940 if (fabs(x1_single[0][kk] - x_est) < min_a) {
1941 min_a = fabs(x1_single[0][kk] - x_est);
1942 rh_segDC1[0][nDC1_segments] = x1_single[0][kk];
1943 rh_sigm_segDC1[0][nDC1_segments] = sigm_x1_single[0][kk];
1944 foundX = true;
1945 }
1946 } // for kk
1947 for (int kk = 0; kk < single_xb1; kk++) {
1948 if (fabs(x1_single[1][kk] - x_est) > 1.5)
1949 continue; //????? 0.5 needs to be reviewed
1950 if (fabs(x1_single[1][kk] - x_est) < min_b) {
1951 min_a = fabs(x1_single[1][kk] - x_est);
1952 rh_segDC1[1][nDC1_segments] = x1_single[1][kk];
1953 rh_sigm_segDC1[1][nDC1_segments] = sigm_x1_single[1][kk];
1954 foundX = true;
1955 }
1956 } // for kk
1957 if (nDC1_segments > 48)
1958 break;
1959 } //! foundX
1960
1961 } // m
1962 if (!foundY && foundX) {
1963 float min_a = 999;
1964 float min_b = 999;
1965 for (int kk = 0; kk < single_ya1; kk++) {
1966 if (fabs(y1_single[0][kk] - y_est) > 1.5)
1967 continue; //????? 0.5 needs to be reviewed
1968 if (fabs(y1_single[0][kk] - y_est) < min_a) {
1969 min_a = fabs(y1_single[0][kk] - y_est);
1970 rh_segDC1[2][nDC1_segments] = y1_single[0][kk];
1971 rh_sigm_segDC1[2][nDC1_segments] = sigm_y1_single[0][kk];
1972 foundY = true;
1973 }
1974 } // for kk
1975 for (int kk = 0; kk < single_yb1; kk++) {
1976 if (fabs(y1_single[1][kk] - y_est) > 1.5)
1977 continue; //????? 0.5 needs to be reviewed
1978 if (fabs(y1_single[1][kk] - y_est) < min_b) {
1979 min_b = fabs(y1_single[1][kk] - y_est);
1980 rh_segDC1[3][nDC1_segments] = y1_single[1][kk];
1981 rh_sigm_segDC1[3][nDC1_segments] = sigm_y1_single[1][kk];
1982 foundY = true;
1983 }
1984 } // for kk
1985 } //! foundY
1986 } //(pair_y1>0)
1987 if (foundX || foundY)
1988 nDC1_segments++;
1989 } // j
1990 } // i
1991 } //(u,v) first
1992
1993 ///-------------DCH2 hit reco------------------///
1994 float par_ab1[4][150];
1995 float par_ab2[4][150];
1996 for (int i = 0; i < 8; i++) {
1997 for (int j = 0; j < 150; j++) {
1998 if (i > 1)
1999 continue;
2000 par_ab1[i][j] = -999;
2001 par_ab2[i][j] = -999;
2002 }
2003 }
2004 int pair_x2 = 0;
2005 int pair_y2 = 0;
2006 int pair_u2 = 0;
2007 int pair_v2 = 0;
2008 int single_xa2 = 0;
2009 int single_ya2 = 0;
2010 int single_ua2 = 0;
2011 int single_va2 = 0;
2012 int single_xb2 = 0;
2013 int single_yb2 = 0;
2014 int single_ub2 = 0;
2015 int single_vb2 = 0;
2016
2017 float xa2_pm[2] = {-999};
2018 float xb2_pm[2] = {-999};
2019 float ya2_pm[2] = {-999};
2020 float yb2_pm[2] = {-999};
2021 float ua2_pm[2] = {-999};
2022 float ub2_pm[2] = {-999};
2023 float va2_pm[2] = {-999};
2024 float vb2_pm[2] = {-999};
2025
2026 // --- X ---
2027
2028 for (Int_t i = 0; i < it_xa2; ++i) {
2029 for (Int_t j = 0; j < it_xb2; ++j) {
2030 if (pair_x2 > 48)
2031 break;
2032 if ((wirenr_xa2[i] != wirenr_xb2[j] && wirenr_xa2[i] != wirenr_xb2[j] + 1))
2033 continue;
2034 int func_nr_a = -1;
2035 int func_nr_b = -1;
2036 for (int t_it = 0; t_it < 4; t_it++) {
2037 if (time_xa2[i] >= t_dc[t_it][8] && time_xa2[i] < t_dc[t_it + 1][8]) {
2038 func_nr_a = t_it;
2039 break;
2040 }
2041 }
2042 double time_xa = time_xa2[i];
2043 double d_a = 0;
2044 double d_b = 0;
2045
2046 if (func_nr_a == 1 || func_nr_a == 2)
2047 d_a =
2048 scale[8]
2049 * (pol_par_dc[1][0][8] + pol_par_dc[1][1][8] * time_xa + pol_par_dc[1][2][8] * time_xa * time_xa
2050 + pol_par_dc[1][3][8] * time_xa * time_xa * time_xa
2051 + pol_par_dc[1][4][8] * time_xa * time_xa * time_xa * time_xa);
2052 if (func_nr_a == 0)
2053 d_a = 0;
2054 if (func_nr_a == 3)
2055 d_a =
2056 scale[8]
2057 * (pol_par_dc[2][0][8] + pol_par_dc[2][1][8] * time_xa + pol_par_dc[2][2][8] * time_xa * time_xa
2058 + pol_par_dc[2][3][8] * time_xa * time_xa * time_xa
2059 + pol_par_dc[2][4][8] * time_xa * time_xa * time_xa * time_xa);
2060
2061 for (int t_it = 0; t_it < 4; t_it++) {
2062 if (time_xb2[j] >= t_dc[t_it][9] && time_xb2[j] < t_dc[t_it + 1][9]) {
2063 func_nr_b = t_it;
2064 break;
2065 }
2066 }
2067 time_xa = time_xb2[j];
2068
2069 if (func_nr_b == 1 || func_nr_b == 2)
2070 d_b =
2071 scale[9]
2072 * (pol_par_dc[1][0][9] + pol_par_dc[1][1][9] * time_xa + pol_par_dc[1][2][9] * time_xa * time_xa
2073 + pol_par_dc[1][3][9] * time_xa * time_xa * time_xa
2074 + pol_par_dc[1][4][9] * time_xa * time_xa * time_xa * time_xa);
2075 if (func_nr_b == 0)
2076 d_b = 0;
2077 if (func_nr_b == 3)
2078 d_b =
2079 scale[9]
2080 * (pol_par_dc[2][0][9] + pol_par_dc[2][1][9] * time_xa + pol_par_dc[2][2][9] * time_xa * time_xa
2081 + pol_par_dc[2][3][9] * time_xa * time_xa * time_xa
2082 + pol_par_dc[2][4][9] * time_xa * time_xa * time_xa * time_xa);
2083
2084 xa2_pm[0] = wirenr_xa2[i] - 119 + d_a;
2085 xa2_pm[1] = wirenr_xa2[i] - 119 - d_a;
2086 xb2_pm[0] = wirenr_xb2[j] - 118.5 + d_b;
2087 xb2_pm[1] = wirenr_xb2[j] - 118.5 - d_b;
2088
2089 double dmin1 = 999;
2090 for (int k = 0; k < 2; k++) {
2091 for (int m = 0; m < 2; m++) {
2092 if (fabs(xa2_pm[k] - xb2_pm[m]) < dmin1) {
2093 dmin1 = fabs(xa2_pm[k] - xb2_pm[m]);
2094 x2_ab[0][pair_x2] = xa2_pm[k];
2095 x2_ab[1][pair_x2] = xb2_pm[m];
2096 }
2097 }
2098 }
2099 if (d_a < 0.02) {
2100 sigm_x2_ab[0][pair_x2] = 0.08 * 0.08;
2101 }
2102 if (d_a >= 0.02 && d_a < 0.1) {
2103 sigm_x2_ab[0][pair_x2] = 0.06 * 0.06;
2104 }
2105 if (d_a >= 0.1 && d_a < 0.4) {
2106 sigm_x2_ab[0][pair_x2] = 0.025 * 0.025;
2107 }
2108 if (d_a >= 0.4 && d_a < 0.41) {
2109 sigm_x2_ab[0][pair_x2] = 0.08 * 0.08;
2110 }
2111 if (d_a >= 0.41) {
2112 sigm_x2_ab[0][pair_x2] = 0.10 * 0.10;
2113 }
2114
2115 if (d_b < 0.02) {
2116 sigm_x2_ab[1][pair_x2] = 0.08 * 0.08;
2117 }
2118 if (d_b >= 0.02 && d_b < 0.1) {
2119 sigm_x2_ab[1][pair_x2] = 0.06 * 0.06;
2120 }
2121 if (d_b >= 0.1 && d_b < 0.4) {
2122 sigm_x2_ab[1][pair_x2] = 0.025 * 0.025;
2123 }
2124 if (d_b >= 0.4 && d_b < 0.41) {
2125 sigm_x2_ab[1][pair_x2] = 0.08 * 0.08;
2126 }
2127 if (d_b >= 0.41) {
2128 sigm_x2_ab[1][pair_x2] = 0.10 * 0.10;
2129 }
2130 pair_x2++;
2131 used_xa2[i] = true;
2132 used_xb2[j] = true;
2133 } // j for2
2134 } // i for1 X.
2135
2136 // reconstruct single X-plane hits
2137
2138 for (Int_t i = 0; i < it_xa2; ++i) {
2139 if (used_xa2[i])
2140 continue;
2141 double d_a = 0;
2142 double d_b = 0;
2143
2144 int func_nr_a = -1;
2145 for (int t_it = 0; t_it < 4; t_it++) {
2146 if (time_xa2[i] >= t_dc[t_it][8] && time_xa2[i] < t_dc[t_it + 1][8]) {
2147 func_nr_a = t_it;
2148 break;
2149 }
2150 }
2151 double time_xa = time_xa2[i];
2152 if (func_nr_a == 1 || func_nr_a == 2)
2153 d_a = scale[8]
2154 * (pol_par_dc[1][0][8] + pol_par_dc[1][1][8] * time_xa + pol_par_dc[1][2][8] * time_xa * time_xa
2155 + pol_par_dc[1][3][8] * time_xa * time_xa * time_xa
2156 + pol_par_dc[1][4][8] * time_xa * time_xa * time_xa * time_xa);
2157 if (func_nr_a == 0)
2158 d_a = 0;
2159 if (func_nr_a == 3)
2160 d_a = scale[8]
2161 * (pol_par_dc[2][0][8] + pol_par_dc[2][1][8] * time_xa + pol_par_dc[2][2][8] * time_xa * time_xa
2162 + pol_par_dc[2][3][8] * time_xa * time_xa * time_xa
2163 + pol_par_dc[2][4][8] * time_xa * time_xa * time_xa * time_xa);
2164
2165 x2_single[0][single_xa2] = wirenr_xa2[i] - 118.5 + d_a;
2166 x2_single[0][single_xa2 + 1] = wirenr_xa2[i] - 118.5 - d_a;
2167 if (d_a < 0.02) {
2168 sigm_x2_single[0][single_xa2] = 0.08 * 0.08;
2169 sigm_x2_single[0][single_xa2 + 1] = 0.08 * 0.08;
2170 }
2171 if (d_a >= 0.02 && d_a < 0.1) {
2172 sigm_x2_single[0][single_xa2] = 0.06 * 0.06;
2173 sigm_x2_single[0][single_xa2 + 1] = 0.06 * 0.06;
2174 }
2175 if (d_a >= 0.1 && d_a < 0.4) {
2176 sigm_x2_single[0][single_xa2] = 0.0250 * 0.0250;
2177 sigm_x2_single[0][single_xa2 + 1] = 0.0250 * 0.0250;
2178 }
2179 if (d_a >= 0.4 && d_a < 0.41) {
2180 sigm_x2_single[0][single_xa2] = 0.08 * 0.08;
2181 sigm_x2_single[0][single_xa2 + 1] = 0.08 * 0.08;
2182 }
2183 if (d_a >= 0.41) {
2184 sigm_x2_single[0][single_xa2] = 0.10 * 0.10;
2185 sigm_x2_single[0][single_xa2 + 1] = 0.10 * 0.10;
2186 }
2187 single_xa2 += 2;
2188 } // for single xa
2189
2190 for (Int_t j = 0; j < it_xb2; ++j) {
2191 if (used_xb2[j])
2192 continue;
2193 int func_nr_b = -1;
2194 double d_a = 0;
2195 double d_b = 0;
2196
2197 for (int t_it = 0; t_it < 4; t_it++) {
2198 if (time_xb2[j] >= t_dc[t_it][9] && time_xb2[j] < t_dc[t_it + 1][9]) {
2199 func_nr_b = t_it;
2200 break;
2201 }
2202 }
2203 double time_xa = time_xb2[j];
2204 if (func_nr_b == 1 || func_nr_b == 2)
2205 d_b = scale[9]
2206 * (pol_par_dc[1][0][9] + pol_par_dc[1][1][9] * time_xa + pol_par_dc[1][2][9] * time_xa * time_xa
2207 + pol_par_dc[1][3][9] * time_xa * time_xa * time_xa
2208 + pol_par_dc[1][4][9] * time_xa * time_xa * time_xa * time_xa);
2209 if (func_nr_b == 0)
2210 d_b = 0;
2211 if (func_nr_b == 3)
2212 d_b = scale[9]
2213 * (pol_par_dc[2][0][9] + pol_par_dc[2][1][9] * time_xa + pol_par_dc[2][2][9] * time_xa * time_xa
2214 + pol_par_dc[2][3][9] * time_xa * time_xa * time_xa
2215 + pol_par_dc[2][4][9] * time_xa * time_xa * time_xa * time_xa);
2216
2217 x2_single[1][single_xb2] = wirenr_xb2[j] - 118.5 + d_b;
2218 x2_single[1][single_xb2 + 1] = wirenr_xb2[j] - 118.5 - d_b;
2219 if (d_b < 0.02) {
2220 sigm_x2_single[1][single_xb2] = 0.08 * 0.08;
2221 sigm_x2_single[1][single_xb2 + 1] = 0.08 * 0.08;
2222 }
2223 if (d_b >= 0.02 && d_b < 0.1) {
2224 sigm_x2_single[1][single_xb2] = 0.06 * 0.06;
2225 sigm_x2_single[1][single_xb2 + 1] = 0.06 * 0.06;
2226 }
2227 if (d_b >= 0.1 && d_b < 0.4) {
2228 sigm_x2_single[1][single_xb2] = 0.0250 * 0.0250;
2229 sigm_x2_single[1][single_xb2 + 1] = 0.0250 * 0.0250;
2230 }
2231 if (d_b >= 0.4 && d_b < 0.41) {
2232 sigm_x2_single[1][single_xb2] = 0.08 * 0.08;
2233 sigm_x2_single[1][single_xb2 + 1] = 0.08 * 0.08;
2234 }
2235 if (d_b >= 0.41) {
2236 sigm_x2_single[1][single_xb2] = 0.1 * 0.1;
2237 sigm_x2_single[1][single_xb2 + 1] = 0.1 * 0.1;
2238 }
2239 single_xb2 += 2;
2240 } // j xb2
2241
2242 //------Y-------
2243
2244 for (Int_t i = 0; i < it_ya2; ++i) {
2245 for (Int_t j = 0; j < it_yb2; ++j) {
2246 if (pair_y2 > 48)
2247 break;
2248 if ((wirenr_ya2[i] != wirenr_yb2[j] && wirenr_ya2[i] != wirenr_yb2[j] + 1))
2249 continue;
2250 double d_a = 0;
2251 double d_b = 0;
2252
2253 int func_nr_a = -1;
2254 int func_nr_b = -1;
2255 for (int t_it = 0; t_it < 4; t_it++) {
2256 if (time_ya2[i] >= t_dc[t_it][10] && time_ya2[i] < t_dc[t_it + 1][10]) {
2257 func_nr_a = t_it;
2258 break;
2259 }
2260 }
2261 if (func_nr_a == 1 || func_nr_a == 2)
2262 d_a = scale[10]
2263 * (pol_par_dc[1][0][10] + pol_par_dc[1][1][10] * time_ya2[i]
2264 + pol_par_dc[1][2][10] * time_ya2[i] * time_ya2[i]
2265 + pol_par_dc[1][3][10] * time_ya2[i] * time_ya2[i] * time_ya2[i]
2266 + pol_par_dc[1][4][10] * time_ya2[i] * time_ya2[i] * time_ya2[i] * time_ya2[i]);
2267
2268 if (func_nr_a == 0)
2269 d_a = 0;
2270 if (func_nr_a == 3)
2271 d_a = scale[10]
2272 * (pol_par_dc[2][0][10] + pol_par_dc[2][1][10] * time_ya2[i]
2273 + pol_par_dc[2][2][10] * time_ya2[i] * time_ya2[i]
2274 + pol_par_dc[2][3][10] * time_ya2[i] * time_ya2[i] * time_ya2[i]
2275 + pol_par_dc[2][4][10] * time_ya2[i] * time_ya2[i] * time_ya2[i] * time_ya2[i]);
2276
2277 for (int t_it = 0; t_it < 4; t_it++) {
2278 if (time_yb2[j] >= t_dc[t_it][11] && time_yb2[j] < t_dc[t_it + 1][11]) {
2279 func_nr_b = t_it;
2280 break;
2281 }
2282 }
2283 if (func_nr_b == 1 || func_nr_b == 2)
2284 d_b = scale[11]
2285 * (pol_par_dc[1][0][11] + pol_par_dc[1][1][11] * time_yb2[j]
2286 + pol_par_dc[1][2][11] * time_yb2[j] * time_yb2[j]
2287 + pol_par_dc[1][3][11] * time_yb2[j] * time_yb2[j] * time_yb2[j]
2288 + pol_par_dc[1][4][11] * time_yb2[j] * time_yb2[j] * time_yb2[j] * time_yb2[j]);
2289 if (func_nr_b == 0)
2290 d_b = 0;
2291 if (func_nr_b == 3)
2292 d_b = scale[11]
2293 * (pol_par_dc[2][0][11] + pol_par_dc[2][1][11] * time_yb2[j]
2294 + pol_par_dc[2][2][11] * time_yb2[j] * time_yb2[j]
2295 + pol_par_dc[2][3][11] * time_yb2[j] * time_yb2[j] * time_yb2[j]
2296 + pol_par_dc[2][4][11] * time_yb2[j] * time_yb2[j] * time_yb2[j] * time_yb2[j]);
2297
2298 ya2_pm[0] = wirenr_ya2[i] - 119 + d_a;
2299 ya2_pm[1] = wirenr_ya2[i] - 119 - d_a;
2300 yb2_pm[0] = wirenr_yb2[j] - 118.5 + d_b;
2301 yb2_pm[1] = wirenr_yb2[j] - 118.5 - d_b;
2302
2303 double dmin1 = 999;
2304 for (int k = 0; k < 2; k++) {
2305 for (int m = 0; m < 2; m++) {
2306 if (fabs(ya2_pm[k] - yb2_pm[m]) < dmin1) {
2307 dmin1 = fabs(ya2_pm[k] - yb2_pm[m]);
2308 y2_ab[0][pair_y2] = ya2_pm[k];
2309 y2_ab[1][pair_y2] = yb2_pm[m];
2310 }
2311 }
2312 }
2313 if (d_a < 0.02) {
2314 sigm_y2_ab[0][pair_y2] = 0.08 * 0.08;
2315 }
2316 if (d_a >= 0.02 && d_a < 0.1) {
2317 sigm_y2_ab[0][pair_y2] = 0.06 * 0.06;
2318 }
2319 if (d_a >= 0.1 && d_a < 0.4) {
2320 sigm_y2_ab[0][pair_y2] = 0.025 * 0.025;
2321 }
2322 if (d_a >= 0.4 && d_a < 0.41) {
2323 sigm_y2_ab[0][pair_y2] = 0.08 * 0.08;
2324 }
2325 if (d_a >= 0.41) {
2326 sigm_y2_ab[0][pair_y2] = 0.10 * 0.10;
2327 }
2328
2329 if (d_b < 0.02) {
2330 sigm_y2_ab[1][pair_y2] = 0.08 * 0.08;
2331 }
2332 if (d_b >= 0.02 && d_b < 0.1) {
2333 sigm_y2_ab[1][pair_y2] = 0.06 * 0.06;
2334 }
2335 if (d_b >= 0.1 && d_b < 0.4) {
2336 sigm_y2_ab[1][pair_y2] = 0.025 * 0.025;
2337 }
2338 if (d_b >= 0.4 && d_b < 0.41) {
2339 sigm_y2_ab[1][pair_y2] = 0.08 * 0.08;
2340 }
2341 if (d_b >= 0.41) {
2342 sigm_y2_ab[1][pair_y2] = 0.10 * 0.10;
2343 }
2344 pair_y2++;
2345 used_ya2[i] = true;
2346 used_yb2[j] = true;
2347
2348 } // j for2
2349 } // i for1 Y.
2350
2351 // reconstruct single Y-plane hits
2352
2353 for (Int_t i = 0; i < it_ya2; ++i) {
2354 if (used_ya2[i])
2355 continue;
2356 double d_a = 0;
2357 double d_b = 0;
2358
2359 int func_nr_a = -1;
2360 for (int t_it = 0; t_it < 4; t_it++) {
2361 if (time_ya2[i] >= t_dc[t_it][10] && time_ya2[i] < t_dc[t_it + 1][10]) {
2362 func_nr_a = t_it;
2363 break;
2364 }
2365 }
2366 if (func_nr_a == 1 || func_nr_a == 2)
2367 d_a = scale[10]
2368 * (pol_par_dc[1][0][10] + pol_par_dc[1][1][10] * time_ya2[i]
2369 + pol_par_dc[1][2][10] * time_ya2[i] * time_ya2[i]
2370 + pol_par_dc[1][3][10] * time_ya2[i] * time_ya2[i] * time_ya2[i]
2371 + pol_par_dc[1][4][10] * time_ya2[i] * time_ya2[i] * time_ya2[i] * time_ya2[i]);
2372
2373 if (func_nr_a == 0)
2374 d_a = 0;
2375 if (func_nr_a == 3)
2376 d_a = scale[10]
2377 * (pol_par_dc[2][0][10] + pol_par_dc[2][1][10] * time_ya2[i]
2378 + pol_par_dc[2][2][10] * time_ya2[i] * time_ya2[i]
2379 + pol_par_dc[2][3][10] * time_ya2[i] * time_ya2[i] * time_ya2[i]
2380 + pol_par_dc[2][4][10] * time_ya2[i] * time_ya2[i] * time_ya2[i] * time_ya2[i]);
2381
2382 y2_single[0][single_ya2] = wirenr_ya2[i] - 119 + d_a;
2383 y2_single[0][single_ya2 + 1] = wirenr_ya2[i] - 119 - d_a;
2384 if (d_a < 0.02) {
2385 sigm_y2_single[0][single_ya2] = 0.08 * 0.08;
2386 sigm_y2_single[0][single_ya2 + 1] = 0.08 * 0.08;
2387 }
2388 if (d_a >= 0.02 && d_a < 0.1) {
2389 sigm_y2_single[0][single_ya2] = 0.06 * 0.06;
2390 sigm_y2_single[0][single_ya2 + 1] = 0.06 * 0.06;
2391 }
2392 if (d_a >= 0.1 && d_a < 0.4) {
2393 sigm_y2_single[0][single_ya2] = 0.025 * 0.025;
2394 sigm_y2_single[0][single_ya2 + 1] = 0.025 * 0.025;
2395 }
2396 if (d_a >= 0.4 && d_a < 0.41) {
2397 sigm_y2_single[0][single_ya2] = 0.08 * 0.08;
2398 sigm_y2_single[0][single_ya2 + 1] = 0.08 * 0.08;
2399 }
2400 if (d_a >= 0.41) {
2401 sigm_y2_single[0][single_ya2] = 0.10 * 0.10;
2402 sigm_y2_single[0][single_ya2 + 1] = 0.10 * 0.10;
2403 }
2404 single_ya2 += 2;
2405 } // for single ya
2406
2407 for (Int_t j = 0; j < it_yb2; ++j) {
2408 if (used_yb2[j])
2409 continue;
2410 double d_a = 0;
2411 double d_b = 0;
2412 int func_nr_b = -1;
2413 for (int t_it = 0; t_it < 4; t_it++) {
2414 if (time_yb2[j] >= t_dc[t_it][11] && time_yb2[j] < t_dc[t_it + 1][11]) {
2415 func_nr_b = t_it;
2416 break;
2417 }
2418 }
2419 if (func_nr_b == 1 || func_nr_b == 2)
2420 d_b = scale[11]
2421 * (pol_par_dc[1][0][11] + pol_par_dc[1][1][11] * time_yb2[j]
2422 + pol_par_dc[1][2][11] * time_yb2[j] * time_yb2[j]
2423 + pol_par_dc[1][3][11] * time_yb2[j] * time_yb2[j] * time_yb2[j]
2424 + pol_par_dc[1][4][11] * time_yb2[j] * time_yb2[j] * time_yb2[j] * time_yb2[j]);
2425 if (func_nr_b == 0)
2426 d_b = 0;
2427 if (func_nr_b == 3)
2428 d_b = scale[11]
2429 * (pol_par_dc[2][0][11] + pol_par_dc[2][1][11] * time_yb2[j]
2430 + pol_par_dc[2][2][11] * time_yb2[j] * time_yb2[j]
2431 + pol_par_dc[2][3][11] * time_yb2[j] * time_yb2[j] * time_yb2[j]
2432 + pol_par_dc[2][4][11] * time_yb2[j] * time_yb2[j] * time_yb2[j] * time_yb2[j]);
2433
2434 y2_single[1][single_yb2] = wirenr_yb2[j] - 118.5 + d_b;
2435 y2_single[1][single_yb2 + 1] = wirenr_yb2[j] - 118.5 - d_b;
2436 if (d_b < 0.02) {
2437 sigm_y2_single[1][single_yb2] = 0.08 * 0.08;
2438 sigm_y2_single[1][single_yb2 + 1] = 0.08 * 0.08;
2439 }
2440 if (d_b >= 0.02 && d_b < 0.1) {
2441 sigm_y2_single[1][single_yb2] = 0.06 * 0.06;
2442 sigm_y2_single[1][single_yb2 + 1] = 0.06 * 0.06;
2443 }
2444 if (d_b >= 0.1 && d_b < 0.4) {
2445 sigm_y2_single[1][single_yb2] = 0.0250 * 0.0250;
2446 sigm_y2_single[1][single_yb2 + 1] = 0.0250 * 0.0250;
2447 }
2448 if (d_b >= 0.4 && d_b < 0.41) {
2449 sigm_y2_single[1][single_yb2] = 0.1 * 0.1;
2450 sigm_y2_single[1][single_yb2 + 1] = 0.1 * 0.1;
2451 }
2452 if (d_b >= 0.41) {
2453 sigm_y2_single[1][single_yb2] = 0.10 * 0.10;
2454 sigm_y2_single[1][single_yb2 + 1] = 0.10 * 0.10;
2455 }
2456 single_yb2 += 2;
2457
2458 } // i for1 Y.
2459
2460 // ---- U ---
2461
2462 for (Int_t i = 0; i < it_ua2; ++i) {
2463 for (Int_t j = 0; j < it_ub2; ++j) {
2464 if (pair_u2 > 48)
2465 break;
2466 if ((wirenr_ua2[i] != wirenr_ub2[j] && wirenr_ua2[i] != wirenr_ub2[j] + 1))
2467 continue;
2468
2469 double d_a = 0;
2470 double d_b = 0;
2471
2472 int func_nr_a = -1;
2473 int func_nr_b = -1;
2474 for (int t_it = 0; t_it < 4; t_it++) {
2475 if (time_ua2[i] >= t_dc[t_it][12] && time_ua2[i] < t_dc[t_it + 1][12]) {
2476 func_nr_a = t_it;
2477 break;
2478 }
2479 }
2480 if (func_nr_a == 1 || func_nr_a == 2)
2481 d_a = scale[12]
2482 * (pol_par_dc[1][0][12] + pol_par_dc[1][1][12] * time_ua2[i]
2483 + pol_par_dc[1][2][12] * time_ua2[i] * time_ua2[i]
2484 + pol_par_dc[1][3][12] * time_ua2[i] * time_ua2[i] * time_ua2[i]
2485 + pol_par_dc[1][4][12] * time_ua2[i] * time_ua2[i] * time_ua2[i] * time_ua2[i]);
2486
2487 if (func_nr_a == 0)
2488 d_a = 0;
2489 if (func_nr_a == 3)
2490 d_a = scale[12]
2491 * (pol_par_dc[1][0][12] + pol_par_dc[1][1][12] * time_ua2[i]
2492 + pol_par_dc[1][2][12] * time_ua2[i] * time_ua2[i]
2493 + pol_par_dc[1][3][12] * time_ua2[i] * time_ua2[i] * time_ua2[i]
2494 + pol_par_dc[1][4][12] * time_ua2[i] * time_ua2[i] * time_ua2[i] * time_ua2[i]);
2495
2496 for (int t_it = 0; t_it < 4; t_it++) {
2497 if (time_ub2[j] >= t_dc[t_it][13] && time_ub2[j] < t_dc[t_it + 1][13]) {
2498 func_nr_b = t_it;
2499 break;
2500 }
2501 }
2502
2503 if (func_nr_b == 0)
2504 d_b = 0;
2505
2506 if (func_nr_b == 1 || func_nr_b == 2)
2507 d_b = scale[13]
2508 * (pol_par_dc[1][0][13] + pol_par_dc[1][1][13] * time_ub2[j]
2509 + pol_par_dc[1][2][13] * time_ub2[j] * time_ub2[j]
2510 + pol_par_dc[1][3][13] * time_ub2[j] * time_ub2[j] * time_ub2[j]
2511 + pol_par_dc[1][4][13] * time_ub2[j] * time_ub2[j] * time_ub2[j] * time_ub2[j]);
2512
2513 if (func_nr_b == 3)
2514 d_b = scale[13]
2515 * (pol_par_dc[1][0][13] + pol_par_dc[1][1][13] * time_ub2[j]
2516 + pol_par_dc[1][2][13] * time_ub2[j] * time_ub2[j]
2517 + pol_par_dc[1][3][13] * time_ub2[j] * time_ub2[j] * time_ub2[j]
2518 + pol_par_dc[1][4][13] * time_ub2[j] * time_ub2[j] * time_ub2[j] * time_ub2[j]);
2519
2520 ua2_pm[0] = wirenr_ua2[i] - 119 + d_a;
2521 ua2_pm[1] = wirenr_ua2[i] - 119 - d_a;
2522 ub2_pm[0] = wirenr_ub2[j] - 118.5 + d_b;
2523 ub2_pm[1] = wirenr_ub2[j] - 118.5 - d_b;
2524
2525 double dmin1 = 999;
2526 for (int k = 0; k < 2; k++) {
2527 for (int m = 0; m < 2; m++) {
2528 if (fabs(ua2_pm[k] - ub2_pm[m]) < dmin1) {
2529 dmin1 = fabs(ua2_pm[k] - ub2_pm[m]);
2530 u2_ab[0][pair_u2] = ua2_pm[k];
2531 u2_ab[1][pair_u2] = ub2_pm[m];
2532 }
2533 }
2534 }
2535 if (d_a < 0.02) {
2536 sigm_u2_ab[0][pair_u2] = 0.08 * 0.08;
2537 }
2538 if (d_a >= 0.02 && d_a < 0.1) {
2539 sigm_u2_ab[0][pair_u2] = 0.06 * 0.06;
2540 }
2541 if (d_a >= 0.1 && d_a < 0.4) {
2542 sigm_u2_ab[0][pair_u2] = 0.025 * 0.025;
2543 }
2544 if (d_a >= 0.4 && d_a < 0.41) {
2545 sigm_u2_ab[0][pair_u2] = 0.08 * 0.08;
2546 }
2547 if (d_a >= 0.41) {
2548 sigm_u2_ab[0][pair_u2] = 0.10 * 0.10;
2549 }
2550
2551 if (d_b < 0.02) {
2552 sigm_u2_ab[1][pair_u2] = 0.08 * 0.08;
2553 }
2554 if (d_b >= 0.02 && d_b < 0.1) {
2555 sigm_u2_ab[1][pair_u2] = 0.06 * 0.06;
2556 }
2557 if (d_b >= 0.1 && d_b < 0.4) {
2558 sigm_u2_ab[1][pair_u2] = 0.025 * 0.025;
2559 }
2560 if (d_b >= 0.4 && d_b < 0.41) {
2561 sigm_u2_ab[1][pair_u2] = 0.08 * 0.08;
2562 }
2563 if (d_b >= 0.41) {
2564 sigm_u2_ab[1][pair_u2] = 0.10 * 0.10;
2565 }
2566 pair_u2++;
2567 used_ua2[i] = true;
2568 used_ub2[j] = true;
2569
2570 } // j for2
2571 } // i for1 U.
2572
2573 // reconstruct single U-plane hits
2574
2575 for (Int_t i = 0; i < it_ua2; ++i) {
2576 if (used_ua2[i])
2577 continue;
2578 double d_a = 0;
2579 double d_b = 0;
2580
2581 int func_nr_a = -1;
2582 for (int t_it = 0; t_it < 4; t_it++) {
2583 if (time_ua2[i] >= t_dc[t_it][12] && time_ua2[i] < t_dc[t_it + 1][12]) {
2584 func_nr_a = t_it;
2585 break;
2586 }
2587 }
2588 if (func_nr_a == 1 || func_nr_a == 2)
2589 d_a = scale[12]
2590 * (pol_par_dc[1][0][12] + pol_par_dc[1][1][12] * time_ua2[i]
2591 + pol_par_dc[1][2][12] * time_ua2[i] * time_ua2[i]
2592 + pol_par_dc[1][3][12] * time_ua2[i] * time_ua2[i] * time_ua2[i]
2593 + pol_par_dc[1][4][12] * time_ua2[i] * time_ua2[i] * time_ua2[i] * time_ua2[i]);
2594 if (func_nr_a == 0)
2595 d_a = 0;
2596 if (func_nr_a == 3)
2597 d_a = scale[12]
2598 * (pol_par_dc[1][0][12] + pol_par_dc[1][1][12] * time_ua2[i]
2599 + pol_par_dc[1][2][12] * time_ua2[i] * time_ua2[i]
2600 + pol_par_dc[1][3][12] * time_ua2[i] * time_ua2[i] * time_ua2[i]
2601 + pol_par_dc[1][4][12] * time_ua2[i] * time_ua2[i] * time_ua2[i] * time_ua2[i]);
2602
2603 u2_single[0][single_ua2] = (wirenr_ua2[i] - 119 + d_a);
2604 u2_single[0][single_ua2 + 1] = (wirenr_ua2[i] - 119 - d_a);
2605 if (d_a < 0.02) {
2606 sigm_u2_single[0][single_ua2] = 0.08 * 0.08;
2607 sigm_u2_single[0][single_ua2 + 1] = 0.08 * 0.08;
2608 }
2609 if (d_a >= 0.02 && d_a < 0.1) {
2610 sigm_u2_single[0][single_ua2] = 0.06 * 0.06;
2611 sigm_u2_single[0][single_ua2 + 1] = 0.06 * 0.06;
2612 }
2613 if (d_a >= 0.1 && d_a < 0.4) {
2614 sigm_u2_single[0][single_ua2] = 0.025 * 0.025;
2615 sigm_u2_single[0][single_ua2 + 1] = 0.025 * 0.025;
2616 }
2617 if (d_a >= 0.4 && d_a < 0.41) {
2618 sigm_u2_single[0][single_ua2] = 0.08 * 0.08;
2619 sigm_u2_single[0][single_ua2 + 1] = 0.08 * 0.08;
2620 }
2621 if (d_a >= 0.41) {
2622 sigm_u2_single[0][single_ua2] = 0.10 * 0.10;
2623 sigm_u2_single[0][single_ua2 + 1] = 0.10 * 0.10;
2624 }
2625 single_ua2 += 2;
2626 } // for single ua
2627
2628 for (Int_t j = 0; j < it_ub2; ++j) {
2629 if (used_ub2[j])
2630 continue;
2631 int func_nr_b = -1;
2632 double d_a = 0;
2633 double d_b = 0;
2634
2635 for (int t_it = 0; t_it < 4; t_it++) {
2636 if (time_ub2[j] >= t_dc[t_it][13] && time_ub2[j] < t_dc[t_it + 1][13]) {
2637 func_nr_b = t_it;
2638 break;
2639 }
2640 }
2641 if (func_nr_b == 0)
2642 d_b = 0;
2643
2644 if (func_nr_b == 1 || func_nr_b == 2)
2645 d_b = scale[13]
2646 * (pol_par_dc[1][0][13] + pol_par_dc[1][1][13] * time_ub2[j]
2647 + pol_par_dc[1][2][13] * time_ub2[j] * time_ub2[j]
2648 + pol_par_dc[1][3][13] * time_ub2[j] * time_ub2[j] * time_ub2[j]
2649 + pol_par_dc[1][4][13] * time_ub2[j] * time_ub2[j] * time_ub2[j] * time_ub2[j]);
2650 if (func_nr_b == 3)
2651 d_b = scale[13]
2652 * (pol_par_dc[1][0][13] + pol_par_dc[1][1][13] * time_ub2[j]
2653 + pol_par_dc[1][2][13] * time_ub2[j] * time_ub2[j]
2654 + pol_par_dc[1][3][13] * time_ub2[j] * time_ub2[j] * time_ub2[j]
2655 + pol_par_dc[1][4][13] * time_ub2[j] * time_ub2[j] * time_ub2[j] * time_ub2[j]);
2656
2657 u2_single[1][single_ub2] = (wirenr_ub2[j] - 118.5 + d_b);
2658 u2_single[1][single_ub2 + 1] = (wirenr_ub2[j] - 118.5 - d_b);
2659 if (d_b < 0.02) {
2660 sigm_u2_single[1][single_ub2] = 0.08 * 0.08;
2661 sigm_u2_single[1][single_ub2 + 1] = 0.08 * 0.08;
2662 }
2663 if (d_b >= 0.02 && d_b < 0.1) {
2664 sigm_u2_single[1][single_ub2] = 0.060 * 0.060;
2665 sigm_u2_single[1][single_ub2 + 1] = 0.060 * 0.060;
2666 }
2667 if (d_b >= 0.1 && d_b < 0.4) {
2668 sigm_u2_single[1][single_ub2] = 0.025 * 0.025;
2669 sigm_u2_single[1][single_ub2 + 1] = 0.0250 * 0.0250;
2670 }
2671 if (d_b >= 0.4 && d_b < 0.41) {
2672 sigm_u2_single[1][single_ub2] = 0.08 * 0.08;
2673 sigm_u2_single[1][single_ub2 + 1] = 0.08 * 0.08;
2674 }
2675 if (d_b >= 0.41) {
2676 sigm_u2_single[1][single_ub2] = 0.10 * 0.10;
2677 sigm_u2_single[1][single_ub2 + 1] = 0.10 * 0.10;
2678 }
2679 single_ub2 += 2;
2680
2681 } // i for1 U.
2682
2683 // ---- V ---
2684
2685 for (Int_t i = 0; i < it_va2; ++i) {
2686 for (Int_t j = 0; j < it_vb2; ++j) {
2687 if (pair_v2 > 48)
2688 break;
2689 if ((wirenr_va2[i] != wirenr_vb2[j] && wirenr_va2[i] != wirenr_vb2[j] + 1))
2690 continue;
2691
2692 double d_a = 0;
2693 double d_b = 0;
2694
2695 int func_nr_a = -1;
2696 int func_nr_b = -1;
2697 for (int t_it = 0; t_it < 4; t_it++) {
2698 if (time_va2[i] >= t_dc[t_it][14] && time_va2[i] < t_dc[t_it + 1][14]) {
2699 func_nr_a = t_it;
2700 break;
2701 }
2702 }
2703 if (func_nr_a == 1 || func_nr_a == 2)
2704 d_a = scale[14]
2705 * (pol_par_dc[1][0][14] + pol_par_dc[1][1][14] * time_va2[i]
2706 + pol_par_dc[1][2][14] * time_va2[i] * time_va2[i]
2707 + pol_par_dc[1][3][14] * time_va2[i] * time_va2[i] * time_va2[i]
2708 + pol_par_dc[1][4][14] * time_va2[i] * time_va2[i] * time_va2[i] * time_va2[i]);
2709
2710 if (func_nr_a == 0)
2711 d_a = 0;
2712 if (func_nr_a == 3)
2713 d_a = scale[14]
2714 * (pol_par_dc[2][0][14] + pol_par_dc[2][1][14] * time_va2[i]
2715 + pol_par_dc[2][2][14] * time_va2[i] * time_va2[i]
2716 + pol_par_dc[2][3][14] * time_va2[i] * time_va2[i] * time_va2[i]
2717 + pol_par_dc[2][4][14] * time_va2[i] * time_va2[i] * time_va2[i] * time_va2[i]);
2718
2719 for (int t_it = 0; t_it < 4; t_it++) {
2720 if (time_vb2[j] >= t_dc[t_it][15] && time_vb2[j] < t_dc[t_it + 1][15]) {
2721 func_nr_b = t_it;
2722 break;
2723 }
2724 }
2725 if (func_nr_b == 1 || func_nr_b == 2)
2726 d_b = scale[15]
2727 * (pol_par_dc[1][0][15] + pol_par_dc[1][1][15] * time_vb2[j]
2728 + pol_par_dc[1][2][15] * time_vb2[j] * time_vb2[j]
2729 + pol_par_dc[1][3][15] * time_vb2[j] * time_vb2[j] * time_vb2[j]
2730 + pol_par_dc[1][4][15] * time_vb2[j] * time_vb2[j] * time_vb2[j] * time_vb2[j]);
2731
2732 if (func_nr_b == 0)
2733 d_b = 0;
2734 if (func_nr_b == 3)
2735 d_b = scale[15]
2736 * (pol_par_dc[2][0][15] + pol_par_dc[2][1][15] * time_vb2[j]
2737 + pol_par_dc[2][2][15] * time_vb2[j] * time_vb2[j]
2738 + pol_par_dc[2][3][15] * time_vb2[j] * time_vb2[j] * time_vb2[j]
2739 + pol_par_dc[2][4][15] * time_vb2[j] * time_vb2[j] * time_vb2[j] * time_vb2[j]);
2740
2741 va2_pm[0] = wirenr_va2[i] - 119 + d_a;
2742 va2_pm[1] = wirenr_va2[i] - 119 - d_a;
2743 vb2_pm[0] = wirenr_vb2[j] - 118.5 + d_b;
2744 vb2_pm[1] = wirenr_vb2[j] - 118.5 - d_b;
2745
2746 double dmin1 = 999;
2747 for (int k = 0; k < 2; k++) {
2748 for (int m = 0; m < 2; m++) {
2749 if (fabs(va2_pm[k] - vb2_pm[m]) < dmin1) {
2750 dmin1 = fabs(va2_pm[k] - vb2_pm[m]);
2751 v2_ab[0][pair_v2] = va2_pm[k];
2752 v2_ab[1][pair_v2] = vb2_pm[m];
2753 }
2754 }
2755 }
2756 if (d_a < 0.02) {
2757 sigm_v2_ab[0][pair_v2] = 0.08 * 0.08;
2758 }
2759 if (d_a >= 0.02 && d_a < 0.1) {
2760 sigm_v2_ab[0][pair_v2] = 0.06 * 0.06;
2761 }
2762 if (d_a >= 0.1 && d_a < 0.4) {
2763 sigm_v2_ab[0][pair_v2] = 0.025 * 0.025;
2764 }
2765 if (d_a >= 0.4 && d_a < 0.41) {
2766 sigm_v2_ab[0][pair_v2] = 0.08 * 0.08;
2767 }
2768 if (d_a >= 0.41) {
2769 sigm_v2_ab[0][pair_v2] = 0.10 * 0.10;
2770 }
2771
2772 if (d_b < 0.02) {
2773 sigm_v2_ab[1][pair_v2] = 0.08 * 0.08;
2774 }
2775 if (d_b >= 0.02 && d_b < 0.1) {
2776 sigm_v2_ab[1][pair_v2] = 0.06 * 0.06;
2777 }
2778 if (d_b >= 0.1 && d_b < 0.4) {
2779 sigm_v2_ab[1][pair_v2] = 0.025 * 0.025;
2780 }
2781 if (d_b >= 0.4 && d_b < 0.41) {
2782 sigm_v2_ab[1][pair_v2] = 0.08 * 0.08;
2783 }
2784 if (d_b >= 0.41) {
2785 sigm_v2_ab[1][pair_v2] = 0.10 * 0.10;
2786 }
2787 pair_v2++;
2788
2789 used_va2[i] = true;
2790 used_vb2[j] = true;
2791
2792 } // j for2
2793 } // i for1 V.
2794
2795 // reconstruct single V-plane hits
2796
2797 for (Int_t i = 0; i < it_va2; ++i) {
2798 if (used_va2[i])
2799 continue;
2800 double d_a = 0;
2801 double d_b = 0;
2802
2803 // added
2804 int func_nr_a = -1;
2805 for (int t_it = 0; t_it < 4; t_it++) {
2806 if (time_va2[i] >= t_dc[t_it][14] && time_va2[i] < t_dc[t_it + 1][14]) {
2807 func_nr_a = t_it;
2808 break;
2809 }
2810 }
2811 if (func_nr_a == 1 || func_nr_a == 2)
2812 d_a = scale[14]
2813 * (pol_par_dc[1][0][14] + pol_par_dc[1][1][14] * time_va2[i]
2814 + pol_par_dc[1][2][14] * time_va2[i] * time_va2[i]
2815 + pol_par_dc[1][3][14] * time_va2[i] * time_va2[i] * time_va2[i]
2816 + pol_par_dc[1][4][14] * time_va2[i] * time_va2[i] * time_va2[i] * time_va2[i]);
2817 if (func_nr_a == 0)
2818 d_a = 0;
2819 if (func_nr_a == 3)
2820 d_a = scale[14]
2821 * (pol_par_dc[2][0][14] + pol_par_dc[2][1][14] * time_va2[i]
2822 + pol_par_dc[2][2][14] * time_va2[i] * time_va2[i]
2823 + pol_par_dc[2][3][14] * time_va2[i] * time_va2[i] * time_va2[i]
2824 + pol_par_dc[2][4][14] * time_va2[i] * time_va2[i] * time_va2[i] * time_va2[i]);
2825
2826 v2_single[0][single_va2] = wirenr_va2[i] - 119 + d_a;
2827 v2_single[0][single_va2 + 1] = wirenr_va2[i] - 119 - d_a;
2828 if (d_a < 0.02) {
2829 sigm_v2_single[0][single_va2] = 0.08 * 0.08;
2830 sigm_v2_single[0][single_va2 + 1] = 0.08 * 0.08;
2831 }
2832 if (d_a >= 0.02 && d_a < 0.1) {
2833 sigm_v2_single[0][single_va2] = 0.06 * 0.06;
2834 sigm_v2_single[0][single_va2 + 1] = 0.06 * 0.06;
2835 }
2836 if (d_a >= 0.1 && d_a < 0.4) {
2837 sigm_v2_single[0][single_va2] = 0.025 * 0.025;
2838 sigm_v2_single[0][single_va2 + 1] = 0.025 * 0.025;
2839 }
2840 if (d_a >= 0.4 && d_a < 0.41) {
2841 sigm_v2_single[0][single_va2] = 0.08 * 0.08;
2842 sigm_v2_single[0][single_va2 + 1] = 0.08 * 0.08;
2843 }
2844 if (d_a >= 0.41) {
2845 sigm_v2_single[0][single_va2] = 0.10 * 0.10;
2846 sigm_v2_single[0][single_va2 + 1] = 0.10 * 0.10;
2847 }
2848
2849 single_va2 += 2;
2850 } // for single va
2851
2852 for (Int_t j = 0; j < it_vb2; ++j) {
2853 if (used_vb2[j])
2854 continue;
2855 double d_a = 0;
2856 double d_b = 0;
2857 int func_nr_b = -1;
2858 for (int t_it = 0; t_it < 4; t_it++) {
2859 if (time_vb2[j] >= t_dc[t_it][15] && time_vb2[j] < t_dc[t_it + 1][15]) {
2860 func_nr_b = t_it;
2861 break;
2862 }
2863 }
2864 if (func_nr_b == 1 || func_nr_b == 2)
2865 d_b = scale[15]
2866 * (pol_par_dc[1][0][15] + pol_par_dc[1][1][15] * time_vb2[j]
2867 + pol_par_dc[1][2][15] * time_vb2[j] * time_vb2[j]
2868 + pol_par_dc[1][3][15] * time_vb2[j] * time_vb2[j] * time_vb2[j]
2869 + pol_par_dc[1][4][15] * time_vb2[j] * time_vb2[j] * time_vb2[j] * time_vb2[j]);
2870 if (func_nr_b == 0)
2871 d_b = 0;
2872 if (func_nr_b == 3)
2873 d_b = scale[15]
2874 * (pol_par_dc[2][0][15] + pol_par_dc[2][1][15] * time_vb2[j]
2875 + pol_par_dc[2][2][15] * time_vb2[j] * time_vb2[j]
2876 + pol_par_dc[2][3][15] * time_vb2[j] * time_vb2[j] * time_vb2[j]
2877 + pol_par_dc[2][4][15] * time_vb2[j] * time_vb2[j] * time_vb2[j] * time_vb2[j]);
2878
2879 v2_single[1][single_vb2] = wirenr_vb2[j] - 118.5 + d_b;
2880 v2_single[1][single_vb2 + 1] = wirenr_vb2[j] - 118.5 - d_b;
2881 if (d_b < 0.02) {
2882 sigm_v2_single[1][single_vb2] = 0.08 * 0.08;
2883 sigm_v2_single[1][single_vb2 + 1] = 0.08 * 0.08;
2884 }
2885 if (d_b >= 0.02 && d_b < 0.1) {
2886 sigm_v2_single[1][single_vb2] = 0.060 * 0.060;
2887 sigm_v2_single[1][single_vb2 + 1] = 0.060 * 0.060;
2888 }
2889 if (d_b >= 0.1 && d_b < 0.4) {
2890 sigm_v2_single[1][single_vb2] = 0.0250 * 0.0250;
2891 sigm_v2_single[1][single_vb2 + 1] = 0.0250 * 0.0250;
2892 }
2893 if (d_b >= 0.4 && d_b < 0.41) {
2894 sigm_v2_single[1][single_vb2] = 0.08 * 0.08;
2895 sigm_v2_single[1][single_vb2 + 1] = 0.08 * 0.08;
2896 }
2897 if (d_b >= 0.41) {
2898 sigm_v2_single[1][single_vb2] = 0.10 * 0.10;
2899 sigm_v2_single[1][single_vb2 + 1] = 0.10 * 0.10;
2900 }
2901 single_vb2 += 2;
2902 } // i for1 V.
2903
2904 // --- DC2 Segment Building ------
2905
2906 found_8p_seg = false;
2907 x_slope = 0.0;
2908 y_slope = 0.0;
2909
2910 if (pair_x2 * pair_y2 > 0) { //(x,y)first
2911
2912 for (int i = 0; i < pair_x2; i++) {
2913 if (nDC2_segments > 48)
2914 break;
2915 float x_coord = (x2_ab[0][i] + x2_ab[1][i]) / 2;
2916 float XU = x_coord + x_slope * 0.5 * (z_loc[5] + z_loc[4] - z_loc[1] - z_loc[0]);
2917 float XV = x_coord + x_slope * 0.5 * (z_loc[7] + z_loc[6] - z_loc[1] - z_loc[0]);
2918
2919 for (int j = 0; j < pair_y2; j++) {
2920 if (nDC2_segments > 48)
2921 break;
2922 float y_coord = (y2_ab[0][j] + y2_ab[1][j]) / 2;
2923 float YU = y_coord + y_slope * 0.5 * (z_loc[5] + z_loc[4] - z_loc[3] - z_loc[2]);
2924 float YV = y_coord + y_slope * 0.5 * (z_loc[7] + z_loc[6] - z_loc[3] - z_loc[2]);
2925 bool foundU = false;
2926 float u_est = isqrt_2 * (YU - XU);
2927 float v_est = isqrt_2 * (YV + XV);
2928 if (pair_u2 > 0) {
2929
2930 double dU_thresh = 1.3;
2931 for (int k = 0; k < pair_u2; k++) {
2932 float u_coord = (u2_ab[0][k] + u2_ab[1][k]) / 2;
2933
2934 if (fabs(u_coord - u_est) > dU_thresh)
2935 continue;
2936 dU_thresh = fabs(u_coord - u_est);
2937
2938 rh_segDC2[0][nDC2_segments] = x2_ab[0][i];
2939 rh_segDC2[1][nDC2_segments] = x2_ab[1][i];
2940 rh_segDC2[2][nDC2_segments] = y2_ab[0][j];
2941 rh_segDC2[3][nDC2_segments] = y2_ab[1][j];
2942 rh_segDC2[4][nDC2_segments] = u2_ab[0][k];
2943 rh_segDC2[5][nDC2_segments] = u2_ab[1][k];
2944 rh_sigm_segDC2[0][nDC2_segments] = sigm_x2_ab[0][i];
2945 rh_sigm_segDC2[1][nDC2_segments] = sigm_x2_ab[1][i];
2946 rh_sigm_segDC2[2][nDC2_segments] = sigm_y2_ab[0][j];
2947 rh_sigm_segDC2[3][nDC2_segments] = sigm_y2_ab[1][j];
2948 rh_sigm_segDC2[4][nDC2_segments] = sigm_u2_ab[0][k];
2949 rh_sigm_segDC2[5][nDC2_segments] = sigm_u2_ab[1][k];
2950
2951 foundU = true;
2952 if (nDC2_segments > 48)
2953 break;
2954 } // k
2955 } //(pair_u2>0)
2956 if (found_8p_seg && !foundU)
2957 continue;
2958 bool foundV = false;
2959 if (pair_v2 > 0) {
2960 double dV_thresh = 1.3;
2961 for (int m = 0; m < pair_v2; m++) {
2962 if (nDC2_segments > 48)
2963 break;
2964 float v_coord = (v2_ab[0][m] + v2_ab[1][m]) / 2;
2965
2966 if (fabs(v_coord - v_est) > dV_thresh)
2967 continue;
2968 dV_thresh = fabs(v_coord - v_est);
2969
2970 foundV = true;
2971 rh_segDC2[0][nDC2_segments] = x2_ab[0][i];
2972 rh_segDC2[1][nDC2_segments] = x2_ab[1][i];
2973 rh_segDC2[2][nDC2_segments] = y2_ab[0][j];
2974 rh_segDC2[3][nDC2_segments] = y2_ab[1][j];
2975 rh_segDC2[6][nDC2_segments] = v2_ab[0][m];
2976 rh_segDC2[7][nDC2_segments] = v2_ab[1][m];
2977 rh_sigm_segDC2[0][nDC2_segments] = sigm_x2_ab[0][i];
2978 rh_sigm_segDC2[1][nDC2_segments] = sigm_x2_ab[1][i];
2979 rh_sigm_segDC2[2][nDC2_segments] = sigm_y2_ab[0][j];
2980 rh_sigm_segDC2[3][nDC2_segments] = sigm_y2_ab[1][j];
2981 rh_sigm_segDC2[6][nDC2_segments] = sigm_v2_ab[0][m];
2982 rh_sigm_segDC2[7][nDC2_segments] = sigm_v2_ab[1][m];
2983 if (!foundU) {
2984 float min_a = 999;
2985 float min_b = 999;
2986 for (int kk = 0; kk < single_ua2; kk++) {
2987 if (fabs(u2_single[0][kk] - u_est) > 1.5)
2988 continue; //????? 0.5 needs to be reviewed
2989 if (fabs(u2_single[0][kk] - u_est) < min_a) {
2990 min_a = fabs(u2_single[0][kk] - u_est);
2991 rh_segDC2[4][nDC2_segments] = u2_single[0][kk];
2992 rh_sigm_segDC2[4][nDC2_segments] = sigm_u2_single[0][kk];
2993 foundU = true;
2994 }
2995 } // for kk
2996 for (int kk = 0; kk < single_ub2; kk++) {
2997 if (fabs(u2_single[1][kk] - u_est) > 1.5)
2998 continue; //????? 0.5 needs to be reviewed
2999 if (fabs(u2_single[1][kk] - u_est) < min_b) {
3000 min_b = fabs(u2_single[1][kk] - u_est);
3001 rh_segDC2[5][nDC2_segments] = u2_single[1][kk];
3002 rh_sigm_segDC2[5][nDC2_segments] = sigm_u2_single[1][kk];
3003 foundU = true;
3004 }
3005 } // for kk
3006 if (nDC2_segments > 48)
3007 break;
3008 } //! foundU
3009
3010 if (nDC2_segments > 48)
3011 break;
3012 } // m
3013 } //(pair_v2>0)
3014 if (!foundV && foundU) {
3015 float min_a = 999;
3016 float min_b = 999;
3017 for (int kk = 0; kk < single_va2; kk++) {
3018 if (fabs(v2_single[0][kk] - v_est) > 1.5)
3019 continue; //????? 0.5 needs to be reviewed
3020 if (fabs(v2_single[0][kk] - v_est) < min_a) {
3021 min_a = fabs(v2_single[0][kk] - v_est);
3022 rh_segDC2[6][nDC2_segments] = v2_single[0][kk];
3023 rh_sigm_segDC2[6][nDC2_segments] = sigm_v2_single[0][kk];
3024 foundV = true;
3025 }
3026 } // for kk
3027 for (int kk = 0; kk < single_vb2; kk++) {
3028 if (fabs(v2_single[1][kk] - v_est) > 1.5)
3029 continue; //????? 0.5 needs to be reviewed
3030 if (fabs(v2_single[1][kk] - v_est) < min_b) {
3031 min_b = fabs(v2_single[1][kk] - v_est);
3032 rh_segDC2[7][nDC2_segments] = v2_single[1][kk];
3033 rh_sigm_segDC2[7][nDC2_segments] = sigm_v2_single[1][kk];
3034 foundV = true;
3035 }
3036 } // for kk
3037 } //! foundV
3038 if (foundV || foundU)
3039 nDC2_segments++;
3040 } // j
3041 } // i
3042 } //(x,y)first
3043
3044 if (pair_u2 * pair_v2 > 0 && !found_8p_seg) { // (u,v) first
3045 for (int i = 0; i < pair_u2; i++) {
3046 if (nDC2_segments > 48)
3047 break;
3048 float u_coord = (u2_ab[0][i] + u2_ab[1][i]) / 2;
3049
3050 for (int j = 0; j < pair_v2; j++) {
3051 if (nDC2_segments > 48)
3052 break;
3053 float v_coord = (v2_ab[0][j] + v2_ab[1][j]) / 2;
3054 bool foundX = false;
3055 float x_est = isqrt_2 * (v_coord - u_coord);
3056 float y_est = isqrt_2 * (u_coord + v_coord);
3057 if (pair_x2 > 0) {
3058 double dX_thresh = 1.5;
3059 for (int k = 0; k < pair_x2; k++) {
3060 float x_coord = (x2_ab[0][k] + x2_ab[1][k]) / 2;
3061 if (nDC2_segments > 48)
3062 break;
3063 if (fabs(x_coord - x_est) > dX_thresh)
3064 continue;
3065 dX_thresh = fabs(x_coord - x_est);
3066
3067 rh_segDC2[0][nDC2_segments] = x2_ab[0][k];
3068 rh_segDC2[1][nDC2_segments] = x2_ab[1][k];
3069 rh_segDC2[4][nDC2_segments] = u2_ab[0][i];
3070 rh_segDC2[5][nDC2_segments] = u2_ab[1][i];
3071 rh_segDC2[6][nDC2_segments] = v2_ab[0][j];
3072 rh_segDC2[7][nDC2_segments] = v2_ab[1][j];
3073 rh_sigm_segDC2[0][nDC2_segments] = sigm_x2_ab[0][k];
3074 rh_sigm_segDC2[1][nDC2_segments] = sigm_x2_ab[1][k];
3075 rh_sigm_segDC2[4][nDC2_segments] = sigm_u2_ab[0][i];
3076 rh_sigm_segDC2[5][nDC2_segments] = sigm_u2_ab[1][i];
3077 rh_sigm_segDC2[6][nDC2_segments] = sigm_v2_ab[0][j];
3078 rh_sigm_segDC2[7][nDC2_segments] = sigm_v2_ab[1][j];
3079
3080 foundX = true;
3081 if (nDC2_segments > 48)
3082 break;
3083 } // k
3084 } //(pair_x2>0)
3085 if (found_8p_seg && !foundX)
3086 continue;
3087 bool foundY = false;
3088 if (pair_y2 > 0) {
3089 double dY_thresh = 1.0;
3090 for (int m = 0; m < pair_y2; m++) {
3091 if (nDC2_segments > 48)
3092 break;
3093 float y_coord = (y2_ab[0][m] + y2_ab[1][m]) / 2;
3094 if (fabs(y_coord - y_est) > dY_thresh)
3095 continue;
3096 dY_thresh = fabs(y_coord - y_est);
3097 foundY = true;
3098 rh_segDC2[2][nDC2_segments] = y2_ab[0][m];
3099 rh_segDC2[3][nDC2_segments] = y2_ab[1][m];
3100 rh_segDC2[4][nDC2_segments] = u2_ab[0][i];
3101 rh_segDC2[5][nDC2_segments] = u2_ab[1][i];
3102 rh_segDC2[6][nDC2_segments] = v2_ab[0][j];
3103 rh_segDC2[7][nDC2_segments] = v2_ab[1][j];
3104 rh_sigm_segDC2[2][nDC2_segments] = sigm_y2_ab[0][m];
3105 rh_sigm_segDC2[3][nDC2_segments] = sigm_y2_ab[1][m];
3106 rh_sigm_segDC2[4][nDC2_segments] = sigm_u2_ab[0][i];
3107 rh_sigm_segDC2[5][nDC2_segments] = sigm_u2_ab[1][i];
3108 rh_sigm_segDC2[6][nDC2_segments] = sigm_v2_ab[0][j];
3109 rh_sigm_segDC2[7][nDC2_segments] = sigm_v2_ab[1][j];
3110 if (!foundX) {
3111 float min_a = 999;
3112 float min_b = 999;
3113 for (int kk = 0; kk < single_xa2; kk++) {
3114 if (fabs(x2_single[1][kk] - x_est) > 1.5)
3115 continue; //????? 0.5 needs to be reviewed
3116
3117 if (fabs(x2_single[0][kk] - x_est) < min_a) {
3118 min_a = fabs(x2_single[0][kk] - x_est);
3119 rh_segDC2[0][nDC2_segments] = x2_single[0][kk];
3120 rh_sigm_segDC2[0][nDC2_segments] = sigm_x2_single[0][kk];
3121 foundX = true;
3122 }
3123 } // for kk
3124 for (int kk = 0; kk < single_xb2; kk++) {
3125
3126 if (fabs(x2_single[1][kk] - x_est) > 1.5)
3127 continue; //????? 0.5 needs to be reviewed
3128 if (fabs(x2_single[1][kk] - x_est) < min_b) {
3129 min_b = fabs(x2_single[1][kk] - x_est);
3130 rh_segDC2[1][nDC2_segments] = x2_single[1][kk];
3131 rh_sigm_segDC2[1][nDC2_segments] = sigm_x2_single[1][kk];
3132 foundX = true;
3133 }
3134 } // for kk
3135 if (nDC2_segments > 48)
3136 break;
3137 } //! foundX
3138 } // m
3139 if (foundX && !foundY) {
3140 float min_a = 999;
3141 float min_b = 999;
3142 for (int kk = 0; kk < single_ya2; kk++) {
3143 if (fabs(y2_single[0][kk] - y_est) > 1.5)
3144 continue; //????? 0.5 needs to be reviewed
3145 if (fabs(y2_single[0][kk] - y_est) < min_a) {
3146 min_a = fabs(y2_single[0][kk] - y_est);
3147 rh_segDC2[2][nDC2_segments] = y2_single[0][kk];
3148 rh_sigm_segDC2[2][nDC2_segments] = sigm_y2_single[0][kk];
3149 foundY = true;
3150 }
3151 } // for kk
3152 for (int kk = 0; kk < single_yb2; kk++) {
3153 if (fabs(y2_single[1][kk] - y_est) > 1.5)
3154 continue; //????? 0.5 needs to be reviewed
3155 if (fabs(y2_single[1][kk] - y_est) < min_b) {
3156 min_b = fabs(y2_single[1][kk] - y_est);
3157 rh_segDC2[3][nDC2_segments] = y2_single[1][kk];
3158 rh_sigm_segDC2[3][nDC2_segments] = sigm_y2_single[1][kk];
3159 foundY = true;
3160 }
3161 } // for kk
3162 }
3163 } //(pair_y2>0)
3164 if (foundX || foundY)
3165 nDC2_segments++;
3166 } // j
3167 } // i
3168 } //(u,v) first
3169
3170 // linear fit for dch1 segs
3171
3172 for (int j = 0; j < nDC1_segments; j++) {
3173 // fit the initial seg and see if the chi2/ndof is too big try to get rid of the hit with the biggest
3174 // deviation from the fit
3175 int worst_hit = -1;
3176 double max_resid = 0;
3177 for (int i = 0; i < 8; i++) {
3178 if (rh_segDC1[i][j] != -999)
3179 size_segDC1[j]++;
3180 }
3181 for (int rej = 0; rej < 2; rej++) { // allow 2 passes max 8->7 & 7->6
3182
3183 float rh_seg[8] = {-999};
3184 float rh_sigm_seg[8] = {-999};
3185 float par_ab[4] = {999, 999, 999, 999};
3186
3187 for (int i = 0; i < 8; i++) {
3188 rh_seg[i] = rh_segDC1[i][j];
3189 rh_sigm_seg[i] = rh_sigm_segDC1[i][j];
3190 }
3191
3192 fit_seg(rh_seg, rh_sigm_seg, par_ab, -1, -1); // usual fit without skipping any plane
3193 for (int i = 0; i < 4; i++) {
3194 par_ab1[i][j] = par_ab[i];
3195 }
3196
3197 chi2_DC1[j] = 0;
3198 float resid = 999;
3199 float dev_coord[8] = {9, 9, 9, 9, 9, 9, 9, 9};
3200
3201 for (int i = 0; i < 8; i++) {
3202 if (i == 0 && rh_segDC1[i][j] != -999) {
3203
3204 resid = rh_segDC1[i][j] - z_loc[i] * par_ab1[0][j] - par_ab1[1][j];
3205
3206 if (fabs(resid) > max_resid) {
3207 worst_hit = i;
3208 max_resid = fabs(resid);
3209 }
3210 chi2_DC1[j] += (resid * resid) / rh_sigm_segDC1[i][j];
3211 }
3212 if (i == 1 && rh_segDC1[i][j] != -999) {
3213 resid = rh_segDC1[i][j] - z_loc[i] * par_ab1[0][j] - par_ab1[1][j];
3214 if (fabs(resid) > max_resid) {
3215 worst_hit = i;
3216 max_resid = fabs(resid);
3217 }
3218 chi2_DC1[j] += (resid * resid) / rh_sigm_segDC1[i][j];
3219 }
3220 if (i == 2 && rh_segDC1[i][j] != -999) {
3221 resid = rh_segDC1[i][j] - z_loc[i] * par_ab1[2][j] - par_ab1[3][j];
3222 if (fabs(resid) > max_resid) {
3223 worst_hit = i;
3224 max_resid = fabs(resid);
3225 }
3226 chi2_DC1[j] += (resid * resid) / rh_sigm_segDC1[i][j];
3227 }
3228 if (i == 3 && rh_segDC1[i][j] != -999) {
3229 resid = rh_segDC1[i][j] - z_loc[i] * par_ab1[2][j] - par_ab1[3][j];
3230
3231 if (fabs(resid) > max_resid) {
3232 worst_hit = i;
3233 max_resid = fabs(resid);
3234 }
3235 chi2_DC1[j] += (resid * resid) / rh_sigm_segDC1[i][j];
3236 }
3237 if (i == 4 && rh_segDC1[i][j] != -999) {
3238 resid = rh_segDC1[i][j] - isqrt_2 * z_loc[i] * (par_ab1[2][j] - par_ab1[0][j])
3239 - isqrt_2 * (par_ab1[3][j] - par_ab1[1][j]);
3240 if (fabs(resid) > max_resid) {
3241 worst_hit = i;
3242 max_resid = fabs(resid);
3243 }
3244 chi2_DC1[j] += (resid * resid) / rh_sigm_segDC1[i][j];
3245 }
3246 if (i == 5 && rh_segDC1[i][j] != -999) {
3247 resid = rh_segDC1[i][j] - isqrt_2 * z_loc[i] * (par_ab1[2][j] - par_ab1[0][j])
3248 - isqrt_2 * (par_ab1[3][j] - par_ab1[1][j]);
3249 if (fabs(resid) > max_resid) {
3250 worst_hit = i;
3251 max_resid = fabs(resid);
3252 }
3253 chi2_DC1[j] += (resid * resid) / rh_sigm_segDC1[i][j];
3254 }
3255 if (i == 6 && rh_segDC1[i][j] != -999) {
3256 resid = rh_segDC1[i][j] - isqrt_2 * z_loc[i] * (par_ab1[2][j] + par_ab1[0][j])
3257 - isqrt_2 * (par_ab1[3][j] + par_ab1[1][j]);
3258 if (fabs(resid) > max_resid) {
3259 worst_hit = i;
3260 max_resid = fabs(resid);
3261 }
3262 chi2_DC1[j] += (resid * resid) / rh_sigm_segDC1[i][j];
3263 }
3264 if (i == 7 && rh_segDC1[i][j] != -999) {
3265
3266 resid = rh_segDC1[i][j] - isqrt_2 * z_loc[i] * (par_ab1[2][j] + par_ab1[0][j])
3267 - isqrt_2 * (par_ab1[3][j] + par_ab1[1][j]);
3268 if (fabs(resid) > max_resid) {
3269 worst_hit = i;
3270 max_resid = fabs(resid);
3271 }
3272 chi2_DC1[j] += (resid * resid) / rh_sigm_segDC1[i][j];
3273 }
3274 } // i
3275
3276 chi2_DC1[j] /= (size_segDC1[j] - 4);
3277
3278 if (chi2_DC1[j] > 30) {
3279 if (size_segDC1[j] == 6) {
3280 chi2_DC1[j] = 999;
3281 break;
3282 } else {
3283 rh_segDC1[worst_hit][j] = -999; // erase worst hit and refit
3284 size_segDC1[j]--;
3285 max_resid = 0;
3286 continue;
3287 }
3288 }
3289 } // rej 0 1 2
3290
3291 // add shifts to slopes and coords
3292 par_ab1[0][j] += x1_slope_sh + x1_slope_sh * par_ab1[0][j] * par_ab1[0][j];
3293 par_ab1[2][j] += y1_slope_sh + y1_slope_sh * par_ab1[2][j] * par_ab1[2][j];
3294 par_ab1[1][j] += x1_sh;
3295 par_ab1[3][j] += y1_sh;
3296
3297 xDC1_glob[j] = par_ab1[0][j] * (99.5) + par_ab1[1][j];
3298 yDC1_glob[j] = par_ab1[2][j] * (99.5) + par_ab1[3][j];
3299 if (size_segDC1[j] > 6)
3300 has7DC1 = true;
3301 } // for DC1 segs
3302
3303 // linear fit for dch1 segs
3304 for (int j = 0; j < nDC2_segments; j++) {
3305 int worst_hit = -1;
3306 double max_resid = 0;
3307 for (int i = 0; i < 8; i++) {
3308 if (rh_segDC2[i][j] != -999)
3309 size_segDC2[j]++;
3310 }
3311 for (int rej = 0; rej < 2; rej++) { // allow 2 passes max 8->7 & 7->6
3312
3313 // end linear fit
3314 float rh_seg[8] = {-999};
3315 float rh_sigm_seg[8] = {-999};
3316 float par_ab[4] = {999, 999, 999, 999};
3317
3318 for (int i = 0; i < 8; i++) {
3319 rh_seg[i] = rh_segDC2[i][j];
3320 rh_sigm_seg[i] = rh_sigm_segDC2[i][j];
3321 }
3322
3323 fit_seg(rh_seg, rh_sigm_seg, par_ab, -1, -1); // usual fit without skipping any plane
3324 for (int i = 0; i < 4; i++) {
3325 par_ab2[i][j] = par_ab[i];
3326 }
3327
3328 chi2_DC2[j] = 0;
3329
3330 float resid = 999;
3331
3332 for (int i = 0; i < 8; i++) {
3333 if (i == 0 && rh_segDC2[i][j] != -999) {
3334 resid = rh_segDC2[i][j] - z_loc[i] * par_ab2[0][j] - par_ab2[1][j];
3335 chi2_DC2[j] += (resid * resid) / rh_sigm_segDC2[i][j];
3336 if (fabs(resid) > max_resid) {
3337 worst_hit = i;
3338 max_resid = fabs(resid);
3339 }
3340 }
3341 if (i == 1 && rh_segDC2[i][j] != -999) {
3342 resid = rh_segDC2[i][j] - z_loc[i] * par_ab2[0][j] - par_ab2[1][j];
3343
3344 chi2_DC2[j] += (resid * resid) / rh_sigm_segDC2[i][j];
3345
3346 if (fabs(resid) > max_resid) {
3347 worst_hit = i;
3348 max_resid = fabs(resid);
3349 }
3350 }
3351 if (i == 2 && rh_segDC2[i][j] != -999) {
3352 resid = rh_segDC2[i][j] - z_loc[i] * par_ab2[2][j] - par_ab2[3][j];
3353
3354 chi2_DC2[j] += (resid * resid) / rh_sigm_segDC2[i][j];
3355
3356 if (fabs(resid) > max_resid) {
3357 worst_hit = i;
3358 max_resid = fabs(resid);
3359 }
3360 }
3361 if (i == 3 && rh_segDC2[i][j] != -999) {
3362 resid = rh_segDC2[i][j] - z_loc[i] * par_ab2[2][j] - par_ab2[3][j];
3363
3364 chi2_DC2[j] += (resid * resid) / rh_sigm_segDC2[i][j];
3365
3366 if (fabs(resid) > max_resid) {
3367 worst_hit = i;
3368 max_resid = fabs(resid);
3369 }
3370 }
3371 if (i == 4 && rh_segDC2[i][j] != -999) {
3372 resid = rh_segDC2[i][j] - isqrt_2 * z_loc[i] * (par_ab2[2][j] - par_ab2[0][j])
3373 - isqrt_2 * (par_ab2[3][j] - par_ab2[1][j]);
3374 chi2_DC2[j] += (resid * resid) / rh_sigm_segDC2[i][j];
3375
3376 if (fabs(resid) > max_resid) {
3377 worst_hit = i;
3378 max_resid = fabs(resid);
3379 }
3380 }
3381 if (i == 5 && rh_segDC2[i][j] != -999) {
3382 resid = rh_segDC2[i][j] - isqrt_2 * z_loc[i] * (par_ab2[2][j] - par_ab2[0][j])
3383 - isqrt_2 * (par_ab2[3][j] - par_ab2[1][j]);
3384 chi2_DC2[j] += (resid * resid) / rh_sigm_segDC2[i][j];
3385 if (fabs(resid) > max_resid) {
3386 worst_hit = i;
3387 max_resid = fabs(resid);
3388 }
3389 }
3390 if (i == 6 && rh_segDC2[i][j] != -999) {
3391 resid = rh_segDC2[i][j] - isqrt_2 * z_loc[i] * (par_ab2[2][j] + par_ab2[0][j])
3392 - isqrt_2 * (par_ab2[3][j] + par_ab2[1][j]);
3393 chi2_DC2[j] += (resid * resid) / rh_sigm_segDC2[i][j];
3394
3395 if (fabs(resid) > max_resid) {
3396 worst_hit = i;
3397 max_resid = fabs(resid);
3398 }
3399 }
3400 if (i == 7 && rh_segDC2[i][j] != -999) {
3401 resid = rh_segDC2[i][j] - isqrt_2 * z_loc[i] * (par_ab2[2][j] + par_ab2[0][j])
3402 - isqrt_2 * (par_ab2[3][j] + par_ab2[1][j]);
3403 chi2_DC2[j] += (resid * resid) / rh_sigm_segDC2[i][j];
3404
3405 if (fabs(resid) > max_resid) {
3406 worst_hit = i;
3407 max_resid = fabs(resid);
3408 }
3409 }
3410 } // i
3411 chi2_DC2[j] /= (size_segDC2[j] - 4);
3412 // if chi2 is big and seg_size = min erase this seg
3413 if (chi2_DC2[j] > 30) {
3414 if (size_segDC2[j] == 6) {
3415 chi2_DC2[j] = 999;
3416 break;
3417 } else {
3418 rh_segDC2[worst_hit][j] = -999; // erase worst hit and refit
3419 size_segDC2[j]--;
3420 max_resid = 0;
3421 continue;
3422 }
3423 }
3424 } // rej 0 1 2
3425
3426 // add shifts to slopes and coords
3427 par_ab2[0][j] += x2_slope_sh + x2_slope_sh * par_ab2[0][j] * par_ab2[0][j];
3428 par_ab2[2][j] += y2_slope_sh + y2_slope_sh * par_ab2[2][j] * par_ab2[2][j];
3429 par_ab2[1][j] += x2_sh;
3430 par_ab2[3][j] += y2_sh;
3431
3432 xDC2_glob[j] = par_ab2[0][j] * (-99.5) + par_ab2[1][j];
3433 yDC2_glob[j] = par_ab2[2][j] * (-99.5) + par_ab2[3][j];
3434 //
3435 if (size_segDC2[j] > 6)
3436 has7DC2 = true;
3437 } // for DC2 segs
3438
3439 // count the number of rh per segment
3440 double x_mid[50]; // x glob of matched segment in the z situated between the two DCH chambers
3441 double y_mid[50]; // y glob of matched segment in the z situated between the two DCH chambers
3442 double a_X[50]; // x slope of the matched segment
3443 double a_Y[50]; // y slope of the matched segment
3444 double imp[50]; // reconstructed particle trajectory momentum
3445 double leng[50]; // the distance from z = 0 to the global point of the matched segment
3446 int seg_it = -1;
3447 int dc1_best[50];
3448 int dc2_best[50];
3449 float Chi2_match[50]; // chi2 of the matched seg
3450 for (int seg = 0; seg < 50; seg++) {
3451 x_mid[seg] = -999;
3452 y_mid[seg] = -999;
3453 a_X[seg] = -999;
3454 a_Y[seg] = -999;
3455 imp[seg] = -999;
3456 leng[seg] = -999;
3457 Chi2_match[seg] = -999;
3458 dc1_best[seg] = 0;
3459 dc2_best[seg] = 0;
3460 }
3461
3462 float xMean = 999;
3463 float yMean = 999;
3464
3465 // leave only longest and best chi2 segments
3466 // dc1
3467 for (int max_size = 8; max_size > 5; max_size--) {
3468 // find longest and best chi2 seg
3469 for (int sizeit1 = 0; sizeit1 < nDC1_segments; sizeit1++) {
3470 if (size_segDC1[sizeit1] != max_size)
3471 continue;
3472 for (int sizeit1_1 = 0; sizeit1_1 < nDC1_segments; sizeit1_1++) {
3473 if (sizeit1_1 == sizeit1)
3474 continue;
3475 for (int hit = 0; hit < 4; hit++) {
3476 if (rh_segDC1[2 * hit][sizeit1] == rh_segDC1[2 * hit][sizeit1_1]
3477 && rh_segDC1[2 * hit + 1][sizeit1] == rh_segDC1[2 * hit + 1][sizeit1_1]
3478 && (chi2_DC1[sizeit1] <= chi2_DC1[sizeit1_1]
3479 || size_segDC1[sizeit1] > size_segDC1[sizeit1_1]))
3480 {
3481 chi2_DC1[sizeit1_1] = 999; // mark seg as bad
3482 break;
3483 }
3484 } // hit
3485 }
3486 }
3487 } // max_size
3488
3489 // dc2
3490 for (int max_size = 8; max_size > 5; max_size--) {
3491 // find longest and best chi2 seg
3492 for (int sizeit2 = 0; sizeit2 < nDC2_segments; sizeit2++) {
3493 if (size_segDC2[sizeit2] != max_size)
3494 continue;
3495 for (int sizeit2_2 = 0; sizeit2_2 < nDC2_segments; sizeit2_2++) {
3496 if (sizeit2_2 == sizeit2)
3497 continue;
3498 for (int hit = 0; hit < 4; hit++) {
3499 if (rh_segDC2[2 * hit][sizeit2] == rh_segDC2[2 * hit][sizeit2_2]
3500 && rh_segDC2[2 * hit + 1][sizeit2] == rh_segDC2[2 * hit + 1][sizeit2_2]
3501 && (chi2_DC2[sizeit2] <= chi2_DC2[sizeit2_2]
3502 || size_segDC2[sizeit2] > size_segDC2[sizeit2_2]))
3503 {
3504 chi2_DC2[sizeit2_2] = 999; // mark seg as bad
3505 break;
3506 }
3507 } // hit
3508 }
3509 }
3510 } // max_size
3511 // fill local segments z,x,y global coords; x-slope; y-slope; Chi2; to be continued...
3512 for (int sizeit1 = 0; sizeit1 < nDC1_segments; sizeit1++) {
3513 if (chi2_DC1[sizeit1] > 50)
3514 continue;
3515 BmnTrack TracksDch1;
3516 FairTrackParam ParamsTrackDch1;
3517 Float_t z0 = Z_dch1;
3518 Float_t x0 = par_ab1[1][sizeit1];
3519 Float_t y0 = par_ab1[3][sizeit1];
3520 ParamsTrackDch1.SetPosition(TVector3(x0, y0, z0));
3521 ParamsTrackDch1.SetTx(par_ab1[0][sizeit1]);
3522 ParamsTrackDch1.SetTy(par_ab1[2][sizeit1]);
3523 // TracksDch1.SetParamFirst(ParamsTrackDch1);
3524 // TracksDch1.SetChi2(chi2_DC1[sizeit1]);
3525 // TracksDch1.SetNDF(Nhits_Dch1[ind_best_Dch1[iBest]]); //WARNING!!! now it is number of hits in matched
3526 // track!
3527 BmnTrack* tr_dch1 = new ((*dchTracks)[dchTracks->GetEntriesFast()]) BmnTrack();
3528 tr_dch1->SetChi2(chi2_DC1[sizeit1]);
3529 tr_dch1->SetNHits(size_segDC1[sizeit1]);
3530 tr_dch1->SetParamFirst(ParamsTrackDch1);
3531 }
3532
3533 for (int sizeit2 = 0; sizeit2 < nDC2_segments; sizeit2++) {
3534 if (chi2_DC2[sizeit2] > 50)
3535 continue;
3536 BmnTrack TracksDch2;
3537 FairTrackParam ParamsTrackDch2;
3538 Float_t z0 = Z_dch2;
3539 Float_t x0 = par_ab2[1][sizeit2];
3540 Float_t y0 = par_ab2[3][sizeit2];
3541 ParamsTrackDch2.SetPosition(TVector3(x0, y0, z0));
3542 ParamsTrackDch2.SetTx(par_ab2[0][sizeit2]);
3543 ParamsTrackDch2.SetTy(par_ab2[2][sizeit2]);
3544 TracksDch2.SetParamFirst(ParamsTrackDch2);
3545 // TracksDch2.SetChi2(chi2_DC2[sizeit2]);
3546 // TracksDch2.SetNDF(Nhits_Dch2[ind_best_Dch2[iBest]]); //WARNING!!! now it is number of hits in matched
3547 // track!
3548
3549 BmnTrack* tr_dch2 = new ((*dchTracks)[dchTracks->GetEntriesFast()]) BmnTrack();
3550 tr_dch2->SetChi2(chi2_DC2[sizeit2]);
3551 tr_dch2->SetNHits(size_segDC2[sizeit2]);
3552 tr_dch2->SetParamFirst(ParamsTrackDch2);
3553 }
3554
3555 // try to match the reconstructed segments from the two chambers
3556 if (has7DC1) {
3557
3558 int match_dc2_seg = -1;
3559 float ax = -999.0;
3560 float ay = -999.0;
3561 float ax1 = -999.0;
3562 float ax2 = -999.0;
3563 float ay1 = -999.0;
3564 float ay2 = -999.0;
3565 float bx1, bx2, by1, by2;
3566 bx1 = bx2 = by1 = by2 = -999.0;
3567
3568 for (int segdc1Nr = 0; segdc1Nr < nDC1_segments; segdc1Nr++) {
3569 if (chi2_DC1[segdc1Nr] > 50 || size_segDC1[segdc1Nr] < 7)
3570 continue; // skip rejected segs with chi2 = 999
3571 float min_distX = 20;
3572 float min_distY = 15;
3573 float min_distSQ = 225;
3574
3575 float dx = -999;
3576 float dy = -999;
3577 float daX = -999;
3578 float daY = -999;
3579 float chi2_match = 0;
3580
3581 for (int segdc2Nr = 0; segdc2Nr < nDC2_segments; segdc2Nr++) {
3582 if (chi2_DC2[segdc2Nr] > 50)
3583 continue; // skip rejected segs with chi2 = 999
3584
3585 float distX = fabs(xDC1_glob[segdc1Nr] - xDC2_glob[segdc2Nr]);
3586 float distY = fabs(yDC1_glob[segdc1Nr] - yDC2_glob[segdc2Nr]);
3587 if (distX < min_distX && distY < min_distY) {
3588 float distSQ = distX * distX + distY * distY;
3589
3590 chi2_match = (fabs(xDC1_glob[segdc1Nr] - xDC2_glob[segdc2Nr])
3591 * fabs(xDC1_glob[segdc1Nr] - xDC2_glob[segdc2Nr]) / 49)
3592 + (distY * distY / 43.56)
3593 + ((par_ab2[0][segdc2Nr] - par_ab1[0][segdc1Nr] * par_ab2[0][segdc2Nr]
3594 - par_ab1[0][segdc1Nr])
3595 / 0.0144)
3596 + (par_ab2[2][segdc2Nr] - par_ab1[2][segdc1Nr] * par_ab2[2][segdc2Nr]
3597 - par_ab1[2][segdc1Nr])
3598 / 0.0225;
3599 if (distSQ < min_distSQ) {
3600
3601 dx = xDC2_glob[segdc2Nr] - xDC1_glob[segdc1Nr];
3602 dy = yDC2_glob[segdc2Nr] - yDC1_glob[segdc1Nr];
3603 xMean = 0.5 * (xDC1_glob[segdc1Nr] + xDC2_glob[segdc2Nr]);
3604 yMean = 0.5 * (yDC1_glob[segdc1Nr] + yDC2_glob[segdc2Nr]);
3605 ax = (par_ab2[1][segdc2Nr] - par_ab1[1][segdc1Nr]) / 199;
3606 ay = (par_ab2[3][segdc2Nr] - par_ab1[3][segdc1Nr]) / 199;
3607
3608 min_distSQ = distSQ;
3609 match_dc2_seg = segdc2Nr; // remember matched dc2 seg
3610 daX = par_ab2[0][segdc2Nr] - par_ab1[0][segdc1Nr];
3611 daY = par_ab2[2][segdc2Nr] - par_ab1[2][segdc1Nr];
3612 ax1 = par_ab1[0][segdc1Nr];
3613 ax2 = par_ab2[0][segdc2Nr];
3614 ay1 = par_ab1[2][segdc1Nr];
3615 ay2 = par_ab2[2][segdc2Nr];
3616 bx2 = par_ab2[1][segdc2Nr];
3617 by2 = par_ab2[3][segdc2Nr];
3618 } // min_dist
3619 } // distX<min_distX && distY<min_distY
3620
3621 } // segdc2Nr
3622
3623 if (min_distSQ == 225)
3624 continue;
3625
3626 //-n if(fabs(dx) < 7 && fabs(dy)<8.2){
3627 if (fabs(dx) < 10 && fabs(dy) < 10) {
3628 if (fabs(daX) < 0.2 && fabs(daY) < 0.2) {
3629 seg_it++;
3630 leng[seg_it] = sqrt(628.65 * 628.65 + (xMean * xMean));
3631 x_mid[seg_it] = xMean;
3632 y_mid[seg_it] = yMean;
3633 // a_X[seg_it] = atan(ax);
3634 // a_Y[seg_it] = atan(ay);
3635 a_X[seg_it] = ax;
3636 a_Y[seg_it] = ay;
3637 imp[seg_it] = -0.4332 / (ax + 0.006774);
3638 Chi2_match[seg_it] = chi2_match;
3639 chi2_DC1[segdc1Nr] = 999; // mark dch1 seg as used for future iterations
3640 chi2_DC2[match_dc2_seg] = 999; // mark dch2 seg as used for future iterations
3641
3642 } // dax day
3643 } // dx<0.1 dy
3644 } // segdc1Nr
3645 } // if(has7DC1) 8p+6p
3646
3647 if (has7DC2) {
3648 int match_dc1_seg = -1;
3649
3650 float ax = -999.0;
3651 float ay = -999.0;
3652 float ax1 = -999.0;
3653 float ax2 = -999.0;
3654 float ay1 = -999.0;
3655 float ay2 = -999.0;
3656
3657 float bx1, bx2, by1, by2;
3658 bx1 = bx2 = by1 = by2 = -999.0;
3659
3660 for (int segdc2Nr = 0; segdc2Nr < nDC2_segments; segdc2Nr++) {
3661
3662 float min_distX = 20;
3663 float min_distY = 15;
3664 float min_distSQ = 225;
3665
3666 float dx = -999;
3667 float dy = -999;
3668 float daX = -999;
3669 float daY = -999;
3670
3671 if (chi2_DC2[segdc2Nr] > 50 || size_segDC2[segdc2Nr] < 7)
3672 continue;
3673
3674 float chi2_match = 0;
3675
3676 for (int segdc1Nr = 0; segdc1Nr < nDC1_segments; segdc1Nr++) {
3677 if (chi2_DC1[segdc1Nr] > 50)
3678 continue;
3679
3680 float distX = fabs(xDC1_glob[segdc1Nr] - xDC2_glob[segdc2Nr]);
3681 float distY = fabs(yDC1_glob[segdc1Nr] - yDC2_glob[segdc2Nr]);
3682
3683 if (distX < min_distX && distY < min_distY) {
3684 float distSQ = distX * distX + distY * distY;
3685
3686 chi2_match = (fabs(xDC1_glob[segdc1Nr] - xDC2_glob[segdc2Nr])
3687 * fabs(xDC1_glob[segdc1Nr] - xDC2_glob[segdc2Nr]) / 49)
3688 + (distY * distY / 43.56)
3689 + ((par_ab2[0][segdc2Nr] - par_ab1[0][segdc1Nr] * par_ab2[0][segdc2Nr]
3690 - par_ab1[0][segdc1Nr])
3691 / 0.0144)
3692 + (par_ab2[2][segdc2Nr] - par_ab1[2][segdc1Nr] * par_ab2[2][segdc2Nr]
3693 - par_ab1[2][segdc1Nr])
3694 / 0.0225;
3695
3696 if (distSQ < min_distSQ) {
3697
3698 dx = xDC2_glob[segdc2Nr] - xDC1_glob[segdc1Nr];
3699 dy = yDC2_glob[segdc2Nr] - yDC1_glob[segdc1Nr];
3700 xMean = 0.5 * (xDC1_glob[segdc1Nr] + xDC2_glob[segdc2Nr]);
3701 yMean = 0.5 * (yDC1_glob[segdc1Nr] + yDC2_glob[segdc2Nr]);
3702 min_distSQ = distSQ;
3703 ax = (par_ab2[1][segdc2Nr] - par_ab1[1][segdc1Nr]) / 199;
3704 ay = (par_ab2[3][segdc2Nr] - par_ab1[3][segdc1Nr]) / 199;
3705 match_dc1_seg = segdc1Nr;
3706 daX = par_ab2[0][segdc2Nr] - par_ab1[0][segdc1Nr];
3707 daY = par_ab2[2][segdc2Nr] - par_ab1[2][segdc1Nr];
3708 } // min_dist
3709 }
3710 } // segdc1Nr
3711
3712 if (min_distSQ >= 225)
3713 continue;
3714
3715 if (fabs(dx) < 10 && fabs(dy) < 10) {
3716 if (fabs(daX) < 0.2 && fabs(daY) < 0.2) {
3717
3718 seg_it++;
3719 leng[seg_it] = sqrt(628.65 * 628.65 + (xMean * xMean));
3720 x_mid[seg_it] = xMean;
3721 y_mid[seg_it] = yMean;
3722 a_X[seg_it] = atan(ax);
3723 a_Y[seg_it] = atan(ay);
3724 imp[seg_it] = -0.4332 / (ax + 0.006774);
3725 Chi2_match[seg_it] = chi2_match;
3726 chi2_DC2[segdc2Nr] = 999; // mark dch2 seg as used for future iterations
3727 chi2_DC1[match_dc1_seg] = 999; // mark dch1 seg as used for future iterations
3728
3729 } // dax day
3730 } // dx dy
3731 } // segdc2Nr
3732 } // 6p+8p has7DC2
3733
3734 // do 6+6p glob segs
3735
3736 int match_dc1_seg = -1;
3737
3738 float ax = -999.0;
3739 float ay = -999.0;
3740 float ax1 = -999.0;
3741 float ax2 = -999.0;
3742 float ay1 = -999.0;
3743 float ay2 = -999.0;
3744
3745 float bx1, bx2, by1, by2;
3746 bx1 = bx2 = by1 = by2 = -999.0;
3747
3748 for (int segdc2Nr = 0; segdc2Nr < nDC2_segments; segdc2Nr++) {
3749 float min_distX = 20;
3750 float min_distY = 15;
3751 float min_distSQ = 225;
3752
3753 float dx = -999;
3754 float dy = -999;
3755 float daX = -999;
3756 float daY = -999;
3757
3758 if (chi2_DC2[segdc2Nr] > 50)
3759 continue;
3760
3761 float chi2_match = 0;
3762 for (int segdc1Nr = 0; segdc1Nr < nDC1_segments; segdc1Nr++) {
3763 if (chi2_DC1[segdc1Nr] > 50)
3764 continue;
3765
3766 float distX = fabs(xDC1_glob[segdc1Nr] - xDC2_glob[segdc2Nr]);
3767 float distY = fabs(yDC1_glob[segdc1Nr] - yDC2_glob[segdc2Nr]);
3768
3769 if (distX < min_distX && distY < min_distY) {
3770 float distSQ = distX * distX + distY * distY;
3771
3772 chi2_match =
3773 (fabs(xDC1_glob[segdc1Nr] - xDC2_glob[segdc2Nr])
3774 * fabs(xDC1_glob[segdc1Nr] - xDC2_glob[segdc2Nr]) / 49)
3775 + (distY * distY / 43.56)
3776 + ((par_ab2[0][segdc2Nr] - par_ab1[0][segdc1Nr] * par_ab2[0][segdc2Nr] - par_ab1[0][segdc1Nr])
3777 / 0.0144)
3778 + (par_ab2[2][segdc2Nr] - par_ab1[2][segdc1Nr] * par_ab2[2][segdc2Nr] - par_ab1[2][segdc1Nr])
3779 / 0.0225;
3780 if (distSQ < min_distSQ) {
3781
3782 dx = xDC2_glob[segdc2Nr] - xDC1_glob[segdc1Nr];
3783 dy = yDC2_glob[segdc2Nr] - yDC1_glob[segdc1Nr];
3784 xMean = 0.5 * (xDC1_glob[segdc1Nr] + xDC2_glob[segdc2Nr]);
3785 yMean = 0.5 * (yDC1_glob[segdc1Nr] + yDC2_glob[segdc2Nr]);
3786 min_distSQ = distSQ;
3787 ax = (par_ab2[1][segdc2Nr] - par_ab1[1][segdc1Nr]) / 199;
3788 ay = (par_ab2[3][segdc2Nr] - par_ab1[3][segdc1Nr]) / 199;
3789 match_dc1_seg = segdc1Nr;
3790 daX = par_ab2[0][segdc2Nr] - par_ab1[0][segdc1Nr];
3791 daY = par_ab2[2][segdc2Nr] - par_ab1[2][segdc1Nr];
3792
3793 } // min_dist
3794 }
3795 } // segdc1Nr
3796
3797 if (min_distSQ >= 225)
3798 continue;
3799
3800 if (fabs(dx) < 10 && fabs(dy) < 10) {
3801 if (fabs(daX) < 0.2 && fabs(daY) < 0.2) {
3802
3803 seg_it++;
3804 leng[seg_it] = sqrt(628.65 * 628.65 + (xMean * xMean));
3805 x_mid[seg_it] = xMean;
3806 y_mid[seg_it] = yMean;
3807 a_X[seg_it] = atan(ax);
3808 a_Y[seg_it] = atan(ay);
3809 imp[seg_it] = -0.4332 / (ax + 0.006774);
3810 Chi2_match[seg_it] = chi2_match;
3811 chi2_DC2[segdc2Nr] = 999; // mark dc2 seg as used for future iterations
3812 chi2_DC1[match_dc1_seg] = 999; // mark dch1 seg as used for future iterations
3813 } // dax day
3814 } // dx dy
3815 } // segdc2Nr
3816
3817 if (seg_it > -1) {
3818 for (int seg = 0; seg < seg_it + 1; seg++) {
3819 // CbmTrack TracksDchMatch;
3820 FairTrackParam ParamsTrackDchMatch;
3821 Float_t z0 = Z_dch_mid;
3822 Float_t x0 = x_mid[seg];
3823 Float_t y0 = y_mid[seg];
3824 ParamsTrackDchMatch.SetPosition(TVector3(x0, y0, z0));
3825 ParamsTrackDchMatch.SetTx(a_X[seg]);
3826 ParamsTrackDchMatch.SetTy(a_Y[seg]);
3827 // TracksDchMatch.SetParamFirst(&ParamsTrackDchMatch);
3828 // TracksDchMatch.SetChiSq(Chi2_match[seg]);
3829 // TracksDch2.SetNDF(Nhits_Dch2[ind_best_Dch2[iBest]]); //WARNING!!! now it is number of hits in matched
3830 // track!
3831 BmnTrack* newTrack = new ((*dchMatchedTracks)[dchMatchedTracks->GetEntriesFast()]) BmnTrack();
3832 newTrack->SetChi2(Chi2_match[seg]);
3833 newTrack->SetParamFirst(ParamsTrackDchMatch);
3834 }
3835 } // seg_it>0
3836
3837 tReco->Fill();
3838
3839 } // event loop
3840
3841 fReco->cd();
3842 tReco->Write();
3843 fReco->Close();
3844}
3845
3846// LocalWords: GetWireNumber X
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
Int_t
int Int_t
Float_t
float Float_t
Short_t
short Short_t
Double_t
double Double_t
gROOT
#define gROOT
Form
char * Form(const char *fmt,...)
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
BmnTrack::SetChi2
void SetChi2(Double_t chi2)
Definition BmnTrack.h:76
BmnTrack::SetNHits
void SetNHits(Int_t n)
Definition BmnTrack.h:80
BmnTrack::SetParamFirst
void SetParamFirst(FairTrackParam &par)
Modifiers.
Definition BmnTrack.h:70
TChain::SetBranchAddress
Int_t SetBranchAddress(const char *bname, T **add, TBranch **ptr=nullptr)
TChain::Add
virtual Int_t Add(const char *name, Long64_t nentries=TTree::kMaxEntries)
TChain::GetEntries
Long64_t GetEntries() const override
TChain::GetEntry
Int_t GetEntry(Long64_t entry=0, Int_t getall=0) override
TClonesArray::Clear
void Clear(Option_t *option="") override
TDirectoryFile::cd
Bool_t cd() override
TFile::Close
void Close(Option_t *option="") override
TObjArray::GetEntriesFast
Int_t GetEntriesFast() const
TObjArray::At
TObject * At(Int_t idx) const override
TString::Format
static TString Format(const char *fmt,...)
TTree::Fill
virtual Int_t Fill()
TTree::Write
Int_t Write(const char *name=nullptr, Int_t option=0, Int_t bufsize=0) const override
TTree::Branch
virtual Int_t Branch(const char *folder, Int_t bufsize=32000, Int_t splitlevel=99)
sqrt
RVec< PromoteType< T > > sqrt(const RVec< T > &v)
atan
RVec< PromoteType< T > > atan(const RVec< T > &v)
v1_slope_sh
const Float_t v1_slope_sh
Definition recoRunDCH_VP.C:123
v2_slope_sh
const Float_t v2_slope_sh
Definition recoRunDCH_VP.C:124
y2_slope_sh
const Float_t y2_slope_sh
Definition recoRunDCH_VP.C:119
sqrt_2
double sqrt_2
Definition recoRunDCH_VP.C:80
y1_slope_sh
const Float_t y1_slope_sh
Definition recoRunDCH_VP.C:117
Z_dch_mid
const Float_t Z_dch_mid
Definition recoRunDCH_VP.C:77
z_glob
const Float_t z_glob[16]
Definition recoRunDCH_VP.C:73
y1_sh
const Float_t y1_sh
Definition recoRunDCH_VP.C:108
u2_sh
const Float_t u2_sh
Definition recoRunDCH_VP.C:112
u1_sh
const Float_t u1_sh
Definition recoRunDCH_VP.C:111
u2_slope_sh
const Float_t u2_slope_sh
Definition recoRunDCH_VP.C:122
fit_seg
void fit_seg(float *rh_seg, float *rh_sigm_seg, float *par_ab, int skip_first, int skip_second)
Definition recoRunDCH_VP.C:126
recoRunDCH_VP
void recoRunDCH_VP(Int_t runId=810)
Definition recoRunDCH_VP.C:312
u1_slope_sh
const Float_t u1_slope_sh
Definition recoRunDCH_VP.C:121
y2_sh
const Float_t y2_sh
Definition recoRunDCH_VP.C:109
intBL
float intBL
Definition recoRunDCH_VP.C:83
v1_sh
const Float_t v1_sh
Definition recoRunDCH_VP.C:113
Z_dch2
const Float_t Z_dch2
Definition recoRunDCH_VP.C:76
Z_dch1
const Float_t Z_dch1
Definition recoRunDCH_VP.C:75
x2_sh
const Float_t x2_sh
Definition recoRunDCH_VP.C:107
x2_slope_sh
const Float_t x2_slope_sh
Definition recoRunDCH_VP.C:118
isqrt_2
double isqrt_2
Definition recoRunDCH_VP.C:81
x1_sh
const Float_t x1_sh
Definition recoRunDCH_VP.C:106
x1_slope_sh
const Float_t x1_slope_sh
Definition recoRunDCH_VP.C:116
v2_sh
const Float_t v2_sh
Definition recoRunDCH_VP.C:114
z_loc
const Float_t z_loc[8]
Definition recoRunDCH_VP.C:72
std
Definition json.hpp:4446