[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)
  ,fNSensInStave(0)
  ,fNStaves(0)
  ,fR(0)
  ,fRMax(0)
  ,fRMin(0)
  ,fZMax(0)
  ,fZMin(0)
  ,fPhiStaMax(0)
  ,fPhiStaMin(0)
  ,fPhiOffs(0)
  ,fSensDZInv(0)
  ,fDPhiStaInv(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)
  ,fNSensInStave(0)
  ,fNStaves(0)
  ,fR(0)
  ,fRMax(0)
  ,fRMin(0)
  ,fZMax(0)
  ,fZMin(0)
  ,fPhiStaMax(0)
  ,fPhiStaMin(0)
  ,fPhiOffs(0)
  ,fSensDZInv(0)
  ,fDPhiStaInv(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[] fPhiStaMax;
  delete[] fPhiStaMin;
  if (GetOwnsClusterArray()) delete fClusters;
}

//______________________________________________________
void AliITSURecoLayer::Print(Option_t* opt) const			      
{
  //print 
  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);
}

//______________________________________________________
void AliITSURecoLayer::Build()
{
  // build internal structures
  const double kSafeR = 0.05; // safety margin for Rmin,Rmax of the layer
  if (fActiveID<0) return;
  fNStaves = fITSGeom->GetNStaves(fActiveID);
  fNSensInStave = fITSGeom->GetNChipsPerModule(fActiveID);
  fNSensors = fNStaves*fNSensInStave;
  fSensors = new AliITSURecoSens*[fNSensors];
  const AliITSsegmentation* kSegm = fITSGeom->GetSegmentation(fActiveID);
  //
  // name layer according its active id, detector type and segmentation tyoe
  TGeoHMatrix mmod;
  const TGeoHMatrix* mt2l;
  fRMin=fZMin=1e9;
  fRMax=fZMax=-1e9;
  double phiTF,rTF, loc[3]={0,0,0},glo[3];
  fNSensors = 0;
  fPhiStaMin = new Double_t[fNStaves];
  fPhiStaMax = new Double_t[fNStaves];
  fSensDZInv = 0;
  fDPhiStaInv = fNStaves/TwoPi();
  //
  for (int ild=0;ild<fNStaves;ild++) {
    fPhiStaMin[ild] = 1e9;
    fPhiStaMax[ild] = -1e9;
    //
    for (int idt=0;idt<fNSensInStave;idt++) {
      AliITSURecoSens* sens = new AliITSURecoSens(fNSensors++);
      fSensors[ild*fNSensInStave+idt] = sens;
      //
      double phiMin=1e9,phiMax=-1e9,zMin=1e9,zMax=-1e9;
      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++) {
	  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];
	  }
	}
      }
      sens->SetBoundaries(phiMin,phiMax,zMin,zMax);
      mt2l = fITSGeom->GetMatrixT2L(fActiveID,ild,idt);
      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);
      //
      if      (fPhiStaMin[ild]>1e8)  fPhiStaMin[ild] = phiMin;
      else if (!OKforPhiMin(fPhiStaMin[ild],phiMin)) fPhiStaMin[ild] = phiMin;
      if      (fPhiStaMax[ild]<-1e8) fPhiStaMax[ild] = phiMax;
      else if (!OKforPhiMax(fPhiStaMax[ild],phiMax)) fPhiStaMax[ild] = phiMax;
      if (fZMin>zMin) fZMin = zMin;
      if (fZMax<zMax) fZMax = zMax;
      //
      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 = fNSensInStave>0 ? fSensDZInv/(fNSensInStave-1)/fNStaves : fZMax-fZMin;
  fSensDZInv = 1./dz;

  const int kNBId[3][3] = { 
    {AliITSURecoSens::kNghbBL,AliITSURecoSens::kNghbB,AliITSURecoSens::kNghbBR},
    {AliITSURecoSens::kNghbL,          -1            ,AliITSURecoSens::kNghbR },
    {AliITSURecoSens::kNghbTL,AliITSURecoSens::kNghbT,AliITSURecoSens::kNghbTR}
  };

  // add neighbours info
  double zTol = 0.45*dz, phiTol = 0.45*TwoPi()/fNStaves;
  for (int ild=0;ild<fNStaves;ild++) {
    for (int idt=0;idt<fNSensInStave;idt++) {
      AliITSURecoSens* sens = GetSensor(ild,idt);
      //
      for (int ils=-1;ils<=1;ils++) {
	int ildN = ild+ils;  // staves of neighbouring sensors
	if (ildN<0) ildN = fNStaves-1; else if (ildN==fNStaves) ildN = 0;
	for (int ids=-1;ids<=1;ids++) {
	  int idtN = idt+ids;
	  if (idtN<0 || idtN==fNSensInStave || (ids==0&&ils==0)) continue;
	  AliITSURecoSens* sensN = GetSensor(ildN,idtN); // potential neighbor
	  int neighbID = ildN*fNSensInStave+idtN;
	  //	  
	  int zType = 1;  // side
	  if (sens->GetZMin()-zTol  > sensN->GetZMax()) continue; // too large distance
	  if (sensN->GetZMin()-zTol > sens->GetZMax() ) continue; // too large distance
	  if      (sens->GetZMin()-zTol>sensN->GetZMin()) zType =  0;     // bottom
	  else if (sensN->GetZMin()-zTol>sens->GetZMin()) zType =  2;     // top
	  //
	  int phiType = 1;

	  double phiTstMn = sensN->GetPhiMin()-phiTol;
	  BringTo02Pi(phiTstMn);
	  if (!OKforPhiMax(sens->GetPhiMax(),phiTstMn)) continue; // too large angle	  
	  double phiTstMx = sensN->GetPhiMax()+phiTol;	  
	  BringTo02Pi(phiTstMx);
	  if (!OKforPhiMin(sens->GetPhiMin(),phiTstMx)) continue; // too large angle
	  //
	  phiTstMn = sensN->GetPhiMin()+phiTol;
	  BringTo02Pi(phiTstMn);
	  phiTstMx = sensN->GetPhiMax()-phiTol;	  
	  BringTo02Pi(phiTstMx);
	  if      (!OKforPhiMax(sens->GetPhiMax(),phiTstMx)) phiType = 0; // left
	  else if (!OKforPhiMin(sens->GetPhiMin(),phiTstMn)) phiType = 2; // right
	  //
	  sens->SetNeighborID(kNBId[zType][phiType], neighbID);
	} // phi scan
      } // z scan
    } // sensors
  } // staves
  //
}

