Retrofeed from the Release developement
authorskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 9 Feb 2005 17:08:43 +0000 (17:08 +0000)
committerskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 9 Feb 2005 17:08:43 +0000 (17:08 +0000)
HBTAN/AliHBTCrab.cxx
HBTAN/AliHBTCrab.h
HBTAN/AliHBTLLWeights.cxx
HBTAN/AliHBTLLWeights.h
HBTAN/AliHBTPositionRandomizer.cxx
HBTAN/AliHBTWeights.h
HBTAN/ltran12.F

index 242b34e..d0623c4 100644 (file)
@@ -113,7 +113,7 @@ Bool_t AliHBTCrab::SetConfig(const AliHBTPair* pair)
 }
 //===================================================================
 
-Double_t AliHBTCrab::GetWeight(const AliHBTPair* partpair)
+Double_t AliHBTCrab::GetWeight(AliHBTPair* partpair)
 {
 //returns the weight
   Double_t qred, r, qdotr, mom;
index 93857fe..b6fe448 100644 (file)
@@ -39,7 +39,7 @@ class AliHBTCrab: public AliHBTWeights
      static AliHBTCrab* Instance();
      void Set();
 
-     Double_t GetWeight(const AliHBTPair* partpair);
+     Double_t GetWeight(AliHBTPair* partpair);
      void Init(Int_t pid1,Int_t pid2); //put the initial values in fortran commons fsiini, led_bldata
      
    private:
index d7a825b..f89a859 100644 (file)
@@ -90,7 +90,7 @@
 //   W=1/.1973D0    ! from fm to 1/GeV
 //   PI=4*DATAN(1.D0)
 //   PI2=2*PI
-//   SPI=DSQRT(PI)
+//   SPI=TMath::Sqrt(PI)
 //   DR=180.D0/PI   ! from radian to degree
 //    _______________________________________________________
 //  !! |Important note: all real quantities are assumed REAL*8 | !!
 # define led_bldata led_bldata_
 # define fsiini fsiini_
 # define ltran12 ltran12_
+# define boosttoprf boosttoprf_
 # define fsiw fsiw_
 # define setpdist setpdist_
 # define type_of_call
 # define led_bldata LED_BLDATA
 # define fsiini FSIINI
 # define ltran12 LTRAN12
+# define boosttoprf BOOSTTOPRF
 # define fsiw FSIW
 # define setpdist SETPDIST
 # define type_of_call _stdcall
 extern "C" void type_of_call led_bldata(); 
 extern "C" void type_of_call fsiini();
 extern "C" void type_of_call ltran12();
+extern "C" void type_of_call boosttoprf();
 extern "C" void type_of_call fsiw();
 extern "C" void type_of_call setpdist(Double_t& r);
 /**************************************************************/
@@ -186,7 +189,7 @@ AliHBTLLWeights::AliHBTLLWeights():
  fColWithResidNuclSwitch(kFALSE),
  fNuclMass(0.0),
  fNuclCharge(0.0),
- fRandomPosition(kFALSE),
+ fRandomPosition(kNone),
  fRadius(0.0),
  fOneMinusLambda(0.0),
  fPID1(0),
@@ -208,7 +211,7 @@ AliHBTLLWeights::AliHBTLLWeights(const AliHBTLLWeights &/*source*/):
  fColWithResidNuclSwitch(kFALSE),
  fNuclMass(0.0),
  fNuclCharge(0.0),
- fRandomPosition(kFALSE),
+ fRandomPosition(kNone),
  fRadius(0.0),
  fOneMinusLambda(0.0),
  fPID1(0),
@@ -259,7 +262,7 @@ void AliHBTLLWeights::Set()
 }
 /************************************************************/
 
-Double_t AliHBTLLWeights::GetWeight(const AliHBTPair* partpair)
+Double_t AliHBTLLWeights::GetWeight(AliHBTPair* partpair)
 {
 // calculates weight for a pair
   static const Double_t kcmtofm = 1.e13;
@@ -298,43 +301,74 @@ Double_t AliHBTLLWeights::GetWeight(const AliHBTPair* partpair)
      if( gRandom->Rndm() < fOneMinusLambda ) return 1.0;
    }
     
