#include "AliTRDgtuTrack.h"
#include "AliTRDrawData.h"
#include "AliTRDdigitsManager.h"
+#include "AliTRDtrackerV1.h"
+#include "AliTRDrecoParam.h"
ClassImp(AliTRDReconstructor)
//
// Reconstruct clusters
//
+
AliInfo("Reconstruct TRD clusters from Digits [Digit TTree -> Cluster TTree]");
AliTRDclusterizer clusterer("clusterer","TRD clusterizer");
// Create a TRD tracker
//
- return new AliTRDtracker(NULL);
+ //return new AliTRDtracker(NULL);
+
+ return new AliTRDtrackerV1(NULL, AliTRDrecoParam::GetLowFluxParam());
}
for (iTime = 0; iTime < nTimeTotal; iTime++) {
// Store the amplitude of the digit if above threshold
- // if (outADC[iTime] > (simParam->GetADCbaseline() + simParam->GetADCthreshold())) {
+ //if (outADC[iTime] > (simParam->GetADCbaseline() + simParam->GetADCthreshold())) {
if (outADC[iTime] != 0 ) { // Now this is enough because there is ZS in raw simulator
nDigits++;
digits->SetDataUnchecked(iRow,iCol,iTime,((Int_t) outADC[iTime]));
for (iTime = 0; iTime < nTimeTotal; iTime++) {
// Store the amplitude of the digit if above threshold
- // if (outADC[iTime] > (adcBaseline + adcThreshold)) {
+ //if (outADC[iTime] > (adcBaseline + adcThreshold)) {
if (outADC[iTime] != 0) { // now this is ok because there is ZS in raw simulation
digitsOut->SetDataUnchecked(iRow,iCol,iTime,((Int_t) outADC[iTime]));
// Copy the dictionary
--- /dev/null
+/**************************************************************************
+ * 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$ */
+
+////////////////////////////////////////////////////////////////////////////
+// //
+// The TRD propagation layer //
+// //
+// Authors: //
+// Marian Ivanov <M.Ivanov@gsi.de> //
+// Alex Bercuci <A.Bercuci@gsi.de> //
+// Markus Fasel <M.Fasel@gsi.de> //
+// //
+////////////////////////////////////////////////////////////////////////////
+
+#include "string.h"
+
+#include "TMath.h"
+
+#include "AliTRDpropagationLayer.h"
+//#include "AliTRDtracker.h"
+#include "AliTRDcluster.h"
+#include "AliTRDgeometry.h"
+
+//_____________________________________________________________________________
+AliTRDpropagationLayer::AliTRDpropagationLayer()
+ :TObject()
+ ,fN(0)
+ ,fSec(0)
+ ,fClusters(NULL)
+ ,fIndex(NULL)
+ ,fX(0.)
+ ,fdX(0.)
+ ,fRho(0.)
+ ,fX0(0.)
+ ,fTimeBinIndex(0)
+ ,fPlane(0)
+ ,fYmax(0)
+ ,fYmaxSensitive(0)
+ ,fHole(kFALSE)
+ ,fHoleZc(0)
+ ,fHoleZmax(0)
+ ,fHoleYc(0)
+ ,fHoleYmax(0)
+ ,fHoleRho(0)
+ ,fHoleX0(0)
+{
+ //
+ // Default constructor
+ //
+
+}
+
+//_____________________________________________________________________________
+AliTRDpropagationLayer::AliTRDpropagationLayer(Double_t x, Double_t dx, Double_t rho
+ , Double_t radLength, Int_t tbIndex, Int_t plane)
+ :TObject()
+ ,fN(0)
+ ,fSec(0)
+ ,fClusters(NULL)
+ ,fIndex(NULL)
+ ,fX(x)
+ ,fdX(dx)
+ ,fRho(rho)
+ ,fX0(radLength)
+ ,fTimeBinIndex(tbIndex)
+ ,fPlane(plane)
+ ,fYmax(0)
+ ,fYmaxSensitive(0)
+ ,fHole(kFALSE)
+ ,fHoleZc(0)
+ ,fHoleZmax(0)
+ ,fHoleYc(0)
+ ,fHoleYmax(0)
+ ,fHoleRho(0)
+ ,fHoleX0(0)
+{
+ //
+ // AliTRDpropagationLayer constructor
+ //
+
+ for (Int_t i = 0; i < (Int_t)kZones; i++) {
+ fZc[i] = 0;
+ fZmax[i] = 0;
+ }
+
+ if (fTimeBinIndex >= 0) {
+ fClusters = new AliTRDcluster*[kMaxClusterPerTimeBin];
+ fIndex = new UInt_t[kMaxClusterPerTimeBin];
+ }
+
+ for (Int_t i = 0; i < 5; i++) {
+ fIsHole[i] = kFALSE;
+ }
+
+}
+
+//_____________________________________________________________________________
+AliTRDpropagationLayer::AliTRDpropagationLayer(const AliTRDpropagationLayer &p)
+ :TObject((TObject&)p)
+ ,fN(p.fN)
+ ,fSec(p.fSec)
+ ,fClusters(0x0)
+ ,fIndex(0x0)
+ ,fX(p.fX)
+ ,fdX(p.fdX)
+ ,fRho(p.fRho)
+ ,fX0(p.fX0)
+ ,fTimeBinIndex(p.fTimeBinIndex)
+ ,fPlane(p.fPlane)
+ ,fYmax(p.fYmax)
+ ,fYmaxSensitive(p.fYmaxSensitive)
+ ,fHole(p.fHole)
+ ,fHoleZc(p.fHoleZc)
+ ,fHoleZmax(p.fHoleZmax)
+ ,fHoleYc(p.fHoleYc)
+ ,fHoleYmax(p.fHoleYmax)
+ ,fHoleRho(p.fHoleRho)
+ ,fHoleX0(p.fHoleX0)
+{
+ //
+ // AliTRDpropagationLayer copy constructor
+ //
+
+ for (Int_t i = 0; i < (Int_t)kZones; i++) {
+ fZc[i] = p.fZc[i];
+ fZmax[i] = p.fZmax[i];
+ fIsHole[i] = p.fIsHole[i];
+ fZmaxSensitive[i] = p.fZmaxSensitive[i];
+ }
+
+ // Make a deep copy of the Clusters array and the Index array unless they are needed in class AliTRDstackLayer
+ Int_t arrsize = (Int_t)kMaxClusterPerTimeBin;
+ if (fTimeBinIndex >= 0) {
+ fClusters = new AliTRDcluster*[arrsize];
+ fIndex = new UInt_t[arrsize];
+ }
+ memset(fIndex, 0, sizeof(UInt_t)*arrsize);
+ memset(fClusters, 0, sizeof(AliTRDcluster *)*arrsize);
+ for(Int_t i = 0; i < arrsize; i++){
+ fClusters[i] = p.fClusters[i];
+ fIndex[i] = p.fIndex[i];
+ }
+}
+
+//_____________________________________________________________________________
+AliTRDpropagationLayer::~AliTRDpropagationLayer()
+{
+ //
+ // Destructor
+ //
+
+ if (fTimeBinIndex >= 0) {
+ delete[] fClusters;
+ delete[] fIndex;
+ }
+
+}
+
+//_____________________________________________________________________________
+void AliTRDpropagationLayer::Copy(TObject &o) const
+{
+ //
+ // Copy function
+ //
+
+ AliTRDpropagationLayer &p = (AliTRDpropagationLayer &)o;
+ p.fN = fN;
+ p.fSec = fSec;
+ p.fX = fX;
+ p.fdX = fdX;
+ p.fRho = fRho;
+ p.fX0 = fX0;
+ p.fTimeBinIndex = fTimeBinIndex;
+ p.fPlane = fPlane;
+ p.fYmax = fYmax;
+ p.fYmaxSensitive = fYmaxSensitive;
+ p.fHole = fHole;
+ p.fHoleZc = fHoleZc;
+ p.fHoleZmax = fHoleZmax;
+ p.fHoleYc = fHoleYc;
+ p.fHoleYmax = fHoleYmax;
+ p.fHoleRho = fHoleRho;
+ p.fHoleX0 = fHoleX0;
+
+ for (Int_t i = 0; i < (Int_t)kZones; i++) {
+ p.fZc[i] = fZc[i];
+ p.fZmax[i] = fZmax[i];
+ p.fIsHole[i] = fIsHole[i];
+ p.fZmaxSensitive[i] = fZmaxSensitive[i];
+ }
+
+ // Make a deep copy of the Clusters array and the Index array
+ // unless they are needed in class AliTRDstackLayer
+ if (fTimeBinIndex >= 0) {
+ if (!p.fClusters)
+ p.fClusters = new AliTRDcluster*[(Int_t)kMaxClusterPerTimeBin];
+ if (!p.fIndex)
+ p.fIndex = new UInt_t[(Int_t)kMaxClusterPerTimeBin];
+ }
+ for (Int_t i = 0; i < (Int_t)kMaxClusterPerTimeBin; i++){
+ //overwrite
+ p.fClusters[i] = fClusters[i];
+ p.fIndex[i] = fIndex[i];
+ }
+
+}
+
+//_____________________________________________________________________________
+void AliTRDpropagationLayer::SetZ(Double_t *center, Double_t *w, Double_t *wsensitive )
+{
+ //
+ // Set centers and the width of sectors
+ //
+
+ for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
+ fZc[icham] = center[icham];
+ fZmax[icham] = w[icham];
+ fZmaxSensitive[icham] = wsensitive[icham];
+ }
+
+}
+
+//_____________________________________________________________________________
+void AliTRDpropagationLayer::SetHoles(Bool_t *holes)
+{
+ //
+ // Set centers and the width of sectors
+ //
+
+ fHole = kFALSE;
+
+ for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
+ fIsHole[icham] = holes[icham];
+ if (holes[icham]) {
+ fHole = kTRUE;
+ }
+ }
+
+}
+
+//_____________________________________________________________________________
+void AliTRDpropagationLayer::InsertCluster(AliTRDcluster *c, UInt_t index)
+{
+ //
+ // Insert cluster in cluster array.
+ // Clusters are sorted according to Y coordinate.
+ //
+
+ if (fTimeBinIndex < 0) {
+ //AliWarning("Attempt to insert cluster into non-sensitive time bin!\n");
+ return;
+ }
+
+ if (fN == (Int_t) kMaxClusterPerTimeBin) {
+ //AliWarning("Too many clusters !\n");
+ return;
+ }
+
+ if (fN == 0) {
+ fIndex[0] = index;
+ fClusters[fN++] = c;
+ return;
+ }
+
+ Int_t i = Find(c->GetY());
+ memmove(fClusters+i+1,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
+ memmove(fIndex +i+1,fIndex +i,(fN-i)*sizeof(UInt_t));
+ fIndex[i] = index;
+ fClusters[i] = c;
+ fN++;
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDpropagationLayer::Find(Float_t y) const
+{
+ //
+ // Returns index of the cluster nearest in Y
+ //
+
+ if (fN <= 0) {
+ return 0;
+ }
+ if (y <= fClusters[0]->GetY()) {
+ return 0;
+ }
+ if (y > fClusters[fN-1]->GetY()) {
+ return fN;
+ }
+
+ Int_t b = 0;
+ Int_t e = fN - 1;
+ Int_t m = (b + e) / 2;
+
+ for ( ; b < e; m = (b + e) / 2) {
+ if (y > fClusters[m]->GetY()) {
+ b = m + 1;
+ }
+ else {
+ e = m;
+ }
+ }
+
+ return m;
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDpropagationLayer::FindNearestCluster(Float_t y, Float_t z
+ , Float_t maxroad
+ , Float_t maxroadz) const
+{
+ //
+ // Returns index of the cluster nearest to the given y,z
+ //
+
+ Int_t index = -1;
+ Int_t maxn = fN;
+ Float_t mindist = maxroad;
+
+ for (Int_t i = Find(y-maxroad); i < maxn; i++) {
+ AliTRDcluster *c = (AliTRDcluster *) (fClusters[i]);
+ Float_t ycl = c->GetY();
+ if (ycl > (y + maxroad)) {
+ break;
+ }
+ if (TMath::Abs(c->GetZ() - z) > maxroadz) {
+ continue;
+ }
+ if (TMath::Abs(ycl - y) < mindist) {
+ mindist = TMath::Abs(ycl - y);
+ index = fIndex[i];
+ }
+ }
+
+ return index;
+
+}
+
--- /dev/null
+#ifndef ALITRDPROPAGATIONLAYER_H
+#define ALITRDPROPAGATIONLAYER_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+
+////////////////////////////////////////////////////////////////////////////
+// //
+// The TRD propagation layer //
+// //
+////////////////////////////////////////////////////////////////////////////
+
+#ifndef ROOT_TObject
+#include "TObject.h"
+#endif
+
+class AliTRDcluster;
+
+class AliTRDpropagationLayer : public TObject
+{
+
+ public:
+
+ enum { kMaxClusterPerTimeBin = 2300
+ , kZones = 5 };
+
+ AliTRDpropagationLayer();
+ AliTRDpropagationLayer(Double_t x, Double_t dx, Double_t rho
+ , Double_t x0, Int_t tbIndex, Int_t plane);
+ AliTRDpropagationLayer(const AliTRDpropagationLayer &/*p*/);
+ ~AliTRDpropagationLayer();
+ AliTRDpropagationLayer &operator=(const AliTRDpropagationLayer &/*p*/) { return *this; }
+
+ operator Int_t() const { return fN; }
+ AliTRDcluster *operator[](Int_t i) { return fClusters[i]; }
+ virtual void Copy(TObject &o) const;
+
+ void SetZmax(Int_t cham, Double_t center, Double_t w) { fZc[cham] = center;
+ fZmax[cham] = w; }
+ void SetYmax(Double_t w, Double_t wsensitive) { fYmax = w;
+ fYmaxSensitive = wsensitive; }
+
+ void SetZ(Double_t* center, Double_t *w, Double_t *wsensitive);
+ void SetHoles(Bool_t* holes);
+ //void SetHole(Double_t Zmax, Double_t Ymax
+ // , Double_t rho = 1.29e-3, Double_t x0 = 36.66
+ // , Double_t Yc = 0.0, Double_t Zc = 0.0);
+ void SetSector(const Int_t sec) { fSec = sec; }
+ void SetX(Double_t x) { fX = x; }
+
+ Double_t GetYmax() const { return fYmax; }
+ Double_t GetZmax(Int_t c) const { return fZmax[c]; }
+ Double_t GetZc(Int_t c) const { return fZc[c]; }
+ UInt_t GetIndex(Int_t i) const { return fIndex[i]; }
+ Double_t GetX() const { return fX; }
+ Double_t GetdX() const { return fdX; }
+ Int_t GetTimeBinIndex() const { return fTimeBinIndex; }
+ Int_t GetPlane() const { return fPlane; }
+ Int_t GetSector() const { return fSec; }
+
+ Bool_t IsHole(Int_t zone) const { return fIsHole[zone]; }
+ Bool_t IsSensitive() const { return (fTimeBinIndex >= 0) ? kTRUE : kFALSE;}
+
+ void Clear(const Option_t * /*o*/) { ; }
+ void Clear() { for (Int_t i = 0; i < fN; i++)
+ fClusters[i] = NULL;
+ fN = 0; }
+
+ void InsertCluster(AliTRDcluster *c, UInt_t index);
+ Int_t Find(Float_t y) const;
+ Int_t FindNearestCluster(Float_t y, Float_t z, Float_t maxroad, Float_t maxroadz) const;
+
+ protected:
+
+ Int_t fN; // Number of clusters
+ Int_t fSec; // Sector mumber
+ AliTRDcluster **fClusters; //[fN] Array of pointers to clusters
+ UInt_t *fIndex; //[fN] Array of cluster indexes
+ Double_t fX; // X coordinate of the middle plane
+ Double_t fdX; // Radial thickness of the time bin
+ Double_t fRho; // Default density of the material
+ Double_t fX0; // Default radiation length
+ Int_t fTimeBinIndex; // Plane * F(local_tb)
+ Int_t fPlane; // Plane number
+ Double_t fZc[kZones]; // Z position of the center for 5 active areas
+ Double_t fZmax[kZones]; // Half of active area length in Z
+ Double_t fZmaxSensitive[kZones]; // Sensitive area for detection Z
+ Bool_t fIsHole[kZones]; // Is hole in given sector
+ Double_t fYmax; // Half of active area length in Y
+ Double_t fYmaxSensitive; // Half of active area length in Y
+ Bool_t fHole; // kTRUE if there is a hole in the layer
+ Double_t fHoleZc; // Z of the center of the hole
+ Double_t fHoleZmax; // Half of the hole length in Z
+ Double_t fHoleYc; // Y of the center of the hole
+ Double_t fHoleYmax; // Half of the hole length in Y
+ Double_t fHoleRho; // Density of the gas in the hole
+ Double_t fHoleX0; // Radiation length of the gas in the hole
+
+ ClassDef(AliTRDpropagationLayer, 1) // Propagation layer for TRD tracker
+
+};
+
+#endif
+
--- /dev/null
+/**************************************************************************
+ * 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$ */
+
+///////////////////////////////////////////////////////////////////////////////
+// //
+// Parameter class for the TRD reconstruction //
+// //
+// Authors: //
+// Alex Bercuci <A.Bercuci@gsi.de> //
+// Markus Fasel <M.Fasel@gsi.de> //
+// //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "AliTRDrecoParam.h"
+
+ClassImp(AliTRDrecoParam)
+
+//______________________________________________________________
+AliTRDrecoParam::AliTRDrecoParam()
+ :AliDetectorRecoParam()
+ ,fkMaxTheta(1.0)
+ ,fkMaxPhi(2.0)
+ ,fkRoad0y(6.0)
+ ,fkRoad0z(8.5)
+ ,fkRoad1y(2.0)
+ ,fkRoad1z(20.0)
+ ,fkRoad2y(3.0)
+ ,fkRoad2z(20.0)
+ ,fkPlaneQualityThreshold(5.0)// 4.2? under Investigation
+ ,fkFindable(.333)
+ ,fkChi2Z(14./*12.5*/)
+ ,fkChi2Y(.25)
+ ,fkTrackLikelihood(-15.)
+{
+ //
+ // Default constructor
+ //
+
+}
+
+//______________________________________________________________
+AliTRDrecoParam *AliTRDrecoParam::GetLowFluxParam()
+{
+ //
+ // Parameters for the low flux environment
+ //
+
+ return new AliTRDrecoParam();
+
+}
+
+//______________________________________________________________
+AliTRDrecoParam *AliTRDrecoParam::GetHighFluxParam()
+{
+ //
+ // Parameters for the high flux environment
+ //
+
+ return new AliTRDrecoParam();
+
+}
--- /dev/null
+#ifndef ALITRDRECOPARAM_H
+#define ALITRDRECOPARAM_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+
+////////////////////////////////////////////////////////////////////////////
+// //
+// Parameter class for the TRD reconstruction //
+// //
+////////////////////////////////////////////////////////////////////////////
+
+#ifndef ALIDETECTORRECOPARAM_H
+#include "AliDetectorRecoParam.h"
+#endif
+
+class AliTRDrecoParam : public AliDetectorRecoParam
+{
+
+ public:
+
+ AliTRDrecoParam();
+ ~AliTRDrecoParam() { }
+
+ Double_t GetChi2Y() const { return fkChi2Y; }
+ Double_t GetChi2Z() const { return fkChi2Z; }
+ Double_t GetFindableClusters() const { return fkFindable; }
+ Double_t GetMaxTheta() const { return fkMaxTheta; }
+ Double_t GetMaxPhi() const { return fkMaxPhi; }
+
+ Double_t GetRoad0y() const { return fkRoad0y; }
+ Double_t GetRoad0z() const { return fkRoad0z; }
+
+ Double_t GetRoad1y() const { return fkRoad1y; }
+ Double_t GetRoad1z() const { return fkRoad1z; }
+
+ Double_t GetRoad2y() const { return fkRoad2y; }
+ Double_t GetRoad2z() const { return fkRoad2z; }
+
+ Double_t GetPlaneQualityThreshold() const { return fkPlaneQualityThreshold; }
+
+ Double_t GetTrackLikelihood() const { return fkTrackLikelihood; }
+
+ static AliTRDrecoParam *GetLowFluxParam();
+ static AliTRDrecoParam *GetHighFluxParam();
+
+ private:
+
+ Double_t fkMaxTheta; // Maximum theta
+ Double_t fkMaxPhi; // Maximum phi
+
+ Double_t fkRoad0y; // Road for middle cluster
+ Double_t fkRoad0z; // Road for middle cluster
+
+ Double_t fkRoad1y; // Road in y for seeded cluster
+ Double_t fkRoad1z; // Road in z for seeded cluster
+
+ Double_t fkRoad2y; // Road in y for extrapolated cluster
+ Double_t fkRoad2z; // Road in z for extrapolated cluster
+
+ Double_t fkPlaneQualityThreshold; // Quality threshold
+ Double_t fkFindable; // Ratio of clusters from a track in one chamber which are at minimum supposed to be found.
+ 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
+
+ ClassDef(AliTRDrecoParam, 1) // Reconstruction parameters for TRD detector
+
+};
+#endif
#include "AliMathBase.h"
#include "AliTRDseed.h"
+#include "AliTRDcalibDB.h"
#include "AliTRDcluster.h"
#include "AliTRDtracker.h"
// Default constructor
//
- for (Int_t i = 0; i < 25; i++) {
+ for (Int_t i = 0; i < knTimebins; i++) {
fX[i] = 0; // x position
fY[i] = 0; // y position
fZ[i] = 0; // z position
// Copy constructor
//
- for (Int_t i = 0; i < 25; i++) {
+ for (Int_t i = 0; i < knTimebins; i++) {
fX[i] = s.fX[i]; // x position
fY[i] = s.fY[i]; // y position
fZ[i] = s.fZ[i]; // z position
}
+//_____________________________________________________________________________
+void AliTRDseed::Copy(TObject &o) const {
+ AliTRDseed &seed = (AliTRDseed &)o;
+ seed.fTilt = fTilt;
+ seed.fPadLength = fPadLength;
+ seed.fX0 = fX0;
+ seed.fSigmaY = fSigmaY;
+ seed.fSigmaY2 = fSigmaY2;
+ seed.fMeanz = fMeanz;
+ seed.fZProb = fZProb;
+ seed.fN = fN;
+ seed.fN2 = fN2;
+ seed.fNUsed = fNUsed;
+ seed.fFreq = fFreq;
+ seed.fNChange = fNChange;
+ seed.fMPads = fMPads;
+ seed.fC = fC;
+ seed.fCC = fCC;
+ seed.fChi2 = fChi2;
+ seed.fChi2Z = fChi2Z;
+ for (Int_t i = 0; i < knTimebins; i++) {
+ seed.fX[i] = fX[i];
+ seed.fY[i] = fY[i];
+ seed.fZ[i] = fZ[i];
+ seed.fIndexes[i] = fIndexes[i];
+ seed.fClusters[i] = fClusters[i];
+ seed.fUsable[i] = fUsable[i];
+ }
+
+ for (Int_t i = 0; i < 2; i++) {
+ seed.fYref[i] = fYref[i];
+ seed.fZref[i] = fZref[i];
+ seed.fYfit[i] = fYfit[i];
+ seed.fYfitR[i] = fYfitR[i];
+ seed.fZfit[i] = fZfit[i];
+ seed.fZfitR[i] = fZfitR[i];
+ seed.fLabels[i] = fLabels[i];
+ }
+
+ TObject::Copy(seed);
+
+}
+
//_____________________________________________________________________________
void AliTRDseed::Reset()
{
// Reset seed
//
- for (Int_t i = 0; i < 25; i++) {
+ for (Int_t i = 0; i < knTimebins; i++) {
fX[i] = 0; // X position
fY[i] = 0; // Y position
fZ[i] = 0; // Z position
fIndexes[i] = 0; // Indexes
- fClusters[i] = 0; // Clusters
+ fClusters[i] = 0x0; // Clusters
fUsable[i] = kFALSE;
}
// Cook 2 labels for seed
//
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ Int_t nTimeBins = cal->GetNumberOfTimeBins();
+
Int_t labels[200];
Int_t out[200];
Int_t nlab = 0;
- for (Int_t i = 0; i < 25; i++) {
+ for (Int_t i = 0; i < nTimeBins+1; i++) {
if (!fClusters[i]) continue;
for (Int_t ilab = 0; ilab < 3; ilab++) {
if (fClusters[i]->GetLabel(ilab) >= 0) {
// Use clusters
//
- for (Int_t i = 0; i < 25; i++) {
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ Int_t nTimeBins = cal->GetNumberOfTimeBins();
+
+ for (Int_t i = 0; i < nTimeBins+1; i++) {
if (!fClusters[i]) continue;
if (!(fClusters[i]->IsUsed())) fClusters[i]->Use();
}
//
const Float_t kRatio = 0.8;
- const Int_t kClmin = 6;
+ const Int_t kClmin = 5;
const Float_t kmaxtan = 2;
- if (TMath::Abs(fYref[1]) > kmaxtan) return; // Track inclined too much
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ Int_t nTimeBins = cal->GetNumberOfTimeBins();
+
+ if (TMath::Abs(fYref[1]) > kmaxtan){
+ //printf("Exit: Abs(fYref[1]) = %3.3f, kmaxtan = %3.3f\n", TMath::Abs(fYref[1]), kmaxtan);
+ return; // Track inclined too much
+ }
Float_t sigmaexp = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction
Float_t ycrosscor = fPadLength * fTilt * 0.5; // Y correction for crossing
Double_t sumwz;
Double_t sumwxz;
- Int_t zints[25]; // Histograming of the z coordinate
+ // Buffering: Leave it constant fot Performance issues
+ Int_t zints[knTimebins]; // Histograming of the z coordinate
// Get 1 and second max probable coodinates in z
- Int_t zouts[50];
- Float_t allowedz[25]; // Allowed z for given time bin
- Float_t yres[25]; // Residuals from reference
+ Int_t zouts[2*knTimebins];
+ Float_t allowedz[knTimebins]; // Allowed z for given time bin
+ Float_t yres[knTimebins]; // Residuals from reference
Float_t anglecor = fTilt * fZref[1]; // Correction to the angle
fN = 0;
fN2 = 0;
- for (Int_t i = 0; i < 25; i++) {
+ for (Int_t i = 0; i < nTimeBins; i++) {
yres[i] = 10000.0;
if (!fClusters[i]) continue;
yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i]; // Residual y
fN++;
}
- if (fN < kClmin) return;
+ if (fN < kClmin){
+ //printf("Exit fN < kClmin: fN = %d\n", fN);
+ return;
+ }
Int_t nz = AliTRDtracker::Freq(fN,zints,zouts,kFALSE);
fZProb = zouts[0];
if (nz <= 1) zouts[3] = 0;
- if (zouts[1] + zouts[3] < kClmin) return;
+ if (zouts[1] + zouts[3] < kClmin) {
+ //printf("Exit zouts[1] = %d, zouts[3] = %d\n",zouts[1],zouts[3]);
+ return;
+ }
// Z distance bigger than pad - length
if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) {
Int_t breaktime = -1;
Bool_t mbefore = kFALSE;
- Int_t cumul[25][2];
+ Int_t cumul[knTimebins][2];
Int_t counts[2] = { 0, 0 };
if (zouts[3] >= 3) {
// with maximal numebr of accepted clusters
//
fNChange = 1;
- for (Int_t i = 0; i < 25; i++) {
+ for (Int_t i = 0; i < nTimeBins; i++) {
cumul[i][0] = counts[0];
cumul[i][1] = counts[1];
if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++;
if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++;
}
Int_t maxcount = 0;
- for (Int_t i = 0; i < 24; i++) {
- Int_t after = cumul[24][0] - cumul[i][0];
+ for (Int_t i = 0; i < nTimeBins; i++) {
+ Int_t after = cumul[nTimeBins][0] - cumul[i][0];
Int_t before = cumul[i][1];
if (after + before > maxcount) {
maxcount = after + before;
breaktime = i;
mbefore = kFALSE;
}
- after = cumul[24][1] - cumul[i][1];
+ after = cumul[nTimeBins-1][1] - cumul[i][1];
before = cumul[i][0];
if (after + before > maxcount) {
maxcount = after + before;
}
- for (Int_t i = 0; i < 25; i++) {
+ for (Int_t i = 0; i < nTimeBins+1; i++) {
if (i > breaktime) allowedz[i] = mbefore ? zouts[2] : zouts[0];
if (i <= breaktime) allowedz[i] = (!mbefore) ? zouts[2] : zouts[0];
}
- if (((allowedz[0] > allowedz[24]) && (fZref[1] < 0)) ||
- ((allowedz[0] < allowedz[24]) && (fZref[1] > 0))) {
+ if (((allowedz[0] > allowedz[nTimeBins]) && (fZref[1] < 0)) ||
+ ((allowedz[0] < allowedz[nTimeBins]) && (fZref[1] > 0))) {
//
// Tracklet z-direction not in correspondance with track z direction
//
fNChange = 0;
- for (Int_t i = 0; i < 25; i++) {
+ for (Int_t i = 0; i < nTimeBins+1; i++) {
allowedz[i] = zouts[0]; // Only longest taken
}
}
//
// Cross pad -row tracklet - take the step change into account
//
- for (Int_t i = 0; i < 25; i++) {
+ for (Int_t i = 0; i < nTimeBins+1; i++) {
if (!fClusters[i]) continue;
if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i]; // Residual y
}
}
- Double_t yres2[25];
+ Double_t yres2[knTimebins];
Double_t mean;
Double_t sigma;
- for (Int_t i = 0; i < 25; i++) {
+ for (Int_t i = 0; i < nTimeBins+1; i++) {
if (!fClusters[i]) continue;
if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
yres2[fN2] = yres[i];
fN2++;
}
if (fN2 < kClmin) {
+ //printf("Exit fN2 < kClmin: fN2 = %d\n", fN2);
fN2 = 0;
return;
}
- AliMathBase::EvaluateUni(fN2,yres2,mean,sigma,Int_t(fN2*kRatio-2));
+ AliMathBase::EvaluateUni(fN2,yres2,mean,sigma, Int_t(fN2*kRatio-2.));
if (sigma < sigmaexp * 0.8) {
sigma = sigmaexp;
}
fMeanz = 0;
fMPads = 0;
- for (Int_t i = 0; i < 25; i++) {
+ for (Int_t i = 0; i < nTimeBins+1; i++) {
fUsable[i] = kFALSE;
if (!fClusters[i]) continue;
}
if (fN2 < kClmin){
+ //printf("Exit fN2 < kClmin(2): fN2 = %d\n",fN2);
fN2 = 0;
return;
}
fYfitR[1] = (sumw * sumwxy - sumwx * sumwy) / det;
fSigmaY2 = 0;
- for (Int_t i = 0; i < 25; i++) {
+ for (Int_t i = 0; i < nTimeBins+1; i++) {
if (!fUsable[i]) continue;
Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i];
fSigmaY2 += delta*delta;
// Update used seed
//
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ Int_t nTimeBins = cal->GetNumberOfTimeBins();
+
fNUsed = 0;
- for (Int_t i = 0; i < 25; i++) {
+ for (Int_t i = 0; i < nTimeBins; i++) {
if (!fClusters[i]) {
continue;
}
// Fitting with tilting pads - kz not fixed
TLinearFitter fitterT2(4,"hyp4");
fitterT2.StoreData(kTRUE);
+
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ Int_t nTimeBins = cal->GetNumberOfTimeBins();
+
Float_t xref2 = (cseed[2].fX0 + cseed[3].fX0) * 0.5; // Reference x0 for z
Int_t npointsT = 0;
if (!cseed[iLayer].IsOK()) continue;
Double_t tilt = cseed[iLayer].fTilt;
- for (Int_t itime = 0; itime < 25; itime++) {
+ for (Int_t itime = 0; itime < nTimeBins+1; itime++) {
if (!cseed[iLayer].fUsable[itime]) continue;
// x relative to the midle chamber
//
Bool_t acceptablez = kTRUE;
for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
- if (cseed[iLayer].IsOK()) {
- Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].fX0 - xref2);
- if (TMath::Abs(cseed[iLayer].fZProb - zT2) > cseed[iLayer].fPadLength * 0.5 + 1.0) {
- acceptablez = kFALSE;
- }
- }
+ if (!cseed[iLayer].IsOK()) continue;
+ Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].fX0 - xref2);
+ if (TMath::Abs(cseed[iLayer].fZProb - zT2) > cseed[iLayer].fPadLength * 0.5 + 1.0) acceptablez = kFALSE;
}
if (!acceptablez) {
Double_t zmf = cseed[2].fZref[0] + cseed[2].fZref[1] * (xref2 - cseed[2].fX0);
params[2] = fitterT2.GetParameter(2);
Double_t curvature = 1.0 + params[1] * params[1] - params[2] * params[0];
+
for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
-
Double_t x = cseed[iLayer].fX0;
Double_t y = 0;
Double_t dy = 0;
if (-params[2]*params[0] + params[1]*params[1] + 1 > 0) {
Double_t rm1 = params[0] / TMath::Sqrt(-params[2]*params[0] + params[1]*params[1] + 1);
if (1.0/(rm1*rm1) - (x-x0) * (x-x0) > 0.0) {
- Double_t res = (x - x0) / TMath::Sqrt(1.0 / (rm1*rm1) - (x-x0)*(x-x0));
- if (params[0] < 0) res *= -1.0;
- dy = res;
+ Double_t res = (x - x0) / TMath::Sqrt(1.0 / (rm1*rm1) - (x-x0)*(x-x0));
+ if (params[0] < 0) res *= -1.0;
+ dy = res;
}
}
z = rpolz0 + rpolz1 * (x - xref2);
class AliTRDseed : public TObject {
public:
+ enum { knTimebins = 35 };
AliTRDseed();
AliTRDseed(const AliTRDseed &s);
AliTRDseed &operator=(const AliTRDseed &s) { *(new(this) AliTRDseed(s));
return *this; }
-
static Float_t FitRiemanTilt(AliTRDseed *seed, Bool_t error);
void UseClusters();
void Update();
void UpdateUsed();
void Reset();
- Bool_t IsOK() const { return fN2 > 8; }
+ Bool_t IsOK() const { return fN2 > 4; }
Bool_t IsUsable(Int_t i) const { return fUsable[i]; }
Float_t GetTilt() const { return fTilt; }
void SetChi2(Float_t chi2) { fChi2 = chi2; }
void SetChi2Z(Float_t chi2z) { fChi2Z = chi2z; }
- private:
+ protected:
+ void Copy(TObject &o) const;
+
Float_t fTilt; // Tilting angle
Float_t fPadLength; // Pad length
Float_t fX0; // X0 position
- Float_t fX[25]; //! X position
- Float_t fY[25]; //! Y position
- Float_t fZ[25]; //! Z position
- Int_t fIndexes[25]; //! Indexes
- AliTRDcluster *fClusters[25]; //! Clusters
- Bool_t fUsable[25]; //! Indication - usable cluster
+ Float_t fX[knTimebins]; //! X position
+ Float_t fY[knTimebins]; //! Y position
+ Float_t fZ[knTimebins]; //! Z position
+ Int_t fIndexes[knTimebins];//! Indexes
+ AliTRDcluster *fClusters[knTimebins]; // Clusters
+ Bool_t fUsable[knTimebins]; //! Indication - usable cluster
Float_t fYref[2]; // Reference y
Float_t fZref[2]; // Reference z
Float_t fYfit[2]; // Y fit position +derivation
--- /dev/null
+/**************************************************************************
+ * 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$ */
+
+////////////////////////////////////////////////////////////////////////////
+// //
+// The TRD track seed //
+// //
+// Authors: //
+// Alex Bercuci <A.Bercuci@gsi.de> //
+// Markus Fasel <M.Fasel@gsi.de> //
+// //
+////////////////////////////////////////////////////////////////////////////
+
+#include "TMath.h"
+#include "TLinearFitter.h"
+
+#include "AliLog.h"
+#include "AliMathBase.h"
+
+#include "AliTRDseedV1.h"
+#include "AliTRDcluster.h"
+#include "AliTRDcalibDB.h"
+#include "AliTRDstackLayer.h"
+#include "AliTRDrecoParam.h"
+
+#define SEED_DEBUG
+
+ClassImp(AliTRDseedV1)
+
+//____________________________________________________________________
+AliTRDseedV1::AliTRDseedV1(Int_t layer, AliTRDrecoParam *p)
+ :AliTRDseed()
+ ,fLayer(layer)
+ ,fTimeBins(0)
+ ,fOwner(kFALSE)
+ ,fRecoParam(p)
+{
+ //
+ // Constructor
+ //
+
+ //AliInfo("");
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ fTimeBins = cal->GetNumberOfTimeBins();
+
+}
+
+//____________________________________________________________________
+AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref, Bool_t owner)
+ :AliTRDseed((AliTRDseed&)ref)
+ ,fLayer(ref.fLayer)
+ ,fTimeBins(ref.fTimeBins)
+ ,fOwner(kFALSE)
+ ,fRecoParam(ref.fRecoParam)
+{
+ //
+ // Copy Constructor performing a deep copy
+ //
+
+ //AliInfo("");
+
+ if(owner){
+ for(int ic=0; ic<fTimeBins; ic++){
+ if(!fClusters[ic]) continue;
+ fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
+ }
+ fOwner = kTRUE;
+ }
+
+}
+
+//____________________________________________________________________
+AliTRDseedV1& AliTRDseedV1::operator=(const AliTRDseedV1 &ref)
+{
+ //
+ // Assignment Operator using the copy function
+ //
+
+ //AliInfo("");
+ if(this != &ref){
+ ref.Copy(*this);
+ }
+ return *this;
+
+}
+
+//____________________________________________________________________
+AliTRDseedV1::~AliTRDseedV1()
+{
+ //
+ // Destructor. The RecoParam object belongs to the underlying tracker.
+ //
+
+ //AliInfo(Form("fOwner[%s]", fOwner?"YES":"NO"));
+
+ if(fOwner) delete [] fClusters;
+}
+
+//____________________________________________________________________
+void AliTRDseedV1::Copy(TObject &ref) const
+{
+ //
+ // Copy function
+ //
+
+ //AliInfo("");
+ AliTRDseedV1 &target = (AliTRDseedV1 &)ref;
+
+ target.fLayer = fLayer;
+ target.fTimeBins = fTimeBins;
+ target.fRecoParam = fRecoParam;
+ AliTRDseed::Copy(target);
+}
+
+//____________________________________________________________________
+Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
+{
+ //
+ // Returns a quality measurement of the current seed
+ //
+
+ Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
+ return .5 * (18.0 - fN2)
+ + 10.* TMath::Abs(fYfit[1] - fYref[1])
+ + 5.* TMath::Abs(fYfit[0] - fYref[0] + zcorr)
+ + 2. * TMath::Abs(fMeanz - fZref[0]) / fPadLength;
+}
+
+//____________________________________________________________________
+Bool_t AliTRDseedV1::AttachClustersIter(AliTRDstackLayer *layer
+ , Float_t quality
+ , Bool_t kZcorr
+ , AliTRDcluster *c)
+{
+ //
+ // Iterative process to register clusters to the seed.
+ // In iteration 0 we try only one pad-row and if quality not
+ // sufficient we try 2 pad-rows (about 5% of tracks cross 2 pad-rows)
+ //
+
+ if(!fRecoParam){
+ AliError("Seed can not be used without a valid RecoParam.");
+ return kFALSE;
+ }
+
+ Float_t tquality;
+ Double_t kroady = fRecoParam->GetRoad1y();
+ Double_t kroadz = fPadLength * .5 + 1.;
+
+ // initialize configuration parameters
+ Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
+ Int_t niter = kZcorr ? 1 : 2;
+
+ Double_t yexp, zexp;
+ Int_t ncl = 0;
+ // start seed update
+ for (Int_t iter = 0; iter < niter; iter++) {
+ //AliInfo(Form("iter = %i", iter));
+ ncl = 0;
+ for (Int_t iTime = 0; iTime < fTimeBins; iTime++) {
+ // define searching configuration
+ Double_t dxlayer = layer[iTime].GetX() - fX0;
+ if(c){
+ zexp = c->GetZ();
+ //Try 2 pad-rows in second iteration
+ if (iter > 0) {
+ zexp = fZref[0] + fZref[1] * dxlayer - zcorr;
+ if (zexp > c->GetZ()) zexp = c->GetZ() + fPadLength*0.5;
+ if (zexp < c->GetZ()) zexp = c->GetZ() - fPadLength*0.5;
+ }
+ } else zexp = fZref[0];
+ yexp = fYref[0] + fYref[1] * dxlayer - zcorr;
+ // get cluster
+// printf("xexp = %3.3f ,yexp = %3.3f, zexp = %3.3f\n",layer[iTime].GetX(),yexp,zexp);
+// printf("layer[%i].GetNClusters() = %i\n", iTime, layer[iTime].GetNClusters());
+ Int_t index = layer[iTime].SearchNearestCluster(yexp, zexp, kroady, kroadz);
+// for(Int_t iclk = 0; iclk < layer[iTime].GetNClusters(); iclk++){
+// AliTRDcluster *testcl = layer[iTime].GetCluster(iclk);
+// printf("Cluster %i: x = %3.3f, y = %3.3f, z = %3.3f\n",iclk,testcl->GetX(), testcl->GetY(), testcl->GetZ());
+// }
+// printf("Index = %i\n",index);
+ if (index < 0) continue;
+
+ // Register cluster
+ AliTRDcluster *cl = (AliTRDcluster*) layer[iTime].GetCluster(index);
+
+ //printf("Cluster %i(0x%x): x = %3.3f, y = %3.3f, z = %3.3f\n", index, cl, cl->GetX(), cl->GetY(), cl->GetZ());
+
+ Int_t GlobalIndex = layer[iTime].GetGlobalIndex(index);
+ fIndexes[iTime] = GlobalIndex;
+ fClusters[iTime] = cl;
+ fX[iTime] = dxlayer;
+ fY[iTime] = cl->GetY();
+ fZ[iTime] = cl->GetZ();
+
+ // Debugging
+ ncl++;
+ }
+
+#ifdef SEED_DEBUG
+// Int_t nclusters = 0;
+// Float_t fD[iter] = 0.;
+// for(int ic=0; ic<fTimeBins+1; ic++){
+// AliTRDcluster *ci = fClusters[ic];
+// if(!ci) continue;
+// for(int jc=ic+1; jc<fTimeBins+1; jc++){
+// AliTRDcluster *cj = fClusters[jc];
+// if(!cj) continue;
+// fD[iter] += TMath::Sqrt((ci->GetY()-cj->GetY())*(ci->GetY()-cj->GetY())+
+// (ci->GetZ()-cj->GetZ())*(ci->GetZ()-cj->GetZ()));
+// nclusters++;
+// }
+// }
+// if(nclusters) fD[iter] /= float(nclusters);
+#endif
+
+ AliTRDseed::Update();
+
+ if(IsOK()){
+ tquality = GetQuality(kZcorr);
+ if(tquality < quality) break;
+ else quality = tquality;
+ }
+ kroadz *= 2.;
+ } // Loop: iter
+ if (!IsOK()) return kFALSE;
+
+ CookLabels();
+ UpdateUsed();
+ return kTRUE;
+}
+
+//____________________________________________________________________
+Bool_t AliTRDseedV1::AttachClustersProj(AliTRDstackLayer *layer
+ , Float_t /*quality*/
+ , Bool_t kZcorr
+ , AliTRDcluster *c)
+{
+ //
+ // Projective algorithm to attach clusters to seeding tracklets
+ //
+ // Parameters
+ //
+ // Output
+ //
+ // Detailed description
+ // 1. Collapse x coordinate for the full detector plane
+ // 2. truncated mean on y (r-phi) direction
+ // 3. purge clusters
+ // 4. truncated mean on z direction
+ // 5. purge clusters
+ // 6. fit tracklet
+ //
+
+ if(!fRecoParam){
+ AliError("Seed can not be used without a valid RecoParam.");
+ return kFALSE;
+ }
+
+ const Int_t knTimeBins = 35;
+ const Int_t kClusterCandidates = 2 * knTimeBins;
+
+ //define roads
+ Double_t kroady = fRecoParam->GetRoad1y();
+ Double_t kroadz = fPadLength * 1.5 + 1.;
+ // correction to y for the tilting angle
+ Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
+
+ // working variables
+ AliTRDcluster *clusters[kClusterCandidates];
+ Double_t cond[4], yexp[knTimeBins], zexp[knTimeBins],
+ yres[kClusterCandidates], zres[kClusterCandidates];
+ Int_t ncl, *index = 0x0, tboundary[knTimeBins];
+
+ // Do cluster projection
+ Int_t nYclusters = 0; Bool_t kEXIT = kFALSE;
+ for (Int_t iTime = 0; iTime < fTimeBins; iTime++) {
+ fX[iTime] = layer[iTime].GetX() - fX0;
+ zexp[iTime] = fZref[0] + fZref[1] * fX[iTime];
+ yexp[iTime] = fYref[0] + fYref[1] * fX[iTime] - zcorr;
+
+ // build condition and process clusters
+ cond[0] = yexp[iTime] - kroady; cond[1] = yexp[iTime] + kroady;
+ cond[2] = zexp[iTime] - kroadz; cond[3] = zexp[iTime] + kroadz;
+ layer[iTime].GetClusters(cond, index, ncl);
+ for(Int_t ic = 0; ic<ncl; ic++){
+ c = layer[iTime].GetCluster(index[ic]);
+ clusters[nYclusters] = c;
+ yres[nYclusters++] = c->GetY() - yexp[iTime];
+ if(nYclusters >= kClusterCandidates) {
+ AliWarning(Form("Cluster candidates reached limit %d. Some may be lost.", kClusterCandidates));
+ kEXIT = kTRUE;
+ break;
+ }
+ }
+ tboundary[iTime] = nYclusters;
+ if(kEXIT) break;
+ }
+
+ // Evaluate truncated mean on the y direction
+ Double_t mean, sigma;
+ AliMathBase::EvaluateUni(nYclusters, yres, mean, sigma, Int_t(nYclusters*.8)-2);
+ //purge cluster candidates
+ Int_t nZclusters = 0;
+ for(Int_t ic = 0; ic<nYclusters; ic++){
+ if(yres[ic] - mean > 4. * sigma){
+ clusters[ic] = 0x0;
+ continue;
+ }
+ zres[nZclusters++] = clusters[ic]->GetZ() - zexp[clusters[ic]->GetLocalTimeBin()];
+ }
+
+ // Evaluate truncated mean on the z direction
+ AliMathBase::EvaluateUni(nZclusters, zres, mean, sigma, Int_t(nZclusters*.8)-2);
+ //purge cluster candidates
+ for(Int_t ic = 0; ic<nZclusters; ic++){
+ if(zres[ic] - mean > 4. * sigma){
+ clusters[ic] = 0x0;
+ continue;
+ }
+ }
+
+
+ // Select only one cluster/TimeBin
+ Int_t lastCluster = 0;
+ fN2 = 0;
+ for (Int_t iTime = 0; iTime < fTimeBins; iTime++) {
+ ncl = tboundary[iTime] - lastCluster;
+ if(!ncl) continue;
+ if(ncl == 1){
+ c = clusters[lastCluster];
+ } else if(ncl > 1){
+ Float_t dold = 9999.; Int_t iptr = lastCluster;
+ for(int ic=lastCluster; ic<tboundary[iTime]; ic++){
+ if(!clusters[ic]) continue;
+ Float_t y = yexp[iTime] - clusters[ic]->GetY();
+ Float_t z = zexp[iTime] - clusters[ic]->GetZ();
+ Float_t d = y * y + z * z;
+ if(d > dold) continue;
+ dold = d;
+ iptr = ic;
+ }
+ c = clusters[iptr];
+ }
+ //Int_t GlobalIndex = layer[iTime].GetGlobalIndex(index);
+ //fIndexes[iTime] = GlobalIndex;
+ fClusters[iTime] = c;
+ fY[iTime] = c->GetY();
+ fZ[iTime] = c->GetZ();
+ lastCluster = tboundary[iTime];
+ fN2++;
+ }
+
+ // number of minimum numbers of clusters expected for the tracklet
+ Int_t kClmin = Int_t(fRecoParam->GetFindableClusters()*fTimeBins);
+ if (fN2 < kClmin){
+ AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
+ fN2 = 0;
+ return kFALSE;
+ }
+ AliTRDseed::Update();
+
+// // fit tracklet and update clusters
+// if(!FitTracklet()) return kFALSE;
+// UpdateUsed();
+ return kTRUE;
+}
+
+//____________________________________________________________________
+Bool_t AliTRDseedV1::FitTracklet()
+{
+ //
+ // Linear fit of the tracklet
+ //
+ // Parameters :
+ //
+ // Output :
+ // True if successful
+ //
+ // Detailed description
+ // 2. Check if tracklet crosses pad row boundary
+ // 1. Calculate residuals in the y (r-phi) direction
+ // 3. Do a Least Square Fit to the data
+ //
+
+ //Float_t sigmaexp = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction
+ Float_t ycrosscor = fPadLength * fTilt * 0.5; // Y correction for crossing
+ Float_t anglecor = fTilt * fZref[1]; // Correction to the angle
+
+ // calculate residuals
+ const Int_t knTimeBins = 35;
+ Float_t yres[knTimeBins]; // y (r-phi) residuals
+ Int_t zint[knTimeBins], // Histograming of the z coordinate
+ zout[2*knTimeBins];//
+
+ fN = 0;
+ for (Int_t iTime = 0; iTime < fTimeBins; iTime++) {
+ if (!fClusters[iTime]) continue;
+ yres[iTime] = fY[iTime] - fYref[0] - (fYref[1] + anglecor) * fX[iTime];
+ zint[fN++] = Int_t(fZ[iTime]);
+ }
+
+ // calculate pad row boundary crosses
+ Int_t kClmin = Int_t(fRecoParam->GetFindableClusters()*fTimeBins);
+ Int_t nz = AliMathBase::Freq(fN, zint, zout, kFALSE);
+ fZProb = zout[0];
+ if(nz <= 1) zout[3] = 0;
+ if(zout[1] + zout[3] < kClmin) {
+ AliWarning(Form("Not enough clusters to fit the cross boundary tracklet %d [%d].", zout[1]+zout[3], kClmin));
+ return kFALSE;
+ }
+ // Z distance bigger than pad - length
+ if (TMath::Abs(zout[0]-zout[2]) > fPadLength) zout[3]=0;
+
+
+ Double_t sumw = 0.,
+ sumwx = 0.,
+ sumwx2 = 0.,
+ sumwy = 0.,
+ sumwxy = 0.,
+ sumwz = 0.,
+ sumwxz = 0.;
+ Int_t npads;
+ fMPads = 0;
+ fMeanz = 0.;
+ for(int iTime=0; iTime<fTimeBins; iTime++){
+ fUsable[iTime] = kFALSE;
+ if (!fClusters[iTime]) continue;
+ npads = fClusters[iTime]->GetNPads();
+
+ fUsable[iTime] = kTRUE;
+ fN2++;
+ fMPads += npads;
+ Float_t weight = 1.0;
+ if(npads > 5) weight = 0.2;
+ else if(npads > 4) weight = 0.5;
+ sumw += weight;
+ sumwx += fX[iTime] * weight;
+ sumwx2 += fX[iTime] * fX[iTime] * weight;
+ sumwy += weight * yres[iTime];
+ sumwxy += weight * yres[iTime] * fX[iTime];
+ sumwz += weight * fZ[iTime];
+ sumwxz += weight * fZ[iTime] * fX[iTime];
+ }
+ if (fN2 < kClmin){
+ AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
+ fN2 = 0;
+ return kFALSE;
+ }
+ fMeanz = sumwz / sumw;
+ fNChange = 0;
+
+ // Tracklet on boundary
+ Float_t correction = 0;
+ if (fNChange > 0) {
+ if (fMeanz < fZProb) correction = ycrosscor;
+ if (fMeanz > fZProb) correction = -ycrosscor;
+ }
+
+ Double_t det = sumw * sumwx2 - sumwx * sumwx;
+ fYfitR[0] = (sumwx2 * sumwy - sumwx * sumwxy) / det;
+ fYfitR[1] = (sumw * sumwxy - sumwx * sumwy) / det;
+
+ fSigmaY2 = 0;
+ for (Int_t i = 0; i < fTimeBins+1; i++) {
+ if (!fUsable[i]) continue;
+ Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i];
+ fSigmaY2 += delta*delta;
+ }
+ fSigmaY2 = TMath::Sqrt(fSigmaY2 / Float_t(fN2-2));
+
+ fZfitR[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det;
+ fZfitR[1] = (sumw * sumwxz - sumwx * sumwz) / det;
+ fZfit[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det;
+ fZfit[1] = (sumw * sumwxz - sumwx * sumwz) / det;
+ fYfitR[0] += fYref[0] + correction;
+ fYfitR[1] += fYref[1];
+ fYfit[0] = fYfitR[0];
+ fYfit[1] = fYfitR[1];
+
+ return kTRUE;
+}
+
+//_____________________________________________________________________________
+Float_t AliTRDseedV1::FitRiemanTilt(AliTRDseedV1 *cseed, Bool_t terror)
+{
+ //
+ // Fit the Rieman tilt
+ //
+
+ // Fitting with tilting pads - kz not fixed
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ Int_t nTimeBins = cal->GetNumberOfTimeBins();
+ TLinearFitter fitterT2(4,"hyp4");
+ fitterT2.StoreData(kTRUE);
+ Float_t xref2 = (cseed[2].fX0 + cseed[3].fX0) * 0.5; // Reference x0 for z
+
+ Int_t npointsT = 0;
+ fitterT2.ClearPoints();
+
+ for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+// printf("\nLayer %d\n", iLayer);
+// cseed[iLayer].Print();
+ if (!cseed[iLayer].IsOK()) continue;
+ Double_t tilt = cseed[iLayer].fTilt;
+
+ for (Int_t itime = 0; itime < nTimeBins+1; itime++) {
+// printf("\ttime %d\n", itime);
+ if (!cseed[iLayer].fUsable[itime]) continue;
+ // x relative to the midle chamber
+ Double_t x = cseed[iLayer].fX[itime] + cseed[iLayer].fX0 - xref2;
+ Double_t y = cseed[iLayer].fY[itime];
+ Double_t z = cseed[iLayer].fZ[itime];
+
+ //
+ // Tilted rieman
+ //
+ Double_t uvt[6];
+ Double_t x2 = cseed[iLayer].fX[itime] + cseed[iLayer].fX0; // Global x
+ Double_t t = 1.0 / (x2*x2 + y*y);
+ uvt[1] = t;
+ uvt[0] = 2.0 * x2 * uvt[1];
+ uvt[2] = 2.0 * tilt * uvt[1];
+ uvt[3] = 2.0 * tilt *uvt[1] * x;
+ uvt[4] = 2.0 * (y + tilt * z) * uvt[1];
+
+ Double_t error = 2.0 * uvt[1];
+ if (terror) {
+ error *= cseed[iLayer].fSigmaY;
+ }
+ else {
+ error *= 0.2; //Default error
+ }
+// printf("\tadd point :\n");
+// for(int i=0; i<5; i++) printf("%f ", uvt[i]);
+// printf("\n");
+ fitterT2.AddPoint(uvt,uvt[4],error);
+ npointsT++;
+
+ }
+
+ }
+ fitterT2.Eval();
+ Double_t rpolz0 = fitterT2.GetParameter(3);
+ Double_t rpolz1 = fitterT2.GetParameter(4);
+
+ //
+ // Linear fitter - not possible to make boundaries
+ // non accept non possible z and dzdx combination
+ //
+ Bool_t acceptablez = kTRUE;
+ for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+ if (cseed[iLayer].IsOK()) {
+ Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].fX0 - xref2);
+ if (TMath::Abs(cseed[iLayer].fZProb - zT2) > cseed[iLayer].fPadLength * 0.5 + 1.0) {
+ acceptablez = kFALSE;
+ }
+ }
+ }
+ if (!acceptablez) {
+ Double_t zmf = cseed[2].fZref[0] + cseed[2].fZref[1] * (xref2 - cseed[2].fX0);
+ Double_t dzmf = (cseed[2].fZref[1] + cseed[3].fZref[1]) * 0.5;
+ fitterT2.FixParameter(3,zmf);
+ fitterT2.FixParameter(4,dzmf);
+ fitterT2.Eval();
+ fitterT2.ReleaseParameter(3);
+ fitterT2.ReleaseParameter(4);
+ rpolz0 = fitterT2.GetParameter(3);
+ rpolz1 = fitterT2.GetParameter(4);
+ }
+
+ Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);
+ Double_t params[3];
+ params[0] = fitterT2.GetParameter(0);
+ params[1] = fitterT2.GetParameter(1);
+ params[2] = fitterT2.GetParameter(2);
+ Double_t curvature = 1.0 + params[1] * params[1] - params[2] * params[0];
+
+ for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+
+ Double_t x = cseed[iLayer].fX0;
+ Double_t y = 0;
+ Double_t dy = 0;
+ Double_t z = 0;
+ Double_t dz = 0;
+
+ // y
+ Double_t res2 = (x * params[0] + params[1]);
+ res2 *= res2;
+ res2 = 1.0 - params[2]*params[0] + params[1]*params[1] - res2;
+ if (res2 >= 0) {
+ res2 = TMath::Sqrt(res2);
+ y = (1.0 - res2) / params[0];
+ }
+
+ //dy
+ Double_t x0 = -params[1] / params[0];
+ if (-params[2]*params[0] + params[1]*params[1] + 1 > 0) {
+ Double_t rm1 = params[0] / TMath::Sqrt(-params[2]*params[0] + params[1]*params[1] + 1);
+ if (1.0/(rm1*rm1) - (x-x0) * (x-x0) > 0.0) {
+ Double_t res = (x - x0) / TMath::Sqrt(1.0 / (rm1*rm1) - (x-x0)*(x-x0));
+ if (params[0] < 0) res *= -1.0;
+ dy = res;
+ }
+ }
+ z = rpolz0 + rpolz1 * (x - xref2);
+ dz = rpolz1;
+ cseed[iLayer].fYref[0] = y;
+ cseed[iLayer].fYref[1] = dy;
+ cseed[iLayer].fZref[0] = z;
+ cseed[iLayer].fZref[1] = dz;
+ cseed[iLayer].fC = curvature;
+
+ }
+
+ return chi2TR;
+
+}
+
+//___________________________________________________________________
+void AliTRDseedV1::Print()
+{
+ //
+ // Printing the seedstatus
+ //
+
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ Int_t nTimeBins = cal->GetNumberOfTimeBins();
+
+ printf("Seed status :\n");
+ printf(" fTilt = %f\n", fTilt);
+ printf(" fPadLength = %f\n", fPadLength);
+ printf(" fX0 = %f\n", fX0);
+ for(int ic=0; ic<nTimeBins; ic++) {
+ const Char_t *isUsable = fUsable[ic]?"Yes":"No";
+ printf(" %d X[%f] Y[%f] Z[%f] Indexes[%d] clusters[%#x] usable[%s]\n"
+ , ic
+ , fX[ic]
+ , fY[ic]
+ , fZ[ic]
+ , fIndexes[ic]
+ , ((Int_t) fClusters[ic])
+ , isUsable);
+ }
+
+ printf(" fYref[0] =%f fYref[1] =%f\n", fYref[0], fYref[1]);
+ printf(" fZref[0] =%f fZref[1] =%f\n", fZref[0], fZref[1]);
+ printf(" fYfit[0] =%f fYfit[1] =%f\n", fYfit[0], fYfit[1]);
+ printf(" fYfitR[0]=%f fYfitR[1]=%f\n", fYfitR[0], fYfitR[1]);
+ printf(" fZfit[0] =%f fZfit[1] =%f\n", fZfit[0], fZfit[1]);
+ printf(" fZfitR[0]=%f fZfitR[1]=%f\n", fZfitR[0], fZfitR[1]);
+ printf(" fSigmaY =%f\n", fSigmaY);
+ printf(" fSigmaY2=%f\n", fSigmaY2);
+ printf(" fMeanz =%f\n", fMeanz);
+ printf(" fZProb =%f\n", fZProb);
+ printf(" fLabels[0]=%d fLabels[1]=%d\n", fLabels[0], fLabels[1]);
+ printf(" fN =%d\n", fN);
+ printf(" fN2 =%d (>8 isOK)\n",fN2);
+ printf(" fNUsed =%d\n", fNUsed);
+ printf(" fFreq =%d\n", fFreq);
+ printf(" fNChange=%d\n", fNChange);
+ printf(" fMPads =%f\n", fMPads);
+
+ printf(" fC =%f\n", fC);
+ printf(" fCC =%f\n",fCC);
+ printf(" fChi2 =%f\n", fChi2);
+ printf(" fChi2Z =%f\n", fChi2Z);
+
+}
--- /dev/null
+#ifndef ALITRDSEEDV1_H
+#define ALITRDSEEDV1_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+
+////////////////////////////////////////////////////////////////////////////
+// //
+// The TRD track seed //
+// //
+////////////////////////////////////////////////////////////////////////////
+
+#ifndef ALITRDSEED_H
+#include "AliTRDseed.h"
+#endif
+
+#ifndef ALIRIEMAN_H
+#include "AliRieman.h"
+#endif
+
+class TTreeSRedirector;
+
+class AliRieman;
+
+class AliTRDstackLayer;
+class AliTRDcluster;
+class AliTRDrecoParam;
+
+class AliTRDseedV1 : public AliTRDseed
+{
+
+ public:
+
+ AliTRDseedV1(Int_t layer = -1, AliTRDrecoParam *p=0x0);
+ ~AliTRDseedV1();
+ AliTRDseedV1(const AliTRDseedV1 &ref, Bool_t owner=kFALSE);
+ AliTRDseedV1& operator=(const AliTRDseedV1 &);
+
+ Bool_t AttachClustersIter(AliTRDstackLayer *layer, Float_t quality, Bool_t kZcorr = kFALSE, AliTRDcluster *c=0x0);
+ Bool_t AttachClustersProj(AliTRDstackLayer *layer, Float_t quality, Bool_t kZcorr = kFALSE, AliTRDcluster *c=0x0);
+ static Float_t FitRiemanTilt(AliTRDseedV1 * cseed, Bool_t terror);
+ Bool_t FitTracklet();
+
+// Bool_t AttachClusters(Double_t *dx, Float_t quality, Bool_t kZcorr=kFALSE, AliTRDcluster *c=0x0);
+ inline Float_t GetChi2Z(const Float_t z = 0.) const;
+ inline Float_t GetChi2Y(const Float_t y = 0.) const;
+ Float_t GetQuality(Bool_t kZcorr) const;
+ Int_t GetLayer() const { return fLayer; }
+
+ inline void Update(const AliRieman *rieman);
+ void Print(Option_t * /*o*/) const { }
+ void Print();
+
+ void SetLayer(Int_t l) { fLayer = l; }
+ void SetNTimeBins(Int_t nTB) { fTimeBins = nTB; }
+ void SetRecoParam(AliTRDrecoParam *p) { fRecoParam = p; }
+
+ protected:
+
+ void Copy(TObject &ref) const;
+
+ private:
+
+ Int_t fLayer; // layer for this seed
+ Int_t fTimeBins; // local copy of the DB info
+ Bool_t fOwner; // owner of the clusters
+ AliTRDrecoParam *fRecoParam; //! local copy of the reco params
+
+ ClassDef(AliTRDseedV1, 1) // New TRD seed
+
+};
+
+//____________________________________________________________
+inline Float_t AliTRDseedV1::GetChi2Z(const Float_t z) const
+{
+ Float_t z1 = (z == 0.) ? fMeanz : z;
+ Float_t chi = fZref[0] - z1;
+ return chi*chi;
+}
+
+//____________________________________________________________
+inline Float_t AliTRDseedV1::GetChi2Y(const Float_t y) const
+{
+ Float_t y1 = (y == 0.) ? fYfitR[0] : y;
+ Float_t chi = fYref[0] - y1;
+ return chi*chi;
+}
+
+//____________________________________________________________
+inline void AliTRDseedV1::Update(const AliRieman *rieman)
+{
+ fZref[0] = rieman->GetZat(fX0);
+ fZref[1] = rieman->GetDZat(fX0);
+ fYref[0] = rieman->GetYat(fX0);
+ fYref[1] = rieman->GetDYat(fX0);
+}
+
+#endif
+
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/* $Id$ */
+
+///////////////////////////////////////////////////////////////////////////////
+// //
+// Organization of clusters at the level of 1 TRD chamber. //
+// The data structure is used for tracking at the stack level. //
+// //
+// Functionalities: //
+// 1. cluster organization and sorting //
+// 2. fast data navigation //
+// //
+// Authors: //
+// Alex Bercuci <A.Bercuci@gsi.de> //
+// Markus Fasel <M.Fasel@gsi.de> //
+// //
+///////////////////////////////////////////////////////////////////////////////
+
+#include <TObject.h>
+#include <TROOT.h>
+#include <TMath.h>
+#include <TStopwatch.h>
+#include <TTreeStream.h>
+
+#include "AliLog.h"
+#include "AliTRDcluster.h"
+
+#include "AliTRDstackLayer.h"
+#include "AliTRDrecoParam.h"
+#include "AliTRDReconstructor.h"
+
+#define DEBUG
+
+ClassImp(AliTRDstackLayer)
+
+//_____________________________________________________________________________
+AliTRDstackLayer::AliTRDstackLayer(Double_t z0, Double_t zLength
+ , UChar_t stackNr, AliTRDrecoParam *p)
+ :AliTRDpropagationLayer()
+ ,fOwner(kFALSE)
+ ,fStackNr(stackNr)
+ ,fNRows(kMaxRows)
+ ,fZ0(z0)
+ ,fZLength(zLength)
+ ,fRecoParam(p)
+ ,fDebugStream(0x0)
+{
+ //
+ // Default constructor (Only provided to use AliTRDstackLayer with arrays)
+ //
+
+ for(int i=0; i<kMaxRows; i++) fPositions[i] = 0;
+}
+
+//_____________________________________________________________________________
+AliTRDstackLayer::AliTRDstackLayer(const AliTRDpropagationLayer &layer, Double_t
+z0, Double_t zLength, UChar_t stackNr, AliTRDrecoParam *p):
+ AliTRDpropagationLayer(layer)
+ ,fOwner(kFALSE)
+ ,fStackNr(stackNr)
+ ,fNRows(kMaxRows)
+ ,fZ0(z0)
+ ,fZLength(zLength)
+ ,fRecoParam(p)
+ ,fDebugStream(0x0)
+{
+// Standard constructor.
+// Initialize also the underlying AliTRDpropagationLayer using the copy constructor.
+
+ for(int i=0; i<kMaxRows; i++) fPositions[i] = 0;
+}
+
+//_____________________________________________________________________________
+AliTRDstackLayer::AliTRDstackLayer(const AliTRDpropagationLayer &layer):
+ AliTRDpropagationLayer(layer)
+ ,fOwner(kFALSE)
+ ,fStackNr(0)
+ ,fNRows(kMaxRows)
+ ,fZ0(0)
+ ,fZLength(0)
+ ,fRecoParam(0x0)
+ ,fDebugStream(0x0)
+{
+// Standard constructor using only AliTRDpropagationLayer.
+
+ for(int i=0; i<kMaxRows; i++) fPositions[i] = 0;
+}
+
+//_____________________________________________________________________________
+AliTRDstackLayer::AliTRDstackLayer(const AliTRDstackLayer &layer):
+ AliTRDpropagationLayer(layer)
+ ,fOwner(layer.fOwner)
+ ,fStackNr(layer.fStackNr)
+ ,fNRows(layer.fNRows)
+ ,fZ0(layer.fZ0)
+ ,fZLength(layer.fZLength)
+ ,fRecoParam(layer.fRecoParam)
+ ,fDebugStream(layer.fDebugStream)
+{
+// Copy Constructor (performs a deep copy)
+
+ for(Int_t i = 0; i < kMaxRows; i++) fPositions[i] = layer.fPositions[i];
+// BuildIndices();
+}
+
+//_____________________________________________________________________________
+AliTRDstackLayer &AliTRDstackLayer::operator=(const AliTRDpropagationLayer &layer)
+{
+// Assignment operator from an AliTRDpropagationLayer
+
+ if (this != &layer) layer.Copy(*this);
+ return *this;
+}
+
+//_____________________________________________________________________________
+AliTRDstackLayer &AliTRDstackLayer::operator=(const AliTRDstackLayer &layer)
+{
+// Assignment operator
+
+ if (this != &layer) layer.Copy(*this);
+ return *this;
+}
+
+//_____________________________________________________________________________
+void AliTRDstackLayer::Copy(TObject &o) const
+{
+// Copy method. Performs a deep copy of all data from this object to object o.
+
+ AliTRDstackLayer &layer = (AliTRDstackLayer &)o;
+ layer.fZ0 = fZ0;
+ layer.fOwner = kFALSE;
+ layer.fNRows = fNRows;
+ layer.fZLength = fZLength;
+ layer.fStackNr = fStackNr;
+ layer.fDebugStream = fDebugStream;
+ layer.fRecoParam = fRecoParam;
+
+ AliTRDpropagationLayer::Copy(layer); // copies everything into layer
+ for(UChar_t i = 0; i < kMaxRows; i++) layer.fPositions[i] = 0;
+// layer.BuildIndices();
+}
+
+//_____________________________________________________________________________
+AliTRDstackLayer::~AliTRDstackLayer()
+{
+// Destructor
+ if(fOwner) for(int ic=0; ic<fN; ic++) delete fClusters[ic];
+}
+
+//_____________________________________________________________________________
+void AliTRDstackLayer::SetRange(const Float_t z0, const Float_t zLength)
+{
+// Sets the range in z-direction
+//
+// Parameters:
+// z0 : starting position of layer in the z direction
+// zLength : length of layer in the z direction
+
+ fZ0 = (z0 <= z0 + zLength) ? z0 : z0 + zLength;
+ fZLength = TMath::Abs(zLength);
+}
+
+//_____________________________________________________________________________
+void AliTRDstackLayer::BuildIndices(Int_t iter)
+{
+// Rearrangement of the clusters belonging to the propagation layer for the stack.
+//
+// Detailed description
+//
+// The array indices of all clusters in one PropagationLayer are stored in
+// array. The array is divided into several bins.
+// The clusters are sorted in increasing order of their y coordinate.
+//
+// Sorting algorithm: TreeSearch
+//
+
+ if(!fN) return;
+
+ // Select clusters that belong to the Stack
+ Int_t nClStack = 0; // Internal counter
+ for(Int_t i = 0; i < fN; i++){
+ Double_t zval = fClusters[i]->GetZ();
+ if(zval < fZ0 || zval > fZ0 + fZLength || fClusters[i]->IsUsed()){
+ fClusters[i] = 0x0;
+ fIndex[i] = 9999;
+ } else nClStack++;
+ }
+ if(nClStack > kMaxClustersLayer) AliWarning(Form("Number of clusters in stack %d exceed buffer size %d", nClStack, kMaxClustersLayer));
+
+ // Nothing in this Stack
+ if(!nClStack){
+ delete fClusters;
+ delete fIndex;
+ fClusters = 0x0;
+ fIndex = 0x0;
+ fN = 0;
+ memset(fPositions, 0, sizeof(UChar_t) * 16);
+ return;
+ }
+
+ // Make a copy
+ AliTRDcluster *helpCL[kMaxClustersLayer];
+ Int_t helpInd[kMaxClustersLayer];
+ nClStack = 0;
+ for(Int_t i = 0; i < TMath::Min(fN, kMaxClustersLayer); i++){
+ if(!fClusters[i]) continue;
+ helpCL[nClStack] = fClusters[i];
+ helpInd[nClStack] = fIndex[i];
+ fClusters[i] = 0x0;
+ fIndex[i] = 9999;
+ nClStack++;
+ }
+
+ // do clusters arrangement
+ fN = nClStack;
+ nClStack = 0;
+ memset(fPositions,0,sizeof(UChar_t)*16); // Reset Positions array
+ for(Int_t i = 0; i < fN; i++){
+ // boundarie check
+ AliTRDcluster *cl = helpCL[i];
+ Double_t zval = cl->GetZ();
+ UChar_t TreeIndex = (UChar_t)(TMath::Abs(fZ0 - zval)/fZLength * fNRows);
+ if(TreeIndex > fNRows - 1) TreeIndex = fNRows - 1;
+ // Insert Leaf
+ Int_t pos = FindYPosition(cl->GetY(), TreeIndex, i);
+ if(pos == -1){ // zbin is empty;
+ Int_t upper = (TreeIndex == fNRows - 1) ? nClStack : fPositions[TreeIndex + 1];
+ memmove(fClusters + upper + 1, fClusters + upper, (sizeof(AliTRDcluster *))*(nClStack-upper));
+ memmove(fIndex + upper + 1, fIndex + upper, (sizeof(UInt_t))*(nClStack-upper));
+ fClusters[upper] = cl;
+ fIndex[upper] = helpInd[i];
+ // Move All pointer one position back
+ for(UChar_t j = TreeIndex + 1; j < fNRows; j++) fPositions[j]++;
+ nClStack++;
+ } else { // zbin not empty
+ memmove(fClusters + pos + 2, fClusters + pos+1, (sizeof(AliTRDcluster *))*(nClStack-(pos+1)));
+ memmove(fIndex + pos + 2, fIndex + pos+1, (sizeof(UInt_t))*(nClStack-(pos+1)));
+ fClusters[pos + 1] = cl; //fIndex[i];
+ fIndex[pos + 1] = helpInd[i];
+ // Move All pointer one position back
+ for(UChar_t j = TreeIndex + 1; j < fNRows; j++) fPositions[j]++;
+ nClStack++;
+ }
+
+
+ // Debug Streaming
+#ifdef DEBUG
+ if(fDebugStream && AliTRDReconstructor::StreamLevel() >= 3){
+ TTreeSRedirector &cstream = *fDebugStream;
+ cstream << "BuildIndices"
+ << "StackNr=" << fStackNr
+ << "SectorNr=" << fSec
+ << "Iter=" << iter
+ << "C.=" << cl
+ << "TreeIdx=" << TreeIndex
+ << "\n";
+ }
+#endif
+ }
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDstackLayer::FindYPosition(Double_t y, UChar_t z, Int_t nClusters) const
+{
+//
+// Tree search Algorithm to find the nearest left cluster for a given
+// y-position in a certain z-bin (in fact AVL-tree).
+// Making use of the fact that clusters are sorted in y-direction.
+//
+// Parameters:
+// y : y position of the reference point in tracking coordinates
+// z : z reference bin.
+// nClusters :
+//
+// Output :
+// Index of the nearest left cluster in the StackLayer indexing (-1 if no clusters are found)
+//
+
+ Int_t start = fPositions[z]; // starting Position of the bin
+ Int_t upper = (Int_t)((z != fNRows - 1) ? fPositions[z+1] : nClusters); // ending Position of the bin
+ Int_t end = upper - 1; // ending Position of the bin
+ if(end < start) return -1; // Bin is empty
+ Int_t middle = static_cast<Int_t>((start + end)/2);
+ // 1st Part: climb down the tree: get the next cluster BEFORE ypos
+ while(start + 1 < end){
+ if(y >= fClusters[middle]->GetY()) start = middle;
+ else end = middle;
+ middle = static_cast<Int_t>((start + end)/2);
+ }
+ if(y > fClusters[end]->GetY()) return end;
+ return start;
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDstackLayer::FindNearestYCluster(Double_t y, UChar_t z) const
+{
+//
+// Tree search Algorithm to find the nearest cluster for a given
+// y-position in a certain z-bin (in fact AVL-tree).
+// Making use of the fact that clusters are sorted in y-direction.
+//
+// Parameters:
+// y : y position of the reference point in tracking coordinates
+// z : z reference bin.
+//
+// Output
+// Index of the nearest cluster in the StackLayer indexing (-1 if no clusters are found)
+//
+
+ Int_t position = FindYPosition(y, z, fN);
+ if(position == -1) return position; // bin empty
+ // FindYPosition always returns the left Neighbor. We don't know if the left or the right Neighbor is nearest
+ // to the Reference y-position, so test both
+ Int_t upper = (Int_t)((z < fNRows-1) ? fPositions[z+1] : fN); // ending Position of the bin
+ if((position + 1) < (upper)){
+ if(TMath::Abs(y - fClusters[position + 1]->GetY()) < TMath::Abs(y - fClusters[position]->GetY())) return position + 1;
+ else return position;
+ }
+ return position;
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDstackLayer::SearchNearestCluster(Double_t y, Double_t z, Double_t maxroady, Double_t maxroadz) const
+{
+//
+// Finds the nearest cluster from a given point in a defined range.
+// Distance is determined in a 2D space by the 2-Norm.
+//
+// Parameters :
+// y : y position of the reference point in tracking coordinates
+// z : z reference bin.
+// maxroady : maximum searching distance in y direction
+// maxroadz : maximum searching distance in z direction
+//
+// Output
+// Index of the nearest cluster in the StackLayer indexing (-1 if no cluster is found).
+// Cluster can be accessed with the operator[] or GetCluster(Int_t index)
+//
+// Detail description
+//
+// The following steps are perfomed:
+// 1. Get the expected z bins inside maxroadz.
+// 2. For each z bin find nearest y cluster.
+// 3. Select best candidate
+//
+ Int_t index = -1;
+ // initial minimal distance will be represented as ellipse: semi-major = z-direction
+ // later 2-Norm will be used
+// Float_t nExcentricity = TMath::Sqrt(maxroadz*maxroadz - maxroad*maxroad)/maxroadz;
+ Float_t mindist = maxroadz;
+
+ // not very nice but unfortunately neccessarry: we have ho check the neighbors in both directions (+ and -) too. How
+ // much neighbors depends on the Quotient maxroadz/fZLength
+ UChar_t maxRows = 3;
+ UChar_t zpos[kMaxRows];
+ // Float_t mindist = TMath::Sqrt(maxroad*maxroad + maxroadz*maxroadz);
+// UChar_t myZbin = FindTreePosition(z, fZ0 + fZLength/2, fZLength/4, 8, 8, kFALSE);
+ UChar_t myZbin = (UChar_t)(TMath::Abs(fZ0 - z)/fZLength * fNRows);
+ if(z < fZ0) myZbin = 0;
+ if(z > fZ0 + fZLength) myZbin = fNRows - 1;
+
+ UChar_t nNeighbors = 0;
+ for(UChar_t i = 0; i < maxRows; i++){
+ if((myZbin - 1 + i) < 0) continue;
+ if((myZbin - 1 + i) > fNRows - 1) break;
+ zpos[nNeighbors] = myZbin - 1 + i;
+ nNeighbors++;
+ }
+ Float_t ycl = 0, zcl = 0;
+ for(UChar_t neighbor = 0; neighbor < nNeighbors; neighbor++){ // Always test the neighbors too
+ Int_t pos = FindNearestYCluster(y,zpos[neighbor]);
+ if(pos == -1) continue; // No cluster in bin
+ AliTRDcluster *c = (AliTRDcluster *) (fClusters[pos]);
+ if(c->IsUsed()) continue; // we are only interested in unused clusters
+ ycl = c->GetY();
+ // Too far away in y-direction (Prearrangement)
+ if (TMath::Abs(ycl - y) > maxroady) continue;
+
+ zcl = c->GetZ();
+ // Too far away in z-Direction
+ // (Prearrangement since we have not so many bins to test)
+ if (TMath::Abs(zcl - z) > maxroadz) continue;
+
+ Float_t dist; // distance defined as 2-Norm
+ // if we havent found a Particle that is in the ellipse around (y,z) with maxroad as semi-minor and
+ // maxroadz as semi-major, we take the radius of the ellipse concerning the cluster as mindist, later we
+ // take the 2-Norm when we found a cluster inside the ellipse (The value 10000 is taken because it is surely
+ // large enough to be usable as an indicator whether we have found a nearer cluster or not)
+// if(mindist > 10000.){
+// Float_t phi = ((zcl - z) == 0) ? TMath::Pi()/2 : TMath::ATan((ycl - y)/(zcl - z));
+// mindist = maxroad/TMath::Sqrt(1 - nExcentricity*nExcentricity * (TMath::Cos(phi))*(TMath::Cos(phi)));
+// }
+ dist = TMath::Max(TMath::Abs(y-ycl),TMath::Abs(z-zcl)); // infinity Norm
+// dist = TMath::Sqrt((ycl - y)*(ycl - y) + (zcl - z)*(zcl - z));
+ if((Int_t)(dist * 100000) < (Int_t)(mindist * 100000)){
+ //if((dist = TMath::Sqrt((ycl - y)*(ycl - y) + (zcl - z)*(zcl - z))) < mindist){
+ mindist = dist;
+ index = pos;
+ }
+ }
+ // This is the Array Position in fIndex2D of the Nearest cluster: if a
+ // cluster is called, then the function has to retrieve the Information
+ // which is Stored in the Array called, the function
+ return index;
+}
+
+//_____________________________________________________________________________
+void AliTRDstackLayer::BuildCond(AliTRDcluster *cl, Double_t *cond, UChar_t Layer, Double_t theta, Double_t phi)
+{
+// Helper function to calculate the area where to expect a cluster in THIS
+// layer.
+//
+// Parameters :
+// cl :
+// cond :
+// Layer :
+// theta :
+// phi :
+//
+// Detail description
+//
+// Helper function to calculate the area where to expect a cluster in THIS
+// layer. by using the information of a former cluster in another layer
+// and the angle in theta- and phi-direction between layer 0 and layer 3.
+// If the layer is zero, initial conditions are calculated. Otherwise a
+// linear interpolation is performed.
+//Begin_Html
+//<img src="gif/build_cond.gif">
+//End_Html
+//
+
+ if(!fRecoParam){
+ AliError("Reconstruction parameters not initialized.");
+ return;
+ }
+
+ if(Layer == 0){
+ cond[0] = cl->GetY(); // center: y-Direction
+ cond[1] = cl->GetZ(); // center: z-Direction
+ cond[2] = fRecoParam->GetMaxPhi() * (cl->GetX() - GetX()) + 1.0; // deviation: y-Direction
+ cond[3] = fRecoParam->GetMaxTheta() * (cl->GetX() - GetX()) + 1.0; // deviation: z-Direction
+ } else {
+ cond[0] = cl->GetY() + phi * (GetX() - cl->GetX());
+ cond[1] = cl->GetZ() + theta * (GetX() - cl->GetX());
+ cond[2] = fRecoParam->GetRoad0y() + phi;
+ cond[3] = fRecoParam->GetRoad0z();
+ }
+}
+
+//_____________________________________________________________________________
+void AliTRDstackLayer::GetClusters(Double_t *cond, Int_t *index, Int_t& ncl, Int_t BufferSize)
+{
+// Finds all clusters situated in this layer inside a rectangle given by the center an ranges.
+//
+// Parameters :
+// cond :
+// index :
+// ncl :
+// BufferSize :
+//
+// Output :
+//
+// Detail description
+//
+// Function returs an array containing the indices in the stacklayer of
+// the clusters found an the number of found clusters in the stacklayer
+
+ ncl = 0;
+ memset(index, 0, BufferSize*sizeof(Int_t));
+ if(fN == 0) return;
+
+ //Boundary checks
+ Double_t zvals[2];
+ zvals[0] = ((cond[1] - cond[3]) < fZ0) ? fZ0 : (cond[1] - cond[3]);
+ zvals[1] = ((cond[1] + cond[3]) < fZ0 + fZLength) ? (cond[1] + cond[3]) : fZ0 + fZLength;
+
+ UChar_t zlo = (fZ0>zvals[0]) ? 0 : (UChar_t)(TMath::Abs(fZ0 - zvals[0])/fZLength * fNRows);
+ UChar_t zhi = (fZ0+fZLength<zvals[1]) ? fNRows - 1 : (UChar_t)(TMath::Abs(fZ0 - zvals[1])/fZLength * fNRows);
+
+ //Preordering in Direction z saves a lot of loops (boundary checked)
+ for(UChar_t z = zlo; z <= zhi; z++){
+ UInt_t upper = (z != fNRows-1) ? fPositions[z+1] : fN;
+ for(Int_t y = fPositions[z]; y < (Int_t)upper; y++){
+ if(ncl == BufferSize){
+ AliWarning("Buffer size riched. Some clusters may be lost.");
+ return; //Buffer filled
+ }
+ if(fClusters[y]->GetY() > (cond[0] + cond[2])) break; // Abbortion conditions!!!
+ if(fClusters[y]->GetY() < (cond[0] - cond[2])) continue; // Too small
+ if(((Int_t)((fClusters[y]->GetZ())*1000) < (Int_t)(zvals[0]*1000)) || ((Int_t)((fClusters[y]->GetZ())*1000) > (Int_t)(zvals[1]*1000))) continue;
+ index[ncl] = y;
+ ncl++;
+ }
+ }
+ if(ncl>fN) AliError(Form("Clusters found %d > %d (clusters in layer)", ncl, fN));
+}
+
+//_____________________________________________________________________________
+AliTRDcluster *AliTRDstackLayer::GetNearestCluster(Double_t *cond)
+{
+// Function returning a pointer to the nearest cluster (nullpointer if not successfull).
+//
+// Parameters :
+// cond :
+//
+// Output :
+// pointer to the nearest cluster (nullpointer if not successfull).
+//
+// Detail description
+//
+// returns a pointer to the nearest cluster (nullpointer if not
+// successfull) by the help of the method FindNearestCluster
+
+
+ Double_t maxroad = fRecoParam->GetRoad2y();
+ Double_t maxroadz = fRecoParam->GetRoad2z();
+
+ Int_t index = SearchNearestCluster(cond[0],cond[1],maxroad,maxroadz);
+ AliTRDcluster *returnCluster = 0x0;
+ if(index != -1) returnCluster = (AliTRDcluster *) fClusters[index];
+ return returnCluster;
+}
+
+//_____________________________________________________________________________
+void AliTRDstackLayer::PrintClusters() const
+{
+// Prints the position of each cluster in the stacklayer on the stdout
+//
+ printf("fDebugStream = %p\n", fDebugStream);
+ printf("fRecoParam = %p\n", fRecoParam);
+
+ for(Int_t i = 0; i < fN; i++){
+ printf("AliTRDstackLayer: index=%i, Cluster: X = %3.3f, Y = %3.3f, Z = %3.3f\n", i, fClusters[i]->GetX(),fClusters[i]->GetY(),fClusters[i]->GetZ());
+ if(fClusters[i]->IsUsed()) printf("cluster allready used. rejected in search algorithm\n");
+ }
+}
--- /dev/null
+#ifndef ALITRDSTACKLAYER_H
+#define ALITRDSTACKLAYER_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+
+////////////////////////////////////////////////////////////////////////////
+// //
+// A TRD layer in a single stack //
+// //
+////////////////////////////////////////////////////////////////////////////
+
+#ifndef ALITRDPROPAGATIONLAYER_H
+#include "AliTRDpropagationLayer.h"
+#endif
+
+#ifndef ALITRDCLUSTER_H
+#include "AliTRDcluster.h"
+#endif
+
+class AliTRDrecoParam;
+
+class AliTRDstackLayer : public AliTRDpropagationLayer
+{
+ public:
+ enum{
+ kMaxClustersLayer = 150,
+ kMaxRows = 16
+ };
+
+ AliTRDstackLayer(Double_t z0 = 0., Double_t zLength = 0., UChar_t stackNr = 0
+ , AliTRDrecoParam *p=0x0);
+ AliTRDstackLayer(const AliTRDpropagationLayer &layer, Double_t z0
+ , Double_t zLength, UChar_t stackNr, AliTRDrecoParam *p = 0x0);
+ AliTRDstackLayer(const AliTRDpropagationLayer &layer);
+ AliTRDstackLayer(const AliTRDstackLayer &layer);
+ ~AliTRDstackLayer();
+ AliTRDstackLayer &operator=(const AliTRDpropagationLayer &myLayer);
+ AliTRDstackLayer &operator=(const AliTRDstackLayer &myLayer);
+ AliTRDcluster *operator[](const Int_t i) {
+ return ((i < fN) && (i >= 0)) ? fClusters[i] : 0x0;
+ }
+
+ void BuildIndices(Int_t iter = 0);
+ void BuildCond(AliTRDcluster *cl, Double_t *cond, UChar_t Layer, Double_t theta=0., Double_t phi=0.);
+ AliTRDcluster* GetCluster(Int_t index) const {return index < fN ? fClusters[index] : 0x0;}
+ Int_t GetGlobalIndex(const Int_t index) const {return ((index < fN) && (index >= 0)) ? fIndex[index] : 0; }
+ void GetClusters(Double_t *cond, Int_t *index, Int_t& ncl, Int_t BufferSize = kMaxClustersLayer);
+ AliTRDcluster* GetNearestCluster(Double_t *cond);
+
+ Double_t GetZ0() const { return fZ0; }
+ Double_t GetDZ0() const { return fZLength; }
+ Int_t GetNClusters() const { return fN; }
+ UInt_t GetStackNr() const { return fStackNr; }
+ void PrintClusters() const;
+ Int_t SearchNearestCluster(const Double_t y, const Double_t z, const Double_t Roady, const Double_t Roadz) const;
+ void SetRange(Float_t z0, Float_t zLength);
+ void SetNRows(const Int_t nRows){ fNRows = nRows; }
+ void SetStackNr(const UInt_t stackNr){ fStackNr = stackNr; }
+ void SetOwner(Bool_t own = kTRUE) {fOwner = own;}
+ void SetClustersArray(AliTRDcluster **cl, Int_t nClusters){fClusters = cl; fN = nClusters;}
+ void SetIndexArray(UInt_t *indexArray){fIndex = indexArray;}
+ void SetDebugStream(TTreeSRedirector *debug) {fDebugStream = debug;}
+ void SetRecoParam(AliTRDrecoParam *p) {fRecoParam = p;}
+
+private:
+ void Copy(TObject &o) const;
+ Int_t FindYPosition(Double_t y, UChar_t z, Int_t nClusters) const;
+ Int_t FindNearestYCluster(Double_t y, UChar_t z) const;
+
+private:
+ Bool_t fOwner; // owner of the clusters
+ UChar_t fStackNr; // stack number in supermodule
+ UChar_t fNRows; // number of pad rows in the chamber
+ UChar_t fPositions[kMaxRows]; // starting index of clusters in pad row
+ Double_t fZ0; // starting position of the layer in Z direction
+ Double_t fZLength; // length of the layer in Z direction
+ AliTRDrecoParam *fRecoParam; //! reconstruction parameters
+ TTreeSRedirector *fDebugStream; //! debug streamer
+
+ ClassDef(AliTRDstackLayer, 1) // stack propagation layer
+
+};
+#endif // ALITRDSTACKLAYER_H_
+
/* $Id$ */
-//#include <Riostream.h>
-
-//#include <TMath.h>
#include <TVector2.h>
#include "AliTracker.h"
#include "AliTRDgeometry.h"
#include "AliTRDcluster.h"
#include "AliTRDtrack.h"
-//#include "AliTRDtracklet.h"
#include "AliTRDcalibDB.h"
#include "Cal/AliTRDCalPID.h"
}
tb = cluster->GetLocalTimeBin();
- if ((tb == 0) || (tb >= ntb)) {
- AliWarning(Form("time bin 0 or > %d in cluster %d", ntb, iClus));
+ if ((tb < 0) || (tb >= ntb)) {
+ AliWarning(Form("time bin < 0 or > %d in cluster %d", ntb, iClus));
AliInfo(Form("dQ/dl %f", fdQdl[iClus]));
continue;
}
return lastplane;
}
-
-//_____________________________________________________________________________
-Int_t AliTRDtracker::Clusters2Tracks(AliESDEvent *event)
-{
- //
- // Finds tracks within the TRD. The ESD event is expected to contain seeds
- // at the outer part of the TRD. The seeds
- // are found within the TRD if fAddTRDseeds is TRUE.
- // The tracks are propagated to the innermost time bin
- // of the TRD and the ESD event is updated
- //
-
- Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
- Float_t foundMin = fgkMinClustersInTrack * timeBins;
- Int_t nseed = 0;
- Int_t found = 0;
- //Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
-
- Int_t n = event->GetNumberOfTracks();
- for (Int_t i = 0; i < n; i++) {
-
- AliESDtrack *seed = event->GetTrack(i);
- ULong_t status = seed->GetStatus();
- if ((status & AliESDtrack::kTRDout) == 0) {
- continue;
- }
- if ((status & AliESDtrack::kTRDin) != 0) {
- continue;
- }
- nseed++;
-
- AliTRDtrack *seed2 = new AliTRDtrack(*seed);
- //seed2->ResetCovariance();
- AliTRDtrack *pt = new AliTRDtrack(*seed2,seed2->GetAlpha());
- AliTRDtrack &t = *pt;
- FollowProlongation(t);
- if (t.GetNumberOfClusters() >= foundMin) {
- UseClusters(&t);
- CookLabel(pt,1 - fgkLabelFraction);
- //t.CookdEdx();
- }
- found++;
-
- Double_t xTPC = 250.0;
- if (PropagateToX(t,xTPC,fgkMaxStep)) {
- seed->UpdateTrackParams(pt, AliESDtrack::kTRDin);
- }
- delete seed2;
- delete pt;
-
- }
-
- AliInfo(Form("Number of loaded seeds: %d",nseed));
- AliInfo(Form("Number of found tracks from loaded seeds: %d",found));
- AliInfo(Form("Total number of found tracks: %d",found));
-
- return 0;
-
-}
//_____________________________________________________________________________
Int_t AliTRDtracker::PropagateBack(AliESDEvent *event)
// by the TPC tracker.
//
- Int_t found = 0; // number of tracks found
- Float_t foundMin = 20.0;
-
- // Sort tracks according to covariance of local Y and Z
- Int_t nSeed = event->GetNumberOfTracks();
- Float_t *quality = new Float_t[nSeed];
- Int_t *index = new Int_t[nSeed];
- for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
- AliESDtrack *seed = event->GetTrack(iSeed);
- Double_t covariance[15];
- seed->GetExternalCovariance(covariance);
- quality[iSeed] = covariance[0] + covariance[2];
- }
- TMath::Sort(nSeed,quality,index,kFALSE);
-
- // Backpropagate all seeds
- for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
-
- // Get the seeds in sorted sequence
- AliESDtrack *seed = event->GetTrack(index[iSeed]);
- fHBackfit->Fill(0); // All seeds
-
- // Check the seed status
- ULong_t status = seed->GetStatus();
- if ((status & AliESDtrack::kTPCout) == 0) {
- fHBackfit->Fill(1); // TPC outer edge reached
- continue;
- }
- if ((status & AliESDtrack::kTRDout) != 0) {
- fHBackfit->Fill(2); // TRD outer edge reached (does this happen ?)
- continue;
- }
-
- // Do the back prolongation
- Int_t lbl = seed->GetLabel();
- AliTRDtrack *track = new AliTRDtrack(*seed);
- track->SetSeedLabel(lbl);
- seed->UpdateTrackParams(track,AliESDtrack::kTRDbackup); // Make backup
- fNseeds++;
- Float_t p4 = track->GetC();
- Int_t expectedClr = FollowBackProlongation(*track);
- fHBackfit->Fill(3); // Back prolongation done
- fHX->Fill(track->GetX());
- // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- if ((TMath::Abs(track->GetC() - p4) / TMath::Abs(p4) < 0.2) ||
- (track->Pt() > 0.8)) {
-
- fHBackfit->Fill(4);
-
- //
- // Make backup for back propagation
- //
-
- Int_t foundClr = track->GetNumberOfClusters();
- if (foundClr >= foundMin) {
- track->CookdEdx();
- track->CookdEdxTimBin(seed->GetID()); // A.Bercuci 25.07.07
- CookLabel(track,1 - fgkLabelFraction);
- if (track->GetBackupTrack()) {
- UseClusters(track->GetBackupTrack());
- }
-
- // Sign only gold tracks
- if (track->GetChi2() / track->GetNumberOfClusters() < 4) {
- if ((seed->GetKinkIndex(0) == 0) &&
- (track->Pt() < 1.5)) {
- UseClusters(track);
- }
- }
- Bool_t isGold = kFALSE;
+ Int_t found = 0; // number of tracks found
+ Float_t foundMin = 20.0;
- // Full gold track
- if (track->GetChi2() / track->GetNumberOfClusters() < 5) {
- //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
- if (track->GetBackupTrack()) {
- seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
-
- }
- isGold = kTRUE;
- //fHBackfit->Fill()
- }
-
- // Almost gold track
- if ((!isGold) &&
- (track->GetNCross() == 0) &&
- (track->GetChi2() / track->GetNumberOfClusters() < 7)) {
- //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
- if (track->GetBackupTrack()) {
- seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
- }
- isGold = kTRUE;
- }
-
- if ((!isGold) &&
- (track->GetBackupTrack())) {
- if ((track->GetBackupTrack()->GetNumberOfClusters() > foundMin) &&
- ((track->GetBackupTrack()->GetChi2()/(track->GetBackupTrack()->GetNumberOfClusters()+1)) < 7)) {
- seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
- isGold = kTRUE;
- }
+ Int_t nSeed = event->GetNumberOfTracks();
+ if(!nSeed){
+ // run stand alone tracking
+ if (AliTRDReconstructor::SeedingOn()) Clusters2Tracks(event);
+ return 0;
}
-
- if ((track->StatusForTOF() > 0) &&
- (track->GetNCross() == 0) &&
- (Float_t(track->GetNumberOfClusters()) / Float_t(track->GetNExpected()) > 0.4)) {
- //seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
+
+ Float_t *quality = new Float_t[nSeed];
+ Int_t *index = new Int_t[nSeed];
+ for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
+ AliESDtrack *seed = event->GetTrack(iSeed);
+ Double_t covariance[15];
+ seed->GetExternalCovariance(covariance);
+ quality[iSeed] = covariance[0] + covariance[2];
}
-
- }
- }
- /**/
-
- /**/
- // 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 (track->GetBackupTrack()) {
- cstream << "Tracks"
- << "EventNrInFile=" << eventNrInFile
- << "ESD.=" << seed
- << "trd.=" << track
- << "trdback.=" << track->GetBackupTrack()
- << "\n";
- }
- else {
- cstream << "Tracks"
- << "EventNrInFile=" << eventNrInFile
- << "ESD.=" << seed
- << "trd.=" << track
- << "trdback.=" << track
- << "\n";
- }
- }
- /**/
-
- // Propagation to the TOF (I.Belikov)
- if (track->GetStop() == kFALSE) {
- fHBackfit->Fill(5);
-
- Double_t xtof = 371.0;
- Double_t xTOF0 = 370.0;
-
- Double_t c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
- if (TMath::Abs(c2) >= 0.99) {
+ // Sort tracks according to covariance of local Y and Z
+ TMath::Sort(nSeed,quality,index,kFALSE);
- fHBackfit->Fill(6);
- delete track;
- continue;
- }
-
- PropagateToX(*track,xTOF0,fgkMaxStep);
-
- // Energy losses taken to the account - check one more time
- c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
- if (TMath::Abs(c2) >= 0.99) {
+ // Backpropagate all seeds
+ for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
- fHBackfit->Fill(7);
- delete track;
- continue;
- }
-
- //if (!PropagateToX(*track,xTOF0,fgkMaxStep)) {
- // fHBackfit->Fill(7);
- //delete track;
- // continue;
- //}
-
- Double_t ymax = xtof * TMath::Tan(0.5 * AliTRDgeometry::GetAlpha());
- Double_t y;
- track->GetYAt(xtof,GetBz(),y);
- if (y > ymax) {
- if (!track->Rotate( AliTRDgeometry::GetAlpha())) {
- fHBackfit->Fill(8);
- delete track;
- continue;
- }
- }
- else if (y < -ymax) {
- if (!track->Rotate(-AliTRDgeometry::GetAlpha())) {
- fHBackfit->Fill(9);
- delete track;
- continue;
- }
- }
-
- if (track->PropagateTo(xtof)) {
- seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
- fHBackfit->Fill(10);
-
- for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
- for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
- seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
- }
- seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
- }
- //seed->SetTRDtrack(new AliTRDtrack(*track));
- if (track->GetNumberOfClusters() > foundMin) {
- fHBackfit->Fill(11);
- found++;
+ // Get the seeds in sorted sequence
+ AliESDtrack *seed = event->GetTrack(index[iSeed]);
+ fHBackfit->Fill(0); // All seeds
+
+ // Check the seed status
+ ULong_t status = seed->GetStatus();
+ if ((status & AliESDtrack::kTPCout) == 0) {
+ fHBackfit->Fill(1); // TPC outer edge reached
+ continue;
+ }
+ if ((status & AliESDtrack::kTRDout) != 0) {
+ fHBackfit->Fill(2); // TRD outer edge reached (does this happen ?)
+ continue;
+ }
+
+ // Do the back prolongation
+ Int_t lbl = seed->GetLabel();
+ AliTRDtrack *track = new AliTRDtrack(*seed);
+ track->SetSeedLabel(lbl);
+ seed->UpdateTrackParams(track,AliESDtrack::kTRDbackup); // Make backup
+ fNseeds++;
+ Float_t p4 = track->GetC();
+ Int_t expectedClr = FollowBackProlongation(*track);
+ fHBackfit->Fill(3); // Back prolongation done
+ fHX->Fill(track->GetX());
+
+ if ((TMath::Abs(track->GetC() - p4) / TMath::Abs(p4) < 0.2) ||
+ (track->Pt() > 0.8)) {
+
+ fHBackfit->Fill(4);
+
+ //
+ // Make backup for back propagation
+ //
+
+ Int_t foundClr = track->GetNumberOfClusters();
+ if (foundClr >= foundMin) {
+ track->CookdEdx();
+ track->CookdEdxTimBin(seed->GetID()); // A.Bercuci 25.07.07
+ CookLabel(track,1 - fgkLabelFraction);
+ if (track->GetBackupTrack()) UseClusters(track->GetBackupTrack());
+
+
+ // Sign only gold tracks
+ if (track->GetChi2() / track->GetNumberOfClusters() < 4) {
+ if ((seed->GetKinkIndex(0) == 0) &&
+ (track->Pt() < 1.5)) {
+ UseClusters(track);
+ }
+ }
+ Bool_t isGold = kFALSE;
+
+ // Full gold track
+ if (track->GetChi2() / track->GetNumberOfClusters() < 5) {
+ //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
+ if (track->GetBackupTrack()) {
+ seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
+ }
+ isGold = kTRUE;
+ //fHBackfit->Fill()
+ }
+
+ // Almost gold track
+ if ((!isGold) && (track->GetNCross() == 0) &&
+ (track->GetChi2() / track->GetNumberOfClusters() < 7)) {
+ //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
+ if (track->GetBackupTrack()) {
+ seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
+ }
+ isGold = kTRUE;
+ }
+
+ if ((!isGold) && (track->GetBackupTrack())) {
+ if ((track->GetBackupTrack()->GetNumberOfClusters() > foundMin) && ((track->GetBackupTrack()->GetChi2()/(track->GetBackupTrack()->GetNumberOfClusters()+1)) < 7)) {
+ seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
+ isGold = kTRUE;
+ }
+ }
+
+ if ((track->StatusForTOF() > 0) && (track->GetNCross() == 0) && (Float_t(track->GetNumberOfClusters()) / Float_t(track->GetNExpected()) > 0.4)) {
+ //seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
+ }
+ }
+ }
+ /**/
+
+ /**/
+ // 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 (track->GetBackupTrack()) {
+ cstream << "Tracks"
+ << "EventNrInFile=" << eventNrInFile
+ << "ESD.=" << seed
+ << "trd.=" << track
+ << "trdback.=" << track->GetBackupTrack()
+ << "\n";
+ }
+ else {
+ cstream << "Tracks"
+ << "EventNrInFile=" << eventNrInFile
+ << "ESD.=" << seed
+ << "trd.=" << track
+ << "trdback.=" << track
+ << "\n";
+ }
+ }
+ /**/
+
+ // Propagation to the TOF (I.Belikov)
+ if (track->GetStop() == kFALSE) {
+ fHBackfit->Fill(5);
+
+ Double_t xtof = 371.0;
+ Double_t xTOF0 = 370.0;
+
+ Double_t c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
+ if (TMath::Abs(c2) >= 0.99) {
+ fHBackfit->Fill(6);
+ delete track;
+ continue;
+ }
+
+ PropagateToX(*track,xTOF0,fgkMaxStep);
+
+ // Energy losses taken to the account - check one more time
+ c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
+ if (TMath::Abs(c2) >= 0.99) {
+ fHBackfit->Fill(7);
+ delete track;
+ continue;
+ }
+
+ //if (!PropagateToX(*track,xTOF0,fgkMaxStep)) {
+ // fHBackfit->Fill(7);
+ //delete track;
+ // continue;
+ //}
+
+ Double_t ymax = xtof * TMath::Tan(0.5 * AliTRDgeometry::GetAlpha());
+ Double_t y;
+ track->GetYAt(xtof,GetBz(),y);
+ if (y > ymax) {
+ if (!track->Rotate( AliTRDgeometry::GetAlpha())) {
+ fHBackfit->Fill(8);
+ delete track;
+ continue;
+ }
+ }
+ else if (y < -ymax) {
+ if (!track->Rotate(-AliTRDgeometry::GetAlpha())) {
+ fHBackfit->Fill(9);
+ delete track;
+ continue;
+ }
+ }
+
+ if (track->PropagateTo(xtof)) {
+ seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
+ fHBackfit->Fill(10);
+
+ for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
+ for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
+ seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
+ }
+ seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
+ }
+ //seed->SetTRDtrack(new AliTRDtrack(*track));
+ if (track->GetNumberOfClusters() > foundMin) {
+ fHBackfit->Fill(11);
+ found++;
+ }
+ }
+ }
+ else {
+ fHBackfit->Fill(12);
+
+ if ((track->GetNumberOfClusters() > 15) &&
+ (track->GetNumberOfClusters() > 0.5*expectedClr)) {
+ seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
+ fHBackfit->Fill(13);
+
+ //seed->SetStatus(AliESDtrack::kTRDStop);
+ for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
+ for (Int_t j = 0; j <AliESDtrack::kNSlice; j++) {
+ seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
+ }
+ seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
+ }
+ //seed->SetTRDtrack(new AliTRDtrack(*track));
+ found++;
+ }
+ }
+
+ seed->SetTRDQuality(track->StatusForTOF());
+ seed->SetTRDBudget(track->GetBudget(0));
+
+ fHBackfit->Fill(14);
+ delete track;
}
- }
-
- }
- else {
-
- fHBackfit->Fill(12);
-
- if ((track->GetNumberOfClusters() > 15) &&
- (track->GetNumberOfClusters() > 0.5*expectedClr)) {
- seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
- fHBackfit->Fill(13);
-
- //seed->SetStatus(AliESDtrack::kTRDStop);
- for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
- for (Int_t j = 0; j <AliESDtrack::kNSlice; j++) {
- seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
- }
- seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
- }
- //seed->SetTRDtrack(new AliTRDtrack(*track));
- found++;
- }
-
- }
-
- seed->SetTRDQuality(track->StatusForTOF());
- seed->SetTRDBudget(track->GetBudget(0));
-
- fHBackfit->Fill(14);
- delete track;
- }
-
- AliInfo(Form("Number of seeds: %d",fNseeds));
- AliInfo(Form("Number of back propagated TRD tracks: %d",found));
-
- // New seeding
- if (AliTRDReconstructor::SeedingOn()) {
- MakeSeedsMI(3,5,event);
- }
-
- fSeeds->Clear();
- fNseeds = 0;
-
- delete [] index;
- delete [] quality;
-
- SaveLogHists();
+ AliInfo(Form("Number of seeds: %d",fNseeds));
+ AliInfo(Form("Number of back propagated TRD tracks: %d",found));
+
+ fSeeds->Clear();
+ fNseeds = 0;
+
+ delete [] index;
+ delete [] quality;
+
+ SaveLogHists();
return 0;
-
}
//_____________________________________________________________________________
// Get material budget
Double_t xyz0[3];
Double_t xyz1[3];
- Double_t param[7];
+ Double_t param[7];
Double_t x;
Double_t y;
Double_t z;
fHClSearch->Fill(hb+6);
continue;
}
- //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
//
// Propagate and update track
//
for (Int_t itime = GetTimeBinsPerPlane()-1; itime >= 0; itime--) {
-
Int_t ilayer = GetGlobalTimeBin(0, iplane,itime);
expectedNumberOfClusters++;
t.SetNExpected(t.GetNExpected() + 1);
x = timeBin.GetX();
if (timeBin) {
-
- if (clusters[ilayer] > 0) {
- index = clusters[ilayer];
- cl = (AliTRDcluster *)GetCluster(index);
- //Double_t h01 = GetTiltFactor(cl); // I.B's fix
- //maxChi2=t.GetPredictedChi2(cl,h01); //
- }
+ if (clusters[ilayer] > 0) {
+ index = clusters[ilayer];
+ cl = (AliTRDcluster *)GetCluster(index);
+ //Double_t h01 = GetTiltFactor(cl); // I.B's fix
+ //maxChi2=t.GetPredictedChi2(cl,h01); //
+ }
if (cl) {
-
- //if (cl->GetNPads() < 5)
- Double_t dxsample = timeBin.GetdX();
- t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
- Double_t h01 = GetTiltFactor(cl);
- Int_t det = cl->GetDetector();
- Int_t plane = fGeom->GetPlane(det);
- if (t.GetX() > 345.0) {
- t.SetNLast(t.GetNLast() + 1);
- t.SetChi2Last(t.GetChi2Last() + maxChi2);
- }
- Double_t xcluster = cl->GetX();
- t.PropagateTo(xcluster,xx0,xrho);
- maxChi2 = t.GetPredictedChi2(cl,h01);
-
- if (maxChi2<1e+10)
- if (!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
- if (!t.Update(cl,maxChi2,index,h01)) {
- // ????
- }
- } // else SetCluster(cl, GetNumberOfClusters()-1); // A.Bercuci 25.07.07
+ //if (cl->GetNPads() < 5)
+ Double_t dxsample = timeBin.GetdX();
+ t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
+ Double_t h01 = GetTiltFactor(cl);
+ Int_t det = cl->GetDetector();
+ Int_t plane = fGeom->GetPlane(det);
+ if (t.GetX() > 345.0) {
+ t.SetNLast(t.GetNLast() + 1);
+ t.SetChi2Last(t.GetChi2Last() + maxChi2);
+ }
+ Double_t xcluster = cl->GetX();
+ t.PropagateTo(xcluster,xx0,xrho);
+ maxChi2 = t.GetPredictedChi2(cl,h01);
+
+ if (maxChi2<1e+10)
+ if (!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
+ if (!t.Update(cl,maxChi2,index,h01)) {
+ // ????
+ }
+ } // else SetCluster(cl, GetNumberOfClusters()-1); // A.Bercuci 25.07.07
if (calibra->GetMITracking()) {
calibra->UpdateHistograms(cl,&t);
}
- // Reset material budget if 2 consecutive gold
- if (plane > 0) {
- if ((t.GetTracklets(plane).GetN() + t.GetTracklets(plane-1).GetN()) > 20) {
- t.SetBudget(2,0.0);
- }
- }
-
- }
+ // Reset material budget if 2 consecutive gold
+ if (plane > 0) {
+ if ((t.GetTracklets(plane).GetN() + t.GetTracklets(plane-1).GetN()) > 20) {
+ t.SetBudget(2,0.0);
+ }
+ }
+
+ }
}
// differs from that of TRD sectors
//
- if (ReadClusters(fClusters,cTree)) {
+
+ if (ReadClusters(fClusters, cTree)) {
AliError("Problem with reading the clusters !");
return 1;
}
}
//_____________________________________________________________________________
-void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/, AliESDEvent *esd)
+Int_t AliTRDtracker::Clusters2Tracks(AliESDEvent *esd)
{
//
// Creates seeds using clusters between position inner plane and outer plane
delete [] pseed;
+ return 0;
}
//_____________________________________________________________________________
// Loop through all entries in the tree
Int_t nEntries = (Int_t) clusterTree->GetEntries();
Int_t nbytes = 0;
- AliTRDcluster *c = 0;
+ AliTRDcluster *c = 0x0;
for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
// Import the tree
// Get the number of points in the detector
Int_t nCluster = clusterArray->GetEntriesFast();
-
// Loop through all TRD digits
for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
- c = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster);
- AliTRDcluster *co = c;
- array->AddLast(co);
+ if(!(c = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster))) continue;
+ array->AddLast(c);
+ //printf("Add cluster 0x%x.\n", c);
clusterArray->RemoveAt(iCluster);
}
}
-
delete clusterArray;
return 0;
}
-//_____________________________________________________________________________
-AliTRDtracker::AliTRDpropagationLayer
- ::AliTRDpropagationLayer(Double_t x, Double_t dx, Double_t rho
- , Double_t radLength, Int_t tbIndex, Int_t plane)
- :fN(0)
- ,fSec(0)
- ,fClusters(NULL)
- ,fIndex(NULL)
- ,fX(x)
- ,fdX(dx)
- ,fRho(rho)
- ,fX0(radLength)
- ,fTimeBinIndex(tbIndex)
- ,fPlane(plane)
- ,fYmax(0)
- ,fYmaxSensitive(0)
- ,fHole(kFALSE)
- ,fHoleZc(0)
- ,fHoleZmax(0)
- ,fHoleYc(0)
- ,fHoleYmax(0)
- ,fHoleRho(0)
- ,fHoleX0(0)
-{
- //
- // AliTRDpropagationLayer constructor
- //
-
- for (Int_t i = 0; i < (Int_t) kZones; i++) {
- fZc[i] = 0;
- fZmax[i] = 0;
- }
-
- if (fTimeBinIndex >= 0) {
- fClusters = new AliTRDcluster*[kMaxClusterPerTimeBin];
- fIndex = new UInt_t[kMaxClusterPerTimeBin];
- }
-
- for (Int_t i = 0; i < 5; i++) {
- fIsHole[i] = kFALSE;
- }
-
-}
-
-//_____________________________________________________________________________
-void AliTRDtracker::AliTRDpropagationLayer
- ::SetHole(Double_t Zmax, Double_t Ymax, Double_t rho
- , Double_t radLength, Double_t Yc, Double_t Zc)
-{
- //
- // Sets hole in the layer
- //
-
- fHole = kTRUE;
- fHoleZc = Zc;
- fHoleZmax = Zmax;
- fHoleYc = Yc;
- fHoleYmax = Ymax;
- fHoleRho = rho;
- fHoleX0 = radLength;
-
-}
//_____________________________________________________________________________
AliTRDtracker::AliTRDtrackingSector
}
-//_____________________________________________________________________________
-void AliTRDtracker::AliTRDpropagationLayer
- ::SetZ(Double_t *center, Double_t *w, Double_t *wsensitive )
-{
- //
- // set centers and the width of sectors
- //
-
- for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
- fZc[icham] = center[icham];
- fZmax[icham] = w[icham];
- fZmaxSensitive[icham] = wsensitive[icham];
- }
-
-}
-
-//_____________________________________________________________________________
-void AliTRDtracker::AliTRDpropagationLayer::SetHoles(Bool_t *holes)
-{
- //
- // set centers and the width of sectors
- //
-
- fHole = kFALSE;
-
- for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
- fIsHole[icham] = holes[icham];
- if (holes[icham]) {
- fHole = kTRUE;
- }
- }
-
-}
-
-//_____________________________________________________________________________
-void AliTRDtracker::AliTRDpropagationLayer
- ::InsertCluster(AliTRDcluster *c, UInt_t index)
-{
- //
- // Insert cluster in cluster array.
- // Clusters are sorted according to Y coordinate.
- //
-
- if (fTimeBinIndex < 0) {
- //AliWarning("Attempt to insert cluster into non-sensitive time bin!\n");
- return;
- }
-
- if (fN == (Int_t) kMaxClusterPerTimeBin) {
- //AliWarning("Too many clusters !\n");
- return;
- }
-
- if (fN == 0) {
- fIndex[0] = index;
- fClusters[fN++] = c;
- return;
- }
-
- Int_t i = Find(c->GetY());
- memmove(fClusters+i+1,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
- memmove(fIndex +i+1,fIndex +i,(fN-i)*sizeof(UInt_t));
- fIndex[i] = index;
- fClusters[i] = c;
- fN++;
-
-}
-
-//_____________________________________________________________________________
-Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Float_t y) const
-{
- //
- // Returns index of the cluster nearest in Y
- //
-
- if (fN <= 0) {
- return 0;
- }
- if (y <= fClusters[0]->GetY()) {
- return 0;
- }
- if (y > fClusters[fN-1]->GetY()) {
- return fN;
- }
-
- Int_t b = 0;
- Int_t e = fN - 1;
- Int_t m = (b + e) / 2;
-
- for ( ; b < e; m = (b + e) / 2) {
- if (y > fClusters[m]->GetY()) {
- b = m + 1;
- }
- else {
- e = m;
- }
- }
-
- return m;
-
-}
-
-//_____________________________________________________________________________
-Int_t AliTRDtracker::AliTRDpropagationLayer
- ::FindNearestCluster(Float_t y, Float_t z, Float_t maxroad
- , Float_t maxroadz) const
-{
- //
- // Returns index of the cluster nearest to the given y,z
- //
-
- Int_t index = -1;
- Int_t maxn = fN;
- Float_t mindist = maxroad;
-
- for (Int_t i = Find(y-maxroad); i < maxn; i++) {
- AliTRDcluster *c = (AliTRDcluster *) (fClusters[i]);
- Float_t ycl = c->GetY();
- if (ycl > (y + maxroad)) {
- break;
- }
- if (TMath::Abs(c->GetZ() - z) > maxroadz) {
- continue;
- }
- if (TMath::Abs(ycl - y) < mindist) {
- mindist = TMath::Abs(ycl - y);
- index = fIndex[i];
- }
- }
-
- return index;
-
-}
-
//_____________________________________________________________________________
Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster *c)
{
}
-
//_____________________________________________________________________________
Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1
, AliTRDtrack *track
AliTRDtrack *AliTRDtracker::RegisterSeed(AliTRDseed *seeds, Double_t *params)
{
//
- // Register a seed
+ // Build a TRD track out of tracklet candidates
+ //
+ // Parameters :
+ // seeds : array of tracklets
+ // params : track parameters (see MakeSeeds() function body for a detailed description)
+ //
+ // Output :
+ // The TRD track.
+ //
+ // Detailed description
+ //
+ // To be discussed with Marian !!
//
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ Int_t nTimeBins = cal->GetNumberOfTimeBins();
+
+
Double_t alpha = AliTRDgeometry::GetAlpha();
Double_t shift = AliTRDgeometry::GetAlpha()/2.0;
Double_t c[15];
AliTRDcluster *cl = 0;
for (Int_t ilayer = 0; ilayer < 6; ilayer++) {
- if (seeds[ilayer].IsOK()) {
- for (Int_t itime = 22; itime > 0; itime--) {
- if (seeds[ilayer].GetIndexes(itime) > 0) {
- index = seeds[ilayer].GetIndexes(itime);
- cl = seeds[ilayer].GetClusters(itime);
- break;
- }
+
+ if (seeds[ilayer].IsOK()) {
+ for (Int_t itime = nTimeBins-1; itime > 0; itime--) {
+ if (seeds[ilayer].GetIndexes(itime) > 0) {
+ index = seeds[ilayer].GetIndexes(itime);
+ cl = seeds[ilayer].GetClusters(itime);
+ //printf("l[%d] index[%d] tb[%d] cptr[%p]\n", ilayer, index, itime, cl);
+ break;
+ }
}
}
if (index > 0) {
// //
////////////////////////////////////////////////////////////////////////////
+#ifndef ROOT_TObjArray
#include "TObjArray.h"
+#endif
+#ifndef ALITRACKER_H
#include "AliTracker.h"
+#endif
+
+#ifndef ALITRDPROPAGATIONLAYER_H
+#include "AliTRDpropagationLayer.h"
+#endif
class TFile;
class TTree;
class AliTRDcluster;
class AliTRDseed;
class AliESDEvent;
+class AliTRDpropagationLayer;
///////////////////////////////////////////////////////////////////////////////
// //
enum { kMaxLayersPerSector = 1000
, kMaxTimeBinIndex = 216
- , kMaxClusterPerTimeBin = 2300
- , kZones = 5
, kTrackingSectors = 18 };
AliTRDtracker();
Int_t ReadClusters(TObjArray *array, TTree *in) const;
AliTRDcluster *GetCluster(AliTRDtrack *track, Int_t plane, Int_t timebin, UInt_t &index);
Int_t FindClusters(Int_t sector, Int_t t0, Int_t t1, AliTRDtrack *track
- , Int_t *clusters, AliTRDtracklet &tracklet);
+ , Int_t *clusters, AliTRDtracklet &tracklet);
protected:
-
- //__________________________________________________________________________________________________
- class AliTRDpropagationLayer {
-
- public:
-
- AliTRDpropagationLayer(Double_t x, Double_t dx, Double_t rho
- , Double_t x0, Int_t tbIndex, Int_t plane);
- AliTRDpropagationLayer(const AliTRDpropagationLayer &/*p*/);
- ~AliTRDpropagationLayer() { if (fTimeBinIndex >= 0) {
- delete[] fClusters;
- delete[] fIndex;
- }
- }
- AliTRDpropagationLayer &operator=(const AliTRDpropagationLayer &/*p*/)
- { return *this; }
-
- operator Int_t() const { return fN; }
- AliTRDcluster *operator[](Int_t i) { return fClusters[i]; }
-
- void SetZmax(Int_t cham, Double_t center, Double_t w) { fZc[cham] = center;
- fZmax[cham] = w; }
- void SetYmax(Double_t w, Double_t wsensitive) { fYmax = w;
- fYmaxSensitive = wsensitive; }
- void SetZ(Double_t* center, Double_t *w, Double_t *wsensitive);
- void SetHoles(Bool_t* holes);
- void SetHole(Double_t Zmax, Double_t Ymax
- , Double_t rho = 1.29e-3, Double_t x0 = 36.66
- , Double_t Yc = 0.0, Double_t Zc = 0.0);
-
- Double_t GetYmax() const { return fYmax; }
- Double_t GetZmax(Int_t c) const { return fZmax[c]; }
- Double_t GetZc(Int_t c) const { return fZc[c]; }
- UInt_t GetIndex(Int_t i) const { return fIndex[i]; }
- Double_t GetX() const { return fX; }
- Double_t GetdX() const { return fdX; }
- Int_t GetTimeBinIndex() const { return fTimeBinIndex; }
- Int_t GetPlane() const { return fPlane; }
- Bool_t IsHole(Int_t zone) const { return fIsHole[zone]; }
- Bool_t IsSensitive() const { return (fTimeBinIndex >= 0) ? kTRUE : kFALSE;}
-
- void Clear() {
- for (Int_t i = 0; i < fN; i++) fClusters[i] = NULL;
- fN = 0;
- }
-
- void InsertCluster(AliTRDcluster *c, UInt_t index);
- Int_t Find(Float_t y) const;
- Int_t FindNearestCluster(Float_t y, Float_t z, Float_t maxroad, Float_t maxroadz) const;
-
- void SetX(Double_t x) { fX = x; }
-
- private:
-
- Int_t fN; // This is fN
- Int_t fSec; // Sector mumber
- AliTRDcluster **fClusters; // Array of pointers to clusters
- UInt_t *fIndex; // Array of cluster indexes
- Double_t fX; // X coordinate of the middle plane
- Double_t fdX; // Radial thickness of the time bin
- Double_t fRho; // Default density of the material
- Double_t fX0; // Default radiation length
- Int_t fTimeBinIndex; // Plane * F(local_tb)
- Int_t fPlane; // Plane number
- Double_t fZc[kZones]; // Z position of the center for 5 active areas
- Double_t fZmax[kZones]; // Half of active area length in Z
- Double_t fZmaxSensitive[kZones]; // Sensitive area for detection Z
- Bool_t fIsHole[kZones]; // Is hole in given sector
- Double_t fYmax; // Half of active area length in Y
- Double_t fYmaxSensitive; // Half of active area length in Y
-
- Bool_t fHole; // kTRUE if there is a hole in the layer
- Double_t fHoleZc; // Z of the center of the hole
- Double_t fHoleZmax; // Half of the hole length in Z
- Double_t fHoleYc; // Y of the center of the hole
- Double_t fHoleYmax; // Half of the hole length in Y
- Double_t fHoleRho; // Density of the gas in the hole
- Double_t fHoleX0; // Radiation length of the gas in the hole
-
- };
+ Bool_t AdjustSector(AliTRDtrack *track);
+ AliTRDtrack *RegisterSeed(AliTRDseed *seeds, Double_t *params);
+ Int_t FollowBackProlongation(AliTRDtrack &t);
+ //void MakeSeedsMI(Int_t inner, Int_t outer, AliESDEvent *esd = 0);
+ protected:
+
//__________________________________________________________________________________________________
class AliTRDtrackingSector {
Int_t GetLayerNumber(Int_t tb) const { return fTimeBinIndex[tb]; }
Double_t GetX(Int_t pl) const { return fLayers[pl]->GetX(); }
AliTRDpropagationLayer* GetLayer(Int_t i) { return fLayers[i]; }
+ Int_t GetSector() const {return fGeomSector;}
void MapTimeBinLayers();
Int_t Find(Double_t x) const;
};
+ protected:
+
AliTRDgeometry *fGeom; // Pointer to TRD geometry
AliTRDtrackingSector *fTrSec[kTrackingSectors]; // Array of tracking sectors;
Int_t fNclusters; // Number of clusters in TRD
TH2D *fHMinD; // QA histogram
TH1D *fHDeltaX; // QA histogram
TH1D *fHXCl; // QA histogram
-
- Bool_t AdjustSector(AliTRDtrack *track);
-
+
private:
-
- AliTRDtrack *RegisterSeed(AliTRDseed *seeds, Double_t *params);
- void MakeSeedsMI(Int_t inner, Int_t outer, AliESDEvent *esd = 0);
-
- Int_t FollowBackProlongation(AliTRDtrack &t);
- Int_t FollowProlongation(AliTRDtrack &t);
- Int_t PropagateToX(AliTRDtrack &t, Double_t xToGo, Double_t maxStep);
- Double_t ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt) const;
- Double_t ExpectedSigmaZ2(Double_t r, Double_t tgl) const;
-
+
+ Int_t FollowProlongation(AliTRDtrack &t);
+ Int_t PropagateToX(AliTRDtrack &t, Double_t xToGo, Double_t maxStep);
+ Double_t ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt) const;
+ Double_t ExpectedSigmaZ2(Double_t r, Double_t tgl) const;
+
+ private:
+
TTreeSRedirector *fDebugStreamer; //!Debug streamer
- ClassDef(AliTRDtracker,3) // TRD tracker
+ ClassDef(AliTRDtracker, 4) // TRD tracker
};
+
+
#endif
--- /dev/null
+/**************************************************************************
+ * 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$ */
+
+///////////////////////////////////////////////////////////////////////////////
+// //
+// Track fitter //
+// //
+// Authors: //
+// Alex Bercuci <A.Bercuci@gsi.de> //
+// Markus Fasel <M.Fasel@gsi.de> //
+// //
+///////////////////////////////////////////////////////////////////////////////
+
+#include <TROOT.h>
+#include <TLinearFitter.h>
+#include <TMath.h>
+#include <TTreeStream.h>
+
+#include "AliLog.h"
+
+#include "AliTRDseed.h"
+#include "AliTRDcalibDB.h"
+#include "AliTRDcluster.h"
+#include "AliTRDReconstructor.h"
+#include "AliTRDseedV1.h"
+#include "AliTRDtrackerFitter.h"
+
+#define DEBUG
+
+ClassImp(AliTRDtrackerFitter)
+
+//_________________________________________________________________________
+AliTRDtrackerFitter::AliTRDtrackerFitter()
+ :TObject()
+ ,fFitterTC(0x0)
+ ,fFitterT2(0x0)
+ ,fRieman1(0x0)
+ ,fRieman2(0x0)
+ ,fChi2TR(0)
+ ,fChi2TC(0)
+ ,fCR(0)
+ ,fCC(0)
+ ,fDca(0)
+ ,fDzmf(0)
+ ,fZmf(0)
+ ,fNlayers(0)
+ ,fDebugStream(0x0)
+{
+ //
+ // Constructor
+ //
+
+ fRieman1 = new AliRieman(1000);
+ fRieman2 = new AliRieman(1000);
+ fFitterTC = new TLinearFitter(2,"hyp2");
+ fFitterT2 = new TLinearFitter(4,"hyp4");
+ fFitterTC->StoreData(kTRUE);
+ fFitterT2->StoreData(kTRUE);
+}
+
+//_________________________________________________________________________
+AliTRDtrackerFitter::AliTRDtrackerFitter(const AliTRDtrackerFitter &f)
+ :TObject(f)
+ ,fFitterTC(0x0)
+ ,fFitterT2(0x0)
+ ,fRieman1(0x0)
+ ,fRieman2(0x0)
+ ,fChi2TR(f.fChi2TR)
+ ,fChi2TC(f.fChi2TC)
+ ,fCR(f.fCR)
+ ,fCC(f.fCC)
+ ,fDca(f.fDca)
+ ,fDzmf(f.fDzmf)
+ ,fZmf(f.fZmf)
+ ,fNlayers(f.fNlayers)
+ ,fDebugStream(0x0)
+{
+ //
+ // Copy Constructor (performs a deep copy)
+ //
+
+ fRieman1 = new AliRieman(*f.fRieman1);
+ fRieman2 = new AliRieman(*f.fRieman2);
+ fFitterTC = new TLinearFitter(*f.fFitterTC);
+ fFitterT2 = new TLinearFitter(*f.fFitterT2);
+ fFitterTC->StoreData(kTRUE);
+ fFitterT2->StoreData(kTRUE);
+}
+
+//_________________________________________________________________________
+AliTRDtrackerFitter::~AliTRDtrackerFitter()
+{
+ //
+ // Destructor
+ //
+
+ delete fRieman1;
+ delete fRieman2;
+ delete fFitterTC;
+ delete fFitterT2;
+}
+
+//_________________________________________________________________________
+AliTRDtrackerFitter &AliTRDtrackerFitter::operator=(const AliTRDtrackerFitter& fitter)
+{
+ //
+ // Assignment operator using the copy Method of this class
+ //
+ if(this != &fitter){
+ fitter.Copy(*this);
+ }
+ return *this;
+}
+
+//_________________________________________________________________________
+void AliTRDtrackerFitter::Copy(TObject &f) const
+{
+ //
+ // Copy method.
+ // Performs a deep copy. Mainly used in the Assignment operator
+ //
+
+ AliTRDtrackerFitter &tf = (AliTRDtrackerFitter &) f;
+
+ if(tf.fFitterTC)
+ delete tf.fFitterTC;
+ tf.fFitterTC = new TLinearFitter(*fFitterTC);
+ if(tf.fFitterT2)
+ delete tf.fFitterT2;
+ tf.fFitterT2 = new TLinearFitter(*fFitterT2);
+ if(tf.fRieman1)
+ delete tf.fRieman1;
+ tf.fRieman1 = new AliRieman(*fRieman1);
+ if(tf.fRieman2)
+ delete fRieman2;
+ tf.fRieman2 = new AliRieman(*fRieman2);
+ tf.fChi2TR = fChi2TR;
+ tf.fChi2TC = fChi2TC;
+ tf.fCR = fCR;
+ tf.fCC = fCC;
+ tf.fDca =fDca;
+ tf.fDzmf = fDzmf;
+ tf.fZmf =fZmf;
+ tf.fNlayers = fNlayers;
+ tf.fDebugStream = 0x0;
+ fFitterTC->StoreData(kTRUE);
+ fFitterT2->StoreData(kTRUE);
+}
+
+//_________________________________________________________________________
+void AliTRDtrackerFitter::FitRieman(AliTRDcluster **cl, Int_t nLayers)
+{
+ //
+ // Performs a Rieman Fit including 4 clusters (one in each layer)
+ //
+ fRieman1->Reset();
+ for(Int_t i = 0; i < nLayers; i++)
+ fRieman1->AddPoint(cl[i]->GetX(), cl[i]->GetY(), cl[i]->GetZ(), 1, 10);
+ fRieman1->Update();
+}
+
+//_________________________________________________________________________
+void AliTRDtrackerFitter::FitRieman(AliTRDseedV1 *cseed, Int_t *planes)
+{
+ //
+ // Performs a Rieman fit using four respectively six seeds
+ // 2 times used: 1st with 4 parameters and Later in the full track
+ // fitting Part with 6 parameters
+ //
+
+ Int_t lplanes[] = {0, 1, 2, 3, 4, 5}, *tplanes;
+ Int_t nplanes;
+ if(planes){
+ nplanes = 4;
+ tplanes = planes;
+ } else {
+ nplanes = 6;
+ tplanes = &lplanes[0];
+ }
+ fRieman1->Reset();
+ for(Int_t iLayer = 0; iLayer < nplanes; iLayer++){
+ if(!cseed[tplanes[iLayer]].IsOK()) continue;
+ fRieman1->AddPoint(cseed[tplanes[iLayer]].GetX0(), cseed[tplanes[iLayer]].GetYfitR(0), cseed[tplanes[iLayer]].GetZProb(), 1, 10);
+ }
+ fRieman1->Update();
+}
+
+//_________________________________________________________________________
+Double_t AliTRDtrackerFitter::FitHyperplane(AliTRDseedV1 *cseed
+ , Double_t chi2ZF
+ , Double_t zFitter)
+{
+ //
+ // performs a hyperplane fit
+ //
+ // Two linear fitters: one in which kz is fixed and one in which kz is a free parameter
+ // checking if values are acceptable and improves the fitresults
+ // Calculates curvature parameters and chisquares
+ // Return: likelihood value as a measurement for the seedquality
+ //
+
+ //const Int_t kMaxTimebins = 100; // Buffer
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ Int_t nTimeBins = cal->GetNumberOfTimeBins();
+ Int_t npointsT = 0;
+ fFitterTC->ClearPoints();
+ fFitterT2->ClearPoints();
+
+ Double_t fHl[6];
+ Double_t xref2 = (cseed[2].GetX0() + cseed[3].GetX0())/2;
+ fRieman2->Reset();
+ for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+ if (!cseed[iLayer].IsOK()) continue;
+ fHl[iLayer] = cseed[iLayer].GetTilt();
+
+ for (Int_t itime = 0; itime < nTimeBins; itime++) {
+ if (!cseed[iLayer].IsUsable(itime)) continue;
+
+ // X relative to the middle chamber
+ Double_t x = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0() - xref2;
+ Double_t y = cseed[iLayer].GetY(itime);
+ Double_t z = cseed[iLayer].GetZ(itime);
+ // ExB correction to the correction
+ // Tilted fRieman1
+ Double_t uvt[6];
+ // Global x
+ Double_t x2 = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0();
+ Double_t t = 1.0 / (x2*x2 + y*y);
+ uvt[1] = t; // t
+ uvt[0] = 2.0 * x2 * uvt[1]; // u
+ uvt[2] = 2.0 * fHl[iLayer] * uvt[1];
+ uvt[3] = 2.0 * fHl[iLayer] * x * uvt[1];
+ uvt[4] = 2.0 * (y + fHl[iLayer]*z) * uvt[1];
+ Double_t error = 2.0 * 0.2 * uvt[1];
+ fFitterT2->AddPoint(uvt,uvt[4],error);
+
+ //Constrained fRieman1
+ z = cseed[iLayer].GetZ(itime);
+ uvt[0] = 2.0 * x2 * t; // u
+ uvt[1] = 2.0 * fHl[iLayer] * x2 * uvt[1];
+ uvt[2] = 2.0 * (y + fHl[iLayer] * (z - zFitter)) * t; // parameter that is most propably wrong
+ fFitterTC->AddPoint(uvt,uvt[2],error);
+ fRieman2->AddPoint(x2,y,z,1,10);
+ npointsT++;
+ }
+ } // Loop: iLayer
+
+ fRieman2->Update();
+ fFitterTC->Eval();
+ fFitterT2->Eval();
+ Double_t rpolz0 = fFitterT2->GetParameter(3);
+ Double_t rpolz1 = fFitterT2->GetParameter(4);
+
+ // Linear fitter - not possible to make boundaries
+ // Do not accept non possible z and dzdx combinations
+ Bool_t acceptablez = kTRUE;
+ for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+ if(!cseed[iLayer].IsOK()) continue;
+ Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].GetX0() - xref2);
+ if (TMath::Abs(cseed[iLayer].GetZProb() - zT2) > cseed[iLayer].GetPadLength() * 0.5 + 1.0) acceptablez = kFALSE;
+ }
+ if (!acceptablez) {
+ fFitterT2->FixParameter(3,fRieman1->GetZat(xref2));
+ fFitterT2->FixParameter(4,fRieman1->GetDZat(xref2));
+ fFitterT2->Eval();
+ fFitterT2->ReleaseParameter(3);
+ fFitterT2->ReleaseParameter(4);
+ rpolz0 = fFitterT2->GetParameter(3);
+ rpolz1 = fFitterT2->GetParameter(4);
+ }
+
+ fChi2TR = fFitterT2->GetChisquare() / Float_t(npointsT);
+ fChi2TC = fFitterTC->GetChisquare() / Float_t(npointsT);
+ Double_t polz1c = fFitterTC->GetParameter(2);
+ Double_t polz0c = polz1c * xref2;
+ Double_t aC = fFitterTC->GetParameter(0);
+ Double_t bC = fFitterTC->GetParameter(1);
+ fCC = aC / TMath::Sqrt(bC * bC + 1.0); // Curvature
+ Double_t aR = fFitterT2->GetParameter(0);
+ Double_t bR = fFitterT2->GetParameter(1);
+ Double_t dR = fFitterT2->GetParameter(2);
+ fCR = 1.0 + bR*bR - dR*aR;
+ fDca = 0.0;
+ if (fCR > 0.0) {
+ fDca = -dR / (TMath::Sqrt(1.0 + bR*bR - dR*aR) + TMath::Sqrt(1.0 + bR*bR));
+ fCR = aR / TMath::Sqrt(fCR);
+ }
+
+
+ Double_t chi2ZT2 = 0.0;
+ Double_t chi2ZTC = 0.0;
+ for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+ if(!cseed[iLayer].IsOK()) continue;
+ Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].GetX0() - xref2);
+ Double_t zTC = polz0c + polz1c * (cseed[iLayer].GetX0() - xref2);
+ chi2ZT2 += TMath::Abs(cseed[iLayer].GetMeanz() - zT2);
+ chi2ZTC += TMath::Abs(cseed[iLayer].GetMeanz() - zTC);
+ }
+ chi2ZT2 /= TMath::Max((fNlayers - 3.0),1.0);
+ chi2ZTC /= TMath::Max((fNlayers - 3.0),1.0);
+
+
+ AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
+ Float_t sumdaf = 0.0;
+ for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+ if(!cseed[iLayer].IsOK()) continue;
+ sumdaf += TMath::Abs((cseed[iLayer].GetYfit(1) - cseed[iLayer].GetYref(1))/ cseed[iLayer].GetSigmaY2());
+ }
+ sumdaf /= Float_t (fNlayers - 2.0);
+
+ // Cook Likelihoods
+ Double_t likezf = TMath::Exp(-chi2ZF * 0.14);
+ Double_t likechi2C = TMath::Exp(-fChi2TC * 0.677);
+ Double_t likechi2TR = TMath::Exp(-fChi2TR * 0.78);
+ Double_t likeaf = TMath::Exp(-sumdaf * 3.23);
+ Double_t likef = likezf * likechi2TR * likeaf;
+
+#ifdef DEBUG
+ if(fDebugStream && AliTRDReconstructor::StreamLevel() >= 2){
+ TTreeSRedirector &treeStreamer = *fDebugStream;
+ treeStreamer << "FitHyperplane"
+ << "seed0.=" << &cseed[0]
+ << "seed1.=" << &cseed[1]
+ << "seed2.=" << &cseed[2]
+ << "seed3.=" << &cseed[3]
+ << "seed4.=" << &cseed[4]
+ << "seed5.=" << &cseed[5]
+ << "chi2TR=" << fChi2TR
+ << "chi2TC=" << fChi2TC
+ << "chi2ZT2=" << chi2ZT2
+ << "chi2ZTC=" << chi2ZTC
+ << "CC=" << fCC
+ << "CR=" << fCR
+ << "DR=" << dR
+ << "DCA=" << fDca
+ << "Polz0=" << polz0c
+ << "Polz1=" << polz1c
+ << "RPolz0=" << rpolz0
+ << "RPolz1=" << rpolz1
+ << "Likechi2C=" << likechi2C
+ << "Likechi2TR=" << likechi2TR
+ << "Likezf=" << likezf
+ << "Likeaf=" << likeaf
+ << "LikeF=" << likef
+ << "Rieman1.=" << fRieman1
+ << "Rieman2.=" << fRieman2
+ << "\n";
+ }
+#endif
+
+ return likef;
+}
+
+//_________________________________________________________________________
+void AliTRDtrackerFitter::GetHyperplaneFitChi2(Double_t *chisquares) const
+{
+ //
+ // Getter method returns the chisquares of the hyperplane fit
+ //
+
+ chisquares[0] = fChi2TR;
+ chisquares[1] = fChi2TC;
+}
+
+//_________________________________________________________________________
+void AliTRDtrackerFitter::GetHyperplaneFitResults(Double_t *params) const
+{
+ //
+ // Getter Method returning the Curvature parameters of the hyperplane fit
+ //
+
+ params[0] = fCC;
+ params[1] = fCR;
+ params[2] = fDca;
+}
+
+//_________________________________________________________________________
+void AliTRDtrackerFitter::Reset()
+{
+ //
+ // Resets the object
+ //
+
+ fRieman1->Reset();
+ fRieman2->Reset();
+ //fFitterTC->Clear();
+ //fFitterT2->Clear();
+}
--- /dev/null
+#ifndef ALITRDTRACKERFITTER_H
+#define ALITRDTRACKERFITTER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+
+///////////////////////////////////////////////////////////////////////////////
+// //
+// Track fitter //
+// //
+// Authors: //
+// Alex Bercuci <A.Bercuci@gsi.de> //
+// Markus Fasel <M.Fasel@gsi.de> //
+// //
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef AliRIEMAN_H
+#include <AliRieman.h>
+#endif
+
+#ifndef ROOT_Rtypes
+#include "Rtypes.h"
+#endif
+
+// implementations in the future
+class TObject;
+class TLinearFitter;
+class TTreeSRedirector;
+
+class AliTRDcluster;
+class AliTRDseedV1;
+
+class AliTRDtrackerFitter : public TObject {
+
+ public:
+
+ AliTRDtrackerFitter();
+ AliTRDtrackerFitter(const AliTRDtrackerFitter &);
+ AliTRDtrackerFitter& operator=(const AliTRDtrackerFitter &);
+ virtual ~AliTRDtrackerFitter();
+
+ void Copy(TObject &f) const;
+ AliRieman *GetRiemanFitter() const { return fRieman1; }
+ void GetHyperplaneFitChi2(Double_t *chisquares) const;
+ void GetHyperplaneFitResults(Double_t *params) const;
+ Double_t GetRiemanCurvature() const { return fRieman1->GetC(); }
+
+ void FitRieman(AliTRDcluster **cl, Int_t nLayers);
+ void FitRieman(AliTRDseedV1 *cseed, Int_t *planes=0x0);
+ Double_t FitHyperplane(AliTRDseedV1 *cseed, Double_t chi2ZF, Double_t zval);
+ void Reset();
+ void SetDebugStream(TTreeSRedirector *debug){fDebugStream = debug;}
+ void SetLayers(const Int_t nlayers) {fNlayers = nlayers;}
+
+ private:
+
+ // Linear fitters in planes
+ TLinearFitter *fFitterTC; // Fitting with tilting pads - kz fixed - kz= Z/x, + vertex const
+ TLinearFitter *fFitterT2; // Fitting with tilting pads - kz not fixed
+ AliRieman *fRieman1; // Rieman fitter
+ AliRieman *fRieman2; // Rieman fitter
+ // Linear Fitter chisquares
+ Double_t fChi2TR; // Chisquared
+ Double_t fChi2TC; // Chisquared
+ // Hyperplane Fit results
+ Double_t fCR; // Track radius
+ Double_t fCC; // Track curvature
+ Double_t fDca; // DCA
+ // Helping Parameters
+ Double_t fDzmf; // Something
+ Double_t fZmf; // Something else
+ Int_t fNlayers; // Some layers
+
+ // Debug
+ TTreeSRedirector *fDebugStream; //! debugger
+
+ ClassDef(AliTRDtrackerFitter,1) // TRD track fitter
+};
+#endif // ALITRDTRACKERFITTER_H
--- /dev/null
+/**************************************************************************
+ * 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$ */
+
+///////////////////////////////////////////////////////////////////////////////
+// //
+// Track finder //
+// //
+// Authors: //
+// Alex Bercuci <A.Bercuci@gsi.de> //
+// Markus Fasel <M.Fasel@gsi.de> //
+// //
+///////////////////////////////////////////////////////////////////////////////
+
+#include <Riostream.h>
+#include <stdio.h>
+#include <string.h>
+
+#include <TBranch.h>
+#include <TFile.h>
+#include <TGraph.h>
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TLinearFitter.h>
+#include <TObjArray.h>
+#include <TROOT.h>
+#include <TTree.h>
+#include <TClonesArray.h>
+#include <TRandom.h>
+#include <TTreeStream.h>
+
+#include "AliLog.h"
+#include "AliESDEvent.h"
+#include "AliAlignObj.h"
+#include "AliRieman.h"
+#include "AliTrackPointArray.h"
+
+#include "AliTRDtracker.h"
+#include "AliTRDtrackerV1.h"
+#include "AliTRDgeometry.h"
+#include "AliTRDpadPlane.h"
+#include "AliTRDgeometry.h"
+#include "AliTRDcluster.h"
+#include "AliTRDtrack.h"
+#include "AliTRDseed.h"
+#include "AliTRDcalibDB.h"
+#include "AliTRDCommonParam.h"
+#include "AliTRDReconstructor.h"
+#include "AliTRDCalibraFillHisto.h"
+#include "AliTRDtrackerFitter.h"
+#include "AliTRDstackLayer.h"
+#include "AliTRDrecoParam.h"
+#include "AliTRDseedV1.h"
+
+#define DEBUG
+
+ClassImp(AliTRDtrackerV1)
+Double_t AliTRDtrackerV1::fTopologicQA[kNConfigs] = {
+ 0.1112, 0.1112, 0.1112, 0.0786, 0.0786,
+ 0.0786, 0.0786, 0.0579, 0.0579, 0.0474,
+ 0.0474, 0.0408, 0.0335, 0.0335, 0.0335
+};
+
+//____________________________________________________________________
+AliTRDtrackerV1::AliTRDtrackerV1(AliTRDrecoParam *p)
+ :AliTRDtracker()
+ ,fSieveSeeding(0)
+ ,fRecoParam(p)
+ ,fFitter(0x0)
+ ,fDebugStreamer(0x0)
+{
+ //
+ // Default constructor. Nothing is initialized.
+ //
+
+}
+
+//____________________________________________________________________
+AliTRDtrackerV1::AliTRDtrackerV1(const TFile *in, AliTRDrecoParam *p)
+ :AliTRDtracker(in)
+ ,fSieveSeeding(0)
+ ,fRecoParam(p)
+ ,fFitter(0x0)
+ ,fDebugStreamer(0x0)
+{
+ //
+ // Standard constructor.
+ // Setting of the geometry file, debug output (if enabled)
+ // and initilize fitter helper.
+ //
+
+ fFitter = new AliTRDtrackerFitter();
+
+#ifdef DEBUG
+ fDebugStreamer = new TTreeSRedirector("TRDdebug.root");
+ fFitter->SetDebugStream(fDebugStreamer);
+#endif
+
+}
+
+//____________________________________________________________________
+AliTRDtrackerV1::~AliTRDtrackerV1()
+{
+ //
+ // Destructor
+ //
+
+ if(fDebugStreamer) delete fDebugStreamer;
+ if(fFitter) delete fFitter;
+ if(fRecoParam) delete fRecoParam;
+
+}
+
+//____________________________________________________________________
+Int_t AliTRDtrackerV1::Clusters2Tracks(AliESDEvent *esd)
+{
+ //
+ // Steering stand alone tracking for full TRD detector
+ //
+ // Parameters :
+ // esd : The ESD event. On output it contains
+ // the ESD tracks found in TRD.
+ //
+ // Output :
+ // Number of tracks found in the TRD detector.
+ //
+ // Detailed description
+ // 1. Launch individual SM trackers.
+ // See AliTRDtrackerV1::Clusters2TracksSM() for details.
+ //
+
+ if(!fRecoParam){
+ AliError("Reconstruction configuration not initialized. Call first AliTRDtrackerV1::SetRecoParam().");
+ return 0;
+ }
+
+ //AliInfo("Start Track Finder ...");
+ Int_t ntracks = 0;
+ for(int ism=0; ism<AliTRDtracker::kTrackingSectors; ism++){
+ //AliInfo(Form("Processing supermodule %i ...", ism));
+ ntracks += Clusters2TracksSM(fTrSec[ism], esd);
+ }
+ AliInfo(Form("Found %d TRD tracks.", ntracks));
+ return ntracks;
+}
+
+//____________________________________________________________________
+Int_t AliTRDtrackerV1::Clusters2TracksSM(AliTRDtracker::AliTRDtrackingSector *sector
+ , AliESDEvent *esd)
+{
+ //
+ // Steer tracking for one SM.
+ //
+ // Parameters :
+ // sector : Array of (SM) propagation layers containing clusters
+ // esd : The current ESD event. On output it contains the also
+ // the ESD (TRD) tracks found in this SM.
+ //
+ // Output :
+ // Number of tracks found in this TRD supermodule.
+ //
+ // Detailed description
+ //
+ // 1. Unpack AliTRDpropagationLayers objects for each stack.
+ // 2. Launch stack tracking.
+ // See AliTRDtrackerV1::Clusters2TracksStack() for details.
+ // 3. Pack results in the ESD event.
+ //
+
+ AliTRDpadPlane *pp = 0x0;
+
+ // allocate space for esd tracks in this SM
+ TClonesArray esdTrackList("AliESDtrack", 2*kMaxTracksStack);
+ esdTrackList.SetOwner();
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ Int_t nTimeBins = cal->GetNumberOfTimeBins();
+ const Int_t kFindable = Int_t(fRecoParam->GetFindableClusters()*6.*nTimeBins);
+
+ Int_t ntracks = 0;
+ Int_t nClStack = 0;
+ for(int istack = 0; istack<AliTRDpropagationLayer::kZones; istack++){
+ AliTRDstackLayer stackLayer[kNPlanes*kNTimeBins];
+
+ nClStack = 0;
+ //AliInfo(Form("Processing stack %i ...",istack));
+ //AliInfo("Building stack propagation layers ...");
+ for(int ilayer=0; ilayer<kNPlanes*nTimeBins; ilayer++){
+ pp = fGeom->GetPadPlane((Int_t)(ilayer/nTimeBins), istack);
+ Double_t stacklength = (pp->GetNrows() - 2) * pp->GetLengthIPad()
+ + 2 * pp->GetLengthOPad() + 2 * pp->GetLengthRim();
+ //Debug
+ Double_t z0 = fGeom->GetRow0((Int_t)(ilayer/nTimeBins),istack,0);
+ const AliTRDpropagationLayer smLayer(*(sector->GetLayer(ilayer)));
+ stackLayer[ilayer] = smLayer;
+#ifdef DEBUG
+ stackLayer[ilayer].SetDebugStream(fDebugStreamer);
+#endif
+ stackLayer[ilayer].SetRange(z0 - stacklength, stacklength);
+ stackLayer[ilayer].SetSector(sector->GetSector());
+ stackLayer[ilayer].SetStackNr(istack);
+ stackLayer[ilayer].SetNRows(pp->GetNrows());
+ stackLayer[ilayer].SetRecoParam(fRecoParam);
+ stackLayer[ilayer].BuildIndices();
+ nClStack += stackLayer[ilayer].GetNClusters();
+ }
+ //AliInfo(Form("Finish building stack propagation layers. nClusters %d.", nClStack));
+ if(nClStack < kFindable) continue;
+ ntracks += Clusters2TracksStack(&stackLayer[0], &esdTrackList);
+ }
+ //AliInfo(Form("Found %d tracks in SM", ntracks));
+
+ for(int itrack=0; itrack<ntracks; itrack++)
+ esd->AddTrack((AliESDtrack*)esdTrackList[itrack]);
+
+ return ntracks;
+}
+
+//____________________________________________________________________
+Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
+ , TClonesArray *esdTrackList)
+{
+ //
+ // Make tracks in one TRD stack.
+ //
+ // Parameters :
+ // layer : Array of stack propagation layers containing clusters
+ // esdTrackList : Array of ESD tracks found by the stand alone tracker.
+ // On exit the tracks found in this stack are appended.
+ //
+ // Output :
+ // Number of tracks found in this stack.
+ //
+ // Detailed description
+ //
+ // 1. Find the 3 most useful seeding chambers. See BuildSeedingConfigs() for details.
+ // 2. Steer AliTRDtrackerV1::MakeSeeds() for 3 seeding layer configurations.
+ // See AliTRDtrackerV1::MakeSeeds() for more details.
+ // 3. Arrange track candidates in decreasing order of their quality
+ // 4. Classify tracks in 5 categories according to:
+ // a) number of layers crossed
+ // b) track quality
+ // 5. Sign clusters by tracks in decreasing order of track quality
+ // 6. Build AliTRDtrack out of seeding tracklets
+ // 7. Cook MC label
+ // 8. Build ESD track and register it to the output list
+ //
+
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ Int_t nTimeBins = cal->GetNumberOfTimeBins();
+ AliTRDseedV1 sseed[kMaxTracksStack*6]; // to be initialized
+ Int_t pars[3]; // MakeSeeds parameters
+
+ //Double_t alpha = AliTRDgeometry::GetAlpha();
+ //Double_t shift = .5 * alpha;
+ Int_t configs[kNConfigs];
+
+ // Build initial seeding configurations
+ Double_t quality = BuildSeedingConfigs(layer, configs);
+#ifdef DEBUG
+ if(AliTRDReconstructor::StreamLevel() > 1)
+ AliInfo(Form("Plane config %d %d %d Quality %f"
+ , configs[0], configs[1], configs[2], quality));
+#endif
+
+ // Initialize contors
+ Int_t ntracks, // number of TRD track candidates
+ ntracks1, // number of registered TRD tracks/iter
+ ntracks2 = 0; // number of all registered TRD tracks in stack
+ fSieveSeeding = 0;
+ do{
+ // Loop over seeding configurations
+ ntracks = 0; ntracks1 = 0;
+ for (Int_t iconf = 0; iconf<3; iconf++) {
+ pars[0] = configs[iconf];
+ pars[1] = layer->GetStackNr();
+ pars[2] = ntracks;
+ ntracks = MakeSeeds(layer, &sseed[6*ntracks], pars);
+ if(ntracks == kMaxTracksStack) break;
+ }
+#ifdef DEBUG
+ if(AliTRDReconstructor::StreamLevel() > 1)
+ AliInfo(Form("Candidate TRD tracks %d in stack %d.", ntracks, pars[1]));
+#endif
+ if(!ntracks) break;
+
+ // Sort the seeds according to their quality
+ Int_t sort[kMaxTracksStack];
+ TMath::Sort(ntracks, fTrackQuality, sort, kTRUE);
+
+ // Initialize number of tracks so far and logic switches
+ Int_t ntracks0 = esdTrackList->GetEntriesFast();
+ Bool_t signedTrack[kMaxTracksStack];
+ Bool_t fakeTrack[kMaxTracksStack];
+ for (Int_t i=0; i<ntracks; i++){
+ signedTrack[i] = kFALSE;
+ fakeTrack[i] = kFALSE;
+ }
+ //AliInfo("Selecting track candidates ...");
+
+ // Sieve clusters in decreasing order of track quality
+ Double_t trackParams[7];
+// AliTRDseedV1 *lseed = 0x0;
+ Int_t jSieve = 0, candidates;
+ do{
+ //AliInfo(Form("\t\tITER = %i ", jSieve));
+
+ // Check track candidates
+ candidates = 0;
+ for (Int_t itrack = 0; itrack < ntracks; itrack++) {
+ Int_t trackIndex = sort[itrack];
+ if (signedTrack[trackIndex] || fakeTrack[trackIndex]) continue;
+
+
+ // Calculate track parameters from tracklets seeds
+ Int_t labelsall[1000];
+ Int_t nlabelsall = 0;
+ Int_t naccepted = 0;
+ Int_t ncl = 0;
+ Int_t nused = 0;
+ Int_t nlayers = 0;
+ Int_t findable = 0;
+ for (Int_t jLayer = 0; jLayer < kNPlanes; jLayer++) {
+ Int_t jseed = kNPlanes*trackIndex+jLayer;
+ if (TMath::Abs(sseed[jseed].GetYref(0) / sseed[jseed].GetX0()) < 0.15)
+ findable++;
+
+ if(!sseed[jseed].IsOK()) continue;
+ sseed[jseed].UpdateUsed();
+ ncl += sseed[jseed].GetN2();
+ nused += sseed[jseed].GetNUsed();
+ nlayers++;
+
+ // Cooking label
+ for (Int_t itime = 0; itime < nTimeBins; itime++) {
+ if(!sseed[jseed].IsUsable(itime)) continue;
+ naccepted++;
+ Int_t tindex = 0, ilab = 0;
+ while(ilab<3 && (tindex = sseed[jseed].GetClusters(itime)->GetLabel(ilab)) >= 0){
+ labelsall[nlabelsall++] = tindex;
+ ilab++;
+ }
+ }
+ }
+ // Filter duplicated tracks
+ if (nused > 30){
+ //printf("Skip nused %d\n", nused);
+ fakeTrack[trackIndex] = kTRUE;
+ continue;
+ }
+ if (Float_t(nused)/ncl >= .25){
+ //printf("Skip nused/ncl >= .25\n");
+ fakeTrack[trackIndex] = kTRUE;
+ continue;
+ }
+
+ // Classify tracks
+ Bool_t skip = kFALSE;
+ switch(jSieve){
+ case 0:
+ if(nlayers < 6) {skip = kTRUE; break;}
+ if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;}
+ break;
+
+ case 1:
+ if(nlayers < findable){skip = kTRUE; break;}
+ if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -4.){skip = kTRUE; break;}
+ break;
+
+ case 2:
+ if ((nlayers == findable) || (nlayers == 6)) { skip = kTRUE; break;}
+ if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -6.0){skip = kTRUE; break;}
+ break;
+
+ case 3:
+ if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;}
+ break;
+
+ case 4:
+ if (nlayers == 3){skip = kTRUE; break;}
+ if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) - nused/(nlayers-3.0) < -15.0){skip = kTRUE; break;}
+ break;
+ }
+ if(skip){
+ candidates++;
+ //printf("REJECTED : %d [%d] nlayers %d trackQuality = %e nused %d\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused);
+ continue;
+ }
+ signedTrack[trackIndex] = kTRUE;
+
+
+ // Build track label - what happens if measured data ???
+ Int_t labels[1000];
+ Int_t outlab[1000];
+ Int_t nlab = 0;
+ for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+ Int_t jseed = kNPlanes*trackIndex+iLayer;
+ if(!sseed[jseed].IsOK()) continue;
+ for(int ilab=0; ilab<2; ilab++){
+ if(sseed[jseed].GetLabels(ilab) < 0) continue;
+ labels[nlab] = sseed[jseed].GetLabels(ilab);
+ nlab++;
+ }
+ }
+ Freq(nlab,labels,outlab,kFALSE);
+ Int_t label = outlab[0];
+ Int_t frequency = outlab[1];
+ Freq(nlabelsall,labelsall,outlab,kFALSE);
+ Int_t label1 = outlab[0];
+ Int_t label2 = outlab[2];
+ Float_t fakeratio = (naccepted - outlab[1]) / Float_t(naccepted);
+
+
+ // Sign clusters
+ AliTRDcluster *cl = 0x0; Int_t clusterIndex = -1;
+ for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
+ Int_t jseed = kNPlanes*trackIndex+jLayer;
+ if(!sseed[jseed].IsOK()) continue;
+ if(TMath::Abs(sseed[jseed].GetYfit(1) - sseed[jseed].GetYfit(1)) >= .2) continue; // check this condition with Marian
+ sseed[jseed].UseClusters();
+ if(!cl){
+ Int_t ic = 0;
+ while(!(cl = sseed[jseed].GetClusters(ic))) ic++;
+ clusterIndex = sseed[jseed].GetIndexes(ic);
+ }
+ }
+ if(!cl) continue;
+
+
+ // Build track parameters
+ AliTRDseedV1 *lseed =&sseed[trackIndex*6];
+ Int_t idx = 0;
+ while(idx<3 && !lseed->IsOK()) {
+ idx++;
+ lseed++;
+ }
+ Double_t cR = lseed->GetC();
+ trackParams[1] = lseed->GetYref(0);
+ trackParams[2] = lseed->GetZref(0);
+ trackParams[3] = lseed->GetX0() * cR - TMath::Sin(TMath::ATan(lseed->GetYref(1)));
+ trackParams[4] = lseed->GetZref(1) / TMath::Sqrt(1. + lseed->GetYref(1) * lseed->GetYref(1));
+ trackParams[5] = cR;
+ trackParams[0] = lseed->GetX0();
+ trackParams[6] = layer[0].GetSector();/* *alpha+shift; // Supermodule*/
+
+#ifdef DEBUG
+ if(AliTRDReconstructor::StreamLevel() > 1) printf("Track %d [%d] nlayers %d trackQuality = %e nused %d, yref = %3.3f\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused, trackParams[1]);
+
+ if(AliTRDReconstructor::StreamLevel() >= 1){
+ Int_t sector = layer[0].GetSector();
+ Int_t nclusters = 0;
+ AliTRDseedV1 *dseed[6];
+ for(int is=0; is<6; is++){
+ dseed[is] = new AliTRDseedV1(sseed[trackIndex*6+is], kTRUE);
+ nclusters += sseed[is].GetN2();
+ //for(int ic=0; ic<30; ic++) if(sseed[trackIndex*6+is].GetClusters(ic)) printf("l[%d] tb[%d] cptr[%p]\n", is, ic, sseed[trackIndex*6+is].GetClusters(ic));
+ }
+ //Int_t eventNrInFile = esd->GetEventNumberInFile();
+ //AliInfo(Form("Number of clusters %d.", nclusters));
+ TTreeSRedirector &cstreamer = *fDebugStreamer;
+ cstreamer << "Clusters2TracksStack"
+ << "Iter=" << fSieveSeeding
+ << "Like=" << fTrackQuality[trackIndex]
+ << "S0.=" << dseed[0]
+ << "S1.=" << dseed[1]
+ << "S2.=" << dseed[2]
+ << "S3.=" << dseed[3]
+ << "S4.=" << dseed[4]
+ << "S5.=" << dseed[5]
+ << "p0=" << trackParams[0]
+ << "p1=" << trackParams[1]
+ << "p2=" << trackParams[2]
+ << "p3=" << trackParams[3]
+ << "p4=" << trackParams[4]
+ << "p5=" << trackParams[5]
+ << "p6=" << trackParams[6]
+ << "Sector=" << sector
+ << "Stack=" << pars[1]
+ << "Label=" << label
+ << "Label1=" << label1
+ << "Label2=" << label2
+ << "FakeRatio=" << fakeratio
+ << "Freq=" << frequency
+ << "Ncl=" << ncl
+ << "NLayers=" << nlayers
+ << "Findable=" << findable
+ << "NUsed=" << nused
+ << "\n";
+ //???for(int is=0; is<6; is++) delete dseed[is];
+ }
+#endif
+
+ AliTRDtrack *track = AliTRDtrackerV1::RegisterSeed(&sseed[trackIndex*kNPlanes], trackParams);
+ if(!track){
+ AliWarning("Fail to build a TRD Track.");
+ continue;
+ }
+ AliESDtrack esdTrack;
+ esdTrack.UpdateTrackParams(track, AliESDtrack::kTRDout);
+ esdTrack.SetLabel(track->GetLabel());
+ new ((*esdTrackList)[ntracks0++]) AliESDtrack(esdTrack);
+ ntracks1++;
+ }
+
+ jSieve++;
+ } while(jSieve<5 && candidates); // end track candidates sieve
+ if(!ntracks1) break;
+
+ // increment counters
+ ntracks2 += ntracks1;
+ fSieveSeeding++;
+
+ // Rebuild plane configurations and indices taking only unused clusters into account
+ quality = BuildSeedingConfigs(layer, configs);
+ //if(quality < fRecoParam->GetPlaneQualityThreshold()) break;
+
+ for(Int_t il = 0; il < kNPlanes * nTimeBins; il++) layer[il].BuildIndices(fSieveSeeding);
+
+#ifdef DEBUG
+ if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Sieve level %d Plane config %d %d %d Quality %f", fSieveSeeding, configs[0], configs[1], configs[2], quality));
+#endif
+ } while(fSieveSeeding<10); // end stack clusters sieve
+
+
+
+ //AliInfo(Form("Registered TRD tracks %d in stack %d.", ntracks2, pars[1]));
+
+ return ntracks2;
+}
+
+//___________________________________________________________________
+Double_t AliTRDtrackerV1::BuildSeedingConfigs(AliTRDstackLayer *layers
+ , Int_t *configs)
+{
+ //
+ // Assign probabilities to chambers according to their
+ // capability of producing seeds.
+ //
+ // Parameters :
+ //
+ // layers : Array of stack propagation layers for all 6 chambers in one stack
+ // configs : On exit array of configuration indexes (see GetSeedingConfig()
+ // for details) in the decreasing order of their seeding probabilities.
+ //
+ // Output :
+ //
+ // Return top configuration quality
+ //
+ // Detailed description:
+ //
+ // To each chamber seeding configuration (see GetSeedingConfig() for
+ // the list of all configurations) one defines 2 quality factors:
+ // - an apriori topological quality (see GetSeedingConfig() for details) and
+ // - a data quality based on the uniformity of the distribution of
+ // clusters over the x range (time bins population). See CookChamberQA() for details.
+ // The overall chamber quality is given by the product of this 2 contributions.
+ //
+
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ Int_t nTimeBins = cal->GetNumberOfTimeBins();
+
+ Double_t chamberQA[kNPlanes];
+ for(int iplane=0; iplane<kNPlanes; iplane++){
+ chamberQA[iplane] = CookPlaneQA(&layers[iplane*nTimeBins]);
+ //printf("chamberQA[%d] = %f\n", iplane, chamberQA[iplane]);
+ }
+
+ Double_t tconfig[kNConfigs];
+ Int_t planes[4];
+ for(int iconf=0; iconf<kNConfigs; iconf++){
+ GetSeedingConfig(iconf, planes);
+ tconfig[iconf] = fTopologicQA[iconf];
+ for(int iplane=0; iplane<4; iplane++) tconfig[iconf] *= chamberQA[planes[iplane]];
+ }
+
+ TMath::Sort(kNConfigs, tconfig, configs, kTRUE);
+ return tconfig[configs[0]];
+}
+
+//____________________________________________________________________
+Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
+ , AliTRDseedV1 *sseed
+ , Int_t *ipar)
+{
+ //
+ // Make tracklet seeds in the TRD stack.
+ //
+ // Parameters :
+ // layers : Array of stack propagation layers containing clusters
+ // sseed : Array of empty tracklet seeds. On exit they are filled.
+ // ipar : Control parameters:
+ // ipar[0] -> seeding chambers configuration
+ // ipar[1] -> stack index
+ // ipar[2] -> number of track candidates found so far
+ //
+ // Output :
+ // Number of tracks candidates found.
+ //
+ // Detailed description
+ //
+ // The following steps are performed:
+ // 1. Select seeding layers from seeding chambers
+ // 2. Select seeding clusters from the seeding AliTRDpropagationLayerStack.
+ // The clusters are taken from layer 3, layer 0, layer 1 and layer 2, in
+ // this order. The parameters controling the range of accepted clusters in
+ // layer 0, 1, and 2 are defined in AliTRDstackLayer::BuildCond().
+ // 3. Helix fit of the cluster set. (see AliTRDtrackerFitter::FitRieman(AliTRDcluster**))
+ // 4. Initialize seeding tracklets in the seeding chambers.
+ // 5. Filter 0.
+ // Chi2 in the Y direction less than threshold ... (1./(3. - sLayer))
+ // Chi2 in the Z direction less than threshold ... (1./(3. - sLayer))
+ // 6. Attach clusters to seeding tracklets and find linear approximation of
+ // the tracklet (see AliTRDseedV1::AttachClustersIter()). The number of used
+ // clusters used by current seeds should not exceed ... (25).
+ // 7. Filter 1.
+ // All 4 seeding tracklets should be correctly constructed (see
+ // AliTRDseedV1::AttachClustersIter())
+ // 8. Helix fit of the seeding tracklets
+ // 9. Filter 2.
+ // Likelihood calculation of the fit. (See AliTRDtrackerV1::CookLikelihood() for details)
+ // 10. Extrapolation of the helix fit to the other 2 chambers:
+ // a) Initialization of extrapolation tracklet with fit parameters
+ // b) Helix fit of tracklets
+ // c) Attach clusters and linear interpolation to extrapolated tracklets
+ // d) Helix fit of tracklets
+ // 11. Improve seeding tracklets quality by reassigning clusters.
+ // See AliTRDtrackerV1::ImproveSeedQuality() for details.
+ // 12. Helix fit of all 6 seeding tracklets and chi2 calculation
+ // 13. Hyperplane fit and track quality calculation. See AliTRDtrackerFitter::FitHyperplane() for details.
+ // 14. Cooking labels for tracklets. Should be done only for MC
+ // 15. Register seeds.
+ //
+
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ Int_t nTimeBins = cal->GetNumberOfTimeBins();
+ AliTRDcluster *c[4] = {0x0, 0x0, 0x0, 0x0}; // initilize seeding clusters
+ AliTRDseedV1 *cseed = &sseed[0]; // initialize tracklets for first track
+ Int_t ncl, mcl; // working variable for looping over clusters
+ Int_t index[AliTRDstackLayer::kMaxClustersLayer], jndex[AliTRDstackLayer::kMaxClustersLayer];
+ // chi2 storage
+ // chi2[0] = tracklet chi2 on the Z direction
+ // chi2[1] = tracklet chi2 on the R direction
+ Double_t chi2[4];
+
+
+ // this should be data member of AliTRDtrack
+ Double_t seedQuality[kMaxTracksStack];
+
+ // unpack control parameters
+ Int_t config = ipar[0];
+ Int_t istack = ipar[1];
+ Int_t ntracks = ipar[2];
+ Int_t planes[kNSeedPlanes]; GetSeedingConfig(config, planes);
+#ifdef DEBUG
+ if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Making seeds Stack[%d] Config[%d] Tracks[%d]...", istack, config, ntracks));
+#endif
+
+ // Init chambers geometry
+ Double_t hL[kNPlanes]; // Tilting angle
+ Float_t padlength[kNPlanes]; // pad lenghts
+ AliTRDpadPlane *pp;
+ for(int il=0; il<kNPlanes; il++){
+ pp = fGeom->GetPadPlane(il, istack); // istack has to be imported
+ hL[il] = TMath::Tan(-TMath::DegToRad()*pp->GetTiltingAngle());
+ padlength[il] = 10.; //pp->GetLengthIPad();
+ }
+
+ Double_t cond0[4], cond1[4], cond2[4];
+ // make seeding layers (to be moved in Clusters2TracksStack)
+ AliTRDstackLayer *layer[] = {0x0, 0x0, 0x0, 0x0};
+ for(int isl=0; isl<kNSeedPlanes; isl++) layer[isl] = MakeSeedingLayer(&layers[planes[isl] * nTimeBins], planes[isl]);
+
+
+ // Start finding seeds
+ Int_t icl = 0;
+ while((c[3] = (*layer[3])[icl++])){
+ if(!c[3]) continue;
+ layer[0]->BuildCond(c[3], cond0, 0);
+ layer[0]->GetClusters(cond0, index, ncl);
+ Int_t jcl = 0;
+ while(jcl<ncl) {
+ c[0] = (*layer[0])[index[jcl++]];
+ if(!c[0]) continue;
+ Double_t dx = c[3]->GetX() - c[0]->GetX();
+ Double_t theta = (c[3]->GetZ() - c[0]->GetZ())/dx;
+ Double_t phi = (c[3]->GetY() - c[0]->GetY())/dx;
+ layer[1]->BuildCond(c[0], cond1, 1, theta, phi);
+ layer[1]->GetClusters(cond1, jndex, mcl);
+
+ Int_t kcl = 0;
+ while(kcl<mcl) {
+ c[1] = (*layer[1])[jndex[kcl++]];
+ if(!c[1]) continue;
+ layer[2]->BuildCond(c[1], cond2, 2, theta, phi);
+ c[2] = layer[2]->GetNearestCluster(cond2);
+ if(!c[2]) continue;
+
+ //AliInfo("Seeding clusters found. Building seeds ...");
+ //for(Int_t i = 0; i < kNSeedPlanes; i++) printf("%i. coordinates: x = %3.3f, y = %3.3f, z = %3.3f\n", i, c[i]->GetX(), c[i]->GetY(), c[i]->GetZ());
+ for (Int_t il = 0; il < 6; il++) cseed[il].Reset();
+
+ fFitter->Reset();
+
+ fFitter->FitRieman(c, kNSeedPlanes);
+
+ chi2[0] = 0.; chi2[1] = 0.;
+ AliTRDseedV1 *tseed = 0x0;
+ for(int iLayer=0; iLayer<kNSeedPlanes; iLayer++){
+ tseed = &cseed[planes[iLayer]];
+ tseed->SetRecoParam(fRecoParam);
+ tseed->SetLayer(planes[iLayer]);
+ tseed->SetTilt(hL[planes[iLayer]]);
+ tseed->SetPadLength(TMath::Sqrt(c[iLayer]->GetSigmaZ2()*12));
+ tseed->SetX0(layer[iLayer]->GetX());
+
+ tseed->Update(fFitter->GetRiemanFitter());
+ chi2[0] += tseed->GetChi2Z(c[iLayer]->GetZ());
+ chi2[1] += tseed->GetChi2Y(c[iLayer]->GetY());
+ }
+
+ Bool_t isFake = kFALSE;
+ 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;
+#ifdef DEBUG
+ if(AliTRDReconstructor::StreamLevel() >= 2){
+ Float_t yref[4], ycluster[4];
+ for(int il=0; il<4; il++){
+ tseed = &cseed[planes[il]];
+ yref[il] = tseed->GetYref(0);
+ ycluster[il] = c[il]->GetY();
+ }
+ Float_t threshold = .5;//1./(3. - sLayer);
+ Int_t ll = c[3]->GetLabel(0);
+ TTreeSRedirector &cs0 = *fDebugStreamer;
+ cs0 << "MakeSeeds0"
+ <<"isFake=" << isFake
+ <<"label=" << ll
+ <<"threshold=" << threshold
+ <<"chi2=" << chi2[1]
+ <<"yref0="<<yref[0]
+ <<"yref1="<<yref[1]
+ <<"yref2="<<yref[2]
+ <<"yref3="<<yref[3]
+ <<"ycluster0="<<ycluster[0]
+ <<"ycluster1="<<ycluster[1]
+ <<"ycluster2="<<ycluster[2]
+ <<"ycluster3="<<ycluster[3]
+ <<"\n";
+ }
+#endif
+
+ if(chi2[0] > fRecoParam->GetChi2Z()/*7./(3. - sLayer)*//*iter*/){
+ //AliInfo(Form("Failed chi2 filter on chi2Z [%f].", chi2[0]));
+ continue;
+ }
+ if(chi2[1] > fRecoParam->GetChi2Y()/*1./(3. - sLayer)*//*iter*/){
+ //AliInfo(Form("Failed chi2 filter on chi2Y [%f].", chi2[1]));
+ continue;
+ }
+ //AliInfo("Passed chi2 filter.");
+
+#ifdef DEBUG
+ if(AliTRDReconstructor::StreamLevel() >= 2){
+ Float_t minmax[2] = { -100.0, 100.0 };
+ for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
+ Float_t max = c[iLayer]->GetZ() + cseed[planes[iLayer]].GetPadLength() * 0.5 + 1.0 - cseed[planes[iLayer]].GetZref(0);
+ if (max < minmax[1]) minmax[1] = max;
+ Float_t min = c[iLayer]->GetZ()-cseed[planes[iLayer]].GetPadLength() * 0.5 - 1.0 - cseed[planes[iLayer]].GetZref(0);
+ if (min > minmax[0]) minmax[0] = min;
+ }
+ Double_t xpos[4];
+ for(Int_t l = 0; l < kNSeedPlanes; l++) xpos[l] = layer[l]->GetX();
+ TTreeSRedirector &cstreamer = *fDebugStreamer;
+ cstreamer << "MakeSeeds1"
+ << "isFake=" << isFake
+ << "config=" << config
+ << "Cl0.=" << c[0]
+ << "Cl1.=" << c[1]
+ << "Cl2.=" << c[2]
+ << "Cl3.=" << c[3]
+ << "X0=" << xpos[0] //layer[sLayer]->GetX()
+ << "X1=" << xpos[1] //layer[sLayer + 1]->GetX()
+ << "X2=" << xpos[2] //layer[sLayer + 2]->GetX()
+ << "X3=" << xpos[3] //layer[sLayer + 3]->GetX()
+ << "Y2exp=" << cond2[0]
+ << "Z2exp=" << cond2[1]
+ << "Chi2R=" << chi2[0]
+ << "Chi2Z=" << chi2[1]
+ << "Seed0.=" << &cseed[planes[0]]
+ << "Seed1.=" << &cseed[planes[1]]
+ << "Seed2.=" << &cseed[planes[2]]
+ << "Seed3.=" << &cseed[planes[3]]
+ << "Zmin=" << minmax[0]
+ << "Zmax=" << minmax[1]
+ << "\n" ;
+ }
+#endif
+ // try attaching clusters to tracklets
+ Int_t nUsedCl = 0;
+ Int_t nlayers = 0;
+ for(int iLayer=0; iLayer<kNSeedPlanes; iLayer++){
+ AliTRDseedV1 tseed = cseed[planes[iLayer]];
+ if(!tseed.AttachClustersIter(&layers[planes[iLayer]*nTimeBins], 5., kFALSE, c[iLayer])) continue;
+ cseed[planes[iLayer]] = tseed;
+ nUsedCl += cseed[planes[iLayer]].GetNUsed();
+ if(nUsedCl > 25) break;
+ nlayers++;
+ }
+ if(nlayers < kNSeedPlanes){
+ //AliInfo("Failed updating all seeds.");
+ continue;
+ }
+ // fit tracklets and cook likelihood
+ chi2[0] = 0.; chi2[1] = 0.;
+ fFitter->FitRieman(&cseed[0], &planes[0]);
+ AliRieman *rim = fFitter->GetRiemanFitter();
+ for(int iLayer=0; iLayer<4; iLayer++){
+ cseed[planes[iLayer]].Update(rim);
+ chi2[0] += (Float_t)cseed[planes[iLayer]].GetChi2Z();
+ chi2[1] += cseed[planes[iLayer]].GetChi2Y();
+ }
+ Double_t chi2r = chi2[1], chi2z = chi2[0];
+ Double_t like = CookLikelihood(&cseed[0], planes, chi2); // to be checked
+ if (TMath::Log(1.E-9 + like) < fRecoParam->GetTrackLikelihood()){
+ //AliInfo(Form("Failed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
+ continue;
+ }
+ //AliInfo(Form("Passed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
+
+
+ // book preliminary results
+ seedQuality[ntracks] = like;
+ fSeedLayer[ntracks] = config;/*sLayer;*/
+
+ // attach clusters to the extrapolation seeds
+ Int_t lextrap[2];
+ GetExtrapolationConfig(config, lextrap);
+ Int_t nusedf = 0; // debug value
+ for(int iLayer=0; iLayer<2; iLayer++){
+ Int_t jLayer = lextrap[iLayer];
+
+ // prepare extrapolated seed
+ cseed[jLayer].Reset();
+ cseed[jLayer].SetRecoParam(fRecoParam);
+ cseed[jLayer].SetLayer(jLayer);
+ cseed[jLayer].SetTilt(hL[jLayer]);
+ cseed[jLayer].SetX0(layers[jLayer * nTimeBins + (nTimeBins/2)].GetX()); // ????????
+ //cseed[jLayer].SetPadLength(??????????);
+ cseed[jLayer].Update(rim);
+
+ AliTRDcluster *cd = FindSeedingCluster(&layers[jLayer*nTimeBins], &cseed[jLayer]);
+ if(cd == 0x0) continue;
+// if(!cd) continue;
+ cseed[jLayer].SetPadLength(TMath::Sqrt(cd->GetSigmaZ2() * 12.));
+ cseed[jLayer].SetX0(cd->GetX()); // reference defined by a seedingLayer which is defined by the x-coordinate of the layers inside
+
+ // fit extrapolated seed
+ AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
+ if ((jLayer == 0) && !(cseed[1].IsOK())) continue;
+ if ((jLayer == 5) && !(cseed[4].IsOK())) continue;
+ AliTRDseedV1 tseed = cseed[jLayer];
+ if(!tseed.AttachClustersIter(&layers[jLayer*nTimeBins], 1000.)) continue;
+ cseed[jLayer] = tseed;
+ nusedf += cseed[jLayer].GetNUsed(); // debug value
+ AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
+ }
+ //AliInfo("Extrapolation done.");
+
+ ImproveSeedQuality(layers, cseed);
+ //AliInfo("Improve seed quality done.");
+
+ nlayers = 0;
+ Int_t nclusters = 0;
+ Int_t findable = 0;
+ for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+ if (TMath::Abs(cseed[iLayer].GetYref(0) / cseed[iLayer].GetX0()) < 0.15) findable++;
+ if (!cseed[iLayer].IsOK()) continue;
+ nclusters += cseed[iLayer].GetN2();
+ nlayers++;
+ }
+ if (nlayers < 3){
+ //AliInfo("Failed quality check on seeds.");
+ continue;
+ }
+
+ // fit full track and cook likelihoods
+ fFitter->FitRieman(&cseed[0]);
+ Double_t chi2ZF = 0., chi2RF = 0.;
+ for(int ilayer=0; ilayer<6; ilayer++){
+ cseed[ilayer].Update(fFitter->GetRiemanFitter());
+ if (!cseed[ilayer].IsOK()) continue;
+ //tchi2 = cseed[ilayer].GetChi2Z();
+ //printf("layer %d chi2 %e\n", ilayer, tchi2);
+ chi2ZF += cseed[ilayer].GetChi2Z();
+ chi2RF += cseed[ilayer].GetChi2Y();
+ }
+ chi2ZF /= TMath::Max((nlayers - 3.), 1.);
+ chi2RF /= TMath::Max((nlayers - 3.), 1.);
+
+ // do the final track fitting
+ fFitter->SetLayers(nlayers);
+#ifdef DEBUG
+ fFitter->SetDebugStream(fDebugStreamer);
+#endif
+ fTrackQuality[ntracks] = fFitter->FitHyperplane(&cseed[0], chi2ZF, GetZ());
+ Double_t param[3];
+ Double_t chi2[2];
+ fFitter->GetHyperplaneFitResults(param);
+ fFitter->GetHyperplaneFitChi2(chi2);
+ //AliInfo("Hyperplane fit done\n");
+
+ // finalize tracklets
+ Int_t labels[12];
+ Int_t outlab[24];
+ Int_t nlab = 0;
+ for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+ if (!cseed[iLayer].IsOK()) continue;
+
+ if (cseed[iLayer].GetLabels(0) >= 0) {
+ labels[nlab] = cseed[iLayer].GetLabels(0);
+ nlab++;
+ }
+
+ if (cseed[iLayer].GetLabels(1) >= 0) {
+ labels[nlab] = cseed[iLayer].GetLabels(1);
+ nlab++;
+ }
+ }
+ Freq(nlab,labels,outlab,kFALSE);
+ Int_t label = outlab[0];
+ Int_t frequency = outlab[1];
+ for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+ cseed[iLayer].SetFreq(frequency);
+ cseed[iLayer].SetC(param[1]/*cR*/);
+ cseed[iLayer].SetCC(param[0]/*cC*/);
+ cseed[iLayer].SetChi2(chi2[0]);
+ cseed[iLayer].SetChi2Z(chi2ZF);
+ }
+
+#ifdef DEBUG
+ if(AliTRDReconstructor::StreamLevel() >= 2){
+ Double_t curv = (fFitter->GetRiemanFitter())->GetC();
+ TTreeSRedirector &cstreamer = *fDebugStreamer;
+ cstreamer << "MakeSeeds2"
+ << "C=" << curv
+ << "Chi2R=" << chi2r
+ << "Chi2Z=" << chi2z
+ << "Chi2TR=" << chi2[0]
+ << "Chi2TC=" << chi2[1]
+ << "Chi2RF=" << chi2RF
+ << "Chi2ZF=" << chi2ZF
+ << "Ncl=" << nclusters
+ << "Nlayers=" << nlayers
+ << "NUsedS=" << nUsedCl
+ << "NUsed=" << nusedf
+ << "Findable" << findable
+ << "Like=" << like
+ << "S0.=" << &cseed[0]
+ << "S1.=" << &cseed[1]
+ << "S2.=" << &cseed[2]
+ << "S3.=" << &cseed[3]
+ << "S4.=" << &cseed[4]
+ << "S5.=" << &cseed[5]
+ << "Label=" << label
+ << "Freq=" << frequency
+ << "\n";
+ }
+#endif
+
+ ntracks++;
+ if(ntracks == kMaxTracksStack){
+ AliWarning(Form("Number of seeds reached maximum allowed (%d) in stack.", kMaxTracksStack));
+ return ntracks;
+ }
+ cseed += 6;
+ }
+ }
+ }
+ for(int isl=0; isl<4; isl++) delete layer[isl];
+
+ return ntracks;
+}
+
+//_____________________________________________________________________________
+AliTRDtrack *AliTRDtrackerV1::RegisterSeed(AliTRDseedV1 *seeds, Double_t *params)
+{
+ //
+ // Build a TRD track out of tracklet candidates
+ //
+ // Parameters :
+ // seeds : array of tracklets
+ // params : track parameters (see MakeSeeds() function body for a detailed description)
+ //
+ // Output :
+ // The TRD track.
+ //
+ // Detailed description
+ //
+ // To be discussed with Marian !!
+ //
+
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ Int_t nTimeBins = cal->GetNumberOfTimeBins();
+
+ Double_t alpha = AliTRDgeometry::GetAlpha();
+ Double_t shift = AliTRDgeometry::GetAlpha()/2.0;
+ Double_t c[15];
+
+ c[ 0] = 0.2;
+ c[ 1] = 0.0; c[ 2] = 2.0;
+ c[ 3] = 0.0; c[ 4] = 0.0; c[ 5] = 0.02;
+ 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;
+
+ Int_t index = 0;
+ AliTRDcluster *cl = 0;
+
+ for (Int_t ilayer = 0; ilayer < 6; ilayer++) {
+ if (seeds[ilayer].IsOK()) {
+ for (Int_t itime = nTimeBins - 1; itime > 0; itime--) {
+ if (seeds[ilayer].GetIndexes(itime) > 0) {
+ index = seeds[ilayer].GetIndexes(itime);
+ cl = seeds[ilayer].GetClusters(itime);
+ break;
+ }
+ }
+ }
+ if (index > 0) {
+ break;
+ }
+ }
+ if (cl == 0) return 0;
+ AliTRDtrack *track = new AliTRDtrack(cl
+ ,index
+ ,¶ms[1]
+ ,c
+ ,params[0]
+ ,params[6]*alpha+shift);
+ // SetCluster(cl, 0); // A. Bercuci
+ track->PropagateTo(params[0]-5.0);
+ track->ResetCovariance(1);
+ Int_t rc = FollowBackProlongation(*track);
+ if (rc < 30) {
+ delete track;
+ track = 0;
+ }
+ else {
+ track->CookdEdx();
+ track->CookdEdxTimBin(-1);
+ CookLabel(track,0.9);
+ }
+
+ return track;
+}
+
+//____________________________________________________________________
+void AliTRDtrackerV1::ImproveSeedQuality(AliTRDstackLayer *layers
+ , AliTRDseedV1 *cseed)
+{
+ //
+ // Sort tracklets according to "quality" and try to "improve" the first 4 worst
+ //
+ // Parameters :
+ // layers : Array of propagation layers for a stack/supermodule
+ // cseed : Array of 6 seeding tracklets which has to be improved
+ //
+ // Output :
+ // cssed : Improved seeds
+ //
+ // Detailed description
+ //
+ // Iterative procedure in which new clusters are searched for each
+ // tracklet seed such that the seed quality (see AliTRDseed::GetQuality())
+ // can be maximized. If some optimization is found the old seeds are replaced.
+ //
+
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ Int_t nTimeBins = cal->GetNumberOfTimeBins();
+
+ // make a local working copy
+ AliTRDseedV1 bseed[6];
+ for (Int_t jLayer = 0; jLayer < 6; jLayer++) bseed[jLayer] = cseed[jLayer];
+
+
+ Float_t lastquality = 10000.0;
+ Float_t lastchi2 = 10000.0;
+ Float_t chi2 = 1000.0;
+
+ for (Int_t iter = 0; iter < 4; iter++) {
+ Float_t sumquality = 0.0;
+ Float_t squality[6];
+ Int_t sortindexes[6];
+
+ for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
+ squality[jLayer] = bseed[jLayer].IsOK() ? bseed[jLayer].GetQuality(kTRUE) : -1.;
+ sumquality +=squality[jLayer];
+ }
+ if ((sumquality >= lastquality) || (chi2 > lastchi2)) break;
+
+
+ lastquality = sumquality;
+ lastchi2 = chi2;
+ if (iter > 0) for (Int_t jLayer = 0; jLayer < 6; jLayer++) cseed[jLayer] = bseed[jLayer];
+
+
+ TMath::Sort(6, squality, sortindexes, kFALSE);
+ for (Int_t jLayer = 5; jLayer > 1; jLayer--) {
+ Int_t bLayer = sortindexes[jLayer];
+ bseed[bLayer].AttachClustersIter(&layers[bLayer*nTimeBins], squality[bLayer], kTRUE);
+ }
+
+ chi2 = AliTRDseedV1::FitRiemanTilt(bseed,kTRUE);
+ } // Loop: iter
+}
+
+//____________________________________________________________________
+Double_t AliTRDtrackerV1::CookPlaneQA(AliTRDstackLayer *layers)
+{
+ //
+ // Calculate plane 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
+ //
+
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ Int_t nTimeBins = cal->GetNumberOfTimeBins();
+
+// Double_t x;
+// TLinearFitter fitter(1, "pol1");
+// fitter.ClearPoints();
+ Int_t ncl = 0;
+ Int_t nused = 0;
+ Int_t nClLayer;
+ for(int itb=0; itb<nTimeBins; itb++){
+ //x = layer[itb].GetX();
+ //printf("x[%d] = %f nCls %d\n", itb, x, layer[itb].GetNClusters());
+ //if(!layer[itb].GetNClusters()) continue;
+ //fitter.AddPoint(&x, layer[itb].GetNClusters(), 1.);
+ nClLayer = layers[itb].GetNClusters();
+ ncl += nClLayer;
+ for(Int_t incl = 0; incl < nClLayer; incl++)
+ if((layers[itb].GetCluster(incl))->IsUsed()) nused++;
+ }
+
+ // calculate the deviation of the mean number of clusters from the
+ // closest integer values
+ Float_t ncl_med = float(ncl-nused)/nTimeBins;
+ Int_t ncli = Int_t(ncl_med);
+ Float_t ncl_dev = TMath::Abs(ncl_med - TMath::Max(ncli, 1));
+ ncl_dev -= (ncl_dev>.5) && ncli ? .5 : 0.;
+ /*Double_t quality = */ return TMath::Exp(-2.*ncl_dev);
+
+// // 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);
+}
+
+//____________________________________________________________________
+Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed
+ , Int_t planes[4]
+ , Double_t *chi2)
+{
+ //
+ // Calculate the probability of this track candidate.
+ //
+ // Parameters :
+ // cseeds : array of candidate tracklets
+ // planes : array of seeding planes (see seeding configuration)
+ // chi2 : chi2 values (on the Z and Y direction) from the rieman fit of the track.
+ //
+ // Output :
+ // likelihood value
+ //
+ // Detailed description
+ //
+ // The track quality is estimated based on the following 4 criteria:
+ // 1. precision of the rieman fit on the Y direction (likea)
+ // 2. chi2 on the Y direction (likechi2y)
+ // 3. chi2 on the Z direction (likechi2z)
+ // 4. number of attached clusters compared to a reference value
+ // (see AliTRDrecoParam::fkFindable) (likeN)
+ //
+ // The distributions for each type of probabilities are given below as of
+ // (date). They have to be checked to assure consistency of estimation.
+ //
+
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ Int_t nTimeBins = cal->GetNumberOfTimeBins();
+ // ratio of the total number of clusters/track which are expected to be found by the tracker.
+ Float_t fgFindable = fRecoParam->GetFindableClusters();
+
+
+ Int_t nclusters = 0;
+ Double_t sumda = 0.;
+ for(UChar_t ilayer = 0; ilayer < 4; ilayer++){
+ Int_t jlayer = planes[ilayer];
+ nclusters += cseed[jlayer].GetN2();
+ sumda += TMath::Abs(cseed[jlayer].GetYfitR(1) - cseed[jlayer].GetYref(1));
+ }
+ Double_t likea = TMath::Exp(-sumda*10.6);
+ Double_t likechi2y = 0.0000000001;
+ if (chi2[1] < 0.5) likechi2y += TMath::Exp(-TMath::Sqrt(chi2[1]) * 7.73);
+ Double_t likechi2z = TMath::Exp(-chi2[0] * 0.088) / TMath::Exp(-chi2[0] * 0.019);
+ Int_t enc = Int_t(fgFindable*4.*nTimeBins); // Expected Number Of Clusters, normally 72
+ Double_t likeN = TMath::Exp(-(enc - nclusters) * 0.19);
+
+ Double_t like = likea * likechi2y * likechi2z * likeN;
+
+#ifdef DEBUG
+ //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){
+ TTreeSRedirector &cstreamer = *fDebugStreamer;
+ cstreamer << "CookLikelihood"
+ << "sumda=" << sumda
+ << "chi0=" << chi2[0]
+ << "chi1=" << chi2[1]
+ << "likea=" << likea
+ << "likechi2y=" << likechi2y
+ << "likechi2z=" << likechi2z
+ << "nclusters=" << nclusters
+ << "likeN=" << likeN
+ << "like=" << like
+ << "\n";
+ }
+#endif
+
+ return like;
+}
+
+//___________________________________________________________________
+void AliTRDtrackerV1::GetMeanCLStack(AliTRDstackLayer *layers
+ , Int_t *planes
+ , Double_t *params)
+{
+ //
+ // Determines the Mean number of clusters per layer.
+ // Needed to determine good Seeding Layers
+ //
+ // Parameters:
+ // - Array of AliTRDstackLayers
+ // - Container for the params
+ //
+ // Detailed description
+ //
+ // Two Iterations:
+ // In the first Iteration the mean is calculted using all layers.
+ // After this, all layers outside the 1-sigma-region are rejected.
+ // Then the mean value and the standard-deviation are calculted a second
+ // time in order to select all layers in the 1-sigma-region as good-candidates.
+ //
+
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ Int_t nTimeBins = cal->GetNumberOfTimeBins();
+
+ Float_t mean = 0, stdev = 0;
+ Double_t ncl[kNTimeBins*kNSeedPlanes], mcl[kNTimeBins*kNSeedPlanes];
+ Int_t position = 0;
+ memset(ncl, 0, sizeof(Int_t)*kNTimeBins*kNSeedPlanes);
+ memset(mcl, 0, sizeof(Int_t)*kNTimeBins*kNSeedPlanes);
+ Int_t nused = 0;
+ for(Int_t ipl = 0; ipl < kNSeedPlanes; ipl++){
+ for(Int_t ils = 0; ils < nTimeBins; ils++){
+ position = planes[ipl]*nTimeBins + ils;
+ ncl[ipl * nTimeBins + ils] = layers[position].GetNClusters();
+ nused = 0;
+ for(Int_t icl = 0; icl < ncl[ipl * nTimeBins + ils]; icl++)
+ if((layers[position].GetCluster(icl))->IsUsed()) nused++;
+ ncl[ipl * nTimeBins + ils] -= nused;
+ }
+ }
+ // Declaration of quartils:
+ //Double_t qvals[3] = {0.0, 0.0, 0.0};
+ //Double_t qprop[3] = {0.16667, 0.5, 0.83333};
+ // Iterations
+ Int_t counter;
+ Double_t *array;
+ Int_t *limit;
+ Int_t nLayers = nTimeBins * kNSeedPlanes;
+ for(Int_t iter = 0; iter < 2; iter++){
+ array = (iter == 0) ? &ncl[0] : &mcl[0];
+ limit = (iter == 0) ? &nLayers : &counter;
+ counter = 0;
+ if(iter == 1){
+ for(Int_t i = 0; i < nTimeBins *kNSeedPlanes; i++){
+ if((ncl[i] > mean + stdev) || (ncl[i] < mean - stdev)) continue; // Outside 1-sigma region
+// if((ncl[i] > qvals[2]) || (ncl[i] < qvals[0])) continue; // Outside 1-sigma region
+ if(ncl[i] == 0) continue; // 0-Layers also rejected
+ mcl[counter] = ncl[i];
+ counter++;
+ }
+ }
+ if(*limit == 0) break;
+ printf("Limit = %d\n", *limit);
+ //using quartils instead of mean and RMS
+// TMath::Quantiles(*limit,3,array,qvals,qprop,kFALSE);
+ mean = TMath::Median(*limit, array, 0x0);
+ stdev = TMath::RMS(*limit, array);
+ }
+// printf("Quantiles: 0.16667 = %3.3f, 0.5 = %3.3f, 0.83333 = %3.3f\n", qvals[0],qvals[1],qvals[2]);
+// memcpy(params,qvals,sizeof(Double_t)*3);
+ params[1] = (Double_t)TMath::Nint(mean);
+ params[0] = (Double_t)TMath::Nint(mean - stdev);
+ params[2] = (Double_t)TMath::Nint(mean + stdev);
+
+}
+
+//___________________________________________________________________
+Int_t AliTRDtrackerV1::GetSeedingLayers(AliTRDstackLayer *layers
+ , Double_t *params)
+{
+ //
+ // Algorithm to find optimal seeding layer
+ // Layers inside one sigma region (given by Quantiles) are sorted
+ // according to their difference.
+ // All layers outside are sorted according t
+ //
+ // Parameters:
+ // - Array of AliTRDstackLayers (in the current plane !!!)
+ // - Container for the Indices of the seeding Layer candidates
+ //
+ // Output:
+ // - Number of Layers inside the 1-sigma-region
+ //
+ // The optimal seeding layer should contain the mean number of
+ // custers in the layers in one chamber.
+ //
+
+ //printf("Params: %3.3f, %3.3f, %3.3f\n", params[0], params[1], params[2]);
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ const Int_t kMaxClustersLayer = AliTRDstackLayer::kMaxClustersLayer;
+ Int_t nTimeBins = cal->GetNumberOfTimeBins();
+ Int_t ncl[kNTimeBins], indices[kNTimeBins], bins[kMaxClustersLayer];
+ memset(ncl, 0, sizeof(Int_t)*kNTimeBins);
+ memset(indices, 0, sizeof(Int_t)*kNTimeBins);
+ memset(bins, 0, sizeof(Int_t)*kMaxClustersLayer);
+ Int_t nused = 0;
+ for(Int_t ils = 0; ils < nTimeBins; ils++){
+ ncl[ils] = layers[ils].GetNClusters();
+ nused = 0;
+ for(Int_t icl = 0; icl < ncl[ils]; icl++)
+ if((layers[ils].GetCluster(icl))->IsUsed()) nused++;
+ ncl[ils] -= nused;
+ }
+
+ Float_t mean = params[1];
+ for(Int_t ils = 0; ils < nTimeBins; ils++){
+ memmove(indices + bins[ncl[ils]+1] + 1, indices + bins[ncl[ils]+1], sizeof(Int_t)*(nTimeBins - ils));
+ indices[bins[ncl[ils]+1]] = ils;
+ for(Int_t i = ncl[ils]+1; i < kMaxClustersLayer; i++)
+ bins[i]++;
+ }
+
+ //for(Int_t i = 0; i < nTimeBins; i++) printf("Bin %d = %d\n", i, bins[i]);
+ Int_t sbin = -1;
+ Int_t nElements;
+ Int_t position = 0;
+ TRandom *r = new TRandom();
+ Int_t iter = 0;
+ while(1){
+ while(sbin < (Int_t)params[0] || sbin > (Int_t)params[2]){
+ // Randomly selecting one bin
+ sbin = (Int_t)r->Poisson(mean);
+ }
+ printf("Bin = %d\n",sbin);
+ //Randomly selecting one Layer in the bin
+ nElements = bins[sbin + 1] - bins[sbin];
+ printf("nElements = %d\n", nElements);
+ if(iter == 5){
+ position = (Int_t)(gRandom->Rndm()*(nTimeBins-1));
+ break;
+ }
+ else if(nElements==0){
+ iter++;
+ continue;
+ }
+ position = (Int_t)(gRandom->Rndm()*(nElements-1)) + bins[sbin];
+ break;
+ }
+ delete r;
+ return indices[position];
+}
+
+//____________________________________________________________________
+AliTRDcluster *AliTRDtrackerV1::FindSeedingCluster(AliTRDstackLayer *layers
+ , AliTRDseedV1/*AliRieman*/ *reference)
+{
+ //
+ // Finds a seeding Cluster for the extrapolation chamber.
+ //
+ // The seeding cluster should be as close as possible to the assumed
+ // track which is represented by a Rieman fit.
+ // Therefore the selecting criterion is the minimum distance between
+ // the best fitting cluster and the Reference which is derived from
+ // the AliTRDseed. Because all layers are assumed to be equally good
+ // a linear search is performed.
+ //
+ // Imput parameters: - layers: array of AliTRDstackLayers (in one chamber!!!)
+ // - sfit: the reference
+ //
+ // Output: - the best seeding cluster
+ //
+
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ Int_t nTimeBins = cal->GetNumberOfTimeBins();
+
+ // distances as squared distances
+ Int_t index = 0;
+ Float_t ypos = 0.0, zpos = 0.0, distance = 0.0, nearest_distance =100000.0;
+ ypos = reference->GetYref(0);
+ zpos = reference->GetZref(0);
+ AliTRDcluster *currentBest = 0x0, *temp = 0x0;
+ for(Int_t ils = 0; ils < nTimeBins; ils++){
+ // Reference positions
+// ypos = reference->GetYat(layers[ils].GetX());
+// zpos = reference->GetZat(layers[ils].GetX());
+ index = layers[ils].SearchNearestCluster(ypos, zpos, fRecoParam->GetRoad2y(), fRecoParam->GetRoad2z());
+ if(index == -1) continue;
+ temp = layers[ils].GetCluster(index);
+ if(!temp) continue;
+ distance = (temp->GetY() - ypos) * (temp->GetY() - ypos) + (temp->GetZ() - zpos) * (temp->GetZ() - zpos);
+ if(distance < nearest_distance){
+ nearest_distance = distance;
+ currentBest = temp;
+ }
+ }
+ return currentBest;
+}
+
+//____________________________________________________________________
+AliTRDstackLayer *AliTRDtrackerV1::MakeSeedingLayer(AliTRDstackLayer *layers
+ , Int_t plane)
+{
+ //
+ // Creates a seeding layer
+ //
+
+ // constants
+ const Int_t kMaxRows = 16;
+ const Int_t kMaxCols = 144;
+ const Int_t kMaxPads = 2304;
+
+ // Get the calculation
+ AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+ Int_t nTimeBins = cal->GetNumberOfTimeBins();
+
+ // Get the geometrical data of the chamber
+ AliTRDpadPlane *pp = fGeom->GetPadPlane(plane, layers[0].GetStackNr());
+ 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());
+ 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 arrpos;
+ Float_t ypos;
+ Int_t irow, 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 < nTimeBins; iTime++){
+ nClusters = layers[iTime].GetNClusters();
+ for(Int_t incl = 0; incl < nClusters; incl++){
+ c = layers[iTime].GetCluster(incl);
+ ypos = c->GetY();
+ if(ypos > ymax && ypos < ymin) continue;
+ irow = pp->GetPadRowNumber(c->GetZ()); // Zbin
+ if(irow < 0)continue;
+ arrpos = static_cast<Int_t>((ypos - ymin)/binlength);
+ if(ypos == ymax) arrpos = nCols - 1;
+ histogram[irow][arrpos]++;
+ sigmas[irow][arrpos] += c->GetSigmaZ2();
+ }
+ }
+
+// Now I have everything in the histogram, do the selection
+// printf("Starting the analysis\n");
+ //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 ti 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[kMaxTracksStack];
+ Float_t sigcands[kMaxTracksStack];
+
+ // helper variables
+ Int_t indices[kMaxPads]; memset(indices, 0, sizeof(Int_t)*kMaxPads);
+ Int_t nCandidates = 0;
+ Float_t norm, cogv;
+ // histogram filled -> Select best bins
+ TMath::Sort(kMaxPads, hvals, indices); // bins storing a 0 should not matter
+ // Set Threshold
+ Int_t maximum = hvals[indices[0]]; // best
+ Int_t threshold = static_cast<UChar_t>(maximum * fRecoParam->GetFindableClusters());
+ Int_t col, row, lower, lower1, upper, upper1;
+ for(Int_t ib = 0; ib < kMaxPads; ib++){
+ if(nCandidates >= kMaxTracksStack){
+ AliWarning(Form("Number of seed candidates %d exceeded maximum allowed per stack %d", nCandidates, 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] = cogv / norm;
+ // passed the filter
+ cand[nCandidates] = row*kMaxCols + 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++;
+ }
+ AliTRDstackLayer *fakeLayer = new AliTRDstackLayer(layers[0].GetZ0(), layers[0].GetDZ0(), layers[0].GetStackNr());
+ fakeLayer->SetX((TMath::Abs(layers[nTimeBins-1].GetX() + layers[0].GetX()))/2);
+ fakeLayer->SetSector(layers[0].GetSector());
+ AliTRDcluster **fakeClusters = 0x0;
+ UInt_t *fakeIndices = 0x0;
+ if(nCandidates){
+ fakeClusters = new AliTRDcluster*[nCandidates];
+ fakeIndices = new UInt_t[nCandidates];
+ UInt_t fakeIndex = 0;
+ for(Int_t ican = 0; ican < nCandidates; ican++){
+ fakeClusters[ican] = new AliTRDcluster();
+ fakeClusters[ican]->SetX(fakeLayer->GetX());
+ fakeClusters[ican]->SetY(cogyvals[cand[ican]]);
+ fakeClusters[ican]->SetZ(cogzvals[cand[ican]]);
+ fakeClusters[ican]->SetSigmaZ2(sigcands[ican]);
+ fakeIndices[ican] = fakeIndex++;// fantasy number
+ }
+ }
+ fakeLayer->SetRecoParam(fRecoParam);
+ fakeLayer->SetClustersArray(fakeClusters, nCandidates);
+ fakeLayer->SetIndexArray(fakeIndices);
+ fakeLayer->SetNRows(nRows);
+ fakeLayer->BuildIndices();
+ //fakeLayer->PrintClusters();
+
+#ifdef DEBUG
+ if(AliTRDReconstructor::StreamLevel() >= 3){
+ TMatrixT<double> 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 = *fDebugStreamer;
+ cstreamer << "MakeSeedingLayer"
+ << "Iteration=" << fSieveSeeding
+ << "plane=" << plane
+ << "ymin=" << ymin
+ << "ymax=" << ymax
+ << "zmin=" << zmin
+ << "zmax=" << zmax
+ << "L.=" << fakeLayer
+ << "Histogram.=" << &hist
+ << "\n";
+ }
+#endif
+ return fakeLayer;
+}
+
+//____________________________________________________________________
+void AliTRDtrackerV1::GetSeedingConfig(Int_t iconfig, Int_t planes[4])
+{
+ //
+ // Map seeding configurations to detector planes.
+ //
+ // Parameters :
+ // iconfig : configuration index
+ // planes : member planes of this configuration. On input empty.
+ //
+ // Output :
+ // planes : contains the planes which are defining the configuration
+ //
+ // Detailed description
+ //
+ // Here is the list of seeding planes configurations together with
+ // their topological classification:
+ //
+ // 0 - 5432 TQ 0
+ // 1 - 4321 TQ 0
+ // 2 - 3210 TQ 0
+ // 3 - 5321 TQ 1
+ // 4 - 4210 TQ 1
+ // 5 - 5431 TQ 1
+ // 6 - 4320 TQ 1
+ // 7 - 5430 TQ 2
+ // 8 - 5210 TQ 2
+ // 9 - 5421 TQ 3
+ // 10 - 4310 TQ 3
+ // 11 - 5410 TQ 4
+ // 12 - 5420 TQ 5
+ // 13 - 5320 TQ 5
+ // 14 - 5310 TQ 5
+ //
+ // The topologic quality is modeled as follows:
+ // 1. The general model is define by the equation:
+ // p(conf) = exp(-conf/2)
+ // 2. According to the topologic classification, configurations from the same
+ // class are assigned the agerage value over the model values.
+ // 3. Quality values are normalized.
+ //
+ // The topologic quality distribution as function of configuration is given below:
+ //Begin_Html
+ // <img src="gif/topologicQA.gif">
+ //End_Html
+ //
+
+ switch(iconfig){
+ case 0: // 5432 TQ 0
+ planes[0] = 2;
+ planes[1] = 3;
+ planes[2] = 4;
+ planes[3] = 5;
+ break;
+ case 1: // 4321 TQ 0
+ planes[0] = 1;
+ planes[1] = 2;
+ planes[2] = 3;
+ planes[3] = 4;
+ break;
+ case 2: // 3210 TQ 0
+ planes[0] = 0;
+ planes[1] = 1;
+ planes[2] = 2;
+ planes[3] = 3;
+ break;
+ case 3: // 5321 TQ 1
+ planes[0] = 1;
+ planes[1] = 2;
+ planes[2] = 3;
+ planes[3] = 5;
+ break;
+ case 4: // 4210 TQ 1
+ planes[0] = 0;
+ planes[1] = 1;
+ planes[2] = 2;
+ planes[3] = 4;
+ break;
+ case 5: // 5431 TQ 1
+ planes[0] = 1;
+ planes[1] = 3;
+ planes[2] = 4;
+ planes[3] = 5;
+ break;
+ case 6: // 4320 TQ 1
+ planes[0] = 0;
+ planes[1] = 2;
+ planes[2] = 3;
+ planes[3] = 4;
+ break;
+ case 7: // 5430 TQ 2
+ planes[0] = 0;
+ planes[1] = 3;
+ planes[2] = 4;
+ planes[3] = 5;
+ break;
+ case 8: // 5210 TQ 2
+ planes[0] = 0;
+ planes[1] = 1;
+ planes[2] = 2;
+ planes[3] = 5;
+ break;
+ case 9: // 5421 TQ 3
+ planes[0] = 1;
+ planes[1] = 2;
+ planes[2] = 4;
+ planes[3] = 5;
+ break;
+ case 10: // 4310 TQ 3
+ planes[0] = 0;
+ planes[1] = 1;
+ planes[2] = 3;
+ planes[3] = 4;
+ break;
+ case 11: // 5410 TQ 4
+ planes[0] = 0;
+ planes[1] = 1;
+ planes[2] = 4;
+ planes[3] = 5;
+ break;
+ case 12: // 5420 TQ 5
+ planes[0] = 0;
+ planes[1] = 2;
+ planes[2] = 4;
+ planes[3] = 5;
+ break;
+ case 13: // 5320 TQ 5
+ planes[0] = 0;
+ planes[1] = 2;
+ planes[2] = 3;
+ planes[3] = 5;
+ break;
+ case 14: // 5310 TQ 5
+ planes[0] = 0;
+ planes[1] = 1;
+ planes[2] = 3;
+ planes[3] = 5;
+ break;
+ }
+}
+
+//____________________________________________________________________
+void AliTRDtrackerV1::GetExtrapolationConfig(Int_t iconfig, Int_t planes[2])
+{
+ //
+ // Returns the extrapolation planes for a seeding configuration.
+ //
+ // Parameters :
+ // iconfig : configuration index
+ // planes : planes which are not in this configuration. On input empty.
+ //
+ // Output :
+ // planes : contains the planes which are not in the configuration
+ //
+ // Detailed description
+ //
+
+ switch(iconfig){
+ case 0: // 5432 TQ 0
+ planes[0] = 1;
+ planes[1] = 0;
+ break;
+ case 1: // 4321 TQ 0
+ planes[0] = 5;
+ planes[1] = 0;
+ break;
+ case 2: // 3210 TQ 0
+ planes[0] = 4;
+ planes[1] = 5;
+ break;
+ case 3: // 5321 TQ 1
+ planes[0] = 4;
+ planes[1] = 0;
+ break;
+ case 4: // 4210 TQ 1
+ planes[0] = 5;
+ planes[1] = 3;
+ break;
+ case 5: // 5431 TQ 1
+ planes[0] = 2;
+ planes[1] = 0;
+ break;
+ case 6: // 4320 TQ 1
+ planes[0] = 5;
+ planes[1] = 1;
+ break;
+ case 7: // 5430 TQ 2
+ planes[0] = 2;
+ planes[1] = 1;
+ break;
+ case 8: // 5210 TQ 2
+ planes[0] = 4;
+ planes[1] = 3;
+ break;
+ case 9: // 5421 TQ 3
+ planes[0] = 3;
+ planes[1] = 0;
+ break;
+ case 10: // 4310 TQ 3
+ planes[0] = 5;
+ planes[1] = 2;
+ break;
+ case 11: // 5410 TQ 4
+ planes[0] = 3;
+ planes[1] = 2;
+ break;
+ case 12: // 5420 TQ 5
+ planes[0] = 3;
+ planes[1] = 1;
+ break;
+ case 13: // 5320 TQ 5
+ planes[0] = 4;
+ planes[1] = 1;
+ break;
+ case 14: // 5310 TQ 5
+ planes[0] = 4;
+ planes[1] = 2;
+ break;
+ }
+}
--- /dev/null
+#ifndef ALITRDTRACKERV1_H
+#define ALITRDTRACKERV1_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+
+////////////////////////////////////////////////////////////////////////////
+// //
+// The TRD tracker //
+// //
+// Authors: //
+// Alex Bercuci <A.Bercuci@gsi.de> //
+// Markus Fasel <M.Fasel@gsi.de> //
+// //
+/////////////////////////////////////////////////////////////////////////////
+
+#ifndef ALITRDTRACKER_H
+#include "AliTRDtracker.h"
+#endif
+
+#define DEBUG
+
+/**************************************************************************
+ * Class Status see source file *
+ **************************************************************************/
+
+class TFile;
+class TTreeSRedirector;
+class TClonesArray;
+
+class AliRieman;
+class AliESDEvent;
+
+class AliTRDseedV1;
+class AliTRDstackLayer;
+class AliTRDtrackerFitter;
+class AliTRDrecoParam;
+
+class AliTRDtrackerV1 : public AliTRDtracker
+{
+
+ public:
+ enum{
+ kNTimeBins = 35,
+ kNPlanes = 6,
+ kNSeedPlanes = 4,
+ kMaxTracksStack = 100,
+ kNConfigs = 15
+ };
+ AliTRDtrackerV1(AliTRDrecoParam *p = 0x0);
+ AliTRDtrackerV1(const TFile *in, AliTRDrecoParam *p);
+ ~AliTRDtrackerV1();
+
+ Int_t Clusters2Tracks(AliESDEvent *esd);
+ void GetSeedingConfig(Int_t iconfig, Int_t planes[4]);
+ void GetExtrapolationConfig(Int_t iconfig, Int_t planes[2]);
+ void SetRecoParam(AliTRDrecoParam *p){fRecoParam = p;}
+
+ protected:
+
+ Double_t BuildSeedingConfigs(AliTRDstackLayer *layer, Int_t *configs);
+ Int_t Clusters2TracksSM(AliTRDtracker::AliTRDtrackingSector *sector, AliESDEvent *esd);
+ Int_t Clusters2TracksStack(AliTRDstackLayer *layer, TClonesArray *esdTrackList);
+ Double_t CookPlaneQA(AliTRDstackLayer *layer);
+ Double_t CookLikelihood(AliTRDseedV1 *cseed, Int_t planes[4], Double_t *chi2);
+ Int_t GetSeedingLayers(AliTRDstackLayer *layers, Double_t *params);
+ void GetMeanCLStack(AliTRDstackLayer *layers, Int_t *planes, Double_t *params);
+ AliTRDcluster *FindSeedingCluster(AliTRDstackLayer *layers, AliTRDseedV1/*AliRieman*/ *sfit);
+ void ImproveSeedQuality(AliTRDstackLayer *layer, AliTRDseedV1 *cseed);
+ Int_t MakeSeeds(AliTRDstackLayer *layers, AliTRDseedV1 *sseed, Int_t *ipar);
+ AliTRDstackLayer *MakeSeedingLayer(AliTRDstackLayer *layers, Int_t Plane);
+ AliTRDtrack* RegisterSeed(AliTRDseedV1 *seeds, Double_t *params);
+
+ private:
+
+ AliTRDtrackerV1(const AliTRDtrackerV1 &tracker);
+ AliTRDtrackerV1 &operator=(const AliTRDtrackerV1 &tracker);
+
+ private:
+
+ static Double_t fTopologicQA[kNConfigs]; // Topologic quality
+ Double_t fTrackQuality[kMaxTracksStack]; // Track quality
+ Int_t fSeedLayer[kMaxTracksStack]; // Seed layer
+ Int_t fSieveSeeding; //! Seeding iterator
+ AliTRDrecoParam *fRecoParam; // Reconstruction parameters
+ AliTRDtrackerFitter *fFitter; //! Fitter class of the tracker
+ TTreeSRedirector *fDebugStreamer; //! Debug stream of the tracker
+
+ ClassDef(AliTRDtrackerV1, 1) // Stand alone tracker development class
+
+};
+#endif
#pragma link C++ class AliTRDtrack+;
#pragma link C++ class AliTRDtracklet+;
#pragma link C++ class AliTRDtracker+;
+#pragma link C++ class AliTRDpropagationLayer+;
#pragma link C++ class AliTRDseed+;
#pragma link C++ class AliTRDpidESD+;
#pragma link C++ class AliTRDtrackingAnalysis+;
+#pragma link C++ class AliTRDrecoParam+;
+#pragma link C++ class AliTRDseedV1+;
+#pragma link C++ class AliTRDtrackerV1+;
+#pragma link C++ class AliTRDstackLayer+;
+#pragma link C++ class AliTRDtrackerFitter+;
+
+
#endif
AliTRDpidESD.cxx \
AliTRDRecParam.cxx \
AliTRDReconstructor.cxx \
- AliTRDtrackingAnalysis.cxx
+ AliTRDtrackingAnalysis.cxx \
+ AliTRDrecoParam.cxx \
+ AliTRDseedV1.cxx \
+ AliTRDtrackerV1.cxx \
+ AliTRDpropagationLayer.cxx \
+ AliTRDstackLayer.cxx \
+ AliTRDtrackerFitter.cxx
-HDRS= $(SRCS:.cxx=.h)
+HDRS= $(SRCS:.cxx=.h)
DHDR= TRDrecLinkDef.h