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