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

/**************************************************************************
 * Copyright(c) 1998-1999, 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: AliTRDgtuParam.cxx 28397 2008-09-02 09:33:00Z cblume $ */

////////////////////////////////////////////////////////////////////////////
//                                                                        //
//  Parameters for GTU simulation                                         //
//                                                                        //
//  Author: J. Klein (Jochen.Klein@cern.ch)                               //
//                                                                        //
////////////////////////////////////////////////////////////////////////////

#include "TROOT.h"
#include "TMath.h"
#include "TMatrix.h"
#include "TDecompLU.h"
#include "TGraphAsymmErrors.h"
#include "TCanvas.h"

#include "AliLog.h"
#include "AliTRDgtuParam.h"
#include "AliTRDgeometry.h"
#include "AliTRDpadPlane.h"

ClassImp(AliTRDgtuParam)

AliTRDgtuParam *AliTRDgtuParam::fgInstance = 0;

// ----- matching windows -----
      Int_t     AliTRDgtuParam::fgDeltaY     = 19;
      Int_t     AliTRDgtuParam::fgDeltaAlpha = 21;

// ----- Bin widths (granularity) -----
const Float_t 	AliTRDgtuParam::fgkBinWidthY  = 160e-4;
const Float_t 	AliTRDgtuParam::fgkBinWidthdY = 140e-4;

// ----- Bit widths (used for internal representation) -----
const Int_t 	AliTRDgtuParam::fgkBitWidthY      = 13;
const Int_t 	AliTRDgtuParam::fgkBitWidthdY     = 7;
const Int_t 	AliTRDgtuParam::fgkBitWidthYProj  = 10;
const Int_t 	AliTRDgtuParam::fgkBitExcessY 	  = 4;
const Int_t 	AliTRDgtuParam::fgkBitExcessAlpha = 10;
const Int_t 	AliTRDgtuParam::fgkBitExcessYProj = 2;

// ----- Tracking parameters -----
/*
const Int_t AliTRDgtuParam::fgkNZChannels = 3; // No. of z-channels
const Int_t AliTRDgtuParam::fgkNLinks = 12;	// No. of links
const Int_t AliTRDgtuParam::fgkFixLayer = 2;	// which layer is fixed for the generation of the z-channel map
const Int_t AliTRDgtuParam::fgkDeltaY = 39;	// accepted deviation in y_proj, default: 9
const Int_t AliTRDgtuParam::fgkDeltaAlpha = 31; // accepted deviation in alpha, default: 11
const Int_t AliTRDgtuParam::fgkNRefLayers = 3;	 // no. of reference layers
*/

AliTRDgtuParam::AliTRDgtuParam() :
  fVertexSize(20.0),
  fCurrTrackletMask(0),
  fRefLayers(0x0),
  fMagField(0.5),
  fGeo(0x0)
{
  // default ctor
  fGeo = new AliTRDgeometry();
  fRefLayers = new Int_t[fgkNRefLayers];
  fRefLayers[0] = 3;
  fRefLayers[1] = 2;
  fRefLayers[2] = 1;
  for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
    fAki[iLayer] = 0.;
    fBki[iLayer] = 0.;
    fCki[iLayer] = 0.;
  }

  GenerateZChannelMap();
}

AliTRDgtuParam::~AliTRDgtuParam()
{
  // dtor

  delete fGeo;
  delete [] fRefLayers;
}

AliTRDgtuParam* AliTRDgtuParam::Instance()
{
  // get (or create) the single instance

  if (fgInstance == 0)
    fgInstance = new AliTRDgtuParam();

  return fgInstance;
}

void AliTRDgtuParam::Terminate()
{
  // destruct the instance

  if (fgInstance != 0) {
    delete fgInstance;
    fgInstance = 0x0;
  }
}

Bool_t AliTRDgtuParam::IsInZChannel(Int_t stack, Int_t layer, Int_t zchannel, Int_t zpos) const
{
  return (fZSubChannel[stack][zchannel][layer][zpos] != 0);
}

Int_t AliTRDgtuParam::GetZSubchannel(Int_t stack, Int_t layer, Int_t zchannel, Int_t zpos) const
{
  return fZSubChannel[stack][zchannel][layer][zpos];
}

Int_t AliTRDgtuParam::GetRefLayer(Int_t refLayerIdx) const
{
  // returns the reference layer indexed by refLayerIdx

  if (refLayerIdx >= 0 && refLayerIdx < fgkNRefLayers)
    return fRefLayers[refLayerIdx];
  else
    return -1;
}

