Initializes the magnetic field component.
This method opens the magnetic field map file and triangular mesh.
312 {
313 if (fieldmapname.empty()) {
314 B2ERROR("The filename for the beamline magnetic field component is empty !");
315 return;
316 }
317 std::string l_fieldmapname = FileSystem::findFile("/data/" + fieldmapname);
318
319 if (interpolname.empty()) {
320 B2ERROR("The filename for the triangulation of the beamline magnetic field component is empty !");
321 return;
322 }
323 std::string l_interpolname = FileSystem::findFile("/data/" + interpolname);
324
325 m_rmax = validRadius;
326
327 B2DEBUG(50, "Delaunay triangulation of the beamline field map: " << l_interpolname);
328 B2DEBUG(50, "Beamline field map: " << l_fieldmapname);
329
330
331 ifstream INd(l_interpolname);
332 int nts; INd >> nts;
333 B2DEBUG(50, "Total number of triangles: " << nts);
334 vector<triangle_t> ts;
335 ts.reserve(nts);
336
337 {
338 triangle_t p;
339 while (INd >> nts >> p.j0 >> p.j1 >> p.j2 >> p.n0 >> p.n1 >> p.n2) ts.push_back(p);
340 }
341
342
343 io::filtering_istream IN;
344 IN.push(io::gzip_decompressor());
345 IN.push(io::file_source(l_fieldmapname));
346
347 int nrphi;
348 IN >> nrphi >> m_rj >> m_nr >> m_nphi;
349 IN >> m_nz >> m_zj >> m_dz0 >> m_dz1;
350 m_idphi = m_nphi / M_PI;
351 m_rj2 = m_rj * m_rj;
352 m_idz0 = 1 / m_dz0;
353 m_idz1 = 1 / m_dz1;
354 int nz0 = 2 * int(m_zj / m_dz0);
355 m_zmax = (m_nz - (nz0 + 1)) / 2 * m_dz1 + m_zj;
356 m_nz1 = (m_nz - (nz0 + 1)) / 2;
357 m_nz2 = m_nz1 + nz0;
358
359
360 struct cs_t {double c, s;};
361 vector<cs_t> cs(nrphi);
362 vector<xy_t> pc;
363 vector<ROOT::Math::XYZVector> tbc;
364 pc.reserve(nrphi);
365 char cbuf[256]; IN.getline(cbuf, 256);
366 double rmax = 0;
367 for (int j = 0; j < nrphi; j++) {
368 IN.getline(cbuf, 256);
369 char* next = cbuf;
370 double r = strtod(next, &next);
371 double phi = strtod(next, &next);
372 strtod(next, &next);
373 double Br = strtod(next, &next);
374 double Bphi = strtod(next, &next);
375 double Bz = strtod(next, nullptr);
376 r *= 100;
377 rmax = std::max(r, rmax);
378 if (phi == 0) {
379 cs[j] = { 1, 0};
380 } else if (phi == 180) {
381 cs[j] = { -1, 0};
382 } else {
383 phi *= M_PI / 180;
384 cs[j] = {cos(phi), sin(phi)};
385 }
386 double x = r * cs[j].c, y = r * cs[j].s;
387 pc.push_back({x, y});
388 if (cs[j].s == 0) Bphi = 0;
389 double Bx = Br * cs[j].c - Bphi * cs[j].s;
390 double By = Br * cs[j].s + Bphi * cs[j].c;
391 tbc.emplace_back(Bx, By, Bz);
392 }
393
394
395 m_idr = m_nr / (rmax - m_rj);
396 m_rmax = std::min(m_rmax, rmax);
397 bool reduce = m_rj > m_rmax;
398
399
400
401
402 vector<bool> ip;
403 if (reduce) {
404 ip.resize(nrphi, false);
405 vector<bool> it(ts.size(), false);
406 auto inside = [this](const xy_t & p)->bool{
407 return p.x * p.x + p.y * p.y < m_rmax * m_rmax;
408 };
409
410 for (int i = 0, imax = ts.size(); i < imax; i++) {
411 const triangle_t& p = ts[i];
412 const xy_t& p0 = pc[p.j0], &p1 = pc[p.j1], &p2 = pc[p.j2];
413 if (inside(p0) || inside(p1) || inside(p2)) {
414 it[i] = true;
415 ip[p.j0] = true;
416 ip[p.j1] = true;
417 ip[p.j2] = true;
418 }
419 }
420
421 vector<short int> pindx(nrphi, -1);
422 int rnp = 0;
423 for (int i = 0, imax = ip.size(); i < imax; i++) {
424 if (ip[i]) pindx[i] = rnp++;
425 }
426 vector<xy_t> rpc;
427 rpc.reserve(rnp);
428 for (int i = 0, imax = pc.size(); i < imax; i++) {
429 if (ip[i]) rpc.push_back(pc[i]);
430 }
431
432 vector<short int> tindx(ts.size(), -1);
433 short int rnt = 0;
434 for (int i = 0, imax = it.size(); i < imax; i++) {
435 if (it[i]) tindx[i] = rnt++;
436 }
437 vector<triangle_t> rts;
438 rts.reserve(rnt);
439 short int nt = ts.size();
440 auto newind = [&nt, &tindx, &rnt](short int n) -> short int {return (n < nt) ? tindx[n] : rnt;};
441 for (int i = 0, imax = ts.size(); i < imax; i++) {
442 if (it[i]) {
443 const triangle_t& t = ts[i];
444 rts.push_back({pindx[t.j0], pindx[t.j1], pindx[t.j2], newind(t.n0), newind(t.n1), newind(t.n2)});
445 }
446 }
447
448 B2DEBUG(50, "Reduce map size to cover only region R<" << m_rmax << " cm: Ntriangles=" << rnt << " Nxypoints = " << rnp <<
449 " Nzslices=" << m_nz << " NBpoints = " << rnp * m_nz);
450 std::swap(rpc, pc);
451 std::swap(rts, ts);
452 } else {
453 ip.resize(nrphi, true);
454 }
455 m_nxy = pc.size();
456
457 m_triInterpol.init(pc, ts, 0.1);
458
459 vector<ROOT::Math::XYZVector> bc(m_nxy * m_nz);
460 unsigned int count = 0;
461 for (int i = 0; i < nrphi; i++) {
462 if (ip[i]) bc[count++] = ROOT::Math::XYZVector(tbc[i]);
463 }
464
465 for (int i = 1; i < m_nz; ++i) {
466 for (int j = 0; j < nrphi; j++) {
467 IN.getline(cbuf, 256);
468 if (!ip[j]) continue;
469 char* next = cbuf;
470 next = strchr(next, ' ');
471 next = strchr(next + 1, ' ');
472 next = strchr(next + 1, ' ');
473 double Br = strtod(next, &next);
474 double Bphi = strtod(next, &next);
475 double Bz = strtod(next, nullptr);
476 if (cs[j].s == 0) Bphi = 0;
477 double Bx = Br * cs[j].c - Bphi * cs[j].s;
478 double By = Br * cs[j].s + Bphi * cs[j].c;
479 bc[count++].SetXYZ(Bx, By, Bz);
480 }
481 }
482 assert(count == bc.size());
483 swap(bc, m_B);
484 }
1.13.2