[u/mrichter/AliRoot.git] / HBTAN / AliHBTCrab.cxx

#include "AliHBTCrab.h"
//______________________________________________________________________
/////////////////////////////////////////////////////////////////////////
//                                                                     //
// AliRoot wrapper to CRAB                                             //
// taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html  //
// written by Scott Pratt                                              //
//                                                                     //
/////////////////////////////////////////////////////////////////////////

#include "AliHBTPair.h"

#include <TMath.h>
#include <TPDGCode.h>

#include "volya_complex.h"

//#include <complex>
//using namespace std;
AliHBTCrab* AliHBTCrab::fgCrab = 0x0;

const Double_t AliHBTCrab::fgkWcons = 1./0.1973;
const Double_t AliHBTCrab::fgkROOT2=1.41421356237309504880;
const double_complex AliHBTCrab::ci(0.0,1.0);

/************************************************************/

AliHBTCrab* AliHBTCrab::Instance()
{
// returns instance of class
 if (fgCrab == 0x0)
  {
    fgCrab = new AliHBTCrab();
  }
 return fgCrab;
}
//===================================================================
void AliHBTCrab::Set()
{
 //sets this as weighitng class
 Info("Set","Setting CRAB as Weighing Class");
 
 if ( fgWeights == 0x0 )  
  {
    fgWeights = AliHBTCrab::Instance();
    return;
  }  
 if ( fgWeights == AliHBTCrab::Instance() ) return;
 delete fgWeights;
 fgWeights = AliHBTCrab::Instance();
}
//===================================================================

AliHBTCrab::AliHBTCrab():
fBreitWigner(kFALSE),
fReducedMom(kTRUE),
fMaxMomentum(100.0)
{
  //ctor
}
//===================================================================
void AliHBTCrab::Init(Int_t pid1,Int_t pid2)
{
  MASS1 = TDatabasePDG::Instance()->GetParticle(pid1)->Mass();
  MASS2 = TDatabasePDG::Instance()->GetParticle(pid2)->Mass();
  INTERACTION_WSYM = 1.0;
  INTERACTION_WANTI = 0.0;
  INTERACTION_WNOSYM = 0.0;
  INTERACTION_DELK  = 1.0;
  INTERACTION_NKMAX = 100;
     
  fPid1 = pid1;   
  fPid2 = pid2;
  
}

Bool_t AliHBTCrab::SetConfig(const AliHBTPair* pair)
{
  Int_t pdg1 = pair->Particle1()->GetPdgCode();
  Int_t pdg2 = pair->Particle2()->GetPdgCode();
  
  if ( ( pdg1 == fPid1) && ( pdg2  == fPid2) ) return kFALSE;
  else Init (pdg1,pdg2);
  
  return kTRUE;
}
//===================================================================

Double_t AliHBTCrab::GetWeight(const AliHBTPair* partpair)
{
  Double_t qred, r, qdotr, mom;
  Int_t test;
  
  SetConfig(partpair);
  
  get_com_quantities(partpair, &qred, &r, &qdotr, &mom, &test);
  if(test==0) 
   {
     Info("GetWeight","Test is 0");
   }
  Double_t corr = corrcalc(qred,qdotr,r);
  
  return corr;
}
//===================================================================

