| Trees | Indices | Help |
|---|
|
|
1 ################################################################################
2 #
3 # Copyright (c) 2009 The MadGraph Development team and Contributors
4 #
5 # This file is a part of the MadGraph 5 project, an application which
6 # automatically generates Feynman diagrams and matrix elements for arbitrary
7 # high-energy processes in the Standard Model and beyond.
8 #
9 # It is subject to the MadGraph license which should accompany this
10 # distribution.
11 #
12 # For more information, please visit: http://madgraph.phys.ucl.ac.be
13 #
14 ################################################################################
15 """Methods and classes to group subprocesses according to initial
16 states, and produce the corresponding grouped subprocess directories."""
17
18 import array
19 import copy
20 import fractions
21 import glob
22 import itertools
23 import logging
24 import os
25 import re
26 import shutil
27 import subprocess
28
29 import madgraph.core.base_objects as base_objects
30 import madgraph.core.diagram_generation as diagram_generation
31 import madgraph.core.helas_objects as helas_objects
32 import madgraph.iolibs.drawing_eps as draw
33 import madgraph.iolibs.files as files
34 import madgraph.iolibs.file_writers as writers
35 import madgraph.iolibs.template_files as template_files
36 import madgraph.iolibs.ufo_expression_parsers as parsers
37
38 import madgraph.various.misc as misc
39
40 import aloha.create_aloha as create_aloha
41
42 import models.sm.write_param_card as write_param_card
43 from madgraph import MG5DIR
44 from madgraph.iolibs.files import cp, ln, mv
45 _file_path = os.path.split(os.path.dirname(os.path.realpath(__file__)))[0] + '/'
46 logger = logging.getLogger('madgraph.group_subprocs')
47
48 #===============================================================================
49 # DiagramTag class to identify matrix elements
50 #===============================================================================
51
52 -class IdentifyConfigTag(diagram_generation.DiagramTag):
53 """DiagramTag daughter class to identify diagrams giving the same
54 config. Need to compare leg number, mass, width, and color."""
55
56 @staticmethod
58 """Returns the end link for a leg needed to identify matrix
59 elements: ((leg numer, state, spin, self_antipart, mass,
60 width, color, decay and is_part), number)."""
61
62 part = model.get_particle(leg.get('id'))
63
64 return [((leg.get('number'),
65 part.get('mass'), part.get('width'), part.get('color')),
66 leg.get('number'))]
67
68 @staticmethod
70 """Returns the info needed to identify matrix elements:
71 interaction color, lorentz, coupling, and wavefunction
72 spin, self_antipart, mass, width, color, decay and
73 is_part. Note that is_part needs to be flipped if we move the
74 final vertex around."""
75
76 inter = model.get_interaction(vertex.get('id'))
77 ret_list = 0
78
79 if last_vertex:
80 return (ret_list,)
81 else:
82 part = model.get_particle(vertex.get('legs')[-1].get('id'))
83 return ((part.get('color'),
84 part.get('mass'), part.get('width')),
85 ret_list)
86
87 @staticmethod
89 """Move the wavefunction part of vertex id appropriately"""
90
91 if len(new_vertex) == 1 and len(old_vertex) == 2:
92 # We go from a last link to next-to-last link - add propagator info
93 return (old_vertex[0],new_vertex[0])
94 elif len(new_vertex) == 2 and len(old_vertex) == 1:
95 # We go from next-to-last link to last link - remove propagator info
96 return (new_vertex[1],)
97 # We should not get here
98 raise diagram_generation.DiagramTag.DiagramTagError, \
99 "Error in IdentifyConfigTag, wrong setup of vertices in link."
100
101 #===============================================================================
102 # SubProcessGroup
103 #===============================================================================
104
105 -class SubProcessGroup(base_objects.PhysicsObject):
106 """Class to group a number of amplitudes with same initial states
107 into a subprocess group"""
108
110 """Define object and give default values"""
111
112 self['number'] = 0
113 self['name'] = ""
114 self['amplitudes'] = diagram_generation.AmplitudeList()
115 self['matrix_elements'] = helas_objects.HelasMatrixElementList()
116 self['mapping_diagrams'] = []
117 self['diagram_maps'] = {}
118 self['diagrams_for_configs'] = []
119 self['amplitude_map'] = {}
120
122 """Filter for valid property values."""
123
124 if name == 'number':
125 if not isinstance(value, int):
126 raise self.PhysicsObjectError, \
127 "%s is not a valid int object" % str(value)
128 if name == 'name':
129 if not isinstance(value, str):
130 raise self.PhysicsObjectError, \
131 "%s is not a valid str object" % str(value)
132 if name == 'amplitudes':
133 if not isinstance(value, diagram_generation.AmplitudeList):
134 raise self.PhysicsObjectError, \
135 "%s is not a valid amplitudelist" % str(value)
136 if name in ['mapping_diagrams', 'diagrams_for_configs']:
137 if not isinstance(value, list):
138 raise self.PhysicsObjectError, \
139 "%s is not a valid list" % str(value)
140 if name == 'diagram_maps':
141 if not isinstance(value, dict):
142 raise self.PhysicsObjectError, \
143 "%s is not a valid dict" % str(value)
144 if name == 'matrix_elements':
145 if not isinstance(value, helas_objects.HelasMatrixElementList):
146 raise self.PhysicsObjectError, \
147 "%s is not a valid HelasMatrixElementList" % str(value)
148
149 if name == 'amplitude_map':
150 if not isinstance(value, dict):
151 raise self.PhysicsObjectError, \
152 "%s is not a valid dict object" % str(value)
153
154 return True
155
157 """Return diagram property names as a nicely sorted list."""
158
159 return ['number', 'name', 'amplitudes', 'mapping_diagrams',
160 'diagram_maps', 'matrix_elements', 'amplitude_map']
161
162 # Enhanced get function
164 """Get the value of the property name."""
165
166 if name == 'matrix_elements' and not self[name]:
167 self.generate_matrix_elements()
168
169 if name in ['mapping_diagrams', 'diagram_maps'] and not self[name]:
170 self.set_mapping_diagrams()
171
172 if name in ['diagrams_for_configs'] and not self[name]:
173 self.set_diagrams_for_configs()
174
175 return super(SubProcessGroup, self).get(name)
176
178 """Set mapping_diagrams and diagram_maps, to prepare for
179 generation of the super-config.inc files."""
180
181 # Find the mapping diagrams
182 mapping_diagrams, diagram_maps = \
183 self.find_mapping_diagrams()
184
185 self.set('mapping_diagrams', mapping_diagrams)
186 self.set('diagram_maps', diagram_maps)
187
188 #===========================================================================
189 # generate_matrix_elements
190 #===========================================================================
192 """Create a HelasMultiProcess corresponding to the amplitudes
193 in self"""
194
195 if not self.get('amplitudes'):
196 raise self.PhysicsObjectError, \
197 "Need amplitudes to generate matrix_elements"
198
199 amplitudes = copy.copy(self.get('amplitudes'))
200
201 self.set('matrix_elements',
202 helas_objects.HelasMultiProcess.\
203 generate_matrix_elements(amplitudes))
204
205 self.set('amplitudes', diagram_generation.AmplitudeList())
206
208 """Generate a convenient name for the group, based on and
209 masses"""
210
211 beam = [l.get('id') for l in process.get('legs') if not l.get('state')]
212 fs = [l.get('id') for l in process.get('legs') if l.get('state')]
213 name = ""
214 for beam in beam:
215 part = process.get('model').get_particle(beam)
216 if part.get('mass').lower() == 'zero' and part.is_fermion() and \
217 part.get('color') != 1:
218 name += "q"
219 elif part.get('mass').lower() == 'zero' and part.is_fermion() and \
220 part.get('color') == 1 and part.get('pdg_code') % 2 == 1:
221 name += "l"
222 elif part.get('mass').lower() == 'zero' and part.is_fermion() and \
223 part.get('color') == 1 and part.get('pdg_code') % 2 == 0:
224 name += "vl"
225 else:
226 name += part.get_name().replace('~', 'x').\
227 replace('+', 'p').replace('-', 'm')
228 name += "_"
229 for fs_part in fs:
230 part = process.get('model').get_particle(fs_part)
231 if part.get('mass').lower() == 'zero' and part.get('color') != 1 \
232 and part.get('spin') == 2:
233 name += "q" # "j"
234 elif part.get('mass').lower() == 'zero' and part.get('color') == 1 \
235 and part.get('spin') == 2:
236 if part.get('charge') == 0:
237 name += "vl"
238 else:
239 name += "l"
240 else:
241 name += part.get_name().replace('~', 'x').\
242 replace('+', 'p').replace('-', 'm')
243
244 for dc in process.get('decay_chains'):
245 name += "_" + self.generate_name(dc)
246
247 return name
248
250 """Get number of external and initial particles for this group"""
251
252 assert self.get('matrix_elements'), \
253 "Need matrix element to call get_nexternal_ninitial"
254
255 return self.get('matrix_elements')[0].\
256 get_nexternal_ninitial()
257
259 """Get number of configs for this group"""
260
261 model = self.get('matrix_elements')[0].get('processes')[0].\
262 get('model')
263
264 next, nini = self.get_nexternal_ninitial()
265
266 return sum([md.get_num_configs(model, nini) for md in
267 self.get('mapping_diagrams')])
268
270 """Find all unique diagrams for all processes in this
271 process class, and the mapping of their diagrams unto this
272 unique diagram."""
273
274 assert self.get('matrix_elements'), \
275 "Need matrix elements to run find_mapping_diagrams"
276
277 matrix_elements = self.get('matrix_elements')
278 model = matrix_elements[0].get('processes')[0].get('model')
279 # mapping_diagrams: The configurations for the non-reducable
280 # diagram topologies
281 mapping_diagrams = []
282 # equiv_diags: Tags identifying diagrams that correspond to
283 # the same configuration
284 equiv_diagrams = []
285 # diagram_maps: A dict from amplitude number to list of
286 # diagram maps, pointing to the mapping_diagrams (starting at
287 # 1). Diagrams with multi-particle vertices will have 0.
288 diagram_maps = {}
289 masswidth_to_pdg = {}
290
291 for ime, me in enumerate(matrix_elements):
292 diagrams = me.get('base_amplitude').get('diagrams')
293 # Check the minimal number of legs we need to include in order
294 # to make sure we'll have some valid configurations
295 max_legs = min([max([len(v.get('legs')) for v in \
296 d.get('vertices') if v.get('id') > 0]) \
297 for d in diagrams])
298 diagram_maps[ime] = []
299 for diagram in diagrams:
300 # Only use diagrams with all vertices == min_legs
301 if any([len(v.get('legs')) > max_legs \
302 for v in diagram.get('vertices') if v.get('id') > 0]):
303 diagram_maps[ime].append(0)
304 continue
305 # Create the equivalent diagram, in the format
306 # [[((ext_number1, mass_width_id1), ..., )],
307 # ...] (for each vertex)
308 equiv_diag = IdentifyConfigTag(diagram, model)
309 try:
310 diagram_maps[ime].append(equiv_diagrams.index(\
311 equiv_diag) + 1)
312 except ValueError:
313 equiv_diagrams.append(equiv_diag)
314 mapping_diagrams.append(diagram)
315 diagram_maps[ime].append(equiv_diagrams.index(\
316 equiv_diag) + 1)
317 return mapping_diagrams, diagram_maps
318
320 """Find the diagrams (number + 1) for all subprocesses
321 corresponding to config number iconfig. Return 0 for subprocesses
322 without corresponding diagram. Note that the iconfig should
323 start at 0."""
324
325 assert self.get('diagram_maps'), \
326 "Need diagram_maps to run get_subproc_diagrams_for_config"
327
328 subproc_diagrams = []
329 for iproc in \
330 range(len(self.get('matrix_elements'))):
331 try:
332 subproc_diagrams.append(self.get('diagram_maps')[iproc].\
333 index(iconfig + 1) + 1)
334 except ValueError:
335 subproc_diagrams.append(0)
336
337 return subproc_diagrams
338
340 """Get a list of all diagrams_for_configs"""
341
342 subproc_diagrams_for_config = []
343 for iconf in range(len(self.get('mapping_diagrams'))):
344 subproc_diagrams_for_config.append(\
345 self.get_subproc_diagrams_for_config(iconf))
346
347 self['diagrams_for_configs'] = subproc_diagrams_for_config
348
349 #===========================================================================
350 # group_amplitudes
351 #===========================================================================
352 @staticmethod
354 """Return a SubProcessGroupList with the amplitudes divided
355 into subprocess groups"""
356
357 assert isinstance(amplitudes, diagram_generation.AmplitudeList), \
358 "Argument to group_amplitudes must be AmplitudeList"
359
360 logger.info("Organizing processes into subprocess groups")
361
362 process_classes = SubProcessGroup.find_process_classes(amplitudes)
363 ret_list = SubProcessGroupList()
364 process_class_numbers = sorted(list(set(process_classes.values())))
365 for num in process_class_numbers:
366 amp_nums = [key for (key, val) in process_classes.items() if \
367 val == num]
368 group = SubProcessGroup()
369 group.set('amplitudes',
370 diagram_generation.AmplitudeList([amplitudes[i] for i in \
371 amp_nums]))
372 group.set('number', group.get('amplitudes')[0].get('process').\
373 get('id'))
374 group.set('name', group.generate_name(\
375 group.get('amplitudes')[0].get('process')))
376 ret_list.append(group)
377
378 return ret_list
379
380 @staticmethod
382 """Find all different process classes, classified according to
383 initial state and final state. For initial state, we
384 differentiate fermions, antifermions, gluons, and masses. For
385 final state, only masses."""
386
387 assert isinstance(amplitudes, diagram_generation.AmplitudeList), \
388 "Argument to find_process_classes must be AmplitudeList"
389 assert amplitudes
390
391 model = amplitudes[0].get('process').get('model')
392 proc_classes = []
393 amplitude_classes = {}
394
395 for iamp, amplitude in enumerate(amplitudes):
396 process = amplitude.get('process')
397 is_parts = [model.get_particle(l.get('id')) for l in \
398 process.get('legs') if not \
399 l.get('state')]
400 fs_parts = [model.get_particle(l.get('id')) for l in \
401 process.get('legs') if l.get('state')]
402 diagrams = amplitude.get('diagrams')
403
404 # This is where the requirements for which particles to
405 # combine are defined. Include p.get('is_part') in
406 # is_parts selection to distinguish between q and qbar,
407 # remove p.get('spin') from fs_parts selection to combine
408 # q and g into "j"
409 proc_class = [ [(p.is_fermion(), ) \
410 for p in is_parts], # p.get('is_part')
411 [(p.get('mass'), p.get('spin'),
412 abs(p.get('color')),l.get('onshell')) for (p, l) \
413 in zip(is_parts + fs_parts, process.get('legs'))],
414 amplitude.get('process').get('id'),
415 process.get('id')]
416 try:
417 amplitude_classes[iamp] = proc_classes.index(proc_class)
418 except ValueError:
419 proc_classes.append(proc_class)
420 amplitude_classes[iamp] = proc_classes.index(proc_class)
421
422 return amplitude_classes
423
424 #===============================================================================
425 # SubProcessGroupList
426 #===============================================================================
427 -class SubProcessGroupList(base_objects.PhysicsObjectList):
428 """List of SubProcessGroup objects"""
429
434
436 """Extract the list of matrix elements"""
437 return helas_objects.HelasMatrixElementList(\
438 sum([group.get('matrix_elements') for group in self], []))
439
441 """Return the list of ALOHA routines used in these matrix elements"""
442
443 return helas_objects.HelasMultiProcess(
444 {'matrix_elements': self.get_matrix_elements()}).get_used_lorentz()
445
447 """Return the list of ALOHA routines used in these matrix elements"""
448
449 return helas_objects.HelasMultiProcess(
450 {'matrix_elements': self.get_matrix_elements()}).get_used_couplings()
451
452 #===============================================================================
453 # DecayChainSubProcessGroup
454 #===============================================================================
455
456 -class DecayChainSubProcessGroup(SubProcessGroup):
457 """Class to keep track of subprocess groups from a list of decay chains"""
458
460 """Define object and give default values"""
461
462 self['core_groups'] = SubProcessGroupList()
463 self['decay_groups'] = DecayChainSubProcessGroupList()
464 # decay_chain_amplitudes is the original DecayChainAmplitudeList
465 self['decay_chain_amplitudes'] = diagram_generation.DecayChainAmplitudeList()
466
468 """Filter for valid property values."""
469
470 if name == 'core_groups':
471 if not isinstance(value, SubProcessGroupList):
472 raise self.PhysicsObjectError, \
473 "%s is not a valid core_groups" % str(value)
474 if name == 'decay_groups':
475 if not isinstance(value, DecayChainSubProcessGroupList):
476 raise self.PhysicsObjectError, \
477 "%s is not a valid decay_groups" % str(value)
478 if name == 'decay_chain_amplitudes':
479 if not isinstance(value, diagram_generation.DecayChainAmplitudeList):
480 raise self.PhysicsObjectError, \
481 "%s is not a valid DecayChainAmplitudeList" % str(value)
482 return True
483
485 """Return diagram property names as a nicely sorted list."""
486
487 return ['core_groups', 'decay_groups', 'decay_chain_amplitudes']
488
490 """Returns a nicely formatted string of the content."""
491
492 mystr = ""
493 for igroup, group in enumerate(self.get('core_groups')):
494 mystr += " " * indent + "Group %d:\n" % (igroup + 1)
495 for amplitude in group.get('amplitudes'):
496 mystr = mystr + amplitude.nice_string(indent + 2) + "\n"
497
498 if self.get('decay_groups'):
499 mystr += " " * indent + "Decay groups:\n"
500 for dec in self.get('decay_groups'):
501 mystr = mystr + dec.nice_string(indent + 2) + "\n"
502
503 return mystr[:-1]
504
505 #===========================================================================
506 # generate_helas_decay_chain_subproc_groups
507 #===========================================================================
509 """Combine core_groups and decay_groups to give
510 HelasDecayChainProcesses and new diagram_maps.
511 """
512
513 # Combine decays
514 matrix_elements = \
515 helas_objects.HelasMultiProcess.generate_matrix_elements(\
516 diagram_generation.AmplitudeList(\
517 self.get('decay_chain_amplitudes')))
518
519 # For each matrix element, check which group it should go into and
520 # calculate diagram_maps
521 me_assignments = {}
522 for me in matrix_elements:
523 group_assignment = self.assign_group_to_decay_process(\
524 me.get('processes')[0])
525 assert group_assignment
526 try:
527 me_assignments[group_assignment].append(me)
528 except KeyError:
529 me_assignments[group_assignment] = [me]
530
531 # Create subprocess groups corresponding to the different
532 # group_assignments
533
534 subproc_groups = SubProcessGroupList()
535 for key in sorted(me_assignments.keys()):
536 group = SubProcessGroup()
537 group.set('matrix_elements', helas_objects.HelasMatrixElementList(\
538 me_assignments[key]))
539 group.set('number', group.get('matrix_elements')[0].\
540 get('processes')[0].get('id'))
541 group.set('name', group.generate_name(\
542 group.get('matrix_elements')[0].\
543 get('processes')[0]))
544 subproc_groups.append(group)
545
546 return subproc_groups
547
549 """Recursively identify which group process belongs to."""
550
551 # Determine properties for the decay chains
552 # The entries of group_assignments are:
553 # [(decay_index, (decay_group_index, ...)),
554 # diagram_map (updated), len(mapping_diagrams)]
555
556 group_assignments = []
557
558 for decay in process.get('decay_chains'):
559 # Find decay group that has this decay in it
560 ids = [l.get('id') for l in decay.get('legs')]
561 decay_groups = [(i, group) for (i, group) in \
562 enumerate(self.get('decay_groups')) \
563 if any([ids in [[l.get('id') for l in \
564 a.get('process').get('legs')] \
565 for a in g.get('amplitudes')] \
566 for g in group.get('core_groups')])]
567
568 for decay_group in decay_groups:
569
570 group_assignment = \
571 decay_group[1].assign_group_to_decay_process(decay)
572
573 if group_assignment:
574 group_assignments.append((decay_group[0], group_assignment))
575
576 if process.get('decay_chains') and not group_assignments:
577 return None
578
579 # Now calculate the corresponding properties for process
580
581 # Find core process group
582 ids = [(l.get('id'),l.get('onshell')) for l in process.get('legs')]
583 core_groups = [(i, group) for (i, group) in \
584 enumerate(self.get('core_groups')) \
585 if ids in [[(l.get('id'),l.get('onshell')) for l in \
586 a.get('process').get('legs')] \
587 for a in group.get('amplitudes')] \
588 and process.get('id') == group.get('number')]
589
590 if not core_groups:
591 return None
592
593 assert len(core_groups) == 1
594
595 core_group = core_groups[0]
596 # This is the first return argument - the chain of group indices
597 group_assignment = (core_group[0],
598 tuple([g for g in group_assignments]))
599
600 if not group_assignments:
601 # No decays - return the values for this process
602 return group_assignment
603
604 return group_assignment
605
606 #===========================================================================
607 # group_amplitudes
608 #===========================================================================
609 @staticmethod
611 """Recursive function. Starting from a DecayChainAmplitude,
612 return a DecayChainSubProcessGroup with the core amplitudes
613 and decay chains divided into subprocess groups"""
614
615 assert isinstance(decay_chain_amps, diagram_generation.DecayChainAmplitudeList), \
616 "Argument to group_amplitudes must be DecayChainAmplitudeList"
617
618 # Collect all amplitudes
619 amplitudes = diagram_generation.AmplitudeList()
620 for amp in decay_chain_amps:
621 amplitudes.extend(amp.get('amplitudes'))
622
623 # Determine core process groups
624 core_groups = SubProcessGroup.group_amplitudes(amplitudes)
625
626 dc_subproc_group = DecayChainSubProcessGroup(\
627 {'core_groups': core_groups,
628 'decay_chain_amplitudes': decay_chain_amps})
629
630 decays = diagram_generation.DecayChainAmplitudeList()
631
632 # Recursively determine decay chain groups
633 for decay_chain_amp in decay_chain_amps:
634 decays.extend(decay_chain_amp.get('decay_chains'))
635
636 if decays:
637 dc_subproc_group.get('decay_groups').append(\
638 DecayChainSubProcessGroup.group_amplitudes(decays))
639
640 return dc_subproc_group
641
642
643
644
645 #===============================================================================
646 # DecayChainSubProcessGroupList
647 #===============================================================================
648 -class DecayChainSubProcessGroupList(base_objects.PhysicsObjectList):
649 """List of DecayChainSubProcessGroup objects"""
650
652 """Test if object obj is a valid element."""
653
654 return isinstance(obj, DecayChainSubProcessGroup)
655
| Trees | Indices | Help |
|---|
| Generated by Epydoc 3.0.1 on Tue Jul 24 08:02:14 2012 | http://epydoc.sourceforge.net |