[u/mrichter/AliRoot.git] / PHOS / reconstruction / shower_functions_v3.f

c ==================================================================

	function coornew(el,er,angl,mmm)
	mmnew=mmconvert(mmm)	
	esum=el+er
	ex=el/esum

	if(esum.le.0.) then
		write (*,*) ' coornew diagnostic bad esum=',esum
		return
	endif
	
	itip=newtipcell(mmnew)
	if(itip.eq.1) then
		xx=xcryexgl(ex)
	elseif(itip.eq.2) then
c		xx=xcryexcr(ex)
		xx=xcryexcrangl(ex,angl)
	endif
	coornew=xx
	return
	end	
	
	function xcryexgl(ar)
c FOR LEAD GLASS
	data p1,p2 /1.480, 1.600/
c ar - energy fraction right from boundary
	

	if(ar.gt.0.5) then
		a=1.-ar
		s=-1.
	else
		a=ar
		s=1.
	endif
	if(a.le.0.) a=0.0001

	y=alog(2.*a)
	x=p2*y**2-p1*y
c lead glass for my geant
c	x=x/0.55
c lead glass for my geant 2*2 cells
	x=x/0.65
c corrections for crystalls 2*2, systematic is better than 0.25 mm for 3Gev geant showers
c94	x=x*0.54-0.25
c end corrections for crystalls
	xcryexgl=s*x
	return
	end

	function xcryexcr(ar)
c FOR CRYSTALS
c	data p1,p2 /1.480, 1.600/
c  very new 26.02.98
	data p1,p2 /1.700, 1.150/
c ar - energy fraction right from boundary
	if(ar.gt.0.5) then
		a=1.-ar
		s=-1.
	else
		a=ar
		s=1.
	endif
	if(a.le.0.) a=0.0001
	y=alog(2.*a)
	x=p2*y**2-p1*y
	xcryexcr=s*x
	return
	end

	function xcryexcrangl(ar,angl)
c FOR CRYSTALS ALICE parametrisation 23.02.98
c	data p1,p2 /1.480, 1.600/
	data p1,p2 /1.700, 1.150/
c ar - energy fraction right from boundary
        aangl=abs(angl)
c        write (*,*) ' ar,angl ',ar,angl
        if(angl.ge.0.) then
                pp1=p1
                pp2=p2-0.7*aangl
        else
                pp1=p1+7.6*aangl
                pp2=p2-0.9*aangl
        endif
c        if(angl.ge.0.) then
c                pp1=p1+7.6*aangl
c                pp2=p2-0.9*aangl
c        else
c                pp1=p1
c                pp2=p2-0.7*aangl
c        endif
	if(ar.gt.0.5) then
		a=1.-ar
		s=-1.
	else
		a=ar
		s=1.
	endif
	if(a.le.0.) a=0.0001
	y=alog(2.*a)
	x=pp2*y**2-pp1*y
	xcryexcrangl=s*x
	return
	end
c ==============================================================
	real function cumulnew(xc,itip)
        real xc
        integer itip

	if(itip.eq.1) then
		a=cucryexgl(xc)
	elseif(itip.eq.2) then
		a=cucryexcr(xc)
	endif
        cumulnew=a
        return
        end
        
	real function dcumulnew(xc,itip)
        real xc
        integer itip

	if(itip.eq.1) then
		a=dcucryexgl(xc)
	elseif(itip.eq.2) then
		a=dcucryexcr(xc)
	endif
        dcumulnew=a
        return
        end
        
	real function cucryexcr(xc)
c fit for ALICE 94 RUN DATA
c	data p1,p2 /2.526, 1.104/
c94	data p1,p2 /2.300, 2.044/
c crystals for my geant
c	data p1,p2 /1.480, 1.600/
c  very new 26.02.98
	data p1,p2 /1.700, 1.150/
	x=abs(xc)
	a=0.5*exp((p1-sqrt(p1**2+4.*p2*x))/2./p2)
	if(xc.ge.0.) then
		cucryexcr=1.-a
	else
		cucryexcr=a
	endif
	return
	end

	real function dcucryexcr(xc)
