~maddevelopers/mg5amcnlo/aMCatNLO : revision 552

1

#! /usr/bin/env python

2

3

from __future__ import division

4

5

################################################################################

6

#

7

8

#

9

# This file is a part of the MadGraph 5 project, an application which

10

# automatically generates Feynman diagrams and matrix elements for arbitrary

11

# high-energy processes in the Standard Model and beyond.

12

#

13

# It is subject to the MadGraph license which should accompany this

14

# distribution.

15

#

16

# For more information, please visit: http://madgraph.phys.ucl.ac.be

17

#

18

################################################################################

19

"""

20

####################################################################

21

#

22

# Routine to decay prodution events in a generic way,

23

# including spin correlation effects

24

#

25

# Ref: S. Frixione, E. Laenen, P. Motylinski, B. R. Webber

26

# JHEP 04 (2007) 081

27

#

28

#

29

#####################################################################

30

"""

31

32

import re

33

import os

34

import shutil

35

import logging

36

import time

37

import cmath

38

pjoin = os.path.join

39

from subprocess import Popen, PIPE, STDOUT

40

41

os.sys.path.append("../.")

42

import string

43

#import madgraph.core.base_objects as base_objects

44

#import madgraph.core.helas_objects as helas_objects

45

#import madgraph.core.diagram_generation as diagram_generation

46

47

import models.import_ufo as import_ufo

48

#import madgraph.various.process_checks as process_checks

49

#from time import strftime

50

#from madgraph.interface.madgraph_interface import MadGraphCmd

51

import madgraph.interface.master_interface as Cmd

52

import madgraph.interface.madevent_interface as me_interface

53

import madgraph.iolibs.files as files

54

import aloha

55

logger = logging.getLogger('decay.stdout') # -> stdout

56

logger_stderr = logging.getLogger('decay.stderr') # ->stderr

57

58

import random

59

import math

60

from madgraph import MG5DIR

61

import madgraph.various.misc as misc

62

#import time

63

64

class Event:

65

""" class to read an event, record the information, write down the event in the lhe format.

66

This class is used both for production and decayed events"""

67

68

def __init__(self, inputfile=None):

69

"""Store the name of the event file """

70

self.inputfile=inputfile

71

self.particle={}

72

73

def give_procdef(self, pid2label):

74

""" Return a string with the process in the format i j > k l ... """

75

proc_line=""

76

initial=0

77

for part in range(1,len(self.particle)+1):

78

if self.particle[part]["istup"]==-1:

79

initial+=1

80

proc_line+=pid2label[self.particle[part]["pid"]]+" "

81

if initial==2 :

82

proc_line+="> "

83

84

if self.particle[part]["istup"]==1:

85

proc_line+=pid2label[self.particle[part]["pid"]]+" "

86

return proc_line

87

88

89

def get_map_resonances(self, branchindex2label,pid2label):

90

"""

91

decay_processes is a dictionary 1 : "A1"

92

2 : "A2"

93

...

94

returns a map {mg index : index of the branch}

95

for each mg index associated with a resonance to be decayed

96

"""

97

# first build a dictionary {}

98

99

dico_map={}

100

dico_label={}

101

102

# dico_symmetry_fac={} # dico {label : symmetry factor}

103

branch_not_yet_found=branchindex2label.keys()

104

branch_found=[]

105

106

107

for index in self.event2mg.keys():

108

if self.event2mg[index]>0:

109

part=self.event2mg[index]

110

pid=self.particle[part]['pid']

111

found=0

112

for index, id in enumerate(branch_not_yet_found):

113

if pid2label[pid]==branchindex2label[id]:

114

found=1

115

dico_map[part]=id

116

label=pid2label[pid]

117

dico_label[part]=label

118

# if not dico_symmetry_fac.has_key(label):

119

# dico_symmetry_fac[label]=1

120

# else:

121

# dico_symmetry_fac[label]+=1

122

123

branch_found.append(id)

124

to_delete=index

125

break

126

if found: del branch_not_yet_found[to_delete]

127

128

if not found:

129

for index in branch_found:

130

if pid2label[pid]==branchindex2label[index]:

131

dico_map[part]=index

132

dico_label[part]=pid2label[pid]

133

134

# now get the symmetry factor:

135

# symm_fac=1.0

136

# print dico_symmetry_fac

137

# for label in dico_symmetry_fac.keys():

138

# for index in range(2,dico_symmetry_fac[label]+1):

139

# symm_fac=symm_fac*(index)

140

141

return dico_map, dico_label #, symm_fac

142

143

144

def give_momenta(self):

145

""" return the set of external momenta of the event,

146

in two different formats:

147

p is list momenta, with each momentum a list [E,px,py,pz]

148

string is a sting

149

"""

150

p=[]

151

string=""

152

for part in range(1,len(self.particle)+1):

153

if self.particle[part]["istup"]<2:

154

mom=[self.particle[part]["momentum"].E, \

155

self.particle[part]["momentum"].px,\

156

self.particle[part]["momentum"].py,\

157

self.particle[part]["momentum"].pz]

158

p.append(mom)

159

string+=str(self.particle[part]["momentum"].E)+" "+\

160

str(self.particle[part]["momentum"].px)\

161

+" "+str(self.particle[part]["momentum"].py)+\

162

" "+str(self.particle[part]["momentum"].pz)+"\n"

163

return p, string

164

165

def string_event_compact(self):

166

""" return a string with the momenta of the event written

167

in an easy readable way

168

"""

169

line=""

170

for part in range(1,len(self.particle)+1):

171

line+=str(self.particle[part]["pid"])+" "

172

line+=str(self.particle[part]["momentum"].px)+" "

173

line+=str(self.particle[part]["momentum"].py)+" "

174

line+=str(self.particle[part]["momentum"].pz)+" "

175

line+=str(self.particle[part]["momentum"].E)+" "

176

line+=str(self.particle[part]["momentum"].m)+" "

177

line+="\n"

178

return line

179

180

def string_event(self):

181

""" return a string with the information of the event written

182

in the lhe format.

183

"""

184

line="<event> \n"

185

line1=' %2d %6d %+13.7e %14.8e %14.8e %14.8e' % \

186

(self.nexternal,self.ievent,self.wgt,self.scale,self.aqed,self.aqcd)

187

line+=line1+"\n"

188

for item in range(1,len(self.event2mg.keys())+1):

189

part=self.event2mg[item]

190

if part>0:

191

particle_line=self.get_particle_line(self.particle[part])

192

else:

193

particle_line=self.get_particle_line(self.resonance[part])

194

line+=particle_line

195

196

if self.diese:

197

line += self.diese

198

if self.rwgt:

199

line += self.rwgt

200

line+="</event> \n"

201

return line

202

203

204

def get_particle_line(self,leg):

205

206

line=" %8d %2d %4d %4d %4d %4d %+13.7e %+13.7e %+13.7e %14.8e %14.8e %10.4e %10.4e" \

207

% (leg["pid"], leg["istup"],leg["mothup1"],leg["mothup2"],\

208

leg["colup1"],leg["colup2"],leg["momentum"].px,leg["momentum"].py,\

209

leg["momentum"].pz,leg["momentum"].E, leg["momentum"].m,\

210

0.0,float(leg["helicity"]) )

211

line+="\n"

212

return line

213

214

def reshuffle_resonances(self,mother):

215

""" reset the momentum of each resonance in the production event

216

to the sum of momenta of the daughters

217

"""

218

219

daughters=[]

220

for part in self.event2mg.keys():

221

index=self.event2mg[part]

222

if index>0:

223

if self.particle[index]["mothup1"]==mother:

224

daughters.append(index)

225

if index<0:

226

if self.resonance[index]["mothup1"]==mother:

227

daughters.append(index)

228

229

if len(daughters)!=2:

230

logger.info("Got more than 2 daughters for one particles")

231

logger.info("in one production event (before decay)")

232

233

234

if daughters[0]>0:

235

momentum_mother=self.particle[daughters[0]]["momentum"].copy()

236

else:

237

momentum_mother=self.resonance[daughters[0]]["momentum"].copy()

238

239

# there might be more than 2 daughters, add all their momentum to get the momentum of the mother

240

for index in range(1,len(daughters)):

241

if daughters[index]>0:

242

momentum_mother=momentum_mother.add(self.particle[daughters[index]]["momentum"])

243

else:

244

momentum_mother=momentum_mother.add(self.resonance[daughters[index]]["momentum"])

245

246

res=self.event2mg[mother]

247

del self.resonance[res]["momentum"]

248

self.resonance[res]["momentum"]=momentum_mother.copy()

249

250

# recurrence:

251

if self.resonance[res]["mothup1"]>2:

252

self.reshuffle_resonances(self.resonance[res]["mothup1"])

253

254

255

def reset_resonances(self):

256

""" re-evaluate the momentum of each resonance, based on the momenta

257

of the external particles

258

"""

259

260

mothers=[]

261

for part in self.particle.keys():

262

if self.particle[part]["mothup1"]>2 and \

263

self.particle[part]["mothup1"] not in mothers :

264

mothers.append(self.particle[part]["mothup1"])

265

self.reshuffle_resonances(self.particle[part]["mothup1"])

266

267

def assign_scale_line(self, line):

268

"""read the line corresponding to global event line

269

format of the line is:

270

Nexternal IEVENT WEIGHT SCALE AEW AS

271

"""

272

line = line.replace('d','e').replace('D','e')

273

inputs = line.split()

274

assert len(inputs) == 6

275

self.nexternal=int(inputs[0])

276

self.ievent=int(inputs[1])

277

self.wgt=float(inputs[2])

278

self.scale=float(inputs[3])

279

self.aqed=float(inputs[4])

280

self.aqcd=float(inputs[5])

281

282

283

284

def get_next_event(self):

285

""" read next event in the lhe event file """

286

line_type = 'none' # support type: init / event / rwgt

287

self.diese = ''

288

for line in self.inputfile:

289

line = line.lower()

290

if line=="":

291

continue

292

# Find special tag in the line

293

if line[0]=="#":

294

self.diese+=line

295

continue

296

if '<event>' in line:

297

#start new_event

298

line_type = 'init'

299

continue

300

elif '<rwgt>' in line:

301

#re-weighting information block

302

line_type = 'rwgt'

303

#No Continue! need to keep track of this line

304

elif '</event>' in line:

305

if not self.particle:

306

continue

307

self.shat=self.particle[1]["momentum"].dot(self.particle[2]["momentum"])

308

return 1

309

310

311

if line_type == 'none':

312

continue

313

# read the line and assign the date accordingly

314

elif line_type == 'init':

315

line_type = 'event'

316

self.assign_scale_line(line)

317

# initialize some local variable

318

index_prod=0

319

index_external=0

320

index_resonance=0

321

self.particle={}

322

self.resonance={}

323

self.max_col=500

324

self.diese=""

325

self.rwgt=""

326

self.event2mg={} # dict. event-like label <-> "madgraph-like" label

327

# in "madgraph-like", resonances are labeled -1, -2, ...

328

elif line_type == 'rwgt': #special aMC@NLO information

329

self.rwgt += line

330

if '</rwgt>' in line:

331

line_type = 'event'

332

elif line_type == 'event':

333

index_prod+=1

334

line=line.replace("\n","")

335

line = line.replace('d','e').replace('D','e')

336

inputs=line.split()

337

pid=int(inputs[0])

338

istup=int(inputs[1])

339

mothup1=int(inputs[2])

340

mothup2=int(inputs[3])

341

colup1=int(inputs[4])

342

if colup1>self.max_col:

343

self.max_col=colup1

344

colup2=int(inputs[5])

345

if colup2>self.max_col:

346

self.max_col=colup2

347

mom=momentum(float(inputs[9]),float(inputs[6]),float(inputs[7]),float(inputs[8]))

348

mass=float(inputs[10])

349

helicity=float(inputs[12])

350

if abs(istup)==1:

351

index_external+=1

352

self.event2mg[index_prod]=index_external

353

self.particle[index_external]={"pid":pid,"istup":istup,"mothup1":mothup1,\

354

"mothup2":mothup2,"colup1":colup1,"colup2":colup2,"momentum":mom,"mass":mass,"helicity":helicity}

355

elif istup==2:

356

index_resonance=index_resonance-1

357

self.event2mg[index_prod]=index_resonance

358

self.resonance[index_resonance]={"pid":pid,"istup":istup,"mothup1":mothup1,\

359

"mothup2":mothup2,"colup1":colup1,"colup2":colup2,"momentum":mom,"mass":mass,"helicity":helicity}

360

else:

361

logger.warning('unknown status in lhe file')

362

return "no_event"

363

364

class pid2label(dict):

365

""" dico pid:label for a given model"""

366

367

def __init__(self,model):

368

369

for particle in model["particles"]:

370

self[particle["pdg_code"]]=particle["name"]

371

self[-particle["pdg_code"]]=particle["antiname"]

372

373

class pid2color(dict):

374

""" dico pid:color rep. for a given model (ex: 5:-3 )"""

375

376

def __init__(self,model):

377

378

for particle in model["particles"]:

379

self[particle["pdg_code"]]=particle["color"]

380

self[-particle["pdg_code"]]=-particle["color"]

381

382

class label2pid(dict):

383

""" dico label:pid for a given model"""

384

385

def __init__(self,model):

386

387

for particle in model["particles"]:

388

self[particle["name"]]=particle.get_pdg_code()

389

self[particle["antiname"]]=-particle.get_pdg_code()

390

if particle['self_antipart']:

391

self[particle["name"]]=abs(self[particle["name"]])

392

self[particle["antiname"]]=abs(self[particle["antiname"]])

393

394

395

class mass_n_width(dict):

396

""" dictionary to extract easily the mass ad the width of the particle in the model

397

{name : {"mass": mass_value, "width": widht_value} }

398

I assume that masses and widths are real

399

"""

400

401

def __init__(self, model, banner):

402

""" Fill the dictionary based on the information in the banner """

403

404

self.model = model

405

self.banner = banner

406

for particle in model["particles"]:

407

self[particle["name"]]={"mass":0.0, "width":0.0}

408

if particle["mass"]!="ZERO":

409

self[particle["name"]]["mass"]=self.get_mass(particle["pdg_code"])

410

if particle["width"]!="ZERO":

411

self[particle["name"]]["width"]=self.get_width(particle["pdg_code"])

412

413

def get_mass(self,pid):

414

""" extract the mass of particle with PDG code=pid from the banner"""

415

try:

416

mass=self.banner.get('param_card', 'mass', pid)

417

except Exception, error:

418

mass=0.0

419

return mass

420

421

422

def get_width(self,pid):

423

""" extract the width of particle with PDG code=pid from the banner"""

424

try:

425

width=self.banner.get('param_card', 'decay', pid)

426

except Exception, error:

427

width=0.0

428

return width

429

430

431

class dc_branch(dict):

432

""" A dictionary to record information necessary to decay particles

433

{ -1 : {"d1": { "label": XX , "nb": YY }, "d2": { "label": XX , "nb": YY } },

434

-2 : {"d1": { "label": XX , "nb": YY }, "d2": { "label": XX , "nb": YY } },

435

....

436

}

437

"""

438

439

def __init__(self, process_line,model,banner,check):

440

441

self.banner=banner

442

self.label2pid=label2pid(model)

443

list_decays=process_line.split(",")

444

self.nb_decays=len(list_decays)

445

self["m_label2index"]={}

446

self["m_index2label"]={}

447

448

for dc_nb, dc in enumerate(list_decays):

449

if dc.find(">")<0: logger.warning('warning: invalid decay chain syntax')

450

mother=dc[:dc.find(">")].replace(" ","")

451

self["m_label2index"][mother]=-dc_nb-1

452

self["m_index2label"][-dc_nb-1]=mother

453

454

self["nexternal"]=self.find_tree(list_decays)

455

if (check): self.check_parameters()

456

# here I check that the relevant masses and widths can

457

# be extracted from the banner

458

459

def check_parameters(self):

460

for res in range(-1,-self.nb_decays-1,-1):

461

d1=self["tree"][res]["d1"]["index"]

462

d2=self["tree"][res]["d2"]["index"]

463

if (d1>0):

464

try:

465

logger.info('Mass of particle with id '\

466

+str(self.label2pid[self["tree"][res]["d1"]["label"]]))

467

logger.info(self.banner.param["Block mass"]\

468

[abs(self.label2pid[self["tree"][res]["d1"]["label"]])])

469

except Exception, error:

470

logger.info('The mass of particle with id '\

471

+str(self.label2pid[self["tree"][res]["d1"]["label"]]))

472

logger.info('was not defined in the param_card.dat')

473

mass=raw_input("Please enter the mass: ")

474

self.banner.param["Block mass"][abs(self.label2pid[self["tree"][res]["d1"]["label"]])]=mass

475

476

if (d2>0):

477

try:

478

logger.info('Mass of particle with id '\

479

+str(self.label2pid[self["tree"][res]["d2"]["label"]]))

480

logger.info( self.banner.param["Block mass"]\

481

[abs(self.label2pid[self["tree"][res]["d2"]["label"]])])

482

except Exception, error:

483

logger.info('The mass of particle with id '\

484

+str(self.label2pid[self["tree"][res]["d2"]["label"]]))

485

logger.info('was not defined in the param_card.dat')

486

mass=raw_input("Please enter the mass: ")

487

self.banner.param["Block mass"][abs(self.label2pid[self["tree"][res]["d2"]["label"]])]=mass

488

489

490

def transpole(self,pole,width):

491

492

""" routine for the generation of a p^2 according to

493

a Breit Wigner distribution

494

the generation window is

495

[ M_pole^2 - 30*M_pole*Gamma , M_pole^2 + 30*M_pole*Gamma ]

496

"""

497

498

zmin = math.atan(-30.0)/width

499

zmax = math.atan(30.0)/width

500

501

z=zmin+(zmax-zmin)*random.random()

502

y = pole+width*math.tan(width*z)

503

504

jac=(width/math.cos(width*z))**2*(zmax-zmin)

505

return y, jac

506

507

def generate_momenta(self,mom_init,ran, pid2width,pid2mass,resonnances,BW_effects):

508

"""Generate the momenta in each decay branch

509

If ran=1: the generation is random, with

510

a. p^2 of each resonance generated according to a BW distribution

511

b. cos(theta) and phi (angles in the rest frame of the decaying particle)

512

are generated according to a flat distribution (no grid)

513

the phase-space weight is also return (up to an overall normalization)

514

since it is needed in the unweighting procedure

515

516

If ran=0: evaluate the momenta based on the previously-generated p^2, cos(theta)

517

and phi in each splitting.

518

This is used in the reshuffling phase (e.g. when we give a mass to gluons

519

in the decay chain )

520

"""

521

index2mom={}

522

# pid2mom={} # a dict { pid : {"status":status, "momentum":momentum} }

523

524

index2mom[-1]={}

525

index2mom[-1]["momentum"]=mom_init

526

if index2mom[-1]['momentum'].m < 1e-3:

527

logger.debug('Decaying particle with m> 1e-3 GeV in generate_momenta')

528

index2mom[-1]["pid"]=self.label2pid[self["m_index2label"][-1]]

529

index2mom[-1]["status"]=2

530

weight=1.0

531

for res in range(-1,-self.nb_decays-1,-1):

532

# Here mA^2 has to be set to p^2:

533

#

534

# IF res=-1:

535

# p^2 has been either fixed to the value in the

536

# production lhe event, or generated according to a Breit-Wigner distr.

537

# during the reshuffling phase of the production event

538

# -> we just need to read the value here

539

# IF res<-1:

540

# p^2 has been generated during the previous iteration of this loop

541

# -> we just need to read the value here

542

543

mA=index2mom[res]["momentum"].m

544

if mA < 1.0e-3:

545

logger.debug('Warning: decaying parting with m<1 MeV in generate_momenta ')

546

#print index2mom[res]["momentum"].px

547

#print index2mom[res]["momentum"].py

548

#print index2mom[res]["momentum"].pz

549

#print index2mom[res]["momentum"].E

550

#print mA

551

#print " "

552

553

d1=self["tree"][res]["d1"]["index"]

554

d2=self["tree"][res]["d2"]["index"]

555

556

# For the daughters, the mass is either generate (intermediate leg + BW mode on)

557

# or set to the pole mass (external leg or BW mode off)

558

# If ran=0, just read the value from the previous generation of momenta

559

# (this is used for reshuffling purposes)

560

if d1>0 or not BW_effects :

561

mB=float(self.banner.get('param_card', 'mass',

562

abs(self.label2pid[self["tree"][res]["d1"]["label"]])).value)

563

elif ran==0: # reshuffling phase

564

mB=self["tree"][res]["d1"]["mass"]

565

else:

566

pid=self.label2pid[self["tree"][res]["d1"]["label"]]

567

# NOTE: here pole and width are normalized by 4.0*mB**2,

568

# Just a convention

569

pole=0.25 #pid2mass[pid]**2/mA**2

570

width=pid2width[pid]*pid2mass[pid]/(4.0*pid2mass[pid]**2) #/mA**2

571

mB, jac=self.transpole(pole,width)

572

mB=math.sqrt(mB*4.0*pid2mass[pid]**2)

573

# record the mass for the reshuffling phase,

574

# in case the point passes the reweighting creteria

575

self["tree"][res]["d1"]["mass"]=mB

576

# update the weigth of the phase-space point

577

weight=weight*jac

578

579

if d2>0 or not BW_effects:

580

mC=float(self.banner.get('param_card', 'mass',

581

abs(self.label2pid[self["tree"][res]["d2"]["label"]])).value)

582

elif ran==0:

583

mC=self["tree"][res]["d2"]["mass"]

584

else:

585

pid=self.label2pid[self["tree"][res]["d2"]["label"]]

586

# NOTE: here pole and width are normalized by 4.0*mC**2,

587

# Just a convention

588

pole=0.25 #pid2mass[pid]**2/mA**2

589

width=pid2width[pid]*pid2mass[pid]/(4.0*pid2mass[pid]**2) #mA**2

590

mC, jac=self.transpole(pole,width)

591

mC=math.sqrt(mC*4.0*pid2mass[pid]**2)

592

# record the mass for the reshuffling phase,

593

# in case the point passes the reweighting creteria

594

self["tree"][res]["d2"]["mass"]=mC

595

# update the weigth of the phase-space point

596

weight=weight*jac

597

598

599

if (mA<mB+mC):

600

logger.debug('mA<mB+mC in generate_momenta')

601

logger.debug('mA = %s' % mA)

602

return 0, 0 # If that happens, throw away the DC phase-space point ...

603

# I don't expect this to be inefficient, since there is a BW cut

604

605

if ran==1:

606

decay_mom=generate_2body_decay(index2mom[res]["momentum"],mA, mB,mC)

607

# record the angles for the reshuffling phase,

608

# in case the point passes the reweighting creteria

609

self["tree"][res]["costh"]=decay_mom.costh

610

self["tree"][res]["sinth"]=decay_mom.sinth

611

self["tree"][res]["cosphi"]=decay_mom.cosphi

612

self["tree"][res]["sinphi"]=decay_mom.sinphi

613

else:

614

# we are in the reshuffling phase,

615

# so we read the angles that have been stored from the

616

# previous phase-space point generation

617

costh=self["tree"][res]["costh"]

618

sinth=self["tree"][res]["sinth"]

619

cosphi=self["tree"][res]["cosphi"]

620

sinphi=self["tree"][res]["sinphi"]

621

decay_mom=generate_2body_decay(index2mom[res]["momentum"],mA, mB,mC,\

622

costh_val=costh, sinth_val=sinth, cosphi_val=cosphi, \

623

sinphi_val=sinphi)

624

625

# record the momenta for later use

626

index2mom[self["tree"][res]["d1"]["index"]]={}

627

index2mom[self["tree"][res]["d1"]["index"]]["momentum"]=decay_mom.momd1

628

index2mom[self["tree"][res]["d1"]["index"]]["pid"]=self.label2pid[self["tree"]\

629

[res]["d1"]["label"]]

630

631

index2mom[self["tree"][res]["d2"]["index"]]={}

632

index2mom[self["tree"][res]["d2"]["index"]]["momentum"]=decay_mom.momd2

633

index2mom[self["tree"][res]["d2"]["index"]]["pid"]=self.label2pid[self["tree"]\

634

[res]["d2"]["label"]]

635

636

if (self["tree"][res]["d1"]["index"]>0):

637

index2mom[self["tree"][res]["d1"]["index"]]["status"]=1

638

else:

639

index2mom[self["tree"][res]["d1"]["index"]]["status"]=2

640

if (self["tree"][res]["d2"]["index"]>0):

641

index2mom[self["tree"][res]["d2"]["index"]]["status"]=1

642

else:

643

index2mom[self["tree"][res]["d2"]["index"]]["status"]=2

644

645

return index2mom, weight

646

647

def find_tree(self,list_decays):

648

"""

649

record the topology of the decay chain in suitable variables

650

This is roughly the equivalent of the configs.inc file in madevent

651

"""

652

self["tree"]={}

653

nexternal=0

654

for mother in range(-1, -len(list_decays)-1, -1):

655

self["tree"][mother]={}

656

self["tree"][mother]["d1"]={}

657

self["tree"][mother]["d2"]={}

658

659

# print "filling the tree"

660

# print "res "+str(mother)

661

dc_nb=-(mother)-1

662

daughters=list_decays[dc_nb][list_decays[dc_nb].find(">")+1:].split()

663

self["tree"][mother]["d1"]["label"]=daughters[0]

664

self["tree"][mother]["d2"]["label"]=daughters[1]

665

666

if self["m_label2index"].has_key(daughters[0]):

667

self["tree"][mother]["d1"]["index"]=self["m_label2index"][daughters[0]]

668

else:

669

nexternal=nexternal+1

670

self["tree"][mother]["d1"]["index"]=nexternal

671

if self["m_label2index"].has_key(daughters[1]):

672

self["tree"][mother]["d2"]["index"]=self["m_label2index"][daughters[1]]

673

else:

674

nexternal=nexternal+1

675

self["tree"][mother]["d2"]["index"]=nexternal

676

return nexternal

677

678

def print_branch(self):

679

"""Print the decay chain structure (for debugging purposes)"""

680

length=len(self["tree"])

681

for res in range(-1,-length-1, -1):

682

logger.info('Decay '+str(res))

683

#print "Mother: "+self["tree"][res]

684

logger.info( 'd1: '+str(self["tree"][res]["d1"]["label"])+\

685

' '+str(self["tree"][res]["d1"]["index"]))

686

logger.info('d2: '+str(self["tree"][res]["d2"]["label"])+\

687

' '+str(self["tree"][res]["d2"]["index"]))

688

689

class momentum:

690

"""A class to handel 4-vectors and the associated operations """

691

def __init__(self,E,px,py,pz):

692

self.px=px

693

self.py=py

694

self.pz=pz

695

self.E=E

696

self.mod2=px**2+py**2+pz**2

697

self.sq=E**2-self.mod2

698

if (self.sq) > self.mod2*1e-10:

699

self.m=math.sqrt(self.sq)

700

# in case we get a very small negative value, set the mass to zero

701

elif (self.sq) > -self.mod2*1e-10: self.m=0.0

702

703

def dot3(self,q):

704

""" return |p|^2 (spatial components only) """

705

return self.px*q.px+self.py*q.py+self.pz*q.pz

706

707

def dot(self,q):

708

""" Minkowski inner product """

709

return self.E*q.E-self.px*q.px-self.py*q.py-self.pz*q.pz

710

711

def subtract(self,q):

712

tot=momentum(self.E-q.E,self.px-q.px,self.py-q.py,self.pz-q.pz)

713

return tot

714

715

def add(self,q):

716

tot=momentum(self.E+q.E,self.px+q.px,self.py+q.py,self.pz+q.pz)

717

return tot

718

719

def nice_string(self):

720

return str(self.E)+" "+str(self.px)+" "+str(self.py)+" "+str(self.pz)

721

722

def boost(self, q):

723

""" boost a vector from a frame where q is at rest to a frame where q is given

724

This routine has been taken from HELAS

725

"""

726

qq = q.mod2

727

728

if (qq > 1E-10*abs(q.E)):

729

pq=self.dot3(q)

730

m=q.m

731

#if (abs(m-self.mA)>1e-6): print "warning: wrong mass"

732

lf=((q.E-m)*pq/qq+self.E)/m

733

pboost=momentum((self.E*q.E+pq)/m, self.px+q.px*lf,\

734

self.py+q.py*lf,self.pz+q.pz*lf)

735

else:

736

pboost=momentum(self.E,self.px,self.py,self.pz)

737

738

return pboost

739

740

def copy(self):

741

copy_mom=momentum(self.E,self.px,self.py,self.pz)

742

return copy_mom

743

744

def invrot(self,q):

745

# inverse of the "rot" operation

746

747

ppE=self.E

748

qt2 = (q.px)**2 + (q.py)**2

749

750

if(qt2==0.0):

751

if ( q.pz>0 ):

752

ppx = self.px

753

ppy = self.py

754

ppz = self.pz

755

else:

756

ppx = -self.px

757

ppy = -self.py

758

ppz = -self.pz

759

else:

760

qq = math.sqrt(qt2+q.pz**2)

761

qt=math.sqrt(qt2)

762

ppy=-q.py/qt*self.px+q.px/qt*self.py

763

if (q.pz==0):

764

ppx=-qq/qt*self.pz

765

if (q.py!=0):

766

ppz=(self.py-q.py*q.pz/qq/qt-q.px/qt*ppy)*qq/q.py

767

else:

768

ppz=(self.px-q.px*q.pz/qq/qt*ppx+q.py/qt*ppy)*qq/q.px

769

else:

770

if (q.py!=0):

771

ppz=(qt**2*self.py+q.py*q.pz*self.pz-q.px*qt*ppy)/qq/q.py

772

else:

773

ppz=(q.px*self.px+q.pz*self.pz)/qq

774

ppx=(-self.pz+q.pz/qq*ppz)*qq/qt

775

pp=momentum(ppE,ppx,ppy,ppz)

776

return pp

777

778

779

def rot(self, q):

780

""" rotate the spatial components of the vector from a frame where q is

781

aligned with the z axis to a frame where the direction of q is specified

782

as an argument

783

Taken from HELAS

784

"""

785

protE =self.E

786

qt2 = (q.px)**2 + (q.py)**2

787

788

if(qt2==0.0):

789

if ( q.pz>0 ):

790

protx = self.px

791

proty = self.py

792

protz = self.pz

793

else:

794

protx = -self.px

795

proty = -self.py

796

protz = -self.pz

797

else:

798

qq = math.sqrt(qt2+q.pz**2)

799

qt = math.sqrt(qt2)

800

protx = q.px*q.pz/qq/qt*self.px -q.py/qt*self.py +q.px/qq*self.pz

801

proty = q.py*q.pz/qq/qt*self.px +q.px/qt*self.py +q.py/qq*self.pz

802

protz = -qt/qq*self.px + q.pz/qq*self.pz

803

804

prot=momentum(protE,protx,proty,protz)

805

return prot

806

807

808

class generate_2body_decay:

809

"""generate momenta for a generic A > B + C decay """

810

811

def __init__(self,p,mA,mB,mC, costh_val=None, sinth_val=None, cosphi_val=None, sinphi_val=None):

812

""" Generate the momentum of B and C in the decay A -> B+C

813

If the angles are given, use them to reconstruct the momenta of B, C

814

in the rest fram of A.

815

If the routine is called without (costh_val, ...), then generate

816

cos(theta) and phi randomly (flat distr.) in the rest frame of A

817

Finally, boost the momenta of B and C in the frame where A has

818

momentum p

819

"""

820

821

self.mA=mA

822

self.mB=mB

823

self.mC=mC

824

825

pmod=self.lambda_fct()/(2.0 * self.mA)

826

if not costh_val:

827

costh=2.0*random.random()-1.0

828

sinth=math.sqrt(1-costh**2)

829

else:

830

costh=costh_val

831

sinth=sinth_val

832

833

if not cosphi_val:

834

phi=random.random()*2.0*math.pi

835

sinphi=math.sin(phi)

836

cosphi=math.cos(phi)

837

else:

838

sinphi=sinphi_val

839

cosphi=cosphi_val

840

841

energyB=math.sqrt(pmod**2+mB**2)

842

energyC=math.sqrt(pmod**2+mC**2)

843

pBrest=momentum(energyB, pmod*cosphi*sinth,pmod*sinphi*sinth, pmod*costh)

844

pCrest=momentum(energyC,-pmod*cosphi*sinth,-pmod*sinphi*sinth, -pmod*costh)

845

self.momd1=pBrest.boost(p)

846

self.momd2=pCrest.boost(p)

847

848

# record costh and phi for later use

849

self.costh=costh

850

self.sinth=sinth

851

self.cosphi=cosphi

852

self.sinphi=sinphi

853

854

def lambda_fct(self):

855

""" The usual lambda function involved in 2-body decay """

856

lam=self.mA**4+self.mB**4+self.mC**4

857

lam=lam-2.0*self.mA**2*self.mB**2-2.0*self.mA**2*self.mC**2\

858

-2.0*self.mC**2*self.mB**2

859

# if lam<0:

860

# print self.mA

861

# print self.mB

862

# print self.mC

863

return math.sqrt(lam)

864

865

866

867

868

class production_topo(dict):

869

""" A dictionnary to record information about a given topology of a production event

870

871

self["branchings"] is a list of the branchings defining the topology (see class branching)

872

self["get_mass2"] is a dictionnary {index -> mass**2 of the corresponding particle}

873

self["get_momentum"] is a dictionnary {index -> momentum of the corresponding particle}

874

self["get_id"] is a dictionnary {index -> pid the corresponding particle}

875

876

Note: index= "madgraph-like" numerotation of the particles

877

"""

878

879

def __init__(self):

880

""" Initialise the dictionaries+list used later on to record the information

881

about the topology of a production event.

882

Note that self["branchings"] is a list,

883

as it makes it easier to scan the topology in the ascendent order

884

"""

885

self["branchings"]=[]

886

self["get_mass2"]={}

887

self["get_momentum"]={}

888

self["get_id"]={}

889

890

def add_one_branching(self,index_propa, index_d1,index_d2,type_propa):

891

""" add the information of one splitting in the topology """

892

branch=branching(index_propa, index_d1,index_d2,type_propa)

893

self["branchings"].append(branch)

894

895

896

def topo2event(self,event,to_decay):

897

"""This routine is typically called and the end of the reshuffling phase.

898

The momenta in the topology were reshuffled in a previous step, and now they are copied

899

back to the production event in this routine

900

"""

901

# start with external legs

902

for part in range(1,len(event.particle)+1):

903

event.particle[part]["momentum"]=self["get_momentum"][part].copy()

904

if part in to_decay:

905

if event.particle[part]["momentum"].m < 1.0e-3:

906

logger.debug('Decaying particle with a mass of less than 1 MeV in topo2event')

907

# print part

908

# print self["get_momentum"][part].nice_string()

909

# print event.particle[part]["momentum"].nice_string()

910

911

def print_topo(self):

912

"""Print the structure of the topology """

913

for branch in self["branchings"]:

914

d1=branch["index_d1"]

915

d2=branch["index_d2"]

916

propa=branch["index_propa"]

917

line=str(propa)+" > "

918

line+=str(d1)+" + "

919

line+=str(d2)+" , type="

920

line+=branch["type"]

921

print line

922

# try:

923

# print "momentum propa"

924

# print self["get_momentum"][propa].nice_string()

925

# print "P^2: "+str(self["get_mass2"][propa])

926

# print "root: "+str(math.sqrt(abs(self["get_mass2"][propa])))

927

# print "momentum d1"

928

# print self["get_momentum"][d1].nice_string()

929

# print "P^2: "+str(self["get_mass2"][d1])

930

# print "root: "+str(math.sqrt(abs(self["get_mass2"][d1])))

931

# print "momentum d2"

932

# print self["get_momentum"][d2].nice_string()

933

# print "P^2: "+str(self["get_mass2"][d2])

934

# print "root: "+str(math.sqrt(abs(self["get_mass2"][d2])))

935

# except Exception, error:

936

# print "topology not yet dressed"