//______________________________________________________
Int_t AliITSURecoLayer::FindSensors(const double* impPar, AliITSURecoSens *sensors[AliITSURecoSens::kNNeighbors])
{
  // 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 0; // outside of Z coverage
  z -= fZMin;
  if (z<-impPar[3]) 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*fDPhiStaInv);  // 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;
  //
  return nsens;
}
//*/
/*
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*fDPhiStaInv);  // 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   = 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 (AliITSURecoSens*)fSensors[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)
	17	:fActiveID(-1)
	18	,fNSensors(0)
852af72e	19	,fNSensInStave(0)
852af72e	20	,fNStaves(0)
a11ef2e4	21	,fR(0)
	22	,fRMax(0)
	23	,fRMin(0)
	24	,fZMax(0)
	25	,fZMin(0)
852af72e	26	,fPhiStaMax(0)
852af72e	27	,fPhiStaMin(0)
a11ef2e4	28	,fPhiOffs(0)
a11ef2e4	29	,fSensDZInv(0)
852af72e	30	,fDPhiStaInv(0)
32d38de2	31	,fMaxStep(0.5)
a11ef2e4	32	,fSensors(0)
a11ef2e4	33	,fITSGeom(0)
32d38de2	34	,fClusters(0)
a11ef2e4	35	{
	36	// def. c-tor
	37	SetNameTitle(name,name);
a11ef2e4	38	}
	39
	40	//______________________________________________________
32d38de2	41	AliITSURecoLayer::AliITSURecoLayer(const char* name, Int_t activeID, AliITSUGeomTGeo* gm)
a11ef2e4	42	:fActiveID(activeID)
a11ef2e4	43	,fNSensors(0)
852af72e	44	,fNSensInStave(0)
852af72e	45	,fNStaves(0)
a11ef2e4	46	,fR(0)
	47	,fRMax(0)
	48	,fRMin(0)
	49	,fZMax(0)
	50	,fZMin(0)
852af72e	51	,fPhiStaMax(0)
852af72e	52	,fPhiStaMin(0)
a11ef2e4	53	,fPhiOffs(0)
a11ef2e4	54	,fSensDZInv(0)
852af72e	55	,fDPhiStaInv(0)
32d38de2	56	,fMaxStep(0.5)
32d38de2	57	,fSensors(0)
a11ef2e4	58	,fITSGeom(gm)
32d38de2	59	,fClusters(0)
a11ef2e4	60	{
	61	// def. c-tor
	62	SetNameTitle(name,name);
32d38de2	63	Build();
a11ef2e4	64	}
	65
	66	//______________________________________________________
	67	AliITSURecoLayer::~AliITSURecoLayer()
	68	{
	69	// def. d-tor
32d38de2	70	delete[] fSensors;
852af72e	71	delete[] fPhiStaMax;
852af72e	72	delete[] fPhiStaMin;
ee52e7b5	73	if (GetOwnsClusterArray()) delete fClusters;
a11ef2e4	74	}
	75
	76	//______________________________________________________
	77	void AliITSURecoLayer::Print(Option_t* opt) const
	78	{
	79	//print
173b3073	80	printf("Lr %-15s %d (act:%+d), NSens: %4d \| MaxStep:%.2f ",GetName(),GetID(),GetActiveID(),GetNSensors(),fMaxStep);
dde91d5d	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);
a11ef2e4	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	//______________________________________________________
32d38de2	87	void AliITSURecoLayer::Build()
a11ef2e4	88	{
a11ef2e4	89	// build internal structures
b515ad7e	90	const double kSafeR = 0.05; // safety margin for Rmin,Rmax of the layer
32d38de2	91	if (fActiveID<0) return;
852af72e	92	fNStaves = fITSGeom->GetNStaves(fActiveID);
	93	fNSensInStave = fITSGeom->GetNChipsPerModule(fActiveID);
	94	fNSensors = fNStaves*fNSensInStave;
32d38de2	95	fSensors = new AliITSURecoSens*[fNSensors];
546d00d8	96	const AliITSsegmentation* kSegm = fITSGeom->GetSegmentation(fActiveID);
a11ef2e4	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;
852af72e	105	fPhiStaMin = new Double_t[fNStaves];
852af72e	106	fPhiStaMax = new Double_t[fNStaves];
a11ef2e4	107	fSensDZInv = 0;
852af72e	108	fDPhiStaInv = fNStaves/TwoPi();
a11ef2e4	109	//
852af72e	110	for (int ild=0;ild<fNStaves;ild++) {
	111	fPhiStaMin[ild] = 1e9;
	112	fPhiStaMax[ild] = -1e9;
a11ef2e4	113	//
852af72e	114	for (int idt=0;idt<fNSensInStave;idt++) {
a11ef2e4	115	AliITSURecoSens* sens = new AliITSURecoSens(fNSensors++);
852af72e	116	fSensors[ild*fNSensInStave+idt] = sens;
a11ef2e4	117	//
a11ef2e4	118	double phiMin=1e9,phiMax=-1e9,zMin=1e9,zMax=-1e9;
b515ad7e	119	mmod = *fITSGeom->GetMatrixSens(fActiveID,ild,idt);
a11ef2e4	120	for (int ix=0;ix<2;ix++) {
546d00d8	121	loc[0] = (ix-0.5)*kSegm->Dx(); // +-DX/2
a11ef2e4	122	for (int iy=0;iy<2;iy++) {
546d00d8	123	loc[1] = (iy-0.5)*kSegm->Dy(); // +-DY/2
a11ef2e4	124	for (int iz=0;iz<2;iz++) {
546d00d8	125	loc[2] = (iz-0.5)*kSegm->Dz(); // +-DZ/2
a11ef2e4	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;
dde91d5d	137	if (glo[2]>zMax) zMax=glo[2];
dde91d5d	138	if (glo[2]<zMin) zMin=glo[2];
a11ef2e4	139	}
	140	}
	141	}
a11ef2e4	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	//
852af72e	154	if (fPhiStaMin[ild]>1e8) fPhiStaMin[ild] = phiMin;
	155	else if (!OKforPhiMin(fPhiStaMin[ild],phiMin)) fPhiStaMin[ild] = phiMin;
	156	if (fPhiStaMax[ild]<-1e8) fPhiStaMax[ild] = phiMax;
	157	else if (!OKforPhiMax(fPhiStaMax[ild],phiMax)) fPhiStaMax[ild] = phiMax;
a11ef2e4	158	if (fZMin>zMin) fZMin = zMin;
	159	if (fZMax<zMax) fZMax = zMax;
	160	//
dde91d5d	161	if (idt>0) fSensDZInv += zMax - GetSensor(ild,idt-1)->GetZMax(); // z interval to previous
a11ef2e4	162	}
	163	}
	164	//
	165	fRMin = Sqrt(fRMin);
	166	fRMax = Sqrt(fRMax);
	167	fR = 0.5*(fRMin+fRMax);
b515ad7e	168	fRMin -= kSafeR;
b515ad7e	169	fRMax += kSafeR;
852af72e	170	double dz = fNSensInStave>0 ? fSensDZInv/(fNSensInStave-1)/fNStaves : fZMax-fZMin;
a11ef2e4	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
852af72e	180	double zTol = 0.45dz, phiTol = 0.45TwoPi()/fNStaves;
	181	for (int ild=0;ild<fNStaves;ild++) {
	182	for (int idt=0;idt<fNSensInStave;idt++) {
a11ef2e4	183	AliITSURecoSens* sens = GetSensor(ild,idt);
	184	//
	185	for (int ils=-1;ils<=1;ils++) {
852af72e	186	int ildN = ild+ils; // staves of neighbouring sensors
852af72e	187	if (ildN<0) ildN = fNStaves-1; else if (ildN==fNStaves) ildN = 0;
a11ef2e4	188	for (int ids=-1;ids<=1;ids++) {
a11ef2e4	189	int idtN = idt+ids;
852af72e	190	if (idtN<0 \|\| idtN==fNSensInStave \|\| (ids==0&&ils==0)) continue;
a11ef2e4	191	AliITSURecoSens* sensN = GetSensor(ildN,idtN); // potential neighbor
852af72e	192	int neighbID = ildN*fNSensInStave+idtN;
a11ef2e4	193	//
a11ef2e4	194	int zType = 1; // side
32d38de2	195	if (sens->GetZMin()-zTol > sensN->GetZMax()) continue; // too large distance
32d38de2	196	if (sensN->GetZMin()-zTol > sens->GetZMax() ) continue; // too large distance
a11ef2e4	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;
32d38de2	201
a11ef2e4	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
852af72e	220	} // staves
a11ef2e4	221	//
a11ef2e4	222	}
08419930	223
a11ef2e4	224	//______________________________________________________
32d38de2	225	Int_t AliITSURecoLayer::FindSensors(const double* impPar, AliITSURecoSens *sensors[AliITSURecoSens::kNNeighbors])
a11ef2e4	226	{
	227	// find sensors having intersection with track
	228	// impPar contains: lab phi of track, dphi, labZ, dz
5785a9d8	229	//
a11ef2e4	230	double z = impPar[2];
0091e9f0	231	if (z>fZMax+impPar[3]) return 0; // outside of Z coverage
a11ef2e4	232	z -= fZMin;
0091e9f0	233	if (z<-impPar[3]) return 0; // outside of Z coverage
852af72e	234	int sensInSta = int(z*fSensDZInv);
	235	if (sensInSta<0) sensInSta = 0;
	236	else if (sensInSta>=fNSensInStave) sensInSta = fNSensInStave-1;
a11ef2e4	237	//
	238	double phi = impPar[0] - fPhiOffs;
	239	BringTo02Pi(phi);
852af72e	240	int staID = int(phi*fDPhiStaInv); // stave id
a11ef2e4	241	int nsens = 0;
a11ef2e4	242	//
852af72e	243	AliITSURecoSens* sensN,sens = GetSensor(staIDfNSensInStave+sensInSta);
a11ef2e4	244	sensors[nsens++] = sens;
5785a9d8	245	//
a11ef2e4	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	//
173b3073	251	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbR)); // neighbor on the right (smaller phi)
173b3073	252	if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax())) sensors[nsens++] = sensN;
a11ef2e4	253	//
173b3073	254	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTR)); // neighbor on the top right (smaller phi, larger Z)
173b3073	255	if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
a11ef2e4	256	//
173b3073	257	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbT)); // neighbor on the top (larger Z)
173b3073	258	if (sensN && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
a11ef2e4	259	//
173b3073	260	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTL)); // neighbor on the top left (larger Z, larger phi)
173b3073	261	if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
a11ef2e4	262	//
173b3073	263	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbL)); // neighbor on the left (larger phi)
173b3073	264	if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin())) sensors[nsens++] = sensN;
a11ef2e4	265	//
173b3073	266	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBL)); // neighbor on the bottom left (smaller Z, larger phi)
173b3073	267	if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
a11ef2e4	268	//
173b3073	269	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbB)); // neighbor on the bottom (smaller Z)
a11ef2e4	270	if (sensN && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
a11ef2e4	271	//
173b3073	272	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBR)); // neighbor on the bottom right (smaller Z, smaller phi)
173b3073	273	if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
a11ef2e4	274	//
	275	return nsens;
	276	}
c03e4f8a	277	//*/
5785a9d8	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) {
852af72e	297	Int_t layS,staS,sensS;
5785a9d8	298	for (int is=0;is<nFnd;is++) {
852af72e	299	fITSGeom->GetChipId(fndSens[is],layS,staS,sensS);
852af72e	300	printf("SNMC#%d(%d): %d %d %d \| ",is,mcLab,layS,staS,sensS); GetSensorFromID(fndSens[is])->Print();
5785a9d8	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	}
852af72e	317	int sensInSta = int(z*fSensDZInv);
	318	if (sensInSta<0) sensInSta = 0;
	319	else if (sensInSta>=fNSensInStave) sensInSta = fNSensInStave-1;
5785a9d8	320	//
	321	double phi = impPar[0] - fPhiOffs;
	322	BringTo02Pi(phi);
852af72e	323	int staID = int(phi*fDPhiStaInv); // stave id
5785a9d8	324	int nsens = 0;
5785a9d8	325	//
852af72e	326	AliITSURecoSens* sensN,sens = GetSensor(staIDfNSensInStave+sensInSta);
5785a9d8	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) {
852af72e	359	Int_t layS,staS,sensS;
5785a9d8	360	printf("Found %d sensors for phi %.3f : %.3f \| Z %.4f %.4f\n", nsens,phiMn,phiMx,zMn,zMx);
5785a9d8	361	for (int is=0;is<nsens;is++) {
852af72e	362	fITSGeom->GetChipId(sensors[is]->GetID()+fITSGeom->GetFirstModIndex(fActiveID),layS,staS,sensS);
852af72e	363	printf("*SNF#%d: %d %d %d \| ",is,layS,staS,sensS); sensors[is]->Print();
5785a9d8	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	*/
32d38de2	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;
0091e9f0	384	for (int i=fNSensors;i--;) GetSensor(i)->SetNClusters(0);
32d38de2	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
852af72e	393	curSens = GetSensor(vID - fITSGeom->GetFirstChipIndex(fActiveID));
32d38de2	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	}
dde91d5d	402
	403	//______________________________________________________
	404	Bool_t AliITSURecoLayer::IsEqual(const TObject* obj) const
	405	{
	406	// check if layers are equal in R
44785f3e	407	const AliITSURecoLayer* lr = (const AliITSURecoLayer*)obj;
dde91d5d	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
44785f3e	415	const AliITSURecoLayer* lr = (const AliITSURecoLayer*)obj;
dde91d5d	416	double dr = GetR() - lr->GetR();
	417	if (Abs(dr)<1e-6) return 0;
	418	return dr>0 ? 1:-1;
	419	//
	420	}
8b16dbae	421
	422	//_________________________________________________________________
	423	AliITSURecoSens* AliITSURecoLayer::GetSensorFromID(Int_t i) const
	424	{
	425	// get sensor from its global id
852af72e	426	i -= fITSGeom->GetFirstChipIndex(fActiveID);
852af72e	427	if (i<0\|\|i>=fNSensors) AliFatal(Form("Sensor with id=%d is not in layer %d",i+fITSGeom->GetFirstChipIndex(fActiveID),fActiveID));
8b16dbae	428	return (AliITSURecoSens*)fSensors[i];
8b16dbae	429	}