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