[u/mrichter/AliRoot.git] / TEvtGen / PHOTOS / phodo.F

      SUBROUTINE PHODO(IP,NCHARB,NEUDAU)
C.----------------------------------------------------------------------
C.
C.    PHOTOS:   PHOton radiation in  decays DOing of KINematics
C.
C.    Purpose:  Starting  from   the  charged  particle energy/momentum,
C.              PNEUTR, photon  energy  fraction and photon  angle  with
C.              respect  to  the axis formed by charged particle energy/
C.              momentum  vector  and PNEUTR, scale the energy/momentum,
C.              keeping the original direction of the neutral system  in
C.              the lab. frame untouched.
C.
C.    Input Parameters:   IP:      Pointer  to   decaying  particle   in
C.                                 /PHOEVT/  and   the   common   itself
C.                        NCHARB:  pointer to the charged radiating
C.                                 daughter in /PHOEVT/.
C.                        NEUDAU:  pointer to the first neutral daughter
C.    Output Parameters:  Common /PHOEVT/, with photon added.
C.
C.    Author(s):  Z. Was, B. van Eijk             Created at:  26/11/89
C.                                                Last Update: 27/05/93
C.
C.----------------------------------------------------------------------
      IMPLICIT NONE
      DOUBLE PRECISION PHOAN1,PHOAN2,ANGLE,FI1,FI3,FI4,FI5,TH1,TH3,TH4
      DOUBLE PRECISION PARNE,QNEW,QOLD,DATA
      INTEGER IP,FI3DUM,I,J,NEUDAU,FIRST,LAST
      INTEGER NCHARB
      REAL*8 EPHOTO,PMAVIR,PHOTRI
      REAL*8 GNEUT,PHORAN,CCOSTH,SSINTH,PVEC(4)
      INTEGER NMXPHO
      PARAMETER (NMXPHO=10000)
      INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
      REAL*8 PPHO,VPHO
      COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
     &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
      DOUBLE PRECISION MCHSQR,MNESQR
      REAL*8 PNEUTR
      COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
      DOUBLE PRECISION COSTHG,SINTHG
      REAL*8 XPHMAX,XPHOTO
      COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
      REAL*8 PI,TWOPI
      COMMON/PHPICO/PI,TWOPI
C--
      EPHOTO=XPHOTO*PPHO(5,IP)/2.D0
      PMAVIR=SQRT(PPHO(5,IP)*(PPHO(5,IP)-2.D0*EPHOTO))
C--
C--   Reconstruct  kinematics  of  charged particle  and  neutral system
      FI1=PHOAN1(PNEUTR(1),PNEUTR(2))
C--
C--   Choose axis along  z of  PNEUTR, calculate  angle  between x and y
C--   components  and z  and x-y plane and  perform Lorentz transform...
      TH1=PHOAN2(PNEUTR(3),SQRT(PNEUTR(1)**2+PNEUTR(2)**2))
      CALL PHORO3(-FI1,PNEUTR(1))
      CALL PHORO2(-TH1,PNEUTR(1))
C--
C--   Take  away  photon energy from charged particle and PNEUTR !  Thus
C--   the onshell charged particle  decays into virtual charged particle
C--   and photon.  The virtual charged  particle mass becomes:
C--   SQRT(PPHO(5,IP)*(PPHO(5,IP)-2*EPHOTO)).  Construct  new PNEUTR mo-
C--   mentum in the rest frame of the parent:
C--   1) Scaling parameters...
      QNEW=PHOTRI(PMAVIR,PNEUTR(5),PPHO(5,NCHARB))
      QOLD=PNEUTR(3)
      GNEUT=(QNEW**2+QOLD**2+MNESQR)/(QNEW*QOLD+SQRT((QNEW**2+MNESQR)*
     &(QOLD**2+MNESQR)))
      IF (GNEUT.LT.1.D0) THEN
        DATA=0.D0
        CALL PHOERR(4,'PHOKIN',DATA)
      ENDIF
      PARNE=GNEUT-SQRT(MAX(GNEUT**2-1.0D0,0.D0))
C--
C--   2) ...reductive boost...
      CALL PHOBO3(PARNE,PNEUTR)
C--
C--   ...calculate photon energy in the reduced system...
      NPHO=NPHO+1
      ISTPHO(NPHO)=1
      IDPHO(NPHO) =22
C--   Photon mother and daughter pointers !
      JMOPHO(1,NPHO)=IP
      JMOPHO(2,NPHO)=0
      JDAPHO(1,NPHO)=0
      JDAPHO(2,NPHO)=0
      PPHO(4,NPHO)=EPHOTO*PPHO(5,IP)/PMAVIR
