Fork me on GitHub

Ignore:
Timestamp:
Apr 16, 2014, 3:56:14 PM (10 years ago)
Author:
Pavel Demin
Message:

switch to a more stable Hector version

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/external/Hector/H_BeamParticle.cc

    r1360 r1365  
    1   /* * * * * * * * * * * * * * * * * * * * * * * * * * * *
    2  *                                                         *
    3 *                   --<--<--  A fast simulator --<--<--     *
    4 *                 / --<--<--     of particle   --<--<--     *
    5 *  ----HECTOR----<                                          *
    6 *                 \ -->-->-- transport through -->-->--     *
    7 *                   -->-->-- generic beamlines -->-->--     *
    8 *                                                           *
    9 * JINST 2:P09005 (2007)                                     *
    10 *      X Rouby, J de Favereau, K Piotrzkowski (CP3)         *
    11 *       http://www.fynu.ucl.ac.be/hector.html               *
    12 *                                                           *
    13 * Center for Cosmology, Particle Physics and Phenomenology  *
    14 *              Universite catholique de Louvain             *
    15 *                 Louvain-la-Neuve, Belgium                 *
    16  *                                                         *
    17    * * * * * * * * * * * * * * * * * * * * * * * * * * * */
     1/*
     2---- Hector the simulator ----
     3   A fast simulator of particles through generic beamlines.
     4   J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be
     5
     6        http://www.fynu.ucl.ac.be/hector.html
     7
     8   Centre de Physique des Particules et de Phénoménologie (CP3)
     9   Université Catholique de Louvain (UCL)
     10*/
    1811
    1912/// \file H_BeamParticle.cc
    … …  
    2619// c++ #includes
    2720#include <iostream>
    28 #include <fstream>
    2921#include <iomanip>
    3022#include <vector>
    … …  
    3527#include "H_Parameters.h"
    3628#include "TRandom.h"
    37 #include "TVectorD.h"
     29//#include "TView.h"
     30//#include "TPolyLine3D.h"
    3831#ifdef _include_pythia_
    3932#include "TPythia6.h"
     33#include "TRandom.h"
    4034#endif
    41 
    4235// local #includes
    4336#include "H_OpticalElement.h"
    … …  
    4740using namespace std;
    4841
    49 void H_BeamParticle::init() { // for more efficiency, put the objects away from init!
     42void H_BeamParticle::init() {
    5043        mp = MP;
    5144        qp = QP;
    … …  
    6558}
    6659
    67 H_BeamParticle::H_BeamParticle() { 
     60H_BeamParticle::H_BeamParticle() {
    6861        init();
    6962}
    7063
    71 H_BeamParticle::H_BeamParticle(const H_BeamParticle& p):
    72         mp(p.mp), qp(p.qp), fs(p.fs), fx(p.fx), fy(p.fy), thx(p.thx), thy(p.thy),
    73         energy(p.energy), hasstopped(p.hasstopped), hasemitted(p.hasemitted),
    74         isphysical(p.isphysical), stop_position(new TVectorD(*(p.stop_position))),
    75         stop_element(0), positions(p.positions) {
     64H_BeamParticle::H_BeamParticle(const H_BeamParticle& p) {
     65        mp = p.mp;
     66        qp = p.qp;
     67        fx = p.fx;
     68        fy = p.fy;
     69        thx = p.thx;
     70        thy = p.thy;
     71        fs = p.fs;
     72        energy = p.energy;
     73        hasstopped = p.hasstopped;
     74        hasemitted = p.hasemitted;
     75        isphysical = p.isphysical;
     76        stop_position = new TVectorD(*(p.stop_position));
    7677        if(p.hasstopped) stop_element = new H_OpticalElement(*(p.stop_element));
     78        positions = p.positions;
    7779}
    7880
    … …  
    99101        stop_position = new TVectorD(*(p.stop_position));
    100102        if(p.hasstopped) stop_element = new H_OpticalElement(*(p.stop_element));
    101         else stop_element = 0;
    102103        positions = p.positions;
    103104        return *this;
    104105}
    105106
    106 const bool H_BeamParticle::stopped(const H_AbstractBeamLine * beamline) {
     107bool H_BeamParticle::stopped(const H_AbstractBeamLine * beamline) {
    107108        vector<TVectorD>::const_iterator position_i;
    108109        for(position_i = positions.begin(); position_i < positions.end()-1; position_i++) {
    … …  
    110111                if(beamline->getElement(pos)->getAperture()->getType()!=NONE) {
    111112                        bool has_passed_entrance = beamline->getElement(pos)->isInside((*position_i)[INDEX_X],(*position_i)[INDEX_Y]);
    112                         bool has_passed_exit     = beamline->getElement(pos)->isInside((*(position_i+1))[INDEX_X],(*(position_i+1))[INDEX_Y]);
    113                         // here we should distinguish between particles passing the input or not.
    114                         //  - particles not passing the input are logically stopped at the input position. period.
    115                         //  - particles passing the input but not the output are stopped somewhere in the element
    116                         //  - particles passing both the input and the output could have been stopped somewhere inside too (less likely)
    117                         if(!has_passed_entrance) {
    118                                 if(VERBOSE) cout<<"Stopped at the entrance of "<<beamline->getElement(pos)->getName();
    119                                 hasstopped=true;
    120                                 stop_element = const_cast<H_OpticalElement*>(beamline->getElement(pos));
    121                                 *stop_position = *position_i;
    122                                 if(VERBOSE) cout<<" at s = "<<(*stop_position)[4]<<endl;
    123                                 return hasstopped;
    124                         } else if (!has_passed_exit) {
    125                                 if(VERBOSE) cout<<"Stopped inside "<<beamline->getElement(pos)->getName();
    126                                 hasstopped=true;
    127                                 stop_element = const_cast<H_OpticalElement*>(beamline->getElement(pos));
    128                                 // this should be computed using the element-based method "H_OpticalElement::getHitPosition" (caution : always nonlinear).
    129                                 *stop_position = beamline->getElement(pos)->getHitPosition(*position_i,BE-energy,mp,qp);
    130                                 if(VERBOSE) cout<<" at s = "<<(*stop_position)[4]<<endl;
    131                                 return hasstopped;
    132                         }
    133 /*
     113                        bool has_passed_exit     = beamline->getElement(pos)->isInside((*(position_i+1))[INDEX_X],(*(position_i+1))[INDEX_Y]);
    134114                        if(!(has_passed_entrance && has_passed_exit)) {
     115        //                      cout << "p =" << (*position_i)[INDEX_X] << "; el=" << beamline->getElement(pos)->getX()<<endl;
     116        //                      beamline->getElement(position_i-positions.begin())->printProperties();
    135117                                if(VERBOSE) cout<<"particle stopped at "<<(beamline->getElement(pos))->getName();
    136118                                if(VERBOSE) cout<<" (s = "<<(*position_i)[4] << ")" << endl;
     119                                if(VERBOSE && !has_passed_exit) cout << "Particle stopped inside the element" << endl;
    137120                                hasstopped=true;
    138121                                stop_element = const_cast<H_OpticalElement*>(beamline->getElement(pos));
    … …  
    140123                                return hasstopped;
    141124                        } // if
    142 */
     125//                      else cout << "outside aperture " << endl;
    143126                } // if
    144127        } // for
    … …  
    159142}
    160143
    161 
    162 
    163144void H_BeamParticle::smearPos(const double dx,const double dy, TRandom* r) {
    164145  // the beam is centered on (fx,fy) at IP
    … …  
    191172}
    192173
    193 
    194174void H_BeamParticle::set4Momentum(const double px, const double py, const double pz, const double ene) {
    195175        /// @param px, py, pz, ene is \f$ (\vec p , E) [GeV]\f$
    … …  
    207187        addPosition(fx,thx,fy,thy,fs);
    208188        return;
    209 }
    210 
    211 void H_BeamParticle::set4Momentum(const TLorentzVector& pmu) {
    212         /// @param pmu is the particle 4-momentum \f$ p^\mu \f$
    213         ///
    214         /// Clears the H_BeamParticle::positions vector.
    215         set4Momentum(pmu.Px(), pmu.Py(), pmu.Pz(), pmu.E());
    216189}
    217190
    … …  
    242215}
    243216
     217void H_BeamParticle::emitGamma(const double gee, const double gq2) {
     218        emitGamma(gee,gq2,0,2*PI);
     219        return;
     220}
     221
    244222void H_BeamParticle::emitGamma(const double gee, const double gq2, const double phimin, const double phimax) {
    245223        /// @param gee = \f$ E_{\gamma} \f$ is the photon energy
    … …  
    271249       
    272250        if( (gq2>q2max) || (gq2<q2min)) {
    273                 if(VERBOSE) cout<<"<H_BeamParticle> WARNING : Non physical particle ! Q2 (" << q2 << ") and E ("<<gee << ") are not compatible." << endl;
     251                if(VERBOSE) cout<<"Non physical particle ! Q2 (" << q2 << ") and E ("<<gee << ") are not compatible." << endl;
    274252                isphysical = false;
    275253        }
    276254
    277         if(hasemitted) { cout<<"<H_BeamParticle> WARNING : particle has already emitted at least one gamma !"<<endl;}
     255        if(hasemitted) { cout<<"particle has already emitted at least one gamma !"<<endl;}
    278256        hasemitted = true;
    279257        energy = energy - gee;
    … …  
    330308}
    331309
    332 std::ostream& operator<< (std::ostream& os, const H_BeamParticle& p) {
    333         os << " M   = " << p.getM()  << "GeV ";
    334         os << " Q   = " << p.getQ()  << "e";
    335         os << " fx  = " << p.getX()  << "m   ";
    336         os << " fy  = " << p.getY()  << "m   ";
    337         os << " thx = " << p.getTX() << "rad ";
    338         os << " thy = " << p.getTY() << "rad ";
    339         os << endl;
    340    return os;
     310void H_BeamParticle::printProperties() const {
     311        cout << " M   = " << getM()  << "GeV ";
     312        cout << " Q   = " << getQ()  << "e";
     313        cout << " fx  = " << getX()  << "m   ";
     314        cout << " fy  = " << getY()  << "m   ";
     315        cout << " thx = " << getTX() << "rad ";
     316        cout << " thy = " << getTY() << "rad ";
     317        cout << endl;
     318        return;
    341319}
    342320
    … …  
    378356                        if((*position_i)[INDEX_S]>=position) {
    379357                                if(position_i==positions.begin()) {
    380                                         if(VERBOSE) cout<<"<H_BeamParticle> ERROR : non reachable value"<<endl;
     358                                        if(VERBOSE) cout<<"ERROR : non reachable value"<<endl;
    381359                                        return;
    382360                                }
    383361                                l = (*position_i)[INDEX_S] - (*(position_i-1))[INDEX_S];
    384362                                if(l==0) {
    385                                         if(VERBOSE) cout<<"<H_BeamParticle> WARNING : no luck in choosing position, no propagation done"<<endl;
     363                                        if(VERBOSE) cout<<"WARNING : no luck in choosing position, no propagation done"<<endl;
    386364                                        return;
    387365                                }
    … …  
    396374                position_i = positions.begin();
    397375                cout << "Desired position is : " << position << " & positions.begin() is " << (*position_i)[INDEX_S] << endl;
    398                 cout<<"<H_BeamParticle> ERROR : position not reachable"<<endl; 
     376                cout<<"ERROR : position not reachable"<<endl;   
    399377                return;
    400378        }
    401379}
    402380
    403 /// Caution : do not use this method (obsolete) !!!
     381/// Caution : do not use this method !!!
    404382void H_BeamParticle::propagate(const H_AbstractBeamLine * beam, const H_OpticalElement * element) {
    405383        TMatrixD X(*getV());
    406         X *= beam->getPartialMatrix(element);
     384        X *= *(beam->getPartialMatrix(element));
    407385        fx = URAD*(X.GetMatrixArray())[0];
    408386        thx = URAD*atan((X.GetMatrixArray())[1]);
    … …  
    412390}
    413391
    414 void H_BeamParticle::propagate(const H_AbstractBeamLine * beam, const string& el_name) {
    415         propagate(beam->getElement(el_name)->getS()+beam->getElement(el_name)->getLength());
     392void H_BeamParticle::propagate(const H_AbstractBeamLine * beam, const string el_name) {
     393        propagate(beam->getElement(el_name)->getS());
    416394        return;
    417395}
    … …  
    419397void H_BeamParticle::propagate(const H_AbstractBeamLine * beam) {
    420398        TMatrixD X(*getV());
    421         X  *= beam->getBeamMatrix();
     399        X  *= *(beam->getBeamMatrix());
    422400        fx  = URAD*(X.GetMatrixArray())[0];
    423401        thx = URAD*atan((X.GetMatrixArray())[1]);
    … …  
    436414        return ;
    437415}
    438 
     416/*
    439417TGraph * H_BeamParticle::getPath(const int x_or_y, const int color) const{
    440418        /// @param x_or_y = 0(1) draws the x(y) component;
    441419
    442420        const int N = (int) positions.size();
    443                 int mycolor = color;
    444         if(N<2) cout<<"<H_BeamParticle> WARNING : particle positions not calculated : please run computePath"<<endl;
     421        int mycolor = color;
     422        if(N<2) cout<<"particle positions not calculated : please run computePath"<<endl;
    445423        double * s = new double[N], * graph = new double[N];
    446424
    … …  
    461439}
    462440
    463 void H_BeamParticle::getPath(const int x_or_y, const string& filename) const{
    464         /// @param x_or_y = 0(1) draws the x(y) component;
    465         ofstream outfile(filename.c_str());
    466 
    467         int index;
    468         if(x_or_y==0) {index = INDEX_X;} else {index = INDEX_Y;}
     441TPolyLine3D *  H_BeamParticle::getPath3D(const H_AbstractBeamLine * beam, const bool isfirst, const int color, const int side) const{
     442        const int N = (int) positions.size();
     443        int mycolor = color;
     444        if(N<2) cout<<"WARNING : particle positions not calculated. Run computePath"<<endl;
     445        double * s = new double[N], * graphx = new double[N], * graphy = new double[N];
     446        int direction = (side<0)?-1:1;
     447
    469448
    470449        vector<TVectorD>::const_iterator position_i;
    471450        for(position_i = positions.begin(); position_i < positions.end(); position_i++) {
    472                 outfile << (*position_i)[index] << "\t" << (*position_i)[INDEX_S] << endl;
     451                graphx[(int)(position_i-positions.begin())] = (*position_i)[INDEX_X];
     452                graphy[(int)(position_i-positions.begin())] = (*position_i)[INDEX_Y];
     453                s[(int)(position_i-positions.begin())] = (*position_i)[INDEX_S]*1000*direction;
    473454        }
    474         outfile.close();
    475 }
    476 
     455
     456        float coi[3] = {beam->getLength()*(-1000),-10000,-5000};
     457        float cof[3] = {beam->getLength()*1000,10000,5000};
     458        TView *view = TView::CreateView(11);
     459        view->SetRange(coi[0],coi[1],coi[2],cof[0],cof[1],cof[2]);
     460        TPolyLine3D* ppath = new TPolyLine3D(N,s,graphx,graphy);
     461        ppath->SetLineColor(mycolor);
     462
     463        if(isfirst) {
     464                ppath->Draw();
     465                view->ShowAxis();
     466        } else {
     467                ppath->Draw("same");
     468        }
     469
     470        delete [] s;
     471        delete [] graphx;
     472        delete [] graphy;
     473        return ppath;
     474} // getPath3D
     475*/
     476void H_BeamParticle::computePath(const H_AbstractBeamLine * beam) {
     477        computePath(beam,true);
     478}
    477479
    478480// should be removed later, to keep only computePath(const H_AbstractBeamLine & , const bool)
    … …  
    505507        for (int i=0; i<N; i++) {
    506508                const unsigned pos = i;
    507                 if(fs < beam->getElement(pos)->getS() && fs > beam->getElement(pos-1)->getS()) {
    508                         if(beam->getElement(pos-1)->getType()!=DRIFT) {
    509                                 cout<<"Path starts inside element "<<beam->getElement(pos-1)->getName()<<endl;
    510                         } else {
    511                                 cout<<"Path starts inside unnamed drift "<<endl;
    512                         }
    513                         H_OpticalElement* temp_el = beam->getElement(pos-1)->clone();
    514                         temp_el->setS(fs);
    515                         temp_el->setLength(beam->getElement(pos)->getS() - fs);
    516                         mat[0][0] = mat[0][0] - temp_el->getX();
    517                         mat[0][1] = mat[0][1] - tan(temp_el->getTX());
    518                         mat[0][2] = mat[0][2] - temp_el->getY();
    519                         mat[0][3] = mat[0][3] - tan(temp_el->getTY());
    520                         mat *= temp_el->getMatrix(energy_loss,mp,qp);
    521                         mat[0][0] = mat[0][0] + temp_el->getX();
    522                         mat[0][1] = mat[0][1] + tan(temp_el->getTX());
    523                         mat[0][2] = mat[0][2] + temp_el->getY();
    524                         mat[0][3] = mat[0][3] + tan(temp_el->getTY());
    525                         xys[0] = mat.GetMatrixArray()[0]*URAD;
    526                         xys[1] = atan(mat.GetMatrixArray()[1])*URAD;
    527                         xys[2] = mat.GetMatrixArray()[2]*URAD;
    528                         xys[3] = atan(mat.GetMatrixArray()[3])*URAD;
    529                         xys[4] = temp_el->getS()+temp_el->getLength();
    530                         addPosition(xys[0],xys[1],xys[2],xys[3],xys[4]);
    531                         if (temp_el) delete temp_el;
    532                 }
    533                 if(fs <= beam->getElement(pos)->getS()) {
    534                         mat[0][0] = mat[0][0] - beam->getElement(pos)->getX();
    535                         mat[0][1] = mat[0][1] - tan(beam->getElement(pos)->getTX()/URAD)*URAD;
    536                         mat[0][2] = mat[0][2] - beam->getElement(pos)->getY();
    537                         mat[0][3] = mat[0][3] - tan(beam->getElement(pos)->getTY()/URAD)*URAD;
    538                         mat *= beam->getElement(pos)->getMatrix(energy_loss,mp,qp);
    539                         mat[0][0] = mat[0][0] + beam->getElement(pos)->getX();
    540                         mat[0][1] = mat[0][1] + tan(beam->getElement(pos)->getTX()/URAD)*URAD;
    541                         mat[0][2] = mat[0][2] + beam->getElement(pos)->getY();
    542                         mat[0][3] = mat[0][3] + tan(beam->getElement(pos)->getTY()/URAD)*URAD;
    543                         xys[0] = mat.GetMatrixArray()[0]*URAD;
    544                         xys[1] = atan(mat.GetMatrixArray()[1])*URAD;
    545                         xys[2] = mat.GetMatrixArray()[2]*URAD;
    546                         xys[3] = atan(mat.GetMatrixArray()[3])*URAD;
    547                         xys[4] = beam->getElement(pos)->getS()+beam->getElement(pos)->getLength();
    548                         addPosition(xys[0],xys[1],xys[2],xys[3],xys[4]);
    549                 }
    550         }
    551 }
    552 
    553 // part about non-ip particle is not ready yet. use the above method in the meantime
     509                mat[0][0] = mat[0][0] - beam->getElement(pos)->getX();
     510                mat[0][1] = mat[0][1] - tan(beam->getElement(pos)->getTX()/URAD)*URAD;
     511                mat[0][2] = mat[0][2] - beam->getElement(pos)->getY();
     512                mat[0][3] = mat[0][3] - tan(beam->getElement(pos)->getTY()/URAD)*URAD;
     513                mat *= beam->getElement(pos)->getMatrix(energy_loss,mp,qp);
     514                mat[0][0] = mat[0][0] + beam->getElement(pos)->getX();
     515                mat[0][1] = mat[0][1] + tan(beam->getElement(pos)->getTX()/URAD)*URAD;
     516                mat[0][2] = mat[0][2] + beam->getElement(pos)->getY();
     517                mat[0][3] = mat[0][3] + tan(beam->getElement(pos)->getTY()/URAD)*URAD;
     518                xys[0] = mat.GetMatrixArray()[0]*URAD;
     519                xys[1] = atan(mat.GetMatrixArray()[1])*URAD;
     520                xys[2] = mat.GetMatrixArray()[2]*URAD;
     521                xys[3] = atan(mat.GetMatrixArray()[3])*URAD;
     522                xys[4] = beam->getElement(pos)->getS()+beam->getElement(pos)->getLength();
     523                addPosition(xys[0],xys[1],xys[2],xys[3],xys[4]);
     524                fx = xys[0];
     525                fy = xys[2];
     526                thx = xys[1];
     527                thy = xys[3];
     528        }
     529}
     530
    554531void H_BeamParticle::computePath(const H_AbstractBeamLine & beam, const bool NonLinear) {
    555532        TMatrixD temp_mat(MDIM,MDIM);
    … …  
    578555        double energy_loss = NonLinear?BE-energy:0;
    579556
    580         // modify here to allow starting at non-IP positions
    581         // s is distance to IP
    582         // initial position is already in positions vector ?
    583557        for (int i=0; i<N; i++) {
    584558                const unsigned pos = i;
    585                 // if we are inside an element, we should start by adding the action
    586                 // of the rest of this element
    587                 if(pos > 0 && fs < beam.getElement(pos)->getS() && fs > beam.getElement(pos-1)->getS()) {
    588                         cout<<"Path starts inside element "<<beam.getElement(pos-1)->getName()<<endl;
    589                         H_OpticalElement* temp_el = new H_OpticalElement(*(beam.getElement(pos-1)));
    590                         temp_el->setS(fs);
    591                         temp_el->setLength(beam.getElement(pos)->getS() - fs);
    592                         mat[0][0] = mat[0][0] - temp_el->getX();
    593                         mat[0][1] = mat[0][1] - tan(temp_el->getTX());
    594                         mat[0][2] = mat[0][2] - temp_el->getY();
    595                         mat[0][3] = mat[0][3] - tan(temp_el->getTY());
    596                         mat *= temp_el->getMatrix(energy_loss,mp,qp);
    597                         mat[0][0] = mat[0][0] + temp_el->getX();
    598                         mat[0][1] = mat[0][1] + tan(temp_el->getTX());
    599                         mat[0][2] = mat[0][2] + temp_el->getY();
    600                         mat[0][3] = mat[0][3] + tan(temp_el->getTY());
    601                 } else if(fs >= beam.getElement(pos)->getS()) {
    602                         mat[0][0] = mat[0][0] - beam.getElement(pos)->getX();
    603                         mat[0][1] = mat[0][1] - tan(beam.getElement(pos)->getTX());
    604                         mat[0][2] = mat[0][2] - beam.getElement(pos)->getY();
    605                         mat[0][3] = mat[0][3] - tan(beam.getElement(pos)->getTY());
    606                         mat *= beam.getElement(pos)->getMatrix(energy_loss,mp,qp);
    607                         mat[0][0] = mat[0][0] + beam.getElement(pos)->getX();
    608                         mat[0][1] = mat[0][1] + tan(beam.getElement(pos)->getTX());
    609                         mat[0][2] = mat[0][2] + beam.getElement(pos)->getY();
    610                         mat[0][3] = mat[0][3] + tan(beam.getElement(pos)->getTY());
    611                 }
    612                 xys[0] = mat.GetMatrixArray()[0]*URAD;
    613                 xys[1] = atan(mat.GetMatrixArray()[1])*URAD;
    614                 xys[2] = mat.GetMatrixArray()[2]*URAD;
    615                 xys[3] = atan(mat.GetMatrixArray()[3])*URAD;
    616                 xys[4] = beam.getElement(pos)->getS()+beam.getElement(pos)->getLength();
    617                 addPosition(xys[0],xys[1],xys[2],xys[3],xys[4]);
    618                 fx = xys[0];
    619                 fy = xys[2];
    620                 thx = xys[1];
    621                 thy = xys[3];
    622         }
     559                mat[0][0] = mat[0][0] - beam.getElement(pos)->getX();
     560                mat[0][1] = mat[0][1] - tan(beam.getElement(pos)->getTX());
     561                mat[0][2] = mat[0][2] - beam.getElement(pos)->getY();
     562                mat[0][3] = mat[0][3] - tan(beam.getElement(pos)->getTY());
     563                mat *= beam.getElement(pos)->getMatrix(energy_loss,mp,qp);
     564                mat[0][0] = mat[0][0] + beam.getElement(pos)->getX();
     565                mat[0][1] = mat[0][1] + tan(beam.getElement(pos)->getTX());
     566                mat[0][2] = mat[0][2] + beam.getElement(pos)->getY();
     567                mat[0][3] = mat[0][3] + tan(beam.getElement(pos)->getTY());
     568                xys[0] = mat.GetMatrixArray()[0]*URAD;
     569                xys[1] = atan(mat.GetMatrixArray()[1])*URAD;
     570                xys[2] = mat.GetMatrixArray()[2]*URAD;
     571                xys[3] = atan(mat.GetMatrixArray()[3])*URAD;
     572                xys[4] = beam.getElement(pos)->getS()+beam.getElement(pos)->getLength();
     573                addPosition(xys[0],xys[1],xys[2],xys[3],xys[4]);
     574                fx = xys[0];
     575                fy = xys[2];
     576                thx = xys[1];
     577                thy = xys[3];
     578        }
    623579}
    624580
Note: See TracChangeset for help on using the changeset viewer.