c fit for ALICE 94 RUN DATA
c	data p1,p2 /2.526, 1.104/
c	data p1,p2 /2.300, 2.044/
c crystals for my geant
c	data p1,p2 /1.480, 1.600/
c  very new 26.02.98
	data p1,p2 /1.700, 1.150/
	x=abs(xc)

	rr=sqrt(p1**2+4.*p2*x)
	da=0.5*exp((p1-rr)/2./p2)/rr
	dcucryexcr=da
	
	return
	end

	real function cucryexgl(xc)
c fit for ALICE 94 RUN DATA
c	data p1,p2 /2.526, 1.104/
c94	data p1,p2 /2.300, 2.044/
c crystals for my geant
	data p1,p2 /1.480, 1.600/
	x=abs(xc)
c lead glass for my geant
	x=x*0.55

	a=0.5*exp((p1-sqrt(p1**2+4.*p2*x))/2./p2)
	if(xc.ge.0.) then
		cucryexgl=1.-a
	else
		cucryexgl=a
	endif
	return
	end

	real function dcucryexgl(xc)
c fit for ALICE 94 RUN DATA
c	data p1,p2 /2.526, 1.104/
c	data p1,p2 /2.300, 2.044/
c crystals for my geant
	data p1,p2 /1.480, 1.600/
	x=abs(xc)
c lead glass for my geant
	x=x*0.55

	rr=sqrt(p1**2+4.*p2*x)
	da=0.5*exp((p1-rr)/2./p2)/rr
c lead glass for my geant
	da=da*0.55
	dcucryexgl=da
	
	return
	end
c ==============================================================

	real function ampcelnew(e,vx,vc,vb,mmm)
	common /comkey/ keykey
	real vx(5),vc(3),vb(3)
        integer mmm

	mmnew=mmconvert(mmm)	
	itip=newtipcell(mmnew)
	xg=vx(1)
	yg=vx(2)
	xc=vc(1)
	yc=vc(2)
	if(keykey.eq.0) then
        	a=aclnew(xg-xc,yg-yc,vb,itip)
	elseif(keykey.eq.1) then
c new with angle
        	a=acl3(e,vx,vc,vb)
	else
		write (*,*) ' Net takoi keykey !! ',keykey
                a=0.
	endif	
        ampcelnew=a
        return
        end

	real function dampcelnew(nd,e,vx,vc,vb,mmm)
	common /comkey/keykey
	real vx(5),vc(3),vb(3)
        integer mmm

	mmnew=mmconvert(mmm)	
	itip=newtipcell(mmnew)
c old no angle
	if(keykey.eq.0) then
		xg=vx(1)
		yg=vx(2)
		xc=vc(1)
		yc=vc(2)
	   if(nd.eq.1) then
		da=dxaclnew(xg-xc,yg-yc,vb,itip)
	   elseif(nd.eq.2) then
		da=dxaclnew(yg-yc,xg-xc,vb,itip)
	   else
		write (*,*) ' net takoi nd ',nd
	   endif
	elseif(keykey.eq.1) then
c new with angle
		da=dacl3(nd,e,vx,vc,vb)
	else
		write (*,*) ' net takoi keykey ',keykey
                da=0.
	endif
        dampcelnew=da

	return
	end

        FUNCTION aclnew(XG,YG,DA,mmn)
C
C     SHOWER FRACTION IN THE CELL
C       XG,YG GAMMA COOR. X=0,Y=0 - CENTER OF THE CELL
C       DA - SIZE OF THE CELL (MM)
C
	real da(3)
        integer mmn
        
        EXTERNAL a2fnew
C       SX=SIGN(1.,XG)
C       SY=SIGN(1.,YG)    FOR ANGLE
        X=-ABS(XG)
        Y=-ABS(YG)
        Dx=DA(1)/2.
	dy=da(2)/2.
        XP=X+Dx
        YP=Y+Dy
        XM=X-Dx
        YM=Y-Dy
        aclnew=a2fnew(XP,YP,mmn)+a2fnew(XM,YM,mmn)-
     -         a2fnew(XP,YM,mmn)-a2fnew(XM,YP,mmn)
