reverting to prev.
[u/mrichter/AliRoot.git] / HBTAN / AliHBTCrab.cxx
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 //                                                                     //
10 /////////////////////////////////////////////////////////////////////////
11  
12 #include "AliHBTPair.h"
13
14 #include <TMath.h>
15
16 #include "volya_complex.h"
17
18 AliHBTCrab* AliHBTCrab::fgCrab = 0x0;
19
20 const Double_t AliHBTCrab::fgkWcons = 1./0.1973;
21 const Double_t AliHBTCrab::fgkROOT2=1.41421356237309504880;
22 #ifdef __DECCXX
23 const complex AliHBTCrab::fgkCI(0.0,1.0);
24 #else
25 const doublecomplex AliHBTCrab::fgkCI(0.0,1.0);
26 #endif
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 //===================================================================
40
41 void AliHBTCrab::Set()
42 {
43 //sets this as weighitng class
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
63   if(fgCrab)
64    {
65      Fatal("AliHBTCrab","Do not use constructor directly. Use Instance() instead.");
66    }
67 }
68 //===================================================================
69
70 AliHBTCrab::AliHBTCrab(const AliHBTCrab &/*source*/):
71 AliHBTWeights(),
72 fBreitWigner(kFALSE),
73 fReducedMom(kTRUE),
74 fMaxMomentum(100.0)
75 {
76   //ctor
77 }
78 //===================================================================
79 AliHBTCrab & AliHBTCrab::operator=(const AliHBTCrab& /*source*/)
80 {
81 //cpy constructor
82  return *AliHBTCrab::Instance();
83 }
84
85 void AliHBTCrab::Init(Int_t pid1,Int_t pid2)
86 {
87 //Initialization method
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;
95      
96   fPid1 = pid1;   
97   fPid2 = pid2;
98 }
99 //===================================================================
100
101 Bool_t AliHBTCrab::SetConfig(const AliHBTPair* pair)
102 {
103 //returns the SetConfig
104   
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 {
117 //returns the weight
118   Double_t qred, r, qdotr, mom;
119   Int_t test;
120   
121   SetConfig(partpair);
122   
123   GetComQuantities(partpair, &qred, &r, &qdotr, &mom, &test);
124   if(test==0) 
125    {
126      Info("GetWeight","Test is 0");
127    }
128   Double_t corr = CorrCalc(qred,qdotr,r);
129   
130   return corr;
131 }
132 //===================================================================
133
134 void AliHBTCrab::GetComQuantities(const AliHBTPair* pair, 
135        double *qred,double *r,double *qdotr,double *mom, int *test)
136 {
137 //************************************
138 //  ALICE //
139
140  double p1[4]; 
141  double p2[4]; 
142  double r1[4]; 
143  double r2[4]; 
144  
145  static const Double_t kCmToFm = 1.e13;
146 // static const Double_t cmtoOneOverGeV = kCmToFm*fgkWcons;
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();
162  r1[1] = part1->Vx()*kCmToFm;
163  r1[2] = part1->Vy()*kCmToFm;
164  r1[3] = part1->Vz()*kCmToFm;
165
166  r2[0] = part2->T();
167  r2[1] = part2->Vx()*kCmToFm;
168  r2[2] = part2->Vy()*kCmToFm;
169  r2[3] = part2->Vz()*kCmToFm;
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;
177   double momtest;
178   if (fReducedMom)
179    {
180      momtest=4.0*fMaxMomentum*fMaxMomentum;
181    }  
182   else
183    {
184      momtest=fMaxMomentum*fMaxMomentum;
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   {
229   //  const double  kdotp=fMass2*fMass2-fMass1*fMass1;
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
276 double  AliHBTCrab::CorrCalc(double trueqred,double trueqdotr,double truer)
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
281   double eta,arg,corr0;
282 //  double xx,xxprime,xxjj,p1,zk;
283 //  int jj,kk,ipart,ipartcount,ispin;
284   int kk;
285   double wsymleftover,wantileftover,wnosymleftover;
286   double qred,qdotr,r;
287 //  const double rmass=fMass1*fMass2/(fMass1+fMass2);
288 #ifdef __DECCXX
289   complex cphi1,cphi2,cphis,cphia;
290 #else
291   doublecomplex cphi1,cphi2,cphis,cphia;
292 #endif
293
294   arg=trueqdotr/197.323-2.0*TMath::Pi()*TMath::Floor(trueqdotr/(197.323*2.0*TMath::Pi()));
295   cphi1=exp(fgkCI*arg);
296   cphis=fgkROOT2*real(cphi1);
297   cphia=fgkCI*fgkROOT2*imag(cphi1);
298   corr0=real(fInteractionWsym*cphis*conj(cphis)
299              +fInteractionWanti*cphia*conj(cphia)
300              +fInteractionWnosym*cphi1*conj(cphi1));
301   goto OUTSIDE_INTERACTION_RANGE;
302
303 #ifdef REDUCED_MOM
304   kk=(int)TMath::Floor(trueqred/fInteractionDelk);
305   qred=(0.5+kk)*fInteractionDelk;
306 #else
307   kk=(int)TMath::Floor(2.0*trueqred/fInteractionDelk);
308   qred=(0.5+kk)*fInteractionDelk/2.0;
309 #endif
310   qdotr=trueqdotr*qred/trueqred;
311   if(kk>=fInteractionNkmax){
312     corr0=1.0;
313     goto OUTSIDE_INTERACTION_RANGE;
314   }
315   r=truer;
316
317   eta=0.0;
318   arg=qdotr/197.323-2.0*TMath::Pi()*TMath::Floor(qdotr/(197.323*2.0*TMath::Pi()));
319   cphi1=exp(fgkCI*arg);
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. */
328   wsymleftover=fInteractionWsym;
329   wantileftover=fInteractionWanti;
330   wnosymleftover=fInteractionWnosym;
331
332   corr0=corr0+real(wsymleftover*cphis*conj(cphis)
333                    +wantileftover*cphia*conj(cphia)
334                    +wnosymleftover*cphi1*conj(cphi1));
335 OUTSIDE_INTERACTION_RANGE:
336
337 #ifdef BREIT_WIGNER
338   corr0=corr0+bwcalc(trueqred,truer);
339 #endif
340
341   return corr0;
342 }
343
344 #ifdef __DECCXX
345 complex AliHBTCrab::CGamma(complex c){
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
350 #else
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. */
354 // This code is written by Scott Pratt
355 // taken from http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html
356 #endif
357 #ifdef __DECCXX
358   complex cg,cphase;
359 #else
360   doublecomplex cg,cphase;
361 #endif
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;
370     if(TMath::Abs(delp)<1E-10) goto CGamma_ESCAPE;
371   }
372   printf("oops not accurate enough, increase jmax\n");
373 CGamma_ESCAPE:
374   phase=phase-2.0*TMath::Pi()*TMath::Floor(phase/(2.0*TMath::Pi()));
375   cphase=exp(fgkCI*phase);
376   cgmag=sqrt(TMath::Pi()*y/sinh(TMath::Pi()*y));
377   mm=(int)TMath::Floor(x+0.5);
378   cg=cgmag*cphase;
379   if(mm<1){
380     for(j=1;j<=-mm+1;j++){
381       cg=cg/(1.0+(double)(-j)+fgkCI*y);
382     }
383   }
384   if(mm>1) {
385     for(j=1;j<=mm-1;j++){
386       cg=cg*((double)(j)+fgkCI*y);
387     }
388   }
389   return cg;
390 }