tpc/generator/Rossegger.h

0001 #ifndef TPCGENERATOR_ROSSEGGER_H
0002 #define TPCGENERATOR_ROSSEGGER_H
0003
0004 //
0005 //  Hello Space Charge Fans:  (although you should hate space charge)
0006 //
0007 //    This is a code that implements the calculations of space charge contributions
0008 //  In a cylindrical TPC.  It uses the solutions discovered by Stefan Rossegger for
0009 //  the ALICE experiment.  These solutions use various Bessel functions as a series
0010 //  solution to the "point charge in a conducting cylinder" problem imposed by
0011 //  The typical configuration of a TPC found at a collider.
0012 //
0013 //    These calculations have only a single dimension [length] and we shall choose
0014 //  the units of cm throughout the calculations.
0015 //
0016 //                                                TKH
0017 //                                                12-3-2015
0018 //
0019
0020 #include <cmath>
0021 #include <format>
0022 #include <iostream>
0023 #include <limits>
0024 #include <map>
0025 #include <string>
0026
0027 class TH2;
0028 class TH3;
0029
0030 #define NumberOfOrders 15  // Convergence problems after 15; Rossegger used 30
0031
0032 class Rossegger
0033 {
0034  public:
0035   explicit Rossegger(const std::string &filename);
0036   Rossegger(double InnerRadius = 30, double OuterRadius = 80, double Rdo_Z = 80, double precision = 1E-4);
0037   virtual ~Rossegger() = default;
0038
0039   void Verbosity(int v)
0040   {
0041     std::cout << std::format("verbosity set to {}.  was {}", v, verbosity) << std::endl;
0042     verbosity = v;
0043     return;
0044   };
0045   double Rmn(int m, int n, double r);               // Rmn function from Rossegger
0046   double Rmn_for_zeroes(int m, double x);           // Rmn function from Rossegger, as used to find Betamn zeroes.
0047   double Rmn1(int m, int n, double r);              // Rmn1 function from Rossegger
0048   double Rmn2(int m, int n, double r);              // Rmn2 function from Rossegger
0049   double RPrime(int m, int n, double ref, double r);  // RPrime function from Rossegger
0050
0051   double Rnk(int n, int k, double r);       // Rnk function from Rossegger
0052   double Rnk_for_zeroes(int n, double mu);  // Rnk function from Rossegger, as used to find munk zeroes.
0053
0054   double Limu(double mu, double x);  // Bessel functions of purely imaginary order
0055   double Kimu(double mu, double x);  // Bessel functions of purely imaginary order
0056
0057   double Ez(double r, double phi, double z, double r1, double phi1, double z1);
0058   double Er(double r, double phi, double z, double r1, double phi1, double z1);
0059   double Ephi(double r, double phi, double z, double r1, double phi1, double z1);
0060
0061   // alternate versions that don't use precalc constants.
0062   double Rmn_(int m, int n, double r);  // Rmn function from Rossegger
0063   // Rmn_for_zeroes doesn't have a way to speed it up with precalcs.
0064   double Rmn1_(int m, int n, double r);              // Rmn1 function from Rossegger
0065   double Rmn2_(int m, int n, double r);              // Rmn2 function from Rossegger
0066   double RPrime_(int m, int n, double ref, double r);  // RPrime function from Rossegger
0067   double Rnk_(int n, int k, double r);               // Rnk function from Rossegger
0068   double Rnk_for_zeroes_(int n, double mu);          // Rnk function from Rossegger, as used to find munk zeroes.
0069   double Ez_(double r, double phi, double z, double r1, double phi1, double z1);
0070   double Er_(double r, double phi, double z, double r1, double phi1, double z1);
0071   double Ephi_(double r, double phi, double z, double r1, double phi1, double z1);
0072
0073  protected:
0074   bool fByFile {false};
0075   double a {std::numeric_limits<double>::quiet_NaN()};
0076   double b {std::numeric_limits<double>::quiet_NaN()};
0077   double L {std::numeric_limits<double>::quiet_NaN()};  //  InnerRadius, OuterRadius, Length of 1/2 the TPC.
0078   int verbosity {0};
0079   double pi {M_PI};
0080   double epsilon {std::numeric_limits<double>::quiet_NaN()};  // precision.
0081
0082   bool tweak {false};
0083
0084   double MinimumDR {std::numeric_limits<double>::quiet_NaN()};
0085   double MinimumDPHI {std::numeric_limits<double>::quiet_NaN()};
0086   double MinimumDZ {std::numeric_limits<double>::quiet_NaN()};
0087
0088   double FindNextZero(double xstart, double epsilon, int order, double (Rossegger::*func)(int, double));  // Routine to find zeroes of func.
0089   void FindBetamn(double epsilon);                                                                        // Routine used to fill the Betamn array with resolution epsilon...
0090   void FindMunk(double epsilon);                                                                          // Routine used to fill the Munk array with resolution epsilon...
0091   bool CheckZeroes(double epsilon);                                                                       // confirm that the zeroes match to the desired precision.
0092
0093   void LoadZeroes(const std::string &destfile);
0094   void SaveZeroes(const std::string &destfile);
0095
0096   double Betamn[NumberOfOrders][NumberOfOrders]{};  //  Betamn array from Rossegger
0097   double N2mn[NumberOfOrders][NumberOfOrders]{};    //  N2mn array from Rossegger
0098   double Munk[NumberOfOrders][NumberOfOrders]{};    //  Munk array from Rossegger
0099   double N2nk[NumberOfOrders][NumberOfOrders]{};    //  N2nk array from Rossegger
0100
0101   void PrecalcFreeConstants();  // Routine used to fill the arrays of other values used repeatedly:
0102   // helpful values to precompute for frequent use:
0103   double BetaN[NumberOfOrders]{};                               //  BetaN=(n+1)*pi/L as used in eg 5.32, .46
0104   double BetaN_a[NumberOfOrders]{};                             //  BetaN*a as used in eg 5.32, .46
0105   double BetaN_b[NumberOfOrders]{};                             //  BetaN*b as used in eg 5.32, .46
0106   double km_BetaN_a[NumberOfOrders][NumberOfOrders]{};          // kn(m,BetaN*a) as used in Rossegger 5.32
0107   double im_BetaN_a[NumberOfOrders][NumberOfOrders]{};          // in(m,BetaN*a) as used in Rossegger 5.32
0108   double km_BetaN_b[NumberOfOrders][NumberOfOrders]{};          // kn(m,BetaN*b) as used in Rossegger 5.33
0109   double im_BetaN_b[NumberOfOrders][NumberOfOrders]{};          // in(m,BetaN*b) as used in Rossegger 5.33
0110   double bessel_denominator[NumberOfOrders][NumberOfOrders]{};  // BesselI(m,BetaN*a)*BesselK(m,BetaN*b)-BesselI(m,BetaN*b)*BesselK(m,BetaN*a) as in Rossegger 5.65
0111
0112   void PrecalcDerivedConstants();                           // Routine used to fill repeated values that depend on the Betamn and Munk zeroes:
0113   double ym_Betamn_a[NumberOfOrders][NumberOfOrders]{};     // yn(m,Betamn[m][n]*a) as used in Rossegger 5.11
0114   double jm_Betamn_a[NumberOfOrders][NumberOfOrders]{};     // jn(m,Betamn[m][n]*a) as used in Rossegger 5.11
0115   double liMunk_BetaN_a[NumberOfOrders][NumberOfOrders]{};  // limu(Munk[n][k],BetaN*a) as used in Rossegger 5.45
0116   double kiMunk_BetaN_a[NumberOfOrders][NumberOfOrders]{};  // kimu(Munk[n][k],BetaN*a) as used in Rossegger 5.45
0117   double sinh_Betamn_L[NumberOfOrders][NumberOfOrders]{};   // sinh(Betamn[m][n]*L)  as in Rossegger 5.64
0118   double sinh_pi_Munk[NumberOfOrders][NumberOfOrders]{};    // sinh(pi*Munk[n][k]) as in Rossegger 5.66
0119
0120   TH2 *Tags {nullptr};
0121   std::map<std::string, TH3 *> Grid;
0122 };
0123
0124 #endif /* __SPACECHARGE_H__ */