~maddevelopers/mg5amcnlo/2.9.4

C Or, because this is very slow, we keep the orignal value. The accuracy returned is -1 and tells the MC that he should not trust the evaluation for checks.

336

CTMODEINIT=CTMODEINIT_BU

337

ENDIF

338

339

IF(DONEHELDOUBLECHECK.AND.(.NOT.HELDOUBLECHECKED)) THEN

340

HELDOUBLECHECKED=.TRUE.

341

DONEHELDOUBLECHECK=.FALSE.

342

ENDIF

343

344

CHECKPHASE=(NTRY.LE.CHECKCYCLE).AND.(((.NOT.FOUNDLOOPFILTER).AND.USELOOPFILTER).OR.(.NOT.FOUNDHELFILTER))

345

346

IF ((.NOT. CHECKPHASE).AND.(.NOT.FOUNDHELFILTER)) THEN

347

OPEN(1, FILE="HelFilter.dat", err=110, status='NEW', action='WRITE')

348

DO I=1,NCOMB

349

WRITE(1,*) GOODHEL(I)

350

ENDDO

351

110 CONTINUE

352

CLOSE(1)

353

FOUNDHELFILTER=.TRUE.

354

ENDIF

355

356

IF ((.NOT. CHECKPHASE).AND.(.NOT.FOUNDLOOPFILTER).AND.USELOOPFILTER) THEN

357

OPEN(1, FILE="LoopFilter.dat", err=111, status='NEW', action='WRITE')

358

DO J=1,NLOOPGROUPS

359

WRITE(1,*) GOODAMP(J)

360

ENDDO

361

111 CONTINUE

362

CLOSE(1)

363

FOUNDLOOPFILTER=.TRUE.

364

ENDIF

365

366

IF (BYPASS_CHECK) THEN

367

OLD_CHECKPHASE = CHECKPHASE

368

OLD_HELDOUBLECHECKED = HELDOUBLECHECKED

369

CHECKPHASE = .FALSE.

370

HELDOUBLECHECKED = .TRUE.

371

DO I=1,NCOMB

372

OLD_GOODHEL(I)=GOODHEL(I)

373

GOODHEL(I)=1

374

ENDDO

375

DO I=1,NLOOPGROUPS

376

OLD_GOODAMP(I)=GOODAMP(I)

377

GOODAMP(I)=.TRUE.

378

ENDDO

379

ENDIF

380

381

IF(CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED)) THEN

382

HELPICKED=1

383

CTMODE=CTMODEINIT

384

ELSE

385

IF (USERHEL.ne.-1) then

386

IF(GOODHEL(USERHEL).eq.-HELOFFSET) THEN

387

ANS(1)=0.0d0

388

ANS(2)=0.0d0

389

ANS(3)=0.0d0

390

goto 9999

391

ENDIF

392

ENDIF

393

HELPICKED=USERHEL

394

IF (CTMODERUN.NE.-1) THEN

395

CTMODE=CTMODERUN

396

ELSE

397

CTMODE=1

398

ENDIF

399

ENDIF

400

401

DO I=1,NEXTERNAL

402

DO J=0,3

403

PS(J,I)=P_USER(J,I)

404

ENDDO

405

ENDDO

406

407

IF (ImprovePSPoint.ge.0) THEN

408

C Make the input PS more precise (exact onshell and energy-momentum conservation)

409

CALL IMPROVE_PS_POINT_PRECISION(PS)

410

ENDIF

411

412

DO I=1,NEXTERNAL

413

DO J=0,3

414

P(J,I)=PS(J,I)

415

ENDDO

416

ENDDO

417

418

%(set_reference)s

419

420

DO K=1, 3

421

BUFFR(K)=0.0d0

422

DO I=1,NCTAMPS

423

AMPL(K,I)=(0.0d0,0.0d0)

424

ENDDO

425

ENDDO

426

427

200 CONTINUE

428

429

IF (.NOT.MP_PS_SET.AND.(CTMODE.EQ.0.OR.CTMODE.GE.4)) THEN

430

CALL SET_MP_PS(P_USER)

431

MP_PS_SET = .TRUE.

432

ENDIF

433

434

LSCALE=DSQRT(ABS((P(0,1)+P(0,2))**2-(P(1,1)+P(1,2))**2-(P(2,1)+P(2,2))**2-(P(3,1)+P(3,2))**2))

435

436

DO I=1,NLOOPS

437

DO J=0,LOOPMAXCOEFS-1

438

LOOPCOEFS(J,I)=(0.0d0,0.0d0)

439

ENDDO

440

