ClassImp(AliTRDReconstructor)
-Bool_t AliTRDReconstructor::fgkSeedingOn = kFALSE;
-Int_t AliTRDReconstructor::fgStreamLevel = 0; // Stream (debug) level
AliTRDrecoParam* AliTRDReconstructor::fgRecoParam = 0x0;
-
//_____________________________________________________________________________
AliTRDReconstructor::~AliTRDReconstructor() {
if(fgRecoParam) delete fgRecoParam;
, esd); }
virtual void FillESD(TTree *digitsTree, TTree *clusterTree, AliESDEvent *esd) const;
- static void SetSeedingOn(Bool_t seeding) { fgkSeedingOn = seeding; }
- static void SetStreamLevel(Int_t level) { fgStreamLevel = level; }
static void SetRecoParam(AliTRDrecoParam *reco) { fgRecoParam = reco; }
- static Bool_t SeedingOn() { return fgkSeedingOn; }
- static Int_t StreamLevel() { return fgStreamLevel; }
private:
- static Bool_t fgkSeedingOn; // Set flag for seeding during reconstruction
- static Int_t fgStreamLevel; // Flag for streaming
static AliTRDrecoParam *fgRecoParam; // Reconstruction parameters
ClassDef(AliTRDReconstructor,0) // Class for the TRD reconstruction
-/**************************************************************************
- * 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: AliTRDchamberTimeBin.cxx 23313 2008-01-11 14:56:43Z cblume $ */
-
-///////////////////////////////////////////////////////////////////////////////
-// //
-// 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 "AliTRDchamberTimeBin.h"
-#include "AliTRDrecoParam.h"
-#include "AliTRDReconstructor.h"
-#include "AliTRDtrackerV1.h"
-
-
-ClassImp(AliTRDchamberTimeBin)
-
-//_____________________________________________________________________________
-AliTRDchamberTimeBin::AliTRDchamberTimeBin(Int_t plane, Int_t stack, Int_t sector, Double_t z0, Double_t zLength)
- :TObject()
- ,fOwner(kFALSE)
- ,fPlane(plane)
- ,fStack(stack)
- ,fSector(sector)
- ,fNRows(kMaxRows)
- ,fN(0)
- ,fX(0.)
- ,fZ0(z0)
- ,fZLength(zLength)
-{
- //
- // Default constructor (Only provided to use AliTRDchamberTimeBin with arrays)
- //
-
- for(int i=0; i<kMaxRows; i++) fPositions[i] = 0xff;
- for(int ic=0; ic<kMaxClustersLayer; ic++){
- fClusters[ic] = 0x0;
- fIndex[ic] = 0xffff;
- }
-}
-
-// //_____________________________________________________________________________
-// AliTRDchamberTimeBin::AliTRDchamberTimeBin(const AliTRDpropagationLayer &layer, Double_t
-// z0, Double_t zLength, UChar_t stackNr):
-// TObject()
-// ,fOwner(kFALSE)
-// ,fPlane(0xff)
-// ,fStack(0xff)
-// ,fSector(0xff)
-// ,fNRows(kMaxRows)
-// ,fN(0)
-// ,fX(0.)
-// ,fZ0(z0)
-// ,fZLength(zLength)
-// {
-// // Standard constructor.
-// // Initialize also the underlying AliTRDpropagationLayer using the copy constructor.
-//
-// SetT0(layer.IsT0());
-// for(int i=0; i<kMaxRows; i++) fPositions[i] = 0xff;
-// for(int ic=0; ic<kMaxClustersLayer; ic++){
-// fClusters[ic] = 0x0;
-// fIndex[ic] = 0xffff;
-// }
-// }
-//
-// //_____________________________________________________________________________
-// AliTRDchamberTimeBin::AliTRDchamberTimeBin(const AliTRDpropagationLayer &layer):
-// TObject()
-// ,fOwner(kFALSE)
-// ,fPlane(0xff)
-// ,fStack(0xff)
-// ,fSector(0xff)
-// ,fNRows(kMaxRows)
-// ,fN(0)
-// ,fX(0.)
-// ,fZ0(0)
-// ,fZLength(0)
-// {
-// // Standard constructor using only AliTRDpropagationLayer.
-//
-// SetT0(layer.IsT0());
-// for(int i=0; i<kMaxRows; i++) fPositions[i] = 0xff;
-// for(int ic=0; ic<kMaxClustersLayer; ic++){
-// fClusters[ic] = 0x0;
-// fIndex[ic] = 0xffff;
-// }
-// }
-// //_____________________________________________________________________________
-// AliTRDchamberTimeBin &AliTRDchamberTimeBin::operator=(const AliTRDpropagationLayer &layer)
-// {
-// // Assignment operator from an AliTRDpropagationLayer
-//
-// if (this != &layer) layer.Copy(*this);
-// return *this;
-// }
-//
-
-//_____________________________________________________________________________
-AliTRDchamberTimeBin::AliTRDchamberTimeBin(const AliTRDchamberTimeBin &layer):
- TObject()
- ,fOwner(layer.fOwner)
- ,fPlane(layer.fPlane)
- ,fStack(layer.fStack)
- ,fSector(layer.fSector)
- ,fNRows(layer.fNRows)
- ,fN(layer.fN)
- ,fX(layer.fX)
- ,fZ0(layer.fZ0)
- ,fZLength(layer.fZLength)
-{
-// Copy Constructor (performs a deep copy)
-
- SetT0(layer.IsT0());
- for(int i=0; i<kMaxRows; i++) fPositions[i] = layer.fPositions[i];
- memcpy(&fClusters[0], &layer.fClusters[0], kMaxClustersLayer*sizeof(UChar_t));
- memcpy(&fIndex[0], &layer.fIndex[0], kMaxClustersLayer*sizeof(UInt_t));
-
-
-// BuildIndices();
-}
-
-//_____________________________________________________________________________
-AliTRDchamberTimeBin &AliTRDchamberTimeBin::operator=(const AliTRDchamberTimeBin &layer)
-{
-// Assignment operator
-
- if (this != &layer) layer.Copy(*this);
- return *this;
-}
-
-//_____________________________________________________________________________
-void AliTRDchamberTimeBin::Copy(TObject &o) const
-{
-// Copy method. Performs a deep copy of all data from this object to object o.
-
- AliTRDchamberTimeBin &layer = (AliTRDchamberTimeBin &)o;
- layer.fOwner = kFALSE;
- layer.fPlane = fPlane;
- layer.fStack = fStack;
- layer.fSector = fSector;
- layer.fNRows = fNRows;
- layer.fN = fN;
- layer.fX = fX;
- layer.fZ0 = fZ0;
- layer.fZLength = fZLength;
- layer.SetT0(IsT0());
-
- for(int i = 0; i < kMaxRows; i++) layer.fPositions[i] = 0;
-
- for(int i=0; i<kMaxRows; i++) layer.fPositions[i] = fPositions[i];
- memcpy(&layer.fClusters[0], &fClusters[0], kMaxClustersLayer*sizeof(UChar_t));
- memcpy(&layer.fIndex[0], &fIndex[0], kMaxClustersLayer*sizeof(UInt_t));
-
- TObject::Copy(layer); // copies everything into layer
-
-// layer.BuildIndices();
-}
-
-//_____________________________________________________________________________
-AliTRDchamberTimeBin::~AliTRDchamberTimeBin()
-{
-// Destructor
- if(fOwner) for(int ic=0; ic<fN; ic++) delete fClusters[ic];
-}
-
-//_____________________________________________________________________________
-void AliTRDchamberTimeBin::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 AliTRDchamberTimeBin::InsertCluster(AliTRDcluster *c, UInt_t index)
-{
- //
- // Insert cluster in cluster array.
- // Clusters are sorted according to Y coordinate.
- //
-
- //if (fTimeBinIndex < 0) {
- //AliWarning("Attempt to insert cluster into non-sensitive time bin!\n");
- //return;
- //}
-
- if (fN == (Int_t) kMaxClustersLayer) {
- //AliWarning("Too many clusters !\n");
- return;
- }
-
- if (fN == 0) {
- fIndex[0] = index;
- fClusters[fN++] = c;
- return;
- }
-
- Int_t i = Find(c->GetY());
- memmove(fClusters+i+1,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
- memmove(fIndex +i+1,fIndex +i,(fN-i)*sizeof(UInt_t));
- fIndex[i] = index;
- fClusters[i] = c;
- fN++;
-
-}
-
-
-//_____________________________________________________________________________
-void AliTRDchamberTimeBin::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++){
- if(fClusters[i]->IsUsed()){
- fClusters[i] = 0x0;
- fIndex[i] = 0xffff;
- } else nClStack++;
- }
- if(nClStack > kMaxClustersLayer) AliWarning(Form("Number of clusters in stack %d exceed buffer size %d. Truncating.", nClStack, kMaxClustersLayer));
-
- // Nothing in this time bin. Reset indexes
- if(!nClStack){
- fN = 0;
- memset(&fPositions[0], 0xff, sizeof(UChar_t) * kMaxRows);
- memset(&fClusters[0], 0x0, sizeof(AliTRDcluster*) * kMaxClustersLayer);
- memset(&fIndex[0], 0xffff, sizeof(UInt_t) * kMaxClustersLayer);
- 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] = 0xffff;
- nClStack++;
- }
-
- // do clusters arrangement
- fX = 0.;
- fN = nClStack;
- nClStack = 0;
- // Reset Positions array
- memset(fPositions, 0, sizeof(UChar_t)*kMaxRows);
- for(Int_t i = 0; i < fN; i++){
- // boundary check
- AliTRDcluster *cl = helpCL[i];
- UChar_t rowIndex = cl->GetPadRow();
- // Insert Leaf
- Int_t pos = FindYPosition(cl->GetY(), rowIndex, i);
- if(pos == -1){ // zbin is empty;
- Int_t upper = (rowIndex == fNRows - 1) ? nClStack : fPositions[rowIndex + 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 = rowIndex + 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 = rowIndex + 1; j < fNRows; j++) fPositions[j]++;
- nClStack++;
- }
-
- // calculate mean x
- fX += cl->GetX();
-
- // Debug Streaming
- if(AliTRDtrackerV1::DebugStreamer() && AliTRDReconstructor::StreamLevel() >= 3){
- TTreeSRedirector &cstream = *AliTRDtrackerV1::DebugStreamer();
- cstream << "BuildIndices"
- << "Plane=" << fPlane
- << "Stack=" << fStack
- << "Sector=" << fSector
- << "Iter=" << iter
- << "C.=" << cl
- << "rowIndex=" << rowIndex
- << "\n";
- }
- }
-
-// AliInfo("Positions");
-// for(int ir=0; ir<fNRows; ir++) printf("pos[%d] %d\n", ir, fPositions[ir]);
-
- fX /= fN;
-}
-
-//_____________________________________________________________________________
-Int_t AliTRDchamberTimeBin::Find(Float_t y) const
-{
- //
- // Returns index of the cluster nearest in Y
- //
-
- if (fN <= 0) return 0;
-
- if (y <= fClusters[0]->GetY()) return 0;
-
- if (y > fClusters[fN-1]->GetY()) return fN;
-
-
- Int_t b = 0;
- Int_t e = fN - 1;
- Int_t m = (b + e) / 2;
-
- for ( ; b < e; m = (b + e) / 2) {
- if (y > fClusters[m]->GetY()) b = m + 1;
- else e = m;
- }
-
- return m;
-}
-
-//_____________________________________________________________________________
-Int_t AliTRDchamberTimeBin::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 AliTRDchamberTimeBin::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 AliTRDchamberTimeBin::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 = fNRows - 1 - (UChar_t)(TMath::Abs(fZ0 - z)/fZLength * fNRows);
- if(z < fZ0) myZbin = fNRows - 1;
- if(z > fZ0 + fZLength) myZbin = 0;
- //printf("\n%f < %f < %f [%d]\n", fZ0, z, fZ0 + fZLength, myZbin);
- //for(int ic=0; ic<fN; ic++) printf("%d z = %f row %d\n", ic, fClusters[ic]->GetZ(), fClusters[ic]->GetPadRow());
-
- 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){
- //printf("y[%f] ycl[%f] roady[%f]\n", y, ycl, 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 AliTRDchamberTimeBin::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(!AliTRDReconstructor::RecoParam()){
- 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] = AliTRDReconstructor::RecoParam()->GetMaxPhi() * (cl->GetX() - GetX()) + 1.0; // deviation: y-Direction
- cond[3] = AliTRDReconstructor::RecoParam()->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] = AliTRDReconstructor::RecoParam()->GetRoad0y() + phi;
- cond[3] = AliTRDReconstructor::RecoParam()->GetRoad0z();
- }
-}
-
-//_____________________________________________________________________________
-void AliTRDchamberTimeBin::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];
- if(cond[1] - cond[3] > fZ0 + fZLength || cond[1] + cond[3] < fZ0) return; // We are outside of the chamvber
- 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 - 1.E-3;
-
- UChar_t zhi = fNRows - 1 - (UChar_t)(TMath::Abs(fZ0 - zvals[0])/fZLength * fNRows);
- UChar_t zlo = fNRows - 1 - (UChar_t)(TMath::Abs(fZ0 - zvals[1])/fZLength * fNRows);
-
-
-// AliInfo(Form("yc[%f] zc[%f] dy[%f] dz[%f]", cond[0], cond[1], cond[2], cond[3]));
-// PrintClusters();
-// AliInfo(Form("zlo[%f] zhi[%f]", zvals[0], zvals[1]));
-// AliInfo(Form("zlo[%d] zhi[%d]", zlo, zhi));
-
- //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;
- //AliInfo(Form("z[%d] y [%d %d]", z, fPositions[z], upper));
- 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))){/*printf("exit z\n"); TODO*/ continue;}
- index[ncl] = y;
- ncl++;
- }
- }
- if(ncl>fN) AliError(Form("Clusters found %d > %d (clusters in layer)", ncl, fN));
-}
-
-//_____________________________________________________________________________
-AliTRDcluster *AliTRDchamberTimeBin::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 = AliTRDReconstructor::RecoParam()->GetRoad2y();
- Double_t maxroadz = AliTRDReconstructor::RecoParam()->GetRoad2z();
-
- Int_t index = SearchNearestCluster(cond[0],cond[1],maxroad,maxroadz);
- AliTRDcluster *returnCluster = 0x0;
- if(index != -1) returnCluster = (AliTRDcluster *) fClusters[index];
- return returnCluster;
-}
-
-//_____________________________________________________________________________
-void AliTRDchamberTimeBin::PrintClusters() const
-{
-// Prints the position of each cluster in the stacklayer on the stdout
-//
- printf("\nnRows = %d\n", fNRows);
- printf("Z0 = %f\n", fZ0);
- printf("Z1 = %f\n", fZ0+fZLength);
- printf("clusters in AliTRDchamberTimeBin %d\n", fN);
- for(Int_t i = 0; i < fN; i++){
- printf("AliTRDchamberTimeBin: index=%i, Cluster: X = %3.3f [%d] Y = %3.3f [%d] Z = %3.3f [%d]\n", i, fClusters[i]->GetX(), fClusters[i]->GetLocalTimeBin(), fClusters[i]->GetY(), fClusters[i]->GetPadCol(), fClusters[i]->GetZ(), fClusters[i]->GetPadRow());
- if(fClusters[i]->IsUsed()) printf("cluster allready used. rejected in search algorithm\n");
- }
-}
+/**************************************************************************\r
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ * *\r
+ * Author: The ALICE Off-line Project. *\r
+ * Contributors are mentioned in the code where appropriate. *\r
+ * *\r
+ * Permission to use, copy, modify and distribute this software and its *\r
+ * documentation strictly for non-commercial purposes is hereby granted *\r
+ * without fee, provided that the above copyright notice appears in all *\r
+ * copies and that both the copyright notice and this permission notice *\r
+ * appear in the supporting documentation. The authors make no claims *\r
+ * about the suitability of this software for any purpose. It is *\r
+ * provided "as is" without express or implied warranty. *\r
+ **************************************************************************/\r
+\r
+/* $Id: AliTRDchamberTimeBin.cxx 23313 2008-01-11 14:56:43Z cblume $ */\r
+\r
+///////////////////////////////////////////////////////////////////////////////\r
+// //\r
+// Organization of clusters at the level of 1 TRD chamber. //\r
+// The data structure is used for tracking at the stack level. //\r
+// //\r
+// Functionalities: //\r
+// 1. cluster organization and sorting //\r
+// 2. fast data navigation //\r
+// //\r
+// Authors: //\r
+// Alex Bercuci <A.Bercuci@gsi.de> //\r
+// Markus Fasel <M.Fasel@gsi.de> //\r
+// //\r
+///////////////////////////////////////////////////////////////////////////////\r
+\r
+#include <TObject.h>\r
+#include <TROOT.h>\r
+#include <TMath.h>\r
+#include <TStopwatch.h>\r
+#include <TTreeStream.h>\r
+\r
+#include "AliLog.h"\r
+\r
+#include "AliTRDcluster.h"\r
+#include "AliTRDchamberTimeBin.h"\r
+#include "AliTRDrecoParam.h"\r
+#include "AliTRDReconstructor.h"\r
+#include "AliTRDtrackerV1.h"\r
+\r
+\r
+ClassImp(AliTRDchamberTimeBin)\r
+\r
+//_____________________________________________________________________________\r
+AliTRDchamberTimeBin::AliTRDchamberTimeBin(Int_t plane, Int_t stack, Int_t sector, Double_t z0, Double_t zLength)\r
+ :TObject()\r
+ ,fOwner(kFALSE)\r
+ ,fPlane(plane)\r
+ ,fStack(stack)\r
+ ,fSector(sector)\r
+ ,fNRows(kMaxRows)\r
+ ,fN(0)\r
+ ,fX(0.)\r
+ ,fZ0(z0)\r
+ ,fZLength(zLength)\r
+{\r
+ //\r
+ // Default constructor (Only provided to use AliTRDchamberTimeBin with arrays)\r
+ //\r
+\r
+ for(int i=0; i<kMaxRows; i++) fPositions[i] = 0xff;\r
+ for(int ic=0; ic<kMaxClustersLayer; ic++){\r
+ fClusters[ic] = 0x0;\r
+ fIndex[ic] = 0xffff;\r
+ }\r
+}\r
+\r
+// //_____________________________________________________________________________\r
+// AliTRDchamberTimeBin::AliTRDchamberTimeBin(const AliTRDpropagationLayer &layer, Double_t\r
+// z0, Double_t zLength, UChar_t stackNr):\r
+// TObject()\r
+// ,fOwner(kFALSE)\r
+// ,fPlane(0xff)\r
+// ,fStack(0xff)\r
+// ,fSector(0xff)\r
+// ,fNRows(kMaxRows)\r
+// ,fN(0)\r
+// ,fX(0.)\r
+// ,fZ0(z0)\r
+// ,fZLength(zLength)\r
+// {\r
+// // Standard constructor.\r
+// // Initialize also the underlying AliTRDpropagationLayer using the copy constructor.\r
+// \r
+// SetT0(layer.IsT0());\r
+// for(int i=0; i<kMaxRows; i++) fPositions[i] = 0xff;\r
+// for(int ic=0; ic<kMaxClustersLayer; ic++){\r
+// fClusters[ic] = 0x0;\r
+// fIndex[ic] = 0xffff;\r
+// }\r
+// }\r
+// \r
+// //_____________________________________________________________________________\r
+// AliTRDchamberTimeBin::AliTRDchamberTimeBin(const AliTRDpropagationLayer &layer):\r
+// TObject()\r
+// ,fOwner(kFALSE)\r
+// ,fPlane(0xff)\r
+// ,fStack(0xff)\r
+// ,fSector(0xff)\r
+// ,fNRows(kMaxRows)\r
+// ,fN(0)\r
+// ,fX(0.)\r
+// ,fZ0(0)\r
+// ,fZLength(0)\r
+// {\r
+// // Standard constructor using only AliTRDpropagationLayer.\r
+// \r
+// SetT0(layer.IsT0());\r
+// for(int i=0; i<kMaxRows; i++) fPositions[i] = 0xff;\r
+// for(int ic=0; ic<kMaxClustersLayer; ic++){\r
+// fClusters[ic] = 0x0;\r
+// fIndex[ic] = 0xffff;\r
+// }\r
+// }\r
+// //_____________________________________________________________________________\r
+// AliTRDchamberTimeBin &AliTRDchamberTimeBin::operator=(const AliTRDpropagationLayer &layer)\r
+// {\r
+// // Assignment operator from an AliTRDpropagationLayer\r
+// \r
+// if (this != &layer) layer.Copy(*this);\r
+// return *this;\r
+// }\r
+// \r
+\r
+//_____________________________________________________________________________\r
+AliTRDchamberTimeBin::AliTRDchamberTimeBin(const AliTRDchamberTimeBin &layer):\r
+ TObject()\r
+ ,fOwner(layer.fOwner)\r
+ ,fPlane(layer.fPlane)\r
+ ,fStack(layer.fStack)\r
+ ,fSector(layer.fSector)\r
+ ,fNRows(layer.fNRows)\r
+ ,fN(layer.fN)\r
+ ,fX(layer.fX)\r
+ ,fZ0(layer.fZ0)\r
+ ,fZLength(layer.fZLength)\r
+{\r
+// Copy Constructor (performs a deep copy)\r
+ \r
+ SetT0(layer.IsT0());\r
+ for(int i=0; i<kMaxRows; i++) fPositions[i] = layer.fPositions[i];\r
+ memcpy(&fClusters[0], &layer.fClusters[0], kMaxClustersLayer*sizeof(UChar_t));\r
+ memcpy(&fIndex[0], &layer.fIndex[0], kMaxClustersLayer*sizeof(UInt_t));\r
+\r
+\r
+// BuildIndices();\r
+}\r
+\r
+//_____________________________________________________________________________\r
+AliTRDchamberTimeBin &AliTRDchamberTimeBin::operator=(const AliTRDchamberTimeBin &layer)\r
+{\r
+// Assignment operator\r
+\r
+ if (this != &layer) layer.Copy(*this);\r
+ return *this;\r
+}\r
+\r
+//_____________________________________________________________________________\r
+void AliTRDchamberTimeBin::Copy(TObject &o) const\r
+{\r
+// Copy method. Performs a deep copy of all data from this object to object o.\r
+ \r
+ AliTRDchamberTimeBin &layer = (AliTRDchamberTimeBin &)o;\r
+ layer.fOwner = kFALSE;\r
+ layer.fPlane = fPlane;\r
+ layer.fStack = fStack;\r
+ layer.fSector = fSector;\r
+ layer.fNRows = fNRows;\r
+ layer.fN = fN;\r
+ layer.fX = fX;\r
+ layer.fZ0 = fZ0;\r
+ layer.fZLength = fZLength;\r
+ layer.SetT0(IsT0());\r
+ \r
+ for(int i = 0; i < kMaxRows; i++) layer.fPositions[i] = 0;\r
+ \r
+ for(int i=0; i<kMaxRows; i++) layer.fPositions[i] = fPositions[i];\r
+ memcpy(&layer.fClusters[0], &fClusters[0], kMaxClustersLayer*sizeof(UChar_t));\r
+ memcpy(&layer.fIndex[0], &fIndex[0], kMaxClustersLayer*sizeof(UInt_t));\r
+ \r
+ TObject::Copy(layer); // copies everything into layer\r
+ \r
+// layer.BuildIndices();\r
+}\r
+\r
+//_____________________________________________________________________________\r
+AliTRDchamberTimeBin::~AliTRDchamberTimeBin()\r
+{\r
+// Destructor\r
+ if(fOwner) for(int ic=0; ic<fN; ic++) delete fClusters[ic];\r
+}\r
+\r
+//_____________________________________________________________________________\r
+void AliTRDchamberTimeBin::SetRange(const Float_t z0, const Float_t zLength)\r
+{\r
+// Sets the range in z-direction\r
+//\r
+// Parameters:\r
+// z0 : starting position of layer in the z direction\r
+// zLength : length of layer in the z direction \r
+\r
+ fZ0 = (z0 <= z0 + zLength) ? z0 : z0 + zLength;\r
+ fZLength = TMath::Abs(zLength);\r
+}\r
+\r
+//_____________________________________________________________________________\r
+void AliTRDchamberTimeBin::InsertCluster(AliTRDcluster *c, UInt_t index) \r
+{\r
+ //\r
+ // Insert cluster in cluster array.\r
+ // Clusters are sorted according to Y coordinate. \r
+ //\r
+\r
+ //if (fTimeBinIndex < 0) { \r
+ //AliWarning("Attempt to insert cluster into non-sensitive time bin!\n");\r
+ //return;\r
+ //}\r
+\r
+ if (fN == (Int_t) kMaxClustersLayer) {\r
+ //AliWarning("Too many clusters !\n"); \r
+ return;\r
+ }\r
+\r
+ if (fN == 0) {\r
+ fIndex[0] = index; \r
+ fClusters[fN++] = c; \r
+ return;\r
+ }\r
+\r
+ Int_t i = Find(c->GetY());\r
+ memmove(fClusters+i+1,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));\r
+ memmove(fIndex +i+1,fIndex +i,(fN-i)*sizeof(UInt_t)); \r
+ fIndex[i] = index; \r
+ fClusters[i] = c; \r
+ fN++;\r
+\r
+} \r
+\r
+\r
+//_____________________________________________________________________________\r
+void AliTRDchamberTimeBin::BuildIndices(Int_t iter)\r
+{\r
+// Rearrangement of the clusters belonging to the propagation layer for the stack.\r
+//\r
+// Detailed description\r
+//\r
+// The array indices of all clusters in one PropagationLayer are stored in\r
+// array. The array is divided into several bins.\r
+// The clusters are sorted in increasing order of their y coordinate.\r
+//\r
+// Sorting algorithm: TreeSearch\r
+//\r
+\r
+ if(!fN) return;\r
+\r
+ // Select clusters that belong to the Stack\r
+ Int_t nClStack = 0; // Internal counter\r
+ for(Int_t i = 0; i < fN; i++){\r
+ if(fClusters[i]->IsUsed()){\r
+ fClusters[i] = 0x0;\r
+ fIndex[i] = 0xffff;\r
+ } else nClStack++;\r
+ }\r
+ if(nClStack > kMaxClustersLayer) AliWarning(Form("Number of clusters in stack %d exceed buffer size %d. Truncating.", nClStack, kMaxClustersLayer));\r
+ \r
+ // Nothing in this time bin. Reset indexes \r
+ if(!nClStack){\r
+ fN = 0;\r
+ memset(&fPositions[0], 0xff, sizeof(UChar_t) * kMaxRows);\r
+ memset(&fClusters[0], 0x0, sizeof(AliTRDcluster*) * kMaxClustersLayer);\r
+ memset(&fIndex[0], 0xffff, sizeof(UInt_t) * kMaxClustersLayer);\r
+ return;\r
+ }\r
+ \r
+ // Make a copy\r
+ AliTRDcluster *helpCL[kMaxClustersLayer];\r
+ Int_t helpInd[kMaxClustersLayer];\r
+ nClStack = 0;\r
+ for(Int_t i = 0; i < TMath::Min(fN, kMaxClustersLayer); i++){\r
+ if(!fClusters[i]) continue;\r
+ helpCL[nClStack] = fClusters[i];\r
+ helpInd[nClStack] = fIndex[i];\r
+ fClusters[i] = 0x0;\r
+ fIndex[i] = 0xffff;\r
+ nClStack++;\r
+ }\r
+ \r
+ // do clusters arrangement\r
+ fX = 0.;\r
+ fN = nClStack;\r
+ nClStack = 0;\r
+ // Reset Positions array\r
+ memset(fPositions, 0, sizeof(UChar_t)*kMaxRows);\r
+ for(Int_t i = 0; i < fN; i++){\r
+ // boundary check\r
+ AliTRDcluster *cl = helpCL[i];\r
+ UChar_t rowIndex = cl->GetPadRow();\r
+ // Insert Leaf\r
+ Int_t pos = FindYPosition(cl->GetY(), rowIndex, i);\r
+ if(pos == -1){ // zbin is empty;\r
+ Int_t upper = (rowIndex == fNRows - 1) ? nClStack : fPositions[rowIndex + 1];\r
+ memmove(fClusters + upper + 1, fClusters + upper, (sizeof(AliTRDcluster *))*(nClStack-upper));\r
+ memmove(fIndex + upper + 1, fIndex + upper, (sizeof(UInt_t))*(nClStack-upper));\r
+ fClusters[upper] = cl;\r
+ fIndex[upper] = helpInd[i]; \r
+ // Move All pointer one position back\r
+ for(UChar_t j = rowIndex + 1; j < fNRows; j++) fPositions[j]++;\r
+ nClStack++;\r
+ } else { // zbin not empty\r
+ memmove(fClusters + pos + 2, fClusters + pos+1, (sizeof(AliTRDcluster *))*(nClStack-(pos+1)));\r
+ memmove(fIndex + pos + 2, fIndex + pos+1, (sizeof(UInt_t))*(nClStack-(pos+1)));\r
+ fClusters[pos + 1] = cl; //fIndex[i];\r
+ fIndex[pos + 1] = helpInd[i];\r
+ // Move All pointer one position back\r
+ for(UChar_t j = rowIndex + 1; j < fNRows; j++) fPositions[j]++; \r
+ nClStack++;\r
+ }\r
+\r
+ // calculate mean x\r
+ fX += cl->GetX(); \r
+ \r
+ // Debug Streaming\r
+ if(AliTRDtrackerV1::DebugStreamer() && AliTRDReconstructor::RecoParam()->GetStreamLevel() >= 3){\r
+ TTreeSRedirector &cstream = *AliTRDtrackerV1::DebugStreamer();\r
+ cstream << "BuildIndices"\r
+ << "Plane=" << fPlane\r
+ << "Stack=" << fStack\r
+ << "Sector=" << fSector\r
+ << "Iter=" << iter\r
+ << "C.=" << cl\r
+ << "rowIndex=" << rowIndex\r
+ << "\n";\r
+ }\r
+ }\r
+\r
+// AliInfo("Positions");\r
+// for(int ir=0; ir<fNRows; ir++) printf("pos[%d] %d\n", ir, fPositions[ir]);\r
+\r
+ fX /= fN;\r
+}\r
+\r
+//_____________________________________________________________________________\r
+Int_t AliTRDchamberTimeBin::Find(Float_t y) const\r
+{\r
+ //\r
+ // Returns index of the cluster nearest in Y \r
+ //\r
+\r
+ if (fN <= 0) return 0;\r
+ \r
+ if (y <= fClusters[0]->GetY()) return 0;\r
+ \r
+ if (y > fClusters[fN-1]->GetY()) return fN;\r
+ \r
+\r
+ Int_t b = 0;\r
+ Int_t e = fN - 1;\r
+ Int_t m = (b + e) / 2;\r
+\r
+ for ( ; b < e; m = (b + e) / 2) {\r
+ if (y > fClusters[m]->GetY()) b = m + 1;\r
+ else e = m;\r
+ }\r
+\r
+ return m;\r
+} \r
+\r
+//_____________________________________________________________________________\r
+Int_t AliTRDchamberTimeBin::FindYPosition(Double_t y, UChar_t z, Int_t nClusters) const\r
+{\r
+//\r
+// Tree search Algorithm to find the nearest left cluster for a given\r
+// y-position in a certain z-bin (in fact AVL-tree). \r
+// Making use of the fact that clusters are sorted in y-direction.\r
+//\r
+// Parameters:\r
+// y : y position of the reference point in tracking coordinates\r
+// z : z reference bin.\r
+// nClusters : \r
+//\r
+// Output :\r
+// Index of the nearest left cluster in the StackLayer indexing (-1 if no clusters are found)\r
+//\r
+\r
+ Int_t start = fPositions[z]; // starting Position of the bin\r
+ Int_t upper = (Int_t)((z != fNRows - 1) ? fPositions[z+1] : nClusters); // ending Position of the bin \r
+ Int_t end = upper - 1; // ending Position of the bin \r
+ if(end < start) return -1; // Bin is empty\r
+ Int_t middle = static_cast<Int_t>((start + end)/2);\r
+ // 1st Part: climb down the tree: get the next cluster BEFORE ypos\r
+ while(start + 1 < end){\r
+ if(y >= fClusters[middle]->GetY()) start = middle;\r
+ else end = middle;\r
+ middle = static_cast<Int_t>((start + end)/2);\r
+ }\r
+ if(y > fClusters[end]->GetY()) return end;\r
+ return start;\r
+}\r
+\r
+//_____________________________________________________________________________\r
+Int_t AliTRDchamberTimeBin::FindNearestYCluster(Double_t y, UChar_t z) const\r
+{\r
+//\r
+// Tree search Algorithm to find the nearest cluster for a given\r
+// y-position in a certain z-bin (in fact AVL-tree). \r
+// Making use of the fact that clusters are sorted in y-direction.\r
+//\r
+// Parameters:\r
+// y : y position of the reference point in tracking coordinates\r
+// z : z reference bin.\r
+//\r
+// Output \r
+// Index of the nearest cluster in the StackLayer indexing (-1 if no clusters are found)\r
+//\r
+\r
+ Int_t position = FindYPosition(y, z, fN);\r
+ if(position == -1) return position; // bin empty\r
+ // FindYPosition always returns the left Neighbor. We don't know if the left or the right Neighbor is nearest\r
+ // to the Reference y-position, so test both\r
+ Int_t upper = (Int_t)((z < fNRows-1) ? fPositions[z+1] : fN); // ending Position of the bin\r
+ if((position + 1) < (upper)){\r
+ if(TMath::Abs(y - fClusters[position + 1]->GetY()) < TMath::Abs(y - fClusters[position]->GetY())) return position + 1;\r
+ else return position;\r
+ }\r
+ return position;\r
+}\r
+\r
+//_____________________________________________________________________________\r
+Int_t AliTRDchamberTimeBin::SearchNearestCluster(Double_t y, Double_t z, Double_t maxroady, Double_t maxroadz) const\r
+{\r
+//\r
+// Finds the nearest cluster from a given point in a defined range.\r
+// Distance is determined in a 2D space by the 2-Norm.\r
+//\r
+// Parameters :\r
+// y : y position of the reference point in tracking coordinates\r
+// z : z reference bin.\r
+// maxroady : maximum searching distance in y direction\r
+// maxroadz : maximum searching distance in z direction\r
+//\r
+// Output \r
+// Index of the nearest cluster in the StackLayer indexing (-1 if no cluster is found).\r
+// Cluster can be accessed with the operator[] or GetCluster(Int_t index)\r
+//\r
+// Detail description\r
+//\r
+// The following steps are perfomed:\r
+// 1. Get the expected z bins inside maxroadz.\r
+// 2. For each z bin find nearest y cluster.\r
+// 3. Select best candidate\r
+//\r
+ Int_t index = -1;\r
+ // initial minimal distance will be represented as ellipse: semi-major = z-direction\r
+ // later 2-Norm will be used \r
+// Float_t nExcentricity = TMath::Sqrt(maxroadz*maxroadz - maxroad*maxroad)/maxroadz;\r
+ Float_t mindist = maxroadz;\r
+ \r
+ // not very nice but unfortunately neccessarry: we have ho check the neighbors in both directions (+ and -) too. How \r
+ // much neighbors depends on the Quotient maxroadz/fZLength \r
+ UChar_t maxRows = 3;\r
+ UChar_t zpos[kMaxRows];\r
+ // Float_t mindist = TMath::Sqrt(maxroad*maxroad + maxroadz*maxroadz);\r
+// UChar_t myZbin = FindTreePosition(z, fZ0 + fZLength/2, fZLength/4, 8, 8, kFALSE);\r
+ UChar_t myZbin = fNRows - 1 - (UChar_t)(TMath::Abs(fZ0 - z)/fZLength * fNRows);\r
+ if(z < fZ0) myZbin = fNRows - 1;\r
+ if(z > fZ0 + fZLength) myZbin = 0;\r
+ //printf("\n%f < %f < %f [%d]\n", fZ0, z, fZ0 + fZLength, myZbin);\r
+ //for(int ic=0; ic<fN; ic++) printf("%d z = %f row %d\n", ic, fClusters[ic]->GetZ(), fClusters[ic]->GetPadRow());\r
+\r
+ UChar_t nNeighbors = 0;\r
+ for(UChar_t i = 0; i < maxRows; i++){\r
+ if((myZbin - 1 + i) < 0) continue;\r
+ if((myZbin - 1 + i) > fNRows - 1) break;\r
+ zpos[nNeighbors] = myZbin - 1 + i;\r
+ nNeighbors++;\r
+ }\r
+ Float_t ycl = 0, zcl = 0;\r
+ for(UChar_t neighbor = 0; neighbor < nNeighbors; neighbor++){ // Always test the neighbors too\r
+ Int_t pos = FindNearestYCluster(y, zpos[neighbor]);\r
+ if(pos == -1) continue; // No cluster in bin\r
+ AliTRDcluster *c = (AliTRDcluster *) (fClusters[pos]);\r
+ if(c->IsUsed()) continue; // we are only interested in unused clusters\r
+ ycl = c->GetY();\r
+ // Too far away in y-direction (Prearrangement)\r
+ if (TMath::Abs(ycl - y) > maxroady){ \r
+ //printf("y[%f] ycl[%f] roady[%f]\n", y, ycl, maxroady);\r
+ continue;\r
+ }\r
+ zcl = c->GetZ();\r
+ // Too far away in z-Direction\r
+ // (Prearrangement since we have not so many bins to test)\r
+ if (TMath::Abs(zcl - z) > maxroadz) continue;\r
+ \r
+ Float_t dist; // distance defined as 2-Norm \r
+ // if we havent found a Particle that is in the ellipse around (y,z) with maxroad as semi-minor and\r
+ // maxroadz as semi-major, we take the radius of the ellipse concerning the cluster as mindist, later we \r
+ // take the 2-Norm when we found a cluster inside the ellipse (The value 10000 is taken because it is surely\r
+ // large enough to be usable as an indicator whether we have found a nearer cluster or not)\r
+// if(mindist > 10000.){\r
+// Float_t phi = ((zcl - z) == 0) ? TMath::Pi()/2 : TMath::ATan((ycl - y)/(zcl - z));\r
+// mindist = maxroad/TMath::Sqrt(1 - nExcentricity*nExcentricity * (TMath::Cos(phi))*(TMath::Cos(phi)));\r
+// }\r
+ dist = TMath::Max(TMath::Abs(y-ycl),TMath::Abs(z-zcl)); // infinity Norm\r
+// dist = TMath::Sqrt((ycl - y)*(ycl - y) + (zcl - z)*(zcl - z));\r
+ if((Int_t)(dist * 100000) < (Int_t)(mindist * 100000)){\r
+ //if((dist = TMath::Sqrt((ycl - y)*(ycl - y) + (zcl - z)*(zcl - z))) < mindist){\r
+ mindist = dist;\r
+ index = pos;\r
+ } \r
+ }\r
+ // This is the Array Position in fIndex2D of the Nearest cluster: if a\r
+ // cluster is called, then the function has to retrieve the Information\r
+ // which is Stored in the Array called, the function\r
+ return index;\r
+}\r
+\r
+//_____________________________________________________________________________\r
+void AliTRDchamberTimeBin::BuildCond(AliTRDcluster *cl, Double_t *cond, UChar_t Layer, Double_t theta, Double_t phi)\r
+{\r
+// Helper function to calculate the area where to expect a cluster in THIS\r
+// layer. \r
+//\r
+// Parameters :\r
+// cl : \r
+// cond :\r
+// Layer : \r
+// theta : \r
+// phi : \r
+//\r
+// Detail description\r
+//\r
+// Helper function to calculate the area where to expect a cluster in THIS\r
+// layer. by using the information of a former cluster in another layer\r
+// and the angle in theta- and phi-direction between layer 0 and layer 3.\r
+// If the layer is zero, initial conditions are calculated. Otherwise a\r
+// linear interpolation is performed. \r
+//Begin_Html\r
+//<img src="gif/build_cond.gif">\r
+//End_Html\r
+//\r
+\r
+ if(!AliTRDReconstructor::RecoParam()){\r
+ AliError("Reconstruction parameters not initialized.");\r
+ return;\r
+ }\r
+ \r
+ if(Layer == 0){\r
+ cond[0] = cl->GetY(); // center: y-Direction\r
+ cond[1] = cl->GetZ(); // center: z-Direction\r
+ cond[2] = AliTRDReconstructor::RecoParam()->GetMaxPhi() * (cl->GetX() - GetX()) + 1.0; // deviation: y-Direction\r
+ cond[3] = AliTRDReconstructor::RecoParam()->GetMaxTheta() * (cl->GetX() - GetX()) + 1.0; // deviation: z-Direction\r
+ } else {\r
+ cond[0] = cl->GetY() + phi * (GetX() - cl->GetX()); \r
+ cond[1] = cl->GetZ() + theta * (GetX() - cl->GetX());\r
+ cond[2] = AliTRDReconstructor::RecoParam()->GetRoad0y() + phi;\r
+ cond[3] = AliTRDReconstructor::RecoParam()->GetRoad0z();\r
+ }\r
+}\r
+\r
+//_____________________________________________________________________________\r
+void AliTRDchamberTimeBin::GetClusters(Double_t *cond, Int_t *index, Int_t& ncl, Int_t BufferSize)\r
+{\r
+// Finds all clusters situated in this layer inside a rectangle given by the center an ranges.\r
+//\r
+// Parameters :\r
+// cond :\r
+// index : \r
+// ncl :\r
+// BufferSize :\r
+//\r
+// Output :\r
+//\r
+// Detail description\r
+//\r
+// Function returs an array containing the indices in the stacklayer of\r
+// the clusters found an the number of found clusters in the stacklayer\r
+\r
+ ncl = 0;\r
+ memset(index, 0, BufferSize*sizeof(Int_t));\r
+ if(fN == 0) return;\r
+ \r
+ //Boundary checks\r
+ Double_t zvals[2];\r
+ if(cond[1] - cond[3] > fZ0 + fZLength || cond[1] + cond[3] < fZ0) return; // We are outside of the chamvber\r
+ zvals[0] = ((cond[1] - cond[3]) < fZ0) ? fZ0 : (cond[1] - cond[3]);\r
+ zvals[1] = ((cond[1] + cond[3]) < fZ0 + fZLength) ? (cond[1] + cond[3]) : fZ0 + fZLength - 1.E-3;\r
+\r
+ UChar_t zhi = fNRows - 1 - (UChar_t)(TMath::Abs(fZ0 - zvals[0])/fZLength * fNRows);\r
+ UChar_t zlo = fNRows - 1 - (UChar_t)(TMath::Abs(fZ0 - zvals[1])/fZLength * fNRows);\r
+ \r
+ \r
+// AliInfo(Form("yc[%f] zc[%f] dy[%f] dz[%f]", cond[0], cond[1], cond[2], cond[3]));\r
+// PrintClusters();\r
+// AliInfo(Form("zlo[%f] zhi[%f]", zvals[0], zvals[1]));\r
+// AliInfo(Form("zlo[%d] zhi[%d]", zlo, zhi));\r
+ \r
+ //Preordering in Direction z saves a lot of loops (boundary checked)\r
+ for(UChar_t z = zlo; z <= zhi; z++){\r
+ UInt_t upper = (z < fNRows-1) ? fPositions[z+1] : fN;\r
+ //AliInfo(Form("z[%d] y [%d %d]", z, fPositions[z], upper));\r
+ for(Int_t y = fPositions[z]; y < (Int_t)upper; y++){\r
+ if(ncl == BufferSize){\r
+ AliWarning("Buffer size riched. Some clusters may be lost.");\r
+ return; //Buffer filled\r
+ }\r
+ \r
+ if(fClusters[y]->GetY() > (cond[0] + cond[2])) break; // Abbortion conditions!!!\r
+ if(fClusters[y]->GetY() < (cond[0] - cond[2])) continue; // Too small\r
+ if(((Int_t)((fClusters[y]->GetZ())*1000) < (Int_t)(zvals[0]*1000)) || ((Int_t)((fClusters[y]->GetZ())*1000) > (Int_t)(zvals[1]*1000))){/*printf("exit z\n"); TODO*/ continue;}\r
+ index[ncl] = y;\r
+ ncl++;\r
+ }\r
+ }\r
+ if(ncl>fN) AliError(Form("Clusters found %d > %d (clusters in layer)", ncl, fN));\r
+}\r
+\r
+//_____________________________________________________________________________\r
+AliTRDcluster *AliTRDchamberTimeBin::GetNearestCluster(Double_t *cond)\r
+{\r
+// Function returning a pointer to the nearest cluster (nullpointer if not successfull).\r
+//\r
+// Parameters :\r
+// cond :\r
+//\r
+// Output :\r
+// pointer to the nearest cluster (nullpointer if not successfull).\r
+// \r
+// Detail description\r
+//\r
+// returns a pointer to the nearest cluster (nullpointer if not\r
+// successfull) by the help of the method FindNearestCluster\r
+ \r
+ \r
+ Double_t maxroad = AliTRDReconstructor::RecoParam()->GetRoad2y();\r
+ Double_t maxroadz = AliTRDReconstructor::RecoParam()->GetRoad2z();\r
+ \r
+ Int_t index = SearchNearestCluster(cond[0],cond[1],maxroad,maxroadz);\r
+ AliTRDcluster *returnCluster = 0x0;\r
+ if(index != -1) returnCluster = (AliTRDcluster *) fClusters[index];\r
+ return returnCluster;\r
+}\r
+\r
+//_____________________________________________________________________________\r
+void AliTRDchamberTimeBin::PrintClusters() const\r
+{\r
+// Prints the position of each cluster in the stacklayer on the stdout\r
+//\r
+ printf("\nnRows = %d\n", fNRows);\r
+ printf("Z0 = %f\n", fZ0);\r
+ printf("Z1 = %f\n", fZ0+fZLength);\r
+ printf("clusters in AliTRDchamberTimeBin %d\n", fN);\r
+ for(Int_t i = 0; i < fN; i++){\r
+ printf("AliTRDchamberTimeBin: index=%i, Cluster: X = %3.3f [%d] Y = %3.3f [%d] Z = %3.3f [%d]\n", i, fClusters[i]->GetX(), fClusters[i]->GetLocalTimeBin(), fClusters[i]->GetY(), fClusters[i]->GetPadCol(), fClusters[i]->GetZ(), fClusters[i]->GetPadRow());\r
+ if(fClusters[i]->IsUsed()) printf("cluster allready used. rejected in search algorithm\n");\r
+ }\r
+}\r
,fkChi2Z(30./*14.*//*12.5*/)
,fkChi2Y(.25)
,fkTrackLikelihood(-15.)
+ ,fkStreamLevel(0)
+ ,fSeedingOn(kFALSE)
+ ,fVertexConstrained(kTRUE)
,fClusMaxThresh(4.5)
,fClusSigThresh(3.5)
,fLUTOn(kTRUE)
AliTRDrawStreamBase::SetRawStreamVersion("TB");
AliTRDrecoParam *par = new AliTRDrecoParam();
par->SetADCbaseline(10);
+ par->SetSeedingOn(kTRUE);
+ par->SetVertexConstrained(kTRUE);
return par;
}
Double_t GetPlaneQualityThreshold() const { return fkPlaneQualityThreshold; }
Double_t GetTrackLikelihood() const { return fkTrackLikelihood; }
+ Int_t GetStreamLevel() const { return fkStreamLevel; }
+ Bool_t SeedingOn() const { return fSeedingOn; }
+ Bool_t IsVertexConstrained() const { return fVertexConstrained; }
+
Double_t GetClusMaxThresh() const { return fClusMaxThresh; };
Double_t GetClusSigThresh() const { return fClusSigThresh; };
Int_t GetTCnexp() const { return fTCnexp; };
Bool_t TCOn() const { return fTCOn; };
Int_t GetADCbaseline() const { return fADCbaseline; };
-
+
static AliTRDrecoParam *GetLowFluxParam();
static AliTRDrecoParam *GetHighFluxParam();
static AliTRDrecoParam *GetCosmicTestParam();
void SetClusterSharing(Bool_t share = kTRUE) { fkClusterSharing = share; };
void SetPIDMethod(Int_t pid = 1) { fkPIDMethod = pid ? 1 : 0; };
+ void SetSeedingOn(Bool_t seedingOn = kTRUE) { fSeedingOn = seedingOn; }
+ void SetVertexConstrained(Bool_t vertexConstrained = kTRUE) { fVertexConstrained = vertexConstrained; }
+ void SetStreamLevel(Int_t streamLevel= 1) { fkStreamLevel = streamLevel; }
void SetLUT(Int_t lutOn = 1) { fLUTOn = lutOn; };
void SetClusMaxThresh(Float_t thresh) { fClusMaxThresh = thresh; };
void SetClusSigThresh(Float_t thresh) { fClusSigThresh = thresh; };
Double_t fkChi2Z; // Max chi2 on the z direction for seeding clusters fit
Double_t fkChi2Y; // Max chi2 on the y direction for seeding clusters Rieman fit
Double_t fkTrackLikelihood; // Track likelihood for tracklets Rieman fit
+ Int_t fkStreamLevel; // Streaming Level in TRD Reconstruction
+
+ Bool_t fSeedingOn; // Do stand alone tracking in the TRD
+ Bool_t fVertexConstrained; // Perform vertex constrained fit
// Clusterization parameter
Double_t fClusMaxThresh; // Threshold value for cluster maximum
}
AliTRDchamberTimeBin *layer = 0x0;
- if(AliTRDReconstructor::StreamLevel()>=7 && c){
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel()>=7 && c){
TClonesArray clusters("AliTRDcluster", 24);
clusters.SetOwner(kTRUE);
AliTRDcluster *cc = 0x0;
fZ[iTime] = cl->GetZ();
ncl++;
}
- if(AliTRDReconstructor::StreamLevel()>=7) AliInfo(Form("iter = %d ncl [%d] = %d", iter, fPlane, ncl));
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel()>=7) AliInfo(Form("iter = %d ncl [%d] = %d", iter, fPlane, ncl));
if(ncl>1){
// calculate length of the time bin (calibration aware)
AliTRDseed::Update();
}
- if(AliTRDReconstructor::StreamLevel()>=7) AliInfo(Form("iter = %d nclFit [%d] = %d", iter, fPlane, fN2));
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel()>=7) AliInfo(Form("iter = %d nclFit [%d] = %d", iter, fPlane, fN2));
if(IsOK()){
tquality = GetQuality(kZcorr);
#include "AliTRDCommonParam.h"
#include "AliTRDtracker.h"
#include "AliTRDReconstructor.h"
+#include "AliTRDrecoParam.h"
#include "AliTRDCalibraFillHisto.h"
ClassImp(AliTRDtracker)
Int_t nSeed = event->GetNumberOfTracks();
if(!nSeed){
// run stand alone tracking
- if (AliTRDReconstructor::SeedingOn()) Clusters2Tracks(event);
+ if (AliTRDReconstructor::RecoParam()->SeedingOn()) Clusters2Tracks(event);
return 0;
}
// Debug part of tracking
TTreeSRedirector &cstream = *fDebugStreamer;
Int_t eventNrInFile = event->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
- if (AliTRDReconstructor::StreamLevel() > 0) {
+ if (AliTRDReconstructor::RecoParam()->GetStreamLevel() > 0) {
if (track->GetBackupTrack()) {
cstream << "Tracks"
<< "EventNrInFile=" << eventNrInFile
// Add TRD track to ESDfriendTrack - maybe this tracks are
// not useful for post-processing - TODO make decision
- if (AliTRDReconstructor::StreamLevel() > 0) {
+ if (AliTRDReconstructor::RecoParam()->GetStreamLevel() > 0) {
seed->AddCalibObject(new AliTRDtrack(*pt2/*, kTRUE*/));
}
delete pt2;
}
// Add TRD track to ESDfriendTrack
- if (AliTRDReconstructor::StreamLevel() > 0) {
+ if (AliTRDReconstructor::RecoParam()->GetStreamLevel() > 0) {
seed->AddCalibObject(new AliTRDtrack(*pt/*, kTRUE*/));
}
delete pt;
isFake = kTRUE;
}
- if (AliTRDReconstructor::StreamLevel() > 0) {
+ if (AliTRDReconstructor::RecoParam()->GetStreamLevel() > 0) {
if ((!isFake) || ((icl3%10) == 0)) { // Debugging print
TTreeSRedirector &cstream = *fDebugStreamer;
cstream << "Seeds0"
if (1 || (!isFake)) {
Float_t zvertex = GetZ();
TTreeSRedirector &cstream = *fDebugStreamer;
- if (AliTRDReconstructor::StreamLevel() > 0) {
+ if (AliTRDReconstructor::RecoParam()->GetStreamLevel() > 0) {
cstream << "Seeds1"
<< "isFake=" << isFake
<< "Vertex=" << zvertex
esdtrack.UpdateTrackParams(track,AliESDtrack::kTRDout);
esdtrack.SetLabel(label);
esd->AddTrack(&esdtrack);
- if (AliTRDReconstructor::StreamLevel() > 0) {
+ if (AliTRDReconstructor::RecoParam()->GetStreamLevel() > 0) {
cstream << "Tracks"
<< "EventNrInFile=" << eventNrInFile
<< "ESD.=" << &esdtrack
}
}
- if (AliTRDReconstructor::StreamLevel() > 0) {
+ if (AliTRDReconstructor::RecoParam()->GetStreamLevel() > 0) {
cstream << "Seeds2"
<< "Iter=" << jter
<< "Track.=" << track
TGraph graphy(t1-t0,x,yt);
TGraph graphz(t1-t0,x,zt);
- if (AliTRDReconstructor::StreamLevel() > 0) {
+ if (AliTRDReconstructor::RecoParam()->GetStreamLevel() > 0) {
cstream << "tracklet"
<< "track.=" << track // Track parameters
<< "tany=" << tany // Tangent of the local track angle
if (!AliTRDcalibDB::Instance()) {
AliFatal("Could not get calibration object");
}
- fgNTimeBins = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
+ if(!fgNTimeBins) fgNTimeBins = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
for (Int_t isector = 0; isector < AliTRDgeometry::kNsect; isector++) new(&fTrSec[isector]) AliTRDtrackingSector(fGeom, isector);
+
+ if (!AliTRDReconstructor::RecoParam()){
+ AliWarning("RecoParams not set in AliTRDReconstructor. Setting to default LowFluxParam.");
+ AliTRDReconstructor::SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
+ }
- if(AliTRDReconstructor::StreamLevel() > 1){
+
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel() > 1){
TDirectory *savedir = gDirectory;
fgDebugStreamer = new TTreeSRedirector("TRD.TrackerDebug.root");
savedir->cd();
Int_t nSeed = event->GetNumberOfTracks();
if(!nSeed){
// run stand alone tracking
- if (AliTRDReconstructor::SeedingOn()) Clusters2Tracks(event);
+ if (AliTRDReconstructor::RecoParam()->SeedingOn()) Clusters2Tracks(event);
return 0;
}
track.UpdateESDtrack(seed);
// Add TRD track to ESDfriendTrack
- if (AliTRDReconstructor::StreamLevel() > 0 /*&& quality TODO*/){
+ if (AliTRDReconstructor::RecoParam()->GetStreamLevel() > 0 /*&& quality TODO*/){
AliTRDtrackV1 *calibTrack = new AliTRDtrackV1(track);
calibTrack->SetOwner();
seed->AddCalibObject(calibTrack);
track.UpdateESDtrack(seed);
// Add TRD track to ESDfriendTrack
-// if (AliTRDReconstructor::StreamLevel() > 0 /*&& quality TODO*/){
+// if (AliTRDReconstructor::RecoParam()->GetStreamLevel() > 0 /*&& quality TODO*/){
// AliTRDtrackV1 *calibTrack = new AliTRDtrackV1(track);
// calibTrack->SetOwner();
// seed->AddCalibObject(calibTrack);
track.UpdateESDtrack(seed);
// Add TRD track to ESDfriendTrack
-// if (AliTRDReconstructor::StreamLevel() > 0 /*&& quality TODO*/){
+// if (AliTRDReconstructor::RecoParam()->GetStreamLevel() > 0 /*&& quality TODO*/){
// AliTRDtrackV1 *calibTrack = new AliTRDtrackV1(track);
// calibTrack->SetOwner();
// seed->AddCalibObject(calibTrack);
}
}
- if(AliTRDReconstructor::StreamLevel() > 1){
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel() > 1){
Int_t index;
for(int iplane=0; iplane<6; iplane++){
AliTRDseedV1 *tracklet = GetTracklet(&t, iplane, index);
} // end planes loop
- if(AliTRDReconstructor::StreamLevel() > 1){
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel() > 1){
TTreeSRedirector &cstreamer = *fgDebugStreamer;
Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
//AliTRDtrackV1 *debugTrack = new AliTRDtrackV1(t);
for(Int_t ip = 0; ip < AliTRDtrackerV1::kNPlanes; ip++)
tracklets[ip].SetCC(curvature);
- if(AliTRDReconstructor::StreamLevel() >= 5){
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel() >= 5){
//Linear Model on z-direction
Double_t xref = CalculateReferenceX(tracklets); // Relative to the middle of the stack
Double_t slope = fitter->GetParameter(2);
tracklets[iLayer].SetChi2(chi2track);
}
- if(AliTRDReconstructor::StreamLevel() >=5){
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel() >=5){
TTreeSRedirector &cstreamer = *fgDebugStreamer;
Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
}
}
- if(AliTRDReconstructor::StreamLevel() >=5){
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel() >=5){
TTreeSRedirector &cstreamer = *fgDebugStreamer;
Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
return new ((*fTracklets)[nentries]) AliTRDseedV1(*tracklet);
}
+//____________________________________________________________________
+AliTRDtrackV1* AliTRDtrackerV1::SetTrack(AliTRDtrackV1 *track)
+{
+ // Add this track to the list of tracks stored in the tracker
+ //
+ // Parameters
+ // - track : pointer to the track to be added to the list
+ //
+ // Output
+ // - the pointer added
+ //
+ // Detailed description
+ // Build the tracks list if it is not yet created (late initialization)
+ // and adds the new track to the list.
+ //
+ if(!fTracks){
+ fTracks = new TClonesArray("AliTRDtrackV1", AliTRDgeometry::Nsect()*kMaxTracksStack);
+ fTracks->SetOwner(kTRUE);
+ }
+ Int_t nentries = fTracks->GetEntriesFast();
+ return new ((*fTracks)[nentries]) AliTRDtrackV1(*track);
+}
+
+
//____________________________________________________________________
Int_t AliTRDtrackerV1::Clusters2TracksSM(Int_t sector, AliESDEvent *esd)
// Build initial seeding configurations
Double_t quality = BuildSeedingConfigs(stack, configs);
- if(AliTRDReconstructor::StreamLevel() > 1){
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel() > 1){
AliInfo(Form("Plane config %d %d %d Quality %f"
, configs[0], configs[1], configs[2], quality));
}
ntracks = MakeSeeds(stack, &sseed[6*ntracks], pars);
if(ntracks == kMaxTracksStack) break;
}
- if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Candidate TRD tracks %d in iteration %d.", ntracks, fSieveSeeding));
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel() > 1) AliInfo(Form("Candidate TRD tracks %d in iteration %d.", ntracks, fSieveSeeding));
if(!ntracks) break;
Int_t ich = 0; while(!(chamber = stack[ich])) ich++;
trackParams[6] = fGeom->GetSector(chamber->GetDetector());/* *alpha+shift; // Supermodule*/
- if(AliTRDReconstructor::StreamLevel() > 1){
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel() > 1){
AliInfo(Form("Track %d [%d] nlayers %d trackQuality = %e nused %d, yref = %3.3f", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused, trackParams[1]));
Int_t nclusters = 0;
esdTrack.SetLabel(track->GetLabel());
track->UpdateESDtrack(&esdTrack);
// write ESD-friends if neccessary
- if (AliTRDReconstructor::StreamLevel() > 0){
+ if (AliTRDReconstructor::RecoParam()->GetStreamLevel() > 0){
//printf("Creating Calibrations Object\n");
AliTRDtrackV1 *calibTrack = new AliTRDtrackV1(*track);
calibTrack->SetOwner();
chamber->Build(fGeom);//Indices(fSieveSeeding);
}
- if(AliTRDReconstructor::StreamLevel() > 1){
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel() > 1){
AliInfo(Form("Sieve level %d Plane config %d %d %d Quality %f", fSieveSeeding, configs[0], configs[1], configs[2], quality));
}
} while(fSieveSeeding<10); // end stack clusters sieve
padlength[iplane] = pp->GetLengthIPad();
}
- if(AliTRDReconstructor::StreamLevel() > 1){
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel() > 1){
AliInfo(Form("Making seeds Stack[%d] Config[%d] Tracks[%d]...", istack, config, ntracks));
}
}
Bool_t isFake = kFALSE;
- if(AliTRDReconstructor::StreamLevel() >= 2){
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel() >= 2){
if (c[0]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
if (c[1]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
if (c[2]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
// AliInfo("Extrapolation done.");
// Debug Stream containing all the 6 tracklets
- if(AliTRDReconstructor::StreamLevel() >= 2){
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel() >= 2){
TTreeSRedirector &cstreamer = *fgDebugStreamer;
TLinearFitter *tiltedRieman = GetTiltedRiemanFitter();
Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
// do the final track fitting (Once with vertex constraint and once without vertex constraint)
Double_t chi2Vals[3];
chi2Vals[0] = FitTiltedRieman(&cseed[0], kFALSE);
- chi2Vals[1] = FitTiltedRiemanConstraint(&cseed[0], GetZ());
+ if(AliTRDReconstructor::RecoParam()->IsVertexConstrained())
+ chi2Vals[1] = FitTiltedRiemanConstraint(&cseed[0], GetZ()); // Do Vertex Constrained fit if desired
+ else
+ chi2Vals[1] = 1.;
chi2Vals[2] = GetChi2Z(&cseed[0]) / TMath::Max((mlayers - 3.), 1.);
// Chi2 definitions in testing stage
//chi2Vals[2] = GetChi2ZTest(&cseed[0]);
cseed[iLayer].SetChi2Z(chi2[1]);
}
- if(AliTRDReconstructor::StreamLevel() >= 2){
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel() >= 2){
TTreeSRedirector &cstreamer = *fgDebugStreamer;
Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
c[ 6] = 0.0; c[ 7] = 0.0; c[ 8] = 0.0; c[ 9] = 0.1;
c[10] = 0.0; c[11] = 0.0; c[12] = 0.0; c[13] = 0.0; c[14] = params[5]*params[5]*0.01;
- AliTRDtrackV1 *track = new AliTRDtrackV1(seeds, ¶ms[1], c, params[0], params[6]*alpha+shift);
- track->PropagateTo(params[0]-5.0);
- track->ResetCovariance(1);
- Int_t nc = FollowBackProlongation(*track);
- if (nc < 30) {
- delete track;
- track = 0x0;
- } else {
- //track->CookdEdx();
- //track->CookdEdxTimBin(-1);
- track->CookLabel(.9);
- // computes PID for track
- track->CookPID();
- // update calibration references using this track
- if(calibra->GetHisto2d()) calibra->UpdateHistogramsV1(track);
- }
-
- return track;
+ AliTRDtrackV1 track(seeds, ¶ms[1], c, params[0], params[6]*alpha+shift);
+ track.PropagateTo(params[0]-5.0);
+ track.ResetCovariance(1);
+ Int_t nc = FollowBackProlongation(track);
+ if (nc < 30) return 0x0;
+
+ AliTRDtrackV1 *ptrTrack = SetTrack(&track);
+ ptrTrack->CookLabel(.9);
+ // computes PID for track
+ ptrTrack->CookPID();
+ // update calibration references using this track
+ if(calibra->GetHisto2d()) calibra->UpdateHistogramsV1(ptrTrack);
+
+ return ptrTrack;
}
}
chi2 = FitTiltedRieman(bseed, kTRUE);
- if(AliTRDReconstructor::StreamLevel() >= 7){
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel() >= 7){
Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
TLinearFitter *tiltedRieman = GetTiltedRiemanFitter();
sumdaf /= Float_t (nLayers - 2.0);
Double_t likeChi2Z = TMath::Exp(-chi2[2] * 0.14); // Chi2Z
- Double_t likeChi2TC = TMath::Exp(-chi2[1] * 0.677); // Constrained Tilted Riemann
+ Double_t likeChi2TC = (AliTRDReconstructor::RecoParam()->IsVertexConstrained()) ?
+ TMath::Exp(-chi2[1] * 0.677) : 1; // Constrained Tilted Riemann
Double_t likeChi2TR = TMath::Exp(-chi2[0] * 0.78); // Non-constrained Tilted Riemann
Double_t likeAF = TMath::Exp(-sumdaf * 3.23);
Double_t trackLikelihood = likeChi2Z * likeChi2TR * likeAF;
- if(AliTRDReconstructor::StreamLevel() >= 2){
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel() >= 2){
Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
TTreeSRedirector &cstreamer = *fgDebugStreamer;
Double_t like = likea * likechi2y * likechi2z * likeN;
// AliInfo(Form("sumda(%f) chi2[0](%f) chi2[1](%f) likea(%f) likechi2y(%f) likechi2z(%f) nclusters(%d) likeN(%f)", sumda, chi2[0], chi2[1], likea, likechi2y, likechi2z, nclusters, likeN));
- if(AliTRDReconstructor::StreamLevel() >= 2){
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel() >= 2){
Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
// The Debug Stream contains the seed
Int_t MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *sseed, Int_t *ipar);
AliTRDtrackV1* MakeTrack(AliTRDseedV1 *seeds, Double_t *params);
static Int_t PropagateToX(AliTRDtrackV1 &t, Double_t xToGo, Double_t maxStep);
+ AliTRDtrackV1* SetTrack(AliTRDtrackV1 *track);
AliTRDseedV1* SetTracklet(AliTRDseedV1 *tracklet);
private:
-/**************************************************************************
- * 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;
-}
-
+/**************************************************************************\r
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ * *\r
+ * Author: The ALICE Off-line Project. *\r
+ * Contributors are mentioned in the code where appropriate. *\r
+ * *\r
+ * Permission to use, copy, modify and distribute this software and its *\r
+ * documentation strictly for non-commercial purposes is hereby granted *\r
+ * without fee, provided that the above copyright notice appears in all *\r
+ * copies and that both the copyright notice and this permission notice *\r
+ * appear in the supporting documentation. The authors make no claims *\r
+ * about the suitability of this software for any purpose. It is *\r
+ * provided "as is" without express or implied warranty. *\r
+ **************************************************************************/\r
+\r
+/* $Id: AliTRDtrackingDebug.cxx 23810 2008-02-08 09:00:27Z hristov $ */\r
+\r
+////////////////////////////////////////////////////////////////////////////\r
+// //\r
+// Tracking in one chamber //\r
+// //\r
+// Authors: //\r
+// Alex Bercuci <A.Bercuci@gsi.de> //\r
+// Markus Fasel <M.Fasel@gsi.de> //\r
+// //\r
+////////////////////////////////////////////////////////////////////////////\r
+\r
+#include "AliTRDtrackingChamber.h"\r
+\r
+#include "TMath.h"\r
+#include "TMatrixTBase.h"\r
+#include <TTreeStream.h>\r
+\r
+#include "AliTRDReconstructor.h"\r
+#include "AliTRDrecoParam.h"\r
+#include "AliTRDtrackerV1.h"\r
+#include "AliTRDgeometry.h"\r
+#include "AliTRDpadPlane.h"\r
+#include "AliTRDcalibDB.h"\r
+\r
+ClassImp(AliTRDtrackingChamber)\r
+\r
+//_______________________________________________________\r
+AliTRDtrackingChamber::AliTRDtrackingChamber(Int_t det) :\r
+ fDetector(det)\r
+ ,fX0(0.)\r
+{} \r
+\r
+//_______________________________________________________\r
+void AliTRDtrackingChamber::Clear(const Option_t *opt)\r
+{\r
+ for(Int_t itb=0; itb<kNTimeBins; itb++) fTB[itb].Clear(opt);\r
+}\r
+\r
+//_______________________________________________________\r
+void AliTRDtrackingChamber::InsertCluster(AliTRDcluster *c, Int_t index)\r
+{\r
+ fTB[c->GetLocalTimeBin()].InsertCluster(c, index);\r
+}\r
+\r
+//_______________________________________________________\r
+Bool_t AliTRDtrackingChamber::Build(AliTRDgeometry *geo)\r
+{\r
+// Init chamber and all time bins (AliTRDchamberTimeBin)\r
+// Calculates radial position of the chamber based on \r
+// radial positions of the time bins (calibration/alignment aware)\r
+//\r
+ Int_t stack = geo->GetChamber(fDetector);\r
+ Int_t plane = geo->GetPlane(fDetector);\r
+ AliTRDpadPlane *pp = geo->GetPadPlane(plane, stack);\r
+ Double_t zl = pp->GetRow0ROC() - pp->GetRowEndROC();\r
+ Double_t z0 = geo->GetRow0(plane, stack, 0) - zl;\r
+ Int_t nrows = pp->GetNrows();\r
+ \r
+ Int_t index[50], jtb = 0;\r
+ for(Int_t itb=0; itb<kNTimeBins; itb++){ \r
+ if(!fTB[itb]) continue;\r
+ fTB[itb].SetRange(z0, zl);\r
+ fTB[itb].SetNRows(nrows);\r
+ fTB[itb].BuildIndices();\r
+ index[jtb++] = itb;\r
+ } \r
+ if(jtb<2) return kFALSE;\r
+ \r
+ \r
+ // ESTIMATE POSITION OF PAD PLANE FOR THIS CHAMBER\r
+ Double_t x0 = fTB[index[0]].GetX();\r
+ Double_t x1 = fTB[index[1]].GetX();\r
+ Double_t dx = (x0 - x1)/(index[1] - index[0]); \r
+ fX0 = x0 + dx*(index[0] - (Int_t)AliTRDcalibDB::Instance()->GetT0Average(fDetector)); \r
+ return kTRUE;\r
+}\r
+ \r
+//_______________________________________________________ \r
+Int_t AliTRDtrackingChamber::GetNClusters() const\r
+{\r
+// Returns number of clusters in chamber\r
+//\r
+ Int_t n = 0;\r
+ for(Int_t itb=0; itb<kNTimeBins; itb++){ \r
+ n += Int_t(fTB[itb]);\r
+ }\r
+ return n; \r
+} \r
+\r
+//_______________________________________________________\r
+Double_t AliTRDtrackingChamber::GetQuality()\r
+{\r
+ //\r
+ // Calculate chamber quality for seeding.\r
+ // \r
+ //\r
+ // Parameters :\r
+ // layers : Array of propagation layers for this plane.\r
+ //\r
+ // Output :\r
+ // plane quality factor for seeding\r
+ // \r
+ // Detailed description\r
+ //\r
+ // The quality of the plane for seeding is higher if:\r
+ // 1. the average timebin population is closer to an integer number\r
+ // 2. the distribution of clusters/timebin is closer to a uniform distribution.\r
+ // - the slope of the first derivative of a parabolic fit is small or\r
+ // - the slope of a linear fit is small\r
+ //\r
+\r
+ Int_t ncl = 0;\r
+ Int_t nused = 0;\r
+ Int_t nClLayer;\r
+ for(int itb=0; itb<kNTimeBins; itb++){\r
+ if(!(nClLayer = fTB[itb].GetNClusters())) continue;\r
+ ncl += nClLayer;\r
+ for(Int_t incl = 0; incl < nClLayer; incl++){\r
+ if((fTB[itb].GetCluster(incl))->IsUsed()) nused++;\r
+ }\r
+ }\r
+ \r
+ // calculate the deviation of the mean number of clusters from the\r
+ // closest integer values\r
+ Float_t nclMed = float(ncl-nused)/AliTRDtrackerV1::GetNTimeBins();\r
+ Int_t ncli = Int_t(nclMed);\r
+ Float_t nclDev = TMath::Abs(nclMed - TMath::Max(ncli, 1));\r
+ nclDev -= (nclDev>.5) && ncli ? 1. : 0.;\r
+ return TMath::Exp(-5.*TMath::Abs(nclDev));\r
+\r
+// // get slope of the derivative\r
+// if(!fitter.Eval()) return quality;\r
+// fitter.PrintResults(3);\r
+// Double_t a = fitter.GetParameter(1);\r
+// \r
+// printf("ncl_dev(%f) a(%f)\n", ncl_dev, a);\r
+// return quality*TMath::Exp(-a);\r
+\r
+}\r
+\r
+\r
+//_______________________________________________________\r
+AliTRDchamberTimeBin *AliTRDtrackingChamber::GetSeedingLayer(AliTRDgeometry *geo)\r
+{\r
+ //\r
+ // Creates a seeding layer\r
+ //\r
+ \r
+ // constants\r
+ const Int_t kMaxRows = 16;\r
+ const Int_t kMaxCols = 144;\r
+ const Int_t kMaxPads = 2304;\r
+ \r
+ // Get the geometrical data of the chamber\r
+ Int_t plane = geo->GetPlane(fDetector);\r
+ Int_t stack = geo->GetChamber(fDetector);\r
+ Int_t sector= geo->GetSector(fDetector);\r
+ AliTRDpadPlane *pp = geo->GetPadPlane(plane, stack);\r
+ Int_t nCols = pp->GetNcols();\r
+ Float_t ymin = TMath::Min(pp->GetCol0(), pp->GetColEnd());\r
+ Float_t ymax = TMath::Max(pp->GetCol0(), pp->GetColEnd());\r
+ Float_t zmin = TMath::Min(pp->GetRow0(), pp->GetRowEnd());\r
+ Float_t zmax = TMath::Max(pp->GetRow0(), pp->GetRowEnd());\r
+ Float_t z0 = -1., zl = -1.;\r
+ Int_t nRows = pp->GetNrows();\r
+ Float_t binlength = (ymax - ymin)/nCols; \r
+ //AliInfo(Form("ymin(%f) ymax(%f) zmin(%f) zmax(%f) nRows(%d) binlength(%f)", ymin, ymax, zmin, zmax, nRows, binlength));\r
+ \r
+ // Fill the histogram\r
+ Int_t nClusters; \r
+ Int_t *histogram[kMaxRows]; // 2D-Histogram\r
+ Int_t hvals[kMaxPads]; memset(hvals, 0, sizeof(Int_t)*kMaxPads); \r
+ Float_t *sigmas[kMaxRows];\r
+ Float_t svals[kMaxPads]; memset(svals, 0, sizeof(Float_t)*kMaxPads); \r
+ AliTRDcluster *c = 0x0;\r
+ for(Int_t irs = 0; irs < kMaxRows; irs++){\r
+ histogram[irs] = &hvals[irs*kMaxCols];\r
+ sigmas[irs] = &svals[irs*kMaxCols];\r
+ }\r
+ for(Int_t iTime = 0; iTime < kNTimeBins; iTime++){\r
+ if(!(nClusters = fTB[iTime].GetNClusters())) continue;\r
+ z0 = fTB[iTime].GetZ0();\r
+ zl = fTB[iTime].GetDZ0();\r
+ for(Int_t incl = 0; incl < nClusters; incl++){\r
+ c = fTB[iTime].GetCluster(incl); \r
+ histogram[c->GetPadRow()][c->GetPadCol()]++;\r
+ sigmas[c->GetPadRow()][c->GetPadCol()] += c->GetSigmaZ2();\r
+ }\r
+ }\r
+ \r
+// Now I have everything in the histogram, do the selection\r
+ //Int_t nPads = nCols * nRows;\r
+ // This is what we are interested in: The center of gravity of the best candidates\r
+ Float_t cogyvals[kMaxPads]; memset(cogyvals, 0, sizeof(Float_t)*kMaxPads);\r
+ Float_t cogzvals[kMaxPads]; memset(cogzvals, 0, sizeof(Float_t)*kMaxPads);\r
+ Float_t *cogy[kMaxRows];\r
+ Float_t *cogz[kMaxRows];\r
+ \r
+ // Lookup-Table storing coordinates according to the bins\r
+ Float_t yLengths[kMaxCols];\r
+ Float_t zLengths[kMaxRows];\r
+ for(Int_t icnt = 0; icnt < nCols; icnt++){\r
+ yLengths[icnt] = pp->GetColPos(nCols - 1 - icnt) + binlength/2;\r
+ }\r
+ for(Int_t icnt = 0; icnt < nRows; icnt++){\r
+ zLengths[icnt] = pp->GetRowPos(icnt) - pp->GetRowSize(icnt)/2;\r
+ }\r
+\r
+ // A bitfield is used to mask the pads as usable\r
+ Short_t mask[kMaxCols]; memset(mask, 0 ,sizeof(Short_t) * kMaxCols);//bool mvals[kMaxPads];\r
+ for(UChar_t icount = 0; icount < nRows; icount++){\r
+ cogy[icount] = &cogyvals[icount*kMaxCols];\r
+ cogz[icount] = &cogzvals[icount*kMaxCols];\r
+ }\r
+ // In this array the array position of the best candidates will be stored\r
+ Int_t cand[AliTRDtrackerV1::kMaxTracksStack];\r
+ Float_t sigcands[AliTRDtrackerV1::kMaxTracksStack];\r
+ \r
+ // helper variables\r
+ Int_t indices[kMaxPads]; memset(indices, -1, sizeof(Int_t)*kMaxPads);\r
+ Int_t nCandidates = 0;\r
+ Float_t norm, cogv;\r
+ // histogram filled -> Select best bins\r
+ Int_t nPads = nCols * nRows;\r
+ TMath::Sort(nPads, hvals, indices); // bins storing a 0 should not matter\r
+ // Set Threshold\r
+ Int_t maximum = hvals[indices[0]]; // best\r
+ Int_t threshold = Int_t(maximum * AliTRDReconstructor::RecoParam()->GetFindableClusters());\r
+ Int_t col, row, lower, lower1, upper, upper1;\r
+ for(Int_t ib = 0; ib < nPads; ib++){\r
+ if(nCandidates >= AliTRDtrackerV1::kMaxTracksStack){\r
+ printf("Number of seed candidates %d exceeded maximum allowed per stack %d", nCandidates, AliTRDtrackerV1::kMaxTracksStack);\r
+ break;\r
+ }\r
+ // Positions\r
+ row = indices[ib]/nCols;\r
+ col = indices[ib]%nCols;\r
+ // here will be the threshold condition:\r
+ if((mask[col] & (1 << row)) != 0) continue; // Pad is masked: continue\r
+ if(histogram[row][col] < TMath::Max(threshold, 1)){ // of course at least one cluster is needed\r
+ break; // number of clusters below threshold: break;\r
+ } \r
+ // passing: Mark the neighbors\r
+ lower = TMath::Max(col - 1, 0); upper = TMath::Min(col + 2, nCols);\r
+ lower1 = TMath::Max(row - 1, 0); upper1 = TMath::Min(row + 2, nCols);\r
+ for(Int_t ic = lower; ic < upper; ++ic)\r
+ for(Int_t ir = lower1; ir < upper1; ++ir){\r
+ if(ic == col && ir == row) continue;\r
+ mask[ic] |= (1 << ir);\r
+ }\r
+ // Storing the position in an array\r
+ // testing for neigboring\r
+ cogv = 0;\r
+ norm = 0;\r
+ lower = TMath::Max(col - 1, 0);\r
+ upper = TMath::Min(col + 2, nCols);\r
+ for(Int_t inb = lower; inb < upper; ++inb){\r
+ cogv += yLengths[inb] * histogram[row][inb];\r
+ norm += histogram[row][inb];\r
+ }\r
+ cogy[row][col] = cogv / norm;\r
+ cogv = 0; norm = 0;\r
+ lower = TMath::Max(row - 1, 0);\r
+ upper = TMath::Min(row + 2, nRows);\r
+ for(Int_t inb = lower; inb < upper; ++inb){\r
+ cogv += zLengths[inb] * histogram[inb][col];\r
+ norm += histogram[inb][col];\r
+ }\r
+ cogz[row][col] = Float_t(cogv) / norm;\r
+ // passed the filter\r
+ cand[nCandidates] = row*nCols + col; // store the position of a passig candidate into an Array\r
+ sigcands[nCandidates] = sigmas[row][col] / histogram[row][col]; // never be a floating point exeption\r
+ // Analysis output\r
+ nCandidates++;\r
+ }\r
+ if(!nCandidates) return 0x0;\r
+ \r
+ Float_t pos[3], sig[2];\r
+ Short_t signal[7]; memset(&signal[0], 0, 7*sizeof(Short_t));\r
+ AliTRDchamberTimeBin *fakeLayer = new AliTRDchamberTimeBin(plane, stack, sector, z0, zl);\r
+ AliTRDcluster *cluster = 0x0;\r
+ if(nCandidates){\r
+ UInt_t fakeIndex = 0;\r
+ for(Int_t ican = 0; ican < nCandidates; ican++){\r
+ row = cand[ican] / nCols;\r
+ col = cand[ican] % nCols;\r
+ //temporary\r
+ Int_t n = 0; Double_t x = 0., y = 0., z = 0.;\r
+ for(int itb=0; itb<kNTimeBins; itb++){\r
+ if(!(nClusters = fTB[itb].GetNClusters())) continue;\r
+ for(Int_t incl = 0; incl < nClusters; incl++){\r
+ c = fTB[itb].GetCluster(incl); \r
+ if(c->GetPadRow() != row) continue;\r
+ if(TMath::Abs(c->GetPadCol() - col) > 2) continue;\r
+ x += c->GetX();\r
+ y += c->GetY();\r
+ z += c->GetZ();\r
+ n++;\r
+ }\r
+ }\r
+ pos[0] = x/n;\r
+ pos[1] = y/n;\r
+ pos[2] = z/n;\r
+ sig[0] = .02;\r
+ sig[1] = sigcands[ican];\r
+ cluster = new AliTRDcluster(fDetector, 0., pos, sig, 0x0, 3, signal, col, row, 0, 0, 0., 0);\r
+ fakeLayer->InsertCluster(cluster, fakeIndex++);\r
+ }\r
+ }\r
+ fakeLayer->SetNRows(nRows);\r
+ fakeLayer->SetOwner();\r
+ fakeLayer->BuildIndices();\r
+ //fakeLayer->PrintClusters();\r
+ \r
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel() >= 3){\r
+ //TMatrixD hist(nRows, nCols);\r
+ //for(Int_t i = 0; i < nRows; i++)\r
+ // for(Int_t j = 0; j < nCols; j++)\r
+ // hist(i,j) = histogram[i][j];\r
+ TTreeSRedirector &cstreamer = *AliTRDtrackerV1::DebugStreamer();\r
+ cstreamer << "GetSeedingLayer"\r
+ << "plane=" << plane\r
+ << "ymin=" << ymin\r
+ << "ymax=" << ymax\r
+ << "zmin=" << zmin\r
+ << "zmax=" << zmax\r
+ << "L.=" << fakeLayer\r
+ //<< "Histogram.=" << &hist\r
+ << "\n";\r
+ }\r
+ \r
+ return fakeLayer;\r
+}\r
+\r
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff