]>
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" |
13 | ||
14 | #include <TMath.h> | |
15 | #include <TPDGCode.h> | |
16 | ||
17 | #include "volya_complex.h" | |
18 | ||
19 | //#include <complex> | |
20 | //using namespace std; | |
21 | AliHBTCrab* AliHBTCrab::fgCrab = 0x0; | |
22 | ||
23 | const Double_t AliHBTCrab::fgkWcons = 1./0.1973; | |
24 | const Double_t AliHBTCrab::fgkROOT2=1.41421356237309504880; | |
8882f543 | 25 | #ifdef __DECCXX |
c114f4f8 | 26 | const complex AliHBTCrab::fgkCI(0.0,1.0); |
8882f543 | 27 | #else |
c114f4f8 | 28 | const double_complex AliHBTCrab::fgkCI(0.0,1.0); |
8882f543 | 29 | #endif |
dd82cadc | 30 | |
31 | /************************************************************/ | |
32 | ||
33 | AliHBTCrab* AliHBTCrab::Instance() | |
34 | { | |
35 | // returns instance of class | |
36 | if (fgCrab == 0x0) | |
37 | { | |
38 | fgCrab = new AliHBTCrab(); | |
39 | } | |
40 | return fgCrab; | |
41 | } | |
42 | //=================================================================== | |
c114f4f8 | 43 | |
dd82cadc | 44 | void AliHBTCrab::Set() |
45 | { | |
c114f4f8 | 46 | //sets this as weighitng class |
dd82cadc | 47 | Info("Set","Setting CRAB as Weighing Class"); |
48 | ||
49 | if ( fgWeights == 0x0 ) | |
50 | { | |
51 | fgWeights = AliHBTCrab::Instance(); | |
52 | return; | |
53 | } | |
54 | if ( fgWeights == AliHBTCrab::Instance() ) return; | |
55 | delete fgWeights; | |
56 | fgWeights = AliHBTCrab::Instance(); | |
57 | } | |
58 | //=================================================================== | |
59 | ||
60 | AliHBTCrab::AliHBTCrab(): | |
61 | fBreitWigner(kFALSE), | |
62 | fReducedMom(kTRUE), | |
63 | fMaxMomentum(100.0) | |
64 | { | |
65 | //ctor | |
66 | } | |
67 | //=================================================================== | |
68 | void AliHBTCrab::Init(Int_t pid1,Int_t pid2) | |
69 | { | |
c114f4f8 | 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; | |
dd82cadc | 77 | |
78 | fPid1 = pid1; | |
79 | fPid2 = pid2; | |
dd82cadc | 80 | } |
81 | ||
82 | Bool_t AliHBTCrab::SetConfig(const AliHBTPair* pair) | |
83 | { | |
c114f4f8 | 84 | //returns the SetConfig |
85 | ||
dd82cadc | 86 | Int_t pdg1 = pair->Particle1()->GetPdgCode(); |
87 | Int_t pdg2 = pair->Particle2()->GetPdgCode(); | |
88 | ||
89 | if ( ( pdg1 == fPid1) && ( pdg2 == fPid2) ) return kFALSE; | |
90 | else Init (pdg1,pdg2); | |
91 | ||
92 | return kTRUE; | |
93 | } | |
94 | //=================================================================== | |
95 | ||
96 | Double_t AliHBTCrab::GetWeight(const AliHBTPair* partpair) | |
97 | { | |
c114f4f8 | 98 | //returns the weight |
dd82cadc | 99 | Double_t qred, r, qdotr, mom; |
100 | Int_t test; | |
101 | ||
102 | SetConfig(partpair); | |
103 | ||
c114f4f8 | 104 | GetComQuantities(partpair, &qred, &r, &qdotr, &mom, &test); |
dd82cadc | 105 | if(test==0) |
106 | { | |
107 | Info("GetWeight","Test is 0"); | |
108 | } | |
c114f4f8 | 109 | Double_t corr = CorrCalc(qred,qdotr,r); |
dd82cadc | 110 | |
111 | return corr; | |
112 | } | |
113 | //=================================================================== | |
114 | ||
c114f4f8 | 115 | void AliHBTCrab::GetComQuantities(const AliHBTPair* pair, |
dd82cadc | 116 | double *qred,double *r,double *qdotr,double *mom, int *test) |
117 | { | |
118 | //************************************ | |
119 | // ALICE // | |
120 | ||
121 | double p1[4]; | |
122 | double p2[4]; | |
123 | double r1[4]; | |
124 | double r2[4]; | |
125 | ||
c114f4f8 | 126 | static const Double_t kCmToFm = 1.e13; |
127 | // static const Double_t cmtoOneOverGeV = kCmToFm*fgkWcons; | |
dd82cadc | 128 | |
129 | AliHBTParticle *part1 = pair->Particle1(); | |
130 | AliHBTParticle *part2 = pair->Particle2(); | |
131 | ||
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; | |
136 | ||
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; | |
141 | ||
142 | r1[0] = part1->T(); | |
c114f4f8 | 143 | r1[1] = part1->Vx()*kCmToFm; |
144 | r1[2] = part1->Vy()*kCmToFm; | |
145 | r1[3] = part1->Vz()*kCmToFm; | |
dd82cadc | 146 | |
147 | r2[0] = part2->T(); | |
c114f4f8 | 148 | r2[1] = part2->Vx()*kCmToFm; |
149 | r2[2] = part2->Vy()*kCmToFm; | |
150 | r2[3] = part2->Vz()*kCmToFm; | |
dd82cadc | 151 | |
152 | // END OF ALICE STUFF | |
153 | ||
154 | // This code is written by Scott Pratt | |
155 | // taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html | |
156 | int alpha; | |
157 | double kdotr; | |
f93c53d8 | 158 | double momtest; |
dd82cadc | 159 | if (fReducedMom) |
160 | { | |
f93c53d8 | 161 | momtest=4.0*fMaxMomentum*fMaxMomentum; |
dd82cadc | 162 | } |
163 | else | |
164 | { | |
f93c53d8 | 165 | momtest=fMaxMomentum*fMaxMomentum; |
dd82cadc | 166 | } |
167 | ||
168 | double ptot2,pdotr,pp,rr; | |
169 | ||
170 | if ( part1->GetPdgCode() == part2->GetPdgCode() ) | |
171 | { | |
172 | *test=1; | |
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]); | |
176 | } | |
177 | //#if ! defined MIXED_PAIRS_FOR_DENOM | |
178 | // if(*mom>momtest){ | |
179 | // *test=0; | |
180 | // return; | |
181 | // } | |
182 | //#endif | |
183 | pp=(p1[0]+p2[0]); | |
184 | rr=(r2[0]-r1[0]); | |
185 | pdotr=pp*rr; | |
186 | kdotr=(p2[0]-p1[0])*rr; | |
187 | ptot2=pp*pp; | |
188 | *r=-rr*rr; | |
189 | for(alpha=1;alpha<4;alpha++){ | |
190 | pp=(p1[alpha]+p2[alpha]); | |
191 | rr=(r2[alpha]-r1[alpha]); | |
192 | pdotr=pdotr-pp*rr; | |
193 | kdotr=kdotr-(p2[alpha]-p1[alpha])*rr; | |
194 | ptot2=ptot2-pp*pp; | |
195 | *r=*r+rr*rr; | |
196 | } | |
197 | *mom=sqrt(*mom); | |
198 | *qred=0.5**mom; | |
199 | ||
200 | if (fReducedMom) | |
201 | { | |
202 | *mom=*qred; | |
203 | } | |
204 | ||
205 | *qdotr=0.5*kdotr; | |
206 | *r=sqrt(*r+pdotr*pdotr/ptot2); | |
207 | } | |
208 | else //identical | |
209 | { | |
c114f4f8 | 210 | // const double kdotp=fMass2*fMass2-fMass1*fMass1; |
dd82cadc | 211 | const double kdotp = part2->GetMass()*part2->GetMass()- part1->GetMass()*part1->GetMass(); |
212 | *test=1; | |
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]); | |
218 | } | |
219 | *mom=*mom+kdotp*kdotp/ptot2; | |
220 | //#if ! defined MIXED_PAIRS_FOR_DENOM | |
221 | // if(*mom>momtest){ | |
222 | // *test=0; | |
223 | // return; | |
224 | // } | |
225 | //#endif | |
226 | pp=(p1[0]+p2[0]); | |
227 | rr=(r2[0]-r1[0]); | |
228 | pdotr=pp*rr; | |
229 | kdotr=(p2[0]-p1[0])*rr; | |
230 | *r=-rr*rr; | |
231 | for(alpha=1;alpha<4;alpha++){ | |
232 | pp=(p1[alpha]+p2[alpha]); | |
233 | rr=(r2[alpha]-r1[alpha]); | |
234 | pdotr=pdotr-pp*rr; | |
235 | kdotr=kdotr-(p2[alpha]-p1[alpha])*rr; | |
236 | *r=*r+rr*rr; | |
237 | } | |
238 | kdotr=(-kdotr+kdotp*pdotr/ptot2); | |
239 | *mom=sqrt(*mom); | |
240 | *qred=0.5**mom; | |
241 | ||
242 | if (fReducedMom) | |
243 | { | |
244 | *mom=*qred; | |
245 | } | |
246 | ||
247 | *qdotr=0.5*kdotr; | |
248 | *r=sqrt(*r+pdotr*pdotr/ptot2); | |
249 | }//not identical | |
250 | ||
251 | return; | |
252 | } | |
253 | ||
254 | ||
255 | //=================================================================== | |
256 | ||
c114f4f8 | 257 | double AliHBTCrab::CorrCalc(double trueqred,double trueqdotr,double truer) |
dd82cadc | 258 | { |
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 | |
f93c53d8 | 262 | double eta,arg,corr0; |
263 | // double xx,xxprime,xxjj,p1,zk; | |
264 | // int jj,kk,ipart,ipartcount,ispin; | |
265 | int kk; | |
266 | double wsym_leftover,wanti_leftover,wnosym_leftover; | |
dd82cadc | 267 | double qred,qdotr,r; |
c114f4f8 | 268 | // const double rmass=fMass1*fMass2/(fMass1+fMass2); |
8882f543 | 269 | #ifdef __DECCXX |
270 | complex cphi1,cphi2,cphis,cphia; | |
271 | #else | |
dd82cadc | 272 | double_complex cphi1,cphi2,cphis,cphia; |
8882f543 | 273 | #endif |
dd82cadc | 274 | |
3294e9ff | 275 | arg=trueqdotr/197.323-2.0*TMath::Pi()*TMath::Floor(trueqdotr/(197.323*2.0*TMath::Pi())); |
c114f4f8 | 276 | cphi1=exp(fgkCI*arg); |
dd82cadc | 277 | cphis=fgkROOT2*real(cphi1); |
c114f4f8 | 278 | cphia=fgkCI*fgkROOT2*imag(cphi1); |
279 | corr0=real(fInteractionWsym*cphis*conj(cphis) | |
280 | +fInteractionWanti*cphia*conj(cphia) | |
281 | +fInteractionWnosym*cphi1*conj(cphi1)); | |
dd82cadc | 282 | goto OUTSIDE_INTERACTION_RANGE; |
283 | ||
284 | #ifdef REDUCED_MOM | |
c114f4f8 | 285 | kk=(int)TMath::Floor(trueqred/fInteractionDelk); |
286 | qred=(0.5+kk)*fInteractionDelk; | |
dd82cadc | 287 | #else |
c114f4f8 | 288 | kk=(int)TMath::Floor(2.0*trueqred/fInteractionDelk); |
289 | qred=(0.5+kk)*fInteractionDelk/2.0; | |
dd82cadc | 290 | #endif |
291 | qdotr=trueqdotr*qred/trueqred; | |
c114f4f8 | 292 | if(kk>=fInteractionNkmax){ |
dd82cadc | 293 | corr0=1.0; |
294 | goto OUTSIDE_INTERACTION_RANGE; | |
295 | } | |
296 | r=truer; | |
297 | ||
298 | eta=0.0; | |
3294e9ff | 299 | arg=qdotr/197.323-2.0*TMath::Pi()*TMath::Floor(qdotr/(197.323*2.0*TMath::Pi())); |
c114f4f8 | 300 | cphi1=exp(fgkCI*arg); |
dd82cadc | 301 | cphi2=conj(cphi1); |
302 | ||
303 | cphis=(cphi1+cphi2)/fgkROOT2; | |
304 | cphia=(cphi1-cphi2)/fgkROOT2; | |
305 | corr0=0.0; | |
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. */ | |
c114f4f8 | 309 | wsym_leftover=fInteractionWsym; |
310 | wanti_leftover=fInteractionWanti; | |
311 | wnosym_leftover=fInteractionWnosym; | |
dd82cadc | 312 | |
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: | |
317 | ||
318 | #ifdef BREIT_WIGNER | |
319 | corr0=corr0+bwcalc(trueqred,truer); | |
320 | #endif | |
321 | ||
322 | return corr0; | |
323 | } | |
324 | ||
8882f543 | 325 | #ifdef __DECCXX |
c114f4f8 | 326 | complex AliHBTCrab::CGamma(complex c){ |
8882f543 | 327 | #else |
c114f4f8 | 328 | double_complex AliHBTCrab::CGamma(double_complex c){ |
8882f543 | 329 | #endif |
dd82cadc | 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 | |
8882f543 | 334 | #ifdef __DECCXX |
335 | complex cg,cphase; | |
336 | #else | |
dd82cadc | 337 | double_complex cg,cphase; |
8882f543 | 338 | #endif |
dd82cadc | 339 | int mm,j; |
340 | double x,y,phase,delp,cgmag; | |
341 | x=real(c); | |
342 | y=imag(c); | |
343 | phase=-TMath::E()*y; | |
344 | for(j=1;j<=100000;j++){ | |
345 | delp=(y/(double)j)-atan(y/(double)j); | |
346 | phase=phase+delp; | |
c114f4f8 | 347 | if(TMath::Abs(delp)<1E-10) goto CGamma_ESCAPE; |
dd82cadc | 348 | } |
349 | printf("oops not accurate enough, increase jmax\n"); | |
c114f4f8 | 350 | CGamma_ESCAPE: |
3294e9ff | 351 | phase=phase-2.0*TMath::Pi()*TMath::Floor(phase/(2.0*TMath::Pi())); |
c114f4f8 | 352 | cphase=exp(fgkCI*phase); |
dd82cadc | 353 | cgmag=sqrt(TMath::Pi()*y/sinh(TMath::Pi()*y)); |
3294e9ff | 354 | mm=(int)TMath::Floor(x+0.5); |
dd82cadc | 355 | cg=cgmag*cphase; |
356 | if(mm<1){ | |
357 | for(j=1;j<=-mm+1;j++){ | |
c114f4f8 | 358 | cg=cg/(1.0+(double)(-j)+fgkCI*y); |
dd82cadc | 359 | } |
360 | } | |
361 | if(mm>1) { | |
362 | for(j=1;j<=mm-1;j++){ | |
c114f4f8 | 363 | cg=cg*((double)(j)+fgkCI*y); |
dd82cadc | 364 | } |
365 | } | |
366 | return cg; | |
367 | } |