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