[u/mrichter/AliRoot.git] / TRD / AliTRDltuParam.cxx

#include <stdio.h>

#include "TMath.h"

#include "AliTRDltuParam.h"

// definition of geometry constants
Float_t AliTRDltuParam::fgZrow[6][5] = {
  {301, 177, 53, -57, -181},
  {301, 177, 53, -57, -181},
  {315, 184, 53, -57, -188},
  {329, 191, 53, -57, -195},
  {343, 198, 53, -57, -202},
  {347, 200, 53, -57, -204}};
Float_t AliTRDltuParam::fgX[6] =
  {300.65, 313.25, 325.85, 338.45, 351.05, 363.65};
Float_t AliTRDltuParam::fgTiltingAngle[6] =
  {-2., 2., -2., 2., -2., 2.};
Int_t   AliTRDltuParam::fgDyMax =  63;
Int_t   AliTRDltuParam::fgDyMin = -64;
Float_t AliTRDltuParam::fgBinDy = 140e-4;
Float_t AliTRDltuParam::fgWidthPad[6] =
  {0.635, 0.665, 0.695, 0.725, 0.755, 0.785};
Float_t AliTRDltuParam::fgLengthInnerPadC1[6] =
  {7.5, 7.5, 8.0, 8.5, 9.0, 9.0};
Float_t AliTRDltuParam::fgLengthOuterPadC1[6] =
  {7.5, 7.5, 7.5, 7.5, 7.5, 8.5};
Float_t AliTRDltuParam::fgLengthInnerPadC0 = 9.0;
Float_t AliTRDltuParam::fgLengthOuterPadC0 = 8.0;
Float_t AliTRDltuParam::fgScalePad = 256. * 32.;
Float_t AliTRDltuParam::fgDriftLength = 3.;

AliTRDltuParam::AliTRDltuParam() :
  TObject(),
  fMagField(0.),
  fOmegaTau(0.),
  fPtMin(0.1),
  fNtimebins(20 << 5),
  fScaleQ0(0),
  fScaleQ1(0),
  fPidTracklengthCorr(kFALSE),
  fTiltCorr(kFALSE),
  fPidGainCorr(kFALSE)
{
  // default constructor
}

AliTRDltuParam::~AliTRDltuParam()
{
  // destructor
}

Int_t AliTRDltuParam::GetDyCorrection(Int_t det, Int_t rob, Int_t mcm) const
{
  // calculate the correction of the deflection
  // i.e. Lorentz angle and tilt correction (if active)

  Int_t layer = det % 6;

  Float_t dyTilt = ( fgDriftLength * TMath::Tan(fgTiltingAngle[layer] * TMath::Pi()/180.) *
  		     GetLocalZ(det, rob, mcm) / fgX[layer] );

  // calculate Lorentz correction
  Float_t dyCorr = - fOmegaTau * fgDriftLength;

  if(fTiltCorr)
    dyCorr += dyTilt; // add tilt correction

  return (int) TMath::Nint(dyCorr * fgScalePad / fgWidthPad[layer]);
}

void AliTRDltuParam::GetDyRange(Int_t det, Int_t rob, Int_t mcm, Int_t ch,
				 Int_t &dyMinInt, Int_t &dyMaxInt) const
{
  // calculate the deflection range in which tracklets are accepted

  dyMinInt = fgDyMin;
  dyMaxInt = fgDyMax;

  // deflection cut is considered for |B| > 0.1 T only
  if (TMath::Abs(fMagField) < 0.1)
    return;

  Float_t e = 0.30;

  Float_t maxDeflTemp = GetPerp(det, rob, mcm, ch)/2. * // Sekante/2 (cm)
    (e * 1e-2 * TMath::Abs(fMagField) / fPtMin);   // 1/R (1/cm)

  Float_t maxDeflAngle = 0.;

  Float_t phi = GetPhi(det, rob, mcm, ch);
  if (maxDeflTemp < TMath::Cos(phi)) {
    maxDeflAngle = TMath::ASin(maxDeflTemp);

    Float_t dyMin = ( fgDriftLength *
		      TMath::Tan(phi - maxDeflAngle) );

    dyMinInt = Int_t(dyMin / fgBinDy);
    // clipping to allowed range
    if (dyMinInt < fgDyMin)
      dyMinInt = fgDyMin;
    else if (dyMinInt > fgDyMax)
      dyMinInt = fgDyMax;

    Float_t dyMax = ( fgDriftLength *
  		      TMath::Tan(phi + maxDeflAngle) );

    dyMaxInt = Int_t(dyMax / fgBinDy);
    // clipping to allowed range
    if (dyMaxInt > fgDyMax)
      dyMaxInt = fgDyMax;
    else if (dyMaxInt < fgDyMin)
      dyMaxInt = fgDyMin;
  }
  else if (maxDeflTemp < 0.) {
    // this must not happen
    printf("Inconsistent calculation of sin(alpha): %f\n", maxDeflTemp);
  }
  else {
    // TRD is not reached at the given pt threshold
    // max range
  }

  if ((dyMaxInt - dyMinInt) <= 0) {
    printf("strange dy range: [%i,%i], using max range now\n", dyMinInt, dyMaxInt);
    dyMaxInt = fgDyMax;
    dyMinInt = fgDyMin;
  }
}

