TPC/AliTPCEfield.cxx

   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
  16
  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
  23
  24 #include "TTreeStream.h"
  25 #include "TMath.h"
  26 #include "TLinearFitter.h"
  27 #include "TRandom.h"
  28 #include "AliTPCEfield.h"
  29
  30 /// \cond CLASSIMP
  31 ClassImp(AliTPCEfield)
  32 /// \endcond
  33
  34 AliTPCEfield* AliTPCEfield::fgInstance=0;
  35
  36 AliTPCEfield::AliTPCEfield():
  37   TNamed(),
  38   fScale(0),
  39   fMaxFreq(0),
  40   fIs2D(kTRUE),
  41   fWorkspace(0)   // file with trees, pictures ...
  42 {
  43   ///
  44
  45   for (Int_t i=0; i<3; i++){
  46     fMin[i]=0; fMax[i]=0;
  47   }
  48   fgInstance =this;
  49 }
  50
  51 AliTPCEfield::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 {
  59   ///
  60
  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
  69 void AliTPCEfield::MakeFitFunctions(Int_t maxFreq){
  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
  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
 124 AliTPCEfield::~AliTPCEfield() {
 125   /// Destructor
 126
 127   if (fWorkspace) delete fWorkspace;
 128 }
 129
 130 void AliTPCEfield::SetRange(Double_t x0, Double_t x1, Double_t y0, Double_t y1, Double_t z0,Double_t z1){
 131   /// Set the ranges - coordinates are rescaled in order to use proper
 132   /// cos,sin expansion in scaled space
 133
 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
 142 void 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){
 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
 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
 169 TTree * AliTPCEfield::GetTree(const char * tname){
 170   ///
 171
 172   return ((*fWorkspace)<<tname).GetTree();
 173 }
 174
 175 Double_t AliTPCEfield::Field(Int_t ftype,  Double_t ifx, Double_t ify, Double_t ifz, Double_t x, Double_t y, Double_t z){
 176   /// Field component in
 177   /// f frequency
 178
 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
 213 Double_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){
 214  /// Field component in
 215  /// f frequency
 216
 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
 260 Double_t AliTPCEfield::EvalField(Int_t ifun, Double_t x, Double_t y, Double_t z, Int_t type){
 261   /// Evaluate function ifun at position gx amd gy
 262   /// type == 0 - field
 263   ///      == 1 - Ex
 264   ///      == 2 - Ey
 265   ///      == 3 - Ez
 266
 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
 278 Double_t AliTPCEfield::Eval(Double_t x, Double_t y, Double_t z, Int_t type){
 279   /// Evaluate function ifun at position gx amd gy
 280   /// type == 0 - field
 281   ///      == 1 - Ex
 282   ///      == 2 - Ey
 283   ///      == 3 - Ez
 284
 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
 298 Double_t AliTPCEfield::EvalS(Double_t x, Double_t y, Double_t z,  Int_t type){
 299   /// static evaluation - possible to use it in the TF1
 300
 301   return fgInstance->Eval(x,y,z,type);
 302 }
 303
 304 void AliTPCEfield::FitField(){
 305   /// Fit the e field
 306   /// Minimize chi2 residuals at the boundary points
 307   /// ?Tempoary sollution - integrals can be calculated analytically -
 308
 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
 365 TMatrixD* AliTPCEfield::MakeCorrelation(TMatrixD &matrix){
 366   ///
 367
 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
 380 void AliTPCEfield::DumpField(Double_t gridSize, Double_t step){
 381   ///
 382
 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
 446 void MakeTPC2DExample(AliTPCEfield *field){
 447   //
 448   /*
 449     .L  $ALICE_ROOT/TPC/AliTPCEfield.cxx++
 450     AliTPCEfield *field =  new AliTPCEfield("field",20, kTRUE,kTRUE);
 451     MakeTPC2DExample(field)
 452     field->FitField()
 453     sqrt(field->fFitter.GetChisquare()/field->fFitter.GetNpoints())
 454
 455     TF2 f2("f2","AliTPCEfield::EvalS(x,y,0,0)",90,245,0,250);
 456     f2->SetNpx(100);     f2->SetNpy(100); f2->Draw("colz");
 457
 458     TF2 f2x("f2x","AliTPCEfield::EvalS(x,y,0,1)",90,240,0,240);
 459     f2x->SetNpx(100);     f2x->SetNpy(100);  f2x->Draw("surf2");
 460
 461     TF2 f2y("f2y","AliTPCEfield::EvalS(x,y,0,2)",90,240,0,240);
 462     f2y->SetNpx(100);     f2y->SetNpy(100);  f2y->Draw("surf2");
 463
 464     field->MakeCorrelation(*(field->fFitCovar)).Print()
 465     Double_t index[100000];
 466     for (Int_t i=0; i<field->fFitCovar->GetNrows(); i++) index[i]=i;
 467     TGraph gr(field->fFitCovar->GetNrows(), index, field->fFitParam->GetMatrixArray());
 468     gr->Draw("alp");
 469
 470     field->GetTree()->Draw("AliTPCEfield::EvalS(x,y,0,0):v");
 471
 472
 473     TF2 f2xdy("f2xdy","AliTPCEfield::EvalS(x,y,0,1)/AliTPCEfield::EvalS(x,y,0,2)",90,240,0,240);
 474     f2xdy->SetNpx(100);     f2xdy->SetNpy(100);  f2xdy->Draw("colz");
 475
 476   */
 477
 478   Double_t p0[4];
 479   Double_t p1[4];
 480   Double_t xmin=85, xmax=245;
 481   Double_t ymin=0, ymax=250, deltaY=0.1*ymax;
 482   Double_t vup=1;
 483   Int_t npoints=1000;
 484   field->SetRange(xmin, xmax,ymin,ymax,0,0);
 485   // upper part
 486   p0[0]=xmin; p0[1]=ymax+deltaY; p0[2]=0; p0[3]=vup;
 487   p1[0]=xmax; p1[1]=ymax-deltaY; p1[2]=0; p1[3]=vup;
 488   field->AddBoundaryLine(p0[0],p0[1], p0[2],p0[3],p1[0],p1[1], p1[2],p1[3],1,npoints);
 489   //left
 490   p0[0]=xmin; p0[1]=ymin+deltaY; p0[2]=0; p0[3]=0;
 491   p1[0]=xmin; p1[1]=ymax+deltaY; p1[2]=0; p1[3]=vup;
 492   field->AddBoundaryLine(p0[0],p0[1], p0[2],p0[3],p1[0],p1[1], p1[2],p1[3],2,npoints);
 493   //right
 494   p0[0]=xmax; p0[1]=ymin-deltaY; p0[2]=0; p0[3]=0;
 495   p1[0]=xmax; p1[1]=ymax-deltaY; p1[2]=0; p1[3]=vup;
 496   field->AddBoundaryLine(p0[0],p0[1], p0[2],p0[3],p1[0],p1[1], p1[2],p1[3],3,npoints);
 497   //
 498   //ROC
 499   p0[0]=xmin; p0[1]=deltaY; p0[2]=0; p0[3]=-0;
 500   p1[0]=xmax; p1[1]=-deltaY; p1[2]=0; p1[3]=0;
 501   field->AddBoundaryLine(p0[0],p0[1], p0[2],p0[3],p1[0],p1[1], p1[2],p1[3],4,npoints);
 502 }
 503
 504
 505