~maddevelopers/mg5amcnlo/2.0.2_onshellsubtract

tests/input_files/IOTestsComparison/IOTestsComparison/short_ML_SMQCD_default/gg_ttx/improve_ps.f

SUBROUTINE IMPROVE_PS_POINT_PRECISION(P)

IMPLICIT NONE

C CONSTANTS

INTEGER NEXTERNAL

PARAMETER (NEXTERNAL=4)

C ARGUMENTS

DOUBLE PRECISION P(0:3,NEXTERNAL)

REAL*16 QP_P(0:3,NEXTERNAL)

C LOCAL VARIABLES

INTEGER I,J

C ----------

C BEGIN CODE

C ----------

DO I=1,NEXTERNAL

DO J=0,3

QP_P(J,I)=P(J,I)

ENDDO

CALL MP_IMPROVE_PS_POINT_PRECISION(QP_P)

DO I=1,NEXTERNAL

DO J=0,3

P(J,I)=QP_P(J,I)

ENDDO

END

SUBROUTINE MP_IMPROVE_PS_POINT_PRECISION(P)

IMPLICIT NONE

C CONSTANTS

INTEGER NEXTERNAL

PARAMETER (NEXTERNAL=4)

C ARGUMENTS

REAL*16 P(0:3,NEXTERNAL)

C LOCAL VARIABLES

INTEGER I,J

INTEGER ERRCODE,ERRCODETMP

REAL*16 NEWP(0:3,NEXTERNAL)

C FUNCTIONS

LOGICAL MP_IS_PHYSICAL

C SAVED VARIABLES

INCLUDE 'MadLoopParams.inc'

C SAVED VARIABLES

INTEGER WARNED

DATA WARNED/0/

LOGICAL TOLD_SUPPRESS

DATA TOLD_SUPPRESS/.FALSE./

C ----------

C BEGIN CODE

C ----------

C ERROR CODES CONVENTION

C 1 :: None physical PS point input

C 100-1000 :: Error in the origianl method for restoring

C precision

C 1000-9999 :: Error when restoring precision ala PSMC

ERRCODETMP=0

ERRCODE=0

DO J=1,NEXTERNAL

DO I=0,3

NEWP(I,J)=P(I,J)

ENDDO

C Check the sanity of the original PS point

IF (.NOT.MP_IS_PHYSICAL(NEWP,WARNED)) THEN

ERRCODE = 1

WRITE(*,*) 'ERROR:: The input PS point is not precise enough.'

GOTO 100

ENDIF

C Now restore the precision

100

IF (IMPROVEPSPOINT.EQ.1) THEN

101

CALL MP_PSMC_IMPROVE_PS_POINT_PRECISION(NEWP,ERRCODE,WARNED)

102

ELSEIF((IMPROVEPSPOINT.EQ.2).OR.(IMPROVEPSPOINT.LE.0)) THEN

103

CALL MP_ORIG_IMPROVE_PS_POINT_PRECISION(NEWP,ERRCODE,WARNED)

104

ENDIF

105

IF (ERRCODE.NE.0) THEN

106

IF (WARNED.LT.20) THEN

107

WRITE(*,*) 'INFO:: Attempting to rescue the precision

108

$ improvement with an alternative method.'

109

WARNED=WARNED+1

110

ENDIF

111

IF (IMPROVEPSPOINT.EQ.1) THEN

112

CALL MP_ORIG_IMPROVE_PS_POINT_PRECISION(NEWP,ERRCODETMP

113

$ ,WARNED)

114

ELSEIF((IMPROVEPSPOINT.EQ.2).OR.(IMPROVEPSPOINT.LE.0)) THEN

115

CALL MP_PSMC_IMPROVE_PS_POINT_PRECISION(NEWP,ERRCODETMP

116

$ ,WARNED)

117

ENDIF

118

IF (ERRCODETMP.NE.0) GOTO 100

119

ENDIF

120

121

C Report to the user or update the PS point.

122

123

GOTO 101

124

100 CONTINUE

125

IF (WARNED.LT.20) THEN

126

WRITE(*,*) 'WARNING:: This PS point could not be improved.

127

$ Error code = ',ERRCODE,ERRCODETMP

128

CALL MP_WRITE_MOM(P)

129

WARNED = WARNED +1

130

ENDIF

131

GOTO 102

132

101 CONTINUE

133

DO J=1,NEXTERNAL

134

DO I=0,3

135

P(I,J)=NEWP(I,J)

136

ENDDO

137

ENDDO

138

102 CONTINUE

139

140

IF (WARNED.GE.20.AND..NOT.TOLD_SUPPRESS) THEN

141

WRITE(*,*) 'INFO:: Further warnings from the improve_ps

142

$ routine will now be supressed.'

143

TOLD_SUPPRESS=.TRUE.

144

ENDIF

145

146

END

147

148

149

FUNCTION MP_IS_CLOSE(P,NEWP,WARNED)

150

IMPLICIT NONE

151

152

C CONSTANTS

153

154

INTEGER NEXTERNAL

155

PARAMETER (NEXTERNAL=4)

156

REAL*16 ZERO

157

PARAMETER (ZERO=0.0E+00_16)

158

REAL*16 THRS_CLOSE

159

PARAMETER (THRS_CLOSE=1.0E-02_16)

160

161

C ARGUMENTS

162

163

REAL*16 P(0:3,NEXTERNAL), NEWP(0:3,NEXTERNAL)

164

LOGICAL MP_IS_CLOSE

165

INTEGER WARNED

166

167

C LOCAL VARIABLES

168

169

INTEGER I,J

170

REAL*16 REF,REF2

171

172

C NOW MAKE SURE THE SHIFTED POINT IS NOT TOO FAR FROM THE ORIGINAL

173

C ONE

174

MP_IS_CLOSE = .TRUE.

175

REF = ZERO

176

REF2 = ZERO

177

DO J=1,NEXTERNAL

178

DO I=0,3

179

REF2 = REF2 + ABS(P(I,J))

180

REF = REF + ABS(P(I,J)-NEWP(I,J))

181

ENDDO

182

ENDDO

183

184

IF ((REF/REF2).GT.THRS_CLOSE) THEN

185

MP_IS_CLOSE = .FALSE.

186

IF (WARNED.LT.20) THEN

187

WRITE(*,*) 'WARNING:: The improved PS point is too far from

188

$ the original one',(REF/REF2)

189

WARNED=WARNED+1

190

ENDIF

191

ENDIF

192

193

END

194

195

FUNCTION MP_IS_PHYSICAL(P,WARNED)

196

IMPLICIT NONE

197

198

C CONSTANTS

199

200

INTEGER NEXTERNAL

201

PARAMETER (NEXTERNAL=4)

202

INTEGER NINITIAL

203

PARAMETER (NINITIAL=2)

204

REAL*16 ZERO

205

PARAMETER (ZERO=0.0E+00_16)

206

REAL*16 MP__ZERO

207

PARAMETER (MP__ZERO=ZERO)

208

REAL*16 ONE

209

PARAMETER (ONE=1.0E+00_16)

210

REAL*16 TWO

211

PARAMETER (TWO=2.0E+00_16)

212

REAL*16 THRES_ONSHELL

213

PARAMETER (THRES_ONSHELL=1.0E-02_16)

214

REAL*16 THRES_FOURMOM

215

PARAMETER (THRES_FOURMOM=1.0E-06_16)

216

217

C ARGUMENTS

218

219

REAL*16 P(0:3,NEXTERNAL)

220

LOGICAL MP_IS_PHYSICAL

221

INTEGER WARNED

222

223

C LOCAL VARIABLES

224

225

INTEGER I,J

226

REAL*16 BUFF,REF

227

REAL*16 MASSES(NEXTERNAL)

228

229

C GLOBAL VARIABLES

230

231

INCLUDE 'mp_coupl.inc'

232

233

MASSES(1)=MP__ZERO

234

MASSES(2)=MP__ZERO

235

MASSES(3)=MP__MDL_MT

236

MASSES(4)=MP__MDL_MT

237

238

C ----------

239

C BEGIN CODE

240

C ----------

241

242

MP_IS_PHYSICAL = .TRUE.

243

244

C WE FIRST CHECK THAT THE INPUT PS POINT IS REASONABLY PHYSICAL

245

C FOR THAT WE NEED A REFERENCE SCALE

246

REF=ZERO

247

DO J=1,NEXTERNAL

248

REF=REF+ABS(P(0,J))

249

ENDDO

250

DO I=0,3

251

BUFF=ZERO

252

DO J=1,NINITIAL

253

BUFF=BUFF-P(I,J)

254

ENDDO

255

DO J=NINITIAL+1,NEXTERNAL

256

BUFF=BUFF+P(I,J)

257

ENDDO

258

IF ((BUFF/REF).GT.THRES_FOURMOM) THEN

259

IF (WARNED.LT.20) THEN

260

WRITE(*,*) 'ERROR:: Four-momentum conservation is not

261

$ accurate enough, ',(BUFF/REF)

262

CALL MP_WRITE_MOM(P)

263

WARNED=WARNED+1

264

ENDIF

265

MP_IS_PHYSICAL = .FALSE.

266

ENDIF

267

ENDDO

268

REF = REF / (ONE*NEXTERNAL)

269

DO I=1,NEXTERNAL

270

REF=ABS(P(0,I))+ABS(P(1,I))+ABS(P(2,I))+ABS(P(3,I))

271

IF ((SQRT(ABS(P(0,I)**2-P(1,I)**2-P(2,I)**2-P(3,I)**2

272

$ -MASSES(I)**2))/REF).GT.THRES_ONSHELL) THEN

273

IF (WARNED.LT.20) THEN

274

WRITE(*,*) 'ERROR:: Onshellness of the momentum of

275

$ particle ',I,' of mass ',MASSES(I),' is not accurate

276

$ enough, ', (SQRT(ABS(P(0,I)**2-P(1,I)**2-P(2,I)**2

277

$ -P(3,I)**2-MASSES(I)**2))/REF)

278

CALL MP_WRITE_MOM(P)

279

WARNED=WARNED+1

280

ENDIF

281

MP_IS_PHYSICAL = .FALSE.

282

ENDIF

283

ENDDO

284

285

END

286

287

SUBROUTINE WRITE_MOM(P)

288

IMPLICIT NONE

289

INTEGER NEXTERNAL

290

PARAMETER (NEXTERNAL=4)

291

INTEGER NINITIAL

292

PARAMETER (NINITIAL=2)

293

DOUBLE PRECISION ZERO

294

PARAMETER (ZERO=0.0D0)

295

DOUBLE PRECISION MDOT

296

297

INTEGER I,J

298

299

300

C ARGUMENTS

301

302

DOUBLE PRECISION P(0:3,NEXTERNAL),PSUM(0:3)

303

DO I=0,3

304

PSUM(I)=ZERO

305

DO J=1,NINITIAL

306

PSUM(I)=PSUM(I)+P(I,J)

307

ENDDO

308

DO J=NINITIAL+1,NEXTERNAL

309

PSUM(I)=PSUM(I)-P(I,J)

310

ENDDO

311

ENDDO

312

WRITE (*,*) ' Phase space point:'

313

WRITE (*,*) ' ---------------------'

314

WRITE (*,*) ' E | px | py | pz | m '

315

DO I=1,NEXTERNAL

316

WRITE (*,'(1x,5e27.17)') P(0,I),P(1,I),P(2,I),P(3,I),SQRT(ABS(M

317

$ DOT(P(0,I),P(0,I))))

318

ENDDO

319

WRITE (*,*) ' Four-momentum conservation sum:'

320

WRITE (*,'(1x,4e27.17)') PSUM(0),PSUM(1),PSUM(2),PSUM(3)

321

WRITE (*,*) ' ---------------------'

322

END

323

324

DOUBLE PRECISION FUNCTION MDOT(P1,P2)

325

IMPLICIT NONE

326

DOUBLE PRECISION P1(0:3),P2(0:3)

327

MDOT=P1(0)*P2(0)-P1(1)*P2(1)-P1(2)*P2(2)-P1(3)*P2(3)

328

RETURN

329

END

330

331

SUBROUTINE MP_WRITE_MOM(P)

332

IMPLICIT NONE

333

INTEGER NEXTERNAL

334

PARAMETER (NEXTERNAL=4)

335

INTEGER NINITIAL

336

PARAMETER (NINITIAL=2)

337

REAL*16 ZERO

338

PARAMETER (ZERO=0.0E+00_16)

339

REAL*16 MP_MDOT

340

341

INTEGER I,J

342

343

344

C ARGUMENTS

345

346

REAL*16 P(0:3,NEXTERNAL),PSUM(0:3),DOT

347

DOUBLE PRECISION DP_P(0:3,NEXTERNAL),DP_PSUM(0:3),DP_DOT

348

349

DO I=0,3

350

PSUM(I)=ZERO

351

DO J=1,NINITIAL

352

PSUM(I)=PSUM(I)+P(I,J)

353

ENDDO

354

DO J=NINITIAL+1,NEXTERNAL

355

PSUM(I)=PSUM(I)-P(I,J)

356

ENDDO

357

ENDDO

358

359

C The GCC4.7 compiler on SLC machines has trouble to write out

360

C quadruple precision variable with the write(*,*) statement. I

361

C therefore perform the cast by hand

362

DO I=0,3

363

DP_PSUM(I)=PSUM(I)

364

DO J=1,NEXTERNAL

365

DP_P(I,J)=P(I,J)

366

ENDDO

367

ENDDO

368

369

WRITE (*,*) ' Phase space point:'

370

WRITE (*,*) ' ---------------------'

371

WRITE (*,*) ' E | px | py | pz | m '

372

DO I=1,NEXTERNAL

373

DOT=SQRT(ABS(MP_MDOT(P(0,I),P(0,I))))

374

DP_DOT=DOT

375

WRITE (*,'(1x,5e27.17)') DP_P(0,I),DP_P(1,I),DP_P(2,I),DP_P(3

376

$ ,I),DP_DOT

377

ENDDO

378

WRITE (*,*) ' Four-momentum conservation sum:'

379

WRITE (*,'(1x,4e27.17)') DP_PSUM(0),DP_PSUM(1),DP_PSUM(2)

380

$ ,DP_PSUM(3)

381

WRITE (*,*) ' ---------------------'

382

END

383

384

REAL*16 FUNCTION MP_MDOT(P1,P2)

385

IMPLICIT NONE

386

REAL*16 P1(0:3),P2(0:3)

387

MP_MDOT=P1(0)*P2(0)-P1(1)*P2(1)-P1(2)*P2(2)-P1(3)*P2(3)

388

RETURN

389

END

390

391

C Rotate_PS rotates the PS point PS (without modifying it)

392

C stores the result in P and for the quadruple precision

393

C version , it also modifies the global variables

394

C PS and MP_DONE accordingly.

395

396

SUBROUTINE ROTATE_PS(P_IN,P,ROTATION)

397

IMPLICIT NONE

398

399

C CONSTANTS

400

401

INTEGER NEXTERNAL

402

PARAMETER (NEXTERNAL=4)

403

404

C ARGUMENTS

405

406

DOUBLE PRECISION P_IN(0:3,NEXTERNAL),P(0:3,NEXTERNAL)

407

INTEGER ROTATION

408

409

C LOCAL VARIABLES

410

411

INTEGER I,J

412

413

C ----------

414

C BEGIN CODE

415

C ----------

416

417

DO I=1,NEXTERNAL

418

C rotation=1 => (xp=z,yp=-x,zp=-y)

419

IF(ROTATION.EQ.1) THEN

420

P(0,I)=P_IN(0,I)

421

P(1,I)=P_IN(3,I)

422

P(2,I)=-P_IN(1,I)

423

P(3,I)=-P_IN(2,I)

424

C rotation=2 => (xp=-z,yp=y,zp=x)

425

ELSEIF(ROTATION.EQ.2) THEN

426

P(0,I)=P_IN(0,I)

427

P(1,I)=-P_IN(3,I)

428

P(2,I)=P_IN(2,I)

429

P(3,I)=P_IN(1,I)

430

ELSE

431

P(0,I)=P_IN(0,I)

432

P(1,I)=P_IN(1,I)

433

P(2,I)=P_IN(2,I)

434

P(3,I)=P_IN(3,I)

435

ENDIF

436

ENDDO

437

438

END

439

440

441

SUBROUTINE MP_ROTATE_PS(P_IN,P,ROTATION)

442

IMPLICIT NONE

443

444

C CONSTANTS

445

446

INTEGER NEXTERNAL

447

PARAMETER (NEXTERNAL=4)

448

449

C ARGUMENTS

450

451

REAL*16 P_IN(0:3,NEXTERNAL),P(0:3,NEXTERNAL)

452

INTEGER ROTATION

453

454

C LOCAL VARIABLES

455

456

INTEGER I,J

457

458

C GLOBAL VARIABLES

459

460

LOGICAL MP_DONE

461

COMMON/MP_DONE/MP_DONE

462

463

C ----------

464

C BEGIN CODE

465

C ----------

466

467

DO I=1,NEXTERNAL

468

C rotation=1 => (xp=z,yp=-x,zp=-y)

469

IF(ROTATION.EQ.1) THEN

470

P(0,I)=P_IN(0,I)

471

P(1,I)=P_IN(3,I)

472

P(2,I)=-P_IN(1,I)

473

P(3,I)=-P_IN(2,I)

474

C rotation=2 => (xp=-z,yp=y,zp=x)

475

ELSEIF(ROTATION.EQ.2) THEN

476

P(0,I)=P_IN(0,I)

477

P(1,I)=-P_IN(3,I)

478

P(2,I)=P_IN(2,I)

479

P(3,I)=P_IN(1,I)

480

ELSE

481

P(0,I)=P_IN(0,I)

482

P(1,I)=P_IN(1,I)

483

P(2,I)=P_IN(2,I)

484

P(3,I)=P_IN(3,I)

485

ENDIF

486

ENDDO

487

488

MP_DONE = .FALSE.

489

490

END

491

492

C *****************************************************************

493

C Beginning of the routine for restoring precision with V.H. method

494

C *****************************************************************

495

496

SUBROUTINE MP_ORIG_IMPROVE_PS_POINT_PRECISION(P,ERRCODE,WARNED)

497

IMPLICIT NONE

498

499

C CONSTANTS

500

501

INTEGER NEXTERNAL

502

PARAMETER (NEXTERNAL=4)

503

INTEGER NINITIAL

504

PARAMETER (NINITIAL=2)

505

REAL*16 ZERO

506

PARAMETER (ZERO=0.0E+00_16)

507

REAL*16 MP__ZERO

508

PARAMETER (MP__ZERO=ZERO)

509

REAL*16 ONE

510

PARAMETER (ONE=1.0E+00_16)

511

REAL*16 TWO

512

PARAMETER (TWO=2.0E+00_16)

513

REAL*16 THRS_TEST

514

PARAMETER (THRS_TEST=1.0E-15_16)

515

516

C ARGUMENTS

517

518

REAL*16 P(0:3,NEXTERNAL)

519

INTEGER ERRCODE, WARNED

520

521

C FUNCTIONS

522

523

LOGICAL MP_IS_CLOSE

524

525

C LOCAL VARIABLES

526

527

INTEGER I,J, P1, P2

528

C PT STANDS FOR PTOT

529

REAL*16 PT(0:3), NEWP(0:3,NEXTERNAL)

530

REAL*16 BUFF,REF,REF2,DISCR

531

REAL*16 MASSES(NEXTERNAL)

532

REAL*16 SHIFTE(2),SHIFTZ(2)

533

534

C GLOBAL VARIABLES

535

536

INCLUDE 'mp_coupl.inc'

537

538

MASSES(1)=MP__ZERO

539

MASSES(2)=MP__ZERO

540

MASSES(3)=MP__MDL_MT

541

MASSES(4)=MP__MDL_MT

542

543

C ----------

544

C BEGIN CODE

545

C ----------

546

ERRCODE = 0

547

548

C NOW WE MAKE SURE THAT THE PS POINT CAN BE IMPROVED BY THE

549

C ALGORITHM

550

REF=ZERO

551

DO J=1,NEXTERNAL

552

REF=REF+ABS(P(0,J))

553

ENDDO

554

555

IF (NINITIAL.NE.2) ERRCODE = 100

556

557

IF (ABS(P(1,1)/REF).GT.THRS_TEST.OR.ABS(P(2,1)/REF).GT.THRS_TEST.

558

$ OR.ABS(P(1,2)/REF).GT.THRS_TEST.OR.ABS(P(2,2)/REF).GT.THRS_TEST

559

$ ) ERRCODE = 200

560

561

IF (MASSES(1).NE.ZERO.OR.MASSES(2).NE.ZERO) ERRCODE = 300

562

563

DO I=1,NEXTERNAL

564

IF (P(0,I).LT.ZERO) ERRCODE = 400 + I

565

ENDDO

566

567

IF (ERRCODE.NE.0) GOTO 100

568

569

C WE FIRST SHIFT ALL THE FINAL STATE PARTICLES TO MAKE THEM

570

C EXACTLY ONSHELL

571

572

DO I=0,3

573

PT(I)=ZERO

574

ENDDO

575

DO I=NINITIAL+1,NEXTERNAL

576

DO J=0,3

577

IF (J.EQ.3) THEN

578

NEWP(3,I)=SIGN(SQRT(ABS(P(0,I)**2-P(1,I)**2-P(2,I)**2

579

$ -MASSES(I)**2)),P(3,I))

580

ELSE

581

NEWP(J,I)=P(J,I)

582

ENDIF

583

PT(J)=PT(J)+NEWP(J,I)

584

ENDDO

585

ENDDO

586

587

C WE CHOOSE P1 IN THE ALGORITHM TO ALWAYS BE THE PARTICLE WITH

588

C POSITIVE PZ

589

IF (P(3,1).GT.ZERO) THEN

590

P1=1

591

P2=2

592

ELSEIF (P(3,2).GT.ZERO) THEN

593

P1=2

594

P2=1

595

ELSE

596

ERRCODE = 500

597

GOTO 100

598

ENDIF

599

600

C Now we calculate the shift to bring to P1 and P2

601

C Mathematica gives

602

C ptotC = {ptotE, ptotX, ptotY, ptotZ};

603

C pm1C = {pm1E + sm1E, pm1X, pm1Y, pm1Z + sm1Z};

604

C {pm0E + sm0E, ptotX - pm1X, ptotY - pm1Y, pm0Z + sm0Z};

605

C sol = Solve[{ptotC[[1]] - pm1C[[1]] - pm0C[[1]] == 0,

606

C ptotC[[4]] - pm1C[[4]] - pm0C[[4]] == 0,

607

C pm1C[[1]]^2 - pm1C[[2]]^2 - pm1C[[3]]^2 - pm1C[[4]]^2 == m1M^2,

608

C pm0C[[1]]^2 - pm0C[[2]]^2 - pm0C[[3]]^2 - pm0C[[4]]^2 == m2M^2},

609

C {sm1E, sm1Z, sm0E, sm0Z}] // FullSimplify;

610

C (solC[[1]] /. {m1M -> 0, m2M -> 0} /. {pm1X -> 0, pm1Y -> 0})

611

C END

612

613

DISCR = -PT(0)**2 + PT(1)**2 + PT(2)**2 + PT(3)**2

614

IF (DISCR.LT.ZERO) DISCR = -DISCR

615

616

SHIFTE(1) = (PT(0)*(-TWO*P(0,P1)*PT(0) + PT(0)**2 + PT(1)**2

617

$ + PT(2)**2) + (TWO*P(0,P1) - PT(0))*PT(3)**2 + PT(3)*DISCR)

618

$ /(TWO*(PT(0) - PT(3))*(PT(0) + PT(3)))

619

SHIFTE(2) = -(PT(0)*(TWO*P(0,P2)*PT(0) - PT(0)**2 + PT(1)**2

620

$ + PT(2)**2) + (-TWO*P(0,P2) + PT(0))*PT(3)**2 + PT(3)*DISCR)

621

$ /(TWO*(PT(0) - PT(3))*(PT(0) + PT(3)))

622

SHIFTZ(1) = (-TWO*P(3,P1)*(PT(0)**2 - PT(3)**2) + PT(3)*(PT(0)*

623

$ *2 + PT(1)**2 + PT(2)**2 - PT(3)**2) + PT(0)*DISCR)/(TWO*(PT(0)

624

$ **2 - PT(3)**2))

625

SHIFTZ(2) = -(TWO*P(3,P2)*(PT(0)**2 - PT(3)**2) + PT(3)*(

626

$ -PT(0)**2 + PT(1)**2 + PT(2)**2 + PT(3)**2) + PT(0)*DISCR)/(TWO

627

$ *(PT(0)**2 - PT(3)**2))

628

NEWP(0,P1) = P(0,P1)+SHIFTE(1)

629

NEWP(3,P1) = P(3,P1)+SHIFTZ(1)

630

NEWP(0,P2) = P(0,P2)+SHIFTE(2)

631

NEWP(3,P2) = P(3,P2)+SHIFTZ(2)

632

NEWP(1,P2) = P(1,P2)

633

NEWP(2,P2) = P(2,P2)

634

DO J=1,2

635

REF=ZERO

636

DO I=NINITIAL+1,NEXTERNAL

637

REF = REF + P(J,I)

638

ENDDO

639

REF = REF - P(J,P2)

640

NEWP(J,P1) = REF

641

ENDDO

642

643

IF (.NOT.MP_IS_CLOSE(P,NEWP,WARNED)) THEN

644

ERRCODE=999

645

GOTO 100

646

ENDIF

647

648

DO J=1,NEXTERNAL

649

DO I=0,3

650

P(I,J)=NEWP(I,J)

651

ENDDO

652

ENDDO

653

654

100 CONTINUE

655

656

END

657

658

C *****************************************************************

659

C Beginning of the routine for restoring precision a la PSMC

660

C *****************************************************************

661

662

SUBROUTINE MP_PSMC_IMPROVE_PS_POINT_PRECISION(P,ERRCODE,WARNED)

663

IMPLICIT NONE

664

665

C CONSTANTS

666

667

INTEGER NEXTERNAL

668

PARAMETER (NEXTERNAL=4)

669

INTEGER NINITIAL

670

PARAMETER (NINITIAL=2)

671

REAL*16 ZERO

672

PARAMETER (ZERO=0.0E+00_16)

673

REAL*16 MP__ZERO

674

PARAMETER (MP__ZERO=ZERO)

675

REAL*16 ONE

676

PARAMETER (ONE=1.0E+00_16)

677

REAL*16 TWO

678

PARAMETER (TWO=2.0E+00_16)

679

REAL*16 CONSISTENCY_THRES

680

PARAMETER (CONSISTENCY_THRES=1.0E-25_16)

681

682

INTEGER NAPPROXZEROS

683

PARAMETER (NAPPROXZEROS=3)

684

685

686

C ARGUMENTS

687

688

REAL*16 P(0:3,NEXTERNAL)

689

INTEGER ERRCODE,ERROR,WARNED

690

691

C FUNCTIONS

692

693

LOGICAL MP_IS_CLOSE

694

695

C LOCAL VARIABLES

696

697

INTEGER I,J, P1, P2

698

REAL*16 NEWP(0:3,NEXTERNAL), PBUFF(0:3)

699

REAL*16 BUFF, BUFF2, XSCALE, APPROX_ZEROS(NAPPROXZEROS)

700

REAL*16 MASSES(NEXTERNAL)

701

702

C GLOBAL VARIABLES

703

704

INCLUDE 'mp_coupl.inc'

705

706

C ----------

707

C BEGIN CODE

708

C ----------

709

710

MASSES(1)=MP__ZERO

711

MASSES(2)=MP__ZERO

712

MASSES(3)=MP__MDL_MT

713

MASSES(4)=MP__MDL_MT

714

715

ERRCODE = 0

716

XSCALE = ONE

717

718

C Define the seeds which should be tried

719

APPROX_ZEROS(1)=1.0E+00_16

720

APPROX_ZEROS(2)=1.1E+00_16

721

APPROX_ZEROS(3)=0.9E+00_16

722

723

C Start by copying the momenta

724

DO I=1,NEXTERNAL

725

DO J=0,3

726

NEWP(J,I)=P(J,I)

727

ENDDO

728

ENDDO

729

730

C First make sur that the space like momentum is exactly conserved

731

DO J=0,3

732

PBUFF(J)=ZERO

733

ENDDO

734

DO I=1,NINITIAL

735

DO J=1,3

736

PBUFF(J)=PBUFF(J)+NEWP(J,I)

737

ENDDO

738

ENDDO

739

DO I=NINITIAL+1,NEXTERNAL-1

740

DO J=1,3

741

PBUFF(J)=PBUFF(J)-NEWP(J,I)

742

ENDDO

743

ENDDO

744

DO J=1,3

745

NEWP(J,NEXTERNAL)=PBUFF(J)

746

ENDDO

747

748

C Now find the 'x' rescaling factor

749

DO I=1,NAPPROXZEROS

750

CALL FINDX(NEWP,APPROX_ZEROS(I),XSCALE,ERROR)

751

IF(ERROR.EQ.0) THEN

752

GOTO 1001

753

ELSE

754

ERRCODE=ERRCODE+(10**(I-1))*ERROR

755

ENDIF

756

ENDDO

757

IF (WARNED.LT.20) THEN

758

WRITE(*,*) 'WARNING:: Could not find the proper rescaling

759

$ factor x. Restoring precision ala PSMC will therefore not be

760

$ used.'

761

WARNED=WARNED+1

762

ENDIF

763

IF (ERRCODE.LT.1000) THEN

764

ERRCODE=ERRCODE+1000

765

ENDIF

766

GOTO 1000

767

1001 CONTINUE

768

ERRCODE = 0

769

770

C Apply the rescaling

771

DO I=1,NEXTERNAL

772

DO J=1,3

773

NEWP(J,I)=NEWP(J,I)*XSCALE

774

ENDDO

775

ENDDO

776

777

C Now restore exact onshellness of the particles.

778

DO I=1,NEXTERNAL

779

BUFF=MASSES(I)**2

780

DO J=1,3

781

BUFF=BUFF+NEWP(J,I)**2

782

ENDDO

783

NEWP(0,I)=SQRT(BUFF)

784

ENDDO

785

786

C Consistency check

787

BUFF=ZERO

788

BUFF2=ZERO

789

DO I=1,NINITIAL

790

BUFF=BUFF-NEWP(0,I)

791

BUFF2=BUFF2+NEWP(0,I)

792

ENDDO

793

DO I=NINITIAL+1,NEXTERNAL

794

BUFF=BUFF+NEWP(0,I)

795

BUFF2=BUFF2+NEWP(0,I)

796

ENDDO

797

IF ((ABS(BUFF)/BUFF2).GT.CONSISTENCY_THRES) THEN

798

IF (WARNED.LT.20) THEN

799

WRITE(*,*) 'WARNING:: The consistency check in the a la PSMC

800

$ precision restoring algorithm failed. The result will

801

$ therefore not be used.'

802

WARNED=WARNED+1

803

ENDIF

804

ERRCODE = 1000

805

GOTO 1000

806

ENDIF

807

808

IF (.NOT.MP_IS_CLOSE(P,NEWP,WARNED)) THEN

809

ERRCODE=999

810

GOTO 1000

811

ENDIF

812

813

DO J=1,NEXTERNAL

814

DO I=0,3

815

P(I,J)=NEWP(I,J)

816

ENDDO

817

ENDDO

818

819

1000 CONTINUE

820

821

END

822

823

824

SUBROUTINE FINDX(P,SEED,XSCALE,ERROR)

825

IMPLICIT NONE

826

827

C CONSTANTS

828

829

INTEGER NEXTERNAL

830

PARAMETER (NEXTERNAL=4)

831

INTEGER NINITIAL

832

PARAMETER (NINITIAL=2)

833

REAL*16 ZERO

834

PARAMETER (ZERO=0.0E+00_16)

835

REAL*16 MP__ZERO

836

PARAMETER (MP__ZERO=ZERO)

837

REAL*16 ONE

838

PARAMETER (ONE=1.0E+00_16)

839

REAL*16 TWO

840

PARAMETER (TWO=2.0E+00_16)

841

INTEGER MAXITERATIONS

842

PARAMETER (MAXITERATIONS=8)

843

REAL*16 CONVERGED

844

PARAMETER (CONVERGED=1.0E-26_16)

845

846

C ARGUMENTS

847

848

REAL*16 P(0:3,NEXTERNAL),SEED,XSCALE

849

INTEGER ERROR

850

851

C LOCAL VARIABLES

852

853

INTEGER I,J,ERR

854

REAL*16 PVECSQ(NEXTERNAL)

855

REAL*16 XN, XNP1,FVAL,DVAL

856

857

C ----------

858

C BEGIN CODE

859

C ----------

860

861

ERROR = 0

862

XSCALE = SEED

863

XN = SEED

864

XNP1 = SEED

865

866

DO I=1,NEXTERNAL

867

PVECSQ(I)=P(1,I)**2+P(2,I)**2+P(3,I)**2

868

ENDDO

869

870

DO I=1,MAXITERATIONS

871

CALL FUNCT(PVECSQ(1),XN,.FALSE.,ERR, FVAL)

872

IF (ERR.NE.0) THEN

873

ERROR=ERR

874

GOTO 710

875

ENDIF

876

CALL FUNCT(PVECSQ(1),XN,.TRUE.,ERR, DVAL)

877

IF (ERR.NE.0) THEN

878

ERROR=ERR

879

GOTO 710

880

ENDIF

881

XNP1=XN-(FVAL/DVAL)

882

IF((ABS(((XNP1-XN)*TWO)/(XNP1+XN))).LT.CONVERGED) THEN

883

XN=XNP1

884

GOTO 700

885

ENDIF

886

XN=XNP1

887

ENDDO

888

ERROR=9

889

GOTO 710

890

891

700 CONTINUE

892

C For good measure, we iterate one last time

893

CALL FUNCT(PVECSQ(1),XN,.FALSE.,ERR, FVAL)

894

IF (ERR.NE.0) THEN

895

ERROR=ERR

896

GOTO 710

897

ENDIF

898

CALL FUNCT(PVECSQ(1),XN,.TRUE.,ERR, DVAL)

899

IF (ERR.NE.0) THEN

900

ERROR=ERR

901

GOTO 710

902

ENDIF

903

904

XSCALE=XN-(FVAL/DVAL)

905

906

710 CONTINUE

907

908

END

909

910

SUBROUTINE FUNCT(PVECSQ,X,DERIVATIVE,ERROR,RES)

911

IMPLICIT NONE

912

913

C CONSTANTS

914

915

INTEGER NEXTERNAL

916

PARAMETER (NEXTERNAL=4)

917

INTEGER NINITIAL

918

PARAMETER (NINITIAL=2)

919

REAL*16 ZERO

920

PARAMETER (ZERO=0.0E+00_16)

921

REAL*16 MP__ZERO

922

PARAMETER (MP__ZERO=ZERO)

923

REAL*16 ONE

924

PARAMETER (ONE=1.0E+00_16)

925

REAL*16 TWO

926

PARAMETER (TWO=2.0E+00_16)

927

928

C ARGUMENTS

929

930

REAL*16 PVECSQ(NEXTERNAL),X,RES

931

INTEGER ERROR

932

LOGICAL DERIVATIVE

933

934

C LOCAL VARIABLES

935

936

INTEGER I,J

937

REAL*16 BUFF,FACTOR

938

REAL*16 MASSES(NEXTERNAL)

939

940

C GLOBAL VARIABLES

941

942

INCLUDE 'mp_coupl.inc'

943

944

C ----------

945

C BEGIN CODE

946

C ----------

947

948

MASSES(1)=MP__ZERO

949

MASSES(2)=MP__ZERO

950

MASSES(3)=MP__MDL_MT

951

MASSES(4)=MP__MDL_MT

952

953

ERROR=0

954

RES=ZERO

955

BUFF=ZERO

956

957

DO I=1,NEXTERNAL

958

IF (I.LE.NINITIAL) THEN

959

FACTOR=-ONE

960

ELSE

961

FACTOR=ONE

962

ENDIF

963

BUFF=MASSES(I)**2+PVECSQ(I)*X**2

964

IF (BUFF.LT.ZERO) THEN

965

RES=ZERO

966

ERROR = 1

967

GOTO 800

968

ENDIF

969

IF (DERIVATIVE) THEN

970

RES=RES + FACTOR*((X*PVECSQ(I))/SQRT(BUFF))

971

ELSE

972

RES=RES + FACTOR*SQRT(BUFF)

973

ENDIF

974

ENDDO

975

976

800 CONTINUE

977

978

END

979

Older »

~maddevelopers/mg5amcnlo/2.0.2_onshellsubtract_newresh