[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONlib.cxx

#include "AliGenMUONlib.h"
#include "AliRun.h"
ClassImp(AliGenMUONlib)
//
//  Pions
Double_t AliGenMUONlib::PtPion(Double_t *px, Double_t *)
{
//
//     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
//     POWER LAW FOR PT > 500 MEV
//     MT SCALING BELOW (T=160 MEV)
//
  const Double_t p0 = 1.3;
  const Double_t xn = 8.28;
  const Double_t xlim=0.5;
  const Double_t t=0.160;
  const Double_t xmpi=0.139;
  const Double_t b=1.;
  Double_t y, y1, xmpi2, ynorm, a;
  Double_t x=*px;
  //
  y1=TMath::Power(p0/(p0+xlim),xn);
  xmpi2=xmpi*xmpi;
  ynorm=b*(TMath::Exp(-sqrt(xlim*xlim+xmpi2)/t));
  a=ynorm/y1;
  if (x > xlim)
    y=a*TMath::Power(p0/(p0+x),xn);
  else
    y=b*TMath::Exp(-sqrt(x*x+xmpi2)/t);
  return y*x;
}
//
// y-distribution
//
Double_t AliGenMUONlib::YPion( Double_t *py, Double_t *)
{
  const Double_t a    = 7000.;
  const Double_t dy   = 4.;

  Double_t y=TMath::Abs(*py);
  //
  Double_t ex = y*y/(2*dy*dy);
  return a*TMath::Exp(-ex);
}
//                 particle composition
//
Int_t AliGenMUONlib::IpPion()
{
    AliMC* pMC = AliMC::GetMC();
    Float_t random[1];
    pMC->Rndm(random,1);
    if (random[0] < 0.5) {
	return  211;
    } else {
	return -211;
    }
}

//____________________________________________________________
//
// Mt-scaling

Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
{
  //    SCALING EN MASSE PAR RAPPORT A PTPI
  //    MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
  const Double_t hm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
  //     VALUE MESON/PI AT 5 GEV
  const Double_t fmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
  np--;
  Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+hm[np]*hm[np])+2.0)),12.3);
  Double_t fmax2=f5/fmax[np];
  // PIONS
  Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
  Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
				 (sqrt(pt*pt+hm[np]*hm[np])+2.0)),12.3)/ fmax2;
  return fmtscal*ptpion;
}
//
// kaon
//
//                pt-distribution
//____________________________________________________________
Double_t AliGenMUONlib::PtKaon( Double_t *px, Double_t *)
{
  return PtScal(*px,2);
}

// y-distribution
//____________________________________________________________
Double_t AliGenMUONlib::YKaon( Double_t *py, Double_t *)
{
  const Double_t a    = 1000.;
  const Double_t dy   = 4.;
  

  Double_t y=TMath::Abs(*py);
  //
  Double_t ex = y*y/(2*dy*dy);
  return a*TMath::Exp(-ex);
}

//                 particle composition
//
Int_t AliGenMUONlib::IpKaon()
{
    AliMC* pMC = AliMC::GetMC();
    Float_t random[1];
    pMC->Rndm(random,1);
    if (random[0] < 0.5) {
	return  321;
    } else {
	return -321;
    }
}

//                    J/Psi 
//
//
//                pt-distribution
//____________________________________________________________
Double_t AliGenMUONlib::PtJpsi( Double_t *px, Double_t *)
{
  const Double_t pt0 = 4.;
  const Double_t xn  = 3.6;
  Double_t x=*px;
  //
  Double_t pass1 = 1.+(x/pt0)*(x/pt0);
  return x/TMath::Power(pass1,xn);
}
//
//               y-distribution
//____________________________________________________________
Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t *)
{
  const Double_t y0 = 4.;
  const Double_t b=1.;
  Double_t yj;
  Double_t y=TMath::Abs(*py);
  //
  if (y < y0)
    yj=b;
  else
    yj=b*TMath::Exp(-(y-y0)*(y-y0)/2);
  return yj;
}
//                 particle composition
//
Int_t AliGenMUONlib::IpJpsi()
{
    return 443;
}

