File indexing completed on 2024-09-27 04:56:57
0001 #include "SimTransport/TotemRPProtonTransportParametrization/interface/LHCOpticsApproximator.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 #include <vector>
0005 #include <iostream>
0006 #include "TROOT.h"
0007 #include "TFile.h"
0008 #include <memory>
0009 #include "TMatrixD.h"
0010 #include "TMath.h"
0011
0012 ClassImp(LHCOpticsApproximator);
0013 ClassImp(LHCApertureApproximator);
0014
0015 void LHCOpticsApproximator::Init() {
0016 out_polynomials.clear();
0017 apertures_.clear();
0018 out_polynomials.push_back(&x_parametrisation);
0019 out_polynomials.push_back(&theta_x_parametrisation);
0020 out_polynomials.push_back(&y_parametrisation);
0021 out_polynomials.push_back(&theta_y_parametrisation);
0022
0023 coord_names.clear();
0024 coord_names.push_back("x");
0025 coord_names.push_back("theta_x");
0026 coord_names.push_back("y");
0027 coord_names.push_back("theta_y");
0028 coord_names.push_back("ksi");
0029
0030 s_begin_ = 0.0;
0031 s_end_ = 0.0;
0032 trained_ = false;
0033 }
0034
0035 LHCOpticsApproximator::LHCOpticsApproximator(std::string name,
0036 std::string title,
0037 TMultiDimFet::EMDFPolyType polynom_type,
0038 std::string beam_direction,
0039 double nominal_beam_momentum)
0040 : x_parametrisation(5, polynom_type, "k"),
0041 theta_x_parametrisation(5, polynom_type, "k"),
0042 y_parametrisation(5, polynom_type, "k"),
0043 theta_y_parametrisation(5, polynom_type, "k") {
0044 this->SetName(name.c_str());
0045 this->SetTitle(title.c_str());
0046 Init();
0047
0048 if (beam_direction == "lhcb1")
0049 beam = lhcb1;
0050 else if (beam_direction == "lhcb2")
0051 beam = lhcb2;
0052 else
0053 beam = lhcb1;
0054
0055 nominal_beam_momentum_ = nominal_beam_momentum;
0056 nominal_beam_energy_ = TMath::Sqrt(nominal_beam_momentum_ * nominal_beam_momentum_ + 0.938272029 * 0.938272029);
0057 }
0058
0059 LHCOpticsApproximator::LHCOpticsApproximator() {
0060 Init();
0061 beam = lhcb1;
0062 nominal_beam_momentum_ = 7000;
0063 nominal_beam_energy_ = TMath::Sqrt(nominal_beam_momentum_ * nominal_beam_momentum_ + 0.938272029 * 0.938272029);
0064 }
0065
0066
0067
0068 double LHCOpticsApproximator::ParameterOutOfRangePenalty(double in[], bool invert_beam_coord_sytems) const {
0069 double in_corrected[5];
0070 if (beam == lhcb1 || !invert_beam_coord_sytems) {
0071 in_corrected[0] = in[0];
0072 in_corrected[1] = in[1];
0073 in_corrected[2] = in[2];
0074 in_corrected[3] = in[3];
0075 in_corrected[4] = in[4];
0076 } else {
0077 in_corrected[0] = -in[0];
0078 in_corrected[1] = -in[1];
0079 in_corrected[2] = in[2];
0080 in_corrected[3] = in[3];
0081 in_corrected[4] = in[4];
0082 }
0083
0084 const TVectorD *min_var = x_parametrisation.GetMinVariables();
0085 const TVectorD *max_var = x_parametrisation.GetMaxVariables();
0086 double res = 0.;
0087
0088 for (int i = 0; i < 5; i++) {
0089 if (in_corrected[i] < (*min_var)(i)) {
0090 double dist = TMath::Abs(((*min_var)(i)-in_corrected[i]) / ((*max_var)(i) - (*min_var)(i)));
0091 res += 8 * (TMath::Exp(dist) - 1.0);
0092 in_corrected[i] = (*min_var)(i);
0093 } else if (in_corrected[i] > (*max_var)(i)) {
0094 double dist = TMath::Abs((in_corrected[i] - (*max_var)(i)) / ((*max_var)(i) - (*min_var)(i)));
0095 res += 8 * (TMath::Exp(dist) - 1.0);
0096 in_corrected[i] = (*max_var)(i);
0097 }
0098 }
0099 return res;
0100 }
0101
0102 bool LHCOpticsApproximator::Transport(const double *in,
0103 double *out,
0104 bool check_apertures,
0105 bool invert_beam_coord_sytems) const {
0106 if (in == nullptr || out == nullptr || !trained_)
0107 return false;
0108
0109 bool res = CheckInputRange(in);
0110 double in_corrected[5];
0111
0112 if (beam == lhcb1 || !invert_beam_coord_sytems) {
0113 in_corrected[0] = in[0];
0114 in_corrected[1] = in[1];
0115 in_corrected[2] = in[2];
0116 in_corrected[3] = in[3];
0117 in_corrected[4] = in[4];
0118 out[0] = x_parametrisation.Eval(in_corrected);
0119 out[1] = theta_x_parametrisation.Eval(in_corrected);
0120 out[2] = y_parametrisation.Eval(in_corrected);
0121 out[3] = theta_y_parametrisation.Eval(in_corrected);
0122 out[4] = in[4];
0123 } else {
0124 in_corrected[0] = -in[0];
0125 in_corrected[1] = -in[1];
0126 in_corrected[2] = in[2];
0127 in_corrected[3] = in[3];
0128 in_corrected[4] = in[4];
0129 out[0] = -x_parametrisation.Eval(in_corrected);
0130 out[1] = -theta_x_parametrisation.Eval(in_corrected);
0131 out[2] = y_parametrisation.Eval(in_corrected);
0132 out[3] = theta_y_parametrisation.Eval(in_corrected);
0133 out[4] = in[4];
0134 }
0135
0136 if (check_apertures) {
0137 for (unsigned int i = 0; i < apertures_.size(); i++) {
0138 res = res && apertures_[i].CheckAperture(in);
0139 }
0140 }
0141 return res;
0142 }
0143
0144 bool LHCOpticsApproximator::Transport2D(const double *in,
0145 double *out,
0146 bool check_apertures,
0147 bool invert_beam_coord_sytems) const
0148
0149 {
0150 if (in == nullptr || out == nullptr || !trained_)
0151 return false;
0152
0153 bool res = CheckInputRange(in);
0154 double in_corrected[5];
0155
0156 if (beam == lhcb1 || !invert_beam_coord_sytems) {
0157 in_corrected[0] = in[0];
0158 in_corrected[1] = in[1];
0159 in_corrected[2] = in[2];
0160 in_corrected[3] = in[3];
0161 in_corrected[4] = in[4];
0162 out[0] = x_parametrisation.Eval(in_corrected);
0163 out[1] = y_parametrisation.Eval(in_corrected);
0164 } else {
0165 in_corrected[0] = -in[0];
0166 in_corrected[1] = -in[1];
0167 in_corrected[2] = in[2];
0168 in_corrected[3] = in[3];
0169 in_corrected[4] = in[4];
0170 out[0] = -x_parametrisation.Eval(in_corrected);
0171 out[1] = y_parametrisation.Eval(in_corrected);
0172 }
0173
0174 if (check_apertures) {
0175 for (unsigned int i = 0; i < apertures_.size(); i++) {
0176 res = res && apertures_[i].CheckAperture(in);
0177 }
0178 }
0179 return res;
0180 }
0181
0182 bool LHCOpticsApproximator::Transport_m_GeV(double in_pos[3],
0183 double in_momentum[3],
0184 double out_pos[3],
0185 double out_momentum[3],
0186 bool check_apertures,
0187 double z2_z1_dist) const {
0188 double in[5];
0189 double out[5];
0190 double part_mom = 0.0;
0191 for (int i = 0; i < 3; ++i)
0192 part_mom += in_momentum[i] * in_momentum[i];
0193
0194 part_mom = std::sqrt(part_mom);
0195
0196 in[0] = in_pos[0];
0197 in[1] = in_momentum[0] / nominal_beam_momentum_;
0198 in[2] = in_pos[1];
0199 in[3] = in_momentum[1] / nominal_beam_momentum_;
0200 in[4] = (part_mom - nominal_beam_momentum_) / nominal_beam_momentum_;
0201
0202 if (!Transport(in, out, check_apertures, true)) {
0203 return false;
0204 }
0205
0206 out_pos[0] = out[0];
0207 out_pos[1] = out[2];
0208 out_pos[2] = in_pos[2] + z2_z1_dist;
0209
0210 out_momentum[0] = out[1] * nominal_beam_momentum_;
0211 out_momentum[1] = out[3] * nominal_beam_momentum_;
0212 double part_out_total_mom = (out[4] + 1) * nominal_beam_momentum_;
0213 out_momentum[2] = std::sqrt(part_out_total_mom * part_out_total_mom - out_momentum[0] * out_momentum[0] -
0214 out_momentum[1] * out_momentum[1]);
0215 out_momentum[2] = TMath::Sign(out_momentum[2], in_momentum[2]);
0216
0217 return true;
0218 }
0219
0220 bool LHCOpticsApproximator::Transport(const MadKinematicDescriptor *in,
0221 MadKinematicDescriptor *out,
0222 bool check_apertures,
0223 bool invert_beam_coord_sytems) const {
0224 if (in == nullptr || out == nullptr || !trained_)
0225 return false;
0226
0227 Double_t input[5];
0228 Double_t output[5];
0229 input[0] = in->x;
0230 input[1] = in->theta_x;
0231 input[2] = in->y;
0232 input[3] = in->theta_y;
0233 input[4] = in->ksi;
0234
0235
0236 if (!Transport(input, output, check_apertures, invert_beam_coord_sytems)) {
0237 return false;
0238 }
0239
0240 out->x = output[0];
0241 out->theta_x = output[1];
0242 out->y = output[2];
0243 out->theta_y = output[3];
0244 out->ksi = output[4];
0245
0246 return true;
0247 }
0248
0249 LHCOpticsApproximator::LHCOpticsApproximator(const LHCOpticsApproximator &org)
0250 : TNamed(org),
0251 x_parametrisation(org.x_parametrisation),
0252 theta_x_parametrisation(org.theta_x_parametrisation),
0253 y_parametrisation(org.y_parametrisation),
0254 theta_y_parametrisation(org.theta_y_parametrisation) {
0255 Init();
0256 s_begin_ = org.s_begin_;
0257 s_end_ = org.s_end_;
0258 trained_ = org.trained_;
0259 apertures_ = org.apertures_;
0260 beam = org.beam;
0261 nominal_beam_energy_ = org.nominal_beam_energy_;
0262 nominal_beam_momentum_ = org.nominal_beam_momentum_;
0263 }
0264
0265 const LHCOpticsApproximator &LHCOpticsApproximator::operator=(const LHCOpticsApproximator &org) {
0266 if (this != &org) {
0267 x_parametrisation = org.x_parametrisation;
0268 theta_x_parametrisation = org.theta_x_parametrisation;
0269 y_parametrisation = org.y_parametrisation;
0270 theta_y_parametrisation = org.theta_y_parametrisation;
0271 Init();
0272
0273 TNamed::operator=(org);
0274 s_begin_ = org.s_begin_;
0275 s_end_ = org.s_end_;
0276 trained_ = org.trained_;
0277
0278 apertures_ = org.apertures_;
0279 beam = org.beam;
0280 nominal_beam_energy_ = org.nominal_beam_energy_;
0281 nominal_beam_momentum_ = org.nominal_beam_momentum_;
0282 }
0283 return org;
0284 }
0285
0286 void LHCOpticsApproximator::Train(TTree *inp_tree,
0287 std::string data_prefix,
0288 polynomials_selection mode,
0289 int max_degree_x,
0290 int max_degree_tx,
0291 int max_degree_y,
0292 int max_degree_ty,
0293 bool common_terms,
0294 double *prec) {
0295 if (inp_tree == nullptr)
0296 return;
0297
0298 InitializeApproximators(mode, max_degree_x, max_degree_tx, max_degree_y, max_degree_ty, common_terms);
0299
0300
0301
0302 double in_var[6];
0303
0304
0305
0306 double out_var[7];
0307
0308
0309 std::string x_in_lab = "x_in";
0310 std::string theta_x_in_lab = "theta_x_in";
0311 std::string y_in_lab = "y_in";
0312 std::string theta_y_in_lab = "theta_y_in";
0313 std::string ksi_in_lab = "ksi_in";
0314 std::string s_in_lab = "s_in";
0315
0316 std::string x_out_lab = data_prefix + "_x_out";
0317 std::string theta_x_out_lab = data_prefix + "_theta_x_out";
0318 std::string y_out_lab = data_prefix + "_y_out";
0319 std::string theta_y_out_lab = data_prefix + "_theta_y_out";
0320 std::string ksi_out_lab = data_prefix + "_ksi_out";
0321 std::string s_out_lab = data_prefix + "_s_out";
0322 std::string valid_out_lab = data_prefix + "_valid_out";
0323
0324
0325 inp_tree->SetBranchStatus("*", false);
0326 inp_tree->SetBranchStatus(x_in_lab.c_str(), true);
0327 inp_tree->SetBranchStatus(theta_x_in_lab.c_str(), true);
0328 inp_tree->SetBranchStatus(y_in_lab.c_str(), true);
0329 inp_tree->SetBranchStatus(theta_y_in_lab.c_str(), true);
0330 inp_tree->SetBranchStatus(ksi_in_lab.c_str(), true);
0331 inp_tree->SetBranchStatus(x_out_lab.c_str(), true);
0332 inp_tree->SetBranchStatus(theta_x_out_lab.c_str(), true);
0333 inp_tree->SetBranchStatus(y_out_lab.c_str(), true);
0334 inp_tree->SetBranchStatus(theta_y_out_lab.c_str(), true);
0335 inp_tree->SetBranchStatus(ksi_out_lab.c_str(), true);
0336 inp_tree->SetBranchStatus(valid_out_lab.c_str(), true);
0337
0338
0339 inp_tree->SetBranchAddress(x_in_lab.c_str(), &(in_var[0]));
0340 inp_tree->SetBranchAddress(theta_x_in_lab.c_str(), &(in_var[1]));
0341 inp_tree->SetBranchAddress(y_in_lab.c_str(), &(in_var[2]));
0342 inp_tree->SetBranchAddress(theta_y_in_lab.c_str(), &(in_var[3]));
0343 inp_tree->SetBranchAddress(ksi_in_lab.c_str(), &(in_var[4]));
0344 inp_tree->SetBranchAddress(s_in_lab.c_str(), &(in_var[5]));
0345
0346
0347 inp_tree->SetBranchAddress(x_out_lab.c_str(), &(out_var[0]));
0348 inp_tree->SetBranchAddress(theta_x_out_lab.c_str(), &(out_var[1]));
0349 inp_tree->SetBranchAddress(y_out_lab.c_str(), &(out_var[2]));
0350 inp_tree->SetBranchAddress(theta_y_out_lab.c_str(), &(out_var[3]));
0351 inp_tree->SetBranchAddress(ksi_out_lab.c_str(), &(out_var[4]));
0352 inp_tree->SetBranchAddress(s_out_lab.c_str(), &(out_var[5]));
0353 inp_tree->SetBranchAddress(valid_out_lab.c_str(), &(out_var[6]));
0354
0355 Long64_t entries = inp_tree->GetEntries();
0356 if (entries > 0) {
0357 inp_tree->SetBranchStatus(s_in_lab.c_str(), true);
0358 inp_tree->SetBranchStatus(s_out_lab.c_str(), true);
0359 inp_tree->GetEntry(0);
0360 s_begin_ = in_var[5];
0361 s_end_ = out_var[5];
0362 inp_tree->SetBranchStatus(s_in_lab.c_str(), false);
0363 inp_tree->SetBranchStatus(s_out_lab.c_str(), false);
0364 }
0365
0366
0367 for (Long64_t i = 0; i < entries; ++i) {
0368 inp_tree->GetEntry(i);
0369 if (out_var[6] != 0)
0370 {
0371 x_parametrisation.AddRow(in_var, out_var[0], 0);
0372 theta_x_parametrisation.AddRow(in_var, out_var[1], 0);
0373 y_parametrisation.AddRow(in_var, out_var[2], 0);
0374 theta_y_parametrisation.AddRow(in_var, out_var[3], 0);
0375 }
0376 }
0377
0378 edm::LogInfo("LHCOpticsApproximator") << "Optical functions parametrizations from " << s_begin_ << " to " << s_end_
0379 << "\n";
0380 PrintInputRange();
0381 for (int i = 0; i < 4; i++) {
0382 double best_precision = 0.0;
0383 if (prec)
0384 best_precision = prec[i];
0385 out_polynomials[i]->FindParameterization(best_precision);
0386 edm::LogInfo("LHCOpticsApproximator") << "Out variable " << coord_names[i] << " polynomial"
0387 << "\n";
0388 out_polynomials[i]->PrintPolynomialsSpecial("M");
0389 edm::LogInfo("LHCOpticsApproximator") << "\n";
0390 }
0391
0392 trained_ = true;
0393 }
0394
0395 void LHCOpticsApproximator::InitializeApproximators(polynomials_selection mode,
0396 int max_degree_x,
0397 int max_degree_tx,
0398 int max_degree_y,
0399 int max_degree_ty,
0400 bool common_terms) {
0401 SetDefaultAproximatorSettings(x_parametrisation, VariableType::X, max_degree_x);
0402 SetDefaultAproximatorSettings(theta_x_parametrisation, VariableType::THETA_X, max_degree_tx);
0403 SetDefaultAproximatorSettings(y_parametrisation, VariableType::Y, max_degree_y);
0404 SetDefaultAproximatorSettings(theta_y_parametrisation, VariableType::THETA_Y, max_degree_ty);
0405
0406 if (mode == PREDEFINED) {
0407 SetTermsManually(x_parametrisation, VariableType::X, max_degree_x, common_terms);
0408 SetTermsManually(theta_x_parametrisation, VariableType::THETA_X, max_degree_tx, common_terms);
0409 SetTermsManually(y_parametrisation, VariableType::Y, max_degree_y, common_terms);
0410 SetTermsManually(theta_y_parametrisation, VariableType::THETA_Y, max_degree_ty, common_terms);
0411 }
0412 }
0413
0414 void LHCOpticsApproximator::SetDefaultAproximatorSettings(TMultiDimFet &approximator,
0415 VariableType var_type,
0416 int max_degree) {
0417 if (max_degree < 1 || max_degree > 20)
0418 max_degree = 10;
0419
0420 if (var_type == VariableType::X || var_type == VariableType::THETA_X) {
0421 Int_t mPowers[] = {1, 1, 0, 0, max_degree};
0422 approximator.SetMaxPowers(mPowers);
0423 approximator.SetMaxFunctions(900);
0424 approximator.SetMaxStudy(1000);
0425 approximator.SetMaxTerms(900);
0426 approximator.SetPowerLimit(1.50);
0427 approximator.SetMinAngle(10);
0428 approximator.SetMaxAngle(10);
0429 approximator.SetMinRelativeError(1e-13);
0430 }
0431
0432 if (var_type == VariableType::Y || var_type == VariableType::THETA_Y) {
0433 Int_t mPowers[] = {0, 0, 1, 1, max_degree};
0434 approximator.SetMaxPowers(mPowers);
0435 approximator.SetMaxFunctions(900);
0436 approximator.SetMaxStudy(1000);
0437 approximator.SetMaxTerms(900);
0438 approximator.SetPowerLimit(1.50);
0439 approximator.SetMinAngle(10);
0440 approximator.SetMaxAngle(10);
0441 approximator.SetMinRelativeError(1e-13);
0442 }
0443 }
0444
0445 void LHCOpticsApproximator::SetTermsManually(TMultiDimFet &approximator,
0446 VariableType variable,
0447 int max_degree,
0448 bool common_terms) {
0449 if (max_degree < 1 || max_degree > 20)
0450 max_degree = 10;
0451
0452
0453
0454
0455
0456 std::vector<Int_t> term_literals;
0457 term_literals.reserve(5000);
0458
0459 if (variable == VariableType::X || variable == VariableType::THETA_X) {
0460
0461 for (int i = 0; i <= max_degree; ++i) {
0462 term_literals.push_back(1);
0463 term_literals.push_back(0);
0464 term_literals.push_back(0);
0465 term_literals.push_back(0);
0466 term_literals.push_back(i);
0467 }
0468
0469 for (int i = 0; i <= max_degree; ++i) {
0470 term_literals.push_back(0);
0471 term_literals.push_back(1);
0472 term_literals.push_back(0);
0473 term_literals.push_back(0);
0474 term_literals.push_back(i);
0475 }
0476
0477 for (int i = 0; i <= max_degree; ++i) {
0478 term_literals.push_back(0);
0479 term_literals.push_back(2);
0480 term_literals.push_back(0);
0481 term_literals.push_back(0);
0482 term_literals.push_back(i);
0483 }
0484
0485 for (int i = 0; i <= max_degree; ++i) {
0486 term_literals.push_back(0);
0487 term_literals.push_back(3);
0488 term_literals.push_back(0);
0489 term_literals.push_back(0);
0490 term_literals.push_back(i);
0491 }
0492
0493 for (int i = 0; i <= max_degree; ++i) {
0494 term_literals.push_back(0);
0495 term_literals.push_back(0);
0496 term_literals.push_back(0);
0497 term_literals.push_back(0);
0498 term_literals.push_back(i);
0499 }
0500 }
0501
0502 if (variable == VariableType::Y || variable == VariableType::THETA_Y) {
0503
0504 for (int i = 0; i <= max_degree; ++i) {
0505 term_literals.push_back(0);
0506 term_literals.push_back(0);
0507 term_literals.push_back(1);
0508 term_literals.push_back(0);
0509 term_literals.push_back(i);
0510 }
0511
0512 for (int i = 0; i <= max_degree; ++i) {
0513 term_literals.push_back(0);
0514 term_literals.push_back(0);
0515 term_literals.push_back(0);
0516 term_literals.push_back(1);
0517 term_literals.push_back(i);
0518 }
0519
0520 for (int i = 0; i <= max_degree; ++i) {
0521 term_literals.push_back(0);
0522 term_literals.push_back(0);
0523 term_literals.push_back(0);
0524 term_literals.push_back(2);
0525 term_literals.push_back(i);
0526 }
0527
0528 for (int i = 0; i <= max_degree; ++i) {
0529 term_literals.push_back(0);
0530 term_literals.push_back(0);
0531 term_literals.push_back(0);
0532 term_literals.push_back(3);
0533 term_literals.push_back(i);
0534 }
0535
0536 for (int i = 0; i <= max_degree; ++i) {
0537 term_literals.push_back(0);
0538 term_literals.push_back(0);
0539 term_literals.push_back(0);
0540 term_literals.push_back(0);
0541 term_literals.push_back(i);
0542 }
0543 }
0544
0545
0546 if (common_terms) {
0547 term_literals.push_back(1), term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(0),
0548 term_literals.push_back(0);
0549 term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(0),
0550 term_literals.push_back(0);
0551 term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(1),
0552 term_literals.push_back(0);
0553 term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(1), term_literals.push_back(0),
0554 term_literals.push_back(0);
0555 term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(1),
0556 term_literals.push_back(0);
0557 term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(1),
0558 term_literals.push_back(0);
0559
0560 term_literals.push_back(1), term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(0),
0561 term_literals.push_back(1);
0562 term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(0),
0563 term_literals.push_back(1);
0564 term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(1),
0565 term_literals.push_back(1);
0566 term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(1), term_literals.push_back(0),
0567 term_literals.push_back(1);
0568 term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(1),
0569 term_literals.push_back(1);
0570 term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(1),
0571 term_literals.push_back(1);
0572
0573 term_literals.push_back(1), term_literals.push_back(2), term_literals.push_back(0), term_literals.push_back(0),
0574 term_literals.push_back(0);
0575 term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(2), term_literals.push_back(0),
0576 term_literals.push_back(0);
0577 term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(2),
0578 term_literals.push_back(0);
0579 term_literals.push_back(2), term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(0),
0580 term_literals.push_back(0);
0581 term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(2), term_literals.push_back(0),
0582 term_literals.push_back(0);
0583 term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(2),
0584 term_literals.push_back(0);
0585 term_literals.push_back(2), term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(0),
0586 term_literals.push_back(0);
0587 term_literals.push_back(0), term_literals.push_back(2), term_literals.push_back(1), term_literals.push_back(0),
0588 term_literals.push_back(0);
0589 term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(2),
0590 term_literals.push_back(0);
0591 term_literals.push_back(2), term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(1),
0592 term_literals.push_back(0);
0593 term_literals.push_back(0), term_literals.push_back(2), term_literals.push_back(0), term_literals.push_back(1),
0594 term_literals.push_back(0);
0595 term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(2), term_literals.push_back(1),
0596 term_literals.push_back(0);
0597
0598 term_literals.push_back(1), term_literals.push_back(2), term_literals.push_back(0), term_literals.push_back(0),
0599 term_literals.push_back(1);
0600 term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(2), term_literals.push_back(0),
0601 term_literals.push_back(1);
0602 term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(2),
0603 term_literals.push_back(1);
0604 term_literals.push_back(2), term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(0),
0605 term_literals.push_back(1);
0606 term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(2), term_literals.push_back(0),
0607 term_literals.push_back(1);
0608 term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(2),
0609 term_literals.push_back(1);
0610 term_literals.push_back(2), term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(0),
0611 term_literals.push_back(1);
0612 term_literals.push_back(0), term_literals.push_back(2), term_literals.push_back(1), term_literals.push_back(0),
0613 term_literals.push_back(1);
0614 term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(2),
0615 term_literals.push_back(1);
0616 term_literals.push_back(2), term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(1),
0617 term_literals.push_back(1);
0618 term_literals.push_back(0), term_literals.push_back(2), term_literals.push_back(0), term_literals.push_back(1),
0619 term_literals.push_back(1);
0620 term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(2), term_literals.push_back(1),
0621 term_literals.push_back(1);
0622 }
0623
0624 std::vector<Int_t> powers;
0625 powers.resize(term_literals.size());
0626
0627 for (unsigned int i = 0; i < term_literals.size(); ++i) {
0628 powers[i] = term_literals[i];
0629 }
0630 approximator.SetPowers(&powers[0], term_literals.size() / 5);
0631 }
0632
0633 void LHCOpticsApproximator::Test(TTree *inp_tree, TFile *f_out, std::string data_prefix, std::string base_out_dir) {
0634 if (inp_tree == nullptr || f_out == nullptr)
0635 return;
0636
0637
0638
0639 double in_var[6];
0640
0641
0642
0643 double out_var[7];
0644
0645
0646 std::string x_in_lab = "x_in";
0647 std::string theta_x_in_lab = "theta_x_in";
0648 std::string y_in_lab = "y_in";
0649 std::string theta_y_in_lab = "theta_y_in";
0650 std::string ksi_in_lab = "ksi_in";
0651 std::string s_in_lab = "s_in";
0652
0653 std::string x_out_lab = data_prefix + "_x_out";
0654 std::string theta_x_out_lab = data_prefix + "_theta_x_out";
0655 std::string y_out_lab = data_prefix + "_y_out";
0656 std::string theta_y_out_lab = data_prefix + "_theta_y_out";
0657 std::string ksi_out_lab = data_prefix + "_ksi_out";
0658 std::string s_out_lab = data_prefix + "_s_out";
0659 std::string valid_out_lab = data_prefix + "_valid_out";
0660
0661
0662 inp_tree->SetBranchStatus("*", false);
0663 inp_tree->SetBranchStatus(x_in_lab.c_str(), true);
0664 inp_tree->SetBranchStatus(theta_x_in_lab.c_str(), true);
0665 inp_tree->SetBranchStatus(y_in_lab.c_str(), true);
0666 inp_tree->SetBranchStatus(theta_y_in_lab.c_str(), true);
0667 inp_tree->SetBranchStatus(ksi_in_lab.c_str(), true);
0668 inp_tree->SetBranchStatus(x_out_lab.c_str(), true);
0669 inp_tree->SetBranchStatus(theta_x_out_lab.c_str(), true);
0670 inp_tree->SetBranchStatus(y_out_lab.c_str(), true);
0671 inp_tree->SetBranchStatus(theta_y_out_lab.c_str(), true);
0672 inp_tree->SetBranchStatus(ksi_out_lab.c_str(), true);
0673 inp_tree->SetBranchStatus(valid_out_lab.c_str(), true);
0674
0675
0676 inp_tree->SetBranchAddress(x_in_lab.c_str(), &(in_var[0]));
0677 inp_tree->SetBranchAddress(theta_x_in_lab.c_str(), &(in_var[1]));
0678 inp_tree->SetBranchAddress(y_in_lab.c_str(), &(in_var[2]));
0679 inp_tree->SetBranchAddress(theta_y_in_lab.c_str(), &(in_var[3]));
0680 inp_tree->SetBranchAddress(ksi_in_lab.c_str(), &(in_var[4]));
0681 inp_tree->SetBranchAddress(s_in_lab.c_str(), &(in_var[5]));
0682
0683
0684 inp_tree->SetBranchAddress(x_out_lab.c_str(), &(out_var[0]));
0685 inp_tree->SetBranchAddress(theta_x_out_lab.c_str(), &(out_var[1]));
0686 inp_tree->SetBranchAddress(y_out_lab.c_str(), &(out_var[2]));
0687 inp_tree->SetBranchAddress(theta_y_out_lab.c_str(), &(out_var[3]));
0688 inp_tree->SetBranchAddress(ksi_out_lab.c_str(), &(out_var[4]));
0689 inp_tree->SetBranchAddress(s_out_lab.c_str(), &(out_var[5]));
0690 inp_tree->SetBranchAddress(valid_out_lab.c_str(), &(out_var[6]));
0691
0692
0693 TH1D *err_hists[4];
0694 TH2D *err_inp_cor_hists[4][5];
0695 TH2D *err_out_cor_hists[4][5];
0696
0697 AllocateErrorHists(err_hists);
0698 AllocateErrorInputCorHists(err_inp_cor_hists);
0699 AllocateErrorOutputCorHists(err_out_cor_hists);
0700
0701 Long64_t entries = inp_tree->GetEntries();
0702
0703 for (Long64_t i = 0; i < entries; ++i) {
0704 double errors[4];
0705 inp_tree->GetEntry(i);
0706
0707 errors[0] = out_var[0] - x_parametrisation.Eval(in_var);
0708 errors[1] = out_var[1] - theta_x_parametrisation.Eval(in_var);
0709 errors[2] = out_var[2] - y_parametrisation.Eval(in_var);
0710 errors[3] = out_var[3] - theta_y_parametrisation.Eval(in_var);
0711
0712 FillErrorHistograms(errors, err_hists);
0713 FillErrorDataCorHistograms(errors, in_var, err_inp_cor_hists);
0714 FillErrorDataCorHistograms(errors, out_var, err_out_cor_hists);
0715 }
0716
0717 WriteHistograms(err_hists, err_inp_cor_hists, err_out_cor_hists, f_out, base_out_dir);
0718
0719 DeleteErrorHists(err_hists);
0720 DeleteErrorCorHistograms(err_inp_cor_hists);
0721 DeleteErrorCorHistograms(err_out_cor_hists);
0722 }
0723
0724 void LHCOpticsApproximator::AllocateErrorHists(TH1D *err_hists[4]) {
0725 std::vector<std::string> error_labels;
0726 error_labels.push_back("x error");
0727 error_labels.push_back("theta_x error");
0728 error_labels.push_back("y error");
0729 error_labels.push_back("theta_y error");
0730
0731 for (int i = 0; i < 4; ++i) {
0732 err_hists[i] = new TH1D(error_labels[i].c_str(), error_labels[i].c_str(), 100, -0.0000000001, 0.0000000001);
0733 err_hists[i]->SetXTitle(error_labels[i].c_str());
0734 err_hists[i]->SetYTitle("counts");
0735 err_hists[i]->SetDirectory(nullptr);
0736 err_hists[i]->SetCanExtend(TH1::kAllAxes);
0737 }
0738 }
0739
0740 void LHCOpticsApproximator::TestAperture(
0741 TTree *inp_tree, TTree *out_tree)
0742 {
0743 if (inp_tree == nullptr || out_tree == nullptr)
0744 return;
0745
0746 Long64_t entries = inp_tree->GetEntries();
0747 double entry[7];
0748 double parametrization_out[5];
0749
0750 inp_tree->SetBranchAddress("x", &(entry[0]));
0751 inp_tree->SetBranchAddress("theta_x", &(entry[1]));
0752 inp_tree->SetBranchAddress("y", &(entry[2]));
0753 inp_tree->SetBranchAddress("theta_y", &(entry[3]));
0754 inp_tree->SetBranchAddress("ksi", &(entry[4]));
0755 inp_tree->SetBranchAddress("mad_accept", &(entry[5]));
0756 inp_tree->SetBranchAddress("par_accept", &(entry[6]));
0757
0758 out_tree->SetBranchAddress("x", &(entry[0]));
0759 out_tree->SetBranchAddress("theta_x", &(entry[1]));
0760 out_tree->SetBranchAddress("y", &(entry[2]));
0761 out_tree->SetBranchAddress("theta_y", &(entry[3]));
0762 out_tree->SetBranchAddress("ksi", &(entry[4]));
0763 out_tree->SetBranchAddress("mad_accept", &(entry[5]));
0764 out_tree->SetBranchAddress("par_accept", &(entry[6]));
0765
0766
0767 for (Long64_t i = 0; i < entries; ++i) {
0768 inp_tree->GetEntry(i);
0769
0770
0771
0772 bool res = Transport(entry, parametrization_out, true, false);
0773
0774 if (res)
0775 entry[6] = 1.0;
0776 else
0777 entry[6] = 0.0;
0778
0779 out_tree->Fill();
0780 }
0781 }
0782
0783 void LHCOpticsApproximator::AllocateErrorInputCorHists(TH2D *err_inp_cor_hists[4][5]) {
0784 std::vector<std::string> error_labels;
0785 std::vector<std::string> data_labels;
0786
0787 error_labels.push_back("x error");
0788 error_labels.push_back("theta_x error");
0789 error_labels.push_back("y error");
0790 error_labels.push_back("theta_y error");
0791
0792 data_labels.push_back("x input");
0793 data_labels.push_back("theta_x input");
0794 data_labels.push_back("y input");
0795 data_labels.push_back("theta_y input");
0796 data_labels.push_back("ksi input");
0797
0798 for (int eri = 0; eri < 4; ++eri) {
0799 for (int dati = 0; dati < 5; ++dati) {
0800 std::string name = error_labels[eri] + " vs. " + data_labels[dati];
0801 const std::string &title = name;
0802 err_inp_cor_hists[eri][dati] =
0803 new TH2D(name.c_str(), title.c_str(), 100, -0.0000000001, 0.0000000001, 100, -0.0000000001, 0.0000000001);
0804 err_inp_cor_hists[eri][dati]->SetXTitle(error_labels[eri].c_str());
0805 err_inp_cor_hists[eri][dati]->SetYTitle(data_labels[dati].c_str());
0806 err_inp_cor_hists[eri][dati]->SetDirectory(nullptr);
0807 err_inp_cor_hists[eri][dati]->SetCanExtend(TH1::kAllAxes);
0808 }
0809 }
0810 }
0811
0812 void LHCOpticsApproximator::AllocateErrorOutputCorHists(TH2D *err_out_cor_hists[4][5]) {
0813 std::vector<std::string> error_labels;
0814 std::vector<std::string> data_labels;
0815
0816 error_labels.push_back("x error");
0817 error_labels.push_back("theta_x error");
0818 error_labels.push_back("y error");
0819 error_labels.push_back("theta_y error");
0820
0821 data_labels.push_back("x output");
0822 data_labels.push_back("theta_x output");
0823 data_labels.push_back("y output");
0824 data_labels.push_back("theta_y output");
0825 data_labels.push_back("ksi output");
0826
0827 for (int eri = 0; eri < 4; ++eri) {
0828 for (int dati = 0; dati < 5; ++dati) {
0829 std::string name = error_labels[eri] + " vs. " + data_labels[dati];
0830 const std::string &title = name;
0831 err_out_cor_hists[eri][dati] =
0832 new TH2D(name.c_str(), title.c_str(), 100, -0.0000000001, 0.0000000001, 100, -0.0000000001, 0.0000000001);
0833 err_out_cor_hists[eri][dati]->SetXTitle(error_labels[eri].c_str());
0834 err_out_cor_hists[eri][dati]->SetYTitle(data_labels[dati].c_str());
0835 err_out_cor_hists[eri][dati]->SetDirectory(nullptr);
0836 err_out_cor_hists[eri][dati]->SetCanExtend(TH1::kAllAxes);
0837 }
0838 }
0839 }
0840
0841 void LHCOpticsApproximator::FillErrorHistograms(double errors[4], TH1D *err_hists[4]) {
0842 for (int i = 0; i < 4; ++i) {
0843 err_hists[i]->Fill(errors[i]);
0844 }
0845 }
0846
0847 void LHCOpticsApproximator::FillErrorDataCorHistograms(double errors[4], double var[5], TH2D *err_cor_hists[4][5]) {
0848 for (int eri = 0; eri < 4; ++eri) {
0849 for (int dati = 0; dati < 5; ++dati) {
0850 err_cor_hists[eri][dati]->Fill(errors[eri], var[dati]);
0851 }
0852 }
0853 }
0854
0855 void LHCOpticsApproximator::DeleteErrorHists(TH1D *err_hists[4]) {
0856 for (int i = 0; i < 4; ++i) {
0857 delete err_hists[i];
0858 }
0859 }
0860
0861 void LHCOpticsApproximator::DeleteErrorCorHistograms(TH2D *err_cor_hists[4][5]) {
0862 for (int eri = 0; eri < 4; ++eri) {
0863 for (int dati = 0; dati < 5; ++dati) {
0864 delete err_cor_hists[eri][dati];
0865 }
0866 }
0867 }
0868
0869 void LHCOpticsApproximator::WriteHistograms(TH1D *err_hists[4],
0870 TH2D *err_inp_cor_hists[4][5],
0871 TH2D *err_out_cor_hists[4][5],
0872 TFile *f_out,
0873 std::string base_out_dir) {
0874 if (f_out == nullptr)
0875 return;
0876
0877 f_out->cd();
0878 if (!gDirectory->cd(base_out_dir.c_str()))
0879 gDirectory->mkdir(base_out_dir.c_str());
0880
0881 gDirectory->cd(base_out_dir.c_str());
0882 gDirectory->mkdir(this->GetName());
0883 gDirectory->cd(this->GetName());
0884 gDirectory->mkdir("x");
0885 gDirectory->mkdir("theta_x");
0886 gDirectory->mkdir("y");
0887 gDirectory->mkdir("theta_y");
0888
0889 gDirectory->cd("x");
0890 err_hists[0]->Write("", TObject::kWriteDelete);
0891 for (int i = 0; i < 5; i++) {
0892 err_inp_cor_hists[0][i]->Write("", TObject::kWriteDelete);
0893 err_out_cor_hists[0][i]->Write("", TObject::kWriteDelete);
0894 }
0895
0896 gDirectory->cd("..");
0897 gDirectory->cd("theta_x");
0898 err_hists[1]->Write("", TObject::kWriteDelete);
0899 for (int i = 0; i < 5; i++) {
0900 err_inp_cor_hists[1][i]->Write("", TObject::kWriteDelete);
0901 err_out_cor_hists[1][i]->Write("", TObject::kWriteDelete);
0902 }
0903
0904 gDirectory->cd("..");
0905 gDirectory->cd("y");
0906 err_hists[2]->Write("", TObject::kWriteDelete);
0907 for (int i = 0; i < 5; i++) {
0908 err_inp_cor_hists[2][i]->Write("", TObject::kWriteDelete);
0909 err_out_cor_hists[2][i]->Write("", TObject::kWriteDelete);
0910 }
0911
0912 gDirectory->cd("..");
0913 gDirectory->cd("theta_y");
0914 err_hists[3]->Write("", TObject::kWriteDelete);
0915 for (int i = 0; i < 5; i++) {
0916 err_inp_cor_hists[3][i]->Write("", TObject::kWriteDelete);
0917 err_out_cor_hists[3][i]->Write("", TObject::kWriteDelete);
0918 }
0919 gDirectory->cd("..");
0920 gDirectory->cd("..");
0921 }
0922
0923 void LHCOpticsApproximator::PrintInputRange() {
0924 const TVectorD *min_var = x_parametrisation.GetMinVariables();
0925 const TVectorD *max_var = x_parametrisation.GetMaxVariables();
0926
0927 edm::LogInfo("LHCOpticsApproximator") << "Covered input parameters range:"
0928 << "\n";
0929 for (int i = 0; i < 5; i++) {
0930 edm::LogInfo("LHCOpticsApproximator") << (*min_var)(i) << " < " << coord_names[i] << " < " << (*max_var)(i) << "\n";
0931 }
0932 edm::LogInfo("LHCOpticsApproximator") << "\n";
0933 }
0934
0935 bool LHCOpticsApproximator::CheckInputRange(const double *in,
0936 bool invert_beam_coord_sytems) const
0937 {
0938 double in_corrected[5];
0939 if (beam == lhcb1 || !invert_beam_coord_sytems) {
0940 in_corrected[0] = in[0];
0941 in_corrected[1] = in[1];
0942 in_corrected[2] = in[2];
0943 in_corrected[3] = in[3];
0944 in_corrected[4] = in[4];
0945 } else {
0946 in_corrected[0] = -in[0];
0947 in_corrected[1] = -in[1];
0948 in_corrected[2] = in[2];
0949 in_corrected[3] = in[3];
0950 in_corrected[4] = in[4];
0951 }
0952
0953 const TVectorD *min_var = x_parametrisation.GetMinVariables();
0954 const TVectorD *max_var = x_parametrisation.GetMaxVariables();
0955 bool res = true;
0956
0957 for (int i = 0; i < 5; i++) {
0958 res = res && in_corrected[i] >= (*min_var)(i) && in_corrected[i] <= (*max_var)(i);
0959 }
0960 return res;
0961 }
0962
0963 void LHCOpticsApproximator::AddRectEllipseAperture(
0964 const LHCOpticsApproximator &in, double rect_x, double rect_y, double r_el_x, double r_el_y) {
0965 apertures_.push_back(
0966 LHCApertureApproximator(in, rect_x, rect_y, r_el_x, r_el_y, LHCApertureApproximator::ApertureType::RECTELLIPSE));
0967 }
0968
0969 LHCApertureApproximator::LHCApertureApproximator() {
0970 rect_x_ = rect_y_ = r_el_x_ = r_el_y_ = 0.0;
0971 ap_type_ = ApertureType::NO_APERTURE;
0972 }
0973
0974 LHCApertureApproximator::LHCApertureApproximator(
0975 const LHCOpticsApproximator &in, double rect_x, double rect_y, double r_el_x, double r_el_y, ApertureType type)
0976 : LHCOpticsApproximator(in) {
0977 rect_x_ = rect_x;
0978 rect_y_ = rect_y;
0979 r_el_x_ = r_el_x;
0980 r_el_y_ = r_el_y;
0981 ap_type_ = type;
0982 }
0983
0984 bool LHCApertureApproximator::CheckAperture(const double *in,
0985 bool invert_beam_coord_sytems) const
0986 {
0987 double out[5];
0988 bool result = Transport(in, out, false, invert_beam_coord_sytems);
0989
0990 if (result && ap_type_ == ApertureType::RECTELLIPSE) {
0991 result = result && out[0] < rect_x_ && out[0] > -rect_x_ && out[2] < rect_y_ && out[2] > -rect_y_ &&
0992 (out[0] * out[0] / (r_el_x_ * r_el_x_) + out[2] * out[2] / (r_el_y_ * r_el_y_) < 1);
0993 }
0994 return result;
0995 }
0996
0997 void LHCOpticsApproximator::PrintOpticalFunctions() {
0998 edm::LogInfo("LHCOpticsApproximator") << std::endl
0999 << "Linear terms of optical functions:"
1000 << "\n";
1001 for (int i = 0; i < 4; i++) {
1002 PrintCoordinateOpticalFunctions(*out_polynomials[i], coord_names[i], coord_names);
1003 }
1004 }
1005
1006 void LHCOpticsApproximator::PrintCoordinateOpticalFunctions(TMultiDimFet ¶metrization,
1007 const std::string &coord_name,
1008 const std::vector<std::string> &input_vars) {
1009 double in[5];
1010
1011 double d_out_d_in[5];
1012 double d_par = 1e-5;
1013 double bias = 0;
1014
1015 for (int j = 0; j < 5; j++)
1016 in[j] = 0.0;
1017
1018 bias = parametrization.Eval(in);
1019
1020 for (int i = 0; i < 5; i++) {
1021 for (int j = 0; j < 5; j++)
1022 in[j] = 0.0;
1023
1024 in[i] = d_par;
1025 d_out_d_in[i] = parametrization.Eval(in);
1026 in[i] = 0.0;
1027 d_out_d_in[i] = d_out_d_in[i] - parametrization.Eval(in);
1028 d_out_d_in[i] = d_out_d_in[i] / d_par;
1029 }
1030 edm::LogInfo("LHCOpticsApproximator") << coord_name << " = " << bias;
1031 for (int i = 0; i < 5; i++) {
1032 edm::LogInfo("LHCOpticsApproximator") << " + " << d_out_d_in[i] << "*" << input_vars[i];
1033 }
1034 edm::LogInfo("LHCOpticsApproximator") << "\n";
1035 }
1036
1037 void LHCOpticsApproximator::GetLinearApproximation(
1038 double atPoint[], double &Cx, double &Lx, double &vx, double &Cy, double &Ly, double &vy, double &D, double ep) {
1039 double out[2];
1040 if (!Transport2D(atPoint, out)) {
1041 return;
1042 }
1043 Cx = out[0];
1044 Cy = out[1];
1045
1046 for (int i = 0; i < 5; ++i) {
1047 atPoint[i] += ep;
1048 if (!Transport2D(atPoint, out)) {
1049 continue;
1050 }
1051 switch (i) {
1052 case 0:
1053 vx = (out[0] - Cx) / ep;
1054 break;
1055 case 1:
1056 Lx = (out[0] - Cx) / ep;
1057 break;
1058 case 2:
1059 vy = (out[1] - Cy) / ep;
1060 break;
1061 case 3:
1062 Ly = (out[1] - Cy) / ep;
1063 break;
1064 case 4:
1065 D = (out[0] - Cx) / ep;
1066 break;
1067 }
1068 atPoint[i] -= ep;
1069 }
1070 }
1071
1072
1073 void LHCOpticsApproximator::GetLineariasedTransportMatrixX(double mad_init_x,
1074 double mad_init_thx,
1075 double mad_init_y,
1076 double mad_init_thy,
1077 double mad_init_xi,
1078 TMatrixD &transp_matrix,
1079 double d_mad_x,
1080 double d_mad_thx) {
1081 double MADX_momentum_correction_factor = 1.0 + mad_init_xi;
1082 transp_matrix.ResizeTo(2, 2);
1083 double in[5];
1084 in[0] = mad_init_x;
1085 in[1] = mad_init_thx;
1086 in[2] = mad_init_y;
1087 in[3] = mad_init_thy;
1088 in[4] = mad_init_xi;
1089
1090 double out[5];
1091
1092 if (!Transport(in, out)) {
1093 return;
1094 };
1095 double x1 = out[0];
1096 double thx1 = out[1];
1097
1098 in[0] = mad_init_x + d_mad_x;
1099 if (!Transport(in, out)) {
1100 return;
1101 }
1102 double x2_dx = out[0];
1103 double thx2_dx = out[1];
1104
1105 in[0] = mad_init_x;
1106 in[1] = mad_init_thx + d_mad_thx;
1107 if (!Transport(in, out)) {
1108 return;
1109 }
1110 double x2_dthx = out[0];
1111 double thx2_dthx = out[1];
1112
1113
1114
1115
1116 transp_matrix(0, 0) = (x2_dx - x1) / d_mad_x;
1117 transp_matrix(1, 0) = (thx2_dx - thx1) / (d_mad_x * MADX_momentum_correction_factor);
1118 transp_matrix(0, 1) = MADX_momentum_correction_factor * (x2_dthx - x1) / d_mad_thx;
1119 transp_matrix(1, 1) = (thx2_dthx - thx1) / d_mad_thx;
1120 }
1121
1122
1123 void LHCOpticsApproximator::GetLineariasedTransportMatrixY(double mad_init_x,
1124 double mad_init_thx,
1125 double mad_init_y,
1126 double mad_init_thy,
1127 double mad_init_xi,
1128 TMatrixD &transp_matrix,
1129 double d_mad_y,
1130 double d_mad_thy) {
1131 double MADX_momentum_correction_factor = 1.0 + mad_init_xi;
1132 transp_matrix.ResizeTo(2, 2);
1133 double in[5];
1134 in[0] = mad_init_x;
1135 in[1] = mad_init_thx;
1136 in[2] = mad_init_y;
1137 in[3] = mad_init_thy;
1138 in[4] = mad_init_xi;
1139
1140 double out[5];
1141
1142 if (!Transport(in, out)) {
1143 return;
1144 };
1145 double y1 = out[2];
1146 double thy1 = out[3];
1147
1148 in[2] = mad_init_y + d_mad_y;
1149 if (!Transport(in, out)) {
1150 return;
1151 }
1152 double y2_dy = out[2];
1153 double thy2_dy = out[3];
1154
1155 in[2] = mad_init_y;
1156 in[3] = mad_init_thy + d_mad_thy;
1157 if (!Transport(in, out)) {
1158 return;
1159 }
1160 double y2_dthy = out[2];
1161 double thy2_dthy = out[3];
1162
1163
1164
1165
1166 transp_matrix(0, 0) = (y2_dy - y1) / d_mad_y;
1167 transp_matrix(1, 0) = (thy2_dy - thy1) / (d_mad_y * MADX_momentum_correction_factor);
1168 transp_matrix(0, 1) = MADX_momentum_correction_factor * (y2_dthy - y1) / d_mad_thy;
1169 transp_matrix(1, 1) = (thy2_dthy - thy1) / d_mad_thy;
1170 }
1171
1172
1173 double LHCOpticsApproximator::GetDx(double mad_init_x,
1174 double mad_init_thx,
1175 double mad_init_y,
1176 double mad_init_thy,
1177 double mad_init_xi,
1178 double d_mad_xi) {
1179 double in[5];
1180 in[0] = mad_init_x;
1181 in[1] = mad_init_thx;
1182 in[2] = mad_init_y;
1183 in[3] = mad_init_thy;
1184 in[4] = mad_init_xi;
1185
1186 double out[5];
1187
1188 if (!Transport(in, out)) {
1189 return 0.0;
1190 }
1191 double x1 = out[0];
1192
1193 in[4] = mad_init_xi + d_mad_xi;
1194 if (!Transport(in, out)) {
1195 return 0.0;
1196 }
1197 double x2_dxi = out[0];
1198 double dispersion = (x2_dxi - x1) / d_mad_xi;
1199
1200 return dispersion;
1201 }
1202
1203
1204
1205 double LHCOpticsApproximator::GetDxds(double mad_init_x,
1206 double mad_init_thx,
1207 double mad_init_y,
1208 double mad_init_thy,
1209 double mad_init_xi,
1210 double d_mad_xi) {
1211 double MADX_momentum_correction_factor = 1.0 + mad_init_xi;
1212 double in[5];
1213 in[0] = mad_init_x;
1214 in[1] = mad_init_thx;
1215 in[2] = mad_init_y;
1216 in[3] = mad_init_thy;
1217 in[4] = mad_init_xi;
1218
1219 double out[5];
1220
1221 if (!Transport(in, out)) {
1222 return 0.0;
1223 }
1224 double thx1 = out[1] / MADX_momentum_correction_factor;
1225
1226 in[4] = mad_init_xi + d_mad_xi;
1227 if (!Transport(in, out)) {
1228 return 0.0;
1229 }
1230 double thx2_dxi = out[1] / MADX_momentum_correction_factor;
1231 double dispersion = (thx2_dxi - thx1) / d_mad_xi;
1232
1233 return dispersion;
1234 }