void AliHBTCrab::get_com_quantities(const AliHBTPair* pair, 
       double *qred,double *r,double *qdotr,double *mom, int *test)
 {
//************************************
//  ALICE //

 double p1[4]; 
 double p2[4]; 
 double r1[4]; 
 double r2[4]; 
 
 static const Double_t cmtofm = 1.e13;
 static const Double_t cmtoOneOverGeV = cmtofm*fgkWcons;
 
 AliHBTParticle *part1 = pair->Particle1();
 AliHBTParticle *part2 = pair->Particle2();

 p1[0] = part1->Energy()*1000.0;
 p1[1] = part1->Px()*1000.0;
 p1[2] = part1->Py()*1000.0;
 p1[3] = part1->Pz()*1000.0;
 
 p2[0] = part2->Energy()*1000.0;
 p2[1] = part2->Px()*1000.0;
 p2[2] = part2->Py()*1000.0;
 p2[3] = part2->Pz()*1000.0;

 r1[0] = part1->T();
 r1[1] = part1->Vx()*cmtofm;
 r1[2] = part1->Vy()*cmtofm;
 r1[3] = part1->Vz()*cmtofm;

 r2[0] = part2->T();
 r2[1] = part2->Vx()*cmtofm;
 r2[2] = part2->Vy()*cmtofm;
 r2[3] = part2->Vz()*cmtofm;

//  END OF ALICE STUFF

// This code is written by Scott Pratt
// taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html
  int alpha;
  double kdotr;
  if (fReducedMom)
   {
     const double momtest=4.0*fMaxMomentum*fMaxMomentum;
   }  
  else
   {
     const double momtest=fMaxMomentum*fMaxMomentum;
   }  
   
  double ptot2,pdotr,pp,rr;
  
  if ( part1->GetPdgCode() == part2->GetPdgCode() )
   {
    *test=1;
    *mom=-(p2[0]-p1[0])*(p2[0]-p1[0]);
    for(alpha=1;alpha<4;alpha++){
      *mom=*mom+(p2[alpha]-p1[alpha])*(p2[alpha]-p1[alpha]);
    }
   //#if ! defined MIXED_PAIRS_FOR_DENOM
   //  if(*mom>momtest){
   //    *test=0;
   //    return;
   //  }
   //#endif
    pp=(p1[0]+p2[0]);
    rr=(r2[0]-r1[0]);
    pdotr=pp*rr;
    kdotr=(p2[0]-p1[0])*rr;
    ptot2=pp*pp;
    *r=-rr*rr;
    for(alpha=1;alpha<4;alpha++){
      pp=(p1[alpha]+p2[alpha]);
      rr=(r2[alpha]-r1[alpha]);
      pdotr=pdotr-pp*rr;
      kdotr=kdotr-(p2[alpha]-p1[alpha])*rr;
      ptot2=ptot2-pp*pp;
      *r=*r+rr*rr;
    }
    *mom=sqrt(*mom);
    *qred=0.5**mom;

    if (fReducedMom)
     {
      *mom=*qred;
     } 

    *qdotr=0.5*kdotr;
    *r=sqrt(*r+pdotr*pdotr/ptot2);
   }  
  else //identical
  {
  //  const double  kdotp=MASS2*MASS2-MASS1*MASS1;
   const double  kdotp = part2->GetMass()*part2->GetMass()- part1->GetMass()*part1->GetMass();
   *test=1;
   *mom=-(p2[0]-p1[0])*(p2[0]-p1[0]);
   ptot2=(p1[0]+p2[0])*(p1[0]+p2[0]);
   for(alpha=1;alpha<4;alpha++){
     *mom=*mom+(p2[alpha]-p1[alpha])*(p2[alpha]-p1[alpha]);
     ptot2=ptot2-(p1[alpha]+p2[alpha])*(p1[alpha]+p2[alpha]);
   }
   *mom=*mom+kdotp*kdotp/ptot2;
  //#if ! defined MIXED_PAIRS_FOR_DENOM
  //  if(*mom>momtest){
  //    *test=0;
  //    return;
  //  }
  //#endif
   pp=(p1[0]+p2[0]);
   rr=(r2[0]-r1[0]);
   pdotr=pp*rr;
   kdotr=(p2[0]-p1[0])*rr;
   *r=-rr*rr;
   for(alpha=1;alpha<4;alpha++){
     pp=(p1[alpha]+p2[alpha]);
     rr=(r2[alpha]-r1[alpha]);
     pdotr=pdotr-pp*rr;
     kdotr=kdotr-(p2[alpha]-p1[alpha])*rr;
     *r=*r+rr*rr;
   }
   kdotr=(-kdotr+kdotp*pdotr/ptot2);
   *mom=sqrt(*mom);
   *qred=0.5**mom;

   if (fReducedMom)
    {
     *mom=*qred;
    } 

   *qdotr=0.5*kdotr;
   *r=sqrt(*r+pdotr*pdotr/ptot2);
  }//not identical  

  return;
}


//===================================================================

