Index: trunk/Utilities/Hector/src/H_Beam.cc
===================================================================
--- trunk/Utilities/Hector/src/H_Beam.cc	(revision 233)
+++ trunk/Utilities/Hector/src/H_Beam.cc	(revision 281)
@@ -1,12 +1,19 @@
-/*
----- Hector the simulator ----
-   A fast simulator of particles through generic beamlines.
-   J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be
-
-        http://www.fynu.ucl.ac.be/hector.html
-
-   Centre de Physique des Particules et de PhÃ©nomÃ©nologie (CP3)
-   UniversitÃ© Catholique de Louvain (UCL)
-*/
+  /* * * * * * * * * * * * * * * * * * * * * * * * * * * *
+ *                                                         *
+*                   --<--<--  A fast simulator --<--<--     *
+*                 / --<--<--     of particle   --<--<--     *
+*  ----HECTOR----<                                          *
+*                 \ -->-->-- transport through -->-->--     *
+*                   -->-->-- generic beamlines -->-->--     *
+*                                                           *
+* JINST 2:P09005 (2007)                                     *
+*      X Rouby, J de Favereau, K Piotrzkowski (CP3)         *
+*       http://www.fynu.ucl.ac.be/hector.html               *
+*                                                           *
+* Center for Cosmology, Particle Physics and Phenomenology  *
+*              Universite catholique de Louvain             *
+*                 Louvain-la-Neuve, Belgium                 *
+ *                                                         *
+   * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
 /// \file H_Beam.cc
@@ -55,9 +62,5 @@
 };
 
-void H_Beam::createBeamParticles(const unsigned int Number_of_particles) {
-	createBeamParticles(Number_of_particles,MP,QP);
-}
-
-void H_Beam::createBeamParticles(const unsigned int Number_of_particles, const double p_mass, const double p_charge) {
+void H_Beam::createBeamParticles(const unsigned int Number_of_particles, const double p_mass, const double p_charge, TRandom* r) {
 	beamParticles.clear();
 	Nparticles = (Number_of_particles<1) ? 1 : Number_of_particles;
@@ -66,8 +69,8 @@
 		p.setPosition(fx_ini,fy_ini,tx_ini,ty_ini,fs_ini);
 		p.setE(fe_ini);
-		p.smearPos(x_disp,y_disp);
-		p.smearAng(tx_disp,ty_disp);
-		p.smearE(e_disp);
-		p.smearS(s_disp);
+		p.smearPos(x_disp,y_disp,r);
+		p.smearAng(tx_disp,ty_disp,r);
+		p.smearE(e_disp,r);
+		p.smearS(s_disp,r);
 		if (VERBOSE) {if (i==0) cout << " x_ini , tx_ini " << p.getX() << " " << p.getTX() << endl;}
 		beamParticles.push_back(p);
@@ -76,5 +79,5 @@
 
 //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) {
-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) {
+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) {
         beamParticles.clear();
         Nparticles = (Number_of_particles<2) ? 2 : Number_of_particles;
@@ -83,17 +86,17 @@
                 H_BeamParticle p(p_mass,p_charge);
                 if (flat) {
-                        gx = gRandom->Uniform(fx_min,fx_max);
-                        gy = gRandom->Uniform(fy_min,fy_max);
-                        gs = gRandom->Uniform(fs_min,fs_max);
-                        gtx = gRandom->Uniform(tx_min,tx_max);
-                        gty = gRandom->Uniform(ty_min,ty_max);
-                        gE = gRandom->Uniform(E_min,E_max);
+                        gx = r->Uniform(fx_min,fx_max);
+                        gy = r->Uniform(fy_min,fy_max);
+                        gs = r->Uniform(fs_min,fs_max);
+                        gtx = r->Uniform(tx_min,tx_max);
+                        gty = r->Uniform(ty_min,ty_max);
+                        gE = r->Uniform(E_min,E_max);
                 } else {
-                        gx = gRandom->Gaus((fx_min+fx_max)/2,(-fx_min+fx_max)/2);
-                        gy = gRandom->Gaus((fy_min+fy_max)/2,(-fy_min+fy_max)/2);
-                        gs = gRandom->Gaus((fs_min+fs_max)/2,(-fs_min+fs_max)/2);
-                        gtx = gRandom->Gaus((tx_min+tx_max)/2,(-tx_min+tx_max)/2);
-                        gty = gRandom->Gaus((ty_min+ty_max)/2,(-ty_min+ty_max)/2);
-                        gE = gRandom->Gaus ((E_min+E_max)/2,(-E_min+E_max)/2);
+                        gx = r->Gaus((fx_min+fx_max)/2,(-fx_min+fx_max)/2);
+                        gy = r->Gaus((fy_min+fy_max)/2,(-fy_min+fy_max)/2);
+                        gs = r->Gaus((fs_min+fs_max)/2,(-fs_min+fs_max)/2);
+                        gtx = r->Gaus((tx_min+tx_max)/2,(-tx_min+tx_max)/2);
+                        gty = r->Gaus((ty_min+ty_max)/2,(-ty_min+ty_max)/2);
+                        gE = r->Gaus ((E_min+E_max)/2,(-E_min+E_max)/2);
                 }
                 p.setPosition(gx,gy,gtx,gty,gs);
@@ -179,8 +182,4 @@
 }
 
-void H_Beam::computePath(const H_AbstractBeamLine * beamline) {
-	computePath(beamline,false);
-}
-
 /// Propagates the beam until a given s
 void H_Beam::propagate(const float position) {
@@ -189,10 +188,4 @@
                 particle_i->propagate(position);
         }
-}
-
-void H_Beam::emitGamma(const double gee, const double gq2) {
-	/// @param gee = \f$ E_{\gamma} \f$ is the photon energy
-	/// @param gq2 = \f$ Q^2 < 0 \f$ is virtuality of photon \f$ Q^{2} = E^{2}-\vec{k}^{2} \f$
-	emitGamma(gee,gq2,0,2*pi);
 }
 
@@ -386,10 +379,11 @@
 }
 
-void H_Beam::printProperties() const {
+std::ostream& operator<< (std::ostream& os, const H_Beam& be) {
 	vector<H_BeamParticle>::const_iterator particle_i;
-	cout << "There are " << Nparticles << " in the beam." << endl;
-	for (particle_i = beamParticles.begin();particle_i < beamParticles.end(); particle_i++) {
-		particle_i->printProperties();
-	}
+	cout << "There are " << be.Nparticles << " in the beam." << endl;
+	for (particle_i = be.beamParticles.begin(); particle_i < be.beamParticles.end(); particle_i++) {
+		cout << *particle_i;
+	}
+   return os;
 }
 
