1 #include "AliHBTLLWeights.h"
2 /**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 //_________________________________________________________________________
18 ///////////////////////////////////////////////////////////////////////////
20 // class AliHBTLLWeights
22 // This class introduces the weight's calculation
23 // according to the Lednicky's algorithm.
28 // Description from fortran code by author R. Lednicky
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
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
49 // NS : approximation used to calculate Bethe-Salpeter amplitude
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,
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
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
103 // AMN = mass of the effective nucleus [GeV/c**2]
104 // CN = charge of the effective nucleus [elem. charge units]
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
119 // and emission times
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
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 //======================================================================
139 /*******************************************************************/
140 /****** ROUTINES USED FOR COMMUNUCATION ********/
141 /******************** WITH FORTRAN ********************/
142 /*******************************************************************/
144 # define led_bldata led_bldata_
145 # define fsiini fsiini_
146 # define ltran12 ltran12_
147 # define boosttoprf boosttoprf_
149 # define setpdist setpdist_
150 # define type_of_call
152 # define led_bldata LED_BLDATA
153 # define fsiini FSIINI
154 # define ltran12 LTRAN12
155 # define boosttoprf BOOSTTOPRF
157 # define setpdist SETPDIST
158 # define type_of_call _stdcall
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 /**************************************************************/
169 #include "AliHBTPair.h"
170 #include "AliVAODParticle.h"
171 #include "WLedCOMMONS.h"
174 #include <TPDGCode.h>
175 #include <TParticlePDG.h>
176 #include <TDatabasePDG.h>
179 ClassImp(AliHBTLLWeights)
181 AliHBTLLWeights* AliHBTLLWeights::fgLLWeights = 0x0;
182 const Double_t AliHBTLLWeights::fgkWcons = 1./0.1973;
184 AliHBTLLWeights::AliHBTLLWeights():
186 fColoumbSwitch(kTRUE),
187 fQuantStatSwitch(kTRUE),
188 fStrongInterSwitch(kTRUE),
189 fColWithResidNuclSwitch(kFALSE),
192 fRandomPosition(kNone),
194 fOneMinusLambda(0.0),
195 fApproximationModel(0),
200 // Default Constructor
202 Fatal("AliHBTLLWeights","LLWeights already instatiated. Use AliHBTLLWeights::Instance()");
204 /**************************************************************/
206 AliHBTLLWeights::AliHBTLLWeights(const AliHBTLLWeights &/*source*/):
209 fColoumbSwitch(kTRUE),
210 fQuantStatSwitch(kTRUE),
211 fStrongInterSwitch(kTRUE),
212 fColWithResidNuclSwitch(kFALSE),
215 fRandomPosition(kNone),
217 fOneMinusLambda(0.0),
218 fApproximationModel(0),
223 //Copy ctor needed by the coding conventions but not used
224 Fatal("AliHBTLLWeights","copy ctor not implemented");
226 /************************************************************/
228 AliHBTLLWeights& AliHBTLLWeights::operator=(const AliHBTLLWeights& /*source*/)
230 //Assignment operator needed by the coding conventions but not used
231 Fatal("AliHBTLLWeights","assignment operator not implemented");
234 /************************************************************/
236 AliHBTLLWeights* AliHBTLLWeights::Instance()
238 // returns instance of class
245 fgLLWeights = new AliHBTLLWeights();
249 /************************************************************/
251 void AliHBTLLWeights::Set()
253 //sets this as weighitng class
254 Info("Set","Setting Lednicky-Lyuboshitz as Weighing Class");
256 if ( fgWeights == 0x0 )
258 fgWeights = AliHBTLLWeights::Instance();
261 if ( fgWeights == AliHBTLLWeights::Instance() ) return;
263 fgWeights = AliHBTLLWeights::Instance();
265 /************************************************************/
267 Double_t AliHBTLLWeights::GetWeight(AliHBTPair* partpair)
269 // calculates weight for a pair
270 static const Double_t kcmtofm = 1.e13;
271 static const Double_t kcmtoOneOverGeV = kcmtofm*fgkWcons;
273 AliVAODParticle *part1 = partpair->Particle1();
274 AliVAODParticle *part2 = partpair->Particle2();
276 if ( (part1 == 0x0) || (part2 == 0x0))
278 Error("GetWeight","Null particle pointer");
282 if ( fPID1 != part1->GetPdgCode() ) return 1.0;
283 if ( fPID2 != part2->GetPdgCode() ) return 1.0;
285 //takes a lot of time
286 if ( (part1->Px() == part2->Px()) &&
287 (part1->Py() == part2->Py()) &&
288 (part1->Pz() == part2->Pz()) )
293 if ((!fRandomPosition) &&
294 (part1->Vx() == part2->Vx()) &&
295 (part1->Vy() == part2->Vy()) &&
296 (part1->Vz() == part2->Vz()) )
301 if(fOneMinusLambda)//implemetation of non-zero intetcept parameter
303 if( gRandom->Rndm() < fOneMinusLambda ) return 1.0;
307 //this must be after ltran12 because it would overwrite what we set below
308 switch (fRandomPosition)
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();
317 FSI_MOM.P2X = part2->Px();
318 FSI_MOM.P2Y = part2->Py();
319 FSI_MOM.P2Z = part2->Pz();
325 //Set particle positions now so they are Gaussian in PRF
326 RandomPairDistances();
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);
335 //boost pair to LCMS and set such a momenta in the common block
336 Int_t retv = SetMomentaInLCMS(part1,part2);
338 //random particle positions/distance so they are Gaussian in LCMS
339 RandomPairDistances();
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();
354 FSI_MOM.P2X = part2->Px();
355 FSI_MOM.P2Y = part2->Py();
356 FSI_MOM.P2Z = part2->Pz();
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;
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;
375 return LEDWEIGHT.WEIN;
377 /************************************************************/
379 void AliHBTLLWeights::Init()
381 //initial parameters of model
383 FSI_NS.NS = fApproximationModel;
385 LEDWEIGHT.ITEST = fTest;
388 FSI_NS.ICH = fColoumbSwitch;
389 FSI_NS.ISI = fStrongInterSwitch;
390 FSI_NS.IQS = fQuantStatSwitch;
391 FSI_NS.I3C = fColWithResidNuclSwitch;
392 LEDWEIGHT.IRANPOS = fRandomPosition;
395 if ( (fPID1 == 0) || (fPID2 == 0) )
397 Fatal("Init","Particles types are not set");
402 FSI_NS.LL = GetPairCode(fPID1,fPID2);
406 Fatal("Init","Particles types are not supported");
410 Info("Init","Setting PIDs %d %d. LL Code is %d",fPID1,fPID2,FSI_NS.LL);
413 TParticlePDG* tpart1 = TDatabasePDG::Instance()->GetParticle(fPID1);
416 Fatal("init","We can not find particle with ID=%d in PDG DataBase",fPID1);
420 FSI_POC.AM1=tpart1->Mass();
421 FSI_POC.C1=tpart1->Charge();
423 TParticlePDG* tpart2 = TDatabasePDG::Instance()->GetParticle(fPID2);
427 Fatal("init","We can not find particle with ID=%d in our DataBase",fPID2);
431 FSI_POC.AM2=tpart2->Mass();
432 FSI_POC.C1=tpart2->Charge();
438 //constants for radii simulation
442 fSigma =TMath::Sqrt(2.)*fRadius;
445 /************************************************************/
447 Int_t AliHBTLLWeights::GetPairCode(const AliHBTPair* partpair)
449 //returns Code corresponding to that pair
450 return GetPairCode(partpair->Particle1()->GetPdgCode(),partpair->Particle2()->GetPdgCode());
452 /************************************************************/
454 Int_t AliHBTLLWeights::GetPairCode(Int_t pid1,Int_t pid2)
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: + + + + + - - - - - - - - - - - - - - - - - - - - - - -
462 //alphas, deuterons and tyts are NOT supported here
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
471 //determine the order of selcetion in switch
472 if (TMath::Abs(pid1) < TMath::Abs(pid2) )
474 if (pid1<0) chargefactor=-1;
475 hpid=pid2*chargefactor;
476 lpid=pid1*chargefactor;
481 if (pid2<0) chargefactor=-1;
482 hpid=pid1*chargefactor;
483 lpid=pid2*chargefactor;
492 //Determine the pair code
493 switch (hpid) //switch on first particle id
499 code = 1; //neutron neutron
503 code = 3; //neutron proton
507 code = 28; //neutron lambda
511 return 0; //given pair not supported
520 code = 2; //proton proton
524 code = 27;//proton lambda
528 return 0; //given pair not supported
539 code = 7; //piplus piplus
543 code = 5; //piplus piminus
547 code = 10; //piplus Kminus
551 code = 11; //piplus Kplus
555 code = 12; //piplus proton
560 return 0; //given pair not supported
572 return 0; //given pair not supported
581 code = 14; //Kplus Kminus
585 code = 15; //Kplus Kplus
589 code = 16; //Kplus proton
593 return 0; //given pair not supported
602 code = 17; //Kminus proton
607 return 0; //given pair not supported
616 code = 2; //Kzero Kzero
620 code = 17; //Kzero KzeroBar
624 return 0; //given pair not supported
633 /************************************************************/
635 void AliHBTLLWeights::SetTest(Bool_t rtest)
640 /************************************************************/
642 void AliHBTLLWeights::SetColoumb(Bool_t col)
644 // (ICH in fortran code) Coulomb interaction between the two particles ON (OFF)
645 fColoumbSwitch = col;
647 /************************************************************/
649 void AliHBTLLWeights::SetQuantumStatistics(Bool_t qss)
651 //IQS: quantum statistics for the two particles ON (OFF)
652 //if non-identical particles automatically off
653 fQuantStatSwitch = qss;
655 /************************************************************/
657 void AliHBTLLWeights::SetStrongInterSwitch(Bool_t sis)
659 //ISI: strong interaction between the two particles ON (OFF)
660 fStrongInterSwitch = sis;
662 /************************************************************/
664 void AliHBTLLWeights::SetColWithResidNuclSwitch(Bool_t crn)
666 //I3C: Coulomb interaction with residual nucleus ON (OFF)
667 fColWithResidNuclSwitch = crn;
669 /************************************************************/
671 void AliHBTLLWeights::SetApproxModel(Int_t ap)
673 //sets Model of Approximation (NS in Fortran code)
674 fApproximationModel=ap;
676 /************************************************************/
678 void AliHBTLLWeights::SetRandomPosition(ERandomizationWay rw)
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;
686 /************************************************************/
688 void AliHBTLLWeights::SetR1dw(Double_t R)
690 //spherical source model radii
693 /************************************************************/
695 void AliHBTLLWeights::SetParticlesTypes(Int_t pid1, Int_t pid2)
697 //set AliRoot particles types
701 /************************************************************/
703 void AliHBTLLWeights::SetNucleusCharge(Double_t ch)
705 // not used now (see comments in fortran code)
708 /************************************************************/
710 void AliHBTLLWeights::SetNucleusMass(Double_t mass)
712 // (see comments in fortran code)
715 /************************************************************/
717 Int_t AliHBTLLWeights::SetMomentaInLCMS(AliVAODParticle* part1, AliVAODParticle* part2)
719 //sets paricle momenta in the common block in LCMS frame
720 //---> Particle energies ---------
722 Double_t p1x = part1->Px();
723 Double_t p1y = part1->Py();
724 Double_t p1z = part1->Pz();
726 Double_t p2x = part2->Px();
727 Double_t p2y = part2->Py();
728 Double_t p2z = part2->Pz();
730 Double_t am1 = part1->Mass();
731 Double_t am2 = part2->Mass();
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
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)
753 Double_t p12ts=p12x*p12x+p12y*p12y;
754 Double_t p12t=TMath::Sqrt(p12ts); //pt
755 //===> azimuthal rotation (pt||x) ============
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);
765 else // rotation impossible
770 //===> co-moving ref. frame ============
773 Boost(p1z,e1,v12z,gamz,p1z,nothing);
774 Boost(p2z,e2,v12z,gamz,p2z,nothing);
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);
792 // Info("Simulate"," %f %f %f ",p1x + p2x, p1y + p2y, p1z + p2z);
796 /**********************************************************/
798 void AliHBTLLWeights::RandomPairDistances()
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();
811 FSI_COOR.X2 = rxcm*fgkWcons;
812 FSI_COOR.Y2 = rycm*fgkWcons;
813 FSI_COOR.Z2 = rzcm*fgkWcons;
819 void AliHBTLLWeights::Rotate(Double_t x, Double_t y, Double_t sphi, Double_t cphi, Double_t& xr, Double_t& yr)
826 void AliHBTLLWeights::Boost(Double_t z, Double_t t, Double_t beta, Double_t gamma, Double_t& zt, Double_t& yt)