Source Code for Module madgraph.various.rambo

  1  from __future__ import division 
  2  import math 
  3  import random 
  4   
  5 -class FortranList(list): 
  6       
  7 -    def __init__(self, min, max=None): 
  8          if max is None: 
  9              self.min = 1 
 10              self.max = min + 1 
 11          else: 
 12              self.min = min 
 13              self.max = max + 1 
 14          list.__init__(self,[0]*(self.max-self.min)) 
 15   
 16           
 17 -    def __getitem__(self, index): 
 18          assert self.min <= index < self.max, 'outside range %s <= %s < %s' % (self.min, index, self.max) 
 19          return list.__getitem__(self, index - self.min) 
 20       
 21 -    def __setitem__(self, index, value): 
 22          assert self.min <= index < self.max 
 23          return list.__setitem__(self, index - self.min , value) 
 24   
 25 -class DoubleFortranList(FortranList): 
 26   
 27 -    def __init__(self, min, max=(None,None)): 
 28           
 29          min1 = min[0] 
 30          max1 = max[0]  
 31          FortranList.__init__(self, min1, max1) 
 32           
 33          min2 = min[1] 
 34          max2 = max[1] 
 35           
 36          for i in range(len(self)): 
 37              list.__setitem__(self,i,FortranList(min2,max2)) 
 38           
 39 -    def __getitem__(self, index): 
 40          var1 = index[0] 
 41          var2 = index[1] 
 42          list1 = FortranList.__getitem__(self, var1) 
 43          return list1.__getitem__(var2) 
 44           
 45 -    def __setitem__(self, index, value): 
 46          var1 = index[0] 
 47          var2 = index[1] 
 48           
 49          list1 = FortranList.__getitem__(self, var1) 
 50          list1.__setitem__(var2, value) 
 51                   
 52   
 53 -class RAMBOError(Exception): 
 54      """ A Error class for RAMBO routine """ 
 55      pass 
 56   
 57 -def RAMBO(N,ET,XM): 
 58      """*********************************************************************** 
 59      *                       RAMBO                                         * 
 60      *    RA(NDOM)  M(OMENTA)  B(EAUTIFULLY)  O(RGANIZED)                  * 
 61      *                                                                     * 
 62      *    A DEMOCRATIC MULTI-PARTICLE PHASE SPACE GENERATOR                * 
 63      *    AUTHORS:  S.D. ELLIS,  R. KLEISS,  W.J. STIRLING                 * 
 64      *    -- ADJUSTED BY HANS KUIJF, WEIGHTS ARE LOGARITHMIC (20-08-90)    * 
 65      *    THIS IS PY VERSION 1.0 -  WRITTEN BY O. MATTELAER                * 
 66      *                                                                     * 
 67      *    N  = NUMBER OF PARTICLES                                         * 
 68      *    ET = TOTAL CENTRE-OF-MASS ENERGY                                 * 
 69      *    XM = PARTICLE MASSES ( DIM=NEXTERNAL-nincoming )                 * 
 70      *  RETURN                                                             * 
 71      *    P  = PARTICLE MOMENTA ( DIM=(4,NEXTERNAL-nincoming) )            * 
 72      *    WT = WEIGHT OF THE EVENT                                         * 
 73      ***********************************************************************""" 
 74      # RUN PARAMETER 
 75      acc = 1e-14 
 76      itmax = 6 
 77      ibegin = 0 
 78      iwarn = FortranList(5) 
 79      Nincoming = 2 
 80       
 81      # Object Initialization 
 82      Z = FortranList(N) 
 83      Q = DoubleFortranList((4,N)) 
 84      P = DoubleFortranList((4,N)) 
 85      R = FortranList(4) 
 86      B = FortranList(3) 
 87      XM2 = FortranList(N) 
 88      P2 = FortranList(N) 
 89      E = FortranList(N) 
 90      V= FortranList(N) 
 91   
 92  # Check input object 
 93      assert isinstance(XM, FortranList) 
 94      assert XM.min == 1 
 95      assert XM.max == N+1  
 96   
 97  # INITIALIZATION STEP: FACTORIALS FOR THE PHASE SPACE WEIGHT                                                                                
 98      if not ibegin: 
 99        ibegin = 1 
