ITS/UPGRADE/AliITSUCATrackingStation.cxx

   1 //-*- Mode: C++ -*-
   2 // **************************************************************************
   3 // This file is property of and copyright by the ALICE ITSU Project         *
   4 // ALICE Experiment at CERN, All rights reserved.                           *
   5 //                                                                          *
   6 // Primary Author: Ruben Shahoyan                                           *
   7 //                                                                          *
   8 // Adapted to ITSU: Maximiliano Puccio <maximiliano.puccio@cern.ch>         *
   9 //                  for the ITS Upgrade project                             *
  10 //                                                                          *
  11 // Permission to use, copy, modify and distribute this software and its     *
  12 // documentation strictly for non-commercial purposes is hereby granted     *
  13 // without fee, provided that the above copyright notice appears in all     *
  14 // copies and that both the copyright notice and this permission notice     *
  15 // appear in the supporting documentation. The authors make no claims       *
  16 // about the suitability of this software for any purpose. It is            *
  17 // provided "as is" without express or implied warranty.                    *
  18 //                                                                          *
  19 //***************************************************************************
  20
  21 #include "AliITSUCATrackingStation.h"
  22 #include <TMath.h>
  23 #include "AliITSUGeomTGeo.h"
  24 #include "AliVertex.h"
  25 #include <TRandom.h>
  26 #include <TStopwatch.h>
  27 #include <TString.h>
  28 #include "AliITSUAux.h"
  29 #include "AliITSURecoSens.h"
  30 #include <Riostream.h>
  31
  32 using AliITSUAux::BringTo02Pi;
  33
  34 //_________________________________________________________________
  35 AliITSUCATrackingStation::AliITSUCATrackingStation()
  36 :fClusters(0x0)
  37 ,fVIDOffset(0)
  38 ,fNClusters(0)
  39 ,fZMin(0)
  40 ,fZMax(0)
  41 ,fDZInv(-1)
  42 ,fDPhiInv(-1)
  43 ,fNZBins(20)
  44 ,fNPhiBins(20)
  45 ,fQueryZBmin(-1)
  46 ,fQueryZBmax(-1)
  47 ,fQueryPhiBmin(-1)
  48 ,fQueryPhiBmax(-1)
  49 ,fBins(0)
  50 ,fOccBins(0)
  51 ,fNOccBins(0)
  52 ,fNFoundClusters(0)
  53 ,fFoundClusterIterator(0)
  54 ,fFoundBinIterator(0)
  55 ,fIndex()
  56 ,fFoundBins(0)
  57 ,fSortedClInfo(0)
  58 {
  59   // def. c-tor
  60 }
  61
  62 //_________________________________________________________________
  63 AliITSUCATrackingStation::AliITSUCATrackingStation(int nzbins,int nphibins,
  64                                                    AliITSURecoLayer *lr, AliITSUGeomTGeo *geo)
  65 :fClusters(0x0)
  66 ,fVIDOffset()
  67 ,fNClusters(0)
  68 ,fZMin(0.f)
  69 ,fZMax(0.f)
  70 ,fDZInv(-1)
  71 ,fDPhiInv(-1)
  72 ,fNZBins(nzbins)
  73 ,fNPhiBins(nphibins)
  74 ,fQueryZBmin(-1)
  75 ,fQueryZBmax(-1)
  76 ,fQueryPhiBmin(-1)
  77 ,fQueryPhiBmax(-1)
  78 ,fBins(0)
  79 ,fOccBins(0)
  80 ,fNOccBins(0)
  81 ,fNFoundClusters(0)
  82 ,fFoundClusterIterator(0)
  83 ,fFoundBinIterator(0)
  84 ,fIndex()
  85 ,fFoundBins(0)
  86 ,fSortedClInfo(0)
  87 {
  88   // c-tor
  89   Init(lr,geo);
  90 }
  91
  92 //_________________________________________________________________
  93 AliITSUCATrackingStation::~AliITSUCATrackingStation()
  94 {
  95   // d-tor
  96   delete[] fBins;
  97   delete[] fOccBins;
  98 //  delete fClusters;
  99 }
 100
 101 //_________________________________________________________________
 102 void AliITSUCATrackingStation::Init(AliITSURecoLayer *lr, AliITSUGeomTGeo *geo)
 103 {
 104   if (fClusters) {
 105     printf("Already initialized\n");
 106     return;
 107   }
 108
 109   fClusters = *lr->GetClustersAddress();
 110   fZMin = lr->GetZMin();
 111   fZMax = lr->GetZMax();
 112   if (fNZBins < 1)   fNZBins = 2;
 113   if (fNPhiBins < 1) fNPhiBins = 1;
 114   fDZInv   = fNZBins / (fZMax - fZMin);
 115   fDPhiInv = fNPhiBins / TMath::TwoPi();
 116   //
 117   fBins = new ClBinInfo_t[fNZBins * fNPhiBins];
 118   fOccBins = new int[fNZBins * fNPhiBins];
 119   fNClusters = fClusters->GetEntriesFast();
 120   fSortedClInfo.reserve(fClusters->GetEntriesFast());
 121   fVIDOffset = ((AliITSUClusterPix*)fClusters->UncheckedAt(0))->GetVolumeId();
 122   //
 123   // prepare detectors info
 124   int detID = -1;
 125   fIndex.resize(lr->GetNSensors(),-1);
 126   fDetectors.reserve(lr->GetNSensors());
 127   for (int iCl = 0; iCl < fClusters->GetEntriesFast(); ++iCl)
 128   { //Fill this layer with detectors
 129     AliITSUClusterPix* c = (AliITSUClusterPix*)fClusters->UncheckedAt(iCl);
 130     if (detID == c->GetVolumeId())
 131     {
 132       continue;
 133     }
 134     detID = c->GetVolumeId();
 135     ITSDetInfo_t det;
 136     det.index = iCl;
 137     //
 138     TGeoHMatrix m;
 139     geo->GetOrigMatrix(detID,m);
 140     //
 141     fIndex[detID - fVIDOffset] = fDetectors.size();
 142     AliITSURecoSens* sens = lr->GetSensorFromID(detID);
 143     const TGeoHMatrix *tm = geo->GetMatrixT2L(detID);
 144     m.Multiply(tm);
 145     double txyz[3] = {0.,0.,0.}, xyz[3] = {0.,0.,0.};
 146     m.LocalToMaster(txyz,xyz);
 147     det.xTF = sens->GetXTF(); // TMath::Sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1]);
 148
 149     det.phiTF = sens->GetPhiTF();//TMath::ATan2(xyz[1],xyz[0]);
 150     det.sinTF = TMath::Sin(det.phiTF);
 151     det.cosTF = TMath::Cos(det.phiTF);
 152     //
 153     // compute the real radius (with misalignment)
 154     TGeoHMatrix mmisal(*(geo->GetMatrix(detID)));
 155     mmisal.Multiply(tm);
 156     xyz[0] = 0.;
 157     xyz[1] = 0.;
 158     xyz[2] = 0.;
 159     mmisal.LocalToMaster(txyz,xyz);
 160     det.xTFmisal = TMath::Sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1]);
 161 //    if (!TMath::AreEqualAbs(det.xTFmisal, sens->GetXTF(), 1e-9)) {
 162 //      cout << "Error: " << lr->GetActiveID() << " " << detID << " " << det.xTFmisal - sens->GetXTF()<< endl;
 163 //    }
 164
 165     fDetectors.push_back(det);
 166   } // end loop on detectors
 167 }
 168
 169 //_________________________________________________________________
 170 void AliITSUCATrackingStation::SortClusters(const AliVertex* vtx)
 171 {
 172   // sort clusters and build fast lookup table
 173   //
 174   ClearSortedInfo();
 175   fSortedClInfo.reserve(fNClusters);
 176   //
 177   ClsInfo_t cl;
 178   for (int icl = fNClusters;icl--;) {
 179     AliITSUClusterPix* cluster = (AliITSUClusterPix*)fClusters->UncheckedAt(icl);
 180     cluster->GetGlobalXYZ( (float*)&cl );
 181     //
 182     if (vtx) { // phi and r will be computed wrt vertex
 183       cl.x -= vtx->GetX();
 184       cl.y -= vtx->GetY();
 185     }
 186     //
 187     cl.r = TMath::Sqrt(cl.x*cl.x + cl.y*cl.y);
 188     cl.phi = TMath::ATan2(cl.y,cl.x);
 189     BringTo02Pi(cl.phi);
 190     cl.index = icl;
 191     cl.zphibin = GetBinIndex(GetZBin(cl.z),GetPhiBin(cl.phi));
 192     cl.detid = cluster->GetVolumeId() - fVIDOffset;
 193     //
 194     fSortedClInfo.push_back( cl );
 195     //
 196   }
 197   sort(fSortedClInfo.begin(), fSortedClInfo.end()); // sort in phi, z
 198   //
 199   // fill cells in phi,z
 200   int currBin = -1;
 201   for (int icl = 0;icl < fNClusters; ++icl)
 202   {
 203     ClsInfo_t &t = fSortedClInfo[icl];
 204     if (t.zphibin > currBin)
 205     { // register new occupied bin
 206       currBin = t.zphibin;
 207       fBins[currBin].first = icl;
 208       fBins[currBin].index = fNOccBins;
 209       fOccBins[fNOccBins++] = currBin;
 210     }
 211     fBins[currBin].ncl++;
 212   }
 213   //  Print("clb"); //RS
 214 }
 215
 216 //_________________________________________________________________
 217 void AliITSUCATrackingStation::Clear(Option_t *)
 218 {
 219   // clear cluster info
 220   ClearSortedInfo();
 221   fIndex.clear();
 222   fNClusters = 0;
 223   if (fClusters) fClusters = 0x0;
 224   //
 225 }
 226
 227 //_________________________________________________________________
 228 void AliITSUCATrackingStation::ClearSortedInfo()
 229 {
 230   // clear cluster info
 231   fSortedClInfo.clear();
 232   memset(fBins,0,fNZBins * fNPhiBins * sizeof(ClBinInfo_t));
 233   memset(fOccBins,0,fNZBins * fNPhiBins * sizeof(int));
 234   fNOccBins = 0;
 235 }
 236
 237 //_________________________________________________________________
 238 void AliITSUCATrackingStation::Print(Option_t *opt) const
 239 {
 240   // dump cluster bins info
 241   TString opts = opt;
 242   opts.ToLower();
 243   printf("Stored %d clusters in %d occupied bins\n",fNClusters,fNOccBins);
 244   //
 245   if (opts.Contains("c")) {
 246     printf("\nCluster info\n");
 247     for (int i = 0; i < fNClusters;i++)
 248     {
 249       const ClsInfo_t &t = fSortedClInfo[i];
 250       printf("#%5d Bin(phi/z):%03d/%03d Z:%+8.3f Phi:%+6.3f R:%7.3f Ind:%d ",
 251              i,t.zphibin/fNZBins,t.zphibin%fNZBins,t.z,t.phi,t.r,t.index);
 252       if (opts.Contains("l")) { // mc labels
 253         AliITSUClusterPix* rp = (AliITSUClusterPix*)fClusters->UncheckedAt(t.index);
 254         for (int l = 0;l < 3; l++) if (rp->GetLabel(l) >= 0) printf("| %d ",rp->GetLabel(l));
 255       }
 256       printf("\n");
 257     }
 258   }
 259   //
 260   if (opts.Contains("b")) {
 261     printf("\nBins info (occupied only)\n");
 262     for (int i=0;i<fNOccBins;i++) {
 263       printf("%4d %5d(phi/z: %03d/%03d) -> %3d cl from %d\n",i,fOccBins[i],fOccBins[i]/fNZBins,fOccBins[i]%fNZBins,
 264              fBins[fOccBins[i]].ncl,fBins[fOccBins[i]].first);
 265     }
 266   }
 267   //
 268 }
 269
 270 //_____________________________________________________________
 271 int AliITSUCATrackingStation::SelectClusters(float zmin,float zmax,float phimin,float phimax)
 272 {
 273   // prepare occupied bins in the requested region
 274   //printf("Select: Z %f %f | Phi: %f %f\n",zmin,zmax,phimin,phimax);
 275   if (!fNOccBins) return 0;
 276   if (zmax < fZMin || zmin > fZMax || zmin > zmax) return 0;
 277   fFoundBins.clear();
 278
 279   fQueryZBmin = GetZBin(zmin);
 280   if (fQueryZBmin < 0) fQueryZBmin = 0;
 281   fQueryZBmax = GetZBin(zmax);
 282   if (fQueryZBmax >= fNZBins) fQueryZBmax = fNZBins - 1;
 283   BringTo02Pi(phimin);
 284   BringTo02Pi(phimax);
 285   fQueryPhiBmin = GetPhiBin(phimin);
 286   fQueryPhiBmax = GetPhiBin(phimax);
 287   int dbz = 0;
 288   fNFoundClusters = 0;
 289   int nbcheck = fQueryPhiBmax - fQueryPhiBmin + 1; //TODO:(MP) check if a circular buffer is feasible
 290   if (nbcheck > 0)
 291   { // no wrapping around 0-2pi, fast case
 292     for (int ip = fQueryPhiBmin;ip <= fQueryPhiBmax;ip++) {
 293       int binID = GetBinIndex(fQueryZBmin,ip);
 294       if ( !(dbz = (fQueryZBmax-fQueryZBmin)) ) { // just one Z bin in the query range
 295         ClBinInfo_t& binInfo = fBins[binID];
 296         if (!binInfo.ncl) continue;
 297         fNFoundClusters += binInfo.ncl;
 298         fFoundBins.push_back(binID);
 299         continue;
 300       }
 301       int binMax = binID + dbz;
 302       for ( ; binID <= binMax; binID++) {
 303         ClBinInfo_t& binInfo = fBins[binID];
 304         if (!binInfo.ncl) continue;
 305         fNFoundClusters += binInfo.ncl;
 306         fFoundBins.push_back(binID);
 307       }
 308     }
 309   }
 310   else
 311   {  // wrapping
 312     nbcheck += fNPhiBins;
 313     for (int ip0 = 0;ip0 <= nbcheck;ip0++) {
 314       int ip = fQueryPhiBmin + ip0;
 315       if (ip >= fNPhiBins) ip -= fNPhiBins;
 316       int binID = GetBinIndex(fQueryZBmin,ip);
 317       if ( !(dbz = (fQueryZBmax - fQueryZBmin)) ) { // just one Z bin in the query range
 318         ClBinInfo_t& binInfo = fBins[binID];
 319         if (!binInfo.ncl) continue;
 320         fNFoundClusters += binInfo.ncl;
 321         fFoundBins.push_back(binID);
 322         continue;
 323       }
 324       int binMax = binID + dbz;
 325       for (;binID <= binMax;binID++) {
 326         ClBinInfo_t& binInfo = fBins[binID];
 327         if (!binInfo.ncl) continue;
 328         fNFoundClusters += binInfo.ncl;
 329         fFoundBins.push_back(binID);
 330       }
 331     }
 332   }
 333   fFoundClusterIterator = fFoundBinIterator = 0;
 334   /*
 335    //printf("Selected -> %d cl in %d bins\n",fNFoundClusters,(int)fFoundBins.size());
 336    for (int i=0;i<(int)fFoundBins.size();i++) {
 337    int bn = fFoundBins[i];
 338    ClBinInfo_t& bin=fBins[bn];
 339    printf("#%d b:%d 1st: %3d Ncl:%d\n",i,bn,bin.first,bin.ncl);
 340    }
 341    printf("\n");
 342    */
 343   return fNFoundClusters;
 344 }
 345
 346 //_____________________________________________________________
 347 int AliITSUCATrackingStation::GetNextClusterInfoID()
 348 {
 349   if (fFoundBinIterator < 0) return 0;
 350   int currBin = fFoundBins[fFoundBinIterator];
 351   if (fFoundClusterIterator < fBins[currBin].ncl) { // same bin
 352     return fBins[currBin].first + fFoundClusterIterator++;
 353   }
 354   if (++fFoundBinIterator < int(fFoundBins.size())) {  // need to change bin
 355     currBin = fFoundBins[fFoundBinIterator];
 356     fFoundClusterIterator = 1;
 357     return fBins[currBin].first;
 358   }
 359   fFoundBinIterator = -1;
 360   return -1;
 361 }
 362
 363 //_____________________________________________________________
 364 void AliITSUCATrackingStation::ResetFoundIterator()
 365 {
 366   // prepare for a new loop over found clusters
 367   if (fNFoundClusters)  fFoundClusterIterator = fFoundBinIterator = 0;
 368 }