[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 "AliVAODParticle.h"
#include "TDatabasePDG.h"
#include <TMath.h>

#include "volya_complex.h"

AliHBTCrab* AliHBTCrab::fgCrab = 0x0;

const Double_t AliHBTCrab::fgkWcons = 1./0.1973;
const Double_t AliHBTCrab::fgkROOT2=1.41421356237309504880;
#ifdef __DECCXX
const complex AliHBTCrab::fgkCI(0.0,1.0);
#else
const doublecomplex AliHBTCrab::fgkCI(0.0,1.0);
#endif

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

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),
  fPid1(0),
  fPid2(0),
  fMass1(0),
  fMass2(0),
  fInteractionWsym(0),
  fInteractionWanti(0),
  fInteractionWnosym(0),
  fInteractionDelk(0),
  fInteractionNkmax(0)
{
  //ctor
  if(fgCrab)
   {
     Fatal("AliHBTCrab","Do not use constructor directly. Use Instance() instead.");
   }
}
//===================================================================

AliHBTCrab::AliHBTCrab(const AliHBTCrab &/*source*/) :
  AliHBTWeights(),
  fBreitWigner(kFALSE),
  fReducedMom(kTRUE),
  fMaxMomentum(100.0),
  fPid1(0),
  fPid2(0),
  fMass1(0),
  fMass2(0),
  fInteractionWsym(0),
  fInteractionWanti(0),
  fInteractionWnosym(0),
  fInteractionDelk(0),
  fInteractionNkmax(0)
{
  //ctor
}
//===================================================================
AliHBTCrab & AliHBTCrab::operator=(const AliHBTCrab& /*source*/)
{
//cpy constructor
 return *AliHBTCrab::Instance();
}

void AliHBTCrab::Init(Int_t pid1,Int_t pid2)
{
//Initialization method
  fMass1 = TDatabasePDG::Instance()->GetParticle(pid1)->Mass();
  fMass2 = TDatabasePDG::Instance()->GetParticle(pid2)->Mass();
  fInteractionWsym = 1.0;
  fInteractionWanti = 0.0;
  fInteractionWnosym = 0.0;
  fInteractionDelk  = 1.0;
  fInteractionNkmax = 100;
     
  fPid1 = pid1;   
  fPid2 = pid2;
}
//===================================================================

