GeFiCa
Germanium detector Field Calculator
Public Member Functions | Public Attributes | Protected Member Functions | Protected Attributes | List of all members
GeFiCa::Grid Class Reference

Data structure of a electric field grid. More...

#include <Grid.h>

Inheritance diagram for GeFiCa::Grid:
GeFiCa::Points GeFiCa::R GeFiCa::Rho GeFiCa::RhoPhi GeFiCa::RhoZ GeFiCa::RTheta GeFiCa::X GeFiCa::XYZ

Public Member Functions

 Grid (size_t n1=0, size_t n2=0, size_t n3=0)
 Default constructor. More...
 
virtual ~Grid ()
 
virtual void SetupWith (Detector &detector)
 Fix potentials on boundaries based on. More...
 
void SuccessiveOverRelax ()
 Successively over-relax potentials on grid points. More...
 
void SolveAnalytically ()
 Solve Poisson's Equation analytically. More...
 
double GetV (double c1, double c2=0, double c3=0) const
 Get potential at (c1,c2,c3) by interpolation. More...
 
double GetE (double c1, double c2=0, double c3=0) const
 
double GetE1 (double c1, double c2=0, double c3=0) const
 
double GetE2 (double c1, double c2=0, double c3=0) const
 
double GetE3 (double c1, double c2=0, double c3=0) const
 
double GetC ()
 Get detector capacitance. More...
 
TTree * GetTree (bool createNew=false)
 Create &/or return a TTree with field data. More...
 
bool IsDepleted ()
 Check if every grid point is depleted. More...
 
Gridoperator*= (double scale)
 Potentials at all points are multiplied by. More...
 
Gridoperator+= (Grid &other)
 Potentials of this grid are summed with those of. More...
 
virtual FieldLineGetFieldLineFrom (double c1, double c2, double c3=0, bool positive=true)
 Propogate a field line from (c1,c2,c3). More...
 
- Public Member Functions inherited from GeFiCa::Points
size_t GetN ()
 total number of points More...
 

Public Attributes

size_t N1
 number of points along the 1st coordinate More...
 
size_t N2
 number of points along the 2nd coordinate More...
 
size_t N3
 number of points along the 3rd coordinate More...
 
std::vector< double > Src
 -(net impurity concentration)x|Qe|/epsilon More...
 
size_t MaxIterations
 maximal iterations of SOR to be performed More...
 
double RelaxationFactor
 within (0,2), used to speed up convergence More...
 
double Tolerance
 SOR stops when error<Tolerance. More...
 
size_t Iterations
 number of iterations of SOR performed More...
 
- Public Attributes inherited from GeFiCa::Points
std::vector< double > C1
 the 1st coordinates of the points More...
 
std::vector< double > C2
 the 2nd coordinates of the points More...
 
std::vector< double > C3
 the 3rd coordinates of the points More...
 
std::vector< double > Vp
 potential at each point More...
 
std::vector< double > Et
 total electric field strength More...
 
std::vector< double > E1
 projection of Et on C1 More...
 
std::vector< double > E2
 projection of Et on C2 More...
 
std::vector< double > E3
 projection of Et on C3 More...
 
std::vector< double > dC1p
 step length to next point alone C1 More...
 
std::vector< double > dC1m
 step length to previous point alone C1 More...
 
std::vector< double > dC2p
 step length to next point along C2 More...
 
std::vector< double > dC2m
 step length to previous point along C2 More...
 
std::vector< double > dC3p
 step length to next point alone C3 More...
 
std::vector< double > dC3m
 step length to previous point alone C3 More...
 

Protected Member Functions

