[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 "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)
#ifdef _ClusterTopology_
  ,fTreeTopology(0)
  ,fFileTopology(0)
  ,fTopology(0)
  ,fMinCol(0)
  ,fMinRow(0)
#endif //_ClusterTopology_
{
  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;
  //
#ifdef _ClusterTopology_
  SaveTopologyTree();
#endif //_ClusterTopology_
}

//______________________________________________________________________________
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);
#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(AliITSUClusterizerClusterCand *cand) 
{
  // fill special tree with cluster topology bit pattern
  //
  enum {kMaxColSpan=30,kMaxRowSpan=30};
  //
  Int_t mnCol=999999,mxCol=0,mnRow=999999,mxRow=0;
  for (AliITSUClusterizerClusterDigit *idigit=cand->fFirstDigit;idigit;idigit=idigit->fNext) {
    AliITSdigit* dig = idigit->fDigit;
    if (mnCol>dig->GetCoord1()) mnCol = dig->GetCoord1();
    if (mxCol<dig->GetCoord1()) mxCol = dig->GetCoord1();
    if (mnRow>dig->GetCoord2()) mnRow = dig->GetCoord2();
    if (mxRow<dig->GetCoord2()) mxRow = dig->GetCoord2();    
  }
  int nColSpan = mxCol-mnCol+1;
  int nRowSpan = mxRow-mnRow+1;
  if (nColSpan<kMaxColSpan && nRowSpan<kMaxRowSpan) {
    fTopology.Clear();
    fTopology.SetBitNumber(nColSpan*nRowSpan-1,kFALSE); // set the size
    fTopology.SetUniqueID( (nRowSpan<<8) + nColSpan );
    for (AliITSUClusterizerClusterDigit *idigit=cand->fFirstDigit;idigit;idigit=idigit->fNext) {
      AliITSdigit* dig = idigit->fDigit;
      fTopology.SetBitNumber( (dig->GetCoord2()-mnRow)*nColSpan + (dig->GetCoord1()-mnCol) );
    }
    if (!fTreeTopology) InitTopologyTree();
    AliInfo(Form("Fill Topology, lr %d",fLayerID));
    fMinCol = mnCol;
    fMinRow = mnRow;
    fTreeTopology->Fill();
  }
  //
}

//______________________________________________________________________
void AliITSUClusterizer::InitTopologyTree() 
{
  // create output tree for the topology stydy
  if (fFileTopology || fTreeTopology) AliFatal("Topology structures are already created");
  TDirectory* cdr = gDirectory;
  fFileTopology = TFile::Open("clusterTopology.root","recreate");
  fFileTopology->cd();
  fTreeTopology = new TTree("clTop","ClusterTopology");
  fTreeTopology->Branch("tp","TBits",&fTopology);
  fTreeTopology->Branch("volID",&fVolID);  
  fTreeTopology->Branch("mnCol",&fMinCol);  
  fTreeTopology->Branch("mnRow",&fMinRow);  
  cdr->cd();
}

//______________________________________________________________________
void AliITSUClusterizer::SaveTopologyTree() 
{
  // write output tree for the topology stydy
  AliInfo("Storing cluster topology tree");
  if (!fFileTopology || !fTreeTopology) {AliError("Topology structures are not created"); return;}
  fFileTopology->cd();
  fTreeTopology->Write();
  delete fTreeTopology;
  fFileTopology->Close();
  delete fFileTopology;
  fTreeTopology = 0;
  fFileTopology = 0;
}

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