ENDDO

441

442

DO K=1,3

443

ANS(K)=0.0d0

444

ENDDO

445

446

C Check if we directly go to multiple precision

447

IF (CTMODE.GE.4) THEN

448

IF (.NOT.MP_DONE) THEN

449

CALL MP_COMPUTE_LOOP_COEFS(MP_P,BUFFR_BIS)

450

C It should be safe to directly set MP_DONE to true already here. But maybe I overlooked something.

451

MP_DONE=.TRUE.

452

ENDIF

453

C Even if MP_DONE is .TRUE. we should anyway skip the

454

C double precision evaluation as it as already been

455

c computed in quadruple precision.

456

goto 300

457

ENDIF

458

459

DO H=1,NCOMB

460

IF ((HELPICKED.EQ.H).OR.((HELPICKED.EQ.-1).AND.(CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED).OR.(GOODHEL(H).GT.-HELOFFSET.AND.GOODHEL(H).NE.0)))) THEN

461

DO I=1,NEXTERNAL

462

NHEL(I)=HELC(I,H)

463

ENDDO

464

%(born_ct_helas_calls)s

465

IF (.NOT.CHECKPHASE.AND.HELDOUBLECHECKED.AND.HELPICKED.EQ.-1) THEN

466

DUMMY=GOODHEL(H)

467

ELSE

468

DUMMY=1

469

ENDIF

470

DO I=1,%(nctamps_or_nloopamps)s

471

DO J=1,%(nbornamps_or_nloopamps)s

472

CFTOT=DCMPLX(CF_N(I,J)/DBLE(ABS(CF_D(I,J))),0.0d0)

473

IF(CF_D(I,J).LT.0) CFTOT=CFTOT*IMAG1

474

%(squaring)s

475

ENDDO

476

ENDDO

477

%(coef_construction)s

478

ENDIF

479

ENDDO

480

481

%(coef_merging)s

482

483

BUFFR_BIS(1)=ANS(1)

484

BUFFR_BIS(2)=ANS(2)

485

BUFFR_BIS(3)=ANS(3)

486

300 CONTINUE

487

ANS(1)=BUFFR_BIS(1)

488

ANS(2)=BUFFR_BIS(2)

489

ANS(3)=BUFFR_BIS(3)

490

491

IF(SKIPLOOPEVAL) THEN

492

GOTO 1226

493

ENDIF

494

495

%(loop_CT_calls)s

496

497

%(actualize_ans)s

498

499

1226 CONTINUE

500

501

IF (CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED)) THEN

502

IF((USERHEL.EQ.-1).OR.(USERHEL.EQ.HELPICKED)) THEN

503

C TO KEEP TRACK OF THE FINAL ANSWER TO BE RETURNED DURING CHECK PHASE

504

BUFFR(1)=BUFFR(1)+ANS(1)

505

BUFFR(2)=BUFFR(2)+ANS(2)

506

BUFFR(3)=BUFFR(3)+ANS(3)

507

ENDIF

508

C SAVE RESULT OF EACH INDEPENDENT HELICITY FOR COMPARISON DURING THE HELICITY FILTER SETUP

509

HELSAVED(1,HELPICKED)=ANS(1)

510

HELSAVED(2,HELPICKED)=ANS(2)

511

HELSAVED(3,HELPICKED)=ANS(3)

512

513

IF (CHECKPHASE) THEN

514

C SET THE HELICITY FILTER

515

IF(.NOT.FOUNDHELFILTER) THEN

516

HEL_INCONSISTENT=.FALSE.

517

IF(ISZERO(ABS(ANS(1))+ABS(ANS(2))+ABS(ANS(3)),REF/DBLE(NCOMB),-1)) THEN

518

IF(NTRY.EQ.1) THEN

519

GOODHEL(HELPICKED)=-HELOFFSET

520

ELSEIF(GOODHEL(HELPICKED).NE.-HELOFFSET) THEN

521

WRITE(*,*) '##W02A WARNING Inconsistent helicity ',HELPICKED

522

IF(HELINITSTARTOVER) THEN

523

WRITE(*,*) '##I01 INFO Initialization starting over because of inconsistency in the helicity filter setup.'

524

NTRY=0

525

ELSE

526

HEL_INCONSISTENT=.TRUE.

527

ENDIF

528

ENDIF

529

ELSE

530

DO H=1,HELPICKED-1

531

IF(GOODHEL(H).GT.-HELOFFSET) THEN

532

C Be looser for helicity check, bring a factor 100

533

DUMMY=ISSAME(ANS,HELSAVED(1,H),REF,.FALSE.)

