[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 "AliITSdigit.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
  if (fInputDigitsReadIndex<fInputDigits->GetEntriesFast()) {
    AliITSdigit *tmp=static_cast<AliITSdigit*>(fInputDigits->UncheckedAt(fInputDigitsReadIndex++));
    AliITSUClusterizerClusterDigit *digit=AllocDigit();
    digit->fDigit=tmp;
    return digit;
  }
  else
    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;
  }
  UChar_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();
  cluster->SetNxNzN(nx,nz,n);
  cluster->SetQ(charge); // note: this is MC info
  //
  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.02) {    // 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"
	12	#include "AliITSUSegmentationPix.h"
	13	#include "AliITSdigit.h"
	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
	124	if (fInputDigitsReadIndex<fInputDigits->GetEntriesFast()) {
	125	AliITSdigit tmp=static_cast<AliITSdigit>(fInputDigits->UncheckedAt(fInputDigitsReadIndex++));
	126	AliITSUClusterizerClusterDigit *digit=AllocDigit();
	127	digit->fDigit=tmp;
b69620f8	128	return digit;
	129	}
	130	else
	131	return 0;
	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	}
5e375bb4	208	UChar_t nx=1,nz=1;
	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);
b69620f8	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();
29ad4146	227	cluster->SetNxNzN(nx,nz,n);
29ad4146	228	cluster->SetQ(charge); // note: this is MC info
5e375bb4	229	//
4bac12be	230	if (!fRawData) {
	231	CheckLabels();
	232	int nl = Min(kMaxLabInCluster,fNLabels);
	233	for (int i=nl;i--;) cluster->SetLabel(fCurrLabels[i],i);
	234	}
	235	//
b69620f8	236	// Set Volume id
	237	cluster->SetVolumeId(fVolID);
	238	// printf("mod %d: (%.4lf,%.4lf)cm\n",fVolID,x,z);
