[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSURecoLayer.cxx

#include <TClonesArray.h>
#include "AliITSURecoLayer.h"
#include "AliITSsegmentation.h"
#include "AliITSUAux.h"
#include "AliITSUClusterPix.h"
#include "AliITSUGeomTGeo.h"
#include "AliLog.h"

using namespace AliITSUAux;
using namespace TMath;

ClassImp(AliITSURecoLayer)


//______________________________________________________
AliITSURecoLayer::AliITSURecoLayer(const char* name)
 :fActiveID(-1)
  ,fNSensors(0)
  ,fNSensorRows(0)
  ,fNSensorsPerRow(0)
  ,fSensVIDtoMatrixID(0)
  ,fR(0)
  ,fRMax(0)
  ,fRMin(0)
  ,fZMax(0)
  ,fZMin(0)
  ,fPhiOffs(0)
  ,fSensDZInv(0)
  ,fSensDPhiInv(0)
  ,fMaxStep(0.5)
  ,fSensors(0)
  ,fITSGeom(0)
  ,fClusters(0)
{
  // def. c-tor
  SetNameTitle(name,name);
}

//______________________________________________________
AliITSURecoLayer::AliITSURecoLayer(const char* name, Int_t activeID, AliITSUGeomTGeo* gm)
  :fActiveID(activeID)
  ,fNSensors(0)
  ,fNSensorRows(0)
  ,fNSensorsPerRow(0)
  ,fSensVIDtoMatrixID(0)
  ,fR(0)
  ,fRMax(0)
  ,fRMin(0)
  ,fZMax(0)
  ,fZMin(0)
  ,fPhiOffs(0)
  ,fSensDZInv(0)
  ,fSensDPhiInv(0)
  ,fMaxStep(0.5)
  ,fSensors(0)
  ,fITSGeom(gm)
  ,fClusters(0)
{
  // def. c-tor
  SetNameTitle(name,name);
  Build();
}

//______________________________________________________
AliITSURecoLayer::~AliITSURecoLayer()
{
  // def. d-tor
  delete fSensors;
  delete[] fSensVIDtoMatrixID;
  if (GetOwnsClusterArray()) delete fClusters;
}

//______________________________________________________
void AliITSURecoLayer::Print(Option_t* opt) const			      
{
  //print 
  printf("Lr %-15s %d (act:%+d), NSens: %4d in %3d rows| MaxStep:%.2f ",
	 GetName(),GetID(),GetActiveID(),GetNSensors(),GetNSensorRows(),fMaxStep);
  printf("%6.3f<R<%6.3f | %+8.3f<Z<%+8.3f dZ:%6.3f dPhi:%6.3f\n",fRMin,fRMax,fZMin,fZMax,
	 fSensDZInv>0 ? 1/fSensDZInv : 0, fSensDPhiInv>0 ? 1/fSensDPhiInv : 0);
  TString opts = opt; opts.ToLower();
  if (opts.Contains("sn")) for (int i=0;i<GetNSensors();i++) GetSensor(i)->Print(opt);
}

//______________________________________________________
void AliITSURecoLayer::Build()
{
  // build internal structures
  const double kSafeR = 0.05; // safety margin for Rmin,Rmax of the layer
  if (fActiveID<0) return;
  //
  int nStaves=fITSGeom->GetNStaves(fActiveID);
  // determine number of sensor rows (sensors aligned at same phi and spanning the Z range of the layer)
  fNSensorRows = nStaves;
  //
  // if the stave has susbtaves, each substave can have multiple rows of sensors (but just 1 row of modules)
  if (fITSGeom->GetNHalfStaves(fActiveID)>0) fNSensorRows *= fITSGeom->GetNHalfStaves(fActiveID);
  //
  // if there are modules defined, the module may have multiple rows of sensors (though not spanning full Z)
  if (fITSGeom->GetNModules(fActiveID)>0) fNSensorRows *= fITSGeom->GetNChipRowsPerModule(fActiveID);
  //
  fNSensors = fITSGeom->GetNChipsPerLayer(fActiveID);
  fNSensorsPerRow = fNSensors/fNSensorRows;
  //
  fSensors = new TObjArray(fNSensors);
  fSensVIDtoMatrixID = new Int_t[fNSensors];
  const AliITSsegmentation* kSegm = fITSGeom->GetSegmentation(fActiveID);
  //
  TGeoHMatrix mmod;
  const TGeoHMatrix* mt2l;
  double phiTF,rTF, loc[3]={0,0,0},glo[3];
  //
  int nSensPerStave = fITSGeom->GetNChipsPerStave(fActiveID);
  for (int staveI=0;staveI<nStaves;staveI++) {
    for (int sensI=0;sensI<nSensPerStave;sensI++) {
      int sID = fITSGeom->GetChipIndex(fActiveID,staveI,sensI);
      AliITSURecoSens* sens = new AliITSURecoSens( sID );
      fSensors->AddLast(sens);
      double phiMin=1e9,phiMax=-1e9,zMin=1e9,zMax=-1e9;
      // this is NOT the sensor matrix, just the ideal chip matrix to get neighbors correct
      fITSGeom->GetOrigMatrix(sID,mmod); 
      //
      for (int ix=0;ix<2;ix++) {       // determine sensor boundaries (ideal)
	loc[0] = (ix-0.5)*kSegm->Dx(); // +-DX/2
	for (int iy=0;iy<2;iy++) {
	  loc[1] = (iy-0.5)*kSegm->Dy(); // +-DY/2
	  for (int iz=0;iz<2;iz++) {
	    loc[2] = (iz-0.5)*kSegm->Dz(); // +-DZ/2
	    mmod.LocalToMaster(loc,glo);
	    double phi = ATan2(glo[1],glo[0]);
	    BringTo02Pi(phi);
	    if      (phiMin>1e8)  phiMin=phi;
	    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];
	    if (glo[2]<zMin) zMin=glo[2];
	  }
	}
      }
      sens->SetBoundaries(phiMin,phiMax,zMin,zMax);
    }
  }
  fSensors->Sort(); // sort sensors to get the neighborhood correct
  //
  // now fill real sensor angles, Z's, accounting for misalignment
  fRMin=fZMin=1e9;
  fRMax=fZMax=-1e9;
  //
  fPhiOffs = 0;
  int firstSensID = fITSGeom->GetFirstChipIndex(fActiveID);
  for (int sensI=0;sensI<fNSensors;sensI++) {
    AliITSURecoSens* sens = GetSensor(sensI);
    mmod = *fITSGeom->GetMatrixSens(sens->GetID());
    fSensVIDtoMatrixID[sens->GetID() - firstSensID] = sensI;
    double phiMin=1e9,phiMax=-1e9,zMin=1e9,zMax=-1e9;
    for (int ix=0;ix<2;ix++) {
      loc[0] = (ix-0.5)*kSegm->Dx(); // +-DX/2
      for (int iy=0;iy<2;iy++) {
	loc[1] = (iy-0.5)*kSegm->Dy(); // +-DY/2
	for (int iz=0;iz<2;iz++) {
	  loc[2] = (iz-0.5)*kSegm->Dz(); // +-DZ/2
	  //
	  mmod.LocalToMaster(loc,glo);
	  double phi = ATan2(glo[1],glo[0]);
	  double r   = glo[0]*glo[0] + glo[1]*glo[1];
	  if (fRMin>r) fRMin = r;
	  if (fRMax<r) fRMax = r;
	  BringTo02Pi(phi);
	  if      (phiMin>1e8) phiMin=phi; 
	  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];
	  if (glo[2]<zMin) zMin=glo[2];
	}
      }
    }
    mt2l = fITSGeom->GetMatrixT2L( sens->GetID() );
    mmod.Multiply(mt2l);	
    loc[0]=loc[1]=loc[2]=0;
    mmod.LocalToMaster(loc,glo);
    rTF   = Sqrt(glo[0]*glo[0] + glo[1]*glo[1]);  //  tracking params (misaligned)
    phiTF = ATan2(glo[1],glo[0]);
    BringTo02Pi(phiTF);
    //
    sens->SetXTF(rTF);
    sens->SetPhiTF(phiTF);
    sens->SetBoundaries(phiMin,phiMax,zMin,zMax);
    if (fZMin>zMin) fZMin = zMin;
    if (fZMax<zMax) fZMax = zMax;
    //
    if (sensI<fNSensorsPerRow) fPhiOffs += MeanPhiSmall(phiMax,phiMin);
  }
  //
  fPhiOffs /= fNSensorsPerRow; // average phi of the 1st row
  fSensDZInv = fNSensorsPerRow/(fZMax-fZMin);
  fSensDPhiInv = fNSensorRows/(2*Pi());
  //
  fRMin = Sqrt(fRMin);
  fRMax = Sqrt(fRMax);
  fR = 0.5*(fRMin+fRMax);
  fRMin -= kSafeR;
  fRMax += kSafeR;
  //
}

