trf++/test/trfzp
diff -N PropZZRK_t.cpp
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ PropZZRK_t.cpp 8 Aug 2011 19:08:37 -0000 1.1
@@ -0,0 +1,440 @@
+// PropZZRK_t.cpp
+
+// Test PropZZRK.
+
+#include "PropZZRK.h"
+#include "PropZZ.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 "trfbase/TrackSurfaceDirection.h"
+#include "SurfZPlane.h"
+
+using std::cout;
+using std::cerr;
+using std::endl;
+using std::string;
+using std::ostringstream;
+using std::istringstream;
+using std::abs;
+using std::max;
+
+using namespace trf;
+
+//**********************************************************************
+
+int main( ) {
+
+ string ok_prefix = "PropZZRK (I): ";
+ string error_prefix = "PropZZRK test (E): ";
+
+ cout << ok_prefix
+ << "-------- Testing component PropZZRK. --------" << 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<PropZZRK>();
+
+ cout << trf_format;
+ cout << ok_prefix << "Test constructor." << endl;
+ PropZZRK prop("constant");
+ cout << prop << endl;
+ // PropZZRK propmc("mc");
+ // cout << prop << endl;
+
+ // Construct equivalent PropZZ propagator.
+
+ PropZZ 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 z1[] = {100.0,-100.0, 100.0,-100.0,100.0,-100.0,100.0,100.0 };
+ double z2[] = {150.0, -50.0, 50.0,-150.0, 50.0, -50.0,150.0, 50.0 };
+ double x[] = { 10.0, 1.0, -1.0, 2.0, 2.0, -2.0, 3.0, 0.0 };
+ double xv[] = { 0.5, 0.03, -0.03, 0.5, -0.5, 1.0, 1.0, 2.0 };
+ double crv[] = { 0.1, -0.1, 0.1, 0.01, 0.01, -1.0, 1.0, 0.02 };
+ double y[] = { 20.0, 0.0, 0.0, 0.0, 0.0, 0.0, 15.0, 0.0 };
+ double yv[] = { -0.5, 0.01,-0.02, 0.0, 0.0, 0.5, -0.5, 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 };
+ bool forw[] = { 1, 0, 0, 1, 0, 1, 1, 1 };
+ double maxdiff = 1.e-7;
+ int ntrk = 8;
+ int i;
+ for ( i=0; i<ntrk; ++i ) {
+ cout << "********** Propagate track " << i << ". **********" << endl;
+ PropStat pstat;
+ SurfZPlane sz1(z1[i]);
+ SurfZPlane sz2(z2[i]);
+ TrackVector vec1;
+ vec1(0) = x[i]; // x
+ vec1(1) = y[i]; // y
+ vec1(2) = xv[i]; // dx/dz
+ vec1(3) = yv[i]; // dy/dz
+ vec1(4) = crv[i]; // q/p
+ TrackSurfaceDirection tdir;
+ if(forw[i])
+ tdir = TSD_FORWARD;
+ else
+ tdir = TSD_BACKWARD;
+ VTrack trv1(SurfacePtr(sz1.new_pure_surface()), vec1, tdir);
+ cout << " starting: " << trv1 << endl;
+ //
+ // Find the direction that will propagate this track from z1 to z2.
+ //
+ Propagator::PropDir dir = Propagator::FORWARD;
+ Propagator::PropDir rdir = Propagator::BACKWARD;
+ if(z2[i] > z1[i] && tdir == TSD_BACKWARD ||
+ z2[i] < z1[i] && tdir == TSD_FORWARD) {
+ dir = Propagator::BACKWARD;
+ rdir = Propagator::FORWARD;
+ }
+
+ //
+ // Propagate.
+ VTrack trv2f = trv1;
+ pstat = prop.vec_dir_prop(trv2f,sz2,dir);
+ assert(pstat.success());
+ cout << " forward: " << trv2f << endl;
+ cout << pstat << endl;
+ if(dir == Propagator::FORWARD)
+ assert( pstat.forward() );
+ if(dir == Propagator::BACKWARD)
+ assert( pstat.backward() );
+
+ //
+ // Propagate using PropZZ and check difference.
+ VTrack trv2f0 = trv1;
+ pstat = prop0.vec_dir_prop(trv2f0,sz2,dir);
+ assert(pstat.success());
+ cout << " forward: " << trv2f0 << endl;
+ cout << pstat << endl;
+ double diff0 =
+ sz2.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,sz1,rdir);
+ assert(pstat.success());
+ 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 =
+ sz1.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,sz1,Propagator::FORWARD);
+ assert(pstat.success());
+ 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,sz1,Propagator::BACKWARD);
+ assert(pstat.success());
+ 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,sz1,Propagator::FORWARD_MOVE);
+ cout << " forward move: " << trv1s << endl;
+ cout << pstat << endl;
+ if ( pstat.forward() ) ++successes;
+ // backward
+ trv1s = trv1;
+ pstat = prop.vec_dir_prop(trv1s,sz1,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,sz1,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,sz2,Propagator::NEAREST, &deriv);
+ assert(pstat.success());
+ VTrack trv1a0 = trv1;
+ TrackDerivative deriv0;
+ pstat = prop0.vec_dir_prop(trv1a0,sz2,Propagator::NEAREST, &deriv0);
+ assert(pstat.success());
+ for(int j=0; j<5; ++j) {
+ double d;
+ if(j < 2)
+ d = 0.25;
+ 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(sz1.new_pure_surface()), vec1b, tdir);
+ VTrack trv1c(SurfacePtr(sz1.new_pure_surface()), vec1c, tdir);
+ pstat = prop.vec_dir_prop(trv1b,sz2,Propagator::NEAREST);
+ assert(pstat.success());
+ pstat = prop.vec_dir_prop(trv1c,sz2,Propagator::NEAREST);
+ assert(pstat.success());
+ for(int i=0; i<5; ++i) {
+ double dij = (trv1b.get_vector(i) - trv1c.get_vector(i))/(2.*d);
+ 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-3);
+ }
+ }
+
+ /*
+ // 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,sz2,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(sz1.new_pure_surface()), vec1bn, tdir);
+ VTrack trv1cn(SurfacePtr(sz1.new_pure_surface()), vec1cn, tdir);
+ pstat = propmc.vec_dir_prop(trv1bn,sz2,Propagator::NEAREST);
+ assert(pstat.success());
+ pstat = propmc.vec_dir_prop(trv1cn,sz2,Propagator::NEAREST);
+ assert(pstat.success());
+ for(int i=0; i<5; ++i) {
+ double dijn = (trv1bn.get_vector(i) - trv1cn.get_vector(i))/(2.*d);
+ 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 < 0.01);
+ }
+ }
+ */
+ }
+
+ //********************************************************************
+
+ // Repeat the above with errors.
+ cout << ok_prefix << "Check reversibility with errors." << endl;
+ double exx[] = { 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1 };
+ double exy[] = { 0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05, 0.05 };
+ double eyy[] = { 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25 };
+ double exxv[]= { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004, 0.004 };
+ double eyxv[]= { 0.02, -0.02, 0.02, -0.02, 0.02, -0.02, 0.02, 0.02 };
+ double exvxv[]={ 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
+ double exyv[]= { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004, 0.004 };
+ double eyyv[]= { 0.04, -0.04, 0.04, -0.04, 0.04, -0.04, 0.04, 0.04 };
+ double exvyv[]={ 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004, 0.004 };
+ double eyvyv[]={ 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02 };
+ double exc[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004, 0.004 };
+ double eyc[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004, 0.004 };
+ double exvc[]= { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004, 0.004 };
+ double eyvc[]= { 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;
+ SurfZPlane dz1(z1[i]);
+ SurfZPlane dz2(z2[i]);
+ TrackVector vec1;
+ vec1(0) = x[i]; // x
+ vec1(1) = y[i]; // y
+ vec1(2) = xv[i]; // dx/dz
+ vec1(3) = yv[i]; // dy/dz
+ vec1(4) = crv[i]; // curvature
+ TrackError err1;
+ err1(0,0) = exx[i];
+ err1(0,1) = exy[i];
+ err1(1,1) = eyy[i];
+ err1(0,2) = exxv[i];
+ err1(1,2) = eyxv[i];
+ err1(2,2) = exvxv[i];
+ err1(0,3) = exyv[i];
+ err1(1,3) = eyyv[i];
+ err1(2,3) = exvyv[i];
+ err1(3,3) = eyvyv[i];
+ err1(0,4) = exc[i];
+ err1(1,4) = eyc[i];
+ err1(2,4) = exvc[i];
+ err1(3,4) = eyvc[i];
+ err1(4,4) = ecc[i];
+ TrackSurfaceDirection tdir;
+ if(forw[i])
+ tdir = TSD_FORWARD;
+ else
+ tdir = TSD_BACKWARD;
+ ETrack tre1(SurfacePtr( dz1.new_pure_surface()), vec1, err1, tdir );
+ 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(z2[i] > z1[i] && tdir == TSD_BACKWARD ||
+ z2[i] < z1[i] && tdir == TSD_FORWARD) {
+ dir = Propagator::BACKWARD;
+ rdir = Propagator::FORWARD;
+ }
+
+ //
+ // Propagate.
+ ETrack tre2f = tre1;
+ pstat = prop.err_dir_prop(tre2f,dz2,dir);
+ assert(pstat.success());
+ cout << " forward: " << tre2f << endl;
+ cout << pstat << endl;
+ if(dir == Propagator::FORWARD)
+ assert( pstat.forward() );
+ if(dir == Propagator::BACKWARD)
+ assert( pstat.backward() );
+
+ //
+ // Propagate using PropZZ and check difference.
+ ETrack tre2f0 = tre1;
+ pstat = prop0.err_dir_prop(tre2f0,dz2,dir);
+ assert(pstat.success());
+ cout << " forward: " << tre2f0 << endl;
+ cout << pstat << endl;
+ double vdiff0 =
+ dz2.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,dz1,rdir);
+ assert(pstat.success());
+ cout << " f return: " << tre2fb << endl;
+ cout << pstat << endl;
+ if(rdir == Propagator::FORWARD)
+ assert( pstat.forward() );
+ if(rdir == Propagator::BACKWARD)
+ assert( pstat.backward() );
+ double difff =
+ dz1.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 PropZZRK("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( "PropZZRK ",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 PropZZRK type=\"constant\" precision=1.e-6 dir=\"forward\" ]";
+ istringstream isstrm(istring);
+ {
+ StdObjStream obstr(isstrm);
+ string name = obstr.read_object();
+ assert( name == "obj2" );
+ const PropZZRK* pobj;
+ ObjTable::get_object(name,pobj);
+ assert( pobj != 0 );
+ assert( pobj->get_type() == PropZZRK::get_static_type() );
+ assert( pobj->get_direction() == Propagator::FORWARD );
+ }
+ }
+
+ //********************************************************************
+
+ cout << ok_prefix
+ << "------------- All tests passed. -------------" << endl;
+ return 0;
+
+ //********************************************************************
+
+}