virtual void OverRelaxAt (size_t idx)
 Over relax potential Vp[. More...
 
size_t GetIdxOfPointToTheRightOf (double c1, size_t begin, size_t end) const
 Get index of point near. More...
 
size_t GetIdxOfPointToTheRightOf (double c1, double c2, size_t begin, size_t end) const
 
size_t GetIdxOfPointToTheRightOf (double c1, double c2, double c3, size_t begin, size_t end) const
 
virtual double GetData (const std::vector< double > &data, double c1, double c2, double c3) const
 Interpolate grid data at (c1,c2,c3). More...
 
size_t GetIdxOfMaxV ()
 Get index of the grid point with max potential. More...
 
size_t GetIdxOfMinV ()
 Get index of the grid point with min potential. More...
 
virtual void CalculateE ()
 Calculate Et, E1, E2, E3 from Vp. More...
 
double twopoint (double dataset[2], double tarlocationset, double pointxset[2]) const
 Calculate interpolation value between two point. More...
 
double threepoint (double dataset[3], double tarlocationset[2], double pointxset[3], double pointyset[3]) const
 Calculate interpolation value between three point (triangle) More...
 
double fourpoint (double dataset[4], double tarlocationset[2], double pointxset[4], double pointyset[4]) const
 Calculate interpolation value between four point (rectangle) More...
 
 ClassDef (Grid, 1)
 

Protected Attributes

std::vector< bool > fIsFixed
 true if field values are fixed More...
 
std::vector< bool > fIsDepleted
 true if a grid point is depleted More...
 
TTree * fTree
 ! ROOT tree to visualize fields More...
 
DetectorfDetector
 ! Pointer to associated detector object More...
 

Detailed Description

Data structure of a electric field grid.

It inherits the flat data structure from Points instead of using the class to save typing. Compare C1[i] with point.C1[i].

Definition at line 50 of file Grid.h.

Constructor & Destructor Documentation

§ Grid()

Grid::Grid ( size_t  n1 = 0,
size_t  n2 = 0,
size_t  n3 = 0 
)

Default constructor.

It also defines ROOT drawing style.

Definition at line 26 of file Grid.cc.

26  : Points(), TNamed("grid","grid"),
27  N1(n1), N2(n2), N3(n3), MaxIterations(5000), RelaxationFactor(1.95),
28  Tolerance(1e-7*volt), Iterations(0), fTree(0), fDetector(0)
29 {
30  // pick up a good style to modify
31  gStyle->SetName("GeFiCa");
32  gStyle->SetLegendBorderSize(0);
33  gStyle->SetLegendFont(132);
34  gStyle->SetLabelFont(132,"XYZ");
35  gStyle->SetTitleFont(132,"XYZ");
36  gStyle->SetLabelSize(0.05,"XYZ");
37  gStyle->SetTitleSize(0.05,"XYZ");
38  gStyle->SetPadTopMargin(0.02);
39  // create a smoother palette than the default one
40  const int nRGBs = 5;
41  const int nCont = 255;
42  double stops[nRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
43  double red[nRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
44  double green[nRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
45  double blue[nRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
46  TColor::CreateGradientColorTable(nRGBs, stops, red, green, blue, nCont);
47  gStyle->SetNumberContours(nCont);
48 }
size_t Iterations
number of iterations of SOR performed
Definition: Grid.h:60
double Tolerance
SOR stops when error<Tolerance.
Definition: Grid.h:59
size_t N2
number of points along the 2nd coordinate
Definition: Grid.h:54
TTree * fTree
! ROOT tree to visualize fields
Definition: Grid.h:130
size_t N1
number of points along the 1st coordinate
Definition: Grid.h:53
A group of discrete points.
Definition: Grid.h:9
size_t MaxIterations
maximal iterations of SOR to be performed
Definition: Grid.h:57
Detector * fDetector
! Pointer to associated detector object
Definition: Grid.h:131
size_t N3
number of points along the 3rd coordinate
Definition: Grid.h:55
static const double volt
Definition: Units.h:20
double RelaxationFactor
within (0,2), used to speed up convergence
Definition: Grid.h:58

§ ~Grid()

virtual GeFiCa::Grid::~Grid ( )
inlinevirtual

Definition at line 66 of file Grid.h.

66 {};

Member Function Documentation

§ CalculateE()

void Grid::CalculateE ( )
protectedvirtual

Calculate Et, E1, E2, E3 from Vp.

Reimplemented in GeFiCa::XYZ, GeFiCa::RhoZ, GeFiCa::RhoPhi, and GeFiCa::RTheta.

Definition at line 551 of file Grid.cc.

552 {
553  for (size_t i=0; i<GetN(); i++) { // deal with E1 only
554  E1[i]=-(Vp[i+1]-Vp[i-1])/(dC1p[i]+dC1m[i]); Et[i]=E1[i];
555  }
556  E1[0]=-(Vp[1]-Vp[0])/dC1p[0]; Et[0]=E1[0];
557  E1[N1-1]=-(Vp[N1-1]-Vp[N1-2])/dC1m[N1-1]; Et[N1-1]=E1[N1-1];
558 }
std::vector< double > dC1m
step length to previous point alone C1
Definition: Grid.h:21
std::vector< double > E1
projection of Et on C1
Definition: Grid.h:17
size_t N1
number of points along the 1st coordinate
Definition: Grid.h:53
std::vector< double > dC1p
step length to next point alone C1
Definition: Grid.h:20
size_t GetN()
total number of points
Definition: Grid.h:26
std::vector< double > Vp
potential at each point
Definition: Grid.h:15
std::vector< double > Et
total electric field strength
Definition: Grid.h:16

§ ClassDef()

GeFiCa::Grid::ClassDef ( Grid  ,
 
)
protected

§ fourpoint()

double Grid::fourpoint ( double  dataset[4],
double  tarlocationset[2],
double  pointxset[4],
double  pointyset[4] 
) const
protected

Calculate interpolation value between four point (rectangle)

Definition at line 604 of file Grid.cc.

606 {
607  //0---------1
608  //| |
609  //| |
610  //| |
611  //2---------3
612  double ab=(pointxset[1]-tarlocationset[0])/(pointxset[1]-pointxset[0]);
613  double aa=1-ab;
614  double ba=(tarlocationset[1]-pointyset[3])/(pointyset[1]-pointyset[3]);
615  double bb=1-ba;
616 
617  double result=(dataset[1]*aa+dataset[0]*ab)*ba+(dataset[2]*ab+dataset[3]*aa)*bb;
618  delete [] dataset;
619  delete [] tarlocationset;
620  delete [] pointxset;
621  delete [] pointyset;
622 
623  if (TMath::Abs(result)<1e-10) return 0;
624  return result;
625 }

§ GetC()

double Grid::GetC ( )

Get detector capacitance.

Calculate C based on \(CV^2/2 = \epsilon \int E^2/2 dx^3\). http://hyperphysics.phy-astr.gsu.edu/hbase/electric/capeng.html

Definition at line 198 of file Grid.cc.

199 {
200  Info("GetC","Start...");
201  SuccessiveOverRelax(); // identify undepleted region
202  std::vector<double>original=Src; // save original rho/epsilon
203  for (size_t i=0; i<GetN(); i++) {
204  Src[i]=0; // set impurity to zero
205  if (fIsDepleted[i]==false) fIsFixed[i]=true; // fix undepleted points
206  }
207  SuccessiveOverRelax(); // calculate potential without impurity
208  Src=original; // set impurity back
209 
210  return 0;
211 }
std::vector< bool > fIsFixed
true if field values are fixed
Definition: Grid.h:128
std::vector< double > Src
-(net impurity concentration)x|Qe|/epsilon
Definition: Grid.h:56
size_t GetN()
total number of points
Definition: Grid.h:26
void SuccessiveOverRelax()
Successively over-relax potentials on grid points.
Definition: Grid.cc:120
std::vector< bool > fIsDepleted
true if a grid point is depleted
Definition: Grid.h:129

§ GetData()

double Grid::GetData ( const std::vector< double > &  data,
double  c1,
double  c2,
double  c3 
) const
protectedvirtual

Interpolate grid data at (c1,c2,c3).

Definition at line 260 of file Grid.cc.

262 {
263  if (N3!=0) { // 3D
264  Warning("GetData", "not yet implemented. Abort."); abort();
265  }
266  if (N2!=0) { // 2D
267  // https://codeplea.com/triangular-interpolation
268  // tmv
269  // +-aa--+--ab--+(C1[idx], C2[idx])
270  // | ^ |
271  // | dyp| ba
272  // | |(x,y) |
273  // +<----+----->+
274  // | | |
275  // | dym| bb
276  // | v |
277  // +-----+------+
278  // bmv
279  size_t idx=GetIdxOfPointToTheRightOf(x,y,0,N2-1); // always exists
280 
281  bool tl=false; // existence of top left grid point
282  bool br=false; // existence of bottom right grid point
283  bool bl=false; // existence of bottom left grid point
284  if (idx%N1!=0) tl=true; // not left boundary
285  if (idx>=N1) br=true; // not bottom boundary
286  if (tl&&bl) bl=true; // neither left nor bottom boundary
287 
288  if(!tl&&!br&&!bl) { // bottom left corner
289  return data[idx];
290  } else if(tl&&!br&&!bl) { // bottom boundary
291  return twopoint(new double[2] {data[idx-1],data[idx]},
292  x,new double[2] {C1[idx-1],C1[idx]});
293  } else if(!tl&&!bl&&br) { // left boundary
294  return twopoint(new double[2]{data[idx-N1],data[idx]},
295  y,new double[2]{C2[idx-N1],C2[idx]});
296  }
297  // no boundary case
298  if(dC1p[idx-1]==dC1m[idx]&&dC1p[idx-N1-1]==dC1m[idx-N1]&&
299  dC2p[idx-N1-1]==dC2m[idx-1]&&dC2p[idx-N1]==dC2m[idx]) {
300  return fourpoint(new double[4]{data[idx-1],data[idx],data[idx-N1-1],data[idx-N1]},new double[2]{x,y},new double[4]{C1[idx-1],C1[idx],C1[idx-N1-1],C1[idx-N1]},new double[4]{C2[idx-1],C2[idx],C2[idx-N1-1],C2[idx-N1]});
301  }
302 
303  //top left case o
304  if(dC1p[idx-1]!=dC1m[idx]&&
305  dC1p[idx-N1-1]==dC1m[idx-N1]&&
306  dC2p[idx-N1-1]!=dC2m[idx-1]&&
307  dC2p[idx-N1]==dC2m[idx])
308  {
309  double xb=dC1m[idx]<dC1p[idx-1] ? C1[idx]-dC1m[idx] : C1[idx-1]+dC1p[idx-1];
310  double yb=dC2m[idx-1]<dC2p[idx-N1-1] ? C2[idx-1]-dC2m[idx-1] : C2[idx-N1-1]+dC2p[idx-N1-1];
311  double bmv=(C1[idx-N1]-xb)/(C1[idx-N1]-C1[idx-N1-1])*data[idx-N1-1]+(xb-C1[idx-N1-1])/(C1[idx-N1]-C1[idx-N1-1])*data[idx-N1];
312  double tmv=fIsFixed[idx]?data[idx]:data[idx-1];
313  double mlv=fIsFixed[idx-1]?data[idx-1]:data[idx-N1-1];
314  double mmv=(C2[idx]-yb)/(C2[idx]-C2[idx-N1])*bmv+(yb-C2[idx-N1])/(C2[idx]-C2[idx-N1])*tmv;
315  if(x>xb)
316  {
317  return fourpoint(new double[4]{tmv,data[idx],bmv,data[idx-N1]},new double[2]{x,y},new double[4]{xb,C1[idx],xb,C1[idx-N1]},new double[4]{C2[idx],C2[idx],C2[idx-N1],C2[idx-N1]});
318  }
319  else if(y<yb)
320  {
321  return fourpoint(new double[4]{mlv,mmv,data[idx-N1-1],bmv},new double[2]{x,y},new double[4]{C1[idx-1],xb,C1[idx-1],xb},new double[4]{yb,yb,C2[idx-N1],C2[idx-N1]});
322  }
323  else if((C2[idx]-yb)/(xb-C1[idx-1]*(x-C1[idx-1])+yb>y))
324  {
325  return threepoint(new double[3]{tmv,mmv,mlv},new double[2]{x,y},new double[3]{xb,xb,C1[idx-1]},new double[3]{C2[idx],yb,yb});
326  }
327  else if((C2[idx]-yb)/(xb-C1[idx-1]*(x-C1[idx-1])+yb<y))
328  {
329  return threepoint(new double[3]{tmv,tmv,mlv},new double[2]{x,y},new double[3]{xb,C1[idx-1],C1[idx-1]},new double[3]{C2[idx],C2[idx],yb});
330  }
331  }
332  //topright case o
333  if(dC1p[idx-1]!=dC1m[idx]&&
334  dC1p[idx-N1-1]==dC1m[idx-N1]&&
335  dC2p[idx-N1-1]==dC2m[idx-1]&&
336  dC2p[idx-N1]!=dC2m[idx])
337  {
338  double xb=dC1m[idx]<dC1p[idx-1] ? C1[idx]-dC1m[idx] : C1[idx-1]+dC1p[idx-1];
339  double yb=dC2m[idx]<dC2p[idx-N1] ? C2[idx]-dC2m[idx] : C2[idx-N1]+dC2p[idx-N1];
340  double bmv=(C1[idx-N1]-xb)/(C1[idx-N1]-C1[idx-N1-1])*data[idx-N1-1]+(xb-C1[idx-N1-1])/(C1[idx-N1]-C1[idx-N1-1])*data[idx-N1];
341  double tmv=fIsFixed[idx]?data[idx]:data[idx-1];
342  double mrv=fIsFixed[idx]?data[idx]:data[idx-N1];
343  double mmv=(C2[idx]-yb)/(C2[idx]-C2[idx-N1])*bmv+(yb-C2[idx-N1])/(C2[idx]-C2[idx-N1])*tmv;
344  if(x<xb)
345  {
346  return fourpoint(new double[4]{data[idx-1],tmv,data[idx-N1-1],bmv},new double[2]{x,y},new double[4]{C1[idx-1],xb,C1[idx-N1-1],xb},new double[4]{C2[idx],C2[idx],C2[idx-N1],C2[idx-N1]});
347  }
348  else if(y<yb)
349  {
350  return fourpoint(new double[4]{mmv,mrv,bmv,data[idx-N1]},new double[2]{x,y},new double[4]{xb,C1[idx],xb,C1[idx-1]},new double[4]{yb,yb,C2[idx-N1],C2[idx-N1]});
351  }
352  else if((C2[idx]-yb)/(C1[idx]-x*(x-C1[idx])+yb>y))
353  {
354  return threepoint(new double[3]{tmv,mmv,mrv},new double[2]{x,y},new double[3]{xb,xb,C1[idx]},new double[3]{C2[idx],yb,yb});
355  }
356  else if((C2[idx]-yb)/(C1[idx]-xb*(x-C1[idx])+yb<y))
357  {
358  return threepoint(new double[3]{tmv,data[idx],mrv},new double[2]{x,y},new double[3]{xb,C1[idx],C1[idx]},new double[3]{C2[idx],C2[idx],yb});
359  }
360  }
361  //bottom left case o
362  if(dC1p[idx-1]==dC1m[idx]&&
363  dC1p[idx-N1-1]!=dC1m[idx-N1]&&
364  dC2p[idx-N1-1]!=dC2m[idx-1]&&
365  dC2p[idx-N1]==dC2m[idx])
366  {
367  double xb=dC1m[idx-N1]<dC1p[idx-N1-1] ? C1[idx-N1]-dC1m[idx-N1] : C1[idx-N1-1]+dC1p[idx-N1-1];
368  double yb=dC2m[idx-1]<dC2p[idx-N1-1] ? C2[idx-1]-dC2m[idx-1] : C2[idx-N1-1]+dC2p[idx-N1-1];
369  double tmv=(C1[idx]-xb)/(C1[idx]-C1[idx-1])*data[idx-1]+(xb-C1[idx-1])/(C1[idx]-C1[idx-1])*data[idx];
370  double bmv=fIsFixed[idx]?data[idx]:data[idx-1];
371  double mlv=fIsFixed[idx-1]?data[idx-1]:data[idx-N1-1];
372  double mmv=(C2[idx]-yb)/(C2[idx]-C2[idx-N1])*bmv+(yb-C2[idx-N1])/(C2[idx]-C2[idx-N1])*tmv;
373  if(x>xb)
374  {
375  return fourpoint(new double[4]{tmv,data[idx],bmv,data[idx-N1]},new double[2]{x,y},new double[4]{xb,C1[idx],xb,C1[idx-N1]},new double[4]{C2[idx],C2[idx],C2[idx-N1],C2[idx-N1]});
376  }
377  else if(y>yb)
378  {
379  return fourpoint(new double[4]{data[idx-N1-1],tmv,mlv,mmv},new double[2]{x,y},new double[4]{C1[idx-1],xb,C1[idx-1],xb},new double[4]{C2[idx],C2[idx],yb,yb});
380  }
381  else if((C2[idx]-yb)/(xb-C1[idx-1]*(x-C1[idx-1])+C2[idx-N1]>y))
382  {
383  return threepoint(new double[3]{bmv,mmv,mlv},new double[2]{x,y},new double[3]{xb,xb,C1[idx-1]},new double[3]{yb,yb,C2[idx-N1]});
384  }
385  else if((C2[idx]-yb)/(xb-C1[idx-1]*(x-C1[idx-1])+C2[idx-N1]<y))
386  {
387  return threepoint(new double[3]{data[idx-N1-1],bmv,mlv},new double[2]{x,y},new double[3]{C1[idx-1],xb,C1[idx-1]},new double[3]{C2[idx],C2[idx],yb});
388  }
389  }
390  //bottom right case o
391  if(dC1p[idx-1]==dC1m[idx]&&
392  dC1p[idx-N1-1]!=dC1m[idx-N1]&&
393  dC2p[idx-N1-1]==dC2m[idx-1]&&
394  dC2p[idx-N1]!=dC2m[idx])
395  {
396  double xb=dC1m[idx-N1]<dC1p[idx-N1-1] ? C1[idx-N1]-dC1m[idx-N1] : C1[idx-N1-1]+dC1p[idx-N1-1];
397  double yb=dC2m[idx]<dC2p[idx-N1] ? C2[idx]-dC2m[idx] : C2[idx-N1]+dC2p[idx-N1];
398  double tmv=(C1[idx]-xb)/(C1[idx]-C1[idx-1])*data[idx-1]+(xb-C1[idx-1])/(C1[idx]-C1[idx-1])*data[idx];
399  double bmv=fIsFixed[idx-N1]?data[idx-N1]:data[idx-N1-1];
400  double mrv=fIsFixed[idx-N1]?data[idx-N1]:data[idx];
401  double mmv=(C2[idx]-yb)/(C2[idx]-C2[idx-N1])*bmv+(yb-C2[idx-N1])/(C2[idx]-C2[idx-N1])*tmv;
402  if(x<xb)
403  {
404  return fourpoint(new double[4]{data[idx-1],tmv,data[idx-N1-1],bmv},new double[2]{x,y},new double[4]{C1[idx-1],xb,C1[idx-N1-1],xb},new double[4]{C2[idx],C2[idx],C2[idx-N1],C2[idx-N1]});
405  }
406  else if(y>yb)
407  {
408  return fourpoint(new double[4]{tmv,data[idx],mmv,mrv},new double[2]{x,y},new double[4]{xb,C1[idx],xb,C1[idx]},new double[4]{C2[idx],C2[idx],yb,yb});
409  }
410  else if((yb-C2[idx-N1])/(C1[idx]-xb*(x-xb)+C2[idx-N1]>y))
411  {
412  return threepoint(new double[3]{bmv,mmv,mrv},new double[2]{x,y},new double[3]{xb,xb,C1[idx]},new double[3]{C2[idx-N1],yb,yb});
413  }
414  else if((yb-C2[idx-N1])/(C1[idx]-xb*(x-xb)+C2[idx-N1]<=y))
415  {
416  return threepoint(new double[3]{bmv,data[idx-N1],mrv},new double[2]{x,y},new double[4]{xb,C1[idx],C1[idx-1]},new double[4]{C2[idx-N1],C2[idx-N1],yb});
417  }
418  }
419  //cross top down case o
420  if(dC1p[idx-1]!=dC1m[idx]&&
421  dC1p[idx-N1-1]!=dC1m[idx-N1]&&
422  dC2p[idx-N1-1]==dC2m[idx-1]&&
423  dC2p[idx-N1]==dC2m[idx])
424  {
425  double topb=dC1m[idx]<dC1p[idx-1] ? C1[idx]-dC1m[idx] : C1[idx-1]+dC1p[idx-1];
426  double bottomb=dC1m[idx-N1]<dC1p[idx-N1-1] ? C1[idx-N1]-dC1m[idx-N1] : C1[idx-N1-1]+dC1p[idx-N1-1];
427  double bmv=fIsFixed[idx-N1]?data[idx-N1]:data[idx-N1-1];
428  double tmv=fIsFixed[idx]?data[idx]:data[idx-1];
429  double TopValueWithBottomCorssPoint=topb>bottomb ? (topb-bottomb)/(topb-C1[idx-1])*data[idx-1]+(bottomb-C1[idx-1])/(topb-C1[idx-1])*tmv :
430  (bottomb-topb)/(C1[idx]-topb)*data[idx]+(C1[idx]-bottomb)/(C1[idx-1]-topb)*tmv ;
431  double BottomValueWithTopCorssPoint=topb>bottomb ? (topb-bottomb)/(C1[idx-N1]-bottomb)*bmv+(C1[idx-N1]-topb)/(C1[idx-N1]-bottomb)*data[idx-N1] :
432  (topb-C1[idx-N1-1])/(bottomb-C1[idx-N1-1])*bmv+(bottomb-topb)/(bottomb-C1[idx-1-N1])*data[idx-N1-1] ;
433  double TopLeft,TopRight,BottomLeft,BottomRight;
434  double TopLeftV,TopRightV,BottomLeftV,BottomRightV;
435  if(topb>bottomb) {
436  TopLeft=bottomb;
437  TopLeftV=TopValueWithBottomCorssPoint;
438  TopRight=topb;
439  TopRightV=tmv;
440  BottomLeft=bottomb;
441  BottomLeftV=bmv;
442  BottomRight=topb;
443  BottomRightV=BottomValueWithTopCorssPoint;
444  } else {
445  TopRight=bottomb;
446  TopRightV=TopValueWithBottomCorssPoint;
447  TopLeft=topb;
448  TopLeftV=tmv;
449  BottomRight=bottomb;
450  BottomRightV=bmv;
451  BottomLeft=topb;
452  BottomLeftV=BottomValueWithTopCorssPoint;
453  }
454  if(x>TopRight) {
455  if(x<1.185&&x>1.175&&y>4.915&&y<4.925)Printf("trv:%f,brv:%f",TopRightV,BottomRightV);
456  return fourpoint(new double[4]{TopRightV,data[idx],BottomRightV,data[idx-N1]},
457  new double[2]{x,y},new double[4]{TopRight,C1[idx],BottomRight,C1[idx-N1]},
458  new double[4]{C2[idx],C2[idx],C2[idx-N1],C2[idx-N1]});
459  } else if(x<TopLeft) {
460  return fourpoint(new double[4]{data[idx-1],TopLeftV,data[idx-N1-1],BottomLeftV},
461  new double[2]{x,y},new double[4]{C1[idx-1],TopLeft,C1[idx-N1-1],BottomLeft},
462  new double[4]{C2[idx-1],C2[idx-1],C2[idx-N1-1],C2[idx-N1-1]});
463  } else if((dC2m[idx])/((topb-bottomb)*(x-(TopRight-TopLeft)/2)+C2[idx-N1]>y-dC2m[idx]/2)) {
464  return topb>bottomb ? threepoint(new double[3]{TopRightV,BottomLeftV,BottomRightV},new double[2]{x,y},new double[3]{TopRight,BottomLeft,BottomRight},new double[3]{C2[idx],C2[idx-N1],C2[idx-N1]}) :
465  threepoint(new double[3]{TopLeft,BottomLeftV,BottomRightV},new double[2]{x,y},new double[3]{TopLeft,BottomLeft,BottomRight},new double[3]{C2[idx],C2[idx-N1],C2[idx-N1]});
466  } else if((dC2m[idx])/((topb-bottomb)*(x-(TopRight-TopLeft)/2)+C2[idx-N1]<=y-dC2m[idx]/2)) {
467  return topb<bottomb ? threepoint(new double[3]{TopRightV,TopLeftV,BottomRightV},new double[2]{x,y},new double[3]{TopRight,TopLeft,BottomRight},new double[3]{C2[idx],C2[idx],C2[idx-N1]}) :
468  threepoint(new double[3]{TopLeft,TopLeftV,BottomLeftV},new double[2]{x,y},new double[3]{TopLeft,BottomLeft,BottomLeft},new double[3]{C2[idx],C2[idx],C2[idx-N1]});
469  }
470  }
471  // cross left right case:
472  // +----+ ... C2[idx]
473  // | |
474  // tlv + trv+ ... t2
475  // |\ |
476  // | \ |
477  // | \ |
478  // | \|
479  // blv + brv+ ... b2
480  // | |
481  // | |
482  // +----+ ... C2[idx-N1]
483  // . .
484  // . .
485  // C1[idx-1] C1[idx]
486  if(dC1p[idx-1]==dC1m[idx]&&
487  dC1p[idx-N1-1]==dC1m[idx-N1]&&
488  dC2p[idx-N1-1]!=dC2m[idx-1]&&
489  dC2p[idx-N1]!=dC2m[idx])
490  {
491  double leftb=dC2m[idx-1]<dC2p[idx-N1-1] ? C2[idx-1]-dC2m[idx-1] : C2[idx-N1-1]+dC2p[idx-N1-1]; // y of left intersection
492  double rightb=dC2m[idx]<dC2p[idx-N1] ? C2[idx]-dC2m[idx] : C2[idx-N1]+dC2p[idx-N1]; // y of right intersection
493  double lmv=fIsFixed[idx-1]?data[idx-1]:data[idx-N1-1]; // voltage of leftb
494  double rmv=fIsFixed[idx]?data[idx]:data[idx-N1]; // voltage of rightb
495  double LeftValueWithRightCorssPoint=leftb>rightb ?
496  (leftb-rightb)/(leftb-C2[idx-N1-1])*data[idx-N1-1]+(rightb-C2[idx-N1-1])/(leftb-C2[idx-N1-1])*lmv :
497  (rightb-leftb)/(C2[idx-1]-leftb)*data[idx-1]+(C2[idx-1]-rightb)/(C2[idx-1]-leftb)*lmv;
498  double RightValueWithLeftCorssPoint=leftb>rightb ?
499  (rightb-leftb)/(C2[idx]-rightb)*data[idx]+(C2[idx]-leftb)/(C2[idx]-rightb)*rmv :
500  (leftb-C2[idx-N1])/(rightb-C2[idx-N1])*rmv+(rightb-leftb)/(rightb-C2[idx-N1])*data[idx-N1];
501  double t2,b2,tlv,trv,blv,brv;
502 
503  if(leftb>rightb) {
504  t2=leftb;
505  tlv=lmv;
506  trv=RightValueWithLeftCorssPoint;
507  b2=rightb;
508  blv=LeftValueWithRightCorssPoint;
509  brv=rmv;
510  } else {
511  trv=rmv;
512  t2=rightb;
513  tlv=LeftValueWithRightCorssPoint;
514  brv=RightValueWithLeftCorssPoint;
515  b2=rightb;
516  blv=lmv;
517  }
518  // find out whether the boundary goes from top-right to bottom-left or
519  // from top-left to bottom-right and then do interpolation
520  if(y>=t2) {
521  return fourpoint(new double[4]{data[idx-1],data[idx],tlv,trv},
522  new double[2]{x,y},new double[4]{C1[idx-1],C1[idx],C1[idx-1],C1[idx]},
523  new double[4]{C2[idx],C2[idx],t2,t2});
524  } else if(y<b2) {
525  return fourpoint(new double[4]{blv,brv,data[idx-N1-1],data[idx-N1]},new double[2]{x,y},new double[4]{C1[idx-1],C1[idx],C1[idx-N1-1],C1[idx-N1]},new double[4]{b2,b2,C2[idx-N1-1],C2[idx-N1]});
526  } else if((rightb-leftb)/(dC1m[idx])*(x-C1[idx-1])+(b2)>y) {
527  return leftb>rightb ? threepoint(new double[3]{trv,blv,brv},new double[2]{x,y},new double[3]{C1[idx],C1[idx-N1],C1[idx-N1]},new double[3]{t2,b2,b2}) :
528  threepoint(new double[3]{tlv,blv,brv},new double[2]{x,y},new double[3]{C1[idx],C1[idx-N1],C1[idx-N1]},new double[3]{t2,b2,b2});
529  } else if((rightb-leftb)/(dC1m[idx])*(x-C1[idx-1])+(b2)<=y) {
530  return leftb>rightb ? threepoint(new double[3]{tlv,blv,brv},new double[2]{x,y},new double[3]{C1[idx-1],C1[idx-1],C1[idx]},new double[3]{t2,b2,b2}) :
531  threepoint(new double[3] {trv,blv,brv},new double[2]{x,y},new double[3]{C1[idx],C1[idx-1],C2[idx]},new double[3]{t2,b2,b2});
532  }
533  }
534  return data[idx];
535  } // 2D
536  // 1D
537  // |<---dC1m[idx]--->|
538  // +---r1---+---r2---+
539  // C1[idx-1] x C1[idx]
540  size_t idx=GetIdxOfPointToTheRightOf(x,0,N1-1);
541  double r2=(C1[idx]-x)/dC1m[idx];
542  double r1=1-r2;
543  double xval=data[idx]*r1+data[idx-1]*r2;
544  if (gDebug>0) Info("GetData","data(x=%.4f)=%.2f, "
545  "C1[%zu]=%.4f, C1[%zu]=%.4f, r1=%.2f, r2=%0.2f",
546  x,xval,idx-1,C1[idx-1],idx,C1[idx],r1,r2);
547  return xval;
548 }
std::vector< double > dC1m
step length to previous point alone C1
Definition: Grid.h:21
double fourpoint(double dataset[4], double tarlocationset[2], double pointxset[4], double pointyset[4]) const
Calculate interpolation value between four point (rectangle)
Definition: Grid.cc:604
std::vector< bool > fIsFixed
true if field values are fixed
Definition: Grid.h:128
double threepoint(double dataset[3], double tarlocationset[2], double pointxset[3], double pointyset[3]) const
Calculate interpolation value between three point (triangle)
Definition: Grid.cc:574
size_t N2
number of points along the 2nd coordinate
Definition: Grid.h:54
size_t N1
number of points along the 1st coordinate
Definition: Grid.h:53
std::vector< double > C2
the 2nd coordinates of the points
Definition: Grid.h:13
double twopoint(double dataset[2], double tarlocationset, double pointxset[2]) const
Calculate interpolation value between two point.
Definition: Grid.cc:561
std::vector< double > dC1p
step length to next point alone C1
Definition: Grid.h:20
size_t GetIdxOfPointToTheRightOf(double c1, size_t begin, size_t end) const
Get index of point near.
Definition: Grid.cc:163
std::vector< double > dC2m
step length to previous point along C2
Definition: Grid.h:23
std::vector< double > dC2p
step length to next point along C2
Definition: Grid.h:22
size_t N3
number of points along the 3rd coordinate
Definition: Grid.h:55
std::vector< double > C1
the 1st coordinates of the points
Definition: Grid.h:12

§ GetE()

double GeFiCa::Grid::GetE ( double  c1,
double  c2 = 0,
double  c3 = 0 
) const
inline

Definition at line 88 of file Grid.h.

89  { return GetData(Et,c1,c2,c3); }
virtual double GetData(const std::vector< double > &data, double c1, double c2, double c3) const
Interpolate grid data at (c1,c2,c3).
Definition: Grid.cc:260
std::vector< double > Et
total electric field strength
Definition: Grid.h:16

§ GetE1()

double GeFiCa::Grid::GetE1 ( double  c1,
double  c2 = 0,
double  c3 = 0 
) const
inline

Definition at line 90 of file Grid.h.

91  { return GetData(E1,c1,c2,c3); }
std::vector< double > E1
projection of Et on C1
Definition: Grid.h:17
virtual double GetData(const std::vector< double > &data, double c1, double c2, double c3) const
Interpolate grid data at (c1,c2,c3).
Definition: Grid.cc:260

§ GetE2()

double GeFiCa::Grid::GetE2 ( double  c1,
double  c2 = 0,
double  c3 = 0 
) const
inline

Definition at line 92 of file Grid.h.

93  { return GetData(E2,c1,c2,c3); }
virtual double GetData(const std::vector< double > &data, double c1, double c2, double c3) const
Interpolate grid data at (c1,c2,c3).
Definition: Grid.cc:260
std::vector< double > E2
projection of Et on C2
Definition: Grid.h:18

§ GetE3()

double GeFiCa::Grid::GetE3 ( double  c1,
double  c2 = 0,
double  c3 = 0 
) const
inline

Definition at line 94 of file Grid.h.

95  { return GetData(E3,c1,c2,c3); }
std::vector< double > E3
projection of Et on C3
Definition: Grid.h:19
virtual double GetData(const std::vector< double > &data, double c1, double c2, double c3) const
Interpolate grid data at (c1,c2,c3).
Definition: Grid.cc:260

§ GetFieldLineFrom()

FieldLine * Grid::GetFieldLineFrom ( double  c1,
double  c2,
double  c3 = 0,
bool  positive = true 
)
virtual

Propogate a field line from (c1,c2,c3).

Examples:
pointContact/drawFields.cc.

Definition at line 629 of file Grid.cc.

630 {
631  FieldLine *line=new FieldLine();
632 
633  int i=0;
634  while (true) { // if (x,y) is in crystal
635  line->C1.push_back(x);
636  line->C2.push_back(y);
637  if (x>=C1[GetN()-1]||x<=C1[0]||y>=C2[GetN()-1]||y<=C2[0]
638  || fIsFixed[GetIdxOfPointToTheRightOf(x,y,0,N2)]) {//out of crystal
639  if (i==0) // initial point is not in crystal
640  Warning("GetFieldLineFrom", "Start point (%.1fcm,%.1fcm)"
641  " is not in crystal! Stop propagating.", x/cm, y/cm);
642  break;
643  }
644 
645  double ex = GetE1(x,y), ey = GetE2(x,y);
646  double et = TMath::Sqrt(ex*ex+ey*ey); // total E
647  if (et==0) {
648  Warning("GetFieldLineFrom", "E@(x=%.1fmm,y=%.1fmm)=%.1fV/cm!",
649  x/mm, y/mm, et/volt*cm);
650  break;
651  }
652  double weight=5/(et/volt*mm); // propagate more in weaker field
653  double dt=10*ns, mu=50000*cm2/volt/sec; // mu is mobility
654  double dx=mu*ex*dt*weight, dy=mu*ey*dt*weight;
655  line->dC1p.push_back(dx);
656  line->dC2p.push_back(dy);
657  if (gDebug>0)
658  Printf("%04d x=%.3fmm, y=%.3fmm, Ex=%.2fV/cm, Ey=%.2fV/cm, "
659  "weight=%.3f, dx=%.3fmm, dy=%.3fmm", i, x/mm, y/mm, ex/volt*cm,
660  ey/volt*cm, weight, dx/mm, dy/mm);
661  if (i>2000) {
662  Info("GetFieldLineFrom", "Propagated more than 2000 steps. Stop");
663  break;
664  }
665  if (positive) { x+=dx; y+=dy; } else { x-=dx; y-=dy; }
666  i++;
667  }
668 
669  return line;
670 }
double GetE1(double c1, double c2=0, double c3=0) const
Definition: Grid.h:90
static const double ns
nano second
Definition: Units.h:7
double GetE2(double c1, double c2=0, double c3=0) const
Definition: Grid.h:92
std::vector< bool > fIsFixed
true if field values are fixed
Definition: Grid.h:128
static const double cm
centimeter
Definition: Units.h:12
size_t N2
number of points along the 2nd coordinate
Definition: Grid.h:54
std::vector< double > C2
the 2nd coordinates of the points
Definition: Grid.h:13
std::vector< double > dC1p
step length to next point alone C1
Definition: Grid.h:20
size_t GetIdxOfPointToTheRightOf(double c1, size_t begin, size_t end) const
Get index of point near.
Definition: Grid.cc:163
size_t GetN()
total number of points
Definition: Grid.h:26
static const double mm
minimeter
Definition: Units.h:15
std::vector< double > dC2p
step length to next point along C2
Definition: Grid.h:22
Electric field line data.
Definition: Grid.h:35
static const double cm2
centimeter squared
Definition: Units.h:13
static const double volt
Definition: Units.h:20
static const double sec
second
Definition: Units.h:10
std::vector< double > C1
the 1st coordinates of the points
Definition: Grid.h:12

§ GetIdxOfMaxV()

size_t Grid::GetIdxOfMaxV ( )
protected

Get index of the grid point with max potential.

Definition at line 73 of file Grid.cc.

74 {
75  double max=Vp[0]; size_t idx=0;
76  for (size_t i=1; i<GetN(); i++) {
77  if (Vp[i]>max) {
78  idx=i;
79  max=Vp[i];
80  }
81  }
82  return idx;
83 }
size_t GetN()
total number of points
Definition: Grid.h:26
std::vector< double > Vp
potential at each point
Definition: Grid.h:15

§ GetIdxOfMinV()

size_t Grid::GetIdxOfMinV ( )
protected

Get index of the grid point with min potential.

Definition at line 86 of file Grid.cc.

87 {
88  double min=Vp[0]; size_t idx=0;
89  for (size_t i=1; i<GetN(); i++) {
90  if (Vp[i]<min) {
91  idx=i;
92  min=Vp[i];
93  }
94  }
95  return idx;
96 }
size_t GetN()
total number of points
Definition: Grid.h:26
std::vector< double > Vp
potential at each point
Definition: Grid.h:15

§ GetIdxOfPointToTheRightOf() [1/3]

size_t Grid::GetIdxOfPointToTheRightOf ( double  c1,
size_t  begin,
size_t  end 
) const
protected

Get index of point near.

Parameters
c1in between
begin&
end.

Definition at line 163 of file Grid.cc.

164 {
165  //if (end==0) end=N1-1;
166  if (begin>=end)return end;
167  size_t mid=(begin+end)/2;
168  if(C1[mid]>=c1)return GetIdxOfPointToTheRightOf(c1,begin,mid);
169  else return GetIdxOfPointToTheRightOf(c1,mid+1,end);
170 }
size_t GetIdxOfPointToTheRightOf(double c1, size_t begin, size_t end) const
Get index of point near.
Definition: Grid.cc:163
std::vector< double > C1
the 1st coordinates of the points
Definition: Grid.h:12

§ GetIdxOfPointToTheRightOf() [2/3]

size_t Grid::GetIdxOfPointToTheRightOf ( double  c1,
double  c2,
size_t  begin,
size_t  end 
) const
protected

Definition at line 173 of file Grid.cc.

175 {
176  //if (end==0) end=N2-1;
177  //search using binary search
178  // if(begin>=end)cout<<"to x"<<begin<<" "<<end<<endl;;
179  if(begin>=end)return GetIdxOfPointToTheRightOf(c1,end*N1,(end+1)*N1-1);
180  size_t mid=((begin+end)/2);
181  if(C2[mid*N1]>=c2) return GetIdxOfPointToTheRightOf(c1,c2,begin,mid);
182  else return GetIdxOfPointToTheRightOf(c1,c2,mid+1,end);
183 }
size_t N1
number of points along the 1st coordinate
Definition: Grid.h:53
std::vector< double > C2
the 2nd coordinates of the points
Definition: Grid.h:13
size_t GetIdxOfPointToTheRightOf(double c1, size_t begin, size_t end) const
Get index of point near.
Definition: Grid.cc:163

§ GetIdxOfPointToTheRightOf() [3/3]

size_t Grid::GetIdxOfPointToTheRightOf ( double  c1,
double  c2,
double  c3,
size_t  begin,
size_t  end 
) const
protected

Definition at line 186 of file Grid.cc.

188 {
189  //if (end==0) end=N3-1;
190  //search using binary search
191  if(begin>=end)return GetIdxOfPointToTheRightOf(c1,c2,begin,begin+N1*N2-1);
192  size_t mid=((begin/(N1*N2)+end/(N1*N2))/2)*N1*N2;
193  if(C3[mid]>=c3)return GetIdxOfPointToTheRightOf(c1,c2,c3,begin,mid);
194  else return GetIdxOfPointToTheRightOf(c1,c2,c3,mid+1,end);
195 }
size_t N2
number of points along the 2nd coordinate
Definition: Grid.h:54
size_t N1
number of points along the 1st coordinate
Definition: Grid.h:53
size_t GetIdxOfPointToTheRightOf(double c1, size_t begin, size_t end) const
Get index of point near.
Definition: Grid.cc:163
std::vector< double > C3
the 3rd coordinates of the points
Definition: Grid.h:14

§ GetTree()

TTree * Grid::GetTree ( bool  createNew = false)

Create &/or return a TTree with field data.

Parameters
[in]createNewis a flag
  • if false (default), the function returns the point of an existing tree, or create one if there is none.
  • if true, the function always create a new tree and delete the old one if it exists.
Examples:
hemispherical/compare2analytic.cc, planar/checkInitialization.cc, planar/compare2analytic.cc, planar/showConvergingSteps.cc, pointContact/checkInitialization.cc, pointContact/drawFields.cc, segmented/drawSliceInPhi.cc, trueCoaxial/checkInitialization.cc, and trueCoaxial/compare2analytic.cc.

Definition at line 214 of file Grid.cc.

215 {
216  if (fTree) { if (createNew) delete fTree; else return fTree; }
217 
218  // define tree
219  bool b,d; double vp,et,e1,e2,e3,c1,c2,c3,p1,p2,p3,m1,m2,m3;
220  fTree = new TTree("t","field data");
221  fTree->SetDirectory(0);
222  fTree->Branch("v",&vp,"v/D");
223  fTree->Branch("e",&et,"e/D");
224  // 1D data
225  fTree->Branch("e1",&e1,"e1/D");
226  fTree->Branch("c1",&c1,"c1/D");
227  fTree->Branch("p1",&p1,"p1/D");
228  fTree->Branch("m1",&m1,"m1/D");
229  if (N2!=0) { // 2D data
230  fTree->Branch("e2",&e2,"e2/D");
231  fTree->Branch("c2",&c2,"c2/D");
232  fTree->Branch("p2",&p2,"p2/D");
233  fTree->Branch("m2",&m2,"m2/D");
234  }
235  if (N3!=0) { // 3D data
236  fTree->Branch("e3",&e3,"e3/D");
237  fTree->Branch("c3",&c3,"c3/D");
238  fTree->Branch("p3",&p3,"p3/D");
239  fTree->Branch("m3",&m3,"m3/D");
240  }
241  fTree->Branch("b",&b,"b/O"); // boundary flag
242  fTree->Branch("d",&d,"d/O"); // depletion flag
243 
244  // fill tree
245  Info("GetTree","%zu entries",GetN());
246  for (size_t i=0; i<GetN(); i++) {
247  e1=E1[i]; c1=C1[i]; p1=dC1p[i]; m1=dC1m[i]; // 1D data
248  if (N2!=0) { e2=E2[i]; c2=C2[i]; p2=dC2p[i]; m2=dC2m[i]; } // 2D data
249  if (N3!=0) { e3=E3[i]; c3=C3[i]; p3=dC3p[i]; m3=dC3m[i]; } // 3D data
250  vp=Vp[i]; et=Et[i]; b=fIsFixed[i]; d=fIsDepleted[i]; // common data
251  fTree->Fill();
252  }
253 
254  fTree->GetListOfBranches()->ls();
255  fTree->ResetBranchAddresses(); // disconnect from local variables
256  return fTree;
257 }
std::vector< double > dC1m
step length to previous point alone C1
Definition: Grid.h:21
std::vector< double > E3
projection of Et on C3
Definition: Grid.h:19
std::vector< bool > fIsFixed
true if field values are fixed
Definition: Grid.h:128
std::vector< double > E1
projection of Et on C1
Definition: Grid.h:17
size_t N2
number of points along the 2nd coordinate
Definition: Grid.h:54
std::vector< double > dC3p
step length to next point alone C3
Definition: Grid.h:24
TTree * fTree
! ROOT tree to visualize fields
Definition: Grid.h:130
std::vector< double > C2
the 2nd coordinates of the points
Definition: Grid.h:13
std::vector< double > dC1p
step length to next point alone C1
Definition: Grid.h:20
size_t GetN()
total number of points
Definition: Grid.h:26
std::vector< double > dC2m
step length to previous point along C2
Definition: Grid.h:23
std::vector< bool > fIsDepleted
true if a grid point is depleted
Definition: Grid.h:129
std::vector< double > dC2p
step length to next point along C2
Definition: Grid.h:22
size_t N3
number of points along the 3rd coordinate
Definition: Grid.h:55
std::vector< double > Vp
potential at each point
Definition: Grid.h:15
std::vector< double > dC3m
step length to previous point alone C3
Definition: Grid.h:25
std::vector< double > Et
total electric field strength
Definition: Grid.h:16
std::vector< double > E2
projection of Et on C2
Definition: Grid.h:18
std::vector< double > C1
the 1st coordinates of the points
Definition: Grid.h:12
std::vector< double > C3
the 3rd coordinates of the points
Definition: Grid.h:14

§ GetV()

double GeFiCa::Grid::GetV ( double  c1,
double  c2 = 0,
double  c3 = 0 
) const
inline

Get potential at (c1,c2,c3) by interpolation.

Examples:
trueCoaxial/verifyCV.cc.

Definition at line 86 of file Grid.h.

87  {return GetData(Vp,c1,c2,c3); }
virtual double GetData(const std::vector< double > &data, double c1, double c2, double c3) const
Interpolate grid data at (c1,c2,c3).
Definition: Grid.cc:260
std::vector< double > Vp
potential at each point
Definition: Grid.h:15

§ IsDepleted()

bool Grid::IsDepleted ( )

Check if every grid point is depleted.

Definition at line 99 of file Grid.cc.

100 {
101  for(size_t i=0;i<GetN();i++) {
102  // calculate one more time in case of adding two fields together,
103  // one is depleted, the other is not
104  OverRelaxAt(i);
105  if (fIsDepleted[i]=false) return false;
106  }
107  return true;
108 }
virtual void OverRelaxAt(size_t idx)
Over relax potential Vp[.
Definition: Grid.h:135
size_t GetN()
total number of points
Definition: Grid.h:26
std::vector< bool > fIsDepleted
true if a grid point is depleted
Definition: Grid.h:129

§ operator*=()

Grid & Grid::operator*= ( double  scale)

Potentials at all points are multiplied by.

Parameters
scale.

Definition at line 66 of file Grid.cc.

67 {
68  for (size_t i=0; i<GetN(); i++) Vp[i]=Vp[i]*scale;
69  return *this;
70 }
size_t GetN()
total number of points
Definition: Grid.h:26
std::vector< double > Vp
potential at each point
Definition: Grid.h:15

§ operator+=()

Grid & Grid::operator+= ( Grid other)

Potentials of this grid are summed with those of.

Parameters
othergrid.

Definition at line 51 of file Grid.cc.

52 {
53  if (GetN()!=other.GetN()) {
54  Warning("+=",
55  "Numbers of points in two grids are different, return *this.");
56  return *this;
57  }
58  for (size_t i=0; i<GetN(); i++) {
59  Vp[i]=Vp[i]+other.Vp[i];
60  Src[i]=Src[i]+other.Src[i];
61  }
62  return *this;
63 }
std::vector< double > Src
-(net impurity concentration)x|Qe|/epsilon
Definition: Grid.h:56
size_t GetN()
total number of points
Definition: Grid.h:26
std::vector< double > Vp
potential at each point
Definition: Grid.h:15

§ OverRelaxAt()

virtual void GeFiCa::Grid::OverRelaxAt ( size_t  idx)
inlineprotectedvirtual

Over relax potential Vp[.

Parameters
idx].

Reimplemented in GeFiCa::XYZ, GeFiCa::R, GeFiCa::Rho, GeFiCa::X, GeFiCa::RhoPhi, GeFiCa::RhoZ, and GeFiCa::RTheta.

Definition at line 135 of file Grid.h.

135 {};

§ SetupWith()

void Grid::SetupWith ( Detector detector)
virtual

Fix potentials on boundaries based on.

Parameters
detectorgeometry. It fills Points data based on
detectorgeometry and N1, N2 and/or N3, and raises the flag fIsFixed for points on/outside electrodes. It has to be called before SuccessiveOverRelax().

Reimplemented in GeFiCa::XYZ, GeFiCa::R, GeFiCa::Rho, GeFiCa::RhoPhi, GeFiCa::RhoZ, GeFiCa::RTheta, and GeFiCa::X.

Definition at line 111 of file Grid.cc.

112 {
113  if (GetN()>0) { // this function can only be called once
114  Warning("SetupWith", "has been called. Do nothing.");
115  return;
116  }
117 }
size_t GetN()
total number of points
Definition: Grid.h:26

§ SolveAnalytically()

void Grid::SolveAnalytically ( )

Solve Poisson's Equation analytically.

It only accepts a constant impurity throughout the grid.

Definition at line 154 of file Grid.cc.

155 {
156  if (Src[0]!=Src[N1-1]) {
157  Error("SolveAnalytically", "can't handle changing impurity.");
158  abort();
159  }
160 }
std::vector< double > Src
-(net impurity concentration)x|Qe|/epsilon
Definition: Grid.h:56
size_t N1
number of points along the 1st coordinate
Definition: Grid.h:53

§ SuccessiveOverRelax()

void Grid::SuccessiveOverRelax ( )

Successively over-relax potentials on grid points.

Examples:
hemispherical/compare2analytic.cc, planar/compare2analytic.cc, planar/optimizeRelaxationFactor.cc, planar/showConvergingSteps.cc, pointContact/calculateFields.cc, pointContact/optimizeRelaxationFactor.cc, segmented/drawSliceInPhi.cc, trueCoaxial/compare2analytic.cc, and trueCoaxial/verifyCV.cc.

Definition at line 120 of file Grid.cc.

121 {
122  if (dC1p.size()<1) {
123  Error("SuccessiveOverRelax", "Grid is not ready. "
124  "Please call SetupWith(Detector&) first.");
125  abort();
126  }
127  Info("SuccessiveOverRelax","Start...");
128  double error=1; Iterations=0;
129  TStopwatch watch; watch.Start();
130  while (Iterations<MaxIterations && error>Tolerance) {
131  if (Iterations%100==0) Printf("%4zu steps, "
132  "error: %.1e (tolerance: %.0e)", Iterations, error, Tolerance);
133  double numerator=0, denominator=0;
134  for (size_t i=0; i<GetN(); i++) {
135  double old=Vp[i]; // save old value of Vp[i]
136  OverRelaxAt(i); // update Vp[i]
137  if(old>0) denominator+=old; else denominator-=old;
138  double diff=Vp[i]-old;
139  if(diff>0) numerator+=(diff); else numerator-=(diff);
140  }
141  if (denominator==0) {
142  Error("SuccessiveOverRelax","Sum of Vs == 0!"); abort();
143  }
144  error = numerator/denominator;
145  Iterations++;
146  }
147  CalculateE();
148  Printf("%4zu steps, error: %.1e (tolerance: %.0e)",
149  Iterations, error, Tolerance);
150  Info("SuccessiveOverRelax", "CPU time: %.1f s", watch.CpuTime());
151 }
size_t Iterations
number of iterations of SOR performed
Definition: Grid.h:60
double Tolerance
SOR stops when error<Tolerance.
Definition: Grid.h:59
virtual void OverRelaxAt(size_t idx)
Over relax potential Vp[.
Definition: Grid.h:135
std::vector< double > dC1p
step length to next point alone C1
Definition: Grid.h:20
size_t GetN()
total number of points
Definition: Grid.h:26
virtual void CalculateE()
Calculate Et, E1, E2, E3 from Vp.
Definition: Grid.cc:551
std::vector< double > Vp
potential at each point
Definition: Grid.h:15

§ threepoint()

double Grid::threepoint ( double  dataset[3],
double  tarlocationset[2],
double  pointxset[3],
double  pointyset[3] 
) const
protected

Calculate interpolation value between three point (triangle)

Definition at line 574 of file Grid.cc.

576 {
577  double x=tarlocationset[0];
578  double x1=pointxset[0];
579  double x2=pointxset[1];
580  double x3=pointxset[2];
581 
582  double y=tarlocationset[1];
583  double y1=pointyset[0];
584  double y2=pointyset[1];
585  double y3=pointyset[2];
586  //x=(1-u-v)*x1+u*x2+v*x3
587  //y=(1-u-v)*y1+u*y2+v*y3
588  double v=((y2-y1)*(x-x1)-(y-y1)*(x2-x1))/((x3-x1*(y2-y1)-(y3-y1)*(x2-x1)));
589  double u=((y3-y1)*(x-x1)-(y-y1)*(x3-x1))/((x2-x1*(y3-y1)-(y2-y1)*(x3-x1)));
590  //double u=(y-y1-v*(y3-y1))/(y2-y1);
591 
592  double result=dataset[0]*(1-u-v)+dataset[1]*u+dataset[2]*v;
593 
594  delete [] dataset;
595  delete [] tarlocationset;
596  delete [] pointxset;
597  delete [] pointyset;
598 
599  if (TMath::Abs(result)<1e-10) return 0;
600  return result;
601 }

§ twopoint()

double Grid::twopoint ( double  dataset[2],
double  tarlocationset,
double  pointxset[2] 
) const
protected

Calculate interpolation value between two point.

Definition at line 561 of file Grid.cc.

563 {
564  double ab=(tarlocationset-pointxset[0])/(pointxset[1]-pointxset[0]);
565 
566  double result=dataset[0]*(1-ab) + dataset[1]*ab;
567  delete [] dataset;
568  delete [] pointxset;
569  if (TMath::Abs(result)<1e-10) return 0;
570  return result;
571 }

Member Data Documentation

§ fDetector

Detector* GeFiCa::Grid::fDetector
protected

! Pointer to associated detector object

Definition at line 131 of file Grid.h.

§ fIsDepleted

std::vector<bool> GeFiCa::Grid::fIsDepleted
protected

true if a grid point is depleted

Definition at line 129 of file Grid.h.

§ fIsFixed

std::vector<bool> GeFiCa::Grid::fIsFixed
protected

true if field values are fixed

Definition at line 128 of file Grid.h.

§ fTree

TTree* GeFiCa::Grid::fTree
protected

! ROOT tree to visualize fields

Definition at line 130 of file Grid.h.

§ Iterations

size_t GeFiCa::Grid::Iterations

number of iterations of SOR performed

Definition at line 60 of file Grid.h.

§ MaxIterations

size_t GeFiCa::Grid::MaxIterations

maximal iterations of SOR to be performed

Examples:
planar/optimizeRelaxationFactor.cc, planar/showConvergingSteps.cc, and pointContact/optimizeRelaxationFactor.cc.

Definition at line 57 of file Grid.h.

§ N1

size_t GeFiCa::Grid::N1

number of points along the 1st coordinate

Definition at line 53 of file Grid.h.

§ N2

size_t GeFiCa::Grid::N2

number of points along the 2nd coordinate

Definition at line 54 of file Grid.h.

§ N3

size_t GeFiCa::Grid::N3

number of points along the 3rd coordinate

Definition at line 55 of file Grid.h.

§ RelaxationFactor

double GeFiCa::Grid::RelaxationFactor

within (0,2), used to speed up convergence

Examples:
planar/optimizeRelaxationFactor.cc, pointContact/calculateFields.cc, pointContact/optimizeRelaxationFactor.cc, and trueCoaxial/verifyCV.cc.

Definition at line 58 of file Grid.h.

§ Src

std::vector<double> GeFiCa::Grid::Src

-(net impurity concentration)x|Qe|/epsilon

Definition at line 56 of file Grid.h.

§ Tolerance

double GeFiCa::Grid::Tolerance

SOR stops when error<Tolerance.

Definition at line 59 of file Grid.h.


The documentation for this class was generated from the following files: