HLT/ITS/trackingSAP/AliITSSAPLayer.cxx

   1 #include <TMath.h>
   2 #include "AliITSSAPLayer.h"
   3 #include "AliITSSAPAux.h"
   4 #include "AliITSRecPoint.h"
   5 #include "AliITSgeomTGeo.h"
   6 #include "AliVertex.h"
   7 #include <TRandom.h>
   8 #include <TStopwatch.h>
   9 #include <TString.h>
  10
  11
  12 using namespace AliITSSAPAux;
  13
  14 //_________________________________________________________________
  15 AliITSSAPLayer::AliITSSAPLayer() :
  16   fClusters(0)
  17   ,fLrID(-1)
  18   ,fVIDOffset(0)
  19   ,fNClusters(0)
  20   ,fZMin(0)
  21   ,fZMax(0)
  22   ,fDZInv(-1)
  23   ,fDPhiInv(-1)
  24   ,fNZBins(20)
  25   ,fNPhiBins(20)
  26   ,fQueryZBmin(-1)
  27   ,fQueryZBmax(-1)
  28   ,fQueryPhiBmin(-1)
  29   ,fQueryPhiBmax(-1)
  30   ,fBins(0)
  31   ,fOccBins(0)
  32   ,fNOccBins(0)
  33   ,fNFoundClusters(0)
  34   ,fFoundClusterIterator(0)
  35   ,fFoundBinIterator(0)
  36   ,fFoundBins(0)
  37   ,fSortedClInfo(0)
  38 {
  39   // def. c-tor
  40 }
  41
  42 //_________________________________________________________________
  43 AliITSSAPLayer::AliITSSAPLayer(int id, float zspan,int nzbins,int nphibins, int buffer) :
  44   fClusters(0)
  45   ,fLrID(id)
  46   ,fVIDOffset((id+1)*2048)
  47   ,fNClusters(0)
  48   ,fZMin(-zspan)
  49   ,fZMax(zspan)
  50   ,fDZInv(-1)
  51   ,fDPhiInv(-1)
  52   ,fNZBins(nzbins)
  53   ,fNPhiBins(nphibins)
  54   ,fQueryZBmin(-1)
  55   ,fQueryZBmax(-1)
  56   ,fQueryPhiBmin(-1)
  57   ,fQueryPhiBmax(-1)
  58   ,fBins(0)
  59   ,fOccBins(0)
  60   ,fNOccBins(0)
  61   ,fNFoundClusters(0)
  62   ,fFoundClusterIterator(0)
  63   ,fFoundBinIterator(0)
  64   ,fFoundBins(0)
  65   ,fSortedClInfo(0)
  66 {
  67   // c-tor
  68   Init(buffer);
  69 }
  70
  71 //_________________________________________________________________
  72 AliITSSAPLayer::~AliITSSAPLayer()
  73 {
  74   // d-tor
  75   delete[] fBins;
  76   delete[] fOccBins;
  77   delete fClusters;
  78 }
  79
  80 //_________________________________________________________________
  81 void AliITSSAPLayer::Init(int buffer)
  82 {
  83   if (fClusters) {
  84     printf("Already initialized\n");
  85     return;
  86   }
  87   if (fNZBins<1)   fNZBins = 2;
  88   if (fNPhiBins<1) fNPhiBins = 1;
  89   fDZInv   = fNZBins/(fZMax-fZMin);
  90   fDPhiInv = fNPhiBins/TMath::TwoPi();
  91   //
  92   fBins = new ClBinInfo_t[fNZBins*fNPhiBins];
  93   fOccBins = new int[fNZBins*fNPhiBins];
  94   if (buffer<100) buffer = 100;
  95   fClusters = new TObjArray(buffer);
  96   fSortedClInfo.reserve(buffer);
  97   //
  98   // prepare detectors info
  99   int id1 = fLrID+1;
 100   Int_t nlad=AliITSgeomTGeo::GetNLadders(id1);
 101   Int_t ndet=AliITSgeomTGeo::GetNDetectors(id1);
 102   int detID = 0;
 103   for (Int_t j=1; j<nlad+1; j++) {
 104     for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
 105       ITSDetInfo_t det;
 106       det.index = detID++;
 107       //
 108       TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(id1,j,k,m);
 109       const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(id1,j,k);
 110       m.Multiply(tm);
 111       Double_t txyz[3] = {0.}, xyz[3] = {0.};
 112       m.LocalToMaster(txyz,xyz);
 113       det.xTF = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 114       det.phiTF = TMath::ATan2(xyz[1],xyz[0]);
 115       //BringTo02Pi(det.phiTF);
 116       det.sinTF = TMath::Sin(det.phiTF);
 117       det.cosTF = TMath::Cos(det.phiTF);
 118       //
 119       // compute the real radius (with misalignment)
 120       TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(id1,j,k)));
 121       mmisal.Multiply(tm);
 122       xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
 123       mmisal.LocalToMaster(txyz,xyz);
 124       det.xTFmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 125       //
 126       fDetectors.push_back(det);
 127     } // end loop on detectors
 128   } // end loop on ladders
 129
 130 }
 131
 132 //_________________________________________________________________
 133 void AliITSSAPLayer::SortClusters(const AliVertex* vtx)
 134 {
 135   // sort clusters and build fast lookup table
 136   //
 137   ClearSortedInfo();
 138   fSortedClInfo.reserve(fNClusters);
 139   //
 140   ClsInfo_t cl;
 141   for (int icl=fNClusters;icl--;) {
 142     AliITSRecPoint* cluster = (AliITSRecPoint*)fClusters->UncheckedAt(icl);
 143     cluster->GetGlobalXYZ( (float*)&cl );
 144     //
 145     if (vtx) { // phi and r will be computed wrt vertex
 146       cl.x -= vtx->GetX();
 147       cl.y -= vtx->GetY();
 148     }
 149     //
 150     cl.r = TMath::Sqrt(cl.x*cl.x + cl.y*cl.y);
 151     cl.phi = TMath::ATan2(cl.y,cl.x);
 152     BringTo02Pi(cl.phi);
 153     cl.index = icl;
 154     cl.zphibin = GetBinIndex(GetZBin(cl.z),GetPhiBin(cl.phi));
 155     cl.detid = cluster->GetVolumeId() - fVIDOffset;
 156     //
 157     fSortedClInfo.push_back( cl );
 158     //
 159   }
 160   sort(fSortedClInfo.begin(), fSortedClInfo.end()); // sort in phi, z
 161   //
 162   // fill cells in phi,z
 163   int currBin = -1;
 164   for (int icl=0;icl<fNClusters;icl++) {
 165     ClsInfo_t &t = fSortedClInfo[icl];
 166     if (t.zphibin>currBin) { // register new occupied bin
 167       currBin = t.zphibin;
 168       fBins[currBin].first = icl;
 169       fBins[currBin].index = fNOccBins;
 170       fOccBins[fNOccBins++] = currBin;
 171     }
 172     fBins[currBin].ncl++;
 173   }
 174   //  Print("clb"); //RS
 175 }
 176
 177 //_________________________________________________________________
 178 void AliITSSAPLayer::Clear(Option_t *)
 179 {
 180   // clear cluster info
 181   ClearSortedInfo();
 182   fNClusters = 0;
 183   if (fClusters) fClusters->Clear();
 184   //
 185 }
 186
 187 //_________________________________________________________________
 188 void AliITSSAPLayer::ClearSortedInfo()
 189 {
 190   // clear cluster info
 191   fSortedClInfo.clear();
 192   memset(fBins,0,fNZBins*fNPhiBins*sizeof(ClBinInfo_t));
 193   memset(fOccBins,0,fNZBins*fNPhiBins*sizeof(int));
 194   fNOccBins = 0;
 195 }
 196
 197 //_________________________________________________________________
 198 void AliITSSAPLayer::Print(Option_t *opt) const
 199 {
 200   // dump cluster bins info
 201   TString opts = opt;
 202   opts.ToLower();
 203   printf("Stored %d clusters in %d occupied bins\n",fNClusters,fNOccBins);
 204   //
 205   if (opts.Contains("c")) {
 206     printf("\nCluster info\n");
 207     for (int i=0;i<fNClusters;i++) {
 208       const ClsInfo_t &t = fSortedClInfo[i];
 209       printf("#%5d Bin(phi/z):%03d/%03d Z:%+8.3f Phi:%+6.3f R:%7.3f Ind:%d ",
 210              i,t.zphibin/fNZBins,t.zphibin%fNZBins,t.z,t.phi,t.r,t.index);
 211       if (opts.Contains("l")) { // mc labels
 212         AliITSRecPoint* rp = (AliITSRecPoint*)fClusters->UncheckedAt(t.index);
 213         for (int l=0;l<3;l++) if (rp->GetLabel(l)>=0) printf("| %d ",rp->GetLabel(l));
 214       }
 215       printf("\n");
 216     }
 217   }
 218   //
 219   if (opts.Contains("b")) {
 220     printf("\nBins info (occupied only)\n");
 221     for (int i=0;i<fNOccBins;i++) {
 222       printf("%4d %5d(phi/z: %03d/%03d) -> %3d cl from %d\n",i,fOccBins[i],fOccBins[i]/fNZBins,fOccBins[i]%fNZBins,
 223              fBins[fOccBins[i]].ncl,fBins[fOccBins[i]].first);
 224     }
 225   }
 226   //
 227 }
 228
 229 //_____________________________________________________________
 230 int AliITSSAPLayer::SelectClusters(float zmin,float zmax,float phimin,float phimax)
 231 {
 232   // prepare occupied bins in the requested region
 233   //printf("Select: Z %f %f | Phi: %f %f\n",zmin,zmax,phimin,phimax);
 234   if (!fNOccBins) return 0;
 235   if (zmax<fZMin || zmin>fZMax || zmin>zmax) return 0;
 236   fFoundBins.clear();
 237
 238   fQueryZBmin = GetZBin(zmin);
 239   if (fQueryZBmin<0) fQueryZBmin = 0;
 240   fQueryZBmax = GetZBin(zmax);
 241   if (fQueryZBmax>=fNZBins) fQueryZBmax = fNZBins-1;
 242   BringTo02Pi(phimin);
 243   BringTo02Pi(phimax);
 244   fQueryPhiBmin = GetPhiBin(phimin);
 245   fQueryPhiBmax = GetPhiBin(phimax);
 246   int dbz=0;
 247   fNFoundClusters = 0;
 248   int nbcheck = fQueryPhiBmax - fQueryPhiBmin + 1;
 249   if (nbcheck>0) { // no wrapping around 0-2pi, fast case
 250     for (int ip=fQueryPhiBmin;ip<=fQueryPhiBmax;ip++) {
 251       int binID = GetBinIndex(fQueryZBmin,ip);
 252       if ( !(dbz=(fQueryZBmax-fQueryZBmin)) ) { // just one Z bin in the query range
 253         ClBinInfo_t& binInfo = fBins[binID];
 254         if (!binInfo.ncl) continue;
 255         fNFoundClusters += binInfo.ncl;
 256         fFoundBins.push_back(binID);
 257         continue;
 258       }
 259       int binMax = binID+dbz;
 260       for (;binID<=binMax;binID++) {
 261         ClBinInfo_t& binInfo = fBins[binID];
 262         if (!binInfo.ncl) continue;
 263         fNFoundClusters += binInfo.ncl;
 264         fFoundBins.push_back(binID);
 265       }
 266     }
 267   }
 268   else {  // wrapping
 269     nbcheck += fNPhiBins;
 270     for (int ip0=0;ip0<=nbcheck;ip0++) {
 271       int ip = fQueryPhiBmin+ip0;
 272       if (ip>=fNPhiBins) ip -= fNPhiBins;
 273       int binID = GetBinIndex(fQueryZBmin,ip);
 274       if ( !(dbz=(fQueryZBmax-fQueryZBmin)) ) { // just one Z bin in the query range
 275         ClBinInfo_t& binInfo = fBins[binID];
 276         if (!binInfo.ncl) continue;
 277         fNFoundClusters += binInfo.ncl;
 278         fFoundBins.push_back(binID);
 279         continue;
 280       }
 281       int binMax = binID+dbz;
 282       for (;binID<=binMax;binID++) {
 283         ClBinInfo_t& binInfo = fBins[binID];
 284         if (!binInfo.ncl) continue;
 285         fNFoundClusters += binInfo.ncl;
 286         fFoundBins.push_back(binID);
 287       }
 288     }
 289   }
 290   fFoundClusterIterator = fFoundBinIterator = 0;
 291   /*
 292   //printf("Selected -> %d cl in %d bins\n",fNFoundClusters,(int)fFoundBins.size());
 293   for (int i=0;i<(int)fFoundBins.size();i++) {
 294   int bn = fFoundBins[i];
 295   ClBinInfo_t& bin=fBins[bn];
 296   printf("#%d b:%d 1st: %3d Ncl:%d\n",i,bn,bin.first,bin.ncl);
 297   }
 298   printf("\n");
 299   */
 300   return fNFoundClusters;
 301 }
 302
 303 //_____________________________________________________________
 304 int AliITSSAPLayer::GetNextClusterInfoID()
 305 {
 306   if (fFoundBinIterator<0) return 0;
 307   int currBin = fFoundBins[fFoundBinIterator];
 308   if (fFoundClusterIterator<fBins[currBin].ncl) { // same bin
 309     return fBins[currBin].first+fFoundClusterIterator++;
 310   }
 311   if (++fFoundBinIterator<int(fFoundBins.size())) {  // need to change bin
 312     currBin = fFoundBins[fFoundBinIterator];
 313     fFoundClusterIterator = 1;
 314     return fBins[currBin].first;
 315   }
 316   fFoundBinIterator = -1;
 317   return -1;
 318 }
 319
 320 //_____________________________________________________________
 321 void AliITSSAPLayer::ResetFoundIterator()
 322 {
 323   // prepare for a new loop over found clusters
 324   if (fNFoundClusters)  fFoundClusterIterator = fFoundBinIterator = 0;
 325 }