TRD/AliTRDtrackingChamber.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: AliTRDtrackingChamber.cxx 23810 2008-02-08 09:00:27Z hristov $ */
  17
  18 ////////////////////////////////////////////////////////////////////////////
  19 //                                                                        //
  20 //  Tracking in one chamber                                               //
  21 //                                                                        //
  22 //  Authors:                                                              //
  23 //    Alex Bercuci <A.Bercuci@gsi.de>                                     //
  24 //    Markus Fasel <M.Fasel@gsi.de>                                       //
  25 //                                                                        //
  26 ////////////////////////////////////////////////////////////////////////////
  27
  28 #include "AliTRDtrackingChamber.h"
  29
  30 #include "TMath.h"
  31 #include "TMatrixTBase.h"
  32 #include <TTreeStream.h>
  33
  34 #include "AliTRDReconstructor.h"
  35 #include "AliTRDrecoParam.h"
  36 #include "AliTRDtrackerV1.h"
  37 #include "AliTRDgeometry.h"
  38 #include "AliTRDpadPlane.h"
  39 #include "AliTRDcalibDB.h"
  40 #include "AliTRDCommonParam.h"
  41 #include "Cal/AliTRDCalDet.h"
  42 #include "Cal/AliTRDCalROC.h"
  43
  44 ClassImp(AliTRDtrackingChamber)
  45
  46 //_______________________________________________________
  47 AliTRDtrackingChamber::AliTRDtrackingChamber()
  48   :TObject()
  49   ,fDetector(-1)
  50   ,fX0(0.)
  51   ,fExB(0.)
  52   ,fVD(0.)
  53   ,fT0(0.)
  54   ,fS2PRF(0.)
  55   ,fDiffL(0.)
  56   ,fDiffT(0.)
  57 {}
  58
  59 //_______________________________________________________
  60 void AliTRDtrackingChamber::Clear(const Option_t *opt)
  61 {
  62   for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++) fTB[itb].Clear(opt);
  63 }
  64
  65 //_______________________________________________________
  66 void AliTRDtrackingChamber::InsertCluster(AliTRDcluster *c, Int_t index, Bool_t hlt)
  67 {
  68 // Add cluster to TB container and recalculate error parameterization (for HLT)
  69
  70   fTB[c->GetPadTime()].InsertCluster(c, index);
  71   if(!hlt) return;
  72
  73   // Define approximate error parameterization for HLT clusters
  74   // if needed the fix values of
  75   //  - drift length of 1.5 can be replaced with c->GetXloc()
  76   //  - pad length can be cached from geometry in Init()
  77   c->SetSigmaY2(fS2PRF, fDiffT, fExB, 1.5);
  78   c->SetSigmaZ2(6.75);
  79 }
  80
  81 //_______________________________________________________
  82 Bool_t AliTRDtrackingChamber::Build(AliTRDgeometry *const geo)
  83 {
  84 // Init chamber and all time bins (AliTRDchamberTimeBin)
  85 // Calculates radial position of the chamber based on
  86 // radial positions of the time bins (calibration/alignment aware)
  87 //
  88   if(fDetector < 0 || fDetector >= AliTRDgeometry::kNdet){
  89     AliWarning(Form("Detector index not set correctly to %d", fDetector));
  90     return kFALSE;
  91   }
  92
  93   Int_t stack = AliTRDgeometry::GetStack(fDetector);
  94   Int_t layer = AliTRDgeometry::GetLayer(fDetector);
  95   AliTRDpadPlane *pp = geo->GetPadPlane(layer, stack);
  96   Double_t zl = pp->GetRow0ROC() - pp->GetRowEndROC();
  97   Double_t z0 = geo->GetRow0(layer, stack, 0) - zl;
  98   Int_t nrows = pp->GetNrows();
  99
 100   Int_t index[50], jtb = 0;
 101   for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++){
 102     if(!fTB[itb]) continue;
 103     fTB[itb].SetRange(z0, zl);
 104     fTB[itb].SetNRows(nrows);
 105     fTB[itb].BuildIndices();
 106     index[jtb++] = itb;
 107   }
 108   if(jtb<2) return kFALSE;
 109
 110   // ESTIMATE POSITION OF PAD PLANE FOR THIS CHAMBER
 111   Int_t t0 = Int_t(fT0);
 112   fTB[t0].SetT0();
 113   Double_t x0 = fTB[index[0]].GetX();
 114   Double_t x1 = fTB[index[1]].GetX();
 115   Double_t dx = (x0 - x1)/(index[1] - index[0]);
 116   fX0 = x0 + dx*(index[0] - t0);
 117   return kTRUE;
 118 }
 119
 120
 121 //_______________________________________________________
 122 void AliTRDtrackingChamber::Init(Int_t det)
 123 {
 124 // Init detector number and Cache calibration parameters
 125
 126   fDetector = det;
 127   AliTRDcalibDB *calib = AliTRDcalibDB::Instance();
 128   fT0    = calib->GetT0Average(fDetector);
 129   fVD    = calib->GetVdriftAverage(fDetector);
 130   fS2PRF = calib->GetPRFROC(fDetector)->GetMean(); fS2PRF *= fS2PRF;
 131   fExB   = AliTRDCommonParam::Instance()->GetOmegaTau(fVD);
 132   AliTRDCommonParam::Instance()->GetDiffCoeff(fDiffL,
 133   fDiffT, fVD);
 134 }
 135
 136 //_______________________________________________________
 137 Int_t AliTRDtrackingChamber::GetNClusters() const
 138 {
 139 // Basic loop method
 140 // Returns number of clusters in chamber
 141 //
 142   Int_t n = 0;
 143   for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++){
 144     n += Int_t(fTB[itb]);
 145   }
 146   return n;
 147 }
 148
 149 //_______________________________________________________
 150 void AliTRDtrackingChamber::Bootstrap(const AliTRDReconstructor *rec)
 151 {
 152 // Basic loop method
 153 // Bootstrap each time bin
 154 //
 155   AliTRDchamberTimeBin *jtb = &fTB[0];
 156   for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++, ++jtb){
 157     (*jtb).Bootstrap(rec, fDetector);
 158   }
 159 }
 160
 161 //_______________________________________________________
 162 void  AliTRDtrackingChamber::SetOwner()
 163 {
 164 // Basic loop method
 165 // Set ownership in time bins
 166 //
 167   AliTRDchamberTimeBin *jtb = &fTB[0];
 168   for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++, ++jtb){
 169     if(!(Int_t(*jtb))) continue;
 170     (*jtb).SetOwner();
 171   }
 172 }
 173
 174 //_______________________________________________________
 175 Double_t AliTRDtrackingChamber::GetQuality()
 176 {
 177   //
 178   // Calculate chamber quality for seeding.
 179   //
 180   //
 181   // Parameters :
 182   //   layers : Array of propagation layers for this plane.
 183   //
 184   // Output :
 185   //   plane quality factor for seeding
 186   //
 187   // Detailed description
 188   //
 189   // The quality of the plane for seeding is higher if:
 190   //  1. the average timebin population is closer to an integer number
 191   //  2. the distribution of clusters/timebin is closer to a uniform distribution.
 192   //    - the slope of the first derivative of a parabolic fit is small or
 193   //    - the slope of a linear fit is small
 194   //
 195
 196   Int_t ncl   = 0;
 197   Int_t nused = 0;
 198   Int_t nClLayer;
 199   for(int itb=0; itb<AliTRDseedV1::kNtb; itb++){
 200     if(!(nClLayer = fTB[itb].GetNClusters())) continue;
 201     ncl += nClLayer;
 202     for(Int_t incl = 0; incl < nClLayer; incl++){
 203       if((fTB[itb].GetCluster(incl))->IsUsed()) nused++;
 204     }
 205   }
 206
 207   // calculate the deviation of the mean number of clusters from the
 208   // closest integer values
 209   Float_t nclMed = float(ncl-nused)/AliTRDtrackerV1::GetNTimeBins();
 210   Int_t ncli = Int_t(nclMed);
 211   Float_t nclDev = TMath::Abs(nclMed - TMath::Max(ncli, 1));
 212   nclDev -= (nclDev>.5) && ncli ? 1. : 0.;
 213   return TMath::Exp(-5.*TMath::Abs(nclDev));
 214
 215 //      // get slope of the derivative
 216 //      if(!fitter.Eval()) return quality;
 217 //      fitter.PrintResults(3);
 218 //      Double_t a = fitter.GetParameter(1);
 219 //
 220 //      printf("ncl_dev(%f)  a(%f)\n", ncl_dev, a);
 221 //      return quality*TMath::Exp(-a);
 222
 223 }
 224
 225
 226 //_______________________________________________________
 227 Bool_t AliTRDtrackingChamber::GetSeedingLayer(AliTRDchamberTimeBin *&fakeLayer, AliTRDgeometry * const geo, const AliTRDReconstructor *rec)
 228 {
 229   //
 230   // Creates a seeding layer
 231   //
 232
 233   // constants
 234   const Int_t kMaxRows = 16;
 235   const Int_t kMaxCols = 144;
 236   const Int_t kMaxPads = 2304;
 237   Int_t timeBinMin = rec->GetRecoParam()->GetNumberOfPresamples();
 238   Int_t timeBinMax = rec->GetRecoParam()->GetNumberOfPostsamples();
 239
 240   // Get the geometrical data of the chamber
 241   Int_t layer = geo->GetLayer(fDetector);
 242   Int_t stack = geo->GetStack(fDetector);
 243   Int_t sector= geo->GetSector(fDetector);
 244   AliTRDpadPlane *pp = geo->GetPadPlane(layer, stack);
 245   Int_t nCols = pp->GetNcols();
 246   Float_t ymin = TMath::Min(pp->GetCol0(), pp->GetColEnd());
 247   Float_t ymax = TMath::Max(pp->GetCol0(), pp->GetColEnd());
 248   Float_t zmin = TMath::Min(pp->GetRow0(), pp->GetRowEnd());
 249   Float_t zmax = TMath::Max(pp->GetRow0(), pp->GetRowEnd());
 250   Float_t z0 = -1., zl = -1.;
 251   Int_t nRows = pp->GetNrows();
 252   Float_t binlength = (ymax - ymin)/nCols;
 253   //AliInfo(Form("ymin(%f) ymax(%f) zmin(%f) zmax(%f) nRows(%d) binlength(%f)", ymin, ymax, zmin, zmax, nRows, binlength));
 254
 255   // Fill the histogram
 256   Int_t nClusters;
 257   Int_t *histogram[kMaxRows];                                                                                   // 2D-Histogram
 258   Int_t hvals[kMaxPads + 1];    memset(hvals, 0, sizeof(Int_t)*kMaxPads);        // one entry in addition for termination flag
 259   Float_t *sigmas[kMaxRows];
 260   Float_t svals[kMaxPads];      memset(svals, 0, sizeof(Float_t)*kMaxPads);
 261   AliTRDcluster *c = NULL;
 262   for(Int_t irs = 0; irs < kMaxRows; irs++){
 263     histogram[irs] = &hvals[irs*kMaxCols];
 264     sigmas[irs] = &svals[irs*kMaxCols];
 265   }
 266   for(Int_t iTime = timeBinMin; iTime < AliTRDseedV1::kNtb-timeBinMax; iTime++){
 267     if(!(nClusters = fTB[iTime].GetNClusters())) continue;
 268     z0 = fTB[iTime].GetZ0();
 269     zl = fTB[iTime].GetDZ0();
 270     for(Int_t incl = 0; incl < nClusters; incl++){
 271       c = fTB[iTime].GetCluster(incl);
 272       histogram[c->GetPadRow()][c->GetPadCol()]++;
 273       sigmas[c->GetPadRow()][c->GetPadCol()] += c->GetSigmaZ2();
 274     }
 275   }
 276
 277 // Now I have everything in the histogram, do the selection
 278   //Int_t nPads = nCols * nRows;
 279   // This is what we are interested in: The center of gravity of the best candidates
 280   Float_t cogyvals[kMaxPads]; memset(cogyvals, 0, sizeof(Float_t)*kMaxPads);
 281   Float_t cogzvals[kMaxPads]; memset(cogzvals, 0, sizeof(Float_t)*kMaxPads);
 282   Float_t *cogy[kMaxRows];
 283   Float_t *cogz[kMaxRows];
 284
 285   // Lookup-Table storing coordinates according to the bins
 286   Float_t yLengths[kMaxCols];
 287   Float_t zLengths[kMaxRows];
 288   for(Int_t icnt = 0; icnt < nCols; icnt++){
 289     yLengths[icnt] = pp->GetColPos(nCols - 1 - icnt) + binlength/2;
 290   }
 291   for(Int_t icnt = 0; icnt < nRows; icnt++){
 292     zLengths[icnt] = pp->GetRowPos(icnt) - pp->GetRowSize(icnt)/2;
 293   }
 294
 295   // A bitfield is used to mask the pads as usable
 296   Short_t mask[kMaxCols]; memset(mask, 0 ,sizeof(Short_t) * kMaxCols);//bool mvals[kMaxPads];
 297   for(UChar_t icount = 0; icount < nRows; icount++){
 298     cogy[icount] = &cogyvals[icount*kMaxCols];
 299     cogz[icount] = &cogzvals[icount*kMaxCols];
 300   }
 301   // In this array the array position of the best candidates will be stored
 302   Int_t   cand[AliTRDtrackerV1::kMaxTracksStack];
 303   Float_t sigcands[AliTRDtrackerV1::kMaxTracksStack];
 304
 305   // helper variables
 306   Int_t indices[kMaxPads]; memset(indices, -1, sizeof(Int_t)*kMaxPads);
 307   Int_t nCandidates = 0;
 308   Float_t norm, cogv;
 309   // histogram filled -> Select best bins
 310   Int_t nPads = nCols * nRows;
 311   // take out all the bins which have less than 3 entries (faster sorting)
 312   Int_t content[kMaxPads], dictionary[kMaxPads], nCont = 0, padnumber = 0;
 313   Int_t *iter = &hvals[0], *citer = &content[0], *diter =  &dictionary[0]; // iterators for preselection
 314   const Int_t threshold = 2;
 315   hvals[nPads] = -1; // termination for iterator
 316   do{
 317     if(*iter > threshold){
 318       *(citer++) = *iter;
 319       *(diter++) = padnumber;
 320       nCont++;
 321     }
 322     padnumber++;
 323   }while(*(++iter) != -1);
 324   TMath::Sort(nCont, content, indices);
 325
 326   Int_t col, row, lower, lower1, upper, upper1;
 327   for(Int_t ib = 0; ib < nCont; ib++){
 328     if(nCandidates >= AliTRDtrackerV1::kMaxTracksStack){
 329       AliWarning(Form("Number of seed candidates %d exceeded maximum allowed per stack %d", nCandidates, AliTRDtrackerV1::kMaxTracksStack));
 330       break;
 331     }
 332     // Positions
 333     row = dictionary[indices[ib]]/nCols;
 334     col = dictionary[indices[ib]]%nCols;
 335     // here will be the threshold condition:
 336     if((mask[col] & (1 << row)) != 0) continue;         // Pad is masked: continue
 337     //  if(histogram[row][col] < TMath::Max(threshold, 1)){     // of course at least one cluster is needed
 338     //          break;                  // number of clusters below threshold: break;
 339     //  }
 340     // passing: Mark the neighbors
 341     lower  = TMath::Max(col - 1, 0); upper  = TMath::Min(col + 2, nCols);
 342     lower1 = TMath::Max(row - 1, 0); upper1 = TMath::Min(row + 2, nCols);
 343     for(Int_t ic = lower; ic < upper; ++ic)
 344       for(Int_t ir = lower1; ir < upper1; ++ir){
 345         if(ic == col && ir == row) continue;
 346         mask[ic] |= (1 << ir);
 347       }
 348     // Storing the position in an array
 349     // testing for neigboring
 350     cogv = 0;
 351     norm = 0;
 352     lower = TMath::Max(col - 1, 0);
 353     upper = TMath::Min(col + 2, nCols);
 354     for(Int_t inb = lower; inb < upper; ++inb){
 355       cogv += yLengths[inb] * histogram[row][inb];
 356       norm += histogram[row][inb];
 357     }
 358     cogy[row][col] = cogv / norm;
 359     cogv = 0; norm = 0;
 360     lower = TMath::Max(row - 1, 0);
 361     upper = TMath::Min(row + 2, nRows);
 362     for(Int_t inb = lower; inb < upper; ++inb){
 363       cogv += zLengths[inb] * histogram[inb][col];
 364       norm += histogram[inb][col];
 365     }
 366     cogz[row][col] = Float_t(cogv) /  norm;
 367     // passed the filter
 368     cand[nCandidates] = row*nCols + col;        // store the position of a passig candidate into an Array
 369     sigcands[nCandidates] = sigmas[row][col] / histogram[row][col]; // never be a floating point exeption
 370     // Analysis output
 371     nCandidates++;
 372   }
 373   if(!nCandidates) return kFALSE;
 374
 375   Float_t pos[3], sig[2];
 376   Short_t signal[7]; memset(&signal[0], 0, 7*sizeof(Short_t));
 377
 378   new(fakeLayer) AliTRDchamberTimeBin(layer, stack, sector, z0, zl);
 379   fakeLayer->SetReconstructor(rec);
 380   fakeLayer->SetNRows(nRows);
 381   fakeLayer->SetOwner(kFALSE);
 382   if(nCandidates){
 383     UInt_t fakeIndex = 0;
 384     for(Int_t ican = 0; ican < nCandidates; ican++){
 385       row = cand[ican] / nCols;
 386       col = cand[ican] % nCols;
 387       //temporary
 388       Int_t n = 0; Double_t x = 0., y = 0., z = 0.;
 389       for(int itb=0; itb<AliTRDseedV1::kNtb; itb++){
 390         if(!(nClusters = fTB[itb].GetNClusters())) continue;
 391         for(Int_t incl = 0; incl < nClusters; incl++){
 392           c = fTB[itb].GetCluster(incl);
 393           if(c->GetPadRow() != row) continue;
 394           if(TMath::Abs(c->GetPadCol() - col) > 2) continue;
 395           x += c->GetX();
 396           y += c->GetY();
 397           z += c->GetZ();
 398           n++;
 399         }
 400       }
 401       pos[0] = x/n;
 402       pos[1] = y/n;
 403       pos[2] = z/n;
 404       sig[0] = .02;
 405       sig[1] = sigcands[ican];
 406       fakeLayer->InsertCluster(new AliTRDcluster(fDetector, 0., pos, sig, NULL, 3, signal, col, row, 0, 0, 0., 0), fakeIndex++);
 407     }
 408   }
 409   fakeLayer->BuildIndices();
 410   //fakeLayer->Print();
 411
 412   if(rec->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) >= 3){
 413     //TMatrixD hist(nRows, nCols);
 414     //for(Int_t i = 0; i < nRows; i++)
 415     //  for(Int_t j = 0; j < nCols; j++)
 416     //          hist(i,j) = histogram[i][j];
 417     TTreeSRedirector &cstreamer = *rec->GetDebugStream(AliTRDrecoParam::kTracker);
 418     cstreamer << "GetSeedingLayer"
 419     << "layer="      << layer
 420     << "ymin="       << ymin
 421     << "ymax="       << ymax
 422     << "zmin="       << zmin
 423     << "zmax="       << zmax
 424     << "L.="         << fakeLayer
 425     //<< "Histogram.=" << &hist
 426     << "\n";
 427   }
 428
 429   return kTRUE;
 430 }
 431
 432
 433 //_______________________________________________________
 434 void AliTRDtrackingChamber::Print(Option_t *opt) const
 435 {
 436   // Print the chamber status
 437   if(!GetNClusters()) return;
 438   AliInfo(Form("fDetector   = %d", fDetector));
 439   AliInfo(Form("fX0         = %7.3f", fX0));
 440   const AliTRDchamberTimeBin *itb = &fTB[0];
 441   for(Int_t jtb=0; jtb<AliTRDseedV1::kNtb; jtb++, itb++) (*itb).Print(opt);
 442 }
 443
 444
 445 //_______________________________________________________
 446 void AliTRDtrackingChamber::Update()
 447 {
 448 // Steer purging of used and shared clusters
 449
 450   AliTRDchamberTimeBin *jtb = &fTB[0];
 451   for(Int_t itb=AliTRDseedV1::kNtb; itb--; ++jtb){
 452     if(!(Int_t(*jtb))) continue;
 453     (*jtb).BuildIndices();
 454   }
 455 }
 456