-  FSI_MOM.P1X = part1->Px();
-  FSI_MOM.P1Y = part1->Py();
-  FSI_MOM.P1Z = part1->Pz();
-      
-  FSI_MOM.P2X = part2->Px();
-  FSI_MOM.P2Y = part2->Py();
-  FSI_MOM.P2Z = part2->Pz();
-
-  FSI_COOR.X1 = part1->Vx()*kcmtoOneOverGeV;
-  FSI_COOR.Y1 = part1->Vy()*kcmtoOneOverGeV;
-  FSI_COOR.Z1 = part1->Vz()*kcmtoOneOverGeV;
-  FSI_COOR.T1 = part1->T();
-
-  FSI_COOR.X2 = part2->Vx()*kcmtoOneOverGeV;
-  FSI_COOR.Y2 = part2->Vy()*kcmtoOneOverGeV;
-  FSI_COOR.Z2 = part2->Vz()*kcmtoOneOverGeV;
-  FSI_COOR.T2 = part2->T();
-  
-  ltran12();
 
   //this must  be after ltran12 because it would overwrite what we set below
-  if (fRandomPosition)
+  switch (fRandomPosition)
    {
-     Double_t rxcm = fSigma*gRandom->Gaus();
-     Double_t rycm = fSigma*gRandom->Gaus();
-     Double_t rzcm = fSigma*gRandom->Gaus();
-
-     FSI_PRF.X=rxcm*fgkWcons;
-     FSI_PRF.Y=rycm*fgkWcons;
-     FSI_PRF.Z=rzcm*fgkWcons;
-     FSI_PRF.T=0.;
-
-     Double_t rps=rxcm*rxcm+rycm*rycm+rzcm*rzcm;
-     Double_t rp=TMath::Sqrt(rps);
-     setpdist(rp);
+    case kGaussianQInv:
+     {
+       // Set Momenta in the common block
+       FSI_MOM.P1X = part1->Px();
+       FSI_MOM.P1Y = part1->Py();
+       FSI_MOM.P1Z = part1->Pz();
+
+       FSI_MOM.P2X = part2->Px();
+       FSI_MOM.P2Y = part2->Py();
+       FSI_MOM.P2Z = part2->Pz();
+      
+       //boost it to PRF
+       ltran12();
+       boosttoprf();
+       
+       //Set particle positions now so they are Gaussian in PRF
+       RandomPairDistances();
+       
+       FSI_PRF.RPS= FSI_COOR.X2*FSI_COOR.X2 + FSI_COOR.Y2*FSI_COOR.Y2 +FSI_COOR.Z2*FSI_COOR.Z2;
+       FSI_PRF.RP=TMath::Sqrt(FSI_PRF.RPS);
+       
+       break;
+      } 
+    case kGaussianOSL:
+      { 
+       //boost pair to LCMS and set such a momenta in the common block
+       Int_t retv = SetMomentaInLCMS(part1,part2);
+       if (retv) return 1;
+       //random particle positions/distance so they are Gaussian in LCMS
+       RandomPairDistances();
+       
+       //Boost to PRF
+       ltran12();
+       boosttoprf();
+    
+       break;
+      }
+    case kNone:
+    default:
+       //set momenta and paricle positions as they are
+       FSI_MOM.P1X = part1->Px();
+       FSI_MOM.P1Y = part1->Py();
+       FSI_MOM.P1Z = part1->Pz();
+
+       FSI_MOM.P2X = part2->Px();
+       FSI_MOM.P2Y = part2->Py();
+       FSI_MOM.P2Z = part2->Pz();
+
+       FSI_COOR.X1 = part1->Vx()*kcmtoOneOverGeV;
+       FSI_COOR.Y1 = part1->Vy()*kcmtoOneOverGeV;
+       FSI_COOR.Z1 = part1->Vz()*kcmtoOneOverGeV;
+       FSI_COOR.T1 = part1->T()*kcmtoOneOverGeV;
+
+       FSI_COOR.X2 = part2->Vx()*kcmtoOneOverGeV;
+       FSI_COOR.Y2 = part2->Vy()*kcmtoOneOverGeV;
+       FSI_COOR.Z2 = part2->Vz()*kcmtoOneOverGeV;
+       FSI_COOR.T2 = part2->T()*kcmtoOneOverGeV;
+
+       ltran12();
+       boosttoprf();
+       
+       break;
    }
