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

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

/* $Id: AliTRDtrackingChamber.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"
#include "AliTRDCommonParam.h"
#include "Cal/AliTRDCalDet.h"
#include "Cal/AliTRDCalROC.h"

ClassImp(AliTRDtrackingChamber)

//_______________________________________________________
AliTRDtrackingChamber::AliTRDtrackingChamber() 
  :TObject()
  ,fDetector(-1)
  ,fX0(0.)
  ,fExB(0.)
  ,fVD(0.)
  ,fT0(0.)
  ,fS2PRF(0.)
  ,fDiffL(0.)
  ,fDiffT(0.)
{}  

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

//_______________________________________________________
void AliTRDtrackingChamber::InsertCluster(AliTRDcluster *c, Int_t index, Bool_t hlt)
{
// Add cluster to TB container and recalculate error parameterization (for HLT)

  fTB[c->GetPadTime()].InsertCluster(c, index);
  if(!hlt) return;

  // Define approximate error parameterization for HLT clusters
  // if needed the fix values of 
  //  - drift length of 1.5 can be replaced with c->GetXloc()
  //  - pad length can be cached from geometry in Init()
  c->SetSigmaY2(fS2PRF, fDiffT, fExB, 1.5);
  c->SetSigmaZ2(6.75);
}

//_______________________________________________________
Bool_t AliTRDtrackingChamber::Build(AliTRDgeometry *const 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)
//
  if(fDetector < 0 || fDetector >= AliTRDgeometry::kNdet){
    AliWarning(Form("Detector index not set correctly to %d", fDetector));
    return kFALSE;
  }

  Int_t stack = AliTRDgeometry::GetStack(fDetector);
  Int_t layer = AliTRDgeometry::GetLayer(fDetector);
  AliTRDpadPlane *pp = geo->GetPadPlane(layer, stack);
  Double_t zl = pp->GetRow0ROC() - pp->GetRowEndROC();
  Double_t z0 = geo->GetRow0(layer, stack, 0) - zl;
  Int_t nrows = pp->GetNrows();
  
  Int_t index[50], jtb = 0;
  for(Int_t itb=0; itb<AliTRDseedV1::kNtb; 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
  Int_t t0 = Int_t(fT0);
  fTB[t0].SetT0();
  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] - t0);	
  return kTRUE;
}
  

//_______________________________________________________
void AliTRDtrackingChamber::Init(Int_t det)
{
// Init detector number and Cache calibration parameters

  fDetector = det;
  AliTRDcalibDB *calib = AliTRDcalibDB::Instance();
  fT0    = calib->GetT0Average(fDetector);
  fVD    = calib->GetVdriftAverage(fDetector);
  fS2PRF = calib->GetPRFROC(fDetector)->GetMean(); fS2PRF *= fS2PRF;
  fExB   = AliTRDCommonParam::Instance()->GetOmegaTau(fVD);
  AliTRDCommonParam::Instance()->GetDiffCoeff(fDiffL,
  fDiffT, fVD);
}

//_______________________________________________________	
Int_t AliTRDtrackingChamber::GetNClusters() const
{
// Basic loop method
// Returns number of clusters in chamber
//
  Int_t n = 0;
  for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++){ 
    n += Int_t(fTB[itb]);
  }
  return n;	
}	

//_______________________________________________________
void AliTRDtrackingChamber::Bootstrap(const AliTRDReconstructor *rec)
{
// Basic loop method
// Bootstrap each time bin
//
  AliTRDchamberTimeBin *jtb = &fTB[0];
  for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++, ++jtb){ 
    (*jtb).Bootstrap(rec, fDetector);
  }
}

//_______________________________________________________
void  AliTRDtrackingChamber::SetOwner()
{
// Basic loop method
// Set ownership in time bins
//
  AliTRDchamberTimeBin *jtb = &fTB[0];
  for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++, ++jtb){ 
    if(!(Int_t(*jtb))) continue;
    (*jtb).SetOwner();
  }
}

//_______________________________________________________
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<AliTRDseedV1::kNtb; 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);

}


//_______________________________________________________
Bool_t AliTRDtrackingChamber::GetSeedingLayer(AliTRDchamberTimeBin *&fakeLayer, AliTRDgeometry * const geo, const AliTRDReconstructor *rec)
{
  //
  // Creates a seeding layer
  //
  
  // constants
  const Int_t kMaxRows = 16;
  const Int_t kMaxCols = 144;
  const Int_t kMaxPads = 2304;
  Int_t timeBinMin = rec->GetRecoParam()->GetNumberOfPresamples();
  Int_t timeBinMax = rec->GetRecoParam()->GetNumberOfPostsamples();

  // Get the geometrical data of the chamber
  Int_t layer = geo->GetLayer(fDetector);
  Int_t stack = geo->GetStack(fDetector);
  Int_t sector= geo->GetSector(fDetector);
  AliTRDpadPlane *pp = geo->GetPadPlane(layer, 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 + 1];	memset(hvals, 0, sizeof(Int_t)*kMaxPads);	 // one entry in addition for termination flag
  Float_t *sigmas[kMaxRows];
  Float_t svals[kMaxPads];	memset(svals, 0, sizeof(Float_t)*kMaxPads);	
  AliTRDcluster *c = NULL;
  for(Int_t irs = 0; irs < kMaxRows; irs++){
    histogram[irs] = &hvals[irs*kMaxCols];
    sigmas[irs] = &svals[irs*kMaxCols];
  }
  for(Int_t iTime = timeBinMin; iTime < AliTRDseedV1::kNtb-timeBinMax; 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;
  // take out all the bins which have less than 3 entries (faster sorting)
  Int_t content[kMaxPads], dictionary[kMaxPads], nCont = 0, padnumber = 0;
  Int_t *iter = &hvals[0], *citer = &content[0], *diter =  &dictionary[0]; // iterators for preselection
  const Int_t threshold = 2;
  hvals[nPads] = -1; // termination for iterator
  do{
    if(*iter > threshold){
      *(citer++) = *iter;
      *(diter++) = padnumber;
      nCont++;
    }
    padnumber++;
  }while(*(++iter) != -1);
  TMath::Sort(nCont, content, indices);		

  Int_t col, row, lower, lower1, upper, upper1;
  for(Int_t ib = 0; ib < nCont; ib++){
    if(nCandidates >= AliTRDtrackerV1::kMaxTracksStack){
      AliWarning(Form("Number of seed candidates %d exceeded maximum allowed per stack %d", nCandidates, AliTRDtrackerV1::kMaxTracksStack));
      break;
    }
    // Positions
    row = dictionary[indices[ib]]/nCols;
    col = dictionary[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 kFALSE;
  
  Float_t pos[3], sig[2];
  Short_t signal[7]; memset(&signal[0], 0, 7*sizeof(Short_t));
  
  new(fakeLayer) AliTRDchamberTimeBin(layer, stack, sector, z0, zl);
  fakeLayer->SetReconstructor(rec);
  fakeLayer->SetNRows(nRows);
  fakeLayer->SetOwner(kFALSE);
  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<AliTRDseedV1::kNtb; 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];
      fakeLayer->InsertCluster(new AliTRDcluster(fDetector, 0., pos, sig, NULL, 3, signal, col, row, 0, 0, 0., 0), fakeIndex++);
    }
  }
  fakeLayer->BuildIndices();
  //fakeLayer->Print();
  
  if(rec->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) >= 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 = *rec->GetDebugStream(AliTRDrecoParam::kTracker);
    cstreamer << "GetSeedingLayer"
    << "layer="      << layer
    << "ymin="       << ymin
    << "ymax="       << ymax
    << "zmin="       << zmin
    << "zmax="       << zmax
    << "L.="         << fakeLayer
    //<< "Histogram.=" << &hist
    << "\n";
  }
  
  return kTRUE;
}


//_______________________________________________________
void AliTRDtrackingChamber::Print(Option_t *opt) const
{
  // Print the chamber status
  if(!GetNClusters()) return;
  AliInfo(Form("fDetector   = %d", fDetector));
  AliInfo(Form("fX0         = %7.3f", fX0));
  const AliTRDchamberTimeBin *itb = &fTB[0];
  for(Int_t jtb=0; jtb<AliTRDseedV1::kNtb; jtb++, itb++) (*itb).Print(opt);
}


//_______________________________________________________
void AliTRDtrackingChamber::Update()
{
// Steer purging of used and shared clusters 

  AliTRDchamberTimeBin *jtb = &fTB[0];
  for(Int_t itb=AliTRDseedV1::kNtb; itb--; ++jtb){ 
    if(!(Int_t(*jtb))) continue;
    (*jtb).BuildIndices();
  }
}
Commit	Line	Data
b0a48c4d	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id: AliTRDtrackingChamber.cxx 23810 2008-02-08 09:00:27Z hristov $ */
	17
	18	////////////////////////////////////////////////////////////////////////////
	19	// //
	20	// Tracking in one chamber //
	21	// //
	22	// Authors: //
	23	// Alex Bercuci <A.Bercuci@gsi.de> //
	24	// Markus Fasel <M.Fasel@gsi.de> //
	25	// //
	26	////////////////////////////////////////////////////////////////////////////
	27
	28	#include "AliTRDtrackingChamber.h"
	29
	30	#include "TMath.h"
	31	#include "TMatrixTBase.h"
	32	#include <TTreeStream.h>
	33
	34	#include "AliTRDReconstructor.h"
	35	#include "AliTRDrecoParam.h"
	36	#include "AliTRDtrackerV1.h"
	37	#include "AliTRDgeometry.h"
	38	#include "AliTRDpadPlane.h"
	39	#include "AliTRDcalibDB.h"
