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 }