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