c79857d5	40	#include "AliTRDCommonParam.h"
b0a48c4d	41	#include "Cal/AliTRDCalDet.h"
	42	#include "Cal/AliTRDCalROC.h"
	43
	44	ClassImp(AliTRDtrackingChamber)
	45
	46	//_______________________________________________________
349f5eeb	47	AliTRDtrackingChamber::AliTRDtrackingChamber()
a8276d32	48	:TObject()
349f5eeb	49	,fDetector(-1)
b0a48c4d	50	,fX0(0.)
c79857d5	51	,fExB(0.)
	52	,fVD(0.)
	53	,fT0(0.)
	54	,fS2PRF(0.)
	55	,fDiffL(0.)
	56	,fDiffT(0.)
b0a48c4d	57	{}
	58
	59	//_______________________________________________________
	60	void AliTRDtrackingChamber::Clear(const Option_t *opt)
	61	{
fac58f00	62	for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++) fTB[itb].Clear(opt);
b0a48c4d	63	}
	64
	65	//_______________________________________________________
c79857d5	66	void AliTRDtrackingChamber::InsertCluster(AliTRDcluster *c, Int_t index, Bool_t hlt)
b0a48c4d	67	{
c79857d5	68	// Add cluster to TB container and recalculate error parameterization (for HLT)
c79857d5	69
b0a48c4d	70	fTB[c->GetPadTime()].InsertCluster(c, index);
c79857d5	71	if(!hlt) return;
	72
	73	// Define approximate error parameterization for HLT clusters
	74	// if needed the fix values of
	75	// - drift length of 1.5 can be replaced with c->GetXloc()
	76	// - pad length can be cached from geometry in Init()
	77	c->SetSigmaY2(fS2PRF, fDiffT, fExB, 1.5);
	78	c->SetSigmaZ2(6.75);
b0a48c4d	79	}
	80
	81	//_______________________________________________________