534

IF(DUMMY.NE.0) THEN

535

IF(NTRY.EQ.1) THEN

536

C Set the matching helicity to be contributing once more

537

GOODHEL(H)=GOODHEL(H)+DUMMY

538

C Use an offset to clearly show it is linked to an other one and to avoid overlap

539

GOODHEL(HELPICKED)=-H-HELOFFSET

540

C Make sure we have paired this hel config to the same one last PS point

541

ELSEIF(GOODHEL(HELPICKED).NE.(-H-HELOFFSET)) THEN

542

WRITE(*,*) '##W02B WARNING Inconsistent helicity ',HELPICKED

543

IF(HELINITSTARTOVER) THEN

544

WRITE(*,*) '##I01 INFO Initialization starting over because of inconsistency in the helicity filter setup.'

545

NTRY=0

546

ELSE

547

HEL_INCONSISTENT=.TRUE.

548

ENDIF

549

ENDIF

550

ENDIF

551

ENDIF

552

ENDDO

553

ENDIF

554

IF(HEL_INCONSISTENT) THEN

555

C This helicity has unstable filter so we will always compute it by itself.

556

C We therefore also need to remove it from the multiplicative factor of the corresponding helicity.

557

IF(GOODHEL(HELPICKED).LT.-HELOFFSET) THEN

558

GOODHEL(-GOODHEL(HELPICKED)-HELOFFSET)=GOODHEL(-GOODHEL(HELPICKED)-HELOFFSET)-1

559

ENDIF

560

C If several helicities were matched to that one, we need to chose another one as reference and redirect the others to this new one

561

C Of course if it is one, then we do not need to do anything (because with HELINITSTARTOVER=.FALSE. we only support exactly identical Hels.)

562

IF(GOODHEL(HELPICKED).GT.-HELOFFSET.AND.GOODHEL(HELPICKED).NE.1) THEN

563

NEWHELREF=-1

564

DO H=1,NCOMB

565

IF (GOODHEL(H).EQ.(-HELOFFSET-HELPICKED)) THEN

566

IF (NEWHELREF.EQ.-1) THEN

567

NEWHELREF=H

568

GOODHEL(H)=GOODHEL(HELPICKED)-1

569

ELSE

570

GOODHEL(H)=-NEWHELREF-HELOFFSET

571

ENDIF

572

ENDIF

573

ENDDO

574

ENDIF

575

C In all cases, from now on this helicity will be computed independantly of the others.

576

C In particular, it is the only thing to do if the helicity was flagged not contributing.

577

GOODHEL(HELPICKED)=1

578

ENDIF

579

ENDIF

580

581

C SET THE LOOP FILTER

582

IF(.NOT.FOUNDLOOPFILTER.AND.USELOOPFILTER) THEN

583

DO I=1,NLOOPGROUPS

584

IF(.NOT.ISZERO(ABS(LOOPRES(1,I))+ABS(LOOPRES(2,I))+ABS(LOOPRES(3,I)),(REF*1.0d-4),I)) THEN

585

IF(NTRY.EQ.1) THEN

586

GOODAMP(I)=.TRUE.

587

LOOPFILTERBUFF(I)=.TRUE.

588

ELSEIF(.NOT.LOOPFILTERBUFF(I)) THEN

589

WRITE(*,*) '##W02 WARNING Inconsistent loop amp ',I,'.'

590

IF(LOOPINITSTARTOVER) THEN

591

WRITE(*,*) '##I01 INFO Initialization starting over because of inconsistency in the loop filter setup.'

592

NTRY=0

593

ELSE

594

GOODAMP(I)=.TRUE.

595

ENDIF

596

ENDIF

597

ENDIF

598

ENDDO

599

ENDIF

600

ELSEIF (.NOT.HELDOUBLECHECKED.AND.NTRY.NE.0)THEN

601

C DOUBLE CHECK THE HELICITY FILTER

602

IF (GOODHEL(HELPICKED).EQ.-HELOFFSET) THEN

603

IF (.NOT.ISZERO(ABS(ANS(1))+ABS(ANS(2))+ABS(ANS(3)),REF/DBLE(NCOMB),-1)) THEN

604

write(*,*) '##W15 Helicity filter could not be successfully double checked.'

605

write(*,*) 'One reason for this is that you might have changed sensible parameters which affected what are the zero helicity configurations.'

606

write(*,*) 'MadLoop will try to reset the Helicity filter with the next PS points it receives.'

607

NTRY=0

608

OPEN(29,FILE='HelFilter.dat',err=348)

609