double  AliHBTCrab::corrcalc(double trueqred,double trueqdotr,double truer)
{
//#define REDUCED_MOM
// This code is written by Scott Pratt
// taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html
  double eta,arg,krmax,corr0;
  double xx,xxprime,xxjj,p1,zk;
  int jj,kk,ipart,ipartcount,ispin;
  double pfactor,wsym_leftover,wanti_leftover,wnosym_leftover;
  double qred,qdotr,r;
  const double rmass=MASS1*MASS2/(MASS1+MASS2);
  double_complex cphi1,cphi2,cphis,cphia;

  arg=trueqdotr/197.323-2.0*TMath::Pi()*floor(trueqdotr/(197.323*2.0*TMath::Pi()));
  cphi1=exp(ci*arg);
  cphis=fgkROOT2*real(cphi1);
  cphia=ci*fgkROOT2*imag(cphi1);
  corr0=real(INTERACTION_WSYM*cphis*conj(cphis)
	     +INTERACTION_WANTI*cphia*conj(cphia)
	     +INTERACTION_WNOSYM*cphi1*conj(cphi1));
  goto OUTSIDE_INTERACTION_RANGE;

#ifdef REDUCED_MOM
  kk=(int)floor(trueqred/INTERACTION_DELK);
  qred=(0.5+kk)*INTERACTION_DELK;
#else
  kk=(int)floor(2.0*trueqred/INTERACTION_DELK);
  qred=(0.5+kk)*INTERACTION_DELK/2.0;
#endif
  qdotr=trueqdotr*qred/trueqred;
  if(kk>=INTERACTION_NKMAX){
    corr0=1.0;
    goto OUTSIDE_INTERACTION_RANGE;
  }
  r=truer;

  eta=0.0;
  arg=qdotr/197.323-2.0*TMath::Pi()*floor(qdotr/(197.323*2.0*TMath::Pi()));
  cphi1=exp(ci*arg);
  cphi2=conj(cphi1);

  cphis=(cphi1+cphi2)/fgkROOT2;
  cphia=(cphi1-cphi2)/fgkROOT2;
  corr0=0.0;
  /* If there are corrections for strong interactions, add the
     change for each partial wave.  If npartial = 0 then there
     are no strong int. corrections. */
  wsym_leftover=INTERACTION_WSYM;
  wanti_leftover=INTERACTION_WANTI;
  wnosym_leftover=INTERACTION_WNOSYM;

  corr0=corr0+real(wsym_leftover*cphis*conj(cphis)
		   +wanti_leftover*cphia*conj(cphia)
		   +wnosym_leftover*cphi1*conj(cphi1));
OUTSIDE_INTERACTION_RANGE:

#ifdef BREIT_WIGNER
  corr0=corr0+bwcalc(trueqred,truer);
#endif

  return corr0;
}

