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"
13 #include "AliVAODParticle.h"
14 #include "TDatabasePDG.h"
17 #include "volya_complex.h"
19 AliHBTCrab* AliHBTCrab::fgCrab = 0x0;
21 const Double_t AliHBTCrab::fgkWcons = 1./0.1973;
22 const Double_t AliHBTCrab::fgkROOT2=1.41421356237309504880;
24 const complex AliHBTCrab::fgkCI(0.0,1.0);
26 const doublecomplex AliHBTCrab::fgkCI(0.0,1.0);
29 /************************************************************/
31 AliHBTCrab* AliHBTCrab::Instance()
33 // returns instance of class
36 fgCrab = new AliHBTCrab();
40 //===================================================================
42 void AliHBTCrab::Set()
44 //sets this as weighitng class
45 Info("Set","Setting CRAB as Weighing Class");
47 if ( fgWeights == 0x0 )
49 fgWeights = AliHBTCrab::Instance();
52 if ( fgWeights == AliHBTCrab::Instance() ) return;
54 fgWeights = AliHBTCrab::Instance();
56 //===================================================================
58 AliHBTCrab::AliHBTCrab() :
68 fInteractionWnosym(0),
75 Fatal("AliHBTCrab","Do not use constructor directly. Use Instance() instead.");
78 //===================================================================
80 AliHBTCrab::AliHBTCrab(const AliHBTCrab &/*source*/) :
91 fInteractionWnosym(0),
97 //===================================================================
98 AliHBTCrab & AliHBTCrab::operator=(const AliHBTCrab& /*source*/)
101 return *AliHBTCrab::Instance();
104 void AliHBTCrab::Init(Int_t pid1,Int_t pid2)
106 //Initialization method
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;
118 //===================================================================
120 Bool_t AliHBTCrab::SetConfig(const AliHBTPair* pair)
122 //returns the SetConfig
124 Int_t pdg1 = pair->Particle1()->GetPdgCode();
125 Int_t pdg2 = pair->Particle2()->GetPdgCode();
127 if ( ( pdg1 == fPid1) && ( pdg2 == fPid2) ) return kFALSE;
128 else Init (pdg1,pdg2);
132 //===================================================================
134 Double_t AliHBTCrab::GetWeight(AliHBTPair* partpair)
137 Double_t qred, r, qdotr, mom;
142 GetComQuantities(partpair, &qred, &r, &qdotr, &mom, &test);
145 Info("GetWeight","Test is 0");
147 Double_t corr = CorrCalc(qred,qdotr,r);
151 //===================================================================
153 void AliHBTCrab::GetComQuantities(const AliHBTPair* pair,
154 double *qred,double *r,double *qdotr,double *mom, int *test)
156 //************************************
164 static const Double_t kCmToFm = 1.e13;
165 // static const Double_t cmtoOneOverGeV = kCmToFm*fgkWcons;
167 AliVAODParticle *part1 = pair->Particle1();
168 AliVAODParticle *part2 = pair->Particle2();
170 p1[0] = part1->E()*1000.0;
171 p1[1] = part1->Px()*1000.0;
172 p1[2] = part1->Py()*1000.0;
173 p1[3] = part1->Pz()*1000.0;
175 p2[0] = part2->E()*1000.0;
176 p2[1] = part2->Px()*1000.0;
177 p2[2] = part2->Py()*1000.0;
178 p2[3] = part2->Pz()*1000.0;
181 r1[1] = part1->Vx()*kCmToFm;
182 r1[2] = part1->Vy()*kCmToFm;
183 r1[3] = part1->Vz()*kCmToFm;
186 r2[1] = part2->Vx()*kCmToFm;
187 r2[2] = part2->Vy()*kCmToFm;
188 r2[3] = part2->Vz()*kCmToFm;
190 // END OF ALICE STUFF
192 // This code is written by Scott Pratt
193 // taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html
199 momtest=4.0*fMaxMomentum*fMaxMomentum;
203 momtest=fMaxMomentum*fMaxMomentum;
206 double ptot2,pdotr,pp,rr;
208 if ( part1->GetPdgCode() == part2->GetPdgCode() )
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]);
215 //#if ! defined MIXED_PAIRS_FOR_DENOM
224 kdotr=(p2[0]-p1[0])*rr;
227 for(alpha=1;alpha<4;alpha++){
228 pp=(p1[alpha]+p2[alpha]);
229 rr=(r2[alpha]-r1[alpha]);
231 kdotr=kdotr-(p2[alpha]-p1[alpha])*rr;
244 *r=sqrt(*r+pdotr*pdotr/ptot2);
248 // const double kdotp=fMass2*fMass2-fMass1*fMass1;
249 const double kdotp = part2->Mass()*part2->Mass()- part1->Mass()*part1->Mass();
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]);
257 *mom=*mom+kdotp*kdotp/ptot2;
258 //#if ! defined MIXED_PAIRS_FOR_DENOM
267 kdotr=(p2[0]-p1[0])*rr;
269 for(alpha=1;alpha<4;alpha++){
270 pp=(p1[alpha]+p2[alpha]);
271 rr=(r2[alpha]-r1[alpha]);
273 kdotr=kdotr-(p2[alpha]-p1[alpha])*rr;
276 kdotr=(-kdotr+kdotp*pdotr/ptot2);
286 *r=sqrt(*r+pdotr*pdotr/ptot2);
293 //===================================================================
295 double AliHBTCrab::CorrCalc(double trueqred,double trueqdotr,double truer)
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
300 double eta,arg,corr0;
301 // double xx,xxprime,xxjj,p1,zk;
302 // int jj,kk,ipart,ipartcount,ispin;
304 double wsymleftover,wantileftover,wnosymleftover;
306 // const double rmass=fMass1*fMass2/(fMass1+fMass2);
308 complex cphi1,cphi2,cphis,cphia;
310 doublecomplex cphi1,cphi2,cphis,cphia;
313 arg=trueqdotr/197.323-2.0*TMath::Pi()*TMath::Floor(trueqdotr/(197.323*2.0*TMath::Pi()));
314 cphi1=exp(fgkCI*arg);
315 cphis=fgkROOT2*real(cphi1);
316 cphia=fgkCI*fgkROOT2*imag(cphi1);
317 corr0=real(fInteractionWsym*cphis*conj(cphis)
318 +fInteractionWanti*cphia*conj(cphia)
319 +fInteractionWnosym*cphi1*conj(cphi1));
320 goto OUTSIDE_INTERACTION_RANGE;
323 kk=(int)TMath::Floor(trueqred/fInteractionDelk);
324 qred=(0.5+kk)*fInteractionDelk;
326 kk=(int)TMath::Floor(2.0*trueqred/fInteractionDelk);
327 qred=(0.5+kk)*fInteractionDelk/2.0;
329 qdotr=trueqdotr*qred/trueqred;
330 if(kk>=fInteractionNkmax){
332 goto OUTSIDE_INTERACTION_RANGE;
337 arg=qdotr/197.323-2.0*TMath::Pi()*TMath::Floor(qdotr/(197.323*2.0*TMath::Pi()));
338 cphi1=exp(fgkCI*arg);
341 cphis=(cphi1+cphi2)/fgkROOT2;
342 cphia=(cphi1-cphi2)/fgkROOT2;
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. */
347 wsymleftover=fInteractionWsym;
348 wantileftover=fInteractionWanti;
349 wnosymleftover=fInteractionWnosym;
351 corr0=corr0+real(wsymleftover*cphis*conj(cphis)
352 +wantileftover*cphia*conj(cphia)
353 +wnosymleftover*cphi1*conj(cphi1));
354 OUTSIDE_INTERACTION_RANGE:
357 corr0=corr0+bwcalc(trueqred,truer);
364 complex AliHBTCrab::CGamma(complex c){
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
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. */
373 // This code is written by Scott Pratt
374 // taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html
379 doublecomplex cg,cphase;
382 double x,y,phase,delp,cgmag;
386 for(j=1;j<=100000;j++){
387 delp=(y/(double)j)-atan(y/(double)j);
389 if(TMath::Abs(delp)<1E-10) goto CGamma_ESCAPE;
391 printf("oops not accurate enough, increase jmax\n");
393 phase=phase-2.0*TMath::Pi()*TMath::Floor(phase/(2.0*TMath::Pi()));
394 cphase=exp(fgkCI*phase);
395 cgmag=sqrt(TMath::Pi()*y/sinh(TMath::Pi()*y));
396 mm=(int)TMath::Floor(x+0.5);
399 for(j=1;j<=-mm+1;j++){
400 cg=cg/(1.0+(double)(-j)+fgkCI*y);
404 for(j=1;j<=mm-1;j++){
405 cg=cg*((double)(j)+fgkCI*y);