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