The present commit corresponds to an important change in the way the
[u/mrichter/AliRoot.git] / TRD / AliTRDgtuParam.cxx
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);
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 }