[u/mrichter/AliRoot.git] / HLT / MUON / OnlineAnalysis / AliHLTMUONCalculations.cxx

/**************************************************************************
 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

// $Id$

////////////////////////////////////////////////////////////////////////////////
//
// Author: Artur Szostak
// Email:  artur@alice.phy.uct.ac.za | artursz@iafrica.com
//
////////////////////////////////////////////////////////////////////////////////

#include "AliHLTMUONCalculations.h"
#include "AliHLTMUONUtils.h"
#include "AliHLTMUONTriggerRecordsBlockStruct.h"
#include <cmath>

AliHLTFloat32_t AliHLTMUONCalculations::fgZf = -975.0;  // cm

AliHLTFloat32_t AliHLTMUONCalculations::fgQBLScaled
	= 3.0 * 2.99792458e8 / 1e9; // T.m.*c/1e9
	
AliHLTMUONParticleSign AliHLTMUONCalculations::fgSign = kSignUnknown;
AliHLTFloat32_t AliHLTMUONCalculations::fgPx = 0;  // GeV/c
AliHLTFloat32_t AliHLTMUONCalculations::fgPy = 0;  // GeV/c
AliHLTFloat32_t AliHLTMUONCalculations::fgPz = 0;  // GeV/c

AliHLTFloat32_t AliHLTMUONCalculations::fgSigmaX2 = 1.;  // cm^2
AliHLTFloat32_t AliHLTMUONCalculations::fgSigmaY2 = 1.;  // cm^2

AliHLTFloat32_t AliHLTMUONCalculations::fgMzx = 0;
AliHLTFloat32_t AliHLTMUONCalculations::fgMzy = 0;
AliHLTFloat32_t AliHLTMUONCalculations::fgCzx = 0;
AliHLTFloat32_t AliHLTMUONCalculations::fgCzy = 0;

AliHLTFloat32_t AliHLTMUONCalculations::fgIdealX1 = 0;  // cm
AliHLTFloat32_t AliHLTMUONCalculations::fgIdealY1 = 0;  // cm
AliHLTFloat32_t AliHLTMUONCalculations::fgIdealZ1 = -1603.5f;  // cm
AliHLTFloat32_t AliHLTMUONCalculations::fgIdealX2 = 0;  // cm
AliHLTFloat32_t AliHLTMUONCalculations::fgIdealY2 = 0;  // cm
AliHLTFloat32_t AliHLTMUONCalculations::fgIdealZ2 = -1703.5f;  // cm


bool AliHLTMUONCalculations::ComputeMomentum(
		AliHLTFloat32_t x1,
		AliHLTFloat32_t y1, AliHLTFloat32_t y2,
		AliHLTFloat32_t z1, AliHLTFloat32_t z2
	)
{
	/// Computes the momentum components based on the equations given in the
	///   ALICE dimuon spectrometer Technical Design Report (TDR-5): trigger section.
	///
	///   Reference: 
	///     "CERN/LHCC 2000-046
	///      Addendum 1 to ALICE TDR 5
	///      15 Dec 2000"
	///     Section 3.1.2 pages 144 and 145.
	///
	/// Input can be in meters, cm or mm. Output is in GeV/c.
	///
	/// \param x1  X coordinate of hit point 1 on the track.
	/// \param y1  Y coordinate of hit point 1 on the track.
	/// \param z1  Z coordinate of hit point 1 on the track.
	/// \param y2  Y coordinate of hit point 2 on the track.
	/// \param z2  Z coordinate of hit point 2 on the track.
	/// \return true if the momentum could be calculated and false otherwise.
	///    If true is returned then the estimated momentum can be fetched by the
	///    method calls: Px(), Py() and Pz() for the individual components.
	
	AliHLTFloat64_t z2mz1 = z2 - z1;
	if (z2mz1 == 0 or z1 == 0)
	{
		fgSign = kSignUnknown;
		fgPx = fgPy = fgPz = 0;
		return false;
	}
	AliHLTFloat64_t thetaTimesZf = (y1*z2 - y2*z1) / z2mz1;
	AliHLTFloat64_t xf = x1 * fgZf / z1;
	AliHLTFloat64_t yf = y2 - ((y2-y1) * (z2-fgZf)) / z2mz1;

	if (thetaTimesZf == 0)
	{
		fgSign = kSignUnknown;
		fgPx = fgPy = fgPz = 0;
		return false;
	}
	AliHLTFloat64_t pDivZf = (fgQBLScaled / thetaTimesZf);
	AliHLTFloat64_t p = pDivZf * fgZf;
	pDivZf = fabs(pDivZf);
	
	if (p < 0)
		fgSign = kSignMinus;
	else if (p > 0)
		fgSign = kSignPlus;
	else
		fgSign = kSignUnknown;
	
	fgPx = AliHLTFloat32_t( pDivZf * xf );
	fgPy = AliHLTFloat32_t( pDivZf * yf );
	AliHLTFloat64_t k = p*p - fgPx*fgPx - fgPy*fgPy;
	if (k > 0)
		fgPz = AliHLTFloat32_t( sqrt(k) );
	else
		fgPz = 0;
	// fgPz must be the same sign as fgZf else it could not have been measured.
	if (fgZf < 0) fgPz = -fgPz;

	return true;
}


AliHLTFloat32_t AliHLTMUONCalculations::QBL()
{
	// We have to convert back into Tesla metres.
	return fgQBLScaled * 1e9 / 2.99792458e8;
}


void AliHLTMUONCalculations::QBL(AliHLTFloat32_t value)
{
	// Note: 2.99792458e8/1e9 is the conversion factor for GeV.
	// It is c/1e9, where c is the speed of light.
	fgQBLScaled = value * 2.99792458e8 / 1e9;
}


AliHLTFloat32_t AliHLTMUONCalculations::ComputeMass(
		AliHLTFloat32_t massA,
		AliHLTFloat32_t pxA,
		AliHLTFloat32_t pyA,
		AliHLTFloat32_t pzA,
		AliHLTFloat32_t massB,
		AliHLTFloat32_t pxB,
		AliHLTFloat32_t pyB,
		AliHLTFloat32_t pzB
	)
{
	/// Calculates the invariant mass for a pair of particles.
	/// \param massA Mmass in GeV/c of particle A.
	/// \param pxA  X component of the momentum in GeV/c for particle A.
	/// \param pyA  Y component of the momentum in GeV/c for particle A.
	/// \param pzA  Z component of the momentum in GeV/c for particle A.
	/// \param massB  Mass in GeV/c of particle B.
	/// \param pxB  X component of the momentum in GeV/c for particle B.
	/// \param pyB  Y component of the momentum in GeV/c for particle B.
	/// \param pzB  Z component of the momentum in GeV/c for particle B.
	/// \return  The invariant mass in GeV/c^2 or -1 if there was a problem
	///          in the calculation due to bad input parameters.
	
	AliHLTFloat32_t massA2 = massA*massA;
	AliHLTFloat32_t massB2 = massB*massB;
	AliHLTFloat32_t energyA = sqrt(massA2 + pxA*pxA + pyA*pyA + pzA*pzA);
	AliHLTFloat32_t energyB = sqrt(massB2 + pxB*pxB + pyB*pyB + pzB*pzB);
	AliHLTFloat32_t mass2 = massA2 + massB2 + 2. * (energyA*energyB - pxA*pxB - pyA*pyB - pzA*pzB);
	if (mass2 < 0.) return -1.;
	return sqrt(mass2);
}


bool AliHLTMUONCalculations::FitLineToTriggerRecord(
		const AliHLTMUONTriggerRecordStruct& trigger
	)
{
	/// Straight lines are fitted in the ZX and ZY planes using a least
	/// squares fit to the coordinates in the trigger record.
	/// http://mathworld.wolfram.com/LeastSquaresFitting.html
	/// If this method returns true, then the fitted parameters can fetched
	/// using the method calls Mzx(), Mzy(), Czx() and Czy(). The lines are
	/// then given by: x = Mzx() * z + Czx() and y = Mzy() * z + Czy()
	/// The ideal coordinates are also calculated and can be fetched with
	/// the method calls: IdealX1(), IdealY1() and IdealZ1() for point on MT1,
	/// and IdealX2(), IdealY2() and IdealZ2() for point on MT2.
	/// \param trigger  The trigger record structure to which we fit a line.
	/// \return  true if the line could be fitted or false otherwise.
	///     The reason for failure could be either too few hits or the slopes
	///     Mzx() or Mzy() would be infinite, implying a line that is
	///     perpendicular to the z axis.
	
	AliHLTMUONParticleSign sign;
	bool hitset[4];
	AliHLTMUONUtils::UnpackTriggerRecordFlags(trigger.fFlags, sign, hitset);
	DebugTrace("hitset = {" << hitset[0] << ", " << hitset[1] << ", "
		<< hitset[2] << ", " << hitset[3] << "}"
	);
	
	return FitLineToTriggerRecord(trigger, hitset);
}


bool AliHLTMUONCalculations::FitLineToTriggerRecord(
		const AliHLTMUONTriggerRecordStruct& trigger,
		const bool hitset[4]
	)
{
	/// Performs a straight line fit like FitLineToTriggerRecord(trigger)
	/// but requires pree-decoded flags indicating which hits were set.
	/// \param trigger  The trigger record structure to which we fit a line.
	/// \param hitset  Flags indicating which hits were set in the trigger record.
	/// \return  true if the line could be fitted or false otherwise.
	
	AliHLTFloat32_t sumX = 0;
	AliHLTFloat32_t sumY = 0;
	AliHLTFloat32_t sumZ = 0;
	int n = 0;
	for (int i = 0; i < 4; i++)
	{
		if (hitset[i])
		{
			sumX += trigger.fHit[i].fX;
			sumY += trigger.fHit[i].fY;
			sumZ += trigger.fHit[i].fZ;
			n++;
		}
	}
	if (n < 2) return false;
	AliHLTFloat32_t meanX = sumX / AliHLTFloat32_t(n);
	AliHLTFloat32_t meanY = sumY / AliHLTFloat32_t(n);
	AliHLTFloat32_t meanZ = sumZ / AliHLTFloat32_t(n);
	
	AliHLTFloat32_t vSSzz = 0;
	AliHLTFloat32_t vSSzx = 0;
	AliHLTFloat32_t vSSzy = 0;
	for (int i = 0; i < 4; i++)
	{
		if (hitset[i])
		{
			vSSzz += (trigger.fHit[i].fZ - meanZ)*(trigger.fHit[i].fZ - meanZ);
			vSSzx += (trigger.fHit[i].fZ - meanZ)*(trigger.fHit[i].fX - meanX);
			vSSzy += (trigger.fHit[i].fZ - meanZ)*(trigger.fHit[i].fY - meanY);
		}
	}
	
	// Calculate params for lines x = fgMzx * z + fgCzx and y = fgMzy * z + fgCzy.
	if (vSSzz == 0) return false;
	fgMzx = vSSzx / vSSzz;
	fgMzy = vSSzy / vSSzz;
	fgCzx = meanX - fgMzx * meanZ;
	fgCzy = meanY - fgMzy * meanZ;
	
	// Calculate ideal points on chambers 11 and 13:
	fgIdealX1 = fgMzx * fgIdealZ1 + fgCzx;
	fgIdealY1 = fgMzy * fgIdealZ1 + fgCzy;
	fgIdealX2 = fgMzx * fgIdealZ2 + fgCzx;
	fgIdealY2 = fgMzy * fgIdealZ2 + fgCzy;
	
	return true;
}


bool AliHLTMUONCalculations::FitLineToData(
		const AliHLTFloat32_t* x, const AliHLTFloat32_t* y,
		const AliHLTFloat32_t* z, AliHLTUInt32_t n
	)
{
	/// Straight lines are fitted in the ZX and ZY planes using a least
	/// squares fit for the (x, y, z) data points.
	/// http://mathworld.wolfram.com/LeastSquaresFitting.html
	/// If this method returns true, then the fitted parameters can fetched
	/// using the method calls Mzx(), Mzy(), Czx() and Czy(). The lines are
	/// then given by: x = Mzx() * z + Czx() and y = Mzy() * z + Czy()
	/// \param x  This must point to the array of x data values.
	/// \param y  This must point to the array of y data values.
	/// \param z  This must point to the array of z data values.
	/// \param n  Specifies the number of data points in the x, y and z arrays.
	/// \return  true if the line could be fitted or false otherwise.
	///     The reason for failure could be either too few data points or the
	///     slopes Mzx() or Mzy() would be infinite, implying a line that is
	///     perpendicular to the z axis.
	
	if (n < 2) return false;
	
	AliHLTFloat32_t sumX = 0;
	AliHLTFloat32_t sumY = 0;
	AliHLTFloat32_t sumZ = 0;
	for (AliHLTUInt32_t i = 0; i < n; i++)
	{
		sumX += x[i];
		sumY += y[i];
		sumZ += z[i];
	}
	AliHLTFloat32_t meanX = sumX / AliHLTFloat32_t(n);
	AliHLTFloat32_t meanY = sumY / AliHLTFloat32_t(n);
	AliHLTFloat32_t meanZ = sumZ / AliHLTFloat32_t(n);
	
	AliHLTFloat32_t vSSzz = 0;
	AliHLTFloat32_t vSSzx = 0;
	AliHLTFloat32_t vSSzy = 0;
	for (AliHLTUInt32_t i = 0; i < n; i++)
	{
		vSSzz += (z[i] - meanZ)*(z[i] - meanZ);
		vSSzx += (z[i] - meanZ)*(x[i] - meanX);
		vSSzy += (z[i] - meanZ)*(y[i] - meanY);
	}
	
	// Calculate params for lines x = fgMzx * z + fgCzx and y = fgMzy * z + fgCzy.
	if (vSSzz == 0) return false;
	fgMzx = vSSzx / vSSzz;
	fgMzy = vSSzy / vSSzz;
	fgCzx = meanX - fgMzx * meanZ;
	fgCzy = meanY - fgMzy * meanZ;
	
	return true;
}


bool AliHLTMUONCalculations::FitLineToData(
		const AliHLTFloat32_t* x, const AliHLTFloat32_t* z, AliHLTUInt32_t n
	)
{
	/// A straight line is fitted in the X, Z data points using a least squares fit.
	/// http://mathworld.wolfram.com/LeastSquaresFitting.html
	/// If this method returns true, then the fitted parameters can fetched using the
	/// method calls Mzx() and Czx(). The line is then given by: x = Mzx() * z + Czx()
	/// \param x  This must point to the array of x data values.
	/// \param z  This must point to the array of z data values.
	/// \param n  Specifies the number of data points in the x and z arrays.
	/// \return  true if the line could be fitted or false otherwise.
	///     The reason for failure could be either too few data points or the slopes
	///     Mzx() would be infinite, implying a line that is perpendicular to the z axis.
	
	if (n < 2) return false;
	
	AliHLTFloat32_t sumX = 0;
	AliHLTFloat32_t sumZ = 0;
	for (AliHLTUInt32_t i = 0; i < n; i++)
	{
		sumX += x[i];
		sumZ += z[i];
	}
	AliHLTFloat32_t meanX = sumX / AliHLTFloat32_t(n);
	AliHLTFloat32_t meanZ = sumZ / AliHLTFloat32_t(n);
	
	AliHLTFloat32_t vSSzz = 0;
	AliHLTFloat32_t vSSzx = 0;
	for (AliHLTUInt32_t i = 0; i < n; i++)
	{
		vSSzz += (z[i] - meanZ)*(z[i] - meanZ);
		vSSzx += (z[i] - meanZ)*(x[i] - meanX);
	}
	
	// Calculate params for line x = fgMzx * z + fgCzx.
	if (vSSzz == 0) return false;
	fgMzx = vSSzx / vSSzz;
	fgCzx = meanX - fgMzx * meanZ;
	
	return true;
}


AliHLTFloat32_t AliHLTMUONCalculations::AliHLTMUONCalculations::ComputeChi2(
		const AliHLTFloat32_t* x, const AliHLTFloat32_t* y,
		const AliHLTFloat32_t* z, AliHLTUInt32_t n
	)
{
	/// Calculates the chi squared value for the set of data points given
	/// the fitted slope and coefficient parameters previously fitted by
	/// LineFit(x, y, z, n);
	/// The fgSigmaX2 and fgSigmaY2 are used as the variance for the X and
	/// Y coordinates respectively. Note we assume that the covariance terms
	/// are zero.
	/// \param x  This must point to the array of x data values.
	/// \param y  This must point to the array of y data values.
	/// \param z  This must point to the array of z data values.
	/// \param n  Specifies the number of data points in the x, y and z arrays.
	/// \return  The chi squared value.
	
	AliHLTFloat32_t chi2 = 0;
	for (AliHLTUInt32_t i = 0; i < n; i++)
	{
		AliHLTFloat32_t residualX = fgMzx * z[i] + fgCzx - x[i];
		AliHLTFloat32_t residualY = fgMzy * z[i] + fgCzy - y[i];
		chi2 += residualX*residualX/fgSigmaX2 + residualY*residualY/fgSigmaY2;
	}
	return chi2;
}
Commit	Line	Data
364df03a	1	/**************************************************************************
	2	* Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
1d8ae082	16	// $Id$
364df03a	17
	18	////////////////////////////////////////////////////////////////////////////////
	19	//
	20	// Author: Artur Szostak
	21	// Email: artur@alice.phy.uct.ac.za \| artursz@iafrica.com
	22	//
	23	////////////////////////////////////////////////////////////////////////////////
	24
	25	#include "AliHLTMUONCalculations.h"
a3d4b6ba	26	#include "AliHLTMUONUtils.h"
a3d4b6ba	27	#include "AliHLTMUONTriggerRecordsBlockStruct.h"
364df03a	28	#include <cmath>
364df03a	29
b92524d0	30	AliHLTFloat32_t AliHLTMUONCalculations::fgZf = -975.0; // cm
	31
	32	AliHLTFloat32_t AliHLTMUONCalculations::fgQBLScaled
	33	= 3.0 * 2.99792458e8 / 1e9; // T.m.*c/1e9
	34
e6357f88	35	AliHLTMUONParticleSign AliHLTMUONCalculations::fgSign = kSignUnknown;
a3d4b6ba	36	AliHLTFloat32_t AliHLTMUONCalculations::fgPx = 0; // GeV/c
	37	AliHLTFloat32_t AliHLTMUONCalculations::fgPy = 0; // GeV/c
	38	AliHLTFloat32_t AliHLTMUONCalculations::fgPz = 0; // GeV/c
	39
	40	AliHLTFloat32_t AliHLTMUONCalculations::fgSigmaX2 = 1.; // cm^2
	41	AliHLTFloat32_t AliHLTMUONCalculations::fgSigmaY2 = 1.; // cm^2
	42
	43	AliHLTFloat32_t AliHLTMUONCalculations::fgMzx = 0;
	44	AliHLTFloat32_t AliHLTMUONCalculations::fgMzy = 0;
	45	AliHLTFloat32_t AliHLTMUONCalculations::fgCzx = 0;
	46	AliHLTFloat32_t AliHLTMUONCalculations::fgCzy = 0;
	47
	48	AliHLTFloat32_t AliHLTMUONCalculations::fgIdealX1 = 0; // cm
	49	AliHLTFloat32_t AliHLTMUONCalculations::fgIdealY1 = 0; // cm
	50	AliHLTFloat32_t AliHLTMUONCalculations::fgIdealZ1 = -1603.5f; // cm
	51	AliHLTFloat32_t AliHLTMUONCalculations::fgIdealX2 = 0; // cm
	52	AliHLTFloat32_t AliHLTMUONCalculations::fgIdealY2 = 0; // cm
	53	AliHLTFloat32_t AliHLTMUONCalculations::fgIdealZ2 = -1703.5f; // cm
e6357f88	54
	55
	56	bool AliHLTMUONCalculations::ComputeMomentum(
	57	AliHLTFloat32_t x1,
	58	AliHLTFloat32_t y1, AliHLTFloat32_t y2,
	59	AliHLTFloat32_t z1, AliHLTFloat32_t z2
	60	)
	61	{
a3d4b6ba	62	/// Computes the momentum components based on the equations given in the
	63	/// ALICE dimuon spectrometer Technical Design Report (TDR-5): trigger section.
	64	///
	65	/// Reference:
	66	/// "CERN/LHCC 2000-046
	67	/// Addendum 1 to ALICE TDR 5
	68	/// 15 Dec 2000"
	69	/// Section 3.1.2 pages 144 and 145.
	70	///
	71	/// Input can be in meters, cm or mm. Output is in GeV/c.
	72	///
	73	/// \param x1 X coordinate of hit point 1 on the track.
	74	/// \param y1 Y coordinate of hit point 1 on the track.
	75	/// \param z1 Z coordinate of hit point 1 on the track.
	76	/// \param y2 Y coordinate of hit point 2 on the track.
	77	/// \param z2 Z coordinate of hit point 2 on the track.
	78	/// \return true if the momentum could be calculated and false otherwise.
	79	/// If true is returned then the estimated momentum can be fetched by the
	80	/// method calls: Px(), Py() and Pz() for the individual components.
	81
e6357f88	82	AliHLTFloat64_t z2mz1 = z2 - z1;
	83	if (z2mz1 == 0 or z1 == 0)
	84	{
	85	fgSign = kSignUnknown;
	86	fgPx = fgPy = fgPz = 0;
	87	return false;
	88	}
	89	AliHLTFloat64_t thetaTimesZf = (y1z2 - y2z1) / z2mz1;
	90	AliHLTFloat64_t xf = x1 * fgZf / z1;
	91	AliHLTFloat64_t yf = y2 - ((y2-y1) * (z2-fgZf)) / z2mz1;
	92
	93	if (thetaTimesZf == 0)
	94	{
	95	fgSign = kSignUnknown;
	96	fgPx = fgPy = fgPz = 0;
	97	return false;
	98	}
b92524d0	99	AliHLTFloat64_t pDivZf = (fgQBLScaled / thetaTimesZf);
e6357f88	100	AliHLTFloat64_t p = pDivZf * fgZf;
	101	pDivZf = fabs(pDivZf);
	102
	103	if (p < 0)
	104	fgSign = kSignMinus;
	105	else if (p > 0)
	106	fgSign = kSignPlus;
	107	else
	108	fgSign = kSignUnknown;
	109
	110	fgPx = AliHLTFloat32_t( pDivZf * xf );
	111	fgPy = AliHLTFloat32_t( pDivZf * yf );
4a9f11d4	112	AliHLTFloat64_t k = pp - fgPxfgPx - fgPy*fgPy;
	113	if (k > 0)
	114	fgPz = AliHLTFloat32_t( sqrt(k) );
	115	else
	116	fgPz = 0;
e6357f88	117	// fgPz must be the same sign as fgZf else it could not have been measured.
	118	if (fgZf < 0) fgPz = -fgPz;
	119
	120	return true;
	121	}
	122
364df03a	123
b92524d0	124	AliHLTFloat32_t AliHLTMUONCalculations::QBL()
364df03a	125	{
b92524d0	126	// We have to convert back into Tesla metres.
	127	return fgQBLScaled * 1e9 / 2.99792458e8;
	128	}
364df03a	129
364df03a	130
b92524d0	131	void AliHLTMUONCalculations::QBL(AliHLTFloat32_t value)
364df03a	132	{
b92524d0	133	// Note: 2.99792458e8/1e9 is the conversion factor for GeV.
	134	// It is c/1e9, where c is the speed of light.
	135	fgQBLScaled = value * 2.99792458e8 / 1e9;
	136	}
c9537879	137
	138
	139	AliHLTFloat32_t AliHLTMUONCalculations::ComputeMass(
	140	AliHLTFloat32_t massA,
	141	AliHLTFloat32_t pxA,
	142	AliHLTFloat32_t pyA,
	143	AliHLTFloat32_t pzA,
	144	AliHLTFloat32_t massB,
	145	AliHLTFloat32_t pxB,
	146	AliHLTFloat32_t pyB,
	147	AliHLTFloat32_t pzB
	148	)
	149	{
	150	/// Calculates the invariant mass for a pair of particles.
	151	/// \param massA Mmass in GeV/c of particle A.
	152	/// \param pxA X component of the momentum in GeV/c for particle A.
	153	/// \param pyA Y component of the momentum in GeV/c for particle A.
	154	/// \param pzA Z component of the momentum in GeV/c for particle A.
	155	/// \param massB Mass in GeV/c of particle B.
	156	/// \param pxB X component of the momentum in GeV/c for particle B.
	157	/// \param pyB Y component of the momentum in GeV/c for particle B.
	158	/// \param pzB Z component of the momentum in GeV/c for particle B.
	159	/// \return The invariant mass in GeV/c^2 or -1 if there was a problem
	160	/// in the calculation due to bad input parameters.
	161
	162	AliHLTFloat32_t massA2 = massA*massA;
	163	AliHLTFloat32_t massB2 = massB*massB;
	164	AliHLTFloat32_t energyA = sqrt(massA2 + pxApxA + pyApyA + pzA*pzA);
	165	AliHLTFloat32_t energyB = sqrt(massB2 + pxBpxB + pyBpyB + pzB*pzB);
	166	AliHLTFloat32_t mass2 = massA2 + massB2 + 2. * (energyAenergyB - pxApxB - pyApyB - pzApzB);
	167	if (mass2 < 0.) return -1.;
	168	return sqrt(mass2);
	169	}
	170
a3d4b6ba	171
	172	bool AliHLTMUONCalculations::FitLineToTriggerRecord(
	173	const AliHLTMUONTriggerRecordStruct& trigger
	174	)
	175	{
	176	/// Straight lines are fitted in the ZX and ZY planes using a least
	177	/// squares fit to the coordinates in the trigger record.
	178	/// http://mathworld.wolfram.com/LeastSquaresFitting.html
	179	/// If this method returns true, then the fitted parameters can fetched
	180	/// using the method calls Mzx(), Mzy(), Czx() and Czy(). The lines are
	181	/// then given by: x = Mzx() * z + Czx() and y = Mzy() * z + Czy()
	182	/// The ideal coordinates are also calculated and can be fetched with
	183	/// the method calls: IdealX1(), IdealY1() and IdealZ1() for point on MT1,
	184	/// and IdealX2(), IdealY2() and IdealZ2() for point on MT2.
	185	/// \param trigger The trigger record structure to which we fit a line.
	186	/// \return true if the line could be fitted or false otherwise.
	187	/// The reason for failure could be either too few hits or the slopes
	188	/// Mzx() or Mzy() would be infinite, implying a line that is
	189	/// perpendicular to the z axis.
	190
	191	AliHLTMUONParticleSign sign;
	192	bool hitset[4];
	193	AliHLTMUONUtils::UnpackTriggerRecordFlags(trigger.fFlags, sign, hitset);
	194	DebugTrace("hitset = {" << hitset[0] << ", " << hitset[1] << ", "
	195	<< hitset[2] << ", " << hitset[3] << "}"
	196	);
	197
	198	return FitLineToTriggerRecord(trigger, hitset);
	199	}
	200
	201
	202	bool AliHLTMUONCalculations::FitLineToTriggerRecord(
	203	const AliHLTMUONTriggerRecordStruct& trigger,
	204	const bool hitset[4]
	205	)
	206	{
	207	/// Performs a straight line fit like FitLineToTriggerRecord(trigger)
	208	/// but requires pree-decoded flags indicating which hits were set.
	209	/// \param trigger The trigger record structure to which we fit a line.
	210	/// \param hitset Flags indicating which hits were set in the trigger record.
	211	/// \return true if the line could be fitted or false otherwise.
	212
	213	AliHLTFloat32_t sumX = 0;
	214	AliHLTFloat32_t sumY = 0;
	215	AliHLTFloat32_t sumZ = 0;
	216	int n = 0;
	217	for (int i = 0; i < 4; i++)
	218	{
	219	if (hitset[i])
	220	{
	221	sumX += trigger.fHit[i].fX;
	222	sumY += trigger.fHit[i].fY;
	223	sumZ += trigger.fHit[i].fZ;
	224	n++;
	225	}
	226	}
	227	if (n < 2) return false;
	228	AliHLTFloat32_t meanX = sumX / AliHLTFloat32_t(n);
	229	AliHLTFloat32_t meanY = sumY / AliHLTFloat32_t(n);
	230	AliHLTFloat32_t meanZ = sumZ / AliHLTFloat32_t(n);
	231
	232	AliHLTFloat32_t vSSzz = 0;
	233	AliHLTFloat32_t vSSzx = 0;
	234	AliHLTFloat32_t vSSzy = 0;
235	for (int i = 0; i < 4; i++)
236	{
237	if (hitset[i])
238	{
239	vSSzz += (trigger.fHit[i].fZ - meanZ)*(trigger.fHit[i].fZ - meanZ);
240	vSSzx += (trigger.fHit[i].fZ - meanZ)*(trigger.fHit[i].fX - meanX);
241	vSSzy += (trigger.fHit[i].fZ - meanZ)*(trigger.fHit[i].fY - meanY);
242	}
243	}
244
245	// Calculate params for lines x = fgMzx * z + fgCzx and y = fgMzy * z + fgCzy.
246	if (vSSzz == 0) return false;
247	fgMzx = vSSzx / vSSzz;
248	fgMzy = vSSzy / vSSzz;
249	fgCzx = meanX - fgMzx * meanZ;
250	fgCzy = meanY - fgMzy * meanZ;
251
252	// Calculate ideal points on chambers 11 and 13:
253	fgIdealX1 = fgMzx * fgIdealZ1 + fgCzx;
254	fgIdealY1 = fgMzy * fgIdealZ1 + fgCzy;
255	fgIdealX2 = fgMzx * fgIdealZ2 + fgCzx;
256	fgIdealY2 = fgMzy * fgIdealZ2 + fgCzy;
257
258	return true;
259	}
260
261
262	bool AliHLTMUONCalculations::FitLineToData(
263	const AliHLTFloat32_t* x, const AliHLTFloat32_t* y,
264	const AliHLTFloat32_t* z, AliHLTUInt32_t n
265	)
266	{
267	/// Straight lines are fitted in the ZX and ZY planes using a least
268	/// squares fit for the (x, y, z) data points.
269	/// http://mathworld.wolfram.com/LeastSquaresFitting.html
270	/// If this method returns true, then the fitted parameters can fetched
271	/// using the method calls Mzx(), Mzy(), Czx() and Czy(). The lines are
272	/// then given by: x = Mzx() * z + Czx() and y = Mzy() * z + Czy()
273	/// \param x This must point to the array of x data values.
274	/// \param y This must point to the array of y data values.
275	/// \param z This must point to the array of z data values.
276	/// \param n Specifies the number of data points in the x, y and z arrays.
277	/// \return true if the line could be fitted or false otherwise.
278	/// The reason for failure could be either too few data points or the
279	/// slopes Mzx() or Mzy() would be infinite, implying a line that is
280	/// perpendicular to the z axis.
281
282	if (n < 2) return false;
283
284	AliHLTFloat32_t sumX = 0;
285	AliHLTFloat32_t sumY = 0;
286	AliHLTFloat32_t sumZ = 0;
287	for (AliHLTUInt32_t i = 0; i < n; i++)
288	{
289	sumX += x[i];
290	sumY += y[i];
291	sumZ += z[i];
292	}
293	AliHLTFloat32_t meanX = sumX / AliHLTFloat32_t(n);
294	AliHLTFloat32_t meanY = sumY / AliHLTFloat32_t(n);
295	AliHLTFloat32_t meanZ = sumZ / AliHLTFloat32_t(n);
296
297	AliHLTFloat32_t vSSzz = 0;
298	AliHLTFloat32_t vSSzx = 0;
299	AliHLTFloat32_t vSSzy = 0;
300	for (AliHLTUInt32_t i = 0; i < n; i++)
301	{
302	vSSzz += (z[i] - meanZ)*(z[i] - meanZ);
303	vSSzx += (z[i] - meanZ)*(x[i] - meanX);
304	vSSzy += (z[i] - meanZ)*(y[i] - meanY);
305	}
306
307	// Calculate params for lines x = fgMzx * z + fgCzx and y = fgMzy * z + fgCzy.
308	if (vSSzz == 0) return false;
309	fgMzx = vSSzx / vSSzz;
310	fgMzy = vSSzy / vSSzz;
311	fgCzx = meanX - fgMzx * meanZ;
312	fgCzy = meanY - fgMzy * meanZ;
313
314	return true;
315	}
316
317
462e3880	318	bool AliHLTMUONCalculations::FitLineToData(
	319	const AliHLTFloat32_t* x, const AliHLTFloat32_t* z, AliHLTUInt32_t n
	320	)
	321	{
	322	/// A straight line is fitted in the X, Z data points using a least squares fit.
	323	/// http://mathworld.wolfram.com/LeastSquaresFitting.html
	324	/// If this method returns true, then the fitted parameters can fetched using the
	325	/// method calls Mzx() and Czx(). The line is then given by: x = Mzx() * z + Czx()
	326	/// \param x This must point to the array of x data values.
	327	/// \param z This must point to the array of z data values.
	328	/// \param n Specifies the number of data points in the x and z arrays.
	329	/// \return true if the line could be fitted or false otherwise.
	330	/// The reason for failure could be either too few data points or the slopes
	331	/// Mzx() would be infinite, implying a line that is perpendicular to the z axis.
	332
	333	if (n < 2) return false;
	334
	335	AliHLTFloat32_t sumX = 0;
	336	AliHLTFloat32_t sumZ = 0;
	337	for (AliHLTUInt32_t i = 0; i < n; i++)
	338	{
	339	sumX += x[i];
	340	sumZ += z[i];
	341	}
	342	AliHLTFloat32_t meanX = sumX / AliHLTFloat32_t(n);
	343	AliHLTFloat32_t meanZ = sumZ / AliHLTFloat32_t(n);
	344
	345	AliHLTFloat32_t vSSzz = 0;
	346	AliHLTFloat32_t vSSzx = 0;
	347	for (AliHLTUInt32_t i = 0; i < n; i++)
	348	{
	349	vSSzz += (z[i] - meanZ)*(z[i] - meanZ);
	350	vSSzx += (z[i] - meanZ)*(x[i] - meanX);
	351	}
	352
	353	// Calculate params for line x = fgMzx * z + fgCzx.
	354	if (vSSzz == 0) return false;
	355	fgMzx = vSSzx / vSSzz;
	356	fgCzx = meanX - fgMzx * meanZ;
	357
	358	return true;
	359	}
	360
	361
a3d4b6ba	362	AliHLTFloat32_t AliHLTMUONCalculations::AliHLTMUONCalculations::ComputeChi2(
	363	const AliHLTFloat32_t* x, const AliHLTFloat32_t* y,
	364	const AliHLTFloat32_t* z, AliHLTUInt32_t n
	365	)
	366	{
	367	/// Calculates the chi squared value for the set of data points given
	368	/// the fitted slope and coefficient parameters previously fitted by
	369	/// LineFit(x, y, z, n);
	370	/// The fgSigmaX2 and fgSigmaY2 are used as the variance for the X and
	371	/// Y coordinates respectively. Note we assume that the covariance terms
	372	/// are zero.
	373	/// \param x This must point to the array of x data values.
	374	/// \param y This must point to the array of y data values.
	375	/// \param z This must point to the array of z data values.
	376	/// \param n Specifies the number of data points in the x, y and z arrays.
	377	/// \return The chi squared value.
	378
	379	AliHLTFloat32_t chi2 = 0;
	380	for (AliHLTUInt32_t i = 0; i < n; i++)
	381	{
	382	AliHLTFloat32_t residualX = fgMzx * z[i] + fgCzx - x[i];
	383	AliHLTFloat32_t residualY = fgMzy * z[i] + fgCzy - y[i];
	384	chi2 += residualXresidualX/fgSigmaX2 + residualYresidualY/fgSigmaY2;
	385	}
	386	return chi2;
	387	}