937

938

def dress_topo_from_event(self,event,to_decay):

939

""" event has been read from the production events file,

940

use these momenta to dress the topology

941

"""

942

943

# start with external legs

944

for part in range(1,len(event.particle)+1):

945

self["get_momentum"][part]=event.particle[part]["momentum"].copy()

946

self["get_mass2"][part]=event.particle[part]["mass"]**2

947

self["get_id"][part]=event.particle[part]["pid"]

948

if part in to_decay :

949

if self["get_momentum"][part].m<1e-3:

950

logger.debug\

951

('decaying particle with m < 1MeV in dress_topo_from_event (1)')

952

if self["get_mass2"][part]<1e-3:

953

logger.debug\

954

('decaying particle with m < 1MeV in dress_topo_from_event (2)')

955

956

# now fill also intermediate legs

957

# Don't care about the pid of intermediate legs

958

for branch in self["branchings"]:

959

part=branch["index_propa"]

960

if branch["type"]=="s":

961

mom_propa=self["get_momentum"][branch["index_d1"]].add(self["get_momentum"][branch["index_d2"]])

962

elif branch["type"]=="t":

963

mom_propa=self["get_momentum"][branch["index_d1"]].subtract(self["get_momentum"][branch["index_d2"]])

964

self["get_momentum"][part]=mom_propa

965

self["get_mass2"][part]=mom_propa.sq

966

967

# Also record shat and rapidity, since the initial momenta will also be reshuffled

968

#

969

p1=self["get_momentum"][1].copy()

970

p2=self["get_momentum"][2].copy()

971

ptot=p1.add(p2)

972

self["shat"]=ptot.sq

973

self["rapidity"]=0.5*math.log((ptot.E+ptot.pz)/(ptot.E-ptot.pz))

974

975

def reshuffle_momenta(self):

976

"""

977

At this stage,

978

- the topo should be dressed with momenta.

979

- the angles should be already extracted.

980

- the masses should be corrected.

981

This routine scan all branchings:

982

- first consider the t-branchings, go to the appropriate frame,

983

modify the three-momenta to account for the corrected masses,

984

- then consider the s-branchings, go to the approriate frame,

985

rescale the three-momenta to account for the corrected masses.

986

"""

987

988

# step number one: need to check if p^2 of each propa

989

# is ok with the new set of masses ...

990

# we first check this for all the s-branchings

991

992

# Close to threshold, problems may occur:

993

# e.g. in A>B+C, one may have mA>mB+mC

994

# Currently, if this happens, I throw away the point

995

# but I am not sure what is the best prescription here ...

996

997

for branch in self["branchings"]:

998

999

# no need to consider t-branching here:

1000

if branch["type"]!="s":

1001

continue

1002

d1=branch["index_d1"]

1003

d2=branch["index_d2"]

1004

propa=branch["index_propa"]

1005

1006

MA=math.sqrt(self["get_mass2"][propa])

1007

MB=math.sqrt(self["get_mass2"][d1])

1008

# print self["get_mass2"][d2]

1009

MC=math.sqrt(self["get_mass2"][d2])

1010

if MA < MB + MC:

1011

# print "WARNING: s-channel propagator with too low p^2 "

1012

# print "in branch "+str(iter)

1013

# print "MA = "+str(MA)

1014

# print "MB = "+str(MB)

1015

# print "MC = "+str(MC)

1016

# print "throw away the point"

1017

1018

# print "increase the value of p^2"

1019

# self["get_mass2"][propa]=(MC+MB+iota)**2

1020

# if iter==len(self["branchings"])-1: # if last branching, needs to

1021

# self["shat"]=self["get_mass2"][propa] # set shat to the new value of MA**2

1022

return 0

1023

1024

1025

1026

# then loop over all t-channels

1027

# and re-genenate the "d2" daughter in each of these branchings

1028

got_a_t_branching=0

1029

for nu, branch in enumerate(self["branchings"]):

1030

1031

# no need to consider the last branching in this loop:

1032

if(nu==len(self["branchings"])-1): break

1033

1034

# no need to scan the s-branching now

1035

if branch["type"]!="t":

1036

continue

1037

1038

got_a_t_branching=1

1039

# t-channel sequence: A+B > 1 + 2, r= pa-p1

1040

ida=branch["index_d1"]

1041

idb=2

1042

id1=branch["index_d2"]

1043

res=branch["index_propa"]

1044

# go to the rest frame of A+B

1045

# set momenta A, B, 1

1046

pa = self["get_momentum"][ida]

1047

pb = self["get_momentum"][idb]

1048

p1 = self["get_momentum"][id1]

1049

# set masses

1050

ma2=self["get_mass2"][ida]

1051

if (self["get_mass2"][id1]>=0):

1052

#print "m1^2, t-branch "+str(iter)

1053

#print self["get_mass2"][id1]

1054

m1=math.sqrt(self["get_mass2"][id1])

1055

else:

1056

# print "WARNING: m1^2 is negative for t-branching "+str(iter)

1057

# print self["get_mass2"][id1]

1058

# print "throw away the point"

1059

return 0

1060

m2=branch["m2"]

1061

t=self["get_mass2"][res]

1062

1063

# express momenta p1 and pa in A+B CMS system

1064

pboost=self["get_momentum"][2].add(self["get_momentum"][branch["index_d1"]])

1065

pboost.px=-pboost.px

1066

pboost.py=-pboost.py

1067

pboost.pz=-pboost.pz

1068

# p1_cms=p1.boost(pboost)

1069

pa_cms=pa.boost(pboost)

1070

1071

# determine the magnitude of p1 in the cms frame

1072

Esum=pboost.sq

1073

1074

if Esum>0 :

1075

Esum=math.sqrt(Esum)

1076

else:

1077

# print "WARNING: (pa+pb)^2 is negative for t-branching "

1078

return 0

1079

md2=(m1+m2)*(m1-m2)

1080

ed=md2/Esum

1081

if (m1*m2==0) :

1082

pp=(Esum-abs(ed))*0.5

1083

else:

1084

pp=(md2/Esum)**2-2.0*(m1**2+m2**2)+Esum**2

1085

if pp>0 :

1086

pp=0.5*math.sqrt(pp)

1087

else:

1088

# print "WARNING: cannot get the momentum of p1 in t-branching "+str(iter)

1089

# print "pp is negative : "+ str(pp)

1090

# print "m1: "+ str(m1)

1091

# print "m2: "+ str(m2)

1092

# print "id1: "+ str(id1)

1093

# print "throw away the point"

1094

return 0

1095

1096

# Now evaluate p1

1097

E_acms=pa_cms.E

1098

p_acms=math.sqrt(pa_cms.mod2)

1099

1100

p1E=(Esum+ed)*0.5

1101

if p1E < m1:

1102

# logger.warning('E1 is smaller than m1 in t-branching')

1103

# logger.warning('Try to reshuffle the momenta once more')

1104

return 0

1105

p1z=-(m1*m1+ma2-t-2.0*p1E*E_acms)/(2.0*p_acms)

1106

ptsq=pp*pp-p1z*p1z

1107

1108

if (ptsq<0 ):

1109

# if pT=0 to begin with, one can get p1z slightly larger

1110

# than pp due to numerical uncertainties.

1111

# In that case, just change slightly the invariant t

1112

if (-ptsq/(pp*pp) <1e-6 ):

1113

oldt=t

1114

p1z=pp

1115

pt=0.0

1116

t=m1*m1+ma2-2.0*p1E*E_acms+2.0*p_acms*p1z

1117

diff_t=abs((t-oldt)/t)*100

1118

if (diff_t>2.0):

1119

logger.warning('t invariant was changed by '+str(diff_t)+' percents')

1120

else:

1121

# print "WARNING: |p|^2 is smaller than p1z^2 in t-branching "+str(iter)

1122

# print "|p| : "+str(pp)

1123

# print "pz : "+str(p1z)

1124

# print "throw away the point"

1125

# print "Set pz^2=|p|^2 and recompute t"

1126

# print "previous t:"+str(-math.sqrt(abs(t)))+"^2"

1127

# p1z=pp

1128

# pt=0.0

1129

# t=m1*m1+ma2-2.0*p1E*E_acms+2.0*p_acms*p1z

1130

# print "new t:"+str(-math.sqrt(abs(t)))+"^2"

1131

return 0

1132

else:

1133

pt=math.sqrt(pp*pp-p1z*p1z)

1134

1135

p1x=pt*branch["cosphi"]

1136

p1y=pt*branch["sinphi"]

1137

1138

p1=momentum(p1E,p1x,p1y,p1z)

1139

1140

p1=p1.rot(pa_cms)

1141

pboost.px=-pboost.px

1142

pboost.py=-pboost.py

1143

pboost.pz=-pboost.pz

1144

p1=p1.boost(pboost)

1145

pr=pa.subtract(p1)

1146

# print " p1 is "

1147

# print p1.nice_string()

1148

#print " pr is "

1149

#print pr.nice_string()

1150

p2=(pa.add(pb)).subtract(p1)

1151

# print " p2 is "

1152

# print p2.nice_string()

1153

# now update momentum

1154

self["get_momentum"][id1]=p1.copy()

1155

self["get_momentum"][res]=pr.copy()

1156

1157

# after we have looped over all t-branchings,

1158

# p2 can be identified with the momentum of the second daughter of the last branching

1159

1160

if got_a_t_branching==1:

1161

pid=self["branchings"][-1]["index_d2"]

1162

self["get_momentum"][pid]=p2.copy()

1163

#else: it means that there were no t-channel at all

1164

# last branching should be associated with shat

1165

# note that the initial momenta will be reshuffled

1166

# at the end of this routine

1167

1168

1169

# Now we can loop over all the s-channel branchings.

1170

# Need to start at the end of the list of branching

1171

for branch in reversed(self["branchings"]):

1172

if branch["type"]!="s":

1173

continue

1174

d1=branch["index_d1"]

1175

d2=branch["index_d2"]

1176

propa=branch["index_propa"]

1177

del self["get_momentum"][d1]

1178

del self["get_momentum"][d2]

1179

mA=math.sqrt(self["get_mass2"][propa])

1180

mB=math.sqrt(self["get_mass2"][d1])

1181

mC=math.sqrt(self["get_mass2"][d2])

1182

mom=self["get_momentum"][propa]

1183

costh=branch["costheta"]

1184

sinth=branch["sintheta"]

1185

cosphi=branch["cosphi"]

1186

sinphi=branch["sinphi"]

1187

decay2body=generate_2body_decay(mom, mA,mB,mC, \

1188

costh_val=costh, sinth_val=sinth, \

1189

cosphi_val=cosphi, sinphi_val=sinphi)

1190

self["get_momentum"][d1]=decay2body.momd1.copy()

1191

self["get_momentum"][d2]=decay2body.momd2.copy()

1192

1193

# Need a special treatment for the 2 > 1 processes:

1194

if len(self["get_mass2"])==3:

1195

self["shat"]=self["get_mass2"][3]

1196

1197

# Finally, compute the initial momenta

1198

# First generate initial momenta in the CMS frame

1199

# Then boost

1200

mB=self["get_mass2"][1]

1201

mC=self["get_mass2"][2]

1202

if mB>0:

1203

mB=math.sqrt(mB)

1204

else:

1205

mB=0.0

1206

if mC>0:

1207

mC=math.sqrt(mC)

1208

else:

1209

mC=0.0

1210

mA=math.sqrt(self["shat"])

1211

Etot=mA*math.cosh(self["rapidity"])

1212

pztot=mA*math.sinh(self["rapidity"])

1213

ptot=momentum(Etot,0.0,0.0,pztot)

1214

1215

decay2body=generate_2body_decay(ptot,mA,mB,mC, \

1216

costh_val=1.0, sinth_val=0.0, \

1217

cosphi_val=0.0, sinphi_val=1.0)

1218

1219

self["get_momentum"][1]=decay2body.momd1

1220

self["get_momentum"][2]=decay2body.momd2

1221

1222

# Need a special treatment for the 2 > 1 processes:

1223

if len(self["get_momentum"])==3:

1224

self["get_momentum"][3]=momentum(Etot,0.0,0.0,pztot)

1225

1226

return 1

1227

1228

def extract_angles(self):

1229

""" the topo should be dressed with momenta at this stage.

1230

Now: extract the angles characterizing each branching.

1231

For t-channel, the equivalent of cos(theta) is the mass of p2

1232

"""

1233

1234

for nu, branch in enumerate(self["branchings"]):

1235

if branch["type"]=="s":

1236

# we have the decay A > B + C

1237

# go to the rest frame of the decaying particle A, and extract cos theta and phi

1238

propa=branch["index_propa"]

1239

# d1=branch["index_d1"]

1240

d2=branch["index_d1"]

1241

pboost=self["get_momentum"][propa].copy()

1242

pboost.px=-pboost.px

1243

pboost.py=-pboost.py

1244

pboost.pz=-pboost.pz

1245

# pb_cms=self["get_momentum"][d1].boost(pboost)

1246

pc_cms=self["get_momentum"][d2].boost(pboost)

1247

mod_pc_cms=math.sqrt(pc_cms.mod2)

1248

branch["costheta"]=pc_cms.pz/mod_pc_cms

1249

branch["sintheta"]=math.sqrt(1.0-branch["costheta"]**2)

1250

branch["cosphi"]=pc_cms.px/mod_pc_cms/branch["sintheta"]

1251

branch["sinphi"]=pc_cms.py/mod_pc_cms/branch["sintheta"]

1252

1253

if branch["type"]=="t":

1254

# we have a t-channel decay A + B > 1 + 2

1255

# go to the rest frame of A+B, and extract phi

1256

1257

# set momenta A, B, 1

1258

pa = self["get_momentum"][branch["index_d1"]]

1259

pb = self["get_momentum"][2]

1260

p1 = self["get_momentum"][branch["index_d2"]]

1261

p2=(pa.add(pb)).subtract(p1)

1262

1263

if (nu==len(self["branchings"])-1):

1264

# last t-channel branching, p2 should be zero

1265

check=p2.E**2+p2.mod2

1266

if check>1e-3:

1267

logger.warning('p2 in the last t-branching is not zero')

1268

# If last t-channel branching, there is no angles to extract

1269

continue

1270

elif (nu==len(self["branchings"])-2):

1271

# here m2 corresponds to the mass of the second daughter in the last splitting

1272

part=self["branchings"][-1]["index_d2"]

1273

if self["get_mass2"][part]<0.0 :

1274

logger.warning('negative mass for particle '+str(part))

1275

logger.warning( self["get_mass2"][part])

1276

branch["m2"]=math.sqrt(self["get_mass2"][part])

1277

elif p2.sq <0 and abs(p2.E)>1e-2:

1278

logger.warning('GET A NEGATIVE M2 MASS IN THE T-BRANCHING '+str(iter))

1279

logger.warning( '(ROUTINE EXTRACT_ANGLES)')

1280

logger.warning( p2.nice_string() )

1281

logger.warning( p2.sq )

1282

else:

1283

branch["m2"]=p2.m

1284

1285

# express momenta p1 and pa in A+B CMS system

1286

pboost=self["get_momentum"][2].add(self["get_momentum"][branch["index_d1"]])

1287

pboost.px=-pboost.px

1288

pboost.py=-pboost.py

1289

pboost.pz=-pboost.pz

1290

1291

if pboost.m<1e-6:

1292

logger.warning('Warning: m=0 in T-BRANCHING '+str(iter))

1293

logger.warning(' pboost.nice_string()')

1294

1295

p1_cms=p1.boost(pboost)

1296

pa_cms=pa.boost(pboost)

1297

1298

# E_acms=pa.E

1299

# mod_acms=math.sqrt(pa_cms.mod2)

1300

#

1301

# now need to go to the frame where pa is aligned along with the z-axis

1302

p_rot=pa_cms.copy()

1303

1304

p1_cmsrot=p1_cms.invrot(p_rot)

1305

# pa_cmsrot=pa_cms.invrot(p_rot)

1306

# now extract cosphi, sinphi

1307

pt=math.sqrt(p1_cmsrot.px**2+p1_cmsrot.py**2)

1308

if (pt>p1_cms.E*10e-10):

1309

branch["cosphi"]=p1_cmsrot.px/pt

1310

branch["sinphi"]=p1_cmsrot.py/pt

1311

else:

1312

branch["cosphi"]=1

1313

branch["sinphi"]=0.0

1314

1315

1316

class branching(dict):

1317

""" A dictionnary to record information about a given branching in a production event

1318

self["type"] is the type of the branching , either "s" for s-channel or "t" for t-channel

1319

self["invariant"] is a real number with the value of p^2 associated with the branching

1320

self["index_d1"] is the mg index of the first daughter

1321

self["index_d2"] is the mg index of the second daughter

1322

self["index_propa"] is the mg index of the propagator

1323

"""

1324

1325

def __init__(self, index_propa, index_d1,index_d2, s_or_t):

1326

self["index_d1"]=index_d1

1327

self["index_d2"]=index_d2

1328

self["index_propa"]=index_propa

1329

self["type"]=s_or_t

1330

1331

class width_estimate:

1332

"""All methods used to calculate branching fractions"""

1333

1334

def __init__(self,resonances,path_me,pid2label_dic,banner,base_model):

1335

1336

self.resonances=resonances

1337

self.path_me=path_me

1338

self.pid2label=pid2label_dic

1339

self.label2pid = self.pid2label

1340

self.model=banner.get('model')

1341

self.banner = banner

1342

self.base_model=base_model

1343

1344

def update_branch(self,branches,to_add):

1345

""" complete the definition of the branch by appending each element of to_add"""

1346

newbranches={}

1347

1348

for item1 in branches.keys():

1349

for item2 in to_add.keys():

1350

tag=item1+item2

1351

newbranches[tag]={}

1352

newbranches[tag]['config']=branches[item1]['config']+to_add[item2]['config']

1353

newbranches[tag]['br']=branches[item1]['br']*to_add[item2]['br']

1354

1355

return newbranches

1356

1357

def get_BR_for_each_decay(self, decay_processes, multiparticles):

1358

""" get the list for possible decays & the associated branching fraction """

1359

1360

model = self.model

1361

base_model = self.base_model

1362

pid2label = self.pid2label

1363

1364

ponctuation=[',','>',')','(']

1365

new_decay_processes={}

1366

1367

for part in decay_processes.keys():

1368

pos_symbol=-1

1369

branch_list=decay_processes[part].split()

1370

new_decay_processes[part]={}

1371

new_decay_processes[part]['']={}

1372

new_decay_processes[part]['']['config']=""

1373

new_decay_processes[part]['']['br']=1.0

1374

1375

initial=""

1376

final=[]

1377

for index, item in enumerate(branch_list):

1378

# First get the symbol at the next position

1379

if index<len(branch_list)-1:

1380

next_symbol=branch_list[index+1]

1381

else:

1382

next_symbol=''

1383

# Then handle the symbol item case by case

1384

if next_symbol=='>': # case1: we have a particle initiating a branching

1385

initial=item

1386

if item not in [ particle['name'] for particle in base_model['particles'] ] \

1387

and item not in [ particle['antiname'] for particle in base_model['particles'] ]:

1388

raise Exception, "No particle "+item+ " in the model "+model

1389

continue

1390

elif item=='>': continue # case 2: we have the > symbole

1391

elif item not in ponctuation : # case 3: we have a particle originating from a branching

1392

final.append(item)

1393

if next_symbol=='' or next_symbol in ponctuation:

1394

#end of a splitting, verify that it exists

1395

if initial not in self.br.keys():

1396

logger.debug('Branching fractions of particle '+initial+' are unknown')

1397

return 0

1398

if len(final)>2:

1399

raise Exception, 'splittings different from A > B +C are currently not implemented '

1400

1401

if final[0] in multiparticles.keys():

1402

set_B=[pid2label[pid] for pid in multiparticles[final[0]]]

1403

else:

1404

if final[0] not in [ particle['name'] for particle in base_model['particles'] ] \

1405

and final[0] not in [ particle['antiname'] for particle in base_model['particles'] ]:

1406

raise Exception, "No particle "+item+ " in the model "

1407

set_B=[final[0]]

1408

if final[1] in multiparticles.keys():

1409

set_C=[pid2label[pid] for pid in multiparticles[final[1]]]

1410

else:

1411

if final[1] not in [ particle['name'] for particle in base_model['particles'] ] \

1412

and final[1] not in [ particle['antiname'] for particle in base_model['particles'] ]:

1413

raise Exception, "No particle "+item+ " in the model "+model

1414

set_C=[final[1]]

1415

1416

splittings={}

1417

counter=0

1418

for chan in range(len(self.br[initial])): # loop over all channels

1419

got_it=0

1420

for d1 in set_B:

1421

for d2 in set_C:

1422

if (d1==self.br[initial][chan]['daughters'][0] and \

1423

d2==self.br[initial][chan]['daughters'][1]) or \

1424

(d2==self.br[initial][chan]['daughters'][0] and \

1425

d1==self.br[initial][chan]['daughters'][1]):

1426

split=" "+initial+" > "+d1+" "+d2+" "

1427

counter+=1

1428

splittings['s'+str(counter)]={}

1429

splittings['s'+str(counter)]['config']=split

1430

splittings['s'+str(counter)]['br']=self.br[initial][chan]['br']

1431

got_it=1

1432

break # to avoid double counting in cases such as w+ > j j

1433

if got_it: break

1434

1435

if len(splittings)==0:

1436

logger.info('Branching '+initial+' > '+final[0]+' '+final[1])

1437

logger.info('is currently unknown')

1438

1439

return 0

1440

else:

1441

new_decay_processes[part]=self.update_branch(new_decay_processes[part],splittings)

1442

1443

inital=""

1444

final=[]

1445

1446

else: # case 4: ponctuation symbol outside a splitting

1447

# just append it to all the current branches

1448

fake_splitting={}

1449

fake_splitting['']={}

1450

fake_splitting['']['br']=1.0

1451

fake_splitting['']['config']=item

1452

new_decay_processes[part]=self.update_branch(new_decay_processes[part],fake_splitting)

1453

1454

return new_decay_processes

1455

1456

def print_branching_fractions(self):

1457

""" print a list of all known branching fractions"""

1458

1459

for res in self.br.keys():

1460

logger.info(' ')

1461

logger.info('decay channels for '+res+' :')

1462

logger.info(' BR d1 d2' )

1463

for decay in self.br[res]:

1464

bran = decay['br']

1465

d1 = decay['daughters'][0]

1466

d2 = decay['daughters'][1]

1467

logger.info(' %e %s %s ' % (bran, d1, d2) )

1468

logger.info(' ')

1469

1470

def print_partial_widths(self):

1471

""" print a list of all known partial widths"""

1472

1473

for res in self.br.keys():

1474

logger.info(' ')

1475

logger.info('decay channels for '+res+' :')

1476

logger.info(' width d1 d2' )

1477

1478

for chan, decay in enumerate(self.br[res]):

1479

width=self.br[res][chan]['width']

1480

d1=self.br[res][chan]['daughters'][0]

1481

d2=self.br[res][chan]['daughters'][1]

1482

logger.info(' %e %s %s ' % (width, d1, d2) )

1483

logger.info(' ')

1484

1485

1486

def extract_br_from_width_evaluation(self):

1487

""" use madgraph to generate me's for res > all all

1488

"""

1489

if os.path.isdir(pjoin(self.path_me,"width_calculator")):

1490

shutil.rmtree(pjoin(self.path_me,"width_calculator"))

1491

1492

assert not os.path.exists(pjoin(self.path_me, "width_calculator"))

1493

1494

path_me = self.path_me

1495

label2pid = self.label2pid

1496

# first build a set resonances with pid>0

1497

1498

particle_set=[]

1499

for part in self.resonances:

1500

if label2pid[part]>0: particle_set.append(part)

1501

for part in self.resonances:

1502

if label2pid[part]<0:

1503

pid_part=-label2pid[part]

1504

if self.pid2label[pid_part] not in particle_set:

1505

particle_set.append(self.pid2label[pid_part])

1506

1507

1508

commandline="import model %s\n" % self.model

1509

commandline+="generate %s > all all \n" % particle_set[0]

1510

commandline+= "set automatic_html_opening False --no_save\n"

1511

if len(particle_set)>1:

1512

for index in range(1,len(particle_set)):

1513

commandline+="add process %s > all all \n" % particle_set[index]

1514

1515

commandline += "output %s/width_calculator -f \n" % path_me

1516

1517

1518

aloha.loop_mode = False

1519

aloha.unitary_gauge = False

1520

cmd = Cmd.MasterCmd()

1521

for line in commandline.split('\n'):

1522

cmd.run_cmd(line)

1523

files.cp(pjoin(path_me, 'Cards', 'param_card.dat'),

1524

pjoin(path_me, 'width_calculator', 'Cards'))

1525

1526

cmd.run_cmd('launch -f')

1527

1528

#me_cmd = me_interface.MadEventCmd(pjoin(path_me,'width_calculator'))

1529

#me_cmd.exec_cmd('set automatic_html_opening False --no_save')

1530

1531

filename=pjoin(path_me,'width_calculator','Events','run_01','param_card.dat')

1532

# misc.sprint(pjoin(path_me,'width_calculator','Events','run_01','param_card.dat'))

1533

self.extract_br_from_card(filename)

1534

1535

def extract_br_for_antiparticle(self):

1536

'''

1537

for each channel with a specific br value,

1538

set the branching fraction of the complex conjugated channel

1539

to the same br value

1540

'''

1541

1542

label2pid = self.label2pid

1543

pid2label = self.label2pid

1544

1545

for res in self.br.keys():

1546

particle=self.base_model.get_particle(label2pid[res])

1547

if particle['self_antipart']:

1548

continue

1549

anti_res=pid2label[-label2pid[res]]

1550

self.br[anti_res] = []

1551

for chan, decay in enumerate(self.br[res]):

1552

self.br[anti_res].append({})

1553

bran=decay['br']

1554

d1=decay['daughters'][0]

1555

d2=decay['daughters'][1]

1556

d1bar=pid2label[-label2pid[d1]]

1557

d2bar=pid2label[-label2pid[d2]]

1558

self.br[anti_res][chan]['br']=bran

1559

self.br[anti_res][chan]['daughters']=[]

1560

self.br[anti_res][chan]['daughters'].append(d1bar)

1561

self.br[anti_res][chan]['daughters'].append(d2bar)

1562

if decay.has_key('width'):

1563

self.br[anti_res][chan]['width']=decay['width']

1564

1565

def launch_width_evaluation(self,resonances, model, mg5cmd):

1566

""" launch the calculation of the partial widths """

1567

1568

label2pid = self.label2pid

1569

pid2label = self.label2pid

1570

# first build a set resonances with pid>0

1571

# since compute_width cannot be used for particle with pid<0

1572

1573

particle_set=[]

1574

for part in resonances:

1575

pid_part = abs(label2pid[part])

1576

if pid_part not in particle_set:

1577

particle_set.append(pid_part)

1578

# erase old info

1579

del self.br

1580

self.br={}

1581

1582

1583

1584

argument = {'particles': particle_set,

1585

'input': pjoin(self.path_me, 'param_card.dat'),

1586

'output': pjoin(self.path_me, 'param_card.dat')}

1587

1588

me_interface.MadEventCmd.compute_widths(model, argument)

1589

self.extract_br_from_card(pjoin(self.path_me, 'param_card.dat'))

1590

return

1591

1592

1593

def extract_br_from_banner(self, banner):

1594

"""get the branching ratio from the banner object:

1595

for each resonance with label 'res', and for each channel with index i,

1596

returns a dictionary branching_fractions[res][i]

1597

with keys

1598

'daughters' : label of the daughters (only 2 body)

1599

'br' : value of the branching fraction"""

1600

1601

self.br = {}

1602

1603

# read the param_card internally to the banner

1604

if not hasattr(banner, 'param_card'):

1605

banner.charge_card('param_card')

1606

param_card = banner.param_card

1607

return self.extract_br_from_card(param_card)

1608

1609

def extract_br_from_card(self, param_card):

1610

"""get the branching ratio from the banner object:

1611

for each resonance with label 'res', and for each channel with index i,

1612

returns a dictionary branching_fractions[res][i]

1613

with keys

1614

'daughters' : label of the daughters (only 2 body)

1615

'br' : value of the branching fraction"""

1616

1617

if isinstance(param_card, str):

1618

import models.check_param_card as check_param_card

1619

param_card = check_param_card.ParamCard(param_card)

1620

1621

if 'decay' not in param_card or not hasattr(param_card['decay'], 'decay_table'):

1622

return self.br

1623

1624

for id, data in param_card['decay'].decay_table.items():

1625

label = self.pid2label[id]

1626

current = [] # tmp name for self.br[label]

1627

for parameter in data:

1628

if parameter.lhacode[0] == 2:

1629

d = [self.pid2label[pid] for pid in parameter.lhacode[1:]]

1630

current.append({'daughters':d, 'br': parameter.value})

1631

self.br[label] = current

1632

1633

self.extract_br_for_antiparticle()

1634

return self.br

1635

1636

1637

1638

1639

1640

1641

1642

1643

class decay_misc:

1644

"""class with various methods for the decay"""

1645

1646

1647

1648

1649

def decay_one_event(self,curr_event,decay_struct,pid2color_dico,\

1650

to_decay, pid2width,pid2mass,resonnances,BW_effects,ran=1):

1651

1652

# Consider the production event recorded in "curr_event", and decay it

1653

# according to the structure recoreded in "decay_struct".

1654

# If ran=1: random decay, phi and cos theta generated according to

1655

# a uniform distribution in the rest frame of the decaying particle

1656

# Ir ran=0 : use the previsously-generated angles and masses to get the momenta

1657

1658

1659

decayed_event=Event()

1660

decayed_event.event2mg={}

1661

1662

if ran==0: # reshuffling phase: the decayed event is about to be written, so we need

1663

# to record some information

1664

decayed_event.ievent=curr_event.ievent

1665

decayed_event.wgt=curr_event.wgt

1666

decayed_event.scale=curr_event.scale

1667

decayed_event.aqed=curr_event.aqed

1668

decayed_event.aqcd=curr_event.aqcd

1669

decayed_event.diese=curr_event.diese

1670

decayed_event.rwgt=curr_event.rwgt

1671

1672

part_number=0

1673

external=0

1674

maxcol=curr_event.max_col

1675

weight=1.0

1676

1677

#event2mg

1678

1679

for index in curr_event.event2mg.keys():

1680

if curr_event.event2mg[index]>0:

1681

part=curr_event.event2mg[index]

1682

if part in to_decay.keys():

1683

mom_init=curr_event.particle[part]["momentum"].copy()

1684

branch_id=to_decay[part]

1685

# sanity check

1686

if mom_init.m<1e-6:

1687

logger.debug('Decaying particle with mass less than 1e-6 GeV in decay_one_event')

1688

decay_products, jac=decay_struct[branch_id].generate_momenta(mom_init,\

1689

ran, pid2width,pid2mass,resonnances,BW_effects)

1690

1691

if ran==1:

1692

if decay_products==0: return 0, 0

1693

weight=weight*jac

1694

1695

# now we need to write the decay products in the event

1696

# follow the decay chain order, so that we can easily keep track of the mother index

1697

for res in range(-1,-len(decay_struct[branch_id]["tree"].keys())-1,-1):

1698

if (res==-1):

1699

part_number+=1

1700

mom=decay_products[res]["momentum"]

1701

pid=decay_products[res]["pid"]

1702

istup=2

1703

mothup1=1

1704

mothup2=2

1705

colup1=curr_event.particle[part]["colup1"]

1706

colup2=curr_event.particle[part]["colup2"]

1707

decay_products[res]["colup1"]=colup1

1708

decay_products[res]["colup2"]=colup2

1709

mass=mom.m

1710

helicity=0.

1711

decayed_event.particle[part_number]={"pid":pid,\

1712

"istup":istup,"mothup1":mothup1,"mothup2":mothup2,\

1713

"colup1":colup1,"colup2":colup2,"momentum":mom,\

1714

"mass":mass,"helicity":helicity}

1715

decayed_event.event2mg[part_number]=part_number

1716

# print part_number

1717

# print pid

1718

# print " "

1719

mothup1=part_number

1720

mothup2=part_number

1721

#

1722

# Extract color information so that we can write the color flow

1723

#

1724

colormother=pid2color_dico[decay_products[res]["pid"]]

1725

colord1=pid2color_dico[decay_products[decay_struct[branch_id]\

1726

["tree"][res]["d1"]["index"]]["pid"]]

1727

colord2=pid2color_dico[decay_products[decay_struct[branch_id]\

1728

["tree"][res]["d2"]["index"]]["pid"]]

1729

1730

colup1=decay_products[res]["colup1"]

1731

colup2=decay_products[res]["colup2"]

1732

1733

# now figure out what is the correct color flow informatio

1734

# Only consider 1,3, 3-bar and 8 color rep.

1735

# Normally, the color flow needs to be determined only

1736

# during the reshuffling phase, but it is currenlty assigned

1737

# for each "trial event"

1738

1739

if abs(colord2)==1:

1740

d1colup1=colup1

1741

d1colup2=colup2

1742

d2colup1=0

1743

d2colup2=0

1744

1745

if abs(colord1)==1:

1746

d2colup1=colup1

1747

d2colup2=colup2

1748

d1colup1=0

1749

d1colup2=0

1750

1751

if colord1==3 and colord2==-3 and colormother ==1:

1752

maxcol+=1

1753

d1colup1=maxcol

1754

d1colup2=0

1755

d2colup1=0

1756

d2colup2=maxcol

1757

1758

if colord1==-3 and colord2==3 and colormother ==1:

1759

maxcol+=1

1760

d1colup1=0

1761

d1colup2=maxcol

1762

d2colup1=maxcol

1763

d2colup2=0

1764

1765

if colord1==3 and colord2==8 and colormother ==3:

1766

maxcol+=1

1767

d2colup1=colup1

1768

d2colup2=maxcol

1769

d1colup1=maxcol

1770

d1colup2=0

1771

1772

if colord2==3 and colord1==8 and colormother ==3:

1773

maxcol+=1

1774

d1colup1=colup1

1775

d1colup2=maxcol

1776

d2colup1=maxcol

1777

d2colup2=0

1778

1779

if colord1==-3 and colord2==8 and colormother ==-3:

1780

maxcol+=1

1781

d2colup2=colup2

1782

d2colup1=maxcol

1783

d1colup2=maxcol

1784

d1colup1=0

1785

1786

if colord2==-3 and colord1==8 and colormother ==-3:

1787

maxcol+=1

1788

d1colup2=colup2

1789

d1colup1=maxcol

1790

d2colup2=maxcol

1791

d2colup1=0

1792

1793

part_number+=1

1794

mom=decay_products[decay_struct[branch_id]\

1795

["tree"][res]["d1"]["index"]]["momentum"]

1796

pid=decay_products[decay_struct[branch_id]\

1797

["tree"][res]["d1"]["index"]]["pid"]

1798

1799

1800

indexd1=decay_struct[branch_id]["tree"][res]["d1"]["index"]

1801

if ( indexd1>0):

1802

istup=1

1803

external+=1

1804

else:

1805

decay_products[indexd1]["colup1"]=d1colup1

1806

decay_products[indexd1]["colup2"]=d1colup2

1807

istup=2

1808

1809

mass=mom.m

1810

helicity=0.

1811

decayed_event.particle[part_number]={"pid":pid,\

1812

"istup":istup,"mothup1":mothup1,"mothup2":mothup2,\

1813

"colup1":d1colup1,"colup2":d1colup2,"momentum":mom,\

1814

"mass":mass,"helicity":helicity}

1815

decayed_event.event2mg[part_number]=part_number

1816

1817

part_number+=1

1818

mom=decay_products[decay_struct[branch_id]["tree"][res]["d2"]\

1819

["index"]]["momentum"]

1820

pid=decay_products[decay_struct[branch_id]["tree"][res]["d2"]\

1821

["index"]]["pid"]

1822

1823

indexd2=decay_struct[branch_id]["tree"][res]["d2"]["index"]

1824

if ( indexd2>0):

1825

istup=1

1826

