[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.

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