//______________________________________________________
Int_t AliITSURecoLayer::FindSensors(const double* impPar, AliITSURecoSens *sensors[kMaxSensMatching])
{
  // find sensors having intersection with track
  // impPar contains: lab phi of track, dphi, labZ, dz
  //
  double zMn=impPar[2]-impPar[3], zMx=impPar[2]+impPar[3]; 
  if (zMn>fZMax) return 0;
  if (zMx<fZMin) return 0;
  //
  int zCenID = int((impPar[2]-fZMin)*fSensDZInv);
  if      (zCenID<0) zCenID = 0;
  else if (zCenID>=fNSensorsPerRow) zCenID = fNSensorsPerRow-1;
  double phiCn = impPar[0] - fPhiOffs;
  //  BringTo02Pi(phiCn); 
  int rowCenID = int(phiCn*fSensDPhiInv);
  //
  // due to the misalignments the actual sensorID's might be shifted
  int res = 0;
  AliITSURecoSens* sensPrev=0, *sens = GetSensor(rowCenID,zCenID);
  //
  //  printf("Guess: Primary Sensor: phiID: %d zID: %d ->",rowCenID,zCenID); sens->Print();
  //
  while ( (res=sens->CheckCoverage(impPar[0], impPar[2])) ) {
    if      (res&AliITSURecoSens::kRight) {if (++rowCenID==fNSensorRows) rowCenID=0;} // neighbor on the right (larger phi)
    else if (res&AliITSURecoSens::kLeft)  {if (--rowCenID<0) rowCenID = fNSensorRows-1;}      // neighbor on the left (smaller phi)
    if      (res&AliITSURecoSens::kUp)    {if (++zCenID==fNSensorsPerRow) zCenID = fNSensorsPerRow-1;}    // neighbor at larger Z (if any)
    else if (res&AliITSURecoSens::kDown)  {if (--zCenID<0) zCenID = 0;}               // neighbor at smaller Z (if any)
    //
    AliITSURecoSens* sensAlt = GetSensor(rowCenID,zCenID);
    if (sensAlt==sens || sensAlt==sensPrev) break;  // there is no better neighbor (z edge) or the point falls in dead area
    //
    sensPrev = sens;
    sens = sensAlt;
  }
  //  printf("Found: Primary Sensor: phiID: %d zID: %d ->",rowCenID,zCenID); sens->Print();
  //
  int nFnd = 0;
  sensors[nFnd++] = sens;
  //
  double phiMn = impPar[0]-impPar[1], phiMx = impPar[0]+impPar[1];
  BringTo02Pi(phiMn);
  BringTo02Pi(phiMx);
  //
  const int kNNeighb = 8;
  const int kCheckNeighb[2][kNNeighb] = { // phi and Z neighbours to check
    { 1, 1, 0,-1,-1,-1, 0, 1},
    { 0, 1, 1, 1, 0,-1,-1,-1}
  };
  //  
  //  printf("Search: %+.4f %+.4f | %+.4f %+.4f\n",phiMn,phiMx, zMn,zMx);
  for (int inb=kNNeighb;inb--;) {
    int idz = kCheckNeighb[1][inb];
    int iz   = zCenID   + idz;
    //    printf("#%d  dp:%+d dz:%+d IZ: %d\n",inb, kCheckNeighb[0][inb], kCheckNeighb[1][inb], iz);
    //
    if (iz<0 || iz>=fNSensorsPerRow) continue;
    int idphi = kCheckNeighb[0][inb];
    int iphi = rowCenID + idphi;
    if      (iphi<0) iphi += fNSensorRows;
    else if (iphi>=fNSensorRows) iphi -= fNSensorRows;
    sens = GetSensor(iphi,iz);
    //
    if      (idz>0) {if (zMx<sens->GetZMin()) continue;}
    else if (idz<0) {if (zMn>sens->GetZMax()) continue;}
    //
    // Z range matches
    if      (idphi>0) {if (!OKforPhiMin(sens->GetPhiMin(),phiMx)) continue;}
    else if (idphi<0) {if (!OKforPhiMax(sens->GetPhiMax(),phiMn)) continue;}
    //
    //    printf("Add %d\n",nFnd);
    sensors[nFnd++] = sens;
    if (nFnd==kMaxSensMatching) break;
  }
  return nFnd;
}