Int_t AliTRDgtuParam::GenerateZChannelMap()
{
  // generate the z-channel map
  // assuming that the tracks come from the vertex
  // +/- fVertexSize in z-direction

  Int_t iSec = 0; // sector is irrelevant
  Bool_t collision = kFALSE;

  for (Int_t iStack = 0; iStack < fGeo->Nstack(); iStack++) {

    Float_t x[6] = { 0 };
    Float_t z[6][16] = {{ 0 }};
    Float_t dZ[6][16] = {{ 0 }};

    for (Int_t iLayer = 0; iLayer < fGeo->Nlayer(); iLayer++) {
      AliTRDpadPlane *pp = fGeo->GetPadPlane(iLayer, iStack);
      x[iLayer]  = fGeo->GetTime0(iLayer) - fGeo->CdrHght(); // ???
      for (Int_t iRow = 0; iRow < fGeo->GetRowMax(iLayer, iStack, iSec); iRow++) {
	z[iLayer][iRow]  = pp->GetRowPos(iRow); // this is the right (pos. z-direction) border of the pad
	dZ[iLayer][iRow] = pp->GetRowSize(iRow); // length of the pad in z-direction
	for (Int_t i = 0; i < fgkNZChannels; i++)
	    fZSubChannel[iStack][i][iLayer][iRow] = 0;
      }
    }

    for (Int_t fixRow = 0; fixRow < fGeo->GetRowMax(fgkFixLayer, iStack, iSec); fixRow++) {

      Double_t fixZmin = z[fgkFixLayer][fixRow] - dZ[fgkFixLayer][fixRow];
      Double_t fixZmax = z[fgkFixLayer][fixRow];
      Double_t fixX    = x[fgkFixLayer] + 1.5; // ??? 1.5 from where?

      for (Int_t iLayer = 0; iLayer < fGeo->Nlayer(); iLayer++) {
	Double_t leftZ, rightZ;

	if (iLayer <= fgkFixLayer) {
	  leftZ  = (fixZmin + fVertexSize) * (x[iLayer] + 1.5) / fixX - fVertexSize;
	  rightZ = (fixZmax - fVertexSize) * (x[iLayer] + 1.5) / fixX + fVertexSize;
	}
	else {
	  leftZ  = (fixZmin - fVertexSize) * (x[iLayer] + 1.5) / fixX + fVertexSize;
	  rightZ = (fixZmax + fVertexSize) * (x[iLayer] + 1.5) / fixX - fVertexSize;
	}

	Double_t epsilon = 0.001;
	for (Int_t iRow = 0; iRow < fGeo->GetRowMax(iLayer, iStack, iSec); iRow++) {
	  if ( (z[iLayer][iRow] )                    > (leftZ  + epsilon) &&
	       (z[iLayer][iRow] - dZ[iLayer][iRow] ) < (rightZ - epsilon) ) {
	    fZChannelMap[iStack][fixRow][iLayer][iRow] = 1;
	    if (fZSubChannel[iStack][fixRow % fgkNZChannels][iLayer][iRow] != 0) {
	      AliError("Collision in Z-Channel assignment occured! No reliable tracking!!!");
	      collision = kTRUE;
	    }
	    else
	      fZSubChannel[iStack][fixRow % fgkNZChannels][iLayer][iRow] = fixRow / fgkNZChannels + 1;
	  }

	}
      }
    }
  }

  return ~collision;
}

Bool_t AliTRDgtuParam::DisplayZChannelMap(Int_t zchannel, Int_t subchannel) const
{
  // display the z-channel map

  if (zchannel >= fgkNZChannels) {
    AliError("Invalid Z channel!");
    return kFALSE;
  }

  Int_t zchmin = zchannel >= 0 ? zchannel : 0;
  Int_t zchmax = zchannel >= 0 ? zchannel + 1 : fgkNZChannels;
  Int_t i = 0;
  Int_t j = 0;
  TCanvas *c = new TCanvas("zchmap", "Z-Chhannel Mapping");
  c->cd();
  TGraph **graphz = new TGraph*[fgkNZChannels];
  for (Int_t zch = zchmin; zch < zchmax; zch++)
    graphz[zch] = new TGraph;
  TGraphAsymmErrors *graph = new TGraphAsymmErrors();
  graph->SetTitle("Z-Channel Map");
  graph->SetPoint(i, 0, 0); // vertex
  graph->SetPointError(i++, 20, 20, 0, 0);
  //  graph->SetRange //????
  for (Int_t iLayer = 0; iLayer < fGeo->Nlayer(); iLayer++) {
    for (Int_t iStack = 0; iStack < fGeo->Nstack(); iStack++) {
      AliTRDpadPlane *pp = fGeo->GetPadPlane(iLayer, iStack);
      for (Int_t iRow = 0; iRow < fGeo->GetRowMax(iLayer, iStack, 0); iRow++) {
	graph->SetPoint(i, pp->GetRowPos(iRow), fGeo->GetTime0(iLayer) - fGeo->CdrHght());
	graph->SetPointError(i++, pp->GetRowSize(iRow), 0, 0, 0);
	for (Int_t zch = zchmin; zch < zchmax; zch++)
	  if (fZSubChannel[iStack][zch][iLayer][iRow] != 0)
	    if (subchannel == 0 || fZSubChannel[iStack][zch][iLayer][iRow] == subchannel)
	      graphz[zch]->SetPoint(j++, pp->GetRowPos(iRow)  - pp->GetRowSize(iRow)/2, fGeo->GetTime0(iLayer) - fGeo->CdrHght());
      }
    }
  }
  graph->SetMarkerStyle(kDot);
  graph->Draw("AP");
  gROOT->Add(graph);
  for (Int_t zch = zchmin; zch < zchmax; zch++) {
    graphz[zch]->SetMarkerStyle(kCircle);
    graphz[zch]->SetMarkerColor(zch+2);
    graphz[zch]->SetMarkerSize(0.3 + zch*0.2);
    graphz[zch]->Draw("P");
    gROOT->Add(graphz[zch]);
  }
  delete [] graphz;
  return kTRUE;
}

