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