[u/mrichter/AliRoot.git] / TRD / AliTRDstackLayer.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//  Organization of clusters at the level of 1 TRD chamber.                  //
//  The data structure is used for tracking at the stack level.              //
//                                                                           //
//  Functionalities:                                                         //
//   1. cluster organization and sorting                                     //
//   2. fast data navigation                                                 //
//                                                                           //
//  Authors:                                                                 //
//    Alex Bercuci <A.Bercuci@gsi.de>                                        //
//    Markus Fasel <M.Fasel@gsi.de>                                          //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////

#include <TObject.h>
#include <TROOT.h>
#include <TMath.h>
#include <TStopwatch.h>
#include <TTreeStream.h>

#include "AliLog.h"
#include "AliTRDcluster.h"

#include "AliTRDstackLayer.h"
#include "AliTRDrecoParam.h"
#include "AliTRDReconstructor.h"

#define DEBUG

ClassImp(AliTRDstackLayer)

//_____________________________________________________________________________
AliTRDstackLayer::AliTRDstackLayer(Double_t z0, Double_t zLength
                                 , UChar_t stackNr, AliTRDrecoParam *p)
  :AliTRDpropagationLayer()
  ,fOwner(kFALSE)
  ,fStackNr(stackNr)
  ,fNRows(kMaxRows)
  ,fZ0(z0)
  ,fZLength(zLength)
  ,fRecoParam(p)
  ,fDebugStream(0x0)
{
  //
  // Default constructor (Only provided to use AliTRDstackLayer with arrays)
  //

	for(int i=0; i<kMaxRows; i++) fPositions[i] = 0;
}

//_____________________________________________________________________________
AliTRDstackLayer::AliTRDstackLayer(const AliTRDpropagationLayer &layer, Double_t
z0, Double_t zLength, UChar_t stackNr, AliTRDrecoParam *p):
	AliTRDpropagationLayer(layer)
	,fOwner(kFALSE)
	,fStackNr(stackNr)
	,fNRows(kMaxRows)
	,fZ0(z0)
	,fZLength(zLength)
	,fRecoParam(p)
	,fDebugStream(0x0)
{
// Standard constructor.
// Initialize also the underlying AliTRDpropagationLayer using the copy constructor.

	SetT0(layer.IsT0());
	for(int i=0; i<kMaxRows; i++) fPositions[i] = 0;
}

//_____________________________________________________________________________
AliTRDstackLayer::AliTRDstackLayer(const AliTRDpropagationLayer &layer):
	AliTRDpropagationLayer(layer)
	,fOwner(kFALSE)
	,fStackNr(0)
	,fNRows(kMaxRows)
	,fZ0(0)
	,fZLength(0)
	,fRecoParam(0x0)
	,fDebugStream(0x0)
{
// Standard constructor using only AliTRDpropagationLayer.
	
	SetT0(layer.IsT0());
	for(int i=0; i<kMaxRows; i++) fPositions[i] = 0;
}

//_____________________________________________________________________________
AliTRDstackLayer::AliTRDstackLayer(const AliTRDstackLayer &layer):
	AliTRDpropagationLayer(layer)
	,fOwner(layer.fOwner)
	,fStackNr(layer.fStackNr)
	,fNRows(layer.fNRows)
	,fZ0(layer.fZ0)
	,fZLength(layer.fZLength)
	,fRecoParam(layer.fRecoParam)
	,fDebugStream(layer.fDebugStream)
{
// Copy Constructor (performs a deep copy)
	
	SetT0(layer.IsT0());
	for(Int_t i = 0; i < kMaxRows; i++) fPositions[i] = layer.fPositions[i];
// 	BuildIndices();
}

//_____________________________________________________________________________
AliTRDstackLayer &AliTRDstackLayer::operator=(const AliTRDpropagationLayer &layer)
{
// Assignment operator from an AliTRDpropagationLayer

	if (this != &layer) layer.Copy(*this);
	return *this;
}

//_____________________________________________________________________________
AliTRDstackLayer &AliTRDstackLayer::operator=(const AliTRDstackLayer &layer)
{
// Assignment operator

	if (this != &layer) layer.Copy(*this);
  return *this;
}

//_____________________________________________________________________________
void AliTRDstackLayer::Copy(TObject &o) const
{
// Copy method. Performs a deep copy of all data from this object to object o.
	
	AliTRDstackLayer &layer = (AliTRDstackLayer &)o;
	layer.fZ0          = fZ0;
	layer.fOwner       = kFALSE;
	layer.fNRows       = fNRows;
	layer.fZLength     = fZLength;
	layer.fStackNr     = fStackNr;
	layer.fDebugStream = fDebugStream;
	layer.fRecoParam   = fRecoParam;
	layer.SetT0(IsT0());
	
	AliTRDpropagationLayer::Copy(layer); // copies everything into layer
	for(UChar_t i = 0; i < kMaxRows; i++) layer.fPositions[i] = 0;
// 	layer.BuildIndices();
}

//_____________________________________________________________________________
AliTRDstackLayer::~AliTRDstackLayer()
{
// Destructor
	if(fOwner) for(int ic=0; ic<fN; ic++) delete fClusters[ic];
}

//_____________________________________________________________________________
void AliTRDstackLayer::SetRange(const Float_t z0, const Float_t zLength)
{
// Sets the range in z-direction
//
// Parameters:
//   z0      : starting position of layer in the z direction
//   zLength : length of layer in the z direction 

	fZ0 = (z0 <= z0 + zLength) ? z0 : z0 + zLength;
	fZLength = TMath::Abs(zLength);
}

//_____________________________________________________________________________
void AliTRDstackLayer::BuildIndices(Int_t iter)
{
// Rearrangement of the clusters belonging to the propagation layer for the stack.
//
// Detailed description
//
// The array indices of all clusters in one PropagationLayer are stored in
// array. The array is divided into several bins.
// The clusters are sorted in increasing order of their y coordinate.
//
// Sorting algorithm: TreeSearch
//

	if(!fN) return;

	// Select clusters that belong to the Stack
	Int_t nClStack = 0;					// Internal counter
	for(Int_t i = 0; i < fN; i++){
		Double_t zval = fClusters[i]->GetZ();
		if(zval < fZ0 || zval > fZ0 + fZLength || fClusters[i]->IsUsed()){
			fClusters[i] = 0x0;
			fIndex[i] = 9999;
		} else nClStack++;
	}
	if(nClStack > kMaxClustersLayer) AliWarning(Form("Number of clusters in stack %d exceed buffer size %d", nClStack, kMaxClustersLayer));
	
	// Nothing in this Stack 
	if(!nClStack){
		delete fClusters;
		delete fIndex;
		fClusters = 0x0;
		fIndex = 0x0;
		fN = 0;
		memset(fPositions, 0, sizeof(UChar_t) * 16);
		return;
	}
	
	// Make a copy
	AliTRDcluster *helpCL[kMaxClustersLayer];
	Int_t helpInd[kMaxClustersLayer];
	nClStack = 0;
	for(Int_t i = 0; i < TMath::Min(fN, kMaxClustersLayer); i++){
		if(!fClusters[i]) continue;
		helpCL[nClStack]  = fClusters[i];
		helpInd[nClStack] = fIndex[i];
		fClusters[i]      = 0x0;
		fIndex[i]         = 9999;
		nClStack++;
	}
	
	// do clusters arrangement
	fN =  nClStack;
	nClStack = 0;
	memset(fPositions,0,sizeof(UChar_t)*16);				// Reset Positions array
	for(Int_t i = 0; i < fN; i++){
		// boundarie check
		AliTRDcluster *cl = helpCL[i];
		Double_t zval = cl->GetZ();
		UChar_t treeIndex = (UChar_t)(TMath::Abs(fZ0 - zval)/fZLength * fNRows);
		if(treeIndex > fNRows - 1) treeIndex = fNRows - 1;
		// Insert Leaf
		Int_t pos = FindYPosition(cl->GetY(), treeIndex, i);
		if(pos == -1){		// zbin is empty;
			Int_t upper = (treeIndex == fNRows - 1) ? nClStack : fPositions[treeIndex + 1];
			memmove(fClusters + upper + 1, fClusters + upper, (sizeof(AliTRDcluster *))*(nClStack-upper));
			memmove(fIndex + upper + 1, fIndex + upper, (sizeof(UInt_t))*(nClStack-upper));
			fClusters[upper] = cl;
			fIndex[upper] = helpInd[i]; 
			// Move All pointer one position back
			for(UChar_t j = treeIndex + 1; j < fNRows; j++) fPositions[j]++;
			nClStack++;
		} else {		// zbin not empty
			memmove(fClusters + pos + 2, fClusters + pos+1, (sizeof(AliTRDcluster *))*(nClStack-(pos+1)));
			memmove(fIndex + pos + 2, fIndex + pos+1, (sizeof(UInt_t))*(nClStack-(pos+1)));
			fClusters[pos + 1] = cl;	//fIndex[i];
			fIndex[pos + 1] = helpInd[i];
			// Move All pointer one position back
			for(UChar_t j = treeIndex + 1; j < fNRows; j++) fPositions[j]++;	
			nClStack++;
		}

		
		// Debug Streaming
#ifdef DEBUG
		if(fDebugStream && AliTRDReconstructor::StreamLevel() >= 3){
			TTreeSRedirector &cstream = *fDebugStream;
			cstream << "BuildIndices"
			<< "StackNr="  << fStackNr
			<< "SectorNr=" << fSec
			<< "Iter="     << iter
			<< "C.="       << cl
			<< "TreeIdx="  << treeIndex
			<< "\n";
		}
#endif	
	}
}

//_____________________________________________________________________________
Int_t AliTRDstackLayer::FindYPosition(Double_t y, UChar_t z, Int_t nClusters) const
{
//
// Tree search Algorithm to find the nearest left cluster for a given
// y-position in a certain z-bin (in fact AVL-tree). 
// Making use of the fact that clusters are sorted in y-direction.
//
// Parameters:
//   y : y position of the reference point in tracking coordinates
//   z : z reference bin.
//   nClusters : 
//
// Output :
// Index of the nearest left cluster in the StackLayer indexing (-1 if no clusters are found)
//

	Int_t start = fPositions[z]; 		// starting Position of the bin
	Int_t upper = (Int_t)((z != fNRows - 1) ? fPositions[z+1] : nClusters); // ending Position of the bin 
	Int_t end = upper - 1; // ending Position of the bin 
	if(end < start) return -1; // Bin is empty
	Int_t middle = static_cast<Int_t>((start + end)/2);
	// 1st Part: climb down the tree: get the next cluster BEFORE ypos
	while(start + 1 < end){
		if(y >= fClusters[middle]->GetY()) start = middle;
		else end = middle;
		middle = static_cast<Int_t>((start + end)/2);
	}
	if(y > fClusters[end]->GetY()) return end;
	return start;
}

//_____________________________________________________________________________
Int_t AliTRDstackLayer::FindNearestYCluster(Double_t y, UChar_t z) const
{
//
// Tree search Algorithm to find the nearest cluster for a given
// y-position in a certain z-bin (in fact AVL-tree). 
// Making use of the fact that clusters are sorted in y-direction.
//
// Parameters:
//   y : y position of the reference point in tracking coordinates
//   z : z reference bin.
//
// Output 
// Index of the nearest cluster in the StackLayer indexing (-1 if no clusters are found)
//

	Int_t position = FindYPosition(y, z, fN);
	if(position == -1) return position; // bin empty
	// FindYPosition always returns the left Neighbor. We don't know if the left or the right Neighbor is nearest
	// to the Reference y-position, so test both
	Int_t upper = (Int_t)((z < fNRows-1) ? fPositions[z+1] : fN); // ending Position of the bin
	if((position + 1) < (upper)){
		if(TMath::Abs(y - fClusters[position + 1]->GetY()) < TMath::Abs(y - fClusters[position]->GetY())) return position + 1;
		else return position;
	}
	return position;
}

