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