TRD/AliTRDtrackingChamber.cxx

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