C--
C--   ...and photon momenta
      CCOSTH=-COSTHG
      SSINTH=SINTHG
      TH3=PHOAN2(CCOSTH,SSINTH)
      FI3=TWOPI*PHORAN(FI3DUM)
      PPHO(1,NPHO)=PPHO(4,NPHO)*SINTHG*COS(FI3)
      PPHO(2,NPHO)=PPHO(4,NPHO)*SINTHG*SIN(FI3)
C--
C--   Minus sign because axis opposite direction of charged particle !
      PPHO(3,NPHO)=-PPHO(4,NPHO)*COSTHG
      PPHO(5,NPHO)=0.D0
C--
C--   Rotate in order to get photon along z-axis
      CALL PHORO3(-FI3,PNEUTR(1))
      CALL PHORO3(-FI3,PPHO(1,NPHO))
      CALL PHORO2(-TH3,PNEUTR(1))
      CALL PHORO2(-TH3,PPHO(1,NPHO))
      ANGLE=EPHOTO/PPHO(4,NPHO)
C--
C--   Boost to the rest frame of decaying particle
      CALL PHOBO3(ANGLE,PNEUTR(1))
      CALL PHOBO3(ANGLE,PPHO(1,NPHO))
C--
C--   Back in the parent rest frame but PNEUTR not yet oriented !
      FI4=PHOAN1(PNEUTR(1),PNEUTR(2))
      TH4=PHOAN2(PNEUTR(3),SQRT(PNEUTR(1)**2+PNEUTR(2)**2))
      CALL PHORO3(FI4,PNEUTR(1))
      CALL PHORO3(FI4,PPHO(1,NPHO))
C--
        DO 60 I=2,4
   60   PVEC(I)=0.D0
        PVEC(1)=1.D0
        CALL PHORO3(-FI3,PVEC)
        CALL PHORO2(-TH3,PVEC)
        CALL PHOBO3(ANGLE,PVEC)
        CALL PHORO3(FI4,PVEC)
        CALL PHORO2(-TH4,PNEUTR)
        CALL PHORO2(-TH4,PPHO(1,NPHO))
        CALL PHORO2(-TH4,PVEC)
        FI5=PHOAN1(PVEC(1),PVEC(2))
C--
C--   Charged particle restores original direction
        CALL PHORO3(-FI5,PNEUTR)
        CALL PHORO3(-FI5,PPHO(1,NPHO))
        CALL PHORO2(TH1,PNEUTR(1))
        CALL PHORO2(TH1,PPHO(1,NPHO))
        CALL PHORO3(FI1,PNEUTR)
        CALL PHORO3(FI1,PPHO(1,NPHO))
C--   See whether neutral system has multiplicity larger than 1...
      IF ((JDAPHO(2,IP)-JDAPHO(1,IP)).GT.1) THEN
C--   Find pointers to components of 'neutral' system
C--
        FIRST=NEUDAU
        LAST=JDAPHO(2,IP)
        DO 70 I=FIRST,LAST
          IF (I.NE.NCHARB.AND.(JMOPHO(1,I).EQ.IP)) THEN
C--
C--   Reconstruct kinematics...
            CALL PHORO3(-FI1,PPHO(1,I))
            CALL PHORO2(-TH1,PPHO(1,I))
C--
C--   ...reductive boost
            CALL PHOBO3(PARNE,PPHO(1,I))
C--
C--   Rotate in order to get photon along z-axis
            CALL PHORO3(-FI3,PPHO(1,I))
            CALL PHORO2(-TH3,PPHO(1,I))
C--
C--   Boost to the rest frame of decaying particle
            CALL PHOBO3(ANGLE,PPHO(1,I))
C--
C--   Back in the parent rest-frame but PNEUTR not yet oriented.
            CALL PHORO3(FI4,PPHO(1,I))
            CALL PHORO2(-TH4,PPHO(1,I))
C--
C--   Charged particle restores original direction
            CALL PHORO3(-FI5,PPHO(1,I))
            CALL PHORO2(TH1,PPHO(1,I))
            CALL PHORO3(FI1,PPHO(1,I))
          ENDIF
   70   CONTINUE
      ELSE
C--
C--   ...only one 'neutral' particle in addition to photon!
        DO 80 J=1,4
   80   PPHO(J,NEUDAU)=PNEUTR(J)
      ENDIF
C--
C--   All 'neutrals' treated, fill /PHOEVT/ for charged particle...
      DO 90 J=1,3
   90 PPHO(J,NCHARB)=-(PPHO(J,NPHO)+PNEUTR(J))
      PPHO(4,NCHARB)=PPHO(5,IP)-(PPHO(4,NPHO)+PNEUTR(4))
C--
      END
Commit	Line	Data
da0e9ce3	1	SUBROUTINE PHODO(IP,NCHARB,NEUDAU)
	2	C.----------------------------------------------------------------------
	3	C.
	4	C. PHOTOS: PHOton radiation in decays DOing of KINematics
	5	C.
	6	C. Purpose: Starting from the charged particle energy/momentum,
	7	C. PNEUTR, photon energy fraction and photon angle with
	8	C. respect to the axis formed by charged particle energy/
	9	C. momentum vector and PNEUTR, scale the energy/momentum,
	10	C. keeping the original direction of the neutral system in
	11	C. the lab. frame untouched.
	12	C.
	13	C. Input Parameters: IP: Pointer to decaying particle in
	14	C. /PHOEVT/ and the common itself
	15	C. NCHARB: pointer to the charged radiating
	16	C. daughter in /PHOEVT/.
	17	C. NEUDAU: pointer to the first neutral daughter
	18	C. Output Parameters: Common /PHOEVT/, with photon added.
	19	C.
	20	C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
	21	C. Last Update: 27/05/93
	22	C.
	23	C.----------------------------------------------------------------------
	24	IMPLICIT NONE
	25	DOUBLE PRECISION PHOAN1,PHOAN2,ANGLE,FI1,FI3,FI4,FI5,TH1,TH3,TH4
	26	DOUBLE PRECISION PARNE,QNEW,QOLD,DATA
	27	INTEGER IP,FI3DUM,I,J,NEUDAU,FIRST,LAST
	28	INTEGER NCHARB
	29	REAL*8 EPHOTO,PMAVIR,PHOTRI
	30	REAL*8 GNEUT,PHORAN,CCOSTH,SSINTH,PVEC(4)
	31	INTEGER NMXPHO
	32	PARAMETER (NMXPHO=10000)
	33	INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
	34	REAL*8 PPHO,VPHO
	35	COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
	36	&JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
	37	DOUBLE PRECISION MCHSQR,MNESQR
	38	REAL*8 PNEUTR
	39	COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
	40	DOUBLE PRECISION COSTHG,SINTHG
	41	REAL*8 XPHMAX,XPHOTO
	42	COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
	43	REAL*8 PI,TWOPI
	44	COMMON/PHPICO/PI,TWOPI
	45	C--
	46	EPHOTO=XPHOTO*PPHO(5,IP)/2.D0
	47	PMAVIR=SQRT(PPHO(5,IP)(PPHO(5,IP)-2.D0EPHOTO))
	48	C--
	49	C-- Reconstruct kinematics of charged particle and neutral system
	50	FI1=PHOAN1(PNEUTR(1),PNEUTR(2))
	51	C--
	52	C-- Choose axis along z of PNEUTR, calculate angle between x and y
	53	C-- components and z and x-y plane and perform Lorentz transform...
	54	TH1=PHOAN2(PNEUTR(3),SQRT(PNEUTR(1)2+PNEUTR(2)2))
	55	CALL PHORO3(-FI1,PNEUTR(1))
	56	CALL PHORO2(-TH1,PNEUTR(1))
	57	C--
	58	C-- Take away photon energy from charged particle and PNEUTR ! Thus
	59	C-- the onshell charged particle decays into virtual charged particle
	60	C-- and photon. The virtual charged particle mass becomes:
	61	C-- SQRT(PPHO(5,IP)(PPHO(5,IP)-2EPHOTO)). Construct new PNEUTR mo-
	62	C-- mentum in the rest frame of the parent:
	63	C-- 1) Scaling parameters...
	64	QNEW=PHOTRI(PMAVIR,PNEUTR(5),PPHO(5,NCHARB))
