Fork me on GitHub

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

switch to a more stable Hector version

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/external/Hector/H_AbstractBeamLine.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*/
     11
    1812
    1913/// \file H_AbstractBeamLine.cc
    … …  
    2721
    2822// ROOT #includes
    29 #include "TPaveLabel.h"
    30 #include "TLine.h"
    31 #include "TGaxis.h"
    32 #include "TLegend.h"
    33 #include "TF1.h"
    34 #include "TROOT.h"
     23//#include "TPaveLabel.h"
     24//#include "TLine.h"
     25//#include "TGaxis.h"
     26//#include "TLegend.h"
     27//#include "TF1.h"
     28//#include "TROOT.h"
    3529
    3630// local #includes
    … …  
    3933#include "H_Drift.h"
    4034#include "H_AbstractBeamLine.h"
     35#include "H_BeamParticle.h"
    4136#include "H_RomanPot.h"
    4237using namespace std;
    4338
    4439void H_AbstractBeamLine::init(const float length) {
    45         beam_mat.ResizeTo(MDIM,MDIM);
    46         beam_mat = driftmat(length);
     40        beam_mat = new TMatrix(MDIM,MDIM);
    4741        beam_length = length;
    4842        H_Drift * drift0 = new H_Drift("Drift0",0.,length);
    … …  
    5145}
    5246
    53 H_AbstractBeamLine::H_AbstractBeamLine(const H_AbstractBeamLine& beamline) :
    54         matrices(beamline.matrices)   {
    55         //elements = beamline.elements; //<-- bad ! the new vector contains the same pointers as the previous one
    56         cloneElements(beamline);
    57         beam_mat.ResizeTo(MDIM,MDIM);
    58         beam_mat = beamline.beam_mat;
     47H_AbstractBeamLine::H_AbstractBeamLine(const H_AbstractBeamLine& beamline) {
     48        elements = beamline.elements;
     49        matrices = beamline.matrices;
     50        beam_mat = new TMatrix(*(beamline.beam_mat));
    5951        beam_length = beamline.beam_length;
    6052}
    … …  
    6254H_AbstractBeamLine& H_AbstractBeamLine::operator=(const H_AbstractBeamLine& beamline) {
    6355        if(this== &beamline) return *this;
    64         //elements = beamline.elements; //<-- bad ! the new vector contains the same pointers as the previous one
    65         cloneElements(beamline);
     56        elements = beamline.elements;
    6657        matrices = beamline.matrices;
    67         beam_mat = beamline.beam_mat;
     58        beam_mat = new TMatrix(*(beamline.beam_mat));
    6859        beam_length = beamline.beam_length;
    6960        return *this;
    70 }
    71 
    72 H_AbstractBeamLine* H_AbstractBeamLine::clone() const {
    73         H_AbstractBeamLine* temp_beam = new H_AbstractBeamLine(beam_length);
    74         vector<H_OpticalElement*>::const_iterator element_i;
    75         for (element_i = elements.begin(); element_i<elements.end(); element_i++) {
    76                 if((*element_i)->getType()!=DRIFT) {
    77                         H_OpticalElement* temp_el = (*element_i)->clone();
    78                         temp_beam->add(temp_el);
    79                 }
    80         }
    81         temp_beam->beam_mat = beam_mat;
    82         temp_beam->matrices = matrices;
    83         return temp_beam;
    8461}
    8562
    … …  
    9168        elements.clear();
    9269        matrices.clear();
    93 }
    94 
    95 void H_AbstractBeamLine::cloneElements(const H_AbstractBeamLine& beam) {
    96         vector<H_OpticalElement*>::const_iterator element_i;
    97         for (element_i = beam.elements.begin(); element_i< beam.elements.end(); element_i++) {
    98                 H_OpticalElement* temp_el = (*element_i)->clone();
    99                 elements.push_back(temp_el);
    100         }
     70        delete beam_mat;
    10171}
    10272
    … …  
    11080        if (a > beam_length)    {
    11181                beam_length = a;
    112                 if(VERBOSE) cout<<"<H_AbstractBeamLine> WARNING : element ("<< newElement->getName()<<") too far away. The beam length has been extended to "<< beam_length << ". "<<endl;
     82                if(VERBOSE) cout<<"WARNING : element ("<< newElement->getName()<<") too far away. The beam length has been extended to "<< beam_length << ". "<<endl;
    11383        }
    11484        calcSequence();
    11585        calcMatrix();
    116         return;
    117 }
    118 
    119 const TMatrix H_AbstractBeamLine::getBeamMatrix() const {
    120         return beam_mat;
    121 }
    122 
    123 const TMatrix H_AbstractBeamLine::getBeamMatrix(const float eloss,const float p_mass, const float p_charge) {
     86}
     87
     88void H_AbstractBeamLine::add(H_OpticalElement & newElement) {
     89        /// @param newElement is added to the beamline
     90//      H_OpticalElement * el = new H_OpticalElement(newElement);
     91//      elements.push_back(el);
     92        elements.push_back(&newElement);
     93        float a = newElement.getS()+newElement.getLength();
     94        if (a > beam_length)    {
     95                beam_length = a;
     96                if(VERBOSE) cout<<"WARNING : element ("<< newElement.getName()<<") too far away. The beam length has been extended to "<< beam_length << ". "<<endl;
     97        }
     98        calcSequence();
     99        calcMatrix();
     100}
     101
     102const TMatrix * H_AbstractBeamLine::getBeamMatrix() const {
     103        TMatrix * mat = new TMatrix(*beam_mat);
     104        return mat;
     105}
     106
     107const TMatrix * H_AbstractBeamLine::getBeamMatrix(const float eloss,const float p_mass, const float p_charge) {
    124108
    125109        vector<H_OpticalElement*>::iterator element_i;
    … …  
    134118                calc_mat *= (*element_i)->getMatrix(eloss,p_mass,p_charge);
    135119        }
    136         return calc_mat;
    137 }
    138 
    139 const TMatrix H_AbstractBeamLine::getPartialMatrix(const string& elname, const float eloss, const float p_mass, const float p_charge) {
     120        const TMatrix* bmat = new TMatrix(calc_mat);
     121        return bmat;
     122}
     123
     124const TMatrix * H_AbstractBeamLine::getPartialMatrix(const string elname, const float eloss, const float p_mass, const float p_charge) {
    140125
    141126        vector<H_OpticalElement*>::iterator element_i;
    … …  
    147132                calc_mat *= (*element_i)->getMatrix(eloss,p_mass,p_charge);
    148133                if(elname==(*element_i)->getName()) {
    149                         return calc_mat;
    150                 }
    151         }
    152         cout<<"<H_AbstractBeamLine> Element "<<elname<<" desn't exist. Returning full beam matrix"<<endl;
    153         return calc_mat;
    154 }
    155 
    156 const TMatrix H_AbstractBeamLine::getPartialMatrix(const unsigned int element_position) const {
     134                        const TMatrix* bmat = new TMatrix(calc_mat);
     135                        return bmat;
     136                }
     137        }
     138        cout<<"Element "<<elname<<" desn't exist. Returning full beam matrix"<<endl;
     139        const TMatrix* bmat = new TMatrix(calc_mat);
     140        return bmat;
     141}
     142
     143const TMatrix * H_AbstractBeamLine::getPartialMatrix(const unsigned int element_position) const {
    157144        //const int N = (element_position<0)?0:(( (element_position)>elements.size()-1)?elements.size()-1:element_position);
    158145        const int N = (element_position>elements.size()-1)?elements.size()-1:element_position;
    159         return *(matrices.begin()+N); // //for optimization of the code :same as return &matrices[N];
    160 }
    161 
    162 const TMatrix H_AbstractBeamLine::getPartialMatrix(const H_OpticalElement * element) const{
     146        return &(*(matrices.begin()+N)); // //for optimization of the code :same as return &matrices[N];
     147}
     148
     149const TMatrix * H_AbstractBeamLine::getPartialMatrix(const H_OpticalElement * element) const{
    163150        // returns the transport matrix to transport until the end of the specified element
    164151        // !!! 2 elements should never have the same name in "elements" !!!
    … …  
    166153        vector<H_OpticalElement*>::const_iterator element_i;
    167154        vector<TMatrix>::const_iterator matrix_i;
    168         TMatrix calc_mat(MDIM,MDIM);
     155        TMatrix * calc_mat = new TMatrix(MDIM,MDIM);
    169156
    170157        // parses the list of optical elements and find the searched one
    … …  
    172159                if(element->getName() == (*element_i)->getName()) {
    173160                        // element has been found
    174                         calc_mat = *matrix_i;
     161                        calc_mat = const_cast<TMatrix*>( &(*matrix_i));
    175162                }
    176163        }
    177         return calc_mat;
     164        return (const TMatrix*) calc_mat;
    178165}
    179166
    … …  
    183170}
    184171
    185 H_OpticalElement * H_AbstractBeamLine::getElement(const unsigned int element_position) const {
     172const H_OpticalElement * H_AbstractBeamLine::getElement(const unsigned int element_position) const {
    186173        const unsigned int N = (element_position>elements.size())?elements.size():element_position;
    187174        return *(elements.begin()+N);//for optimization of the code :same as return &elements[N];
    … …  
    189176
    190177
    191 H_OpticalElement * H_AbstractBeamLine::getElement(const string& el_name) {
     178H_OpticalElement * H_AbstractBeamLine::getElement(const string el_name) {
    192179        for(unsigned int i=0; i < elements.size(); i++) {
    193180                if( (*(elements.begin()+i))->getName() == el_name )
    194181                        return *(elements.begin()+i);
    195182        } // if found -> return ; else : not found at all !     
    196         cout<<"<H_AbstractBeamLine> Element "<<el_name<<" not found"<<endl;
     183        cout<<"Element "<<el_name<<" not found"<<endl;
    197184        return *(elements.begin()+1);
    198185}
    199186
    200 H_OpticalElement * H_AbstractBeamLine::getElement(const string& el_name) const {
     187const H_OpticalElement * H_AbstractBeamLine::getElement(const string el_name) const {
    201188        for(unsigned int i=0; i < elements.size(); i++) {
    202189                if( (*(elements.begin()+i))->getName() == el_name)
    203190                        return *(elements.begin()+i);
    204191        } // if found -> return ; else : not found at all !
    205         cout<<"<H_AbstractBeamLine> Element "<<el_name<<" not found"<<endl;
     192        cout<<"Element "<<el_name<<" not found"<<endl;
    206193        return *(elements.begin()+1);
    207194}
    208195
    209 std::ostream& operator<< (std::ostream& os, const H_AbstractBeamLine& be) {
    210         vector<H_OpticalElement*>::const_iterator element_i;
    211         os << "Beamline content" << endl;
    212         for (element_i = be.elements.begin(); element_i < be.elements.end(); element_i++) {
    213                 os << (int)(element_i - be.elements.begin()) << "\t" << (*element_i)->getName() << "\t" << (*element_i)->getS() << endl;
    214         }
    215   return os;
    216 }
    217 
    218 void H_AbstractBeamLine::printProperties() const { cout << *this; return; }
     196void H_AbstractBeamLine::printProperties() const {
     197        vector<H_OpticalElement*>::const_iterator element_i;
     198        cout << "Pointeurs des elements du faisceau" << endl;
     199        for (element_i = elements.begin(); element_i < elements.end(); element_i++) {
     200                cout << (int)(element_i-elements.begin()) << "\t" << (*element_i)->getName() << "\t" << (*element_i)->getS() << endl;   
     201        }
     202        return;
     203}
    219204
    220205void H_AbstractBeamLine::showElements() const{
    … …  
    254239                temp = *matrix_i;
    255240                cout << "Matrix for transport until s=" << (*element_i)->getS() + (*element_i)->getLength() << "m (" << (*element_i)->getName() << "). " << endl;
    256                 printMatrix(temp);
     241                printMatrix(&temp);
    257242                cout << endl;
    258243        }
    … …  
    271256        // getting rid of drifts before calculating
    272257        for(element_i = elements.begin(); element_i < elements.end(); element_i++) {
    273                 if((*element_i)->getType() == DRIFT) {delete (*element_i); elements.erase(element_i); }
     258                if((*element_i)->getType() == DRIFT) {elements.erase(element_i); }
    274259        }
    275260
    … …  
    296281                        temp_elements.push_back(dr);
    297282        }
    298        
    299         // cleaning : avoid some memory leaks
    300283        elements.clear();
    301 
    302284        for(element_i=temp_elements.begin(); element_i < temp_elements.end(); element_i++) {
    303285                elements.push_back(*element_i);
    … …  
    321303        }
    322304
    323         beam_mat.ResizeTo(MDIM,MDIM);
    324         beam_mat = calc_mat;
     305        *beam_mat = calc_mat;
    325306        return;
    326307}
    … …  
    338319}
    339320
    340 void H_AbstractBeamLine::draw(const float xmin, const float ymin, const float xmax, const float ymax) const{
    341         gROOT->SetStyle("Plain");
    342         TLegend* leg = new TLegend(xmin,ymin,xmax,ymax,"");
     321void H_AbstractBeamLine::draw() const{
     322/*      gROOT->SetStyle("Plain");
     323        TLegend* leg = new TLegend(0.85,0.50,1,1,"Legend");
    343324        leg->SetBorderSize(1);
    344         leg->SetFillColor(kWhite);
    345325        TBox* b1 = new TBox();
    346326        TBox* b2 = new TBox(0,0,10,10);
    … …  
    365345        leg->AddEntry(b7,"RCollimator");
    366346        leg->Draw();
     347*/
    367348/*      TLine* l1 = new TLine(0.05,0.5,0.95,0.5);
    368349        TLine* l2 = new TLine(0.1,0.1,0.1,0.9);
    … …  
    413394}
    414395
    415 void H_AbstractBeamLine::drawX(const float a_min, const float a_max, const float scale) const{
     396void H_AbstractBeamLine::drawX(const float a_min, const float a_max) const{
    416397        /// @param a_min defines the size of the drawing
    417398        /// @param a_max defines the size of the drawing
    418         /// @param scale allows to multiply the drawing, i.e. changing the units
    419399        const int N = getNumberOfElements();
    420400        for(int i=0;i<N;i++) {
    … …  
    422402                float meight = fabs(a_min);
    423403                float size = (height>meight)?meight:height;
    424                 float middle = getElement(i)->getX()*URAD*scale;
     404                float middle = getElement(i)->getX()*URAD;
    425405                if(getElement(i)->getType()!=DRIFT) getElement(i)->draw(middle+size/2.,middle-size/2.);
    426406        }
    … …  
    440420}
    441421
    442 void H_AbstractBeamLine::moveElement(const string& name, const float new_s) {
     422void H_AbstractBeamLine::moveElement(const string name, const float new_s) {
    443423        /// @param name identifies the element to move
    444424        /// @param new_s is where to put it
    … …  
    453433}
    454434
    455 void H_AbstractBeamLine::alignElement(const string& name, const float disp_x, const float disp_y) {
     435void H_AbstractBeamLine::alignElement(const string name, const float disp_x, const float disp_y) {
    456436        /// @param name identifies the element to move
    457437        /// @param disp_x identifies the displacement to add in x [\f$ \mu m \f$]
    … …  
    465445                        }
    466446                }
    467                 cout<<"<H_AbstractBeamLine> WARNING : Element "<<name<<" not found."<<endl;
     447                cout<<"Element "<<name<<" not found."<<endl;
     448                if(VERBOSE) cout<<"Element "<<name<<" not found."<<endl;
    468449                return;
    469450}
    470451
    471 void H_AbstractBeamLine::tiltElement(const string& name, const float ang_x, const float ang_y) {
     452void H_AbstractBeamLine::tiltElement(const string name, const float ang_x, const float ang_y) {
    472453        /// @param name identifies the element to move
    473454        /// @param ang_x identifies the angle to add in x
    … …  
    481462                        }
    482463                }
    483                 cout<<"<H_AbstractBeamLine> WARNING : Element "<<name<<" not found."<<endl;
     464                if(VERBOSE) cout<<"Element "<<name<<" not found."<<endl;
    484465                return;
    485466}
    … …  
    500481}
    501482
     483/*
    502484TGraph * H_AbstractBeamLine::getBetaX() const{
    503485        const int N = elements.size();
    … …  
    648630        return rely;
    649631}
    650 
     632*/
Note: See TracChangeset for help on using the changeset viewer.