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