c	write (*,*) ' xg,yg,dx,dy,acl ',xg,yg,dx,dy,acl
        RETURN
        END
C
        FUNCTION dxaclnew(XG,YG,DA,mmn)
C
C       DXACL=D(ACL(X,Y,DA))/D(X)
C
        integer mmn
        EXTERNAL dxa2fnew
        SX=SIGN(1.,XG)
C       SY=SIGN(1.,YG)    FOR ANGLE
        X=-ABS(XG)
        Y=-ABS(YG)
        D=DA/2.
        XP=X+D
        YP=Y+D
        XM=X-D
        YM=Y-D
        A=dxa2fnew(XP,YP,mmn)+dxa2fnew(XM,YM,mmn)-
     -    dxa2fnew(XP,YM,mmn)-dxa2fnew(XM,YP,mmn)
        dxaclnew=-SX*A
        RETURN
        END

        FUNCTION a2fnew(X,Y,mmn)
        integer mmn
        EXTERNAL afnew
        AX=afnew(X,mmn)
        AY=afnew(Y,mmn)
	a2fnew=aa2fz(ax,ay)
	return
        END

        FUNCTION dxa2fnew(X,Y,mmn)
        integer mmn
        EXTERNAL afnew,dafnew
        AX=afnew(X,mmn)
        AY=afnew(Y,mmn)
	dxa2fnew=daa2fz(ax,ay)*dafnew(x,mmn)
	return
        END

        FUNCTION afnew(x,mmn)
	afnew=cumulnew(x,mmn)
	return
        end
        
        FUNCTION dafnew(x,mmn)
	dafnew=dcumulnew(x,mmn)
	return
        end

	function sigwlgam0(ewl,egam)
c        write (*,*) ' ewl,egam ',ewl,egam 
	if(egam.le.0.) then
		sigwlgam0=0.
	else
		sqwl=sqrt(ewl)
		fr=ewl/egam
		if(fr.gt.0.8) fr=0.8
c		sigwlgam0=0.36*(1.-fr/0.88)*sqwl
        xp=abs(fr-0.05)
        sfun=0.5*sin(3.14/0.7*xp)+0.5*sin(3.14/0.76*sqrt(xp*0.7))
c        sfun=2.5*egam*sfun+0.13-0.13*xp
c        fm=(1.+exp(-egam/0.2))
c        if(fm.lt.1.) fm=1.
c        if(fm.gt.2.) fm=2.
        fm=2.5
        sfun=egam*fm*sfun+0.13-0.13*xp
c        sigwlgam0=sqrt(egam)*sfun*sqwl
        sigwlgam0=0.1/sqrt(egam)*sfun*sqwl
c        write (*,*) ' fr sfun ',fr,sfun 
c        write (*,*) ' sigwlgam0 ',sigwlgam0,ewl,egam 

	endif
	return
	end	

	function dispwlm(ewl,m)

        common /comsigma/sigmaph,sigmapd,sigphsq,sigpdsq


	dispph=sigphsq*ewl
	disppd=sigpdsq	
	dispwlm=dispph+disppd

	return
	end

        real function zmidlshower(ee,m)
	include 'comgam.for'
        zmidlshower = 46.
c        zmidlshower = 50.
        return
        end
