7 static double const TWO_PI =
M_PI * 2;
8 while ( phi < -
M_PI ) phi += TWO_PI;
9 while ( phi >
M_PI ) phi -= TWO_PI;
140 double *
obj,
double *rhs,
double *lbound,
double *ubound,
142 double *diag,
double *odiag,
146 double *dxs,
double *dxsn,
double *
up,
150 double *ddspr,
double *ddsprn,
double *dsup,
double *ddsup,
155 double *ddv,
double *ddvn,
double *prinf,
157 double *upinf,
double *duinf,
double *
scale,
167 int *ecolpnt,
int *
count,
int *vcstat,
int *pivots,
168 int *invprm,
int *snhead,
int *nodtyp,
int *inta1,
178 int *iter,
int *corect,
int *fixn,
int *dropn,
int *fnzmax,
179 int *fnzmin,
double *addobj,
double *bigbou,
189 #pragma clang diagnostic push
190 #pragma clang diagnostic ignored "-Wunused-const-variable"
194 static const int minimize = 1;
195 static const int maximize = -1;
196 static const char less_equal =
'L';
197 static const char equal =
'E';
198 static const char greater_equal =
'G';
199 static const char range =
'R';
200 #pragma clang diagnostic pop
203 bool _variable_named;
205 std::vector<double> _rhs;
206 std::vector<char> _constraint_sense;
207 std::vector<double> _range_value;
208 std::vector<char *> _row_name;
210 std::vector<double> _objective;
211 std::vector<double> _lower_bound;
212 std::vector<double> _upper_bound;
213 std::vector<char *> _column_name;
215 std::vector<int> _coefficient_row;
216 std::vector<int> _coefficient_column;
217 std::vector<double> _coefficient_value;
222 _constraint_sense.clear();
223 for (std::vector<char *>::const_iterator iterator =
225 iterator != _row_name.end(); iterator++) {
230 void clear_column(
void)
233 _lower_bound.clear();
234 _upper_bound.clear();
235 for (std::vector<char *>::const_iterator iterator =
236 _column_name.begin();
237 iterator != _column_name.end(); iterator++) {
240 _column_name.clear();
242 void clear_coefficient(
void)
244 _coefficient_row.clear();
245 _coefficient_column.clear();
246 _coefficient_value.clear();
248 void to_csr(std::vector<int> &column_pointer,
249 std::vector<int> &row_index,
250 std::vector<double> &nonzero_value,
251 const int index_start = 1)
256 std::vector<std::vector<std::pair<int, double> > >
257 column_index(_objective.size(),
258 std::vector<std::pair<int, double> >());
260 std::vector<int>::const_iterator iterator_row =
261 _coefficient_row.begin();
262 std::vector<int>::const_iterator iterator_column =
263 _coefficient_column.begin();
264 std::vector<double>::const_iterator iterator_value =
265 _coefficient_value.begin();
267 for (; iterator_value != _coefficient_value.end();
268 iterator_row++, iterator_column++, iterator_value++) {
269 column_index[*iterator_column].push_back(
270 std::pair<int, double>(
272 *iterator_row + index_start, *iterator_value));
275 for (std::vector<std::vector<std::pair<int, double> > >::
276 iterator iterator_outer = column_index.begin();
277 iterator_outer != column_index.end(); iterator_outer++) {
279 column_pointer.push_back(row_index.size() + index_start);
280 std::sort(iterator_outer->begin(), iterator_outer->end());
281 for (std::vector<std::pair<int, double> >::const_iterator
282 iterator_inner = iterator_outer->begin();
283 iterator_inner != iterator_outer->end();
285 row_index.push_back(iterator_inner->first);
286 nonzero_value.push_back(iterator_inner->second);
290 column_pointer.push_back(row_index.size() + index_start);
293 problem_t(
const bool variable_named)
294 : _optimized(
false), _variable_named(variable_named)
307 virtual int populate(
void) = 0;
308 size_t nrow(
void)
const
312 size_t ncolumn(
void)
const
314 return _objective.size();
316 void push_back_row(
const char constraint_sense,
320 _constraint_sense.push_back(constraint_sense);
322 if (_variable_named) {
326 snprintf(name, name_length,
"c%llu",
327 static_cast<unsigned long long>(_rhs.size()));
328 _row_name.push_back(name);
331 void push_back_row(
const std::string &constraint_sense,
336 if (constraint_sense ==
"<=") {
338 push_back_row(rhs, cplex_sense);
340 else if (constraint_sense ==
"==") {
342 push_back_row(rhs, cplex_sense);
344 else if (constraint_sense ==
">=") {
346 push_back_row(rhs, cplex_sense);
349 fprintf(stderr,
"%s:%d: illegal sense (`%s')\n",
350 __FILE__, __LINE__, constraint_sense.c_str());
353 void push_back_column(
const double objective,
354 const double lower_bound,
355 const double upper_bound)
357 _objective.push_back(objective);
358 _lower_bound.push_back(lower_bound);
359 _upper_bound.push_back(upper_bound);
361 if (_variable_named) {
362 static const size_t name_length = 24;
365 snprintf(name, name_length,
"x%llu",
366 static_cast<unsigned long long>(
368 _column_name.push_back(name);
371 void push_back_coefficient(
372 const int row,
const int column,
const double value)
374 _coefficient_row.push_back(row);
375 _coefficient_column.push_back(column);
376 _coefficient_value.push_back(value);
378 virtual int optimize(
void) = 0;
379 int optimize_primal_simplex(
void)
383 int optimize_dual_simplex(
void)
387 int optimize_barrier(
void)
391 int optimize_network(
void)
395 int optimize_sifting(
void)
399 int optimize_concurrent(
void)
403 int optimize_deterministic_concurrent(
void)
409 int &solution_status,
double &objective_value,
410 std::vector<double> &variable_primal,
411 std::vector<double> &variable_dual,
412 std::vector<double> &variable_slack_surplus,
413 std::vector<double> &reduced_cost) = 0;
415 double &objective_value,
416 std::vector<double> &variable_primal,
417 std::vector<double> &variable_dual,
418 std::vector<double> &variable_slack_surplus,
419 std::vector<double> &reduced_cost)
423 return solve(solution_status, objective_value,
424 variable_primal, variable_dual,
425 variable_slack_surplus, reduced_cost);
428 std::vector<double> &variable_primal,
429 std::vector<double> &variable_dual,
430 std::vector<double> &variable_slack_surplus,
431 std::vector<double> &reduced_cost)
434 double objective_value;
436 return solve(solution_status, objective_value,
437 variable_primal, variable_dual,
438 variable_slack_surplus, reduced_cost);
441 int &solution_status,
double &objective_value,
442 std::vector<double> &variable_primal,
443 std::vector<double> &variable_dual)
445 std::vector<double> variable_slack_surplus;
446 std::vector<double> reduced_cost;
448 return solve(solution_status, objective_value,
449 variable_primal, variable_dual,
450 variable_slack_surplus, reduced_cost);
453 std::vector<double> &variable_primal,
454 std::vector<double> &variable_dual)
457 double objective_value;
459 return solve(solution_status, objective_value, variable_primal,
463 std::vector<double> &variable_primal)
465 std::vector<double> variable_dual;
467 return solve(variable_primal, variable_dual);
471 class environment_t {
473 std::vector<problem_t *> _problem;
476 #ifndef BPMPD_INFINITY_BOUND
477 #define BPMPD_INFINITY_BOUND 1e+30
478 #endif // BPMPD_INFINITY_BOUND
480 #ifndef BPMPD_VERSION_STRING
481 #define BPMPD_VERSION_STRING "Version 2.11 (1996 December)"
482 #endif // BPMPD_VERSION_STRING
484 class bpmpd_environment_t;
486 class bpmpd_problem_t :
public problem_t {
488 #pragma clang diagnostic push
489 #pragma clang diagnostic ignored "-Wunused-const-variable"
491 #pragma clang diagnostic pop
493 double _objective_sense;
494 double _objective_value;
495 std::vector<double> _variable_primal;
496 std::vector<double> _variable_dual;
497 int _solution_status;
499 int ampl_solution_status(
const int termination_code)
501 switch (termination_code) {
515 fprintf(stderr,
"%s:%d: error: unknown termination code "
516 "%d\n", __FILE__, __LINE__, termination_code);
519 fprintf(stderr,
"%s:%d: %d\n", __FILE__, __LINE__,
522 void set_bpmpd_parameter(
void)
607 void append_constraint_sense_bound(
608 std::vector<double> &lower_bound_bpmpd,
609 std::vector<double> &upper_bound_bpmpd)
611 for (std::vector<char>::const_iterator iterator =
612 _constraint_sense.begin();
613 iterator != _constraint_sense.end(); iterator++) {
617 upper_bound_bpmpd.push_back(0);
620 lower_bound_bpmpd.push_back(0);
621 upper_bound_bpmpd.push_back(0);
624 lower_bound_bpmpd.push_back(0);
630 iterator - _constraint_sense.begin();
632 lower_bound_bpmpd.push_back(0);
633 upper_bound_bpmpd.push_back(_range_value[index] -
639 lower_bound_bpmpd.reserve(
dims_.
mn);
640 lower_bound_bpmpd.insert(lower_bound_bpmpd.end(),
641 _lower_bound.begin(),
643 upper_bound_bpmpd.reserve(
dims_.
mn);
644 upper_bound_bpmpd.insert(upper_bound_bpmpd.end(),
645 _upper_bound.begin(),
654 bpmpd_problem_t(
void)
655 : problem_t(
false), _objective_sense(1.0),
656 _objective_value(NAN)
659 inline size_t nrow(
void)
const
663 inline size_t ncolumn(
void)
const
667 inline void set_objective_sense(
int sense)
669 _objective_sense = sense;
678 void push_back_column(
const double objective,
679 const double lower_bound = 0,
682 problem_t::push_back_column(objective, lower_bound,
685 void set_size(
const size_t nonzero_value_size,
686 const double mf = 6.0,
const size_t ma = 0)
688 dims_.
n = _objective.size();
702 std::max(static_cast<size_t>(0), ma);
709 std::vector<int> column_pointer;
711 std::vector<int> row_index;
713 std::vector<double> nonzero_value;
715 to_csr(column_pointer, row_index, nonzero_value);
718 std::bind1st(std::multiplies<double>(),
722 for (
size_t i = 0;
i < 3;
i++) {
723 set_size(nonzero_value.size(), 6.0 * (1 <<
i));
728 set_bpmpd_parameter();
730 std::vector<double> diag(
dims_.
mn, 0);
731 std::vector<double> odiag(
dims_.
mn, 0);
732 std::vector<double> dxs(
dims_.
mn, 0);
733 std::vector<double> dxsn(
dims_.
mn, 0);
735 std::vector<double> dual_residual(
dims_.
mn, 0);
736 std::vector<double> ddspr(
dims_.
mn, 0);
737 std::vector<double> ddsprn(
dims_.
mn, 0);
738 std::vector<double> dsup(
dims_.
mn, 0);
739 std::vector<double> ddsup(
dims_.
mn, 0);
740 std::vector<double> ddsupn(
dims_.
mn, 0);
741 std::vector<double> ddv(
dims_.
m, 0);
742 std::vector<double> ddvn(
dims_.
m, 0);
743 std::vector<double> prinf(
dims_.
m, 0);
744 std::vector<double> upinf(
dims_.
mn, 0);
745 std::vector<double> duinf(
dims_.
mn, 0);
747 std::vector<int> vartyp(
dims_.
n, 0);
748 std::vector<int> slktyp(
dims_.
m, 0);
749 std::vector<int> colpnt(
dims_.
n1, 0);
750 std::vector<int> ecolpnt(
dims_.
mn, 0);
752 std::vector<int> vcstat(
dims_.
mn, 0);
753 std::vector<int> pivots(
dims_.
mn, 0);
754 std::vector<int> invprm(
dims_.
mn, 0);
755 std::vector<int> snhead(
dims_.
mn, 0);
756 std::vector<int> nodtyp(
dims_.
mn, 0);
757 std::vector<int> inta1(
dims_.
mn, 0);
758 std::vector<int> prehis(
dims_.
mn, 0);
760 int termination_code;
769 double bigbou = 1
e+15;
773 _variable_primal.resize(
dims_.
mn);
774 _variable_dual.resize(
dims_.
m);
776 std::vector<double> lower_bound_bpmpd = _lower_bound;
777 std::vector<double> upper_bound_bpmpd = _upper_bound;
779 append_constraint_sense_bound(lower_bound_bpmpd,
782 solver_(&_objective[0], &_rhs[0], &lower_bound_bpmpd[0],
783 &upper_bound_bpmpd[0], &diag[0], &odiag[0],
784 &_variable_primal[0], &dxs[0], &dxsn[0], &
up[0],
785 &dual_residual[0], &ddspr[0], &ddsprn[0],
786 &dsup[0], &ddsup[0], &ddsupn[0],
787 &_variable_dual[0], &ddv[0], &ddvn[0], &prinf[0],
788 &upinf[0], &duinf[0], &
scale[0],
789 &nonzero_value[0], &vartyp[0], &slktyp[0],
790 &column_pointer[0], &ecolpnt[0], &
count[0],
791 &vcstat[0], &pivots[0], &invprm[0], &snhead[0],
792 &nodtyp[0], &inta1[0], &prehis[0], &row_index[0],
793 &rindex[0], &termination_code, &_objective_value,
794 &iter, &correct, &fixn, &dropn, &fnzmax, &fnzmin,
795 &addobj, &bigbou, &infinity_copy, &ft);
797 _objective_value *= _objective_sense;
798 _solution_status = ampl_solution_status(termination_code);
799 if (termination_code != -2) {
805 return _solution_status == 0 ? 0 : 1;
808 int &solution_status,
double &objective_value,
809 std::vector<double> &variable_primal,
810 std::vector<double> &variable_dual,
811 std::vector<double> &variable_slack_surplus,
812 std::vector<double> &reduced_cost)
820 int &solution_status,
double &objective_value,
821 std::vector<double> &variable_primal,
822 std::vector<double> &variable_dual)
824 variable_primal = std::vector<double>(
825 _variable_primal.begin(),
826 _variable_primal.begin() + _objective.size());
827 variable_dual = _variable_dual;
831 friend class bpmpd_environment_t;
834 class bpmpd_environment_t :
public environment_t {
836 bpmpd_environment_t(
void)
839 ~bpmpd_environment_t(
void)
842 int set_verbose(
void);
843 int set_data_checking(
void);
850 return bpmpd_problem_t();
858 double angular_range_reduce(
const double x)
860 if (!std::isfinite(x)) {
864 static const double cody_waite_x_max = 1608.4954386379741381;
865 static const double two_pi_0 = 6.2831853071795649157;
866 static const double two_pi_1 = 2.1561211432631314669e-14;
867 static const double two_pi_2 = 1.1615423895917441336e-27;
870 if (x >= -cody_waite_x_max && x <= cody_waite_x_max) {
871 static const double inverse_two_pi =
872 0.15915494309189534197;
873 const double k = rint(x * inverse_two_pi);
874 ret = ((x - (k * two_pi_0)) - k * two_pi_1) -
887 double hermite_h_normalized(
const size_t n,
const double x)
892 case 3: y = -3.913998411780905*x + 2.6093322745206033*
std::pow(x,3);
break;
893 case 5: y = 4.931174490213579*x - 6.574899320284771*
std::pow(x,3) + 1.3149798640569543*
std::pow(x,5);
break;
894 case 7: y = -5.773117374387059*x + 11.546234748774118*
std::pow(x,3) - 4.618493899509647*
std::pow(x,5) + 0.43985656185806166*
std::pow(x,7);
break;
895 case 9: y = 6.507479403136423*x - 17.353278408363792*
std::pow(x,3) + 10.411967045018276*
std::pow(x,5) - 1.9832318180987192*
std::pow(x,7) + 0.11017954544992885*
std::pow(x,9);
break;
896 case 11: y = -7.167191940825306*x + 23.89063980275102*
std::pow(x,3) - 19.112511842200817*
std::pow(x,5) + 5.460717669200234*
std::pow(x,7) - 0.6067464076889149*
std::pow(x,9) + 0.02206350573414236*
std::pow(x,11);
break;
897 case 13: y = 7.771206704387521*x - 31.084826817550084*
std::pow(x,3) + 31.084826817550084*
std::pow(x,5) - 11.841838787638126*
std::pow(x,7) + 1.9736397979396878*
std::pow(x,9) - 0.14353743985015913*
std::pow(x,11) + 0.0036804471756451056*
std::pow(x,13);
break;
898 case 15: y = -8.331608118589472*x + 38.88083788675087*
std::pow(x,3) - 46.65700546410104*
std::pow(x,5) + 22.217621649571925*
std::pow(x,7) - 4.9372492554604275*
std::pow(x,9) + 0.5386090096865921*
std::pow(x,11) - 0.027620974855722673*
std::pow(x,13) + 0.00052611380677567*
std::pow(x,15);
break;
899 case 17: y = 8.856659222944476*x - 47.23551585570387*
std::pow(x,3) + 66.12972219798543*
std::pow(x,5) - 37.7884126845631*
std::pow(x,7) + 10.496781301267527*
std::pow(x,9) - 1.5268045529116403*
std::pow(x,11) + 0.11744650407012618*
std::pow(x,13) - 0.004474152536004807*
std::pow(x,15) + 0.0000657963608236001*
std::pow(x,17);
break;
910 static const size_t ncms_hcal_edge_pseudorapidity = 82 + 1;
911 static const double cms_hcal_edge_pseudorapidity[
912 ncms_hcal_edge_pseudorapidity] = {
913 -5.191, -4.889, -4.716, -4.538, -4.363, -4.191, -4.013,
914 -3.839, -3.664, -3.489, -3.314, -3.139, -2.964, -2.853,
915 -2.650, -2.500, -2.322, -2.172, -2.043, -1.930, -1.830,
916 -1.740, -1.653, -1.566, -1.479, -1.392, -1.305, -1.218,
917 -1.131, -1.044, -0.957, -0.879, -0.783, -0.696, -0.609,
918 -0.522, -0.435, -0.348, -0.261, -0.174, -0.087,
920 0.087, 0.174, 0.261, 0.348, 0.435, 0.522, 0.609,
921 0.696, 0.783, 0.879, 0.957, 1.044, 1.131, 1.218,
922 1.305, 1.392, 1.479, 1.566, 1.653, 1.740, 1.830,
923 1.930, 2.043, 2.172, 2.322, 2.500, 2.650, 2.853,
924 2.964, 3.139, 3.314, 3.489, 3.664, 3.839, 4.013,
925 4.191, 4.363, 4.538, 4.716, 4.889, 5.191
929 cms_hcal_edge_pseudorapidity,
930 cms_hcal_edge_pseudorapidity +
931 ncms_hcal_edge_pseudorapidity);
933 static const size_t ncms_ecal_edge_pseudorapidity = 344 + 1;
935 for (
size_t i = 0;
i < ncms_ecal_edge_pseudorapidity;
i++) {
938 (ncms_ecal_edge_pseudorapidity - 1)) -
959 for (std::vector<particle_t>::const_iterator iterator =
961 iterator !=
_event.end(); iterator++) {
962 const unsigned int reduced_id =
963 iterator->reduced_particle_flow_id;
967 if (iterator->momentum.Eta() >=
969 iterator->momentum.Eta() <
971 const double azimuth =
972 iterator->momentum.Phi();
975 (*_perp_fourier)[
k - 1][reduced_id]
977 iterator->momentum.Pt() *
979 (*_perp_fourier)[
k - 1][reduced_id]
981 iterator->momentum.Pt() *
990 const size_t nfeature = 2 *
nfourier - 1;
1003 scale[0] = 1.0 / 5400.0;
1006 scale[1] = 1.0 / 130.0;
1009 scale[2] = 1.0 / 220.0;
1012 const size_t index_edge_end =
1018 ((*_perp_fourier)[0 ][
j][0][0] +
1019 (*_perp_fourier)[index_edge_end][
j][0][0]);
1026 ((*_perp_fourier)[0 ][
j][
k][0] +
1027 (*_perp_fourier)[index_edge_end][
j][
k][0]);
1032 ((*_perp_fourier)[0 ][
j][
k][1] +
1033 (*_perp_fourier)[index_edge_end][
j][
k][1]);
1052 #ifdef HAVE_SPARSEHASH
1054 const voronoi_diagram_t::Face face_empty;
1055 google::dense_hash_map<voronoi_diagram_t::Face_handle,
1056 size_t, hash<voronoi_diagram_t::Face_handle> >
1059 face_index.set_empty_key(face_empty);
1060 #else // HAVE_SPARSEHASH
1061 std::map<voronoi_diagram_t::Face_handle, size_t>
1063 #endif // HAVE_SPARSEHASH
1065 for (std::vector<particle_t>::const_iterator iterator =
1067 iterator !=
_event.end(); iterator++) {
1071 for (
int k = -1;
k <= 1;
k++) {
1073 iterator->momentum.Eta(),
1074 iterator->momentum.Phi() +
1076 const voronoi_diagram_t::Face_handle
handle =
1086 for (std::vector<particle_t>::iterator iterator =
1088 iterator !=
_event.end(); iterator++) {
1089 const voronoi_diagram_t::Locate_result
result =
1090 diagram.locate(*iterator);
1091 const voronoi_diagram_t::Face_handle *face =
1092 boost::get<voronoi_diagram_t::Face_handle>(
1094 double polygon_area;
1097 voronoi_diagram_t::Ccb_halfedge_circulator
1098 circulator_start = (*face)->outer_ccb();
1099 bool unbounded =
false;
1102 voronoi_diagram_t::Ccb_halfedge_circulator
1103 circulator = circulator_start;
1108 if (circulator->has_target()) {
1110 circulator->target()->point());
1111 _event[face_index[*face]].incident.
1114 face_index[circulator->twin()->
1122 while (++circulator != circulator_start);
1124 polygon_area = INFINITY;
1127 polygon_area = polygon.area();
1133 iterator->area = fabs(polygon_area);
1138 for (std::vector<particle_t>::iterator iterator =
1140 iterator !=
_event.end(); iterator++) {
1141 int predictor_index = -1;
1142 int interpolation_index = -1;
1147 if (iterator->momentum.Eta() >=
1149 iterator->momentum.Eta() <
1151 predictor_index =
l - 1;
1160 if (iterator->momentum.Eta() >=
1162 iterator->momentum.Eta() <
1164 interpolation_index =
l - 1;
1172 if (iterator->momentum.Eta() >=
1174 iterator->momentum.Eta() <
1176 interpolation_index =
l - 1;
1181 if (predictor_index >= 0 && interpolation_index >= 0) {
1185 const double azimuth = iterator->momentum.Phi();
1186 const float (*
p)[2][82] =
1188 ue_predictor_pf[
j][predictor_index]
1191 #endif // STANDALONE
1196 const size_t norder =
l == 0 ? 9 : 1;
1198 for (
size_t m = 0;
m < 2;
m++) {
1199 float u =
p[
l][
m][0];
1201 for (
size_t n = 0;
n < 2 * nfourier - 1;
n++) {
1202 if ((
l == 0 &&
n == 0) || (
l == 2 && (
n == 3 ||
n == 4))) {
1204 for (
size_t o = 2;
o < norder + 1;
o++) {
1205 u +=
p[
l][
m][9 *
n +
o] *
1206 hermite_h_normalized(
1207 2 *
o - 1, _feature[
n]) *
1208 exp(-_feature[n] * _feature[n]);
1213 pred += u * (
l == 0 ? 1.0 : 2.0) *
1214 (
m == 0 ?
cos(l * azimuth) :
1216 if (l == 0 &&
m == 0) {
1230 ue_interpolation_pf0[predictor_index][
1231 interpolation_index];
1235 ue_interpolation_pf1[predictor_index][
1236 interpolation_index];
1240 ue_interpolation_pf2[predictor_index][
1241 interpolation_index];
1247 interpolation_index];
1252 interpolation_index];
1257 interpolation_index];
1259 #endif // STANDALONE
1271 if (std::isfinite(iterator->area) && density >= 0) {
1274 iterator->momentum_perp_subtracted =
1275 iterator->momentum.Pt() -
1276 density * iterator->area;
1279 iterator->momentum_perp_subtracted =
1280 iterator->momentum.Pt();
1282 iterator->momentum_perp_subtracted_unequalized =
1283 iterator->momentum_perp_subtracted;
1288 boost::multi_array<double, 2> radial_distance_square(
1291 for (std::vector<particle_t>::const_iterator
1292 iterator_outer =
_event.begin();
1293 iterator_outer !=
_event.end(); iterator_outer++) {
1294 radial_distance_square
1295 [iterator_outer -
_event.begin()]
1296 [iterator_outer -
_event.begin()] = 0;
1298 for (std::vector<particle_t>::const_iterator
1299 iterator_inner =
_event.begin();
1300 iterator_inner != iterator_outer;
1302 const double deta = iterator_outer->momentum.Eta() -
1303 iterator_inner->momentum.Eta();
1304 const double dphi = angular_range_reduce(
1305 iterator_outer->momentum.Phi() -
1306 iterator_inner->momentum.Phi());
1308 radial_distance_square
1309 [iterator_outer -
_event.begin()]
1310 [iterator_inner -
_event.begin()] =
1311 deta * deta + dphi * dphi;
1312 radial_distance_square
1313 [iterator_inner -
_event.begin()]
1314 [iterator_outer -
_event.begin()] =
1315 radial_distance_square
1316 [iterator_outer -
_event.begin()]
1317 [iterator_inner -
_event.begin()];
1323 for (std::vector<particle_t>::const_iterator
1324 iterator_outer =
_event.begin();
1325 iterator_outer !=
_event.end(); iterator_outer++) {
1326 double incident_area_sum = iterator_outer->area;
1328 for (std::set<std::vector<particle_t>::iterator>::
1329 const_iterator iterator_inner =
1330 iterator_outer->incident.begin();
1332 iterator_outer->incident.end();
1334 incident_area_sum += (*iterator_inner)->area;
1336 _active.push_back(incident_area_sum < 2.0);
1341 _event.size(), std::vector<size_t>());
1343 _event.size(), std::vector<size_t>());
1349 static const size_t npair_max = 36;
1351 for (
size_t i = 0;
i <
_event.size();
i++) {
1352 for (
size_t j = 0;
j <
_event.size();
j++) {
1354 const size_t incident_count =
1359 (radial_distance_square[
i][
j] <
1361 incident_count > 0)) {
1366 if (
_event[
i].momentum_perp_subtracted < 0) {
1367 std::vector<double> radial_distance_square_list;
1369 for (std::vector<size_t>::const_iterator iterator =
1373 const size_t j = *iterator;
1375 if (
_event[j].momentum_perp_subtracted > 0) {
1376 radial_distance_square_list.push_back(
1377 radial_distance_square[
i][j]);
1381 double radial_distance_square_max_equalization_cut =
1384 if (radial_distance_square_list.size() >= npair_max) {
1385 std::sort(radial_distance_square_list.begin(),
1386 radial_distance_square_list.end());
1387 radial_distance_square_max_equalization_cut =
1388 radial_distance_square_list[npair_max - 1];
1391 for (std::vector<size_t>::const_iterator iterator =
1395 const size_t j = *iterator;
1397 if (
_event[j].momentum_perp_subtracted > 0 &&
1398 radial_distance_square[
i][j] <
1399 radial_distance_square_max_equalization_cut) {
1405 std::pair<size_t, size_t>(
i, j));
1407 radial_distance_square[
i][j] /
1416 bpmpd_problem_t *
p =
reinterpret_cast<bpmpd_problem_t *
>(lp_problem);
1437 p->set_objective_sense(bpmpd_problem_t::minimize);
1441 static const size_t nsector_azimuth = 12;
1446 static const size_t ncms_hcal_edge_pseudorapidity = 19 + 1;
1447 static const double cms_hcal_edge_pseudorapidity[
1448 ncms_hcal_edge_pseudorapidity] = {
1449 -5.191, -4.538, -4.013,
1450 -3.489, -2.853, -2.322, -1.830, -1.305, -0.783, -0.261,
1451 0.261, 0.783, 1.305, 1.830, 2.322, 2.853, 3.489,
1456 const double *edge_pseudorapidity;
1458 nedge_pseudorapidity = ncms_hcal_edge_pseudorapidity;
1459 edge_pseudorapidity = cms_hcal_edge_pseudorapidity;
1461 const size_t nsuperblock = (nedge_pseudorapidity - 2) *
1464 size_t index_row = 0;
1465 for (
size_t index_pseudorapidity = 0;
1466 index_pseudorapidity < nedge_pseudorapidity - 2;
1467 index_pseudorapidity++) {
1468 for (
size_t index_azimuth = 0;
1469 index_azimuth < nsector_azimuth - 1;
1471 const size_t index_column =
1472 index_pseudorapidity * nsector_azimuth +
1475 bpmpd_problem_t::greater_equal, 0);
1476 p->push_back_coefficient(
1477 index_row, index_column, 1);
1478 p->push_back_coefficient(
1479 index_row, nsuperblock + index_column, -1);
1482 bpmpd_problem_t::greater_equal, 0);
1483 p->push_back_coefficient(
1484 index_row, index_column, 1);
1485 p->push_back_coefficient(
1486 index_row, nsuperblock + index_column + 1, -1);
1489 bpmpd_problem_t::greater_equal, 0);
1490 p->push_back_coefficient(
1491 index_row, index_column, 1);
1492 p->push_back_coefficient(
1494 nsuperblock + index_column + nsector_azimuth, -1);
1497 bpmpd_problem_t::greater_equal, 0);
1498 p->push_back_coefficient(
1499 index_row, index_column, 1);
1500 p->push_back_coefficient(
1502 nsuperblock + index_column + nsector_azimuth + 1,
1506 const size_t index_column =
1507 index_pseudorapidity * nsector_azimuth +
1508 nsector_azimuth - 1;
1510 bpmpd_problem_t::greater_equal, 0);
1511 p->push_back_coefficient(
1512 index_row, index_column, 1);
1513 p->push_back_coefficient(
1514 index_row, nsuperblock + index_column, -1);
1517 bpmpd_problem_t::greater_equal, 0);
1518 p->push_back_coefficient(
1519 index_row, index_column, 1);
1520 p->push_back_coefficient(
1522 nsuperblock + index_column - (nsector_azimuth - 1),
1526 bpmpd_problem_t::greater_equal, 0);
1527 p->push_back_coefficient(
1528 index_row, index_column, 1);
1529 p->push_back_coefficient(
1531 nsuperblock + index_column + nsector_azimuth, -1);
1534 bpmpd_problem_t::greater_equal, 0);
1535 p->push_back_coefficient(
1536 index_row, index_column, 1);
1537 p->push_back_coefficient(
1539 nsuperblock + index_column + nsector_azimuth -
1540 (nsector_azimuth - 1),
1545 const size_t nstaggered_block =
1546 (nedge_pseudorapidity - 1) * nsector_azimuth;
1547 const size_t nblock = nsuperblock + 2 * nstaggered_block;
1553 size_t positive_count = 0;
1555 for (std::vector<particle_t>::const_iterator iterator =
1557 iterator !=
_event.end(); iterator++) {
1558 if (iterator->momentum_perp_subtracted >= 0) {
1559 positive_index[iterator -
_event.begin()] =
1565 _ncost = nblock + positive_count;
1572 std::vector<particle_t>::const_iterator
1573 iterator_particle =
_event.begin();
1574 std::vector<bool>::const_iterator iterator_active =
1576 std::vector<std::vector<size_t> >::const_iterator
1577 iterator_recombine_index_outer =
1579 std::vector<std::vector<size_t> >::const_iterator
1580 iterator_recombine_unsigned_outer =
1582 size_t index_column_max =
_ncost - 1;
1583 for (; iterator_particle !=
_event.end();
1584 iterator_particle++, iterator_active++,
1585 iterator_recombine_index_outer++,
1586 iterator_recombine_unsigned_outer++) {
1587 if (*iterator_active) {
1588 int index_pseudorapidity = -1;
1592 if (iterator_particle->momentum.Eta() >= edge_pseudorapidity[
i - 1] &&
1593 iterator_particle->momentum.Eta() < edge_pseudorapidity[
i]) {
1594 index_pseudorapidity =
i - 1;
1598 const int index_azimuth = floor(
1599 (iterator_particle->momentum.Phi() +
M_PI) *
1600 ((nsector_azimuth >> 1) /
M_PI));
1602 if (index_pseudorapidity != -1) {
1610 iterator_particle->momentum_perp_subtracted >= 0 ?
1612 bpmpd_problem_t::greater_equal,
1613 iterator_particle->momentum_perp_subtracted);
1616 const double sign = iterator_particle->momentum_perp_subtracted >= 0 ? 1 : -1;
1617 const size_t index_column_block_subtract =
1619 (nedge_pseudorapidity - 1) * nsector_azimuth +
1620 index_pseudorapidity * nsector_azimuth +
1624 index_column_block_subtract;
1626 if (iterator_particle->momentum_perp_subtracted >= 0) {
1627 const size_t index_column_cost =
1628 nblock + positive_index[iterator_particle -
_event.begin()];
1630 p->push_back_coefficient(
1631 index_row, index_column_cost, 1);
1633 std::max(index_column_max, index_column_cost);
1635 p->push_back_coefficient(
1636 index_row, index_column_block_subtract, 1);
1638 std::max(index_column_max, index_column_block_subtract);
1640 for (std::vector<size_t>::const_iterator
1641 iterator_recombine_index_inner =
1642 iterator_recombine_index_outer->begin();
1643 iterator_recombine_index_inner !=
1644 iterator_recombine_index_outer->end();
1645 iterator_recombine_index_inner++) {
1646 const size_t index_column =
1647 *iterator_recombine_index_inner +
1650 p->push_back_coefficient(
1651 index_row, index_column, sign);
1653 std::max(index_column_max, index_column);
1657 const size_t index_column_block =
1659 index_pseudorapidity * nsector_azimuth +
1667 double sum_unequalized;
1669 sum_unequalized = 0;
1670 for (std::vector<size_t>::const_iterator
1671 iterator_recombine_unsigned_inner =
1672 iterator_recombine_unsigned_outer->begin();
1673 iterator_recombine_unsigned_inner !=
1674 iterator_recombine_unsigned_outer->end();
1675 iterator_recombine_unsigned_inner++) {
1677 _event[*iterator_recombine_unsigned_inner].momentum_perp_subtracted;
1679 sum_unequalized =
std::max(0.0, sum_unequalized);
1681 if (sum_unequalized >= sum_unequalized_3 ||
1682 (sum_unequalized >= sum_unequalized_2 &&
1683 (iterator_particle -
_event.begin()) % 2 == 0) ||
1684 (sum_unequalized >= sum_unequalized_1 &&
1685 (iterator_particle -
_event.begin()) % 4 == 0) ||
1686 (sum_unequalized >= sum_unequalized_0 &&
1687 (iterator_particle -
_event.begin()) % 8 == 0)) {
1689 const double weight = sum_unequalized *
1691 iterator_particle->area));
1695 bpmpd_problem_t::greater_equal,
1698 p->push_back_coefficient(
1699 index_row, index_column_block, 1.0 / weight);
1701 for (std::vector<size_t>::const_iterator
1702 iterator_recombine_unsigned_inner =
1703 iterator_recombine_unsigned_outer->begin();
1704 iterator_recombine_unsigned_inner !=
1705 iterator_recombine_unsigned_outer->end();
1706 iterator_recombine_unsigned_inner++) {
1707 if (
_event[*iterator_recombine_unsigned_inner].momentum_perp_subtracted >= 0) {
1708 const size_t index_column_cost =
1710 positive_index[*iterator_recombine_unsigned_inner];
1712 p->push_back_coefficient(
1713 index_row, index_column_cost, 1);
1715 std::max(index_column_max, index_column_cost);
1721 bpmpd_problem_t::greater_equal,
1724 p->push_back_coefficient(
1727 for (std::vector<size_t>::const_iterator iterator_recombine_unsigned_inner = iterator_recombine_unsigned_outer->begin();
1728 iterator_recombine_unsigned_inner != iterator_recombine_unsigned_outer->end();
1729 iterator_recombine_unsigned_inner++) {
1730 if (
_event[*iterator_recombine_unsigned_inner].momentum_perp_subtracted >= 0) {
1731 const size_t index_column_cost =
1733 positive_index[*iterator_recombine_unsigned_inner];
1735 p->push_back_coefficient(
1736 index_row, index_column_cost, -1);
1738 std::max(index_column_max, index_column_cost);
1752 static const double epsilon_degeneracy = 1
e-2;
1758 for (
size_t i = 0;
i < nsuperblock;
i++) {
1759 p->push_back_column(
1762 for (
size_t i = nsuperblock;
i < nsuperblock + nstaggered_block;
i++) {
1763 p->push_back_column(
1766 for (
size_t i = nsuperblock + nstaggered_block;
i < nsuperblock + 2 * nstaggered_block;
i++) {
1767 p->push_back_column(
1770 for (
size_t i = nsuperblock + 2 * nstaggered_block;
i <
_ncost;
i++) {
1771 p->push_back_column(
1777 for (
size_t i = _ncost;
i <= index_column_max;
i++) {
1778 p->push_back_column(
1785 bpmpd_problem_t lp_problem =
reinterpret_cast<bpmpd_environment_t *
>(
_lp_environment)->problem();
1789 lp_problem.optimize();
1791 int solution_status;
1792 double objective_value;
1793 std::vector<double>
x;
1794 std::vector<double>
pi;
1796 lp_problem.solve(solution_status, objective_value,
1799 for (
size_t k =
_ncost;
k < x.size();
k++) {
1801 momentum_perp_subtracted < 0 &&
1803 momentum_perp_subtracted >= 0 && x[
k] >= 0) {
1805 momentum_perp_subtracted += x[
k];
1807 momentum_perp_subtracted -= x[
k];
1810 for (
size_t k = 0;
k <
_event.size();
k++) {
1813 _event[
k].momentum_perp_subtracted -=
1820 for (std::vector<particle_t>::iterator iterator =
1822 iterator !=
_event.end(); iterator++) {
1823 iterator->momentum_perp_subtracted =
std::max(
1824 0.0, iterator->momentum_perp_subtracted);
1844 const double dr_max,
1845 const std::pair<double, double> equalization_threshold,
1846 const bool remove_nonpositive)
1847 : _remove_nonpositive(remove_nonpositive),
1848 _equalization_threshold(equalization_threshold),
1849 _radial_distance_square_max(dr_max * dr_max),
1850 _positive_bound_scale(0.2),
1858 -5.191, -2.650, -2.043, -1.740, -1.479, -1.131, -0.783, -0.522, 0.522, 0.783, 1.131, 1.479, 1.740, 2.043, 2.650, 5.191
1862 edge_pseudorapidity,
1882 const double perp,
const double pseudorapidity,
1883 const double azimuth,
1884 const unsigned int reduced_particle_flow_id)
1908 std::vector<double>
ret;
1910 for (std::vector<particle_t>::const_iterator iterator =
1912 iterator !=
_event.end(); iterator++) {
1913 ret.push_back(iterator->momentum_perp_subtracted);
1922 std::vector<double>
ret;
1924 for (std::vector<particle_t>::const_iterator iterator =
1926 iterator !=
_event.end(); iterator++) {
1927 ret.push_back(iterator->momentum_perp_subtracted_unequalized);
1942 std::vector<double>
ret;
1944 for (std::vector<particle_t>::const_iterator iterator =
1946 iterator !=
_event.end(); iterator++) {
1947 ret.push_back(iterator->area);
1964 std::vector<std::set<size_t> >
ret;
1966 for (std::vector<particle_t>::const_iterator
1967 iterator_outer =
_event.begin();
1968 iterator_outer !=
_event.end(); iterator_outer++) {
1971 for (std::set<std::vector<particle_t>::iterator>::
1972 const_iterator iterator_inner =
1973 iterator_outer->incident.begin();
1974 iterator_inner != iterator_outer->incident.begin();
1976 e.insert(*iterator_inner -
_event.begin());
1987 return std::vector<double>(
std::vector< bool > _active
std::vector< std::vector< size_t > > _recombine_index
tuple ret
prodAgent to be discontinued
std::vector< double > perp_fourier(void)
void subtract_momentum(void)
Sin< T >::type sin(const T &t)
CGAL::Voronoi_diagram_2< delaunay_triangulation_t, CGAL::Delaunay_triangulation_adaptation_traits_2< delaunay_triangulation_t >, CGAL::Delaunay_triangulation_caching_degeneracy_removal_policy_2< delaunay_triangulation_t > > voronoi_diagram_t
double _positive_bound_scale
void remove_nonpositive(void)
std::vector< double > _edge_pseudorapidity
void initialize_geometry(void)
std::vector< double > _recombine_tie
bool equal(const T &first, const T &second)
int solver_(double *obj, double *rhs, double *lbound, double *ubound, double *diag, double *odiag, double *xs, double *dxs, double *dxsn, double *up, double *dspr, double *ddspr, double *ddsprn, double *dsup, double *ddsup, double *ddsupn, double *dv, double *ddv, double *ddvn, double *prinf, double *upinf, double *duinf, double *scale, double *nonzeros, int *vartyp, int *slktyp, int *colpnt, int *ecolpnt, int *count, int *vcstat, int *pivots, int *invprm, int *snhead, int *nodtyp, int *inta1, int *prehis, int *rowidx, int *rindex, int *code, double *opt, int *iter, int *corect, int *fixn, int *dropn, int *fnzmax, int *fnzmin, double *addobj, double *bigbou, double *big, int *ft)
U second(std::pair< T, U > const &p)
T x() const
Cartesian x coordinate.
delaunay_triangulation_t::Point point_2d_t
std::vector< double > _cms_ecal_edge_pseudorapidity
std::vector< double > particle_area(void)
void lp_populate(void *lp_problem)
void subtract_if_necessary(void)
void push_back_particle(const double perp, const double pseudorapidity, const double azimuth, const unsigned int reduced_particle_flow_id)
std::vector< std::pair< size_t, size_t > > _recombine
#define BPMPD_VERSION_STRING
void clear(CLHEP::HepGenMatrix &m)
Helper function: Reset all elements of a matrix to 0.
static const size_t nreduced_particle_flow_id
std::vector< std::set< size_t > > particle_incident(void)
void recombine_link(void)
Cos< T >::type cos(const T &t)
static const size_t nfourier
std::vector< size_t > _nblock_subtract
std::vector< double > _feature
bool insert(Storage &iStorage, ItemType *iItem, const IdTag &iIdTag)
double _radial_distance_square_max
boost::multi_array< double, 4 > * _perp_fourier
PtEtaPhiELorentzVectorD PtEtaPhiELorentzVector
Lorentz vector with cartesian internal representation.
std::vector< std::vector< size_t > > _recombine_unsigned
float ue_interpolation_pf0[15][344]
void feature_extract(void)
Geom::Phi< T > phi() const
VoronoiAlgorithm(const UECalibration *ue, const double dr_max, const std::pair< double, double > equalization_threshold=std::pair< double, double >(5.0, 35.0), const bool remove_nonpositive=true)
T perp() const
Magnitude of transverse component.
void voronoi_area_incident(void)
std::pair< double, double > _equalization_threshold
size_t nedge_pseudorapidity(void) const
float ue_predictor_pf[3][15][5][2][82]
std::vector< double > _cms_hcal_edge_pseudorapidity
volatile std::atomic< bool > shutdown_flag false
CGAL::Polygon_2< CGAL::Exact_predicates_inexact_constructions_kernel > polygon_t
std::vector< double > subtracted_unequalized_perp(void)
float ue_interpolation_pf1[15][344]
float ue_interpolation_pf2[15][82]
std::vector< double > subtracted_equalized_perp(void)
Power< A, B >::type pow(const A &a, const B &b)
#define BPMPD_INFINITY_BOUND