348 CONTINUE

610

CLOSE(29,STATUS='delete')

611

ENDIF

612

ENDIF

613

IF (GOODHEL(HELPICKED).LT.-HELOFFSET.AND.NTRY.NE.0) THEN

614

IF(ISSAME(ANS,HELSAVED(1,ABS(GOODHEL(HELPICKED)+HELOFFSET)),REF,.TRUE.).EQ.0) THEN

615

write(*,*) '##W15 Helicity filter could not be successfully double checked.'

616

write(*,*) 'One reason for this is that you might have changed sensible parameters which affected the helicity dependance relations.'

617

write(*,*) 'MadLoop will try to reset the Helicity filter with the next PS points it receives.'

618

NTRY=0

619

OPEN(30,FILE='HelFilter.dat',err=349)

620

349 CONTINUE

621

CLOSE(30,STATUS='delete')

622

ENDIF

623

ENDIF

624

C SET HELDOUBLECHECKED TO .TRUE. WHEN DONE

625

C even if it failed we do not want to redo the check afterwards if HELINITSTARTOVER=.FALSE.

626

IF (HELPICKED.EQ.NCOMB.AND.(NTRY.NE.0.OR..NOT.HELINITSTARTOVER)) THEN

627

DONEHELDOUBLECHECK=.TRUE.

628

ENDIF

629

ENDIF

630

631

C GOTO NEXT HELICITY OR FINISH

632

IF(HELPICKED.NE.NCOMB) THEN

633

HELPICKED=HELPICKED+1

634

MP_DONE=.FALSE.

635

goto 200

636

ELSE

637

C Useful printout

638

c do I=1,NCOMB

639

c write(*,*) 'HELSAVED(1,',I,')=',HELSAVED(1,I)

640

c write(*,*) 'HELSAVED(2,',I,')=',HELSAVED(2,I)

641

c write(*,*) 'HELSAVED(3,',I,')=',HELSAVED(3,I)

642

c write(*,*) ' GOODHEL(',I,')=',GOODHEL(I)

643

c ENDDO

644

ANS(1)=BUFFR(1)

645

ANS(2)=BUFFR(2)

646

ANS(3)=BUFFR(3)

647

IF(NTRY.EQ.0) THEN

648

NATTEMPTS=NATTEMPTS+1

649

IF(NATTEMPTS.EQ.MAXATTEMPTS) THEN

650

WRITE(*,*) '##E01 ERROR Could not initialize the filters in ',MAXATTEMPTS,' trials'

651

STOP 1

652

ENDIF

653

ENDIF

654

ENDIF

655

656

ENDIF

657

658

DO K=1,3

659

ANS(K)=ANS(K)/DBLE(IDEN)

660

IF (USERHEL.NE.-1) THEN

661

ANS(K)=ANS(K)*HELAVGFACTOR

662

ENDIF

663

ENDDO

664

665

IF(.NOT.CHECKPHASE.AND.HELDOUBLECHECKED.AND.(CTMODERUN.EQ.-1)) THEN

666

STAB_INDEX=STAB_INDEX+1

667

IF(DOING_QP_EVALS) THEN

668

QP_RES(1,STAB_INDEX)=ANS(1)

669

QP_RES(2,STAB_INDEX)=ANS(2)

670

QP_RES(3,STAB_INDEX)=ANS(3)

671

ELSE

672

DP_RES(1,STAB_INDEX)=ANS(1)

673

DP_RES(2,STAB_INDEX)=ANS(2)

674

DP_RES(3,STAB_INDEX)=ANS(3)

675

ENDIF

676

677

IF(DOING_QP_EVALS) THEN

678

BASIC_CT_MODE=4

679

ELSE

680

BASIC_CT_MODE=1

681

ENDIF

682

683

C BEGINNING OF THE DEFINITIONS OF THE DIFFERENT EVALUATION METHODS

684

685

IF(.NOT.EVAL_DONE(2)) THEN

686

EVAL_DONE(2)=.TRUE.

687

CTMODE=BASIC_CT_MODE+1

688

goto 300

689

ENDIF

690

691

CTMODE=BASIC_CT_MODE

692

693

IF(.NOT.EVAL_DONE(3).AND. ((DOING_QP_EVALS.AND.NRotations_QP.GE.1).OR.((.NOT.DOING_QP_EVALS).AND.NRotations_DP.GE.1)) ) THEN

694

EVAL_DONE(3)=.TRUE.

695

CALL ROTATE_PS(PS,P,1)

696

IF (DOING_QP_EVALS) CALL MP_ROTATE_PS(MP_PS,MP_P,1)

