Merge pull request #35 from MultithreadCorner/develop

Develop

Author: AAAlvesJr

examples/phase_space/phsp_chain.inl

@@ -218,12 +218,12 @@ int main(int argv, char** argc)
 		auto start = std::chrono::high_resolution_clock::now();
 
 		//generate the final state particles for B0 -> K pi J/psi
-		phsp1.Generate(B0, Chain_d.GetDecay(_0).begin(), Chain_d.GetDecay(_0).end());
+		phsp1.Generate(B0, Chain_d.GetDecays(_0).begin(), Chain_d.GetDecays(_0).end());
 
 		//pass the list of J/psi to generate the final
 		//state particles for J/psi -> mu+ mu-
-		phsp2.Generate(Chain_d.GetDecay(_0).GetDaughters(0).begin(), Chain_d.GetDecay(_0).GetDaughters(0).end(),
-				Chain_d.GetDecay(_1).begin());
+		phsp2.Generate(Chain_d.GetDecays(_0).GetDaughters(0).begin(), Chain_d.GetDecays(_0).GetDaughters(0).end(),
+				Chain_d.GetDecays(_1).begin());
 
 		auto end = std::chrono::high_resolution_clock::now();
 

examples/root_macros/particle_dispersion.C

@@ -47,92 +47,129 @@
 #include <hydra/multiarray.h>
 #include <hydra/ForEach.h>
 
