From: cblume Date: Mon, 2 Jun 2008 15:36:18 +0000 (+0000) Subject: Enlarging the recoParam functionality and populating the track list (needed by HLT) X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=d20df6fcfae893dc566bc6e74585136316d7c1d7;p=u%2Fmrichter%2FAliRoot.git Enlarging the recoParam functionality and populating the track list (needed by HLT) --- diff --git a/TRD/AliTRDReconstructor.cxx b/TRD/AliTRDReconstructor.cxx index 4769387e31c..8f39bb9f056 100644 --- a/TRD/AliTRDReconstructor.cxx +++ b/TRD/AliTRDReconstructor.cxx @@ -41,11 +41,8 @@ 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; diff --git a/TRD/AliTRDReconstructor.h b/TRD/AliTRDReconstructor.h index 710314188ff..b887ecf4e9d 100644 --- a/TRD/AliTRDReconstructor.h +++ b/TRD/AliTRDReconstructor.h @@ -40,17 +40,11 @@ class AliTRDReconstructor: public AliReconstructor , 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 diff --git a/TRD/AliTRDchamberTimeBin.cxx b/TRD/AliTRDchamberTimeBin.cxx index 7e0fced095b..31b4f0801bb 100644 --- a/TRD/AliTRDchamberTimeBin.cxx +++ b/TRD/AliTRDchamberTimeBin.cxx @@ -1,662 +1,662 @@ -/************************************************************************** - * 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 // -// Markus Fasel // -// // -/////////////////////////////////////////////////////////////////////////////// - -#include -#include -#include -#include -#include - -#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; iGetY()); - 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; irGetY()) 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((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((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; icGetZ(), 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 -// -//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"); - } -} +/************************************************************************** + * 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 // +// Markus Fasel // +// // +/////////////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#include +#include + +#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; iGetY()); + 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::RecoParam()->GetStreamLevel() >= 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; irGetY()) 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((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((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; icGetZ(), 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 +// +//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"); + } +} diff --git a/TRD/AliTRDrecoParam.cxx b/TRD/AliTRDrecoParam.cxx index 87a9410dda2..d24992efda5 100644 --- a/TRD/AliTRDrecoParam.cxx +++ b/TRD/AliTRDrecoParam.cxx @@ -48,6 +48,9 @@ AliTRDrecoParam::AliTRDrecoParam() ,fkChi2Z(30./*14.*//*12.5*/) ,fkChi2Y(.25) ,fkTrackLikelihood(-15.) + ,fkStreamLevel(0) + ,fSeedingOn(kFALSE) + ,fVertexConstrained(kTRUE) ,fClusMaxThresh(4.5) ,fClusSigThresh(3.5) ,fLUTOn(kTRUE) @@ -93,6 +96,8 @@ AliTRDrecoParam *AliTRDrecoParam::GetCosmicTestParam() AliTRDrawStreamBase::SetRawStreamVersion("TB"); AliTRDrecoParam *par = new AliTRDrecoParam(); par->SetADCbaseline(10); + par->SetSeedingOn(kTRUE); + par->SetVertexConstrained(kTRUE); return par; } diff --git a/TRD/AliTRDrecoParam.h b/TRD/AliTRDrecoParam.h index fbc5c062e8a..69d15750293 100644 --- a/TRD/AliTRDrecoParam.h +++ b/TRD/AliTRDrecoParam.h @@ -46,7 +46,11 @@ class AliTRDrecoParam : public AliDetectorRecoParam 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; }; @@ -54,13 +58,16 @@ class AliTRDrecoParam : public AliDetectorRecoParam 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; }; @@ -90,6 +97,10 @@ class AliTRDrecoParam : public AliDetectorRecoParam 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 diff --git a/TRD/AliTRDseedV1.cxx b/TRD/AliTRDseedV1.cxx index 04af0f3a5c4..3f45e8a660d 100644 --- a/TRD/AliTRDseedV1.cxx +++ b/TRD/AliTRDseedV1.cxx @@ -355,7 +355,7 @@ Bool_t AliTRDseedV1::AttachClustersIter(AliTRDtrackingChamber *chamber, Float_t } AliTRDchamberTimeBin *layer = 0x0; - if(AliTRDReconstructor::StreamLevel()>=7 && c){ + if(AliTRDReconstructor::RecoParam()->GetStreamLevel()>=7 && c){ TClonesArray clusters("AliTRDcluster", 24); clusters.SetOwner(kTRUE); AliTRDcluster *cc = 0x0; @@ -423,7 +423,7 @@ Bool_t AliTRDseedV1::AttachClustersIter(AliTRDtrackingChamber *chamber, Float_t 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) @@ -464,7 +464,7 @@ Bool_t AliTRDseedV1::AttachClustersIter(AliTRDtrackingChamber *chamber, Float_t 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); diff --git a/TRD/AliTRDtracker.cxx b/TRD/AliTRDtracker.cxx index b0926794f44..e0d182b2947 100644 --- a/TRD/AliTRDtracker.cxx +++ b/TRD/AliTRDtracker.cxx @@ -50,6 +50,7 @@ #include "AliTRDCommonParam.h" #include "AliTRDtracker.h" #include "AliTRDReconstructor.h" +#include "AliTRDrecoParam.h" #include "AliTRDCalibraFillHisto.h" ClassImp(AliTRDtracker) @@ -417,7 +418,7 @@ Int_t AliTRDtracker::PropagateBack(AliESDEvent *event) Int_t nSeed = event->GetNumberOfTracks(); if(!nSeed){ // run stand alone tracking - if (AliTRDReconstructor::SeedingOn()) Clusters2Tracks(event); + if (AliTRDReconstructor::RecoParam()->SeedingOn()) Clusters2Tracks(event); return 0; } @@ -525,7 +526,7 @@ Int_t AliTRDtracker::PropagateBack(AliESDEvent *event) // 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 @@ -773,7 +774,7 @@ Int_t AliTRDtracker::RefitInward(AliESDEvent *event) // 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; @@ -781,7 +782,7 @@ Int_t AliTRDtracker::RefitInward(AliESDEvent *event) } // Add TRD track to ESDfriendTrack - if (AliTRDReconstructor::StreamLevel() > 0) { + if (AliTRDReconstructor::RecoParam()->GetStreamLevel() > 0) { seed->AddCalibObject(new AliTRDtrack(*pt/*, kTRUE*/)); } delete pt; @@ -1559,7 +1560,7 @@ Int_t AliTRDtracker::Clusters2Tracks(AliESDEvent *esd) isFake = kTRUE; } - if (AliTRDReconstructor::StreamLevel() > 0) { + if (AliTRDReconstructor::RecoParam()->GetStreamLevel() > 0) { if ((!isFake) || ((icl3%10) == 0)) { // Debugging print TTreeSRedirector &cstream = *fDebugStreamer; cstream << "Seeds0" @@ -2197,7 +2198,7 @@ Int_t AliTRDtracker::Clusters2Tracks(AliESDEvent *esd) 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 @@ -2418,7 +2419,7 @@ Int_t AliTRDtracker::Clusters2Tracks(AliESDEvent *esd) 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 @@ -2428,7 +2429,7 @@ Int_t AliTRDtracker::Clusters2Tracks(AliESDEvent *esd) } } - if (AliTRDReconstructor::StreamLevel() > 0) { + if (AliTRDReconstructor::RecoParam()->GetStreamLevel() > 0) { cstream << "Seeds2" << "Iter=" << jter << "Track.=" << track @@ -3661,7 +3662,7 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 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 diff --git a/TRD/AliTRDtrackerV1.cxx b/TRD/AliTRDtrackerV1.cxx index 2768098f345..632aaec98d6 100644 --- a/TRD/AliTRDtrackerV1.cxx +++ b/TRD/AliTRDtrackerV1.cxx @@ -94,11 +94,17 @@ AliTRDtrackerV1::AliTRDtrackerV1() 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(); @@ -252,7 +258,7 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event) Int_t nSeed = event->GetNumberOfTracks(); if(!nSeed){ // run stand alone tracking - if (AliTRDReconstructor::SeedingOn()) Clusters2Tracks(event); + if (AliTRDReconstructor::RecoParam()->SeedingOn()) Clusters2Tracks(event); return 0; } @@ -299,7 +305,7 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event) 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); @@ -388,7 +394,7 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event) 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); @@ -402,7 +408,7 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event) 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); @@ -556,7 +562,7 @@ Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t) } } - 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); @@ -725,7 +731,7 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t) } // 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); @@ -879,7 +885,7 @@ Float_t AliTRDtrackerV1::FitTiltedRiemanConstraint(AliTRDseedV1 *tracklets, Doub 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); @@ -1051,7 +1057,7 @@ Float_t AliTRDtrackerV1::FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigErro 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(); @@ -1283,7 +1289,7 @@ Double_t AliTRDtrackerV1::FitRiemanTilt(AliTRDtrackV1 *track, AliTRDseedV1 *trac } } - if(AliTRDReconstructor::StreamLevel() >=5){ + if(AliTRDReconstructor::RecoParam()->GetStreamLevel() >=5){ TTreeSRedirector &cstreamer = *fgDebugStreamer; Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber(); @@ -1568,6 +1574,30 @@ AliTRDseedV1* AliTRDtrackerV1::SetTracklet(AliTRDseedV1 *tracklet) 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) @@ -1661,7 +1691,7 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClon // 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)); } @@ -1680,7 +1710,7 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClon 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; @@ -1844,7 +1874,7 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClon 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; @@ -1903,7 +1933,7 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClon 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(); @@ -1931,7 +1961,7 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClon 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 @@ -2076,7 +2106,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss 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)); } @@ -2138,7 +2168,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss } 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; @@ -2256,7 +2286,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss // 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(); @@ -2288,7 +2318,10 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss // 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]); @@ -2320,7 +2353,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss 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(); @@ -2395,24 +2428,20 @@ AliTRDtrackV1* AliTRDtrackerV1::MakeTrack(AliTRDseedV1 *seeds, Double_t *params) 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; } @@ -2473,7 +2502,7 @@ Int_t AliTRDtrackerV1::ImproveSeedQuality(AliTRDtrackingChamber **stack, AliTRDs } 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(); @@ -2526,12 +2555,13 @@ Double_t AliTRDtrackerV1::CalculateTrackLikelihood(AliTRDseedV1 *tracklets, Doub 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; @@ -2598,7 +2628,7 @@ Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed, Int_t planes[4] 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 diff --git a/TRD/AliTRDtrackerV1.h b/TRD/AliTRDtrackerV1.h index 042362d42a2..7fed733919e 100644 --- a/TRD/AliTRDtrackerV1.h +++ b/TRD/AliTRDtrackerV1.h @@ -130,6 +130,7 @@ protected: 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: diff --git a/TRD/AliTRDtrackingChamber.cxx b/TRD/AliTRDtrackingChamber.cxx index 22629309509..5d3d1aeed42 100644 --- a/TRD/AliTRDtrackingChamber.cxx +++ b/TRD/AliTRDtrackingChamber.cxx @@ -1,350 +1,350 @@ -/************************************************************************** - * 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 // -// Markus Fasel // -// // -//////////////////////////////////////////////////////////////////////////// - -#include "AliTRDtrackingChamber.h" - -#include "TMath.h" -#include "TMatrixTBase.h" -#include - -#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; itbGetLocalTimeBin()].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; itbGetT0Average(fDetector)); - return kTRUE; -} - -//_______________________________________________________ -Int_t AliTRDtrackingChamber::GetNClusters() const -{ -// Returns number of clusters in chamber -// - Int_t n = 0; - for(Int_t itb=0; itbIsUsed()) 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; itbGetPadRow() != 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; -} - +/************************************************************************** + * 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 // +// Markus Fasel // +// // +//////////////////////////////////////////////////////////////////////////// + +#include "AliTRDtrackingChamber.h" + +#include "TMath.h" +#include "TMatrixTBase.h" +#include + +#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; itbGetLocalTimeBin()].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; itbGetT0Average(fDetector)); + return kTRUE; +} + +//_______________________________________________________ +Int_t AliTRDtrackingChamber::GetNClusters() const +{ +// Returns number of clusters in chamber +// + Int_t n = 0; + for(Int_t itb=0; itbIsUsed()) 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; itbGetPadRow() != 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::RecoParam()->GetStreamLevel() >= 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; +} +