Retrofeed from the Release developement
[u/mrichter/AliRoot.git] / HBTAN / AliHBTLLWeights.cxx
CommitLineData
4fdf4eb3 1#include "AliHBTLLWeights.h"
88cb7938 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)
c81f9591 86// and ializes various parameters.
88cb7938 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
1a1e58ac 93// SPI=TMath::Sqrt(PI)
88cb7938 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
4fdf4eb3 138
7f92929e 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_
1a1e58ac 147# define boosttoprf boosttoprf_
7f92929e 148# define fsiw fsiw_
88cb7938 149# define setpdist setpdist_
7f92929e 150# define type_of_call
151#else
152# define led_bldata LED_BLDATA
153# define fsiini FSIINI
154# define ltran12 LTRAN12
1a1e58ac 155# define boosttoprf BOOSTTOPRF
7f92929e 156# define fsiw FSIW
88cb7938 157# define setpdist SETPDIST
7f92929e 158# define type_of_call _stdcall
159#endif
160/****************************************************************/
161extern "C" void type_of_call led_bldata();
162extern "C" void type_of_call fsiini();
163extern "C" void type_of_call ltran12();
1a1e58ac 164extern "C" void type_of_call boosttoprf();
7f92929e 165extern "C" void type_of_call fsiw();
88cb7938 166extern "C" void type_of_call setpdist(Double_t& r);
7f92929e 167/**************************************************************/
168
88cb7938 169#include "AliHBTPair.h"
78d7c6d3 170#include "AliVAODParticle.h"
88cb7938 171#include "WLedCOMMONS.h"
88cb7938 172#include <TRandom.h>
173#include <TMath.h>
174#include <TPDGCode.h>
78d7c6d3 175#include <TParticlePDG.h>
176#include <TDatabasePDG.h>
88cb7938 177
178
7f92929e 179ClassImp(AliHBTLLWeights)
7f92929e 180
88cb7938 181AliHBTLLWeights* AliHBTLLWeights::fgLLWeights = 0x0;
182const Double_t AliHBTLLWeights::fgkWcons = 1./0.1973;
7f92929e 183
88cb7938 184AliHBTLLWeights::AliHBTLLWeights():
185 fTest(kTRUE),
186 fColoumbSwitch(kTRUE),
187 fQuantStatSwitch(kTRUE),
188 fStrongInterSwitch(kTRUE),
111e505b 189 fColWithResidNuclSwitch(kFALSE),
88cb7938 190 fNuclMass(0.0),
191 fNuclCharge(0.0),
1a1e58ac 192 fRandomPosition(kNone),
88cb7938 193 fRadius(0.0),
526c2bd5 194 fOneMinusLambda(0.0),
88cb7938 195 fPID1(0),
196 fPID2(0),
197 fSigma(0.0)
7f92929e 198{
88cb7938 199// Default Constructor
ae4a4473 200 if (fgLLWeights)
201 Fatal("AliHBTLLWeights","LLWeights already instatiated. Use AliHBTLLWeights::Instance()");
7f92929e 202}
88cb7938 203/**************************************************************/
7f92929e 204
111e505b 205AliHBTLLWeights::AliHBTLLWeights(const AliHBTLLWeights &/*source*/):
dd82cadc 206 AliHBTWeights(),
111e505b 207 fTest(kTRUE),
208 fColoumbSwitch(kTRUE),
209 fQuantStatSwitch(kTRUE),
210 fStrongInterSwitch(kTRUE),
211 fColWithResidNuclSwitch(kFALSE),
212 fNuclMass(0.0),
213 fNuclCharge(0.0),
1a1e58ac 214 fRandomPosition(kNone),
111e505b 215 fRadius(0.0),
526c2bd5 216 fOneMinusLambda(0.0),
111e505b 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
226AliHBTLLWeights& 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
4fdf4eb3 234AliHBTLLWeights* AliHBTLLWeights::Instance()
88cb7938 235{
236// returns instance of class
237 if (fgLLWeights)
238 {
4fdf4eb3 239 return fgLLWeights;
88cb7938 240 }
241 else
242 {
243 fgLLWeights = new AliHBTLLWeights();
244 return fgLLWeights;
245 }
246}
526c2bd5 247/************************************************************/
7f92929e 248
dd82cadc 249void 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
1a1e58ac 265Double_t AliHBTLLWeights::GetWeight(AliHBTPair* partpair)
7f92929e 266{
88cb7938 267// calculates weight for a pair
37e71815 268 static const Double_t kcmtofm = 1.e13;
269 static const Double_t kcmtoOneOverGeV = kcmtofm*fgkWcons;
88cb7938 270
78d7c6d3 271 AliVAODParticle *part1 = partpair->Particle1();
272 AliVAODParticle *part2 = partpair->Particle2();
2f8eea63 273
4fdf4eb3 274 if ( (part1 == 0x0) || (part2 == 0x0))
88cb7938 275 {
276 Error("GetWeight","Null particle pointer");
277 return 0.0;
278 }
279
72e72d00 280 if ( fPID1 != part1->GetPdgCode() ) return 1.0;
281 if ( fPID2 != part2->GetPdgCode() ) return 1.0;
88cb7938 282
526c2bd5 283//takes a lot of time
4fdf4eb3 284 if ( (part1->Px() == part2->Px()) &&
285 (part1->Py() == part2->Py()) &&
286 (part1->Pz() == part2->Pz()) )
88cb7938 287 {
288 return 0.0;
289 }
290
4fdf4eb3 291 if ((!fRandomPosition) &&
88cb7938 292 (part1->Vx() == part2->Vx()) &&
293 (part1->Vy() == part2->Vy()) &&
294 (part1->Vz() == part2->Vz()) )
2f8eea63 295 {
4fdf4eb3 296 return 0.0;
2f8eea63 297 }
88cb7938 298
526c2bd5 299 if(fOneMinusLambda)//implemetation of non-zero intetcept parameter
300 {
301 if( gRandom->Rndm() < fOneMinusLambda ) return 1.0;
302 }
303
7f92929e 304
88cb7938 305 //this must be after ltran12 because it would overwrite what we set below
1a1e58ac 306 switch (fRandomPosition)
88cb7938 307 {
1a1e58ac 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;
88cb7938 370 }
1a1e58ac 371
88cb7938 372 fsiw();
4fdf4eb3 373 return LEDWEIGHT.WEIN;
374}
7f92929e 375/************************************************************/
88cb7938 376
7f92929e 377void AliHBTLLWeights::Init()
4fdf4eb3 378{
88cb7938 379//initial parameters of model
380
4fdf4eb3 381 FSI_NS.NS = fApproximationModel;
382
88cb7938 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
4fdf4eb3 393 if ( (fPID1 == 0) || (fPID2 == 0) )
88cb7938 394 {
395 Fatal("Init","Particles types are not set");
396 return;//pro forma
397 }
c81f9591 398
399
4fdf4eb3 400 FSI_NS.LL = GetPairCode(fPID1,fPID2);
88cb7938 401
4fdf4eb3 402 if (FSI_NS.LL == 0)
88cb7938 403 {
404 Fatal("Init","Particles types are not supported");
405 return;//pro forma
406 }
407
c81f9591 408 Info("Init","Setting PIDs %d %d. LL Code is %d",fPID1,fPID2,FSI_NS.LL);
409
88cb7938 410
4fdf4eb3 411 TParticlePDG* tpart1 = TDatabasePDG::Instance()->GetParticle(fPID1);
412 if (tpart1 == 0x0)
88cb7938 413 {
414 Fatal("init","We can not find particle with ID=%d in PDG DataBase",fPID1);
415 return;
416 }
417
4fdf4eb3 418 FSI_POC.AM1=tpart1->Mass();
419 FSI_POC.C1=tpart1->Charge();
88cb7938 420
4fdf4eb3 421 TParticlePDG* tpart2 = TDatabasePDG::Instance()->GetParticle(fPID2);
88cb7938 422//lv
4fdf4eb3 423 if (tpart2 == 0x0)
88cb7938 424 {
425 Fatal("init","We can not find particle with ID=%d in our DataBase",fPID2);
426 return;
427 }
428
4fdf4eb3 429 FSI_POC.AM2=tpart2->Mass();
430 FSI_POC.C1=tpart2->Charge();
88cb7938 431
4fdf4eb3 432 led_bldata();
433 fsiini();
88cb7938 434
435
436//constants for radii simulation
437
438 if(fRandomPosition)
439 {
440 fSigma =TMath::Sqrt(2.)*fRadius;
441 }
7f92929e 442}
88cb7938 443/************************************************************/
7f92929e 444
445Int_t AliHBTLLWeights::GetPairCode(const AliHBTPair* partpair)
446{
88cb7938 447//returns Code corresponding to that pair
448 return GetPairCode(partpair->Particle1()->GetPdgCode(),partpair->Particle2()->GetPdgCode());
7f92929e 449}
88cb7938 450/************************************************************/
7f92929e 451
452Int_t AliHBTLLWeights::GetPairCode(Int_t pid1,Int_t pid2)
453{
88cb7938 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
7f92929e 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
88cb7938 469//determine the order of selcetion in switch
7f92929e 470 if (TMath::Abs(pid1) < TMath::Abs(pid2) )
88cb7938 471 {
472 if (pid1<0) chargefactor=-1;
473 hpid=pid2*chargefactor;
474 lpid=pid1*chargefactor;
475 swap = kFALSE;
476 }
7f92929e 477 else
88cb7938 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
7f92929e 491 switch (hpid) //switch on first particle id
88cb7938 492 {
493 case kNeutron:
7f92929e 494 switch (lpid)
88cb7938 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 }
7f92929e 512 break;
88cb7938 513
514 case kProton:
7f92929e 515 switch (lpid)
88cb7938 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 }
7f92929e 530 break;
88cb7938 531
532 case kPiPlus:
533
7f92929e 534 switch (lpid)
88cb7938 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 }
7f92929e 561 break;
88cb7938 562 case kPi0:
7f92929e 563 switch (lpid)
88cb7938 564 {
565 case kPi0:
566 code = 6;
567 break;
568
569 default:
570 return 0; //given pair not supported
571 break;
572 }
7f92929e 573 break;
574
88cb7938 575 case kKPlus:
7f92929e 576 switch (lpid)
88cb7938 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 }
7f92929e 594 break;
595
88cb7938 596 case kKMinus:
7f92929e 597 switch (lpid)
88cb7938 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 }
7f92929e 608 break;
609
88cb7938 610 case kK0:
7f92929e 611 switch (lpid)
88cb7938 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 }
7f92929e 625 break;
88cb7938 626
627 default: return 0;
628 }
7f92929e 629 return code;
630}
88cb7938 631/************************************************************/
632
633void AliHBTLLWeights::SetTest(Bool_t rtest)
634{
635 //Sets fTest member
636 fTest = rtest;
637}
638/************************************************************/
639
640void 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
647void 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
655void AliHBTLLWeights::SetStrongInterSwitch(Bool_t sis)
656{
657 //ISI: strong interaction between the two particles ON (OFF)
658 fStrongInterSwitch = sis;
659}
660/************************************************************/
7f92929e 661
88cb7938 662void AliHBTLLWeights::SetColWithResidNuclSwitch(Bool_t crn)
663{
664 //I3C: Coulomb interaction with residual nucleus ON (OFF)
665 fColWithResidNuclSwitch = crn;
666}
667/************************************************************/
668
669void AliHBTLLWeights::SetApproxModel(Int_t ap)
670{
671 //sets Model of Approximation (NS in Fortran code)
672 fApproximationModel=ap;
673}
674/************************************************************/
675
1a1e58ac 676void AliHBTLLWeights::SetRandomPosition(ERandomizationWay rw)
88cb7938 677{
1a1e58ac 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;
88cb7938 683}
684/************************************************************/
685
686void AliHBTLLWeights::SetR1dw(Double_t R)
687{
688 //spherical source model radii
689 fRadius=R;
690}
691/************************************************************/
692
693void AliHBTLLWeights::SetParticlesTypes(Int_t pid1, Int_t pid2)
694{
695 //set AliRoot particles types
696 fPID1 = pid1;
697 fPID2 = pid2;
698}
699/************************************************************/
700
701void AliHBTLLWeights::SetNucleusCharge(Double_t ch)
702{
703 // not used now (see comments in fortran code)
704 fNuclCharge=ch;
705}
706/************************************************************/
707
708void AliHBTLLWeights::SetNucleusMass(Double_t mass)
709{
710 // (see comments in fortran code)
711 fNuclMass=mass;
712}
1a1e58ac 713/************************************************************/
714
715Int_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
796void 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
817void 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
824void 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}