// GtrTrfPropagator.cpp

#include <iostream>
#include "GtrTrfPropagator.hpp"
#include "trfbase/PropDispatch.h"
#include "trfbase/PropStat.h"
#include "trfbase/PropNull.h"
#include "trfcyl/SurfCylinder.h"
#include "trfcyl/PropCyl.h"
#include "trfcyl/PropJoinCyl.h"
#include "trfzp/SurfZPlane.h"
#include "trfzp/PropZZ.h"
#include "trfxyp/SurfXYPlane.h"
#include "trfxyp/PropXYXY.h"
#include "trfcylplane/PropZXY.h"
#include "trfcylplane/PropXYZ.h"
#include "trfcylplane/PropZCyl.h"
#include "trfcylplane/PropXYCyl.h"
#include "trfcylplane/PropCylZ.h"
#include "trfcylplane/PropCylXY.h"
#include "trfdca/SurfDCA.h"
#include "trfdca/PropCylDCA.h"
#include "trfdca/PropDCACyl.h"
#include "trfbase/ObjDoublePtr.h"
#include "gtrbase/GtrFieldPkg.hpp"
#include "objstream/ObjTable.hpp"

using namespace trf;
using std::cerr;
using std::endl;

//**********************************************************************

// Default Constructor.

GtrTrfPropagator::GtrTrfPropagator() {

  // Extract the types.
  TrfObject::Type tcyl = SurfCylinder::get_static_type();
  TrfObject::Type tzpl = SurfZPlane::get_static_type();
  TrfObject::Type txyp = SurfXYPlane::get_static_type();
  TrfObject::Type tdca = SurfDCA::get_static_type();

  // Get global magnetic field value.
  GtrFieldPkg::initialize();
  ObjDoublePtr bp;
  ObjTable::Status stat = ObjTable::get_object("bfield", bp);
  if(stat != ObjTable::OK) {
    cerr << "GtrTrfPropagator: Can't get global bfield value." << endl;
    abort();
  }

  // Construct the constituent propagators.
  double rmin = 1.0;
  double rfac = 1.1;
  PropagatorPtr pcylcyl( new PropCyl(bp));
  PropagatorPtr pzplzpl( new PropZZ(bp));
  PropagatorPtr pxypxyp( new PropXYXY(bp));
  PropagatorPtr pzplxyp( new PropZXY(bp));
  PropagatorPtr pxypzpl( new PropXYZ(bp));
  PropagatorPtr pcylzpl( new PropCylZ(bp));
  PropagatorPtr pcylxyp( new PropCylXY(bp));
  PropagatorPtr pzplcyl( new PropZCyl(bp));
  PropagatorPtr pxypcyl( new PropXYCyl(bp));
  PropagatorPtr pcyldca( new PropCylDCA(bp));
  PropagatorPtr pdcacyl( new PropDCACyl(bp));
  PropagatorPtr pzpldca( new PropJoinCyl(rmin,1.0,pzplcyl,pcyldca) );
  PropagatorPtr pxypdca( new PropJoinCyl(rmin,1.0,pxypcyl,pcyldca) );
  PropagatorPtr pdcazpl( new PropJoinCyl(rmin,rfac,pdcacyl,pcylzpl) );
  PropagatorPtr pdcaxyp( new PropJoinCyl(rmin,rfac,pdcacyl,pcylxyp) );
  PropagatorPtr pdcadca( new PropNull );

  // Construct the dispatch propagator.
  PropDispatch* _pdispatch = new PropDispatch;
  _pprop = _pdispatch;
  _pdispatch->add_propagator(tcyl,tcyl,pcylcyl);
  _pdispatch->add_propagator(tzpl,tzpl,pzplzpl);
  _pdispatch->add_propagator(txyp,txyp,pxypxyp);
  _pdispatch->add_propagator(tzpl,txyp,pzplxyp);
  _pdispatch->add_propagator(txyp,tzpl,pxypzpl);
  _pdispatch->add_propagator(tcyl,tzpl,pcylzpl);
  _pdispatch->add_propagator(tcyl,txyp,pcylxyp);
  _pdispatch->add_propagator(tzpl,tcyl,pzplcyl);
  _pdispatch->add_propagator(txyp,tcyl,pxypcyl);
  _pdispatch->add_propagator(tcyl,tdca,pcyldca);
  _pdispatch->add_propagator(tdca,tcyl,pdcacyl);
  _pdispatch->add_propagator(tzpl,tdca,pzpldca);
  _pdispatch->add_propagator(txyp,tdca,pxypdca);
  _pdispatch->add_propagator(tdca,tzpl,pdcazpl);
  _pdispatch->add_propagator(tdca,txyp,pdcaxyp);
  _pdispatch->add_propagator(tdca,tdca,pdcadca);

}

//**********************************************************************