#ifdef TRIPLETS_TRIPLET_CU_H_
#error "this file can only be included once"
#endif

#define TRIPLETS_TRIPLET_CU_H_

#include "data.h"
#include "layout-ptr.h"
#include "triplet.h"
#include "util.h"
#include "vec.h"

template<class Tl, class T>
hostdevice__ void ParamTripletBase_flex<Tl, T>::set_pivot_hit
(const Hit &pivot_hit, layptr<TlTube, Tube> tubes) {
	this->pivot_time = this->max_hit_time = this->min_hit_time 
		= pivot_hit.t;
	//sum_hit_time = pivot_hit.t;
	//oldest_tube_id = pivot_hit.tube_id;
	this->pivot_tube_id = pivot_hit.tube_id;
	this->nhits = 1;
	//sum_hit_pos = tubes[pivot_hit.tube_id].pos;
	((T *)this)->update_cms(pivot_hit, tubes);
	//is_active = true;
}  // set_pivot_hit

template<class Tl>
hostdevice__ void Triplet_flex<Tl>::update_cms
(const Hit &hit, layptr<TlTube, Tube> tubes) {
	Tube tube = tubes[hit.tube_id];
	fvec3 tube_pos = tube.pos;
	this->sum_hit_pos += tube_pos;
}  // update_cms

template<class Tl>
hostdevice__ void HalfSkewlet_flex<Tl>::update_cms
(const Hit &hit, layptr<TlTube, Tube> tubes) {
	Tube tube = tubes[hit.tube_id];
	fvec3 tube_pos = tube.pos;
	float pz = tube_pos.z;
	if(pz > 36.0f || pz < 34.0f) {
		// skewed tube, special handling required
		// TODO: for performance optimization, make a specialization for 
		// non-skewed case
		float dz = tube.dir.z, k = (35.0f - pz) / dz;
		float half_length = tube.half_length;
		zmin = max(zmin, pz - dz * half_length);
		zmax = min(zmax, pz + dz * half_length);
		tube_pos += k * tube.dir;
	} 
	this->sum_hit_pos += tube_pos;
}  // update_cms

template<class Tl, class T>
hostdevice__ void ParamTripletBase_flex<Tl, T>::add_hit
(const Hit &hit, layptr<TlTube, Tube> tubes) {
	//sum_hit_time += hit.t;
	this->nhits++;
	//sum_hit_pos += tubes[hit.tube_id].pos;
	((T *)this)->update_cms(hit, tubes);
	if(hit.t < this->min_hit_time) {
		this->min_hit_time = hit.t;
		//oldest_tube_id = hit.tube_id;
	} else if(hit.t > this->max_hit_time)
		this->max_hit_time = hit.t;
}  // add_hit

template<class Tl, class T>
hostdevice__ void ParamTripletBase_flex<Tl, T>::add_neighbors
	(const int *nneighbors, const int *neighbor_tubes, int neighbor_pitch, 
	 layptr<TlHit, Hit> tube_hits, const int *tube_hit_starts, 
	 layptr<TlTube, Tube> tubes) {
	int itube = this->pivot_tube_id;
	for(int jneighbor = 0; jneighbor < nneighbors[itube]; jneighbor++) {
		int jtube = neighbor_tubes[jneighbor * neighbor_pitch + itube];
		int jhit_start = tube_hit_starts[jtube];
		int jhit_end = tube_hit_starts[jtube + 1];
		for(int jhit = jhit_start; jhit < jhit_end; jhit++) {
			if(this->accepts_hit(tube_hits[jhit]))
				add_hit(tube_hits[jhit], tubes);
		}
	}  // for(each neighbor)
}  // add_neighbors

template<class Tl>
hostdevice__ bool TripletBase_flex<Tl>::accepts_hit(const Hit &hit) const {
	return 
		pivot_time - DRIFT_TIME < hit.t && hit.t < pivot_time + DRIFT_TIME && 
		min_accept_window() <= hit.t && hit.t <= max_accept_window();
}  // accepts_hit

template<class Tl>
hostdevice__ void Skewlet_flex<Tl>::add_triplet(const HalfSkewlet &triplet) {
	//sum_hit_time += triplet.sum_hit_time;
	this->nhits += triplet.nhits;
	// TODO: ensure that CMS is updated correctly, and that 3D CMS update is not
	// required when combining half-skewlets
	this->sum_hit_pos += triplet.sum_hit_pos;
	if(triplet.min_hit_time < this->min_hit_time) {
		this->min_hit_time = triplet.min_hit_time;
		//oldest_tube_id = triplet.oldest_tube_id;
	}
	if(triplet.max_hit_time > this->max_hit_time)
		this->max_hit_time = triplet.max_hit_time;
}  // add_triplet

inline __host__ __device__ fvec3 find_crossing
(const fvec3 &p1, const fvec3 &u1, const fvec3 &p2, const fvec3 &u2) {
	fvec3 p21 = p2 - p1;
	if(len2(p21) < 1e-12f)
		return p1;
	fvec3 m = cross(u2, u1);
	float mlen2 = len2(m);
	if(mlen2 < 1e-12f)
		return p1 + 0.5f * p21;
	float inv_mlen2 = 1.0f / mlen2;
	//fvec3 mnorm = inv_mlen2 * m, r = cross(p21, mnorm);
	fvec3 r = cross(p21, m);
	float t1 = dot(r, u2), t2 = dot(r, u1);
	t1 *= inv_mlen2;
	t2 *= inv_mlen2;
	fvec3 q1 = p1 + t1 * u1, q2 = p2 + t2 * u2; 
	//fvec3 q21 = q2 - q1;
	//return q1 + 0.5f * q21;
	return 0.5f * (q1 + q2);
}  // find_crossing

inline __host__ __device__ fvec3 find_crossing_fast
(const fvec3 &p1, const fvec3 &u1, const fvec3 &p2, const fvec3 &u2) {
	fvec2 p21 = p2.xy() - p1.xy();
	fvec2 m = cross(u2, u1).xy();
	float mlen2 = len2(m);
	float inv_mlen2 = 1.0f / mlen2;
	float r = cross(p21, m) * inv_mlen2;
	fvec2 q21xy = 0.5f * (p1.xy() + p2.xy());
	float q21z = p1.z + r * u1.z * u2.z;
	return fvec3::from(q21xy.x, q21xy.y, q21z);
}  // find_crossing