New version of the code, see history.txt for more details
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 7 Jun 2000 18:44:11 +0000 (18:44 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 7 Jun 2000 18:44:11 +0000 (18:44 +0000)
RALICE/AliMath.cxx
RALICE/AliMath.h
RALICE/AliVertex.cxx
RALICE/AliVertex.h

index 0826a0a..5c78bd6 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.3  1999/11/03 14:23:17  fca
+New version of RALICE introduced
+
 Revision 1.2  1999/09/29 09:24:28  fca
 Introduction of the Copyright and cvs Log
 
@@ -36,6 +39,7 @@ Introduction of the Copyright and cvs Log
 //                             // correct model
 //
 //--- Author: Nick van Eijndhoven 14-nov-1998 UU-SAP Utrecht
+//- Modified: NvE 12-mar-2000 UU-SAP Utrecht
 ///////////////////////////////////////////////////////////////////////////
 
 #include "AliMath.h"
@@ -52,7 +56,7 @@ AliMath::~AliMath()
 // Destructor
 }
 ///////////////////////////////////////////////////////////////////////////
-Float_t AliMath::Gamma(Float_t z)
+Double_t AliMath::Gamma(Double_t z)
 {
 // Computation of gamma(z) for all z>0.
 //
@@ -69,11 +73,11 @@ Float_t AliMath::Gamma(Float_t z)
   return 0;
  }
  
Float_t v=LnGamma(z);
Double_t v=LnGamma(z);
  return exp(v);
 }
 ///////////////////////////////////////////////////////////////////////////
-Float_t AliMath::Gamma(Float_t a,Float_t x)
+Double_t AliMath::Gamma(Double_t a,Double_t x)
 {
 // Computation of the incomplete gamma function P(a,x)
 //
@@ -104,7 +108,7 @@ Float_t AliMath::Gamma(Float_t a,Float_t x)
  }
 }
 ///////////////////////////////////////////////////////////////////////////
-Float_t AliMath::LnGamma(Float_t z)
+Double_t AliMath::LnGamma(Double_t z)
 {
 // Computation of ln[gamma(z)] for all z>0.
 //
@@ -143,11 +147,11 @@ Float_t AliMath::LnGamma(Float_t z)
   y+=1.;
   ser+=c[i]/y;
  }
Float_t v=tmp+log(c[0]*ser/x);
Double_t v=tmp+log(c[0]*ser/x);
  return v;
 }
 ///////////////////////////////////////////////////////////////////////////
-Float_t AliMath::GamSer(Float_t a,Float_t x)
+Double_t AliMath::GamSer(Double_t a,Double_t x)
 {
 // Computation of the incomplete gamma function P(a,x)
 // via its series representation.
@@ -158,7 +162,7 @@ Float_t AliMath::GamSer(Float_t a,Float_t x)
 //--- Nve 14-nov-1998 UU-SAP Utrecht
  
  Int_t itmax=100;   // Maximum number of iterations
Float_t eps=3.e-7; // Relative accuracy
Double_t eps=3.e-7; // Relative accuracy
  
  if (a<=0.)
  {
@@ -172,10 +176,10 @@ Float_t AliMath::GamSer(Float_t a,Float_t x)
   return 0;
  }
  
Float_t gln=LnGamma(a);
Float_t ap=a;
Float_t sum=1./a;
Float_t del=sum;
Double_t gln=LnGamma(a);
Double_t ap=a;
Double_t sum=1./a;
Double_t del=sum;
  for (Int_t n=1; n<=itmax; n++)
  {
   ap+=1.;
@@ -184,11 +188,11 @@ Float_t AliMath::GamSer(Float_t a,Float_t x)
   if (fabs(del)<fabs(sum*eps)) break;
   if (n==itmax) cout << "*GamSer(a,x)* a too large or itmax too small" << endl;
  }
Float_t v=sum*exp(-x+a*log(x)-gln);
Double_t v=sum*exp(-x+a*log(x)-gln);
  return v;
 }
 ///////////////////////////////////////////////////////////////////////////
