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