#include #include #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); }