#include <TGeoVolume.h>
#include <TGeoTube.h>
#include <TGeoManager.h>
#include <TTree.h>
#include "AliITSURecoDet.h"
#include "AliITSUGeomTGeo.h"
#include "AliITSsegmentation.h"
#include "AliITSUClusterPix.h"
#include "AliITSUReconstructor.h"
ClassImp(AliITSURecoDet)
const Char_t* AliITSURecoDet::fgkBeamPipeVolName = "IP_PIPE";
AliITSURecoDet::AliITSURecoDet()
: fNLayers(0)
,fNLayersActive(0)
,fRMax(-1)
,fRMin(-1)
,fRITSTPCRef(-1)
,fLayers(0)
,fLayersActive(0)
,fGeom(0)
{
}
AliITSURecoDet::AliITSURecoDet(AliITSUGeomTGeo* geom, const char* name)
: fNLayers(0)
,fNLayersActive(0)
,fRMax(-1)
,fRMin(-1)
,fRITSTPCRef(-1)
,fLayers(0)
,fLayersActive(0)
,fGeom(geom)
{
SetNameTitle(name,name);
fLayers.SetOwner(kTRUE);
fLayersActive.SetOwner(kFALSE);
Build();
}
AliITSURecoDet::~AliITSURecoDet()
{
fLayersActive.Clear();
fLayers.Clear();
}
void AliITSURecoDet::Print(Option_t* opt) const
{
printf("Detector %s, %d layers, %d active layers\n",GetName(),GetNLayers(),GetNLayersActive());
TString opts = opt; opts.ToLower();
if (opts.Contains("lr")) {
for (int i=0;i<GetNLayers();i++) GetLayer(i)->Print(opt);
printf("ITS-TPC matching reference R: %.3f\n",fRITSTPCRef);
}
}
void AliITSURecoDet::AddLayer(const AliITSURecoLayer* lr)
{
fLayers.AddLast((TObject*)lr);
fNLayers++;
if (lr->IsActive()) {
fLayersActive.AddLast((TObject*)lr);
fNLayersActive++;
}
}
Bool_t AliITSURecoDet::Build()
{
if (!fGeom) AliFatal("Geometry interface is not set");
int nlr = fGeom->GetNLayers();
if (!nlr) AliFatal("No geometry loaded");
for (int ilr=0;ilr<nlr;ilr++) {
int lrTyp = fGeom->GetLayerChipTypeID(ilr);
AliITSURecoLayer* lra = new AliITSURecoLayer(Form("Lr%d%s%d",ilr,fGeom->GetChipTypeName(lrTyp),
lrTyp%AliITSUGeomTGeo::kMaxSegmPerChipType),
ilr,fGeom);
lra->SetPassive(kFALSE);
AddLayer(lra);
}
double rMin,rMax,zMin,zMax;
TGeoVolume *v = gGeoManager->GetVolume(fgkBeamPipeVolName);
AliITSURecoLayer* lrp = 0;
if (!v) AliWarning("No beam pipe found in geometry");
else {
TGeoTube *t=(TGeoTube*)v->GetShape();
rMin = t->GetRmin();
rMax = t->GetRmax();
zMin =-t->GetDz();
zMax = t->GetDz();
lrp = new AliITSURecoLayer("BeamPipe");
lrp->SetRMin(rMin);
lrp->SetRMax(rMax);
lrp->SetR(0.5*(rMin+rMax));
lrp->SetZMin(zMin);
lrp->SetZMax(zMax);
lrp->SetPassive(kTRUE);
AddLayer(lrp);
}
const AliITSURecoParam* recopar = AliITSUReconstructor::GetRecoParam();
if (recopar) {
lrp = new AliITSURecoLayer("TPC-ITSwall");
lrp->SetRMin(AliITSUReconstructor::GetRecoParam()->GetTPCITSWallRMin());
lrp->SetRMax(AliITSUReconstructor::GetRecoParam()->GetTPCITSWallRMax());
lrp->SetR(0.5*(lrp->GetRMin()+lrp->GetRMax()));
lrp->SetZMin(-AliITSUReconstructor::GetRecoParam()->GetTPCITSWallZSpanH());
lrp->SetZMax( AliITSUReconstructor::GetRecoParam()->GetTPCITSWallZSpanH());
lrp->SetMaxStep( AliITSUReconstructor::GetRecoParam()->GetTPCITSWallMaxStep());
lrp->SetPassive(kTRUE);
AddLayer(lrp);
}
else {
AliWarning("RecoParam is not available, TPC-ITS wall is not set");
}
IndexLayers();
Print("lr");
return kTRUE;
}
void AliITSURecoDet::IndexLayers()
{
const Double_t kRMargin = 1e-2;
fLayersActive.Sort();
for (int i=0;i<fNLayersActive;i++) GetLayerActive(i)->SetActiveID(i);
fLayers.Sort();
for (int i=0;i<fNLayers;i++) GetLayer(i)->SetID(i);
if (fNLayers>0) {
SetRMin(GetLayer(0)->GetRMin()-kRMargin);
SetRMax(GetLayer(fNLayers-1)->GetRMax()+kRMargin);
}
const double kOffsLastActR = 5.;
int lastActive = GetLrIDActive(fNLayersActive-1);
AliITSURecoLayer* lrA = GetLayer(lastActive);
double rref = lrA->GetRMax() + kOffsLastActR;
if (lastActive < fNLayers-1) {
AliITSURecoLayer* lrL = GetLayer(lastActive+1);
if (lrL->GetRMin()<=rref) rref = lrL->GetRMin();
if (rref - lrA->GetRMax()<kOffsLastActR) {
AliError(Form("The ITS-TPC matching reference R=%.2f is too close to last active R=%.3f",rref,lrA->GetRMax()));
}
}
SetRITSTPCRef(rref);
}
Int_t AliITSURecoDet::FindLastLayerID(Double_t r, int dir) const
{
int ilr;
if (dir>0) {
for (ilr=0;ilr<fNLayers;ilr++) {
AliITSURecoLayer* lr = GetLayer(ilr);
if ( r<lr->GetR(-dir) ) break;
}
return --ilr;
}
else {
for (ilr=fNLayers;ilr--;) {
AliITSURecoLayer* lr = GetLayer(ilr);
if ( r>lr->GetR(-dir) ) break;
}
ilr++;
return ilr<fNLayers ? ilr:-1;
}
}
Int_t AliITSURecoDet::FindFirstLayerID(Double_t r, int dir) const
{
int ilr;
if (dir>0) {
for (ilr=0;ilr<fNLayers;ilr++) {
AliITSURecoLayer* lr = GetLayer(ilr);
if ( r<lr->GetR(dir) ) break;
}
return ilr<fNLayers ? ilr:-1;
}
else {
for (ilr=fNLayers;ilr--;) {
AliITSURecoLayer* lr = GetLayer(ilr);
if ( r>lr->GetR(dir) ) break;
}
return ilr;
}
}
void AliITSURecoDet::CreateClusterArrays()
{
for (int ilr=0;ilr<fNLayersActive;ilr++) {
AliITSURecoLayer* lr = GetLayerActive(ilr);
lr->SetOwnsClusterArray(kTRUE);
int tpDet = fGeom->GetLayerChipTypeID(ilr)/AliITSUGeomTGeo::kMaxSegmPerChipType;
if (tpDet == AliITSUGeomTGeo::kChipTypePix) {
lr->SetClusters(new TClonesArray(AliITSUClusterPix::Class()));
}
else {
AliFatal(Form("Unknown detector type %d",tpDet));
}
}
}
Int_t AliITSURecoDet::LoadClusters(TTree* treeRP)
{
if (!treeRP) return 0;
for (int ilr=fNLayersActive;ilr--;) {
TBranch* br = treeRP->GetBranch(Form("ITSRecPoints%d",ilr));
if (!br) AliFatal(Form("Provided cluster tree does not contain branch for layer %d",ilr));
br->SetAddress( GetLayerActive(ilr)->GetClustersAddress() );
}
return treeRP->GetEntry(0);
}
Int_t FindIndex(AliITSUClusterPix **parr, Int_t i, AliITSUClusterPix *c) {
Int_t b=0;
AliITSUClusterPix *cc=parr[b];
if (c->Compare(cc) <= 0) return b;
Int_t e=b+i-1;
cc=parr[e];
if (c->Compare(cc) > 0) return e+1;
Int_t m=(b+e)/2;
for (; b<e; m=(b+e)/2) {
cc=parr[m];
if (c->Compare(cc) > 0) b=m+1;
else e=m;
}
return m;
}
void AliITSURecoDet::SortClusters(AliITSUClusterPix::SortMode_t mode)
{
AliITSUClusterPix::SetSortMode( mode );
for (int ilr=fNLayersActive;ilr--;) {
TClonesArray *arr=GetLayerActive(ilr)->GetClusters();
const Int_t kNcl=arr->GetEntriesFast();
AliITSUClusterPix **parr=new AliITSUClusterPix*[kNcl];
parr[0]=(AliITSUClusterPix *)arr->UncheckedAt(0);
for (Int_t j=1; j<kNcl; j++) {
AliITSUClusterPix *c=(AliITSUClusterPix*)arr->UncheckedAt(j);
Int_t i=FindIndex(parr,j,c);
Int_t k=j-i;
memmove(parr+i+1 ,parr+i,k*sizeof(AliITSUClusterPix*));
parr[i]=c;
}
TObject **cont=arr->GetObjectRef();
for (Int_t i=0; i<kNcl; i++) cont[i]=parr[i];
delete[] parr;
}
}