#include <iostream>
#include <TLorentzVector.h>
#include "const.h"
#include "spin.h"

using namespace std;

Spin::Spin() {
  spin_transport.ResizeTo(3,3);
  cosy_file = "";
}

Spin::~Spin() {
}

void Spin::print() {
  spin_transport.Print();
}

void Spin::load_cosy_file(TString filename) {
  if (filename != cosy_file) {
    // code which reads in the COSY spin file.
    char coef1s[100];
    char coef2s[100];
    char coef3s[100];
    char junkstring[100];
    char p1s[20];
    char p2s[20];
    char p3s[20];
    char p4s[20];
    char p5s[20];
    char zip[20];
     
    int numlines[3];
    numlines[0]=252;
    numlines[1]=247;
    numlines[2]=252;
    
    cosy_file = filename;
    //char loadfile[20] = cosy_file;
    FILE *spfi;
    spfi = fopen (cosy_file.Data(), "r");
    if (spfi == NULL) perror ("Error Opening COSY File");
    else {
      cout << "Opened COSY file successfully ... " << endl;
      for (int j = 0; j<3; j++) {
        for (int i = 0; i < numlines[j]; i++) {
          Int_t k = j*252 + i;
          fgets(coef1s,22,spfi);
          fgets(coef2s,21,spfi);
          fgets(coef3s,21,spfi);
          fgets(zip,2,spfi);
          fgets(p1s,2,spfi);
          fgets(p2s,2,spfi);
          fgets(p3s,2,spfi);
          fgets(p4s,2,spfi);
          fgets(zip,2,spfi);
          fgets(p5s,2,spfi);
          //fgets(zip,3,spfi);
          zip2 = fgetc(spfi); // not used
          coef1[k] = strtod(coef1s,NULL);
          coef2[k] = strtod(coef2s,NULL);
          coef3[k] = strtod(coef3s,NULL);
	  //cout << " =========== reading cosy file ========== ";
	  //cout << k << " " << coef1[k] << " " << coef2[k] << " " << coef3[k] << endl;
          p1[k] = atoi(p1s);
          p2[k] = atoi(p2s);
          p3[k] = atoi(p3s);
          p4[k] = atoi(p4s);
          p5[k] = atoi(p5s);
	  //cout << p1[k] << " " << p2[k] << " " << p3[k] << " " << p4[k] << " " << p5[k] << endl;
        }
        fgets(junkstring,100, spfi);
      }
    }
    fclose(spfi);
  }
}

// ----------------------------------------------------------------------------

void Spin::set_to_unit_matrix() {
  spin_transport.UnitMatrix();
}

// ----------------------------------------------------------------------------

void Spin::lab_to_tra_EP(Double_t thLeft, Double_t Th_tgt_p, Double_t Ph_tgt_p) {
  Double_t Psi_tgt_p = atan(tan(Th_tgt_p) * cos(Ph_tgt_p));
  Double_t cos_phi_p = cos(Ph_tgt_p);
  Double_t sin_phi_p = sin(Ph_tgt_p);
  Double_t cos_psi_p = cos(-Psi_tgt_p);
  Double_t sin_psi_p = sin(-Psi_tgt_p);

  // The unit vectors are constructed in the proton-arm transport system
  //    k : incident electron beam direction
  //    q : proton momentum (= three momentum-transfer vector) direction

  TVector3 k(0, -sin(thLeft), cos(thLeft));
  TVector3 q(sin_psi_p, cos_psi_p * sin_phi_p, cos_psi_p * cos_phi_p);
  TVector3 y = (q.Cross(k)).Unit();
  TVector3 x = y.Cross(q);

  TMatrixD rot(3,3);

  rot(0,0) = x(0); rot(0,1) = y(0); rot(0,2) = q(0);
  rot(1,0) = x(1); rot(1,1) = y(1); rot(1,2) = q(1);
  rot(2,0) = x(2); rot(2,1) = y(2); rot(2,2) = q(2);
//
// For testing purposes
//
//  rot(0,0) = 0; rot(0,1) = -1; rot(0,2) = 0;
//  rot(1,0) = 1; rot(1,1) = 0; rot(1,2) = 0;
//  rot(2,0) = 0; rot(2,1) = 0; rot(2,2) = 1;

//  cout << "==== EP scattering ========" << endl;
//  for(int i=0;i<3;i++){
// 	for(int j=0;j<3;j++){
//		cout << rot(i,j) << "    ";
//	}
//	cout << endl;
//  }
//  cout << "==== end of EP scattering  ========" << endl;

  rot *= spin_transport;
  spin_transport = rot;
}