-Float_t AliMath::GamCf(Float_t a,Float_t x)
+Double_t AliMath::GamCf(Double_t a,Double_t x)
 {
 // Computation of the incomplete gamma function P(a,x)
 // via its continued fraction representation.
@@ -199,8 +203,8 @@ Float_t AliMath::GamCf(Float_t a,Float_t x)
 //--- Nve 14-nov-1998 UU-SAP Utrecht
  
  Int_t itmax=100;      // Maximum number of iterations
Float_t eps=3.e-7;    // Relative accuracy
Float_t fpmin=1.e-30; // Smallest Float_t value allowed here
Double_t eps=3.e-7;    // Relative accuracy
Double_t fpmin=1.e-30; // Smallest Double_t value allowed here
  
  if (a<=0.)
  {
@@ -214,15 +218,15 @@ Float_t AliMath::GamCf(Float_t a,Float_t x)
   return 0;
  }
  
Float_t gln=LnGamma(a);
Float_t b=x+1.-a;
Float_t c=1./fpmin;
Float_t d=1./b;
Float_t h=d;
Float_t an,del;
Double_t gln=LnGamma(a);
Double_t b=x+1.-a;
Double_t c=1./fpmin;
Double_t d=1./b;
Double_t h=d;
Double_t an,del;
  for (Int_t i=1; i<=itmax; i++)
  {
-  an=float(-i)*(float(i)-a);
+  an=double(-i)*(double(i)-a);
   b+=2.;
   d=an*d+b;
   if (fabs(d)<fpmin) d=fpmin;
@@ -234,11 +238,11 @@ Float_t AliMath::GamCf(Float_t a,Float_t x)
   if (fabs(del-1.)<eps) break;
   if (i==itmax) cout << "*GamCf(a,x)* a too large or itmax too small" << endl;
  }
Float_t v=exp(-x+a*log(x)-gln)*h;
Double_t v=exp(-x+a*log(x)-gln)*h;
  return (1.-v);
 }
 ///////////////////////////////////////////////////////////////////////////
-Float_t AliMath::Erf(Float_t x)
+Double_t AliMath::Erf(Double_t x)
 {
 // Computation of the error function erf(x).
 //
@@ -247,7 +251,7 @@ Float_t AliMath::Erf(Float_t x)
  return (1.-Erfc(x));
 }
 ///////////////////////////////////////////////////////////////////////////
-Float_t AliMath::Erfc(Float_t x)
+Double_t AliMath::Erfc(Double_t x)
 {
 // Computation of the complementary error function erfc(x).
 //
@@ -259,19 +263,19 @@ Float_t AliMath::Erfc(Float_t x)
 //--- Nve 14-nov-1998 UU-SAP Utrecht
  
  // The parameters of the Chebyshev fit
- const Float_t a1=-1.26551223,  a2=1.00002368,
-               a3= 0.37409196,  a4=0.09678418,
-               a5=-0.18628806,  a6=0.27886807,
-               a7=-1.13520398,  a8=1.48851587,
-               a9=-0.82215223, a10=0.17087277;
+ const Double_t a1=-1.26551223,  a2=1.00002368,
+                a3= 0.37409196,  a4=0.09678418,
+                a5=-0.18628806,  a6=0.27886807,
+                a7=-1.13520398,  a8=1.48851587,
+                a9=-0.82215223, a10=0.17087277;
  
Float_t v=1.; // The return value
Double_t v=1.; // The return value
  
Float_t z=fabs(x);
Double_t z=fabs(x);
  
  if (z <= 0.) return v; // erfc(0)=1
  
Float_t t=1./(1.+0.5*z);
Double_t t=1./(1.+0.5*z);
  
  v=t*exp((-z*z)
          +a1+t*(a2+t*(a3+t*(a4+t*(a5+t*(a6+t*(a7+t*(a8+t*(a9+t*a10)))))))));
@@ -281,7 +285,7 @@ Float_t AliMath::Erfc(Float_t x)
  return v;
 }
 ///////////////////////////////////////////////////////////////////////////
-Float_t AliMath::Prob(Float_t chi2,Int_t ndf)
+Double_t AliMath::Prob(Double_t chi2,Int_t ndf)
 {
 // Computation of the probability for a certain Chi-squared (chi2)
 // and number of degrees of freedom (ndf).
@@ -316,22 +320,272 @@ Float_t AliMath::Prob(Float_t chi2,Int_t ndf)
 // This code may be activated in case the gamma function gives problems
 // if (ndf==1)
 // {
-//  Float_t v=1.-Erf(sqrt(chi2)/sqrt(2.));
+//  Double_t v=1.-Erf(sqrt(chi2)/sqrt(2.));
 //  return v;
 // }
  
 // Gaussian approximation for large ndf
 // This code may be activated in case the gamma function shows a problem
-// Float_t q=sqrt(2.*chi2)-sqrt(float(2*ndf-1));
+// Double_t q=sqrt(2.*chi2)-sqrt(double(2*ndf-1));
 // if (n>30 && q>0.)
 // {
-//  Float_t v=0.5*(1.-Erf(q/sqrt(2.)));
+//  Double_t v=0.5*(1.-Erf(q/sqrt(2.)));
 //  return v;
 // }
  
  // Evaluate the incomplete gamma function
Float_t a=float(ndf)/2.;
Float_t x=chi2/2.;
Double_t a=double(ndf)/2.;
Double_t x=chi2/2.;
  return (1.-Gamma(a,x));
 }
 ///////////////////////////////////////////////////////////////////////////
+Double_t AliMath::BesselI0(Double_t x)
+{
+// Computation of the modified Bessel function I_0(x) for any real x.
+//
+// The algorithm is based on the article by Abramowitz and Stegun [1]
+// as denoted in Numerical Recipes 2nd ed. on p. 230 (W.H.Press et al.).
+//
+// [1] M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
+//     Applied Mathematics Series vol. 55 (1964), Washington.  
+//
+//--- NvE 12-mar-2000 UU-SAP Utrecht
+
+ // Parameters of the polynomial approximation  
+ const Double_t p1=1.0,          p2=3.5156229,    p3=3.0899424,
+                p4=1.2067492,    p5=0.2659732,    p6=3.60768e-2,  p7=4.5813e-3;
+
+ const Double_t q1= 0.39894228,  q2= 1.328592e-2, q3= 2.25319e-3,
+                q4=-1.57565e-3,  q5= 9.16281e-3,  q6=-2.057706e-2,
+                q7= 2.635537e-2, q8=-1.647633e-2, q9= 3.92377e-3; 
+
+ Double_t ax=fabs(x);
+
+ Double_t y=0,result=0;
+
+ if (ax < 3.75)
+ {
+  y=pow(x/3.75,2);
+  result=p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))));
+ }
+ else
+ {
+  y=3.75/ax;
+  result=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9))))))));
+ }
+
+ return result;
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t AliMath::BesselK0(Double_t x)
+{
+// Computation of the modified Bessel function K_0(x) for positive real x.
+//
+// The algorithm is based on the article by Abramowitz and Stegun [1]
+// as denoted in Numerical Recipes 2nd ed. on p. 230 (W.H.Press et al.).
+//
+// [1] M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
+//     Applied Mathematics Series vol. 55 (1964), Washington.  
+//
+//--- NvE 12-mar-2000 UU-SAP Utrecht
+
+ // Parameters of the polynomial approximation  
+ const Double_t p1=-0.57721566,  p2=0.42278420,   p3=0.23069756,
+                p4= 3.488590e-2, p5=2.62698e-3,   p6=1.0750e-4,    p7=7.4e-5;
+
+ const Double_t q1= 1.25331414,  q2=-7.832358e-2, q3= 2.189568e-2,
+                q4=-1.062446e-2, q5= 5.87872e-3,  q6=-2.51540e-3,  q7=5.3208e-4;
+
+ if (x <= 0)
+ {
+  cout << " *BesselK0* Invalid argument x = " << x << endl;
+  return 0;
+ }
+
+ Double_t y=0,result=0;
+
+ if (x <= 2)
+ {
+  y=x*x/4.;
+  result=(-log(x/2.)*BesselI0(x))+(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))));
+ }
+ else
+ {
+  y=2./x;
+  result=(exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*q7))))));
+ }
+
+ return result;
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t AliMath::BesselI1(Double_t x)
+{
+// Computation of the modified Bessel function I_1(x) for any real x.
+//
+// The algorithm is based on the article by Abramowitz and Stegun [1]
+// as denoted in Numerical Recipes 2nd ed. on p. 230 (W.H.Press et al.).
+//
+// [1] M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
+//     Applied Mathematics Series vol. 55 (1964), Washington.  
+//
+//--- NvE 12-mar-2000 UU-SAP Utrecht
+
+ // Parameters of the polynomial approximation  
+ const Double_t p1=0.5,          p2=0.87890594,   p3=0.51498869,
+                p4=0.15084934,   p5=2.658733e-2,  p6=3.01532e-3,  p7=3.2411e-4;
+
+ const Double_t q1= 0.39894228,  q2=-3.988024e-2, q3=-3.62018e-3,
+                q4= 1.63801e-3,  q5=-1.031555e-2, q6= 2.282967e-2,
+                q7=-2.895312e-2, q8= 1.787654e-2, q9=-4.20059e-3; 
+
+ Double_t ax=fabs(x);
+
+ Double_t y=0,result=0;
+
+ if (ax < 3.75)
+ {
+  y=pow(x/3.75,2);
+  result=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))));
+ }
+ else
+ {
+  y=3.75/ax;
+  result=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9))))))));
+  if (x < 0) result=-result;
+ }
+
+ return result;
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t AliMath::BesselK1(Double_t x)
+{
+// Computation of the modified Bessel function K_1(x) for positive real x.
+//
+// The algorithm is based on the article by Abramowitz and Stegun [1]
+// as denoted in Numerical Recipes 2nd ed. on p. 230 (W.H.Press et al.).
+//
+// [1] M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
+//     Applied Mathematics Series vol. 55 (1964), Washington.  
+//
+//--- NvE 12-mar-2000 UU-SAP Utrecht
+
+ // Parameters of the polynomial approximation  
+ const Double_t p1= 1.,          p2= 0.15443144,  p3=-0.67278579,
+                p4=-0.18156897,  p5=-1.919402e-2, p6=-1.10404e-3,  p7=-4.686e-5;
+
+ const Double_t q1= 1.25331414,  q2= 0.23498619,  q3=-3.655620e-2,
+                q4= 1.504268e-2, q5=-7.80353e-3,  q6= 3.25614e-3,  q7=-6.8245e-4;
+
+ if (x <= 0)
+ {
+  cout << " *BesselK1* Invalid argument x = " << x << endl;
+  return 0;
+ }
+
+ Double_t y=0,result=0;
+
+ if (x <= 2)
+ {
+  y=x*x/4.;
+  result=(log(x/2.)*BesselI1(x))+(1./x)*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))));
+ }
+ else
+ {
+  y=2./x;
+  result=(exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*q7))))));
+ }
+
+ return result;
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t AliMath::BesselK(Int_t n,Double_t x)
+{
+// Computation of the Integer Order Modified Bessel function K_n(x)
+// for n=0,1,2,... and positive real x.
+//
+// The algorithm uses the recurrence relation
+//
+//               K_n+1(x) = (2n/x)*K_n(x) + K_n-1(x) 
+//
+// as denoted in Numerical Recipes 2nd ed. on p. 232 (W.H.Press et al.).
+//
+//--- NvE 12-mar-2000 UU-SAP Utrecht
+
+ if (x <= 0 || n < 0)
+ {
+  cout << " *BesselK* Invalid argument(s) (n,x) = (" << n << " , " << x << ")" << endl;
+  return 0;
+ }
+
+ if (n==0) return BesselK0(x);
+
+ if (n==1) return BesselK1(x);
+
+ // Perform upward recurrence for all x
+ Double_t tox=2./x;
+ Double_t bkm=BesselK0(x);
+ Double_t bk=BesselK1(x);
+ Double_t bkp=0;
+ for (Int_t j=1; j<n; j++)
+ {
+  bkp=bkm+double(j)*tox*bk;
+  bkm=bk;
+  bk=bkp;
+ }
+
+ return bk;
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t AliMath::BesselI(Int_t n,Double_t x)
+{
+// Computation of the Integer Order Modified Bessel function I_n(x)
+// for n=0,1,2,... and any real x.
+//
+// The algorithm uses the recurrence relation
+//
+//               I_n+1(x) = (-2n/x)*I_n(x) + I_n-1(x) 
+//
+// as denoted in Numerical Recipes 2nd ed. on p. 232 (W.H.Press et al.).
+//
+//--- NvE 12-mar-2000 UU-SAP Utrecht
+
+ Int_t iacc=40; // Increase to enhance accuracy
+ Double_t bigno=1.e10, bigni=1.e-10;
+
+ if (n < 0)
+ {
+  cout << " *BesselI* Invalid argument (n,x) = (" << n << " , " << x << ")" << endl;
+  return 0;
+ }
+
+ if (n==0) return BesselI0(x);
+
+ if (n==1) return BesselI1(x);
+
+ if (fabs(x) < 1.e-10) return 0;
+
+ Double_t tox=2./fabs(x);
+ Double_t bip=0,bim=0;
+ Double_t bi=1;
+ Double_t result=0;
+ Int_t m=2*((n+int(sqrt(float(iacc*n))))); // Downward recurrence from even m
+ for (Int_t j=m; j<=1; j--)
+ {
+  bim=bip+double(j)*tox*bi;
+  bip=bi;
+  bi=bim;
+  if (fabs(bi) > bigno) // Renormalise to prevent overflows
+  {
+   result*=bigni;
+   bi*=bigni;
+   bip*=bigni;
+  }
+  if (j==n) result=bip;
+ }
+
+ result*=BesselI0(x)/bi; // Normalise with I0(x)
+ if ((x < 0) && (n%2 == 1)) result=-result;
+
+ return result;
+}
+///////////////////////////////////////////////////////////////////////////
index 4b3485d..8560b5c 100644 (file)
 class AliMath : public TObject
 {
  public:
-  AliMath();                            // Default constructor
-  ~AliMath();                           // Destructor
-  Float_t Gamma(Float_t z);             // Standard gamma function Gamma(z)
-  Float_t Gamma(Float_t a,Float_t x);   // Incomplete gamma function P(a,x)
-  Float_t LnGamma(Float_t z);           // Compute ln[Gamma(z)]
-  Float_t Erf(Float_t x);               // Error function erf(x)
-  Float_t Erfc(Float_t x);              // Complementary error function erfc(x)
-  Float_t Prob(Float_t chi2,Int_t ndf); // Compute Chi-squared probability
+  AliMath();                              // Default constructor
+  ~AliMath();                             // Destructor
+  Double_t Gamma(Double_t z);             // Standard gamma function Gamma(z)
+  Double_t Gamma(Double_t a,Double_t x);  // Incomplete gamma function P(a,x)
+  Double_t LnGamma(Double_t z);           // Compute ln[Gamma(z)]
+  Double_t Erf(Double_t x);               // Error function erf(x)
+  Double_t Erfc(Double_t x);              // Complementary error function erfc(x)
+  Double_t Prob(Double_t chi2,Int_t ndf); // Compute Chi-squared probability
+  Double_t BesselI(Int_t n,Double_t x);   // Compute integer order modified Bessel function I_n(x)
+  Double_t BesselK(Int_t n,Double_t x);   // Compute integer order modified Bessel function K_n(x)
  
  protected:
-  Float_t GamSer(Float_t a,Float_t x);  // Compute P(a,x) via serial representation
-  Float_t GamCf(Float_t a,Float_t x);   // Compute P(a,x) via continued fractions
+  Double_t GamSer(Double_t a,Double_t x); // Compute P(a,x) via serial representation
+  Double_t GamCf(Double_t a,Double_t x);  // Compute P(a,x) via continued fractions
+  Double_t BesselI0(Double_t x);          // Compute modified Bessel function I_0(x)
+  Double_t BesselK0(Double_t x);          // Compute modified Bessel function K_0(x)
+  Double_t BesselI1(Double_t x);          // Compute modified Bessel function I_1(x)
+  Double_t BesselK1(Double_t x);          // Compute modified Bessel function K_1(x)
  
  ClassDef(AliMath,1) // Various mathematical tools for physics analysis.
  
index a670151..41b61d4 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.4  1999/11/03 14:23:18  fca
+New version of RALICE introduced
+
 Revision 1.3  1999/09/29 09:24:28  fca
 Introduction of the Copyright and cvs Log
 
@@ -123,6 +126,7 @@ AliVertex::AliVertex()
 // Initial maximum number of sec. vertices is set to the default value.
  fNvmax=0;
  fVertices=0;
+ fConnects=0;
  Reset();
  SetNtinit();
  SetNvmax();
@@ -134,6 +138,7 @@ AliVertex::AliVertex(Int_t n)
 // All variables initialised to 0
  fNvmax=0;
  fVertices=0;
+ fConnects=0;
  Reset();
  if (n > 0)
  {
@@ -155,6 +160,12 @@ AliVertex::~AliVertex()
 // Default destructor
  if (fVertices) delete fVertices;
  fVertices=0;
+ if (fConnects)
+ {
+  fConnects->Delete();
+  delete fConnects;
+  fConnects=0;
+ }
 }
 ///////////////////////////////////////////////////////////////////////////
 void AliVertex::SetNvmax(Int_t n)
@@ -182,6 +193,12 @@ void AliVertex::Reset()
 
  fNvtx=0;
  if (fNvmax>0) SetNvmax(fNvmax);
+ if (fConnects)
+ {
+  fConnects->Delete();
+  delete fConnects;
+  fConnects=0;
+ }
 }
 ///////////////////////////////////////////////////////////////////////////
 void AliVertex::Add(AliJet& j)
@@ -195,26 +212,54 @@ void AliVertex::Add(AliJet& j)
  }
 }
 ///////////////////////////////////////////////////////////////////////////
-void AliVertex::Add(AliVertex& v)
+void AliVertex::Add(AliVertex& v,Int_t connect)
 {
 // Add a (secondary) vertex to the current vertex.
 // In case the maximum number of (secondary) vertices has been reached,
 // the array space will be extended automatically
 //
-// Note : The 4-momentum of the current (primary) vertex
-//        is updated automatically, but the track connecting
-//        both vertices has to be entered separately by the user.
+// Note : By default the 4-momentum and charge of the current (primary) vertex
+//        are updated by automatically creating the track connecting
+//        both vertices. The track parameters are taken from the
+//        4-momentum and charge of the secondary vertex.
+//        The automatic creation of the connecting track and updating
+//        of the (primary) vertex 4-momentum and charge can be suppressed
+//        by specifying connect=0. In this case, however, the user
+//        has to introduce the connecting track lateron by hand
+//        explicitly in order to match the kinematics and charge.
 //
  if (fNvtx == fNvmax) // Check if maximum vertex number is reached
  {
   fNvmax++;
   fVertices->Expand(fNvmax);
  }
- // Update 4-momentum for current vertex
+
+ // Add the linked (secondary) vertex to the list 
  fNvtx++;
  fVertices->Add(&v);
- (Ali4Vector)(*this)+=v;
+
+ // Create connecting track and update 4-momentum and charge for current vertex
+ if (connect)
+ {
+  AliPosition r1=GetPosition();
+  AliPosition r2=v.GetPosition();
+  Float_t q=v.GetCharge();
+  Ali3Vector p=v.Get3Momentum();
+  Double_t v2=v.GetInvariant();
+  Double_t dv2=v.Ali4Vector::GetResultError();
+
+  AliTrack* t=new AliTrack;
+  t->SetBeginPoint(r1);
+  t->SetEndPoint(r2);
+  t->SetCharge(q);
+  t->Set3Momentum(p);
+  t->SetInvariant(v2,dv2);
+
+  AliJet::Add(t);
+
+  if (!fConnects) fConnects=new TObjArray(fNvmax);
+  fConnects->Add(t);
+ }
 }
 ///////////////////////////////////////////////////////////////////////////
 void AliVertex::Info(TString f)