65	QOLD=PNEUTR(3)
66	GNEUT=(QNEW2+QOLD2+MNESQR)/(QNEWQOLD+SQRT((QNEW2+MNESQR)
67	&(QOLD**2+MNESQR)))
68	IF (GNEUT.LT.1.D0) THEN
69	DATA=0.D0
70	CALL PHOERR(4,'PHOKIN',DATA)
71	ENDIF
72	PARNE=GNEUT-SQRT(MAX(GNEUT**2-1.0D0,0.D0))
73	C--
74	C-- 2) ...reductive boost...
75	CALL PHOBO3(PARNE,PNEUTR)
76	C--
77	C-- ...calculate photon energy in the reduced system...
78	NPHO=NPHO+1
79	ISTPHO(NPHO)=1
80	IDPHO(NPHO) =22
81	C-- Photon mother and daughter pointers !
82	JMOPHO(1,NPHO)=IP
83	JMOPHO(2,NPHO)=0
84	JDAPHO(1,NPHO)=0
85	JDAPHO(2,NPHO)=0
86	PPHO(4,NPHO)=EPHOTO*PPHO(5,IP)/PMAVIR
87	C--
88	C-- ...and photon momenta
89	CCOSTH=-COSTHG
90	SSINTH=SINTHG
91	TH3=PHOAN2(CCOSTH,SSINTH)
92	FI3=TWOPI*PHORAN(FI3DUM)
93	PPHO(1,NPHO)=PPHO(4,NPHO)SINTHGCOS(FI3)
94	PPHO(2,NPHO)=PPHO(4,NPHO)SINTHGSIN(FI3)
95	C--
96	C-- Minus sign because axis opposite direction of charged particle !
97	PPHO(3,NPHO)=-PPHO(4,NPHO)*COSTHG
98	PPHO(5,NPHO)=0.D0
99	C--
100	C-- Rotate in order to get photon along z-axis
101	CALL PHORO3(-FI3,PNEUTR(1))
102	CALL PHORO3(-FI3,PPHO(1,NPHO))
103	CALL PHORO2(-TH3,PNEUTR(1))
104	CALL PHORO2(-TH3,PPHO(1,NPHO))
105	ANGLE=EPHOTO/PPHO(4,NPHO)
106	C--
107	C-- Boost to the rest frame of decaying particle
108	CALL PHOBO3(ANGLE,PNEUTR(1))
109	CALL PHOBO3(ANGLE,PPHO(1,NPHO))
110	C--
111	C-- Back in the parent rest frame but PNEUTR not yet oriented !
112	FI4=PHOAN1(PNEUTR(1),PNEUTR(2))
113	TH4=PHOAN2(PNEUTR(3),SQRT(PNEUTR(1)2+PNEUTR(2)2))
114	CALL PHORO3(FI4,PNEUTR(1))
115	CALL PHORO3(FI4,PPHO(1,NPHO))
116	C--
117	DO 60 I=2,4
118	60 PVEC(I)=0.D0
119	PVEC(1)=1.D0
120	CALL PHORO3(-FI3,PVEC)
121	CALL PHORO2(-TH3,PVEC)
122	CALL PHOBO3(ANGLE,PVEC)
123	CALL PHORO3(FI4,PVEC)
124	CALL PHORO2(-TH4,PNEUTR)
125	CALL PHORO2(-TH4,PPHO(1,NPHO))
126	CALL PHORO2(-TH4,PVEC)
127	FI5=PHOAN1(PVEC(1),PVEC(2))
128	C--
129	C-- Charged particle restores original direction
130	CALL PHORO3(-FI5,PNEUTR)
131	CALL PHORO3(-FI5,PPHO(1,NPHO))
132	CALL PHORO2(TH1,PNEUTR(1))
133	CALL PHORO2(TH1,PPHO(1,NPHO))
134	CALL PHORO3(FI1,PNEUTR)
135	CALL PHORO3(FI1,PPHO(1,NPHO))
136	C-- See whether neutral system has multiplicity larger than 1...
137	IF ((JDAPHO(2,IP)-JDAPHO(1,IP)).GT.1) THEN
138	C-- Find pointers to components of 'neutral' system
139	C--
140	FIRST=NEUDAU
141	LAST=JDAPHO(2,IP)
142	DO 70 I=FIRST,LAST
143	IF (I.NE.NCHARB.AND.(JMOPHO(1,I).EQ.IP)) THEN
144	C--
145	C-- Reconstruct kinematics...
146	CALL PHORO3(-FI1,PPHO(1,I))
147	CALL PHORO2(-TH1,PPHO(1,I))
148	C--
149	C-- ...reductive boost
150	CALL PHOBO3(PARNE,PPHO(1,I))
151	C--
152	C-- Rotate in order to get photon along z-axis
153	CALL PHORO3(-FI3,PPHO(1,I))
154	CALL PHORO2(-TH3,PPHO(1,I))
155	C--
156	C-- Boost to the rest frame of decaying particle
157	CALL PHOBO3(ANGLE,PPHO(1,I))
158	C--
159	C-- Back in the parent rest-frame but PNEUTR not yet oriented.
160	CALL PHORO3(FI4,PPHO(1,I))
161	CALL PHORO2(-TH4,PPHO(1,I))
162	C--
163	C-- Charged particle restores original direction
164	CALL PHORO3(-FI5,PPHO(1,I))
165	CALL PHORO2(TH1,PPHO(1,I))
166	CALL PHORO3(FI1,PPHO(1,I))
167	ENDIF
168	70 CONTINUE
169	ELSE
170	C--
171	C-- ...only one 'neutral' particle in addition to photon!
172	DO 80 J=1,4
173	80 PPHO(J,NEUDAU)=PNEUTR(J)
174	ENDIF
175	C--
176	C-- All 'neutrals' treated, fill /PHOEVT/ for charged particle...
177	DO 90 J=1,3
178	90 PPHO(J,NCHARB)=-(PPHO(J,NPHO)+PNEUTR(J))
179	PPHO(4,NCHARB)=PPHO(5,IP)-(PPHO(4,NPHO)+PNEUTR(4))
180	C--
181	END