[u/mrichter/AliRoot.git] / TRD / AliTRDsim.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/*
$Log$
*/

#include <stdlib.h>

#include "TH1.h"
#include "TRandom.h"
#include "TMath.h"

#include "AliTRDsim.h"
#include "AliTRDconst.h"

ClassImp(AliTRDsim)


const Float_t kD1 = kPeThick / kRaFoils;
const Float_t kD2 = kRaThick / (kRaFoils + 1);

//Root specials, to be removed

static TH1F *h100, *h101, *h102;

AliTRDsim::AliTRDsim()
{
  fNj=200;
  fIrst=1;
}

AliTRDsim::AliTRDsim(AliModule* mod, Int_t foil, Int_t gas)
{
  Float_t a1, z1, ro1, rad, abs;
  Float_t a2, z2, ro2;
  char * name[21];
  mod->AliGetMaterial(foil, name, a1, z1, ro1, rad, abs);
  mod->AliGetMaterial(gas, name, a2, z2, ro2, rad, abs);
  fOmega1 = 28.8*TMath::Sqrt(ro1*z1/a1);
  fOmega2 = 28.8*TMath::Sqrt(ro2*z2/a2);
}

void AliTRDsim::trd_sim()
{

  const Float_t amass[4] = { 5.11e-4,.13957,.4937,.10566 };
  const Double_t of[4] = { 20.9,24.4,14.27,26.9 };
  const Double_t og[4] = { .28,.7,.74,.74 };
  Int_t ifl = 0;
  Int_t ig = 1;
  Int_t nev = 1000;
  Double_t gamma = -10.;
  
  /* Local variables */
  static Float_t temp, pres;
  static Int_t i, j;
  static Float_t o, sigma[200];
  static Float_t trEn[10];
  static Double_t omega1, omega2;
  static Float_t am;
  static Int_t np;
  static Int_t ipa;
  
  /* ***********************************************************************
   */
  /*  TRD simulation - multimodule (regular rad.) */
  /*     after: M. CASTELLANO et al., */
  /*   COMP. PHYS. COMM. 51 (1988) 431 + COMP. PHYS. COMM. 61 (1990) 395 */
  
  /*   17.07.1998 - A.Andronic */
  /*   08.12.1998 - simplified version */
  
  ipa = 0;
  /* that's electron */
  am = amass[ipa];
  omega1 = of[ifl];
  /* plasma frequency: foil and gap */
  omega2 = og[ig - 1];
  if (gamma < -1e5) printf("*** Momentum steps !!! ***\n");
  if (gamma < 0. && gamma >= -1e5) {
    gamma = sqrt(gamma * gamma + am * am) / am;
    printf("*** Gamma (electron) = %f\n",gamma);
  }
  temp = 20.;
  pres = 1.;
  fBin = 100. /  fNj;
  /* binsize */
  fL = 1. - fBin / 2.;
  fU = fL + 100.;
  /*  setting the stage ................................... */
  for (j = 0; j < fNj; ++j) {
    /* getting the sigma values - for fixed energy values */
    o = fBin *  j + 1.;
    /* omega in keV */
    /* abs. in rad. (1 foi */
    sigma[j] = fsigmaRad(ifl, ig, o);
  }
  printf(" Working...\n");
  /*  sampling over some events ........................... */
  for (i = 0; i < nev; ++i) {
    xtr(gamma, omega1, omega2, ro1, ro2, sigma, np, trEn);
    /* TR: n, E */
    h101->Fill(np);
    /* sample nTR distr. */
    for (j = 0; j < np; ++j) {
      h102->Fill(trEn[j], 1. / fBin);
      /* sample the TR en. distr. */
    }
  }
  /* ------------------------------------------------------------------- */
  /*      else  !ns steps */
  /*      enddo  !imod */
  /* events */
  h100->Draw();
  h101->Draw();
  h102->Draw();
} /* trd_sim__ */