Bool_t AliHBTCrab::SetConfig(const AliHBTPair* pair)
{
//returns the SetConfig
  
  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(AliHBTPair* partpair)
{
//returns the weight
  Double_t qred, r, qdotr, mom;
  Int_t test;
  
  SetConfig(partpair);
  
  GetComQuantities(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::GetComQuantities(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 kCmToFm = 1.e13;
// static const Double_t cmtoOneOverGeV = kCmToFm*fgkWcons;
 
 AliVAODParticle *part1 = pair->Particle1();
 AliVAODParticle *part2 = pair->Particle2();

 p1[0] = part1->E()*1000.0;
 p1[1] = part1->Px()*1000.0;
 p1[2] = part1->Py()*1000.0;
 p1[3] = part1->Pz()*1000.0;
 
 p2[0] = part2->E()*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()*kCmToFm;
 r1[2] = part1->Vy()*kCmToFm;
 r1[3] = part1->Vz()*kCmToFm;

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

//  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;
  double momtest;
  if (fReducedMom)
   {
     momtest=4.0*fMaxMomentum*fMaxMomentum;
   }  
  else
   {
     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=fMass2*fMass2-fMass1*fMass1;
   const double  kdotp = part2->Mass()*part2->Mass()- part1->Mass()*part1->Mass();
   *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,corr0;
//  double xx,xxprime,xxjj,p1,zk;
//  int jj,kk,ipart,ipartcount,ispin;
  int kk;
  double wsymleftover,wantileftover,wnosymleftover;
  double qred,qdotr,r;
//  const double rmass=fMass1*fMass2/(fMass1+fMass2);
#ifdef __DECCXX
  complex cphi1,cphi2,cphis,cphia;
#else
  doublecomplex cphi1,cphi2,cphis,cphia;
#endif

  arg=trueqdotr/197.323-2.0*TMath::Pi()*TMath::Floor(trueqdotr/(197.323*2.0*TMath::Pi()));
  cphi1=exp(fgkCI*arg);
  cphis=fgkROOT2*real(cphi1);
  cphia=fgkCI*fgkROOT2*imag(cphi1);
  corr0=real(fInteractionWsym*cphis*conj(cphis)
	     +fInteractionWanti*cphia*conj(cphia)
	     +fInteractionWnosym*cphi1*conj(cphi1));
  goto OUTSIDE_INTERACTION_RANGE;

#ifdef REDUCED_MOM
  kk=(int)TMath::Floor(trueqred/fInteractionDelk);
  qred=(0.5+kk)*fInteractionDelk;
#else
  kk=(int)TMath::Floor(2.0*trueqred/fInteractionDelk);
  qred=(0.5+kk)*fInteractionDelk/2.0;
#endif
  qdotr=trueqdotr*qred/trueqred;
  if(kk>=fInteractionNkmax){
    corr0=1.0;
    goto OUTSIDE_INTERACTION_RANGE;
  }
  r=truer;

  eta=0.0;
  arg=qdotr/197.323-2.0*TMath::Pi()*TMath::Floor(qdotr/(197.323*2.0*TMath::Pi()));
  cphi1=exp(fgkCI*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. */
  wsymleftover=fInteractionWsym;
  wantileftover=fInteractionWanti;
  wnosymleftover=fInteractionWnosym;

  corr0=corr0+real(wsymleftover*cphis*conj(cphis)
		   +wantileftover*cphia*conj(cphia)
		   +wnosymleftover*cphi1*conj(cphi1));
OUTSIDE_INTERACTION_RANGE:

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

  return corr0;
}

#ifdef __DECCXX
complex AliHBTCrab::CGamma(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
#else
doublecomplex AliHBTCrab::CGamma(doublecomplex 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
#endif
#ifdef __DECCXX
  complex cg,cphase;
#else
  doublecomplex cg,cphase;
#endif
  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(TMath::Abs(delp)<1E-10) goto CGamma_ESCAPE;
  }
  printf("oops not accurate enough, increase jmax\n");
CGamma_ESCAPE:
  phase=phase-2.0*TMath::Pi()*TMath::Floor(phase/(2.0*TMath::Pi()));
  cphase=exp(fgkCI*phase);
  cgmag=sqrt(TMath::Pi()*y/sinh(TMath::Pi()*y));
  mm=(int)TMath::Floor(x+0.5);
  cg=cgmag*cphase;
  if(mm<1){
    for(j=1;j<=-mm+1;j++){
      cg=cg/(1.0+(double)(-j)+fgkCI*y);
    }
  }
  if(mm>1) {
    for(j=1;j<=mm-1;j++){
      cg=cg*((double)(j)+fgkCI*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	// //
c114f4f8	9	// //
dd82cadc	10	/////////////////////////////////////////////////////////////////////////
f93c53d8	11
dd82cadc	12	#include "AliHBTPair.h"
78d7c6d3	13	#include "AliVAODParticle.h"
78d7c6d3	14	#include "TDatabasePDG.h"
dd82cadc	15	#include <TMath.h>
dd82cadc	16
	17	#include "volya_complex.h"
	18
dd82cadc	19	AliHBTCrab* AliHBTCrab::fgCrab = 0x0;
	20
	21	const Double_t AliHBTCrab::fgkWcons = 1./0.1973;
	22	const Double_t AliHBTCrab::fgkROOT2=1.41421356237309504880;
8882f543	23	#ifdef __DECCXX
c114f4f8	24	const complex AliHBTCrab::fgkCI(0.0,1.0);
8882f543	25	#else
c031a417	26	const doublecomplex AliHBTCrab::fgkCI(0.0,1.0);
8882f543	27	#endif
dd82cadc	28
	29	/************************************************************/
	30
	31	AliHBTCrab* AliHBTCrab::Instance()
	32	{
	33	// returns instance of class
	34	if (fgCrab == 0x0)
	35	{
	36	fgCrab = new AliHBTCrab();
	37	}
	38	return fgCrab;
	39	}
	40	//===================================================================
c114f4f8	41
dd82cadc	42	void AliHBTCrab::Set()
dd82cadc	43	{
c114f4f8	44	//sets this as weighitng class
dd82cadc	45	Info("Set","Setting CRAB as Weighing Class");
	46
	47	if ( fgWeights == 0x0 )
	48	{
	49	fgWeights = AliHBTCrab::Instance();
	50	return;
	51	}
	52	if ( fgWeights == AliHBTCrab::Instance() ) return;
	53	delete fgWeights;
	54	fgWeights = AliHBTCrab::Instance();
	55	}
	56	//===================================================================
	57
4b1c9620	58	AliHBTCrab::AliHBTCrab() :
	59	fBreitWigner(kFALSE),
	60	fReducedMom(kTRUE),
	61	fMaxMomentum(100.0),
	62	fPid1(0),
	63	fPid2(0),
	64	fMass1(0),
	65	fMass2(0),
	66	fInteractionWsym(0),
	67	fInteractionWanti(0),
	68	fInteractionWnosym(0),
	69	fInteractionDelk(0),
	70	fInteractionNkmax(0)
dd82cadc	71	{
dd82cadc	72	//ctor
8ef3740b	73	if(fgCrab)
	74	{
	75	Fatal("AliHBTCrab","Do not use constructor directly. Use Instance() instead.");
	76	}
dd82cadc	77	}
dd82cadc	78	//===================================================================
c031a417	79
4b1c9620	80	AliHBTCrab::AliHBTCrab(const AliHBTCrab &/source/) :
	81	AliHBTWeights(),
	82	fBreitWigner(kFALSE),
	83	fReducedMom(kTRUE),
	84	fMaxMomentum(100.0),
	85	fPid1(0),
	86	fPid2(0),
	87	fMass1(0),
	88	fMass2(0),
	89	fInteractionWsym(0),
	90	fInteractionWanti(0),
	91	fInteractionWnosym(0),
	92	fInteractionDelk(0),
	93	fInteractionNkmax(0)
c031a417	94	{
	95	//ctor
	96	}
	97	//===================================================================
34914285	98	AliHBTCrab & AliHBTCrab::operator=(const AliHBTCrab& /source/)
c031a417	99	{
	100	//cpy constructor
	101	return *AliHBTCrab::Instance();
	102	}
	103
dd82cadc	104	void AliHBTCrab::Init(Int_t pid1,Int_t pid2)
dd82cadc	105	{
c031a417	106	//Initialization method
c114f4f8	107	fMass1 = TDatabasePDG::Instance()->GetParticle(pid1)->Mass();
	108	fMass2 = TDatabasePDG::Instance()->GetParticle(pid2)->Mass();
	109	fInteractionWsym = 1.0;
	110	fInteractionWanti = 0.0;
	111	fInteractionWnosym = 0.0;
	112	fInteractionDelk = 1.0;
	113	fInteractionNkmax = 100;
dd82cadc	114
	115	fPid1 = pid1;
	116	fPid2 = pid2;
dd82cadc	117	}
c031a417	118	//===================================================================
dd82cadc	119
	120	Bool_t AliHBTCrab::SetConfig(const AliHBTPair* pair)
	121	{
c114f4f8	122	//returns the SetConfig
c114f4f8	123
dd82cadc	124	Int_t pdg1 = pair->Particle1()->GetPdgCode();
	125	Int_t pdg2 = pair->Particle2()->GetPdgCode();
	126
	127	if ( ( pdg1 == fPid1) && ( pdg2 == fPid2) ) return kFALSE;
	128	else Init (pdg1,pdg2);
	129
	130	return kTRUE;
	131	}
	132	//===================================================================
	133
1a1e58ac	134	Double_t AliHBTCrab::GetWeight(AliHBTPair* partpair)
dd82cadc	135	{
c114f4f8	136	//returns the weight
dd82cadc	137	Double_t qred, r, qdotr, mom;
	138	Int_t test;
	139
	140	SetConfig(partpair);
	141
c114f4f8	142	GetComQuantities(partpair, &qred, &r, &qdotr, &mom, &test);
dd82cadc	143	if(test==0)
	144	{
	145	Info("GetWeight","Test is 0");
	146	}
c114f4f8	147	Double_t corr = CorrCalc(qred,qdotr,r);
dd82cadc	148
	149	return corr;
	150	}
	151	//===================================================================
	152
c114f4f8	153	void AliHBTCrab::GetComQuantities(const AliHBTPair* pair,
dd82cadc	154	double qred,double r,double qdotr,double mom, int *test)
c031a417	155	{
dd82cadc	156	//************************************
	157	// ALICE //
	158
	159	double p1[4];
	160	double p2[4];
	161	double r1[4];
	162	double r2[4];
	163
c114f4f8	164	static const Double_t kCmToFm = 1.e13;
c114f4f8	165	// static const Double_t cmtoOneOverGeV = kCmToFm*fgkWcons;
dd82cadc	166
78d7c6d3	167	AliVAODParticle *part1 = pair->Particle1();
78d7c6d3	168	AliVAODParticle *part2 = pair->Particle2();
dd82cadc	169
78d7c6d3	170	p1[0] = part1->E()*1000.0;
dd82cadc	171	p1[1] = part1->Px()*1000.0;
	172	p1[2] = part1->Py()*1000.0;
	173	p1[3] = part1->Pz()*1000.0;
	174
78d7c6d3	175	p2[0] = part2->E()*1000.0;
dd82cadc	176	p2[1] = part2->Px()*1000.0;
	177	p2[2] = part2->Py()*1000.0;
	178	p2[3] = part2->Pz()*1000.0;
	179
	180	r1[0] = part1->T();
c114f4f8	181	r1[1] = part1->Vx()*kCmToFm;
	182	r1[2] = part1->Vy()*kCmToFm;
	183	r1[3] = part1->Vz()*kCmToFm;
dd82cadc	184
dd82cadc	185	r2[0] = part2->T();
c114f4f8	186	r2[1] = part2->Vx()*kCmToFm;
	187	r2[2] = part2->Vy()*kCmToFm;
	188	r2[3] = part2->Vz()*kCmToFm;
dd82cadc	189
	190	// END OF ALICE STUFF
	191
	192	// This code is written by Scott Pratt
	193	// taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html
	194	int alpha;
	195	double kdotr;
f93c53d8	196	double momtest;
dd82cadc	197	if (fReducedMom)
dd82cadc	198	{
f93c53d8	199	momtest=4.0fMaxMomentumfMaxMomentum;
dd82cadc	200	}
	201	else
	202	{
f93c53d8	203	momtest=fMaxMomentum*fMaxMomentum;
dd82cadc	204	}
	205
	206	double ptot2,pdotr,pp,rr;
	207
	208	if ( part1->GetPdgCode() == part2->GetPdgCode() )
	209	{
	210	*test=1;
	211	mom=-(p2[0]-p1[0])(p2[0]-p1[0]);
	212	for(alpha=1;alpha<4;alpha++){
	213	mom=mom+(p2[alpha]-p1[alpha])*(p2[alpha]-p1[alpha]);
	214	}
	215	//#if ! defined MIXED_PAIRS_FOR_DENOM
	216	// if(*mom>momtest){
	217	// *test=0;
	218	// return;
	219	// }
	220	//#endif
	221	pp=(p1[0]+p2[0]);
	222	rr=(r2[0]-r1[0]);
	223	pdotr=pp*rr;
	224	kdotr=(p2[0]-p1[0])*rr;
	225	ptot2=pp*pp;
	226	r=-rrrr;
	227	for(alpha=1;alpha<4;alpha++){
	228	pp=(p1[alpha]+p2[alpha]);
	229	rr=(r2[alpha]-r1[alpha]);
	230	pdotr=pdotr-pp*rr;
	231	kdotr=kdotr-(p2[alpha]-p1[alpha])*rr;
	232	ptot2=ptot2-pp*pp;
	233	r=r+rr*rr;
	234	}
	235	mom=sqrt(mom);
	236	qred=0.5*mom;
	237
	238	if (fReducedMom)
	239	{
	240	mom=qred;
	241	}
	242
	243	qdotr=0.5kdotr;
	244	r=sqrt(r+pdotr*pdotr/ptot2);
	245	}
	246	else //identical
	247	{
c114f4f8	248	// const double kdotp=fMass2fMass2-fMass1fMass1;
78d7c6d3	249	const double kdotp = part2->Mass()part2->Mass()- part1->Mass()part1->Mass();
dd82cadc	250	*test=1;
	251	mom=-(p2[0]-p1[0])(p2[0]-p1[0]);
	252	ptot2=(p1[0]+p2[0])*(p1[0]+p2[0]);
	253	for(alpha=1;alpha<4;alpha++){
	254	mom=mom+(p2[alpha]-p1[alpha])*(p2[alpha]-p1[alpha]);
	255	ptot2=ptot2-(p1[alpha]+p2[alpha])*(p1[alpha]+p2[alpha]);
	256	}
	257	mom=mom+kdotp*kdotp/ptot2;
	258	//#if ! defined MIXED_PAIRS_FOR_DENOM
	259	// if(*mom>momtest){
	260	// *test=0;
	261	// return;
	262	// }
	263	//#endif
	264	pp=(p1[0]+p2[0]);
	265	rr=(r2[0]-r1[0]);
	266	pdotr=pp*rr;
	267	kdotr=(p2[0]-p1[0])*rr;
	268	r=-rrrr;
	269	for(alpha=1;alpha<4;alpha++){
	270	pp=(p1[alpha]+p2[alpha]);
	271	rr=(r2[alpha]-r1[alpha]);
	272	pdotr=pdotr-pp*rr;
	273	kdotr=kdotr-(p2[alpha]-p1[alpha])*rr;
	274	r=r+rr*rr;
	275	}
	276	kdotr=(-kdotr+kdotp*pdotr/ptot2);
	277	mom=sqrt(mom);
	278	qred=0.5*mom;
	279
	280	if (fReducedMom)
	281	{
	282	mom=qred;
	283	}
	284
	285	qdotr=0.5kdotr;
	286	r=sqrt(r+pdotr*pdotr/ptot2);
	287	}//not identical
	288
	289	return;
	290	}
	291
	292
	293	//===================================================================
	294
c114f4f8	295	double AliHBTCrab::CorrCalc(double trueqred,double trueqdotr,double truer)
dd82cadc	296	{
	297	//#define REDUCED_MOM
	298	// This code is written by Scott Pratt
	299	// taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html
f93c53d8	300	double eta,arg,corr0;
	301	// double xx,xxprime,xxjj,p1,zk;
	302	// int jj,kk,ipart,ipartcount,ispin;
	303	int kk;
c031a417	304	double wsymleftover,wantileftover,wnosymleftover;
dd82cadc	305	double qred,qdotr,r;
c114f4f8	306	// const double rmass=fMass1*fMass2/(fMass1+fMass2);
8882f543	307	#ifdef __DECCXX
	308	complex cphi1,cphi2,cphis,cphia;
	309	#else
c031a417	310	doublecomplex cphi1,cphi2,cphis,cphia;
8882f543	311	#endif
dd82cadc	312
3294e9ff	313	arg=trueqdotr/197.323-2.0TMath::Pi()TMath::Floor(trueqdotr/(197.3232.0TMath::Pi()));
c114f4f8	314	cphi1=exp(fgkCI*arg);
dd82cadc	315	cphis=fgkROOT2*real(cphi1);
c114f4f8	316	cphia=fgkCIfgkROOT2imag(cphi1);
	317	corr0=real(fInteractionWsymcphisconj(cphis)
	318	+fInteractionWanticphiaconj(cphia)
	319	+fInteractionWnosymcphi1conj(cphi1));
dd82cadc	320	goto OUTSIDE_INTERACTION_RANGE;
	321
	322	#ifdef REDUCED_MOM
c114f4f8	323	kk=(int)TMath::Floor(trueqred/fInteractionDelk);
c114f4f8	324	qred=(0.5+kk)*fInteractionDelk;
dd82cadc	325	#else
c114f4f8	326	kk=(int)TMath::Floor(2.0*trueqred/fInteractionDelk);
c114f4f8	327	qred=(0.5+kk)*fInteractionDelk/2.0;
dd82cadc	328	#endif
dd82cadc	329	qdotr=trueqdotr*qred/trueqred;
c114f4f8	330	if(kk>=fInteractionNkmax){
dd82cadc	331	corr0=1.0;
	332	goto OUTSIDE_INTERACTION_RANGE;
	333	}
	334	r=truer;
	335
	336	eta=0.0;
3294e9ff	337	arg=qdotr/197.323-2.0TMath::Pi()TMath::Floor(qdotr/(197.3232.0TMath::Pi()));
c114f4f8	338	cphi1=exp(fgkCI*arg);
dd82cadc	339	cphi2=conj(cphi1);
	340
	341	cphis=(cphi1+cphi2)/fgkROOT2;
	342	cphia=(cphi1-cphi2)/fgkROOT2;
	343	corr0=0.0;
	344	/* If there are corrections for strong interactions, add the
	345	change for each partial wave. If npartial = 0 then there
	346	are no strong int. corrections. */
c031a417	347	wsymleftover=fInteractionWsym;
	348	wantileftover=fInteractionWanti;
	349	wnosymleftover=fInteractionWnosym;
dd82cadc	350
c031a417	351	corr0=corr0+real(wsymleftovercphisconj(cphis)
	352	+wantileftovercphiaconj(cphia)
	353	+wnosymleftovercphi1conj(cphi1));
dd82cadc	354	OUTSIDE_INTERACTION_RANGE:
	355
	356	#ifdef BREIT_WIGNER
	357	corr0=corr0+bwcalc(trueqred,truer);
	358	#endif
	359
	360	return corr0;
	361	}
	362
8882f543	363	#ifdef __DECCXX
c114f4f8	364	complex AliHBTCrab::CGamma(complex c){
c031a417	365	/* This calc.s gamma functions which are in the form gamma(n+i*y)
	366	where n is an int and y is real. */
	367	// This code is written by Scott Pratt
	368	// taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html
8882f543	369	#else
c031a417	370	doublecomplex AliHBTCrab::CGamma(doublecomplex c){
	371	/* This calc.s gamma functions which are in the form gamma(n+i*y)
	372	where n is an int and y is real. */
dd82cadc	373	// This code is written by Scott Pratt
dd82cadc	374	// taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html
c031a417	375	#endif
8882f543	376	#ifdef __DECCXX
	377	complex cg,cphase;
	378	#else
c031a417	379	doublecomplex cg,cphase;
8882f543	380	#endif
dd82cadc	381	int mm,j;
	382	double x,y,phase,delp,cgmag;
	383	x=real(c);
	384	y=imag(c);
	385	phase=-TMath::E()*y;
	386	for(j=1;j<=100000;j++){
	387	delp=(y/(double)j)-atan(y/(double)j);
	388	phase=phase+delp;
c114f4f8	389	if(TMath::Abs(delp)<1E-10) goto CGamma_ESCAPE;
dd82cadc	390	}
dd82cadc	391	printf("oops not accurate enough, increase jmax\n");
c114f4f8	392	CGamma_ESCAPE:
3294e9ff	393	phase=phase-2.0TMath::Pi()TMath::Floor(phase/(2.0*TMath::Pi()));
c114f4f8	394	cphase=exp(fgkCI*phase);
dd82cadc	395	cgmag=sqrt(TMath::Pi()y/sinh(TMath::Pi()y));
3294e9ff	396	mm=(int)TMath::Floor(x+0.5);
dd82cadc	397	cg=cgmag*cphase;
	398	if(mm<1){
	399	for(j=1;j<=-mm+1;j++){
c114f4f8	400	cg=cg/(1.0+(double)(-j)+fgkCI*y);
dd82cadc	401	}
	402	}
	403	if(mm>1) {
	404	for(j=1;j<=mm-1;j++){
c114f4f8	405	cg=cg((double)(j)+fgkCIy);
dd82cadc	406	}
	407	}
	408	return cg;
	409	}