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