Introduction of the Copyright and cvs Log
[u/mrichter/AliRoot.git] / TRD / AliTRDsim.cxx
CommitLineData
4c039060 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/*
17$Log$
18*/
19
99d5402e 20#include <stdlib.h>
21
22#include "TH1.h"
23#include "TRandom.h"
24#include "TMath.h"
25
26#include "AliTRDsim.h"
27#include "AliTRDconst.h"
28
29ClassImp(AliTRDsim)
30
31
32const Float_t kD1 = kPeThick / kRaFoils;
33const Float_t kD2 = kRaThick / (kRaFoils + 1);
34
35//Root specials, to be removed
36
37static TH1F *h100, *h101, *h102;
38
39AliTRDsim::AliTRDsim()
40{
41 fNj=200;
42 fIrst=1;
43}
44
45AliTRDsim::AliTRDsim(AliModule* mod, Int_t foil, Int_t gas)
46{
47 Float_t a1, z1, ro1, rad, abs;
48 Float_t a2, z2, ro2;
49 char * name[21];
50 mod->AliGetMaterial(foil, name, a1, z1, ro1, rad, abs);
51 mod->AliGetMaterial(gas, name, a2, z2, ro2, rad, abs);
52 fOmega1 = 28.8*TMath::Sqrt(ro1*z1/a1);
53 fOmega2 = 28.8*TMath::Sqrt(ro2*z2/a2);
54}
55
56void AliTRDsim::trd_sim()
57{
58
59 const Float_t amass[4] = { 5.11e-4,.13957,.4937,.10566 };
60 const Double_t of[4] = { 20.9,24.4,14.27,26.9 };
61 const Double_t og[4] = { .28,.7,.74,.74 };
62 Int_t ifl = 0;
63 Int_t ig = 1;
64 Int_t nev = 1000;
65 Double_t gamma = -10.;
66
67 /* Local variables */
68 static Float_t temp, pres;
69 static Int_t i, j;
70 static Float_t o, sigma[200];
71 static Float_t trEn[10];
72 static Double_t omega1, omega2;
73 static Float_t am;
74 static Int_t np;
75 static Int_t ipa;
76
77 /* ***********************************************************************
78 */
79 /* TRD simulation - multimodule (regular rad.) */
80 /* after: M. CASTELLANO et al., */
81 /* COMP. PHYS. COMM. 51 (1988) 431 + COMP. PHYS. COMM. 61 (1990) 395 */
82
83 /* 17.07.1998 - A.Andronic */
84 /* 08.12.1998 - simplified version */
85
86 ipa = 0;
87 /* that's electron */
88 am = amass[ipa];
89 omega1 = of[ifl];
90 /* plasma frequency: foil and gap */
91 omega2 = og[ig - 1];
92 if (gamma < -1e5) printf("*** Momentum steps !!! ***\n");
93 if (gamma < 0. && gamma >= -1e5) {
94 gamma = sqrt(gamma * gamma + am * am) / am;
95 printf("*** Gamma (electron) = %f\n",gamma);
96 }
97 temp = 20.;
98 pres = 1.;
99 fBin = 100. / fNj;
100 /* binsize */
101 fL = 1. - fBin / 2.;
102 fU = fL + 100.;
103 /* setting the stage ................................... */
104 for (j = 0; j < fNj; ++j) {
105 /* getting the sigma values - for fixed energy values */
106 o = fBin * j + 1.;
107 /* omega in keV */
108 /* abs. in rad. (1 foi */
109 sigma[j] = fsigmaRad(ifl, ig, o);
110 }
111 printf(" Working...\n");
112 /* sampling over some events ........................... */
113 for (i = 0; i < nev; ++i) {
114 xtr(gamma, omega1, omega2, ro1, ro2, sigma, np, trEn);
115 /* TR: n, E */
116 h101->Fill(np);
117 /* sample nTR distr. */
118 for (j = 0; j < np; ++j) {
119 h102->Fill(trEn[j], 1. / fBin);
120 /* sample the TR en. distr. */
121 }
122 }
123 /* ------------------------------------------------------------------- */
124 /* else !ns steps */
125 /* enddo !imod */
126 /* events */
127 h100->Draw();
128 h101->Draw();
129 h102->Draw();
130} /* trd_sim__ */
131
132void AliTRDsim::xtr(Double_t gamma, Double_t omega1, Double_t omega2, Double_t ro1,
133 Double_t ro2,
134 Float_t *sigmaRad, Int_t &np, Float_t *trEn)
135{
136 /* Initialized data */
137
138 static Double_t alfa = .0072973;
139 static Double_t pi = 3.14159265;
140
141 /* Local variables */
142 static Double_t conv, a;
143 static Int_t i, j;
144 static Float_t o, w[200], omega;
145 static Double_t tetan, stemp;
146 static Float_t om;
147 static Double_t sk;
148 static Float_t wn[200];
149 static Double_t cs1, cs2;
150 static Double_t ro11, ro22, aux;
151 static Float_t ntr;
152 static Double_t sum;
153
154 /************************************************************************
155 ******/
156 /* TR: number and energy distr. */
157
158 /* Function Body */
159 sk = kD2 / kD1;
160 /* -------------- starts with the TR spectrum ------------- */
161
162 stemp = 0.;
163 for (j = 0; j < fNj; ++j) {
164 /* TR spectrum */
165 omega = (fBin * j + 1.) * 1e3;
166 /* keV->eV */
167 cs1 = omega1 / omega;
168 cs2 = omega2 / omega;
169 ro11 = omega * kD1 * 2.5 * (1. / (gamma * gamma) + cs1*cs1);
170 ro22 = omega * kD1 * 2.5 * (1. / (gamma * gamma) + cs2*cs2);
171 sum = 0.;
172 for (i = 0; i < 10; ++i) {
173/* 30 - it matters a bit */
174 tetan = (pi * 2. * (i+1) - (ro11 + sk * ro22)) / (sk + 1.);
175 if (tetan < 0.) {
176 tetan = 0.;
177 }
178 aux = 1. / (ro11 + tetan) - 1. / (ro22 + tetan);
179 a = tetan * (aux * aux) * (1. - cos(ro11 + tetan));
180 sum += a;
181 }
182 o = omega * .001;
183 /* eV->keV */
184 conv = 1. - exp(-kRaFoils * sigmaRad[j]);
185 w[j] = alfa * 4. / (sigmaRad[j] * (sk + 1.)) * conv * sum;
186 /* dW/domega */
187 wn[j] = w[j] / o;
188 /* dN/domega */
189 stemp += wn[j];
190 if (fIrst == 1) {
191 h100->Fill(o, wn[j]);
192 /* double precision not accepted */
193 }
194 }
195 /* -------------- done with the spectrum ------------- */
196 /* j (omega spectrum) */
197 ntr = stemp * fBin;
198 /* <nTR> (binsize corr.) */
199 om = h100->GetMean();
200 /* <Etr> */
201 if (fIrst == 1) {
202 /* prints the production */
203 printf(" Produced TR - <n>, <E>: %5.2f %6.2f KeV\n",ntr,om);
204 fIrst = 0;
205 }
206 /* prob. distr. */
207 np = gRandom->Poisson(ntr);
208 /* Np TR photons Poiss distr. from mean */
209 for (j = 0; j < np; ++j) {
210 /* TR energy (binsize corr.) */
211 trEn[j] = hisran(wn, fNj, fL, fBin);
212 /* their energy */
213 }
214}
215
216Float_t AliTRDsim::fsigmaRad(Float_t ro1, Float_t ro2, Int_t ig, Float_t o)
217{
218
219 /* Local variables */
220 static Float_t pres;
221 static Double_t mumu;
222 static Int_t j;
223 static Double_t t;
224 static Int_t i1, i2;
225 static Double_t x1;
226 static Double_t mu1, mu2, deo, omf[36], omg[36], muf[36], mug[36];
227
228 static Bool_t first = kTRUE;
229
230 /* cccccccccccccccccccccccccccccccccccccccccccc */
231 /* calculates sigma for radiator - one foil+one gap */
232
233 if(first) {
234 FILE* inp = fopen("po.x","r");
235 for (j=0;j<36;++j) {
236 fscanf(inp,"%lf %lf %lf",&omf[j],&muf[j],&mumu);
237 }
238 fclose(inp);
239 inp = fopen("he.x","r");
240 for (j=0;j<36;++j) {
241 fscanf(inp,"%lf %lf %lf",&omg[j],&mug[j],&mumu);
242 }
243 fclose(inp);
244 first=kFALSE;
245 }
246 /* first */
247 x1 = o * .001;
248 /* keV->MeV */
249 if (x1 >= .001) {
250 locate(omf, 36, x1, i1, deo);
251 mu1 = muf[i1] - deo * (muf[i1] - muf[i1+1]) / (omf[i1+1] - omf[i1]);
252 locate(omg, 36, x1, i2, deo);
253 mu2 = mug[i2] - deo * (mug[i2] - mug[i2+1]) / (omg[i2+1] - omg[i2]);
254 t = 273.16;
255 /* gases at 0 C */
256 return (mu1*ro1*kD1+mu2*293.16/t * ro2*kD2)/1e4;
257 /* mu */
258 } else {
259 return 1e6;
260 }
261}
262
263Int_t AliTRDsim::locate(Double_t *xv, Int_t n, Double_t xval,
264 Int_t &kl, Double_t &dx)
265{
266 /* -------------------------------------------------------------- */
267 /* locates a point (xval) in a 1-dim grid (xv(n)) --> iloc,dx,ier */
268 /* -------------------------------------------------------------- */
269 /* Function Body */
270 if (xval >= xv[n-1]) return 1;
271 if (xval < xv[0]) return -1;
272 Int_t km,kh=n-1;
273 kl=0;
274 while(kh-kl>1) if(xval<xv[km=(kl+kh)/2]) kh=km; else kl=km;
275 if(xval<xv[kl] || xval > xv[kl+1] || kl >= n-1) {
276 printf("locate failed xv[%d] %f xval %f xv[%d] %f!!!\n",
277 kl,xv[kl],xval,kl+1,xv[kl+1]);
278 exit(1);
279 }
280 dx=xval-xv[kl];
281 return 0;
282}
283
284Float_t AliTRDsim::hisran(Float_t *y, Int_t n, Float_t xlo, Float_t xwid)
285{
286 /* Local variables */
287 Float_t yinv, ytot=0;
288 Int_t i;
289 Float_t yr;
290
291/* SUBROUTINE TO GENERATE RANDOM NUMBERS */
292/* ACCORDING TO AN EMPIRICAL DISTRIBUTION */
293/* SUPPLIED BY THE USER IN THE FORM OF A HISTOGRAM */
294/* F. JAMES, MAY, 1976 */
295
296 if (y[n-1] != 1.) {
297
298/* INITIALIZE HISTOGRAM TO FORM CUMULATIVE DISTRIBUTION */
299
300 ytot = 0.;
301 for (i = 0; i < n; ++i) {
302 if (y[i] < 0.) {
303 printf("hisran: found value y[%d] = %f\n",i,y[i]);
304 exit(1);
305 }
306 ytot += y[i];
307 y[i] = ytot;
308 }
309 if (ytot <= 0.) {
310 printf("hisran: total probability %f < 0\n",ytot);
311 exit(1);
312 }
313 yinv = 1. / ytot;
314 for (i = 0; i < n-1; ++i) {
315 y[i] *= yinv;
316 }
317 y[n-1] = 1.;
318 }
319/* NOW GENERATE RANDOM NUMBER BETWEEN 0 AND ONE */
320 yr = gRandom->Rndm();
321/* AND TRANSFORM IT INTO THE CORRESPONDING X-VALUE */
322 if(yr<=y[0]) return xlo + xwid * (yr / y[0]);
323 else {
324 Int_t km,kl=0,kh=n-1;
325 while(kh-kl>1) if(yr<y[km=(kl+kh)/2]) kh=km; else kl=km;
326 return xlo + xwid * (kl + (yr - y[kl]) / (y[kl + 1] - y[kl]));
327 }
328}
329