CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
LHCOpticsApproximator.cc
Go to the documentation of this file.
1 #include "SimTransport/TotemRPProtonTransportParametrization/interface/LHCOpticsApproximator.h"
2 #include "FWCore/MessageLogger/interface/MessageLogger.h"
3 
4 #include <vector>
5 #include <iostream>
6 #include "TROOT.h"
7 #include "TFile.h"
8 #include <memory>
9 #include "TMatrixD.h"
10 #include "TMath.h"
11 
12 ClassImp(LHCOpticsApproximator);
13 ClassImp(LHCApertureApproximator);
14 
15 void LHCOpticsApproximator::Init() {
16  out_polynomials.clear();
17  apertures_.clear();
18  out_polynomials.push_back(&x_parametrisation);
19  out_polynomials.push_back(&theta_x_parametrisation);
20  out_polynomials.push_back(&y_parametrisation);
21  out_polynomials.push_back(&theta_y_parametrisation);
22 
23  coord_names.clear();
24  coord_names.push_back("x");
25  coord_names.push_back("theta_x");
26  coord_names.push_back("y");
27  coord_names.push_back("theta_y");
28  coord_names.push_back("ksi");
29 
30  s_begin_ = 0.0;
31  s_end_ = 0.0;
32  trained_ = false;
33 }
34 
35 LHCOpticsApproximator::LHCOpticsApproximator(std::string name,
36  std::string title,
37  TMultiDimFet::EMDFPolyType polynom_type,
38  std::string beam_direction,
39  double nominal_beam_momentum)
40  : x_parametrisation(5, polynom_type, "k"),
41  theta_x_parametrisation(5, polynom_type, "k"),
42  y_parametrisation(5, polynom_type, "k"),
43  theta_y_parametrisation(5, polynom_type, "k") {
44  this->SetName(name.c_str());
45  this->SetTitle(title.c_str());
46  Init();
47 
48  if (beam_direction == "lhcb1")
49  beam = lhcb1;
50  else if (beam_direction == "lhcb2")
51  beam = lhcb2;
52  else
53  beam = lhcb1;
54 
55  nominal_beam_momentum_ = nominal_beam_momentum;
56  nominal_beam_energy_ = TMath::Sqrt(nominal_beam_momentum_ * nominal_beam_momentum_ + 0.938272029 * 0.938272029);
57 }
58 
59 LHCOpticsApproximator::LHCOpticsApproximator() {
60  Init();
61  beam = lhcb1;
62  nominal_beam_momentum_ = 7000;
63  nominal_beam_energy_ = TMath::Sqrt(nominal_beam_momentum_ * nominal_beam_momentum_ + 0.938272029 * 0.938272029);
64 }
65 
66 //MADX canonical variables
67 //(x, theta_x, y, theta_y, ksi) [m, rad, m, rad, -1..0]
68 double LHCOpticsApproximator::ParameterOutOfRangePenalty(double in[], bool invert_beam_coord_sytems) const {
69  double in_corrected[5];
70  if (beam == lhcb1 || !invert_beam_coord_sytems) {
71  in_corrected[0] = in[0];
72  in_corrected[1] = in[1];
73  in_corrected[2] = in[2];
74  in_corrected[3] = in[3];
75  in_corrected[4] = in[4];
76  } else {
77  in_corrected[0] = -in[0];
78  in_corrected[1] = -in[1];
79  in_corrected[2] = in[2];
80  in_corrected[3] = in[3];
81  in_corrected[4] = in[4];
82  }
83 
84  const TVectorD *min_var = x_parametrisation.GetMinVariables();
85  const TVectorD *max_var = x_parametrisation.GetMaxVariables();
86  double res = 0.;
87 
88  for (int i = 0; i < 5; i++) {
89  if (in_corrected[i] < (*min_var)(i)) {
90  double dist = TMath::Abs(((*min_var)(i)-in_corrected[i]) / ((*max_var)(i) - (*min_var)(i)));
91  res += 8 * (TMath::Exp(dist) - 1.0);
92  in_corrected[i] = (*min_var)(i);
93  } else if (in_corrected[i] > (*max_var)(i)) {
94  double dist = TMath::Abs((in_corrected[i] - (*max_var)(i)) / ((*max_var)(i) - (*min_var)(i)));
95  res += 8 * (TMath::Exp(dist) - 1.0);
96  in_corrected[i] = (*max_var)(i);
97  }
98  }
99  return res;
100 }
101 
102 bool LHCOpticsApproximator::Transport(const double *in,
103  double *out,
104  bool check_apertures,
105  bool invert_beam_coord_sytems) const {
106  if (in == nullptr || out == nullptr || !trained_)
107  return false;
108 
109  bool res = CheckInputRange(in);
110  double in_corrected[5];
111 
112  if (beam == lhcb1 || !invert_beam_coord_sytems) {
113  in_corrected[0] = in[0];
114  in_corrected[1] = in[1];
115  in_corrected[2] = in[2];
116  in_corrected[3] = in[3];
117  in_corrected[4] = in[4];
118  out[0] = x_parametrisation.Eval(in_corrected);
119  out[1] = theta_x_parametrisation.Eval(in_corrected);
120  out[2] = y_parametrisation.Eval(in_corrected);
121  out[3] = theta_y_parametrisation.Eval(in_corrected);
122  out[4] = in[4];
123  } else {
124  in_corrected[0] = -in[0];
125  in_corrected[1] = -in[1];
126  in_corrected[2] = in[2];
127  in_corrected[3] = in[3];
128  in_corrected[4] = in[4];
129  out[0] = -x_parametrisation.Eval(in_corrected);
130  out[1] = -theta_x_parametrisation.Eval(in_corrected);
131  out[2] = y_parametrisation.Eval(in_corrected);
132  out[3] = theta_y_parametrisation.Eval(in_corrected);
133  out[4] = in[4];
134  }
135 
136  if (check_apertures) {
137  for (unsigned int i = 0; i < apertures_.size(); i++) {
138  res = res && apertures_[i].CheckAperture(in);
139  }
140  }
141  return res;
142 }
143 
144 bool LHCOpticsApproximator::Transport2D(const double *in,
145  double *out,
146  bool check_apertures,
147  bool invert_beam_coord_sytems) const
148 //return true if transport possible, double x, y
149 {
150  if (in == nullptr || out == nullptr || !trained_)
151  return false;
152 
153  bool res = CheckInputRange(in);
154  double in_corrected[5];
155 
156  if (beam == lhcb1 || !invert_beam_coord_sytems) {
157  in_corrected[0] = in[0];
158  in_corrected[1] = in[1];
159  in_corrected[2] = in[2];
160  in_corrected[3] = in[3];
161  in_corrected[4] = in[4];
162  out[0] = x_parametrisation.Eval(in_corrected);
163  out[1] = y_parametrisation.Eval(in_corrected);
164  } else {
165  in_corrected[0] = -in[0];
166  in_corrected[1] = -in[1];
167  in_corrected[2] = in[2];
168  in_corrected[3] = in[3];
169  in_corrected[4] = in[4];
170  out[0] = -x_parametrisation.Eval(in_corrected);
171  out[1] = y_parametrisation.Eval(in_corrected);
172  }
173 
174  if (check_apertures) {
175  for (unsigned int i = 0; i < apertures_.size(); i++) {
176  res = res && apertures_[i].CheckAperture(in);
177  }
178  }
179  return res;
180 }
181 
182 bool LHCOpticsApproximator::Transport_m_GeV(double in_pos[3],
183  double in_momentum[3],
184  double out_pos[3],
185  double out_momentum[3],
186  bool check_apertures,
187  double z2_z1_dist) const {
188  double in[5];
189  double out[5];
190  double part_mom = 0.0;
191  for (int i = 0; i < 3; i++)
192  part_mom += in_momentum[i] * in_momentum[i];
193 
194  part_mom = TMath::Sqrt(part_mom);
195 
196  in[0] = in_pos[0];
197  in[1] = in_momentum[0] / nominal_beam_momentum_;
198  in[2] = in_pos[1];
199  in[3] = in_momentum[1] / nominal_beam_momentum_;
200  in[4] = (part_mom - nominal_beam_momentum_) / nominal_beam_momentum_;
201 
202  bool res = Transport(in, out, check_apertures, true);
203 
204  out_pos[0] = out[0];
205  out_pos[1] = out[2];
206  out_pos[2] = in_pos[2] + z2_z1_dist;
207 
208  out_momentum[0] = out[1] * nominal_beam_momentum_;
209  out_momentum[1] = out[3] * nominal_beam_momentum_;
210  double part_out_total_mom = (out[4] + 1) * nominal_beam_momentum_;
211  out_momentum[2] = TMath::Sqrt(part_out_total_mom * part_out_total_mom - out_momentum[0] * out_momentum[0] -
212  out_momentum[1] * out_momentum[1]);
213  out_momentum[2] = TMath::Sign(out_momentum[2], in_momentum[2]);
214 
215  return res;
216 }
217 
218 bool LHCOpticsApproximator::Transport(const MadKinematicDescriptor *in,
219  MadKinematicDescriptor *out,
220  bool check_apertures,
221  bool invert_beam_coord_sytems) const {
222  if (in == nullptr || out == nullptr || !trained_)
223  return false;
224 
225  Double_t input[5];
226  Double_t output[5];
227  input[0] = in->x;
228  input[1] = in->theta_x;
229  input[2] = in->y;
230  input[3] = in->theta_y;
231  input[4] = in->ksi;
232 
233  //transport inverts the coordinate systems
234  bool res = Transport(input, output, check_apertures, invert_beam_coord_sytems);
235 
236  out->x = output[0];
237  out->theta_x = output[1];
238  out->y = output[2];
239  out->theta_y = output[3];
240  out->ksi = output[4];
241 
242  return res;
243 }
244 
245 LHCOpticsApproximator::LHCOpticsApproximator(const LHCOpticsApproximator &org)
246  : TNamed(org),
247  x_parametrisation(org.x_parametrisation),
248  theta_x_parametrisation(org.theta_x_parametrisation),
249  y_parametrisation(org.y_parametrisation),
250  theta_y_parametrisation(org.theta_y_parametrisation) {
251  Init();
252  s_begin_ = org.s_begin_;
253  s_end_ = org.s_end_;
254  trained_ = org.trained_;
255  apertures_ = org.apertures_;
256  beam = org.beam;
257  nominal_beam_energy_ = org.nominal_beam_energy_;
258  nominal_beam_momentum_ = org.nominal_beam_momentum_;
259 }
260 
261 const LHCOpticsApproximator &LHCOpticsApproximator::operator=(const LHCOpticsApproximator &org) {
262  if (this != &org) {
263  x_parametrisation = org.x_parametrisation;
264  theta_x_parametrisation = org.theta_x_parametrisation;
265  y_parametrisation = org.y_parametrisation;
266  theta_y_parametrisation = org.theta_y_parametrisation;
267  Init();
268 
269  TNamed::operator=(org);
270  s_begin_ = org.s_begin_;
271  s_end_ = org.s_end_;
272  trained_ = org.trained_;
273 
274  apertures_ = org.apertures_;
275  beam = org.beam;
276  nominal_beam_energy_ = org.nominal_beam_energy_;
277  nominal_beam_momentum_ = org.nominal_beam_momentum_;
278  }
279  return org;
280 }
281 
282 void LHCOpticsApproximator::Train(TTree *inp_tree,
283  std::string data_prefix,
284  polynomials_selection mode,
285  int max_degree_x,
286  int max_degree_tx,
287  int max_degree_y,
288  int max_degree_ty,
289  bool common_terms,
290  double *prec) {
291  if (inp_tree == nullptr)
292  return;
293 
294  InitializeApproximators(mode, max_degree_x, max_degree_tx, max_degree_y, max_degree_ty, common_terms);
295 
296  //in-variables
297  //x_in, theta_x_in, y_in, theta_y_in, ksi_in, s_in
298  double in_var[6];
299 
300  //out-variables
301  //x_out, theta_x_out, y_out, theta_y_out, ksi_out, s_out, valid_out;
302  double out_var[7];
303 
304  //in- out-lables
305  std::string x_in_lab = "x_in";
306  std::string theta_x_in_lab = "theta_x_in";
307  std::string y_in_lab = "y_in";
308  std::string theta_y_in_lab = "theta_y_in";
309  std::string ksi_in_lab = "ksi_in";
310  std::string s_in_lab = "s_in";
311 
312  std::string x_out_lab = data_prefix + "_x_out";
313  std::string theta_x_out_lab = data_prefix + "_theta_x_out";
314  std::string y_out_lab = data_prefix + "_y_out";
315  std::string theta_y_out_lab = data_prefix + "_theta_y_out";
316  std::string ksi_out_lab = data_prefix + "_ksi_out";
317  std::string s_out_lab = data_prefix + "_s_out";
318  std::string valid_out_lab = data_prefix + "_valid_out";
319 
320  //disable not needed branches to speed up the readin
321  inp_tree->SetBranchStatus("*", false); //disable all branches
322  inp_tree->SetBranchStatus(x_in_lab.c_str(), true);
323  inp_tree->SetBranchStatus(theta_x_in_lab.c_str(), true);
324  inp_tree->SetBranchStatus(y_in_lab.c_str(), true);
325  inp_tree->SetBranchStatus(theta_y_in_lab.c_str(), true);
326  inp_tree->SetBranchStatus(ksi_in_lab.c_str(), true);
327  inp_tree->SetBranchStatus(x_out_lab.c_str(), true);
328  inp_tree->SetBranchStatus(theta_x_out_lab.c_str(), true);
329  inp_tree->SetBranchStatus(y_out_lab.c_str(), true);
330  inp_tree->SetBranchStatus(theta_y_out_lab.c_str(), true);
331  inp_tree->SetBranchStatus(ksi_out_lab.c_str(), true);
332  inp_tree->SetBranchStatus(valid_out_lab.c_str(), true);
333 
334  //set input data adresses
335  inp_tree->SetBranchAddress(x_in_lab.c_str(), &(in_var[0]));
336  inp_tree->SetBranchAddress(theta_x_in_lab.c_str(), &(in_var[1]));
337  inp_tree->SetBranchAddress(y_in_lab.c_str(), &(in_var[2]));
338  inp_tree->SetBranchAddress(theta_y_in_lab.c_str(), &(in_var[3]));
339  inp_tree->SetBranchAddress(ksi_in_lab.c_str(), &(in_var[4]));
340  inp_tree->SetBranchAddress(s_in_lab.c_str(), &(in_var[5]));
341 
342  //set output data adresses
343  inp_tree->SetBranchAddress(x_out_lab.c_str(), &(out_var[0]));
344  inp_tree->SetBranchAddress(theta_x_out_lab.c_str(), &(out_var[1]));
345  inp_tree->SetBranchAddress(y_out_lab.c_str(), &(out_var[2]));
346  inp_tree->SetBranchAddress(theta_y_out_lab.c_str(), &(out_var[3]));
347  inp_tree->SetBranchAddress(ksi_out_lab.c_str(), &(out_var[4]));
348  inp_tree->SetBranchAddress(s_out_lab.c_str(), &(out_var[5]));
349  inp_tree->SetBranchAddress(valid_out_lab.c_str(), &(out_var[6]));
350 
351  Long64_t entries = inp_tree->GetEntries();
352  if (entries > 0) {
353  inp_tree->SetBranchStatus(s_in_lab.c_str(), true);
354  inp_tree->SetBranchStatus(s_out_lab.c_str(), true);
355  inp_tree->GetEntry(0);
356  s_begin_ = in_var[5];
357  s_end_ = out_var[5];
358  inp_tree->SetBranchStatus(s_in_lab.c_str(), false);
359  inp_tree->SetBranchStatus(s_out_lab.c_str(), false);
360  }
361 
362  //set input and output variables for fitting
363  for (Long64_t i = 0; i < entries; ++i) {
364  inp_tree->GetEntry(i);
365  if (out_var[6] != 0) //if out data valid
366  {
367  x_parametrisation.AddRow(in_var, out_var[0], 0);
368  theta_x_parametrisation.AddRow(in_var, out_var[1], 0);
369  y_parametrisation.AddRow(in_var, out_var[2], 0);
370  theta_y_parametrisation.AddRow(in_var, out_var[3], 0);
371  }
372  }
373 
374  edm::LogInfo("LHCOpticsApproximator") << "Optical functions parametrizations from " << s_begin_ << " to " << s_end_
375  << "\n";
376  PrintInputRange();
377  for (int i = 0; i < 4; i++) {
378  double best_precision = 0.0;
379  if (prec)
380  best_precision = prec[i];
381  out_polynomials[i]->FindParameterization(best_precision);
382  edm::LogInfo("LHCOpticsApproximator") << "Out variable " << coord_names[i] << " polynomial"
383  << "\n";
384  out_polynomials[i]->PrintPolynomialsSpecial("M");
385  edm::LogInfo("LHCOpticsApproximator") << "\n";
386  }
387 
388  trained_ = true;
389 }
390 
391 void LHCOpticsApproximator::InitializeApproximators(polynomials_selection mode,
392  int max_degree_x,
393  int max_degree_tx,
394  int max_degree_y,
395  int max_degree_ty,
396  bool common_terms) {
397  SetDefaultAproximatorSettings(x_parametrisation, VariableType::X, max_degree_x);
398  SetDefaultAproximatorSettings(theta_x_parametrisation, VariableType::THETA_X, max_degree_tx);
399  SetDefaultAproximatorSettings(y_parametrisation, VariableType::Y, max_degree_y);
400  SetDefaultAproximatorSettings(theta_y_parametrisation, VariableType::THETA_Y, max_degree_ty);
401 
402  if (mode == PREDEFINED) {
403  SetTermsManually(x_parametrisation, VariableType::X, max_degree_x, common_terms);
404  SetTermsManually(theta_x_parametrisation, VariableType::THETA_X, max_degree_tx, common_terms);
405  SetTermsManually(y_parametrisation, VariableType::Y, max_degree_y, common_terms);
406  SetTermsManually(theta_y_parametrisation, VariableType::THETA_Y, max_degree_ty, common_terms);
407  }
408 }
409 
410 void LHCOpticsApproximator::SetDefaultAproximatorSettings(TMultiDimFet &approximator,
411  VariableType var_type,
412  int max_degree) {
413  if (max_degree < 1 || max_degree > 20)
414  max_degree = 10;
415 
416  if (var_type == VariableType::X || var_type == VariableType::THETA_X) {
417  Int_t mPowers[] = {1, 1, 0, 0, max_degree};
418  approximator.SetMaxPowers(mPowers);
419  approximator.SetMaxFunctions(900);
420  approximator.SetMaxStudy(1000);
421  approximator.SetMaxTerms(900);
422  approximator.SetPowerLimit(1.50);
423  approximator.SetMinAngle(10);
424  approximator.SetMaxAngle(10);
425  approximator.SetMinRelativeError(1e-13);
426  }
427 
428  if (var_type == VariableType::Y || var_type == VariableType::THETA_Y) {
429  Int_t mPowers[] = {0, 0, 1, 1, max_degree};
430  approximator.SetMaxPowers(mPowers);
431  approximator.SetMaxFunctions(900);
432  approximator.SetMaxStudy(1000);
433  approximator.SetMaxTerms(900);
434  approximator.SetPowerLimit(1.50);
435  approximator.SetMinAngle(10);
436  approximator.SetMaxAngle(10);
437  approximator.SetMinRelativeError(1e-13);
438  }
439 }
440 
441 void LHCOpticsApproximator::SetTermsManually(TMultiDimFet &approximator,
442  VariableType variable,
443  int max_degree,
444  bool common_terms) {
445  if (max_degree < 1 || max_degree > 20)
446  max_degree = 10;
447 
448  //put terms of shape:
449  //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
450  //t: 0,1,...,max_degree
451 
452  std::vector<Int_t> term_literals;
453  term_literals.reserve(5000);
454 
455  if (variable == VariableType::X || variable == VariableType::THETA_X) {
456  //1,0,0,0,t
457  for (int i = 0; i <= max_degree; ++i) {
458  term_literals.push_back(1);
459  term_literals.push_back(0);
460  term_literals.push_back(0);
461  term_literals.push_back(0);
462  term_literals.push_back(i);
463  }
464  //0,1,0,0,t
465  for (int i = 0; i <= max_degree; ++i) {
466  term_literals.push_back(0);
467  term_literals.push_back(1);
468  term_literals.push_back(0);
469  term_literals.push_back(0);
470  term_literals.push_back(i);
471  }
472  //0,2,0,0,t
473  for (int i = 0; i <= max_degree; ++i) {
474  term_literals.push_back(0);
475  term_literals.push_back(2);
476  term_literals.push_back(0);
477  term_literals.push_back(0);
478  term_literals.push_back(i);
479  }
480  //0,3,0,0,t
481  for (int i = 0; i <= max_degree; ++i) {
482  term_literals.push_back(0);
483  term_literals.push_back(3);
484  term_literals.push_back(0);
485  term_literals.push_back(0);
486  term_literals.push_back(i);
487  }
488  //0,0,0,0,t
489  for (int i = 0; i <= max_degree; ++i) {
490  term_literals.push_back(0);
491  term_literals.push_back(0);
492  term_literals.push_back(0);
493  term_literals.push_back(0);
494  term_literals.push_back(i);
495  }
496  }
497 
498  if (variable == VariableType::Y || variable == VariableType::THETA_Y) {
499  //0,0,1,0,t
500  for (int i = 0; i <= max_degree; ++i) {
501  term_literals.push_back(0);
502  term_literals.push_back(0);
503  term_literals.push_back(1);
504  term_literals.push_back(0);
505  term_literals.push_back(i);
506  }
507  //0,0,0,1,t
508  for (int i = 0; i <= max_degree; ++i) {
509  term_literals.push_back(0);
510  term_literals.push_back(0);
511  term_literals.push_back(0);
512  term_literals.push_back(1);
513  term_literals.push_back(i);
514  }
515  //0,0,0,2,t
516  for (int i = 0; i <= max_degree; ++i) {
517  term_literals.push_back(0);
518  term_literals.push_back(0);
519  term_literals.push_back(0);
520  term_literals.push_back(2);
521  term_literals.push_back(i);
522  }
523  //0,0,0,3,t
524  for (int i = 0; i <= max_degree; ++i) {
525  term_literals.push_back(0);
526  term_literals.push_back(0);
527  term_literals.push_back(0);
528  term_literals.push_back(3);
529  term_literals.push_back(i);
530  }
531  //0,0,0,0,t
532  for (int i = 0; i <= max_degree; ++i) {
533  term_literals.push_back(0);
534  term_literals.push_back(0);
535  term_literals.push_back(0);
536  term_literals.push_back(0);
537  term_literals.push_back(i);
538  }
539  }
540 
541  //push common terms
542  if (common_terms) {
543  term_literals.push_back(1), term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(0),
544  term_literals.push_back(0);
545  term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(0),
546  term_literals.push_back(0);
547  term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(1),
548  term_literals.push_back(0);
549  term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(1), term_literals.push_back(0),
550  term_literals.push_back(0);
551  term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(1),
552  term_literals.push_back(0);
553  term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(1),
554  term_literals.push_back(0);
555 
556  term_literals.push_back(1), term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(0),
557  term_literals.push_back(1);
558  term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(0),
559  term_literals.push_back(1);
560  term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(1),
561  term_literals.push_back(1);
562  term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(1), term_literals.push_back(0),
563  term_literals.push_back(1);
564  term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(1),
565  term_literals.push_back(1);
566  term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(1),
567  term_literals.push_back(1);
568 
569  term_literals.push_back(1), term_literals.push_back(2), term_literals.push_back(0), term_literals.push_back(0),
570  term_literals.push_back(0);
571  term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(2), term_literals.push_back(0),
572  term_literals.push_back(0);
573  term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(2),
574  term_literals.push_back(0);
575  term_literals.push_back(2), term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(0),
576  term_literals.push_back(0);
577  term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(2), term_literals.push_back(0),
578  term_literals.push_back(0);
579  term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(2),
580  term_literals.push_back(0);
581  term_literals.push_back(2), term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(0),
582  term_literals.push_back(0);
583  term_literals.push_back(0), term_literals.push_back(2), term_literals.push_back(1), term_literals.push_back(0),
584  term_literals.push_back(0);
585  term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(2),
586  term_literals.push_back(0);
587  term_literals.push_back(2), term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(1),
588  term_literals.push_back(0);
589  term_literals.push_back(0), term_literals.push_back(2), term_literals.push_back(0), term_literals.push_back(1),
590  term_literals.push_back(0);
591  term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(2), term_literals.push_back(1),
592  term_literals.push_back(0);
593 
594  term_literals.push_back(1), term_literals.push_back(2), term_literals.push_back(0), term_literals.push_back(0),
595  term_literals.push_back(1);
596  term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(2), term_literals.push_back(0),
597  term_literals.push_back(1);
598  term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(2),
599  term_literals.push_back(1);
600  term_literals.push_back(2), term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(0),
601  term_literals.push_back(1);
602  term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(2), term_literals.push_back(0),
603  term_literals.push_back(1);
604  term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(0), term_literals.push_back(2),
605  term_literals.push_back(1);
606  term_literals.push_back(2), term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(0),
607  term_literals.push_back(1);
608  term_literals.push_back(0), term_literals.push_back(2), term_literals.push_back(1), term_literals.push_back(0),
609  term_literals.push_back(1);
610  term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(1), term_literals.push_back(2),
611  term_literals.push_back(1);
612  term_literals.push_back(2), term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(1),
613  term_literals.push_back(1);
614  term_literals.push_back(0), term_literals.push_back(2), term_literals.push_back(0), term_literals.push_back(1),
615  term_literals.push_back(1);
616  term_literals.push_back(0), term_literals.push_back(0), term_literals.push_back(2), term_literals.push_back(1),
617  term_literals.push_back(1);
618  }
619 
620  std::vector<Int_t> powers;
621  powers.resize(term_literals.size());
622 
623  for (unsigned int i = 0; i < term_literals.size(); ++i) {
624  powers[i] = term_literals[i];
625  }
626  approximator.SetPowers(&powers[0], term_literals.size() / 5);
627 }
628 
629 void LHCOpticsApproximator::Test(TTree *inp_tree, TFile *f_out, std::string data_prefix, std::string base_out_dir) {
630  if (inp_tree == nullptr || f_out == nullptr)
631  return;
632 
633  //in-variables
634  //x_in, theta_x_in, y_in, theta_y_in, ksi_in, s_in
635  double in_var[6];
636 
637  //out-variables
638  //x_out, theta_x_out, y_out, theta_y_out, ksi_out, s_out, valid_out;
639  double out_var[7];
640 
641  //in- out-lables
642  std::string x_in_lab = "x_in";
643  std::string theta_x_in_lab = "theta_x_in";
644  std::string y_in_lab = "y_in";
645  std::string theta_y_in_lab = "theta_y_in";
646  std::string ksi_in_lab = "ksi_in";
647  std::string s_in_lab = "s_in";
648 
649  std::string x_out_lab = data_prefix + "_x_out";
650  std::string theta_x_out_lab = data_prefix + "_theta_x_out";
651  std::string y_out_lab = data_prefix + "_y_out";
652  std::string theta_y_out_lab = data_prefix + "_theta_y_out";
653  std::string ksi_out_lab = data_prefix + "_ksi_out";
654  std::string s_out_lab = data_prefix + "_s_out";
655  std::string valid_out_lab = data_prefix + "_valid_out";
656 
657  //disable not needed branches to speed up the readin
658  inp_tree->SetBranchStatus("*", false); //disable all branches
659  inp_tree->SetBranchStatus(x_in_lab.c_str(), true);
660  inp_tree->SetBranchStatus(theta_x_in_lab.c_str(), true);
661  inp_tree->SetBranchStatus(y_in_lab.c_str(), true);
662  inp_tree->SetBranchStatus(theta_y_in_lab.c_str(), true);
663  inp_tree->SetBranchStatus(ksi_in_lab.c_str(), true);
664  inp_tree->SetBranchStatus(x_out_lab.c_str(), true);
665  inp_tree->SetBranchStatus(theta_x_out_lab.c_str(), true);
666  inp_tree->SetBranchStatus(y_out_lab.c_str(), true);
667  inp_tree->SetBranchStatus(theta_y_out_lab.c_str(), true);
668  inp_tree->SetBranchStatus(ksi_out_lab.c_str(), true);
669  inp_tree->SetBranchStatus(valid_out_lab.c_str(), true);
670 
671  //set input data adresses
672  inp_tree->SetBranchAddress(x_in_lab.c_str(), &(in_var[0]));
673  inp_tree->SetBranchAddress(theta_x_in_lab.c_str(), &(in_var[1]));
674  inp_tree->SetBranchAddress(y_in_lab.c_str(), &(in_var[2]));
675  inp_tree->SetBranchAddress(theta_y_in_lab.c_str(), &(in_var[3]));
676  inp_tree->SetBranchAddress(ksi_in_lab.c_str(), &(in_var[4]));
677  inp_tree->SetBranchAddress(s_in_lab.c_str(), &(in_var[5]));
678 
679  //set output data adresses
680  inp_tree->SetBranchAddress(x_out_lab.c_str(), &(out_var[0]));
681  inp_tree->SetBranchAddress(theta_x_out_lab.c_str(), &(out_var[1]));
682  inp_tree->SetBranchAddress(y_out_lab.c_str(), &(out_var[2]));
683  inp_tree->SetBranchAddress(theta_y_out_lab.c_str(), &(out_var[3]));
684  inp_tree->SetBranchAddress(ksi_out_lab.c_str(), &(out_var[4]));
685  inp_tree->SetBranchAddress(s_out_lab.c_str(), &(out_var[5]));
686  inp_tree->SetBranchAddress(valid_out_lab.c_str(), &(out_var[6]));
687 
688  //test histogramms
689  TH1D *err_hists[4];
690  TH2D *err_inp_cor_hists[4][5];
691  TH2D *err_out_cor_hists[4][5];
692 
693  AllocateErrorHists(err_hists);
694  AllocateErrorInputCorHists(err_inp_cor_hists);
695  AllocateErrorOutputCorHists(err_out_cor_hists);
696 
697  Long64_t entries = inp_tree->GetEntries();
698  //set input and output variables for fitting
699  for (Long64_t i = 0; i < entries; ++i) {
700  double errors[4];
701  inp_tree->GetEntry(i);
702 
703  errors[0] = out_var[0] - x_parametrisation.Eval(in_var);
704  errors[1] = out_var[1] - theta_x_parametrisation.Eval(in_var);
705  errors[2] = out_var[2] - y_parametrisation.Eval(in_var);
706  errors[3] = out_var[3] - theta_y_parametrisation.Eval(in_var);
707 
708  FillErrorHistograms(errors, err_hists);
709  FillErrorDataCorHistograms(errors, in_var, err_inp_cor_hists);
710  FillErrorDataCorHistograms(errors, out_var, err_out_cor_hists);
711  }
712 
713  WriteHistograms(err_hists, err_inp_cor_hists, err_out_cor_hists, f_out, base_out_dir);
714 
715  DeleteErrorHists(err_hists);
716  DeleteErrorCorHistograms(err_inp_cor_hists);
717  DeleteErrorCorHistograms(err_out_cor_hists);
718 }
719 
720 void LHCOpticsApproximator::AllocateErrorHists(TH1D *err_hists[4]) {
721  std::vector<std::string> error_labels;
722  error_labels.push_back("x error");
723  error_labels.push_back("theta_x error");
724  error_labels.push_back("y error");
725  error_labels.push_back("theta_y error");
726 
727  for (int i = 0; i < 4; ++i) {
728  err_hists[i] = new TH1D(error_labels[i].c_str(), error_labels[i].c_str(), 100, -0.0000000001, 0.0000000001);
729  err_hists[i]->SetXTitle(error_labels[i].c_str());
730  err_hists[i]->SetYTitle("counts");
731  err_hists[i]->SetDirectory(nullptr);
732  err_hists[i]->SetCanExtend(TH1::kAllAxes);
733  }
734 }
735 
736 void LHCOpticsApproximator::TestAperture(
737  TTree *inp_tree, TTree *out_tree) //x, theta_x, y, theta_y, ksi, mad_accepted, parametriz_accepted
738 {
739  if (inp_tree == nullptr || out_tree == nullptr)
740  return;
741 
742  Long64_t entries = inp_tree->GetEntries();
743  double entry[7];
744  double parametrization_out[5];
745 
746  inp_tree->SetBranchAddress("x", &(entry[0]));
747  inp_tree->SetBranchAddress("theta_x", &(entry[1]));
748  inp_tree->SetBranchAddress("y", &(entry[2]));
749  inp_tree->SetBranchAddress("theta_y", &(entry[3]));
750  inp_tree->SetBranchAddress("ksi", &(entry[4]));
751  inp_tree->SetBranchAddress("mad_accept", &(entry[5]));
752  inp_tree->SetBranchAddress("par_accept", &(entry[6]));
753 
754  out_tree->SetBranchAddress("x", &(entry[0]));
755  out_tree->SetBranchAddress("theta_x", &(entry[1]));
756  out_tree->SetBranchAddress("y", &(entry[2]));
757  out_tree->SetBranchAddress("theta_y", &(entry[3]));
758  out_tree->SetBranchAddress("ksi", &(entry[4]));
759  out_tree->SetBranchAddress("mad_accept", &(entry[5]));
760  out_tree->SetBranchAddress("par_accept", &(entry[6]));
761 
762  // int ind=0;
763  for (Long64_t i = 0; i < entries; i++) {
764  inp_tree->GetEntry(i);
765 
766  //Don't invert the coordinate systems, appertures are defined in the
767  //coordinate system of the beam - perhaps to be changed
768  bool res = Transport(entry, parametrization_out, true, false);
769 
770  if (res)
771  entry[6] = 1.0;
772  else
773  entry[6] = 0.0;
774 
775  out_tree->Fill();
776  }
777 }
778 
779 void LHCOpticsApproximator::AllocateErrorInputCorHists(TH2D *err_inp_cor_hists[4][5]) {
780  std::vector<std::string> error_labels;
781  std::vector<std::string> data_labels;
782 
783  error_labels.push_back("x error");
784  error_labels.push_back("theta_x error");
785  error_labels.push_back("y error");
786  error_labels.push_back("theta_y error");
787 
788  data_labels.push_back("x input");
789  data_labels.push_back("theta_x input");
790  data_labels.push_back("y input");
791  data_labels.push_back("theta_y input");
792  data_labels.push_back("ksi input");
793 
794  for (int eri = 0; eri < 4; ++eri) {
795  for (int dati = 0; dati < 5; ++dati) {
796  std::string name = error_labels[eri] + " vs. " + data_labels[dati];
797  const std::string &title = name;
798  err_inp_cor_hists[eri][dati] =
799  new TH2D(name.c_str(), title.c_str(), 100, -0.0000000001, 0.0000000001, 100, -0.0000000001, 0.0000000001);
800  err_inp_cor_hists[eri][dati]->SetXTitle(error_labels[eri].c_str());
801  err_inp_cor_hists[eri][dati]->SetYTitle(data_labels[dati].c_str());
802  err_inp_cor_hists[eri][dati]->SetDirectory(nullptr);
803  err_inp_cor_hists[eri][dati]->SetCanExtend(TH1::kAllAxes);
804  }
805  }
806 }
807 
808 void LHCOpticsApproximator::AllocateErrorOutputCorHists(TH2D *err_out_cor_hists[4][5]) {
809  std::vector<std::string> error_labels;
810  std::vector<std::string> data_labels;
811 
812  error_labels.push_back("x error");
813  error_labels.push_back("theta_x error");
814  error_labels.push_back("y error");
815  error_labels.push_back("theta_y error");
816 
817  data_labels.push_back("x output");
818  data_labels.push_back("theta_x output");
819  data_labels.push_back("y output");
820  data_labels.push_back("theta_y output");
821  data_labels.push_back("ksi output");
822 
823  for (int eri = 0; eri < 4; ++eri) {
824  for (int dati = 0; dati < 5; ++dati) {
825  std::string name = error_labels[eri] + " vs. " + data_labels[dati];
826  const std::string &title = name;
827  err_out_cor_hists[eri][dati] =
828  new TH2D(name.c_str(), title.c_str(), 100, -0.0000000001, 0.0000000001, 100, -0.0000000001, 0.0000000001);
829  err_out_cor_hists[eri][dati]->SetXTitle(error_labels[eri].c_str());
830  err_out_cor_hists[eri][dati]->SetYTitle(data_labels[dati].c_str());
831  err_out_cor_hists[eri][dati]->SetDirectory(nullptr);
832  err_out_cor_hists[eri][dati]->SetCanExtend(TH1::kAllAxes);
833  }
834  }
835 }
836 
837 void LHCOpticsApproximator::FillErrorHistograms(double errors[4], TH1D *err_hists[4]) {
838  for (int i = 0; i < 4; ++i) {
839  err_hists[i]->Fill(errors[i]);
840  }
841 }
842 
843 void LHCOpticsApproximator::FillErrorDataCorHistograms(double errors[4], double var[5], TH2D *err_cor_hists[4][5]) {
844  for (int eri = 0; eri < 4; ++eri) {
845  for (int dati = 0; dati < 5; ++dati) {
846  err_cor_hists[eri][dati]->Fill(errors[eri], var[dati]);
847  }
848  }
849 }
850 
851 void LHCOpticsApproximator::DeleteErrorHists(TH1D *err_hists[4]) {
852  for (int i = 0; i < 4; ++i) {
853  delete err_hists[i];
854  }
855 }
856 
857 void LHCOpticsApproximator::DeleteErrorCorHistograms(TH2D *err_cor_hists[4][5]) {
858  for (int eri = 0; eri < 4; ++eri) {
859  for (int dati = 0; dati < 5; ++dati) {
860  delete err_cor_hists[eri][dati];
861  }
862  }
863 }
864 
865 void LHCOpticsApproximator::WriteHistograms(TH1D *err_hists[4],
866  TH2D *err_inp_cor_hists[4][5],
867  TH2D *err_out_cor_hists[4][5],
868  TFile *f_out,
869  std::string base_out_dir) {
870  if (f_out == nullptr)
871  return;
872 
873  f_out->cd();
874  if (!gDirectory->cd(base_out_dir.c_str()))
875  gDirectory->mkdir(base_out_dir.c_str());
876 
877  gDirectory->cd(base_out_dir.c_str());
878  gDirectory->mkdir(this->GetName());
879  gDirectory->cd(this->GetName());
880  gDirectory->mkdir("x");
881  gDirectory->mkdir("theta_x");
882  gDirectory->mkdir("y");
883  gDirectory->mkdir("theta_y");
884 
885  gDirectory->cd("x");
886  err_hists[0]->Write("", TObject::kWriteDelete);
887  for (int i = 0; i < 5; i++) {
888  err_inp_cor_hists[0][i]->Write("", TObject::kWriteDelete);
889  err_out_cor_hists[0][i]->Write("", TObject::kWriteDelete);
890  }
891 
892  gDirectory->cd("..");
893  gDirectory->cd("theta_x");
894  err_hists[1]->Write("", TObject::kWriteDelete);
895  for (int i = 0; i < 5; i++) {
896  err_inp_cor_hists[1][i]->Write("", TObject::kWriteDelete);
897  err_out_cor_hists[1][i]->Write("", TObject::kWriteDelete);
898  }
899 
900  gDirectory->cd("..");
901  gDirectory->cd("y");
902  err_hists[2]->Write("", TObject::kWriteDelete);
903  for (int i = 0; i < 5; i++) {
904  err_inp_cor_hists[2][i]->Write("", TObject::kWriteDelete);
905  err_out_cor_hists[2][i]->Write("", TObject::kWriteDelete);
906  }
907 
908  gDirectory->cd("..");
909  gDirectory->cd("theta_y");
910  err_hists[3]->Write("", TObject::kWriteDelete);
911  for (int i = 0; i < 5; i++) {
912  err_inp_cor_hists[3][i]->Write("", TObject::kWriteDelete);
913  err_out_cor_hists[3][i]->Write("", TObject::kWriteDelete);
914  }
915  gDirectory->cd("..");
916  gDirectory->cd("..");
917 }
918 
919 void LHCOpticsApproximator::PrintInputRange() {
920  const TVectorD *min_var = x_parametrisation.GetMinVariables();
921  const TVectorD *max_var = x_parametrisation.GetMaxVariables();
922 
923  edm::LogInfo("LHCOpticsApproximator") << "Covered input parameters range:"
924  << "\n";
925  for (int i = 0; i < 5; i++) {
926  edm::LogInfo("LHCOpticsApproximator") << (*min_var)(i) << " < " << coord_names[i] << " < " << (*max_var)(i) << "\n";
927  }
928  edm::LogInfo("LHCOpticsApproximator") << "\n";
929 }
930 
931 bool LHCOpticsApproximator::CheckInputRange(const double *in,
932  bool invert_beam_coord_sytems) const //x, thx, y, thy, ksi
933 {
934  double in_corrected[5];
935  if (beam == lhcb1 || !invert_beam_coord_sytems) {
936  in_corrected[0] = in[0];
937  in_corrected[1] = in[1];
938  in_corrected[2] = in[2];
939  in_corrected[3] = in[3];
940  in_corrected[4] = in[4];
941  } else {
942  in_corrected[0] = -in[0];
943  in_corrected[1] = -in[1];
944  in_corrected[2] = in[2];
945  in_corrected[3] = in[3];
946  in_corrected[4] = in[4];
947  }
948 
949  const TVectorD *min_var = x_parametrisation.GetMinVariables();
950  const TVectorD *max_var = x_parametrisation.GetMaxVariables();
951  bool res = true;
952 
953  for (int i = 0; i < 5; i++) {
954  res = res && in_corrected[i] >= (*min_var)(i) && in_corrected[i] <= (*max_var)(i);
955  }
956  return res;
957 }
958 
959 void LHCOpticsApproximator::AddRectEllipseAperture(
960  const LHCOpticsApproximator &in, double rect_x, double rect_y, double r_el_x, double r_el_y) {
961  apertures_.push_back(
962  LHCApertureApproximator(in, rect_x, rect_y, r_el_x, r_el_y, LHCApertureApproximator::ApertureType::RECTELLIPSE));
963 }
964 
965 LHCApertureApproximator::LHCApertureApproximator() {
966  rect_x_ = rect_y_ = r_el_x_ = r_el_y_ = 0.0;
967  ap_type_ = ApertureType::NO_APERTURE;
968 }
969 
970 LHCApertureApproximator::LHCApertureApproximator(
971  const LHCOpticsApproximator &in, double rect_x, double rect_y, double r_el_x, double r_el_y, ApertureType type)
972  : LHCOpticsApproximator(in) {
973  rect_x_ = rect_x;
974  rect_y_ = rect_y;
975  r_el_x_ = r_el_x;
976  r_el_y_ = r_el_y;
977  ap_type_ = type;
978 }
979 
980 bool LHCApertureApproximator::CheckAperture(const double *in,
981  bool invert_beam_coord_sytems) const //x, thx. y, thy, ksi
982 {
983  double out[5];
984  bool result = Transport(in, out, false, invert_beam_coord_sytems);
985 
986  if (ap_type_ == ApertureType::RECTELLIPSE) {
987  result = result && out[0] < rect_x_ && out[0] > -rect_x_ && out[2] < rect_y_ && out[2] > -rect_y_ &&
988  (out[0] * out[0] / (r_el_x_ * r_el_x_) + out[2] * out[2] / (r_el_y_ * r_el_y_) < 1);
989  }
990  return result;
991 }
992 
993 void LHCOpticsApproximator::PrintOpticalFunctions() {
994  edm::LogInfo("LHCOpticsApproximator") << std::endl
995  << "Linear terms of optical functions:"
996  << "\n";
997  for (int i = 0; i < 4; i++) {
998  PrintCoordinateOpticalFunctions(*out_polynomials[i], coord_names[i], coord_names);
999  }
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 }
LHCOpticsApproximator::InitializeApproximators
void InitializeApproximators(polynomials_selection mode, int max_degree_x, int max_degree_tx, int max_degree_y, int max_degree_ty, bool common_terms)
Definition: LHCOpticsApproximator.cc:391
LHCOpticsApproximator::FillErrorDataCorHistograms
void FillErrorDataCorHistograms(double errors[4], double var[5], TH2D *err_cor_hists[4][5])
Definition: LHCOpticsApproximator.cc:843
LHCOpticsApproximator::DeleteErrorHists
void DeleteErrorHists(TH1D *err_hists[4])
Definition: LHCOpticsApproximator.cc:851
LHCOpticsApproximator::TestAperture
void TestAperture(TTree *in_tree, TTree *out_tree)
x, theta_x, y, theta_y, ksi, mad_accepted, parametriz_accepted
Definition: LHCOpticsApproximator.cc:736
mps_fire.i
i
Definition: mps_fire.py:428
LHCOpticsApproximator::Transport2D
bool Transport2D(const double *in, double *out, bool check_apertures=false, bool invert_beam_coord_sytems=true) const
Definition: LHCOpticsApproximator.cc:144
TMultiDimFet::SetPowerLimit
void SetPowerLimit(Double_t limit=1e-3)
Definition: TMultiDimFet.cc:2047
LHCOpticsApproximator::Transport
bool Transport(const double *in, double *out, bool check_apertures=false, bool invert_beam_coord_sytems=true) const
Definition: LHCOpticsApproximator.cc:102
LHCOpticsApproximator::apertures_
std::vector< LHCApertureApproximator > apertures_
apertures on the way
Definition: LHCOpticsApproximator.h:170
Sign
T Sign(T A, T B)
Definition: MathUtil.h:54
LHCOpticsApproximator::theta_x_parametrisation
TMultiDimFet theta_x_parametrisation
polynomial approximation for theta_x
Definition: LHCOpticsApproximator.h:177
TMultiDimFet::SetMaxPowers
void SetMaxPowers(const Int_t *powers)
Definition: TMultiDimFet.cc:2056
TMultiDimFet::GetMinVariables
const TVectorD * GetMinVariables() const
Definition: TMultiDimFet.h:166
TMultiDimFet::SetMaxTerms
void SetMaxTerms(Int_t terms)
Definition: TMultiDimFet.h:206
LHCOpticsApproximator::nominal_beam_momentum_
double nominal_beam_momentum_
GeV/c.
Definition: LHCOpticsApproximator.h:166
LHCOpticsApproximator::AllocateErrorOutputCorHists
void AllocateErrorOutputCorHists(TH2D *err_out_cor_hists[4][5])
Definition: LHCOpticsApproximator.cc:808
LHCOpticsApproximator::Transport_m_GeV
bool Transport_m_GeV(double in_pos[3], double in_momentum[3], double out_pos[3], double out_momentum[3], bool check_apertures, double z2_z1_dist) const
pos, momentum: x,y,z; pos in m, momentum in GeV/c
Definition: LHCOpticsApproximator.cc:182
LHCOpticsApproximator::CheckInputRange
bool CheckInputRange(const double *in, bool invert_beam_coord_sytems=true) const
Definition: LHCOpticsApproximator.cc:931
LHCOpticsApproximator::SetTermsManually
void SetTermsManually(TMultiDimFet &approximator, VariableType variable, int max_degree, bool common_terms)
Definition: LHCOpticsApproximator.cc:441
LHCOpticsApproximator::GetLineariasedTransportMatrixY
void GetLineariasedTransportMatrixY(double mad_init_x, double mad_init_thx, double mad_init_y, double mad_init_thy, double mad_init_xi, TMatrixD &tr_matrix, double d_mad_y=10e-6, double d_mad_thy=10e-6)
returns linearised transport matrix for y projection | dy_out/dy_in dy_out/dthy_in | | dthy_out/dy_in...
Definition: LHCOpticsApproximator.cc:1109
LHCOpticsApproximator::AddRectEllipseAperture
void AddRectEllipseAperture(const LHCOpticsApproximator &in, double rect_x, double rect_y, double r_el_x, double r_el_y)
Definition: LHCOpticsApproximator.cc:959
LHCOpticsApproximator::ParameterOutOfRangePenalty
double ParameterOutOfRangePenalty(double par_m[], bool invert_beam_coord_sytems=true) const
Definition: LHCOpticsApproximator.cc:68
LHCOpticsApproximator::s_begin_
double s_begin_
begin of transport along the reference orbit
Definition: LHCOpticsApproximator.h:162
TMultiDimFet::SetMinRelativeError
void SetMinRelativeError(Double_t error)
Definition: TMultiDimFet.cc:2068
type
type
Definition: SiPixelVCal_PayloadInspector.cc:39
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
input
static std::string const input
Definition: EdmProvDump.cc:47
mps_fire.result
tuple result
Definition: mps_fire.py:311
TMultiDimFet::SetPowers
virtual void SetPowers(const Int_t *powers, Int_t terms)
Definition: TMultiDimFet.cc:2027
TMultiDimFet::SetMaxStudy
void SetMaxStudy(Int_t n)
Definition: TMultiDimFet.h:205
LHCOpticsApproximator::theta_y_parametrisation
TMultiDimFet theta_y_parametrisation
polynomial approximation for theta_y
Definition: LHCOpticsApproximator.h:179
LHCOpticsApproximator::AllocateErrorHists
void AllocateErrorHists(TH1D *err_hists[4])
Definition: LHCOpticsApproximator.cc:720
EcalCondDBWriter_cfi.beam
string beam
Definition: EcalCondDBWriter_cfi.py:45
LHCOpticsApproximator::GetLinearApproximation
void GetLinearApproximation(double atPoint[], double &Cx, double &Lx, double &vx, double &Cy, double &Ly, double &vy, double &D, double ep=1E-5)
Definition: LHCOpticsApproximator.cc:1033
LHCOpticsApproximator::trained_
bool trained_
trained polynomials
Definition: LHCOpticsApproximator.h:167
LHCOpticsApproximator::PrintCoordinateOpticalFunctions
void PrintCoordinateOpticalFunctions(TMultiDimFet &parametrization, const std::string &coord_name, const std::vector< std::string > &input_vars)
Definition: LHCOpticsApproximator.cc:1002
LHCOpticsApproximator::out_polynomials
std::vector< TMultiDimFet * > out_polynomials
Definition: LHCOpticsApproximator.h:168
isotrackApplyRegressor.var
list var
if using global norm cols_to_minmax = [&#39;t_delta&#39;, &#39;t_hmaxNearP&#39;,&#39;t_emaxNearP&#39;, &#39;t_hAnnular&#39;, &#39;t_eAnnular&#39;,&#39;t_pt&#39;,&#39;t_nVtx&#39;,&#39;t_ieta&#39;,&#39;t_eHcal10&#39;, &#39;t_eHcal30&#39;,&#39;t_rhoh&#39;,&#39;t_eHcal&#39;] df[cols_to_minmax] = df[cols_to_minmax].apply(lambda x: (x - x.min()) / (x.max() - x.min()) if (x.max() - x.min() &gt; 0) else 1.0/200.0)
Definition: isotrackApplyRegressor.py:79
Abs
T Abs(T a)
Definition: MathUtil.h:49
ClassImp
ClassImp(AliDaqEventHeader)
LHCOpticsApproximator::GetDxds
double GetDxds(double mad_init_x, double mad_init_thx, double mad_init_y, double mad_init_thy, double mad_init_xi, double d_mad_xi=0.001)
Definition: LHCOpticsApproximator.cc:1181
LHCOpticsApproximator::Test
void Test(TTree *inp_tree, TFile *f_out, std::string data_prefix=std::string("def"), std::string base_out_dir=std::string(""))
Definition: LHCOpticsApproximator.cc:629
TMultiDimFet::GetMaxVariables
const TVectorD * GetMaxVariables() const
Definition: TMultiDimFet.h:160
LHCOpticsApproximator
Class finds the parametrisation of MADX proton transport and transports the protons according to it 5...
Definition: LHCOpticsApproximator.h:29
TMultiDimFet::AddRow
virtual void AddRow(const Double_t *x, Double_t D, Double_t E=0)
Definition: TMultiDimFet.cc:236
LHCOpticsApproximator::GetLineariasedTransportMatrixX
void GetLineariasedTransportMatrixX(double mad_init_x, double mad_init_thx, double mad_init_y, double mad_init_thy, double mad_init_xi, TMatrixD &tr_matrix, double d_mad_x=10e-6, double d_mad_thx=10e-6)
returns linearised transport matrix for x projection | dx_out/dx_in dx_out/dthx_in | | dthx_out/dx_in...
Definition: LHCOpticsApproximator.cc:1065
LHCOpticsApproximator::DeleteErrorCorHistograms
void DeleteErrorCorHistograms(TH2D *err_cor_hists[4][5])
Definition: LHCOpticsApproximator.cc:857
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:125
TMultiDimFet::SetMaxFunctions
void SetMaxFunctions(Int_t n)
Definition: TMultiDimFet.h:203
convertSQLitetoXML_cfg.output
tuple output
Definition: convertSQLitetoXML_cfg.py:72
LHCOpticsApproximator::SetDefaultAproximatorSettings
void SetDefaultAproximatorSettings(TMultiDimFet &approximator, VariableType var_type, int max_degree)
Definition: LHCOpticsApproximator.cc:410
funct::D
DecomposeProduct< arg, typename Div::arg > D
Definition: Factorize.h:141
LHCOpticsApproximator::s_end_
double s_end_
end of transport along the reference orbit
Definition: LHCOpticsApproximator.h:163
LHCOpticsApproximator::FillErrorHistograms
void FillErrorHistograms(double errors[4], TH1D *err_hists[4])
Definition: LHCOpticsApproximator.cc:837
LHCOpticsApproximator::GetDx
double GetDx(double mad_init_x, double mad_init_thx, double mad_init_y, double mad_init_thy, double mad_init_xi, double d_mad_xi=0.001)
Definition: LHCOpticsApproximator.cc:1153
LHCOpticsApproximator::AllocateErrorInputCorHists
void AllocateErrorInputCorHists(TH2D *err_inp_cor_hists[4][5])
Definition: LHCOpticsApproximator.cc:779
TMultiDimFet::SetMinAngle
void SetMinAngle(Double_t angle=1)
Definition: TMultiDimFet.cc:2013
mps_splice.entry
list entry
Definition: mps_splice.py:68
LHCOpticsApproximator::coord_names
std::vector< std::string > coord_names
pointers to polynomials
Definition: LHCOpticsApproximator.h:169
LHCOpticsApproximator::y_parametrisation
TMultiDimFet y_parametrisation
polynomial approximation for y
Definition: LHCOpticsApproximator.h:178
LHCOpticsApproximator::x_parametrisation
TMultiDimFet x_parametrisation
polynomial approximation for x
Definition: LHCOpticsApproximator.h:176
LHCOpticsApproximator::Train
void Train(TTree *inp_tree, std::string data_prefix=std::string("def"), polynomials_selection mode=PREDEFINED, int max_degree_x=10, int max_degree_tx=10, int max_degree_y=10, int max_degree_ty=10, bool common_terms=false, double *prec=nullptr)
Definition: LHCOpticsApproximator.cc:282
LHCOpticsApproximator::WriteHistograms
void WriteHistograms(TH1D *err_hists[4], TH2D *err_inp_cor_hists[4][5], TH2D *err_out_cor_hists[4][5], TFile *f_out, std::string base_out_dir)
Definition: LHCOpticsApproximator.cc:865
TMultiDimFet::SetMaxAngle
void SetMaxAngle(Double_t angle=0)
Definition: TMultiDimFet.cc:1998
LHCApertureApproximator::CheckAperture
bool CheckAperture(const double *in, bool invert_beam_coord_sytems=true) const
Definition: LHCOpticsApproximator.cc:980
LHCOpticsApproximator::operator=
const LHCOpticsApproximator & operator=(const LHCOpticsApproximator &org)
Definition: LHCOpticsApproximator.cc:261
TMultiDimFet::Eval
virtual Double_t Eval(const Double_t *x, const Double_t *coeff=nullptr) const
Definition: TMultiDimFet.cc:442
Generated for CMSSW Reference Manual by doxygen 1.8.5