[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;
  int charge=0;
  for (AliITSUClusterizerClusterDigit *idigit=cand->fFirstDigit;idigit;idigit=idigit->fNext) {
    AliITSdigit* dig = idigit->fDigit;
    fSegm->DetToLocal(dig->GetCoord2(),dig->GetCoord1(),cx,cz);
    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;
    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(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
    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;
29ad4146	160	int charge=0;
b69620f8	161	for (AliITSUClusterizerClusterDigit *idigit=cand->fFirstDigit;idigit;idigit=idigit->fNext) {
4bac12be	162	AliITSdigit* dig = idigit->fDigit;
0e84ce67	163	fSegm->DetToLocal(dig->GetCoord2(),dig->GetCoord1(),cx,cz);
29ad4146	164	charge += dig->GetSignal();
5e375bb4	165	x += cx;
	166	z += cz;
	167	if (cx<xmn) xmn=cx;
	168	if (cx>xmx) xmx=cx;
	169	if (cz<zmn) zmn=cz;
	170	if (cz>zmx) zmx=cz;
4bac12be	171	px += fSegm->Dpx(dig->GetCoord2());
	172	pz += fSegm->Dpz(dig->GetCoord1());
	173	//
	174	if (!fRawData) {
	175	for(Int_t dlab=0;dlab<maxLbinDigit;dlab++){
	176	Int_t digitlab = (dig->GetTracks())[dlab];
	177	if(digitlab<0) continue;
	178	AddLabel(digitlab);
	179	}
	180	}
	181	//
b69620f8	182	++n;
b69620f8	183	}
5e375bb4	184	UChar_t nx=1,nz=1;
	185	double dx = xmx-xmn, dz = zmx-zmn;
	186	if (n>1) {
	187	double fac=1./n;
	188	x *= fac; // mean coordinates
	189	z *= fac;
	190	px *= fac; // mean pitch
	191	pz *= fac;
	192	nx = 1+Nint(dx/px);
	193	nz = 1+Nint(dz/pz);
b69620f8	194	}
889b1493	195	x -= fLorAngCorrection; // LorentzAngle correction
b69620f8	196	cluster->SetX(x);
b69620f8	197	cluster->SetZ(z);
5e375bb4	198	cluster->SetY(0);
c8d1f258	199	cluster->SetSigmaZ2(nz>1 ? dzdzk1to12 : pzpzk1to12);
c8d1f258	200	cluster->SetSigmaY2(nx>1 ? dxdxk1to12 : pxpxk1to12);
5e375bb4	201	cluster->SetSigmaYZ(0);
5e375bb4	202	cluster->SetFrameLoc();
29ad4146	203	cluster->SetNxNzN(nx,nz,n);
29ad4146	204	cluster->SetQ(charge); // note: this is MC info
5e375bb4	205	//
4bac12be	206	if (!fRawData) {
	207	CheckLabels();
	208	int nl = Min(kMaxLabInCluster,fNLabels);
	209	for (int i=nl;i--;) cluster->SetLabel(fCurrLabels[i],i);
	210	}
	211	//
b69620f8	212	// Set Volume id
	213	cluster->SetVolumeId(fVolID);
	214	// printf("mod %d: (%.4lf,%.4lf)cm\n",fVolID,x,z);
	215	}
	216
	217	//______________________________________________________________________________
	218	void AliITSUClusterizer::CloseCand(AliITSUClusterizerClusterCand *cand)
	219	{
	220	// finish cluster
5e375bb4	221	AliITSUClusterPix cluster = (AliITSUClusterPix)NextCluster();
b69620f8	222	Transform(cluster,cand);
	223	DeallocDigits(cand->fFirstDigit,cand->fLastDigit);
	224	DeallocCand(cand);
	225	}
	226
	227	//______________________________________________________________________________
	228	void AliITSUClusterizer::ClosePart(AliITSUClusterizerClusterPart *part)
	229	{
	230	// finish cluster part
	231	AliITSUClusterizerClusterCand *cand=part->fParent;
	232	DetachPartFromCand(cand,part);
	233	DeallocPart(part);
	234	if (!cand->fFirstPart) CloseCand(cand);
	235	}
	236
	237	//______________________________________________________________________________
	238	void AliITSUClusterizer::CloseRemainingParts(AliITSUClusterizerClusterPart *part)
	239	{
	240	// finish what is left
	241	while (part) {
	242	AliITSUClusterizerClusterPart *next=part->fNextInRow;
	243	ClosePart(part);
	244	part=next;
	245	}
	246	}
	247
	248	//______________________________________________________________________________
	249	void AliITSUClusterizer::Clusterize()
	250	{
4bac12be	251	// main algo for MC clustererization
	252	SetRawData(kFALSE);
	253	//
b69620f8	254	AliITSUClusterizerClusterDigit *iDigit=NextDigit();
	255	AliITSUClusterizerClusterPart *iPrevRowBegin=0;
	256	AliITSUClusterizerClusterPart *iNextRowBegin=0;
	257	AliITSUClusterizerClusterPart *iPrevRow=0;
	258	AliITSUClusterizerClusterPart *iNextRow=0;
	259	Int_t lastV=0;
	260	while (iDigit) {
4bac12be	261	if (iDigit->fDigit->GetCoord2()!=lastV) {
b69620f8	262	// NEW ROW
	263	if (iNextRow) iNextRow->fNextInRow=0;
	264	if (iPrevRowBegin) CloseRemainingParts(iPrevRowBegin);
4bac12be	265	if (iDigit->fDigit->GetCoord2()==lastV+1) {
b69620f8	266	iPrevRowBegin=iNextRowBegin;
	267	iPrevRow =iNextRowBegin;
	268	}
	269	else {
	270	// there was an empty row
	271	CloseRemainingParts(iNextRowBegin);
	272	iPrevRowBegin=0;
	273	iPrevRow =0;
	274	}
	275	iNextRowBegin=0;
	276	iNextRow =0;
4bac12be	277	lastV=iDigit->fDigit->GetCoord2();
b69620f8	278	}
b69620f8	279	// skip cluster parts before this digit
4bac12be	280	while (iPrevRow && iPrevRow->fUEnd<iDigit->fDigit->GetCoord1()) {
b69620f8	281	iPrevRow=iPrevRow->fNextInRow;
	282	}
	283	// find the longest continous line of digits [iDigit,pDigit]=[iDigit,jDigit)
	284	AliITSUClusterizerClusterCand *cand=AllocCand();
	285	AliITSUClusterizerClusterDigit *pDigit=iDigit;
	286	AliITSUClusterizerClusterDigit *jDigit=NextDigit();
	287	cand->fFirstPart=0;
	288	cand->fFirstDigit=cand->fLastDigit=iDigit; // NB: first diggit is attached differently
	289	iDigit->fNext=0;
4bac12be	290	Int_t lastU =iDigit->fDigit->GetCoord1();
b69620f8	291	Int_t lastU1=lastU+1;
4bac12be	292	while (jDigit && jDigit->fDigit->GetCoord1()==lastU1 && jDigit->fDigit->GetCoord2()==lastV) {
b69620f8	293	pDigit=jDigit;
	294	jDigit=NextDigit();
	295	AttachDigitToCand(cand,pDigit);
	296	++lastU1;
	297	}
	298	--lastU1;
	299	AliITSUClusterizerClusterPart *part=AllocPart();
	300	part->fUBegin=lastU ;
	301	part->fUEnd =lastU1;
	302	AttachPartToCand(cand,part);
	303	// merge all cluster candidates of the previous line touching this one,
	304	// advance to the last one, but keep that one the next active one
	305	AliITSUClusterizerClusterPart *jPrevRow=iPrevRow;
	306	while (jPrevRow && jPrevRow->fUBegin<=lastU1) {
	307	if (jPrevRow->fParent!=cand) {
	308	MergeCands(jPrevRow->fParent,cand);
	309	DeallocCand(cand);
	310	cand=jPrevRow->fParent;
	311	}
	312	iPrevRow=jPrevRow;
	313	jPrevRow=jPrevRow->fNextInRow;
	314	}
	315	if (iNextRow)
	316	iNextRow->fNextInRow=part;
	317	else
	318	iNextRowBegin=part;
	319	iNextRow=part;
	320	iDigit=jDigit;
	321	}
	322	// remove remaining cluster parts
	323	CloseRemainingParts(iPrevRowBegin);
	324	if (iNextRow) iNextRow->fNextInRow=0;
	325	CloseRemainingParts(iNextRowBegin);
	326	return;
	327	}
889b1493	328
	329	//______________________________________________________________________________
	330	void AliITSUClusterizer::PrepareLorentzAngleCorrection(Double_t bz)
	331	{
	332	// calculate parameters for Lorentz Angle correction. Must be called
	333	// after setting segmentation and recoparams
	334	fLorAngCorrection = 0.5fRecoParam->GetTanLorentzAngle(fLayerID)bz/kNominalBz*fSegm->Dy();
	335	}
4bac12be	336
	337	//______________________________________________________________________
	338	void AliITSUClusterizer::CheckLabels()
	339	{
	340	// Tries to find mother's labels
	341	//
	342	if (fNLabels<1) return;
	343	AliRunLoader *rl = AliRunLoader::Instance();
	344	if(!rl) return;
	345	TTree trK=(TTree)rl->TreeK();
	346	if (!trK) return;
	347	//
	348	static int labS[kMaxLabels];
	349	static float kine[kMaxLabels];
	350	Int_t nlabels = fNLabels;
	351	Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
	352	for (Int_t i=fNLabels;i--;) labS[i] = fCurrLabels[i];
	353	//
	354	for (Int_t i=0;i<nlabels;i++) {
	355	Int_t label = labS[i];
	356	if (label>=ntracks) continue;
	357	TParticle part=(TParticle)gAlice->GetMCApp()->Particle(label);
	358	kine[i] = part->Energy() - part->GetCalcMass(); // kinetic energy
	359	if (kine[i] < 0.02) { // reduce soft particles from the same cluster
	360	Int_t m=part->GetFirstMother();
	361	if (m<0) continue; // primary
	362	//
	363	if (part->GetStatusCode()>0) continue;
	364	//
	365	// if the parent is within the same cluster, assign parent's label
	366	for (int j=0;j<nlabels;j++) if (labS[j]==m) { labS[i] = m; break;}
	367	}
	368	}
	369	//
	370	if (nlabels>kMaxLabInCluster) { // only 3 labels are stored in cluster, sort in decreasing momentum
	371	static int ind[kMaxLabels],labSS[kMaxLabels];
	372	TMath::Sort(nlabels,kine,ind);
	373	for (int i=nlabels;i--;) labSS[i] = labS[i];
	374	for (int i=nlabels;i--;) labS[i] = labSS[ind[i]];
	375	}
	376	//
	377	//compress labels -- if multi-times the same
	378	for (Int_t i=0;i<nlabels;i++) fCurrLabels[i]=-2;
	379	fNLabels = 0;
	380	int j=0;
	381	for (int i=0;i<nlabels;i++) {
	382	for (j=fNLabels;j--;) if (labS[i]==fCurrLabels[j]) break; // the label already there
	383	if (j<0) fCurrLabels[fNLabels++] = labS[i];
	384	}
	385	//
	386	}