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  fApproximationModel(0),
 196  fPID1(0),
 197  fPID2(0),
 198  fSigma(0.0)
 199 {
 200 // Default Constructor
 201   if (fgLLWeights)
 202    Fatal("AliHBTLLWeights","LLWeights already instatiated. Use AliHBTLLWeights::Instance()");
 203 }
 204 /**************************************************************/
 205
 206 AliHBTLLWeights::AliHBTLLWeights(const AliHBTLLWeights &/*source*/):
 207  AliHBTWeights(),
 208  fTest(kTRUE),
 209  fColoumbSwitch(kTRUE),
 210  fQuantStatSwitch(kTRUE),
 211  fStrongInterSwitch(kTRUE),
 212  fColWithResidNuclSwitch(kFALSE),
 213  fNuclMass(0.0),
 214  fNuclCharge(0.0),
 215  fRandomPosition(kNone),
 216  fRadius(0.0),
 217  fOneMinusLambda(0.0),
 218  fApproximationModel(0),
 219  fPID1(0),
 220  fPID2(0),
 221  fSigma(0.0)
 222 {
 223   //Copy ctor needed by the coding conventions but not used
 224   Fatal("AliHBTLLWeights","copy ctor not implemented");
 225 }
 226 /************************************************************/
 227
 228 AliHBTLLWeights& AliHBTLLWeights::operator=(const AliHBTLLWeights& /*source*/)
 229 {
 230   //Assignment operator needed by the coding conventions but not used
 231   Fatal("AliHBTLLWeights","assignment operator not implemented");
 232   return * this;
 233 }
 234 /************************************************************/
 235
 236 AliHBTLLWeights* AliHBTLLWeights::Instance()
 237 {
 238 // returns instance of class
 239  if (fgLLWeights)
 240   {
 241     return fgLLWeights;
 242   }
 243  else
 244   {
 245    fgLLWeights = new AliHBTLLWeights();
 246    return fgLLWeights;
 247   }
 248 }
 249 /************************************************************/
 250
 251 void AliHBTLLWeights::Set()
 252 {
 253  //sets this as weighitng class
 254  Info("Set","Setting Lednicky-Lyuboshitz as Weighing Class");
 255
 256  if ( fgWeights == 0x0 )
 257   {
 258     fgWeights = AliHBTLLWeights::Instance();
 259     return;
 260   }
 261  if ( fgWeights == AliHBTLLWeights::Instance() ) return;
 262  delete fgWeights;
 263  fgWeights = AliHBTLLWeights::Instance();
 264 }
 265 /************************************************************/
 266
 267 Double_t AliHBTLLWeights::GetWeight(AliHBTPair* partpair)
 268 {
 269 // calculates weight for a pair
 270   static const Double_t kcmtofm = 1.e13;
 271   static const Double_t kcmtoOneOverGeV = kcmtofm*fgkWcons;
 272
 273   AliVAODParticle *part1 = partpair->Particle1();
 274   AliVAODParticle *part2 = partpair->Particle2();
 275
 276   if ( (part1 == 0x0) || (part2 == 0x0))
 277    {
 278      Error("GetWeight","Null particle pointer");
 279      return 0.0;
 280    }
 281
 282   if ( fPID1 != part1->GetPdgCode() ) return 1.0;
 283   if ( fPID2 != part2->GetPdgCode() ) return 1.0;
 284
 285 //takes a lot of time
 286   if ( (part1->Px() == part2->Px()) &&
 287        (part1->Py() == part2->Py()) &&
 288        (part1->Pz() == part2->Pz()) )
 289    {
 290      return 0.0;
 291    }
 292
 293   if ((!fRandomPosition) &&
 294       (part1->Vx()  == part2->Vx()) &&
 295       (part1->Vy()  == part2->Vy()) &&
 296       (part1->Vz()  == part2->Vz()) )
 297     {
 298       return 0.0;
 299     }
 300
 301   if(fOneMinusLambda)//implemetation of non-zero intetcept parameter
 302    {
 303      if( gRandom->Rndm() < fOneMinusLambda ) return 1.0;
 304    }
 305
 306
 307   //this must  be after ltran12 because it would overwrite what we set below
 308   switch (fRandomPosition)
 309    {
 310     case kGaussianQInv:
 311      {
 312        // Set Momenta in the common block
 313        FSI_MOM.P1X = part1->Px();
 314        FSI_MOM.P1Y = part1->Py();
 315        FSI_MOM.P1Z = part1->Pz();
 316
 317        FSI_MOM.P2X = part2->Px();
 318        FSI_MOM.P2Y = part2->Py();
 319        FSI_MOM.P2Z = part2->Pz();
 320
 321        //boost it to PRF
 322        ltran12();
 323        boosttoprf();
 324
 325        //Set particle positions now so they are Gaussian in PRF
 326        RandomPairDistances();
 327
 328        FSI_PRF.RPS= FSI_COOR.X2*FSI_COOR.X2 + FSI_COOR.Y2*FSI_COOR.Y2 +FSI_COOR.Z2*FSI_COOR.Z2;
 329        FSI_PRF.RP=TMath::Sqrt(FSI_PRF.RPS);
 330
 331        break;
 332       }
 333     case kGaussianOSL:
 334       {
 335        //boost pair to LCMS and set such a momenta in the common block
 336        Int_t retv = SetMomentaInLCMS(part1,part2);
 337        if (retv) return 1;
 338        //random particle positions/distance so they are Gaussian in LCMS
 339        RandomPairDistances();
 340
 341        //Boost to PRF
 342        ltran12();
 343        boosttoprf();
 344
 345        break;
 346       }
 347     case kNone:
 348     default:
 349        //set momenta and paricle positions as they are
 350        FSI_MOM.P1X = part1->Px();
 351        FSI_MOM.P1Y = part1->Py();
 352        FSI_MOM.P1Z = part1->Pz();
 353
 354        FSI_MOM.P2X = part2->Px();
 355        FSI_MOM.P2Y = part2->Py();
 356        FSI_MOM.P2Z = part2->Pz();
 357
 358        FSI_COOR.X1 = part1->Vx()*kcmtoOneOverGeV;
 359        FSI_COOR.Y1 = part1->Vy()*kcmtoOneOverGeV;
 360        FSI_COOR.Z1 = part1->Vz()*kcmtoOneOverGeV;
 361        FSI_COOR.T1 = part1->T()*kcmtoOneOverGeV;
 362
 363        FSI_COOR.X2 = part2->Vx()*kcmtoOneOverGeV;
 364        FSI_COOR.Y2 = part2->Vy()*kcmtoOneOverGeV;
 365        FSI_COOR.Z2 = part2->Vz()*kcmtoOneOverGeV;
 366        FSI_COOR.T2 = part2->T()*kcmtoOneOverGeV;
 367
 368        ltran12();
 369        boosttoprf();
 370
 371        break;
 372    }
 373
 374   fsiw();
 375   return LEDWEIGHT.WEIN;
 376 }
 377 /************************************************************/
 378
 379 void AliHBTLLWeights::Init()
 380 {
 381 //initial parameters of model
 382
 383   FSI_NS.NS = fApproximationModel;
 384
 385   LEDWEIGHT.ITEST = fTest;
 386   if(fTest)
 387    {
 388      FSI_NS.ICH = fColoumbSwitch;
 389      FSI_NS.ISI = fStrongInterSwitch;
 390      FSI_NS.IQS = fQuantStatSwitch;
 391      FSI_NS.I3C = fColWithResidNuclSwitch;
 392      LEDWEIGHT.IRANPOS = fRandomPosition;
 393    }
 394
 395   if ( (fPID1 == 0) || (fPID2 == 0) )
 396    {
 397      Fatal("Init","Particles types are not set");
 398      return;//pro forma
 399    }
 400
 401
 402   FSI_NS.LL = GetPairCode(fPID1,fPID2);
 403
 404   if (FSI_NS.LL == 0)
 405    {
 406      Fatal("Init","Particles types are not supported");
 407      return;//pro forma
 408    }
 409
 410   Info("Init","Setting PIDs %d %d. LL Code is %d",fPID1,fPID2,FSI_NS.LL);
 411
 412
 413   TParticlePDG* tpart1 = TDatabasePDG::Instance()->GetParticle(fPID1);
 414   if (tpart1 == 0x0)
 415    {
 416      Fatal("init","We can not find particle with ID=%d in PDG DataBase",fPID1);
 417      return;
 418    }
 419
 420   FSI_POC.AM1=tpart1->Mass();
 421   FSI_POC.C1=tpart1->Charge();
 422
 423   TParticlePDG* tpart2 = TDatabasePDG::Instance()->GetParticle(fPID2);
 424 //lv
 425   if (tpart2 == 0x0)
 426    {
 427      Fatal("init","We can not find particle with ID=%d in our DataBase",fPID2);
 428      return;
 429    }
 430
 431   FSI_POC.AM2=tpart2->Mass();
 432   FSI_POC.C1=tpart2->Charge();
 433
 434   led_bldata();
 435   fsiini();
 436
 437
 438 //constants for radii simulation
 439
 440   if(fRandomPosition)
 441    {
 442      fSigma =TMath::Sqrt(2.)*fRadius;
 443    }
 444 }
 445 /************************************************************/
 446
 447 Int_t AliHBTLLWeights::GetPairCode(const AliHBTPair* partpair)
 448 {
 449 //returns Code corresponding to that pair
 450  return GetPairCode(partpair->Particle1()->GetPdgCode(),partpair->Particle2()->GetPdgCode());
 451 }
 452 /************************************************************/
 453
 454 Int_t AliHBTLLWeights::GetPairCode(Int_t pid1,Int_t pid2)
 455 {
 456 // returns code corresponding to the pair of PIDs
 457 //   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
 458 //   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
 459 //   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
 460 //   NS=1 y/n:  +  +  +  +   +    -   -  -  -  -   -   -   -  -  -  -  -   -  -    -    -    -  -   - -  -      -      -
 461
 462 //alphas, deuterons and tyts are NOT supported here
 463
 464   Int_t chargefactor = 1;
 465   Int_t hpid; //pid in higher row
 466   Int_t lpid; //pid in lower row
 467   Int_t code; //pairCode
 468
 469   Bool_t swap;
 470
 471 //determine the order of selcetion in switch
 472   if (TMath::Abs(pid1) < TMath::Abs(pid2) )
 473    {
 474     if (pid1<0) chargefactor=-1;
 475     hpid=pid2*chargefactor;
 476     lpid=pid1*chargefactor;
 477     swap = kFALSE;
 478    }
 479   else
 480    {
 481     if (pid2<0) chargefactor=-1;
 482     hpid=pid1*chargefactor;
 483     lpid=pid2*chargefactor;
 484     swap = kTRUE;
 485    }
 486
 487 //mlv
 488    hpid=pid1;
 489    lpid=pid2;
 490
 491
 492 //Determine the pair code
 493   switch (hpid) //switch on first  particle id
 494    {
 495      case kNeutron:
 496       switch (lpid)
 497        {
 498          case kNeutron:
 499            code = 1;  //neutron neutron
 500            break;
 501
 502          case kProton:
 503            code = 3;  //neutron proton
 504            break;
 505
 506          case kLambda0:
 507            code = 28;  //neutron lambda
 508            break;
 509
 510          default:
 511            return 0; //given pair not supported
 512            break;
 513        }
 514       break;
 515
 516      case kProton:
 517       switch (lpid)
 518        {
 519          case kProton:
 520            code = 2; //proton proton
 521            break;
 522
 523          case kLambda0:
 524            code = 27;//proton lambda
 525            break;
 526
 527          default:
 528            return 0; //given pair not supported
 529            break;
 530
 531        }
 532       break;
 533
 534      case kPiPlus:
 535
 536       switch (lpid)
 537        {
 538          case kPiPlus:
 539            code = 7; //piplus piplus
 540            break;
 541
 542          case kPiMinus:
 543            code = 5; //piplus piminus
 544            break;
 545
 546          case kKMinus:
 547            code = 10; //piplus Kminus
 548            break;
 549
 550          case kKPlus:
 551            code = 11; //piplus Kplus
 552            break;
 553
 554          case kProton:
 555            code = 12; //piplus proton
 556            chargefactor*=-1;
 557            break;
 558
 559          default:
 560            return 0; //given pair not supported
 561            break;
 562        }
 563       break;
 564      case kPi0:
 565       switch (lpid)
 566        {
 567          case kPi0:
 568            code = 6;
 569            break;
 570
 571          default:
 572            return 0; //given pair not supported
 573            break;
 574        }
 575       break;
 576
 577      case kKPlus:
 578       switch (lpid)
 579        {
 580          case kKMinus:
 581            code = 14; //Kplus Kminus
 582            break;
 583
 584          case kKPlus:
 585            code = 15; //Kplus Kplus
 586            break;
 587
 588          case kProton:
 589            code = 16; //Kplus proton
 590            break;
 591
 592          default:
 593            return 0; //given pair not supported
 594            break;
 595        }
 596       break;
 597
 598      case kKMinus:
 599       switch (lpid)
 600        {
 601          case kProton:
 602            code = 17; //Kminus proton
 603            chargefactor*=1;
 604            break;
 605
 606          default:
 607            return 0; //given pair not supported
 608            break;
 609        }
 610       break;
 611
 612      case kK0:
 613       switch (lpid)
 614        {
 615          case kK0:
 616            code = 2; //Kzero Kzero
 617            break;
 618
 619          case kK0Bar:
 620            code = 17; //Kzero KzeroBar
 621            break;
 622
 623          default:
 624            return 0; //given pair not supported
 625            break;
 626        }
 627       break;
 628
 629      default: return 0;
 630    }
 631   return code;
 632 }
 633 /************************************************************/
 634
 635 void AliHBTLLWeights::SetTest(Bool_t rtest)
 636 {
 637   //Sets fTest member
 638   fTest = rtest;
 639 }
 640 /************************************************************/
 641
 642 void AliHBTLLWeights::SetColoumb(Bool_t col)
 643 {
 644   // (ICH in fortran code) Coulomb interaction between the two particles ON (OFF)
 645   fColoumbSwitch = col;
 646 }
 647 /************************************************************/
 648
 649 void AliHBTLLWeights::SetQuantumStatistics(Bool_t qss)
 650 {
 651   //IQS: quantum statistics for the two particles ON (OFF)
 652   //if non-identical particles automatically off
 653   fQuantStatSwitch = qss;
 654 }
 655 /************************************************************/
 656
 657 void AliHBTLLWeights::SetStrongInterSwitch(Bool_t sis)
 658 {
 659   //ISI: strong interaction between the two particles ON (OFF)
 660   fStrongInterSwitch = sis;
 661 }
 662 /************************************************************/
 663
 664 void AliHBTLLWeights::SetColWithResidNuclSwitch(Bool_t crn)
 665 {
 666   //I3C: Coulomb interaction with residual nucleus ON (OFF)
 667   fColWithResidNuclSwitch = crn;
 668 }
 669 /************************************************************/
 670
 671 void AliHBTLLWeights::SetApproxModel(Int_t ap)
 672 {
 673   //sets  Model of Approximation (NS in Fortran code)
 674   fApproximationModel=ap;
 675 }
 676 /************************************************************/
 677
 678 void AliHBTLLWeights::SetRandomPosition(ERandomizationWay rw)
 679 {
 680 // rw can be: kGaussianQInv - so 1D Qinv correlation function has a Gaussian shape
 681 //                            (and also Qlong and Qside, but not Qout)
 682 //            kGaussianOSL - so 3D Qout-Qside-Qlong correlation function has a Gaussian shape
 683 //            kNone  - no randomization performed (DEFAULT)
 684  fRandomPosition = rw;
 685 }
 686 /************************************************************/
 687
 688 void AliHBTLLWeights::SetR1dw(Double_t R)
 689 {
 690   //spherical source model radii
 691   fRadius=R;
 692 }
 693 /************************************************************/
 694
 695 void AliHBTLLWeights::SetParticlesTypes(Int_t pid1, Int_t pid2)
 696 {
 697   //set AliRoot particles types
 698   fPID1 = pid1;
 699   fPID2 = pid2;
 700 }
 701 /************************************************************/
 702
 703 void AliHBTLLWeights::SetNucleusCharge(Double_t ch)
 704 {
 705   // not used now  (see comments in fortran code)
 706   fNuclCharge=ch;
 707 }
 708 /************************************************************/
 709
 710 void AliHBTLLWeights::SetNucleusMass(Double_t mass)
 711 {
 712   // (see comments in fortran code)
 713   fNuclMass=mass;
 714 }
 715 /************************************************************/
 716
 717 Int_t AliHBTLLWeights::SetMomentaInLCMS(AliVAODParticle* part1, AliVAODParticle* part2)
 718 {
 719 //sets paricle momenta in the common block in LCMS frame
 720 //---> Particle energies ---------
 721
 722   Double_t p1x = part1->Px();
 723   Double_t p1y = part1->Py();
 724   Double_t p1z = part1->Pz();
 725
 726   Double_t p2x = part2->Px();
 727   Double_t p2y = part2->Py();
 728   Double_t p2z = part2->Pz();
 729
 730   Double_t am1 = part1->Mass();
 731   Double_t am2 = part2->Mass();
 732
 733   Double_t p1s=p1x*p1x+p1y*p1y+p1z*p1z;
 734   Double_t p2s=p2x*p2x+p2y*p2y+p2z*p2z;
 735   Double_t e1=TMath::Sqrt(am1*am1+p1s);
 736   Double_t e2=TMath::Sqrt(am2*am2+p2s);
 737 //---> pair parameters -----------
 738   Double_t e12=e1+e2;       // energy
 739   Double_t p12x=p1x+p2x;    // px
 740   Double_t p12y=p1y+p2y;    // py
 741   Double_t p12z=p1z+p2z;    // pz
 742   Double_t p12s=p12x*p12x+p12y*p12y+p12z*p12z;
 743   Double_t p12 =TMath::Sqrt(p12s);// momentum
 744   Double_t v12 =p12/e12;    // velocity
 745
 746   Double_t cth =p12z/p12;   // cos(theta)
 747 //  Double_t sth =TMath::Sqrt(1. - cth*cth); //sin
 748   Double_t v12z=v12*cth;    // longit. v
 749   Double_t gamz=1./TMath::Sqrt(1.-v12z*v12z);
 750 //--      v12t=v12*sth    // transv. v in cms (not needed)
 751
 752
 753   Double_t p12ts=p12x*p12x+p12y*p12y;
 754   Double_t p12t=TMath::Sqrt(p12ts); //pt
 755 //===> azimuthal rotation (pt||x) ============
 756    if(p12t != 0.0)
 757     {
 758      Double_t cphi=p12x/p12t;  // cos(phi)
 759      Double_t sphi=p12y/p12t;  // sin(phi)
 760      Rotate(p1x,p1y,sphi,cphi,p1x,p1y);
 761      Rotate(p2x,p2y,sphi,cphi,p2x,p2y);
 762 //     Rotate(x1,y1,sphi,cphi,x1,y1);
 763 //     Rotate(x2,y2,sphi,cphi,x2,y2);
 764      }
 765    else // rotation impossible
 766     {
 767       return 1;
 768     }
 769
 770 //===> co-moving ref. frame       ============
 771
 772    Double_t nothing;
 773    Boost(p1z,e1,v12z,gamz,p1z,nothing);
 774    Boost(p2z,e2,v12z,gamz,p2z,nothing);
 775
 776 //   p1s=p1x*p1x+p1y*p1y+p1z*p1z;
 777 //   p2s=p2x*p2x+p2y*p2y+p2z*p2z;
 778 //   e1=TMath::Sqrt(am1*am1+p1s);
 779 //   e2=TMath::Sqrt(am2*am2+p2s);
 780 //   Boost(z1,t1,v12z,gamz,z1,t1);
 781 //   Boost(z2,t2,v12z,gamz,z2,t2);
 782
 783
 784      FSI_MOM.P1X = p1x;
 785      FSI_MOM.P1Y = p1y;
 786      FSI_MOM.P1Z = p1z;
 787
 788      FSI_MOM.P2X = p2x;
 789      FSI_MOM.P2Y = p2y;
 790      FSI_MOM.P2Z = p2z;
 791
 792 //     Info("Simulate"," %f %f %f ",p1x + p2x, p1y + p2y, p1z + p2z);
 793      return 0;
 794 }
 795
 796 /**********************************************************/
 797
 798 void AliHBTLLWeights::RandomPairDistances()
 799 {
 800  //randomizes pair distances
 801   Double_t rxcm = fSigma*gRandom->Gaus();
 802   Double_t rycm = fSigma*gRandom->Gaus();
 803   Double_t rzcm = fSigma*gRandom->Gaus();
 804
 805
 806   FSI_COOR.X1 = 0;
 807   FSI_COOR.Y1 = 0;
 808   FSI_COOR.Z1 = 0;
 809   FSI_COOR.T1 = 0;
 810
 811   FSI_COOR.X2 = rxcm*fgkWcons;
 812   FSI_COOR.Y2 = rycm*fgkWcons;
 813   FSI_COOR.Z2 = rzcm*fgkWcons;
 814   FSI_COOR.T2 = 0;
 815
 816 }
 817
 818
 819 void AliHBTLLWeights::Rotate(Double_t x, Double_t y, Double_t sphi, Double_t cphi, Double_t& xr, Double_t& yr)
 820 {
 821 //rotates
 822  yr=y*cphi-x*sphi;
 823  xr=x*cphi+y*sphi;
 824 }
 825
 826 void AliHBTLLWeights::Boost(Double_t z, Double_t t, Double_t beta, Double_t gamma, Double_t& zt, Double_t& yt)
 827 {
 828 //boosts
 829   zt=gamma*(z-beta*t);
 830   yt=gamma*(t-beta*z);
 831 }