305 {
306 map<int, G4ThreeVector> v;
307
308 double wh = h * h, wH = H * H, wnorm = wh + wH;
309 double tn = ((b - a) / 2 / h * wh + (B - A) / 2 / H * wH) / wnorm, d = h * tn, D = H * tn;
310
311 double m = (a + b) * 0.5, M = (A + B) * 0.5;
312
313 const double eps = 0.5e-3;
314 if (fabs(a - (m - d)) > eps || fabs(b - (m + d)) > eps || fabs(A - (M - D)) > eps || fabs(B - (M + D)) > eps) {
315 double alfa =
atan(tn) * 180 / M_PI;
316 B2WARNING("Cannot make parallel sides better than 0.5 mcm: alpha =" << alpha << " alpha from sides = " << alfa << " da = " << a -
317 (m - d) << " db = " << b - (m + d) << " dA = " << A - (M - D) << " dB = " << B - (M + D));
318 }
319
320 v[1] = G4ThreeVector(-(m - d) / 2, h / 2, -150);
321 v[2] = G4ThreeVector(-(m + d) / 2, -h / 2, -150);
322 v[3] = G4ThreeVector((m + d) / 2, -h / 2, -150);
323 v[4] = G4ThreeVector((m - d) / 2, h / 2, -150);
324 v[5] = G4ThreeVector(-(M - D) / 2, H / 2, 150);
325 v[6] = G4ThreeVector(-(M + D) / 2, -H / 2, 150);
326 v[7] = G4ThreeVector((M + D) / 2, -H / 2, 150);
327 v[8] = G4ThreeVector((M - D) / 2, H / 2, 150);
328
329 if (wrapthick != 0) {
330 map<int, G4ThreeVector> nv;
331 nv[1] = newvertex(wrapthick, v[1], v[5], v[2], v[4]);
332 nv[2] = newvertex(wrapthick, v[2], v[6], v[3], v[1]);
333 nv[3] = newvertex(wrapthick, v[3], v[7], v[4], v[2]);
334 nv[4] = newvertex(wrapthick, v[4], v[8], v[1], v[3]);
335 nv[5] = newvertex(wrapthick, v[5], v[1], v[8], v[6]);
336 nv[6] = newvertex(wrapthick, v[6], v[2], v[5], v[7]);
337 nv[7] = newvertex(wrapthick, v[7], v[3], v[6], v[8]);
338 nv[8] = newvertex(wrapthick, v[8], v[4], v[7], v[5]);
339 std::swap(nv, v);
340 }
341
342
343
344 return v;
345 }
double atan(double a)
atan for double