//                      Upsilon
//
//
//                  pt-distribution
//____________________________________________________________
Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t * )
{
  const Double_t pt0 = 5.3;
  const Double_t xn  = 2.5;
  Double_t x=*px;
  //
  Double_t pass1 = 1.+(x/pt0)*(x/pt0);
  return x/TMath::Power(pass1,xn);
}
//
//                    y-distribution
//
//____________________________________________________________
Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t *)
{
  const Double_t y0 = 3.;
  const Double_t b=1.;
  Double_t yu;
  Double_t y=TMath::Abs(*py);
  //
  if (y < y0)
    yu=b;
  else
    yu=b*TMath::Exp(-(y-y0)*(y-y0)/2);
  return yu;
}
//                 particle composition
//
Int_t AliGenMUONlib::IpUpsilon()
{
    return 553;
}

//
//                        Phi
//
//
//    pt-distribution (by scaling of pion distribution)
//____________________________________________________________
Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t *)
{
  return PtScal(*px,7);
}
//    y-distribution
Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t *)
{
    Double_t *dummy=0;
    return YJpsi(px,dummy);
}
//                 particle composition
//
Int_t AliGenMUONlib::IpPhi()
{
    return 41;
}

//
//                        Charm
//
//
//                    pt-distribution
//____________________________________________________________
Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t *)
{
  const Double_t pt0 = 4.08;
  const Double_t xn  = 9.40;
  Double_t x=*px;
  //
  Double_t pass1 = 1.+(x/pt0)*(x/pt0);
  return x/TMath::Power(pass1,xn);
}
//                  y-distribution
Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t *)
{
    Double_t *dummy=0;
    return YJpsi(px,dummy);
}

Int_t AliGenMUONlib::IpCharm()
{  
    AliMC* pMC = AliMC::GetMC();
    Float_t random[2];
    Int_t ip;
//    411,421,431,4122
    pMC->Rndm(random,2);
    if (random[0] < 0.5) {
	ip=411;
    } else if (random[0] < 0.75) {
	ip=421;
    } else if (random[0] < 0.90) {
	ip=431;
    } else {
	ip=4122;
    }
    if (random[1] < 0.5) {ip=-ip;}
    
    return ip;
}


//
//                        Beauty
//
//
//                    pt-distribution
//____________________________________________________________
Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t *)
{
  const Double_t pt0 = 4.;
  const Double_t xn  = 3.6;
  Double_t x=*px;
  //
  Double_t pass1 = 1.+(x/pt0)*(x/pt0);
  return x/TMath::Power(pass1,xn);
}
//                     y-distribution
Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t *)
{
    Double_t *dummy=0;
    return YJpsi(px,dummy);
}

Int_t AliGenMUONlib::IpBeauty()
{  
    AliMC* pMC = AliMC::GetMC();
    Float_t random[2];
    Int_t ip;
    pMC->Rndm(random,2);
    if (random[0] < 0.5) {
	ip=511;
    } else if (random[0] < 0.75) {
	ip=521;
    } else if (random[0] < 0.90) {
	ip=531;
    } else {
	ip=5122;
    }
    if (random[1] < 0.5) {ip=-ip;}
    
    return ip;
}

typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
GenFunc AliGenMUONlib::GetPt(Param_t param)
{
    GenFunc func;
    switch (param) 
    {
    case phi_p:
	func=PtPhi;
	break;
    case jpsi_p:
	func=PtJpsi;
	break;
    case upsilon_p:
	func=PtUpsilon;
	break;
    case charm_p:
	func=PtCharm;
	break;
    case beauty_p:
	func=PtBeauty;
	break;
    case pion_p:
	func=PtPion;
	break;
    case kaon_p:
	func=PtKaon;
	break;
    default:
        func=0;
        printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
    }
    return func;
}

