Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:31:09

0001 #ifndef TotemRPProtonTransportParametrization_LHC_OPTICS_APPROXIMATOR_H
0002 #define TotemRPProtonTransportParametrization_LHC_OPTICS_APPROXIMATOR_H
0003 
0004 #include <string>
0005 #include "TNamed.h"
0006 #include "TTree.h"
0007 #include "TH1D.h"
0008 #include "TH2D.h"
0009 #include <memory>
0010 #include "TMatrixD.h"
0011 
0012 #include "SimTransport/TotemRPProtonTransportParametrization/interface/TMultiDimFet.h"
0013 
0014 struct MadKinematicDescriptor {
0015   double x;        ///< m
0016   double theta_x;  ///< rad
0017   double y;        ///< m
0018   double theta_y;  ///< rad
0019   double ksi;      ///< 1
0020 };
0021 
0022 class LHCApertureApproximator;
0023 
0024 /**
0025  *\brief Class finds the parametrisation of MADX proton transport and transports the protons according to it
0026  * 5 phase space variables are taken in to configuration: x, y, theta_x, theta_y, xi
0027  * xi < 0 for momentum losses (that is for diffractive protons)
0028 **/
0029 class LHCOpticsApproximator : public TNamed {
0030 public:
0031   LHCOpticsApproximator();
0032   /// begin and end position along the beam of the particle to transport, training_tree, prefix of data branch in the tree
0033   LHCOpticsApproximator(std::string name,
0034                         std::string title,
0035                         TMultiDimFet::EMDFPolyType polynom_type,
0036                         std::string beam_direction,
0037                         double nominal_beam_momentum);
0038   LHCOpticsApproximator(const LHCOpticsApproximator &org);
0039   const LHCOpticsApproximator &operator=(const LHCOpticsApproximator &org);
0040 
0041   enum polynomials_selection { AUTOMATIC, PREDEFINED };
0042   enum beam_type { lhcb1, lhcb2 };
0043   void Train(TTree *inp_tree,
0044              std::string data_prefix = std::string("def"),
0045              polynomials_selection mode = PREDEFINED,
0046              int max_degree_x = 10,
0047              int max_degree_tx = 10,
0048              int max_degree_y = 10,
0049              int max_degree_ty = 10,
0050              bool common_terms = false,
0051              double *prec = nullptr);
0052   void Test(TTree *inp_tree,
0053             TFile *f_out,
0054             std::string data_prefix = std::string("def"),
0055             std::string base_out_dir = std::string(""));
0056   void TestAperture(TTree *in_tree,
0057                     TTree *out_tree);  ///< x, theta_x, y, theta_y, ksi, mad_accepted, parametriz_accepted
0058 
0059   double ParameterOutOfRangePenalty(double par_m[], bool invert_beam_coord_sytems = true) const;
0060 
0061   /// Basic 3D transport method
0062   /// MADX canonical variables
0063   /// IN/OUT: (x, theta_x, y, theta_y, xi) [m, rad, m, rad, 1]
0064   /// returns true if transport possible
0065   /// if theta is calculated from momentum p, use theta_x = p.x() / p.mag() and theta_y = p.y() / p.mag()
0066   bool Transport(const double *in,
0067                  double *out,
0068                  bool check_apertures = false,
0069                  bool invert_beam_coord_sytems = true) const;
0070   bool Transport(const MadKinematicDescriptor *in,
0071                  MadKinematicDescriptor *out,
0072                  bool check_apertures = false,
0073                  bool invert_beam_coord_sytems = true) const;  //return true if transport possible
0074 
0075   /// Basic 2D transport method
0076   /// MADX canonical variables
0077   /// IN : (x, theta_x, y, theta_y, xi) [m, rad, m, rad, 1]
0078   /// OUT : (x, y) [m, m]
0079   /// returns true if transport possible
0080   bool Transport2D(const double *in,
0081                    double *out,
0082                    bool check_apertures = false,
0083                    bool invert_beam_coord_sytems = true) const;
0084 
0085   bool Transport_m_GeV(double in_pos[3],
0086                        double in_momentum[3],
0087                        double out_pos[3],
0088                        double out_momentum[3],
0089                        bool check_apertures,
0090                        double z2_z1_dist) const;  ///< pos, momentum: x,y,z;  pos in m, momentum in GeV/c
0091 
0092   void PrintInputRange();
0093   bool CheckInputRange(const double *in, bool invert_beam_coord_sytems = true) const;
0094   void AddRectEllipseAperture(
0095       const LHCOpticsApproximator &in, double rect_x, double rect_y, double r_el_x, double r_el_y);
0096   void PrintOpticalFunctions();
0097   void PrintCoordinateOpticalFunctions(TMultiDimFet &parametrization,
0098                                        const std::string &coord_name,
0099                                        const std::vector<std::string> &input_vars);
0100 
0101   /**
0102      *\brief returns linearised transport matrix for x projection
0103      * |  dx_out/dx_in    dx_out/dthx_in   |
0104      * | dthx_out/dx_in  dthx_out/dthx_in  |
0105      *
0106      * input:  [m], [rad], xi:-1...0
0107      */
0108   void GetLineariasedTransportMatrixX(double mad_init_x,
0109                                       double mad_init_thx,
0110                                       double mad_init_y,
0111                                       double mad_init_thy,
0112                                       double mad_init_xi,
0113                                       TMatrixD &tr_matrix,
0114                                       double d_mad_x = 10e-6,
0115                                       double d_mad_thx = 10e-6);
0116 
0117   /**
0118      *\brief returns linearised transport matrix for y projection
0119      * |  dy_out/dy_in    dy_out/dthy_in   |
0120      * | dthy_out/dy_in  dthy_out/dthy_in  |
0121      *
0122      * input:  [m], [rad], xi:-1...0
0123      */
0124   void GetLineariasedTransportMatrixY(double mad_init_x,
0125                                       double mad_init_thx,
0126                                       double mad_init_y,
0127                                       double mad_init_thy,
0128                                       double mad_init_xi,
0129                                       TMatrixD &tr_matrix,
0130                                       double d_mad_y = 10e-6,
0131                                       double d_mad_thy = 10e-6);
0132 
0133   double GetDx(double mad_init_x,
0134                double mad_init_thx,
0135                double mad_init_y,
0136                double mad_init_thy,
0137                double mad_init_xi,
0138                double d_mad_xi = 0.001);
0139   double GetDxds(double mad_init_x,
0140                  double mad_init_thx,
0141                  double mad_init_y,
0142                  double mad_init_thy,
0143                  double mad_init_xi,
0144                  double d_mad_xi = 0.001);
0145   inline beam_type GetBeamType() const { return beam; }
0146 
0147   /// returns linear approximation of the transport parameterization
0148   /// takes numerical derivatives (see parameter ep) around point `atPoint' (this array has the same structure as `in' parameter in Transport method)
0149   /// the linearized transport: x = Cx + Lx*theta_x + vx*x_star
0150   void GetLinearApproximation(double atPoint[],
0151                               double &Cx,
0152                               double &Lx,
0153                               double &vx,
0154                               double &Cy,
0155                               double &Ly,
0156                               double &vy,
0157                               double &D,
0158                               double ep = 1E-5);
0159 
0160 private:
0161   void Init();
0162   double s_begin_;  ///< begin of transport along the reference orbit
0163   double s_end_;    ///< end of transport along the reference orbit
0164   beam_type beam;
0165   double nominal_beam_energy_;                  ///< GeV
0166   double nominal_beam_momentum_;                ///< GeV/c
0167   bool trained_;                                ///< trained polynomials
0168   std::vector<TMultiDimFet *> out_polynomials;  //! pointers to polynomials
0169   std::vector<std::string> coord_names;
0170   std::vector<LHCApertureApproximator> apertures_;  ///< apertures on the way
0171 
0172 #ifndef __CINT_
0173   friend class ProtonTransportFunctionsESSource;
0174 #endif  // __CINT__
0175 
0176   TMultiDimFet x_parametrisation;        ///< polynomial approximation for x
0177   TMultiDimFet theta_x_parametrisation;  ///< polynomial approximation for theta_x
0178   TMultiDimFet y_parametrisation;        ///< polynomial approximation for y
0179   TMultiDimFet theta_y_parametrisation;  ///< polynomial approximation for theta_y
0180 
0181   //train_mode mode_;  //polynomial selection mode - selection done by fitting function or selection from the list according to the specified order
0182   enum class VariableType { X, THETA_X, Y, THETA_Y };
0183   //internal methods
0184   void InitializeApproximators(polynomials_selection mode,
0185                                int max_degree_x,
0186                                int max_degree_tx,
0187                                int max_degree_y,
0188                                int max_degree_ty,
0189                                bool common_terms);
0190   void SetDefaultAproximatorSettings(TMultiDimFet &approximator, VariableType var_type, int max_degree);
0191   void SetTermsManually(TMultiDimFet &approximator, VariableType variable, int max_degree, bool common_terms);
0192 
0193   void AllocateErrorHists(TH1D *err_hists[4]);
0194   void AllocateErrorInputCorHists(TH2D *err_inp_cor_hists[4][5]);
0195   void AllocateErrorOutputCorHists(TH2D *err_out_cor_hists[4][5]);
0196 
0197   void DeleteErrorHists(TH1D *err_hists[4]);
0198   void DeleteErrorCorHistograms(TH2D *err_cor_hists[4][5]);
0199 
0200   void FillErrorHistograms(double errors[4], TH1D *err_hists[4]);
0201   void FillErrorDataCorHistograms(double errors[4], double var[5], TH2D *err_cor_hists[4][5]);
0202 
0203   void WriteHistograms(TH1D *err_hists[4],
0204                        TH2D *err_inp_cor_hists[4][5],
0205                        TH2D *err_out_cor_hists[4][5],
0206                        TFile *f_out,
0207                        std::string base_out_dir);
0208 
0209   ClassDef(LHCOpticsApproximator, 1)  // Proton transport approximator
0210 };
0211 
0212 class LHCApertureApproximator : public LHCOpticsApproximator {
0213 public:
0214   enum class ApertureType { NO_APERTURE, RECTELLIPSE };
0215 
0216   LHCApertureApproximator();
0217   LHCApertureApproximator(const LHCOpticsApproximator &in,
0218                           double rect_x,
0219                           double rect_y,
0220                           double r_el_x,
0221                           double r_el_y,
0222                           ApertureType type = ApertureType::RECTELLIPSE);
0223 
0224   bool CheckAperture(const double *in, bool invert_beam_coord_sytems = true) const;  //x, thx. y, thy, ksi
0225   //bool CheckAperture(MadKinematicDescriptor *in);  //x, thx. y, thy, ksi
0226 private:
0227   double rect_x_, rect_y_, r_el_x_, r_el_y_;
0228   ApertureType ap_type_;
0229 
0230   ClassDef(LHCApertureApproximator, 1)  // Aperture approximator
0231 };
0232 #endif  //TotemRPProtonTransportParametrization_LHC_OPTICS_APPROXIMATOR_H