propagate cluster error parametrization
[u/mrichter/AliRoot.git] / TRD / AliTRDgtuParam.cxx
CommitLineData
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
37ClassImp(AliTRDgtuParam)
38
39AliTRDgtuParam *AliTRDgtuParam::fgInstance = 0;
40
41// ----- Bin widths (granularity) -----
42const Float_t AliTRDgtuParam::fgkBinWidthY = 160e-4;
43const Float_t AliTRDgtuParam::fgkBinWidthdY = 140e-4;
44
45// ----- Bit widths (used for internal representation) -----
46const Int_t AliTRDgtuParam::fgkBitWidthY = 13;
47const Int_t AliTRDgtuParam::fgkBitWidthdY = 7;
48const Int_t AliTRDgtuParam::fgkBitWidthYProj = 10;
49const Int_t AliTRDgtuParam::fgkBitExcessY = 4;
50const Int_t AliTRDgtuParam::fgkBitExcessAlpha = 10;
51const Int_t AliTRDgtuParam::fgkBitExcessYProj = 2;
52
53// ----- Tracking parameters -----
54/*
55const Int_t AliTRDgtuParam::fgkNZChannels = 3; // No. of z-channels
56const Int_t AliTRDgtuParam::fgkNLinks = 12; // No. of links
57const Int_t AliTRDgtuParam::fgkFixLayer = 2; // which layer is fixed for the generation of the z-channel map
58const Int_t AliTRDgtuParam::fgkDeltaY = 39; // accepted deviation in y_proj, default: 9
59const Int_t AliTRDgtuParam::fgkDeltaAlpha = 31; // accepted deviation in alpha, default: 11
60const Int_t AliTRDgtuParam::fgkNRefLayers = 3; // no. of reference layers
61*/
62
63AliTRDgtuParam::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
78AliTRDgtuParam::~AliTRDgtuParam()
79{
80 // dtor
81
82 delete fGeo;
83 delete [] fRefLayers;
84}
85
86AliTRDgtuParam* 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
96void AliTRDgtuParam::Terminate()
97{
98 // destruct the instance
99
100 if (fgInstance != 0) {
101 delete fgInstance;
102 fgInstance = 0x0;
103 }
104}
105
106Bool_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
111Int_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
116Int_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
124Int_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
189Bool_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
236Int_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
244Int_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
253Int_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
260Bool_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
294Float_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
303Float_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
313Float_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/*
324Float_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
346Bool_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
372Bool_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
394Float_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);
400 printf("c1: %f\n", c1);
401 Float_t r = (375. / 10000.) * c1 * 256 / (a >> 1);
402 return r;
403
404 Float_t y1 = a + b*x1;
405 Float_t y2 = a + b*x2;
406 Float_t alpha = TMath::Abs( TMath::ATan(y2/x2) - TMath::ATan(y1/x1) );
407 d = TMath::Sqrt( TMath::Power(x2-x1, 2) + TMath::Power(y2-y1, 2) );
408 r = d / 2. / TMath::Sin(alpha);
409 return r;
410}