18#include "FairLogger.h"
19#include "FairMCPoint.h"
20#include "FairRootManager.h"
27static Double_t workTime = 0.0;
32 if (!IsActive())
return;
34 clock_t tStart = clock();
35 PrepareArraysToProcessEvent();
36 if (fDebug) cout <<
"\n======================== Silicon track finder exec started ===================\n" << endl;
37 if (fDebug) cout <<
"Event number: " << fEventNo << endl;
41 if(expData) StripsReading(DigitsArrayX,DigitsArrayXp);
43 if(fDebug) cout<<
"------------Clustering----------------------"<<endl;
44 if(expData) Clustering(DigitsArrayX, DigitsArrayXp, XClust, XClust_width, sumQx, NclustX, XpClust, XpClust_width, sumQxp, NclustXp);
46 if(fDebug) cout<<
"------------Coordinate Calculation----------"<<endl;
47 CoordinateCalculation(NclustX, NclustXp, XCoord, XpCoord, SigmaX, SigmaXp, XspCoord, XpspCoord, YspCoord, SigmspX, SigmspXp, SigmspY, NhitsXYmod,
48 XClust, XClust_width, sumQx,XpClust, XpClust_width, sumQxp,
49 XspClust, XspClust_width, sumQxsp,XpspClust, XpspClust_width, sumQxpsp, vec_points, Sp_pdg);
53 if(expData) CoordinateAlignment(NhitsXYmod, XspCoord, XpspCoord, YspCoord, NclustX, NclustXp, XCoord, XpCoord);
57 if(fDebug) cout<<
"------------Case1: spatial track------------"<<endl;
58 CountSpatialPoints(NhitsXYmod,Nsp_st);
59 if ( Nsp_st[1] > 0 && Nsp_st[2] > 0 && Nsp_st[3] > 0 ){
60 Case1(NhitsXYmod, XspCoord, XpspCoord, YspCoord, SigmspX, SigmspXp, SigmspY, vec_tracks,
61 XspClust, XspClust_width, sumQxsp,XpspClust, XpspClust_width, sumQxpsp);
62 RecordingTracksAfterSpatialPoints(vec_tracks,CleanTr);
64 CheckPoints(NclustX, NclustXp, XCoord, XpCoord, SigmaX, SigmaXp, XspCoord, XpspCoord, YspCoord, SigmspX, SigmspXp, SigmspY, NhitsXYmod, CleanTr,
65 Nleftoversp, NleftoverX, NleftoverXp,leftoverXsp,leftoverXpsp,leftoverYsp,leftoverXsigsp,leftoverXpsigsp,leftoverYsigsp,leftoverX,leftoverXp,leftoverXsig,leftoverXpsig);
66 if(fDebug) cout<<
"------------Case1.--------------------------"<<endl;
68 if(fDebug) cout<<
"------------Case2: 2spatial points + XorX'---"<<endl;
69 CountSpatialPoints(Nleftoversp,Nsp_st);
70 if ( Nsp_st[3] > 0 && (Nsp_st[1] + Nsp_st[2]) > 0 ) {
71 for (Int_t Istn = 1; Istn < 3; ++Istn) {
72 Case2( Istn, Nleftoversp, leftoverXsp, leftoverXpsp, leftoverYsp, leftoverXsigsp, leftoverXpsigsp, leftoverYsigsp, NleftoverX, NleftoverXp,
73 leftoverX, leftoverXp, leftoverXsig, leftoverXpsig, vec_tracks,
74 XClust, XClust_width, sumQx,XpspClust, XpClust_width, sumQxp,
75 XspClust, XspClust_width, sumQxsp,XpspClust, XpspClust_width, sumQxpsp);
76 RecordingTracks(vec_tracks,CleanTr);
77 CheckPoints(NclustX, NclustXp, XCoord, XpCoord, SigmaX, SigmaXp, XspCoord, XpspCoord, YspCoord, SigmspX, SigmspXp, SigmspY, NhitsXYmod, CleanTr,
78 Nleftoversp, NleftoverX, NleftoverXp,leftoverXsp,leftoverXpsp,leftoverYsp,leftoverXsigsp,leftoverXpsigsp,leftoverYsigsp,leftoverX,leftoverXp,leftoverXsig,leftoverXpsig);
81 if(fDebug) cout<<
"------------Case2.--------------------------"<<endl;
83 if(fDebug) cout<<
"------------Case3:2spatial points + check to the vertex--"<<endl;
84 CountSpatialPoints(Nleftoversp,Nsp_st);
85 if ( Nsp_st[1] > 0 && Nsp_st[2] > 0 ) {
86 Case3( Nleftoversp, leftoverXsp, leftoverXpsp, leftoverYsp, leftoverXsigsp, leftoverXpsigsp, leftoverYsigsp,
87 NleftoverX, NleftoverXp, leftoverX, leftoverXp, leftoverXsig, leftoverXpsig, vec_tracks,
88 XClust, XClust_width, sumQx,XpspClust, XpClust_width, sumQxp,
89 XspClust, XspClust_width, sumQxsp,XpspClust, XpspClust_width, sumQxpsp);
90 RecordingTracks_case3(vec_tracks,CleanTr);
91 CheckPoints(NclustX, NclustXp, XCoord, XpCoord, SigmaX, SigmaXp, XspCoord, XpspCoord, YspCoord, SigmspX, SigmspXp, SigmspY, NhitsXYmod, CleanTr,
92 Nleftoversp, NleftoverX, NleftoverXp,leftoverXsp,leftoverXpsp,leftoverYsp,leftoverXsigsp,leftoverXpsigsp,leftoverYsigsp,leftoverX,leftoverXp,leftoverXsig,leftoverXpsig);
94 if(fDebug) cout<<
"------------Case3.--------------------------"<<endl;
96 if(fDebug) cout<<
"------------Case4:1spatial point in 3st + X/X' in st 1&2-"<<endl;
97 CountSpatialPoints(Nleftoversp,Nsp_st);
99 for (Int_t Istn = 1; Istn < 3; ++Istn) {
100 Case4( Istn, Nleftoversp, leftoverXsp, leftoverXpsp, leftoverYsp, leftoverXsigsp, leftoverXpsigsp, leftoverYsigsp, NleftoverX,
101 NleftoverXp, leftoverX, leftoverXp, leftoverXsig, leftoverXpsig, vec_tracks,
102 XClust, XClust_width, sumQx,XpspClust, XpClust_width, sumQxp,
103 XspClust, XspClust_width, sumQxsp,XpspClust, XpspClust_width, sumQxpsp);
104 RecordingTracks(vec_tracks,CleanTr);
105 CheckPoints(NclustX, NclustXp, XCoord, XpCoord, SigmaX, SigmaXp, XspCoord, XpspCoord, YspCoord, SigmspX, SigmspXp, SigmspY, NhitsXYmod, CleanTr,
106 Nleftoversp, NleftoverX, NleftoverXp,leftoverXsp,leftoverXpsp,leftoverYsp,leftoverXsigsp,leftoverXpsigsp,leftoverYsigsp,leftoverX,leftoverXp,leftoverXsig,leftoverXpsig);
109 if(fDebug)cout<<
"-------------Case4.--------------------------"<<endl;
110 CheckLeftover(Nleftoversp, NleftoverX, NleftoverXp,leftoverXsp,leftoverXpsp,leftoverYsp,leftoverXsigsp,leftoverXpsigsp,leftoverYsigsp,leftoverX,leftoverXp,leftoverXsig,leftoverXpsig, CleanTr,
111 NhitsXYmod, XspCoord, YspCoord, SigmspX, SigmspY);
117 if(fDebug) PrintAllTracks(CleanTr);
119 if(!expData) MCefficiencyCalculation(vec_points, CleanTr);
120 if(fDebug)cout<<
"---------------------------------------------"<<endl;
122 WriteSiliconHits(NclustX, NclustXp, NhitsXYmod, XCoord, XClust_width, XClust, sumQx, XpCoord, XpClust_width, XpClust, sumQxp, XspCoord, YspCoord, XspClust_width, XpspClust_width,
123 XspClust, XpspClust, sumQxsp, sumQxpsp, SigmaX, SigmaXp,SigmspX, SigmspXp, SigmspY, Sp_pdg);
124 WriteSiliconTracks(CleanTr);
126 if (fDebug) cout << endl;
127 if (fDebug) cout <<
"======================== Silicon track finder exec finished ===================" << endl;
128 clock_t tFinish = clock();
129 workTime += ((Double_t) (tFinish - tStart)) / CLOCKS_PER_SEC;
134void BmnSiliconTrackFinder::WriteSiliconTracks(vector<tracksX>& Tr){
135 if (fDebug) cout <<
"--WriteSiliconTracks--"<<endl;
137 Int_t Track_counter = 0;
140 for (
size_t InIter = 0; InIter < Tr.size(); ++InIter) {
145 Double_t z1 = Tr.at(InIter).CoordZ[1] + Zcentr;
146 Double_t x1 = Tr.at(InIter).param[0] * Tr.at(InIter).CoordZ[1] + Tr.at(InIter).param[1];
147 Double_t y1 = Tr.at(InIter).param[2] * Tr.at(InIter).CoordZ[1] + Tr.at(InIter).param[3];
148 Double_t z3 = Tr.at(InIter).CoordZ[3] + Zcentr;
149 Double_t x3 = Tr.at(InIter).param[0] * Tr.at(InIter).CoordZ[3] + Tr.at(InIter).param[1];
150 Double_t y3 = Tr.at(InIter).param[2] * Tr.at(InIter).CoordZ[3] + Tr.at(InIter).param[3];
154 for (Int_t st = 1; st < fNstations; ++st) {
155 if (Tr.at(InIter).CoordX[st] < -900. && Tr.at(InIter).CoordXp[st] > -900.){
156 Tr.at(InIter).CoordY[st] = Tr.at(InIter).param[2] * Tr.at(InIter).CoordZ[st] + Tr.at(InIter).param[3];
157 Tr.at(InIter).SigmaY[st] = 0.105581;
161 BmnSiliconHit* hit =
new ((*fBmnSiliconHitsArray)[Track_counter])
162 BmnSiliconHit(0, TVector3(Tr.at(InIter).CoordX[st]+Shift_toCenterOfMagnetX, Tr.at(InIter).CoordY[st]+Shift_toCenterOfMagnetY, Tr.at(InIter).CoordZ[st]+ Zcentr),
163 TVector3(Tr.at(InIter).SigmaX[st], Tr.at(InIter).SigmaY[st], -1), InIter);
166 hit->
SetModule(Tr.at(InIter).ModNum[st]);
177 FairTrackParam trackParamf;
178 trackParamf.SetPosition(TVector3(x1+Shift_toCenterOfMagnetX, y1+Shift_toCenterOfMagnetY, z1));
179 trackParamf.SetTx(Tr.at(InIter).param[0]+ Shift_toCenterOfMagnetAX + Shift_toCenterOfMagnetAX * Tr.at(InIter).param[0]* Tr.at(InIter).param[0]);
180 trackParamf.SetTy(Tr.at(InIter).param[2]+ Shift_toCenterOfMagnetAY + Shift_toCenterOfMagnetAY * Tr.at(InIter).param[2]* Tr.at(InIter).param[2]);
182 FairTrackParam trackParaml;
183 trackParaml.SetPosition(TVector3(x3+Shift_toCenterOfMagnetX, y3+Shift_toCenterOfMagnetY, z3));
184 trackParaml.SetTx(Tr.at(InIter).param[0]+ Shift_toCenterOfMagnetAX + Shift_toCenterOfMagnetAX * Tr.at(InIter).param[0]* Tr.at(InIter).param[0]);
185 trackParaml.SetTy(Tr.at(InIter).param[2]+ Shift_toCenterOfMagnetAY + Shift_toCenterOfMagnetAY * Tr.at(InIter).param[2]* Tr.at(InIter).param[2]);
188 BmnTrack *track =
new ((*fSiTracks)[fSiTracks->GetEntriesFast()])
BmnTrack();
189 track->
SetChi2(Tr.at(InIter).Chi2);
190 track->
SetNHits(Tr.at(InIter).Nhits);
192 track->
SetNDF(Tr.at(InIter).Pdg);
201void BmnSiliconTrackFinder::WriteSiliconHits(Int_t **NclustX_, Int_t **NclustXp_, Int_t **NhitsXYmod_,
202 Double_t ***XCoord_, Double_t ***XClust_width_, Double_t ***XClust_, Double_t ***sumQx_,
203 Double_t ***XpCoord_, Double_t ***XpClust_width_, Double_t ***XpClust_, Double_t ***sumQxp_,
204 Double_t ***XspCoord_, Double_t ***YspCoord_, Double_t ***XspClust_width_, Double_t ***XpspClust_width_,
205 Double_t ***XspClust_, Double_t ***XpspClust_, Double_t ***sumQxsp_, Double_t ***sumQxpsp_,
206 Double_t ***SigmaX_, Double_t ***SigmaXp_,Double_t ***SigmspX_, Double_t ***SigmspXp_, Double_t ***SigmspY_, Int_t *** Sp_pdg_ ){
208 for (Int_t ist = 1; ist < fNstations; ist++) {
209 for (Int_t imod = 0; imod < fNmodules; imod++) {
210 if (ist == 1 && imod > 3)
continue;
211 if (ist == 2 && imod > 1)
continue;
212 for (Int_t cl = 0; cl < NclustX_[ist][imod]; cl++) {
214 if (XCoord_[ist][imod][cl] < -900.)
continue;
217 BmnSiliconHit(0, TVector3(XCoord_[ist][imod][cl]+Shift_toCenterOfMagnetX, -100., Zstation[ist][imod]),
218 TVector3(SigmaX_[ist][imod][cl], -100, -1), cl);
235 for (Int_t ist = 1; ist < fNstations; ist++) {
236 for (Int_t imod = 0; imod < fNmodules; imod++) {
237 if (ist == 1 && imod > 3)
continue;
238 if (ist == 2 && imod > 1)
continue;
239 for (Int_t cl = 0; cl < NclustXp_[ist][imod]; cl++) {
241 if (XpCoord_[ist][imod][cl] < -900.)
continue;
243 BmnSiliconHit* hit =
new ((*fBmnSiliconHitsXpArray)[countXp])
244 BmnSiliconHit(0, TVector3(XpCoord_[ist][imod][cl]+Shift_toCenterOfMagnetX, -100., Zstation[ist][imod]),
245 TVector3(SigmaXp_[ist][imod][cl], -100, -1), cl);
262 for (Int_t ist = 1; ist < fNstations; ist++) {
263 for (Int_t imod = 0; imod < fNmodules; imod++) {
264 if (ist == 1 && imod > 3)
continue;
265 if (ist == 2 && imod > 1)
continue;
266 for (Int_t cl = 0; cl < NhitsXYmod_[ist][imod]; cl++) {
268 if(fDebug)cout<<
" st "<<ist<<
" hit: X "<<XspCoord_[ist][imod][cl] <<
" Y "<<YspCoord_[ist][imod][cl]<<
" pdg "<<Sp_pdg_[ist][imod][cl]<<endl;
269 if (XspCoord_[ist][imod][cl] < -900.)
continue;
270 BmnSiliconHit* hit =
new ((*fBmnSiliconHitsXYArray)[countXY])
271 BmnSiliconHit(0, TVector3(XspCoord_[ist][imod][cl]+Shift_toCenterOfMagnetX, YspCoord_[ist][imod][cl]+Shift_toCenterOfMagnetY, Zstation[ist][imod]),
272 TVector3(SigmspX_[ist][imod][cl], SigmspY_[ist][imod][cl], -1), cl);
276 hit->
SetPdg(Sp_pdg_[ist][imod][cl]);
290void BmnSiliconTrackFinder::MCefficiencyCalculation(vector<MC_points>& vec, vector<tracksX>& Tr){
292 if (fDebug) cout<<
" ---Silicon MC tracks association--"<<endl;
294 Double_t delt[4] = {-999.,-999.,-999.,-999.};
295 Double_t sig[4] = {0.04, 0.08, 0.05, 0.08};
297 Double_t dmatch = 0.;
298 Double_t dmc_match[kMaxMC];
299 Int_t mc_tr_assoc[kMaxMC];
301 for (Int_t j = 0; j < kMaxMC; j++) {
302 dmc_match[j] = 1000.;
306 Double_t dax = -999.;
307 Double_t day = -999.;
311 if (fDebug) cout<<
" Nmc "<<vec.size()<<
" Nsi "<<Tr.size()<<endl;
313 Int_t Ngood_mc_tracks = 0;
314 Int_t Ngood_reco_tracks = 0;
318 for (
size_t itr = 0; itr < vec.size(); itr++) {
322 if (fDebug && vec.at(itr).Np >= 3 && vec.at(itr).wo3st == 0){
323 hDen_mctrSi->Fill(0);
325 cout<<
" SiDen "<<endl;
326 if (vec.at(itr).Pdg == PDG_Li7) wasLi7 = 1;
327 if (vec.at(itr).Pdg == PDG_He4) wasHe4 = 1;
330 for (
size_t InIter = 0; InIter < Tr.size(); InIter++){
332 delt[0] = vec.at(itr).param[0] - Tr.at(InIter).param[0];
333 delt[1] = vec.at(itr).param[1] - Tr.at(InIter).param[1];
334 delt[2] = vec.at(itr).param[2] - Tr.at(InIter).param[2];
335 delt[3] = vec.at(itr).param[3] - Tr.at(InIter).param[3];
338 if (delt[0] > -900.) hdAx_MC_tr_comb->Fill(dax);
339 if (delt[1] > -900.) hdX_MC_tr_comb ->Fill(dx);
340 if (delt[2] > -900.) hdAy_MC_tr_comb->Fill(day);
341 if (delt[3] > -900.) hdY_MC_tr_comb ->Fill(dy);
345 dmatch = (delt[0]*delt[0])/(sig[0]*sig[0])+ (delt[1]*delt[1])/(sig[1]*sig[1])+
346 (delt[2]*delt[2])/(sig[2]*sig[2])+ (delt[3]*delt[3])/(sig[3]*sig[3]);
348 if ( dmc_match[itr] > dmatch){
349 dmc_match[itr] = dmatch;
350 mc_tr_assoc[itr] = InIter;
358 if (mc_tr_assoc[itr] > -1){
361 if (fDebug) cout <<
" itr "<<itr<<
" Np "<<vec.at(itr).Np<<
" mc_Id "<<vec.at(itr).Id<<
362 " ax_mc "<<vec.at(itr).param[0]<<
" reco_ind "<<mc_tr_assoc[itr]<<
" ax "<<Tr.at(mc_tr_assoc[itr]).param[0]<<
363 " dmc_match "<<dmc_match[itr]<<endl;
365 if (dax > -900.) hdAx_MC_tr->Fill(dax);
366 if (dx > -900.) hdX_MC_tr ->Fill(dx);
367 if (day > -900.) hdAy_MC_tr->Fill(day);
368 if (dy > -900.) hdY_MC_tr ->Fill(dy);
373 if (wasLi7 && wasHe4) {
374 if (fDebug) cout <<
" wasLi7 && wasHe4 "<<wasLi7<<
" "<<wasHe4<<endl;
375 hDen_mcreaction->Fill(0);
379 if (fDebug) cout<<
"Si reject poorly chosen association segments "<<endl;
380 for (
size_t itr = 0; itr < vec.size(); itr++) {
381 if (mc_tr_assoc[itr] == -1)
continue;
383 for (
size_t itr2 = 0; itr2 < vec.size(); itr2++) {
384 if (itr2 == itr)
continue;
385 if (mc_tr_assoc[itr2] == -1)
continue;
387 if (mc_tr_assoc[itr] == mc_tr_assoc[itr2]){
388 if (dmc_match[itr2] > 2*dmc_match[itr] ) mc_tr_assoc[itr2] = -1;
397 if (fDebug) cout<<
" mc_Id "<<vec.at(itr).Id<<
" assoc "<<mc_tr_assoc[itr]<<
" pdg "<<vec.at(itr).Pdg<<endl;
398 if (fDebug && mc_tr_assoc[itr] > -1 && vec.at(itr).Np >= 3 && vec.at(itr).wo3st == 0){
401 hNum_mctrSi->Fill(0);
406 if (fDebug) cout<<
" Nassoc "<<Nassoc<<endl;
407 if (Nassoc == 2 && wasLi7 && wasHe4){
408 if (fDebug) cout<<
" hNum_mcreaction "<<endl;
409 hNum_mcreaction->Fill(0);
412 hNtrsi_mc ->Fill(Ngood_mc_tracks);
413 hNtrsi_reco->Fill(Ngood_reco_tracks);
414 hNtrsi_mc_vs_reco ->Fill(Ngood_reco_tracks,Ngood_mc_tracks);
417 for (
size_t InIter = 0; InIter < Tr.size(); InIter++){
418 for (
size_t itr = 0; itr < vec.size(); itr++) {
419 if ( mc_tr_assoc[itr] == (Int_t)InIter){
420 Tr.at(InIter).Pdg = vec.at(itr).Pdg;
424 for (
size_t InIter = 0; InIter < Tr.size(); InIter++){
425 if (fDebug) cout<<
" InIter "<<InIter<<
" Pdg "<<Tr.at(InIter).Pdg<<endl;
429void BmnSiliconTrackFinder::StripsReading(Double_t ***DigitsArrayX_, Double_t ***DigitsArrayXp_){
431 Int_t stat, mod, layer, strip, new_strip; Double_t Ampl;
432 for (Int_t iDig = 0; iDig < fBmnSiDigitsArray->GetEntriesFast(); ++iDig) {
438 if ((strip = digit ->GetStripNumber() - 2) < 0)
continue;
443 DigitsArrayX_[stat][mod][strip] = Ampl;
447 if (fRunPeriod == 7 && stat == 1) {
448 new_strip = nstripXpsm - strip;
449 if (new_strip < 0)
continue;
450 DigitsArrayXp_[stat][mod][new_strip] = Ampl;
452 new_strip = nstripXp - strip;
453 if (new_strip < 0)
continue;
454 DigitsArrayXp_[stat][mod][new_strip] = Ampl;
467void BmnSiliconTrackFinder::GetXYspatial(Int_t **NClX, Int_t **NClXp, Double_t ***XCoor, Double_t ***XpCoor, Double_t ***Xsp,
468 Double_t ***Xpsp, Double_t ***Ysp, Int_t **NXYsp, Double_t ***SigmX, Double_t ***SigmXp, Double_t ***SigspX, Double_t ***SigspXp, Double_t ***SigspY,
469 Double_t ***XCl, Double_t ***XCl_width, Double_t ***sQx, Double_t ***XpCl, Double_t ***XpCl_width, Double_t ***sQxp,
470 Double_t ***XspCl, Double_t ***XspCl_width, Double_t ***sQxsp, Double_t ***XpspCl, Double_t ***XpspCl_width, Double_t ***sQxpsp){
472 SensitiveAreaY =
new Double_t[fNstations];
474 for (Int_t istat = 1; istat < fNstations; istat++) {
475 SensitiveAreaY[istat] = 2 * 6.0705 + 0.1148 + 0.2303;
477 if (fRunPeriod == 7){
478 SensitiveAreaY[1] = 6.0705 + 0.1148;
480 SensitiveAreaY[0] = 0.;
481 Double_t YCoor_cand = -1.;
483 for (Int_t istat = 1; istat < fNstations; istat++) {
484 for (Int_t imod = 0; imod < fNmodules; imod++) {
485 NXYsp[istat][imod] = 0;
487 for (Int_t clx = 0; clx < NClX[istat][imod]; ++clx) {
488 for (Int_t clxp = 0; clxp < NClXp[istat][imod]; ++clxp) {
489 if (XCoor[istat][imod][clx] > -900. && XpCoor[istat][imod][clxp] > -900.) {
492 YCoor_cand = (XpCoor[istat][imod][clxp] - XCoor[istat][imod][clx]) / tg2_5;
496 if (YCoor_cand >= 0.1148 && YCoor_cand <= SensitiveAreaY[istat]) {
498 if(fRunPeriod == 7 && istat > 1 && YCoor_cand > 6.0705 + 0.1148 && YCoor_cand < 6.0705 + 0.1148 + 0.2303 )
continue;
501 if (fRunPeriod == 7 && istat > 1 && YCoor_cand <= half_module) {
502 YCoor_cand -= shift_after_zigzag;
505 Xsp[istat][imod][NXYsp[istat][imod]] = XCoor[istat][imod][clx];
506 SigspX[istat][imod][NXYsp[istat][imod]] = SigmX[istat][imod][clx];
507 XspCl[istat][imod][NXYsp[istat][imod]] = XCl[istat][imod][clx];
508 XspCl_width[istat][imod][NXYsp[istat][imod]] = XCl_width[istat][imod][clx];
509 sQxsp[istat][imod][NXYsp[istat][imod]] = sQx[istat][imod][clx];
511 Xpsp[istat][imod][NXYsp[istat][imod]] = XpCoor[istat][imod][clxp];
512 SigspXp[istat][imod][NXYsp[istat][imod]] = SigmXp[istat][imod][clxp];
513 XpspCl[istat][imod][NXYsp[istat][imod]] = XpCl[istat][imod][clxp];
514 XpspCl_width[istat][imod][NXYsp[istat][imod]]= XpCl_width[istat][imod][clxp];
515 sQxpsp[istat][imod][NXYsp[istat][imod]] = sQxp[istat][imod][clxp];
517 Ysp[istat][imod][NXYsp[istat][imod]] = YCoor_cand;
518 SigspY[istat][imod][NXYsp[istat][imod]] =
sqrt( pow (SigmX[istat][imod][clx],2) + pow(SigmXp[istat][imod][clxp],2) ) / tg2_5;
519 NXYsp[istat][imod]++;
538void BmnSiliconTrackFinder::Clustering( Double_t*** DigitsArrayX_, Double_t *** DigitsArrayXp_,
539 Double_t *** XClust_, Double_t *** XClust_width_, Double_t *** sumQx_, Int_t **NclustX_,
540 Double_t *** XpClust_, Double_t *** XpClust_width_, Double_t *** sumQxp_, Int_t **NclustXp_){
559 for (Int_t istat = 1; istat < fNstations; istat++) {
560 for (Int_t imod = 0; imod < fNmodules; imod++) {
561 for (Int_t istr = 0; istr < fNstrips; istr++) {
562 if (istat == 1 && imod > 3)
continue;
563 if (istat == 2 && imod > 1)
continue;
565 if (fDebug && istat == 1 && DigitsArrayX_[istat][imod][istr] > 0.) hoccupancyX1.at(imod) ->Fill(istr);
566 if (fDebug && istat == 2 && DigitsArrayX_[istat][imod][istr] > 0.) hoccupancyX2.at(imod) ->Fill(istr);
567 if (fDebug && istat == 3 && DigitsArrayX_[istat][imod][istr] > 0.) hoccupancyX3.at(imod) ->Fill(istr);
569 if (fDebug && istat == 1 && DigitsArrayXp_[istat][imod][istr] > 0.) hoccupancyXp1.at(imod) ->Fill(istr);
570 if (fDebug && istat == 2 && DigitsArrayXp_[istat][imod][istr] > 0.) hoccupancyXp2.at(imod) ->Fill(istr);
571 if (fDebug && istat == 3 && DigitsArrayXp_[istat][imod][istr] > 0.) hoccupancyXp3.at(imod) ->Fill(istr);
573 if ( DigitsArrayX_[istat][imod][istr] == 0. && DigitsArrayX_[istat][imod][istr-1] > Cut_AmplX && DigitsArrayX_[istat][imod][istr + 1] > Cut_AmplX)
574 DigitsArrayX_[istat][imod][istr] = 0.5 * ( DigitsArrayX_[istat][imod][istr-1] + DigitsArrayX_[istat][imod][istr + 1] ) ;
576 if (istat == 1 && imod == 0 ) {
578 if ( DigitsArrayX_[istat][imod][404] > Cut_AmplX ) DigitsArrayX_[istat][imod][405] = 0.5* DigitsArrayX_[istat][imod][404] ;
579 if ( DigitsArrayX_[istat][imod][413] > Cut_AmplX ) DigitsArrayX_[istat][imod][412] = 0.5* DigitsArrayX_[istat][imod][413] ;
582 if ( DigitsArrayXp_[istat][imod][istr] == 0. && DigitsArrayXp_[istat][imod][istr-1] > 120. && DigitsArrayXp_[istat][imod][istr + 1] > 120.)
583 DigitsArrayXp_[istat][imod][istr] = 0.5 * ( DigitsArrayXp_[istat][imod][istr-1] + DigitsArrayXp_[istat][imod][istr + 1] );
588 Int_t NfirstX, NlastX, NfirstXp, NlastXp, ClusterSizeX, ClusterSizeXp;
589 Double_t SumAmpl, SumAmplXp;
591 for (Int_t istat = 1; istat < fNstations; istat++) {
593 for (Int_t imod = 0; imod < fNmodules; imod++) {
601 for (Int_t istr = 0; istr < fNstrips; istr++){
606 if (NfirstX < 0 && DigitsArrayX_[istat][imod][istr] == 0.)
continue;
607 if (NfirstX < 0 && DigitsArrayX_[istat][imod][istr] > 0.) NfirstX = istr;
608 if (NfirstX >= 0 && DigitsArrayX_[istat][imod][istr + 1] == 0. ) NlastX = istr;
610 if (NfirstX >= 0 && NlastX > 0) {
612 Int_t max_ampl_strip = -1,
613 num_max_ampl_strip = -1,
616 num_centr_strip = -1 ;
618 ClusterSizeX = NlastX - NfirstX + 1;
621 Int_t Counter_overflow = 0;
623 for (Int_t is = NfirstX; is < NlastX + 1; is++) {
624 Ampl_strX[is - NfirstX] = DigitsArrayX_[istat][imod][is];
626 if ( DigitsArrayX_[istat][imod][is] > Cut_overflow) Counter_overflow++;
628 if ( DigitsArrayX_[istat][imod][is] > max_ampl_strip ) {
629 max_ampl_strip = DigitsArrayX_[istat][imod][is];
630 num_max_ampl_strip = is;
633 CoG = FindClusterCenter(Ampl_strX, ClusterSizeX, SumAmpl);
637 for (Int_t is = NfirstX; is < NlastX + 1; is++) {
638 if ( max_ampl_strip > 0) {
641 Int_t dif_strip = is - num_max_ampl_strip;
642 if ( dif_strip == 0) num_centr_strip = is;
650 if ( num_centr_strip > 0 && DigitsArrayX_[istat][imod][num_centr_strip] < Cut_AmplStripX) {
657 if ( SumAmpl < Cut_AmplX ) {
665 XClust_[istat][imod][NclustX_[istat][imod]] = NfirstX + CoG;
666 XClust_width_[istat][imod][NclustX_[istat][imod]] = ClusterSizeX;
667 sumQx_[istat][imod][NclustX_[istat][imod]] = SumAmpl;
670 NclustX_[istat][imod] ++;
680 if (NclustX_[istat][imod] == 0)
continue;
687 for (Int_t istr = 0; istr < fNstrips; istr++) {
692 if (NfirstXp < 0 && DigitsArrayXp_[istat][imod][istr] == 0.)
continue;
693 if (NfirstXp < 0 && DigitsArrayXp_[istat][imod][istr] > 0.) NfirstXp = istr;
694 if (NfirstXp >= 0 && DigitsArrayXp_[istat][imod][istr + 1] == 0.) NlastXp = istr;
696 if (NfirstXp >= 0 && NlastXp > 0) {
698 Int_t max_ampl_strip = -1,
699 num_max_ampl_strip = -1,
702 num_centr_strip = -1 ;
704 ClusterSizeXp = NlastXp - NfirstXp + 1;
707 Int_t Counter_overflow = 0;
709 for (Int_t is = NfirstXp; is < NlastXp + 1; is++) {
710 Ampl_strXp[is - NfirstXp] = DigitsArrayXp_[istat][imod][is];
712 if ( DigitsArrayXp_[istat][imod][is] > Cut_overflow) Counter_overflow++;
714 if ( DigitsArrayXp_[istat][imod][is] > max_ampl_strip ) {
715 max_ampl_strip = DigitsArrayXp_[istat][imod][is];
716 num_max_ampl_strip = is;
720 CoG = FindClusterCenter(Ampl_strXp, ClusterSizeXp, SumAmplXp);
723 for (Int_t is = NfirstX; is < NlastX + 1; is++) {
724 if ( max_ampl_strip > 0) {
726 Int_t dif_strip = is - num_max_ampl_strip;
727 if ( dif_strip == 0) num_centr_strip = is;
736 if ( num_centr_strip > 0 && DigitsArrayXp_[istat][imod][num_centr_strip] < Cut_AmplStripXp) {
743 if ( SumAmplXp < Cut_AmplXp) {
756 XpClust_[istat][imod][NclustXp_[istat][imod]] = NfirstXp + CoG;
757 XpClust_width_[istat][imod][NclustXp_[istat][imod]] = ClusterSizeXp;
758 sumQxp_[istat][imod][NclustXp_[istat][imod]] = SumAmplXp;
759 NclustXp_[istat][imod]++;
776void BmnSiliconTrackFinder::CoordinateCalculation(Int_t **NclustX_, Int_t **NclustXp_, Double_t ***XCoord_, Double_t ***XpCoord_, Double_t ***SigmaX_, Double_t ***SigmaXp_,
777 Double_t ***XspCoord_, Double_t ***XpspCoord_, Double_t ***YspCoord_, Double_t ***SigmspX_, Double_t ***SigmspXp_, Double_t ***SigmspY_,
779 Double_t ***XClust_, Double_t ***XClust_width_, Double_t ***sumQx_, Double_t ***XpClust_, Double_t ***XpClust_width_, Double_t ***sumQxp_,
780 Double_t ***XspClust_, Double_t ***XspClust_width_, Double_t ***sumQxsp_, Double_t ***XpspClust_, Double_t ***XpspClust_width_, Double_t ***sumQxpsp_,
781 vector<MC_points>& vec, Int_t *** Sp_pdg_ ){
785 for (Int_t ist = 1; ist < fNstations; ist++) {
786 for (Int_t imod = 0; imod < fNmodules; imod++) {
787 if (ist == 1 && imod > 3)
continue;
788 if (ist == 2 && imod > 1)
continue;
790 for (Int_t cl = 0; cl < NclustX_[ist][imod]; cl++) {
791 XCoord_[ist][imod][cl] = deltaX * XClust_[ist][imod][cl];
792 if (XClust_width[ist][imod][cl] > 1 ) SigmaX_[ist][imod][cl] = (XClust_width_[ist][imod][cl] * deltaX) / (2.*sq12);
793 else SigmaX_[ist][imod][cl] = deltaX / sq12;
800 for (Int_t cl = 0; cl < NclustXp_[ist][imod]; cl++) {
801 XpCoord_[ist][imod][cl] = deltaXp * XpClust_[ist][imod][cl];
802 if (XpClust_width[ist][imod][cl] > 1 ) SigmaXp_[ist][imod][cl] = (XpClust_width_[ist][imod][cl] * deltaXp) / (2.*sq12);
803 else SigmaXp_[ist][imod][cl] = deltaXp / sq12;
836 GetXYspatial(NclustX_, NclustXp_, XCoord_, XpCoord_, XspCoord_, XpspCoord_, YspCoord_, NhitsXYmod_, SigmaX_, SigmaXp_, SigmspX_, SigmspXp_, SigmspY_,
837 XClust_, XClust_width_, sumQx_, XpClust_, XpClust_width_, sumQxp_, XspClust_, XspClust_width_, sumQxsp_, XpspClust_, XpspClust_width_, sumQxpsp_);
886 if (fDebug) cout<<
" Sim: "<<endl;
887 for (Int_t iMC = 0; iMC < fBmnHitsArray2->GetEntriesFast(); ++iMC) {
890 Double_t x_MC = hit->GetX();
891 Double_t xp_MC = hit->GetY();
892 Double_t z_MC = hit->GetZ();
895 Double_t sigx_MC = hit->GetDx();
896 Double_t sigxp_MC = hit->GetDy();
898 if (sigx_MC < 0 ) sigx_MC = -sigx_MC;
899 if (sigxp_MC < 0 ) sigxp_MC= -sigxp_MC;
901 if (z_MC > -320.){ st_mc = 3;
if (fDebug) hSi_st3mc->Fill(x_MC,(xp_MC - x_MC)/tg2_5 );}
902 if (z_MC < -434. && z_MC > -436.) {st_mc = 2;
if (fDebug) hSi_st2mc->Fill(x_MC,(xp_MC - x_MC)/tg2_5 );}
903 if (z_MC < -438.) {st_mc = 1;
if (fDebug) hSi_st1mc->Fill(x_MC,(xp_MC - x_MC)/tg2_5 );}
905 if (st_mc == 2) imod = 1;
907 XCoord_[st_mc][imod][NclustX[st_mc][imod]] = x_MC;
908 XpCoord_[st_mc][imod][NclustXp[st_mc][imod]] = xp_MC;
909 SigmaX_[st_mc][imod][NclustXp[st_mc][imod]] = sigx_MC;
910 SigmaXp_[st_mc][imod][NclustXp[st_mc][imod]] = sigxp_MC;
911 NclustX_[st_mc][imod]++;
912 NclustXp_[st_mc][imod]++;
914 XspCoord_[st_mc][imod][NhitsXYmod[st_mc][imod]] = x_MC;
915 if (xp_MC > -900.) XpspCoord_[st_mc][imod][NhitsXYmod[st_mc][imod]] = xp_MC;
916 YspCoord_[st_mc][imod][NhitsXYmod[st_mc][imod]] = (xp_MC - x_MC)/tg2_5;
917 Sp_pdg_[st_mc][imod][NhitsXYmod[st_mc][imod]] = hit->
GetType();
918 SigmspX_[st_mc][imod][NhitsXYmod[st_mc][imod]] = sigx_MC;
919 SigmspXp_[st_mc][imod][NhitsXYmod[st_mc][imod]] = sigxp_MC;
920 SigmspY_[st_mc][imod][NhitsXYmod[st_mc][imod]] =
sqrt( pow (sigx_MC,2) + pow(sigxp_MC,2) ) / tg2_5;
921 if (fDebug) cout<<
" Id "<<Id<<
" x "<<x_MC<<
" xp "<<xp_MC<<
" y "<<YspCoord_[st_mc][imod][NhitsXYmod[st_mc][imod]] <<
" z_MC "<<z_MC<<
" st_mc "<<st_mc<<
" pdg "<<hit->
GetType()<<endl;
923 NhitsXYmod_[st_mc][imod]++;
927 if (fDebug) cout<<
"True MC BmnSiliconHitClean"<<endl;
932 int Np_in_3st[kMaxMC];
933 int mcTracksArray[kMaxMC];
934 int mcPdgCode[kMaxMC];
935 double Xmc[kMaxMC][fNstations];
936 double Ymc[kMaxMC][fNstations];
937 double Zmc[kMaxMC][fNstations];
938 for (Int_t Id= 0; Id < kMaxMC; Id++) {
941 mcTracksArray[Id] = -1;
943 for (Int_t
i = 0;
i < fNstations;
i++) {
950 for (Int_t iMC = 0; iMC < fBmnHitsArray->GetEntriesFast(); ++iMC) {
953 Double_t x_MC = hit->GetX();
954 Double_t y_MC = hit->GetY();
955 Double_t z_MC = hit->GetZ();
958 Int_t pdgCode = hit->
GetType();
960 if (tr_before != trackId_MC) {
962 mcTracksArray[Nmc_tracks] = trackId_MC;
964 tr_before = trackId_MC;
966 if (z_MC > -320.) {st_mc = 3; Np_in_3st[Nmc_tracks]++; }
967 if (z_MC < -434. && z_MC > -436. ) st_mc = 2;
968 if (z_MC < -438.) st_mc = 1;
969 if (fDebug) cout<<
" pdgId "<<trackId_MC<<
" pdgCode "<<pdgCode<<
" x "<<x_MC<<
" y "<<y_MC<<
" z "<<z_MC<<
" st_mc "<<st_mc<<endl;
971 Xmc[Nmc_tracks][st_mc] = x_MC;
972 Ymc[Nmc_tracks][st_mc] = y_MC;
973 Zmc[Nmc_tracks][st_mc] = z_MC;
974 mcPdgCode[Nmc_tracks] = pdgCode;
975 Nst_mc[Nmc_tracks]++;
980 for (Int_t Id= 0; Id < kMaxMC; Id++) {
981 if (Nst_mc[Id] > 0 ){
982 for (Int_t
i = 0;
i < fNstations;
i++) {
983 tmpTr.x[
i] = Xmc[Id][
i];
984 tmpTr.y[
i] = Ymc[Id][
i];
985 tmpTr.z[
i] = Zmc[Id][
i];
987 tmpTr.Id = mcTracksArray[Id];
988 tmpTr.Np = Nst_mc[Id];
989 tmpTr.Pdg= mcPdgCode[Id];
990 if (Nst_mc[Id] >= 2 && Np_in_3st[Id] > 0 ) vec.push_back(tmpTr);
994 if (fDebug) cout<<
" vec_points.size() "<<vec.size()<<endl;
996 for (
size_t itr = 0; itr < vec.size(); itr++) {
998 if (vec.at(itr).x[3] < -900.) vec.at(itr).wo3st = 1;
999 if (vec.at(itr).x[1] > -900.){
1000 vec.at(itr).param[0] = (vec.at(itr).x[1] - vec.at(itr).x[3])/ (vec.at(itr).z[1] - vec.at(itr).z[3]);
1001 vec.at(itr).param[1] = (vec.at(itr).x[1] + vec.at(itr).x[3])*0.5;
1002 vec.at(itr).param[2] = (vec.at(itr).y[1] - vec.at(itr).y[3])/ (vec.at(itr).z[1] - vec.at(itr).z[3]);
1003 vec.at(itr).param[3] = (vec.at(itr).y[1] + vec.at(itr).y[3])*0.5;
1005 vec.at(itr).param[0] = (vec.at(itr).x[2] - vec.at(itr).x[3])/ (vec.at(itr).z[2] - vec.at(itr).z[3]);
1006 vec.at(itr).param[1] = (vec.at(itr).x[2] + vec.at(itr).x[3])*0.5;
1007 vec.at(itr).param[2] = (vec.at(itr).y[2] - vec.at(itr).y[3])/ (vec.at(itr).z[2] - vec.at(itr).z[3]);
1008 vec.at(itr).param[3] = (vec.at(itr).y[2] + vec.at(itr).y[3])*0.5;
1010 if (fDebug) cout<<
" mcId "<<vec.at(itr).Id<<
" Pdg "<<vec.at(itr).Pdg<<
" Ax "<<vec.at(itr).param[0]<<
" bx "<<vec.at(itr).param[1]<<
" Ay "<<vec.at(itr).param[2]<<
" by "<<vec.at(itr).param[3]<<
" Np "<<vec.at(itr).Np<<endl;
1011 if (fDebug) cout<<
" Xt "<<vec.at(itr).param[0]*( Z0_SRC_target - Zcentr ) + vec.at(itr).param[1]<<
1012 " Yt "<<vec.at(itr).param[2]*( Z0_SRC_target - Zcentr ) + vec.at(itr).param[3]<<endl;
1014 if (fDebug) hAxsi_mctrue->Fill(vec.at(itr).param[0]);
1015 if (fDebug) hAysi_mctrue->Fill(vec.at(itr).param[2]);
1016 if (fDebug) hBxsi_mctrue->Fill(vec.at(itr).param[1]);
1017 if (fDebug) hBysi_mctrue->Fill(vec.at(itr).param[3]);
1028void BmnSiliconTrackFinder::CoordinateAlignment( Int_t **NhitsXYmod_, Double_t ***XspCoord_, Double_t ***XpspCoord_, Double_t ***YspCoord_,
1029 Int_t **NclustX_, Int_t **NclustXp_, Double_t ***XCoord_, Double_t ***XpCoord_){
1033 for (Int_t ist = 1; ist < fNstations; ist++) {
1035 for (Int_t imod = 0; imod < fNmodules; imod++) {
1037 if (ist == 1 && imod > 3 )
continue;
1038 if (ist == 2 && imod > 1 )
continue;
1043 if (fRunPeriod == 7){
1045 for(Int_t ist = 1; ist < fNstations; ist++) {
1046 for (Int_t imod = 0; imod < fNmodules; imod++) {
1047 if (ist == 1 && imod > 3 )
continue;
1048 if (ist == 2 && imod > 1 )
continue;
1049 for (Int_t cl = 0; cl < NhitsXYmod_[ist][imod]; cl++) {
1051 if( (ist == 2 || ist == 3)) {
1052 if( YspCoord_[ist][imod][cl] > 0. && YspCoord_[ist][imod][cl] <= 0.1151)
continue;
1053 if( YspCoord_[ist][imod][cl] < 0. && YspCoord_[ist][imod][cl] >= -0.1151)
continue;
1054 if( YspCoord_[ist][imod][cl] > 0. ) YspCoord_[ist][imod][cl] += 0.1151;
1058 (ist == 1 && (imod == 1 || imod == 3)) ||
1059 (ist == 2 && imod == 1) ||
1060 (ist == 3 && (imod == 2 || imod ==3 || imod ==6 || imod ==7) )
1063 XspCoord_[ist][imod][cl] = -XspCoord_[ist][imod][cl];
1064 XpspCoord_[ist][imod][cl] = -XpspCoord_[ist][imod][cl];
1072 (ist == 1 && (imod == 2 || imod == 3) )||
1073 (ist == 3 && imod > 3 )
1075 YspCoord_[ist][imod][cl] = -YspCoord_[ist][imod][cl];
1083 XspCoord_[ist][imod][cl] += shiftX[ist][imod];
1084 YspCoord_[ist][imod][cl] += shiftY[ist][imod];
1087 XspCoord_[ist][imod][cl] = -XspCoord_[ist][imod][cl];
1092 XpspCoord_[ist][imod][cl] += shiftY[ist][imod]*(-Angle(ist,imod))*tg2_5 + shiftX[ist][imod];
1094 XpspCoord_[ist][imod][cl] = -XpspCoord_[ist][imod][cl];
1096 XpspCoord_[ist][imod][cl] += shiftStYtoGlob[ist]*(Angle(ist,imod))*tg2_5 + shiftStXtoGlob[ist];
1099 XspCoord_[ist][imod][cl] += shiftStXtoGlob[ist];
1100 YspCoord_[ist][imod][cl] += shiftStYtoGlob[ist];
1112 for (Int_t ist = 1; ist < fNstations; ist++) {
1113 for (Int_t imod = 0; imod < fNmodules; imod++) {
1114 if (ist == 1 && imod > 3 )
continue;
1115 if (ist == 2 && imod > 1 )
continue;
1117 for (Int_t cl = 0; cl < NclustX_[ist][imod]; cl++) {
1119 (ist == 1 && (imod == 1 || imod == 3)) ||
1120 (ist == 2 && imod == 1) ||
1121 (ist == 3 && (imod == 2 || imod ==3 || imod ==6 || imod ==7) )
1123 XCoord_[ist][imod][cl]= -XCoord_[ist][imod][cl];
1125 XCoord_[ist][imod][cl] += shiftX[ist][imod];
1128 XCoord_[ist][imod][cl] = -XCoord_[ist][imod][cl];
1129 XCoord_[ist][imod][cl] += shiftStXtoGlob[ist];
1132 for (Int_t cl = 0; cl < NclustXp_[ist][imod]; cl++) {
1135 (ist == 1 && (imod == 1 || imod == 3)) ||
1136 (ist == 2 && imod == 1) ||
1137 (ist == 3 && (imod == 2 || imod ==3 || imod ==6 || imod ==7) )
1139 XpCoord_[ist][imod][cl]= -XpCoord_[ist][imod][cl];
1143 XpCoord_[ist][imod][cl] += shiftY[ist][imod]*(-Angle(ist,imod))*tg2_5 + shiftX[ist][imod];
1148 XpCoord_[ist][imod][cl] = -XpCoord_[ist][imod][cl];
1151 XpCoord_[ist][imod][cl] += shiftStYtoGlob[ist]*(Angle(ist,imod))*tg2_5 + shiftStXtoGlob[ist];
1158 if (fRunPeriod == 8){
1161 for(Int_t ist = 1; ist < fNstations; ist++) {
1162 for (Int_t imod = 0; imod < fNmodules; imod++) {
1163 for (Int_t cl = 0; cl < NhitsXYmod_[ist][imod]; cl++) {
1167 XspCoord_[ist][imod][cl] = -XspCoord_[ist][imod][cl];
1168 XpspCoord_[ist][imod][cl] = -XpspCoord_[ist][imod][cl];
1174 YspCoord_[ist][imod][cl] = -YspCoord_[ist][imod][cl];
1180 XspCoord_[ist][imod][cl] += shiftX[ist][imod];
1181 YspCoord_[ist][imod][cl] += shiftY[ist][imod];
1182 XpspCoord_[ist][imod][cl] += shiftY[ist][imod]*(-Angle(ist,imod))*tg2_5 + shiftX[ist][imod];
1186 if( ist == 1 || ist == 3){
1187 Xtmp = XspCoord_[ist][imod][cl];
1188 XspCoord_[ist][imod][cl] = YspCoord_[ist][imod][cl];
1189 YspCoord_[ist][imod][cl] = Xtmp;
1191 if( ist == 2 || ist == 4){
1192 XspCoord_[ist][imod][cl] = -XspCoord_[ist][imod][cl];
1193 YspCoord_[ist][imod][cl] = -YspCoord_[ist][imod][cl];
1197 XspCoord_[ist][imod][cl] += shiftStXtoGlob[ist];
1198 YspCoord_[ist][imod][cl] += shiftStYtoGlob[ist];
1199 XpspCoord_[ist][imod][cl] += shiftStYtoGlob[ist]*(Angle(ist,imod))*tg2_5 + shiftStXtoGlob[ist];
1211 for (Int_t ist = 1; ist < fNstations; ist++) {
1212 for (Int_t imod = 0; imod < fNmodules; imod++) {
1213 for (Int_t cl = 0; cl < NclustX_[ist][imod]; cl++) {
1215 XCoord_[ist][imod][cl]= -XCoord_[ist][imod][cl];
1217 XCoord_[ist][imod][cl] += shiftX[ist][imod];
1219 if( ist == 2 || ist == 4){
1220 XCoord_[ist][imod][cl] = - XCoord_[ist][imod][cl];
1222 XCoord_[ist][imod][cl] += shiftStXtoGlob[ist];
1225 for (Int_t cl = 0; cl < NclustXp_[ist][imod]; cl++) {
1228 XpCoord_[ist][imod][cl]= -XpCoord_[ist][imod][cl];
1232 XpCoord_[ist][imod][cl] += shiftY[ist][imod]*(-Angle(ist,imod))*tg2_5 + shiftX[ist][imod];
1234 if( ist == 2 || ist == 4){
1235 XpCoord_[ist][imod][cl] = - XpCoord_[ist][imod][cl];
1239 XpCoord_[ist][imod][cl] += shiftStYtoGlob[ist]*(Angle(ist,imod))*tg2_5 + shiftStXtoGlob[ist];
1246 if(fDebug) cout<<
" --after shift-- "<<endl;
1248 for (Int_t ist = 1; ist < fNstations; ist++) {
1249 if (fDebug) cout <<
" XCoord_[" << ist <<
"] ";
1250 for (Int_t imod = 0; imod < fNmodules; imod++) {
1251 for (Int_t cl = 0; cl < NclustX_[ist][imod]; cl++) {
1252 if (fDebug) cout << XCoord_[ist][imod][cl] <<
" ";
1255 if (fDebug) cout << endl;
1257 if (fDebug) cout << endl;
1258 for (Int_t ist = 1; ist < fNstations; ist++) {
1259 if (fDebug) cout <<
" XpCoord_[" << ist <<
"] ";
1260 for (Int_t imod = 0; imod < fNmodules; imod++) {
1261 for (Int_t cl = 0; cl < NclustXp_[ist][imod]; cl++) {
1262 if (fDebug) cout << XpCoord_[ist][imod][cl] <<
" ";
1265 if (fDebug) cout << endl;
1267 if(fDebug) cout<<
" --spatial-- "<<endl;
1269 for(Int_t ist = 1; ist < fNstations; ist++) {
1270 if(fDebug) cout<<
" Xspatial["<<ist<<
"] ";
1271 for (Int_t imod = 0; imod < fNmodules; imod++) {
1272 if (ist == 1 && imod > 3 )
continue;
1273 if (ist == 2 && imod > 1 )
continue;
1275 for (Int_t cl = 0; cl < NhitsXYmod_[ist][imod]; cl++) {
1276 if(fDebug) cout<<XspCoord_[ist][imod][cl]<<
" ";
1279 if(fDebug) cout<<endl;
1281 if(fDebug) cout<<endl;
1283 for(Int_t ist = 1; ist < fNstations; ist++) {
1284 if(fDebug) cout<<
" Xpspatial["<<ist<<
"] ";
1285 for (Int_t imod = 0; imod < fNmodules; imod++) {
1286 if (ist == 1 && imod > 3 )
continue;
1287 if (ist == 2 && imod > 1 )
continue;
1289 for (Int_t cl = 0; cl < NhitsXYmod_[ist][imod]; cl++) {
1290 if(fDebug) cout<<XpspCoord_[ist][imod][cl]<<
" ";
1293 if(fDebug) cout<<endl;
1295 if(fDebug) cout<<endl;
1297 for(Int_t ist = 1; ist < fNstations; ist++) {
1298 if(fDebug) cout<<
" Yspatial["<<ist<<
"] ";
1299 for (Int_t imod = 0; imod < fNmodules; imod++) {
1300 if (ist == 1 && imod > 3 )
continue;
1301 if (ist == 2 && imod > 1 )
continue;
1302 for (Int_t cl = 0; cl < NhitsXYmod_[ist][imod]; cl++) {
1303 if(fDebug) cout<<YspCoord_[ist][imod][cl]<<
" ";
1306 if(fDebug) cout<<endl;
1308 if(fDebug) cout<<endl;
1310 for (Int_t ist = 1; ist < fNstations ; ist++) {
1311 for (Int_t imod = 0; imod < fNmodules; imod++) {
1312 if (ist == 1 && imod > 3 )
continue;
1313 if (ist == 2 && imod > 1 )
continue;
1314 for (Int_t cl = 0; cl < NhitsXYmod_[ist][imod]; cl++) {
1315 if(fDebug && ist == 1) hY_XSisp1bef->Fill(XspCoord_[ist][imod][cl],YspCoord_[ist][imod][cl]);
1316 if(fDebug && ist == 2) hY_XSisp2bef->Fill(XspCoord_[ist][imod][cl],YspCoord_[ist][imod][cl]);
1327void BmnSiliconTrackFinder::CountSpatialPoints(Int_t **Nleftoversp_, Int_t *Nsp_st_){
1328 if(fDebug) cout<<
" CountSpatialPoints "<<endl;
1330 for (Int_t st = 1; st < fNstations; ++st) {
1332 for (Int_t imod = 0; imod < fNmodules; imod++) {
1334 Nsp_st_[st] += Nleftoversp_[st][imod];
1337 if ( fDebug) cout<<
" next case? Nsp_st[1] "<<Nsp_st[1]<<
" Nsp_st[2] "<<Nsp_st[2]<<
" Nsp_st[3] "<<Nsp_st[3]<<endl;
1343Bool_t BmnSiliconTrackFinder::SkipEvent(Int_t **NclustX_, Int_t **NclustXp_, Double_t ***XCoord_, Double_t ***XpCoord_){
1346 Int_t NfithitsX[fNstations], NfithitsXp[fNstations];
1347 for (Int_t ist = 0; ist < fNstations; ist++) {
1349 NfithitsXp[ist] = 0;
1352 for (Int_t ist = 1; ist < fNstations; ist++) {
1353 for (Int_t imod = 0; imod < fNmodules; imod++) {
1354 if (ist == 1 && imod > 3 )
continue;
1355 if (ist == 2 && imod > 1 )
continue;
1357 for (Int_t cl = 0; cl < NclustX_[ist][imod]; cl++) {
1361 for (Int_t cl = 0; cl < NclustXp_[ist][imod]; cl++) {
1368 cout<<
" Number of hits: "<<endl;
1369 for (Int_t ist = 1; ist < fNstations; ist++) {
1370 cout<<
" NfithitsX["<<ist<<
"] "<<NfithitsX[ist]<<
" NfithitsXp["<<ist<<
"] "<<NfithitsXp[ist]<<endl;
1376 if ( (NfithitsX[1] + NfithitsX[2] + NfithitsX[3]) > 39 && (NfithitsXp[1] + NfithitsXp[2] + NfithitsXp[3]) > 36 ){
1377 if(fDebug) cout<<
" --Skip event!!!!! with high Si-hits multiplicity --"<<endl;
1380 if ( (NfithitsX[1] + NfithitsX[2] + NfithitsX[3]) > 51 ){
1381 if(fDebug) cout<<
" --Skip event!!!!! with high Si-hits multiplicity --"<<endl;
1391void BmnSiliconTrackFinder::RecordingTracksAfterSpatialPoints(vector<tracksX> & vec_tracks_, vector<tracksX> & CleanTr_){
1393 if(fDebug) cout<<endl;
1394 if(fDebug) cout<<
" RecordingTracksAfterSpatialPoints size:"<<vec_tracks_.size()<<endl;
1395 for(
size_t itr = 0; itr< vec_tracks_.size(); ++itr) {
1396 CleanTr_.push_back(vec_tracks_.at(itr));
1398 if (fDebug) cout<<
" itr "<<itr<<
" Chi2 "<<CleanTr_.at(itr).Chi2<<
" Nhits "<<CleanTr_.at(itr).Nhits<<
" Xv "<<CleanTr_.at(itr).Xv<<
" Yv "<<CleanTr_.at(itr).Yv<<endl;
1399 for (Int_t st = 1; st < fNstations; ++st) {
1400 if (fDebug) cout<<
" CleanTr_.at("<<itr<<
").CoordX["<<st<<
"] "<<CleanTr_.at(itr).CoordX[st]<<
" Xp "<<CleanTr_.at(itr).CoordXp[st]<<endl;
1403 if(fDebug) cout<<
" CleanTr_.size() "<<CleanTr_.size()<<endl;
1404 if(fDebug) cout<<endl;
1408 if (fDebug) cout <<
" sorting vector"<<endl;
1409 for (
size_t itr = 0; itr < CleanTr_.size(); ++itr) {
1410 if (fDebug) cout <<
" itr " << itr <<
" Chi2 " << CleanTr_.at(itr).Chi2 <<endl;
1412 for (
size_t itr = 0; itr < CleanTr_.size(); ++itr) {
1413 for (Int_t st = 1; st < fNstations; ++st) {
1414 if (fDebug) cout<<
" st "<<st<<
" Z "<<CleanTr_.at(itr).CoordZ[st]<<endl;
1415 if (fDebug) cout<<
" CleanTr.CoordX["<<st <<
"] "<<CleanTr_.at(itr).CoordX[st]<<
" strip "<<CleanTr_.at(itr).StripNumX[st]<<
" Clsize "<<CleanTr_.at(itr).ClSizeX[st]<<
" Q "<<CleanTr_.at(itr).SumQX[st]<<endl;
1416 if (fDebug) cout<<
" CleanTr.CoordXp["<<st <<
"] "<<CleanTr_.at(itr).CoordXp[st]<<
" strip "<<CleanTr_.at(itr).StripNumXp[st]<<
" Clsize "<<CleanTr_.at(itr).ClSizeXp[st]<<
" Q "<<CleanTr_.at(itr).SumQXp[st]<<endl;
1419 vec_tracks_.clear();
1425void BmnSiliconTrackFinder::RecordingTracks(vector<tracksX> & vec_tracks_, vector<tracksX> & CleanTr_){
1426 if(fDebug) cout<<endl;
1427 if(fDebug) cout<<
" RecordingTracks size:"<<vec_tracks_.size()<<endl;
1428 if(fDebug) cout<<
" CleanTr size: "<<CleanTr_.size()<<endl;
1429 Double_t ZposC[fNstations];
1431 ZposC[1] = 213.999 + 0.03 -6.9;
1432 ZposC[2] = 213.999 + 0.03 -6.9;
1433 ZposC[3] = 327.130 + 0.03;
1435 for (Int_t is = 1; is < fNstations; is++) {
1436 ZposC[is] += Z0_SRC_target;
1438 for(
size_t itr = 0; itr< vec_tracks_.size(); ++itr) {
1441 Double_t X1_cand = vec_tracks_.at(itr).param[0] * (ZposC[1] - Zcentr) + vec_tracks_.at(itr).param[1];
1442 Double_t Y1_cand = vec_tracks_.at(itr).param[2] * (ZposC[1] - Zcentr) + vec_tracks_.at(itr).param[3];
1444 Double_t X3_cand = vec_tracks_.at(itr).param[0] * (ZposC[3] - Zcentr) + vec_tracks_.at(itr).param[1];
1445 Double_t Y3_cand = vec_tracks_.at(itr).param[2] * (ZposC[3] - Zcentr) + vec_tracks_.at(itr).param[3];
1448 Bool_t Cross_6p = 0;
1450 for(
size_t InIter = 0; InIter < CleanTr_.size(); ++InIter) {
1452 Double_t X1_tr = CleanTr_.at(InIter).param[0] * (ZposC[1] - Zcentr) + CleanTr_.at(InIter).param[1];
1453 Double_t Y1_tr = CleanTr_.at(InIter).param[2] * (ZposC[1] - Zcentr) + CleanTr_.at(InIter).param[3];
1455 Double_t X3_tr = CleanTr_.at(InIter).param[0] * (ZposC[3] - Zcentr) + CleanTr_.at(InIter).param[1];
1456 Double_t Y3_tr = CleanTr_.at(InIter).param[2] * (ZposC[3] - Zcentr) + CleanTr_.at(InIter).param[3];
1459 if ( X1_cand < X1_tr && X3_cand > X3_tr) {
1463 if ( Y1_cand < Y1_tr && Y3_cand > Y3_tr) {
1467 if ( X1_cand > X1_tr && X3_cand < X3_tr) {
1471 if ( Y1_cand > Y1_tr && Y3_cand < Y3_tr) {
1478 if (fDebug) cout<<
" itr "<<itr<<
" CleanTr size: "<<CleanTr_.size()<<endl;
1479 CleanTr_.push_back(vec_tracks_.at(itr));
1481 if (fDebug) cout<<
" itr "<<itr<<
" Chi2 "<<vec_tracks_.at(itr).Chi2<<
" Nhits "<<vec_tracks_.at(itr).Nhits<<
" Xv "<<vec_tracks_.at(itr).Xv<<
" Yv "<<vec_tracks_.at(itr).Yv<<endl;
1482 for (Int_t st = 1; st < fNstations; ++st) {
1484 cout<<
" st "<<st<<
" Z "<<vec_tracks_.at(itr).CoordZ[st]<<endl;
1485 cout<<
" CleanTr.CoordX["<<st <<
"] "<<vec_tracks_.at(itr).CoordX[st]<<
" strip "<<vec_tracks_.at(itr).StripNumX[st]<<
" Clsize "<<vec_tracks_.at(itr).ClSizeX[st]<<
" Q "<<vec_tracks_.at(itr).SumQX[st]<<endl;
1486 cout<<
" CleanTr.CoordXp["<<st <<
"] "<<vec_tracks_.at(itr).CoordXp[st]<<
" strip "<<vec_tracks_.at(itr).StripNumXp[st]<<
" Clsize "<<vec_tracks_.at(itr).ClSizeXp[st]<<
" Q "<<vec_tracks_.at(itr).SumQXp[st]<<endl;
1491 if(fDebug) cout<<endl;
1492 if(fDebug) cout<<
" now CleanTr.size() "<<CleanTr_.size()<<endl;
1493 if(fDebug) cout<<endl;
1494 vec_tracks_.clear();
1500void BmnSiliconTrackFinder::RecordingTracks_case3(vector<tracksX> & vec_tracks_, vector<tracksX> & CleanTr_){
1501 if(fDebug) cout<<endl;
1502 if(fDebug) cout<<
" RecordingTracks size:"<<vec_tracks_.size()<<endl;
1503 Double_t ZposC[fNstations];
1505 ZposC[1] = 213.999 + 0.03 -6.9;
1506 ZposC[2] = 213.999 + 0.03 -6.9;
1507 ZposC[3] = 327.130 + 0.03;
1509 if (fDebug) cout<<
" after case 3 write cand: "<<vec_tracks_.size()<<endl;
1512 Double_t sigma_y = 0.1;
1513 for(
size_t itr = 0; itr< vec_tracks_.size(); ++itr){
1515 Double_t X1_cand = vec_tracks_.at(itr).param[0] * (ZposC[1] - Zcentr) + vec_tracks_.at(itr).param[1];
1516 Double_t Y1_cand = vec_tracks_.at(itr).param[2] * (ZposC[1] - Zcentr) + vec_tracks_.at(itr).param[3];
1518 Double_t X3_cand = vec_tracks_.at(itr).param[0] * (ZposC[3] - Zcentr) + vec_tracks_.at(itr).param[1];
1519 Double_t Y3_cand = vec_tracks_.at(itr).param[2] * (ZposC[3] - Zcentr) + vec_tracks_.at(itr).param[3];
1524 Bool_t Cross_6p = 0, rej_cand = 0;
1526 for(
size_t InIter = 0; InIter < CleanTr_.size(); ++InIter){
1528 Double_t X1_tr = CleanTr_.at(InIter).param[0] * (ZposC[1] - Zcentr) + CleanTr_.at(InIter).param[1];
1529 Double_t Y1_tr = CleanTr_.at(InIter).param[2] * (ZposC[1] - Zcentr) + CleanTr_.at(InIter).param[3];
1531 Double_t X3_tr = CleanTr_.at(InIter).param[0] * (ZposC[3] - Zcentr) + CleanTr_.at(InIter).param[1];
1532 Double_t Y3_tr = CleanTr_.at(InIter).param[2] * (ZposC[3] - Zcentr) + CleanTr_.at(InIter).param[3];
1537 if ( X1_cand < X1_tr && X3_cand > X3_tr) {Cross_6p = 1; }
1538 if ( Y1_cand < Y1_tr && (Y3_cand - sigma_y) > Y3_tr) {Cross_6p = 1; }
1540 if ( X1_cand > X1_tr && X3_cand < X3_tr) {Cross_6p = 1; }
1541 if ( Y1_cand > Y1_tr && (Y3_cand + sigma_y) < Y3_tr) {Cross_6p = 1; }
1545 if ( CleanTr_.at(InIter).Nhits > vec_tracks_.at(itr).Nhits ) {
1550 if ( !rej_cand && CleanTr_.at(InIter).Nhits < vec_tracks_.at(itr).Nhits ){
1553 CleanTr_.erase(CleanTr_.begin()+ InIter);
1556 if( !rej_cand && Cross_6p ){
1557 if (CleanTr_.at(InIter).Nhits == vec_tracks_.at(itr).Nhits && CleanTr_.at(InIter).Nhits == 4 ){
1563 if ( CleanTr_.at(InIter).Nhits == vec_tracks_.at(itr).Nhits && CleanTr_.at(InIter).Nhits > 4 ){
1564 if ( CleanTr_.at(InIter).Chi2 < vec_tracks_.at(itr).Chi2 ) rej_cand = 1;
1565 if ( CleanTr_.at(InIter).Chi2 > vec_tracks_.at(itr).Chi2 ) {
1566 CleanTr_.erase(CleanTr_.begin()+ InIter);
1579 if ( !Cross_6p && !rej_cand ) {
1580 if (fDebug) cout<<
" write good itr "<<itr<<endl;
1581 CleanTr_.push_back(vec_tracks_.at(itr));
1582 if (fDebug) cout<<
" push good itr "<<itr<<endl;
1584 if (fDebug) cout<<
" itr "<<itr<<
" Chi2 "<<vec_tracks_.at(itr).Chi2<<
" Nhits "<<vec_tracks_.at(itr).Nhits<<
" Xv "<<vec_tracks_.at(itr).Xv<<
" Yv "<<vec_tracks_.at(itr).Xv<<endl;
1585 for (Int_t st = 1; st < fNstations; ++st) {
1586 if (fDebug) cout<<
" st "<<st<<
" Z "<<CleanTr_.at(itr).CoordZ[st]<<endl;
1587 if (fDebug) cout<<
" CleanTr.CoordX["<<st <<
"] "<<vec_tracks_.at(itr).CoordX[st]<<
" strip "<<vec_tracks_.at(itr).StripNumX[st]<<
" Clsize "<<vec_tracks_.at(itr).ClSizeX[st]<<
" Q "<<vec_tracks_.at(itr).SumQX[st]<<endl;
1588 if (fDebug) cout<<
" CleanTr.CoordXp["<<st <<
"] "<<vec_tracks_.at(itr).CoordXp[st]<<
" strip "<<vec_tracks_.at(itr).StripNumXp[st]<<
" Clsize "<<vec_tracks_.at(itr).ClSizeXp[st]<<
" Q "<<vec_tracks_.at(itr).SumQXp[st]<<endl;
1593 if(fDebug) cout<<endl;
1594 if(fDebug) cout<<
" now CleanTr_.size() "<<CleanTr_.size()<<endl;
1595 if(fDebug) cout<<endl;
1596 vec_tracks_.clear();
1604void BmnSiliconTrackFinder::Case1(Int_t **NhitsXYmod_, Double_t ***XspCoord_, Double_t ***XpspCoord_, Double_t ***YspCoord_, Double_t ***SigmspX_, Double_t ***SigmspXp_, Double_t ***SigmspY_,
1605 vector<tracksX> & out_tracks,
1606 Double_t ***XspClust_, Double_t ***XspClust_width_, Double_t ***sumQxsp_, Double_t ***XpspClust_, Double_t ***XpspClust_width_, Double_t ***sumQxpsp_){
1608 Int_t NpX_in_line = 0;
1609 Double_t dminx1_x2 = 999.;
1615 for ( Int_t
i = 1;
i < fNstations; ++
i){
1620 Double_t bx_candidate = 999., by_candidate = 999.;
1621 Int_t iClustXMod_candidate = -1, iModX_candidate = -1;
1622 Double_t dxV_thr = .2, dyV_thr = .4;
1624 if(fDebug) cout<<endl;
1626 if(fDebug) cout<<
" case1: Spatial track "<<endl;
1628 for (Int_t imod1 = 0; imod1 < fNmodules1; imod1++) {
1629 if(fDebug && NhitsXYmod_[1][imod1] > 0 ) cout<<
" NhitsXYmod_[1]["<<imod1<<
"] "<<NhitsXYmod_[1][imod1]<<endl;
1630 for (Int_t cl1 = 0; cl1 < NhitsXYmod_[1][imod1]; cl1++){
1633 for (Int_t imod2 = 0; imod2 < fNmodules2; imod2++) {
1634 if(fDebug && NhitsXYmod_[2][imod2]) cout<<
" NhitsXYmod_[2]["<<imod2<<
"] "<<NhitsXYmod_[2][imod2]<<endl;
1635 for (Int_t cl2 = 0; cl2 < NhitsXYmod_[2][imod2]; cl2++){
1645 if (
fabs( XspCoord_[1][imod1][cl1] - XspCoord_[2][imod2][cl2] ) > half_roadX1_X2)
continue;
1646 if (
fabs( XspCoord_[1][imod1][cl1] - XpspCoord_[1][imod1][cl1] ) > half_roadX1_Xp1)
continue;
1647 if (
fabs( XspCoord_[2][imod2][cl2] - XpspCoord_[2][imod2][cl2] ) > half_roadX2_Xp2)
continue;
1648 if (
fabs( YspCoord_[1][imod1][cl1] - YspCoord_[2][imod2][cl2] ) > half_roadY1_Y2)
continue;
1653 if ( dminx1_x2 >
fabs( XspCoord_[1][imod1][cl1] - XspCoord_[2][imod2][cl2] ) ){
1654 dminx1_x2 =
fabs( XspCoord_[1][imod1][cl1] - XspCoord_[2][imod2][cl2] );
1655 iClustXMod[1] = cl1;
1656 iClustXMod[2] = cl2;
1662 if ( NpX_in_line < 2)
continue;
1664 Bool_t was_candidate = 0;
1665 bx_candidate = 999., by_candidate = 999.;
1669 for (Int_t imod3 = 0; imod3 < fNmodules3; imod3++) {
1670 if(fDebug && NhitsXYmod_[3][imod3] > 0 ) cout<<
" NhitsXYmod_[3]["<<imod3<<
"] "<<NhitsXYmod_[3][imod3]<<endl;
1671 for (Int_t cl3 = 0; cl3 < NhitsXYmod_[3][imod3]; cl3++){
1677 if (
fabs( XspCoord_[3][imod3][cl3] - XpspCoord_[3][imod3][cl3] ) > half_roadX3_Xp3)
continue;
1690 Double_t ax = (XspCoord_[3][imod3][cl3] - XspCoord_[1][imod1][cl1] ) / (Zstation[3][imod3] - Zstation[1][imod1]);
1691 Double_t bx_target_region = ax * ( Z0_SRC_target - Zstation[3][imod3] ) + XspCoord_[3][imod3][cl3];
1693 Double_t ay = (YspCoord_[3][imod3][cl3] - YspCoord_[1][imod1][cl1] ) / (Zstation[3][imod3] - Zstation[1][imod1]);
1694 Double_t by_target_region = ay * ( Z0_SRC_target - Zstation[3][imod3] ) + YspCoord_[3][imod3][cl3];
1696 if(fDebug) cout<<
" X1 "<< XspCoord_[1][imod1][cl1]<<
" X2 "<< XspCoord_[2][imod2][cl2]<<
" X3 "<< XspCoord_[3][imod3][cl3]<<
" was_candidate "<<was_candidate<<
1697 " ax "<<ax<<
" ay "<<ay<<
" Z1 "<<Zstation[1][imod1]<<
" Z3 "<<Zstation[3][imod1]<<
1702 if(fDebug) cout<<
" bx_target_region "<<bx_target_region<<
" half_target_regionX "<<half_target_regionX<<endl;
1703 if (
fabs(kX_target - bx_target_region) > half_target_regionX)
continue;
1704 if(fDebug) cout<<
" by_target_region "<<by_target_region<<
" half_target_regionY "<<half_target_regionY<<endl;
1705 if (
fabs(kY_target - by_target_region) > half_target_regionY)
continue;
1707 if ( was_candidate == 0 ) {
1708 bx_candidate = bx_target_region;
1709 by_candidate = by_target_region;
1710 iModX_candidate = imod3;
1711 iClustXMod_candidate = cl3;
1713 if (fDebug) cout<<
" write candidate: bx_target_region "<<bx_target_region<<
" bx_candidate "<<bx_candidate<<
" by_target_region "<<by_target_region<<
" by_candidate "<<by_candidate<<endl;
1716 if ( iModX_candidate == imod3 && iClustXMod_candidate == cl3 ) {
if (fDebug) cout<<
" 1"<<endl;
continue;}
1717 if (
fabs(kX_target - bx_target_region) > (
fabs(kX_target - bx_candidate) + dxV_thr) ) {
if (fDebug) cout<<
" 2"<<endl;
continue;}
1719 if (
fabs(kY_target - by_target_region) > (
fabs(kY_target - by_candidate) + dyV_thr) ) {
if (fDebug) cout<<
" 3"<<endl;
continue;}
1723 if ((
fabs(kX_target - bx_target_region) + 0.25*
fabs(kY_target - by_target_region) ) > (
fabs(kX_target - bx_candidate) + 0.25*
fabs(kY_target - by_candidate)) )
1724 {
if (fDebug) cout<<
" 4"<<endl;
continue;}
1727 bx_candidate = bx_target_region;
1728 by_candidate = by_target_region;
1729 iModX_candidate = imod3;
1730 iClustXMod_candidate = cl3;
1731 if (fDebug) cout<<
" re - write "<<endl;
1737 if (bx_candidate > 900.)
continue;
1739 if ( tmpTr.
CoordX[3] == XspCoord_[3][iModX_candidate][iClustXMod_candidate] )
continue;
1740 if ( tmpTr.
CoordXp[3] == XpspCoord_[3][iModX_candidate][iClustXMod_candidate] )
continue;
1743 if ( iModX_candidate > -1){
1744 iModX[3] = iModX_candidate;
1745 iClustXMod[3] = iClustXMod_candidate;
1752 for (Int_t ist = 1; ist < fNstations; ++ist) {
1753 for (Int_t imod = 0; imod < fNmodules; imod++) {
1754 Xforglfit[ist][imod] = -999.;
1755 Xpforglfit[ist][imod] = -999.;
1756 SigmXforglfit[ist][imod] = -999.;
1757 SigmXpforglfit[ist][imod] = -999.;
1761 for (Int_t idx = 1; idx < fNstations; ++idx) {
1762 if ( iClustXMod[idx] < 0 )
continue;
1765 Xforglfit[idx][iModX[idx]] = XspCoord_[idx][iModX[idx]][iClustXMod[idx]];
1766 SigmXforglfit[idx][iModX[idx]] = SigmspX_[idx][iModX[idx]][iClustXMod[idx]];
1768 Xpforglfit[idx][iModX[idx]] = XpspCoord_[idx][iModX[idx]][iClustXMod[idx]];
1769 SigmXpforglfit[idx][iModX[idx]]= SigmspXp_[idx][iModX[idx]][iClustXMod[idx]];
1771 Yforglfit[idx][iModX[idx]] = YspCoord_[idx][iModX[idx]][iClustXMod[idx]];
1772 SigmYforglfit[idx][iModX[idx]] = SigmspY_[idx][iModX[idx]][iClustXMod[idx]];
1776 Int_t TotalNumberOfHits = 0;
1777 Double_t ChiSquareNdf = 0.;
1785 calculationOfChi2(Xforglfit, Xpforglfit, SigmXforglfit, SigmXpforglfit, iClustXMod, iModX, ChiSquare, line);
1787 TotalNumberOfHits = 6;
1788 ChiSquareNdf = ChiSquare/(TotalNumberOfHits - 4);
1790 if ( fDebug) cout<<
" before cut ChiSquareNDF "<<ChiSquareNdf<<
" Chi2_Globcut "<<Chi2_Globcut<<
" TotalNumberOfHits "<<TotalNumberOfHits <<endl;
1791 if (fDebug) cout<<endl;
1793 if ( ChiSquareNdf > Chi2_Globcut )
continue;
1801 if (TotalNumberOfHits == 6 ){
1803 for(
size_t itr = 0; itr< out_tracks.size(); ++itr){
1806 (out_tracks.at(itr).CoordX[1] == Xforglfit[1][iModX[1]] && out_tracks.at(itr).CoordXp[1] == Xpforglfit[1][iModX[1]] )
1808 (out_tracks.at(itr).CoordX[2] == Xforglfit[2][iModX[2]] && out_tracks.at(itr).CoordXp[2] == Xpforglfit[2][iModX[2]] )
1810 (out_tracks.at(itr).CoordX[3] == Xforglfit[3][iModX[3]] && out_tracks.at(itr).CoordXp[3] == Xpforglfit[3][iModX[3]] )){
1811 if ( out_tracks.at(itr).Chi2 > ChiSquareNdf) out_tracks.erase(out_tracks.begin()+ itr);
1812 else ChiSquareNdf = 999.;
1816 if ( ChiSquareNdf > Chi2_Globcut )
continue;
1818 tmpTr.
Chi2 = ChiSquareNdf ;
1819 tmpTr.
Nhits = TotalNumberOfHits;
1820 tmpTr.
Xv = bx_candidate;
1821 tmpTr.
Yv = by_candidate;
1822 tmpTr.
param[0] = line[0];
1823 tmpTr.
param[1] = line[1];
1824 tmpTr.
param[2] = line[2];
1825 tmpTr.
param[3] = line[3];
1827 if (fDebug) cout<<
" Xv_candidate "<<bx_candidate<<
" Yv_candidate "<<by_candidate<<endl;
1828 for (Int_t st = 1; st < fNstations; ++st) {
1829 if (iModX[st] < 0)
continue;
1831 tmpTr.
ModNum[st] = iModX[st];
1832 tmpTr.
CoordX[st] = Xforglfit[st][iModX[st]];
1833 tmpTr.
CoordY[st] = Yforglfit[st][iModX[st]];
1834 tmpTr.
CoordXp[st] = Xpforglfit[st][iModX[st]];
1835 tmpTr.
SigmaX[st] = SigmXforglfit[st][iModX[st]];
1836 tmpTr.
SigmaY[st] = SigmYforglfit[st][iModX[st]];
1837 tmpTr.
SigmaXp[st] = SigmXpforglfit[st][iModX[st]];
1838 tmpTr.
CoordZ[st] = ZstnCforglfit[st][iModX[st]];
1840 tmpTr.
StripNumX[st] = XspClust_[st][iModX[st]][iClustXMod[st]];
1841 tmpTr.
StripNumXp[st]= XpspClust_[st][iModX[st]][iClustXMod[st]];
1842 tmpTr.
ClSizeX[st] = XspClust_width_[st][iModX[st]][iClustXMod[st]];
1843 tmpTr.
ClSizeXp[st] = XpspClust_width_[st][iModX[st]][iClustXMod[st]];
1844 tmpTr.
SumQX[st] = sumQxsp_[st][iModX[st]][iClustXMod[st]];
1845 tmpTr.
SumQXp[st] = sumQxpsp_[st][iModX[st]][iClustXMod[st]];
1847 if (fDebug) cout<<
" st "<<st<<
" Z "<<tmpTr.
CoordZ[st]<<endl;
1848 if (fDebug) cout<<
" CoordX["<<st <<
"] "<<tmpTr.
CoordX[st]<<
" strip "<<tmpTr.
StripNumX[st]<<
" Clsize "<<tmpTr.
ClSizeX[st]<<
" Q "<<tmpTr.
SumQX[st]<<endl;
1849 if (fDebug) cout<<
" CoordXp["<<st <<
"] "<<tmpTr.
CoordXp[st]<<
" strip "<<tmpTr.
StripNumXp[st]<<
" Clsize "<<tmpTr.
ClSizeXp[st]<<
" Q "<<tmpTr.
SumQXp[st]<<endl;
1852 out_tracks.push_back(tmpTr);
1854 if (fDebug) cout<<
" ax "<<tmpTr.
param[0]<<
" bx "<<tmpTr.
param[1]<<
" ay "<<tmpTr.
param[2]<<
" by "<<tmpTr.
param[3]<<
" Chi2 "<<tmpTr.
Chi2<<endl;
1855 if (fDebug) cout<<
" write "<<out_tracks.size()<<endl;
1867void BmnSiliconTrackFinder::CheckPoints(Int_t **NclustX_, Int_t **NclustXp_,
1868 Double_t ***XCoord_, Double_t ***XpCoord_, Double_t ***SigmaX_, Double_t ***SigmaXp_,Double_t ***XspCoord_, Double_t ***XpspCoord_, Double_t ***YspCoord_, Double_t ***SigmspX_, Double_t ***SigmspXp_, Double_t ***SigmspY_,
1869 Int_t **NhitsXYmod_, vector<tracksX> & tracks,Int_t **Nspatial, Int_t **NpointsX,Int_t **NpointsXp,
1870 Double_t ***pointsXsp, Double_t ***pointsXpsp,Double_t ***pointsYsp,Double_t ***pointsXsigsp,Double_t ***pointsXpsigsp, Double_t ***pointsYsigsp,
1871 Double_t ***pointsX, Double_t ***pointsXp, Double_t ***pointsXsig,Double_t ***pointsXpsig){
1873 if(fDebug) cout<<
" CheckPoints "<<endl;
1875 for (Int_t st = 0; st < fNstations; ++st) {
1876 for (Int_t imod = 0; imod < fNmodules; imod++) {
1877 Nspatial[st][imod] = 0;
1878 NpointsX[st][imod] = 0;
1879 NpointsXp[st][imod] = 0;
1880 for (Int_t j = 0; j < kBig; j++) {
1881 pointsXsp[st][imod][j] = -999.;
1882 pointsXpsp[st][imod][j] = -999.;
1883 pointsYsp[st][imod][j] = -999.;
1884 pointsXsigsp[st][imod][j] = -999.;
1885 pointsXpsigsp[st][imod][j] = -999.;
1886 pointsYsigsp[st][imod][j] = -999.;
1887 pointsX[st][imod][j] = -999.;
1888 pointsXp[st][imod][j] = -999.;
1889 pointsXsig[st][imod][j] = -999.;
1890 pointsXpsig[st][imod][j] = -999.;
1895 for (Int_t st = 1; st < fNstations; ++st){
1896 for (Int_t imod = 0; imod < fNmodules; imod++) {
1898 for (Int_t cl = 0; cl < NhitsXYmod_[st][imod]; cl++){
1900 for (
size_t itr = 0 ; itr < tracks.size(); ++itr){
1902 if ( XspCoord[st][imod][cl] == tracks.at(itr).CoordX[st]){
1906 if ( XpspCoord[st][imod][cl] == tracks.at(itr).CoordXp[st]){
1911 if ( clux > -1)
continue;
1912 pointsXsp[st][imod][Nspatial[st][imod]] = XspCoord[st][imod][cl] ;
1913 pointsXpsp[st][imod][Nspatial[st][imod]] = XpspCoord[st][imod][cl] ;
1914 pointsYsp[st][imod][Nspatial[st][imod]] = YspCoord[st][imod][cl] ;
1915 pointsXsigsp[st][imod][Nspatial[st][imod]] = SigmspX[st][imod][cl];
1916 pointsXpsigsp[st][imod][Nspatial[st][imod]]= SigmspXp[st][imod][cl];
1917 pointsYsigsp[st][imod][Nspatial[st][imod]] = SigmspY[st][imod][cl];
1918 Nspatial[st][imod]++;
1924 for (Int_t st = 1; st < fNstations; ++st) {
1925 for (Int_t imod = 0; imod < fNmodules; imod++) {
1926 for (Int_t cl = 0; cl < NclustX[st][imod]; cl++){
1928 for (
size_t itr = 0 ; itr < tracks.size(); ++itr){
1929 if ( XCoord[st][imod][cl] == tracks.at(itr).CoordX[st]){
1934 if ( clux > -1)
continue;
1935 pointsX[st][imod][NpointsX[st][imod]] = XCoord[st][imod][cl] ;
1936 pointsXsig[st][imod][NpointsX[st][imod]]= SigmaX[st][imod][cl] ;
1937 NpointsX[st][imod]++;
1940 for (Int_t cl = 0; cl < NclustXp[st][imod]; cl++){
1942 for (
size_t itr = 0 ; itr < tracks.size(); ++itr){
1943 if ( XpCoord[st][imod][cl] == tracks.at(itr).CoordXp[st]){
1948 if ( clux > -1)
continue;
1949 pointsXp[st][imod][NpointsXp[st][imod]] = XpCoord[st][imod][cl] ;
1950 pointsXpsig[st][imod][NpointsX[st][imod]] = SigmaXp[st][imod][cl] ;
1951 NpointsXp[st][imod]++;
1956 if(fDebug) cout<<
" after rejection "<<endl;
1958 for(Int_t ist = 1; ist < fNstations; ist++){
1959 if(fDebug) cout<<
" pointsXsp["<<ist<<
"] ";
1960 for (Int_t imod = 0; imod < fNmodules; imod++) {
1961 for (Int_t cl = 0; cl < Nspatial[ist][imod]; cl++){
1962 if(fDebug) cout<<pointsXsp[ist][imod][cl]<<
" ";
1965 if(fDebug) cout<<endl;
1967 if(fDebug) cout<<endl;
1969 for(Int_t ist = 1; ist < fNstations; ist++) {
1970 if(fDebug) cout<<
" pointsX["<<ist<<
"] ";
1971 for (Int_t imod = 0; imod < fNmodules; imod++) {
1972 for (Int_t cl = 0; cl < NpointsX[ist][imod]; cl++){
1973 if(fDebug) cout<<pointsX[ist][imod][cl]<<
" ";
1976 if(fDebug) cout<<endl;
1978 if(fDebug) cout<<endl;
1980 for(Int_t ist = 1; ist < fNstations; ist++) {
1981 if(fDebug) cout<<
" pointsXp["<<ist<<
"] ";
1982 for (Int_t imod = 0; imod < fNmodules; imod++) {
1983 for (Int_t cl = 0; cl < NpointsXp[ist][imod]; cl++){
1984 if(fDebug) cout<<pointsXp[ist][imod][cl]<<
" ";
1987 if(fDebug) cout<<endl;
1989 if(fDebug) cout<<endl;
1995void BmnSiliconTrackFinder::CheckLeftover(Int_t **Nspatial, Int_t **NpointsX,Int_t **NpointsXp,
1996 Double_t ***pointsXsp, Double_t ***pointsXpsp,Double_t ***pointsYsp,
1997 Double_t ***pointsXsigsp,Double_t ***pointsXpsigsp, Double_t ***pointsYsigsp,
1998 Double_t ***pointsX, Double_t ***pointsXp, Double_t ***pointsXsig,
1999 Double_t ***pointsXpsig, vector<tracksX> & CleanTr_,
2000 Int_t **NhitsXYmod_, Double_t ***XspCoord_, Double_t ***YspCoord_, Double_t ***SigmspX_, Double_t ***SigmspY_){
2002 if(fDebug) cout<<
" CheckLeftover"<<endl;
2003 bool woSp[fNstations];
2004 for (Int_t st = 0; st < fNstations ; ++st) {
2008 if(fDebug) cout<<endl;
2009 if (fDebug && CleanTr_.size() > 0) cout<<
"track was but: "<<endl;
2011 for (
size_t InIter = 0; InIter < CleanTr_.size(); ++InIter){
2013 if (CleanTr_.at(InIter).Nhits < 5 || CleanTr_.at(InIter).Nhits == 6)
continue;
2015 if(fDebug) cout<<
" itr "<<InIter<<
" if Nhits = 5: "<<CleanTr_.at(InIter).Nhits<<endl;
2016 for (Int_t st = 1; st < fNstations -1 ; ++st) {
2020 if ( CleanTr_.at(InIter).CoordX[st] < -900.){
2022 if (fDebug) cout<<
"no X["<<st<<
"]: "<< CleanTr_.at(InIter).CoordX[st]<<endl;
2023 if (fDebug) cout<<
" was sp points? "<<endl;
2024 for (Int_t imod = 0; imod < fNmodules; imod++) {
2025 if ( st == 1 && imod > 3)
continue;
2026 if ( st == 2 && imod > 1)
continue;
2028 if (fDebug) cout<<
" Nspatial["<<st<<
"]= "<<NhitsXYmod_[st][imod]<<endl;
2030 if (NhitsXYmod_[st][imod] > 0 ) woSp[st] = 1;
2032 for (Int_t cl = 0; cl < NhitsXYmod_[st][imod]; cl++){
2033 if(fDebug) cout<<
"st["<<st<<
"] "<<
" X "<<XspCoord_[st][imod][cl]<<endl;
2037 Double_t FitX = CleanTr_.at(InIter).param[0] * CleanTr_.at(InIter).CoordZ[st] + CleanTr_.at(InIter).param[1];
2038 Double_t FitY = CleanTr_.at(InIter).param[2] * CleanTr_.at(InIter).CoordZ[st] + CleanTr_.at(InIter).param[3];
2041 if ( woSp[st] == 0){
2042 if(fDebug) cout<<
" woSp["<<st<<
"] "<<woSp[st]<<
"-> FitX "<<FitX<<
" FitY "<<FitY<<endl;
2043 XspCoord_[st][1][NhitsXYmod_[st][1]] = FitX;
2044 YspCoord_[st][1][NhitsXYmod_[st][1]] = FitY;
2045 SigmspX_[st][1][NhitsXYmod_[st][1]] = 0.0060;
2046 SigmspY_[st][1][NhitsXYmod_[st][1]] = 0.2;
2047 NhitsXYmod_[st][1]++;
2052 if ( CleanTr_.at(InIter).CoordY[st] < -900.){
2053 if (fDebug) cout<<
"no Y["<<st<<
"] "<< CleanTr_.at(InIter).CoordY[st]<<endl;
2054 for (Int_t imod = 0; imod < fNmodules; imod++) {
2055 if ( st == 1 && imod > 3)
continue;
2056 if ( st == 2 && imod > 1)
continue;
2058 if (fDebug) cout<<
" Nspatial = "<<NhitsXYmod_[st][imod]<<endl;
2060 if (NhitsXYmod_[st][imod] > 0 ) woSp[st] = 1;
2062 for (Int_t cl = 0; cl < NhitsXYmod_[st][imod]; cl++){
2063 if(fDebug) cout<<
"st["<<st<<
"] "<<
" Y "<<YspCoord_[st][imod][cl]<<endl;
2066 Double_t FitX = CleanTr_.at(InIter).param[0] * CleanTr_.at(InIter).CoordZ[st] + CleanTr_.at(InIter).param[1];
2067 Double_t FitY = CleanTr_.at(InIter).param[2] * CleanTr_.at(InIter).CoordZ[st] + CleanTr_.at(InIter).param[3];
2069 if ( woSp[st] == 0){
2070 if(fDebug) cout<<
" woSp["<<st<<
"] "<<woSp[st]<<
"-> FitX "<<FitX<<
" FitY "<<FitY<<endl;
2071 XspCoord_[st][1][NhitsXYmod_[st][1]] = FitX;
2072 YspCoord_[st][1][NhitsXYmod_[st][1]] = FitY;
2073 SigmspX_[st][1][NhitsXYmod_[st][1]] = 0.0060;
2074 SigmspY_[st][1][NhitsXYmod_[st][1]] = 0.2;
2075 NhitsXYmod_[st][1]++;
2086 for(Int_t ist = 1; ist < fNstations; ist++){
2087 if(fDebug) cout<<
" pointsXsp["<<ist<<
"] ";
2088 for (Int_t imod = 0; imod < fNmodules; imod++) {
2089 for (Int_t cl = 0; cl < Nspatial[ist][imod]; cl++){
2090 if(fDebug) cout<<pointsXsp[ist][imod][cl]<<
" ";
2093 if(fDebug) cout<<endl;
2095 if(fDebug) cout<<endl;
2097 for(Int_t ist = 1; ist < fNstations; ist++){
2098 if(fDebug) cout<<
" pointsXpsp["<<ist<<
"] ";
2099 for (Int_t imod = 0; imod < fNmodules; imod++) {
2100 for (Int_t cl = 0; cl < Nspatial[ist][imod]; cl++){
2101 if(fDebug) cout<<pointsXpsp[ist][imod][cl]<<
" ";
2104 if(fDebug) cout<<endl;
2106 if(fDebug) cout<<endl;
2107 for(Int_t ist = 1; ist < fNstations; ist++) {
2108 if(fDebug) cout<<
" pointsX["<<ist<<
"] ";
2109 for (Int_t imod = 0; imod < fNmodules; imod++) {
2110 for (Int_t cl = 0; cl < NpointsX[ist][imod]; cl++){
2111 if(fDebug) cout<<pointsX[ist][imod][cl]<<
" ";
2114 if(fDebug) cout<<endl;
2116 for(Int_t ist = 1; ist < fNstations; ist++) {
2117 if(fDebug) cout<<
" pointsXp["<<ist<<
"] ";
2118 for (Int_t imod = 0; imod < fNmodules; imod++) {
2119 for (Int_t cl = 0; cl < NpointsXp[ist][imod]; cl++){
2120 if(fDebug) cout<<pointsXp[ist][imod][cl]<<
" ";
2123 if(fDebug) cout<<endl;
2127 for(Int_t ist = 1; ist < fNstations -1; ist++){
2128 for (Int_t imod = 0; imod < fNmodules; imod++) {
2129 if ( ist == 1 && imod > 3)
continue;
2130 if ( ist == 2 && imod > 1)
continue;
2133 if (NhitsXYmod_[ist][imod] > 0 ) woSp[ist] = 1;
2135 for (Int_t imod = 0; imod < fNmodules; imod++) {
2136 if (woSp[ist] == 1)
continue;
2137 if (fDebug) cout<<
" woSp["<<ist<<
"] "<<woSp[ist]<<endl;
2139 for (Int_t cl = 0; cl < NpointsX[ist][imod]; cl++){
2140 for(Int_t is = 1; is < fNstations -1 ; is++) {
2141 for (Int_t im = 0; im < fNmodules; im++) {
2142 if ( is == 1 && im > 3)
continue;
2143 if ( is == 2 && im > 1)
continue;
2145 for (Int_t clp = 0; clp < NpointsXp[is][im]; clp++){
2146 if (ist == 1 && is == 2){
2147 if(fDebug) cout<<
" X["<<ist<<
"] "<<pointsX[ist][imod][cl]<<
" Xp["<<is<<
"] "<< pointsXp[is][im][clp]<<
" XXp "<<pointsX[ist][imod][cl] - pointsXp[is][im][clp]<<endl;
2148 if (fDebug &&
fabs(pointsX[ist][imod][cl] - pointsXp[is][im][clp]) > 0.) hXXp12CheckLeftover->Fill(pointsX[ist][imod][cl] - pointsXp[is][im][clp]);
2149 if (fDebug &&
fabs(pointsX[ist][imod][cl] - pointsXp[is][im][clp]) > 0. &&
fabs(pointsX[ist][imod][cl] - pointsXp[is][im][clp]) < 0.2 ){
2150 hXXp12CheckLeftover03->Fill(pointsX[ist][imod][cl] - pointsXp[is][im][clp]);
2152 XspCoord_[ist][imod][NhitsXYmod_[ist][imod]] = pointsX[ist][imod][cl];
2153 YspCoord_[ist][imod][NhitsXYmod_[ist][imod]] =(pointsX[ist][imod][cl] - pointsXp[is][im][clp]) / tg2_5;
2154 SigmspX_[ist][imod][NhitsXYmod_[ist][imod]] = pointsXsig[ist][imod][cl];
2155 SigmspY_[ist][imod][NhitsXYmod_[ist][imod]] = 0.2;
2156 if(fDebug) cout<<
" Ysp "<<(pointsX[ist][imod][cl] - pointsXp[is][im][clp]) / tg2_5<<endl;
2157 NhitsXYmod_[ist][imod]++;
2166 if(fDebug) cout<<endl;
2169 for(Int_t ist = 1; ist < fNstations -1 ; ist++) {
2170 for (Int_t imod = 0; imod < fNmodules; imod++) {
2171 if ( ist == 1 && imod > 3)
continue;
2172 if ( ist == 2 && imod > 1)
continue;
2176 if (NhitsXYmod_[ist][imod] > 0 ) woSp[ist] = 1;
2178 for (Int_t imod = 0; imod < fNmodules; imod++) {
2179 if (woSp[ist] == 1)
continue;
2180 if (fDebug) cout<<
" woSp["<<ist<<
"] "<<woSp[ist]<<endl;
2182 for (Int_t cl = 0; cl < NpointsXp[ist][imod]; cl++){
2183 for(Int_t is = 1; is < fNstations -1; is++) {
2184 for (Int_t im = 0; im < fNmodules; im++) {
2185 if ( is == 1 && im > 3)
continue;
2186 if ( is == 2 && im > 1)
continue;
2188 for (Int_t clp = 0; clp < NpointsX[is][im]; clp++){
2189 if (ist == 1 && is == 2){
2190 if (fDebug) cout<<
" Xp["<<ist<<
"] "<<pointsXp[ist][imod][cl]<<
" X["<<is<<
"] "<< pointsX[is][im][clp]<<
" XXp "<<pointsXp[ist][imod][cl] - pointsX[is][im][clp]<<endl;
2191 if (fDebug &&
fabs(pointsXp[ist][imod][cl] - pointsX[is][im][clp]) > 0.) hXXp12CheckLeftover->Fill(pointsXp[ist][imod][cl] - pointsX[is][im][clp]);
2192 if (fDebug &&
fabs(pointsXp[ist][imod][cl] - pointsX[is][im][clp]) > 0. &&
fabs(pointsXp[ist][imod][cl] - pointsX[is][im][clp]) < 0.3 ){
2193 hXXp12CheckLeftover03->Fill(pointsXp[ist][imod][cl] - pointsX[is][im][clp]);
2195 XspCoord_[ist][imod][NhitsXYmod_[ist][imod]] = pointsX[is][im][clp];
2196 YspCoord_[ist][imod][NhitsXYmod_[ist][imod]] =(pointsX[is][im][clp] - pointsXp[ist][imod][cl]) / tg2_5;
2197 SigmspX_[ist][imod][NhitsXYmod_[ist][imod]] = pointsXsig[is][im][clp];
2198 SigmspY_[ist][imod][NhitsXYmod_[ist][imod]] = 0.2;
2199 if(fDebug) cout<<
" Ysp "<<(pointsX[is][im][clp] - pointsXp[ist][imod][cl]) / tg2_5<<endl;
2200 NhitsXYmod_[ist][imod]++;
2211 for(Int_t is = 1; is < fNstations -1; is++) {
2212 for (Int_t im = 0; im < fNmodules; im++) {
2213 if ( is == 1 && im > 3)
continue;
2214 if ( is == 2 && im > 1)
continue;
2216 for(Int_t icl = 0; icl < NhitsXYmod_[is][im]; icl++) {
2222 for (Int_t ist = 1; ist < fNstations ; ist++) {
2223 for (Int_t imod = 0; imod < fNmodules; imod++) {
2224 if (ist == 1 && imod > 3 )
continue;
2225 if (ist == 2 && imod > 1 )
continue;
2226 for (Int_t cl = 0; cl < NhitsXYmod_[ist][imod]; cl++) {
2227 if(fDebug && ist == 1) hY_XSisp1->Fill(XspCoord_[ist][imod][cl],YspCoord_[ist][imod][cl]);
2228 if(fDebug && ist == 2) hY_XSisp2->Fill(XspCoord_[ist][imod][cl],YspCoord_[ist][imod][cl]);
2239void BmnSiliconTrackFinder::PrintAllTracks(vector<tracksX> & CleanTr_){
2241 cout<<
"-----end of track finding ------ Ntr "<< CleanTr_.size()<<endl;
2242 Double_t ver_averx = 0., ver_avery = 0., sum_xv = 0., sum_yv = 0.;
2243 hNtracks->Fill(CleanTr_.size());
2244 for (
size_t InIter = 0; InIter < CleanTr_.size(); ++InIter){
2246 hNpoint->Fill(CleanTr_.at(InIter).Nhits );
2247 hAxglob->Fill(CleanTr_.at(InIter).param[0]);
2248 hBxglob->Fill(CleanTr_.at(InIter).param[1]);
2249 hAyglob->Fill(CleanTr_.at(InIter).param[2]);
2250 hByglob->Fill(CleanTr_.at(InIter).param[3]);
2254 hChiSquareNdf->Fill(CleanTr_.at(InIter).Chi2);
2255 if (InIter == 0) hAx_first_tr->Fill(CleanTr_.at(InIter).param[0]);
2256 else hAx_more_first_tr->Fill(CleanTr_.at(InIter).param[0]);
2257 hvertexXY->Fill(CleanTr_.at(InIter).Xv, CleanTr_.at(InIter).Yv);
2259 sum_xv += CleanTr_.at(InIter).Xv;
2260 sum_yv += CleanTr_.at(InIter).Yv;
2262 Double_t Z1 = Z0_SRC_target + 126.;
2263 Double_t Z2 = Z0_SRC_target + 176.;
2265 Double_t X1 = CleanTr_.at(InIter).param[0] * ( Z1 -Zcentr) + CleanTr_.at(InIter).param[1];
2266 Double_t Y1 = CleanTr_.at(InIter).param[2] * ( Z1 -Zcentr) + CleanTr_.at(InIter).param[3];
2268 Double_t X2 = CleanTr_.at(InIter).param[0] * ( Z2 -Zcentr) + CleanTr_.at(InIter).param[1];
2269 Double_t Y2 = CleanTr_.at(InIter).param[2] * ( Z2 -Zcentr) + CleanTr_.at(InIter).param[3];
2271 hYXSi1_run8->Fill(X1,Y1);
2272 hYXSi2_run8->Fill(X2,Y2);
2273 hXSi1_run8->Fill(X1);
2274 hYSi1_run8->Fill(Y1);
2275 hXSi2_run8->Fill(X2);
2276 hYSi2_run8->Fill(Y2);
2279 cout<<
" ---iTr "<<InIter<<
" Nhits "<<CleanTr_.at(InIter).Nhits<<
2280 " ax "<<CleanTr_.at(InIter).param[0]<<
" bx "<<CleanTr_.at(InIter).param[1]<<
" ay "<<CleanTr_.at(InIter).param[2] <<
" by "<<CleanTr_.at(InIter).param[3]<<
2281 " Chi2 "<<CleanTr_.at(InIter).Chi2<<endl;
2282 cout<<
" Xv "<<CleanTr_.at(InIter).Xv<<
" Yv "<<CleanTr_.at(InIter).Yv<<endl;
2283 cout<<
" --X--"<<endl;
2284 Double_t X1_0 = 999., X1_1= 999., X1_2= 999., X1_3= 999., X2_0= 999.,X2_1= 999.,X3_1= 999.,X3_2= 999.;
2286 for (Int_t st = 1; st < fNstations; ++st) {
2287 if ( CleanTr_.at(InIter).CoordX[st] < -900.)
continue;
2288 Double_t
Hit = CleanTr_.at(InIter).CoordX[st];
2289 Double_t Fit = CleanTr_.at(InIter).param[0] * CleanTr_.at(InIter).CoordZ[st] + CleanTr_.at(InIter).param[1];
2290 cout<<
" st "<<st<<
" mod "<<CleanTr_.at(InIter).ModNum[st]<<
" Xhit "<<
Hit<<
" Xfit "<< Fit<<
" dX "<<
Hit - Fit<<
" Z "<<CleanTr_.at(InIter).CoordZ[st]<<
" strip "<<CleanTr_.at(InIter).StripNumX[st]<<
" Clsize "<<CleanTr_.at(InIter).ClSizeX[st]<<
" Q "<<CleanTr_.at(InIter).SumQX[st]<<endl;
2293 hX_st_1 -> Fill(CleanTr_.at(InIter).CoordX[st]);
2294 hdXvsXst_1 -> Fill(Hit,Hit - Fit);
2295 hdX_st_1-> Fill(Hit - Fit);
2296 if (CleanTr_.at(InIter).ModNum[st] == 0 ) {
2297 X1_0 = CleanTr_.at(InIter).CoordX[1];
2298 hdXvsXst1_0 -> Fill(CleanTr_.at(InIter).CoordX[1],CleanTr_.at(InIter).CoordX[1] - Fit);
2300 if (CleanTr_.at(InIter).ModNum[st] == 1 ){
2301 X1_1 = CleanTr_.at(InIter).CoordX[1];
2302 hdXvsXst1_1 -> Fill(CleanTr_.at(InIter).CoordX[1],CleanTr_.at(InIter).CoordX[1] - Fit);
2304 if (CleanTr_.at(InIter).ModNum[st] == 2 ){
2305 X1_2 = CleanTr_.at(InIter).CoordX[1];
2306 hdXvsXst1_2 -> Fill(CleanTr_.at(InIter).CoordX[1],CleanTr_.at(InIter).CoordX[1] - Fit);
2308 if (CleanTr_.at(InIter).ModNum[st] == 3 ){
2309 X1_3 = CleanTr_.at(InIter).CoordX[1];
2310 hdXvsXst1_3 -> Fill(CleanTr_.at(InIter).CoordX[1],CleanTr_.at(InIter).CoordX[1] - Fit);
2314 hX_st_2 -> Fill(CleanTr_.at(InIter).CoordX[st]);
2315 hdXvsXst_2 -> Fill(Hit,Hit - Fit);
2316 hdX_st_2-> Fill(Hit - Fit);
2317 if (CleanTr_.at(InIter).ModNum[2] == 0 ){
2318 X2_0 = CleanTr_.at(InIter).CoordX[2];
2319 hdXvsXst2_0 -> Fill(CleanTr_.at(InIter).CoordX[2],CleanTr_.at(InIter).CoordX[2] - Fit);
2321 if (CleanTr_.at(InIter).ModNum[2] == 1 ){
2322 X2_1 = CleanTr_.at(InIter).CoordX[2];
2323 hdXvsXst2_1 -> Fill(CleanTr_.at(InIter).CoordX[2],CleanTr_.at(InIter).CoordX[2] - Fit);
2328 hX_st_3 -> Fill(CleanTr_.at(InIter).CoordX[st]);
2329 hdXvsXst_3 -> Fill(Hit,Hit - Fit);
2330 hdX_st_3-> Fill(Hit - Fit);
2331 if (CleanTr_.at(InIter).ModNum[3] == 1 ){
2332 X3_1 = CleanTr_.at(InIter).CoordX[3];
2333 hdXvsXst3_1 -> Fill(CleanTr_.at(InIter).CoordX[3],CleanTr_.at(InIter).CoordX[3] - Fit);
2334 hdXXp3_mod1 -> Fill(CleanTr_.at(InIter).CoordX[st] - CleanTr_.at(InIter).CoordXp[st]);
2336 if (CleanTr_.at(InIter).ModNum[3] == 2 ){
2337 X3_2 = CleanTr_.at(InIter).CoordX[3];
2338 hdXvsXst3_2 -> Fill(CleanTr_.at(InIter).CoordX[3],CleanTr_.at(InIter).CoordX[3] - Fit);
2339 hdXXp3_mod2 -> Fill(CleanTr_.at(InIter).CoordX[st] - CleanTr_.at(InIter).CoordXp[st]);
2343 if (X1_0 < 900. && X2_0 < 900.) hdXst1_0_st2_0->Fill(X1_0 - X2_0);
2344 if (X1_0 < 900. && X2_1 < 900.) hdXst1_0_st2_1->Fill(X1_0 - X2_1);
2345 if (X1_1 < 900. && X2_0 < 900.) hdXst1_1_st2_0->Fill(X1_1 - X2_0);
2346 if (X1_1 < 900. && X2_1 < 900.) hdXst1_1_st2_1->Fill(X1_1 - X2_1);
2347 if (X1_2 < 900. && X2_0 < 900.) hdXst1_2_st2_0->Fill(X1_2 - X2_0);
2348 if (X1_2 < 900. && X2_1 < 900.) hdXst1_2_st2_1->Fill(X1_2 - X2_1);
2349 if (X1_3 < 900. && X2_0 < 900.) hdXst1_3_st2_0->Fill(X1_3 - X2_0);
2350 if (X1_3 < 900. && X2_1 < 900.) hdXst1_3_st2_1->Fill(X1_3 - X2_1);
2351 if (X2_0 < 900. && X3_1 < 900.) hdXst2_0_st3_1->Fill(X2_0 - X3_1);
2352 if (X2_0 < 900. && X3_2 < 900.) hdXst2_0_st3_2->Fill(X2_0 - X3_2);
2353 if (X2_1 < 900. && X3_1 < 900.) hdXst2_1_st3_1->Fill(X2_1 - X3_1);
2354 if (X2_1 < 900. && X3_2 < 900.) hdXst2_1_st3_2->Fill(X2_1 - X3_2);
2356 if(CleanTr_.at(InIter).Nhits == 6 ){
2357 Double_t Ax = (CleanTr_.at(InIter).CoordX[1] -CleanTr_.at(InIter).CoordX[3])/(CleanTr_.at(InIter).CoordZ[1] -CleanTr_.at(InIter).CoordZ[3]);
2358 if (X2_0 < 900.) hX13_X2_m0->Fill( X2_0 - (Ax*(CleanTr_.at(InIter).CoordZ[2] - CleanTr_.at(InIter).CoordZ[1]) + CleanTr_.at(InIter).CoordX[1]) );
2359 if (X2_1 < 900.) hX13_X2_m1->Fill( X2_1 - (Ax*(CleanTr_.at(InIter).CoordZ[2] - CleanTr_.at(InIter).CoordZ[1]) + CleanTr_.at(InIter).CoordX[1]) );
2361 hdX_st1_st2->Fill( CleanTr_.at(InIter).CoordX[1] -CleanTr_.at(InIter).CoordX[2] );
2362 hdX_st2_st3->Fill( CleanTr_.at(InIter).CoordX[2] -CleanTr_.at(InIter).CoordX[3] );
2363 hdX_st1_st3->Fill( CleanTr_.at(InIter).CoordX[1] -CleanTr_.at(InIter).CoordX[3] );
2366 cout<<
" --Xp--"<<endl;
2367 Double_t Xp2_0= 999.,Xp2_1= 999.;
2368 for (Int_t st = 1; st < fNstations; ++st) {
2369 if ( CleanTr_.at(InIter).CoordXp[st] < -900.)
continue;
2370 Double_t
Hit = CleanTr_.at(InIter).CoordXp[st];
2371 Double_t Fit = ( CleanTr_.at(InIter).param[0] * CleanTr_.at(InIter).CoordZ[st] + CleanTr_.at(InIter).param[1] ) +
2372 ( CleanTr_.at(InIter).param[2] * CleanTr_.at(InIter).CoordZ[st] + CleanTr_.at(InIter).param[3] ) * Angle(st,CleanTr_.at(InIter).ModNum[st]) * tg2_5 ;
2374 hdXp_st_1-> Fill(Hit - Fit);
2375 cout<<
" st "<<st<<
" Xphit "<<
Hit<<
" Xpfit "<< Fit<<
" dXp "<<
Hit - Fit<<
" Z "<<CleanTr_.at(InIter).CoordZ[st]<<
" strip "<<CleanTr_.at(InIter).StripNumXp[st]<<
" Clsize "<<CleanTr_.at(InIter).ClSizeXp[st]<<
" Q "<<CleanTr_.at(InIter).SumQXp[st]<<endl;
2377 hdXp_st_2-> Fill(Hit - Fit);
2378 if (CleanTr_.at(InIter).ModNum[2] == 0 ){
2379 Xp2_0 = CleanTr_.at(InIter).CoordXp[2];
2381 if (CleanTr_.at(InIter).ModNum[2] == 1 ){
2382 Xp2_1 = CleanTr_.at(InIter).CoordXp[2];
2386 hdXp_st_3-> Fill(Hit - Fit);
2387 if ( CleanTr_.at(InIter).ModNum[3] == 1 ){
2388 hdXp3_mod1 -> Fill(Hit - Fit);
2390 if ( CleanTr_.at(InIter).ModNum[3] == 2 ){
2391 hdXp3_mod2 -> Fill(Hit - Fit);
2395 if(CleanTr_.at(InIter).Nhits == 6 ){
2396 Double_t Axp = (CleanTr_.at(InIter).CoordXp[1] -CleanTr_.at(InIter).CoordXp[3])/(CleanTr_.at(InIter).CoordZ[1] -CleanTr_.at(InIter).CoordZ[3]);
2397 if (Xp2_0 < 900.) hXp13_Xp2_m0->Fill( Xp2_0 - (Axp*(CleanTr_.at(InIter).CoordZ[2] - CleanTr_.at(InIter).CoordZ[1]) + CleanTr_.at(InIter).CoordXp[1]) );
2398 if (Xp2_1 < 900.) hXp13_Xp2_m1->Fill( Xp2_1 - (Axp*(CleanTr_.at(InIter).CoordZ[2] - CleanTr_.at(InIter).CoordZ[1]) + CleanTr_.at(InIter).CoordXp[1]) );
2401 if (fDebug) cout<<
" --Y--"<<endl;
2402 Double_t Y1_0 = 999., Y1_1= 999., Y1_2= 999., Y1_3= 999., Y2_0= 999.,Y2_1= 999.,Y3_1= 999.,Y3_2= 999.;
2404 for (Int_t st = 1; st < fNstations; ++st) {
2405 if ( CleanTr_.at(InIter).CoordY[st] < -900. )
continue;
2406 Double_t
Hit = CleanTr_.at(InIter).CoordY[st];
2407 Double_t Fit = CleanTr_.at(InIter).param[2] * CleanTr_.at(InIter).CoordZ[st] + CleanTr_.at(InIter).param[3];
2408 cout<<
" st "<<st<<
" Yhit "<<
Hit<<
" Yfit "<< Fit<<
" dY "<<
Hit - Fit<<
" Z "<<CleanTr_.at(InIter).CoordZ[st]<<endl;
2410 if (st == 1 ) {hY_st_1 -> Fill(CleanTr_.at(InIter).CoordY[st]);hdYvsYst_1 -> Fill(Hit,Hit - Fit); hdY_st_1-> Fill(Hit - Fit);}
2411 if (st == 2 ) {hY_st_2 -> Fill(CleanTr_.at(InIter).CoordY[st]);hdYvsYst_2 -> Fill(Hit,Hit - Fit); hdY_st_2-> Fill(Hit - Fit);}
2412 if (st == 3 ) {hY_st_3 -> Fill(CleanTr_.at(InIter).CoordY[st]);hdYvsYst_3 -> Fill(CleanTr_.at(InIter).CoordY[3],CleanTr_.at(InIter).CoordY[3] - Fit);
2413 hdY_st_3-> Fill(CleanTr_.at(InIter).CoordY[3] - Fit);
2414 if (
fabs(Hit - Fit) > 0.1) cout<<
" > 0.1 "<<endl;
2419 if ( CleanTr_.at(InIter).ModNum[1] == 0 ){
2420 hdYvsYst1_mod0 -> Fill(Hit,Hit - Fit);
2421 Y1_0 = CleanTr_.at(InIter).CoordY[st];
2423 if ( CleanTr_.at(InIter).ModNum[1] == 1 ){
2424 hdYvsYst1_mod1 -> Fill(Hit,Hit - Fit);
2425 Y1_1 = CleanTr_.at(InIter).CoordY[st];
2427 if ( CleanTr_.at(InIter).ModNum[1] == 2 ){
2428 hdYvsYst1_mod2 -> Fill(Hit,Hit - Fit);
2429 Y1_2 = CleanTr_.at(InIter).CoordY[st];
2431 if ( CleanTr_.at(InIter).ModNum[1] == 3 ){
2432 hdYvsYst1_mod3 -> Fill(Hit,Hit - Fit);
2433 Y1_3 = CleanTr_.at(InIter).CoordY[st];
2437 if (CleanTr_.at(InIter).ModNum[2] == 0 ){
2438 hdYvsYst2_mod0 -> Fill(Hit,Hit - Fit);
2439 Y2_0 = CleanTr_.at(InIter).CoordY[st];
2441 if (CleanTr_.at(InIter).ModNum[2] == 1 ){
2442 hdYvsYst2_mod1 -> Fill(Hit,Hit - Fit);
2443 Y2_1 = CleanTr_.at(InIter).CoordY[st];
2448 if ( CleanTr_.at(InIter).ModNum[3] == 1 ){
2449 hdYvsYst3_mod1 -> Fill(Hit,Hit - Fit);
2450 Y3_1 = CleanTr_.at(InIter).CoordY[st];
2452 if ( CleanTr_.at(InIter).ModNum[3] == 2 ){
2453 hdYvsYst3_mod2 -> Fill(Hit,Hit - Fit);
2454 Y3_2 = CleanTr_.at(InIter).CoordY[st];
2459 double Ay_pr= -999.;
double Y_pr = -999.;
2460 if (CleanTr_.at(InIter).CoordY[1] > -900. && CleanTr_.at(InIter).CoordY[3] > -900.){
2461 Ay_pr = (CleanTr_.at(InIter).CoordY[1] - CleanTr_.at(InIter).Yv)/(CleanTr_.at(InIter).CoordZ[1] - Z0_SRC_target);
2462 Y_pr = Ay_pr* (CleanTr_.at(InIter).CoordZ[3] - CleanTr_.at(InIter).CoordZ[1]) + CleanTr_.at(InIter).CoordY[1];
2464 if (CleanTr_.at(InIter).CoordY[2] > -900. && CleanTr_.at(InIter).CoordY[3] > -900.){
2465 Ay_pr = (CleanTr_.at(InIter).CoordY[2] - CleanTr_.at(InIter).Yv)/(CleanTr_.at(InIter).CoordZ[2] - Z0_SRC_target);
2466 Y_pr = Ay_pr* (CleanTr_.at(InIter).CoordZ[3] - CleanTr_.at(InIter).CoordZ[2]) + CleanTr_.at(InIter).CoordY[2];
2468 cout<<
" Y3 predicted "<<Y_pr<<
" Y "<<CleanTr_.at(InIter).CoordY[3]<<endl;
2470 hdY_st1_st2->Fill( CleanTr_.at(InIter).CoordY[1] -CleanTr_.at(InIter).CoordY[2] );
2471 hdY_st2_st3->Fill( CleanTr_.at(InIter).CoordY[2] -CleanTr_.at(InIter).CoordY[3] );
2472 hdY_st1_st3->Fill( CleanTr_.at(InIter).CoordY[1] -CleanTr_.at(InIter).CoordY[3] );
2475 if (Y1_0 < 900. && Y2_0 < 900.) hdYst1_0_st2_0->Fill(Y1_0 - Y2_0);
2476 if (Y1_0 < 900. && Y2_1 < 900.) hdYst1_0_st2_1->Fill(Y1_0 - Y2_1);
2477 if (Y1_1 < 900. && Y2_0 < 900.) hdYst1_1_st2_0->Fill(Y1_1 - Y2_0);
2478 if (Y1_1 < 900. && Y2_1 < 900.) hdYst1_1_st2_1->Fill(Y1_1 - Y2_1);
2479 if (Y1_2 < 900. && Y2_0 < 900.) hdYst1_2_st2_0->Fill(Y1_2 - Y2_0);
2480 if (Y1_2 < 900. && Y2_1 < 900.) hdYst1_2_st2_1->Fill(Y1_2 - Y2_1);
2481 if (Y1_3 < 900. && Y2_0 < 900.) hdYst1_3_st2_0->Fill(Y1_3 - Y2_0);
2482 if (Y1_3 < 900. && Y2_1 < 900.) hdYst1_3_st2_1->Fill(Y1_3 - Y2_1);
2483 if (Y2_0 < 900. && Y3_1 < 900.) hdYst2_0_st3_1->Fill(Y2_0 - Y3_1);
2484 if (Y2_0 < 900. && Y3_2 < 900.) hdYst2_0_st3_2->Fill(Y2_0 - Y3_2);
2485 if (Y2_1 < 900. && Y3_1 < 900.) hdYst2_1_st3_1->Fill(Y2_1 - Y3_1);
2486 if (Y2_1 < 900. && Y3_2 < 900.) hdYst2_1_st3_2->Fill(Y2_1 - Y3_2);
2488 if(CleanTr_.at(InIter).Nhits == 6 ){
2489 Double_t Ax = (CleanTr_.at(InIter).CoordY[1] -CleanTr_.at(InIter).CoordY[3])/(CleanTr_.at(InIter).CoordZ[1] -CleanTr_.at(InIter).CoordZ[3]);
2490 Double_t Ax23 = (CleanTr_.at(InIter).CoordY[2] -CleanTr_.at(InIter).CoordY[3])/(CleanTr_.at(InIter).CoordZ[2] -CleanTr_.at(InIter).CoordZ[3]);
2492 if (Y1_0 < 900.) hY1m0_Y23->Fill( Y1_0 - (Ax23*(CleanTr_.at(InIter).CoordZ[1] - CleanTr_.at(InIter).CoordZ[2]) + CleanTr_.at(InIter).CoordY[2]) );
2493 if (Y1_1 < 900.) hY1m1_Y23->Fill( Y1_1 - (Ax23*(CleanTr_.at(InIter).CoordZ[1] - CleanTr_.at(InIter).CoordZ[2]) + CleanTr_.at(InIter).CoordY[2]) );
2494 if (Y1_2 < 900.) hY1m2_Y23->Fill( Y1_2 - (Ax23*(CleanTr_.at(InIter).CoordZ[1] - CleanTr_.at(InIter).CoordZ[2]) + CleanTr_.at(InIter).CoordY[2]) );
2495 if (Y1_3 < 900.) hY1m3_Y23->Fill( Y1_3 - (Ax23*(CleanTr_.at(InIter).CoordZ[1] - CleanTr_.at(InIter).CoordZ[2]) + CleanTr_.at(InIter).CoordY[2]) );
2497 if (Y2_0 < 900.) hY13_Y2_m0->Fill( Y2_0 - (Ax*(CleanTr_.at(InIter).CoordZ[2] - CleanTr_.at(InIter).CoordZ[1]) + CleanTr_.at(InIter).CoordY[1]) );
2498 if (Y2_1 < 900.) hY13_Y2_m1->Fill( Y2_1 - (Ax*(CleanTr_.at(InIter).CoordZ[2] - CleanTr_.at(InIter).CoordZ[1]) + CleanTr_.at(InIter).CoordY[1]) );
2502 if (CleanTr_.at(InIter).CoordX[1] > -900. || CleanTr_.at(InIter).CoordXp[1] > -900.)
2503 hprofile_beam_z1->Fill(CleanTr_.at(InIter).param[0] * (CleanTr_.at(InIter).CoordZ[1] - z_cnt) + CleanTr_.at(InIter).param[1],CleanTr_.at(InIter).param[2] * (CleanTr_.at(InIter).CoordZ[1]- z_cnt) + CleanTr_.at(InIter).param[3] );
2504 if (CleanTr_.at(InIter).CoordX[2] > -900. || CleanTr_.at(InIter).CoordXp[2] > -900.)
2505 hprofile_beam_z2->Fill(CleanTr_.at(InIter).param[0] * (CleanTr_.at(InIter).CoordZ[2] - z_cnt) + CleanTr_.at(InIter).param[1],CleanTr_.at(InIter).param[2] * (CleanTr_.at(InIter).CoordZ[2]- z_cnt) + CleanTr_.at(InIter).param[3] );
2506 hprofile_beam_z3->Fill(CleanTr_.at(InIter).param[0] * (CleanTr_.at(InIter).CoordZ[3] - z_cnt) + CleanTr_.at(InIter).param[1],CleanTr_.at(InIter).param[2] * (CleanTr_.at(InIter).CoordZ[3]- z_cnt) + CleanTr_.at(InIter).param[3] );
2509 for (Int_t st = 1; st < fNstations; ++st) {
2511 if (CleanTr_.at(InIter).Nhits == N_min_points) {
2512 if ( CleanTr_.at(InIter).CoordX[st] < -900. ) D_eff->Fill(2*st-1);
2513 if ( CleanTr_.at(InIter).CoordXp[st] < -900. ) D_eff->Fill(2*st);
2515 if (CleanTr_.at(InIter).Nhits > N_min_points) {
2516 D_eff->Fill(2*st-1);
2518 if ( CleanTr_.at(InIter).CoordX[st] > -900. ) N_eff->Fill(2*st-1);
2519 if ( CleanTr_.at(InIter).CoordXp[st] > -900. ) N_eff->Fill(2*st);
2525 if (CleanTr_.size() > 0) {
2526 ver_averx = sum_xv/CleanTr_.size();
2527 ver_avery = sum_yv/CleanTr_.size();
2528 hvertex_aver_XY->Fill(ver_averx, ver_avery);
2530 cout<<
"-----end find tracks ------ Ntr "<< CleanTr_.size()<<endl;
2538void BmnSiliconTrackFinder::Case2( Int_t & Fi_St,Int_t **Nleftoversp_,Double_t ***leftoverXsp_, Double_t ***leftoverXpsp_, Double_t ***leftoverYsp_,
2539 Double_t ***leftoverXsigsp_, Double_t ***leftoverXpsigsp_, Double_t ***leftoverYsigsp_,
2540 Int_t **NleftoverX_, Int_t **NleftoverXp_,Double_t ***leftoverX_, Double_t ***leftoverXp_,Double_t ***leftoverXsig_, Double_t ***leftoverXpsig_,
2541 vector<tracksX> & vec_tracks_,
2542 Double_t ***XClust_, Double_t ***XClust_width_, Double_t ***sumQx_, Double_t ***XpClust_, Double_t ***XpClust_width_, Double_t ***sumQxp_,
2543 Double_t ***XspClust_, Double_t ***XspClust_width_, Double_t ***sumQxsp_, Double_t ***XpspClust_, Double_t ***XpspClust_width_, Double_t ***sumQxpsp_){
2545 if(fDebug) cout<<
" case2:"<<Fi_St<<
"time "<<endl;
2549 Double_t xv_candidate = 999., yv_candidate = 999., ax_candidate = 999., ay_candidate = 999.;
2550 Int_t iClustXMod_candidate = -1, iModX_candidate = -1;
2551 Double_t dxV_thr = .2, dyV_thr = .4;
2553 Int_t TotalNumberOfHits = 0;
2554 Double_t ChiSquareNdf = 0.;
2556 Bool_t x5_exist = 0;
2557 Bool_t xp5_exist = 0;
2558 Int_t Nmod_Fi = fNmodules1;
2559 Int_t Nmod_Sec = fNmodules2;
2563 Nmod_Fi = fNmodules2;
2564 Nmod_Sec = fNmodules1;
2568 for (Int_t imod1 = 0; imod1 < Nmod_Fi; imod1++) {
2569 for (Int_t cl1 = 0; cl1 < Nleftoversp_[Fi_St][imod1]; cl1++){
2571 Double_t ax_target_region = -999., bx_target_region= -999., ay_target_region = -999., by_target_region= -999.;
2572 Bool_t was_candidate = 0; ax_candidate = 999., ay_candidate = 999., xv_candidate = 999., yv_candidate = 999.;
2578 for (Int_t imod3 = 0; imod3 < fNmodules3; imod3++) {
2579 for (Int_t cl3 = 0; cl3 < Nleftoversp_[3][imod3]; cl3++){
2581 for ( Int_t
i = 1;
i < fNstations; ++
i){
2586 ax_target_region = (leftoverXsp_[3][imod3][cl3] - leftoverXsp_[Fi_St][imod1][cl1] ) / (Zstation[3][imod3] - Zstation[Fi_St][imod1] );
2587 bx_target_region = -ax_target_region * ( Zstation[3][imod3] - Z0_SRC_target ) + leftoverXsp_[3][imod3][cl3];
2589 ay_target_region = (leftoverYsp_[3][imod3][cl3] - leftoverYsp_[Fi_St][imod1][cl1] ) / (Zstation[3][imod3] - Zstation[Fi_St][imod1] );
2590 by_target_region = -ay_target_region * ( Zstation[3][imod3] - Z0_SRC_target ) + leftoverYsp_[3][imod3][cl3];
2592 if (fDebug) cout<<
" bx_target_region "<<bx_target_region<<
" by_target_region "<<by_target_region<<endl;
2594 if (
fabs(bx_target_region) > half_target_regionX)
continue;
2595 if (
fabs(by_target_region) > half_target_regionY)
continue;
2597 if ( was_candidate == 0 ) {
2599 xv_candidate = bx_target_region;
2600 yv_candidate = by_target_region;
2601 ax_candidate = ax_target_region;
2602 ay_candidate = ay_target_region;
2603 iModX_candidate = imod3;
2604 iClustXMod_candidate = cl3;
2608 if (
fabs(bx_target_region) > (
fabs(xv_candidate) + dxV_thr) )
continue;
2609 if (
fabs(by_target_region) > (
fabs(yv_candidate) + dyV_thr) )
continue;
2610 if ( (
fabs(bx_target_region) + 0.5*
fabs(by_target_region) ) > (
fabs(xv_candidate) + 0.25*
fabs(yv_candidate) ) )
continue;
2612 xv_candidate = bx_target_region;
2613 yv_candidate = by_target_region;
2614 ax_candidate = ax_target_region;
2615 ay_candidate = ay_target_region;
2617 iModX_candidate = imod3;
2618 iClustXMod_candidate = cl3;
2628 if (fDebug) cout<<
" aft.st3 cycle: xv_candidate "<<xv_candidate<<
" by_cand "<<yv_candidate<<
" was_candidate "<<was_candidate<<endl;
2629 if (xv_candidate > 900.)
continue;
2631 for (Int_t ist = 1; ist < fNstations; ++ist) {
2632 for (Int_t imod = 0; imod < fNmodules; imod++) {
2633 Xforglfit[ist][imod] = -999.;
2634 Yforglfit[ist][imod] = -999.;
2635 Xpforglfit[ist][imod] = -999.;
2636 SigmYforglfit[ist][imod] = -999.;
2637 SigmXforglfit[ist][imod] = -999.;
2638 SigmXpforglfit[ist][imod] = -999.;
2642 if ( iModX_candidate > -1){
2643 iModX[Fi_St] = imod1;
2644 iModX[3] = iModX_candidate;
2645 iClustXMod[Fi_St] = cl1;
2646 iClustXMod[3] = iClustXMod_candidate;
2649 for (Int_t ist = 1; ist < fNstations; ++ist) {
2650 if ( ist == Sec_St )
continue;
2651 Xforglfit[ist][iModX[ist]] = leftoverXsp_[ist][iModX[ist]][iClustXMod[ist]];
2652 SigmXforglfit[ist][iModX[ist]] = leftoverXsigsp_[ist][iModX[ist]][iClustXMod[ist]];
2653 Xpforglfit[ist][iModX[ist]] = leftoverXpsp_[ist][iModX[ist]][iClustXMod[ist]];
2654 SigmXpforglfit[ist][iModX[ist]] = leftoverXpsigsp_[ist][iModX[ist]][iClustXMod[ist]];
2655 Yforglfit[ist][iModX[ist]] = leftoverYsp_[ist][iModX[ist]][iClustXMod[ist]];
2656 SigmYforglfit[ist][iModX[ist]] = leftoverYsigsp_[ist][iModX[ist]][iClustXMod[ist]];
2660 if ( fDebug) cout<<
" try to search 5th p. Sec_St "<<Sec_St<<
" Nmod_Sec "<<Nmod_Sec<<endl;
2661 for (Int_t imod2 = 0; imod2 < Nmod_Sec; imod2++) {
2664 if ( x5_exist ) { Npoint = 5;
break;}
2666 for (Int_t clx2 = 0; clx2 < NleftoverX_[Sec_St][imod2]; clx2++){
2670 iClustXMod[Sec_St] = -1;
2671 Xforglfit[Sec_St][iModX[Sec_St]] = -999.;
2672 SigmXforglfit[Sec_St][iModX[Sec_St]] = -999.;
2673 Xpforglfit[Sec_St][iModX[Sec_St]] = -999.;
2674 SigmXpforglfit[Sec_St][iModX[Sec_St]] = -999.;
2676 if ( fDebug) cout<<
"search 5p in X. tmp x1 "<<leftoverXsp_[Fi_St][imod1][cl1] <<
" x2 "<<leftoverX_[Sec_St][imod2][clx2]<<
" hrd "<<half_roadX1_X2<<endl;
2677 if (
fabs(leftoverXsp_[Fi_St][imod1][cl1] - leftoverX_[Sec_St][imod2][clx2]) > half_roadX1_X2 )
continue;
2678 if ( fDebug) cout<<
"found 5th p in X"<<Sec_St<<
" "<<leftoverX_[Sec_St][imod2][clx2]<<endl;
2682 if ( clx2 < (NleftoverX_[Sec_St][imod2] -1) ) {
2683 if (
fabs(leftoverXsp_[Fi_St][imod1][cl1]-leftoverX_[Sec_St][imod2][clx2]) >
fabs(leftoverXsp_[Fi_St][imod1][cl1]-leftoverX_[Sec_St][imod2][clx2+1]) )
2687 iModX[Sec_St] = imod2;
2688 iClustXMod[Sec_St] = clx2;
2689 Xforglfit[Sec_St][iModX[Sec_St]] = leftoverX_[Sec_St][iModX[Sec_St]][clx2];
2690 SigmXforglfit[Sec_St][iModX[Sec_St]] = leftoverXsig_[Sec_St][iModX[Sec_St]][clx2];
2691 Xpforglfit[Sec_St][iModX[Sec_St]] = -999.;
2692 SigmXpforglfit[Sec_St][iModX[Sec_St]] = -999.;
2695 TotalNumberOfHits = 0;
2699 for (Int_t coeff = 0; coeff < 4; ++coeff) {
2705 calculationOfChi2(Xforglfit, Xpforglfit, SigmXforglfit, SigmXpforglfit, iClustXMod, iModX, ChiSquare, line);
2707 TotalNumberOfHits = Npoint;
2709 ChiSquareNdf = ChiSquare/(TotalNumberOfHits - 4);
2710 else ChiSquareNdf = ChiSquare;
2712 if ( fDebug) cout<<
" before cut ChiSquareNDF "<<ChiSquareNdf<<
" TotalNumberOfHits "<<TotalNumberOfHits <<endl;
2713 if (fDebug) cout<<endl;
2715 if ( ChiSquareNdf > Chi2_Globcut ) {
2718 iClustXMod[Sec_St] = -1;
2719 Xforglfit[Sec_St][iModX[Sec_St]] = -999.;
2720 SigmXforglfit[Sec_St][iModX[Sec_St]] = -999.;
2721 Xpforglfit[Sec_St][imod2] = -999.;
2722 SigmXpforglfit[Sec_St][imod2] = -999.;
2728 if (fDebug) cout<<
" case2: x5_exist "<<endl;
2729 if (fDebug) cout<<
" X "<<Xforglfit[Sec_St][iModX[Sec_St]]<<
" strip "<<XClust_[Sec_St][iModX[Sec_St]][iClustXMod[Sec_St]]<<
" Q "<<sumQx_[Sec_St][iModX[Sec_St]][iClustXMod[Sec_St]]<<endl;
2735 if ( x5_exist )
continue;
2738 for (Int_t clxp2 = 0; clxp2 < NleftoverXp_[Sec_St][imod2]; clxp2++){
2739 if ( xp5_exist ) { Npoint = 5;
break;}
2742 iClustXMod[Sec_St] = -1;
2743 Xforglfit[Sec_St][iModX[Sec_St]] = -999.;
2744 SigmXforglfit[Sec_St][iModX[Sec_St]] = -999.;
2745 Xpforglfit[Sec_St][imod2] = -999.;
2746 SigmXpforglfit[Sec_St][imod2] = -999.;
2748 Int_t st1 = Fi_St; Int_t st2 = Sec_St;
2749 Int_t sign1 = Angle(st1,imod1);
2750 Int_t sign2 = Angle(st2,imod2);
2751 if(fDebug) cout<<
" sign1 "<<sign1<<
" sign2 "<<sign2<<endl;
2753 if ( (sign1 == sign2) && (
fabs(leftoverXpsp_[Fi_St][imod1][cl1] - leftoverXp_[Sec_St][imod2][clxp2] ) > half_roadXp1_Xp2) )
continue;
2754 if ( (sign1 != sign2) && (
fabs(leftoverXpsp_[Fi_St][imod1][cl1] - leftoverXp_[Sec_St][imod2][clxp2] ) > half_roadXp1_Xp2_diff_sign) )
continue;
2755 if ( fDebug) cout<<
" --Xp1 "<<leftoverXpsp_[Fi_St][imod1][cl1]<<
" Xp2 "<<leftoverXp_[Sec_St][imod2][clxp2]<<
" dXp "<<leftoverXpsp_[Fi_St][imod1][cl1] - leftoverXp_[Sec_St][imod2][clxp2]<<endl;
2759 iModX[Sec_St] = imod2;
2760 iClustXMod[Sec_St] = clxp2;
2761 Xforglfit[Sec_St][iModX[Sec_St]] = -999.;
2762 SigmXforglfit[Sec_St][iModX[Sec_St]] = -999.;
2763 Xpforglfit[Sec_St][imod2] = leftoverXp_[Sec_St][iModX[Sec_St]][clxp2];
2764 SigmXpforglfit[Sec_St][imod2] = leftoverXpsig_[Sec_St][iModX[Sec_St]][clxp2];
2767 TotalNumberOfHits = 0;
2771 for (Int_t coeff = 0; coeff < 4; ++coeff) {
2775 calculationOfChi2(Xforglfit, Xpforglfit, SigmXforglfit, SigmXpforglfit, iClustXMod, iModX, ChiSquare, line);
2777 TotalNumberOfHits = Npoint;
2779 ChiSquareNdf = ChiSquare/(TotalNumberOfHits - 4);
2780 else ChiSquareNdf = ChiSquare;
2782 if ( fDebug) cout<<
" before cut ChiSquareNDF "<<ChiSquareNdf<<
" TotalNumberOfHits "<<TotalNumberOfHits <<endl;
2783 if (fDebug) cout<<endl;
2785 if ( ChiSquareNdf > Chi2_Globcut ){
2788 iClustXMod[Sec_St] = -1;
2789 Xforglfit[Sec_St][iModX[Sec_St]] = -999.;
2790 SigmXforglfit[Sec_St][iModX[Sec_St]] = -999.;
2791 Xpforglfit[Sec_St][imod2] = -999.;
2792 SigmXpforglfit[Sec_St][imod2] = -999.;
2798 if (fDebug) cout<<
" case2: xp5_exist "<<xp5_exist<<endl;
2799 if (fDebug) cout<<
" Xp "<<Xpforglfit[Sec_St][iModX[Sec_St]]<<
" strip "<<XpClust_[Sec_St][iModX[Sec_St]][iClustXMod[Sec_St]]<<
" Q "<<sumQxp_[Sec_St][iModX[Sec_St]][iClustXMod[Sec_St]]<<endl;
2803 if ( xp5_exist )
break;
2807 if (fDebug) cout<<
" case2: end search ------------"<<endl;
2808 if (fDebug) cout<<
" Npoint "<<Npoint<<
" ChiSquareNdf "<<ChiSquareNdf<<
" x5_exist "<<x5_exist<<
" xp5_exist "<<xp5_exist<<endl;
2810 TotalNumberOfHits = Npoint;
2811 Int_t Mod_cand[fNstations];
2812 Double_t X_cand[fNstations], Xp_cand[fNstations], Y_cand[fNstations],
2813 SigX_cand[fNstations], SigXp_cand[fNstations], SigY_cand[fNstations],
2814 Chi2_cand = 999., line_cand[4];
2816 for (Int_t coeff = 0; coeff < 4; ++coeff) {
2817 line_cand[coeff] = 0.;
2820 for (Int_t st = 1; st < fNstations; ++st) {
2823 Xp_cand[st] = -999.;
2825 SigX_cand[st] = -999.;
2826 SigXp_cand[st] = -999.;
2827 SigY_cand[st] = -999.;
2830 if (TotalNumberOfHits > 4 && ChiSquareNdf < Chi2_cand && ChiSquareNdf < Chi2_Globcut ){
2832 Chi2_cand = ChiSquareNdf;
2833 line_cand[0] = line[0];
2834 line_cand[1] = line[1];
2835 line_cand[2] = line[2];
2836 line_cand[3] = line[3];
2838 if (fDebug) cout<<
" case2 5p cand Chi2Ndf "<<Chi2_cand<<endl;
2840 for (Int_t st = 1; st < fNstations; ++st) {
2841 if (iModX[st] < 0)
continue;
2843 Mod_cand[st] = iModX[st];
2844 X_cand[st] = Xforglfit[st][iModX[st]];
2845 Y_cand[st] = Yforglfit[st][iModX[st]];
2846 Xp_cand[st] = Xpforglfit[st][iModX[st]];
2847 SigX_cand[st] = SigmXforglfit[st][iModX[st]];
2848 SigY_cand[st] = SigmYforglfit[st][iModX[st]];
2849 SigXp_cand[st] = SigmXpforglfit[st][iModX[st]];
2851 if (fDebug) cout<<
" X_cand["<<st<<
"] "<< X_cand[st] <<
" Xp_cand["<<st<<
"] "<<Xp_cand[st] <<
" iModX "<<iModX[st]<<endl;
2856 if ( TotalNumberOfHits == 4) {
2857 if (fDebug) cout<<
"TotalNumberOfHits == 4"<<endl;
2861 line_cand[0] = ax_candidate;
2862 line_cand[1] = xv_candidate + ax_candidate * ( Zcentr - Z0_SRC_target );
2863 line_cand[2] = ay_candidate;
2864 line_cand[3] = yv_candidate + ay_candidate * ( Zcentr - Z0_SRC_target );
2866 if (fDebug) cout<<
" case2 4p cand Xv "<<xv_candidate<<
" Yv "<<yv_candidate<<endl;
2867 if (fDebug) cout<<
" ax "<<line_cand[0]<<
" bx "<<line_cand[1]<<
" ay "<<line_cand[2]<<
" by "<<line_cand[3]<<endl;
2869 for (Int_t st = 1; st < fNstations; ++st) {
2870 if (iModX[st] < 0)
continue;
2872 Mod_cand[st] = iModX[st];
2873 X_cand[st] = Xforglfit[st][iModX[st]];
2874 Y_cand[st] = Yforglfit[st][iModX[st]];
2875 Xp_cand[st] = Xpforglfit[st][iModX[st]];
2876 SigX_cand[st] = SigmXforglfit[st][iModX[st]];
2877 SigY_cand[st] = SigmYforglfit[st][iModX[st]];
2878 SigXp_cand[st] = SigmXpforglfit[st][iModX[st]];
2880 if (fDebug) cout<<
" X_cand["<<st<<
"] "<< X_cand[st] <<
" Xp_cand["<<st<<
"] "<<Xp_cand[st] <<
" iModX "<<iModX[st]<<endl;
2886 if (TotalNumberOfHits >= 4 && Chi2_cand < Chi2_Globcut ){
2889 cout<<
"case2: before writing: "<<endl;
2890 for (Int_t st = 1; st < fNstations; ++st) {
2891 cout<<
" st "<<st<<
" X_cand "<< X_cand[st] <<
" Xp_cand "<< Xp_cand[st] <<endl;
2895 if (
fabs(line_cand[2] * ( Z0_SRC_target - Zcentr ) + line_cand[3]) > half_target_regionY){
2896 if (fDebug) cout<<
" case2:--------half_target_regionY cut: "<<half_target_regionY<<
" Yv "<<
fabs(line_cand[2] * ( Z0_SRC_target - Zcentr ) + line_cand[3])<<endl;
2897 ChiSquareNdf = 999.;
2902 for(
size_t itr = 0; itr< vec_tracks_.size(); ++itr){
2905 for (Int_t st = 1; st < fNstations; ++st) {
2906 cout<<
" itr "<<itr<<
" st "<<st<<
" X_tr "<<vec_tracks_.at(itr).CoordX[st]<<
" Xp_tr "<<vec_tracks_.at(itr).CoordXp[st]<<endl;
2910 ( (vec_tracks_.at(itr).CoordX[1] == X_cand[1] && X_cand[1] > -900.) || (vec_tracks_.at(itr).CoordXp[1] == Xp_cand[1] && Xp_cand[1] > -900.) )
2912 ( (vec_tracks_.at(itr).CoordX[2] == X_cand[2] && X_cand[2] > -900.) || (vec_tracks_.at(itr).CoordXp[2] == Xp_cand[2] && Xp_cand[2] > -900.) )
2914 ( (vec_tracks_.at(itr).CoordX[3] == X_cand[3] && X_cand[3] > -900.) || (vec_tracks_.at(itr).CoordXp[3] == Xp_cand[3] && Xp_cand[3] > -900.) )
2917 if(fDebug ) cout<<
" TotalNumberOfHits "<<TotalNumberOfHits<<
" Nhits["<<itr<<
"] "<<vec_tracks_.at(itr).Nhits<<endl;
2918 if ( vec_tracks_.at(itr).Nhits < TotalNumberOfHits) { vec_tracks_.erase(vec_tracks_.begin()+ itr);
continue;}
2919 if ( vec_tracks_.at(itr).Nhits > TotalNumberOfHits) {
2920 if(fDebug ) cout<<
" Nhits > TotalNumberOfHits . Reject cand "<<endl;
2921 Chi2_cand = 999.;
continue;}
2923 if ( vec_tracks_.at(itr).Nhits == TotalNumberOfHits && TotalNumberOfHits == 4 ){
2924 if(fDebug ) cout<<
" Nhits "<<vec_tracks_.at(itr).Nhits <<
" X1 "<< vec_tracks_.at(itr).CoordX[1]<<
" X2 "<<vec_tracks_.at(itr).CoordX[2] <<endl;
2926 if ( (
fabs(vec_tracks_.at(itr).Xv) + 0.25*
fabs(vec_tracks_.at(itr).Yv) ) > (
fabs(xv_candidate) + 0.25*
fabs(yv_candidate) ) ){
2927 vec_tracks_.erase(vec_tracks_.begin()+ itr);
2931 if(fDebug ) cout<<
" Nhits == TotalNumberOfHits == 4 . Reject cand "<<endl;
2937 if ( vec_tracks_.at(itr).Nhits == TotalNumberOfHits && TotalNumberOfHits == 5 ){
2938 if ( vec_tracks_.at(itr).Chi2 > Chi2_cand ) {vec_tracks_.erase(vec_tracks_.begin()+ itr);
continue;}
2940 if(fDebug ) cout<<
" Nhits == TotalNumberOfHits == 5 . Reject cand ??? "<<endl;
2945 if ( vec_tracks_.at(itr).Nhits > TotalNumberOfHits){
2946 if(fDebug ) cout<<
" Nhits=5 > TotalNumberOfHits=4 . Reject cand ??? "<<endl;
2952 if(fDebug ) cout<<
" Chi2_cand "<<Chi2_cand <<endl;
2953 if ( Chi2_cand > Chi2_Globcut)
continue;
2955 tmpTr.
Chi2 = Chi2_cand ;
2956 tmpTr.
Nhits = TotalNumberOfHits;
2957 tmpTr.
Xv = xv_candidate;
2958 tmpTr.
Yv = yv_candidate;
2959 tmpTr.
param[0] = line_cand[0];
2960 tmpTr.
param[1] = line_cand[1];
2961 tmpTr.
param[2] = line_cand[2];
2962 tmpTr.
param[3] = line_cand[3];
2964 if (fDebug) cout<<
" Xv_candidate "<<xv_candidate<<
" Yv_candidate "<<yv_candidate<<
" ChiNdf "<<Chi2_cand<<endl;
2965 for (Int_t st = 1; st < fNstations; ++st) {
2966 if (Mod_cand[st] < 0)
continue;
2968 tmpTr.
ModNum[st] = iModX[st];
2969 tmpTr.
CoordX[st] = X_cand[st];
2970 tmpTr.
CoordXp[st] = Xp_cand[st];
2971 tmpTr.
SigmaX[st] = SigX_cand[st];
2972 tmpTr.
SigmaXp[st] = SigXp_cand[st];
2973 tmpTr.
CoordZ[st] = ZstnCforglfit[st][Mod_cand[st]];
2974 tmpTr.
CoordY[st] = Y_cand[st];
2975 tmpTr.
SigmaY[st] = SigY_cand[st];
2978 tmpTr.
StripNumX[st] = XClust_[st][iModX[st]][iClustXMod[st]];
2979 tmpTr.
StripNumXp[st]= XpClust_[st][iModX[st]][iClustXMod[st]];
2980 tmpTr.
ClSizeX[st] = XClust_width_[st][iModX[st]][iClustXMod[st]];
2981 tmpTr.
ClSizeXp[st] = XpClust_width_[st][iModX[st]][iClustXMod[st]];
2982 tmpTr.
SumQX[st] = sumQx_[st][iModX[st]][iClustXMod[st]];
2983 tmpTr.
SumQXp[st] = sumQxp_[st][iModX[st]][iClustXMod[st]];
2986 tmpTr.
StripNumX[st] = XspClust_[st][iModX[st]][iClustXMod[st]];
2987 tmpTr.
StripNumXp[st]= XpspClust_[st][iModX[st]][iClustXMod[st]];
2988 tmpTr.
ClSizeX[st] = XspClust_width_[st][iModX[st]][iClustXMod[st]];
2989 tmpTr.
ClSizeXp[st] = XpspClust_width_[st][iModX[st]][iClustXMod[st]];
2990 tmpTr.
SumQX[st] = sumQxsp_[st][iModX[st]][iClustXMod[st]];
2991 tmpTr.
SumQXp[st] = sumQxpsp_[st][iModX[st]][iClustXMod[st]];
2993 if (fDebug) cout<<
" st "<<st<<
" Z "<<tmpTr.
CoordZ[st]<<endl;
2994 if (fDebug) cout<<
" tmpTr.CoordX["<<st <<
"] "<<tmpTr.
CoordX[st]<<
" strip "<<tmpTr.
StripNumX[st]<<
" Clsize "<<tmpTr.
ClSizeX[st]<<
" Q "<<tmpTr.
SumQX[st]<<endl;
2995 if (fDebug) cout<<
" tmpTr.CoordXp["<<st <<
"] "<<tmpTr.
CoordXp[st]<<
" strip "<<tmpTr.
StripNumXp[st]<<
" Clsize "<<tmpTr.
ClSizeXp[st]<<
" Q "<<tmpTr.
SumQXp[st]<<endl;
2998 vec_tracks_.push_back(tmpTr);
2999 if (fDebug) cout<<
" tmpTr.Xv "<<tmpTr.
Xv<<
" tmpTr.Yv "<<tmpTr.
Yv<<endl;
3011void BmnSiliconTrackFinder::Case3( Int_t **Nleftoversp_,Double_t ***leftoverXsp_, Double_t ***leftoverXpsp_, Double_t ***leftoverYsp_,Double_t ***leftoverXsigsp_, Double_t ***leftoverXpsigsp_, Double_t ***leftoverYsigsp_,
3012 Int_t **NleftoverX_, Int_t **NleftoverXp_,Double_t ***leftoverX_, Double_t ***leftoverXp_,Double_t ***leftoverXsig_, Double_t ***leftoverXpsig_,
3013 vector<tracksX> & vec_tracks_,
3014 Double_t ***XClust_, Double_t ***XClust_width_, Double_t ***sumQx_, Double_t ***XpClust_, Double_t ***XpClust_width_, Double_t ***sumQxp_,
3015 Double_t ***XspClust_, Double_t ***XspClust_width_, Double_t ***sumQxsp_, Double_t ***XpspClust_, Double_t ***XpspClust_width_, Double_t ***sumQxpsp_){
3018 if(fDebug) cout<<
" case3 "<<endl;
3020 Bool_t x5_exist = 0;
3022 Double_t TotalNumberOfHits = 0;
3023 Double_t ChiSquareNdf = 999.;
3024 Double_t hrd_case3 = 0.7;
3026 for (Int_t imod1 = 0; imod1 < fNmodules1; imod1++) {
3027 for (Int_t cl1 = 0; cl1 < Nleftoversp_[1][imod1]; cl1++){
3029 for (Int_t imod2 = 0; imod2 < fNmodules2; imod2++) {
3030 for (Int_t cl2 = 0; cl2 < Nleftoversp_[2][imod2]; cl2++){
3032 Bool_t was_candidate = 0;
3033 Int_t iClustXMod_candidate = -1, iModX_candidate = -1;
3034 Double_t ax12_from_target = 999., x_ex3 = 999., ay12_from_target = 999., y_ex3 = 999.;
3035 Double_t ay1_2_cand = 999;
3036 Double_t bx_candidate = 999., by_candidate = 999.;
3038 for ( Int_t
i = 1;
i < fNstations; ++
i){
3043 if (
fabs(leftoverXsp_[1][imod1][cl1] - leftoverXsp_[2][imod2][cl2] ) > half_roadX1_X2 )
continue;
3044 if (
fabs(leftoverXsp_[1][imod1][cl1] - leftoverXpsp_[1][imod1][cl1]) > half_roadX1_Xp1 )
continue;
3045 if (
fabs(leftoverXsp_[2][imod2][cl2] - leftoverXpsp_[2][imod2][cl2]) > half_roadX2_Xp2 )
continue;
3046 if (
fabs(leftoverYsp_[1][imod1][cl1] - leftoverYsp_[2][imod2][cl2] ) > half_roadY1_Y2)
continue;
3053 if ( fDebug) cout<<
"case3: iModX[1] "<<iModX[1]<<
" iCl "<<iClustXMod[1]<<
" iModX[2] "<<iModX[2]<<
" iCl "<<iClustXMod[2]<<endl;
3055 Double_t z1_2 = 0.5 * ( Zstation[1][imod1] + Zstation[2][imod2] );
3056 Double_t y1_2 = 0.5 * ( leftoverYsp_[1][imod1][cl1] +leftoverYsp_[2][imod2][cl2] );
3057 ax12_from_target = (Xv_av - leftoverXsp_[1][imod1][cl1] )/( Z0_SRC_target - Zstation[1][imod1] );
3058 ay12_from_target = (Yv_av - y1_2 )/( Z0_SRC_target - z1_2 );
3060 if ( fDebug) cout<<
" ax12_from_target "<<ax12_from_target<<
" Xv_av "<<Xv_av<<
" Yv_av "<<Yv_av<<endl;
3062 for (Int_t ist = 1; ist < fNstations; ++ist) {
3063 for (Int_t imod = 0; imod < fNmodules; imod++) {
3064 Xforglfit[ist][imod] = -999.;
3065 Yforglfit[ist][imod] = -999.;
3066 Xpforglfit[ist][imod] = -999.;
3067 SigmYforglfit[ist][imod] = -999.;
3068 SigmXforglfit[ist][imod] = -999.;
3069 SigmXpforglfit[ist][imod] = -999.;
3074 for (Int_t ist = 1; ist < fNstations; ++ist) {
3075 if ( ist > 2 )
continue;
3076 Xforglfit[ist][iModX[ist]] = leftoverXsp_[ist][iModX[ist]][iClustXMod[ist]];
3077 SigmXforglfit[ist][iModX[ist]] = leftoverXsigsp_[ist][iModX[ist]][iClustXMod[ist]];
3078 Xpforglfit[ist][iModX[ist]] = leftoverXpsp_[ist][iModX[ist]][iClustXMod[ist]];
3079 SigmXpforglfit[ist][iModX[ist]] = leftoverXpsigsp_[ist][iModX[ist]][iClustXMod[ist]];
3080 Yforglfit[ist][iModX[ist]] = leftoverYsp_[ist][iModX[ist]][iClustXMod[ist]];
3081 SigmYforglfit[ist][iModX[ist]] = leftoverYsigsp_[ist][iModX[ist]][iClustXMod[ist]];
3084 for (Int_t imod3 = 0; imod3 < fNmodules3; imod3++) {
3085 for (Int_t clx3 = 0; clx3 < NleftoverX_[3][imod3]; clx3++){
3087 x_ex3 = ax12_from_target*(Zstation[3][imod3] - Z0_SRC_target ) + Xv_av;
3088 y_ex3 = ay12_from_target*(Zstation[3][imod3] - Z0_SRC_target ) + Yv_av;
3090 if ( fDebug) cout<<
" X1 "<<leftoverXsp_[1][imod1][cl1]<<
" X2 "<<leftoverXsp_[2][imod2][cl2]<<
" x_ex3 "<<x_ex3<<endl;
3091 if ( fDebug) cout<<
" Y1 "<<leftoverYsp_[1][imod1][cl1]<<
" Y2 "<<leftoverYsp_[2][imod2][cl2]<<
" y_ex3 "<<y_ex3<<endl;
3093 if (
fabs(x_ex3 - leftoverX_[3][imod3][clx3]) > hrd_case3 )
continue;
3095 if ( !was_candidate ) {
3097 iModX_candidate = imod3;
3098 iClustXMod_candidate = clx3;
3102 ay1_2_cand = ay12_from_target;
3103 if (fDebug) cout<<
" write candidate3 x3 = "<<leftoverX_[3][imod3][clx3]<<endl;
3105 if ( (x_ex3 - leftoverX_[3][imod3][clx3]) < (x_ex3 - leftoverX_[3][iModX_candidate][iClustXMod_candidate]) ) {
3106 iModX_candidate = imod3;
3107 iClustXMod_candidate = clx3;
3109 ay1_2_cand = ay12_from_target;
3110 if (fDebug) cout<<
" re-write cand3 x3 = "<<leftoverX_[3][imod3][clx3]<<endl;
3114 if ( was_candidate ) Npoint = 5;
3119 if ( was_candidate ){
3120 iModX[3] = iModX_candidate;
3121 iClustXMod[3] = iClustXMod_candidate;
3122 Xforglfit[3][iModX_candidate] = leftoverX_[3][iModX_candidate][iClustXMod_candidate];
3123 SigmXforglfit[3][iModX_candidate] = leftoverXsig[3][iModX_candidate][iClustXMod_candidate];
3124 Xpforglfit[3][iModX_candidate] = -999.;
3125 SigmXpforglfit[3][iModX_candidate] = -999.;
3126 Yforglfit[3][iModX_candidate] = -999.;
3127 SigmYforglfit[3][iModX_candidate] = -999.;
3131 for (Int_t coeff = 0; coeff < 4; ++coeff) {
3135 calculationOfChi2(Xforglfit, Xpforglfit, SigmXforglfit, SigmXpforglfit, iClustXMod, iModX, ChiSquare, line);
3137 TotalNumberOfHits = Npoint;
3139 ChiSquareNdf = ChiSquare/(TotalNumberOfHits - 4);
3140 else ChiSquareNdf = ChiSquare;
3142 if ( fDebug) cout<<
" case3 x: before cut ChiSquareNDF "<<ChiSquareNdf<<
" TotalNumberOfHits "<<TotalNumberOfHits
3143 <<
" fitted line par "<<line[0]<<
" "<<line[1]<<
" "<<line[2]<<
" "<<line[3]<<endl;
3144 if (fDebug) cout<<endl;
3147 if ( ChiSquareNdf < Chi2_Globcut ){
3150 line[2] = ay1_2_cand;
3151 line[3] = Yv_av - line[2] * ( Z0_SRC_target - Zcentr );
3152 if ( fDebug) cout<<
" replace to Yv-Y12 line par "<<line[2]<<
" "<<line[3]<<endl;
3158 Xforglfit[3][iModX[3]] = -999.;
3159 SigmXforglfit[3][iModX[3]] = -999.;
3160 Xpforglfit[3][iModX[3]] = -999.;
3161 SigmXpforglfit[3][iModX[3]] = -999.;
3162 iClustXMod_candidate = -1;
3165 if ( fDebug) cout<<
"---x.----"<<endl;
3166 if ( fDebug) cout<<
"case3: end 3x st: iModX[1] "<<iModX[1]<<
" iCl "<<iClustXMod[1]<<
" iModX[2] "<<iModX[2]<<
" iCl "<<iClustXMod[2]<<endl;
3173 for (Int_t imod3 = 0; imod3 < fNmodules3; imod3++) {
3174 for (Int_t clxp3 = 0; clxp3 < NleftoverXp_[3][imod3]; clxp3++){
3178 Xforglfit[3][iModX[3]] = -999.;
3179 SigmXforglfit[3][iModX[3]] = -999.;
3180 Xpforglfit[3][iModX[3]] = -999.;
3181 SigmXpforglfit[3][iModX[3]] = -999.;
3183 x_ex3 = ax12_from_target*(Zstation[3][imod3] - Z0_SRC_target ) + Xv_av;
3184 y_ex3 = ay12_from_target*(Zstation[3][imod3] - Z0_SRC_target ) + Yv_av;
3186 if ( fDebug) cout<<
" X1 "<<leftoverXsp_[1][imod1][cl1]<<
" X2 "<<leftoverXsp_[2][imod2][cl2]<<
" x_ex3 "<<x_ex3<<endl;
3187 if ( fDebug) cout<<
" Y1 "<<leftoverYsp_[1][imod1][cl1]<<
" Y2 "<<leftoverYsp_[2][imod2][cl2]<<
" y_ex3 "<<y_ex3<<endl;
3190 Double_t xp_ex3 = y_ex3*(Angle(st3,imod3))*tg2_5 + x_ex3;
3191 if (
fabs(xp_ex3 - leftoverXp_[3][imod3][clxp3]) > hrd_case3 )
continue;
3193 if ( !was_candidate ) {
3194 iModX_candidate = imod3;
3195 iClustXMod_candidate = clxp3;
3197 if (fDebug) cout<<
" write candidate3 xp3 = "<<leftoverXp_[3][imod3][clxp3]<<
" xp_ex3 "<<xp_ex3<<endl;
3199 if (
fabs(xp_ex3 - leftoverXp_[3][imod3][clxp3]) <
fabs(xp_ex3 - leftoverXp_[3][iModX_candidate][iClustXMod_candidate]) ) {
3201 iModX_candidate = imod3;
3202 iClustXMod_candidate = clxp3;
3203 if (fDebug) cout<<
" re-write cand3 xp3 = "<<leftoverXp_[3][imod3][clxp3]<<
" xp_ex3 "<<xp_ex3<<endl;
3208 if (was_candidate) Npoint = 5;
3213 if ( was_candidate ){
3214 iModX[3] = iModX_candidate;
3215 iClustXMod[3] = iClustXMod_candidate;
3216 Xforglfit[3][iModX_candidate] = -999.;
3217 SigmXforglfit[3][iModX_candidate] = -999.;
3218 Xpforglfit[3][iModX_candidate] = leftoverXp_[3][iModX_candidate][iClustXMod_candidate];
3219 SigmXpforglfit[3][iModX_candidate] = leftoverXpsig[3][iModX_candidate][iClustXMod_candidate];
3220 Yforglfit[3][iModX_candidate] = -999.;
3221 SigmYforglfit[3][iModX_candidate] = -999.;
3224 TotalNumberOfHits = 0;
3225 ChiSquareNdf = 999.;
3228 for (Int_t coeff = 0; coeff < 4; ++coeff) {
3232 calculationOfChi2(Xforglfit, Xpforglfit, SigmXforglfit, SigmXpforglfit, iClustXMod, iModX, ChiSquare, line);
3234 TotalNumberOfHits = Npoint;
3236 ChiSquareNdf = ChiSquare/(TotalNumberOfHits - 4);
3237 else ChiSquareNdf = ChiSquare;
3239 if ( fDebug) cout<<
" case3 xp: before cut ChiSquareNDF "<<ChiSquareNdf<<
" TotalNumberOfHits "<<TotalNumberOfHits <<
" fitted line par "<<line[0]<<
" "<<line[1]<<
" "<<line[2]<<
" "<<line[3]<<endl;
3240 if (fDebug) cout<<endl;
3241 if ( fDebug) cout<<
" ChiSquareNdf "<<ChiSquareNdf<<endl;
3244 if ( fDebug) cout<<
" ChiSquareNdf "<<ChiSquareNdf<<endl;
3246 if ( ChiSquareNdf < Chi2_Globcut ){
3248 if (
fabs(line[2] * ( Z0_SRC_target - Zcentr ) + line[3]) > half_target_regionY ){
3249 line[2] = ay12_from_target;
3250 line[3] = Yv_av - line[2] * ( Z0_SRC_target - Zcentr );
3251 if ( fDebug) cout<<
" xp replace to Yv-Y12 line par "<<line[2]<<
" "<<line[3]<<endl;
3255 if ( fDebug) cout<<
"xp5_exist"<<endl;
3262 Xforglfit[3][iModX[3]] = -999.;
3263 SigmXforglfit[3][iModX[3]] = -999.;
3264 Xpforglfit[3][iModX[3]] = -999.;
3265 SigmXpforglfit[3][iModX[3]] = -999.;
3270 if (fDebug) cout<<
"---xp.----"<<endl;
3271 if (fDebug) cout<<
" was_candidate "<<was_candidate<<endl;
3272 if (fDebug) cout<<
"case3: end 3xp st: iModX[1] "<<iModX[1]<<
" iCl "<<iClustXMod[1]<<
" iModX[2] "<<iModX[2]<<
" iCl "<<iClustXMod[2]<<endl;
3273 if (fDebug) cout<<
" TotalNumberOfHits "<<TotalNumberOfHits<<
" ChiSquareNdf "<<ChiSquareNdf<<endl;
3276 Int_t Mod_cand[fNstations], iCl_cand[fNstations];
3277 Double_t X_cand[fNstations], Xp_cand[fNstations], Y_cand[fNstations],
3278 SigX_cand[fNstations], SigXp_cand[fNstations], SigY_cand[fNstations],
3279 Chi2_cand = 999., line_cand[4];
3281 for (Int_t coeff = 0; coeff < 4; ++coeff) {
3282 line_cand[coeff] = 0.;
3285 for (Int_t st = 1; st < fNstations; ++st) {
3289 Xp_cand[st] = -999.;
3291 SigX_cand[st] = -999.;
3292 SigXp_cand[st] = -999.;
3293 SigY_cand[st] = -999.;
3296 if (fDebug) cout<<
" TotalNumberOfHits "<<TotalNumberOfHits<<
" ChiSquareNdf "<<ChiSquareNdf<<
" Chi2_Globcut "<<Chi2_Globcut<<endl;
3299 if (TotalNumberOfHits > 4){
3300 if (
fabs(line[0] * ( Z0_SRC_target - Zcentr ) + line[1]) > half_target_regionX){
3301 if (fDebug) cout<<
" half_target_regionX cut: "<<half_target_regionX<<
" Xv "<<
fabs(line[0] * ( Z0_SRC_target - Zcentr ) + line[1])<<endl;
3302 ChiSquareNdf = 999.;
3307 if (TotalNumberOfHits > 4 && ChiSquareNdf < Chi2_cand && ChiSquareNdf < Chi2_Globcut && was_candidate ){
3308 if (fDebug) cout<<
" TotalNumberOfHits "<<TotalNumberOfHits<<
" ChiSquareNdf "<<ChiSquareNdf<<endl;
3309 Chi2_cand = ChiSquareNdf;
3310 line_cand[0] = line[0];
3311 line_cand[1] = line[1];
3312 line_cand[2] = line[2];
3313 line_cand[3] = line[3];
3315 for (Int_t st = 1; st < fNstations; ++st) {
3316 if (iModX[st] < 0)
continue;
3318 Mod_cand[st] = iModX[st];
3319 iCl_cand[st] = iClustXMod[st];
3320 X_cand[st] = Xforglfit[st][iModX[st]];
3321 Y_cand[st] = Yforglfit[st][iModX[st]];
3322 Xp_cand[st] = Xpforglfit[st][iModX[st]];
3323 SigX_cand[st] = SigmXforglfit[st][iModX[st]];
3324 SigY_cand[st] = SigmYforglfit[st][iModX[st]];
3325 SigXp_cand[st] = SigmXpforglfit[st][iModX[st]];
3327 if (fDebug) cout<<
" X_cand["<<st<<
"] "<< X_cand[st] <<
" Xp_cand["<<st<<
"] "<<Xp_cand[st] <<
" iMod "<<iModX[st]<<
" iCl "<<iClustXMod[st]<<
" ClSizeXp "<<XpspClust_width_[st][Mod_cand[st]][iClustXMod[st]]<<endl;
3332 if (TotalNumberOfHits > 4 && Chi2_cand < Chi2_Globcut && was_candidate ){
3334 for(
size_t itr = 0; itr< vec_tracks_.size(); ++itr){
3337 ( (vec_tracks_.at(itr).CoordX[1] == X_cand[1] && X_cand[1] > -900.) || (vec_tracks_.at(itr).CoordXp[1] == Xp_cand[1] && Xp_cand[1] > -900.) )
3339 ( (vec_tracks_.at(itr).CoordX[2] == X_cand[2] && X_cand[2] > -900.) || (vec_tracks_.at(itr).CoordXp[2] == Xp_cand[2] && Xp_cand[2] > -900.) )
3341 ( (vec_tracks_.at(itr).CoordX[3] == X_cand[3] && X_cand[3] > -900.) || (vec_tracks_.at(itr).CoordXp[3] == Xp_cand[3] && Xp_cand[3] > -900.) )
3344 if ( vec_tracks_.at(itr).Chi2 > Chi2_cand ){
3345 vec_tracks_.erase(vec_tracks_.begin()+ itr);
3354 if ( Chi2_cand > Chi2_Globcut)
continue;
3356 tmpTr.
Chi2 = Chi2_cand;
3357 tmpTr.
Nhits = TotalNumberOfHits;
3358 bx_candidate = line_cand[0] * ( Z0_SRC_target -Zcentr) + line_cand[1];
3359 by_candidate = line_cand[2] * ( Z0_SRC_target -Zcentr) + line_cand[3];
3360 tmpTr.
Xv = bx_candidate;
3361 tmpTr.
Yv = by_candidate;
3362 tmpTr.
param[0] = line_cand[0];
3363 tmpTr.
param[1] = line_cand[1];
3364 tmpTr.
param[2] = line_cand[2];
3365 tmpTr.
param[3] = line_cand[3];
3367 if (fDebug) cout<<
" Xv_candidate "<<bx_candidate<<
" Yv_candidate "<<by_candidate<<
" ChiSquareNdf "<<ChiSquareNdf<<endl;
3368 for (Int_t st = 1; st < fNstations; ++st) {
3369 if (Mod_cand[st] < 0)
continue;
3371 tmpTr.
ModNum[st] = iModX[st];
3372 tmpTr.
CoordX[st] = X_cand[st];
3373 tmpTr.
CoordXp[st] = Xp_cand[st];
3374 tmpTr.
SigmaX[st] = SigX_cand[st];
3375 tmpTr.
SigmaXp[st] = SigXp_cand[st];
3376 tmpTr.
CoordZ[st] = ZstnCforglfit[st][Mod_cand[st]];
3377 tmpTr.
CoordY[st] = Y_cand[st];
3378 tmpTr.
SigmaY[st] = SigY_cand[st];
3381 tmpTr.
StripNumX[st] = XspClust_[st][Mod_cand[st]][iCl_cand[st]];
3382 tmpTr.
StripNumXp[st]= XpspClust_[st][Mod_cand[st]][iCl_cand[st]];
3383 tmpTr.
ClSizeX[st] = XspClust_width_[st][Mod_cand[st]][iCl_cand[st]];
3384 tmpTr.
ClSizeXp[st] = XpspClust_width_[st][Mod_cand[st]][iCl_cand[st]];
3385 tmpTr.
SumQX[st] = sumQxsp_[st][Mod_cand[st]][iCl_cand[st]];
3386 tmpTr.
SumQXp[st] = sumQxpsp_[st][Mod_cand[st]][iCl_cand[st]];
3388 tmpTr.
StripNumX[st] = XClust_[st][iModX[st]][iClustXMod_candidate];
3389 tmpTr.
StripNumXp[st]= XpClust_[st][iModX[st]][iClustXMod_candidate];
3390 tmpTr.
ClSizeX[st] = XClust_width_[st][iModX[st]][iClustXMod_candidate];
3391 tmpTr.
ClSizeXp[st] = XpClust_width_[st][iModX[st]][iClustXMod_candidate];
3392 tmpTr.
SumQX[st] = sumQx_[st][iModX[st]][iClustXMod_candidate];
3393 tmpTr.
SumQXp[st] = sumQxp_[st][iModX[st]][iClustXMod_candidate];
3396 if (fDebug) cout<<
" st "<<st<<
" Z "<<tmpTr.
CoordZ[st]<<endl;
3397 if (fDebug) cout<<
" tmpTr.CoordX["<<st <<
"] "<<tmpTr.
CoordX[st]<<
" strip "<<tmpTr.
StripNumX[st]<<
" Clsize "<<tmpTr.
ClSizeX[st]<<
" Q "<<tmpTr.
SumQX[st]<<endl;
3398 if (fDebug) cout<<
" tmpTr.CoordXp["<<st <<
"] "<<tmpTr.
CoordXp[st]<<
" strip "<<tmpTr.
StripNumXp[st]<<
" Clsize "<<tmpTr.
ClSizeXp[st]<<
" Q "<<tmpTr.
SumQXp[st]<<endl;
3401 vec_tracks_.push_back(tmpTr);
3408 if (fDebug) cout<<endl;
3409 if (fDebug) cout<<
" after second case 3 write cand: "<<vec_tracks_.size()<<endl;
3416void BmnSiliconTrackFinder::Case4( Int_t & Fi_St,Int_t **Nleftoversp_,Double_t ***leftoverXsp_, Double_t ***leftoverXpsp_, Double_t ***leftoverYsp_,Double_t ***leftoverXsigsp_, Double_t ***leftoverXpsigsp_, Double_t ***leftoverYsigsp_,
3417 Int_t **NleftoverX_, Int_t **NleftoverXp_,Double_t ***leftoverX_, Double_t ***leftoverXp_,Double_t ***leftoverXsig_, Double_t ***leftoverXpsig_,
3418 vector<tracksX> & vec_tracks_,
3419 Double_t ***XClust_, Double_t ***XClust_width_, Double_t ***sumQx_, Double_t ***XpClust_, Double_t ***XpClust_width_, Double_t ***sumQxp_,
3420 Double_t ***XspClust_, Double_t ***XspClust_width_, Double_t ***sumQxsp_, Double_t ***XpspClust_, Double_t ***XpspClust_width_, Double_t ***sumQxpsp_){
3424 Int_t Nmod_Fi = fNmodules1;
3425 Int_t Nmod_Sec = fNmodules2;
3429 Nmod_Fi = fNmodules2;
3430 Nmod_Sec = fNmodules1;
3434 Int_t Mod_cand[fNstations], iCl_cand[fNstations];
3435 Double_t X_cand[fNstations], Xp_cand[fNstations], Y_cand[fNstations],
3436 SigX_cand[fNstations], SigXp_cand[fNstations], SigY_cand[fNstations],
3437 Chi2_cand = 999., line_cand[4];
3439 for (Int_t coeff = 0; coeff < 4; ++coeff) {
3440 line_cand[coeff] = 0.;
3443 for (Int_t st = 1; st < fNstations; ++st) {
3447 Xp_cand[st] = -999.;
3449 SigX_cand[st] = -999.;
3450 SigXp_cand[st] = -999.;
3451 SigY_cand[st] = -999.;
3454 Double_t bx_candidate = 999., by_candidate = 999.;
3457 for (Int_t imod3 = 0; imod3 < fNmodules3; imod3++) {
3458 for (Int_t cl3 = 0; cl3 < Nleftoversp_[3][imod3]; cl3++){
3460 Double_t Xv_min = 2.5, Yv_min = 2.5;
3463 for (Int_t imod1 = 0; imod1 < Nmod_Fi; imod1++) {
3464 for (Int_t cl1 = 0; cl1 < NleftoverX[Fi_St][imod1]; cl1++){
3466 Double_t ax_target_region = -999., ay_target_region = -999.;
3469 for (Int_t imod2 = 0; imod2 < Nmod_Sec; imod2++) {
3470 for (Int_t cl2 = 0; cl2 < NleftoverXp_[Sec_St][imod2]; cl2++){
3472 if (
fabs(leftoverX_[Fi_St][imod1][cl1] - leftoverXp_[Sec_St][imod2][cl2]) > half_roadX1_Xp2 )
continue;
3474 Double_t z1_2 = 0.5 * ( Zstation[1][imod1] + Zstation[2][imod2] );
3475 Double_t Y_1_2 = GetXYspatial_station1_2(Fi_St, imod2, leftoverX_[Fi_St][imod1][cl1] , leftoverXp_[Sec_St][imod2][cl2]);
3477 if(fDebug) cout<<
" Fi_St "<<Fi_St<<
" X "<<leftoverX_[Fi_St][imod1][cl1]<<
" Xp "<<leftoverXp_[Sec_St][imod2][cl2]<<
" new Y "<<Y_1_2<<endl;
3479 ax_target_region = (leftoverX_[Fi_St][imod1][cl1] - leftoverXsp_[3][imod3][cl3]) / (Zstation[Fi_St][imod1] - Zstation[3][imod3]);
3480 ay_target_region = (Y_1_2 - leftoverYsp_[3][imod3][cl3]) / (z1_2 - Zstation[3][imod3]);
3482 Double_t Xv_cand = ax_target_region * (- z1_2 + Z0_SRC_target )+ leftoverX_[Fi_St][imod1][cl1];
3483 Double_t Yv_cand = ay_target_region * ( - z1_2 + Z0_SRC_target )+ Y_1_2;
3485 if (
fabs(Xv_cand) > half_target_regionX)
continue;
3486 if (
fabs(Yv_cand) > half_target_regionY)
continue;
3488 if ( (
fabs(Xv_min) + 0.25*
fabs(Yv_min) ) < (
fabs(Xv_cand) + 0.25*
fabs(Yv_cand) ) )
continue;
3490 if(fDebug) cout<<
" Xv_cand "<<Xv_cand<<
" Yv_cand "<<Yv_cand<<endl;
3495 bx_candidate = Xv_cand;
3496 by_candidate = Yv_cand;
3498 line_cand[0] = ax_target_region;
3499 line_cand[1] = Xv_cand + ax_target_region * ( Zcentr - Z0_SRC_target );
3500 line_cand[2] = ay_target_region;
3501 line_cand[3] = Yv_cand + ay_target_region * ( Zcentr - Z0_SRC_target );
3503 Mod_cand[Fi_St] = imod1;
3504 iCl_cand[Fi_St] = cl1;
3505 X_cand[Fi_St] = leftoverX_[Fi_St][imod1][cl1];
3506 Xp_cand[Fi_St] = -999.;
3507 Y_cand[Fi_St] = -999.;
3508 SigX_cand[Fi_St] = 1.;
3509 SigXp_cand[Fi_St] = 1.;
3510 SigY_cand[Fi_St] = 1.;
3512 Mod_cand[Sec_St] = imod2;
3513 iCl_cand[Sec_St] = cl2;
3514 X_cand[Sec_St] = -999.;
3515 Xp_cand[Sec_St] = leftoverXp_[Sec_St][imod2][cl2];
3516 Y_cand[Sec_St] = -999.;
3517 SigX_cand[Sec_St] = 1.;
3518 SigXp_cand[Sec_St] = 1.;
3519 SigY_cand[Sec_St] = 1.;
3521 Mod_cand[3] = imod3;
3523 X_cand[3] = leftoverXsp_[3][imod3][cl3];
3524 Xp_cand[3] = leftoverXpsp_[3][imod3][cl3];
3525 Y_cand[3] = leftoverYsp_[3][imod3][cl3];
3526 SigX_cand[3] = leftoverXsigsp_[3][imod3][cl3];
3527 SigXp_cand[3] = leftoverXpsigsp_[3][imod3][cl3];
3528 SigY_cand[3] = leftoverYsigsp_[3][imod3][cl3];
3534 for(
size_t itr = 0; itr< vec_tracks.size(); ++itr){
3537 ( (vec_tracks.at(itr).CoordX[1] == X_cand[1] && X_cand[1] > -900.) || (vec_tracks.at(itr).CoordXp[1] == Xp_cand[1] && Xp_cand[1] > -900.) )
3539 ( (vec_tracks.at(itr).CoordX[2] == X_cand[2] && X_cand[2] > -900.) || (vec_tracks.at(itr).CoordXp[2] == Xp_cand[2] && Xp_cand[2] > -900.) )
3541 ( (vec_tracks.at(itr).CoordX[3] == X_cand[3] && X_cand[3] > -900.) || (vec_tracks.at(itr).CoordXp[3] == Xp_cand[3] && Xp_cand[3] > -900.) )
3543 if ( (
fabs(vec_tracks.at(itr).Xv) + 0.25*
fabs(vec_tracks.at(itr).Yv) ) > (
fabs(bx_candidate) + 0.25*
fabs(by_candidate) ) ){
3544 vec_tracks.erase(vec_tracks.begin()+ itr);
3554 if ( Chi2_cand > Chi2_Globcut)
continue;
3556 tmpTr.
Chi2 = Chi2_cand ;
3558 tmpTr.
Xv = bx_candidate;
3559 tmpTr.
Yv = by_candidate;
3560 tmpTr.
param[0] = line_cand[0];
3561 tmpTr.
param[1] = line_cand[1];
3562 tmpTr.
param[2] = line_cand[2];
3563 tmpTr.
param[3] = line_cand[3];
3565 if (fDebug) cout<<endl;
3566 if (fDebug) cout<<
" Xv_candidate "<<bx_candidate<<
" Yv_candidate "<<by_candidate<<
" Nhits "<<tmpTr.
Nhits<<endl;
3567 for (Int_t st = 1; st < fNstations; ++st) {
3568 if (Mod_cand[st] < 0)
continue;
3570 tmpTr.
ModNum[st] = iModX[st];
3571 tmpTr.
CoordX[st] = X_cand[st];
3572 tmpTr.
CoordXp[st] = Xp_cand[st];
3573 tmpTr.
SigmaX[st] = SigX_cand[st];
3574 tmpTr.
SigmaXp[st] = SigXp_cand[st];
3575 tmpTr.
CoordZ[st] = ZstnCforglfit[st][Mod_cand[st]];
3576 tmpTr.
CoordY[st] = Y_cand[st];
3577 tmpTr.
SigmaY[st] = SigY_cand[st];
3580 tmpTr.
StripNumX[st] = XClust_[st][Mod_cand[st]][iCl_cand[st]];
3581 tmpTr.
StripNumXp[st]= XpClust_[st][Mod_cand[st]][iCl_cand[st]];
3582 tmpTr.
ClSizeX[st] = XClust_width_[st][Mod_cand[st]][iCl_cand[st]];
3583 tmpTr.
ClSizeXp[st] = XpClust_width_[st][Mod_cand[st]][iCl_cand[st]];
3584 tmpTr.
SumQX[st] = sumQx_[st][Mod_cand[st]][iCl_cand[st]];
3585 tmpTr.
SumQXp[st] = sumQxp_[st][Mod_cand[st]][iCl_cand[st]];
3587 tmpTr.
StripNumX[st] = XspClust_[st][Mod_cand[st]][iCl_cand[st]];
3588 tmpTr.
StripNumXp[st]= XpspClust_[st][Mod_cand[st]][iCl_cand[st]];
3589 tmpTr.
ClSizeX[st] = XspClust_width_[st][Mod_cand[st]][iCl_cand[st]];
3590 tmpTr.
ClSizeXp[st] = XpspClust_width_[st][Mod_cand[st]][iCl_cand[st]];
3591 tmpTr.
SumQX[st] = sumQxsp_[st][Mod_cand[st]][iCl_cand[st]];
3592 tmpTr.
SumQXp[st] = sumQxpsp_[st][Mod_cand[st]][iCl_cand[st]];
3595 if (fDebug) cout<<
" st "<<st<<
" Z "<<tmpTr.
CoordZ[st]<<endl;
3596 if (fDebug) cout<<
" tmpTr.CoordX["<<st <<
"] "<<tmpTr.
CoordX[st]<<
" strip "<<tmpTr.
StripNumX[st]<<
" Clsize "<<tmpTr.
ClSizeX[st]<<
" Q "<<tmpTr.
SumQX[st]<<endl;
3597 if (fDebug) cout<<
" tmpTr.CoordXp["<<st <<
"] "<<tmpTr.
CoordXp[st]<<
" strip "<<tmpTr.
StripNumXp[st]<<
" Clsize "<<tmpTr.
ClSizeXp[st]<<
" Q "<<tmpTr.
SumQXp[st]<<endl;
3601 vec_tracks.push_back(tmpTr);
3602 if (fDebug) cout<<
" tmpTr.Xv "<<tmpTr.
Xv<<
" tmpTr.Yv "<<tmpTr.
Yv<<endl;
3603 if (fDebug) cout<<endl;
3613Double_t BmnSiliconTrackFinder::GetXYspatial_station1_2(Int_t &Fi_St, Int_t &imod, Double_t & XCoor, Double_t &XpCoor ) {
3618 Double_t Ysp_1_2 = -999.0;
3619 Double_t YCoor_cand;
3627 if ( XCoor > -900. && XpCoor > -900.){
3628 YCoor_cand = ( XpCoor - XCoor ) / tg2_5;
3630 if (YCoor_cand >= shiftY[Sec_St][imod] - 0.5 && YCoor_cand <= shiftY[Sec_St][imod] + 12. +0.5 ) {
3631 Ysp_1_2 = YCoor_cand;
3638void BmnSiliconTrackFinder::calculationOfChi2(Double_t **Xforglfit_, Double_t **Xpforglfit_, Double_t **SigmXforglfit_, Double_t **SigmXpforglfit_,
3639 Int_t *iClustXMod_, Int_t *iModX_, Double_t & ChiSquare_, Double_t *line_){
3643 if (fDebug) cout<<endl;
3644 if (fDebug) cout<<
"--GlobalFit--"<<endl;
3647 for (Int_t
i = 0;
i < 4; ++
i) {
3648 for (Int_t j = 0; j < 4; ++j) {
3652 for (Int_t
i = 0;
i < 4; ++
i) {
3656 for (Int_t
i = 0;
i < 16; ++
i) {
3659 for (Int_t
i = 0;
i < 16; ++
i) {
3660 AmatrInverted[
i] = 0;
3663 GlobalFit( Xforglfit_, Xpforglfit_, ZstnCforglfit, SigmXforglfit_, SigmXpforglfit_, Amatr, rhs);
3667 for (Int_t row = 0; row < 4; ++row) {
3668 for (Int_t col = 0; col < 4; ++col) {
3669 AmatrArray[4*col + row] = Amatr[row][col];
3673 bool ifSuccess = InvertMatrix(AmatrArray, AmatrInverted);
3675 if (!ifSuccess) {cout<<
"???? InvertMatrix is unsuccessfull ????"<<endl; }
3676 if (!ifSuccess) {cout<<
"???? InvertMatrix is unsuccessfull ????"<<endl;
return;}
3678 for (Int_t coeff = 0; coeff < 4; ++coeff) {
3679 for (Int_t iCol = 0; iCol < 4; ++iCol) {
3680 line_[coeff] += AmatrInverted[4*iCol + coeff]*rhs[iCol];
3684 if ( fDebug ) cout<<
" ax "<<line_[0]<<
" bx "<<line_[1]<<
" ay "<<line_[2]<<
" by "<<line_[3]<<endl;
3687 if ( fDebug) cout<<
"ChiSquare_ "<<ChiSquare_<<endl;
3688 for (Int_t ist = 1; ist < fNstations; ++ist) {
3689 if (Xforglfit_[ist][iModX_[ist]] < -900. )
continue;
3691 Double_t hit = Xforglfit_[ist][iModX_[ist]] ;
3692 Double_t fit = line_[0]*ZstnCforglfit[ist][iModX_[ist]] + line_[1];
3694 Double_t sigmaX_2 = SigmXforglfit_[ist][iModX_[ist]]*SigmXforglfit_[ist][iModX_[ist]];
3695 ChiSquare_ += ((hit-fit)*(hit-fit))/(sigmaX_2);
3696 if ( fDebug) cout<<
" ist "<<ist<<
" hitx "<<hit<<
" dx "<<hit-fit<<
" Z "<<ZstnCforglfit[ist][iModX_[ist]]<<
" sigma "<<SigmXforglfit_[ist][iModX_[ist]]<<
" ChiSquare_ "<<ChiSquare_<<endl;
3698 for (Int_t ist = 1; ist < fNstations; ++ist) {
3699 if (Xpforglfit_[ist][iModX_[ist]] < -900. )
continue;
3701 Double_t hit = Xpforglfit_[ist][iModX_[ist]];
3702 Double_t fit = (line_[0]*ZstnCforglfit[ist][iModX_[ist]] + line_[1]) + (line_[2]*ZstnCforglfit[ist][iModX_[ist]] + line_[3]) * Angle(ist,iModX_[ist]) * tg2_5 ;
3704 Double_t sigmaXp_2 = SigmXpforglfit_[ist][iModX_[ist]]*SigmXpforglfit_[ist][iModX_[ist]];
3705 ChiSquare_ += ((hit-fit)*(hit-fit))/(sigmaXp_2);
3706 if ( fDebug) cout<<
" ist "<<ist<<
" hitxp "<<hit<<
" dxp "<<hit-fit<<
" Z "<<ZstnCforglfit[ist][iModX_[ist]]<<
" sigma "<<SigmXpforglfit_[ist][iModX_[ist]]<<
" ChiSquare_ "<<ChiSquare_<<endl;
3715Int_t BmnSiliconTrackFinder::Angle(Int_t &st, Int_t &Mod) {
3718 if (st == 1 && Mod == 0) sign = -1;
3719 if (st == 1 && Mod == 1) sign = 1;
3720 if (st == 1 && Mod == 2) sign = 1;
3721 if (st == 1 && Mod == 3) sign = -1;
3723 if (st == 2 && Mod == 0) sign = -1;
3724 if (st == 2 && Mod == 1) sign = 1;
3726 if (st == 3 && Mod == 0) sign = -1;
3727 if (st == 3 && Mod == 1) sign = -1;
3728 if (st == 3 && Mod == 2) sign = 1;
3729 if (st == 3 && Mod == 3) sign = 1;
3730 if (st == 3 && Mod == 4) sign = 1;
3731 if (st == 3 && Mod == 5) sign = 1;
3732 if (st == 3 && Mod == 6) sign = -1;
3733 if (st == 3 && Mod == 7) sign = -1;
3741bool BmnSiliconTrackFinder::InvertMatrix(Double_t *
m, Double_t *invOut) {
3742 Double_t inv[16], det;
3745 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];
3746 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];
3747 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];
3748 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];
3749 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];
3750 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];
3751 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];
3752 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];
3753 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];
3754 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];
3755 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];
3756 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];
3757 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];
3758 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];
3759 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];
3760 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];
3761 det =
m[0] * inv[0] +
m[1] * inv[4] +
m[2] * inv[8] +
m[3] * inv[12];
3766 for (
i = 0;
i < 16;
i++)
3767 invOut[
i] = inv[
i] * det;
3774void BmnSiliconTrackFinder::GlobalFit(Double_t **XHits, Double_t **XpHits, Double_t **ZHits, Double_t **SigmaXHits, Double_t **SigmaXpHits, Double_t **Amatr_, Double_t *rhs_) {
3782 Double_t sigmaX_2 = 0.;
3783 Double_t sigmaU_2 = 0.;
3784 Double_t sigmaV_2 = 0.;
3786 for (Int_t iSt = 1; iSt < fNstations; ++iSt) {
3787 for (Int_t imod = 0; imod < fNmodules; ++imod) {
3789 if (iSt == 1 && imod > 3)
continue;
3790 if (iSt == 2 && imod > 1)
continue;
3792 Z = ZHits[iSt][imod];
3795 if (XHits[iSt][imod] > -900.) {
3797 sigmaX_2 = SigmaXHits[iSt][imod] * SigmaXHits[iSt][imod];
3798 X = XHits[iSt][imod];
3801 Amatr_[0][0] += 2 * Z * Z / sigmaX_2;
3802 Amatr_[0][1] += 2 * Z / sigmaX_2;
3803 rhs_[0] += 2 * Z * X / sigmaX_2;
3804 Amatr_[1][0] += 2 * Z / sigmaX_2;
3805 Amatr_[1][1] += 2 / sigmaX_2;
3806 rhs_[1] += 2 * X / sigmaX_2;
3809 if (XpHits[iSt][imod] > -900.) {
3811 if (Angle(iSt, imod) > 0) {
3814 sigmaU_2 = SigmaXpHits[iSt][imod] * SigmaXpHits[iSt][imod];
3815 U = XpHits[iSt][imod];
3818 Amatr_[0][0] += 2 * Z * Z * cos2_5 * cos2_5 / sigmaU_2;
3819 Amatr_[0][1] += 2 * Z * cos2_5 * cos2_5 / sigmaU_2;
3820 Amatr_[0][2] += 2 * Z * Z * cos2_5 * sin2_5 / sigmaU_2;
3821 Amatr_[0][3] += 2 * Z * cos2_5 * sin2_5 / sigmaU_2;
3822 rhs_[0] += 2 * Z * cos2_5 * U / sigmaU_2;
3824 Amatr_[1][0] += 2 * Z * cos2_5 * cos2_5 / sigmaU_2;
3825 Amatr_[1][1] += 2 * cos2_5 * cos2_5 / sigmaU_2;
3826 Amatr_[1][2] += 2 * Z * cos2_5 * sin2_5 / sigmaU_2;
3827 Amatr_[1][3] += 2 * cos2_5 * sin2_5 / sigmaU_2;
3828 rhs_[1] += 2 * cos2_5 * U / sigmaU_2;
3830 Amatr_[2][0] += 2 * Z * Z * cos2_5 * sin2_5 / sigmaU_2;
3831 Amatr_[2][1] += 2 * Z * cos2_5 * sin2_5 / sigmaU_2;
3832 Amatr_[2][2] += 2 * Z * Z * sin2_5 * sin2_5 / sigmaU_2;
3833 Amatr_[2][3] += 2 * Z * sin2_5 * sin2_5 / sigmaU_2;
3834 rhs_[2] += 2 * Z * sin2_5 * U / sigmaU_2;
3836 Amatr_[3][0] += 2 * Z * cos2_5 * sin2_5 / sigmaU_2;
3837 Amatr_[3][1] += 2 * cos2_5 * sin2_5 / sigmaU_2;
3838 Amatr_[3][2] += 2 * Z * sin2_5 * sin2_5 / sigmaU_2;
3839 Amatr_[3][3] += 2 * sin2_5 * sin2_5 / sigmaU_2;
3840 rhs_[3] += 2 * sin2_5 * U / sigmaU_2;
3846 sigmaV_2 = SigmaXpHits[iSt][imod] * SigmaXpHits[iSt][imod];
3848 V = XpHits[iSt][imod];
3851 Amatr_[0][0] += 2 * Z * Z * cos2_5 * cos2_5 / sigmaV_2;
3852 Amatr_[0][1] += 2 * Z * cos2_5 * cos2_5 / sigmaV_2;
3853 Amatr_[0][2] += -2 * Z * Z * cos2_5 * sin2_5 / sigmaV_2;
3854 Amatr_[0][3] += -2 * Z * cos2_5 * sin2_5 / sigmaV_2;
3855 rhs_[0] += 2 * Z * cos2_5 * V / sigmaV_2;
3857 Amatr_[1][0] += 2 * Z * cos2_5 * cos2_5 / sigmaV_2;
3858 Amatr_[1][1] += 2 * cos2_5 * cos2_5 / sigmaV_2;
3859 Amatr_[1][2] += -2 * Z * cos2_5 * sin2_5 / sigmaV_2;
3860 Amatr_[1][3] += -2 * cos2_5 * sin2_5 / sigmaV_2;
3861 rhs_[1] += 2 * cos2_5 * V / sigmaV_2;
3863 Amatr_[2][0] += -2 * Z * Z * cos2_5 * sin2_5 / sigmaV_2;
3864 Amatr_[2][1] += -2 * Z * cos2_5 * sin2_5 / sigmaV_2;
3865 Amatr_[2][2] += 2 * Z * Z * sin2_5 * sin2_5 / sigmaV_2;
3866 Amatr_[2][3] += 2 * Z * sin2_5 * sin2_5 / sigmaV_2;
3867 rhs_[2] += -2 * Z * sin2_5 * V / sigmaV_2;
3869 Amatr_[3][0] += -2 * Z * cos2_5 * sin2_5 / sigmaV_2;
3870 Amatr_[3][1] += -2 * cos2_5 * sin2_5 / sigmaV_2;
3871 Amatr_[3][2] += 2 * Z * sin2_5 * sin2_5 / sigmaV_2;
3872 Amatr_[3][3] += 2 * sin2_5 * sin2_5 / sigmaV_2;
3873 rhs_[3] += -2 * sin2_5 * V / sigmaV_2;
3893Double_t BmnSiliconTrackFinder::FindClusterCenter(Double_t* Ampl, Int_t nElements, Double_t &SumAmpll) {
3897 for (Int_t
i = 0;
i < nElements; ++
i) {
3900 XQ += Double_t(
i + 0.5) * Ampl[
i];
3903 Double_t CoG = XQ / Q;
3914 fInputBranchName =
"BmnSiliconHitClean";
3915 fInputBranchName2 =
"BmnSiliconHitSim";
3919 if (fRunPeriod == 7){
3926 if (fRunPeriod == 8){
3937 hNpoint =
new TH1D(
"Npoint",
"Npoints: X+ X'; Npoints; Events", 7,1,8);fList.Add(hNpoint);
3938 hChiSquareNdf =
new TH1D(
"ChiSquareNdf",
"ChiSquareNdf", 200, 0, Chi2_Globcut);fList.Add(hChiSquareNdf);
3939 hNtracks =
new TH1D(
"Ntracks",
"Ntracks; Ntracks; Events", 10,0,10);fList.Add(hNtracks);
3940 hAxglob =
new TH1D(
"Axglob",
"AngleX; tgX[rad]; Events", 200, -.1,.1);fList.Add(hAxglob);
3941 hBxglob =
new TH1D(
"Bxglob",
"Bx; [cm]; Events", 200, -10,10);fList.Add(hBxglob);
3942 hAyglob =
new TH1D(
"Ayglob",
"AngleY; tgY[rad]; Events", 200, -.1,.1);fList.Add(hAyglob);
3943 hByglob =
new TH1D(
"Byglob",
"By; [cm]; Events", 200, -10,10);fList.Add(hByglob);
3944 hvertexXY =
new TH2D(
"vertexXY",
"vertexXY", 100, -5, 5, 100, -10, 5);fList.Add(hvertexXY);
3945 hvertex_aver_XY =
new TH2D(
"vertex_aver_XY",
"vertex_aver_XY", 100, -5, 5, 100, -10, 5);fList.Add(hvertex_aver_XY);
3946 hprofile_beam_z1 =
new TH2D(
"profile_beam_z1",
"beam_profile_Siz1; [cm];[cm]", 100, -10, 10, 100, -15, 5); fList.Add(hprofile_beam_z1);
3947 hprofile_beam_z2 =
new TH2D(
"profile_beam_z2",
"beam_profile_Siz2; [cm];[cm]", 100, -10, 10, 100, -15, 5);fList.Add(hprofile_beam_z2);
3948 hprofile_beam_z3 =
new TH2D(
"profile_beam_z3",
"beam_profile_Siz3; [cm];[cm]", 100, -15, 15, 100, -20, 5);fList.Add(hprofile_beam_z3);
3949 hAx_first_tr =
new TH1D(
"Ax_first_tr",
"Ax_first_tr; AngleX[mrad]; Events", 200, -.1,.1);fList.Add(hAx_first_tr);
3950 hAx_more_first_tr=
new TH1D(
"Ax_more_first_tr",
"Ax_more_first_tr; AngleX[mrad]; Events", 200, -.1,.1); fList.Add(hAx_more_first_tr);
3951 hYXSi1_run8 =
new TH2D(
"hYXSi1_run8",
"YXSi1_run8", 100, -15,15, 100, -15,15); fList.Add(hYXSi1_run8);
3952 hYXSi2_run8 =
new TH2D(
"hYXSi2_run8",
"YXSi2_run8", 100, -15,15, 100, -15,15); fList.Add(hYXSi2_run8);
3953 hXSi1_run8 =
new TH1D(
"XSi1_run8",
"XSi1; [cm]; Events", 200, -10,10);fList.Add(hXSi1_run8);
3954 hYSi1_run8 =
new TH1D(
"YSi1_run8",
"YSi1; [cm]; Events", 200, -10,10);fList.Add(hYSi1_run8);
3955 hXSi2_run8 =
new TH1D(
"XSi2_run8",
"XSi2; [cm]; Events", 200, -10,10);fList.Add(hXSi2_run8);
3956 hYSi2_run8 =
new TH1D(
"YSi2_run8",
"YSi2; [cm]; Events", 200, -10,10);fList.Add(hYSi2_run8);
3958 hdX_st1_st2 =
new TH1D(
"X_st1_st2",
"X (st1-st2)",100, -0.1, 0.1);fList.Add(hdX_st1_st2);
3959 hdX_st2_st3 =
new TH1D(
"X_st2_st3",
"X (st2-st3)",100, -5., 5.);fList.Add(hdX_st2_st3);
3960 hdX_st1_st3 =
new TH1D(
"X_st1_st3",
"X (st1-st3)",100, -5., 5.);fList.Add(hdX_st1_st3);
3962 hdY_st1_st2 =
new TH1D(
"Y_st1_st2",
"Y (st1-st2)",100, -1, 1);fList.Add(hdY_st1_st2);
3963 hdY_st2_st3 =
new TH1D(
"Y_st2_st3",
"Y (st2-st3)",100, -5., 5.);fList.Add(hdY_st2_st3);
3964 hdY_st1_st3 =
new TH1D(
"Y_st1_st3",
"Y (st1-st3)",100, -5., 5.);fList.Add(hdY_st1_st3);
3967 hX_st_1 =
new TH1D(
"X_st_1",
"X_st_1", 100, -6, 6);fList.Add(hX_st_1);
3968 hX_st_2 =
new TH1D(
"X_st_2",
"X_st_2", 100, -6, 6);fList.Add(hX_st_2);
3969 hX_st_3 =
new TH1D(
"X_st_3",
"X_st_3", 100, -6, 6);fList.Add(hX_st_3);
3971 hdX_st_1 =
new TH1D(
"dX_st_1",
"dX_st_1", 100, -0.02, 0.02);fList.Add(hdX_st_1);
3972 hdX_st_2 =
new TH1D(
"dX_st_2",
"dX_st_2", 100, -0.02, 0.02);fList.Add(hdX_st_2);
3973 hdX_st_3 =
new TH1D(
"dX_st_3",
"dX_st_3", 100, -0.002, 0.002);fList.Add(hdX_st_3);
3975 hdXp_st_1 =
new TH1D(
"dXp_st_1",
"dXp_st_1", 100, -0.02, 0.02);fList.Add(hdXp_st_1);
3976 hdXp_st_2 =
new TH1D(
"dXp_st_2",
"dXp_st_2", 100, -0.02, 0.02);fList.Add(hdXp_st_2);
3977 hdXp_st_3 =
new TH1D(
"dXp_st_3",
"dXp_st_3", 100, -0.02, 0.02);fList.Add(hdXp_st_3);
3979 hdXvsXst_1 =
new TH2D(
"dXvsXst_1",
"dXvsXst_1;XSi[cm];hit-fit[cm]", 100, -6, 6, 100, -0.03, 0.03);fList.Add(hdXvsXst_1);
3980 hdXvsXst_2 =
new TH2D(
"dXvsXst_2",
"dXvsXst_2;XSi[cm];hit-fit[cm]", 100, -6, 6, 100, -0.03, 0.03);fList.Add(hdXvsXst_2);
3981 hdXvsXst_3 =
new TH2D(
"dXvsXst_3",
"dXvsXst_3;XSi[cm];hit-fit[cm]", 100, -6, 6, 100, -0.0005, 0.0005);fList.Add(hdXvsXst_3);
3983 hX13_X2_m0 =
new TH1D(
"X13_X2_m0",
"X2_X13_m0", 200, -0.05, 0.05);fList.Add(hX13_X2_m0);
3984 hX13_X2_m1 =
new TH1D(
"X13_X2_m1",
"X2_X13_m1", 200, -0.05, 0.05);fList.Add(hX13_X2_m1);
3986 hXp13_Xp2_m0 =
new TH1D(
"Xp13_Xp2_m0",
"Xp2_Xp13__m0", 200, -0.4, 0.4);fList.Add(hXp13_Xp2_m0);
3987 hXp13_Xp2_m1 =
new TH1D(
"Xp13_Xp2_m1",
"Xp2_Xp13__m1", 200, -0.4, 0.4);fList.Add(hXp13_Xp2_m1);
3989 hY13_Y2_m0 =
new TH1D(
"Y13_Y2_m0",
"Y2_Y13_m0", 200, -1, 1);fList.Add(hY13_Y2_m0);
3990 hY13_Y2_m1 =
new TH1D(
"Y13_Y2_m1",
"Y2_Y13_m1", 200, -1, 1);fList.Add(hY13_Y2_m1);
3992 hY1m0_Y23 =
new TH1D(
"Y1m0_Y23",
"Y1m0_Y23", 200, -1, 1);fList.Add(hY1m0_Y23);
3993 hY1m1_Y23 =
new TH1D(
"Y1m1_Y23",
"Y1m1_Y23", 200, -1, 1);fList.Add(hY1m1_Y23);
3994 hY1m2_Y23 =
new TH1D(
"Y1m2_Y23",
"Y1m2_Y23", 200, -1, 1);fList.Add(hY1m2_Y23);
3995 hY1m3_Y23 =
new TH1D(
"Y1m3_Y23",
"Y1m3_Y23", 200, -1, 1);fList.Add(hY1m3_Y23);
3997 hY_st_1 =
new TH1D(
"Y_st_1",
"Y_st_1", 100, -6, 6);fList.Add(hY_st_1);
3998 hY_st_2 =
new TH1D(
"Y_st_2",
"Y_st_2", 100, -6, 6);fList.Add(hY_st_2);
3999 hY_st_3 =
new TH1D(
"Y_st_3",
"Y_st_3", 100, -6, 6);fList.Add(hY_st_3);
4001 hdY_st_1 =
new TH1D(
"dY_st_1",
"dY_st_1", 100, -0.6, 0.6);fList.Add(hdY_st_1);
4002 hdY_st_2 =
new TH1D(
"dY_st_2",
"dY_st_2", 100, -0.6, 0.6);fList.Add(hdY_st_2);
4003 hdY_st_3 =
new TH1D(
"dY_st_3",
"dY_st_3", 100, -0.6, 0.6);fList.Add(hdY_st_3);
4005 hdYvsYst_1 =
new TH2D(
"dYvsYst_1",
"dYvsYst_1;YSi[cm];hit-fit[cm]", 100, -6, 6, 100, -0.6, 0.6);fList.Add(hdYvsYst_1);
4006 hdYvsYst_2 =
new TH2D(
"dYvsYst_2",
"dYvsYst_2;YSi[cm];hit-fit[cm]", 100, -6, 6, 100, -0.6, 0.6);fList.Add(hdYvsYst_2);
4007 hdYvsYst_3 =
new TH2D(
"dYvsYst_3",
"dYvsYst_3;YSi[cm];hit-fit[cm]", 100, -6, 6, 100, -0.4, 0.4);fList.Add(hdYvsYst_3);
4009 hdXst1_0_st2_0=
new TH1D(
"dXst1_0_st2_0",
"dXst1_0_st2_0",100, -0.1, 0.1);
4010 hdXst1_0_st2_1=
new TH1D(
"dXst1_0_st2_1",
"dXst1_0_st2_1",100, -0.1, 0.1);
4011 hdXst1_1_st2_0=
new TH1D(
"dXst1_1_st2_0",
"dXst1_1_st2_0",100, -0.1, 0.1);
4012 hdXst1_1_st2_1=
new TH1D(
"dXst1_1_st2_1",
"dXst1_1_st2_1",100, -0.1, 0.1);
4013 hdXst1_2_st2_0=
new TH1D(
"dXst1_2_st2_0",
"dXst1_2_st2_0",100, -0.1, 0.1);
4014 hdXst1_2_st2_1=
new TH1D(
"dXst1_2_st2_1",
"dXst1_2_st2_1",100, -0.1, 0.1);
4015 hdXst1_3_st2_0=
new TH1D(
"dXst1_3_st2_0",
"dXst1_3_st2_0",100, -0.1, 0.1);
4016 hdXst1_3_st2_1=
new TH1D(
"dXst1_3_st2_1",
"dXst1_3_st2_1",100, -0.1, 0.1);
4017 hdXst2_0_st3_1=
new TH1D(
"dXst2_0_st3_1",
"dXst2_0_st3_1",100, -3, 3);
4018 hdXst2_0_st3_2=
new TH1D(
"dXst2_0_st3_2",
"dXst2_0_st3_2",100, -3, 3);
4019 hdXst2_1_st3_1=
new TH1D(
"dXst2_1_st3_1",
"dXst2_1_st3_1",100, -3, 3);
4020 hdXst2_1_st3_2=
new TH1D(
"dXst2_1_st3_2",
"dXst2_1_st3_2",100, -3, 3);
4021 fList.Add(hdXst1_0_st2_0);
4022 fList.Add(hdXst1_1_st2_0);
4023 fList.Add(hdXst1_2_st2_0);
4024 fList.Add(hdXst1_3_st2_0);
4025 fList.Add(hdXst1_0_st2_1);
4026 fList.Add(hdXst1_1_st2_1);
4027 fList.Add(hdXst1_2_st2_1);
4028 fList.Add(hdXst1_3_st2_1);
4029 fList.Add(hdXst2_0_st3_1);
4030 fList.Add(hdXst2_0_st3_2);
4031 fList.Add(hdXst2_1_st3_1);
4032 fList.Add(hdXst2_1_st3_2);
4034 hdXp3_mod1 =
new TH1D(
"hdXp3_mod1",
"hdXp3_mod1",100, -0.02, 0.02);fList.Add(hdXp3_mod1);
4035 hdXp3_mod2 =
new TH1D(
"hdXp3_mod2",
"hdXp3_mod2",100, -0.02, 0.02);fList.Add(hdXp3_mod2);
4036 hdXXp3_mod1 =
new TH1D(
"hdXXp3_mod1",
"hdXXp3_mod1",100, -0.5, 0.5);fList.Add(hdXXp3_mod1);
4037 hdXXp3_mod2 =
new TH1D(
"hdXXp3_mod2",
"hdXXp3_mod2",100, -0.5, 0.5);fList.Add(hdXXp3_mod2);
4038 hXXp12CheckLeftover =
new TH1D(
"hXXp12CheckLeftover",
"hXXp12CheckLeftover",200, -1, 1);fList.Add(hXXp12CheckLeftover);
4039 hXXp12CheckLeftover03 =
new TH1D(
"hXXp12CheckLeftover03",
"hXXp12CheckLeftover03",900, -0.15, 0.15);fList.Add(hXXp12CheckLeftover03);
4041 hdXvsXst1_0 =
new TH2D(
"dXvsXst1_0",
"dXvsXst1_0", 100, -6, 6, 100, -0.05, 0.05);fList.Add(hdXvsXst1_0);
4042 hdXvsXst1_1 =
new TH2D(
"dXvsXst1_1",
"dXvsXst1_1", 100, -6, 6, 100, -0.05, 0.05);fList.Add(hdXvsXst1_1);
4043 hdXvsXst1_2 =
new TH2D(
"dXvsXst1_2",
"dXvsXst1_2", 100, -6, 6, 100, -0.05, 0.05);fList.Add(hdXvsXst1_2);
4044 hdXvsXst1_3 =
new TH2D(
"dXvsXst1_3",
"dXvsXst1_3", 100, -6, 6, 100, -0.05, 0.05);fList.Add(hdXvsXst1_3);
4045 hdXvsXst2_0 =
new TH2D(
"dXvsXst2_0",
"dXvsXst2_0", 100, -6, 6, 100, -0.05, 0.05);fList.Add(hdXvsXst2_0);
4046 hdXvsXst2_1 =
new TH2D(
"dXvsXst2_1",
"dXvsXst2_1", 100, -6, 6, 100, -0.05, 0.05);fList.Add(hdXvsXst2_1);
4047 hdXvsXst3_1 =
new TH2D(
"dXvsXst3_1",
"dXvsXst3_1", 100, -6, 6, 100, -0.001, 0.001);fList.Add(hdXvsXst3_1);
4048 hdXvsXst3_2 =
new TH2D(
"dXvsXst3_2",
"dXvsXst3_2", 100, -6, 6, 100, -0.001, 0.001);fList.Add(hdXvsXst3_2);
4050 hdYvsYst1_mod0 =
new TH2D(
"dYvsYst1_mod0",
"dYvsYst1_mod0", 100, -6, 6, 100, -0.6, 0.6);fList.Add(hdYvsYst1_mod0);
4051 hdYvsYst1_mod1 =
new TH2D(
"dYvsYst1_mod1",
"dYvsYst1_mod1", 100, -6, 6, 100, -0.6, 0.6);fList.Add(hdYvsYst1_mod1);
4052 hdYvsYst1_mod2 =
new TH2D(
"dYvsYst1_mod2",
"dYvsYst1_mod2", 100, -6, 6, 100, -0.6, 0.6);fList.Add(hdYvsYst1_mod2);
4053 hdYvsYst1_mod3 =
new TH2D(
"dYvsYst1_mod3",
"dYvsYst1_mod3", 100, -6, 6, 100, -0.6, 0.6);fList.Add(hdYvsYst1_mod3);
4054 hdYvsYst2_mod0 =
new TH2D(
"dYvsYst2_mod0",
"dYvsYst2_mod0", 100, -6, 6, 100, -0.6, 0.6);fList.Add(hdYvsYst2_mod0);
4055 hdYvsYst2_mod1 =
new TH2D(
"dYvsYst2_mod1",
"dYvsYst2_mod1", 100, -6, 6, 100, -0.6, 0.6);fList.Add(hdYvsYst2_mod1);
4056 hdYvsYst3_mod1 =
new TH2D(
"dYvsYst3_mod1",
"dYvsYst3_mod1", 100, -6, 6, 100, -0.2, 0.2);fList.Add(hdYvsYst3_mod1);
4057 hdYvsYst3_mod2 =
new TH2D(
"dYvsYst3_mod2",
"dYvsYst3_mod2", 100, -6, 6, 100, -0.2, 0.2);fList.Add(hdYvsYst3_mod2);
4059 hdYst1_0_st2_0=
new TH1D(
"dYst1_0_st2_0",
"dYst1_0_st2_0",100, -1, 1);
4060 hdYst1_0_st2_1=
new TH1D(
"dYst1_0_st2_1",
"dYst1_0_st2_1",100, -1, 1);
4061 hdYst1_1_st2_0=
new TH1D(
"dYst1_1_st2_0",
"dYst1_1_st2_0",100, -1, 1);
4062 hdYst1_1_st2_1=
new TH1D(
"dYst1_1_st2_1",
"dYst1_1_st2_1",100, -1, 1);
4063 hdYst1_2_st2_0=
new TH1D(
"dYst1_2_st2_0",
"dYst1_2_st2_0",100, -1, 1);
4064 hdYst1_2_st2_1=
new TH1D(
"dYst1_2_st2_1",
"dYst1_2_st2_1",100, -1, 1);
4065 hdYst1_3_st2_0=
new TH1D(
"dYst1_3_st2_0",
"dYst1_3_st2_0",100, -1, 1);
4066 hdYst1_3_st2_1=
new TH1D(
"dYst1_3_st2_1",
"dYst1_3_st2_1",100, -1, 1);
4067 hdYst2_0_st3_1=
new TH1D(
"dYst2_0_st3_1",
"dYst2_0_st3_1",100, -3, 3);
4068 hdYst2_0_st3_2=
new TH1D(
"dYst2_0_st3_2",
"dYst2_0_st3_2",100, -3, 3);
4069 hdYst2_1_st3_1=
new TH1D(
"dYst2_1_st3_1",
"dYst2_1_st3_1",100, -3, 3);
4070 hdYst2_1_st3_2=
new TH1D(
"dYst2_1_st3_2",
"dYst2_1_st3_2",100, -3, 3);
4071 fList.Add(hdYst1_0_st2_0);
4072 fList.Add(hdYst1_1_st2_0);
4073 fList.Add(hdYst1_2_st2_0);
4074 fList.Add(hdYst1_3_st2_0);
4075 fList.Add(hdYst1_0_st2_1);
4076 fList.Add(hdYst1_1_st2_1);
4077 fList.Add(hdYst1_2_st2_1);
4078 fList.Add(hdYst1_3_st2_1);
4079 fList.Add(hdYst2_0_st3_1);
4080 fList.Add(hdYst2_0_st3_2);
4081 fList.Add(hdYst2_1_st3_1);
4082 fList.Add(hdYst2_1_st3_2);
4084 hY_XSisp1=
new TH2D(
"hY_XSisp1",
"hY_XSisp1", 100, -10, 10, 100, -10, 10);fList.Add(hY_XSisp1);
4085 hY_XSisp2=
new TH2D(
"hY_XSisp2",
"hY_XSisp2", 100, -10, 10, 100, -10, 10);fList.Add(hY_XSisp2);
4086 hY_XSisp3=
new TH2D(
"hY_XSisp3",
"hY_XSisp3", 100, -10, 10, 100, -20, 0);fList.Add(hY_XSisp3);
4088 hY_XSisp1bef=
new TH2D(
"hY_XSisp1bef",
"hY_XSisp1bef", 100, -10, 10, 100, -10, 10);fList.Add(hY_XSisp1bef);
4089 hY_XSisp2bef=
new TH2D(
"hY_XSisp2bef",
"hY_XSisp2bef", 100, -10, 10, 100, -10, 10);fList.Add(hY_XSisp2bef);
4091 N_eff =
new TH1D(
"N_eff",
"N_eff", 6,1,7);fList.Add(N_eff);
4092 D_eff =
new TH1D(
"D_eff",
"D_eff", 6,1,7);fList.Add(D_eff);
4093 E_eff =
new TH1D(
"E_eff",
"Hit Efficiency per Layer on Track", 6,1,7);fList.Add(E_eff);
4096 hdAx_MC_tr =
new TH1D(
"dAx_MC_tr",
"dAx_MC_tr;[rad]",200, -.001,.001);fList.Add(hdAx_MC_tr);
4097 hdAy_MC_tr =
new TH1D(
"dAy_MC_tr",
"dAy_MC_tr;[rad]",200, -.01,.01);fList.Add(hdAy_MC_tr);
4098 hdX_MC_tr =
new TH1D(
"dX_MC_tr",
"dX_MC_tr;[cm]",200, -2,2);fList.Add(hdX_MC_tr);
4099 hdY_MC_tr =
new TH1D(
"dY_MC_tr",
"dY_MC_tr;[cm]",200, -2,2);fList.Add(hdY_MC_tr);
4101 hAxsi_mctrue =
new TH1D(
"Axsi_mctrue",
"Ax_mctrueX; tgX[rad]; Events", 200, -.1,.1);fList.Add(hAxsi_mctrue);
4102 hBxsi_mctrue =
new TH1D(
"Bxsi_mctrue",
"Bx_mctrue; [cm]; Events", 200, -10,10);fList.Add(hBxsi_mctrue);
4103 hAysi_mctrue =
new TH1D(
"Aysi_mctrue",
"Ay_mctrue; tgY[rad]; Events", 200, -.1,.1);fList.Add(hAysi_mctrue);
4104 hBysi_mctrue =
new TH1D(
"Bysi_mctrue",
"By_mctrue; [cm]; Events", 200, -10,10);fList.Add(hBysi_mctrue);
4106 hSi_st1mc =
new TH2D(
"Si_st1mc",
"beam_profileSi_st1mc smeared; [cm];[cm]", 100, -10, 10, 100, -15, 5); fList.Add(hSi_st1mc);
4107 hSi_st2mc =
new TH2D(
"Si_st2mc",
"beam_profileSi_st2mc smeared; [cm];[cm]", 100, -10, 10, 100, -15, 5); fList.Add(hSi_st2mc);
4108 hSi_st3mc =
new TH2D(
"Si_st3mc",
"beam_profileSi_st3mc smeared; [cm];[cm]", 100, -15, 15, 100, -25, 5); fList.Add(hSi_st3mc);
4110 hdAx_MC_tr_comb =
new TH1D(
"dAx_MC_tr_comb",
"dAx_MC_tr_comb;[rad]",200, -.005,.005);fList.Add(hdAx_MC_tr_comb);
4111 hdAy_MC_tr_comb =
new TH1D(
"dAy_MC_tr_comb",
"dAy_MC_tr_comb;[rad]",200, -.005,.005);fList.Add(hdAy_MC_tr_comb);
4112 hdX_MC_tr_comb =
new TH1D(
"dX_MC_tr_comb",
"dX_MC_tr_comb;[cm]",200, -2,2);fList.Add(hdX_MC_tr_comb);
4113 hdY_MC_tr_comb =
new TH1D(
"dY_MC_tr_comb",
"dY_MC_tr_comb;[cm]",200, -2,2);fList.Add(hdY_MC_tr_comb);
4115 hDen_mctrSi =
new TH1D(
"hDen_mctrSi",
"Den_mctrSi", 1, 0, 1);fList.Add(hDen_mctrSi);
4116 hNum_mctrSi =
new TH1D(
"hNum_mctrSi",
"Num_mctrSi", 1, 0, 1);fList.Add(hNum_mctrSi);
4117 hEff_mctrSi =
new TH1D(
"hEff_mctrSi",
"Eff_mctrSi", 1, 0, 1);fList.Add(hEff_mctrSi);
4119 hDen_mcreaction=
new TH1D(
"hDen_mcreaction",
"hDen_mcreaction", 1, 0, 1);fList.Add(hDen_mcreaction);
4120 hNum_mcreaction=
new TH1D(
"hNum_mcreaction",
"hNum_mcreaction", 1, 0, 1);fList.Add(hNum_mcreaction);
4121 hEff_mcreaction=
new TH1D(
"hEff_mcreaction",
"hEff_mcreaction", 1, 0, 1);fList.Add(hEff_mcreaction);
4123 hNtrsi_mc =
new TH1D(
"Ntrsi_mc",
";N of Si-tr mc", 10, 0,10); fList.Add(hNtrsi_mc);
4124 hNtrsi_reco =
new TH1D(
"Ntrsi_reco",
";N of Si-tr reco", 10, 0,10); fList.Add(hNtrsi_reco);
4125 hNtrsi_mc_vs_reco=
new TH2D(
"hNtr_mc_vs_reco",
"hNtr_mc_vs_reco", 10, 0,10, 10, 0,10); fList.Add(hNtrsi_mc_vs_reco);
4129 for (Int_t imod = 0; imod < fNmodules; imod++) {
4130 if (imod > 3)
continue;
4131 h =
new TH1D(Form(
"hoccupancyX1_mod%d",imod), Form(
"hoccupancyX1_mod%d",imod), fNstrips+1, 0.,fNstrips+1);fList.Add(h); hoccupancyX1.push_back(h);
4133 for (Int_t imod = 0; imod < fNmodules; imod++) {
4134 if (imod > 1)
continue;
4135 h =
new TH1D(Form(
"hoccupancyX2_mod%d",imod), Form(
"hoccupancyX2_mod%d",imod), fNstrips+1, 0.,fNstrips+1);fList.Add(h); hoccupancyX2.push_back(h);
4137 for (Int_t imod = 0; imod < fNmodules; imod++) {
4138 h =
new TH1D(Form(
"hoccupancyX3_mod%d",imod), Form(
"hoccupancyX3_mod%d",imod), fNstrips+1, 0.,fNstrips+1);fList.Add(h); hoccupancyX3.push_back(h);
4141 for (Int_t imod = 0; imod < fNmodules; imod++) {
4142 if (imod > 3)
continue;
4143 h =
new TH1D(Form(
"hoccupancyXp1_mod%d",imod), Form(
"hoccupancyXp1_mod%d",imod), fNstrips+1, 0.,fNstrips+1);fList.Add(h); hoccupancyXp1.push_back(h);
4145 for (Int_t imod = 0; imod < fNmodules; imod++) {
4146 if (imod > 1)
continue;
4147 h =
new TH1D(Form(
"hoccupancyXp2_mod%d",imod), Form(
"hoccupancyXp2_mod%d",imod), fNstrips+1, 0.,fNstrips+1);fList.Add(h); hoccupancyXp2.push_back(h);
4149 for (Int_t imod = 0; imod < fNmodules; imod++) {
4150 h =
new TH1D(Form(
"hoccupancyXp3_mod%d",imod), Form(
"hoccupancyXp3_mod%d",imod), fNstrips+1, 0.,fNstrips+1);fList.Add(h); hoccupancyXp3.push_back(h);
4156 AmatrInverted =
new Double_t[16];
4157 rhs =
new Double_t[4];
4158 AmatrArray =
new Double_t[16];
4159 Amatr =
new Double_t*[4];
4160 line =
new Double_t[4];
4161 for (Int_t
i = 0;
i < 4;
i++) {
4162 Amatr[
i] =
new Double_t[4];
4164 Ampl_strX =
new Double_t[fNstrips];
4165 Ampl_strXp =
new Double_t[fNstrips];
4166 Nsp_st =
new Int_t[fNstations];
4167 Nleftoversp =
new Int_t*[fNstations];
4168 NleftoverX =
new Int_t*[fNstations];
4169 NleftoverXp =
new Int_t*[fNstations];
4170 iClustXMod =
new Int_t[fNstations];
4171 iModX =
new Int_t[fNstations];
4172 shiftStXtoGlob =
new Double_t[fNstations];
4173 shiftStYtoGlob =
new Double_t[fNstations];
4174 shiftX =
new Double_t*[fNstations];
4175 shiftY =
new Double_t*[fNstations];
4176 shiftYbelow =
new Double_t*[fNstations];
4177 Zstation =
new Double_t*[fNstations];
4178 Zposition =
new Double_t[fNstations];
4179 NclustX =
new Int_t*[fNstations];
4180 NclustXp =
new Int_t*[fNstations];
4181 NhitsXYmod =
new Int_t*[fNstations];
4182 Xforglfit =
new Double_t*[fNstations];
4183 Xpforglfit =
new Double_t*[fNstations];
4184 Yforglfit =
new Double_t*[fNstations];
4185 SigmXforglfit =
new Double_t*[fNstations];
4186 SigmXpforglfit =
new Double_t*[fNstations];
4187 SigmYforglfit =
new Double_t*[fNstations];
4188 ZstnCforglfit =
new Double_t*[fNstations];
4189 XCoord =
new Double_t**[fNstations];
4190 XpCoord =
new Double_t**[fNstations];
4191 XspCoord =
new Double_t**[fNstations];
4192 XpspCoord =
new Double_t**[fNstations];
4193 YspCoord =
new Double_t**[fNstations];
4194 Sp_pdg =
new Int_t**[fNstations];
4195 SigmaX =
new Double_t**[fNstations];
4196 SigmaXp =
new Double_t**[fNstations];
4197 SigmspX =
new Double_t**[fNstations];
4198 SigmspXp =
new Double_t**[fNstations];
4199 SigmspY =
new Double_t**[fNstations];
4200 DigitsArrayX =
new Double_t**[fNstations];
4201 DigitsArrayXp =
new Double_t**[fNstations];
4202 XClust =
new Double_t**[fNstations];
4203 XpClust =
new Double_t**[fNstations];
4204 XClust_width =
new Double_t**[fNstations];
4205 XpClust_width =
new Double_t**[fNstations];
4206 sumQx =
new Double_t**[fNstations];
4207 sumQxp =
new Double_t**[fNstations];
4208 XspClust =
new Double_t**[fNstations];
4209 XpspClust =
new Double_t**[fNstations];
4210 XspClust_width =
new Double_t**[fNstations];
4211 XpspClust_width=
new Double_t**[fNstations];
4212 sumQxsp =
new Double_t**[fNstations];
4213 sumQxpsp =
new Double_t**[fNstations];
4214 leftoverXsp =
new Double_t**[fNstations];
4215 leftoverXpsp =
new Double_t**[fNstations];
4216 leftoverYsp =
new Double_t**[fNstations];
4217 leftoverXsigsp =
new Double_t**[fNstations];
4218 leftoverXpsigsp=
new Double_t**[fNstations];
4219 leftoverYsigsp =
new Double_t**[fNstations];
4220 leftoverX =
new Double_t**[fNstations];
4221 leftoverXp =
new Double_t**[fNstations];
4222 leftoverXsig =
new Double_t**[fNstations];
4223 leftoverXpsig =
new Double_t**[fNstations];
4225 for (Int_t istat = 0; istat < fNstations; istat++) {
4226 Nleftoversp[istat] =
new Int_t[fNmodules];
4227 NleftoverX[istat] =
new Int_t[fNmodules];
4228 NleftoverXp[istat] =
new Int_t[fNmodules];
4229 leftoverXsp[istat] =
new Double_t*[fNmodules];
4230 leftoverXpsp[istat] =
new Double_t*[fNmodules];
4231 leftoverYsp[istat] =
new Double_t*[fNmodules];
4232 leftoverXsigsp[istat] =
new Double_t*[fNmodules];
4233 leftoverXpsigsp[istat]=
new Double_t*[fNmodules];
4234 leftoverYsigsp[istat] =
new Double_t*[fNmodules];
4235 leftoverX[istat] =
new Double_t*[fNmodules];
4236 leftoverXp[istat] =
new Double_t*[fNmodules];
4237 leftoverXsig[istat] =
new Double_t*[fNmodules];
4238 leftoverXpsig[istat] =
new Double_t*[fNmodules];
4239 shiftX[istat] =
new Double_t[fNmodules];
4240 shiftY[istat] =
new Double_t[fNmodules];
4241 shiftYbelow[istat] =
new Double_t[fNmodules];
4242 Zstation[istat] =
new Double_t[fNmodules];
4243 NclustX[istat] =
new Int_t[fNmodules];
4244 NclustXp[istat] =
new Int_t[fNmodules];
4245 NhitsXYmod[istat] =
new Int_t[fNmodules];
4246 SigmXforglfit[istat] =
new Double_t[fNmodules];
4247 SigmXpforglfit[istat] =
new Double_t[fNmodules];
4248 SigmYforglfit[istat] =
new Double_t[fNmodules];
4249 ZstnCforglfit[istat] =
new Double_t[fNmodules];
4250 Xforglfit[istat] =
new Double_t[fNmodules];
4251 Xpforglfit[istat] =
new Double_t[fNmodules];
4252 Yforglfit[istat] =
new Double_t[fNmodules];
4253 XCoord[istat] =
new Double_t*[fNmodules];
4254 XpCoord[istat] =
new Double_t*[fNmodules];
4255 XspCoord[istat] =
new Double_t*[fNmodules];
4256 XpspCoord[istat] =
new Double_t*[fNmodules];
4257 YspCoord[istat] =
new Double_t*[fNmodules];
4258 Sp_pdg[istat] =
new Int_t*[fNmodules];
4259 SigmaX[istat] =
new Double_t*[fNmodules];
4260 SigmaXp[istat] =
new Double_t*[fNmodules];
4261 SigmspX[istat] =
new Double_t*[fNmodules];
4262 SigmspXp[istat] =
new Double_t*[fNmodules];
4263 SigmspY[istat] =
new Double_t*[fNmodules];
4264 DigitsArrayX[istat] =
new Double_t*[fNmodules];
4265 DigitsArrayXp[istat] =
new Double_t*[fNmodules];
4266 XClust[istat] =
new Double_t*[fNmodules];
4267 XpClust[istat] =
new Double_t*[fNmodules];
4268 XClust_width[istat] =
new Double_t*[fNmodules];
4269 XpClust_width[istat] =
new Double_t*[fNmodules];
4270 sumQx[istat] =
new Double_t*[fNmodules];
4271 sumQxp[istat] =
new Double_t*[fNmodules];
4272 XspClust[istat] =
new Double_t*[fNmodules];
4273 XpspClust[istat] =
new Double_t*[fNmodules];
4274 XspClust_width[istat] =
new Double_t*[fNmodules];
4275 XpspClust_width[istat]=
new Double_t*[fNmodules];
4276 sumQxsp[istat] =
new Double_t*[fNmodules];
4277 sumQxpsp[istat] =
new Double_t*[fNmodules];
4279 for (Int_t im = 0; im < fNmodules; im++) {
4280 leftoverXsp[istat][im] =
new Double_t[kBig];
4281 leftoverXpsp[istat][im] =
new Double_t[kBig];
4282 leftoverYsp[istat][im] =
new Double_t[kBig];
4283 leftoverXsigsp[istat][im] =
new Double_t[kBig];
4284 leftoverXpsigsp[istat][im] =
new Double_t[kBig];
4285 leftoverYsigsp[istat][im] =
new Double_t[kBig];
4286 leftoverX[istat][im] =
new Double_t[kBig];
4287 leftoverXp[istat][im] =
new Double_t[kBig];
4288 leftoverXsig[istat][im] =
new Double_t[kBig];
4289 leftoverXpsig[istat][im] =
new Double_t[kBig];
4290 XCoord[istat][im] =
new Double_t[kBig];
4291 XpCoord[istat][im] =
new Double_t[kBig];
4292 XspCoord[istat][im] =
new Double_t[kBig];
4293 XpspCoord[istat][im] =
new Double_t[kBig];
4294 YspCoord[istat][im] =
new Double_t[kBig];
4295 Sp_pdg[istat][im] =
new Int_t[kBig];
4296 SigmaX[istat][im] =
new Double_t[kBig];
4297 SigmaXp[istat][im] =
new Double_t[kBig];
4298 SigmspX[istat][im] =
new Double_t[kBig];
4299 SigmspXp[istat][im] =
new Double_t[kBig];
4300 SigmspY[istat][im] =
new Double_t[kBig];
4301 DigitsArrayX[istat][im] =
new Double_t[fNstrips];
4302 DigitsArrayXp[istat][im] =
new Double_t[fNstrips];
4303 XClust[istat][im] =
new Double_t[fNstrips];
4304 XpClust[istat][im] =
new Double_t[fNstrips];
4305 XClust_width[istat][im] =
new Double_t[fNstrips];
4306 XpClust_width[istat][im] =
new Double_t[fNstrips];
4307 sumQx[istat][im] =
new Double_t[fNstrips];
4308 sumQxp[istat][im] =
new Double_t[fNstrips];
4309 XspClust[istat][im] =
new Double_t[fNstrips];
4310 XpspClust[istat][im] =
new Double_t[fNstrips];
4311 XspClust_width[istat][im] =
new Double_t[fNstrips];
4312 XpspClust_width[istat][im] =
new Double_t[fNstrips];
4313 sumQxsp[istat][im] =
new Double_t[fNstrips];
4314 sumQxpsp[istat][im] =
new Double_t[fNstrips];
4318 if (fVerbose) cout <<
"BmnSiliconTrackFinder::Init()" << endl;
4319 FairRootManager* ioman = FairRootManager::Instance();
4322 if (fDebug) cout<<
" expData "<<endl;
4323 fBmnSiDigitsArray = (TClonesArray*) ioman->GetObject(
"SILICON");
4324 if (!fBmnSiDigitsArray) {
4325 if (fDebug) cout<<
"fBmnSiDigitsArray::Init(): branch "<<
" not found! Task will be deactivated"<<endl;
4330 if (fDebug) cout<<
" !expData BmnSiliconTrackFinder::Init()" << endl;
4331 fBmnHitsArray = (TClonesArray*)ioman->GetObject(fInputBranchName);
4332 if (fDebug) cout <<
"fInputBranchName = " << fInputBranchName <<
"\n";
4333 if (!fBmnHitsArray) {
4334 cout <<
"BmnSiliconTrackFinder::Init(): branch " <<fInputBranchName<<
" not found! Task will be deactivated" << endl;
4338 fBmnHitsArray2 = (TClonesArray*)ioman->GetObject(fInputBranchName2);
4339 if (fDebug) cout <<
"fInputBranchName2 = " << fInputBranchName2 <<
"\n";
4340 if (!fBmnHitsArray2) {
4341 cout <<
"BmnSiliconTrackFinder::Init(): branch " <<fInputBranchName2<<
" not found! Task will be deactivated" << endl;
4349 half_target_regionX = 3.;
4350 half_target_regionY = 5.5;
4357 Cut_overflow = 1800.;
4359 Cut_AmplStripX = 120.;
4360 Cut_AmplStripXp = 120.;
4364 if (fRunPeriod == 7){
4367 Z0_SRC_target = -645.191;
4369 Shift_toCenterOfMagnetX = 0.54;
4370 Shift_toCenterOfMagnetY = -3.43;
4371 Shift_toCenterOfMagnetAX = 0.;
4372 Shift_toCenterOfMagnetAY = 0.0024;
4374 shiftX[1][0] = -6.4 + 0.1148 + 0.0026; shiftX[2][0] = -5.9 + .1148 + 0.0075; shiftX[3][0] = -12.15 + 0.1148; shiftX[3][4] = -12.15 + 0.1148;
4375 shiftX[1][1] = 5.9 - 0.1148 + 0.0085; shiftX[2][1] = 6.4 - .1148 + 0.009; shiftX[3][1] = -6.15 + 0.1148; shiftX[3][5] = -6.15 + 0.1148;
4376 shiftX[1][2] = -5.9 + 0.1148 + 0.011; shiftX[3][2] = 6.15 - 0.1148; shiftX[3][6] = 6.15 - 0.1148;
4377 shiftX[1][3] = 6.4 - 0.1148 + 0.019; shiftX[3][3] = 12.15 - 0.1148; shiftX[3][7] = 12.15 - 0.1148;
4380 shiftY[1][0] = -0.4 -0.08; shiftY[2][0] = -6.22 -0.02 ; shiftY[3][0] = -0.15; shiftY[3][4] = 0.15;
4381 shiftY[1][1] = 0.1 -0.06; shiftY[2][1] = -6.72 ; shiftY[3][1] = -0.15 -0.25; shiftY[3][5] = 0.15;
4382 shiftY[1][2] = -0.2 -0.05; shiftY[3][2] = -0.15 -0.4; shiftY[3][6] = 0.15;
4383 shiftY[1][3] = 0.4 -0.16; shiftY[3][3] = -0.15; shiftY[3][7] = 0.15;
4385 Zstation[1][0] = 1.44 + 0.5 * 0.003; Zstation[2][0] = 4.71 + 0.5 * 0.003; Zstation[3][0] = 3.51 + 0.5 * 0.003; Zstation[3][4] = 2.25 + 0.5 * 0.003;
4386 Zstation[1][1] = 2.17 + 0.5 * 0.003; Zstation[2][1] = 5.44 + 0.5 * 0.003; Zstation[3][1] = 1.98 + 0.5 * 0.003 - .08; Zstation[3][5] = 3.78 + 0.5 * 0.003;
4387 Zstation[1][2] = 1.44 + 0.5 * 0.003; Zstation[3][2] = 3.29 + 0.5 * 0.003 + .08; Zstation[3][6] = 2.47 + 0.5 * 0.003;
4388 Zstation[1][3] = 2.17 + 0.5 * 0.003; Zstation[3][3] = 1.76 + 0.5 * 0.003; Zstation[3][7] = 4.00 + 0.5 * 0.003;
4391 shiftStXtoGlob[0] = 0.; shiftStYtoGlob[0] = 0.; Zposition[0] = 0.;
4392 shiftStXtoGlob[1] = 0.149 ; shiftStYtoGlob[1] = 0.096 + .191; Zposition[1] = 213.999 -6.87;
4393 shiftStXtoGlob[2] = 0.149 ; shiftStYtoGlob[2] = 0.096 + .151; Zposition[2] = 213.999 -6.87;
4394 shiftStXtoGlob[3] = 0.047 ; shiftStYtoGlob[3] = -6.274 + .174; Zposition[3] = 327.130 + 0.03;
4397 if (fRunPeriod == 8){
4400 Z0_SRC_target = -574.91;
4404 shiftX[1][0] = -6.150 + 0.1148; shiftX[2][0] = -6.130 + 0.1148; shiftX[3][0] = -6.160 + 0.1148; shiftX[4][0] = -6.158 + 0.1148;
4405 shiftX[1][1] = 6.470 - 0.1148; shiftX[2][1] = 6.170 - 0.1148; shiftX[3][1] = 6.122 - 0.1148; shiftX[4][1] = 6.129 - 0.1148;
4409 shiftY[1][0] = 6.203; shiftY[2][0] = 6.218; shiftY[3][0] = 6.206; shiftY[4][0] = 6.209;
4410 shiftY[1][1] = 6.193; shiftY[2][1] = 6.199; shiftY[3][1] = 6.197; shiftY[4][1] = 6.213;
4412 shiftYbelow[1][0] =-0.098; shiftYbelow[2][0] =-0.084; shiftYbelow[3][0] =-0.094; shiftYbelow[4][0] =-0.091;
4413 shiftYbelow[1][1] =-0.104; shiftYbelow[2][1] =-0.097; shiftYbelow[3][1] =-0.102; shiftYbelow[4][0] =-0.089;
4415 Zstation[1][0] =-1.14 -1.3; Zstation[2][0] = 1.145 + 1.3; Zstation[3][0] =-1.095 - 1.3; Zstation[4][0] = 1.140 + 1.3;
4416 Zstation[1][1] =-1.85 -1.3; Zstation[2][1] = 1.870 + 1.3; Zstation[3][1] =-1.850 - 1.3; Zstation[4][1] = 1.875 + 1.3;
4419 shiftStXtoGlob[1] = 0.193; shiftStYtoGlob[1] = 0.153; Zposition[1] = 124.779;
4420 shiftStXtoGlob[2] = 0.193; shiftStYtoGlob[2] = 0.153; Zposition[2] = 126.679;
4421 shiftStXtoGlob[3] = 0.193; shiftStYtoGlob[3] = 0.153; Zposition[3] = 176.679;
4422 shiftStXtoGlob[4] = 0.193; shiftStYtoGlob[4] = 0.153; Zposition[4] = 178.679;
4424 Shift_toCenterOfMagnetX = 0.;
4425 Shift_toCenterOfMagnetY = 0.;
4426 Shift_toCenterOfMagnetAX = 0.;
4427 Shift_toCenterOfMagnetAY = 0.;
4438 Cut_overflow = 1800.;
4440 Cut_AmplStripX = 0.;
4441 Cut_AmplStripXp = 0.;
4445 Shift_toCenterOfMagnetX = 0.;
4446 Shift_toCenterOfMagnetY = 0.;
4447 Shift_toCenterOfMagnetAX = 0.;
4448 Shift_toCenterOfMagnetAY = 0.;
4451 shiftX[1][0] = -6.4 + 0.1148; shiftX[2][0] = -5.9 + .1148; shiftX[3][0] = -12.15 + 0.1148; shiftX[3][4] = -12.15 + 0.1148;
4452 shiftX[1][1] = 5.9 - 0.1148; shiftX[2][1] = 6.4 - .1148; shiftX[3][1] = -6.15 + 0.1148; shiftX[3][5] = -6.15 + 0.1148;
4453 shiftX[1][2] = -5.9 + 0.1148; shiftX[3][2] = 6.15 - 0.1148; shiftX[3][6] = 6.15 - 0.1148;
4454 shiftX[1][3] = 6.4 - 0.1148; shiftX[3][3] = 12.15 - 0.1148; shiftX[3][7] = 12.15 - 0.1148;
4457 shiftY[1][0] = -0.4 ; shiftY[2][0] = -6.22; shiftY[3][0] = -0.15; shiftY[3][4] = 0.15;
4458 shiftY[1][1] = 0.1 ; shiftY[2][1] = -6.72; shiftY[3][1] = -0.15; shiftY[3][5] = 0.15;
4459 shiftY[1][2] = -0.2 ; shiftY[3][2] = -0.15; shiftY[3][6] = 0.15;
4460 shiftY[1][3] = 0.4 ; shiftY[3][3] = -0.15; shiftY[3][7] = 0.15;
4462 shiftStXtoGlob[0] = 0.; shiftStYtoGlob[0] = 0.;
4463 shiftStXtoGlob[1] = 0.; shiftStYtoGlob[1] = 0.;
4464 shiftStXtoGlob[2] = 0.; shiftStYtoGlob[2] = 0.;
4465 shiftStXtoGlob[3] = 0.; shiftStYtoGlob[3] = 0.;
4468 Zposition[1] = 213.999 -6.9 + 0.03;
4469 Zposition[2] = 213.999 -6.9 + 0.03;
4470 Zposition[3] = 327.130 + 0.03;
4472 Zstation[1][0] = 1.44 + 0.5 * 0.003; Zstation[2][0] = 4.71 + 0.5 * 0.003; Zstation[3][0] = 3.51 + 0.5 * 0.003; Zstation[3][4] = 2.25 + 0.5 * 0.003;
4473 Zstation[1][1] = 2.17 + 0.5 * 0.003; Zstation[2][1] = 5.44 + 0.5 * 0.003; Zstation[3][1] = 1.98 + 0.5 * 0.003 - .08; Zstation[3][5] = 3.78 + 0.5 * 0.003;
4474 Zstation[1][2] = 1.44 + 0.5 * 0.003; Zstation[3][2] = 3.29 + 0.5 * 0.003 + .08; Zstation[3][6] = 2.47 + 0.5 * 0.003;
4475 Zstation[1][3] = 2.17 + 0.5 * 0.003; Zstation[3][3] = 1.76 + 0.5 * 0.003; Zstation[3][7] = 4.00 + 0.5 * 0.003;
4479 for (Int_t is = 1; is < fNstations; is++) {
4480 for (
int im = 0; im < fNmodules; im++) {
4481 if ( fRunPeriod == 7 && is == 1 && im > 3)
continue;
4482 if ( fRunPeriod == 7 && is == 2 && im > 1)
continue;
4483 Zstation[is][im] += Zposition[is];
4484 Zstation[is][im] += Z0_SRC_target;
4487 ZstnCforglfit[is][im] = Zstation[is][im] - Zcentr;
4492 for (Int_t ist = 1; ist < fNstations; ist++) {
4494 for (Int_t imod = 0; imod < fNmodules; imod++) {
4495 if (ist == 1 && imod > 3 )
continue;
4496 if (ist == 2 && imod > 1 )
continue;
4503 fBmnSiliconHitsArray =
new TClonesArray(fOutputHitsBranchName);
4504 ioman->Register(
"BmnSiliconHits",
"SI", fBmnSiliconHitsArray, kTRUE);
4506 fBmnSiliconHitsXArray =
new TClonesArray(fOutputHitsBranchName);
4507 ioman->Register(
"BmnSiliconHitsX",
"SI", fBmnSiliconHitsXArray, kTRUE);
4509 fBmnSiliconHitsXpArray =
new TClonesArray(fOutputHitsBranchName);
4510 ioman->Register(
"BmnSiliconHitsXp",
"SI", fBmnSiliconHitsXpArray, kTRUE);
4512 fBmnSiliconHitsXYArray =
new TClonesArray(fOutputHitsBranchName);
4513 ioman->Register(
"BmnSiliconHitsXY",
"SI", fBmnSiliconHitsXYArray, kTRUE);
4516 fSiTracks =
new TClonesArray(
"BmnSiliconTrack");
4517 ioman->Register(
"BmnSiliconTrack",
"SI", fSiTracks, kTRUE);
4525void BmnSiliconTrackFinder::PrepareArraysToProcessEvent() {
4527 fBmnSiliconHitsArray->Delete();
4528 fBmnSiliconHitsXArray->Delete();
4529 fBmnSiliconHitsXpArray->Delete();
4530 fBmnSiliconHitsXYArray->Delete();
4535 for (Int_t istat = 0; istat < fNstations; istat++) {
4538 for (Int_t im = 0; im < fNmodules; im++) {
4540 Nleftoversp[istat][im] = 0;
4541 NleftoverX[istat][im] = 0;
4542 NleftoverXp[istat][im] = 0;
4543 NclustX[istat][im] = 0;
4544 NclustXp[istat][im] = 0;
4545 NhitsXYmod[istat][im] = 0;
4546 Xforglfit[istat][im] = -999.;
4547 Xpforglfit[istat][im] = -999.;
4548 Yforglfit[istat][im] = -999.;
4549 SigmXforglfit[istat][im] = -999.;
4550 SigmXpforglfit[istat][im] = -999.;
4551 SigmYforglfit[istat][im] = -999.;
4553 for (Int_t ii = 0; ii < fNstrips; ii++) {
4555 Ampl_strXp[ii] = 0.;
4556 DigitsArrayX[istat][im][ii] = 0.;
4557 DigitsArrayXp[istat][im][ii] = 0.;
4560 for (Int_t ii = 0; ii < kBig; ii++) {
4561 leftoverXsp[istat][im][ii] = -999.;
4562 leftoverXpsp[istat][im][ii] = -999.;
4563 leftoverYsp[istat][im][ii] = -999.;
4564 leftoverXsigsp[istat][im][ii] = -999.;
4565 leftoverXpsigsp[istat][im][ii] = -999.;
4566 leftoverYsigsp[istat][im][ii] = -999.;
4567 leftoverX[istat][im][ii] = -999.;
4568 leftoverXp[istat][im][ii] = -999.;
4569 leftoverXsig[istat][im][ii] = -999.;
4570 leftoverXpsig[istat][im][ii] = -999.;
4571 XCoord[istat][im][ii] = -999.;
4572 XpCoord[istat][im][ii] = -999.;
4573 XspCoord[istat][im][ii] = -999.;
4574 XpspCoord[istat][im][ii] = -999.;
4575 YspCoord[istat][im][ii] = -999.;
4576 Sp_pdg[istat][im][ii] = -1;
4577 SigmaX[istat][im][ii] = -999.;
4578 SigmaXp[istat][im][ii] = -999.;
4579 SigmspX[istat][im][ii] = -999.;
4580 SigmspXp[istat][im][ii] = -999.;
4581 SigmspY[istat][im][ii] = -999.;
4582 XClust[istat][im][ii] = 0.;
4583 XpClust[istat][im][ii] = 0.;
4584 XClust_width[istat][im][ii] = 0.;
4585 XpClust_width[istat][im][ii] = 0.;
4586 sumQx[istat][im][ii] = 0.;
4587 sumQxp[istat][im][ii] = 0.;
4588 XspClust[istat][im][ii] = 0.;
4589 XpspClust[istat][im][ii] = 0.;
4590 XspClust_width[istat][im][ii] = 0.;
4591 XpspClust_width[istat][im][ii] = 0.;
4592 sumQxsp[istat][im][ii] = 0.;
4593 sumQxpsp[istat][im][ii] = 0.;
4604 fOutputHitsBranchName =
"BmnSiliconHit";
4610 fRunNumber = runNumber;
4612 fRunPeriod = RunPeriod;
4613 fOutputHitsBranchName =
"BmnSiliconHit";
4621 for (Int_t istat = 0; istat < fNstations; istat++) {
4623 for (Int_t im = 0; im < fNmodules; im++) {
4624 delete[] XCoord[istat][im];
4625 delete[] XpCoord[istat][im];
4626 delete[] XspCoord[istat][im];
4627 delete[] XpspCoord[istat][im];
4628 delete[] YspCoord[istat][im];
4629 delete[] Sp_pdg[istat][im];
4630 delete[] SigmaX[istat][im];
4631 delete[] SigmaXp[istat][im];
4632 delete[] SigmspX[istat][im];
4633 delete[] SigmspXp[istat][im];
4634 delete[] SigmspY[istat][im];
4635 delete[] leftoverXsp[istat][im];
4636 delete[] leftoverXpsp[istat][im];
4637 delete[] leftoverYsp[istat][im];
4638 delete[] leftoverXsigsp[istat][im];
4639 delete[] leftoverXpsigsp[istat][im];
4640 delete[] leftoverYsigsp[istat][im];
4641 delete[] leftoverX[istat][im];
4642 delete[] leftoverXp[istat][im];
4643 delete[] leftoverXsig[istat][im];
4644 delete[] leftoverXpsig[istat][im];
4645 delete[] DigitsArrayX[istat][im];
4646 delete[] DigitsArrayXp[istat][im];
4647 delete[] XClust[istat][im];
4648 delete[] XpClust[istat][im];
4649 delete[] XClust_width[istat][im];
4650 delete[] XpClust_width[istat][im];
4651 delete[] sumQx[istat][im];
4652 delete[] sumQxp[istat][im];
4653 delete[] XspClust[istat][im];
4654 delete[] XpspClust[istat][im];
4655 delete[] XspClust_width[istat][im];
4656 delete[] XpspClust_width[istat][im];
4657 delete[] sumQxsp[istat][im];
4658 delete[] sumQxpsp[istat][im];
4661 delete[] NclustX[istat];
4662 delete[] NclustXp[istat];
4663 delete[] NhitsXYmod[istat];
4664 delete[] XCoord[istat];
4665 delete[] XpCoord[istat];
4666 delete[] XspCoord[istat];
4667 delete[] XpspCoord[istat];
4668 delete[] YspCoord[istat];
4669 delete[] Sp_pdg[istat];
4670 delete[] SigmaX[istat];
4671 delete[] SigmaXp[istat];
4672 delete[] SigmspX[istat];
4673 delete[] SigmspXp[istat];
4674 delete[] SigmspY[istat];
4675 delete[] Xforglfit[istat];
4676 delete[] Xpforglfit[istat];
4677 delete[] Yforglfit[istat];
4678 delete[] SigmXforglfit[istat];
4679 delete[] SigmXpforglfit[istat];
4680 delete[] SigmYforglfit[istat];
4681 delete[] ZstnCforglfit[istat];
4682 delete[] shiftX[istat];
4683 delete[] shiftY[istat];
4684 delete[] shiftYbelow[istat];
4685 delete[] Zstation[istat];
4686 delete[] Nleftoversp[istat];
4687 delete[] NleftoverX[istat];
4688 delete[] NleftoverXp[istat];
4689 delete[] leftoverXsp[istat];
4690 delete[] leftoverXpsp[istat];
4691 delete[] leftoverYsp[istat];
4692 delete[] leftoverXsigsp[istat];
4693 delete[] leftoverXpsigsp[istat];
4694 delete[] leftoverYsigsp[istat];
4695 delete[] leftoverX[istat];
4696 delete[] leftoverXp[istat];
4697 delete[] leftoverXsig[istat];
4698 delete[] leftoverXpsig[istat];
4700 delete[] DigitsArrayX[istat];
4701 delete[] DigitsArrayXp[istat];
4702 delete[] XClust[istat];
4703 delete[] XpClust[istat];
4704 delete[] XClust_width[istat];
4705 delete[] XpClust_width[istat];
4706 delete[] sumQx[istat];
4707 delete[] sumQxp[istat];
4708 delete[] XspClust[istat];
4709 delete[] XpspClust[istat];
4710 delete[] XspClust_width[istat];
4711 delete[] XpspClust_width[istat];
4712 delete[] sumQxsp[istat];
4713 delete[] sumQxpsp[istat];
4718 delete[] NhitsXYmod;
4731 delete[] Xpforglfit;
4733 delete[] SigmXforglfit;
4734 delete[] SigmXpforglfit;
4735 delete[] SigmYforglfit;
4736 delete[] ZstnCforglfit;
4737 delete[] shiftStXtoGlob;
4738 delete[] shiftStYtoGlob;
4741 delete[] shiftYbelow;
4744 delete[] Nleftoversp;
4745 delete[] NleftoverX;
4746 delete[] NleftoverXp;
4747 delete[] leftoverXsp;
4748 delete[] leftoverXpsp;
4749 delete[] leftoverYsp;
4750 delete[] leftoverXsigsp;
4751 delete[] leftoverXpsigsp;
4752 delete[] leftoverYsigsp;
4754 delete[] leftoverXp;
4755 delete[] leftoverXsig;
4756 delete[] leftoverXpsig;
4758 delete[] iClustXMod;
4761 delete[] Ampl_strXp;
4762 delete[] DigitsArrayX;
4763 delete[] DigitsArrayXp;
4766 delete[] XClust_width;
4767 delete[] XpClust_width;
4772 delete[] XspClust_width;
4773 delete[] XpspClust_width;
4777 for (Int_t
i = 0;
i < 4;
i++) {
4780 delete[] AmatrInverted;
4782 delete[] AmatrArray;
4796 E_eff->Divide(N_eff,D_eff,1,1);
4797 hEff_mctrSi->Divide(hNum_mctrSi,hDen_mctrSi,1,1);
4798 hEff_mcreaction->Divide(hNum_mcreaction,hDen_mcreaction,1,1);
4800 printf(
"BmnSiliconTrackFinder: write hists to file... ");
4801 fOutputFileName = Form(
"hSiliconTracks_run%d.root", fRunNumber);
4802 cout<< fOutputFileName <<endl;
4803 TFile file(fOutputFileName,
"RECREATE");
4808 cout <<
"Work time of the SI track finder: " << workTime <<
" s" << endl;
friend F32vec4 sqrt(const F32vec4 &a)
friend F32vec4 fabs(const F32vec4 &a)
void SetModule(Int_t mod)
void SetStation(Short_t st)
void SetStripTotalSignalInUpperLayer(Double_t sig)
void SetClusterSizeInUpperLayer(Int_t csize)
void SetClusterSizeInLowerLayer(Int_t csize)
void SetStripPositionInLowerLayer(Double_t spos)
void SetStripTotalSignalInLowerLayer(Double_t sig)
void SetStripPositionInUpperLayer(Double_t spos)
virtual void Exec(Option_t *opt)
virtual ~BmnSiliconTrackFinder()
virtual InitStatus Init()
static bool compareSegments(const tracksX &a, const tracksX &b)
Double_t GetStripSignal()
void SetChi2(Double_t chi2)
void SetParamLast(FairTrackParam &par)
void SetParamFirst(FairTrackParam &par)