]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCEfield.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
doxy: TPC/DA
[u/mrichter/AliRoot.git] / TPC / AliTPCEfield.cxx
CommitLineData
b1f0a2a5 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
b1f0a2a5 16
0fd1fc0b 17/// \class AliTPCEfield
18///
19/// Solution of Laplace equation in cartesian system, with boundary condition.
20///
21/// See details:
22/// http://web.mit.edu/6.013_book/www/chapter5/5.10.html
b1f0a2a5 23
24#include "TTreeStream.h"
25#include "TMath.h"
26#include "TLinearFitter.h"
27#include "TRandom.h"
28#include "AliTPCEfield.h"
29
0fd1fc0b 30/// \cond CLASSIMP
b1f0a2a5 31ClassImp(AliTPCEfield)
0fd1fc0b 32/// \endcond
b1f0a2a5 33
34AliTPCEfield* AliTPCEfield::fgInstance=0;
35
36AliTPCEfield::AliTPCEfield():
37 TNamed(),
38 fScale(0),
39 fMaxFreq(0),
40 fIs2D(kTRUE),
41 fWorkspace(0) // file with trees, pictures ...
42{
3a4edebe 43 ///
0fd1fc0b 44
b1f0a2a5 45 for (Int_t i=0; i<3; i++){
46 fMin[i]=0; fMax[i]=0;
47 }
48 fgInstance =this;
49}
50
51AliTPCEfield::AliTPCEfield(const char* name, Int_t maxFreq, Bool_t is2D, Bool_t useLinear):
52 TNamed(name,name),
53 fScale(0),
54 fMaxFreq(maxFreq),
55 fIs2D(is2D),
56 fUseLinear(useLinear),
57 fWorkspace(0) // file with trees, pictures ...
58{
3a4edebe 59 ///
0fd1fc0b 60
b1f0a2a5 61 for (Int_t i=0; i<3; i++){
62 fMin[i]=0; fMax[i]=0;
63 }
64 fWorkspace=new TTreeSRedirector(Form("%s.root",name));
65 MakeFitFunctions(maxFreq);
66 fgInstance =this;
67}
68
69void AliTPCEfield::MakeFitFunctions(Int_t maxFreq){
0fd1fc0b 70 /// fit functions = \f$ f(x,y,z) = fx(x)*fy(y)*fz(z) \f$
71 /// function can be of following types:
72 /// 0 - constant
73 /// 1 - linear
74 /// 2 - hypx
75 /// 3 - hypy
76 /// 4 - hypz
77
b1f0a2a5 78 Int_t nfunctions=0;
79 if (fIs2D) nfunctions = 1+(maxFreq)*8;
80 if (!fIs2D) nfunctions = 1+(maxFreq)*8;
81 if (fUseLinear) nfunctions+=2;
82 fFitFunctions = new TMatrixD(nfunctions,4);
83 fFitParam = new TVectorD(nfunctions);
84 fFitCovar = new TMatrixD(nfunctions,nfunctions);
85 TMatrixD &fitF = *fFitFunctions;
86 // constant function
87 Int_t counter=1;
88 fitF(0,0)=0;
89 //
90 // linear functions in one cordinate - constan in others
91 if (fUseLinear)
92 for (Int_t ifx=0; ifx<=1; ifx++)
93 for (Int_t ify=0; ify<=1; ify++)
94 for (Int_t ifz=0; ifz<=1; ifz++){
95 if (fIs2D && ifz>0) continue;
96 if ((ifx+ify+ifz)==0) continue;
97 if ((ifx+ify+ifz)>1) continue;
98 fitF(counter,0)= 1; // function type
99 fitF(counter,1)= ifx; // type of x function
100 fitF(counter,2)= ify; // type of y function
101 fitF(counter,3)= ifz; // type of z function
102 counter++;
103 }
104 //
105 if (fIs2D){
106 for (Int_t ihyp=0; ihyp<2; ihyp++)
107 for (Int_t ifx=-maxFreq; ifx<=maxFreq; ifx++)
108 for (Int_t ify=-maxFreq; ify<=maxFreq; ify++){
109 if (TMath::Abs(ify)!=TMath::Abs(ifx)) continue;
110 if (ifx==0) continue;
111 if (ify==0) continue;
112 fitF(counter,0)= 2+ihyp; // function type
113 fitF(counter,1)= ifx; // type of x function - + sinus - cosin
114 fitF(counter,2)= ify; // type of y function
115 fitF(counter,3)= 0; // type of y function
116 counter++;
117 }
118 }
119
120}
121
122
123
124AliTPCEfield::~AliTPCEfield() {
0fd1fc0b 125 /// Destructor
126
b1f0a2a5 127 if (fWorkspace) delete fWorkspace;
128}
129
130void AliTPCEfield::SetRange(Double_t x0, Double_t x1, Double_t y0, Double_t y1, Double_t z0,Double_t z1){
0fd1fc0b 131 /// Set the ranges - coordinates are rescaled in order to use proper
132 /// cos,sin expansion in scaled space
133
b1f0a2a5 134 fMin[0]=x0; fMax[0]=x1;
135 fMin[1]=y0; fMax[1]=y1;
136 fMin[2]=z0; fMax[2]=z1;
137 if (fIs2D) fScale=0.5*(TMath::Abs(x1-x0)+TMath::Abs(y1-y0));
138 if (!fIs2D) fScale=0.5*(TMath::Abs(x1-x0)+TMath::Abs(y1-y0)+TMath::Abs(z1-z0));
139}
140
141
142void AliTPCEfield::AddBoundaryLine(Double_t x0,Double_t y0,Double_t z0, Double_t v0, Double_t x1, Double_t y1, Double_t z1,Double_t v1, Int_t id, Int_t npoints){
0fd1fc0b 143 /// Add a e field boundary line
144 /// From point (x0,y0) to point (x1,y1)
145 /// Linear decrease of potential is assumed
146 /// Boundary can be identified using boundary ID
147 /// The line is written into tree Boundary
148
b1f0a2a5 149 Double_t deltaX = (x1-x0);
150 Double_t deltaY = (y1-y0);
151 Double_t deltaZ = (z1-z0);
152 Double_t deltaV = (v1-v0);
153 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
154 Double_t bpoint=gRandom->Rndm();
155 Double_t x = x0+deltaX*bpoint;
156 Double_t y = y0+deltaY*bpoint;
157 Double_t z = z0+deltaZ*bpoint;
158 Double_t v = v0+deltaV*bpoint;
159 (*fWorkspace)<<"Boundary"<<
160 "x="<<x<< // x coordinate
161 "y="<<y<< // y coordinate
162 "z="<<z<< // z coordinate
163 "v="<<v<< // potential
164 "id="<<id<< // boundary ID
165 "\n";
166 }
167}
168
169TTree * AliTPCEfield::GetTree(const char * tname){
3a4edebe 170 ///
0fd1fc0b 171
b1f0a2a5 172 return ((*fWorkspace)<<tname).GetTree();
173}
174
175Double_t AliTPCEfield::Field(Int_t ftype, Double_t ifx, Double_t ify, Double_t ifz, Double_t x, Double_t y, Double_t z){
0fd1fc0b 176 /// Field component in
177 /// f frequency
178
b1f0a2a5 179 Double_t fx=1,fy=1,fz=1;
180 const Double_t kEps=0.01;
181 //
182 if (ftype==0) return 1;
183 if (ftype==1) {
184 if (TMath::Nint(ifx)==1) return x;
185 if (TMath::Nint(ify)==1) return y;
186 if (TMath::Nint(ifz)==1) return z;
187 }
188 Double_t pi = TMath::Pi();
189 if (ifx>kEps) fx = (ftype==2) ? SinHNorm(ifx*pi*x,TMath::Abs(ifx*pi)): TMath::Sin(ifx*pi*x);
190 if (ifx<-kEps) fx = (ftype==2) ? CosHNorm(ifx*pi*x,TMath::Abs(ifx*pi)): TMath::Cos(ifx*pi*x);
191 //
192 if (ify>kEps) fy = (ftype==3) ? SinHNorm(ify*pi*y,TMath::Abs(ify*pi)): TMath::Sin(ify*pi*y);
193 if (ify<-kEps) fy = (ftype==3) ? CosHNorm(ify*pi*y,TMath::Abs(ify*pi)): TMath::Cos(ify*pi*y);
194 //
195 if (ifz>kEps) fz = (ftype==4) ? SinHNorm(ifz*pi*z,TMath::Abs(ifz*pi)): TMath::Sin(ifz*pi*z);
196 if (ifz<-kEps) fz = (ftype==4) ? CosHNorm(ifz*pi*z,TMath::Abs(ifz*pi)): TMath::Cos(ifz*pi*z);
197 (*fWorkspace)<<"eval"<<
198 "x="<<x<<
199 "y="<<y<<
200 "z="<<z<<
201 "ifx="<<ifx<<
202 "ify="<<ify<<
203 "ifz="<<ifz<<
204 "fx="<<fx<<
205 "fy="<<fy<<
206 "fz="<<fz<<
207 "ftype="<<ftype<<
208 "\n";
209 return fx*fy*fz;
210}
211
212
213Double_t AliTPCEfield::FieldDn(Int_t ftype, Double_t ifx, Double_t ify, Double_t ifz, Int_t dn, Double_t x, Double_t y, Double_t z){
0fd1fc0b 214 /// Field component in
215 /// f frequency
216
b1f0a2a5 217 Double_t fx=1,fy=1,fz=1;
218 const Double_t kEps=0.01;
219 //
220 if (ftype==0) return 0.;
221 if (ftype==1) {
222 Double_t value=0;
223 if (TMath::Nint(ifx)==1 &&dn==0) value=1.;
224 if (TMath::Nint(ify)==1 &&dn==1) value=1.;
225 if (TMath::Nint(ifz)==1 &&dn==2) value=1.;
226 return value;
227 }
228 Double_t pi = TMath::Pi();
229 if (ifx>kEps) fx = (ftype==2) ? SinHNorm(ifx*pi*x,TMath::Abs(ifx*pi)): TMath::Sin(ifx*pi*x);
230 if (ifx<-kEps) fx = (ftype==2) ? CosHNorm(ifx*pi*x,TMath::Abs(ifx*pi)): TMath::Cos(ifx*pi*x);
231 //
232 if (ify>kEps) fy = (ftype==3) ? SinHNorm(ify*pi*y,TMath::Abs(ify*pi)): TMath::Sin(ify*pi*y);
233 if (ify<-kEps) fy = (ftype==3) ? CosHNorm(ify*pi*y,TMath::Abs(ify*pi)): TMath::Cos(ify*pi*y);
234 //
235 if (ifz>kEps) fz = (ftype==4) ? SinHNorm(ifz*pi*z,TMath::Abs(ifz*pi)): TMath::Sin(ifz*pi*z);
236 if (ifz<-kEps) fz = (ftype==4) ? CosHNorm(ifz*pi*z,-TMath::Abs(ifz*pi)): TMath::Cos(ifz*pi*z);
237
238 if (dn==0){
239 if (ifx>kEps) fx = (ftype==2) ? CosHNorm(ifx*pi*x,TMath::Abs(ifx*pi)): TMath::Cos(ifx*pi*x);
240 if (ifx<-kEps) fx = (ftype==2) ? SinHNorm(ifx*pi*x,TMath::Abs(ifx*pi)): -TMath::Sin(ifx*pi*x);
241 fx*=ifx*pi;
242 }
243 if (dn==1){
244 if (ify>kEps) fy = (ftype==3) ? CosHNorm(ify*pi*y,TMath::Abs(ify*pi)): TMath::Cos(ify*pi*y);
245 if (ify<-kEps) fy = (ftype==3) ? SinHNorm(ify*pi*y,TMath::Abs(ify*pi)): -TMath::Sin(ify*pi*y);
246 fy*=ify*pi;
247 }
248 if (dn==2){
249 if (ifz>kEps) fz = (ftype==4) ? CosHNorm(ifz*pi*z,TMath::Abs(ifz*pi)): TMath::Cos(ifz*pi*z);
250 if (ifz<-kEps) fz = (ftype==4) ? SinHNorm(ifz*pi*z,TMath::Abs(ifz*pi)): -TMath::Sin(ifz*pi*z);
251 fz*=ifz*pi;
252 }
253
254 return fx*fy*fz;
255}
256
257
258
259
260Double_t AliTPCEfield::EvalField(Int_t ifun, Double_t x, Double_t y, Double_t z, Int_t type){
0fd1fc0b 261 /// Evaluate function ifun at position gx amd gy
262 /// type == 0 - field
263 /// == 1 - Ex
264 /// == 2 - Ey
265 /// == 3 - Ez
266
b1f0a2a5 267 TMatrixD &mat = *fFitFunctions;
268 Int_t fid = TMath::Nint(mat(ifun,0));
269 Double_t ifx = (mat(ifun,1));
270 Double_t ify = (mat(ifun,2));
271 Double_t ifz = (mat(ifun,3));
272 //
273 if (type==0) return Field(fid,ifx,ify,ifz, x, y,z);
274 if (type>0) return FieldDn(fid,ifx,ify,ifz,type-1, x, y,z);
275 return 0;
276}
277
278Double_t AliTPCEfield::Eval(Double_t x, Double_t y, Double_t z, Int_t type){
0fd1fc0b 279 /// Evaluate function ifun at position gx amd gy
280 /// type == 0 - field
281 /// == 1 - Ex
282 /// == 2 - Ey
283 /// == 3 - Ez
284
b1f0a2a5 285 Double_t value=0;
286 Double_t lx= 2.*(x-(fMin[0]+fMax[0])*0.5)/fScale;
287 Double_t ly= 2.*(y-(fMin[1]+fMax[1])*0.5)/fScale;
288 Double_t lz= 2.*(z-(fMin[2]+fMax[2])*0.5)/fScale;
289 //
290 Int_t nfun=fFitFunctions->GetNrows();
291 for (Int_t ifun=0; ifun<nfun; ifun++){
292 if (type==0) value+=(*fFitParam)[ifun]*EvalField(ifun,lx,ly,lz,type);
293 if (type>0) value+=2*(*fFitParam)[ifun]*EvalField(ifun,lx,ly,lz,type)/fScale;
294 }
295 return value;
296}
297
298Double_t AliTPCEfield::EvalS(Double_t x, Double_t y, Double_t z, Int_t type){
0fd1fc0b 299 /// static evaluation - possible to use it in the TF1
300
b1f0a2a5 301 return fgInstance->Eval(x,y,z,type);
302}
303
304void AliTPCEfield::FitField(){
0fd1fc0b 305 /// Fit the e field
306 /// Minimize chi2 residuals at the boundary points
307 /// ?Tempoary sollution - integrals can be calculated analytically -
308
b1f0a2a5 309 Int_t nfun=fFitFunctions->GetNrows();
310 Double_t *fun =new Double_t[nfun];
311 fFitter= new TLinearFitter(nfun, Form("hyp%d", nfun-1));
312 //
313 TTree * tree = GetTree("Boundary");
314 Int_t npoints = tree->GetEntries();
315 Int_t *indexes = new Int_t[npoints];
316 Double_t *rindex = new Double_t[npoints];
317 //
318
319 Double_t x=0, y=0, z=0, v=0;
320 tree->SetBranchAddress("x",&x);
321 tree->SetBranchAddress("y",&y);
322 tree->SetBranchAddress("z",&z);
323 tree->SetBranchAddress("v",&v);
324 TMatrixD valMatrix(npoints,4);
325 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
326 tree->GetEntry(ipoint);
327 valMatrix(ipoint,0)=2.*(x-(fMin[0]+fMax[0])*0.5)/fScale;
328 valMatrix(ipoint,1)=2.*(y-(fMin[1]+fMax[1])*0.5)/fScale;
329 valMatrix(ipoint,2)=2.*(z-(fMin[2]+fMax[2])*0.5)/fScale;
330 valMatrix(ipoint,3)=v;
331 rindex[ipoint]=gRandom->Rndm();
332 }
333 TMath::Sort(npoints,rindex,indexes);
334 for (Int_t jpoint=0; jpoint<npoints; jpoint++){
335 Int_t ipoint=indexes[jpoint];
336 //
337 Double_t lx= valMatrix(ipoint,0);
338 Double_t ly= valMatrix(ipoint,1);
339 Double_t lz= valMatrix(ipoint,2);
340 Double_t value = valMatrix(ipoint,3);
341 for (Int_t ifun=1; ifun<nfun; ifun++){
342 Double_t ffun=EvalField(ifun,lx,ly,lz,0);
343 fun[ifun-1]=ffun;
344 }
345 fFitter->AddPoint(fun,value,1);
346 }
347 fFitter->Eval();
348 fFitter->GetCovarianceMatrix(*fFitCovar);
349 fFitter->GetParameters(*fFitParam);
350
351 (*fWorkspace)<<"fitGlobal"<<
352 "covar.="<<fFitCovar<< // covariance matrix
353 "param.="<<fFitParam<< // fit parameters
354 "fun.="<<fFitFunctions<< // type of fit functions
355 "\n";
356
357 for (Int_t ifun=1; ifun<nfun; ifun++){
358
359 }
360
361
362}
363
364
365TMatrixD* AliTPCEfield::MakeCorrelation(TMatrixD &matrix){
3a4edebe 366 ///
0fd1fc0b 367
b1f0a2a5 368 Int_t nrows = matrix.GetNrows();
369 TMatrixD * mat = new TMatrixD(nrows,nrows);
370 for (Int_t irow=0; irow<nrows; irow++)
371 for (Int_t icol=0; icol<nrows; icol++){
372 (*mat)(irow,icol)= matrix(irow,icol)/TMath::Sqrt(matrix(irow,irow)*matrix(icol,icol));
373 }
374 return mat;
375}
376
377
378
379
380void AliTPCEfield::DumpField(Double_t gridSize, Double_t step){
3a4edebe 381 ///
0fd1fc0b 382
b1f0a2a5 383 Double_t stepSize=0.001*fScale/fMaxFreq;
384 //
385 for (Double_t x = fMin[0]+stepSize; x<=fMax[0]-stepSize; x+=gridSize){
386 for (Double_t y = fMin[1]+stepSize; y<=fMax[1]-stepSize; y+=gridSize)
387 for (Double_t z = fMin[2]; z<=fMax[2]; z+=gridSize){
388 //
389 //
390 Double_t v = Eval(x,y,z,0);
391 Double_t ex = Eval(x,y,z,1);
392 Double_t ey = Eval(x,y,z,2);
393 Double_t ez = Eval(x,y,z,3);
394 Double_t dexdx = (Eval(x,y,z,1)-Eval(x-stepSize,y,z,1))/stepSize; // numerical derivative
395 Double_t deydy = (Eval(x,y,z,2)-Eval(x,y-stepSize,z,2))/stepSize;
396 Double_t dezdz = (Eval(x,y,z,3)-Eval(x,y,z-stepSize,3))/stepSize;
397 (*fWorkspace)<<"dumpField"<<
398 "x="<<x<< // position
399 "y="<<y<<
400 "z="<<z<<
401 "v="<<v<< // potential
402 "ex="<<ex<< // Efield
403 "ey="<<ey<<
404 "ez="<<ez<<
405 "dexdx="<<dexdx<< //gradient of e field
406 "deydy="<<deydy<<
407 "dezdz="<<dezdz<<
408 "\n";
409 }
410 }
411 //
412 //
413 for (Double_t x = fMin[0]+stepSize; x<=fMax[0]-stepSize; x+=gridSize){
414 Double_t sumEx=0;
415 Double_t sumEy=0;
416 Double_t sumExEy=0;
417 for (Double_t y = fMin[1]+0.2; y<=fMax[1]-stepSize; y+=step)
418 for (Double_t z = fMin[2]; z<=fMax[2]; z+=gridSize){
419 //
420 //
421 Double_t v = Eval(x,y+step*0.5,z,0);
422 Double_t ex = Eval(x,y+step*0.5,z,1);
423 Double_t ey = Eval(x,y+step*0.5,z,2);
424 Double_t ez = Eval(x,y+step*0.5,z,3);
425 sumEx+=ex*step;
426 sumEy+=ey*step;
427 sumExEy+=(ex/ey)*step;
428 (*fWorkspace)<<"dumpDistortion"<<
429 "x="<<x<< // position
430 "y="<<y<<
431 "z="<<z<<
432 "v="<<v<< // potential
433 "ex="<<ex<< // Efield
434 "ey="<<ey<<
435 "ez="<<ez<<
436 //
437 "sEx="<<sumEx<< // field x integral
438 "sEy="<<sumEy<< // field y integral
439 "sExEy="<<sumExEy<< // tan integral
440 "\n";
441 }
442 }
443}
444
445
446void MakeTPC2DExample(AliTPCEfield *field){
3a4edebe 447 ///
448
b1f0a2a5 449 /*
450 .L $ALICE_ROOT/TPC/AliTPCEfield.cxx++
451 AliTPCEfield *field = new AliTPCEfield("field",20, kTRUE,kTRUE);
452 MakeTPC2DExample(field)
453 field->FitField()
454 sqrt(field->fFitter.GetChisquare()/field->fFitter.GetNpoints())
455
456 TF2 f2("f2","AliTPCEfield::EvalS(x,y,0,0)",90,245,0,250);
457 f2->SetNpx(100); f2->SetNpy(100); f2->Draw("colz");
458
459 TF2 f2x("f2x","AliTPCEfield::EvalS(x,y,0,1)",90,240,0,240);
460 f2x->SetNpx(100); f2x->SetNpy(100); f2x->Draw("surf2");
461
462 TF2 f2y("f2y","AliTPCEfield::EvalS(x,y,0,2)",90,240,0,240);
463 f2y->SetNpx(100); f2y->SetNpy(100); f2y->Draw("surf2");
464
465 field->MakeCorrelation(*(field->fFitCovar)).Print()
466 Double_t index[100000];
467 for (Int_t i=0; i<field->fFitCovar->GetNrows(); i++) index[i]=i;
468 TGraph gr(field->fFitCovar->GetNrows(), index, field->fFitParam->GetMatrixArray());
469 gr->Draw("alp");
470
471 field->GetTree()->Draw("AliTPCEfield::EvalS(x,y,0,0):v");
472
473
474 TF2 f2xdy("f2xdy","AliTPCEfield::EvalS(x,y,0,1)/AliTPCEfield::EvalS(x,y,0,2)",90,240,0,240);
475 f2xdy->SetNpx(100); f2xdy->SetNpy(100); f2xdy->Draw("colz");
476
477 */
478
479 Double_t p0[4];
480 Double_t p1[4];
481 Double_t xmin=85, xmax=245;
482 Double_t ymin=0, ymax=250, deltaY=0.1*ymax;
483 Double_t vup=1;
484 Int_t npoints=1000;
485 field->SetRange(xmin, xmax,ymin,ymax,0,0);
486 // upper part
487 p0[0]=xmin; p0[1]=ymax+deltaY; p0[2]=0; p0[3]=vup;
488 p1[0]=xmax; p1[1]=ymax-deltaY; p1[2]=0; p1[3]=vup;
489 field->AddBoundaryLine(p0[0],p0[1], p0[2],p0[3],p1[0],p1[1], p1[2],p1[3],1,npoints);
490 //left
491 p0[0]=xmin; p0[1]=ymin+deltaY; p0[2]=0; p0[3]=0;
492 p1[0]=xmin; p1[1]=ymax+deltaY; p1[2]=0; p1[3]=vup;
493 field->AddBoundaryLine(p0[0],p0[1], p0[2],p0[3],p1[0],p1[1], p1[2],p1[3],2,npoints);
494 //right
495 p0[0]=xmax; p0[1]=ymin-deltaY; p0[2]=0; p0[3]=0;
496 p1[0]=xmax; p1[1]=ymax-deltaY; p1[2]=0; p1[3]=vup;
497 field->AddBoundaryLine(p0[0],p0[1], p0[2],p0[3],p1[0],p1[1], p1[2],p1[3],3,npoints);
498 //
499 //ROC
500 p0[0]=xmin; p0[1]=deltaY; p0[2]=0; p0[3]=-0;
501 p1[0]=xmax; p1[1]=-deltaY; p1[2]=0; p1[3]=0;
502 field->AddBoundaryLine(p0[0],p0[1], p0[2],p0[3],p1[0],p1[1], p1[2],p1[3],4,npoints);
503}
504
505
506