// ----------------------------------------------------------------------------

void Spin::lab_to_tra_EEP(Double_t e0, Double_t p_e,
                          Double_t thRight, Double_t thLeft,
                          Double_t Th_tgt_e, Double_t Ph_tgt_e) {

  Double_t Psi_tgt_e = atan(tan(Th_tgt_e) * cos(Ph_tgt_e));
  Double_t cos_phi_e = cos(Ph_tgt_e);
  Double_t sin_phi_e = sin(Ph_tgt_e);
  Double_t cos_psi_e = cos(-Psi_tgt_e);
  Double_t sin_psi_e = sin(-Psi_tgt_e);

  // ki and kf are the initial and final electron-momentum vectors
  // in the electron arm transport frame

  TVector3 ki(0, -e0*sin(thRight), e0*cos(thRight));
  TVector3 kf(p_e*sin_psi_e, -p_e*cos_psi_e * sin_phi_e, p_e*cos_psi_e * cos_phi_e);
  TVector3 y = (ki.Cross(kf)).Unit();
  TVector3 z = (ki - kf).Unit();
  TVector3 x = y.Cross(z);

  TMatrixD rot(3,3);

  rot(0,0) = x(0); rot(0,1) = y(0); rot(0,2) = z(0);
  rot(1,0) = x(1); rot(1,1) = y(1); rot(1,2) = z(1);
  rot(2,0) = x(2); rot(2,1) = y(2); rot(2,2) = z(2);

//  cout << "==== first one ========" << endl;
//  for(int i=0;i<3;i++){
//	for(int j=0;j<3;j++){
//		cout << rot(i,j) << "    ";
//	}
//	cout << endl;
//  }
//  cout << "==== second one ========" << endl;

  rot *= spin_transport;
  spin_transport = rot;

  Double_t c = cos(thLeft-thRight);
  Double_t s = sin(thLeft-thRight);

  rot(0,0) = 1.; rot(0,1) = 0; rot(0,2) =  0;
  rot(1,0) = 0.; rot(1,1) = c; rot(1,2) = -s;
  rot(2,0) = 0.; rot(2,1) = s; rot(2,2) =  c;

//  for(int i=0;i<3;i++){
//	for(int j=0;j<3;j++){
//		cout << rot(i,j) << "    ";
//	}
//	cout << endl;
//  }
//  cout << "===== third one =======" << endl;

  rot *= spin_transport;
  spin_transport = rot;

//  for(int i=0;i<3;i++){
//	for(int j=0;j<3;j++){
//		cout << rot(i,j) << "    ";
//	}
//	cout << endl;
//  }
//  cout << "====================" << endl;

}

// ----------------------------------------------------------------------------

void Spin::spin_dipole(Double_t precession_angle) {
    TMatrixD rot(3,3);

    Double_t c = cos(precession_angle);
    Double_t s = sin(precession_angle);

//    rot(0,0) =  -1; rot(0,1) = 0; rot(0,2) =  0;
//    rot(1,0) =  0; rot(1,1) = c; rot(1,2) =  -s;
//    rot(2,0) =  0; rot(2,1) = s; rot(2,2) = c;

//
// The following are some matrices for testing purposes
// 
// Hall A 
//
//    rot(0,0) =  0; rot(0,1) = 1; rot(0,2) =  0;
//    rot(1,0) =  c; rot(1,1) = 0; rot(1,2) = -s;
//    rot(2,0) =  s; rot(2,1) = 0; rot(2,2) =  c;
//
// Hall C
//
    rot(0,0) =  c; rot(0,1) = 0; rot(0,2) =  s;
    rot(1,0) =  0; rot(1,1) = 1; rot(1,2) =  0;
    rot(2,0) = -s; rot(2,1) = 0; rot(2,2) =  c;
//
//    rot(0,0) =  1; rot(0,1) = 0; rot(0,2) =  0;
//    rot(1,0) =  0; rot(1,1) = 0.99999; rot(1,2) =  0.00001;
//    rot(2,0) =  0; rot(2,1) = 0.00001; rot(2,2) =  0.99999;
//
//    rot(0,0) =  1.; rot(0,1) = 0.; rot(0,2) =  0.;
//    rot(1,0) =  0.; rot(1,1) = 0.707; rot(1,2) =  -0.707;
//    rot(2,0) =  0.; rot(2,1) = 0.707; rot(2,2) =  0.707;
//
//    rot(0,0) =  1; rot(0,1) = 0; rot(0,2) =  0;
//    rot(1,0) =  0; rot(1,1) = 0.017; rot(1,2) =  -0.999;
//    rot(2,0) =  0; rot(2,1) = 0.999; rot(2,2) =  0.017;
  
//  for (int i=0;i<3;i++){
//	for(int j=0;j<3;j++){
//		cout << rot(i,j) << "    ";
//	}
//	cout << endl;
//  }
//  cout << "====================" << endl;

  rot *= spin_transport;
  spin_transport = rot;
}

