Back to home page

Project CMSSW displayed by LXR
 
 

    

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 //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 = 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   //transport inverts the coordinate systems
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   //in-variables
0297   //x_in, theta_x_in, y_in, theta_y_in, ksi_in, s_in
0298   double in_var[6];
0299 
0300   //out-variables
0301   //x_out, theta_x_out, y_out, theta_y_out, ksi_out, s_out, valid_out;
0302   double out_var[7];
0303 
0304   //in- out-lables
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   //disable not needed branches to speed up the readin
0321   inp_tree->SetBranchStatus("*", false);  //disable all branches
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   //set input data adresses
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   //set output data adresses
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   //set input and output variables for fitting
0363   for (Long64_t i = 0; i < entries; ++i) {
0364     inp_tree->GetEntry(i);
0365     if (out_var[6] != 0)  //if out data valid
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   //put terms of shape:
0449   //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
0450   //t: 0,1,...,max_degree
0451 
0452   std::vector<Int_t> term_literals;
0453   term_literals.reserve(5000);
0454 
0455   if (variable == VariableType::X || variable == VariableType::THETA_X) {
0456     //1,0,0,0,t
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     //0,1,0,0,t
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     //0,2,0,0,t
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     //0,3,0,0,t
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     //0,0,0,0,t
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     //0,0,1,0,t
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     //0,0,0,1,t
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     //0,0,0,2,t
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     //0,0,0,3,t
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     //0,0,0,0,t
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   //push common terms
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   //in-variables
0634   //x_in, theta_x_in, y_in, theta_y_in, ksi_in, s_in
0635   double in_var[6];
0636 
0637   //out-variables
0638   //x_out, theta_x_out, y_out, theta_y_out, ksi_out, s_out, valid_out;
0639   double out_var[7];
0640 
0641   //in- out-lables
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   //disable not needed branches to speed up the readin
0658   inp_tree->SetBranchStatus("*", false);  //disable all branches
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   //set input data adresses
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   //set output data adresses
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   //test histogramms
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   //set input and output variables for fitting
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)  //x, theta_x, y, theta_y, ksi, mad_accepted, parametriz_accepted
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   //  int ind=0;
0763   for (Long64_t i = 0; i < entries; i++) {
0764     inp_tree->GetEntry(i);
0765 
0766     //Don't invert the coordinate systems, appertures are defined in the
0767     //coordinate system of the beam - perhaps to be changed
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  //x, thx, y, thy, ksi
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  //x, thx. y, thy, ksi
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 &parametrization,
1003                                                             const std::string &coord_name,
1004                                                             const std::vector<std::string> &input_vars) {
1005   double in[5];
1006   //  double out;
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 //real angles in the matrix, MADX convention used only for input
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   //  | dx/dx,   dx/dthx    |
1100   //  | dthx/dx, dtchx/dthx |
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 //real angles in the matrix, MADX convention used only for input
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   //  | dy/dy,   dy/dthy    |
1144   //  | dthy/dy, dtchy/dthy |
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 //MADX convention used only for input
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 //MADX convention used only for input
1180 //angular dispersion
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 }
[ 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!