Commit	Line	Data
fe4da5cc	1	c ==================================================================
	2
	3	function coornew(el,er,angl,mmm)
	4	mmnew=mmconvert(mmm)
	5	esum=el+er
	6	ex=el/esum
	7
	8	if(esum.le.0.) then
	9	write (,) ' coornew diagnostic bad esum=',esum
	10	return
	11	endif
	12
	13	itip=newtipcell(mmnew)
	14	if(itip.eq.1) then
	15	xx=xcryexgl(ex)
	16	elseif(itip.eq.2) then
	17	c xx=xcryexcr(ex)
	18	xx=xcryexcrangl(ex,angl)
	19	endif
	20	coornew=xx
	21	return
	22	end
	23
	24	function xcryexgl(ar)
	25	c FOR LEAD GLASS
	26	data p1,p2 /1.480, 1.600/
	27	c ar - energy fraction right from boundary
	28
	29
	30	if(ar.gt.0.5) then
	31	a=1.-ar
	32	s=-1.
	33	else
	34	a=ar
	35	s=1.
	36	endif
	37	if(a.le.0.) a=0.0001
	38
	39	y=alog(2.*a)
	40	x=p2y2-p1y
	41	c lead glass for my geant
	42	c x=x/0.55
	43	c lead glass for my geant 2*2 cells
	44	x=x/0.65
	45	c corrections for crystalls 2*2, systematic is better than 0.25 mm for 3Gev geant showers
	46	c94 x=x*0.54-0.25
	47	c end corrections for crystalls
	48	xcryexgl=s*x
	49	return
	50	end
	51
	52	function xcryexcr(ar)
	53	c FOR CRYSTALS
	54	c data p1,p2 /1.480, 1.600/
	55	c very new 26.02.98
	56	data p1,p2 /1.700, 1.150/
	57	c ar - energy fraction right from boundary
	58	if(ar.gt.0.5) then
	59	a=1.-ar
	60	s=-1.
	61	else
	62	a=ar
	63	s=1.
	64	endif
