trf++/test/trfcyl
diff -N PropCylRK_t.cpp
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ PropCylRK_t.cpp 8 Aug 2011 19:05:50 -0000 1.1
@@ -0,0 +1,423 @@
+// PropCylRK_t.cpp
+
+// Test PropCylRK.
+
+#include "PropCylRK.h"
+#include "PropCyl.h"
+#include <iostream>
+#include <string>
+#include <cassert>
+#include <sstream>
+#include <cmath>
+#include "objstream/StdObjStream.hpp"
+#include "trfutil/trfstream.h"
+#include "trfutil/TRFMath.h"
+#include "trfbase/PropStat.h"
+#include "trfbase/TrackVector.h"
+#include "trfbase/VTrack.h"
+#include "trfbase/ETrack.h"
+#include "SurfCylinder.h"
+#include "PropCylFieldConst.h"
+
+using std::cout;
+using std::cerr;
+using std::endl;
+using std::string;
+using std::ostringstream;
+using std::istringstream;
+using std::abs;
+
+using namespace trf;
+
+//**********************************************************************
+
+int main( ) {
+
+ string ok_prefix = "PropCylRK (I): ";
+ string error_prefix = "PropCylRK test (E): ";
+
+ cout << ok_prefix
+ << "-------- Testing component PropCylRK. --------" << endl;
+
+ // Make sure assert is enabled.
+ bool assert_flag = false;
+ assert ( ( assert_flag = true, assert_flag ) );
+ if ( ! assert_flag ) {
+ cerr << "Assert is disabled" << endl;
+ return 1;
+ }
+
+ //********************************************************************
+
+ ObjTable::register_type<PropCylFieldConst>();
+ ObjTable::register_type<PropCylRK>();
+
+ cout << trf_format;
+ cout << ok_prefix << "Test constructor." << endl;
+ PropCylRK prop("constant");
+ cout << prop << endl;
+ PropCylRK propmc("mc");
+ cout << prop << endl;
+
+ // Construct equivalent PropCyl propagator.
+
+ PropCyl prop0(ObjDoublePtr(new ObjDouble(-2.0)));
+
+ //********************************************************************
+
+ // Here we propagate some tracks both forward and backward and then
+ // each back to the original track. We check that the returned
+ // track parameters match those of the original track.
+ cout << ok_prefix << "Check reversibility." << endl;
+ double r1[] = { 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0 };
+ double r2[] = { 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0 };
+ double phi1[] = { 0.0, 1.0, -1.0, 2.0, 2.0, -2.0, 3.0, 0.0 };
+ double alf1[] = { 0.0, 0.1, -0.1, 0.9, -0.9, 0.9, 0.9, 3.0 };
+ double crv[] = { 0.0, 0.0, 0.01, -0.1, 0.1, -0.05, 0.1, 0.0 };
+ double z[] = {100.0, 100.0, 100.0, 100.0,100.0, 100.0,115.0,100.0 };
+ double tlm[] = { 0.0, -0.1, 0.5, -1.0, 2.0, -1.2, 0.35, 0.0 };
+ double fbdf[] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
+ double bfdf[] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
+ double maxdiff = 1.e-7;
+ int ntrk = 8;
+ int i;
+ for ( i=0; i<ntrk; ++i ) {
+ cout << "********** Propagate track " << i << ". **********" << endl;
+ PropStat pstat;
+ SurfCylinder scy1(r1[i]);
+ SurfCylinder scy2(r2[i]);
+ TrackVector vec1;
+ vec1(0) = phi1[i]; // phi
+ vec1(1) = z[i]; // z
+ vec1(2) = alf1[i]; // alpha
+ vec1(3) = tlm[i]; // tlam
+ vec1(4) = crv[i]; // curv
+ VTrack trv1(SurfacePtr(scy1.new_pure_surface()),vec1);
+ cout << " starting: " << trv1 << endl;
+ //
+ // Find the direction that will propagate this track from r1 to r2.
+ //
+ Propagator::PropDir dir = Propagator::FORWARD;
+ Propagator::PropDir rdir = Propagator::BACKWARD;
+ if(r2[i] > r1[i] && abs(alf1[i]) > PI2 ||
+ r2[i] < r1[i] && abs(alf1[i]) <= PI2) {
+ dir = Propagator::BACKWARD;
+ rdir = Propagator::FORWARD;
+ }
+
+ //
+ // Propagate.
+ VTrack trv2f = trv1;
+ pstat = prop.vec_dir_prop(trv2f,scy2,dir);
+ cout << " forward: " << trv2f << endl;
+ cout << pstat << endl;
+ if(dir == Propagator::FORWARD)
+ assert( pstat.forward() );
+ if(dir == Propagator::BACKWARD)
+ assert( pstat.backward() );
+
+ //
+ // Propagate using PropCyl and check difference.
+ VTrack trv2f0 = trv1;
+ pstat = prop0.vec_dir_prop(trv2f0,scy2,dir);
+ cout << " forward: " << trv2f0 << endl;
+ cout << pstat << endl;
+ double diff0 =
+ scy2.vec_diff(trv2f.get_vector(),trv2f0.get_vector()).amax();
+ cout << "diff: " << diff0 << endl;
+ assert( diff0<maxdiff );
+
+ //
+ // Propagate in reverse direction.
+ VTrack trv2fb = trv2f;
+ pstat = prop.vec_dir_prop(trv2fb,scy1,rdir);
+ cout << " f return: " << trv2fb << endl;
+ cout << pstat << endl;
+ if(rdir == Propagator::FORWARD)
+ assert( pstat.forward() );
+ if(rdir == Propagator::BACKWARD)
+ assert( pstat.backward() );
+ // Check the return differences.
+ double difff =
+ scy1.vec_diff(trv2fb.get_vector(),trv1.get_vector()).amax();
+ cout << "diff: " << difff << endl;
+ assert( difff<maxdiff );
+ //
+ // Check no-move forward propagation to the same surface.
+ VTrack trv1s = trv1;
+ pstat = prop.vec_dir_prop(trv1s,scy1,Propagator::FORWARD);
+ cout << " same f: " << trv1s << endl;
+ cout << pstat << endl;
+ assert( pstat.same() );
+ assert( pstat.path_distance() == 0 );
+ //
+ // Check no-move backward propagation to the same surface.
+ trv1s = trv1;
+ pstat = prop.vec_dir_prop(trv1s,scy1,Propagator::BACKWARD);
+ cout << " same b: " << trv1s << endl;
+ cout << pstat << endl;
+ assert( pstat.same() );
+ assert( pstat.path_distance() == 0 );
+ //
+ // Check move propagation to the same surface.
+ //
+ // forward
+ int successes = 0;
+ trv1s = trv1;
+ pstat = prop.vec_dir_prop(trv1s,scy1,Propagator::FORWARD_MOVE);
+ cout << " forward move: " << trv1s << endl;
+ cout << pstat << endl;
+ if ( pstat.forward() ) ++successes;
+ // backward
+ trv1s = trv1;
+ pstat = prop.vec_dir_prop(trv1s,scy1,Propagator::BACKWARD_MOVE);
+ cout << " backward move: " << trv1s << endl;
+ cout << pstat << endl;
+ if ( pstat.backward() ) ++successes;
+ // Neither of these should have succeeded.
+ assert( successes == 0 );
+ //
+ // nearest
+ trv1s = trv1;
+ pstat = prop.vec_dir_prop(trv1s,scy1,Propagator::NEAREST_MOVE);
+ cout << " nearest move: " << trv1s << endl;
+ cout << pstat << endl;
+ assert ( !pstat.success() );
+
+ // Test derivatives numerically using uniform field.
+ cout << "Testing derivatives with uniform field" << endl;
+ VTrack trv1a = trv1;
+ TrackDerivative deriv;
+ pstat = prop.vec_dir_prop(trv1a,scy2,Propagator::NEAREST, &deriv);
+ assert(pstat.success());
+ VTrack trv1a0 = trv1;
+ TrackDerivative deriv0;
+ pstat = prop0.vec_dir_prop(trv1a0,scy2,Propagator::NEAREST, &deriv0);
+ assert(pstat.success());
+ for(int j=0; j<5; ++j) {
+ double d;
+ if(j < 2)
+ d = 0.1;
+ else if(j < 4)
+ d = 0.01;
+ else
+ d = 0.001;
+ TrackVector vec1b = vec1;
+ TrackVector vec1c = vec1;
+ vec1b(j) = vec1(j) + d;
+ vec1c(j) = vec1(j) - d;
+ VTrack trv1b(SurfacePtr(scy1.new_pure_surface()), vec1b);
+ VTrack trv1c(SurfacePtr(scy1.new_pure_surface()), vec1c);
+ pstat = prop.vec_dir_prop(trv1b,scy2,Propagator::NEAREST);
+ assert(pstat.success());
+ pstat = prop.vec_dir_prop(trv1c,scy2,Propagator::NEAREST);
+ assert(pstat.success());
+ for(int i=0; i<5; ++i) {
+ double dij = (trv1b.get_vector(i) - trv1c.get_vector(i))/(2.*d);
+ if(i == 0)
+ dij = fmod2(dij, TWOPI);
+ cout << i << ", " << j << '\t' << dij << '\t' << deriv(i,j)
+ << '\t' << deriv0(i,j) << endl;
+ double scale = max(abs(deriv(i,j)), 1.);
+ assert(abs(deriv(i,j)-deriv0(i,j))/scale < maxdiff);
+ assert(abs(dij-deriv(i,j))/scale < 1.e-4);
+ }
+ }
+
+ // Test derivatives numerically using non-uniform field.
+ cout << "Testing derivatives with non-uniform field" << endl;
+ VTrack trv1an = trv1;
+ TrackDerivative derivn;
+ pstat = propmc.vec_dir_prop(trv1an,scy2,Propagator::NEAREST, &derivn);
+ assert(pstat.success());
+ for(int j=0; j<5; ++j) {
+ double d;
+ if(j < 2)
+ d = 0.1;
+ else if(j < 4)
+ d = 0.01;
+ else
+ d = 0.001;
+ TrackVector vec1bn = vec1;
+ TrackVector vec1cn = vec1;
+ vec1bn(j) = vec1(j) + d;
+ vec1cn(j) = vec1(j) - d;
+ VTrack trv1bn(SurfacePtr(scy1.new_pure_surface()), vec1bn);
+ VTrack trv1cn(SurfacePtr(scy1.new_pure_surface()), vec1cn);
+ pstat = propmc.vec_dir_prop(trv1bn,scy2,Propagator::NEAREST);
+ assert(pstat.success());
+ pstat = propmc.vec_dir_prop(trv1cn,scy2,Propagator::NEAREST);
+ assert(pstat.success());
+ for(int i=0; i<5; ++i) {
+ double dijn = (trv1bn.get_vector(i) - trv1cn.get_vector(i))/(2.*d);
+ if(i == 0)
+ dijn = fmod2(dijn, TWOPI);
+ cout << i << ", " << j << '\t' << dijn << '\t' << derivn(i,j) << endl;
+ double scale = max(abs(derivn(i,j)), 1.);
+ assert(abs(dijn-derivn(i,j))/scale < 2.e-3);
+ }
+ }
+ }
+
+ //********************************************************************
+
+ // Repeat the above with errors.
+ cout << ok_prefix << "Check reversibility with errors." << endl;
+ double epp[] = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
+ double epz[] = { 0.01, -0.01, 0.01, -0.01, 0.01, -0.01, 0.01, 0.01 };
+ double ezz[] = { 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25 };
+ double epa[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004, 0.004 };
+ double eza[] = { 0.04, -0.04, 0.04, -0.04, 0.04, -0.04, 0.04, 0.04 };
+ double eaa[] = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
+ double epl[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004, 0.004 };
+ double eal[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004, 0.004 };
+ double ezl[] = { 0.04, -0.04, 0.04, -0.04, 0.04, -0.04, 0.04, 0.04 };
+ double ell[] = { 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02 };
+ double epc[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004, 0.004 };
+ double ezc[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004, 0.004 };
+ double eac[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004, 0.004 };
+ double elc[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004, 0.004 };
+ double ecc[] = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
+ for ( i=0; i<ntrk; ++i ) {
+ cout << "********** Propagate track " << i << ". **********" << endl;
+ PropStat pstat;
+ SurfCylinder scy1(r1[i]);
+ SurfCylinder scy2(r2[i]);
+ TrackVector vec1;
+ vec1(0) = phi1[i]; // phi
+ vec1(1) = z[i]; // z
+ vec1(2) = alf1[i]; // alpha
+ vec1(3) = tlm[i]; // tlam
+ vec1(4) = crv[i]; // curv
+ TrackError err1;
+ err1(0,0) = epp[i];
+ err1(0,1) = epz[i];
+ err1(1,1) = ezz[i];
+ err1(0,2) = epa[i];
+ err1(1,2) = eza[i];
+ err1(2,2) = eaa[i];
+ err1(0,3) = epl[i];
+ err1(1,3) = ezl[i];
+ err1(2,3) = eal[i];
+ err1(3,3) = ell[i];
+ err1(0,4) = epc[i];
+ err1(1,4) = ezc[i];
+ err1(2,4) = eac[i];
+ err1(3,4) = elc[i];
+ err1(4,4) = ecc[i];
+ ETrack tre1(SurfacePtr( scy1.new_pure_surface()), vec1, err1 );
+ cout << " starting: " << tre1 << endl;
+ //
+ // Find the direction that will propagate this track from r1 to r2.
+ //
+ Propagator::PropDir dir = Propagator::FORWARD;
+ Propagator::PropDir rdir = Propagator::BACKWARD;
+ if(r2[i] > r1[i] && abs(alf1[i]) > PI2 ||
+ r2[i] < r1[i] && abs(alf1[i]) <= PI2) {
+ dir = Propagator::BACKWARD;
+ rdir = Propagator::FORWARD;
+ }
+
+ //
+ // Propagate.
+ ETrack tre2f = tre1;
+ pstat = prop.err_dir_prop(tre2f,scy2,dir);
+ cout << " forward: " << tre2f << endl;
+ cout << pstat << endl;
+ if(dir == Propagator::FORWARD)
+ assert( pstat.forward() );
+ if(dir == Propagator::BACKWARD)
+ assert( pstat.backward() );
+
+ //
+ // Propagate using PropCyl and check difference.
+ ETrack tre2f0 = tre1;
+ pstat = prop0.err_dir_prop(tre2f0,scy2,dir);
+ cout << " forward: " << tre2f0 << endl;
+ cout << pstat << endl;
+ double vdiff0 =
+ scy2.vec_diff(tre2f.get_vector(),tre2f0.get_vector()).amax();
+ cout << "vec diff: " << vdiff0 << endl;
+ assert( vdiff0<maxdiff );
+ TrackError df0 = tre2f.get_error() - tre2f0.get_error();
+ double ediff0 = df0.amax();
+ cout << "err diff: " << ediff0 << endl;
+ assert( ediff0<maxdiff );
+
+ //
+ // Propagate backward
+ ETrack tre2fb = tre2f;
+ pstat = prop.err_dir_prop(tre2fb,scy1,rdir);
+ cout << " f return: " << tre2fb << endl;
+ cout << pstat << endl;
+ if(rdir == Propagator::FORWARD)
+ assert( pstat.forward() );
+ if(rdir == Propagator::BACKWARD)
+ assert( pstat.backward() );
+ double difff =
+ scy1.vec_diff(tre2fb.get_vector(),tre1.get_vector()).amax();
+ cout << "vec diff: " << difff << endl;
+ assert( difff<maxdiff );
+ TrackError dfb = tre2fb.get_error() - tre1.get_error();
+ double edifff = dfb.amax();
+ cout << "err diff: " << edifff << endl;
+ assert( edifff<maxdiff );
+ }
+
+ //********************************************************************
+
+ cout << ok_prefix << "Test cloning." << endl;
+ assert( prop.new_propagator() != 0 );
+
+ //********************************************************************
+
+ cout << ok_prefix << "Test the field." << endl;
+ assert( prop.get_bfield() != 0 );
+ assert( &*prop.get_bfield_map() != 0 );
+
+ //********************************************************************
+
+ cout << ok_prefix << "Write object data." << endl;
+ {
+ ObjPtr pobj( new PropCylRK("constant") );
+ ostringstream mystream;
+ StdObjStream objstream(mystream);
+ objstream.write_object("obj1",pobj);
+ cout << mystream.str() << endl;
+ assert( ObjTable::has_object_name("obj1") );
+ string::size_type pos = 0;
+ assert( (pos=mystream.str().find( "obj1 ",pos)) != string::npos );
+ assert( (pos=mystream.str().find( "PropCylRK ",pos)) != string::npos );
+ assert( (pos=mystream.str().find( "type",pos)) != string::npos );
+ assert( (pos=mystream.str().find( "precision",pos)) != string::npos );
+ assert( (pos=mystream.str().find( "scale1",pos)) != string::npos );
+ assert( (pos=mystream.str().find( "scale2",pos)) != string::npos );
+ }
+
+ //********************************************************************
+
+ cout << ok_prefix << "Read object data." << endl;
+ {
+ string istring = "[ obj2 PropCylRK type=\"constant\" precision=1.e-6 dir=\"forward\" ]";
+ istringstream isstrm(istring);
+ {
+ StdObjStream obstr(isstrm);
+ string name = obstr.read_object();
+ assert( name == "obj2" );
+ const PropCylRK* pobj;
+ ObjTable::get_object(name,pobj);
+ assert( pobj != 0 );
+ assert( pobj->get_type() == PropCylRK::get_static_type() );
+ assert( pobj->get_direction() == Propagator::FORWARD );
+ }
+ }
+
+ //********************************************************************
+
+ cout << ok_prefix
+ << "------------- All tests passed. -------------" << endl;
+ return 0;
+
+ //********************************************************************
+
+}