BmnRoot
Loading...
Searching...
No Matches
BmnUpstreamTracking.cxx
Go to the documentation of this file.
1// Author: Vasilisa Lenivenko (VBLHEP) <vasilisa@jinr.ru> 2019-07-18
2
4// //
5// BmnUpstreamTracking //
6// //
7// Implementation of an algorithm developed by //
8// V.Lenivenko and V.Paltchik //
9// to the BmnRoot software //
10// //
11// The algorithm serves for searching for tracks //
12// in the Silicon detectors & MWPCs //
13// of the configuration SRC of BM@N experiment //
14// //
16
17#include "BmnUpstreamTracking.h"
18
19#include "FairLogger.h"
20#include "FairRootManager.h"
21
22#include <TMath.h>
23#include "TFile.h"
24
25#include <vector>
26
27static Double_t workTime = 0.0;
28TList fhList;
29
30//----------------------------------------------------------------------
31void BmnUpstreamTracking::Exec(Option_t* opt) {
32 if (!IsActive()) return;
33
34 clock_t tStart = clock();
35 PrepareArraysToProcessEvent();
36 if (fDebug) cout << "\n======================== Upstream track finder exec started ===================\n" << endl;
37 if (fDebug) printf("Event number: %d\n", fEventNo);
38 fEventNo++;
39
40 if (fDebug) cout<<"--------------Read-------------"<<endl;
41 ReadSiliconTracks(par_ab_SiTr,par_SiTr_z,NSiTracks, vec_points);
42 ReadSiliconHits(SiXYhits,NSiXYhits);
43
44 if (fDebug) cout<<"--------------Reco: MWPCSegments-----"<<endl;
45 ReadMWPCSegments(par_ab_Ch,par_Seg_z,Nseg_Ch, vec_points);
46 if (fDebug) cout<<"--------------Reco: MWPCTracks"<<endl;
47 ReadMWPCTracks(par_ab_tr,par_ab_trz,NPCTracks);
48
49 //-----------------Matching-------------------------------------------
50 if (fDebug) cout<<"--------------Matching Si-tracks & PC-tracks-----"<<endl;
51 PCTracksSiTracksMatching(par_ab_SiTr,par_SiTr_z,NSiTracks,par_ab_tr,par_ab_trz,NPCTracks,vtmpSeg,OutVector);//Step1
52 if (fDebug) cout<<"--------------Matching Si-tracks & PC-segments---"<<endl;
53 PCSegmentsSiTracksMatching(par_ab_SiTr,par_SiTr_z,NSiTracks,par_ab_Ch,par_Seg_z,Nseg_Ch,vtmpSeg,OutVector); //Step2
54 if (fDebug) cout<<"--------------Adding hits-------------------------"<<endl;
55 GetAddSiHits(OutVector,SiXYhits,NSiXYhits);
56 GetHitsToPoints(OutVector,par_ab_tr,NPCTracks,par_ab_Ch,Nseg_Ch,Points,NumPoints);
57
58 if (fDebug) cout<<"--------------Fit--------------------------------"<<endl;
59 CalculationOfChi2(Points, NumPoints, vecUpTracks);
60
61 if (fDebug) PrintAllTracks(vecUpTracks);
62 if (!expData) MCefficiencyCalculation(vec_points, vecUpTracks);
63
64 if (fDebug) cout<<"----TrackRecording--------------------------------"<<endl;
65 TrackRecording(vecUpTracks);
66
67 if (fDebug) cout << endl;
68 if (fDebug) cout << "======================== Upstream track finder exec finished ===================" << endl;
69 clock_t tFinish = clock();
70 workTime += ((Double_t) (tFinish - tStart)) / CLOCKS_PER_SEC;
71}
72//----------------------------------------------------------------------
73
74void BmnUpstreamTracking::MCefficiencyCalculation(vector<MC_points> &points, vector<UpTracks> &vecUp){
75
76 if (fDebug) cout<<" ---UpMC tracks association--"<<endl;
77 // ax, bx, ay, by
78 Double_t delta2[4] = {-999.,-999.,-999.,-999.};
79 Double_t sig[4] = {0.03, 4., 0.03, 4.};
80
81 Double_t dmatch = 0.;
82 Double_t dmc_match[kMaxMC];
83 Int_t mc_tr_assoc[kMaxMC];
84
85 for (Int_t j = 0; j < kMaxMC; j++) {
86 dmc_match[j] = 1000.;
87 mc_tr_assoc[j] = -1;
88 }
89
90 Double_t dax = -999.;
91 Double_t day = -999.;
92 Double_t dx = -999.;
93 Double_t dy = -999.;
94
95 if (fDebug) cout<<" MC "<<points.size()<<" UpTracks "<<vecUp.size()<<endl;
96
97 Int_t Ngood_mc_tracks = 0;
98 Int_t Ngood_reco_tracks = 0;
99 for (size_t itr = 0; itr < points.size(); itr++) {//mc_tr
100
101 //---MC Eff ---
102 //---Den
103 if (fDebug ) cout<<" Np "<<points.at(itr).Np<<" pdg "<<points.at(itr).Pdg<<endl;
104 if (fDebug && points.at(itr).Np >=7 //>=10
105 && points.at(itr).np_3si == 0
106 && (points.at(itr).Np_ch2 > 3 || points.at(itr).Np_ch3 > 3) ){
107 // if ( points.at(itr).xWas3 && points.at(itr).uWas3 && points.at(itr).vWas3 &&
108 // points.at(itr).xWas2 && points.at(itr).uWas2 && points.at(itr).vWas2 )
109 hDen_mcuptr->Fill(1);
110 hDen_mcuptr->Fill(0);
111 hY_vs_Xmctrue->Fill(points.at(itr).param[1],points.at(itr).param[3]);
112
113 Ngood_mc_tracks++;
114 cout<<" itr "<<itr<<" Ngood_mc_tracks "<<Ngood_mc_tracks<<endl;
115 cout<<"UpDen"<<endl;
116 }
117
118 for(size_t iup= 0; iup < vecUp.size(); iup ++) {//reco
119 delta2[0] = points.at(itr).param[0] - vecUp.at(iup).param[0];
120 delta2[1] = points.at(itr).param[1] - vecUp.at(iup).param[1];
121 delta2[2] = points.at(itr).param[2] - vecUp.at(iup).param[2];
122 delta2[3] = points.at(itr).param[3] - vecUp.at(iup).param[3];
123
124 dmatch = 0.;
125 dmatch = (delta2[0]*delta2[0])/(sig[0]*sig[0])+
126 (delta2[1]*delta2[1])/(sig[1]*sig[1])+
127 (delta2[2]*delta2[2])/(sig[2]*sig[2])+
128 (delta2[3]*delta2[3])/(sig[3]*sig[3]);
129
130 //combinatorics
131 if (fDebug && delta2[0] > -900.) hdAx_tr_mc_comb->Fill(delta2[0]);
132 if (fDebug && delta2[1] > -900.) hdX_tr_mc_comb ->Fill(delta2[1]);
133 if (fDebug && delta2[2] > -900.) hdAy_tr_mc_comb->Fill(delta2[2]);
134 if (fDebug && delta2[3] > -900.) hdY_tr_mc_comb ->Fill(delta2[3]);
135
136 if (fDebug) cout<<" imc "<<itr<<" iup "<<iup<<" Nhits "<<vecUp.at(iup).Nhits<<" ax_mc "<<points.at(itr).param[0]<<" ay_mc "<<points.at(itr).param[2]
137 <<" ax_up "<<vecUp.at(iup).param[0]<<" ay_up "<<vecUp.at(iup).param[2]<<
138 " dmc_match "<< dmc_match[itr] << " "<<dmatch<<" delta2[0] "<<delta2[0]<<" delta2[1] "<<delta2[1]<<" delta2[2] "<<delta2[2]<<" delta2[3] "<<delta2[3]<<endl;
139
140 if ( dmc_match[itr] > dmatch){
141 dmc_match[itr] = dmatch;
142 mc_tr_assoc[itr] = iup;
143 //vecUp.at(iup).Pdg = points.at(itr).Pdg;
144 dax = delta2[0];
145 dx = delta2[1];
146 day = delta2[2];
147 dy = delta2[3];
148 }
149
150 }//vecUp.size()
151 if (fDebug){
152 // cout<<" itr "<<itr<<" mc_Id "<<points.at(itr).Id<<
153 // " ax_mc "<<points.at(itr).paramSi[0]<<" reco_ind "<<mc_tr_assoc[itr]<<" ax "<< vecUp.at(mc_tr_assoc[itr]).param[0] <<
154 // " dmc_match "<<dmc_match[itr]<<endl;
155 if (mc_tr_assoc[itr] > -1){
156 cout<<" mc itr "<<itr<<" rec "<<mc_tr_assoc[itr]<<endl;
157 if (dax > -900.) hdAx_uptr_mc->Fill(dax);
158 if (dx > -900.) hdX_uptr_mc ->Fill(dx);
159 if (day > -900.) hdAy_uptr_mc->Fill(day);
160 if (dy > -900.) hdY_uptr_mc ->Fill(dy);
161
162 }//if mc_tr_assoc[itr] > -1
163 }//if
164 }//points.size
165
166 if (fDebug) cout<<"reject poorly chosen association segments "<<endl;
167 for (size_t itr = 0; itr < points.size(); itr++) {//mc_tr
168 if (mc_tr_assoc[itr] == -1) continue;
169
170 for (size_t itr2 = 0; itr2 < points.size(); itr2++) {//mc_tr
171 if (itr2 == itr) continue;
172 if (mc_tr_assoc[itr2] == -1) continue;
173
174 if (mc_tr_assoc[itr] == mc_tr_assoc[itr2]){
175 if (dmc_match[itr2] > dmc_match[itr] ) mc_tr_assoc[itr2] = -1;
176 else {
177 mc_tr_assoc[itr] = -1;
178 break;
179 }
180 }
181 }//itr2
182 //---MC Eff ---
183 //---Num
184 if (fDebug) cout<<" mc_Id "<<points.at(itr).Id<<" pdg "<<points.at(itr).Pdg<<" assoc "<<mc_tr_assoc[itr]<<endl;
185 if (fDebug && mc_tr_assoc[itr] > -1 && fDebug && points.at(itr).Np >=7 //>= 10
186 && points.at(itr).np_3si == 0
187 && (points.at(itr).Np_ch2 > 3 || points.at(itr).Np_ch3 > 3)){
188 if ( points.at(itr).xWas3 && points.at(itr).uWas3 && points.at(itr).vWas3 &&
189 points.at(itr).xWas2 && points.at(itr).uWas2 && points.at(itr).vWas2 )hNum_mcuptr->Fill(1);
190
191 hNum_mcuptr->Fill(0);
192 Ngood_reco_tracks++;
193 //cout<<" itr "<<itr<<" Ngood_reco_tracks "<<Ngood_reco_tracks<<endl;
194 cout<<"UpNum"<<endl;
195 }
196 }//itr
197 if (fDebug) cout<<" Ngood_mc_tracks "<<Ngood_mc_tracks<<" Ngood_reco_tracks "<<Ngood_reco_tracks<<endl;
198 if (fDebug) {
199 hNtr_reco->Fill(Ngood_reco_tracks);
200 hNtr_mc ->Fill(Ngood_mc_tracks);
201 hNtr_mc_vs_reco ->Fill(Ngood_reco_tracks,Ngood_mc_tracks);
202 if (Ngood_mc_tracks == 2) hNrecoTrif2mc->Fill(Ngood_reco_tracks);
203 }
204 //if (fDebug && points.size() == 2) hNrecoTrif2mc->Fill(vecUp.size());
205 for (size_t itr = 0; itr < points.size(); itr++) {//mc_tr
206 if (mc_tr_assoc[itr] > -1) vecUp.at(mc_tr_assoc[itr]).Pdg = points.at(itr).Pdg;
207 }
208
209 twotracks tt;
210 Bool_t track1_was = 0;
211 Bool_t ifNmc2 = 0;
212 Bool_t ifSelectedCases = 0;
213 if (fDebug) cout<<" reco tracks: "<<endl;
214
215 if (Ngood_mc_tracks == 2 ) ifNmc2 = 1;
216
217 if ( vecUp.size() == 2 && Ngood_mc_tracks == 2){
218 for (size_t iup= 0; iup < vecUp.size(); iup ++) {//reco
219 if (fDebug) cout<<" iup "<<iup<< " found code pdg "<<vecUp.at(iup).Pdg <<endl;
220 if (track1_was == 0){
221 track1_was =1;
222 tt.track1_pdg = vecUp.at(0).Pdg;
223 tt.track1[0] = vecUp.at(iup).param[0];
224 tt.track1[1] = vecUp.at(iup).param[1];
225 tt.track1[2] = vecUp.at(iup).param[2];
226 tt.track1[3] = vecUp.at(iup).param[3];
227 }
228 if (track1_was == 1){
229 tt.track2_pdg = vecUp.at(1).Pdg;
230 tt.track2[0] = vecUp.at(iup).param[0];
231 tt.track2[1] = vecUp.at(iup).param[1];
232 tt.track2[2] = vecUp.at(iup).param[2];
233 tt.track2[3] = vecUp.at(iup).param[3];
234 }
235 }//reco
236 if (track1_was == 1) vec_twotracks.push_back(tt);
237 if (fDebug) cout<<" vec_twotracks "<<vec_twotracks.size()<<endl;
238
239 //for(int itr = 0; itr< vec_twotracks.size(); ++itr){
240 if (vec_twotracks.size() == 1){
241 Int_t itr = 0;
242 if (
243 (vec_twotracks.at(itr).track1_pdg == PDG_Be7 && vec_twotracks.at(itr).track2_pdg == PDG_p) || (vec_twotracks.at(itr).track2_pdg == PDG_Be7 && vec_twotracks.at(itr).track1_pdg == PDG_p) ||
244 (vec_twotracks.at(itr).track1_pdg == PDG_Be7 && vec_twotracks.at(itr).track2_pdg == PDG_H2) || (vec_twotracks.at(itr).track2_pdg == PDG_Be7 && vec_twotracks.at(itr).track1_pdg == PDG_H2) ||
245 (vec_twotracks.at(itr).track1_pdg == PDG_Be7 && vec_twotracks.at(itr).track2_pdg == PDG_He4) || (vec_twotracks.at(itr).track2_pdg == PDG_Be7 && vec_twotracks.at(itr).track1_pdg == PDG_He4)||
246 (vec_twotracks.at(itr).track1_pdg == PDG_Be7 && vec_twotracks.at(itr).track2_pdg == PDG_He3) || (vec_twotracks.at(itr).track2_pdg == PDG_Be7 && vec_twotracks.at(itr).track1_pdg == PDG_He3)||
247 (vec_twotracks.at(itr).track1_pdg == PDG_Li7 && vec_twotracks.at(itr).track2_pdg == PDG_He4) || (vec_twotracks.at(itr).track2_pdg == PDG_Li7 && vec_twotracks.at(itr).track1_pdg == PDG_He4)||
248 (vec_twotracks.at(itr).track1_pdg == PDG_Li7 && vec_twotracks.at(itr).track2_pdg == PDG_He3) || (vec_twotracks.at(itr).track2_pdg == PDG_Li7 && vec_twotracks.at(itr).track1_pdg == PDG_He3)||
249 (vec_twotracks.at(itr).track1_pdg == PDG_Li6 && vec_twotracks.at(itr).track2_pdg == PDG_He4) || (vec_twotracks.at(itr).track2_pdg == PDG_Li6 && vec_twotracks.at(itr).track1_pdg == PDG_He4)||
250 (vec_twotracks.at(itr).track1_pdg == PDG_Li6 && vec_twotracks.at(itr).track2_pdg == PDG_He3) || (vec_twotracks.at(itr).track2_pdg == PDG_Li6 && vec_twotracks.at(itr).track1_pdg == PDG_He3)||
251 (vec_twotracks.at(itr).track1_pdg == PDG_He4 && vec_twotracks.at(itr).track2_pdg == PDG_He4)
252
253 ) ifSelectedCases = 1;
254
255 if(fDebug) cout<<" pdg1 "<<vec_twotracks.at(itr).track1_pdg <<" pdg2 "<<vec_twotracks.at(itr).track2_pdg<<" ifSelectedCases "<<ifSelectedCases<<endl;
256 Double_t CosA = (vec_twotracks.at(itr).track1[0]*vec_twotracks.at(itr).track2[0] + vec_twotracks.at(itr).track1[2]*vec_twotracks.at(itr).track2[2] + 1)/
257 (sqrt((vec_twotracks.at(itr).track1[0]*vec_twotracks.at(itr).track1[0] + vec_twotracks.at(itr).track1[2]*vec_twotracks.at(itr).track1[2] + 1)*
258 (vec_twotracks.at(itr).track2[0]*vec_twotracks.at(itr).track2[0] + vec_twotracks.at(itr).track2[2]*vec_twotracks.at(itr).track2[2] + 1)));
259 Double_t alpha = TMath::ACos(CosA);
260 Double_t alphaDeg = TMath::RadToDeg()* alpha;
261 if(fDebug) cout<<" alpha "<<alpha<<" CosA "<<CosA<<" ToDeg "<<alphaDeg<<endl;
262 hAngle_reco->Fill(alphaDeg);
263 if (ifNmc2) hAngle_recoifNmc2->Fill(alphaDeg);
264 if (ifNmc2 && ifSelectedCases) hAngle_recoifNmc2Cases->Fill(alphaDeg);
265 }
266
267 }//Ngood_mc_tracks == 2
268
269
270}//MCefficiencyCalculation
271
272
273//-----------------------------ReadSiliconTracks------------------------
274void BmnUpstreamTracking::ReadSiliconTracks(Double_t** par_ab, Double_t* par_z, Int_t & NTracks, vector<MC_points> & vec){
275 MC_points tmpTr;
276 if (!expData){
277 if (fDebug) cout<<" MC_SiliconTracks:"<<endl;
278 int tr_before = -1;
279 int Nmc_tracks = -1;
280
281 int Nst_mc[kMaxMC];
282 int Np_in_3st[kMaxMC];
283 int mcTracksArray[kMaxMC];
284 double Xmc[kMaxMC][fNstations];
285 double Ymc[kMaxMC][fNstations];
286 double Zmc[kMaxMC][fNstations];
287 for (Int_t Id= 0; Id < kMaxMC; Id++) {
288 Nst_mc[Id] = 0;
289 Np_in_3st[Id] = 0;
290 mcTracksArray[Id] = -1;
291 for (Int_t i = 0; i < fNstations; i++) {
292 Xmc[Id][i] = -999.;
293 Ymc[Id][i] = -999.;
294 Zmc[Id][i] = -999.;
295 }
296 }
297
298 for (Int_t iMC = 0; iMC < fSiTracksSim->GetEntriesFast(); ++iMC) {
299 BmnSiliconHit* hit = (BmnSiliconHit*)fSiTracksSim->UncheckedAt(iMC);
300
301 Double_t x_MC = hit->GetX();
302 Double_t y_MC = hit->GetY();
303 Double_t z_MC = hit->GetZ();
304 Int_t trackId_MC= hit->GetIndex();
305 Int_t st_mc = -1;
306
307 if (tr_before != trackId_MC) {
308 Nmc_tracks++;
309 mcTracksArray[Nmc_tracks] = trackId_MC;
310 }
311 tr_before = trackId_MC;
312
313 if (z_MC > -320.) {st_mc = 3; Np_in_3st[Nmc_tracks]++; }
314 if (z_MC < -434. && z_MC > -436. ) st_mc = 2;
315 if (z_MC < -438.) st_mc = 1;
316 if (fDebug) cout<<" Id "<<trackId_MC<<" x "<<x_MC<<" y "<<y_MC<<" z "<<z_MC<<" st_mc "<<st_mc<<endl;
317
318 Xmc[Nmc_tracks][st_mc] = x_MC;
319 Ymc[Nmc_tracks][st_mc] = y_MC;
320 Zmc[Nmc_tracks][st_mc] = z_MC;
321 Nst_mc[Nmc_tracks]++;
322
323 }//iMC
324 Nmc_tracks++;
325
326 for (Int_t Id= 0; Id < kMaxMC; Id++) {
327
328 for (Int_t i = 0; i < fNstations; i++) {
329 tmpTr.x_si[i] = Xmc[Id][i];
330 tmpTr.y_si[i] = Ymc[Id][i];
331 tmpTr.z_si[i] = Zmc[Id][i];
332 }
333
334 tmpTr.Id = mcTracksArray[Id];
335 tmpTr.Np_Si = Nst_mc[Id];
336 if (Np_in_3st[Id] == tmpTr.Np_Si ) tmpTr.np_3si = 1;
337 if (fDebug && Np_in_3st[Id] > 0 ) cout<<" Np_in_3st "<<Np_in_3st[Id] <<" tmpTr.Np_Si "<<tmpTr.Np_Si<<" tmpTr.np_3si "<<tmpTr.np_3si<<endl;
338 if (Nst_mc[Id] >= 2 && Np_in_3st[Id] > 0) vec.push_back(tmpTr);
339 }
340
341 if (fDebug) cout<<" vec_points.size() "<<vec.size()<<endl;
342
343 }// !expData
344 if ( fDebug ) cout<<endl;
345
346
347 if (fDebug) cout<<"----Reco_SiliconTracks:"<<endl;
348 for (Int_t iTrack = 0; iTrack < fSiTracks->GetEntriesFast(); iTrack++) {
349 BmnTrack* track = (BmnTrack*) fSiTracks ->At(iTrack);
350
351 FairTrackParam* par = track->GetParamLast();
352 //FairTrackParam* par = track->GetParamFirst();
353
354 if ( fDebug ) cout<<" SiTr "<<iTrack<<" SiTx "<<par->GetTx()<<" SiX "<<par->GetX()<<" SiTy "<<par->GetTy()<<" SiY "<<par->GetY()<<" SiZ "<<par->GetZ()<<" Nhits "<<track->GetNHits()<<endl;
355
356 par_ab[0][NTracks]= par->GetTx();
357 par_ab[1][NTracks]= par->GetX();
358 par_ab[2][NTracks]= par->GetTy();
359 par_ab[3][NTracks]= par->GetY();
360 par_ab[4][NTracks]= track->GetNHits();
361 par_z[NTracks] = par->GetZ();
362 NTracks++;
363 }
364
365 if ( fDebug ) cout<<endl;
366
367}//ReadSiliconTracks
368//----------------------------------------------------------------------
369
370//----------------------------------------------------------------------
371void BmnUpstreamTracking::ReadSiliconHits(Double_t*** hits, Int_t* NSihits){
372 //-------------Read Silicon Hits------------------------------------
373 if ( fDebug ) cout<<" Hits from Tracks "<<endl;
374 for (Int_t iSihit = 0; iSihit < fSilHits->GetEntriesFast(); iSihit++) {
375 BmnSiliconHit* hit = (BmnSiliconHit*) fSilHits ->At(iSihit);
376
377 // if ( fDebug ) cout<<" iSihit "<<iSihit<<endl;
378 Int_t st = hit->GetStation();
379
380 if ( fDebug ) cout<<" st "<<st<<" hit: X "<<hit->GetX()<<" Y "<<hit->GetY()<<" Sigx "<<hit->GetDx()<<" Sigy "<<hit->GetDy()<<endl;
381 hits[st][NSihits[st]][0] = hit->GetX();
382 hits[st][NSihits[st]][1] = hit->GetY();
383 hits[st][NSihits[st]][2] = hit->GetZ();
384 hits[st][NSihits[st]][3] = hit->GetDx();
385 hits[st][NSihits[st]][4] = hit->GetDy();
386 hits[st][NSihits[st]][5] = hit->GetModule();
387
388 NSihits[st]++;
389 }
390}//ReadSiliconHits
391//----------------------------------------------------------------------
392
393//----------------------------------------------------------------------
394void BmnUpstreamTracking::ReadMWPCSegments(Double_t*** par_ab, Double_t** par_z, Int_t* Nseg, vector<MC_points> & vec){
395 for (Int_t iSeg = 0; iSeg < fMWPCSegments->GetEntries() ; ++iSeg) {
396 BmnTrack* segment = (BmnTrack*) fMWPCSegments ->At(iSeg) ;
397
398 Int_t ich = -1;
399 Double_t Z_mwpc = segment->GetParamFirst()->GetZ();
400 //Double_t Chi2 = segment->GetChi2();
401
402 if ( Z_mwpc >= -859. && Z_mwpc <= -857.) ich = 0;
403 if ( Z_mwpc >= -760. && Z_mwpc <= -758.) ich = 1;
404 if ( Z_mwpc >= -365. && Z_mwpc <= -363.) ich = 2;
405 if ( Z_mwpc >= -215. && Z_mwpc <= -213.) ich = 3;
406
407 if ( ich < 0) continue;
408
409 if ( fDebug ) cout<<" MWPCSeg: ich "<<ich<<" Tx "<<segment->GetParamFirst()->GetTx()<<" X "<<segment->GetParamFirst()->GetX()<<" Ty "<<segment->GetParamFirst()->GetTy()<<" Y "<<segment->GetParamFirst()->GetY()<< " Z "<<Z_mwpc<<" Nhits "<<segment->GetNHits()<<endl;
410 par_ab[ich][0][Nseg[ich]]= segment->GetParamFirst()->GetTx();
411 par_ab[ich][1][Nseg[ich]]= segment->GetParamFirst()->GetX();
412 par_ab[ich][2][Nseg[ich]]= segment->GetParamFirst()->GetTy();
413 par_ab[ich][3][Nseg[ich]]= segment->GetParamFirst()->GetY();
414 par_ab[ich][4][Nseg[ich]]= segment->GetNHits();
415 par_ab[ich][5][Nseg[ich]]= segment->GetParamFirst()->GetCovariance(0, 0);
416 par_ab[ich][6][Nseg[ich]]= segment->GetParamFirst()->GetCovariance(1, 1);
417 par_ab[ich][7][Nseg[ich]]= segment->GetParamFirst()->GetCovariance(2, 2);
418 par_ab[ich][8][Nseg[ich]]= segment->GetParamFirst()->GetCovariance(3, 3);
419 par_z[ich][Nseg[ich]] = segment->GetParamFirst()->GetZ();
420
421 if ( ich > 1 ){
422 par_ab[ich][1][Nseg[ich]]+= -X_shift - X_shift_seg[ich];
423 par_ab[ich][3][Nseg[ich]]+= -Y_shift;
424 }
425 Nseg[ich]++;
426
427 }
428
429 if (!expData){
430 if (fDebug)cout<<"----MC: MWPCSegments"<<endl;
431
432 Int_t mcTracksArray[kMaxMC];
433 Int_t mcTracksPdg[kMaxMC];
434 Double_t X2mc[kMaxMC][kNPlanes];
435 //Double_t Y2mc[kMaxMC][kNPlanes];
436 //Double_t Z2mc[kMaxMC][kNPlanes];
437 Double_t X3mc[kMaxMC][kNPlanes];
438 //Double_t Y3mc[kMaxMC][kNPlanes];
439 //Double_t Z3mc[kMaxMC][kNPlanes];
440 int Npl_MC2[kMaxMC]; int Npl_MC3[kMaxMC];
441
442 for (Int_t Id= 0; Id < kMaxMC; Id++) {
443 Npl_MC2[Id] = 0;
444 Npl_MC3[Id] = 0;
445 mcTracksArray[Id] = -1;
446 mcTracksPdg[Id] = -1;
447 for (Int_t i = 0; i < kNPlanes; i++) {
448 X2mc[Id][i] = -999.;
449 //Y2mc[Id][i] = -999.;
450 //Z2mc[Id][i] = -999.;
451 X3mc[Id][i] = -999.;
452 //Y3mc[Id][i] = -999.;
453 //Z3mc[Id][i] = -999.;
454 }
455 }
456
457 int tr_before = -1;
458 int Nmc_tracks = -1;
459
460
461 for (Int_t iMC = 0; iMC < fBmnHitsArray->GetEntriesFast(); ++iMC) {
462 BmnMwpcHit* hit = (BmnMwpcHit*)fBmnHitsArray->UncheckedAt(iMC);
463
464 Int_t st_MC = hit->GetMwpcId();
465 Int_t trackId_MC = hit->GetHitId();
466 //Int_t pdg_MC = hit->GetPdgId();
467 Int_t pl_MC = hit->GetPlaneId();
468 //Short_t wire_MC = hit->GetWireNumber();
469 //Double_t time_MC = hit->GetWireTime();
470
471 //if (fDebug)cout<<" st_MC "<<st_MC<<" trackId_MC "<<trackId_MC<<" pl_MC "<<pl_MC<<" X "<<hit->GetX()<<" wire_MC "<<wire_MC<<endl;
472
473 if (tr_before != trackId_MC) {
474 Nmc_tracks++;
475 mcTracksArray[Nmc_tracks] = hit->GetHitId();
476 mcTracksPdg[Nmc_tracks] = hit->GetPdgId();
477 }
478 tr_before = trackId_MC;
479
480 if (st_MC == 2){
481 X2mc[Nmc_tracks][pl_MC] = hit->GetX();
482 //Y2mc[Nmc_tracks][pl_MC] = hit->GetY();
483 //Z2mc[Nmc_tracks][pl_MC] = hit->GetZ();
484 Npl_MC2[Nmc_tracks]++;
485 }
486 if (st_MC == 3){
487 X3mc[Nmc_tracks][pl_MC] = hit->GetX();
488 //Y3mc[Nmc_tracks][pl_MC] = hit->GetY();
489 //Z3mc[Nmc_tracks][pl_MC] = hit->GetZ();
490 Npl_MC3[Nmc_tracks]++;
491
492 }
493 //if (fDebug)cout<<" X2["<<Nmc_tracks<<"]["<<pl_MC<<"] "<<X2mc[Nmc_tracks][pl_MC]<<" Npl_MC2 "<< Npl_MC2[Nmc_tracks]<<" X3["<<Nmc_tracks<<"]["<<pl_MC<<"] "<<X3mc[Nmc_tracks][pl_MC]<<" Npl_MC3 "<< Npl_MC3[Nmc_tracks]<<" z3 "<<Z3mc[Nmc_tracks][pl_MC]<<endl;
494
495 }//iMC
496 Nmc_tracks++;
497
498
499 if (fDebug)cout<<" Nmc_tracks "<<Nmc_tracks<<" MC vec_points.size() "<<vec.size()<<endl;
500 MC_points tmpTr;
501 //bool mc_new = 0;
502 for (Int_t id = 0; id < Nmc_tracks; id++) { //mc_tr_curr
503 if (fDebug)cout<<" id "<<id<<" Id_mc "<< mcTracksArray[id]<<" pdg "<<mcTracksPdg[id]<<" Npl2 "<<Npl_MC2[id]<<" Npl3 "<<Npl_MC3[id]<<endl;
504
505 for (size_t itr = 0; itr < vec.size(); itr++) {// mc_tr_before
506
507 if (Npl_MC2[id] >= 4 || Npl_MC3[id] >= 4){
508 if (vec.at(itr).Id == mcTracksArray[id]) {
509 vec.at(itr).Pdg = mcTracksPdg[id];
510 vec.at(itr).Np_ch2 = Npl_MC2[id];
511 vec.at(itr).Np_ch3 = Npl_MC3[id];
512 vec.at(itr).Np_P = Npl_MC2[id] + Npl_MC3[id];
513 vec.at(itr).Np = vec.at(itr).Np_P + vec.at(itr).Np_Si;
514
515 int i3 = -1;
516 if (X3mc[id][5] > -900.) i3 =5;
517 if (i3 < 0 && X3mc[id][4] > -900.) i3 =4;
518 if (i3 < 0 && X3mc[id][3] > -900.) i3 =3;
519
520 int isi = -1;
521 if (isi < 0 && vec.at(itr).x_si[1] > -900.) isi = 1;
522 if (isi < 0 && vec.at(itr).x_si[2] > -900.) isi = 2;
523
524
525 vec.at(itr).param[0] = (vec.at(itr).x_si[isi] - vec.at(itr).x_si[3])/(vec.at(itr).z_si[isi] - vec.at(itr).z_si[3]);
526 vec.at(itr).param[1] = vec.at(itr).param[0]*(-350 - vec.at(itr).z_si[3]) + vec.at(itr).x_si[3];
527 vec.at(itr).param[2] = (vec.at(itr).y_si[isi] - vec.at(itr).y_si[3])/(vec.at(itr).z_si[isi] - vec.at(itr).z_si[3]);
528 vec.at(itr).param[3] = vec.at(itr).param[2]*(-350 - vec.at(itr).z_si[3]) + vec.at(itr).y_si[3];
529
530 if (fDebug)cout<<" si mc: x3 "<<vec.at(itr).x_si[3]<<" z3 "<<vec.at(itr).z_si[3]<<" UpX "<<vec.at(itr).param[1]<<endl;
531
532 if (X2mc[id][0] > -900. || X2mc[id][3] > -900.) vec.at(itr).xWas2 = 1;
533 if (X2mc[id][1] > -900. || X2mc[id][4] > -900.) vec.at(itr).vWas2 = 1;
534 if (X2mc[id][2] > -900. || X2mc[id][5] > -900.) vec.at(itr).uWas2 = 1;
535
536 if (X3mc[id][0] > -900. || X3mc[id][3] > -900.) vec.at(itr).xWas3 = 1;
537 if (X3mc[id][1] > -900. || X3mc[id][4] > -900.) vec.at(itr).vWas3 = 1;
538 if (X3mc[id][2] > -900. || X3mc[id][5] > -900.) vec.at(itr).uWas3 = 1;
539
540
541 }//vec_points.at(itr).Id == mcTracksArray[id]
542 }//Npl_MC2[id] >= 4 || Npl_MC3[id] >= 4
543
544 }//vec_points.size
545 }//Nmc_tracksId = -1;
546
547 if (fDebug){
548 cout<<endl;
549 cout<<"----MCtrue all systems, size: "<<vec.size()<<endl;
550 for (size_t itr = 0; itr < vec.size(); itr++) {
551 cout<<" itr "<<itr<<" Id "<<vec.at(itr).Id<<" Np_Si "<<vec.at(itr).Np_Si<<" Np_P "<<vec.at(itr).Np_P<<" Np_ch2 "<<vec.at(itr).Np_ch2<<" Np_ch3 "<<vec.at(itr).Np_ch3<<endl;
552 cout<<" itr "<<itr<<" Id "<<vec.at(itr).Id<<" Up: ax "<<vec.at(itr).param[0]<<" bx "<<vec.at(itr).param[1]<<
553 " ay "<<vec.at(itr).param[2]<<" by "<<vec.at(itr).param[3]<<endl;
554
555 cout<<" Xt "<<vec.at(itr).param[0]*( kZ_target - Zcentr ) + vec.at(itr).param[1]<<
556 " Yt "<<vec.at(itr).param[2]*( kZ_target - Zcentr ) + vec.at(itr).param[3]<<endl;
557
558 hAx_upmc->Fill(vec.at(itr).param[0]);
559 hX_upmc ->Fill(vec.at(itr).param[1]);
560 hAy_upmc->Fill(vec.at(itr).param[2]);
561 hY_upmc ->Fill(vec.at(itr).param[3]);
562 }
563 }//if
564
565 }// !expData
566 if (fDebug)cout<<endl;
567}//ReadMWPCSegments
568//----------------------------------------------------------------------
569
570void BmnUpstreamTracking::ReadMWPCTracks(Double_t*** par_ab, Double_t** par_z, Int_t* NTracks){
571 for (Int_t iTr = 0; iTr < fMWPCTracks->GetEntriesFast() ; ++iTr) {
572 BmnTrack* track = (BmnTrack*) fMWPCTracks->At(iTr) ;
573
574 Int_t ip = -1;
575 Double_t Z = track->GetParamFirst()->GetZ();
576
577 if ( Z >= -809. && Z <= -807.) ip = 0;
578 if ( Z >= -290. && Z <= -288.) ip = 1;
579 if ( ip < 0) continue;
580
581 if ( fDebug ) cout<<" MWPCTr: ip "<<iTr<<" Tx "<<track->GetParamFirst()->GetTx()<<" X "<<track->GetParamFirst()->GetX()<<" Ty "<<track->GetParamFirst()->GetTy()<<" Y "<<track->GetParamFirst()->GetY()<<" Z "<<Z<<endl;
582
583 par_ab[ip][0][NTracks[ip]] = track->GetParamFirst()->GetTx() ;
584 par_ab[ip][1][NTracks[ip]] = track->GetParamFirst()->GetX() ;
585 par_ab[ip][2][NTracks[ip]] = track->GetParamFirst()->GetTy() ;
586 par_ab[ip][3][NTracks[ip]] = track->GetParamFirst()->GetY() ;
587 par_ab[ip][4][NTracks[ip]] = track->GetNHits();
588 par_ab[ip][5][NTracks[ip]] = track->GetParamFirst()->GetCovariance(0, 0);
589 par_ab[ip][6][NTracks[ip]] = track->GetParamFirst()->GetCovariance(1, 1);
590 par_ab[ip][7][NTracks[ip]] = track->GetParamFirst()->GetCovariance(2, 2);
591 par_ab[ip][8][NTracks[ip]] = track->GetParamFirst()->GetCovariance(3, 3);
592 par_z[ip][NTracks[ip]] = track->GetParamFirst()->GetZ();
593
594 par_ab[1][1][NTracks[ip]]+= -X_shift;
595 par_ab[1][3][NTracks[ip]]+= -Y_shift;
596 /*
597 hAxPCtr->Fill(par_ab[1][0][NTracks[ip]]);
598 hxPCtr->Fill(par_ab[1][1][NTracks[ip]]);
599 hAyPCtr->Fill(par_ab[1][2][NTracks[ip]]);
600 hyPCtr->Fill(par_ab[1][3][NTracks[ip]]);
601 hbeam_profile_pair_[ip]->Fill(track->GetParamFirst()->GetX(),track->GetParamFirst()->GetY());
602 */
603 NTracks[ip]++;
604
605 }//itr
606
607}//ReadMWPCTracks
608//----------------------------------------------------------------------
609
610
611//----------------------------------------------------------------------
612void BmnUpstreamTracking::PCTracksSiTracksMatching(Double_t** par_Si, Double_t* par_Siz, Int_t & NTracksSi,
613Double_t*** par_PC, Double_t** par_PCz, Int_t* NTracksPC, vector<Smatch> & vec_m, vector<Smatch> & Out){
614
615 Smatch tmpSeg;
616 vec_m.clear();
617
618 if (fDebug) cout<<" NTracksSi "<<NTracksSi<<" NPCTracks[1] "<<NPCTracks[1]<<endl;
619
620 for (Int_t iSi=0; iSi < NTracksSi; ++iSi){
621 for (Int_t iPCtr = 0; iPCtr < NPCTracks[1]; ++iPCtr){
622
623 Double_t dz = par_PCz[1][iPCtr] - par_Siz[iSi];
624
625 Double_t dAx = par_Si[0][iSi] - par_PC[1][0][iPCtr] ;
626 Double_t dAy = par_Si[2][iSi] - par_PC[1][2][iPCtr] ;
627 Double_t dX = (par_Si[0][iSi]* dz + par_Si[1][iSi]) - par_PC[1][1][iPCtr];
628 Double_t dY = (par_Si[2][iSi]* dz + par_Si[3][iSi]) - par_PC[1][3][iPCtr];
629 /*
630 hdX_MWPCpair1_Sitr->Fill(dX);
631 hdAx_MWPCpair1_Sitr->Fill(dAx);
632 hdY_MWPCpair1_Sitr->Fill(dY);
633 hdAy_MWPCpair1_Sitr->Fill(dAy);
634 */
635 if (fDebug) cout<<" Si("<<iSi<<")-P1("<<iPCtr<<"): dAx("<<cutAx<<") "<<dAx<<" dAy("<<cutAy<<") "<<dAy<<" dX("<<cutX<<") "<<dX<<" dY("<<cutY<<") "<<dY<<endl;
636 if (fabs(dAx) < cutAx && fabs(dAy) < cutAy && fabs(dX) < cutX && fabs(dY) < cutY ){
637 if (fDebug) cout<<"after cuts: Si "<<iSi<<"- P1 "<<iPCtr<<endl;
638 Double_t Chi2_match = dX*dX + 0.25*dY*dY;
639
640 tmpSeg.Ind1 = iSi;
641 tmpSeg.Ind2 = iPCtr;
642 tmpSeg.Inds2 = -1;
643 tmpSeg.Inds3 = -1;
644 tmpSeg.Nhits1 = par_Si[4][iSi];
645 tmpSeg.Nhits2 = par_PC[1][4][iPCtr];
646 tmpSeg.Z1 = par_Siz[iSi];
647 tmpSeg.Z2 = par_PCz[1][iPCtr];
648 tmpSeg.distX = dX;
649 tmpSeg.distY = dY;
650 tmpSeg.Chi2m = Chi2_match;
651 // if (fDebug) cout<<" iSi "<<iSi<<" iPCtr "<<iPCtr<<" dX "<<dx<<" dY "<<dy<<" Chi2_match "<<Chi2_match<<endl;
652
653 for(int ipar = 0; ipar < 4; ipar++){
654 tmpSeg.param1[ipar] = par_Si[ipar][iSi];
655 tmpSeg.param2[ipar] = par_PC[1][ipar][iPCtr];
656 }
657 vec_m.push_back(tmpSeg);
658
659 /*
660 hdX_MWPCpair1_Sitr_cut->Fill(dX);
661 hdAx_MWPCpair1_Sitr_cut->Fill(dAx);
662 hdY_MWPCpair1_Sitr_cut->Fill(dY);
663 hdAy_MWPCpair1_Sitr_cut->Fill(dAy);
664 */
665 }//if
666 }//iPCtr
667 }//iSi
668 if (fDebug) cout<<" --End PCtracks: size "<<vtmpSeg.size()<<endl;
669
670 if (vtmpSeg.size() > 0){
671 Int_t PCtrack = 1;
672 RecordingBest(PCtrack, vec_m, Out);
673 }
674
675 if (fDebug) cout<<"PCTracksSiTracksMatching: OutVector.size() "<<Out.size()<<endl;
676 for (size_t iter = 0; iter < Out.size(); ++iter) {
677 if (fDebug) printf(" OutVector.at(%zu): %8.4f | Si %d - PC %d Seg2 %d - Seg3 %d\n", iter, Out.at(iter).Chi2m, Out.at(iter).Ind1, Out.at(iter).Ind2, Out.at(iter).Inds2, Out.at(iter).Inds3 );
678 }
679}//PCTracksSiTracksMatching
680//----------------------------------------------------------------------
681
682//----------------------------------------------------------------------
683void BmnUpstreamTracking::PCSegmentsSiTracksMatching(Double_t** par_Si, Double_t* par_Siz, Int_t & NTracksSi,
684Double_t*** par_PC, Double_t** par_PCz, Int_t* NSegPC, vector<Smatch> & vec_m, vector<Smatch> & Out){
685
686 vec_m.clear();
687 if ( fDebug ) cout<<endl;
688 if ( fDebug ) cout<<" ----Matching: PC segments Ch2 ---"<<endl;
689 if ( fDebug ) cout<<" NTracksSi "<<NTracksSi<<" NSegPC[2] "<<NSegPC[2]<<endl;
690 for (int iSi=0; iSi < NTracksSi; ++iSi){
691 for(int iseg2 = 0; iseg2 < NSegPC[2]; ++iseg2){
692
693 Double_t dz2 = par_PCz[2][iseg2] - par_Siz[iSi];
694 // if ( fDebug ) cout<<" Zpc "<<par_PCz[2][iseg2]<<" Zsi "<<par_Siz[iSi]<<endl;
695 Double_t X_extr2 = par_Si[0][iSi]* dz2 + par_Si[1][iSi];
696 Double_t Y_extr2 = par_Si[2][iSi]* dz2 + par_Si[3][iSi];
697
698 Double_t dX = X_extr2 - par_PC[2][1][iseg2];
699 Double_t dY = Y_extr2 - par_PC[2][3][iseg2];
700
701 //hdX_MWPCSeg2_Sitr->Fill(dX);
702 //hdY_MWPCSeg2_Sitr->Fill(dY);
703
704 if (fDebug) cout<<" Si["<<iSi<<"]-Seg2["<<iseg2<<"] : dX("<<cutX<<") "<<dX<<" dY("<<cutY<<") "<<dY<<endl;
705 if (fabs(dX) < cutX && fabs(dY) < cutY ){
706
707 Int_t iseg3 = -1;
708 if (fDebug) cout<<" cuts ok:iseg2 "<<iseg2<<endl;
709 RecordCandidateSeg(iSi, iseg2, iseg3, par_Si, par_PC, par_Siz, par_PCz, dX, dY, vec_m);
710
711 }//if
712 }//iseg2
713 }//iSi
714
715 if (fDebug) cout<<" --End PCiseg2 : now vec_m.size "<<vec_m.size()<<endl;
716 if (fDebug) cout<<endl;
717 //-----------------------------------------------------------------
718 if ( fDebug ) cout<<endl;
719 if ( fDebug ) cout<<" ----Matching: PC segments Ch3 ---"<<endl;
720 if ( fDebug ) cout<<" NTracksSi "<<NTracksSi<<" NSegPC[3] "<<NSegPC[3]<<endl;
721 for (int iSi=0; iSi < NSiTracks; ++iSi){
722 for(int iseg3 = 0; iseg3 < NSegPC[3]; ++iseg3){
723
724 Double_t dz3 = par_PCz[3][iseg3] - par_Siz[iSi];
725 Double_t X_extr3 = par_Si[0][iSi]* dz3 + par_Si[1][iSi];
726 Double_t Y_extr3 = par_Si[2][iSi]* dz3 + par_Si[3][iSi];
727
728 Double_t dX = X_extr3 - par_PC[3][1][iseg3];
729 Double_t dY = Y_extr3 - par_PC[3][3][iseg3];
730 //hdX_MWPCSeg3_Sitr->Fill(dX);
731 //hdY_MWPCSeg3_Sitr->Fill(dY);
732
733 if (fDebug) cout<<" Si["<<iSi<<"]-Seg3["<<iseg3<<"] : dX("<<cutX<<") "<<dX<<" dY("<<cutY<<") "<<dY<<endl;
734 if (fabs(dX) < cutX && fabs(dY) < cutY ){
735
736 Int_t iseg2 = -1;
737 RecordCandidateSeg(iSi, iseg2, iseg3, par_Si, par_PC, par_Siz, par_PCz, dX, dY, vec_m);
738
739 }//if
740 }//iseg3
741 }//iSi
742
743 if (fDebug) cout<<" --End PCiseg3 : now vec_m.size "<<vec_m.size()<<endl;
744
745 if (vec_m.size() > 0){
746 Int_t PCtrack = -1;
747 RecordingBest(PCtrack, vec_m, Out);
748 }
749 if (fDebug) cout<<" vec_m.size() "<<vec_m.size()<<" Out.size() "<<Out.size()<<endl;
750 if (fDebug) cout<<endl;
751}//PCSegmentsSiTracksMatching
752//----------------------------------------------------------------------
753
754//----------------------------------------------------------------------
755void BmnUpstreamTracking::RecordCandidateSeg(Int_t & iSi, Int_t & is2, Int_t & is3,
756Double_t ** par_Si, Double_t*** par_PC,
757Double_t* par_Siz, Double_t** ZPC,
758Double_t & dx, Double_t & dy, vector<Smatch>& vtmp){
759
760 if ( fDebug ) cout<<" RecordCandidateSeg: Inds: seg2 "<<is2<<" seg3 "<<is3<<endl;
761 Smatch tmpSeg;
762 Int_t iseg = is2;
763 Int_t ich = 2;
764 if ( is2 < 0 ) {
765 iseg = is3;
766 ich = 3;
767 }
768
769 Double_t Chi2_match = dx*dx + 0.25*dy*dy;
770
771 tmpSeg.Ind1 = iSi;
772 tmpSeg.Ind2 = -1;
773 tmpSeg.Inds2 = is2;
774 tmpSeg.Inds3 = is3;
775 tmpSeg.Z1 = par_Siz[iSi];
776 tmpSeg.Zs2 = ZPC[ich][iseg];
777 tmpSeg.distX = dx;
778 tmpSeg.distY = dy;
779 tmpSeg.Chi2m = Chi2_match;
780 if ( fDebug ) cout<<" iSi "<<iSi<<" iseg2 "<<is2<<" iseg3 "<<is3<<" dX "<<dx<<" dY "<<dy<<" Chi2_match "<<Chi2_match<<endl;
781
782 for(int ipar = 0; ipar < 4; ipar++) {
783 tmpSeg.param1[ipar] = par_Si[ipar][iSi];
784 tmpSeg.params2[ipar] = par_PC[ich][ipar][iseg];
785 }
786 vtmp.push_back(tmpSeg);
787}
788//----------------------------------------------------------------------
789
790
791//----------------------------------------------------------------------
792void BmnUpstreamTracking::RecordingBest(Int_t & iPCtrack, vector<Smatch> & vec, vector<Smatch> & OutVec){
793 if ( iPCtrack > -1 ){
794 if (fDebug) cout<<endl;
795 if (fDebug) cout<<"--RecordingBest: tracks --"<<endl;
796 // vector sorting
797 sort(vec.begin(), vec.end(), compareSegments);
798
799 //first best
800 OutVec.push_back(vec.at(0));
801 if (fDebug) cout<<"--isMatch: Chi2 "<<vec.at(0).Chi2m<<endl;
802
803 //reject repeat index
804 Bool_t isMatch;
805 for (size_t iter = 0; iter < vec.size(); ++iter) {
806 if (fDebug) printf(" vtmpSeg.at(%zu): %8.4f | %d - %d\n", iter, vec.at(iter).Chi2m, vec.at(iter).Ind1, vec.at(iter).Ind2 );
807 isMatch = 0;
808 for (size_t InIter = 0; InIter < OutVec.size(); ++InIter) {
809 if(vec.at(iter).Ind1 == OutVec.at(InIter).Ind1 || vec.at(iter).Ind2 == OutVec.at(InIter).Ind2) {
810 isMatch = 1;
811 continue;
812 }
813 }
814 //writing unique index
815 if (isMatch == 0) {
816 if (fDebug) cout<<"--isMatch: Chi2 "<<vec.at(iter).Chi2m<<endl;
817 OutVec.push_back(vec.at(iter));
818 }
819 }//iter
820
821 if (fDebug) cout<<" Now OutVec.size "<<OutVec.size()<<endl;
822 vec.clear();
823 }else{//---segs--
824 if (fDebug) cout<<"--RecordingBest: segs --"<<endl;
825 // vector sorting
826 sort(vec.begin(), vec.end(), compareSegments);
827
828 //reject repeat index
829 Bool_t isMatch;
830 for (size_t iter = 0; iter < vec.size(); ++iter) {
831 if (fDebug) printf(" vec.at(%zu): %8.4f | Si %d - Seg2 %d - Seg3 %d\n", iter, vec.at(iter).Chi2m, vec.at(iter).Ind1, vec.at(iter).Inds2, vec.at(iter).Inds3 );
832 isMatch = 0;
833 for (size_t InIter = 0; InIter < OutVec.size(); ++InIter) {
834 if ( vec.at(iter).Ind1 == OutVec.at(InIter).Ind1 && OutVec.at(InIter).Ind2 > -1 ) {
835 isMatch = 1;
836 continue;
837 }
838 if((vec.at(iter).Ind1 == OutVec.at(InIter).Ind1 && OutVec.at(InIter).Inds3 > -1) || (vec.at(iter).Ind1 == OutVec.at(InIter).Ind1 && vec.at(iter).Inds2 == OutVec.at(InIter).Inds2) ){
839 isMatch = 1;
840 continue;
841 }
842 if((vec.at(iter).Ind1 == OutVec.at(InIter).Ind1 && OutVec.at(InIter).Inds2 > -1) || (vec.at(iter).Ind1 == OutVec.at(InIter).Ind1 && vec.at(iter).Inds3 == OutVec.at(InIter).Inds3) ){
843 isMatch = 1;
844 continue;
845 }
846 if( (OutVec.at(InIter).Inds2 > -1 && vec.at(iter).Inds2 == OutVec.at(InIter).Inds2 ) || ( OutVec.at(InIter).Inds3 > -1 && vec.at(iter).Inds3 == OutVec.at(InIter).Inds3) ){
847 isMatch = 1;
848 continue;
849 }
850 }
851 if (fDebug) cout<<" writing unique index "<<endl;
852 if (isMatch == 0) {
853 if (fDebug) cout<<"--isMatch: Chi2 "<<vec.at(iter).Chi2m<<" Inds2 "<<vec.at(iter).Inds2<<" Inds3 "<<vec.at(iter).Inds3<<endl;
854 OutVec.push_back(vec.at(iter));
855 }
856 }//iter
857
858 if (fDebug) cout<<" RecordingBest: Now OutVec.size "<<OutVec.size()<<endl;
859 for (size_t iter = 0; iter < OutVec.size(); ++iter) {
860 if (fDebug) cout<<" iter "<<iter<<" index: Si:"<<OutVec.at(iter).Ind1<<" PCtr: "<< OutVec.at(iter).Ind2<<" PCseg2: "<<OutVec.at(iter).Inds2<<" PCseg3: "<<OutVec.at(iter).Inds3<<endl;
861 }
862
863 vec.clear();
864 }
865
866}//RecordingBest
867//----------------------------------------------------------------------
868
869
870//----------------------------------------------------------------------
871void BmnUpstreamTracking::GetAddSiHits(vector<Smatch> & Out, Double_t*** hits, Int_t* Nh){
872Bool_t was_ind;
873Double_t tmp_array[kNumStSi][6];
874if (fDebug) cout<<"--GetAddSiHits--"<<endl;
875
876 for (size_t iter = 0; iter < Out.size(); ++iter) {
877 if (fDebug) printf(" Out.at(%zu): %8.4f | Si %d - PC %d Seg2 %d - Seg3 %d\n", iter, Out.at(iter).Chi2m, Out.at(iter).Ind1, Out.at(iter).Ind2, Out.at(iter).Inds2, Out.at(iter).Inds3 );
878
879 was_ind = 0;
880 for (int st=1; st < kNumStSi; ++st){
881 for (int i=0; i < 5; ++i){
882 tmp_array[st][i] = -999.;
883 }
884 for (int ihit=0; ihit < Nh[st]; ++ihit){
885 if ( Out.at(iter).Ind1 == ihit ){
886 was_ind = 1;
887 tmp_array[st][0] = hits[st][ihit][0];
888 tmp_array[st][1] = hits[st][ihit][1];
889 tmp_array[st][2] = hits[st][ihit][2];
890 tmp_array[st][3] = hits[st][ihit][3];
891 tmp_array[st][4] = hits[st][ihit][4];
892 tmp_array[st][5] = hits[st][ihit][5];
893 }
894 }
895 }//st
896 if (was_ind){
897 for(int st = 1; st < kNumStSi; ++st){
898 Out.at(iter).SiCoordX[st] = tmp_array[st][0];
899 Out.at(iter).SiCoordY[st] = tmp_array[st][1];
900 Out.at(iter).SiCoordZ[st] = tmp_array[st][2];
901 Out.at(iter).SiSigmaX[st] = tmp_array[st][3];
902 Out.at(iter).SiSigmaY[st] = tmp_array[st][4];
903 Out.at(iter).SiMod[st] = tmp_array[st][5];
904 if (fDebug) cout<<" st "<<st<<" X "<<tmp_array[st][0]<<" Y "<<tmp_array[st][1]<<" Z "<<tmp_array[st][2]<<" mod "<<tmp_array[st][5]<<endl;
905 }
906 }
907 }//iter
908
909}//GetAddSiHits
910//----------------------------------------------------------------------
911
912
913//----------------------------------------------------------------------
914void BmnUpstreamTracking::GetHitsToPoints(vector<Smatch> & Out,Double_t*** ab_tr, Int_t* NPCTracks_, Double_t*** ab_Ch,
915Int_t* Nseg_ch_, Double_t*** Points_, Int_t & NumPoints_){
916 if (fDebug) cout<<" ---GetHitsToPoints---"<<endl;
917 for (size_t iter = 0; iter < Out.size(); iter++) {
918 if ( fDebug ) cout<<" silicon points "<<endl;
919 for (int ist=1; ist < kNumStSi; ++ist){
920 if (ist < 3){
921 Points_[iter][ist-1][0] = Out.at(iter).SiCoordX[ist];//x
922 Points_[iter][ist-1][1] = Out.at(iter).SiCoordY[ist];//y
923 Points_[iter][ist-1][2] = Out.at(iter).SiCoordZ[ist];//z
924 Points_[iter][ist-1][3] = Out.at(iter).SiSigmaX[ist];//sigx
925 Points_[iter][ist-1][4] = Out.at(iter).SiSigmaY[ist];//sigy
926 if ( fDebug ) cout<<" X["<<ist<<"] "<<Points_[iter][ist-1][0] <<" Y["<<ist<<"] "<<Points_[ist-1][ist][1] <<endl;
927 }else{
928 Points_[iter][ist][0] = Out.at(iter).SiCoordX[ist];//x
929 Points_[iter][ist][1] = Out.at(iter).SiCoordY[ist];//y
930 Points_[iter][ist][2] = Out.at(iter).SiCoordZ[ist];//z
931 Points_[iter][ist][3] = Out.at(iter).SiSigmaX[ist];//sigx
932 Points_[iter][ist][4] = Out.at(iter).SiSigmaY[ist];//sigy
933 if ( fDebug ) cout<<" X["<<ist<<"] "<<Points_[iter][ist][0] <<" Y["<<ist<<"] "<<Points_[iter][ist][1] <<endl;
934 }
935 }
936 for (int ipctr=0; ipctr < NPCTracks_[1]; ++ipctr){
937 // if was MWPC-track
938 if ( Out.at(iter).Ind2 == ipctr){
939 // x = kx*(dz) + bx //Z_pair[kNumPairs] = {-808.816, -289.401}; //Z_Chamber[kNumChambers] = {-858., -761., -364., -214.};
940 Double_t deltaAx = ab_tr[1][7][ipctr];
941 Double_t deltaAy = ab_tr[1][8][ipctr];
942
943 Points_[iter][2][0] = Out.at(iter).param2[0]*(Z_Chamber[2]-Z_pair[1]) + Out.at(iter).param2[1];//x3
944 Points_[iter][4][0] = Out.at(iter).param2[0]*(Z_Chamber[3]-Z_pair[1]) + Out.at(iter).param2[1];//x4 //(Z_Chamber[3]-Z_pair[1])
945
946 Points_[iter][2][1] = Out.at(iter).param2[2]*(Z_Chamber[2]-Z_pair[1]) + Out.at(iter).param2[3];//y3
947 Points_[iter][4][1] = Out.at(iter).param2[2]*(Z_Chamber[3]-Z_pair[1]) + Out.at(iter).param2[3];//y4 //Z_pair[1] - Z_Chamber[3]
948 // Ax*deltaX + X*deltaAx + deltaX
949 Points_[iter][2][3] = Out.at(iter).param2[0]*SigmXmwpc + Out.at(iter).param2[1]*deltaAx + SigmXmwpc;//sigx
950 Points_[iter][2][4] = Out.at(iter).param2[2]*SigmYmwpc + Out.at(iter).param2[3]*deltaAy + SigmYmwpc;//sigy
951 Points_[iter][4][3] = SigmXmwpc;//sigx// Points_[iter][3][3];
952 Points_[iter][4][4] = SigmYmwpc;//sigy// Points_[iter][3][4];
953 if ( fDebug ){
954 cout<<" MWPCs from track"<<endl;
955 cout<<" X3 "<<Points_[iter][2][0] <<" X4 "<<Points_[iter][4][0] <<endl;//" sigX "<<Points_[iter][3][3]<<endl;
956 cout<<" Y3 "<<Points_[iter][2][1] <<" X4 "<<Points_[iter][4][1] <<endl;//" sigX "<<Points_[iter][3][3]<<endl;
957 //" sigXcov "<<ab_tr[1][5][ipctr]<<endl;
958 }
959 }
960 }//NPCTracks_[1]
961 //--if was segment--
962 for (int iseg2=0; iseg2 < Nseg_ch_[2]; ++iseg2){
963 if ( Out.at(iter).Inds2 == iseg2){
964 Points_[iter][2][0] = ab_Ch[2][1][iseg2];//x
965 Points_[iter][2][1] = ab_Ch[2][3][iseg2];//y
966 Points_[iter][2][3] = SigmXmwpc;//sigx
967 Points_[iter][2][4] = SigmYmwpc;//sigy
968 if ( fDebug ){
969 cout<<" MWPCs from segment"<<endl;
970 cout<<" X3 "<<Points_[iter][2][0] <<" Y3 "<<Points_[iter][2][1]<<endl;//" sigX "<<Points_[iter][3][3]<<endl;
971 //" sigXcov "<<ab_tr[1][5][ipctr]<<endl;
972 }
973 }
974 }//Nseg_ch_[2]
975 for (int iseg3=0; iseg3 < Nseg_ch_[3]; ++iseg3){
976 if ( Out.at(iter).Inds3 == iseg3){
977 Points_[iter][4][0] = ab_Ch[3][1][iseg3];//x
978 Points_[iter][4][1] = ab_Ch[3][3][iseg3];//y
979 Points_[iter][4][3] = SigmXmwpc;//sigx
980 Points_[iter][4][4] = SigmYmwpc;//sigy
981 if ( fDebug ){
982 cout<<" MWPCs from segment"<<endl;
983 cout<<" X4 "<<Points_[iter][4][0] <<" Y4 "<<Points_[iter][4][1]<<endl;
984 }
985 }
986 }//Nseg_ch_[2]
987
988
989 Points_[iter][2][2] = Z_Chamber[2];
990 Points_[iter][4][2] = Z_Chamber[3];
991 NumPoints_ = iter+1;
992 }//Out.size()
993 //????????
994}
995//----------------------------------------------------------------------
996
997
998//----------------------------------------------------------------------
999void BmnUpstreamTracking::CalculationOfChi2(Double_t*** Points_, Int_t & NumPoints_, vector<UpTracks> & vec ){
1000 if (fDebug) cout<<" Check points before fit "<<" Zcentr "<<Zcentr<<endl;
1001 Int_t Nhits_before;
1002 for (int ip =0; ip < NumPoints_; ++ip){
1003 Nhits_before = 0;
1004 for (int ist =0; ist < kPoints; ++ist){
1005 Points_[ip][ist][0] += X_shiftUp[ist];
1006 Points_[ip][ist][1] += Y_shiftUp[ist];
1007 Points_[ip][ist][2] -= Zcentr;
1008 // if (fDebug) cout<<"ip "<<ip<<" X["<<ist<<"] "<< Points_[ip][ist][0]<<" Y "<<Points_[ip][ist][1]<<" Z "<<Points_[ip][ist][2]<<" sigx "<<Points_[ip][ist][3]<<" sigy "<<Points_[ip][ist][4]<<endl;
1009 if (Points_[ip][ist][0] > -900.) Nhits_before++;
1010 if (Points_[ip][ist][1] > -900.) Nhits_before++;
1011 }
1012 }
1013 if (fDebug) cout<<" Nhits_before "<<Nhits_before<<endl;
1014 UpTracks tmptr;
1015
1016 for (Int_t i = 0; i < kPoints; ++i) {
1017 for (Int_t j = 0; j < 5; ++j) {
1018 Xforglfit[i][j] = -999.;
1019 }
1020 }
1021 Double_t ChiSquare_ ;
1022 for (Int_t itr = 0; itr < NumPoints_; ++itr) {
1023
1024 if (Nhits_before < 6 ) continue;
1025
1026 ChiSquare_ = 0.;
1027 if(fDebug) cout<<endl;
1028 if(fDebug) cout<<"--GlobalFit--"<<endl;
1029 //setting global fit matrices to zero
1030
1031 for (Int_t i = 0; i < 4; ++i){
1032 for (Int_t j = 0; j < 4; ++j){
1033 Amatr[i][j] = 0.; }}
1034 for (Int_t i = 0; i < 4; ++i) { rhs[i] = 0; }
1035 for (Int_t i = 0; i < 16; ++i) { AmatrArray[i] = 0; }
1036 for (Int_t i = 0; i < 16; ++i) { AmatrInverted[i] = 0; }
1037
1038 for (Int_t i = 0; i < kPoints; ++i) {
1039 for (Int_t j = 0; j < kPoin_Par; ++j) {
1040 Xforglfit[i][j] = Points_[itr][i][j];
1041 //if ( fDebug && j == 0) cout<<" x ["<<i<<"]["<<j<<"] "<<Xforglfit[i][j]<<endl;
1042 }
1043 for (Int_t j = 0; j < kPoin_Par; ++j) {
1044 // if ( fDebug && j == 1) cout<<" y ["<<i<<"]["<<j<<"] "<<Xforglfit[i][j]<<endl;
1045 }
1046 }
1047
1048 GlobalFit(Xforglfit,Amatr, rhs);
1049
1050 //here we do smth with Amatr & rhs (solving)
1051 // Amatr is row-major
1052 for (Int_t row = 0; row < 4; ++row) {
1053 for (Int_t col = 0; col < 4; ++col) {
1054 AmatrArray[4*col + row] = Amatr[row][col];
1055 }
1056 }//AmatrArray is col-major
1057
1058 bool ifSuccess = InvertMatrix(AmatrArray, AmatrInverted);
1059
1060 if (!ifSuccess) {cout<< "???? InvertMatrix is unsuccessfull ????"<<endl; }
1061 if (!ifSuccess) {cout<< "???? InvertMatrix is unsuccessfull ????"<<endl; return;}
1062 //line_ = { ax, bx, ay, by}
1063 Double_t line_[4] = {0.,0.,0.,0.};
1064 for (Int_t coeff = 0; coeff < 4; ++coeff) {
1065 for (Int_t iCol = 0; iCol < 4; ++iCol) {
1066 // if ( fDebug ) cout<<" AmatrInverted["<<4*iCol + coeff<<"] "<<AmatrInverted[4*iCol + coeff]<<" rhs["<<iCol<<"] "<<rhs[iCol]<<endl;
1067 line_[coeff] += AmatrInverted[4*iCol + coeff]*rhs[iCol];
1068 }
1069 }
1070
1071 if ( fDebug ) cout<<" ax "<<line_[0]<<" bx "<<line_[1]<<" ay "<<line_[2]<<" by "<<line_[3]<<endl;
1072 if ( fDebug ) cout<<endl;
1073
1074 //---------------------------Chi2---------------------------------
1075 Int_t nhits = 0;
1076 for (Int_t ist = 0; ist < kPoints; ++ist) {
1077 if (Xforglfit[ist][0] < -900. ) continue;//x
1078
1079 Double_t hit = Xforglfit[ist][0] ;
1080 Double_t fit = line_[0]*Xforglfit[ist][2] + line_[1];
1081 nhits++;
1082
1083 Double_t sigmaX_2 = Xforglfit[ist][3]*Xforglfit[ist][3];
1084 ChiSquare_ += ((hit-fit)*(hit-fit))/(sigmaX_2);
1085 //if ( fDebug )
1086 // if ((ist == 2 || ist == 4) && fabs(hit-fit) > 0.2 ) cout<<" fEventNo "<<fEventNo<<" PC ist "<<ist<<" hitx "<<hit<<" fit "<<fit<<" dx "<<hit-fit<<" Z "<<Xforglfit[ist][2]<<" sigma "<<Xforglfit[ist][3]<<" ChiSquare_ "<<ChiSquare_<<endl;
1087 // if ((ist == 0 || ist == 1 || ist == 3) && fabs(hit-fit) > 0.02 ) cout<<" fEventNo "<<fEventNo<<" Si ist "<<ist<<" hitx "<<hit<<" fit "<<fit<<" dx "<<hit-fit<<" Z "<<Xforglfit[ist][2]<<" sigma "<<Xforglfit[ist][3]<<" ChiSquare_ "<<ChiSquare_<<endl;
1088 if ( fDebug ) cout<<" ist "<<ist<<" hitx "<<hit<<" fit "<<fit<<" dx "<<hit-fit<<" Z "<<Xforglfit[ist][2]<<" sigma "<<Xforglfit[ist][3]<<" ChiSquare_ "<<ChiSquare_<<endl;
1089
1090 }
1091 for (Int_t ist = 0; ist < kPoints; ++ist) {
1092 if (Xforglfit[ist][1] < -900. ) continue;//y
1093
1094 Double_t hit = Xforglfit[ist][1];
1095 Double_t fit = line_[2]*Xforglfit[ist][2] + line_[3];
1096 nhits++;
1097
1098 Double_t sigmaY_2 = Xforglfit[ist][4]*Xforglfit[ist][4];
1099 ChiSquare_ += ((hit-fit)*(hit-fit))/(sigmaY_2);
1100 //if ( fDebug)
1101 //if (fabs(hit-fit) > 0.2 )
1102 if ( fDebug) cout<<" ist "<<ist<<" hity "<<hit<<" fit "<<fit<<" dy "<<hit-fit<<" Z "<<Xforglfit[ist][2]<<" sigma "<<Xforglfit[ist][4]<<" ChiSquare_ "<<ChiSquare_<<endl;
1103 }//ist
1104
1105 Double_t ChiSquareNdf = ChiSquare_/(nhits - 4);
1106
1107 if (ChiSquareNdf < kChi2_Max && nhits > 5){
1108
1109 tmptr.Chi2 = ChiSquareNdf;
1110
1111 if ( fDebug) cout<<" Chi2 "<<tmptr.Chi2<<" nhits "<<nhits<<endl;
1112 if ( fDebug) cout<<" ------"<<endl;
1113
1114 for (Int_t iCol = 0; iCol < 4; ++iCol) {
1115 tmptr.param[iCol] = line_[iCol];
1116 }
1117 // if ( fDebug) cout<<" Ax "<<tmptr.param[0]<<" bx "<<tmptr.param[1]<<" Ay "<<tmptr.param[2]<<" by "<<tmptr.param[3]<<endl;
1118 // if ( fDebug) cout<<" X coor "<<endl;
1119 for (Int_t ist = 0; ist < kPoints; ++ist) {
1120 tmptr.CoordZ[ist] = Xforglfit[ist][2];
1121 tmptr.CoordX[ist] = Xforglfit[ist][0];
1122 if (tmptr.CoordX[ist]< -900. ) continue;
1123 // if ( fDebug) cout<<" ist "<<ist<<" hitx "<<tmptr.CoordX[ist]<<endl;
1124 }
1125 //if ( fDebug) cout<<" Y coor "<<endl;
1126 for (Int_t ist = 0; ist < kPoints; ++ist) {
1127 tmptr.CoordY[ist] = Xforglfit[ist][1];
1128 if (tmptr.CoordY[ist]< -900. ) continue;
1129 // if ( fDebug) cout<<" ist "<<ist<<" hity "<<tmptr.CoordY[ist]<<endl;
1130 }
1131 tmptr.Nhits = nhits;
1132 vec.push_back(tmptr);
1133 }
1134 }//itr
1135
1136
1137}//calculationOfChi2
1138//----------------------------------------------------------------------
1139
1140//----------------------------------------------------------------------
1141void BmnUpstreamTracking::PrintAllTracks(vector<UpTracks> & vecUp){
1142
1143 cout<<" ------PrintAllTracks: "<<vecUp.size()<<endl;
1144 for (size_t itr =0; itr < vecUp.size(); ++itr){
1145 cout<<endl;
1146 cout<<" itr "<<itr<<" Chi2 "<<vecUp.at(itr).Chi2<<" Ax "<<vecUp.at(itr).param[0]<<" bx "<<vecUp.at(itr).param[1]<<" Ay "<<vecUp.at(itr).param[2]<<" by "<<vecUp.at(itr).param[3]<<" Nhits "<<vecUp.at(itr).Nhits<<endl;
1147
1148 cout<<"--X--"<<endl;
1149 for (Int_t st = 0; st < kPoints; ++st) {
1150 if ( vecUp.at(itr).CoordX[st] < -900.) continue;
1151 Double_t Hit = vecUp.at(itr).CoordX[st];
1152 Double_t Fit = vecUp.at(itr).param[0] * vecUp.at(itr).CoordZ[st] + vecUp.at(itr).param[1];
1153 hResXst[st]->Fill(Hit - Fit);
1154 cout<<" st "<<st<<" Xhit "<<Hit<<" Xfit "<< Fit<<" Z "<<vecUp.at(itr).CoordZ[st]<<endl;
1155 }
1156 cout<<"--Y--"<<endl;
1157 for (Int_t st = 0; st < kPoints; ++st) {
1158 if ( vecUp.at(itr).CoordY[st] < -900.) continue;
1159 Double_t Hit = vecUp.at(itr).CoordY[st];
1160 Double_t Fit = vecUp.at(itr).param[2] * vecUp.at(itr).CoordZ[st] + vecUp.at(itr).param[3];
1161 hResYst[st]->Fill(Hit - Fit);
1162 cout<<" st "<<st<<" Yhit "<<Hit<<" Yfit "<< Fit<<" Z "<<vecUp.at(itr).CoordZ[st]<<endl;
1163 }
1164 Double_t x_target = vecUp.at(itr).param[0]*( kZ_target - Zcentr ) + vecUp.at(itr).param[1];
1165 Double_t y_target = vecUp.at(itr).param[2]*( kZ_target - Zcentr ) + vecUp.at(itr).param[3];
1166
1167 hAx_fitUp->Fill(vecUp.at(itr).param[0]);
1168 hAy_fitUp->Fill(vecUp.at(itr).param[2]);
1169 hx_fitUp->Fill(vecUp.at(itr).param[1]);
1170 hy_fitUp->Fill(vecUp.at(itr).param[3]);
1171 hvertexXYUp->Fill(x_target,y_target);
1172 hy_vs_x_Up->Fill(vecUp.at(itr).param[1],vecUp.at(itr).param[3]);
1173 hchi2_fitUp->Fill(vecUp.at(itr).Chi2);
1174 hNhitsUp->Fill(vecUp.at(itr).Nhits);
1175 hTyTx_Up->Fill(vecUp.at(itr).param[0],vecUp.at(itr).param[2]);
1176
1177 }
1178 cout<<"---"<<endl;
1179}
1180//----------------------------------------------------------------------
1181
1182//----------------------------------------------------------------------
1183void BmnUpstreamTracking::TrackRecording(vector<UpTracks> & vecUp)
1184{
1185 for (size_t InIter = 0; InIter < vecUp.size(); ++InIter) {
1186/*
1187 Double_t z0 = vecUp.at(InIter).CoordZ[1] + Zcentr;
1188 Double_t x0 = vecUp.at(InIter).param[0] * vecUp.at(InIter).CoordZ[0] + vecUp.at(InIter).param[1]; //ax(z)+bx
1189 Double_t y0 = vecUp.at(InIter).param[2] * vecUp.at(InIter).CoordZ[0] + vecUp.at(InIter).param[3]; //y
1190 Double_t z4 = vecUp.at(InIter).CoordZ[4] + Zcentr;
1191 Double_t x4 = vecUp.at(InIter).param[0] * vecUp.at(InIter).CoordZ[4] + vecUp.at(InIter).param[1]; //ax(z)+bx
1192 Double_t y4 = vecUp.at(InIter).param[2] * vecUp.at(InIter).CoordZ[4] + vecUp.at(InIter).param[3]; //y
1193
1194 FairTrackParam trackParamf;
1195 trackParamf.SetPosition(TVector3(x0, y0, z0));
1196 trackParamf.SetTx(vecUp.at(InIter).param[0]);
1197 trackParamf.SetTy(vecUp.at(InIter).param[2]);
1198
1199 FairTrackParam trackParaml;
1200 trackParaml.SetPosition(TVector3(x4, y4, z4));
1201 trackParaml.SetTx(vecUp.at(InIter).param[0]);
1202 trackParaml.SetTy(vecUp.at(InIter).param[2]);
1203*/
1204
1205 FairTrackParam trackParaml;
1206 trackParaml.SetPosition(TVector3(vecUp.at(InIter).param[1]+ Shift_toCenterOfMagnetX, vecUp.at(InIter).param[3]+ Shift_toCenterOfMagnetY,Zcentr));
1207 trackParaml.SetTx(vecUp.at(InIter).param[0]+Shift_toCenterOfMagnetAX);
1208 trackParaml.SetTy(vecUp.at(InIter).param[2]+Shift_toCenterOfMagnetAY);
1209
1210 BmnTrack *track = new ((*fBmnUpstreamTracksArray)[fBmnUpstreamTracksArray->GetEntriesFast()]) BmnTrack();
1211 track->SetChi2(vecUp.at(InIter).Chi2);
1212 track->SetNHits(vecUp.at(InIter).Nhits);
1213 track->SetFlag(InIter);
1214 //track->SetParamFirst(trackParamf);
1215 track->SetParamLast(trackParaml);
1216
1217 }//InIter
1218
1219}
1220//----------------------------------------------------------------------
1221
1222//----------------------------------------------------------------------
1223bool BmnUpstreamTracking::InvertMatrix(Double_t *m, Double_t *invOut) {
1224 Double_t inv[16], det;
1225 int i;
1226
1227 inv[0] = m[5] * m[10] * m[15] -m[5] * m[11] * m[14] -m[9] * m[6] * m[15] +m[9] * m[7] * m[14] +m[13] * m[6] * m[11] -m[13] * m[7] * m[10];
1228 inv[4] = -m[4] * m[10] * m[15] +m[4] * m[11] * m[14] +m[8] * m[6] * m[15] -m[8] * m[7] * m[14] -m[12] * m[6] * m[11] +m[12] * m[7] * m[10];
1229 inv[8] = m[4] * m[9] * m[15] -m[4] * m[11] * m[13] -m[8] * m[5] * m[15] +m[8] * m[7] * m[13] +m[12] * m[5] * m[11] -m[12] * m[7] * m[9];
1230 inv[12] = -m[4] * m[9] * m[14] +m[4] * m[10] * m[13] +m[8] * m[5] * m[14] -m[8] * m[6] * m[13] -m[12] * m[5] * m[10] +m[12] * m[6] * m[9];
1231 inv[1] = -m[1] * m[10] * m[15] +m[1] * m[11] * m[14] +m[9] * m[2] * m[15] -m[9] * m[3] * m[14] -m[13] * m[2] * m[11] +m[13] * m[3] * m[10];
1232 inv[5] = m[0] * m[10] * m[15] -m[0] * m[11] * m[14] -m[8] * m[2] * m[15] +m[8] * m[3] * m[14] +m[12] * m[2] * m[11] -m[12] * m[3] * m[10];
1233 inv[9] = -m[0] * m[9] * m[15] +m[0] * m[11] * m[13] +m[8] * m[1] * m[15] -m[8] * m[3] * m[13] -m[12] * m[1] * m[11] +m[12] * m[3] * m[9];
1234 inv[13] = m[0] * m[9] * m[14] -m[0] * m[10] * m[13] -m[8] * m[1] * m[14] +m[8] * m[2] * m[13] +m[12] * m[1] * m[10] -m[12] * m[2] * m[9];
1235 inv[2] = m[1] * m[6] * m[15] -m[1] * m[7] * m[14] -m[5] * m[2] * m[15] + m[5] * m[3] * m[14] +m[13] * m[2] * m[7] -m[13] * m[3] * m[6];
1236 inv[6] = -m[0] * m[6] * m[15] +m[0] * m[7] * m[14] +m[4] * m[2] * m[15] -m[4] * m[3] * m[14] -m[12] * m[2] * m[7] +m[12] * m[3] * m[6];
1237 inv[10] = m[0] * m[5] * m[15] -m[0] * m[7] * m[13] -m[4] * m[1] * m[15] +m[4] * m[3] * m[13] + m[12] * m[1] * m[7] - m[12] * m[3] * m[5];
1238 inv[14] = -m[0] * m[5] * m[14] +m[0] * m[6] * m[13] +m[4] * m[1] * m[14] -m[4] * m[2] * m[13] -m[12] * m[1] * m[6] +m[12] * m[2] * m[5];
1239 inv[3] = -m[1] * m[6] * m[11] +m[1] * m[7] * m[10] + m[5] * m[2] * m[11] -m[5] * m[3] * m[10] - m[9] * m[2] * m[7] +m[9] * m[3] * m[6];
1240 inv[7] = m[0] * m[6] * m[11] -m[0] * m[7] * m[10] -m[4] * m[2] * m[11] +m[4] * m[3] * m[10] +m[8] * m[2] * m[7] -m[8] * m[3] * m[6];
1241 inv[11] = -m[0] * m[5] * m[11] +m[0] * m[7] * m[9] +m[4] * m[1] * m[11] -m[4] * m[3] * m[9] -m[8] * m[1] * m[7] +m[8] * m[3] * m[5];
1242 inv[15] = m[0] * m[5] * m[10] -m[0] * m[6] * m[9] -m[4] * m[1] * m[10] +m[4] * m[2] * m[9] +m[8] * m[1] * m[6] -m[8] * m[2] * m[5];
1243 det = m[0] * inv[0] + m[1] * inv[4] + m[2] * inv[8] + m[3] * inv[12];
1244
1245 if (det == 0)
1246 return false;
1247 det = 1.0 / det;
1248 for (i = 0; i < 16; i++)
1249 invOut[i] = inv[i] * det;
1250 return true;
1251}
1252//----------------------------------------------------------------------
1253
1254//----------------------------------------------------------------------
1255void BmnUpstreamTracking::GlobalFit(Double_t** XYZsigXY, Double_t** Amatr_, Double_t* rhs_){
1256
1257 Double_t sigmaX_2 = 999.;
1258 Double_t sigmaY_2 = 999.;
1259 Double_t X = 999., Y = 999., Z =999.;
1260
1261 for (Int_t iSt = 0; iSt < kPoints; ++iSt) {
1262
1263 Z = XYZsigXY[iSt][2];
1264
1265 if (XYZsigXY[iSt][0] > -900.){
1266 X = XYZsigXY[iSt][0];
1267 sigmaX_2 = XYZsigXY[iSt][3]* XYZsigXY[iSt][3];
1268
1269 Amatr_[0][0] += 2*Z*Z/sigmaX_2;
1270 Amatr_[0][1] += 2*Z/sigmaX_2;
1271 rhs_[0] += 2*Z*X/sigmaX_2;
1272
1273 Amatr_[1][0] += 2*Z/sigmaX_2;
1274 Amatr_[1][1] += 2/sigmaX_2;
1275 rhs_[1] += 2*X/sigmaX_2;
1276 }
1277
1278 if (XYZsigXY[iSt][1] > -900.){
1279 Y = XYZsigXY[iSt][1];
1280 sigmaY_2 = XYZsigXY[iSt][4]* XYZsigXY[iSt][4];
1281
1282 Amatr_[2][2] += 2*Z*Z/sigmaY_2;
1283 Amatr_[2][3] += 2*Z/sigmaY_2;
1284 rhs_[2] += 2*Z*Y/sigmaY_2;
1285
1286 Amatr_[3][2] += 2*Z/sigmaY_2;
1287 Amatr_[3][3] += 2/sigmaY_2;
1288 rhs_[3] += 2*Y/sigmaY_2;
1289 }
1290
1291 for (Int_t coeff = 0; coeff < 4; ++coeff) {
1292 for (Int_t iCol = 0; iCol < 4; ++iCol) {
1293 //if ( fDebug ) cout << "GLOBALFIT: station = " << iSt << " " << "Amatr_[" << coeff << "][" << iCol << "] = " << Amatr_[coeff][iCol] << "\n";
1294 }
1295 // if ( fDebug ) cout<<" rhs_["<<coeff<<"] "<< rhs_[coeff]<<endl;
1296 }
1297 }
1298
1299}//GlobalFit
1300//----------------------------------------------------------------------
1301
1302
1303//----------------------------------------------------------------------
1304InitStatus BmnUpstreamTracking::Init() {
1305
1306 if (fVerbose) cout << "BmnUpstreamTracking == Init ==" << endl;
1307 FairRootManager* ioman = FairRootManager::Instance();
1308
1309 //----Data & MC---
1310 fSiTracks = (TClonesArray*) ioman->GetObject(fInputBranchName1);
1311 if (!fSiTracks)
1312 {
1313 cout<<"BmnUpstreamTracking::Init(): branch "<<fInputBranchName1<<" not found! Task will be deactivated"<<endl;
1314 SetActive(kFALSE);
1315 return kERROR;
1316 }
1317
1318 fSilHits = (TClonesArray*) ioman->GetObject(fInputBranchHits);
1319 if (!fSilHits){
1320 cout<<"BmnUpstreamTracking::Init(): branch "<<fInputBranchHits<<" not found! Task will be deactivated"<<endl;
1321 SetActive(kFALSE);
1322 return kERROR;
1323 }
1324
1325 fMWPCTracks = (TClonesArray*) ioman->GetObject(fInputBranchName2);
1326 if (!fMWPCTracks){
1327 cout<<"BmnUpstreamTracking::Init(): branch "<<fInputBranchName2<<" not found! Task will be deactivated"<<endl;
1328 SetActive(kFALSE);
1329 return kERROR;
1330 }
1331 fMWPCSegments = (TClonesArray*) ioman->GetObject(fInputBranchName3);
1332 if (!fMWPCSegments)
1333 {
1334 cout<<"BmnUpstreamTracking::Init(): branch "<<fInputBranchName3<<" not found! Task will be deactivated"<<endl;
1335 SetActive(kFALSE);
1336 return kERROR;
1337 }
1338
1339 //----MC true---
1340 if (!expData ){
1341 cout<<" !expData "<<endl;
1342 fSiTracksSim = (TClonesArray*) ioman->GetObject(fInputBranchNameSimTrue);
1343 if (!fSiTracksSim){
1344 cout<<"BmnUpstreamTracking::Init(): branch "<<fInputBranchNameSimTrue<<" not found! Task will be deactivated"<<endl;
1345 SetActive(kFALSE);
1346 return kERROR;
1347 }
1348 fBmnHitsArray = (TClonesArray*)ioman->GetObject(fInputBranchName);
1349 cout << "fBmnHitsArray = " << fInputBranchName << "\n";
1350 if (!fBmnHitsArray) {
1351 cout << "BmnMwpcHitFinder::Init(): branch " << fInputBranchName << " not found! Task will be deactivated" << endl;
1352 SetActive(kFALSE);
1353 return kERROR;
1354 }
1355 }// !expData
1356
1357 // Create and register tracks arrays
1358 fBmnUpstreamTracksArray = new TClonesArray(fOutputTracksBranchName);
1359 ioman->Register("BmnUpstreamTrack", "UpstreamTrack", fBmnUpstreamTracksArray, kTRUE);
1360
1361 par_ab_Ch = new Double_t**[kNumChambers];
1362 par_ab_tr = new Double_t**[kNumPairs];
1363 SiXYhits = new Double_t**[kNumStSi];
1364 Points = new Double_t**[kBig];
1365 par_ab_SiTr = new Double_t*[kNumPars];
1366 par_ab_trz = new Double_t*[kNumPairs];
1367 par_Seg_z = new Double_t*[kNumChambers];
1368 Xforglfit = new Double_t*[kPoints];
1369 Amatr = new Double_t*[4];
1370 par_SiTr_z = new Double_t[kBig];
1371 X_shift_seg = new Double_t[kNumChambers];
1372 Z_pair = new Double_t[kNumPairs];
1373 Z_Chamber = new Double_t[kNumChambers];
1374 AmatrInverted = new Double_t[16];
1375 rhs = new Double_t[4];
1376 AmatrArray = new Double_t[16];
1377 line = new Double_t[4];
1378 Nseg_Ch = new Int_t[kNumChambers];
1379 NPCTracks = new Int_t[kNumPairs];
1380 NSiXYhits = new Int_t[kNumStSi];
1381 X_shiftUp = new Double_t[kPoints];
1382 Y_shiftUp = new Double_t[kPoints];
1383
1384 for (Int_t i = 0; i < 4; i++) {
1385 Amatr[i] = new Double_t[4];
1386 }
1387
1388 for (Int_t iPars = 0; iPars < kNumPars; iPars++) {
1389 par_ab_SiTr[iPars] = new Double_t[kBig];
1390 }
1391
1392 for (Int_t iCh = 0; iCh < kNumChambers; iCh++) {
1393 par_Seg_z[iCh] = new Double_t[kBig];
1394 par_ab_Ch[iCh] = new Double_t*[kNumPars];
1395 for (Int_t iPars = 0; iPars < kNumPars; iPars++) {
1396 par_ab_Ch[iCh][iPars] = new Double_t[kBig];
1397 }
1398 }
1399 for (Int_t ip = 0; ip < kNumPairs; ip++) {
1400 par_ab_tr[ip] = new Double_t*[kNumPars];
1401 par_ab_trz[ip] = new Double_t[kBig];
1402 for (Int_t iPars = 0; iPars < kNumPars; iPars++) {
1403 par_ab_tr[ip][iPars] = new Double_t[kBig];
1404 }
1405 }
1406 for (Int_t ist = 0; ist < kNumStSi; ist++) {
1407 SiXYhits[ist] = new Double_t*[kBig];
1408 for (Int_t ib = 0; ib < kBig; ib++) {
1409 SiXYhits[ist][ib] = new Double_t[kNumPars];
1410 }
1411 }
1412 for (Int_t ib = 0; ib < kBig; ib++) {
1413 Points[ib] = new Double_t*[kPoints];
1414 for (Int_t ip = 0; ip < kPoints; ip++) {
1415 Points[ib][ip] = new Double_t[kPoin_Par];
1416 }
1417 }
1418 for (Int_t ip = 0; ip < kPoints; ip++) {
1419 Xforglfit[ip] = new Double_t[kPoin_Par];
1420 }
1421
1422 //--some hists--
1423 if (fDebug){
1424
1425 hAx_fitUp = new TH1D("Ax_fitUp","Tx_fitUp;[rad]", 200, -0.03, 0.03);fList.Add(hAx_fitUp);
1426 hAy_fitUp = new TH1D("Ay_fitUp","Ty_fitUp;[rad]", 200, -0.05, 0.05); fList.Add(hAy_fitUp);
1427 hx_fitUp = new TH1D("x_fitUp","x_fitUp;[cm]", 200, -10, 10);fList.Add(hx_fitUp);
1428 hy_fitUp = new TH1D("y_fitUp","y_fitUp;[cm]", 200, -15, 5);fList.Add(hy_fitUp);
1429 hy_vs_x_Up = new TH2D("y_vs_x_Up",";X[cm];Y[cm]", 200, -10, 10, 200, -15, 5);fList.Add(hy_vs_x_Up);
1430 hvertexXYUp = new TH2D("vertexXYUp","vertexUp;X[cm];Y[cm]", 200, -10, 10, 200, -10, 10);fList.Add(hvertexXYUp);
1431 hTyTx_Up= new TH2D("hTyTx_Up","hTyTx_Up", 100, -0.04, 0.04, 100, -0.04, 0.04);fList.Add(hTyTx_Up);
1432
1433 hchi2_fitUp = new TH1D("chi2_fitUp","Chi2", 100, 0, 100);fList.Add(hchi2_fitUp);
1434 hNhitsUp = new TH1D("NhitsUp","Number of hits", 10, 0.5, 10.5);fList.Add(hNhitsUp);
1435
1436 hAx_upmc = new TH1D("Ax_upmc","Ax_upmc;[rad]", 100, -0.05, 0.05);fList.Add(hAx_upmc);
1437 hAy_upmc = new TH1D("Ay_upmc","Ay_upmc;[rad]", 100, -0.05, 0.05);fList.Add(hAy_upmc);
1438 hX_upmc = new TH1D("X_upmc","X_upmc;[rad]", 100, -10, 10);fList.Add(hX_upmc);
1439 hY_upmc = new TH1D("Y_upmc","Y_upmc;[rad]", 100, -10, 10);fList.Add(hY_upmc);
1440
1441 hdAx_tr_mc_comb = new TH1D("dAx_tr_mc_comb","dAx_tr_mc_comb", 100, -0.005, 0.005);fList.Add(hdAx_tr_mc_comb);
1442 hdAy_tr_mc_comb = new TH1D("dAy_tr_mc_comb","dAy_tr_mc_comb", 100, -0.005, 0.005);fList.Add(hdAy_tr_mc_comb);
1443 hdX_tr_mc_comb = new TH1D("dX_tr_mc_comb","dX_tr_mc_comb", 100, -10, 10); fList.Add(hdX_tr_mc_comb);
1444 hdY_tr_mc_comb = new TH1D("dY_tr_mc_comb","dY_tr_mc_comb", 100, -10, 10); fList.Add(hdY_tr_mc_comb);
1445
1446 hdAx_uptr_mc = new TH1D("dAx_uptr_mc","dAx_uptr_mc", 200, -0.002, 0.002);fList.Add(hdAx_uptr_mc);
1447 hdAy_uptr_mc = new TH1D("dAy_uptr_mc","dAy_uptr_mc", 200, -0.005, 0.005);fList.Add(hdAy_uptr_mc);
1448 hdX_uptr_mc = new TH1D("dX_uptr_mc","dX_uptr_mc", 200, -2, 2); fList.Add(hdX_uptr_mc);
1449 hdY_uptr_mc = new TH1D("dY_uptr_mc","dY_uptr_mc", 200, -2, 2); fList.Add(hdY_uptr_mc);
1450 hY_vs_Xmctrue= new TH2D("Y_vs_Xmctrue","mctrue;X[cm];Y[cm]", 200, -10, 10, 200, -15, 5);fList.Add(hY_vs_Xmctrue);
1451
1452
1453 hNtr_reco = new TH1D("Ntr_reco","N up tr reco", 10, 0,10); fList.Add(hNtr_reco);
1454 hNtr_mc = new TH1D("Ntr_mc","N up tr mc", 10, 0,10); fList.Add(hNtr_mc);
1455 hNtr_mc_vs_reco = new TH2D("Ntr_mc_vs_reco","; Nupreco; Nupmctrue", 10, 0,10, 10, 0,10);fList.Add(hNtr_mc_vs_reco);
1456 hNrecoTrif2mc = new TH1D("NrecoTrif2mc","N Up tracks if Nmc=2", 3, 0,3); fList.Add(hNrecoTrif2mc);
1457 hAngle_reco = new TH1D("Angle_reco","Angle_reco if Nreco = 2 ;[deg];N/0.02^{o}", 200, 0., 4.); fList.Add(hAngle_reco);
1458 hAngle_recoifNmc2 = new TH1D("Angle_recoifNmc2","Angle_reco if Nreco = 2 & Nmc = 2;[deg];N/0.02^{o}", 200, 0., 4.); fList.Add(hAngle_recoifNmc2);
1459 hAngle_recoifNmc2Cases = new TH1D("Angle_recoifNmc2Cases","Angle_reco if Nreco = 2 & Nmc = 2 & SelectedCases;[deg];N/0.02^{o}", 200, 0., 4.); fList.Add(hAngle_recoifNmc2Cases);
1460
1461 hNum_mcuptr = new TH1D("hNum_mcuptr","Num_mcuptr", 2, 0, 2);fList.Add(hNum_mcuptr);
1462 hDen_mcuptr = new TH1D("hDen_mcuptr","Den_mcuptr", 2, 0, 2);fList.Add(hDen_mcuptr);
1463 hEff_mcuptr = new TH1D("hEff_mcuptr","Eff_mc_uptr", 2, 0, 2);fList.Add(hEff_mcuptr);
1464 hEff_mcuptr->GetXaxis()->SetBinLabel(1,"2+3 coord");
1465 hEff_mcuptr->GetXaxis()->SetBinLabel(2,"only 3coord");
1466
1467 TH1D *h;
1468 for (Int_t i = 0; i < kPoints; ++i) {
1469 if (i == 2 || i == 4) h = new TH1D(Form("ResXst%d", i), Form("ResXst%d", i), 100, -0.5, 0.5);
1470 else h = new TH1D(Form("ResXst%d", i), Form("ResXst%d", i), 150, -0.015, 0.015);
1471 fList.Add(h);
1472 hResXst.push_back(h);
1473 }
1474 for (Int_t i = 0; i < kPoints; ++i) {
1475 if (i == 2 || i == 4) h = new TH1D(Form("ResYst%d", i), Form("ResYst%d", i), 100, -0.25, 0.25);
1476 //else if(i == 4) h = new TH1D(Form("ResYst%d", i), Form("ResYst%d", i), 100, -0.2, 0.2);
1477 else h = new TH1D(Form("ResYst%d", i), Form("ResYst%d", i), 100, -0.4, 0.4);
1478 fList.Add(h);
1479 hResYst.push_back(h);
1480 }
1481 }
1482
1483 Z_pair[0] = -808.816;
1484 Z_pair[1] = -289.401;
1485 Z_Chamber[0] = -858.155579;
1486 Z_Chamber[1] = -761.476562;
1487 Z_Chamber[2] = -364.255554;
1488 Z_Chamber[3] = -214.545547;
1489
1490 //--Shift--
1491 if (expData){
1492 X_shift_seg[0] = 0.;
1493 X_shift_seg[1] = 0.;
1494 X_shift_seg[2] = -0.04;
1495 X_shift_seg[3] = -0.05;
1496 Shift_toCenterOfMagnetX = 0.039;
1497 Shift_toCenterOfMagnetY = 0.9;//0.13;// drawing
1498 Shift_toCenterOfMagnetAX = 0.;
1499 Shift_toCenterOfMagnetAY = 0.0019;//
1500 X_shift = 0.02;
1501 Y_shift = 0.3;
1502
1503 }else{
1504 X_shift_seg[0] = 0.;
1505 X_shift_seg[1] = 0.;
1506 X_shift_seg[2] = 0.;
1507 X_shift_seg[3] = 0.;
1508 Shift_toCenterOfMagnetX = 0.;
1509 Shift_toCenterOfMagnetY = 0.;
1510 Shift_toCenterOfMagnetAX = 0.;
1511 Shift_toCenterOfMagnetAY = 0.;
1512 X_shift = 0.;
1513 Y_shift = 0.;
1514 }
1515
1516 //-- 5 detectors Shift---
1517 X_shiftUp[0] = 0.;//
1518 X_shiftUp[1] = 0.;
1519 X_shiftUp[2] = 0.;
1520 X_shiftUp[3] = 0.;
1521 X_shiftUp[4] = 0.;//-0.015;
1522
1523 Y_shiftUp[0] = 0.;//-0.005 + 0.02-0.005-0.008+0.03;
1524 Y_shiftUp[1] = 0.;//+0.03;
1525 Y_shiftUp[2] = 0.;//+0.006;
1526 Y_shiftUp[3] = 0.;//+.1;
1527 Y_shiftUp[4] = 0.;//+0.002;
1528
1529
1530 return kSUCCESS;
1531}
1532//----------------------------------------------------------------------
1533
1534
1535//---------------------Initialisation-----------------------------------
1536void BmnUpstreamTracking::PrepareArraysToProcessEvent() {
1537 fBmnUpstreamTracksArray->Clear();
1538 OutVector.clear();
1539 vecUpTracks.clear();
1540 vec_points.clear();
1541 vec_twotracks.clear();
1542
1543 NSiTracks = 0;
1544 for (Int_t iPars = 0; iPars < kNumPars; iPars++) {
1545 for (Int_t j = 0; j < kBig; j++) {
1546 par_ab_SiTr[iPars][j] = -999.;
1547 par_SiTr_z[j] = -999.;
1548 }
1549 }
1550
1551 for (Int_t iCh = 0; iCh < kNumChambers; iCh++) {
1552 Nseg_Ch[iCh] = 0;
1553
1554 for (Int_t iPars = 0; iPars < kNumPars; iPars++) {
1555 for (Int_t j = 0; j < kBig; j++) {
1556 par_ab_Ch[iCh][iPars][j] = -999.;
1557 par_Seg_z[iCh][j] = -999.;
1558 }
1559 }
1560 }
1561 for (Int_t ip = 0; ip < kNumPairs; ip++) {
1562 NPCTracks[ip] = 0;
1563 for (Int_t iPars = 0; iPars < kNumPars; iPars++) {
1564 for (Int_t j = 0; j < kBig; j++) {
1565 par_ab_tr[ip][iPars][j] = -999.;
1566 }
1567 }
1568 }
1569 for (Int_t ist = 0; ist < kNumStSi; ist++) {
1570 NSiXYhits[ist] = 0;
1571 for (Int_t ib = 0; ib < kBig; ib++) {
1572 for (Int_t iPars = 0; iPars < kNumPars; iPars++) {
1573 SiXYhits[ist][ib][iPars] = -999.;
1574 }
1575 }
1576 }
1577 for (Int_t ib = 0; ib < kBig; ib++) {
1578 for (Int_t ip = 0; ip < kPoints; ip++) {
1579 for (Int_t i = 0; i < kPoin_Par; i++) {
1580 Points[ib][ip][i] = -999.;
1581 }
1582 }
1583 }
1584
1585}
1586//----------------------------------------------------------------------
1587
1588BmnUpstreamTracking::BmnUpstreamTracking() :
1589 hAx_fitUp(nullptr),
1590 hAy_fitUp(nullptr),
1591 hx_fitUp(nullptr),
1592 hy_fitUp(nullptr),
1593 hchi2_fitUp(nullptr),
1594 hNhitsUp(nullptr),
1595 hdAx_tr_mc_comb(nullptr),
1596 hdX_tr_mc_comb(nullptr),
1597 hdAy_tr_mc_comb(nullptr),
1598 hdY_tr_mc_comb(nullptr),
1599 hdAx_uptr_mc(nullptr),
1600 hdX_uptr_mc(nullptr),
1601 hdAy_uptr_mc(nullptr),
1602 hdY_uptr_mc(nullptr),
1603 hDen_mcuptr(nullptr),
1604 hNum_mcuptr(nullptr),
1605 hEff_mcuptr(nullptr),
1606 hAx_upmc(nullptr),
1607 hAy_upmc(nullptr),
1608 hX_upmc(nullptr),
1609 hY_upmc(nullptr),
1610 hNtr_reco(nullptr),
1611 hNtr_mc(nullptr),
1612 hNrecoTrif2mc(nullptr),
1613 hAngle_reco(nullptr),
1614 hAngle_recoifNmc2(nullptr),
1615 hAngle_recoifNmc2Cases(nullptr),
1616 hNtr_mc_vs_reco(nullptr),
1617 hy_vs_x_Up(nullptr),
1618 hY_vs_Xmctrue(nullptr),
1619 hvertexXYUp(nullptr),
1620 hTyTx_Up(nullptr)
1621{}
1622
1623//-------------------------constructor----------------------------------
1624BmnUpstreamTracking::BmnUpstreamTracking(Bool_t isExp, Int_t runNumber)
1625{
1626 fRunNumber = runNumber;
1627 expData = isExp;
1628 fInputBranchName1 = "BmnSiliconTrack";
1629 fInputBranchName2 = "BmnMwpcTrack";
1630 fInputBranchName3 = "BmnMwpcSegment";
1631 fInputBranchHits = "BmnSiliconHits";
1632 fOutputTracksBranchName = "BmnTrack";
1633 if (!expData){
1634 fRunNumber = 0001;
1635 fInputBranchNameSimTrue = "BmnSiliconHitClean";
1636 fInputBranchName = "BmnMwpcHit";
1637 }
1638
1639}
1640//----------------------------------------------------------------------
1641
1642
1643//-------------------Destructor-----------------------------------------
1644BmnUpstreamTracking::~BmnUpstreamTracking() {
1645
1646 for (Int_t iPars = 0; iPars < kNumPars; iPars++) {
1647 delete[] par_ab_SiTr[iPars];
1648 }
1649 for (Int_t iCh = 0; iCh < kNumChambers; iCh++) {
1650 for (Int_t iPars = 0; iPars < kNumPars; iPars++) {
1651 delete[] par_ab_Ch[iCh][iPars];
1652 }
1653 delete[] par_ab_Ch[iCh];
1654 delete[] par_Seg_z[iCh];
1655 }
1656 for (Int_t ip = 0; ip < kNumPairs; ip++) {
1657 for (Int_t iPars = 0; iPars < kNumPars; iPars++) {
1658 delete[] par_ab_tr[ip][iPars];
1659 }
1660 delete[] par_ab_tr[ip];
1661 delete[] par_ab_trz[ip];
1662 }
1663 for (Int_t ist = 0; ist < kNumStSi; ist++) {
1664 for (Int_t ib = 0; ib < kBig; ib++) {
1665 delete[] SiXYhits[ist][ib];
1666 }
1667 delete[] SiXYhits[ist];
1668 }
1669 for (Int_t ib = 0; ib < kBig; ib++) {
1670 for (Int_t ip = 0; ip < kPoints; ip++) {
1671 delete[] Points[ib][ip];
1672 }
1673 delete[] Points[ib];
1674 }
1675 for (Int_t ip = 0; ip < kPoints; ip++) {
1676 delete[] Xforglfit[ip];
1677 }
1678 for (Int_t i = 0; i < 4; i++) {
1679 delete[] Amatr[i];
1680 }
1681 delete[] X_shiftUp;
1682 delete[] AmatrInverted;
1683 delete[] rhs;
1684 delete[] AmatrArray;
1685 delete[] Amatr;
1686 delete[] line;
1687 delete[] Xforglfit;
1688 delete[] Points;
1689 delete[] NSiXYhits;
1690 delete[] SiXYhits;
1691 delete[] par_ab_trz;
1692 delete[] par_ab_tr;
1693 delete[] par_ab_SiTr;
1694 delete[] par_SiTr_z;
1695 delete[] par_ab_Ch;
1696 delete[] par_Seg_z;
1697 delete[] Nseg_Ch;
1698 delete[] X_shift_seg;
1699 delete[] Z_pair;
1700 delete[] Z_Chamber;
1701 delete[] NPCTracks;
1702 delete[] Y_shiftUp;
1703}
1704//----------------------------------------------------------------------
1705
1706
1707
1708//----------------------------------------------------------------------
1709void BmnUpstreamTracking::Finish() {
1710
1711 if (fDebug) {
1712
1713 hEff_mcuptr->Divide(hNum_mcuptr,hDen_mcuptr,1,1);
1714 printf("BmnUpstreamTracking: write hists to file... ");
1715 fOutputFileName = Form("hUpstreamTracks_run%d.root", fRunNumber);
1716 cout<< fOutputFileName <<endl;
1717 TFile file(fOutputFileName, "RECREATE");
1718 fList.Write();
1719 file.Close();
1720 }
1721 cout << "Work time of the UpstreamTrack finder: " << workTime << " s" << endl;
1722}
1723//----------------------------------------------------------------------
fhList
TList fhList
Definition BmnUpstreamTracking.cxx:28
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition P4_F32vec4.h:34
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
i
int i
Definition P4_F32vec4.h:22
m
__m128 m
Definition P4_F32vec4.h:27
BmnHit::GetIndex
Int_t GetIndex() const
Definition BmnHit.h:37
BmnHit::GetModule
Int_t GetModule()
Definition BmnHit.h:77
BmnHit::GetStation
Short_t GetStation() const
Definition BmnHit.h:45
BmnMwpcHit::GetPdgId
Int_t GetPdgId()
Definition BmnMwpcHit.h:77
BmnMwpcHit::GetMwpcId
Short_t GetMwpcId() const
Definition BmnMwpcHit.h:37
BmnMwpcHit::GetPlaneId
Int_t GetPlaneId()
Definition BmnMwpcHit.h:76
BmnMwpcHit::GetHitId
Int_t GetHitId() const
Definition BmnMwpcHit.h:41
BmnTrack::SetChi2
void SetChi2(Double_t chi2)
Definition BmnTrack.h:97
BmnTrack::SetParamLast
void SetParamLast(FairTrackParam &par)
Definition BmnTrack.h:89
BmnTrack::SetFlag
void SetFlag(Int_t flag)
Definition BmnTrack.h:93
BmnTrack::GetParamFirst
FairTrackParam * GetParamFirst()
Definition BmnTrack.h:72
BmnTrack::GetNHits
Int_t GetNHits() const
Definition BmnTrack.h:44
BmnTrack::SetNHits
void SetNHits(Int_t n)
Definition BmnTrack.h:105
BmnTrack::GetParamLast
FairTrackParam * GetParamLast()
Definition BmnTrack.h:76
BmnUpstreamTracking::compareSegments
static bool compareSegments(const Smatch &a, const Smatch &b)
Definition BmnUpstreamTracking.h:95
BmnUpstreamTracking::Exec
virtual void Exec(Option_t *opt)
Definition BmnUpstreamTracking.cxx:31
BmnUpstreamTracking::Init
virtual InitStatus Init()
Definition BmnUpstreamTracking.cxx:1304
Smatch::params2
Double_t params2[4]
Definition BmnUpstreamTracking.h:61
Smatch::param1
Double_t param1[4]
Definition BmnUpstreamTracking.h:59
Smatch::param2
Double_t param2[4]
Definition BmnUpstreamTracking.h:60
Smatch::Chi2m
Double_t Chi2m
Definition BmnUpstreamTracking.h:58
Smatch::distX
Double_t distX
Definition BmnUpstreamTracking.h:56
Smatch::distY
Double_t distY
Definition BmnUpstreamTracking.h:57
UpTracks::CoordX
Double_t CoordX[5]
Definition BmnUpstreamTracking.h:74
UpTracks::CoordY
Double_t CoordY[5]
Definition BmnUpstreamTracking.h:75
UpTracks::CoordZ
Double_t CoordZ[5]
Definition BmnUpstreamTracking.h:76
UpTracks::param
Double_t param[4]
Definition BmnUpstreamTracking.h:73