[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 %d 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[kNNeighbors])
{
  // find sensors having intersection with track
  // impPar contains: lab phi of track, dphi, labZ, dz
  //
  return 0;
  /*
  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 sensInRow = int(z*fSensDZInv);
  if      (sensInRow<0) sensInRow = 0;
  else if (sensInRow>=fNSensorsPerRow) sensInRow = fNSensorsPerRow-1;
  //
  double phi = impPar[0] - fPhiOffs;
  BringTo02Pi(phi);
  int rowID = int(phi*fSensDPhiInv);  // stave id
  //
  int sensID = rowID*fNSensorsPerRow + sensInRow; // most probable candidate
  int nsens = 0;
  //
  AliITSURecoSens* sensN,*sens = GetSensor(sesnID);
  //
  // make sure this is best matching sensor
  if (sens->GetZMin()<impPar[2] && sens->GetZMax()>impPar[2] && 
      OKforPhiMin(sens->GetPhiMin(),impPar[0]) && OKforPhiMax(sens->GetPhiMax(),impPar[0]) ) 

  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*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   = 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)
	17	:fActiveID(-1)
	18	,fNSensors(0)
5009207f	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)
5009207f	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)
5009207f	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)
5009207f	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
5009207f	68	delete fSensors;
5009207f	69	delete[] fSensVIDtoMatrixID;
ee52e7b5	70	if (GetOwnsClusterArray()) delete fClusters;
a11ef2e4	71	}
	72
	73	//______________________________________________________
	74	void AliITSURecoLayer::Print(Option_t* opt) const
	75	{
	76	//print
5009207f	77	printf("Lr %-15s %d (act:%+d), NSens: %4d in %d rows\| MaxStep:%.2f ",
	78	GetName(),GetID(),GetActiveID(),GetNSensors(),GetNSensorRows(),fMaxStep);
	79	printf("%6.3f<R<%6.3f \| %+8.3f<Z<%+8.3f dZ:%6.3f dPhi:%6.3f\n",fRMin,fRMax,fZMin,fZMax,
	80	fSensDZInv>0 ? 1/fSensDZInv : 0, fSensDPhiInv>0 ? 1/fSensDPhiInv : 0);
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;
5009207f	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];
5009207f	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;
5009207f	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);
5009207f	132	if (phiMin>1e8) phiMin=phi;
a11ef2e4	133	else if (!OKforPhiMin(phiMin,phi)) phiMin=phi;
5009207f	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	}
5009207f	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	//______________________________________________________
5009207f	209	Int_t AliITSURecoLayer::FindSensors(const double* impPar, AliITSURecoSens *sensors[kNNeighbors])
a11ef2e4	210	{
	211	// find sensors having intersection with track
	212	// impPar contains: lab phi of track, dphi, labZ, dz
5785a9d8	213	//
5009207f	214	return 0;
5009207f	215	/*
a11ef2e4	216	double z = impPar[2];
0091e9f0	217	if (z>fZMax+impPar[3]) return 0; // outside of Z coverage
a11ef2e4	218	z -= fZMin;
0091e9f0	219	if (z<-impPar[3]) return 0; // outside of Z coverage
5009207f	220	int sensInRow = int(z*fSensDZInv);
	221	if (sensInRow<0) sensInRow = 0;
	222	else if (sensInRow>=fNSensorsPerRow) sensInRow = fNSensorsPerRow-1;
a11ef2e4	223	//
	224	double phi = impPar[0] - fPhiOffs;
	225	BringTo02Pi(phi);
5009207f	226	int rowID = int(phi*fSensDPhiInv); // stave id
	227	//
	228	int sensID = rowID*fNSensorsPerRow + sensInRow; // most probable candidate
a11ef2e4	229	int nsens = 0;
a11ef2e4	230	//
5009207f	231	AliITSURecoSens* sensN,*sens = GetSensor(sesnID);
	232	//
	233	// make sure this is best matching sensor
	234	if (sens->GetZMin()<impPar[2] && sens->GetZMax()>impPar[2] &&
	235	OKforPhiMin(sens->GetPhiMin(),impPar[0]) && OKforPhiMax(sens->GetPhiMax(),impPar[0]) )
	236
a11ef2e4	237	sensors[nsens++] = sens;
5785a9d8	238	//
a11ef2e4	239	// check neighbours
	240	double zMn=impPar[2]-impPar[3], zMx=impPar[2]+impPar[3], phiMn=impPar[0]-impPar[1], phiMx=impPar[0]+impPar[1];
	241	BringTo02Pi(phiMn);
	242	BringTo02Pi(phiMx);
	243	//
173b3073	244	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbR)); // neighbor on the right (smaller phi)
173b3073	245	if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax())) sensors[nsens++] = sensN;
a11ef2e4	246	//
173b3073	247	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTR)); // neighbor on the top right (smaller phi, larger Z)
173b3073	248	if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
a11ef2e4	249	//
173b3073	250	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbT)); // neighbor on the top (larger Z)
173b3073	251	if (sensN && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
a11ef2e4	252	//
173b3073	253	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTL)); // neighbor on the top left (larger Z, larger phi)
173b3073	254	if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
a11ef2e4	255	//
173b3073	256	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbL)); // neighbor on the left (larger phi)
173b3073	257	if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin())) sensors[nsens++] = sensN;
a11ef2e4	258	//
173b3073	259	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBL)); // neighbor on the bottom left (smaller Z, larger phi)
173b3073	260	if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
a11ef2e4	261	//
173b3073	262	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbB)); // neighbor on the bottom (smaller Z)
a11ef2e4	263	if (sensN && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
a11ef2e4	264	//
173b3073	265	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBR)); // neighbor on the bottom right (smaller Z, smaller phi)
173b3073	266	if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
a11ef2e4	267	//
a11ef2e4	268	return nsens;
5009207f	269	*/
a11ef2e4	270	}
5009207f	271
c03e4f8a	272	//*/
5785a9d8	273	/*
	274	Int_t AliITSURecoLayer::FindSensors(const double* impPar, AliITSURecoSens *sensors[AliITSURecoSens::kNNeighbors],int mcLab)
	275	{
	276	// find sensors having intersection with track
	277	// impPar contains: lab phi of track, dphi, labZ, dz
	278
	279	//tmp>>>
	280	int nFnd = 0;
	281	int fndSens[50];
	282	if (mcLab>=0) { // find correct sensors from MC info
	283	int ncl = GetNClusters();
	284	for (int icl=ncl;icl--;) {
	285	AliCluster* cl = GetCluster(icl);
	286	for (int ilb=0;ilb<3;ilb++) {
	287	if (cl->GetLabel(ilb)<0) break;
	288	if (cl->GetLabel(ilb)==mcLab) {fndSens[nFnd++] = cl->GetVolumeId(); break;}
	289	}
	290	}
	291	if (nFnd>0) {
852af72e	292	Int_t layS,staS,sensS;
5785a9d8	293	for (int is=0;is<nFnd;is++) {
852af72e	294	fITSGeom->GetChipId(fndSens[is],layS,staS,sensS);
852af72e	295	printf("SNMC#%d(%d): %d %d %d \| ",is,mcLab,layS,staS,sensS); GetSensorFromID(fndSens[is])->Print();
5785a9d8	296	}
	297	}
	298	}
	299
	300	//tmp<<<
	301
	302	double z = impPar[2];
	303	if (z>fZMax+impPar[3]) {
	304	if (nFnd>0) printf("MissedSens!!!\n");
	305	return 0; // outside of Z coverage
	306	}
	307	z -= fZMin;
	308	if (z<-impPar[3]) {
	309	if (nFnd>0) printf("MissedSens!!!\n");
	310	return 0; // outside of Z coverage
	311	}
852af72e	312	int sensInSta = int(z*fSensDZInv);
	313	if (sensInSta<0) sensInSta = 0;
	314	else if (sensInSta>=fNSensInStave) sensInSta = fNSensInStave-1;
5785a9d8	315	//
	316	double phi = impPar[0] - fPhiOffs;
	317	BringTo02Pi(phi);
5009207f	318	int staID = int(phi*fSensDPhiInv); // stave id
5785a9d8	319	int nsens = 0;
5785a9d8	320	//
852af72e	321	AliITSURecoSens* sensN,sens = GetSensor(staIDfNSensInStave+sensInSta);
5785a9d8	322	sensors[nsens++] = sens;
	323	//
	324	// check neighbours
	325	double zMn=impPar[2]-impPar[3], zMx=impPar[2]+impPar[3], phiMn=impPar[0]-impPar[1], phiMx=impPar[0]+impPar[1];
	326	BringTo02Pi(phiMn);
	327	BringTo02Pi(phiMx);
	328	//
	329	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbR)); // neighbor on the right (smaller phi)
	330	if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax())) sensors[nsens++] = sensN;
	331	//
	332	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTR)); // neighbor on the top right (smaller phi, larger Z)
	333	if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
	334	//
	335	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbT)); // neighbor on the top (larger Z)
	336	if (sensN && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
	337	//
	338	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbTL)); // neighbor on the top left (larger Z, larger phi)
	339	if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMin()<zMx) sensors[nsens++] = sensN;
	340	//
	341	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbL)); // neighbor on the left (larger phi)
	342	if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin())) sensors[nsens++] = sensN;
	343	//
	344	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBL)); // neighbor on the bottom left (smaller Z, larger phi)
	345	if (sensN && OKforPhiMax(phiMx,sensN->GetPhiMin()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
	346	//
	347	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbB)); // neighbor on the bottom (smaller Z)
	348	if (sensN && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
	349	//
	350	sensN = GetSensor(sens->GetNeighborID(AliITSURecoSens::kNghbBR)); // neighbor on the bottom right (smaller Z, smaller phi)
	351	if (sensN && OKforPhiMin(phiMn,sensN->GetPhiMax()) && sensN->GetZMax()>zMn) sensors[nsens++] = sensN;
	352	//
	353	if (mcLab>=0) {
852af72e	354	Int_t layS,staS,sensS;
5785a9d8	355	printf("Found %d sensors for phi %.3f : %.3f \| Z %.4f %.4f\n", nsens,phiMn,phiMx,zMn,zMx);
5785a9d8	356	for (int is=0;is<nsens;is++) {
852af72e	357	fITSGeom->GetChipId(sensors[is]->GetID()+fITSGeom->GetFirstModIndex(fActiveID),layS,staS,sensS);
852af72e	358	printf("*SNF#%d: %d %d %d \| ",is,layS,staS,sensS); sensors[is]->Print();
5785a9d8	359	}
	360	for (int ism=0;ism<nFnd;ism++) {
	361	AliITSURecoSens* snMC = GetSensorFromID(fndSens[ism]);
	362	Bool_t ok=kFALSE;
	363	for (int isf=0;isf<nsens;isf++) {
	364	if (snMC==sensors[isf]) {ok=kTRUE;break;}
	365	}
	366	if (!ok) printf("MissedSens %d!!!\n",ism);
	367	}
	368	}
	369	return nsens;
	370	}
	371	*/
32d38de2	372	//______________________________________________________
	373	void AliITSURecoLayer::ProcessClusters(Int_t mode)
	374	{
	375	// register in each sensor of the layer its cluster.
	376	// the clusters of the layer must be sorted per sensor
	377	int ncl = fClusters->GetEntriesFast();
	378	int curSensID = -1;
0091e9f0	379	for (int i=fNSensors;i--;) GetSensor(i)->SetNClusters(0);
32d38de2	380	AliITSURecoSens* curSens = 0;
	381	for (int icl=0;icl<ncl;icl++) {
	382	AliITSUClusterPix* cl = (AliITSUClusterPix*) fClusters->UncheckedAt(icl);
	383	cl->GoToFrameTrk();
	384	int vID = cl->GetVolumeId();
	385	if (vID<curSensID) {AliFatal("Clusters are not sorted in increasing sensorID");}
	386	if (vID>curSensID) {
	387	if (curSens) curSens->ProcessClusters(mode); // prepare clusters for reconstruction
852af72e	388	curSens = GetSensor(vID - fITSGeom->GetFirstChipIndex(fActiveID));
32d38de2	389	curSensID = vID;
	390	curSens->SetFirstClusterId(icl);
	391	}
	392	curSens->IncNClusters();
	393	}
	394	if (curSens) curSens->ProcessClusters(mode); // last sensor was not processed yet
	395	//
	396	}
dde91d5d	397
	398	//______________________________________________________
	399	Bool_t AliITSURecoLayer::IsEqual(const TObject* obj) const
	400	{
	401	// check if layers are equal in R
44785f3e	402	const AliITSURecoLayer* lr = (const AliITSURecoLayer*)obj;
dde91d5d	403	return Abs(lr->GetR()-GetR())<1e-6 ? kTRUE : kFALSE;
	404	}
	405
	406	//______________________________________________________
	407	Int_t AliITSURecoLayer::Compare(const TObject* obj) const
	408	{
	409	// compare two layers
44785f3e	410	const AliITSURecoLayer* lr = (const AliITSURecoLayer*)obj;
dde91d5d	411	double dr = GetR() - lr->GetR();
	412	if (Abs(dr)<1e-6) return 0;
	413	return dr>0 ? 1:-1;
	414	//
	415	}
8b16dbae	416
	417	//_________________________________________________________________
	418	AliITSURecoSens* AliITSURecoLayer::GetSensorFromID(Int_t i) const
	419	{
	420	// get sensor from its global id
852af72e	421	i -= fITSGeom->GetFirstChipIndex(fActiveID);
852af72e	422	if (i<0\|\|i>=fNSensors) AliFatal(Form("Sensor with id=%d is not in layer %d",i+fITSGeom->GetFirstChipIndex(fActiveID),fActiveID));
5009207f	423	return GetSensor(SensVIDtoMatrixID(i));
8b16dbae	424	}