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
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_
148 # define setpdist setpdist_
149 # define type_of_call
151 # define led_bldata LED_BLDATA
152 # define fsiini FSIINI
153 # define ltran12 LTRAN12
155 # define setpdist SETPDIST
156 # define type_of_call _stdcall
158 /****************************************************************/
159 extern "C" void type_of_call led_bldata();
160 extern "C" void type_of_call fsiini();
161 extern "C" void type_of_call ltran12();
162 extern "C" void type_of_call fsiw();
163 extern "C" void type_of_call setpdist(Double_t& r);
164 /**************************************************************/
166 #include "AliHBTPair.h"
167 #include "AliVAODParticle.h"
168 #include "WLedCOMMONS.h"
171 #include <TPDGCode.h>
172 #include <TParticlePDG.h>
173 #include <TDatabasePDG.h>
176 ClassImp(AliHBTLLWeights)
178 AliHBTLLWeights* AliHBTLLWeights::fgLLWeights = 0x0;
179 const Double_t AliHBTLLWeights::fgkWcons = 1./0.1973;
181 AliHBTLLWeights::AliHBTLLWeights():
183 fColoumbSwitch(kTRUE),
184 fQuantStatSwitch(kTRUE),
185 fStrongInterSwitch(kTRUE),
186 fColWithResidNuclSwitch(kFALSE),
189 fRandomPosition(kFALSE),
191 fOneMinusLambda(0.0),
196 // Default Constructor
198 Fatal("AliHBTLLWeights","LLWeights already instatiated. Use AliHBTLLWeights::Instance()");
200 /**************************************************************/
202 AliHBTLLWeights::AliHBTLLWeights(const AliHBTLLWeights &/*source*/):
205 fColoumbSwitch(kTRUE),
206 fQuantStatSwitch(kTRUE),
207 fStrongInterSwitch(kTRUE),
208 fColWithResidNuclSwitch(kFALSE),
211 fRandomPosition(kFALSE),
213 fOneMinusLambda(0.0),
218 //Copy ctor needed by the coding conventions but not used
219 Fatal("AliHBTLLWeights","copy ctor not implemented");
221 /************************************************************/
223 AliHBTLLWeights& AliHBTLLWeights::operator=(const AliHBTLLWeights& /*source*/)
225 //Assignment operator needed by the coding conventions but not used
226 Fatal("AliHBTLLWeights","assignment operator not implemented");
229 /************************************************************/
231 AliHBTLLWeights* AliHBTLLWeights::Instance()
233 // returns instance of class
240 fgLLWeights = new AliHBTLLWeights();
244 /************************************************************/
246 void AliHBTLLWeights::Set()
248 //sets this as weighitng class
249 Info("Set","Setting Lednicky-Lyuboshitz as Weighing Class");
251 if ( fgWeights == 0x0 )
253 fgWeights = AliHBTLLWeights::Instance();
256 if ( fgWeights == AliHBTLLWeights::Instance() ) return;
258 fgWeights = AliHBTLLWeights::Instance();
260 /************************************************************/
262 Double_t AliHBTLLWeights::GetWeight(const AliHBTPair* partpair)
264 // calculates weight for a pair
265 static const Double_t kcmtofm = 1.e13;
266 static const Double_t kcmtoOneOverGeV = kcmtofm*fgkWcons;
268 AliVAODParticle *part1 = partpair->Particle1();
269 AliVAODParticle *part2 = partpair->Particle2();
271 if ( (part1 == 0x0) || (part2 == 0x0))
273 Error("GetWeight","Null particle pointer");
277 if ( fPID1 != part1->GetPdgCode() ) return 1.0;
278 if ( fPID2 != part2->GetPdgCode() ) return 1.0;
280 //takes a lot of time
281 if ( (part1->Px() == part2->Px()) &&
282 (part1->Py() == part2->Py()) &&
283 (part1->Pz() == part2->Pz()) )
288 if ((!fRandomPosition) &&
289 (part1->Vx() == part2->Vx()) &&
290 (part1->Vy() == part2->Vy()) &&
291 (part1->Vz() == part2->Vz()) )
296 if(fOneMinusLambda)//implemetation of non-zero intetcept parameter
298 if( gRandom->Rndm() < fOneMinusLambda ) return 1.0;
301 FSI_MOM.P1X = part1->Px();
302 FSI_MOM.P1Y = part1->Py();
303 FSI_MOM.P1Z = part1->Pz();
305 FSI_MOM.P2X = part2->Px();
306 FSI_MOM.P2Y = part2->Py();
307 FSI_MOM.P2Z = part2->Pz();
309 FSI_COOR.X1 = part1->Vx()*kcmtoOneOverGeV;
310 FSI_COOR.Y1 = part1->Vy()*kcmtoOneOverGeV;
311 FSI_COOR.Z1 = part1->Vz()*kcmtoOneOverGeV;
312 FSI_COOR.T1 = part1->T();
314 FSI_COOR.X2 = part2->Vx()*kcmtoOneOverGeV;
315 FSI_COOR.Y2 = part2->Vy()*kcmtoOneOverGeV;
316 FSI_COOR.Z2 = part2->Vz()*kcmtoOneOverGeV;
317 FSI_COOR.T2 = part2->T();
321 //this must be after ltran12 because it would overwrite what we set below
324 Double_t rxcm = fSigma*gRandom->Gaus();
325 Double_t rycm = fSigma*gRandom->Gaus();
326 Double_t rzcm = fSigma*gRandom->Gaus();
328 FSI_PRF.X=rxcm*fgkWcons;
329 FSI_PRF.Y=rycm*fgkWcons;
330 FSI_PRF.Z=rzcm*fgkWcons;
333 Double_t rps=rxcm*rxcm+rycm*rycm+rzcm*rzcm;
334 Double_t rp=TMath::Sqrt(rps);
339 return LEDWEIGHT.WEIN;
341 /************************************************************/
343 void AliHBTLLWeights::Init()
345 //initial parameters of model
347 FSI_NS.NS = fApproximationModel;
349 LEDWEIGHT.ITEST = fTest;
352 FSI_NS.ICH = fColoumbSwitch;
353 FSI_NS.ISI = fStrongInterSwitch;
354 FSI_NS.IQS = fQuantStatSwitch;
355 FSI_NS.I3C = fColWithResidNuclSwitch;
356 LEDWEIGHT.IRANPOS = fRandomPosition;
359 if ( (fPID1 == 0) || (fPID2 == 0) )
361 Fatal("Init","Particles types are not set");
366 FSI_NS.LL = GetPairCode(fPID1,fPID2);
370 Fatal("Init","Particles types are not supported");
374 Info("Init","Setting PIDs %d %d. LL Code is %d",fPID1,fPID2,FSI_NS.LL);
377 TParticlePDG* tpart1 = TDatabasePDG::Instance()->GetParticle(fPID1);
380 Fatal("init","We can not find particle with ID=%d in PDG DataBase",fPID1);
384 FSI_POC.AM1=tpart1->Mass();
385 FSI_POC.C1=tpart1->Charge();
387 TParticlePDG* tpart2 = TDatabasePDG::Instance()->GetParticle(fPID2);
391 Fatal("init","We can not find particle with ID=%d in our DataBase",fPID2);
395 FSI_POC.AM2=tpart2->Mass();
396 FSI_POC.C1=tpart2->Charge();
402 //constants for radii simulation
406 fSigma =TMath::Sqrt(2.)*fRadius;
409 /************************************************************/
411 Int_t AliHBTLLWeights::GetPairCode(const AliHBTPair* partpair)
413 //returns Code corresponding to that pair
414 return GetPairCode(partpair->Particle1()->GetPdgCode(),partpair->Particle2()->GetPdgCode());
416 /************************************************************/
418 Int_t AliHBTLLWeights::GetPairCode(Int_t pid1,Int_t pid2)
420 // returns code corresponding to the pair of PIDs
421 // 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
422 // 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
423 // 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
424 // NS=1 y/n: + + + + + - - - - - - - - - - - - - - - - - - - - - - -
426 //alphas, deuterons and tyts are NOT supported here
428 Int_t chargefactor = 1;
429 Int_t hpid; //pid in higher row
430 Int_t lpid; //pid in lower row
431 Int_t code; //pairCode
435 //determine the order of selcetion in switch
436 if (TMath::Abs(pid1) < TMath::Abs(pid2) )
438 if (pid1<0) chargefactor=-1;
439 hpid=pid2*chargefactor;
440 lpid=pid1*chargefactor;
445 if (pid2<0) chargefactor=-1;
446 hpid=pid1*chargefactor;
447 lpid=pid2*chargefactor;
456 //Determine the pair code
457 switch (hpid) //switch on first particle id
463 code = 1; //neutron neutron
467 code = 3; //neutron proton
471 code = 28; //neutron lambda
475 return 0; //given pair not supported
484 code = 2; //proton proton
488 code = 27;//proton lambda
492 return 0; //given pair not supported
503 code = 7; //piplus piplus
507 code = 5; //piplus piminus
511 code = 10; //piplus Kminus
515 code = 11; //piplus Kplus
519 code = 12; //piplus proton
524 return 0; //given pair not supported
536 return 0; //given pair not supported
545 code = 14; //Kplus Kminus
549 code = 15; //Kplus Kplus
553 code = 16; //Kplus proton
557 return 0; //given pair not supported
566 code = 17; //Kminus proton
571 return 0; //given pair not supported
580 code = 2; //Kzero Kzero
584 code = 17; //Kzero KzeroBar
588 return 0; //given pair not supported
597 /************************************************************/
599 void AliHBTLLWeights::SetTest(Bool_t rtest)
604 /************************************************************/
606 void AliHBTLLWeights::SetColoumb(Bool_t col)
608 // (ICH in fortran code) Coulomb interaction between the two particles ON (OFF)
609 fColoumbSwitch = col;
611 /************************************************************/
613 void AliHBTLLWeights::SetQuantumStatistics(Bool_t qss)
615 //IQS: quantum statistics for the two particles ON (OFF)
616 //if non-identical particles automatically off
617 fQuantStatSwitch = qss;
619 /************************************************************/
621 void AliHBTLLWeights::SetStrongInterSwitch(Bool_t sis)
623 //ISI: strong interaction between the two particles ON (OFF)
624 fStrongInterSwitch = sis;
626 /************************************************************/
628 void AliHBTLLWeights::SetColWithResidNuclSwitch(Bool_t crn)
630 //I3C: Coulomb interaction with residual nucleus ON (OFF)
631 fColWithResidNuclSwitch = crn;
633 /************************************************************/
635 void AliHBTLLWeights::SetApproxModel(Int_t ap)
637 //sets Model of Approximation (NS in Fortran code)
638 fApproximationModel=ap;
640 /************************************************************/
642 void AliHBTLLWeights::SetRandomPosition(Bool_t rp)
644 //ON=kTRUE(OFF=kFALSE)
645 //ON -- calculation of the Gauss source radii
646 //if the generator don't allows the source generation (for example MeVSim)
647 //if ON the following parameters are requested:
648 fRandomPosition = rp;
650 /************************************************************/
652 void AliHBTLLWeights::SetR1dw(Double_t R)
654 //spherical source model radii
657 /************************************************************/
659 void AliHBTLLWeights::SetParticlesTypes(Int_t pid1, Int_t pid2)
661 //set AliRoot particles types
665 /************************************************************/
667 void AliHBTLLWeights::SetNucleusCharge(Double_t ch)
669 // not used now (see comments in fortran code)
672 /************************************************************/
674 void AliHBTLLWeights::SetNucleusMass(Double_t mass)
676 // (see comments in fortran code)