GenFunc AliGenMUONlib::GetY(Param_t param)
{
    GenFunc func;
    switch (param) 
    {
    case phi_p:
	func=YPhi;
	break;
    case jpsi_p:
	func=YJpsi;
	break;
    case upsilon_p:
	func=YUpsilon;
	break;
    case charm_p:
	func=YCharm;
	break;
    case beauty_p:
	func=YBeauty;
	break;
    case pion_p:
	func=YPion;
	break;
    case kaon_p:
	func=YKaon;
	break;
    default:
        func=0;
        printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
    }
    return func;
}
typedef Int_t (*GenFuncIp) ();
GenFuncIp AliGenMUONlib::GetIp(Param_t param)
{
    GenFuncIp func;
    switch (param) 
    {
    case phi_p:
	func=IpPhi;
	break;
    case jpsi_p:
	func=IpJpsi;
	break;
    case upsilon_p:
	func=IpUpsilon;
	break;
    case charm_p:
	func=IpCharm;
	break;
    case beauty_p:
	func=IpBeauty;
	break;
    case pion_p:
	func=IpPion;
	break;
    case kaon_p:
	func=IpKaon;
	break;
    default:
        func=0;
        printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
    }
    return func;
}
Commit	Line	Data
fe4da5cc	1	#include "AliGenMUONlib.h"
	2	#include "AliRun.h"
	3	ClassImp(AliGenMUONlib)
	4	//
	5	// Pions
	6	Double_t AliGenMUONlib::PtPion(Double_t px, Double_t )
	7	{
	8	//
	9	// PT-PARAMETERIZATION CDF, PRL 61(88) 1819
	10	// POWER LAW FOR PT > 500 MEV
	11	// MT SCALING BELOW (T=160 MEV)
	12	//
	13	const Double_t p0 = 1.3;
	14	const Double_t xn = 8.28;
	15	const Double_t xlim=0.5;
	16	const Double_t t=0.160;
	17	const Double_t xmpi=0.139;
	18	const Double_t b=1.;
	19	Double_t y, y1, xmpi2, ynorm, a;
	20	Double_t x=*px;
	21	//
	22	y1=TMath::Power(p0/(p0+xlim),xn);
	23	xmpi2=xmpi*xmpi;
	24	ynorm=b(TMath::Exp(-sqrt(xlimxlim+xmpi2)/t));
	25	a=ynorm/y1;
	26	if (x > xlim)
	27	y=a*TMath::Power(p0/(p0+x),xn);
	28	else
	29	y=bTMath::Exp(-sqrt(xx+xmpi2)/t);
	30	return y*x;
	31	}
753690b0	32	//
	33	// y-distribution
	34	//
	35	Double_t AliGenMUONlib::YPion( Double_t py, Double_t )
	36	{
	37	const Double_t a = 7000.;
	38	const Double_t dy = 4.;
	39
	40	Double_t y=TMath::Abs(*py);
	41	//
	42	Double_t ex = yy/(2dy*dy);
	43	return a*TMath::Exp(-ex);
	44	}
	45	// particle composition
	46	//
	47	Int_t AliGenMUONlib::IpPion()
	48	{
	49	AliMC* pMC = AliMC::GetMC();
	50	Float_t random[1];
	51	pMC->Rndm(random,1);
	52	if (random[0] < 0.5) {
	53	return 211;
	54	} else {
	55	return -211;
	56	}
	57	}
fe4da5cc	58
	59	//____________________________________________________________
	60	//
	61	// Mt-scaling
	62
	63	Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
	64	{
	65	// SCALING EN MASSE PAR RAPPORT A PTPI
	66	// MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
	67	const Double_t hm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
	68	// VALUE MESON/PI AT 5 GEV
	69	const Double_t fmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
	70	np--;
	71	Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+hm[np]*hm[np])+2.0)),12.3);
	72	Double_t fmax2=f5/fmax[np];
	73	// PIONS
	74	Double_t ptpion=100.PtPion(&pt, (Double_t) 0);
	75	Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
	76	(sqrt(ptpt+hm[np]hm[np])+2.0)),12.3)/ fmax2;
	77	return fmtscal*ptpion;
	78	}
	79	//
753690b0	80	// kaon
	81	//
	82	// pt-distribution
	83	//____________________________________________________________
	84	Double_t AliGenMUONlib::PtKaon( Double_t px, Double_t )
	85	{
	86	return PtScal(*px,2);
	87	}
	88
	89	// y-distribution
