5 #include <TStopwatch.h> 15 Warning(
"GetGraph",
"No data for graph. Return 0!");
18 if (fGl) {
return fGl; }
19 fGl =
new TGraph(
GetN(),
C1.data(),
C2.data());
20 fGl->SetName(GetName()); fGl->SetTitle(GetTitle());
21 fGl->SetMarkerStyle(8); fGl->SetMarkerSize(0.4); fGl->SetMarkerColor(12);
27 N1(n1), N2(n2), N3(n3), MaxIterations(5000), RelaxationFactor(1.95),
28 Tolerance(1e-7*
volt), Iterations(0), fTree(0), fDetector(0)
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);
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);
55 "Numbers of points in two grids are different, return *this.");
58 for (
size_t i=0; i<
GetN(); i++) {
68 for (
size_t i=0; i<
GetN(); i++)
Vp[i]=
Vp[i]*scale;
75 double max=
Vp[0];
size_t idx=0;
76 for (
size_t i=1; i<
GetN(); i++) {
88 double min=
Vp[0];
size_t idx=0;
89 for (
size_t i=1; i<
GetN(); i++) {
101 for(
size_t i=0;i<
GetN();i++) {
114 Warning(
"SetupWith",
"has been called. Do nothing.");
123 Error(
"SuccessiveOverRelax",
"Grid is not ready. " 124 "Please call SetupWith(Detector&) first.");
127 Info(
"SuccessiveOverRelax",
"Start...");
129 TStopwatch watch; watch.Start();
130 while (Iterations<MaxIterations && error>
Tolerance) {
132 "error: %.1e (tolerance: %.0e)",
Iterations, error, Tolerance);
133 double numerator=0, denominator=0;
134 for (
size_t i=0; i<
GetN(); 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);
141 if (denominator==0) {
142 Error(
"SuccessiveOverRelax",
"Sum of Vs == 0!"); abort();
144 error = numerator/denominator;
148 Printf(
"%4zu steps, error: %.1e (tolerance: %.0e)",
150 Info(
"SuccessiveOverRelax",
"CPU time: %.1f s", watch.CpuTime());
157 Error(
"SolveAnalytically",
"can't handle changing impurity.");
166 if (begin>=end)
return end;
167 size_t mid=(begin+end)/2;
174 size_t begin,
size_t end)
const 180 size_t mid=((begin+end)/2);
187 size_t begin,
size_t end)
const 200 Info(
"GetC",
"Start...");
202 std::vector<double>original=
Src;
203 for (
size_t i=0; i<
GetN(); i++) {
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");
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");
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");
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");
241 fTree->Branch(
"b",&b,
"b/O");
242 fTree->Branch(
"d",&d,
"d/O");
245 Info(
"GetTree",
"%zu entries",
GetN());
246 for (
size_t i=0; i<
GetN(); i++) {
254 fTree->GetListOfBranches()->ls();
255 fTree->ResetBranchAddresses();
261 double x,
double y,
double z)
const 264 Warning(
"GetData",
"not yet implemented. Abort."); abort();
284 if (idx%
N1!=0) tl=
true;
285 if (idx>=
N1) br=
true;
290 }
else if(tl&&!br&&!bl) {
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) {
294 return twopoint(
new double[2]{data[idx-
N1],data[idx]},
295 y,
new double[2]{
C2[idx-
N1],
C2[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]});
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;
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]});
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]});
323 else if((
C2[idx]-yb)/(xb-
C1[idx-1]*(x-
C1[idx-1])+yb>y))
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});
327 else if((
C2[idx]-yb)/(xb-
C1[idx-1]*(x-
C1[idx-1])+yb<y))
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});
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;
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]});
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]});
352 else if((
C2[idx]-yb)/(
C1[idx]-x*(x-
C1[idx])+yb>y))
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});
356 else if((
C2[idx]-yb)/(
C1[idx]-xb*(x-
C1[idx])+yb<y))
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});
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;
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]});
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});
381 else if((
C2[idx]-yb)/(xb-
C1[idx-1]*(x-
C1[idx-1])+
C2[idx-
N1]>y))
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]});
385 else if((
C2[idx]-yb)/(xb-
C1[idx-1]*(x-
C1[idx-1])+
C2[idx-
N1]<y))
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});
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];
401 double mmv=(
C2[idx]-yb)/(
C2[idx]-
C2[idx-
N1])*bmv+(yb-
C2[idx-
N1])/(
C2[idx]-
C2[idx-
N1])*tmv;
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]});
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});
410 else if((yb-
C2[idx-
N1])/(
C1[idx]-xb*(x-xb)+
C2[idx-
N1]>y))
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});
414 else if((yb-
C2[idx-
N1])/(
C1[idx]-xb*(x-xb)+
C2[idx-
N1]<=y))
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});
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;
437 TopLeftV=TopValueWithBottomCorssPoint;
443 BottomRightV=BottomValueWithTopCorssPoint;
446 TopRightV=TopValueWithBottomCorssPoint;
452 BottomLeftV=BottomValueWithTopCorssPoint;
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]});
493 double lmv=
fIsFixed[idx-1]?data[idx-1]:data[idx-
N1-1];
494 double rmv=
fIsFixed[idx]?data[idx]:data[idx-
N1];
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;
506 trv=RightValueWithLeftCorssPoint;
508 blv=LeftValueWithRightCorssPoint;
513 tlv=LeftValueWithRightCorssPoint;
514 brv=RightValueWithLeftCorssPoint;
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});
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});
541 double r2=(
C1[idx]-x)/
dC1m[idx];
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);
553 for (
size_t i=0; i<
GetN(); i++) {
562 double pointxset[2])
const 564 double ab=(tarlocationset-pointxset[0])/(pointxset[1]-pointxset[0]);
566 double result=dataset[0]*(1-ab) + dataset[1]*ab;
569 if (TMath::Abs(result)<1e-10)
return 0;
575 double pointxset[3],
double pointyset[3])
const 577 double x=tarlocationset[0];
578 double x1=pointxset[0];
579 double x2=pointxset[1];
580 double x3=pointxset[2];
582 double y=tarlocationset[1];
583 double y1=pointyset[0];
584 double y2=pointyset[1];
585 double y3=pointyset[2];
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)));
592 double result=dataset[0]*(1-u-v)+dataset[1]*u+dataset[2]*v;
595 delete [] tarlocationset;
599 if (TMath::Abs(result)<1e-10)
return 0;
605 double pointxset[4],
double pointyset[4])
const 612 double ab=(pointxset[1]-tarlocationset[0])/(pointxset[1]-pointxset[0]);
614 double ba=(tarlocationset[1]-pointyset[3])/(pointyset[1]-pointyset[3]);
617 double result=(dataset[1]*aa+dataset[0]*ab)*ba+(dataset[2]*ab+dataset[3]*aa)*bb;
619 delete [] tarlocationset;
623 if (TMath::Abs(result)<1e-10)
return 0;
635 line->
C1.push_back(x);
636 line->
C2.push_back(y);
640 Warning(
"GetFieldLineFrom",
"Start point (%.1fcm,%.1fcm)" 641 " is not in crystal! Stop propagating.", x/
cm, y/
cm);
646 double et = TMath::Sqrt(ex*ex+ey*ey);
648 Warning(
"GetFieldLineFrom",
"E@(x=%.1fmm,y=%.1fmm)=%.1fV/cm!",
652 double weight=5/(et/
volt*
mm);
654 double dx=mu*ex*dt*weight, dy=mu*ey*dt*weight;
655 line->
dC1p.push_back(dx);
656 line->
dC2p.push_back(dy);
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,
662 Info(
"GetFieldLineFrom",
"Propagated more than 2000 steps. Stop");
665 if (positive) { x+=dx; y+=dy; }
else { x-=dx; y-=dy; }
679 line->
C1.push_back(x);
680 line->
C2.push_back(y);
681 line->
C3.push_back(z);
686 Warning(
"GetFieldLineFrom",
"Start point (%.1fcm,%.1fcm)" 687 " is not in crystal! Stop propagating.", x/
cm, y/
cm);
694 double et = TMath::Sqrt(ex*ex+ey*ey+ez*ez);
696 Warning(
"GetFieldLineFrom",
"E@(x=%.1fmm,y=%.1fmm,z=%.1f)=%.1fV/cm!",
700 double weight=5/(et/
volt*
mm);
702 double dx=mu*ex*dt*weight, dy=mu*ey*dt*weight, dz=mu*ez*dt*weight;
703 line->
dC1p.push_back(dx);
704 line->
dC2p.push_back(dy);
705 line->
dC3p.push_back(dz);
707 Printf(
"%04d x=%.3fmm, y=%.3fmm, z=%.3fmm, Ex=%.2fV/cm, Ey=%.2fV/cm, " 708 "Ez=%.2fV/cm, weight=%.3f, dx=%.3fmm, dy=%.3fmm", i, x/
mm, y/
mm,
711 Info(
"GetFieldLineFrom",
"Propagated more than 2000 steps. Stop");
714 if (positive) { x+=dx; y+=dy; z+=dz; }
else { x-=dx; y-=dy; z-=dz;}
std::vector< double > dC1m
step length to previous point alone C1
size_t Iterations
number of iterations of SOR performed
double GetE1(double c1, double c2=0, double c3=0) const
Grid & operator*=(double scale)
Potentials at all points are multiplied by.
std::vector< double > E3
projection of Et on C3
size_t GetIdxOfMinV()
Get index of the grid point with min potential.
double fourpoint(double dataset[4], double tarlocationset[2], double pointxset[4], double pointyset[4]) const
Calculate interpolation value between four point (rectangle)
Grid & operator+=(Grid &other)
Potentials of this grid are summed with those of.
Data structure of a electric field grid.
static const double ns
nano second
double GetE2(double c1, double c2=0, double c3=0) const
Grid(size_t n1=0, size_t n2=0, size_t n3=0)
Default constructor.
std::vector< bool > fIsFixed
true if field values are fixed
double threepoint(double dataset[3], double tarlocationset[2], double pointxset[3], double pointyset[3]) const
Calculate interpolation value between three point (triangle)
double Tolerance
SOR stops when error<Tolerance.
static const double cm
centimeter
std::vector< double > E1
projection of Et on C1
size_t N2
number of points along the 2nd coordinate
std::vector< double > dC3p
step length to next point alone C3
virtual FieldLine * GetFieldLineFrom(double c1, double c2, bool positive=true)
Propogate a field line from (c1,c2,c3).
std::vector< double > Src
-(net impurity concentration)x|Qe|/epsilon
TTree * fTree
! ROOT tree to visualize fields
size_t N1
number of points along the 1st coordinate
std::vector< double > C2
the 2nd coordinates of the points
void SolveAnalytically()
Solve Poisson's Equation analytically.
virtual void SetupWith(Detector &detector)
Fix potentials on boundaries based on.
virtual void OverRelaxAt(size_t idx)
Over relax potential Vp[.
double twopoint(double dataset[2], double tarlocationset, double pointxset[2]) const
Calculate interpolation value between two point.
Detector & crystal properties.
std::vector< double > dC1p
step length to next point alone C1
A group of discrete points.
size_t GetIdxOfPointToTheRightOf(double c1, size_t begin, size_t end) const
Get index of point near.
virtual double GetData(const std::vector< double > &data, double c1, double c2, double c3) const
Interpolate grid data at (c1,c2,c3).
size_t GetN()
total number of points
std::vector< double > dC2m
step length to previous point along C2
void SuccessiveOverRelax()
Successively over-relax potentials on grid points.
double GetE3(double c1, double c2=0, double c3=0) const
static const double mm
minimeter
std::vector< bool > fIsDepleted
true if a grid point is depleted
std::vector< double > dC2p
step length to next point along C2
A file defining commonly used units & constants.
Electric field line data.
size_t N3
number of points along the 3rd coordinate
static const double cm2
centimeter squared
TTree * GetTree(bool createNew=false)
Create &/or return a TTree with field data.
size_t GetIdxOfMaxV()
Get index of the grid point with max potential.
virtual void CalculateE()
Calculate Et, E1, E2, E3 from Vp.
double GetC()
Get detector capacitance.
std::vector< double > Vp
potential at each point
static const double sec
second
std::vector< double > dC3m
step length to previous point alone C3
The only namespace in GeFiCa.
bool IsDepleted()
Check if every grid point is depleted.
std::vector< double > Et
total electric field strength
std::vector< double > E2
projection of Et on C2
std::vector< double > C1
the 1st coordinates of the points
std::vector< double > C3
the 3rd coordinates of the points