tpc/generator/AnalyticFieldModel.cc

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