//_____________________________________________________________________________
Int_t AliTRDstackLayer::SearchNearestCluster(Double_t y, Double_t z, Double_t maxroady, Double_t maxroadz) const
{
//
// Finds the nearest cluster from a given point in a defined range.
// Distance is determined in a 2D space by the 2-Norm.
//
// Parameters :
//   y : y position of the reference point in tracking coordinates
//   z : z reference bin.
//   maxroady : maximum searching distance in y direction
//   maxroadz : maximum searching distance in z direction
//
// Output 
// Index of the nearest cluster in the StackLayer indexing (-1 if no cluster is found).
// Cluster can be accessed with the operator[] or GetCluster(Int_t index)
//
// Detail description
//
// The following steps are perfomed:
// 1. Get the expected z bins inside maxroadz.
// 2. For each z bin find nearest y cluster.
// 3. Select best candidate
//
	Int_t   index   = -1;
	// initial minimal distance will be represented as ellipse: semi-major = z-direction
	// later 2-Norm will be used  
// 	Float_t nExcentricity = TMath::Sqrt(maxroadz*maxroadz - maxroad*maxroad)/maxroadz;
	Float_t mindist = maxroadz;
	
	// not very nice but unfortunately neccessarry: we have ho check the neighbors in both directions (+ and -) too. How 
	// much neighbors depends on the Quotient maxroadz/fZLength   
	UChar_t maxRows = 3;
	UChar_t zpos[kMaxRows];
  // Float_t mindist = TMath::Sqrt(maxroad*maxroad + maxroadz*maxroadz);
// 	UChar_t myZbin = FindTreePosition(z, fZ0 + fZLength/2, fZLength/4, 8, 8, kFALSE);
	UChar_t myZbin = (UChar_t)(TMath::Abs(fZ0 - z)/fZLength * fNRows);
	if(z < fZ0) myZbin = 0;
	if(z > fZ0 + fZLength) myZbin = fNRows - 1;

	UChar_t nNeighbors = 0;
	for(UChar_t i = 0; i < maxRows; i++){
		if((myZbin - 1 + i) < 0) continue;
		if((myZbin - 1 + i) > fNRows - 1) break;
		zpos[nNeighbors] = myZbin - 1 + i;
		nNeighbors++;
	}
	Float_t ycl = 0, zcl = 0;
	for(UChar_t neighbor = 0; neighbor < nNeighbors; neighbor++){	// Always test the neighbors too
		Int_t pos = FindNearestYCluster(y,zpos[neighbor]);
		if(pos == -1) continue;	// No cluster in bin
		AliTRDcluster *c = (AliTRDcluster *) (fClusters[pos]);
		if(c->IsUsed()) continue;		// we are only interested in unused clusters
		ycl = c->GetY();
		// Too far away in y-direction (Prearrangement)
		if (TMath::Abs(ycl - y) > maxroady) continue;

		zcl = c->GetZ();
		// Too far away in z-Direction
		// (Prearrangement since we have not so many bins to test)
		if (TMath::Abs(zcl - z) > maxroadz) continue;
		
		Float_t dist;		// distance defined as 2-Norm	
		// if we havent found a Particle that is in the ellipse around (y,z) with maxroad as semi-minor and
		// maxroadz as semi-major, we take the radius of the ellipse concerning the cluster as mindist, later we 
		// take the 2-Norm when we found a cluster inside the ellipse (The value 10000 is taken because it is surely
		// large enough to be usable as an indicator whether we have found a nearer cluster or not)
// 		if(mindist > 10000.){
// 			Float_t phi = ((zcl - z) == 0) ? TMath::Pi()/2 : TMath::ATan((ycl - y)/(zcl - z));
// 			mindist = maxroad/TMath::Sqrt(1 - nExcentricity*nExcentricity * (TMath::Cos(phi))*(TMath::Cos(phi)));
// 		}
		dist = TMath::Max(TMath::Abs(y-ycl),TMath::Abs(z-zcl));	// infinity Norm
// 		dist = TMath::Sqrt((ycl - y)*(ycl - y) + (zcl - z)*(zcl - z));
		if((Int_t)(dist * 100000) < (Int_t)(mindist * 100000)){
		//if((dist = TMath::Sqrt((ycl - y)*(ycl - y) + (zcl - z)*(zcl - z))) < mindist){
			mindist = dist;
			index   = pos;
		}	
	}
	// This is the Array Position in fIndex2D of the Nearest cluster: if a
	// cluster is called, then the function has to retrieve the Information
	// which is Stored in the Array called, the function
	return index;
}

//_____________________________________________________________________________
void AliTRDstackLayer::BuildCond(AliTRDcluster *cl, Double_t *cond, UChar_t Layer, Double_t theta, Double_t phi)
{
// Helper function to calculate the area where to expect a cluster in THIS
// layer. 
//
// Parameters :
//   cl    : 
//   cond  :
//   Layer : 
//   theta : 
//   phi   : 
//
// Detail description
//
// Helper function to calculate the area where to expect a cluster in THIS
// layer. by using the information of a former cluster in another layer
// and the angle in theta- and phi-direction between layer 0 and layer 3.
// If the layer is zero, initial conditions are calculated. Otherwise a
// linear interpolation is performed. 
//Begin_Html
//<img src="gif/build_cond.gif">
//End_Html
//

	if(!fRecoParam){
		AliError("Reconstruction parameters not initialized.");
		return;
	}
	
	if(Layer == 0){
		cond[0] = cl->GetY();			// center: y-Direction
		cond[1] = cl->GetZ();			// center: z-Direction
		cond[2] = fRecoParam->GetMaxPhi()   * (cl->GetX() - GetX()) + 1.0;			// deviation: y-Direction
		cond[3] = fRecoParam->GetMaxTheta() * (cl->GetX() - GetX()) + 1.0;			// deviation: z-Direction
	} else {
		cond[0] = cl->GetY() + phi   * (GetX() - cl->GetX()); 
		cond[1] = cl->GetZ() + theta * (GetX() - cl->GetX());
		cond[2] = fRecoParam->GetRoad0y() + phi;
		cond[3] = fRecoParam->GetRoad0z();
	}
}

