Index: external/Hector/H_Beam.cc
===================================================================
--- external/Hector/H_Beam.cc	(revision 5b822e567e81b047f743b282d78eaae9b7c4fb20)
+++ external/Hector/H_Beam.cc	(revision 3c40083cb5c3b145665840d29f1189d1de1c3a0d)
@@ -1,19 +1,12 @@
-  /* * * * * * * * * * * * * * * * * * * * * * * * * * * *
- *                                                         *
-*                   --<--<--  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                 *
- *                                                         *
-   * * * * * * * * * * * * * * * * * * * * * * * * * * * */
+/*
+---- 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)
+*/
 
 /// \file H_Beam.cc
@@ -22,5 +15,5 @@
 
 // ROOT #includes
-#include "TGraph.h"
+//#include "TGraph.h"
 #include "TRandom.h"
 
@@ -62,5 +55,9 @@
 };
 
-void H_Beam::createBeamParticles(const unsigned int Number_of_particles, const double p_mass, const double p_charge, TRandom* r) {
+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) {
 	beamParticles.clear();
 	Nparticles = (Number_of_particles<1) ? 1 : Number_of_particles;
@@ -69,42 +66,12 @@
 		p.setPosition(fx_ini,fy_ini,tx_ini,ty_ini,fs_ini);
 		p.setE(fe_ini);
-		p.smearPos(x_disp,y_disp,r);
-		p.smearAng(tx_disp,ty_disp,r);
-		p.smearE(e_disp,r);
-		p.smearS(s_disp,r);
+		p.smearPos(x_disp,y_disp);
+		p.smearAng(tx_disp,ty_disp);
+		p.smearE(e_disp);
+		p.smearS(s_disp);
 		if (VERBOSE) {if (i==0) cout << " x_ini , tx_ini " << p.getX() << " " << p.getTX() << endl;}
 		beamParticles.push_back(p);
 	}
 }
-
-//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, TRandom* r) {
-        beamParticles.clear();
-        Nparticles = (Number_of_particles<2) ? 2 : Number_of_particles;
-        float gx,gy,gs,gtx,gty,gE;
-        for (unsigned int i=0; i<Nparticles; i++) {
-                H_BeamParticle p(p_mass,p_charge);
-                if (flat) {
-                        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 = 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);
-                p.setE(gE);
-                beamParticles.push_back(p);
-        }
-        return;
-}
-
 
 void H_Beam::createXScanningBeamParticles(const unsigned int Number_of_particles, const float fx_max) {
@@ -182,4 +149,8 @@
 }
 
+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) {
@@ -188,4 +159,10 @@
                 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);
 }
 
@@ -237,5 +214,5 @@
 	return EY2;
 }
-
+/*
 TGraphErrors * H_Beam::getBetaX(const float length, const unsigned int number_of_points) {
 	/// @param length [m]
@@ -277,4 +254,5 @@
 	return betay;
 }
+*/
 
 float H_Beam::getX(const float s, float& error_on_posx) {
@@ -313,7 +291,5 @@
 }
 
-vector<TVectorD> H_Beam::getStoppingElements(const H_AbstractBeamLine * beamline, vector<H_OpticalElement>& list, vector<int>& numb) {
-
-		vector<TVectorD> stop_positions;
+void H_Beam::getStoppingElements(const H_AbstractBeamLine * beamline, vector<H_OpticalElement>& list, vector<int>& numb) {
         vector<H_BeamParticle>::iterator particle_i;
         vector<H_OpticalElement>::iterator element_i;
@@ -331,10 +307,8 @@
                 if(particle_i->stopped(beamline)) {
                         temp_el = *(particle_i->getStoppingElement());
-						stop_positions.push_back(*(particle_i->getStopPosition()));
                         if(list.size()==0) {
                                 number=1;
                                 list.push_back(temp_el);
                                 numb.push_back(number);
-
                         } else {
                                 for (element_i = list.begin(), n_i = numb.begin(); element_i < list.end(); element_i++, n_i++) {
@@ -356,5 +330,4 @@
                 } // if particle_i->stopped
         }// for particle_i
-		return stop_positions;
 } // H_Beam::getStoppingElements
 
@@ -379,11 +352,10 @@
 }
 
-std::ostream& operator<< (std::ostream& os, const H_Beam& be) {
+void H_Beam::printProperties() const {
 	vector<H_BeamParticle>::const_iterator particle_i;
-	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;
+	cout << "There are " << Nparticles << " in the beam." << endl;
+	for (particle_i = beamParticles.begin();particle_i < beamParticles.end(); particle_i++) {
+		particle_i->printProperties();
+	}
 }
 
@@ -400,53 +372,5 @@
 	}
 } // H_Beam::printStoppingElements
-
-TH2F *  H_Beam::drawAngleProfile(const float s) {
-        /// not a const method because does a propagate to s!
-        char title[50];
-        sprintf(title,"Beam profile at %.2f m",s);
-        vector<H_BeamParticle>::iterator particle_i;
-        float xmax, xmin, ymax, ymin;
-        float xx, yy, xborder, yborder;
-
-        particle_i=beamParticles.begin();
-        xmin = particle_i->getTX();
-        xmax = particle_i->getTX();
-        ymin = particle_i->getTY();
-        ymax = particle_i->getTY();
-
-        for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
-                particle_i->propagate(s);
-                xx = particle_i->getTX();
-                yy = particle_i->getTY();
-
-                xmax = xx>xmax ? xx : xmax;
-                ymax = yy>ymax ? yy : ymax;
-                xmin = xx<xmin ? xx : xmin;
-                ymin = yy<ymin ? yy : ymin;
-        }
-
-        // in order to avoid some drawing problems, when the beam divergence is null
-        if(!(xmax || xmin)) xmax +=0.1;
-        if(!(ymax || ymin)) xmax +=0.1;
-
-        if(xmax == xmin) xmax *= 1.1;
-        if(ymax == ymin) ymax *= 1.1;
-
-        xborder = (xmax-xmin)*0.2;
-        yborder = (ymax-ymin)*0.2;
-
-        xmax += xborder;
-        xmin -= xborder;
-        ymax += yborder;
-        ymin -= yborder;
-
-        TH2F * profile = new TH2F("profile",title,10000,xmin,xmax,1000,ymin,ymax);
-        for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
-                profile->Fill(particle_i->getTX(), particle_i->getTY());
-        }
-        return profile;
-}
-
-
+/*
 TH2F *  H_Beam::drawProfile(const float s) {
 	/// not a const method because does a propagate to s!
@@ -475,7 +399,4 @@
 
 	// in order to avoid some drawing problems, when the beam divergence is null
-        if(!(xmax || xmin)) xmax +=0.1;
-        if(!(ymax || ymin)) xmax +=0.1;
-
 	if(xmax == xmin) xmax += 0.1;
 	if(ymax == ymin) ymax += 0.1;
@@ -494,6 +415,6 @@
 	}
 	return profile;
-}
-
+}*/
+/*
 TMultiGraph * H_Beam::drawBeamX(const int color) const {
 	int mycolor = color;
@@ -519,2 +440,3 @@
         return beam_profile_y;
 }
+*/
