]>
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" |
78d7c6d3 | 13 | #include "AliVAODParticle.h" |
14 | #include "TDatabasePDG.h" | |
dd82cadc | 15 | #include <TMath.h> |
dd82cadc | 16 | |
17 | #include "volya_complex.h" | |
18 | ||
dd82cadc | 19 | AliHBTCrab* AliHBTCrab::fgCrab = 0x0; |
20 | ||
21 | const Double_t AliHBTCrab::fgkWcons = 1./0.1973; | |
22 | const Double_t AliHBTCrab::fgkROOT2=1.41421356237309504880; | |
8882f543 | 23 | #ifdef __DECCXX |
c114f4f8 | 24 | const complex AliHBTCrab::fgkCI(0.0,1.0); |
8882f543 | 25 | #else |
c031a417 | 26 | const doublecomplex AliHBTCrab::fgkCI(0.0,1.0); |
8882f543 | 27 | #endif |
dd82cadc | 28 | |
29 | /************************************************************/ | |
30 | ||
31 | AliHBTCrab* 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 | 42 | void 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 | ||
58 | AliHBTCrab::AliHBTCrab(): | |
59 | fBreitWigner(kFALSE), | |
60 | fReducedMom(kTRUE), | |
61 | fMaxMomentum(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 | |
71 | AliHBTCrab::AliHBTCrab(const AliHBTCrab &/*source*/): | |
72 | AliHBTWeights(), | |
73 | fBreitWigner(kFALSE), | |
74 | fReducedMom(kTRUE), | |
75 | fMaxMomentum(100.0) | |
76 | { | |
77 | //ctor | |
78 | } | |
79 | //=================================================================== | |
34914285 | 80 | AliHBTCrab & AliHBTCrab::operator=(const AliHBTCrab& /*source*/) |
c031a417 | 81 | { |
82 | //cpy constructor | |
83 | return *AliHBTCrab::Instance(); | |
84 | } | |
85 | ||
dd82cadc | 86 | void 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 | |
102 | Bool_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 | ||
116 | Double_t AliHBTCrab::GetWeight(const AliHBTPair* partpair) | |
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 | 135 | void 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 | 277 | double 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 | 336 | OUTSIDE_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 | 346 | complex 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 | 352 | doublecomplex 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 | 374 | CGamma_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 | } |