From f38e508ebaefeee7863e72e9e8893857040906df Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o?= Date: Fri, 7 Sep 2018 12:28:52 +0200 Subject: [PATCH 01/19] Add data.generate and fix label for MeV/c --- p2sat/_Data.py | 51 +++++++++++++++++++++++++++++++++++++++++++++++--- p2sat/_Plot.py | 4 ++-- 2 files changed, 50 insertions(+), 5 deletions(-) diff --git a/p2sat/_Data.py b/p2sat/_Data.py index 820e2a6..cc562a2 100644 --- a/p2sat/_Data.py +++ b/p2sat/_Data.py @@ -124,19 +124,64 @@ def update(self,w,x,y,z,px,py,pz,t,verbose=True): mass = self._ps.mass if mass == 0: self.ekin = self.p - self.gamma = self.ekin # FIXME : vérifier + self.gamma = np.inf else: self.ekin = (np.sqrt((self.p/mass)**2 + 1) - 1) * mass self.gamma = self.ekin/mass + 1. if verbose: print("Done !") - def generate(self,**kargs): + def generate(self,Nconf,Npart,ekin,theta,phi,pos=None,time=None,verbose=True): """ Generate a particle phase space from given laws TODO """ - pass + # Print a starting message + if verbose: + print("Generate particle phase-space for \"%s\" ekin law, \"%s\" theta law, \"%s\" phi law ..." + %(ekin["law"],theta["law"],phi["law"])) + + # Ensure that Nconf is of type int (for values such as 1e6) + Nconf = int(Nconf) + # Generate weights + wnorm = float(Npart)/Nconf + w = np.array([wnorm] * Nconf) + + # Generate theta angle + if theta["law"]=="iso": + theta0 = np.random.uniform(0.,np.pi,Nconf) + elif theta["law"]=="gauss": + pass + # Generate phi angle + if phi["law"]=="iso": + phi0 = np.random.uniform(0,np.pi,Nconf) + elif phi["law"]=="gauss": + pass + # Generate energy + if ekin["law"]=="mono": + ekin0 = np.array([ekin["E"]] * Nconf) + elif ekin["law"]=="exp": + pass + + # Reconstruct momentum from energy and angle distributions + mass = self._ps.mass + px = np.sqrt((ekin0/mass)**2/(1. + np.tan(theta0)**2 + np.tan(phi0)**2)) * mass + py = px * np.tan(theta0) + pz = px * np.tan(phi0) + + # Generate position + if pos is None: + x = np.array([0.] * Nconf) + y = np.array([0.] * Nconf) + z = np.array([0.] * Nconf) + + # Generate time + if time is None: + t = np.array([0.] * Nconf) + + if verbose: print("Done !") + # Update current object + self.update(w,x,y,z,px,py,pz,t,verbose=verbose) def select(self,axis,faxis,frange,fpp=1e-7): """ diff --git a/p2sat/_Plot.py b/p2sat/_Plot.py index fa4e939..1f59316 100644 --- a/p2sat/_Plot.py +++ b/p2sat/_Plot.py @@ -39,11 +39,11 @@ def get_labels(self,axes,wnorm): unit = self._ps.data.units[ax] if unit is not None: labels.append("{} ({})".format(name,unit)) - if wnorm is None:res += "/{}".format(unit) + if wnorm is None:res += "{} ".format(unit) else: labels.append(name) - labels.append("Number{}".format(res)) + labels.append("Number/({})".format(res[:-1])) return labels def figure(self,number=None,clear=True): From 276f3682d1df223768ebcf8ab50a66773a3d921a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o?= Date: Fri, 7 Sep 2018 14:37:25 +0200 Subject: [PATCH 02/19] Add doc on generate --- p2sat/_Data.py | 71 ++++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 60 insertions(+), 11 deletions(-) diff --git a/p2sat/_Data.py b/p2sat/_Data.py index cc562a2..dc3de6d 100644 --- a/p2sat/_Data.py +++ b/p2sat/_Data.py @@ -44,7 +44,7 @@ class _Data(object): - gamma is defined as - :math:`E_{kin}/m_e c^2 + 1` for massive species - - :math:`...` otherwise + - :math:`+\infty` otherwise Details of the calculations can be found at ... TODO """ @@ -132,9 +132,46 @@ def update(self,w,x,y,z,px,py,pz,t,verbose=True): def generate(self,Nconf,Npart,ekin,theta,phi,pos=None,time=None,verbose=True): """ - Generate a particle phase space from given laws - - TODO + Generate a particle phase space from given laws. + + Parameters + ---------- + Nconf : int + total number of configurations + Npart : float + total number of particles + ekin : dict + parameters to generate kinetic energy + theta : dict + parameters to generate theta angle distribution + phi : dict + parameters to generate phi angle distribution + pos : dict, optional + parameters to generate position distribution. Default is 0 for x,y,z + time : dict + parameters to generate time distribution. Default is 0 + + Notes + ----- + The dictionnaries must each time at least contain the key 'law' with a value + depending on which law are available for each physical quantity + + For dict `ekin`, available laws are : + - 'mono', for a mono-energetic source. Energy must be given as a value of keyword 'energy' + - 'exp', for exponential energy. Temperature must be given as a value of keyword 'T' + + For dict `theta` and `phi`, available laws are : + - 'mono', for a directional source. Angle must be given as a value of keyword 'angle' + - 'iso', for an isotropic source. An optional keyword 'max' can be given to specify a maximum angle + - 'gauss', for a gaussian spreading. Center of the distribution must be given with keyword 'mu', and standard deviantion with keyword 'sigma' + + Examples + -------- + Assuming a `PhaseSpace` object is instanciated as `eps`, you can generate + a mono-energetic source in isotropic direction for 1e12 particles represented + by 1e6 configurations as follows + + >>> eps.data.generate(Nconf=1e6,Npart=1e12,ekin={"law":"mono","E":20.0},theta={"law":"iso"},phi={"law":"iso"}) """ # Print a starting message if verbose: @@ -148,20 +185,32 @@ def generate(self,Nconf,Npart,ekin,theta,phi,pos=None,time=None,verbose=True): w = np.array([wnorm] * Nconf) # Generate theta angle - if theta["law"]=="iso": - theta0 = np.random.uniform(0.,np.pi,Nconf) + if theta["law"]=="mono": + theta0 = np.array([theta["angle"]] * Nconf) + elif theta["law"]=="iso": + try: + mangle = theta["max"]*np.pi/180. + except KeyError: + mangle = np.pi + theta0 = np.random.uniform(0.,mangle,Nconf) elif theta["law"]=="gauss": - pass + theta0 = np.random.normal(theta["mu"],theta["sigma"],Nconf) # Generate phi angle - if phi["law"]=="iso": - phi0 = np.random.uniform(0,np.pi,Nconf) + if phi["law"]=="mono": + phi0 = np.array([phi["angle"]] * Nconf) + elif phi["law"]=="iso": + try: + mangle = theta["max"]*np.pi/180. + except KeyError: + mangle = np.pi + phi0 = np.random.uniform(0,mangle,Nconf) elif phi["law"]=="gauss": - pass + phi0 = np.random.normal(phi["mu"],phi["sigma"],Nconf) # Generate energy if ekin["law"]=="mono": ekin0 = np.array([ekin["E"]] * Nconf) elif ekin["law"]=="exp": - pass + ekin0 = np.random.exponential(ekin["T"],Nconf) # Reconstruct momentum from energy and angle distributions mass = self._ps.mass From 33281d8f30dccdb2b1e57fa23d4932f988b32ff0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o?= Date: Fri, 7 Sep 2018 15:51:10 +0200 Subject: [PATCH 03/19] lorrentz transform prototype + fix plot against t --- p2sat/_Data.py | 77 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 77 insertions(+) diff --git a/p2sat/_Data.py b/p2sat/_Data.py index dc3de6d..00a8829 100644 --- a/p2sat/_Data.py +++ b/p2sat/_Data.py @@ -65,6 +65,7 @@ def __init__(self,PhaseSpace): 'px' : '$p_x$', 'py' : '$p_y$', 'pz' : '$p_z$', + 't' : '$t$', 'r' : '$r$', 'p' : '$p$', @@ -82,6 +83,7 @@ def __init__(self,PhaseSpace): 'px' : 'MeV/c', 'py' : 'MeV/c', 'pz' : 'MeV/c', + 't' : 'fs', 'r' : 'um', 'p' : 'MeV/c', @@ -313,6 +315,81 @@ def propagate(self): """ pass + def lorrentz(self,beta_CM,pos_CM=None): + """ + Lorrentz-transformate the particle phase-space with given speed of the center of mass. + """ + # If no position given, default is the origin + if pos_CM is None: + pos_CM = [0.,0.,0.] + # Calculate gamma factor from beta_CM + # beta_CM = np.array(beta_CM) + # v_CM = beta_CM * 2.99792458e8 * 1e6/1e15 # speed of CM in um/fs + # gamma_CM = 1./np.sqrt(1-beta_CM**2) + # # Time transformation + # t = gamma_CM*(self.t + # # Position transformation + # # x = (self.x - pos_CM[0])*gamma_CM[0] + # # y = (self.y - pos_CM[1])*gamma_CM[1] + # # z = (self.z - pos_CM[2])*gamma_CM[2] + # x = self.x + # y = self.y + # z = self.z + # # Momentum transformation + # px = self.px*gamma_CM[0] + # py = self.px*gamma_CM[1] + # pz = self.px*gamma_CM[2] + # # No transformations on statistical weights + # w = self.w + beta_CM = np.array(beta_CM) + b = np.dot(beta_CM,beta_CM) + n = beta_CM/b + c = 2.99792458e8 * 1e6/1e15 # speed of CM in um/fs + v = b * c + gamma = 1./np.sqrt(1-beta_CM**2) + g = 1./np.sqrt(1-b**2) + + G = np.matrix([ + [g ,-g*b*n[0] ,-g*b*n[1] ,-g*b*n[2] ], + [-g*b*n[0] ,1+(g-1)*n[0]**2 ,(g-1)*n[0]*n[1] ,(g-1)*n[0]*n[2]], + [-g*b*n[1] ,(g-1)*n[1]*n[0] ,1+(g-1)*n[1]**2 ,(g-1)*n[1]*n[2]], + [-g*b*n[2] ,(g-1)*n[2]*n[0] ,(g-1)*n[2]*n[1] ,1+(g-1)*n[2]**2] + ]) + print(G) + # Time transformation + q1 = [c*self.t.copy(),self.x.copy(),self.y.copy(),self.z.copy()] + p1 = [self.ekin.copy()/c,self.px.copy(),self.py.copy(),self.pz.copy()] + + q2 = np.dot(G,q1) + p2 = np.dot(G,p1) + # print(q2) + # print() + # print(p2) + w = self.w + t = np.array(q2[0]/c)[0] + x = np.array(q2[1])[0] + y = np.array(q2[2])[0] + z = np.array(q2[3])[0] + px = np.array(p1[1])[0] + py = np.array(p1[2])[0] + pz = np.array(p1[3])[0] + # print(x,y) + # self.update(w,x[0],y[0],z[0],px[0],py[0],pz[0],t[0]) + self.update(w,x,y,z,px,py,pz,t) + # # Position transformation + # # x = (self.x - pos_CM[0])*gamma_CM[0] + # # y = (self.y - pos_CM[1])*gamma_CM[1] + # # z = (self.z - pos_CM[2])*gamma_CM[2] + # x = self.x + # y = self.y + # z = self.z + # # Momentum transformation + # px = self.px*gamma_CM[0] + # py = self.px*gamma_CM[1] + # pz = self.px*gamma_CM[2] + # # No transformations on statistical weights + # w = self.w + def discretize(self,with_time=True,verbose=True,**kargs): """ Discretize the particles phase space in a 6 or 7 D histogram. From 00c9fd0a9598786255fe6caca68e6c7078106eff Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o?= Date: Fri, 7 Sep 2018 16:25:25 +0200 Subject: [PATCH 04/19] Add __str__ method to PhaseSpace object --- p2sat/PhaseSpace.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/p2sat/PhaseSpace.py b/p2sat/PhaseSpace.py index b8f50b0..d11ec76 100644 --- a/p2sat/PhaseSpace.py +++ b/p2sat/PhaseSpace.py @@ -73,6 +73,28 @@ def __add__(self,other): return ps + def __str__(self): + txt = "" + txt += "p2sat PhaseSpace instance located at %s\n\n"%hex(id(self)) + txt += "Specie : %s\n"%self.specie + txt += "Number of configurations : %i\n"%len(self) + txt += "Total number of particles : %.4E\n\n"%sum(self.data.w) + txt += "Statistics : ( min , max , mean , std ) unit\n" + txt += " x : (% .3E, % .3E, % .3E, % .3E) um\n"%(min(self.data.x),max(self.data.x),self.data.x.mean(),self.data.x.std()) + txt += " y : (% .3E, % .3E, % .3E, % .3E) um\n"%(min(self.data.y),max(self.data.y),self.data.y.mean(),self.data.y.std()) + txt += " z : (% .3E, % .3E, % .3E, % .3E) um\n"%(min(self.data.z),max(self.data.z),self.data.z.mean(),self.data.z.std()) + txt += " r : (% .3E, % .3E, % .3E, % .3E) um\n\n"%(min(self.data.r),max(self.data.r),self.data.r.mean(),self.data.r.std()) + txt += " px : (% .3E, % .3E, % .3E, % .3E) MeV/c\n"%(min(self.data.px),max(self.data.px),self.data.px.mean(),self.data.px.std()) + txt += " py : (% .3E, % .3E, % .3E, % .3E) MeV/c\n"%(min(self.data.py),max(self.data.py),self.data.py.mean(),self.data.py.std()) + txt += " pz : (% .3E, % .3E, % .3E, % .3E) MeV/c\n"%(min(self.data.pz),max(self.data.pz),self.data.pz.mean(),self.data.pz.std()) + txt += " p : (% .3E, % .3E, % .3E, % .3E) MeV/c\n\n"%(min(self.data.p),max(self.data.p),self.data.p.mean(),self.data.p.std()) + txt += " ekin : (% .3E, % .3E, % .3E, % .3E) MeV\n"%(min(self.data.ekin),max(self.data.ekin),self.data.ekin.mean(),self.data.ekin.std()) + txt += " theta : (% .3E, % .3E, % .3E, % .3E) deg\n"%(min(self.data.theta),max(self.data.theta),self.data.theta.mean(),self.data.theta.std()) + txt += " phi : (% .3E, % .3E, % .3E, % .3E) deg\n"%(min(self.data.phi),max(self.data.phi),self.data.phi.mean(),self.data.phi.std()) + txt += "" + + return txt + def __len__(self): """ Return the total number of configurations. From ab4c71d82005fda5197be2e9136cc164596418bc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o?= Date: Fri, 7 Sep 2018 16:58:58 +0200 Subject: [PATCH 05/19] Fix generate : right energy --- p2sat/_Data.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/p2sat/_Data.py b/p2sat/_Data.py index 00a8829..b346abc 100644 --- a/p2sat/_Data.py +++ b/p2sat/_Data.py @@ -128,8 +128,10 @@ def update(self,w,x,y,z,px,py,pz,t,verbose=True): self.ekin = self.p self.gamma = np.inf else: - self.ekin = (np.sqrt((self.p/mass)**2 + 1) - 1) * mass - self.gamma = self.ekin/mass + 1. + # self.ekin = (np.sqrt((self.p/mass)**2 + 1) - 1) * mass + # self.gamma = self.ekin/mass + 1. + self.gamma = np.sqrt((self.p/mass)**2 + 1.) + self.ekin = (self.gamma - 1) * mass if verbose: print("Done !") def generate(self,Nconf,Npart,ekin,theta,phi,pos=None,time=None,verbose=True): @@ -216,7 +218,8 @@ def generate(self,Nconf,Npart,ekin,theta,phi,pos=None,time=None,verbose=True): # Reconstruct momentum from energy and angle distributions mass = self._ps.mass - px = np.sqrt((ekin0/mass)**2/(1. + np.tan(theta0)**2 + np.tan(phi0)**2)) * mass + etot = ekin0 + mass + px = np.sqrt((etot**2 - mass**2)/(1. + np.tan(theta0)**2 + np.tan(phi0)**2)) py = px * np.tan(theta0) pz = px * np.tan(phi0) From 2542db4aedf6847a8c699c66dac8cc504e9e2c31 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o?= Date: Fri, 7 Sep 2018 17:14:22 +0200 Subject: [PATCH 06/19] WIP: Update lorrentz --- p2sat/_Data.py | 90 +++++++++++++++++++++++++------------------------- 1 file changed, 45 insertions(+), 45 deletions(-) diff --git a/p2sat/_Data.py b/p2sat/_Data.py index b346abc..a14718f 100644 --- a/p2sat/_Data.py +++ b/p2sat/_Data.py @@ -128,10 +128,8 @@ def update(self,w,x,y,z,px,py,pz,t,verbose=True): self.ekin = self.p self.gamma = np.inf else: - # self.ekin = (np.sqrt((self.p/mass)**2 + 1) - 1) * mass - # self.gamma = self.ekin/mass + 1. - self.gamma = np.sqrt((self.p/mass)**2 + 1.) - self.ekin = (self.gamma - 1) * mass + self.ekin = (np.sqrt((self.p/mass)**2 + 1) - 1) * mass + self.gamma = self.ekin/mass + 1. if verbose: print("Done !") def generate(self,Nconf,Npart,ekin,theta,phi,pos=None,time=None,verbose=True): @@ -344,54 +342,56 @@ def lorrentz(self,beta_CM,pos_CM=None): # pz = self.px*gamma_CM[2] # # No transformations on statistical weights # w = self.w + + # beta_CM = np.array(beta_CM) + # b = np.dot(beta_CM,beta_CM) + # n = beta_CM/b + # c = 2.99792458e8 * 1e6/1e15 # speed of CM in um/fs + # v = b * c + # gamma = 1./np.sqrt(1-beta_CM**2) + # g = 1./np.sqrt(1-b**2) + # G = np.matrix([ + # [g ,-g*b*n[0] ,-g*b*n[1] ,-g*b*n[2] ], + # [-g*b*n[0] ,1+(g-1)*n[0]**2 ,(g-1)*n[0]*n[1] ,(g-1)*n[0]*n[2]], + # [-g*b*n[1] ,(g-1)*n[1]*n[0] ,1+(g-1)*n[1]**2 ,(g-1)*n[1]*n[2]], + # [-g*b*n[2] ,(g-1)*n[2]*n[0] ,(g-1)*n[2]*n[1] ,1+(g-1)*n[2]**2] + # ]) + # print(G) + # # Time transformation + # q1 = [c*self.t.copy(),self.x.copy(),self.y.copy(),self.z.copy()] + # p1 = [self.ekin.copy()/c,self.px.copy(),self.py.copy(),self.pz.copy()] + # q2 = np.dot(G,q1) + # p2 = np.dot(G,p1) + # w = self.w + # t = np.array(q2[0]/c)[0] + # x = np.array(q2[1])[0] + # y = np.array(q2[2])[0] + # z = np.array(q2[3])[0] + # px = np.array(p1[1])[0] + # py = np.array(p1[2])[0] + # pz = np.array(p1[3])[0] + # self.update(w,x,y,z,px,py,pz,t) + beta_CM = np.array(beta_CM) b = np.dot(beta_CM,beta_CM) n = beta_CM/b c = 2.99792458e8 * 1e6/1e15 # speed of CM in um/fs - v = b * c + v = beta_CM * c gamma = 1./np.sqrt(1-beta_CM**2) g = 1./np.sqrt(1-b**2) - G = np.matrix([ - [g ,-g*b*n[0] ,-g*b*n[1] ,-g*b*n[2] ], - [-g*b*n[0] ,1+(g-1)*n[0]**2 ,(g-1)*n[0]*n[1] ,(g-1)*n[0]*n[2]], - [-g*b*n[1] ,(g-1)*n[1]*n[0] ,1+(g-1)*n[1]**2 ,(g-1)*n[1]*n[2]], - [-g*b*n[2] ,(g-1)*n[2]*n[0] ,(g-1)*n[2]*n[1] ,1+(g-1)*n[2]**2] - ]) - print(G) - # Time transformation - q1 = [c*self.t.copy(),self.x.copy(),self.y.copy(),self.z.copy()] - p1 = [self.ekin.copy()/c,self.px.copy(),self.py.copy(),self.pz.copy()] - - q2 = np.dot(G,q1) - p2 = np.dot(G,p1) - # print(q2) - # print() - # print(p2) - w = self.w - t = np.array(q2[0]/c)[0] - x = np.array(q2[1])[0] - y = np.array(q2[2])[0] - z = np.array(q2[3])[0] - px = np.array(p1[1])[0] - py = np.array(p1[2])[0] - pz = np.array(p1[3])[0] - # print(x,y) - # self.update(w,x[0],y[0],z[0],px[0],py[0],pz[0],t[0]) - self.update(w,x,y,z,px,py,pz,t) - # # Position transformation - # # x = (self.x - pos_CM[0])*gamma_CM[0] - # # y = (self.y - pos_CM[1])*gamma_CM[1] - # # z = (self.z - pos_CM[2])*gamma_CM[2] - # x = self.x - # y = self.y - # z = self.z - # # Momentum transformation - # px = self.px*gamma_CM[0] - # py = self.px*gamma_CM[1] - # pz = self.px*gamma_CM[2] - # # No transformations on statistical weights - # w = self.w + w1 = self.w + t1 = self.t + r1 = np.array([self.x,self.y,self.z]) + p1 = np.array([self.px,self.py,self.pz]) + + w2 = w1 + t2 = g*(t1 + sum(r1*v)/c**2) + r2 = r1 + ((g-1)/np.dot(v,v) * sum(r1*v) + g*t1)*v + p2 = p1 + + self.update(w2,r2[0],r2[1],r2[2],p2[0],p2[1],p2[2],t2) + def discretize(self,with_time=True,verbose=True,**kargs): """ From 963969d5d2a3ef94a8e434bb0e2ba63f4822c636 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o?= Date: Mon, 10 Sep 2018 11:46:15 +0200 Subject: [PATCH 07/19] Update generate doc & fix max iso for phi --- p2sat/_Data.py | 99 +++++++++++++++++++------------------------------- 1 file changed, 37 insertions(+), 62 deletions(-) diff --git a/p2sat/_Data.py b/p2sat/_Data.py index a14718f..cc6f5e0 100644 --- a/p2sat/_Data.py +++ b/p2sat/_Data.py @@ -167,6 +167,21 @@ def generate(self,Nconf,Npart,ekin,theta,phi,pos=None,time=None,verbose=True): - 'iso', for an isotropic source. An optional keyword 'max' can be given to specify a maximum angle - 'gauss', for a gaussian spreading. Center of the distribution must be given with keyword 'mu', and standard deviantion with keyword 'sigma' + Details of the calculations : + + Considering :math:`E_T = E_k + E_m` being the total energy, with + :math:`E_k` the kinetic energy and :math:`E_m` the rest mass energy, we have + :math:`E_T^2 = p^2 + E_m^2` and :math:`p^2=p_x^2 + p_y^2 + p_z^2` so + :math:`p_x^2+p_y^2+p_z^2=E_T^2 - E_m^2`. + + Assuming :math:`\\tan{\\theta}=\\frac{p_y}{p_x}` and + :math:`\\tan{\\phi}=\\frac{p_z}{p_x}` we get + + :math:`p_x = \sqrt{\\frac{E_T^2 - E_m^2}{1 + \\tan{\\theta}^2 + \\tan{\phi}^2}}` + :math:`p_y = p_x \\tan{\\theta}` + :math:`p_z = p_x \\tan{\phi}` + + Examples -------- Assuming a `PhaseSpace` object is instanciated as `eps`, you can generate @@ -202,7 +217,7 @@ def generate(self,Nconf,Npart,ekin,theta,phi,pos=None,time=None,verbose=True): phi0 = np.array([phi["angle"]] * Nconf) elif phi["law"]=="iso": try: - mangle = theta["max"]*np.pi/180. + mangle = phi["max"]*np.pi/180. except KeyError: mangle = np.pi phi0 = np.random.uniform(0,mangle,Nconf) @@ -323,75 +338,35 @@ def lorrentz(self,beta_CM,pos_CM=None): # If no position given, default is the origin if pos_CM is None: pos_CM = [0.,0.,0.] - # Calculate gamma factor from beta_CM - # beta_CM = np.array(beta_CM) - # v_CM = beta_CM * 2.99792458e8 * 1e6/1e15 # speed of CM in um/fs - # gamma_CM = 1./np.sqrt(1-beta_CM**2) - # # Time transformation - # t = gamma_CM*(self.t - # # Position transformation - # # x = (self.x - pos_CM[0])*gamma_CM[0] - # # y = (self.y - pos_CM[1])*gamma_CM[1] - # # z = (self.z - pos_CM[2])*gamma_CM[2] - # x = self.x - # y = self.y - # z = self.z - # # Momentum transformation - # px = self.px*gamma_CM[0] - # py = self.px*gamma_CM[1] - # pz = self.px*gamma_CM[2] - # # No transformations on statistical weights - # w = self.w - - # beta_CM = np.array(beta_CM) - # b = np.dot(beta_CM,beta_CM) - # n = beta_CM/b - # c = 2.99792458e8 * 1e6/1e15 # speed of CM in um/fs - # v = b * c - # gamma = 1./np.sqrt(1-beta_CM**2) - # g = 1./np.sqrt(1-b**2) - # G = np.matrix([ - # [g ,-g*b*n[0] ,-g*b*n[1] ,-g*b*n[2] ], - # [-g*b*n[0] ,1+(g-1)*n[0]**2 ,(g-1)*n[0]*n[1] ,(g-1)*n[0]*n[2]], - # [-g*b*n[1] ,(g-1)*n[1]*n[0] ,1+(g-1)*n[1]**2 ,(g-1)*n[1]*n[2]], - # [-g*b*n[2] ,(g-1)*n[2]*n[0] ,(g-1)*n[2]*n[1] ,1+(g-1)*n[2]**2] - # ]) - # print(G) - # # Time transformation - # q1 = [c*self.t.copy(),self.x.copy(),self.y.copy(),self.z.copy()] - # p1 = [self.ekin.copy()/c,self.px.copy(),self.py.copy(),self.pz.copy()] - # q2 = np.dot(G,q1) - # p2 = np.dot(G,p1) - # w = self.w - # t = np.array(q2[0]/c)[0] - # x = np.array(q2[1])[0] - # y = np.array(q2[2])[0] - # z = np.array(q2[3])[0] - # px = np.array(p1[1])[0] - # py = np.array(p1[2])[0] - # pz = np.array(p1[3])[0] - # self.update(w,x,y,z,px,py,pz,t) beta_CM = np.array(beta_CM) b = np.dot(beta_CM,beta_CM) n = beta_CM/b c = 2.99792458e8 * 1e6/1e15 # speed of CM in um/fs - v = beta_CM * c + v = b * c gamma = 1./np.sqrt(1-beta_CM**2) g = 1./np.sqrt(1-b**2) + G = np.matrix([ + [g ,-g*b*n[0] ,-g*b*n[1] ,-g*b*n[2] ], + [-g*b*n[0] ,1+(g-1)*n[0]**2 ,(g-1)*n[0]*n[1] ,(g-1)*n[0]*n[2]], + [-g*b*n[1] ,(g-1)*n[1]*n[0] ,1+(g-1)*n[1]**2 ,(g-1)*n[1]*n[2]], + [-g*b*n[2] ,(g-1)*n[2]*n[0] ,(g-1)*n[2]*n[1] ,1+(g-1)*n[2]**2] + ]) - w1 = self.w - t1 = self.t - r1 = np.array([self.x,self.y,self.z]) - p1 = np.array([self.px,self.py,self.pz]) - - w2 = w1 - t2 = g*(t1 + sum(r1*v)/c**2) - r2 = r1 + ((g-1)/np.dot(v,v) * sum(r1*v) + g*t1)*v - p2 = p1 - - self.update(w2,r2[0],r2[1],r2[2],p2[0],p2[1],p2[2],t2) - + # # Time transformation + q1 = [c*self.t.copy(),self.x.copy(),self.y.copy(),self.z.copy()] + p1 = [self.ekin.copy()/c,self.px.copy(),self.py.copy(),self.pz.copy()] + ct,x,y,z = np.dot(G,q1) + ekin_c,px,py,pz = np.dot(G,p1) + w = self.w + t = np.array(ct)[0]/c + x = np.array(x)[0] + y = np.array(y)[0] + z = np.array(z)[0] + px = np.array(px)[0] + py = np.array(py)[0] + pz = np.array(pz)[0] + self.update(w,x,y,z,px,py,pz,t) def discretize(self,with_time=True,verbose=True,**kargs): """ From e976990f7f92a582efb3e55928b9d6c15f4706ed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o?= Date: Mon, 10 Sep 2018 12:29:35 +0200 Subject: [PATCH 08/19] Fix: add plt.show at the end of plots functions --- p2sat/_Plot.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/p2sat/_Plot.py b/p2sat/_Plot.py index 1f59316..81fd3dd 100644 --- a/p2sat/_Plot.py +++ b/p2sat/_Plot.py @@ -147,6 +147,8 @@ def h1(self,axis,where='post',log=False,polar=False,reverse=False,**kargs): a.grid(True) + plt.show() + return a def f1(self,axis,func_name,log=False,polar=False,reverse=False,**kargs): @@ -205,6 +207,8 @@ def f1(self,axis,func_name,log=False,polar=False,reverse=False,**kargs): a.grid(True) + plt.show() + return a def h2(self,axis1,axis2,log=False,polar=False,**kargs): @@ -253,6 +257,8 @@ def h2(self,axis1,axis2,log=False,polar=False,**kargs): a.grid(True) plt.colorbar(a2,label=labels[2]) + plt.show() + return a def c2(self,axis1,axis2,log=False,polar=False,gfilter=0.0,**kargs): @@ -309,6 +315,8 @@ def c2(self,axis1,axis2,log=False,polar=False,gfilter=0.0,**kargs): plt.clabel(a2, inline=1, fontsize=10 ,fmt='%1.1e') + plt.show() + return a def s2(self,axis1,axis2,log=False,polar=False,select=None): @@ -362,6 +370,8 @@ def s2(self,axis1,axis2,log=False,polar=False,select=None): plt.colorbar(a2,label=labels[2]) + plt.show() + return a def h2h1(self,axis1,axis2,log=False,**kargs): @@ -403,6 +413,8 @@ def h2h1(self,axis1,axis2,log=False,**kargs): self.autoclear=tmp + plt.show() + def s2h1(self,axis1,axis2,log=False): """ """ @@ -441,3 +453,5 @@ def h3(self,axis1,axis2,axis3, tmp[i1][i2][i3]=e3 a.scatter3D(g1,g2,g3,s=snorm*tmp,cmap='hot') + + plt.show() From 0eb0ce0fa1206a8b494ff346ce9646181238e10a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o?= Date: Mon, 10 Sep 2018 12:38:38 +0200 Subject: [PATCH 09/19] Fix?: lorrentz --- p2sat/_Data.py | 89 ++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 64 insertions(+), 25 deletions(-) diff --git a/p2sat/_Data.py b/p2sat/_Data.py index cc6f5e0..9413bd8 100644 --- a/p2sat/_Data.py +++ b/p2sat/_Data.py @@ -168,7 +168,7 @@ def generate(self,Nconf,Npart,ekin,theta,phi,pos=None,time=None,verbose=True): - 'gauss', for a gaussian spreading. Center of the distribution must be given with keyword 'mu', and standard deviantion with keyword 'sigma' Details of the calculations : - + Considering :math:`E_T = E_k + E_m` being the total energy, with :math:`E_k` the kinetic energy and :math:`E_m` the rest mass energy, we have :math:`E_T^2 = p^2 + E_m^2` and :math:`p^2=p_x^2 + p_y^2 + p_z^2` so @@ -334,39 +334,78 @@ def propagate(self): def lorrentz(self,beta_CM,pos_CM=None): """ Lorrentz-transformate the particle phase-space with given speed of the center of mass. + + Notes + ----- + https://en.wikipedia.org/wiki/Lorentz_transformation#Transformation_of_other_quantities """ # If no position given, default is the origin if pos_CM is None: pos_CM = [0.,0.,0.] - beta_CM = np.array(beta_CM) - b = np.dot(beta_CM,beta_CM) - n = beta_CM/b - c = 2.99792458e8 * 1e6/1e15 # speed of CM in um/fs + # lowercase : scalar, caption : vector + + B = np.array(beta_CM) + b = np.dot(B,B) + N = B/b + c = 2.99792458e8 * 1e6/1e15 # speed of light in um/fs v = b * c - gamma = 1./np.sqrt(1-beta_CM**2) g = 1./np.sqrt(1-b**2) - G = np.matrix([ - [g ,-g*b*n[0] ,-g*b*n[1] ,-g*b*n[2] ], - [-g*b*n[0] ,1+(g-1)*n[0]**2 ,(g-1)*n[0]*n[1] ,(g-1)*n[0]*n[2]], - [-g*b*n[1] ,(g-1)*n[1]*n[0] ,1+(g-1)*n[1]**2 ,(g-1)*n[1]*n[2]], - [-g*b*n[2] ,(g-1)*n[2]*n[0] ,(g-1)*n[2]*n[1] ,1+(g-1)*n[2]**2] - ]) - - # # Time transformation - q1 = [c*self.t.copy(),self.x.copy(),self.y.copy(),self.z.copy()] - p1 = [self.ekin.copy()/c,self.px.copy(),self.py.copy(),self.pz.copy()] - ct,x,y,z = np.dot(G,q1) - ekin_c,px,py,pz = np.dot(G,p1) + + # Four position + a1 = c*self.t + Z1 = np.array([self.x,self.y,self.z]).T + + a2,Z2 =[],[] + for i in range(len(a1)): + a2.append(g*(a1[i] - (v*np.dot(N,Z1[i]))/c)) + Z2.append(Z1[i] + (g-1)*np.dot(Z1[i],N)*N - g*(a1[i]*v*N)/c) + + t = np.array(a2)/c + x,y,z = np.array(Z2).T + + # Four momentum + a1 = self.ekin/c + Z1 = np.array([self.px,self.py,self.pz]).T + + a2,Z2 =[],[] + for i in range(len(a1)): + a2.append(g*(a1[i] - (v*np.dot(N,Z1[i]))/c)) + Z2.append(Z1[i] + (g-1)*np.dot(Z1[i],N)*N - g*(a1[i]*v*N)/c) + + px,py,pz = np.array(Z2).T + w = self.w - t = np.array(ct)[0]/c - x = np.array(x)[0] - y = np.array(y)[0] - z = np.array(z)[0] - px = np.array(px)[0] - py = np.array(py)[0] - pz = np.array(pz)[0] self.update(w,x,y,z,px,py,pz,t) + # + # beta_CM = np.array(beta_CM) + # b = np.dot(beta_CM,beta_CM) + # n = beta_CM/b + # c = 2.99792458e8 * 1e6/1e15 # speed of CM in um/fs + # v = b * c + # gamma = 1./np.sqrt(1-beta_CM**2) + # g = 1./np.sqrt(1-b**2) + # G = np.matrix([ + # [g ,-g*b*n[0] ,-g*b*n[1] ,-g*b*n[2] ], + # [-g*b*n[0] ,1+(g-1)*n[0]**2 ,(g-1)*n[0]*n[1] ,(g-1)*n[0]*n[2]], + # [-g*b*n[1] ,(g-1)*n[1]*n[0] ,1+(g-1)*n[1]**2 ,(g-1)*n[1]*n[2]], + # [-g*b*n[2] ,(g-1)*n[2]*n[0] ,(g-1)*n[2]*n[1] ,1+(g-1)*n[2]**2] + # ]) + # + # # # Time transformation + # q1 = [c*self.t.copy(),self.x.copy(),self.y.copy(),self.z.copy()] + # p1 = [self.ekin.copy()/c,self.px.copy(),self.py.copy(),self.pz.copy()] + # ct,x,y,z = np.dot(G,q1) + # ekin_c,px,py,pz = np.dot(G,p1) + # w = self.w + # t = np.array(ct)[0]/c + # x = np.array(x)[0] + # y = np.array(y)[0] + # z = np.array(z)[0] + # px = np.array(px)[0] + # py = np.array(py)[0] + # pz = np.array(pz)[0] + # self.update(w,x,y,z,px,py,pz,t) def discretize(self,with_time=True,verbose=True,**kargs): """ From 68a7269c206eb89f828bcb30d3c8b9c6352a78e1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o?= Date: Mon, 10 Sep 2018 15:17:29 +0200 Subject: [PATCH 10/19] Update get_labels method --- p2sat/PhaseSpace.py | 14 +++++++++++--- p2sat/_Data.py | 9 +++++++-- p2sat/_Plot.py | 33 +++++++++++++++++++++++++-------- 3 files changed, 43 insertions(+), 13 deletions(-) diff --git a/p2sat/PhaseSpace.py b/p2sat/PhaseSpace.py index d11ec76..da6cdc6 100644 --- a/p2sat/PhaseSpace.py +++ b/p2sat/PhaseSpace.py @@ -38,12 +38,20 @@ class PhaseSpace(object): See sub-objects documentation for more informations """ def __init__(self,specie): - self.specie = specie if specie in ("gamma","g"): + self.specie = "gamma" self.mass = 0 - elif specie in ("positron","electron","e-","e+","e"): + elif specie in ("positron","e+"): + self.specie = "e+" self.mass = 0.511 - elif specie in ("muon","muon+","muon-","mu+","mu-","mu"): + elif specie in ("electron","e-"): + self.specie = "e-" + self.mass = 0.511 + elif specie in ("muon+","mu+"): + self.specie = "mu+" + self.mass = 105.6 + elif specie in ("muon-","mu-"): + self.specie = "mu-" self.mass = 105.6 else: raise NameError("Unknown particle specie.") diff --git a/p2sat/_Data.py b/p2sat/_Data.py index 9413bd8..f2b2af7 100644 --- a/p2sat/_Data.py +++ b/p2sat/_Data.py @@ -73,7 +73,12 @@ def __init__(self,PhaseSpace): 'gamma' : '$\gamma$', 'theta' : '$\\theta$', - 'phi' : '$\\phi$' + 'phi' : '$\\phi$', + + 'e-' : '$e^-$', + 'e+' : '$e^+$', + 'mu-' : '$\mu^-$', + 'mu+' : '$\mu^+$' } self.units = { @@ -372,7 +377,7 @@ def lorrentz(self,beta_CM,pos_CM=None): for i in range(len(a1)): a2.append(g*(a1[i] - (v*np.dot(N,Z1[i]))/c)) Z2.append(Z1[i] + (g-1)*np.dot(Z1[i],N)*N - g*(a1[i]*v*N)/c) - + px,py,pz = np.array(Z2).T w = self.w diff --git a/p2sat/_Plot.py b/p2sat/_Plot.py index 81fd3dd..d3fe8f1 100644 --- a/p2sat/_Plot.py +++ b/p2sat/_Plot.py @@ -29,21 +29,38 @@ def get_labels(self,axes,wnorm): labels : list of str Labels of given axes and label of weight """ + # Initialization + specie = self._ps.specie labels=[] - res="" + N_name="" + N_unit="" + # Loop over all the axes for ax in axes: + # No label if ax is not a str if type(ax) is not str: labels.append("") else: - name = self._ps.data.labels[ax] - unit = self._ps.data.units[ax] - if unit is not None: - labels.append("{} ({})".format(name,unit)) - if wnorm is None:res += "{} ".format(unit) + # Else get axis name and unit from dicts + ax_name = self._ps.data.labels[ax] + ax_unit = self._ps.data.units[ax] + # ekin is E_\gamma when specie is gamma + if ax=="ekin" and specie =="gamma": + ax_name = "$E_\gamma$" + # Format the label for axis and unit of N + if ax_unit is not None: + labels.append("{} ({})".format(ax_name,ax_unit)) + if wnorm is None:N_name += "d%s "%(ax_name[1:-1]) # LaTeX without '$' + if wnorm is None:N_unit += "%s "%(ax_unit) else: - labels.append(name) + labels.append(ax_name) + + # Format number name + specie_name = self._ps.data.labels[specie] + if N_unit =="": + labels.append("$N_{%s}$"%(specie_name[1:-1])) + else: + labels.append("$\\frac{d N_{%s}}{%s}$ (%s)$^{-1}$"%(specie_name[1:-1],N_name[:-1],N_unit[:-1])) - labels.append("Number/({})".format(res[:-1])) return labels def figure(self,number=None,clear=True): From 142f5aadcad77005b70bacd7f8f2f50bb8ecc887 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o?= Date: Mon, 10 Sep 2018 15:34:29 +0200 Subject: [PATCH 11/19] Change wnorm -> normed --- p2sat/_Data.py | 10 +++++----- p2sat/_Hist.py | 31 +++++++++++++++++++------------ p2sat/_Plot.py | 20 ++++++++++---------- 3 files changed, 34 insertions(+), 27 deletions(-) diff --git a/p2sat/_Data.py b/p2sat/_Data.py index f2b2af7..2c124f9 100644 --- a/p2sat/_Data.py +++ b/p2sat/_Data.py @@ -203,8 +203,8 @@ def generate(self,Nconf,Npart,ekin,theta,phi,pos=None,time=None,verbose=True): # Ensure that Nconf is of type int (for values such as 1e6) Nconf = int(Nconf) # Generate weights - wnorm = float(Npart)/Nconf - w = np.array([wnorm] * Nconf) + weight = float(Npart)/Nconf + w = np.array([weight] * Nconf) # Generate theta angle if theta["law"]=="mono": @@ -437,12 +437,12 @@ def discretize(self,with_time=True,verbose=True,**kargs): hn=self._ps.hist.hn if verbose : print('Data discretization ...') - #bi,hi=hn(['x','y','z','px','py','pz'],bwidth=bwidth,brange=brange,wnorm=[1.0]*6,select=select) + if with_time: - bi,hi=hn([self.x,self.y,self.z,self.px,self.py,self.pz,self.t],wnorm=[1.0]*7,**kargs) + bi,hi=hn([self.x,self.y,self.z,self.px,self.py,self.pz,self.t],normed=False,**kargs) bx,by,bz,bpx,bpy,bpz,bt=bi else: - bi,hi=hn([self.x,self.y,self.z,self.px,self.py,self.pz],wnorm=[1.0]*6,**kargs) + bi,hi=hn([self.x,self.y,self.z,self.px,self.py,self.pz],normed=False,**kargs) bx,by,bz,bpx,bpy,bpz=bi if verbose : print('Done !') w = [] diff --git a/p2sat/_Hist.py b/p2sat/_Hist.py index bacb625..9ef4e09 100644 --- a/p2sat/_Hist.py +++ b/p2sat/_Hist.py @@ -8,7 +8,7 @@ class _Hist(object): def __init__(self,PhaseSpace): self._ps=PhaseSpace - def hn(self,axis,blen=None,bwidth=None,brange=None,wnorm=None,select=None): + def hn(self,axis,blen=None,bwidth=None,brange=None,normed=True,select=None): """ Create and return the n-dimensional histo of axis list. @@ -22,8 +22,8 @@ def hn(self,axis,blen=None,bwidth=None,brange=None,wnorm=None,select=None): list of bin width. If a bwidth element is None, a calculation is done to have 10 bins in the correspondant axis brange : list of list of 2 float, optional list of bin minimum and maximum. If a brange element is None, the minimum/maximum of the axis is taken - wnorm : list of float, optional - weight normalization. If a wnorm element is None, the bin width is taken + normed : list of bool, optional + weight normalization. If a normed element is True, the bin width is taken select : dict, optional filtering dictionary @@ -44,10 +44,10 @@ def hn(self,axis,blen=None,bwidth=None,brange=None,wnorm=None,select=None): Examples -------- - >>> hn(['x'],bwidth=[50],brange=[[0,1000]],wnorm=[1.0],select={'ekin':(0.511,None)}) + >>> hn(['x'],bwidth=[50],brange=[[0,1000]],normed=[True],select={'ekin':(0.511,None)}) returns the number of particles with :math:`ekin \in [0.511, +\infty] MeV` in function of x - wnorm=[1.0] to not divide nb of particles by bin width (otherwise number per um) + normed=[True] to not divide nb of particles by bin width (otherwise number per um) >>> hn(['r','ekin'],bwidth=[10.0,0.1],brange=[[0,1000],[0.1,50.0]],select={'x':150}) @@ -70,7 +70,12 @@ def hn(self,axis,blen=None,bwidth=None,brange=None,wnorm=None,select=None): axis[i]=d.select(ax,faxis=select.keys(),frange=select.values()) # Define default bin range - if wnorm is None : wnorm=[None]*len(axis) + if type(normed) is bool: + if normed: + normed=[True]*len(axis) + else: + normed=[False]*len(axis) + if brange is None : brange=[[None,None]]*len(axis) if bwidth is None : bwidth=[None]*len(axis) if blen is None : blen=[None]*len(axis) @@ -89,17 +94,19 @@ def hn(self,axis,blen=None,bwidth=None,brange=None,wnorm=None,select=None): bwidth[i] = (brange[i][1] - brange[i][0])/blen[i] else: bwidth[i] = (brange[i][1] - brange[i][0])/blen[i] - if wnorm[i] is None: wnorm[i]=bwidth[i] + + wnorm = 1. + if normed[i] is True: wnorm*=bwidth[i] # bins.append(np.linspace(brange[i][0],brange[i][1]+bwidth[i],blen[i])) bins.append(np.linspace(brange[i][0],brange[i][1],blen[i])) # Calculate the multi dimensional histo, normalized by wnorm - h,b=np.histogramdd(axis,weights=w/np.product(wnorm),bins=bins) + h,b=np.histogramdd(axis,weights=w/wnorm,bins=bins) # Return the bins and histo return b,h - def h1(self,axis,bwidth=None,brange=None,wnorm=None,select=None): + def h1(self,axis,bwidth=None,brange=None,normed=True,select=None): """ Create and return the 1 dimensional histogram of given axis. @@ -131,11 +138,11 @@ def h1(self,axis,bwidth=None,brange=None,wnorm=None,select=None): """ if not brange : brange = [None,None] - b,h=self.hn([axis],bwidth=[bwidth],brange=[brange],wnorm=[wnorm],select=select) + b,h=self.hn([axis],bwidth=[bwidth],brange=[brange],normed=normed,select=select) return b[0],h - def h2(self,axis1,axis2,bwidth1=None,bwidth2=None,brange1=None,brange2=None,wnorm=None,select=None): + def h2(self,axis1,axis2,bwidth1=None,bwidth2=None,brange1=None,brange2=None,normed=True,select=None): """ Create and return the 2 dimensional histogram of given axis. @@ -168,7 +175,7 @@ def h2(self,axis1,axis2,bwidth1=None,bwidth2=None,brange1=None,brange2=None,wnor if not brange1 : brange1 = [None,None] if not brange2 : brange2 = [None,None] - b,h=self.hn([axis1,axis2],bwidth=[bwidth1,bwidth2],brange=[brange1,brange2],wnorm=wnorm,select=select) + b,h=self.hn([axis1,axis2],bwidth=[bwidth1,bwidth2],brange=[brange1,brange2],normed=normed,select=select) return b[0],b[1],h diff --git a/p2sat/_Plot.py b/p2sat/_Plot.py index d3fe8f1..81c7d49 100644 --- a/p2sat/_Plot.py +++ b/p2sat/_Plot.py @@ -13,7 +13,7 @@ def __init__(self,PhaseSpace): self.cmap="viridis" plt.ion() - def get_labels(self,axes,wnorm): + def get_labels(self,axes,normed): """ Returns the labels of given axes. @@ -21,7 +21,7 @@ def get_labels(self,axes,wnorm): ---------- axes : list of str Names of the axes - wnorm : float or None + normed : float or None Weight normalization. If None, the last labels element is \"Number\", otherwise it is \"Number/unit1/unit2/...\" Returns @@ -35,7 +35,7 @@ def get_labels(self,axes,wnorm): N_name="" N_unit="" # Loop over all the axes - for ax in axes: + for i,ax in enumerate(axes): # No label if ax is not a str if type(ax) is not str: labels.append("") @@ -49,8 +49,8 @@ def get_labels(self,axes,wnorm): # Format the label for axis and unit of N if ax_unit is not None: labels.append("{} ({})".format(ax_name,ax_unit)) - if wnorm is None:N_name += "d%s "%(ax_name[1:-1]) # LaTeX without '$' - if wnorm is None:N_unit += "%s "%(ax_unit) + if normed[i]: N_name += "d%s "%(ax_name[1:-1]) # LaTeX without '$' + if normed[i]: N_unit += "%s "%(ax_unit) else: labels.append(ax_name) @@ -129,7 +129,7 @@ def h1(self,axis,where='post',log=False,polar=False,reverse=False,**kargs): else: a=plt.gca() - labels=self.get_labels([axis],kargs.get('wnorm',None)) + labels=self.get_labels([axis],kargs.get('normed',[True])) b,h=self._h.h1(axis,**kargs) @@ -197,7 +197,7 @@ def f1(self,axis,func_name,log=False,polar=False,reverse=False,**kargs): else: a=plt.gca() - labels=self.get_labels([axis],kargs.get('wnorm',None)) + labels=self.get_labels([axis],kargs.get('normed',[True,True])) b,h=self._h.f1(axis,func_name,return_fit=True,**kargs) @@ -252,7 +252,7 @@ def h2(self,axis1,axis2,log=False,polar=False,**kargs): a=plt.gca(polar=True) else: a=plt.gca() - labels=self.get_labels([axis1,axis2],kargs.get('wnorm',None)) + labels=self.get_labels([axis1,axis2],kargs.get('normed',[True,True])) b1,b2,h=self._h.h2(axis1,axis2,**kargs) g1,g2=np.meshgrid(b1,b2,indexing='ij') @@ -304,7 +304,7 @@ def c2(self,axis1,axis2,log=False,polar=False,gfilter=0.0,**kargs): a=plt.gca(polar=True) else: a=plt.gca() - labels=self.get_labels([axis1,axis2],kargs.get('wnorm',None)) + labels=self.get_labels([axis1,axis2],kargs.get('normed',[True,True])) if log: from matplotlib.colors import LogNorm @@ -356,7 +356,7 @@ def s2(self,axis1,axis2,log=False,polar=False,select=None): a=plt.gca(polar=True) else: a=plt.gca() - labels=self.get_labels([axis1,axis2],wnorm=1) + labels=self.get_labels([axis1,axis2],normed=[True,True]) d = self._ps.data From d378ad00ffe7fe8176319314955d0396db3d30c9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o?= Date: Mon, 10 Sep 2018 15:53:47 +0200 Subject: [PATCH 12/19] Comment and simplify hn --- p2sat/_Hist.py | 63 +++++++++++++++++++++++++------------------------- 1 file changed, 32 insertions(+), 31 deletions(-) diff --git a/p2sat/_Hist.py b/p2sat/_Hist.py index 9ef4e09..833d4f8 100644 --- a/p2sat/_Hist.py +++ b/p2sat/_Hist.py @@ -8,7 +8,7 @@ class _Hist(object): def __init__(self,PhaseSpace): self._ps=PhaseSpace - def hn(self,axis,blen=None,bwidth=None,brange=None,normed=True,select=None): + def hn(self,axis,bwidth=None,brange=None,normed=True,select=None): """ Create and return the n-dimensional histo of axis list. @@ -16,14 +16,12 @@ def hn(self,axis,blen=None,bwidth=None,brange=None,normed=True,select=None): ---------- axis : list of str/np.array list of axis to hist - blen : list of int, optional - list of number of bins. If a blen element is None, the default value is 10 bwidth : list of float, optional list of bin width. If a bwidth element is None, a calculation is done to have 10 bins in the correspondant axis brange : list of list of 2 float, optional list of bin minimum and maximum. If a brange element is None, the minimum/maximum of the axis is taken - normed : list of bool, optional - weight normalization. If a normed element is True, the bin width is taken + normed : bool or list of bool, optional + weight normalization. If a normed element is True, the bin width is taken as weight normalization select : dict, optional filtering dictionary @@ -44,18 +42,19 @@ def hn(self,axis,blen=None,bwidth=None,brange=None,normed=True,select=None): Examples -------- - >>> hn(['x'],bwidth=[50],brange=[[0,1000]],normed=[True],select={'ekin':(0.511,None)}) + >>> hn(['x'],bwidth=[50],brange=[[0,1000]],normed=[False],select={'ekin':(0.511,None)}) returns the number of particles with :math:`ekin \in [0.511, +\infty] MeV` in function of x - normed=[True] to not divide nb of particles by bin width (otherwise number per um) + normed=[False] to not divide nb of particles by bin width (otherwise number per um) >>> hn(['r','ekin'],bwidth=[10.0,0.1],brange=[[0,1000],[0.1,50.0]],select={'x':150}) returns a number of e- per um per MeV at x=150 um """ + # Get a shortcut to data object d=self._ps.data - # Get a copy the axis from a str if needed + # Get a copy of the axes for i,ax in enumerate(axis): if type(ax) is str:ax = eval("d.%s"%ax) axis[i]=np.array(ax) @@ -69,37 +68,39 @@ def hn(self,axis,blen=None,bwidth=None,brange=None,normed=True,select=None): for i,ax in enumerate(axis): axis[i]=d.select(ax,faxis=select.keys(),frange=select.values()) - # Define default bin range - if type(normed) is bool: - if normed: - normed=[True]*len(axis) - else: - normed=[False]*len(axis) - + # Define bin range if brange is None : brange=[[None,None]]*len(axis) - if bwidth is None : bwidth=[None]*len(axis) - if blen is None : blen=[None]*len(axis) - bins=[] for i,_ in enumerate(axis): if brange[i][0] is None:brange[i][0]=min(axis[i]) if brange[i][1] is None:brange[i][1]=max(axis[i]) - # if brange[i][0]==brange[i][1]: - # brange[i][1]=2 - # blen[i]=2 - if blen[i] is None: - if bwidth[i] is not None: - blen[i] = int(np.ceil((brange[i][1] + bwidth[i] - brange[i][0])/bwidth[i])) - else: - blen[i] = 10 - bwidth[i] = (brange[i][1] - brange[i][0])/blen[i] + + # Define bin width + if bwidth is None : bwidth=[None]*len(axis) + blen=[None]*len(axis) + for i,_ in enumerate(axis): + if bwidth[i] is not None: + # Bin width is over-estimated to fit with bin range + blen[i] = int(np.ceil((brange[i][1] + bwidth[i] - brange[i][0])/bwidth[i])) else: - bwidth[i] = (brange[i][1] - brange[i][0])/blen[i] + # Default is 10 bins + blen[i] = 10 + bwidth[i] = float(brange[i][1] - brange[i][0])/blen[i] - wnorm = 1. - if normed[i] is True: wnorm*=bwidth[i] - # bins.append(np.linspace(brange[i][0],brange[i][1]+bwidth[i],blen[i])) + # Calculate bins + bins=[] + for i,_ in enumerate(axis): bins.append(np.linspace(brange[i][0],brange[i][1],blen[i])) + # Get weight normalization + if type(normed) is bool: + if normed: + normed=[True]*len(axis) + else: + normed=[False]*len(axis) + wnorm = 1. + for i,_ in enumerate(axis): + if normed[i]: wnorm*=bwidth[i] + # Calculate the multi dimensional histo, normalized by wnorm h,b=np.histogramdd(axis,weights=w/wnorm,bins=bins) From 224229112eebfd55b1d18c2bdf000b0c6d20d19a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o?= Date: Mon, 10 Sep 2018 15:58:15 +0200 Subject: [PATCH 13/19] Fix: normed in plot --- p2sat/_Plot.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/p2sat/_Plot.py b/p2sat/_Plot.py index 81c7d49..07a7348 100644 --- a/p2sat/_Plot.py +++ b/p2sat/_Plot.py @@ -34,6 +34,11 @@ def get_labels(self,axes,normed): labels=[] N_name="" N_unit="" + if type(normed) is bool: + if normed: + normed=[True]*len(axes) + else: + normed=[False]*len(axes) # Loop over all the axes for i,ax in enumerate(axes): # No label if ax is not a str @@ -129,7 +134,7 @@ def h1(self,axis,where='post',log=False,polar=False,reverse=False,**kargs): else: a=plt.gca() - labels=self.get_labels([axis],kargs.get('normed',[True])) + labels=self.get_labels([axis],kargs.get('normed',True)) b,h=self._h.h1(axis,**kargs) @@ -197,7 +202,7 @@ def f1(self,axis,func_name,log=False,polar=False,reverse=False,**kargs): else: a=plt.gca() - labels=self.get_labels([axis],kargs.get('normed',[True,True])) + labels=self.get_labels([axis],kargs.get('normed',True)) b,h=self._h.f1(axis,func_name,return_fit=True,**kargs) @@ -252,7 +257,7 @@ def h2(self,axis1,axis2,log=False,polar=False,**kargs): a=plt.gca(polar=True) else: a=plt.gca() - labels=self.get_labels([axis1,axis2],kargs.get('normed',[True,True])) + labels=self.get_labels([axis1,axis2],kargs.get('normed',True)) b1,b2,h=self._h.h2(axis1,axis2,**kargs) g1,g2=np.meshgrid(b1,b2,indexing='ij') @@ -304,7 +309,7 @@ def c2(self,axis1,axis2,log=False,polar=False,gfilter=0.0,**kargs): a=plt.gca(polar=True) else: a=plt.gca() - labels=self.get_labels([axis1,axis2],kargs.get('normed',[True,True])) + labels=self.get_labels([axis1,axis2],kargs.get('normed',True)) if log: from matplotlib.colors import LogNorm @@ -356,7 +361,7 @@ def s2(self,axis1,axis2,log=False,polar=False,select=None): a=plt.gca(polar=True) else: a=plt.gca() - labels=self.get_labels([axis1,axis2],normed=[True,True]) + labels=self.get_labels([axis1,axis2],normed=True) d = self._ps.data From ae999b868356fb44be2d2495c7612c43ed774014 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o?= Date: Mon, 10 Sep 2018 16:37:43 +0200 Subject: [PATCH 14/19] Add beta and v --- p2sat/_Data.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/p2sat/_Data.py b/p2sat/_Data.py index 2c124f9..ffac069 100644 --- a/p2sat/_Data.py +++ b/p2sat/_Data.py @@ -71,6 +71,8 @@ def __init__(self,PhaseSpace): 'p' : '$p$', 'ekin' : '$E_{kin}$', 'gamma' : '$\gamma$', + 'beta' : '$\\beta$', + 'v' : '$v$', 'theta' : '$\\theta$', 'phi' : '$\\phi$', @@ -129,12 +131,17 @@ def update(self,w,x,y,z,px,py,pz,t,verbose=True): self.theta = np.degrees(np.arctan(self.py/self.px)) self.phi = np.degrees(np.arctan(self.pz/self.px)) # Geometrical effect ? change -> 0 pi mass = self._ps.mass + c = 2.99792458e8 * 1e6/1e15 # speed of light in um/fs if mass == 0: self.ekin = self.p self.gamma = np.inf + self.beta = 1 + self.v = c else: self.ekin = (np.sqrt((self.p/mass)**2 + 1) - 1) * mass self.gamma = self.ekin/mass + 1. + self.beta = np.sqrt(1.-1/self.gamma**2) + self.v = self.beta * c if verbose: print("Done !") def generate(self,Nconf,Npart,ekin,theta,phi,pos=None,time=None,verbose=True): From 696e260a000e6839bb349dcded35cc99902bdccb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o?= Date: Mon, 10 Sep 2018 17:11:51 +0200 Subject: [PATCH 15/19] Add propagate --- p2sat/_Data.py | 50 ++++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 44 insertions(+), 6 deletions(-) diff --git a/p2sat/_Data.py b/p2sat/_Data.py index ffac069..988744c 100644 --- a/p2sat/_Data.py +++ b/p2sat/_Data.py @@ -96,6 +96,8 @@ def __init__(self,PhaseSpace): 'p' : 'MeV/c', 'ekin' : 'MeV', 'gamma' : None, + 'beta' : None, + 'v' : 'um/fs', 'theta' : 'deg', 'phi' : 'deg' @@ -134,9 +136,9 @@ def update(self,w,x,y,z,px,py,pz,t,verbose=True): c = 2.99792458e8 * 1e6/1e15 # speed of light in um/fs if mass == 0: self.ekin = self.p - self.gamma = np.inf - self.beta = 1 - self.v = c + self.gamma = np.array([np.inf]*len(w)) + self.beta = np.array([1]*len(w)) + self.v = np.array([c]*len(w)) else: self.ekin = (np.sqrt((self.p/mass)**2 + 1) - 1) * mass self.gamma = self.ekin/mass + 1. @@ -335,13 +337,49 @@ def transformate(self,translation=(0.0,0.0,0.0),rotation=(0.0,0.0)): """ pass - def propagate(self): + def propagate(self,x_pos=None,time=None,verbose=True): """ Propagate the phase space to a given position or time. - TODO + Parameters + ---------- + x_pos : float, optional + propagate the phase-space untill x = x_pos. Default is None (no propagation) + time : float, optional + propagate the phase-space untill t = time. Default is None (no propagation) + + Notes + ----- + x_pos and time can not be defined simultaneously. """ - pass + if time is None and x_pos is None: + raise ValueError("You must define time or x_pos.") + if time is not None and x_pos is not None: + raise ValueError("time and x_pos can not be defined simultaneously.") + + w = self.w + px = self.px + py = self.py + pz = self.pz + + if time is not None: + if verbose: print("Propagate %s phase-space to time = %.4E fs."%(self._ps.specie,time)) + t = np.array([time]*len(w)) + Dt = t - self.t + x = self.x + (self.px/self.p)*self.v*Dt + y = self.y + (self.py/self.p)*self.v*Dt + z = self.z + (self.pz/self.p)*self.v*Dt + + if x_pos is not None: + if verbose: print("Propagate %s phase-space to x = %.4E um."%(self._ps.specie,x_pos)) + x = np.array([x_pos]*len(w)) + Dt = (x - self.x)/self.v + t = self.t + Dt + y = self.y + (self.py/self.p)*self.v*Dt + z = self.z + (self.pz/self.p)*self.v*Dt + + self.update(w,x,y,z,px,py,pz,t,verbose=verbose) + def lorrentz(self,beta_CM,pos_CM=None): """ From 82f3e5d63a482066151ffca4980aee16acd4df6e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o?= Date: Mon, 10 Sep 2018 17:21:21 +0200 Subject: [PATCH 16/19] Update doc & PhaseSpace.__str__ --- p2sat/PhaseSpace.py | 11 ++++++++++- p2sat/_Data.py | 44 ++++++-------------------------------------- 2 files changed, 16 insertions(+), 39 deletions(-) diff --git a/p2sat/PhaseSpace.py b/p2sat/PhaseSpace.py index da6cdc6..5f2fb01 100644 --- a/p2sat/PhaseSpace.py +++ b/p2sat/PhaseSpace.py @@ -82,21 +82,30 @@ def __add__(self,other): return ps def __str__(self): - txt = "" + txt = "\n" txt += "p2sat PhaseSpace instance located at %s\n\n"%hex(id(self)) txt += "Specie : %s\n"%self.specie txt += "Number of configurations : %i\n"%len(self) txt += "Total number of particles : %.4E\n\n"%sum(self.data.w) txt += "Statistics : ( min , max , mean , std ) unit\n" + txt += " w : (% .3E, % .3E, % .3E, % .3E) \n"%(min(self.data.w),max(self.data.w),self.data.w.mean(),self.data.w.std()) + txt += " t : (% .3E, % .3E, % .3E, % .3E) fs\n\n"%(min(self.data.t),max(self.data.t),self.data.t.mean(),self.data.t.std()) + txt += " x : (% .3E, % .3E, % .3E, % .3E) um\n"%(min(self.data.x),max(self.data.x),self.data.x.mean(),self.data.x.std()) txt += " y : (% .3E, % .3E, % .3E, % .3E) um\n"%(min(self.data.y),max(self.data.y),self.data.y.mean(),self.data.y.std()) txt += " z : (% .3E, % .3E, % .3E, % .3E) um\n"%(min(self.data.z),max(self.data.z),self.data.z.mean(),self.data.z.std()) txt += " r : (% .3E, % .3E, % .3E, % .3E) um\n\n"%(min(self.data.r),max(self.data.r),self.data.r.mean(),self.data.r.std()) + txt += " px : (% .3E, % .3E, % .3E, % .3E) MeV/c\n"%(min(self.data.px),max(self.data.px),self.data.px.mean(),self.data.px.std()) txt += " py : (% .3E, % .3E, % .3E, % .3E) MeV/c\n"%(min(self.data.py),max(self.data.py),self.data.py.mean(),self.data.py.std()) txt += " pz : (% .3E, % .3E, % .3E, % .3E) MeV/c\n"%(min(self.data.pz),max(self.data.pz),self.data.pz.mean(),self.data.pz.std()) txt += " p : (% .3E, % .3E, % .3E, % .3E) MeV/c\n\n"%(min(self.data.p),max(self.data.p),self.data.p.mean(),self.data.p.std()) + txt += " ekin : (% .3E, % .3E, % .3E, % .3E) MeV\n"%(min(self.data.ekin),max(self.data.ekin),self.data.ekin.mean(),self.data.ekin.std()) + txt += " gamma : (% .3E, % .3E, % .3E, % .3E) \n"%(min(self.data.gamma),max(self.data.gamma),self.data.gamma.mean(),self.data.gamma.std()) + txt += " beta : (% .3E, % .3E, % .3E, % .3E) \n"%(min(self.data.beta),max(self.data.beta),self.data.beta.mean(),self.data.beta.std()) + txt += " v : (% .3E, % .3E, % .3E, % .3E) um/fs\n\n"%(min(self.data.v),max(self.data.v),self.data.v.mean(),self.data.v.std()) + txt += " theta : (% .3E, % .3E, % .3E, % .3E) deg\n"%(min(self.data.theta),max(self.data.theta),self.data.theta.mean(),self.data.theta.std()) txt += " phi : (% .3E, % .3E, % .3E, % .3E) deg\n"%(min(self.data.phi),max(self.data.phi),self.data.phi.mean(),self.data.phi.std()) txt += "" diff --git a/p2sat/_Data.py b/p2sat/_Data.py index 988744c..1a2508c 100644 --- a/p2sat/_Data.py +++ b/p2sat/_Data.py @@ -206,8 +206,8 @@ def generate(self,Nconf,Npart,ekin,theta,phi,pos=None,time=None,verbose=True): """ # Print a starting message if verbose: - print("Generate particle phase-space for \"%s\" ekin law, \"%s\" theta law, \"%s\" phi law ..." - %(ekin["law"],theta["law"],phi["law"])) + print("Generate %s phase-space for \"%s\" ekin law, \"%s\" theta law, \"%s\" phi law ..." + %(self._ps.specie,ekin["law"],theta["law"],phi["law"])) # Ensure that Nconf is of type int (for values such as 1e6) Nconf = int(Nconf) @@ -381,21 +381,18 @@ def propagate(self,x_pos=None,time=None,verbose=True): self.update(w,x,y,z,px,py,pz,t,verbose=verbose) - def lorrentz(self,beta_CM,pos_CM=None): + def lorentz(self,beta_CM,verbose=True): """ - Lorrentz-transformate the particle phase-space with given speed of the center of mass. + Lorentz-transformate the particle phase-space with given speed of the center of mass. Notes ----- https://en.wikipedia.org/wiki/Lorentz_transformation#Transformation_of_other_quantities """ - # If no position given, default is the origin - if pos_CM is None: - pos_CM = [0.,0.,0.] - + if verbose: print("Lorentz-transform %s phase-space with center of mass frame moving at %s c."%(self._ps.specie,beta_CM)) # lowercase : scalar, caption : vector - B = np.array(beta_CM) + B = -np.array(beta_CM) b = np.dot(B,B) N = B/b c = 2.99792458e8 * 1e6/1e15 # speed of light in um/fs @@ -427,35 +424,6 @@ def lorrentz(self,beta_CM,pos_CM=None): w = self.w self.update(w,x,y,z,px,py,pz,t) - # - # beta_CM = np.array(beta_CM) - # b = np.dot(beta_CM,beta_CM) - # n = beta_CM/b - # c = 2.99792458e8 * 1e6/1e15 # speed of CM in um/fs - # v = b * c - # gamma = 1./np.sqrt(1-beta_CM**2) - # g = 1./np.sqrt(1-b**2) - # G = np.matrix([ - # [g ,-g*b*n[0] ,-g*b*n[1] ,-g*b*n[2] ], - # [-g*b*n[0] ,1+(g-1)*n[0]**2 ,(g-1)*n[0]*n[1] ,(g-1)*n[0]*n[2]], - # [-g*b*n[1] ,(g-1)*n[1]*n[0] ,1+(g-1)*n[1]**2 ,(g-1)*n[1]*n[2]], - # [-g*b*n[2] ,(g-1)*n[2]*n[0] ,(g-1)*n[2]*n[1] ,1+(g-1)*n[2]**2] - # ]) - # - # # # Time transformation - # q1 = [c*self.t.copy(),self.x.copy(),self.y.copy(),self.z.copy()] - # p1 = [self.ekin.copy()/c,self.px.copy(),self.py.copy(),self.pz.copy()] - # ct,x,y,z = np.dot(G,q1) - # ekin_c,px,py,pz = np.dot(G,p1) - # w = self.w - # t = np.array(ct)[0]/c - # x = np.array(x)[0] - # y = np.array(y)[0] - # z = np.array(z)[0] - # px = np.array(px)[0] - # py = np.array(py)[0] - # pz = np.array(pz)[0] - # self.update(w,x,y,z,px,py,pz,t) def discretize(self,with_time=True,verbose=True,**kargs): """ From 60c409f8c4892faf3a8017025aff7b76c9b0dbb7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o?= Date: Mon, 10 Sep 2018 17:37:20 +0200 Subject: [PATCH 17/19] Minor synthax changes --- p2sat/_Data.py | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/p2sat/_Data.py b/p2sat/_Data.py index 1a2508c..be5f959 100644 --- a/p2sat/_Data.py +++ b/p2sat/_Data.py @@ -135,15 +135,15 @@ def update(self,w,x,y,z,px,py,pz,t,verbose=True): mass = self._ps.mass c = 2.99792458e8 * 1e6/1e15 # speed of light in um/fs if mass == 0: - self.ekin = self.p - self.gamma = np.array([np.inf]*len(w)) - self.beta = np.array([1]*len(w)) - self.v = np.array([c]*len(w)) + self.ekin = self.p + self.gamma = np.array([np.inf]*len(w)) + self.beta = np.array([1.]*len(w)) + self.v = np.array([c]*len(w)) else: - self.ekin = (np.sqrt((self.p/mass)**2 + 1) - 1) * mass - self.gamma = self.ekin/mass + 1. - self.beta = np.sqrt(1.-1/self.gamma**2) - self.v = self.beta * c + self.ekin = (np.sqrt((self.p/mass)**2 + 1) - 1) * mass + self.gamma = self.ekin/mass + 1. + self.beta = np.sqrt(1.-1/self.gamma**2) + self.v = self.beta * c if verbose: print("Done !") def generate(self,Nconf,Npart,ekin,theta,phi,pos=None,time=None,verbose=True): @@ -234,7 +234,7 @@ def generate(self,Nconf,Npart,ekin,theta,phi,pos=None,time=None,verbose=True): mangle = phi["max"]*np.pi/180. except KeyError: mangle = np.pi - phi0 = np.random.uniform(0,mangle,Nconf) + phi0 = np.random.uniform(0.,mangle,Nconf) elif phi["law"]=="gauss": phi0 = np.random.normal(phi["mu"],phi["sigma"],Nconf) # Generate energy @@ -380,7 +380,6 @@ def propagate(self,x_pos=None,time=None,verbose=True): self.update(w,x,y,z,px,py,pz,t,verbose=verbose) - def lorentz(self,beta_CM,verbose=True): """ Lorentz-transformate the particle phase-space with given speed of the center of mass. From d060895817e7636cff71856cab03ba5c734bae16 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o?= Date: Thu, 13 Sep 2018 17:01:33 +0200 Subject: [PATCH 18/19] Fix: number of particles in hist.f1 --- p2sat/_Hist.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/p2sat/_Hist.py b/p2sat/_Hist.py index 833d4f8..dd27d7b 100644 --- a/p2sat/_Hist.py +++ b/p2sat/_Hist.py @@ -253,15 +253,17 @@ def f1(self,axis,func_name,return_fit=False,verbose=True,**kargs): and returns fit parameters A,sigma,mu. """ # Get the hist data + if verbose:print("Processing 1D \"{}\" fit of \"{}\" ...".format(func_name,str(axis))) x,w = self._ps.hist.h1(axis,**kargs) + Ntot = np.trapz(w,x[:-1]) # Define fit function and default values for fit parameters if func_name=="exp": f = lambda x,A,T: A*np.exp(-x/T)/T - p0 = [sum(w),1] + p0 = [Ntot,1] elif func_name=="gauss": f = lambda x,A,sigma,mu: A/(np.sqrt(2*np.pi) * sigma) * np.exp(-(x-mu)**2/(2*sigma**2)) - p0 = [sum(w),x.std(),0] + p0 = [Ntot,x.std(),0] else: raise NameError("Unknown func_name.") @@ -271,9 +273,9 @@ def f1(self,axis,func_name,return_fit=False,verbose=True,**kargs): # Print error on number of particles if verbose: - print('Parameters are {}'.format(popt)) - diff = (popt[0]-sum(w))/sum(w) * 100 - print('Error on number of particles for \"{}\" fit : {:.2F} %'.format(func_name,diff)) + print('Parameters : {}'.format(popt)) + diff = (popt[0]-Ntot)/Ntot * 100 + print('Error on number of particles : {:.2F} %'.format(diff)) # Choice of the return values if return_fit: From f66ca562093cee23e878684e92ff93a24ff4c073 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o?= Date: Thu, 13 Sep 2018 17:06:55 +0200 Subject: [PATCH 19/19] Add example to generate phase space --- examples/generate.py | 42 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) create mode 100644 examples/generate.py diff --git a/examples/generate.py b/examples/generate.py new file mode 100644 index 0000000..16dc2b1 --- /dev/null +++ b/examples/generate.py @@ -0,0 +1,42 @@ +#coding:utf8 + +""" +This is an example of how to use the `PhaseSpace.stat` object of p2sat. + +It allows to make statistics on phase space data in a very simple way +""" + +# Import p2sat +p2sat_path="../" +import sys +if p2sat_path not in sys.path:sys.path.append(p2sat_path) +import p2sat + +# Boolean to export or not the generated phase space +export = False + +# Instanciate a PhaseSpace object for electron specie +eps = p2sat.PhaseSpace(specie="electron") + +# Define energy and angle parameters +ekin_dict = {"law":"exp","T":1.0} +# theta_dict = {"law":"gauss","mu":0.,"sigma":5.} +theta_dict = {"law":"iso","max":20.} +phi_dict = {"law":"iso"} + +# Generate particle phase space +eps.data.generate(Nconf = 1e5, Npart = 1e12, ekin=ekin_dict, theta=theta_dict, phi=phi_dict) + +# Look at the consistency of phase space generation +print(eps) +eps.plot.figure(0) +eps.plot.h1('ekin',bwidth=.1,log=True) +eps.plot.f1('ekin',func_name="exp",log=True,bwidth=.1) + +eps.plot.figure(1) +eps.plot.h2('theta','phi',log=True, + bwidth1=.5,bwidth2=1.) + +# Export phase space if needed +if export: + eps.export.txt("test_eps.csv",sep=",")