// ----------------------------------------------------------------------------

void Spin::spin_cosy(Double_t pLeft, Double_t y_tgt_p,
                     Double_t Th_tgt_p, Double_t Ph_tgt_p,
                     Double_t p_p) {

  TMatrixD rot(3,3);
  TMatrixD rottemp(3,3);

  // here goes the code which calculates the matrix elements of the COSY spin-transport matrix

  //need momentum and mass in MeV;
  
  // for testing purposes
  //p_p = pLeft*1.03;
  //Ph_tgt_p = 0.01;
  //Th_tgt_p = 0.01;
  //y_tgt_p = 1.00;
 
  Double_t p_pMeV = p_p*1000;
  Double_t m_pMeV = m_p*1000;
  Double_t pLeftMeV = pLeft*1000;

  Double_t tan_ph = tan(Ph_tgt_p);
  Double_t tan_th = tan(Th_tgt_p);

  //other fill variables

  Double_t cosy_vect[5];
  Double_t fillarray[9];

  // Setting up vector to change fpp variables to cosy variable

  // x     -->     x
  // th    -->     a
  // y     -->     y
  // ph    -->     b
  // z     -->     tof (always zero)
  // delta -->     kin

  cosy_vect[0] = 0;
  cosy_vect[1] = (tan_th*p_pMeV)/(pLeftMeV*sqrt(tan_th*tan_th + tan_ph*tan_ph +1));
  cosy_vect[2] = y_tgt_p/100.0;
  cosy_vect[3] = (tan_ph*p_pMeV)/(pLeftMeV*sqrt(tan_th*tan_th + tan_ph*tan_ph +1));

  double tkin = sqrt(p_pMeV*p_pMeV + m_pMeV*m_pMeV) - m_pMeV;
  double tkin_cent = sqrt(pLeftMeV*pLeftMeV + m_pMeV*m_pMeV) - m_pMeV;

  cosy_vect[4] = (tkin - tkin_cent)/tkin_cent;

  // for testing purposes
  //cosy_vect[0] = 0;
  //cosy_vect[1] = 0;
  //cosy_vect[2] = 0;
  //cosy_vect[3] = 0;
  //cosy_vect[4] = 0;

  
//  cout << "COSY variables: " << tkin << " " << tkin_cent << endl;
//  cout << "COSY variables: " << p_pMeV << " " << pLeftMeV << endl;
//  cout << "COSY variables: " << m_pMeV << " " << y_tgt_p << endl;
//  cout << "COSY variables: " << Ph_tgt_p << " " << Th_tgt_p << endl;
//  cout << "COSY variables: " << tan_ph << " " << tan_th << endl;

  //cout << "COSY Vector" << endl;
  // for (int i=0; i<5; i++){cout << i << " " << cosy_vect[i] << endl;};

  // read in cosy matrix file
  //
  int numlines[3];
  numlines[0]=252;
  numlines[1]=247;
  numlines[2]=252;
  
  for (Int_t j =0; j <3; j++){
    fillarray[j] = 0;
    fillarray[j+3] = 0;
    fillarray[j+6] = 0;
    for (Int_t i = 0; i < numlines[j]; i++){
      Int_t k2 = j * 252 + i;
      //cout << "========= fill array ==========" << endl;
      fillarray[j] += coef1[k2] * pow(cosy_vect[0],p1[k2]) * pow(cosy_vect[1],p2[k2]) *
          pow(cosy_vect[2],p3[k2]) * pow(cosy_vect[3],p4[k2]) * pow(cosy_vect[4],p5[k2]);
        
      fillarray[j+3] += coef2[k2] * pow(cosy_vect[0],p1[k2]) * pow(cosy_vect[1],p2[k2]) *
          pow(cosy_vect[2],p3[k2]) * pow(cosy_vect[3],p4[k2]) * pow(cosy_vect[4],p5[k2]);
        
      fillarray[j+6] += coef3[k2] * pow(cosy_vect[0],p1[k2]) * pow(cosy_vect[1],p2[k2]) *
          pow(cosy_vect[2],p3[k2]) * pow(cosy_vect[3],p4[k2]) * pow(cosy_vect[4],p5[k2]);
      
      //cout << coef1[k2] << " " << coef2[k2] << " " << coef3[k2] << endl;
      //cout << p1[k2] << " " << p2[k2] << " " << p3[k2] << " " <<
//	     p4[k2] << " " << p5[k2] << endl;
      //cout << "term1 = " << coef1[k2] * pow(cosy_vect[0],p1[k2]) * pow(cosy_vect[1],p2[k2]) * pow(cosy_vect[2],p3[k2]) * pow(cosy_vect[3],p4[k2]) * pow(cosy_vect[4],p5[k2]) << endl;  
      //cout << "term2 = " << coef2[k2] * pow(cosy_vect[0],p1[k2]) * pow(cosy_vect[1],p2[k2]) * pow(cosy_vect[2],p3[k2]) * pow(cosy_vect[3],p4[k2]) * pow(cosy_vect[4],p5[k2]) << endl;  
      //cout << "term3 = " << coef3[k2] * pow(cosy_vect[0],p1[k2]) * pow(cosy_vect[1],p2[k2]) * pow(cosy_vect[2],p3[k2]) * pow(cosy_vect[3],p4[k2]) * pow(cosy_vect[4],p5[k2]) << endl;  
    }
  }
  //for(int j=0;j<3;j++){
//	cout << fillarray[j] << "  " << fillarray[j+3] << "  " << fillarray[j+6] << endl;
 // }	
//  cout << "====================" << endl;


  // ...

//    double c=cos(104.5*3.14159265/180.0);
//    double s=sin(104.5*3.14159265/180.0);

//    rot(0,0) =  0; rot(0,1) = 1; rot(0,2) =  0;
//    rot(1,0) =  c; rot(1,1) = 0; rot(1,2) = -s;
//    rot(2,0) =  s; rot(2,1) = 0; rot(2,2) =  c;
 

  // "F"  : column major (Fortran) m[i][j] = array[i+j*fNrows]
  // else : row major    (C)       m[i][j] = array[i*fNcols+j] (default)
  rot.SetMatrixArray(fillarray,"F");

  //If you wish to check the determinant, enable the four lines below, and make public det_cosy_mat.

  //Double_t deta,detb;
  //TDecompLU lu(rot);
  //lu.Det(deta,detb);
  //det_cosy_mat = deta*pow(2,detb);
      
  //rottemp(0,0)=rot(1,0); rot(1,0)=rot(0,0);rot(0,0)=rottemp(0,0); 
  //rottemp(0,1)=rot(1,1); rot(1,1)=rot(0,1);rot(0,1)=rottemp(0,1); 
  //rottemp(0,2)=rot(1,2); rot(1,2)=rot(0,2);rot(0,2)=rottemp(0,2); 
  rot(0,2)=-rot(0,2);
  rot(2,0)=-rot(2,0);
  rottemp(1,0)=rot(2,1); rot(2,1)=rot(0,1); rot(0,1)=rottemp(1,0); 
  rottemp(1,0)=rot(0,1); rot(0,1)=rot(1,0); rot(1,0)=rottemp(1,0); 
  rottemp(2,1)=rot(1,2); rot(1,2)=rot(2,1); rot(2,1)=rottemp(2,1); 
  rot(0,1)=-rot(0,1);

//  for(int i=0;i<3;i++){
//	for(int j=0;j<3;j++){
//		cout << rot(i,j) << "    ";
//	}
//	cout << endl;
//  }
//  cout << "====================" << endl;

  rot *= spin_transport;
  spin_transport = rot;
}