//*/
/*
Int_t AliITSURecoLayer::FindSensors(const double* impPar, AliITSURecoSens *sensors[AliITSURecoSens::kNNeighbors],int mcLab)
{
  // find sensors having intersection with track
  // impPar contains: lab phi of track, dphi, labZ, dz

  //tmp>>>
  int nFnd = 0;
  int fndSens[50];
  if (mcLab>=0) { // find correct sensors from MC info
    int ncl = GetNClusters();
    for (int icl=ncl;icl--;) {
      AliCluster* cl = GetCluster(icl);
      for (int ilb=0;ilb<3;ilb++) {
	if (cl->GetLabel(ilb)<0) break;
	if (cl->GetLabel(ilb)==mcLab) {fndSens[nFnd++] = cl->GetVolumeId(); break;}
      }
    }
    if (nFnd>0) {
      Int_t layS,staS,sensS;
      for (int is=0;is<nFnd;is++) {
	fITSGeom->GetChipId(fndSens[is],layS,staS,sensS);
	printf("SNMC#%d(%d): %d %d %d | ",is,mcLab,layS,staS,sensS); GetSensorFromID(fndSens[is])->Print();
      }
    }
  }

  //tmp<<<

  double z = impPar[2];
  if (z>fZMax+impPar[3]) {
    if (nFnd>0) printf("MissedSens!!!\n");
    return 0; // outside of Z coverage
  }
  z -= fZMin;
  if (z<-impPar[3]) {
    if (nFnd>0) printf("MissedSens!!!\n");
    return 0; // outside of Z coverage
  }
  int sensInSta = int(z*fSensDZInv);
  if      (sensInSta<0) sensInSta = 0;
  else if (sensInSta>=fNSensInStave) sensInSta = fNSensInStave-1;
  //
  double phi = impPar[0] - fPhiOffs;
  BringTo02Pi(phi);
  int staID = int(phi*fSensDPhiInv);  // stave id
  int nsens = 0;
  //
  AliITSURecoSens* sensN,*sens = GetSensor(staID*fNSensInStave+sensInSta);
  sensors[nsens++] = sens;
  //
  // check neighbours
  double zMn=impPar[2]-impPar[3], zMx=impPar[2]+impPar[3], phiMn=impPar[0]-impPar[1], phiMx=impPar[0]+impPar[1];
  BringTo02Pi(phiMn);
  BringTo02Pi(phiMx);
  //
  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 (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 (larger Z)
  if (sensN && 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 (larger phi)
  if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin())) 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 (smaller Z)
  if (sensN && 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;
  //
  if (mcLab>=0) {
    Int_t layS,staS,sensS;
    printf("Found %d sensors for phi %.3f : %.3f | Z %.4f %.4f\n", nsens,phiMn,phiMx,zMn,zMx); 
    for (int is=0;is<nsens;is++) {
      fITSGeom->GetChipId(sensors[is]->GetID()+fITSGeom->GetFirstModIndex(fActiveID),layS,staS,sensS);
      printf("*SNF#%d: %d %d %d | ",is,layS,staS,sensS); sensors[is]->Print();
    }
    for (int ism=0;ism<nFnd;ism++) {
      AliITSURecoSens* snMC = GetSensorFromID(fndSens[ism]);
      Bool_t ok=kFALSE;
      for (int isf=0;isf<nsens;isf++) {
	if (snMC==sensors[isf]) {ok=kTRUE;break;}
      }
      if (!ok) printf("MissedSens %d!!!\n",ism);
    }
  }
  return nsens;
}
*/
//______________________________________________________
void AliITSURecoLayer::ProcessClusters(Int_t mode)
{
  // register in each sensor of the layer its cluster.
  // 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);
    cl->GoToFrameTrk();
    int vID = cl->GetVolumeId();
    if (vID<curSensID) {AliFatal("Clusters are not sorted in increasing sensorID");}
    if (vID>curSensID) {
      if (curSens) curSens->ProcessClusters(mode);    // prepare clusters for reconstruction
      curSens   = GetSensorFromID(vID); //GetSensor(vID - fITSGeom->GetFirstChipIndex(fActiveID));
      curSensID = vID;
      curSens->SetFirstClusterId(icl);
    }
    curSens->IncNClusters();
  }
  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->GetFirstChipIndex(fActiveID);
  if (i<0||i>=fNSensors) AliFatal(Form("Sensor with id=%d is not in layer %d",i+fITSGeom->GetFirstChipIndex(fActiveID),fActiveID));
  return GetSensor(SensVIDtoMatrixID(i));
}
Commit	Line	Data
32d38de2	1	#include <TClonesArray.h>
a11ef2e4	2	#include "AliITSURecoLayer.h"
a11ef2e4	3	#include "AliITSsegmentation.h"
a11ef2e4	4	#include "AliITSUAux.h"
32d38de2	5	#include "AliITSUClusterPix.h"
8b16dbae	6	#include "AliITSUGeomTGeo.h"
8b16dbae	7	#include "AliLog.h"
a11ef2e4	8
	9	using namespace AliITSUAux;
	10	using namespace TMath;
	11
	12	ClassImp(AliITSURecoLayer)
	13
	14
	15	//______________________________________________________
	16	AliITSURecoLayer::AliITSURecoLayer(const char* name)