65	if(a.le.0.) a=0.0001
66	y=alog(2.*a)
67	x=p2y2-p1y
68	xcryexcr=s*x
69	return
70	end
71
72	function xcryexcrangl(ar,angl)
73	c FOR CRYSTALS ALICE parametrisation 23.02.98
74	c data p1,p2 /1.480, 1.600/
75	data p1,p2 /1.700, 1.150/
76	c ar - energy fraction right from boundary
77	aangl=abs(angl)
78	c write (,) ' ar,angl ',ar,angl
79	if(angl.ge.0.) then
80	pp1=p1
81	pp2=p2-0.7*aangl
82	else
83	pp1=p1+7.6*aangl
84	pp2=p2-0.9*aangl
85	endif
86	c if(angl.ge.0.) then
87	c pp1=p1+7.6*aangl
88	c pp2=p2-0.9*aangl
89	c else
90	c pp1=p1
91	c pp2=p2-0.7*aangl
92	c endif
93	if(ar.gt.0.5) then
94	a=1.-ar
95	s=-1.
96	else
97	a=ar
98	s=1.
99	endif
100	if(a.le.0.) a=0.0001
101	y=alog(2.*a)
102	x=pp2y2-pp1y
103	xcryexcrangl=s*x
104	return
105	end
106	c ==============================================================
107	real function cumulnew(xc,itip)
108	real xc
109	integer itip
110
111	if(itip.eq.1) then
112	a=cucryexgl(xc)
113	elseif(itip.eq.2) then
114	a=cucryexcr(xc)
115	endif
116	cumulnew=a
117	return
118	end
119
120	real function dcumulnew(xc,itip)
121	real xc
122	integer itip
123
124	if(itip.eq.1) then
125	a=dcucryexgl(xc)
126	elseif(itip.eq.2) then
127	a=dcucryexcr(xc)
128	endif
129	dcumulnew=a
130	return
131	end
132
133	real function cucryexcr(xc)
134	c fit for ALICE 94 RUN DATA
135	c data p1,p2 /2.526, 1.104/
136	c94 data p1,p2 /2.300, 2.044/
137	c crystals for my geant
138	c data p1,p2 /1.480, 1.600/
139	c very new 26.02.98
140	data p1,p2 /1.700, 1.150/
141	x=abs(xc)
142	a=0.5exp((p1-sqrt(p12+4.p2*x))/2./p2)
143	if(xc.ge.0.) then
144	cucryexcr=1.-a
145	else
146	cucryexcr=a
147	endif
148	return
149	end
150
151	real function dcucryexcr(xc)
152	c fit for ALICE 94 RUN DATA
153	c data p1,p2 /2.526, 1.104/
154	c data p1,p2 /2.300, 2.044/
155	c crystals for my geant
156	c data p1,p2 /1.480, 1.600/
157	c very new 26.02.98
158	data p1,p2 /1.700, 1.150/
159	x=abs(xc)
160
161	rr=sqrt(p1*2+4.p2*x)
162	da=0.5*exp((p1-rr)/2./p2)/rr
163	dcucryexcr=da
164
165	return
166	end
167
168	real function cucryexgl(xc)
169	c fit for ALICE 94 RUN DATA
170	c data p1,p2 /2.526, 1.104/
171	c94 data p1,p2 /2.300, 2.044/
172	c crystals for my geant
173	data p1,p2 /1.480, 1.600/
174	x=abs(xc)
175	c lead glass for my geant
176	x=x*0.55
177
178	a=0.5exp((p1-sqrt(p12+4.p2*x))/2./p2)
179	if(xc.ge.0.) then
180	cucryexgl=1.-a
181	else
182	cucryexgl=a
183	endif
184	return
185	end
186
187	real function dcucryexgl(xc)
188	c fit for ALICE 94 RUN DATA
189	c data p1,p2 /2.526, 1.104/
190	c data p1,p2 /2.300, 2.044/
191	c crystals for my geant
192	data p1,p2 /1.480, 1.600/
193	x=abs(xc)
194	c lead glass for my geant
195	x=x*0.55
196
197	rr=sqrt(p1*2+4.p2*x)
198	da=0.5*exp((p1-rr)/2./p2)/rr
199	c lead glass for my geant
200	da=da*0.55
201	dcucryexgl=da
202
203	return
204	end
205	c ==============================================================
206
207	real function ampcelnew(e,vx,vc,vb,mmm)
208	common /comkey/ keykey
209	real vx(5),vc(3),vb(3)
210	integer mmm
211
212	mmnew=mmconvert(mmm)
213	itip=newtipcell(mmnew)
214	xg=vx(1)
215	yg=vx(2)
216	xc=vc(1)
217	yc=vc(2)
218	if(keykey.eq.0) then
219	a=aclnew(xg-xc,yg-yc,vb,itip)
220	elseif(keykey.eq.1) then
221	c new with angle
222	a=acl3(e,vx,vc,vb)
223	else
224	write (,) ' Net takoi keykey !! ',keykey
225	a=0.
226	endif
227	ampcelnew=a
228	return
229	end
230
231	real function dampcelnew(nd,e,vx,vc,vb,mmm)
232	common /comkey/keykey
233	real vx(5),vc(3),vb(3)
234	integer mmm
235
236	mmnew=mmconvert(mmm)
237	itip=newtipcell(mmnew)
238	c old no angle
239	if(keykey.eq.0) then
240	xg=vx(1)
241	yg=vx(2)
242	xc=vc(1)
243	yc=vc(2)
244	if(nd.eq.1) then
245	da=dxaclnew(xg-xc,yg-yc,vb,itip)
246	elseif(nd.eq.2) then
247	da=dxaclnew(yg-yc,xg-xc,vb,itip)
248	else
249	write (,) ' net takoi nd ',nd
250	endif
251	elseif(keykey.eq.1) then
252	c new with angle
253	da=dacl3(nd,e,vx,vc,vb)
254	else
255	write (,) ' net takoi keykey ',keykey
256	da=0.
257	endif
258	dampcelnew=da
259
260	return
261	end
262
263	FUNCTION aclnew(XG,YG,DA,mmn)
264	C
265	C SHOWER FRACTION IN THE CELL
266	C XG,YG GAMMA COOR. X=0,Y=0 - CENTER OF THE CELL
267	C DA - SIZE OF THE CELL (MM)
268	C
269	real da(3)
270	integer mmn
271
272	EXTERNAL a2fnew
273	C SX=SIGN(1.,XG)
274	C SY=SIGN(1.,YG) FOR ANGLE
275	X=-ABS(XG)
276	Y=-ABS(YG)
277	Dx=DA(1)/2.
278	dy=da(2)/2.
279	XP=X+Dx
280	YP=Y+Dy
281	XM=X-Dx
282	YM=Y-Dy
283	aclnew=a2fnew(XP,YP,mmn)+a2fnew(XM,YM,mmn)-
284	- a2fnew(XP,YM,mmn)-a2fnew(XM,YP,mmn)
285	c write (,) ' xg,yg,dx,dy,acl ',xg,yg,dx,dy,acl
286	RETURN
287	END
288	C
289	FUNCTION dxaclnew(XG,YG,DA,mmn)
290	C
291	C DXACL=D(ACL(X,Y,DA))/D(X)
292	C
293	integer mmn
294	EXTERNAL dxa2fnew
295	SX=SIGN(1.,XG)
296	C SY=SIGN(1.,YG) FOR ANGLE
297	X=-ABS(XG)
298	Y=-ABS(YG)
299	D=DA/2.
300	XP=X+D
301	YP=Y+D
302	XM=X-D
303	YM=Y-D
304	A=dxa2fnew(XP,YP,mmn)+dxa2fnew(XM,YM,mmn)-
305	- dxa2fnew(XP,YM,mmn)-dxa2fnew(XM,YP,mmn)
306	dxaclnew=-SX*A
307	RETURN
308	END
309
310	FUNCTION a2fnew(X,Y,mmn)
311	integer mmn
312	EXTERNAL afnew
313	AX=afnew(X,mmn)
314	AY=afnew(Y,mmn)
315	a2fnew=aa2fz(ax,ay)
316	return
317	END
318
319	FUNCTION dxa2fnew(X,Y,mmn)
320	integer mmn
321	EXTERNAL afnew,dafnew
322	AX=afnew(X,mmn)
323	AY=afnew(Y,mmn)
324	dxa2fnew=daa2fz(ax,ay)*dafnew(x,mmn)
325	return
326	END
327
328	FUNCTION afnew(x,mmn)
329	afnew=cumulnew(x,mmn)
330	return
331	end
332
333	FUNCTION dafnew(x,mmn)
334	dafnew=dcumulnew(x,mmn)
335	return
336	end
337
338	function sigwlgam0(ewl,egam)
339	c write (,) ' ewl,egam ',ewl,egam
340	if(egam.le.0.) then
341	sigwlgam0=0.
342	else
343	sqwl=sqrt(ewl)
344	fr=ewl/egam
345	if(fr.gt.0.8) fr=0.8
346	c sigwlgam0=0.36(1.-fr/0.88)sqwl
347	xp=abs(fr-0.05)
348	sfun=0.5sin(3.14/0.7xp)+0.5sin(3.14/0.76sqrt(xp*0.7))
349	c sfun=2.5egamsfun+0.13-0.13*xp
350	c fm=(1.+exp(-egam/0.2))
351	c if(fm.lt.1.) fm=1.
352	c if(fm.gt.2.) fm=2.
353	fm=2.5
354	sfun=egamfmsfun+0.13-0.13*xp
355	c sigwlgam0=sqrt(egam)sfunsqwl
356	sigwlgam0=0.1/sqrt(egam)sfunsqwl
357	c write (,) ' fr sfun ',fr,sfun
358	c write (,) ' sigwlgam0 ',sigwlgam0,ewl,egam
359
360	endif
361	return
362	end
363
364	function dispwlm(ewl,m)
365
366	common /comsigma/sigmaph,sigmapd,sigphsq,sigpdsq
367
368
369	dispph=sigphsq*ewl
370	disppd=sigpdsq
371	dispwlm=dispph+disppd
372
373	return
374	end
375
376	real function zmidlshower(ee,m)
377	include 'comgam.for'
378	zmidlshower = 46.
379	c zmidlshower = 50.
380	return
381	end
382
383