6cae87c5	239	//
6cae87c5	240	#ifdef _ClusterTopology_
6f54924d	241	FillClusterTopology(cand,cluster);
6cae87c5	242	#endif //_ClusterTopology_
6cae87c5	243	//
b69620f8	244	}
	245
	246	//______________________________________________________________________________
	247	void AliITSUClusterizer::CloseCand(AliITSUClusterizerClusterCand *cand)
	248	{
	249	// finish cluster
5e375bb4	250	AliITSUClusterPix cluster = (AliITSUClusterPix)NextCluster();
b69620f8	251	Transform(cluster,cand);
	252	DeallocDigits(cand->fFirstDigit,cand->fLastDigit);
	253	DeallocCand(cand);
	254	}
	255
	256	//______________________________________________________________________________
	257	void AliITSUClusterizer::ClosePart(AliITSUClusterizerClusterPart *part)
	258	{
	259	// finish cluster part
	260	AliITSUClusterizerClusterCand *cand=part->fParent;
	261	DetachPartFromCand(cand,part);
	262	DeallocPart(part);
	263	if (!cand->fFirstPart) CloseCand(cand);
	264	}
	265
	266	//______________________________________________________________________________
	267	void AliITSUClusterizer::CloseRemainingParts(AliITSUClusterizerClusterPart *part)
	268	{
	269	// finish what is left
	270	while (part) {
	271	AliITSUClusterizerClusterPart *next=part->fNextInRow;
	272	ClosePart(part);
	273	part=next;
	274	}
	275	}
	276
	277	//______________________________________________________________________________
	278	void AliITSUClusterizer::Clusterize()
	279	{
4bac12be	280	// main algo for MC clustererization
	281	SetRawData(kFALSE);
	282	//
b69620f8	283	AliITSUClusterizerClusterDigit *iDigit=NextDigit();
	284	AliITSUClusterizerClusterPart *iPrevRowBegin=0;
	285	AliITSUClusterizerClusterPart *iNextRowBegin=0;
	286	AliITSUClusterizerClusterPart *iPrevRow=0;
	287	AliITSUClusterizerClusterPart *iNextRow=0;
	288	Int_t lastV=0;
6cae87c5	289	//
b69620f8	290	while (iDigit) {
4bac12be	291	if (iDigit->fDigit->GetCoord2()!=lastV) {
b69620f8	292	// NEW ROW
	293	if (iNextRow) iNextRow->fNextInRow=0;
	294	if (iPrevRowBegin) CloseRemainingParts(iPrevRowBegin);
4bac12be	295	if (iDigit->fDigit->GetCoord2()==lastV+1) {
b69620f8	296	iPrevRowBegin=iNextRowBegin;
	297	iPrevRow =iNextRowBegin;
	298	}
	299	else {
	300	// there was an empty row
	301	CloseRemainingParts(iNextRowBegin);
	302	iPrevRowBegin=0;
	303	iPrevRow =0;
	304	}
	305	iNextRowBegin=0;
	306	iNextRow =0;
4bac12be	307	lastV=iDigit->fDigit->GetCoord2();
b69620f8	308	}
b69620f8	309	// skip cluster parts before this digit
ee58ce21	310	int limCol = iDigit->fDigit->GetCoord1()-fAllowDiagonalClusterization;
ee58ce21	311	while (iPrevRow && iPrevRow->fUEnd < limCol) {
b69620f8	312	iPrevRow=iPrevRow->fNextInRow;
	313	}
	314	// find the longest continous line of digits [iDigit,pDigit]=[iDigit,jDigit)
	315	AliITSUClusterizerClusterCand *cand=AllocCand();
	316	AliITSUClusterizerClusterDigit *pDigit=iDigit;
	317	AliITSUClusterizerClusterDigit *jDigit=NextDigit();
	318	cand->fFirstPart=0;
	319	cand->fFirstDigit=cand->fLastDigit=iDigit; // NB: first diggit is attached differently
	320	iDigit->fNext=0;
4bac12be	321	Int_t lastU =iDigit->fDigit->GetCoord1();
b69620f8	322	Int_t lastU1=lastU+1;
4bac12be	323	while (jDigit && jDigit->fDigit->GetCoord1()==lastU1 && jDigit->fDigit->GetCoord2()==lastV) {
b69620f8	324	pDigit=jDigit;
	325	jDigit=NextDigit();
	326	AttachDigitToCand(cand,pDigit);
	327	++lastU1;
	328	}
ee58ce21	329	if (!fAllowDiagonalClusterization) --lastU1;
b69620f8	330	AliITSUClusterizerClusterPart *part=AllocPart();
	331	part->fUBegin=lastU ;
	332	part->fUEnd =lastU1;
	333	AttachPartToCand(cand,part);
	334	// merge all cluster candidates of the previous line touching this one,
	335	// advance to the last one, but keep that one the next active one
	336	AliITSUClusterizerClusterPart *jPrevRow=iPrevRow;
	337	while (jPrevRow && jPrevRow->fUBegin<=lastU1) {
	338	if (jPrevRow->fParent!=cand) {
	339	MergeCands(jPrevRow->fParent,cand);
	340	DeallocCand(cand);
	341	cand=jPrevRow->fParent;
	342	}
	343	iPrevRow=jPrevRow;
	344	jPrevRow=jPrevRow->fNextInRow;
	345	}
	346	if (iNextRow)
	347	iNextRow->fNextInRow=part;
	348	else
	349	iNextRowBegin=part;
	350	iNextRow=part;
	351	iDigit=jDigit;
	352	}
	353	// remove remaining cluster parts
	354	CloseRemainingParts(iPrevRowBegin);
	355	if (iNextRow) iNextRow->fNextInRow=0;
	356	CloseRemainingParts(iNextRowBegin);
	357	return;
	358	}
889b1493	359
	360	//______________________________________________________________________________
	361	void AliITSUClusterizer::PrepareLorentzAngleCorrection(Double_t bz)
	362	{
	363	// calculate parameters for Lorentz Angle correction. Must be called
	364	// after setting segmentation and recoparams
	365	fLorAngCorrection = 0.5fRecoParam->GetTanLorentzAngle(fLayerID)bz/kNominalBz*fSegm->Dy();
	366	}
