tpc/generator/AnalyticFieldModel.cc

0001
0002 #include "AnalyticFieldModel.h"
0003
0004 #include <TFormula.h>
0005 #include <TVector3.h>
0006
0007 #include <format>
0008 #include <iostream>
0009 #include <string>
0010
0011 AnalyticFieldModel::AnalyticFieldModel(float _ifc_radius, float _ofc_radius, float _z_max, float scalefactor)
0012   : vTestFunction1(new TFormula("f1", "[0]*(x^4 - [3] *x^3 + [4] * x^2)*cos([1]* y)^2*exp(-1* [2] * z^2)"))
0013   , rhoTestFunction1(new TFormula("ff1", "[0]*(((16.0 * x^2 - 9.0 * [3] * x + 4.0*[4]) *cos([1] * y)^2 * exp(-1 *[2]*z^2)) - ((x^2 -  [3] * x + [4]) * 2 * [1]^2 * cos(2 * [1] * y) * exp(-1 *[2]*z^2)) + ((x^4 -  [3] * x^3 + [4] * x^2) * cos([1] * y)^2 * (4*[2]^2*z^2 - 2 * [2]) * exp(-1 *[2]*z^2)))"))
0014   , erTestFunction1(new TFormula("er", " [0]*(4*x^3 - 3 * [3] *x^2 + 2 * [4] * x)*cos([1]* y)^2*exp(-1* [2] * z^2)"))
0015   , ePhiTestFunction1(new TFormula("ePhi", "  [0]*(x^3 - [3] *x^2 +  [4] * x)* -1  * [1] * sin(2 * [1]* y)*exp(-1* [2] * z^2)"))
0016   , ezTestFunction1(new TFormula("ez", " [0]*(x^4 - [3] *x^3 + [4] * x^2)*cos([1]* y)^2*-1*2*[2]*z*exp(-1* [2] * z^2)"))
0017   , intErDzTestFunction1(new TFormula("intErDz", " [0]*(4*x^3 - 3 * [3] *x^2 + 2 * [4] * x)*cos([1]* y)^2*((sqrt(pi)*TMath::Erf(sqrt([2]) * z))/(2 * sqrt([2]))) "))
0018   , intEPhiDzTestFunction1(new TFormula("intEPhiDz", "[0]* (x^3 - [3] *x^2 +  [4] * x)* -1  * [1] * sin(2 * [1]* y)*((sqrt(pi)*TMath::Erf(sqrt([2]) * z))/(2 * sqrt([2])))"))
0019   , intEzDzTestFunction1(new TFormula("intEzDz", "[0]* (x^4 - [3] *x^3 + [4] * x^2)*cos([1]* y)^2*exp(-1* [2] * z^2)"))
0020 {
0021   double ifc_radius = _ifc_radius;
0022   double ofc_radius = _ofc_radius;
0023   double tpc_halfz = _z_max;
0024
0025   double sum = ifc_radius + ofc_radius;   // 338 in ALICE, [3] in args
0026   double prod = ifc_radius * ofc_radius;  // 21250.75 in ALICE [4] in args
0027   double diff = ofc_radius - ifc_radius;
0028
0029   double a = ofc_radius * ofc_radius;
0030   a *= (diff);
0031   a *= (diff);
0032   // a =  a/1000.0;
0033   a = 1 / a;
0034   a *= scalefactor;
0035   double b = 0.5;
0036   double c = 1.0 / (((tpc_halfz) / 2.0) * ((tpc_halfz) / 2.0));
0037   double d = sum;
0038   double e = prod;
0039
0040   std::cout << "Setting Analytic Formula, variables:" << std::endl;
0041   std::cout << std::format("ifc={:.6f}\tofc={:.6f}\tdelz={:.6f}\ndiff={:.6f}\tscale={:.6f}",
0042                            ifc_radius, ofc_radius, tpc_halfz, diff, scalefactor)
0043             << std::endl;
0044   std::cout << std::format("a={:E}\nb={:E}\nc={:E}\nd={:.6f}\ne={:.6f}", a, b, c, d, e) << std::endl;
0045
0046   vTestFunction1->SetParameters(a, b, c, d, e);
0047   rhoTestFunction1->SetParameters(a, b, c, d, e);
0048
0049   erTestFunction1->SetParameters(-a, b, c, d, e);
0050   ePhiTestFunction1->SetParameters(-a, b, c, d, e);
0051   ezTestFunction1->SetParameters(-a, b, c, d, e);
0052   intErDzTestFunction1->SetParameters(-a, b, c, d, e);
0053   intEPhiDzTestFunction1->SetParameters(-a, b, c, d, e);
0054   intEzDzTestFunction1->SetParameters(-a, b, c, d, e);
0055   return;
0056 }
0057
0058 TVector3 AnalyticFieldModel::E(const TVector3& pos)
0059 {  // field as a function of position
0060   // in rhat phihat zhat coordinates: (at phi=0, phi is the +Y position, Perp is the +X direction and Z is Z)
0061   TVector3 ret(erTestFunction1->Eval(pos.Perp(), pos.Phi(), pos.Z()),
0062                ePhiTestFunction1->Eval(pos.Perp(), pos.Phi(), pos.Z()),
0063                ezTestFunction1->Eval(pos.Perp(), pos.Phi(), pos.Z()));
0064   // now rotate this to the position we evaluated it at, to match the global coordinate system.
0065   ret.RotateZ(pos.Phi());
0066   return ret;
0067 }
0068
0069 double AnalyticFieldModel::Rho(const TVector3& pos)
0070 {                                             // charge density as a function of position
0071   const double alice_chargescale = 8.85e-14;  // their rho has charge density in units of C/cm^3 /eps0.  This is eps0 in (V*cm)/C units so that I can multiple by the volume in cm^3 to get Q in C.
0072   // at phi=0, phi is the +Y position, Perp is the +X direction and Z is Z.
0073   return alice_chargescale * rhoTestFunction1->Eval(pos.Perp(), pos.Phi(), pos.Z());
0074 }
0075
0076 TVector3 AnalyticFieldModel::Eint(float zfinal, const TVector3& pos)
0077 {  // field integral from 'pos' to z-position zfinal.
0078   // in rhat phihat zhat coordinates: (at phi=0, phi is the +Y position, Perp is the +X direction and Z is Z)
0079   TVector3 eintI;
0080   TVector3 eintF;
0081   eintI.SetXYZ(intErDzTestFunction1->Eval(pos.Perp(), pos.Phi(), pos.Z()),
0082                intEPhiDzTestFunction1->Eval(pos.Perp(), pos.Phi(), pos.Z()),
0083                intEzDzTestFunction1->Eval(pos.Perp(), pos.Phi(), pos.Z()));
0084   eintF.SetXYZ(intErDzTestFunction1->Eval(pos.Perp(), pos.Phi(), zfinal),
0085                intEPhiDzTestFunction1->Eval(pos.Perp(), pos.Phi(), zfinal),
0086                intEzDzTestFunction1->Eval(pos.Perp(), pos.Phi(), zfinal));
0087
0088   TVector3 ret = eintF - eintI;
0089   // std::cout <<  boost::str(boost::format("Integrating z=%E to z=%E, delz=%E.  Before rotation, field integrals: (xyz)  (%E,%E,%E) to (%E,%E,%E), diff=(%E,%E,%E)") %pos.Z() %zfinal %zfinal-pos.Z() %eintI.X() %eintI.Y() %eintI.Z() %eintF.X() %eintF.Y() %eintF.Z() %ret.X() %ret.Y() %ret.Z()) << std::endl;
0090   // now rotate this to the position we evaluated it at, to match the global coordinate system.
0091   ret.RotateZ(pos.Phi());
0092   return ret;
0093 }