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