FairRoot/PandaRoot
visualize_fieldmaps.cxx
Go to the documentation of this file.
1 // visualization of field lines in the Panda field map
2 
3 #include<TString.h>
4 #include<TStopwatch.h>
5 #include<TROOT.h>
6 #include<FairModule.h>
7 #include<PndCave.h>
8 #include<PndMagnet.h>
9 #include<PndPipe.h>
10 #include<PndMvdDetector.h>
11 #include<FairRunSim.h>
12 #include<PndLmdDetector.h>
13 #include<FairPrimaryGenerator.h>
14 #include<FairBoxGenerator.h>
15 #include<PndMultiField.h>
16 #include<PndTransMap.h>
17 #include<PndDipoleMap.h>
18 #include<PndSolenoidMap.h>
19 #include<FairParRootFileIo.h>
20 #include<PndMultiFieldPar.h>
21 #include<map>
22 #include<vector>
23 #include<fstream>
24 #include<iostream>
25 #include<iomanip>
26 #include<sstream>
27 #include<TH1.h>
28 #include<TTree.h>
29 #include<TFile.h>
30 #include<TCanvas.h>
31 #include<TVector3.h>
32 #include "TPolyLine3D.h"
33 #include "TView.h"
34 #include "TCanvas.h"
35 #include "TRandom.h"
36 #include<TApplication.h>
37 #include "TGLViewer.h"
38 #include <TGeoManager.h>
39 #include "TGeoMaterial.h"
40 #include "TGeoMedium.h"
41 #include "TGeoVolume.h"
42 #include "TString.h"
43 #include "TSystem.h"
44 #include "FairGeoLoader.h"
45 #include "FairGeoInterface.h"
46 #include "FairGeoMedia.h"
47 #include "FairGeoBuilder.h"
48 #include "FairModule.h"
49 
50 using namespace std;
51 
52 vector<string> fieldmapnames;
53 
54 // read fem calculated fieldmap and transform it
55 // fieldtype 1: no symetry
56 // fieldtype 2: Solenoidfieled
57 // fieldtype 3: Dipolefield
58 // fieldtype 4: Transition region
59 // filename_out is the pandaroot filename without the extenstion
60 // in contrast to filename_in having the extension
61 void visualize_fieldmaps() {
62  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
63 
64  cout << " setup of visualization " << endl;
65  TApplication* myapp = NULL;
66  TCanvas *canvas_fieldlines = NULL;
67  myapp = new TApplication("myapp",0,0);
68 
69 
70  TString vmcWorkdir = gSystem->Getenv("VMCWORKDIR");
71  // materials and media
72  FairGeoLoader* geoLoad = new FairGeoLoader("TGeo", "FairGeoLoader");
73  FairGeoInterface* geoFace = geoLoad->getGeoInterface();
74  geoFace->setMediaFile(vmcWorkdir + "/geometry/media_pnd.geo");
75  geoFace->readMedia();
76  geoFace->print();
77  FairGeoMedia* geoMedia = geoFace->getMedia();
78  FairGeoMedium *FairMediumAir = geoMedia->getMedium("air");
79  FairGeoBuilder* geoBuild = geoLoad->getGeoBuilder();
80  Int_t nmed=geoBuild->createMedium(FairMediumAir);
81  TGeoVolume *top_vol = gGeoManager->MakeBox("cave", gGeoManager->GetMedium("air"), 400, 400, 1000);
82  top_vol->SetVisibility(0);
83 
84  TFile* solenoid_file = new TFile("$VMCWORKDIR/geometry/FullSolenoid.root" ,"OPEN");
85  TGeoVolume * solenoid = (TGeoVolume*) solenoid_file->Get("topNode");
86  if (!solenoid){
87  cout << " 'Error: no solenoid file found to display " << endl;
88  } else {
89  solenoid->SetLineColor(12);
90  solenoid->SetTransparency(50);
91  //top_vol->AddNode(solenoid, 1);
92  }
93  TFile* beampipe_file = new TFile("$VMCWORKDIR/geometry/beampipe_201309.root" ,"OPEN");
94  TGeoVolume * beampipe = (TGeoVolume*) beampipe_file->Get("pipeassembly");
95  if (!beampipe){
96  cout << " 'Error: no beam pipe file found to display " << endl;
97  } else {
98  beampipe->SetTransparency(50);
99  top_vol->AddNode(beampipe, 1);
100  }
101  TFile* lmd_setup_file = new TFile("$VMCWORKDIR/geometry/Dipole.root" ,"OPEN");
102  TGeoVolume * lmd_setup = (TGeoVolume*) lmd_setup_file->Get("cave");
103  if (!lmd_setup){
104  cout << " 'Error: no lmd setup found to display " << endl;
105  } else {
106  //lmd_setup->SetLineColor(12);
107  lmd_setup->SetTransparency(50);
108  //top_vol->AddNode(lmd_setup, 1);
109  }
110  gGeoManager->SetTopVolume(top_vol);
111  top_vol->Draw("ogl");
112 
113  cout << " loading panda field map for 1.5 GeV/c beam momentum" << endl;
114 
115  PndMultiField *fField= new PndMultiField();
116 
117  PndTransMap *map_t = new PndTransMap("TransMap", "R", 1.5);
118  PndDipoleMap *map_d1 = new PndDipoleMap("DipoleMap1", "R", 1.5);
119  PndDipoleMap *map_d2 = new PndDipoleMap("DipoleMap2", "R", 1.5);
120  PndSolenoidMap *map_s1 = new PndSolenoidMap("SolenoidMap1", "R");
121  PndSolenoidMap *map_s2 = new PndSolenoidMap("SolenoidMap2", "R");
122  PndSolenoidMap *map_s3 = new PndSolenoidMap("SolenoidMap3", "R");
123  PndSolenoidMap *map_s4 = new PndSolenoidMap("SolenoidMap4", "R");
124 
125  fField->AddField(map_t);
126  fField->AddField(map_d1);
127  fField->AddField(map_d2);
128  fField->AddField(map_s1);
129  fField->AddField(map_s2);
130  fField->AddField(map_s3);
131  fField->AddField(map_s4);
132 
133  fField->Init();
134 
135  // define the valid area of the field map
136  // currently the dipole area is the smallest one
137  const double x_min = -30.;
138  const double x_max = 30.;
139  const double y_min = -30.;
140  const double y_max = 30.;
141  const double z_min = -172.;
142  const double z_max = 660.;
143 
144  if (0){
145  cout << " drawing the field bounding box" << endl;
146  double x[2];
147  double y[2];
148  double z[2];
149  TPolyLine3D* line;
150 
151  for (int ix1 = 0; ix1 < 2; ix1++){
152  if (ix1) x[0] = x_min; else x[0] = x_max;
153  for (int ix2 = 0; ix2 < 2; ix2++){
154  if (ix2) x[1] = x_min; else x[1] = x_max;
155  for (int iy1 = 0; iy1 < 2; iy1++){
156  if (iy1) y[0] = y_min; else y[0] = y_max;
157  for (int iy2 = 0; iy2 < 2; iy2++){
158  if (iy2) y[1] = y_min; else y[1] = y_max;
159  for (int iz1 = 0; iz1 < 2; iz1++){
160  if (iz1) z[0] = z_min; else z[0] = z_max;
161  for (int iz2 = 0; iz2 < 2; iz2++){
162  if (iz2) z[1] = z_min; else z[1] = z_max;
163  line = new TPolyLine3D(2, x, y, z);
164  line->SetLineColor(3);
165  line->Draw("same");
166  }
167  }
168  }
169  }
170  }
171  }
172  }
173 
174  // store the starting points
175  vector<double>xstart_prev;
176  vector<double>ystart_prev;
177  vector<double>zstart_prev;
178  vector<int> line_color;
179  // Initialize with a starting grid
180  // field lines will be drawn into both directions of the grid
181  if (xstart_prev.size()==0){
182  //const Int_t n = 100;
183  //create a raster for one plane in the middle of the solenoid
184  //to start the field lines from
185  // Check here the field strength in order to have same distances
186  // for same field strength
187  double pval[3] = {0.,0.,0.};
188  double Bval[3] = {0.,0.,0.};
189  // get the field strength perpendicular to the
190  // plane in the following loop
191  fField->GetFieldValue(pval, Bval);
192  const double B_strength_init = fabs(Bval[2]);
193  cout << " the solenoid field strength is " << B_strength_init;
194  // that strength has to to correspond to a step size of 0.2 cm
195  const double norm_step_size = 0.2;
196  // span a x-y-plane starting from the point with
197  // the strongest field
198  // a rotational symmetry of the field map is assumed in all 4 quarters
199  cout << " initializing solenoid grid for field lines " << endl;
200  for (double iy = -2.; iy < 2.; iy += norm_step_size){
201  for (double ix = -2.; ix < 2.; ix += norm_step_size){
202  xstart_prev.push_back(ix);
203  ystart_prev.push_back(iy);
204  zstart_prev.push_back(0.);
205  int color = 4;//12; // + floor(7./4.*r2);
206  line_color.push_back(color);
207  /*pval[0] = ix;
208  pval[1] = 0.;
209  pval[2] = 0.;
210  fField->GetFieldValue(pval, Bval);
211  double B_strength = fabs(Bval[2]);
212  int mod = B_strength/B_strength_init;
213  */
214  }
215  }
216  // take the middle of the dipole
217  const double dipole_middle_z = 475.;
218  const double iy = 0.;
219  // create starting points for the dipole
220  // start in the middle of the dipole and follow x
221  // assume in the first loop a y-z mirror plane symmetry
222  cout << " creating variables " << endl;
223  bool print_info = true;
224  bool** found_mod;
225  const int max_mods = 250;
226  found_mod = new bool*[max_mods];
227  for (int i = 0; i < max_mods; i++)
228  found_mod[i] = new bool[max_mods];
229  cout << " initializing dipole grid for field lines " << endl;
230  for (int dir_x = -1; dir_x < 2; dir_x += 2){
231  for (int dir_z = -1; dir_z < 2; dir_z += 2){
232  for (int i = 0; i < max_mods; i++){
233  for (int j = 0; j < max_mods; j++){
234  found_mod[i][j] = false;
235  }
236  }
237  for (double _ix = 0.; _ix < 200.; _ix += .1){
238  for (double _iz = 0.; _iz < 500.; _iz += .1){
239  // skip double entries for the lines on the symmetry axis
240  if ((_ix == 0 && dir_x == -1) || (_iz == 0 && dir_z == -1)) continue;
241  double ix = _ix*dir_x;
242  double iz = _iz*dir_z + dipole_middle_z;
243  pval[0] = ix;
244  pval[1] = iy;
245  pval[2] = iz;
246  // check the validity of the field map range
247  //cout << x_min << " " << pval[0] << " " << x_max << endl;
248  //cout << y_min << " " << pval[1] << " " << y_max << endl;
249  //cout << z_min << " " << pval[2] << " " << z_max << endl << endl;
250  if (!(x_min < pval[0] && pval[0] < x_max)||
251  !(y_min < pval[1] && pval[1] < y_max)||
252  !(z_min < pval[2] && pval[2] < z_max)){
253  continue;
254  }
255  // get the field strength perpendicular to the
256  // plane in the following loop
257  fField->GetFieldValue(pval, Bval);
258  if (print_info){
259  print_info = false;
260  cout << " dipole field strength to normalize to is " << Bval[1] << endl;
261  }
262  int mod_x = (int) floor(fabs(_ix/norm_step_size*Bval[1]/B_strength_init));
263  int mod_z = (int) floor(fabs(_iz/norm_step_size*Bval[1]/B_strength_init));
264  if (mod_x >= max_mods || mod_z >= max_mods || mod_x < 0 || mod_z < 0){
265  cout << " Grid is out of bounds for mod x " << mod_x << " mod z " << mod_z << endl;
266  continue;
267  }
268  if (!found_mod[mod_x][mod_z]){
269  found_mod[mod_x][mod_z] = true;
270  xstart_prev.push_back(ix);
271  ystart_prev.push_back(iy);
272  zstart_prev.push_back(iz);
273  int color = 4;
274  line_color.push_back(color);
275  }
276  }// loop over iz
277  }// loop over ix
278  }// loop over dir_z
279  }// loop over dir_x
280  cout << " deleting " << endl;
281  for (int i = 0; i < max_mods; i++)
282  delete[] found_mod[i];
283  delete[] found_mod;
284  }
285  // define here the smallest grid spacing
286  double d_x = 1.;
287  double d_y = 1.;
288  double d_z = 1.;
289  cout << " drawing field lines " << endl;
290  // draw field lines following the path of the field vectors
291  // follow the field lines in both directions
292  for (int direction = 1; direction > -2; direction -= 2){
293  for (unsigned int istartpoint = 0; istartpoint < xstart_prev.size(); istartpoint++){
294  int npoints = 10000;
295  double stepsize = d_x;
296  if (d_y < stepsize) stepsize = d_y;
297  if (d_z < stepsize) stepsize = d_z;
298  vector<double> xs;
299  vector<double> ys;
300  vector<double> zs;
301  xs.push_back(xstart_prev[istartpoint]);
302  ys.push_back(ystart_prev[istartpoint]);
303  zs.push_back(zstart_prev[istartpoint]);
304  //cout << endl;
305  for (int ipoint = 1; ipoint < npoints; ipoint++){
306  // get the previous position
307  double pval[3] = {xs[ipoint-1],ys[ipoint-1],zs[ipoint-1]};
308  // check the validity of the field map range
309  //cout << x_min << " " << pval[0] << " " << x_max << endl;
310  //cout << y_min << " " << pval[1] << " " << y_max << endl;
311  //cout << z_min << " " << pval[2] << " " << z_max << endl << endl;
312  if (!(x_min < pval[0] && pval[0] < x_max)||
313  !(y_min < pval[1] && pval[1] < y_max)||
314  !(z_min < pval[2] && pval[2] < z_max)){
315  //if (xs.size()>2)
316  //cout << " break due to out of range after " << ipoint << " points " << endl;
317  break;
318  }
319  // get the corresponding field value
320  // what is in fact the direction of the field line
321  double Bval[3];
322  fField->GetFieldValue(pval, Bval);
323  // scale to the stepsize and add to the position
324  double fieldstrength = sqrt(Bval[0]*Bval[0]+Bval[1]*Bval[1]+Bval[2]*Bval[2]);
325  pval[0] = pval[0] + stepsize/fieldstrength*Bval[0]*direction;
326  pval[1] = pval[1] + stepsize/fieldstrength*Bval[1]*direction;
327  pval[2] = pval[2] + stepsize/fieldstrength*Bval[2]*direction;
328  // store the next point
329  xs.push_back(pval[0]);
330  ys.push_back(pval[1]);
331  zs.push_back(pval[2]);
332  }
333  // create a line out of the points
334  if (xs.size() > 2){
335  //cout << " creating polyline " << endl;
336  TPolyLine3D* fieldline = new TPolyLine3D(xs.size(), &xs[0], &ys[0], &zs[0]);
337  fieldline->SetLineColor(line_color[istartpoint]);
338  fieldline->Draw("same");
339  // mirror it
340  for (unsigned int ipoint = 0; ipoint < xs.size(); ipoint++){
341  xs[ipoint] = -xs[ipoint];
342  }
343  }
344  }
345  }
346  myapp->Run();
347 
348 }
349 
350 int main(int argc, char **argv) {
351  visualize_fieldmaps();
352  return 0;
353 }
354 
fField
PndMultiField * fField
Definition: sim_emc_apd.C:97
main
int main(int argc, char **argv)
Definition: MvdMQ/run/runFileSampler.cxx:29
geoLoad
FairGeoLoader * geoLoad
Definition: createRootGeoFile.C:16
i
Int_t i
Definition: run_full.C:25
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
PndMultiField::GetFieldValue
void GetFieldValue(const Double_t point[3], Double_t *bField)
Definition: PndMultiField.cxx:257
map_d2
PndDipoleMap * map_d2
Definition: sim_hit_emc.C:111
map_s1
PndSolenoidMap * map_s1
Definition: sim_hit_emc.C:112
gGeoManager
TGeoManager * gGeoManager
Definition: createRootGeoFile.C:97
FairMediumAir
FairGeoMedium * FairMediumAir
Definition: createRootGeoFile.C:85
geoMedia
FairGeoMedia * geoMedia
Definition: createRootGeoFile.C:22
map_s3
PndSolenoidMap * map_s3
Definition: sim_hit_emc.C:114
z
Double_t z
Definition: createRootFscGeometryFile.C:146
visualize_fieldmaps
void visualize_fieldmaps()
Definition: visualize_fieldmaps.cxx:61
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
map_s2
PndSolenoidMap * map_s2
Definition: sim_hit_emc.C:113
x
Double_t x
Definition: createRootFscGeometryFile.C:159
PndMultiField::AddField
void AddField(FairField *field)
Definition: PndMultiField.cxx:199
map_d1
PndDipoleMap * map_d1
Definition: sim_hit_emc.C:110
map_t
PndTransMap * map_t
Definition: sim_hit_emc.C:109
geoFace
FairGeoInterface * geoFace
Definition: createRootGeoFile.C:17
y
Double_t y
Definition: createRootFscGeometryFile.C:159
fieldmapnames
vector< string > fieldmapnames
Definition: interpolate_fieldmaps.cxx:41
geoBuild
FairGeoBuilder * geoBuild
Definition: createRootGeoFile.C:23
map_s4
PndSolenoidMap * map_s4
Definition: sim_hit_emc.C:115
Generated on Wed Apr 3 2019 08:44:14 for FairRoot/PandaRoot by   doxygen 1.8.5