Float_t AliTRDltuParam::GetElongation(Int_t det, Int_t rob, Int_t mcm, Int_t ch) const
{
  // calculate the ratio of the distance to the primary vertex and the
  // distance in x-direction for the given ADC channel

  Int_t layer = det % 6;

  Float_t elongation = TMath::Abs(GetDist(det, rob, mcm, ch) / fgX[layer]);

  // sanity check
  if(elongation<0.001) {
    elongation=1.;
  }
  return elongation;
}

void AliTRDltuParam::GetCorrectionFactors(Int_t det, Int_t rob, Int_t mcm, Int_t ch,
					  UInt_t &cor0, UInt_t &cor1, Float_t gain) const
{
  // calculate the gain correction factors for the given ADC channel

  if (fPidGainCorr==kFALSE)
    gain=1;

  if (fPidTracklengthCorr == kTRUE ) {
    cor0 = UInt_t ((1.0*fScaleQ0* (1./GetElongation(det, rob, mcm, ch))) / gain );
    cor1 = UInt_t ((1.0*fScaleQ1* (1./GetElongation(det, rob, mcm, ch))) / gain );
  }
  else {
    cor0 = UInt_t (fScaleQ0 / gain);
    cor1 = UInt_t ( fScaleQ1 / gain);
  }
}

Int_t AliTRDltuParam::GetNtimebins() const
{
  // return the number of timebins used

  return fNtimebins;
}

Float_t AliTRDltuParam::GetX(Int_t det, Int_t /* rob */, Int_t /* mcm */) const
{
  // return the distance to the beam axis in x-direction

  Int_t layer = det%6;
  return fgX[layer];
}

Float_t AliTRDltuParam::GetLocalY(Int_t det, Int_t rob, Int_t mcm, Int_t ch) const
{
  // get local y-position (r-phi) w.r.t. the chamber centre

  Int_t layer = det%6;
  // calculate the pad position as in the TRAP
  Float_t ypos = (-4 + 1 + (rob&0x1) * 4 + (mcm&0x3)) * 18 - ch - 0.5; // y position in bins of pad widths
  return ypos*fgWidthPad[layer];
}

Float_t AliTRDltuParam::GetLocalZ(Int_t det, Int_t rob, Int_t mcm) const
{
  // get local z-position w.r.t. to the chamber boundary

  Int_t stack = (det%30) / 6;
  Int_t layer = det % 6;
  Int_t row   = (rob/2) * 4 + mcm/4;

  if (stack == 2) {
    if (row == 0)
      return (fgZrow[layer][stack] - 0.5 * fgLengthOuterPadC0);
    else if (row == 11)
      return (fgZrow[layer][stack] - 1.5 * fgLengthOuterPadC0 - (row - 1) * fgLengthInnerPadC0);
    else
      return (fgZrow[layer][stack] - fgLengthOuterPadC0 - (row - 0.5) * fgLengthInnerPadC0);
  }
  else {
    if (row == 0)
      return (fgZrow[layer][stack] - 0.5 * fgLengthOuterPadC1[layer]);
    else if (row == 15)
      return (fgZrow[layer][stack] - 1.5 * fgLengthOuterPadC1[layer] - (row - 1) * fgLengthInnerPadC1[layer]);
    else
      return (fgZrow[layer][stack] - fgLengthOuterPadC1[layer] - (row - 0.5) * fgLengthInnerPadC1[layer]);
  }
}

Float_t AliTRDltuParam::GetPerp(Int_t det, Int_t rob, Int_t mcm, Int_t ch) const
{
  // get transverse distance to the beam axis

  return TMath::Sqrt(GetLocalY(det, rob, mcm, ch)*GetLocalY(det, rob, mcm, ch) +
		     GetX(det, rob, mcm)*GetX(det, rob, mcm) );
}