void AliTRDsim::xtr(Double_t gamma, Double_t omega1, Double_t omega2, Double_t ro1,
		    Double_t ro2,
		     Float_t *sigmaRad, Int_t &np, Float_t *trEn)
{
    /* Initialized data */

  static Double_t alfa = .0072973;
  static Double_t pi = 3.14159265;
  
  /* Local variables */
  static Double_t conv, a;
  static Int_t i, j;
  static Float_t o, w[200], omega;
  static Double_t tetan, stemp;
  static Float_t om;
  static Double_t sk;
  static Float_t wn[200];
  static Double_t cs1, cs2;
  static Double_t ro11, ro22, aux;
  static Float_t ntr;
  static Double_t sum;
  
  /************************************************************************
   ******/
  /*   TR: number and energy distr. */
  
  /* Function Body */
  sk = kD2 / kD1;
  /* -------------- starts with the TR spectrum ------------- */
  
  stemp = 0.;
  for (j = 0; j < fNj; ++j) {
    /* TR spectrum */
    omega = (fBin *  j + 1.) * 1e3;
    /* keV->eV */
    cs1 = omega1 / omega;
    cs2 = omega2 / omega;
    ro11 = omega * kD1 * 2.5 * (1. / (gamma * gamma) + cs1*cs1);
    ro22 = omega * kD1 * 2.5 * (1. / (gamma * gamma) + cs2*cs2);
    sum = 0.;
    for (i = 0; i < 10; ++i) {
/* 30 - it matters a bit */
      tetan = (pi * 2. *  (i+1) - (ro11 + sk * ro22)) / (sk + 1.);
      if (tetan < 0.) {
	tetan = 0.;
      }
      aux = 1. / (ro11 + tetan) - 1. / (ro22 + tetan);
      a = tetan * (aux * aux) * (1. - cos(ro11 + tetan));
      sum += a;
    }
    o = omega * .001;
    /* eV->keV */
    conv = 1. - exp(-kRaFoils * sigmaRad[j]);
    w[j] = alfa * 4. / (sigmaRad[j] * (sk + 1.)) * conv * sum;
    /* dW/domega */
    wn[j] = w[j] / o;
    /* dN/domega */
    stemp += wn[j];
    if (fIrst == 1) {
      h100->Fill(o, wn[j]);
      /* double precision not accepted */
    }
  }
  /* -------------- done with the spectrum ------------- */
  /* j (omega spectrum) */
  ntr = stemp * fBin;
  /* <nTR> (binsize corr.) */
  om = h100->GetMean();
  /* <Etr> */
  if (fIrst == 1) {
    /* prints the production */
    printf(" Produced TR - <n>, <E>: %5.2f  %6.2f  KeV\n",ntr,om);
    fIrst = 0;
  }
  /* prob. distr. */
  np = gRandom->Poisson(ntr);
  /* Np TR photons Poiss distr. from mean */
  for (j = 0; j < np; ++j) {
    /* TR energy (binsize corr.) */
    trEn[j] = hisran(wn, fNj, fL, fBin);
    /* their energy */
  }
}

