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