68b7631d	17	:fActiveID(-1)
a11ef2e4	18	,fNSensors(0)
68b7631d	19	,fNSensorRows(0)
	20	,fNSensorsPerRow(0)
	21	,fSensVIDtoMatrixID(0)
a11ef2e4	22	,fR(0)
	23	,fRMax(0)
	24	,fRMin(0)
	25	,fZMax(0)
	26	,fZMin(0)
a11ef2e4	27	,fPhiOffs(0)
a11ef2e4	28	,fSensDZInv(0)
68b7631d	29	,fSensDPhiInv(0)
32d38de2	30	,fMaxStep(0.5)
a11ef2e4	31	,fSensors(0)
a11ef2e4	32	,fITSGeom(0)
32d38de2	33	,fClusters(0)
a11ef2e4	34	{
	35	// def. c-tor
	36	SetNameTitle(name,name);
a11ef2e4	37	}
	38
	39	//______________________________________________________
32d38de2	40	AliITSURecoLayer::AliITSURecoLayer(const char* name, Int_t activeID, AliITSUGeomTGeo* gm)
a11ef2e4	41	:fActiveID(activeID)
a11ef2e4	42	,fNSensors(0)
68b7631d	43	,fNSensorRows(0)
	44	,fNSensorsPerRow(0)
	45	,fSensVIDtoMatrixID(0)
a11ef2e4	46	,fR(0)
	47	,fRMax(0)
	48	,fRMin(0)
	49	,fZMax(0)
	50	,fZMin(0)
a11ef2e4	51	,fPhiOffs(0)
a11ef2e4	52	,fSensDZInv(0)
68b7631d	53	,fSensDPhiInv(0)
32d38de2	54	,fMaxStep(0.5)
32d38de2	55	,fSensors(0)
a11ef2e4	56	,fITSGeom(gm)
32d38de2	57	,fClusters(0)
a11ef2e4	58	{
	59	// def. c-tor
	60	SetNameTitle(name,name);
32d38de2	61	Build();
a11ef2e4	62	}
	63
	64	//______________________________________________________
	65	AliITSURecoLayer::~AliITSURecoLayer()
	66	{
	67	// def. d-tor
68b7631d	68	delete fSensors;
68b7631d	69	delete[] fSensVIDtoMatrixID;
ee52e7b5	70	if (GetOwnsClusterArray()) delete fClusters;
a11ef2e4	71	}
	72
	73	//______________________________________________________
	74	void AliITSURecoLayer::Print(Option_t* opt) const
	75	{
	76	//print
68b7631d	77	printf("Lr %-15s %d (act:%+d), NSens: %4d in %3d 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);
a11ef2e4	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	//______________________________________________________
32d38de2	86	void AliITSURecoLayer::Build()
a11ef2e4	87	{
a11ef2e4	88	// build internal structures
b515ad7e	89	const double kSafeR = 0.05; // safety margin for Rmin,Rmax of the layer
32d38de2	90	if (fActiveID<0) return;
68b7631d	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];
546d00d8	107	const AliITSsegmentation* kSegm = fITSGeom->GetSegmentation(fActiveID);
a11ef2e4	108	//
a11ef2e4	109	TGeoHMatrix mmod;
a11ef2e4	110	const TGeoHMatrix* mt2l;
a11ef2e4	111	double phiTF,rTF, loc[3]={0,0,0},glo[3];
68b7631d	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);
a11ef2e4	119	double phiMin=1e9,phiMax=-1e9,zMin=1e9,zMax=-1e9;
68b7631d	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)
546d00d8	124	loc[0] = (ix-0.5)*kSegm->Dx(); // +-DX/2
a11ef2e4	125	for (int iy=0;iy<2;iy++) {
546d00d8	126	loc[1] = (iy-0.5)*kSegm->Dy(); // +-DY/2
a11ef2e4	127	for (int iz=0;iz<2;iz++) {
546d00d8	128	loc[2] = (iz-0.5)*kSegm->Dz(); // +-DZ/2
a11ef2e4	129	mmod.LocalToMaster(loc,glo);
a11ef2e4	130	double phi = ATan2(glo[1],glo[0]);
a11ef2e4	131	BringTo02Pi(phi);
68b7631d	132	if (phiMin>1e8) phiMin=phi;
a11ef2e4	133	else if (!OKforPhiMin(phiMin,phi)) phiMin=phi;
68b7631d	134	if (phiMax<-1e8) phiMax=phi;
a11ef2e4	135	else if (!OKforPhiMax(phiMax,phi)) phiMax=phi;
dde91d5d	136	if (glo[2]>zMax) zMax=glo[2];
dde91d5d	137	if (glo[2]<zMin) zMin=glo[2];
a11ef2e4	138	}
	139	}
	140	}
a11ef2e4	141	sens->SetBoundaries(phiMin,phiMax,zMin,zMax);
a11ef2e4	142	}
a11ef2e4	143	}
68b7631d	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());
a11ef2e4	199	//
	200	fRMin = Sqrt(fRMin);
	201	fRMax = Sqrt(fRMax);
	202	fR = 0.5*(fRMin+fRMax);
b515ad7e	203	fRMin -= kSafeR;
b515ad7e	204	fRMax += kSafeR;
a11ef2e4	205	//
a11ef2e4	206	}
08419930	207
a11ef2e4	208	//______________________________________________________
68b7631d	209	Int_t AliITSURecoLayer::FindSensors(const double* impPar, AliITSURecoSens *sensors[kMaxSensMatching])
a11ef2e4	210	{
	211	// find sensors having intersection with track
	212	// impPar contains: lab phi of track, dphi, labZ, dz
5785a9d8	213	//
68b7631d	214	double zMn=impPar[2]-impPar[3], zMx=impPar[2]+impPar[3];
	215	if (zMn>fZMax) return 0;
	216	if (zMx<fZMin) return 0;
a11ef2e4	217	//
68b7631d	218	int zCenID = int((impPar[2]-fZMin)*fSensDZInv);
	219	if (zCenID<0) zCenID = 0;
	220	else if (zCenID>=fNSensorsPerRow) zCenID = fNSensorsPerRow-1;
	221	double phiCn = impPar[0] - fPhiOffs;
	222	// BringTo02Pi(phiCn);
	223	int rowCenID = int(phiCn*fSensDPhiInv);
a11ef2e4	224	//
68b7631d	225	// due to the misalignments the actual sensorID's might be shifted
	226	int res = 0;
	227	AliITSURecoSens* sensPrev=0, *sens = GetSensor(rowCenID,zCenID);
5785a9d8	228	//
68b7631d	229	// printf("Guess: Primary Sensor: phiID: %d zID: %d ->",rowCenID,zCenID); sens->Print();
a11ef2e4	230	//
68b7631d	231	while ( (res=sens->CheckCoverage(impPar[0], impPar[2])) ) {
	232	if (res&AliITSURecoSens::kRight) {if (++rowCenID==fNSensorRows) rowCenID=0;} // neighbor on the right (larger phi)
	233	else if (res&AliITSURecoSens::kLeft) {if (--rowCenID<0) rowCenID = fNSensorRows-1;} // neighbor on the left (smaller phi)
	234	if (res&AliITSURecoSens::kUp) {if (++zCenID==fNSensorsPerRow) zCenID = fNSensorsPerRow-1;} // neighbor at larger Z (if any)
	235	else if (res&AliITSURecoSens::kDown) {if (--zCenID<0) zCenID = 0;} // neighbor at smaller Z (if any)
	236	//
	237	AliITSURecoSens* sensAlt = GetSensor(rowCenID,zCenID);
	238	if (sensAlt==sens \|\| sensAlt==sensPrev) break; // there is no better neighbor (z edge) or the point falls in dead area
	239	//
	240	sensPrev = sens;
	241	sens = sensAlt;
	242	}
	243	// printf("Found: Primary Sensor: phiID: %d zID: %d ->",rowCenID,zCenID); sens->Print();
a11ef2e4	244	//
68b7631d	245	int nFnd = 0;
68b7631d	246	sensors[nFnd++] = sens;
a11ef2e4	247	//
68b7631d	248	double phiMn = impPar[0]-impPar[1], phiMx = impPar[0]+impPar[1];
	249	BringTo02Pi(phiMn);
	250	BringTo02Pi(phiMx);
a11ef2e4	251	//
68b7631d	252	const int kNNeighb = 8;
	253	const int kCheckNeighb[2][kNNeighb] = { // phi and Z neighbours to check
	254	{ 1, 1, 0,-1,-1,-1, 0, 1},
	255	{ 0, 1, 1, 1, 0,-1,-1,-1}
	256	};
	257	//
	258	// printf("Search: %+.4f %+.4f \| %+.4f %+.4f\n",phiMn,phiMx, zMn,zMx);
	259	for (int inb=kNNeighb;inb--;) {
	260	int idz = kCheckNeighb[1][inb];
	261	int iz = zCenID + idz;
	262	// printf("#%d dp:%+d dz:%+d IZ: %d\n",inb, kCheckNeighb[0][inb], kCheckNeighb[1][inb], iz);
	263	//
	264	if (iz<0 \|\| iz>=fNSensorsPerRow) continue;
	265	int idphi = kCheckNeighb[0][inb];
	266	int iphi = rowCenID + idphi;
	267	if (iphi<0) iphi += fNSensorRows;
	268	else if (iphi>=fNSensorRows) iphi -= fNSensorRows;
	269	sens = GetSensor(iphi,iz);
	270	//
	271	if (idz>0) {if (zMx<sens->GetZMin()) continue;}
	272	else if (idz<0) {if (zMn>sens->GetZMax()) continue;}
	273	//
	274	// Z range matches
	275	if (idphi>0) {if (!OKforPhiMin(sens->GetPhiMin(),phiMx)) continue;}
	276	else if (idphi<0) {if (!OKforPhiMax(sens->GetPhiMax(),phiMn)) continue;}
	277	//
	278	// printf("Add %d\n",nFnd);
	279	sensors[nFnd++] = sens;
	280	if (nFnd==kMaxSensMatching) break;
	281	}
	282	return nFnd;
a11ef2e4	283	}
68b7631d	284
c03e4f8a	285	//*/
5785a9d8	286	/*
	287	Int_t AliITSURecoLayer::FindSensors(const double* impPar, AliITSURecoSens *sensors[AliITSURecoSens::kNNeighbors],int mcLab)
	288	{
	289	// find sensors having intersection with track
	290	// impPar contains: lab phi of track, dphi, labZ, dz
	291
	292	//tmp>>>
	293	int nFnd = 0;
	294	int fndSens[50];
	295	if (mcLab>=0) { // find correct sensors from MC info
	296	int ncl = GetNClusters();
	297	for (int icl=ncl;icl--;) {
	298	AliCluster* cl = GetCluster(icl);
	299	for (int ilb=0;ilb<3;ilb++) {
	300	if (cl->GetLabel(ilb)<0) break;
	301	if (cl->GetLabel(ilb)==mcLab) {fndSens[nFnd++] = cl->GetVolumeId(); break;}
	302	}
	303	}
	304	if (nFnd>0) {
852af72e	305	Int_t layS,staS,sensS;
5785a9d8	306	for (int is=0;is<nFnd;is++) {
852af72e	307	fITSGeom->GetChipId(fndSens[is],layS,staS,sensS);
852af72e	308	printf("SNMC#%d(%d): %d %d %d \| ",is,mcLab,layS,staS,sensS); GetSensorFromID(fndSens[is])->Print();
5785a9d8	309	}
	310	}
	311	}
	312
	313	//tmp<<<
	314
	315	double z = impPar[2];
	316	if (z>fZMax+impPar[3]) {
	317	if (nFnd>0) printf("MissedSens!!!\n");
	318	return 0; // outside of Z coverage
	319	}
	320	z -= fZMin;
	321	if (z<-impPar[3]) {
	322	if (nFnd>0) printf("MissedSens!!!\n");
	323	return 0; // outside of Z coverage
	324	}
852af72e	325	int sensInSta = int(z*fSensDZInv);
	326	if (sensInSta<0) sensInSta = 0;
	327	else if (sensInSta>=fNSensInStave) sensInSta = fNSensInStave-1;
5785a9d8	328	//
	329	double phi = impPar[0] - fPhiOffs;
	330	BringTo02Pi(phi);
68b7631d	331	int staID = int(phi*fSensDPhiInv); // stave id
5785a9d8	332	int nsens = 0;
5785a9d8	333	//
852af72e	334	AliITSURecoSens* sensN,sens = GetSensor(staIDfNSensInStave+sensInSta);
5785a9d8	335	sensors[nsens++] = sens;
	336	//
	337	// check neighbours
	338	double zMn=impPar[2]-impPar[3], zMx=impPar[2]+impPar[3], phiMn=impPar[0]-impPar[1], phiMx=impPar[0]+impPar[1];
	339	BringTo02Pi(phiMn);
	340	BringTo02Pi(phiMx);
	341	//
	342	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbR)); // neighbor on the right (smaller phi)
	343	if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax())) sensors[nsens++] = sensN;
	344	//
	345	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTR)); // neighbor on the top right (smaller phi, larger Z)
	346	if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
	347	//
	348	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbT)); // neighbor on the top (larger Z)
	349	if (sensN && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
	350	//
	351	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTL)); // neighbor on the top left (larger Z, larger phi)
	352	if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
	353	//
	354	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbL)); // neighbor on the left (larger phi)
	355	if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin())) sensors[nsens++] = sensN;
	356	//
	357	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBL)); // neighbor on the bottom left (smaller Z, larger phi)
	358	if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
	359	//
	360	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbB)); // neighbor on the bottom (smaller Z)
	361	if (sensN && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
	362	//
	363	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBR)); // neighbor on the bottom right (smaller Z, smaller phi)
	364	if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
	365	//
	366	if (mcLab>=0) {
852af72e	367	Int_t layS,staS,sensS;
5785a9d8	368	printf("Found %d sensors for phi %.3f : %.3f \| Z %.4f %.4f\n", nsens,phiMn,phiMx,zMn,zMx);
5785a9d8	369	for (int is=0;is<nsens;is++) {
852af72e	370	fITSGeom->GetChipId(sensors[is]->GetID()+fITSGeom->GetFirstModIndex(fActiveID),layS,staS,sensS);
852af72e	371	printf("*SNF#%d: %d %d %d \| ",is,layS,staS,sensS); sensors[is]->Print();
5785a9d8	372	}
	373	for (int ism=0;ism<nFnd;ism++) {
	374	AliITSURecoSens* snMC = GetSensorFromID(fndSens[ism]);
	375	Bool_t ok=kFALSE;
	376	for (int isf=0;isf<nsens;isf++) {
	377	if (snMC==sensors[isf]) {ok=kTRUE;break;}
	378	}
	379	if (!ok) printf("MissedSens %d!!!\n",ism);
	380	}
	381	}
	382	return nsens;
	383	}
	384	*/
32d38de2	385	//______________________________________________________
	386	void AliITSURecoLayer::ProcessClusters(Int_t mode)
	387	{
	388	// register in each sensor of the layer its cluster.
	389	// the clusters of the layer must be sorted per sensor
	390	int ncl = fClusters->GetEntriesFast();
	391	int curSensID = -1;
0091e9f0	392	for (int i=fNSensors;i--;) GetSensor(i)->SetNClusters(0);
32d38de2	393	AliITSURecoSens* curSens = 0;
	394	for (int icl=0;icl<ncl;icl++) {
	395	AliITSUClusterPix* cl = (AliITSUClusterPix*) fClusters->UncheckedAt(icl);
	396	cl->GoToFrameTrk();
	397	int vID = cl->GetVolumeId();
	398	if (vID<curSensID) {AliFatal("Clusters are not sorted in increasing sensorID");}
	399	if (vID>curSensID) {
	400	if (curSens) curSens->ProcessClusters(mode); // prepare clusters for reconstruction
6a6a65db	401	curSens = GetSensorFromID(vID); //GetSensor(vID - fITSGeom->GetFirstChipIndex(fActiveID));
32d38de2	402	curSensID = vID;
	403	curSens->SetFirstClusterId(icl);
	404	}
	405	curSens->IncNClusters();
	406	}
	407	if (curSens) curSens->ProcessClusters(mode); // last sensor was not processed yet
	408	//
	409	}