-               
+
   fsiw();
   return LEDWEIGHT.WEIN;
 }
@@ -639,13 +673,13 @@ void AliHBTLLWeights::SetApproxModel(Int_t ap)
 }
 /************************************************************/
      
-void AliHBTLLWeights::SetRandomPosition(Bool_t rp)
+void AliHBTLLWeights::SetRandomPosition(ERandomizationWay rw)
 { 
- //ON=kTRUE(OFF=kFALSE)
- //ON -- calculation of the Gauss source radii 
- //if the generator don't allows the source generation (for example MeVSim)
- //if ON the following parameters are requested:
- fRandomPosition = rp;
+// rw can be: kGaussianQInv - so 1D Qinv correlation function has a Gaussian shape 
+//                            (and also Qlong and Qside, but not Qout)
+//            kGaussianOSL - so 3D Qout-Qside-Qlong correlation function has a Gaussian shape
+//            kNone  - no randomization performed (DEFAULT)
+ fRandomPosition = rw;
 }
 /************************************************************/
 
@@ -676,3 +710,120 @@ void AliHBTLLWeights::SetNucleusMass(Double_t mass)
   // (see comments in fortran code)
   fNuclMass=mass;
 }
+/************************************************************/
+
+Int_t AliHBTLLWeights::SetMomentaInLCMS(AliVAODParticle* part1, AliVAODParticle* part2)
+{
+//sets paricle momenta in the common block in LCMS frame
+//---> Particle energies ---------
+
+  Double_t p1x = part1->Px();
+  Double_t p1y = part1->Py();
+  Double_t p1z = part1->Pz();
+      
+  Double_t p2x = part2->Px();
+  Double_t p2y = part2->Py();
+  Double_t p2z = part2->Pz();
+
+  Double_t am1 = part1->Mass();
+  Double_t am2 = part2->Mass();
+  
+  Double_t p1s=p1x*p1x+p1y*p1y+p1z*p1z;
+  Double_t p2s=p2x*p2x+p2y*p2y+p2z*p2z;
+  Double_t e1=TMath::Sqrt(am1*am1+p1s);
+  Double_t e2=TMath::Sqrt(am2*am2+p2s);
+//---> pair parameters -----------
+  Double_t e12=e1+e2;       // energy
+  Double_t p12x=p1x+p2x;    // px
+  Double_t p12y=p1y+p2y;    // py
+  Double_t p12z=p1z+p2z;    // pz
+  Double_t p12s=p12x*p12x+p12y*p12y+p12z*p12z;
+  Double_t p12 =TMath::Sqrt(p12s);// momentum
+  Double_t v12 =p12/e12;    // velocity
+
+  Double_t cth =p12z/p12;   // cos(theta)
+//  Double_t sth =TMath::Sqrt(1. - cth*cth); //sin
+  Double_t v12z=v12*cth;    // longit. v
+  Double_t gamz=1./TMath::Sqrt(1.-v12z*v12z);
+//--      v12t=v12*sth    // transv. v in cms (not needed)
+
+
+  Double_t p12ts=p12x*p12x+p12y*p12y;
+  Double_t p12t=TMath::Sqrt(p12ts); //pt
+//===> azimuthal rotation (pt||x) ============
+   if(p12t != 0.0) 
+    {
+     Double_t cphi=p12x/p12t;  // cos(phi)
+     Double_t sphi=p12y/p12t;  // sin(phi)
+     Rotate(p1x,p1y,sphi,cphi,p1x,p1y);       
+     Rotate(p2x,p2y,sphi,cphi,p2x,p2y);
+//     Rotate(x1,y1,sphi,cphi,x1,y1);
+//     Rotate(x2,y2,sphi,cphi,x2,y2);
+     }    
+   else // rotation impossible 
+    {
+      return 1;
+    }
+
+//===> co-moving ref. frame       ============
+   
+   Double_t nothing;
+   Boost(p1z,e1,v12z,gamz,p1z,nothing);
+   Boost(p2z,e2,v12z,gamz,p2z,nothing);
+
+//   p1s=p1x*p1x+p1y*p1y+p1z*p1z;
+//   p2s=p2x*p2x+p2y*p2y+p2z*p2z;
+//   e1=TMath::Sqrt(am1*am1+p1s);
+//   e2=TMath::Sqrt(am2*am2+p2s);
+//   Boost(z1,t1,v12z,gamz,z1,t1);
+//   Boost(z2,t2,v12z,gamz,z2,t2);
+
+
+     FSI_MOM.P1X = p1x;
+     FSI_MOM.P1Y = p1y;
+     FSI_MOM.P1Z = p1z;
+
+     FSI_MOM.P2X = p2x;
+     FSI_MOM.P2Y = p2y;
+     FSI_MOM.P2Z = p2z;
+     
+//     Info("Simulate"," %f %f %f ",p1x + p2x, p1y + p2y, p1z + p2z);
+     return 0;
+}
+
+/**********************************************************/  
+
+void AliHBTLLWeights::RandomPairDistances()
+{
+ //randomizes pair distances
+  Double_t rxcm = fSigma*gRandom->Gaus();
+  Double_t rycm = fSigma*gRandom->Gaus();
+  Double_t rzcm = fSigma*gRandom->Gaus();
+
+
+  FSI_COOR.X1 = 0;
+  FSI_COOR.Y1 = 0;
+  FSI_COOR.Z1 = 0;
+  FSI_COOR.T1 = 0;
+
+  FSI_COOR.X2 = rxcm*fgkWcons;
+  FSI_COOR.Y2 = rycm*fgkWcons;
+  FSI_COOR.Z2 = rzcm*fgkWcons;
+  FSI_COOR.T2 = 0;
+
+}
+
+
+void AliHBTLLWeights::Rotate(Double_t x, Double_t y, Double_t sphi, Double_t cphi, Double_t& xr, Double_t& yr)
+{
+//rotates 
+ yr=y*cphi-x*sphi;
+ xr=x*cphi+y*sphi;
+}
+
+void AliHBTLLWeights::Boost(Double_t z, Double_t t, Double_t beta, Double_t gamma, Double_t& zt, Double_t& yt)
+{
+//boosts
+  zt=gamma*(z-beta*t);
+  yt=gamma*(t-beta*z);
+}
index fe1ff8c..47561f7 100644 (file)
@@ -1,10 +1,28 @@
 /* $Id$ */
 