Float_t AliTRDltuParam::GetPhi(Int_t det, Int_t rob, Int_t mcm, Int_t ch) const
{
  // calculate the azimuthal angle for the given ADC channel

  return TMath::ATan2(GetLocalY(det, rob, mcm, ch), GetX(det, rob, mcm));
}

Float_t AliTRDltuParam::GetDist(Int_t det, Int_t rob, Int_t mcm, Int_t ch) const
{
  // calculate the distance from the origin for the given ADC channel

  return TMath::Sqrt(GetLocalY(det, rob, mcm, ch)*GetLocalY(det, rob, mcm, ch) +
		     GetX(det, rob, mcm)*GetX(det, rob, mcm) +
		     GetLocalZ(det, rob, mcm)*GetLocalZ(det, rob, mcm) );
}
Commit	Line	Data
5e86ff99	1	#include <stdio.h>
	2
	3	#include "TMath.h"
	4
	5	#include "AliTRDltuParam.h"
	6
	7	// definition of geometry constants
	8	Float_t AliTRDltuParam::fgZrow[6][5] = {
	9	{301, 177, 53, -57, -181},
	10	{301, 177, 53, -57, -181},
	11	{315, 184, 53, -57, -188},
	12	{329, 191, 53, -57, -195},
	13	{343, 198, 53, -57, -202},
	14	{347, 200, 53, -57, -204}};
	15	Float_t AliTRDltuParam::fgX[6] =
	16	{300.65, 313.25, 325.85, 338.45, 351.05, 363.65};
	17	Float_t AliTRDltuParam::fgTiltingAngle[6] =
	18	{-2., 2., -2., 2., -2., 2.};
	19	Int_t AliTRDltuParam::fgDyMax = 63;
	20	Int_t AliTRDltuParam::fgDyMin = -64;
1fa993c2	21	Float_t AliTRDltuParam::fgBinDy = 140e-4;
5e86ff99	22	Float_t AliTRDltuParam::fgWidthPad[6] =
	23	{0.635, 0.665, 0.695, 0.725, 0.755, 0.785};
	24	Float_t AliTRDltuParam::fgLengthInnerPadC1[6] =
	25	{7.5, 7.5, 8.0, 8.5, 9.0, 9.0};
	26	Float_t AliTRDltuParam::fgLengthOuterPadC1[6] =
	27	{7.5, 7.5, 7.5, 7.5, 7.5, 8.5};
	28	Float_t AliTRDltuParam::fgLengthInnerPadC0 = 9.0;
	29	Float_t AliTRDltuParam::fgLengthOuterPadC0 = 8.0;
	30	Float_t AliTRDltuParam::fgScalePad = 256. * 32.;
918cceaa	31	Float_t AliTRDltuParam::fgDriftLength = 3.;
5e86ff99	32
	33	AliTRDltuParam::AliTRDltuParam() :
	34	TObject(),
	35	fMagField(0.),
	36	fOmegaTau(0.),
	37	fPtMin(0.1),
	38	fNtimebins(20 << 5),