697

goto 200

698

ENDIF

699

700

IF(.NOT.EVAL_DONE(4).AND. ((DOING_QP_EVALS.AND.NRotations_QP.GE.2).OR.((.NOT.DOING_QP_EVALS).AND.NRotations_DP.GE.2)) ) THEN

701

EVAL_DONE(4)=.TRUE.

702

CALL ROTATE_PS(PS,P,2)

703

IF (DOING_QP_EVALS) CALL MP_ROTATE_PS(MP_PS,MP_P,2)

704

goto 200

705

ENDIF

706

707

CALL ROTATE_PS(PS,P,0)

708

IF (DOING_QP_EVALS) CALL MP_ROTATE_PS(MP_PS,MP_P,0)

709

710

C END OF THE DEFINITIONS OF THE DIFFERENT EVALUATION METHODS

711

712

IF(DOING_QP_EVALS) THEN

713

CALL COMPUTE_ACCURACY(QP_RES,N_QP_EVAL,ACC,ANS)

714

ACCURACY=ACC

715

RET_CODE_H=3

716

IF(ACC.GE.MLSTABTHRES) THEN

717

RET_CODE_H=4

718

NEPS=NEPS+1

719

CALL COMPUTE_ACCURACY(DP_RES,N_DP_EVAL,TEMP1,TEMP)

720

WRITE(*,*) '##W03 WARNING An exceptional PS point was', ' detected.'

721

WRITE(*,*) '(DP,QP) accuracies : (',TEMP1,',',ACC,')'

722

WRITE(*,*) 'Best estimate (fin,1eps,2eps) :',(ANS(I),I=1,3)

723

IF(NEPS.LE.10) THEN

724

WRITE(*,*) 'Double precision evaluations :',(DP_RES(1,I),I=1,N_DP_EVAL)

725

WRITE(*,*) 'Quad precision evaluations :',(QP_RES(1,I),I=1,N_QP_EVAL)

726

WRITE(*,*) 'PS point specification :'

727

WRITE(*,*) 'Renormalization scale MU_R=',MU_R

728

DO I=1,NEXTERNAL

729

WRITE (*,'(i2,1x,4e27.17)') i, P(0,i),P(1,i),P(2,i),P(3,i)

730

ENDDO

731

ENDIF

732

IF(NEPS.EQ.10) THEN

733

WRITE(*,*) 'Further output of the details of these unstable PS points will now be suppressed.'

734

ENDIF

735

ENDIF

736

ELSE

737

CALL COMPUTE_ACCURACY(DP_RES,N_DP_EVAL,ACC,ANS)

738

IF(ACC.GE.MLSTABTHRES) THEN

739

DOING_QP_EVALS=.TRUE.

740

EVAL_DONE(1)=.TRUE.

741

DO I=2,MAXSTABILITYLENGTH

742

EVAL_DONE(I)=.FALSE.

743

ENDDO

744

STAB_INDEX=0

745

CTMODE=4

746

goto 200

747

ELSE

748

RET_CODE_H=2

749

ACCURACY=ACC

750

ENDIF

751

ENDIF

752

ELSE

753

RET_CODE_H=1

754

ACCURACY=-1.0d0

755

ENDIF

756

757

9999 CONTINUE

758

759

C Finalize the return code

760

IF (MP_DONE_ONCE) THEN

761

RET_CODE_T=2

762

ELSE

763

RET_CODE_T=1

764

ENDIF

765

IF(CHECKPHASE.OR..NOT.HELDOUBLECHECKED) THEN

766

RET_CODE_H=1

767

RET_CODE_T=RET_CODE_T+2

768

ACCURACY=-1.0d0

769

ENDIF

770

RET_CODE_U=0

771

772

C Reinitialize the default threshold if it was specified by the user

773

IF (USER_STAB_PREC.GT.0.0d0) THEN

774

MLSTABTHRES=MLSTABTHRES_BU

775

CTMODEINIT=CTMODEINIT_BU

776

ENDIF

777

778

C Reinitialize the check phase logicals and the filters if check bypassed

779

IF (BYPASS_CHECK) THEN

780

CHECKPHASE = OLD_CHECKPHASE

781

HELDOUBLECHECKED = OLD_HELDOUBLECHECKED

782

DO I=1,NCOMB

783

GOODHEL(I)=OLD_GOODHEL(I)

784

ENDDO

785

DO I=1,NLOOPGROUPS

786

GOODAMP(I)=OLD_GOODAMP(I)

787

ENDDO

788

ENDIF

789

END

790

791

792

logical function IsZero(toTest, reference_value, loop)

793

IMPLICIT NONE

794

795

C CONSTANTS

796

797

INTEGER NLOOPGROUPS

798

PARAMETER (NLOOPGROUPS=%(nloop_groups)d)

799

800

C ARGUMENTS

801

802

%(real_dp_format)s toTest, reference_value

803

integer loop

804

805

C GLOBAL

806

807

INCLUDE 'MadLoopParams.inc'

808

%(complex_dp_format)s LOOPRES(3,NLOOPGROUPS)

809

LOGICAL S(NLOOPGROUPS)

810

common/LOOPRES/LOOPRES,S

811

C ----------

812

C BEGIN CODE

813

C ----------

814

IF(abs(reference_value).eq.0.0d0) then

815

IsZero=.FALSE.

816

write(*,*) '##E02 ERRROR Reference value for comparison is zero.'

817

STOP 1

818

else

819

IsZero=((abs(toTest)/abs(reference_value)).lt.ZEROTHRES)

820

endif

821

822

IF(loop.NE.-1) THEN

823

IF((.NOT.ISZERO).AND.(.NOT.S(loop))) THEN

824

write(*,*) '##W01 WARNING Contribution ',loop,' is detected as contributing with CR=',(abs(toTest)/abs(reference_value)),' but is unstable.'

825

ENDIF

826

ENDIF

827

828

end

829

830

integer function ISSAME(resA,resB,ref,useMax)

831

IMPLICIT NONE

832

C This function compares the result from two different helicity configuration A and B

833

C It returns 0 if they are not related and (+/-wgt) if A=(+/-wgt)*B.

834

C For now, the only wgt implemented is the integer 1 or -1.

835

C If useMax is .TRUE., it uses all implemented weights no matter what is HELINITSTARTOVER

836

837

C CONSTANTS

838

839

integer MAX_WGT_TO_TRY

840

parameter (MAX_WGT_TO_TRY=2)

841

842

C ARGUMENTS

843

844

%(real_dp_format)s resA(3), resB(3)

845

%(real_dp_format)s ref

846

logical useMax

847

848

C LOCAL VARIABLES

849

850

LOGICAL ISZERO

851

integer I,J

852

integer N_WGT_TO_TRY

853

integer WGT_TO_TRY(MAX_WGT_TO_TRY)

854

data WGT_TO_TRY/1,-1/

855

856

C INCLUDES

857

858

include 'MadLoopParams.inc'

859

C ----------

860

C BEGIN CODE

861

C ----------

862

ISSAME=0

863

864

C If the helicity can be constructed progressively while allowing inconsistency, then we only allow for weight one comparisons.

865

IF (.NOT.HELINITSTARTOVER.AND..NOT.USEMAX) THEN

866

N_WGT_TO_TRY=1

867

ELSE

868

N_WGT_TO_TRY=MAX_WGT_TO_TRY

869

ENDIF

870

871

DO I=1,N_WGT_TO_TRY

872

DO J=1,3

873

IF (IsZero(ABS(resB(J)),ref,-1)) THEN

874

IF(.NOT.IsZero(ABS(resB(J))+ABS(resA(J)),ref,-1)) THEN

875

GOTO 1231

876

ENDIF

877

C Be loser for helicity comparison, so bring a factor 100

878

ELSEIF(.NOT.IsZero((resA(J)/resB(J))-DBLE(WGT_TO_TRY(I)),ref*100.0d0,-1)) THEN

879

GOTO 1231

880

ENDIF

881

ENDDO

882

ISSAME = WGT_TO_TRY(I)

883

RETURN

884

1231 CONTINUE

885

ENDDO

886

END

887

888

SUBROUTINE compute_accuracy(fulllist, length, acc, estimate)

889

implicit none

890

891

C PARAMETERS

892

893

integer maxstabilitylength

894

common/stability_tests/maxstabilitylength

895

896

C ARGUMENTS

897

898

real*8 fulllist(3,maxstabilitylength)

899

integer length

900

real*8 acc, estimate(3)

901

902

C LOCAL VARIABLES

903

904

logical mask(maxstabilitylength)

905

logical mask3(3)

906

data mask3/.TRUE.,.TRUE.,.TRUE./

907

integer i,j

908

real*8 avg

909

real*8 diff

910

real*8 accuracies(3)

911

real*8 list(maxstabilitylength)

912

913

C ----------

914

C BEGIN CODE

915

C ----------

916

do i=1,length

917

mask(i)=.TRUE.

918

enddo

919

do i=length+1,maxstabilitylength

920