Int_t AliTRDgtuParam::GetCiAlpha(Int_t layer) const
{
  // get the constant for the calculation of alpha

  Int_t ci = (Int_t) (GetChamberThickness() / fGeo->GetTime0(layer) * GetBinWidthY() / GetBinWidthdY() * (1 << (GetBitExcessAlpha() + GetBitExcessY() + 1)) );
  return ci;
}

Int_t AliTRDgtuParam::GetCiYProj(Int_t layer) const
{
  // get the constant for the calculation of y_proj

  Float_t xmid = (fGeo->GetTime0(0) + fGeo->GetTime0(fGeo->Nlayer()-1)) / 2.;
  Int_t ci = (Int_t) (- (fGeo->GetTime0(layer) - xmid) / GetChamberThickness() * GetBinWidthdY() / GetBinWidthY() * (1 << GetBitExcessYProj()) );
  return ci;
}

Int_t AliTRDgtuParam::GetYt(Int_t stack, Int_t layer, Int_t zrow) const
{
    return (Int_t) (- ( (layer % 2 ? 1 : -1) *
			(GetGeo()->GetPadPlane(layer, stack)->GetRowPos(zrow) - GetGeo()->GetPadPlane(layer, stack)->GetRowSize(zrow) / 2) *
			TMath::Tan(- 2.0 / 180.0 * TMath::Pi()) ) / 0.016 );
}

Bool_t AliTRDgtuParam::GenerateRecoCoefficients(Int_t trackletMask)
{
  // calculate the coefficients for the straight line fit
  // depending on the mask of contributing tracklets

  fCurrTrackletMask = trackletMask;

  TMatrix a(GetNLayers(), 3);
  TMatrix b(3, GetNLayers());
  TMatrix c(3, 3);

  for (Int_t layer = 0; layer < GetNLayers(); layer++) {
      if ( (trackletMask & (1 << layer)) == 0) {
	  a(layer, 0) = 0;
	  a(layer, 1) = 0;
	  a(layer, 2) = 0;
      }
      else {
	  a(layer, 0) = 1;
	  a(layer, 1) = fGeo->GetTime0(layer);
	  a(layer, 2) = (layer % 2 ? 1 : -1) * fGeo->GetTime0(layer);
      }
  }

  b.Transpose(a);
  c = b * a;
  c.InvertFast();
  b = c * b;

  for (Int_t layer = 0; layer < GetNLayers(); layer++) {
      fAki[layer] = b.GetMatrixArray()[layer];
      fBki[layer] = b.GetMatrixArray()[GetNLayers() + layer];
      fCki[layer] = b.GetMatrixArray()[2 * GetNLayers() + layer];
    }
  return kTRUE;
}

Float_t AliTRDgtuParam::GetAki(Int_t k, Int_t i)
{
  // get A_ki for the calculation of the tracking parameters
  if (fCurrTrackletMask != k)
    GenerateRecoCoefficients(k);

  return fAki[i];
}

Float_t AliTRDgtuParam::GetBki(Int_t k, Int_t i)
{
  // get B_ki for the calculation of the tracking parameters

  if (fCurrTrackletMask != k)
    GenerateRecoCoefficients(k);

  return fBki[i];
}

Float_t AliTRDgtuParam::GetCki(Int_t k, Int_t i)
{
  // get B_ki for the calculation of the tracking parameters

  if (fCurrTrackletMask != k)
    GenerateRecoCoefficients(k);

  return fCki[i];
}

