Finding/documenting common geometric calculations

[drbenmorgan]

This is mostly a question for @yramachers, but feel free to comment if you have anything additional!

In working through #149, I found large and relatively kludgy blocks of vector/plane geometry calculations, e.g.

Falaise/source/falaise/snemo/simulation/gg_step_hit_processor.cc

Lines 478 to 739 in 3fa8064

	/* Scheme of the current step hit (XY view in DCCF)
	* and the generation of randomly distributed ion/pair
	* creation vertex along the hit segment:
	*
	* step hit segment
	*
	* M (minimum approach)
	* start stop
	* --------+==:=:==:=====*====:==:===>+----> track direction
	* t1 t : / t t2
	* : /
	* : / d: shortest drift distance
	* :/
	* +
	* C (the cell center == origin in DCCF)
	*/
	// loop on ionizations:
	for (size_t ion = 0; ion < number_of_ionizations; ++ion) {
	// shoot time and vertex:
	const double ran = get_rng().uniform();
	const double ran_ionization_time = hit_time_start + ran * (hit_time_stop - hit_time_start);
	const geomtools::vector_3d ran_ionization_pos = cell_hit_pos_start + ran * step_dir;
	const double hit_x = ran_ionization_pos.x();
	const double hit_y = ran_ionization_pos.y();
	const double hit_z = ran_ionization_pos.z();
	// compute the drift distance to the anodic wire:
	const double drift_distance = hypot(hit_x, hit_y);
	// compute the impact of the Geiger avalanche on the anodic wire (DCCF):
	geomtools::vector_3d avalanche_impact_cell_pos(0.0, 0.0, hit_z);
	// Check if the ionization occurs in a fiducial cylinder (if requested):
	// In the drift cell XY plane:
	if (datatools::is_valid(fiducialCellRadius_)) {
	// not in the fiducial drift X-Y circular region:
	if (drift_distance > fiducialCellRadius_) {
	// drop this ion/electron pair which will not
	// produce a Geiger avalanche and thus is sterile:
	continue;
	}
	}
	// Along the drift cell Z axis:
	if (datatools::is_valid(fiducialCellLength_)) {
	// not in the fiducial drift Z longitudinal region:
	// drop this ion/electron pair which will not
	// produce a Geiger avalanche and thus is sterile.
	if (std::abs(hit_z) > 0.5 * fiducialCellLength_) {
	continue;
	}
	}
	if (!datatools::is_valid(min_drift_distance)) {
	// The first processed ion/electron pair is accepted as
	// the closest to the anodic wire (so far):
	min_drift_distance = drift_distance;
	// record the time of the pair creation:
	ionization_time_discrete = ran_ionization_time;
	// record its position:
	ionization_world_pos_discrete =
	world_hit_pos_start + ran * (world_hit_pos_stop - world_hit_pos_start);
	// 2011-12-08 FM: added
	ionization_world_momentum_discrete = world_hit_momentum_start;
	// compute the impact of the Geiger avalanche
	// on the anodic wire (WCS):
	cell_world_plcmt.child_to_mother(avalanche_impact_cell_pos,
	avalanche_impact_world_pos_discrete);
	} else {
	if (drift_distance < min_drift_distance) {
	// This ion/electron pair is closest than the
	// current best candidate:
	min_drift_distance = drift_distance;
	// So we update now the ionization time and position
	// with this new pair:
	ionization_time_discrete = ran_ionization_time;
	ionization_world_pos_discrete =
	world_hit_pos_start + ran * (world_hit_pos_stop - world_hit_pos_start);
	// 2011-12-08 FM: added
	ionization_world_momentum_discrete = world_hit_momentum_start;
	// compute the impact of the Geiger avalanche
	// on the anodic wire (WCS):
	cell_world_plcmt.child_to_mother(avalanche_impact_cell_pos,
	avalanche_impact_world_pos_discrete);
	}
	}
	} // for
	// at this stage, the space/time coordinates of the best ion/electron
	// pair is available through:
	// 'ionization_time'
	// 'ionization_world_pos'
	// 'avalanche_impact_world_pos'
	// unless no first ionization ion/electron pair was produced for the
	// current step hit.
	} // end of 'if (algo_discrete)'
	if (algo_continuous) {
	/* Algorithm (b)
	*
	* For highly ionizing particle (alpha), ionization is considered as
	* a continuous process with a large amount of electron/ion pairs
	* created along the step:
	*/
	// initial boundaries of the step hit (DCCF):
	geomtools::vector_2d s1(cell_hit_pos_start.x(), cell_hit_pos_start.y());
	geomtools::vector_2d s2(cell_hit_pos_stop.x(), cell_hit_pos_stop.y());
	/* Scheme of a typical step hit (XY view in DCCF):
	*
	* step hit segment
	* S1 S2
	* --------+===============>+--------> track direction
	* t1 ~~~~~~~~~~~~~~~ t2
	* ~~~~~~~~~~~~
	* ~~~~~~~~~ <-- the Geiger avalanche
	* ~~~~~ region
	* ~~
	* +
	* C (the cell center == origin in DCCF)
	*/
	// If fiducial drift Z is defined:
	if (datatools::is_valid(fiducialCellLength_)) {
	// approximation (could be refined but at the price of tricky geometrics):
	const double mean_hit_z = 0.5 * (cell_hit_pos_start.z() + cell_hit_pos_stop.z());
	// not in the fiducial drift Z longitudinal region:
	if (std::abs(mean_hit_z) > 0.5 * fiducialCellLength_) {
	// consider this step as sterile
	// => no ionization is generated
	continue; // skip to next step hit
	}
	}
	// If fiducial drift radius is defined (XY plane),
	if (datatools::is_valid(fiducialCellRadius_)) {
	/* We recompute the effective boundaries of the step hit.
	* We must find the active part/segment of the
	* step hit segment that belongs to the circular fiducial
	* drift region.
	*/
	// compute start=S1 and stop=S2 of the step inside
	// the fiducial drift cylindric volume:
	geomtools::vector_2d effective_s1, effective_s2;
	geomtools::vector_2d cell_center(0., 0.);
	if (geomtools::intersection::find_intersection_segment_disk_2d(
	s1, s2, cell_center, fiducialCellRadius_, effective_s1, effective_s2)) {
	// update the fiducial part of the step hit segment:
	// this remove the part(s) of the step hit that lies
	// outside the drift region:
	s1 = effective_s1;
	s2 = effective_s2;
	} else {
	// step segment (s1,s2) is totally outside the fiducial
	// drift volume (XY plane)
	// => no ionization is generated
	continue; // skip to next step hit
	}
	}
	/* compute the coordinates of the point H along the (S1,S2) axis
	* that is the closest to the anodic wire (cell center C in DCCF):
	*/
	const geomtools::vector_2d s1s2 = s2 - s1;
	const geomtools::vector_2d u = s1s2.unit();
	double a1, b1, c1;
	double a2, b2, c2;
	double xh, yh;
	a1 = u.x();
	b1 = u.y();
	c1 = 0.0;
	a2 = u.y();
	b2 = -u.x();
	c2 = s1.x() * u.y() - s1.y() * u.x();
	if (mygsl::linear_system_2x2_solve(a1, b1, c1, a2, b2, c2, xh, yh)) {
	const geomtools::vector_2d h(xh, yh);
	const geomtools::vector_2d s1h = h - s1;
	if (((s1h.dot(u) > 0.0) && (s1h.mag() > s1s2.mag())) || (s1h.dot(u) < 0.0)) {
	// H does not belong to the step hit segment (S1,S2):
	const double d1 = s1.mag();
	const double d2 = s2.mag();
	if (d1 < d2) {
	/* Step start is the nearest:
	*
	* H S1 S2
	* ----+----+================+------------>
	* : / t1 t2
	* : /
	* : / d1
	* :/
	* C +
	*/
	ionization_time_continuous = hit_time_start;
	ionization_world_pos_continuous = world_hit_pos_start;
	// 2011-12-08 FM: added
	ionization_world_momentum_continuous = world_hit_momentum_start;
	const geomtools::vector_3d avalanche_impact_cell_pos(0.0, 0.0, cell_hit_pos_start.z());
	cell_world_plcmt.child_to_mother(avalanche_impact_cell_pos,
	avalanche_impact_world_pos_continuous);
	} else {
	/* Step stop is the nearest:
	*
	* S1 S2 H
	* --------+================+----+---------->
	* t1 t2 \ :
	* \ :
	* d2 \ :
	* \:
	* + C
	*/
	ionization_time_continuous = hit_time_stop;
	ionization_world_pos_continuous = world_hit_pos_stop;
	// 2011-12-08 FM: added
	ionization_world_momentum_continuous = world_hit_momentum_start;
	const geomtools::vector_3d avalanche_impact_cell_pos(0.0, 0.0, cell_hit_pos_stop.z());
	cell_world_plcmt.child_to_mother(avalanche_impact_cell_pos,
	avalanche_impact_world_pos_continuous);
	}
	} else {
	/* H is the nearest:
	*
	* step prop
	* <---------->
	* S1 H S2
	* --------+==========+======+-------->
	* t1 : t t2
	* :
	* d :
	* :
	* + C
	*/
	const double step_prop = s1h.mag() / s1s2.mag();
	// compute the time of nearest ionization:
	ionization_time_continuous =
	hit_time_start + step_prop * (hit_time_stop - hit_time_start);
	// compute the position of nearest ionization in (WCF):
	ionization_world_pos_continuous =
	world_hit_pos_start + step_prop * (world_hit_pos_stop - world_hit_pos_start);
	// 2011-12-08 FM: added
	ionization_world_momentum_continuous = world_hit_momentum_start;
	const double zh =
	cell_hit_pos_start.z() + step_prop * (cell_hit_pos_stop.z() - cell_hit_pos_start.z());
	const geomtools::vector_3d avalanche_impact_cell_pos(0.0, 0.0, zh);
	cell_world_plcmt.child_to_mother(avalanche_impact_cell_pos,
	avalanche_impact_world_pos_continuous);
	}
	}
	if (computeMinimumApproachPosition_) {
	minimum_approach_time = ionization_time_continuous;
	minimum_approach_world_pos = ionization_world_pos_continuous;
	minimum_approach_world_distance =
	(avalanche_impact_world_pos_continuous - ionization_world_pos_continuous).mag();
	}
	} // end of if (algo_continuous)

It would be good to document any existing functionality we have for such calculations (distance to plane, point of closest approach etc), or add any code needed to implement anything missing. @yramachers, I know you've done quite a bit here in reconstruction, so could you point me to the code if it's on GitHub please?