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