photos_hepevt_example.f
1 C********************************************************
2 C Example of use of Photos C++ interface.
3 C D+ D- -> tau+ tau- HEPEVT events are constructed.
4 C Taus are subsequently decayed via Photos++.
5 C
6 C @author Tomasz Przedzinski
7 C @date 21 August 2013
8 C********************************************************
9  PROGRAM ph_hepvt_test
10 
11 C INITIALIZE PHOTOS++
12  CALL photos_init
13 
14 C PREPARE AN EVENT in HEPEVT
15 C TWO TEST EVENTS HAVE BEEN PROVIDED WITH THIS EXAMPLE
16  CALL simple_ee_to_z_to_tautau_event
17 c CALL EVEN_SIMPLER_Z_TO_EE_EVENT
18 
19  WRITE(*,*) "##############"
20  WRITE(*,*) "PHODMP: BEFORE"
21  WRITE(*,*) "##############"
22  CALL phodmp
23 
24  CALL photos_process
25 C CALL PHOTOS_PROCESS_PARTICLE(4)
26 C 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
35 C serve as an exaple; to be replaced with event generated by user
36  REAL*8 amell
37 
38  amell = 0.0005111;
39 
40 C CREATE SIMPLE EVENT: e+ e- -> Z -> tau+ tau- -> pi nu_tau, pi nu_tau
41 C ---------------------------------------------------------------------
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
63 C serve as an exaple; to be replaced with event generated by user
64  REAL*8 amell
65 
66  amell = 0.0005111;
67 
68 C CREATE SIMPLE EVENT: Z -> e+ e-
69 C ---------------------------------------------------------------------
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
81 C-----------------------------------------------------------------------------
82 C COMMON HEPEVT
83 C
84 C IMPORTANT: The definition of HEPEVT is also present in:
85 C TAUOLA/src/eventRecords/PhotosHEPEVTEvent.cxx
86 C If the definition changes (eg. different value of NMXHEP or REAL
87 C instead of REAL*8) it has to be updated in that file as well
88 C and Tauola++ has to be recompiled.
89 C-----------------------------------------------------------------------------
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 
105 C
106 C WARNING: note that common PHOQED is missing in this example.
107 C it is because its content, so far, is not used in C++
108 C implementations.
109 C Its function is taken by other methods, which are probably more
110 C comfortable and elegant.
111 C
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
129 C.----------------------------------------------------------------------
130 C.
131 C. PHOTOS: PHOton radiation in decays event DuMP routine
132 C.
133 C. Purpose: Print event record.
134 C.
135 C. Input Parameters: Common /HEPEVT/
136 C.
137 C. Output Parameters: None
138 C.
139 C. Author(s): B. van Eijk Created at: 05/06/90
140 C. Last Update: 05/06/90
141 C.
142 C.----------------------------------------------------------------------
143 C IMPLICIT NONE
144  DOUBLE PRECISION sumvec(5)
145  INTEGER i,j
146 C 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.
166 C--
167 C-- 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
173 C--
174 C-- 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