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