dde91d5d	410
	411	//______________________________________________________
	412	Bool_t AliITSURecoLayer::IsEqual(const TObject* obj) const
	413	{
	414	// check if layers are equal in R
44785f3e	415	const AliITSURecoLayer* lr = (const AliITSURecoLayer*)obj;
dde91d5d	416	return Abs(lr->GetR()-GetR())<1e-6 ? kTRUE : kFALSE;
	417	}
	418
	419	//______________________________________________________
	420	Int_t AliITSURecoLayer::Compare(const TObject* obj) const
	421	{
	422	// compare two layers
44785f3e	423	const AliITSURecoLayer* lr = (const AliITSURecoLayer*)obj;
dde91d5d	424	double dr = GetR() - lr->GetR();
	425	if (Abs(dr)<1e-6) return 0;
	426	return dr>0 ? 1:-1;
	427	//
	428	}
8b16dbae	429
	430	//_________________________________________________________________
	431	AliITSURecoSens* AliITSURecoLayer::GetSensorFromID(Int_t i) const
	432	{
	433	// get sensor from its global id
852af72e	434	i -= fITSGeom->GetFirstChipIndex(fActiveID);
852af72e	435	if (i<0\|\|i>=fNSensors) AliFatal(Form("Sensor with id=%d is not in layer %d",i+fITSGeom->GetFirstChipIndex(fActiveID),fActiveID));
68b7631d	436	return GetSensor(SensVIDtoMatrixID(i));
8b16dbae	437	}