/*
Float_t AliTRDgtuParam::GetD(Int_t k) const
{
  // get the determinant for the calculation of the tracking parameters

  TMatrix t(3, 3);
  for (Int_t i = 0; i < GetNLayers(); i++) {
    if ( !((k >> i) & 0x1) )
      continue;
    Float_t xi = fGeo->GetTime0(i);
    t(0,0) += 1;
    t(1,0) += xi;
    t(2,0) += TMath::Power(-1, i) * xi;
    t(0,1) += xi;
    t(1,1) += TMath::Power(xi, 2);
    t(2,1) += TMath::Power(-1, i) * TMath::Power(xi, 2);
    t(0,2) += TMath::Power(-1, i) * xi;
    t(1,2) += TMath::Power(-1, i) * TMath::Power(xi, 2);
    t(2,2) += TMath::Power(xi, 2);
  }
  return t.Determinant();
}

Bool_t AliTRDgtuParam::GetFitParams(TVectorD& rhs, Int_t k)
{
  // calculate the fitting parameters
  // will be changed!

  TMatrix t(3,3);
  for (Int_t i = 0; i < GetNLayers(); i++) {
    if ( !((k >> i) & 0x1) )
      continue;
    Float_t xi = fGeo->GetTime0(i);
    t(0,0) += 1;
    t(1,0) += xi;
    t(2,0) += TMath::Power(-1, i) * xi;
    t(0,1) += xi;
    t(1,1) += TMath::Power(xi, 2);
    t(2,1) += TMath::Power(-1, i) * TMath::Power(xi, 2);
    t(0,2) -= TMath::Power(-1, i) * xi;
    t(1,2) -= TMath::Power(-1, i) * TMath::Power(xi, 2);
    t(2,2) -= TMath::Power(xi, 2);
  }
  TDecompLU lr(t);
  lr.Solve(rhs);
  return lr.Decompose();
}
*/

Bool_t AliTRDgtuParam::GetIntersectionPoints(Int_t k, Float_t &x1, Float_t &x2)
{
  // get the x-coord. of the assumed circle/straight line intersection points

  Int_t l1 = -1;
  Int_t l2 = -1;
  Int_t nHits = 0;
  for (Int_t layer = 0; layer < GetNLayers(); layer++) {
    if ( (k >> layer) & 0x1 ) {
      if (l1 < 0)
	l1 = layer;
      l2 = layer;
      nHits++;
    }
  }

  if ( (l1 >= 0) && (l2 >= 0) ) {
    x1 = fGeo->GetTime0(l1) + 10./6 * (nHits -1);
    x2 = fGeo->GetTime0(l2) - 10./6 * (nHits -1);
    return kTRUE;
  }
  else
    return kFALSE;
}

Float_t AliTRDgtuParam::GetPt(Int_t a, Float_t /* b */, Float_t x1, Float_t x2) const
{
  // returns 0.3 * B * 1/a (1/128 GeV/c)
  // a : offset, b : slope (not used)

  Float_t c1 = x1 * x2 / 2. / 10000.; // conversion cm to m
  Float_t r = 0;
  if ( (a >> 1) != 0)
    r = (0.3 * fMagField / 2. / (fgkBinWidthY/100.)) * (((Int_t) c1) << 8) / (a >> 1); //??? why shift of a?
  return r;
}
Commit	Line	Data
52c19022	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id: AliTRDgtuParam.cxx 28397 2008-09-02 09:33:00Z cblume $ */
	17
	18	////////////////////////////////////////////////////////////////////////////
	19	// //
	20	// Parameters for GTU simulation //
	21	// //
	22	// Author: J. Klein (Jochen.Klein@cern.ch) //
	23	// //
	24	////////////////////////////////////////////////////////////////////////////
	25
637666cd	26	#include "TROOT.h"
52c19022	27	#include "TMath.h"
	28	#include "TMatrix.h"
	29	#include "TDecompLU.h"
	30	#include "TGraphAsymmErrors.h"
	31	#include "TCanvas.h"
	32
	33	#include "AliLog.h"
	34	#include "AliTRDgtuParam.h"
	35	#include "AliTRDgeometry.h"
	36	#include "AliTRDpadPlane.h"
	37
	38	ClassImp(AliTRDgtuParam)
	39
	40	AliTRDgtuParam *AliTRDgtuParam::fgInstance = 0;
	41
2cf67435	42	// ----- matching windows -----
	43	Int_t AliTRDgtuParam::fgDeltaY = 19;
	44	Int_t AliTRDgtuParam::fgDeltaAlpha = 21;
	45
52c19022	46	// ----- Bin widths (granularity) -----
	47	const Float_t AliTRDgtuParam::fgkBinWidthY = 160e-4;
	48	const Float_t AliTRDgtuParam::fgkBinWidthdY = 140e-4;
	49
	50	// ----- Bit widths (used for internal representation) -----
	51	const Int_t AliTRDgtuParam::fgkBitWidthY = 13;
5f006bd7	52	const Int_t AliTRDgtuParam::fgkBitWidthdY = 7;
52c19022	53	const Int_t AliTRDgtuParam::fgkBitWidthYProj = 10;
5f006bd7	54	const Int_t AliTRDgtuParam::fgkBitExcessY = 4;
	55	const Int_t AliTRDgtuParam::fgkBitExcessAlpha = 10;
	56	const Int_t AliTRDgtuParam::fgkBitExcessYProj = 2;
