ITS/UPGRADE/AliITSURecoLayer.cxx

   1 #include <TClonesArray.h>
   2 #include "AliITSURecoLayer.h"
   3 #include "AliITSsegmentation.h"
   4 #include "AliITSUAux.h"
   5 #include "AliITSUClusterPix.h"
   6 #include "AliITSUGeomTGeo.h"
   7 #include "AliLog.h"
   8
   9 using namespace AliITSUAux;
  10 using namespace TMath;
  11
  12 ClassImp(AliITSURecoLayer)
  13
  14
  15 //______________________________________________________
  16 AliITSURecoLayer::AliITSURecoLayer(const char* name)
  17   :fActiveID(-1)
  18   ,fNSensors(0)
  19   ,fNSensInLadder(0)
  20   ,fNLadders(0)
  21   ,fR(0)
  22   ,fRMax(0)
  23   ,fRMin(0)
  24   ,fZMax(0)
  25   ,fZMin(0)
  26   ,fPhiLadMax(0)
  27   ,fPhiLadMin(0)
  28   ,fPhiOffs(0)
  29   ,fSensDZInv(0)
  30   ,fDPhiLadInv(0)
  31   ,fMaxStep(0.5)
  32   ,fSensors(0)
  33   ,fITSGeom(0)
  34   ,fClusters(0)
  35 {
  36   // def. c-tor
  37   SetNameTitle(name,name);
  38 }
  39
  40 //______________________________________________________
  41 AliITSURecoLayer::AliITSURecoLayer(const char* name, Int_t activeID, AliITSUGeomTGeo* gm)
  42   :fActiveID(activeID)
  43   ,fNSensors(0)
  44   ,fNSensInLadder(0)
  45   ,fNLadders(0)
  46   ,fR(0)
  47   ,fRMax(0)
  48   ,fRMin(0)
  49   ,fZMax(0)
  50   ,fZMin(0)
  51   ,fPhiLadMax(0)
  52   ,fPhiLadMin(0)
  53   ,fPhiOffs(0)
  54   ,fSensDZInv(0)
  55   ,fDPhiLadInv(0)
  56   ,fMaxStep(0.5)
  57   ,fSensors(0)
  58   ,fITSGeom(gm)
  59   ,fClusters(0)
  60 {
  61   // def. c-tor
  62   SetNameTitle(name,name);
  63   Build();
  64 }
  65
  66 //______________________________________________________
  67 AliITSURecoLayer::~AliITSURecoLayer()
  68 {
  69   // def. d-tor
  70   delete[] fSensors;
  71   delete[] fPhiLadMax;
  72   delete[] fPhiLadMin;
  73   if (GetOwnsClusterArray()) delete fClusters;
  74 }
  75
  76 //______________________________________________________
  77 void AliITSURecoLayer::Print(Option_t* opt) const
  78 {
  79   //print
  80   printf("Lr %-15s %d (act:%+d), NSens: %4d | MaxStep:%.2f ",GetName(),GetID(),GetActiveID(),GetNSensors(),fMaxStep);
  81   printf("%6.3f<R<%6.3f | %+8.3f<Z<%+8.3f dZ:%6.3f\n",fRMin,fRMax,fZMin,fZMax,fSensDZInv>0 ? 1/fSensDZInv : 0);
  82   TString opts = opt; opts.ToLower();
  83   if (opts.Contains("sn")) for (int i=0;i<GetNSensors();i++) GetSensor(i)->Print(opt);
  84 }
  85
  86 //______________________________________________________
  87 void AliITSURecoLayer::Build()
  88 {
  89   // build internal structures
  90   const double kSafeR = 0.05; // safety margin for Rmin,Rmax of the layer
  91   if (fActiveID<0) return;
  92   fNLadders = fITSGeom->GetNLadders(fActiveID);
  93   fNSensInLadder = fITSGeom->GetNDetectors(fActiveID);
  94   fNSensors = fNLadders*fNSensInLadder;
  95   fSensors = new AliITSURecoSens*[fNSensors];
  96   const AliITSsegmentation* kSegm = fITSGeom->GetSegmentation(fActiveID);
  97   //
  98   // name layer according its active id, detector type and segmentation tyoe
  99   TGeoHMatrix mmod;
 100   const TGeoHMatrix* mt2l;
 101   fRMin=fZMin=1e9;
 102   fRMax=fZMax=-1e9;
 103   double phiTF,rTF, loc[3]={0,0,0},glo[3];
 104   fNSensors = 0;
 105   fPhiLadMin = new Double_t[fNLadders];
 106   fPhiLadMax = new Double_t[fNLadders];
 107   fSensDZInv = 0;
 108   fDPhiLadInv = fNLadders/TwoPi();
 109   //
 110   for (int ild=0;ild<fNLadders;ild++) {
 111     fPhiLadMin[ild] = 1e9;
 112     fPhiLadMax[ild] = -1e9;
 113     //
 114     for (int idt=0;idt<fNSensInLadder;idt++) {
 115       AliITSURecoSens* sens = new AliITSURecoSens(fNSensors++);
 116       fSensors[ild*fNSensInLadder+idt] = sens;
 117       //
 118       double phiMin=1e9,phiMax=-1e9,zMin=1e9,zMax=-1e9;
 119       mmod = *fITSGeom->GetMatrixSens(fActiveID,ild,idt);
 120       for (int ix=0;ix<2;ix++) {
 121         loc[0] = (ix-0.5)*kSegm->Dx(); // +-DX/2
 122         for (int iy=0;iy<2;iy++) {
 123           loc[1] = (iy-0.5)*kSegm->Dy(); // +-DY/2
 124           for (int iz=0;iz<2;iz++) {
 125             loc[2] = (iz-0.5)*kSegm->Dz(); // +-DZ/2
 126             //
 127             mmod.LocalToMaster(loc,glo);
 128             double phi = ATan2(glo[1],glo[0]);
 129             double r   = glo[0]*glo[0] + glo[1]*glo[1];
 130             if (fRMin>r) fRMin = r;
 131             if (fRMax<r) fRMax = r;
 132             BringTo02Pi(phi);
 133             if      (phiMin>1e8) phiMin=phi;
 134             else if (!OKforPhiMin(phiMin,phi)) phiMin=phi;
 135             if      (phiMax<-1e8) phiMax=phi;
 136             else if (!OKforPhiMax(phiMax,phi)) phiMax=phi;
 137             if (glo[2]>zMax) zMax=glo[2];
 138             if (glo[2]<zMin) zMin=glo[2];
 139           }
 140         }
 141       }
 142       sens->SetBoundaries(phiMin,phiMax,zMin,zMax);
 143       mt2l = fITSGeom->GetMatrixT2L(fActiveID,ild,idt);
 144       mmod.Multiply(mt2l);
 145       loc[0]=loc[1]=loc[2]=0;
 146       mmod.LocalToMaster(loc,glo);
 147       rTF   = Sqrt(glo[0]*glo[0] + glo[1]*glo[1]);  //  tracking params (misaligned)
 148       phiTF = ATan2(glo[1],glo[0]);
 149       BringTo02Pi(phiTF);
 150       //
 151       sens->SetXTF(rTF);
 152       sens->SetPhiTF(phiTF);
 153       //
 154       if      (fPhiLadMin[ild]>1e8)  fPhiLadMin[ild] = phiMin;
 155       else if (!OKforPhiMin(fPhiLadMin[ild],phiMin)) fPhiLadMin[ild] = phiMin;
 156       if      (fPhiLadMax[ild]<-1e8) fPhiLadMax[ild] = phiMax;
 157       else if (!OKforPhiMax(fPhiLadMax[ild],phiMax)) fPhiLadMax[ild] = phiMax;
 158       if (fZMin>zMin) fZMin = zMin;
 159       if (fZMax<zMax) fZMax = zMax;
 160       //
 161       if (idt>0) fSensDZInv += zMax - GetSensor(ild,idt-1)->GetZMax(); // z interval to previous
 162     }
 163   }
 164   //
 165   fRMin = Sqrt(fRMin);
 166   fRMax = Sqrt(fRMax);
 167   fR = 0.5*(fRMin+fRMax);
 168   fRMin -= kSafeR;
 169   fRMax += kSafeR;
 170   double dz = fNSensInLadder>0 ? fSensDZInv/(fNSensInLadder-1)/fNLadders : fZMax-fZMin;
 171   fSensDZInv = 1./dz;
 172
 173   const int kNBId[3][3] = {
 174     {AliITSURecoSens::kNghbBL,AliITSURecoSens::kNghbB,AliITSURecoSens::kNghbBR},
 175     {AliITSURecoSens::kNghbL,          -1            ,AliITSURecoSens::kNghbR },
 176     {AliITSURecoSens::kNghbTL,AliITSURecoSens::kNghbT,AliITSURecoSens::kNghbTR}
 177   };
 178
 179   // add neighbours info
 180   double zTol = 0.45*dz, phiTol = 0.45*TwoPi()/fNLadders;
 181   for (int ild=0;ild<fNLadders;ild++) {
 182     for (int idt=0;idt<fNSensInLadder;idt++) {
 183       AliITSURecoSens* sens = GetSensor(ild,idt);
 184       //
 185       for (int ils=-1;ils<=1;ils++) {
 186         int ildN = ild+ils;  // ladders of neighbouring sensors
 187         if (ildN<0) ildN = fNLadders-1; else if (ildN==fNLadders) ildN = 0;
 188         for (int ids=-1;ids<=1;ids++) {
 189           int idtN = idt+ids;
 190           if (idtN<0 || idtN==fNSensInLadder || (ids==0&&ils==0)) continue;
 191           AliITSURecoSens* sensN = GetSensor(ildN,idtN); // potential neighbor
 192           int neighbID = ildN*fNSensInLadder+idtN;
 193           //
 194           int zType = 1;  // side
 195           if (sens->GetZMin()-zTol  > sensN->GetZMax()) continue; // too large distance
 196           if (sensN->GetZMin()-zTol > sens->GetZMax() ) continue; // too large distance
 197           if      (sens->GetZMin()-zTol>sensN->GetZMin()) zType =  0;     // bottom
 198           else if (sensN->GetZMin()-zTol>sens->GetZMin()) zType =  2;     // top
 199           //
 200           int phiType = 1;
 201
 202           double phiTstMn = sensN->GetPhiMin()-phiTol;
 203           BringTo02Pi(phiTstMn);
 204           if (!OKforPhiMax(sens->GetPhiMax(),phiTstMn)) continue; // too large angle
 205           double phiTstMx = sensN->GetPhiMax()+phiTol;
 206           BringTo02Pi(phiTstMx);
 207           if (!OKforPhiMin(sens->GetPhiMin(),phiTstMx)) continue; // too large angle
 208           //
 209           phiTstMn = sensN->GetPhiMin()+phiTol;
 210           BringTo02Pi(phiTstMn);
 211           phiTstMx = sensN->GetPhiMax()-phiTol;
 212           BringTo02Pi(phiTstMx);
 213           if      (!OKforPhiMax(sens->GetPhiMax(),phiTstMx)) phiType = 0; // left
 214           else if (!OKforPhiMin(sens->GetPhiMin(),phiTstMn)) phiType = 2; // right
 215           //
 216           sens->SetNeighborID(kNBId[zType][phiType], neighbID);
 217         } // phi scan
 218       } // z scan
 219     } // sensors
 220   } // ladders
 221   //
 222 }
 223 ///*
 224 //______________________________________________________
 225 Int_t AliITSURecoLayer::FindSensors(const double* impPar, AliITSURecoSens *sensors[AliITSURecoSens::kNNeighbors])
 226 {
 227   // find sensors having intersection with track
 228   // impPar contains: lab phi of track, dphi, labZ, dz
 229   //
 230   double z = impPar[2];
 231   if (z>fZMax+impPar[3]) return 0; // outside of Z coverage
 232   z -= fZMin;
 233   if (z<-impPar[3]) return 0; // outside of Z coverage
 234   int sensInLad = int(z*fSensDZInv);
 235   if      (sensInLad<0) sensInLad = 0;
 236   else if (sensInLad>=fNSensInLadder) sensInLad = fNSensInLadder-1;
 237   //
 238   double phi = impPar[0] - fPhiOffs;
 239   BringTo02Pi(phi);
 240   int ladID = int(phi*fDPhiLadInv);  // ladder id
 241   int nsens = 0;
 242   //
 243   AliITSURecoSens* sensN,*sens = GetSensor(ladID*fNSensInLadder+sensInLad);
 244   sensors[nsens++] = sens;
 245   //
 246   // check neighbours
 247   double zMn=impPar[2]-impPar[3], zMx=impPar[2]+impPar[3], phiMn=impPar[0]-impPar[1], phiMx=impPar[0]+impPar[1];
 248   BringTo02Pi(phiMn);
 249   BringTo02Pi(phiMx);
 250   //
 251   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbR)); // neighbor on the right (smaller phi)
 252   if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax())) sensors[nsens++] = sensN;
 253   //
 254   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTR)); // neighbor on the top right (smaller phi, larger Z)
 255   if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
 256   //
 257   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbT)); // neighbor on the top (larger Z)
 258   if (sensN && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
 259   //
 260   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTL)); // neighbor on the top left (larger Z, larger phi)
 261   if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
 262   //
 263   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbL)); // neighbor on the left (larger phi)
 264   if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin())) sensors[nsens++] = sensN;
 265   //
 266   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBL)); // neighbor on the bottom left (smaller Z, larger phi)
 267   if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
 268   //
 269   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbB));  // neighbor on the bottom (smaller Z)
 270   if (sensN && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
 271   //
 272   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBR)); // neighbor on the bottom right (smaller Z, smaller phi)
 273   if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
 274   //
 275   return nsens;
 276 }
 277 //*/
 278 /*
 279 Int_t AliITSURecoLayer::FindSensors(const double* impPar, AliITSURecoSens *sensors[AliITSURecoSens::kNNeighbors],int mcLab)
 280 {
 281   // find sensors having intersection with track
 282   // impPar contains: lab phi of track, dphi, labZ, dz
 283
 284   //tmp>>>
 285   int nFnd = 0;
 286   int fndSens[50];
 287   if (mcLab>=0) { // find correct sensors from MC info
 288     int ncl = GetNClusters();
 289     for (int icl=ncl;icl--;) {
 290       AliCluster* cl = GetCluster(icl);
 291       for (int ilb=0;ilb<3;ilb++) {
 292         if (cl->GetLabel(ilb)<0) break;
 293         if (cl->GetLabel(ilb)==mcLab) {fndSens[nFnd++] = cl->GetVolumeId(); break;}
 294       }
 295     }
 296     if (nFnd>0) {
 297       Int_t layS,ladS,sensS;
 298       for (int is=0;is<nFnd;is++) {
 299         fITSGeom->GetModuleId(fndSens[is],layS,ladS,sensS);
 300         printf("SNMC#%d(%d): %d %d %d | ",is,mcLab,layS,ladS,sensS); GetSensorFromID(fndSens[is])->Print();
 301       }
 302     }
 303   }
 304
 305   //tmp<<<
 306
 307   double z = impPar[2];
 308   if (z>fZMax+impPar[3]) {
 309     if (nFnd>0) printf("MissedSens!!!\n");
 310     return 0; // outside of Z coverage
 311   }
 312   z -= fZMin;
 313   if (z<-impPar[3]) {
 314     if (nFnd>0) printf("MissedSens!!!\n");
 315     return 0; // outside of Z coverage
 316   }
 317   int sensInLad = int(z*fSensDZInv);
 318   if      (sensInLad<0) sensInLad = 0;
 319   else if (sensInLad>=fNSensInLadder) sensInLad = fNSensInLadder-1;
 320   //
 321   double phi = impPar[0] - fPhiOffs;
 322   BringTo02Pi(phi);
 323   int ladID = int(phi*fDPhiLadInv);  // ladder id
 324   int nsens = 0;
 325   //
 326   AliITSURecoSens* sensN,*sens = GetSensor(ladID*fNSensInLadder+sensInLad);
 327   sensors[nsens++] = sens;
 328   //
 329   // check neighbours
 330   double zMn=impPar[2]-impPar[3], zMx=impPar[2]+impPar[3], phiMn=impPar[0]-impPar[1], phiMx=impPar[0]+impPar[1];
 331   BringTo02Pi(phiMn);
 332   BringTo02Pi(phiMx);
 333   //
 334   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbR)); // neighbor on the right (smaller phi)
 335   if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax())) sensors[nsens++] = sensN;
 336   //
 337   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTR)); // neighbor on the top right (smaller phi, larger Z)
 338   if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
 339   //
 340   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbT)); // neighbor on the top (larger Z)
 341   if (sensN && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
 342   //
 343   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTL)); // neighbor on the top left (larger Z, larger phi)
 344   if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
 345   //
 346   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbL)); // neighbor on the left (larger phi)
 347   if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin())) sensors[nsens++] = sensN;
 348   //
 349   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBL)); // neighbor on the bottom left (smaller Z, larger phi)
 350   if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
 351   //
 352   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbB));  // neighbor on the bottom (smaller Z)
 353   if (sensN && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
 354   //
 355   sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBR)); // neighbor on the bottom right (smaller Z, smaller phi)
 356   if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
 357   //
 358   if (mcLab>=0) {
 359     Int_t layS,ladS,sensS;
 360     printf("Found %d sensors for phi %.3f : %.3f | Z %.4f %.4f\n", nsens,phiMn,phiMx,zMn,zMx);
 361     for (int is=0;is<nsens;is++) {
 362       fITSGeom->GetModuleId(sensors[is]->GetID()+fITSGeom->GetFirstModIndex(fActiveID),layS,ladS,sensS);
 363       printf("*SNF#%d: %d %d %d | ",is,layS,ladS,sensS); sensors[is]->Print();
 364     }
 365     for (int ism=0;ism<nFnd;ism++) {
 366       AliITSURecoSens* snMC = GetSensorFromID(fndSens[ism]);
 367       Bool_t ok=kFALSE;
 368       for (int isf=0;isf<nsens;isf++) {
 369         if (snMC==sensors[isf]) {ok=kTRUE;break;}
 370       }
 371       if (!ok) printf("MissedSens %d!!!\n",ism);
 372     }
 373   }
 374   return nsens;
 375 }
 376 */
 377 //______________________________________________________
 378 void AliITSURecoLayer::ProcessClusters(Int_t mode)
 379 {
 380   // register in each sensor of the layer its cluster.
 381   // the clusters of the layer must be sorted per sensor
 382   int ncl = fClusters->GetEntriesFast();
 383   int curSensID = -1;
 384   for (int i=fNSensors;i--;) GetSensor(i)->SetNClusters(0);
 385   AliITSURecoSens* curSens = 0;
 386   for (int icl=0;icl<ncl;icl++) {
 387     AliITSUClusterPix* cl = (AliITSUClusterPix*) fClusters->UncheckedAt(icl);
 388     cl->GoToFrameTrk();
 389     int vID = cl->GetVolumeId();
 390     if (vID<curSensID) {AliFatal("Clusters are not sorted in increasing sensorID");}
 391     if (vID>curSensID) {
 392       if (curSens) curSens->ProcessClusters(mode);    // prepare clusters for reconstruction
 393       curSens   = GetSensor(vID - fITSGeom->GetFirstModIndex(fActiveID));
 394       curSensID = vID;
 395       curSens->SetFirstClusterId(icl);
 396     }
 397     curSens->IncNClusters();
 398   }
 399   if (curSens) curSens->ProcessClusters(mode); // last sensor was not processed yet
 400   //
 401 }
 402
 403 //______________________________________________________
 404 Bool_t AliITSURecoLayer::IsEqual(const TObject* obj) const
 405 {
 406   // check if layers are equal in R
 407   const AliITSURecoLayer* lr = (const AliITSURecoLayer*)obj;
 408   return Abs(lr->GetR()-GetR())<1e-6 ? kTRUE : kFALSE;
 409 }
 410
 411 //______________________________________________________
 412 Int_t  AliITSURecoLayer::Compare(const TObject* obj) const
 413 {
 414   // compare two layers
 415   const AliITSURecoLayer* lr = (const AliITSURecoLayer*)obj;
 416   double dr = GetR() - lr->GetR();
 417   if (Abs(dr)<1e-6) return 0;
 418   return dr>0 ? 1:-1;
 419   //
 420 }
 421
 422 //_________________________________________________________________
 423 AliITSURecoSens* AliITSURecoLayer::GetSensorFromID(Int_t i) const
 424 {
 425   // get sensor from its global id
 426   i -= fITSGeom->GetFirstModIndex(fActiveID);
 427   if (i<0||i>=fNSensors) AliFatal(Form("Sensor with id=%d is not in layer %d",i+fITSGeom->GetFirstModIndex(fActiveID),fActiveID));
 428   return (AliITSURecoSens*)fSensors[i];
 429 }