]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HBTAN/AliHBTLLWeights.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Working version (P.Skowronski)
[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)
86// and initializes 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=DSQRT(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
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_
147# define fsiw fsiw_
88cb7938 148# define setpdist setpdist_
7f92929e 149# define type_of_call
150#else
151# define led_bldata LED_BLDATA
152# define fsiini FSIINI
153# define ltran12 LTRAN12
154# define fsiw FSIW
88cb7938 155# define setpdist SETPDIST
7f92929e 156# define type_of_call _stdcall
157#endif
158/****************************************************************/
159extern "C" void type_of_call led_bldata();
160extern "C" void type_of_call fsiini();
161extern "C" void type_of_call ltran12();
162extern "C" void type_of_call fsiw();
88cb7938 163extern "C" void type_of_call setpdist(Double_t& r);
7f92929e 164/**************************************************************/
165
88cb7938 166#include "AliHBTPair.h"
167#include "AliHBTParticle.h"
168#include "WLedCOMMONS.h"
169#include <TList.h>
170#include <TRandom.h>
171#include <TMath.h>
172#include <TPDGCode.h>
173
174
7f92929e 175ClassImp(AliHBTLLWeights)
7f92929e 176
88cb7938 177AliHBTLLWeights* AliHBTLLWeights::fgLLWeights = 0x0;
178const Double_t AliHBTLLWeights::fgkWcons = 1./0.1973;
7f92929e 179
88cb7938 180AliHBTLLWeights::AliHBTLLWeights():
181 fTest(kTRUE),
182 fColoumbSwitch(kTRUE),
183 fQuantStatSwitch(kTRUE),
184 fStrongInterSwitch(kTRUE),
111e505b 185 fColWithResidNuclSwitch(kFALSE),
88cb7938 186 fNuclMass(0.0),
187 fNuclCharge(0.0),
188 fRandomPosition(kFALSE),
189 fRadius(0.0),
190 fPID1(0),
191 fPID2(0),
192 fSigma(0.0)
7f92929e 193{
88cb7938 194// Default Constructor
7f92929e 195}
88cb7938 196/**************************************************************/
7f92929e 197
111e505b 198AliHBTLLWeights::AliHBTLLWeights(const AliHBTLLWeights &/*source*/):
199 TObject(),
200 fTest(kTRUE),
201 fColoumbSwitch(kTRUE),
202 fQuantStatSwitch(kTRUE),
203 fStrongInterSwitch(kTRUE),
204 fColWithResidNuclSwitch(kFALSE),
205 fNuclMass(0.0),
206 fNuclCharge(0.0),
207 fRandomPosition(kFALSE),
208 fRadius(0.0),
209 fPID1(0),
210 fPID2(0),
211 fSigma(0.0)
212{
213 //Copy ctor needed by the coding conventions but not used
214 Fatal("AliHBTLLWeights","copy ctor not implemented");
215}
216/************************************************************/
217
218AliHBTLLWeights& AliHBTLLWeights::operator=(const AliHBTLLWeights& /*source*/)
219{
220 //Assignment operator needed by the coding conventions but not used
221 Fatal("AliHBTLLWeights","assignment operator not implemented");
222 return * this;
223}
224/************************************************************/
225
4fdf4eb3 226AliHBTLLWeights* AliHBTLLWeights::Instance()
88cb7938 227{
228// returns instance of class
229 if (fgLLWeights)
230 {
4fdf4eb3 231 return fgLLWeights;
88cb7938 232 }
233 else
234 {
235 fgLLWeights = new AliHBTLLWeights();
236 return fgLLWeights;
237 }
238}
239
7f92929e 240
241Double_t AliHBTLLWeights::GetWeight(const AliHBTPair* partpair)
242{
88cb7938 243// calculates weight for a pair
244 static const Double_t cmtofm = 1.e13;
111e505b 245 static const Double_t cmtoOneOverGeV = cmtofm*fgkWcons;
88cb7938 246
4fdf4eb3 247 AliHBTParticle *part1 = partpair->Particle1();
248 AliHBTParticle *part2 = partpair->Particle2();
2f8eea63 249
4fdf4eb3 250 if ( (part1 == 0x0) || (part2 == 0x0))
88cb7938 251 {
252 Error("GetWeight","Null particle pointer");
253 return 0.0;
254 }
255
256
257//eats a lot of time
4fdf4eb3 258 if ( (part1->Px() == part2->Px()) &&
259 (part1->Py() == part2->Py()) &&
260 (part1->Pz() == part2->Pz()) )
88cb7938 261 {
262 return 0.0;
263 }
264
4fdf4eb3 265 if ((!fRandomPosition) &&
88cb7938 266 (part1->Vx() == part2->Vx()) &&
267 (part1->Vy() == part2->Vy()) &&
268 (part1->Vz() == part2->Vz()) )
2f8eea63 269 {
4fdf4eb3 270 return 0.0;
2f8eea63 271 }
88cb7938 272
273 FSI_MOM.P1X = part1->Px();
274 FSI_MOM.P1Y = part1->Py();
275 FSI_MOM.P1Z = part1->Pz();
276
277 FSI_MOM.P2X = part2->Px();
278 FSI_MOM.P2Y = part2->Py();
279 FSI_MOM.P2Z = part2->Pz();
280
111e505b 281 FSI_COOR.X1 = part1->Vx()*cmtoOneOverGeV;
282 FSI_COOR.Y1 = part1->Vy()*cmtoOneOverGeV;
283 FSI_COOR.Z1 = part1->Vz()*cmtoOneOverGeV;
88cb7938 284 FSI_COOR.T1 = part1->T();
285
111e505b 286 FSI_COOR.X2 = part2->Vx()*cmtoOneOverGeV;
287 FSI_COOR.Y2 = part2->Vy()*cmtoOneOverGeV;
288 FSI_COOR.Z2 = part2->Vz()*cmtoOneOverGeV;
88cb7938 289 FSI_COOR.T2 = part2->T();
4fdf4eb3 290
291 ltran12();
7f92929e 292
88cb7938 293 //this must be after ltran12 because it would overwrite what we set below
294 if (fRandomPosition)
295 {
296 Double_t rxcm = fSigma*gRandom->Gaus();
297 Double_t rycm = fSigma*gRandom->Gaus();
298 Double_t rzcm = fSigma*gRandom->Gaus();
299
300 FSI_PRF.X=rxcm*fgkWcons;
301 FSI_PRF.Y=rycm*fgkWcons;
302 FSI_PRF.Z=rzcm*fgkWcons;
303 FSI_PRF.T=0.;
304
305 Double_t rps=rxcm*rxcm+rycm*rycm+rzcm*rzcm;
306 Double_t rp=TMath::Sqrt(rps);
307 setpdist(rp);
308 }
309
310 fsiw();
4fdf4eb3 311 return LEDWEIGHT.WEIN;
312}
7f92929e 313/************************************************************/
88cb7938 314
7f92929e 315void AliHBTLLWeights::Init()
4fdf4eb3 316{
88cb7938 317//initial parameters of model
318
4fdf4eb3 319 FSI_NS.NS = fApproximationModel;
320
88cb7938 321 LEDWEIGHT.ITEST = fTest;
322 if(fTest)
323 {
324 FSI_NS.ICH = fColoumbSwitch;
325 FSI_NS.ISI = fStrongInterSwitch;
326 FSI_NS.IQS = fQuantStatSwitch;
327 FSI_NS.I3C = fColWithResidNuclSwitch;
328 LEDWEIGHT.IRANPOS = fRandomPosition;
329 }
330
4fdf4eb3 331 if ( (fPID1 == 0) || (fPID2 == 0) )
88cb7938 332 {
333 Fatal("Init","Particles types are not set");
334 return;//pro forma
335 }
4fdf4eb3 336 FSI_NS.LL = GetPairCode(fPID1,fPID2);
88cb7938 337
4fdf4eb3 338 if (FSI_NS.LL == 0)
88cb7938 339 {
340 Fatal("Init","Particles types are not supported");
341 return;//pro forma
342 }
343
344
4fdf4eb3 345 TParticlePDG* tpart1 = TDatabasePDG::Instance()->GetParticle(fPID1);
346 if (tpart1 == 0x0)
88cb7938 347 {
348 Fatal("init","We can not find particle with ID=%d in PDG DataBase",fPID1);
349 return;
350 }
351
4fdf4eb3 352 FSI_POC.AM1=tpart1->Mass();
353 FSI_POC.C1=tpart1->Charge();
88cb7938 354
4fdf4eb3 355 TParticlePDG* tpart2 = TDatabasePDG::Instance()->GetParticle(fPID2);
88cb7938 356//lv
4fdf4eb3 357 if (tpart2 == 0x0)
88cb7938 358 {
359 Fatal("init","We can not find particle with ID=%d in our DataBase",fPID2);
360 return;
361 }
362
4fdf4eb3 363 FSI_POC.AM2=tpart2->Mass();
364 FSI_POC.C1=tpart2->Charge();
88cb7938 365
4fdf4eb3 366 led_bldata();
367 fsiini();
88cb7938 368
369
370//constants for radii simulation
371
372 if(fRandomPosition)
373 {
374 fSigma =TMath::Sqrt(2.)*fRadius;
375 }
7f92929e 376}
88cb7938 377/************************************************************/
7f92929e 378
379Int_t AliHBTLLWeights::GetPairCode(const AliHBTPair* partpair)
380{
88cb7938 381//returns Code corresponding to that pair
382 return GetPairCode(partpair->Particle1()->GetPdgCode(),partpair->Particle2()->GetPdgCode());
7f92929e 383}
88cb7938 384/************************************************************/
7f92929e 385
386Int_t AliHBTLLWeights::GetPairCode(Int_t pid1,Int_t pid2)
387{
88cb7938 388// returns code corresponding to the pair of PIDs
389// 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
390// 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
391// 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
392// NS=1 y/n: + + + + + - - - - - - - - - - - - - - - - - - - - - - -
393
394//alphas, deuterons and tyts are NOT supported here
395
7f92929e 396 Int_t chargefactor = 1;
397 Int_t hpid; //pid in higher row
398 Int_t lpid; //pid in lower row
399 Int_t code; //pairCode
400
401 Bool_t swap;
402
88cb7938 403//determine the order of selcetion in switch
7f92929e 404 if (TMath::Abs(pid1) < TMath::Abs(pid2) )
88cb7938 405 {
406 if (pid1<0) chargefactor=-1;
407 hpid=pid2*chargefactor;
408 lpid=pid1*chargefactor;
409 swap = kFALSE;
410 }
7f92929e 411 else
88cb7938 412 {
413 if (pid2<0) chargefactor=-1;
414 hpid=pid1*chargefactor;
415 lpid=pid2*chargefactor;
416 swap = kTRUE;
417 }
418
419//mlv
420 hpid=pid1;
421 lpid=pid2;
422
423
424//Determine the pair code
7f92929e 425 switch (hpid) //switch on first particle id
88cb7938 426 {
427 case kNeutron:
7f92929e 428 switch (lpid)
88cb7938 429 {
430 case kNeutron:
431 code = 1; //neutron neutron
432 break;
433
434 case kProton:
435 code = 3; //neutron proton
436 break;
437
438 case kLambda0:
439 code = 28; //neutron lambda
440 break;
441
442 default:
443 return 0; //given pair not supported
444 break;
445 }
7f92929e 446 break;
88cb7938 447
448 case kProton:
7f92929e 449 switch (lpid)
88cb7938 450 {
451 case kProton:
452 code = 2; //proton proton
453 break;
454
455 case kLambda0:
456 code = 27;//proton lambda
457 break;
458
459 default:
460 return 0; //given pair not supported
461 break;
462
463 }
7f92929e 464 break;
88cb7938 465
466 case kPiPlus:
467
7f92929e 468 switch (lpid)
88cb7938 469 {
470 case kPiPlus:
471 code = 7; //piplus piplus
472 break;
473
474 case kPiMinus:
475 code = 5; //piplus piminus
476 break;
477
478 case kKMinus:
479 code = 10; //piplus Kminus
480 break;
481
482 case kKPlus:
483 code = 11; //piplus Kplus
484 break;
485
486 case kProton:
487 code = 12; //piplus proton
488 chargefactor*=-1;
489 break;
490
491 default:
492 return 0; //given pair not supported
493 break;
494 }
7f92929e 495 break;
88cb7938 496 case kPi0:
7f92929e 497 switch (lpid)
88cb7938 498 {
499 case kPi0:
500 code = 6;
501 break;
502
503 default:
504 return 0; //given pair not supported
505 break;
506 }
7f92929e 507 break;
508
88cb7938 509 case kKPlus:
7f92929e 510 switch (lpid)
88cb7938 511 {
512 case kKMinus:
513 code = 14; //Kplus Kminus
514 break;
515
516 case kKPlus:
517 code = 15; //Kplus Kplus
518 break;
519
520 case kProton:
521 code = 16; //Kplus proton
522 break;
523
524 default:
525 return 0; //given pair not supported
526 break;
527 }
7f92929e 528 break;
529
88cb7938 530 case kKMinus:
7f92929e 531 switch (lpid)
88cb7938 532 {
533 case kProton:
534 code = 17; //Kminus proton
535 chargefactor*=1;
536 break;
537
538 default:
539 return 0; //given pair not supported
540 break;
541 }
7f92929e 542 break;
543
88cb7938 544 case kK0:
7f92929e 545 switch (lpid)
88cb7938 546 {
547 case kK0:
548 code = 2; //Kzero Kzero
549 break;
550
551 case kK0Bar:
552 code = 17; //Kzero KzeroBar
553 break;
554
555 default:
556 return 0; //given pair not supported
557 break;
558 }
7f92929e 559 break;
88cb7938 560
561 default: return 0;
562 }
7f92929e 563 return code;
564}
88cb7938 565/************************************************************/
566
567void AliHBTLLWeights::SetTest(Bool_t rtest)
568{
569 //Sets fTest member
570 fTest = rtest;
571}
572/************************************************************/
573
574void AliHBTLLWeights::SetColoumb(Bool_t col)
575{
576 // (ICH in fortran code) Coulomb interaction between the two particles ON (OFF)
577 fColoumbSwitch = col;
578}
579/************************************************************/
580
581void AliHBTLLWeights::SetQuantumStatistics(Bool_t qss)
582{
583 //IQS: quantum statistics for the two particles ON (OFF)
584 //if non-identical particles automatically off
585 fQuantStatSwitch = qss;
586}
587/************************************************************/
588
589void AliHBTLLWeights::SetStrongInterSwitch(Bool_t sis)
590{
591 //ISI: strong interaction between the two particles ON (OFF)
592 fStrongInterSwitch = sis;
593}
594/************************************************************/
7f92929e 595
88cb7938 596void AliHBTLLWeights::SetColWithResidNuclSwitch(Bool_t crn)
597{
598 //I3C: Coulomb interaction with residual nucleus ON (OFF)
599 fColWithResidNuclSwitch = crn;
600}
601/************************************************************/
602
603void AliHBTLLWeights::SetApproxModel(Int_t ap)
604{
605 //sets Model of Approximation (NS in Fortran code)
606 fApproximationModel=ap;
607}
608/************************************************************/
609
610void AliHBTLLWeights::SetRandomPosition(Bool_t rp)
611{
612 //ON=kTRUE(OFF=kFALSE)
613 //ON -- calculation of the Gauss source radii
614 //if the generator don't allows the source generation (for example MeVSim)
615 //if ON the following parameters are requested:
616 fRandomPosition = rp;
617}
618/************************************************************/
619
620void AliHBTLLWeights::SetR1dw(Double_t R)
621{
622 //spherical source model radii
623 fRadius=R;
624}
625/************************************************************/
626
627void AliHBTLLWeights::SetParticlesTypes(Int_t pid1, Int_t pid2)
628{
629 //set AliRoot particles types
630 fPID1 = pid1;
631 fPID2 = pid2;
632}
633/************************************************************/
634
635void AliHBTLLWeights::SetNucleusCharge(Double_t ch)
636{
637 // not used now (see comments in fortran code)
638 fNuclCharge=ch;
639}
640/************************************************************/
641
642void AliHBTLLWeights::SetNucleusMass(Double_t mass)
643{
644 // (see comments in fortran code)
645 fNuclMass=mass;
646}