From 5443e65ea0e4ba7bbf3dd08d9893a31c9dd8a539 Mon Sep 17 00:00:00 2001 From: cblume Date: Wed, 12 Jun 2002 09:54:36 +0000 Subject: [PATCH] Update of tracking code provided by Sergei --- TRD/AliTRD.cxx | 29 +- TRD/AliTRD.h | 5 +- TRD/AliTRDcluster.cxx | 115 +- TRD/AliTRDcluster.h | 84 +- TRD/AliTRDclusterizer.cxx | 7 + TRD/AliTRDclusterizer.h | 37 +- TRD/AliTRDclusterizerV0.cxx | 39 +- TRD/AliTRDclusterizerV1.cxx | 49 +- TRD/AliTRDclusterizerV1.h | 7 +- TRD/AliTRDgeometry.cxx | 185 +-- TRD/AliTRDgeometry.h | 58 +- TRD/AliTRDgeometryDetail.cxx | 24 +- TRD/AliTRDparameter.cxx | 6 +- TRD/AliTRDpid.cxx | 17 +- TRD/AliTRDrecPoint.cxx | 64 +- TRD/AliTRDsimpleMC.cxx | 15 +- TRD/AliTRDsimpleMC.h | 44 +- TRD/AliTRDtrack.cxx | 291 ++-- TRD/AliTRDtrack.h | 75 +- TRD/AliTRDtracker.cxx | 2534 ++++++++++++++++++++++++++-------- TRD/AliTRDtracker.h | 271 +++- TRD/AliTRDtrackingSector.cxx | 32 +- TRD/AliTRDtrackingSector.h | 7 +- TRD/Makefile | 2 - TRD/TRDLinkDef.h | 4 +- 25 files changed, 2667 insertions(+), 1334 deletions(-) diff --git a/TRD/AliTRD.cxx b/TRD/AliTRD.cxx index 4762f1f5310..94971d745fc 100644 --- a/TRD/AliTRD.cxx +++ b/TRD/AliTRD.cxx @@ -15,6 +15,9 @@ /* $Log$ +Revision 1.38 2002/03/28 14:59:07 cblume +Coding conventions + Revision 1.37 2002/03/25 20:01:49 cblume Introduce parameter class @@ -300,37 +303,23 @@ AliTRD::~AliTRD() } //_____________________________________________________________________________ -void AliTRD::AddCluster(Float_t *pos, Int_t *digits, Int_t det, Float_t amp - , Int_t *tracks, Float_t sigmaY2, Int_t iType) +void AliTRD::AddCluster(Float_t *pos, Int_t det, Float_t amp + , Int_t *tracks, Float_t *sig, Int_t iType) { // // Add a cluster for the TRD // - Int_t pla = fGeometry->GetPlane(det); - Int_t cha = fGeometry->GetChamber(det); - Int_t sec = fGeometry->GetSector(det); - - Float_t padRow = pos[0]; - Float_t padCol = pos[1]; - Float_t row0 = fGeometry->GetRow0(pla,cha,sec); - Float_t col0 = fGeometry->GetCol0(pla); - Float_t rowSize = fGeometry->GetRowPadSize(pla,cha,sec); - Float_t colSize = fGeometry->GetColPadSize(pla); - AliTRDcluster *c = new AliTRDcluster(); c->SetDetector(det); c->AddTrackIndex(tracks); - c->SetQ(amp); - + c->SetY(pos[0]); + c->SetZ(pos[1]); + c->SetSigmaY2(sig[0]); + c->SetSigmaZ2(sig[1]); c->SetLocalTimeBin(((Int_t) pos[2])); - c->SetY(col0 + padCol * colSize); - c->SetZ(row0 + padRow * rowSize); - - c->SetSigmaY2((sigmaY2 + 1./12.) * colSize*colSize); - c->SetSigmaZ2(rowSize * rowSize / 12.); switch (iType) { case 0: diff --git a/TRD/AliTRD.h b/TRD/AliTRD.h index f8796b1b873..e6ee05be339 100644 --- a/TRD/AliTRD.h +++ b/TRD/AliTRD.h @@ -34,9 +34,8 @@ class AliTRD : public AliDetector { AliTRD &operator=(const AliTRD &trd); virtual void AddHit(Int_t track, Int_t det, Float_t *hits, Int_t q, Bool_t inDrift); - virtual void AddCluster(Float_t *pos, Int_t *digits - , Int_t det, Float_t amp, Int_t *tracks - , Float_t sigmaY2, Int_t iType); + virtual void AddCluster(Float_t *pos, Int_t det, Float_t amp, Int_t *tracks + , Float_t *sig, Int_t iType); virtual void BuildGeometry(); virtual void Copy(TObject &trd); virtual void CreateGeometry(); diff --git a/TRD/AliTRDcluster.cxx b/TRD/AliTRDcluster.cxx index 2b0191652a5..5bfc692e6e4 100644 --- a/TRD/AliTRDcluster.cxx +++ b/TRD/AliTRDcluster.cxx @@ -28,40 +28,16 @@ Revision 1.1.2.1 2000/09/22 14:47:52 cblume Add the tracking code */ - -/////////////////////////////////////////////////////////////////////////////// -// // -// TRD cluster // -// // -/////////////////////////////////////////////////////////////////////////////// #include "AliTRDcluster.h" +#include "AliTRDgeometry.h" #include "AliTRDrecPoint.h" ClassImp(AliTRDcluster) -//_____________________________________________________________________________ -AliTRDcluster::AliTRDcluster() -{ - // - // Default constructor - // - - fDetector = 0; - fTimeBin = 0; - fTracks[0] = 0; - fTracks[1] = 0; - fTracks[2] = 0; - fY = 0; - fZ = 0; - fQ = 0; - fSigmaY2 = 0; - fSigmaZ2 = 0; - -} //_____________________________________________________________________________ -AliTRDcluster::AliTRDcluster(const AliTRDrecPoint &p) + AliTRDcluster::AliTRDcluster(const AliTRDrecPoint &p):AliCluster() { // // Constructor from AliTRDrecPoint @@ -87,61 +63,27 @@ AliTRDcluster::AliTRDcluster(const AliTRDrecPoint &p) } //_____________________________________________________________________________ -AliTRDcluster::AliTRDcluster(const AliTRDcluster &c) +AliTRDcluster::AliTRDcluster(const AliTRDcluster &c):AliCluster() { // // Copy constructor // - ((AliTRDcluster &) c).Copy(*this); - -} - -//_____________________________________________________________________________ -AliTRDcluster::~AliTRDcluster() -{ - // - // AliTRDcluster destructor - // - -} - -//_____________________________________________________________________________ -AliTRDcluster &AliTRDcluster::operator=(const AliTRDcluster &c) -{ - // - // Assignment operator - // - - if (this != &c) ((AliTRDcluster &) c).Copy(*this); - return *this; + fTracks[0] = c.GetLabel(0); + fTracks[1] = c.GetLabel(1); + fTracks[2] = c.GetLabel(2); -} + fY = c.GetY(); + fZ = c.GetZ(); + fSigmaY2 = c.GetSigmaY2(); + fSigmaZ2 = c.GetSigmaZ2(); -//_____________________________________________________________________________ -void AliTRDcluster::Copy(TObject &c) -{ - // - // Copy function - // - - ((AliTRDcluster &) c).fDetector = fDetector; - ((AliTRDcluster &) c).fTimeBin = fTimeBin; - - ((AliTRDcluster &) c).fTracks[0] = fTracks[0]; - ((AliTRDcluster &) c).fTracks[1] = fTracks[1]; - ((AliTRDcluster &) c).fTracks[2] = fTracks[2]; - - ((AliTRDcluster &) c).fQ = fQ; - - ((AliTRDcluster &) c).fY = fY; - ((AliTRDcluster &) c).fZ = fZ; - ((AliTRDcluster &) c).fSigmaY2 = fSigmaY2; - ((AliTRDcluster &) c).fSigmaZ2 = fSigmaZ2; + fDetector = c.GetDetector(); + fTimeBin = c.GetLocalTimeBin(); + fQ = c.GetQ(); } - //_____________________________________________________________________________ void AliTRDcluster::AddTrackIndex(Int_t *track) { @@ -154,38 +96,37 @@ void AliTRDcluster::AddTrackIndex(Int_t *track) // ones are stored first; // - const Int_t kSize = 9; + const Int_t size = 9; - Int_t entries[kSize][2], i, j, index; + Int_t entries[size][2], i, j, index; - Bool_t indexAdded; + Bool_t index_added; - for (i=0; i= 0) { - while ( (!indexAdded) && ( j < kSize ) ) { + while ( (!index_added) && ( j < size ) ) { if ((entries[j][0]==index) || (entries[j][1]==0)) { entries[j][0]=index; entries[j][1]=entries[j][1]+1; - indexAdded=kTRUE; + index_added=kTRUE; } j++; } } - } + } // sort by number of appearances and index value Int_t swap=1, tmp0, tmp1; while ( swap > 0) { swap=0; - for(i=0; i<(kSize-1); i++) { + for(i=0; i<(size-1); i++) { if ((entries[i][0] >= 0) && (entries[i+1][0] >= 0)) { if ((entries[i][1] < entries[i+1][1]) || ((entries[i][1] == entries[i+1][1]) && @@ -200,14 +141,12 @@ void AliTRDcluster::AddTrackIndex(Int_t *track) } } } - } + } // set track indexes - for(i=0; i<3; i++) { - fTracks[i] = entries[i][0]; - } + for(i=0; i<3; i++) SetLabel(entries[i][0],i); return; -} +} diff --git a/TRD/AliTRDcluster.h b/TRD/AliTRDcluster.h index cc16efc9216..451df071615 100644 --- a/TRD/AliTRDcluster.h +++ b/TRD/AliTRDcluster.h @@ -5,64 +5,46 @@ /* $Id$ */ -#include - -/////////////////////////////////////////////////////////////////////////////// -// // -// TRD cluster // -// // -/////////////////////////////////////////////////////////////////////////////// + +#include "AliCluster.h" class AliTRDrecPoint; -class AliTRDcluster : public TObject { +class AliTRDcluster : public AliCluster { public: - AliTRDcluster(); + AliTRDcluster() : AliCluster() { fQ=0; fTimeBin=0; fDetector=0; } AliTRDcluster(const AliTRDcluster &c); AliTRDcluster(const AliTRDrecPoint &p); - virtual ~AliTRDcluster(); - AliTRDcluster &operator=(const AliTRDcluster &c); - virtual void Copy(TObject &c); - virtual void AddTrackIndex(Int_t *i); - - Int_t GetDetector() const { return fDetector; } - Int_t GetLocalTimeBin() const { return fTimeBin; } + virtual void AddTrackIndex(Int_t *i); - Float_t GetSigmaY2() const { return fSigmaY2; } - Float_t GetSigmaZ2() const { return fSigmaZ2; } - Float_t GetY() const { return fY; } - Float_t GetZ() const { return fZ; } - Float_t GetQ() const { return fQ; } - Int_t IsUsed() const { return (fQ < 0) ? 1 : 0; } - void Use() { fQ = -fQ; } - Int_t GetTrackIndex(Int_t i) const { return fTracks[i]; } - - Bool_t From2pad() const { return TestBit(k2pad); } - Bool_t From3pad() const { return TestBit(k3pad); } - Bool_t From4pad() const { return TestBit(k4pad); } - Bool_t From5pad() const { return TestBit(k5pad); } - Bool_t FromLarge() const { return TestBit(kLarge); } - Bool_t Isolated() const { return (TestBit(k2pad) || TestBit(k3pad)); } + Int_t IsUsed() const { return (fQ < 0) ? 1 : 0; } + void Use() { fQ = -fQ; } + + Bool_t From2pad() const { return TestBit(k2pad); } + Bool_t From3pad() const { return TestBit(k3pad); } + Bool_t From4pad() const { return TestBit(k4pad); } + Bool_t From5pad() const { return TestBit(k5pad); } + Bool_t FromLarge() const { return TestBit(kLarge);} + Bool_t Isolated() const { return (TestBit(k2pad) || TestBit(k3pad)); } - void SetDetector(Int_t d) { fDetector = d; } - void SetLocalTimeBin(Int_t t) { fTimeBin = t; } - void SetQ(Float_t q) { fQ = q; } - void SetY(Float_t y) { fY = y; } - void SetZ(Float_t z) { fZ = z; } - void SetTrackIndex(Int_t i, Int_t t) { fTracks[i] = t; } - void SetSigmaY2(Float_t s) { fSigmaY2 = s; } - void SetSigmaZ2(Float_t s) { fSigmaZ2 = s; } - - void Set2pad() { SetBit(k2pad); } - void Set3pad() { SetBit(k3pad); } - void Set4pad() { SetBit(k4pad); } - void Set5pad() { SetBit(k5pad); } - void SetLarge() { SetBit(kLarge); } - + void SetDetector(Int_t d) { fDetector = d; } + void SetLocalTimeBin(Int_t t) { fTimeBin = t; } + void SetQ(Float_t q) { fQ = q; } + + Int_t GetDetector() const { return fDetector; } + Int_t GetLocalTimeBin() const { return fTimeBin; } + Float_t GetQ() const { return fQ; } + + void Set2pad() { SetBit(k2pad); } + void Set3pad() { SetBit(k3pad); } + void Set4pad() { SetBit(k4pad); } + void Set5pad() { SetBit(k5pad); } + void SetLarge() { SetBit(kLarge); } + protected: enum { @@ -72,17 +54,11 @@ class AliTRDcluster : public TObject { k5pad = 0x00000008, // 5 pad cluster kLarge = 0x00000016 // Large cluster }; - + Int_t fDetector; // TRD detector number Int_t fTimeBin; // Time bin number within the detector - - Int_t fTracks[3]; // labels of overlapped tracks Float_t fQ; // amplitude - Float_t fY; // local Rphi coordinate (cm) within tracking sector - Float_t fZ; // local Z coordinate (cm) within tracking sector - Float_t fSigmaY2; // Y variance (cm) - Float_t fSigmaZ2; // Z variance (cm) - + ClassDef(AliTRDcluster,1) // Cluster for the TRD }; diff --git a/TRD/AliTRDclusterizer.cxx b/TRD/AliTRDclusterizer.cxx index d185bb79c88..ce1d21b3894 100644 --- a/TRD/AliTRDclusterizer.cxx +++ b/TRD/AliTRDclusterizer.cxx @@ -15,6 +15,9 @@ /* $Log$ +Revision 1.11 2001/11/27 08:50:33 hristov +BranchOld replaced by Branch + Revision 1.10 2001/11/14 10:50:45 cblume Changes in digits IO. Add merging of summable digits @@ -91,6 +94,7 @@ Add new TRD classes #include "AliTRDcluster.h" #include "AliTRDrecPoint.h" #include "AliTRDgeometry.h" +#include "AliTRDparameter.h" ClassImp(AliTRDclusterizer) @@ -107,6 +111,7 @@ AliTRDclusterizer::AliTRDclusterizer():TNamed() fTRD = 0; fEvent = 0; fVerbose = 0; + fPar = 0; } @@ -123,6 +128,7 @@ AliTRDclusterizer::AliTRDclusterizer(const Text_t* name, const Text_t* title) fClusterTree = NULL; fEvent = 0; fVerbose = 0; + fPar = 0; } @@ -184,6 +190,7 @@ void AliTRDclusterizer::Copy(TObject &c) ((AliTRDclusterizer &) c).fClusterTree = NULL; ((AliTRDclusterizer &) c).fEvent = 0; ((AliTRDclusterizer &) c).fVerbose = fVerbose; + ((AliTRDclusterizer &) c).fPar = 0; } diff --git a/TRD/AliTRDclusterizer.h b/TRD/AliTRDclusterizer.h index 94a1532d19a..451d59fee72 100644 --- a/TRD/AliTRDclusterizer.h +++ b/TRD/AliTRDclusterizer.h @@ -9,6 +9,8 @@ class TFile; +class AliTRDparameter; + /////////////////////////////////////////////////////// // Finds and handles cluster // /////////////////////////////////////////////////////// @@ -23,27 +25,32 @@ class AliTRDclusterizer : public TNamed { virtual ~AliTRDclusterizer(); AliTRDclusterizer &operator=(const AliTRDclusterizer &c); - virtual void Copy(TObject &c); - virtual Bool_t Open(const Char_t *name, Int_t nEvent = 0); - virtual Bool_t Open(const Char_t *inname, const Char_t *outname, Int_t nEvent = 0); - virtual Bool_t OpenInput(const Char_t *name, Int_t nEvent = 0); - virtual Bool_t OpenOutput(const Char_t *name); - virtual Bool_t MakeClusters() = 0; - virtual Bool_t WriteClusters(Int_t det); + virtual void Copy(TObject &c); + virtual Bool_t Open(const Char_t *name, Int_t nEvent = 0); + virtual Bool_t Open(const Char_t *inname, const Char_t *outname, Int_t nEvent = 0); + virtual Bool_t OpenInput(const Char_t *name, Int_t nEvent = 0); + virtual Bool_t OpenOutput(const Char_t *name); + virtual Bool_t MakeClusters() = 0; + virtual Bool_t WriteClusters(Int_t det); + + void SetVerbose(Int_t v = 1) { fVerbose = v; }; + + virtual void SetParameter(AliTRDparameter *par) { fPar = par; }; - void SetVerbose(Int_t v = 1) { fVerbose = v; }; + AliTRDparameter *GetParameter() const { return fPar; }; protected: - TFile *fInputFile; //! AliROOT input file - TFile *fOutputFile; //! AliROOT output file - TTree *fClusterTree; //! Tree with the cluster - AliTRD *fTRD; //! The TRD object + TFile *fInputFile; //! AliROOT input file + TFile *fOutputFile; //! AliROOT output file + TTree *fClusterTree; //! Tree with the cluster + AliTRD *fTRD; //! The TRD object + AliTRDparameter *fPar; // TRD digitization parameter object - Int_t fEvent; //! Event number - Int_t fVerbose; // Sets the verbose level + Int_t fEvent; //! Event number + Int_t fVerbose; // Sets the verbose level - ClassDef(AliTRDclusterizer,2) // TRD-Cluster manager base class + ClassDef(AliTRDclusterizer,3) // TRD-Cluster manager base class }; diff --git a/TRD/AliTRDclusterizerV0.cxx b/TRD/AliTRDclusterizerV0.cxx index 9e28834f38c..178e5e7a347 100644 --- a/TRD/AliTRDclusterizerV0.cxx +++ b/TRD/AliTRDclusterizerV0.cxx @@ -15,6 +15,9 @@ /* $Log$ +Revision 1.10 2001/12/05 15:04:34 hristov +Changes related to the corrections of AliRecPoint + Revision 1.9 2001/05/28 17:07:58 hristov Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of material in G10 TEC frames and material between TEC planes (C.Blume,S.Sedykh) @@ -98,6 +101,7 @@ Add new TRD classes #include "AliTRDhit.h" #include "AliTRDgeometry.h" #include "AliTRDrecPoint.h" +#include "AliTRDparameter.h" ClassImp(AliTRDclusterizerV0) @@ -160,6 +164,15 @@ Bool_t AliTRDclusterizerV0::MakeClusters() // Get the geometry AliTRDgeometry *geo = fTRD->GetGeometry(); + + // Create a default parameter class if none is defined + if (!fPar) { + fPar = new AliTRDparameter("TRDparameter","Standard TRD parameter"); + if (fVerbose > 0) { + printf(" "); + printf("Create the default parameter object.\n"); + } + } printf("AliTRDclusterizerV0::MakeCluster -- "); printf("Start creating cluster.\n"); @@ -184,14 +197,14 @@ Bool_t AliTRDclusterizerV0::MakeClusters() for (Int_t iplan = 0; iplan < AliTRDgeometry::Nplan(); iplan++) { for (Int_t isect = 0; isect < AliTRDgeometry::Nsect(); isect++) { - Int_t nColMax = geo->GetColMax(iplan); - Float_t row0 = geo->GetRow0(iplan,icham,isect); - Float_t col0 = geo->GetCol0(iplan); - Float_t time0 = geo->GetTime0(iplan); + Int_t nColMax = fPar->GetColMax(iplan); + Float_t row0 = fPar->GetRow0(iplan,icham,isect); + Float_t col0 = fPar->GetCol0(iplan); + Float_t time0 = fPar->GetTime0(iplan); - Float_t rowPadSize = geo->GetRowPadSize(iplan,icham,isect); - Float_t colPadSize = geo->GetColPadSize(iplan); - Float_t timeBinSize = geo->GetTimeBinSize(); + Float_t rowPadSize = fPar->GetRowPadSize(iplan,icham,isect); + Float_t colPadSize = fPar->GetColPadSize(iplan); + Float_t timeBinSize = fPar->GetTimeBinSize(); // Loop through all entries in the tree for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) { @@ -320,10 +333,14 @@ Bool_t AliTRDclusterizerV0::MakeClusters() smear[2] = (Int_t) ((time0 - smear[2]) / timeBinSize); // Add the smeared cluster to the output array - Int_t detector = recPoint1->GetDetector(); - Int_t digits[3] = {0}; - Int_t tr[9] = {-1}; - fTRD->AddCluster(smear,digits,detector,0.0,tr,0,0); + Int_t detector = recPoint1->GetDetector(); + Int_t tr[9] = { -1 }; + Float_t pos[3]; + Float_t sigma[2] = { 0.0 }; + pos[0] = smear[1]; + pos[1] = smear[0]; + pos[2] = (time0 - smear[2]) / timeBinSize; + fTRD->AddCluster(pos,detector,0.0,tr,sigma,0); } diff --git a/TRD/AliTRDclusterizerV1.cxx b/TRD/AliTRDclusterizerV1.cxx index 75cfe837054..e62926bdfd3 100644 --- a/TRD/AliTRDclusterizerV1.cxx +++ b/TRD/AliTRDclusterizerV1.cxx @@ -15,6 +15,9 @@ /* $Log$ +Revision 1.16 2002/03/25 20:01:30 cblume +Introduce parameter class + Revision 1.15 2001/11/14 12:09:11 cblume Use correct name for digitizer @@ -107,7 +110,6 @@ AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer() // fDigitsManager = 0; - fPar = 0; } @@ -121,7 +123,6 @@ AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title fDigitsManager = new AliTRDdigitsManager(); fDigitsManager->CreateArrays(); - fPar = 0; } @@ -170,7 +171,6 @@ void AliTRDclusterizerV1::Copy(TObject &c) // ((AliTRDclusterizerV1 &) c).fDigitsManager = 0; - ((AliTRDclusterizerV1 &) c).fPar = 0; AliTRDclusterizer::Copy(c); @@ -217,10 +217,8 @@ Bool_t AliTRDclusterizerV1::MakeClusters() // Create a default parameter class if none is defined if (!fPar) { fPar = new AliTRDparameter("TRDparameter","Standard TRD parameter"); - if (fVerbose > 0) { - printf(" "); - printf("Create the default parameter object.\n"); - } + printf(" "); + printf("Create the default parameter object.\n"); } Float_t timeBinSize = fPar->GetTimeBinSize(); @@ -312,10 +310,15 @@ Bool_t AliTRDclusterizerV1::MakeClusters() ,icham,iplan,isect); } - Int_t nRowMax = fPar->GetRowMax(iplan,icham,isect); - Int_t nColMax = fPar->GetColMax(iplan); - Int_t nTimeBefore = fPar->GetTimeBefore(); - Int_t nTimeTotal = fPar->GetTimeTotal(); + Int_t nRowMax = fPar->GetRowMax(iplan,icham,isect); + Int_t nColMax = fPar->GetColMax(iplan); + Int_t nTimeBefore = fPar->GetTimeBefore(); + Int_t nTimeTotal = fPar->GetTimeTotal(); + + Float_t row0 = fPar->GetRow0(iplan,icham,isect); + Float_t col0 = fPar->GetCol0(iplan); + Float_t rowSize = fPar->GetRowPadSize(iplan,icham,isect); + Float_t colSize = fPar->GetColPadSize(iplan); // Get the digits digits = fDigitsManager->GetDigits(idet); @@ -471,7 +474,8 @@ Bool_t AliTRDclusterizerV1::MakeClusters() // Calculate the position of the cluster by using the // lookup table method - clusterPads[1] = fPar->LUTposition(iplan,clusterSignal[0] + clusterPads[1] = col + 0.5 + + fPar->LUTposition(iplan,clusterSignal[0] ,clusterSignal[1] ,clusterSignal[2]); @@ -486,8 +490,11 @@ Bool_t AliTRDclusterizerV1::MakeClusters() } - Float_t clusterSigmaY2 = (clusterSignal[2] + clusterSignal[0]) / clusterCharge - - (clusterPads[1]-col-0.5) * (clusterPads[1]-col-0.5); + Float_t q0 = clusterSignal[0]; + Float_t q1 = clusterSignal[1]; + Float_t q2 = clusterSignal[2]; + Float_t clusterSigmaY2 = (q1*(q0+q2)+4*q0*q2) / + (clusterCharge*clusterCharge); // Correct for ExB displacement if (fPar->ExBOn()) { @@ -520,13 +527,21 @@ Bool_t AliTRDclusterizerV1::MakeClusters() printf("Type = %d, Number of pads = %d\n",iType,nPadCount); } + // Calculate the position and the error + Float_t clusterPos[3]; + clusterPos[0] = clusterPads[1] * colSize + col0; + clusterPos[1] = clusterPads[0] * rowSize + row0; + clusterPos[2] = clusterPads[2]; + Float_t clusterSig[2]; + clusterSig[0] = (clusterSigmaY2 + 1./12.) * colSize*colSize; + clusterSig[1] = rowSize * rowSize / 12.; + // Add the cluster to the output array - fTRD->AddCluster(clusterPads - ,clusterDigit + fTRD->AddCluster(clusterPos ,idet ,clusterCharge ,clusterTracks - ,clusterSigmaY2 + ,clusterSig ,iType); } diff --git a/TRD/AliTRDclusterizerV1.h b/TRD/AliTRDclusterizerV1.h index 01f168c6dd5..da30c6ea06e 100644 --- a/TRD/AliTRDclusterizerV1.h +++ b/TRD/AliTRDclusterizerV1.h @@ -12,7 +12,6 @@ /////////////////////////////////////////////////////// class AliTRDdigitsManager; -class AliTRDparameter; class AliTRDclusterizerV1 : public AliTRDclusterizer { @@ -27,20 +26,16 @@ class AliTRDclusterizerV1 : public AliTRDclusterizer { virtual void Copy(TObject &c); virtual Bool_t MakeClusters(); virtual Bool_t ReadDigits(); - virtual void SetParameter(AliTRDparameter *par) { fPar = par; }; - - AliTRDparameter *GetParameter() const { return fPar; }; protected: AliTRDdigitsManager *fDigitsManager; //! TRD digits manager - AliTRDparameter *fPar; // TRD digitization parameter object private: virtual Float_t Unfold(Float_t eps, Int_t plane, Float_t *padSignal); - ClassDef(AliTRDclusterizerV1,4) // TRD-Cluster finder, slow simulator + ClassDef(AliTRDclusterizerV1,5) // TRD-Cluster finder, slow simulator }; diff --git a/TRD/AliTRDgeometry.cxx b/TRD/AliTRDgeometry.cxx index b007e431496..e12cc6b66bf 100644 --- a/TRD/AliTRDgeometry.cxx +++ b/TRD/AliTRDgeometry.cxx @@ -15,6 +15,9 @@ /* $Log$ +Revision 1.17 2002/04/05 13:20:12 cblume +Remove const for CreateGeometry + Revision 1.16 2002/03/28 14:59:07 cblume Coding conventions @@ -107,7 +110,7 @@ Add new TRD classes #include "AliMC.h" #include "AliTRDgeometry.h" -#include "AliTRDrecPoint.h" +#include "AliTRDparameter.h" #include "AliMC.h" ClassImp(AliTRDgeometry) @@ -287,28 +290,10 @@ void AliTRDgeometry::Init() } } - // The pad size in column direction (rphi-direction) - SetColPadSize(0,0.65); - SetColPadSize(1,0.68); - SetColPadSize(2,0.71); - SetColPadSize(3,0.74); - SetColPadSize(4,0.77); - SetColPadSize(5,0.80); - - // The pad row (z-direction) - SetNRowPad(); - - // The number of time bins. Default is 100 ns timbin size - SetNTimeBin(15); - - // Additional time bins before and after the drift region. - // Default is to only sample the drift region - SetExpandTimeBin(0,0); - // The rotation matrix elements Float_t phi = 0; for (isect = 0; isect < fgkNsect; isect++) { - phi = -2.0 * kPI / (Float_t) fgkNsect * ((Float_t) isect + 0.5); + phi = -2.0 * TMath::Pi() / (Float_t) fgkNsect * ((Float_t) isect + 0.5); fRotA11[isect] = TMath::Cos(phi); fRotA12[isect] = TMath::Sin(phi); fRotA21[isect] = TMath::Sin(phi); @@ -322,97 +307,6 @@ void AliTRDgeometry::Init() } -//_____________________________________________________________________________ -void AliTRDgeometry::SetNRowPad(const Int_t p, const Int_t c, const Int_t npad) -{ - // - // Redefines the number of pads in raw direction for - // a given plane and chamber number - // - - for (Int_t isect = 0; isect < fgkNsect; isect++) { - - fRowMax[p][c][isect] = npad; - - fRowPadSize[p][c][isect] = (fClength[p][c] - 2.*fgkRpadW) - / ((Float_t) npad); - - } - -} - -//_____________________________________________________________________________ -void AliTRDgeometry::SetNRowPad() -{ - // - // Defines the number of pads in row direction - // - - Int_t isect; - Int_t icham; - Int_t iplan; - - Int_t rowMax[kNplan][kNcham] = { { 16, 16, 12, 16, 16 } - , { 16, 16, 12, 16, 16 } - , { 16, 16, 12, 16, 16 } - , { 16, 16, 12, 16, 16 } - , { 14, 16, 12, 16, 14 } - , { 13, 16, 12, 16, 13 } }; - - for (isect = 0; isect < kNsect; isect++) { - for (icham = 0; icham < kNcham; icham++) { - for (iplan = 0; iplan < kNplan; iplan++) { - - fRowMax[iplan][icham][isect] = rowMax[iplan][icham]; - - fRowPadSize[iplan][icham][isect] = (fClength[iplan][icham] - 2.*fgkRpadW) - / ((Float_t) rowMax[iplan][icham]); - - Float_t row0 = fgkRpadW - fClength[iplan][0] - - fClength[iplan][1] - - fClength[iplan][2]/2.; - for (Int_t ic = 0; ic < icham; ic++) { - row0 += fClength[iplan][ic]; - } - fRow0[iplan][icham][isect] = row0; - - } - } - } - -} - -//_____________________________________________________________________________ -void AliTRDgeometry::SetColPadSize(const Int_t p, const Float_t s) -{ - // - // Redefines the pad size in column direction - // - - fColPadSize[p] = s; - fCol0[p] = - fCwidth[p]/2. + fgkCpadW; - fColMax[p] = ((Int_t) ((fCwidth[p] - 2. * fgkCpadW) / s)); - -} - -//_____________________________________________________________________________ -void AliTRDgeometry::SetNTimeBin(const Int_t nbin) -{ - // - // Redefines the number of time bins in the drift region. - // The time bin width is defined by the length of the - // drift region divided by . - // - - fTimeMax = nbin; - fTimeBinSize = fgkDrThick / ((Float_t) fTimeMax); - for (Int_t iplan = 0; iplan < fgkNplan; iplan++) { - fTime0[iplan] = fgkRmin + fgkCraH + fgkCdrH - + iplan * (fgkCH + fgkVspace); - } - -} - //_____________________________________________________________________________ void AliTRDgeometry::CreateGeometry(Int_t *idtmed) { @@ -423,7 +317,9 @@ void AliTRDgeometry::CreateGeometry(Int_t *idtmed) } //_____________________________________________________________________________ -Bool_t AliTRDgeometry::Local2Global(Int_t idet, Float_t *local, Float_t *global) const +Bool_t AliTRDgeometry::Local2Global(Int_t idet, Float_t *local + , Float_t *global + , AliTRDparameter *par) const { // // Converts local pad-coordinates (row,col,time) into @@ -434,35 +330,44 @@ Bool_t AliTRDgeometry::Local2Global(Int_t idet, Float_t *local, Float_t *global) Int_t isect = GetSector(idet); // Sector info (0-17) Int_t iplan = GetPlane(idet); // Plane info (0-5) - return Local2Global(iplan,icham,isect,local,global); + return Local2Global(iplan,icham,isect,local,global,par); } //_____________________________________________________________________________ Bool_t AliTRDgeometry::Local2Global(Int_t iplan, Int_t icham, Int_t isect - , Float_t *local, Float_t *global) const + , Float_t *local, Float_t *global + , AliTRDparameter *par) const { // // Converts local pad-coordinates (row,col,time) into // global ALICE reference frame coordinates (x,y,z) // + if (!par) { + Error("Local2Global","No parameter defined\n"); + return kFALSE; + } + Int_t idet = GetDetector(iplan,icham,isect); // Detector number Float_t padRow = local[0]+0.5; // Pad Row position Float_t padCol = local[1]+0.5; // Pad Column position Float_t timeSlice = local[2]+0.5; // Time "position" - Float_t row0 = GetRow0(iplan,icham,isect); - Float_t col0 = GetCol0(iplan); - Float_t time0 = GetTime0(iplan); + Float_t row0 = par->GetRow0(iplan,icham,isect); + Float_t col0 = par->GetCol0(iplan); + Float_t time0 = par->GetTime0(iplan); Float_t rot[3]; // calculate (x,y,z) position in rotated chamber - rot[0] = time0 - (timeSlice - fTimeBefore) * fTimeBinSize; - rot[1] = col0 + padCol * fColPadSize[iplan]; - rot[2] = row0 + padRow * fRowPadSize[iplan][icham][isect]; + rot[0] = time0 - (timeSlice - par->GetTimeBefore()) + * par->GetTimeBinSize(); + rot[1] = col0 + padCol + * par->GetColPadSize(iplan); + rot[2] = row0 + padRow + * par->GetRowPadSize(iplan,icham,isect); // Rotate back to original position return RotateBack(idet,rot,global); @@ -562,43 +467,3 @@ Int_t AliTRDgeometry::GetSector(const Int_t d) const } -//_____________________________________________________________________________ -void AliTRDgeometry::GetGlobal(const AliRecPoint *p, TVector3 &pos - , TMatrix &mat) const -{ - // - // Returns the global coordinate and error matrix of a AliTRDrecPoint - // - - GetGlobal(p,pos); - mat.Zero(); - -} - -//_____________________________________________________________________________ -void AliTRDgeometry::GetGlobal(const AliRecPoint *p, TVector3 &pos) const -{ - // - // Returns the global coordinate and error matrix of a AliTRDrecPoint - // - - Int_t detector = ((AliTRDrecPoint *) p)->GetDetector(); - - Float_t global[3]; - Float_t local[3]; - local[0] = ((AliTRDrecPoint *) p)->GetLocalRow(); - local[1] = ((AliTRDrecPoint *) p)->GetLocalCol(); - local[2] = ((AliTRDrecPoint *) p)->GetLocalTime(); - - if (Local2Global(detector,local,global)) { - pos.SetX(global[0]); - pos.SetY(global[1]); - pos.SetZ(global[2]); - } - else { - pos.SetX(0.0); - pos.SetY(0.0); - pos.SetZ(0.0); - } - -} diff --git a/TRD/AliTRDgeometry.h b/TRD/AliTRDgeometry.h index 541beeefdc1..0f1b754aff2 100644 --- a/TRD/AliTRDgeometry.h +++ b/TRD/AliTRDgeometry.h @@ -13,6 +13,8 @@ #include "AliGeometry.h" +class AliTRDparameter; + class AliTRDgeometry : public AliGeometry { public: @@ -25,8 +27,9 @@ class AliTRDgeometry : public AliGeometry { virtual void CreateGeometry(Int_t *idtmed); virtual Int_t IsVersion() const = 0; virtual void Init(); - virtual Bool_t Local2Global(Int_t d, Float_t *local, Float_t *global) const; - virtual Bool_t Local2Global(Int_t p, Int_t c, Int_t s, Float_t *local, Float_t *global) const; + virtual Bool_t Impact(const TParticle * particle) const { return kTRUE; }; + virtual Bool_t Local2Global(Int_t d, Float_t *local, Float_t *global, AliTRDparameter *par) const; + virtual Bool_t Local2Global(Int_t p, Int_t c, Int_t s, Float_t *local, Float_t *global, AliTRDparameter *par) const; virtual Bool_t Rotate(Int_t d, Float_t *pos, Float_t *rot) const; virtual Bool_t RotateBack(Int_t d, Float_t *rot, Float_t *pos) const; @@ -58,14 +61,6 @@ class AliTRDgeometry : public AliGeometry { virtual void SetPHOShole() = 0; virtual void SetRICHhole() = 0; - virtual void SetNRowPad(); - virtual void SetNRowPad(const Int_t p, const Int_t c, const Int_t npad); - virtual void SetColPadSize(const Int_t p, const Float_t s); - virtual void SetNTimeBin(const Int_t nbin); - virtual void SetExpandTimeBin(const Int_t nbefore, const Int_t nafter) - { fTimeBefore = nbefore; - fTimeAfter = nafter; }; - virtual Bool_t GetPHOShole() const = 0; virtual Bool_t GetRICHhole() const = 0; @@ -78,30 +73,9 @@ class AliTRDgeometry : public AliGeometry { Float_t GetChamberWidth(const Int_t p) const { return fCwidth[p]; }; Float_t GetChamberLength(const Int_t p, const Int_t c) const { return fClength[p][c]; }; - - Int_t GetRowMax(const Int_t p, const Int_t c, const Int_t s) - const { return fRowMax[p][c][s]; }; - Int_t GetColMax(const Int_t p) const { return fColMax[p]; }; - Int_t GetTimeMax() const { return fTimeMax; }; - Int_t GetTimeBefore() const { return fTimeBefore; }; - Int_t GetTimeAfter() const { return fTimeAfter; }; - Int_t GetTimeTotal() const { return fTimeMax - + fTimeBefore - + fTimeAfter; }; - - Float_t GetRow0(const Int_t p, const Int_t c, const Int_t s) - const { return fRow0[p][c][s]; }; - Float_t GetCol0(const Int_t p) const { return fCol0[p]; }; - Float_t GetTime0(const Int_t p) const { return fTime0[p]; }; - - Float_t GetRowPadSize(const Int_t p, const Int_t c, const Int_t s) - const { return fRowPadSize[p][c][s]; }; - Float_t GetColPadSize(const Int_t p) const { return fColPadSize[p]; }; - Float_t GetTimeBinSize() const { return fTimeBinSize; }; - - virtual void GetGlobal(const AliRecPoint *p, TVector3 &pos, TMatrix &mat) const; - virtual void GetGlobal(const AliRecPoint *p, TVector3 &pos) const; - virtual Bool_t Impact(const TParticle * particle) const {return kTRUE;} + virtual void GetGlobal(const AliRecPoint *p, TVector3 &pos, TMatrix &mat) const { }; + virtual void GetGlobal(const AliRecPoint *p, TVector3 &pos) const { }; + static Double_t GetAlpha() { return 2 * 3.14159265358979323846 / fgkNsect; }; protected: @@ -165,25 +139,11 @@ class AliTRDgeometry : public AliGeometry { static const Float_t fgkCoZpos; // Position of the PE of the cooling device static const Float_t fgkWaZpos; // Position of the colling water - Int_t fRowMax[kNplan][kNcham][kNsect]; // Number of pad-rows - Int_t fColMax[kNplan]; // Number of pad-columns - Int_t fTimeMax; // Number of timebins in the drift region - Int_t fTimeBefore; // Number of timebins before the drift region - Int_t fTimeAfter; // Number of timebins after the drift region - Float_t fCwidth[kNplan]; // Outer widths of the chambers Float_t fClength[kNplan][kNcham]; // Outer lengths of the chambers Float_t fClengthPH[kNplan][kNcham]; // For sectors with holes for the PHOS Float_t fClengthRH[kNplan][kNcham]; // For sectors with holes for the RICH - Float_t fRow0[kNplan][kNcham][kNsect]; // Row-position of pad 0 - Float_t fCol0[kNplan]; // Column-position of pad 0 - Float_t fTime0[kNplan]; // Time-position of pad 0 - - Float_t fRowPadSize[kNplan][kNcham][kNsect]; // Pad size in z-direction - Float_t fColPadSize[kNplan]; // Pad size in rphi-direction - Float_t fTimeBinSize; // Size of the time buckets - Float_t fRotA11[kNsect]; // Matrix elements for the rotation Float_t fRotA12[kNsect]; // Matrix elements for the rotation Float_t fRotA21[kNsect]; // Matrix elements for the rotation @@ -194,7 +154,7 @@ class AliTRDgeometry : public AliGeometry { Float_t fRotB21[kNsect]; // Matrix elements for the backward rotation Float_t fRotB22[kNsect]; // Matrix elements for the backward rotation - ClassDef(AliTRDgeometry,4) // TRD geometry base class + ClassDef(AliTRDgeometry,5) // TRD geometry base class }; diff --git a/TRD/AliTRDgeometryDetail.cxx b/TRD/AliTRDgeometryDetail.cxx index c8c2597f21f..c20b7d8ace6 100644 --- a/TRD/AliTRDgeometryDetail.cxx +++ b/TRD/AliTRDgeometryDetail.cxx @@ -15,6 +15,9 @@ /* $Log$ +Revision 1.4 2002/03/28 14:59:07 cblume +Coding conventions + Revision 1.3 2002/02/11 14:21:16 cblume Update of the geometry. Get rid of MANY @@ -35,6 +38,7 @@ Add detailed geometry and simple simulator #include "AliMC.h" #include "AliTRDgeometryDetail.h" +#include "AliTRDparameter.h" ClassImp(AliTRDgeometryDetail) @@ -428,15 +432,17 @@ void AliTRDgeometryDetail::PositionReadout(Int_t ipla, Int_t icha) const Int_t kNmcmChannel = 18; - Int_t nMCMrow = GetRowMax(ipla,icha,0); - Int_t nMCMcol = GetColMax(ipla) / kNmcmChannel; + AliTRDparameter *parameter = new AliTRDparameter(); + + Int_t nMCMrow = parameter->GetRowMax(ipla,icha,0); + Int_t nMCMcol = parameter->GetColMax(ipla) / kNmcmChannel; Float_t xSize = (GetChamberWidth(ipla) - 2.*fgkCpadW) / ((Float_t) nMCMcol); Float_t ySize = (GetChamberLength(ipla,icha) - 2.*fgkRpadW) / ((Float_t) nMCMrow); - Float_t x0 = GetCol0(ipla); - Float_t y0 = GetRow0(ipla,icha,0); + Float_t x0 = parameter->GetCol0(ipla); + Float_t y0 = parameter->GetRow0(ipla,icha,0); Int_t iCopy = GetDetector(ipla,icha,0) * 1000; for (Int_t iMCMrow = 0; iMCMrow < nMCMrow; iMCMrow++) { @@ -450,6 +456,8 @@ void AliTRDgeometryDetail::PositionReadout(Int_t ipla, Int_t icha) } } + delete parameter; + } //_____________________________________________________________________________ @@ -500,12 +508,14 @@ void AliTRDgeometryDetail::PositionCooling(Int_t ipla, Int_t icha, Int_t idrotm) Float_t ypos; Float_t zpos; + AliTRDparameter *parameter = new AliTRDparameter(); + Int_t iCopy = GetDetector(ipla,icha,0) * 100; - Int_t nMCMrow = GetRowMax(ipla,icha,0); + Int_t nMCMrow = parameter->GetRowMax(ipla,icha,0); Float_t ySize = (GetChamberLength(ipla,icha) - 2.*fgkRpadW) / ((Float_t) nMCMrow); - Float_t y0 = GetRow0(ipla,icha,0); + Float_t y0 = parameter->GetRow0(ipla,icha,0); // Position the cooling pipes for (Int_t iMCMrow = 0; iMCMrow < nMCMrow; iMCMrow++) { @@ -522,4 +532,6 @@ void AliTRDgeometryDetail::PositionCooling(Int_t ipla, Int_t icha, Int_t idrotm) } + delete parameter; + } diff --git a/TRD/AliTRDparameter.cxx b/TRD/AliTRDparameter.cxx index f2d2e46cbd7..0ed581725a6 100644 --- a/TRD/AliTRDparameter.cxx +++ b/TRD/AliTRDparameter.cxx @@ -15,6 +15,9 @@ /* $Log$ +Revision 1.5 2002/04/30 08:30:40 cblume +gAlice now only read by AliRunDigitizer. Therefore it is just deleted in AliTRDmerge.C + Revision 1.4 2002/04/12 12:13:23 cblume Add Jiris changes @@ -382,7 +385,8 @@ void AliTRDparameter::Init() fTimeCoupling = 0.4; // The tilting angle for the readout pads - SetTiltingAngle(5.0); + //SetTiltingAngle(5.0); + SetTiltingAngle(0.0); // The magnetic field strength in Tesla //fField = 0.2 * gAlice->Field()->Factor(); diff --git a/TRD/AliTRDpid.cxx b/TRD/AliTRDpid.cxx index 6f06f7138ac..edad9fa9008 100644 --- a/TRD/AliTRDpid.cxx +++ b/TRD/AliTRDpid.cxx @@ -15,6 +15,9 @@ /* $Log$ +Revision 1.4 2001/11/19 08:44:08 cblume +Fix bugs reported by Rene + Revision 1.3 2001/11/14 10:50:46 cblume Changes in digits IO. Add merging of summable digits @@ -361,7 +364,7 @@ Bool_t AliTRDpid::ReadCluster(const Char_t *name) printf("AliTRDpid::ReadCluster -- "); printf("Open file %s\n",name); - AliTRDtracker *tracker = new AliTRDtracker("dummy","dummy"); + AliTRDtracker *tracker = new AliTRDtracker(); tracker->ReadClusters(fClusterArray,name); if (!fClusterArray) { @@ -458,7 +461,7 @@ Int_t AliTRDpid::MCpid(const AliTRDtrack *t) } // Loop through all clusters associated to this track - Int_t nCluster = t->GetNclusters(); + Int_t nCluster = t->GetNumberOfClusters(); for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { // Get a cluster @@ -468,8 +471,8 @@ Int_t AliTRDpid::MCpid(const AliTRDtrack *t) } // Get the first two MC track indices - Int_t track0 = cluster->GetTrackIndex(0); - Int_t track1 = cluster->GetTrackIndex(1); + Int_t track0 = cluster->GetLabel(0); + Int_t track1 = cluster->GetLabel(1); // Check on the track index to find the right primaries if ((track0 > fPIDindexMin) && @@ -574,7 +577,7 @@ Int_t AliTRDpid::MCpid(const AliTRDtrack *t, Int_t *pdg } // Loop through all clusters associated to this track - Int_t nCluster = t->GetNclusters(); + Int_t nCluster = t->GetNumberOfClusters(); for (iCluster = 0; iCluster < nCluster; iCluster++) { // Get a cluster @@ -586,7 +589,7 @@ Int_t AliTRDpid::MCpid(const AliTRDtrack *t, Int_t *pdg // Get the MC track indices for (iTrack = 0; iTrack < kNtrack; iTrack++) { - Int_t trackIndex = cluster->GetTrackIndex(iTrack); + Int_t trackIndex = cluster->GetLabel(iTrack); if (trackIndex >= 0) { particle = gAlice->Particle(trackIndex); Int_t pdgCode = particle->GetPdgCode(); @@ -655,7 +658,7 @@ Bool_t AliTRDpid::SumCharge(const AliTRDtrack *t } // Loop through all clusters associated to this track - Int_t nClus = t->GetNclusters(); + Int_t nClus = t->GetNumberOfClusters(); for (Int_t iClus = 0; iClus < nClus; iClus++) { // Get a cluster diff --git a/TRD/AliTRDrecPoint.cxx b/TRD/AliTRDrecPoint.cxx index 0eee224fe96..5d153363da6 100644 --- a/TRD/AliTRDrecPoint.cxx +++ b/TRD/AliTRDrecPoint.cxx @@ -15,6 +15,9 @@ /* $Log$ +Revision 1.8 2002/03/28 14:59:07 cblume +Coding conventions + Revision 1.7 2001/12/05 15:04:34 hristov Changes related to the corrections of AliRecPoint @@ -146,25 +149,25 @@ void AliTRDrecPoint::SetLocalPosition(TVector3 &pos) // system. // - const Float_t kSq12 = 3.464101615; + //const Float_t kSq12 = 3.464101615; // Set the position - fLocPos = pos; + //fLocPos = pos; // Set the error matrix // row: pad-size / sqrt(12) // col: not defined yet // time: bin-size / sqrt(12) - Int_t plane = ((AliTRDgeometry *) fGeom)->GetPlane(fDetector); - Int_t chamber = ((AliTRDgeometry *) fGeom)->GetChamber(fDetector); - Int_t sector = ((AliTRDgeometry *) fGeom)->GetSector(fDetector); - fLocPosM->operator()(0,0) = ((AliTRDgeometry *) fGeom)->GetRowPadSize(plane - ,chamber - ,sector) - / kSq12; - fLocPosM->operator()(1,1) = 0.0; - fLocPosM->operator()(2,2) = ((AliTRDgeometry *) fGeom)->GetTimeBinSize() - / kSq12; + //Int_t plane = ((AliTRDgeometry *) fGeom)->GetPlane(fDetector); + //Int_t chamber = ((AliTRDgeometry *) fGeom)->GetChamber(fDetector); + //Int_t sector = ((AliTRDgeometry *) fGeom)->GetSector(fDetector); + //fLocPosM->operator()(0,0) = ((AliTRDgeometry *) fGeom)->GetRowPadSize(plane + // ,chamber + // ,sector) + // / kSq12; + //fLocPosM->operator()(1,1) = 0.0; + //fLocPosM->operator()(2,2) = ((AliTRDgeometry *) fGeom)->GetTimeBinSize() + // / kSq12; // printf("rec. point: row = %f, col = %f, time = %f \n", // fLocPos[0],fLocPos[1],fLocPos[2]); @@ -179,34 +182,34 @@ void AliTRDrecPoint::SetTrackingYZ(Float_t sigmaY, Float_t sigmaZ) // of tracking sector // - Int_t plane = ((AliTRDgeometry *) fGeom)->GetPlane(fDetector); - Int_t chamber = ((AliTRDgeometry *) fGeom)->GetChamber(fDetector); - Int_t sector = ((AliTRDgeometry *) fGeom)->GetSector(fDetector); + //Int_t plane = ((AliTRDgeometry *) fGeom)->GetPlane(fDetector); + //Int_t chamber = ((AliTRDgeometry *) fGeom)->GetChamber(fDetector); + //Int_t sector = ((AliTRDgeometry *) fGeom)->GetSector(fDetector); // Set the position - Float_t padRow = fLocPos[0]; // Pad Row position - Float_t padCol = fLocPos[1]; // Pad Column position + //Float_t padRow = fLocPos[0]; // Pad Row position + //Float_t padCol = fLocPos[1]; // Pad Column position - Float_t col0 = ((AliTRDgeometry *) fGeom)->GetCol0(plane); - Float_t row0 = ((AliTRDgeometry *) fGeom)->GetRow0(plane,chamber,sector); + //Float_t col0 = ((AliTRDgeometry *) fGeom)->GetCol0(plane); + //Float_t row0 = ((AliTRDgeometry *) fGeom)->GetRow0(plane,chamber,sector); // Float_t offset = 0.5 * ((AliTRDgeometry *) fGeom)->GetChamberWidth(plane); - fY = - (col0 + padCol * ((AliTRDgeometry *) fGeom)->GetColPadSize(plane)); - fZ = row0 + padRow * ((AliTRDgeometry *) fGeom)->GetRowPadSize(plane - ,chamber - ,sector); + //fY = - (col0 + padCol * ((AliTRDgeometry *) fGeom)->GetColPadSize(plane)); + //fZ = row0 + padRow * ((AliTRDgeometry *) fGeom)->GetRowPadSize(plane + // ,chamber + // ,sector); // fSigmaY = sigmaY * sigmaY; // fSigmaZ = sigmaZ * sigmaZ; - fSigmaY2 = 0.05 * 0.05; +//fSigmaY2 = 0.05 * 0.05; - fSigmaZ2 = ((AliTRDgeometry *) fGeom)->GetRowPadSize(plane,chamber,sector) - * ((AliTRDgeometry *) fGeom)->GetRowPadSize(plane,chamber,sector) - / 12.; +//fSigmaZ2 = ((AliTRDgeometry *) fGeom)->GetRowPadSize(plane,chamber,sector) +// * ((AliTRDgeometry *) fGeom)->GetRowPadSize(plane,chamber,sector) +// / 12.; } @@ -277,3 +280,10 @@ void AliTRDrecPoint::AddTrackIndex(Int_t *track) return; } + + + + + + + diff --git a/TRD/AliTRDsimpleMC.cxx b/TRD/AliTRDsimpleMC.cxx index 36e01a1c6e0..2805dc16585 100644 --- a/TRD/AliTRDsimpleMC.cxx +++ b/TRD/AliTRDsimpleMC.cxx @@ -14,7 +14,10 @@ **************************************************************************/ /* -$Log$ +$Log$ +Revision 1.1 2001/11/06 17:19:41 cblume +Add detailed geometry and simple simulator + */ /////////////////////////////////////////////////////////////////////////////// @@ -32,6 +35,7 @@ $Log$ #include "AliTRDsimpleMC.h" #include "AliTRDgeometry.h" #include "AliTRDv1.h" +#include "AliTRDparameter.h" ClassImp(AliTRDsimpleMC) @@ -60,6 +64,7 @@ AliTRDsimpleMC::AliTRDsimpleMC():AliMC() fTrackEntering = kFALSE; fTRD = NULL; + fPar = NULL; } @@ -89,6 +94,7 @@ AliTRDsimpleMC::AliTRDsimpleMC(const char *name, const char *title) fTrackEntering = kFALSE; fTRD = NULL; + fPar = NULL; } @@ -158,10 +164,13 @@ void AliTRDsimpleMC::NewTrack(Int_t iTrack, Int_t pdg // Starts a new track. // + if (!fPar) { + fPar = new AliTRDparameter("TRDparameter","Standard TRD parameter"); + } + if (!fTRD) { fTRD = (AliTRDv1 *) gAlice->GetDetector("TRD"); - AliTRDgeometry *geometry = fTRD->GetGeometry(); - fX0 = geometry->GetTime0(0) - AliTRDgeometry::DrThick(); + fX0 = fPar->GetTime0(0) - AliTRDgeometry::DrThick(); } fTRD->ResetHits(); diff --git a/TRD/AliTRDsimpleMC.h b/TRD/AliTRDsimpleMC.h index 9ed9163aec6..ba6a2352ed4 100644 --- a/TRD/AliTRDsimpleMC.h +++ b/TRD/AliTRDsimpleMC.h @@ -15,6 +15,7 @@ #include "AliMCProcess.h" class AliTRDv1; +class AliTRDparameter; class AliTRDsimpleMC : public AliMC { @@ -216,27 +217,28 @@ class AliTRDsimpleMC : public AliMC { , kVolDrCh }; - Float_t fMaxStep; // Maximum step size - Int_t fNStep; // Number of steps - Int_t fTrack; // Track number - Double_t fTrackPx; // Track px - Double_t fTrackPy; // Track py - Double_t fTrackPz; // Track pz - Double_t fTrackPtot; // Track total momentum - Double_t fTrackEtot; // Track total energy - Double_t fTrackX; // Track x position - Double_t fTrackY; // Track y position - Double_t fTrackZ; // Track z position - Double_t fX0; // X position of the beginning of the chamber - Double_t fTrackStep; // Track step size - Int_t fTrackPid; // Track PID - Float_t fTrackCharge; // Track charge - Float_t fTrackMass; // Track particle mass - Bool_t fTrackEntering; // Track entering chamber - - AliTRDv1 *fTRD; //! TRD detector object - - ClassDef(AliTRDsimpleMC,1) // Simple TRD Monte Carlo class + Float_t fMaxStep; // Maximum step size + Int_t fNStep; // Number of steps + Int_t fTrack; // Track number + Double_t fTrackPx; // Track px + Double_t fTrackPy; // Track py + Double_t fTrackPz; // Track pz + Double_t fTrackPtot; // Track total momentum + Double_t fTrackEtot; // Track total energy + Double_t fTrackX; // Track x position + Double_t fTrackY; // Track y position + Double_t fTrackZ; // Track z position + Double_t fX0; // X position of the beginning of the chamber + Double_t fTrackStep; // Track step size + Int_t fTrackPid; // Track PID + Float_t fTrackCharge; // Track charge + Float_t fTrackMass; // Track particle mass + Bool_t fTrackEntering; // Track entering chamber + + AliTRDv1 *fTRD; //! TRD detector object + AliTRDparameter *fPar; //! TRD parameter object + + ClassDef(AliTRDsimpleMC,2) // Simple TRD Monte Carlo class }; #endif diff --git a/TRD/AliTRDtrack.cxx b/TRD/AliTRDtrack.cxx index ca389636cc3..b3d050e286c 100644 --- a/TRD/AliTRDtrack.cxx +++ b/TRD/AliTRDtrack.cxx @@ -35,20 +35,13 @@ Add the tracking code */ -///////////////////////////////////////////////////////////////////////////// -// // -// TRD track object // -// // -///////////////////////////////////////////////////////////////////////////// - #include +#include -#include - -#include "AliTRD.h" #include "AliTRDgeometry.h" #include "AliTRDcluster.h" #include "AliTRDtrack.h" +#include "../TPC/AliTPCtrack.h" ClassImp(AliTRDtrack) @@ -56,15 +49,18 @@ ClassImp(AliTRDtrack) //_____________________________________________________________________________ AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, UInt_t index, -const Double_t xx[5], const Double_t cc[15], Double_t xref, Double_t alpha) { + const Double_t xx[5], const Double_t cc[15], + Double_t xref, Double_t alpha) : AliKalmanTrack() { //----------------------------------------------------------------- // This is the main track constructor. //----------------------------------------------------------------- - fLab=-1; - fChi2=0.; - fdEdx=0.; + + fSeedLab = -1; fAlpha=alpha; + if (fAlpha<-TMath::Pi()) fAlpha += 2*TMath::Pi(); + if (fAlpha>=TMath::Pi()) fAlpha -= 2*TMath::Pi(); + fX=xref; fY=xx[0]; fZ=xx[1]; fC=xx[2]; fE=xx[3]; fT=xx[4]; @@ -75,24 +71,28 @@ const Double_t xx[5], const Double_t cc[15], Double_t xref, Double_t alpha) { fCey=cc[6]; fCez=cc[7]; fCec=cc[8]; fCee=cc[9]; fCty=cc[10]; fCtz=cc[11]; fCtc=cc[12]; fCte=cc[13]; fCtt=cc[14]; - fN=0; - fIndex[fN]=index; + fIndex[0]=index; + SetNumberOfClusters(1); + + fdEdx=0.; - Float_t q = c->GetQ(); + Double_t q = TMath::Abs(c->GetQ()); Double_t s = fX*fC - fE, t=fT; - q *= TMath::Sqrt((1-s*s)/(1+t*t)); - fdQdl[fN++] = q; + if(s*s < 1) q *= TMath::Sqrt((1-s*s)/(1+t*t)); + + fdQdl[0] = q; } //_____________________________________________________________________________ -AliTRDtrack::AliTRDtrack(const AliTRDtrack& t) { +AliTRDtrack::AliTRDtrack(const AliTRDtrack& t) : AliKalmanTrack(t) { // // Copy constructor. // - fLab=t.fLab; + SetLabel(t.GetLabel()); + fSeedLab=t.GetSeedLabel(); - fChi2=t.fChi2; + SetChi2(t.GetChi2()); fdEdx=t.fdEdx; fAlpha=t.fAlpha; @@ -106,12 +106,87 @@ AliTRDtrack::AliTRDtrack(const AliTRDtrack& t) { fCey=t.fCey; fCez=t.fCez; fCec=t.fCec; fCee=t.fCee; fCty=t.fCty; fCtz=t.fCtz; fCtc=t.fCtc; fCte=t.fCte; fCtt=t.fCtt; - fN=t.fN; - for (Int_t i=0; i= TMath::Pi()) fAlpha -= 2*TMath::Pi(); + + Double_t x, p[5]; t.GetExternalParameters(x,p); + + fX=x; + + x = GetConvConst(); + + fY=p[0]; fZ=p[1]; fC=p[4]/x; + fE=fX*fC-p[2]; fT=p[3]; + + //Conversion of the covariance matrix + Double_t c[15]; t.GetExternalCovariance(c); + + c[10]/=x; c[11]/=x; c[12]/=x; c[13]/=x; c[14]/=x*x; + + Double_t c22=fX*fX*c[14] - 2*fX*c[12] + c[5]; + Double_t c32=fX*c[13] - c[8]; + Double_t c20=fX*c[10] - c[3], c21=fX*c[11] - c[4], c42=fX*c[14] - c[12]; + + fCyy=c[0 ]; + fCzy=c[1 ]; fCzz=c[2 ]; + fCcy=c[10]; fCcz=c[11]; fCcc=c[14]; + fCey=c20; fCez=c21; fCec=c42; fCee=c22; + fCty=c[6 ]; fCtz=c[7 ]; fCtc=c[13]; fCte=c32; fCtt=c[9 ]; + +} + +//____________________________________________________________________________ +void AliTRDtrack::GetExternalParameters(Double_t& xr, Double_t x[5]) const { + // + // This function returns external TRD track representation + // + xr=fX; + x[0]=GetY(); + x[1]=GetZ(); + x[2]=GetSnp(); + x[3]=GetTgl(); + x[4]=fC*GetConvConst(); +} + +//_____________________________________________________________________________ +void AliTRDtrack::GetExternalCovariance(Double_t cc[15]) const { + // + // This function returns external representation of the covriance matrix. + // + Double_t a=GetConvConst(); + + Double_t c22=fX*fX*fCcc-2*fX*fCec+fCee; + Double_t c32=fX*fCtc-fCte; + Double_t c20=fX*fCcy-fCey, c21=fX*fCcz-fCez, c42=fX*fCcc-fCec; + + cc[0 ]=fCyy; + cc[1 ]=fCzy; cc[2 ]=fCzz; + cc[3 ]=c20; cc[4 ]=c21; cc[5 ]=c22; + cc[6 ]=fCty; cc[7 ]=fCtz; cc[8 ]=c32; cc[9 ]=fCtt; + cc[10]=fCcy*a; cc[11]=fCcz*a; cc[12]=c42*a; cc[13]=fCtc*a; cc[14]=fCcc*a*a; +} + //_____________________________________________________________________________ void AliTRDtrack::GetCovariance(Double_t cc[15]) const { @@ -125,7 +200,7 @@ void AliTRDtrack::GetCovariance(Double_t cc[15]) const { //_____________________________________________________________________________ Int_t AliTRDtrack::Compare(const TObject *o) const { -// Compares tracks according to their Y2 +// Compares tracks according to their Y2 or curvature AliTRDtrack *t=(AliTRDtrack*)o; // Double_t co=t->GetSigmaY2(); @@ -144,13 +219,17 @@ void AliTRDtrack::CookdEdx(Double_t low, Double_t up) { //----------------------------------------------------------------- // Calculates dE/dX within the "low" and "up" cuts. //----------------------------------------------------------------- + Int_t i; - Int_t nc=GetNclusters(); + Int_t nc=GetNumberOfClusters(); - Float_t sorted[200]; - for (i=0; i<200; i++) sorted[i]=fdQdl[i]; + Float_t sorted[kMAX_CLUSTERS_PER_TRACK]; + for (i=0; i < nc; i++) { + sorted[i]=fdQdl[i]; + } Int_t swap; + do { swap=0; for (i=0; i= 0.99999) { - if (fN>4) cerr<4) cerr< 1) return 0; + Double_t r1=sqrt(1.- c1*c1); + Double_t c2=fC*x2 - fE; + if((c2*c2) > 1) return 0; + Double_t r2=sqrt(1.- c2*c2); fY += dx*(c1+c2)/(r1+r2); fZ += dx*(c1+c2)/(c1*r2 + c2*r1)*fT; @@ -225,6 +311,7 @@ Int_t AliTRDtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm) Double_t d=sqrt((x1-fX)*(x1-fX)+(y1-fY)*(y1-fY)+(z1-fZ)*(z1-fZ)); Double_t p2=GetPt()*GetPt()*(1.+fT*fT); + p2 = TMath::Min(p2,1e+08); // to avoid division by (1-1) for stiff tracks Double_t beta2=p2/(p2 + pm*pm); Double_t ey=fC*fX - fE, ez=fT; @@ -240,19 +327,20 @@ Int_t AliTRDtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm) //Energy losses************************ - + if (x1 < x2) d=-d; Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d*rho; - if (x1 < x2) dE=-dE; + SetLength(GetLength()+d); + + cc = fC; fC*=(1.- sqrt(p2+pm*pm)/p2*dE); - //fE*=(1.- sqrt(p2+pm*pm)/p2*dE); + fE+=fX*(fC-cc); return 1; } - //_____________________________________________________________________________ -void AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, UInt_t index) +Int_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, UInt_t index) { // Assignes found cluster to the track and updates track information @@ -270,8 +358,9 @@ void AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, UInt_t index) Double_t dy=c->GetY() - fY, dz=c->GetZ() - fZ; Double_t cur=fC + k20*dy + k21*dz, eta=fE + k30*dy + k31*dz; if (TMath::Abs(cur*fX-eta) >= 0.99999) { - if (fN>4) cerr<4) cerr<GetQ(); - Double_t s = fX*fC - fE, t=fT; - q *= TMath::Sqrt((1-s*s)/(1+t*t)); - fdQdl[fN++] = q; - - fChi2 += chisq; + Int_t n=GetNumberOfClusters(); + fIndex[n]=index; + SetNumberOfClusters(n+1); + SetChi2(GetChi2()+chisq); // cerr<<"in update: fIndex["<= 0.99999) { + Int_t n=GetNumberOfClusters(); + if (n>4) cerr<= 0.99999) { - if (fN>4) cerr<4) cerr<= 0.) { - if (fN>4) cerr<4) cerr< 1) { py = pt; px = 0; } + else if(r < -1) { py = -pt; px = 0; } + else { + y0=fY + sqrt(1.- r*r)/fC; + px=-pt*(fY-y0)*fC; //cos(phi); + py=-pt*(fE-fX*fC); //sin(phi); + } pz=pt*fT; Double_t tmp=px*TMath::Cos(fAlpha) - py*TMath::Sin(fAlpha); py=px*TMath::Sin(fAlpha) + py*TMath::Cos(fAlpha); @@ -416,74 +517,28 @@ void AliTRDtrack::GetPxPyPz(Double_t& px, Double_t& py, Double_t& pz) const } -//____________________________________________________________________________ -void AliTRDtrack::Streamer(TBuffer &R__b) +//_________________________________________________________________________ +void AliTRDtrack::GetGlobalXYZ(Double_t& x, Double_t& y, Double_t& z) const { - Int_t i; + // Returns reconstructed track coordinates in the global system. + + x = fX; y = fY; z = fZ; + Double_t tmp=x*TMath::Cos(fAlpha) - y*TMath::Sin(fAlpha); + y=x*TMath::Sin(fAlpha) + y*TMath::Cos(fAlpha); + x=tmp; + +} - if (R__b.IsReading()) { - Version_t R__v = R__b.ReadVersion(); if (R__v) { } - TObject::Streamer(R__b); - R__b >> fLab; - R__b >> fChi2; - R__b >> fdEdx; - R__b >> fAlpha; - R__b >> fX; - R__b >> fY; - R__b >> fZ; - R__b >> fC; - R__b >> fE; - R__b >> fT; - R__b >> fCyy; - R__b >> fCzy; - R__b >> fCzz; - R__b >> fCcy; - R__b >> fCcz; - R__b >> fCcc; - R__b >> fCey; - R__b >> fCez; - R__b >> fCec; - R__b >> fCee; - R__b >> fCty; - R__b >> fCtz; - R__b >> fCtc; - R__b >> fCte; - R__b >> fCtt; - R__b >> fN; - for (i=0; i> fIndex[i]; - for (i=0; i> fdQdl[i]; - } else { - R__b.WriteVersion(AliTRDtrack::IsA()); - TObject::Streamer(R__b); - R__b << fLab; - R__b << fChi2; - R__b << fdEdx; - R__b << fAlpha; - R__b << fX; - R__b << fY; - R__b << fZ; - R__b << fC; - R__b << fE; - R__b << fT; - R__b << fCyy; - R__b << fCzy; - R__b << fCzz; - R__b << fCcy; - R__b << fCcz; - R__b << fCcc; - R__b << fCey; - R__b << fCez; - R__b << fCec; - R__b << fCee; - R__b << fCty; - R__b << fCtz; - R__b << fCtc; - R__b << fCte; - R__b << fCtt; - R__b << fN; - for (i=0; i +#include +#include class AliTRDcluster; +class AliTPCtrack; + +const unsigned kMAX_CLUSTERS_PER_TRACK=210; -class AliTRDtrack : public TObject { +class AliTRDtrack : public AliKalmanTrack { // Represents reconstructed TRD track public: - AliTRDtrack() { fN=0;} + AliTRDtrack():AliKalmanTrack(){} AliTRDtrack(const AliTRDcluster *c, UInt_t index, const Double_t xx[5], const Double_t cc[15], Double_t xr, Double_t alpha); AliTRDtrack(const AliTRDtrack& t); + AliTRDtrack(const AliKalmanTrack& t, Double_t alpha); Int_t Compare(const TObject *o) const; void CookdEdx(Double_t low=0.05, Double_t up=0.70); @@ -32,48 +31,58 @@ public: Double_t GetC() const {return fC;} Int_t GetClusterIndex(Int_t i) const {return fIndex[i];} Float_t GetClusterdQdl(Int_t i) const {return fdQdl[i];} - Double_t GetChi2() const {return fChi2;} - Double_t GetZchi2() const {return fZchi2;} + void GetCovariance(Double_t cov[15]) const; - Double_t GetdEdx() const {return fdEdx;} + Float_t GetdEdx() const {return fdEdx;} Double_t GetEta() const {return fE;} - Int_t GetLabel() const {return fLab;} - Int_t GetNclusters() const {return fN;} + + void GetExternalCovariance(Double_t cov[15]) const ; + void GetExternalParameters(Double_t& xr, Double_t x[5]) const ; + + Double_t GetLikelihoodElectron() const { return fLhElectron; }; + + Double_t Get1Pt() const {return TMath::Abs(fC*GetConvConst());} Double_t GetP() const { return TMath::Abs(GetPt())*sqrt(1.+GetTgl()*GetTgl()); } - Double_t GetPredictedChi2(const AliTRDcluster *c) const ; - Double_t GetPt() const {return 0.299792458*0.4/GetC()/100;} + Double_t GetPredictedChi2(const AliTRDcluster*) const ; + Double_t GetPt() const {return 1./Get1Pt();} void GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const ; + void GetGlobalXYZ(Double_t &x, Double_t &y, Double_t &z) const ; + Int_t GetSeedLabel() const { return fSeedLab; } Double_t GetSigmaC2() const {return fCcc;} Double_t GetSigmaTgl2() const {return fCtt;} Double_t GetSigmaY2() const {return fCyy;} Double_t GetSigmaZ2() const {return fCzz;} - Double_t GetTgl() const {return fT;} - Double_t GetX() const {return fX;} - Double_t GetY() const {return fY;} // returns running Y - Double_t GetZ() const {return fZ;} - - Float_t GetLikelihoodElectron() const { return fLhElectron; }; - - Bool_t IsSortable() const {return kTRUE;} + Double_t GetSnp() const {return fX*fC - fE;} + Double_t GetTgl() const {return fT;} + Double_t GetX() const {return fX;} + Double_t GetY() const {return fY;} + Double_t GetZ() const {return fZ;} Int_t PropagateTo(Double_t xr, Double_t x0=8.72,Double_t rho=5.86e-3,Double_t pm=0.139); + void ResetCovariance(); Int_t Rotate(Double_t angle); - void SetLabel(Int_t lab=0) {fLab=lab;} void SetdEdx(Float_t dedx) {fdEdx=dedx;} + void SetLikelihoodElectron(Float_t l) { fLhElectron = l; }; + + void SetSampledEdx(Float_t q, Int_t i) { + Double_t s=GetSnp(), t=GetTgl(); + q*= TMath::Sqrt((1-s*s)/(1+t*t)); + fdQdl[i]=q; + } + + void SetSeedLabel(Int_t lab) { fSeedLab=lab; } + + Int_t Update(const AliTRDcluster* c, Double_t chi2, UInt_t i); - void Update(const AliTRDcluster* c, Double_t chi2, UInt_t i); - void SetLikelihoodElectron(Float_t l) { fLhElectron = l; }; protected: - Int_t fLab; // track label - Double_t fChi2; // total chi2 value for R*phi measurements - Double_t fZchi2; // total chi2 value for Z measurements + Int_t fSeedLab; // track label taken from seeding Float_t fdEdx; // dE/dx Double_t fAlpha; // rotation angle @@ -91,9 +100,9 @@ protected: Double_t fCey, fCez, fCec, fCee; // track Double_t fCty, fCtz, fCtc, fCte, fCtt; // parameters - Short_t fN; // number of clusters associated with the track - UInt_t fIndex[200]; // global indexes of these clusters - Float_t fdQdl[200]; // cluster amplitudes corrected for track angles + UInt_t fIndex[kMAX_CLUSTERS_PER_TRACK]; // global indexes of clusters + Float_t fdQdl[kMAX_CLUSTERS_PER_TRACK]; // cluster amplitudes corrected + // for track angles Float_t fLhElectron; // Likelihood to be an electron diff --git a/TRD/AliTRDtracker.cxx b/TRD/AliTRDtracker.cxx index 16b15f15264..65f605ce4cb 100644 --- a/TRD/AliTRDtracker.cxx +++ b/TRD/AliTRDtracker.cxx @@ -53,189 +53,156 @@ Add the tracking code */ -//////////////////////////////////////////////////////////////////////////////// -// // -// The TRD tracker // -// // -//////////////////////////////////////////////////////////////////////////////// - #include #include -#include #include -#include +#include +#include -#include "AliRun.h" -#include "AliTRD.h" #include "AliTRDgeometry.h" -#include "AliTRDrecPoint.h" +#include "AliTRDparameter.h" +#include "AliTRDgeometryDetail.h" #include "AliTRDcluster.h" #include "AliTRDtrack.h" -#include "AliTRDtrackingSector.h" -#include "AliTRDtimeBin.h" +#include "../TPC/AliTPCtrack.h" #include "AliTRDtracker.h" ClassImp(AliTRDtracker) - const Float_t AliTRDtracker::fgkSeedDepth = 0.5; - const Float_t AliTRDtracker::fgkSeedStep = 0.05; - const Float_t AliTRDtracker::fgkSeedGap = 0.25; + const Float_t AliTRDtracker::fSeedDepth = 0.5; + const Float_t AliTRDtracker::fSeedStep = 0.10; + const Float_t AliTRDtracker::fSeedGap = 0.25; + + const Float_t AliTRDtracker::fMaxSeedDeltaZ12 = 40.; + const Float_t AliTRDtracker::fMaxSeedDeltaZ = 25.; + const Float_t AliTRDtracker::fMaxSeedC = 0.0052; + const Float_t AliTRDtracker::fMaxSeedTan = 1.2; + const Float_t AliTRDtracker::fMaxSeedVertexZ = 150.; - const Float_t AliTRDtracker::fgkMaxSeedDeltaZ12 = 40.; - const Float_t AliTRDtracker::fgkMaxSeedDeltaZ = 25.; - const Float_t AliTRDtracker::fgkMaxSeedC = 0.0052; - const Float_t AliTRDtracker::fgkMaxSeedTan = 1.2; - const Float_t AliTRDtracker::fgkMaxSeedVertexZ = 150.; + const Double_t AliTRDtracker::fSeedErrorSY = 0.2; + const Double_t AliTRDtracker::fSeedErrorSY3 = 2.5; + const Double_t AliTRDtracker::fSeedErrorSZ = 0.1; - const Double_t AliTRDtracker::fgkSeedErrorSY = 0.2; - const Double_t AliTRDtracker::fgkSeedErrorSY3 = 2.5; - const Double_t AliTRDtracker::fgkSeedErrorSZ = 0.1; + const Float_t AliTRDtracker::fMinClustersInSeed = 0.7; - const Float_t AliTRDtracker::fgkMinClustersInSeed = 0.7; + const Float_t AliTRDtracker::fMinClustersInTrack = 0.5; + const Float_t AliTRDtracker::fMinFractionOfFoundClusters = 0.8; - const Float_t AliTRDtracker::fgkMinClustersInTrack = 0.5; - const Float_t AliTRDtracker::fgkSkipDepth = 0.05; - const Float_t AliTRDtracker::fgkLabelFraction = 0.5; - const Float_t AliTRDtracker::fgkWideRoad = 20.; + const Float_t AliTRDtracker::fSkipDepth = 0.05; + const Float_t AliTRDtracker::fLabelFraction = 0.8; + const Float_t AliTRDtracker::fWideRoad = 20.; + + const Double_t AliTRDtracker::fMaxChi2 = 24.; - const Double_t AliTRDtracker::fgkMaxChi2 = 24.; //____________________________________________________________________ -AliTRDtracker::AliTRDtracker() +AliTRDtracker::AliTRDtracker(const TFile *geomfile) { - // - // Default constructor - // + // + // Main constructor + // - fEvent = 0; - fGeom = NULL; + Float_t fTzero = 0; + Int_t version = 2; + + fAddTRDseeds = kFALSE; - fNclusters = 0; - fClusters = NULL; - fNseeds = 0; - fSeeds = NULL; - fNtracks = 0; - fTracks = NULL; + fGeom = NULL; + + TDirectory *savedir=gDirectory; + TFile *in=(TFile*)geomfile; + if (!in->IsOpen()) { + printf("AliTRDtracker::AliTRDtracker(): geometry file is not open!\n"); + printf(" DETAIL TRD geometry and DEFAULT TRD parameter will be used\n"); + } + else { + in->cd(); + in->ls(); + fGeom = (AliTRDgeometry*) in->Get("TRDgeometry"); + fPar = (AliTRDparameter*) in->Get("TRDparameter"); + fGeom->Dump(); + } - fSY2corr = 0.025; -} + if(fGeom) { + // fTzero = geo->GetT0(); + fTzero = 0.; + version = fGeom->IsVersion(); + printf("Found geometry version %d on file \n", version); + } + else { + printf("AliTRDtracker::AliTRDtracker(): cann't find TRD geometry!\n"); + printf(" DETAIL TRD geometry and DEFAULT TRD parameter will be used\n"); + fGeom = new AliTRDgeometryDetail(); + fPar = new AliTRDparameter(); + } + + savedir->cd(); -//____________________________________________________________________ -AliTRDtracker::AliTRDtracker(const Text_t* name, const Text_t* title) - :TNamed(name, title) -{ - // - // TRD tracker contructor - // + + // fGeom->SetT0(fTzero); fEvent = 0; - fGeom = NULL; fNclusters = 0; fClusters = new TObjArray(2000); fNseeds = 0; - fSeeds = new TObjArray(20000); + fSeeds = new TObjArray(2000); fNtracks = 0; - fTracks = new TObjArray(10000); - - fSY2corr = 0.025; -} - -//___________________________________________________________________ -AliTRDtracker::~AliTRDtracker() -{ - // - // Destructor - // - - delete fClusters; - delete fTracks; - delete fSeeds; - delete fGeom; - -} - -//___________________________________________________________________ -void AliTRDtracker::Clusters2Tracks(TH1F *hs, TH1F *hd) -{ - // - // Do the trackfinding - // - - Int_t i; - - Int_t inner, outer; - Int_t fTotalNofTimeBins = fGeom->GetTimeMax() * AliTRDgeometry::Nplan(); - Int_t nSteps = (Int_t) (fgkSeedDepth / fgkSeedStep); - Int_t gap = (Int_t) (fTotalNofTimeBins * fgkSeedGap); - Int_t step = (Int_t) (fTotalNofTimeBins * fgkSeedStep); - + fTracks = new TObjArray(1000); - // nSteps = 1; - - if (!fClusters) return; - - AliTRDtrackingSector fTrSec[AliTRDgeometry::kNsect]; - SetUpSectors(fTrSec); + for(Int_t geom_s = 0; geom_s < kTRACKING_SECTORS; geom_s++) { + Int_t tr_s = CookSectorIndex(geom_s); + fTrSec[tr_s] = new AliTRDtrackingSector(fGeom, geom_s, fPar); + } - // make a first turn looking for seed ends in the same (n,n) - // and in the adjacent sectors (n,n+1) + fSY2corr = 0.025; + fSZ2corr = 12.; - for(i=0; iCamHght(); // Amplification region + Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region - for(i=0; iGetTimeBinSize(); + Int_t tbAmp = fPar->GetTimeBefore(); + Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx); + Int_t tbDrift = fPar->GetTimeMax(); + Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx); -} + tbDrift = TMath::Min(tbDrift,maxDrift); + tbAmp = TMath::Min(tbAmp,maxAmp); -//_____________________________________________________________________ -Double_t AliTRDtracker::ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt) const -{ - // - // Parametrised "expected" error of the cluster reconstruction in Y - // + fTimeBinsPerPlane = tbAmp + tbDrift; + fMaxGap = (Int_t) (fTimeBinsPerPlane * fGeom->Nplan() * fSkipDepth); - Double_t s = 0.08 * 0.08; - return s; + fVocal = kFALSE; -} +} -//_____________________________________________________________________ -Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t r, Double_t tgl) const +//___________________________________________________________________ +AliTRDtracker::~AliTRDtracker() { - // - // Parametrised "expected" error of the cluster reconstruction in Z - // - - Double_t s = 6 * 6 /12.; - return s; + delete fClusters; + delete fTracks; + delete fSeeds; + delete fGeom; + delete fPar; -} + for(Int_t geom_s = 0; geom_s < kTRACKING_SECTORS; geom_s++) { + delete fTrSec[geom_s]; + } +} //_____________________________________________________________________ -Double_t f1trd(Double_t x1,Double_t y1, - Double_t x2,Double_t y2, - Double_t x3,Double_t y3) +inline Double_t f1trd(Double_t x1,Double_t y1, + Double_t x2,Double_t y2, + Double_t x3,Double_t y3) { // // Initial approximation of the track curvature - // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch // - Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1); Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)- (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2)); @@ -245,17 +212,16 @@ Double_t f1trd(Double_t x1,Double_t y1, Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b); return -xr*yr/sqrt(xr*xr+yr*yr); - } //_____________________________________________________________________ -Double_t f2trd(Double_t x1,Double_t y1, - Double_t x2,Double_t y2, - Double_t x3,Double_t y3) +inline Double_t f2trd(Double_t x1,Double_t y1, + Double_t x2,Double_t y2, + Double_t x3,Double_t y3) { // - // Initial approximation of the track curvature times center of curvature - // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch + // Initial approximation of the track curvature times X coordinate + // of the center of curvature // Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1); @@ -267,472 +233,1167 @@ Double_t f2trd(Double_t x1,Double_t y1, Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b); return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr); - } //_____________________________________________________________________ -Double_t f3trd(Double_t x1,Double_t y1, - Double_t x2,Double_t y2, - Double_t z1,Double_t z2) +inline Double_t f3trd(Double_t x1,Double_t y1, + Double_t x2,Double_t y2, + Double_t z1,Double_t z2) { // // Initial approximation of the tangent of the track dip angle - // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch // return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)); - } - //___________________________________________________________________ - -Int_t AliTRDtracker::FindProlongation(AliTRDtrack& t, AliTRDtrackingSector *sec, - Int_t s, Int_t rf, Int_t matchedIndex, - TH1F *hs, TH1F *hd) +Int_t AliTRDtracker::Clusters2Tracks(const TFile *inp, TFile *out) { // - // Starting from current position on track=t this function tries - // to extrapolate the track up to timeBin=rf and to confirm prolongation - // if a close cluster is found. *sec is a pointer to allocated - // array of sectors, in which the initial sector has index=s. + // Finds tracks within the TRD. File is expected to contain seeds + // at the outer part of the TRD. If is NULL, the seeds + // are found within the TRD if fAddTRDseeds is TRUE. + // The tracks are propagated to the innermost time bin + // of the TRD and stored in file . // - // TH1F *hsame = hs; - // TH1F *hdiff = hd; + LoadEvent(); + + TDirectory *savedir=gDirectory; - // Bool_t good_match; + char tname[100]; - const Int_t kTimeBinsToSkip=Int_t(fgkSkipDepth*sec->GetNtimeBins()); - Int_t tryAgain=kTimeBinsToSkip; + if (!out->IsOpen()) { + cerr<<"AliTRDtracker::Clusters2Tracks(): output file is not open !\n"; + return 1; + } - Double_t alpha=AliTRDgeometry::GetAlpha(); + sprintf(tname,"seedTRDtoTPC_%d",fEvent); + TTree tpc_tree(tname,"Tree with seeds from TRD at outer TPC pad row"); + AliTPCtrack *iotrack=0; + tpc_tree.Branch("tracks","AliTPCtrack",&iotrack,32000,0); + + sprintf(tname,"TreeT%d_TRD",fEvent); + TTree trd_tree(tname,"TRD tracks at inner TRD time bin"); + AliTRDtrack *iotrack_trd=0; + trd_tree.Branch("tracks","AliTRDtrack",&iotrack_trd,32000,0); + + Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins(); + Float_t foundMin = fMinClustersInTrack * timeBins; + + if (inp) { + TFile *in=(TFile*)inp; + if (!in->IsOpen()) { + cerr<<"AliTRDtracker::Clusters2Tracks(): file with seeds is not open !\n"; + cerr<<" ... going for seeds finding inside the TRD\n"; + } + else { + in->cd(); + sprintf(tname,"TRDb_%d",fEvent); + TTree *seedTree=(TTree*)in->Get(tname); + if (!seedTree) { + cerr<<"AliTRDtracker::Clusters2Tracks(): "; + cerr<<"can't get a tree with track seeds !\n"; + return 3; + } + AliTRDtrack *seed=new AliTRDtrack; + seedTree->SetBranchAddress("tracks",&seed); + + Int_t n=(Int_t)seedTree->GetEntries(); + for (Int_t i=0; iGetEvent(i); + seed->ResetCovariance(); + fSeeds->AddLast(new AliTRDtrack(*seed)); + fNseeds++; + } + delete seed; + } + } - Int_t ns=Int_t(2*TMath::Pi()/alpha+0.5); + out->cd(); - Double_t x0, rho; + // find tracks from loaded seeds - for (Int_t nr=sec->GetTimeBinNumber(t.GetX())-1; nr>=rf; nr--) { + Int_t nseed=fSeeds->GetEntriesFast(); + Int_t i, found = 0; + Int_t innerTB = fTrSec[0]->GetInnerTimeBin(); + + for (i=0; iUncheckedAt(i), &t=*pt; + FollowProlongation(t, innerTB); + if (t.GetNumberOfClusters() >= foundMin) { + UseClusters(&t); + CookLabel(pt, 1-fLabelFraction); + // t.CookdEdx(); + } + iotrack_trd = pt; + trd_tree.Fill(); + cerr<GetAlpha()); + iotrack = tpc; + tpc_tree.Fill(); + delete tpc; + } + delete fSeeds->RemoveAt(i); + fNseeds--; + } - Double_t x=sec->GetX(nr); - Double_t ymax=x*TMath::Tan(0.5*alpha); + cout<<"Number of loaded seeds: "<TECframe(nr,t.GetY(),t.GetZ())) { - rho = 1.7; x0 = 33.0; // G10 frame - } - if(TMath::Abs(x - t.GetX()) > 3) { - rho = 0.0559; x0 = 55.6; // radiator - } - if (!t.PropagateTo(x,x0,rho,0.139)) { - cerr<<"Can't propagate to x = "<GetNumberOfTimeBins(); + + Int_t nSteps = (Int_t) (fSeedDepth / fSeedStep); + Int_t gap = (Int_t) (timeBins * fSeedGap); + Int_t step = (Int_t) (timeBins * fSeedStep); + + // make a first turn with tight cut on initial curvature + for(Int_t turn = 1; turn <= 2; turn++) { + if(turn == 2) { + nSteps = (Int_t) (fSeedDepth / (3*fSeedStep)); + step = (Int_t) (timeBins * (3*fSeedStep)); + } + for(Int_t i=0; ifgkWideRoad) { - if (t.GetNclusters() > 4) { - cerr<GetEntriesFast(); + printf("\n initial number of seeds %d\n", nseed); + + MakeSeeds(inner, outer, turn); + + nseed=fSeeds->GetEntriesFast(); + printf("\n number of seeds after MakeSeeds %d\n", nseed); + + + for (Int_t i=0; iUncheckedAt(i), &t=*pt; + FollowProlongation(t,innerTB); + if (t.GetNumberOfClusters() >= foundMin) { + UseClusters(&t); + CookLabel(pt, 1-fLabelFraction); + t.CookdEdx(); + cerr<GetAlpha()); + iotrack = tpc; + tpc_tree.Fill(); + delete tpc; + } + } + delete fSeeds->RemoveAt(i); + fNseeds--; + } } - return 0; - } + } + } + tpc_tree.Write(); + trd_tree.Write(); + + cout<<"Total number of found tracks: "<cd(); + + return 0; +} + + + +//_____________________________________________________________________________ +Int_t AliTRDtracker::PropagateBack(const TFile *inp, TFile *out) { + // + // Reads seeds from file . The seeds are AliTPCtrack's found and + // backpropagated by the TPC tracker. Each seed is first propagated + // to the TRD, and then its prolongation is searched in the TRD. + // If sufficiently long continuation of the track is found in the TRD + // the track is updated, otherwise it's stored as originaly defined + // by the TPC tracker. + // + + + LoadEvent(); + + TDirectory *savedir=gDirectory; + + TFile *in=(TFile*)inp; + + if (!in->IsOpen()) { + cerr<<"AliTRDtracker::PropagateBack(): "; + cerr<<"file with back propagated TPC tracks is not open !\n"; + return 1; + } + + if (!out->IsOpen()) { + cerr<<"AliTRDtracker::PropagateBack(): "; + cerr<<"file for back propagated TRD tracks is not open !\n"; + return 2; + } + + in->cd(); + char tname[100]; + sprintf(tname,"seedsTPCtoTRD_%d",fEvent); + TTree *seedTree=(TTree*)in->Get(tname); + if (!seedTree) { + cerr<<"AliTRDtracker::PropagateBack(): "; + cerr<<"can't get a tree with seeds from TPC !\n"; + cerr<<"check if your version of TPC tracker creates tree "<SetBranchAddress("tracks",&seed); + + Int_t n=(Int_t)seedTree->GetEntries(); + for (Int_t i=0; iGetEvent(i); + Int_t lbl = seed->GetLabel(); + AliTRDtrack *tr = new AliTRDtrack(*seed,seed->GetAlpha()); + tr->SetSeedLabel(lbl); + fSeeds->AddLast(tr); + fNseeds++; + } - Double_t maxChi2=fgkMaxChi2; + delete seed; + delete seedTree; - if (timeBin) { + out->cd(); - for (Int_t i=timeBin.Find(y-road); iGetTrackIndex(0) == matchedIndex) || - // (c->GetTrackIndex(1) == matchedIndex) || - // (c->GetTrackIndex(2) == matchedIndex)) good_match = kTRUE; - // if(hsame) hsame->Fill(TMath::Abs(c->GetY()-y)/road); - // if(hdiff) hdiff->Fill(road); + sprintf(tname,"seedsTRDtoTOF1_%d",fEvent); + TTree tofTree1(tname,"Tracks back propagated through TPC and TRD"); + tofTree1.Branch("tracks","AliTPCtrack",&otrack,32000,0); - if (c->GetY() > y+road) break; - if (c->IsUsed() > 0) continue; + sprintf(tname,"seedsTRDtoTOF2_%d",fEvent); + TTree tofTree2(tname,"Tracks back propagated through TPC and TRD"); + tofTree2.Branch("tracks","AliTPCtrack",&otrack,32000,0); - // if(good_match) hsame->Fill(TMath::Abs(c->GetZ()-z)); - // else hdiff->Fill(TMath::Abs(c->GetZ()-z)); + sprintf(tname,"seedsTRDtoPHOS_%d",fEvent); + TTree phosTree(tname,"Tracks back propagated through TPC and TRD"); + phosTree.Branch("tracks","AliTPCtrack",&otrack,32000,0); - // if(!good_match) continue; + sprintf(tname,"seedsTRDtoRICH_%d",fEvent); + TTree richTree(tname,"Tracks back propagated through TPC and TRD"); + richTree.Branch("tracks","AliTPCtrack",&otrack,32000,0); - if((c->GetZ()-z)*(c->GetZ()-z) > 3 * 12 * sz2) continue; + sprintf(tname,"TRDb_%d",fEvent); + TTree trdTree(tname,"Back propagated TRD tracks at outer TRD time bin"); + AliTRDtrack *otrack_trd=0; + trdTree.Branch("tracks","AliTRDtrack",&otrack_trd,32000,0); + + Int_t found=0; - Double_t chi2=t.GetPredictedChi2(c); + Int_t nseed=fSeeds->GetEntriesFast(); - // if((c->GetTrackIndex(0) == matchedIndex) || - // (c->GetTrackIndex(1) == matchedIndex) || - // (c->GetTrackIndex(2) == matchedIndex)) - // hdiff->Fill(chi2); + Float_t foundMin = fMinClustersInTrack * fTimeBinsPerPlane * fGeom->Nplan(); - // ncl++; + Int_t outermost_tb = fTrSec[0]->GetOuterTimeBin(); - if (chi2 > maxChi2) continue; - maxChi2=chi2; - cl=c; - index=timeBin.GetIndex(i); - } - - if(!cl) { + for (Int_t i=0; iUncheckedAt(i), &s=*ps; + Int_t expectedClr = FollowBackProlongation(s); + Int_t foundClr = s.GetNumberOfClusters(); + Int_t last_tb = fTrSec[0]->GetLayerNumber(s.GetX()); - if (c->GetY() > y+road) break; - if (c->IsUsed() > 0) continue; - if((c->GetZ()-z)*(c->GetZ()-z) > 3.25 * 12 * sz2) continue; - - Double_t chi2=t.GetPredictedChi2(c); - - // ncl++; + printf("seed %d: found %d out of %d expected clusters, Min is %f\n", + i, foundClr, expectedClr, foundMin); - if (chi2 > maxChi2) continue; - maxChi2=chi2; - cl=c; - index=timeBin.GetIndex(i); - } - } - + if (foundClr >= foundMin) { + s.CookdEdx(); + CookLabel(ps, 1-fLabelFraction); + UseClusters(ps); + otrack_trd=ps; + trdTree.Fill(); + cout< ymax) { - // cerr<<"y > ymax: "< "<GetListOfFiles()->FindObject(hitfile); - if (!fInputFile) { - printf("AliTRDtracker::Open -- "); - printf("Open the ALIROOT-file %s.\n",hitfile); - fInputFile = new TFile(hitfile,"UPDATE"); - } - else { - printf("AliTRDtracker::Open -- "); - printf("%s is already open.\n",hitfile); - } + Float_t wIndex, wTB, wChi2; + Float_t wYrt, wYclosest, wYcorrect, wYwindow; + Float_t wZrt, wZclosest, wZcorrect, wZwindow; + Float_t wPx, wPy, wPz, wC; + Double_t Px, Py, Pz; + Float_t wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2; - // Get AliRun object from file or create it if not on file + Int_t trackIndex = t.GetLabel(); - gAlice = (AliRun*) fInputFile->Get("gAlice"); - if (gAlice) { - printf("AliTRDtracker::GetEvent -- "); - printf("AliRun object found on file.\n"); - } - else { - printf("AliTRDtracker::GetEvent -- "); - printf("Could not find AliRun object.\n"); - } + Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5); - /* - // Import the Trees for the event nEvent in the file - Int_t nparticles = gAlice->GetEvent(fEvent); - cerr<<"nparticles = "<GetDetector("TRD"); - fGeom = trd->GetGeometry(); + if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); + if (alpha < 0. ) alpha += 2.*TMath::Pi(); -} + Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect; + Double_t x0, rho, x, dx, y, ymax, z; -//_____________________________________________________________________________ -void AliTRDtracker::SetUpSectors(AliTRDtrackingSector *sec) -{ - // - // Fills clusters into TRD tracking_sectors - // Note that the numbering scheme for the TRD tracking_sectors - // differs from that of TRD sectors - // + Int_t expectedNumberOfClusters = 0; + Bool_t lookForCluster; - for (Int_t i=0; iGetLayerNumber(t.GetX()); nr>rf; nr--) { + + y = t.GetY(); z = t.GetZ(); + + // first propagate to the inner surface of the current time bin + fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster); + x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; y = t.GetY(); z = t.GetZ(); + if(!t.PropagateTo(x,x0,rho,0.139)) break; + y = t.GetY(); + ymax = x*TMath::Tan(0.5*alpha); + if (y > ymax) { + s = (s+1) % ns; + if (!t.Rotate(alpha)) break; + if(!t.PropagateTo(x,x0,rho,0.139)) break; + } else if (y <-ymax) { + s = (s-1+ns) % ns; + if (!t.Rotate(-alpha)) break; + if(!t.PropagateTo(x,x0,rho,0.139)) break; + } - cerr<<"SetUpSectors: sorting clusters"<GetEntriesFast(); - UInt_t index; - while (ncl--) { - printf("\r %d left ",ncl); - AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl); - Int_t detector=c->GetDetector(), localTimeBin=c->GetLocalTimeBin(); - Int_t sector=fGeom->GetSector(detector); + y = t.GetY(); z = t.GetZ(); + + // now propagate to the middle plane of the next time bin + fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster); + x = fTrSec[s]->GetLayer(nr-1)->GetX(); y = t.GetY(); z = t.GetZ(); + if(!t.PropagateTo(x,x0,rho,0.139)) break; + y = t.GetY(); + ymax = x*TMath::Tan(0.5*alpha); + if (y > ymax) { + s = (s+1) % ns; + if (!t.Rotate(alpha)) break; + if(!t.PropagateTo(x,x0,rho,0.139)) break; + } else if (y <-ymax) { + s = (s-1+ns) % ns; + if (!t.Rotate(-alpha)) break; + if(!t.PropagateTo(x,x0,rho,0.139)) break; + } - Int_t trackingSector = AliTRDgeometry::kNsect - sector - 1; - Int_t tb=sec[sector].GetTimeBin(detector,localTimeBin); - index=ncl; - sec[trackingSector][tb].InsertCluster(c,index); + if(lookForCluster) { - } - printf("\r\n"); -} + expectedNumberOfClusters++; -//_____________________________________________________________________________ -void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer, - AliTRDtrackingSector* fTrSec, Int_t turn, - TH1F *hs, TH1F *hd) -{ - // - // Creates track seeds using clusters in timeBins=i1,i2 - // + wIndex = (Float_t) t.GetLabel(); + wTB = nr; - Int_t i2 = inner, i1 = outer; - Int_t ti[3], to[3]; - Int_t nprim = 85600/2; + AliTRDpropagationLayer& time_bin=*(fTrSec[s]->GetLayer(nr-1)); - TH1F *hsame = hs; - TH1F *hdiff = hd; - Bool_t match = false; - Int_t matchedIndex; - // find seeds + Double_t sy2=ExpectedSigmaY2(x,t.GetTgl(),t.GetPt()); - Double_t x[5], c[15]; - Int_t maxSec=AliTRDgeometry::kNsect; - Double_t alpha=AliTRDgeometry::GetAlpha(); - Double_t shift=AliTRDgeometry::GetAlpha()/2.; - Double_t cs=cos(alpha), sn=sin(alpha); - Double_t cs2=cos(2.*alpha), sn2=sin(2.*alpha); + Double_t sz2=ExpectedSigmaZ2(x,t.GetTgl()); - Double_t x1 =fTrSec[0].GetX(i1); - Double_t xx2=fTrSec[0].GetX(i2); + Double_t road; + if((t.GetSigmaY2() + sy2) > 0) road=10.*sqrt(t.GetSigmaY2() + sy2); + else return expectedNumberOfClusters; + + wYrt = (Float_t) y; + wZrt = (Float_t) z; + wYwindow = (Float_t) road; + t.GetPxPyPz(Px,Py,Pz); + wPx = (Float_t) Px; + wPy = (Float_t) Py; + wPz = (Float_t) Pz; + wC = (Float_t) t.GetC(); + wSigmaC2 = (Float_t) t.GetSigmaC2(); + wSigmaTgl2 = (Float_t) t.GetSigmaTgl2(); + wSigmaY2 = (Float_t) t.GetSigmaY2(); + wSigmaZ2 = (Float_t) t.GetSigmaZ2(); + wChi2 = -1; + + if (road>fWideRoad) { + if (t.GetNumberOfClusters()>4) + cerr<GetLabel(0) != trackIndex) && + (c->GetLabel(1) != trackIndex) && + (c->GetLabel(2) != trackIndex)) continue; + if(TMath::Abs(c->GetY() - y) > minDY) continue; + minDY = TMath::Abs(c->GetY() - y); + wYcorrect = c->GetY(); + wZcorrect = c->GetZ(); + wChi2 = t.GetPredictedChi2(c); + } + } + // Now go for the real cluster search - printf("\n"); + if (time_bin) { - if((turn != 1)&&(turn != 2)) { - printf("*** Error in MakeSeeds: unexpected turn = %d \n", turn); - return; - } + for (Int_t i=time_bin.Find(y-road); iGetY() > y+road) break; + if (c->IsUsed() > 0) continue; + if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue; + Double_t chi2=t.GetPredictedChi2(c); + + if (chi2 > max_chi2) continue; + max_chi2=chi2; + cl=c; + index=time_bin.GetIndex(i); + } + + if(!cl) { + + for (Int_t i=time_bin.Find(y-road); iGetY() > y+road) break; + if (c->IsUsed() > 0) continue; + if((c->GetZ()-z)*(c->GetZ()-z) > 2.25 * 12 * sz2) continue; + + Double_t chi2=t.GetPredictedChi2(c); + + if (chi2 > max_chi2) continue; + max_chi2=chi2; + cl=c; + index=time_bin.GetIndex(i); + } + } + + if (cl) { + wYclosest = cl->GetY(); + wZclosest = cl->GetZ(); - for (Int_t ns=0; nsGetQ()/dx,t.GetNumberOfClusters()); + if(!t.Update(cl,max_chi2,index)) { + if(!try_again--) return 0; + } + else try_again=fMaxGap; + } + else { + if (try_again==0) break; + try_again--; + } - printf("\n MakeSeeds: sector %d \n", ns); + /* + if((((Int_t) wTB)%15 == 0) || (((Int_t) wTB)%15 == 14)) { + + printf(" %f", wIndex); //1 + printf(" %f", wTB); //2 + printf(" %f", wYrt); //3 + printf(" %f", wYclosest); //4 + printf(" %f", wYcorrect); //5 + printf(" %f", wYwindow); //6 + printf(" %f", wZrt); //7 + printf(" %f", wZclosest); //8 + printf(" %f", wZcorrect); //9 + printf(" %f", wZwindow); //10 + printf(" %f", wPx); //11 + printf(" %f", wPy); //12 + printf(" %f", wPz); //13 + printf(" %f", wSigmaC2*1000000); //14 + printf(" %f", wSigmaTgl2*1000); //15 + printf(" %f", wSigmaY2); //16 + // printf(" %f", wSigmaZ2); //17 + printf(" %f", wChi2); //17 + printf(" %f", wC); //18 + printf("\n"); + } + */ + } + } + } + return expectedNumberOfClusters; + + +} - Int_t nl2=fTrSec[(ns-2+maxSec)%maxSec][i2]; - Int_t nl=fTrSec[(ns-1+maxSec)%maxSec][i2]; - Int_t nm=fTrSec[ns][i2]; - Int_t nu=fTrSec[(ns+1)%maxSec][i2]; - Int_t nu2=fTrSec[(ns+2)%maxSec][i2]; +//___________________________________________________________________ - AliTRDtimeBin& r1=fTrSec[ns][i1]; +Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t) +{ + // Starting from current radial position of track this function + // extrapolates the track up to outer timebin and in the sensitive + // layers confirms prolongation if a close cluster is found. + // Returns the number of clusters expected to be found in sensitive layers - for (Int_t is=0; is < r1; is++) { - Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ(); - for(Int_t ii=0; ii<3; ii++) to[ii] = r1[is]->GetTrackIndex(ii); + Float_t wIndex, wTB, wChi2; + Float_t wYrt, wYclosest, wYcorrect, wYwindow; + Float_t wZrt, wZclosest, wZcorrect, wZwindow; + Float_t wPx, wPy, wPz, wC; + Double_t Px, Py, Pz; + Float_t wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2; - for (Int_t js=0; js < nl2+nl+nm+nu+nu2; js++) { - - const AliTRDcluster *cl; - Double_t x2, y2, z2; - Double_t x3=0., y3=0.; + Int_t trackIndex = t.GetLabel(); - if (jsGetY(); z2=cl->GetZ(); - for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii); + Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5); - x2= xx2*cs2+y2*sn2; - y2=-xx2*sn2+y2*cs2; - } - else if (jsGetY(); z2=cl->GetZ(); - for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii); + Int_t try_again=fMaxGap; - x2= xx2*cs+y2*sn; - y2=-xx2*sn+y2*cs; + Double_t alpha=t.GetAlpha(); - } - else if (jsGetY(); z2=cl->GetZ(); - for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii); - } - else if (jsGetY(); z2=cl->GetZ(); - for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii); + if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); + if (alpha < 0. ) alpha += 2.*TMath::Pi(); - x2=xx2*cs-y2*sn; - y2=xx2*sn+y2*cs; + Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect; - } - else { - if(turn != 2) continue; - AliTRDtimeBin& r2=fTrSec[(ns+2)%maxSec][i2]; - cl=r2[js-nl2-nl-nm-nu]; - y2=cl->GetY(); z2=cl->GetZ(); - for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii); - - x2=xx2*cs2-y2*sn2; - y2=xx2*sn2+y2*cs2; - } - + Int_t outerTB = fTrSec[0]->GetOuterTimeBin(); - match = false; - matchedIndex = -1; - for (Int_t ii=0; ii<3; ii++) { - // cerr<<"ti["< fgkMaxSeedDeltaZ12) continue; + Double_t x0, rho, x, dx, y, ymax, z; - Double_t zz=z1 - z1/x1*(x1-x2); + Bool_t lookForCluster; - if (TMath::Abs(zz-z2)>fgkMaxSeedDeltaZ) continue; + Int_t expectedNumberOfClusters = 0; + x = t.GetX(); - Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1); - if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;} + alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning - x[0]=y1; - x[1]=z1; - x[2]=f1trd(x1,y1,x2,y2,x3,y3); - if (TMath::Abs(x[2]) > fgkMaxSeedC) continue; + for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr= 0.99999) continue; + // first propagate to the outer surface of the current time bin - x[4]=f3trd(x1,y1,x2,y2,z1,z2); + fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster); + x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2; y = t.GetY(); z = t.GetZ(); - if (TMath::Abs(x[4]) > fgkMaxSeedTan) continue; + if(!t.PropagateTo(x,x0,rho,0.139)) break; - Double_t a=asin(x[3]); - Double_t zv=z1 - x[4]/x[2]*(a+asin(x[2]*x1-x[3])); + y = t.GetY(); + ymax = x*TMath::Tan(0.5*alpha); - if (TMath::Abs(zv)>fgkMaxSeedVertexZ) continue; + if (y > ymax) { + s = (s+1) % ns; + if (!t.Rotate(alpha)) break; + if(!t.PropagateTo(x,x0,rho,0.139)) break; + } else if (y <-ymax) { + s = (s-1+ns) % ns; + if (!t.Rotate(-alpha)) break; + if(!t.PropagateTo(x,x0,rho,0.139)) break; + } + y = t.GetY(); z = t.GetZ(); - Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2(); - Double_t sy2=cl->GetSigmaY2(), sz2=cl->GetSigmaZ2(); - Double_t sy3=fgkSeedErrorSY3, sy=fgkSeedErrorSY, sz=fgkSeedErrorSZ; + // now propagate to the middle plane of the next time bin + fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster); - Double_t f20=(f1trd(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy; - Double_t f22=(f1trd(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy; - Double_t f24=(f1trd(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy; - Double_t f30=(f2trd(x1,y1+sy,x2,y2,x3,y3)-x[3])/sy; - Double_t f32=(f2trd(x1,y1,x2,y2+sy,x3,y3)-x[3])/sy; - Double_t f34=(f2trd(x1,y1,x2,y2,x3,y3+sy)-x[3])/sy; - Double_t f40=(f3trd(x1,y1+sy,x2,y2,z1,z2)-x[4])/sy; - Double_t f41=(f3trd(x1,y1,x2,y2,z1+sz,z2)-x[4])/sz; - Double_t f42=(f3trd(x1,y1,x2,y2+sy,z1,z2)-x[4])/sy; - Double_t f43=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[4])/sz; + x = fTrSec[s]->GetLayer(nr+1)->GetX(); y = t.GetY(); z = t.GetZ(); - c[0]=sy1; - c[1]=0.; c[2]=sz1; - c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f24*sy3*f24; - c[6]=f30*sy1; c[7]=0.; c[8]=f30*sy1*f20+f32*sy2*f22+f34*sy3*f24; - c[9]=f30*sy1*f30+f32*sy2*f32+f34*sy3*f34; - c[10]=f40*sy1; c[11]=f41*sz1; c[12]=f40*sy1*f20+f42*sy2*f22; - c[13]=f40*sy1*f30+f42*sy2*f32; - c[14]=f40*sy1*f40+f41*sz1*f41+f42*sy2*f42+f43*sz2*f43; + if(!t.PropagateTo(x,x0,rho,0.139)) break; - UInt_t index=r1.GetIndex(is); - - AliTRDtrack *track=new AliTRDtrack(r1[is],index,x,c,x1,ns*alpha+shift); + y = t.GetY(); - Int_t rc=FindProlongation(*track,fTrSec,ns,i2,matchedIndex,hsame,hdiff); + ymax = x*TMath::Tan(0.5*alpha); - // if (match) hsame->Fill((Float_t) track->GetNclusters()); - // else hdiff->Fill((Float_t) track->GetNclusters()); - // delete track; - // continue; + if(fVocal) printf("nr+1=%d, x %f, z %f, y %f, ymax %f\n",nr+1,x,z,y,ymax); + + if (y > ymax) { + s = (s+1) % ns; + if (!t.Rotate(alpha)) break; + if(!t.PropagateTo(x,x0,rho,0.139)) break; + } else if (y <-ymax) { + s = (s-1+ns) % ns; + if (!t.Rotate(-alpha)) break; + if(!t.PropagateTo(x,x0,rho,0.139)) break; + } - if ((rc < 0) || - (track->GetNclusters() < (i1-i2)*fgkMinClustersInSeed)) delete track; - else { - fSeeds->AddLast(track); fNseeds++; - printf("\r found seed %d ", fNseeds); + // printf("label %d, pl %d, lookForCluster %d \n", + // trackIndex, nr+1, lookForCluster); + + if(lookForCluster) { + expectedNumberOfClusters++; + + wIndex = (Float_t) t.GetLabel(); + wTB = fTrSec[s]->GetLayer(nr+1)->GetTimeBinIndex(); + + AliTRDpropagationLayer& time_bin=*(fTrSec[s]->GetLayer(nr+1)); + Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt()); + Double_t sz2=ExpectedSigmaZ2(t.GetX(),t.GetTgl()); + if((t.GetSigmaY2() + sy2) < 0) break; + Double_t road = 10.*sqrt(t.GetSigmaY2() + sy2); + Double_t y=t.GetY(), z=t.GetZ(); + + wYrt = (Float_t) y; + wZrt = (Float_t) z; + wYwindow = (Float_t) road; + t.GetPxPyPz(Px,Py,Pz); + wPx = (Float_t) Px; + wPy = (Float_t) Py; + wPz = (Float_t) Pz; + wC = (Float_t) t.GetC(); + wSigmaC2 = (Float_t) t.GetSigmaC2(); + wSigmaTgl2 = (Float_t) t.GetSigmaTgl2(); + wSigmaY2 = (Float_t) t.GetSigmaY2(); + wSigmaZ2 = (Float_t) t.GetSigmaZ2(); + wChi2 = -1; + + if (road>fWideRoad) { + if (t.GetNumberOfClusters()>4) + cerr<GetLabel(0) != trackIndex) && + (c->GetLabel(1) != trackIndex) && + (c->GetLabel(2) != trackIndex)) continue; + if(TMath::Abs(c->GetY() - y) > minDY) continue; + minDY = TMath::Abs(c->GetY() - y); + wYcorrect = c->GetY(); + wZcorrect = c->GetZ(); + wChi2 = t.GetPredictedChi2(c); } - } - } - } + } - fSeeds->Sort(); -} + // Now go for the real cluster search -//_____________________________________________________________________________ -void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename) -{ - // + if (time_bin) { + + for (Int_t i=time_bin.Find(y-road); iGetY() > y+road) break; + if (c->IsUsed() > 0) continue; + if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue; + Double_t chi2=t.GetPredictedChi2(c); + + if (chi2 > max_chi2) continue; + max_chi2=chi2; + cl=c; + index=time_bin.GetIndex(i); + } + + if(!cl) { + + for (Int_t i=time_bin.Find(y-road); iGetY() > y+road) break; + if (c->IsUsed() > 0) continue; + if((c->GetZ()-z)*(c->GetZ()-z) > 2.25 * 12 * sz2) continue; + + Double_t chi2=t.GetPredictedChi2(c); + + if (chi2 > max_chi2) continue; + max_chi2=chi2; + cl=c; + index=time_bin.GetIndex(i); + } + } + + if (cl) { + wYclosest = cl->GetY(); + wZclosest = cl->GetZ(); + + t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters()); + if(!t.Update(cl,max_chi2,index)) { + if(!try_again--) return 0; + } + else try_again=fMaxGap; + } + else { + if (try_again==0) break; + try_again--; + } + + /* + if((((Int_t) wTB)%15 == 0) || (((Int_t) wTB)%15 == 14)) { + + printf(" %f", wIndex); //1 + printf(" %f", wTB); //2 + printf(" %f", wYrt); //3 + printf(" %f", wYclosest); //4 + printf(" %f", wYcorrect); //5 + printf(" %f", wYwindow); //6 + printf(" %f", wZrt); //7 + printf(" %f", wZclosest); //8 + printf(" %f", wZcorrect); //9 + printf(" %f", wZwindow); //10 + printf(" %f", wPx); //11 + printf(" %f", wPy); //12 + printf(" %f", wPz); //13 + printf(" %f", wSigmaC2*1000000); //14 + printf(" %f", wSigmaTgl2*1000); //15 + printf(" %f", wSigmaY2); //16 + // printf(" %f", wSigmaZ2); //17 + printf(" %f", wChi2); //17 + printf(" %f", wC); //18 + printf("\n"); + } + */ + } + } + } + return expectedNumberOfClusters; +} + +//___________________________________________________________________ + +Int_t AliTRDtracker::PropagateToOuterPlane(AliTRDtrack& t, Double_t xToGo) +{ + // Starting from current radial position of track this function + // extrapolates the track up to radial position . + // Returns 1 if track reaches the plane, and 0 otherwise + + Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5); + + Double_t alpha=t.GetAlpha(); + + if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); + if (alpha < 0. ) alpha += 2.*TMath::Pi(); + + Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect; + + Bool_t lookForCluster; + Double_t x0, rho, x, dx, y, ymax, z; + + x = t.GetX(); + + alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning + + Int_t plToGo = fTrSec[0]->GetLayerNumber(xToGo); + + for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nrGetLayer(nr)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster); + x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2; y = t.GetY(); z = t.GetZ(); + if(!t.PropagateTo(x,x0,rho,0.139)) return 0; + y = t.GetY(); + ymax = x*TMath::Tan(0.5*alpha); + if (y > ymax) { + s = (s+1) % ns; + if (!t.Rotate(alpha)) return 0; + } else if (y <-ymax) { + s = (s-1+ns) % ns; + if (!t.Rotate(-alpha)) return 0; + } + if(!t.PropagateTo(x,x0,rho,0.139)) return 0; + + y = t.GetY(); z = t.GetZ(); + + // now propagate to the middle plane of the next time bin + fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster); + x = fTrSec[s]->GetLayer(nr+1)->GetX(); y = t.GetY(); z = t.GetZ(); + if(!t.PropagateTo(x,x0,rho,0.139)) return 0; + y = t.GetY(); + ymax = x*TMath::Tan(0.5*alpha); + if (y > ymax) { + s = (s+1) % ns; + if (!t.Rotate(alpha)) return 0; + } else if (y <-ymax) { + s = (s-1+ns) % ns; + if (!t.Rotate(-alpha)) return 0; + } + if(!t.PropagateTo(x,x0,rho,0.139)) return 0; + } + return 1; +} + +//___________________________________________________________________ + +Int_t AliTRDtracker::PropagateToTPC(AliTRDtrack& t) +{ + // Starting from current radial position of track this function + // extrapolates the track up to radial position of the outermost + // padrow of the TPC. + // Returns 1 if track reaches the TPC, and 0 otherwise + + Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5); + + Double_t alpha=t.GetAlpha(); + + if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); + if (alpha < 0. ) alpha += 2.*TMath::Pi(); + + Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect; + + Bool_t lookForCluster; + Double_t x0, rho, x, dx, y, ymax, z; + + x = t.GetX(); + + alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning + + Int_t plTPC = fTrSec[0]->GetLayerNumber(246.055); + + for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr>plTPC; nr--) { + + y = t.GetY(); z = t.GetZ(); + + // first propagate to the outer surface of the current time bin + fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster); + x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; y = t.GetY(); z = t.GetZ(); + if(!t.PropagateTo(x,x0,rho,0.139)) return 0; + y = t.GetY(); + ymax = x*TMath::Tan(0.5*alpha); + if (y > ymax) { + s = (s+1) % ns; + if (!t.Rotate(alpha)) return 0; + } else if (y <-ymax) { + s = (s-1+ns) % ns; + if (!t.Rotate(-alpha)) return 0; + } + if(!t.PropagateTo(x,x0,rho,0.139)) return 0; + + y = t.GetY(); z = t.GetZ(); + + // now propagate to the middle plane of the next time bin + fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster); + x = fTrSec[s]->GetLayer(nr-1)->GetX(); y = t.GetY(); z = t.GetZ(); + if(!t.PropagateTo(x,x0,rho,0.139)) return 0; + y = t.GetY(); + ymax = x*TMath::Tan(0.5*alpha); + if (y > ymax) { + s = (s+1) % ns; + if (!t.Rotate(alpha)) return 0; + } else if (y <-ymax) { + s = (s-1+ns) % ns; + if (!t.Rotate(-alpha)) return 0; + } + if(!t.PropagateTo(x,x0,rho,0.139)) return 0; + } + return 1; +} + + +//_____________________________________________________________________________ +void AliTRDtracker::LoadEvent() +{ + // Fills clusters into TRD tracking_sectors + // Note that the numbering scheme for the TRD tracking_sectors + // differs from that of TRD sectors + + ReadClusters(fClusters); + Int_t ncl=fClusters->GetEntriesFast(); + cout<<"\n LoadSectors: sorting "<UncheckedAt(ncl); + Int_t detector=c->GetDetector(), local_time_bin=c->GetLocalTimeBin(); + Int_t sector=fGeom->GetSector(detector); + Int_t plane=fGeom->GetPlane(detector); + + Int_t tracking_sector = CookSectorIndex(sector); + + Int_t gtb = fTrSec[tracking_sector]->CookTimeBinIndex(plane,local_time_bin); + Int_t layer = fTrSec[tracking_sector]->GetLayerNumber(gtb); + + index=ncl; + fTrSec[tracking_sector]->GetLayer(layer)->InsertCluster(c,index); + } + printf("\r\n"); + +} + +//_____________________________________________________________________________ +void AliTRDtracker::UnloadEvent() +{ + // + // Clears the arrays of clusters and tracks. Resets sectors and timebins + // + + Int_t i, nentr; + + nentr = fClusters->GetEntriesFast(); + for (i = 0; i < nentr; i++) delete fClusters->RemoveAt(i); + + nentr = fSeeds->GetEntriesFast(); + for (i = 0; i < nentr; i++) delete fSeeds->RemoveAt(i); + + nentr = fTracks->GetEntriesFast(); + for (i = 0; i < nentr; i++) delete fTracks->RemoveAt(i); + + Int_t nsec = AliTRDgeometry::kNsect; + + for (i = 0; i < nsec; i++) { + for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) { + fTrSec[i]->GetLayer(pl)->Clear(); + } + } + +} + +//__________________________________________________________________________ +void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer, Int_t turn) +{ + // Creates track seeds using clusters in timeBins=i1,i2 + + if(turn > 2) { + cerr<<"MakeSeeds: turn "<GetLayerNumber(inner); + Int_t i1 = fTrSec[0]->GetLayerNumber(outer); + + Double_t x1 =fTrSec[0]->GetX(i1); + Double_t xx2=fTrSec[0]->GetX(i2); + + for (Int_t ns=0; nsGetLayer(i2)); + Int_t nl=(*fTrSec[(ns-1+max_sec)%max_sec]->GetLayer(i2)); + Int_t nm=(*fTrSec[ns]->GetLayer(i2)); + Int_t nu=(*fTrSec[(ns+1)%max_sec]->GetLayer(i2)); + Int_t nu2=(*fTrSec[(ns+2)%max_sec]->GetLayer(i2)); + + AliTRDpropagationLayer& r1=*(fTrSec[ns]->GetLayer(i1)); + + for (Int_t is=0; is < r1; is++) { + Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ(); + + for (Int_t js=0; js < nl2+nl+nm+nu+nu2; js++) { + + const AliTRDcluster *cl; + Double_t x2, y2, z2; + Double_t x3=0., y3=0.; + + if (jsGetLayer(i2)); + cl=r2[js]; + y2=cl->GetY(); z2=cl->GetZ(); + + x2= xx2*cs2+y2*sn2; + y2=-xx2*sn2+y2*cs2; + } + else if (jsGetLayer(i2)); + cl=r2[js-nl2]; + y2=cl->GetY(); z2=cl->GetZ(); + + x2= xx2*cs+y2*sn; + y2=-xx2*sn+y2*cs; + } + else if (jsGetLayer(i2)); + cl=r2[js-nl2-nl]; + x2=xx2; y2=cl->GetY(); z2=cl->GetZ(); + } + else if (jsGetLayer(i2)); + cl=r2[js-nl2-nl-nm]; + y2=cl->GetY(); z2=cl->GetZ(); + + x2=xx2*cs-y2*sn; + y2=xx2*sn+y2*cs; + } + else { + if(turn != 2) continue; + AliTRDpropagationLayer& r2=*(fTrSec[(ns+2)%max_sec]->GetLayer(i2)); + cl=r2[js-nl2-nl-nm-nu]; + y2=cl->GetY(); z2=cl->GetZ(); + + x2=xx2*cs2-y2*sn2; + y2=xx2*sn2+y2*cs2; + } + + if(TMath::Abs(z1-z2) > fMaxSeedDeltaZ12) continue; + + Double_t zz=z1 - z1/x1*(x1-x2); + + if (TMath::Abs(zz-z2)>fMaxSeedDeltaZ) continue; + + Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1); + if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;} + + x[0]=y1; + x[1]=z1; + x[2]=f1trd(x1,y1,x2,y2,x3,y3); + + if (TMath::Abs(x[2]) > fMaxSeedC) continue; + + x[3]=f2trd(x1,y1,x2,y2,x3,y3); + + if (TMath::Abs(x[2]*x1-x[3]) >= 0.99999) continue; + + x[4]=f3trd(x1,y1,x2,y2,z1,z2); + + if (TMath::Abs(x[4]) > fMaxSeedTan) continue; + + Double_t a=asin(x[3]); + Double_t zv=z1 - x[4]/x[2]*(a+asin(x[2]*x1-x[3])); + + if (TMath::Abs(zv)>fMaxSeedVertexZ) continue; + + Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2(); + Double_t sy2=cl->GetSigmaY2(), sz2=cl->GetSigmaZ2(); + Double_t sy3=fSeedErrorSY3, sy=fSeedErrorSY, sz=fSeedErrorSZ; + + Double_t f20=(f1trd(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy; + Double_t f22=(f1trd(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy; + Double_t f24=(f1trd(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy; + Double_t f30=(f2trd(x1,y1+sy,x2,y2,x3,y3)-x[3])/sy; + Double_t f32=(f2trd(x1,y1,x2,y2+sy,x3,y3)-x[3])/sy; + Double_t f34=(f2trd(x1,y1,x2,y2,x3,y3+sy)-x[3])/sy; + Double_t f40=(f3trd(x1,y1+sy,x2,y2,z1,z2)-x[4])/sy; + Double_t f41=(f3trd(x1,y1,x2,y2,z1+sz,z2)-x[4])/sz; + Double_t f42=(f3trd(x1,y1,x2,y2+sy,z1,z2)-x[4])/sy; + Double_t f43=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[4])/sz; + + c[0]=sy1; + c[1]=0.; c[2]=sz1; + c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f24*sy3*f24; + c[6]=f30*sy1; c[7]=0.; c[8]=f30*sy1*f20+f32*sy2*f22+f34*sy3*f24; + c[9]=f30*sy1*f30+f32*sy2*f32+f34*sy3*f34; + c[10]=f40*sy1; c[11]=f41*sz1; c[12]=f40*sy1*f20+f42*sy2*f22; + c[13]=f40*sy1*f30+f42*sy2*f32; + c[14]=f40*sy1*f40+f41*sz1*f41+f42*sy2*f42+f43*sz2*f43; + + UInt_t index=r1.GetIndex(is); + + AliTRDtrack *track=new AliTRDtrack(r1[is],index,x,c,x1,ns*alpha+shift); + + Int_t rc=FollowProlongation(*track, i2); + + if ((rc < 1) || + (track->GetNumberOfClusters() < + (outer-inner)*fMinClustersInSeed)) delete track; + else { + fSeeds->AddLast(track); fNseeds++; + cerr<<"\r found seed "<= 0) or AliTRDrecPoints (option < 0) // from the file. The names of the cluster tree and branches // should match the ones used in AliTRDclusterizer::WriteClusters() @@ -740,155 +1401,169 @@ void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename) TDirectory *savedir=gDirectory; - TFile *file = TFile::Open(filename); - if (!file->IsOpen()) {printf("Can't open file %s !\n",filename); return;} + if (inp) { + TFile *in=(TFile*)inp; + if (!in->IsOpen()) { + cerr<<"AliTRDtracker::ReadClusters(): input file is not open !\n"; + return; + } + else{ + in->cd(); + } + } Char_t treeName[12]; sprintf(treeName,"TreeR%d_TRD",fEvent); - TTree *clusterTree = (TTree*) file->Get(treeName); - - TObjArray *clusterArray = new TObjArray(400); - - clusterTree->GetBranch("TRDcluster")->SetAddress(&clusterArray); + TTree *ClusterTree = (TTree*) gDirectory->Get(treeName); + + TObjArray *ClusterArray = new TObjArray(400); + + ClusterTree->GetBranch("TRDcluster")->SetAddress(&ClusterArray); + + Int_t nEntries = (Int_t) ClusterTree->GetEntries(); + printf("found %d entries in %s.\n",nEntries,ClusterTree->GetName()); - Int_t nEntries = (Int_t) clusterTree->GetEntries(); - printf("found %d entries in %s.\n",nEntries,clusterTree->GetName()); - // Loop through all entries in the tree Int_t nbytes; AliTRDcluster *c = 0; + printf("\n"); for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) { // Import the tree - nbytes += clusterTree->GetEvent(iEntry); - + nbytes += ClusterTree->GetEvent(iEntry); + // Get the number of points in the detector - Int_t nCluster = clusterArray->GetEntriesFast(); - printf("Read %d clusters from entry %d \n", nCluster, iEntry); - + Int_t nCluster = ClusterArray->GetEntriesFast(); + printf("\r Read %d clusters from entry %d", nCluster, iEntry); + // Loop through all TRD digits for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { - c = (AliTRDcluster*)clusterArray->UncheckedAt(iCluster); + c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster); AliTRDcluster *co = new AliTRDcluster(*c); co->SetSigmaY2(c->GetSigmaY2() * fSY2corr); + Int_t ltb = co->GetLocalTimeBin(); + if(ltb != 0) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr); + array->AddLast(co); - delete clusterArray->RemoveAt(iCluster); + delete ClusterArray->RemoveAt(iCluster); } } - file->Close(); - delete clusterArray; - savedir->cd(); - + delete ClusterArray; + savedir->cd(); + } -//___________________________________________________________________ -void AliTRDtracker::FindTracks(AliTRDtrackingSector* fTrSec, TH1F *hs, TH1F *hd) +//______________________________________________________________________ +void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename) { // - // Finds tracks in TRD + // Reads AliTRDclusters from file . The names of the cluster + // tree and branches should match the ones used in + // AliTRDclusterizer::WriteClusters() + // if == 0, clusters are added into AliTRDtracker fCluster array // - TH1F *hsame = hs; - TH1F *hdiff = hd; + TDirectory *savedir=gDirectory; - Int_t numOfTimeBins = fTrSec[0].GetNtimeBins(); - Int_t nseed=fSeeds->GetEntriesFast(); + TFile *file = TFile::Open(filename); + if (!file->IsOpen()) { + cerr<<"Can't open file with TRD clusters"<Get(treeName); - AliTRDtrack& t=*((AliTRDtrack*)fSeeds->UncheckedAt(i)); + if (!ClusterTree) { + cerr<<"AliTRDtracker::ReadClusters(): "; + cerr<<"can't get a tree with clusters !\n"; + return; + } - nSeedClusters = t.GetNclusters(); - Double_t alpha=t.GetAlpha(); + TObjArray *ClusterArray = new TObjArray(400); - if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); - if (alpha < 0. ) alpha += 2.*TMath::Pi(); - Int_t ns=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect; + ClusterTree->GetBranch("TRDcluster")->SetAddress(&ClusterArray); - Int_t label = GetTrackLabel(t); + Int_t nEntries = (Int_t) ClusterTree->GetEntries(); + cout<<"found "<RemoveAt(i); - fNseeds--; - } -} + for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) { -//__________________________________________________________________ -void AliTRDtracker::UseClusters(AliTRDtrack t) -{ - // - // Mark used cluster - // + // Import the tree + nbytes += ClusterTree->GetEvent(iEntry); - Int_t ncl=t.GetNclusters(); - for (Int_t i=0; iUncheckedAt(index); - c->Use(); + // Get the number of points in the detector + Int_t nCluster = ClusterArray->GetEntriesFast(); + printf("\r Read %d clusters from entry %d", nCluster, iEntry); + + // Loop through all TRD digits + for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { + c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster); + AliTRDcluster *co = new AliTRDcluster(*c); + co->SetSigmaY2(c->GetSigmaY2() * fSY2corr); + Int_t ltb = co->GetLocalTimeBin(); + if(ltb != 0) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr); + array->AddLast(co); + delete ClusterArray->RemoveAt(iCluster); + } } -} + file->Close(); + delete ClusterArray; + savedir->cd(); + +} + //__________________________________________________________________ -Int_t AliTRDtracker::GetTrackLabel(AliTRDtrack t) -{ - // - // Get MC label - // +void AliTRDtracker::CookLabel(AliKalmanTrack* pt, Float_t wrong) const { Int_t label=123456789, index, i, j; - Int_t ncl=t.GetNclusters(); - const Int_t kRange = AliTRDgeometry::kNplan * fGeom->GetTimeMax(); - Bool_t labelAdded; + Int_t ncl=pt->GetNumberOfClusters(); + const Int_t range = fTrSec[0]->GetOuterTimeBin()+1; + + Bool_t label_added; - // Int_t s[kRange][2]; - Int_t **s = new Int_t* [kRange]; - for (i=0; iGetClusterIndex(i); AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index); - t0=c->GetTrackIndex(0); - t1=c->GetTrackIndex(1); - t2=c->GetTrackIndex(2); + t0=c->GetLabel(0); + t1=c->GetLabel(1); + t2=c->GetLabel(2); } for (i=0; iGetClusterIndex(i); AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index); for (Int_t k=0; k<3; k++) { - label=c->GetTrackIndex(k); - labelAdded=kFALSE; j=0; + label=c->GetLabel(k); + label_added=kFALSE; j=0; if (label >= 0) { - while ( (!labelAdded) && ( j < kRange ) ) { + while ( (!label_added) && ( j < range ) ) { if (s[j][0]==label || s[j][1]==0) { s[j][0]=label; s[j][1]=s[j][1]+1; - labelAdded=kTRUE; + label_added=kTRUE; } j++; } @@ -899,60 +1574,701 @@ Int_t AliTRDtracker::GetTrackLabel(AliTRDtrack t) Int_t max=0; label = -123456789; - for (i=0; imax) { max=s[i][1]; label=s[i][0]; } } + + for (i=0; i ncl*fgkLabelFraction) return label; - else return -1; + + if ((1.- Float_t(max)/ncl) > wrong) label=-label; + + pt->SetLabel(label); + } -//___________________________________________________________________ -Int_t AliTRDtracker::WriteTracks(const Char_t *filename) + +//__________________________________________________________________ +void AliTRDtracker::UseClusters(const AliKalmanTrack* t, Int_t from = 0) const { + Int_t ncl=t->GetNumberOfClusters(); + for (Int_t i=from; iGetClusterIndex(i); + AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index); + c->Use(); + } +} + + +//_____________________________________________________________________ +Double_t AliTRDtracker::ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt) +{ + // Parametrised "expected" error of the cluster reconstruction in Y + + Double_t s = 0.08 * 0.08; + return s; +} + +//_____________________________________________________________________ +Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t r, Double_t tgl) { + // Parametrised "expected" error of the cluster reconstruction in Z + + Double_t s = 6 * 6 /12.; + return s; +} + + +//_____________________________________________________________________ +Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t local_tb) const +{ + // + // Returns radial position which corresponds to time bin + // in tracking sector and plane + // + + Int_t index = fTrSec[sector]->CookTimeBinIndex(plane, local_tb); + Int_t pl = fTrSec[sector]->GetLayerNumber(index); + return fTrSec[sector]->GetLayer(pl)->GetX(); + +} + + +//_______________________________________________________ +AliTRDtracker::AliTRDpropagationLayer::AliTRDpropagationLayer(Double_t x, + Double_t dx, Double_t rho, Double_t x0, Int_t tb_index) +{ // - // Write the tracks to the output file + // AliTRDpropagationLayer constructor // - TDirectory *savedir=gDirectory; + fN = 0; fX = x; fdX = dx; fRho = rho; fX0 = x0; + fClusters = NULL; fIndex = NULL; fTimeBinIndex = tb_index; + - //TFile *out=TFile::Open(filename,"RECREATE"); - TFile *out = (TFile*) gROOT->GetListOfFiles()->FindObject(filename); - if (!out) { - printf("AliTRDtracker::Open -- "); - printf("Open the ALIROOT-file %s.\n",filename); - out = new TFile(filename,"RECREATE"); + for(Int_t i=0; i < (Int_t) kZONES; i++) { + fZc[i]=0; fZmax[i] = 0; } - else { - printf("AliTRDtracker::Open -- "); - printf("%s is already open.\n",filename); + + fYmax = 0; + + if(fTimeBinIndex >= 0) { + fClusters = new (AliTRDcluster*)[kMAX_CLUSTER_PER_TIME_BIN]; + fIndex = new (UInt_t)[kMAX_CLUSTER_PER_TIME_BIN]; } - Char_t treeName[12]; - sprintf(treeName,"TreeT%d_TRD",fEvent); - TTree tracktree(treeName,"Tree with TRD tracks"); + fHole = kFALSE; + fHoleZc = 0; + fHoleZmax = 0; + fHoleYc = 0; + fHoleYmax = 0; + fHoleRho = 0; + fHoleX0 = 0; + +} + +//_______________________________________________________ +void AliTRDtracker::AliTRDpropagationLayer::SetHole( + Double_t Zmax, Double_t Ymax, Double_t rho, + Double_t x0, Double_t Yc, Double_t Zc) +{ + // + // Sets hole in the layer + // + + fHole = kTRUE; + fHoleZc = Zc; + fHoleZmax = Zmax; + fHoleYc = Yc; + fHoleYmax = Ymax; + fHoleRho = rho; + fHoleX0 = x0; +} + - AliTRDtrack *iotrack=0; - tracktree.Branch("tracks","AliTRDtrack",&iotrack,32000,0); +//_______________________________________________________ +AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, Int_t gs, AliTRDparameter* par) +{ + // + // AliTRDtrackingSector Constructor + // + + fGeom = geo; + fPar = par; + fGeomSector = gs; + fTzeroShift = 0.13; + fN = 0; + + for(UInt_t i=0; i < kMAX_TIME_BIN_INDEX; i++) fTimeBinIndex[i] = -1; + + + AliTRDpropagationLayer* ppl; + + Double_t x, xin, xout, dx, rho, x0; + Int_t steps; - Int_t ntracks=fTracks->GetEntriesFast(); + // set time bins in the gas of the TPC - for (Int_t i=0; iUncheckedAt(i); - iotrack=pt; - tracktree.Fill(); - cerr<<"WriteTracks: put track "<Close(); + // set time bins in the outer field cage vessel - savedir->cd(); + dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; x0 = 44.77; // Tedlar + ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1); + InsertLayer(ppl); - cerr<<"WriteTracks: done"<Rmin(); + + // add layers between TPC and TRD (Air temporarily) + xin = xout; xout = xtrd; + steps = 50; dx = (xout - xin)/steps; + rho = 1.2e-3; x0 = 36.66; + + for(Int_t i=0; iCroHght(); // Rohacell + Double_t dxSpace = (Double_t) fGeom->Cspace(); // Spacing between planes + Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region + Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region + Double_t dxRad = (Double_t) fGeom->CraHght(); // Radiator + Double_t dxTEC = dxRad + dxDrift + dxAmp + dxRo; + Double_t dxPlane = dxTEC + dxSpace; + + Int_t tbBefore = (Int_t) (dxAmp/fPar->GetTimeBinSize()); + if(tbBefore > fPar->GetTimeBefore()) tbBefore = fPar->GetTimeBefore(); + + Int_t tb, tb_index; + const Int_t nChambers = AliTRDgeometry::Ncham(); + Double_t Ymax = 0, holeYmax = 0, Zc[nChambers], Zmax[nChambers]; + Double_t holeZmax = 1000.; // the whole sector is missing + + + for(Int_t plane = 0; plane < AliTRDgeometry::Nplan(); plane++) { + + // Radiator + xin = xtrd + plane * dxPlane; xout = xin + dxRad; + steps = 12; dx = (xout - xin)/steps; rho = 0.074; x0 = 40.6; + for(Int_t i=0; iGetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) { + holeYmax = x*TMath::Tan(0.5*alpha); + ppl->SetHole(holeYmax, holeZmax); + } + if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) { + holeYmax = x*TMath::Tan(0.5*alpha); + ppl->SetHole(holeYmax, holeZmax); + } + InsertLayer(ppl); + } + + Ymax = fGeom->GetChamberWidth(plane)/2; + for(Int_t ch = 0; ch < nChambers; ch++) { + Zmax[ch] = fGeom->GetChamberLength(plane,ch)/2; + Float_t pad = fPar->GetRowPadSize(plane,ch,0); + Float_t row0 = fPar->GetRow0(plane,ch,0); + Int_t nPads = fPar->GetRowMax(plane,ch,0); + Zc[ch] = (pad * nPads)/2 + row0 - pad/2; + } + + dx = fPar->GetTimeBinSize(); + rho = 0.00295 * 0.85; x0 = 11.0; + + Double_t x0 = (Double_t) fPar->GetTime0(plane); + Double_t xbottom = x0 - dxDrift; + Double_t xtop = x0 + dxAmp; + + // Amplification region + + steps = (Int_t) (dxAmp/dx); + + for(tb = 0; tb < steps; tb++) { + x = x0 + tb * dx + dx/2; + tb_index = CookTimeBinIndex(plane, -tb-1); + ppl = new AliTRDpropagationLayer(x,dx,rho,x0,tb_index); + ppl->SetYmax(Ymax); + for(Int_t ch = 0; ch < nChambers; ch++) { + ppl->SetZmax(ch, Zc[ch], Zmax[ch]); + } + if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) { + holeYmax = x*TMath::Tan(0.5*alpha); + ppl->SetHole(holeYmax, holeZmax); + } + if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) { + holeYmax = x*TMath::Tan(0.5*alpha); + ppl->SetHole(holeYmax, holeZmax); + } + InsertLayer(ppl); + } + tb_index = CookTimeBinIndex(plane, -steps); + x = (x + dx/2 + xtop)/2; + dx = 2*(xtop-x); + ppl = new AliTRDpropagationLayer(x,dx,rho,x0,tb_index); + ppl->SetYmax(Ymax); + for(Int_t ch = 0; ch < nChambers; ch++) { + ppl->SetZmax(ch, Zc[ch], Zmax[ch]); + } + if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) { + holeYmax = x*TMath::Tan(0.5*alpha); + ppl->SetHole(holeYmax, holeZmax); + } + if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) { + holeYmax = x*TMath::Tan(0.5*alpha); + ppl->SetHole(holeYmax, holeZmax); + } + InsertLayer(ppl); + + // Drift region + dx = fPar->GetTimeBinSize(); + steps = (Int_t) (dxDrift/dx); + + for(tb = 0; tb < steps; tb++) { + x = x0 - tb * dx - dx/2; + tb_index = CookTimeBinIndex(plane, tb); + + ppl = new AliTRDpropagationLayer(x,dx,rho,x0,tb_index); + ppl->SetYmax(Ymax); + for(Int_t ch = 0; ch < nChambers; ch++) { + ppl->SetZmax(ch, Zc[ch], Zmax[ch]); + } + if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) { + holeYmax = x*TMath::Tan(0.5*alpha); + ppl->SetHole(holeYmax, holeZmax); + } + if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) { + holeYmax = x*TMath::Tan(0.5*alpha); + ppl->SetHole(holeYmax, holeZmax); + } + InsertLayer(ppl); + } + tb_index = CookTimeBinIndex(plane, steps); + x = (x - dx/2 + xbottom)/2; + dx = 2*(x-xbottom); + ppl = new AliTRDpropagationLayer(x,dx,rho,x0,tb_index); + ppl->SetYmax(Ymax); + for(Int_t ch = 0; ch < nChambers; ch++) { + ppl->SetZmax(ch, Zc[ch], Zmax[ch]); + } + if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) { + holeYmax = x*TMath::Tan(0.5*alpha); + ppl->SetHole(holeYmax, holeZmax); + } + if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) { + holeYmax = x*TMath::Tan(0.5*alpha); + ppl->SetHole(holeYmax, holeZmax); + } + InsertLayer(ppl); + + // Pad Plane + xin = xtop; dx = 0.025; xout = xin + dx; rho = 1.7; x0 = 33.0; + ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1); + if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) { + holeYmax = (xin+dx/2)*TMath::Tan(0.5*alpha); + ppl->SetHole(holeYmax, holeZmax); + } + if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) { + holeYmax = (xin+dx/2)*TMath::Tan(0.5*alpha); + ppl->SetHole(holeYmax, holeZmax); + } + InsertLayer(ppl); + + // Rohacell + xin = xout; xout = xtrd + (plane + 1) * dxPlane - dxSpace; + steps = 5; dx = (xout - xin)/steps; rho = 0.074; x0 = 40.6; + for(Int_t i=0; iGetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) { + holeYmax = x*TMath::Tan(0.5*alpha); + ppl->SetHole(holeYmax, holeZmax); + } + if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) { + holeYmax = x*TMath::Tan(0.5*alpha); + ppl->SetHole(holeYmax, holeZmax); + } + InsertLayer(ppl); + } + + // Space between the chambers, air + xin = xout; xout = xtrd + (plane + 1) * dxPlane; + steps = 5; dx = (xout - xin)/steps; rho = 1.29e-3; x0 = 36.66; + for(Int_t i=0; iGetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) { + holeYmax = x*TMath::Tan(0.5*alpha); + ppl->SetHole(holeYmax, holeZmax); + } + if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) { + holeYmax = x*TMath::Tan(0.5*alpha); + ppl->SetHole(holeYmax, holeZmax); + } + InsertLayer(ppl); + } + } + + // Space between the TRD and RICH + Double_t xRICH = 500.; + xin = xout; xout = xRICH; + steps = 200; dx = (xout - xin)/steps; rho = 1.29e-3; x0 = 36.66; + for(Int_t i=0; i in plane + // + + Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region + Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region + Double_t dx = (Double_t) fPar->GetTimeBinSize(); + + Int_t tbAmp = fPar->GetTimeBefore(); + Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx); + Int_t tbDrift = fPar->GetTimeMax(); + Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx); + + Int_t tb_per_plane = TMath::Min(tbAmp,maxAmp) + TMath::Min(tbDrift,maxDrift); + + Int_t gtb = (plane+1) * tb_per_plane - local_tb - 1; + + if((local_tb < 0) && + (TMath::Abs(local_tb) > TMath::Min(tbAmp,maxAmp))) return -1; + if(local_tb >= TMath::Min(tbDrift,maxDrift)) return -1; + + return gtb; + + +} + +//______________________________________________________ + +void AliTRDtracker::AliTRDtrackingSector::MapTimeBinLayers() +{ + // + // For all sensitive time bins sets corresponding layer index + // in the array fTimeBins + // + + Int_t index; + + for(Int_t i = 0; i < fN; i++) { + index = fLayers[i]->GetTimeBinIndex(); + + // printf("gtb %d -> pl %d -> x %f \n", index, i, fLayers[i]->GetX()); + + if(index < 0) continue; + if(index >= (Int_t) kMAX_TIME_BIN_INDEX) { + printf("*** AliTRDtracker::MapTimeBinLayers: \n"); + printf(" index %d exceeds allowed maximum of %d!\n", + index, kMAX_TIME_BIN_INDEX-1); + continue; + } + fTimeBinIndex[index] = i; + } + + Double_t x1, dx1, x2, dx2, gap; + + for(Int_t i = 0; i < fN-1; i++) { + x1 = fLayers[i]->GetX(); + dx1 = fLayers[i]->GetdX(); + x2 = fLayers[i+1]->GetX(); + dx2 = fLayers[i+1]->GetdX(); + gap = (x2 - dx2/2) - (x1 + dx1/2); + if(gap < -0.01) { + printf("*** warning: layers %d and %d are overlayed:\n",i,i+1); + printf(" %f + %f + %f > %f\n", x1, dx1/2, dx2/2, x2); + } + if(gap > 0.01) { + printf("*** warning: layers %d and %d have a large gap:\n",i,i+1); + printf(" (%f - %f) - (%f + %f) = %f\n", + x2, dx2/2, x1, dx1, gap); + } + } +} + + +//______________________________________________________ + + +Int_t AliTRDtracker::AliTRDtrackingSector::GetLayerNumber(Double_t x) const +{ + // + // Returns the number of time bin which in radial position is closest to + // + + if(x >= fLayers[fN-1]->GetX()) return fN-1; + if(x <= fLayers[0]->GetX()) return 0; + + Int_t b=0, e=fN-1, m=(b+e)/2; + for (; b fLayers[m]->GetX()) b=m+1; + else e=m; + } + if(TMath::Abs(x - fLayers[m]->GetX()) > + TMath::Abs(x - fLayers[m+1]->GetX())) return m+1; + else return m; + +} + +//______________________________________________________ + +Int_t AliTRDtracker::AliTRDtrackingSector::GetInnerTimeBin() const +{ + // + // Returns number of the innermost SENSITIVE propagation layer + // + + return GetLayerNumber(0); +} + +//______________________________________________________ + +Int_t AliTRDtracker::AliTRDtrackingSector::GetOuterTimeBin() const +{ + // + // Returns number of the outermost SENSITIVE time bin + // + + return GetLayerNumber(GetNumberOfTimeBins() - 1); } +//______________________________________________________ + +Int_t AliTRDtracker::AliTRDtrackingSector::GetNumberOfTimeBins() const +{ + // + // Returns number of SENSITIVE time bins + // + + Int_t tb, layer; + for(tb = kMAX_TIME_BIN_INDEX-1; tb >=0; tb--) { + layer = GetLayerNumber(tb); + if(layer>=0) break; + } + return tb+1; +} + +//______________________________________________________ + +void AliTRDtracker::AliTRDtrackingSector::InsertLayer(AliTRDpropagationLayer* pl) +{ + // + // Insert layer in fLayers array. + // Layers are sorted according to X coordinate. + + if ( fN == ((Int_t) kMAX_LAYERS_PER_SECTOR)) { + printf("AliTRDtrackingSector::InsertLayer(): Too many layers !\n"); + return; + } + if (fN==0) {fLayers[fN++] = pl; return;} + Int_t i=Find(pl->GetX()); + + memmove(fLayers+i+1 ,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*)); + fLayers[i]=pl; fN++; + +} + +//______________________________________________________ + +Int_t AliTRDtracker::AliTRDtrackingSector::Find(Double_t x) const +{ + // + // Returns index of the propagation layer nearest to X + // + + if (x <= fLayers[0]->GetX()) return 0; + if (x > fLayers[fN-1]->GetX()) return fN; + Int_t b=0, e=fN-1, m=(b+e)/2; + for (; b fLayers[m]->GetX()) b=m+1; + else e=m; + } + return m; +} + +//______________________________________________________ + +void AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters( + Double_t y, Double_t z, Double_t &dx, Double_t &rho, Double_t &x0, + Bool_t &lookForCluster) const +{ + // + // Returns radial step , density , rad. length , + // and sensitivity in point + // + + dx = fdX; + rho = fRho; + x0 = fX0; + lookForCluster = kFALSE; + + // check dead regions + if(fTimeBinIndex >= 0) { + for(Int_t ch = 0; ch < (Int_t) kZONES; ch++) { + if(TMath::Abs(z - fZc[ch]) < fZmax[ch]) + lookForCluster = kTRUE; + } + if(TMath::Abs(y) > fYmax) lookForCluster = kFALSE; + if(!lookForCluster) { + // rho = 1.7; x0 = 33.0; // G10 + } + } + + // check hole + if(fHole && (TMath::Abs(y - fHoleYc) < fHoleYmax) && + (TMath::Abs(z - fHoleZc) < fHoleZmax)) { + lookForCluster = kFALSE; + rho = fHoleRho; + x0 = fHoleX0; + } + + return; +} + +//______________________________________________________ + +void AliTRDtracker::AliTRDpropagationLayer::InsertCluster(AliTRDcluster* c, + UInt_t index) { + +// Insert cluster in cluster array. +// Clusters are sorted according to Y coordinate. + + if(fTimeBinIndex < 0) { + printf("*** attempt to insert cluster into non-sensitive time bin!\n"); + return; + } + + if (fN== (Int_t) kMAX_CLUSTER_PER_TIME_BIN) { + printf("AliTRDpropagationLayer::InsertCluster(): 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(Double_t y) const { + +// Returns index of the cluster nearest in Y + + if (y <= fClusters[0]->GetY()) return 0; + if (y > fClusters[fN-1]->GetY()) return fN; + Int_t b=0, e=fN-1, m=(b+e)/2; + for (; b fClusters[m]->GetY()) b=m+1; + else e=m; + } + return m; +} + + + + + + + + + diff --git a/TRD/AliTRDtracker.h b/TRD/AliTRDtracker.h index 68e91fef04a..4dbc91d36ac 100644 --- a/TRD/AliTRDtracker.h +++ b/TRD/AliTRDtracker.h @@ -4,14 +4,7 @@ /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * See cxx source for full Copyright notice */ -//////////////////////////////////////////////////////////////////////////////// -// // -// The TRD tracker // -// // -//////////////////////////////////////////////////////////////////////////////// - -#include -#include +#include "AliTracker.h" class TFile; class TParticle; @@ -19,60 +12,187 @@ class TParticlePDG; class TObjArray; class AliTRDgeometry; +class AliTRDparameter; class AliTRDtrack; +class AliTRDcluster; class AliTRDmcTrack; -class AliTRDtrackingSector; -class AliTRDtracker : public TNamed { +const unsigned kMAX_LAYERS_PER_SECTOR = 1000; +const unsigned kMAX_TIME_BIN_INDEX = 216; // (30 drift + 6 ampl) * 6 planes +const unsigned kMAX_CLUSTER_PER_TIME_BIN = 3500; +const unsigned kZONES = 5; +const Int_t kTRACKING_SECTORS = 18; + +class AliTRDtracker : public AliTracker { public: - AliTRDtracker(); - AliTRDtracker(const Text_t* name, const Text_t* title); - virtual ~AliTRDtracker(); - - virtual void Clusters2Tracks(TH1F *hs, TH1F *hd); - 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 FindProlongation(AliTRDtrack& t, AliTRDtrackingSector *sec, - Int_t s, Int_t rf=0, Int_t matched_index = -1, - TH1F *hs=0, TH1F *hd=0); - void GetEvent(const Char_t *hitfile, const Char_t *clusterfile); - void SetUpSectors(AliTRDtrackingSector *sec); - virtual void MakeSeeds(Int_t inner, Int_t outer, AliTRDtrackingSector *sec, - Int_t turn, TH1F *hs, TH1F *hd); - virtual void FindTracks(AliTRDtrackingSector *sec, TH1F *hs, TH1F *hd); - virtual void UseClusters(AliTRDtrack t); - virtual Int_t GetTrackLabel(AliTRDtrack t); - Int_t WriteTracks(const Char_t *filename); - void ReadClusters(TObjArray *array, const Char_t *filename); - - Float_t GetSeedGap() const {return fgkSeedGap;} - Float_t GetSeedStep() const {return fgkSeedStep;} - Float_t GetSeedDepth() const {return fgkSeedDepth;} - Float_t GetSkipDepth() const {return fgkSkipDepth;} - Double_t GetMaxChi2() const {return fgkMaxChi2;} - Float_t GetMaxSeedC() const {return fgkMaxSeedC;} - Float_t GetMaxSeedTan() const {return fgkMaxSeedTan;} - Double_t GetSeedErrorSY() const {return fgkSeedErrorSY;} - Double_t GetSeedErrorSY3() const {return fgkSeedErrorSY3;} - Double_t GetSeedErrorSZ() const {return fgkSeedErrorSZ;} - Float_t GetLabelFraction() const {return fgkLabelFraction;} - Float_t GetWideRoad() const {return fgkWideRoad;} - - Float_t GetMinClustersInTrack() const {return fgkMinClustersInTrack;} - Float_t GetMinClustersInSeed() const {return fgkMinClustersInSeed;} - Float_t GetMaxSeedDeltaZ() const {return fgkMaxSeedDeltaZ;} - Float_t GetMaxSeedVertexZ() const {return fgkMaxSeedVertexZ;} - - void SetSY2corr(Float_t w) {fSY2corr = w;} + AliTRDtracker():AliTracker() {} + AliTRDtracker(const TFile *in); + ~AliTRDtracker(); + + Int_t Clusters2Tracks(const TFile *in, TFile *out); + Int_t PropagateBack(const TFile *in, TFile *out); + AliCluster *GetCluster(Int_t index) const { return NULL; }; + virtual void CookLabel(AliKalmanTrack *t,Float_t wrong) const; + virtual void UseClusters(const AliKalmanTrack *t, Int_t from=0) const; + + void SetEventNumber(Int_t event) { fEvent = event; } + void SetAddTRDseeds() { fAddTRDseeds = kTRUE; } + + void ReadClusters(TObjArray *array, const Char_t *filename); + Int_t CookSectorIndex(Int_t gs) { return kTRACKING_SECTORS - 1 - gs; } + + + Float_t GetSeedGap() const {return fSeedGap;} + Int_t GetMaxGap() const {return fMaxGap;} + Int_t GetTimeBinsPerPlane() const {return fTimeBinsPerPlane;} + Float_t GetSeedStep() const {return fSeedStep;} + Float_t GetSeedDepth() const {return fSeedDepth;} + Float_t GetSkipDepth() const {return fSkipDepth;} + Double_t GetMaxChi2() const {return fMaxChi2;} + Float_t GetMaxSeedC() const {return fMaxSeedC;} + Float_t GetMaxSeedTan() const {return fMaxSeedTan;} + Double_t GetSeedErrorSY() const {return fSeedErrorSY;} + Double_t GetSeedErrorSY3() const {return fSeedErrorSY3;} + Double_t GetSeedErrorSZ() const {return fSeedErrorSZ;} + Float_t GetLabelFraction() const {return fLabelFraction;} + Float_t GetWideRoad() const {return fWideRoad;} + + Float_t GetMinClustersInTrack() const {return fMinClustersInTrack;} + Float_t GetMinClustersInSeed() const {return fMinClustersInSeed;} + Float_t GetMaxSeedDeltaZ() const {return fMaxSeedDeltaZ;} + Float_t GetMaxSeedVertexZ() const {return fMaxSeedVertexZ;} + + // x <-> timebin conversions useful in analysis macros + Double_t GetX(Int_t sec, Int_t plane, Int_t local_tb) const; + Double_t GetX(Int_t sec, Int_t pl) const { + return fTrSec[sec]->GetLayer(pl)->GetX(); } + Int_t GetGlobalTimeBin(Int_t sec, Int_t plane, Int_t local_tb) const { + return fTrSec[sec]->CookTimeBinIndex(plane,local_tb); } + Double_t GetLayerNumber(Int_t sec, Double_t x) const { + return fTrSec[sec]->GetLayerNumber(x); } + + public: + class AliTRDpropagationLayer { + // ***************** internal class ******************* + public: + AliTRDpropagationLayer(Double_t x, Double_t dx, Double_t rho, + Double_t x0, Int_t tb_index); + + ~AliTRDpropagationLayer() { + if(fTimeBinIndex >= 0) { delete[] fClusters; delete[] fIndex; } + } + void InsertCluster(AliTRDcluster*, UInt_t); + operator Int_t() const {return fN;} + AliTRDcluster* operator[](Int_t i) {return fClusters[i];} + UInt_t GetIndex(Int_t i) const {return fIndex[i];} + Double_t GetX() const { return fX; } + Double_t GetdX() const { return fdX; } + Double_t GetRho() const { return fRho; } + Double_t GetX0() const { return fX0; } + Int_t GetTimeBinIndex() const { return fTimeBinIndex; } + void GetPropagationParameters(Double_t y, Double_t z, + Double_t &dx, Double_t &rho, Double_t &x0, + Bool_t &lookForCluster) const; + Int_t Find(Double_t y) const; + void SetZmax(Int_t cham, Double_t center, Double_t w) + { fZc[cham] = center; fZmax[cham] = w; } + void SetYmax(Double_t w) { fYmax = w; } + 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]; } + + void SetHole(Double_t Zmax, Double_t Ymax, + Double_t rho = 1.29e-3, Double_t x0 = 36.66, + Double_t Yc = 0, Double_t Zc = 0); + + void Clear() {for(Int_t i=0; iGetX(); } + void MapTimeBinLayers(); + Int_t GetLayerNumber(Double_t x) const; + Int_t GetInnerTimeBin() const; + Int_t GetOuterTimeBin() const; + Int_t GetLayerNumber(Int_t tb) const {return fTimeBinIndex[tb];} + Float_t GetTzeroShift() const { return fTzeroShift; } + Int_t Find(Double_t x) const; + void InsertLayer(AliTRDpropagationLayer* pl); + // AliTRDpropagationLayer* operator[](Int_t i) { return fLayers[i]; } + AliTRDpropagationLayer* GetLayer(Int_t i) { return fLayers[i]; } + Int_t CookTimeBinIndex(Int_t plane, Int_t local_tb) const; + + private: + Int_t fN; // total number of layers + AliTRDgeometry *fGeom; + AliTRDparameter *fPar; + AliTRDpropagationLayer *fLayers[kMAX_LAYERS_PER_SECTOR]; + Int_t fTimeBinIndex[kMAX_TIME_BIN_INDEX]; + Float_t fTzeroShift; // T0 shift in cm + Int_t fGeomSector; // sector # in AliTRDgeometry + }; + + private: + + void LoadEvent(); + void UnloadEvent(); + + virtual void MakeSeeds(Int_t inner, Int_t outer, Int_t turn); + + Int_t FollowProlongation(AliTRDtrack& t, Int_t rf); + Int_t FollowBackProlongation(AliTRDtrack& t); + + Int_t PropagateToTPC(AliTRDtrack& t); + Int_t PropagateToOuterPlane(AliTRDtrack& t, Double_t x); + + Int_t WriteTracks(); + void ReadClusters(TObjArray *array, const TFile *in = 0); + + void SetSY2corr(Float_t w) {fSY2corr = w;} + void SetSZ2corr(Float_t w) {fSZ2corr = w;} + Double_t ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt); + Double_t ExpectedSigmaZ2(Double_t r, Double_t tgl); + protected: - Int_t fEvent; // Event number + Int_t fEvent; // Event number - AliTRDgeometry *fGeom; // Pointer to TRD geometry + AliTRDgeometry *fGeom; // Pointer to TRD geometry + AliTRDparameter *fPar; + AliTRDtrackingSector *fTrSec[kTRACKING_SECTORS]; + // array of tracking sectors; + Int_t fNclusters; // Number of clusters in TRD TObjArray *fClusters; // List of clusters for all sectors @@ -85,28 +205,43 @@ class AliTRDtracker : public TNamed { Float_t fSY2corr; // Correction coefficient for // cluster SigmaY2 - static const Float_t fgkSeedGap; // Distance between inner and outer + Float_t fSZ2corr; // Correction coefficient for + // cluster SigmaZ2 + + static const Float_t fSeedGap; // Distance between inner and outer // time bin in seeding // (fraction of all time bins) - static const Float_t fgkSeedStep; // Step in iterations - static const Float_t fgkSeedDepth; // Fraction of TRD allocated for seeding - static const Float_t fgkSkipDepth; // Fraction of TRD which can be skipped + static const Float_t fSeedStep; // Step in iterations + static const Float_t fSeedDepth; // Fraction of TRD allocated for seeding + static const Float_t fSkipDepth; // Fraction of TRD which can be skipped // in track prolongation - static const Double_t fgkMaxChi2; // max increment in track chi2 + Int_t fTimeBinsPerPlane; // number of sensitive timebins per plane + Int_t fMaxGap; // max gap (in time bins) in the track + // in track prolongation + + static const Double_t fMaxChi2; // max increment in track chi2 - static const Float_t fgkMinClustersInTrack; // min fraction of clusters in track - static const Float_t fgkMinClustersInSeed; // min fraction of clusters in seed - static const Float_t fgkMaxSeedDeltaZ; // max dZ in MakeSeeds - static const Float_t fgkMaxSeedDeltaZ12; // max abs(z1-z2) in MakeSeeds - static const Float_t fgkMaxSeedC; // max initial curvature in MakeSeeds - static const Float_t fgkMaxSeedTan; // max initial Tangens(lambda) in MakeSeeds - static const Float_t fgkMaxSeedVertexZ; // max vertex Z in MakeSeeds - static const Double_t fgkSeedErrorSY; // sy parameter in MakeSeeds - static const Double_t fgkSeedErrorSY3; // sy3 parameter in MakeSeeds - static const Double_t fgkSeedErrorSZ; // sz parameter in MakeSeeds - static const Float_t fgkLabelFraction; // min fraction of clusters in GetTrackLabel - static const Float_t fgkWideRoad; // max road width in FindProlongation + static const Float_t fMinClustersInTrack; // min number of clusters in track + // out of total timebins + + static const Float_t fMinFractionOfFoundClusters; // min found clusters + // out of expected + + static const Float_t fMinClustersInSeed; // min fraction of clusters in seed + static const Float_t fMaxSeedDeltaZ; // max dZ in MakeSeeds + static const Float_t fMaxSeedDeltaZ12; // max abs(z1-z2) in MakeSeeds + static const Float_t fMaxSeedC; // max initial curvature in MakeSeeds + static const Float_t fMaxSeedTan; // max initial Tangens(lambda) in MakeSeeds + static const Float_t fMaxSeedVertexZ; // max vertex Z in MakeSeeds + static const Double_t fSeedErrorSY; // sy parameter in MakeSeeds + static const Double_t fSeedErrorSY3; // sy3 parameter in MakeSeeds + static const Double_t fSeedErrorSZ; // sz parameter in MakeSeeds + static const Float_t fLabelFraction; // min fraction of same label + static const Float_t fWideRoad; // max road width in FindProlongation + + Bool_t fVocal; + Bool_t fAddTRDseeds; ClassDef(AliTRDtracker,1) // manager base class diff --git a/TRD/AliTRDtrackingSector.cxx b/TRD/AliTRDtrackingSector.cxx index 22521eaeedd..7c663266854 100644 --- a/TRD/AliTRDtrackingSector.cxx +++ b/TRD/AliTRDtrackingSector.cxx @@ -15,6 +15,9 @@ /* $Log$ +Revision 1.8 2002/03/28 14:59:07 cblume +Coding conventions + Revision 1.7 2001/11/19 08:44:08 cblume Fix bugs reported by Rene @@ -56,6 +59,7 @@ Add the tracking code #include "AliTRDcluster.h" #include "AliTRDtimeBin.h" #include "AliTRDtrackingSector.h" +#include "AliTRDparameter.h" ClassImp(AliTRDtrackingSector) @@ -91,13 +95,17 @@ void AliTRDtrackingSector::SetUp() AliTRD *trd = (AliTRD*) gAlice->GetDetector("TRD"); fGeom = trd->GetGeometry(); - fTimeBinSize = fGeom->GetTimeBinSize(); - fN = AliTRDgeometry::Nplan() * (Int_t(AliTRDgeometry::DrThick() /fTimeBinSize) + 1); fTimeBin = new AliTRDtimeBin[fN]; + if (!fPar) { + fPar = new AliTRDparameter("TRDparameter","Standard TRD parameter"); + } + + fTimeBinSize = fPar->GetTimeBinSize(); + } //______________________________________________________ @@ -117,7 +125,7 @@ Double_t AliTRDtrackingSector::GetX(Int_t tb) const Int_t localTb = tbPerPlane - tb%tbPerPlane - 1; Int_t plane = tb/tbPerPlane; - Float_t t0 = fGeom->GetTime0(plane); + Float_t t0 = fPar->GetTime0(plane); Double_t x = t0 - (localTb + 0.5) * fTimeBinSize; return x; @@ -133,8 +141,8 @@ Int_t AliTRDtrackingSector::GetTimeBinNumber(Double_t x) const // Returns the time bin number // - Float_t rOut = fGeom->GetTime0(AliTRDgeometry::Nplan()-1); - Float_t rIn = fGeom->GetTime0(0) - AliTRDgeometry::DrThick(); + Float_t rOut = fPar->GetTime0(AliTRDgeometry::Nplan()-1); + Float_t rIn = fPar->GetTime0(0) - AliTRDgeometry::DrThick(); if(x >= rOut) return fN-1; @@ -142,11 +150,11 @@ Int_t AliTRDtrackingSector::GetTimeBinNumber(Double_t x) const Int_t plane; for (plane = AliTRDgeometry::Nplan()-1; plane >= 0; plane--) { - if(x > (fGeom->GetTime0(plane) - AliTRDgeometry::DrThick())) break; + if(x > (fPar->GetTime0(plane) - AliTRDgeometry::DrThick())) break; } Int_t tbPerPlane = fN/AliTRDgeometry::Nplan(); - Int_t localTb = Int_t((fGeom->GetTime0(plane)-x)/fTimeBinSize); + Int_t localTb = Int_t((fPar->GetTime0(plane)-x)/fTimeBinSize); if((localTb < 0) || (localTb >= tbPerPlane)) { printf("AliTRDtrackingSector::GetTimeBinNumber: \n"); @@ -204,13 +212,13 @@ Bool_t AliTRDtrackingSector::TECframe(Int_t tb, Double_t y, Double_t z) const for(Int_t iCha = 1; iCha < AliTRDgeometry::Ncham(); iCha++) { - fRow0 = fGeom->GetRow0(plane,iCha-1,0); - fRowPadSize = fGeom->GetRowPadSize(plane,iCha-1,0); - nPadRows = fGeom->GetRowMax(plane,iCha-1,0); + fRow0 = fPar->GetRow0(plane,iCha-1,0); + fRowPadSize = fPar->GetRowPadSize(plane,iCha-1,0); + nPadRows = fPar->GetRowMax(plane,iCha-1,0); zmin = fRow0 - fRowPadSize/2 + fRowPadSize * nPadRows; - fRow0 = fGeom->GetRow0(plane,iCha,0); - fRowPadSize = fGeom->GetRowPadSize(plane,iCha,0); + fRow0 = fPar->GetRow0(plane,iCha,0); + fRowPadSize = fPar->GetRowPadSize(plane,iCha,0); zmax = fRow0 - fRowPadSize/2; if((z > zmin) && (z < zmax)) return kTRUE; diff --git a/TRD/AliTRDtrackingSector.h b/TRD/AliTRDtrackingSector.h index 2ab16a47fda..49996d4aa03 100644 --- a/TRD/AliTRDtrackingSector.h +++ b/TRD/AliTRDtrackingSector.h @@ -16,6 +16,7 @@ class AliTRDtimeBin; class AliTRDgeometry; +class AliTRDparameter; class AliTRDtrackingSector : public TObject { @@ -34,14 +35,18 @@ public: Int_t GetTimeBin(Int_t det, Int_t local_tb) const; Bool_t TECframe(Int_t tb, Double_t y, Double_t z) const; + virtual void SetParameter(AliTRDparameter *par) { fPar = par; }; + AliTRDparameter *GetParameter() const { return fPar; }; + protected: Int_t fN; // ??????? AliTRDgeometry *fGeom; // Pointer to TRD geometry AliTRDtimeBin *fTimeBin; // Pointer to array of AliTRDtimeBin Float_t fTimeBinSize; // Time bin size in cm + AliTRDparameter *fPar; // TRD parameter - ClassDef(AliTRDtrackingSector,1) // Provides tools to address clusters which lay within one sector + ClassDef(AliTRDtrackingSector,2) // Provides tools to address clusters which lay within one sector }; diff --git a/TRD/Makefile b/TRD/Makefile index 5da3919d1ed..ed1a6ed53be 100644 --- a/TRD/Makefile +++ b/TRD/Makefile @@ -39,8 +39,6 @@ SRCS = AliTRD.cxx \ AliTRDarrayI.cxx \ AliTRDarrayF.cxx \ AliTRDpoints.cxx \ - AliTRDtimeBin.cxx \ - AliTRDtrackingSector.cxx \ AliTRDtrackHits.cxx \ AliTRDtrack.cxx \ AliTRDtracker.cxx \ diff --git a/TRD/TRDLinkDef.h b/TRD/TRDLinkDef.h index 0307727b201..b518bf61a74 100644 --- a/TRD/TRDLinkDef.h +++ b/TRD/TRDLinkDef.h @@ -32,9 +32,7 @@ #pragma link C++ class AliTRDdataArrayF+; #pragma link C++ class AliTRDsim+; #pragma link C++ class AliTRDpoints+; -#pragma link C++ class AliTRDtimeBin+; -#pragma link C++ class AliTRDtrackingSector+; -#pragma link C++ class AliTRDtrack-; +#pragma link C++ class AliTRDtrack+; #pragma link C++ class AliTRDtracker+; #pragma link C++ class AliTRDtrackHits+; #pragma link C++ class AliTRDcluster+; -- 2.43.0