100        twopi = 8 * math.atan(1) 
101        po2log = math.log(twopi/4) 
102        Z[2] = po2log 
103        for k in range(3, N+1): 
104            Z[k] = Z[k-1] + po2log - 2.*math.log(k-2) - math.log(k-1) 
105             
106  # CHECK ON THE NUMBER OF PARTICLES 
107        assert 1 < N < 101 
108   
109  # CHECK WHETHER TOTAL ENERGY IS SUFFICIENT; COUNT NONZERO MASSES 
110      xmt =  0 
111      nm = 0 
112      for i in range(1,N+1): 
113          if XM[i] != 0: 
114              nm +=1 
115          xmt += abs(XM[i]) 
116           
117      if xmt > ET: 
118          raise RAMBOError, ' Not enough energy in this case' 
119   
120  #                                                                                                                                           
121  # THE PARAMETER VALUES ARE NOW ACCEPTED                                                                                                     
122  #                                                                                                                                           
123  # GENERATE N MASSLESS MOMENTA IN INFINITE PHASE SPACE                                                                                       
124      for i in range(1,N+1): 
125          r1=random_nb(1) 
126          c = 2 * r1 -1 
127          s = math.sqrt(1 - c**2) 
128          f = twopi * random_nb(2) 
129          r1 = random_nb(3) 
130          r2 = random_nb(4) 
131           
132          Q[(4,i)]=-math.log(r1*r2) 
133          Q[(3,i)]= Q[(4,i)]*c 
134          Q[(2,i)]=Q[(4,i)]*s*math.cos(f) 
135          Q[(1,i)]=Q[(4,i)]*s*math.sin(f) 
136           
137  # CALCULATE THE PARAMETERS OF THE CONFORMAL TRANSFORMATION                                                                                  
138      for i in range(1, N+1): 
139          for k in range(1,5): 
140              R[k] = R[k] + Q[(k,i)] 
141      rmas = math.sqrt(R[4]**2-R[3]**2-R[2]**2-R[1]**2) 
142      for k in range(1,4): 
143          B[k] = - R[k]/rmas 
144   
145      g = R[4] / rmas 
146      a = 1.0 / (1+g) 
147      x = ET / rmas 
148   
149  # TRANSFORM THE Q'S CONFORMALLY INTO THE P'S 
150      for i in range(1, N+1): 
151          bq = B[1]*Q[(1,i)]+B[2]*Q[(2,i)]+B[3]*Q[(3,i)] 
152          for k in range(1,4): 
153              P[k,i] = x*(Q[(k,i)]+B[k]*(Q[(4,i)]+a*bq))  
154          P[(4,i)] = x*(g*Q[(4,i)]+bq) 
155           
156  # CALCULATE WEIGHT AND POSSIBLE WARNINGS                                                                                                    
157      wt = po2log 
158      if N != 2: 
159          wt = (2 * N-4) * math.log(ET) + Z[N] 
160      if wt < -180 and IWARN[1] < 5: 
161          print "RAMBO WARNS: WEIGHT = EXP(%f20.9) MAY UNDERFLOW" % wt 
162          IWARN[1] += 1 
163      if wt > 174 and IWARN[2] < 5:       
164          print " RAMBO WARNS: WEIGHT = EXP(%f20.9) MAY  OVERFLOW" % wt 
165          IWARN[2] += 1 
166   
167                                                                                                                                             
168  # RETURN FOR WEIGHTED MASSLESS MOMENTA                                                                                                      
169      if nm == 0:  
170          return P, wt 
171   
172   
173  # MASSIVE PARTICLES: RESCALE THE MOMENTA BY A FACTOR X 
174      xmax = math.sqrt(1-(xmt/ET)**2) 
175      for i in range(1,N+1): 
176          XM2[i] = XM[i] **2 
177          P2[i] = P[(4,i)]**2 
178      iter = 0 
179      x= xmax 
180      accu = ET * acc 
181     
182      while 1: 
183          f0 = -ET 
184          g0 = 0 
185          x2 = x**2 
186          for i in range(1, N+1): 
187              E[i] = math.sqrt(XM2[i]+x2*P2[i]) 
188              f0 += E[i] 
189              g0 += P2[i]/E[i]  
190          if abs(f0) <= accu: 
191              break 
192          iter += 1 
193          if iter > itmax: 
194              print "RAMBO WARNS: %s ITERATIONS DID NOT GIVE THE DESIRED ACCURACY = %s" \ 
195                      %(iter, f0) 
196              break 
197          x=x-f0/(x*g0) 
198      for i in range(1, N+1): 
199          V[i] = x * P[(4,i)] 
200          for k in range(1,4): 
201              P[(k,i)] = x * P[(k,i)] 
202          P[(4,i)] = E[i] 
203           
204  # CALCULATE THE MASS-EFFECT WEIGHT FACTOR                                                                                                   
205      wt2 = 1. 
206      wt3 = 0. 
207      for i in range(1, N+1): 
208          wt2 *= V[i]/E[i] 
209          wt3 += V[i]**2/E[i] 
210      wtm = (2.*N-3.)*math.log(x)+math.log(wt2/wt3*ET) 
211   
212  # RETURN FOR  WEIGHTED MASSIVE MOMENTA                                                                                                      
213      wt += wtm 
214      if(wt < -180 and iwarn[3] < 5): 
215          print " RAMBO WARNS: WEIGHT = EXP(%s) MAY UNDERFLOW" % wt 
216          iwarn[3] += 1 
217      if(wt > 174  and iwarn[4] > 5): 
218          print " RAMBO WARNS: WEIGHT = EXP(%s) MAY OVERFLOW" % wt 
219          iwarn[4] += 1 
220   
221  # RETURN 
222      return P, wt 
223   
224 -def random_nb(value): 
225      """ random number """ 
226      output = 0 
227      while output < 1e-16: 
228          output= random.uniform(0,1) 
229      return output 
230