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