4bac12be	367
	368	//______________________________________________________________________
	369	void AliITSUClusterizer::CheckLabels()
	370	{
	371	// Tries to find mother's labels
	372	//
	373	if (fNLabels<1) return;
	374	AliRunLoader *rl = AliRunLoader::Instance();
	375	if(!rl) return;
	376	TTree trK=(TTree)rl->TreeK();
	377	if (!trK) return;
	378	//
	379	static int labS[kMaxLabels];
	380	static float kine[kMaxLabels];
	381	Int_t nlabels = fNLabels;
	382	Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
	383	for (Int_t i=fNLabels;i--;) labS[i] = fCurrLabels[i];
	384	//
	385	for (Int_t i=0;i<nlabels;i++) {
	386	Int_t label = labS[i];
	387	if (label>=ntracks) continue;
	388	TParticle part=(TParticle)gAlice->GetMCApp()->Particle(label);
	389	kine[i] = part->Energy() - part->GetCalcMass(); // kinetic energy
	390	if (kine[i] < 0.02) { // reduce soft particles from the same cluster
	391	Int_t m=part->GetFirstMother();
	392	if (m<0) continue; // primary
	393	//
	394	if (part->GetStatusCode()>0) continue;
	395	//
	396	// if the parent is within the same cluster, assign parent's label
	397	for (int j=0;j<nlabels;j++) if (labS[j]==m) { labS[i] = m; break;}
	398	}
	399	}
	400	//
	401	if (nlabels>kMaxLabInCluster) { // only 3 labels are stored in cluster, sort in decreasing momentum
	402	static int ind[kMaxLabels],labSS[kMaxLabels];
	403	TMath::Sort(nlabels,kine,ind);
	404	for (int i=nlabels;i--;) labSS[i] = labS[i];
	405	for (int i=nlabels;i--;) labS[i] = labSS[ind[i]];
	406	}
	407	//
	408	//compress labels -- if multi-times the same
	409	for (Int_t i=0;i<nlabels;i++) fCurrLabels[i]=-2;
	410	fNLabels = 0;
	411	int j=0;
	412	for (int i=0;i<nlabels;i++) {
	413	for (j=fNLabels;j--;) if (labS[i]==fCurrLabels[j]) break; // the label already there
	414	if (j<0) fCurrLabels[fNLabels++] = labS[i];
	415	}
	416	//
	417	}
6cae87c5	418
	419	//------------------------------------------------------------------------
	420
	421	#ifdef _ClusterTopology_
	422	//
	423	//______________________________________________________________________
6f54924d	424	void AliITSUClusterizer::FillClusterTopology(const AliITSUClusterizerClusterCand cand, AliITSUClusterPix cl) const
6cae87c5	425	{
6f54924d	426	// fill cluster topology bit pattern
6cae87c5	427	//
6cae87c5	428	//
6f54924d	429	UShort_t mnCol=0xffff,mxCol=0,mnRow=0xffff,mxRow=0;
6cae87c5	430	for (AliITSUClusterizerClusterDigit *idigit=cand->fFirstDigit;idigit;idigit=idigit->fNext) {
6cae87c5	431	AliITSdigit* dig = idigit->fDigit;
6f54924d	432	if (mnCol>dig->GetCoord1()) mnCol = (UShort_t)dig->GetCoord1();
	433	if (mxCol<dig->GetCoord1()) mxCol = (UShort_t)dig->GetCoord1();
	434	if (mnRow>dig->GetCoord2()) mnRow = (UShort_t)dig->GetCoord2();
	435	if (mxRow<dig->GetCoord2()) mxRow = (UShort_t)dig->GetCoord2();
6cae87c5	436	}
6f54924d	437	UShort_t nColSpan = UShort_t(mxCol-mnCol+1);
	438	UShort_t nRowSpan = UShort_t(mxRow-mnRow+1);
	439	UShort_t nColSpanW = nColSpan;
	440	UShort_t nRowSpanW = nRowSpan;
	441	if (nColSpan*nRowSpan>AliITSUClusterPix::kMaxPatternBits) { // need to store partial info
	442	// will crtail largest dimension
	443	if (nColSpan>nRowSpan) {
	444	if ( (nColSpanW=AliITSUClusterPix::kMaxPatternBits/nRowSpan)==0 ) {
	445	nColSpanW = 1;
	446	nRowSpanW = AliITSUClusterPix::kMaxPatternBits;
	447	}
	448	}
	449	else {
	450	if ( (nRowSpanW=AliITSUClusterPix::kMaxPatternBits/nColSpan)==0 ) {
	451	nRowSpanW = 1;
	452	nColSpanW = AliITSUClusterPix::kMaxPatternBits;
	453	}
6cae87c5	454	}
6f54924d	455	}
	456	cl->ResetPattern();
	457	//
	458	cl->SetPatternRowSpan(nRowSpanW,nRowSpanW<nRowSpan);
	459	cl->SetPatternColSpan(nColSpanW,nColSpanW<nColSpan);
	460	cl->SetPatternMinRow(mnRow);
	461	cl->SetPatternMinCol(mnCol);
	462	//
	463	for (AliITSUClusterizerClusterDigit *idigit=cand->fFirstDigit;idigit;idigit=idigit->fNext) {
	464	AliITSdigit* dig = idigit->fDigit;
	465	UInt_t ir = dig->GetCoord2()-mnRow;
	466	UInt_t ic = dig->GetCoord1()-mnCol;
	467	if (ir<nRowSpanW && ic<nColSpanW) cl->SetPixel(ir,ic);
6cae87c5	468	}
6cae87c5	469	//
6cae87c5	470	}
	471
	472	#endif //_ClusterTopology_