photos_hepevt_example.f
1C********************************************************
2C Example of use of Photos C++ interface.
3C D+ D- -> tau+ tau- HEPEVT events are constructed.
4C Taus are subsequently decayed via Photos++.
5C
6C @author Tomasz Przedzinski
7C @date 21 August 2013
8C********************************************************
9 PROGRAM ph_hepvt_test
10
11C INITIALIZE PHOTOS++
12 CALL photos_init
13
14C PREPARE AN EVENT in HEPEVT
15C TWO TEST EVENTS HAVE BEEN PROVIDED WITH THIS EXAMPLE
16 CALL simple_ee_to_z_to_tautau_event
17c CALL EVEN_SIMPLER_Z_TO_EE_EVENT
18
19 WRITE(*,*) "##############"
20 WRITE(*,*) "PHODMP: BEFORE"
21 WRITE(*,*) "##############"
22 CALL phodmp
23
24 CALL photos_process
25C CALL PHOTOS_PROCESS_PARTICLE(4)
26C CALL PHOTOS_PROCESS_BRANCH(4)
27
28 WRITE(*,*) "##############"
29 WRITE(*,*) "PHODMP: AFTER"
30 WRITE(*,*) "##############"
31 CALL phodmp
32 END
33
34 SUBROUTINE simple_ee_to_z_to_tautau_event
35C serve as an exaple; to be replaced with event generated by user
36 real*8 amell
37
38 amell = 0.0005111;
39
40C CREATE SIMPLE EVENT: e+ e- -> Z -> tau+ tau- -> pi nu_tau, pi nu_tau
41C ---------------------------------------------------------------------
42 CALL add_particle( 11, 6, 1.7763568394002505d-15, -3.5565894425761324d-15, 4.5521681043409913d+01, 4.5521681043409934d+01,
43 $ amell, 0, 0, 3, 3)
44 CALL add_particle( -11, 6, -1.7763568394002505d-15, 3.5488352204797800d-15, -4.5584999071936601d+01, 4.5584999071936622d+01,
45 $ amell, 0, 0, 3, 3)
46 CALL add_particle( 23, 5, 0d0, 0d0, -6.3318028526687442d-02, 9.1106680115346506d+01,
47 $ 9.1106658112716090d+01, 1, 2, 4, 5)
48 CALL add_particle( 15, 2, -2.3191595992562256d+01, -2.6310500920665142d+01, -2.9046412466624929d+01, 4.5573504956498098d+01,
49 $ 1.7769900000002097d+00, 3, 0, 6, 7)
50 CALL add_particle( -15, 2, 2.3191595992562256d+01, 2.6310500920665142d+01, 2.8983094438098242d+01, 4.5533175158848429d+01,
51 $ 1.7769900000000818d+00, 3, 0, 8, 9)
52 CALL add_particle( 16, 1, -1.2566536214715378d+00, -1.7970251138317268d+00, -1.3801323581022720d+00, 2.5910119010468553d+00,
53 $ 9.9872238934040070d-03, 4, 0, 0, 0)
54 CALL add_particle(-211, 1, -2.1935073012334062d+01, -2.4513624017269400d+01, -2.7666443730700312d+01, 4.2982749776866747d+01,
55 $ 1.3956783711910248d-01, 4, 0, 0, 0)
56 CALL add_particle( -16, 1, 8.4364531743909055d+00, 8.3202830831667836d+00, 9.6202800273055971d+00, 1.5262723881157640d+01,
57 $ 9.9829332903027534d-03, 5, 0, 0, 0)
58 CALL add_particle( 211, 1, 1.4755273459419701d+01, 1.7990366047940022d+01, 1.9362977676297948d+01, 3.0270707771933196d+01,
59 $ 1.3956753909587860d-01, 5, 0, 0, 0)
60 END
61
62 SUBROUTINE even_simpler_z_to_ee_event
63C serve as an exaple; to be replaced with event generated by user
64 real*8 amell
65
66 amell = 0.0005111;
67
68C CREATE SIMPLE EVENT: Z -> e+ e-
69C ---------------------------------------------------------------------
70 CALL add_particle( 23, 2, 0.000000000000000d+00, 0.000000000000000d+00,
71 $ -4.821915360414504d+02, 4.900549668552968d+02, 9.118760000000000d+01, 0, 0, 2, 3)
72 CALL add_particle( 11, 1, -3.018114335125028d+01, -1.204071905454172d+01,
73 $ -7.717282429699088d+01, 8.373485020774415d+01, 5.109989200000000d-04, 1, 0, 0, 0)
74 CALL add_particle(-11, 1, 3.018114335125028d+01, 1.204071905454172d+01,
75 $ -4.050187117444596d+02, 4.063201166475526d+02, 5.109989200000000d-04, 1, 0, 0, 0)
76 END
77
78 SUBROUTINE add_particle(ID,STATUS,PX,PY,PZ,E,M,MOTHER1,MOTHER2,DAUGHTER1,DAUGHTER2)
79 INTEGER ID,STATUS,MOTHER1,MOTHER2,DAUGHTER1,DAUGHTER2
80 real*8 px,py,pz,e,m
81C-----------------------------------------------------------------------------
82C COMMON HEPEVT
83C
84C IMPORTANT: The definition of HEPEVT is also present in:
85C TAUOLA/src/eventRecords/PhotosHEPEVTEvent.cxx
86C If the definition changes (eg. different value of NMXHEP or REAL
87C instead of REAL*8) it has to be updated in that file as well
88C and Tauola++ has to be recompiled.
89C-----------------------------------------------------------------------------
90 INTEGER NMXHEP
91 parameter(nmxhep=10000)
92 real*8 phep, vhep ! to be real*4/ *8 depending on host
93 INTEGER nevhep,nhep,isthep,idhep,jmohep,
94 $ jdahep
95 COMMON /hepevt/
96 $ nevhep, ! serial number
97 $ nhep, ! number of particles
98 $ isthep(nmxhep), ! status code
99 $ idhep(nmxhep), ! particle ident KF
100 $ jmohep(2,nmxhep), ! parent particles
101 $ jdahep(2,nmxhep), ! childreen particles
102 $ phep(5,nmxhep), ! four-momentum, mass [GeV]
103 $ vhep(4,nmxhep) ! vertex [mm]
104
105C
106C WARNING: note that common PHOQED is missing in this example.
107C it is because its content, so far, is not used in C++
108C implementations.
109C Its function is taken by other methods, which are probably more
110C comfortable and elegant.
111C
112
113 SAVE hepevt
114 nhep=nhep+1
115 idhep(nhep) =id
116 isthep(nhep)=status
117 phep(1,nhep)=px
118 phep(2,nhep)=py
119 phep(3,nhep)=pz
120 phep(4,nhep)=e
121 phep(5,nhep)=m
122 jmohep(1,nhep)=mother1
123 jmohep(2,nhep)=mother2
124 jdahep(1,nhep)=daughter1
125 jdahep(2,nhep)=daughter2
126 END
127
128 SUBROUTINE phodmp
129C.----------------------------------------------------------------------
130C.
131C. PHOTOS: PHOton radiation in decays event DuMP routine
132C.
133C. Purpose: Print event record.
134C.
135C. Input Parameters: Common /HEPEVT/
136C.
137C. Output Parameters: None
138C.
139C. Author(s): B. van Eijk Created at: 05/06/90
140C. Last Update: 05/06/90
141C.
142C.----------------------------------------------------------------------
143C IMPLICIT NONE
144 DOUBLE PRECISION SUMVEC(5)
145 INTEGER I,J
146C this is the hepevt class in old style. No d_h_ class prD-name
147 INTEGER NMXHEP
148 parameter(nmxhep=10000)
149 real*8 phep, vhep ! to be real*4/ *8 depending on host
150 INTEGER nevhep,nhep,isthep,idhep,jmohep,
151 $ jdahep
152 COMMON /hepevt/
153 $ nevhep, ! serial number
154 $ nhep, ! number of particles
155 $ isthep(nmxhep), ! status code
156 $ idhep(nmxhep), ! particle ident KF
157 $ jmohep(2,nmxhep), ! parent particles
158 $ jdahep(2,nmxhep), ! childreen particles
159 $ phep(5,nmxhep), ! four-momentum, mass [GeV]
160 $ vhep(4,nmxhep) ! vertex [mm]
161
162 INTEGER PHLUN
163 common/pholun/phlun
164 DO 10 i=1,5
165 10 sumvec(i)=0.
166C--
167C-- Print event number...
168 WRITE(phlun,9000)
169 WRITE(phlun,9010) nevhep
170 WRITE(phlun,9080)
171 WRITE(phlun,9020)
172 DO 30 i=1,nhep
173C--
174C-- For 'stable particle' calculate vector momentum sum
175 IF (jdahep(1,i).EQ.0) THEN
176 DO 20 j=1,4
177 20 sumvec(j)=sumvec(j)+phep(j,i)
178 IF (jmohep(2,i).EQ.0) THEN
179 WRITE(phlun,9030) i,idhep(i),jmohep(1,i),(phep(j,i),j=1,5)
180 ELSE
181 WRITE(phlun,9040) i,idhep(i),jmohep(1,i),jmohep(2,i),(phep
182 & (j,i),j=1,5)
183 ENDIF
184 ELSE
185 IF (jmohep(2,i).EQ.0) THEN
186 WRITE(phlun,9050) i,idhep(i),jmohep(1,i),jdahep(1,i),
187 & jdahep(2,i),(phep(j,i),j=1,5)
188 ELSE
189 WRITE(phlun,9060) i,idhep(i),jmohep(1,i),jmohep(2,i),
190 & jdahep(1,i),jdahep(2,i),(phep(j,i),j=1,5)
191 ENDIF
192 ENDIF
193 30 CONTINUE
194 sumvec(5)=sqrt(sumvec(4)**2-sumvec(1)**2-sumvec(2)**2-
195 &sumvec(3)**2)
196 WRITE(phlun,9070) (sumvec(j),j=1,5)
197 RETURN
198 9000 FORMAT(1h0,80('='))
199 9010 FORMAT(1h ,29x,'Event No.:',i10)
200 9020 FORMAT(1h0,1x,'Nr',3x,'Type',3x,'Parent(s)',2x,'Daughter(s)',6x,
201 &'Px',7x,'Py',7x,'Pz',7x,'E',4x,'Inv. M.')
202 9030 FORMAT(1h ,i4,i7,3x,i4,9x,'Stable',2x,5f9.2)
203 9040 FORMAT(1h ,i4,i7,i4,' - ',i4,5x,'Stable',2x,5f9.2)
204 9050 FORMAT(1h ,i4,i7,3x,i4,6x,i4,' - ',i4,5f9.2)
205 9060 FORMAT(1h ,i4,i7,i4,' - ',i4,2x,i4,' - ',i4,5f9.2)
206 9070 FORMAT(1h0,23x,'Vector Sum: ', 5f9.2)
207 9080 FORMAT(1h0,6x,'Particle Parameters')
208 END