c79857d5	82	Bool_t AliTRDtrackingChamber::Build(AliTRDgeometry *const geo)
b0a48c4d	83	{
	84	// Init chamber and all time bins (AliTRDchamberTimeBin)
	85	// Calculates radial position of the chamber based on
	86	// radial positions of the time bins (calibration/alignment aware)
	87	//
349f5eeb	88	if(fDetector < 0 \|\| fDetector >= AliTRDgeometry::kNdet){
	89	AliWarning(Form("Detector index not set correctly to %d", fDetector));
	90	return kFALSE;
	91	}
	92
c79857d5	93	Int_t stack = AliTRDgeometry::GetStack(fDetector);
c79857d5	94	Int_t layer = AliTRDgeometry::GetLayer(fDetector);
b0a48c4d	95	AliTRDpadPlane *pp = geo->GetPadPlane(layer, stack);
	96	Double_t zl = pp->GetRow0ROC() - pp->GetRowEndROC();
	97	Double_t z0 = geo->GetRow0(layer, stack, 0) - zl;
	98	Int_t nrows = pp->GetNrows();
	99
	100	Int_t index[50], jtb = 0;
fac58f00	101	for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++){
b0a48c4d	102	if(!fTB[itb]) continue;
	103	fTB[itb].SetRange(z0, zl);
	104	fTB[itb].SetNRows(nrows);
	105	fTB[itb].BuildIndices();
	106	index[jtb++] = itb;
	107	}
	108	if(jtb<2) return kFALSE;
	109
