ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AnalyticFieldModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file AnalyticFieldModel.cc
1 
2 #include "AnalyticFieldModel.h"
3 #include "TFormula.h"
4 #include "TVector3.h"
5 
6 #include <cstdio>
7 
8 AnalyticFieldModel::AnalyticFieldModel(float _ifc_radius, float _ofc_radius, float _z_max, float scalefactor)
9 {
10  double ifc_radius = _ifc_radius;
11  double ofc_radius = _ofc_radius;
12  double tpc_halfz = _z_max;
13 
14  double sum = ifc_radius + ofc_radius; //338 in ALICE, [3] in args
15  double prod = ifc_radius * ofc_radius; //21250.75 in ALICE [4] in args
16  double diff = ofc_radius - ifc_radius;
17 
18  double a = ofc_radius * ofc_radius;
19  a *= (diff);
20  a *= (diff);
21  //a = a/1000.0;
22  a = 1 / a;
23  a *= scalefactor;
24  double b = 0.5;
25  double c = 1.0 / (((tpc_halfz) / 2.0) * ((tpc_halfz) / 2.0));
26  double d = sum;
27  double e = prod;
28 
29  vTestFunction1 = new TFormula("f1", "[0]*(x^4 - [3] *x^3 + [4] * x^2)*cos([1]* y)^2*exp(-1* [2] * z^2)");
30  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)))");
31 
32  erTestFunction1 = new TFormula("er", " [0]*(4*x^3 - 3 * [3] *x^2 + 2 * [4] * x)*cos([1]* y)^2*exp(-1* [2] * z^2)");
33  ePhiTestFunction1 = new TFormula("ePhi",
34  " [0]*(x^3 - [3] *x^2 + [4] * x)* -1 * [1] * sin(2 * [1]* y)*exp(-1* [2] * z^2)");
35  ezTestFunction1 = new TFormula("ez",
36  " [0]*(x^4 - [3] *x^3 + [4] * x^2)*cos([1]* y)^2*-1*2*[2]*z*exp(-1* [2] * z^2)");
37 
38  intErDzTestFunction1 = new TFormula("intErDz",
39  " [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]))) ");
40  intEPhiDzTestFunction1 = new TFormula("intEPhiDz",
41  "[0]* (x^3 - [3] *x^2 + [4] * x)* -1 * [1] * sin(2 * [1]* y)*((sqrt(pi)*TMath::Erf(sqrt([2]) * z))/(2 * sqrt([2])))");
42  intEzDzTestFunction1 = new TFormula("intEzDz",
43  "[0]* (x^4 - [3] *x^3 + [4] * x^2)*cos([1]* y)^2*exp(-1* [2] * z^2)");
44 
45  printf("Setting Analytic Formula, variables:\n");
46  printf("ifc=%f\tofc=%f\tdelz=%f\ndiff=%f\tscale=%f\n", ifc_radius, ofc_radius, tpc_halfz, diff, scalefactor);
47  printf("a=%E\nb=%E\nc=%E\nd=%f\ne=%f\n", a, b, c, d, e);
48 
49  vTestFunction1->SetParameters(a, b, c, d, e);
50  rhoTestFunction1->SetParameters(a, b, c, d, e);
51  //printf("rho value at (rmid,1,zmid)=%f\n",
52 
53  erTestFunction1->SetParameters(-a, b, c, d, e);
54  ePhiTestFunction1->SetParameters(-a, b, c, d, e);
55  ezTestFunction1->SetParameters(-a, b, c, d, e);
56  intErDzTestFunction1->SetParameters(-a, b, c, d, e);
57  intEPhiDzTestFunction1->SetParameters(-a, b, c, d, e);
58  intEzDzTestFunction1->SetParameters(-a, b, c, d, e);
59  return;
60 }
61 
62 TVector3 AnalyticFieldModel::E(TVector3 pos)
63 { //field as a function of position
64  //in rhat phihat zhat coordinates: (at phi=0, phi is the +Y position, Perp is the +X direction and Z is Z)
65  TVector3 ret(erTestFunction1->Eval(pos.Perp(), pos.Phi(), pos.Z()),
66  ePhiTestFunction1->Eval(pos.Perp(), pos.Phi(), pos.Z()),
67  ezTestFunction1->Eval(pos.Perp(), pos.Phi(), pos.Z()));
68  //now rotate this to the position we evaluated it at, to match the global coordinate system.
69  ret.RotateZ(pos.Phi());
70  return ret;
71 }
72 
73 double AnalyticFieldModel::Rho(TVector3 pos)
74 { //charge density as a function of position
75  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.
76  //at phi=0, phi is the +Y position, Perp is the +X direction and Z is Z.
77  return alice_chargescale * rhoTestFunction1->Eval(pos.Perp(), pos.Phi(), pos.Z());
78 }
79 
80 TVector3 AnalyticFieldModel::Eint(float zfinal, TVector3 pos)
81 { //field integral from 'pos' to z-position zfinal.
82  //in rhat phihat zhat coordinates: (at phi=0, phi is the +Y position, Perp is the +X direction and Z is Z)
83  TVector3 eintI, eintF;
84  eintI.SetXYZ(intErDzTestFunction1->Eval(pos.Perp(), pos.Phi(), pos.Z()),
85  intEPhiDzTestFunction1->Eval(pos.Perp(), pos.Phi(), pos.Z()),
86  intEzDzTestFunction1->Eval(pos.Perp(), pos.Phi(), pos.Z()));
87  eintF.SetXYZ(intErDzTestFunction1->Eval(pos.Perp(), pos.Phi(), zfinal),
88  intEPhiDzTestFunction1->Eval(pos.Perp(), pos.Phi(), zfinal),
89  intEzDzTestFunction1->Eval(pos.Perp(), pos.Phi(), zfinal));
90 
91  TVector3 ret = eintF - eintI;
92  // printf("Integrating z=%E to z=%E, delz=%E. Before rotation, field integrals: (xyz) (%E,%E,%E) to (%E,%E,%E), diff=(%E,%E,%E)\n",pos.Z(),zfinal,zfinal-pos.Z(),eintI.X(),eintI.Y(),eintI.Z(),eintF.X(),eintF.Y(),eintF.Z(),ret.X(),ret.Y(),ret.Z());
93  //now rotate this to the position we evaluated it at, to match the global coordinate system.
94  ret.RotateZ(pos.Phi());
95  return ret;
96 }