+//_________________________________________________________________
+///////////////////////////////////////////////////////////////////////
+//
+// class AliHBTLLWeights
+//
 // This class introduces the weight's calculation 
 // according to the Lednicky's algorithm.
 // The detailed description of the algorithm can be found 
 // in comments to fortran code:
 // fsiw.f, fsiini.f  
+//
+// Source simulation: Particle position/distance randomization
+// This class has functionality of randomizing particle positions / pair distances
+// so obtained correlation function follows a given shape, 
+// depending on a value of fRandomPosition member
+// By default this feature is switched off: fRandomPosition == AliHBTLLWeights::kNone
+// It can be: kGaussianQInv - so 1D Qinv correlation function has a Gaussian shape 
+//                            (and also Qlong and Qside, but not Qout)
+//            kGaussianOSL - so 3D Qout-Qside-Qlong correlation function has a Gaussian shape
+//
+// Ludmila Malinina, Piotr Krzysztof Skowronski
+//
+///////////////////////////////////////////////////////////////////////
 
 #ifndef ALIHBTLLWEIGHTS_H
 #define ALIHBTLLWEIGHTS_H
@@ -12,6 +30,7 @@
 #include "AliHBTWeights.h"
 
 class AliHBTPair;
+class AliVAODParticle;
 class AliHBTLLWeights: public AliHBTWeights
  {
    public:
@@ -21,7 +40,7 @@ class AliHBTLLWeights: public AliHBTWeights
      
      void Set();
      
-     Double_t GetWeight(const AliHBTPair* partpair); //get weight calculated by Lednicky's algorithm
+     Double_t GetWeight(AliHBTPair* partpair); //get weight calculated by Lednicky's algorithm
 
      void Init(); //put the initial values in fortran commons fsiini, led_bldata
      void SetTest(Bool_t rtest = kTRUE);//Sets fTest member
@@ -32,12 +51,15 @@ class AliHBTLLWeights: public AliHBTWeights
      void SetColWithResidNuclSwitch(Bool_t crn = kTRUE);//I3C: Coulomb interaction with residual nucleus ON (OFF)  
      void SetLambda(Double_t la){fOneMinusLambda=1.-la;}  //lambda=haoticity
      void SetApproxModel(Int_t ap);//sets  Model of Approximation (NS in Fortran code)
-     void SetRandomPosition(Bool_t rp = kTRUE); //ON=kTRUE(OFF=kFALSE)
-     void SetR1dw(Double_t R);   //spherical source model radii                                                                                                                                                                                                
      void SetParticlesTypes(Int_t pid1, Int_t pid2); //set AliRoot particles types   
      void SetNucleusCharge(Double_t ch); // not used now  (see comments in fortran code)
      void SetNucleusMass(Double_t mass); // (see comments in fortran code)
      
+     enum ERandomizationWay {kNone = 0, kGaussianQInv, kGaussianOSL};
+
+     void SetRandomPosition(ERandomizationWay rw = kNone); //Choose if, and if yes, how to randomize distance between particles
+     void SetR1dw(Double_t R);   //spherical source model radii                                                                                                                                                                                                
+     
    protected:
      
      Bool_t fTest;           //flag indicating if parameters listed below are to be calculated automatically (0)
@@ -50,7 +72,7 @@ class AliHBTLLWeights: public AliHBTWeights
      Double_t fNuclMass;   //mass of nucleus
      Double_t fNuclCharge; //charge of nucleus
      
-     Bool_t   fRandomPosition;//flag indicating if positions of particles shuold be randomized according to parameters below
+     ERandomizationWay   fRandomPosition;//flag indicating if positions of particles shuold be randomized according to parameters below
      Double_t fRadius;        //radius used for randomizing particle vertex position
      
      Double_t fOneMinusLambda; //1 - intercept parameter
@@ -71,7 +93,13 @@ class AliHBTLLWeights: public AliHBTWeights
      AliHBTLLWeights(const AliHBTLLWeights &/*source*/);
      AliHBTLLWeights & operator=(const AliHBTLLWeights& /*source*/);
      
-     ClassDef(AliHBTLLWeights,1)
+     Int_t SetMomentaInLCMS(AliVAODParticle* part1, AliVAODParticle* part2);
+     void  RandomPairDistances();
+     
+          
+     void Rotate(Double_t x, Double_t y, Double_t sphi, Double_t cphi, Double_t& xr, Double_t& yr);
+     void Boost (Double_t z, Double_t t, Double_t beta, Double_t gamma, Double_t& zt, Double_t& yt);
+     ClassDef(AliHBTLLWeights,2)
  };
 
 #endif
