7 #include <GaudiKernel/IMessageSvc.h>
38 if (
stationName.substr(0,3)==
"BML" && stationEta==7) ilay+=2;
42 double R=std::sqrt(gPosX*gPosX+gPosY*gPosY);
43 const double Phi=std::atan2(gPosY,gPosX);
62 std::array<std::reference_wrapper<double>, 3>& result_aw,
63 std::array<std::reference_wrapper<double>, 3>& result_bw,
69 if(rpc_x.at(6).size()+rpc_x.at(7).size() >0) layer_end = 7;
72 std::bitset<8> result_pat{};
73 double result_dMM{9999}, result_dMO{9999};
75 std::array<size_t, 8> result_index{};
78 for (
int l_start=0; l_start<layer_end; l_start++){
80 for (
size_t i_start = 0; i_start < rpc_x.at(l_start).
size(); ++i_start){
84 std::bitset<8>
pat (1<<l_start);
85 std::array<size_t, 8>
index {};
86 index[l_start] = i_start;
90 int current_l = l_start;
91 double current_x = rpc_x.at(l_start).at(i_start);
95 for (
int l_test=l_start+1; l_test<=layer_end; l_test++){
96 double min_delta {999};
97 const std::vector<double>& test_layer_hits {rpc_x.at(l_test)};
99 for (
size_t i_test = 0; i_test < test_layer_hits.size(); ++i_test){
102 if (deltaOK(current_l,l_test,current_x, test_layer_hits.at(i_test),
false,delta)){
104 if (delta < min_delta) {
106 index[l_test] = i_test;
110 if (min_delta < 998){
112 current_x = test_layer_hits.at(
index[l_test]);
115 dMO = (l_start<4 and l_test>=4) ?
std::min(dMO, min_delta) : dMO;
116 dMM = (l_start<2 and l_test>=2 and l_test<4) or (l_start>=4 and l_start<5 and l_test>=6) ?
std::min(dMM, min_delta) : dMM;
121 if (nHits>nHits_pat) {
127 }
else if (nHits==nHits_pat) {
129 if (dMM<result_dMM or (dMM==result_dMM and dMO<result_dMO)){
140 abcal(result_pat, result_index, result_aw, result_bw, rpcLayerHits);
142 std::ostringstream ossR, ossZ;
144 for (
int i=0;
i<8; ++
i){
145 if(result_pat.test(
i)){
156 std::ostringstream oss;
157 std::copy(result_index.begin(), result_index.end(), std::ostream_iterator<int>(oss,
" "));
158 ATH_MSG_VERBOSE(
"patfinder: BEST pat= " << result_pat <<
" nHit: " << nHits_pat <<
" Idx: " << oss.str()
159 <<
" dMM= "<<result_dMM <<
" dMO= "<<result_dMO <<
" R_hits: " << ossR.str() <<
" Z_hits: " << ossZ.str()
160 <<
" Slopes: " << result_aw[0] <<
"," << result_aw[1] <<
"," << result_aw[2] <<
" Offsets: " << result_bw[0] <<
"," << result_bw[1] <<
"," << result_bw[2]);
172 const int N_layers=8;
174 const std::vector<std::vector<double>>& rpc_phi {rpcLayerHits.
hits_in_layer_phi};
177 if (rpc_phi.at(6).size()+rpc_phi.at(7).size()>0) l_start_max=5;
182 double result_dMM{9999}, result_dMO{9999};
186 for (
int l_start=0; l_start<=l_start_max; l_start++){
188 for (
const double& phi_start : rpc_phi.at(l_start)){
194 int current_l = l_start;
195 double current_phi = phi_start;
199 for (
int l_test=l_start+1; l_test<N_layers; l_test++){
200 double min_delta {999};
201 double layer_phi {0.};
203 for (
const double& phi_test : rpc_phi.at(l_test)){
206 if (deltaOK(current_l,l_test,current_phi, phi_test,
true,delta)){
208 if (delta < min_delta) {
210 layer_phi = phi_test;
214 if (min_delta < 998){
216 current_phi = layer_phi;
218 dMO = (l_start<4 and l_test>=4) ?
std::min(dMO, min_delta) : dMO;
219 dMM = (l_start<2 and l_test>=2 and l_test<4) or (l_start>=4 and l_start<5 and l_test>=6) ?
std::min(dMM, min_delta) : dMM;
224 if (nHits>nHits_pat and current_l > 1) {
228 phi_middle=phi_start;
229 phi_outer=current_phi;
230 }
else if (nHits==nHits_pat and current_l > 1) {
232 if (dMM<result_dMM or (dMM==result_dMM and dMO<result_dMO)){
235 phi_middle=phi_start;
236 phi_outer=current_phi;
243 ATH_MSG_DEBUG(
"patfinder: BEST phi path dMM= "<<result_dMM <<
" dMO= "<<result_dMO
244 <<
" phi_middle= "<<phi_middle <<
" phi_outer= "<<phi_outer);
271 const double delta_gasgap_eta = 0.01;
272 const double delta_lowpt_eta = 0.05;
273 const double delta_highpt_eta = 0.1;
274 const double delta_feet_eta = 0.05;
276 const double delta_gasgap_phi = 0.01;
277 const double delta_lowpt_phi = 0.03;
278 const double delta_highpt_phi = 0.04;
279 const double delta_feet_phi = 0.03;
283 else delta=std::abs(
x2-
x1);
294 delta_max=delta_gasgap_phi;
295 }
else if (
l1<2&&
l2<4) {
296 delta_max=delta_lowpt_phi;
297 }
else if (l1<4&&l2>=4) {
298 delta_max=delta_highpt_phi;
299 }
else if (l1<6&&l1>=4&&
l2>=6) {
300 delta_max=delta_feet_phi;
304 delta_max=delta_gasgap_eta;
305 }
else if (l1<2&&l2>=2&&
l2<4) {
306 delta_max=delta_lowpt_eta;
307 }
else if (l1<4&&l2>=4) {
308 delta_max=delta_highpt_eta;
309 }
else if (l1<6&&l1>=4&&
l2>=6) {
310 delta_max=delta_feet_eta;
317 if (delta<delta_max) OK=
true;
327 double DeltaPhi, temp_phi;
331 DeltaPhi= 999; temp_phi=9999;
332 for(
int inum=0;inum < 8;inum++){
333 temp_phi = std::abs((inum *
M_PI/4.0 )- calibPhi);
334 if(temp_phi < DeltaPhi) DeltaPhi = temp_phi;
336 }
else if(std::string::npos !=
stationName.rfind(
'S') ||
339 DeltaPhi= 999; temp_phi=9999;
341 for(
int inum=0;inum < 8;inum++){
342 temp_phi = std::abs(inum *(
M_PI/4.0 )+(
M_PI/8.0) - calibPhi);
343 if(temp_phi < DeltaPhi) DeltaPhi = temp_phi;
347 double calibR = R *
std::cos(DeltaPhi);
356 const std::array<size_t, 8>&
index,
357 std::array<std::reference_wrapper<double>, 3>& aw,
358 std::array<std::reference_wrapper<double>, 3>& bw,
361 const std::vector<std::vector<double>>& rpc_R {rpcLayerHits.
hits_in_layer_R};
362 const std::vector<std::vector<double>>& rpc_Z {rpcLayerHits.
hits_in_layer_Z};
365 auto getCompanion = [](
const int&
l) ->
int {
366 return (
l%2==0) ?
l+1 :
l-1;
369 auto getAvgRZ= [&rpc_R, &rpc_Z, &result_pat, &
index]
370 (
double& R,
double& Z,
const int&
l,
const int& companion) ->
void {
372 if (result_pat.test(companion)){
373 R = (rpc_R.at(
l).at(
index[
l]) + rpc_R.at(companion).at(
index[companion])) / 2.0;
374 Z = (rpc_Z.at(
l).at(
index[
l]) + rpc_Z.at(companion).at(
index[companion])) / 2.0;
383 double R1{0}, R2{0}, Z1{0}, Z2{0};
384 for(
int i=0;
i<4;
i++){
385 if(!result_pat.test(
i))
continue;
389 const int comp1{getCompanion(
l1)}, comp2{getCompanion(
l2)};
390 if (
l1 != 4 and
l1 !=
l2 and
l2 != comp1){
392 getAvgRZ(R1, Z1,
l1, comp1);
393 getAvgRZ(R2, Z2,
l2, comp2);
396 if(((result_pat & std::bitset<8>(
"00001111")).
count() > 3) and std::abs(Z2-Z1) >
ZERO_LIMIT){
398 double theta_m {std::atan2(R1, Z1)};
399 double theta_t {std::atan2(R2-R1, Z2-Z1)};
400 double theta_f {(theta_m+theta_t)/2.};
403 bw[0].get() = R1 - Z1*aw[0].get();
412 aw[1].get() = (R2-R1)/(Z2-Z1);
413 bw[1].get() = R2 - Z2*aw[1].get();
421 getAvgRZ(R1, Z1,
l1, comp1);
424 aw[1].get() = aw[0].get();
431 double R3{0}, R4{0}, Z3{0}, Z4{0};
432 for(
int i=2;
i<8;
i++){
433 if(!result_pat.test(
i))
continue;
437 const int comp3{getCompanion(l3)}, comp4{getCompanion(l4)};
438 if ( l3 != 8 and l3 != l4 and l4 != comp3){
440 getAvgRZ(R3, Z3, l3, comp3);
441 getAvgRZ(R4, Z4, l4, comp4);
444 aw[2].get() = (R4-R3)/(Z4-Z3);
445 bw[2].get() = R4 - Z4*aw[2].get();
453 getAvgRZ(R3, Z3, l3, comp3);
459 getAvgRZ(R3, Z3, l3, comp3);
461 aw[1].get() = aw[0].get();
464 getAvgRZ(R2, Z2,
l2, comp2);