51380642	39	fScaleQ0(0),
51380642	40	fScaleQ1(0),
5e86ff99	41	fPidTracklengthCorr(kFALSE),
5f7c3c48	42	fTiltCorr(kFALSE),
5f7c3c48	43	fPidGainCorr(kFALSE)
5e86ff99	44	{
6419bebb	45	// default constructor
5e86ff99	46	}
	47
	48	AliTRDltuParam::~AliTRDltuParam()
	49	{
6419bebb	50	// destructor
5e86ff99	51	}
	52
	53	Int_t AliTRDltuParam::GetDyCorrection(Int_t det, Int_t rob, Int_t mcm) const
	54	{
	55	// calculate the correction of the deflection
	56	// i.e. Lorentz angle and tilt correction (if active)
	57
	58	Int_t layer = det % 6;
	59
16477104	60	Float_t dyTilt = ( fgDriftLength * TMath::Tan(fgTiltingAngle[layer] * TMath::Pi()/180.) *
5e86ff99	61	GetLocalZ(det, rob, mcm) / fgX[layer] );
	62
	63	// calculate Lorentz correction
	64	Float_t dyCorr = - fOmegaTau * fgDriftLength;
	65
	66	if(fTiltCorr)
	67	dyCorr += dyTilt; // add tilt correction
	68
	69	return (int) TMath::Nint(dyCorr * fgScalePad / fgWidthPad[layer]);
	70	}
	71
	72	void AliTRDltuParam::GetDyRange(Int_t det, Int_t rob, Int_t mcm, Int_t ch,
	73	Int_t &dyMinInt, Int_t &dyMaxInt) const
	74	{
	75	// calculate the deflection range in which tracklets are accepted
	76
	77	dyMinInt = fgDyMin;
	78	dyMaxInt = fgDyMax;
	79
1fa993c2	80	// deflection cut is considered for \|B\| > 0.1 T only
39e4e4cf	81	if (TMath::Abs(fMagField) < 0.1)
5e86ff99	82	return;
	83
	84	Float_t e = 0.30;
	85
1fa993c2	86	Float_t maxDeflTemp = GetPerp(det, rob, mcm, ch)/2. * // Sekante/2 (cm)
1fa993c2	87	(e * 1e-2 * TMath::Abs(fMagField) / fPtMin); // 1/R (1/cm)
5e86ff99	88
	89	Float_t maxDeflAngle = 0.;
	90
1fa993c2	91	Float_t phi = GetPhi(det, rob, mcm, ch);
1fa993c2	92	if (maxDeflTemp < TMath::Cos(phi)) {
5e86ff99	93	maxDeflAngle = TMath::ASin(maxDeflTemp);
	94
	95	Float_t dyMin = ( fgDriftLength *
1fa993c2	96	TMath::Tan(phi - maxDeflAngle) );
5e86ff99	97
5e86ff99	98	dyMinInt = Int_t(dyMin / fgBinDy);
1fa993c2	99	// clipping to allowed range
5e86ff99	100	if (dyMinInt < fgDyMin)
5e86ff99	101	dyMinInt = fgDyMin;
1fa993c2	102	else if (dyMinInt > fgDyMax)
1fa993c2	103	dyMinInt = fgDyMax;
5e86ff99	104
5e86ff99	105	Float_t dyMax = ( fgDriftLength *
1fa993c2	106	TMath::Tan(phi + maxDeflAngle) );
5e86ff99	107
5e86ff99	108	dyMaxInt = Int_t(dyMax / fgBinDy);
1fa993c2	109	// clipping to allowed range
5e86ff99	110	if (dyMaxInt > fgDyMax)
5e86ff99	111	dyMaxInt = fgDyMax;
1fa993c2	112	else if (dyMaxInt < fgDyMin)
1fa993c2	113	dyMaxInt = fgDyMin;
5e86ff99	114	}
1fa993c2	115	else if (maxDeflTemp < 0.) {
	116	// this must not happen
	117	printf("Inconsistent calculation of sin(alpha): %f\n", maxDeflTemp);
	118	}
	119	else {
	120	// TRD is not reached at the given pt threshold
	121	// max range
	122	}
	123
5e86ff99	124	if ((dyMaxInt - dyMinInt) <= 0) {
1fa993c2	125	printf("strange dy range: [%i,%i], using max range now\n", dyMinInt, dyMaxInt);
	126	dyMaxInt = fgDyMax;
	127	dyMinInt = fgDyMin;
5e86ff99	128	}
	129	}
	130
	131	Float_t AliTRDltuParam::GetElongation(Int_t det, Int_t rob, Int_t mcm, Int_t ch) const
	132	{
6419bebb	133	// calculate the ratio of the distance to the primary vertex and the
	134	// distance in x-direction for the given ADC channel
	135
5e86ff99	136	Int_t layer = det % 6;
	137
	138	Float_t elongation = TMath::Abs(GetDist(det, rob, mcm, ch) / fgX[layer]);
	139
	140	// sanity check
	141	if(elongation<0.001) {
	142	elongation=1.;
	143	}
	144	return elongation;
	145	}
	146
	147	void AliTRDltuParam::GetCorrectionFactors(Int_t det, Int_t rob, Int_t mcm, Int_t ch,
5f7c3c48	148	UInt_t &cor0, UInt_t &cor1, Float_t gain) const
5e86ff99	149	{
6419bebb	150	// calculate the gain correction factors for the given ADC channel
6419bebb	151
5f7c3c48	152	if (fPidGainCorr==kFALSE)
	153	gain=1;
	154
5e86ff99	155	if (fPidTracklengthCorr == kTRUE ) {
5f7c3c48	156	cor0 = UInt_t ((1.0fScaleQ0 (1./GetElongation(det, rob, mcm, ch))) / gain );
5f7c3c48	157	cor1 = UInt_t ((1.0fScaleQ1 (1./GetElongation(det, rob, mcm, ch))) / gain );
5e86ff99	158	}
5e86ff99	159	else {
5f7c3c48	160	cor0 = UInt_t (fScaleQ0 / gain);
5f7c3c48	161	cor1 = UInt_t ( fScaleQ1 / gain);
5e86ff99	162	}
	163	}
	164
	165	Int_t AliTRDltuParam::GetNtimebins() const
	166	{
6419bebb	167	// return the number of timebins used
6419bebb	168
5e86ff99	169	return fNtimebins;
	170	}
	171
	172	Float_t AliTRDltuParam::GetX(Int_t det, Int_t /* rob /, Int_t / mcm */) const
	173	{
6419bebb	174	// return the distance to the beam axis in x-direction
6419bebb	175
5e86ff99	176	Int_t layer = det%6;
	177	return fgX[layer];
	178	}
	179
	180	Float_t AliTRDltuParam::GetLocalY(Int_t det, Int_t rob, Int_t mcm, Int_t ch) const
	181	{
6419bebb	182	// get local y-position (r-phi) w.r.t. the chamber centre
6419bebb	183
5e86ff99	184	Int_t layer = det%6;
5e86ff99	185	// calculate the pad position as in the TRAP
39e4e4cf	186	Float_t ypos = (-4 + 1 + (rob&0x1) * 4 + (mcm&0x3)) * 18 - ch - 0.5; // y position in bins of pad widths
5e86ff99	187	return ypos*fgWidthPad[layer];
	188	}
	189
	190	Float_t AliTRDltuParam::GetLocalZ(Int_t det, Int_t rob, Int_t mcm) const
	191	{
6419bebb	192	// get local z-position w.r.t. to the chamber boundary
6419bebb	193
5e86ff99	194	Int_t stack = (det%30) / 6;
	195	Int_t layer = det % 6;
	196	Int_t row = (rob/2) * 4 + mcm/4;
	197
	198	if (stack == 2) {
	199	if (row == 0)
	200	return (fgZrow[layer][stack] - 0.5 * fgLengthOuterPadC0);
	201	else if (row == 11)
	202	return (fgZrow[layer][stack] - 1.5 * fgLengthOuterPadC0 - (row - 1) * fgLengthInnerPadC0);
	203	else
	204	return (fgZrow[layer][stack] - fgLengthOuterPadC0 - (row - 0.5) * fgLengthInnerPadC0);
	205	}
	206	else {
	207	if (row == 0)
	208	return (fgZrow[layer][stack] - 0.5 * fgLengthOuterPadC1[layer]);
	209	else if (row == 15)
	210	return (fgZrow[layer][stack] - 1.5 * fgLengthOuterPadC1[layer] - (row - 1) * fgLengthInnerPadC1[layer]);
	211	else
	212	return (fgZrow[layer][stack] - fgLengthOuterPadC1[layer] - (row - 0.5) * fgLengthInnerPadC1[layer]);
	213	}
	214	}
	215
	216	Float_t AliTRDltuParam::GetPerp(Int_t det, Int_t rob, Int_t mcm, Int_t ch) const
	217	{
6419bebb	218	// get transverse distance to the beam axis
6419bebb	219
5e86ff99	220	return TMath::Sqrt(GetLocalY(det, rob, mcm, ch)*GetLocalY(det, rob, mcm, ch) +
	221	GetX(det, rob, mcm)*GetX(det, rob, mcm) );
	222	}
	223
	224	Float_t AliTRDltuParam::GetPhi(Int_t det, Int_t rob, Int_t mcm, Int_t ch) const
	225	{
6419bebb	226	// calculate the azimuthal angle for the given ADC channel
6419bebb	227
5e86ff99	228	return TMath::ATan2(GetLocalY(det, rob, mcm, ch), GetX(det, rob, mcm));
	229	}
	230
	231	Float_t AliTRDltuParam::GetDist(Int_t det, Int_t rob, Int_t mcm, Int_t ch) const
	232	{
6419bebb	233	// calculate the distance from the origin for the given ADC channel
6419bebb	234
5e86ff99	235	return TMath::Sqrt(GetLocalY(det, rob, mcm, ch)*GetLocalY(det, rob, mcm, ch) +
	236	GetX(det, rob, mcm)*GetX(det, rob, mcm) +
	237	GetLocalZ(det, rob, mcm)*GetLocalZ(det, rob, mcm) );
	238	}