Fix in layers rmin/rmax calculation and their accounting in track navigation.
[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSURecoLayer.cxx
index 7d33104..7b2e12e 100644 (file)
@@ -1,9 +1,10 @@
 #include <TClonesArray.h>
 #include "AliITSURecoLayer.h"
-#include "AliITSUGeomTGeo.h"
 #include "AliITSsegmentation.h"
 #include "AliITSUAux.h"
 #include "AliITSUClusterPix.h"
+#include "AliITSUGeomTGeo.h"
+#include "AliLog.h"
 
 using namespace AliITSUAux;
 using namespace TMath;
@@ -69,13 +70,15 @@ AliITSURecoLayer::~AliITSURecoLayer()
   delete[] fSensors;
   delete[] fPhiLadMax;
   delete[] fPhiLadMin;
+  if (GetOwnsClusterArray()) delete fClusters;
 }
 
 //______________________________________________________
 void AliITSURecoLayer::Print(Option_t* opt) const                            
 {
   //print 
-  printf("Layer %s %d (active:%+d), NSensors: %d\n",GetName(),GetID(),GetActiveID(),GetNSensors());
+  printf("Lr %-15s %d (act:%+d), NSens: %4d | MaxStep:%.2f ",GetName(),GetID(),GetActiveID(),GetNSensors(),fMaxStep);
+  printf("%6.3f<R<%6.3f | %+8.3f<Z<%+8.3f dZ:%6.3f\n",fRMin,fRMax,fZMin,fZMax,fSensDZInv>0 ? 1/fSensDZInv : 0);
   TString opts = opt; opts.ToLower();
   if (opts.Contains("sn")) for (int i=0;i<GetNSensors();i++) GetSensor(i)->Print(opt);
 }
@@ -84,6 +87,7 @@ void AliITSURecoLayer::Print(Option_t* opt) const
 void AliITSURecoLayer::Build()
 {
   // build internal structures
+  const double kSafeR = 0.05; // safety margin for Rmin,Rmax of the layer
   if (fActiveID<0) return;
   fNLadders = fITSGeom->GetNLadders(fActiveID);
   fNSensInLadder = fITSGeom->GetNDetectors(fActiveID);
@@ -101,6 +105,7 @@ void AliITSURecoLayer::Build()
   fPhiLadMin = new Double_t[fNLadders];
   fPhiLadMax = new Double_t[fNLadders];
   fSensDZInv = 0;
+  fDPhiLadInv = fNLadders/TwoPi();
   //
   for (int ild=0;ild<fNLadders;ild++) {
     fPhiLadMin[ild] = 1e9;
@@ -108,10 +113,10 @@ void AliITSURecoLayer::Build()
     //
     for (int idt=0;idt<fNSensInLadder;idt++) {
       AliITSURecoSens* sens = new AliITSURecoSens(fNSensors++);
-      fSensors[idt] = sens;
+      fSensors[ild*fNSensInLadder+idt] = sens;
       //
       double phiMin=1e9,phiMax=-1e9,zMin=1e9,zMax=-1e9;
-      mmod = *fITSGeom->GetMatrix(fActiveID,ild,idt);
+      mmod = *fITSGeom->GetMatrixSens(fActiveID,ild,idt);
       for (int ix=0;ix<2;ix++) {
        loc[0] = (ix-0.5)*kSegm->Dx(); // +-DX/2
        for (int iy=0;iy<2;iy++) {
@@ -129,12 +134,11 @@ void AliITSURecoLayer::Build()
            else if (!OKforPhiMin(phiMin,phi)) phiMin=phi;
            if      (phiMax<-1e8) phiMax=phi;
            else if (!OKforPhiMax(phiMax,phi)) phiMax=phi;            
-           if (glo[2]>zMax) zMax=glo[2]; else if (glo[2]<zMin) zMin=glo[2];
+           if (glo[2]>zMax) zMax=glo[2];
+           if (glo[2]<zMin) zMin=glo[2];
          }
        }
       }
-      printf("%d %d (%d)  %f %f\n",ild,idt, ild*fNSensInLadder+idt,phiMin,phiMax);
-
       sens->SetBoundaries(phiMin,phiMax,zMin,zMax);
       mt2l = fITSGeom->GetMatrixT2L(fActiveID,ild,idt);
       mmod.Multiply(mt2l);     
@@ -154,13 +158,15 @@ void AliITSURecoLayer::Build()
       if (fZMin>zMin) fZMin = zMin;
       if (fZMax<zMax) fZMax = zMax;
       //
-      if (idt>0) fSensDZInv += zMax - GetSensor(fNSensors-2)->GetZMax(); // z interval to previoud
+      if (idt>0) fSensDZInv += zMax - GetSensor(ild,idt-1)->GetZMax(); // z interval to previous
     }
   }
   //
   fRMin = Sqrt(fRMin);
   fRMax = Sqrt(fRMax);
   fR = 0.5*(fRMin+fRMax);
