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