external+=1

1827

else:

1828

istup=2

1829

decay_products[indexd2]["colup1"]=d2colup1

1830

decay_products[indexd2]["colup2"]=d2colup2

1831

1832

mothup1=part_number-2

1833

mothup2=part_number-2

1834

mass=mom.m

1835

helicity=0.

1836

decayed_event.particle[part_number]={"pid":pid,"istup":istup,\

1837

"mothup1":mothup1,"mothup2":mothup2,"colup1":d2colup1,\

1838

"colup2":d2colup2,\

1839

"momentum":mom,"mass":mass,"helicity":helicity}

1840

1841

decayed_event.event2mg[part_number]=part_number

1842

1843

else:

1844

external+=1

1845

part_number+=1

1846

decayed_event.particle[part_number]=curr_event.particle[part]

1847

decayed_event.event2mg[part_number]=part_number

1848

1849

else: # resonance in the production event

1850

if (ran==0): # write resonances in the prod. event ONLY if the

1851

# decayed event is ready to be written down

1852

part=curr_event.event2mg[index]

1853

part_number+=1

1854

decayed_event.particle[part_number]=curr_event.resonance[part]

1855

decayed_event.event2mg[part_number]=part_number

1856

# Here I need to check that the daughters still have the correct mothup1 and mothup2

1857

for part in curr_event.resonance.keys():

1858

mothup1=curr_event.resonance[part]["mothup1"]

1859

mothup2=curr_event.resonance[part]["mothup2"]

1860

if mothup1==index:

1861

if mothup2!=index: print "Warning: mothup1!=mothup2"

1862

curr_event.resonance[part]["mothup1"]=part_number

1863

curr_event.resonance[part]["mothup2"]=part_number

1864

for part in curr_event.particle.keys():

1865

mothup1=curr_event.particle[part]["mothup1"]

1866

mothup2=curr_event.particle[part]["mothup2"]

1867

if mothup1==index:

1868

if mothup2!=index: print "Warning: mothup1!=mothup2"

1869

curr_event.particle[part]["mothup1"]=part_number

1870

curr_event.particle[part]["mothup2"]=part_number

1871

1872

decayed_event.nexternal=part_number

1873

1874

return decayed_event, weight

1875

1876

1877

def get_topologies(self,matrix_element):

1878

"""Extraction of the phase-space topologies from mg5 matrix elements

1879

This is used for the production matrix element only.

1880

1881

the routine is essentially equivalent to write_configs_file_from_diagrams

1882

except that I don't write the topology in a file,

1883

I record it in an object production_topo (the class is defined above in this file)

1884

"""

1885

1886

# Extract number of external particles

1887

( nexternal, ninitial) = matrix_element.get_nexternal_ninitial()

1888

1889

del nexternal

1890

preconfigs = [(i+1, d) for i,d in enumerate(matrix_element.get('diagrams'))]

1891

mapconfigs = [c[0] for c in preconfigs]

1892

configs=[[c[1]] for c in preconfigs]

1893

model = matrix_element.get('processes')[0].get('model')

1894

1895

1896

topologies ={} # dictionnary {mapconfig number -> production_topology}

1897

# this is the object to be returned at the end of this routine

1898

1899

s_and_t_channels = []

1900

1901

minvert = min([max([d for d in config if d][0].get_vertex_leg_numbers()) \

1902

for config in configs])

1903

1904

# Number of subprocesses

1905

# nsubprocs = len(configs[0])

1906

1907

nconfigs = 0

1908

1909

new_pdg = model.get_first_non_pdg()

1910

1911

for iconfig, helas_diags in enumerate(configs):

1912

if any([vert > minvert for vert in

1913

[d for d in helas_diags if d][0].get_vertex_leg_numbers()]):

1914

# Only 3-vertices allowed in configs.inc

1915

continue

1916

nconfigs += 1

1917

1918

# Need s- and t-channels for all subprocesses, including

1919

# those that don't contribute to this config

1920

empty_verts = []

1921

stchannels = []

1922

for h in helas_diags:

1923

if h:

1924

# get_s_and_t_channels gives vertices starting from

1925

# final state external particles and working inwards

1926

stchannels.append(h.get('amplitudes')[0].\

1927

get_s_and_t_channels(ninitial, new_pdg))

1928

else:

1929

stchannels.append((empty_verts, None))

1930

1931

# For t-channels, just need the first non-empty one

1932

tchannels = [t for s,t in stchannels if t != None][0]

1933

1934

# For s_and_t_channels (to be used later) use only first config

1935

s_and_t_channels.append([[s for s,t in stchannels if t != None][0],

1936

tchannels])

1937

1938

1939

1940

# Make sure empty_verts is same length as real vertices

1941

if any([s for s,t in stchannels]):

1942

empty_verts[:] = [None]*max([len(s) for s,t in stchannels])

1943

1944

# Reorganize s-channel vertices to get a list of all

1945

# subprocesses for each vertex

1946

schannels = zip(*[s for s,t in stchannels])

1947

else:

1948

schannels = []

1949

1950

allchannels = schannels

1951

if len(tchannels) > 1:

1952

# Write out tchannels only if there are any non-trivial ones

1953

allchannels = schannels + tchannels

1954

1955

# Write out propagators for s-channel and t-channel vertices

1956

1957

# use the AMP2 index to label the topologies

1958

tag_topo=mapconfigs[iconfig]

1959

topologies[tag_topo]=production_topo()

1960

1961

for verts in allchannels:

1962

if verts in schannels:

1963

vert = [v for v in verts if v][0]

1964

else:

1965

vert = verts

1966

daughters = [leg.get('number') for leg in vert.get('legs')[:-1]]

1967

last_leg = vert.get('legs')[-1]

1968

1969

1970

if verts in schannels:

1971

type_propa="s"

1972

elif verts in tchannels[:-1]:

1973

type_propa="t"

1974

1975

1976

if (type_propa):

1977

topologies[tag_topo].add_one_branching(last_leg.get('number'),\

1978

daughters[0],daughters[1],type_propa)

1979

1980

return topologies

1981

1982

@misc.mute_logger()

1983

def generate_fortran_me(self,processes,base_model,mode, mgcmd,path_me):

1984

"""Given a process and a model, use the standanlone module of mg5

1985

to generate a fortran executable for the evaluation of the

1986

corresponding matrix element

1987

mode=0 : production part

1988

mode=1 : process fully decayed

1989

"""

1990

1991

commandline="import model "+base_model

1992

mgcmd.exec_cmd(commandline)

1993

1994

mgcmd.exec_cmd("set group_subprocesses False")

1995

1996

commandline="generate "+processes[0]

1997

mgcmd.exec_cmd(commandline)

1998

1999

# output the result in Fortran format:

2000

if mode==0: # production process

2001

mgcmd.exec_cmd("output standalone_ms %s -f" % pjoin(path_me,'production_me') )

2002

2003

elif mode==1: # full process

2004

mgcmd.exec_cmd("output standalone_ms %s -f" % pjoin(path_me,'full_me'))

2005

2006

2007

2008

2009

# now extract the information about the topologies

2010

if mode==0:

2011

me_list=mgcmd._curr_matrix_elements.get_matrix_elements()

2012

if(len(me_list)!=1):

2013

logger.warning('WARNING: unexpected number of matrix elements')

2014

topo=self.get_topologies(me_list[0])

2015

return topo

2016

2017

2018

def get_resonances(self,decay_processes):

2019

""" return a list of labels of each resonance involved in the decay chain """

2020

2021

resonances=[]

2022

for line in decay_processes:

2023

line=line.replace(">", " > ")

2024

line=line.replace("(", " ( ")

2025

line=line.replace(")", " ) ")

2026

list_proc=line.split()

2027

for i , item in enumerate(list_proc):

2028

if item ==">":

2029

resonances.append(list_proc[i-1])

2030

return resonances

2031

2032

2033

def get_full_process_structure(self,decay_processes,line_prod_proc, base_model, banner, proc_option, check=0):

2034

""" return a string with the definition of the process fully decayed

2035

and also a list of dc_branch objects with all infomation about the topology

2036

of each decay branch

2037

"""

2038

2039

decay_struct={}

2040

full_proc_line=line_prod_proc+" "+proc_option+" , "

2041

for proc_index in decay_processes.keys():

2042

if ',' in decay_processes[proc_index]:

2043

current_branch=decay_processes[proc_index]

2044

# list_branch=current_branch.split(",")

2045

# current_nb_decays=len(list_branch)

2046

# for nu in range(current_nb_decays):

2047

# if nu >0 and nu < current_nb_decays-1: list_branch[nu]=" ( "+list_branch[nu]

2048

# for nu in range(current_nb_decays-2):

2049

# list_branch[current_nb_decays-1]=list_branch[current_nb_decays-1]+" ) "

2050

# current_branch=""

2051

# for nu in range(current_nb_decays-1):

2052

# current_branch+=list_branch[nu]+ " , "

2053

# current_branch+=list_branch[current_nb_decays-1]

2054

full_proc_line=full_proc_line+" ( "+current_branch+ " ) , "

2055

else:

2056

full_proc_line=full_proc_line+" "+\

2057

decay_processes[proc_index]+ " , "

2058

decay_proc=decay_processes[proc_index]

2059

decay_proc=decay_proc.replace("\n","")

2060

decay_proc=decay_proc.replace("("," ")

2061

decay_proc=decay_proc.replace(")"," ")

2062

2063

decay_struct[proc_index]=dc_branch(decay_proc,base_model,banner,check)

2064

# decay_struct[proc_index].print_branch()

2065

# decay_struct[proc_index].print_branch()

2066

full_proc_line=full_proc_line[:-3]

2067

#print full_proc_line

2068

2069

return full_proc_line, decay_struct

2070

2071

def compile_fortran_me_production(self, path_me):

2072

""" Compile the fortran executables associated with the evalutation of the

2073

matrix elements (production process)

2074

Returns the path to the fortran executable

2075

"""

2076

2077

list_prod=os.listdir(pjoin(path_me,"production_me/SubProcesses"))

2078

counter=0

2079

logger.debug("""Finalizing production me's """)

2080

2081

2082

for direc in list_prod:

2083

if direc[0]=="P":

2084

counter+=1

2085

prod_name=direc[string.find(direc,"_")+1:]

2086

2087

old_path=pjoin(path_me,'production_me','SubProcesses',direc)

2088

new_path=pjoin(path_me,'production_me','SubProcesses',prod_name)

2089

if os.path.isdir(new_path): shutil.rmtree(new_path)

2090

os.rename(old_path, new_path)

2091

2092

file_madspin=pjoin(MG5DIR, 'MadSpin', 'src', 'driver_prod.f')

2093

shutil.copyfile(file_madspin, pjoin(new_path,"check_sa.f"))

2094

2095

file_madspin=pjoin(MG5DIR, 'MadSpin', 'src', 'makefile_ms')

2096

shutil.copyfile(file_madspin, pjoin(new_path,"makefile") )

2097

2098

file=pjoin(path_me, 'param_card.dat')

2099

shutil.copyfile(file,pjoin(path_me,"production_me","Cards","param_card.dat"))

2100

2101

# files to produce the parameters:

2102

file_madspin=pjoin(MG5DIR, 'MadSpin', 'src', 'initialize.f')

2103

shutil.copyfile(file_madspin,pjoin(new_path,"initialize.f"))

2104

2105

file_madspin=pjoin(MG5DIR, 'MadSpin', 'src', 'lha_read_ms.f')

2106

shutil.copyfile(file_madspin, pjoin(path_me,"production_me","Source","MODEL","lha_read.f" ))

2107

shutil.copyfile(pjoin(path_me,'production_me','Source','MODEL','input.inc'),pjoin(new_path,'input.inc'))

2108

2109

2110

# COMPILATION

2111

# in case there are new DHELAS routines, we need to recompile

2112

misc.compile(arg=['clean'], cwd=pjoin(path_me,"production_me","Source", "DHELAS"), mode='fortran')

2113

misc.compile( cwd=pjoin(path_me,"production_me","Source","DHELAS"), mode='fortran')

2114

2115

misc.compile(arg=['clean'], cwd=pjoin(path_me,"production_me","Source", "MODEL"), mode='fortran')

2116

misc.compile( cwd=pjoin(path_me,"production_me","Source","MODEL"), mode='fortran')

2117

2118

#os.chdir(new_path)

2119

misc.compile(arg=['clean'], cwd=new_path, mode='fortran')

2120

misc.compile(arg=['init'],cwd=new_path,mode='fortran')

2121

misc.call('./init', cwd=new_path)

2122

2123

shutil.copyfile(pjoin(new_path,'parameters.inc'),

2124

pjoin(new_path,os.path.pardir, 'parameters.inc'))

2125

os.chdir(path_me)

2126

2127

misc.compile(cwd=new_path, mode='fortran')

2128

2129

if(os.path.getsize(pjoin(path_me,'production_me','SubProcesses', 'parameters.inc'))<10):

2130

raise Exception, "Parameters of the model were not written correctly ! "

2131

return prod_name

2132

2133

2134

def compile_fortran_me_full(self,path_me):

2135

""" Compile the fortran executables associated with the evalutation of the

2136

matrix elements (full process)

2137

Returns the path to the fortran executable

2138

"""

2139

2140

2141

list_full=os.listdir(pjoin(path_me,"full_me","SubProcesses"))

2142

2143

logger.debug("""Finalizing decay chain me's """)

2144

for direc in list_full:

2145

if direc[0]=="P":

2146

2147

decay_name=direc[string.find(direc,"_")+1:]

2148

2149

old_path=pjoin(path_me,'full_me','SubProcesses',direc)

2150

new_path=pjoin(path_me, 'full_me','SubProcesses',decay_name)

2151

2152

2153

if os.path.isdir(new_path): shutil.rmtree(new_path)

2154

os.rename(old_path, new_path)

2155

2156

2157

file_madspin=pjoin(MG5DIR, 'MadSpin', 'src', 'driver_full.f')

2158

shutil.copyfile(file_madspin, pjoin(new_path,"check_sa.f") )

2159

2160

2161

file_madspin=pjoin(MG5DIR, 'MadSpin', 'src', 'makefile_ms')

2162

shutil.copyfile(file_madspin, pjoin(new_path,"makefile") )

2163

2164

shutil.copyfile(pjoin(path_me,'full_me','Source','MODEL','input.inc'),pjoin(new_path,'input.inc'))

2165

2166

# write all the parameters:

2167

file_madspin=pjoin(MG5DIR, 'MadSpin', 'src', 'initialize.f')

2168

shutil.copyfile(file_madspin,pjoin(new_path,"initialize.f"))

2169

2170

file_madspin=pjoin(MG5DIR, 'MadSpin', 'src', 'lha_read_ms.f')

2171

shutil.copyfile(file_madspin, pjoin(path_me,"full_me","Source","MODEL","lha_read.f" ))

2172

2173

file=pjoin(path_me, 'param_card.dat')

2174

shutil.copyfile(file,pjoin(path_me,"full_me","Cards","param_card.dat"))

2175

2176

# BEGIN COMPILATION

2177

# in case there are new DHELAS routines, we need to recompile

2178

misc.compile(arg=['clean'], cwd=pjoin(path_me,"full_me","Source","DHELAS"), mode='fortran')

2179

misc.compile( cwd=pjoin(path_me,"full_me","Source","DHELAS"), mode='fortran')

2180

2181

misc.compile(arg=['clean'], cwd=pjoin(path_me,"full_me","Source","MODEL"), mode='fortran')

2182

misc.compile( cwd=pjoin(path_me,"full_me","Source","MODEL"), mode='fortran')

2183

2184

os.chdir(new_path)

2185

misc.compile(arg=['clean'], cwd=new_path, mode='fortran')

2186

misc.compile(arg=['init'],cwd=new_path,mode='fortran')

2187

misc.call('./init')

2188

shutil.copyfile('parameters.inc', '../parameters.inc')

2189

os.chdir(path_me)

2190

2191

# now we can compile check

2192

misc.compile(arg=['check'], cwd=new_path, mode='fortran')

2193

# END COMPILATION

2194

2195

2196

if(os.path.getsize(pjoin(path_me,'full_me','SubProcesses', 'parameters.inc'))<10):

2197

raise Exception, "Parameters of the model were not written correctly ! "

2198

2199

#decay_pattern=direc[string.find(direc,"_")+1:]

2200

#decay_pattern=decay_pattern[string.find(decay_pattern,"_")+1:]

2201

#decay_pattern=decay_pattern[string.find(decay_pattern,"_")+1:]

2202

2203

os.chdir(path_me)

2204

return decay_name

2205

2206

def restore_light_parton_masses(self,topo,event):

2207

""" masses of light partons were set to zero for

2208

the evaluation of the matrix elements

2209

now we need to restore the initial masses before

2210

writting the decayed event in the lhe file

2211

"""

2212

2213

light_partons=[21,1,2,3]

2214

for part in topo["get_id"].keys():

2215

if abs(topo["get_id"][part]) in light_partons:

2216

topo["get_mass2"][part]=event.particle[part]["mass"]**2

2217

2218

# need to check if last branch is a t-branching. If it is,

2219

# we need to update the value of branch["m2"]

2220

# since this will be used in the reshuffling procedure

2221

if len(topo["branchings"])>0: # Exclude 2>1 topologies

2222

if topo["branchings"][-1]["type"]=="t":

2223

if topo["branchings"][-2]["type"]!="t":

2224

logger.warning('last branching is t-channel')

2225

logger.warning('but last-but-one branching is not t-channel')

2226

else:

2227

part=topo["branchings"][-1]["index_d2"]

2228

if part >0: # reset the mass only if "part" is an external particle

2229

topo["branchings"][-2]["m2"]=math.sqrt(topo["get_mass2"][part])

2230

2231

def reorder_branch(self,branch):

2232

""" branch is a string with the definition of a decay chain

2233

If branch contains " A > B C , B > ... "

2234

reorder into " A > C B , B > ... "

2235

2236

"""

2237

branch=branch.replace(',', ' , ')

2238

branch=branch.replace(')', ' ) ')

2239

branch=branch.replace('(', ' ( ')

2240

list_branch=branch.split(" ")

2241

