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 "Cal/AliTRDCalDet.h"
  41 #include "Cal/AliTRDCalROC.h"
  42
  43 ClassImp(AliTRDtrackingChamber)
  44
  45 //_______________________________________________________
  46 AliTRDtrackingChamber::AliTRDtrackingChamber(Int_t det) :
  47   fDetector(det)
  48   ,fX0(0.)
  49 {}
  50
  51 //_______________________________________________________
  52 void AliTRDtrackingChamber::Clear(const Option_t *opt)
  53 {
  54   for(Int_t itb=0; itb<kNTimeBins; itb++) fTB[itb].Clear(opt);
  55 }
  56
  57 //_______________________________________________________
  58 void AliTRDtrackingChamber::InsertCluster(AliTRDcluster *c, Int_t index)
  59 {
  60   fTB[c->GetPadTime()].InsertCluster(c, index);
  61 }
  62
  63 //_______________________________________________________
  64 Bool_t AliTRDtrackingChamber::Build(AliTRDgeometry *geo, const AliTRDCalDet *cal, Bool_t hlt)
  65 {
  66 // Init chamber and all time bins (AliTRDchamberTimeBin)
  67 // Calculates radial position of the chamber based on
  68 // radial positions of the time bins (calibration/alignment aware)
  69 //
  70   Int_t stack = geo->GetStack(fDetector);
  71   Int_t layer = geo->GetLayer(fDetector);
  72   AliTRDpadPlane *pp = geo->GetPadPlane(layer, stack);
  73   Double_t zl = pp->GetRow0ROC() - pp->GetRowEndROC();
  74   Double_t z0 = geo->GetRow0(layer, stack, 0) - zl;
  75   Int_t nrows = pp->GetNrows();
  76
  77   Int_t index[50], jtb = 0;
  78   for(Int_t itb=0; itb<kNTimeBins; itb++){
  79     if(!fTB[itb]) continue;
  80     fTB[itb].SetRange(z0, zl);
  81     fTB[itb].SetNRows(nrows);
  82     fTB[itb].BuildIndices();
  83     index[jtb++] = itb;
  84   }
  85   if(jtb<2) return kFALSE;
  86
  87
  88   // ESTIMATE POSITION OF PAD PLANE FOR THIS CHAMBER
  89   Double_t x0 = fTB[index[0]].GetX();
  90   Double_t x1 = fTB[index[1]].GetX();
  91   Double_t dx = (x0 - x1)/(index[1] - index[0]);
  92
  93   Int_t t0 = (Int_t)cal->GetValue(fDetector);
  94   if(!hlt){
  95     Double_t mean = 0.0;
  96     AliTRDCalROC *roc = AliTRDcalibDB::Instance()->GetT0ROC(fDetector);
  97     for(Int_t k = 0; k<roc->GetNchannels(); k++) mean += roc->GetValue(k);
  98     mean /= roc->GetNchannels();
  99     t0 = (Int_t)(cal->GetValue(fDetector) + mean);
 100   }
 101
 102   fX0 = x0 + dx*(index[0] - t0);
 103   return kTRUE;
 104 }
 105
 106 //_______________________________________________________
 107 Int_t AliTRDtrackingChamber::GetNClusters() const
 108 {
 109 // Returns number of clusters in chamber
 110 //
 111   Int_t n = 0;
 112   for(Int_t itb=0; itb<kNTimeBins; itb++){
 113     n += Int_t(fTB[itb]);
 114   }
 115   return n;
 116 }
 117
 118 //_______________________________________________________
 119 Double_t AliTRDtrackingChamber::GetQuality()
 120 {
 121   //
 122   // Calculate chamber quality for seeding.
 123   //
 124   //
 125   // Parameters :
 126   //   layers : Array of propagation layers for this plane.
 127   //
 128   // Output :
 129   //   plane quality factor for seeding
 130   //
 131   // Detailed description
 132   //
 133   // The quality of the plane for seeding is higher if:
 134   //  1. the average timebin population is closer to an integer number
 135   //  2. the distribution of clusters/timebin is closer to a uniform distribution.
 136   //    - the slope of the first derivative of a parabolic fit is small or
 137   //    - the slope of a linear fit is small
 138   //
 139
 140   Int_t ncl   = 0;
 141   Int_t nused = 0;
 142   Int_t nClLayer;
 143   for(int itb=0; itb<kNTimeBins; itb++){
 144     if(!(nClLayer = fTB[itb].GetNClusters())) continue;
 145     ncl += nClLayer;
 146     for(Int_t incl = 0; incl < nClLayer; incl++){
 147       if((fTB[itb].GetCluster(incl))->IsUsed()) nused++;
 148     }
 149   }
 150
 151   // calculate the deviation of the mean number of clusters from the
 152   // closest integer values
 153   Float_t nclMed = float(ncl-nused)/AliTRDtrackerV1::GetNTimeBins();
 154   Int_t ncli = Int_t(nclMed);
 155   Float_t nclDev = TMath::Abs(nclMed - TMath::Max(ncli, 1));
 156   nclDev -= (nclDev>.5) && ncli ? 1. : 0.;
 157   return TMath::Exp(-5.*TMath::Abs(nclDev));
 158
 159 //      // get slope of the derivative
 160 //      if(!fitter.Eval()) return quality;
 161 //      fitter.PrintResults(3);
 162 //      Double_t a = fitter.GetParameter(1);
 163 //
 164 //      printf("ncl_dev(%f)  a(%f)\n", ncl_dev, a);
 165 //      return quality*TMath::Exp(-a);
 166
 167 }
 168
 169
 170 //_______________________________________________________
 171 Bool_t AliTRDtrackingChamber::GetSeedingLayer(AliTRDchamberTimeBin *&fakeLayer, AliTRDgeometry *geo, const AliTRDReconstructor *rec)
 172 {
 173   //
 174   // Creates a seeding layer
 175   //
 176
 177   // constants
 178   const Int_t kMaxRows = 16;
 179   const Int_t kMaxCols = 144;
 180   const Int_t kMaxPads = 2304;
 181   Int_t timeBinMin = rec->GetRecoParam()->GetNumberOfPresamples();
 182   Int_t timeBinMax = rec->GetRecoParam()->GetNumberOfPostsamples();
 183
 184   // Get the geometrical data of the chamber
 185   Int_t layer = geo->GetLayer(fDetector);
 186   Int_t stack = geo->GetStack(fDetector);
 187   Int_t sector= geo->GetSector(fDetector);
 188   AliTRDpadPlane *pp = geo->GetPadPlane(layer, stack);
 189   Int_t nCols = pp->GetNcols();
 190   Float_t ymin = TMath::Min(pp->GetCol0(), pp->GetColEnd());
 191   Float_t ymax = TMath::Max(pp->GetCol0(), pp->GetColEnd());
 192   Float_t zmin = TMath::Min(pp->GetRow0(), pp->GetRowEnd());
 193   Float_t zmax = TMath::Max(pp->GetRow0(), pp->GetRowEnd());
 194   Float_t z0 = -1., zl = -1.;
 195   Int_t nRows = pp->GetNrows();
 196   Float_t binlength = (ymax - ymin)/nCols;
 197   //AliInfo(Form("ymin(%f) ymax(%f) zmin(%f) zmax(%f) nRows(%d) binlength(%f)", ymin, ymax, zmin, zmax, nRows, binlength));
 198
 199   // Fill the histogram
 200   Int_t nClusters;
 201   Int_t *histogram[kMaxRows];                                                                                   // 2D-Histogram
 202   Int_t hvals[kMaxPads + 1];    memset(hvals, 0, sizeof(Int_t)*kMaxPads);        // one entry in addition for termination flag
 203   Float_t *sigmas[kMaxRows];
 204   Float_t svals[kMaxPads];      memset(svals, 0, sizeof(Float_t)*kMaxPads);
 205   AliTRDcluster *c = 0x0;
 206   for(Int_t irs = 0; irs < kMaxRows; irs++){
 207     histogram[irs] = &hvals[irs*kMaxCols];
 208     sigmas[irs] = &svals[irs*kMaxCols];
 209   }
 210   for(Int_t iTime = timeBinMin; iTime < kNTimeBins-timeBinMax; iTime++){
 211     if(!(nClusters = fTB[iTime].GetNClusters())) continue;
 212     z0 = fTB[iTime].GetZ0();
 213     zl = fTB[iTime].GetDZ0();
 214     for(Int_t incl = 0; incl < nClusters; incl++){
 215       c = fTB[iTime].GetCluster(incl);
 216       histogram[c->GetPadRow()][c->GetPadCol()]++;
 217       sigmas[c->GetPadRow()][c->GetPadCol()] += c->GetSigmaZ2();
 218     }
 219   }
 220
 221 // Now I have everything in the histogram, do the selection
 222   //Int_t nPads = nCols * nRows;
 223   // This is what we are interested in: The center of gravity of the best candidates
 224   Float_t cogyvals[kMaxPads]; memset(cogyvals, 0, sizeof(Float_t)*kMaxPads);
 225   Float_t cogzvals[kMaxPads]; memset(cogzvals, 0, sizeof(Float_t)*kMaxPads);
 226   Float_t *cogy[kMaxRows];
 227   Float_t *cogz[kMaxRows];
 228
 229   // Lookup-Table storing coordinates according to the bins
 230   Float_t yLengths[kMaxCols];
 231   Float_t zLengths[kMaxRows];
 232   for(Int_t icnt = 0; icnt < nCols; icnt++){
 233     yLengths[icnt] = pp->GetColPos(nCols - 1 - icnt) + binlength/2;
 234   }
 235   for(Int_t icnt = 0; icnt < nRows; icnt++){
 236     zLengths[icnt] = pp->GetRowPos(icnt) - pp->GetRowSize(icnt)/2;
 237   }
 238
 239   // A bitfield is used to mask the pads as usable
 240   Short_t mask[kMaxCols]; memset(mask, 0 ,sizeof(Short_t) * kMaxCols);//bool mvals[kMaxPads];
 241   for(UChar_t icount = 0; icount < nRows; icount++){
 242     cogy[icount] = &cogyvals[icount*kMaxCols];
 243     cogz[icount] = &cogzvals[icount*kMaxCols];
 244   }
 245   // In this array the array position of the best candidates will be stored
 246   Int_t   cand[AliTRDtrackerV1::kMaxTracksStack];
 247   Float_t sigcands[AliTRDtrackerV1::kMaxTracksStack];
 248
 249   // helper variables
 250   Int_t indices[kMaxPads]; memset(indices, -1, sizeof(Int_t)*kMaxPads);
 251   Int_t nCandidates = 0;
 252   Float_t norm, cogv;
 253   // histogram filled -> Select best bins
 254   Int_t nPads = nCols * nRows;
 255   // take out all the bins which have less than 3 entries (faster sorting)
 256   Int_t content[kMaxPads], dictionary[kMaxPads], nCont = 0, padnumber = 0;
 257   Int_t *iter = &hvals[0], *citer = &content[0], *diter =  &dictionary[0]; // iterators for preselection
 258   const Int_t threshold = 2;
 259   hvals[nPads] = -1; // termination for iterator
 260   do{
 261     if(*iter > threshold){
 262       *(citer++) = *iter;
 263       *(diter++) = padnumber;
 264       nCont++;
 265     }
 266     padnumber++;
 267   }while(*(++iter) != -1);
 268   TMath::Sort(nCont, content, indices);
 269
 270   Int_t col, row, lower, lower1, upper, upper1;
 271   for(Int_t ib = 0; ib < nCont; ib++){
 272     if(nCandidates >= AliTRDtrackerV1::kMaxTracksStack){
 273       printf("Number of seed candidates %d exceeded maximum allowed per stack %d", nCandidates, AliTRDtrackerV1::kMaxTracksStack);
 274       break;
 275     }
 276     // Positions
 277     row = dictionary[indices[ib]]/nCols;
 278     col = dictionary[indices[ib]]%nCols;
 279     // here will be the threshold condition:
 280     if((mask[col] & (1 << row)) != 0) continue;         // Pad is masked: continue
 281     //  if(histogram[row][col] < TMath::Max(threshold, 1)){     // of course at least one cluster is needed
 282     //          break;                  // number of clusters below threshold: break;
 283     //  }
 284     // passing: Mark the neighbors
 285     lower  = TMath::Max(col - 1, 0); upper  = TMath::Min(col + 2, nCols);
 286     lower1 = TMath::Max(row - 1, 0); upper1 = TMath::Min(row + 2, nCols);
 287     for(Int_t ic = lower; ic < upper; ++ic)
 288       for(Int_t ir = lower1; ir < upper1; ++ir){
 289         if(ic == col && ir == row) continue;
 290         mask[ic] |= (1 << ir);
 291       }
 292     // Storing the position in an array
 293     // testing for neigboring
 294     cogv = 0;
 295     norm = 0;
 296     lower = TMath::Max(col - 1, 0);
 297     upper = TMath::Min(col + 2, nCols);
 298     for(Int_t inb = lower; inb < upper; ++inb){
 299       cogv += yLengths[inb] * histogram[row][inb];
 300       norm += histogram[row][inb];
 301     }
 302     cogy[row][col] = cogv / norm;
 303     cogv = 0; norm = 0;
 304     lower = TMath::Max(row - 1, 0);
 305     upper = TMath::Min(row + 2, nRows);
 306     for(Int_t inb = lower; inb < upper; ++inb){
 307       cogv += zLengths[inb] * histogram[inb][col];
 308       norm += histogram[inb][col];
 309     }
 310     cogz[row][col] = Float_t(cogv) /  norm;
 311     // passed the filter
 312     cand[nCandidates] = row*nCols + col;        // store the position of a passig candidate into an Array
 313     sigcands[nCandidates] = sigmas[row][col] / histogram[row][col]; // never be a floating point exeption
 314     // Analysis output
 315     nCandidates++;
 316   }
 317   if(!nCandidates) return kFALSE;
 318
 319   Float_t pos[3], sig[2];
 320   Short_t signal[7]; memset(&signal[0], 0, 7*sizeof(Short_t));
 321
 322   new(fakeLayer) AliTRDchamberTimeBin(layer, stack, sector, z0, zl);
 323   fakeLayer->SetReconstructor(rec);
 324   AliTRDcluster *cluster = 0x0;
 325   if(nCandidates){
 326     UInt_t fakeIndex = 0;
 327     for(Int_t ican = 0; ican < nCandidates; ican++){
 328       row = cand[ican] / nCols;
 329       col = cand[ican] % nCols;
 330       //temporary
 331       Int_t n = 0; Double_t x = 0., y = 0., z = 0.;
 332       for(int itb=0; itb<kNTimeBins; itb++){
 333         if(!(nClusters = fTB[itb].GetNClusters())) continue;
 334         for(Int_t incl = 0; incl < nClusters; incl++){
 335           c = fTB[itb].GetCluster(incl);
 336           if(c->GetPadRow() != row) continue;
 337           if(TMath::Abs(c->GetPadCol() - col) > 2) continue;
 338           x += c->GetX();
 339           y += c->GetY();
 340           z += c->GetZ();
 341           n++;
 342         }
 343       }
 344       pos[0] = x/n;
 345       pos[1] = y/n;
 346       pos[2] = z/n;
 347       sig[0] = .02;
 348       sig[1] = sigcands[ican];
 349       cluster = new AliTRDcluster(fDetector, 0., pos, sig, 0x0, 3, signal, col, row, 0, 0, 0., 0);
 350       fakeLayer->InsertCluster(cluster, fakeIndex++);
 351     }
 352   }
 353   fakeLayer->SetNRows(nRows);
 354   fakeLayer->SetOwner();
 355   fakeLayer->BuildIndices();
 356   //fakeLayer->PrintClusters();
 357
 358   if(rec->GetStreamLevel(AliTRDReconstructor::kTracker) >= 3){
 359     //TMatrixD hist(nRows, nCols);
 360     //for(Int_t i = 0; i < nRows; i++)
 361     //  for(Int_t j = 0; j < nCols; j++)
 362     //          hist(i,j) = histogram[i][j];
 363     TTreeSRedirector &cstreamer = *AliTRDtrackerV1::DebugStreamer();
 364     cstreamer << "GetSeedingLayer"
 365     << "layer="      << layer
 366     << "ymin="       << ymin
 367     << "ymax="       << ymax
 368     << "zmin="       << zmin
 369     << "zmax="       << zmax
 370     << "L.="         << fakeLayer
 371     //<< "Histogram.=" << &hist
 372     << "\n";
 373   }
 374
 375   return kTRUE;
 376 }
 377