b0a48c4d	110	// ESTIMATE POSITION OF PAD PLANE FOR THIS CHAMBER
c79857d5	111	Int_t t0 = Int_t(fT0);
c79857d5	112	fTB[t0].SetT0();
b0a48c4d	113	Double_t x0 = fTB[index[0]].GetX();
	114	Double_t x1 = fTB[index[1]].GetX();
	115	Double_t dx = (x0 - x1)/(index[1] - index[0]);
b0a48c4d	116	fX0 = x0 + dx*(index[0] - t0);
	117	return kTRUE;
	118	}
	119
c79857d5	120
	121	//_______________________________________________________
	122	void AliTRDtrackingChamber::Init(Int_t det)
	123	{
	124	// Init detector number and Cache calibration parameters
	125
	126	fDetector = det;
	127	AliTRDcalibDB *calib = AliTRDcalibDB::Instance();
	128	fT0 = calib->GetT0Average(fDetector);
	129	fVD = calib->GetVdriftAverage(fDetector);
	130	fS2PRF = calib->GetPRFROC(fDetector)->GetMean(); fS2PRF *= fS2PRF;
	131	fExB = AliTRDCommonParam::Instance()->GetOmegaTau(fVD);
	132	AliTRDCommonParam::Instance()->GetDiffCoeff(fDiffL,
	133	fDiffT, fVD);
	134	}
	135
b0a48c4d	136	//_______________________________________________________
	137	Int_t AliTRDtrackingChamber::GetNClusters() const
	138	{
ac1dec3b	139	// Basic loop method
b0a48c4d	140	// Returns number of clusters in chamber
	141	//
	142	Int_t n = 0;
fac58f00	143	for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++){
b0a48c4d	144	n += Int_t(fTB[itb]);
	145	}
	146	return n;
	147	}
	148
ac1dec3b	149	//_______________________________________________________
	150	void AliTRDtrackingChamber::Bootstrap(const AliTRDReconstructor *rec)
	151	{
	152	// Basic loop method
	153	// Bootstrap each time bin
	154	//
	155	AliTRDchamberTimeBin *jtb = &fTB[0];
fac58f00	156	for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++, ++jtb){
ac1dec3b	157	(*jtb).Bootstrap(rec, fDetector);
	158	}
	159	}
	160
	161	//_______________________________________________________
	162	void AliTRDtrackingChamber::SetOwner()
	163	{
	164	// Basic loop method
	165	// Set ownership in time bins
	166	//
	167	AliTRDchamberTimeBin *jtb = &fTB[0];
fac58f00	168	for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++, ++jtb){
ac1dec3b	169	if(!(Int_t(*jtb))) continue;
	170	(*jtb).SetOwner();
	171	}
	172	}
	173
b0a48c4d	174	//_______________________________________________________
	175	Double_t AliTRDtrackingChamber::GetQuality()
	176	{
	177	//
	178	// Calculate chamber quality for seeding.
	179	//
	180	//
	181	// Parameters :
	182	// layers : Array of propagation layers for this plane.
	183	//
	184	// Output :
	185	// plane quality factor for seeding
	186	//
	187	// Detailed description
	188	//
	189	// The quality of the plane for seeding is higher if:
	190	// 1. the average timebin population is closer to an integer number
	191	// 2. the distribution of clusters/timebin is closer to a uniform distribution.
	192	// - the slope of the first derivative of a parabolic fit is small or
	193	// - the slope of a linear fit is small
	194	//
	195
	196	Int_t ncl = 0;
	197	Int_t nused = 0;
	198	Int_t nClLayer;
fac58f00	199	for(int itb=0; itb<AliTRDseedV1::kNtb; itb++){
b0a48c4d	200	if(!(nClLayer = fTB[itb].GetNClusters())) continue;
	201	ncl += nClLayer;
	202	for(Int_t incl = 0; incl < nClLayer; incl++){
	203	if((fTB[itb].GetCluster(incl))->IsUsed()) nused++;
	204	}
	205	}
	206
	207	// calculate the deviation of the mean number of clusters from the
	208	// closest integer values
	209	Float_t nclMed = float(ncl-nused)/AliTRDtrackerV1::GetNTimeBins();
	210	Int_t ncli = Int_t(nclMed);
	211	Float_t nclDev = TMath::Abs(nclMed - TMath::Max(ncli, 1));
	212	nclDev -= (nclDev>.5) && ncli ? 1. : 0.;
	213	return TMath::Exp(-5.*TMath::Abs(nclDev));
	214
	215	// // get slope of the derivative
	216	// if(!fitter.Eval()) return quality;
	217	// fitter.PrintResults(3);
	218	// Double_t a = fitter.GetParameter(1);
	219	//
	220	// printf("ncl_dev(%f) a(%f)\n", ncl_dev, a);
	221	// return quality*TMath::Exp(-a);
	222
	223	}
	224
	225
	226	//_______________________________________________________
