QwAnalysis
QwMagneticFieldConvert.cc
Go to the documentation of this file.
1 // Standard C and C++ headers
2 #include <iostream>
3 #include <string>
4 #include <sys/time.h>
5 
6 // ROOT headers
7 #include <TRandom3.h>
8 
9 // Qweak headers
10 #include "QwLog.h"
11 #include "QwOptions.h"
12 #include "QwParameterFile.h"
13 #include "QwMagneticField.h"
14 
15 
16 #define NSAMPLES 1000000
17 
18 int main (int argc, char* argv[])
19 {
20  // Search path for parameter files
21  QwParameterFile::AppendToSearchPath(getenv_safe_string("QW_PRMINPUT"));
22  QwParameterFile::AppendToSearchPath(getenv_safe_string("QWANALYSIS") + "/Tracking/prminput");
23 
24  // Option definitions of standard objects
25  QwMagneticField::DefineOptions(gQwOptions);
26 
27  // Custom options of this executable
28  gQwOptions.AddOptions("Magnetic field map actions")
29  ("convert",po::value<bool>()->default_bool_value(false),
30  "Convert field map from ASCII to binary");
31  gQwOptions.AddOptions("Magnetic field map actions")
32  ("timing",po::value<bool>()->default_bool_value(false),
33  "Time magnetic field reading routines");
34 
35  // Set command line
36  gQwOptions.SetCommandLine(argc, argv);
37 
38  // Get options
39  bool timing = gQwOptions.GetValue<bool>("timing");
40  bool convert = gQwOptions.GetValue<bool>("convert");
41  // Note: mapname is the mapfile without extension
42  std::string mapfile = gQwOptions.GetValue<std::string>("QwMagneticField.mapfile");
43  std::string mapname = mapfile.substr(0,mapfile.find_last_of("."));
44 
45  timeval time_start, time_finish;
46 
47 
48  // Initialize the field map
49  gettimeofday(&time_start, 0);
50  QwMagneticField* magneticfield = new QwMagneticField(gQwOptions,true);
51  gettimeofday(&time_finish, 0);
52  if (timing) {
53  int time_initialize_sec = ((int) time_finish.tv_sec - (int) time_start.tv_sec);
54  int time_initialize_usec = ((int) time_finish.tv_usec - (int) time_start.tv_usec);
55  if (time_initialize_usec < 0) { time_initialize_usec += 1000000; time_initialize_sec--; }
56  QwMessage << "Initialization: "
57  << time_initialize_sec << " sec and "
58  << time_initialize_usec << " usec" << QwLog::endl;
59  }
60 
61  // Set data reduction factor
62  //magneticfield->SetDataReductionFactor(2);
63  //QwMessage << "Data reduction factor set to "
64  // << magneticfield->GetDataReductionFactor() << QwLog::endl;
65 
66 
67  // Initialize common variables for all tests
68  double mean = 30.0;
69  double sigma = 3.0;
70  double field[3];
71  TRandom3 random;
72  int time_sampling_sec;
73  int time_sampling_usec;
74  double time_sampling;
75 
76  // Test point
77  double point[3];
78  // r = 100 cm, phi = 67.5 degrees, z = 100 cm
79  //point[0] = 38.2683 * Qw::cm; point[1] = 92.388 * Qw::cm; point[2] = 100.0 * Qw::cm;
80  // r = 100 cm, phi = 0 degrees, z = 100 cm
81  point[0] = 100.0 * Qw::cm; point[1] = 0.0 * Qw::cm; point[2] = 100.0 * Qw::cm;
82 
83  // Read the text field map
84  gettimeofday(&time_start, 0);
85  magneticfield->SetFilename(mapfile);
86  magneticfield->ReadFieldMap();
87  gettimeofday(&time_finish, 0);
88  if (timing) {
89  int time_reading_sec = ((int) time_finish.tv_sec - (int) time_start.tv_sec);
90  int time_reading_usec = ((int) time_finish.tv_usec - (int) time_start.tv_usec);
91  if (time_reading_usec < 0) { time_reading_usec += 1000000; time_reading_sec--; }
92  QwMessage << "Reading field map (text): "
93  << time_reading_sec << " sec and "
94  << time_reading_usec << " usec" << QwLog::endl;
95  }
96  magneticfield->GetCartesianFieldValue(point, field);
97  QwMessage << "Field at " << point[0]/Qw::cm << "," << point[1]/Qw::cm << "," << point[2]/Qw::cm << " cm "
98  << "is " << field[0]/Qw::kG << "," << field[1]/Qw::kG << "," << field[2]/Qw::kG << " kG"
99  << QwLog::endl;
100 
101  // Conversion from text field map to binary field map
102  if (convert) {
103 
104  // Write the binary field map
105  gettimeofday(&time_start, 0);
106  magneticfield->WriteBinaryFile(mapname + ".bin");
107  gettimeofday(&time_finish, 0);
108  if (timing) {
109  int time_writing_sec = ((int) time_finish.tv_sec - (int) time_start.tv_sec);
110  int time_writing_usec = ((int) time_finish.tv_usec - (int) time_start.tv_usec);
111  if (time_writing_usec < 0) { time_writing_usec += 1000000; time_writing_sec--; }
112  QwMessage << "Writing field map (binary): "
113  << time_writing_sec << " sec and "
114  << time_writing_usec << " usec" << QwLog::endl;
115  }
116 
117  // Delete old field map
118  gettimeofday(&time_start, 0);
119  delete magneticfield;
120  gettimeofday(&time_finish, 0);
121  if (timing) {
122  int time_destruction_sec = ((int) time_finish.tv_sec - (int) time_start.tv_sec);
123  int time_destruction_usec = ((int) time_finish.tv_usec - (int) time_start.tv_usec);
124  if (time_destruction_usec < 0) { time_destruction_usec += 1000000; time_destruction_sec--; }
125  QwMessage << "Destruction: "
126  << time_destruction_sec << " sec and "
127  << time_destruction_usec << " usec" << QwLog::endl;
128  }
129 
130  // Create new field map
131  gettimeofday(&time_start, 0);
132  magneticfield = new QwMagneticField(gQwOptions,true);
133  gettimeofday(&time_finish, 0);
134  if (timing) {
135  int time_initialize2_sec = ((int) time_finish.tv_sec - (int) time_start.tv_sec);
136  int time_initialize2_usec = ((int) time_finish.tv_usec - (int) time_start.tv_usec);
137  if (time_initialize2_usec < 0) { time_initialize2_usec += 1000000; time_initialize2_sec--; }
138  QwMessage << "Initialization: "
139  << time_initialize2_sec << " sec and "
140  << time_initialize2_usec << " usec" << QwLog::endl;
141  }
142 
143  // Read the binary field map
144  gettimeofday(&time_start, 0);
145  magneticfield->SetFilename(mapname + ".bin");
146  magneticfield->ReadFieldMap();
147  gettimeofday(&time_finish, 0);
148  if (timing) {
149  int time_reading2_sec = ((int) time_finish.tv_sec - (int) time_start.tv_sec);
150  int time_reading2_usec = ((int) time_finish.tv_usec - (int) time_start.tv_usec);
151  if (time_reading2_usec < 0) { time_reading2_usec += 1000000; time_reading2_sec--; }
152  QwMessage << "Reading field map (binary): "
153  << time_reading2_sec << " sec and "
154  << time_reading2_usec << " usec" << QwLog::endl;
155  }
156 
157  } // end of conversion
158 
159  // Multilinear interpolation
160  magneticfield->SetInterpolationMethod(kMultiLinear);
161  QwMessage << "Interpolation method set to multilinear, method "
162  << magneticfield->GetInterpolationMethod() << QwLog::endl;
163 
164  magneticfield->GetCartesianFieldValue(point, field);
165  QwMessage << "Field at " << point[0]/Qw::cm << "," << point[1]/Qw::cm << "," << point[2]/Qw::cm << " cm "
166  << "is " << field[0]/Qw::kG << "," << field[1]/Qw::kG << "," << field[2]/Qw::kG << " kG"
167  << QwLog::endl;
168 
169  if (timing) {
170  gettimeofday(&time_start, 0);
171  double rpoint[3];
172  for (int i = 0; i < NSAMPLES; i++) {
173  rpoint[0] = point[0] + random.Gaus(mean,sigma);
174  rpoint[1] = point[1] + random.Gaus(mean,sigma);
175  rpoint[2] = point[2] + random.Gaus(mean,sigma);
176  magneticfield->GetCartesianFieldValue(rpoint, field);
177  }
178  gettimeofday(&time_finish, 0);
179  time_sampling_sec = ((int) time_finish.tv_sec - (int) time_start.tv_sec);
180  time_sampling_usec = ((int) time_finish.tv_usec - (int) time_start.tv_usec);
181  if (time_sampling_usec < 0) { time_sampling_usec += 1000000; time_sampling_sec--; }
182  time_sampling = 1000000 * time_sampling_sec + time_sampling_usec;
183  QwMessage << "Sampling (total of " << NSAMPLES << " samples): "
184  << time_sampling_sec << " sec and "
185  << time_sampling_usec << " usec" << QwLog::endl;
186  QwMessage << "Sampling (per sample): " << time_sampling / NSAMPLES << " usec" << QwLog::endl;
187  }
188 
189  // Nearest-neighbor interpolation
190  magneticfield->SetInterpolationMethod(kNearestNeighbor);
191  QwMessage << "Interpolation method set to nearest-neighbor, method "
192  << magneticfield->GetInterpolationMethod() << QwLog::endl;
193 
194  magneticfield->GetCartesianFieldValue(point, field);
195  QwMessage << "Field at " << point[0]/Qw::cm << "," << point[1]/Qw::cm << "," << point[2]/Qw::cm << " cm "
196  << "is " << field[0]/Qw::kG << "," << field[1]/Qw::kG << "," << field[2]/Qw::kG << " kG"
197  << QwLog::endl;
198 
199  if (timing) {
200  gettimeofday(&time_start, 0);
201  double rpoint[3];
202  for (int i = 0; i < NSAMPLES; i++) {
203  rpoint[0] = point[0] + random.Gaus(mean,sigma);
204  rpoint[1] = point[1] + random.Gaus(mean,sigma);
205  rpoint[2] = point[2] + random.Gaus(mean,sigma);
206  magneticfield->GetCartesianFieldValue(rpoint,field);
207  }
208  gettimeofday(&time_finish, 0);
209  time_sampling_sec = ((int) time_finish.tv_sec - (int) time_start.tv_sec);
210  time_sampling_usec = ((int) time_finish.tv_usec - (int) time_start.tv_usec);
211  if (time_sampling_usec < 0) { time_sampling_usec += 1000000; time_sampling_sec--; }
212  time_sampling = 1000000 * time_sampling_sec + time_sampling_usec;
213  QwMessage << "Sampling (total of " << NSAMPLES << " samples): "
214  << time_sampling_sec << " sec and "
215  << time_sampling_usec << " usec" << QwLog::endl;
216  QwMessage << "Sampling (per sample): " << time_sampling / NSAMPLES << " usec" << QwLog::endl;
217  }
218 
219  // Delete the field map
220  delete magneticfield;
221 
222  return 0;
223 }
QwMessage
#define QwMessage
Predefined log drain for regular messages.
Definition: QwLog.h:50
QwParameterFile::AppendToSearchPath
static void AppendToSearchPath(const TString &searchdir)
Add a directory to the search path.
Definition: QwParameterFile.cc:35
QwMagneticField
Magnetic field map object.
Definition: QwMagneticField.h:50
default_bool_value
#define default_bool_value(b)
Definition: QwOptions.h:51
QwMagneticField::SetFilename
void SetFilename(const std::string &filename)
Set the filename.
Definition: QwMagneticField.h:77
gQwOptions
QwOptions gQwOptions
Definition: QwOptions.cc:27
QwOptions::AddOptions
po::options_description_easy_init AddOptions(const std::string &blockname="Specialized options")
Add an option to a named block or create new block.
Definition: QwOptions.h:164
Qw::kG
static const double kG
Definition: QwUnits.h:113
QwOptions::GetValue
T GetValue(const std::string &key)
Get a templated value.
Definition: QwOptions.h:240
QwMagneticField::SetInterpolationMethod
void SetInterpolationMethod(const EQwInterpolationMethod method)
Set interpolation method.
Definition: QwMagneticField.h:188
QwMagneticField::GetCartesianFieldValue
void GetCartesianFieldValue(const double point_xyz[3], double field_xyz[3]) const
Get the cartesian components of the field value.
Definition: QwMagneticField.h:124
QwLog.h
A logfile class, based on an identical class in the Hermes analyzer.
NSAMPLES
#define NSAMPLES
Definition: QwMagneticFieldConvert.cc:16
QwLog::endl
static std::ostream & endl(std::ostream &)
End of the line.
Definition: QwLog.cc:299
getenv_safe_string
const std::string getenv_safe_string(const char *name)
Definition: QwOptions.h:37
QwOptions.h
An options class which parses command line, config file and environment.
QwMagneticField::DefineOptions
static void DefineOptions(QwOptions &options)
Define command line and config file options.
Definition: QwMagneticField.cc:79
QwMagneticField::WriteBinaryFile
bool WriteBinaryFile(const std::string &fieldmap) const
Write a binary field map.
Definition: QwMagneticField.h:176
QwMagneticField::GetInterpolationMethod
EQwInterpolationMethod GetInterpolationMethod() const
Get interpolation method.
Definition: QwMagneticField.h:192
main
int main(int argc, char **argv)
Definition: QwRoot.cc:20
QwOptions::SetCommandLine
void SetCommandLine(int argc, char *argv[], bool default_config_file=true)
Set the command line arguments.
Definition: QwOptions.cc:112
Qw::cm
static const double cm
Length units: base unit is mm.
Definition: QwUnits.h:61
QwMagneticField::ReadFieldMap
bool ReadFieldMap()
Read a field map.
Definition: QwMagneticField.cc:176