52c19022	57
	58	// ----- Tracking parameters -----
	59	/*
	60	const Int_t AliTRDgtuParam::fgkNZChannels = 3; // No. of z-channels
	61	const Int_t AliTRDgtuParam::fgkNLinks = 12; // No. of links
	62	const Int_t AliTRDgtuParam::fgkFixLayer = 2; // which layer is fixed for the generation of the z-channel map
	63	const Int_t AliTRDgtuParam::fgkDeltaY = 39; // accepted deviation in y_proj, default: 9
	64	const Int_t AliTRDgtuParam::fgkDeltaAlpha = 31; // accepted deviation in alpha, default: 11
	65	const Int_t AliTRDgtuParam::fgkNRefLayers = 3; // no. of reference layers
	66	*/
	67
	68	AliTRDgtuParam::AliTRDgtuParam() :
	69	fVertexSize(20.0),
	70	fCurrTrackletMask(0),
	71	fRefLayers(0x0),
b491d23b	72	fMagField(0.5),
52c19022	73	fGeo(0x0)
	74	{
	75	// default ctor
	76	fGeo = new AliTRDgeometry();
	77	fRefLayers = new Int_t[fgkNRefLayers];
	78	fRefLayers[0] = 3;
	79	fRefLayers[1] = 2;
	80	fRefLayers[2] = 1;
637666cd	81	for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
	82	fAki[iLayer] = 0.;
	83	fBki[iLayer] = 0.;
	84	fCki[iLayer] = 0.;
	85	}
	86
5f006bd7	87	GenerateZChannelMap();
52c19022	88	}
52c19022	89
5f006bd7	90	AliTRDgtuParam::~AliTRDgtuParam()
52c19022	91	{
	92	// dtor
	93
	94	delete fGeo;
	95	delete [] fRefLayers;
	96	}
	97
5f006bd7	98	AliTRDgtuParam* AliTRDgtuParam::Instance()
52c19022	99	{
	100	// get (or create) the single instance
	101
5f006bd7	102	if (fgInstance == 0)
52c19022	103	fgInstance = new AliTRDgtuParam();
	104
	105	return fgInstance;
	106	}
	107
5f006bd7	108	void AliTRDgtuParam::Terminate()
52c19022	109	{
	110	// destruct the instance
	111
	112	if (fgInstance != 0) {
	113	delete fgInstance;
	114	fgInstance = 0x0;
	115	}
	116	}
	117
5f006bd7	118	Bool_t AliTRDgtuParam::IsInZChannel(Int_t stack, Int_t layer, Int_t zchannel, Int_t zpos) const
52c19022	119	{
	120	return (fZSubChannel[stack][zchannel][layer][zpos] != 0);
	121	}
	122
	123	Int_t AliTRDgtuParam::GetZSubchannel(Int_t stack, Int_t layer, Int_t zchannel, Int_t zpos) const
	124	{
	125	return fZSubChannel[stack][zchannel][layer][zpos];
	126	}
	127
5f006bd7	128	Int_t AliTRDgtuParam::GetRefLayer(Int_t refLayerIdx) const
52c19022	129	{
36dc3337	130	// returns the reference layer indexed by refLayerIdx
36dc3337	131
52c19022	132	if (refLayerIdx >= 0 && refLayerIdx < fgkNRefLayers)
52c19022	133	return fRefLayers[refLayerIdx];
5f006bd7	134	else
52c19022	135	return -1;
	136	}
	137
