GeFiCa
Germanium detector Field Calculator
Grid.cc
Go to the documentation of this file.
1 #include <TMath.h>
2 #include <TTree.h>
3 #include <TStyle.h>
4 #include <TGraph.h>
5 #include <TStopwatch.h>
6 
7 #include "Grid.h"
8 #include "Units.h"
9 #include "Detector.h"
10 using namespace GeFiCa;
11 
13 {
14  if (GetN()<1) {
15  Warning("GetGraph", "No data for graph. Return 0!");
16  return 0;
17  }
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);
22  return fGl;
23 }
24 //______________________________________________________________________________
25 //
26 Grid::Grid(size_t n1, size_t n2, size_t n3) : 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 }
49 //______________________________________________________________________________
50 //
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 }
64 //______________________________________________________________________________
65 //
66 Grid& Grid::operator*=(double scale)
67 {
68  for (size_t i=0; i<GetN(); i++) Vp[i]=Vp[i]*scale;
69  return *this;
70 }
71 //______________________________________________________________________________
72 //
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 }
84 //______________________________________________________________________________
85 //
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 }
97 //______________________________________________________________________________
98 //
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 }
109 //______________________________________________________________________________
110 //
111 void Grid::SetupWith(Detector &detector)
112 {
113  if (GetN()>0) { // this function can only be called once
114  Warning("SetupWith", "has been called. Do nothing.");
115  return;
116  }
117 }
118 //______________________________________________________________________________
119 //
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 }
152 //______________________________________________________________________________
153 //
155 {
156  if (Src[0]!=Src[N1-1]) {
157  Error("SolveAnalytically", "can't handle changing impurity.");
158  abort();
159  }
160 }
161 //______________________________________________________________________________
162 //
163 size_t Grid::GetIdxOfPointToTheRightOf(double c1,size_t begin,size_t end) const
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 }
171 //______________________________________________________________________________
172 //
173 size_t Grid::GetIdxOfPointToTheRightOf(double c1,double c2,
174  size_t begin,size_t end) const
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 }
184 //______________________________________________________________________________
185 //
186 size_t Grid::GetIdxOfPointToTheRightOf(double c1, double c2,double c3,
187  size_t begin,size_t end) const
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 }
196 //______________________________________________________________________________
197 //
198 double Grid::GetC()
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 }
212 //______________________________________________________________________________
213 //
214 TTree* Grid::GetTree(bool createNew)
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 }
258 //______________________________________________________________________________
259 //
260 double Grid::GetData(const std::vector<double> &data,
261  double x, double y, double z) const
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 }
549 //______________________________________________________________________________
550 //
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 }
559 //______________________________________________________________________________
560 //
561 double Grid::twopoint(double dataset[2],double tarlocationset,
562  double pointxset[2])const
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 }
572 //______________________________________________________________________________
573 //
574 double Grid::threepoint(double dataset[3],double tarlocationset[2],
575  double pointxset[3],double pointyset[3])const
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 }
602 //______________________________________________________________________________
603 //
604 double Grid::fourpoint(double dataset[4],double tarlocationset[2],
605  double pointxset[4],double pointyset[4])const
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 }
626 
627 //______________________________________________________________________________
628 //
629 FieldLine* Grid::GetFieldLineFrom(double x, double y, bool positive)
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 }
671 //______________________________________________________________________________
672 //
673 FieldLine* Grid::GetFieldLineFrom(double x, double y, double z, bool positive)
674 {
675  FieldLine *line=new FieldLine();
676 
677  int i=0;
678  while (true) { // if (x,y,z) is in crystal
679  line->C1.push_back(x);
680  line->C2.push_back(y);
681  line->C3.push_back(z);
682  if (x>=C1[GetN()-1]||x<=C1[0]||y>=C2[GetN()-1]||y<=C2[0]
683  || z>=C3[GetN()-1]||z<=C3[0]
684  || fIsFixed[GetIdxOfPointToTheRightOf(x,y,z,N2)]) {//out of crystal
685  if (i==0) // initial point is not in crystal
686  Warning("GetFieldLineFrom", "Start point (%.1fcm,%.1fcm)"
687  " is not in crystal! Stop propagating.", x/cm, y/cm);
688  break;
689  }
690 
691  Printf("here");
692  double ex = GetE1(x,y,z), ey = GetE2(x,y,z), ez = GetE3(x,y,z);
693  Printf("here");
694  double et = TMath::Sqrt(ex*ex+ey*ey+ez*ez); // total E
695  if (et==0) {
696  Warning("GetFieldLineFrom", "E@(x=%.1fmm,y=%.1fmm,z=%.1f)=%.1fV/cm!",
697  x/mm, y/mm, z/mm, et/volt*cm);
698  break;
699  }
700  double weight=5/(et/volt*mm); // propagate more in weaker field
701  double dt=10*ns, mu=50000*cm2/volt/sec; // mu is mobility
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);
706  if (gDebug>0)
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,
709  z/mm, ex/volt*cm, ey/volt*cm, ez/volt*cm, weight, dx/mm, dy/mm);
710  if (i>2000) {
711  Info("GetFieldLineFrom", "Propagated more than 2000 steps. Stop");
712  break;
713  }
714  if (positive) { x+=dx; y+=dy; z+=dz; } else { x-=dx; y-=dy; z-=dz;}
715  i++;
716  }
717 
718  return line;
719 }
std::vector< double > dC1m
step length to previous point alone C1
Definition: Grid.h:21
size_t Iterations
number of iterations of SOR performed
Definition: Grid.h:60
double GetE1(double c1, double c2=0, double c3=0) const
Definition: Grid.h:90
Grid & operator*=(double scale)
Potentials at all points are multiplied by.
Definition: Grid.cc:66
std::vector< double > E3
projection of Et on C3
Definition: Grid.h:19
TGraph * GetGraph()
Definition: Grid.cc:12
size_t GetIdxOfMinV()
Get index of the grid point with min potential.
Definition: Grid.cc:86
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
Grid & operator+=(Grid &other)
Potentials of this grid are summed with those of.
Definition: Grid.cc:51
Data structure of a electric field grid.
Definition: Grid.h:50
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
Grid(size_t n1=0, size_t n2=0, size_t n3=0)
Default constructor.
Definition: Grid.cc:26
std::vector< bool > fIsFixed
true if field values are fixed
Definition: Grid.h:129
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
double Tolerance
SOR stops when error<Tolerance.
Definition: Grid.h:59
static const double cm
centimeter
Definition: Units.h:12
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
virtual FieldLine * GetFieldLineFrom(double c1, double c2, bool positive=true)
Propogate a field line from (c1,c2,c3).
Definition: Grid.cc:629
std::vector< double > Src
-(net impurity concentration)x|Qe|/epsilon
Definition: Grid.h:56
TTree * fTree
! ROOT tree to visualize fields
Definition: Grid.h:131
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
void SolveAnalytically()
Solve Poisson&#39;s Equation analytically.
Definition: Grid.cc:154
virtual void SetupWith(Detector &detector)
Fix potentials on boundaries based on.
Definition: Grid.cc:111
virtual void OverRelaxAt(size_t idx)
Over relax potential Vp[.
Definition: Grid.h:136
double twopoint(double dataset[2], double tarlocationset, double pointxset[2]) const
Calculate interpolation value between two point.
Definition: Grid.cc:561
Detector & crystal properties.
Definition: Detector.h:32
std::vector< double > dC1p
step length to next point alone C1
Definition: Grid.h:20
A group of discrete points.
Definition: Grid.h:9
size_t GetIdxOfPointToTheRightOf(double c1, size_t begin, size_t end) const
Get index of point near.
Definition: Grid.cc:163
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
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
void SuccessiveOverRelax()
Successively over-relax potentials on grid points.
Definition: Grid.cc:120
double GetE3(double c1, double c2=0, double c3=0) const
Definition: Grid.h:94
static const double mm
minimeter
Definition: Units.h:15
std::vector< bool > fIsDepleted
true if a grid point is depleted
Definition: Grid.h:130
std::vector< double > dC2p
step length to next point along C2
Definition: Grid.h:22
A file defining commonly used units & constants.
Electric field line data.
Definition: Grid.h:35
size_t N3
number of points along the 3rd coordinate
Definition: Grid.h:55
static const double cm2
centimeter squared
Definition: Units.h:13
static const double volt
Definition: Units.h:20
TTree * GetTree(bool createNew=false)
Create &/or return a TTree with field data.
Definition: Grid.cc:214
size_t GetIdxOfMaxV()
Get index of the grid point with max potential.
Definition: Grid.cc:73
virtual void CalculateE()
Calculate Et, E1, E2, E3 from Vp.
Definition: Grid.cc:551
double GetC()
Get detector capacitance.
Definition: Grid.cc:198
std::vector< double > Vp
potential at each point
Definition: Grid.h:15
static const double sec
second
Definition: Units.h:10
std::vector< double > dC3m
step length to previous point alone C3
Definition: Grid.h:25
The only namespace in GeFiCa.
Definition: Detector.h:6
bool IsDepleted()
Check if every grid point is depleted.
Definition: Grid.cc:99
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