[u/mrichter/AliRoot.git] / TRD / AliTRDtrackingChamber.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: AliTRDtrackingDebug.cxx 23810 2008-02-08 09:00:27Z hristov $ */

////////////////////////////////////////////////////////////////////////////
//                                                                        //
//  Tracking in one chamber                                               //
//                                                                        //
//  Authors:                                                              //
//    Alex Bercuci <A.Bercuci@gsi.de>                                     //
//    Markus Fasel <M.Fasel@gsi.de>                                       //
//                                                                        //
////////////////////////////////////////////////////////////////////////////

#include "AliTRDtrackingChamber.h"

#include "TMath.h"
#include "TMatrixTBase.h"
#include <TTreeStream.h>

#include "AliTRDReconstructor.h"
#include "AliTRDrecoParam.h"
#include "AliTRDtrackerV1.h"
#include "AliTRDgeometry.h"
#include "AliTRDpadPlane.h"
#include "AliTRDcalibDB.h"

ClassImp(AliTRDtrackingChamber)

//_______________________________________________________
AliTRDtrackingChamber::AliTRDtrackingChamber(Int_t det) :
	fDetector(det)
	,fX0(0.)
{}  

//_______________________________________________________
void AliTRDtrackingChamber::Clear(const Option_t *opt)
{
	for(Int_t itb=0; itb<kNTimeBins; itb++) fTB[itb].Clear(opt);
}

//_______________________________________________________
void AliTRDtrackingChamber::InsertCluster(AliTRDcluster *c, Int_t index)
{
	fTB[c->GetLocalTimeBin()].InsertCluster(c, index);
}

//_______________________________________________________
Bool_t AliTRDtrackingChamber::Build(AliTRDgeometry *geo)
{
// Init chamber and all time bins (AliTRDchamberTimeBin)
// Calculates radial position of the chamber based on 
// radial positions of the time bins (calibration/alignment aware)
//
	Int_t stack = geo->GetChamber(fDetector);
	Int_t plane = geo->GetPlane(fDetector);
	AliTRDpadPlane *pp = geo->GetPadPlane(plane, stack);
	Double_t zl = pp->GetRow0ROC() - pp->GetRowEndROC();
	Double_t z0 = geo->GetRow0(plane, stack, 0) - zl;
	Int_t nrows = pp->GetNrows();
	
	Int_t index[50], jtb = 0;
	for(Int_t itb=0; itb<kNTimeBins; itb++){ 
		if(!fTB[itb]) continue;
		fTB[itb].SetRange(z0, zl);
		fTB[itb].SetNRows(nrows);
		fTB[itb].BuildIndices();
		index[jtb++] = itb;
	}	
	if(jtb<2) return kFALSE;
	
	
	// ESTIMATE POSITION OF PAD PLANE FOR THIS CHAMBER
	Double_t x0 = fTB[index[0]].GetX();
	Double_t x1 = fTB[index[1]].GetX();
	Double_t dx = (x0 - x1)/(index[1] - index[0]); 
	fX0 = x0 + dx*(index[0] - (Int_t)AliTRDcalibDB::Instance()->GetT0Average(fDetector));	
	return kTRUE;
}
	
//_______________________________________________________	
Int_t AliTRDtrackingChamber::GetNClusters() const
{
// Returns number of clusters in chamber
//
	Int_t n = 0;
	for(Int_t itb=0; itb<kNTimeBins; itb++){ 
		n += Int_t(fTB[itb]);
	}
	return n;	
}	

