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