4d6aee34	227	Bool_t AliTRDtrackingChamber::GetSeedingLayer(AliTRDchamberTimeBin &fakeLayer, AliTRDgeometry const geo, const AliTRDReconstructor *rec)
b0a48c4d	228	{
	229	//
	230	// Creates a seeding layer
	231	//
	232
	233	// constants
	234	const Int_t kMaxRows = 16;
	235	const Int_t kMaxCols = 144;
	236	const Int_t kMaxPads = 2304;
	237	Int_t timeBinMin = rec->GetRecoParam()->GetNumberOfPresamples();
	238	Int_t timeBinMax = rec->GetRecoParam()->GetNumberOfPostsamples();
	239
	240	// Get the geometrical data of the chamber
	241	Int_t layer = geo->GetLayer(fDetector);
	242	Int_t stack = geo->GetStack(fDetector);
	243	Int_t sector= geo->GetSector(fDetector);
	244	AliTRDpadPlane *pp = geo->GetPadPlane(layer, stack);
	245	Int_t nCols = pp->GetNcols();
	246	Float_t ymin = TMath::Min(pp->GetCol0(), pp->GetColEnd());
	247	Float_t ymax = TMath::Max(pp->GetCol0(), pp->GetColEnd());
	248	Float_t zmin = TMath::Min(pp->GetRow0(), pp->GetRowEnd());
	249	Float_t zmax = TMath::Max(pp->GetRow0(), pp->GetRowEnd());
	250	Float_t z0 = -1., zl = -1.;
	251	Int_t nRows = pp->GetNrows();
	252	Float_t binlength = (ymax - ymin)/nCols;
	253	//AliInfo(Form("ymin(%f) ymax(%f) zmin(%f) zmax(%f) nRows(%d) binlength(%f)", ymin, ymax, zmin, zmax, nRows, binlength));
	254
	255	// Fill the histogram
	256	Int_t nClusters;
	257	Int_t *histogram[kMaxRows]; // 2D-Histogram
	258	Int_t hvals[kMaxPads + 1]; memset(hvals, 0, sizeof(Int_t)*kMaxPads); // one entry in addition for termination flag
	259	Float_t *sigmas[kMaxRows];
	260	Float_t svals[kMaxPads]; memset(svals, 0, sizeof(Float_t)*kMaxPads);
4d6aee34	261	AliTRDcluster *c = NULL;
b0a48c4d	262	for(Int_t irs = 0; irs < kMaxRows; irs++){
	263	histogram[irs] = &hvals[irs*kMaxCols];
	264	sigmas[irs] = &svals[irs*kMaxCols];
	265	}
fac58f00	266	for(Int_t iTime = timeBinMin; iTime < AliTRDseedV1::kNtb-timeBinMax; iTime++){
b0a48c4d	267	if(!(nClusters = fTB[iTime].GetNClusters())) continue;
	268	z0 = fTB[iTime].GetZ0();
	269	zl = fTB[iTime].GetDZ0();
	270	for(Int_t incl = 0; incl < nClusters; incl++){
	271	c = fTB[iTime].GetCluster(incl);
	272	histogram[c->GetPadRow()][c->GetPadCol()]++;
	273	sigmas[c->GetPadRow()][c->GetPadCol()] += c->GetSigmaZ2();
	274	}
	275	}
	276
	277	// Now I have everything in the histogram, do the selection
	278	//Int_t nPads = nCols * nRows;
	279	// This is what we are interested in: The center of gravity of the best candidates
	280	Float_t cogyvals[kMaxPads]; memset(cogyvals, 0, sizeof(Float_t)*kMaxPads);
	281	Float_t cogzvals[kMaxPads]; memset(cogzvals, 0, sizeof(Float_t)*kMaxPads);
	282	Float_t *cogy[kMaxRows];
	283	Float_t *cogz[kMaxRows];
	284
	285	// Lookup-Table storing coordinates according to the bins
	286	Float_t yLengths[kMaxCols];
	287	Float_t zLengths[kMaxRows];
	288	for(Int_t icnt = 0; icnt < nCols; icnt++){
	289	yLengths[icnt] = pp->GetColPos(nCols - 1 - icnt) + binlength/2;
	290	}
	291	for(Int_t icnt = 0; icnt < nRows; icnt++){
	292	zLengths[icnt] = pp->GetRowPos(icnt) - pp->GetRowSize(icnt)/2;
	293	}
	294
	295	// A bitfield is used to mask the pads as usable
	296	Short_t mask[kMaxCols]; memset(mask, 0 ,sizeof(Short_t) * kMaxCols);//bool mvals[kMaxPads];
	297	for(UChar_t icount = 0; icount < nRows; icount++){
	298	cogy[icount] = &cogyvals[icount*kMaxCols];
	299	cogz[icount] = &cogzvals[icount*kMaxCols];
	300	}
	301	// In this array the array position of the best candidates will be stored
	302	Int_t cand[AliTRDtrackerV1::kMaxTracksStack];
	303	Float_t sigcands[AliTRDtrackerV1::kMaxTracksStack];
	304
	305	// helper variables
	306	Int_t indices[kMaxPads]; memset(indices, -1, sizeof(Int_t)*kMaxPads);
	307	Int_t nCandidates = 0;
	308	Float_t norm, cogv;
	309	// histogram filled -> Select best bins
	310	Int_t nPads = nCols * nRows;
	311	// take out all the bins which have less than 3 entries (faster sorting)
	312	Int_t content[kMaxPads], dictionary[kMaxPads], nCont = 0, padnumber = 0;
	313	Int_t iter = &hvals[0], citer = &content[0], *diter = &dictionary[0]; // iterators for preselection
	314	const Int_t threshold = 2;
	315	hvals[nPads] = -1; // termination for iterator
	316	do{
	317	if(*iter > threshold){
	318	(citer++) = iter;
	319	*(diter++) = padnumber;
	320	nCont++;
	321	}
	322	padnumber++;
	323	}while(*(++iter) != -1);
	324	TMath::Sort(nCont, content, indices);
	325
	326	Int_t col, row, lower, lower1, upper, upper1;
	327	for(Int_t ib = 0; ib < nCont; ib++){
	328	if(nCandidates >= AliTRDtrackerV1::kMaxTracksStack){
2a766b25	329	AliWarning(Form("Number of seed candidates %d exceeded maximum allowed per stack %d", nCandidates, AliTRDtrackerV1::kMaxTracksStack));
b0a48c4d	330	break;
	331	}
	332	// Positions
	333	row = dictionary[indices[ib]]/nCols;
	334	col = dictionary[indices[ib]]%nCols;
	335	// here will be the threshold condition:
	336	if((mask[col] & (1 << row)) != 0) continue; // Pad is masked: continue
	337	// if(histogram[row][col] < TMath::Max(threshold, 1)){ // of course at least one cluster is needed
	338	// break; // number of clusters below threshold: break;
	339	// }
	340	// passing: Mark the neighbors
	341	lower = TMath::Max(col - 1, 0); upper = TMath::Min(col + 2, nCols);
	342	lower1 = TMath::Max(row - 1, 0); upper1 = TMath::Min(row + 2, nCols);
	343	for(Int_t ic = lower; ic < upper; ++ic)
	344	for(Int_t ir = lower1; ir < upper1; ++ir){
	345	if(ic == col && ir == row) continue;
	346	mask[ic] \|= (1 << ir);
	347	}
	348	// Storing the position in an array
	349	// testing for neigboring
	350	cogv = 0;
	351	norm = 0;
	352	lower = TMath::Max(col - 1, 0);
	353	upper = TMath::Min(col + 2, nCols);
	354	for(Int_t inb = lower; inb < upper; ++inb){
	355	cogv += yLengths[inb] * histogram[row][inb];
	356	norm += histogram[row][inb];
	357	}
	358	cogy[row][col] = cogv / norm;
	359	cogv = 0; norm = 0;
	360	lower = TMath::Max(row - 1, 0);
	361	upper = TMath::Min(row + 2, nRows);
	362	for(Int_t inb = lower; inb < upper; ++inb){
	363	cogv += zLengths[inb] * histogram[inb][col];
	364	norm += histogram[inb][col];
	365	}
	366	cogz[row][col] = Float_t(cogv) / norm;
	367	// passed the filter
	368	cand[nCandidates] = row*nCols + col; // store the position of a passig candidate into an Array
	369	sigcands[nCandidates] = sigmas[row][col] / histogram[row][col]; // never be a floating point exeption
	370	// Analysis output
	371	nCandidates++;
	372	}
	373	if(!nCandidates) return kFALSE;
	374
	375	Float_t pos[3], sig[2];
	376	Short_t signal[7]; memset(&signal[0], 0, 7*sizeof(Short_t));
	377
	378	new(fakeLayer) AliTRDchamberTimeBin(layer, stack, sector, z0, zl);
	379	fakeLayer->SetReconstructor(rec);
ac1dec3b	380	fakeLayer->SetNRows(nRows);
ac1dec3b	381	fakeLayer->SetOwner(kFALSE);
b0a48c4d	382	if(nCandidates){
	383	UInt_t fakeIndex = 0;
	384	for(Int_t ican = 0; ican < nCandidates; ican++){
	385	row = cand[ican] / nCols;
	386	col = cand[ican] % nCols;
	387	//temporary
	388	Int_t n = 0; Double_t x = 0., y = 0., z = 0.;
fac58f00	389	for(int itb=0; itb<AliTRDseedV1::kNtb; itb++){
b0a48c4d	390	if(!(nClusters = fTB[itb].GetNClusters())) continue;
	391	for(Int_t incl = 0; incl < nClusters; incl++){
	392	c = fTB[itb].GetCluster(incl);
	393	if(c->GetPadRow() != row) continue;
	394	if(TMath::Abs(c->GetPadCol() - col) > 2) continue;
	395	x += c->GetX();
	396	y += c->GetY();
	397	z += c->GetZ();
	398	n++;
	399	}
	400	}
	401	pos[0] = x/n;
	402	pos[1] = y/n;
	403	pos[2] = z/n;
	404	sig[0] = .02;
	405	sig[1] = sigcands[ican];
4d6aee34	406	fakeLayer->InsertCluster(new AliTRDcluster(fDetector, 0., pos, sig, NULL, 3, signal, col, row, 0, 0, 0., 0), fakeIndex++);
b0a48c4d	407	}
b0a48c4d	408	}
b0a48c4d	409	fakeLayer->BuildIndices();
ac1dec3b	410	//fakeLayer->Print();
b0a48c4d	411
a2fbb6ec	412	if(rec->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) >= 3){
b0a48c4d	413	//TMatrixD hist(nRows, nCols);
	414	//for(Int_t i = 0; i < nRows; i++)
	415	// for(Int_t j = 0; j < nCols; j++)
	416	// hist(i,j) = histogram[i][j];
a2fbb6ec	417	TTreeSRedirector &cstreamer = *rec->GetDebugStream(AliTRDrecoParam::kTracker);
b0a48c4d	418	cstreamer << "GetSeedingLayer"
	419	<< "layer=" << layer
	420	<< "ymin=" << ymin
	421	<< "ymax=" << ymax
	422	<< "zmin=" << zmin
	423	<< "zmax=" << zmax
	424	<< "L.=" << fakeLayer
	425	//<< "Histogram.=" << &hist
	426	<< "\n";
	427	}
	428
	429	return kTRUE;
	430	}
	431
6c207d50	432
	433	//_______________________________________________________
	434	void AliTRDtrackingChamber::Print(Option_t *opt) const
	435	{
4d6aee34	436	// Print the chamber status
6c207d50	437	if(!GetNClusters()) return;
	438	AliInfo(Form("fDetector = %d", fDetector));
	439	AliInfo(Form("fX0 = %7.3f", fX0));
	440	const AliTRDchamberTimeBin *itb = &fTB[0];
fac58f00	441	for(Int_t jtb=0; jtb<AliTRDseedV1::kNtb; jtb++, itb++) (*itb).Print(opt);
6c207d50	442	}
fac58f00	443
	444
	445	//_______________________________________________________
	446	void AliTRDtrackingChamber::Update()
	447	{
	448	// Steer purging of used and shared clusters
	449
	450	AliTRDchamberTimeBin *jtb = &fTB[0];
	451	for(Int_t itb=AliTRDseedV1::kNtb; itb--; ++jtb){
	452	if(!(Int_t(*jtb))) continue;
	453	(*jtb).BuildIndices();
	454	}
	455	}
	456