//_______________________________________________________
Double_t AliTRDtrackingChamber::GetQuality()
{
  //
  // Calculate chamber quality for seeding.
  // 
  //
  // Parameters :
  //   layers : Array of propagation layers for this plane.
  //
  // Output :
  //   plane quality factor for seeding
  // 
  // Detailed description
  //
  // The quality of the plane for seeding is higher if:
  //  1. the average timebin population is closer to an integer number
  //  2. the distribution of clusters/timebin is closer to a uniform distribution.
  //    - the slope of the first derivative of a parabolic fit is small or
  //    - the slope of a linear fit is small
  //

	Int_t ncl   = 0;
	Int_t nused = 0;
	Int_t nClLayer;
	for(int itb=0; itb<kNTimeBins; itb++){
		if(!(nClLayer = fTB[itb].GetNClusters())) continue;
		ncl += nClLayer;
		for(Int_t incl = 0; incl < nClLayer; incl++){
			if((fTB[itb].GetCluster(incl))->IsUsed()) nused++;
		}
	}
	
	// calculate the deviation of the mean number of clusters from the
	// closest integer values
	Float_t nclMed = float(ncl-nused)/AliTRDtrackerV1::GetNTimeBins();
	Int_t ncli = Int_t(nclMed);
	Float_t nclDev = TMath::Abs(nclMed - TMath::Max(ncli, 1));
	nclDev -= (nclDev>.5) && ncli ? 1. : 0.;
	return TMath::Exp(-5.*TMath::Abs(nclDev));

// 	// get slope of the derivative
// 	if(!fitter.Eval()) return quality;
// 	fitter.PrintResults(3);
// 	Double_t a = fitter.GetParameter(1);
// 
// 	printf("ncl_dev(%f)  a(%f)\n", ncl_dev, a);
// 	return quality*TMath::Exp(-a);

}