5f006bd7	138	Int_t AliTRDgtuParam::GenerateZChannelMap()
52c19022	139	{
52c19022	140	// generate the z-channel map
5f006bd7	141	// assuming that the tracks come from the vertex
52c19022	142	// +/- fVertexSize in z-direction
	143
	144	Int_t iSec = 0; // sector is irrelevant
	145	Bool_t collision = kFALSE;
	146
	147	for (Int_t iStack = 0; iStack < fGeo->Nstack(); iStack++) {
	148
36dc3337	149	Float_t x[6] = { 0 };
36dc3337	150	Float_t z[6][16] = {{ 0 }};
52c19022	151	Float_t dZ[6][16] = {{ 0 }};
5f006bd7	152
52c19022	153	for (Int_t iLayer = 0; iLayer < fGeo->Nlayer(); iLayer++) {
52c19022	154	AliTRDpadPlane *pp = fGeo->GetPadPlane(iLayer, iStack);
36dc3337	155	x[iLayer] = fGeo->GetTime0(iLayer) - fGeo->CdrHght(); // ???
52c19022	156	for (Int_t iRow = 0; iRow < fGeo->GetRowMax(iLayer, iStack, iSec); iRow++) {
36dc3337	157	z[iLayer][iRow] = pp->GetRowPos(iRow); // this is the right (pos. z-direction) border of the pad
52c19022	158	dZ[iLayer][iRow] = pp->GetRowSize(iRow); // length of the pad in z-direction
5f006bd7	159	for (Int_t i = 0; i < fgkNZChannels; i++)
52c19022	160	fZSubChannel[iStack][i][iLayer][iRow] = 0;
	161	}
	162	}
	163
	164	for (Int_t fixRow = 0; fixRow < fGeo->GetRowMax(fgkFixLayer, iStack, iSec); fixRow++) {
5f006bd7	165
5f006bd7	166	Double_t fixZmin = z[fgkFixLayer][fixRow] - dZ[fgkFixLayer][fixRow];
36dc3337	167	Double_t fixZmax = z[fgkFixLayer][fixRow];
5f006bd7	168	Double_t fixX = x[fgkFixLayer] + 1.5; // ??? 1.5 from where?
52c19022	169
	170	for (Int_t iLayer = 0; iLayer < fGeo->Nlayer(); iLayer++) {
	171	Double_t leftZ, rightZ;
5f006bd7	172
52c19022	173	if (iLayer <= fgkFixLayer) {
36dc3337	174	leftZ = (fixZmin + fVertexSize) * (x[iLayer] + 1.5) / fixX - fVertexSize;
36dc3337	175	rightZ = (fixZmax - fVertexSize) * (x[iLayer] + 1.5) / fixX + fVertexSize;
52c19022	176	}
52c19022	177	else {
36dc3337	178	leftZ = (fixZmin - fVertexSize) * (x[iLayer] + 1.5) / fixX + fVertexSize;
36dc3337	179	rightZ = (fixZmax + fVertexSize) * (x[iLayer] + 1.5) / fixX - fVertexSize;
52c19022	180	}
5f006bd7	181
52c19022	182	Double_t epsilon = 0.001;
52c19022	183	for (Int_t iRow = 0; iRow < fGeo->GetRowMax(iLayer, iStack, iSec); iRow++) {
5f006bd7	184	if ( (z[iLayer][iRow] ) > (leftZ + epsilon) &&
36dc3337	185	(z[iLayer][iRow] - dZ[iLayer][iRow] ) < (rightZ - epsilon) ) {
52c19022	186	fZChannelMap[iStack][fixRow][iLayer][iRow] = 1;
	187	if (fZSubChannel[iStack][fixRow % fgkNZChannels][iLayer][iRow] != 0) {
	188	AliError("Collision in Z-Channel assignment occured! No reliable tracking!!!");
	189	collision = kTRUE;
	190	}
5f006bd7	191	else
52c19022	192	fZSubChannel[iStack][fixRow % fgkNZChannels][iLayer][iRow] = fixRow / fgkNZChannels + 1;
	193	}
	194
	195	}
	196	}
	197	}
	198	}
	199
	200	return ~collision;
	201	}
	202
5f006bd7	203	Bool_t AliTRDgtuParam::DisplayZChannelMap(Int_t zchannel, Int_t subchannel) const
52c19022	204	{
5f006bd7	205	// display the z-channel map
52c19022	206
637666cd	207	if (zchannel >= fgkNZChannels) {
52c19022	208	AliError("Invalid Z channel!");
	209	return kFALSE;
	210	}
	211
	212	Int_t zchmin = zchannel >= 0 ? zchannel : 0;
	213	Int_t zchmax = zchannel >= 0 ? zchannel + 1 : fgkNZChannels;
	214	Int_t i = 0;
	215	Int_t j = 0;
	216	TCanvas *c = new TCanvas("zchmap", "Z-Chhannel Mapping");
	217	c->cd();
	218	TGraph *graphz = new TGraph[fgkNZChannels];
5f006bd7	219	for (Int_t zch = zchmin; zch < zchmax; zch++)
52c19022	220	graphz[zch] = new TGraph;
	221	TGraphAsymmErrors *graph = new TGraphAsymmErrors();
	222	graph->SetTitle("Z-Channel Map");
	223	graph->SetPoint(i, 0, 0); // vertex
	224	graph->SetPointError(i++, 20, 20, 0, 0);
	225	// graph->SetRange //????
	226	for (Int_t iLayer = 0; iLayer < fGeo->Nlayer(); iLayer++) {
	227	for (Int_t iStack = 0; iStack < fGeo->Nstack(); iStack++) {
	228	AliTRDpadPlane *pp = fGeo->GetPadPlane(iLayer, iStack);
	229	for (Int_t iRow = 0; iRow < fGeo->GetRowMax(iLayer, iStack, 0); iRow++) {
	230	graph->SetPoint(i, pp->GetRowPos(iRow), fGeo->GetTime0(iLayer) - fGeo->CdrHght());
	231	graph->SetPointError(i++, pp->GetRowSize(iRow), 0, 0, 0);
	232	for (Int_t zch = zchmin; zch < zchmax; zch++)
	233	if (fZSubChannel[iStack][zch][iLayer][iRow] != 0)
	234	if (subchannel == 0 \|\| fZSubChannel[iStack][zch][iLayer][iRow] == subchannel)
	235	graphz[zch]->SetPoint(j++, pp->GetRowPos(iRow) - pp->GetRowSize(iRow)/2, fGeo->GetTime0(iLayer) - fGeo->CdrHght());
	236	}
	237	}
	238	}
	239	graph->SetMarkerStyle(kDot);
	240	graph->Draw("AP");
637666cd	241	gROOT->Add(graph);
52c19022	242	for (Int_t zch = zchmin; zch < zchmax; zch++) {
	243	graphz[zch]->SetMarkerStyle(kCircle);
	244	graphz[zch]->SetMarkerColor(zch+2);
	245	graphz[zch]->SetMarkerSize(0.3 + zch*0.2);
	246	graphz[zch]->Draw("P");
637666cd	247	gROOT->Add(graphz[zch]);
52c19022	248	}
54d34aac	249	delete [] graphz;
52c19022	250	return kTRUE;
	251	}
	252
