Fork me on GitHub

Ignore:
Timestamp:
Mar 1, 2009, 3:56:11 PM (15 years ago)
Author:
Xavier Rouby
Message:

new Hector version

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/Utilities/Hector/src/H_Beam.cc

    r233 r281  
    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 */
     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   * * * * * * * * * * * * * * * * * * * * * * * * * * * */
    1118
    1219/// \file H_Beam.cc
    … …  
    5562};
    5663
    57 void H_Beam::createBeamParticles(const unsigned int Number_of_particles) {
    58         createBeamParticles(Number_of_particles,MP,QP);
    59 }
    60 
    61 void H_Beam::createBeamParticles(const unsigned int Number_of_particles, const double p_mass, const double p_charge) {
     64void H_Beam::createBeamParticles(const unsigned int Number_of_particles, const double p_mass, const double p_charge, TRandom* r) {
    6265        beamParticles.clear();
    6366        Nparticles = (Number_of_particles<1) ? 1 : Number_of_particles;
    … …  
    6669                p.setPosition(fx_ini,fy_ini,tx_ini,ty_ini,fs_ini);
    6770                p.setE(fe_ini);
    68                 p.smearPos(x_disp,y_disp);
    69                 p.smearAng(tx_disp,ty_disp);
    70                 p.smearE(e_disp);
    71                 p.smearS(s_disp);
     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);
    7275                if (VERBOSE) {if (i==0) cout << " x_ini , tx_ini " << p.getX() << " " << p.getTX() << endl;}
    7376                beamParticles.push_back(p);
    … …  
    7679
    7780//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) {
    78 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) {
     81void 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) {
    7982        beamParticles.clear();
    8083        Nparticles = (Number_of_particles<2) ? 2 : Number_of_particles;
    … …  
    8386                H_BeamParticle p(p_mass,p_charge);
    8487                if (flat) {
    85                         gx = gRandom->Uniform(fx_min,fx_max);
    86                         gy = gRandom->Uniform(fy_min,fy_max);
    87                         gs = gRandom->Uniform(fs_min,fs_max);
    88                         gtx = gRandom->Uniform(tx_min,tx_max);
    89                         gty = gRandom->Uniform(ty_min,ty_max);
    90                         gE = gRandom->Uniform(E_min,E_max);
     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);
    9194                } else {
    92                         gx = gRandom->Gaus((fx_min+fx_max)/2,(-fx_min+fx_max)/2);
    93                         gy = gRandom->Gaus((fy_min+fy_max)/2,(-fy_min+fy_max)/2);
    94                         gs = gRandom->Gaus((fs_min+fs_max)/2,(-fs_min+fs_max)/2);
    95                         gtx = gRandom->Gaus((tx_min+tx_max)/2,(-tx_min+tx_max)/2);
    96                         gty = gRandom->Gaus((ty_min+ty_max)/2,(-ty_min+ty_max)/2);
    97                         gE = gRandom->Gaus ((E_min+E_max)/2,(-E_min+E_max)/2);
     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);
    98101                }
    99102                p.setPosition(gx,gy,gtx,gty,gs);
    … …  
    179182}
    180183
    181 void H_Beam::computePath(const H_AbstractBeamLine * beamline) {
    182         computePath(beamline,false);
    183 }
    184 
    185184/// Propagates the beam until a given s
    186185void H_Beam::propagate(const float position) {
    … …  
    189188                particle_i->propagate(position);
    190189        }
    191 }
    192 
    193 void H_Beam::emitGamma(const double gee, const double gq2) {
    194         /// @param gee = \f$ E_{\gamma} \f$ is the photon energy
    195         /// @param gq2 = \f$ Q^2 < 0 \f$ is virtuality of photon \f$ Q^{2} = E^{2}-\vec{k}^{2} \f$
    196         emitGamma(gee,gq2,0,2*pi);
    197190}
    198191
    … …  
    386379}
    387380
    388 void H_Beam::printProperties() const {
     381std::ostream& operator<< (std::ostream& os, const H_Beam& be) {
    389382        vector<H_BeamParticle>::const_iterator particle_i;
    390         cout << "There are " << Nparticles << " in the beam." << endl;
    391         for (particle_i = beamParticles.begin();particle_i < beamParticles.end(); particle_i++) {
    392                 particle_i->printProperties();
    393         }
     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;
    394388}
    395389
Note: See TracChangeset for help on using the changeset viewer.