mask(i)=.FALSE.

921

C For some architectures, it is necessary to initialize all the elements of fulllist(i,j)

922

C Beware that if the length provided is incorrect, then this can corrup the fulllist given in argument.

923

do j=1,3

924

fulllist(j,i)=0.0d0

925

enddo

926

enddo

927

928

do i=1,3

929

do j=1,maxstabilitylength

930

list(j)=fulllist(i,j)

931

enddo

932

diff=maxval(list,1,mask)-minval(list,1,mask)

933

avg=(maxval(list,1,mask)+minval(list,1,mask))/2.0d0

934

estimate(i)=avg

935

if (avg.eq.0.0d0) then

936

accuracies(i)=diff

937

else

938

accuracies(i)=diff/abs(avg)

939

endif

940

enddo

941

942

C The technique below is too sensitive, typically to

943

C unstablities in very small poles

944

C ACC=MAXVAL(ACCURACIES,1,MASK3)

945

C The following is used instead

946

ACC = 0.0d0

947

AVG = 0.0d0

948

DO I=1,3

949

ACC = ACC + ACCURACIES(I)*ABS(ESTIMATE(I))

950

AVG = AVG + ESTIMATE(I)

951

ENDDO

952

ACC = ACC / ( ABS(AVG) / 3.0d0)

953

954

end

955

956

SUBROUTINE SET_N_EVALS(N_DP_EVALS,N_QP_EVALS)

957

958

IMPLICIT NONE

959

INTEGER N_DP_EVALS, N_QP_EVALS

960

961

include 'MadLoopParams.inc'

962

963

IF(CTMODERUN.LE.-1) THEN

964

N_DP_EVALS=2+NRotations_DP

965

N_QP_EVALS=2+NRotations_QP

966

ELSE

967

N_DP_EVALS=1

968

N_QP_EVALS=1

969

ENDIF

970

971

IF(N_DP_EVALS.GT.20.OR.N_QP_EVALS.GT.20) THEN

972

WRITE(*,*) 'ERROR:: Increase hardcoded maxstabilitylength.'

973

stop 1

974

ENDIF

975

976

END

977

978

C THIS SUBROUTINE SIMPLY SET THE GLOBAL PS CONFIGURATION GLOBAL VARIABLES FROM A GIVEN VARIABLE IN DOUBLE PRECISION

979

SUBROUTINE SET_MP_PS(P)

980

981

INTEGER NEXTERNAL

982

PARAMETER (NEXTERNAL=%(nexternal)d)

983

%(real_mp_format)s MP_PS(0:3,NEXTERNAL),MP_P(0:3,NEXTERNAL)

984

common/MP_PSPOINT/MP_PS,MP_P

985

%(real_dp_format)s P(0:3,NEXTERNAL)

986

987

DO I=1,NEXTERNAL

988

DO J=0,3

989

MP_PS(J,I)=P(J,I)

990

ENDDO

991

ENDDO

992

CALL MP_IMPROVE_PS_POINT_PRECISION(MP_PS)

993

DO I=1,NEXTERNAL

994

DO J=0,3

995

MP_P(J,I)=MP_PS(J,I)

996

ENDDO

997

ENDDO

998

999

END

1000

1001

SUBROUTINE FORCE_STABILITY_CHECK(ONOFF)

1002

1003

C This function can be called by the MadLoop user so as to always have stability

1004

C checked, even during initialisation, when calling the *_thres routines.

1005

1006

LOGICAL ONOFF

1007

1008

LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY

1009

DATA BYPASS_CHECK, ALWAYS_TEST_STABILITY /.FALSE.,.FALSE./

1010

COMMON/BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY

1011

1012

ALWAYS_TEST_STABILITY = ONOFF

1013

1014

END

1015

1016

SUBROUTINE SLOOPMATRIXHEL_THRES(P,HEL,ANS,PREC_ASKED,PREC_FOUND,RET_CODE)

1017

IMPLICIT NONE

1018

1019

C CONSTANTS

1020

1021

INTEGER NEXTERNAL

1022

PARAMETER (NEXTERNAL=%(nexternal)d)

1023

1024

C ARGUMENTS

1025

1026

%(real_dp_format)s P(0:3,NEXTERNAL)

1027

%(real_dp_format)s ANS(3)

1028

INTEGER HEL,RET_CODE

1029

%(real_dp_format)s PREC_ASKED,PREC_FOUND

1030

1031

C GLOBAL VARIABLES

1032

1033

%(real_dp_format)s USER_STAB_PREC

1034

COMMON/USER_STAB_PREC/USER_STAB_PREC

