Fork me on GitHub

Changeset 1365 in svn for trunk/external/Hector/H_Beam.cc


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_Beam.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_Beam.cc
    … …  
    2215
    2316// ROOT #includes
    24 #include "TGraph.h"
     17//#include "TGraph.h"
    2518#include "TRandom.h"
    2619
    … …  
    6255};
    6356
    64 void H_Beam::createBeamParticles(const unsigned int Number_of_particles, const double p_mass, const double p_charge, TRandom* r) {
     57void H_Beam::createBeamParticles(const unsigned int Number_of_particles) {
     58        createBeamParticles(Number_of_particles,MP,QP);
     59}
     60
     61void H_Beam::createBeamParticles(const unsigned int Number_of_particles, const double p_mass, const double p_charge) {
    6562        beamParticles.clear();
    6663        Nparticles = (Number_of_particles<1) ? 1 : Number_of_particles;
    … …  
    6966                p.setPosition(fx_ini,fy_ini,tx_ini,ty_ini,fs_ini);
    7067                p.setE(fe_ini);
    71                 p.smearPos(x_disp,y_disp,r);
    72                 p.smearAng(tx_disp,ty_disp,r);
    73                 p.smearE(e_disp,r);
    74                 p.smearS(s_disp,r);
     68                p.smearPos(x_disp,y_disp);
     69                p.smearAng(tx_disp,ty_disp);
     70                p.smearE(e_disp);
     71                p.smearS(s_disp);
    7572                if (VERBOSE) {if (i==0) cout << " x_ini , tx_ini " << p.getX() << " " << p.getTX() << endl;}
    7673                beamParticles.push_back(p);
    7774        }
    7875}
    79 
    80 //void H_Beam::particleGun(const unsigned int Number_of_particles, const float E_min=BE, const float E_max=BE, const float fs_min=0, const float fs_max=0, const float fx_min=0, const float fx_max=0, const float fy_min=0, const float fy_max=0, const float tx_min=-PI/2., const float tx_max=PI/2., const float ty_min=-PI/2., const float ty_max=PI/2., const float p_mass=MP, const double p_charge=QP) {
    81 void H_Beam::particleGun(const unsigned int Number_of_particles, const float E_min, const float E_max, const float fs_min, const float fs_max, const float fx_min, const float fx_max, const float fy_min, const float fy_max, const float tx_min, const float tx_max, const float ty_min, const float ty_max, const float p_mass, const double p_charge, const bool flat, TRandom* r) {
    82         beamParticles.clear();
    83         Nparticles = (Number_of_particles<2) ? 2 : Number_of_particles;
    84         float gx,gy,gs,gtx,gty,gE;
    85         for (unsigned int i=0; i<Nparticles; i++) {
    86                 H_BeamParticle p(p_mass,p_charge);
    87                 if (flat) {
    88                         gx = r->Uniform(fx_min,fx_max);
    89                         gy = r->Uniform(fy_min,fy_max);
    90                         gs = r->Uniform(fs_min,fs_max);
    91                         gtx = r->Uniform(tx_min,tx_max);
    92                         gty = r->Uniform(ty_min,ty_max);
    93                         gE = r->Uniform(E_min,E_max);
    94                 } else {
    95                         gx = r->Gaus((fx_min+fx_max)/2,(-fx_min+fx_max)/2);
    96                         gy = r->Gaus((fy_min+fy_max)/2,(-fy_min+fy_max)/2);
    97                         gs = r->Gaus((fs_min+fs_max)/2,(-fs_min+fs_max)/2);
    98                         gtx = r->Gaus((tx_min+tx_max)/2,(-tx_min+tx_max)/2);
    99                         gty = r->Gaus((ty_min+ty_max)/2,(-ty_min+ty_max)/2);
    100                         gE = r->Gaus ((E_min+E_max)/2,(-E_min+E_max)/2);
    101                 }
    102                 p.setPosition(gx,gy,gtx,gty,gs);
    103                 p.setE(gE);
    104                 beamParticles.push_back(p);
    105         }
    106         return;
    107 }
    108 
    10976
    11077void H_Beam::createXScanningBeamParticles(const unsigned int Number_of_particles, const float fx_max) {
    … …  
    182149}
    183150
     151void H_Beam::computePath(const H_AbstractBeamLine * beamline) {
     152        computePath(beamline,false);
     153}
     154
    184155/// Propagates the beam until a given s
    185156void H_Beam::propagate(const float position) {
    … …  
    188159                particle_i->propagate(position);
    189160        }
     161}
     162
     163void H_Beam::emitGamma(const double gee, const double gq2) {
     164        /// @param gee = \f$ E_{\gamma} \f$ is the photon energy
     165        /// @param gq2 = \f$ Q^2 < 0 \f$ is virtuality of photon \f$ Q^{2} = E^{2}-\vec{k}^{2} \f$
     166        emitGamma(gee,gq2,0,2*PI);
    190167}
    191168
    … …  
    237214        return EY2;
    238215}
    239 
     216/*
    240217TGraphErrors * H_Beam::getBetaX(const float length, const unsigned int number_of_points) {
    241218        /// @param length [m]
    … …  
    277254        return betay;
    278255}
     256*/
    279257
    280258float H_Beam::getX(const float s, float& error_on_posx) {
    … …  
    313291}
    314292
    315 vector<TVectorD> H_Beam::getStoppingElements(const H_AbstractBeamLine * beamline, vector<H_OpticalElement>& list, vector<int>& numb) {
    316 
    317                 vector<TVectorD> stop_positions;
     293void H_Beam::getStoppingElements(const H_AbstractBeamLine * beamline, vector<H_OpticalElement>& list, vector<int>& numb) {
    318294        vector<H_BeamParticle>::iterator particle_i;
    319295        vector<H_OpticalElement>::iterator element_i;
    … …  
    331307                if(particle_i->stopped(beamline)) {
    332308                        temp_el = *(particle_i->getStoppingElement());
    333                                                 stop_positions.push_back(*(particle_i->getStopPosition()));
    334309                        if(list.size()==0) {
    335310                                number=1;
    336311                                list.push_back(temp_el);
    337312                                numb.push_back(number);
    338 
    339313                        } else {
    340314                                for (element_i = list.begin(), n_i = numb.begin(); element_i < list.end(); element_i++, n_i++) {
    … …  
    356330                } // if particle_i->stopped
    357331        }// for particle_i
    358                 return stop_positions;
    359332} // H_Beam::getStoppingElements
    360333
    … …  
    379352}
    380353
    381 std::ostream& operator<< (std::ostream& os, const H_Beam& be) {
     354void H_Beam::printProperties() const {
    382355        vector<H_BeamParticle>::const_iterator particle_i;
    383         cout << "There are " << be.Nparticles << " in the beam." << endl;
    384         for (particle_i = be.beamParticles.begin(); particle_i < be.beamParticles.end(); particle_i++) {
    385                 cout << *particle_i;
    386         }
    387    return os;
     356        cout << "There are " << Nparticles << " in the beam." << endl;
     357        for (particle_i = beamParticles.begin();particle_i < beamParticles.end(); particle_i++) {
     358                particle_i->printProperties();
     359        }
    388360}
    389361
    … …  
    400372        }
    401373} // H_Beam::printStoppingElements
    402 
    403 TH2F *  H_Beam::drawAngleProfile(const float s) {
    404         /// not a const method because does a propagate to s!
    405         char title[50];
    406         sprintf(title,"Beam profile at %.2f m",s);
    407         vector<H_BeamParticle>::iterator particle_i;
    408         float xmax, xmin, ymax, ymin;
    409         float xx, yy, xborder, yborder;
    410 
    411         particle_i=beamParticles.begin();
    412         xmin = particle_i->getTX();
    413         xmax = particle_i->getTX();
    414         ymin = particle_i->getTY();
    415         ymax = particle_i->getTY();
    416 
    417         for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
    418                 particle_i->propagate(s);
    419                 xx = particle_i->getTX();
    420                 yy = particle_i->getTY();
    421 
    422                 xmax = xx>xmax ? xx : xmax;
    423                 ymax = yy>ymax ? yy : ymax;
    424                 xmin = xx<xmin ? xx : xmin;
    425                 ymin = yy<ymin ? yy : ymin;
    426         }
    427 
    428         // in order to avoid some drawing problems, when the beam divergence is null
    429         if(!(xmax || xmin)) xmax +=0.1;
    430         if(!(ymax || ymin)) xmax +=0.1;
    431 
    432         if(xmax == xmin) xmax *= 1.1;
    433         if(ymax == ymin) ymax *= 1.1;
    434 
    435         xborder = (xmax-xmin)*0.2;
    436         yborder = (ymax-ymin)*0.2;
    437 
    438         xmax += xborder;
    439         xmin -= xborder;
    440         ymax += yborder;
    441         ymin -= yborder;
    442 
    443         TH2F * profile = new TH2F("profile",title,10000,xmin,xmax,1000,ymin,ymax);
    444         for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
    445                 profile->Fill(particle_i->getTX(), particle_i->getTY());
    446         }
    447         return profile;
    448 }
    449 
    450 
     374/*
    451375TH2F *  H_Beam::drawProfile(const float s) {
    452376        /// not a const method because does a propagate to s!
    … …  
    475399
    476400        // in order to avoid some drawing problems, when the beam divergence is null
    477         if(!(xmax || xmin)) xmax +=0.1;
    478         if(!(ymax || ymin)) xmax +=0.1;
    479 
    480401        if(xmax == xmin) xmax += 0.1;
    481402        if(ymax == ymin) ymax += 0.1;
    … …  
    494415        }
    495416        return profile;
    496 }
    497 
     417}*/
     418/*
    498419TMultiGraph * H_Beam::drawBeamX(const int color) const {
    499420        int mycolor = color;
    … …  
    519440        return beam_profile_y;
    520441}
     442*/
Note: See TracChangeset for help on using the changeset viewer.