double_complex AliHBTCrab::cgamma(double_complex c){
  /* This calc.s gamma functions which are in the form gamma(n+i*y)
     where n is an int and y is real. */
// This code is written by Scott Pratt
// taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html
  double_complex cg,cphase;
  int mm,j;
  double x,y,phase,delp,cgmag;
  x=real(c);
  y=imag(c);
  phase=-TMath::E()*y;
  for(j=1;j<=100000;j++){
    delp=(y/(double)j)-atan(y/(double)j);
    phase=phase+delp;
    if(fabs(delp)<1E-10) goto CGAMMA_ESCAPE;
  }
  printf("oops not accurate enough, increase jmax\n");
CGAMMA_ESCAPE:
  phase=phase-2.0*TMath::Pi()*floor(phase/(2.0*TMath::Pi()));
  cphase=exp(ci*phase);
  cgmag=sqrt(TMath::Pi()*y/sinh(TMath::Pi()*y));
  mm=(int)floor(x+0.5);
  cg=cgmag*cphase;
  if(mm<1){
    for(j=1;j<=-mm+1;j++){
      cg=cg/(1.0+(double)(-j)+ci*y);
    }
  }
  if(mm>1) {
    for(j=1;j<=mm-1;j++){
      cg=cg*((double)(j)+ci*y);
    }
  }
  return cg;
}
//===================================================================
Commit	Line	Data
dd82cadc	1	#include "AliHBTCrab.h"
	2	//______________________________________________________________________
	3	/////////////////////////////////////////////////////////////////////////
	4	// //
	5	// AliRoot wrapper to CRAB //
	6	// taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html //
	7	// written by Scott Pratt //
	8	// //
	9	/////////////////////////////////////////////////////////////////////////
	10
	11	#include "AliHBTPair.h"
	12
	13	#include <TMath.h>
	14	#include <TPDGCode.h>
	15
	16	#include "volya_complex.h"
	17
	18	//#include <complex>
	19	//using namespace std;
	20	AliHBTCrab* AliHBTCrab::fgCrab = 0x0;
	21
	22	const Double_t AliHBTCrab::fgkWcons = 1./0.1973;
	23	const Double_t AliHBTCrab::fgkROOT2=1.41421356237309504880;
	24	const double_complex AliHBTCrab::ci(0.0,1.0);
	25
	26	/************************************************************/
	27
	28	AliHBTCrab* AliHBTCrab::Instance()
	29	{
	30	// returns instance of class
	31	if (fgCrab == 0x0)
	32	{
	33	fgCrab = new AliHBTCrab();
	34	}
	35	return fgCrab;
	36	}
	37	//===================================================================
	38	void AliHBTCrab::Set()
	39	{
	40	//sets this as weighitng class
	41	Info("Set","Setting CRAB as Weighing Class");
	42
	43	if ( fgWeights == 0x0 )
	44	{
	45	fgWeights = AliHBTCrab::Instance();
	46	return;
	47	}
	48	if ( fgWeights == AliHBTCrab::Instance() ) return;
	49	delete fgWeights;
	50	fgWeights = AliHBTCrab::Instance();
	51	}
	52	//===================================================================
	53
	54	AliHBTCrab::AliHBTCrab():
	55	fBreitWigner(kFALSE),
	56	fReducedMom(kTRUE),
	57	fMaxMomentum(100.0)
	58	{
	59	//ctor
	60	}
	61	//===================================================================
	62	void AliHBTCrab::Init(Int_t pid1,Int_t pid2)
	63	{
	64	MASS1 = TDatabasePDG::Instance()->GetParticle(pid1)->Mass();
65	MASS2 = TDatabasePDG::Instance()->GetParticle(pid2)->Mass();
66	INTERACTION_WSYM = 1.0;
67	INTERACTION_WANTI = 0.0;
68	INTERACTION_WNOSYM = 0.0;
69	INTERACTION_DELK = 1.0;
70	INTERACTION_NKMAX = 100;
71
72	fPid1 = pid1;
73	fPid2 = pid2;
74
75	}
76
77	Bool_t AliHBTCrab::SetConfig(const AliHBTPair* pair)
78	{
79	Int_t pdg1 = pair->Particle1()->GetPdgCode();
80	Int_t pdg2 = pair->Particle2()->GetPdgCode();
81
82	if ( ( pdg1 == fPid1) && ( pdg2 == fPid2) ) return kFALSE;
83	else Init (pdg1,pdg2);
84
85	return kTRUE;
86	}
87	//===================================================================
88
89	Double_t AliHBTCrab::GetWeight(const AliHBTPair* partpair)
90	{
91	Double_t qred, r, qdotr, mom;
92	Int_t test;
93
94	SetConfig(partpair);
95
96	get_com_quantities(partpair, &qred, &r, &qdotr, &mom, &test);
97	if(test==0)
98	{
99	Info("GetWeight","Test is 0");
100	}
101	Double_t corr = corrcalc(qred,qdotr,r);
102
103	return corr;
104	}
105	//===================================================================
106
107	void AliHBTCrab::get_com_quantities(const AliHBTPair* pair,
108	double qred,double r,double qdotr,double mom, int *test)
109	{
110	//************************************
111	// ALICE //
112
113	double p1[4];
114	double p2[4];
115	double r1[4];
116	double r2[4];
117
118	static const Double_t cmtofm = 1.e13;
119	static const Double_t cmtoOneOverGeV = cmtofm*fgkWcons;
120
121	AliHBTParticle *part1 = pair->Particle1();
122	AliHBTParticle *part2 = pair->Particle2();
123
124	p1[0] = part1->Energy()*1000.0;
125	p1[1] = part1->Px()*1000.0;
126	p1[2] = part1->Py()*1000.0;
127	p1[3] = part1->Pz()*1000.0;
128
129	p2[0] = part2->Energy()*1000.0;
130	p2[1] = part2->Px()*1000.0;
131	p2[2] = part2->Py()*1000.0;
132	p2[3] = part2->Pz()*1000.0;
133
134	r1[0] = part1->T();
135	r1[1] = part1->Vx()*cmtofm;
136	r1[2] = part1->Vy()*cmtofm;
137	r1[3] = part1->Vz()*cmtofm;
138
139	r2[0] = part2->T();
140	r2[1] = part2->Vx()*cmtofm;
141	r2[2] = part2->Vy()*cmtofm;
142	r2[3] = part2->Vz()*cmtofm;
143
144	// END OF ALICE STUFF
145
146	// This code is written by Scott Pratt
147	// taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html
148	int alpha;
149	double kdotr;
150	if (fReducedMom)
151	{
152	const double momtest=4.0fMaxMomentumfMaxMomentum;
153	}
154	else
155	{
156	const double momtest=fMaxMomentum*fMaxMomentum;
157	}
158
159	double ptot2,pdotr,pp,rr;
160
161	if ( part1->GetPdgCode() == part2->GetPdgCode() )
162	{
163	*test=1;
164	mom=-(p2[0]-p1[0])(p2[0]-p1[0]);
165	for(alpha=1;alpha<4;alpha++){
166	mom=mom+(p2[alpha]-p1[alpha])*(p2[alpha]-p1[alpha]);
167	}
168	//#if ! defined MIXED_PAIRS_FOR_DENOM
169	// if(*mom>momtest){
170	// *test=0;
171	// return;
172	// }
173	//#endif
174	pp=(p1[0]+p2[0]);
175	rr=(r2[0]-r1[0]);
176	pdotr=pp*rr;
177	kdotr=(p2[0]-p1[0])*rr;
178	ptot2=pp*pp;
179	r=-rrrr;
180	for(alpha=1;alpha<4;alpha++){
181	pp=(p1[alpha]+p2[alpha]);
182	rr=(r2[alpha]-r1[alpha]);
183	pdotr=pdotr-pp*rr;
184	kdotr=kdotr-(p2[alpha]-p1[alpha])*rr;
185	ptot2=ptot2-pp*pp;
186	r=r+rr*rr;
187	}
188	mom=sqrt(mom);
189	qred=0.5*mom;
190
191	if (fReducedMom)
192	{
193	mom=qred;
194	}
195
196	qdotr=0.5kdotr;
197	r=sqrt(r+pdotr*pdotr/ptot2);
198	}
199	else //identical
200	{
201	// const double kdotp=MASS2MASS2-MASS1MASS1;
202	const double kdotp = part2->GetMass()part2->GetMass()- part1->GetMass()part1->GetMass();
203	*test=1;
204	mom=-(p2[0]-p1[0])(p2[0]-p1[0]);
205	ptot2=(p1[0]+p2[0])*(p1[0]+p2[0]);
206	for(alpha=1;alpha<4;alpha++){
207	mom=mom+(p2[alpha]-p1[alpha])*(p2[alpha]-p1[alpha]);
208	ptot2=ptot2-(p1[alpha]+p2[alpha])*(p1[alpha]+p2[alpha]);
209	}
210	mom=mom+kdotp*kdotp/ptot2;
211	//#if ! defined MIXED_PAIRS_FOR_DENOM
212	// if(*mom>momtest){
213	// *test=0;
214	// return;
215	// }
216	//#endif
217	pp=(p1[0]+p2[0]);
218	rr=(r2[0]-r1[0]);
219	pdotr=pp*rr;
220	kdotr=(p2[0]-p1[0])*rr;
221	r=-rrrr;
222	for(alpha=1;alpha<4;alpha++){
223	pp=(p1[alpha]+p2[alpha]);
224	rr=(r2[alpha]-r1[alpha]);
225	pdotr=pdotr-pp*rr;
226	kdotr=kdotr-(p2[alpha]-p1[alpha])*rr;
227	r=r+rr*rr;
228	}
229	kdotr=(-kdotr+kdotp*pdotr/ptot2);
230	mom=sqrt(mom);
231	qred=0.5*mom;
232
233	if (fReducedMom)
234	{
235	mom=qred;
236	}
237
238	qdotr=0.5kdotr;
239	r=sqrt(r+pdotr*pdotr/ptot2);
240	}//not identical
241
242	return;
243	}
244
245
246	//===================================================================
247
248	double AliHBTCrab::corrcalc(double trueqred,double trueqdotr,double truer)
249	{
250	//#define REDUCED_MOM
251	// This code is written by Scott Pratt
252	// taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html
253	double eta,arg,krmax,corr0;
254	double xx,xxprime,xxjj,p1,zk;
255	int jj,kk,ipart,ipartcount,ispin;
256	double pfactor,wsym_leftover,wanti_leftover,wnosym_leftover;
257	double qred,qdotr,r;
258	const double rmass=MASS1*MASS2/(MASS1+MASS2);
259	double_complex cphi1,cphi2,cphis,cphia;
260
261	arg=trueqdotr/197.323-2.0TMath::Pi()floor(trueqdotr/(197.3232.0TMath::Pi()));
262	cphi1=exp(ci*arg);
263	cphis=fgkROOT2*real(cphi1);
264	cphia=cifgkROOT2imag(cphi1);
265	corr0=real(INTERACTION_WSYMcphisconj(cphis)
266	+INTERACTION_WANTIcphiaconj(cphia)
267	+INTERACTION_WNOSYMcphi1conj(cphi1));
268	goto OUTSIDE_INTERACTION_RANGE;
269
270	#ifdef REDUCED_MOM
271	kk=(int)floor(trueqred/INTERACTION_DELK);
272	qred=(0.5+kk)*INTERACTION_DELK;
273	#else
274	kk=(int)floor(2.0*trueqred/INTERACTION_DELK);
275	qred=(0.5+kk)*INTERACTION_DELK/2.0;
276	#endif
277	qdotr=trueqdotr*qred/trueqred;
278	if(kk>=INTERACTION_NKMAX){
279	corr0=1.0;
280	goto OUTSIDE_INTERACTION_RANGE;
281	}
282	r=truer;
283
284	eta=0.0;
285	arg=qdotr/197.323-2.0TMath::Pi()floor(qdotr/(197.3232.0TMath::Pi()));
286	cphi1=exp(ci*arg);
287	cphi2=conj(cphi1);
288
289	cphis=(cphi1+cphi2)/fgkROOT2;
290	cphia=(cphi1-cphi2)/fgkROOT2;
291	corr0=0.0;
292	/* If there are corrections for strong interactions, add the
293	change for each partial wave. If npartial = 0 then there
294	are no strong int. corrections. */
295	wsym_leftover=INTERACTION_WSYM;
296	wanti_leftover=INTERACTION_WANTI;
297	wnosym_leftover=INTERACTION_WNOSYM;
298
299	corr0=corr0+real(wsym_leftovercphisconj(cphis)
300	+wanti_leftovercphiaconj(cphia)
301	+wnosym_leftovercphi1conj(cphi1));
302	OUTSIDE_INTERACTION_RANGE:
303
304	#ifdef BREIT_WIGNER
305	corr0=corr0+bwcalc(trueqred,truer);
306	#endif
307
308	return corr0;
309	}
310
311	double_complex AliHBTCrab::cgamma(double_complex c){
312	/* This calc.s gamma functions which are in the form gamma(n+i*y)
313	where n is an int and y is real. */
314	// This code is written by Scott Pratt
315	// taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html
316	double_complex cg,cphase;
317	int mm,j;
318	double x,y,phase,delp,cgmag;
319	x=real(c);
320	y=imag(c);
321	phase=-TMath::E()*y;
322	for(j=1;j<=100000;j++){
323	delp=(y/(double)j)-atan(y/(double)j);
324	phase=phase+delp;
325	if(fabs(delp)<1E-10) goto CGAMMA_ESCAPE;
326	}
327	printf("oops not accurate enough, increase jmax\n");
328	CGAMMA_ESCAPE:
329	phase=phase-2.0TMath::Pi()floor(phase/(2.0*TMath::Pi()));
330	cphase=exp(ci*phase);
331	cgmag=sqrt(TMath::Pi()y/sinh(TMath::Pi()y));
332	mm=(int)floor(x+0.5);
333	cg=cgmag*cphase;
334	if(mm<1){
335	for(j=1;j<=-mm+1;j++){
336	cg=cg/(1.0+(double)(-j)+ci*y);
337	}
338	}
339	if(mm>1) {
340	for(j=1;j<=mm-1;j++){
341	cg=cg((double)(j)+ciy);
342	}
343	}
344	return cg;
345	}
346	//===================================================================
347