1035

1036

INTEGER H,T,U

1037

%(real_dp_format)s ACCURACY

1038

common/ACC/ACCURACY,H,T,U

1039

1040

LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY

1041

COMMON/BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY

1042

1043

C ----------

1044

C BEGIN CODE

1045

C ----------

1046

USER_STAB_PREC = PREC_ASKED

1047

1048

CALL SLOOPMATRIXHEL(P,HEL,ANS)

1049

IF(ALWAYS_TEST_STABILITY.AND.(H.eq.1.OR.ACCURACY.lt.0.0d0)) THEN

1050

BYPASS_CHECK = .TRUE.

1051

CALL SLOOPMATRIXHEL(P,HEL,ANS)

1052

BYPASS_CHECK = .FALSE.

1053

C Make sure we correctly return an initialization-type T code

1054

IF (T.eq.2) T=4

1055

IF (T.eq.1) T=3

1056

ENDIF

1057

1058

C Reset it to default value not to affect next runs

1059

USER_STAB_PREC = -1.0d0

1060

PREC_FOUND=ACCURACY

1061

RET_CODE=100*H+10*T+U

1062

1063

END

1064

1065

SUBROUTINE SLOOPMATRIX_THRES(P,ANS,PREC_ASKED,PREC_FOUND,RET_CODE)

1066

1067

C Inputs are:

1068

C P(0:3, Nexternal) double :: Kinematic configuration (E,px,py,pz)

1069

C PEC_ASKED double :: Target relative accuracy, -1 for default

1070

1071

C Outputs are:

1072

C ANS(3) double :: Result (finite, single pole, double pole)

1073

C PREC_FOUND double :: Relative accuracy estimated for the result

1074

C Returns -1 if no stab test could be performed.

1075

C RET_CODE integer :: Return code. See below for details

1076

1077

C Return code conventions: RET_CODE = H*100 + T*10 + U

1078

1079

C H == 1

1080

C Stability unknown.

1081

C H == 2

1082

C Stable PS (SPS) point.

1083

C No stability rescue was necessary.

1084

C H == 3

1085

C Unstable PS (UPS) point.

1086

C Stability rescue necessary, and successful.

1087

C H == 4

1088

C Exceptional PS (EPS) point.

1089

C Stability rescue attempted, but unsuccessful.

1090

1091

C T == 1

1092

C Default computation (double prec.) was performed.

1093

C T == 2

1094

C Quadruple precision was used for this PS point.

1095

C T == 3

1096

C MadLoop in initialization phase. Only double precision used.

1097

C T == 4

1098

C MadLoop in initialization phase. Quadruple precision used.

1099

1100

C U is a number left for future use (always set to 0 for now).

1101

C example: TIR vs OPP usage.

1102

1103

IMPLICIT NONE

1104

1105

C CONSTANTS

1106

1107

INTEGER NEXTERNAL

1108

PARAMETER (NEXTERNAL=%(nexternal)d)

1109

1110

C ARGUMENTS

1111

1112

%(real_dp_format)s P(0:3,NEXTERNAL)

1113

%(real_dp_format)s ANS(3)

1114

%(real_dp_format)s PREC_ASKED,PREC_FOUND

1115

INTEGER RET_CODE

1116

1117

C GLOBAL VARIABLES

1118

1119

%(real_dp_format)s USER_STAB_PREC

1120

COMMON/USER_STAB_PREC/USER_STAB_PREC

1121

1122

INTEGER H,T,U

1123

%(real_dp_format)s ACCURACY

1124

common/ACC/ACCURACY,H,T,U

1125

1126

LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY

1127

COMMON/BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY

1128

1129

C ----------

1130

C BEGIN CODE

1131

C ----------

1132

USER_STAB_PREC = PREC_ASKED

1133

1134

CALL SLOOPMATRIX(P,ANS)

1135

IF(ALWAYS_TEST_STABILITY.AND.(H.eq.1.OR.ACCURACY.lt.0.0d0)) THEN

1136

BYPASS_CHECK = .TRUE.

1137

CALL SLOOPMATRIX(P,ANS)

1138

BYPASS_CHECK = .FALSE.

1139

C Make sure we correctly return an initialization-type T code

1140

IF (T.eq.2) T=4

1141

IF (T.eq.1) T=3

1142

ENDIF

1143

1144

C Reset it to default value not to affect next runs

1145

USER_STAB_PREC = -1.0d0

1146

PREC_FOUND=ACCURACY

1147

RET_CODE=100*H+10*T+U

1148

1149

END

Older »