//_____________________________________________________________________________
void AliTRDstackLayer::GetClusters(Double_t *cond, Int_t *index, Int_t& ncl, Int_t BufferSize)
{
// Finds all clusters situated in this layer inside a rectangle  given by the center an ranges.
//
// Parameters :
//   cond  :
//   index : 
//   ncl :
//   BufferSize   :
//
// Output :
//
// Detail description
//
// Function returs an array containing the indices in the stacklayer of
// the clusters found an  the number of found clusters in the stacklayer

	ncl = 0;
	memset(index, 0, BufferSize*sizeof(Int_t));
	if(fN == 0) return;
	
	//Boundary checks
	Double_t zvals[2];
	zvals[0] = ((cond[1] - cond[3]) < fZ0) ? fZ0 : (cond[1] - cond[3]);
	zvals[1] = ((cond[1] + cond[3]) < fZ0 + fZLength) ? (cond[1] + cond[3]) : fZ0 + fZLength;

	UChar_t zlo = (fZ0>zvals[0]) ? 0 : (UChar_t)(TMath::Abs(fZ0 - zvals[0])/fZLength * fNRows);
	UChar_t zhi = (fZ0+fZLength<zvals[1]) ? fNRows - 1 : (UChar_t)(TMath::Abs(fZ0 - zvals[1])/fZLength * fNRows);

	//Preordering in Direction z saves a lot of loops (boundary checked)
	for(UChar_t z = zlo; z <= zhi; z++){
		UInt_t upper = (z < fNRows-1) ? fPositions[z+1] : fN;
		for(Int_t y = fPositions[z]; y < (Int_t)upper; y++){
			if(ncl == BufferSize){
				AliWarning("Buffer size riched. Some clusters may be lost.");
				return;	//Buffer filled
			}
			if(fClusters[y]->GetY() > (cond[0] + cond[2])) break;			// Abbortion conditions!!!
			if(fClusters[y]->GetY() < (cond[0] - cond[2])) continue;	// Too small
			if(((Int_t)((fClusters[y]->GetZ())*1000) < (Int_t)(zvals[0]*1000)) || ((Int_t)((fClusters[y]->GetZ())*1000) > (Int_t)(zvals[1]*1000))) continue;
			index[ncl] = y;
			ncl++;
		}
	}
	if(ncl>fN) AliError(Form("Clusters found %d > %d (clusters in layer)", ncl, fN));
}

//_____________________________________________________________________________
AliTRDcluster *AliTRDstackLayer::GetNearestCluster(Double_t *cond)
{
// Function returning a pointer to the nearest cluster (nullpointer if not successfull).
//
// Parameters :
//   cond  :
//
// Output :
//   pointer to the nearest cluster (nullpointer if not successfull).
// 
// Detail description
//
// returns a pointer to the nearest cluster (nullpointer if not
// successfull) by the help of the method FindNearestCluster
	
	
	Double_t maxroad  = fRecoParam->GetRoad2y();
	Double_t maxroadz = fRecoParam->GetRoad2z();
	
	Int_t index = SearchNearestCluster(cond[0],cond[1],maxroad,maxroadz);
	AliTRDcluster *returnCluster = 0x0;
	if(index != -1) returnCluster = (AliTRDcluster *) fClusters[index];
	return returnCluster;
}

