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