fe4da5cc	90	//____________________________________________________________
753690b0	91	Double_t AliGenMUONlib::YKaon( Double_t py, Double_t )
fe4da5cc	92	{
753690b0	93	const Double_t a = 1000.;
	94	const Double_t dy = 4.;
	95
	96
fe4da5cc	97	Double_t y=TMath::Abs(*py);
fe4da5cc	98	//
753690b0	99	Double_t ex = yy/(2dy*dy);
	100	return a*TMath::Exp(-ex);
	101	}
	102
	103	// particle composition
	104	//
	105	Int_t AliGenMUONlib::IpKaon()
	106	{
	107	AliMC* pMC = AliMC::GetMC();
	108	Float_t random[1];
	109	pMC->Rndm(random,1);
	110	if (random[0] < 0.5) {
	111	return 321;
	112	} else {
	113	return -321;
	114	}
fe4da5cc	115	}
753690b0	116
fe4da5cc	117	// J/Psi
	118	//
	119	//
	120	// pt-distribution
	121	//____________________________________________________________
	122	Double_t AliGenMUONlib::PtJpsi( Double_t px, Double_t )
	123	{
	124	const Double_t pt0 = 4.;
	125	const Double_t xn = 3.6;
	126	Double_t x=*px;
	127	//
	128	Double_t pass1 = 1.+(x/pt0)*(x/pt0);
	129	return x/TMath::Power(pass1,xn);
	130	}
	131	//
	132	// y-distribution
	133	//____________________________________________________________
	134	Double_t AliGenMUONlib::YJpsi(Double_t py, Double_t )
	135	{
	136	const Double_t y0 = 4.;
	137	const Double_t b=1.;
	138	Double_t yj;
	139	Double_t y=TMath::Abs(*py);
	140	//
	141	if (y < y0)
	142	yj=b;
	143	else
	144	yj=bTMath::Exp(-(y-y0)(y-y0)/2);
	145	return yj;
	146	}
	147	// particle composition
	148	//
	149	Int_t AliGenMUONlib::IpJpsi()
	150	{
	151	return 443;
	152	}
	153
	154	// Upsilon
	155	//
	156	//
	157	// pt-distribution
	158	//____________________________________________________________
	159	Double_t AliGenMUONlib::PtUpsilon( Double_t px, Double_t )
	160	{
	161	const Double_t pt0 = 5.3;
	162	const Double_t xn = 2.5;
	163	Double_t x=*px;
	164	//
	165	Double_t pass1 = 1.+(x/pt0)*(x/pt0);
	166	return x/TMath::Power(pass1,xn);
	167	}
	168	//
	169	// y-distribution
	170	//
	171	//____________________________________________________________
	172	Double_t AliGenMUONlib::YUpsilon(Double_t py, Double_t )
	173	{
	174	const Double_t y0 = 3.;
	175	const Double_t b=1.;
	176	Double_t yu;
	177	Double_t y=TMath::Abs(*py);
	178	//
	179	if (y < y0)
	180	yu=b;
181	else
182	yu=bTMath::Exp(-(y-y0)(y-y0)/2);
183	return yu;
184	}
185	// particle composition
186	//
187	Int_t AliGenMUONlib::IpUpsilon()
188	{
189	return 553;
190	}
191
192	//
193	// Phi
194	//
195	//
196	// pt-distribution (by scaling of pion distribution)
197	//____________________________________________________________
198	Double_t AliGenMUONlib::PtPhi( Double_t px, Double_t )
199	{
200	return PtScal(*px,7);
201	}
202	// y-distribution
203	Double_t AliGenMUONlib::YPhi( Double_t px, Double_t )
204	{
ccdc344b	205	Double_t *dummy=0;
fe4da5cc	206	return YJpsi(px,dummy);
	207	}
	208	// particle composition
	209	//
	210	Int_t AliGenMUONlib::IpPhi()
	211	{
	212	return 41;
	213	}
	214
	215	//
	216	// Charm
	217	//
	218	//
	219	// pt-distribution
	220	//____________________________________________________________
	221	Double_t AliGenMUONlib::PtCharm( Double_t px, Double_t )
	222	{
	223	const Double_t pt0 = 4.08;
	224	const Double_t xn = 9.40;
	225	Double_t x=*px;
	226	//
	227	Double_t pass1 = 1.+(x/pt0)*(x/pt0);
	228	return x/TMath::Power(pass1,xn);
	229	}
	230	// y-distribution
	231	Double_t AliGenMUONlib::YCharm( Double_t px, Double_t )
	232	{
ccdc344b	233	Double_t *dummy=0;
fe4da5cc	234	return YJpsi(px,dummy);
	235	}
	236
	237	Int_t AliGenMUONlib::IpCharm()
	238	{
	239	AliMC* pMC = AliMC::GetMC();
	240	Float_t random[2];
	241	Int_t ip;
	242	// 411,421,431,4122
	243	pMC->Rndm(random,2);
	244	if (random[0] < 0.5) {
	245	ip=411;
	246	} else if (random[0] < 0.75) {
	247	ip=421;
	248	} else if (random[0] < 0.90) {
	249	ip=431;
	250	} else {
	251	ip=4122;
	252	}
	253	if (random[1] < 0.5) {ip=-ip;}
	254
	255	return ip;
	256	}
	257
	258
	259	//
	260	// Beauty
	261	//
	262	//
	263	// pt-distribution
	264	//____________________________________________________________
	265	Double_t AliGenMUONlib::PtBeauty( Double_t px, Double_t )
	266	{
	267	const Double_t pt0 = 4.;
	268	const Double_t xn = 3.6;
	269	Double_t x=*px;
	270	//
	271	Double_t pass1 = 1.+(x/pt0)*(x/pt0);
	272	return x/TMath::Power(pass1,xn);
	273	}
	274	// y-distribution
	275	Double_t AliGenMUONlib::YBeauty( Double_t px, Double_t )
	276	{
ccdc344b	277	Double_t *dummy=0;
fe4da5cc	278	return YJpsi(px,dummy);
	279	}
	280
	281	Int_t AliGenMUONlib::IpBeauty()
	282	{
	283	AliMC* pMC = AliMC::GetMC();
	284	Float_t random[2];
	285	Int_t ip;
	286	pMC->Rndm(random,2);
	287	if (random[0] < 0.5) {
	288	ip=511;
	289	} else if (random[0] < 0.75) {
	290	ip=521;
	291	} else if (random[0] < 0.90) {
	292	ip=531;
	293	} else {
	294	ip=5122;
	295	}
	296	if (random[1] < 0.5) {ip=-ip;}
	297
	298	return ip;
	299	}
	300
	301	typedef Double_t (GenFunc) (Double_t, Double_t*);
