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