[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSUClusterizer.cxx

#include <TTree.h>
#include <TFile.h>
#include <TObjArray.h>
#include <TParticle.h>
#include <TMath.h>
#include "AliRun.h"
#include "AliLog.h"
#include "AliMC.h"
#include <AliITSUSegmentationPix.h>
#include "AliITSUClusterizer.h"
#include "AliITSUGeomTGeo.h"
#include "AliITSUSegmentationPix.h"
#include "AliITSUDigitPix.h"
#include "AliITSURecoParam.h"
#include "AliITSUAux.h"
using namespace TMath;
using namespace AliITSUAux;

ClassImp(AliITSUClusterizer)

//______________________________________________________________________________
AliITSUClusterizer::AliITSUClusterizer(Int_t initNRow) 
:  fVolID(-1)
  ,fAllowDiagonalClusterization(kFALSE)
  ,fSegm(0)
  ,fRecoParam(0)
  ,fInputDigits(0)
  ,fInputDigitsReadIndex(0)
  ,fLayerID(0)
  ,fNLabels(0)
  ,fRawData(kFALSE)
  ,fLorAngCorrection(0)
  ,fOutputClusters(0)
  ,fDigitFreelist(0)
  ,fPartFreelist(0)
  ,fCandFreelist(0)
  ,fDigitFreelistBptrFirst(0)
  ,fDigitFreelistBptrLast(0)
  ,fPartFreelistBptr(0)
  ,fCandFreelistBptr(0)
#ifdef _ClusterTopology_
  ,fTopology(0)
  ,fMinCol(0)
  ,fMinRow(0)
#endif //_ClusterTopology_
{
  SetUniqueID(0);
  // c-tor
  SetNRow(initNRow);
//
#ifdef _ClusterTopology_
  AliInfo("*********************************************************************");
  AliInfo("ATTENTION: YOU ARE RUNNING IN SPECIAL MODE OF STORING CLUSTER PATTERN");
  AliInfo("*********************************************************************");
#endif //_ClusterTopology_

//
}

//______________________________________________________________________________
void AliITSUClusterizer::SetSegmentation(const AliITSUSegmentationPix *segm) 
{
  // attach segmentation, if needed, reinitialize array
  fSegm = segm;
  SetNRow(fSegm->GetNRow()); // reinitialize if needed

}

//______________________________________________________________________________
void AliITSUClusterizer::SetNRow(Int_t nr)
{
  // update buffers
  int nrOld = GetUniqueID();
  if (nrOld>=nr) return;
  SetUniqueID(nr);
  while (fDigitFreelistBptrFirst) {
    AliITSUClusterizerClusterDigit *next = fDigitFreelistBptrFirst[kDigitChunkSize-1].fNext;
    delete[] fDigitFreelistBptrFirst;
    fDigitFreelistBptrFirst=next;
  }
  delete[] fPartFreelistBptr;
  delete[] fCandFreelistBptr;
  //
  fPartFreelist=fPartFreelistBptr = new AliITSUClusterizerClusterPart[nr+1];
  fCandFreelist=fCandFreelistBptr = new AliITSUClusterizerClusterCand[nr+1];
  for (int i=0;i<nr;++i) {
    fPartFreelistBptr[i].fNextInRow = &fPartFreelistBptr[i+1];
    fCandFreelistBptr[i].fNext      = &fCandFreelistBptr[i+1];
  }  
}

//______________________________________________________________________________
AliITSUClusterizer::~AliITSUClusterizer() 
{
  // d-tor
  while (fDigitFreelistBptrFirst) {
    AliITSUClusterizerClusterDigit *next = fDigitFreelistBptrFirst[kDigitChunkSize-1].fNext;
    delete[] fDigitFreelistBptrFirst;
    fDigitFreelistBptrFirst=next;
  }
  delete[] fPartFreelistBptr;
  delete[] fCandFreelistBptr;
  //
}

//______________________________________________________________________________
AliITSUClusterizer::AliITSUClusterizerClusterDigit* AliITSUClusterizer::AllocDigitFreelist() 
{
  // allocate aux space
  AliITSUClusterizerClusterDigit *tmp = new AliITSUClusterizerClusterDigit[kDigitChunkSize];
  for (int i=0;i<kDigitChunkSize-2;++i) tmp[i].fNext=&tmp[i+1];
  tmp[kDigitChunkSize-2].fNext=0;
  tmp[kDigitChunkSize-1].fNext=0;
  if (!fDigitFreelistBptrFirst) fDigitFreelistBptrFirst=tmp;
  else                          fDigitFreelistBptrLast[kDigitChunkSize-1].fNext=tmp;
  fDigitFreelistBptrLast=tmp;
  return tmp;
}

//______________________________________________________________________________
AliITSUClusterizer::AliITSUClusterizerClusterDigit* AliITSUClusterizer::NextDigit() 
{
  // get next digit
  while (fInputDigitsReadIndex<fInputDigits->GetEntriesFast()) {
    AliITSUDigitPix *tmp=static_cast<AliITSUDigitPix*>(fInputDigits->UncheckedAt(fInputDigitsReadIndex++));
    if (TMath::Abs(tmp->GetROCycle())>=fRecoParam->GetMaxROCycle()) continue;
    AliITSUClusterizerClusterDigit *digit=AllocDigit();
    digit->fDigit=tmp;
    return digit;
  }
  return 0;
}

//______________________________________________________________________________
void AliITSUClusterizer::AttachPartToCand(AliITSUClusterizerClusterCand *cand,AliITSUClusterizerClusterPart *part) 
{
  // attach part
  part->fParent = cand;
  part->fPrevInCluster = 0;
  part->fNextInCluster = cand->fFirstPart;
  if (cand->fFirstPart) cand->fFirstPart->fPrevInCluster = part;
  cand->fFirstPart=part;
}

//______________________________________________________________________________
void AliITSUClusterizer::MergeCands(AliITSUClusterizerClusterCand *a,AliITSUClusterizerClusterCand *b) 
{
  // merge cluster parts
  AliITSUClusterizerClusterPart *ipart=b->fFirstPart; 
  AliITSUClusterizerClusterPart *jpart; 
  do {
    jpart=ipart;
    jpart->fParent=a;
  } while ((ipart=ipart->fNextInCluster));
  jpart->fNextInCluster=a->fFirstPart;
  jpart->fNextInCluster->fPrevInCluster=jpart;
  a->fFirstPart=b->fFirstPart;
  // merge digits
  b->fLastDigit->fNext=a->fFirstDigit;
  a->fFirstDigit=b->fFirstDigit;
  //  DeallocCand(b);
}

//______________________________________________________________________________
void AliITSUClusterizer::Transform(AliITSUClusterPix *cluster,AliITSUClusterizerClusterCand *cand) 
{
  // convert set of digits to cluster data in LOCAL frame
  const double k1to12 = 1./12;
  static int maxLbinDigit = AliITSdigit::GetNTracks();
  //
  fNLabels = 0;
  Int_t n=0;
  cand->fLastDigit->fNext=0;
  double x=0,z=0,xmn=1e9,xmx=-1e9,zmn=1e9,zmx=-1e9,px=0,pz=0;
  float  cx,cz;
  int charge=0;
  //
  for (AliITSUClusterizerClusterDigit *idigit=cand->fFirstDigit;idigit;idigit=idigit->fNext) {
    AliITSdigit* dig = idigit->fDigit;
    fSegm->DetToLocal(dig->GetCoord2(),dig->GetCoord1(),cx,cz); // center of pixel
    //
    // account for possible diod shift
    double ddx,ddz, dx=fSegm->Dpx(dig->GetCoord2()), dz=fSegm->Dpz(dig->GetCoord1());
    fSegm->GetDiodShift(dig->GetCoord2(),dig->GetCoord1(),ddx,ddz);
    //
    charge += dig->GetSignal();
    x += cx;
    z += cz;
    if (cx<xmn) xmn=cx;
    if (cx>xmx) xmx=cx;
    if (cz<zmn) zmn=cz;
    if (cz>zmx) zmx=cz;
    x += ddx*dx;
    z += ddz*dz;
    px += dx;
    pz += dz;
    //
    if (!fRawData) {
      for(Int_t dlab=0;dlab<maxLbinDigit;dlab++){
	Int_t digitlab = (dig->GetTracks())[dlab];
	if(digitlab<0) continue;
	AddLabel(digitlab);
      }
    }
    //
    ++n;
  }
  Int_t nx=1,nz=1;
  double dx = xmx-xmn, dz = zmx-zmn;
  if (n>1) {
    double fac=1./n;
    x  *= fac;  // mean coordinates
    z  *= fac;
    px *= fac;  // mean pitch
    pz *= fac; 
    nx = 1+Nint(dx/px);
    nz = 1+Nint(dz/pz);
  }
  x -= fLorAngCorrection;  // LorentzAngle correction
  cluster->SetX(x);
  cluster->SetZ(z);  
  cluster->SetY(0);
  cluster->SetSigmaZ2(nz>1 ? dz*dz*k1to12 : pz*pz*k1to12);
  cluster->SetSigmaY2(nx>1 ? dx*dx*k1to12 : px*px*k1to12);
  cluster->SetSigmaYZ(0);
  cluster->SetFrameLoc();
  if (n>AliITSUClusterPix::kMaskNPix) n = AliITSUClusterPix::kMaskNPix;
  cluster->SetNxNzN(nx>AliITSUClusterPix::kMaskNX ? AliITSUClusterPix::kMaskNX : nx,
		    nz>AliITSUClusterPix::kMaskNZ ? AliITSUClusterPix::kMaskNZ : nz,
		    n>AliITSUClusterPix::kMaskNPix ? AliITSUClusterPix::kMaskNPix : n);		    
  cluster->SetQ(charge); // note: this is MC info
  //  if (cluster->GetNPix()!=n || cluster->GetNx()!=nx || cluster->GetNz()!=nz) {
  //  printf("PROBLEM: %d %d %d -> %d %d %d\n",nx,nz,n,
  //	   cluster->GetNx(),cluster->GetNz(),cluster->GetNPix());
  //  cluster->Dump();
  //  }
  //
  if (!fRawData) {
    CheckLabels();
    int nl = Min(kMaxLabInCluster,fNLabels);
    for (int i=nl;i--;) cluster->SetLabel(fCurrLabels[i],i);
  }
  //
  // Set Volume id
  cluster->SetVolumeId(fVolID);
  //    printf("mod %d: (%.4lf,%.4lf)cm\n",fVolID,x,z);
  //
#ifdef _ClusterTopology_
  FillClusterTopology(cand,cluster);
#endif //_ClusterTopology_
  //
}

//______________________________________________________________________________
void AliITSUClusterizer::CloseCand(AliITSUClusterizerClusterCand *cand) 
{
  // finish cluster
  AliITSUClusterPix *cluster = (AliITSUClusterPix*)NextCluster();
  Transform(cluster,cand);
  DeallocDigits(cand->fFirstDigit,cand->fLastDigit);
  DeallocCand(cand);
}

//______________________________________________________________________________
void AliITSUClusterizer::ClosePart(AliITSUClusterizerClusterPart *part) 
{
  // finish cluster part
  AliITSUClusterizerClusterCand *cand=part->fParent;
  DetachPartFromCand(cand,part);
  DeallocPart(part);
  if (!cand->fFirstPart) CloseCand(cand); 
}

//______________________________________________________________________________
void AliITSUClusterizer::CloseRemainingParts(AliITSUClusterizerClusterPart *part) 
{
  // finish what is left
  while (part) {
    AliITSUClusterizerClusterPart *next=part->fNextInRow;
    ClosePart(part);
    part=next;
  } 
}

//______________________________________________________________________________
void AliITSUClusterizer::Clusterize() 
{
  // main algo for MC clustererization
  SetRawData(kFALSE);
  //
  AliITSUClusterizerClusterDigit *iDigit=NextDigit();
  AliITSUClusterizerClusterPart *iPrevRowBegin=0;
  AliITSUClusterizerClusterPart *iNextRowBegin=0;
  AliITSUClusterizerClusterPart *iPrevRow=0;
  AliITSUClusterizerClusterPart *iNextRow=0;
  Int_t lastV=0;
  //
  while (iDigit) {
    if (iDigit->fDigit->GetCoord2()!=lastV) {
      // NEW ROW
      if (iNextRow) iNextRow->fNextInRow=0;
      if (iPrevRowBegin) CloseRemainingParts(iPrevRowBegin);
      if (iDigit->fDigit->GetCoord2()==lastV+1) {
	iPrevRowBegin=iNextRowBegin;
	iPrevRow     =iNextRowBegin;
      }
      else {
	// there was an empty row
	CloseRemainingParts(iNextRowBegin);
	iPrevRowBegin=0;
	iPrevRow     =0;
      }
      iNextRowBegin=0;
      iNextRow     =0;
      lastV=iDigit->fDigit->GetCoord2(); 
    }
    // skip cluster parts before this digit
    int limCol = iDigit->fDigit->GetCoord1()-fAllowDiagonalClusterization;
    while (iPrevRow && iPrevRow->fUEnd < limCol) {
      iPrevRow=iPrevRow->fNextInRow;
    }
    // find the longest continous line of digits [iDigit,pDigit]=[iDigit,jDigit)
    AliITSUClusterizerClusterCand *cand=AllocCand(); 
    AliITSUClusterizerClusterDigit *pDigit=iDigit;
    AliITSUClusterizerClusterDigit *jDigit=NextDigit();
    cand->fFirstPart=0;
    cand->fFirstDigit=cand->fLastDigit=iDigit; // NB: first diggit is attached differently
    iDigit->fNext=0;
    Int_t lastU =iDigit->fDigit->GetCoord1();
    Int_t lastU1=lastU+1;
    while (jDigit && jDigit->fDigit->GetCoord1()==lastU1 && jDigit->fDigit->GetCoord2()==lastV) {
      pDigit=jDigit;
      jDigit=NextDigit();
      AttachDigitToCand(cand,pDigit);
      ++lastU1;
    }
    if (!fAllowDiagonalClusterization) --lastU1;
    AliITSUClusterizerClusterPart *part=AllocPart();
    part->fUBegin=lastU ;
    part->fUEnd  =lastU1;
    AttachPartToCand(cand,part);
    // merge all cluster candidates of the previous line touching this one,
    // advance to the last one, but keep that one the next active one
    AliITSUClusterizerClusterPart *jPrevRow=iPrevRow;
    while (jPrevRow && jPrevRow->fUBegin<=lastU1) {
      if (jPrevRow->fParent!=cand) {
	MergeCands(jPrevRow->fParent,cand);
	DeallocCand(cand);
	cand=jPrevRow->fParent;
      }
      iPrevRow=jPrevRow;
      jPrevRow=jPrevRow->fNextInRow;
    }
    if (iNextRow)
      iNextRow->fNextInRow=part;
    else
      iNextRowBegin=part;
    iNextRow=part;
    iDigit=jDigit;
  }
  // remove remaining cluster parts
  CloseRemainingParts(iPrevRowBegin);
  if (iNextRow) iNextRow->fNextInRow=0;
  CloseRemainingParts(iNextRowBegin);
  return;
}

//______________________________________________________________________________
void AliITSUClusterizer::PrepareLorentzAngleCorrection(Double_t bz)
{
  // calculate parameters for Lorentz Angle correction. Must be called 
  // after setting segmentation and recoparams
  fLorAngCorrection = 0.5*fRecoParam->GetTanLorentzAngle(fLayerID)*bz/kNominalBz*fSegm->Dy();
}

//______________________________________________________________________
void AliITSUClusterizer::CheckLabels() 
{
  // Tries to find mother's labels
  //
  if (fNLabels<1) return;
  AliRunLoader *rl = AliRunLoader::Instance();
  if(!rl) return;
  TTree *trK=(TTree*)rl->TreeK();
  if (!trK) return;
  //
  static int   labS[kMaxLabels];
  static float kine[kMaxLabels];
  Int_t nlabels = fNLabels; 
  Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
  for (Int_t i=fNLabels;i--;) labS[i] = fCurrLabels[i];
  //
  for (Int_t i=0;i<nlabels;i++) {
    Int_t label = labS[i];
    if (label>=ntracks) continue;
    TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
    kine[i] = part->Energy() - part->GetCalcMass(); // kinetic energy 
    if (kine[i] < 0.001) {    // reduce soft particles from the same cluster
      Int_t m=part->GetFirstMother();
      if (m<0) continue; // primary
      //
      if (part->GetStatusCode()>0) continue;
      //
      // if the parent is within the same cluster, assign parent's label
      for (int j=0;j<nlabels;j++) if (labS[j]==m) { labS[i] = m; break;}
    }
  } 
  //
  if (nlabels>kMaxLabInCluster) { // only 3 labels are stored in cluster, sort in decreasing momentum
    static int ind[kMaxLabels],labSS[kMaxLabels];
    TMath::Sort(nlabels,kine,ind);
    for (int i=nlabels;i--;) labSS[i] = labS[i];
    for (int i=nlabels;i--;) labS[i] = labSS[ind[i]]; 
  }
  //
  //compress labels -- if multi-times the same
  for (Int_t i=0;i<nlabels;i++) fCurrLabels[i]=-2;
  fNLabels = 0;
  int j=0;
  for (int i=0;i<nlabels;i++) {
    for (j=fNLabels;j--;) if (labS[i]==fCurrLabels[j]) break; // the label already there
    if (j<0) fCurrLabels[fNLabels++] = labS[i];
  }
  //
}

//------------------------------------------------------------------------

#ifdef _ClusterTopology_
//
//______________________________________________________________________
void AliITSUClusterizer::FillClusterTopology(const AliITSUClusterizerClusterCand *cand, AliITSUClusterPix* cl) const 
{
  // fill cluster topology bit pattern
  //
  //
  UShort_t mnCol=0xffff,mxCol=0,mnRow=0xffff,mxRow=0;
  for (AliITSUClusterizerClusterDigit *idigit=cand->fFirstDigit;idigit;idigit=idigit->fNext) {
    AliITSdigit* dig = idigit->fDigit;
    if (mnCol>dig->GetCoord1()) mnCol = (UShort_t)dig->GetCoord1();
    if (mxCol<dig->GetCoord1()) mxCol = (UShort_t)dig->GetCoord1();
    if (mnRow>dig->GetCoord2()) mnRow = (UShort_t)dig->GetCoord2();
    if (mxRow<dig->GetCoord2()) mxRow = (UShort_t)dig->GetCoord2();    
  }
  UShort_t nColSpan = UShort_t(mxCol-mnCol+1);
  UShort_t nRowSpan = UShort_t(mxRow-mnRow+1);
  UShort_t nColSpanW = nColSpan;
  UShort_t nRowSpanW = nRowSpan;
  if (nColSpan*nRowSpan>AliITSUClusterPix::kMaxPatternBits) { // need to store partial info
    // will crtail largest dimension
    if (nColSpan>nRowSpan) {
      if ( (nColSpanW=AliITSUClusterPix::kMaxPatternBits/nRowSpan)==0 ) {
	nColSpanW = 1;
	nRowSpanW = AliITSUClusterPix::kMaxPatternBits;
      }
    }
    else {
      if ( (nRowSpanW=AliITSUClusterPix::kMaxPatternBits/nColSpan)==0 ) {
	nRowSpanW = 1;
	nColSpanW = AliITSUClusterPix::kMaxPatternBits;
      }
    }
  }
  cl->ResetPattern();
  //
  cl->SetPatternRowSpan(nRowSpanW,nRowSpanW<nRowSpan);
  cl->SetPatternColSpan(nColSpanW,nColSpanW<nColSpan);
  cl->SetPatternMinRow(mnRow);
  cl->SetPatternMinCol(mnCol);
  //
  for (AliITSUClusterizerClusterDigit *idigit=cand->fFirstDigit;idigit;idigit=idigit->fNext) {
    AliITSdigit* dig = idigit->fDigit;
    UInt_t ir = dig->GetCoord2()-mnRow;
    UInt_t ic = dig->GetCoord1()-mnCol;
    if (ir<nRowSpanW && ic<nColSpanW) cl->SetPixel(ir,ic);
  }
  //
}

#endif //_ClusterTopology_
Commit	Line	Data
5e375bb4	1	#include <TTree.h>
6cae87c5	2	#include <TFile.h>
5e375bb4	3	#include <TObjArray.h>
4bac12be	4	#include <TParticle.h>
5e375bb4	5	#include <TMath.h>
4bac12be	6	#include "AliRun.h"
6f54924d	7	#include "AliLog.h"
4bac12be	8	#include "AliMC.h"
5e375bb4	9	#include <AliITSUSegmentationPix.h>
b69620f8	10	#include "AliITSUClusterizer.h"
5e375bb4	11	#include "AliITSUGeomTGeo.h"
5e375bb4	12	#include "AliITSUSegmentationPix.h"
c7811d12	13	#include "AliITSUDigitPix.h"
5e375bb4	14	#include "AliITSURecoParam.h"
889b1493	15	#include "AliITSUAux.h"
5e375bb4	16	using namespace TMath;
889b1493	17	using namespace AliITSUAux;
b69620f8	18
	19	ClassImp(AliITSUClusterizer)
	20
	21	//______________________________________________________________________________
	22	AliITSUClusterizer::AliITSUClusterizer(Int_t initNRow)
	23	: fVolID(-1)
ee58ce21	24	,fAllowDiagonalClusterization(kFALSE)
b69620f8	25	,fSegm(0)
5e375bb4	26	,fRecoParam(0)
b69620f8	27	,fInputDigits(0)
b69620f8	28	,fInputDigitsReadIndex(0)
889b1493	29	,fLayerID(0)
4bac12be	30	,fNLabels(0)
4bac12be	31	,fRawData(kFALSE)
889b1493	32	,fLorAngCorrection(0)
b69620f8	33	,fOutputClusters(0)
	34	,fDigitFreelist(0)
	35	,fPartFreelist(0)
	36	,fCandFreelist(0)
	37	,fDigitFreelistBptrFirst(0)
	38	,fDigitFreelistBptrLast(0)
	39	,fPartFreelistBptr(0)
	40	,fCandFreelistBptr(0)
6cae87c5	41	#ifdef _ClusterTopology_
6cae87c5	42	,fTopology(0)
	43	,fMinCol(0)
	44	,fMinRow(0)
	45	#endif //_ClusterTopology_
b69620f8	46	{
	47	SetUniqueID(0);
	48	// c-tor
	49	SetNRow(initNRow);
6f54924d	50	//
	51	#ifdef _ClusterTopology_
	52	AliInfo("*********************************************************************");
	53	AliInfo("ATTENTION: YOU ARE RUNNING IN SPECIAL MODE OF STORING CLUSTER PATTERN");
	54	AliInfo("*********************************************************************");
	55	#endif //_ClusterTopology_
	56
	57	//
b69620f8	58	}
	59
	60	//______________________________________________________________________________
	61	void AliITSUClusterizer::SetSegmentation(const AliITSUSegmentationPix *segm)
	62	{
	63	// attach segmentation, if needed, reinitialize array
	64	fSegm = segm;
	65	SetNRow(fSegm->GetNRow()); // reinitialize if needed
	66
	67	}
	68
	69	//______________________________________________________________________________
	70	void AliITSUClusterizer::SetNRow(Int_t nr)
	71	{
	72	// update buffers
	73	int nrOld = GetUniqueID();
	74	if (nrOld>=nr) return;
	75	SetUniqueID(nr);
	76	while (fDigitFreelistBptrFirst) {
	77	AliITSUClusterizerClusterDigit *next = fDigitFreelistBptrFirst[kDigitChunkSize-1].fNext;
	78	delete[] fDigitFreelistBptrFirst;
	79	fDigitFreelistBptrFirst=next;
	80	}
	81	delete[] fPartFreelistBptr;
	82	delete[] fCandFreelistBptr;
	83	//
	84	fPartFreelist=fPartFreelistBptr = new AliITSUClusterizerClusterPart[nr+1];
	85	fCandFreelist=fCandFreelistBptr = new AliITSUClusterizerClusterCand[nr+1];
	86	for (int i=0;i<nr;++i) {
	87	fPartFreelistBptr[i].fNextInRow = &fPartFreelistBptr[i+1];
	88	fCandFreelistBptr[i].fNext = &fCandFreelistBptr[i+1];
	89	}
	90	}
	91
	92	//______________________________________________________________________________
	93	AliITSUClusterizer::~AliITSUClusterizer()
	94	{
	95	// d-tor
	96	while (fDigitFreelistBptrFirst) {
	97	AliITSUClusterizerClusterDigit *next = fDigitFreelistBptrFirst[kDigitChunkSize-1].fNext;
	98	delete[] fDigitFreelistBptrFirst;
	99	fDigitFreelistBptrFirst=next;
	100	}
	101	delete[] fPartFreelistBptr;
	102	delete[] fCandFreelistBptr;
6cae87c5	103	//
b69620f8	104	}
	105
	106	//______________________________________________________________________________
	107	AliITSUClusterizer::AliITSUClusterizerClusterDigit* AliITSUClusterizer::AllocDigitFreelist()
	108	{
	109	// allocate aux space
	110	AliITSUClusterizerClusterDigit *tmp = new AliITSUClusterizerClusterDigit[kDigitChunkSize];
	111	for (int i=0;i<kDigitChunkSize-2;++i) tmp[i].fNext=&tmp[i+1];
	112	tmp[kDigitChunkSize-2].fNext=0;
	113	tmp[kDigitChunkSize-1].fNext=0;
	114	if (!fDigitFreelistBptrFirst) fDigitFreelistBptrFirst=tmp;
	115	else fDigitFreelistBptrLast[kDigitChunkSize-1].fNext=tmp;
	116	fDigitFreelistBptrLast=tmp;
	117	return tmp;
	118	}
	119
	120	//______________________________________________________________________________
	121	AliITSUClusterizer::AliITSUClusterizerClusterDigit* AliITSUClusterizer::NextDigit()
	122	{
	123	// get next digit
c7811d12	124	while (fInputDigitsReadIndex<fInputDigits->GetEntriesFast()) {
	125	AliITSUDigitPix tmp=static_cast<AliITSUDigitPix>(fInputDigits->UncheckedAt(fInputDigitsReadIndex++));
	126	if (TMath::Abs(tmp->GetROCycle())>=fRecoParam->GetMaxROCycle()) continue;
b69620f8	127	AliITSUClusterizerClusterDigit *digit=AllocDigit();
b69620f8	128	digit->fDigit=tmp;
b69620f8	129	return digit;
b69620f8	130	}
c7811d12	131	return 0;
b69620f8	132	}
	133
	134	//______________________________________________________________________________
	135	void AliITSUClusterizer::AttachPartToCand(AliITSUClusterizerClusterCand cand,AliITSUClusterizerClusterPart part)
	136	{
	137	// attach part
	138	part->fParent = cand;
	139	part->fPrevInCluster = 0;
	140	part->fNextInCluster = cand->fFirstPart;
	141	if (cand->fFirstPart) cand->fFirstPart->fPrevInCluster = part;
	142	cand->fFirstPart=part;
	143	}
	144
	145	//______________________________________________________________________________
	146	void AliITSUClusterizer::MergeCands(AliITSUClusterizerClusterCand a,AliITSUClusterizerClusterCand b)
	147	{
	148	// merge cluster parts
	149	AliITSUClusterizerClusterPart *ipart=b->fFirstPart;
	150	AliITSUClusterizerClusterPart *jpart;
	151	do {
	152	jpart=ipart;
	153	jpart->fParent=a;
	154	} while ((ipart=ipart->fNextInCluster));
	155	jpart->fNextInCluster=a->fFirstPart;
	156	jpart->fNextInCluster->fPrevInCluster=jpart;
	157	a->fFirstPart=b->fFirstPart;
	158	// merge digits
	159	b->fLastDigit->fNext=a->fFirstDigit;
	160	a->fFirstDigit=b->fFirstDigit;
	161	// DeallocCand(b);
	162	}
	163
	164	//______________________________________________________________________________
5e375bb4	165	void AliITSUClusterizer::Transform(AliITSUClusterPix cluster,AliITSUClusterizerClusterCand cand)
b69620f8	166	{
5e375bb4	167	// convert set of digits to cluster data in LOCAL frame
5e375bb4	168	const double k1to12 = 1./12;
4bac12be	169	static int maxLbinDigit = AliITSdigit::GetNTracks();
5e375bb4	170	//
4bac12be	171	fNLabels = 0;
b69620f8	172	Int_t n=0;
b69620f8	173	cand->fLastDigit->fNext=0;
5e375bb4	174	double x=0,z=0,xmn=1e9,xmx=-1e9,zmn=1e9,zmx=-1e9,px=0,pz=0;
5e375bb4	175	float cx,cz;
29ad4146	176	int charge=0;
6cae87c5	177	//
b69620f8	178	for (AliITSUClusterizerClusterDigit *idigit=cand->fFirstDigit;idigit;idigit=idigit->fNext) {
4bac12be	179	AliITSdigit* dig = idigit->fDigit;
ee58ce21	180	fSegm->DetToLocal(dig->GetCoord2(),dig->GetCoord1(),cx,cz); // center of pixel
	181	//
	182	// account for possible diod shift
	183	double ddx,ddz, dx=fSegm->Dpx(dig->GetCoord2()), dz=fSegm->Dpz(dig->GetCoord1());
	184	fSegm->GetDiodShift(dig->GetCoord2(),dig->GetCoord1(),ddx,ddz);
	185	//
29ad4146	186	charge += dig->GetSignal();
5e375bb4	187	x += cx;
	188	z += cz;
	189	if (cx<xmn) xmn=cx;
	190	if (cx>xmx) xmx=cx;
	191	if (cz<zmn) zmn=cz;
	192	if (cz>zmx) zmx=cz;
ee58ce21	193	x += ddx*dx;
	194	z += ddz*dz;
	195	px += dx;
	196	pz += dz;
4bac12be	197	//
	198	if (!fRawData) {
	199	for(Int_t dlab=0;dlab<maxLbinDigit;dlab++){
	200	Int_t digitlab = (dig->GetTracks())[dlab];
	201	if(digitlab<0) continue;
	202	AddLabel(digitlab);
	203	}
	204	}
	205	//
b69620f8	206	++n;
b69620f8	207	}
b7e5fec5	208	Int_t nx=1,nz=1;
5e375bb4	209	double dx = xmx-xmn, dz = zmx-zmn;
	210	if (n>1) {
	211	double fac=1./n;
	212	x *= fac; // mean coordinates
	213	z *= fac;
	214	px *= fac; // mean pitch
	215	pz *= fac;
	216	nx = 1+Nint(dx/px);
	217	nz = 1+Nint(dz/pz);
b69620f8	218	}
889b1493	219	x -= fLorAngCorrection; // LorentzAngle correction
b69620f8	220	cluster->SetX(x);
b7e5fec5	221	cluster->SetZ(z);
5e375bb4	222	cluster->SetY(0);
c8d1f258	223	cluster->SetSigmaZ2(nz>1 ? dzdzk1to12 : pzpzk1to12);
c8d1f258	224	cluster->SetSigmaY2(nx>1 ? dxdxk1to12 : pxpxk1to12);
5e375bb4	225	cluster->SetSigmaYZ(0);
5e375bb4	226	cluster->SetFrameLoc();
b7e5fec5	227	if (n>AliITSUClusterPix::kMaskNPix) n = AliITSUClusterPix::kMaskNPix;
	228	cluster->SetNxNzN(nx>AliITSUClusterPix::kMaskNX ? AliITSUClusterPix::kMaskNX : nx,
	229	nz>AliITSUClusterPix::kMaskNZ ? AliITSUClusterPix::kMaskNZ : nz,
	230	n>AliITSUClusterPix::kMaskNPix ? AliITSUClusterPix::kMaskNPix : n);
29ad4146	231	cluster->SetQ(charge); // note: this is MC info
b7e5fec5	232	// if (cluster->GetNPix()!=n \|\| cluster->GetNx()!=nx \|\| cluster->GetNz()!=nz) {
	233	// printf("PROBLEM: %d %d %d -> %d %d %d\n",nx,nz,n,
	234	// cluster->GetNx(),cluster->GetNz(),cluster->GetNPix());
	235	// cluster->Dump();
	236	// }
5e375bb4	237	//
4bac12be	238	if (!fRawData) {
	239	CheckLabels();
	240	int nl = Min(kMaxLabInCluster,fNLabels);
	241	for (int i=nl;i--;) cluster->SetLabel(fCurrLabels[i],i);
	242	}
	243	//
b69620f8	244	// Set Volume id
	245	cluster->SetVolumeId(fVolID);
	246	// printf("mod %d: (%.4lf,%.4lf)cm\n",fVolID,x,z);
6cae87c5	247	//
6cae87c5	248	#ifdef _ClusterTopology_
6f54924d	249	FillClusterTopology(cand,cluster);
6cae87c5	250	#endif //_ClusterTopology_
6cae87c5	251	//
b69620f8	252	}
	253
	254	//______________________________________________________________________________
	255	void AliITSUClusterizer::CloseCand(AliITSUClusterizerClusterCand *cand)
	256	{
	257	// finish cluster
5e375bb4	258	AliITSUClusterPix cluster = (AliITSUClusterPix)NextCluster();
b69620f8	259	Transform(cluster,cand);
	260	DeallocDigits(cand->fFirstDigit,cand->fLastDigit);
	261	DeallocCand(cand);
	262	}
	263
	264	//______________________________________________________________________________
	265	void AliITSUClusterizer::ClosePart(AliITSUClusterizerClusterPart *part)
	266	{
	267	// finish cluster part
	268	AliITSUClusterizerClusterCand *cand=part->fParent;
	269	DetachPartFromCand(cand,part);
	270	DeallocPart(part);
	271	if (!cand->fFirstPart) CloseCand(cand);
	272	}
	273
	274	//______________________________________________________________________________
	275	void AliITSUClusterizer::CloseRemainingParts(AliITSUClusterizerClusterPart *part)
	276	{
	277	// finish what is left
	278	while (part) {
	279	AliITSUClusterizerClusterPart *next=part->fNextInRow;
	280	ClosePart(part);
	281	part=next;
	282	}
	283	}
	284
	285	//______________________________________________________________________________
	286	void AliITSUClusterizer::Clusterize()
	287	{
4bac12be	288	// main algo for MC clustererization
	289	SetRawData(kFALSE);
	290	//
b69620f8	291	AliITSUClusterizerClusterDigit *iDigit=NextDigit();
	292	AliITSUClusterizerClusterPart *iPrevRowBegin=0;
	293	AliITSUClusterizerClusterPart *iNextRowBegin=0;
	294	AliITSUClusterizerClusterPart *iPrevRow=0;
	295	AliITSUClusterizerClusterPart *iNextRow=0;
	296	Int_t lastV=0;
6cae87c5	297	//
b69620f8	298	while (iDigit) {
4bac12be	299	if (iDigit->fDigit->GetCoord2()!=lastV) {
b69620f8	300	// NEW ROW
	301	if (iNextRow) iNextRow->fNextInRow=0;
	302	if (iPrevRowBegin) CloseRemainingParts(iPrevRowBegin);
4bac12be	303	if (iDigit->fDigit->GetCoord2()==lastV+1) {
b69620f8	304	iPrevRowBegin=iNextRowBegin;
	305	iPrevRow =iNextRowBegin;
	306	}
	307	else {
	308	// there was an empty row
	309	CloseRemainingParts(iNextRowBegin);
	310	iPrevRowBegin=0;
	311	iPrevRow =0;
	312	}
	313	iNextRowBegin=0;
	314	iNextRow =0;
4bac12be	315	lastV=iDigit->fDigit->GetCoord2();
b69620f8	316	}
b69620f8	317	// skip cluster parts before this digit
ee58ce21	318	int limCol = iDigit->fDigit->GetCoord1()-fAllowDiagonalClusterization;
ee58ce21	319	while (iPrevRow && iPrevRow->fUEnd < limCol) {
b69620f8	320	iPrevRow=iPrevRow->fNextInRow;
	321	}
	322	// find the longest continous line of digits [iDigit,pDigit]=[iDigit,jDigit)
	323	AliITSUClusterizerClusterCand *cand=AllocCand();
	324	AliITSUClusterizerClusterDigit *pDigit=iDigit;
	325	AliITSUClusterizerClusterDigit *jDigit=NextDigit();
	326	cand->fFirstPart=0;
	327	cand->fFirstDigit=cand->fLastDigit=iDigit; // NB: first diggit is attached differently
	328	iDigit->fNext=0;
4bac12be	329	Int_t lastU =iDigit->fDigit->GetCoord1();
b69620f8	330	Int_t lastU1=lastU+1;
4bac12be	331	while (jDigit && jDigit->fDigit->GetCoord1()==lastU1 && jDigit->fDigit->GetCoord2()==lastV) {
b69620f8	332	pDigit=jDigit;
	333	jDigit=NextDigit();
	334	AttachDigitToCand(cand,pDigit);
	335	++lastU1;
	336	}
ee58ce21	337	if (!fAllowDiagonalClusterization) --lastU1;
b69620f8	338	AliITSUClusterizerClusterPart *part=AllocPart();
	339	part->fUBegin=lastU ;
	340	part->fUEnd =lastU1;
	341	AttachPartToCand(cand,part);
	342	// merge all cluster candidates of the previous line touching this one,
	343	// advance to the last one, but keep that one the next active one
	344	AliITSUClusterizerClusterPart *jPrevRow=iPrevRow;
	345	while (jPrevRow && jPrevRow->fUBegin<=lastU1) {
	346	if (jPrevRow->fParent!=cand) {
	347	MergeCands(jPrevRow->fParent,cand);
	348	DeallocCand(cand);
	349	cand=jPrevRow->fParent;
	350	}
	351	iPrevRow=jPrevRow;
	352	jPrevRow=jPrevRow->fNextInRow;
	353	}
	354	if (iNextRow)
	355	iNextRow->fNextInRow=part;
	356	else
	357	iNextRowBegin=part;
	358	iNextRow=part;
	359	iDigit=jDigit;
	360	}
	361	// remove remaining cluster parts
	362	CloseRemainingParts(iPrevRowBegin);
	363	if (iNextRow) iNextRow->fNextInRow=0;
	364	CloseRemainingParts(iNextRowBegin);
	365	return;
	366	}
889b1493	367
	368	//______________________________________________________________________________
	369	void AliITSUClusterizer::PrepareLorentzAngleCorrection(Double_t bz)
	370	{
	371	// calculate parameters for Lorentz Angle correction. Must be called
	372	// after setting segmentation and recoparams
	373	fLorAngCorrection = 0.5fRecoParam->GetTanLorentzAngle(fLayerID)bz/kNominalBz*fSegm->Dy();
	374	}
4bac12be	375
	376	//______________________________________________________________________
	377	void AliITSUClusterizer::CheckLabels()
	378	{
	379	// Tries to find mother's labels
	380	//
	381	if (fNLabels<1) return;
	382	AliRunLoader *rl = AliRunLoader::Instance();
	383	if(!rl) return;
	384	TTree trK=(TTree)rl->TreeK();
	385	if (!trK) return;
	386	//
	387	static int labS[kMaxLabels];
	388	static float kine[kMaxLabels];
	389	Int_t nlabels = fNLabels;
	390	Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
	391	for (Int_t i=fNLabels;i--;) labS[i] = fCurrLabels[i];
	392	//
	393	for (Int_t i=0;i<nlabels;i++) {
	394	Int_t label = labS[i];
	395	if (label>=ntracks) continue;
	396	TParticle part=(TParticle)gAlice->GetMCApp()->Particle(label);
	397	kine[i] = part->Energy() - part->GetCalcMass(); // kinetic energy
81a0c667	398	if (kine[i] < 0.001) { // reduce soft particles from the same cluster
4bac12be	399	Int_t m=part->GetFirstMother();
	400	if (m<0) continue; // primary
	401	//
	402	if (part->GetStatusCode()>0) continue;
	403	//
	404	// if the parent is within the same cluster, assign parent's label
	405	for (int j=0;j<nlabels;j++) if (labS[j]==m) { labS[i] = m; break;}
	406	}
	407	}
	408	//
	409	if (nlabels>kMaxLabInCluster) { // only 3 labels are stored in cluster, sort in decreasing momentum
	410	static int ind[kMaxLabels],labSS[kMaxLabels];
	411	TMath::Sort(nlabels,kine,ind);
	412	for (int i=nlabels;i--;) labSS[i] = labS[i];
	413	for (int i=nlabels;i--;) labS[i] = labSS[ind[i]];
	414	}
	415	//
	416	//compress labels -- if multi-times the same
	417	for (Int_t i=0;i<nlabels;i++) fCurrLabels[i]=-2;
	418	fNLabels = 0;
	419	int j=0;
	420	for (int i=0;i<nlabels;i++) {
	421	for (j=fNLabels;j--;) if (labS[i]==fCurrLabels[j]) break; // the label already there
	422	if (j<0) fCurrLabels[fNLabels++] = labS[i];
	423	}
	424	//
	425	}
6cae87c5	426
	427	//------------------------------------------------------------------------
	428
	429	#ifdef _ClusterTopology_
	430	//
	431	//______________________________________________________________________
6f54924d	432	void AliITSUClusterizer::FillClusterTopology(const AliITSUClusterizerClusterCand cand, AliITSUClusterPix cl) const
6cae87c5	433	{
6f54924d	434	// fill cluster topology bit pattern
6cae87c5	435	//
6cae87c5	436	//
6f54924d	437	UShort_t mnCol=0xffff,mxCol=0,mnRow=0xffff,mxRow=0;
6cae87c5	438	for (AliITSUClusterizerClusterDigit *idigit=cand->fFirstDigit;idigit;idigit=idigit->fNext) {
6cae87c5	439	AliITSdigit* dig = idigit->fDigit;
6f54924d	440	if (mnCol>dig->GetCoord1()) mnCol = (UShort_t)dig->GetCoord1();
	441	if (mxCol<dig->GetCoord1()) mxCol = (UShort_t)dig->GetCoord1();
	442	if (mnRow>dig->GetCoord2()) mnRow = (UShort_t)dig->GetCoord2();
	443	if (mxRow<dig->GetCoord2()) mxRow = (UShort_t)dig->GetCoord2();
6cae87c5	444	}
6f54924d	445	UShort_t nColSpan = UShort_t(mxCol-mnCol+1);
	446	UShort_t nRowSpan = UShort_t(mxRow-mnRow+1);
	447	UShort_t nColSpanW = nColSpan;
	448	UShort_t nRowSpanW = nRowSpan;
	449	if (nColSpan*nRowSpan>AliITSUClusterPix::kMaxPatternBits) { // need to store partial info
	450	// will crtail largest dimension
	451	if (nColSpan>nRowSpan) {
	452	if ( (nColSpanW=AliITSUClusterPix::kMaxPatternBits/nRowSpan)==0 ) {
	453	nColSpanW = 1;
	454	nRowSpanW = AliITSUClusterPix::kMaxPatternBits;
	455	}
	456	}
	457	else {
	458	if ( (nRowSpanW=AliITSUClusterPix::kMaxPatternBits/nColSpan)==0 ) {
	459	nRowSpanW = 1;
	460	nColSpanW = AliITSUClusterPix::kMaxPatternBits;
	461	}
6cae87c5	462	}
6f54924d	463	}
	464	cl->ResetPattern();
	465	//
	466	cl->SetPatternRowSpan(nRowSpanW,nRowSpanW<nRowSpan);
	467	cl->SetPatternColSpan(nColSpanW,nColSpanW<nColSpan);
	468	cl->SetPatternMinRow(mnRow);
	469	cl->SetPatternMinCol(mnCol);
	470	//
	471	for (AliITSUClusterizerClusterDigit *idigit=cand->fFirstDigit;idigit;idigit=idigit->fNext) {
	472	AliITSdigit* dig = idigit->fDigit;
	473	UInt_t ir = dig->GetCoord2()-mnRow;
	474	UInt_t ic = dig->GetCoord1()-mnCol;
	475	if (ir<nRowSpanW && ic<nColSpanW) cl->SetPixel(ir,ic);
6cae87c5	476	}
6cae87c5	477	//
6cae87c5	478	}
	479
	480	#endif //_ClusterTopology_