index 3c36fbf..3d5d463 100644 (file)
 class AliVertex : public AliJet,public AliPosition
 {
  public:
-  AliVertex();                        // Default constructor
-  AliVertex(Int_t n);                 // Create a vertex to hold initially n tracks
-  ~AliVertex();                       // Default destructor
-  void Reset();                       // Reset all values
-  void Add(AliJet& j);                // Add a jet of tracks to the vertex
-  void Add(AliVertex& v);             // Add a (secondary) vertex to the current vertex
+  AliVertex();                            // Default constructor
+  AliVertex(Int_t n);                     // Create a vertex to hold initially n tracks
+  ~AliVertex();                           // Default destructor
+  void Reset();                           // Reset all values
+  void Add(AliJet& j);                    // Add a jet of tracks to the vertex
+  void Add(AliVertex& v,Int_t connect=1); // Add (and connect) a (sec.) vertex to the current vertex
   void Add(AliJet* j)    { Add(*j); }
-  void Add(AliVertex* v) { Add(*v); }
-  void Info(TString f="car");         // Print the vertex info within coordinate frame f
-  void List(TString f="car");         // Print vertex prim. track information for coord. frame f
-  void ListAll(TString f="car");      // Print prim. + sec. vertex full track info for coord. frame f
-  Int_t GetNvertices();               // Return the number of (secondary) vertices
-  AliVertex* GetVertex(Int_t i);      // Provide i-th (secondary) vertex
-  void SetNvmax(Int_t n=2);           // Set the initial max. number of (secondary) vertices
+  void Add(AliVertex* v,Int_t connect=1) { Add(*v,connect); }
+  void Info(TString f="car");             // Print the vertex info within coordinate frame f
+  void List(TString f="car");             // Print vertex prim. track information for coord. frame f
+  void ListAll(TString f="car");          // Print prim. + sec. vertex full track info for coord. frame f
+  Int_t GetNvertices();                   // Return the number of (secondary) vertices
+  AliVertex* GetVertex(Int_t i);          // Provide i-th (secondary) vertex
+  void SetNvmax(Int_t n=2);               // Set the initial max. number of (secondary) vertices
 
  protected:
   Int_t fNvmax;         // The maximum number of (secondary) vertices
   Int_t fNvtx;          // The number of (secondary) vertices
   TObjArray* fVertices; // Array to hold the pointers to the (secondary) vertices
+  TObjArray* fConnects; // Array to hold the pointers to the auto-generated connecting tracks
 
  private:
   void Dump(AliVertex* v,Int_t n,TString f); // Recursively print all sec. vertices