1 #include "AliHBTCrab.h"
2 //______________________________________________________________________
3 /////////////////////////////////////////////////////////////////////////
5 // AliRoot wrapper to CRAB //
6 // taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html //
7 // written by Scott Pratt //
10 /////////////////////////////////////////////////////////////////////////
12 #include "AliHBTPair.h"
17 #include "volya_complex.h"
20 //using namespace std;
21 AliHBTCrab* AliHBTCrab::fgCrab = 0x0;
23 const Double_t AliHBTCrab::fgkWcons = 1./0.1973;
24 const Double_t AliHBTCrab::fgkROOT2=1.41421356237309504880;
26 const complex AliHBTCrab::fgkCI(0.0,1.0);
28 const double_complex AliHBTCrab::fgkCI(0.0,1.0);
31 /************************************************************/
33 AliHBTCrab* AliHBTCrab::Instance()
35 // returns instance of class
38 fgCrab = new AliHBTCrab();
42 //===================================================================
44 void AliHBTCrab::Set()
46 //sets this as weighitng class
47 Info("Set","Setting CRAB as Weighing Class");
49 if ( fgWeights == 0x0 )
51 fgWeights = AliHBTCrab::Instance();
54 if ( fgWeights == AliHBTCrab::Instance() ) return;
56 fgWeights = AliHBTCrab::Instance();
58 //===================================================================
60 AliHBTCrab::AliHBTCrab():
67 //===================================================================
68 void AliHBTCrab::Init(Int_t pid1,Int_t pid2)
70 fMass1 = TDatabasePDG::Instance()->GetParticle(pid1)->Mass();
71 fMass2 = TDatabasePDG::Instance()->GetParticle(pid2)->Mass();
72 fInteractionWsym = 1.0;
73 fInteractionWanti = 0.0;
74 fInteractionWnosym = 0.0;
75 fInteractionDelk = 1.0;
76 fInteractionNkmax = 100;
82 Bool_t AliHBTCrab::SetConfig(const AliHBTPair* pair)
84 //returns the SetConfig
86 Int_t pdg1 = pair->Particle1()->GetPdgCode();
87 Int_t pdg2 = pair->Particle2()->GetPdgCode();
89 if ( ( pdg1 == fPid1) && ( pdg2 == fPid2) ) return kFALSE;
90 else Init (pdg1,pdg2);
94 //===================================================================
96 Double_t AliHBTCrab::GetWeight(const AliHBTPair* partpair)
99 Double_t qred, r, qdotr, mom;
104 GetComQuantities(partpair, &qred, &r, &qdotr, &mom, &test);
107 Info("GetWeight","Test is 0");
109 Double_t corr = CorrCalc(qred,qdotr,r);
113 //===================================================================
115 void AliHBTCrab::GetComQuantities(const AliHBTPair* pair,
116 double *qred,double *r,double *qdotr,double *mom, int *test)
118 //************************************
126 static const Double_t kCmToFm = 1.e13;
127 // static const Double_t cmtoOneOverGeV = kCmToFm*fgkWcons;
129 AliHBTParticle *part1 = pair->Particle1();
130 AliHBTParticle *part2 = pair->Particle2();
132 p1[0] = part1->Energy()*1000.0;
133 p1[1] = part1->Px()*1000.0;
134 p1[2] = part1->Py()*1000.0;
135 p1[3] = part1->Pz()*1000.0;
137 p2[0] = part2->Energy()*1000.0;
138 p2[1] = part2->Px()*1000.0;
139 p2[2] = part2->Py()*1000.0;
140 p2[3] = part2->Pz()*1000.0;
143 r1[1] = part1->Vx()*kCmToFm;
144 r1[2] = part1->Vy()*kCmToFm;
145 r1[3] = part1->Vz()*kCmToFm;
148 r2[1] = part2->Vx()*kCmToFm;
149 r2[2] = part2->Vy()*kCmToFm;
150 r2[3] = part2->Vz()*kCmToFm;
152 // END OF ALICE STUFF
154 // This code is written by Scott Pratt
155 // taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html
161 momtest=4.0*fMaxMomentum*fMaxMomentum;
165 momtest=fMaxMomentum*fMaxMomentum;
168 double ptot2,pdotr,pp,rr;
170 if ( part1->GetPdgCode() == part2->GetPdgCode() )
173 *mom=-(p2[0]-p1[0])*(p2[0]-p1[0]);
174 for(alpha=1;alpha<4;alpha++){
175 *mom=*mom+(p2[alpha]-p1[alpha])*(p2[alpha]-p1[alpha]);
177 //#if ! defined MIXED_PAIRS_FOR_DENOM
186 kdotr=(p2[0]-p1[0])*rr;
189 for(alpha=1;alpha<4;alpha++){
190 pp=(p1[alpha]+p2[alpha]);
191 rr=(r2[alpha]-r1[alpha]);
193 kdotr=kdotr-(p2[alpha]-p1[alpha])*rr;
206 *r=sqrt(*r+pdotr*pdotr/ptot2);
210 // const double kdotp=fMass2*fMass2-fMass1*fMass1;
211 const double kdotp = part2->GetMass()*part2->GetMass()- part1->GetMass()*part1->GetMass();
213 *mom=-(p2[0]-p1[0])*(p2[0]-p1[0]);
214 ptot2=(p1[0]+p2[0])*(p1[0]+p2[0]);
215 for(alpha=1;alpha<4;alpha++){
216 *mom=*mom+(p2[alpha]-p1[alpha])*(p2[alpha]-p1[alpha]);
217 ptot2=ptot2-(p1[alpha]+p2[alpha])*(p1[alpha]+p2[alpha]);
219 *mom=*mom+kdotp*kdotp/ptot2;
220 //#if ! defined MIXED_PAIRS_FOR_DENOM
229 kdotr=(p2[0]-p1[0])*rr;
231 for(alpha=1;alpha<4;alpha++){
232 pp=(p1[alpha]+p2[alpha]);
233 rr=(r2[alpha]-r1[alpha]);
235 kdotr=kdotr-(p2[alpha]-p1[alpha])*rr;
238 kdotr=(-kdotr+kdotp*pdotr/ptot2);
248 *r=sqrt(*r+pdotr*pdotr/ptot2);
255 //===================================================================
257 double AliHBTCrab::CorrCalc(double trueqred,double trueqdotr,double truer)
259 //#define REDUCED_MOM
260 // This code is written by Scott Pratt
261 // taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html
262 double eta,arg,corr0;
263 // double xx,xxprime,xxjj,p1,zk;
264 // int jj,kk,ipart,ipartcount,ispin;
266 double wsym_leftover,wanti_leftover,wnosym_leftover;
268 // const double rmass=fMass1*fMass2/(fMass1+fMass2);
270 complex cphi1,cphi2,cphis,cphia;
272 double_complex cphi1,cphi2,cphis,cphia;
275 arg=trueqdotr/197.323-2.0*TMath::Pi()*TMath::Floor(trueqdotr/(197.323*2.0*TMath::Pi()));
276 cphi1=exp(fgkCI*arg);
277 cphis=fgkROOT2*real(cphi1);
278 cphia=fgkCI*fgkROOT2*imag(cphi1);
279 corr0=real(fInteractionWsym*cphis*conj(cphis)
280 +fInteractionWanti*cphia*conj(cphia)
281 +fInteractionWnosym*cphi1*conj(cphi1));
282 goto OUTSIDE_INTERACTION_RANGE;
285 kk=(int)TMath::Floor(trueqred/fInteractionDelk);
286 qred=(0.5+kk)*fInteractionDelk;
288 kk=(int)TMath::Floor(2.0*trueqred/fInteractionDelk);
289 qred=(0.5+kk)*fInteractionDelk/2.0;
291 qdotr=trueqdotr*qred/trueqred;
292 if(kk>=fInteractionNkmax){
294 goto OUTSIDE_INTERACTION_RANGE;
299 arg=qdotr/197.323-2.0*TMath::Pi()*TMath::Floor(qdotr/(197.323*2.0*TMath::Pi()));
300 cphi1=exp(fgkCI*arg);
303 cphis=(cphi1+cphi2)/fgkROOT2;
304 cphia=(cphi1-cphi2)/fgkROOT2;
306 /* If there are corrections for strong interactions, add the
307 change for each partial wave. If npartial = 0 then there
308 are no strong int. corrections. */
309 wsym_leftover=fInteractionWsym;
310 wanti_leftover=fInteractionWanti;
311 wnosym_leftover=fInteractionWnosym;
313 corr0=corr0+real(wsym_leftover*cphis*conj(cphis)
314 +wanti_leftover*cphia*conj(cphia)
315 +wnosym_leftover*cphi1*conj(cphi1));
316 OUTSIDE_INTERACTION_RANGE:
319 corr0=corr0+bwcalc(trueqred,truer);
326 complex AliHBTCrab::CGamma(complex c){
328 double_complex AliHBTCrab::CGamma(double_complex c){
330 /* This calc.s gamma functions which are in the form gamma(n+i*y)
331 where n is an int and y is real. */
332 // This code is written by Scott Pratt
333 // taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html
337 double_complex cg,cphase;
340 double x,y,phase,delp,cgmag;
344 for(j=1;j<=100000;j++){
345 delp=(y/(double)j)-atan(y/(double)j);
347 if(TMath::Abs(delp)<1E-10) goto CGamma_ESCAPE;
349 printf("oops not accurate enough, increase jmax\n");
351 phase=phase-2.0*TMath::Pi()*TMath::Floor(phase/(2.0*TMath::Pi()));
352 cphase=exp(fgkCI*phase);
353 cgmag=sqrt(TMath::Pi()*y/sinh(TMath::Pi()*y));
354 mm=(int)TMath::Floor(x+0.5);
357 for(j=1;j<=-mm+1;j++){
358 cg=cg/(1.0+(double)(-j)+fgkCI*y);
362 for(j=1;j<=mm-1;j++){
363 cg=cg*((double)(j)+fgkCI*y);