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