+ ny = TMath::Max(Int_t(fNdiv*TMath::Exp(-(y-fCurrentY)*(y-fCurrentY)/(2*fOrigSigmaY*fOrigSigmaY))),5);
+ ndy = fOrigSigmaY/Double_t(ny);
+ if (ndy>(y2-y-dy)) {
+ ndy =y2-y-dy;
+ if (ndy<kprec) ndy=2*kprec; //calculate new delta y
+ }
+ //
+ Double_t sumch=0;
+ //calculation of x borders and initial step
+ Double_t deltay = (y-y1chev);
+
+ Double_t xp1 = x0+deltay*kx;
+ //x begining of pad at position y
+ Double_t xp2 =xp1+fWidth; //x end of pad at position y
+ Double_t xp3 =xp1+kx*dy; //...at position y+dy
+ Double_t xp4 =xp2+kx*dy; //..
+
+ Double_t x1 = TMath::Min(xp1,xp3);
+ x1 = TMath::Max(xp1,xch-3.5*fOrigSigmaX); //beging of integration
+ Double_t x2 = TMath::Max(xp2,xp4);
+ x2 = TMath::Min(xp2+dy*kx,xch+3.5*fOrigSigmaX); //end of integration
+
+ Int_t nx = TMath::Max(Int_t(fNdiv*TMath::Exp(-(x1-xch)*(x1-xch)/(2*fOrigSigmaX*fOrigSigmaX))*
+ TMath::Exp(-(y1-fCurrentY)*(y1-fCurrentY)/(2*fOrigSigmaY*fOrigSigmaY))),2);
+ Double_t dx = TMath::Min(fOrigSigmaX/Double_t(nx),x2-x1)/5.; //on the border more iteration
+ Double_t ndx=dx;
+
+ if (x2>(x1+kprec)) {
+ for (Double_t x = x1; x<x2+kprec ;){
+ //new step SIZE
+ nx = TMath::Max(Int_t(fNdiv*TMath::Exp(-(x-xch)*(x-xch)/(2*fOrigSigmaX*fOrigSigmaX))),3);
+ ndx = fOrigSigmaX/Double_t(nx);
+ if (ndx>(x2-x-dx)) {
+ ndx =x2-x-dx;
+ }
+ if ( ( (x+dx+ndx)<TMath::Max(xp3,xp1)) || ( (x+dx+ndx)>TMath::Min(xp4,xp2))) {
+ ndx/=5.;
+ }
+ if (ndx<kprec) ndx=2*kprec;
+ //INTEGRAL APROXIMATION
+ Double_t ddx,ddy,dddx,dddy;
+ ddx = xch-(x+dx/2.);
+ ddy = fCurrentY-(y+dy/2.);
+ dddx = cos*ddx-sin*ddy;
+ dddy = sin*ddx+cos*ddy;
+ Double_t z0=fGRF->Eval(dddx,dddy); //middle point
+
+ ddx = xch-(x+dx/2.);
+ ddy = fCurrentY-(y);
+ dddx = cos*ddx-sin*ddy;
+ dddy = sin*ddx+cos*ddy;
+ Double_t z1=fGRF->Eval(dddx,dddy); //point down
+
+ ddx = xch-(x+dx/2.);
+ ddy = fCurrentY-(y+dy);
+ dddx = cos*ddx-sin*ddy;
+ dddy = sin*ddx+cos*ddy;
+ Double_t z3=fGRF->Eval(dddx,dddy); //point up
+
+ ddx = xch-(x);
+ ddy = fCurrentY-(y+dy/2.);
+ dddx = cos*ddx-sin*ddy;
+ dddy = sin*ddx+cos*ddy;
+ Double_t z2=fGRF->Eval(dddx,dddy); //point left
+
+ ddx = xch-(x+dx);
+ ddy = fCurrentY-(y+dy/2.);
+ dddx = cos*ddx-sin*ddy;
+ dddy = sin*ddx+cos*ddy;
+ Double_t z4=fGRF->Eval(dddx,dddy); //point right
+
+
+ if (z0<0) {z0=0;z1=0;z2=0;z3=0;z4=0;}
+
+ Double_t f2x= (z3+z1-2*z0)*4.;//second derivation in y
+ Double_t f2y= (z2+z4-2*z0)*4.;//second derivation in x
+ Double_t f1y= (z3-z1);
+ Double_t z ;
+ z = (z0+f2x/6.+f2y/6.);//second order aproxiation of integral
+ if (kx>kprec){ //positive derivation
+ if (x<(xp1+dy*kx)){ //calculate volume at left border
+ Double_t xx1 = x;
+ Double_t xx2 = TMath::Min(x+dx,xp1+dy*kx);
+ Double_t yy1 = y+(xx1-xp1)/kx;
+ Double_t yy2 = TMath::Min(y+(xx2-xp1)/kx,y+dy);
+ z=z0;
+ if (yy2<y+dy) {
+ z-= z0*(y+dy-yy2)/dy; //constant part rectangle
+ z-= f1y*(xx2-xx1)*(y+dy-yy2)*(y+dy-yy2)/(2.*dx*dy);
+ }
+ z-=z0*(xx2-xx1)*(yy2-yy1)/(2*dx*dy); //constant part rectangle
+
+ }
+ if (x>xp2){ //calculate volume at right border
+ Double_t xx1 = x;
+ Double_t xx2 = x+dx;
+ Double_t yy1 = y+(xx1-xp2)/kx;
+ Double_t yy2 = y+(xx2-xp2)/kx;
+ z=z0;
+ //rectangle part
+ z-=z0*(yy1-y)/dy; //constant part
+ z-=f1y*(xx2-xx1)*(yy1-y)*(yy1-y)/(2*dx*dy);
+ //triangle part
+ z-=z0*(xx2-xx1)*(yy2-yy1)/(2*dx*dy); //constant part
+ }
+ }
+ if (kx<-kprec){ //negative derivation
+ if (x<(xp1+dy*kx)){ //calculate volume at left border
+ Double_t xx1 = x;
+ Double_t xx2 = TMath::Min(x+dx,xp3-dy/kx);
+ Double_t yy1 = y+(xx1-xp1)/kx;
+ Double_t yy2 = TMath::Max(y,yy1+(xx2-xx1)/kx); //yy2<yy1
+ z = z0;
+ z-= z0*(yy2-y)/dy; // constant part rectangle
+ z-= f1y*(xx2-xx1)*(yy2-y)*(yy2-y)/(2.*dx*dy);
+ z-=z0*(xx2-xx1)*(yy1-yy2)/(2*dx*dy); //constant part triangle