RALICE/AliBoost.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 // $Id$
  17
  18 ///////////////////////////////////////////////////////////////////////////
  19 // Class AliBoost
  20 // Perform various Lorentz transformations.
  21 //
  22 // Example :
  23 // =========
  24 //
  25 // Float_t a[3]={0.1,0.2,0.3};
  26 // Float_t ea[3]={0.01,0.02,0.03};
  27 // Ali3Vector beta;
  28 // beta.SetVector(a,"car");
  29 // beta.SetErrors(ea,"car");
  30 //
  31 // AliBoost b1;
  32 // b1.SetBeta(beta);
  33 // b1.Data();
  34 //
  35 // Float_t b[4]={14,1,2,3};
  36 // Float_t eb[4]={1.4,0.1,0.2,0.3};
  37 // Ali4Vector p;
  38 // p.SetVector(b,"car");
  39 // p.SetErrors(eb,"car");
  40 // Ali4Vector pprim=b1.Boost(p);
  41 // p.Data();
  42 // pprim.Data();
  43 //
  44 // p=b1.Inverse(pprim);
  45 // pprim.Data();
  46 // p.Data();
  47 //
  48 // Float_t c[4]={5,0,0,4};
  49 // Float_t ec[4]={0.5,0,0,0.4};
  50 // Ali4Vector q;
  51 // q.SetVector(c,"car");
  52 // q.SetErrors(ec,"car");
  53 //
  54 // AliBoost b2;
  55 // b2.Set4Momentum(q);
  56 // b2.Data("sph");
  57 //
  58 //--- Author: Nick van Eijndhoven 14-may-1996 UU-SAP Utrecht
  59 //- Modified: NvE $Date$ UU-SAP Utrecht
  60 ///////////////////////////////////////////////////////////////////////////
  61
  62 #include "AliBoost.h"
  63 #include "Riostream.h"
  64
  65 ClassImp(AliBoost) // Class implementation to enable ROOT I/O
  66
  67 AliBoost::AliBoost() : TObject()
  68 {
  69 // Creation of a Lorentz boost object and initialisation of parameters.
  70 // Beta is set to (0,0,0) and consequently Gamma=1.
  71 // All errors are initialised to 0.
  72  Double_t a[3]={0,0,0};
  73  fBeta.SetVector(a,"sph");
  74  fGamma=1;
  75  fDgamma=0;
  76  fDresult=0;
  77 }
  78 ///////////////////////////////////////////////////////////////////////////
  79 AliBoost::~AliBoost()
  80 {
  81 // Default destructor.
  82 }
  83 ///////////////////////////////////////////////////////////////////////////
  84 AliBoost::AliBoost(AliBoost& b) : TObject(b)
  85 {
  86 // Copy constructor
  87  fBeta=b.fBeta;
  88  fGamma=b.fGamma;
  89  fDgamma=b.fDgamma;
  90  fDresult=b.fDresult;
  91 }
  92 ///////////////////////////////////////////////////////////////////////////
  93 void AliBoost::SetBeta(Ali3Vector& b)
  94 {
  95 // Setting of boost parameters on basis of beta 3-vector.
  96 // The errors on the beta 3-vector are taken from the input 3-vector.
  97 // The gamma value and its error are calculated accordingly.
  98  fBeta=b;
  99  Double_t beta2=fBeta.Dot(fBeta);
 100  Double_t dbeta2=fBeta.GetResultError();
 101
 102  if (beta2 > 1.)
 103  {
 104   cout << " *AliBoost::SetBeta* beta > 1." << endl;
 105  }
 106  fGamma=0;
 107  fDgamma=0;
 108  Double_t temp=1.-beta2;
 109  if (temp > 0.)
 110  {
 111   fGamma=sqrt(1./temp);
 112   fDgamma=fabs(dbeta2/(2.*pow(temp,1.5)));
 113  }
 114 }
 115 ///////////////////////////////////////////////////////////////////////////
 116 void AliBoost::Set4Momentum(Ali4Vector& p)
 117 {
 118 // Setting of boost parameters on basis of momentum 4-vector data.
 119 // The errors of the input 4-vector are used to calculate the
 120 // errors on the beta 3-vector and the gamma factor.
 121  Double_t E=p.GetScalar();
 122  Double_t dE=p.GetResultError();
 123  if (E <= 0.)
 124  {
 125   cout << " *AliBoost::Set4Momentum* Unphysical situation." << endl;
 126   p.Data();
 127  }
 128  else
 129  {
 130   Ali3Vector b=p.Get3Vector();
 131   Double_t vb[3],eb[3];
 132   b.GetVector(vb,"car");
 133   b.GetErrors(eb,"car");
 134   b=b/E;
 135   for (Int_t i=0; i<3; i++)
 136   {
 137    eb[i]=sqrt(pow(eb[i]/E,2)+pow(vb[i]*dE/(E*E),2));
 138   }
 139   b.SetErrors(eb,"car");
 140   SetBeta(b);
 141  }
 142 }
 143 ///////////////////////////////////////////////////////////////////////////
 144 Ali3Vector AliBoost::GetBetaVector()
 145 {
 146 // Provide the beta 3-vector.
 147  return fBeta;
 148 }
 149 ///////////////////////////////////////////////////////////////////////////
 150 Double_t AliBoost::GetBeta()
 151 {
 152 // Provide the norm of the beta 3-vector.
 153 // The error on the value can be obtained via GetResultError().
 154  Double_t norm=fBeta.GetNorm();
 155  fDresult=fBeta.GetResultError();
 156  return norm;
 157 }
 158 ///////////////////////////////////////////////////////////////////////////
 159 Double_t AliBoost::GetGamma()
 160 {
 161 // Provide the gamma factor.
 162 // The error on the value can be obtained via GetResultError().
 163  fDresult=fDgamma;
 164  return fGamma;
 165 }
 166 ///////////////////////////////////////////////////////////////////////////
 167 Double_t AliBoost::GetResultError()
 168 {
 169 // Provide the error on the result of an operation yielding a scalar.
 170 // E.g. GetBeta() or GetGamma()
 171  return fDresult;
 172 }
 173 ///////////////////////////////////////////////////////////////////////////
 174 void AliBoost::Data(TString f)
 175 {
 176 // Printing of the boost parameter info in coordinate frame f.
 177  Double_t beta=fBeta.GetNorm();
 178  Double_t dbeta=fBeta.GetResultError();
 179  cout << " *AliBoost::Data* beta : " << beta << " error : " << dbeta
 180       << " gamma : " << fGamma << " error : " << fDgamma << endl;
 181  cout << "  Beta";
 182  fBeta.Data(f);
 183 }
 184 ///////////////////////////////////////////////////////////////////////////
 185 Ali4Vector AliBoost::Boost(Ali4Vector& v)
 186 {
 187 // Perform the Lorentz boost on the 4-vector v.
 188 // Error propagation is performed automatically.
 189 // Note : As an approximation Beta and p.Dot(Beta) are considered as
 190 //        independent quantities.
 191
 192  Double_t beta=fBeta.GetNorm();
 193  Double_t dbeta=fBeta.GetResultError();
 194
 195  Double_t beta2=pow(beta,2);
 196
 197  if (beta > 1.e-10)
 198  {
 199   Double_t E=v.GetScalar();
 200   Double_t dE=v.GetResultError();
 201
 202   Ali3Vector p=v.Get3Vector();
 203
 204   Double_t pdotbeta=p.Dot(fBeta);
 205   Double_t dpdotbeta=p.GetResultError();
 206
 207   // Determine the new vector components
 208   Double_t Eprim=fGamma*(E-pdotbeta);
 209
 210   Double_t z=((fGamma-1.)*pdotbeta/beta2)-fGamma*E;
 211   Ali3Vector add=fBeta*z;
 212
 213   // Determine errors on the new vector components
 214   Double_t dEprim=sqrt(pow((E-pdotbeta)*fDgamma,2)+pow(fGamma*dE,2)
 215                       +pow(fGamma*dpdotbeta,2));
 216   Double_t dz=sqrt( pow(((fGamma-1.)/beta2)*dpdotbeta,2) + pow(fGamma*dE,2)
 217                    +pow((
 218     ((2./beta)-(4.*pow(beta,3)-6.*pow(beta,5))/(2.*pow((pow(beta,4)-pow(beta,6)),1.5)))*pdotbeta
 219     +beta*E/pow(fGamma,3))*dbeta,2) );
 220
 221   Double_t vb[3],eb[3];
 222   fBeta.GetVector(vb,"car");
 223   fBeta.GetErrors(eb,"car");
 224   for (Int_t i=0; i<3; i++)
 225   {
 226    eb[i]=sqrt(pow(z*eb[i],2)+pow(vb[i]*dz,2));
 227   }
 228   add.SetErrors(eb,"car");
 229
 230   // Calculate the new 3-vector
 231   Ali3Vector pprim=p+add;
 232
 233   // Set the components and errors of the new 4-vector
 234   Ali4Vector w;
 235   w.SetVector(Eprim,pprim);
 236   w.SetScalarError(dEprim);
 237
 238   return w;
 239  }
 240  else
 241  {
 242   return v;
 243  }
 244 }
 245 ///////////////////////////////////////////////////////////////////////////
 246 Ali4Vector AliBoost::Inverse(Ali4Vector& vprim)
 247 {
 248 // Perform the inverse Lorentz boost on the 4-vector vprim.
 249 // Error propagation is performed automatically.
 250 // Note : As an approximation Beta and pprim.Dot(Beta) are considered as
 251 //        independent quantities.
 252
 253  Double_t beta=fBeta.GetNorm();
 254  Double_t dbeta=fBeta.GetResultError();
 255
 256  Double_t beta2=pow(beta,2);
 257
 258  if (beta > 1.e-10)
 259  {
 260   Double_t Eprim=vprim.GetScalar();
 261   Double_t dEprim=vprim.GetResultError();
 262
 263   Ali3Vector pprim=vprim.Get3Vector();
 264
 265   Double_t pprimdotbeta=pprim.Dot(fBeta);
 266   Double_t dpprimdotbeta=pprim.GetResultError();
 267
 268   // Determine the new vector components
 269   Double_t E=fGamma*(Eprim+pprimdotbeta);
 270
 271   Double_t z=((fGamma-1.)*pprimdotbeta/beta2)+fGamma*Eprim;
 272   Ali3Vector add=fBeta*z;
 273
 274   // Determine errors on the prime-vector components
 275   Double_t dE=sqrt(pow((Eprim+pprimdotbeta)*fDgamma,2)+pow(fGamma*dEprim,2)
 276                       +pow(fGamma*dpprimdotbeta,2));
 277   Double_t dz=sqrt( pow(((fGamma-1.)/beta2)*dpprimdotbeta,2) + pow(fGamma*dEprim,2)
 278                    +pow((
 279     ((2./beta)-(4.*pow(beta,3)-6.*pow(beta,5))/(2.*pow((pow(beta,4)-pow(beta,6)),1.5)))*pprimdotbeta
 280     -beta*Eprim/pow(fGamma,3))*dbeta,2) );
 281
 282   Double_t vb[3],eb[3];
 283   fBeta.GetVector(vb,"car");
 284   fBeta.GetErrors(eb,"car");
 285   for (Int_t i=0; i<3; i++)
 286   {
 287    eb[i]=sqrt(pow(z*eb[i],2)+pow(vb[i]*dz,2));
 288   }
 289   add.SetErrors(eb,"car");
 290
 291   // Calculate the new 3-vector
 292   Ali3Vector p=pprim+add;
 293
 294   // Set the components and errors of the new 4-vector
 295   Ali4Vector w;
 296   w.SetVector(E,p);
 297   w.SetScalarError(dE);
 298
 299   return w;
 300  }
 301  else
 302  {
 303   return vprim;
 304  }
 305 }
 306 ///////////////////////////////////////////////////////////////////////////