#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)
{
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)
{
SetNameTitle(name,name);
Build();
}
AliITSURecoLayer::~AliITSURecoLayer()
{
delete fSensors;
delete[] fSensVIDtoMatrixID;
if (GetOwnsClusterArray()) delete fClusters;
}
void AliITSURecoLayer::Print(Option_t* opt) const
{
printf("Lr %-15s %d (act:%+d), NSens: %4d in %3d rows| MaxStep:%.2f ",
GetName(),GetID(),GetActiveID(),GetNSensors(),GetNSensorRows(),fMaxStep);
printf("%6.3f<R<%6.3f | %+8.3f<Z<%+8.3f dZ:%6.3f dPhi:%6.3f\n",fRMin,fRMax,fZMin,fZMax,
fSensDZInv>0 ? 1/fSensDZInv : 0, fSensDPhiInv>0 ? 1/fSensDPhiInv : 0);
TString opts = opt; opts.ToLower();
if (opts.Contains("sn")) for (int i=0;i<GetNSensors();i++) GetSensor(i)->Print(opt);
}
void AliITSURecoLayer::Build()
{
const double kSafeR = 0.05;
if (fActiveID<0) return;
int nStaves=fITSGeom->GetNStaves(fActiveID);
fNSensorRows = nStaves;
if (fITSGeom->GetNHalfStaves(fActiveID)>0) fNSensorRows *= fITSGeom->GetNHalfStaves(fActiveID);
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;
fITSGeom->GetOrigMatrix(sID,mmod);
for (int ix=0;ix<2;ix++) {
loc[0] = (ix-0.5)*kSegm->Dx();
for (int iy=0;iy<2;iy++) {
loc[1] = (iy-0.5)*kSegm->Dy();
for (int iz=0;iz<2;iz++) {
loc[2] = (iz-0.5)*kSegm->Dz();
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();
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();
for (int iy=0;iy<2;iy++) {
loc[1] = (iy-0.5)*kSegm->Dy();
for (int iz=0;iz<2;iz++) {
loc[2] = (iz-0.5)*kSegm->Dz();
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]);
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;
fSensDZInv = fNSensorsPerRow/(fZMax-fZMin);
fSensDPhiInv = fNSensorRows/(2*Pi());
fRMin = Sqrt(fRMin);
fRMax = Sqrt(fRMax);
fR = 0.5*(fRMin+fRMax);
fRMin -= kSafeR;
fRMax += kSafeR;
}
Int_t AliITSURecoLayer::FindSensors(const double* impPar, AliITSURecoSens *sensors[kMaxSensMatching])
{
double zMn=impPar[2]-impPar[3], zMx=impPar[2]+impPar[3];
if (zMn>fZMax) return 0;
if (zMx<fZMin) return 0;
int zCenID = int((impPar[2]-fZMin)*fSensDZInv);
if (zCenID<0) zCenID = 0;
else if (zCenID>=fNSensorsPerRow) zCenID = fNSensorsPerRow-1;
double phiCn = impPar[0] - fPhiOffs;
int rowCenID = int(phiCn*fSensDPhiInv);
int res = 0;
AliITSURecoSens* sensPrev=0, *sens = GetSensor(rowCenID,zCenID);
while ( (res=sens->CheckCoverage(impPar[0], impPar[2])) ) {
if (res&AliITSURecoSens::kRight) {if (++rowCenID==fNSensorRows) rowCenID=0;}
else if (res&AliITSURecoSens::kLeft) {if (--rowCenID<0) rowCenID = fNSensorRows-1;}
if (res&AliITSURecoSens::kUp) {if (++zCenID==fNSensorsPerRow) zCenID = fNSensorsPerRow-1;}
else if (res&AliITSURecoSens::kDown) {if (--zCenID<0) zCenID = 0;}
AliITSURecoSens* sensAlt = GetSensor(rowCenID,zCenID);
if (sensAlt==sens || sensAlt==sensPrev) break;
sensPrev = sens;
sens = sensAlt;
}
int nFnd = 0;
sensors[nFnd++] = sens;
double phiMn = impPar[0]-impPar[1], phiMx = impPar[0]+impPar[1];
BringTo02Pi(phiMn);
BringTo02Pi(phiMx);
const int kNNeighb = 8;
const int kCheckNeighb[2][kNNeighb] = {
{ 1, 1, 0,-1,-1,-1, 0, 1},
{ 0, 1, 1, 1, 0,-1,-1,-1}
};
for (int inb=kNNeighb;inb--;) {
int idz = kCheckNeighb[1][inb];
int iz = zCenID + idz;
if (iz<0 || iz>=fNSensorsPerRow) continue;
int idphi = kCheckNeighb[0][inb];
int iphi = rowCenID + idphi;
if (iphi<0) iphi += fNSensorRows;
else if (iphi>=fNSensorRows) iphi -= fNSensorRows;
sens = GetSensor(iphi,iz);
if (idz>0) {if (zMx<sens->GetZMin()) continue;}
else if (idz<0) {if (zMn>sens->GetZMax()) continue;}
if (idphi>0) {if (!OKforPhiMin(sens->GetPhiMin(),phiMx)) continue;}
else if (idphi<0) {if (!OKforPhiMax(sens->GetPhiMax(),phiMn)) continue;}
sensors[nFnd++] = sens;
if (nFnd==kMaxSensMatching) break;
}
return nFnd;
}
void AliITSURecoLayer::ProcessClusters(Int_t mode)
{
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->SetRecoInfo(0);
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);
curSens = GetSensorFromID(vID);
curSensID = vID;
curSens->SetFirstClusterId(icl);
}
curSens->IncNClusters();
}
if (curSens) curSens->ProcessClusters(mode);
}
Bool_t AliITSURecoLayer::IsEqual(const TObject* obj) const
{
const AliITSURecoLayer* lr = (const AliITSURecoLayer*)obj;
return Abs(lr->GetR()-GetR())<1e-6 ? kTRUE : kFALSE;
}
Int_t AliITSURecoLayer::Compare(const TObject* obj) const
{
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
{
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));
}