From b5a40a1637a5c8c928a744d71cdf8664e549fa0b Mon Sep 17 00:00:00 2001 From: Alex Buts Date: Thu, 10 Oct 2013 12:24:25 +0100 Subject: [PATCH] refs #8095 presumably fixed issue but the unit tests still need some work on --- .../inc/MantidDataHandling/FindDetectorsPar.h | 90 +++-- .../DataHandling/src/FindDetectorsPar.cpp | 373 ++++++++++-------- 2 files changed, 273 insertions(+), 190 deletions(-) diff --git a/Code/Mantid/Framework/DataHandling/inc/MantidDataHandling/FindDetectorsPar.h b/Code/Mantid/Framework/DataHandling/inc/MantidDataHandling/FindDetectorsPar.h index 3f424b214a9e..b494964644c3 100644 --- a/Code/Mantid/Framework/DataHandling/inc/MantidDataHandling/FindDetectorsPar.h +++ b/Code/Mantid/Framework/DataHandling/inc/MantidDataHandling/FindDetectorsPar.h @@ -112,6 +112,9 @@ namespace Mantid File change history is stored at: . Code Documentation is available at: */ + // predefine class, used to cashe precalculated detector's parameters + class DetParameters; + class DLLExport FindDetectorsPar : public API::Algorithm { public: @@ -127,12 +130,12 @@ namespace Mantid /// the accessors, used to return algorithm results when called as Child Algorithm, without setting the properties; std::vectorconst & getAzimuthal()const{return azimuthal;} std::vectorconst & getPolar()const{return polar;} - std::vectorconst & getAzimWidth()const{return azimuthal_width;} - std::vectorconst & getPolarWidth()const{return polar_width;} - std::vectorconst & getFlightPath()const{return secondary_flightpath;} - std::vectorconst & getDetID()const{return det_ID;} + std::vectorconst & getAzimWidth()const{return azimuthalWidth;} + std::vectorconst & getPolarWidth()const{return polarWidth;} + std::vectorconst & getFlightPath()const{return secondaryFlightpath;} + std::vectorconst& getDetID()const{return detID;} /// number of real detectors, calculated by algorithm - size_t getNDetectors()const{return nDetectors;} + size_t getNDetectors()const{return m_nDetectors;} private: /// Sets documentation strings for this algorithm virtual void initDocs(); @@ -141,36 +144,29 @@ namespace Mantid void exec(); /** the variable defines if algorithm needs to calculate linear ranges for the detectors (dX,dY) * instead of azimuthal_width and polar_width */ - bool return_linear_ranges; - // number of real (not monitors)detectors, processed by algorithm; - size_t nDetectors; + bool m_SizesAreLinear; + + // numner of real(valid and non-monitor) detectors calculated by the alvorithm + size_t m_nDetectors; + // the vectors which represent detector's parameters as linear structure std::vector azimuthal; std::vector polar; - std::vector azimuthal_width; - std::vector polar_width; - std::vector secondary_flightpath; - std::vector width; - std::vector height; - std::vector det_ID; + std::vector azimuthalWidth; + std::vector polarWidth; + std::vector secondaryFlightpath; + std::vector detID; /// logger -> to provide logging, for MD workspaces static Kernel::Logger& g_log; - /// calculates par values for a detectors ring; - void calc_cylDetPar(const Geometry::IDetector_const_sptr spDet, - const Geometry::IObjComponent_const_sptr sample, - const Kernel::V3D &groupCentre, - double &azim, double &polar, double &azim_width, double &polar_width,double &dist); - /// calculates par values for a detectors block or a detector; - void calc_rectDetPar(const API::MatrixWorkspace_sptr inputWS, - const Geometry::IDetector_const_sptr spDet, - const Geometry::IObjComponent_const_sptr sample, - const Kernel::V3D &groupCentre, - double &azim, double &polar, double &azim_width, double &polar_width,double &dist); + // calculate generic detectors parameters: + void calcDetPar(const Geometry::IDetector_const_sptr &spDet,const Kernel::V3D &GroupCenter,DetParameters &Detector); /// if ASCII file is selected as the datasource, this structure describes the type of this file. FileTypeDescriptor current_ASCII_file; /// internal function which sets the output table according to the algorithms properties - void set_output_table(); + void setOutputTable(); + /// extract valid detectors parameters into vectors above + void extractAndLinearize(const std::vector &detPar); /// functions used to populate data from the phx or par file void populate_values_from_file(const API::MatrixWorkspace_sptr & inputWS); /// load data from par or phx file; @@ -187,6 +183,48 @@ namespace Mantid void load_plain(std::ifstream &stream,std::vector &Data,FileTypeDescriptor const &FILE_TYPE); }; +/**Small helper class used to precalculate the detectors parameters */ +class DetParameters +{ +public: + /// azimuthal detector's angle in spherical coordinate system alighned with the beam + double azimutAngle; + /// polar detector's angle in spherical coordinate system alighned with the beam + double polarAngle; + /// scattering source -- detector' distance + double secondaryFlightPath; + /// linear or angular size of the bounding box encapsulating detector and alighned tangentially to the constant scattering angle circle + double azimWidth,polarWidth; + /// the detector's ID + int64_t detID; + // default detector ID -- -1 means undefined + DetParameters():detID(-1){} + +}; + +/** helper class-collection to keep together the parameters, which characterize average composite detector + and help to calculate these parameters*/ +class AvrgDetector +{ + double m_AzimutSum; + double m_PolarSum; + double m_FlightPathSum; + double m_AzimMin,m_PolarMin,m_AzimMax,m_PolarMax; + // if azimuthal and polar sizes expressed in angular or linear units + bool m_useSphericalSizes; + /// numbr of primary detectors, contributing into this detector + size_t m_nComponents; +public: + AvrgDetector():m_AzimutSum(0),m_PolarSum(0),m_FlightPathSum(0), + m_useSphericalSizes(false), + m_AzimMin(FLT_MAX),m_PolarMin(FLT_MAX),m_AzimMax(-FLT_MAX),m_PolarMax(-FLT_MAX), + m_nComponents(0) + {} + void addDetInfo(const Geometry::IDetector_const_sptr &spDet,const Kernel::V3D &Observer); + void returnAvrgDetPar(DetParameters &det); + +}; + } //end namespace DataHandling } //end namespace Mandid diff --git a/Code/Mantid/Framework/DataHandling/src/FindDetectorsPar.cpp b/Code/Mantid/Framework/DataHandling/src/FindDetectorsPar.cpp index 453bca6200dd..998e576c4839 100644 --- a/Code/Mantid/Framework/DataHandling/src/FindDetectorsPar.cpp +++ b/Code/Mantid/Framework/DataHandling/src/FindDetectorsPar.cpp @@ -67,7 +67,7 @@ using namespace Kernel; using namespace API; // nothing here according to mantid FindDetectorsPar::FindDetectorsPar(): -return_linear_ranges(false) +m_SizesAreLinear(false) {}; FindDetectorsPar::~FindDetectorsPar(){}; @@ -108,7 +108,7 @@ FindDetectorsPar::exec() throw(Kernel::Exception::NotFoundError("can not obtain InoputWorkspace for the algorithm to work","")); } // Number of spectra - const size_t nHist = inputWS->getNumberHistograms(); + const int64_t nHist = (int64_t)inputWS->getNumberHistograms(); // try to load par file if one is availible std::string fileName = this->getProperty("ParFile"); @@ -118,9 +118,9 @@ FindDetectorsPar::exec() throw(Kernel::Exception::FileError(" file not exist",fileName)); } size_t nPars = loadParFile(fileName); - if(nPars == nHist){ + if(nPars ==(size_t)nHist){ this->populate_values_from_file(inputWS); - this->set_output_table(); + this->setOutputTable(); return; }else{ g_log.warning()<<" number of parameters in the file: "<getProperty("ReturnLinearRanges"); + m_SizesAreLinear = this->getProperty("ReturnLinearRanges"); - // Get a pointer to the sample - Geometry::IObjComponent_const_sptr sample =inputWS->getInstrument()->getSample(); - - azimuthal.assign(nHist,std::numeric_limits::quiet_NaN()); - polar.assign(nHist,std::numeric_limits::quiet_NaN()); - azimuthal_width.assign(nHist,std::numeric_limits::quiet_NaN()); - polar_width.assign(nHist,std::numeric_limits::quiet_NaN()); - secondary_flightpath.assign(nHist,std::numeric_limits::quiet_NaN()); - det_ID.assign(nHist,std::numeric_limits::quiet_NaN()); - this->nDetectors = 0; + std::vector Detectors(nHist); + DetParameters AverageDetector; + this->m_nDetectors = 0; Progress progress(this,0,1,100); const int progStep = (int)(ceil(double(nHist)/100.0)); + + // define the centre of coordinates: + Kernel::V3D Observer =inputWS->getInstrument()->getSample()->getPos(); + // Loop over the spectra - size_t ic(0); - for (size_t i = 0; i < nHist; i++) + //PRAGMA_OMP + for (int64_t i = 0; i < nHist; i++) { Geometry::IDetector_const_sptr spDet; try{ @@ -160,34 +157,30 @@ FindDetectorsPar::exec() // Check that we aren't writing a monitor... if (spDet->isMonitor())continue; - det_ID[ic] = spDet->getID(); - - Kernel::V3D groupCentre; - Geometry::det_topology group_shape= spDet->getTopology(groupCentre); - if(group_shape == Geometry::cyl){ // we have a ring; - calc_cylDetPar(spDet,sample,groupCentre,azimuthal[ic], polar[ic], - azimuthal_width[ic], polar_width[ic],secondary_flightpath[ic]); - }else{ // we have a detector or a rectangular shape - calc_rectDetPar(inputWS,spDet,sample,groupCentre,azimuthal[ic],polar[ic], - azimuthal_width[ic],polar_width[ic],secondary_flightpath[ic]); - } - ic++ ; + + // valid detector has valid detID + Detectors[i].detID = spDet->getID(); + + // calculate all parameters for current composite detector + calcDetPar(spDet,Observer,Detectors[i]); + // make regular progress reports and check for canceling the algorithm if ( i % progStep == 0 ){ progress.report(); } } - nDetectors = ic; - this->set_output_table(); + this->extractAndLinearize(Detectors); + // if necessary set up table workspace with detectors parameters. + this->setOutputTable(); } // Constant for converting Radians to Degrees const double rad2deg = 180.0 / M_PI; + // functions defines the ouptput table with parameters -void -FindDetectorsPar::set_output_table() +void FindDetectorsPar::setOutputTable() { std::string output = getProperty("OutputParTable"); if(output.empty())return; @@ -206,7 +199,7 @@ FindDetectorsPar::set_output_table() m_result->addColumn("double","twoTheta"); m_result->addColumn("double","azimuthal"); m_result->addColumn("double","secondary_flightpath"); - if(return_linear_ranges){ + if(m_SizesAreLinear){ m_result->addColumn("double","det_width"); m_result->addColumn("double","det_height"); }else{ @@ -215,171 +208,228 @@ FindDetectorsPar::set_output_table() } m_result->addColumn("long64","detID"); - for(size_t i=0;iappendRow(); - row << polar[i] << azimuthal[i] << secondary_flightpath[i] << polar_width[i] << azimuthal_width[i] << (int64_t)det_ID[i]; + row << polar[i] << azimuthal[i] << secondaryFlightpath[i] << polarWidth[i] << azimuthalWidth[i] <getPos(); + Kernel::V3D toDet = (detPos - Observer); + + double dist2Det,Polar,Azimut; + // identify the detector' position in the beam coordinate system: + toDet.getSpherical(dist2Det,Polar,Azimut); + m_FlightPathSum +=dist2Det; + m_PolarSum +=Polar; + m_AzimutSum +=Azimut; + + + // centre of the azimuthal ring (the ring detectors form around the beam) + Kernel::V3D ringCentre(0,0,toDet.Z()); + + // Get the bounding box + Geometry::BoundingBox bbox; + std::vector coord(3); + + Kernel::V3D er(0,1,0),e_th,ez(0,0,1); //ez along beamline, which is always oz; (can be amended) + if(dist2Det)er = toDet/dist2Det; // direction to the detector + Kernel::V3D e_tg = er.cross_prod(ez); // tangential to the ring and anticloakwise; + + coord[0]=er; + coord[1]=ez; + coord[2]=e_tg ; + bbox.setBoxAlignment(ringCentre,coord); + + spDet->getBoundingBox(bbox); + + // linear extensions of the bounding box orientied tangentially to the equal scattering angle circle + double azimMin = bbox.xMin(); + double azimMax = bbox.xMax(); + double polarMin = bbox.zMin(); // bounding box has been rotated according to coord above, so z is along e_tg + double polarMax = bbox.zMax(); + + + if (m_useSphericalSizes) + { + if (dist2Det==0) dist2Det = 1; + + // convert to angular units + double polarHalfSize = rad2deg*atan2(0.5*(polarMax-polarMin), dist2Det); + double azimHalfSize = rad2deg*atan2(0.5*(polarMax-azimMax), dist2Det); + + double polarMin = Polar -polarHalfSize; + double polarMax = Polar +polarHalfSize; + double azimMin = Azimut -azimHalfSize; + double azimMax = Azimut +azimHalfSize; + + } + if (m_AzimMin>azimMin)m_AzimMin=azimMin; + if (m_AzimMaxpolarMin)m_PolarMin=polarMin; + if (m_PolarMax coord(3); - - // get vector leading from the sample to the ring centre - Kernel::V3D Observer = sample->getPos(); - coord[1] = (GroupCenter-Observer); - double d0 = coord[1].norm(); - coord[1] /= d0; - // access contributed detectors; - Geometry::DetectorGroup_const_sptr pDetGroup = boost::dynamic_pointer_cast(spDet); - if(!pDetGroup){ + + // return undefined detector parameters if no average detector is defined; + if(m_nComponents==0)return; + + avrgDet.azimutAngle = m_AzimutSum/m_nComponents; + avrgDet.polarAngle = m_PolarSum/m_nComponents; + avrgDet.secondaryFlightPath = m_FlightPathSum/m_nComponents; + + avrgDet.azimWidth =(m_AzimMax-m_AzimMin); + avrgDet.polarWidth=(m_PolarMax-m_PolarMin); + +} + +void FindDetectorsPar::calcDetPar(const Geometry::IDetector_const_sptr &spDet,const Kernel::V3D &Observer, + DetParameters &Detector) +{ + + // get number of basic detectors within the composit detector + size_t nDetectors = spDet->nDets(); + // define summator + AvrgDetector detSum; + // do we want spherical or linear box sizes? + detSum.m_useSphericalSizes = !m_SizesAreLinear; + + if( nDetectors == 1) + { + detSum.addDetInfo(spDet,Observer); + } + else + { + // access contributing detectors; + Geometry::DetectorGroup_const_sptr spDetGroup = boost::dynamic_pointer_cast(spDet); + if(!spDetGroup){ g_log.error()<<"calc_cylDetPar: can not downcast IDetector_sptr to detector group for det->ID: "<getID()< pDets = pDetGroup->getDetectors(); - Geometry::BoundingBox bbox; - - // loop through all detectors in the group - for(size_t i=0;igetPos(); - coord[0] = center-GroupCenter; - double d1 = coord[0].norm(); - coord[0] /= d1; - coord[2] = coord[0].cross_prod(coord[1]); - - // obtain the bounding box, aligned accordingly to the coordinates; - bbox.nullify(); - bbox.setBoxAlignment(center,coord); - pDets[i]->getBoundingBox(bbox); - - double d_min = d1+bbox.xMin(); if(d_mind1_max)d1_max = d_max; - double d_azim = (bbox.zMax()-bbox.zMin())/d1; - azim_width+=d_azim; - - d1_sum +=d1; - dist_sum +=d1*d1+d0*d0; + auto detectors = spDetGroup->getDetectors(); + auto it = detectors.begin(); + auto it_end = detectors.end(); + for(;it!=it_end;it++) + { + detSum.addDetInfo(*it,Observer); } - double dNdet= double(pDets.size()); - dist = sqrt(dist_sum/dNdet); - if(return_linear_ranges){ - polar_width = d1_max-d1_min; // the width of the detector's ring; - azim_width = 2*M_PI*(d1_sum/dNdet); // the length of the detector's ring - }else{ - polar_width = (atan2(d1_max,d0)-atan2(d1_min,d0))*rad2deg; - polar = atan2(d1_sum/dNdet,d0)*rad2deg; - - azim_width *= rad2deg; - } + } + // calculate averages and return the detector parameters + detSum.returnAvrgDetPar(Detector); + } -void -FindDetectorsPar::calc_rectDetPar(const API::MatrixWorkspace_sptr inputWS, - const Geometry::IDetector_const_sptr spDet,const Geometry::IObjComponent_const_sptr sample, - const Kernel::V3D &GroupCentre, - double &azim, double &polar, double &azim_width, double &polar_width,double &dist) +void FindDetectorsPar::extractAndLinearize(const std::vector &detPar) { - // Get Sample->Detector distance - dist = spDet->getDistance(*sample); - polar = inputWS->detectorTwoTheta(spDet)*rad2deg; - azim = spDet->getPhi()*rad2deg; - // Now let's work out the detector widths on basis of bounding box tangential to the 2Theta=const ring; - Kernel::V3D beamDetVector(GroupCentre.X(),GroupCentre.Y(),0); // group centre minus the projection of this centre to the beamline - beamDetVector.normalize(); - std::vector coord(3); - coord[0] = beamDetVector; - coord[1] = Kernel::V3D(0,0,1); // along beamline, which is always oz; (can be amended) - coord[2] = coord[0].cross_prod(coord[1]); // tangential to the ring and anticloakwise; - + size_t nDetectors; + + // provisional number + nDetectors = detPar.size(); + + this->azimuthal.resize(nDetectors); + this->polar.resize(nDetectors); + this->azimuthalWidth.resize(nDetectors); + this->polarWidth.resize(nDetectors); + this->secondaryFlightpath.resize(nDetectors); + this->detID.resize(nDetectors); + + nDetectors = 0; + for(size_t i=0;i(detPar[i].detID); + nDetectors++; + } + // store caluclated value + m_nDetectors = nDetectors; - // Get the bounding box - Geometry::BoundingBox bbox; - bbox.setBoxAlignment(GroupCentre,coord); + // resize to actual detector's number + this->azimuthal.resize(nDetectors); + this->polar.resize(nDetectors); + this->azimuthalWidth.resize(nDetectors); + this->polarWidth.resize(nDetectors); + this->secondaryFlightpath.resize(nDetectors); + this->detID.resize(nDetectors); - spDet->getBoundingBox(bbox); - double xsize = bbox.xMax() - bbox.xMin(); - double ysize = bbox.zMax() - bbox.zMin(); // bounding box has been rotated according to coord above, so z is along coord[2] - if(return_linear_ranges){ - polar_width = xsize; // width - azim_width = ysize; // height - }else{ - polar_width = 2*rad2deg*atan2((xsize/2.0), dist); - azim_width = 2*rad2deg*atan2((ysize/2.0), dist); - } } - + // -size_t -FindDetectorsPar::loadParFile(const std::string &fileName){ +size_t FindDetectorsPar::loadParFile(const std::string &fileName){ // load ASCII par or phx file std::ifstream dataStream; std::vector result; this->current_ASCII_file = get_ASCII_header(fileName,dataStream); load_plain(dataStream,result,current_ASCII_file); - this->nDetectors = current_ASCII_file.nData_records; + m_nDetectors = current_ASCII_file.nData_records; dataStream.close(); // transfer par data into internal algorithm parameters; - azimuthal.resize(nDetectors); - polar.resize(nDetectors); - det_ID.resize(nDetectors); + azimuthal.resize(m_nDetectors); + polar.resize(m_nDetectors); + detID.resize(m_nDetectors); int Block_size,shift; - if(current_ASCII_file.Type==PAR_type){ + if(current_ASCII_file.Type==PAR_type) + { + m_SizesAreLinear = true; Block_size = 5; // this value coinside with the value defined in load_plain shift = 0; - width.resize(nDetectors); - height.resize(nDetectors); - secondary_flightpath.resize(nDetectors,std::numeric_limits::quiet_NaN()); + azimuthalWidth.resize(m_nDetectors); + polarWidth.resize(m_nDetectors); + secondaryFlightpath.resize(m_nDetectors,std::numeric_limits::quiet_NaN()); - for(size_t i=0;igetNumberHistograms(); - if(this->current_ASCII_file.Type == PAR_type){ + if(this->current_ASCII_file.Type == PAR_type) + { // in this case data in azimuthal width and polar width are in fact real sizes in meters; have to transform it in into angular values - azimuthal_width.resize(nHist); - polar_width.resize(nHist); for (size_t i = 0; i < nHist; i++){ - if((azimuthal[i]>-45&&azimuthal[i]<45)||(azimuthal[i]>135)||(azimuthal[i]<-135)){ - azimuthal_width[i]=atan2(height[i],secondary_flightpath[i])*rad2deg; - polar_width[i] =atan2(width[i],secondary_flightpath[i])*rad2deg; - }else{ - azimuthal_width[i]=atan2(width[i],secondary_flightpath[i])*rad2deg; - polar_width[i] =atan2(height[i],secondary_flightpath[i])*rad2deg; - } + azimuthalWidth[i]=atan2(azimuthalWidth[i],secondaryFlightpath[i])*rad2deg; + polarWidth[i] =atan2(polarWidth[i],secondaryFlightpath[i])*rad2deg; } - height.resize(0); - width.resize(0); - }else{ + m_SizesAreLinear = false; + } + else + { Geometry::IObjComponent_const_sptr sample =inputWS->getInstrument()->getSample(); - secondary_flightpath.resize(nHist); + secondaryFlightpath.resize(nHist); // Loop over the spectra for (size_t i = 0; i < nHist; i++){ Geometry::IDetector_const_sptr spDet; @@ -417,7 +462,7 @@ FindDetectorsPar::populate_values_from_file(const API::MatrixWorkspace_sptr & in // Check that we aren't writing a monitor... if (spDet->isMonitor())continue; /// this is the only value, which is not defined in phx file, so we calculate it - secondary_flightpath[i] = spDet->getDistance(*sample); + secondaryFlightpath[i] = spDet->getDistance(*sample); } }