Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-04-25 02:43:25 

0001 #! /usr/bin/env python
0002 
0003 import sys
0004 import os
0005 import math
0006 import random
0007 from scipy import optimize
0008 from scipy import interpolate
0009 import numpy
0010 import copy
0011 
0012 Q2max=1.0 # 1 GeV^2 as the maximally allowed Q2
0013 ion_Form=1 # Form1: Q**2=kT**2+(mn*x)**2, Qmin**2=(mn*x)**2; 
0014            # Form2: Q**2=kT**2/(1-x)+(mn*x)**2/(1-x), Qmin**2=(mn*x)**2/(1-x)
0015 
0016 files=[arg for arg in sys.argv[1:] if arg.startswith('--file=')]
0017 nuclei=[arg for arg in sys.argv[1:] if arg.startswith('--beams=')]
0018 if not files or not nuclei:
0019     raise Exception( "The usage of it should be e.g., ./lhe_ktsmearing_UPC --beams='Pb208 Pb208' --file='/PATH/TO/file.lhe' --out='ktsmearing.lhe4upc' ")
0020 files=files[0]
0021 files=files.replace('--file=','')
0022 #files=[file.lower() for file in files.split(' ')]
0023 files=[file for file in files.split(' ')]
0024 files=[files[0]]
0025 nuclei=nuclei[0]
0026 nuclei=nuclei.replace('--beams=','')
0027 nuclei=[nucleus.rstrip().lstrip() for nucleus in nuclei.split(' ')]
0028 
0029 # name:(RA,aA,wA), RA and aA are in fm, need divide by GeVm12fm to get GeV-1
0030 GeVm12fm=0.1973
0031 WoodsSaxon={'H2':(0.01,0.5882,0),'Li7':(1.77,0.327,0),'Be9':(1.791,0.611,0),'B10':(1.71,0.837,0),'B11':(1.69,0.811,0),\
0032                 'C13':(1.635,1.403,0),'C14':(1.73,1.38,0),'N14':(2.570,0.5052,-0.180),'N15':(2.334,0.498,0.139),'O16':(2.608,0.513,-0.051),'Ne20':(2.791,0.698,-0.168),\
0033                 'Mg24':(3.108,0.607,-0.163),'Mg25':(3.22,0.58,-0.236),'Al27':(3.07,0.519,0),'Si28':(3.340,0.580,-0.233),'Si29':(3.338,0.547,-0.203),'Si30':(3.338,0.547,-0.203),\
0034                 'P31':(3.369,0.582,-0.173),'Cl35':(3.476,0.599,-0.10),'Cl37':(3.554,0.588,-0.13),'Ar40':(3.766,0.586,-0.161),'K39':(3.743,0.595,-0.201),'Ca40':(3.766,0.586,-0.161),\
0035                 'Ca48':(3.7369,0.5245,-0.030),'Ni58':(4.3092,0.5169,-0.1308),'Ni60':(4.4891,0.5369,-0.2668),'Ni61':(4.4024,0.5401,-0.1983),'Ni62':(4.4425,0.5386,-0.2090),'Ni64':(4.5211,0.5278,-0.2284),\
0036                 'Cu63':(4.214,0.586,0),'Kr78':(4.5,0.5,0),'Ag110':(5.33,0.535,0),'Sb122':(5.32,0.57,0),'Xe129':(5.36,0.59,0),'Xe132':(5.4,0.61,0),\
0037                 'Nd142':(5.6135,0.5868,0.096),'Er166':(5.98,0.446,0.19),'W186':(6.58,0.480,0),'Au197':(6.38,0.535,0),'Pb207':(6.62,0.546,0),'Pb208':(6.624,0.549,0)}
0038 
0039 if nuclei[0] != 'p' and nuclei[0] not in WoodsSaxon.keys():
0040     raise ValueError('do not know the first beam type = %s'%nuclei[0])
0041 
0042 if nuclei[1] != 'p' and nuclei[1] not in WoodsSaxon.keys():
0043     raise ValueError('do not know the second beam type = %s'%nuclei[1])
0044 
0045 outfile=[arg for arg in sys.argv[1:] if arg.startswith('--out=')]
0046 if not outfile:
0047     outfile=['ktsmearing.lhe4upc']
0048 
0049 outfile=outfile[0]
0050 outfile=outfile.replace('--out=','')
0051 
0052 currentdir=os.getcwd()
0053 
0054 p_Q2max_save=1
0055 p_x_array=None # an array of log10(1/x)
0056 p_xmax_array=None # an array of maximal function value at logQ2/Q02, where Q02=0.71
0057 p_fmax_array=None # an array of maximal function value
0058 p_xmax_interp=None
0059 p_fmax_interp=None
0060 
0061 offset=100
0062 
0063 def generate_Q2_epa_proton(x,Q2max):
0064     if x >= 1.0 or x <= 0:
0065         raise ValueError( "x >= 1 or x <= 0")
0066     mp=0.938272081 # proton mass in unit of GeV
0067     mupomuN=2.793
0068     Q02=0.71  # in unit of GeV**2
0069     mp2=mp**2
0070     Q2min=mp2*x**2/(1-x)
0071 
0072     def xmaxvalue(Q2MAX):
0073         val=(math.sqrt(Q2MAX*(4*mp2+Q2MAX))-Q2MAX)/(2*mp2)
0074         return val
0075 
0076     global p_x_array
0077     global p_Q2max_save
0078     global p_xmax_array
0079     global p_fmax_array
0080     global p_xmax_interp
0081     global p_fmax_interp
0082 
0083     if Q2max <= Q2min or x >= xmaxvalue(Q2max) : return Q2max
0084 
0085     logQ2oQ02max = math.log(Q2max/Q02)
0086     logQ2oQ02min = math.log(Q2min/Q02)
0087 
0088     def distfun(xx,logQ2oQ02):
0089         exp=math.exp(logQ2oQ02)
0090         funvalue=(-8*mp2**2*xx**2+exp**2*mupomuN**2*Q02**2*\
0091                        (2-2*xx+xx**2)+2*exp*mp2*Q02*(4-4*xx+mupomuN**2*xx**2))\
0092                        /(2*exp*(1+exp)**4*Q02*(4*mp2+exp*Q02))
0093         return funvalue
0094 
0095     if p_x_array is None or (p_Q2max_save != Q2max):
0096         # we need to generate the grid first
0097         p_Q2max_save = Q2max
0098         xmaxQ2max=xmaxvalue(Q2max)
0099         log10xmaxQ2maxm1=math.log10(1/xmaxQ2max)
0100         p_x_array=[]
0101         p_xmax_array=[]
0102         p_fmax_array=[]
0103         for log10xm1 in range(10):
0104             for j in range(10):
0105                 tlog10xm1=log10xmaxQ2maxm1+0.1*j+log10xm1
0106                 p_x_array.append(tlog10xm1)
0107                 xx=10**(-tlog10xm1)
0108                 if log10xm1 == 0 and j == 0:
0109                     max_Q2 = logQ2oQ02max
0110                     max_fun = distfun(xx,max_Q2)
0111                     p_xmax_array.append(max_Q2)
0112                     p_fmax_array.append(max_fun)
0113                 else:
0114                     max_Q2 = optimize.fmin(lambda x0: -distfun(xx,x0),\
0115                                                     (logQ2oQ02max+logQ2oQ02min)/2,\
0116                                                full_output=False,disp=False)
0117                     max_fun = distfun(xx,max_Q2[0])
0118                     p_xmax_array.append(max_Q2[0])
0119                     p_fmax_array.append(max_fun)
0120         p_x_array=numpy.array(p_x_array)
0121         p_xmax_array=numpy.array(p_xmax_array)
0122         p_fmax_array=numpy.array(p_fmax_array)
0123         p_xmax_interp=interpolate.interp1d(p_x_array,p_xmax_array)
0124         p_fmax_interp=interpolate.interp1d(p_x_array,p_fmax_array)
0125     log10xm1=math.log10(1/x)
0126     max_x = p_xmax_interp(log10xm1)
0127     max_fun = p_fmax_interp(log10xm1)
0128     logQ2oQ02now=logQ2oQ02min
0129     while True:
0130         r1=random.random() # a random float number between 0 and 1
0131         logQ2oQ02now=(logQ2oQ02max-logQ2oQ02min)*r1+logQ2oQ02min
0132         w=distfun(x,logQ2oQ02now)/max_fun
0133         r2=random.random() # a random float number between 0 and 1
0134         if r2 <= w: break
0135     Q2v=math.exp(logQ2oQ02now)*Q02
0136     return Q2v
0137 
0138 A_Q2max_save=[1,1]
0139 A_x_array=[None,None]  # an array of log10(1/x)
0140 A_xmax_array=[None,None] # an array of maximal function value at logQ2/Q02, where Q02=0.71
0141 A_fmax_array=[None,None] # an array of maximal function value
0142 A_xmax_interp=[None,None]
0143 A_fmax_interp=[None,None]
0144 
0145 # first beam: ibeam=0; second beam: ibeam=1
0146 def generate_Q2_epa_ion(ibeam,x,Q2max,RA,aA,wA):
0147     if x >= 1.0 or x <= 0:
0148         raise ValueError( "x >= 1 or x <= 0")
0149     if ibeam not in [0,1]:
0150         raise ValueError( "ibeam != 0,1")
0151     mn=0.9315 # averaged nucleon mass in unit of GeV
0152     Q02=0.71
0153     mn2=mn**2
0154     if ion_Form == 2:
0155         Q2min=mn2*x**2/(1-x)
0156     else:
0157         Q2min=mn2*x**2
0158     RAA=RA/GeVm12fm # from fm to GeV-1
0159     aAA=aA/GeVm12fm # from fm to GeV-1
0160     
0161     
0162     def xmaxvalue(Q2MAX):
0163         val=(math.sqrt(Q2MAX*(4*mn2+Q2MAX))-Q2MAX)/(2*mn2)
0164         return val
0165 
0166     global A_x_array
0167     global A_Q2max_save
0168     global A_xmax_array
0169     global A_fmax_array
0170     global A_xmax_interp
0171     global A_fmax_interp
0172 
0173     if Q2max <= Q2min or x >= xmaxvalue(Q2max) : return Q2max
0174 
0175     logQ2oQ02max = math.log(Q2max/Q02)
0176     logQ2oQ02min = math.log(Q2min/Q02)
0177 
0178     # set rhoA0=1 (irrelvant for this global factor)
0179     def FchA1(q):
0180         piqaA=math.pi*q*aAA
0181         funval=4*math.pi**4*aAA**3/(piqaA**2*math.sinh(piqaA)**2)*\
0182             (piqaA*math.cosh(piqaA)*math.sin(q*RAA)*(1-wA*aAA**2/RAA**2*\
0183             (6*math.pi**2/math.sinh(piqaA)**2+math.pi**2-3*RAA**2/aAA**2))\
0184             -q*RAA*math.sinh(piqaA)*math.cos(q*RAA)*(1-wA*aAA**2/RAA**2*\
0185             (6*math.pi**2/math.sinh(piqaA)**2+3*math.pi**2-RAA**2/aAA**2)))
0186         return funval
0187 
0188     # set rhoA0=1 (irrelvant for this global factor
0189     def FchA2(q):
0190         funval=0
0191         # only keep the first two terms
0192         for n in range(1,3):
0193             funval=funval+(-1)**(n-1)*n*math.exp(-n*RAA/aAA)/(n**2+q**2*aAA**2)**2*\
0194                 (1+12*wA*aAA**2/RAA**2*(n**2-q**2*aAA**2)/(n**2+q**2*aAA**2)**2)
0195         funval=funval*8*math.pi*aAA**3
0196         return funval
0197 
0198     def distfun(xx,logQ2oQ02):
0199         exp=math.exp(logQ2oQ02)*Q02
0200         if ion_Form == 2:
0201             FchA=FchA1(math.sqrt((1-xx)*exp))+FchA2(math.sqrt((1-xx)*exp))
0202         else:
0203             FchA=FchA1(math.sqrt(exp))+FchA2(math.sqrt(exp))
0204         funvalue=(1-Q2min/exp)*FchA**2
0205         return funvalue
0206 
0207     if A_x_array[ibeam] == None or (A_Q2max_save[ibeam] != Q2max):
0208         # we need to generate the grid first
0209         A_Q2max_save[ibeam] = Q2max
0210         xmaxQ2max=xmaxvalue(Q2max)
0211         log10xmaxQ2maxm1=math.log10(1/xmaxQ2max)
0212         A_x_array[ibeam]=[]
0213         A_xmax_array[ibeam]=[]
0214         A_fmax_array[ibeam]=[]
0215         for log10xm1 in range(10):
0216             for j in range(10):
0217                 tlog10xm1=log10xmaxQ2maxm1+0.1*j+log10xm1
0218                 A_x_array[ibeam].append(tlog10xm1)
0219                 xx=10**(-tlog10xm1)
0220                 if log10xm1 == 0 and j == 0:
0221                     max_Q2 = logQ2oQ02max
0222                     max_fun = distfun(xx,max_Q2)
0223                     A_xmax_array[ibeam].append(max_Q2)
0224                     A_fmax_array[ibeam].append(max_fun)
0225                 else:
0226                     max_Q2 = optimize.fmin(lambda x0: -distfun(xx,x0),\
0227                                                     (logQ2oQ02max+logQ2oQ02min)/2,\
0228                                                full_output=False,disp=False)
0229                     max_fun = distfun(xx,max_Q2[0])
0230                     A_xmax_array[ibeam].append(max_Q2[0])
0231                     A_fmax_array[ibeam].append(max_fun)
0232         A_x_array[ibeam]=numpy.array(A_x_array[ibeam])
0233         A_xmax_array[ibeam]=numpy.array(A_xmax_array[ibeam])
0234         A_fmax_array[ibeam]=numpy.array(A_fmax_array[ibeam])
0235         A_xmax_interp[ibeam]=interpolate.interp1d(A_x_array[ibeam],A_xmax_array[ibeam])
0236         A_fmax_interp[ibeam]=interpolate.interp1d(A_x_array[ibeam],A_fmax_array[ibeam])
0237     log10xm1=math.log10(1/x)
0238     max_x = A_xmax_interp[ibeam](log10xm1)
0239     max_fun = A_fmax_interp[ibeam](log10xm1)
0240     logQ2oQ02now=logQ2oQ02min
0241     while True:
0242         r1=random.random() # a random float number between 0 and 1
0243         logQ2oQ02now=(logQ2oQ02max-logQ2oQ02min)*r1+logQ2oQ02min
0244         w=distfun(x,logQ2oQ02now)/max_fun
0245         r2=random.random() # a random float number between 0 and 1
0246         if r2 <= w: break
0247     Q2v=math.exp(logQ2oQ02now)*Q02
0248     return Q2v
0249 
0250 #stream=open("Q2.dat",'w')
0251 #for i in range(100000):
0252 #    Q2v=generate_Q2_epa_ion(1,1e-1,1.0,WoodsSaxon['Pb208'][0],\
0253 #                                WoodsSaxon['Pb208'][1],WoodsSaxon['Pb208'][2])
0254 #    stream.write('%12.7e\n'%Q2v)
0255 #stream.close()
0256 
0257 def boostl(Q,PBOO,P):
0258     """Boost P via PBOO with PBOO^2=Q^2 to PLB"""
0259     # it boosts P from (Q,0,0,0) to PBOO
0260     # if P=(PBOO[0],-PBOO[1],-PBOO[2],-PBOO[3])
0261     # it will boost P to (Q,0,0,0)
0262     PLB=[0,0,0,0] # energy, px, py, pz in unit of GeV
0263     PLB[0]=(PBOO[0]*P[0]+PBOO[3]*P[3]+PBOO[2]*P[2]+PBOO[1]*P[1])/Q
0264     FACT=(PLB[0]+P[0])/(Q+PBOO[0])
0265     for j in range(1,4):
0266         PLB[j]=P[j]+FACT*PBOO[j]
0267     return PLB
0268 
0269 def boostl2(Q,PBOO1,PBOO2,P):
0270     """Boost P from PBOO1 (PBOO1^2=Q^2) to PBOO2 (PBOO2^2=Q^2) frame"""
0271     PBOO10=[PBOO1[0],-PBOO1[1],-PBOO1[2],-PBOO1[3]]
0272     PRES=boostl(Q,PBOO10,P) # PRES is in (Q,0,0,0) frame
0273     PLB=boostl(Q,PBOO2,PRES) # PLB is in PBOO2 frame
0274     return PLB
0275 
0276 def boostToEcm(E1,E2,pext):
0277     Ecm=2*math.sqrt(E1*E2)
0278     PBOO=[E1+E2,0,0,E2-E1]
0279     pext2=copy.deepcopy(pext)
0280     for j in range(len(pext)):
0281         pext2[j]=boostl(Ecm,PBOO,pext[j])
0282     return pext2
0283 
0284 def boostFromEcm(E1,E2,pext):
0285     Ecm=2*math.sqrt(E1*E2)
0286     PBOO=[E1+E2,0,0,E1-E2]
0287     pext2=copy.deepcopy(pext)
0288     for j in range(len(pext)):
0289         pext2[j]=boostl(Ecm,PBOO,pext[j])
0290     return pext2
0291 
0292 def InitialMomentumReshuffle(Ecm,x1,x2,Q1,Q2,pext):
0293     r1=random.random() # a random float number between 0 and 1
0294     r2=random.random() # a random float number between 0 and 1
0295     ph1=2*math.pi*r1
0296     ph2=2*math.pi*r2
0297     Kperp2=Q1**2+Q2**2+2*Q1*Q2*math.cos(ph1-ph2)
0298     Kperp2max=Ecm**2*(min(1,x1/x2,x2/x1)-x1*x2)
0299     if Kperp2 >= Kperp2max:
0300         return None
0301     x1bar=math.sqrt(x1/x2*Kperp2/Ecm**2+x1**2)
0302     x2bar=math.sqrt(x2/x1*Kperp2/Ecm**2+x2**2)
0303     if x1bar >= 1.0 or x2bar >= 1.0: return None
0304     pext2=copy.deepcopy(pext)
0305     # new initial state
0306     pext2[0][0]=Ecm/2*x1bar
0307     pext2[0][1]=Q1*math.cos(ph1)
0308     pext2[0][2]=Q1*math.sin(ph1)
0309     pext2[0][3]=Ecm/2*x1bar
0310     pext2[1][0]=Ecm/2*x2bar
0311     pext2[1][1]=Q2*math.cos(ph2)
0312     pext2[1][2]=Q2*math.sin(ph2)
0313     pext2[1][3]=-Ecm/2*x2bar
0314     # new final state
0315     PBOO1=[0,0,0,0]
0316     PBOO2=[0,0,0,0]
0317     for j in range(4):
0318         PBOO1[j]=pext[0][j]+pext[1][j]
0319         PBOO2[j]=pext2[0][j]+pext2[1][j]
0320     Q=math.sqrt(x1*x2)*Ecm
0321     for j in range(2,len(pext)):
0322         pext2[j]=boostl2(Q,PBOO1,PBOO2,pext[j])
0323     return pext2
0324 
0325 
0326 headers=[]
0327 inits=[]
0328 events=[]
0329 ninit0=0
0330 ninit1=0
0331 firstinit=""
0332 E_beam1=0
0333 E_beam2=0
0334 PID_beam1=0
0335 PID_beam2=0
0336 
0337 nevent=0
0338 
0339 ilil=0
0340 for i,file in enumerate(files):
0341     stream=open(file,'r')
0342     headQ=True
0343     initQ=False
0344     iinit=-1
0345     ievent=-1
0346     eventQ=False
0347     this_event=[]
0348     n_particles=0
0349     rwgtQ=False
0350     mgrwtQ=False
0351     procid=None
0352     proc_dict={}
0353     for line in stream:
0354         sline=line.replace('\n','')
0355         if "<init>" in line or "<init " in line:
0356             initQ=True
0357             headQ=False
0358             iinit=iinit+1
0359             if i==0: inits.append(sline)
0360         elif headQ and i == 0:
0361             headers.append(sline)
0362         elif "</init>" in line or "</init " in line:
0363             initQ=False
0364             iinit=-1
0365             if i==0: inits.append(sline)
0366         elif initQ:
0367             iinit=iinit+1
0368             if iinit == 1:
0369                 if i == 0:
0370                     firstinit=sline
0371                     ninit0=len(inits)
0372                     inits.append(sline)
0373                     firstinit=' '.join(firstinit.split()[:-1])
0374                     ff=firstinit.strip().split()
0375                     PID_beam1=int(ff[0])
0376                     PID_beam2=int(ff[1])
0377                     E_beam1=float(ff[2])
0378                     E_beam2=float(ff[3])
0379                     if abs(PID_beam1) != 2212 or abs(PID_beam2) != 2212:
0380                         raise ValueError( "Not a proton-proton collider")
0381                     ninit1=int(sline.split()[-1])
0382                 else:
0383                     ninit1=ninit1+int(sline.split()[-1])
0384                     sline=' '.join(sline.split()[:-1])
0385                     if not sline == firstinit:
0386                         raise Exception( "the beam information of the LHE files is not identical")
0387             elif iinit == 2:
0388                 procid=sline.split()[-1]
0389                 ilil=ilil+1
0390                 sline=' '.join(sline.split()[:-1])+(' %d'%(offset+ilil))
0391                 proc_dict[procid]=offset+ilil
0392                 if i == 0:
0393                     inits.append(sline)
0394                 else:
0395                     inits.insert(-1,sline)
0396             elif iinit >= 3:
0397                 if i == 0:
0398                     inits.append(sline)
0399                 else:
0400                     inits.insert(-1,sline)
0401             else:
0402                 raise Exception( "should not reach here. Do not understand the <init> block")
0403         elif "<event>" in line or "<event " in line:
0404             eventQ=True
0405             ievent=ievent+1
0406             events.append(sline)
0407         elif "</event>" in line or "</event " in line:
0408             nevent=nevent+1
0409             eventQ=False
0410             rwgtQ=False
0411             mgrwtQ=False
0412             ievent=-1
0413             this_event=[]
0414             n_particles=0
0415             events.append(sline)
0416             #if nevent >= 10: break
0417         elif eventQ:
0418             ievent=ievent+1
0419             if ievent == 1:
0420                 found=False
0421                 for procid,new_procid in proc_dict.items():
0422                     if ' '+procid+' ' not in sline: continue
0423                     procpos=sline.index(' '+procid+' ')
0424                     found=True
0425                     sline=sline[:procpos]+(' %d'%(new_procid))+sline[procpos+len(' '+procid):]
0426                     break
0427                 if not found: raise Exception( "do not find the correct proc id !")
0428                 n_particles=int(sline.split()[0])
0429                 #procpos=sline.index(' '+procid)
0430                 #sline=sline[:procpos]+(' %d'%(1+i))+sline[procpos+len(' '+procid):]
0431             elif "<rwgt" in sline:
0432                 rwgtQ=True
0433             elif "</rwgt" in sline:
0434                 rwgtQ=False
0435             elif "<mgrwt" in sline:
0436                 mgrwtQ=True
0437             elif "</mgrwt" in sline:
0438                 mgrwtQ=False                
0439             elif not rwgtQ and not mgrwtQ:
0440                 sline2=sline.split()
0441                 particle=[int(sline2[0]),int(sline2[1]),int(sline2[2]),int(sline2[3]),\
0442                               int(sline2[4]),int(sline2[5]),float(sline2[6]),float(sline2[7]),\
0443                               float(sline2[8]),float(sline2[9]),float(sline2[10]),\
0444                               float(sline2[11]),float(sline2[12])]
0445                 this_event.append(particle)
0446                 if ievent == n_particles+1:
0447                     # get the momenta and masses
0448                     x1=this_event[0][9]/E_beam1
0449                     x2=this_event[1][9]/E_beam2
0450                     pext=[]
0451                     mass=[]
0452                     for j in range(n_particles):
0453                         pext.append([this_event[j][9],this_event[j][6],\
0454                                          this_event[j][7],this_event[j][8]])
0455                         mass.append(this_event[j][10])
0456                     # first we need to boost from antisymmetric beams to symmetric beams
0457                     if E_beam1 != E_beam2:
0458                         pext=boostToEcm(E_beam1,E_beam2,pext)
0459                     Ecm=2*math.sqrt(E_beam1*E_beam2)
0460                     pext_new = None
0461                     Q1=0
0462                     Q2=0
0463                     while pext_new == None:
0464                         # generate Q1 and Q2
0465                         if nuclei[0] == 'p':
0466                             Q12=generate_Q2_epa_proton(x1,Q2max)
0467                         else:
0468                             RA,aA,wA=WoodsSaxon[nuclei[0]]
0469                             Q12=generate_Q2_epa_ion(0,x1,Q2max,RA,aA,wA)
0470                         if nuclei[1] == 'p':
0471                             Q22=generate_Q2_epa_proton(x2,Q2max)
0472                         else:
0473                             if nuclei[0] == nuclei[1]:
0474                                 RA,aA,wA=WoodsSaxon[nuclei[0]]
0475                                 Q22=generate_Q2_epa_ion(0,x2,Q2max,RA,aA,wA)
0476                             else:
0477                                 RA,aA,wA=WoodsSaxon[nuclei[1]]
0478                                 Q22=generate_Q2_epa_ion(1,x2,Q2max,RA,aA,wA)
0479                         Q1=math.sqrt(Q12)
0480                         Q2=math.sqrt(Q22)
0481                         # perform the initial momentum reshuffling
0482                         pext_new=InitialMomentumReshuffle(Ecm,x1,x2,Q1,Q2,pext)
0483 
0484                     if E_beam1 != E_beam2:
0485                         # boost back from the symmetric beams to antisymmetric beams
0486                         pext_new=boostFromEcm(E_beam1,E_beam2,pext_new)
0487                     # update the event information
0488                     # negative invariant mass means negative invariant mass square (-Q**2, spacelike)
0489                     this_event[0][10]=-Q1
0490                     this_event[1][10]=-Q2
0491                     for j in range(n_particles):
0492                         this_event[j][9]=pext_new[j][0]
0493                         this_event[j][6]=pext_new[j][1]
0494                         this_event[j][7]=pext_new[j][2]
0495                         this_event[j][8]=pext_new[j][3]
0496                         newsline="      %d    %d     %d    %d    %d    %d  %12.7e  %12.7e  %12.7e  %12.7e  %12.7e  %12.7e  %12.7e"%tuple(this_event[j])
0497                         events.append(newsline)
0498                 continue
0499             events.append(sline)
0500     stream.close()
0501 
0502 # modify the number of process information
0503 firstinit=firstinit+(' %d'%ninit1)
0504 inits[ninit0]=firstinit
0505 
0506 text='\n'.join(headers)+'\n'
0507 text=text+'\n'.join(inits)+'\n'
0508 text=text+'\n'.join(events)
0509 if '<LesHouchesEvents' in headers[0]: text=text+'\n</LesHouchesEvents>\n'
0510 
0511 stream=open(outfile,'w')
0512 stream.write(text)
0513 stream.close()
0514 print ("INFO: The final produced lhe file is %s"%outfile)
0515
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!