for index in range(len(list_branch)-1,-1,-1):

2242

if list_branch[index]==' ' or list_branch[index]=='': del list_branch[index]

2243

#print list_branch

2244

for index, item in enumerate(list_branch):

2245

if item =="," and list_branch[index+1]!="(":

2246

if list_branch[index-2]==list_branch[index+1]:

2247

# swap the two particles before the comma:

2248

temp=list_branch[index-2]

2249

list_branch[index-2]=list_branch[index-1]

2250

list_branch[index-1]=temp

2251

if item =="," and list_branch[index+1]=="(":

2252

if list_branch[index-2]==list_branch[index+2]:

2253

# swap the two particles before the comma:

2254

temp=list_branch[index-2]

2255

list_branch[index-2]=list_branch[index-1]

2256

list_branch[index-1]=temp

2257

2258

2259

new_branch=""

2260

for item in list_branch:

2261

new_branch+=item+" "

2262

2263

return new_branch, list_branch[0]

2264

2265

2266

def get_symm_fac(self, decay_processes):

2267

""" get the symmetry factor associated with associated resonance """

2268

2269

dico_initpart={}

2270

2271

for index in decay_processes.keys():

2272

line=decay_processes[index]

2273

list_obj=line.split()

2274

res=list_obj[0]

2275

if not dico_initpart.has_key(res):

2276

dico_initpart[res]=[]

2277

if line not in dico_initpart[res]:

2278

dico_initpart[res].append(line)

2279

2280

2281

symm_fac=1

2282

for res in dico_initpart.keys():

2283

nb_decay_channels=len(dico_initpart[res])

2284

for index in range(1,nb_decay_channels+1):

2285

symm_fac=symm_fac*index

2286

2287

return symm_fac

2288

2289

def set_light_parton_massless(self,topo):

2290

""" masses of light partons are set to zero for

2291

the evaluation of the matrix elements

2292

"""

2293

2294

light_partons=[21,1,2,3]

2295

for part in topo["get_id"].keys():

2296

if abs(topo["get_id"][part]) in light_partons :

2297

topo["get_mass2"][part]=0.0

2298

2299

# need to check if last branch is a t-branching. If it is,

2300

# we need to update the value of branch["m2"]

2301

# since this will be used in the reshuffling procedure

2302

if len(topo["branchings"])>0: # Exclude 2>1 topologies

2303

if topo["branchings"][-1]["type"]=="t":

2304

if topo["branchings"][-2]["type"]!="t":

2305

logger.info('last branching is t-channel')

2306

logger.info('but last-but-one branching is not t-channel')

2307

else:

2308

part=topo["branchings"][-1]["index_d2"]

2309

if part >0: # reset the mass only if "part" is an external particle

2310

topo["branchings"][-2]["m2"]=math.sqrt(topo["get_mass2"][part])

2311

2312

2313

def transpole(self,pole,width):

2314

2315

""" routine for the generation of a p^2 according to

2316

a Breit Wigner distribution

2317

the generation window is

2318

[ M_pole^2 - 30*M_pole*Gamma , M_pole^2 + 30*M_pole*Gamma ]

2319

"""

2320

2321

zmin = math.atan(-30.0)/width

2322

zmax = math.atan(30.0)/width

2323

2324

z=zmin+(zmax-zmin)*random.random()

2325

y = pole+width*math.tan(width*z)

2326

2327

jac=(width/math.cos(width*z))**2*(zmax-zmin)

2328

return y, jac

2329

2330

2331

def generate_BW_masses (self,topo, to_decay, pid2name, pid2width, pid2mass,s):

2332

"""Generate the BW masses of the particles to be decayed in the production event """

2333

2334

weight=1.0

2335

for part in topo["get_id"].keys():

2336

pid=topo["get_id"][part]

2337

if pid2name[pid] in to_decay:

2338

mass=0.25

2339

width=pid2width[pid]*pid2mass[pid]/(2.0*pid2mass[pid])**2

2340

virtualmass2, jac=self.transpole(mass,width)

2341

virtualmass2=virtualmass2*(2.0*pid2mass[pid])**2

2342

weight=weight*jac

2343

#print "need to generate BW mass of "+str(part)

2344

##print "mass: "+str(pid2mass[pid])

2345

#print "width: "+str(pid2width[pid])

2346

#print "virtual mass: "+str(math.sqrt(virtualmass2))

2347

#print "jac: "+str(jac)

2348

old_mass=topo["get_mass2"][part]

2349

topo["get_mass2"][part]=virtualmass2

2350

# sanity check

2351

if pid2mass[pid]<1e-6:

2352

logger.debug('A decaying particle has a mass of less than 1e-6 GeV')

2353

# for debugg purposes:

2354

if abs((pid2mass[pid]-math.sqrt(topo["get_mass2"][part]))/pid2mass[pid])>1.0 :

2355

logger.debug('Mass after BW smearing affected by more than 100 % (1)')

2356

logger.debug('Pole mass: '+str(pid2mass[pid]))

2357

logger.debug('Virtual mass: '+str(math.sqrt(topo["get_mass2"][part])))

2358

2359

#print topo["get_mass2"]

2360

2361

# need to check if last branch is a t-branching. If it is,

2362

# we need to update the value of branch["m2"]

2363

2364

if len(topo["branchings"])>0: # Exclude 2>1 topologies

2365

if topo["branchings"][-1]["type"]=="t":

2366

if topo["branchings"][-2]["type"]!="t":

2367

logger.debug('last branching is t-channel')

2368

logger.debug('but last-but-one branching is not t-channel')

2369

else:

2370

part=topo["branchings"][-1]["index_d2"]

2371

if part >0: # reset the mass only if "part" refers to an external particle

2372

old_mass=topo["branchings"][-2]["m2"]

2373

topo["branchings"][-2]["m2"]=math.sqrt(topo["get_mass2"][part])

2374

#sanity check

2375

2376

if abs(old_mass-topo["branchings"][-2]["m2"])>1e-10:

2377

if abs((old_mass-topo["branchings"][-2]["m2"])/old_mass)>1.0 :

2378

logger.debug('Mass after BW smearing affected by more than 100 % (2)')

2379

logger.debug('Previous value: '+ str(old_mass))

2380

logger.debug('New mass: '+ str((topo["branchings"][-2]["m2"])))

2381

try:

2382

pid=topo["get_id"][part]

2383

logger.debug('pole mass: %s' % pid2mass[pid])

2384

except Exception:

2385

pass

2386

return weight

2387

2388

def modify_param_card(self,pid2widths,path_me):

2389

"""Modify the param_card w/r to what is read from the banner:

2390

if the value of a width is set to zero in the banner,

2391

it is automatically set to its default value in this code

2392

"""

2393

2394

param_card=open(pjoin(path_me,'param_card.dat'), 'r')

2395

new_param_card=""

2396

while 1:

2397

line=param_card.readline()

2398

if line =="": break

2399

list_line=line.split()

2400

if len(list_line)>2:

2401

if list_line[0]=="DECAY" and int(list_line[1]) in pid2widths.keys():

2402

list_line[2]=str(pid2widths[int(list_line[1])])

2403

line=""

2404

for item in list_line:

2405

line+=item+ " "

2406

line+="\n"

2407

new_param_card+=line

2408

2409

param_card.close()

2410

param_card=open(pjoin(path_me, 'param_card.dat'), 'w')

2411

param_card.write(new_param_card)

2412

param_card.close()

2413

2414

2415

def select_one_topo(self,prod_values):

2416

#

2417

# prod_values[0] is the value of |M_prod|^2

2418

# prod_values[1], prod_values[2], ... are the values of individual diagrams

2419

# (called AMP2 in mg)

2420

#

2421

2422

# first evaluate the sum of all single diagram values

2423

total=0.0

2424

cumul=[0.0]

2425

for i in range(1,len(prod_values)):

2426

cumul.append(cumul[i-1]+float(prod_values[i]))

2427

total+=float(prod_values[i])

2428

2429

for i in range(len(cumul)): cumul[i]=cumul[i]/total

2430

2431

#print "Cumulative AMP2 values: "

2432

#print cumul

2433

select_topo=random.random()

2434

2435

for i in range(1,len(cumul)):

2436

if select_topo>cumul[i-1] and select_topo<cumul[i]:

2437

good_topo=i

2438

break

2439

2440

#print "Selected topology"

2441

#print good_topo

2442

return good_topo, cumul

2443

2444

def find_resonances(self,proc, branches):

2445

""" restore the resonances in a production process

2446

the expected decay chains (variable branches)

2447

were extracted from the banner

2448

"""

2449

2450

pos1=proc.find(">")

2451

list_finalstate=proc[pos1+1:].split()

2452

2453

for res in branches.keys():

2454

resFS=[ item for item in branches[res]["finalstate"]]

2455

for index in range(len(list_finalstate)-1,-1,-1):

2456

if list_finalstate[index] in resFS:

2457

pos2 = resFS.index(list_finalstate[index])

2458

del resFS[pos2]

2459

del list_finalstate[index]

2460

if resFS!=[]:

2461

