00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #include "PSCalib/GeometryObject.h"
00017
00018
00019
00020
00021 #include <iostream>
00022 #include <sstream>
00023 #include <iomanip>
00024 #include <cmath>
00025 #include <cstring>
00026
00027
00028
00029
00030 #include "MsgLogger/MsgLogger.h"
00031 #include "PSCalib/GlobalMethods.h"
00032
00033
00034
00035
00036
00037 using namespace std;
00038
00039
00040
00041
00042
00043 namespace PSCalib {
00044
00045
00046
00047
00048 GeometryObject::GeometryObject ( std::string pname,
00049 unsigned pindex,
00050 std::string oname,
00051 unsigned oindex,
00052 double x0,
00053 double y0,
00054 double z0,
00055 double rot_z,
00056 double rot_y,
00057 double rot_x,
00058 double tilt_z,
00059 double tilt_y,
00060 double tilt_x
00061 )
00062 : m_pname (pname)
00063 , m_pindex (pindex)
00064 , m_oname (oname)
00065 , m_oindex (oindex)
00066 , m_x0 (x0)
00067 , m_y0 (y0)
00068 , m_z0 (z0)
00069 , m_rot_z (rot_z)
00070 , m_rot_y (rot_y)
00071 , m_rot_x (rot_x)
00072 , m_tilt_z (tilt_z)
00073 , m_tilt_y (tilt_y)
00074 , m_tilt_x (tilt_x)
00075 , m_do_tilt(true)
00076 , m_mbits(0377)
00077 , m_size(0)
00078 , p_xarr(0)
00079 , p_yarr(0)
00080 , p_zarr(0)
00081 , p_aarr(0)
00082 , p_marr(0)
00083 {
00084 const unsigned print_bits=0;
00085 m_seggeom = PSCalib::SegGeometryStore::Create(m_oname, print_bits);
00086 m_parent = shpGO();
00087 v_list_of_children.clear();
00088 }
00089
00090
00091
00092
00093 GeometryObject::~GeometryObject ()
00094 {
00095 if (m_seggeom) delete m_seggeom;
00096
00097 deallocate_memory();
00098 }
00099
00100
00101 void GeometryObject::deallocate_memory()
00102 {
00103 if (p_xarr) delete [] p_xarr; p_xarr=0;
00104 if (p_yarr) delete [] p_yarr; p_yarr=0;
00105 if (p_zarr) delete [] p_zarr; p_zarr=0;
00106 if (p_aarr) delete [] p_aarr; p_aarr=0;
00107 if (p_marr) delete [] p_marr; p_marr=0;
00108 }
00109
00110
00111 std::string GeometryObject::string_geo()
00112 {
00113 std::stringstream ss;
00114 ss << "parent:" << std::setw(10) << std::right << m_pname
00115 << " pind:" << std::setw(2) << m_pindex
00116 << " geo:" << std::setw(10) << m_oname
00117 << " oind:" << std::setw(2) << m_oindex << std::right << std::fixed
00118 << " x0:" << std::setw(8) << std::setprecision(0) << m_x0
00119 << " y0:" << std::setw(8) << std::setprecision(0) << m_y0
00120 << " z0:" << std::setw(8) << std::setprecision(0) << m_z0
00121 << " rot_z:" << std::setw(6) << std::setprecision(1) << m_rot_z
00122 << " rot_y:" << std::setw(6) << std::setprecision(1) << m_rot_y
00123 << " rot_x:" << std::setw(6) << std::setprecision(1) << m_rot_x
00124 << " tilt_z:" << std::setw(8) << std::setprecision(5) << m_tilt_z
00125 << " tilt_y:" << std::setw(8) << std::setprecision(5) << m_tilt_y
00126 << " tilt_x:" << std::setw(8) << std::setprecision(5) << m_tilt_x;
00127 return ss.str();
00128 }
00129
00130
00131 std::string GeometryObject::str_data()
00132 {
00133 std::stringstream ss;
00134 ss << std::setw(11) << std::left << m_pname
00135 << " " << std::setw(3) << std::right << m_pindex
00136 << " " << std::setw(11) << std::left << m_oname
00137 << " " << std::setw(3) << std::right << m_oindex << std::fixed
00138 << " " << std::setw(8) << std::setprecision(0) << m_x0
00139 << " " << std::setw(8) << std::setprecision(0) << m_y0
00140 << " " << std::setw(8) << std::setprecision(0) << m_z0
00141 << " " << std::setw(6) << std::setprecision(1) << m_rot_z
00142 << " " << std::setw(6) << std::setprecision(1) << m_rot_y
00143 << " " << std::setw(6) << std::setprecision(1) << m_rot_x
00144 << " " << std::setw(8) << std::setprecision(5) << m_tilt_z
00145 << " " << std::setw(8) << std::setprecision(5) << m_tilt_y
00146 << " " << std::setw(8) << std::setprecision(5) << m_tilt_x;
00147 return ss.str();
00148 }
00149
00150
00151 void GeometryObject::print_geo()
00152 {
00153
00154 MsgLog(name(), info, string_geo());
00155 }
00156
00157
00158 std::string GeometryObject::string_geo_children()
00159 {
00160 std::stringstream ss;
00161 ss << "parent:" << std::setw(10) << std::left << m_pname
00162 << " i:" << std::setw(2) << m_pindex
00163 << " geo:" << std::setw(10) << m_oname
00164 << " i:" << std::setw(2) << m_oindex
00165 << " #children:" << v_list_of_children.size();
00166 for(std::vector<shpGO>::iterator it = v_list_of_children.begin();
00167 it != v_list_of_children.end(); ++ it) {
00168 ss << " " << (*it)->get_geo_name() << " " << (*it)->get_geo_index();
00169 }
00170 return ss.str();
00171 }
00172
00173
00174 void GeometryObject::print_geo_children()
00175 {
00176
00177 MsgLog(name(), info, string_geo_children());
00178 }
00179
00180
00181 void GeometryObject::transform_geo_coord_arrays(const double* X,
00182 const double* Y,
00183 const double* Z,
00184 const unsigned size,
00185 double* Xt,
00186 double* Yt,
00187 double* Zt,
00188 const bool do_tilt)
00189 {
00190
00191 double angle_x = (do_tilt) ? m_rot_x + m_tilt_x : m_rot_x;
00192 double angle_y = (do_tilt) ? m_rot_y + m_tilt_y : m_rot_y;
00193 double angle_z = (do_tilt) ? m_rot_z + m_tilt_z : m_rot_z;
00194
00195
00196 double* X1 = new double [size];
00197 double* Y1 = new double [size];
00198 double* Z2 = new double [size];
00199
00200
00201 rotation(X, Y, size, angle_z, X1, Y1);
00202 rotation(Z, X1, size, angle_y, Z2, Xt);
00203 rotation(Y1, Z2, size, angle_x, Yt, Zt);
00204
00205
00206 for(unsigned i=0; i<size; ++i) {
00207 Xt[i] += m_x0;
00208 Yt[i] += m_y0;
00209 Zt[i] += m_z0;
00210 }
00211
00212
00213 delete [] X1;
00214 delete [] Y1;
00215 delete [] Z2;
00216 }
00217
00218
00219 unsigned GeometryObject::get_size_geo_array()
00220 {
00221 if(m_seggeom) return m_seggeom -> size();
00222
00223 unsigned size=0;
00224 for(std::vector<shpGO>::iterator it = v_list_of_children.begin();
00225 it != v_list_of_children.end(); ++it) {
00226 size += (*it)->get_size_geo_array();
00227 }
00228 return size;
00229 }
00230
00231
00232 double GeometryObject::get_pixel_scale_size()
00233 {
00234 if(m_seggeom) return m_seggeom -> pixel_scale_size();
00235
00236 double pixel_scale_size=1;
00237
00238 for(std::vector<shpGO>::iterator it = v_list_of_children.begin();
00239 it != v_list_of_children.end(); ++it) {
00240 pixel_scale_size = (*it)->get_pixel_scale_size();
00241 break;
00242 }
00243 return pixel_scale_size;
00244 }
00245
00246
00247 void GeometryObject::get_pixel_coords(const double*& X, const double*& Y, const double*& Z, unsigned& size,
00248 const bool do_tilt, const bool do_eval)
00249 {
00250
00251 if(p_xarr==0 || do_tilt != m_do_tilt || do_eval) evaluate_pixel_coords(do_tilt, do_eval);
00252 X = p_xarr;
00253 Y = p_yarr;
00254 Z = p_zarr;
00255 size = m_size;
00256 }
00257
00258
00259 void GeometryObject::get_pixel_areas(const double*& areas, unsigned& size)
00260 {
00261 if(p_aarr==0) evaluate_pixel_coords();
00262 areas = p_aarr;
00263 size = m_size;
00264 }
00265
00266
00267 void GeometryObject::get_pixel_mask(const int*& mask, unsigned& size, const unsigned& mbits)
00268 {
00269 if(mbits != m_mbits or p_marr==0) { m_mbits = mbits; evaluate_pixel_coords(); }
00270 mask = p_marr;
00271 size = m_size;
00272 }
00273
00274
00275 void GeometryObject::evaluate_pixel_coords(const bool do_tilt, const bool do_eval)
00276 {
00277 m_do_tilt = do_tilt;
00278
00279 unsigned size = get_size_geo_array();
00280
00281 if(p_xarr==0 || size != m_size) {
00282
00283 m_size = size;
00284
00285 this->deallocate_memory();
00286
00287 p_xarr = new double [m_size];
00288 p_yarr = new double [m_size];
00289 p_zarr = new double [m_size];
00290 p_aarr = new double [m_size];
00291 p_marr = new int [m_size];
00292 }
00293
00294 if(m_seggeom) {
00295
00296 const double* x_arr = m_seggeom -> pixel_coord_array (SG::AXIS_X);
00297 const double* y_arr = m_seggeom -> pixel_coord_array (SG::AXIS_Y);
00298 const double* z_arr = m_seggeom -> pixel_coord_array (SG::AXIS_Z);
00299 const double* a_arr = m_seggeom -> pixel_area_array();
00300 const int* m_arr = m_seggeom -> pixel_mask_array(m_mbits);
00301
00302 transform_geo_coord_arrays(x_arr, y_arr, z_arr, m_size, p_xarr, p_yarr, p_zarr, do_tilt);
00303 std::memcpy(&p_aarr[0], a_arr, m_size*sizeof(double));
00304 std::memcpy(&p_marr[0], m_arr, m_size*sizeof(int));
00305 return;
00306 }
00307
00308 unsigned ibase=0;
00309 unsigned ind=0;
00310 for(std::vector<shpGO>::iterator it = v_list_of_children.begin();
00311 it != v_list_of_children.end(); ++it, ++ind) {
00312
00313 if((*it)->get_geo_index() != ind) {
00314 std::stringstream ss;
00315 ss << "WARNING! Geometry object:" << (*it)->get_geo_name() << ":" << (*it)->get_geo_index()
00316 << " has non-consequtive index at reconstruction: " << ind << '\n';
00317
00318 MsgLog(name(), warning, ss.str());
00319 }
00320
00321 const double* pXch;
00322 const double* pYch;
00323 const double* pZch;
00324 const double* pAch;
00325 const int* pMch;
00326 unsigned sizech;
00327
00328 (*it)->get_pixel_coords(pXch, pYch, pZch, sizech, do_tilt, do_eval);
00329 transform_geo_coord_arrays(pXch, pYch, pZch, sizech, &p_xarr[ibase], &p_yarr[ibase], &p_zarr[ibase], do_tilt);
00330
00331 (*it)->get_pixel_areas(pAch, sizech);
00332 std::memcpy(&p_aarr[ibase], pAch, sizech*sizeof(double));
00333
00334 (*it)->get_pixel_mask(pMch, sizech, m_mbits);
00335 std::memcpy(&p_marr[ibase], pMch, sizech*sizeof(int));
00336
00337 ibase += sizech;
00338
00339 (*it)->deallocate_memory();
00340 }
00341
00342 if(ibase == PSCalib::SIZE2X2 && m_oname == "CSPAD2X2:V1") {
00343
00344
00345
00346 two2x1ToData2x2<double>(p_xarr);
00347 two2x1ToData2x2<double>(p_yarr);
00348 two2x1ToData2x2<double>(p_zarr);
00349 two2x1ToData2x2<double>(p_aarr);
00350 two2x1ToData2x2<int> (p_marr);
00351 }
00352 }
00353
00354
00355
00356 void GeometryObject::get_geo_pars( double& x0,
00357 double& y0,
00358 double& z0,
00359 double& rot_z,
00360 double& rot_y,
00361 double& rot_x,
00362 double& tilt_z,
00363 double& tilt_y,
00364 double& tilt_x
00365 )
00366 {
00367 x0 = m_x0;
00368 y0 = m_y0;
00369 z0 = m_z0;
00370 rot_z = m_rot_z;
00371 rot_y = m_rot_y;
00372 rot_x = m_rot_x;
00373 tilt_z = m_tilt_z;
00374 tilt_y = m_tilt_y;
00375 tilt_x = m_tilt_x;
00376 }
00377
00378
00379
00380 void GeometryObject::set_geo_pars( const double& x0,
00381 const double& y0,
00382 const double& z0,
00383 const double& rot_z,
00384 const double& rot_y,
00385 const double& rot_x,
00386 const double& tilt_z,
00387 const double& tilt_y,
00388 const double& tilt_x
00389 )
00390 {
00391 m_x0 = x0;
00392 m_y0 = y0;
00393 m_z0 = z0;
00394 m_rot_z = rot_z;
00395 m_rot_y = rot_y;
00396 m_rot_x = rot_x;
00397 m_tilt_z = tilt_z;
00398 m_tilt_y = tilt_y;
00399 m_tilt_x = tilt_x;
00400 }
00401
00402
00403
00404 void GeometryObject::move_geo( const double& dx,
00405 const double& dy,
00406 const double& dz
00407 )
00408 {
00409 m_x0 += dx;
00410 m_y0 += dy;
00411 m_z0 += dz;
00412 }
00413
00414
00415
00416 void GeometryObject::tilt_geo( const double& dt_x,
00417 const double& dt_y,
00418 const double& dt_z
00419 )
00420 {
00421 m_tilt_z += dt_z;
00422 m_tilt_y += dt_y;
00423 m_tilt_x += dt_x;
00424 }
00425
00426
00427
00428
00429
00430
00431 void GeometryObject::rotation(const double* X, const double* Y, const unsigned size,
00432 const double C, const double S,
00433 double* Xrot, double* Yrot)
00434 {
00435 for(unsigned i=0; i<size; ++i) {
00436 Xrot[i] = X[i]*C - Y[i]*S;
00437 Yrot[i] = Y[i]*C + X[i]*S;
00438 }
00439 }
00440
00441
00442
00443 void GeometryObject::rotation(const double* X, const double* Y, const unsigned size, const double angle_deg,
00444 double* Xrot, double* Yrot)
00445 {
00446 const double angle_rad = angle_deg * DEG_TO_RAD;
00447 const double C = cos(angle_rad);
00448 const double S = sin(angle_rad);
00449 rotation(X, Y, size, C, S, Xrot, Yrot);
00450 }
00451
00452
00453
00454 }
00455
00456
00457
00458
00459