Fork me on GitHub

Ignore:
Timestamp:
Apr 16, 2014, 3:56:14 PM (11 years ago)
Author:
pavel <pavel@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
64a4950
Parents:
f6b9fec
Message:

switch to a more stable Hector version

File:
1 edited

Legend:

Unmodified
Added
Removed

  • external/Hector/H_OpticalElement.cc

    rf6b9fec r3c40083  
    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_OpticalElement.cc
    … …  
    2518
    2619// ROOT #includes
    27 #include "TPaveLabel.h"
     20//#include "TPaveLabel.h"
    2821
    2922// local #includes
    30 #include "H_Parameters.h"
    3123#include "H_TransportMatrices.h"
    3224#include "H_Aperture.h"
    … …  
    3527
    3628/// called by the constructors
    37 void H_OpticalElement::init(const string& nameE, const int typeE, const double s, const double k, const double l) {
     29void H_OpticalElement::init(const string nameE, const int typeE, const double s, const double k, const double l, H_Aperture* the_app) {
    3830        // this is called by the constructors
    3931        // must be in public section !
    … …  
    4840        element_length = l;
    4941        type = typeE;
    50         element_mat.ResizeTo(MDIM,MDIM);
    51         element_mat = driftmat(l);
    52 
    53         // do NOT use setAperture for the initialisation ! there are protections there
     42        element_mat = new TMatrix(MDIM,MDIM);
     43        setAperture(the_app);
    5444       
    55         if(element_length<0)  { if(VERBOSE) cout<<"<H_OpticalElement> ERROR : Interpenetration of elements !"<<endl; }
    56         if(element_length==0) { if(VERBOSE) cout<<"<H_OpticalElement> WARNING : 0-length element ! (" << name << ") " << " at " << fs << endl; }
     45        if(element_length<0)  { if(VERBOSE) cout<<"\t ERROR : Interpenetration of elements !"<<endl; }
     46        if(element_length==0) { if(VERBOSE) cout<<"\t WARNING : 0-length element ! (" << name << ") " << " at " << fs << endl; }
    5747        betax =0;
    5848        betay =0;
    5949}
    6050
    61 H_OpticalElement::H_OpticalElement() : element_aperture(new H_Aperture()) {
    62         init("",DRIFT,0.,0.,0.1);
     51H_OpticalElement::H_OpticalElement() {
     52        H_Aperture* the_app = new H_Aperture();
     53        init("",DRIFT,0.,0.,0.1,the_app);
    6354}
    6455
    65 H_OpticalElement::H_OpticalElement(const string& nameE, const int typeE, const double s, const double k, const double l) : element_aperture(new H_Aperture()) {
    66         init(nameE,typeE,s,k,l);
     56H_OpticalElement::H_OpticalElement(const string nameE, const int typeE, const double s, const double k, const double l) {
     57        H_Aperture* the_app = new H_Aperture();
     58        init(nameE,typeE,s,k,l,the_app);
    6759}
    6860
    69 H_OpticalElement::H_OpticalElement(const string& nameE, const int typeE, const double s, const double k, const double l, H_Aperture* the_app) : element_aperture(the_app->clone()) {
    70         init(nameE,typeE,s,k,l);
     61H_OpticalElement::H_OpticalElement(const string nameE, const int typeE, const double s, const double k, const double l, H_Aperture* the_app) {
     62        init(nameE,typeE,s,k,l,the_app);
    7163}
    7264
    73 H_OpticalElement::H_OpticalElement(const int typeE, const double s, const double k, const double l, H_Aperture* the_app) : element_aperture(the_app->clone()) {
    74         init("",typeE,s,k,l);
     65H_OpticalElement::H_OpticalElement(const int typeE, const double s, const double k, const double l, H_Aperture* the_app) {
     66        init("",typeE,s,k,l,the_app);
    7567}
    7668
    77 H_OpticalElement::H_OpticalElement(const int typeE, const double s, const double k, const double l) : element_aperture(new H_Aperture()) {
    78         init("",typeE,s,k,l);
     69H_OpticalElement::H_OpticalElement(const int typeE, const double s, const double k, const double l) {
     70        H_Aperture* the_app = new H_Aperture();
     71        init("",typeE,s,k,l,the_app);
    7972}
    8073
    … …  
    9285        name = el.name;
    9386        typestring = el.typestring;
    94         element_mat.ResizeTo(MDIM,MDIM);
    95         element_mat = el.element_mat;
    96         element_aperture = el.element_aperture->clone();
     87        element_mat = new TMatrix(*(el.element_mat));
     88        element_aperture = new H_Aperture(*(el.element_aperture));
    9789}
    9890
    … …  
    10496        xpos = el.xpos;
    10597        ypos = el.ypos;
    106         txpos = el.txpos;
    107         typos = el.typos;
     98                txpos = el.txpos;
     99                typos = el.typos;
    108100        betax = el.betax;
    109101        betay = el.betay;
    … …  
    111103        name = el.name;
    112104        typestring = el.typestring;
    113         element_mat.ResizeTo(MDIM,MDIM);
    114         element_mat = el.element_mat;
    115         element_aperture = el.element_aperture->clone();
     105        delete element_mat;
     106        delete element_aperture;
     107        element_mat = new TMatrix(*(el.element_mat));
     108        element_aperture = new H_Aperture(*(el.element_aperture));
    116109        return *this;
    117110}
    118111
    119 void H_OpticalElement::setAperture(const H_Aperture* ap) {
    120         // do NOT use setAperture in your constructor, as element_aperture is not initialized
    121         // this function do not take into account ap if ap=0
    122         // do nothing if element_mat = ap
    123                 if (!ap) {cout << "<H_OpticalElement> Trying to set an empty pointer for the aperture ! Nothing done.\n"; return;}
    124                 if (element_aperture != ap) {
    125                         delete element_aperture;
    126                         element_aperture = ap->clone();
    127                 }
     112void H_OpticalElement::setAperture(H_Aperture * ap) {
     113//      element_aperture = const_cast<H_Aperture*>(ap);
     114        element_aperture = ap;
    128115        return;
    129116}
    130 /*
    131 void H_OpticalElement::setAperture(H_Aperture* ap) {
    132         // do NOT use setAperture in your constructor, as element_aperture is not initialized
    133         // this function do not take into account ap if ap=0
    134         // do nothing if element_mat = ap
    135                 if (!ap) {cout << "<H_OpticalElement> Trying to set an empty pointer for the aperture ! Nothing done.\n"; return;}
    136                 if (element_aperture != ap) {
    137                         delete element_aperture;
    138                         element_aperture = ap; //->clone();
    139                 }
    140         return;
    141 }
    142 */
    143117
    144 std::ostream& operator<< (std::ostream& os, const H_OpticalElement& el) {
    145         os << el.typestring << el.name << "\t at s = " << el.fs << "\t length = "<< el.element_length;
    146         if(el.fk!=0) os <<"\t strength = " << el.fk;
    147         if(el.element_aperture->getType()!=NONE) {
    148                 os << *(el.element_aperture) << endl;
    149         }
    150         os<<endl;
    151         if(el.element_length<0 && VERBOSE)
    152                 os <<"<H_OpticalElement> ERROR : Interpenetration of elements !"<<endl;
    153         else if(el.element_length==0 && VERBOSE)
    154                 os <<"<H_OpticalElement> WARNING : 0-length "<< el.typestring << " !" << endl;
    155   return os;
     118void H_OpticalElement::printProperties() const {
     119                cout << typestring;
     120                cout << name;
     121                cout <<"\t at s = " << fs;
     122                cout <<"\t length = "<< element_length;
     123                if(fk!=0) cout <<"\t strength = " << fk;
     124                if(element_aperture->getType()!=NONE) {
     125                        cout <<"\t aperture type = " << element_aperture->getTypeString();
     126                        element_aperture->printProperties();                   
     127                }
     128       
     129                cout<<endl;
     130                if(element_length<0)  { if(VERBOSE) cout<<"\t ERROR : Interpenetration of elements !"<<endl; }
     131                if(element_length==0) { if(VERBOSE) cout<<"\t WARNING : 0-length "<< typestring << " !" << endl; }
     132
     133 return;
    156134}
    157135
    158136void H_OpticalElement::showMatrix() const {
    159         printMatrix(element_mat);
     137                printMatrix(element_mat);
    160138        return;
    161139}
    162140
    163141void H_OpticalElement::drawAperture() const {
    164         element_aperture->draw(1);
     142        element_aperture->draw();
    165143        return;
    166144}
    167145
    168 TMatrix H_OpticalElement::getMatrix() {
    169         setMatrix(0,MP,1);
    170         return element_mat;
     146TMatrix H_OpticalElement::getMatrix() const {
     147        return *element_mat;
    171148}
    172149
    173 TMatrix H_OpticalElement::getMatrix(const float eloss, const float p_mass, const float p_charge) {
     150TMatrix H_OpticalElement::getMatrix(const float eloss, const float p_mass, const float p_charge) const {
    174151        setMatrix(eloss,p_mass,p_charge);
    175         return element_mat;
     152        return *element_mat;
    176153}
    177154
    178155void H_OpticalElement::draw(const float meight, const float height) const{
    179         /// @param meight is the minimal extend of the graph
     156/*      /// @param meight is the minimal extend of the graph
    180157        /// @param height is the maximal extend of the graph
    181158        float x1 = getS();
    … …  
    188165        cur_box->SetFillColor((int)getType());
    189166        cur_box->Draw();
     167*/
    190168}
    191169
    192 TVectorD H_OpticalElement::getHitPosition(const TVectorD& init_pos, const double energy_loss, const double mp, const double qp) {
    193         if(!element_length) {
    194                 // cout<<"O-length element ("<<getName()<<"), should not appear here !"<<endl;
    195                 return init_pos;
    196         }
    197         // some declarations
    198         bool inside = false;
    199         double vec[MDIM] = {init_pos[INDEX_X]/URAD,  tan(init_pos[INDEX_TX]/URAD),
    200                             init_pos[INDEX_Y]/URAD,  tan(init_pos[INDEX_TY]/URAD),
    201                             -energy_loss,            1};
    202         TMatrixD mat_init(1,MDIM,vec);
    203         TMatrixD mat_min(1,MDIM,vec);
    204         TMatrixD mat_max(1,MDIM,vec);
    205         TMatrixD mat_stop(1,MDIM,vec);
    206         H_OpticalElement* temp_el = clone();
    207         // initialasing boundaries
    208         double min_pos = 0;
    209         double max_pos = element_length/2.;
    210         double max_old = element_length/2.;
    211         // number of iterations
    212         // (idea : fix precision instead of number of iterations + add security)
    213         // (idea : interpolate between max and min and give the error)
    214         const int N = 10;
    215         // starting search loop
    216         for(int i = 0; i < N; i++) {
    217                 // fixing position to be investigated
    218                 temp_el->setLength(max_pos);
    219                 // initialising the vector at the initial vector + possible shift/tilt
    220                 mat_max[0][0] = mat_init[0][0] - temp_el->getX();
    221                 mat_max[0][1] = mat_init[0][1] - tan(temp_el->getTX());
    222                 mat_max[0][2] = mat_init[0][2] - temp_el->getY();
    223                 mat_max[0][3] = mat_init[0][3] - tan(temp_el->getTY());
    224                 // propagating
    225                 mat_max *= temp_el->getMatrix(energy_loss,mp,qp);
    226                 // compensating for the previously stated shifts/tilts
    227                 mat_max[0][0] = mat_max[0][0] + temp_el->getX();
    228                 mat_max[0][1] = mat_max[0][1] + tan(temp_el->getTX());
    229                 mat_max[0][2] = mat_max[0][2] + temp_el->getY();
    230                 mat_max[0][3] = mat_max[0][3] + tan(temp_el->getTY());
    231                 // fixing new boundaries
    232                 if(temp_el->isInside(mat_max.GetMatrixArray()[0]*URAD,mat_max.GetMatrixArray()[2]*URAD)) {
    233                         max_old = max_pos;
    234                         max_pos = max_pos + (max_pos - min_pos)/2.;
    235                         min_pos = max_old;
    236                         inside = true;
    237                 } else {
    238                         max_pos = min_pos + (max_pos - min_pos)/2.;
    239                         inside = false;
    240                 }
    241                 // end of loop
    242         }
    243         // if it passes at the last iteration, choosing the other range²
    244         if(inside) min_pos = max_old;
    245         // here the interpolation method : now we are sure that the intercept is between min_pos and max_pos
    246         // getting vector at min_pos (for first boundary) :
    247         bool precision_estimate = false;
    248         if(precision_estimate) {
    249                 temp_el->setLength(min_pos);
    250                 mat_min[0][0] = mat_init[0][0] - temp_el->getX();
    251                 mat_min[0][1] = mat_init[0][1] - tan(temp_el->getTX());
    252                 mat_min[0][2] = mat_init[0][2] - temp_el->getY();
    253                 mat_min[0][3] = mat_init[0][3] - tan(temp_el->getTY());
    254                 mat_min *= temp_el->getMatrix(energy_loss,mp,qp);
    255                 mat_min[0][0] = mat_min[0][0] + temp_el->getX();
    256                 mat_min[0][1] = mat_min[0][1] + tan(temp_el->getTX());
    257                 mat_min[0][2] = mat_min[0][2] + temp_el->getY();
    258                 mat_min[0][3] = mat_min[0][3] + tan(temp_el->getTY());
    259                 mat_min[0][4] = min_pos + init_pos[4];
    260                 // getting vector at max_pos (for second boundary) :
    261                 temp_el->setLength(max_pos);
    262                 mat_max[0][0] = mat_init[0][0] - temp_el->getX();
    263                 mat_max[0][1] = mat_init[0][1] - tan(temp_el->getTX());
    264                 mat_max[0][2] = mat_init[0][2] - temp_el->getY();
    265                 mat_max[0][3] = mat_init[0][3] - tan(temp_el->getTY());
    266                 mat_max *= temp_el->getMatrix(energy_loss,mp,qp);
    267                 mat_max[0][0] = mat_max[0][0] + temp_el->getX();
    268                 mat_max[0][1] = mat_max[0][1] + tan(temp_el->getTX());
    269                 mat_max[0][2] = mat_max[0][2] + temp_el->getY();
    270                 mat_max[0][3] = mat_max[0][3] + tan(temp_el->getTY());
    271                 mat_max[0][4] = max_pos + init_pos[4];
    272         }
    273         // getting vector in the middle (for estimate) :
    274         temp_el->setLength((max_pos+min_pos)/2.);
    275         mat_stop[0][0] = mat_init[0][0] - temp_el->getX();
    276         mat_stop[0][1] = mat_init[0][1] - tan(temp_el->getTX());
    277         mat_stop[0][2] = mat_init[0][2] - temp_el->getY();
    278         mat_stop[0][3] = mat_init[0][3] - tan(temp_el->getTY());
    279         mat_stop *= temp_el->getMatrix(energy_loss,mp,qp);
    280         mat_stop[0][0] = mat_stop[0][0] + temp_el->getX();
    281         mat_stop[0][1] = mat_stop[0][1] + tan(temp_el->getTX());
    282         mat_stop[0][2] = mat_stop[0][2] + temp_el->getY();
    283         mat_stop[0][3] = mat_stop[0][3] + tan(temp_el->getTY());
    284         mat_stop[0][4] = (max_pos+min_pos)/2. + init_pos[4];
    285 
    286         double xys[LENGTH_VEC];
    287         xys[INDEX_X]=  mat_stop[0][0]*URAD;
    288         xys[INDEX_TX]= atan(mat_stop[0][1])*URAD;
    289         xys[INDEX_Y]=  mat_stop[0][2]*URAD;
    290         xys[INDEX_TY]= atan(mat_stop[0][3])*URAD;
    291         xys[INDEX_S]=  mat_stop[0][4] ;
    292         TVectorD temp_vec(LENGTH_VEC,xys);
    293 
    294         if(precision_estimate) {
    295                 cout<<"--- Results and precision estimates ---"<<endl;
    296                 cout<<"\t Stopping element : "<<getName()<<endl;
    297                 cout<<"\t hit point s  : "<<mat_stop[0][4]<<" m +- "<<(mat_max[0][4]-mat_stop[0][4])*1000.<<" mm"<<endl;
    298                 cout<<"\t hit point x  : "<<mat_stop[0][0]*URAD;
    299                 cout<<" + "<<fabs(((mat_min[0][0]<mat_max[0][0])?(mat_max[0][0]-mat_stop[0][0])*URAD:(mat_min[0][0]-mat_stop[0][0])*URAD));
    300                 cout<<" - "<<fabs(((mat_min[0][0]<mat_max[0][0])?(mat_stop[0][0]-mat_min[0][0])*URAD:(mat_stop[0][0]-mat_max[0][0])*URAD))<<" µm"<<endl;
    301                 cout<<"\t hit point y  : "<<mat_stop[0][2]*URAD;
    302                 cout<<" + "<<fabs(((mat_min[0][0]<mat_max[0][2])?(mat_max[0][2]-mat_stop[0][2])*URAD:(mat_min[0][2]-mat_stop[0][2])*URAD));
    303                 cout<<" - "<<fabs(((mat_min[0][0]<mat_max[0][2])?(mat_stop[0][2]-mat_min[0][2])*URAD:(mat_stop[0][2]-mat_max[0][2])*URAD))<<" µm"<<endl;
    304                 cout<<"\t hit point tx : "<<mat_stop[0][1]*URAD;
    305                 cout<<" + "<<fabs(((mat_min[0][0]<mat_max[0][1])?(mat_max[0][1]-mat_stop[0][1])*URAD:(mat_min[0][1]-mat_stop[0][1])*URAD));
    306                 cout<<" - "<<fabs(((mat_min[0][0]<mat_max[0][1])?(mat_stop[0][1]-mat_min[0][1])*URAD:(mat_stop[0][1]-mat_max[0][1])*URAD))<<" µrad"<<endl;
    307                 cout<<"\t hit point ty : "<<mat_stop[0][3]*URAD;
    308                 cout<<" + "<<fabs(((mat_min[0][0]<mat_max[0][3])?(mat_max[0][3]-mat_stop[0][3])*URAD:(mat_min[0][3]-mat_stop[0][3])*URAD));
    309                 cout<<" - "<<fabs(((mat_min[0][0]<mat_max[0][3])?(mat_stop[0][3]-mat_min[0][3])*URAD:(mat_stop[0][3]-mat_max[0][3])*URAD))<<" µrad"<<endl;
    310 
    311         }
    312 
    313         delete temp_el;
    314         return temp_vec;
    315 }
Note: See TracChangeset for help on using the changeset viewer.