// ----------------------------------------------------------------------------

void Spin::spin_pentchev(Double_t p_p, Double_t y_tgt_p,
                         Double_t Th_tgt_p, Double_t Ph_tgt_p,
                         Double_t Th_tra_p, Double_t Ph_tra_p,
                         Int_t mode) {
  /*
     This rotation is a replacement of the COSY rotation,
     based on an approximation by Lubomir Pentchev  (see
     JLAB TN-03-024).
  */

  double gamma = sqrt(p_p*p_p + m_p*m_p)/m_p;

  double fcPhTgtP = Ph_tgt_p;

  const double phi_y = -0.317;
  const double phi_phi = -0.066;
  const double dphi_y = 0.218;
  const double dphi_phi = -0.287;

  double PhiD = phi_y*y_tgt_p/1000. + phi_phi*fcPhTgtP;
  double deltaPhiD = dphi_y*y_tgt_p/1000. + dphi_phi*fcPhTgtP;

  double ChiTheta =  (BendAngle*DEGTORAD + Th_tgt_p - Th_tra_p)
    * kappa * gamma;
  double ChiPhi = gamma*kappa*((Ph_tra_p*cos(ChiTheta)-fcPhTgtP)
                               + PhiD*(1-cos(ChiTheta))
               - deltaPhiD*(1-2*sin(ChiTheta)/ChiTheta + cos(ChiTheta)));
  double ChiPhiPrime = gamma*kappa*(fcPhTgtP*sin(ChiTheta) - PhiD*sin(ChiTheta)
                                    + deltaPhiD*(2*(1-cos(ChiTheta))/ChiTheta - sin(ChiTheta)));

// Phil: why create the matrices in such an obfuscated way?
// TODO: rewrite this after I figure out what cases 1 and 2 are for
// cout << "====== start =============" << endl;
// for(int i=0;i<3;i++){
//	for(int j=0;j<3;j++){
//		cout << spin_transport(i,j) << "    ";
//	}
//	cout << endl;
//  }
//  cout << "====== next ==========" << endl;
  
  switch (mode) {
    case 0: {
      Double_t cChiPhi = cos(ChiPhi);
      Double_t sChiPhi = sin(ChiPhi);
      Double_t cChiPhiPrime = cos(ChiPhiPrime);
      Double_t sChiPhiPrime = sin(ChiPhiPrime);
      Double_t cChiTheta = cos(ChiTheta);
      Double_t sChiTheta = sin(ChiTheta);
  
      TMatrixD rot(3,3);
      rot(0,0) = cChiPhiPrime * cChiTheta;
      rot(1,0) = -sChiPhiPrime;
      rot(2,0) = -cChiPhiPrime * sChiTheta;
      //rot(2,0) = cChiPhiPrime * sChiTheta;
  
      rot(0,1) = -sChiPhi * cChiPhiPrime * sChiTheta + cChiPhi * sChiPhiPrime * cChiTheta;
      rot(1,1) = cChiPhi * cChiPhiPrime;
      rot(2,1) = -cChiPhi * sChiPhiPrime * sChiTheta - sChiPhi * cChiPhiPrime * cChiTheta;
  
      rot(0,2) = cChiPhi * sChiTheta;
      //rot(0,2) = -cChiPhi * sChiTheta;
      rot(1,2) = sChiPhi;
      rot(2,2) = cChiPhi * cChiTheta;

//      for(int i=0;i<3;i++){
//	for(int j=0;j<3;j++){
//		cout << rot(i,j) << "    ";
//	}
//	cout << endl;
//      }
//      cout << "====================" << endl;
      
      rot *= spin_transport;
      spin_transport = rot;
      
//      for(int i=0;i<3;i++){
//	for(int j=0;j<3;j++){
//		cout << spin_transport(i,j) << "    ";
//	}
//	cout << endl;
//      }
//      cout << "====================" << endl;
  
      break;
    }
    case 1:
    case 2: {
      TMatrixD rot1(3,3);
      TMatrixD rot2(3,3);
      TMatrixD rot(3,3);
      TMatrixD rottemp(3,3);
      rot(0,0) = 1.; rot(0,1) = 0; rot(0,2) =  0;
      rot(1,0) = 0.; rot(1,1) = 1; rot(1,2) =  0;
      rot(2,0) = 0.; rot(2,1) = 0; rot(2,2) =  1;
      Double_t c = cos(-ChiPhi);
      Double_t s = sin(-ChiPhi);
      rot1(0,0) = 1.; rot1(0,1) = 0; rot1(0,2) =  0;
      rot1(1,0) = 0.; rot1(1,1) = c; rot1(1,2) = -s;
      rot1(2,0) = 0.; rot1(2,1) = s; rot1(2,2) =  c;
      c = cos(ChiTheta);
      s = sin(ChiTheta);
      //rot2(0,0) =  c; rot2(0,1) = 0; rot2(0,2) =  s;
      rot2(0,0) =  c; rot2(0,1) = 0; rot2(0,2) =  -s;
      rot2(1,0) =  0; rot2(1,1) = 1; rot2(1,2) =  0;
      //rot2(2,0) = -s; rot2(2,1) = 0; rot2(2,2) =  c;
      rot2(2,0) = s; rot2(2,1) = 0; rot2(2,2) =  c;
      if (mode == 1) {
        rot1 *= rot;
        rot = rot1;
        rot2 *= rot;
        rot = rot2;
      } else {
        rot2 *= rot;
        rot = rot2;
        rot1 *= rot;
        rot = rot1;
      }

//      rottemp(0,0)=rot(1,0); rot(1,0)=rot(0,0);rot(0,0)=rottemp(0,0); 
//      rottemp(0,1)=rot(1,1); rot(1,1)=rot(0,1);rot(0,1)=rottemp(0,1); 
//      rottemp(0,2)=rot(1,2); rot(1,2)=rot(0,2);rot(0,2)=rottemp(0,2); 

//      for(int i=0;i<3;i++){
//	for(int j=0;j<3;j++){
//		cout << rot(i,j) << "    ";
//	}
//	cout << endl;
//      }
//      cout << "======= end ===========" << endl;
  
      rot *= spin_transport;
      spin_transport = rot;

//      for(int i=0;i<3;i++){
//	for(int j=0;j<3;j++){
//		cout << spin_transport(i,j) << "    ";
//	}
//	cout << endl;
//      }
//      cout << "====================" << endl;
  
      break;
    }
    default:
      cout << " " << endl;
      exit(1);
  }
}

