← Back to branch summary

~maddevelopers/mg5amcnlo/WWW5_caching

« back to all changes in this revision

Viewing changes to users/mardelcourt/PROC_407857/PROC_407857/Source/CERNLIB/dlsqp2.f

  • Committer: John Doe
  • Date: 2013-03-25 20:27:02 UTC
  • Revision ID: john.doe@gmail.com-20130325202702-5sk3t1r8h33ca4p4
first clean version

Show diffs side-by-side

added added

removed removed

Lines of Context:
users/mardelcourt/PROC_407857/PROC_407857/Source/CERNLIB/dlsqp2.f
 
1
*
 
2
* $Id: dlsqp2.f,v 1.1 2009/07/30 22:46:16 madgraph Exp $
 
3
*
 
4
* $Log: dlsqp2.f,v $
 
5
* Revision 1.1  2009/07/30 22:46:16  madgraph
 
6
* JA: Implemented CKKW-style matching with Pythia pT-ordered showers
 
7
*
 
8
* Revision 1.1.1.1  1996/04/01 15:02:24  mclareni
 
9
* Mathlib gen
 
10
*
 
11
*
 
12
      SUBROUTINE DLSQP2(N,X,Y,A0,A1,A2,SD,IFAIL)
 
13
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 
14
 
 
15
      DIMENSION X(*),Y(*)
 
16
 
 
17
      PARAMETER (R0 = 0)
 
18
 
 
19
      A0=0
 
20
      A1=0
 
21
      A2=0
 
22
      SD=0
 
23
      IF(N .LE. 2) THEN
 
24
       IFAIL=1
 
25
      ELSE
 
26
       FN=N
 
27
       XM=0
 
28
       DO 1 K = 1,N
 
29
       XM=XM+X(K)
 
30
    1  CONTINUE
 
31
       XM=XM/FN
 
32
       SX=0
 
33
       SXX=0
 
34
       SXXX=0
 
35
       SXXXX=0
 
36
       SY=0
 
37
       SYY=0
 
38
       SXY=0
 
39
       SXXY=0
 
40
       DO 2 K = 1,N
 
41
       XK=X(K)-XM
 
42
       YK=Y(K)
 
43
       XK2=XK**2
 
44
       SX=SX+XK
 
45
       SXX=SXX+XK2
 
46
       SXXX=SXXX+XK2*XK
 
47
       SXXXX=SXXXX+XK2**2
 
48
       SY=SY+YK
 
49
       SYY=SYY+YK**2
 
50
       SXY=SXY+XK*YK
 
51
       SXXY=SXXY+XK2*YK
 
52
    2  CONTINUE
 
53
       DET=(FN*SXXXX-SXX**2)*SXX-FN*SXXX**2
 
54
       IF(DET .GT. 0) THEN
 
55
        A2=(SXX*(FN*SXXY-SXX*SY)-FN*SXXX*SXY)/DET
 
56
        A1=(SXY-SXXX*A2)/SXX
 
57
        A0=(SY-SXX*A2)/FN
 
58
        IFAIL=0
 
59
       ELSE
 
60
        IFAIL=-1
 
61
       ENDIF
 
62
      ENDIF
 
63
      IF(IFAIL .EQ. 0 .AND. N .GT. 3)
 
64
     1   SD=SQRT(MAX(R0,SYY-A0*SY-A1*SXY-A2*SXXY)/(N-3))
 
65
      A0=A0+XM*(XM*A2-A1)
 
66
      A1=A1-2*XM*A2
 
67
      RETURN
 
68
      END
 
69