5f006bd7	253	Int_t AliTRDgtuParam::GetCiAlpha(Int_t layer) const
52c19022	254	{
	255	// get the constant for the calculation of alpha
	256
36dc3337	257	Int_t ci = (Int_t) (GetChamberThickness() / fGeo->GetTime0(layer) * GetBinWidthY() / GetBinWidthdY() * (1 << (GetBitExcessAlpha() + GetBitExcessY() + 1)) );
36dc3337	258	return ci;
52c19022	259	}
52c19022	260
5f006bd7	261	Int_t AliTRDgtuParam::GetCiYProj(Int_t layer) const
52c19022	262	{
	263	// get the constant for the calculation of y_proj
	264
5f006bd7	265	Float_t xmid = (fGeo->GetTime0(0) + fGeo->GetTime0(fGeo->Nlayer()-1)) / 2.;
36dc3337	266	Int_t ci = (Int_t) (- (fGeo->GetTime0(layer) - xmid) / GetChamberThickness() * GetBinWidthdY() / GetBinWidthY() * (1 << GetBitExcessYProj()) );
36dc3337	267	return ci;
52c19022	268	}
	269
	270	Int_t AliTRDgtuParam::GetYt(Int_t stack, Int_t layer, Int_t zrow) const
	271	{
5f006bd7	272	return (Int_t) (- ( (layer % 2 ? 1 : -1) *
5f006bd7	273	(GetGeo()->GetPadPlane(layer, stack)->GetRowPos(zrow) - GetGeo()->GetPadPlane(layer, stack)->GetRowSize(zrow) / 2) *
52c19022	274	TMath::Tan(- 2.0 / 180.0 * TMath::Pi()) ) / 0.016 );
	275	}
	276
5f006bd7	277	Bool_t AliTRDgtuParam::GenerateRecoCoefficients(Int_t trackletMask)
52c19022	278	{
5f006bd7	279	// calculate the coefficients for the straight line fit
36dc3337	280	// depending on the mask of contributing tracklets
36dc3337	281
52c19022	282	fCurrTrackletMask = trackletMask;
	283
	284	TMatrix a(GetNLayers(), 3);
	285	TMatrix b(3, GetNLayers());
	286	TMatrix c(3, 3);
	287
	288	for (Int_t layer = 0; layer < GetNLayers(); layer++) {
	289	if ( (trackletMask & (1 << layer)) == 0) {
	290	a(layer, 0) = 0;
	291	a(layer, 1) = 0;
	292	a(layer, 2) = 0;
5f006bd7	293	}
52c19022	294	else {
	295	a(layer, 0) = 1;
	296	a(layer, 1) = fGeo->GetTime0(layer);
	297	a(layer, 2) = (layer % 2 ? 1 : -1) * fGeo->GetTime0(layer);
	298	}
	299	}
	300
	301	b.Transpose(a);
	302	c = b * a;
	303	c.InvertFast();
	304	b = c * b;
	305
	306	for (Int_t layer = 0; layer < GetNLayers(); layer++) {
	307	fAki[layer] = b.GetMatrixArray()[layer];
	308	fBki[layer] = b.GetMatrixArray()[GetNLayers() + layer];
	309	fCki[layer] = b.GetMatrixArray()[2 * GetNLayers() + layer];
	310	}
	311	return kTRUE;
	312	}
	313
5f006bd7	314	Float_t AliTRDgtuParam::GetAki(Int_t k, Int_t i)
52c19022	315	{
	316	// get A_ki for the calculation of the tracking parameters
	317	if (fCurrTrackletMask != k)
	318	GenerateRecoCoefficients(k);
	319
	320	return fAki[i];
	321	}
	322