753690b0	302	GenFunc AliGenMUONlib::GetPt(Param_t param)
fe4da5cc	303	{
fe4da5cc	304	GenFunc func;
753690b0	305	switch (param)
fe4da5cc	306	{
753690b0	307	case phi_p:
fe4da5cc	308	func=PtPhi;
fe4da5cc	309	break;
753690b0	310	case jpsi_p:
fe4da5cc	311	func=PtJpsi;
fe4da5cc	312	break;
753690b0	313	case upsilon_p:
fe4da5cc	314	func=PtUpsilon;
fe4da5cc	315	break;
753690b0	316	case charm_p:
fe4da5cc	317	func=PtCharm;
fe4da5cc	318	break;
753690b0	319	case beauty_p:
fe4da5cc	320	func=PtBeauty;
fe4da5cc	321	break;
753690b0	322	case pion_p:
	323	func=PtPion;
	324	break;
	325	case kaon_p:
	326	func=PtKaon;
	327	break;
119b35c7	328	default:
	329	func=0;
	330	printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
fe4da5cc	331	}
	332	return func;
	333	}
	334
753690b0	335	GenFunc AliGenMUONlib::GetY(Param_t param)
fe4da5cc	336	{
fe4da5cc	337	GenFunc func;
753690b0	338	switch (param)
fe4da5cc	339	{
753690b0	340	case phi_p:
fe4da5cc	341	func=YPhi;
fe4da5cc	342	break;
753690b0	343	case jpsi_p:
fe4da5cc	344	func=YJpsi;
fe4da5cc	345	break;
753690b0	346	case upsilon_p:
fe4da5cc	347	func=YUpsilon;
fe4da5cc	348	break;
753690b0	349	case charm_p:
fe4da5cc	350	func=YCharm;
fe4da5cc	351	break;
753690b0	352	case beauty_p:
fe4da5cc	353	func=YBeauty;
fe4da5cc	354	break;
753690b0	355	case pion_p:
	356	func=YPion;
	357	break;
	358	case kaon_p:
	359	func=YKaon;
	360	break;
119b35c7	361	default:
	362	func=0;
	363	printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
fe4da5cc	364	}
	365	return func;
	366	}
	367	typedef Int_t (*GenFuncIp) ();
753690b0	368	GenFuncIp AliGenMUONlib::GetIp(Param_t param)
fe4da5cc	369	{
fe4da5cc	370	GenFuncIp func;
753690b0	371	switch (param)
fe4da5cc	372	{
753690b0	373	case phi_p:
fe4da5cc	374	func=IpPhi;
fe4da5cc	375	break;
753690b0	376	case jpsi_p:
fe4da5cc	377	func=IpJpsi;
fe4da5cc	378	break;
753690b0	379	case upsilon_p:
fe4da5cc	380	func=IpUpsilon;
fe4da5cc	381	break;
753690b0	382	case charm_p:
fe4da5cc	383	func=IpCharm;
fe4da5cc	384	break;
753690b0	385	case beauty_p:
fe4da5cc	386	func=IpBeauty;
fe4da5cc	387	break;
753690b0	388	case pion_p:
	389	func=IpPion;
	390	break;
	391	case kaon_p:
	392	func=IpKaon;
	393	break;
119b35c7	394	default:
	395	func=0;
	396	printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
fe4da5cc	397	}
	398	return func;
	399	}
	400
	401
753690b0	402
753690b0	403