[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 "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;

// ----- 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),
  fGeo(0x0)
{
  // default ctor
  fGeo = new AliTRDgeometry();
  fRefLayers = new Int_t[fgkNRefLayers];
  fRefLayers[0] = 3;
  fRefLayers[1] = 2;
  fRefLayers[2] = 1;
  zChannelGen(); 
}

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 
{
  if (refLayerIdx >= 0 && refLayerIdx < fgkNRefLayers)
    return fRefLayers[refLayerIdx];
  else 
    return -1;
}

Int_t AliTRDgtuParam::zChannelGen() 
{
  // 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");
  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");
  }
  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) 
{
  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++;
    }
  }

  x1 = fGeo->GetTime0(l1) + 10./6 * (nHits -1);
  x2 = fGeo->GetTime0(l2) - 10./6 * (nHits -1);

  return ( (l1 >= 0) && (l2 >= 0) );
}

Float_t AliTRDgtuParam::GetRadius(Int_t a, Float_t b, Float_t x1, Float_t x2) 
{
  // get the radius for the track
  Float_t d = (1 + b * b /2 ) * (x2 - x1);
  Float_t c1 = x1 * x2 / 2;
//  Float_t c2 = (x1 + x2) / (x1 * x2);
//  printf("c1: %f\n", c1);
  Float_t r = 0;
  if ( (a >> 1) != 0)
    r = (375. / 10000.) * c1 * 256 / (a >> 1);
  return r;

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