Fork me on GitHub

Changeset 3c40083 in git for external/Hector/H_BeamLine.cc


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_BeamLine.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_BeamLine.cc
    … …  
    2417#include <fstream>
    2518#include <sstream>
    26 #include <cstdlib>
    2719
    2820// local #includes
    … …  
    4739// are opposite to Wille's. See fill() for more details.
    4840
    49 H_BeamLine::H_BeamLine(): H_AbstractBeamLine(),
    50         direction(1), ips(0), ipx(0), ipy(0), iptx(0), ipty(0) {
    51 }
    52 
    53 H_BeamLine::H_BeamLine(const int si, const float length) : H_AbstractBeamLine(length),
    54         direction((si >= abs(si)) ? 1 : -1), ips(0), ipx(0), ipy(0), iptx(0), ipty(0) {
    55 }
    56 
    57 H_BeamLine::H_BeamLine(const H_BeamLine& beam) : H_AbstractBeamLine(beam),
    58         direction(beam.direction), ips(beam.ips), ipx(beam.ipx), ipy(beam.ipy), iptx(beam.iptx), ipty(beam.ipty) {
    59 }
    60 
    61 H_BeamLine& H_BeamLine::operator=(const H_BeamLine& beam) {
    62         if(this==&beam) return *this;
    63         H_AbstractBeamLine::operator=(beam); // call the mother's operator=
     41H_BeamLine::H_BeamLine(const int si, const float length) : H_AbstractBeamLine(length){
     42        direction = (si >= abs(si)) ? 1 : -1;
     43        ips=0;
     44        ipx=0;
     45        ipy=0;
     46        iptx=0;
     47        ipty=0;
     48}
     49
     50H_BeamLine::H_BeamLine(const H_BeamLine& beam) : H_AbstractBeamLine(beam) {
    6451        direction = beam.direction;
    6552        ips = beam.ips;
    … …  
    6855        iptx = beam.iptx;
    6956        ipty = beam.ipty;
     57}
     58
     59H_BeamLine& H_BeamLine::operator=(const H_BeamLine& beam) {
     60        if(this==&beam) return *this;
     61        direction = beam.direction;
     62        ips = beam.ips;
     63        ipx = beam.ipx;
     64    ipy = beam.ipy;
     65    iptx = beam.iptx;
     66    ipty = beam.ipty;
    7067        return *this;
    7168}
    7269
    73 void H_BeamLine::findIP(const string& filename, const string& ipname) {
     70void H_BeamLine::findIP(const string filename) {
     71        findIP(filename,"IP5");
     72        return;
     73}
     74
     75void H_BeamLine::findIP(const string filename, const string ipname) {
    7476        // searches for the IP position in the extended table.
    7577        ifstream tabfile(filename.c_str());
    76                 if (! tabfile.is_open()) cout << "<H_BeamLine> ERROR: I Can't open \"" << filename << "\"" << endl;
     78                if (! tabfile.is_open()) cout << "\t ERROR: I Can't open \"" << filename << "\"" << endl;
    7779        bool found = false;
    7880        int N_col=0;
    … …  
    112114        } // while (!eof)
    113115
    114         if (!found) cout << "<H_BeamLine> ERROR ! IP not found." << endl;
     116        if (!found) cout << "\t ERROR ! IP not found." << endl;
    115117        tabfile.close();
    116118}
    … …  
    125127}
    126128
    127 
    128 void H_BeamLine::fill(const string& filename, const int dir, const string& ipname) {
     129void H_BeamLine::fill(const string filename) {
     130        fill(filename,1,"IP5");
     131        return;
     132}
     133
     134
     135void H_BeamLine::fill(const string filename, const int dir, const string ipname) {
    129136        string headers[40];  // table of all column headers
    130137        int header_type[40]; // table of all column types
    131138        findIP(filename,ipname);
    132139        ifstream tabfile(filename.c_str());
    133                 if (! tabfile.is_open()) cout << "<H_BeamLine> ERROR: I Can't open \"" << filename << "\"" << endl;
     140                if (! tabfile.is_open()) cout << "\t ERROR: I Can't open \"" << filename << "\"" << endl;
    134141        if(VERBOSE) {
    135142                cout<<"Using file : "<< filename <<" in the "<<((direction>0)?"positive":"negative")<<" direction."<<endl;
    … …  
    187194                        }
    188195
     196                        ap = 0; //init
     197                        if(e.aper_1 !=0 || e.aper_2 !=0 || e.aper_3 !=0 || e.aper_4 != 0) {
     198                                e.aper_1 *= URAD;
     199                                e.aper_2 *= URAD;
     200                                e.aper_3 *= URAD;
     201                                e.aper_4 *= URAD; // in [m] in the tables !
     202
     203                                if(strstr(e.apertype.c_str(),"RECTELLIPSE")) ap = new H_RectEllipticAperture(e.aper_1,e.aper_2,e.aper_3,e.aper_4,0,0);
     204                                else if(strstr(e.apertype.c_str(),"CIRCLE")) ap = new H_CircularAperture(e.aper_1,0,0);
     205                                else if(strstr(e.apertype.c_str(),"RECTANGLE")) ap = new H_RectangularAperture(e.aper_1,e.aper_2,0,0);
     206                                else if(strstr(e.apertype.c_str(),"ELLIPSE")) ap = new H_EllipticAperture(e.aper_1,e.aper_2,0,0);
     207                        }
     208
    189209                        if(el!=0) {
    190                                 ap = 0; //init
    191                                 if(e.aper_1 !=0 || e.aper_2 !=0 || e.aper_3 !=0 || e.aper_4 != 0) {
    192                                         e.aper_1 *= URAD;
    193                                         e.aper_2 *= URAD;
    194                                         e.aper_3 *= URAD;
    195                                         e.aper_4 *= URAD; // in [m] in the tables !
    196 
    197                                         if(strstr(e.apertype.c_str(),"RECTELLIPSE"))
    198                                                 ap = new H_RectEllipticAperture(e.aper_1,e.aper_2,e.aper_3,e.aper_4,0,0);
    199                                         else if(strstr(e.apertype.c_str(),"CIRCLE"))
    200                                                 ap = new H_CircularAperture(e.aper_1,0,0);
    201                                         else if(strstr(e.apertype.c_str(),"RECTANGLE"))
    202                                                 ap = new H_RectangularAperture(e.aper_1,e.aper_2,0,0);
    203                                         else if(strstr(e.apertype.c_str(),"ELLIPSE"))
    204                                                 ap = new H_EllipticAperture(e.aper_1,e.aper_2,0,0);
     210                                if (direction<0) {
     211                                        el->setBetaX(e.betx);
     212                                        el->setBetaY(e.bety);
     213                                        el->setDX(e.dx);
     214                                        el->setDY(e.dy);
     215                                        el->setRelX(e.x);
     216                                        el->setRelY(e.y);
     217                                } else {
     218                                        el->setBetaX(previous_betax);
     219                                        el->setBetaY(previous_betay);
     220                                        el->setDX(previous_dx);
     221                                        el->setDY(previous_dy);
     222                                        el->setRelX(previous_x);
     223                                        el->setRelY(previous_y);
     224                                }
     225                                if(ap!=0) {
     226                                        el->setAperture(ap);
     227                                //              delete ap; // ap deleted in H_AbstractBeamLine::~H_AbstractBeamLine
    205228                                }
    206229       
    207                                 if (direction<0) {
    208                                         el->setBetaX(e.betx);   el->setBetaY(e.bety);
    209                                         el->setDX(e.dx);        el->setDY(e.dy);
    210                                         el->setRelX(e.x);       el->setRelY(e.y);
    211                                 } else {
    212                                         el->setBetaX(previous_betax);   el->setBetaY(previous_betay);
    213                                         el->setDX(previous_dx);         el->setDY(previous_dy);
    214                                         el->setRelX(previous_x);        el->setRelY(previous_y);
    215                                 }
    216                                 if(ap!=0) {
    217                                         el->setAperture(ap);
    218                                         delete ap;
    219                                         // if memory is allocated, i.e. if ap!=0,
    220                                         // it wont be deallocated in H_OpticalElement::~H_OpticalElement
    221                                         // ap should then be deleted here
    222                                 } // memory leak-free!
    223        
    224                                 /// Parses all the elements from the beginning of the file,
    225                                 /// but only keeps the ones from the IP till the desired length
    226                                 if(e.s>=0 && e.s<beam_length) add(el);
    227                                 else { delete el;}
    228                                 // NB : if "el" is added to the beamline, it will be
    229                                 // deleted by H_AbstractBeamLine::~H_AbstractBeamLine
    230                                 // Otherwise, it should be deleted explicitly
    231 
     230                                /// Parses all the elements, but only keeps the ones from the IP till the desired length
     231                                if(e.s>=0 && e.s<beam_length) add(*el);
     232
     233                        // delete el; // el deleted in H_AbstractBeamLine::~H_AbstractBeamLine
    232234                        } // if(el!=0)
    233235                } // if (found)
Note: See TracChangeset for help on using the changeset viewer.