index 3a27165..08579bb 100644 (file)
@@ -310,7 +310,7 @@ void AliHBTRndmGaussBall::Randomize(Double_t& x,Double_t& y,Double_t&z,Double_t&
   
   if (fTsigma == 0.0)
    {
-     t = 0.0;
+     t = fTmean;
      return;
    }
   
index 7d3af8d..6a212d4 100644 (file)
@@ -9,8 +9,8 @@ class AliHBTWeights: public TObject
  {
    public:
       virtual ~AliHBTWeights();
-      static Double_t Weight(const AliHBTPair* partpair){return (fgWeights)?fgWeights->GetWeight(partpair):0.0;}
-      virtual Double_t GetWeight(const AliHBTPair* partpair) = 0;
+      static Double_t Weight(AliHBTPair* partpair){return (fgWeights)?fgWeights->GetWeight(partpair):0.0;}
+      virtual Double_t GetWeight(AliHBTPair* partpair) = 0;
       virtual void Set() = 0;
       static AliHBTWeights* Instance() {return fgWeights;}
       
index c3fa195..caa8952 100644 (file)
@@ -1,4 +1,3 @@
-
        SUBROUTINE ltran12
 C
 
@@ -18,67 +17,26 @@ C-FSI ***************************************************
       COMMON/FSI_CVK/V,CVK
 C-FSI ***************************************************
       COMMON /PAIR/P12T,V12Z,GAMZ,V12T,CPHI,SPHI
-      
-      icrf=1
-      irot=1
-      
-C---> Particle energies ---------
+C   calculating Ri, Pi and Ei
+      R1=DSQRT(X1*X1+Y1*Y1+Z1*Z1)
+      R2=DSQRT(X2*X2+Y2*Y2+Z2*Z2)
       P1S=P1X*P1X+P1Y*P1Y+P1Z*P1Z
       P2S=P2X*P2X+P2Y*P2Y+P2Z*P2Z
-      E1=DSQRT(AM1*AM1+P1S)
-      E2=DSQRT(AM2*AM2+P2S)
-C---> Pair parameters -----------
-      E12=E1+E2       ! Energy
-      P12X=P1X+P2X    ! Px
-      P12Y=P1Y+P2Y    ! Py
-      P12Z=P1Z+P2Z    ! Pz
-      P12S=P12X**2+P12Y**2+P12Z**2
-      P12 =DSQRT(P12S)! Momentum
-      V12 =P12/E12    ! Velocity
-      CTH =P12Z/P12   ! cos(theta)
-      STH =DSQRT(1.D0-CTH**2) !sin
-      V12Z=V12*CTH    ! Longit. V
-      GAMZ=1.D0/DSQRT(1.D0-V12Z**2)
-C--      V12T=V12*STH    ! Transv. V in CMS (not needed)
-      P12TS=P12X*P12X+P12Y*P12Y
-      P12T=DSQRT(P12TS) !Pt
-C===> Azimuthal rotation (Pt||X) ============
-      IF(V12T.NE.0.D0) THEN
-        CPHI=P12X/P12T  ! cos(phi)
-        SPHI=P12Y/P12T  ! sin(phi)
-        IF((irot.eq.1)) THEN 
-          CALL ROT8(P1X,P1Y,SPHI,CPHI,P1X,P1Y)       
-          CALL ROT8(P2X,P2Y,SPHI,CPHI,P2X,P2Y)
-          CALL ROT8(X1,Y1,SPHI,CPHI,X1,Y1)       
-          CALL ROT8(X2,Y2,SPHI,CPHI,X2,Y2) 
-        END IF             
-      ELSE ! Rotation impossible 
-       CPHI=2.D0 ! to avoid
-       SPHI=2.D0 ! using it !  
-      END IF             
-C===> Co-moving ref. frame       ============
-      IF(icrf.gt.0) THEN
-        CALL LTR8(P1Z,E1,V12Z,GAMZ,P1Z,E1a)
-        CALL LTR8(P2Z,E2,V12Z,GAMZ,P2Z,E2a)
-        P1S=P1X*P1X+P1Y*P1Y+P1Z*P1Z
-        P2S=P2X*P2X+P2Y*P2Y+P2Z*P2Z
-        E1=DSQRT(AM1*AM1+P1S)
-        E2=DSQRT(AM2*AM2+P2S)
-        CALL LTR8(Z1,T1,V12Z,GAMZ,Z1,T1)
-        CALL LTR8(Z2,T2,V12Z,GAMZ,Z2,T2)
-      END IF        
-C===> Pair reference frame       ============
       P1=DSQRT(P1S)
       P2=DSQRT(P2S)
+      E1=DSQRT(AM1*AM1+P1S)
+      E2=DSQRT(AM2*AM2+P2S)
+C-----------------------------------------------------------------------
       E12=E1+E2
       P12X=P1X+P2X
       P12Y=P1Y+P2Y
       P12Z=P1Z+P2Z
       P12S=P12X**2+P12Y**2+P12Z**2
-      P12 =DSQRT(P12S)
-      AM12S=E12*E12-P12S
-      AM12=DSQRT(AM12S)
+      AM12=DSQRT(E12**2-P12S)
       EPM=E12+AM12
+      P12=DSQRT(P12S)
       P112=P1X*P12X+P1Y*P12Y+P1Z*P12Z
       H1=(P112/EPM-E1)/AM12
       PPX=P1X+P12X*H1
@@ -87,28 +45,47 @@ C===> Pair reference frame       ============
       EE=(E12*E1-P112)/AM12
       AKS=EE**2-AM1**2
       AK=DSQRT(AKS)
-      CVK=(P12X*PPX+P12Y*PPY+P12Z*PPZ)/(P12*AK)
-      V=P12/E12 
-      V12T=P12T/SQRT(AM12S+P12TS) ! transverse velocity in LCMS   
-C---> Coordinates -----------------------------
+
+      RETURN
+      END
+C====       
+C==== ===============================================================
+C==== ===============================================================
+C==== ===============================================================
+
+      subroutine BoostToPrf()
+      IMPLICIT REAL*8 (A-H,O-Z)
+      COMMON/FSI_CVK/V,CVK
+      COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1,  !part. momenta in NRF
+     1               P2X,P2Y,P2Z,E2,P2
+      COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS,
+     1               X,Y,Z,T,RP,RPS
+      COMMON/FSI_COOR/X1,Y1,Z1,T1,R1, ! 4-coord. of emis. points in NRF
+     1                X2,Y2,Z2,T2,R2
+      COMMON/FSI_P12/P12X,P12Y,P12Z,E12,P12,AM12,EPM
+
       XS=X1-X2
       YS=Y1-Y2
       ZS=Z1-Z2
-      TS=T1-T2  
+      TS=T1-T2
       RS12=XS*P12X+YS*P12Y+ZS*P12Z
       H1=(RS12/EPM-TS)/AM12
-
       X=XS+P12X*H1
       Y=YS+P12Y*H1
       Z=ZS+P12Z*H1
       T=(E12*TS-RS12)/AM12
       RPS=X*X+Y*Y+Z*Z
       RP=DSQRT(RPS)
+CW      WRITE(6,38)'RP ',RP,'X ',X,Y,Z,T
+38    FORMAT(A7,E11.4,A7,4E11.4)
+      CVK=(P12X*PPX+P12Y*PPY+P12Z*PPZ)/(P12*AK)
+      V=P12/E12
+      return 
+      end
+C==== ===============================================================
+C==== ===============================================================
 
-
-      RETURN
-      END
-C====       
       SUBROUTINE LTR8(Z,T,BETA,GAMMA,ZT,TT)
 C===> Lorentz Transf. of Z(Pz) and T(E) to moving ref. frame.(REAL*8) 
 CInp: Z,T-Zcoord,Time before tr., BETA,GAMMA- velocity, Lor.fact.
@@ -150,19 +127,7 @@ C=====Just sets distance between particles
       IMPLICIT REAL*8 (A-H,O-Z)
       COMMON/FSI_PRF/PPX,PPY,PPZ,AK,AKS, ! momenta in PRF
      1               X,Y,Z,T,RP,RPS
-      RPS=R
-      RP=R*R
+      RP=R
+      RPS=R*R
       RETURN
       END
-
-
-
-
-
-
-
-
-
-
-