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