// ----------------------------------------------------------------------------

void Spin::tra_to_fpp(Double_t Th_tra_p, Double_t Ph_tra_p) {
  // this is equal to PALM
  Double_t Psi_tra_p = atan(tan(Th_tra_p) * cos(Ph_tra_p));

  TMatrixD rot(3,3);

  Double_t c = cos(Ph_tra_p);
  Double_t s = sin(Ph_tra_p);

  rot(0,0) = 1.; rot(0,1) = 0; rot(0,2) =  0;
  rot(1,0) = 0.; rot(1,1) = c; rot(1,2) = -s;
  rot(2,0) = 0.; rot(2,1) = s; rot(2,2) =  c;
//
// For testing purposes
//
//  rot(0,0) = 1.; rot(0,1) = 0; rot(0,2) =  0;
//  rot(1,0) = 0.; rot(1,1) = 1; rot(1,2) =  0;
//  rot(2,0) = 0.; rot(2,1) = 0; rot(2,2) =  1;

  rot *= spin_transport;
  spin_transport = rot;

  c = cos(-Psi_tra_p);
  s = sin(-Psi_tra_p);

  rot(0,0) =  c; rot(0,1) = 0; rot(0,2) =  s;
  rot(1,0) =  0; rot(1,1) = 1; rot(1,2) =  0;
  rot(2,0) = -s; rot(2,1) = 0; rot(2,2) =  c;
//
// For testing purposes
//
//  rot(0,0) =  1; rot(0,1) = 0; rot(0,2) =  0;
//  rot(1,0) =  0; rot(1,1) = 1; rot(1,2) =  0;
//  rot(2,0) =  0; rot(2,1) = 0; rot(2,2) =  1;
 
  rot *= spin_transport;
  spin_transport = rot;
  
//  cout << "======== tra_to_fpp =======" << endl;
//  for(int i=0;i<3;i++){
//	for(int j=0;j<3;j++){
//  	cout << spin_transport(i,j) << "        ";
//	}
//	cout << endl;
//  }
//  cout << "======== end tra_to_fpp =======" << endl;

}

// ----------------------------------------------------------------------------

Double_t &Spin::operator()(int rown, int coln) {
  return spin_transport(rown, coln);
}