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 }