HBTAN/AliHBTLLWeights.cxx

   1 #include "AliHBTLLWeights.h"
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 //_________________________________________________________________________
  18 ///////////////////////////////////////////////////////////////////////////
  19 //
  20 //  class AliHBTLLWeights
  21 //
  22 //  This class introduces the weight's calculation
  23 //  according to the Lednicky's algorithm.
  24 //
  25 //
  26 //  fsiw.f, fsiini.f
  27 //
  28 //  Description from fortran code by author R. Lednicky
  29 //
  30 //  Calculates final state interaction (FSI) weights
  31 //  WEIF = weight due to particle - (effective) nucleus FSI (p-N)
  32 //  WEI  = weight due to p-p-N FSI
  33 //  WEIN = weight due to p-p FSI; note that WEIN=WEI if I3C=0;
  34 //                                note that if I3C=1 the calculation of
  35 //                                WEIN can be skipped by putting J=0
  36 //.......................................................................
  37 //  Correlation Functions:
  38 //  CF(p-p-N)   = sum(WEI)/sum(WEIF)
  39 //  CF(p-p)     = sum(WEIN)/sum(1); here the nucleus is completely
  40 //                                  inactive
  41 //  CF(p-p-"N") = sum(WEIN*WEIF')/sum(WEIF'), where WEIN and WEIF'
  42 //                are not correlated (calculated at different emission
  43 //                points, e.g., for different events);
  44 //                thus here the nucleus affects one-particle
  45 //                spectra but not the correlation
  46 //.......................................................................
  47 //  User must supply data file <fn> on unit NUNIT (e.g. =11) specifying
  48 //  LL   : particle pair
  49 //  NS   : approximation used to calculate Bethe-Salpeter amplitude
  50 //  ITEST: test switch
  51 //         If ITEST=1 then also following parameters are required
  52 //  ICH  : 1(0) Coulomb interaction between the two particles ON (OFF)
  53 //  IQS  : 1(0) quantum statistics for the two particles ON (OFF)
  54 //  ISI  : 1(0) strong interaction between the two particles ON (OFF)
  55 //  I3C  : 1(0) Coulomb interaction with residual nucleus ON (OFF)
  56 //  This data file can contain other information useful for the user.
  57 //  It is read by subroutines READINT4 and READREA8(4) (or READ_FILE).
  58 //  -------------------------------------------------------------------
  59 //-   LL       1  2  3  4   5    6   7  8  9 10  11  12  13  14 15 16 17
  60 //-   part. 1: n  p  n alfa pi+ pi0 pi+ n  p pi+ pi+ pi+ pi- K+ K+ K+ K-
  61 //-   part. 2: n  p  p alfa pi- pi0 pi+ d  d  K-  K+  p   p  K- K+ p  p
  62 //   NS=1 y/n: +  +  +  +   +    -   -  -  -  -   -   -   -  -  -  -  -
  63 //  -------------------------------------------------------------------
  64 //-   LL       18  19 20  21  22 23 24 25 26    27     28
  65 //-   part. 1:  d  d   t  t   K0 K0  d p  p      p      n
  66 //-   part. 2:  d alfa t alfa K0 K0b t t alfa lambda lambda
  67 //   NS=1 y/n:  -  -   -  -   -  -   - -  -      +      +
  68 //  -------------------------------------------------------------------
  69 //   NS=1  Square well potential,
  70 //   NS=3  not used
  71 //   NS=4  scattered wave approximated by the spherical wave,
  72 //   NS=2  same as NS=4 but the approx. of equal emission times in PRF
  73 //       not required (t=0 approx. used in all other cases).
  74 //   Note: if NS=2,4, the B-S amplitude diverges at zero distance r* in
  75 //       the two-particle c.m.s.; user can specify a cutoff AA in
  76 //       SUBROUTINE FSIINI, for example:
  77 //       IF(NS.EQ.2.OR.NS.EQ.4)AA=5.D0 !! in 1/GeV --> AA=1. fm
  78 //  ------------------------------------------------------------------
  79 //  ITEST=1 any values of parameters ICH, IQS, ISI, I3C are allowed
  80 //          and should be given in data file <fn>
  81 //  ITEST=0 physical values of these parameters are put automatically
  82 //          in FSIINI (their values are not required in data file)
  83 //=====================================================================
  84 //  At the beginning of calculation user should call FSIINI,
  85 //  which reads LL, NS, ITEST (and eventually ICH, IQS, ISI, I3C)
  86 //  and ializes various parameters.
  87 //  In particular the constants in
  88 //    COMMON/FSI_CONS/PI,PI2,SPI,DR,W
  89 //  may be useful for the user:
  90 //   W=1/.1973D0    ! from fm to 1/GeV
  91 //   PI=4*DATAN(1.D0)
  92 //   PI2=2*PI
  93 //   SPI=TMath::Sqrt(PI)
  94 //   DR=180.D0/PI   ! from radian to degree
  95 //    _______________________________________________________
  96 //  !! |Important note: all real quantities are assumed REAL*8 | !!
  97 //    -------------------------------------------------------
  98 //  For each event user should fill in the following information
  99 //  in COMMONs (all COMMONs in FSI calculation start with FSI_):
 100 //  ...................................................................
 101 //   COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
 102 //  Only
 103 //       AMN  = mass of the effective nucleus   [GeV/c**2]
 104 //       CN   = charge of the effective nucleus [elem. charge units]
 105 //  are required
 106 //  ...................................................................
 107 //   COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, !part. momenta in the rest frame
 108 //  1               P2X,P2Y,P2Z,E2,P2  !of effective nucleus (NRF)
 109 //  Only the components
 110 //                      PiX,PiY,PiZ  [GeV/c]
 111 //  in NRF are required.
 112 //  To make the corresponding Lorentz transformation user can use the
 113 //  subroutines LTRAN and LTRANB
 114 //  ...................................................................
 115 //  COMMON/FSI_COOR/X1,Y1,Z1,T1,R1,     ! 4-coord. of emission
 116 //  1               X2,Y2,Z2,T2,R2      ! points in NRF
 117 //  The componets
 118 //                     Xi,Yi,Zi  [fm]
 119 //  and emission times
 120 //                        Ti   [fm/c]
 121 //  should be given in NRF with the origin assumed at the center
 122 //  of the effective nucleus. If the effect of residual nucleus is
 123 //  not calculated within FSIW, the NRF can be any fixed frame.
 124 //  --------------------------------------------------------------------
 125 //  Before calling FSIW the user must call
 126 //   CALL LTRAN12
 127 //  Besides Lorentz transformation to pair rest frame:
 128 //  (p1-p2)/2 --> k* it also transforms 4-coordinates of
 129 //  emission points from fm to 1/GeV and calculates Ei,Pi and Ri.
 130 //  Note that |k*|=AK in COMMON/FSI_PRF/
 131 //  --------------------------------------------------------------------
 132 //  After making some additional filtering using k* (say k* < k*max)
 133 //  or direction of vector k*,
 134 //  user can finally call FSIW to calculate the FSI weights
 135 //  to be used to construct the correlation function
 136 //======================================================================
 137
 138
 139 /*******************************************************************/
 140 /******      ROUTINES    USED    FOR     COMMUNUCATION      ********/
 141 /********************     WITH      FORTRAN     ********************/
 142 /*******************************************************************/
 143 #ifndef WIN32
 144 # define led_bldata led_bldata_
 145 # define fsiini fsiini_
 146 # define ltran12 ltran12_
 147 # define boosttoprf boosttoprf_
 148 # define fsiw fsiw_
 149 # define setpdist setpdist_
 150 # define type_of_call
 151 #else
 152 # define led_bldata LED_BLDATA
 153 # define fsiini FSIINI
 154 # define ltran12 LTRAN12
 155 # define boosttoprf BOOSTTOPRF
 156 # define fsiw FSIW
 157 # define setpdist SETPDIST
 158 # define type_of_call _stdcall
 159 #endif
 160 /****************************************************************/
 161 extern "C" void type_of_call led_bldata();
 162 extern "C" void type_of_call fsiini();
 163 extern "C" void type_of_call ltran12();
 164 extern "C" void type_of_call boosttoprf();
 165 extern "C" void type_of_call fsiw();
 166 extern "C" void type_of_call setpdist(Double_t& r);
 167 /**************************************************************/
 168
 169 #include "AliHBTPair.h"
 170 #include "AliVAODParticle.h"
 171 #include "WLedCOMMONS.h"
 172 #include <TRandom.h>
 173 #include <TMath.h>
 174 #include <TPDGCode.h>
 175 #include <TParticlePDG.h>
 176 #include <TDatabasePDG.h>
 177
 178
 179 ClassImp(AliHBTLLWeights)
 180
 181 AliHBTLLWeights* AliHBTLLWeights::fgLLWeights = 0x0;
 182 const Double_t AliHBTLLWeights::fgkWcons = 1./0.1973;
 183
 184 AliHBTLLWeights::AliHBTLLWeights():
 185  fTest(kTRUE),
 186  fColoumbSwitch(kTRUE),
 187  fQuantStatSwitch(kTRUE),
 188  fStrongInterSwitch(kTRUE),
 189  fColWithResidNuclSwitch(kFALSE),
 190  fNuclMass(0.0),
 191  fNuclCharge(0.0),
 192  fRandomPosition(kNone),
 193  fRadius(0.0),
 194  fOneMinusLambda(0.0),
 195  fPID1(0),
 196  fPID2(0),
 197  fSigma(0.0)
 198 {
 199 // Default Constructor
 200   if (fgLLWeights)
 201    Fatal("AliHBTLLWeights","LLWeights already instatiated. Use AliHBTLLWeights::Instance()");
 202 }
 203 /**************************************************************/
 204
 205 AliHBTLLWeights::AliHBTLLWeights(const AliHBTLLWeights &/*source*/):
 206  AliHBTWeights(),
 207  fTest(kTRUE),
 208  fColoumbSwitch(kTRUE),
 209  fQuantStatSwitch(kTRUE),
 210  fStrongInterSwitch(kTRUE),
 211  fColWithResidNuclSwitch(kFALSE),
 212  fNuclMass(0.0),
 213  fNuclCharge(0.0),
 214  fRandomPosition(kNone),
 215  fRadius(0.0),
 216  fOneMinusLambda(0.0),
 217  fPID1(0),
 218  fPID2(0),
 219  fSigma(0.0)
 220 {
 221   //Copy ctor needed by the coding conventions but not used
 222   Fatal("AliHBTLLWeights","copy ctor not implemented");
 223 }
 224 /************************************************************/
 225
 226 AliHBTLLWeights& AliHBTLLWeights::operator=(const AliHBTLLWeights& /*source*/)
 227 {
 228   //Assignment operator needed by the coding conventions but not used
 229   Fatal("AliHBTLLWeights","assignment operator not implemented");
 230   return * this;
 231 }
 232 /************************************************************/
 233
 234 AliHBTLLWeights* AliHBTLLWeights::Instance()
 235 {
 236 // returns instance of class
 237  if (fgLLWeights)
 238   {
 239     return fgLLWeights;
 240   }
 241  else
 242   {
 243    fgLLWeights = new AliHBTLLWeights();
 244    return fgLLWeights;
 245   }
 246 }
 247 /************************************************************/
 248
 249 void AliHBTLLWeights::Set()
 250 {
 251  //sets this as weighitng class
 252  Info("Set","Setting Lednicky-Lyuboshitz as Weighing Class");
 253
 254  if ( fgWeights == 0x0 )
 255   {
 256     fgWeights = AliHBTLLWeights::Instance();
 257     return;
 258   }
 259  if ( fgWeights == AliHBTLLWeights::Instance() ) return;
 260  delete fgWeights;
 261  fgWeights = AliHBTLLWeights::Instance();
 262 }
 263 /************************************************************/
 264
 265 Double_t AliHBTLLWeights::GetWeight(AliHBTPair* partpair)
 266 {
 267 // calculates weight for a pair
 268   static const Double_t kcmtofm = 1.e13;
 269   static const Double_t kcmtoOneOverGeV = kcmtofm*fgkWcons;
 270
 271   AliVAODParticle *part1 = partpair->Particle1();
 272   AliVAODParticle *part2 = partpair->Particle2();
 273
 274   if ( (part1 == 0x0) || (part2 == 0x0))
 275    {
 276      Error("GetWeight","Null particle pointer");
 277      return 0.0;
 278    }
 279
 280   if ( fPID1 != part1->GetPdgCode() ) return 1.0;
 281   if ( fPID2 != part2->GetPdgCode() ) return 1.0;
 282
 283 //takes a lot of time
 284   if ( (part1->Px() == part2->Px()) &&
 285        (part1->Py() == part2->Py()) &&
 286        (part1->Pz() == part2->Pz()) )
 287    {
 288      return 0.0;
 289    }
 290
 291   if ((!fRandomPosition) &&
 292       (part1->Vx()  == part2->Vx()) &&
 293       (part1->Vy()  == part2->Vy()) &&
 294       (part1->Vz()  == part2->Vz()) )
 295     {
 296       return 0.0;
 297     }
 298
 299   if(fOneMinusLambda)//implemetation of non-zero intetcept parameter
 300    {
 301      if( gRandom->Rndm() < fOneMinusLambda ) return 1.0;
 302    }
 303
 304
 305   //this must  be after ltran12 because it would overwrite what we set below
 306   switch (fRandomPosition)
 307    {
 308     case kGaussianQInv:
 309      {
 310        // Set Momenta in the common block
 311        FSI_MOM.P1X = part1->Px();
 312        FSI_MOM.P1Y = part1->Py();
 313        FSI_MOM.P1Z = part1->Pz();
 314
 315        FSI_MOM.P2X = part2->Px();
 316        FSI_MOM.P2Y = part2->Py();
 317        FSI_MOM.P2Z = part2->Pz();
 318
 319        //boost it to PRF
 320        ltran12();
 321        boosttoprf();
 322
 323        //Set particle positions now so they are Gaussian in PRF
 324        RandomPairDistances();
 325
 326        FSI_PRF.RPS= FSI_COOR.X2*FSI_COOR.X2 + FSI_COOR.Y2*FSI_COOR.Y2 +FSI_COOR.Z2*FSI_COOR.Z2;
 327        FSI_PRF.RP=TMath::Sqrt(FSI_PRF.RPS);
 328
 329        break;
 330       }
 331     case kGaussianOSL:
 332       {
 333        //boost pair to LCMS and set such a momenta in the common block
 334        Int_t retv = SetMomentaInLCMS(part1,part2);
 335        if (retv) return 1;
 336        //random particle positions/distance so they are Gaussian in LCMS
 337        RandomPairDistances();
 338
 339        //Boost to PRF
 340        ltran12();
 341        boosttoprf();
 342
 343        break;
 344       }
 345     case kNone:
 346     default:
 347        //set momenta and paricle positions as they are
 348        FSI_MOM.P1X = part1->Px();
 349        FSI_MOM.P1Y = part1->Py();
 350        FSI_MOM.P1Z = part1->Pz();
 351
 352        FSI_MOM.P2X = part2->Px();
 353        FSI_MOM.P2Y = part2->Py();
 354        FSI_MOM.P2Z = part2->Pz();
 355
 356        FSI_COOR.X1 = part1->Vx()*kcmtoOneOverGeV;
 357        FSI_COOR.Y1 = part1->Vy()*kcmtoOneOverGeV;
 358        FSI_COOR.Z1 = part1->Vz()*kcmtoOneOverGeV;
 359        FSI_COOR.T1 = part1->T()*kcmtoOneOverGeV;
 360
 361        FSI_COOR.X2 = part2->Vx()*kcmtoOneOverGeV;
 362        FSI_COOR.Y2 = part2->Vy()*kcmtoOneOverGeV;
 363        FSI_COOR.Z2 = part2->Vz()*kcmtoOneOverGeV;
 364        FSI_COOR.T2 = part2->T()*kcmtoOneOverGeV;
 365
 366        ltran12();
 367        boosttoprf();
 368
 369        break;
 370    }
 371
 372   fsiw();
 373   return LEDWEIGHT.WEIN;
 374 }
 375 /************************************************************/
 376
 377 void AliHBTLLWeights::Init()
 378 {
 379 //initial parameters of model
 380
 381   FSI_NS.NS = fApproximationModel;
 382
 383   LEDWEIGHT.ITEST = fTest;
 384   if(fTest)
 385    {
 386      FSI_NS.ICH = fColoumbSwitch;
 387      FSI_NS.ISI = fStrongInterSwitch;
 388      FSI_NS.IQS = fQuantStatSwitch;
 389      FSI_NS.I3C = fColWithResidNuclSwitch;
 390      LEDWEIGHT.IRANPOS = fRandomPosition;
 391    }
 392
 393   if ( (fPID1 == 0) || (fPID2 == 0) )
 394    {
 395      Fatal("Init","Particles types are not set");
 396      return;//pro forma
 397    }
 398
 399
 400   FSI_NS.LL = GetPairCode(fPID1,fPID2);
 401
 402   if (FSI_NS.LL == 0)
 403    {
 404      Fatal("Init","Particles types are not supported");
 405      return;//pro forma
 406    }
 407
 408   Info("Init","Setting PIDs %d %d. LL Code is %d",fPID1,fPID2,FSI_NS.LL);
 409
 410
 411   TParticlePDG* tpart1 = TDatabasePDG::Instance()->GetParticle(fPID1);
 412   if (tpart1 == 0x0)
 413    {
 414      Fatal("init","We can not find particle with ID=%d in PDG DataBase",fPID1);
 415      return;
 416    }
 417
 418   FSI_POC.AM1=tpart1->Mass();
 419   FSI_POC.C1=tpart1->Charge();
 420
 421   TParticlePDG* tpart2 = TDatabasePDG::Instance()->GetParticle(fPID2);
 422 //lv
 423   if (tpart2 == 0x0)
 424    {
 425      Fatal("init","We can not find particle with ID=%d in our DataBase",fPID2);
 426      return;
 427    }
 428
 429   FSI_POC.AM2=tpart2->Mass();
 430   FSI_POC.C1=tpart2->Charge();
 431
 432   led_bldata();
 433   fsiini();
 434
 435
 436 //constants for radii simulation
 437
 438   if(fRandomPosition)
 439    {
 440      fSigma =TMath::Sqrt(2.)*fRadius;
 441    }
 442 }
 443 /************************************************************/
 444
 445 Int_t AliHBTLLWeights::GetPairCode(const AliHBTPair* partpair)
 446 {
 447 //returns Code corresponding to that pair
 448  return GetPairCode(partpair->Particle1()->GetPdgCode(),partpair->Particle2()->GetPdgCode());
 449 }
 450 /************************************************************/
 451
 452 Int_t AliHBTLLWeights::GetPairCode(Int_t pid1,Int_t pid2)
 453 {
 454 // returns code corresponding to the pair of PIDs
 455 //   pairCode   1  2  3  4   5    6   7  8  9 10  11  12  13  14 15 16 17 18  19  20   21   22  23 24 25 26    27     28
 456 //   hpid:      n  p  n alfa pi+ pi0 pi+ n  p pi+ pi+ pi+ pi- K+ K+ K+ K-  d  d    t   t    K0  K0  d p  p      p      n
 457 //   lpid:      n  p  p alfa pi- pi0 pi+ d  d  K-  K+  p   p  K- K+ p  p   d alfa  t  alfa  K0  K0b t t alfa lambda lambda
 458 //   NS=1 y/n:  +  +  +  +   +    -   -  -  -  -   -   -   -  -  -  -  -   -  -    -    -    -  -   - -  -      -      -
 459
 460 //alphas, deuterons and tyts are NOT supported here
 461
 462   Int_t chargefactor = 1;
 463   Int_t hpid; //pid in higher row
 464   Int_t lpid; //pid in lower row
 465   Int_t code; //pairCode
 466
 467   Bool_t swap;
 468
 469 //determine the order of selcetion in switch
 470   if (TMath::Abs(pid1) < TMath::Abs(pid2) )
 471    {
 472     if (pid1<0) chargefactor=-1;
 473     hpid=pid2*chargefactor;
 474     lpid=pid1*chargefactor;
 475     swap = kFALSE;
 476    }
 477   else
 478    {
 479     if (pid2<0) chargefactor=-1;
 480     hpid=pid1*chargefactor;
 481     lpid=pid2*chargefactor;
 482     swap = kTRUE;
 483    }
 484
 485 //mlv
 486    hpid=pid1;
 487    lpid=pid2;
 488
 489
 490 //Determine the pair code
 491   switch (hpid) //switch on first  particle id
 492    {
 493      case kNeutron:
 494       switch (lpid)
 495        {
 496          case kNeutron:
 497            code = 1;  //neutron neutron
 498            break;
 499
 500          case kProton:
 501            code = 3;  //neutron proton
 502            break;
 503
 504          case kLambda0:
 505            code = 28;  //neutron lambda
 506            break;
 507
 508          default:
 509            return 0; //given pair not supported
 510            break;
 511        }
 512       break;
 513
 514      case kProton:
 515       switch (lpid)
 516        {
 517          case kProton:
 518            code = 2; //proton proton
 519            break;
 520
 521          case kLambda0:
 522            code = 27;//proton lambda
 523            break;
 524
 525          default:
 526            return 0; //given pair not supported
 527            break;
 528
 529        }
 530       break;
 531
 532      case kPiPlus:
 533
 534       switch (lpid)
 535        {
 536          case kPiPlus:
 537            code = 7; //piplus piplus
 538            break;
 539
 540          case kPiMinus:
 541            code = 5; //piplus piminus
 542            break;
 543
 544          case kKMinus:
 545            code = 10; //piplus Kminus
 546            break;
 547
 548          case kKPlus:
 549            code = 11; //piplus Kplus
 550            break;
 551
 552          case kProton:
 553            code = 12; //piplus proton
 554            chargefactor*=-1;
 555            break;
 556
 557          default:
 558            return 0; //given pair not supported
 559            break;
 560        }
 561       break;
 562      case kPi0:
 563       switch (lpid)
 564        {
 565          case kPi0:
 566            code = 6;
 567            break;
 568
 569          default:
 570            return 0; //given pair not supported
 571            break;
 572        }
 573       break;
 574
 575      case kKPlus:
 576       switch (lpid)
 577        {
 578          case kKMinus:
 579            code = 14; //Kplus Kminus
 580            break;
 581
 582          case kKPlus:
 583            code = 15; //Kplus Kplus
 584            break;
 585
 586          case kProton:
 587            code = 16; //Kplus proton
 588            break;
 589
 590          default:
 591            return 0; //given pair not supported
 592            break;
 593        }
 594       break;
 595
 596      case kKMinus:
 597       switch (lpid)
 598        {
 599          case kProton:
 600            code = 17; //Kminus proton
 601            chargefactor*=1;
 602            break;
 603
 604          default:
 605            return 0; //given pair not supported
 606            break;
 607        }
 608       break;
 609
 610      case kK0:
 611       switch (lpid)
 612        {
 613          case kK0:
 614            code = 2; //Kzero Kzero
 615            break;
 616
 617          case kK0Bar:
 618            code = 17; //Kzero KzeroBar
 619            break;
 620
 621          default:
 622            return 0; //given pair not supported
 623            break;
 624        }
 625       break;
 626
 627      default: return 0;
 628    }
 629   return code;
 630 }
 631 /************************************************************/
 632
 633 void AliHBTLLWeights::SetTest(Bool_t rtest)
 634 {
 635   //Sets fTest member
 636   fTest = rtest;
 637 }
 638 /************************************************************/
 639
 640 void AliHBTLLWeights::SetColoumb(Bool_t col)
 641 {
 642   // (ICH in fortran code) Coulomb interaction between the two particles ON (OFF)
 643   fColoumbSwitch = col;
 644 }
 645 /************************************************************/
 646
 647 void AliHBTLLWeights::SetQuantumStatistics(Bool_t qss)
 648 {
 649   //IQS: quantum statistics for the two particles ON (OFF)
 650   //if non-identical particles automatically off
 651   fQuantStatSwitch = qss;
 652 }
 653 /************************************************************/
 654
 655 void AliHBTLLWeights::SetStrongInterSwitch(Bool_t sis)
 656 {
 657   //ISI: strong interaction between the two particles ON (OFF)
 658   fStrongInterSwitch = sis;
 659 }
 660 /************************************************************/
 661
 662 void AliHBTLLWeights::SetColWithResidNuclSwitch(Bool_t crn)
 663 {
 664   //I3C: Coulomb interaction with residual nucleus ON (OFF)
 665   fColWithResidNuclSwitch = crn;
 666 }
 667 /************************************************************/
 668
 669 void AliHBTLLWeights::SetApproxModel(Int_t ap)
 670 {
 671   //sets  Model of Approximation (NS in Fortran code)
 672   fApproximationModel=ap;
 673 }
 674 /************************************************************/
 675
 676 void AliHBTLLWeights::SetRandomPosition(ERandomizationWay rw)
 677 {
 678 // rw can be: kGaussianQInv - so 1D Qinv correlation function has a Gaussian shape
 679 //                            (and also Qlong and Qside, but not Qout)
 680 //            kGaussianOSL - so 3D Qout-Qside-Qlong correlation function has a Gaussian shape
 681 //            kNone  - no randomization performed (DEFAULT)
 682  fRandomPosition = rw;
 683 }
 684 /************************************************************/
 685
 686 void AliHBTLLWeights::SetR1dw(Double_t R)
 687 {
 688   //spherical source model radii
 689   fRadius=R;
 690 }
 691 /************************************************************/
 692
 693 void AliHBTLLWeights::SetParticlesTypes(Int_t pid1, Int_t pid2)
 694 {
 695   //set AliRoot particles types
 696   fPID1 = pid1;
 697   fPID2 = pid2;
 698 }
 699 /************************************************************/
 700
 701 void AliHBTLLWeights::SetNucleusCharge(Double_t ch)
 702 {
 703   // not used now  (see comments in fortran code)
 704   fNuclCharge=ch;
 705 }
 706 /************************************************************/
 707
 708 void AliHBTLLWeights::SetNucleusMass(Double_t mass)
 709 {
 710   // (see comments in fortran code)
 711   fNuclMass=mass;
 712 }
 713 /************************************************************/
 714
 715 Int_t AliHBTLLWeights::SetMomentaInLCMS(AliVAODParticle* part1, AliVAODParticle* part2)
 716 {
 717 //sets paricle momenta in the common block in LCMS frame
 718 //---> Particle energies ---------
 719
 720   Double_t p1x = part1->Px();
 721   Double_t p1y = part1->Py();
 722   Double_t p1z = part1->Pz();
 723
 724   Double_t p2x = part2->Px();
 725   Double_t p2y = part2->Py();
 726   Double_t p2z = part2->Pz();
 727
 728   Double_t am1 = part1->Mass();
 729   Double_t am2 = part2->Mass();
 730
 731   Double_t p1s=p1x*p1x+p1y*p1y+p1z*p1z;
 732   Double_t p2s=p2x*p2x+p2y*p2y+p2z*p2z;
 733   Double_t e1=TMath::Sqrt(am1*am1+p1s);
 734   Double_t e2=TMath::Sqrt(am2*am2+p2s);
 735 //---> pair parameters -----------
 736   Double_t e12=e1+e2;       // energy
 737   Double_t p12x=p1x+p2x;    // px
 738   Double_t p12y=p1y+p2y;    // py
 739   Double_t p12z=p1z+p2z;    // pz
 740   Double_t p12s=p12x*p12x+p12y*p12y+p12z*p12z;
 741   Double_t p12 =TMath::Sqrt(p12s);// momentum
 742   Double_t v12 =p12/e12;    // velocity
 743
 744   Double_t cth =p12z/p12;   // cos(theta)
 745 //  Double_t sth =TMath::Sqrt(1. - cth*cth); //sin
 746   Double_t v12z=v12*cth;    // longit. v
 747   Double_t gamz=1./TMath::Sqrt(1.-v12z*v12z);
 748 //--      v12t=v12*sth    // transv. v in cms (not needed)
 749
 750
 751   Double_t p12ts=p12x*p12x+p12y*p12y;
 752   Double_t p12t=TMath::Sqrt(p12ts); //pt
 753 //===> azimuthal rotation (pt||x) ============
 754    if(p12t != 0.0)
 755     {
 756      Double_t cphi=p12x/p12t;  // cos(phi)
 757      Double_t sphi=p12y/p12t;  // sin(phi)
 758      Rotate(p1x,p1y,sphi,cphi,p1x,p1y);
 759      Rotate(p2x,p2y,sphi,cphi,p2x,p2y);
 760 //     Rotate(x1,y1,sphi,cphi,x1,y1);
 761 //     Rotate(x2,y2,sphi,cphi,x2,y2);
 762      }
 763    else // rotation impossible
 764     {
 765       return 1;
 766     }
 767
 768 //===> co-moving ref. frame       ============
 769
 770    Double_t nothing;
 771    Boost(p1z,e1,v12z,gamz,p1z,nothing);
 772    Boost(p2z,e2,v12z,gamz,p2z,nothing);
 773
 774 //   p1s=p1x*p1x+p1y*p1y+p1z*p1z;
 775 //   p2s=p2x*p2x+p2y*p2y+p2z*p2z;
 776 //   e1=TMath::Sqrt(am1*am1+p1s);
 777 //   e2=TMath::Sqrt(am2*am2+p2s);
 778 //   Boost(z1,t1,v12z,gamz,z1,t1);
 779 //   Boost(z2,t2,v12z,gamz,z2,t2);
 780
 781
 782      FSI_MOM.P1X = p1x;
 783      FSI_MOM.P1Y = p1y;
 784      FSI_MOM.P1Z = p1z;
 785
 786      FSI_MOM.P2X = p2x;
 787      FSI_MOM.P2Y = p2y;
 788      FSI_MOM.P2Z = p2z;
 789
 790 //     Info("Simulate"," %f %f %f ",p1x + p2x, p1y + p2y, p1z + p2z);
 791      return 0;
 792 }
 793
 794 /**********************************************************/
 795
 796 void AliHBTLLWeights::RandomPairDistances()
 797 {
 798  //randomizes pair distances
 799   Double_t rxcm = fSigma*gRandom->Gaus();
 800   Double_t rycm = fSigma*gRandom->Gaus();
 801   Double_t rzcm = fSigma*gRandom->Gaus();
 802
 803
 804   FSI_COOR.X1 = 0;
 805   FSI_COOR.Y1 = 0;
 806   FSI_COOR.Z1 = 0;
 807   FSI_COOR.T1 = 0;
 808
 809   FSI_COOR.X2 = rxcm*fgkWcons;
 810   FSI_COOR.Y2 = rycm*fgkWcons;
 811   FSI_COOR.Z2 = rzcm*fgkWcons;
 812   FSI_COOR.T2 = 0;
 813
 814 }
 815
 816
 817 void AliHBTLLWeights::Rotate(Double_t x, Double_t y, Double_t sphi, Double_t cphi, Double_t& xr, Double_t& yr)
 818 {
 819 //rotates
 820  yr=y*cphi-x*sphi;
 821  xr=x*cphi+y*sphi;
 822 }
 823
 824 void AliHBTLLWeights::Boost(Double_t z, Double_t t, Double_t beta, Double_t gamma, Double_t& zt, Double_t& yt)
 825 {
 826 //boosts
 827   zt=gamma*(z-beta*t);
 828   yt=gamma*(t-beta*z);
 829 }