+#include "TCanvas.h"
+#include "TGraphTime.h"
+#include "TROOT.h"
+#include "TMarker.h"
 
 struct Evolve
 {
 	Evolve() = delete;
 
-	Evolve(double time, double step):
-		fTime(time),
-		fStep(step)
+	Evolve(double ispeed, double time, double slope):
+		fInitialSpeed(ispeed),
+		fSlope(slope),
+		fTime(time)
 	{}
 
+	__hydra_dual__
 	Evolve(Evolve const& other):
-		fTime(other.fTime),
-		fStep(other.fStep)
+		fInitialSpeed(other.fInitialSpeed),
+		fSlope(other.fSlope),
+		fTime(other.fTime)
 	{}
 
-	__hydra_dual__
 	template<typename Particle>
+	__hydra_dual__
 	inline void operator()(Particle p) {
+
 		double theta =  hydra::get<3>(p);
-		hydra::get<0>(p) += delta_x(theta);
-		hydra::get<1>(p) += delta_y(theta);
-		hydra::get<2>(p) += delta_z(theta);
+		double phi   =  hydra::get<4>(p);
+
+		hydra::get<0>(p) += delta_x(theta, phi);
+		hydra::get<1>(p) += delta_y(theta, phi);
+		hydra::get<2>(p) += delta_z(theta, phi);
 	}
 
 	__hydra_dual__
-	inline double delta_x(double theta){
-		return ::exp(-fTime)*::cos(theta)*fStep;
+	inline double delta_x(double theta, double phi){
+		return displacement()*::cos(theta)*sin(phi);
 	}
 
 	__hydra_dual__
-	inline double delta_y(double theta){
-		return ::exp(-fTime)*::sin(theta)*fStep;
+	inline double delta_y(double theta, double phi){
+		return displacement()*::sin(theta)*sin(phi);
 	}
 
 	__hydra_dual__
-	inline double delta_z(double theta)	{
-		return ::exp(-fTime)*fStep;
+	inline double delta_z(double theta, double phi)	{
+		return displacement()*sin(phi);
 	}
 
+
+	__hydra_dual__
+	inline double displacement(){
+		return ::exp(-fSlope*fTime)*(fSlope*fTime*fTime + fInitialSpeed) ;
+	}
+
+
+	double fInitialSpeed;
+	double fSlope;
 	double fTime;
-	double fStep;
+
 };
 
 
-typedef hydra::multiarray<double, 4, hydra::device::sys_t> StateDevice_t;
-typedef hydra::multiarray<double, 4, hydra::host::sys_t>   StateHost_t;
+typedef hydra::multiarray<double, 5, hydra::device::sys_t> StateDevice_t;
+typedef hydra::multiarray<double, 5, hydra::host::sys_t>   StateHost_t;
 typedef std::vector<StateHost_t> Ensamble_t;
 
-void particle_dispersion( size_t nparicles, size_t nsteps, double epsilon ){
+void particle_dispersion( size_t nparicles, size_t nsteps, double initial_speed, double slope ){
 
-	Ensamble_t ensamble(nsteps, StateHost_t(nparicles, hydra::make_tuple(0.0, 0.0, 0.0, 0.0) ));
+	Ensamble_t ensamble(nsteps, StateHost_t(nparicles, hydra::make_tuple(0.0, 0.0, 0.0, 0.0, 0.0) ));
 
 
 	{
 		//initial state
-		StateDevice_t initial_state(nparicles, hydra::make_tuple(0.0, 0.0, 0.0, 0.0));
+		StateDevice_t initial_state(nparicles, hydra::make_tuple(0.0, 0.0, 0.0, 0.0, 0.0));
 		//distribute particles in the surface
 		hydra::Random<> RND( std::chrono::system_clock::now().time_since_epoch().count());
 		// set x-coordinate
 		RND.SetSeed(159);
-		RND.Uniform(-1.0, 1.0, initial_state.begin(0), initial_state.end(0));
+		RND.Uniform(-0.1, 0.1, initial_state.begin(0), initial_state.end(0));
 		// set y-coordinate
 		RND.SetSeed(753);
-		RND.Uniform(-1.0, 1.0, initial_state.begin(1), initial_state.end(1));
+		RND.Uniform(-2.0, 2.0, initial_state.begin(1), initial_state.end(1));
 
 		size_t i=1;
 
 		for(auto& final_state:ensamble){
 
-			hydra::Random<> RND(std::rand());
+			RND.SetSeed(std::rand());
 			RND.Uniform(0.0, 2.0*PI, initial_state.begin(3), initial_state.end(3));
-			hydra::for_each( initial_state.begin(), initial_state.end(), Evolve(i*epsilon, 2.0)  );
+			RND.SetSeed(std::rand());
+			RND.Uniform(0.0, 2.0*PI, initial_state.begin(4), initial_state.end(4));
+			hydra::for_each( initial_state.begin(), initial_state.end(), Evolve(initial_speed , double(i)/nsteps, slope)  );
 			hydra::copy(initial_state.begin(), initial_state.end(), final_state.begin());
 			i++;
 		}
 
 	}
 
-	for(auto& state:ensamble){
-		std::cout<< "=====================================" << std::endl;
-		std::cout<< "=====================================" << std::endl;
-		std::cout<< std::endl	<< std::endl	<< std::endl;
-		for(auto particle:state)
-			std::cout<< particle << std::endl;
+
+	TCanvas canvas("canvas", "", 500, 500);
+	TGraphTime *g = new TGraphTime(nsteps,-10.0, -10.0, 10.0, 10.0);
+
+	int step=0;
+	for(auto final_state:ensamble){
+
+		for( auto particle:final_state){
+
+			double x = hydra::get<0>(particle);
+			double y = hydra::get<1>(particle);
+
+			TMarker *m = new TMarker(x,y,20);
+			m->SetMarkerColor(kBlack);
+			m->SetMarkerSize(0.5);
+			g->Add(m,step);
+		}
+
+		step++;
 	}
 
+
+	g->Draw();
+
 }
 
 

hydra/Chains.h

@@ -231,13 +231,13 @@ class Chains<Decays<N, hydra::detail::BackendPolicy<BACKEND> >...> {
 
 	template<unsigned int I>
     typename HYDRA_EXTERNAL_NS::thrust::tuple_element<I, decays_type >::type&
-	GetDecay(placeholders::placeholder<I> ) {
+	GetDecays(placeholders::placeholder<I> ) {
 		return HYDRA_EXTERNAL_NS::thrust::get<I>(this->fDecays);
 	}
 
 	template<unsigned int I>
 	typename HYDRA_EXTERNAL_NS::thrust::tuple_element<I, decays_type >::type const&
-	GetDecay(placeholders::placeholder<I> ) const {
+	GetDecays(placeholders::placeholder<I> ) const {
 		return HYDRA_EXTERNAL_NS::thrust::get<I>(this->fDecays);
 	}
 

hydra/Vector4R.h

@@ -97,6 +97,10 @@ class Vector4R
 			const Vector4R& v);
 	__hydra_host__      __hydra_device__       inline GReal_t mass2() const;
 	__hydra_host__      __hydra_device__       inline GReal_t mass() const;
+	__hydra_host__      __hydra_device__       inline GReal_t p2() const;
+	__hydra_host__      __hydra_device__       inline GReal_t p() const;
+
+
 	__hydra_host__ __hydra_device__ inline void applyRotateEuler(GReal_t alpha,
 			GReal_t beta, GReal_t gamma);
 	__hydra_host__ __hydra_device__ inline void applyBoostTo(const Vector4R& p4,

hydra/detail/Vector4R.inl

@@ -385,7 +385,32 @@ inline Vector3R Vector4R::vector3()
 
 	return temp;
 }
+__hydra_host__ __hydra_device__
+inline GReal_t Vector4R::p() const
+
+// returns the 3 momentum mag.
+{
+	GReal_t temp;
+
+	temp = v[1] * v[1] + v[2] * v[2] + v[3] * v[3];
+
+	temp = sqrt(temp);
 
+	return temp;
+} //
+
+__hydra_host__ __hydra_device__
+inline GReal_t Vector4R::p2() const
+
+// returns the 3 momentum mag.
+{
+	GReal_t temp;
+
+	temp = v[1] * v[1] + v[2] * v[2] + v[3] * v[3];
+
+
+	return temp;
+} //
 
 __hydra_host__ __hydra_device__
 inline GReal_t Vector4R::d3mag() const
@@ -400,6 +425,8 @@ inline GReal_t Vector4R::d3mag() const
 
 	return temp;
 } // r3mag
+
+
 __hydra_host__ __hydra_device__
 inline GReal_t Vector4R::dot(const Vector4R& p2) const
 {