[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)
  ,fNSensInLadder(0)
  ,fNLadders(0)
  ,fR(0)
  ,fRMax(0)
  ,fRMin(0)
  ,fZMax(0)
  ,fZMin(0)
  ,fPhiLadMax(0)
  ,fPhiLadMin(0)
  ,fPhiOffs(0)
  ,fSensDZInv(0)
  ,fDPhiLadInv(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)
  ,fNSensInLadder(0)
  ,fNLadders(0)
  ,fR(0)
  ,fRMax(0)
  ,fRMin(0)
  ,fZMax(0)
  ,fZMin(0)
  ,fPhiLadMax(0)
  ,fPhiLadMin(0)
  ,fPhiOffs(0)
  ,fSensDZInv(0)
  ,fDPhiLadInv(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[] fPhiLadMax;
  delete[] fPhiLadMin;
  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;
  fNLadders = fITSGeom->GetNLadders(fActiveID);
  fNSensInLadder = fITSGeom->GetNDetectors(fActiveID);
  fNSensors = fNLadders*fNSensInLadder;
  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;
  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;
    fPhiLadMax[ild] = -1e9;
    //
    for (int idt=0;idt<fNSensInLadder;idt++) {
      AliITSURecoSens* sens = new AliITSURecoSens(fNSensors++);
      fSensors[ild*fNSensInLadder+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      (fPhiLadMin[ild]>1e8)  fPhiLadMin[ild] = phiMin;
      else if (!OKforPhiMin(fPhiLadMin[ild],phiMin)) fPhiLadMin[ild] = phiMin;
      if      (fPhiLadMax[ild]<-1e8) fPhiLadMax[ild] = phiMax;
      else if (!OKforPhiMax(fPhiLadMax[ild],phiMax)) fPhiLadMax[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 = fNSensInLadder>0 ? fSensDZInv/(fNSensInLadder-1)/fNLadders : 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()/fNLadders;
  for (int ild=0;ild<fNLadders;ild++) {
    for (int idt=0;idt<fNSensInLadder;idt++) {
      AliITSURecoSens* sens = GetSensor(ild,idt);
      //
      for (int ils=-1;ils<=1;ils++) {
	int ildN = ild+ils;  // ladders of neighbouring sensors
	if (ildN<0) ildN = fNLadders-1; else if (ildN==fNLadders) ildN = 0;
	for (int ids=-1;ids<=1;ids++) {
	  int idtN = idt+ids;
	  if (idtN<0 || idtN==fNSensInLadder || (ids==0&&ils==0)) continue;
	  AliITSURecoSens* sensN = GetSensor(ildN,idtN); // potential neighbor
	  int neighbID = ildN*fNSensInLadder+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
  } // ladders
  //
}

//______________________________________________________
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 sensInLad = int(z*fSensDZInv);
  if      (sensInLad<0) sensInLad = 0;
  else if (sensInLad>=fNSensInLadder) sensInLad = fNSensInLadder-1;
  //
  double phi = impPar[0] - fPhiOffs;
  BringTo02Pi(phi);
  int ladID = int(phi*fDPhiLadInv);  // ladder id
  int nsens = 0;
  //
  AliITSURecoSens* sensN,*sens = GetSensor(ladID*fNSensInLadder+sensInLad);
  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,ladS,sensS;
      for (int is=0;is<nFnd;is++) {
	fITSGeom->GetModuleId(fndSens[is],layS,ladS,sensS);
	printf("SNMC#%d(%d): %d %d %d | ",is,mcLab,layS,ladS,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 sensInLad = int(z*fSensDZInv);
  if      (sensInLad<0) sensInLad = 0;
  else if (sensInLad>=fNSensInLadder) sensInLad = fNSensInLadder-1;
  //
  double phi = impPar[0] - fPhiOffs;
  BringTo02Pi(phi);
  int ladID = int(phi*fDPhiLadInv);  // ladder id
  int nsens = 0;
  //
  AliITSURecoSens* sensN,*sens = GetSensor(ladID*fNSensInLadder+sensInLad);
  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,ladS,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->GetModuleId(sensors[is]->GetID()+fITSGeom->GetFirstModIndex(fActiveID),layS,ladS,sensS);
      printf("*SNF#%d: %d %d %d | ",is,layS,ladS,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->GetFirstModIndex(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->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];
}
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)
	19	,fNSensInLadder(0)
	20	,fNLadders(0)
	21	,fR(0)
	22	,fRMax(0)
	23	,fRMin(0)
	24	,fZMax(0)
	25	,fZMin(0)
	26	,fPhiLadMax(0)
	27	,fPhiLadMin(0)
	28	,fPhiOffs(0)
	29	,fSensDZInv(0)
	30	,fDPhiLadInv(0)
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)
	43	,fNSensors(0)
	44	,fNSensInLadder(0)
	45	,fNLadders(0)
	46	,fR(0)
	47	,fRMax(0)
	48	,fRMin(0)
	49	,fZMax(0)
	50	,fZMin(0)
	51	,fPhiLadMax(0)
	52	,fPhiLadMin(0)
	53	,fPhiOffs(0)
	54	,fSensDZInv(0)
	55	,fDPhiLadInv(0)
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;
	71	delete[] fPhiLadMax;
	72	delete[] fPhiLadMin;
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;
a11ef2e4	92	fNLadders = fITSGeom->GetNLadders(fActiveID);
a11ef2e4	93	fNSensInLadder = fITSGeom->GetNDetectors(fActiveID);
32d38de2	94	fNSensors = fNLadders*fNSensInLadder;
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;
	105	fPhiLadMin = new Double_t[fNLadders];
	106	fPhiLadMax = new Double_t[fNLadders];
	107	fSensDZInv = 0;
173b3073	108	fDPhiLadInv = fNLadders/TwoPi();
a11ef2e4	109	//
	110	for (int ild=0;ild<fNLadders;ild++) {
	111	fPhiLadMin[ild] = 1e9;
	112	fPhiLadMax[ild] = -1e9;
	113	//
	114	for (int idt=0;idt<fNSensInLadder;idt++) {
	115	AliITSURecoSens* sens = new AliITSURecoSens(fNSensors++);
dde91d5d	116	fSensors[ild*fNSensInLadder+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	//
	154	if (fPhiLadMin[ild]>1e8) fPhiLadMin[ild] = phiMin;
	155	else if (!OKforPhiMin(fPhiLadMin[ild],phiMin)) fPhiLadMin[ild] = phiMin;
	156	if (fPhiLadMax[ild]<-1e8) fPhiLadMax[ild] = phiMax;
	157	else if (!OKforPhiMax(fPhiLadMax[ild],phiMax)) fPhiLadMax[ild] = phiMax;
	158	if (fZMin>zMin) fZMin = zMin;
	159	if (fZMax<zMax) fZMax = zMax;
	160	//
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;
a11ef2e4	170	double dz = fNSensInLadder>0 ? fSensDZInv/(fNSensInLadder-1)/fNLadders : fZMax-fZMin;
	171	fSensDZInv = 1./dz;
	172
	173	const int kNBId[3][3] = {
	174	{AliITSURecoSens::kNghbBL,AliITSURecoSens::kNghbB,AliITSURecoSens::kNghbBR},
	175	{AliITSURecoSens::kNghbL, -1 ,AliITSURecoSens::kNghbR },
	176	{AliITSURecoSens::kNghbTL,AliITSURecoSens::kNghbT,AliITSURecoSens::kNghbTR}
	177	};
	178
	179	// add neighbours info
	180	double zTol = 0.45dz, phiTol = 0.45TwoPi()/fNLadders;
	181	for (int ild=0;ild<fNLadders;ild++) {
	182	for (int idt=0;idt<fNSensInLadder;idt++) {
	183	AliITSURecoSens* sens = GetSensor(ild,idt);
	184	//
	185	for (int ils=-1;ils<=1;ils++) {
	186	int ildN = ild+ils; // ladders of neighbouring sensors
	187	if (ildN<0) ildN = fNLadders-1; else if (ildN==fNLadders) ildN = 0;
	188	for (int ids=-1;ids<=1;ids++) {
	189	int idtN = idt+ids;
	190	if (idtN<0 \|\| idtN==fNSensInLadder \|\| (ids==0&&ils==0)) continue;
	191	AliITSURecoSens* sensN = GetSensor(ildN,idtN); // potential neighbor
	192	int neighbID = ildN*fNSensInLadder+idtN;
	193	//
	194	int zType = 1; // side
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
	220	} // ladders
	221	//
a11ef2e4	222	}
	223
	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
a11ef2e4	234	int sensInLad = int(z*fSensDZInv);
0091e9f0	235	if (sensInLad<0) sensInLad = 0;
0091e9f0	236	else if (sensInLad>=fNSensInLadder) sensInLad = fNSensInLadder-1;
a11ef2e4	237	//
	238	double phi = impPar[0] - fPhiOffs;
	239	BringTo02Pi(phi);
	240	int ladID = int(phi*fDPhiLadInv); // ladder id
	241	int nsens = 0;
	242	//
	243	AliITSURecoSens* sensN,sens = GetSensor(ladIDfNSensInLadder+sensInLad);
	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	}
	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) {
	297	Int_t layS,ladS,sensS;
	298	for (int is=0;is<nFnd;is++) {
	299	fITSGeom->GetModuleId(fndSens[is],layS,ladS,sensS);
	300	printf("SNMC#%d(%d): %d %d %d \| ",is,mcLab,layS,ladS,sensS); GetSensorFromID(fndSens[is])->Print();
	301	}
	302	}
	303	}
	304
	305	//tmp<<<
	306
	307	double z = impPar[2];
	308	if (z>fZMax+impPar[3]) {
	309	if (nFnd>0) printf("MissedSens!!!\n");
	310	return 0; // outside of Z coverage
	311	}
	312	z -= fZMin;
	313	if (z<-impPar[3]) {
	314	if (nFnd>0) printf("MissedSens!!!\n");
	315	return 0; // outside of Z coverage
	316	}
	317	int sensInLad = int(z*fSensDZInv);
	318	if (sensInLad<0) sensInLad = 0;
	319	else if (sensInLad>=fNSensInLadder) sensInLad = fNSensInLadder-1;
	320	//
	321	double phi = impPar[0] - fPhiOffs;
	322	BringTo02Pi(phi);
	323	int ladID = int(phi*fDPhiLadInv); // ladder id
	324	int nsens = 0;
	325	//
	326	AliITSURecoSens* sensN,sens = GetSensor(ladIDfNSensInLadder+sensInLad);
	327	sensors[nsens++] = sens;
	328	//
	329	// check neighbours
	330	double zMn=impPar[2]-impPar[3], zMx=impPar[2]+impPar[3], phiMn=impPar[0]-impPar[1], phiMx=impPar[0]+impPar[1];
	331	BringTo02Pi(phiMn);
	332	BringTo02Pi(phiMx);
	333	//
	334	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbR)); // neighbor on the right (smaller phi)
	335	if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax())) sensors[nsens++] = sensN;
	336	//
	337	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTR)); // neighbor on the top right (smaller phi, larger Z)
	338	if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
	339	//
	340	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbT)); // neighbor on the top (larger Z)
	341	if (sensN && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
342	//
343	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTL)); // neighbor on the top left (larger Z, larger phi)
344	if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
345	//
346	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbL)); // neighbor on the left (larger phi)
347	if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin())) sensors[nsens++] = sensN;
348	//
349	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBL)); // neighbor on the bottom left (smaller Z, larger phi)
350	if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
351	//
352	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbB)); // neighbor on the bottom (smaller Z)
353	if (sensN && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
354	//
355	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBR)); // neighbor on the bottom right (smaller Z, smaller phi)
356	if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
357	//
358	if (mcLab>=0) {
359	Int_t layS,ladS,sensS;
360	printf("Found %d sensors for phi %.3f : %.3f \| Z %.4f %.4f\n", nsens,phiMn,phiMx,zMn,zMx);
361	for (int is=0;is<nsens;is++) {
362	fITSGeom->GetModuleId(sensors[is]->GetID()+fITSGeom->GetFirstModIndex(fActiveID),layS,ladS,sensS);
363	printf("*SNF#%d: %d %d %d \| ",is,layS,ladS,sensS); sensors[is]->Print();
364	}
365	for (int ism=0;ism<nFnd;ism++) {
366	AliITSURecoSens* snMC = GetSensorFromID(fndSens[ism]);
367	Bool_t ok=kFALSE;
368	for (int isf=0;isf<nsens;isf++) {
369	if (snMC==sensors[isf]) {ok=kTRUE;break;}
370	}
371	if (!ok) printf("MissedSens %d!!!\n",ism);
372	}
373	}
374	return nsens;
375	}
376	*/
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
	393	curSens = GetSensor(vID - fITSGeom->GetFirstModIndex(fActiveID));
	394	curSensID = vID;
	395	curSens->SetFirstClusterId(icl);
	396	}
	397	curSens->IncNClusters();
	398	}
	399	if (curSens) curSens->ProcessClusters(mode); // last sensor was not processed yet
	400	//
	401	}
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
	426	i -= fITSGeom->GetFirstModIndex(fActiveID);
	427	if (i<0\|\|i>=fNSensors) AliFatal(Form("Sensor with id=%d is not in layer %d",i+fITSGeom->GetFirstModIndex(fActiveID),fActiveID));
	428	return (AliITSURecoSens*)fSensors[i];
	429	}