~chirality-flow/chiralityflow/ChiralityFlowMG

logical is_LC ! for not leading color bypass the writing of intermediate particle since the diagram is a very good candididate (and that it leads to issue)

double precision pt

integer get_color,elim_indices,set_colmp,fix_tchannel_color,combid

real ran1

external pt,ran1,get_color,elim_indices,set_colmp,fix_tchannel_color

if (first_time) then

include 'props.inc'

first_time=.false.

endif

npart = nexternal

buff = ' '

100

101

do i=-nexternal+2,nexternal

102

icolalt(1,i)=0

103

icolalt(2,i)=0

104

enddo

105

106

107

c Choose the config (diagram) which was actually used to produce the event

108

109

c ...unless the diagram is passed in igraphs(1); then use that diagram

110

lconfig=iconfig

111

if (ickkw.gt.0) then

112

if (btest(mlevel,3)) then

113

write(*,*)'unwgt.f: write out diagram ',igraphs(1)

114

endif

115

lconfig=igraphs(1)

116

endif

117

118

119

c Choose a color flow which is certain to work with the propagator

120

c structure of the chosen diagram and use that as an alternative

121

122

123

nc = int(jamp2(0))

124

is_LC = .true.

125

maxcolor=0

126

if(nc.gt.0)then

127

if(icolamp(1,%(iconfig)s,iproc)) then

128

targetamp(1)=jamp2(1)

129

c print *,'Color flow 1 allowed for config ',lconfig

130

else

131

targetamp(1)=0d0

132

endif

133

do ic =2,nc

134

if(icolamp(ic,%(iconfig)s,iproc))then

135

targetamp(ic) = jamp2(ic)+targetamp(ic-1)

136

c print *,'Color flow ',ic,' allowed for config ',lconfig,targetamp(ic)

137

else

138

targetamp(ic)=targetamp(ic-1)

139

endif

140

enddo

141

c ensure that at least one leading color is different of zero if not allow

142

c all subleading color.

143

if (targetamp(nc).eq.0)then

144

is_LC = .false.

145

targetamp(1)=jamp2(1)

146

do ic =2,nc

147

targetamp(ic) = jamp2(ic)+targetamp(ic-1)

148

enddo

149

endif

150

151

152

xtarget=ran1(iseed)*targetamp(nc)

153

154

ic = 1

155

do while (targetamp(ic) .lt. xtarget .and. ic .lt. nc)

156

ic=ic+1

157

enddo

158

if(targetamp(nc).eq.0) ic=0

159

c print *,'Chose color flow ',ic

160

do i=1,nexternal

161

if(ic.gt.0) then

162

icolalt(1,isym(i,jsym))=icolup(1,i,ic,numproc)

163

icolalt(2,isym(i,jsym))=icolup(2,i,ic,numproc)

164

c write(*,*) i, icolalt(1,isym(i,jsym)), icolalt(2,isym(i,jsym))

165

if (abs(icolup(1,i,ic, numproc)).gt.maxcolor) maxcolor=icolup(1,i,ic, numproc)

166

if (abs(icolup(2,i,ic, numproc)).gt.maxcolor) maxcolor=icolup(2,i,ic, numproc)

167

endif

168

enddo

169

else ! nc.gt.0

170

171

do i=1,nexternal

172

icolalt(1,i)=0

173

icolalt(2,i)=0

174

enddo

175

176

endif ! nc.gt.0

177

178

c Store original maxcolor to know if we have epsilon vertices

179

maxorg=maxcolor

180

c Keep track of IS colors that go through to final state

181

c (since we shouldn't replace pop-up indices with those)

182

do i=1,nincoming

183

do j=1,2

184

is_colors(j,i)=0

185

do k=3,nexternal

186

if (iabs(icolalt(j,i)).eq.iabs(icolalt(j,k))) then

187

c This color is going through to FS

188

is_colors(j,i)=iabs(icolalt(j,i))

189

endif

190

enddo

191

enddo

192

enddo

193

194

c Get mother information from chosen graph

195

196

197

c Set idij for external particles (needed to keep track of BWs)

198

if(ickkw.gt.0) then

199

do i=1,nexternal

200

idij(i)=ishft(1,i-1)

201

enddo

202

endif

203

204

c First check number of resonant s-channel propagators

205

ns=0

206

nres=0

207

tchannel=.false.

208

c Ensure that mother-daughter information starts from 0

209

do i=-nexternal+3,0

210

jpart(2,i) = 0

211

jpart(3,i) = 0

212

enddo

213

214

215

c Loop over propagators to find mother-daughter information

216

do i=-1,-nexternal+2,-1

217

c Daughters

218

if(i.gt.-nexternal+2)then

219

ida(1)=iforest(1,i,lconfig)

220

ida(2)=iforest(2,i,lconfig)

221

do j=1,2

222

if(ida(j).gt.0) ida(j)=isym(ida(j),jsym)

223

enddo

224

c Set idij (needed to keep track of BWs)

225

if(ickkw.gt.0) idij(i)=combid(idij(ida(1)),idij(ida(2)))

226

endif

227

c Decide s- or t-channel (for not LC -> set to none

228

if(i.gt.-nexternal+2.and.is_LC.and.

229

$ iabs(sprop(numproc,i,lconfig)).gt.0) then ! s-channel propagator

230

jpart(1,i)=sprop(numproc,i,lconfig)

231

ns=ns+1

232

else if(nres.gt.0.and.maxcolor.gt.maxorg.and.is_LC) then

233

c For t-channel propagators, just check that the colors are ok

234

if(i.eq.-nexternal+2) then

235

c This is the final t-channel, combining with leg 2

236

mo_color=0

237

if(.not.tchannel)then

238

c There are no previous t-channels, so this is a combination of

239

c 2, 1 and the last s-channel

240

ida(1)=1

241

ida(2)=i+1

242

else

243

c The daughter info is in iforest

244

ida(1)=iforest(1,i,lconfig)

245

ida(2)=iforest(2,i,lconfig)

246

endif

247

c Reverse colors of t-channels to get right color ordering

248

ncolmp=0

249

ncolmp=set_colmp(ncolmp,icolmp,2,jpart,

250

$ iforest(1,-max_branch,lconfig),icolalt,

251

$ icolalt(2,2),icolalt(1,2))

252

else

253

jpart(1,i)=tprid(i,lconfig)

254

mo_color=get_color(jpart(1,i))

255

ncolmp=0

256

endif

257

if(mo_color.gt.1.and.

258

$ mo_color.ne.3.and.mo_color.ne.8)then

259

da_color(1)=get_color(jpart(1,ida(1)))

260

da_color(2)=get_color(jpart(1,ida(2)))

261

call write_error(da_color(1), da_color(2), mo_color)

262

endif

263

c Set icolmp for daughters

264

ncolmp=set_colmp(ncolmp,icolmp,ida(2),jpart,

265

$ iforest(1,-max_branch,lconfig),icolalt,

266

$ icolalt(1,ida(2)),icolalt(2,ida(2)))

267

c Reverse colors of t-channels to get right color ordering

268

ncolmp=set_colmp(ncolmp,icolmp,ida(1),jpart,

269

$ iforest(1,-max_branch,lconfig),icolalt,

270

$ icolalt(2,ida(1)),icolalt(1,ida(1)))

271

c Fix t-channel color

272

c print *,'t-channel: ',i,ida(1),ida(2),mo_color

273

c print *,'colors: ',((icolmp(j,k),j=1,2),k=1,ncolmp)

274

maxcolor=fix_tchannel_color(mo_color,maxcolor,

275

$ ncolmp,icolmp,i,icolalt,is_colors)

276

tchannel=.true.

277

cycle

278

else

279

goto 100

280

endif

281

c Set status codes for propagator

282

c if((igscl(0).ne.0.and.

283

c $ (iabs(jpart(1,i)).gt.5.and.iabs(jpart(1,i)).lt.11).or.

284

c $ (iabs(jpart(1,i)).gt.16.and.iabs(jpart(1,i)).ne.21)).or.

285

c $ (igscl(0).eq.0.and.OnBW(i))) then

286

if(ickkw.eq.0.and.OnBW(i))then

287

c Resonance whose mass should be preserved

288

jpart(6,i)=2

289

nres=nres+1

290

else if (ickkw.gt.0) then

291

if(isbw(idij(i))) then

292

c Resonance whose mass should be preserved

293

jpart(6,i)=2

294

nres=nres+1

295

else

296

jpart(6,i)=3

297

endif

298

else

299

c Propagator for documentation only - not included

300

jpart(6,i)=3

301

endif

302

c Calculate momentum (p1+p2 for s-channel, p2-p1 for t-channel)

303

do j=0,3

304

pb(j,i)=pb(j,ida(1))+pb(j,ida(2))

305

enddo

306

pb(4,i)=sqrt(max(0d0,pb(0,i)**2-pb(1,i)**2-pb(2,i)**2-pb(3,i)**2))

307

c if(jpart(6,i).eq.2.and.

308

c $ abs(pb(4,i)-prmass(i,lconfig)).gt.5d0*prwidth(i,lconfig)) then

309

c jpart(6,i)=3

310

c nres=nres-1

311

c endif

312

c Set color info for all s-channels

313

mo_color = get_color(jpart(1,i))

314

c If inside multipart. vertex (indicated by color 2) cycle

315

c Set tentative mothers

316

jpart(2,i) = 1

317

jpart(3,i) = 2

318

c Set mother info for daughters

319

do j=1,2

320

jpart(2,ida(j)) = i

321

jpart(3,ida(j)) = i

322

enddo

323

if(mo_color.eq.2) cycle

324

c Reset list of color indices

325

ncolmp=0

326

c Add new color indices to list of color indices

327

do j=1,2

328

ncolmp=set_colmp(ncolmp,icolmp,ida(j),jpart,

329

$ iforest(1,-max_branch,lconfig),icolalt,

330

$ icolalt(1,ida(j)),icolalt(2,ida(j)))

331

enddo

332

c print *,'s-channel: ',i,mo_color,ida(1),ida(2)

333

c print *,'colors: ',((icolmp(j,k),j=1,2),k=1,ncolmp)

334

if(is_LC)then

335

if (icolmp(1,1).eq.1000.or.icolmp(2,1).eq.1000)then

336

if (jpart(6,i).eq.2)then

337

stop 5

338

endif

339

elseif(mo_color.eq.1) then ! color singlet

340

maxcolor=elim_indices(0,0,ncolmp,icolmp,i,icolalt,

341

$ is_colors,maxcolor)

342

elseif(mo_color.eq.-3) then ! color anti-triplet

343

maxcolor=elim_indices(0,1,ncolmp,icolmp,i,icolalt,

344

$ is_colors,maxcolor)

345

elseif(mo_color.eq.3) then ! color triplet

346

maxcolor=elim_indices(1,0,ncolmp,icolmp,i,icolalt,

347

$ is_colors,maxcolor)

348

elseif(mo_color.eq.-6) then ! color anti-sextet

349

maxcolor=elim_indices(0,2,ncolmp,icolmp,i,icolalt,

350

$ is_colors,maxcolor)

351

elseif(mo_color.eq.6) then ! color sextet

352

maxcolor=elim_indices(2,0,ncolmp,icolmp,i,icolalt,

353

$ is_colors,maxcolor)

354

elseif(mo_color.eq.8) then ! color octet

355

maxcolor=elim_indices(1,1,ncolmp,icolmp,i,icolalt,

356

$ is_colors,maxcolor)

357

else ! 2 indicates multipart. vertex

358

da_color(1) = get_color(jpart(1,ida(1)))

359

da_color(2) = get_color(jpart(1,ida(2)))

360

call write_error(da_color(1), da_color(2), mo_color)

361

endif

362

endif !end of check on LC

363

364

c Just zero helicity info for intermediate states

365

jpart(7,i) = 0

366

enddo ! do i

367

100 continue

368

if (is_LC) call check_pure_internal_flow(icolalt,jpart, maxcolor)

369

370

c Remove non-resonant mothers, set position of particles

371

ires=0

372

do i=-ns,nexternal

373

jpart(4,i)=icolalt(1,i)

374

jpart(5,i)=icolalt(2,i)

375

if(i.eq.1.or.i.eq.2) then

376

ito(i)=i ! initial state particle

377

else if(i.ge.3) then

378

ito(i)=i+nres ! final state particle

379

else if(i.le.-1.and.jpart(6,i).eq.2) then

380

ires=ires+1

381

ito(i)=2+ires ! s-channel resonances

382

else

383

ito(i)=0

384

if(i.eq.0) cycle

385

endif

386

if(jpart(2,i).lt.0.and.jpart(6,jpart(2,i)).ne.2) then

387

jpart(2,i)=jpart(2,jpart(2,i))

388

jpart(3,i)=jpart(3,jpart(3,i))

389

endif

390

enddo

391

392

c Shift particles to right place and set mothers of particles

393

394

do i=nexternal,-ns,-1

395

if(ito(i).le.0) cycle

396

do j=1,7

397

jpart(j,ito(i))=jpart(j,i)

398

enddo

399

if(jpart(2,ito(i)).lt.0) then

400

jpart(2,ito(i))=ito(jpart(2,ito(i)))

401

jpart(3,ito(i))=ito(jpart(3,ito(i)))

402

endif

403

do j=0,4

404

pb(j,ito(i))=pb(j,i)

405

enddo

406

enddo

407

408

c Set correct mother number for clustering info

409

410

if (icluster(1,1).ne.0) then

411

do i=1,nexternal-2

412

if(icluster(4,i).gt.0)then

413

icluster(4,i)=ito(icluster(4,i))

414

else

415

icluster(4,i)=-1

416

endif

417

if(icluster(3,i).eq.0)then

418

icluster(3,i)=-1

419

endif

420

if(ito(icluster(1,i)).gt.0)

421

$ icluster(1,i)=ito(icluster(1,i))

422

if(ito(icluster(2,i)).gt.0)

423

$ icluster(2,i)=ito(icluster(2,i))

424

if(flip)then

425

if(icluster(1,i).le.2)

426

$ icluster(1,i)=3-icluster(1,i)

427

if(icluster(2,i).le.2)

428

$ icluster(2,i)=3-icluster(2,i)

429

if(icluster(3,i).ge.1.and.icluster(3,i).le.2)

430

$ icluster(3,i)=3-icluster(3,i)

431

endif

432

enddo

433

endif

434

435

if (flip) then

436

c Need to flip initial state color, since might be overwritten

437

do i=1,7

438

j=jpart(i,1)

439

jpart(i,1)=jpart(i,2)

440

jpart(i,2)=j

441

enddo

442

endif

443

444

if(ickkw.gt.0) then

445

if (lhe_version.lt.3d0) then

446

write(cform,'(a4,i2,a6)') '(a1,',max(nexternal,10),'e15.7)'

447

write(buff,cform) '#',(ptclus(i),i=3,nexternal)

448

else if(nexternal.gt.2)then

449

temp0='<scales '

450

temp=''

451

do i=3,nexternal

452

integfour=''

453

float=''

454

Write(float,'(f16.5)') ptclus(i)

455

write(integfour,'(i4)') ito(i)

456

temp=trim(temp)//' pt_clust_'//trim(adjustl(integfour))//'="'//trim(adjustl(float))//'"'

457

enddo

458

ptclusstring=trim(adjustl(temp0//trim(temp)//'></scales>'))

459

c write(*,*)'WRITING THE ptclusscale:',trim(adjustl(ptclusstring))

460

write(buff,'(a)') trim(adjustl(ptclusstring))

461

endif

462

endif

463

464

npart = nexternal+nres

465

466

return

467

end

468

469

c *************************************

470

subroutine write_error(ida1,ida2,imo)

471

c *************************************

472

implicit none

473

integer ida1,ida2,imo

474

475

open(unit=26,file='../../../error',status='unknown',err=999)

476

if (ida1.eq.1000)then

477

write(26,*) 'Error: too many particles in multipart. vertex,',

478

$ ' please increase maxcolmp in addmothers.f'

479

write(*,*) 'Error: too many particles in multipart. vertex,',

480

$ ' please increase maxcolmp in addmothers.f'

481

stop

482

endif

483

if (ida1.eq.1001)then

484

write(26,*) 'Error: failed to reduce to color indices: ',ida2,imo

485

write(*,*) 'Error: failed to reduce to color indices: ',ida2,imo

486

stop

487

endif

488

write(26,*) 'Error: Color combination ',ida1,ida2,

489

$ '->',imo,' not implemented in addmothers.f'

490

write(*,*) 'Error: Color combination ',ida1,ida2,

491

$ '->',imo,' not implemented in addmothers.f'

492

stop

493

494

999 write(*,*) 'error'

495

end

496

497

c *********************************************************************

498

function set_colmp(ncolmp,icolmp,npart,jpart,forest,icol,icol1,icol2)

499

c *********************************************************************

500

implicit none

501

integer maxcolmp

502

parameter(maxcolmp=20)

503

include 'nexternal.inc'

504

include 'genps.inc'

505

c Arguments

506

integer set_colmp

507

integer ncolmp,icolmp(2,*),npart,icol1,icol2

508

integer icol(2,-nexternal+2:2*nexternal-3)

509

integer jpart(7,-nexternal+3:2*nexternal-3)

510

integer forest(2,-max_branch:-1)

511

c Local

512

integer da_color(2),itmp,ida(2),icolor,ipart,i,j

513

integer get_color,set1colmp

514

515

set_colmp=ncolmp

516

icolor=get_color(jpart(1,npart))

517

if(icolor.eq.1) return

518

if(icolor.eq.2) then

519

c Multiparticle vertex - need to go through daughters and collect all colors

520

ipart=npart

521

do while(icolor.eq.2)

522

ida(1)=forest(1,ipart)

523

ida(2)=forest(2,ipart)

524

da_color(1)=get_color(jpart(1,ida(1)))

525

da_color(2)=get_color(jpart(1,ida(2)))

526

c print *,'iforest: ',ipart,ida(1),ida(2),da_color(1),da_color(2)

527

if(da_color(1).ne.2.and.da_color(2).lt.da_color(1).or.

528

$ da_color(2).eq.2)then

529

c Order daughters according to color, but always color 2 first

530

itmp=ida(1)

531

ida(1)=ida(2)

532

ida(2)=itmp

533

itmp=da_color(1)

534

da_color(1)=da_color(2)

535

da_color(2)=itmp

536

endif

537

do i=1,2

538

if(da_color(i).ne.1.and.da_color(i).ne.2)then

539

ncolmp=set1colmp(ncolmp,icolmp,icol(1,ida(i)),

540

$ icol(2,ida(i)))

541

endif

542

enddo

543

icolor=da_color(1)

544

ipart=ida(1)

545

enddo

546

else

547

ncolmp=set1colmp(ncolmp,icolmp,icol1,icol2)

548

endif

549

set_colmp=ncolmp

550

return

551

end

552

553

c ******************************************************

554

function set1colmp(ncolmp,icolmp,icol1,icol2)

555

c ******************************************************

556

implicit none

557

c Arguments

558

integer maxcolmp

559

parameter(maxcolmp=20)

560

integer set1colmp

561

integer ncolmp,icolmp(2,*),icol1,icol2,i,j

562

563

c print *,'icol1,icol2: ',icol1,icol2

564

565

ncolmp=ncolmp+1

566

icolmp(1,ncolmp)=icol1

567

icolmp(2,ncolmp)=icol2

568

c Avoid color sextet-type negative indices

569

if(icolmp(1,ncolmp).lt.0)then

570

ncolmp=ncolmp+1

571

icolmp(2,ncolmp)=-icolmp(1,ncolmp-1)

572

icolmp(1,ncolmp-1)=0

573

icolmp(1,ncolmp)=0

574

elseif(icolmp(2,ncolmp).lt.0)then

575

ncolmp=ncolmp+1

576

icolmp(1,ncolmp)=-icolmp(2,ncolmp-1)

577

icolmp(2,ncolmp-1)=0

578

icolmp(2,ncolmp)=0

579

endif

580

c print *,'icolmp: ',((icolmp(i,j),i=1,2),j=1,ncolmp)

581

if(ncolmp.gt.maxcolmp)

582

$ call write_error(1000,ncolmp,maxcolmp)

583

set1colmp=ncolmp

584

return

585

end

586

587

c********************************************************************

588

function fix_tchannel_color(mo_color,maxcolor,ncolmp,icolmp,ires,

589

$ icol,is_colors)

590

c********************************************************************

591

c Successively eliminate identical pairwise color indices from the

592

c icolmp list, until only (max) one triplet and one antitriplet remains

593

594

595

implicit none

596

include 'nexternal.inc'

597

integer fix_tchannel_color

598

integer mo_color,maxcolor,ncolmp,icolmp(2,*)

599

integer ires,icol(2,-nexternal+2:2*nexternal-3)

600

integer is_colors(2,nincoming)

601

integer i,j,i3,i3bar,max3,max3bar,min3,min3bar,maxcol,mincol

602

integer count

603

604

c Successively eliminate color indices in pairs until only the wanted

605

c indices remain

606

do i=1,ncolmp

607

do j=1,ncolmp

608

if(icolmp(1,i).ne.0.and.icolmp(1,i).eq.icolmp(2,j)) then

609

icolmp(1,i)=0

610

icolmp(2,j)=0

611

endif

612

enddo

613

enddo

614

615

i3=0

616

i3bar=0

617

icol(1,ires)=0

618

icol(2,ires)=0

619

do i=1,ncolmp

620

if(icolmp(1,i).gt.0)then

621

i3=i3+1

622

c color for t-channels needs to be reversed

623

if(i3.eq.1) icol(2,ires)=icolmp(1,i)

624

if(i3.eq.2) icol(1,ires)=-icolmp(1,i)

625

endif

626

if(icolmp(2,i).gt.0)then

627

i3bar=i3bar+1

628

c color for t-channels needs to be reversed

629

if(i3bar.eq.1) icol(1,ires)=icolmp(2,i)

630

if(i3bar.eq.2) icol(2,ires)=-icolmp(2,i)

631

endif

632

enddo

633

634

if(mo_color.eq.0)then

635

icol(1,ires)=0

636

icol(2,ires)=0

637

endif

638

639

fix_tchannel_color=maxcolor

640

if(mo_color.le.1.and.i3.eq.0.and.i3bar.eq.0) return

641

if(mo_color.eq.3.and.(i3.eq.1.and.i3bar.eq.0

642

$ .or.i3bar.eq.1.and.i3.eq.0)) return

643

if(mo_color.eq.8.and.i3.eq.1.and.i3bar.eq.1) return

644

645

c Make sure that max and min don't come from the same octet

646

call find_max_min(icolmp,ncolmp,max3,min3,max3bar,min3bar,

647

$ i3,i3bar,is_colors)

648

c print *,'After finding: ',ncolmp,((icolmp(j,i),j=1,2),i=1,ncolmp)

649

c print *,'mo_color = ',mo_color

650

651

if(mo_color.le.1.and.i3-i3bar.eq.2.or.

652

$ mo_color.le.1.and.i3bar-i3.eq.2.or.

653

$ mo_color.le.1.and.i3.eq.1.and.i3bar.eq.1) then

654

c Replace the maximum index with the minimum one everywhere

655

maxcol=max(max3,max3bar)

656

if(maxcol.eq.max3) then

657

mincol=min3bar

658

else

659

mincol=min3

660

endif

661

do i=ires+1,-1

662

do j=1,2

663

if(icol(j,i).eq.maxcol)

664

$ icol(j,i)=mincol

665

enddo

666

enddo

667

c print *,'Replaced ',maxcol,' by ',mincol

668

elseif(mo_color.le.1.and.i3.eq.2.and.i3bar.eq.2) then

669

c Ensure that max > min

670

if(max3bar.lt.min3)then

671

i=min3

672

min3=max3bar

673

max3bar=i

674

endif

675

if(max3.lt.min3bar)then

676

i=min3bar

677

min3bar=max3

678

max3=i

679

endif

680

c Replace the maximum indices with the minimum ones everywhere

681

do i=ires+1,-1

682

do j=1,2

683

if(icol(j,i).eq.max3bar)

684

$ icol(j,i)=min3

685

if(icol(j,i).eq.max3)

686

$ icol(j,i)=min3bar

687

enddo

688

enddo

689

c print *,'Replaced ',max3bar,' by ',min3,' and ',max3,' by ',min3bar

690

elseif(mo_color.le.1.and.mod(i3,3).eq.0.and.mod(i3bar,3).eq.0)then

691

c This is epsilon index - do nothing

692

continue

693

else if(mo_color.eq.3.and.mod(i3-i3bar,3).eq.2) then

694

c This is an epsilon index

695

maxcolor=maxcolor+1

696

icol(1,ires)=maxcolor

697

icol(2,ires)=0

698

c print *,'Set mother color for ',ires,' to ',(icol(j,ires),j=1,2)

699

else if(mo_color.eq.3.and.mod(i3bar-i3,3).eq.2) then

700

c This is an epsilon index

701

maxcolor=maxcolor+1

702

icol(1,ires)=0

703

icol(2,ires)=maxcolor

704

c print *,'Set mother color for ',ires,' to ',(icol(j,ires),j=1,2)

705

else if(mo_color.eq.3.and.(i3-i3bar.eq.1.or.i3bar-i3.eq.1).or.

706

$ mo_color.eq.8.and.i3.eq.2.and.i3bar.eq.2) then

707

c Replace the maximum index with the minimum one everywhere

708

c (we don't know if we should replace i3 with i3bar or vice versa)

709

c Actually we know if one of the index is repeated (we do not want to replace that one)

710

maxcol=max(max3,max3bar)

711

if(maxcol.eq.max3) then

712

mincol=min3bar

713

else

714

mincol=min3

715

endif

716

do i=ires+1,-1

717

do j=1,2

718

if(icol(j,i).eq.maxcol)

719

$ icol(j,i)=mincol

720

enddo

721

enddo

722

723

if (mincol.eq.0.and.mo_color.eq.3) then

724

c situation like (possible if they are epsilon in the gluon decay

725

c (503,0)----------+ggggggggggggg (509,508)

726

c |

727

c |(x,y)

728

c in this case maxcol=509 and mincol=0

729

c The correct solution in this case is:

730

c (503,0)----------+ggggggggggggg (503,508)

731

c |

732

c |(0,508)

733

if(icolmp(2,1).eq.0)then

734

maxcol = icolmp(1,2)

735

mincol = icolmp(1,1)

736

icol(1,ires) = 0

737

icol(2,ires) = icolmp(2,2)

738

elseif(icolmp(1,1).eq.0)then

739

maxcol = icolmp(2,2)

740

mincol = icolmp(2,1)

741

icol(1,ires) = icolmp(1,2)

742

icol(2,ires) = 0

743

elseif(icolmp(2,2).eq.0)then

744

maxcol = icolmp(1,1)

745

mincol = icolmp(1,2)

746

icol(1,ires) = 0

747

icol(2,ires) = icolmp(2,1)

748

elseif(icolmp(1,2).eq.0)then

749

maxcol = icolmp(2,1)

750

mincol = icolmp(2,2)

751

icol(1,ires) = icolmp(1,1)

752

icol(2,ires) = 0

753

else

754

!should not happen

755

write(*,*) "weird color combination in addmothers.f"

756

write(*,*) icolmp(1,1), icolmp(2,1), icolmp(1,2), icolmp(2,2)

757

call write_error(1001,0,0)

758

endif

759

c now maxcol=509 and mincol=503 so replace all occurence of 509-> 503

760

c print *,'Replaced ',maxcol,' by ',mincol

761

do i=ires+1,nexternal

762

do j=1,2

763

if(icol(j,i).eq.maxcol)

764

$ icol(j,i)=mincol

765

enddo

766

enddo

767

else

768

c standard case

769

c First check that mincol is not a going trough index. If it is it

770

C should not be assign to one of the temporary index

771

count=0

772

do i=1,nexternal

773

do j=1,2

774

if (icol(j,i).eq.mincol) count = count +1

775

enddo

776

enddo

777

778

if(count.eq.2)then

779

c we do not want to use that index pass to the other one

780

if (mincol.eq.min3)then

781

mincol = min3bar

782

maxcol = max3

783

else

784

mincol = min3

785

maxcol = max3bar

786

endif

787

endif

788

789

c Fix the color for ires (remember 3<->3bar for t-channels)

790

icol(1,ires)=0

791

icol(2,ires)=0

792

c print *,'Replaced ',maxcol,' by ',mincol

793

if(max3.eq.maxcol)then

794

if(i3-i3bar.ge.0) icol(2,ires)=min3

795

if(i3bar-i3.ge.0) icol(1,ires)=max3bar

796

else

797

if(i3-i3bar.ge.0) icol(2,ires)=max3

798

if(i3bar-i3.ge.0) icol(1,ires)=min3bar

799

endif

800

endif

801

c print *,'Set mother color for ',ires,' to ',(icol(j,ires),j=1,2)

802

else

803

c Don't know how to deal with this

804

call write_error(i3,i3bar,mo_color)

805

endif

806

fix_tchannel_color=maxcolor

807

808

return

809

end

810

811

c*******************************************************************

812

function elim_indices(n3,n3bar,ncolmp,icolmp,ires,icol,

813

$ is_colors,maxcolor)

814

c*******************************************************************

815

c Successively eliminate identical pairwise color indices from the

816

c icolmp list, until only the wanted indices remain

817

c n3 gives the number of triplet indices, n3bar number of antitriplets

818

c n3=1 for triplet, n3bar=1 for antitriplet,

819

c (n3,n3bar)=(1,1) for octet,

820

c n3=2 for sextet, n3bar=2 for antisextet

821

c If there are epsilon^{ijk} or epsilonbar color couplings, we

822

c need to introduce new index based on maxcolor.

823

824

825

implicit none

826

include 'nexternal.inc'

827

integer elim_indices

828

integer n3,n3bar,ncolmp,icolmp(2,*),maxcolor

829

integer ires,icol(2,-nexternal+2:2*nexternal-3)

830

integer is_colors(2,nincoming)

831

integer i,j,i3,i3bar

832

833

c Successively eliminate color indices in pairs until only the wanted

834

c indices remain

835

do i=1,ncolmp

836

do j=1,ncolmp

837

if(icolmp(1,i).ne.0.and.icolmp(1,i).eq.icolmp(2,j)) then

838

icolmp(1,i)=0

839

icolmp(2,j)=0

840

endif

841

enddo

842

enddo

843

844

i3=0

845

i3bar=0

846

icol(1,ires)=0

847

icol(2,ires)=0

848

do i=1,ncolmp

849

if(icolmp(1,i).gt.0)then

850

i3=i3+1

851

if(i3.eq.1) icol(1,ires)=icolmp(1,i)

852

if(i3.eq.2) icol(2,ires)=-icolmp(1,i)

853

endif

854

if(icolmp(2,i).gt.0)then

855

i3bar=i3bar+1

856

if(i3bar.eq.1) icol(2,ires)=icolmp(2,i)

857

if(i3bar.eq.2) icol(1,ires)=-icolmp(2,i)

858

endif

859

enddo

860

861

c print *,'i3,n3,i3bar,n3bar: ',i3,n3,i3bar,n3bar

862

c print *,'icol(1,ires),icol(2,ires): ',icol(1,ires),icol(2,ires)

863

864

if(n3bar.le.1.and.n3.eq.0) icol(1,ires)=0

865

if(n3.le.1.and.n3bar.eq.0) icol(2,ires)=0

866

867

if(i3.ne.n3.or.i3bar.ne.n3bar) then

868

if(n3.gt.0.and.n3bar.eq.0.and.mod(i3bar+n3,3).eq.0.and.i3.eq.0)then

869

c This is an epsilon index interaction

870

c write(*,*) i3, n3, i3bar, n3bar, ires

871

maxcolor=maxcolor+1

872

icol(1,ires)=maxcolor

873

if(n3.eq.2)then

874

maxcolor=maxcolor+1

875

icol(2,ires)=-maxcolor

876

endif

877

elseif(n3bar.gt.0.and.n3.eq.0.and.mod(i3+n3bar,3).eq.0.and.i3bar.eq.0)then

878

c This is an epsilonbar index interaction

879

c write(*,*) i3, n3, i3bar, n3bar, ires

880

maxcolor=maxcolor+1

881

icol(2,ires)=maxcolor

882

if(n3.eq.2)then

883

maxcolor=maxcolor+1

884

icol(1,ires)=-maxcolor

885

endif

886

elseif(n3.gt.0.and.n3bar.eq.0.and.i3-i3bar.eq.n3.or.

887

$ n3bar.gt.0.and.n3.eq.0.and.i3bar-i3.eq.n3bar.or.

888

$ n3.eq.1.and.n3bar.eq.1.and.i3-i3bar.eq.0.or.

889

$ n3.eq.0.and.n3bar.eq.0.and.i3-i3bar.eq.0.or.

890

$ n3bar.gt.0.and.n3.eq.0.and.mod(i3+n3bar,3).eq.0.and.i3bar.ne.0.or.

891

$ n3.gt.0.and.n3bar.eq.0.and.mod(i3bar+n3,3).eq.0.and.i3.ne.0)then

892

c We have a previous epsilon which gives the wrong pop-up index

893

call fix_s_color_indices(n3,n3bar,i3,i3bar,ncolmp,icolmp,

894

$ ires,icol,is_colors)

895

else

896

c Don't know how to deal with this

897

call write_error(1001,n3,n3bar)

898

endif

899

endif

900

901

elim_indices=maxcolor

902

903

return

904

end

905

906

c*******************************************************************

907

subroutine fix_s_color_indices(n3,n3bar,i3,i3bar,ncolmp,icolmp,

908

$ ires,icol,is_colors)

909

c*******************************************************************

910

911

c Fix color flow if some particle has got the wrong pop-up color

912

c due to epsilon-ijk vertices

913

914

915

implicit none

916

include 'nexternal.inc'

917

integer n3,n3bar,ncolmp,icolmp(2,*),maxcolor

918

integer ires,icol(2,-nexternal+2:2*nexternal-3)

919

integer is_colors(2,nincoming)

920

integer i,j,i3,i3bar

921

integer max_n3,max_n3bar,min_n3,min_n3bar,maxcol,mincol

922

923

icol(1,ires)=0

924

icol(2,ires)=0

925

926

c print *,'Colors: ',ncolmp,((icolmp(j,i),j=1,2),i=1,ncolmp)

927

c print *,'n3,n3bar,i3,i3bar: ',n3,n3bar,i3,i3bar

928

929

c Make sure that max and min don't come from the same octet

930

call find_max_min(icolmp,ncolmp,max_n3,min_n3,

931

$ max_n3bar,min_n3bar,i3,i3bar,is_colors)

932

c print *,'max3,min3bar,min3,max3bar: ',max_n3,min_n3bar,min_n3,max_n3bar

933

934

if(n3.eq.1.and.n3bar.eq.0.and.i3-i3bar.eq.n3.or.

935

$ n3bar.eq.1.and.n3.eq.0.and.i3bar-i3.eq.n3bar.or.

936

$ n3bar.eq.1.and.n3.eq.1.and.i3bar-i3.eq.0.or.

937

$ n3bar.eq.0.and.n3.eq.0.and.i3bar-i3.eq.0.or.

938

$ n3bar.gt.0.and.n3.eq.0.and.mod(i3+n3bar,3).eq.0.and.i3bar.ne.0.or.

939

$ n3.gt.0.and.n3bar.eq.0.and.mod(i3bar+n3,3).eq.0.and.i3.ne.0)then

940

941

if ((i3.eq.2.or.i3bar.eq.2).and.(n3bar+n3.eq.1))then

942

c Special case:

943

c -------------------- (0,503)

944

c g

945

c g

946

c g (504,505)

947

948

c need to correct to

949

c ------------------------- (0,503)

950

c (0,505) g

951

c g

952

c g (503,505)

953

if (i3.eq.2) then

954

icol(1,ires) = max(icolmp(1,1), icolmp(1,2))

955

icol(2,ires) = 0

956

maxcol = max(icolmp(2,1), icolmp(2,2))

957

mincol = min(icolmp(1,1), icolmp(1,2))

958

c replace maxcol by mincol

959

elseif (i3bar.eq.2) then

960

icol(1,ires) = 0

961

icol(2,ires) = max(icolmp(2,1), icolmp(2,2))

962

maxcol = max(icolmp(1,1), icolmp(1,2))

963

mincol = min(icolmp(2,1), icolmp(2,2))

964

endif

965

c write(*,*) "replace ", maxcol,"by",mincol

966

do i=ires+1,nexternal

967

do j=1,2

968

if(icol(j,i).eq.maxcol)

969

$ icol(j,i)=mincol

970

enddo

971

enddo

972

973

974

else

975

c Replace the highest 3bar-index with the lowest 3-index,

976

c or vice versa

977

978

979

maxcol=max(max_n3,max_n3bar)

980

if(maxcol.eq.max_n3) then

981

mincol=min_n3bar

982

else

983

mincol=min_n3

984

endif

985

do i=ires,-1

986

do j=1,2

987

if(icol(j,i).eq.maxcol)

988

$ icol(j,i)=mincol

989

enddo

990

enddo

991

c print *,'Replaced ',maxcol,' with ',mincol

992

if(max_n3.eq.maxcol)then

993

if(n3.eq.1) icol(1,ires)=min_n3

994

if(n3bar.eq.1) icol(2,ires)=max_n3bar

995

else

996

if(n3.eq.1) icol(1,ires)=max_n3

997

if(n3bar.eq.1) icol(2,ires)=min_n3bar

998

endif

999

c print *,'Set mother color for ',ires,' to ',(icol(j,ires),j=1,2)

1000

endif

1001

else

1002

c Don't know how to deal with this

1003

call write_error(1001,n3,n3bar)

1004

endif

1005

return

1006

end

1007

1008

1009

subroutine check_pure_internal_flow(icol,jpart, maxcolor)

1010

1011

implicit none

1012

include 'nexternal.inc'

1013

integer jpart(7,-nexternal+3:2*nexternal-3)

1014

integer icol(2,-nexternal+2:2*nexternal-3)

1015

integer maxcolor

1016

integer i,j,k,l

1017

logical found

1018

1019

c do i=-nexternal+3,nexternal

1020

c write(*,*) i, icol(1,i), icol(2,i),(jpart(j,i) , j=1,3)

1021

c enddo

1022

do i=-nexternal+3,-1

1023

if (jpart(2,i).eq.0.or.jpart(3,i).eq.0) goto 20 ! not define mother -> continue

1024

if (icol(1,i).eq.1000.or.icol(2,i).eq.1000) goto 20 ! special color value -> continue

1025

do k = 1,2

1026

found=.false.

1027

do j=1,nexternal

1028

if(icol(k,i).eq.icol(1,j).or.icol(k,i).eq.icol(2,j))then

1029

found=.true.

1030

goto 10 ! break

1031

endif

1032

enddo

1033

10 continue

1034

if (.not.found)then

1035

call correct_external_flow_epsilon(icol, jpart, maxcolor,

1036

& icol(k,i))

1037

endif

1038

enddo

1039

20 continue

1040

enddo

1041

return

1042

end

1043

1044

1045

1046

subroutine correct_external_flow_epsilon(icol, jpart, maxcolor, mincol)

1047

1048

c for avoiding double epsilon problem

1049

1050

implicit none

1051

include 'nexternal.inc'

1052

integer jpart(7,-nexternal+3:2*nexternal-3)

1053

integer maxcolor

1054

integer icol(2,-nexternal+2:2*nexternal-3)

1055

integer i,j,i3

1056

integer mincol ! the potential propagator linked to the two epsilon.

1057

integer k,l

1058

integer potential_index(2)

1059

integer epsilon_index(4)

1060

integer mothers(2*nexternal-3)

1061

logical to_change

1062

1063

C In presence of two epsilon_ijk connected by a flow we need to ensure that the

1064

C the index of the non summed indices do not repeat each other

1065

l=0

1066

do i=-nexternal+3,2*nexternal-3

1067

if (icol(1,i).eq.mincol.or.icol(2,i).eq.mincol)then

1068

potential_index(1)=0

1069

c write(*,*) "particle",i,"has color index", mincol

1070

k=0 !index to see how many child we found so far

1071

do j=-nexternal+3,2*nexternal-3

1072

if (jpart(2,j).eq.i.or.jpart(3,j).eq.i)then

1073

c write(*,*) "find the child", j

1074

if (icol(1,j).eq.mincol.or.icol(2,j).eq.mincol)then

1075

potential_index(1)=0

1076

c write(*,*) "the color", mincol,

1077

c & "is pass to one of the children ->no epsilon at this stage"

1078

c the color flow is pass to a child so no need to do anything for this part/junction

1079

goto 10 ! break

1080

elseif(icol(1,j).ne.0) then

1081

c write(*,*) "child has not colour", mincol, "add", icol(1,j)

1082

k = k+1

1083

potential_index(k) = icol(1,j)

1084

mothers(1) = i

1085

elseif(icol(2,j).ne.0)then

1086

c write(*,*) "child has not colour", mincol, "add", icol(2,j)

1087

k = k+1

1088

potential_index(k) = icol(2,j)

1089

mothers(1) = i

1090

else

1091

call write_error(1001,0,0)

1092

endif

1093

endif

1094

enddo

1095

10 continue

1096

c store the index of the final junction for this color

1097

c write(*,*) "found", potential_index

1098

if (potential_index(1).ne.0) then

1099

l = l+1

1100

epsilon_index(l) = potential_index(1)

1101

l = l+1

1102

epsilon_index(l) = potential_index(2)

1103

endif

1104

endif

1105

enddo

1106

C Remove the duplication index if any

1107

mothers(2) = 0

1108

c check the mother of mothers1 and change the color index

1109

c epsilon_index(3) -> maxcolor+1, epsilon_index(4) -> maxcolor+2

1110

c then add info on mothers to recursively change

1111

c Firs check if we have to change

1112

to_change = .false.

1113

do i=3,4

1114

do j=1,2

1115

if (epsilon_index(i).eq.epsilon_index(j))then

1116

C ` The index is duplicated need to change one

1117

to_change = .true.

1118

endif

1119

enddo

1120

enddo

1121

if (epsilon_index(4).eq.0) to_change = .false.

1122

if (to_change)then

1123

k=1

1124

do i =1, 2*nexternal-3

1125

if (mothers(i).eq.0)then

1126

goto 20 !break

1127

endif

1128

do j=mothers(i)+1,2*nexternal-3

1129

if (jpart(2,j).eq.mothers(i).or.jpart(3,j).eq.mothers(i))then

1130

if (icol(1,j).eq.epsilon_index(3))then

1131

icol(1,j) = maxcolor + 1

1132

k = k+1

1133

mothers(k) = j

1134

mothers(k+1) = 0

1135

elseif (icol(2,j).eq.epsilon_index(3))then

1136

icol(2,j) = maxcolor + 1

1137

k = k+1

1138

mothers(k) = j

1139

mothers(k+1) = 0

1140

elseif (icol(1,j).eq.epsilon_index(4))then

1141

icol(1,j) = maxcolor + 2

1142

k = k+1

1143

mothers(k) = j

1144

mothers(k+1) = 0

1145

elseif (icol(2,j).eq.epsilon_index(4))then

1146

icol(2,j) = maxcolor + 2

1147

k = k+1

1148

mothers(k) = j

1149

mothers(k+1) = 0

1150

endif

1151

endif

1152

enddo

1153

enddo

1154

20 continue

1155

maxcolor = maxcolor +2

1156

endif

1157

end

1158

1159

c*******************************************************************************

1160

subroutine find_max_min(icolmp,ncolmp,max3,min3,max3bar,min3bar,

1161

$ i3,i3bar,is_colors)

1162

c*******************************************************************************

1163

implicit none

1164

include 'nexternal.inc'

1165

integer ncolmp,icolmp(2,*)

1166

integer is_colors(2,nincoming)

1167

integer i,j,k,max3,max3bar,min3,min3bar,i3,i3bar,i3now,i3barnow

1168

integer allpairs(2,nexternal),npairs,maxcol,mincol

1169

1170

i3now=i3

1171

i3barnow=i3bar

1172

1173

c Create all possible pairs (3,3bar) that

1174

c 1. come from different octets, 2. are different,

1175

c 3. don't contain any color lines passing through the event

1176

npairs = 0

1177

do 20 i=1,ncolmp

1178

if(icolmp(1,i).eq.0) goto 20

1179

do k=1,nincoming

1180

if(icolmp(1,i).eq.is_colors(1,k)) goto 20

1181

enddo

1182

do 10 j=1,ncolmp

1183

if(i.eq.j.or.icolmp(2,j).eq.0.or.

1184

$ icolmp(1,i).eq.icolmp(2,j)) goto 10

1185

do k=1,nincoming

1186

if(icolmp(2,j).eq.is_colors(2,k)) goto 10

1187

enddo

1188

npairs=npairs+1

1189

allpairs(1,npairs)=icolmp(1,i)

1190

allpairs(2,npairs)=icolmp(2,j)

1191

10 enddo

1192

20 enddo

1193

1194

c print *,'is_colors: ',((is_colors(i,j),i=1,2),j=1,nincoming)

1195

c print *,'pairs: ',((allpairs(i,j),i=1,2),j=1,npairs)

1196

1197

c Find the pairs with maximum 3 and 3bar

1198

min3=1000

1199

min3bar=1000

1200

max3=0

1201

max3bar=0

1202

do i=1,npairs

1203

if(allpairs(1,i).gt.max3.and.

1204

$ (allpairs(2,i).lt.max3bar.or.

1205

$ allpairs(1,i).gt.allpairs(2,i)))then

1206

max3=allpairs(1,i)

1207

min3bar=allpairs(2,i)

1208

else if(allpairs(2,i).gt.max3bar.and.

1209

$ (allpairs(1,i).lt.max3.or.

1210

$ allpairs(2,i).gt.allpairs(1,i)))then

1211

max3bar=allpairs(2,i)

1212

min3=allpairs(1,i)

1213

endif

1214

enddo

1215

1216

c Find "maximum" pairs with minimum 3 and 3bar

1217

do i=1,npairs

1218

if(allpairs(1,i).eq.max3.and.

1219

$ allpairs(2,i).lt.min3bar.and.

1220

$ allpairs(2,i).ne.max3bar)

1221

$ min3bar=allpairs(2,i)

1222

if(allpairs(2,i).eq.max3bar.and.

1223

$ allpairs(1,i).lt.min3.and.

1224

$ allpairs(1,i).ne.max3)

1225

$ min3=allpairs(1,i)

1226

enddo

1227

1228

c Check that the pair are indeed different. Might not be the case if

1229

c The process contains some epsilon_ijk somewhere else.

1230

if (i3bar.gt.1.and.i3.gt.1)then

1231

if (max3bar.eq.min3bar)then

1232

c try to change min3bar

1233

min3bar=10000

1234

do i=1,npairs

1235

c search a new pair but with a different index!

1236

if(allpairs(1,i).eq.max3.and.

1237

$ allpairs(2,i).lt.min3bar.and.

1238

$ allpairs(2,i).ne.max3bar)then

1239

min3bar=allpairs(2,i)

1240

endif

1241

enddo

1242

c check if we found a new one. If not try to change the other index (max3bar)

1243

if (min3bar.eq.10000)then

1244

min3bar = max3bar

1245

max3bar = 0

1246

c search a new pair but with a different index!

1247

do i=1,npairs

1248

if(allpairs(1,i).eq.min3.and.

1249

$ allpairs(2,i).gt.max3bar.and.

1250

$ allpairs(2,i).ne.min3bar) then

1251

max3bar = allpairs(2,i)

1252

endif

1253

enddo

1254

c This should not happen.

1255

if (max3bar.eq.0)then

1256

write(*,*) "colorflow problem"

1257

call write_error(1001,0,0)

1258

endif

1259

endif

1260

endif

1261

c Doing the same but for the color index.

1262

if (max3.eq.min3)then

1263

c try to change min3

1264

min3=10000

1265

do i=1,npairs

1266

if(allpairs(2,i).eq.max3bar.and.

1267

$ allpairs(1,i).lt.min3.and.

1268

$ allpairs(1,i).ne.max3)then

1269

min3=allpairs(1,i)

1270

endif

1271

enddo

1272

if (min3.eq.10000)then

1273

min3 = max3 ! restore value

1274

max3 = 0

1275

c try to change max3

1276

do i=1,npairs

1277

if(allpairs(2,i).eq.min3bar.and.

1278

$ allpairs(1,i).gt.max3.and.

1279

$ allpairs(1,i).ne.min3) then

1280

max3 = allpairs(1,i)

1281

endif

1282

enddo

1283

if (max3.eq.0)then

1284

write(*,*) "colorflow problem"

1285

stop

1286

endif

1287

endif

1288

endif

1289

endif

1290

1291

if (max3.gt.0.and.max3bar.gt.0) then

1292

c We have found our two pairs, so we're done

1293

return

1294

endif

1295

1296

if(max3.gt.0.or.max3bar.gt.0)then

1297

i3now=i3now-1

1298

i3barnow=i3barnow-1

1299

endif

1300

1301

c Find pair for non-maximum (where we allow all colors)

1302

maxcol=max(max3,max3bar)

1303

if(maxcol.eq.max3) then

1304

mincol=min3bar

1305

else

1306

mincol=min3

1307

endif

1308

1309

npairs=0

1310

do i=1,ncolmp

1311

if(icolmp(1,i).eq.0.and.i3now.gt.0) cycle

1312

if(icolmp(1,i).eq.maxcol.or.icolmp(1,i).eq.mincol)

1313

$ cycle

1314

do j=1,ncolmp

1315

if(icolmp(2,j).eq.0.and.i3barnow.gt.0) cycle

1316

if(i.eq.j.or.icolmp(1,i).eq.icolmp(2,j)) cycle

1317

if(icolmp(2,j).eq.maxcol.or.icolmp(2,j).eq.mincol)

1318

$ cycle

1319

npairs=npairs+1

1320

allpairs(1,npairs)=icolmp(1,i)

1321

allpairs(2,npairs)=icolmp(2,j)

1322

enddo

1323

enddo

1324

if(npairs.ge.1)then

1325

if(maxcol.eq.max3)then

1326

min3=allpairs(1,1)

1327

max3bar=allpairs(2,1)

1328

else

1329

max3=allpairs(1,1)

1330

min3bar=allpairs(2,1)

1331

endif

1332

endif

1333

1334

c print *,'allpairs: ',((allpairs(i,j),i=1,2),j=1,npairs)

1335

1336

end

Older »