5f006bd7	323	Float_t AliTRDgtuParam::GetBki(Int_t k, Int_t i)
52c19022	324	{
	325	// get B_ki for the calculation of the tracking parameters
	326
	327	if (fCurrTrackletMask != k)
	328	GenerateRecoCoefficients(k);
	329
	330	return fBki[i];
	331	}
	332
5f006bd7	333	Float_t AliTRDgtuParam::GetCki(Int_t k, Int_t i)
52c19022	334	{
	335	// get B_ki for the calculation of the tracking parameters
	336
	337	if (fCurrTrackletMask != k)
	338	GenerateRecoCoefficients(k);
	339
	340	return fCki[i];
	341	}
	342
	343	/*
5f006bd7	344	Float_t AliTRDgtuParam::GetD(Int_t k) const
52c19022	345	{
	346	// get the determinant for the calculation of the tracking parameters
	347
	348	TMatrix t(3, 3);
	349	for (Int_t i = 0; i < GetNLayers(); i++) {
	350	if ( !((k >> i) & 0x1) )
	351	continue;
	352	Float_t xi = fGeo->GetTime0(i);
	353	t(0,0) += 1;
	354	t(1,0) += xi;
	355	t(2,0) += TMath::Power(-1, i) * xi;
	356	t(0,1) += xi;
	357	t(1,1) += TMath::Power(xi, 2);
	358	t(2,1) += TMath::Power(-1, i) * TMath::Power(xi, 2);
	359	t(0,2) += TMath::Power(-1, i) * xi;
	360	t(1,2) += TMath::Power(-1, i) * TMath::Power(xi, 2);
	361	t(2,2) += TMath::Power(xi, 2);
	362	}
	363	return t.Determinant();
	364	}
	365
5f006bd7	366	Bool_t AliTRDgtuParam::GetFitParams(TVectorD& rhs, Int_t k)
52c19022	367	{
	368	// calculate the fitting parameters
	369	// will be changed!
	370
	371	TMatrix t(3,3);
	372	for (Int_t i = 0; i < GetNLayers(); i++) {
	373	if ( !((k >> i) & 0x1) )
	374	continue;
	375	Float_t xi = fGeo->GetTime0(i);
	376	t(0,0) += 1;
	377	t(1,0) += xi;
	378	t(2,0) += TMath::Power(-1, i) * xi;
	379	t(0,1) += xi;
	380	t(1,1) += TMath::Power(xi, 2);
	381	t(2,1) += TMath::Power(-1, i) * TMath::Power(xi, 2);
	382	t(0,2) -= TMath::Power(-1, i) * xi;
	383	t(1,2) -= TMath::Power(-1, i) * TMath::Power(xi, 2);
	384	t(2,2) -= TMath::Power(xi, 2);
	385	}
	386	TDecompLU lr(t);
	387	lr.Solve(rhs);
	388	return lr.Decompose();
	389	}
	390	*/
	391
5f006bd7	392	Bool_t AliTRDgtuParam::GetIntersectionPoints(Int_t k, Float_t &x1, Float_t &x2)
52c19022	393	{
	394	// get the x-coord. of the assumed circle/straight line intersection points
	395
	396	Int_t l1 = -1;
	397	Int_t l2 = -1;
	398	Int_t nHits = 0;
	399	for (Int_t layer = 0; layer < GetNLayers(); layer++) {
	400	if ( (k >> layer) & 0x1 ) {
5f006bd7	401	if (l1 < 0)
52c19022	402	l1 = layer;
	403	l2 = layer;
	404	nHits++;
	405	}
	406	}
	407
637666cd	408	if ( (l1 >= 0) && (l2 >= 0) ) {
	409	x1 = fGeo->GetTime0(l1) + 10./6 * (nHits -1);
	410	x2 = fGeo->GetTime0(l2) - 10./6 * (nHits -1);
	411	return kTRUE;
	412	}
	413	else
	414	return kFALSE;
52c19022	415	}
52c19022	416
b491d23b	417	Float_t AliTRDgtuParam::GetPt(Int_t a, Float_t /* b */, Float_t x1, Float_t x2) const
52c19022	418	{
b491d23b	419	// returns 0.3 * B * 1/a (1/128 GeV/c)
	420	// a : offset, b : slope (not used)
	421
	422	Float_t c1 = x1 * x2 / 2. / 10000.; // conversion cm to m
80f93426	423	Float_t r = 0;
80f93426	424	if ( (a >> 1) != 0)
b491d23b	425	r = (0.3 * fMagField / 2. / (fgkBinWidthY/100.)) * (((Int_t) c1) << 8) / (a >> 1); //??? why shift of a?
52c19022	426	return r;
52c19022	427	}