//_____________________________________________________________________________
void AliTRDstackLayer::PrintClusters() const
{
// Prints the position of each cluster in the stacklayer on the stdout
//
  //printf("fDebugStream = %#o\n", ((Int_t) fDebugStream));
  //printf("fRecoParam   = %#o\n", ((Int_t) fRecoParam));
	
	for(Int_t i = 0; i < fN; i++){
		printf("AliTRDstackLayer: index=%i, Cluster: X = %3.3f, Y = %3.3f, Z = %3.3f\n", i,  fClusters[i]->GetX(),fClusters[i]->GetY(),fClusters[i]->GetZ());
		if(fClusters[i]->IsUsed()) printf("cluster allready used. rejected in search algorithm\n");
	}
}
Commit	Line	Data
e4f2f73d	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id$ */
	17
	18	///////////////////////////////////////////////////////////////////////////////
	19	// //
	20	// Organization of clusters at the level of 1 TRD chamber. //
	21	// The data structure is used for tracking at the stack level. //
	22	// //
	23	// Functionalities: //
	24	// 1. cluster organization and sorting //
	25	// 2. fast data navigation //
	26	// //
	27	// Authors: //
	28	// Alex Bercuci <A.Bercuci@gsi.de> //
	29	// Markus Fasel <M.Fasel@gsi.de> //
	30	// //
	31	///////////////////////////////////////////////////////////////////////////////
	32
	33	#include <TObject.h>
	34	#include <TROOT.h>
	35	#include <TMath.h>
	36	#include <TStopwatch.h>
	37	#include <TTreeStream.h>
	38
	39	#include "AliLog.h"
	40	#include "AliTRDcluster.h"
	41
	42	#include "AliTRDstackLayer.h"
	43	#include "AliTRDrecoParam.h"
	44	#include "AliTRDReconstructor.h"
	45
	46	#define DEBUG
	47
	48	ClassImp(AliTRDstackLayer)
	49
	50	//_____________________________________________________________________________
	51	AliTRDstackLayer::AliTRDstackLayer(Double_t z0, Double_t zLength
	52	, UChar_t stackNr, AliTRDrecoParam *p)
	53	:AliTRDpropagationLayer()
	54	,fOwner(kFALSE)
	55	,fStackNr(stackNr)
	56	,fNRows(kMaxRows)
	57	,fZ0(z0)
	58	,fZLength(zLength)
	59	,fRecoParam(p)
	60	,fDebugStream(0x0)
	61	{
	62	//
	63	// Default constructor (Only provided to use AliTRDstackLayer with arrays)
	64	//
65
66	for(int i=0; i<kMaxRows; i++) fPositions[i] = 0;
67	}
68
69	//_____________________________________________________________________________
70	AliTRDstackLayer::AliTRDstackLayer(const AliTRDpropagationLayer &layer, Double_t
71	z0, Double_t zLength, UChar_t stackNr, AliTRDrecoParam *p):
72	AliTRDpropagationLayer(layer)
73	,fOwner(kFALSE)
74	,fStackNr(stackNr)
75	,fNRows(kMaxRows)
76	,fZ0(z0)
77	,fZLength(zLength)
78	,fRecoParam(p)
79	,fDebugStream(0x0)
80	{
81	// Standard constructor.
82	// Initialize also the underlying AliTRDpropagationLayer using the copy constructor.
83
bcb6fb78	84	SetT0(layer.IsT0());
e4f2f73d	85	for(int i=0; i<kMaxRows; i++) fPositions[i] = 0;
	86	}
	87
	88	//_____________________________________________________________________________
	89	AliTRDstackLayer::AliTRDstackLayer(const AliTRDpropagationLayer &layer):
	90	AliTRDpropagationLayer(layer)
	91	,fOwner(kFALSE)
	92	,fStackNr(0)
	93	,fNRows(kMaxRows)
	94	,fZ0(0)
	95	,fZLength(0)
	96	,fRecoParam(0x0)
	97	,fDebugStream(0x0)
	98	{
	99	// Standard constructor using only AliTRDpropagationLayer.
	100
bcb6fb78	101	SetT0(layer.IsT0());
e4f2f73d	102	for(int i=0; i<kMaxRows; i++) fPositions[i] = 0;
	103	}
	104
	105	//_____________________________________________________________________________
	106	AliTRDstackLayer::AliTRDstackLayer(const AliTRDstackLayer &layer):
	107	AliTRDpropagationLayer(layer)
	108	,fOwner(layer.fOwner)
	109	,fStackNr(layer.fStackNr)
	110	,fNRows(layer.fNRows)
	111	,fZ0(layer.fZ0)
	112	,fZLength(layer.fZLength)
	113	,fRecoParam(layer.fRecoParam)
	114	,fDebugStream(layer.fDebugStream)
	115	{
	116	// Copy Constructor (performs a deep copy)
	117
bcb6fb78	118	SetT0(layer.IsT0());
e4f2f73d	119	for(Int_t i = 0; i < kMaxRows; i++) fPositions[i] = layer.fPositions[i];
	120	// BuildIndices();
	121	}
	122
	123	//_____________________________________________________________________________
	124	AliTRDstackLayer &AliTRDstackLayer::operator=(const AliTRDpropagationLayer &layer)
	125	{
	126	// Assignment operator from an AliTRDpropagationLayer
	127
	128	if (this != &layer) layer.Copy(*this);
	129	return *this;
	130	}
	131
	132	//_____________________________________________________________________________
	133	AliTRDstackLayer &AliTRDstackLayer::operator=(const AliTRDstackLayer &layer)
	134	{
	135	// Assignment operator
	136
	137	if (this != &layer) layer.Copy(*this);
	138	return *this;
	139	}
	140
	141	//_____________________________________________________________________________
	142	void AliTRDstackLayer::Copy(TObject &o) const
	143	{
	144	// Copy method. Performs a deep copy of all data from this object to object o.
	145
	146	AliTRDstackLayer &layer = (AliTRDstackLayer &)o;
	147	layer.fZ0 = fZ0;
	148	layer.fOwner = kFALSE;
	149	layer.fNRows = fNRows;
	150	layer.fZLength = fZLength;
	151	layer.fStackNr = fStackNr;
	152	layer.fDebugStream = fDebugStream;
	153	layer.fRecoParam = fRecoParam;
bcb6fb78	154	layer.SetT0(IsT0());
e4f2f73d	155
	156	AliTRDpropagationLayer::Copy(layer); // copies everything into layer
	157	for(UChar_t i = 0; i < kMaxRows; i++) layer.fPositions[i] = 0;
	158	// layer.BuildIndices();
	159	}
	160
	161	//_____________________________________________________________________________
	162	AliTRDstackLayer::~AliTRDstackLayer()
	163	{
	164	// Destructor
	165	if(fOwner) for(int ic=0; ic<fN; ic++) delete fClusters[ic];
	166	}
	167
	168	//_____________________________________________________________________________
	169	void AliTRDstackLayer::SetRange(const Float_t z0, const Float_t zLength)
	170	{
	171	// Sets the range in z-direction
	172	//
	173	// Parameters:
	174	// z0 : starting position of layer in the z direction
	175	// zLength : length of layer in the z direction
	176
	177	fZ0 = (z0 <= z0 + zLength) ? z0 : z0 + zLength;
	178	fZLength = TMath::Abs(zLength);
	179	}
	180
	181	//_____________________________________________________________________________
	182	void AliTRDstackLayer::BuildIndices(Int_t iter)
	183	{
	184	// Rearrangement of the clusters belonging to the propagation layer for the stack.
	185	//
	186	// Detailed description
	187	//
	188	// The array indices of all clusters in one PropagationLayer are stored in
	189	// array. The array is divided into several bins.
	190	// The clusters are sorted in increasing order of their y coordinate.
	191	//
	192	// Sorting algorithm: TreeSearch
	193	//
	194
	195	if(!fN) return;
	196
	197	// Select clusters that belong to the Stack
	198	Int_t nClStack = 0; // Internal counter
	199	for(Int_t i = 0; i < fN; i++){
	200	Double_t zval = fClusters[i]->GetZ();
	201	if(zval < fZ0 \|\| zval > fZ0 + fZLength \|\| fClusters[i]->IsUsed()){
	202	fClusters[i] = 0x0;
	203	fIndex[i] = 9999;
	204	} else nClStack++;
	205	}
	206	if(nClStack > kMaxClustersLayer) AliWarning(Form("Number of clusters in stack %d exceed buffer size %d", nClStack, kMaxClustersLayer));
	207
	208	// Nothing in this Stack
	209	if(!nClStack){
	210	delete fClusters;
	211	delete fIndex;
	212	fClusters = 0x0;
	213	fIndex = 0x0;
	214	fN = 0;
	215	memset(fPositions, 0, sizeof(UChar_t) * 16);
	216	return;
	217	}
	218
219	// Make a copy
220	AliTRDcluster *helpCL[kMaxClustersLayer];
221	Int_t helpInd[kMaxClustersLayer];
222	nClStack = 0;
223	for(Int_t i = 0; i < TMath::Min(fN, kMaxClustersLayer); i++){
224	if(!fClusters[i]) continue;
225	helpCL[nClStack] = fClusters[i];
226	helpInd[nClStack] = fIndex[i];
227	fClusters[i] = 0x0;
228	fIndex[i] = 9999;
229	nClStack++;
230	}
231
232	// do clusters arrangement
233	fN = nClStack;
234	nClStack = 0;
235	memset(fPositions,0,sizeof(UChar_t)*16); // Reset Positions array
236	for(Int_t i = 0; i < fN; i++){
237	// boundarie check
238	AliTRDcluster *cl = helpCL[i];
239	Double_t zval = cl->GetZ();
d76231c8	240	UChar_t treeIndex = (UChar_t)(TMath::Abs(fZ0 - zval)/fZLength * fNRows);
d76231c8	241	if(treeIndex > fNRows - 1) treeIndex = fNRows - 1;
e4f2f73d	242	// Insert Leaf
d76231c8	243	Int_t pos = FindYPosition(cl->GetY(), treeIndex, i);
e4f2f73d	244	if(pos == -1){ // zbin is empty;
d76231c8	245	Int_t upper = (treeIndex == fNRows - 1) ? nClStack : fPositions[treeIndex + 1];
e4f2f73d	246	memmove(fClusters + upper + 1, fClusters + upper, (sizeof(AliTRDcluster ))(nClStack-upper));
	247	memmove(fIndex + upper + 1, fIndex + upper, (sizeof(UInt_t))*(nClStack-upper));
	248	fClusters[upper] = cl;
	249	fIndex[upper] = helpInd[i];
	250	// Move All pointer one position back
d76231c8	251	for(UChar_t j = treeIndex + 1; j < fNRows; j++) fPositions[j]++;
e4f2f73d	252	nClStack++;
	253	} else { // zbin not empty
	254	memmove(fClusters + pos + 2, fClusters + pos+1, (sizeof(AliTRDcluster ))(nClStack-(pos+1)));
	255	memmove(fIndex + pos + 2, fIndex + pos+1, (sizeof(UInt_t))*(nClStack-(pos+1)));
	256	fClusters[pos + 1] = cl; //fIndex[i];
	257	fIndex[pos + 1] = helpInd[i];
	258	// Move All pointer one position back
d76231c8	259	for(UChar_t j = treeIndex + 1; j < fNRows; j++) fPositions[j]++;
e4f2f73d	260	nClStack++;
	261	}
	262
	263
	264	// Debug Streaming
	265	#ifdef DEBUG
	266	if(fDebugStream && AliTRDReconstructor::StreamLevel() >= 3){
	267	TTreeSRedirector &cstream = *fDebugStream;
	268	cstream << "BuildIndices"
	269	<< "StackNr=" << fStackNr
	270	<< "SectorNr=" << fSec
	271	<< "Iter=" << iter
	272	<< "C.=" << cl
d76231c8	273	<< "TreeIdx=" << treeIndex
e4f2f73d	274	<< "\n";
	275	}
	276	#endif
	277	}
	278	}
	279
	280	//_____________________________________________________________________________
	281	Int_t AliTRDstackLayer::FindYPosition(Double_t y, UChar_t z, Int_t nClusters) const
	282	{
	283	//
	284	// Tree search Algorithm to find the nearest left cluster for a given
	285	// y-position in a certain z-bin (in fact AVL-tree).
	286	// Making use of the fact that clusters are sorted in y-direction.
	287	//
	288	// Parameters:
	289	// y : y position of the reference point in tracking coordinates
	290	// z : z reference bin.
	291	// nClusters :
	292	//
	293	// Output :
	294	// Index of the nearest left cluster in the StackLayer indexing (-1 if no clusters are found)
	295	//
	296
	297	Int_t start = fPositions[z]; // starting Position of the bin
	298	Int_t upper = (Int_t)((z != fNRows - 1) ? fPositions[z+1] : nClusters); // ending Position of the bin
	299	Int_t end = upper - 1; // ending Position of the bin
	300	if(end < start) return -1; // Bin is empty
	301	Int_t middle = static_cast<Int_t>((start + end)/2);
	302	// 1st Part: climb down the tree: get the next cluster BEFORE ypos
	303	while(start + 1 < end){
	304	if(y >= fClusters[middle]->GetY()) start = middle;
	305	else end = middle;
	306	middle = static_cast<Int_t>((start + end)/2);
	307	}
	308	if(y > fClusters[end]->GetY()) return end;
	309	return start;
	310	}
	311
	312	//_____________________________________________________________________________
	313	Int_t AliTRDstackLayer::FindNearestYCluster(Double_t y, UChar_t z) const
	314	{
	315	//
	316	// Tree search Algorithm to find the nearest cluster for a given
	317	// y-position in a certain z-bin (in fact AVL-tree).
	318	// Making use of the fact that clusters are sorted in y-direction.
	319	//
	320	// Parameters:
	321	// y : y position of the reference point in tracking coordinates
	322	// z : z reference bin.
	323	//
	324	// Output
	325	// Index of the nearest cluster in the StackLayer indexing (-1 if no clusters are found)
	326	//
	327
	328	Int_t position = FindYPosition(y, z, fN);
	329	if(position == -1) return position; // bin empty
	330	// FindYPosition always returns the left Neighbor. We don't know if the left or the right Neighbor is nearest
	331	// to the Reference y-position, so test both
	332	Int_t upper = (Int_t)((z < fNRows-1) ? fPositions[z+1] : fN); // ending Position of the bin
	333	if((position + 1) < (upper)){
	334	if(TMath::Abs(y - fClusters[position + 1]->GetY()) < TMath::Abs(y - fClusters[position]->GetY())) return position + 1;
	335	else return position;
	336	}
	337	return position;
338	}
339
340	//_____________________________________________________________________________
341	Int_t AliTRDstackLayer::SearchNearestCluster(Double_t y, Double_t z, Double_t maxroady, Double_t maxroadz) const
342	{
343	//
344	// Finds the nearest cluster from a given point in a defined range.
345	// Distance is determined in a 2D space by the 2-Norm.
346	//
347	// Parameters :
348	// y : y position of the reference point in tracking coordinates
349	// z : z reference bin.
350	// maxroady : maximum searching distance in y direction
351	// maxroadz : maximum searching distance in z direction
352	//
353	// Output
354	// Index of the nearest cluster in the StackLayer indexing (-1 if no cluster is found).
355	// Cluster can be accessed with the operator[] or GetCluster(Int_t index)
356	//
357	// Detail description
358	//
359	// The following steps are perfomed:
360	// 1. Get the expected z bins inside maxroadz.
361	// 2. For each z bin find nearest y cluster.
362	// 3. Select best candidate
363	//
364	Int_t index = -1;
365	// initial minimal distance will be represented as ellipse: semi-major = z-direction
366	// later 2-Norm will be used
367	// Float_t nExcentricity = TMath::Sqrt(maxroadzmaxroadz - maxroadmaxroad)/maxroadz;
368	Float_t mindist = maxroadz;
369
370	// not very nice but unfortunately neccessarry: we have ho check the neighbors in both directions (+ and -) too. How
371	// much neighbors depends on the Quotient maxroadz/fZLength
372	UChar_t maxRows = 3;
373	UChar_t zpos[kMaxRows];
374	// Float_t mindist = TMath::Sqrt(maxroadmaxroad + maxroadzmaxroadz);
375	// UChar_t myZbin = FindTreePosition(z, fZ0 + fZLength/2, fZLength/4, 8, 8, kFALSE);
376	UChar_t myZbin = (UChar_t)(TMath::Abs(fZ0 - z)/fZLength * fNRows);
377	if(z < fZ0) myZbin = 0;
378	if(z > fZ0 + fZLength) myZbin = fNRows - 1;
379
380	UChar_t nNeighbors = 0;
381	for(UChar_t i = 0; i < maxRows; i++){
382	if((myZbin - 1 + i) < 0) continue;
383	if((myZbin - 1 + i) > fNRows - 1) break;
384	zpos[nNeighbors] = myZbin - 1 + i;
385	nNeighbors++;
386	}
387	Float_t ycl = 0, zcl = 0;
388	for(UChar_t neighbor = 0; neighbor < nNeighbors; neighbor++){ // Always test the neighbors too
389	Int_t pos = FindNearestYCluster(y,zpos[neighbor]);
390	if(pos == -1) continue; // No cluster in bin
391	AliTRDcluster c = (AliTRDcluster ) (fClusters[pos]);
392	if(c->IsUsed()) continue; // we are only interested in unused clusters
393	ycl = c->GetY();
394	// Too far away in y-direction (Prearrangement)
395	if (TMath::Abs(ycl - y) > maxroady) continue;
396
397	zcl = c->GetZ();
398	// Too far away in z-Direction
399	// (Prearrangement since we have not so many bins to test)
400	if (TMath::Abs(zcl - z) > maxroadz) continue;
401
402	Float_t dist; // distance defined as 2-Norm
403	// if we havent found a Particle that is in the ellipse around (y,z) with maxroad as semi-minor and
404	// maxroadz as semi-major, we take the radius of the ellipse concerning the cluster as mindist, later we
405	// take the 2-Norm when we found a cluster inside the ellipse (The value 10000 is taken because it is surely
406	// large enough to be usable as an indicator whether we have found a nearer cluster or not)
407	// if(mindist > 10000.){
408	// Float_t phi = ((zcl - z) == 0) ? TMath::Pi()/2 : TMath::ATan((ycl - y)/(zcl - z));
409	// mindist = maxroad/TMath::Sqrt(1 - nExcentricitynExcentricity (TMath::Cos(phi))*(TMath::Cos(phi)));
410	// }
411	dist = TMath::Max(TMath::Abs(y-ycl),TMath::Abs(z-zcl)); // infinity Norm
412	// dist = TMath::Sqrt((ycl - y)(ycl - y) + (zcl - z)(zcl - z));
413	if((Int_t)(dist * 100000) < (Int_t)(mindist * 100000)){
414	//if((dist = TMath::Sqrt((ycl - y)(ycl - y) + (zcl - z)(zcl - z))) < mindist){
415	mindist = dist;
416	index = pos;
417	}
418	}
419	// This is the Array Position in fIndex2D of the Nearest cluster: if a
420	// cluster is called, then the function has to retrieve the Information
421	// which is Stored in the Array called, the function
422	return index;
423	}
424
425	//_____________________________________________________________________________
426	void AliTRDstackLayer::BuildCond(AliTRDcluster cl, Double_t cond, UChar_t Layer, Double_t theta, Double_t phi)
427	{
428	// Helper function to calculate the area where to expect a cluster in THIS
429	// layer.
430	//
431	// Parameters :
432	// cl :
433	// cond :
434	// Layer :
435	// theta :
436	// phi :
437	//
438	// Detail description
439	//
440	// Helper function to calculate the area where to expect a cluster in THIS
441	// layer. by using the information of a former cluster in another layer
442	// and the angle in theta- and phi-direction between layer 0 and layer 3.
443	// If the layer is zero, initial conditions are calculated. Otherwise a
444	// linear interpolation is performed.
445	//Begin_Html
446	//<img src="gif/build_cond.gif">
447	//End_Html
448	//
449
450	if(!fRecoParam){
451	AliError("Reconstruction parameters not initialized.");
452	return;
453	}
454
455	if(Layer == 0){
456	cond[0] = cl->GetY(); // center: y-Direction
457	cond[1] = cl->GetZ(); // center: z-Direction
458	cond[2] = fRecoParam->GetMaxPhi() * (cl->GetX() - GetX()) + 1.0; // deviation: y-Direction
459	cond[3] = fRecoParam->GetMaxTheta() * (cl->GetX() - GetX()) + 1.0; // deviation: z-Direction
460	} else {
461	cond[0] = cl->GetY() + phi * (GetX() - cl->GetX());
462	cond[1] = cl->GetZ() + theta * (GetX() - cl->GetX());
463	cond[2] = fRecoParam->GetRoad0y() + phi;
464	cond[3] = fRecoParam->GetRoad0z();
465	}
466	}
467
468	//_____________________________________________________________________________
469	void AliTRDstackLayer::GetClusters(Double_t cond, Int_t index, Int_t& ncl, Int_t BufferSize)
470	{
471	// Finds all clusters situated in this layer inside a rectangle given by the center an ranges.
472	//
473	// Parameters :
474	// cond :
475	// index :
476	// ncl :
477	// BufferSize :
478	//
479	// Output :
480	//
481	// Detail description
482	//
483	// Function returs an array containing the indices in the stacklayer of
484	// the clusters found an the number of found clusters in the stacklayer
485
486	ncl = 0;
487	memset(index, 0, BufferSize*sizeof(Int_t));
488	if(fN == 0) return;
489
490	//Boundary checks
491	Double_t zvals[2];
492	zvals[0] = ((cond[1] - cond[3]) < fZ0) ? fZ0 : (cond[1] - cond[3]);
493	zvals[1] = ((cond[1] + cond[3]) < fZ0 + fZLength) ? (cond[1] + cond[3]) : fZ0 + fZLength;
494
495	UChar_t zlo = (fZ0>zvals[0]) ? 0 : (UChar_t)(TMath::Abs(fZ0 - zvals[0])/fZLength * fNRows);
496	UChar_t zhi = (fZ0+fZLength<zvals[1]) ? fNRows - 1 : (UChar_t)(TMath::Abs(fZ0 - zvals[1])/fZLength * fNRows);
497
498	//Preordering in Direction z saves a lot of loops (boundary checked)
499	for(UChar_t z = zlo; z <= zhi; z++){
1e8d1d74	500	UInt_t upper = (z < fNRows-1) ? fPositions[z+1] : fN;
e4f2f73d	501	for(Int_t y = fPositions[z]; y < (Int_t)upper; y++){
	502	if(ncl == BufferSize){
	503	AliWarning("Buffer size riched. Some clusters may be lost.");
	504	return; //Buffer filled
	505	}
	506	if(fClusters[y]->GetY() > (cond[0] + cond[2])) break; // Abbortion conditions!!!
	507	if(fClusters[y]->GetY() < (cond[0] - cond[2])) continue; // Too small
	508	if(((Int_t)((fClusters[y]->GetZ())1000) < (Int_t)(zvals[0]1000)) \|\| ((Int_t)((fClusters[y]->GetZ())1000) > (Int_t)(zvals[1]1000))) continue;
	509	index[ncl] = y;
	510	ncl++;
	511	}
	512	}
	513	if(ncl>fN) AliError(Form("Clusters found %d > %d (clusters in layer)", ncl, fN));
	514	}
	515
	516	//_____________________________________________________________________________
	517	AliTRDcluster AliTRDstackLayer::GetNearestCluster(Double_t cond)
	518	{
	519	// Function returning a pointer to the nearest cluster (nullpointer if not successfull).
	520	//
	521	// Parameters :
	522	// cond :
	523	//
	524	// Output :
	525	// pointer to the nearest cluster (nullpointer if not successfull).
	526	//
	527	// Detail description
	528	//
	529	// returns a pointer to the nearest cluster (nullpointer if not
	530	// successfull) by the help of the method FindNearestCluster
	531
	532
	533	Double_t maxroad = fRecoParam->GetRoad2y();
	534	Double_t maxroadz = fRecoParam->GetRoad2z();
	535
	536	Int_t index = SearchNearestCluster(cond[0],cond[1],maxroad,maxroadz);
	537	AliTRDcluster *returnCluster = 0x0;
	538	if(index != -1) returnCluster = (AliTRDcluster *) fClusters[index];
	539	return returnCluster;
	540	}
	541
	542	//_____________________________________________________________________________
	543	void AliTRDstackLayer::PrintClusters() const
	544	{
	545	// Prints the position of each cluster in the stacklayer on the stdout
	546	//
669ade02	547	//printf("fDebugStream = %#o\n", ((Int_t) fDebugStream));
669ade02	548	//printf("fRecoParam = %#o\n", ((Int_t) fRecoParam));
e4f2f73d	549
	550	for(Int_t i = 0; i < fN; i++){
	551	printf("AliTRDstackLayer: index=%i, Cluster: X = %3.3f, Y = %3.3f, Z = %3.3f\n", i, fClusters[i]->GetX(),fClusters[i]->GetY(),fClusters[i]->GetZ());
	552	if(fClusters[i]->IsUsed()) printf("cluster allready used. rejected in search algorithm\n");
	553	}
	554	}