ionhit.cpp

Go to the documentation of this file.

/*
 * ionhit.cpp - Ion event data class
 * Copyright (C) 2014 D Haley
 * 
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#include "atomprobe/primitives/ionHit.h"
#include "atomprobe/primitives/boundcube.h"

#include "atomprobe/helper/aptAssert.h"
#include "atomprobe/helper/endianTest.h"

#ifdef _OPENMP
#include <omp.h>
#endif
using std::vector;
using AtomProbe::BoundCube;

#include <fstream>
#include <limits>
#include <utility>

namespace AtomProbe
{

IonHit::IonHit() 
{
    //At this point I deliberately don't initialise the point class
    //as in DEBUG mode, the point class will catch failure to init
}

IonHit::IonHit(float *buffer) : massToCharge(buffer[3]), pos(buffer)
{
}

IonHit::IonHit(float x, float y, float z, float m) : massToCharge(m), pos(x,y,z)
{
}

IonHit::IonHit(const IonHit &obj2) : massToCharge(obj2.massToCharge), pos(obj2.pos)
{
}

IonHit::IonHit(const Point3D &p, float newMass) : massToCharge(newMass), pos(p)
{
}

void IonHit::setMassToCharge(float newMass)
{
    massToCharge=newMass;
}

float IonHit::getMassToCharge() const
{
    return massToCharge;
}   

void IonHit::setPos(unsigned int idx, float value)
{
    ASSERT(idx < 3);
    pos[idx]=value;
}

void IonHit::setPos(const Point3D &p)
{
    pos=p;
}

#ifdef __LITTLE_ENDIAN__
void IonHit::switchEndian()
{
    
    pos.switchEndian();
    floatSwapBytes(&(massToCharge));
}
#endif

const IonHit &IonHit::operator=(const IonHit &obj) 
{
    massToCharge=obj.massToCharge;
    pos = obj.pos;

    return *this;
}

bool IonHit::operator==(const IonHit &obj)  const
{
    return massToCharge==obj.massToCharge && pos == obj.pos;

}
IonHit IonHit::operator+(const Point3D &obj)
{
    pos+=obj;   
    return *this;
}

void IonHit::getPoints(const vector<IonHit> &ions, vector<Point3D> &p)
{
    p.resize(ions.size());
#pragma omp parallel for
    for(size_t ui=0;ui<ions.size();ui++)
        p[ui] = ions[ui].getPosRef();
}

Point3D IonHit::getPos() const
{
    return pos;
}   

bool IonHit::hasNaN()
{
    return (std::isnan(massToCharge) || std::isnan(pos[0]) || 
                std::isnan(pos[1]) || std::isnan(pos[2]));
}

void IonHit::getCentroid(const std::vector<IonHit> &points,Point3D &centroid)
{
    centroid=Point3D(0,0,0);
    size_t nPoints=points.size();
#ifdef _OPENMP
    
    //Parallel version
    //--
    vector<Point3D> centroids(omp_get_max_threads(),Point3D(0,0,0));
#pragma omp parallel for 
    for(size_t ui=0;ui<nPoints;ui++)
        centroids[omp_get_thread_num()]+=points[ui].getPos();

    for(size_t ui=0;ui<centroids.size();ui++)
        centroid+=centroids[ui];
    //--

#else
    for(unsigned int ui=0;ui<nPoints;ui++)
        centroid+=points[ui].getPos();
#endif
    
    centroid*=1.0f/(float)nPoints;
}

void IonHit::getBoundCube(const std::vector<IonHit> &points,BoundCube &b)
{
    ASSERT(points.size());

    float bounds[3][2];
    for(unsigned int ui=0;ui<3;ui++)
    {
        bounds[ui][0]=std::numeric_limits<float>::max();
        bounds[ui][1]=-std::numeric_limits<float>::max();
    }
    
    for(unsigned int ui=0; ui<points.size(); ui++)
    {
        Point3D p;
        p=points[ui].getPos();
        for(unsigned int uj=0; uj<3; uj++)
        {
            bounds[uj][0] = std::min(p.getValue(uj),bounds[uj][0]);
            bounds[uj][1] = std::max(p.getValue(uj),bounds[uj][1]);
        }
    }

    b.setBounds(bounds[0][0],bounds[1][0],
            bounds[2][0],bounds[0][1],
            bounds[1][1],bounds[2][1]);

}

std::ostream& operator<<(std::ostream &stream, const IonHit &ion)
{
       stream << "(" << ion.getPos()[0] << "," << ion.getPos()[1] << "," << ion.getPos()[2] << "," << ion.getMassToCharge() << ")";
       return stream;
}


BoundCube IonHit::getIonDataLimits(const std::vector<IonHit> &points)
{
    ASSERT(points.size());

    BoundCube b;    
    b.setInverseLimits();   
#ifndef OPENMP
    float bounds[3][2];
    for(unsigned int ui=0;ui<3;ui++)
    {
        bounds[ui][0]=std::numeric_limits<float>::max();
        bounds[ui][1]=-std::numeric_limits<float>::max();
    }
    
    for(unsigned int ui=0; ui<points.size(); ui++)
    {

        for(unsigned int uj=0; uj<3; uj++)
        {
            Point3D p;
            p=points[ui].getPos();
            if(p.getValue(uj) < bounds[uj][0])
                bounds[uj][0] = p.getValue(uj);
            
            if(p.getValue(uj) > bounds[uj][1])
                bounds[uj][1] = p.getValue(uj);
        }
    }

    b.setBounds(bounds[0][0],bounds[1][0],
            bounds[2][0],bounds[0][1],
            bounds[1][1],bounds[2][1]);
#else
    // parallel version
    vector<BoundCube> cubes;

    unsigned int nT=omp_get_max_threads();
    cubes.resize(nT);
    for(unsigned int ui=0;ui<cubes.size();ui++)
        cube[ui].setInverseLimits();

    unsigned int tCount=1;
    #pragma omp parallel for reduction(tCount|+)
    for(unsigned int ui=0;ui<points.size();ui++)
    {
        Point3D p;
        p=points[ui].getPos();
        for(unsigned int uj=0;uj<3;uj++)
        {
            b.setBounds(uj,0,std::min(b.getBound(uj,0),p[uj]));
            b.setBounds(uj,1,std::min(b.getBound(uj,0),p[uj]));
        }
    }

    for(unsigned int ui=0;ui<std::min(tCount,nT);ui++)
        b.expand(cubes[ui]);

#endif

    return b;
}


float IonHit::operator[](unsigned int ui) const
{
    ASSERT(ui <= 3);
    if (ui < 3)
        return pos[ui];
    else
        return massToCharge;
}

float &IonHit::operator[](unsigned int ui)
{
    ASSERT(ui <=3);
    if (ui < 3) 
        return pos[ui];
    else 
        return massToCharge;
}



bool EPOS_ENTRY::operator==(const EPOS_ENTRY &obj)  const
{
    if(massToCharge != obj.massToCharge)
        return false;
    if(x!=obj.x || y != obj.y || z !=obj.z)
        return false;
    if(timeOfFlight!= obj.timeOfFlight)
        return false;
    if(voltDC != obj.voltDC)
        return false;
    if(voltPulse!= obj.voltPulse)
        return false;
    if(xDetector!= obj.xDetector || yDetector!=obj.yDetector)
        return false;
    if(deltaPulse!= obj.deltaPulse)
        return false;
    if(hitMultiplicity!= obj.hitMultiplicity)
        return false;

    return true;

}



#ifdef DEBUG

bool testIonHit()
{
    //test the boundcube function
    vector<IonHit> h;
    IonHit hit;
    hit.setMassToCharge(1);
    
    for(size_t ui=0;ui<8;ui++)
    {
        hit.setPos(Point3D(ui&4 >> 2,ui&2 >> 1,ui&1));
        h.push_back(hit);
    }

    BoundCube bc;
    IonHit::getBoundCube(h,bc);
    TEST(bc.isValid(),"check boundcube");


    return true;
}



#endif
};