REST-for-Physics  v2.3
Rare Event Searches ToolKit for Physics
TRestAxionMagneticField.cxx
1 /******************** REST disclaimer ***********************************
2  * This file is part of the REST software framework. *
3  * *
4  * Copyright (C) 2016 GIFNA/TREX (University of Zaragoza) *
5  * For more information see http://gifna.unizar.es/trex *
6  * *
7  * REST is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * REST is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have a copy of the GNU General Public License along with *
18  * REST in $REST_PATH/LICENSE. *
19  * If not, see http://www.gnu.org/licenses/. *
20  * For the list of contributors see $REST_PATH/CREDITS. *
21  *************************************************************************/
22 
218 
219 #include "TRestAxionMagneticField.h"
220 
221 using namespace std;
222 
223 #include "TRestPhysics.h"
224 using namespace REST_Physics;
225 
226 #include "TGraph.h"
227 ClassImp(TRestAxionMagneticField);
228 
233 
248 TRestAxionMagneticField::TRestAxionMagneticField(const char* cfgFileName, string name)
249  : TRestMetadata(cfgFileName) {
250  RESTDebug << "Entering TRestAxionMagneticField constructor( cfgFileName, name )" << RESTendl;
251 
252  Initialize();
253 
255 
257 }
258 
263  RESTDebug << "Entering ... TRestAxionMagneticField() destructor." << RESTendl;
264 }
265 
270  SetSectionName(this->ClassName());
271  SetLibraryVersion(LIBRARY_VERSION);
272 
273  fCanvas = NULL;
274  fHisto = NULL;
275 }
276 
321 TCanvas* TRestAxionMagneticField::DrawHistogram(TString projection, TString Bcomp, Int_t volIndex,
322  Double_t step, TString style, Double_t depth) {
323  Double_t step_x, step_y, step_z;
325 
326  if (fCanvas != NULL) {
327  delete fCanvas;
328  fCanvas = NULL;
329  }
330 
331  if (fHisto != NULL) {
332  delete fHisto;
333  fHisto = NULL;
334  }
335 
336  if (volIndex < 0) volIndex = 0;
337 
338  if (volIndex >= GetNumberOfVolumes()) {
339  RESTError << volIndex << " corresponds to none volume index " << RESTendl;
340  RESTError << "Total number of volumes : " << GetNumberOfVolumes() << RESTendl;
341  RESTError << "Setting volIndex to the first volume" << RESTendl;
342  volIndex = 0;
343  }
344 
345  if (step <= 0) {
346  step_x = fMeshSize[volIndex].X();
347  step_y = fMeshSize[volIndex].Y();
348  step_z = fMeshSize[volIndex].Z();
349  } else
350  step_x = step_y = step_z = step;
351 
352  MagneticFieldVolume* vol = GetMagneticVolume(volIndex);
353  if (!vol) return fCanvas;
354 
355  if (!(projection == "XY" || projection == "XZ" || projection == "YZ")) {
356  RESTError << "You entered : " << projection << " as a projection but you have to choose XY, XZ or YZ"
357  << RESTendl;
358  return fCanvas;
359  }
360 
361  Double_t centerX = fPositions[volIndex][0];
362  Double_t centerY = fPositions[volIndex][1];
363  Double_t centerZ = fPositions[volIndex][2];
364 
365  Double_t halfSizeX = vol->mesh.GetNetSizeX() / 2.;
366  Double_t halfSizeY = vol->mesh.GetNetSizeY() / 2.;
367  Double_t halfSizeZ = vol->mesh.GetNetSizeZ() / 2.;
368 
369  Double_t xMin = centerX - halfSizeX;
370  Double_t yMin = centerY - halfSizeY;
371  Double_t zMin = centerZ - halfSizeZ;
372 
373  Double_t xMax = centerX + halfSizeX;
374  Double_t yMax = centerY + halfSizeY;
375  Double_t zMax = centerZ + halfSizeZ;
376 
377  Int_t nBinsX = (xMax - xMin) / step_x;
378  Int_t nBinsY = (yMax - yMin) / step_y;
379  Int_t nBinsZ = (zMax - zMin) / step_z;
380 
381  Double_t x, y, z;
382  Double_t B;
383  TVector3 Bvec;
384 
385  if (projection == "XY") {
386  fCanvas = new TCanvas("fCanvas", "");
387  fHisto = new TH2D("", "", nBinsX, xMin, xMax, nBinsY, yMin, yMax);
388 
389  if (depth < -100000.0)
390  z = (zMin + zMax) / 2.0;
391  else if ((depth >= zMin) && (depth <= zMax))
392  z = depth;
393  else
394  RESTError << "You entered depth = " << depth << ", but you have to choose depth between " << zMin
395  << " and " << zMax << RESTendl;
396  x = xMin;
397 
398  for (Int_t i = 0; i < nBinsX; i++) {
399  y = yMin;
400  for (Int_t j = 0; j < nBinsY; j++) {
401  Bvec = GetMagneticField(TVector3(x, y, z), false);
402  if (Bcomp == "X")
403  B = Bvec[0];
404  else {
405  if (Bcomp == "Y")
406  B = Bvec[1];
407  else {
408  if (Bcomp == "Z")
409  B = Bvec[2];
410  else
411  RESTError << "You entered : " << Bcomp
412  << " as a B component but you have to choose X, Y or Z" << RESTendl;
413  }
414  }
415  fHisto->Fill(x, y, B);
416  y = y + step_y;
417  }
418  x = x + step_x;
419  }
420 
421  fCanvas->cd();
422  fHisto->SetBit(TH1::kNoStats);
423  fHisto->GetXaxis()->SetTitle("x (mm)");
424  fHisto->GetYaxis()->SetTitle("y (mm)");
425 
426  if (Bcomp == "X") {
427  TString title = Form("B_{x} against x and y @ z = %.1f", z);
428  fHisto->SetTitle(title);
429  fCanvas->SetTitle(title);
430  }
431  if (Bcomp == "Y") {
432  TString title = Form("B_{y} against x and y @ z = %.1f", z);
433  fHisto->SetTitle(title);
434  fCanvas->SetTitle(title);
435  }
436  if (Bcomp == "Z") {
437  TString title = Form("B_{z} against x and y @ z = %.1f", z);
438  fHisto->SetTitle(title);
439  fCanvas->SetTitle(title);
440  }
441 
442  fHisto->Draw(style);
443  return fCanvas;
444  }
445 
446  if (projection == "XZ") {
447  TCanvas* fCanvas = new TCanvas("fCanvas", "");
448  fHisto = new TH2D("", "", nBinsX, xMin, xMax, nBinsZ, zMin, zMax);
449 
450  if (depth < -100000.0)
451  y = (yMin + yMax) / 2.0;
452  else if ((depth >= yMin) && (depth <= yMax))
453  y = depth;
454  else
455  RESTError << "You entered depth = " << depth << ", but you have to choose depth between " << yMin
456  << " and " << yMax << RESTendl;
457  x = xMin;
458 
459  for (Int_t i = 0; i < nBinsX; i++) {
460  z = zMin;
461  for (Int_t j = 0; j < nBinsZ; j++) {
462  Bvec = GetMagneticField(TVector3(x, y, z), false);
463  if (Bcomp == "X")
464  B = Bvec[0];
465  else {
466  if (Bcomp == "Y")
467  B = Bvec[1];
468  else {
469  if (Bcomp == "Z")
470  B = Bvec[2];
471  else
472  RESTError << "You entered : " << Bcomp
473  << " as a B component but you have to choose X, Y or Z" << RESTendl;
474  }
475  }
476  fHisto->Fill(x, z, B);
477  z = z + step_z;
478  }
479  x = x + step_x;
480  }
481 
482  fCanvas->cd();
483  fHisto->SetBit(TH1::kNoStats);
484  fHisto->GetXaxis()->SetTitle("x (mm)");
485  fHisto->GetYaxis()->SetTitle("z (mm)");
486 
487  if (Bcomp == "X") {
488  TString title = Form("B_{x} against x and z @ y = %.1f", y);
489  fHisto->SetTitle(title);
490  fCanvas->SetTitle(title);
491  }
492  if (Bcomp == "Y") {
493  TString title = Form("B_{y} against x and z @ y = %.1f", y);
494  fHisto->SetTitle(title);
495  fCanvas->SetTitle(title);
496  }
497  if (Bcomp == "Z") {
498  TString title = Form("B_{z} against x and z @ y = %.1f", y);
499  fHisto->SetTitle(title);
500  fCanvas->SetTitle(title);
501  }
502 
503  fHisto->Draw(style);
504  return fCanvas;
505  }
506 
507  if (projection == "YZ") {
508  TCanvas* fCanvas = new TCanvas("fCanvas", "");
509  fHisto = new TH2D("", "", nBinsY, yMin, yMax, nBinsZ, zMin, zMax);
510 
511  if (depth < -100000.0)
512  x = (xMin + xMax) / 2.0;
513  else if ((depth >= xMin) && (depth <= xMax))
514  x = depth;
515  else
516  RESTError << "You entered depth = " << depth << ", but you have to choose depth between " << xMin
517  << " and " << xMax << RESTendl;
518  y = yMin;
519 
520  for (Int_t i = 0; i < nBinsY; i++) {
521  z = zMin;
522  for (Int_t j = 0; j < nBinsZ; j++) {
523  Bvec = GetMagneticField(TVector3(x, y, z), false);
524  if (Bcomp == "X")
525  B = Bvec[0];
526  else {
527  if (Bcomp == "Y")
528  B = Bvec[1];
529  else {
530  if (Bcomp == "Z")
531  B = Bvec[2];
532  else
533  RESTError << "You entered : " << Bcomp
534  << " as a B component but you have to choose X, Y or Z" << RESTendl;
535  }
536  }
537  fHisto->Fill(y, z, B);
538  z = z + step_z;
539  }
540  y = y + step_y;
541  }
542 
543  fCanvas->cd();
544  fHisto->SetBit(TH1::kNoStats);
545  fHisto->GetXaxis()->SetTitle("y (mm)");
546  fHisto->GetYaxis()->SetTitle("z (mm)");
547 
548  if (Bcomp == "X") {
549  TString title = Form("B_{x} against y and z @ x = %.1f", x);
550  fHisto->SetTitle(title);
551  fCanvas->SetTitle(title);
552  }
553  if (Bcomp == "Y") {
554  TString title = Form("B_{y} against y and z @ x = %.1f", x);
555  fHisto->SetTitle(title);
556  fCanvas->SetTitle(title);
557  }
558  if (Bcomp == "Z") {
559  TString title = Form("B_{z} against y and z @ x = %.1f", x);
560  fHisto->SetTitle(title);
561  fCanvas->SetTitle(title);
562  }
563 
564  fHisto->Draw(style);
565  return fCanvas;
566  }
567 
568  return fCanvas;
569 }
570 
575 TCanvas* TRestAxionMagneticField::DrawTracks(TVector3 vanishingPoint, Int_t divisions, Int_t volId) {
576  if (fCanvas != NULL) {
577  delete fCanvas;
578  fCanvas = NULL;
579  }
580  fCanvas = new TCanvas("fCanvas", "", 1600, 600);
581 
582  TPad* pad1 = new TPad("pad1", "This is pad1", 0.01, 0.02, 0.99, 0.97);
583  pad1->Divide(2, 1);
584  pad1->Draw();
585 
586  pad1->cd(1);
587 
588  Double_t genSizeY = fBoundMax[volId].Y() * 3.;
589  Double_t genPositionZ = fPositions[volId][2] - fBoundMax[volId].Z() - 2000;
590  Double_t finalPositionZ = fPositions[volId][2] + fBoundMax[volId].Z() + 2000;
591 
592  // We generate the particle at different Y-highs
593  TGraph* bBox = new TGraph();
594  bBox->SetPoint(0, fPositions[volId][2] - fBoundMax[volId].Z(),
595  fPositions[volId][1] - fBoundMax[volId].Y());
596  bBox->SetPoint(1, fPositions[volId][2] - fBoundMax[volId].Z(),
597  fPositions[volId][1] + fBoundMax[volId].Y());
598  bBox->SetPoint(2, fPositions[volId][2] + fBoundMax[volId].Z(),
599  fPositions[volId][1] + fBoundMax[volId].Y());
600  bBox->SetPoint(3, fPositions[volId][2] + fBoundMax[volId].Z(),
601  fPositions[volId][1] - fBoundMax[volId].Y());
602  bBox->SetPoint(4, fPositions[volId][2] - fBoundMax[volId].Z(),
603  fPositions[volId][1] - fBoundMax[volId].Y());
604 
605  RESTDebug << "Gen position : " << genPositionZ << RESTendl;
606 
607  bBox->GetXaxis()->SetLimits(genPositionZ - 500, finalPositionZ + 500);
608  bBox->GetHistogram()->SetMaximum(genSizeY + 100);
609  bBox->GetHistogram()->SetMinimum(-genSizeY - 100);
610 
611  bBox->GetXaxis()->SetTitle("Z [mm]");
612  bBox->GetXaxis()->SetTitleSize(0.05);
613  bBox->GetXaxis()->SetLabelSize(0.05);
614  bBox->GetYaxis()->SetTitle("Y [mm]");
615  bBox->GetYaxis()->SetTitleOffset(1.3);
616  bBox->GetYaxis()->SetTitleSize(0.05);
617  bBox->GetYaxis()->SetLabelSize(0.05);
618  bBox->SetLineWidth(2);
619  bBox->Draw("AL");
620 
621  Int_t n = 0;
622  for (Double_t y = genSizeY; y >= -genSizeY; y -= fBoundMax[volId].Y()) {
623  TVector3 position(0, y, genPositionZ);
624  TVector3 direction = (vanishingPoint - position).Unit();
625 
626  std::vector<TVector3> trackBounds = this->GetFieldBoundaries(position, direction);
627  TGraph* grB = new TGraph();
628  grB->SetPoint(0, trackBounds[0].Z(), trackBounds[0].Y());
629  grB->SetPoint(1, trackBounds[1].Z(), trackBounds[1].Y());
630 
631  RESTDebug << "Initial" << RESTendl;
632  RESTDebug << "-------" << RESTendl;
634 
635  RESTDebug << RESTendl;
636  RESTDebug << "Start moving along" << RESTendl;
637  RESTDebug << "++++++++++++++++++" << RESTendl;
638 
639  TGraph* fieldGr = new TGraph();
640  Double_t posZ = fPositions[volId][2] - fBoundMax[volId].Z() - 10;
641  Double_t delta = fBoundMax[volId][2] * 2. / divisions;
642 
643  Double_t field = this->GetTransversalComponent(position, direction);
644  fieldGr->SetPoint(fieldGr->GetN(), posZ, field);
645 
646  while (posZ <= fPositions[volId][2] + fBoundMax[volId].Z()) {
647  TVector3 posAlongAxis = TVector3(fPositions[volId][0], fPositions[volId][1], posZ);
648 
649  position = MoveToPlane(position, direction, TVector3(0, 0, 1), posAlongAxis);
650  Double_t field = this->GetTransversalComponent(position, direction);
651 
652  fieldGr->SetPoint(fieldGr->GetN(), posZ, field);
653 
654  posZ += delta;
655  }
656 
657  TVector3 posAlongAxis = TVector3(fPositions[volId][0], fPositions[volId][1], posZ + 10);
658  position = MoveToPlane(position, direction, TVector3(0, 0, 1), posAlongAxis);
659 
660  Double_t field2 = this->GetTransversalComponent(position, direction);
661  fieldGr->SetPoint(fieldGr->GetN(), posZ, field2);
662 
663  pad1->cd(2);
664  fieldGr->SetLineWidth(3);
665  fieldGr->SetLineColor(38 + n);
666  fieldGr->GetXaxis()->SetLimits(genPositionZ - 500, finalPositionZ + 500);
667  fieldGr->GetHistogram()->SetMaximum(2.5);
668  fieldGr->GetHistogram()->SetMinimum(0);
669  fieldGr->GetXaxis()->SetTitle("Z [mm]");
670  fieldGr->GetXaxis()->SetTitleSize(0.05);
671  fieldGr->GetXaxis()->SetLabelSize(0.05);
672  fieldGr->GetYaxis()->SetTitle("B [T]");
673  fieldGr->GetYaxis()->SetTitleOffset(1.3);
674  fieldGr->GetYaxis()->SetTitleSize(0.05);
675  fieldGr->GetYaxis()->SetLabelSize(0.05);
676  if (y == genSizeY)
677  fieldGr->Draw("AL");
678  else
679  fieldGr->Draw("L");
680  pad1->cd(1);
681 
682  TGraph* gr = new TGraph();
683  gr->SetPoint(0, genPositionZ, y);
684  position = MoveToPlane(position, direction, TVector3(0, 0, 1), TVector3(0, 0, finalPositionZ));
685  gr->SetPoint(1, finalPositionZ, position.Y());
686 
687  gr->SetLineWidth(1.5);
688  gr->Draw("L");
689  grB->SetLineColor(38 + n);
690  n++;
691  grB->SetLineWidth(3);
692  grB->Draw("L");
693  }
694  bBox->Draw("L");
695 
696  return fCanvas;
697 }
698 
705 void TRestAxionMagneticField::LoadMagneticFieldData(MagneticFieldVolume& mVol,
706  std::vector<std::vector<Float_t>> data) {
707  mVol.field.resize(mVol.mesh.GetNodesX());
708  for (int n = 0; n < mVol.field.size(); n++) {
709  mVol.field[n].resize(mVol.mesh.GetNodesY());
710  for (int m = 0; m < mVol.field[n].size(); m++) mVol.field[n][m].resize(mVol.mesh.GetNodesZ());
711  }
712 
713  RESTDebug << "TRestAxionMagneticField::LoadMagneticFieldData. Printing first 5 data rows" << RESTendl;
714  for (Int_t n = 0; n < data.size(); n++) {
715  // The magnetic field map is centered at zero.
716  // But the mesh definition contains the offset position
717  // We shift the data to match the mesh node network.
718  Int_t nX = mVol.mesh.GetNodeX((Int_t)(data[n][0] + mVol.mesh.GetNetSizeX() / 2.), true);
719  Int_t nY = mVol.mesh.GetNodeY((Int_t)(data[n][1] + mVol.mesh.GetNetSizeY() / 2.), true);
720  Int_t nZ = mVol.mesh.GetNodeZ((Int_t)(data[n][2] + mVol.mesh.GetNetSizeZ() / 2.), true);
721 
722  if (n < 5) {
723  RESTDebug << "X: " << data[n][0] << " Y: " << data[n][1] << " Z: " << data[n][2] << RESTendl;
724  RESTDebug << "absX: " << data[n][0] + mVol.position.X()
725  << " absY: " << data[n][1] + mVol.position.Y()
726  << " absZ: " << data[n][2] + mVol.position.Z() << RESTendl;
727  RESTDebug << "nX: " << nX << " nY: " << nY << " nZ: " << nZ << RESTendl;
728  RESTDebug << "Bx: " << data[n][3] << " By: " << data[n][4] << " Bz: " << data[n][5] << RESTendl;
730  }
731 
732  if (mVol.field[nX][nY][nZ] != TVector3(0.0, 0.0, 0.0)) {
733  RESTWarning << "X: " << data[n][0] << " Y: " << data[n][1] << " Z: " << data[n][2] << RESTendl;
734  RESTWarning << "nX: " << nX << " nY: " << nY << " nZ: " << nZ << RESTendl;
735  RESTWarning << "WARNING: field[nX][nY][nZ] element not equal to initial value (0, 0, 0) !!"
736  << RESTendl;
737  RESTWarning << "It has value: "
738  << "mVol.field[" << nX << "][" << nY << "][" << nZ << "] = ("
739  << mVol.field[nX][nY][nZ].X() << " , " << mVol.field[nX][nY][nZ].Y() << " , "
740  << mVol.field[nX][nY][nZ].Z() << ")" << RESTendl;
741  RESTWarning << "Values to write: "
742  << "Bx: " << data[n][3] << " By: " << data[n][4] << " Bz: " << data[n][5] << RESTendl
743  << RESTendl;
744 
745  this->SetError("There was a problem assigning the field matrix!");
747  }
748 
749  mVol.field[nX][nY][nZ] = TVector3(data[n][3], data[n][4], data[n][5]);
750  }
751 }
752 
759  for (unsigned int n = 0; n < fPositions.size(); n++) {
760  string fullPathName = SearchFile((string)fFileNames[n]);
761  RESTDebug << "Reading file : " << fFileNames[n] << RESTendl;
762  RESTDebug << "Full path : " << fullPathName << RESTendl;
763 
764  if (fFileNames[n] != "none" && fullPathName == "") {
765  RESTError << "TRestAxionMagneticField::LoadMagneticVolumes. File " << fFileNames[n]
766  << " not found!" << RESTendl;
767  RESTError
768  << "REST will look for this file at any path given by <searchPath at globals definitions"
769  << RESTendl;
770  exit(5);
771  }
772 
773  std::vector<std::vector<Float_t>> fieldData;
774  if (fFileNames[n] != "none")
775  if (fullPathName.find(".dat") != string::npos) {
776  RESTDebug << "Reading ASCII format" << RESTendl;
777  if (!TRestTools::ReadASCIITable(fullPathName, fieldData)) {
778  RESTError << "Problem reading file : " << fullPathName << RESTendl;
779  exit(1);
780  }
781  } else {
782  if (fullPathName.find(".bin") != string::npos) {
783  RESTDebug << "Reading binary format" << RESTendl;
784  if (!TRestTools::ReadBinaryTable(fullPathName, fieldData, 6)) {
785  RESTError << "Problem reading file : " << fullPathName << RESTendl;
786  exit(2);
787  }
788  }
789  }
790  else if (fFileNames[n] != "none") {
791  RESTError << "Filename : " << fullPathName << RESTendl;
792  RESTError << "File format not recognized!" << RESTendl;
793  exit(3);
794  }
795 
796  if (fFileNames[n] != "none" && fieldData.size() < 2) {
797  RESTError << "Field data size is no more than 2 grid points!" << RESTendl;
798  RESTError << "Filename : " << fullPathName << RESTendl;
799  RESTError << "Probably something went wrong loading the file" << RESTendl;
800  exit(4);
801  }
802 
803  Float_t xMin = -fBoundMax[n].X(), yMin = -fBoundMax[n].Y(), zMin = -fBoundMax[n].Z();
804  Float_t xMax = fBoundMax[n].X(), yMax = fBoundMax[n].Y(), zMax = fBoundMax[n].Z();
805  Float_t meshSizeX = fMeshSize[n].X(), meshSizeY = fMeshSize[n].Y(), meshSizeZ = fMeshSize[n].Z();
806 
807  // If a field map is defined we get the boundaries, and mesh size from the volume
808  if (fieldData.size() > 0) {
809  RESTDebug << "Reading max boundary values" << RESTendl;
810  xMax = TRestTools::GetMaxValueFromTable(fieldData, 0);
811  yMax = TRestTools::GetMaxValueFromTable(fieldData, 1);
812  zMax = TRestTools::GetMaxValueFromTable(fieldData, 2);
813 
814  if (fBoundMax[n] != TVector3(0, 0, 0)) {
815  if (fBoundMax[n] != TVector3(xMax, yMax, zMax)) {
816  RESTWarning << "Volume : " << n << RESTendl;
817  RESTWarning << "boundMax was defined in RML but does not match the field map boundaries!"
818  << RESTendl;
819  RESTWarning << "Max. Field map boundaries : (" << xMax << ", " << yMax << ", " << zMax
820  << ")" << RESTendl;
821  }
822  }
823 
824  RESTDebug << "Reading min boundary values" << RESTendl;
825  xMin = TRestTools::GetMinValueFromTable(fieldData, 0);
826  yMin = TRestTools::GetMinValueFromTable(fieldData, 1);
827  zMin = TRestTools::GetMinValueFromTable(fieldData, 2);
828 
829  if (fBoundMax[n] != TVector3(0, 0, 0)) {
830  if (-fBoundMax[n] != TVector3(xMin, yMin, zMin)) {
831  RESTWarning << "Volume : " << n << RESTendl;
832  RESTWarning << "boundMax was defined in RML but does not match the field map boundaries"
833  << RESTendl;
834  RESTWarning << "Min. Field map boundaries : (" << xMin << ", " << yMin << ", " << zMin
835  << ")" << RESTendl;
836  }
837  }
838  fBoundMax[n] = TVector3(xMax, yMax, zMax);
839 
840  RESTDebug << "Reading mesh size" << RESTendl;
841  meshSizeX = TRestTools::GetLowestIncreaseFromTable(fieldData, 0);
842  meshSizeY = TRestTools::GetLowestIncreaseFromTable(fieldData, 1);
843  meshSizeZ = TRestTools::GetLowestIncreaseFromTable(fieldData, 2);
844 
845  if (fMeshSize[n] != TVector3(0, 0, 0)) {
846  if (fMeshSize[n] != TVector3(meshSizeX, meshSizeY, meshSizeZ)) {
847  RESTWarning << "Volume : " << n << RESTendl;
848  RESTWarning
849  << "MeshSize was defined in RML but does not match the mesh size deduced from "
850  "field map"
851  << RESTendl;
852  RESTWarning << "Mesh size : (" << meshSizeX << ", " << meshSizeY << ", " << meshSizeZ
853  << ")" << RESTendl;
854  }
855  }
856  fMeshSize[n] = TVector3(meshSizeX, meshSizeY, meshSizeZ);
857  }
858 
860  RESTDebug << "Reading magnetic field map" << RESTendl;
861  RESTDebug << "--------------------------" << RESTendl;
862 
863  RESTDebug << "Full path : " << fullPathName << RESTendl;
864 
865  RESTDebug << "Boundaries" << RESTendl;
866  RESTDebug << "xMin: " << xMin << " yMin: " << yMin << " zMin: " << zMin << RESTendl;
867  RESTDebug << "xMax: " << xMax << " yMax: " << yMax << " zMax: " << zMax << RESTendl;
868  RESTDebug << "Mesh size" << RESTendl;
869 
870  RESTDebug << "sX: " << meshSizeX << " sY: " << meshSizeY << " sZ: " << meshSizeZ << RESTendl;
871 
872  if (fieldData.size() > 4) {
873  RESTDebug << "Printing beginning of magnetic file table : " << fieldData.size() << RESTendl;
874  TRestTools::PrintTable(fieldData, 0, 5);
875  } else {
876  RESTDebug << "The data file contains no field map" << RESTendl;
877  }
878  }
880 
881  // Number of nodes
882  Int_t nx = (Int_t)(2 * xMax / meshSizeX) + 1;
883  Int_t ny = (Int_t)(2 * yMax / meshSizeY) + 1;
884  Int_t nz = (Int_t)(2 * zMax / meshSizeZ) + 1;
885 
886  // We create an auxiliar mesh helping to initialize the fieldMap
887  // The mesh is centered at zero. Absolute position is defined in the Magnetic volume
888  // TODO It would be interesting that TRestMesh could be used in cylindrical coordinates.
889  TRestMesh restMesh;
890  restMesh.SetSize(2 * xMax, 2 * yMax, 2 * zMax);
891  restMesh.SetOrigin(fPositions[n] - TVector3(xMax, yMax, zMax));
892  restMesh.SetNodes(nx, ny, nz);
893  if (fMeshType[n] == "cylinder")
894  restMesh.SetCylindrical(true);
895  else
896  restMesh.SetCylindrical(false);
897 
898  MagneticFieldVolume mVolume;
899  mVolume.position = fPositions[n];
900  mVolume.mesh = restMesh;
901 
902  if (fieldData.size() > 0) LoadMagneticFieldData(mVolume, fieldData);
903 
904  if (fBoundMax[n] == TVector3(0, 0, 0)) {
905  RESTError << "The bounding box was not defined for volume " << n << "!" << RESTendl;
906  RESTError << "Please review RML configuration for TRestAxionMagneticField" << RESTendl;
907  exit(22);
908  } else if (fMeshSize[n] == TVector3(0, 0, 0)) {
909  RESTError << "The mesh grid size was not defined for volume " << n << "!" << RESTendl;
910  RESTError << "Please review RML configuration for TRestAxionMagneticField" << RESTendl;
911  exit(22);
912  }
913  fMagneticFieldVolumes.push_back(mVolume);
914  }
915 
916  if (CheckOverlaps()) {
917  RESTError << "TRestAxionMagneticField::LoadMagneticVolumes. Volumes overlap!" << RESTendl;
918  exit(1);
919  }
920  RESTDebug << "Finished loading magnetic volumes" << RESTendl;
921 }
922 
926 TVector3 TRestAxionMagneticField::GetMagneticField(Double_t x, Double_t y, Double_t z) {
927  return GetMagneticField(TVector3(x, y, z));
928 }
929 
938 TVector3 TRestAxionMagneticField::GetMagneticField(TVector3 pos, Bool_t showWarning) {
939  Int_t id = GetVolumeIndex(pos);
940 
941  if (id < 0) {
942  if (showWarning)
943  RESTWarning << "TRestAxionMagneticField::GetMagneticField position (" << pos.X() << ", "
944  << pos.Y() << ", " << pos.Z() << ") is outside any volume" << RESTendl;
945  return TVector3(0, 0, 0);
946  } else {
947  if (IsFieldConstant(id)) return fConstantField[id];
948  TVector3 node = GetMagneticVolumeNode(fMagneticFieldVolumes[id], pos);
949  Int_t nX = node.X();
950  Int_t nY = node.Y();
951  Int_t nZ = node.Z();
952 
953  Int_t nX_1, nY_1, nZ_1;
954 
955  if ((nX + 1) < fMagneticFieldVolumes[id].mesh.GetNodesX())
956  nX_1 = nX + 1;
957  else
958  nX_1 = nX;
959  if ((nY + 1) < fMagneticFieldVolumes[id].mesh.GetNodesY())
960  nY_1 = nY + 1;
961  else
962  nY_1 = nY;
963  if ((nZ + 1) < fMagneticFieldVolumes[id].mesh.GetNodesZ())
964  nZ_1 = nZ + 1;
965  else
966  nZ_1 = nZ;
967 
968  TVector3 C000 = fMagneticFieldVolumes[id].field[nX][nY][nZ] + fConstantField[id];
969  TVector3 C100 = fMagneticFieldVolumes[id].field[nX_1][nY][nZ] + fConstantField[id];
970  TVector3 C010 = fMagneticFieldVolumes[id].field[nX][nY_1][nZ] + fConstantField[id];
971  TVector3 C110 = fMagneticFieldVolumes[id].field[nX_1][nY_1][nZ] + fConstantField[id];
972  TVector3 C001 = fMagneticFieldVolumes[id].field[nX][nY][nZ_1] + fConstantField[id];
973  TVector3 C101 = fMagneticFieldVolumes[id].field[nX_1][nY][nZ_1] + fConstantField[id];
974  TVector3 C011 = fMagneticFieldVolumes[id].field[nX][nY_1][nZ_1] + fConstantField[id];
975  TVector3 C111 = fMagneticFieldVolumes[id].field[nX_1][nY_1][nZ_1] + fConstantField[id];
976 
977  Double_t x0 = fMagneticFieldVolumes[id].mesh.GetX(nX);
978  Double_t x1 = fMagneticFieldVolumes[id].mesh.GetX(nX_1);
979  Double_t xd;
980  if (x0 == x1)
981  xd = 0;
982  else
983  xd = (pos.X() - x0) / (x1 - x0);
984  if ((xd < -0.00001) || (xd > 1.00001))
985  RESTWarning << "TRestAxionMagneticField::GetMagneticField Error: xd NOT between 0 and 1"
986  << RESTendl;
987 
988  Double_t y0 = fMagneticFieldVolumes[id].mesh.GetY(nY);
989  Double_t y1 = fMagneticFieldVolumes[id].mesh.GetY(nY_1);
990  Double_t yd;
991  if (y0 == y1)
992  yd = 0;
993  else
994  yd = (pos.Y() - y0) / (y1 - y0);
995  if ((yd < -0.00001) || (yd > 1.00001))
996  RESTWarning << "TRestAxionMagneticField::GetMagneticField Error: yd NOT between 0 and 1"
997  << RESTendl;
998 
999  Double_t z0 = fMagneticFieldVolumes[id].mesh.GetZ(nZ);
1000  Double_t z1 = fMagneticFieldVolumes[id].mesh.GetZ(nZ_1);
1001  Double_t zd;
1002  if (z0 == z1)
1003  zd = 0;
1004  else
1005  zd = (pos.Z() - z0) / (z1 - z0);
1006  if ((zd < -0.00001) || (zd > 1.00001))
1007  RESTWarning << "TRestAxionMagneticField::GetMagneticField Error: zd NOT between 0 and 1"
1008  << RESTendl;
1009 
1010  // first we interpolate along x-axis
1011  TVector3 C00 = C000 * (1.0 - xd) + C100 * xd;
1012  TVector3 C01 = C001 * (1.0 - xd) + C101 * xd;
1013  TVector3 C10 = C010 * (1.0 - xd) + C110 * xd;
1014  TVector3 C11 = C011 * (1.0 - xd) + C111 * xd;
1015 
1016  // then we interpolate along y-axis
1017  TVector3 C0 = C00 * (1.0 - yd) + C10 * yd;
1018  TVector3 C1 = C01 * (1.0 - yd) + C11 * yd;
1019 
1020  // finally we interpolate along z-axis
1021  TVector3 C = C0 * (1.0 - zd) + C1 * zd;
1022 
1023  RESTDebug << "position = (" << pos.X() << ", " << pos.Y() << ", " << pos.Z() << ") ";
1024  RESTDebug << "nX = " << nX << " nY = " << nY << " nZ = " << nZ << " nX_1 = " << nX_1
1025  << " nY_1 = " << nY_1 << " nZ_1 = " << nZ_1 << RESTendl << RESTendl;
1026  RESTDebug << "C000 = (" << C000.X() << ", " << C000.Y() << ", " << C000.Z() << ")" << RESTendl
1027  << RESTendl;
1028  RESTDebug << "C100 = (" << C100.X() << ", " << C100.Y() << ", " << C100.Z() << ")" << RESTendl
1029  << RESTendl;
1030  RESTDebug << "C010 = (" << C010.X() << ", " << C010.Y() << ", " << C010.Z() << ")" << RESTendl
1031  << RESTendl;
1032  RESTDebug << "C110 = (" << C110.X() << ", " << C110.Y() << ", " << C110.Z() << ")" << RESTendl
1033  << RESTendl;
1034  RESTDebug << "C001 = (" << C001.X() << ", " << C001.Y() << ", " << C001.Z() << ")" << RESTendl
1035  << RESTendl;
1036  RESTDebug << "C101 = (" << C101.X() << ", " << C101.Y() << ", " << C101.Z() << ")" << RESTendl
1037  << RESTendl;
1038  RESTDebug << "C011 = (" << C011.X() << ", " << C011.Y() << ", " << C011.Z() << ")" << RESTendl
1039  << RESTendl;
1040  RESTDebug << "C111 = (" << C111.X() << ", " << C111.Y() << ", " << C111.Z() << ")" << RESTendl
1041  << RESTendl;
1042  RESTDebug << " -------------------------------------------------------" << RESTendl;
1043  RESTDebug << "C00 = (" << C00.X() << ", " << C00.Y() << ", " << C00.Z() << ")" << RESTendl
1044  << RESTendl;
1045  RESTDebug << "C01 = (" << C01.X() << ", " << C01.Y() << ", " << C01.Z() << ")" << RESTendl
1046  << RESTendl;
1047  RESTDebug << "C10 = (" << C10.X() << ", " << C10.Y() << ", " << C10.Z() << ")" << RESTendl
1048  << RESTendl;
1049  RESTDebug << "C11 = (" << C11.X() << ", " << C11.Y() << ", " << C11.Z() << ")" << RESTendl
1050  << RESTendl;
1051  RESTDebug << " -------------------------------------------------------" << RESTendl;
1052  RESTDebug << "C0 = (" << C0.X() << ", " << C0.Y() << ", " << C0.Z() << ")" << RESTendl << RESTendl;
1053  RESTDebug << "C1 = (" << C1.X() << ", " << C1.Y() << ", " << C1.Z() << ")" << RESTendl << RESTendl;
1054  RESTDebug << " -------------------------------------------------------" << RESTendl;
1055  RESTDebug << "C = (" << C.X() << ", " << C.Y() << ", " << C.Z() << ")" << RESTendl << RESTendl;
1056 
1057  return C;
1058  }
1059 }
1060 
1066  if (!FieldLoaded()) LoadMagneticVolumes();
1067 
1068  for (int n = 0; n < fMagneticFieldVolumes.size(); n++) {
1069  if (fMagneticFieldVolumes[n].mesh.IsInside(pos)) return n;
1070  }
1071  return -1;
1072 }
1073 
1078  if (GetVolumeIndex(pos) >= 0) return true;
1079  return false;
1080 }
1081 
1086 
1091  if (GetNumberOfVolumes() > id)
1092  return fPositions[id];
1093  else {
1094  RESTWarning << "TRestAxionMagneticField::GetVolumePosition. Id : " << id << " out of range!"
1095  << RESTendl;
1096  RESTWarning << "Number of volumes defined : " << GetNumberOfVolumes() << RESTendl;
1097  return TVector3(0, 0, 0);
1098  }
1099 }
1100 
1105 Double_t TRestAxionMagneticField::GetTransversalComponent(TVector3 position, TVector3 direction) {
1106  return abs(GetMagneticField(position, false).Perp(direction));
1107 }
1108 
1120 std::vector<Double_t> TRestAxionMagneticField::GetTransversalComponentAlongPath(TVector3 from, TVector3 to,
1121  Double_t dl, Int_t Nmax) {
1122  Double_t length = (to - from).Mag();
1123 
1124  Double_t diff = dl;
1125  if (Nmax > 0) {
1126  if (length / dl > Nmax) {
1127  diff = length / Nmax;
1128  RESTWarning << "TRestAxionMagneticField::GetTransversalComponentAlongPath. Nmax reached!"
1129  << RESTendl;
1130  RESTWarning << "Nmax = " << Nmax << RESTendl;
1131  RESTWarning << "Adjusting differential step to : " << diff << " mm" << RESTendl;
1132  }
1133  }
1134 
1135  TVector3 direction = (to - from).Unit();
1136 
1137  std::vector<Double_t> Bt;
1138  for (Double_t d = 0; d < length; d += diff) {
1139  Bt.push_back(GetTransversalComponent(from + d * direction, direction));
1140  }
1141 
1142  return Bt;
1143 }
1144 
1155 Double_t TRestAxionMagneticField::GetTransversalFieldAverage(TVector3 from, TVector3 to, Double_t dl,
1156  Int_t Nmax) {
1157  Double_t length = (to - from).Mag();
1158 
1159  Double_t Bavg = 0.;
1160  std::vector<Double_t> Bt = GetTransversalComponentAlongPath(from, to, dl, Nmax);
1161  for (auto& b : Bt) Bavg += b;
1162 
1163  if (length > 0) return Bavg / length;
1164 
1165  RESTError << "TRestAxionMagneticField::GetTransversalFieldAverage. Lenght is zero!" << RESTendl;
1166  return 0.;
1167 }
1168 
1180 TVector3 TRestAxionMagneticField::GetFieldAverageTransverseVector(TVector3 from, TVector3 to, Double_t dl,
1181  Int_t Nmax) {
1182  Double_t length = (to - from).Mag();
1183 
1184  Double_t diff = dl;
1185  if (Nmax > 0) {
1186  if (length / dl > Nmax) {
1187  diff = length / Nmax;
1188  RESTWarning << "TRestAxionMagneticField::GetFieldAverageTransverseVector Nmax reached!"
1189  << RESTendl;
1190  RESTWarning << "Nmax = " << Nmax << RESTendl;
1191  RESTWarning << "Adjusting differential step to : " << diff << " mm" << RESTendl;
1192  }
1193  }
1194 
1195  TVector3 direction = (to - from).Unit();
1196  TVector3 Bavg = TVector3(0.0, 0.0, 0.0);
1197  TVector3 BTavg = TVector3(0.0, 0.0, 0.0);
1198  Int_t numberofpoints = 0;
1199 
1200  for (Double_t d = 0; d <= length; d += diff) {
1201  Bavg = Bavg + GetMagneticField(from + d * direction);
1202  numberofpoints = numberofpoints + 1;
1203  }
1204 
1205  if ((length > 0) && (numberofpoints > 0)) {
1206  Bavg = Bavg * (1.0 / numberofpoints); // calculates the average magnetic field vector
1207  BTavg =
1208  Bavg - (Bavg * direction) *
1209  direction; // calculates the transverse component of the average magnetic field vector
1210  RESTDebug << "B average vector = (" << Bavg.x() << ", " << Bavg.y() << ", " << Bavg.z() << ")"
1211  << RESTendl;
1212  RESTDebug << "Transverse B average vector = (" << BTavg.x() << ", " << BTavg.y() << ", " << BTavg.z()
1213  << ")" << RESTendl;
1214  return BTavg;
1215  }
1216  RESTError << "TRestAxionMagneticField::GetTransversalFieldAverage. Lenght is zero!" << RESTendl;
1217  return TVector3(0.0, 0.0, 0.0);
1218 }
1219 
1228 TVector3 TRestAxionMagneticField::GetMagneticVolumeNode(MagneticFieldVolume mVol, TVector3 pos) {
1229  Int_t nx = mVol.mesh.GetNodeX(pos.X());
1230  Int_t ny = mVol.mesh.GetNodeY(pos.Y());
1231  Int_t nz = mVol.mesh.GetNodeZ(pos.Z());
1232  return TVector3(nx, ny, nz);
1233 }
1234 
1239  RESTDebug << "Checking overlaps" << RESTendl;
1240  for (int n = 0; n < GetNumberOfVolumes(); n++) {
1241  for (int m = 0; m < GetNumberOfVolumes(); m++) {
1242  if (m == n) continue;
1243  RESTDebug << "Volume : " << m << RESTendl;
1244 
1245  TVector3 b = GetMagneticVolume(m)->mesh.GetVertex(0);
1246  RESTDebug << "Relative bottom vertex : (" << b.X() << ", " << b.Y() << ", " << b.Z() << ")"
1247  << RESTendl;
1248  if (GetMagneticVolume(n)->mesh.IsInsideBoundingBox(b)) return true;
1249 
1250  TVector3 t = GetMagneticVolume(m)->mesh.GetVertex(1);
1251  RESTDebug << "Relative top vertex : (" << t.X() << ", " << t.Y() << ", " << t.Z() << ")"
1252  << RESTendl;
1253  if (GetMagneticVolume(n)->mesh.IsInsideBoundingBox(t)) return true;
1254  }
1255  }
1256  return false;
1257 }
1258 
1272 std::vector<TVector3> TRestAxionMagneticField::GetVolumeBoundaries(TVector3 pos, TVector3 dir, Int_t id) {
1273  MagneticFieldVolume* vol = GetMagneticVolume(id);
1274 
1275  std::vector<TVector3> boundaries;
1276 
1277  if (vol) boundaries = vol->mesh.GetTrackBoundaries(pos, dir);
1278 
1279  return boundaries;
1280 }
1281 
1292 std::vector<TVector3> TRestAxionMagneticField::GetFieldBoundaries(TVector3 pos, TVector3 dir,
1293  Double_t precision, Int_t id) {
1294  std::vector<TVector3> volumeBoundaries = GetVolumeBoundaries(pos, dir, id);
1295  if (volumeBoundaries.size() != 2) return volumeBoundaries;
1296 
1297  if (IsFieldConstant(id)) return volumeBoundaries;
1298 
1299  MagneticFieldVolume* vol = GetMagneticVolume(id);
1300  if (!vol) return volumeBoundaries;
1301 
1302  if (precision == 0) precision = min(fMeshSize[id].X(), min(fMeshSize[id].Y(), fMeshSize[id].Z())) / 2.;
1303 
1304  TVector3 unit = dir.Unit();
1305  std::vector<TVector3> fieldBoundaries;
1306 
1307  TVector3 in = volumeBoundaries[0];
1308  while ((((volumeBoundaries[1] - in) * dir) > 0) && (GetTransversalComponent(in, dir) == 0))
1309  in = MoveByDistanceFast(in, unit, precision);
1310  if (((volumeBoundaries[1] - in) * dir) > 0)
1311  fieldBoundaries.push_back(in);
1312  else
1313  return fieldBoundaries;
1314 
1315  TVector3 out = volumeBoundaries[1];
1316  while ((((volumeBoundaries[0] - out) * dir) < 0) && (GetTransversalComponent(out, -dir) == 0) &&
1317  (((out - in) * dir) > 0))
1318  out = MoveByDistanceFast(out, -unit, precision);
1319  if ((((volumeBoundaries[0] - out) * dir) < 0) && (((out - in) * dir) > 0))
1320  fieldBoundaries.push_back(out);
1321  else
1322  return fieldBoundaries;
1323 
1324  return fieldBoundaries;
1325 }
1326 
1331  this->Initialize();
1332 
1333  auto magVolumeDef = GetElement("addMagneticVolume");
1334  while (magVolumeDef) {
1335  string filename = GetFieldValue("fileName", magVolumeDef);
1336  if (filename == "Not defined")
1337  fFileNames.push_back("none");
1338  else
1339  fFileNames.push_back(filename);
1340 
1341  TVector3 position = Get3DVectorParameterWithUnits("position", magVolumeDef);
1342  if (position == TVector3(-1, -1, -1))
1343  fPositions.push_back(TVector3(0, 0, 0));
1344  else
1345  fPositions.push_back(position);
1346 
1347  TVector3 field = Get3DVectorParameterWithUnits("field", magVolumeDef);
1348  if (field == TVector3(-1, -1, -1))
1349  fConstantField.push_back(TVector3(0, 0, 0));
1350  else
1351  fConstantField.push_back(field);
1352 
1353  TVector3 boundMax = Get3DVectorParameterWithUnits("boundMax", magVolumeDef);
1354  if (boundMax == TVector3(-1, -1, -1))
1355  fBoundMax.push_back(TVector3(0, 0, 0));
1356  else
1357  fBoundMax.push_back(boundMax);
1358 
1359  TVector3 meshSize = Get3DVectorParameterWithUnits("meshSize", magVolumeDef);
1360  if (meshSize == TVector3(-1, -1, -1))
1361  fMeshSize.push_back(TVector3(0, 0, 0));
1362  else
1363  fMeshSize.push_back(meshSize);
1364 
1365  string type = GetParameter("meshType", magVolumeDef);
1366  if (type == "NO_SUCH_PARA")
1367  fMeshType.push_back("cylinder");
1368  else
1369  fMeshType.push_back(type);
1370 
1371  // TRestMesh will only consider the first bounding component anyway
1372  if (fMeshType.back() == "cylinder" && fBoundMax.back().X() != fBoundMax.back().Y()) {
1373  RESTWarning << "Mesh type is cylinder. But X and Y inside boundMax are not the same!" << RESTendl;
1374  RESTWarning << "Making second bound component Y equal to the X bound component!" << RESTendl;
1375  fBoundMax.back().SetY(fBoundMax.back().X());
1376  }
1377 
1378  RESTDebug << "Reading new magnetic volume" << RESTendl;
1379  RESTDebug << "-----" << RESTendl;
1380  RESTDebug << "Filename : " << filename << RESTendl;
1381  RESTDebug << "Position: ( " << position.X() << ", " << position.Y() << ", " << position.Z() << ") mm"
1382  << RESTendl;
1383  RESTDebug << "Field: ( " << field.X() << ", " << field.Y() << ", " << field.Z() << ") T" << RESTendl;
1384  RESTDebug << "Max bounding box ( " << boundMax.X() << ", " << boundMax.Y() << ", " << boundMax.Z()
1385  << ")" << RESTendl;
1386  RESTDebug << "Mesh size ( " << meshSize.X() << ", " << meshSize.Y() << ", " << meshSize.Z() << ")"
1387  << RESTendl;
1388  RESTDebug << "----" << RESTendl;
1389 
1390  magVolumeDef = GetNextElement(magVolumeDef);
1391  }
1392 
1394 
1395  // TODO we should check that the volumes do not overlap
1396 }
1397 
1403 
1404  RESTMetadata << " - Number of magnetic volumes : " << GetNumberOfVolumes() << RESTendl;
1405  RESTMetadata << " ------------------------------------------------ " << RESTendl;
1406  for (int p = 0; p < GetNumberOfVolumes(); p++) {
1407  if (p > 0) RESTMetadata << " ------------------------------------------------ " << RESTendl;
1408 
1409  Double_t centerX = fPositions[p][0];
1410  Double_t centerY = fPositions[p][1];
1411  Double_t centerZ = fPositions[p][2];
1412 
1413  Double_t halfSizeX = fBoundMax[p].X();
1414  Double_t halfSizeY = fBoundMax[p].Y();
1415  Double_t halfSizeZ = fBoundMax[p].Z();
1416 
1417  Double_t xMin = centerX - halfSizeX;
1418  Double_t yMin = centerY - halfSizeY;
1419  Double_t zMin = centerZ - halfSizeZ;
1420 
1421  Double_t xMax = centerX + halfSizeX;
1422  Double_t yMax = centerY + halfSizeY;
1423  Double_t zMax = centerZ + halfSizeZ;
1424 
1425  RESTMetadata << "* Volume " << p << " centered at (" << centerX << "," << centerY << "," << centerZ
1426  << ") mm" << RESTendl;
1427  RESTMetadata << " - Grid mesh element size. X: " << fMeshSize[p].X() << "mm "
1428  << " Y: " << fMeshSize[p].Y() << "mm "
1429  << " Z: " << fMeshSize[p].Z() << "mm " << RESTendl;
1430  RESTMetadata << " - Offset field [T] : (" << fConstantField[p].X() << ", " << fConstantField[p].Y()
1431  << ", " << fConstantField[p].Z() << ")" << RESTendl;
1432  RESTMetadata << " - File loaded : " << fFileNames[p] << RESTendl;
1433  RESTMetadata << " " << RESTendl;
1434  RESTMetadata << " - Bounds : " << RESTendl;
1435  RESTMetadata << " xmin : " << xMin << " mm , xmax : " << xMax << " mm" << RESTendl;
1436  RESTMetadata << " ymin : " << yMin << " mm, ymax : " << yMax << " mm" << RESTendl;
1437  RESTMetadata << " zmin : " << zMin << " mm, zmax : " << zMax << " mm" << RESTendl;
1438  RESTMetadata << " " << RESTendl;
1439  }
1440  RESTMetadata << "+++++++++++++++++++++++++++++++++++++++++++++++++" << RESTendl;
1441 }
std::vector< TVector3 > fMeshSize
The size of a grid element from the mesh in mm.
TVector3 MoveToPlane(TVector3 pos, TVector3 dir, TVector3 n, TVector3 a)
This method will translate the vector with direction dir starting at position pos to the plane define...
std::vector< TVector3 > fBoundMax
A vector to store the maximum bounding box values.
Int_t GetVolumeIndex(TVector3 pos)
It returns the corresponding volume index at the given position. If not found it will return -1...
std::vector< TVector3 > GetFieldBoundaries(TVector3 pos, TVector3 dir, Double_t precision=0, Int_t id=0)
Finds the in/out particle trajectory boundaries for a particular magnetic region, similar to the meth...
void SetLibraryVersion(TString version)
Set the library version of this metadata class.
A class to load magnetic field maps and evaluate the field on those maps including interpolation...
void SetSectionName(std::string sName)
set the section name, clear the section content
std::vector< TString > fMeshType
The type of the mesh used (default is cylindrical)
virtual void PrintMetadata()
Implemented it in the derived metadata class to print out specific metadata information.
TVector3 GetVolumePosition(Int_t id)
it returns the volume position (or center) for the given volume id.
static T GetMinValueFromTable(const std::vector< std::vector< T >> &data, Int_t column=-1)
It returns the minimum value for a particular column from the table given by argument. If no column is specified in the arguments, then it gets the minimum from the full table.
Definition: TRestTools.cxx:394
TVector3 GetVolumeCenter(Int_t id)
it returns the volume position (or center) for the given volume id.
void SetCylindrical(Bool_t v)
Sets the coordinate system to cylindrical.
Definition: TRestMesh.h:136
void LoadMagneticFieldData(MagneticFieldVolume &mVol, std::vector< std::vector< Float_t >> data)
A method to help loading magnetic field data, as x,y,z,Bx,By,Bz into a magnetic volume definition usi...
void LoadMagneticVolumes()
It will load the magnetic field data from the data filenames specified at the RML definition...
TiXmlElement * GetElement(std::string eleDeclare, TiXmlElement *e=nullptr)
Get an xml element from a given parent element, according to its declaration.
static int ReadASCIITable(std::string fName, std::vector< std::vector< Double_t >> &data, Int_t skipLines=0)
Reads an ASCII file containing a table with values.
Definition: TRestTools.cxx:511
STL namespace.
TVector3 GetMagneticVolumeNode(MagneticFieldVolume mVol, TVector3 pos)
It returns the corresponding mesh node in the magnetic volume.
void SetSize(Double_t sX, Double_t sY, Double_t sZ)
Sets the origin of the bounding-box and initializes the nodes vector to zero.
Definition: TRestMesh.cxx:527
std::vector< TVector3 > fConstantField
A constant field component that will be added to the field map.
TCanvas * DrawTracks(TVector3 vanishingPoint, Int_t divisions, Int_t volId=0)
A method that creates a canvas where tracks traversing the magnetic volume are drawm together with th...
TVector3 GetMagneticField(Double_t x, Double_t y, Double_t z)
It returns the magnetic field vector at x,y,z.
static int ReadBinaryTable(std::string fName, std::vector< std::vector< T >> &data, Int_t columns=-1)
Reads a binary file containing a fixed-columns table with values.
Definition: TRestTools.cxx:242
static T GetMaxValueFromTable(const std::vector< std::vector< T >> &data, Int_t column=-1)
It returns the maximum value for a particular column from the table given by argument. If no column is specified in the arguments, then it gets the maximum from the full table.
Definition: TRestTools.cxx:359
std::vector< TVector3 > GetVolumeBoundaries(TVector3 pos, TVector3 dir, Int_t id=0)
Finds the in/out particle trajectory boundaries for a particular magnetic region bounding box...
std::string fConfigFileName
Full name of the rml file.
Int_t GetNumberOfVolumes()
The number of magnetic volumes loaded into the object.
+show most of the information for each steps
TCanvas * fCanvas
A canvas to insert the histogram drawing.
std::vector< TVector3 > fPositions
The absolute position of each of the magnetic volumes defined in this class.
void PrintMetadata()
Prints on screen the information about the metadata members of TRestAxionMagneticField.
+show the defined debug messages
static T GetLowestIncreaseFromTable(std::vector< std::vector< T >> data, Int_t column)
It returns the lowest increase, different from zero, between the elements of a particular column from...
Definition: TRestTools.cxx:431
std::vector< std::string > fFileNames
The name of the filenames containing the field data.
void SetNodes(Int_t nX, Int_t nY, Int_t nZ)
Sets the number of nodes and initializes the nodes vector to zero.
Definition: TRestMesh.cxx:539
TRestStringOutput::REST_Verbose_Level GetVerboseLevel()
returns the verboselevel in type of REST_Verbose_Level enumerator
Bool_t CheckOverlaps()
It will return true if the magnetic the regions overlap.
Double_t GetTransversalComponent(TVector3 position, TVector3 direction)
It returns the intensity of the transversal magnetic field component for the defined propagation dire...
TiXmlElement * GetNextElement(TiXmlElement *e)
Get the next sibling xml element of this element, with same eleDeclare.
TRestAxionMagneticField()
Default constructor.
Bool_t IsInside(TVector3 pos)
It returns true if the given position is found inside a magnetic volume. False otherwise.
Double_t GetTransversalFieldAverage(TVector3 from, TVector3 to, Double_t dl=1., Int_t Nmax=0)
It returns the average of the transversal magnetic field intensity between the 3-d coordinates from a...
Int_t GetChar(std::string hint="Press a KEY to continue ...")
Helps to pause the program, printing a message before pausing.
Bool_t IsFieldConstant(Int_t id)
It returns true if no magnetic field map was loaded for that volume.
static int PrintTable(std::vector< std::vector< T >> data, Int_t start=0, Int_t end=0)
Prints the contents of the vector table given as argument in screen. Allowed types are Int_t...
Definition: TRestTools.cxx:154
std::string GetFieldValue(std::string parName, TiXmlElement *e)
Returns the field value of an xml element which has the specified name.
void SetOrigin(Double_t oX, Double_t oY, Double_t oZ)
Sets the origin of the bounding-box and initializes the nodes vector to zero.
Definition: TRestMesh.cxx:507
std::string SearchFile(std::string filename)
Search files in current directory and directories specified in "searchPath" section.
void Initialize()
Initialization of TRestAxionMagneticField members.
Bool_t FieldLoaded()
This private method returns true if the magnetic field volumes loaded are the same as the volumes def...
Int_t LoadConfigFromFile(const std::string &configFilename, const std::string &sectionName="")
Give the file name, find out the corresponding section. Then call the main starter.
std::vector< MagneticFieldVolume > fMagneticFieldVolumes
A magnetic field volume structure to store field data and mesh.
TVector3 GetFieldAverageTransverseVector(TVector3 from, TVector3 to, Double_t dl=10., Int_t Nmax=0)
It returns the transverse component of the average magnetic field vector calculated along the line th...
TH2D * fHisto
A helper histogram to plot the field.
void SetError(std::string message="", bool print=true)
A metadata class may use this method to signal that something went wrong.
A base class for any REST metadata class.
Definition: TRestMetadata.h:70
endl_t RESTendl
Termination flag object for TRestStringOutput.
std::vector< Double_t > GetTransversalComponentAlongPath(TVector3 from, TVector3 to, Double_t dl=1., Int_t Nmax=0)
It returns a vector describing the transversal magnetic field component between from and to positions...
~TRestAxionMagneticField()
Default destructor.
void InitFromConfigFile()
Initialization of TRestAxionMagnetic field members through a RML file.
A basic class inheriting from TObject to help creating a node grid definition.
Definition: TRestMesh.h:38
MagneticFieldVolume * GetMagneticVolume(Int_t id)
It returns a pointer to the corresponding magnetic volume id.
TVector3 MoveByDistanceFast(TVector3 pos, TVector3 dir, Double_t d)
This method transports a position pos by a distance d in the direction defined by dir...
std::string GetParameter(std::string parName, TiXmlElement *e, TString defaultValue=PARAMETER_NOT_FOUND_STR)
Returns the value for the parameter named parName in the given section.
TCanvas * DrawHistogram(TString projection, TString Bcomp, Int_t volIndex=-1, Double_t step=-1, TString style="COLZ0", Double_t depth=-100010.0)
A method that creates a canvas where magnetic field map is drawn.
This namespace serves to define physics constants and other basic physical operations.
Definition: TRestPhysics.h:34