Introduction of the Copyright and cvs Log
[u/mrichter/AliRoot.git] / TRD / AliTRDsim.cxx
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
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
29 ClassImp(AliTRDsim)
30
31
32 const Float_t kD1 = kPeThick / kRaFoils;
33 const Float_t kD2 = kRaThick / (kRaFoils + 1);
34
35 //Root specials, to be removed
36
37 static TH1F *h100, *h101, *h102;
38
39 AliTRDsim::AliTRDsim()
40 {
41   fNj=200;
42   fIrst=1;
43 }
44
45 AliTRDsim::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
56 void 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
132 void 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
216 Float_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
263 Int_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
284 Float_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