+  fRMin -= kSafeR;
+  fRMax += kSafeR;
   double dz = fNSensInLadder>0 ? fSensDZInv/(fNSensInLadder-1)/fNLadders : fZMax-fZMin;
   fSensDZInv = 1./dz;
 
@@ -213,7 +219,6 @@ void AliITSURecoLayer::Build()
     } // sensors
   } // ladders
   //
-  //
 }
 
 //______________________________________________________
@@ -222,10 +227,12 @@ Int_t AliITSURecoLayer::FindSensors(const double* impPar, AliITSURecoSens *senso
   // find sensors having intersection with track
   // impPar contains: lab phi of track, dphi, labZ, dz
   double z = impPar[2];
-  if (z>fZMax+impPar[3]) return -1; // outside of Z coverage
+  if (z>fZMax+impPar[3]) return 0; // outside of Z coverage
   z -= fZMin;
-  if (z<-impPar[3]) return -1; // outside of Z coverage
+  if (z<-impPar[3]) return 0; // outside of Z coverage
   int sensInLad = int(z*fSensDZInv);
+  if      (sensInLad<0) sensInLad = 0;
+  else if (sensInLad>=fNSensInLadder) sensInLad = fNSensInLadder-1;
   //
   double phi = impPar[0] - fPhiOffs;
   BringTo02Pi(phi);
@@ -239,29 +246,29 @@ Int_t AliITSURecoLayer::FindSensors(const double* impPar, AliITSURecoSens *senso
   BringTo02Pi(phiMn);
   BringTo02Pi(phiMx);
   //
-  sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbR)); // neighbor on the right
-  if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin())) sensors[nsens++] = sensN;
+  sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbR)); // neighbor on the right (smaller phi)
+  if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax())) sensors[nsens++] = sensN;
   //
-  sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTR)); // neighbor on the top right
-  if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
+  sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTR)); // neighbor on the top right (smaller phi, larger Z)
+  if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
   //
-  sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbT)); // neighbor on the top
-  if (sensN && sensN->GetZMin()>zMx) sensors[nsens++] = sensN;
+  sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbT)); // neighbor on the top (larger Z)
+  if (sensN && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
   //
-  sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTL)); // neighbor on the top left
-  if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
+  sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTL)); // neighbor on the top left (larger Z, larger phi)
+  if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
   //
-  sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbL)); // neighbor on the left
-  if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax())) sensors[nsens++] = sensN;
+  sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbL)); // neighbor on the left (larger phi)
+  if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin())) sensors[nsens++] = sensN;
   //
-  sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBL)); // neighbor on the bottom left
-  if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
+  sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBL)); // neighbor on the bottom left (smaller Z, larger phi)
+  if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
   //
-  sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbB));  // neighbor on the bottom
+  sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbB));  // neighbor on the bottom (smaller Z)
   if (sensN && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
   //
-  sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBR)); // neighbor on the bottom right
-  if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
+  sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBR)); // neighbor on the bottom right (smaller Z, smaller phi)
+  if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
   //
   return nsens;
 }
@@ -273,6 +280,7 @@ void AliITSURecoLayer::ProcessClusters(Int_t mode)
   // the clusters of the layer must be sorted per sensor
   int ncl = fClusters->GetEntriesFast();
   int curSensID = -1;
+  for (int i=fNSensors;i--;) GetSensor(i)->SetNClusters(0);
   AliITSURecoSens* curSens = 0;
   for (int icl=0;icl<ncl;icl++) {
     AliITSUClusterPix* cl = (AliITSUClusterPix*) fClusters->UncheckedAt(icl);
@@ -290,3 +298,31 @@ void AliITSURecoLayer::ProcessClusters(Int_t mode)
   if (curSens) curSens->ProcessClusters(mode); // last sensor was not processed yet
   //
 }
+
+//______________________________________________________
+Bool_t AliITSURecoLayer::IsEqual(const TObject* obj) const
+{
+  // check if layers are equal in R
+  const AliITSURecoLayer* lr = (const AliITSURecoLayer*)obj;
+  return Abs(lr->GetR()-GetR())<1e-6 ? kTRUE : kFALSE;
+}
+
+//______________________________________________________
+Int_t  AliITSURecoLayer::Compare(const TObject* obj) const
+{
+  // compare two layers
+  const AliITSURecoLayer* lr = (const AliITSURecoLayer*)obj;
+  double dr = GetR() - lr->GetR();
+  if (Abs(dr)<1e-6) return 0;
+  return dr>0 ? 1:-1;
+  //      
+}
+
+//_________________________________________________________________
+AliITSURecoSens* AliITSURecoLayer::GetSensorFromID(Int_t i) const 
+{
+  // get sensor from its global id
+  i -= fITSGeom->GetFirstModIndex(fActiveID);
+  if (i<0||i>=fNSensors) AliFatal(Form("Sensor with id=%d is not in layer %d",i+fITSGeom->GetFirstModIndex(fActiveID),fActiveID));
+  return (AliITSURecoSens*)fSensors[i];
+}