//_______________________________________________________
AliTRDchamberTimeBin *AliTRDtrackingChamber::GetSeedingLayer(AliTRDgeometry *geo)
{
  //
  // Creates a seeding layer
  //
	
	// constants
	const Int_t kMaxRows = 16;
	const Int_t kMaxCols = 144;
	const Int_t kMaxPads = 2304;
		
	// Get the geometrical data of the chamber
	Int_t plane = geo->GetPlane(fDetector);
	Int_t stack = geo->GetChamber(fDetector);
	Int_t sector= geo->GetSector(fDetector);
	AliTRDpadPlane *pp = geo->GetPadPlane(plane, stack);
	Int_t nCols = pp->GetNcols();
	Float_t ymin = TMath::Min(pp->GetCol0(), pp->GetColEnd());
	Float_t ymax = TMath::Max(pp->GetCol0(), pp->GetColEnd());
	Float_t zmin = TMath::Min(pp->GetRow0(), pp->GetRowEnd());
	Float_t zmax = TMath::Max(pp->GetRow0(), pp->GetRowEnd());
	Float_t z0 = -1., zl = -1.;
	Int_t nRows = pp->GetNrows();
	Float_t binlength = (ymax - ymin)/nCols; 
	//AliInfo(Form("ymin(%f) ymax(%f) zmin(%f) zmax(%f) nRows(%d) binlength(%f)", ymin, ymax, zmin, zmax, nRows, binlength));
	
	// Fill the histogram
	Int_t nClusters;	
	Int_t *histogram[kMaxRows];											// 2D-Histogram
	Int_t hvals[kMaxPads];	memset(hvals, 0, sizeof(Int_t)*kMaxPads);	
	Float_t *sigmas[kMaxRows];
	Float_t svals[kMaxPads];	memset(svals, 0, sizeof(Float_t)*kMaxPads);	
	AliTRDcluster *c = 0x0;
	for(Int_t irs = 0; irs < kMaxRows; irs++){
		histogram[irs] = &hvals[irs*kMaxCols];
		sigmas[irs] = &svals[irs*kMaxCols];
	}
	for(Int_t iTime = 0; iTime < kNTimeBins; iTime++){
		if(!(nClusters = fTB[iTime].GetNClusters())) continue;
		z0 = fTB[iTime].GetZ0();
		zl = fTB[iTime].GetDZ0();
		for(Int_t incl = 0; incl < nClusters; incl++){
			c = fTB[iTime].GetCluster(incl);	
			histogram[c->GetPadRow()][c->GetPadCol()]++;
			sigmas[c->GetPadRow()][c->GetPadCol()] += c->GetSigmaZ2();
		}
	}
	
// Now I have everything in the histogram, do the selection
	//Int_t nPads = nCols * nRows;
	// This is what we are interested in: The center of gravity of the best candidates
	Float_t cogyvals[kMaxPads]; memset(cogyvals, 0, sizeof(Float_t)*kMaxPads);
	Float_t cogzvals[kMaxPads]; memset(cogzvals, 0, sizeof(Float_t)*kMaxPads);
	Float_t *cogy[kMaxRows];
	Float_t *cogz[kMaxRows];
	
	// Lookup-Table storing coordinates according to the bins
	Float_t yLengths[kMaxCols];
	Float_t zLengths[kMaxRows];
	for(Int_t icnt = 0; icnt < nCols; icnt++){
 		yLengths[icnt] = pp->GetColPos(nCols - 1 - icnt) + binlength/2;
	}
	for(Int_t icnt = 0; icnt < nRows; icnt++){
		zLengths[icnt] = pp->GetRowPos(icnt) - pp->GetRowSize(icnt)/2;
	}

	// A bitfield is used to mask the pads as usable
	Short_t mask[kMaxCols]; memset(mask, 0 ,sizeof(Short_t) * kMaxCols);//bool mvals[kMaxPads];
	for(UChar_t icount = 0; icount < nRows; icount++){
		cogy[icount] = &cogyvals[icount*kMaxCols];
		cogz[icount] = &cogzvals[icount*kMaxCols];
	}
	// In this array the array position of the best candidates will be stored
	Int_t   cand[AliTRDtrackerV1::kMaxTracksStack];
	Float_t sigcands[AliTRDtrackerV1::kMaxTracksStack];
 	
	// helper variables
	Int_t indices[kMaxPads]; memset(indices, -1, sizeof(Int_t)*kMaxPads);
	Int_t nCandidates = 0;
 	Float_t norm, cogv;
	// histogram filled -> Select best bins
	Int_t nPads = nCols * nRows;
	TMath::Sort(nPads, hvals, indices);			// bins storing a 0 should not matter
	// Set Threshold
	Int_t maximum = hvals[indices[0]];	// best
	Int_t threshold = Int_t(maximum * AliTRDReconstructor::RecoParam()->GetFindableClusters());
	Int_t col, row, lower, lower1, upper, upper1;
	for(Int_t ib = 0; ib < nPads; ib++){
		if(nCandidates >= AliTRDtrackerV1::kMaxTracksStack){
			printf("Number of seed candidates %d exceeded maximum allowed per stack %d", nCandidates, AliTRDtrackerV1::kMaxTracksStack);
			break;
		}
		// Positions
		row = indices[ib]/nCols;
		col = indices[ib]%nCols;
		// here will be the threshold condition:
		if((mask[col] & (1 << row)) != 0) continue;		// Pad is masked: continue
		if(histogram[row][col] < TMath::Max(threshold, 1)){	// of course at least one cluster is needed
			break;			// number of clusters below threshold: break;
		} 
		// passing: Mark the neighbors
		lower  = TMath::Max(col - 1, 0); upper  = TMath::Min(col + 2, nCols);
		lower1 = TMath::Max(row - 1, 0); upper1 = TMath::Min(row + 2, nCols);
		for(Int_t ic = lower; ic < upper; ++ic)
			for(Int_t ir = lower1; ir < upper1; ++ir){
				if(ic == col && ir == row) continue;
				mask[ic] |= (1 << ir);
			}
		// Storing the position in an array
		// testing for neigboring
		cogv = 0;
		norm = 0;
		lower = TMath::Max(col - 1, 0);
		upper = TMath::Min(col + 2, nCols);
		for(Int_t inb = lower; inb < upper; ++inb){
			cogv += yLengths[inb] * histogram[row][inb];
			norm += histogram[row][inb];
		}
		cogy[row][col] = cogv / norm;
		cogv = 0; norm = 0;
		lower = TMath::Max(row - 1, 0);
		upper = TMath::Min(row + 2, nRows);
		for(Int_t inb = lower; inb < upper; ++inb){
			cogv += zLengths[inb] * histogram[inb][col];
			norm += histogram[inb][col];
		}
		cogz[row][col] = Float_t(cogv) /  norm;
		// passed the filter
		cand[nCandidates] = row*nCols + col;	// store the position of a passig candidate into an Array
		sigcands[nCandidates] = sigmas[row][col] / histogram[row][col]; // never be a floating point exeption
		// Analysis output
		nCandidates++;
	}
	if(!nCandidates) return 0x0;
	
	Float_t pos[3], sig[2];
	Short_t signal[7]; memset(&signal[0], 0, 7*sizeof(Short_t));
	AliTRDchamberTimeBin *fakeLayer = new AliTRDchamberTimeBin(plane, stack, sector, z0, zl);
	AliTRDcluster *cluster = 0x0;
	if(nCandidates){
		UInt_t fakeIndex = 0;
		for(Int_t ican = 0; ican < nCandidates; ican++){
			row = cand[ican] / nCols;
			col = cand[ican] % nCols;
			//temporary
			Int_t n = 0; Double_t x = 0., y = 0., z = 0.;
			for(int itb=0; itb<kNTimeBins; itb++){
				if(!(nClusters = fTB[itb].GetNClusters())) continue;
				for(Int_t incl = 0; incl < nClusters; incl++){
					c = fTB[itb].GetCluster(incl);	
					if(c->GetPadRow() != row) continue;
					if(TMath::Abs(c->GetPadCol() - col) > 2) continue;
					x += c->GetX();
					y += c->GetY();
					z += c->GetZ();
					n++;
				}
			}
			pos[0] = x/n;
			pos[1] = y/n;
			pos[2] = z/n;
			sig[0] = .02;
			sig[1] = sigcands[ican];
			cluster = new AliTRDcluster(fDetector, 0., pos, sig, 0x0, 3, signal, col, row, 0, 0, 0., 0);
			fakeLayer->InsertCluster(cluster, fakeIndex++);
		}
	}
	fakeLayer->SetNRows(nRows);
	fakeLayer->SetOwner();
	fakeLayer->BuildIndices();
	//fakeLayer->PrintClusters();
	
	if(AliTRDReconstructor::StreamLevel() >= 3){
		//TMatrixD hist(nRows, nCols);
		//for(Int_t i = 0; i < nRows; i++)
		//	for(Int_t j = 0; j < nCols; j++)
		//		hist(i,j) = histogram[i][j];
		TTreeSRedirector &cstreamer = *AliTRDtrackerV1::DebugStreamer();
		cstreamer << "GetSeedingLayer"
		<< "plane="      << plane
		<< "ymin="       << ymin
		<< "ymax="       << ymax
		<< "zmin="       << zmin
		<< "zmax="       << zmax
		<< "L.="         << fakeLayer
		//<< "Histogram.=" << &hist
		<< "\n";
	}
	
	return fakeLayer;
}
Commit	Line	Data
eb38ed55	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: AliTRDtrackingDebug.cxx 23810 2008-02-08 09:00:27Z hristov $ */
	17
	18	////////////////////////////////////////////////////////////////////////////
	19	// //
	20	// Tracking in one chamber //
	21	// //
	22	// Authors: //
	23	// Alex Bercuci <A.Bercuci@gsi.de> //
	24	// Markus Fasel <M.Fasel@gsi.de> //
	25	// //
	26	////////////////////////////////////////////////////////////////////////////
	27
	28	#include "AliTRDtrackingChamber.h"
	29
	30	#include "TMath.h"
	31	#include "TMatrixTBase.h"
	32	#include <TTreeStream.h>
	33
	34	#include "AliTRDReconstructor.h"
	35	#include "AliTRDrecoParam.h"
	36	#include "AliTRDtrackerV1.h"
	37	#include "AliTRDgeometry.h"
	38	#include "AliTRDpadPlane.h"
	39	#include "AliTRDcalibDB.h"
	40
	41	ClassImp(AliTRDtrackingChamber)
	42
	43	//_______________________________________________________
	44	AliTRDtrackingChamber::AliTRDtrackingChamber(Int_t det) :
	45	fDetector(det)
	46	,fX0(0.)
	47	{}
	48
	49	//_______________________________________________________
	50	void AliTRDtrackingChamber::Clear(const Option_t *opt)
	51	{
	52	for(Int_t itb=0; itb<kNTimeBins; itb++) fTB[itb].Clear(opt);
	53	}
	54
	55	//_______________________________________________________
	56	void AliTRDtrackingChamber::InsertCluster(AliTRDcluster *c, Int_t index)
	57	{
	58	fTB[c->GetLocalTimeBin()].InsertCluster(c, index);
	59	}
	60
	61	//_______________________________________________________
	62	Bool_t AliTRDtrackingChamber::Build(AliTRDgeometry *geo)
	63	{
	64	// Init chamber and all time bins (AliTRDchamberTimeBin)
65	// Calculates radial position of the chamber based on
66	// radial positions of the time bins (calibration/alignment aware)
67	//
68	Int_t stack = geo->GetChamber(fDetector);
69	Int_t plane = geo->GetPlane(fDetector);
70	AliTRDpadPlane *pp = geo->GetPadPlane(plane, stack);
71	Double_t zl = pp->GetRow0ROC() - pp->GetRowEndROC();
72	Double_t z0 = geo->GetRow0(plane, stack, 0) - zl;
73	Int_t nrows = pp->GetNrows();
74
75	Int_t index[50], jtb = 0;
76	for(Int_t itb=0; itb<kNTimeBins; itb++){
77	if(!fTB[itb]) continue;
78	fTB[itb].SetRange(z0, zl);
79	fTB[itb].SetNRows(nrows);
80	fTB[itb].BuildIndices();
81	index[jtb++] = itb;
82	}
83	if(jtb<2) return kFALSE;
84
85
86	// ESTIMATE POSITION OF PAD PLANE FOR THIS CHAMBER
87	Double_t x0 = fTB[index[0]].GetX();
88	Double_t x1 = fTB[index[1]].GetX();
89	Double_t dx = (x0 - x1)/(index[1] - index[0]);
90	fX0 = x0 + dx*(index[0] - (Int_t)AliTRDcalibDB::Instance()->GetT0Average(fDetector));
91	return kTRUE;
92	}
93
94	//_______________________________________________________
95	Int_t AliTRDtrackingChamber::GetNClusters() const
96	{
97	// Returns number of clusters in chamber
98	//
99	Int_t n = 0;
100	for(Int_t itb=0; itb<kNTimeBins; itb++){
101	n += Int_t(fTB[itb]);
102	}
103	return n;
104	}
105
106	//_______________________________________________________
2985ffcb	107	Double_t AliTRDtrackingChamber::GetQuality()
eb38ed55	108	{
	109	//
	110	// Calculate chamber quality for seeding.
	111	//
	112	//
	113	// Parameters :
	114	// layers : Array of propagation layers for this plane.
	115	//
	116	// Output :
	117	// plane quality factor for seeding
	118	//
	119	// Detailed description
	120	//
	121	// The quality of the plane for seeding is higher if:
	122	// 1. the average timebin population is closer to an integer number
	123	// 2. the distribution of clusters/timebin is closer to a uniform distribution.
	124	// - the slope of the first derivative of a parabolic fit is small or
	125	// - the slope of a linear fit is small
	126	//
	127
	128	Int_t ncl = 0;
	129	Int_t nused = 0;
	130	Int_t nClLayer;
	131	for(int itb=0; itb<kNTimeBins; itb++){
	132	if(!(nClLayer = fTB[itb].GetNClusters())) continue;
	133	ncl += nClLayer;
	134	for(Int_t incl = 0; incl < nClLayer; incl++){
	135	if((fTB[itb].GetCluster(incl))->IsUsed()) nused++;
	136	}
	137	}
	138
	139	// calculate the deviation of the mean number of clusters from the
	140	// closest integer values
2985ffcb	141	Float_t nclMed = float(ncl-nused)/AliTRDtrackerV1::GetNTimeBins();
eb38ed55	142	Int_t ncli = Int_t(nclMed);
	143	Float_t nclDev = TMath::Abs(nclMed - TMath::Max(ncli, 1));
	144	nclDev -= (nclDev>.5) && ncli ? 1. : 0.;
	145	return TMath::Exp(-5.*TMath::Abs(nclDev));
	146
	147	// // get slope of the derivative
	148	// if(!fitter.Eval()) return quality;
	149	// fitter.PrintResults(3);
	150	// Double_t a = fitter.GetParameter(1);
	151	//
	152	// printf("ncl_dev(%f) a(%f)\n", ncl_dev, a);
	153	// return quality*TMath::Exp(-a);
	154
	155	}
	156
	157
	158	//_______________________________________________________
	159	AliTRDchamberTimeBin AliTRDtrackingChamber::GetSeedingLayer(AliTRDgeometry geo)
	160	{
	161	//
	162	// Creates a seeding layer
	163	//
	164
	165	// constants
	166	const Int_t kMaxRows = 16;
	167	const Int_t kMaxCols = 144;
	168	const Int_t kMaxPads = 2304;
	169
	170	// Get the geometrical data of the chamber
	171	Int_t plane = geo->GetPlane(fDetector);
	172	Int_t stack = geo->GetChamber(fDetector);
	173	Int_t sector= geo->GetSector(fDetector);
	174	AliTRDpadPlane *pp = geo->GetPadPlane(plane, stack);
	175	Int_t nCols = pp->GetNcols();
	176	Float_t ymin = TMath::Min(pp->GetCol0(), pp->GetColEnd());
	177	Float_t ymax = TMath::Max(pp->GetCol0(), pp->GetColEnd());
	178	Float_t zmin = TMath::Min(pp->GetRow0(), pp->GetRowEnd());
	179	Float_t zmax = TMath::Max(pp->GetRow0(), pp->GetRowEnd());
	180	Float_t z0 = -1., zl = -1.;
	181	Int_t nRows = pp->GetNrows();
	182	Float_t binlength = (ymax - ymin)/nCols;
	183	//AliInfo(Form("ymin(%f) ymax(%f) zmin(%f) zmax(%f) nRows(%d) binlength(%f)", ymin, ymax, zmin, zmax, nRows, binlength));
	184
	185	// Fill the histogram
	186	Int_t nClusters;
	187	Int_t *histogram[kMaxRows]; // 2D-Histogram
	188	Int_t hvals[kMaxPads]; memset(hvals, 0, sizeof(Int_t)*kMaxPads);
	189	Float_t *sigmas[kMaxRows];
	190	Float_t svals[kMaxPads]; memset(svals, 0, sizeof(Float_t)*kMaxPads);
	191	AliTRDcluster *c = 0x0;
	192	for(Int_t irs = 0; irs < kMaxRows; irs++){
	193	histogram[irs] = &hvals[irs*kMaxCols];
	194	sigmas[irs] = &svals[irs*kMaxCols];
	195	}
	196	for(Int_t iTime = 0; iTime < kNTimeBins; iTime++){
	197	if(!(nClusters = fTB[iTime].GetNClusters())) continue;
	198	z0 = fTB[iTime].GetZ0();
	199	zl = fTB[iTime].GetDZ0();
	200	for(Int_t incl = 0; incl < nClusters; incl++){
	201	c = fTB[iTime].GetCluster(incl);
	202	histogram[c->GetPadRow()][c->GetPadCol()]++;
	203	sigmas[c->GetPadRow()][c->GetPadCol()] += c->GetSigmaZ2();
	204	}
	205	}
206
207	// Now I have everything in the histogram, do the selection
208	//Int_t nPads = nCols * nRows;
209	// This is what we are interested in: The center of gravity of the best candidates
210	Float_t cogyvals[kMaxPads]; memset(cogyvals, 0, sizeof(Float_t)*kMaxPads);
211	Float_t cogzvals[kMaxPads]; memset(cogzvals, 0, sizeof(Float_t)*kMaxPads);
212	Float_t *cogy[kMaxRows];
213	Float_t *cogz[kMaxRows];
214
215	// Lookup-Table storing coordinates according to the bins
216	Float_t yLengths[kMaxCols];
217	Float_t zLengths[kMaxRows];
218	for(Int_t icnt = 0; icnt < nCols; icnt++){
219	yLengths[icnt] = pp->GetColPos(nCols - 1 - icnt) + binlength/2;
220	}
221	for(Int_t icnt = 0; icnt < nRows; icnt++){
222	zLengths[icnt] = pp->GetRowPos(icnt) - pp->GetRowSize(icnt)/2;
223	}
224
225	// A bitfield is used to mask the pads as usable
226	Short_t mask[kMaxCols]; memset(mask, 0 ,sizeof(Short_t) * kMaxCols);//bool mvals[kMaxPads];
227	for(UChar_t icount = 0; icount < nRows; icount++){
228	cogy[icount] = &cogyvals[icount*kMaxCols];
229	cogz[icount] = &cogzvals[icount*kMaxCols];
230	}
231	// In this array the array position of the best candidates will be stored
232	Int_t cand[AliTRDtrackerV1::kMaxTracksStack];
233	Float_t sigcands[AliTRDtrackerV1::kMaxTracksStack];
234
235	// helper variables
236	Int_t indices[kMaxPads]; memset(indices, -1, sizeof(Int_t)*kMaxPads);
237	Int_t nCandidates = 0;
238	Float_t norm, cogv;
239	// histogram filled -> Select best bins
240	Int_t nPads = nCols * nRows;
241	TMath::Sort(nPads, hvals, indices); // bins storing a 0 should not matter
242	// Set Threshold
243	Int_t maximum = hvals[indices[0]]; // best
244	Int_t threshold = Int_t(maximum * AliTRDReconstructor::RecoParam()->GetFindableClusters());
245	Int_t col, row, lower, lower1, upper, upper1;
246	for(Int_t ib = 0; ib < nPads; ib++){
247	if(nCandidates >= AliTRDtrackerV1::kMaxTracksStack){
248	printf("Number of seed candidates %d exceeded maximum allowed per stack %d", nCandidates, AliTRDtrackerV1::kMaxTracksStack);
249	break;
250	}
251	// Positions
252	row = indices[ib]/nCols;
253	col = indices[ib]%nCols;
254	// here will be the threshold condition:
255	if((mask[col] & (1 << row)) != 0) continue; // Pad is masked: continue
256	if(histogram[row][col] < TMath::Max(threshold, 1)){ // of course at least one cluster is needed
257	break; // number of clusters below threshold: break;
258	}
259	// passing: Mark the neighbors
260	lower = TMath::Max(col - 1, 0); upper = TMath::Min(col + 2, nCols);
261	lower1 = TMath::Max(row - 1, 0); upper1 = TMath::Min(row + 2, nCols);
262	for(Int_t ic = lower; ic < upper; ++ic)
263	for(Int_t ir = lower1; ir < upper1; ++ir){
264	if(ic == col && ir == row) continue;
265	mask[ic] \|= (1 << ir);
266	}
267	// Storing the position in an array
268	// testing for neigboring
269	cogv = 0;
270	norm = 0;
271	lower = TMath::Max(col - 1, 0);
272	upper = TMath::Min(col + 2, nCols);
273	for(Int_t inb = lower; inb < upper; ++inb){
274	cogv += yLengths[inb] * histogram[row][inb];
275	norm += histogram[row][inb];
276	}
277	cogy[row][col] = cogv / norm;
278	cogv = 0; norm = 0;
279	lower = TMath::Max(row - 1, 0);
280	upper = TMath::Min(row + 2, nRows);
281	for(Int_t inb = lower; inb < upper; ++inb){
282	cogv += zLengths[inb] * histogram[inb][col];
283	norm += histogram[inb][col];
284	}
285	cogz[row][col] = Float_t(cogv) / norm;
286	// passed the filter
287	cand[nCandidates] = row*nCols + col; // store the position of a passig candidate into an Array
288	sigcands[nCandidates] = sigmas[row][col] / histogram[row][col]; // never be a floating point exeption
289	// Analysis output
290	nCandidates++;
291	}
292	if(!nCandidates) return 0x0;
293
294	Float_t pos[3], sig[2];
295	Short_t signal[7]; memset(&signal[0], 0, 7*sizeof(Short_t));
296	AliTRDchamberTimeBin *fakeLayer = new AliTRDchamberTimeBin(plane, stack, sector, z0, zl);
297	AliTRDcluster *cluster = 0x0;
298	if(nCandidates){
299	UInt_t fakeIndex = 0;
300	for(Int_t ican = 0; ican < nCandidates; ican++){
301	row = cand[ican] / nCols;
302	col = cand[ican] % nCols;
303	//temporary
304	Int_t n = 0; Double_t x = 0., y = 0., z = 0.;
305	for(int itb=0; itb<kNTimeBins; itb++){
306	if(!(nClusters = fTB[itb].GetNClusters())) continue;
307	for(Int_t incl = 0; incl < nClusters; incl++){
308	c = fTB[itb].GetCluster(incl);
309	if(c->GetPadRow() != row) continue;
310	if(TMath::Abs(c->GetPadCol() - col) > 2) continue;
311	x += c->GetX();
312	y += c->GetY();
313	z += c->GetZ();
314	n++;
315	}
316	}
317	pos[0] = x/n;
318	pos[1] = y/n;
319	pos[2] = z/n;
320	sig[0] = .02;
321	sig[1] = sigcands[ican];
322	cluster = new AliTRDcluster(fDetector, 0., pos, sig, 0x0, 3, signal, col, row, 0, 0, 0., 0);
323	fakeLayer->InsertCluster(cluster, fakeIndex++);
324	}
325	}
326	fakeLayer->SetNRows(nRows);
327	fakeLayer->SetOwner();
328	fakeLayer->BuildIndices();
329	//fakeLayer->PrintClusters();
330
331	if(AliTRDReconstructor::StreamLevel() >= 3){
332	//TMatrixD hist(nRows, nCols);
333	//for(Int_t i = 0; i < nRows; i++)
334	// for(Int_t j = 0; j < nCols; j++)
335	// hist(i,j) = histogram[i][j];
336	TTreeSRedirector &cstreamer = *AliTRDtrackerV1::DebugStreamer();
337	cstreamer << "GetSeedingLayer"
338	<< "plane=" << plane
339	<< "ymin=" << ymin
340	<< "ymax=" << ymax
341	<< "zmin=" << zmin
342	<< "zmax=" << zmax
343	<< "L.=" << fakeLayer
344	//<< "Histogram.=" << &hist
345	<< "\n";
346	}
347
348	return fakeLayer;
349	}
350