Float_t AliTRDsim::fsigmaRad(Float_t ro1, Float_t ro2, Int_t ig, Float_t o)
{

  /* Local variables */
  static Float_t pres;
  static Double_t mumu;
  static Int_t j;
  static Double_t t;
  static Int_t i1, i2;
  static Double_t x1;
  static Double_t mu1, mu2, deo, omf[36], omg[36], muf[36], mug[36];
  
  static Bool_t first = kTRUE;
  
  /* cccccccccccccccccccccccccccccccccccccccccccc */
  /*    calculates sigma for radiator - one foil+one gap */
  
  if(first) {
    FILE* inp = fopen("po.x","r");
    for (j=0;j<36;++j) {
      fscanf(inp,"%lf %lf %lf",&omf[j],&muf[j],&mumu);
    }
    fclose(inp);
    inp = fopen("he.x","r");
    for (j=0;j<36;++j) {
      fscanf(inp,"%lf %lf %lf",&omg[j],&mug[j],&mumu);
    }
    fclose(inp);
    first=kFALSE;
  }
  /* first */
  x1 = o * .001;
  /* keV->MeV */
  if (x1 >= .001) {
    locate(omf, 36, x1, i1, deo);
    mu1 = muf[i1] - deo * (muf[i1] - muf[i1+1]) / (omf[i1+1] - omf[i1]);
    locate(omg, 36, x1, i2, deo);
    mu2 = mug[i2] - deo * (mug[i2] - mug[i2+1]) / (omg[i2+1] - omg[i2]);
    t = 273.16;
    /* gases at 0 C */
    return (mu1*ro1*kD1+mu2*293.16/t * ro2*kD2)/1e4;
    /* mu */
  } else {
    return 1e6;
  }
} 

Int_t AliTRDsim::locate(Double_t *xv, Int_t n, Double_t xval, 
		       Int_t &kl, Double_t &dx)
{
  /* -------------------------------------------------------------- */
  /*  locates a point (xval) in a 1-dim grid (xv(n)) --> iloc,dx,ier */
  /* -------------------------------------------------------------- */
  /* Function Body */
  if (xval >= xv[n-1]) return 1;
  if (xval < xv[0]) return -1;
  Int_t km,kh=n-1;
  kl=0;
  while(kh-kl>1) if(xval<xv[km=(kl+kh)/2]) kh=km; else kl=km;
  if(xval<xv[kl] || xval > xv[kl+1] || kl >= n-1) {
    printf("locate failed xv[%d] %f xval %f xv[%d] %f!!!\n",
	   kl,xv[kl],xval,kl+1,xv[kl+1]);
    exit(1);
  }
  dx=xval-xv[kl];
  return 0;
}

Float_t AliTRDsim::hisran(Float_t *y, Int_t n, Float_t xlo, Float_t xwid)
{
    /* Local variables */
    Float_t yinv, ytot=0;
    Int_t i;
    Float_t yr;

/*         SUBROUTINE TO GENERATE RANDOM NUMBERS */
/*         ACCORDING TO AN EMPIRICAL DISTRIBUTION */
/*         SUPPLIED BY THE USER IN THE FORM OF A HISTOGRAM */
/*         F. JAMES,    MAY, 1976 */

    if (y[n-1] != 1.) {

/*         INITIALIZE HISTOGRAM TO FORM CUMULATIVE DISTRIBUTION */

      ytot = 0.;
      for (i = 0; i < n; ++i) {
	if (y[i] < 0.) {
	  printf("hisran: found value y[%d] = %f\n",i,y[i]);
	  exit(1);
	}
	ytot += y[i];
	y[i] = ytot;
      }
      if (ytot <= 0.) {
	printf("hisran: total probability %f < 0\n",ytot);
	exit(1);
      }
      yinv = 1. / ytot;
      for (i = 0; i < n-1; ++i) {
	y[i] *= yinv;
      }
      y[n-1] = 1.;
    }
/*         NOW GENERATE RANDOM NUMBER BETWEEN 0 AND ONE */
    yr = gRandom->Rndm();
/*         AND TRANSFORM IT INTO THE CORRESPONDING X-VALUE */
    if(yr<=y[0]) return xlo + xwid * (yr / y[0]);
    else {
      Int_t km,kl=0,kh=n-1;
      while(kh-kl>1) if(yr<y[km=(kl+kh)/2]) kh=km; else kl=km;
      return xlo + xwid * (kl + (yr - y[kl]) / (y[kl + 1] - y[kl]));
    }
}
Commit	Line	Data
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
	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.8TMath::Sqrt(ro1z1/a1);
	53	fOmega2 = 28.8TMath::Sqrt(ro2z2/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 (mu1ro1kD1+mu2293.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