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

#include <TTree.h>
#include <TObjArray.h>
#include <TParticle.h>
#include <TMath.h>
#include "AliRun.h"
#include "AliMC.h"
#include <AliITSUSegmentationPix.h>
#include "AliITSUClusterizer.h"
#include "AliITSUClusterPix.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)
{
  SetUniqueID(0);
  // c-tor
  SetNRow(initNRow);
}

//______________________________________________________________________________
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);
}

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