Back to home page

Project CMSSW displayed by LXR
 
 

    

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 //MADX canonical variables
0067 //(x, theta_x, y, theta_y, ksi) [m, rad, m, rad, -1..0]
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 //return true if transport possible, double x, y
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   //transport inverts the coordinate systems
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   //in-variables
0301   //x_in, theta_x_in, y_in, theta_y_in, ksi_in, s_in
0302   double in_var[6];
0303 
0304   //out-variables
0305   //x_out, theta_x_out, y_out, theta_y_out, ksi_out, s_out, valid_out;
0306   double out_var[7];
0307 
0308   //in- out-lables
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   //disable not needed branches to speed up the readin
0325   inp_tree->SetBranchStatus("*", false);  //disable all branches
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   //set input data adresses
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   //set output data adresses
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   //set input and output variables for fitting
0367   for (Long64_t i = 0; i < entries; ++i) {
0368     inp_tree->GetEntry(i);
0369     if (out_var[6] != 0)  //if out data valid
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   //put terms of shape:
0453   //1,0,0,0,t    0,1,0,0,t    0,2,0,0,t    0,3,0,0,t    0,0,0,0,t
0454   //t: 0,1,...,max_degree
0455 
0456   std::vector<Int_t> term_literals;
0457   term_literals.reserve(5000);
0458 
0459   if (variable == VariableType::X || variable == VariableType::THETA_X) {
0460     //1,0,0,0,t
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     //0,1,0,0,t
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     //0,2,0,0,t
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     //0,3,0,0,t
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     //0,0,0,0,t
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     //0,0,1,0,t
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     //0,0,0,1,t
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     //0,0,0,2,t
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     //0,0,0,3,t
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     //0,0,0,0,t
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   //push common terms
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   //in-variables
0638   //x_in, theta_x_in, y_in, theta_y_in, ksi_in, s_in
0639   double in_var[6];
0640 
0641   //out-variables
0642   //x_out, theta_x_out, y_out, theta_y_out, ksi_out, s_out, valid_out;
0643   double out_var[7];
0644 
0645   //in- out-lables
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   //disable not needed branches to speed up the readin
0662   inp_tree->SetBranchStatus("*", false);  //disable all branches
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   //set input data adresses
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   //set output data adresses
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   //test histogramms
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   //set input and output variables for fitting
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)  //x, theta_x, y, theta_y, ksi, mad_accepted, parametriz_accepted
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   //  int ind=0;
0767   for (Long64_t i = 0; i < entries; ++i) {
0768     inp_tree->GetEntry(i);
0769 
0770     //Don't invert the coordinate systems, appertures are defined in the
0771     //coordinate system of the beam - perhaps to be changed
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  //x, thx, y, thy, ksi
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  //x, thx. y, thy, ksi
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 &parametrization,
1007                                                             const std::string &coord_name,
1008                                                             const std::vector<std::string> &input_vars) {
1009   double in[5];
1010   //  double out;
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 //real angles in the matrix, MADX convention used only for input
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   //  | dx/dx,   dx/dthx    |
1114   //  | dthx/dx, dtchx/dthx |
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 //real angles in the matrix, MADX convention used only for input
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   //  | dy/dy,   dy/dthy    |
1164   //  | dthy/dy, dtchy/dthy |
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 //MADX convention used only for input
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 //MADX convention used only for input
1204 //angular dispersion
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!