logger.warning('CANNOT RECOGNIZE THE EXPECTED DECAY \

2462

CHAIN STRUCTURE IN PRODUCTION EVENT')

2463

return proc

2464

else:

2465

list_finalstate.append(res)

2466

2467

initstate=proc[:pos1]

2468

finalstate=""

2469

for part in list_finalstate:

2470

finalstate+=" "+part

2471

for res in branches.keys():

2472

finalstate+=" , "+branches[res]["branch"]

2473

newproc=initstate+" > "+finalstate+" "

2474

return newproc

2475

2476

2477

2478

def get_final_state_compact(self,final_state_full):

2479

2480

# prod_resonances={}

2481

dc_pos=final_state_full.find(",")

2482

2483

if dc_pos>0:

2484

branch_list=final_state_full.split(",")

2485

del branch_list[0]

2486

list_obj=final_state_full.split()

2487

final_state_compact=""

2488

to_be_deleted=[]

2489

for index, obj in enumerate(list_obj):

2490

if obj==">":

2491

to_be_deleted.append(list_obj[index-1])

2492

2493

for obj in list_obj:

2494

if obj!=">" and obj!="," and obj not in to_be_deleted:

2495

final_state_compact+=obj+" "

2496

2497

branches={}

2498

for branch in branch_list:

2499

list_part=branch.split()

2500

branches[list_part[0]]={"finalstate":list_part[2:]}

2501

branches[list_part[0]]["branch"], dummy= self.reorder_branch(branch)

2502

else:

2503

final_state_compact=final_state_full

2504

branches={}

2505

2506

return final_state_compact, branches

2507

2508

def get_banner(self):

2509

pass

2510

2511

2512

def set_cumul_proba_for_tag_decay(self,multi_decay_processes,decay_tags):

2513

"""

2514

"""

2515

2516

sum_br=0.0

2517

for index, tag_decay in enumerate(decay_tags):

2518

sum_br+=multi_decay_processes[tag_decay]['br']

2519

2520

cumul=0.0

2521

for index, tag_decay in enumerate(decay_tags):

2522

cumul+=multi_decay_processes[tag_decay]['br']/sum_br

2523

multi_decay_processes[tag_decay]['cumul_br']=cumul

2524

return sum_br

2525

2526

2527

2528

2529

def generate_tag_decay(self,multi_decay_processes,decay_tags ):

2530

""" generate randomly the tag for the decay config. using a probability law

2531

based on branching fractions

2532

"""

2533

2534

r=random.random()

2535

2536

for index, tag_decay in enumerate(decay_tags):

2537

if r < multi_decay_processes[tag_decay]['cumul_br']:

2538

return tag_decay

2539

if r==1.0:

2540

return decay_tags[-1]

2541

2542

2543

2544

def update_tag_decays(self, decay_tags, branch_tags):

2545

""" update the set of tags to identify the final state content of a decay channel """

2546

2547

new_decay_tags=[]

2548

for item1 in branch_tags:

2549

if len(decay_tags)==0:

2550

new_decay_tags.append( (item1,) )

2551

else:

2552

for item2 in decay_tags:

2553

new_decay_tags.append(item2+(item1,))

2554

return new_decay_tags

2555

2556

def get_mean_sd(self,list_obj):

2557

""" return the mean value and the standard deviation of a list of reals """

2558

sum=0.0

2559

N=float(len(list_obj))

2560

for item in list_obj:

2561

sum+=item

2562

mean=sum/N

2563

sd=0.0

2564

for item in list_obj:

2565

sd+=(item-mean)**2

2566

sd=sd/(N-1.0)

2567

return mean, sd

2568

2569

2570

def check_decay_tag(self, tag,tag_list,check_weights):

2571

""" check is check_weights[tag] is proportional to

2572

check_weights[x] for one x in tag_list """

2573

2574

for x in tag_list:

2575

# build a list of [el1/el2]

2576

ratio=[ check_weights[x][index]/item for index, item in enumerate(check_weights[tag])]

2577

mean, sd = self.get_mean_sd(ratio)

2578

2579

if sd<mean*1e-6: return x

2580

return 0

2581

2582

def check_param_card(self, param_card):

2583

raise Exception

2584

list_line=param_card.split('\n')

2585

output=""

2586

loop_block=0

2587

for line in list_line:

2588

if line=="": continue

2589

current_line=line.split()

2590

2591

if loop_block==1 and line[0]!='#': continue

2592

if loop_block==1 and line[0]=='#':

2593

loop_block=0

2594

continue

2595

2596

if len(current_line)<2:

2597

output+=line+"\n"

2598

elif current_line[0]!='Block' or current_line[1]!='loop':

2599

output+=line+"\n"

2600

else:

2601

loop_block=1

2602

2603

return output

2604

2605

def get_dico_branchindex2label(self, decay_processes):

2606

""" return a dictionary {index branch : label of decaying particle}"""

2607

2608

dico={}

2609

for part in decay_processes:

2610

line_proc=decay_processes[part].split()

2611

dico[part]=line_proc[0]

2612

2613

return dico

2614

2615

def process_decay_syntax(self,decay_processes):

2616

""" add spaces to avoid any confusion in the decay chain syntax """

2617

2618

for part in decay_processes:

2619

decay_processes[part]=decay_processes[part].replace(',',' , ')

2620

decay_processes[part]=decay_processes[part].replace('>',' > ')

2621

decay_processes[part]=decay_processes[part].replace(')',' ) ')

2622

decay_processes[part]=decay_processes[part].replace('(',' ( ')

2623

2624

2625

class decay_all_events:

2626

2627

@misc.mute_logger()

2628

def __init__(self, ms_interface, inputfile, mybanner, decay_processes,\

2629

prod_branches, proc_option, options):

2630

2631

self.calculator = {}

2632

self.calculator_nbcall = {}

2633

self.options = options

2634

max_weight_arg = options['max_weight']

2635

BW_effects = options['BW_effect']

2636

path_me = os.path.realpath(options['curr_dir'])

2637

self.path_me = path_me

2638

2639

# need to unbuffer all I/O in fortran, otherwise

2640

# the values of matrix elements are not passed to the Python script

2641

os.environ['GFORTRAN_UNBUFFERED_ALL']='y'

2642

# os.system(" export GFORTRAN_UNBUFFERED_ALL ")

2643

2644

2645

mgcmd = ms_interface.mg5cmd

2646

base_model = ms_interface.model

2647

2648

curr_dir=os.getcwd()

2649

2650

# Remove old stuff from previous runs

2651

# so that the current run is not confused

2652

if os.path.isfile(pjoin(path_me,"param_card.dat")):

2653

os.remove(pjoin(path_me,"param_card.dat"))

2654

2655

if os.path.isdir(pjoin(path_me,"production_me")):

2656

shutil.rmtree(pjoin(path_me,"production_me"))

2657

2658

if os.path.isdir(pjoin(path_me,"full_me")):

2659

shutil.rmtree(pjoin(path_me,"full_me"))

2660

decay_tools=decay_misc()

2661

2662

2663

# process a bit the decay chain strings, so that the code is not confused by the syntax

2664

decay_tools.process_decay_syntax(decay_processes)

2665

2666

2667

# we will need a dictionary pid > label

2668

pid2label_dict=pid2label(base_model)

2669

2670

mybanner.check_pid(pid2label_dict)

2671

pid2label_dict.update(label2pid(base_model))

2672

2673

# we will also need a dictionary pid > color_rep

2674

pid2color_dict=pid2color(base_model)

2675

2676

2677

# now overwrite the param_card.dat in Cards:

2678

param_card=mybanner['slha']

2679

#param_card=decay_tools.check_param_card( param_card)

2680

2681

# now we can write the param_card.dat:

2682

# Note that the width of each resonance in the

2683

# decay chain should be >0 , we will check that later on

2684

param=open(pjoin(path_me,'param_card.dat'),"w")

2685

param.write(param_card)

2686

param.close()

2687

2688

# extract all resonances in the decay:

2689

resonances=decay_tools.get_resonances(decay_processes.values())

2690

logger.info('List of resonances:')

2691

logger.info(resonances)

2692

2693

label2width={}

2694

label2mass={}

2695

pid2width={}

2696

pid2mass={}

2697

# now extract the width of the resonances:

2698

for particle_label in resonances:

2699

try:

2700

part=pid2label_dict[particle_label]

2701

mass = mybanner.get('param_card','mass', abs(part))

2702

width = mybanner.get('param_card','decay', abs(part))

2703

except ValueError, error:

2704

continue

2705

else:

2706

label2width[particle_label]=float(width.value)

2707

label2mass[particle_label]=float(mass.value)

2708

pid2mass[part]=label2mass[particle_label]

2709

pid2width[part]=label2width[particle_label]

2710

if label2width[particle_label]==0.0:

2711

for param in base_model["parameters"][('external',)]:

2712

if param.lhablock=="DECAY" and param.lhacode==[abs(part)]:

2713

label2width[particle_label]=param.value

2714

pid2width[part]=label2width[particle_label]

2715

logger.warning('ATTENTION')

2716

logger.warning('Found a zero width in the param_card for particle '\

2717

+str(particle_label))

2718

logger.warning('Use instead the default value '\

2719

+str(label2width[particle_label]))

2720

2721

2722

# now we need to modify the values of the width

2723

# in param_card.dat, since this is where the input

2724

# parameters will be read when evaluating matrix elements

2725

decay_tools.modify_param_card(pid2width,path_me)

2726

2727

2728

# now we need to evaluate the branching fractions:

2729

# =================================================

2730

logger.info('Determining partial decay widths...')

2731

2732

calculate_br = width_estimate(resonances, path_me, pid2label_dict,

2733

mybanner, base_model)

2734

# Maybe the branching fractions are already given in the banner:

2735

calculate_br.extract_br_from_banner(mybanner)

2736

calculate_br.print_branching_fractions()

2737

#

2738

# now we check that we have all needed pieces of info regarding the branching fraction:

2739

multiparticles = mgcmd._multiparticles

2740

branching_per_channel=calculate_br.get_BR_for_each_decay(decay_processes,

2741

multiparticles)

2742

2743

2744

# check that we get branching fractions for all resonances to be decayed:

2745

2746

if branching_per_channel == 0:

2747

logger.debug('We need to recalculate the branching fractions')

2748

if hasattr(base_model.get('particles')[0], 'partial_widths'):

2749

logger.debug('using the compute_width module of madevent')

2750

calculate_br.launch_width_evaluation(resonances,base_model, mgcmd) # use FR to get all partial widths # set the br to partial_width/total_width

2751

else:

2752

logger.debug('compute_width module not available, use numerical estimates instead ')

2753

calculate_br.extract_br_from_width_evaluation()

2754

calculate_br.print_branching_fractions()

2755

branching_per_channel=calculate_br.get_BR_for_each_decay(decay_processes,multiparticles)

2756

2757

if branching_per_channel==0:

2758

raise Exception, 'Failed to extract the branching fraction associated with each decay channel'

2759

2760

2761

2762

2763

2764

# now we need to sort all the different decay configurations, and get the br for each of them

2765

# ===========================================================================================

2766

2767

# 1. first create a list of branches that ordered according to the canonical numeratation 1, 2, ...:

2768

list_particle_to_decay=decay_processes.keys()

2769

list_particle_to_decay.sort()

2770

2771

# 2. then use a tuple to identify a decay channel

2772

#

2773

# (tag1 , tag2 , tag3, ...)

2774

#

2775

# where the number of entries is the number of branches,

2776

# tag1 is the tag that identifies the final state of branch 1,

2777

# tag2 is the tag that identifies the final state of branch 2,

2778

# etc ...

2779

2780

# 2.a. first get decay_tags = the list of all the tuples

2781

2782

decay_tags=[]

2783

for part in list_particle_to_decay: # loop over particle to decay in the production process

2784

branch_tags=[ fs for fs in branching_per_channel[part]] # list of tags in a given branch

2785

decay_tags=decay_tools.update_tag_decays(decay_tags, branch_tags)

2786

2787

# 2.b. then build the dictionary multi_decay_processes = multi dico

2788

# first key = a tag in decay_tags

2789

# second key : ['br'] = float with the branching fraction for the FS associated with 'tag'

2790

# ['config'] = a list of strings, each of then giving the definition of a branch

2791

2792

multi_decay_processes={}

2793

for tag in decay_tags:

2794

# compute br + get the congis

2795

br=1.0

2796

list_branches={}

2797

for index, part in enumerate(list_particle_to_decay):

2798

br=br*branching_per_channel[part][tag[index]]['br']

2799

list_branches[part]=branching_per_channel[part][tag[index]]['config']

2800

multi_decay_processes[tag]={}

2801

multi_decay_processes[tag]['br']=br

2802

multi_decay_processes[tag]['config']=list_branches

2803

2804

# compute the cumulative probabilities associated with the branching fractions

2805

# keep track of sum of br over the decay channels at work, since we need to rescale

2806

# the event weight by this number

2807

sum_br=decay_tools.set_cumul_proba_for_tag_decay(multi_decay_processes,decay_tags)

2808

2809

2810

2811

# Now we have a dictionary (multi_decay_processes) of which values

2812

# identifies all the decay channels to be considered, and the associated branching fractions.

2813

2814

# The only difference with the one-channel implementation resides in the fact that

2815

# decay_processes is replaced by multi_decay_processes

2816

2817

2818

# A few initialisations:

2819

#========================

2820

# consider the possibility of several production process

2821

set_of_processes=[]

2822

decay_struct={}

2823

full_proc_line={}

2824

decay_path={} # dictionary to record the name of the directory with decay fortran me

2825

production_path={} # dictionary to record the name of the directory with production fortran me

2826

# also for production matrix elements,

2827

# we need to keep track of the topologies

2828

topologies={}

2829

2830

dico_branchindex2label=decay_tools.get_dico_branchindex2label(decay_processes)

2831

symm_fac=decay_tools.get_symm_fac(decay_processes)

2832

logger.info('Symmetry factor: '+str(symm_fac))

2833

2834

2835

#Next step: we need to determine which matrix elements are really necessary

2836

#==========================================================================

2837

2838

# create a dictionnary that maps a 'tag_decay' onto 'tag_me_decay' = tag for the matrix elements

2839

# call this map 'map_decay_me'

2840

logger.info('Checking the equality of matrix elements for different decay channels ... ')

2841

check_weights={}

2842

curr_event=Event(inputfile)

2843

curr_event.get_next_event()

2844

to_decay_map, to_decay_label =curr_event.get_map_resonances(dico_branchindex2label, pid2label_dict)

2845

2846

# print to_decay_map,

2847

# print to_decay_label

2848

2849

prod_process=curr_event.give_procdef(pid2label_dict)

2850

extended_prod_process=decay_tools.find_resonances(prod_process, prod_branches)

2851

2852

tag_production=prod_process.replace(">", "_")

2853

tag_production=tag_production.replace(" ","")

2854

tag_production=tag_production.replace("~","x")

2855

set_of_processes.append(tag_production)

2856

# generate fortran me for production only

2857

logger.debug(tag_production)

2858

topologies[tag_production]=\

2859

decay_tools.generate_fortran_me([extended_prod_process+proc_option],\

2860

mybanner.get("model"), 0,mgcmd,path_me)

2861

prod_name=decay_tools.compile_fortran_me_production(path_me)

2862

production_path[tag_production]=prod_name

2863

2864

# generate fortran me for the whole process

2865

decay_struct[tag_production]={}

2866

decay_path[tag_production]={}

2867

for tag_decay in decay_tags:

2868

check_weights[tag_decay]=[]

2869

new_full_proc_line, new_decay_struct=\

2870

decay_tools.get_full_process_structure(multi_decay_processes[tag_decay]['config'],\

2871

extended_prod_process, base_model, mybanner, proc_option)

2872

decay_struct[tag_production][tag_decay]=new_decay_struct

2873

#

2874

decay_tools.generate_fortran_me([new_full_proc_line+proc_option],\

2875

mybanner.get("model"), 1,mgcmd,path_me)

2876

2877

decay_name=decay_tools.compile_fortran_me_full(path_me)

2878

decay_path[tag_production][tag_decay]=decay_name

2879

2880

# select a topology for the reshuffling

2881

p, p_str=curr_event.give_momenta()

2882

prod_values = self.calculate_matrix_element('prod',

2883

production_path[tag_production], p_str)

2884

prod_values=prod_values.replace("\n", "")

2885

prod_values=prod_values.split()

2886

mg5_me_prod = float(prod_values[0])

2887

tag_topo, cumul_proba = decay_tools.select_one_topo(prod_values)

2888

2889

# extract the canonical phase-space variable based on that topology

2890

topologies[tag_production][tag_topo].dress_topo_from_event(curr_event,to_decay_label)

2891

topologies[tag_production][tag_topo].extract_angles()

2892

2893

# Breit-Wigner effects + reshuffling

2894

decay_tools.set_light_parton_massless(topologies[tag_production][tag_topo])

2895

for dec in range(100):

2896

try_reshuffle=0

2897

while 1:

2898

try_reshuffle+=1

2899

if try_reshuffle > 10:

2900

logger.debug('current %s' % try_reshuffle)

2901

BW_weight_prod=decay_tools.generate_BW_masses(\

2902

topologies[tag_production][tag_topo], \

2903

to_decay_label.values(),pid2label_dict, pid2width,pid2mass, \

2904

curr_event.shat)

2905

2906

succeed=topologies[tag_production][tag_topo].reshuffle_momenta()

2907

# sanlity check

2908

for part in topologies[tag_production][tag_topo]['get_momentum'].keys():

2909

if part in to_decay_map and \

2910

topologies[tag_production][tag_topo]['get_momentum'][part].m<1.0:

2911

logger.debug('Mass of a particle to decay is less than 1 GeV')

2912

logger.debug('in reshuffling loop')

2913

# end sanity check

2914

if succeed:

2915

if try_reshuffle > 10:

2916

logger.debug('pass at %s' % try_reshuffle)

2917

break

2918

if try_reshuffle % 10 == 0:

2919

logger.debug( 'tried %ix to reshuffle the momenta, failed'% try_reshuffle)

2920

logger.debug( ' So let us try with another topology')

2921

tag_topo, cumul_proba=decay_tools.select_one_topo(prod_values)

2922

topologies[tag_production][tag_topo].dress_topo_from_event(\

2923

curr_event,to_decay_label)

2924

topologies[tag_production][tag_topo].extract_angles()

2925

2926

# sometimes not possible to set the masses of the external partons to zero,

2927

# keep the original masses

2928

# decay_tools.set_light_parton_massless(topologies\

2929

# [tag_production][tag_topo])

2930

continue

2931

#try_reshuffle=0

2932

if try_reshuffle >100:

2933

misc.sprint('fail 100 times')

2934

break

2935

2936

2937

decayed_event, BW_weight_decay=decay_tools.decay_one_event(\

2938

curr_event,decay_struct[tag_production][decay_tags[0]], \

2939

pid2color_dict, to_decay_map, pid2width, \

2940

pid2mass, resonances,BW_effects)

2941

2942

if decayed_event==0:

2943

logger.warning('failed to decay event properly')

2944

continue

2945

for tag_decay in decay_tags:

2946

# set the momenta for the production event and the decayed event:

2947

p, p_str=curr_event.give_momenta()

2948

p_full, p_full_str=decayed_event.give_momenta()

2949

2950

# Evaluate matrix elements

2951

# start with production weight:

2952

prod_values =self.calculate_matrix_element('prod',

2953

production_path[tag_production], p_str)

2954

2955

prod_values=prod_values.replace("\n", "")

2956

prod_values=prod_values.split()

2957

mg5_me_prod = float(prod_values[0])

2958

# then decayed weight:

2959

full_values =self.calculate_matrix_element('full',

2960

decay_path[tag_production][tag_decay], p_full_str)

2961

2962

mg5_me_full = float(full_values)

2963

#mg5_me_full = float(external.communicate(input=p_full_str)[0])

2964

mg5_me_full=mg5_me_full*BW_weight_prod*BW_weight_decay

2965

os.chdir(curr_dir)

2966

if(not mg5_me_full>0 or not mg5_me_prod >0 ):

2967

logger.warning('WARNING: NEGATIVE MATRIX ELEMENT !!')

2968

weight=mg5_me_full/mg5_me_prod

2969

check_weights[tag_decay].append(weight)

2970

2971

# verify if some matrix elements are identical up to an overal factor

2972

map_decay_me={}

2973

decay_me_tags=[]

2974

for tag_decay in decay_tags:

2975

tag=decay_tools.check_decay_tag(tag_decay,map_decay_me.values(),check_weights)

2976

if tag==0:

2977

map_decay_me[tag_decay]=tag_decay

2978

decay_me_tags.append(tag_decay)

2979

else:

2980

map_decay_me[tag_decay]=tag

2981

self.terminate_fortran_executables(decay_path[tag_production][tag_decay])

2982

shutil.rmtree(pjoin(self.path_me,'full_me', 'SubProcesses',decay_path[tag_production][tag_decay]))

2983

2984

logger.info('Out of %d decay channels, %s matrix elements are independent ' % (len(decay_tags), len(decay_me_tags)))

2985

# for tag in decay_me_tags:

2986

# logger.info(multi_decay_processes[tag])

2987

del check_weights

2988

# Estimation of the maximum weight

2989

#=================================

2990

inputfile.seek(0)

2991

del curr_event

2992

os.system("date")

2993

# Now we are ready to start the evaluation of the maximum weight

2994

if max_weight_arg>0:

2995

max_weight={}

2996

for tag_decay in decay_me_tags: max_weight[tag_decay]=max_weight_arg

2997

else:

2998

numberev=20 # number of events

2999

numberps=10000 # number of phase pace points per event

3000

3001

logger.info(' ')

3002

logger.info(' Estimating the maximum weight ')

3003

logger.info(' ***************************** ')

3004

logger.info(' Probing the first '+str(numberev)+' events')

3005

logger.info(' at '+str(numberps)+' phase space points')

3006

logger.info(' ')

3007

3008

curr_event=Event(inputfile)

3009

probe_weight=[]

3010

3011

starttime = time.time()

3012

for ev in range(numberev):

3013

probe_weight.append({})

3014

for tag_decay in decay_me_tags:

3015

probe_weight[ev][tag_decay]=0.0

3016

curr_event.get_next_event()

3017

to_decay_map, to_decay_label =\

3018

curr_event.get_map_resonances(dico_branchindex2label, pid2label_dict)

3019

# check if we have a new production process, in which case

3020

# we need to generate the corresponding matrix elements

3021

prod_process=curr_event.give_procdef(pid2label_dict)

3022

extended_prod_process=decay_tools.find_resonances(prod_process, prod_branches)

3023

3024

tag_production=prod_process.replace(">", "_")

3025

tag_production=tag_production.replace(" ","")

3026

tag_production=tag_production.replace("~","x")

3027

if tag_production not in set_of_processes: # we need to generate new matrix elements

3028

logger.info('Found a new process: '+extended_prod_process+proc_option)

3029

# logger.info(prod_process)

3030

# logger.info('Re-interpreted as ')

3031

# logger.info(extended_prod_process+proc_option)

3032

# logger.info( tag_production)

3033

logger.debug( ' -> need to generate the corresponding fortran matrix element ... ')

3034

set_of_processes.append(tag_production)

3035

3036

# generate fortran me for production only

3037

topologies[tag_production]=\

3038

decay_tools.generate_fortran_me([extended_prod_process+proc_option],\

3039

mybanner.get("model"), 0,mgcmd,path_me)

3040

3041

prod_name=decay_tools.compile_fortran_me_production(path_me)

3042

production_path[tag_production]=prod_name

3043

3044

# for the decay, we need to keep track of all possibilities for the decay final state:

3045

decay_struct[tag_production]={}

3046

decay_path[tag_production]={}

3047

for tag_decay in decay_tags:

3048

new_full_proc_line, new_decay_struct=\

3049

decay_tools.get_full_process_structure(multi_decay_processes[tag_decay]['config'],\

3050

extended_prod_process, base_model, mybanner,proc_option)

3051

decay_struct[tag_production][tag_decay]=new_decay_struct

3052

if tag_decay in decay_me_tags:

3053

full_proc_line[tag_decay]=new_full_proc_line

3054

#

3055

for tag_decay in decay_me_tags:

3056

new_full_proc_line=full_proc_line[tag_decay]

3057

decay_tools.generate_fortran_me([new_full_proc_line+proc_option],\

3058

mybanner.get("model"), 1,mgcmd,path_me)

3059

3060

decay_name=decay_tools.compile_fortran_me_full(path_me)

3061

decay_path[tag_production][tag_decay]=decay_name

3062

3063

logger.debug('Done.')

3064

3065

# Now the relevant matrix elements for the current event are there.

3066

# But we still need to select a production topology

3067

# So first we evaluate the production me's only,

3068

# and select one production topolgy randomly,

3069

# with each topology weighted by the corresponding diagram squared.

3070

3071

3072

p, p_str=curr_event.give_momenta()

3073

# Note here that no momentum reshuffling is done,

3074

# since we don't know yet which topology should be used.

3075

# so light quarks and gluons in the final state may have a small mass

3076

3077

3078

prod_values = self.calculate_matrix_element('prod',

3079

production_path[tag_production], p_str)

3080

prod_values=prod_values.replace("\n", "")

3081

prod_values=prod_values.split()

3082

mg5_me_prod = float(prod_values[0])

3083

3084

tag_topo, cumul_proba = decay_tools.select_one_topo(prod_values)

3085

3086

# logger.debug(' ')

3087

running_time = misc.format_timer(time.time()-starttime)

3088

logger.debug('Event %s %s: ' % (ev+1, running_time))

3089

# print "Shat: "+str(curr_event.shat)

3090

logger.debug('Number of production topologies : '\

3091

+str(len(topologies[tag_production].keys())))

3092

# print "Cumulative probabilities : "+str(cumul_proba)

3093

logger.debug('Selected topology : '\

3094

+str(tag_topo))

3095

# topologies[tag_production][tag_topo].print_topo()

3096

3097

# print "Event before reshuffling:"

3098

# print curr_event.string_event_compact()

3099

3100

# We have our topology now. So we can

3101

# 1. dress the topology with the momenta

3102

# in the production event,

3103

# 2. set the masses the light partons to zero,

3104

# 3. generate the BW masses,

3105

# 4. reshuffle the momenta in the production event

3106

# 5. pass the info to curr_event

3107

3108

topologies[tag_production][tag_topo].dress_topo_from_event(curr_event,to_decay_label)

3109

topologies[tag_production][tag_topo].extract_angles()

3110

3111

if BW_effects:

3112

decay_tools.set_light_parton_massless(topologies[tag_production][tag_topo])

3113

3114

for dec in range(numberps):

3115

3116

# try to reshuffle the momenta

3117

if BW_effects:

3118

try_reshuffle=0

3119

while 1:

3120

try_reshuffle+=1

3121

BW_weight_prod=decay_tools.generate_BW_masses(\

3122

topologies[tag_production][tag_topo], \

3123

to_decay_label.values(),pid2label_dict, pid2width,pid2mass, \

3124

curr_event.shat)

3125

3126

succeed=topologies[tag_production][tag_topo].reshuffle_momenta()

3127

# sanlity check

3128

for part in topologies[tag_production][tag_topo]['get_momentum'].keys():

3129

if part in to_decay_map and \

3130

topologies[tag_production][tag_topo]['get_momentum'][part].m<1.0:

3131

logger.debug('Mass of a particle to decay is less than 1 GeV')

3132

logger.debug('in reshuffling loop')

3133

# end sanity check

3134

if succeed: break

3135

if try_reshuffle==10:

3136

logger.debug( 'tried 10x to reshuffle the momenta, failed')

3137

logger.debug( ' So let us try with another topology')

3138

tag_topo, cumul_proba=decay_tools.select_one_topo(prod_values)

3139

topologies[tag_production][tag_topo].dress_topo_from_event(\

3140

curr_event,to_decay_label)

3141

topologies[tag_production][tag_topo].extract_angles()

3142

# keep original masses in the event

3143

# decay_tools.set_light_parton_massless(topologies\

3144

# [tag_production][tag_topo])

3145

try_reshuffle=0

3146

else:

3147

BW_weight_prod=1.0

3148

3149

topologies[tag_production][tag_topo].topo2event(curr_event,to_decay_label)

3150

3151

# if dec==0:

3152

# print "Event after reshuffling:"

3153

# print curr_event.string_event_compact()

3154

3155

# Here the production event has been reshuffled,

3156

# now we need to decay it.

3157

# There might be several decay channels -> loop over them

3158

3159

for tag_decay in decay_me_tags:

3160

3161

decayed_event, BW_weight_decay=decay_tools.decay_one_event(\

3162

curr_event,decay_struct[tag_production][tag_decay], \

3163

pid2color_dict, to_decay_map, pid2width, \

3164

pid2mass, resonances,BW_effects)

3165

3166

if decayed_event==0:

3167

logger.warning('failed to decay event properly')

3168

continue

3169

# set the momenta for the production event and the decayed event:

3170

p, p_str=curr_event.give_momenta()

3171

p_full, p_full_str=decayed_event.give_momenta()

3172

3173

# start with production weight:

3174

prod_values =self.calculate_matrix_element('prod',

3175

production_path[tag_production], p_str)

3176

3177

prod_values=prod_values.replace("\n", "")

3178

prod_values=prod_values.split()

3179

mg5_me_prod = float(prod_values[0])

3180

# then decayed weight:

3181

full_values =self.calculate_matrix_element('full',

3182

decay_path[tag_production][tag_decay], p_full_str)

3183

3184

mg5_me_full = float(full_values)

3185

#mg5_me_full = float(external.communicate(input=p_full_str)[0])

3186

mg5_me_full=mg5_me_full*BW_weight_prod*BW_weight_decay

3187

os.chdir(curr_dir)

3188

3189

if(not mg5_me_full>0 or not mg5_me_prod >0 ):

3190

logger.warning('WARNING: NEGATIVE MATRIX ELEMENT !!')

3191

3192

weight=mg5_me_full/mg5_me_prod

3193

if (weight>probe_weight[ev][tag_decay]): probe_weight[ev][tag_decay]=weight

3194

3195

for index,tag_decay in enumerate(decay_me_tags):

3196

logger.info('Event '+str(ev+1)+\

3197

' , decay config '+str(index+1)+' : '+str(probe_weight[ev][tag_decay])+' %s' % running_time)

3198

3199

# Computation of the maximum weight used in the unweighting procedure

3200

max_weight={}

3201

3202

for index,tag_decay in enumerate(decay_me_tags):

3203

weights=[]

3204

for ev in range(numberev):

3205

weights.append(probe_weight[ev][tag_decay])

3206

ave_weight, std_weight=decay_tools.get_mean_sd(weights)

3207

std_weight=math.sqrt(std_weight)

3208

logger.info(' ')

3209

logger.info(' Decay channel '+str(index+1))

3210

for part in multi_decay_processes[tag_decay]['config'].keys():

3211

logger.info(multi_decay_processes[tag_decay]['config'][part])

3212

# logger.info(' average maximum weight that we got is '+str(ave_weight))

3213

# logger.info(' with a standard deviation of '+str(std_weight))

3214

# logger.info(' -> W_max = average + 4 * standard deviation')

3215

max_weight[tag_decay]=ave_weight+4.0*std_weight

3216

logger.info(' Using maximum weight '+str(max_weight[tag_decay]))

3217

del curr_event

3218

3219

3220

3221

#

3222

os.system("date")

3223

logger.info(' ' )

3224

logger.info('Decaying the events... ')

3225

inputfile.seek(0)

3226

outputfile = open(pjoin(path_me,'decayed_events.lhe'), 'w')

3227

curr_event=Event(inputfile)

3228

ms_banner = ""

3229

total_br = 0

3230

for index,tag_decay in enumerate(decay_tags):

3231

ms_banner+="# Decay channel "+str(index+1)+"\n"

3232

for part in multi_decay_processes[tag_decay]['config'].keys():

3233

ms_banner+="# "+multi_decay_processes[tag_decay]['config'][part]+"\n"

3234

ms_banner+="# branching fraction: "+str(multi_decay_processes[tag_decay]['br']) + "\n"

3235

ms_banner+="# estimate of the maximum weight: "+str(max_weight[map_decay_me[tag_decay]]) + "\n"

3236

total_br += multi_decay_processes[tag_decay]['br']

3237

self.branching_ratio = total_br

3238

mybanner['madspin'] += ms_banner

3239

# Update cross-section in the banner

3240

if 'mggenerationinfo' in mybanner:

3241

mg_info = mybanner['mggenerationinfo'].split('\n')

3242

for i,line in enumerate(mg_info):

3243

if 'Events' in line:

3244

continue

3245

if ':' not in line:

3246

continue

3247

info, value = line.rsplit(':',1)

3248

try:

3249

value = float(value)

3250

except:

3251

continue

3252

mg_info[i] = '%s : %s' % (info, value * total_br)

3253

mybanner['mggenerationinfo'] = '\n'.join(mg_info)

3254

if 'init' in mybanner:

3255

new_init =''

3256

for line in mybanner['init'].split('\n'):

3257

if len(line.split()) != 4:

3258

new_init += '%s\n' % line

3259

else:

3260

data = [float(nb) for nb in line.split()]

3261

data[:3] = [ data[i] *total_br for i in range(3)]

3262

new_init += ' %.12E %.12E %.12E %i\n' % tuple(data)

3263

mybanner['init'] = new_init

3264

mybanner.write(outputfile, close_tag=False)

3265

3266

event_nb=0

3267

trial_nb_all_events=0

3268

starttime = time.time()

3269

while 1:

3270

if (curr_event.get_next_event()=="no_event"): break

3271

event_nb+=1

3272

if (event_nb % max(int(10**int(math.log10(float(event_nb)))),10)==0):

3273

running_time = misc.format_timer(time.time()-starttime)

3274

logger.info('Event nb %s %s' % (event_nb, running_time))

3275

trial_nb=0

3276

to_decay_map, to_decay_label =\

3277

curr_event.get_map_resonances(dico_branchindex2label, pid2label_dict)

3278

3279

# if event_nb>10: break

3280

# check if we have a new production process, in which case we need to generate the correspomding matrix elements

3281

prod_process=curr_event.give_procdef(pid2label_dict)

3282

extended_prod_process=decay_tools.find_resonances(prod_process, prod_branches)

3283

3284

tag_production=prod_process.replace(">", "_")

3285

tag_production=tag_production.replace(" ","")

3286

tag_production=tag_production.replace("~","x")

3287

if tag_production not in set_of_processes: # we need to generate new matrix elements

3288

# logger.info(' ')

3289

logger.info('Found a new process: '+extended_prod_process+proc_option)

3290

# logger.info(prod_process)

3291

# logger.info('Re-interpreted as ')

3292

# logger.info(extended_prod_process+proc_option)

3293

logger.debug( tag_production)

3294

logger.debug( ' -> need to generate the corresponding fortran matrix element ... ')

3295

set_of_processes.append(tag_production)

3296

3297

# generate fortran me for production only

3298

topologies[tag_production]=\

3299

decay_tools.generate_fortran_me([extended_prod_process+proc_option],\

3300

mybanner.get("model"), 0,mgcmd,path_me)

3301

prod_name=decay_tools.compile_fortran_me_production(path_me)

3302

production_path[tag_production]=prod_name

3303

3304

# for the decay, we need to keep track of all possibilities for the decay final state:

3305

decay_struct[tag_production]={}

3306

decay_path[tag_production]={}

3307

for tag_decay in decay_tags:

3308

new_full_proc_line, new_decay_struct=\

3309

decay_tools.get_full_process_structure(multi_decay_processes[tag_decay]['config'],\

3310

extended_prod_process, base_model, mybanner,proc_option)

3311

decay_struct[tag_production][tag_decay]=new_decay_struct

3312

if tag_decay in decay_me_tags:

3313

full_proc_line[tag_decay]=new_full_proc_line

3314

#

3315

for tag_decay in decay_me_tags:

3316

new_full_proc_line=full_proc_line[tag_decay]

3317

decay_tools.generate_fortran_me([new_full_proc_line+proc_option],\

3318

mybanner.get("model"), 1,mgcmd,path_me)

3319

3320

decay_name=decay_tools.compile_fortran_me_full(path_me)

3321

decay_path[tag_production][tag_decay]=decay_name

3322

3323

logger.debug('Done.')

3324

3325

3326

# First evaluate production matrix element

3327

p, p_str=curr_event.give_momenta()

3328

prod_values = self.calculate_matrix_element('prod',

3329

production_path[tag_production], p_str)

3330

prod_values=prod_values.replace("\n", "")

3331

prod_values=prod_values.split()

3332

mg5_me_prod = float(prod_values[0])

3333

3334

# select topology based on sigle-diagram weights

3335

tag_topo, cumul_proba=decay_tools.select_one_topo(prod_values)

3336

3337

# dress the topology with momenta and extract the canonical numbers

3338

topologies[tag_production][tag_topo].dress_topo_from_event(curr_event,to_decay_label)

3339

topologies[tag_production][tag_topo].extract_angles()

3340

3341

if BW_effects:

3342

decay_tools.set_light_parton_massless(topologies[tag_production][tag_topo])

3343

3344

while 1:

3345

trial_nb+=1

3346

3347

# try to reshuffle the event:

3348

if BW_effects:

3349

try_reshuffle=0

3350

while 1:

3351

try_reshuffle+=1

3352

BW_weight_prod=decay_tools.generate_BW_masses(topologies[tag_production][tag_topo], \

3353

to_decay_label.values(),pid2label_dict, pid2width,pid2mass, \

3354

curr_event.shat)

3355

succeed=topologies[tag_production][tag_topo].reshuffle_momenta()

3356

if succeed: break

3357

if try_reshuffle==10:

3358

logger.debug('WARNING: tried 10x to reshuffle the momenta, failed')

3359

logger.debug(' So let us try with another topology')

3360

# print "Event: "+str(event_nb)

3361

# topologies[tag_production][tag_topo].print_topo()

3362

tag_topo, cumul_proba=decay_tools.select_one_topo(prod_values)

3363

topologies[tag_production][tag_topo].dress_topo_from_event(curr_event,to_decay_label)

3364

topologies[tag_production][tag_topo].extract_angles()

3365

3366

# keep original masses in the production event

3367

# decay_tools.set_light_parton_massless(topologies[tag_production][tag_topo])

3368

try_reshuffle=0

3369

3370

else:

3371

BW_weight_prod=1.0

3372

3373

# Here we need to select a decay configuration on a random basis:

3374

tag_decay=decay_tools.generate_tag_decay(multi_decay_processes,decay_tags)

3375

3376

topologies[tag_production][tag_topo].topo2event(curr_event,to_decay_label)

3377

decayed_event, BW_weight_decay=decay_tools.decay_one_event(curr_event,decay_struct[tag_production][tag_decay], \

3378

pid2color_dict, to_decay_map, pid2width, pid2mass, resonances,BW_effects)

3379

3380

if decayed_event==0:

3381

logger.info('failed to decay one event properly')

3382

continue # means we had mA<mB+mC in one splitting A->B+C

3383

p, p_str=curr_event.give_momenta()

3384

p_full, p_full_str=decayed_event.give_momenta()

3385

3386

# Now evaluate the matrix elements ...

3387

# start with production weight:

3388

prod_values = self.calculate_matrix_element('prod',

3389

production_path[tag_production], p_str)

3390

prod_values=prod_values.replace("\n", "")

3391

prod_values=prod_values.split()

3392

mg5_me_prod = float(prod_values[0])

3393

# then decayed weight:

3394

full_value = self.calculate_matrix_element('full',

3395

decay_path[tag_production][map_decay_me[tag_decay]], p_full_str)

3396

mg5_me_full = float(full_value)

3397

temp=mg5_me_full

3398

mg5_me_full=mg5_me_full*BW_weight_prod*BW_weight_decay

3399

3400

if(not mg5_me_full>0 or not mg5_me_prod >0 ):

3401

logger.warning('NEGATIVE MATRIX ELEMENT !!')

3402

3403

# mg5_me_prod, amp2_prod = evaluator.evaluate_matrix_element(me_prod[tag_production],p)

3404

# mg5_me_full, amp_full = evaluator.evaluate_matrix_element(me_full[tag_production],p_full)

3405

3406

weight=mg5_me_full/mg5_me_prod

3407

if weight>max_weight[map_decay_me[tag_decay]]:

3408

logger.info('warning: got a larger weight than max_weight estimate')

3409

logger.info('the ratio with the max_weight estimate is '+str(weight/max_weight[map_decay_me[tag_decay]]))

3410

logger.info('decay channel ')

3411

logger.info(multi_decay_processes[tag_decay])

3412

3413

if (weight/max_weight[map_decay_me[tag_decay]]> random.random()):

3414

3415

# Here we need to restore the masses of the light partons

3416

# initially found in the lhe production event

3417

decay_tools.restore_light_parton_masses(topologies[tag_production][tag_topo],curr_event)

3418

succeed=topologies[tag_production][tag_topo].reshuffle_momenta()

3419

if not succeed:

3420

logger.info('Warning: unable to restore masses of light partons')

3421

else:

3422

topologies[tag_production][tag_topo].topo2event(curr_event,to_decay_label)

3423

curr_event.reset_resonances() # re-evaluate the momentum of each resonance in prod. event

3424

decayed_event, BW_weight_decay=decay_tools.decay_one_event(curr_event,decay_struct[tag_production][tag_decay], \

3425

pid2color_dict, to_decay_map, pid2width, pid2mass, resonances,BW_effects,ran=0)

3426

decayed_event.wgt=decayed_event.wgt*sum_br*float(symm_fac)

3427

outputfile.write(decayed_event.string_event())

3428

# print "number of trials: "+str(trial_nb)

3429

trial_nb_all_events+=trial_nb

3430

break

3431

3432

os.system("date")

3433

outputfile.write('</LesHouchesEvents>\n')

3434

inputfile.close()

3435

outputfile.close()

3436

3437

logger.info('Total number of events: '+str(event_nb))

3438

logger.info('Average number of trial points per production event: '\

3439

+str(float(trial_nb_all_events)/float(event_nb)))

3440

logger.info('Number of subprocesses '+str(len(decay_path)))

3441

self.terminate_fortran_executables()

3442

shutil.rmtree(pjoin(path_me,'production_me'))

3443

shutil.rmtree(pjoin(path_me,'full_me'))

3444

3445

# set the environment variable GFORTRAN_UNBUFFERED_ALL

3446

# to its original value

3447

os.environ['GFORTRAN_UNBUFFERED_ALL']='n'

3448

3449

3450

3451

def terminate_fortran_executables(self, path_to_decay=0 ):

3452

"""routine to terminate all fortran executables"""

3453

3454

if not path_to_decay:

3455

for (mode, production) in self.calculator:

3456

external = self.calculator[(mode, production)]

3457

external.terminate()

3458

else:

3459

try:

3460

external = self.calculator[('full', path_to_decay)]

3461

except Exception:

3462

pass

3463

else:

3464

external.terminate()

3465

3466

def calculate_matrix_element(self, mode, production, stdin_text):

3467

"""routine to return the matrix element"""

3468

tmpdir = ''

3469

if (mode, production) in self.calculator:

3470

external = self.calculator[(mode, production)]

3471

self.calculator_nbcall[(mode, production)] += 1

3472

else:

3473

logger.debug('we have %s calculator ready' % len(self.calculator))

3474

if mode == 'prod':

3475

tmpdir = pjoin(self.path_me,'production_me', 'SubProcesses',

3476

production)

3477

else:

3478

tmpdir = pjoin(self.path_me,'full_me', 'SubProcesses',

3479

production)

3480

executable_prod="./check"

3481

external = Popen(executable_prod, stdout=PIPE, stdin=PIPE,

3482

stderr=STDOUT, cwd=tmpdir)

3483

self.calculator[(mode, production)] = external

3484

self.calculator_nbcall[(mode, production)] = 1

3485

3486

3487

external.stdin.write(stdin_text)

3488

if mode == 'prod':

3489

info = int(external.stdout.readline())

3490

nb_output = abs(info)+1

3491

else:

3492

info = 1

3493

nb_output = 1

3494

3495

3496

3497

prod_values = ' '.join([external.stdout.readline() for i in range(nb_output)])

3498

if info < 0:

3499

print 'ZERO DETECTED'

3500

print prod_values

3501

print stdin_text

3502

os.system('lsof -p %s' % external.pid)

3503

return ' '.join(prod_values.split()[-1*(nb_output-1):])

3504

3505

if len(self.calculator) > 100:

3506

logger.debug('more than 100 calculator. Perform cleaning')

3507

nb_calls = self.calculator_nbcall.values()

3508

nb_calls.sort()

3509

cut = max([nb_calls[len(nb_calls)//2], 0.001 * nb_calls[-1]])

3510

for key, external in list(self.calculator.items()):

3511

nb = self.calculator_nbcall[key]

3512

if nb < cut:

3513

external.stdin.close()

3514

external.stdout.close()

3515

external.terminate()

3516

del self.calculator[key]

3517

del self.calculator_nbcall[key]

3518

else:

3519

self.calculator_nbcall[key] = self.calculator_nbcall[key] //10

3520

3521

3522

return prod_values

3523

3524

3525