src/graphics/DrawBField.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : M.Frank
0011 //
0012 //==========================================================================
0013
0014 // Framework includes
0015 #include <DD4hep/Printout.h>
0016 #include <DD4hep/FieldTypes.h>
0017 #include <DD4hep/MatrixHelpers.h>
0018 #include <DD4hep/DetFactoryHelper.h>
0019 #include <cmath>
0020
0021 using namespace dd4hep;
0022
0023 #include <TH2.h>
0024 #include <TPad.h>
0025 #include <TRint.h>
0026 #include <TAxis.h>
0027 #include <TStyle.h>
0028 #include <TArrow.h>
0029 #include <TCanvas.h>
0030 #include <THistPainter.h>
0031
0032 #include "Hoption.h"
0033 #include "Hparam.h"
0034
0035 extern Hoption_t Hoption;
0036 extern Hparam_t  Hparam;
0037
0038
0039 namespace  {
0040   void help_bfield(int argc, char** argv)   {
0041     std::cout <<
0042       "Usage: DD4hep_DrawBField -arg [-arg]                                 \n"
0043       "  The plugin implicitly assumes that the geometry is already loaded. \n"
0044       "     -area  <string>  Define area for which the B-field will be drawn\n"
0045       "     -nx:   <number>  Number of bins in X direction                  \n"
0046       "     -ny:   <number>  Number of bins in Y direction                  \n"
0047       "     -dx:   <range>   Definition of the range in X                   \n"
0048       "     -dy:   <range>   Definition of the range in Y                   \n"
0049       "     -dz:   <range>   Definition of the range in Z                   \n"
0050       "     Range definition as tuple: min,max                              \n"
0051       "                 min:  lower boundary                                \n"
0052       "                 max:  upper boundary                                \n"
0053       "\tArguments given: " << arguments(argc,argv) << std::endl << std::flush;
0054     ::exit(EINVAL);
0055   }
0056
0057   struct range_t  {
0058     double rmin=0e0, rmax=0e0;
0059   };
0060   struct field_t  {
0061     Position  position;
0062     Direction bfield;
0063   };
0064
0065   range_t get_range(const std::string& s, int argc, char** argv)   {
0066     int ival=0;
0067     range_t range;
0068     std::string val;
0069     for(std::size_t i=0; i<s.length()+1; ++i)   {
0070       char c = s.c_str()[i];
0071       if ( c == ',' || i == s.length() )  {
0072     switch(ival)   {
0073     case 0:
0074       range.rmin = _toDouble(val);
0075       break;
0076     case 1:
0077       range.rmax = _toDouble(val);
0078       break;
0079     default:
0080       except("","+++ Too many value for range descriptor provided: %s", s.c_str());
0081       break;
0082     }
0083     val = "";
0084     ++ival;
0085     continue;
0086       }
0087       val += c;
0088     }
0089     if ( ival != 2 )   {
0090       help_bfield(argc, argv);
0091     }
0092     return range;
0093   }
0094
0095   class MyHistPainter : public THistPainter {
0096   public:
0097     using THistPainter::THistPainter;
0098     TH2* paintArrows(TH2* histo, const std::vector<field_t>& points)   {
0099       fYaxis->SetTitle("Y coordinate [cm]");
0100       fXaxis->SetTitle("X coordinate [cm]");
0101       for(const auto& p : points)    {
0102     double strength2 = p.bfield.mag2();
0103     if ( strength2 > 1e-8 )   {
0104       double strength = sqrt(p.bfield.X()*p.bfield.X()+p.bfield.Y()*p.bfield.Y());
0105       int ix = fXaxis->FindBin(p.position.X());
0106       int iy = fYaxis->FindBin(p.position.Y());
0107       double xlo = fXaxis->GetBinLowEdge(ix);
0108       double xhi = fXaxis->GetBinLowEdge(ix+1);
0109       double ylo = fYaxis->GetBinLowEdge(iy);
0110       double yhi = fYaxis->GetBinLowEdge(iy+1);
0111       auto norm_field = p.bfield * ((xhi-xlo) / strength);
0112       double x1  = xlo + (xhi-xlo)/2.0;
0113       double x2  = x1 + norm_field.X();
0114       double y1  = ylo + (yhi-ylo)/2.0;
0115       double y2  = y1 + norm_field.Y();
0116       auto* arrow = new TArrow(x1, y1, x2, y2, 0.015, "|>");
0117       arrow->SetAngle(25);
0118       arrow->SetFillColor(kRed);
0119       arrow->SetLineColor(kRed);
0120       arrow->Draw();
0121       histo->Fill(p.position.X(), p.position.X(), strength);
0122       ::fprintf(stdout, "Arrow %+15.8e  %+15.8e  %+15.8e  %+15.8e    %+15.8e  %+15.8e\n",
0123             x1, y1, x2, y2, x2-x1, y2-y1);
0124       //delete arrow;
0125     }
0126       }
0127       PaintPalette();
0128       return histo;
0129     }
0130   };
0131
0132   /// Basic entry point to interprete an XML document
0133   /**
0134    *  - The file URI to be opened
0135    *    is passed as first argument to the call.
0136    *  - The processing hint (build type) is passed as optional
0137    *    second argument.
0138    *
0139    *  Factory: DD4hep_CompactLoader
0140    *
0141    *  \author  M.Frank
0142    *  \version 1.0
0143    *  \date    01/04/2014
0144    */
0145   static long draw_bfield(Detector& description, int argc, char** argv) {
0146     std::vector<std::string> srange_x, srange_y, srange_z;
0147     std::size_t nbin_x = 100, nbin_y = 100, nbin_z = 1;
0148     std::string output;
0149     double z_value = 0e0;
0150     bool draw = false;
0151
0152     printout(WARNING,"DrawBField","This is a TEST. It does not function properly!");
0153
0154     for(int i = 0; i < argc && argv[i]; ++i)  {
0155       if ( 0 == ::strncmp("-rx",argv[i],4) )
0156     srange_x.push_back(argv[++i]);
0157       else if ( 0 == ::strncmp("-ry",argv[i],3) )
0158     srange_y.push_back(argv[++i]);
0159       else if ( 0 == ::strncmp("-rz",argv[i],3) )
0160     srange_z.push_back(argv[++i]);
0161       else if ( 0 == ::strncmp("-nx",argv[i],3) )
0162     nbin_x = _toULong(argv[++i]);
0163       else if ( 0 == ::strncmp("-ny",argv[i],3) )
0164     nbin_y = _toULong(argv[++i]);
0165       else if ( 0 == ::strncmp("-nz",argv[i],3) )
0166     nbin_z = _toULong(argv[++i]);
0167       else if ( 0 == ::strncmp("-z",argv[i],2) )
0168     z_value = _toDouble(argv[++i]);
0169       else if ( 0 == ::strncmp("-draw",argv[i],4) )
0170         draw = true;
0171       else if ( 0 == ::strncmp("-output",argv[i],4) )
0172         output = argv[++i];
0173       else if ( 0 == ::strncmp("-help",argv[i],2) )
0174     help_bfield(argc, argv);
0175     }
0176     if ( srange_x.empty() || srange_y.empty() )   {
0177       help_bfield(argc, argv);
0178     }
0179     std::vector<range_t> range_x, range_y, range_z;
0180     for(const auto& r : srange_x)  {
0181       range_t rg = get_range(r, argc, argv);
0182       range_x.push_back(rg);
0183     }
0184     for(const auto& r : srange_y)  {
0185       range_t rg = get_range(r, argc, argv);
0186       range_y.push_back(rg);
0187     }
0188     for(const auto& r : srange_z)  {
0189       range_t rg = get_range(r, argc, argv);
0190       range_z.push_back(rg);
0191     }
0192
0193     double ma = std::numeric_limits<double>::max();
0194     range_t envelope_x  { ma, -ma };
0195     range_t envelope_y  { ma, -ma };
0196     range_t envelope_z  { ma, -ma };
0197     for( const auto& range : range_x )   {
0198       envelope_x.rmin = std::min(range.rmin, envelope_x.rmin);
0199       envelope_x.rmax = std::max(range.rmax, envelope_x.rmax);
0200     }
0201     for( const auto& range : range_y )   {
0202       envelope_y.rmin = std::min(range.rmin, envelope_y.rmin);
0203       envelope_y.rmax = std::max(range.rmax, envelope_y.rmax);
0204     }
0205     for( const auto& range : range_z )   {
0206       envelope_z.rmin = std::min(range.rmin, envelope_z.rmin);
0207       envelope_z.rmax = std::max(range.rmax, envelope_z.rmax);
0208     }
0209
0210     double dx = (envelope_x.rmax - envelope_x.rmin) / double(nbin_x);
0211     double dy = (envelope_y.rmax - envelope_y.rmin) / double(nbin_y);
0212     double dz = nbin_z == 1 ? 0e0 : (envelope_z.rmax - envelope_z.rmin) / double(nbin_z);
0213     printf("Range(x) min:%4ld bins %+15.8e cm max:%+15.8e cm dx:%+15.8e cm\n",
0214        nbin_x, envelope_x.rmin/cm, envelope_x.rmax/cm, dx/cm);
0215     printf("Range(y) min:%4ld bins %+15.8e cm max:%+15.8e cm dx:%+15.8e cm\n",
0216        nbin_y, envelope_y.rmin/cm, envelope_y.rmax/cm, dy/cm);
0217     if ( nbin_z > 1 ) printf("Range(z) min:%4ld bins %+15.8e cm max:%+15.8e cm dx:%+15.8e cm\n",
0218                  nbin_z, envelope_z.rmin/cm, envelope_z.rmax/cm, dz/cm);
0219
0220     FILE* out_file = ::fopen(output.empty() ? "/dev/null" : output.c_str(), "w");
0221     ::fprintf(out_file,"#######################################################################################################\n");
0222     ::fprintf(out_file,"      x[cm]            y[cm]            z[cm]          Bx[Tesla]        By[Tesla]        Bz[Tesla]     \n");
0223     std::vector<field_t> field_values;
0224     for( std::size_t i = 0; i < nbin_x; ++i )   {
0225       float x = envelope_x.rmin + double(i)*dx + dx/2e0;
0226       for( std::size_t j = 0; j < nbin_y; ++j )   {
0227     float y = envelope_y.rmin + double(j)*dy + dy/2e0;
0228     for( std::size_t k = 0; k < nbin_z; ++k )   {
0229       float z = nbin_z == 1 ? z_value : envelope_z.rmin + double(k)*dz + dz/2e0;
0230       field_t value;
0231       value.position = { x, y, z };
0232       value.bfield   = { 0e0, 0e0, 0e0 };
0233       value.bfield = description.field().magneticField(value.position);
0234       ::fprintf(out_file, " %+15.8e  %+15.8e  %+15.8e  %+15.8e  %+15.8e  %+15.8e\n",
0235          value.position.X()/cm, value.position.Y()/cm,  value.position.Z()/cm,
0236          value.bfield.X()/dd4hep::tesla, value.bfield.Y()/dd4hep::tesla, value.bfield.Z()/dd4hep::tesla);
0237       field_values.emplace_back(value);
0238     }
0239       }
0240     }
0241     ::fclose(out_file);
0242     if ( draw )   {
0243       if ( 0 == gApplication )  {
0244     std::pair<int, char**> a(argc,argv);
0245     gApplication = new TRint("DD4hepRint", &a.first, a.second);
0246     printout(INFO,"DD4hepRint","++ Created ROOT interpreter instance for DD4hepUI.");
0247       }
0248       auto* histo = new TH2F("Bfield", "B-Field strength in Tesla",
0249                 nbin_x, envelope_x.rmin, envelope_x.rmax,
0250                 nbin_y, envelope_y.rmin, envelope_y.rmax);
0251       MyHistPainter paint;
0252       paint.SetHistogram(histo);
0253       TCanvas* c1 = new TCanvas("B-Field");
0254       c1->SetWindowSize(1000,1000);
0255       //paint.Paint();
0256       histo->SetStats(kFALSE);
0257       //histo->Draw();
0258       paint.paintArrows(histo, field_values);
0259       histo->Draw("COLZ SAME");
0260       gPad->SetGridx();
0261       gPad->SetGridy();
0262       if ( !gApplication->IsRunning() )  {
0263     gApplication->Run();
0264       }
0265     }
0266     return 1;
0267   }
0268 }
0269 DECLARE_APPLY(DD4hep_DrawBField,draw_bfield)