which included commits to RCS files with non-trunk default branches.
--- /dev/null
+#include "Ali3Vector.h"
+
+ClassImp(Ali3Vector) // Class implementation to enable ROOT I/O
+
+Ali3Vector::Ali3Vector()
+{
+// Creation of an Ali3Vector object and initialisation of parameters
+ fV=0;
+ fTheta=0;
+ fPhi=0;
+}
+///////////////////////////////////////////////////////////////////////////
+Ali3Vector::~Ali3Vector()
+{
+// Destructor to delete dynamically allocated memory
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali3Vector::SetVector(Double_t* v,TString f)
+{
+// Store vector according to reference frame f
+ Double_t pi=acos(-1.);
+ Int_t frame=0;
+ if (f == "car") frame=1;
+ if (f == "sph") frame=2;
+ if (f == "cyl") frame=3;
+
+ Double_t x,y,z,rho,phi;
+
+ switch (frame)
+ {
+ case 1: // Cartesian coordinates
+ x=v[0];
+ y=v[1];
+ z=v[2];
+ fV=sqrt(x*x+y*y+z*z);
+ fTheta=0;
+ if (fV && fabs(z/fV)<=1.)
+ {
+ fTheta=acos(z/fV);
+ }
+ else
+ {
+ if (z<0.) fTheta=pi;
+ }
+ if (fTheta<0.) fTheta+=2.*pi;
+ fPhi=0;
+ if (x || y) fPhi=atan2(y,x);
+ if (fPhi<0.) fPhi+=2.*pi;
+ break;
+
+ case 2: // Spherical coordinates
+ fV=v[0];
+ fTheta=v[1];
+ fPhi=v[2];
+ break;
+
+ case 3: // Cylindrical coordinates
+ rho=v[0];
+ phi=v[1];
+ z=v[2];
+ fV=sqrt(rho*rho+z*z);
+ fPhi=phi;
+ if (fPhi<0.) fPhi+=2.*pi;
+ fTheta=0;
+ if (fV && fabs(z/fV)<=1.)
+ {
+ fTheta=acos(z/fV);
+ }
+ else
+ {
+ if (z<0.) fTheta=pi;
+ }
+ if (fTheta<0.) fTheta+=2.*pi;
+ break;
+
+ default: // Unsupported reference frame
+ cout << "*Ali3Vector::SetVector* Unsupported frame : " << f << endl
+ << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
+ fV=0;
+ fTheta=0;
+ fPhi=0;
+ break;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali3Vector::GetVector(Double_t* v,TString f)
+{
+// Provide vector according to reference frame f
+ Int_t frame=0;
+ if (f == "car") frame=1;
+ if (f == "sph") frame=2;
+ if (f == "cyl") frame=3;
+
+ switch (frame)
+ {
+ case 1: // Cartesian coordinates
+ v[0]=fV*sin(fTheta)*cos(fPhi);
+ v[1]=fV*sin(fTheta)*sin(fPhi);
+ v[2]=fV*cos(fTheta);
+ break;
+
+ case 2: // Spherical coordinates
+ v[0]=fV;
+ v[1]=fTheta;
+ v[2]=fPhi;
+ break;
+
+ case 3: // Cylindrical coordinates
+ v[0]=fV*sin(fTheta);
+ v[1]=fPhi;
+ v[2]=fV*cos(fTheta);
+ break;
+
+ default: // Unsupported reference frame
+ cout << "*Ali3Vector::GetVector* Unsupported frame : " << f << endl
+ << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
+ for (Int_t i=0; i<3; i++)
+ {
+ v[i]=0;
+ }
+ break;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali3Vector::SetVector(Float_t* v,TString f)
+{
+// Store vector according to reference frame f
+ Double_t vec[3];
+ for (Int_t i=0; i<3; i++)
+ {
+ vec[i]=v[i];
+ }
+ SetVector(vec,f);
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali3Vector::GetVector(Float_t* v,TString f)
+{
+// Provide vector according to reference frame f
+ Double_t vec[3];
+ GetVector(vec,f);
+ for (Int_t i=0; i<3; i++)
+ {
+ v[i]=vec[i];
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali3Vector::Info(TString f)
+{
+// Print vector components according to reference frame f
+ if (f=="car" || f=="sph" || f=="cyl")
+ {
+ Double_t vec[3];
+ GetVector(vec,f);
+ cout << " Vector in " << f << " coordinates : "
+ << vec[0] << " " << vec[1] << " " << vec[2] << endl;
+ }
+ else
+ {
+ cout << " *Ali3Vector::Info* Unsupported frame : " << f << endl
+ << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t Ali3Vector::GetNorm()
+{
+ return fV;
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t Ali3Vector::Dot(Ali3Vector& q)
+{
+// Provide the dot product of the current vector with vector q
+ Double_t a[3],b[3];
+ Double_t dotpro=0;
+
+ GetVector(a,"car");
+ q.GetVector(b,"car");
+ for (Int_t i=0; i<3; i++)
+ {
+ dotpro+=a[i]*b[i];
+ }
+
+ return dotpro;
+}
+///////////////////////////////////////////////////////////////////////////
+Ali3Vector Ali3Vector::Cross(Ali3Vector& q)
+{
+// Provide the cross product of the current vector with vector q
+ Double_t a[3],b[3],c[3];
+
+ GetVector(a,"car");
+ q.GetVector(b,"car");
+
+ c[0]=a[1]*b[2]-a[2]*b[1];
+ c[1]=a[2]*b[0]-a[0]*b[2];
+ c[2]=a[0]*b[1]-a[1]*b[0];
+
+ Ali3Vector v;
+ v.SetVector(c,"car");
+
+ return v;
+}
+///////////////////////////////////////////////////////////////////////////
+Ali3Vector Ali3Vector::operator+(Ali3Vector& q)
+{
+// Add vector q to the current vector
+ Double_t a[3],b[3];
+
+ GetVector(a,"car");
+ q.GetVector(b,"car");
+
+ for (Int_t i=0; i<3; i++)
+ {
+ a[i]+=b[i];
+ }
+
+ Ali3Vector v;
+ v.SetVector(a,"car");
+
+ return v;
+}
+///////////////////////////////////////////////////////////////////////////
+Ali3Vector Ali3Vector::operator-(Ali3Vector& q)
+{
+// Subtract vector q from the current vector
+ Double_t a[3],b[3];
+
+ GetVector(a,"car");
+ q.GetVector(b,"car");
+
+ for (Int_t i=0; i<3; i++)
+ {
+ a[i]-=b[i];
+ }
+
+ Ali3Vector v;
+ v.SetVector(a,"car");
+
+ return v;
+}
+///////////////////////////////////////////////////////////////////////////
+Ali3Vector Ali3Vector::operator*(Double_t s)
+{
+// Multiply the current vector with a scalar s
+ Double_t a[3];
+
+ GetVector(a,"car");
+
+ for (Int_t i=0; i<3; i++)
+ {
+ a[i]*=s;
+ }
+
+ Ali3Vector v;
+ v.SetVector(a,"car");
+
+ return v;
+}
+///////////////////////////////////////////////////////////////////////////
+Ali3Vector Ali3Vector::operator/(Double_t s)
+{
+// Divide the current vector by a scalar s
+
+ if (fabs(s)<1.e-20) // Protect against division by 0
+ {
+ cout << " *Ali3Vector::/* Division by 0 detected. No action taken." << endl;
+ return *this;
+ }
+ else
+ {
+ Double_t a[3];
+
+ GetVector(a,"car");
+
+ for (Int_t i=0; i<3; i++)
+ {
+ a[i]/=s;
+ }
+
+ Ali3Vector v;
+ v.SetVector(a,"car");
+
+ return v;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Ali3Vector& Ali3Vector::operator+=(Ali3Vector& q)
+{
+// Add vector q to the current vector
+ Double_t a[3],b[3];
+
+ GetVector(a,"car");
+ q.GetVector(b,"car");
+
+ for (Int_t i=0; i<3; i++)
+ {
+ a[i]+=b[i];
+ }
+
+ SetVector(a,"car");
+
+ return *this;
+}
+///////////////////////////////////////////////////////////////////////////
+Ali3Vector& Ali3Vector::operator-=(Ali3Vector& q)
+{
+// Subtract vector q from the current vector
+ Double_t a[3],b[3];
+
+ GetVector(a,"car");
+ q.GetVector(b,"car");
+
+ for (Int_t i=0; i<3; i++)
+ {
+ a[i]-=b[i];
+ }
+
+ SetVector(a,"car");
+
+ return *this;
+}
+///////////////////////////////////////////////////////////////////////////
+Ali3Vector& Ali3Vector::operator*=(Double_t s)
+{
+// Multiply the current vector with a scalar s
+ Double_t a[3];
+
+ GetVector(a,"car");
+
+ for (Int_t i=0; i<3; i++)
+ {
+ a[i]*=s;
+ }
+
+ SetVector(a,"car");
+
+ return *this;
+}
+///////////////////////////////////////////////////////////////////////////
+Ali3Vector& Ali3Vector::operator/=(Double_t s)
+{
+// Divide the current vector by a scalar s
+
+ if (fabs(s)<1.e-20) // Protect against division by 0
+ {
+ cout << " *Ali3Vector::/=* Division by 0 detected. No action taken." << endl;
+ return *this;
+ }
+ else
+ {
+ Double_t a[3];
+
+ GetVector(a,"car");
+
+ for (Int_t i=0; i<3; i++)
+ {
+ a[i]/=s;
+ }
+
+ SetVector(a,"car");
+
+ return *this;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
--- /dev/null
+#ifndef ALI3VECTOR_H
+#define ALI3VECTOR_H
+///////////////////////////////////////////////////////////////////////////
+// Class Ali3Vector
+// Handling of 3-vectors in various reference frames.
+//
+// This class is meant to serve as a base class for ALICE objects
+// that have 3-dimensional vector characteristics.
+//
+// Note :
+// ------
+// Vectors (v) and reference frames (f) are specified via
+// SetVector(Float_t* v,TString f) under the following conventions :
+//
+// f="car" ==> v in Cartesian coordinates (x,y,z)
+// f="sph" ==> v in Spherical coordinates (r,theta,phi)
+// f="cyl" ==> v in Cylindrical coordinates (rho,phi,z)
+//
+// All angles are in radians.
+//
+// Example :
+// ---------
+//
+// Ali3Vector a;
+// Float_t v[3]={-1,25,7};
+// a.SetVector(v,"car");
+// a.Info();
+//
+// Float_t vec[3];
+// a.GetVector(vec,"sph");
+//
+// Ali3Vector b;
+// Float_t v2[3]={6,-18,33};
+// b.SetVector(v2,"car");
+//
+// Float_t dotpro=a.Dot(b);
+//
+// Ali3Vector c=a.Cross(b);
+// c.Info("sph");
+// c.GetVector(vec,"cyl");
+// Float_t norm=c.GetNorm();
+// c=a+b;
+// c=a-b;
+// c=a*5;
+//
+//--- NvE 30-mar-1999 UU-SAP Utrecht
+///////////////////////////////////////////////////////////////////////////
+
+#include <iostream.h>
+#include <math.h>
+
+#include "TObject.h"
+#include "TString.h"
+
+class Ali3Vector
+{
+ public:
+ Ali3Vector(); // Default constructor
+ virtual ~Ali3Vector(); // Destructor
+ virtual void SetVector(Double_t* v,TString f); // Store vector v in frame f
+ virtual void GetVector(Double_t* v,TString f); // Provide vector v in frame f
+ virtual void SetVector(Float_t* v,TString f); // Store vector v in frame f
+ virtual void GetVector(Float_t* v,TString f); // Provide vector v in frame f
+ virtual void Info(TString f="car"); // Print vector components in frame f
+ Double_t GetNorm(); // Provide norm of the vector
+ Double_t Dot(Ali3Vector& q); // Provide dot product with q
+ Ali3Vector Cross(Ali3Vector& q); // Provide cross product with q
+ Ali3Vector operator+(Ali3Vector& q); // Add vector q
+ Ali3Vector operator-(Ali3Vector& q); // Subtract vector q
+ Ali3Vector operator*(Double_t s); // Multiply vector with scalar s
+ Ali3Vector operator/(Double_t s); // Divide vector by scalar s
+ Ali3Vector& operator+=(Ali3Vector& q); // Add vector q
+ Ali3Vector& operator-=(Ali3Vector& q); // Subtract vector q
+ Ali3Vector& operator*=(Double_t s); // Multiply with scalar s
+ Ali3Vector& operator/=(Double_t s); // Divide by scalar s
+
+ protected:
+ Double_t fV,fTheta,fPhi; // Vector in spherical coordinates
+
+ ClassDef(Ali3Vector,1) // Class definition to enable ROOT I/O
+};
+#endif
--- /dev/null
+#include "Ali4Vector.h"
+
+ClassImp(Ali4Vector) // Class implementation to enable ROOT I/O
+
+Ali4Vector::Ali4Vector()
+{
+// Creation of a contravariant 4-vector and initialisation of parameters
+ fV0=0;
+ Double_t a[3]={0,0,0};
+ fV.SetVector(a,"sph");
+}
+///////////////////////////////////////////////////////////////////////////
+Ali4Vector::~Ali4Vector()
+{
+// Destructor to delete dynamically allocated memory
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali4Vector::SetVector(Double_t v0,Ali3Vector v)
+{
+// Store contravariant vector
+ fV0=v0;
+ fV=v;
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali4Vector::SetVector(Double_t* v,TString f)
+{
+// Store vector according to reference frame f
+ fV0=v[0];
+ Double_t a[3];
+ for (Int_t i=0; i<3; i++)
+ {
+ a[i]=v[i+1];
+ }
+ fV.SetVector(a,f);
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali4Vector::GetVector(Double_t* v,TString f)
+{
+// Provide vector according to reference frame f
+ v[0]=fV0;
+ Double_t a[3];
+ fV.GetVector(a,f);
+ for (Int_t i=0; i<3; i++)
+ {
+ v[i+1]=a[i];
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali4Vector::SetVector(Float_t* v,TString f)
+{
+// Store vector according to reference frame f
+ Double_t vec[4];
+ for (Int_t i=0; i<4; i++)
+ {
+ vec[i]=v[i];
+ }
+ SetVector(vec,f);
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali4Vector::GetVector(Float_t* v,TString f)
+{
+// Provide vector according to reference frame f
+ Double_t vec[4];
+ GetVector(vec,f);
+ for (Int_t i=0; i<4; i++)
+ {
+ v[i]=vec[i];
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t Ali4Vector::GetScalar()
+{
+// Provide the scalar part
+ return fV0;
+}
+///////////////////////////////////////////////////////////////////////////
+Ali3Vector Ali4Vector::Get3Vector()
+{
+// Provide the 3-vector part
+ return fV;
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali4Vector::Info(TString f)
+{
+// Print contravariant vector components according to reference frame f
+ if (f=="car" || f=="sph" || f=="cyl")
+ {
+ Double_t vec[4];
+ GetVector(vec,f);
+ cout << " Contravariant vector in " << f << " coordinates : "
+ << vec[0] << " " << vec[1] << " " << vec[2] << " " << vec[3] << endl;
+ }
+ else
+ {
+ cout << " *Ali4Vector::Info* Unsupported frame : " << f << endl
+ << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t Ali4Vector::Dot(Ali4Vector& q)
+{
+// Provide the dot product of the current vector with vector q
+ Double_t a[4],b[4];
+ Double_t dotpro;
+
+ GetVector(a,"car");
+ q.GetVector(b,"car");
+
+ dotpro=a[0]*b[0];
+ for (Int_t i=1; i<4; i++)
+ {
+ dotpro-=a[i]*b[i];
+ }
+
+ return dotpro;
+}
+///////////////////////////////////////////////////////////////////////////
+Ali4Vector Ali4Vector::operator+(Ali4Vector& q)
+{
+// Add vector q to the current vector
+ Double_t a[4],b[4];
+
+ GetVector(a,"car");
+ q.GetVector(b,"car");
+
+ for (Int_t i=0; i<4; i++)
+ {
+ a[i]+=b[i];
+ }
+
+ Ali4Vector v;
+ v.SetVector(a,"car");
+
+ return v;
+}
+///////////////////////////////////////////////////////////////////////////
+Ali4Vector Ali4Vector::operator-(Ali4Vector& q)
+{
+// Subtract vector q from the current vector
+ Double_t a[4],b[4];
+
+ GetVector(a,"car");
+ q.GetVector(b,"car");
+
+ for (Int_t i=0; i<4; i++)
+ {
+ a[i]-=b[i];
+ }
+
+ Ali4Vector v;
+ v.SetVector(a,"car");
+
+ return v;
+}
+///////////////////////////////////////////////////////////////////////////
+Ali4Vector Ali4Vector::operator*(Double_t s)
+{
+// Multiply the current vector with a scalar s
+ Double_t a[4];
+
+ GetVector(a,"car");
+
+ for (Int_t i=0; i<4; i++)
+ {
+ a[i]*=s;
+ }
+
+ Ali4Vector v;
+ v.SetVector(a,"car");
+
+ return v;
+}
+///////////////////////////////////////////////////////////////////////////
+Ali4Vector Ali4Vector::operator/(Double_t s)
+{
+// Divide the current vector by a scalar s
+
+ if (fabs(s)<1.e-20) // Protect against division by 0
+ {
+ cout << " *Ali4Vector::/* Division by 0 detected. No action taken." << endl;
+ return *this;
+ }
+ else
+ {
+ Double_t a[4];
+
+ GetVector(a,"car");
+
+ for (Int_t i=0; i<4; i++)
+ {
+ a[i]/=s;
+ }
+
+ Ali4Vector v;
+ v.SetVector(a,"car");
+
+ return v;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Ali4Vector& Ali4Vector::operator+=(Ali4Vector& q)
+{
+// Add vector q to the current vector
+ Double_t a[4],b[4];
+
+ GetVector(a,"car");
+ q.GetVector(b,"car");
+
+ for (Int_t i=0; i<4; i++)
+ {
+ a[i]+=b[i];
+ }
+
+ SetVector(a,"car");
+
+ return *this;
+}
+///////////////////////////////////////////////////////////////////////////
+Ali4Vector& Ali4Vector::operator-=(Ali4Vector& q)
+{
+// Subtract vector q from the current vector
+ Double_t a[4],b[4];
+
+ GetVector(a,"car");
+ q.GetVector(b,"car");
+
+ for (Int_t i=0; i<4; i++)
+ {
+ a[i]-=b[i];
+ }
+
+ SetVector(a,"car");
+
+ return *this;
+}
+///////////////////////////////////////////////////////////////////////////
+Ali4Vector& Ali4Vector::operator*=(Double_t s)
+{
+// Multiply the current vector with a scalar s
+ Double_t a[4];
+
+ GetVector(a,"car");
+
+ for (Int_t i=0; i<4; i++)
+ {
+ a[i]*=s;
+ }
+
+ SetVector(a,"car");
+
+ return *this;
+}
+///////////////////////////////////////////////////////////////////////////
+Ali4Vector& Ali4Vector::operator/=(Double_t s)
+{
+// Divide the current vector by a scalar s
+
+ if (fabs(s)<1.e-20) // Protect against division by 0
+ {
+ cout << " *Ali4Vector::/=* Division by 0 detected. No action taken." << endl;
+ return *this;
+ }
+ else
+ {
+ Double_t a[4];
+
+ GetVector(a,"car");
+
+ for (Int_t i=0; i<4; i++)
+ {
+ a[i]/=s;
+ }
+
+ SetVector(a,"car");
+
+ return *this;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
--- /dev/null
+#ifndef ALI4VECTOR_H
+#define ALI4VECTOR_H
+///////////////////////////////////////////////////////////////////////////
+// Class Ali4Vector
+// Handling of Lorentz 4-vectors in various reference frames.
+//
+// This class is meant to serve as a base class for ALICE objects
+// that have Lorentz 4-vector characteristics.
+//
+// All 4-vectors are treated in the contravariant form and the convention
+// for the metric and the 4-vector components is according to the one
+// used in the book "Classical Electrodynamics" by J.D. Jackson.
+//
+// The dotproduct is defined such that p.Dot(p) yields the Lorentz invariant
+// scalar of the 4-vector p (i.e. m**2 in case p is a 4-momentum).
+//
+// Note :
+// ------
+// Vectors (v) and reference frames (f) are specified via
+// SetVector(Float_t* v,TString f) under the following conventions :
+//
+// f="car" ==> 3-vector part of v in Cartesian coordinates (x,y,z)
+// f="sph" ==> 3-vector part of v in Spherical coordinates (r,theta,phi)
+// f="cyl" ==> 3-vector part of v in Cylindrical coordinates (rho,phi,z)
+//
+// All angles are in radians.
+//
+// Example :
+// ---------
+//
+// Ali4Vector a;
+//
+// Float_t v[4]={25,-1,3,7};
+// a.SetVector(v,"car");
+//
+// Float_t vec[4];
+// a.GetVector(vec,"sph");
+//
+// Ali4Vector b;
+// Float_t v2[4]={33,6,-18,2};
+// b.SetVector(v2,"car");
+//
+// Float_t dotpro=a.Dot(b);
+//
+// Float_t x0=16;
+// Ali3Vector x;
+// Float_t vec2[3]={1,2,3};
+// x.SetVector(vec2,"car");
+//
+// Ali4Vector c;
+// c.SetVector(x0,x);
+// c.GetVector(vec,"car");
+// c.Info("cyl");
+// c=a+b;
+// c=a-b;
+// c=a*5;
+//
+//--- NvE 01-apr-1999 UU-SAP Utrecht
+///////////////////////////////////////////////////////////////////////////
+
+#include <iostream.h>
+#include <math.h>
+
+#include "Ali3Vector.h"
+
+class Ali4Vector
+{
+ public:
+ Ali4Vector(); // Default constructor for contravariant vector
+ virtual ~Ali4Vector(); // Destructor
+ virtual void SetVector(Double_t v0,Ali3Vector v); // Store contravariant vector
+ virtual void SetVector(Double_t* v,TString f); // Store contravariant vector v^i in frame f
+ virtual void GetVector(Double_t* v,TString f); // Provide contravariant vector v^i in frame f
+ virtual void SetVector(Float_t* v,TString f); // Store contravariant vector v^i in frame f
+ virtual void GetVector(Float_t* v,TString f); // Provide contravariant vector v^i in frame f
+ Double_t GetScalar(); // Provide the scalar part of v
+ Ali3Vector Get3Vector(); // Provide the 3-vector part of v
+ virtual void Info(TString f="car"); // Print contravariant components in frame f
+ Double_t Dot(Ali4Vector& q); // Provide dot product v^i*q_i
+ Ali4Vector operator+(Ali4Vector& q); // Add contravariant vector q
+ Ali4Vector operator-(Ali4Vector& q); // Subtract contravariant vector q
+ Ali4Vector operator*(Double_t s); // Multiply contravariant vector with scalar s
+ Ali4Vector operator/(Double_t s); // Divide contravariant vector by scalar s
+ Ali4Vector& operator+=(Ali4Vector& q); // Add contravariant vector q
+ Ali4Vector& operator-=(Ali4Vector& q); // Subtract contravariant vector q
+ Ali4Vector& operator*=(Double_t s); // Multiply with scalar s
+ Ali4Vector& operator/=(Double_t s); // Divide by scalar s
+
+ protected:
+ Double_t fV0; // The scalar part
+ Ali3Vector fV; // The 3-vector part
+
+ ClassDef(Ali4Vector,1) // Class definition to enable ROOT I/O
+};
+#endif
--- /dev/null
+#include "AliBoost.h"
+
+ClassImp(AliBoost) // Class implementation to enable ROOT I/O
+
+AliBoost::AliBoost()
+{
+// Creation of a Lorentz boost object and initialisation of parameters
+ fGamma=1;
+ fBeta2=0;
+ Double_t a[3]={0,0,0};
+ fBeta.SetVector(a,"sph");
+}
+///////////////////////////////////////////////////////////////////////////
+AliBoost::~AliBoost()
+{
+// Default destructor
+}
+///////////////////////////////////////////////////////////////////////////
+void AliBoost::SetBeta(Ali3Vector b)
+{
+// Setting of boost parameters on basis of beta 3-vector
+ fBeta2=b.Dot(b);
+ fBeta=b;
+
+ if (fBeta2 > 1.)
+ {
+ cout << " *AliBoost::SetBeta* beta > 1." << endl;
+ }
+ Double_t test=1.-fBeta2;
+ fGamma=0;
+ if (test > 0.) fGamma=sqrt(1./test);
+}
+///////////////////////////////////////////////////////////////////////////
+void AliBoost::SetGamma(Double_t g,Ali3Vector v)
+{
+// Setting of boost parameters on basis of gamma and direction 3-vector
+ if (g >= 1.)
+ {
+ fGamma=g;
+ fBeta2=1.-(1./(fGamma*fGamma));
+ fBeta=v*sqrt(fBeta2);
+ }
+ else
+ {
+ cout << " *AliBoost::SetGamma* Invalid input gamma = " << g << endl;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliBoost::Set4Momentum(Ali4Vector& p)
+{
+// Setting of boost parameters on basis of momentum 4-vector data
+ Double_t E=p.GetScalar();
+ if (E <= 0.)
+ {
+ cout << " *AliBoost::Set4Momentum* Unphysical situation." << endl;
+ p.Info();
+ }
+ else
+ {
+ Ali3Vector b=p.Get3Vector();
+ b=b/E;
+ SetBeta(b);
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Ali3Vector AliBoost::GetBetaVector()
+{
+// Provide the the beta 3-vector
+ return fBeta;
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t AliBoost::GetBeta()
+{
+// Provide the norm of the beta 3-vector
+ return sqrt(fBeta2);
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t AliBoost::GetGamma()
+{
+// Provide the gamma factor
+ return fGamma;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliBoost::Info(TString f)
+{
+// Printing of the boost parameter info in coordinate frame f
+
+ cout << " *AliBoost::Info* beta = " << sqrt(fBeta2) << " gamma = " << fGamma << endl
+ << " Beta";
+ fBeta.Info(f);
+}
+///////////////////////////////////////////////////////////////////////////
+Ali4Vector AliBoost::Boost(Ali4Vector& v)
+{
+// Perform the Lorentz boost on the 4-vector v
+ if (fBeta2 > 1.e-20)
+ {
+ Double_t E=v.GetScalar();
+ Ali3Vector p=v.Get3Vector();
+ Double_t pdotbeta=p.Dot(fBeta);
+
+ Double_t Eprim;
+ Eprim=fGamma*(E-pdotbeta);
+
+ Ali3Vector term1,term2,pprim;
+ term1=fBeta*((fGamma-1.)*pdotbeta/fBeta2);
+ term2=fBeta*(fGamma*E);
+ pprim=p+term1-term2;
+
+ Ali4Vector w;
+ w.SetVector(Eprim,pprim);
+
+ return w;
+ }
+ else
+ {
+ return v;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Ali4Vector AliBoost::Inverse(Ali4Vector& vprim)
+{
+// Perform the inverse Lorentz boost on the 4-vector vprim
+ if (fBeta2 > 1.e-20)
+ {
+ Double_t Eprim=vprim.GetScalar();
+ Ali3Vector pprim=vprim.Get3Vector();
+ Double_t pprimdotbeta=pprim.Dot(fBeta);
+
+ Double_t E;
+ E=fGamma*(Eprim+pprimdotbeta);
+
+ Ali3Vector term1,term2,p;
+ term1=fBeta*((fGamma-1.)*pprimdotbeta/fBeta2);
+ term2=fBeta*(fGamma*Eprim);
+ p=pprim+term1+term2;
+
+ Ali4Vector w;
+ w.SetVector(E,p);
+
+ return w;
+ }
+ else
+ {
+ return vprim;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
--- /dev/null
+#ifndef ALIBOOST_H
+#define ALIBOOST_H
+///////////////////////////////////////////////////////////////////////////
+// Class AliBoost
+// Perform various Lorentz transformations.
+//
+// Example :
+// =========
+//
+// Float_t a[3]={0.1,0.2,0.3};
+// Ali3Vector beta;
+// beta.SetVector(a,"car");
+//
+// AliBoost b1;
+// b1.SetBeta(beta);
+// cout << " Boost b1 :" << endl;
+// b1.Info();
+//
+// Float_t b[4]={14,1,2,3};
+// Ali4Vector p;
+// p.SetVector(b,"car");
+// Ali4Vector pprim=b1.Boost(p);
+// cout << " Boost b1 result p,pprim :" << endl;
+// p.Info();
+// pprim.Info();
+//
+// p=b1.Inverse(pprim);
+// cout << " Inverse b1 result pprim,p :" << endl;
+// pprim.Info();
+// p.Info();
+//
+// Float_t c[4]={5,0,0,4};
+// Ali4Vector q;
+// q.SetVector(c,"car");
+//
+// AliBoost b2;
+// b2.Set4Momentum(q);
+// cout << " Lorbo b2 : " << endl;
+// b2.Info("sph");
+//
+//--- NvE 14-may-1996 UU-SAP Utrecht
+//--- Modified : NvE 01-apr-1999 UU-SAP Utrecht using Ali3Vector/Ali4Vector
+///////////////////////////////////////////////////////////////////////////
+
+#include <iostream.h>
+#include <math.h>
+
+#include "TObject.h"
+
+#include "Ali4Vector.h"
+
+class AliBoost : public TObject
+{
+ public:
+ AliBoost(); // Default constructor
+ ~AliBoost(); // Default destructor
+ void SetBeta(Ali3Vector b); // Set boost parameters by beta 3-vector
+ void SetGamma(Double_t g,Ali3Vector v); // Set boost parameters by gamma and direction 3-vector
+ void Set4Momentum(Ali4Vector& p); // Set boost parameters by 4-momentum
+ Ali3Vector GetBetaVector(); // Provide the beta 3-vector
+ Double_t GetBeta(); // Provide norm of beta 3-vector
+ Double_t GetGamma(); // Provide gamma value
+ void Info(TString f="car"); // Print boost parameter info in coord. frame f
+ Ali4Vector Boost(Ali4Vector& v); // Perform Lorentz boost on 4-vector v
+ Ali4Vector Inverse(Ali4Vector& v); // Perform inverse Lorentz boost on 4-vector v
+
+ protected:
+ Ali3Vector fBeta; // The beta 3-vector
+ Double_t fBeta2; // beta**2
+ Double_t fGamma; // The gamma factor
+
+ ClassDef(AliBoost,1) // Class definition to enable ROOT I/O
+};
+#endif
--- /dev/null
+#include "AliCalcluster.h"
+
+ClassImp(AliCalcluster) // Class implementation to enable ROOT I/O
+
+AliCalcluster::AliCalcluster()
+{
+// Default constructer, all data is set to 0
+ fCenter=0;
+ fSig=0.;
+ fNmods=0;
+ fSig11=0.;
+ fSig33=0.;
+ fSig55=0.;
+ fRowdisp=0.;
+ fColdisp=0.;
+ fNvetos=0;
+ fVetos=0;
+}
+///////////////////////////////////////////////////////////////////////////
+AliCalcluster::~AliCalcluster()
+{
+// Destructor to delete dynamically allocated memory
+ if (fVetos)
+ {
+ fVetos->Delete();
+ delete fVetos;
+ fVetos=0;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+AliCalcluster::AliCalcluster(AliCalmodule& m)
+{
+// Cluster constructor with module m as center.
+// Module data is only entered for a module which contains
+// a signal and has not been used in a cluster yet.
+// Modules situated at a detector edge are not allowed to start a cluster.
+
+ Float_t pos[3]={0,0,0};
+
+ if ((m.GetClusteredSignal() > 0.) && (m.GetEdgeValue() == 0))
+ {
+ fCenter=&m;
+ m.GetPosition(pos,"sph");
+ SetPosition(pos,"sph");
+ fSig=m.GetClusteredSignal();
+ fNmods=1;
+ fSig11=m.GetClusteredSignal();
+ fSig33=m.GetClusteredSignal();
+ fSig55=m.GetClusteredSignal();
+ fRowdisp=0.;
+ fColdisp=0.;
+ m.SetClusteredSignal(0.); // mark module as used in cluster
+ fNvetos=0;
+ fVetos=0;
+ }
+ else
+ {
+ fCenter=0;
+ SetPosition(pos,"sph");
+ fSig=0.;
+ fNmods=0;
+ fSig11=0.;
+ fSig33=0.;
+ fSig55=0.;
+ fRowdisp=0.;
+ fColdisp=0.;
+ fNvetos=0;
+ fVetos=0;
+}
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliCalcluster::GetRow()
+{
+// Provide the row number of the cluster center
+ if (fCenter)
+ {
+ return fCenter->GetRow();
+ }
+ else
+ {
+ return 0;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliCalcluster::GetColumn()
+{
+// Provide the column number of the cluster center
+ if (fCenter)
+ {
+ return fCenter->GetColumn();
+ }
+ else
+ {
+ return 0;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliCalcluster::GetSignal(Int_t n)
+{
+// Provide the total signal of a module matrix of n modules around
+// the cluster center.
+// Example : n=9 --> total signal in the 3x3 matrix
+// 1 --> signal of central module
+// Note : n=0 provides the total cluster signal (Default)
+
+ Float_t signal=-1;
+
+ if (n == 0) signal=fSig;
+ if (n == 1) signal=fSig11;
+ if (n == 9) signal=fSig33;
+ if (n == 25) signal=fSig55;
+
+ if (signal > 0.)
+ {
+ return signal;
+ }
+ else
+ {
+ cout << " *AliCalcluster::GetSignal* Invalid argument n = " << n << endl;
+ return 0;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliCalcluster::GetNmodules()
+{
+// Provide the number of modules in the cluster
+ return fNmods;
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliCalcluster::GetRowDispersion()
+{
+// Provide the normalised row dispersion of the cluster
+ if (fSig > 0.)
+ {
+ return fRowdisp/fSig;
+ }
+ else
+ {
+ return 0.;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliCalcluster::GetColumnDispersion()
+{
+// Provide the normalised column dispersion of the cluster
+ if (fSig > 0.)
+ {
+ return fColdisp/fSig;
+ }
+ else
+ {
+ return 0.;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalcluster::Start(AliCalmodule& m)
+{
+// Reset cluster data and start with module m
+// A module can only start a cluster when it contains
+// a signal, has not been used in a cluster yet and is not
+// situated at a detector edge
+
+ Float_t pos[3]={0,0,0};
+
+ if ((m.GetClusteredSignal() > 0.) && (m.GetEdgeValue() == 0))
+ {
+ fCenter=&m;
+ m.GetPosition(pos,"sph");
+ SetPosition(pos,"sph");
+ fSig=m.GetSignal();
+ fNmods=1;
+ fSig11=m.GetSignal();
+ fSig33=m.GetSignal();
+ fSig55=m.GetSignal();
+ fRowdisp=0.;
+ fColdisp=0.;
+ m.SetClusteredSignal(0.); // mark module as used in cluster
+ }
+ else
+ {
+ fCenter=0;
+ SetPosition(pos,"sph");
+ fSig=0.;
+ fNmods=0;
+ fSig11=0.;
+ fSig33=0.;
+ fSig55=0.;
+ fRowdisp=0.;
+ fColdisp=0.;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalcluster::Add(AliCalmodule& m)
+{
+// Add module data to the cluster
+
+ Float_t signal=m.GetClusteredSignal();
+
+ if (signal > 0.) // only add unused modules
+ {
+ fSig+=signal;
+ fNmods+=1;
+ Int_t drow=int(fabs(double(GetRow()-m.GetRow()))); // row distance to center
+ Int_t dcol=int(fabs(double(GetColumn()-m.GetColumn()))); // column distance to center
+ if ((drow < 2) && (dcol < 2)) fSig33+=signal;
+ if ((drow < 3) && (dcol < 3)) fSig55+=signal;
+ fRowdisp+=signal*float(drow*drow);
+ fColdisp+=signal*float(dcol*dcol);
+ m.SetClusteredSignal(0.); // mark module as used in cluster
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalcluster::AddVetoSignal(Float_t* r,TString f,Float_t s)
+{
+// Associate an (extrapolated) AliSignal at location r as veto to the cluster
+// Note : The default signal value (s) is 0
+ if (!fVetos)
+ {
+ fNvetos=0;
+ fVetos=new TObjArray();
+ }
+
+ fVetos->Add(new AliSignal);
+ fNvetos++;
+
+ ((AliSignal*)fVetos->At(fNvetos-1))->SetPosition(r,f);
+ ((AliSignal*)fVetos->At(fNvetos-1))->SetSignal(s);
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliCalcluster::GetNvetos()
+{
+// Provide the number of veto signals associated to the cluster
+ return fNvetos;
+}
+///////////////////////////////////////////////////////////////////////////
+AliSignal* AliCalcluster::GetVetoSignal(Int_t i)
+{
+// Provide access to the i-th veto signal of this cluster
+// Note : The first hit corresponds to i=1
+
+ if (i>0 && i<=fNvetos)
+ {
+ return (AliSignal*)fVetos->At(i-1);
+ }
+ else
+ {
+ cout << " *AliCalcluster::GetVetoSignal* Signal number " << i
+ << " out of range." << endl;
+ cout << " --- First signal (if any) returned." << endl;
+ return (AliSignal*)fVetos->At(0);
+ }
+}
+///////////////////////////////////////////////////////////////////////////
--- /dev/null
+#ifndef ALICALCLUSTER_H
+#define ALICALCLUSTER_H
+///////////////////////////////////////////////////////////////////////////
+// Class AliCalcluster
+// Description of a cluster of calorimeter modules.
+// A matrix geometry is assumed in which a cluster center
+// is identified by (row,col) and contains sig as signal
+// being the signal of the complete cluster.
+// Some info about cluster topology is provided in order
+// to enable EM or hadronic cluster identification
+//
+//--- NvE 13-jun-1997 UU-SAP Utrecht
+///////////////////////////////////////////////////////////////////////////
+
+#include <iostream.h>
+#include <math.h>
+
+#include "TObject.h"
+#include "TObjArray.h"
+#include "TString.h"
+
+#include "AliCalmodule.h"
+
+class AliCalcluster : public TObject,public AliPosition
+{
+ public:
+ AliCalcluster(); // Default constructor, all data initialised to 0
+ ~AliCalcluster(); // Default destructor
+ AliCalcluster(AliCalmodule& m); // Create new cluster starting at module m
+ Int_t GetRow(); // Return row number of cluster center
+ Int_t GetColumn(); // Return column number of cluster center
+ Float_t GetSignal(Int_t n=0); // Return signal of nxn matrix around the center
+ Int_t GetNmodules(); // Return number of modules in cluster
+ Float_t GetRowDispersion(); // Return normalised row dispersion of cluster
+ Float_t GetColumnDispersion(); // Return normalised column dispersion of cluster
+ void Start(AliCalmodule& m); // Reset cluster data to start with module m
+ void Add(AliCalmodule& m); // Add module data to cluster
+ void AddVetoSignal(Float_t* r,TString f,Float_t s=0); // Associate (extrapolated) signal
+ AliSignal* GetVetoSignal(Int_t j); // Access to veto signal number j
+ Int_t GetNvetos(); // Provide the number of veto signals
+
+ protected:
+ AliCalmodule* fCenter; // Pointer to the central module of the cluster
+ Float_t fSig; // The total signal value of the cluster
+ Int_t fNmods; // The number of modules in the cluster
+ Float_t fSig11; // Cluster signal of the central module
+ Float_t fSig33; // Cluster signal in 3x3 matrix around the center
+ Float_t fSig55; // Cluster signal in 5x5 matrix around the center
+ Float_t fRowdisp; // Row dispersion of cluster (not normalised)
+ Float_t fColdisp; // Column dispersion of cluster (not normalised)
+ Int_t fNvetos; // The number of associated veto signals
+ TObjArray* fVetos; // The array of associated veto signals
+
+ ClassDef(AliCalcluster,1) // Class definition to enable ROOT I/O
+};
+#endif
--- /dev/null
+#include "AliCalmodule.h"
+
+ClassImp(AliCalmodule) // Class implementation to enable ROOT I/O
+
+AliCalmodule::AliCalmodule()
+{
+// Default constructor, all module data is set to 0
+ fRow=0;
+ fCol=0;
+ fSigc=0;
+ fEdge=0;
+ fDead=0;
+ fGain=1;
+}
+///////////////////////////////////////////////////////////////////////////
+AliCalmodule::~AliCalmodule()
+{
+// Default destructor
+}
+///////////////////////////////////////////////////////////////////////////
+AliCalmodule::AliCalmodule(Int_t row,Int_t col,Float_t sig)
+{
+// Module constructor with initialisation of module data
+ fRow=row;
+ fCol=col;
+ fSignal=sig;
+ fSigc=sig;
+ fEdge=0;
+ fDead=0;
+ fGain=1;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalmodule::SetRow(Int_t i)
+{
+// Set the row number for this module
+ fRow=i;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalmodule::SetColumn(Int_t i)
+{
+// Set the column number for this module
+ fCol=i;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalmodule::SetSignal(Int_t row,Int_t col,Float_t sig)
+{
+// Set or change the data of the module
+ fRow=row;
+ fCol=col;
+ fSignal=sig;
+ fSigc=sig;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalmodule::AddSignal(Int_t row,Int_t col,Float_t sig)
+{
+// Add or change the data of the module
+ fRow=row;
+ fCol=col;
+ fSignal+=sig;
+ fSigc+=sig;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalmodule::SetClusteredSignal(Float_t sig)
+{
+// Set or change the signal of the module after clustering
+ fSigc=sig;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalmodule::SetEdgeOn()
+{
+// Indicate the module as edge module
+ fEdge=1;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalmodule::SetEdgeOff()
+{
+// Indicate the module as non-edge module
+ fEdge=0;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalmodule::EdgeUp()
+{
+// Increase the edge value by 1
+// This simplifies treatment of edge modules around temp. dead modules
+ fEdge+=1;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalmodule::EdgeDown()
+{
+// Decrease the edge value by 1
+// This simplifies treatment of edge modules around temp. dead modules
+ if (fEdge > 0) fEdge-=1;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalmodule::SetDead()
+{
+// Indicate the module as dead
+ fDead=1;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalmodule::SetAlive()
+{
+// Indicate the module as dead
+ fDead=0;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalmodule::SetGain(Float_t gain)
+{
+// Set the gain value of the readout system
+ fGain=gain;
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliCalmodule::GetRow()
+{
+// Provide the row number of the module
+ return fRow;
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliCalmodule::GetColumn()
+{
+// Provide the column number of the module
+ return fCol;
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliCalmodule::GetSignal()
+{
+// Provide the signal value of the module
+ if (!fDead)
+ {
+ return fSignal;
+ }
+ else
+ {
+ return 0;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliCalmodule::GetClusteredSignal()
+{
+// Provide the signal of the module after clustering
+ if (!fDead)
+ {
+ return fSigc;
+ }
+ else
+ {
+ return 0;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliCalmodule::GetEdgeValue()
+{
+// Provide the value of the edge indicator (1=edge 0=no-edge)
+ return fEdge;
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliCalmodule::GetDeadValue()
+{
+// Provide the value of the dead indicator (1=dead 0=alive)
+ return fDead;
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliCalmodule::GetGain()
+{
+// Provide the gain value of the readout system
+ return fGain;
+}
+///////////////////////////////////////////////////////////////////////////
--- /dev/null
+#ifndef ALICALMODULE_H
+#define ALICALMODULE_H
+///////////////////////////////////////////////////////////////////////////
+// Class AliCalmodule
+// Description of a module in a calorimeter system.
+// A matrix geometry is assumed, such that a module
+// is identified by (row,col) and contains a certain signal
+// Note : row and col start counting at 1
+//
+//--- NvE 13-jun-1997 UU-SAP Utrecht
+///////////////////////////////////////////////////////////////////////////
+
+#include <iostream.h>
+
+#include "AliSignal.h"
+
+class AliCalmodule : public AliSignal
+{
+ public:
+ AliCalmodule(); // Default constructor
+ ~AliCalmodule(); // Default destructor
+ AliCalmodule(Int_t row,Int_t col,Float_t sig); // Create a module and initialise data
+ void SetSignal(Int_t row,Int_t col,Float_t sig); // Set or change data for certain module
+ void AddSignal(Int_t row,Int_t col,Float_t sig); // Add signal to a certain module
+ void SetRow(Int_t i); // Set the row number of the module
+ void SetColumn(Int_t i); // Set the column number of the module
+ Int_t GetRow(); // Return the row number of the module
+ Int_t GetColumn(); // Return the column number of the module
+ Float_t GetSignal(); // Return the signal value of the module
+ void SetClusteredSignal(Float_t val); // Set the signal of the module after clustering
+ Float_t GetClusteredSignal(); // Return module signal after clustering
+ void SetEdgeOn(); // Set flag to indicate modules at edges
+ void SetEdgeOff(); // Set flag to indicate modules not at edges
+ void EdgeUp(); // Increase edge value by 1
+ void EdgeDown(); // Decrease edge value by 1
+ Int_t GetEdgeValue(); // Return the value of the edge indicator
+ void SetDead(); // Set flag to indicate dead modules
+ void SetAlive(); // Set flag to indicate active modules
+ Int_t GetDeadValue(); // Return the value of the dead module indicator
+ void SetGain(Float_t gain); // Set gain of the module's readout system
+ Float_t GetGain(); // Return the gain value
+
+ protected:
+ Int_t fRow; // The current row number
+ Int_t fCol; // The current column number
+ Float_t fSigc; // The signal after clustering
+ Int_t fEdge; // Flag to indicate edge module (>0=edge 0=no edge)
+ Int_t fDead; // Flag to indicate dead module (1=dead 0=alive)
+ Float_t fGain; // Gain of the module's readout system
+
+ ClassDef(AliCalmodule,1) // Class definition to enable ROOT I/O
+};
+#endif
--- /dev/null
+#include "AliCalorimeter.h"
+
+ClassImp(AliCalorimeter) // Class implementation to enable ROOT I/O
+
+AliCalorimeter::AliCalorimeter()
+{
+// Default constructor, all parameters set to 0
+ fNrows=0;
+ fNcolumns=0;
+ fNsignals=0;
+ fNclusters=0;
+ fMatrix=0;
+ fClusters=0;
+ fModules=0;
+ fHmodules=0;
+ fHclusters=0;
+ fNvetos=0;
+ fVetos=0;
+}
+///////////////////////////////////////////////////////////////////////////
+AliCalorimeter::~AliCalorimeter()
+{
+// Destructor to delete memory allocated for matrix and cluster array
+ if (fMatrix)
+ {
+ for (Int_t i=0; i<fNrows; i++)
+ {
+ delete [] fMatrix[i];
+ }
+ delete [] fMatrix;
+ }
+ fMatrix=0;
+ if (fModules) delete fModules;
+ fModules=0;
+ if (fClusters)
+ {
+ fClusters->Delete();
+ delete fClusters;
+ fClusters=0;
+ }
+ if (fHmodules) delete fHmodules;
+ fHmodules=0;
+ if (fHclusters) delete fHclusters;
+ fHclusters=0;
+ if (fVetos)
+ {
+ fVetos->Delete();
+ delete fVetos;
+ fVetos=0;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+AliCalorimeter::AliCalorimeter(Int_t nrow, Int_t ncol)
+{
+// Create a calorimeter module matrix
+ fNrows=nrow;
+ fNcolumns=ncol;
+ fNsignals=0;
+ fNclusters=0;
+ fClusters=0;
+ fMatrix=new AliCalmodule*[nrow];
+ for (Int_t i=0; i<nrow; i++)
+ {
+ fMatrix[i]=new AliCalmodule[ncol];
+ }
+ // Mark the edge modules
+ for (Int_t j=0; j<ncol; j++)
+ {
+ fMatrix[0][j].SetEdgeOn();
+ fMatrix[nrow-1][j].SetEdgeOn();
+ }
+ for (Int_t k=0; k<nrow; k++)
+ {
+ fMatrix[k][0].SetEdgeOn();
+ fMatrix[k][ncol-1].SetEdgeOn();
+ }
+
+ fModules=new TObjArray(); // Default size, expanded automatically
+
+ fHmodules=0;
+ fHclusters=0;
+
+ fNvetos=0;
+ fVetos=0;
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliCalorimeter::GetNrows()
+{
+// Provide the number of rows for the calorimeter module matrix
+ return fNrows;
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliCalorimeter::GetNcolumns()
+{
+// Provide the number of columns for the calorimeter module matrix
+ return fNcolumns;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalorimeter::SetSignal(Int_t row, Int_t col, Float_t sig)
+{
+// Set the signal for a certain calorimeter module
+
+ if (!fMatrix) LoadMatrix(); // Restore matrix data in case of reading input
+
+ if (row>0 && row<=fNrows && col>0 && col<=fNcolumns)
+ {
+ if (fMatrix[row-1][col-1].GetSignal() <= 0.) // only count new modules
+ {
+ fNsignals++;
+ fModules->Add(&(fMatrix[row-1][col-1]));
+ }
+ fMatrix[row-1][col-1].SetSignal(row,col,sig);
+ }
+ else
+ {
+ cout << " *AliCalorimeter::SetSignal* row,col : " << row << "," << col
+ << " out of range." << endl;
+ cout << " Nrows,Ncols = " << fNrows << "," << fNcolumns << endl;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalorimeter::AddSignal(Int_t row, Int_t col, Float_t sig)
+{
+// Add the signal to a certain calorimeter module
+
+ if (!fMatrix) LoadMatrix(); // Restore matrix data in case of reading input
+
+ if (row>0 && row<=fNrows && col>0 && col<=fNcolumns)
+ {
+ if (fMatrix[row-1][col-1].GetSignal() <= 0.) // only count new modules
+ {
+ fNsignals++;
+ fModules->Add(&(fMatrix[row-1][col-1]));
+ }
+ fMatrix[row-1][col-1].AddSignal(row,col,sig);
+ }
+ else
+ {
+ cout << " *AliCalorimeter::AddSignal* row,col : " << row << "," << col
+ << " out of range." << endl;
+ cout << " Nrows,Ncols = " << fNrows << "," << fNcolumns << endl;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalorimeter::Reset(Int_t row, Int_t col)
+{
+// Reset the signal for a certain calorimeter module
+
+ if (!fMatrix) LoadMatrix(); // Restore matrix data in case of reading input
+
+ if (row>0 && row<=fNrows && col>0 && col<=fNcolumns)
+ {
+ fMatrix[row-1][col-1].SetSignal(row,col,0);
+ fNsignals--;
+ fModules->Remove(&(fMatrix[row-1][col-1]));
+ }
+ else
+ {
+ cout << " *AliCalorimeter::Reset* row,col : " << row << "," << col
+ << " out of range." << endl;
+ cout << " Nrows,Ncols = " << fNrows << "," << fNcolumns << endl;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalorimeter::Reset()
+{
+// Reset the signals for the complete calorimeter
+// Normally this is done to prepare for the data of the next event
+// Note : Module gains, edge and dead flags remain unchanged
+
+ if (!fMatrix) LoadMatrix(); // Restore matrix data in case of reading input
+
+ fNsignals=0;
+ for (Int_t i=0; i<fNrows; i++)
+ {
+ for (Int_t j=0; j<fNcolumns; j++)
+ {
+ fMatrix[i][j].SetSignal(i+1,j+1,0);
+ }
+ }
+ if (fModules) fModules->Clear();
+
+ fNclusters=0;
+ if (fClusters)
+ {
+ fClusters->Delete();
+ delete fClusters;
+ fClusters=0;
+ }
+
+ fNvetos=0;
+ if (fVetos)
+ {
+ fVetos->Delete();
+ delete fVetos;
+ fVetos=0;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliCalorimeter::GetSignal(Int_t row, Int_t col)
+{
+// Provide the signal of a certain calorimeter module
+
+ if (!fMatrix) LoadMatrix(); // Restore matrix data in case of reading input
+
+ if (row>0 && row<=fNrows && col>0 && col<=fNcolumns)
+ {
+ return fMatrix[row-1][col-1].GetSignal();
+ }
+ else
+ {
+ cout << " *AliCalorimeter::GetSignal* row,col : " << row << "," << col
+ << " out of range." << endl;
+ cout << " Nrows,Ncols = " << fNrows << "," << fNcolumns << endl;
+ return 0;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalorimeter::SetEdgeOn(Int_t row, Int_t col)
+{
+// Indicate a certain calorimeter module as 'edge module'
+
+ if (!fMatrix) LoadMatrix(); // Restore matrix data in case of reading input
+
+ if (row>0 && row<=fNrows && col>0 && col<=fNcolumns)
+ {
+ fMatrix[row-1][col-1].SetEdgeOn();
+ }
+ else
+ {
+ cout << " *AliCalorimeter::SetEdgeOn* row,col : " << row << "," << col
+ << " out of range." << endl;
+ cout << " Nrows,Ncols = " << fNrows << "," << fNcolumns << endl;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalorimeter::SetEdgeOff(Int_t row, Int_t col)
+{
+// Indicate a certain calorimeter module as 'non-edge module'
+
+ if (!fMatrix) LoadMatrix(); // Restore matrix data in case of reading input
+
+ if (row>0 && row<=fNrows && col>0 && col<=fNcolumns)
+ {
+ fMatrix[row-1][col-1].SetEdgeOff();
+ }
+ else
+ {
+ cout << " *AliCalorimeter::SetEdgeOff* row,col : " << row << "," << col
+ << " out of range." << endl;
+ cout << " Nrows,Ncols = " << fNrows << "," << fNcolumns << endl;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalorimeter::SetDead(Int_t row, Int_t col)
+{
+// Indicate a certain calorimeter module as 'dead module'
+
+ if (!fMatrix) LoadMatrix(); // Restore matrix data in case of reading input
+
+ if (row>0 && row<=fNrows && col>0 && col<=fNcolumns)
+ {
+ fMatrix[row-1][col-1].SetDead();
+
+ // Increase the 'edge value' of surrounding modules
+ Int_t rlow=row-1;
+ Int_t rup=row+1;
+ Int_t clow=col-1;
+ Int_t cup=col+1;
+
+ if (rlow < 1) rlow=row;
+ if (rup > fNrows) rup=fNrows;
+ if (clow < 1) clow=col;
+ if (cup > fNcolumns) cup=fNcolumns;
+
+ for (Int_t i=rlow; i<=rup; i++)
+ {
+ for (Int_t j=clow; j<=cup; j++)
+ {
+ fMatrix[i-1][j-1].EdgeUp();
+ }
+ }
+
+ // Correct the 'edge value' for the dead module itself
+ fMatrix[row-1][col-1].EdgeDown();
+ }
+ else
+ {
+ cout << " *AliCalorimeter::SetDead* row,col : " << row << "," << col
+ << " out of range." << endl;
+ cout << " Nrows,Ncols = " << fNrows << "," << fNcolumns << endl;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalorimeter::SetAlive(Int_t row, Int_t col)
+{
+// Indicate a certain calorimeter module as 'active module'
+
+ if (!fMatrix) LoadMatrix(); // Restore matrix data in case of reading input
+
+ if (row>0 && row<=fNrows && col>0 && col<=fNcolumns)
+ {
+ fMatrix[row-1][col-1].SetAlive();
+
+ // Decrease the 'edge value' of surrounding modules
+ Int_t rlow=row-1;
+ Int_t rup=row+1;
+ Int_t clow=col-1;
+ Int_t cup=col+1;
+
+ if (rlow < 1) rlow=row;
+ if (rup > fNrows) rup=fNrows;
+ if (clow < 1) clow=col;
+ if (cup > fNcolumns) cup=fNcolumns;
+
+ for (Int_t i=rlow; i<=rup; i++)
+ {
+ for (Int_t j=clow; j<=cup; j++)
+ {
+ fMatrix[i-1][j-1].EdgeDown();
+ }
+ }
+
+ // Correct the 'edge value' for the dead module itself
+ fMatrix[row-1][col-1].EdgeUp();
+ }
+ else
+ {
+ cout << " *AliCalorimeter::SetAlive* row,col : " << row << "," << col
+ << " out of range." << endl;
+ cout << " Nrows,Ncols = " << fNrows << "," << fNcolumns << endl;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalorimeter::SetGain(Int_t row, Int_t col, Float_t gain)
+{
+// Set the gain value for a certain calorimeter module
+
+ if (!fMatrix) LoadMatrix(); // Restore matrix data in case of reading input
+
+ if (row>0 && row<=fNrows && col>0 && col<=fNcolumns)
+ {
+ fMatrix[row-1][col-1].SetGain(gain);
+ }
+ else
+ {
+ cout << " *AliCalorimeter::SetGain* row,col : " << row << "," << col
+ << " out of range." << endl;
+ cout << " Nrows,Ncols = " << fNrows << "," << fNcolumns << endl;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalorimeter::SetPosition(Int_t row,Int_t col,Float_t* vec,TString f)
+{
+// Set the position in user coordinates for a certain calorimeter module
+
+ if (!fMatrix) LoadMatrix(); // Restore matrix data in case of reading input
+
+ if (row>0 && row<=fNrows && col>0 && col<=fNcolumns)
+ {
+ fMatrix[row-1][col-1].SetPosition(vec,f);
+ }
+ else
+ {
+ cout << " *AliCalorimeter::SetPosition* row,col : " << row << "," << col
+ << " out of range." << endl;
+ cout << " Nrows,Ncols = " << fNrows << "," << fNcolumns << endl;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliCalorimeter::GetEdgeValue(Int_t row, Int_t col)
+{
+// Provide the value of the edge flag of a certain module
+
+ if (!fMatrix) LoadMatrix(); // Restore matrix data in case of reading input
+
+ if (row>0 && row<=fNrows && col>0 && col<=fNcolumns)
+ {
+ return fMatrix[row-1][col-1].GetEdgeValue();
+ }
+ else
+ {
+ cout << " *AliCalorimeter::GetEdgeValue* row,col : " << row << "," << col
+ << " out of range." << endl;
+ cout << " Nrows,Ncols = " << fNrows << "," << fNcolumns << endl;
+ return 0;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliCalorimeter::GetDeadValue(Int_t row, Int_t col)
+{
+// Provide the value of the dead flag of a certain module
+
+ if (!fMatrix) LoadMatrix(); // Restore matrix data in case of reading input
+
+ if (row>0 && row<=fNrows && col>0 && col<=fNcolumns)
+ {
+ return fMatrix[row-1][col-1].GetDeadValue();
+ }
+ else
+ {
+ cout << " *AliCalorimeter::GetDeadValue* row,col : " << row << "," << col
+ << " out of range." << endl;
+ cout << " Nrows,Ncols = " << fNrows << "," << fNcolumns << endl;
+ return 0;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliCalorimeter::GetGain(Int_t row, Int_t col)
+{
+// Provide the gain value of a certain module
+
+ if (!fMatrix) LoadMatrix(); // Restore matrix data in case of reading input
+
+ if (row>0 && row<=fNrows && col>0 && col<=fNcolumns)
+ {
+ return fMatrix[row-1][col-1].GetGain();
+ }
+ else
+ {
+ cout << " *AliCalorimeter::GetGain* row,col : " << row << "," << col
+ << " out of range." << endl;
+ cout << " Nrows,Ncols = " << fNrows << "," << fNcolumns << endl;
+ return 0;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalorimeter::GetPosition(Int_t row,Int_t col,Float_t* vec,TString f)
+{
+// Return the position in user coordinates for a certain calorimeter module
+
+ if (!fMatrix) LoadMatrix(); // Restore matrix data in case of reading input
+
+ if (row>0 && row<=fNrows && col>0 && col<=fNcolumns)
+ {
+ fMatrix[row-1][col-1].GetPosition(vec,f);
+ }
+ else
+ {
+ cout << " *AliCalorimeter::GetPosition* row,col : " << row << "," << col
+ << " out of range." << endl;
+ cout << " Nrows,Ncols = " << fNrows << "," << fNcolumns << endl;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliCalorimeter::GetClusteredSignal(Int_t row, Int_t col)
+{
+// Provide the module signal after clustering
+
+ if (!fMatrix) LoadMatrix(); // Restore matrix data in case of reading input
+
+ if (row>0 && row<=fNrows && col>0 && col<=fNcolumns)
+ {
+ return fMatrix[row-1][col-1].GetClusteredSignal();
+ }
+ else
+ {
+ cout << " *AliCalorimeter::GetClusteredSignal* row,col : " << row << "," << col
+ << " out of range." << endl;
+ cout << " Nrows,Ncols = " << fNrows << "," << fNcolumns << endl;
+ return 0;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliCalorimeter::GetNsignals()
+{
+// Provide the number of modules that contain a signal
+// Note : The number of modules marked 'dead' but which had a signal
+// are included.
+ return fNsignals;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalorimeter::Group(Int_t n)
+{
+// Group the individual modules into clusters
+// Module signals of n rings around the central module will be grouped
+
+ if (fNsignals > 0) // Directly return if no modules fired
+ {
+ if (!fMatrix) LoadMatrix(); // Restore matrix data in case of reading input
+
+ if (fNclusters > 0) Ungroup(); // Restore unclustered situation if needed
+
+ // Order the modules with decreasing signal
+ AliCalmodule* ordered=new AliCalmodule[fNsignals]; // temp. array for ordered modules
+ Sortm(ordered);
+
+ // Clustering of modules. Start with the highest signal.
+ if (fClusters)
+ {
+ fClusters->Delete();
+ delete fClusters;
+ fClusters=0;
+ }
+ fClusters=new TObjArray();
+ fNclusters=0;
+ Int_t row=0;
+ Int_t col=0;
+ AliCalcluster* c=0;
+ for (Int_t i=0; i<fNsignals; i++)
+ {
+ row=ordered[i].GetRow(); // row number of cluster center
+ col=ordered[i].GetColumn(); // column number of cluster center
+ if (row>0 && row<=fNrows && col>0 && col<=fNcolumns)
+ {
+ // only use modules not yet used in a cluster
+ if (fMatrix[row-1][col-1].GetClusteredSignal() > 0.)
+ {
+ c=new AliCalcluster;
+ c->Start(fMatrix[row-1][col-1]); // module to start the cluster
+ if (c->GetNmodules() > 0) // cluster started successfully (no edge)
+ {
+ fClusters->Add(c);
+ fNclusters++; // update cluster counter
+ AddRing(row,col,n); // add signals of n rings around the center
+ }
+ else
+ {
+ if (c) delete c;
+ c=0;
+ }
+ }
+ }
+ }
+
+ // Delete the temp. array
+ delete [] ordered;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalorimeter::Sortm(AliCalmodule* ordered)
+{
+// Order the modules with decreasing signal
+
+ Int_t nord=0;
+ for (Int_t i=0; i<fNrows; i++) // loop over all modules of the matrix
+ {
+ for (Int_t ii=0; ii<fNcolumns; ii++)
+ {
+ if (fMatrix[i][ii].GetSignal() <= 0.) continue; // only take modules with a signal
+
+ if (nord == 0) // store the first module with a signal at the first ordered position
+ {
+ nord++;
+ ordered[nord-1]=fMatrix[i][ii];
+ continue;
+ }
+
+ for (Int_t j=0; j<=nord; j++) // put module in the right ordered position
+ {
+ if (j == nord) // module has smallest signal seen so far
+ {
+ nord++;
+ ordered[j]=fMatrix[i][ii]; // add module at the end
+ break; // go for next matrix module
+ }
+
+ if (fMatrix[i][ii].GetSignal() < ordered[j].GetSignal()) continue;
+
+ nord++;
+ for (Int_t k=nord-1; k>j; k--) {ordered[k]=ordered[k-1];} // create empty position
+ ordered[j]=fMatrix[i][ii]; // put module at empty position
+ break; // go for next matrix module
+ }
+ }
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalorimeter::AddRing(Int_t row, Int_t col, Int_t n)
+{
+// Add module signals of 1 ring around (row,col) to current cluster
+// n denotes the maximum number of rings around cluster center
+// Note : This function is used recursively
+
+ if (n >= 1) // Check if any rings left for recursive calls
+ {
+ Float_t signal=GetSignal(row,col); // signal of (row,col) module
+
+ Int_t lrow=row-1; if (lrow < 1) lrow=1; // row lowerbound for ring
+ Int_t urow=row+1; if (urow > fNrows) urow=fNrows; // row upperbound for ring
+ Int_t lcol=col-1; if (lcol < 1) lcol=1; // col lowerbound for ring
+ Int_t ucol=col+1; if (ucol > fNcolumns) ucol=fNcolumns; // row upperbound for ring
+
+ for (Int_t i=lrow; i<=urow; i++)
+ {
+ for (Int_t j=lcol; j<=ucol; j++)
+ {
+ // add module(i,j) to cluster if the signal <= signal(row,col)
+ if (fMatrix[i-1][j-1].GetSignal() <= signal)
+ {
+ ((AliCalcluster*)fClusters->At(fNclusters-1))->Add(fMatrix[i-1][j-1]);
+ }
+ AddRing(i,j,n-1); // Go for ring of modules around this (i,j) one
+ }
+ }
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliCalorimeter::GetNclusters()
+{
+// Provide the number of clusters
+ return fNclusters;
+}
+///////////////////////////////////////////////////////////////////////////
+AliCalcluster* AliCalorimeter::GetCluster(Int_t j)
+{
+// Provide cluster number j
+// Note : j=1 denotes the first cluster
+ if ((j >= 1) && (j <= fNclusters))
+ {
+ return (AliCalcluster*)fClusters->At(j-1);
+ }
+ else
+ {
+ cout << " *AliCalorimeter::GetCluster* cluster number : " << j
+ << " out of range." << endl;
+ cout << " -- Cluster number 1 (if any) returned " << endl;
+ return (AliCalcluster*)fClusters->At(0);
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+AliCalmodule* AliCalorimeter::GetModule(Int_t j)
+{
+// Provide 'fired' module number j
+// Note : j=1 denotes the first 'fired' module
+ if ((j >= 1) && (j <= fNsignals))
+ {
+ return (AliCalmodule*)fModules->At(j-1);
+ }
+ else
+ {
+ cout << " *AliCalorimeter::GetModule* module number : " << j
+ << " out of range." << endl;
+ cout << " -- Fired module number 1 (if any) returned " << endl;
+ return (AliCalmodule*)fModules->At(0);
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+TH2F* AliCalorimeter::DrawModules()
+{
+// Provide a lego plot of the module signals
+
+ if (fHmodules)
+ {
+ fHmodules->Reset();
+ }
+ else
+ {
+ fHmodules=new TH2F("fHmodules","Module signals",
+ fNcolumns,0.5,float(fNcolumns)+0.5,fNrows,0.5,float(fNrows)+0.5);
+
+ fHmodules->SetDirectory(0); // Suppress global character of histo pointer
+ }
+
+ AliCalmodule* m;
+ Float_t row,col,signal;
+ for (Int_t i=0; i<fNsignals; i++)
+ {
+ m=(AliCalmodule*)fModules->At(i);
+ if (m)
+ {
+ row=float(m->GetRow());
+ col=float(m->GetColumn());
+ signal=m->GetSignal();
+ if (signal>0.) fHmodules->Fill(col,row,signal);
+ }
+ }
+
+ fHmodules->Draw("lego");
+ return fHmodules;
+}
+///////////////////////////////////////////////////////////////////////////
+TH2F* AliCalorimeter::DrawClusters()
+{
+// Provide a lego plot of the cluster signals
+
+ if (fHclusters)
+ {
+ fHclusters->Reset();
+ }
+ else
+ {
+ fHclusters=new TH2F("fHclusters","Cluster signals",
+ fNcolumns,0.5,float(fNcolumns)+0.5,fNrows,0.5,float(fNrows)+0.5);
+
+ fHclusters->SetDirectory(0); // Suppress global character of histo pointer
+ }
+
+ AliCalcluster* c;
+ Float_t row,col,signal;
+ for (Int_t i=0; i<fNclusters; i++)
+ {
+ c=(AliCalcluster*)fClusters->At(i);
+ if (c)
+ {
+ row=float(c->GetRow());
+ col=float(c->GetColumn());
+ signal=c->GetSignal();
+ if (signal>0.) fHclusters->Fill(col,row,signal);
+ }
+ }
+
+ fHclusters->Draw("lego");
+ return fHclusters;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalorimeter::LoadMatrix()
+{
+// Load the Calorimeter module matrix data back from the TObjArray
+
+ // Create the module matrix space
+ if (fMatrix)
+ {
+ for (Int_t k=0; k<fNrows; k++)
+ {
+ delete [] fMatrix[k];
+ }
+ delete [] fMatrix;
+ }
+ fMatrix=new AliCalmodule*[fNrows];
+ for (Int_t i=0; i<fNrows; i++)
+ {
+ fMatrix[i]=new AliCalmodule[fNcolumns];
+ }
+
+ // Copy the module data back into the matrix
+ AliCalmodule* m;
+ Int_t row;
+ Int_t col;
+ for (Int_t j=0; j<fNsignals; j++)
+ {
+ m=(AliCalmodule*)fModules->At(j);
+ row=m->GetRow();
+ col=m->GetColumn();
+ fMatrix[row-1][col-1]=*m;
+ fModules->AddAt(&(fMatrix[row-1][col-1]),j); // Store new pointer
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalorimeter::Ungroup()
+{
+// Set the module signals back to the non-clustered situation
+
+ if (!fMatrix) LoadMatrix(); // Restore matrix data in case of reading input
+
+ Float_t signal=0;
+ for (Int_t i=0; i<fNrows; i++)
+ {
+ for (Int_t j=0; j<fNcolumns; j++)
+ {
+ signal=fMatrix[i][j].GetSignal();
+ fMatrix[i][j].SetClusteredSignal(signal);
+ }
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliCalorimeter::AddVetoSignal(Float_t* r,TString f,Float_t s)
+{
+// Associate an (extrapolated) AliSignal at location r as veto to the cal.
+// Note : The default signal value (s) is 0
+ if (!fVetos)
+ {
+ fNvetos=0;
+ fVetos=new TObjArray();
+ }
+
+ fVetos->Add(new AliSignal);
+ fNvetos++;
+
+ ((AliSignal*)fVetos->At(fNvetos-1))->SetPosition(r,f);
+ ((AliSignal*)fVetos->At(fNvetos-1))->SetSignal(s);
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliCalorimeter::GetNvetos()
+{
+// Provide the number of veto signals associated to the calorimeter
+ return fNvetos;
+}
+///////////////////////////////////////////////////////////////////////////
+AliSignal* AliCalorimeter::GetVetoSignal(Int_t i)
+{
+// Provide access to the i-th veto signal of this calorimeter
+// Note : The first hit corresponds to i=1
+
+ if (i>0 && i<=fNvetos)
+ {
+ return (AliSignal*)fVetos->At(i-1);
+ }
+ else
+ {
+ cout << " *AliCalorimeter::GetVetoSignal* Signal number " << i
+ << " out of range." << endl;
+ cout << " --- First signal (if any) returned." << endl;
+ return (AliSignal*)fVetos->At(0);
+ }
+}
+///////////////////////////////////////////////////////////////////////////
--- /dev/null
+#ifndef ALICALORIMETER_H
+#define ALICALORIMETER_H
+///////////////////////////////////////////////////////////////////////////
+// Class AliCalorimeter
+// Description of a modular calorimeter system.
+// A matrix geometry is used in which a module is identified by (row,col).
+// Note : First module is identified as (1,1).
+//
+// This is the way to define and enter signals into a calorimeter :
+//
+// AliCalorimeter cal(10,15); // Calorimeter of 10x15 modules
+// // All module signals set to 0.
+// cal.AddSignal(5,7,85.4);
+// cal.AddSignal(5,7,25.9);
+// cal.AddSignal(3,5,1000);
+// cal.SetSignal(5,7,10.3);
+// cal.Reset(3,5); // Reset module (3,5) as being 'not fired'
+// // All module data are re-initialised.
+// cal.SetEdgeOn(1,1); // Declare module (1,1) as an 'edge module'
+// cal.SetDead(8,3);
+// cal.SetGain(2,8,3.2);
+//
+// Float_t vec[3]={6,1,20};
+// cal.SetPosition(2,8,vec,"car");
+//
+// Float_t loc[3]={-1,12,3};
+// cal.AddVetoSignal(loc,"car"); // Associate (extrapolated) position as a veto
+//
+// cal.Group(2); // Group 'fired' modules into clusters
+// // Perform grouping over 2 rings around the center
+// cal.Reset(); // Reset the complete calorimeter
+// // Normally to prepare for the next event data
+// // Note : Module gain, edge and dead flags remain
+//
+//--- NvE 13-jun-1997 UU-SAP Utrecht
+///////////////////////////////////////////////////////////////////////////
+
+#include <iostream.h>
+#include <math.h>
+
+#include "TObject.h"
+#include "TObjArray.h"
+#include "TH2.h"
+#include "TString.h"
+
+#include "AliDetector.h"
+
+#include "AliCalmodule.h"
+#include "AliCalcluster.h"
+#include "AliSignal.h"
+
+class AliCalorimeter : public AliDetector
+{
+ public:
+ AliCalorimeter(); // Default constructor
+ AliCalorimeter(Int_t nrow,Int_t ncol); // Create a calorimeter matrix
+ ~AliCalorimeter(); // Destructor
+ Int_t GetNrows(); // Return number of rows of the matrix
+ Int_t GetNcolumns(); // Return number of columns of the matrix
+ void SetSignal(Int_t row,Int_t col,Float_t s); // Set signal for a certain module
+ void AddSignal(Int_t row,Int_t col,Float_t s); // Add signal to a certain module
+ void Reset(Int_t row,Int_t col); // Reset signal for a certain module
+ void Reset(); // Reset the complete calorimeter
+ Float_t GetSignal(Int_t row,Int_t col); // Provide signal of a certain module
+ Int_t GetNsignals(); // Return number of modules with a signal
+ void Group(Int_t n); // Group modules into clusters (n rings)
+ Int_t GetNclusters(); // Return number of clusters
+ Float_t GetClusteredSignal(Int_t row,Int_t col); // Provide module signal after clustering
+ AliCalcluster* GetCluster(Int_t j); // Access to cluster number j
+ AliCalmodule* GetModule(Int_t j); // Access to 'fired' module number j
+ void SetEdgeOn(Int_t row,Int_t col); // Indicate module as 'edge module'
+ void SetEdgeOff(Int_t row,Int_t col); // Indicate module as 'non-edge module'
+ Int_t GetEdgeValue(Int_t row,Int_t col); // Provide the edge flag of a module
+ void SetDead(Int_t row,Int_t col); // Indicate module as 'dead module'
+ void SetAlive(Int_t row,Int_t col); // Indicate module as 'active module'
+ Int_t GetDeadValue(Int_t row,Int_t col); // Provide the dead flag of a module
+ void SetGain(Int_t row,Int_t col,Float_t g); // Set the gain value for a module
+ Float_t GetGain(Int_t row,Int_t col); // Provide the gain value of a module
+ void SetPosition(Int_t row,Int_t col,Float_t* r,TString f); // Set module position
+ void GetPosition(Int_t row,Int_t col,Float_t* r,TString f); // Return module position
+ TH2F* DrawModules(); // Draw lego plot of module signals
+ TH2F* DrawClusters(); // Draw lego plot of cluster signals
+ void AddVetoSignal(Float_t* r,TString f,Float_t s=0); // Associate (extrapolated) signal
+ AliSignal* GetVetoSignal(Int_t j); // Access to veto signal number j
+ Int_t GetNvetos(); // Provide the number of veto signals
+
+ protected:
+ Int_t fNrows; // The number of rows
+ Int_t fNcolumns; // The number of columns
+ Int_t fNsignals; // The number of modules with a signal
+ Int_t fNclusters; // The number of clusters
+ AliCalmodule** fMatrix; //! The matrix of modules for internal use
+ void Sortm(AliCalmodule*); // Order the modules with decreasing signal
+ TObjArray* fClusters; // The array of clusters
+ void AddRing(Int_t row,Int_t col,Int_t n); // add signals of n rings around cluster center
+ TObjArray* fModules; // The array of modules for output
+ void LoadMatrix(); // Load calorimeter matrix data from input
+ void Ungroup(); // Restore module matrix as before clustering
+ TH2F* fHmodules; //! The module 2-D histogram
+ TH2F* fHclusters; //! The cluster 2-D histogram
+ Int_t fNvetos; // The number of associated veto signals
+ TObjArray* fVetos; // The array of associated (extrapolated) veto signals
+
+ ClassDef(AliCalorimeter,1) // Class definition to enable ROOT I/O
+};
+#endif
--- /dev/null
+#include "AliInvmass.h"
+
+ClassImp(AliInvmass) // Class implementation to enable ROOT I/O
+
+AliInvmass::AliInvmass()
+{
+// Creation of an AliInvmass object and initialisation of parameters
+ fPi=acos(-1.);
+ fMode=2;
+ fBkg=0;
+ fNewtheta=1;
+ fNewphi=1;
+ fMinv=0;
+ fMbkg=0;
+}
+////////////////////////////////////////////////////////////////////////////////
+AliInvmass::~AliInvmass()
+{
+// Destructor to delete dynamically allocated memory
+ if (fMinv)
+ {
+ fMinv->Delete();
+ delete fMinv;
+ fMinv=0;
+ }
+
+ if (fMbkg)
+ {
+ fMbkg->Delete();
+ delete fMbkg;
+ fMbkg=0;
+ }
+}
+////////////////////////////////////////////////////////////////////////////////
+void AliInvmass::SetStorageMode(Int_t m)
+{
+// Set storage mode for the result arrays for inv. mass and comb. background
+ fMode=2;
+ if (m==1) fMode=1;
+}
+////////////////////////////////////////////////////////////////////////////////
+void AliInvmass::SetThetaSwitch(Int_t i)
+{
+// Enable/Disable (1/0) switching of theta angle in comb. bkg. reconstruction.
+// Default : Switching of theta is enabled.
+ fNewtheta=1;
+ if (i==0) fNewtheta=0;
+}
+////////////////////////////////////////////////////////////////////////////////
+void AliInvmass::SetPhiSwitch(Int_t i)
+{
+// Enable/Disable (1/0) switching of phi angle in comb. bkg. reconstruction.
+// Default : Switching of phi is enabled.
+ fNewphi=1;
+ if (i==0) fNewphi=0;
+}
+////////////////////////////////////////////////////////////////////////////////
+Int_t AliInvmass::GetStorageMode()
+{
+// Provide mode of storage for the result arrays for inv. mass and comb. background
+ return fMode;
+}
+////////////////////////////////////////////////////////////////////////////////
+Int_t AliInvmass::GetThetaSwitch()
+{
+// Provide the theta switching flag
+ return fNewtheta;
+}
+////////////////////////////////////////////////////////////////////////////////
+Int_t AliInvmass::GetPhiSwitch()
+{
+// Provide the phi switching flag
+ return fNewphi;
+}
+////////////////////////////////////////////////////////////////////////////////
+TObjArray* AliInvmass::Invmass(TObjArray* a1,TObjArray* a2)
+{
+// Perform two-particle invariant mass reconstruction
+ fBkg=0;
+ Combine(a1,a2);
+ return fMinv;
+}
+////////////////////////////////////////////////////////////////////////////////
+TObjArray* AliInvmass::CombBkg(TObjArray* a1,TObjArray* a2)
+{
+// Perform two-particle combinatorial background reconstruction
+ fBkg=1;
+ Combine(a1,a2);
+ if (fMode != 1)
+ {
+ return fMbkg;
+ }
+ else
+ {
+ return fMinv;
+ }
+}
+////////////////////////////////////////////////////////////////////////////////
+void AliInvmass::Combine(TObjArray* a1,TObjArray* a2)
+{
+// Perform two-particle invariant mass reconstruction
+
+ if ((!fBkg || fMode==1) && fMinv)
+ {
+ fMinv->Delete();
+ delete fMinv;
+ fMinv=0;
+ }
+
+ if (fBkg && (fMode !=1) && fMbkg)
+ {
+ fMbkg->Delete();
+ delete fMbkg;
+ fMbkg=0;
+ }
+
+ Int_t isame; // Indicates whether both lists are identical
+ isame=0;
+ if (a1==a2) isame=1;
+
+ // Index i must loop over the shortest of a1 and a2
+ TObjArray* listi=a1;
+ TObjArray* listj=a2;
+ Int_t ni=a1->GetEntries();
+ Int_t nj=a2->GetEntries();
+ if (nj < ni)
+ {
+ ni=a2->GetEntries();
+ nj=a1->GetEntries();
+ listi=a2;
+ listj=a1;
+ }
+
+ AliTrack* p1=0;
+ AliTrack* p2=0;
+ AliTrack* px=0;
+ Ali4Vector ptot;
+ AliTrack* t=0;
+ Double_t v2[4],vx[4];
+ Float_t q1,q2;
+
+ Int_t jmin; // Start index for list j
+ Int_t jx; // Index for randomly picked particle for comb. bkg. reconstruction
+
+ for (Int_t i=0; i<ni; i++) // Select first a particle from list i
+ {
+ p1=(AliTrack*)listi->At(i);
+ p2=0;
+
+ if (!p1) continue;
+
+ jmin=0;
+ if (isame) jmin=i+1;
+ for (Int_t j=jmin; j<nj; j++) // Select also a particle from list j
+ {
+ p2=(AliTrack*)listj->At(j);
+ if (p1==p2) p2=0; // Don't combine particle with itself
+
+ if (!p2) continue;
+
+ p2->GetVector(v2,"sph");
+
+ // Take theta and phi from randomly chosen other list j particle for bkg. reconstr.
+ if (fBkg)
+ {
+ px=0;
+ if ((!isame && nj>1) || (isame && nj>2))
+ {
+ jx=int(fRndm.Uniform(0,float(nj)));
+ px=(AliTrack*)listj->At(jx);
+
+ while (!px || px==p2 || px==p1)
+ {
+ jx++;
+ if (jx >= nj) jx=0;
+ px=(AliTrack*)listj->At(jx);
+ }
+
+ px->GetVector(vx,"sph");
+ if (fNewtheta) v2[2]=vx[2]; // Replace the theta angle in the v2 vector
+ if (fNewphi) v2[3]=vx[3]; // Replace the phi angle in the v2 vector
+ }
+ }
+
+ if ((!fBkg && p2) || (fBkg && px))
+ {
+ // Store the data of this two-particle combination
+ ptot.SetVector(v2,"sph");
+ ptot=(Ali4Vector)(ptot+(*p1));
+ q1=p1->GetCharge();
+ q2=p2->GetCharge();
+ t=new AliTrack;
+ t->Set4Momentum(ptot);
+ t->SetCharge(q1+q2);
+ if (!fBkg || fMode==1)
+ {
+ if (!fMinv) fMinv=new TObjArray();
+ fMinv->Add(t);
+ }
+ else
+ {
+ if (!fMbkg) fMbkg=new TObjArray();
+ fMbkg->Add(t);
+ }
+ }
+
+ } // End of second particle loop
+
+ } // End of first particle loop
+}
+////////////////////////////////////////////////////////////////////////////////
--- /dev/null
+#ifndef ALIINVMASS_H
+#define ALIINVMASS_H
+////////////////////////////////////////////////////////////////////////////////
+// Class AliInvmass
+// Construction of invariant mass and combinatorial background.
+//
+// Example :
+// ---------
+//
+// TObjArray* photons=new TObjArray(); // Array with photon tracks for pi0 rec.
+//
+// // Code to create some photon tracks from pi0 decays
+// Int_t ntracks=200;
+// for (Int_t i=0; i<ntracks; i++)
+// {
+// photons->Add(new Alitrack);
+// ...
+// ...
+// ...
+// }
+//
+// // Perform the invariant mass and comb. bkg. reconstruction
+//
+// TObjArray* allm=q.Invmass(photons,photons); // All reconstructed invariant masses
+//
+// TH1F* hall=new TH1F("hall","hall",200,0,2); // Histo with M_inv of all combinations
+//
+// Int_t nall=0;
+// if (allm) nall=allm->GetEntries();
+//
+// AliTrack* t;
+// Float_t minv;
+// for (Int_t j=0; j<nall; j++)
+// {
+// t=(AliTrack*)allm->At(j);
+// if (t)
+// {
+// minv=t->GetMass();
+// hall->Fill(minv);
+// }
+// }
+//
+// TObjArray* bkgm=q.CombBkg(photons,photons); // Reconstructed comb. background
+//
+// TH1F* hbkg=new TH1F("hbkg","hbkg",200,0,2); // Histo with M_inv. of comb. background
+//
+// Int_t nbkg=0;
+// if (bkgm) nbkg=bkgm->GetEntries();
+//
+// for (Int_t j=0; j<nbkg; j++)
+// {
+// t=(AliTrack*)bkgm->At(j);
+// if (t)
+// {
+// minv=t->GetMass();
+// hbkg->Fill(minv);
+// }
+// }
+//
+// TH1F* hsig=new TH1F("sig","sig",200,0,2); // Histo with the bkg. subtracted signal
+// hsig->Sumw2();
+// hsig->Add(hall,hbkg,1,-1);
+//
+//
+// Note : By default the storage of the reconstructed information is performed
+// in separate TObjArrays for the signal and comb. background resp.
+// In order to limit the memory usage, AliInvmass::SetStorageMode(1) may be
+// used to activate only a single TObjArray to store the reconstructed information.
+// Consequently, the following statements
+//
+// TObjArray* allm=q.Invmass(photons,photons);
+// TObjArray* bkgm=q.CombBkg(photons,photons);
+//
+// will result in the fact that after he invokation of CombBkg
+// the information of "allm" is lost due to the fact that the storage is
+// is re-used for "bkgm" in case the "single storage" option has been selected.
+// Usage of the, in that case invalid, pointer "allm" may cause your
+// program to crash.
+//
+// * Thus : In case of single storage usage, all invokations of the returned
+// array pointer have to be completed before invoking any memberfunction
+// of the same AliInvmass object again.
+//
+//
+//
+//--- NvE 12-apr-1999 UU-SAP Utrecht
+////////////////////////////////////////////////////////////////////////////////
+
+#include <iostream.h>
+#include <math.h>
+
+#include "TObject.h"
+#include "TObjArray.h"
+
+#include "AliRandom.h"
+#include "AliTrack.h"
+
+class AliInvmass : public TObject
+{
+ public:
+ AliInvmass(); // Default constructor
+ ~AliInvmass(); // Destructor
+ void SetStorageMode(Int_t m); // Set storage mode (1=single, 2=multiple)
+ void SetThetaSwitch(Int_t i=1); // Enable (1/0) new theta for comb. bkg. reco.
+ void SetPhiSwitch(Int_t i=1); // Enable (1/0) new phi for comb. bkg. reco.
+ Int_t GetStorageMode(); // Provide storage mode
+ Int_t GetThetaSwitch(); // Provide theta switch flag
+ Int_t GetPhiSwitch(); // Provide phi switch flag
+ TObjArray* Invmass(TObjArray* a1,TObjArray* a2); // Two-particle inv. mass reco.
+ TObjArray* CombBkg(TObjArray* a1,TObjArray* a2); // Two-particle comb. background reco.
+
+ protected:
+ Double_t fPi; // Value of pi
+ Int_t fMode; // Storage mode for signal and bkg. results (2=separate arrays)
+ Int_t fBkg; // Flag to denote comb. background processing
+ AliRandom fRndm; // The random number generator for the comb. bkg. reconstruction
+ Int_t fNewtheta; // Flag to denote enabling of switching theta for comb. bkg. reco.
+ Int_t fNewphi; // Flag to denote enabling of switching phi for comb. bkg. reco.
+ TObjArray* fMinv; // Array with reconstructed invariant mass 'tracks'
+ TObjArray* fMbkg; // Array with reconstructed comb. background 'tracks'
+
+ private:
+ void Combine(TObjArray* a1,TObjArray* a2); // Make two-particle combinations
+
+ ClassDef(AliInvmass,1) // Class definition to enable ROOT I/O
+};
+#endif
--- /dev/null
+#include "AliJet.h"
+
+ClassImp(AliJet) // Class implementation to enable ROOT I/O
+
+AliJet::AliJet()
+{
+// Default constructor
+// All variables initialised to 0
+// Initial maximum number of tracks is set to the default value
+ fTracks=0;
+ fNtinit=0;
+ Reset();
+ SetNtinit();
+}
+///////////////////////////////////////////////////////////////////////////
+AliJet::AliJet(Int_t n)
+{
+// Create a jet to hold initially a maximum of n tracks
+// All variables initialised to 0
+ fTracks=0;
+ fNtinit=0;
+ Reset();
+ if (n > 0)
+ {
+ SetNtinit(n);
+ }
+ else
+ {
+ cout << endl;
+ cout << " *AliJet* Initial max. number of tracks entered : " << n << endl;
+ cout << " This is invalid. Default initial maximum will be used." << endl;
+ cout << endl;
+ SetNtinit();
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+AliJet::~AliJet()
+{
+// Default destructor
+ if (fTracks) delete fTracks;
+ fTracks=0;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliJet::SetNtinit(Int_t n)
+{
+// Set the initial maximum number of tracks for this jet
+ fNtinit=n;
+ fNtmax=n;
+ if (fTracks) delete fTracks;
+ fTracks=new TObjArray(fNtmax);
+}
+///////////////////////////////////////////////////////////////////////////
+void AliJet::Reset()
+{
+// Reset all variables to 0
+// The max. number of tracks is set to the initial value again
+ fNtrk=0;
+ fQ=0;
+ Double_t a[4]={0,0,0,0};
+ SetVector(a,"sph");
+ if (fNtinit > 0) SetNtinit(fNtinit);
+}
+///////////////////////////////////////////////////////////////////////////
+void AliJet::Add(AliTrack& t)
+{
+// Add a track to the jet
+// In case the maximum number of tracks has been reached
+// space will be extended to hold an additional amount of tracks as
+// was initially reserved
+ if (fNtrk == fNtmax) // Check if maximum track number is reached
+ {
+ fNtmax+=fNtinit;
+ fTracks->Expand(fNtmax);
+ }
+
+ // Add the track to this jet
+ fNtrk++;
+ fTracks->Add(&t);
+ (*this)+=(Ali4Vector&)t;
+ fQ+=t.GetCharge();
+}
+///////////////////////////////////////////////////////////////////////////
+void AliJet::Info(TString f)
+{
+// Provide jet information within the coordinate frame f
+ cout << " *AliJet::Info* Invmass : " << GetInvmass() << " Charge : " << fQ
+ << " Momentum : " << GetMomentum() << " Ntracks : " << fNtrk << endl;
+ cout << " ";
+ Ali4Vector::Info(f);
+}
+///////////////////////////////////////////////////////////////////////////
+void AliJet::List(TString f)
+{
+// Provide jet and primary track information within the coordinate frame f
+
+ Info(f); // Information of the current jet
+
+ // The tracks of this jet
+ AliTrack* t;
+ for (Int_t it=1; it<=fNtrk; it++)
+ {
+ t=GetTrack(it);
+ if (t)
+ {
+ cout << " ---Track no. " << it << endl;
+ cout << " ";
+ t->Info(f);
+ }
+ else
+ {
+ cout << " *AliJet::List* Error : No track present." << endl;
+ }
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliJet::ListAll(TString f)
+{
+// Provide jet and prim.+sec. track information within the coordinate frame f
+
+ Info(f); // Information of the current jet
+
+ // The tracks of this jet
+ AliTrack* t;
+ for (Int_t it=1; it<=fNtrk; it++)
+ {
+ t=GetTrack(it);
+ if (t)
+ {
+ cout << " ---Track no. " << it << endl;
+ cout << " ";
+ t->ListAll(f);
+ }
+ else
+ {
+ cout << " *AliJet::List* Error : No track present." << endl;
+ }
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliJet::GetNtracks()
+{
+// Return the current number of tracks of this jet
+ return fNtrk;
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t AliJet::GetEnergy()
+{
+// Return the total energy of the jet
+ return GetScalar();
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t AliJet::GetMomentum()
+{
+// Return the value of the total jet 3-momentum
+ Ali3Vector p=Get3Vector();
+ Double_t p2=p.Dot(p);
+ return sqrt(p2);
+}
+///////////////////////////////////////////////////////////////////////////
+Ali3Vector AliJet::Get3Momentum()
+{
+// Return the the total jet 3-momentum
+ Ali3Vector p=Get3Vector();
+ return p;
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t AliJet::GetInvmass()
+{
+// Return the invariant mass of the jet
+ Double_t m2=Dot(*this);
+ if (m2>0)
+ {
+ return sqrt(m2);
+ }
+ else
+ {
+ return 0;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliJet::GetCharge()
+{
+// Return the total charge of the jet
+ return fQ;
+}
+///////////////////////////////////////////////////////////////////////////
+AliTrack* AliJet::GetTrack(Int_t i)
+{
+// Return the i-th track of this jet
+ return (AliTrack*)fTracks->At(i-1);
+}
+///////////////////////////////////////////////////////////////////////////
--- /dev/null
+#ifndef ALIJET_H
+#define ALIJET_H
+///////////////////////////////////////////////////////////////////////////
+// Class AliJet
+// Creation and investigation of a jet of particle tracks.
+// An AliJet can be constructed by adding AliTracks.
+//
+// Coding example to make 2 jets j1 and j2.
+// ----------------------------------------
+// j1 contains the AliTracks 1 and 2
+// j2 contains the AliTracks 3 and 4
+//
+// AliTrack t1,t2,t3,t4;
+// ...
+// ... // code to fill the AliTrack data
+// ...
+// AliJet j1(5);
+// AliJet j2(12);
+// j1.Add(t1);
+// j1.Add(t2);
+// j2.Add(t3);
+// j2.Add(t4);
+//
+// j1.Info();
+// j2.Info("sph");
+//
+// Float_t e1=j1.GetEnergy();
+// Float_t pnorm=j1->GetMomentum();
+// Ali3Vector p=j1->Get3Momentum();
+// Float_t m=j1.GetInvmass();
+// Int_t ntk=j1.GetNtracks();
+// AliTrack* tj=j1.GetTrack(1);
+//
+// Note : All quantities are in GeV, GeV/c or GeV/c**2
+//
+//--- NvE 10-jul-1997 UU-SAP Utrecht
+//--- Modified : NvE 06-apr-1999 UU-SAP Utrecht to inherit from Ali4Vector
+///////////////////////////////////////////////////////////////////////////
+
+#include <iostream.h>
+#include <math.h>
+
+#include "TObject.h"
+#include "TObjArray.h"
+
+#include "Ali4Vector.h"
+#include "AliTrack.h"
+
+class AliJet : public TObject,public Ali4Vector
+{
+ public:
+ AliJet(); // Default constructor
+ AliJet(Int_t n); // Create a Jet to hold initially n Tracks
+ ~AliJet(); // Default destructor
+ void Reset(); // Reset all values
+ virtual void Add(AliTrack& t); // Add a track to the jet
+ virtual void Add(AliTrack* t) { Add(*t); }
+ void Info(TString f); // Print jet information in coordinate frame f
+ void List(TString f="car"); // Print jet prim. track information for coord. frame f
+ void ListAll(TString f="car"); // Print jet prim. and decay track information for coord. frame f
+ Double_t GetEnergy(); // Provide the total jet energy
+ Double_t GetMomentum(); // Provide the value of the total jet 3-momentum
+ Ali3Vector Get3Momentum(); // Provide the total jet 3-momentum
+ Double_t GetInvmass(); // Provide the invariant mass
+ Float_t GetCharge(); // Provide the total charge of the jet
+ Int_t GetNtracks(); // Return the number of tracks in the jet
+ AliTrack* GetTrack(Int_t i); // Provide i-th track of the jet (1=first track)
+
+ protected:
+ void SetNtinit(Int_t n=2); // Set the initial max. number of tracks for this Jet
+ Int_t fNtinit; // The initial max. number of tracks for this jet
+ Int_t fNtmax; // The maximum number of tracks for this Jet
+ Float_t fQ; // The total charge of the jet
+ Int_t fNtrk; // The number of tracks in the jet
+ TObjArray* fTracks; // Array to hold the pointers to the tracks of the jet
+
+ ClassDef(AliJet,1) // Class definition to enable ROOT I/O
+};
+#endif
--- /dev/null
+#include "AliMath.h"
+
+ClassImp(AliMath) // Class implementation to enable ROOT I/O
+
+AliMath::AliMath()
+{
+// Default constructor
+}
+///////////////////////////////////////////////////////////////////////////
+AliMath::~AliMath()
+{
+// Destructor
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliMath::Gamma(Float_t z)
+{
+// Computation of gamma(z) for all z>0.
+//
+// The algorithm is based on the article by C.Lanczos [1] as denoted in
+// Numerical Recipes 2nd ed. on p. 207 (W.H.Press et al.).
+//
+// [1] C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86.
+//
+//--- Nve 14-nov-1998 UU-SAP Utrecht
+
+ if (z<=0.)
+ {
+ cout << "*Gamma(z)* Wrong argument z = " << z << endl;
+ return 0;
+ }
+
+ Float_t v=LnGamma(z);
+ return exp(v);
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliMath::Gamma(Float_t a,Float_t x)
+{
+// Computation of the incomplete gamma function P(a,x)
+//
+// The algorithm is based on the formulas and code as denoted in
+// Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
+//
+//--- Nve 14-nov-1998 UU-SAP Utrecht
+
+ if (a<=0.)
+ {
+ cout << "*Gamma(a,x)* Invalid argument a = " << a << endl;
+ return 0;
+ }
+
+ if (x<=0.)
+ {
+ if (x<0) cout << "*Gamma(a,x)* Invalid argument x = " << x << endl;
+ return 0;
+ }
+
+ if (x<(a+1.))
+ {
+ return GamSer(a,x);
+ }
+ else
+ {
+ return GamCf(a,x);
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliMath::LnGamma(Float_t z)
+{
+// Computation of ln[gamma(z)] for all z>0.
+//
+// The algorithm is based on the article by C.Lanczos [1] as denoted in
+// Numerical Recipes 2nd ed. on p. 207 (W.H.Press et al.).
+//
+// [1] C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86.
+//
+// The accuracy of the result is better than 2e-10.
+//
+//--- Nve 14-nov-1998 UU-SAP Utrecht
+
+ if (z<=0.)
+ {
+ cout << "*LnGamma(z)* Wrong argument z = " << z << endl;
+ return 0;
+ }
+
+ // Coefficients for the series expansion
+ Double_t c[7];
+ c[0]= 2.5066282746310005;
+ c[1]= 76.18009172947146;
+ c[2]=-86.50532032941677;
+ c[3]= 24.01409824083091;
+ c[4]= -1.231739572450155;
+ c[5]= 0.1208650973866179e-2;
+ c[6]= -0.5395239384953e-5;
+
+ Double_t x=z;
+ Double_t y=x;
+ Double_t tmp=x+5.5;
+ tmp=(x+0.5)*log(tmp)-tmp;
+ Double_t ser=1.000000000190015;
+ for (Int_t i=1; i<7; i++)
+ {
+ y+=1.;
+ ser+=c[i]/y;
+ }
+ Float_t v=tmp+log(c[0]*ser/x);
+ return v;
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliMath::GamSer(Float_t a,Float_t x)
+{
+// Computation of the incomplete gamma function P(a,x)
+// via its series representation.
+//
+// The algorithm is based on the formulas and code as denoted in
+// Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
+//
+//--- Nve 14-nov-1998 UU-SAP Utrecht
+
+ Int_t itmax=100; // Maximum number of iterations
+ Float_t eps=3.e-7; // Relative accuracy
+
+ if (a<=0.)
+ {
+ cout << "*GamSer(a,x)* Invalid argument a = " << a << endl;
+ return 0;
+ }
+
+ if (x<=0.)
+ {
+ if (x<0) cout << "*GamSer(a,x)* Invalid argument x = " << x << endl;
+ return 0;
+ }
+
+ Float_t gln=LnGamma(a);
+ Float_t ap=a;
+ Float_t sum=1./a;
+ Float_t del=sum;
+ for (Int_t n=1; n<=itmax; n++)
+ {
+ ap+=1.;
+ del=del*x/ap;
+ sum+=del;
+ 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);
+ return v;
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliMath::GamCf(Float_t a,Float_t x)
+{
+// Computation of the incomplete gamma function P(a,x)
+// via its continued fraction representation.
+//
+// The algorithm is based on the formulas and code as denoted in
+// Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
+//
+//--- 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
+
+ if (a<=0.)
+ {
+ cout << "*GamCf(a,x)* Invalid argument a = " << a << endl;
+ return 0;
+ }
+
+ if (x<=0.)
+ {
+ if (x<0) cout << "*GamCf(a,x)* Invalid argument x = " << x << endl;
+ 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;
+ for (Int_t i=1; i<=itmax; i++)
+ {
+ an=float(-i)*(float(i)-a);
+ b+=2.;
+ d=an*d+b;
+ if (fabs(d)<fpmin) d=fpmin;
+ c=b+an/c;
+ if (fabs(c)<fpmin) c=fpmin;
+ d=1./d;
+ del=d*c;
+ h=h*del;
+ 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;
+ return (1.-v);
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliMath::Erf(Float_t x)
+{
+// Computation of the error function erf(x).
+//
+//--- NvE 14-nov-1998 UU-SAP Utrecht
+
+ return (1.-Erfc(x));
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliMath::Erfc(Float_t x)
+{
+// Computation of the complementary error function erfc(x).
+//
+// The algorithm is based on a Chebyshev fit as denoted in
+// Numerical Recipes 2nd ed. on p. 214 (W.H.Press et al.).
+//
+// The fractional error is always less than 1.2e-7.
+//
+//--- 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;
+
+ Float_t v=1.; // The return value
+
+ Float_t z=fabs(x);
+
+ if (z <= 0.) return v; // erfc(0)=1
+
+ Float_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)))))))));
+
+ if (x < 0.) v=2.-v; // erfc(-x)=2-erfc(x)
+
+ return v;
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliMath::Prob(Float_t chi2,Int_t ndf)
+{
+// Computation of the probability for a certain Chi-squared (chi2)
+// and number of degrees of freedom (ndf).
+//
+// Calculations are based on the incomplete gamma function P(a,x),
+// where a=ndf/2 and x=chi2/2.
+//
+// P(a,x) represents the probability that the observed Chi-squared
+// for a correct model should be less than the value chi2.
+//
+// The returned probability corresponds to 1-P(a,x),
+// which denotes the probability that an observed Chi-squared exceeds
+// the value chi2 by chance, even for a correct model.
+//
+//--- NvE 14-nov-1998 UU-SAP Utrecht
+
+ if (ndf <= 0) return 0; // Set CL to zero in case ndf<=0
+
+ if (chi2 <= 0.)
+ {
+ if (chi2 < 0.)
+ {
+ return 0;
+ }
+ else
+ {
+ return 1;
+ }
+ }
+
+// Alternative which is exact
+// This code may be activated in case the gamma function gives problems
+// if (ndf==1)
+// {
+// Float_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));
+// if (n>30 && q>0.)
+// {
+// Float_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.;
+ return (1.-Gamma(a,x));
+}
+///////////////////////////////////////////////////////////////////////////
--- /dev/null
+#ifndef ALIMATH_H
+#define ALIMATH_H
+///////////////////////////////////////////////////////////////////////////
+// Class AliMath
+// Various mathematical tools which may be very convenient while
+// performing physics analysis.
+//
+// Example : Probability of a Chi-squared value
+// =========
+//
+// AliMath M;
+// Float_t chi2=20 ; // The chi-squared value
+// Int_t ndf=12; // The number of degrees of freedom
+// Float_t p=M.Prob(chi2,ndf); // The probability that at least a Chi-squared
+// // value of chi2 will be observed, even for a
+// // correct model
+//
+//--- NvE 14-nov-1998 UU-SAP Utrecht
+///////////////////////////////////////////////////////////////////////////
+
+#include <iostream.h>
+#include <math.h>
+
+#include "TObject.h"
+
+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
+
+ 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
+
+ ClassDef(AliMath,1) // Class definition to enable ROOT I/O
+
+};
+#endif
--- /dev/null
+#include "AliPosition.h"
+
+ClassImp(AliPosition) // Class implementation to enable ROOT I/O
+
+AliPosition::AliPosition()
+{
+// Creation of an AliPosition object and initialisation of parameters
+}
+///////////////////////////////////////////////////////////////////////////
+AliPosition::~AliPosition()
+{
+// Destructor to delete dynamically allocated memory
+}
+///////////////////////////////////////////////////////////////////////////
+void AliPosition::SetPosition(Double_t* r,TString f)
+{
+// Store position according to reference frame f
+ SetVector(r,f);
+}
+///////////////////////////////////////////////////////////////////////////
+void AliPosition::GetPosition(Double_t* r,TString f)
+{
+// Provide position according to reference frame f
+ GetVector(r,f);
+}
+///////////////////////////////////////////////////////////////////////////
+void AliPosition::SetPosition(Float_t* r,TString f)
+{
+// Store position according to reference frame f
+ SetVector(r,f);
+}
+///////////////////////////////////////////////////////////////////////////
+void AliPosition::GetPosition(Float_t* r,TString f)
+{
+// Provide position according to reference frame f
+ GetVector(r,f);
+}
+///////////////////////////////////////////////////////////////////////////
+AliPosition& AliPosition::GetPosition()
+{
+// Provide position
+ return (*this);
+}
+///////////////////////////////////////////////////////////////////////////
+void AliPosition::SetPosition(Ali3Vector& r)
+{
+// Set position
+ Double_t a[3];
+ r.GetVector(a,"sph");
+ SetVector(a,"sph");
+}
+///////////////////////////////////////////////////////////////////////////
--- /dev/null
+#ifndef ALIPOSITION_H
+#define ALIPOSITION_H
+///////////////////////////////////////////////////////////////////////////
+// Class AliPosition
+// Handling of positions in various reference frames.
+//
+// This class is meant to serve as a base class for ALICE objects
+// that have a unique position in 3-dimensional space.
+//
+// Note :
+// ------
+// Positions (r) and reference frames (f) are specified via
+// SetPosition(Float_t* r,TString f) under the following conventions :
+//
+// f="car" ==> r in Cartesian coordinates (x,y,z)
+// f="sph" ==> r in Spherical coordinates (r,theta,phi)
+// f="cyl" ==> r in Cylindrical coordinates (rho,phi,z)
+//
+// All angles are in radians.
+//
+// Example :
+// ---------
+//
+// AliPosition q;
+// Float_t pos[3]={-1,25,7};
+// q.SetPosition(pos,"car");
+// Float_t loc[3];
+// q.GetPosition(loc,"sph");
+//
+//--- NvE 06-feb-1999 UU-SAP Utrecht
+///////////////////////////////////////////////////////////////////////////
+
+#include <iostream.h>
+#include <math.h>
+
+#include "TObject.h"
+#include "TString.h"
+
+#include "Ali3Vector.h"
+
+class AliPosition : public Ali3Vector
+{
+ public:
+ AliPosition(); // Default constructor
+ virtual ~AliPosition(); // Destructor
+ virtual void SetPosition(Double_t* r,TString f); // Store position r in frame f
+ virtual void GetPosition(Double_t* r,TString f); // Provide position r in frame f
+ virtual void SetPosition(Float_t* r,TString f); // Store position r in frame f
+ virtual void GetPosition(Float_t* r,TString f); // Provide position r in frame f
+ AliPosition& GetPosition(); // Provide position
+ virtual void SetPosition(Ali3Vector& r); // Store position r
+
+ ClassDef(AliPosition,1) // Class definition to enable ROOT I/O
+};
+#endif
--- /dev/null
+#include "AliRandom.h"
+
+ClassImp(AliRandom) // Class implementation to enable ROOT I/O
+
+AliRandom::AliRandom()
+{
+// Creation of an AliRandom object and default initialisation
+//
+// A seed is used to create the initial u[97] table.
+// This seed is converted into four startup parameters i j k and l
+// (see member function "unpack").
+//
+// Suggested test values : i=12 j=34 k=56 l=78 (see article)
+// which corresponds to : seed = 53310452
+
+ Int_t seed=53310452; // Default seed
+ Start(seed,0,0); // Start the sequence for this seed from scratch
+}
+///////////////////////////////////////////////////////////////////////////
+AliRandom::AliRandom(Int_t seed)
+{
+// Creation of an AliRandom object and user defined initialisation
+
+ Start(seed,0,0); // Start the sequence for this seed from scratch
+}
+///////////////////////////////////////////////////////////////////////////
+AliRandom::AliRandom(Int_t seed,Int_t cnt1,Int_t cnt2)
+{
+// Creation of an AliRandom object and user defined initialisation
+//
+// seed is the seed to create the initial u[97] table.
+// The range of the seed is : 0 <= seed <= 921350143
+//
+// cnt1 and cnt2 are the parameters for the counting system
+// to enable a start of the sequence at a certain point.
+// The current values of seed, cnt1 and cnt2 can be obtained
+// via the member functions "GetSeed", "GetCnt1" and "GetCnt2" resp.
+// To start from scratch one should select : cnt1=0 and cnt2=0
+
+ Start(seed,cnt1,cnt2); // Start the sequence from a user defined point
+}
+///////////////////////////////////////////////////////////////////////////
+AliRandom::~AliRandom()
+{
+// Destructor to delete memory allocated for the area function arrays
+ if (fXa) delete [] fXa;
+ fXa=0;
+ if (fYa) delete [] fYa;
+ fYa=0;
+ if (fIbins) delete [] fIbins;
+ fIbins=0;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliRandom::Start(Int_t seed,Int_t cnt1,Int_t cnt2)
+{
+// Start a certain sequence from scratch or from a user defined point
+//
+// The algorithm to start from scratch is based on the routine RSTART
+// as described in the report by G.Marsaglia and A.Zaman
+// (FSU-SCRI-87-50 Florida State University 1987).
+//
+// seed is the seed to create the initial u[97] table.
+// This seed is converted into four startup parameters i j k and l
+// (see the member function "unpack").
+//
+// The range of the seed is : 0 <= seed <= 921350143
+//
+// Suggested test values : i=12 j=34 k=56 l=78 (see article)
+// which corresponds to : seed = 53310452
+//
+// cnt1 and cnt2 are the parameters for the counting system
+// to enable a start of the sequence at a certain point.
+// The current values of seed, cnt1 and cnt2 can be obtained
+// via the member functions "GetSeed", "GetCnt1" and "GetCnt2" resp.
+// To start from scratch one should select : cnt1=0 and cnt2=0
+
+// Reset the area function
+ fNa=0;
+ fXa=0;
+ fYa=0;
+ fIbins=0;
+
+// Clipping parameter to prevent overflow of the counting system
+ fClip=1000000;
+
+// Set the lags for the Fibonacci sequence of the first part
+// The sequence is set to F(97,33,*) (see article)
+ fI=97;
+ fJ=33;
+
+// Unpack the seed value and determine i, j, k and l
+ fSeed=seed;
+ Int_t i,j,k,l;
+ Unpack(seed,i,j,k,l);
+
+// Reset counters
+ fCnt1=0;
+ fCnt2=0;
+
+// Fill the starting table for the first part of the combination
+ Float_t s,t;
+ Int_t m;
+ for (Int_t ii=0; ii<97; ii++)
+ {
+ s=0.;
+ t=0.5;
+
+ for (Int_t jj=1; jj<=24; jj++)
+ {
+ m=(((i*j)%179)*k)%179;
+ i=j;
+ j=k;
+ k=m;
+ l=((53*l)+1)%169;
+ if ((l*m)%64 >= 32) s+=t;
+ t=0.5*t;
+ }
+ fU[ii]=s;
+ }
+
+// Initialise the second part of the combination
+ fC=362436./16777216.;
+ fCd=7654321./16777216.;
+ fCm=16777213./16777216.;
+
+// Generate random numbers upto the user selected starting point
+// on basis of the counting system
+ if (cnt1 > 0) Uniform(cnt1);
+ if (cnt2 > 0)
+ {
+ for (Int_t n=1; n<=cnt2; n++)
+ {
+ Uniform(fClip);
+ }
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliRandom::Unpack(Int_t seed,Int_t& i,Int_t& j,Int_t& k,Int_t& l)
+{
+// Unpack the seed into the four startup parameters i,j,k and l
+//
+// The range of the seed is : 0 <= seed <= 921350143
+//
+// From the article :
+// The i,j and k values may be chosen in the interval : 1 <= n <= 178
+// The l value may be in the interval : 0 <= l <= 168
+//
+// Taking into account the period of the 3-lagged Fibonacci sequence
+// The following "bad" combinations have to be ruled out :
+//
+// i j k l period
+// 1 1 1 X 1
+// 178 1 1 X 4
+// 1 178 1 X 2
+// 1 1 178 X 4
+// 178 178 1 X 4
+// 178 1 178 X 2
+// 1 178 178 X 4
+// 178 178 178 X 1
+//
+// To rule these "bad" combinations out all together, we choose
+// the following allowed ranges :
+// The i,j and k values may be chosen in the interval : 2 <= n <= 177
+// The l value may be in the interval : 0 <= l <= 168
+//
+// and use the formula :
+// seed = (i-2)*176*176*169 + (j-2)*176*169 + (k-2)*169 + l
+
+ if ((seed >= 0) && (seed <= 921350143)) // Check seed value
+ {
+ Int_t idum=seed;
+ Int_t imin2=idum/(176*176*169);
+ idum=idum%(176*176*169);
+ Int_t jmin2=idum/(176*169);
+ idum=idum%(176*169);
+ Int_t kmin2=idum/169;
+
+ i=imin2+2;
+ j=jmin2+2;
+ k=kmin2+2;
+ l=seed%169;
+ }
+ else
+ {
+ cout << " *AliRandom::unpack()* Unallowed seed value encountered."
+ << " seed = " << seed << endl;
+ cout << " Seed will be set to default value." << endl;
+
+ seed=53310452; // Default seed
+ Start(seed,0,0); // Start the sequence for this seed from scratch
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliRandom::GetSeed()
+{
+// Provide the current seed value
+ return fSeed;
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliRandom::GetCnt1()
+{
+// Provide the current value of the counter cnt1
+ return fCnt1;
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliRandom::GetCnt2()
+{
+// Provide the current value of the counter cnt2
+ return fCnt2;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliRandom::Info()
+{
+// Print the current seed, cnt1 and cnt2 values
+ cout << " *Random* seed = " << fSeed
+ << " cnt1 = " << fCnt1 << " cnt2 = " << fCnt2 << endl;
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliRandom::Uniform()
+{
+// Generate uniform random numbers in the interval <0,1>
+//
+// The algorithm is based on lagged Fibonacci sequences (first part)
+// combined with a congruential method (second part)
+// as described in the report by G.Marsaglia and A.Zaman
+// (FSU-SCRI-87-50 Florida State University 1987).
+//
+// Features :
+// 1) Period = 2**144
+// 2) Same sequence of 24-bit real numbers on all common machines
+
+// First part of the combination : F(97,33,*) (see article for explanation)
+ Float_t unirnu=fU[fI-1]-fU[fJ-1];
+ if (unirnu < 0) unirnu+=1.;
+ fU[fI-1]=unirnu;
+ fI-=1;
+ if (fI == 0) fI=97;
+ fJ-=1;
+ if (fJ == 0) fJ=97;
+
+// Second part of the combination (see article for explanation)
+ fC-=fCd;
+ if (fC < 0.) fC+=fCm;
+ unirnu-=fC;
+ if (unirnu < 0.) unirnu+=1.;
+
+// Update the counting system to enable sequence continuation
+// at an arbitrary starting position.
+// Two counters have been introduced to avoid overflow
+// fCnt1 : Counter which goes up to fClip
+// and is reset when fClip is reached
+// fCnt2 : Counts the number of times fClip has been reached
+ fCnt1+=1;
+ if (fCnt1 >= fClip)
+ {
+ fCnt1=0;
+ fCnt2+=1;
+ }
+
+ if (unirnu <= 0.) unirnu=Uniform(); // Exclude 0. from the range
+
+ return unirnu;
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliRandom::Uniform(Float_t a,Float_t b)
+{
+// Generate uniform random numbers in the interval <a,b>
+ Float_t rmin=a;
+ if (a > b) rmin=b;
+
+ Float_t rndm=Uniform();
+ rndm=rmin+fabs(rndm*(a-b));
+
+ return rndm;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliRandom::Uniform(Float_t* vec,Int_t n,Float_t a,Float_t b)
+{
+// Generate a vector of uniform random numbers in the interval <a,b>
+// This saves lots of (member)function calls in case many random
+// numbers are needed at once.
+//
+// n = The number of random numbers to be generated
+//
+// The algorithm is based on lagged Fibonacci sequences (first part)
+// combined with a congruential method (second part)
+// as described in the report by G.Marsaglia and A.Zaman
+// (FSU-SCRI-87-50 Florida State University 1987).
+//
+// Features :
+// 1) Period = 2**144
+// 2) Same sequence of 24-bit real numbers on all common machines
+
+// Determine the minimum of a and b
+ Float_t rmin=a;
+ if (a > b) rmin=b;
+
+// First generate random numbers within <0,1>
+ if (n > 0) // Check n value
+ {
+ for (Int_t jvec=0; jvec<n; jvec++)
+ {
+ // First part of the combination : F(97,33,*)
+ Float_t unirnu=fU[fI-1]-fU[fJ-1];
+ if (unirnu < 0) unirnu+=1.;
+ fU[fI-1]=unirnu;
+ fI-=1;
+ if (fI == 0) fI=97;
+ fJ-=1;
+ if (fJ == 0) fJ=97;
+
+ // Second part of the combination
+ fC-=fCd;
+ if (fC < 0.) fC+=fCm;
+ unirnu-=fC;
+ if (unirnu < 0.) unirnu+=1.;
+
+ // Update the counting system to enable sequence continuation
+ // at an arbitrary starting position.
+ // Two counters have been introduced to avoid overflow
+ // fCnt1 : Counter which goes up to fClip
+ // and is reset when fClip is reached
+ // fCnt2 : Counts the number of times fClip has been reached
+ fCnt1+=1;
+ if (fCnt1 >= fClip)
+ {
+ fCnt1=0;
+ fCnt2+=1;
+ }
+
+ if (unirnu <= 0.) unirnu=Uniform(); // Exclude 0. from the range
+
+ // Fill the vector within the selected range
+ vec[jvec]=rmin+fabs(unirnu*(a-b));
+ }
+ }
+ else
+ {
+ cout << " *AliRandom::Uniform* Invalid value n = " << n << endl;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliRandom::Uniform(Float_t* vec,Int_t n)
+{
+// Generate a vector of uniform random numbers in the interval <0,1>
+// This saves lots of (member)function calls in case many random
+// numbers are needed at once.
+//
+// n = The number of random numbers to be generated
+
+ Uniform(vec,n,0.,1.);
+}
+///////////////////////////////////////////////////////////////////////////
+void AliRandom::Uniform(Int_t n)
+{
+// Generate n uniform random numbers in in one go.
+// This saves lots of (member)function calls in case one needs to skip
+// to a certain point in a sequence.
+//
+// n = The number of random numbers to be generated
+//
+// Note : No check is made here to exclude 0 from the range.
+// It's only the number of generated randoms that counts
+//
+// The algorithm is based on lagged Fibonacci sequences (first part)
+// combined with a congruential method (second part)
+// as described in the report by G.Marsaglia and A.Zaman
+// (FSU-SCRI-87-50 Florida State University 1987).
+//
+// Features :
+// 1) Period = 2**144
+// 2) Same sequence of 24-bit real numbers on all common machines
+
+ if (n > 0) // Check n value
+ {
+ for (Int_t jvec=0; jvec<n; jvec++)
+ {
+ // First part of the combination : F(97,33,*)
+ Float_t unirnu=fU[fI-1]-fU[fJ-1];
+ if (unirnu < 0) unirnu+=1.;
+ fU[fI-1]=unirnu;
+ fI-=1;
+ if (fI == 0) fI=97;
+ fJ-=1;
+ if (fJ == 0) fJ=97;
+
+ // Second part of the combination
+ fC-=fCd;
+ if (fC < 0.) fC+=fCm;
+ unirnu-=fC;
+ if (unirnu < 0.) unirnu+=1.;
+
+ // Update the counting system to enable sequence continuation
+ // at an arbitrary starting position.
+ // Two counters have been introduced to avoid overflow
+ // fCnt1 : Counter which goes up to fClip
+ // and is reset when fClip is reached
+ // fCnt2 : Counts the number of times fClip has been reached
+ fCnt1+=1;
+ if (fCnt1 >= fClip)
+ {
+ fCnt1=0;
+ fCnt2+=1;
+ }
+ }
+ }
+ else
+ {
+ cout << " *AliRandom::Uniform* Invalid value n = " << n << endl;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliRandom::Gauss(Float_t mean,Float_t sigma)
+{
+// Generate gaussian distributed random numbers with certain mean and sigma
+//
+// Method :
+// P(x) = The gaussian distribution function
+// --> ln(P) provides an expression for (x-xmean)**2 from which
+// the following expression for x can be obtained
+//
+// x = xmean +/- sigma * sqrt(-2*ln(q))
+//
+// in which q is an expression in terms of pi, sigma and p and lies within
+// the interval <0,1>.
+//
+// The gaussian distributed x values are obtained as follows :
+//
+// 1) Two uniform random numbers q1 and q2 in <0,1> are generated.
+// 2) q1 is in fact a uniform generated value for P which is substituted
+// directly in the formula above.
+// 3) The value of q2 determines whether we use the + or - sign.
+
+// Generate the two uniform random numbers q1 and q2 in <0,1>
+ Float_t q1,q2;
+ q1=Uniform();
+ q2=Uniform();
+
+// Use q1 and q2 to get the gaussian distributed random number
+ Float_t pi=acos(-1.);
+ Float_t unirng=mean+cos(2.*pi*q2)*sigma*sqrt(-2.*log(q1));
+
+ return unirng;
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliRandom::Gauss()
+{
+// Generate gaussian distributed random numbers with mean=0 and sigma=1
+
+ return Gauss(0.,1.);
+}
+///////////////////////////////////////////////////////////////////////////
+void AliRandom::Gauss(Float_t* vec,Int_t n,Float_t mean,Float_t sigma)
+{
+// Generate a vector of gaussian random numbers with certain mean and sigma
+// This saves lots of (member)function calls in case many random
+// numbers are needed at once.
+//
+// n = The number of random numbers to be generated
+
+ if (n > 0) // Check n value
+ {
+ // The vector to hold the q1 and q2 values.
+ // Two subsequent q[] values are used for q1 and q2
+ // in order to obtain identical random numbers in the vector
+ // as when generating n single ones.
+ Int_t m=2*n;
+ Float_t* q=new Float_t[m];
+ Uniform(q,m);
+
+ // Fill the vector with gaussian random numbers
+ Float_t pi=acos(-1.);
+ Float_t q1,q2;
+ for (Int_t jvec=0; jvec<n; jvec++)
+ {
+ q1=q[jvec*2]; // use two subsequent q[] values
+ q2=q[(jvec*2)+1];
+ vec[jvec]=mean+cos(2.*pi*q2)*sigma*sqrt(-2.*log(q1));
+ }
+ delete [] q;
+ }
+ else
+ {
+ cout << " *AliRandom::Gauss* Invalid value n = " << n << endl;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliRandom::Gauss(Float_t* vec,Int_t n)
+{
+// Generate a vector of gaussian random numbers with mean=0 and sigma=1
+// This saves lots of (member)function calls in case many random
+// numbers are needed at once.
+//
+// n = The number of random numbers to be generated
+
+ Gauss(vec,n,0.,1.);
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliRandom::Poisson(Float_t mean)
+{
+// Generate Poisson distributed random numbers with certain mean
+//
+// Method :
+//
+// P(n) = exp(-mean)*mean**n/n! is the Poisson distribution function
+//
+// with : n = 0,1,2,3,... and mean > 0
+//
+// To generate the distribution, the "sum trick" is used as mentioned
+// in "Formulae and Methods in Experimental data Evaluation Vol. 1"
+
+ Float_t unirnp=0.; // Initialise the random number value
+
+ if (mean <= 0.) // Check the mean value
+ {
+ cout << " *AliRandom::Poisson* Invalid value mean = " << mean
+ << " Should be positive." << endl;
+ }
+ else
+ {
+ if (mean > 80.) // Use gaussian dist. for high mean values to save time
+ {
+ Float_t grndm=Gauss();
+ Float_t rpois=mean+grndm*sqrt(mean);
+ Int_t npois=int(rpois);
+ if ((rpois-float(npois)) > 0.5) npois++;
+ unirnp=float(npois);
+ }
+ else // Construct a Poisson random number from a uniform one
+ {
+ Int_t npois=-1;
+ Float_t expxm=exp(-mean);
+ Float_t poitst=1.;
+ while (poitst > expxm)
+ {
+ Float_t rndm=Uniform();
+ npois++;
+ poitst=poitst*rndm;
+ }
+ unirnp=float(npois);
+ } // end of check for using Gauss method
+ } // end of mean validity checkn
+ return unirnp;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliRandom::Poisson(Float_t* vec,Int_t n,Float_t mean)
+{
+// Generate a vector of Poisson dist. random numbers with certain mean
+// This saves lots of (member)function calls in case many random
+// numbers are needed at once.
+//
+// n = The number of random numbers to be generated
+//
+// Method :
+//
+// P(n) = exp(-mean)*mean**n/n! is the Poisson distribution function
+//
+// with : n = 0,1,2,3,... and mean > 0
+//
+// To generate the distribution, the "sum trick" is used as mentioned
+// in "Formulae and Methods in Experimental data Evaluation Vol. 1"
+
+ if (n <= 0) // Check n value
+ {
+ cout << " *AliRandom::Poisson* Invalid value n = " << n << endl;
+ }
+ else
+ {
+ if (mean <= 0.) // Check the mean value
+ {
+ cout << " *AliRandom::Poisson* Invalid value mean = " << mean
+ << " Should be positive." << endl;
+ }
+ else
+ {
+ if (mean > 80.) // Use gaussian dist. for high mean values to save time
+ {
+ Float_t* grndm=new Float_t[n];
+ Gauss(grndm,n);
+ Int_t npois;
+ Float_t rpois;
+ for (Int_t jvec=0; jvec<n; jvec++)
+ {
+ rpois=mean+grndm[jvec]*sqrt(mean);
+ npois=int(rpois);
+ if ((rpois-float(npois)) > 0.5) npois++;
+ vec[jvec]=float(npois);
+ }
+ delete [] grndm;
+ }
+ else // Construct Poisson random numbers from a uniform ones
+ {
+ Float_t expxm=exp(-mean);
+ Int_t npois;
+ Float_t poitst;
+ for (Int_t jvec=0; jvec<n; jvec++)
+ {
+ npois=-1;
+ poitst=1.;
+ while (poitst > expxm)
+ {
+ Float_t rndm=Uniform();
+ npois++;
+ poitst=poitst*rndm;
+ }
+ vec[jvec]=float(npois);
+ }
+ } // end of check for using Gauss method
+ } // end of mean validity check
+ } // end of n validity check
+}
+///////////////////////////////////////////////////////////////////////////
+void AliRandom::SetUser(Float_t a,Float_t b,Int_t n,Float_t (*f)(Float_t))
+{
+// Determine the area under f(x) as a function of x
+// This is called the "area function" and serves as a basis to
+// provide random numbers in [a,b] according to the user defined
+// distribution f(x).
+// The area function is normalised such that the most extreme
+// value is 1 or -1.
+
+ fNa=n+1; // The number of bins for the area function
+ fXa=new Float_t[fNa]; // The binned x values of the area function
+ fYa=new Float_t[fNa]; // The corresponding summed f(x) values
+ fIbins=new Int_t[fNa]; // The bin numbers of the random x candidates
+
+ Float_t xmin=a;
+ if (a > b) xmin=b;
+ Float_t step=fabs(a-b)/float(n);
+
+ Float_t x;
+ Float_t extreme=0;
+ for (Int_t i=0; i<fNa; i++) // Fill bins of the area function
+ {
+ x=xmin+float(i)*step;
+ fXa[i]=x;
+ fYa[i]=f(x);
+ if (i > 0) fYa[i]+=fYa[i-1];
+ if (fabs(fYa[i]) > extreme) extreme=fabs(fYa[i]);
+ }
+ fYamin=fYa[0]/extreme;
+ fYamax=fYa[0]/extreme;
+ for (Int_t j=0; j<fNa; j++) // Normalise the area function
+ {
+ fYa[j]=fYa[j]/extreme;
+ if (fYa[j] < fYamin) fYamin=fYa[j];
+ if (fYa[j] > fYamax) fYamax=fYa[j];
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliRandom::SetUser(Float_t* x,Float_t* y,Int_t n)
+{
+// Determine the area under y[i] as a function of x[i]
+// This is called the "area function" and serves as a basis to
+// provide random numbers in x according to the user provided
+// distribution (x[i],y[i]).
+// The area function is normalised such that the most extreme
+// value is 1 or -1.
+
+ fNa=n; // The number of bins for the area function
+ fXa=new Float_t[fNa]; // The binned x values of the area function
+ fYa=new Float_t[fNa]; // The corresponding summed y values
+ fIbins=new Int_t[fNa]; // The bin numbers of the random x candidates
+
+// Order input data with increasing x
+ fXa[0]=x[0];
+ fYa[0]=y[0];
+ for (Int_t i=1; i<fNa; i++) // Loop over x[]
+ {
+ for (Int_t j=0; j<i; j++) // Loop over xa[]
+ {
+ if (x[i] < fXa[j])
+ {
+ for (Int_t k=i; k>=j; k--) // Create empty position
+ {
+ fXa[k+1]=fXa[k];
+ fYa[k+1]=fYa[k];
+ }
+ fXa[j]=x[i]; // Put x[] value in empty position
+ fYa[j]=y[i]; // Put y[] value in empty position
+ break; // Go for next x[] value
+ }
+ if (j == i-1) // This x[] value is the largest so far
+ {
+ fXa[i]=x[i]; // Put x[] value at the end of x[]
+ fYa[i]=y[i]; // Put y[] value at the end of y[]
+ }
+ } // End loop over fXa[]
+ } // End loop over x[]
+
+ Float_t extreme=0;
+ for (Int_t l=0; l<fNa; l++) // Fill area function
+ {
+ if (l > 0) fYa[l]+=fYa[l-1];
+ if (fabs(fYa[l]) > extreme) extreme=fabs(fYa[l]);
+ }
+ fYamin=fYa[0]/extreme;
+ fYamax=fYa[0]/extreme;
+ for (Int_t m=0; m<fNa; m++) // Normalise the area function
+ {
+ fYa[m]=fYa[m]/extreme;
+ if (fYa[m] < fYamin) fYamin=fYa[m];
+ if (fYa[m] > fYamax) fYamax=fYa[m];
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliRandom::User()
+{
+// Provide a random number according to the user defined distribution
+//
+// Method :
+// --------
+// Select by a uniform random number a certain area fraction (from fYa[])
+// of the area function.
+// The required random number is given by the corresponding x value (fXa[])
+// of the area function.
+// In case of more than one x value candidate, select randomly one of them.
+
+ Float_t unirnf=0;
+
+ Float_t ra=Uniform(fYamin,fYamax); // Random value of the area function
+ Float_t dmin=100.*fabs(fYamax-fYamin); // Init. the min. distance encountered
+ Int_t ncand=0;
+ Float_t dist;
+ for (Int_t i=0; i<fNa; i++) // Search for fYa[] value(s) closest to ra
+ {
+ dist=fabs(ra-fYa[i]);
+ if (dist <= dmin) // fYa[i] within smallest distance --> extra candidate
+ {
+ ncand++;
+ if (dist < dmin) ncand=1; // Smallest distance so far --> THE candidate
+ dmin=dist;
+ fIbins[ncand-1]=i+1;
+ }
+ }
+
+ Int_t jbin=0; // The bin number to hold the required x value
+ if (ncand == 1) jbin=fIbins[0];
+
+ if (ncand > 1) // Multiple x value candidates --> pick one randomly
+ {
+ Float_t cand=Uniform(1.,float(ncand));
+ Int_t jcand=int(cand);
+ if ((cand-float(jcand)) > 0.5) jcand++;
+ jbin=fIbins[jcand-1];
+ }
+
+ if (jbin > 0) // Pick randomly a value in this x-bin
+ {
+ Float_t xlow=fXa[jbin-1];
+ if (jbin > 1) xlow=fXa[jbin-2];
+ Float_t xup=fXa[jbin-1];
+ if (jbin < fNa-1) xup=fXa[jbin];
+ unirnf=Uniform(xlow,xup);
+ }
+
+ if (ncand == 0) cout << " *AliRandom::User* No candidate found." << endl;
+
+ return unirnf;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliRandom::User(Float_t* vec,Int_t n)
+{
+// Generate a vector of random numbers according to a user defined dist.
+// This saves lots of (member)function calls in case many random
+// numbers are needed at once.
+//
+// n = The number of random numbers to be generated
+//
+// Method :
+// --------
+// Select by a uniform random number a certain area fraction (from fYa[])
+// of the area function.
+// The required random number is given by the corresponding x value (fXa[])
+// of the area function.
+// In case of more than one x value candidate, select randomly one of them.
+
+ Float_t unirnf,ra,dmin,dist;
+ Int_t ncand,jbin;
+ for (Int_t jvec=0; jvec<n; jvec++)
+ {
+ unirnf=0;
+ ra=Uniform(fYamin,fYamax); // Random value of the area function
+ dmin=100.*fabs(fYamax-fYamin); // Init. the min. distance encountered
+ ncand=0;
+ for (Int_t i=0; i<fNa; i++) // Search for fYa[] value(s) closest to ra
+ {
+ dist=fabs(ra-fYa[i]);
+ if (dist <= dmin) // fYa[i] within smallest distance --> extra candidate
+ {
+ ncand++;
+ if (dist < dmin) ncand=1; // Smallest distance so far --> THE candidate
+ dmin=dist;
+ fIbins[ncand-1]=i+1;
+ }
+ }
+
+ jbin=0; // The bin number to hold the required x value
+ if (ncand == 1) jbin=fIbins[0];
+
+ if (ncand > 1) // Multiple x value candidates --> pick one randomly
+ {
+ Float_t cand=Uniform(1.,float(ncand));
+ Int_t jcand=int(cand);
+ if ((cand-float(jcand)) > 0.5) jcand++;
+ jbin=fIbins[jcand-1];
+ }
+
+ if (jbin > 0) // Pick randomly a value in this x-bin
+ {
+ Float_t xlow=fXa[jbin-1];
+ if (jbin > 1) xlow=fXa[jbin-2];
+ Float_t xup=fXa[jbin-1];
+ if (jbin < fNa-1) xup=fXa[jbin];
+ unirnf=Uniform(xlow,xup);
+ }
+
+ if (ncand == 0) cout << " *AliRandom::User* No candidate found." << endl;
+
+ vec[jvec]=unirnf;
+
+ }
+}
+///////////////////////////////////////////////////////////////////////////
--- /dev/null
+#ifndef ALIRANDOM_H
+#define ALIRANDOM_H
+///////////////////////////////////////////////////////////////////////////
+// Class AliRandom
+// Generate universal random numbers on all common machines.
+// Available distributions : Uniform, Gaussian, Poisson and
+// User defined function
+//
+// Features :
+// ----------
+// 1) Period = 2**144
+// 2) Same sequence of 24-bit real numbers on all common machines
+//
+// Reference :
+// -----------
+// G.Marsaglia and A.Zaman, FSU-SCRI-87-50, Florida State University, 1987.
+//
+// Coding example :
+// ----------------
+//
+// Float_t rndm; // Variable to hold a single random number
+// const Int_t n=1000;
+// Float_t rvec[n]; // Vector to hold n random numbers
+//
+// AliRandom r; // Create a Random object with default sequence
+//
+// rndm=r.Uniform(); // Provide a uniform random number in <0,1>
+// Float_t a=3.;
+// Float_t b=5.;
+// rndm=r.Uniform(a,b); // Provide a uniform random number in <a,b>
+// r.Uniform(rvec,n); // Provide n uniform randoms in <0,1> in rvec
+// r.Uniform(rvec,n,a,b); // Provide n uniform randoms in <a,b> in rvec
+//
+// rndm=r.Gauss(); // Provide a Gaussian random number with
+// // mean=0 and sigma=1
+// Float_t mean=25.;
+// Float_t sigma=5.;
+// rndm=r.Gauss(mean,sigma); // Provide a Gaussian random number
+// // with specified mean and sigma
+// r.Gauss(rvec,n); // n Gaussian randoms mean=0 sigma=1
+// r.Gauss(rvec,n,mean,sigma); // n Gaussian randoms with specified
+// // mean and sigma
+//
+// rndm=r.Poisson(mean); // Provide a Poisson random number with
+// // specified mean
+// r.Poisson(rvec,nmean); // n Poisson randoms with specified mean
+//
+// Int_t seed=1837724
+// AliRandom p(seed); // Create a Random object with specified seed.
+// // The sequence is started from scratch.
+// Int_t cnt1=25;
+// Int_t cnt2=8;
+// AliRandom q(seed,cnt1,cnt2); // Create a Random object with specified seed
+// // The sequence is started at the location
+// // denoted by the counters cnt1 and cnt2.
+//
+// q.Info(); // Print the current seed, cnt1 and cnt2 values.
+// q.GetSeed(); // Provide the current seed value.
+// q.GetCnt1(); // Provide the current cnt1 value.
+// q.GetCnt2(); // Provide the current cnt2 value.
+//
+// Float_t udist(Float_t x) // A user defined distribution
+// {
+// return x*x-4.*x;
+// }
+//
+// Int_t nbins=100;
+// q.SetUser(a,b,nbins,udist); // Initialise generator for udist distribution
+// q.User(); // Provide a random number according to the udist distribution
+// q.User(rvec,n); // Provide n randoms according to the udist distribution
+//
+// Float_t* x=new Float_t[nbins];
+// Float_t* y=new Float_t[nbins];
+//
+// ... code to fill x[] and y[] ..
+//
+// AliRandom s;
+// s.SetUser(x,y,nbins); // Initialise generator for (x[i],y[i]) distribution
+// s.User(); // Provide a random number according to the user distribution
+// s.User(rvec,n); // Provide n randoms according to the user distribution
+//
+// Notes :
+// -------
+// 1) Allowed seed values : 0 <= seed <= 921350143
+// Default seed = 53310452
+// 2) To ensure a unique sequence for each run, one can automatically
+// construct a seed value by e.g. using the date and time.
+// 3) Using the rvec facility saves a lot of CPU time for large n values.
+//
+//--- NvE 11-oct-1997 UU-SAP Utrecht
+///////////////////////////////////////////////////////////////////////////
+
+#include <iostream.h>
+#include <math.h>
+
+#include "TObject.h"
+
+class AliRandom : public TObject
+{
+ public:
+ AliRandom(); // Constructor with default sequence
+ AliRandom(Int_t seed); // Constructor with user defined seed
+ AliRandom(Int_t seed,Int_t cnt1,Int_t cnt2); // User defined starting point
+ ~AliRandom(); // Destructor
+ Int_t GetSeed(); // Provide current seed value
+ Int_t GetCnt1(); // Provide current counter value cnt1
+ Int_t GetCnt2(); // Provide current counter value cnt2
+ void Info(); // Print current seed, cnt1 and cnt2
+ Float_t Uniform(); // Uniform dist. within <0,1>
+ Float_t Uniform(Float_t a,Float_t b); // Uniform dist. within <a,b>
+ void Uniform(Float_t* vec,Int_t n); // n uniform randoms in <0,1>
+ void Uniform(Float_t* vec,Int_t n,Float_t a,Float_t b); // see above
+ Float_t Gauss(); // Gaussian dist. with mean=0 sigma=1
+ Float_t Gauss(Float_t mean,Float_t sigma); // Gaussian dist. with mean and sigma
+ void Gauss(Float_t* vec,Int_t n); // n Gaussian randoms mean=0 sigma=1
+ void Gauss(Float_t* vec,Int_t n,Float_t mean,Float_t sigma); // see above
+ Float_t Poisson(Float_t mean); // Poisson dist. with certain mean
+ void Poisson(Float_t* vec,Int_t n,Float_t mean); // n Poisson randoms with mean
+ void SetUser(Float_t a,Float_t b,Int_t n,Float_t (*f)(Float_t)); // User dist. f(x)
+ void SetUser(Float_t* x,Float_t* y,Int_t n); // User dist. arrays
+ Float_t User(); // Provide random in [a,b] according to user distribution
+ void User(Float_t* vec,Int_t n); // n randoms in [a,b] from user dist.
+
+ private:
+ Int_t fI,fJ,fSeed,fCnt1,fCnt2,fClip; // Indices, seed and counters
+ Float_t fU[97],fC,fCd,fCm; // The Fibonacci parameters
+ void Start(Int_t seed,Int_t cnt1,Int_t cnt2); // Start at certain point
+ void Unpack(Int_t seed,Int_t& i,Int_t& j,Int_t& k,Int_t& l); // Unpack the seed
+ void Uniform(Int_t n); // n uniform randoms for quick skipping
+ Int_t fNa; //! The number of bins of the area function
+ Float_t* fXa; //! The binned x values of the area function
+ Float_t* fYa; //! The corresponding y values of the area function
+ Float_t fYamin,fYamax; //! The min. and max. y values of the area function
+ Int_t* fIbins; //! The bin numbers of the random x candidates
+
+ ClassDef(AliRandom,1) // Class definition to enable ROOT I/O
+};
+#endif
--- /dev/null
+#include "AliSample.h"
+
+AliSample::AliSample()
+{
+// Creation of an Aliample object and resetting the statistics values
+// The dimension is initialised to maximum
+ fDim=fMaxdim;
+ fNames[0]='X';
+ fNames[1]='Y';
+ fNames[2]='Z';
+ fN=0;
+ Reset();
+}
+///////////////////////////////////////////////////////////////////////////
+AliSample::~AliSample()
+{
+// Default destructor
+}
+///////////////////////////////////////////////////////////////////////////
+void AliSample::Reset()
+{
+// Resetting the statistics values for a certain Sample object
+// Dimension is NOT changed
+ fN=0;
+ for (Int_t i=0; i<fDim; i++)
+ {
+ fSum[i]=0.;
+ fMean[i]=0.;
+ fVar[i]=0.;
+ fSigma[i]=0.;
+ for (Int_t j=0; j<fDim; j++)
+ {
+ fSum2[i][j]=0.;
+ fCov[i][j]=0.;
+ fCor[i][j]=0.;
+ }
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliSample::Enter(Float_t x)
+{
+// Entering a value into a 1-dim. sample
+// In case of first entry the dimension is set to 1
+ if (fN == 0)
+ {
+ fDim=1;
+ fNames[0]='X';
+ fNames[1]='-';
+ fNames[2]='-';
+ }
+ if (fDim != 1)
+ {
+ cout << " *AliSample::enter* Error : Not a 1-dim sample." << endl;
+ }
+ else
+ {
+ fN+=1;
+ fSum[0]+=x;
+ fSum2[0][0]+=x*x;
+ Compute();
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliSample::Remove(Float_t x)
+{
+// Removing a value from a 1-dim. sample
+ if (fDim != 1)
+ {
+ cout << " *AliSample::remove* Error : Not a 1-dim sample." << endl;
+ }
+ else
+ {
+ fN-=1;
+ fSum[0]-=x;
+ fSum2[0][0]-=x*x;
+ Compute();
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliSample::Enter(Float_t x,Float_t y)
+{
+// Entering a pair (x,y) into a 2-dim. sample
+// In case of first entry the dimension is set to 2
+ if (fN == 0)
+ {
+ fDim=2;
+ fNames[0]='X';
+ fNames[1]='Y';
+ fNames[2]='-';
+ }
+ if (fDim != 2)
+ {
+ cout << " *AliSample::enter* Error : Not a 2-dim sample." << endl;
+ }
+ else
+ {
+ fN+=1;
+ fSum[0]+=x;
+ fSum[1]+=y;
+ fSum2[0][0]+=x*x;
+ fSum2[0][1]+=x*y;
+ fSum2[1][0]+=y*x;
+ fSum2[1][1]+=y*y;
+ Compute();
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliSample::Remove(Float_t x,Float_t y)
+{
+// Removing a pair (x,y) from a 2-dim. sample
+ if (fDim != 2)
+ {
+ cout << " *AliSample::remove* Error : Not a 2-dim sample." << endl;
+ }
+ else
+ {
+ fN-=1;
+ fSum[0]-=x;
+ fSum[1]-=y;
+ fSum2[0][0]-=x*x;
+ fSum2[0][1]-=x*y;
+ fSum2[1][0]-=y*x;
+ fSum2[1][1]-=y*y;
+ Compute();
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliSample::Enter(Float_t x,Float_t y,Float_t z)
+{
+// Entering a set (x,y,z) into a 3-dim. sample
+// In case of first entry the dimension is set to 3
+ if (fN == 0)
+ {
+ fDim=3;
+ fNames[0]='X';
+ fNames[1]='Y';
+ fNames[2]='Z';
+ }
+ if (fDim != 3)
+ {
+ cout << " *AliSample::enter* Error : Not a 3-dim sample." << endl;
+ }
+ else
+ {
+ fN+=1;
+ fSum[0]+=x;
+ fSum[1]+=y;
+ fSum[2]+=z;
+ fSum2[0][0]+=x*x;
+ fSum2[0][1]+=x*y;
+ fSum2[0][2]+=x*z;
+ fSum2[1][0]+=y*x;
+ fSum2[1][1]+=y*y;
+ fSum2[1][2]+=y*z;
+ fSum2[2][0]+=z*x;
+ fSum2[2][1]+=z*y;
+ fSum2[2][2]+=z*z;
+ Compute();
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliSample::Remove(Float_t x,Float_t y,Float_t z)
+{
+// Removing a set (x,y,z) from a 3-dim. sample
+ if (fDim != 3)
+ {
+ cout << " *AliSample::remove* Error : Not a 3-dim sample." << endl;
+ }
+ else
+ {
+ fN-=1;
+ fSum[0]-=x;
+ fSum[1]-=y;
+ fSum[2]-=z;
+ fSum2[0][0]-=x*x;
+ fSum2[0][1]-=x*y;
+ fSum2[0][2]-=x*z;
+ fSum2[1][0]-=y*x;
+ fSum2[1][1]-=y*y;
+ fSum2[1][2]-=y*z;
+ fSum2[2][0]-=z*x;
+ fSum2[2][1]-=z*y;
+ fSum2[2][2]-=z*z;
+ Compute();
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliSample::Compute()
+{
+// Computation of the various statistical values
+// after each entering or removing action on a certain sample
+ Float_t rn=fN;
+ for (Int_t k=0; k<fDim; k++)
+ {
+ fMean[k]=fSum[k]/rn;
+ fVar[k]=(fSum2[k][k]/rn)-(fMean[k]*fMean[k]);
+ if (fVar[k] < 0.) fVar[k]=0.;
+ fSigma[k]=sqrt(fVar[k]);
+ }
+ for (Int_t i=0; i<fDim; i++)
+ {
+ for (Int_t j=0; j<fDim; j++)
+ {
+ fCov[i][j]=(fSum2[i][j]/rn)-(fMean[i]*fMean[j]);
+ Float_t sigij=fSigma[i]*fSigma[j];
+ if (sigij != 0.) fCor[i][j]=fCov[i][j]/sigij;
+ }
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliSample::GetDimension()
+{
+// Provide the dimension of a certain sample
+ return fDim;
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliSample::GetN()
+{
+// Provide the number of entries of a certain sample
+ return fN;
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliSample::GetSum(Int_t i)
+{
+// Provide the sum of a certain variable
+ if (fDim < i)
+ {
+ cout << " *AliSample::sum* Error : Dimension less than " << i << endl;
+ return 0.;
+ }
+ else
+ {
+ return fSum[i-1];
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliSample::GetMean(Int_t i)
+{
+// Provide the mean of a certain variable
+ if (fDim < i)
+ {
+ cout << " *AliSample::mean* Error : Dimension less than " << i << endl;
+ return 0.;
+ }
+ else
+ {
+ return fMean[i-1];
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliSample::GetVar(Int_t i)
+{
+// Provide the variance of a certain variable
+ if (fDim < i)
+ {
+ cout << " *AliSample::var* Error : Dimension less than " << i << endl;
+ return 0.;
+ }
+ else
+ {
+ return fVar[i-1];
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliSample::GetSigma(Int_t i)
+{
+// Provide the standard deviation of a certain variable
+ if (fDim < i)
+ {
+ cout << " *AliSample::sigma* Error : Dimension less than " << i << endl;
+ return 0.;
+ }
+ else
+ {
+ return fSigma[i-1];
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliSample::GetCov(Int_t i,Int_t j)
+{
+// Provide the covariance between variables i and j
+ if ((fDim < i) || (fDim < j))
+ {
+ Int_t k=i;
+ if (j > i) k=j;
+ cout << " *AliSample::cov* Error : Dimension less than " << k << endl;
+ return 0.;
+ }
+ else
+ {
+ return fCov[i-1][j-1];
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliSample::GetCor(Int_t i,Int_t j)
+{
+// Provide the correlation between variables i and j
+ if ((fDim < i) || (fDim < j))
+ {
+ Int_t k=i;
+ if (j > i) k=j;
+ cout << " *AliSample::cor* Error : Dimension less than " << k << endl;
+ return 0.;
+ }
+ else
+ {
+ return fCor[i-1][j-1];
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliSample::Info()
+{
+// Printing of statistics of all variables
+ for (Int_t i=0; i<fDim; i++)
+ {
+ cout << " " << fNames[i] << " : N = " << fN;
+ cout << " Sum = " << fSum[i] << " Mean = " << fMean[i];
+ cout << " Var = " << fVar[i] << " Sigma = " << fSigma[i] << endl;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliSample::Info(Int_t i)
+{
+// Printing of statistics of ith variable
+ if (fDim < i)
+ {
+ cout << " *AliSample::Info(i)* Error : Dimension less than " << i << endl;
+ }
+ else
+ {
+ cout << " " << fNames[i-1] << " : N = " << fN;
+ cout << " Sum = " << fSum[i-1] << " Mean = " << fMean[i-1];
+ cout << " Var = " << fVar[i-1] << " Sigma = " << fSigma[i-1] << endl;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliSample::Info(Int_t i,Int_t j)
+{
+// Printing of covariance and correlation between variables i and j
+ if ((fDim < i) || (fDim < j))
+ {
+ Int_t k=i;
+ if (j > i) k=j;
+ cout << " *AliSample::Info(i,j)* Error : Dimension less than " << k << endl;
+ }
+ else
+ {
+ cout << " " << fNames[i-1] << "-" << fNames[j-1] << " correlation :";
+ cout << " Cov. = " << fCov[i-1][j-1] << " Cor. = " << fCor[i-1][j-1] << endl;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
--- /dev/null
+#ifndef ALISAMPLE_H
+#define ALISAMPLE_H
+///////////////////////////////////////////////////////////////////////////
+// Class AliSample
+// Perform statistics on various multi-dimensional data samples
+// A data sample can be filled using the "Enter" and/or "Remove" functions,
+// whereas the "Reset" function resets the complete sample to 'empty'.
+// The info which can be extracted from a certain data sample are the
+// sum, mean, variance, sigma, covariance and correlation.
+// The "Info" function provides all statistics data for a certain sample.
+// The variables for which these stat. parameters have to be calculated
+// are indicated by the index of the variable which is passed as an
+// argument to the various member functions.
+// The index convention for a data point (x,y) is : x=1 y=2
+//
+// Example :
+// ---------
+// For an AliSample s a data point (x,y) can be entered as s.Enter(x,y) and
+// the mean_x can be obtained as s.GetMean(1) whereas the mean_y is obtained
+// via s.GetMean(2).
+// The correlation between x and y is available via s.GetCor(1,2).
+// The x-statistics are obtained via s.Info(1), y-statistics via s.Info(2),
+// and the covariance and correlation between x and y via s.Info(1,2).
+// All statistics of a sample are obtained via s.Info().
+//
+//--- NvE 30-mar-1996 CERN Geneva
+///////////////////////////////////////////////////////////////////////////
+
+#include <math.h>
+#include <iostream.h>
+
+#include "Rtypes.h"
+
+class AliSample
+{
+ public:
+ AliSample(); // Default constructor
+ virtual ~AliSample(); // Default destructor
+ void Reset(); // Reset complete statistics
+ void Enter(Float_t x); // Enter value for 1-dim. sample
+ void Remove(Float_t x); // Remove value from 1-dim. sample
+ void Enter(Float_t x, Float_t y); // Enter value for 2-dim. sample
+ void Remove(Float_t x, Float_t y); // Remove value from 2-dim. sample
+ void Enter(Float_t x, Float_t y, Float_t z); // Enter value for 3-dim. sample
+ void Remove(Float_t x, Float_t y, Float_t z); // Remove value from 3-dim. sample
+ Int_t GetDimension(); // Provide dimension of the sample
+ Int_t GetN(); // Provide the number of entries
+ Float_t GetSum(Int_t i); // Provide sum for i-th variable
+ Float_t GetMean(Int_t i); // Provide mean for i-th variable
+ Float_t GetVar(Int_t i); // Provide variance for i-th variable
+ Float_t GetSigma(Int_t i); // Standard deviation for i-th variable
+ Float_t GetCov(Int_t i, Int_t j); // Covariance for i-th and j-th variable
+ Float_t GetCor(Int_t i, Int_t j); // Correlation for i-th and j-th variable
+ void Info(); // Stat. info for the complete sample
+ void Info(Int_t i); // Stat. info for the i-th variable
+ void Info(Int_t i, Int_t j); // Stat. info for i-th and j-th variable
+
+ private:
+ Int_t fDim; // Dimension of the sample
+ Int_t fN; // Number of entries of the sample
+ enum {fMaxdim=3}; // Maximum supported dimension
+ char fNames[fMaxdim]; // Variable names i.e. X,Y,Z
+ Float_t fSum[fMaxdim]; // Total sum for each variable
+ Float_t fSum2[fMaxdim][fMaxdim]; // Total sum**2 for each variable
+ void Compute(); // Compute the various quantities
+ Float_t fMean[fMaxdim]; // Mean for each variable
+ Float_t fVar[fMaxdim]; // Variation for each variable
+ Float_t fSigma[fMaxdim]; // Standard deviation for each variable
+ Float_t fCov[fMaxdim][fMaxdim]; // Covariances of pairs of variables
+ Float_t fCor[fMaxdim][fMaxdim]; // Correlations of pairs of variables
+};
+#endif
--- /dev/null
+#include "AliSignal.h"
+
+ClassImp(AliSignal) // Class implementation to enable ROOT I/O
+
+AliSignal::AliSignal()
+{
+// Creation of an AliSignal object and initialisation of parameters
+ Reset();
+}
+///////////////////////////////////////////////////////////////////////////
+AliSignal::~AliSignal()
+{
+// Destructor to delete dynamically allocated memory
+}
+///////////////////////////////////////////////////////////////////////////
+void AliSignal::Reset()
+{
+// Reset all values
+ Float_t r[3]={0,0,0};
+ SetPosition(r,"sph");
+ fSignal=0;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliSignal::SetSignal(Float_t sig)
+{
+// Store signal value
+ fSignal=sig;
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliSignal::GetSignal()
+{
+// Provide signal value
+ return fSignal;
+}
+///////////////////////////////////////////////////////////////////////////
--- /dev/null
+#ifndef ALISIGNAL_H
+#define ALISIGNAL_H
+///////////////////////////////////////////////////////////////////////////
+// Class AliSignal
+// Handling of ALICE (extrapolated) signals.
+//
+// Note :
+// ------
+// Signal positions (r) and reference frames (f) are specified via
+// SetPosition(r,f) under the following conventions :
+//
+// f="car" ==> r is Cartesian (x,y,z)
+// f="sph" ==> r is Spherical (r,theta,phi)
+// f="cyl" ==> r is Cylindrical (rho,phi,z)
+//
+// All angles are in radians.
+//
+// Example :
+// ---------
+//
+// AliSignal s;
+// Float_t pos[3]={-1,25,7};
+// Float_t signal=120.8;
+// s.SetPosition(pos,"car");
+// s.SetSignal(signal);
+// Float_t loc[3];
+// s.GetPosition(loc,"sph");
+// Float_t adc=s.GetSignal();
+//
+//--- NvE 23-jan-1999 UU-SAP Utrecht
+///////////////////////////////////////////////////////////////////////////
+
+#include "AliPosition.h"
+
+class AliSignal : public TObject,public AliPosition
+{
+ public:
+ AliSignal(); // Default constructor
+ ~AliSignal(); // Destructor
+ virtual void SetSignal(Float_t sig); // Store signal value
+ virtual Float_t GetSignal(); // Provide signal value
+ virtual void Reset(); // Reset all values to 0
+
+ protected:
+ Float_t fSignal; // Signal value
+
+ ClassDef(AliSignal,1) // Class definition to enable ROOT I/O
+};
+#endif
--- /dev/null
+#include "AliTrack.h"
+
+ClassImp(AliTrack) // Class implementation to enable ROOT I/O
+
+AliTrack::AliTrack()
+{
+// Default constructor
+// All variables initialised to 0
+ fDecays=0;
+ Reset();
+}
+///////////////////////////////////////////////////////////////////////////
+AliTrack::~AliTrack()
+{
+// Destructor to delete memory allocated for decay tracks array
+ if (fDecays)
+ {
+ fDecays->Delete();
+ delete fDecays;
+ fDecays=0;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliTrack::Reset()
+{
+// Reset all variables to 0
+ fM=0;
+ fQ=0;
+ fNdec=0;
+ Double_t a[4]={0,0,0,0};
+ SetVector(a,"sph");
+ if (fDecays)
+ {
+ fDecays->Delete();
+ delete fDecays;
+ fDecays=0;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliTrack::Set3Momentum(Ali3Vector& p)
+{
+// Set the track parameters according to the 3-momentum p
+ Double_t E=sqrt(p.Dot(p)+fM*fM);
+ SetVector(E,p);
+}
+///////////////////////////////////////////////////////////////////////////
+void AliTrack::Set4Momentum(Ali4Vector& p)
+{
+// Set the track parameters according to the 4-momentum p
+ Double_t E=p.GetScalar();
+ Ali3Vector pv=p.Get3Vector();
+ SetVector(E,pv);
+
+ Double_t m2=p.Dot(p);
+ fM=0;
+ if (m2 > 0.) fM=sqrt(m2);
+}
+///////////////////////////////////////////////////////////////////////////
+void AliTrack::SetMass(Double_t m)
+{
+// Set the particle mass
+ fM=m;
+ Ali3Vector p=Get3Vector();
+ Double_t E=sqrt(p.Dot(p)+fM*fM);
+ SetVector(E,p);
+}
+///////////////////////////////////////////////////////////////////////////
+void AliTrack::SetCharge(Float_t q)
+{
+// Set the particle charge
+ fQ=q;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliTrack::Info(TString f)
+{
+// Provide track information within the coordinate frame f
+ cout << " *AliTrack::Info* Mass : " << fM << " Charge : " << fQ
+ << " Momentum : " << GetMomentum() << " Ntracks : " << fNdec << endl;
+ cout << " ";
+ Ali4Vector::Info(f);
+}
+///////////////////////////////////////////////////////////////////////////
+void AliTrack::List(TString f)
+{
+// Provide current track and decay level 1 information within coordinate frame f
+
+ Info(f); // Information of the current track
+
+ // Decay products of this track
+ AliTrack* td;
+ for (Int_t id=1; id<=fNdec; id++)
+ {
+ td=GetDecayTrack(id);
+ if (td)
+ {
+ cout << " ---Level 1 sec. track no. " << id << endl;
+ cout << " ";
+ td->Info(f);
+ }
+ else
+ {
+ cout << " *AliTrack::List* Error : No decay track present." << endl;
+ }
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliTrack::ListAll(TString f)
+{
+// Provide complete track and decay information within the coordinate frame f
+
+ Info(f); // Information of the current track
+
+ AliTrack* t=this;
+ Dump(t,1,f); // Information of all decay products
+}
+//////////////////////////////////////////////////////////////////////////
+void AliTrack::Dump(AliTrack* t,Int_t n,TString f)
+{
+// Recursively provide the info of all decay levels of this track
+ AliTrack* td;
+ for (Int_t id=1; id<=t->GetNdecay(); id++)
+ {
+ td=t->GetDecayTrack(id);
+ if (td)
+ {
+ cout << " ---Level " << n << " sec. track no. " << id << endl;
+ cout << " ";
+ td->Info(f);
+
+ // Go for next decay level of this decay track recursively
+ Dump(td,n+1,f);
+ }
+ else
+ {
+ cout << " *AliTrack::Dump* Error : No decay track present." << endl;
+ }
+ }
+}
+//////////////////////////////////////////////////////////////////////////
+Double_t AliTrack::GetMomentum()
+{
+// Provide the value of the track 3-momentum
+ Ali3Vector p=Get3Vector();
+ return sqrt(p.Dot(p));
+}
+///////////////////////////////////////////////////////////////////////////
+Ali3Vector AliTrack::Get3Momentum()
+{
+// Provide the track 3-momentum
+ return (Ali3Vector)Get3Vector();
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t AliTrack::GetMass()
+{
+// Provide the particle mass
+ return fM;
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliTrack::GetCharge()
+{
+// Provide the particle charge
+ return fQ;
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t AliTrack::GetEnergy()
+{
+// Provide the particle's energy
+ return GetScalar();
+}
+///////////////////////////////////////////////////////////////////////////
+void AliTrack::Decay(Double_t m1,Double_t m2,Double_t thcms,Double_t phicms)
+{
+// Perform 2-body decay of current track
+// m1 : mass of decay product 1
+// m2 : mass of decay product 2
+// thcms : cms theta decay angle (in rad.) of m1
+// phicms : cms phi decay angle (in rad.) of m1
+
+ fNdec=2; // it's a 2-body decay
+
+// Compute the 4-momenta of the decay products in the cms
+// Note : p2=p1=pnorm for a 2-body decay
+ Double_t e1=((fM*fM)+(m1*m1)-(m2*m2))/(2.*fM);
+ Double_t e2=((fM*fM)+(m2*m2)-(m1*m1))/(2.*fM);
+ Double_t pnorm=(e1*e1)-(m1*m1);
+ if (pnorm>0.)
+ {
+ pnorm=sqrt(pnorm);
+ }
+ else
+ {
+ pnorm=0;
+ }
+
+ Double_t a[3];
+ a[0]=pnorm;
+ a[1]=thcms;
+ a[2]=phicms;
+ Ali3Vector p;
+ p.SetVector(a,"sph");
+
+ Ali4Vector pprim1;
+ pprim1.SetVector(e1,p);
+
+ Ali4Vector pprim2;
+ p*=-1;
+ pprim2.SetVector(e2,p);
+
+ // Determine boost parameters from the parent particle
+ Double_t E=GetScalar();
+ p=Get3Vector();
+ Ali4Vector pmu;
+ pmu.SetVector(E,p);
+
+ AliBoost q;
+ q.Set4Momentum(pmu);
+
+ Ali4Vector p1=q.Inverse(pprim1); // Boost decay product 1
+ Ali4Vector p2=q.Inverse(pprim2); // Boost decay product 2
+
+ // Enter the boosted data into the decay tracks array
+ if (fDecays)
+ {
+ fDecays->Delete();
+ delete fDecays;
+ }
+ fDecays=new TObjArray();
+
+ fDecays->Add(new AliTrack);
+ ((AliTrack*)fDecays->At(0))->Set4Momentum(p1);
+ fDecays->Add(new AliTrack);
+ ((AliTrack*)fDecays->At(1))->Set4Momentum(p2);
+
+// Set the mass values to m1 and m2 to omit roundoff errors
+ ((AliTrack*)fDecays->At(0))->SetMass(m1);
+ ((AliTrack*)fDecays->At(1))->SetMass(m2);
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliTrack::GetNdecay()
+{
+// Provide the number of decay produced tracks
+ return fNdec;
+}
+///////////////////////////////////////////////////////////////////////////
+AliTrack* AliTrack::GetDecayTrack(Int_t j)
+{
+// Provide decay produced track number j
+// Note : j=1 denotes the first decay track
+ if ((j >= 1) && (j <= fNdec))
+ {
+ return (AliTrack*)fDecays->At(j-1);
+ }
+ else
+ {
+ cout << " *AliTrack* decay track number : " << j << " out of range." << endl;
+ cout << " -- Decay track number 1 (if any) returned." << endl;
+ return (AliTrack*)fDecays->At(0);
+ }
+}
+///////////////////////////////////////////////////////////////////////////
--- /dev/null
+#ifndef ALITRACK_H
+#define ALITRACK_H
+///////////////////////////////////////////////////////////////////////////
+// Class AliTrack
+// Handling of the attributes of a reconstructed particle track.
+//
+// Coding example :
+// ----------------
+//
+// Float_t a[4]={195.,1.2,-0.04,8.5};
+// Ali4Vector pmu;
+// pmu.SetVector(a,"car");
+// AliTrack t1;
+// t1.Set4Momentum(pmu);
+//
+// Float_t b[3]={1.2,-0.04,8.5};
+// Ali3Vector p;
+// p.SetVector(b,"car");
+// AliTrack t2;
+// t2.Set3Momentum(p);
+// t2.SetCharge(0);
+// t2.SetMass(1.115);
+//
+// t1.Info();
+// t2.Info();
+//
+// Float_t pi=acos(-1.);
+// Float_t thcms=0.2*pi; // decay theta angle in cms
+// Float_t phicms=pi/4.; // decay theta angle in cms
+// Float_t m1=0.938;
+// Float_t m2=0.140;
+// t2.Decay(m1,m2,thcms,phicms); // Track t2 decay : Lambda -> proton + pion
+//
+// t2.List();
+//
+// Int_t ndec=t2.GetNdecay();
+// AliTrack* d1=t2.GetDecayTrack(1); // Access to decay track number 1
+// AliTrack* d2=t2.GetDecayTrack(2); // Access to decay track number 2
+//
+// Note : All quantities are in GeV, GeV/c or GeV/c**2
+//
+//--- NvE 10-jul-1997 UU-SAP Utrecht
+//--- Modified : NvE 06-apr-1999 UU-SAP Utrecht to inherit from Ali4Vector
+///////////////////////////////////////////////////////////////////////////
+
+#include "TObject.h"
+#include "TObjArray.h"
+
+#include "AliBoost.h"
+
+class AliTrack : public TObject,public Ali4Vector
+{
+ public:
+ AliTrack(); // Default constructor
+ ~AliTrack(); // Destructor
+ void Reset(); // Reset all values to 0
+ void Set4Momentum(Ali4Vector& p); // Set track 4-momentum
+ void Set3Momentum(Ali3Vector& p); // Set track 3-momentum
+ void SetMass(Double_t m); // Set particle mass
+ void SetCharge(Float_t q); // Set particle charge
+ void Info(TString f="car"); // Print track information for coord. frame f
+ void List(TString f="car"); // Print track and decay level 1 information for coord. frame f
+ void ListAll(TString f="car"); // Print track and all decay level information for coord. frame f
+ Ali3Vector Get3Momentum(); // Provide track 3-momentum
+ Double_t GetMomentum(); // Provide value of track 3-momentum
+ Double_t GetMass(); // Provide particle mass
+ Float_t GetCharge(); // Provide particle charge
+ Double_t GetEnergy(); // Provide particle total energy
+ void Decay(Double_t m1,Double_t m2,Double_t thcms,Double_t phicms); // Perform 2-body decay
+ Int_t GetNdecay(); // Provide number of decay products
+ AliTrack* GetDecayTrack(Int_t j); // Access to decay produced track number j
+
+ protected:
+ Double_t fM; // The mass of the particle
+ Float_t fQ; // The charge of the particle
+ Int_t fNdec; // The number of decay products
+ TObjArray* fDecays; // The array of decay produced tracks for output
+
+ private:
+ void Dump(AliTrack* t,Int_t n,TString f); // Recursively print all decay levels
+
+ ClassDef(AliTrack,1) // Class definition to enable ROOT I/O
+};
+#endif
--- /dev/null
+#include "AliVertex.h"
+
+ClassImp(AliVertex) // Class implementation to enable ROOT I/O
+
+AliVertex::AliVertex()
+{
+// Default constructor
+// All variables initialised to 0
+// Initial maximum number of tracks is set to the default value
+// Initial maximum number of sec. vertices is set to the default value
+ fNvmax=0;
+ fVertices=0;
+ Reset();
+ SetNtinit();
+ SetNvmax();
+}
+///////////////////////////////////////////////////////////////////////////
+AliVertex::AliVertex(Int_t n)
+{
+// Create a vertex to hold initially a maximum of n tracks
+// All variables initialised to 0
+ fNvmax=0;
+ fVertices=0;
+ Reset();
+ if (n > 0)
+ {
+ SetNtinit(n);
+ }
+ else
+ {
+ cout << endl;
+ cout << " *AliVertex* Initial max. number of tracks entered : " << n << endl;
+ cout << " This is invalid. Default initial maximum will be used." << endl;
+ cout << endl;
+ SetNtinit();
+ }
+ SetNvmax();
+}
+///////////////////////////////////////////////////////////////////////////
+AliVertex::~AliVertex()
+{
+// Default destructor
+ if (fVertices) delete fVertices;
+ fVertices=0;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliVertex::SetNvmax(Int_t n)
+{
+// Set the initial maximum number of (secondary) vertices
+ if (n > 0)
+ {
+ fNvmax=n;
+ }
+ else
+ {
+ fNvmax=1;
+ }
+ if (fVertices) delete fVertices;
+ fVertices=new TObjArray(fNvmax);
+}
+///////////////////////////////////////////////////////////////////////////
+void AliVertex::Reset()
+{
+// Reset all variables to 0
+// The max. number of tracks is set to the initial value again
+// The max. number of vertices is set to the default value again
+
+ AliJet::Reset();
+
+ fNvtx=0;
+ if (fNvmax>0) SetNvmax(fNvmax);
+}
+///////////////////////////////////////////////////////////////////////////
+void AliVertex::Add(AliJet& j)
+{
+// Add the tracks of a jet to the vertex
+ AliTrack* tj;
+ for (Int_t i=1; i<=j.GetNtracks(); i++)
+ {
+ tj=j.GetTrack(i);
+ AliJet::Add(tj);
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliVertex::Add(AliVertex& v)
+{
+// 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.
+//
+ if (fNvtx == fNvmax) // Check if maximum vertex number is reached
+ {
+ fNvmax++;
+ fVertices->Expand(fNvmax);
+ }
+
+ // Update 4-momentum for current vertex
+ fNvtx++;
+ fVertices->Add(&v);
+ (Ali4Vector)(*this)+=v;
+}
+///////////////////////////////////////////////////////////////////////////
+void AliVertex::Info(TString f)
+{
+// Provide vertex information within the coordinate frame f
+ cout << " *AliVertex::Info* Invmass : " << GetInvmass()
+ << " Charge : " << GetCharge() << " Momentum : " << GetMomentum()
+ << " Ntracks : " << GetNtracks() << " Nvertices : " << fNvtx << endl;
+ cout << " ";
+ Ali4Vector::Info(f);
+ cout << " Position";
+ AliPosition::Info(f);
+}
+///////////////////////////////////////////////////////////////////////////
+void AliVertex::List(TString f)
+{
+// Provide primary track and sec. vertex information within the coordinate frame f
+
+ Info(f); // Information of the current vertex
+
+ // The tracks of this vertex
+ AliTrack* t;
+ for (Int_t it=1; it<=GetNtracks(); it++)
+ {
+ t=GetTrack(it);
+ if (t)
+ {
+ cout << " ---Track no. " << it << endl;
+ cout << " ";
+ t->Info(f);
+ }
+ else
+ {
+ cout << " *AliVertex::List* Error : No track present." << endl;
+ }
+ }
+
+ // The secondary vertices of this vertex
+ AliVertex* v;
+ for (Int_t iv=1; iv<=GetNvertices(); iv++)
+ {
+ v=GetVertex(iv);
+ if (v)
+ {
+ cout << " ---Level 1 sec. vertex no. " << iv << endl;
+ cout << " ";
+ v->Info(f);
+ }
+ else
+ {
+ cout << " *AliVertex::List* Error : No sec. vertex present." << endl;
+ }
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void AliVertex::ListAll(TString f)
+{
+// Provide complete (sec) vertex and (decay) track info within the coordinate frame f
+
+ Info(f); // Information of the current vertex
+
+ // The tracks of this vertex
+ AliTrack* t;
+ for (Int_t it=1; it<=GetNtracks(); it++)
+ {
+ t=GetTrack(it);
+ if (t)
+ {
+ cout << " ---Track no. " << it << endl;
+ cout << " ";
+ t->ListAll(f);
+ }
+ else
+ {
+ cout << " *AliVertex::ListAll* Error : No track present." << endl;
+ }
+ }
+
+ AliVertex* v=this;
+ Dump(v,1,f); // Information of all sec. vertices
+}
+//////////////////////////////////////////////////////////////////////////
+void AliVertex::Dump(AliVertex* v,Int_t n,TString f)
+{
+// Recursively provide the info of all secondary vertices of this vertex
+ AliVertex* vs;
+ for (Int_t iv=1; iv<=v->GetNvertices(); iv++)
+ {
+ vs=v->GetVertex(iv);
+ if (vs)
+ {
+ cout << " ---Level " << n << " sec. vertex no. " << iv << endl;
+ cout << " ";
+ vs->Info(f);
+
+ // The tracks of this vertex
+ AliTrack* t;
+ for (Int_t it=1; it<=vs->GetNtracks(); it++)
+ {
+ t=vs->GetTrack(it);
+ if (t)
+ {
+ cout << " ---Track no. " << it << endl;
+ cout << " ";
+ t->ListAll(f);
+ }
+ else
+ {
+ cout << " *AliVertex::Dump* Error : No track present." << endl;
+ }
+ }
+
+ // Go for next sec. vertex level of this sec. vertex recursively
+ Dump(vs,n+1,f);
+ }
+ else
+ {
+ cout << " *AliVertex::Dump* Error : No sec. vertex present." << endl;
+ }
+ }
+}
+//////////////////////////////////////////////////////////////////////////
+Int_t AliVertex::GetNvertices()
+{
+// Return the current number of (secondary) vertices
+ return fNvtx;
+}
+///////////////////////////////////////////////////////////////////////////
+AliVertex* AliVertex::GetVertex(Int_t i)
+{
+// Return the i-th (secondary) vertex of the current vertex
+ return (AliVertex*)fVertices->At(i-1);
+}
+///////////////////////////////////////////////////////////////////////////
--- /dev/null
+#ifndef ALIVERTEX_H
+#define ALIVERTEX_H
+///////////////////////////////////////////////////////////////////////////
+// Class AliVertex
+// Creation and investigation of an AliVertex.
+// An AliVertex can be constructed by adding AliTracks and/or AliJets.
+//
+// Note : Also (secondary) vertices can be added to a vertex.
+//
+// Coding example to make 3 vertices v1, v2 and v3.
+// ------------------------------------------------
+// v1 contains the tracks 1,2,3 and 4
+// v2 contains the tracks 5,6 and 7
+// v3 contains the jets 1 and 2
+//
+// AliTrack t1,t2,t3,t4,t5,t6,t7;
+// ...
+// ... // code to fill the track data
+// ...
+//
+// AliJet j1,j2;
+// ...
+// ... // code to fill the jet data
+// ...
+//
+// AliVertex v1(5);
+//
+// v1.Add(t1);
+// v1.Add(t2);
+// v1.Add(t3);
+// v1.Add(t4);
+//
+// Float_t r1[3]={2.4,0.1,-8.5};
+// v1.SetPosition(r1,"car");
+//
+// AliVertex v2(2);
+// v2.Add(t5);
+// v2.Add(t6);
+// v2.Add(t7);
+//
+// Float_t r2[3]={1.6,-3.2,5.7};
+// v2.SetPosition(r2,"car");
+//
+// AliVertex v3;
+//
+// v3.Add(j1);
+// v3.Add(j2);
+//
+// Float_t r3[3]={6.2,4.8,1.3};
+// v3.SetPosition(r3,"car");
+//
+// v1.Info("sph");
+// v2.ListAll();
+// v3.List("cyl");
+//
+// Float_t e1=v1.GetEnergy();
+// Ali3Vector p1=v1.Get3Momentum();
+// Float_t loc[3];
+// v1.GetPosition(loc,"sph");
+// AliPosition r=v2.GetPosition();
+// r.Info();
+// Int_t nt=v2.GetNtracks();
+// AliTrack* tv=v2.GetTrack(1); // Access track number 1 of Vertex v2
+//
+// Specify the vertices v2 and v3 as secondary vertices of v1
+//
+// v1.Add(v2);
+// v1.Add(v3);
+//
+// v1.List();
+//
+// Int_t nv=v1.GetNvtx();
+// AliVertex* vx=v1.GetVertex(1); // Access 1st secondary vertex of v1
+// Float_t e=vx->GetEnergy();
+//
+// Float_t M=v1.GetInvmass();
+//
+// Reconstruct Vertex v1 from scratch
+//
+// v1.Reset();
+// v1.SetNvmax(25); // Increase initial no. of sec. vertices
+// v1.Add(t3);
+// v1.Add(t7);
+// v1.Add(j2);
+// Float_t pos[3]={7,9,4};
+// v1.SetPosition(pos,"car");
+//
+// Note : All quantities are in GeV, GeV/c or GeV/c**2
+//
+//--- NvE 04-apr-1998 UU-SAP Utrecht
+//--- Modified : NvE 08-apr-1999 UU-SAP Utrecht to inherit from AliJet
+///////////////////////////////////////////////////////////////////////////
+
+#include <iostream.h>
+#include <math.h>
+
+#include "TObject.h"
+#include "TObjArray.h"
+
+#include "AliJet.h"
+#include "AliPosition.h"
+
+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
+ 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
+
+ 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
+
+ private:
+ void Dump(AliVertex* v,Int_t n,TString f); // Recursively print all sec. vertices
+
+ ClassDef(AliVertex,1) // Class definition to enable ROOT I/O
+};
+#endif
--- /dev/null
+############################### TPC Makefile ##################################
+
+# Include machine specific definitions
+
+include $(ALICE_ROOT)/conf/GeneralDef
+include $(ALICE_ROOT)/conf/MachineDef.$(ALICE_TARGET)
+
+PACKAGE = RALICE
+
+# C++ sources
+
+SRCS = Ali3Vector.cxx Ali4Vector.cxx AliBoost.cxx AliCalcluster.cxx \
+ AliCalmodule.cxx AliCalorimeter.cxx AliInvmass.cxx AliJet.cxx \
+ AliMath.cxx AliPosition.cxx AliRandom.cxx AliSample.cxx AliSignal.cxx \
+ AliTrack.cxx AliVertex.cxx
+
+
+# C++ Headers
+
+HDRS = $(SRCS:.cxx=.h) RALICELinkDef.h
+
+# Library dictionary
+
+DICT = RALICECint.cxx
+DICTH = $(DICT:.cxx=.h)
+DICTO = $(patsubst %.cxx,$(ALICE_TARGET)/%.o,$(DICT))
+
+# FORTRAN Objectrs
+
+FOBJS = $(FSRCS:.f=.o)
+
+# C Objects
+
+COBJS = $(CSRCS:.c=.o)
+
+# C++ Objects
+
+OBJS = $(patsubst %.cxx,$(ALICE_TARGET)/%.o,$(SRCS)) $(DICTO)
+
+# C++ compilation flags
+
+CXXFLAGS = $(CXXOPTS) -I$(ROOTSYS)/include -I. -I$(ALICE_ROOT)/include/
+
+# FORTRAN compilation flags
+
+FFLAGS = $(FOPT)
+
+##### TARGETS #####
+
+# Target
+
+SLIBRARY = $(LIBDIR)/libRALICE.$(SL)
+
+default: $(SLIBRARY)
+
+$(LIBDIR)/libRALICE.$(SL): $(OBJS)
+
+$(DICT): $(HDRS)
+
+depend: $(SRCS)
+
+TOCLEAN = $(OBJS) *Cint.h *Cint.cxx
+
+############################### General Macros ################################
+
+include $(ALICE_ROOT)/conf/GeneralMacros
+
+############################ Dependencies #####################################
+
+-include Make-depend
--- /dev/null
+////////////////////////////////////////////////////////////////////////////////
+// All classes of RALICE
+// This class list is used to create the RALICE dictionary via rootcint
+// in the automatic installation of the ROOT loadable library.
+//
+// Note : Headers have also to be entered into the list in file allhead.h
+//
+//--- NvE 12-apr-1998 UU-SAP Utrecht
+////////////////////////////////////////////////////////////////////////////////
+
+#ifdef __CINT__
+ #pragma link off all globals;
+ #pragma link off all classes;
+ #pragma link off all functions;
+
+ #pragma link C++ class AliMath;
+ #pragma link C++ class AliSample;
+ #pragma link C++ class AliRandom;
+ #pragma link C++ class Ali3Vector;
+ #pragma link C++ class Ali4Vector;
+ #pragma link C++ class AliBoost;
+ #pragma link C++ class AliPosition;
+ #pragma link C++ class AliSignal;
+ #pragma link C++ class AliCalorimeter;
+ #pragma link C++ class AliCalmodule;
+ #pragma link C++ class AliCalcluster;
+ #pragma link C++ class AliTrack;
+ #pragma link C++ class AliJet;
+ #pragma link C++ class AliVertex;
+ #pragma link C++ class AliInvmass;
+#endif
+
--- /dev/null
+////////////////////////////////////////////////////////////////////////////////
+// All headers of RALICE
+// This header list is used to create the RALICE dictionary via rootcint
+// in the automatic installation of the ROOT loadable library.
+//
+// Note : Class names have also to be entered into the list in file linkdef.h
+//
+//--- NvE 12-apr-1998 UU-SAP Utrecht
+////////////////////////////////////////////////////////////////////////////////
+
+#include "AliMath.h"
+#include "AliRandom.h"
+#include "AliSample.h"
+#include "Ali3Vector.h"
+#include "Ali4Vector.h"
+#include "AliBoost.h"
+#include "AliPosition.h"
+#include "AliSignal.h"
+#include "AliCalorimeter.h"
+#include "AliCalmodule.h"
+#include "AliCalcluster.h"
+#include "AliTrack.h"
+#include "AliJet.h"
+#include "AliVertex.h"
+#include "AliInvmass.h"
--- /dev/null
+//*CMZ : 1.00/01 04/07/97 16.43.13 by Nick van Eijndhoven (UU/CERN)
+//*-- Author : Nick van Eijndhoven (UU/CERN) 04/07/97
+
+This RALICE cmz file is an attempt to provide an Object Oriented framework,
+consisting of C++ classes, in which event reconstruction of the ALICE
+detector data can be performed.
+In switching to Object Oriented programming, I myself have started to
+perform the WA93 and WA98 data analysis within the ROOT [1] framework.
+Having seen the great advantages of this system I have started to make my
+C++ classes more general in order to use them as an onset for an ALICE
+reconstruction (and physics analysis) framework.
+The RALICE package can be compiled on all platforms using the GNU G++
+compiler and the various classes can be used in standalone mode.
+However, running the programs within the ROOT framework greatly enlarges
+the analysis capabilities (e.g. histogramming, fitting, graphics etc...).
+In addition the high level of interactivity of the ROOT/CINT system allows
+program development without the time consuming compile/link/load/execute cycle,
+whereas also the ROOT tree output format provides a completely machine
+independent data format providing efficient and easy to use data access
+capable to cope with the most complex data analyses programs.
+
+Only the (proposed) C++ ANSI standard is used in the source code and as
+such is fully compatible with all standard C++ compilers as well as with
+the ROOT/CINT interpreting system.
+
+The comments in the source code are placed in the locations proposed
+in the ROOT manual pages [1] such that the automatic source code
+documentation system of ROOT can be used directly from the source code.
+This has turned out to be very convenient (time saving) and guarantees
+always updated documentation compatible with the current source code.
+
+Coding conventions :
+--------------------
+In order not to clash with the (class) names of the ROOT framework
+and (future) packages of other groups, a few rules concerning names
+of classes, (member)functions and variables have to be obeyed.
+The rules are the following :
+
+ 1) Only (proposed) ANSI standard C++ is allowed, with an even stricter
+ requirement that the code should compile without any warnings
+ under the GNU g++, msvc++ and the native C++ compilers of HP
+ and DECAlpha machines.
+ This will assure the programs to run on all standard ALICE platforms.
+ 2) Class names start with "Ali" followed by an uppercase character,
+ all other characters are lowercase.
+ Example : AliCalorimeter
+ In this way the RALICE class names will NEVER clash with the ones
+ of ROOT whereas the probability of a clash with the class names of
+ other group's code (e.g. ATLAS, CDF, PHENIX etc...) is minimised.
+ To prevent name clashes within the various (future) ALICE packages,
+ please refer to the general note at the end.
+ 3) Names of detector specific classes should start with "Ali" followed
+ by the detector name in uppercase, all other characters are lowercase
+ except the first character following the detector name, which has to
+ be uppercase..
+ Example : AliTPCSegment or AliPPCTiming.
+ These detector specific classes should only be introduced when there
+ is really a need for it.
+ E.g. when a track segment of the TPC and ITS have a lot in common
+ it would be better to introduce a general AliTracksegment class
+ instead of AliTPCSegment and AliITSSegment classes.
+ 4) Class headers should be under the control of "#ifndef" and the name
+ should consist of "CLASSNAME_H" (i.e. the classname in uppercase).
+ Example : #ifndef ALITRACK_H
+ #define ALITRACK_H
+ In this way also the ifdefs will be unique and prevents the danger
+ of having the name of an ifdef being the same as a Classname.
+ 5) The private area in the class header has to be defined as the last item.
+ Macros, like the ROOT ClassDef() statement (if needed) must be put
+ appear at the right location, i.e. just before the "};" of the
+ class definition.
+ 6) Names of member functions should start with a capital character
+ and should NOT contain underscores (which drop out with HTML).
+ From the name it should be clear what the functionality is and
+ capital characters should be used to indicate various "words".
+ Example : AliTrack::Set3Momentum(float* p)
+ 7) Names of datamembers of a class should start with a lowercase "f"
+ and the next character has to be uppercase.
+ Example : float fEnergy
+ This will allow directly identification of datamembers in the code.
+ The names of all other local variables may be chosen freely by the
+ author.
+ Note : It is recommended to use only lowercase characters
+ for local variables.
+ 8) Names of global variables should start with "gAli" and the next
+ character has to be uppercase.
+ Example : gAliRun
+ This will allow directly identification of global variables in the
+ code and will not clash with the existing ROOT globals like
+ for instance gDirectory etc...
+ Note : Usage of global variables should be avoided as much as
+ possible. Most of the data transfer should go via the classes
+ and their member functions (data hiding).
+ 9) Comments should be placed at the positions as outlined in the ROOT docs.
+ This will enable the automatic HTML machinery of ROOT.
+10) Each class header should contain a short description of the class
+ functionality including some examples.
+
+General note :
+--------------
+Within the ALICE software pool it may happen that e.g. in simulation
+applications one wants to define for instance a Track class which
+contains as data members some additional information (e.g. which was
+the corresponding parent particle) compared to the AliTrack class.
+Since objects reconstructed from real data will always contain the
+minimal amount of information compared to e.g. objects from simulation,
+it is in the above case then necessary to introduce a new class
+AliSTrack (simulation track).
+Obviously such a newly defined object (AliSTrack) can be derived from
+the reconstruction object (AliTrack) and just have some data members
+and/or memberfunctions added to it.
+In such a way maximum flexibility is provided within every (future)
+ALICE project, whereas all produced data can always be analysed using
+the RALICE tools.
+In view of this it might even be preferred to impose as a convention
+for future projects to adopt a unique prefix for their specific classes.
+For example the prefixes "AliS" and "AliD" could be used to indicate
+the simulation and DAQ specific classes respectively.
+
+Installation :
+--------------
+The RALICE library can be automatically installed using the automatic CMZ
+installation procedure.
+
+Available installation modes :
+------------------------------
+cmz -install ralice --> GNU G++ loadable libralice.a
+cmz -install ralice shared --> GNU G++ loadable shared library ralice.sl
+cmz -install ralice root --> ROOT (G++ based) loadable library ralice.sl
+cmz -install ralice - msvc --> MSVC++ loadable library
+cmz -install ralice shared msvc --> MSVC++ loadable shared library ralice.dll
+cmz -install ralice root msvc --> ROOT (MSVC++ based) loadable library ralice.dll
+cmz -install ralice - hpcc --> HP CC loadable library
+cmz -install ralice shared hpcc --> HP CC loadable shared library ralice.sl
+cmz -install ralice root hpcc --> ROOT (HP CC based) loadable library ralice.sl
+
+
+[1] http://root.cern.ch
+
+
+
+ Nick van Eijndhoven
+ Subatomic Physics Dept.
+ Utrecht University/NIKHEF
+ P.O. Box 80.000
+ NL-3508 TA Utrecht
+ The Netherlands
+ Email: nick@phys.uu.nl
+ WWW: http://www.phys.uu.nl/~nick
--- /dev/null
+//*CMZ : 1.00/02 04/07/97 16.48.20 by Nick van Eijndhoven (UU/CERN)
+//*-- Author : Nick van Eijndhoven (UU/CERN) 04/07/97
+//////////////////////////////////////////////////////////////////////////
+// History of updates //
+//////////////////////////////////////////////////////////////////////////
+04-jul-1997 NvE First release of the package
+11-jul-1997 NvE Class Invmass added to /GENUTILS and directories /PHYSICS
+ and /TRACKING created to hold Jet and Track classes resp.
+15-aug-1997 NvE Pilot patch *RALICE removed and automatic installation
+ now completely performed by pure CMZ commands in the
+ INSTALL macro.
+12-sep-1997 NvE Signal of central cluster module and position in user
+ coordinates for all modules introduced in the Calorimetry
+ package.
+ Also all "float vec[]" replaced by "float* vec" in the
+ Lorbo class for stylistic reasons.
+20-sep-1997 NvE Name of member function "getc()" of Calorimeter changed to
+ "get_c()" to avoid clash with native C function.
+ Also protection added in Calorimeter::group() against cases
+ where no modules fired.
+23-sep-1997 NvE Macro "install" modified to also enable automatic installation
+ of a G++ or ROOT loadable shared library 'ralice.sl'.
+24-sep-1997 NvE Class "Calorimeter" optimized by introduction of pointers
+ for the Module matrix instead of an indexed array.
+26-sep-1997 NvE Directory /ROOT added to hold the universal ROOT dictionaries
+ (dict.h and dict.cc) as generated by "rootcint" for all the
+ classes of this RALICE package.
+ This will allow automatic installation of the RALICE library
+ without running "rootcint".
+30-sep-1997 NvE 2-body decay added to class Track and invokation of "delete []"
+ only performed in case the object exists in the destructors
+ of the Calorimeter, Jet and Track classes.
+02-oct-1997 NvE Installation for MSVC++ added into macro "install" and several
+ "return 0.;" statements added in some "Sample" member functions.
+ The latter caused MSVC++ compiler errors.
+ Also 'int ..=fabs(int..)' changed to 'int ..=int(fabs(int..))'
+ in "Cluster::Add()" to prevent G++ warning.
+03-oct-1997 NvE The above code is now ANSI compatible, but MSVC++ 5.0 still
+ gives an error when comp. extensions are switched off (flag /Za).
+ The reason is probably an incorrect implementation of the ANSI
+ type conversion rules in MSVC++. To satisfy MSVC++ for the moment
+ the code in in "Cluster::Add()" has been changed to
+ 'int ..=int(fabs(double(int..)))' so that at least a Microsoft
+ compatible RALICE.DLL can be produced. As soon as G++ is able
+ to produce Microsoft compatible dll's the support for MSVC++
+ will be dropped for ralice.cmz.
+ Message to all : USE G++ AS COMPILATION STANDARD !
+ Also CMZ standard adopted to include a "$" in the DECK names
+ which contain KEEP definitions.
+13-oct-1997 NvE Compiler directive for g++ changed in "flogon" and "install"
+ to provide warning for each non-ANSI C++ statement.
+ Also ROOT libs linked with g++ in "install" in case of
+ ROOT loadable shared library creation.
+ Unused variables removed in "Calorimeter::add_ring" and
+ "Track::decay".
+14-oct-1997 NvE Strict ANSI requirements dropped for ROOT loadable lib creation
+ to avoid many harmless warnings from rootcint generated code.
+ Class "Random" added to /GENUTILS and /ROOT updated accordingly.
+14-oct-1997 NvE Class "Random" extended with user definable distribution and
+ /ROOT updated accordingly.
+19-oct-1997 NvE Class "Random" extended with vector facility for user
+ defined dist. and also bug fixed. /ROOT updated accordingly.
+27-oct-1997 NvE Class "Random" optimised for user definable distribution and
+ /ROOT updated accordingly.
+10-dec-1997 NvE Creation of MSVC++ loadable library introduced in "install".
+12-dec-1997 NvE Keep and Deck in /ROOT given unique names to enable loading
+ of this package together with other packages.
+20-jan-1998 NvE New ralicedict produced because of new ROOT 1.03/09 version.
+23-mar-1998 NvE Macro "install" modified to support HP-CC.
+ Multiple usage of for-loop index names removed to satisfy HP-CC
+ compiler (although the original code was ANSI C++).
+ Cluster data just reset to 'empty cluster' in case of cluster
+ 'creation' or 'start' for edge module or no-signal module.
+31-mar-1998 NvE Array index bug fixed in class Jet.
+03-apr-1998 NvE Jet and Invmass classes optimised by using pointers and
+ compiler options /GD and /DWIN32 introduced for msvc++
+ in order to optimize for windows DLL's and ensure correct
+ pre-processing in RConfig.h for ROOT library creation.
+13-apr-1998 NvE ROOTCINT processing implemented in $KUMACS/INSTALL and
+ contents of /ROOT updated accordingly.
+ Class Vertex added to /TRACKING and classes Track and Jet
+ protected against modifying of returned addresses.
+14-apr-1998 NvE "rm .def" removed for Unix systems in $KUMACS/INSTALL since
+ on Unix no .def file is created by ROOTCINT.
+17-apr-1998 NvE Return type Track& changed to Track* in the get_trk()
+ member functions of the Jet and Vertex classes.
+21-apr-1998 NvE Classes Track, Vertex, Jet, Lorbo and Invmass derived from
+ TObject to enable ROOT I/O.
+ $KUMACS/FLOGON and $KUMACS/INSTALL updated accordingly.
+26-apr-1998 NvE Classes Module, Cluster and Calorimeter derived from TObject
+ to enable ROOT I/O.
+ Functionality of class Calorimeter has also been extended to
+ enable re-clustering of original Module signals.
+28-apr-1998 NvE Decay produced Tracks in class Track stored also into TObjArray
+ to allow direct I/O.
+03-jun-1998 NvE MSVC++ compiler option /DWIN32 removed in $KUMACS/INSTALL
+ because of bug fix in ROOT version 2.00.
+ Printout of version number introduced in $KUMACS/FLOGON and
+ DOC section updated for new email and www address.
+17-aug-1998 NvE New ALICE coding conventions adopted and automatic marking
+ of 'edge modules' introduced in class AliCalorimeter.
+ The automatic treatment of edge modules required new member-
+ functions EdgeUp() and EdgeDown() in class AliModule.
+19-aug-1998 NvE Call to LoadMatrix() added in AliCalorimeter::Group()
+ to correctly treat different event sizes (i.e. empty events)
+ when reading data from input file(s).
+25-aug-1998 NvE AliCalorimeter::GetModule() introduced to provide quick
+ access to all 'fired' modules.
+27-aug-1998 NvE GetSignal() and GetClusteredSignal() of AliModule now
+ return 0 in case the module was marked as dead.
+ Because of the above, protection against row=col=0
+ introduced in AliCalorimeter::Group().
+ To enable access of modified calorimeter module data via
+ AliCalorimeter::GetModule() after reading in, the "fModules"
+ elements are updated in AliCalorimeter::LoadMatrix().
+16-sep-1998 NvE Correct deletion of the module matrix introduced in the
+ AliCalorimeter destructor and AliCalorimeter::LoadMatrix().
+16-oct-1998 NvE New ALICE coding conventions adopted. This resulted in
+ modifications in /$KUMACS and /ROOT (.cxx extension instead of .cc)
+ and /HEADERS (public class area in front of private).
+ #include <iostream.h> added in AliCluster (thanks to A.Zvyagin).
+19-oct-1998 NvE Protections agains crazy row,col numbers introduced in
+ AliCalorimeter (thanks to E. van der Pijll).
+22-oct-1998 NvE Bug fixed in setting edge values in SetDead() and SetAlive()
+ of AliCalorimeter (thanks to E. van der Pijll).
+ Also row and column dispersion introduced in AliCluster.
+24-oct-1998 NvE Memberfunctions DrawModules() and DrawClusters() introduced
+ in AliCalorimeter.
+ Also AliRandom derived from TObject to enable ROOT I/O.
+26-oct-1998 NvE DrawModules() and DrawClusters() histogram axes defined such
+ that the X and Y axis corresponds to columns and rows resp.
+ in AliCalorimeter. This view matches the 'common sense'.
+06-nov-1998 NvE SetDirectory(0) invoked for the histos created in
+ DrawModules() and DrawClusters() to suppress global character
+ of the histo pointer.
+16-nov-1998 NvE New class AliMath implemented in /GENUTILS and /ROOT
+ updated accordingly.
+23-jan-1999 NvE New class AliSignal implemented and AliCalorimeter extended
+ to contain veto signal information.
+25-jan-1999 NvE Most of the private areas modified into protected areas to enable
+ the use of inheritance.
+ Also support for veto signal information introduced in AliCluster.
+04-feb-1999 NvE Cluster array for internal use removed in AliCalorimeter; now only
+ the "TObjArray* fClusters" is used for performance optimalisation.
+06-feb-1999 NvE New class AliPosition implemented for universal position handling
+ and AliSignal, AliCluster and AliVertex all derived from AliPosition.
+ AliModule inherits from AliSignal and arguments of SetPosition()
+ and GetPosition() in AliCalorimeter updated to match the AliSignal
+ inheritance of AliModule.
+11-mar-1999 NvE All types changed to Int_t, Float_t etc... in order to guarantee
+ portability of all data within the ROOT framework.
+12-mar-1999 NvE Default destructors added in AliInvmass, AliLorbo and AliSample.
+30-mar-1999 NvE New base class Ali3Vector introduced for universal 3D vector handling
+ and AliPosition derived from it.
+ In AliSignal::Reset() memberfunction call introduced instead of setting
+ datamembers; this enhances flexibility for base class development.
+01-apr-1999 NvE "car" introduced as default frame for Ali3Vector::Print().
+03-apr-1999 NvE New base class Ali4Vector introduced for universal handling of
+ Lorentz 4-vectors.
+ Also class AliBoost introduced to perform Lorentz boosts on Ali4Vector
+ objects. This class will replace the old AliLorbo.
+ Memberfunction Print() renamed to Info() for Ali3Vector, Ali4Vector
+ and AliBoost to prevent clash with TObject::Print().
+05-apr-1999 NvE Functionality of Ali3Vector and Ali4Vector extended with += etc...
+06-apr-1999 NvE Class AliTrack derived from Ali4Vector and updated accordingly.
+ Also class AliJet modified by deriving it from AliTrack.
+08-apr-1999 NvE Class AliVertex derived from AliJet and updated accordingly.
+ Also recursive printing facilities List() and ListAll() introduced
+ for AliTrack, AliJet and AliVertex.
+09-apr-1999 NvE Obsolete classes AliLorbo and AliInvmass removed.
+ A new AliInvmass class is planned to provide invariant mass reconstruction
+ functionality to serve e.g. di-lepton and pi0 studies.
+16-apr-1999 NvE New class AliInvmass introduced for invariant mass and comb. bkg. reconstruction.
+21-apr-1999 NvE AliBoost::SetGamma() introduced for enhanced accuracy in boost parameters.
+27-apr-1999 NvE AliModule and AliCluster classes renamed to AliCalmodule and AliCalcluster resp.
+ in view of future general Aliroot base classes AliModule and AliCluster.
+28-jun-1999 NvE Ali4Vector casting changed in AliVertex::Add() and some explicit type conversions
+ added and an unused variable removed in AliInvmass to prevent g++ errors/warnings
+ (thanks to Eugene van der Pijll).
--- /dev/null
+#!/bin/sh
+### Shell script to create a ROOT loadable HP-CC shared lib out of .cxx source code
+###
+### NvE 28-jun-1999 UU-SAP Utrecht
+#
+### The option strings for HP-CC shared lib compilation and linking ***
+### For the HP-CC ROOT loadable shared lib the strict requirements are ***
+### dropped to avoid many warnings from the rootcint generated code ***
+
+hpcomp="-c -s -z +z +a1 +w +DAportable -I$ROOTSYS/include"
+hproot="-c -s -z +z +a1 +DAportable -I$ROOTSYS/include"
+hplink="-L$ROOTSYS/lib/ -l*.sl -lm"
+
+rootcint zzzralicedict.cxx -c allhead.h linkdef.h
+
+CC $hproot *.cxx
+
+CC -b -o ralice.sl *.o
+
+rm zzzralicedict.*
+rm *.o
+
+echo '*** hpcclib done. Result in ralice.sl'
--- /dev/null
+Ali3Vector.cxx
+Ali4Vector.cxx
+AliBoost.cxx
+AliCalcluster.cxx
+AliCalmodule.cxx
+AliCalorimeter.cxx
+AliInvmass.cxx
+AliJet.cxx
+AliMath.cxx
+AliPosition.cxx
+AliRandom.cxx
+AliSample.cxx
+AliSignal.cxx
+AliTrack.cxx
+AliVertex.cxx
--- /dev/null
+////////////////////////////////////////////////////////////////////////////////
+// All classes of RALICE
+// This class list is used to create the RALICE dictionary via rootcint
+// in the automatic installation of the ROOT loadable library.
+//
+// Note : Headers have also to be entered into the list in file allhead.h
+//
+//--- NvE 12-apr-1998 UU-SAP Utrecht
+////////////////////////////////////////////////////////////////////////////////
+
+#ifdef __CINT__
+ #pragma link off all globals;
+ #pragma link off all classes;
+ #pragma link off all functions;
+
+ #pragma link C++ class AliMath;
+ #pragma link C++ class AliSample;
+ #pragma link C++ class AliRandom;
+ #pragma link C++ class Ali3Vector;
+ #pragma link C++ class Ali4Vector;
+ #pragma link C++ class AliBoost;
+ #pragma link C++ class AliPosition;
+ #pragma link C++ class AliSignal;
+ #pragma link C++ class AliCalorimeter;
+ #pragma link C++ class AliCalmodule;
+ #pragma link C++ class AliCalcluster;
+ #pragma link C++ class AliTrack;
+ #pragma link C++ class AliJet;
+ #pragma link C++ class AliVertex;
+ #pragma link C++ class AliInvmass;
+#endif
+
--- /dev/null
+@echo off
+rem ****************************************************************************
+rem * Script to create a relocatable MSVC++ DLL from *.cxx files
+rem *
+rem * Usage : mkdll
+rem *
+rem * This will create ralice.dll from all .h and .cxx files in the current dir
+rem *
+rem * In view of the ROOTCINT processing, the following two standard files
+rem * are always required :
+rem *
+rem * allhead.h ==> containing an include of all .h files
+rem * linkdef.h ==> containing the #pragma's to define all classes
+rem *
+rem ****************************************************************************
+rem *
+echo .
+echo === Automatic ROOT DLL production of file ralice.dll ===
+echo .
+rem *
+rem --- The option strings for MSVC++ DLL compilation and linking ***
+set mscomp=/nologo /c /TP /Za /MD /I%ROOTSYS%\include
+set msdll=/nologo /TP /Za /MD /LD /GD /I%ROOTSYS%\include
+set mslink=/ENTRY:_DllMainCRTStartup@12 %ROOTSYS%\lib\*.lib %MYLIBS%\*.lib
+rem *
+rootcint zzzralicedict.cxx -c allhead.h linkdef.h
+rem *
+cl %msdll% *.cxx /link %mslink% /OUT:ralice.dll
+rem *
+rem --- Delete all intermediate files
+del .def
+del containing
+del zzzralicedict.h
+del zzzralicedict.cxx
+del *.obj
+rem *
+echo *** mkdll done.
--- /dev/null
+@echo off
+rem ****************************************************************************
+rem * Script to create an MSVC++ LIB from *.cxx files
+rem *
+rem * Usage : mklib
+rem *
+rem * This will create ralice.lib from all .h and .cxx files in the current dir
+rem *
+rem * In view of the ROOTCINT processing, the following two standard files
+rem * are always required :
+rem *
+rem * allhead.h ==> containing an include of all .h files
+rem * linkdef.h ==> containing the #pragma's to define all classes
+rem *
+rem ****************************************************************************
+rem *
+echo .
+echo === Automatic ROOT LIB production of file ralice.lib ===
+echo .
+rem *
+rem --- The option strings for MSVC++ DLL compilation and linking ***
+set mscomp=/nologo /c /TP /Za /MD /I%ROOTSYS%\include
+set msdll=/nologo /TP /Za /MD /LD /GD /I%ROOTSYS%\include
+set mslink=/ENTRY:_DllMainCRTStartup@12 %ROOTSYS%\lib\*.lib %MYLIBS%\*.lib
+rem *
+rootcint zzzralicedict.cxx -c allhead.h linkdef.h
+rem *
+cl %mscomp% *.cxx
+lib /nologo *.obj /OUT:ralice.lib
+rem *
+rem --- Delete all intermediate files
+del .def
+del containing
+del zzzralicedict.h
+del zzzralicedict.cxx
+del *.obj
+rem *
+echo *** mklib done.
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff