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