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